A time synchronization method based on BeiDou satellite

Abstract

This invention discloses a time synchronization method based on BeiDou satellites, comprising the following steps: S1: acquiring the ephemeris, satellite observation data, and B2b satellite correction data messages of BeiDou-3 satellites; S2: performing anomaly processing on the satellite observation data; S3: combining the ephemeris and the anomaly-processed satellite observation data, using a PPP floating-point solution function model for pseudorange single-point positioning; S4: calculating the satellite's position coordinates based on the ephemeris and B2b satellite correction data messages; S5: inputting the pseudorange single-point positioning result and the satellite's position coordinates into a Kalman filter for filtering, followed by PPP precise single-point positioning processing to obtain the receiver clock bias; S6: inputting the receiver clock bias data into a clock control module to adjust the clock and complete synchronization. In this invention, each device acts as a master station, achieving high-precision positioning of each device through the PPP-B2b model, combining positioning technology with time synchronization technology, thereby improving the synchronization accuracy between devices.

Description

Technical Field

[0001] This invention belongs to the field of satellite time synchronization technology, specifically relating to a time synchronization method based on the BeiDou satellite. Background Technology

[0002] Since the completion and global service launch of the BeiDou-3 satellite system in August 2020, it has provided global positioning, navigation, and timing (PNT) services, short message communication, international search and rescue, satellite-based augmentation system (SBAS), and precision point positioning (PPP) services. The PPP service is achieved through the PPP-B2b signal broadcast by three GEO satellites, and currently relies on the precise correction products of 26 BDS3 satellites broadcasting the PPP-B2b signal. This enables real-time dynamic decimeter-level and static centimeter-level positioning accuracy.

[0003] However, the commonly used technology at present is RTK technology to achieve high-precision time synchronization between devices. But using RTK technology requires setting up an observation base station, and then the observation base station transmits the time deviation to each sub-station. The transmitted time deviation is used to compensate for the phase of the sub-station, thereby ensuring time synchronization between the sub-stations. If the observation base station broadcasts abnormally, the time synchronization of the sub-stations cannot be guaranteed.

[0004] To address the aforementioned issues, this system proposes a high-precision time synchronization technology based on PPP-B2b. This technology utilizes the PPP-B2b service signal broadcast by the BeiDou-3 satellite to achieve time synchronization between devices using PPP-B2b precise single-point positioning GNSS. This method can achieve synchronization accuracy between devices at the 3ns level. This technology meets the dual requirements of time synchronization accuracy and flexibility in the national defense and military fields, and also provides a precise means for cross-regional precise time synchronization, rapid backup of time scales, and long baseline measurement. Summary of the Invention

[0005] To address the problems mentioned in the background section, this invention provides a time synchronization method based on the BeiDou satellite system, which solves the problem in existing traditional RTK models that require the establishment of fixed base stations and cannot guarantee time synchronization when observing abnormal broadcasts from the base stations.

[0006] To achieve the above objectives, the present invention provides the following technical solution:

[0007] A time synchronization method based on BeiDou satellites includes the following steps:

[0008] S1: Obtain the ephemeris, satellite observation data, and B2b satellite correction data messages of BeiDou-3 satellites through the receiver;

[0009] S2: Perform anomaly handling, which includes a preprocessing algorithm and a cycle slip detection algorithm. The preprocessing algorithm is used to remove satellite observation data that is below a predetermined value, and the cycle slip detection algorithm is used to remove abnormal clock error data in the satellite observation data.

[0010] S3: Combining ephemeris and anomaly-processed satellite observation data, pseudorange single-point positioning is performed using the PPP floating-point solution function model;

[0011] S4: Based on the ephemeris and B2b satellite correction data, determine the position and clock bias of the B2b satellite, recover the precise ephemeris of the B2b satellite through the clock bias correction algorithm, and calculate the position coordinates of the satellite.

[0012] S5: Input the pseudorange point positioning result and the satellite's position coordinates into a Kalman filter for filtering. Then, perform PPP precise point positioning on the filtered data to calculate the receiver's position and clock error.

[0013] S6: Input the receiver clock difference data obtained in S5 into the clock control module to periodically adjust the frequency and phase of the clock to complete time synchronization.

[0014] Preferably, the preprocessing algorithm in S2 is as follows:

[0015] Suppose the satellite number received by the receiver is , ,…… The corresponding satellite signal-to-noise ratio is , ,…… The corresponding satellite elevation angle is , ,…… The preprocessing method involves eliminating satellites that do not meet the requirements based on the signal-to-noise ratio and elevation angle. The specific formula is as follows:

[0016] ;

[0017] ;

[0018] The sum of the two formulas above gives the calculated satellite mass value, i.e.:

[0019] ;

[0020] If the current satellite number Calculation results If the value is less than 60, the satellite is marked as unhealthy and removed from the list.

[0021] Preferably, the cycle slip detection algorithm in S2 is a phase carrier wide-lane combination of combined observations, and the observation equation is:

[0022] ;

[0023] ;

[0024] ;

[0025] in , Indicates the current frequency. , This represents the corresponding carrier phase observation value. Indicates the time of data collection. This indicates the cycle slip detection threshold, and E represents the satellite elevation angle. Indicates the wavelength of the wide lane;

[0026] pass If the collected observation value is found to be within the cycle slip detection threshold through iterative calculation, the data is retained for subsequent calculations.

[0027] Preferably, in S3, the PPP floating-point solution function model simultaneously uses approximate range and carrier phase observations for parameter estimation, and constructs observation equations to describe the functional relationship between the observations and various estimated parameters and correction errors. Let the frequency be... The pseudorange and carrier phase observation equations are as follows:

[0028] Formula A: ;

[0029] Among them, superscript Indicates satellite, subscript Indicates receiver, subscript Indicates the frequency band of the observed value. Represents pseudorange observations, in units of This represents the carrier phase observation value, in cycles. This represents the geometric distance between the phase centers of the receiver antenna and the satellite antenna, in units of... and Clock bias indicates the time difference between the receiver and satellite clocks relative to a time reference, and is measured in units of 1. The unit representing the speed of light in a vacuum is _____. and The unit representing the frequency-dependent pseudorange hardware delay at the receiver and satellite ends is . This represents the frequency-dependent hardware delay of the receiver and satellite carrier phase, measured in cycles. This represents the frequency-dependent ionospheric delay, measured in meters (m). This represents the tropospheric delayed wet component mapping function related to the satellite elevation angle. The tropospheric delay moisture component in the zenith direction of the observation station is shown in meters. Represents frequency The corresponding carrier wavelength, in meters (m). This represents the carrier phase integer ambiguity, measured in cycles. and The unmodeling error noise represents the pseudorange and carrier phase observations, in units of ;

[0030] The formula for calculating the first-order ionospheric delay error of pseudorange and carrier phase observations is as follows:

[0031] ;

[0032] ;

[0033] In the formula, and The ionospheric delay error, representing frequency-dependent pseudorange and carrier phase observations, is expressed in units of . , Indicates the signal propagation path, in meters (m). Represents electron density, with units of electron number. ;

[0034] As shown in the formula for calculating the first-order ionospheric delay error, the ionospheric delay error is related to the signal frequency and the ionospheric electron density. The first-order ionospheric delay errors of the pseudorange and carrier phase are equal in value but opposite in sign, and the frequency... The ionospheric delay of the dual-frequency observations satisfies:

[0035] ;

[0036] For an observation model without an ionospheric array, the pseudorange and carrier phase observation equations can be expressed as:

[0037] Formula B: ;

[0038] In the formula, Let be the combination coefficients, and satisfy . ;

[0039] Combining equations A and B yields the simplified equation for the deionized system:

[0040] ;

[0041] In the formula, This represents the carrier wavelength corresponding to the ionospheric combination observation, in units of... This represents the phase ambiguity of the ionospherically unbounded combined carrier wave, in cycles. and This represents the pseudorange hardware delay at the receiver and satellite ends without an ionospheric combination, in units of... and The hardware delay of the carrier phase at the receiver and satellite ends without ionospheric combination is expressed in weeks. The actual estimated parameters of PPP include: receiver coordinates, receiver clock error, wet component of tropospheric delay in the zenith direction of the station, and ambiguity parameters.

[0042] Preferably, the specific steps of S4 are as follows:

[0043] S4.1: Based on the data from the ephemeris and B2b satellite correction data messages, determine whether there are B2b correction data. If so, further determine whether the correction data has exceeded its validity period. If it is available, proceed to the next step.

[0044] S4.2: Check if the IODC of B2b correction information type 2 and information type 4, i.e., the version numbers of the track correction and clock error correction, match. If they match, proceed to the next step.

[0045] S4.3: Based on the IODN matching number in information type 2, determine whether the satellite has an IODC number and IODN that match in the broadcast IODN. If so, proceed to the next step.

[0046] S4.4: Based on the satellite's ephemeris data, determine the satellite's position and clock bias. Using the clock bias correction algorithm, recover the B2b precise ephemeris and calculate the satellite's position coordinates.

[0047] S4.5: Based on the URA index in information type 2, determine the satellite positioning accuracy of the satellite and assign corresponding weights to the satellite.

[0048] Compared with the prior art, the beneficial effects of the present invention are:

[0049] 1. Compared with the traditional RTK model, this invention does not require the establishment of fixed base stations. Each device in this invention is a master station, and high-precision positioning of each device is achieved through the PPP-B2b model.

[0050] 2. By calculating high-quality positioning accuracy, the timing accuracy of the equipment is improved. Combining positioning technology with timing technology improves the synchronization accuracy between equipment. Attached Figure Description

[0051] Figure 1 This is a schematic diagram of the method flow of the present invention;

[0052] Figure 2 This is a flowchart for removing abnormal clock bias data during anomaly handling;

[0053] Figure 3 This is a flowchart of S4;

[0054] Figure 4 This is a flowchart of S5. Detailed Implementation

[0055] To facilitate understanding of the technical content of this invention by those skilled in the art, the invention will be further described in detail below with reference to the accompanying drawings and specific examples. It should be understood that the specific examples described herein are merely illustrative and not intended to limit the scope of the invention.

[0056] Example 1

[0057] like Figure 1 , Figure 4 As shown, a time synchronization method based on BeiDou satellites includes the following steps:

[0058] S1: Obtain the ephemeris, satellite observation data, and B2b satellite correction data messages of BeiDou-3 satellites through the receiver;

[0059] S2: Due to interference from external conditions such as the environment, the collected satellite observation data may contain pseudorange noise and carrier cycle slips. Using these abnormal observation data for parameter estimation will affect the solution accuracy of PPP. Therefore, it is necessary to perform anomaly processing on the collected satellite observation data. Anomaly processing includes preprocessing algorithms and cycle slip detection algorithms. The preprocessing algorithm is used to remove satellite observation data with observation quality lower than the predetermined value, and the cycle slip detection algorithm is used to remove abnormal clock error data in the satellite observation data.

[0060] S3: Combining ephemeris and anomaly-processed satellite observation data, pseudorange single-point positioning is performed using the PPP floating-point solution function model;

[0061] S4: Based on the ephemeris and B2b satellite correction data, determine the satellite's position and clock bias. Using a clock bias correction algorithm, recover the precise ephemeris of the B2b satellite and calculate its position coordinates.

[0062] S5: Input the pseudorange point positioning result and the satellite's position coordinates into a Kalman filter for filtering. Then, perform PPP precise point positioning on the filtered data to calculate the receiver's position and clock error.

[0063] S6: Input the receiver clock difference data obtained in S5 into the clock control module to periodically adjust the frequency and phase of the clock to complete time synchronization.

[0064] In this embodiment, compared with the traditional RTK model, the present invention does not require the establishment of a fixed base station. Each device in the present invention is a master station. High-precision positioning of each device is achieved through the PPP-B2b model, and the timing accuracy of the devices is improved by calculating high-quality positioning accuracy. The combination of positioning technology and timing technology improves the synchronization accuracy between devices.

[0065] Example 2

[0066] The difference between this embodiment and embodiment 1 is that the preprocessing algorithm in S2 is specifically as follows:

[0067] Suppose the satellite number received by the receiver is , ,…… The corresponding satellite signal-to-noise ratio is , ,…… The corresponding satellite elevation angle is , ,…… The preprocessing method involves eliminating satellites that do not meet the requirements based on the signal-to-noise ratio and elevation angle. The specific formula is as follows:

[0068] ;

[0069] ;

[0070] The sum of the two formulas above gives the calculated satellite mass value, i.e.:

[0071] ;

[0072] If the current satellite number Calculation results If the value is less than 60, the satellite is marked as unhealthy and removed from the list.

[0073] Example 3

[0074] The difference between this embodiment and Embodiment 1 is that, as Figure 2 As shown, the cycle slip detection algorithm in S2 uses a phase carrier wide-lane combination of combined observations, and the observation equation is:

[0075] ;

[0076] ;

[0077] ;

[0078] in , Indicates the current frequency. , This represents the corresponding carrier phase observation value. Indicates the time of data collection. This indicates the cycle slip detection threshold, and E represents the satellite elevation angle. Indicates the wavelength of the wide lane;

[0079] pass If the collected observation value is found to be within the cycle slip detection threshold through iterative calculation, the data is retained for subsequent calculations.

[0080] Example 4

[0081] The difference between this embodiment and Embodiment 1 is that in S3, the PPP floating-point solution function model uses both approximate range and carrier phase observations for parameter estimation, and constructs observation equations to describe the functional relationship between the observations and various estimated parameters and correction errors. Let the frequency be... The pseudorange and carrier phase observation equations are as follows:

[0082] Formula A: ;

[0083] Among them, superscript Indicates satellite, subscript Indicates receiver, subscript Indicates the frequency band of the observed value. Represents pseudorange observations, in units of This represents the carrier phase observation value, in cycles. This represents the geometric distance between the phase centers of the receiver antenna and the satellite antenna, in units of... and Clock bias indicates the time difference between the receiver and satellite clocks relative to a time reference, and is measured in units of 1. The unit representing the speed of light in a vacuum is _____. and The unit representing the frequency-dependent pseudorange hardware delay at the receiver and satellite ends is . This represents the frequency-dependent hardware delay of the receiver and satellite carrier phase, measured in cycles. This represents the frequency-dependent ionospheric delay, measured in meters (m). This represents the tropospheric delayed wet component mapping function related to the satellite elevation angle. The tropospheric delay moisture component in the zenith direction of the observation station is shown in meters. Represents frequency The corresponding carrier wavelength, in meters (m). This represents the carrier phase integer ambiguity, measured in cycles. and The unmodeling error noise represents the pseudorange and carrier phase observations, in units of Furthermore, errors in the observation equations, such as receiver and satellite antenna phase center deviation and variation, relativistic effects, Earth rotation correction, and tropospheric delay dry component, have been corrected in advance using the model.

[0084] According to radio propagation theory, pseudorange propagates in the ionosphere at group velocity, while carrier phase propagates at phase velocity. Therefore, the distance measured by pseudorange observations is longer than the actual distance, while the distance measured by carrier phase observations is shorter than the actual distance. The formulas for calculating the first-order ionospheric delay error of pseudorange and carrier phase observations are as follows:

[0085] ;

[0086] ;

[0087] In the formula, and The ionospheric delay error, representing frequency-dependent pseudorange and carrier phase observations, is expressed in units of . , Indicates the signal propagation path, in meters (m). Represents electron density, with units of electron number. ;

[0088] As shown in the formula for calculating the first-order ionospheric delay error, the ionospheric delay error is related to the signal frequency and the ionospheric electron density. The first-order ionospheric delay errors of the pseudorange and carrier phase are equal in value but opposite in sign, and the frequency... The ionospheric delay of the dual-frequency observations satisfies:

[0089] ;

[0090] Based on the specific numerical relationship between carrier frequency and ionospheric delay, users can eliminate low-order ionospheric delay errors using a linear combination of dual-frequency carrier observations and pseudorange measurements. This combination of observations is also known as the ionosphere-free (IF) combined model. The ionosphere-free combined model is the earliest and most widely used observation model in PPP applications, and its pseudorange and carrier phase observation equations can be expressed as:

[0091] Formula B: ;

[0092] In the formula, Let be the combination coefficients, and satisfy . ;

[0093] Combining equations A and B yields the simplified equation for the deionized system:

[0094] ;

[0095] In the formula, This represents the carrier wavelength corresponding to the ionospheric combination observation, in units of... This represents the phase ambiguity of the ionospherically unbounded combined carrier wave, in cycles. and This represents the pseudorange hardware delay at the receiver and satellite ends without an ionospheric combination, in units of... and The hardware delay of the carrier phase at the receiver and satellite ends without ionospheric combination is expressed in cycles. Since the clock bias parameter and the hardware delay are linearly related, they cannot be directly separated during parameter estimation. The actual estimated clock bias is the sum of the true clock bias and the corresponding hardware delay error. Therefore, the actual estimated parameters of PPP include: receiver coordinates, receiver clock bias, wet component of tropospheric delay in the zenith direction of the station, and ambiguity parameters.

[0096] Example 5

[0097] The difference between this embodiment and Embodiment 1 is that, as Figure 3 As shown, the specific steps of S4 are as follows:

[0098] S4.1: Based on the data from the ephemeris and B2b satellite correction data messages, determine whether there are B2b correction data. If so, further determine whether the correction data has exceeded its validity period. If it is available, proceed to the next step.

[0099] S4.2: Check if the IODC of B2b correction information type 2 and information type 4, i.e., the version numbers of the track correction and clock error correction, match. If they match, proceed to the next step.

[0100] S4.3: Based on the IODN matching number in information type 2, determine whether the satellite has an IODC number and IODN that match in the broadcast IODN. If so, proceed to the next step.

[0101] S4.4: Based on the satellite's ephemeris data, determine the satellite's position and clock bias. Using the clock bias correction algorithm, recover the B2b precise ephemeris and calculate the satellite's position coordinates.

[0102] S4.5: Based on the URA index in information type 2, determine the satellite positioning accuracy of the satellite and assign corresponding weights to the satellite.

[0103] Test results: Assuming the calculated current clock bias is... The time is The calculated clock difference for the next second is The time is The clock difference calculated at the nth second is... The time is The average clock difference is calculated using a 60-second sliding window as Feq( The average clock difference of the next sliding window is Feq( The average clock difference for the Nth iteration is Feq( If the frequency control word is ), then the frequency control word is .

[0104] The synchronization accuracy between devices controlled in this way can reach within 3ns.

Claims

1. A time synchronization method based on BeiDou satellites, characterized in that, Includes the following steps: S1: Obtain the ephemeris, satellite observation data, and B2b satellite correction data messages of BeiDou-3 satellites through the receiver; S2: Perform anomaly handling, which includes a preprocessing algorithm and a cycle slip detection algorithm. The preprocessing algorithm is used to remove satellite observation data that is below a predetermined value, and the cycle slip detection algorithm is used to remove abnormal clock error data in the satellite observation data. The preprocessing algorithm is as follows: Suppose the satellite number received by the receiver is , ,…… The corresponding satellite signal-to-noise ratio is , ,…… The corresponding satellite elevation angle is , ,…… The preprocessing method involves eliminating satellites that do not meet the requirements based on the signal-to-noise ratio and elevation angle. The specific formula is as follows: ； ； The sum of the two formulas above gives the calculated satellite mass value, i.e.: ； If the current satellite number Calculation results If the value is less than 60, the satellite is marked as unhealthy and removed from the list. The cycle slip detection algorithm uses a phase carrier wide-lane combination of combined observations, and the observation equation is: ； ； ； in , Indicates the current frequency. , This represents the corresponding carrier phase observation value. Indicates the time of data collection. This indicates the cycle slip detection threshold, and E represents the satellite elevation angle. Indicates the wavelength of the wide lane; pass If the collected observation value is found to be within the cycle slip detection threshold through iterative calculation, then the data is retained for subsequent calculations. S3: Combining ephemeris and anomaly-processed satellite observation data, pseudorange single-point positioning is performed using the PPP floating-point solution function model; S4: Based on the ephemeris and B2b satellite correction data, determine the position and clock bias of the B2b satellite, recover the precise ephemeris of the B2b satellite through the clock bias correction algorithm, and calculate the position coordinates of the satellite. S5: Input the pseudorange point positioning result and the satellite's position coordinates into a Kalman filter for filtering. Then, perform PPP precise point positioning on the filtered data to calculate the receiver's position and clock error. S6: Input the receiver clock difference data obtained in S5 into the clock control module to periodically adjust the frequency and phase of the clock to complete time synchronization.

2. The time synchronization method based on BeiDou satellite according to claim 1, characterized in that, In S3, the PPP floating-point solution function model uses both analog range and carrier phase observations for parameter estimation, and constructs observation equations to describe the functional relationship between the observations and various estimated parameters and correction errors. Let the frequency be... The pseudorange and carrier phase observation equations are as follows: Formula A: ; Among them, superscript Indicates satellite, subscript Indicates receiver, subscript Indicates the frequency band of the observed value. Represents pseudorange observations, in units of This represents the carrier phase observation value, in cycles. This represents the geometric distance between the phase centers of the receiver antenna and the satellite antenna, in units of... and Clock bias indicates the time difference between the receiver and satellite clocks relative to a time reference, and is measured in units of 1. The unit representing the speed of light in a vacuum is _____. , and The unit representing the frequency-dependent pseudorange hardware delay at the receiver and satellite ends is . , This represents the frequency-dependent hardware delay of the receiver and satellite carrier phase, measured in cycles. This represents the frequency-dependent ionospheric delay, measured in meters (m). This represents the tropospheric delayed wet component mapping function related to the satellite elevation angle. The tropospheric delay moisture component in the zenith direction of the observation station is shown in meters. Represents frequency The corresponding carrier wavelength, in meters (m). This represents the carrier phase integer ambiguity, measured in cycles. and The unmodeling error noise represents the pseudorange and carrier phase observations, in units of ; The formula for calculating the first-order ionospheric delay error of pseudorange and carrier phase observations is as follows: ； ； In the formula, and The ionospheric delay error, representing frequency-dependent pseudorange and carrier phase observations, is expressed in units of . , Indicates the signal propagation path, in meters (m). Represents electron density, with units of electron number. ; As shown in the formula for calculating the first-order ionospheric delay error, the ionospheric delay error is related to the signal frequency and the ionospheric electron density. The first-order ionospheric delay errors of the pseudorange and carrier phase are equal in value but opposite in sign, and the frequency... The ionospheric delay of the dual-frequency observations satisfies: ； For an observation model without an ionospheric array, the pseudorange and carrier phase observation equations can be expressed as: Formula B: ; In the formula, Let be the combination coefficients, and satisfy . ; Combining equations A and B yields the simplified equation for the deionized system: ； In the formula, This represents the carrier wavelength corresponding to the ionospheric combination observation, in units of... This represents the phase ambiguity of the ionospherically unbounded combined carrier wave, in cycles. and This represents the pseudorange hardware delay at the receiver and satellite ends without an ionospheric combination, in units of... and The hardware delay of the carrier phase at the receiver and satellite ends without ionospheric combination is expressed in weeks. The actual estimated parameters of PPP include: receiver coordinates, receiver clock error, wet component of tropospheric delay in the zenith direction of the station, and ambiguity parameters.

3. The time synchronization method based on BeiDou satellite according to claim 1, characterized in that, The specific steps of S4 are as follows: S4.1: Based on the data from the ephemeris and B2b satellite correction data messages, determine whether there are B2b correction data. If so, further determine whether the correction data has exceeded its validity period. If it is available, proceed to the next step. S4.2: Check if the IODC of B2b correction information type 2 and information type 4, i.e., the version numbers of the track correction and clock error correction, match. If they match, proceed to the next step. S4.3: Based on the IODN matching number in information type 2, determine whether the satellite has an IODC number and IODN that match in the broadcast IODN. If so, proceed to the next step. S4.4: Based on the satellite's ephemeris data, determine the satellite's position and clock bias. Using the clock bias correction algorithm, recover the B2b precise ephemeris and calculate the satellite's position coordinates. S4.5: Based on the URA index in information type 2, determine the satellite positioning accuracy of the satellite and assign corresponding weights to the satellite.

Images