CN101383653B - Automatic time synchronizing method on satellite - Google Patents

Abstract

An on-board autonomous time statistics method, which endows the existing time management unit with long-term, high-precision time keeping capabilities and autonomous phase adjustment capabilities. The time difference processor uses the time difference observations to estimate the periodicity of the latest filtering results in real time with a time difference EKF filter The time difference processor broadcasts the forecast time difference to the user on the one hand, and at the same time judges whether the forecast time difference at the current moment enters the phase modulation interval. If it enters the phase modulation interval, it calculates the phase modulation amount and the nominal adjustment The phase residual is used to drive the second pulse generator to adjust the corresponding phase according to the calculated phasing amount, and replace the initial forecast time difference and the a priori estimated value of the time difference in the on-board time difference prediction model with the nominal phase modulation residual, thus realizing The autonomous time on the star is unified. The invention can reduce the dependence on the ground station, help to reduce the requirement on the performance of the on-board clock, save the investment and maintenance cost of the constellation system, and improve the reliability and safety of the constellation system.

Description

A Method of Autonomous Time System on Satellite

technical field

The invention relates to an on-board time system method, in particular to an autonomous synchronization method of an on-board clock and a ground system clock, and is compatible with traditional ground station intervention time system methods, and is especially suitable for autonomous time system processing of satellite constellation systems.

Background technique

The output frequency of the clock (or frequency source) has main error factors such as frequency calibration deviation (called frequency difference or accuracy), aging phenomenon and random noise. Therefore, various pulse signals ( represents a time interval) and the corresponding ideal value always has some degree of phase deviation (or time difference), and as time goes by, this phase deviation becomes larger and larger. In order to realize some interactive functions, the phase between the pulses output by each frequency source set in different systems must maintain a certain precision synchronization, or, it is not necessary to adjust the real phase difference, but it is necessary to broadcast the pulse phase difference to the relevant users Adjust accordingly. Among them, the adjustment of the pulse phase is a physical process, called physical synchronization; without phase adjustment, just broadcasting the phase difference for synchronization is a mathematical process, called mathematical synchronization. Physical synchronization and mathematical synchronization constitute the two major contents of time synchronization.

The key to maintaining high-precision synchronization between the on-board clock and the ground-based system clock lies in the accurate tracking of the satellite-ground time difference. The methods of tracking satellite-ground time difference mainly include: 1) time difference observation, and the typical three types of methods are: the ground station establishes a ranging link with the transiting satellite to measure the satellite-ground time difference; the two-way time comparison measurement based on the inter-satellite link The inter-satellite time difference, supplemented by the satellite-ground time difference observation, can be converted to obtain the time difference observation of each constellation star relative to the ground system clock; based on the time service of the GPS global navigation system, the ground system clock and the satellite clock relative to the navigation system clock The deviation is converted into the observed amount of time difference between the satellite and the earth. 2) According to the understanding of the change law of the satellite-ground time difference, the time difference prediction model is used for forecasting. For the high stability crystal oscillator (QZ clock), the long-term forecast should adopt the secondary time difference prediction model including the influence of the aging rate, while for the atomic clock ( For example, RB or CS clock), it is not necessary to consider the influence of aging rate, and only use a linear time difference prediction model. The above two types of tracking methods have their own advantages and disadvantages: the former is limited by the observation area and time (for example, the distribution of ground stations is limited; the frequency of satellite-to-earth measurement is limited), and there are observation errors in time-difference observations, etc.; the latter It cannot take into account the random changing part of the clock (random noise), and relies heavily on the estimation accuracy of model variables (initial forecast time difference, frequency difference and aging rate), and the long-term time difference forecast error will grow rapidly at any time.

The main disadvantage of the time statistics method based on the intervention of the ground station is that the ability to accurately track the time difference between the satellite and the ground is weak, which leads to excessive requirements on hardware performance and relies heavily on the intervention of the ground station. Due to the limited distribution area of my country's ground stations, the satellite transit time is short, and the time in the invisible interval is very long, the ground station can only observe the satellite-ground time difference during the short period of satellite transit, and estimate the satellite time difference forecast accordingly In the non-visible interval, the long-term time difference forecast can only be performed by the time difference forecasting model. Therefore, the trend of rapid increase of the time difference forecast error at any time cannot be avoided. In order to meet the requirements of high-precision time synchronization indicators, the traditional method relies heavily on high-performance clock products, and the ground intervention time system processing cannot be carried out independently from the time difference observation in time, but can only be performed by the ground when the satellite is in transit. The station issues instructions such as time service, centralized time correction, uniform time correction or second pulse phase adjustment to complete the corresponding satellite-ground time synchronization process.

The patent application number is 200810074975.3, which introduces "a satellite constellation high-precision time unification method", which provides a solution for the time synchronization of a small formation satellite system composed of three stars (1 master and 2 slaves), representing my country's The latest application research results in constellation timing technology. The main purpose of this invention is to realize the time synchronization between the three satellites. Therefore, although it is proposed to use the domestic ground station to observe the time difference between the three satellites, the purpose is to control the time difference of each satellite within the required range of the index. This method is still a time system method that mainly relies on the intervention of the ground station. The reason is that, as the key equipment for on-board time system processing, the on-board time management unit in this invention (as shown in Figure 1, includes frequency source, frequency Synthetic splitter, time difference processor, second pulse generator and communication controller) only play a role in cooperating with the management of the ground station, but do not have the ability of autonomous precise timekeeping and autonomous phase modulation on the star, but can only passively accept The variable estimation (for example, frequency difference) uploaded by the ground station after processing, in order to correct the time difference forecast on the satellite, rely on the ground station to issue a time difference adjustment command for centralized time correction or uniform time correction processing, in order to adjust the second pulse phase on the satellite, rely on the ground station The phasing command issued.

In addition, this invention also proposes to use the two inter-satellite links between the primary star and the secondary star to observe the inter-satellite time difference, and use the GPS navigation system to provide time to each star to observe the inter-satellite time difference, and use these time difference observations to correct the star time difference. Method for uploading time difference forecasts (for broadcasting). These two methods have a certain ability of autonomous mathematical synchronization on the star, but the geometric distance between the primary and secondary stars is only a few hundred kilometers, and the inter-satellite time difference observations are directly used (without any filtering processing) to correct the inter-satellite forecast time difference The method is likely not to meet the synchronization requirements of constellation systems separated by more than thousands of kilometers. Moreover, these two methods do not have the capability of autonomous physical synchronization on the star. More importantly, these two methods cannot complete the autonomous synchronization of each satellite clock and the ground system clock.

At present, the timing tasks of the well-known foreign satellite constellation systems all rely on the ground station systems deployed around the world to complete. For example, the Iridium system with the inter-satellite link has 11 globally distributed ground stations; the Globalstar system needs more than 300 ground stations because there is no inter-satellite link; the global navigation system GPS also relies on the ground station system deployed around the world, but, GPS The system has installed inter-satellite links starting from Block II R satellites, and will form a certain inter-satellite autonomous time system capability after 2010. However, according to the current data, the system also does not have the ability of autonomous physical synchronization, and as a navigation system, there is no need for precise satellite-to-ground timing.

The on-board time system method based on ground station control will cause a considerable workload and high cost burden on the ground station for a constellation system with dozens of satellites. Compared with the big countries in Europe and the United States, my country has limited funds for developing constellation systems and the deployment of ground stations is restricted by regions. Therefore, it is very necessary to have independent time synchronization technology on the satellite.

Contents of the invention

The technology of the present invention solves the problem: overcomes the deficiencies of the prior art, and provides an autonomous time system method on the star. The present invention can independently realize the mathematical synchronization and physical synchronization on the star, which helps to reduce the dependence on the ground station, reduces the The requirements for on-board clock performance save the investment and maintenance costs of the constellation system, and at the same time improve the reliability, safety and operability of the constellation system.

Technical solution of the present invention: a kind of autonomous time system method on the star, comprises the following steps:

(A) the time difference processor judges the time difference adjustment mark of the current moment, if the time difference adjustment mark is "0", then enters step (B) and carries out self-timekeeping; if the time difference adjustment mark is "1", then enters the ground timing processing flow; If the time difference adjustment flag is "2", enter the process of centralized time adjustment on the ground;

(B) The time difference processor judges the on-board time-difference prediction model variable replacement flag, if the replacement flag is "1", then use the latest real-time estimation result of the on-board EKF filter to replace the corresponding variable of the current secondary time-difference prediction model, and use the The current moment is the initial forecasting time, and the time difference forecast model variable replacement sign is reset to "0", and then enters step (C); if the replacement sign is "0", then directly proceeds to step (C);

(C) The time difference processor judges the ground intervention phase modulation flag, if the phase modulation flag is "1", then set the current phase modulation flag representing the phase modulation process at the current moment to "1", and reset the ground intervention phase modulation flag Set to "0", then go to step (E); if the phase modulation flag is "0", then go to step (D);

(D) The time difference processor judges the on-board autonomous phasing flag, if the on-satellite autonomous phasing flag is "0", then set the current phasing flag to "0", and then go to step (G); If the phase modulation flag is "1", set the current phase modulation flag to "1", reset the autonomous phase modulation flag to "0", and then go to step (E);

(E) Calculate the phasing amount and the nominal phasing residual according to the current forecast time difference, and then turn to step (F);

(F) The time-of-flight processor issues a second-pulse phase modulation command to drive the second-pulse generator to perform corresponding phase adjustments according to the calculated phase-modulation amount, and unconditionally replace the initial time-difference amount in the secondary time-difference prediction model with the nominal phase-modulation residual , and take the current phasing moment as the initial forecasting moment, and then turn to step (G);

(G) Based on the secondary time difference prediction model on the star, the time difference between the satellite and the ground is forecasted according to the required time on the star, and the forecast time difference is broadcast to the user on the one hand, and on the other hand is used for the judging step (H) of the timing of the autonomous phase adjustment;

(H) the time difference processor judges whether the time difference forecast at the next moment is in the phase modulation interval, if "Yes", then the self-phase modulation timing sign is set to "1", then proceeds to step (1), otherwise, directly proceeds to Enter step (1);

(1) the time difference processor utilizes the EKF filter to carry out state update and time update to the time difference amount, frequency difference and aging rate of the secondary time difference forecasting model, then proceed to step (J);

(J) The time difference processor judges whether the current moment needs to utilize the state update and the time update result of the EKF filter to replace the variable estimation in the secondary time difference forecasting model according to the replacement period, if the current moment is the replacement time, then the replacement sign is set to "1", otherwise set the replacement flag to "0", and then go to step (A).

The processing flow of the time service in the step (A) is: when performing the time service operation, the time difference processor unconditionally accepts the satellite time information in the time service data as the current star time, and accepts the time difference amount in the time service data as the current moment forecast time difference.

The processing flow of the centralized time correction in the step (A) is: when the centralized time correction operation is performed, the time difference processor adjusts the star time according to the star time adjustment amount in the centralized time correction data at one time, and adjusts the time according to the centralized time correction data. The time difference adjustment amount in is to adjust the forecast time difference calculated at the current moment at one time.

In the step (B), the time difference processor utilizes the latest real-time estimation result of the EKF filter on the star to replace the corresponding variable of the current secondary time difference prediction model. Relative to the time difference observation sequence of the reference clock, the EKF filter algorithm is used to give the latest estimated value of the quadratic time difference forecast model variable in real time. The time difference processor judges whether the current time is the replacement time according to the replacement period of the model variable. If it is the variable replacement time, Then use the latest estimated value to replace the estimated value of the quadratic time difference forecast model variable estimated value originally stored on the star; for:

$\mathrm{<span\; class="notranslate">\; \&Delta;</span>}{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; \&Delta;</span>}{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}$

Among them: Δα _k is the forecast time difference at the required time t _k , Δα ₀ is the estimated time difference at the initial forecast time t ₀ , α ₁ is the frequency difference estimate of the on-board clock relative to the ground system clock, and α ₂ is the on-board clock Aging rate estimates relative to ground system clocks.

The calculation method of the phasing amount and the nominal phasing residual in the step (E) is: the phasing amount of the second pulse is: Δα _{k, remove} =-1/f ₀ round(Δα _k /(1/f ₀ ) ), the corresponding nominal phasing residual is: Δα _{k, res} = Δα _k + Δα _{k, remove} ;

Among them, round(·) indicates the rounding operation according to the rounding principle, and f ₀ is the fundamental frequency.

The method of judging whether the time difference forecast at the next moment in the step (H) is within the phase modulation interval is as follows: set the phase modulation index as _Δlim , and the phase modulation threshold as _Δgate , if the forecast time difference Δα _k at the next moment satisfies the following condition:

‖Δα _k ‖≤Δ _lim -Δ _gate

Then it is considered that this moment is the moment of autonomous phase modulation.

In the described step (1), the time difference processor utilizes the EKF filter to carry out the process of status update and time update to the time difference amount, frequency difference and aging rate of the secondary time difference prediction model:

(1) The EKF filter reads the current phase modulation identification, if the current phase modulation identification is "0", then enter step (3); if the current phase modulation identification is "1", then enter step (2);

(2) The EKF filter utilizes the nominal phase modulation residual to replace the prior estimate of the time difference amount given at the previous moment of the EKF filter, and then enters step (3);

(3) The time difference processor judges whether there is a time difference observation sequence of the local satellite time relative to the reference clock transmitted from the communication controller at the current moment, if the time difference observation sequence is entered into step (4), if there is no time difference observation sequence, then directly Go to step (5);

(4) Use the satellite-ground time difference observations given by the satellite-ground time difference observation to improve the prior estimation of the current time difference, frequency difference and aging rate, and then enter step (5);

(5) Using the dynamic change law of time difference, frequency difference and aging rate variables, predict the posterior estimated value of these variables, and give the prior estimate of the next moment;

(6) Finally, the time difference processor judges whether the current moment needs to use the state update and time update results of the EKF filter to replace the variable estimates in the quadratic time difference forecasting model according to the replacement period.

Compared with the prior art, the present invention has the advantages that: the present invention inherits the hardware device and main function division of the existing mature time management unit, and only needs to properly modify the software processing function of the time difference processor and the data entry and exit of the communication controller , can realize independent mathematical synchronization and physical synchronization operation, can greatly reduce the workload of the ground station, and ensure the ability of the system to work normally within a certain period of time when the ground station fails or is destroyed, especially suitable for the time-based task of the constellation system. This autonomous time system method endows the time management unit with long-term, high-precision punctuality and autonomous phase adjustment capabilities. In order to ensure the accuracy of the long-term time difference forecast on the satellite, the time difference processor uses the time difference between the on-board clock and the reference clock as much as possible. Observations, using an EKF filter to estimate the time difference, frequency difference and aging rate in real time, the latest filtering results of the time difference prediction model variables periodically replace the previously stored variable estimates, thus ensuring long-term and accurate tracking of the on-board time difference forecast Real time difference capability. The time difference processor judges the forecast time difference at the current moment. When the time difference enters the phase modulation interval, it is considered as the moment of autonomous phase modulation at the current moment. According to the principle of minimum phase modulation residual error, the phase modulation amount and the nominal phase modulation residual are calculated to drive the second pulse The generator performs corresponding phase adjustment according to the calculated phasing amount, and replaces the initial predicted time difference in the on-board time difference prediction model and the prior mean value estimate of the time difference component given by the EKF filter with the nominal phasing residual. The invention has carried out applied research and mathematical simulation in the satellite-earth time synchronization task of the 24-star Walker constellation. A total of two different satellite-ground time difference observation schemes are considered: one is to use the ground station and the inter-satellite link; the other is to use the GPS system to provide time service to the ground station and the satellites. The simulation results in these two different cases show that even if the on-board clock uses a typical high-temperature crystal oscillator, the autonomous timing method can strictly control the satellite-ground second pulse deviation within the required range of the index, which fully demonstrates the autonomous timing method. feasibility, effectiveness and superiority. The invention is also compatible with the autonomous time system method and the traditional ground station intervention time system method, so that the time system method provided by the invention has good reliability, safety and maneuverability.

Description of drawings

Fig. 1 is the structural representation of the existing on-board time management unit that the present invention adopts;

Fig. 2 is that the present invention utilizes the on-board time management unit to carry out the main flowchart of time system processing;

Fig. 3 is the flow chart that the present invention utilizes on-board time management unit to carry out EKF filtering;

Fig. 4 is the star-earth time difference estimation error result on a certain representative star in embodiment 1;

Fig. 5 is the time difference prediction error result on a certain representative star in embodiment 1;

Fig. 6 is the real second pulse phase deviation of the satellite on a certain representative star in embodiment 1: Fig. 6 A shows the trend of the phase deviation increasing with time when the phase modulation operation is not taken; Fig. 6 B represents the phase when taking the autonomous phase modulation operation The tendency of the deviation to increase over time;

Fig. 7 is the time difference estimation error result of the ground system clock relative to the navigation system clock that the EKF filter on the ground provides in embodiment 2;

Fig. 8 is the satellite-ground time difference estimation error result that the EKF filter on a certain representative star provides in embodiment 2;

Fig. 9 is the time difference prediction error result on a representative star in embodiment 2;

Fig. 10 is the real satellite-ground second pulse phase deviation of a certain representative star in embodiment 2; Fig. 10A shows the trend that the phase deviation increases at any time when the phase modulation operation is not taken; Fig. 10B shows that the phase deviation increases at any time when the phase modulation operation is taken the trend of.

Detailed ways

Below in conjunction with accompanying drawing and specific embodiment the present invention is described in further detail:

As shown in Figure 1, the existing time management unit used in the present invention is composed of five functional modules: a high-precision frequency source, a frequency synthesis splitter, a time difference processor, a second pulse generator and a communication controller.

High-precision frequency source: The stability of the frequency generated by the high-precision frequency source will affect the time uniformity of each satellite in the constellation, which will also affect the performance of the entire constellation. Therefore, the high-precision frequency source must be able to provide high-precision and stable frequencies as the reference frequency of the entire satellite.

Frequency synthesis splitter: The frequency synthesis splitter receives the reference frequency output by the high-precision frequency source, and decomposes it into multiple reference frequencies as needed to provide to the user of the high-precision time difference clock on the satellite, and provides 1 reference frequency for the second pulse generation The frequency synthesizer is used to generate multiple 1PPS reference pulses per second, and the frequency synthesis splitter should ensure the phase consistency between the decomposed multiple reference frequencies and the stability of each reference frequency.

Second pulse generator: The second pulse generator receives a reference frequency generated by the branch of the frequency synthesis splitter, and is responsible for converting it into multiple channels of 1PPS reference second pulse, and outputs it to the user of the high-precision time difference clock on the satellite. The second pulse generator can receive the second pulse phase adjustment command issued by the time difference processor, and adjust the phases of the output multi-channel reference second pulse according to the command. The second pulse generator shall ensure that the phase of the output reference second pulse is consistent with the input reference frequency.

Time difference processor: The time difference processor is the core module of the time management unit, which is responsible for solving and managing the time difference between the star time and the reference time. The time difference processor generates high-precision time at the beginning of the time management unit and transmits it to the time difference clock through the communication controller When the receiving device and the time difference clock user act as the local star, the high-precision time is generated by the crystal oscillator inside the time difference processor. The time difference processor receives the time difference information including the star time and the reference time from the time difference collection and sending device through the communication controller, and judges the validity of the received time difference information. For the median calculation, the initial time difference Δ ₀ is obtained by sorting the uploaded 5 time difference observations and then taking the median value. The time difference processor can calculate the standard time difference Δ according to the following formula:

Δ＝Δ ₀ +K(TT ₀ )

In the formula, T is the local satellite time at the moment of calculating the standard time difference Δ, and K is the frequency difference.

The time difference processor packs the current local star time plus the standard time difference Δ to form high-precision time information, and sends it to the time difference clock receiving device and the time difference clock user through the communication controller. The time difference processor receives the time difference clock adjustment command on the ground through the communication controller to adjust the standard time difference so that the local star time is adjusted to the reference time.

Communication controller: The communication controller is the interface module between the time management unit and the on-board bus, responsible for receiving the time difference transmitted from the on-board bus and the time difference clock adjustment command sent by the ground, and responsible for broadcasting the time difference processor calculated by the on-board bus The deviation in second pulses between local star time and the reference time.

The working process of the time management unit is: the time difference processor generates high-precision time at the beginning of the time management unit and transmits it to the time difference clock receiving device and the time difference clock user as the star through the communication controller. The crystal oscillator is generated, and the high-precision frequency source generates the reference frequency and sends it to the frequency synthesis splitter. The frequency synthesis splitter decomposes the reference frequency into n+1 reference frequencies, n≥1, and one of the reference frequencies is sent to the second pulse generator Generating n channels of reference second pulse signals, the time difference processor receives the time difference through the communication controller to collect the time difference information transmitted by the sending device, and the time difference processor calculates and processes the time difference information to generate the standard time difference between the star time and the reference time, time difference processing The device receives the time difference clock adjustment instructions on the ground through the communication controller to adjust the standard time difference, and the other n channels of reference frequency, the adjusted n channels of reference second pulse, standard time difference and high-precision time constitute the high-precision time information together to form the time difference The clock is transmitted to the time difference clock receiving device and the time difference clock user through the communication controller, so as to adjust the local satellite time to the reference time.

This invention adopts the same hardware configuration and main function division: the frequency source adopts a high-stable crystal oscillator or an atomic clock to provide a high-frequency (usually 10MHz) reference frequency signal as a time reference for on-board mission execution; frequency synthesis The splitter divides the base frequency signal into multiple channels for relevant users; the second pulse generator receives one reference frequency from the frequency synthesis splitter, and is responsible for converting it into multiple second pulse signals, which are output to the high-speed satellite. The precision time users use. After receiving the phase modulation command, the second pulse generator makes corresponding adjustments according to the phase adjustment amount; the communication controller is the interface module between the time management unit and the on-board bus, and is responsible for receiving various commands transmitted from the on-board bus (ground or Other stars send) time difference information, and are responsible for broadcasting the time difference information of this star calculated by the time difference processor to the bus on the star; the time difference processor is the core of the time management unit, and this invention further improves its main functions and software content Expansion, not only manages the generation and broadcasting of high-precision time information (combination of local star time and forecast time difference), but also has considerable command judgment and decision-making capabilities, and can complete autonomous timekeeping and phase modulation decision-making capabilities on the star. In addition, it also has With considerable data processing capabilities, it can complete the task of estimating the variables of the time difference forecast model.

The autonomous time system method of the present invention is as shown in Figure 2, and concrete realization steps are as follows:

(1) The time difference processor judges the time difference adjustment mark at the current moment, if the time difference adjustment mark is "0", then enters the autonomous time keeping step (2); if the time difference adjustment mark is "1", then enters the traditional ground timing processing flow ; If the time difference adjustment flag is "2", enter the traditional ground intervention centralized time adjustment process;

(2) The time difference processor judges the variable replacement flag of the on-board time difference prediction model. If the variable replacement flag is "1", the latest real-time estimation result of the on-board EKF filter is used to replace the corresponding variable of the current secondary time difference prediction model, and This moment is the initial prediction moment, and the variable replacement flag needs to be reset to "0", and then enter step (3); if the variable substitution flag is "0", then directly go to step (3);

(3) The time difference processor judges the ground intervention phasing flag, if the ground intervention phasing flag is "1", then the current phase modulation flag representing the phase modulation process at the current moment is set to "1", and the ground intervention phasing The flag is reset to "0", and then go to step (5); if the ground intervention phase modulation flag is "0", then go to step (4);

(4) The time difference processor judges the on-board autonomous phasing flag, if the on-satellite autonomous phasing flag is "0", then set the current phasing flag to "0", and then go to step (7); If the phase modulation flag is "1", then the current phase modulation flag is set to "1", and the autonomous phase modulation flag is set to "0", then go to step (5);

(5) According to the principle of minimum residual error after phase modulation, calculate the phase modulation amount and the nominal phase modulation residual error according to the current forecast time difference, and then turn to step (6);

(6) The time difference processor issues a second pulse phase modulation command to drive the second pulse generator to perform corresponding phase adjustment according to the calculated phase modulation amount, and unconditionally replace the initial time difference amount in the time difference prediction model with the nominal phase modulation residual, and use The current phasing moment is the initial forecast moment, then proceed to step (7);

(7) Based on the secondary time difference prediction model on the satellite, the satellite-ground time difference is predicted according to the required time on the satellite. On the one hand, the forecast time difference is used for broadcasting, and on the other hand, it is also used in the judging step (8) of the timing of autonomous phase adjustment;

(8) The time difference processor judges whether the time difference forecast at the next moment is in the phase modulation interval, if "Yes", then the self-phase modulation timing flag is set to "1", and then goes to step (9), otherwise, directly goes to Enter step (9);

(9) the time difference processor utilizes an EKF filter to carry out state update and time update processing to these variable estimates of the time difference amount, frequency difference and aging rate of the satellite-ground time difference prediction model, and then proceed to step (10);

(10) The time difference processor judges whether the current moment needs to utilize the result of the EKF filter to replace the relevant variable estimation in the quadratic time difference forecasting model according to the replacement period, if the current moment is the replacement moment, then the replacement flag is set to "1", Otherwise, set the replacement flag to "0", and then go to step (1).

The time difference adjustment mark in step (1) is a parameter stored on the satellite, and the time difference adjustment mark is modified by a time difference adjustment command sent from the ground. The content of the time difference adjustment in this invention is adjusted relative to the corresponding content of the patent application number 200810074975.3: the time difference processor of the present invention introduces an autonomous timekeeping method, and the method of ground intervention phase modulation is removed from these mathematical synchronization processing contents Removed, but also abandoned the uniform timing method. This invention defines the default value of the time difference adjustment flag as "0", indicating that the satellite adopts an autonomous timekeeping method; Just change it to "1" or "2". After the ground intervention is over, the flag should be reset to "0".

The time service or centralized time correction in step (1) is a mathematical synchronization method for changing the on-board time difference forecast. When the timing operation is performed, the time difference processor unconditionally accepts the satellite time information in the timing data as the local satellite time, and accepts the time difference in the timing data as the forecast time difference at the current moment. The effect of this operation on the time difference forecast is actually to change the initial time difference amount in the time difference forecast model. The on-board time difference prediction after the time service is completed by taking the time difference amount at the time service moment as the initial time difference amount, and the time difference prediction model composed of the original frequency difference and aging rate estimation. The time difference processor combines the time service satellite time and the forecast time difference, and broadcasts to relevant users through the communication controller. When the centralized time adjustment operation is performed, the time difference processor adjusts the star time at one time according to the star time adjustment amount in the centralized time adjustment data, and adjusts the forecast time difference calculated at the current time at one time according to the time difference adjustment amount in the centralized time adjustment data . The jet-difference forecast after the centralized time adjustment is completed by the jet-difference forecast model composed of the jet-difference amount at the time-adjusted moment as the initial jet-difference amount, and the original frequency difference and aging rate estimation. The time difference processor broadcasts the local satellite time and forecast time difference after centralized time correction to relevant users through the communication controller.

The on-board time difference forecast model variable replacement flag in step (2) is a parameter stored on the star, and the logo is modified by the replacement command sent by the time difference processor. If the flag is "0", it means that the real-time estimation result of the on-board EKF filter does not need to replace the original variable estimation of the quadratic time difference prediction model on the satellite; if the flag is "1", the replacement operation is performed.

The process of using the real-time estimation results of the on-board EKF filter to replace the original variable estimation process of the on-board quadratic time-difference prediction model is as follows: the time-difference processor bases the time-difference observation sequence of the local satellite time relative to the reference clock transmitted from the communication controller , using the EKF filter algorithm to give the latest estimated value of the quadratic time difference forecast model variable in real time, the time difference processor judges whether the current time is the replacement time according to the replacement period of the model variable, if it is the variable replacement time, use the latest estimated value to replace the on-board The previously stored quadratic time difference forecast model variable estimated value; otherwise, the next time difference forecast still adopts the previously stored quadratic time difference forecast model variable estimated value on the star. The quadratic time difference forecast model is:

$\mathrm{<span\; class="notranslate">\; \&Delta;</span>}{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; \&Delta;</span>}{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}<span\; class="notranslate">\; )</span>}^{\mathrm{<span\; class="notranslate">\; 2</span>}}$

Among them: Δα _k is the predicted time difference at the required time t _k , Δα ₀ is the estimated time difference at the initial forecast time t ₀ , α ₁ and α ₂ represent the frequency difference and aging rate estimation of the on-board clock relative to the ground system clock value.

The ground intervention phasing flag in step (3) is a parameter stored on the satellite, and the flag is modified by the ground intervention phasing command sent by the ground control unit. Ground intervention phasing is a physical synchronization operation separated from the original time difference adjustment content. This invention defines the default value of the ground intervention phasing flag as "0", which means the on-board autonomous phasing solution. When the ground decides to intervene During the phasing operation on the satellite, it is enough to issue an instruction to change the ground intervention phasing flag to "1". After the ground intervention is over, the flag will be restored to the default "0" value.

The on-board autonomous phasing flag and the current phasing flag in step (4) are two parameters stored on the star, both of which indicate that the current moment is the phasing moment, but the former indicates that the current moment is the autonomous phasing moment , and the latter may also indicate that the current moment is the moment of ground intervention phase modulation. The default values of the autonomous phasing flag and the current phasing flag are both "0", which means that the current moment is a non-phasing moment, and "1" means the phasing moment. After the relevant processing is completed, these two types of flags should be restored again to the default "0" value.

Calculation of phasing amount and nominal phasing residual in step (5). The phase modulation strategy is to make the nominal phase modulation residual error as small as possible. Since the phase modulation amount can only be an integer multiple of 1/f ₀ (f ₀ refers to the fundamental frequency), the phase modulation amount should be:

Δα _{k, remove} = -1/f ₀ round(Δα _k /(1/f ₀ ))

Among them, round(·) indicates the rounding operation according to the rounding principle. The corresponding nominal phasing residual is

Δα _{k, res} = Δα _k + Δα _{k, remove}

In step (7), in order to ensure long-term accurate forecasting of time difference, the time difference processor needs to be replaced periodically by the latest estimated value of the EKF filter. In addition, ground timing, centralized timing correction or phase modulation processing will change the time difference forecast after taking these measures, and the time difference processor also needs to automatically change some model variables and initial forecast time, and the relevant replacement methods have been mentioned above.

The next moment in step (8) is the judgment of the timing of the autonomous phase adjustment. At the current moment, the time difference processor judges whether the predicted time difference at the next time has entered the phase modulation interval, and if it has entered the interval, it will set the autonomous phase modulation flag to "1". The phase modulation interval is jointly determined by the phase modulation index _Δlim and the phase modulation threshold _Δgate . The judgment conditions for an autonomous phasing timing are described as follows:

‖Δα _k ‖≤Δ _lim -Δ _gate

If the above conditions are met, it is considered that the time for autonomous phase adjustment has arrived. Among them, _Δlim is the system time synchronization index, and the initial phase modulation threshold _Δgate should be selected according to the time difference estimation accuracy, but in order to avoid misjudgment, it should be finally determined by simulation.

The EKF filtering process of the time difference forecast model variable in step (9) is shown in Figure 3. Compared with other general EKF filtering processes, step (A) is added. This step is to eliminate the adaptation of the influence of phase modulation processing on the estimation of EKF filter sexual countermeasures. Due to the phase modulation, the real time difference has changed greatly instantaneously, and the time difference component estimate given by the general time difference EKF filter cannot predict this large change. Therefore, the EKF filter may cause the filter to Entering a longer perturbation state may even cause the filter to diverge (expressed as inconsistent estimation). Since the nominal phase modulation residual is known, using this known quantity to replace the prior mean value estimate of the time difference component given by the EKF filter at the previous moment can ensure the smooth working performance of the filter. After that, enter the general EKF filtering process. Including an observation update step: using possible satellite-ground time difference observation opportunities (for example, satellite-ground time difference observation based on ground stations, satellite-ground time difference observations based on ground stations and inter-satellite links, satellite-ground time difference observations based on navigation systems) to To improve the a priori estimates of the current time difference, frequency difference and aging rate, and a time update step: using the dynamic change law of the time difference, frequency difference and aging rate variables, the posterior estimation of these variables The prior estimate is used to predict the value and give the prior estimate of the next moment. The specific mathematical description is as follows:

The state quantity is defined as

x _k = [Δα _k α _{1, k} α _{2, k} ] ^T

Among them, the first component represents the amount of time difference at the current moment; the latter two components represent the frequency difference and the aging rate respectively. The general EKF filter requires a filter initialization step, that is, it needs to be given an initial prior state mean estimate and covariance estimation P _0|-1 , as well as system noise variance Q _k and observation noise variance R _k , then the following recursive estimation process can be started:

● State estimation observation update (with observation time):

${\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; =</span>{\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}{\mathrm{<span\; class="notranslate">\; \&upsi;</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}$

${\mathrm{<span\; class="notranslate">\; P</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; P</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}{\mathrm{<span\; class="notranslate">\; P</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}^{\mathrm{<span\; class="notranslate">\; \&upsi;\&upsi;</span>}}{\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}$

Among them, the new information and its variance matrix are

${\mathrm{<span\; class="notranslate">\; \&upsi;</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; h</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}{\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}$

${\mathrm{<span\; class="notranslate">\; P</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}^{\mathrm{<span\; class="notranslate">\; \&upsi;\&upsi;</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; h</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}{\mathrm{<span\; class="notranslate">\; P</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; h</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; +</span>{\mathrm{<span\; class="notranslate">\; R</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}$

The covariance matrix is

${\mathrm{<span\; class="notranslate">\; P</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}^{\mathrm{<span\; class="notranslate">\; xz</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; P</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\mathrm{<span\; class="notranslate">\; h</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}$

The gain matrix is

${\mathrm{<span\; class="notranslate">\; K</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; P</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}^{\mathrm{<span\; class="notranslate">\; xz</span>}}{<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; P</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}^{\mathrm{<span\; class="notranslate">\; \&upsi;\&upsi;</span>}}<span\; class="notranslate">\; )</span>}^{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}$

● State estimation time update:

${\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; k</span>}}{\stackrel{<span\; class="notranslate">\; ~</span>}{\mathrm{<span\; class="notranslate">\; x</span>}}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; k</span>}}$

${\mathrm{<span\; class="notranslate">\; P</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; k</span>}}{\mathrm{<span\; class="notranslate">\; P</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; |</span>\mathrm{<span\; class="notranslate">\; k</span>}}{\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; k</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; \&Gamma;</span>}{\mathrm{<span\; class="notranslate">\; Q</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}{\mathrm{<span\; class="notranslate">\; \&Gamma;</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}$

The specific description of several parameters in the above general time difference EKF filter formula is as follows: the state transition matrix is

${\mathrm{<span\; class="notranslate">\; \&Phi;</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; =</span>\left[\begin{array}{ccc}\mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; \&Delta;\&tau;</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; /</span>\mathrm{<span\; class="notranslate">\; 2</span>}\mathrm{<span\; class="notranslate">\; \&Delta;</span>}{\mathrm{<span\; class="notranslate">\; \&tau;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}\\ \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; \&Delta;\&tau;</span>}\\ \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}\end{array}\right]$

Δτ=t _k -t _k-1 represents a discrete period;

The discrete system noise w _k ＝[ε _k η _k ξ _k ] ^T is determined by the performance of the clock product itself, and the process noise variance matrix Q _k is a quantity related to Δτ:

${\mathrm{<span\; class="notranslate">\; Q</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; =</span>\left[\begin{array}{ccc}{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; \&alpha;\&alpha;</span>}}& {\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; \&alpha;\&alpha;</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}& {\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; \&alpha;\&alpha;</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}\\ {\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; \&alpha;\&alpha;</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}& {\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}\mathrm{<span\; class="notranslate">\; 1</span>}\mathrm{<span\; class="notranslate">\; \&alpha;</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}& {\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}\mathrm{<span\; class="notranslate">\; 1</span>}\mathrm{<span\; class="notranslate">\; \&alpha;</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}\\ {\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; \&alpha;\&alpha;</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}& {\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}\mathrm{<span\; class="notranslate">\; 1</span>}\mathrm{<span\; class="notranslate">\; \&alpha;</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}& {\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}\mathrm{<span\; class="notranslate">\; 2</span>}\mathrm{<span\; class="notranslate">\; \&alpha;</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}\end{array}\right]$

Its components are expressed as follows:

${\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; \&alpha;\&alpha;</span>}}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; <figure-callout\; id="0"\; label="h"\; filenames="H2008102251757E00021.png,H2008102251757E00031.png"\; state="\{\{state\}\}">h</figure-callout></span>}}_{\mathrm{<span\; class="notranslate">\; 0</span>}}}{\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; \&Delta;\&tau;</span>}<span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; \&pi;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; h</span>}}_{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 2</span>}}}{\mathrm{<span\; class="notranslate">\; 3</span>}}\mathrm{<span\; class="notranslate">\; \&Delta;</span>}{\mathrm{<span\; class="notranslate">\; \&tau;</span>}}^{\mathrm{<span\; class="notranslate">\; 3</span>}}<span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; \&pi;</span>}}^{\mathrm{<span\; class="notranslate">\; 4</span>}}{\mathrm{<span\; class="notranslate">\; h</span>}}_{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 4</span>}}}{\mathrm{<span\; class="notranslate">\; 5</span>}}\mathrm{<span\; class="notranslate">\; \&Delta;</span>}{\mathrm{<span\; class="notranslate">\; \&tau;</span>}}^{\mathrm{<span\; class="notranslate">\; 5</span>}}$

${\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; \&alpha;</span>}\mathrm{<span\; class="notranslate">\; 1</span>}\mathrm{<span\; class="notranslate">\; \&alpha;</span>}\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; \&pi;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; h</span>}}_{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 2</span>}}\mathrm{<span\; class="notranslate">\; \&Delta;\&tau;</span>}<span\; class="notranslate">\; +</span>\frac{\mathrm{<span\; class="notranslate">\; 8</span>}{\mathrm{<span\; class="notranslate">\; \&pi;</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; h</span>}}_{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 4</span>}}}{\mathrm{<span\; class="notranslate">\; 3</span>}}\mathrm{<span\; class="notranslate">\; \&Delta;</span>}{\mathrm{<span\; class="notranslate">\; \&tau;</span>}}^{\mathrm{<span\; class="notranslate">\; 3</span>}}$

q _α2α2 = 8π ⁴ h _-4 Δτ

q _αα1 = π ² h _-2 Δτ ² + π ⁴ h _-4 Δτ ⁴

${\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; \&alpha;\&alpha;</span>}\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; =</span>\frac{\mathrm{<span\; class="notranslate">\; 4</span>}{\mathrm{<span\; class="notranslate">\; \&pi;</span>}}^{\mathrm{<span\; class="notranslate">\; 4</span>}}{\mathrm{<span\; class="notranslate">\; h</span>}}_{<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 4</span>}}\mathrm{<span\; class="notranslate">\; \&Delta;</span>}{\mathrm{<span\; class="notranslate">\; \&tau;</span>}}^{\mathrm{<span\; class="notranslate">\; 3</span>}}}{\mathrm{<span\; class="notranslate">\; 3</span>}}$

q _α1α2 = 4π ⁴ h _-4 Δτ ²

h _a in the above formula represents the noise intensity coefficient of the energy spectrum whose noise index is a; the system noise distribution matrix is

$\mathrm{<span\; class="notranslate">\; \&Gamma;</span>}<span\; class="notranslate">\; =</span>\left[\begin{array}{ccc}\mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}\\ \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}\\ \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 0</span>}& \mathrm{<span\; class="notranslate">\; 1</span>}\end{array}\right]$

The observation matrix is

${\mathrm{<span\; class="notranslate">\; h</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; =</span>\frac{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; h</span>}}{<span\; class="notranslate">\; \&PartialD;</span>\mathrm{<span\; class="notranslate">\; x</span>}}{<span\; class="notranslate">\; |</span>}_{\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; x</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}}$

The variance matrix of the observation noise v _k is R _k .

Selection of the model variable replacement period in step (10). When the variable replacement period is equal to the observation update period of the EKF filter (that is, the state posterior estimation at each observation update moment is used to replace the corresponding variables of the quadratic time difference forecast model), the forecast accuracy of the quadratic time difference forecast model is equal to that of the EKF filter The prior estimation accuracy of the time difference component of the sensor. If the variable replacement period is greater than the observation update period, it will lead to a decline in the accuracy of the quadratic time difference forecast model. Considering the difficulty of engineering realization, the two cycles can be different, of course, should be properly selected according to the simulation.

Example 1

Use the ground station and the inter-satellite link to realize the autonomous synchronization processing between the clocks on each satellite and the ground system clock in the satellite constellation (24-star Walker constellation) (for example, the goal is to control the second pulse synchronization error between the clocks on each satellite and the ground system clock at all times In the range of 1us.) The clock on the star uses a typical high-stable crystal oscillator QZ clock, and the ground uses a CS clock. In each time period (assuming 5s), the ground station measures the time difference between the transiting satellite clock and the ground system clock (assuming the accuracy is 30ns), and the inter-satellite time difference is measured between adjacent satellites through the inter-satellite link (assuming the accuracy for 10ns). Based on the satellite-ground time difference observation data of the uploaded transit satellites, combined with the inter-satellite time measurement data, the transit star is used as the starting star, and the satellite-ground time difference observations of each star relative to the ground station are converted according to the principle of the shortest satellite-ground path. The processor uses the on-board EKF filter to estimate the satellite-ground time difference prediction model variables. If the current moment is the time to replace the variables of the secondary time difference forecast model (assuming the replacement period is 10s), the time difference processor will issue a model variable replacement command to replace the original variable estimates with the latest estimated results given by the filter. Time management The unit performs time difference forecast based on the new quadratic time difference forecast model, and uses the forecast time difference for broadcasting and judging the timing of autonomous phasing. When autonomous phasing is required, calculate the amount of phasing and nominal phasing residual, and issue instructions The driving second pulse generator is adjusted accordingly according to the calculated phasing amount. The specific description of the above-mentioned on-board autonomous time system realization process is as follows:

(1) the time difference processor judges the time difference adjustment sign of the current moment, and finds that it is "0", then enters the autonomous time-keeping step (2);

(2) Determine the variable replacement flag of the on-board time difference forecast model, if it is "1", use the latest real-time estimation result of the on-board EKF filter to replace the corresponding variables of the current secondary time difference forecast model, and use this moment as the starting forecast At the moment, the variable replacement flag is reset to "0", and then enters step (3); if the flag is "0", then directly proceeds to step (3);

(3) The time difference processor judges the ground intervention phase modulation flag, finds that it is "0", and directly transfers to step (4);

(4) The time difference processor judges the autonomous phasing flag on the star. If it is "0", set the current phasing flag to "0", and then go to step (7); if it is "1", set the current phasing flag to "0". The phase flag is set to "1", and the autonomous phase modulation flag is set to "0", and then go to step (5);

(5) According to the principle of minimum residual error after phase modulation, calculate the phase modulation amount and the nominal phase modulation residual error according to the current forecast time difference, and then turn to step (6);

(6) The time difference processor issues a second pulse phase modulation command to drive the second pulse generator to perform corresponding phase adjustment according to the calculated phase modulation amount, and unconditionally replace the initial time difference amount in the time difference prediction model with the nominal phase modulation residual, and use The current phasing moment is the initial forecast moment, then proceed to step (7);

(7) Based on the secondary time difference prediction model on the satellite, the satellite-ground time difference is forecasted according to the needs on the satellite. On the one hand, the forecast time difference is used for broadcasting, and on the other hand, it is also used in the judging step (8) of the timing of autonomous phase adjustment;

(8) The time difference processor judges whether the time difference forecast at the next moment is in the phase modulation interval, if "Yes", then the self-phase modulation timing flag is set to "1", and then goes to step (9), otherwise, directly goes to Enter step (9);

(9) The time difference processor uses the EKF filter to perform status update and time update processing on the estimates of variables such as time difference, frequency difference and aging rate of the quadratic time difference forecast model. If the current phase modulation flag is "1", an additional important operation needs to be performed first, that is, the calculated nominal phase modulation residual is used to replace the prior estimate of the mean value of the time difference component. Then enter the usual EKF filtering process. If there is satellite-earth time difference observation data (after conversion) at the current moment, then the prior state estimation is updated by using the observation data, and the posteriori state estimation result is given. Regardless of whether there is time difference observation data, the conventional time update process is performed, and the prior state estimation of the next moment is given. Go to step (10) at last;

(10) The time difference processor judges whether the relevant variables in the quadratic time difference forecasting model need to be replaced by the results of the EKF filter at the current moment according to the replacement period. If the current moment is the replacement moment, then the replacement identifier is set to "1", otherwise the replacement identifier is set to "0", then proceed to step (1);

A star is selected as a representative star, and a set of mathematical simulation results after the above-mentioned time system processing are shown in Figures 4, 5, and 6. Fig. 4 is the estimated error result of the satellite-ground time difference given by the EKF filter on a certain representative star in embodiment 1, the result shows that the EKF filter estimate is a consistent estimate, and the phase modulation process disturbs the EKF filter on the star The influence is also well suppressed, and the steady-state estimation accuracy of the time difference is better than 10 ns; Fig. 5 is the result of the time difference prediction error of the quadratic time difference prediction model on a certain representative star in Example 1, and the result shows After the estimated result of the above EKF filter replaces the original model variable estimation, the predicted time difference generally tracks the time difference estimation result of the EKF filter on the star; Deviation: Figure 6A shows that the phase deviation increases at any time when the phase modulation operation is not adopted, and soon exceeds the index requirement of 1us; The expected synchronization indicator requirements are met.

Example 2

Use the GPS navigation system to realize the autonomous synchronization processing between the clocks on each satellite constellation (24-star Walker constellation) and the ground system clock (for example, the second pulse synchronization error between the clocks on each satellite and the ground system clock is always controlled within 1us.) A typical QZ clock is used as the clock on the star, and a CS clock is used on the ground. In each time period (assuming 5s), the ground station uses the GPS navigation system to serve time, and measures the system time difference observation sequence of the ground system clock relative to the navigation system clock (assuming that the timing accuracy is 100ns, but the actual accuracy is higher than this!), A ground EKF filter similar to the on-board EKF filter, based on this type of observation data, gives a posteriori estimates of the time difference, frequency difference and aging rate of the system clock, and passes through the satellite-ground link at regular intervals (assuming 1000s) The road and the inter-satellite link are transmitted to each on-board time management unit. In addition, the ground does not intervene in the on-board time system processing. Through the GPS navigation system timing on the star, the time difference observation sequence of the on-board clock relative to the navigation system clock is measured (assuming that the timing accuracy is 100ns), and the time difference prediction model variable of the ground system clock relative to the navigation system clock stored on the star is uploaded periodically The system time difference at the current moment is predicted, and the two types of time difference data are converted into the observed amount of the satellite-ground time difference at the current moment, and the EKF filter on the satellite estimates the satellite-ground time difference prediction model variables accordingly. If the current moment is the time to replace the variable estimation of the secondary time difference forecast model on the satellite (assuming that the replacement cycle is 10s), the time difference processor issues a model variable replacement command to replace the original estimated result with the latest estimated result given by the EKF filter on the satellite. Variable estimation, the time management unit forecasts the time difference based on the new quadratic time difference prediction model, and uses it to judge the timing of broadcasting and autonomous phasing. The instruction drives the second pulse generator to make corresponding adjustments according to the calculated phasing amount. The specific description of the above-mentioned on-board autonomous time system realization process is as follows:

(1) the time difference processor judges the time difference adjustment sign of the current moment, and finds that it is "0", then enters the autonomous time-keeping step (2);

(2) Determine the variable replacement flag of the on-board time difference forecast model, if it is "1", use the latest real-time estimation result of the on-board EKF filter to replace the corresponding variables of the current secondary time difference forecast model, and use this moment as the starting forecast At this time, it is also necessary to reset the variable replacement flag to "0", and then enter step (3); if the flag is "0", then directly go to step (3);

(3) The time difference processor judges the ground intervention phase modulation flag, finds that it is "0", and turns to step (4);

(4) The time difference processor judges the autonomous phasing flag on the star. If it is "0", set the current phasing flag to "0", and then go to step (7); if it is "1", set the current phasing flag to "0". The phase flag is set to "1", and the autonomous phase modulation flag is set to "0", and then go to step (5);

(5) According to the principle of minimum residual error after phase modulation, calculate the phase modulation amount and the nominal phase modulation residual error according to the current forecast time difference, and then turn to step (6);

(6) The time difference processor issues a second pulse phase modulation command to drive the second pulse generator to perform corresponding phase adjustment according to the calculated phase modulation amount, and unconditionally replace the initial time difference amount in the time difference prediction model with the nominal phase modulation residual, and use The current phasing moment is the initial forecast moment, then proceed to step (7);

(7) Based on the secondary time difference prediction model on the satellite, the satellite-ground time difference is forecasted according to the needs on the satellite. On the one hand, the forecast time difference is used for broadcasting, and on the other hand, it is also used in the judging step (8) of the timing of autonomous phase adjustment;

(8) The time difference processor judges whether the time difference forecast at the next moment is in the phase modulation interval, if "Yes", then the self-phase modulation timing flag is set to "1", and then goes to step (9), otherwise, directly goes to Enter step (9);

(9) The time difference processor uses the EKF filter to perform status update and time update processing on the estimates of variables such as time difference, frequency difference and aging rate of the quadratic time difference forecast model. If the current phase modulation flag is "1", an additional important operation needs to be performed first, that is, the calculated nominal phase modulation residual is used to replace the prior estimate of the mean value of the time difference component. Then enter the usual EKF filtering process. If there is satellite-earth time difference observation data (after conversion) at the current moment, then the prior state estimation is updated by using the observation data, and the posteriori state estimation result is given. Regardless of whether there is time difference observation data, the conventional time update process is performed, and the prior state estimation of the next moment is given. Go to step (10) at last;

(10) The time difference processor judges whether the relevant variables in the quadratic time difference forecasting model need to be replaced by the results of the EKF filter at the current moment according to the replacement cycle. If the current moment is the replacement moment, then the replacement identifier is set to "1", otherwise the replacement identifier is set to "0", then proceed to step (1);

Choose a star as a representative star, and a group of mathematical simulation results after the above-mentioned time system processing are shown in Figures 7, 8, 9, and 10: Figure 7 is the ground system provided by the ground EKF filter in Embodiment 2 The time difference estimation error result of the clock relative to the navigation system clock, the result shows that the EKF filter estimation is a consistent estimation, and its steady-state estimation accuracy reaches 1 ns; Fig. 8 is given by the EKF filter on a certain representative star in embodiment 2 The estimation error results of the satellite-ground time difference, the results show that the EKF filter estimation is a consistent estimation, and the disturbance effect of the phase modulation process on the EKF filter on the satellite is also well suppressed, and its steady-state estimation accuracy is about 10 ns; Fig. 9 It is the time difference prediction error result of the secondary time difference prediction model on a representative star in Example 2. The result shows that the estimated results of the EKF filter on the satellite are used to replace the original model variables at intervals of one time period, and the forecast time difference is roughly tracked The time difference estimation result of the EKF filter on the star; Fig. 10 is the real satellite-ground second pulse phase deviation of a certain representative star in embodiment 2: Fig. 10A shows that the phase deviation increases at any time when the phase modulation operation is not taken, and soon The index requirement of 1us is exceeded; Fig. 10B shows that the phase deviation after adopting the autonomous timing system has been controlled within the range of 1us, which meets the expected synchronization index requirement.

It is not difficult to see that the content of autonomous time statistics in Embodiment 1 and Embodiment 2 is almost the same. This shows that the present invention has good versatility and adaptability. The simulation result of embodiment 1 and embodiment 2 illustrates, under the situation that adopts QZ clock on the star, this autonomous time system method still strictly controls the satellite-ground second pulse phase difference in the scope of index requirement (absolute value is no more than 1us) , which cannot be achieved by traditional methods, which fully demonstrates the feasibility, effectiveness and superiority of this method.

The content not described in detail in the present invention is well known to those skilled in the art.

Claims (7)

1. A kind of autonomous time system method on the star, it is characterized in that comprising the following steps: (A) the time difference processor judges the time difference adjustment mark of the current moment, if the time difference adjustment mark is "0", then enters step (B) and carries out self-timekeeping; if the time difference adjustment mark is "1", then enters the ground timing processing flow; If the time difference adjustment flag is "2", enter the process of centralized time adjustment on the ground; (B) The time difference processor judges the on-board quadratic time-difference forecast model variable replacement flag, if the variable substitution flag is "1", then use the latest real-time estimated results of the on-board EKF filter to replace the corresponding variables of the current quadratic time-difference prediction model, And take this current moment as the initial forecasting time, and reset the variable replacement sign of the secondary time difference forecast model to "0", and then enter step (C); if the variable replacement sign is "0", then directly go to step ( C); (C) The time difference processor judges the ground intervention phasing flag, if the ground intervention phasing flag is "1", then the current phase modulation flag representing the phase modulation process at the current moment is set to "1", and the ground intervention phasing The flag is reset to "0", and then go to step (E); if the ground intervention phase modulation flag is "0", go to step (D); (D) The time difference processor judges the on-board autonomous phasing flag, if the on-satellite autonomous phasing flag is "0", then set the current phasing flag to "0", and then go to step (G); If the phase modulation flag is "1", set the current phase modulation flag to "1", reset the autonomous phase modulation flag to "0", and then go to step (E); (E) Calculate the phasing amount and the nominal phasing residual according to the current forecast satellite-ground time difference, and then turn to step (F); (F) The time-of-flight processor issues a second-pulse phase modulation command to drive the second-pulse generator to perform corresponding phase adjustments according to the calculated phase-modulation amount, and unconditionally replace the initial time-difference amount in the secondary time-difference prediction model with the nominal phase-modulation residual , and take the current phasing moment as the initial forecasting moment, and then turn to step (G); (G) Based on the secondary time difference prediction model on the satellite, the satellite-ground time difference is forecasted according to the required time on the star, and the satellite-ground time difference is broadcast to the user on the one hand, and is used for the judging step (H) of the timing of the autonomous phase adjustment on the other hand; (H) the time difference processor judges whether the satellite-earth time difference forecast of the next moment is in the phase modulation interval, if "Yes", then the self-phase modulation timing mark is set to "1", then proceeds to step (1), otherwise, Go directly to step (I); (1) the time difference processor utilizes the EKF filter to carry out state update and time update to the time difference amount, frequency difference and aging rate of the secondary time difference forecasting model, then proceed to step (J); (J) The time difference processor judges whether the current moment needs to utilize the state update and the time update result of the EKF filter to replace the variable estimation in the secondary time difference forecasting model according to the replacement period, if the current moment is the replacement time, then the replacement sign is set to "1", otherwise set the replacement flag to "0", and then go to step (A). 2. The autonomous time system method on a kind of star according to claim 1, is characterized in that: the processing flow of time service in the described step (A) is: when carrying out time service operation, time difference processor unconditionally will be in the time service data The stellar time information is accepted as the star time, and the time difference in the timing data is accepted as the forecast star-ground time difference at the current moment. 3. The autonomous time system method on a kind of star according to claim 1, it is characterized in that: the processing flow of centralized time correction in the described step (A) is: when carrying out centralized time correction operation, time difference processor one-time Adjust the local star time according to the star time adjustment amount in the centralized time correction data, and adjust the forecast star-earth time difference calculated at the current moment at one time according to the time difference adjustment amount in the centralized time correction data. 4. a kind of on-star autonomous time statistics method according to claim 1 is characterized in that: in the described step (B), the time difference processor utilizes the latest real-time estimation result of the EKF filter on the star to replace the current secondary time difference prediction model The realization process of the corresponding variable is as follows: the time difference processor uses the EKF filter algorithm to give the latest estimated value of the secondary time difference forecast model variables in real time according to the time difference observation sequence of the local satellite time relative to the reference clock transmitted from the communication controller, and the time difference The processor judges whether the current time is the replacement time according to the replacement period of the model variable. If it is the variable replacement time, the latest estimated value is used to replace the estimated value of the secondary time difference forecast model variable previously stored on the star; otherwise, the next satellite-ground time difference forecast Still adopting the variable estimated value of the quadratic time difference forecast model originally stored on the star, the described quadratic time difference forecast model is:

Among them: Δα _k is the predicted satellite-ground time difference at the required time t _k , Δα ₀ is the time difference estimate at the initial forecast time t ₀ , α ₁ is the frequency difference estimate of the on-board clock relative to the ground system clock, and α ₂ is the satellite-ground time difference estimate. Estimation of the aging rate of the upper clock relative to the ground system clock. 5. a kind of on-satellite autonomous time statistics method according to claim 1 is characterized in that: in the described step (E), the calculation method of phasing amount and nominal phasing residual is: the phasing amount of second pulse is : Δα _{k, remove} = -1/f ₀ round(Δα _k /(1/f ₀ )), the corresponding nominal phase modulation residual is: Δα _{k, res} = Δα _k + Δα _{k, remove} ; Among them, Δα _k is the predicted satellite-ground time difference at the required time t _k , round(·) means rounding operation according to the rounding principle, and f ₀ is the fundamental frequency. 6. The autonomous time statistics method on a kind of star according to claim 1, it is characterized in that: in the described step (H), whether the satellite-earth time difference forecast of next moment is in the judging method in the phase modulation interval is: set modulation The phase index is Δ _lim , and the phase modulation threshold is Δ _gate . If the predicted satellite-ground time difference Δα _k at the next moment satisfies the following conditions: ||Δα _k ||≤Δ _lim -Δ _gate Then it is considered that this moment is the moment of autonomous phase modulation. 7. a kind of autonomous time statistics method on the star according to claim 1, is characterized in that: in described step (1), time difference processor utilizes EKF filter to the time difference amount, frequency difference and aging of secondary time difference prediction model The process of status update and time update is as follows: (1) The EKF filter reads the current phase modulation identification, if the current phase modulation identification is "0", then enter step (3); if the current phase modulation identification is "1", then enter step (2); (2) The EKF filter utilizes the nominal phase modulation residual to replace the prior estimate of the time difference amount given at the previous moment of the EKF filter, and then enters step (3); (3) The time difference processor judges whether there is a time difference observation sequence of the local satellite time relative to the reference clock transmitted from the communication controller at the current moment, if the time difference observation sequence is entered into step (4), if there is no time difference observation sequence, then directly Go to step (5); (4) Use the satellite-ground time difference observations given by the satellite-ground time difference observation to improve the prior estimation of the current time difference, frequency difference and aging rate, and then enter step (5); (5) Using the dynamic change law of time difference, frequency difference and aging rate variables, predict the posterior estimated value of these variables, and give the prior estimate of the next moment. the

Images