Method for autonomous integrity monitoring of navigation satellite messages

Abstract

The application provides a navigation satellite message autonomous integrity monitoring method, which comprises the following steps: calculating first state information at a predetermined time through extrapolation based on historical navigation message parameters of an on-board message, and calculating second state information according to navigation message parameters of a new on-board message on the predetermined time, calculating a first difference between the two; obtaining regional system signals and global system signals on the same time, and calculating corresponding state information and a second difference between the state information; judging whether the first difference and the second difference exceed corresponding alarm thresholds, and if both exceed, generating an alarm information of the message integrity. Thus, the application can realize on-board message autonomous integrity monitoring, can monitor and alarm abnormalities of the message caused by uplink injection, malicious tampering and the like, and prevent error positioning results caused by the error message broadcast to the ground users.

Description

Technical Field

[0001] This invention relates to the field of navigation satellite technology, and in particular to a method for autonomous monitoring of the integrity of navigation satellite messages. Background Technology

[0002] Integrity refers to the system's ability to promptly alert the user and terminate the signal when the navigation system's error exceeds the permissible limit and it fails to perform its intended navigation function. Integrity is a crucial indicator of a navigation system's performance, and its performance is vital for users seeking high security. Organizations such as the International Civil Aviation Organization (ICAO) have established strict regulations regarding the integrity indicators of satellite navigation systems.

[0003] Generally, the ways to achieve operational integrity include: ground control system monitoring, satellite autonomous integrity monitoring, and receiver autonomous integrity monitoring. Due to the limitation of global ground station deployment, the BeiDou Navigation Satellite System has limited capabilities for monitoring the integrity of overseas satellites. Therefore, to achieve the system's integrity indicators, it cannot rely solely on ground system monitoring; it requires a combination of satellite autonomous integrity monitoring and ground monitoring to ensure the BeiDou system's basic global operational integrity service indicators.

[0004] Currently, satellite autonomous integrity monitoring does not include monitoring the correctness of messages. That is, if there are no errors in the message but errors occur in its own physical quantities, the onboard integrity monitoring unit cannot detect or alarm. This will create an integrity risk. When the ground user terminal receives the erroneous message, it will generate incorrect positioning results, affecting navigation services.

[0005] Therefore, there is an urgent need for a method that can effectively address the aforementioned problems in the existing technology. Summary of the Invention

[0006] The purpose of this invention is to provide a method for autonomous integrity monitoring of navigation satellite messages, which can monitor and alert on anomalies caused by uplink injection, malicious tampering, etc., and prevent incorrect positioning results from being caused by erroneous messages being broadcast to ground users.

[0007] To achieve the above objectives, the present invention provides a method for autonomous integrity monitoring of navigation satellite messages, comprising the following steps:

[0008] Obtain the first navigation message parameters from the first message, and extrapolate the first state information at a predetermined time based on the first navigation message parameters. Calculate the first difference between the first state information and the second state information; wherein the second state information is calculated based on the second navigation message parameters in the second message annotated at the predetermined time; the first message is a historical annotated message relative to the second message, and the first state information and the second state information are satellite positions or clock biases;

[0009] The system acquires regional system signals and global system signals transmitted at the same time, and calculates corresponding third and fourth state information based on the message parameters in the regional system signals and global system signals, respectively, and calculates a second difference between the third and fourth state information; wherein the third and fourth state information are satellite positions or clock biases.

[0010] Determine whether the first difference exceeds the first alarm threshold and whether the second difference exceeds the second alarm threshold;

[0011] If the first difference exceeds the first alarm threshold and the second difference exceeds the second alarm threshold, an incompleteness alarm message is generated.

[0012] Optionally, the first status information and the second status information are satellite positions;

[0013] The step of obtaining the first navigation message parameters in the first message, extrapolating and calculating the first state information at a predetermined time based on the first navigation message parameters, and calculating the first difference between the first state information and the second state information includes:

[0014] Obtain the first ephemeris parameter from the first message;

[0015] Calculate the extrapolated duration between the ephemeris reference time and the predetermined time in the first ephemeris parameter;

[0016] The orbital elements of the satellite are extrapolated based on the extrapolation duration.

[0017] Calculate the three-dimensional coordinates of the first satellite at the predetermined time based on the extrapolated orbital elements;

[0018] The three-dimensional coordinates of the second satellite are calculated based on the second ephemeris parameters in the second message sent at the predetermined time.

[0019] Calculate the position difference between the three-dimensional coordinates of the first satellite and the three-dimensional coordinates of the second satellite.

[0020] Optionally, the extrapolation duration is calculated based on the following formula:

[0021] t k =tt oe Where t is the predetermined time, t oe The ephemeris reference time is mentioned above;

[0022] The extrapolation calculation of the orbital elements of the satellite based on the extrapolation duration includes:

[0023] According to the formula A0 = A ref+ΔA, calculate the semi-major axis; where ΔA is the deviation of the semi-major axis from the reference value, A red At a depth of 27,906,100 meters in Earth's orbit, A ref The inclination geosynchronous orbit or geostationary orbit has a depth of 42,162,200m.

[0024] According to the formula Calculate the satellite's average angular velocity; where μ is the Earth's gravitational constant;

[0025] According to the formula Calculate the deviation of the satellite's average angular velocity; where Δn0 is the difference between the satellite's average velocity and the calculated value. The rate of change of the difference between the satellite's average velocity and the calculated value;

[0026] According to formula n A =n0+Δn A Calculate the satellite's average angular velocity after deviation correction;

[0027] According to formula M k =M0+n A t k Calculate the mean anomaly angle; where M0 is the mean anomaly angle at the ephemeris reference time;

[0028] According to formula M k =E k -esinE k The eccentricity angle is calculated iteratively; where e is the eccentricity.

[0029] According to the formula Calculate the true anterior angle;

[0030] According to the formula φ k =v k +ω calculates the latitude argument; ω is the perigee deflection.

[0031] According to the formula Calculate the latitude argument correction term, radial distance correction term, and orbital inclination correction term separately; where δu k For the dimension argument correction term, δr k For radial distance correction, δi k For the track inclination correction term, C us C is the amplitude of the sinusoidal harmonic correction term for the latitude argument. uc C is the amplitude of the cosine harmonic correction term for the latitude argument. rs C is the amplitude of the sinusoidal harmonic correction term for the orbital radius. rc C is the amplitude of the cosine harmonic correction term for the orbital radius. is C is the amplitude of the sinusoidal harmonic correction term for the orbital inclination angle. icThe amplitude of the cosine harmonic correction term for the orbital inclination angle;

[0032] According to formula r k =A k (1-ecosE k )+δr k Calculate the corrected radial distance; where A k For the major half-axis;

[0033] According to the formula Calculate the satellite's coordinates in its orbital plane;

[0034] According to the formula Calculate the corrected longitude of the ascending node, where Ω0 represents the Earth's rotation rate, and Ω0 represents the longitude of the reference ascending node. The rate of change of longitude at the ascending node;

[0035] According to the formula Calculate the orbital inclination at the ephemeris reference time; where i0 is the orbital inclination at the reference time. This represents the rate of change of the orbital inclination angle.

[0036] Optionally, calculating the three-dimensional coordinates of the first satellite at the predetermined time based on the extrapolated orbital elements includes:

[0037] According to the formula Calculate the three-dimensional coordinates of the first satellite at the predetermined time;

[0038] The calculation of the three-dimensional coordinates of the second satellite based on the second ephemeris parameters in the second message sent at the predetermined time includes:

[0039] According to the second ephemeris parameters and formula in the second message sent at the predetermined time. And let t k =0 to calculate the three-dimensional coordinates of the second satellite.

[0040] Optionally, the first state information and the second state information are clock differences;

[0041] The step of obtaining the first navigation message parameters in the first message, extrapolating and calculating the first state information at a predetermined time based on the first navigation message parameters, and calculating the first difference between the first state information and the second state information includes:

[0042] Obtain the first clock bias parameter from the first message; wherein the first clock bias parameter includes the satellite clock deviation coefficient, drift coefficient, drift rate coefficient, and clock bias parameter reference time;

[0043] Based on the first clock difference parameter, extrapolate and calculate the first clock difference at the predetermined time;

[0044] The second clock error is calculated based on the second clock error parameters in the second message recorded at the predetermined time; wherein the second clock error parameters include the satellite clock bias coefficient, drift coefficient, drift rate coefficient, and clock error parameter reference time;

[0045] Calculate the time difference between the first clock difference and the second clock difference.

[0046] Optionally, the first clock bias is calculated based on the following formula:

[0047] Δt 外推 =a0+a1(tt) oc )+a2(tt oc ) 2 +Δt r ;

[0048]

[0049] Where, Δt r This is a relativistic correction term, where e is the satellite orbital eccentricity. E is the square root of the semi-major axis of the satellite's orbit. k Let a0 be the satellite orbital perihelion angle, a1 be the satellite clock deviation coefficient in the first message, a2 be the drift coefficient in the first message, and t be the drift rate coefficient in the first message. oc The reference time for the clock difference parameter in the first message; Where μ is the gravitational constant and C is the speed of light.

[0050] Optionally, the second clock bias is calculated based on the following formula:

[0051] Δt=a0′+Δt r ;

[0052] Where a0′ is the satellite clock offset coefficient in the second message.

[0053] Optionally, before determining whether the first difference exceeds the first alarm threshold and whether the second difference exceeds the second alarm threshold, the method further includes:

[0054] Multiple sets of the first difference and the second difference are obtained as sample values;

[0055] Statistical analysis was performed on multiple sets of the first difference and the second difference to obtain the first standard deviation and the mean;

[0056] Subtract the mean from the sample value and divide by the first standard deviation to obtain the normalized parameter value;

[0057] Calculate the probability distribution and second standard deviation of the normalized parameter;

[0058] Multiply the second standard deviation by the standard deviation inflation factor to obtain the standard deviation;

[0059] Based on the standard deviation and integrity risk requirements, the first alarm threshold and the second alarm threshold are calculated.

[0060] Optionally, if the first difference exceeds a first alarm threshold and the second difference exceeds a second alarm threshold, then an incompleteness alarm is generated.

[0061] If the first difference exceeds the first alarm threshold and the second difference exceeds the second alarm threshold, then the corresponding frequency ephemeris or clock error incompleteness flag is set.

[0062] The navigation satellite message autonomous integrity monitoring method of this invention calculates first state information at a predetermined time based on extrapolation of navigation message parameters from historically injected messages, and calculates second state information based on navigation message parameters of newly injected messages at the predetermined time, and calculates a first difference between the two; it also acquires regional system signals and global system signals injected at the same time, and calculates their corresponding state information and a second difference between them respectively; it determines whether both the first and second differences exceed the corresponding alarm thresholds. If both exceed the thresholds, an integrity alarm is generated. Thus, this invention fills the gap in on-board autonomous integrity monitoring, realizing on-board message autonomous integrity monitoring, and can monitor and alarm for anomalies in messages caused by uplink injection, malicious tampering, etc., preventing erroneous messages from being broadcast to ground users and leading to incorrect positioning results. Attached Figure Description

[0063] Figure 1 This diagram illustrates a flowchart of the autonomous integrity monitoring method for navigation satellite messages provided in an embodiment of the present invention.

[0064] Figure 2 The flowchart of the autonomous integrity monitoring method for navigation satellite messages described in this invention is shown. Detailed Implementation

[0065] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the invention.

[0066] It should be noted that references to "an embodiment," "embodiment," "example embodiment," etc., in this specification refer to the described embodiment including specific features, structures, or characteristics, but not every embodiment must include these specific features, structures, or characteristics. Furthermore, such expressions do not refer to the same embodiment. Moreover, when describing specific features, structures, or characteristics in conjunction with embodiments, whether or not explicitly described, it is indicated that incorporating such features, structures, or characteristics into other embodiments is within the knowledge of those skilled in the art.

[0067] Furthermore, certain terms are used in the specification and subsequent claims to refer to specific components or parts. Those skilled in the art will understand that manufacturers may use different names or terms to refer to the same component or part. This specification and subsequent claims do not distinguish components or parts by differences in name, but rather by differences in function. The terms "comprising" and "including" used throughout the specification and subsequent claims are open-ended and should be interpreted as "including but not limited to." Additionally, the term "connection" here includes any direct and indirect electrical connection means. Indirect electrical connection means include connections made through other means.

[0068] To address the technical challenge that current satellite autonomous integrity monitoring does not include monitoring the correctness of messages (i.e., when there are no errors in the message but errors occur in its own physical quantities), and the onboard integrity monitoring unit cannot monitor and alarm, this invention provides a navigation satellite message autonomous integrity monitoring method to solve the problem of the inability to perform onboard autonomous monitoring of messages.

[0069] Figure 1 This invention illustrates a method for monitoring the autonomous integrity of navigation satellite messages according to an embodiment of the present invention. This method is mainly applied to monitoring the autonomous integrity of onboard navigation satellite messages and includes the following steps:

[0070] S101: Obtain the first navigation message parameters in the first message, and extrapolate the first state information at a predetermined time based on the first navigation message parameters, and calculate the first difference between the first state information and the second state information; wherein, the second state information is calculated based on the second navigation message parameters in the second message at the predetermined time; the first message is a historical message relative to the second message, and the first state information and the second state information are satellite positions or clock biases; that is, in this embodiment, in step S101, the first state information at the predetermined time is extrapolated based on the first navigation message parameters carried in the historical message. Extrapolation calculation refers to extrapolating the first state information at the predetermined time based on a known data sequence or functional relationship, beyond the known data range. Extrapolation is the process of predicting values ​​within a given interval; that is, using the first navigation message parameters of historically recorded messages to predict a future moment to obtain the predicted state of the satellite. Specifically, in aerospace engineering, extrapolation is the process of extrapolating the orbital parameters (such as position, velocity, clock bias, etc.) of a satellite or spacecraft within a known time period to predict its future orbital position or clock bias and other satellite state information. After obtaining the first state information at a predetermined moment through extrapolation, a new recorded message is obtained at that predetermined moment, and the second state information is calculated based on the second navigation message parameters in the new recorded message. Then, the first state information is compared with the second state information to obtain the first difference between the two.

[0071] In this embodiment, step S101 specifically involves extrapolating historical messages at the same frequency point to a predetermined time to obtain the predicted first state information. Then, the second state information at the predetermined time is calculated using the newly added message. The two states calculated at the same frequency point are compared, i.e., the difference between the predicted result and the actual measurement result is compared. In this embodiment, the predetermined time is specifically the current time, the first message refers to messages added before the current time, and the second message is the new message added at the current time.

[0072] S102: Acquire the regional system signal and global system signal transmitted at the same time, and calculate the corresponding third and fourth state information based on the message parameters in the regional system signal and global system signal, respectively, and calculate the second difference between the third and fourth state information; wherein, the third and fourth state information are satellite positions or clock biases; for satellites of the BeiDou-3 system that simultaneously broadcast regional system signals (smooth transition signals B1I, B3I) and global system signals (global signals B1C, B2a), the regional system messages (B1I, B3I) and global messages (B1C, B2a) are transmitted separately through different information types. After receiving the data, the satellites broadcast it in the corresponding signal components; the regional signal messages and global signal messages are calculated and generated separately using monitoring quantities of different signals, and are transmitted through different information types during the transmission process, and can be considered to be independent of each other; that is, in this embodiment, by acquiring the message parameters received simultaneously by signal branches at different frequency points, calculating their respective corresponding state information and comparing them, the difference between the two is obtained.

[0073] S103: Determine whether the first difference exceeds the first alarm threshold and whether the second difference exceeds the second alarm threshold.

[0074] S104: If the first difference exceeds the first alarm threshold and the second difference exceeds the second alarm threshold, an incompleteness alarm message is generated. That is, if the first difference calculated in step S101 and the second difference calculated in step S102 both exceed their respective alarm thresholds, an alarm is issued.

[0075] See Figure 2 This embodiment combines the comparison of historical and newly injected messages within the same frequency point with the comparison of different signal types between different frequency points to determine whether an alarm is needed based on the comparison results.

[0076] In this embodiment, all state information can be either satellite position or clock bias; the following will explain these two cases respectively:

[0077] In an optional implementation, where all the state information to be calculated is satellite position, step S101 specifically includes:

[0078] 1. Obtain the first ephemeris parameter from the first message; in practice, the satellite's ephemeris is uploaded via information type 1, with an uploading strategy of once per hour, allowing only parameters for the current day and the next two days to be uploaded. Normal startup mode is generally used, meaning the uploading is completed one minute before startup at the exact hour of the next hour after receiving this type of information. The satellite can obtain the ephemeris by receiving information type 1.

[0079] The satellite ephemeris in information type 1 referred to in this embodiment consists of 18 orbital parameters, plus 8 bits of ephemeris data age and 2 bits of satellite orbit type (SatType), for a total of 45 bits. The definition and characteristics of the ephemeris parameters are shown in Table 1 below.

[0080] Table 1:

[0081]

[0082]

[0083]

[0084] 2. Calculate the extrapolated duration between the ephemeris reference time and the predetermined time in the first ephemeris parameter; that is, using the ephemeris reference time t... oe As a benchmark, the difference t between the extrapolated time (i.e., the predetermined time) and the reference time is first calculated. k For normal startup mode, the extrapolated duration is generally the hour of the next hour; for real-time startup mode, it is the actual activation time. Specifically, the extrapolated duration is calculated based on the following formula: t k =tt oe Where t is the predetermined time, t oe For ephemeris reference time.

[0085] 3. Extrapolate and calculate the orbital elements of the satellite based on the extrapolated duration.

[0086] In practice, the extrapolation duration t can be calculated according to the following steps. k orbital elements:

[0087] According to the formula A0 = A ref +ΔA, calculate the semi-major axis; where ΔA is the deviation of the semi-major axis from the reference value, A ref At a depth of 27,906,100 meters in Earth's orbit, A ref The inclination is 42162200m in inclined geosynchronous orbit or geostationary orbit.

[0088] According to the formula Calculate the satellite's average angular velocity; where μ is the Earth's gravitational constant, μ = 3.986004418 × 10⁻⁶. 14 m 3 / s 2 .

[0089] According to the formula Calculate the deviation of the satellite's average angular velocity; where Δn0 is the difference between the satellite's average velocity and the calculated value. The rate of change of the difference between the satellite's average velocity and the calculated value.

[0090] According to formula n A =n0+Δn A Calculate the satellite's average angular velocity after deviation correction.

[0091] According to formula M k =M0+n A t k Calculate the mean aperimeter angle; where M0 is the mean aperimeter angle at the ephemeris reference time.

[0092] According to formula M k =E k -esinE k The eccentricity angle is calculated iteratively; where e is the eccentricity.

[0093] According to the formula Calculate the true anterior angle.

[0094] According to the formula φ k =v k +ω calculates the latitude argument; ω is the perigee deflection.

[0095] According to the formula Calculate the latitude argument correction term, radial distance correction term, and orbital inclination correction term separately; where δu k For the dimension argument correction term, δr k For radial distance correction, δi k For the track inclination correction term, C us C is the amplitude of the sinusoidal harmonic correction term for the latitude argument. uc C is the amplitude of the cosine harmonic correction term for the latitude argument. rs C is the amplitude of the sinusoidal harmonic correction term for the orbital radius. rc C is the amplitude of the cosine harmonic correction term for the orbital radius. is C is the amplitude of the sinusoidal harmonic correction term for the orbital inclination angle. ic The amplitude of the cosine harmonic correction term for the orbital inclination angle.

[0096] According to formula r k =A k (1-ecosE k )+δr k Calculate the corrected radial distance; where A k For the major half-axis;

[0097] According to the formula Calculate the satellite's coordinates in its orbital plane;

[0098] According to the formula Calculate the corrected longitude of the ascending node, where The Earth's rotation rate, Ω0 is the reference longitude of the ascending node. The rate of change of longitude at the ascending node;

[0099] According to the formula Calculate the orbital inclination at the ephemeris reference time; where i0 is the orbital inclination at the reference time. This represents the rate of change of the orbital inclination angle.

[0100] 4. Calculate the three-dimensional coordinates of the first satellite at the predetermined time based on the extrapolated orbital elements. Specifically, calculate the three-dimensional coordinates of the first satellite at the predetermined time based on the calculated orbital elements, i.e.:

[0101] According to the formula Calculate the three-dimensional coordinates (X, Y, Z) of the first satellite at the predetermined time. k Y k Z k );

[0102] 5. Calculate the three-dimensional coordinates of the second satellite based on the second ephemeris parameters in the second message sent at the predetermined time;

[0103] Specifically, the above-mentioned methods for calculating the number of orbital elements can be used, and let t k =0 to calculate the three-dimensional coordinates of the second satellite, that is, according to the second ephemeris parameters and formula in the second message sent at the predetermined time. And let t k =0 to calculate the three-dimensional coordinates (X0, Y0, Z0) of the second satellite.

[0104] 6. Calculate the position difference between the three-dimensional coordinates of the first satellite and the three-dimensional coordinates of the second satellite. Combining the above calculation results, the position difference (ΔX, ΔY, ΔZ) calculated using the old ephemeris extrapolation and the new ephemeris can be calculated as follows: (X...) k Y k Z k )-(X0, Y0, Z0).

[0105] Furthermore, in step S102, the regional signal message and the global signal message are calculated and generated separately using different signal monitoring quantities. They are uploaded using different information types during the uploading process and can be considered to be independent of each other. The corresponding ephemeris parameters in the two messages are compared.

[0106] The ephemeris format in another information type 62 is shown in Table 2 below. Table 2 is specifically used to show the basic navigation information of the regional system.

[0107] Table 2:

[0108]

[0109]

[0110]

[0111]

[0112] By comparing the message format with that in Table 1, it can be seen that due to the smooth transition, the differences between the satellite's semi-major axis, the rate of change of right ascension of the ascending node, and the satellite's average motion speed and the calculated values ​​in the global signal are slightly different. Based on the parameters in Table 2, the third and fourth state information corresponding to the regional system signal and the global system signal uploaded at the same time can be calculated according to the following steps:

[0113] 1. Calculate the major semi-axis

[0114] 2. Calculate the satellite's average angular velocity. Among them, the gravitational constant μ = 3.986004418 × 10 14 m 3 / s 2 .

[0115] 3. Calculate the time difference t between the observed epoch and the reference epoch. k =tt oa ; where t oa For reference time;

[0116] 4. Calculate the angle M at the point of approach. k =M0+n0t k Where n0 is the average angular velocity of the satellite;

[0117] 5. Iteratively calculate the near-point angle E k M k =E k -esinE k .

[0118] 6. Calculate the true anterior angle v. k ,

[0119] 7. Calculate the argument φ. k =v k +ω, where ω is the perigee deflection angle;

[0120] 8. According to the formula Calculate the latitude argument correction term, radial distance correction term, and orbital inclination correction term separately; where δu k For the dimension argument correction term, δr k For radial distance correction, δi k For the track inclination correction term, C usC is the amplitude of the sinusoidal harmonic correction term for the latitude argument. uc C is the amplitude of the cosine harmonic correction term for the latitude argument. rs C is the amplitude of the sinusoidal harmonic correction term for the orbital radius. rc C is the amplitude of the cosine harmonic correction term for the orbital radius. is C is the amplitude of the sinusoidal harmonic correction term for the orbital inclination angle. ic The amplitude of the cosine harmonic correction term for the orbital inclination angle.

[0121] 9. According to the formula r k =A k (1-ecosE k )+δr k Calculate the corrected radial distance; where A k For the major half-axis;

[0122] 10. Calculate the satellite's coordinates in the orbital plane.

[0123] 11. Calculate the corrected longitude of the ascending node. in The Earth's rotation rate, Ω0 is the reference longitude of the ascending node. The rate of change of longitude at the ascending node;

[0124] 11. The orbital inclination at the reference time is i = i0 + IDOT·t k +δi k ;

[0125] 12. Calculate the satellite's coordinates (X, Y, F, Z) in the BDCS (BeiDou) coordinate system using the following formula. k ′, Y k ′, Z k ′):

[0126]

[0127] In this embodiment, "the same time" in step S102 specifically refers to the predetermined time. Therefore, the satellite position difference (ΔX′, ΔY′, ΔZ′) obtained using the regional system signal and the global system signal is calculated as follows: k ′, Y k ′, Z k ′)-(X0, Y0, Z0).

[0128] The first difference (ΔX, ΔY, ΔZ) and the second difference (ΔX′, ΔY′, ΔZ′) calculated above are compared with the corresponding first alarm threshold and second alarm threshold respectively. It is determined whether the first difference (ΔX, ΔY, ΔZ) exceeds the preset first alarm threshold and whether the second difference (ΔX′, ΔY′, ΔZ′) exceeds the preset second alarm threshold. If both exceed the preset threshold, an ephemeris parameter alarm is issued.

[0129] In another optional implementation, where the state information to be calculated is a clock difference, step S101 specifically includes:

[0130] Obtain the first clock bias parameter from the first message; wherein the first clock bias parameter includes the satellite clock bias coefficient, drift coefficient, drift rate coefficient, and clock bias parameter reference time.

[0131] Based on the first clock difference parameter, the first clock difference at the predetermined time is extrapolated and calculated; specifically, the first clock difference is calculated based on the following formula:

[0132] Δt 外推 =a0+a1(tt) oc )+a2(tt oc ) 2 +Δt r ;

[0133]

[0134] Where, Δt r This is a relativistic correction term, where e is the satellite orbital eccentricity. E is the square root of the semi-major axis of the satellite's orbit. k Let a0 be the satellite orbital perihelion angle, a1 be the satellite clock deviation coefficient in the first message, a2 be the drift coefficient in the first message, and t be the drift rate coefficient in the first message. oc The reference time for the clock difference parameter in the first message; Where μ is the gravitational constant and C is the speed of light.

[0135] The second clock error is calculated based on the second clock error parameters in the second message injected at the predetermined time. These parameters include the satellite clock bias coefficient, drift coefficient, drift rate coefficient, and clock error parameter reference time. Specifically, the real-time clock error Δt is calculated using the corresponding satellite clock bias coefficient a0′, drift coefficient a1′, and drift rate coefficient a2′ in the newly injected message. Since t = T, the actual clock error Δt is calculated at this point. oc Therefore, substituting into the above formula, we get ΔT=a0′+ΔT r That is, the second clock bias is calculated based on the following formula:

[0136] ΔT=a0′+Δt r ;

[0137] Where a0′ is the satellite clock offset coefficient in the second message.

[0138] Calculate the time difference between the first clock difference and the second clock difference. Then, calculate Δt. 外推 The first difference can be obtained by subtracting Δt from the first difference.

[0139] Furthermore, for satellites of the BeiDou-3 system that simultaneously broadcast regional system signals (smooth transition signals B1I and B3I) and global system signals (global signals B1C and B2a), the regional system messages (B1I and B3I) and global messages (B1C and B2a) are respectively uploaded using separate information types. After receiving the signals, the satellites broadcast them in the corresponding signal components. The regional signal messages and global signal messages are calculated and generated separately using different signal monitoring quantities, and are uploaded using different information types during the uploading process; they can be considered independent of each other. Step S102 may specifically include:

[0140] The real clock difference Δt′ is calculated based on the satellite clock deviation coefficient, drift coefficient, and drift rate coefficient in the smooth transition signal.

[0141] The real clock difference Δt is calculated based on the satellite clock bias coefficient, drift coefficient, and drift rate coefficient in the global system signal.

[0142] The second difference can be obtained by subtracting Δt′ from Δt.

[0143] The first and second differences calculated above are compared with the corresponding first and second alarm thresholds to determine whether the first difference exceeds the preset first alarm threshold and whether the second difference exceeds the preset second alarm threshold. If both exceed the preset first alarm threshold, a clock difference alarm is issued.

[0144] Optionally, before step S103, the following steps are also included:

[0145] Multiple sets of the first difference and the second difference are obtained as sample values; specifically, the aforementioned method can be used to calculate the first difference between the extrapolated ephemeris and clock difference of multiple sets of the same frequency point and the newly injected ephemeris and clock difference, as well as the second difference between the real-time ephemeris and clock difference of different frequency points, thereby obtaining multiple sets of sample data.

[0146] Statistical analysis was performed on multiple groups of the first and second differences to obtain the first standard deviation and mean.

[0147] The sample value is subtracted from the mean and divided by the first standard deviation to obtain the normalized parameter value.

[0148] Calculate the probability distribution and second standard deviation of the normalized parameter.

[0149] The second standard deviation is multiplied by the standard deviation inflation factor (f) to obtain the standard deviation. That is, in this embodiment, the standard deviation inflation factor (f) can be used to determine the final variance. The original standard deviation σ is then multiplied by f, and the final standard deviation is fσ.

[0150] Based on the aforementioned standard deviation and integrity risk requirements, the first alarm threshold and the second alarm threshold are calculated. Since the parameter values ​​follow the formula X ~ N(μ, f 2 σ 2 If the data follows a Gaussian distribution, the probability of exceeding the threshold is Pr(|x-μ|>fσ). For multi-source data, since the data from each source conforms to a Gaussian distribution, the probability distribution of the data from multiple sources can be calculated using a multivariate Gaussian distribution:

[0151]

[0152] For both same-frequency and different-frequency data sources, the probability of simultaneously exceeding the threshold is Pr(x1-μ1|>f1σ1&x2-μ2|>f2σ2), which is the integrity risk corresponding to the threshold Th. Based on the integrity risk requirement, the corresponding thresholds T1=μ1+σ1 and T2=μ2+σ2 are obtained. For the ephemeris and clock bias parameter distributions, two thresholds are obtained, namely… and

[0153] In this embodiment, after the satellite is in orbit, the difference between the newly received ephemeris and the historical ephemeris at different frequencies is used. When both differences simultaneously exceed a threshold, the system will determine the appropriate threshold. If an alarm is triggered, a flag indicating incomplete ephemeris information for the corresponding frequency point can be set.

[0154] It is also possible to use the calculated new receive clock bias and historical clock bias, the difference between the new receive clock bias at different frequencies, and the fact that both differences simultaneously exceed a threshold. If an alarm is triggered, a clock error flag can be set for the corresponding frequency point.

[0155] Of course, in other embodiments, the present invention can also be applied to the integrity monitoring of other status information of satellite navigation messages; that is, the present invention is scalable and can be applied to the monitoring of other data in the message, without being limited by the type of data.

[0156] In summary, the autonomous integrity monitoring method for navigation satellite messages described in this invention can simultaneously monitor the orbit and clock bias within the message, solving the problem of the inability to perform autonomous onboard monitoring of messages. This invention combines data from both co-frequency and inter-frequency sources, performing joint comparison and judgment through multi-source data, thus improving the accuracy of anomaly detection; it avoids situations where the ephemeris itself is normal, but anomalies occur during data transmission, leading to erroneous alarms (false alarms). Furthermore, it fully utilizes ground-based ephemeris data used during satellite operation, without adding new data, and features low implementation cost and on-orbit upgradeability.

[0157] Of course, the present invention may have other various embodiments. Without departing from the spirit and essence of the present invention, those skilled in the art can make various corresponding changes and modifications according to the present invention, but these corresponding changes and modifications should all fall within the protection scope of the appended claims.

Claims

1. A method for autonomous integrity monitoring of navigation satellite messages, characterized in that, Including the following steps: Obtain the first navigation message parameters from the first message, and extrapolate the first state information at a predetermined time based on the first navigation message parameters. Calculate the first difference between the first state information and the second state information; wherein the second state information is calculated based on the second navigation message parameters in the second message annotated at the predetermined time; the first message is a historical annotated message relative to the second message, and the first state information and the second state information are satellite positions or clock biases; The system acquires regional system signals and global system signals transmitted at the same time, and calculates corresponding third and fourth state information based on the message parameters in the regional system signals and global system signals, respectively, and calculates a second difference between the third and fourth state information; wherein the third and fourth state information are satellite positions or clock biases. Determine whether the first difference exceeds the first alarm threshold and whether the second difference exceeds the second alarm threshold; If the first difference exceeds the first alarm threshold and the second difference exceeds the second alarm threshold, an incompleteness alarm message is generated.

2. The autonomous integrity monitoring method for navigation satellite messages according to claim 1, characterized in that, The first status information and the second status information are satellite positions; The step of obtaining the first navigation message parameters in the first message, extrapolating and calculating the first state information at a predetermined time based on the first navigation message parameters, and calculating the first difference between the first state information and the second state information includes: Obtain the first ephemeris parameter from the first message; Calculate the extrapolated duration between the ephemeris reference time and the predetermined time in the first ephemeris parameter; The orbital elements of the satellite are extrapolated based on the extrapolation duration. Calculate the three-dimensional coordinates of the first satellite at the predetermined time based on the extrapolated orbital elements; The three-dimensional coordinates of the second satellite are calculated based on the second ephemeris parameters in the second message sent at the predetermined time. Calculate the position difference between the three-dimensional coordinates of the first satellite and the three-dimensional coordinates of the second satellite.

3. The autonomous integrity monitoring method for navigation satellite messages according to claim 2, characterized in that, The extrapolation duration is calculated based on the following formula: t k =tt oe Where t is the predetermined time, t oe The ephemeris reference time is as stated; The extrapolation calculation of the orbital elements of the satellite based on the extrapolation duration includes: According to the formula A0 = A ref +ΔA, calculate the semi-major axis; where ΔA is the deviation of the semi-major axis from the reference value, A ref At a depth of 27,906,100 meters in Earth's orbit, A ref The inclination geosynchronous orbit or geostationary orbit has a depth of 42,162,200m. According to the formula Calculate the satellite's average angular velocity; where μ is the Earth's gravitational constant; According to the formula Calculate the deviation of the satellite's average angular velocity; where Δn0 is the difference between the satellite's average velocity and the calculated value. The rate of change of the difference between the satellite's average velocity and the calculated value; According to formula n A =n0+Δn A Calculate the satellite's average angular velocity after deviation correction; According to formula M k =M0+n A t k Calculate the mean anomaly angle; where M0 is the mean anomaly angle at the ephemeris reference time; According to formula M k =E k -esinE k The eccentricity angle is calculated iteratively; where e is the eccentricity. According to the formula Calculate the true anterior angle; According to the formula φ k =v k +ω calculates the latitude argument; ω is the perigee deflection. According to the formula Calculate the latitude argument correction term, radial distance correction term, and orbital inclination correction term separately; where δu k For the latitude argument correction term, δr k For radial distance correction, δi k For the track inclination correction term, C us C is the amplitude of the sinusoidal harmonic correction term for the latitude argument. uc C is the amplitude of the cosine harmonic correction term for the latitude argument. rs C is the amplitude of the sinusoidal harmonic correction term for the orbital radius. rc C is the amplitude of the cosine harmonic correction term for the orbital radius. is C is the amplitude of the sinusoidal harmonic correction term for the orbital inclination angle. ic The amplitude of the cosine harmonic correction term for the orbital inclination angle; According to formula r k =A k (1-ecosE k )+δr k Calculate the corrected radial distance; where A k For the major half-axis; According to the formula Calculate the satellite's coordinates in its orbital plane; According to the formula Calculate the corrected longitude of the ascending node, where Ω0 represents the Earth's rotation rate, and Ω0 represents the longitude of the reference ascending node. The rate of change of longitude at the ascending node; According to the formula Calculate the orbital inclination at the ephemeris reference time; where i0 is the orbital inclination at the reference time. This represents the rate of change of the orbital inclination angle.

4. The autonomous integrity monitoring method for navigation satellite messages according to claim 3, characterized in that, The three-dimensional coordinates of the first satellite at the predetermined time, calculated based on the extrapolated orbital elements, include: According to the formula Calculate the three-dimensional coordinates of the first satellite at the predetermined time; The calculation of the three-dimensional coordinates of the second satellite based on the second ephemeris parameters in the second message sent at the predetermined time includes: According to the second ephemeris parameters and formula in the second message sent at the predetermined time. And let t k =0 to calculate the three-dimensional coordinates of the second satellite.

5. The autonomous integrity monitoring method for navigation satellite messages according to claim 1, characterized in that, The first state information and the second state information are clock differences; The step of obtaining the first navigation message parameters in the first message, extrapolating and calculating the first state information at a predetermined time based on the first navigation message parameters, and calculating the first difference between the first state information and the second state information includes: Obtain the first clock bias parameter from the first message; wherein the first clock bias parameter includes the satellite clock deviation coefficient, drift coefficient, drift rate coefficient, and clock bias parameter reference time; Based on the first clock difference parameter, extrapolate and calculate the first clock difference at the predetermined time; The second clock error is calculated based on the second clock error parameters in the second message recorded at the predetermined time; wherein the second clock error parameters include the satellite clock bias coefficient, drift coefficient, drift rate coefficient, and clock error parameter reference time; Calculate the time difference between the first clock difference and the second clock difference.

6. The autonomous integrity monitoring method for navigation satellite messages according to claim 5, characterized in that, The first clock bias is calculated based on the following formula: Δt 外推 =a0+a1(t-t oc )+a2(t-t oc ) 2 +Δt r ; Where, Δt r This is a relativistic correction term, where e is the satellite orbital eccentricity. E is the square root of the semi-major axis of the satellite's orbit. k Let a0 be the satellite orbital perihelion angle, a1 be the satellite clock deviation coefficient in the first message, a2 be the drift coefficient in the first message, and t be the drift rate coefficient in the first message. oc The reference time for the clock difference parameter in the first message; Where μ is the gravitational constant and C is the speed of light.

7. The autonomous integrity monitoring method for navigation satellite messages according to claim 6, characterized in that, The second clock bias is calculated based on the following formula: Δt=a0′+Δt r ; Where a0′ is the satellite clock offset coefficient in the second message.

8. The method for autonomous integrity monitoring of navigation satellite messages according to claim 1, characterized in that, Before determining whether the first difference exceeds the first alarm threshold and whether the second difference exceeds the second alarm threshold, the method further includes: Multiple sets of the first difference and the second difference are obtained as sample values; Statistical analysis was performed on multiple sets of the first difference and the second difference to obtain the first standard deviation and the mean; Subtract the mean from the sample value and divide by the first standard deviation to obtain the normalized parameter value; Calculate the probability distribution and second standard deviation of the normalized parameter; Multiply the second standard deviation by the standard deviation inflation factor to obtain the standard deviation; Based on the standard deviation and integrity risk requirements, the first alarm threshold and the second alarm threshold are calculated.

9. The method for autonomous integrity monitoring of navigation satellite messages according to claim 1, characterized in that, If the first difference exceeds the first alarm threshold and the second difference exceeds the second alarm threshold, then an incompleteness alarm is generated. If the first difference exceeds the first alarm threshold and the second difference exceeds the second alarm threshold, then the corresponding frequency ephemeris or clock error incompleteness flag is set.

Images