Autonomous intact monitoring method for navigation satellite messages

Abstract

The invention provides an autonomous intact monitoring method for a navigation satellite message, which comprises the following steps: calculating first state information at a predetermined moment by extrapolating navigation message parameters of a historical top-injection message, calculating second state information according to navigation message parameters of a new top-injection message top-injection at the predetermined moment, and calculating a first difference value between the first state information and the second state information; acquiring a regional system signal and a global system signal which are uploaded at the same moment, and respectively calculating state information corresponding to the regional system signal and the global system signal and a second difference value between the state information; and judging whether the first difference value and the second difference value both exceed corresponding alarm thresholds or not, and if yes, generating incomplete alarm information. Thus, autonomous integrity monitoring of the on-satellite telegraph text can be achieved, monitoring and warning can be conducted on abnormity of the telegraph text caused by uplink injection, malicious tampering and the like, and the situation that a wrong positioning result is caused after a wrong telegraph text is broadcasted to a ground user is prevented.

Description

Technical Field

[0001] The present invention relates to the field of navigation satellite technology, and in particular to a navigation satellite message autonomous integrity monitoring method. Background Art

[0002] Integrity refers to the ability of the system to promptly alarm and notify the user to terminate the signal when the error of the navigation system exceeds the allowed limit and cannot perform the specified navigation work. Integrity is one of the important indicators of the navigation system. The integrity performance index is crucial to high-security users. Organizations such as ICAO have strictly regulated the integrity indicators of satellite navigation systems.

[0003] Generally speaking, the ways to achieve integrity functions include: ground operation and control system monitoring, satellite autonomous integrity monitoring, and receiver autonomous integrity monitoring. The BeiDou satellite navigation system is affected by the fact that the ground system cannot be deployed globally, and its integrity monitoring capabilities for overseas satellites are limited. Therefore, in order to achieve the integrity indicators of the system, it cannot rely solely on ground system monitoring. It is necessary to use satellite autonomous integrity monitoring and ground monitoring to jointly ensure the global basic integrity service indicators of the BeiDou system.

[0004] At present, the autonomous integrity of satellites does not include the monitoring of the correctness of telegrams. That is, if there is no bit error in the telegram, but there is an error in the physical quantity itself, the onboard integrity monitoring unit cannot monitor and alarm. This will create an integrity risk. After the ground user terminal receives the erroneous telegram, it will generate an erroneous positioning result, affecting the navigation service.

[0005] Therefore, there is an urgent need for a method that can effectively solve the above-mentioned prior art problems. Summary of the invention

[0006] The purpose of the present invention is to provide a navigation satellite message autonomous integrity monitoring method for realizing on-board message autonomous integrity monitoring, which can monitor and alarm the message anomalies caused by uplink injection, malicious tampering, etc., to prevent erroneous message broadcasting to ground users and cause erroneous positioning results.

[0007] In order to achieve the above object, the present invention provides a method for autonomous integrity monitoring of navigation satellite messages, comprising the steps of:

[0008] Acquire a first navigation message parameter in a first message, and extrapolate and calculate first state information at a predetermined time based on the first navigation message parameter, and calculate a first difference between the first state information and second state information; wherein the second state information is calculated based on the second navigation message parameter in the second message annotated at the predetermined time; the first message is a historically annotated message relative to the second message, and the first state information and the second state information are satellite positions or clock errors;

[0009] Acquire a regional system signal and a global system signal registered at the same time, and calculate corresponding third state information and fourth state information based on message parameters in the regional system signal and the global system signal, and calculate a second difference between the third state information and the fourth state information; wherein the third state information and the fourth state information are satellite positions or clock errors;

[0010] Determining whether the first difference exceeds a first alarm threshold, and whether the second difference exceeds a second alarm threshold;

[0011] If the first difference exceeds a first alarm threshold and the second difference exceeds a 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 acquiring of the first navigation message parameter in the first message, and extrapolating and calculating the first state information at a predetermined time based on the first navigation message parameter, and calculating the first difference between the first state information and the second state information comprises:

[0014] Obtaining first ephemeris parameters in a first message;

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

[0016] Performing extrapolation calculation on each orbital element of the satellite based on the extrapolation duration;

[0017] Calculate the three-dimensional coordinates of the first satellite at the predetermined time according to the orbital elements obtained by extrapolation;

[0018] Calculating the three-dimensional coordinates of the second satellite according to the second ephemeris parameters in the second message annotated at the predetermined time;

[0019] The position difference between the first satellite three-dimensional coordinates and the second satellite three-dimensional coordinates is calculated.

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

[0021] t k =tt oe ; Wherein, t is the predetermined time, t oe is the ephemeris reference time;

[0022] The extrapolating calculation of each orbital element of the satellite based on the extrapolation duration includes:

[0023] According to the formula A0=A ref+ΔA, calculate the major semi-axis; where ΔA is the deviation of the major semi-axis from the reference value, A red In medium earth orbit, it is equal to 27906100m, A ref In inclined geosynchronous orbit or geostationary orbit equal to 42162200m;

[0024] According to the formula Calculate the average angular velocity of the satellite; where μ is the 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 motion rate and the calculated value, The rate of change of the difference between the satellite's average motion velocity and the calculated value;

[0026] According to the formula n A =n0+Δn A , calculate the average angular velocity of the satellite after bias correction;

[0027] According to the 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 the formula M k =E k -esinE k , iteratively calculate the eccentric anomaly angle; where e is the eccentricity;

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

[0030] According to the formula φ k =v k +ω, calculate the latitude argument; ω is the perigee declination;

[0031] According to the formula Calculate the latitude argument correction, radial distance correction and orbital inclination correction respectively; where δu k is the dimension argument correction term, δr k is the radial distance correction term, δi k is the orbit inclination correction term, C us is the amplitude of the sine harmonic correction term of the latitude argument, C uc is the amplitude of the cosine harmonic correction term of the latitude argument, C rs is the amplitude of the sinusoidal harmonic correction term of the orbital radius, C rc is the amplitude of the cosine harmonic correction term of the orbital radius, C is is the amplitude of the sinusoidal harmonic correction term of the orbital inclination, C icis the amplitude of the cosine harmonic correction term of the orbital inclination;

[0032] According to the formula r k =A k (1-ecosE k )+δr k , calculate the corrected radial distance; where A k is the major semi-axis;

[0033] According to the formula Calculate the coordinates of the satellite in the orbital plane;

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

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

[0036] Optionally, calculating the three-dimensional coordinates of the first satellite at the predetermined time according to the orbital elements obtained by extrapolation includes:

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

[0038] The step of calculating the second satellite three-dimensional coordinates according to the second ephemeris parameter in the second message annotated at the predetermined time includes:

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

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

[0041] The acquiring of the first navigation message parameter in the first message, and extrapolating and calculating the first state information at a predetermined time based on the first navigation message parameter, and calculating the first difference between the first state information and the second state information comprises:

[0042] Acquire a first clock error parameter in the first message; wherein the first clock error parameter includes a satellite clock deviation coefficient, a drift coefficient, a drift rate coefficient, and a clock error parameter reference time;

[0043] extrapolating and calculating a first clock difference at a predetermined time according to the first clock difference parameter;

[0044] Calculate the second clock error according to the second clock error parameter in the second message annotated at the predetermined time; wherein the second clock error parameter includes a satellite clock deviation coefficient, a drift coefficient, a drift rate coefficient and a clock error parameter reference time;

[0045] The time difference between the first clock difference and the second clock difference is calculated.

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

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

[0048]

[0049] Among them, Δt r is the relativistic correction term, e is the satellite orbit eccentricity, is the square root of the semi-major axis of the satellite orbit, E k is the satellite orbit anomaly angle, a0 is the satellite clock bias coefficient in the first message, a1 is the drift coefficient in the first message, a2 is the drift rate coefficient in the first message, t oc The reference time of the clock error parameter in the first message; Where μ is the gravitational constant and C is the speed of light.

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

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

[0052] Wherein, a0′ is the satellite clock deviation coefficient in the second message.

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

[0054] Acquire multiple groups of the first difference values ​​and the second difference values ​​as sample values;

[0055] Performing statistics on multiple groups of the first differences and the second differences to obtain a first standard deviation and a mean;

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

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

[0058] Multiplying the second standard deviation by the standard deviation expansion factor to obtain a standard deviation;

[0059] The first alarm threshold and the second alarm threshold are calculated based on the standard deviation and the integrity risk requirement.

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

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

[0062] The method for autonomous integrity monitoring of navigation satellite messages described in the present invention calculates the first state information at a predetermined time by extrapolating the navigation message parameters of the historically injected messages, and calculates the second state information according to the navigation message parameters of the newly injected messages injected at the predetermined time, and calculates the first difference between the two; and obtains the regional system signal and the global system signal injected at the same time, and respectively calculates the state information corresponding to the two and the second difference between the state information; determines whether the first difference and the second difference both exceed the corresponding alarm threshold, and if both exceed, generates an alarm message of incompleteness. In this way, the present invention fills the gap in autonomous integrity monitoring on satellites, realizes autonomous integrity monitoring of satellite messages, and can monitor and alarm abnormalities of messages caused by uplink injection, malicious tampering, etc., to prevent erroneous messages from being broadcast to ground users and causing erroneous positioning results. BRIEF DESCRIPTION OF THE DRAWINGS

[0063] Figure 1 A flowchart showing the steps of the method for autonomous integrity monitoring of navigation satellite messages provided by an embodiment of the present invention is shown;

[0064] Figure 2 A flowchart of the method for autonomous integrity monitoring of navigation satellite messages according to the present invention is shown. DETAILED DESCRIPTION

[0065] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention.

[0066] It should be noted that references to "one embodiment", "embodiment", "example embodiment", etc. in this specification refer to the embodiment described, which may include specific features, structures or characteristics, but not every embodiment must include these specific features, structures or characteristics. In addition, such expressions do not refer to the same embodiment. Furthermore, when describing specific features, structures or characteristics in conjunction with an embodiment, whether or not there is an explicit description, it has been shown that it is within the knowledge of those skilled in the art to combine such features, structures or characteristics into other embodiments.

[0067] In addition, certain words are used in the specification and subsequent claims to refer to specific components or parts. Those with ordinary knowledge in the relevant field should understand that manufacturers can use different nouns or terms to refer to the same component or part. This specification and subsequent claims do not use differences in names as a way to distinguish components or parts, but use differences in the functions of components or parts as the criteria for distinction. "Including" and "including" mentioned throughout the specification and subsequent claims are open-ended terms and should be interpreted as "including but not limited to". In addition, the word "connected" here includes any direct and indirect electrical connection means. Indirect electrical connection means include connection through other devices.

[0068] Aiming at the technical problem that the current satellite autonomous integrity does not include the monitoring of the correctness of the telegram, that is, the telegram has no bit errors, but the physical quantity itself has errors, and the on-board integrity monitoring unit cannot monitor and alarm, the present invention provides a navigation satellite telegram autonomous integrity monitoring method to solve the problem that the telegram cannot be autonomously monitored on the satellite.

[0069] Figure 1 A method for autonomous integrity monitoring of navigation satellite messages provided by an embodiment of the present invention is shown. The method is mainly applied to autonomous integrity monitoring of messages on navigation satellites, and includes the following steps:

[0070] S101: Acquire first navigation message parameters in the first message, and extrapolate and calculate first state information at a predetermined time based on the first navigation message parameters, and calculate a 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 historically annotated message relative to the second message, and the first state information and the second state information are satellite positions or clock differences; that is, in this embodiment, in step S101, extrapolate and calculate the first state information at a predetermined time based on the first navigation message parameters carried in the historically annotated message, and the extrapolation calculation refers to calculating the first state information at a predetermined time based on a known data sequence or a functional relationship, when the first state information exceeds the known data range. The process of predicting a numerical value in a range of; that is, predicting a future moment based on the first navigation message parameters of the historical injection message, and obtaining the predicted state of the satellite; specifically, in aerospace engineering, extrapolation calculation is to perform extrapolation calculation based on the orbital parameters (such as position, speed, clock error, etc.) of the satellite or spacecraft in a known time period to predict its future orbital position or clock error and other satellite state information; after the first state information at a predetermined moment is obtained by extrapolation calculation, a new injection message at the predetermined moment is further obtained, and the second state information is calculated based on the second navigation message parameters in the new injection message, and then the first state information is compared with the second state information to obtain a first difference between the two.

[0071] That is, step S101 of this embodiment specifically obtains the predicted first state information by extrapolating the historically injected telegrams at the same frequency point to the predetermined time, and then uses the newly injected telegrams to calculate the second state information at the predetermined time, and compares the two state information calculated at the same frequency point, that is, compares the difference between the predicted result and the actual measurement result. The predetermined time of this embodiment is specifically the current time, the first telegram refers to the telegram injected before the current time, and the second telegram is the new telegram injected at the current time.

[0072] S102: Acquire the regional system signal and the global system signal that are uploaded at the same time, and calculate the corresponding third state information and fourth state information based on the message parameters in the regional system signal and the global system signal, and calculate the second difference between the third state information and the fourth state information; wherein the third state information and the fourth state information are satellite positions or clock differences; for the satellites of the BeiDou-3 system that simultaneously broadcast the regional system signal (smooth transition signal B1I, B3I) and the global system signal (global signal B1C, B2a), the regional system message (B1I, B3I) and the global message (B1C, B2a) are uploaded respectively through different information types. After receiving, the satellite broadcasts the corresponding signal component; the regional signal message and the global signal message are calculated and generated separately using different signal monitoring quantities, and are uploaded through different information types during the upload process, and can be considered to be independent of each other; that is, this embodiment obtains the message parameters received by the signal branches of different frequency points at the same time, calculates the corresponding state information and compares them, so as to obtain the corresponding difference between the two.

[0073] S103: Determine whether the first difference exceeds a first alarm threshold, and whether the second difference exceeds a 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 corresponding alarm thresholds, an alarm prompt is issued.

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

[0076] The status information of this embodiment may be all satellite positions or all clock errors; the following will describe these two situations respectively:

[0077] In an optional implementation manner, each state information to be calculated is a satellite position, and step S101 specifically includes:

[0078] 1. Obtain the first ephemeris parameter in the first message; in the specific implementation process, the satellite's ephemeris is injected through the injection information type 1, and its injection strategy is to inject once every hour, and only the parameters of the current day and the next two days can be injected. Generally, the normal startup mode is adopted, that is, the injection is completed at the next full hour after receiving this type of information, and 1 minute before startup. The satellite can obtain the ephemeris by receiving the information type 1.

[0079] The satellite ephemeris in the information type 1 referred to in this embodiment is composed of 18 orbit parameters, plus 8 bits of ephemeris data age, and 2 bits of satellite orbit type (SatType), 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 time from the ephemeris reference time in the first ephemeris parameter to the predetermined time; that is, take the ephemeris reference time t oe As a benchmark, first calculate the difference t between the extrapolated time (i.e., the predetermined time) and the reference time k For normal startup mode, it is generally the next hour, and for real-time startup mode, it is the actual activation time. Specifically, the extrapolated duration is calculated based on the following formula: k =tt oe ; where t is the scheduled time, t oe is the ephemeris reference time.

[0085] 3. Perform extrapolation calculation on each orbital element of the satellite based on the extrapolation duration.

[0086] In specific implementation, the extrapolation time t can be calculated according to the following steps: k The number of orbital elements:

[0087] According to the formula A0=A ref +ΔA, calculate the major semi-axis; where ΔA is the deviation of the major semi-axis from the reference value, A ref In medium earth orbit, it is equal to 27906100m, A ref In inclined geosynchronous orbit or geostationary orbit is equal to 42162200m.

[0088] According to the formula Calculate the average angular velocity of the satellite; where μ is the 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 motion rate and the calculated value, The rate of change of the difference between the satellite's average motion rate and the calculated value.

[0090] According to the formula n A =n0+Δn A , calculate the bias-corrected average angular velocity of the satellite.

[0091] According to the 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.

[0092] According to the formula M k =E k -esinE k , iteratively calculate the eccentric anomaly angle; where e is the eccentricity.

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

[0094] According to the formula φ k =v k +ω, calculate the latitude argument; ω is the perigee declination;

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

[0096] According to the formula r k =A k (1-ecosE k )+δr k , calculate the corrected radial distance; where A k is the major semi-axis;

[0097] According to the formula Calculate the coordinates of the satellite in the orbital plane;

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

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

[0100] 4. Calculate the three-dimensional coordinates of the first satellite at the predetermined time according to the orbital elements obtained by extrapolation. Specifically, the three-dimensional coordinates of the first satellite at the predetermined time are calculated according to the orbital elements calculated above, that is:

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

[0102] 5. Calculate the three-dimensional coordinates of the second satellite according to the second ephemeris parameters in the second message annotated at the predetermined time;

[0103] Specifically, the calculation method of the above orbital elements can be used, and t k = 0 to calculate the second satellite three-dimensional coordinates, that is, according to the second ephemeris parameters and formula in the second message annotated 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 first satellite 3D coordinates and the second satellite 3D coordinates. Combining the above calculation results, the position difference calculated using the old ephemeris extrapolation and the new ephemeris can be calculated as (ΔX, ΔY, ΔZ) = (X k , Y k , Z k )-(X0, Y0, Z0).

[0105] Further, in step S102, the regional signal message and the global signal message are calculated and generated separately by monitoring different signal monitoring quantities. Different information types are added during the annotation process and can be considered 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 the following Table 2, which is specifically used to show basic navigation information of a regional system.

[0107] Table 2:

[0108]

[0109]

[0110]

[0111]

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

[0113] 1. Calculate the major semiaxis

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

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

[0116] 4. Calculate the mean anomaly angle M k =M0+n0t k , where n0 is the average angular velocity of the satellite;

[0117] 5. Iterative calculation of the eccentric angle E k , M k =E k -esinE k .

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

[0119] 7. Calculate the dimension angle φ k =v k +ω, ω is the perigee declination angle;

[0120] 8. According to the formula Calculate the latitude argument correction, radial distance correction and orbital inclination correction respectively; where δu k is the dimension argument correction term, δr k is the radial distance correction term, δi k is the orbit inclination correction term, C usis the amplitude of the sine harmonic correction term of the latitude argument, C uc is the amplitude of the cosine harmonic correction term of the latitude argument, C rs is the amplitude of the sinusoidal harmonic correction term of the orbital radius, C rc is the amplitude of the cosine harmonic correction term of the orbital radius, C is is the amplitude of the sinusoidal harmonic correction term of the orbital inclination, C ic is the amplitude of the cosine harmonic correction term of the orbital inclination.

[0121] 9. According to formula r k =A k (1-ecosE k )+δr k , calculate the corrected radial distance; where A k is the major semi-axis;

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

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

[0124] 11. Orbital inclination at reference time i = i0 + IDOT·t k +δi k ;

[0125] 12. Calculate the satellite coordinates (X k ′,Y k ′,Z k ′):

[0126]

[0127] In this embodiment, the “same time” in step S102 specifically refers to the predetermined time, so the satellite position difference (ΔX′, ΔY′, ΔZ′) obtained by obtaining the regional system signal and the global system signal is calculated = (X 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 to determine 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, an ephemeris parameter alarm prompt is issued.

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

[0130] Obtain a first clock error parameter in a first message; wherein the first clock error parameter includes a satellite clock deviation coefficient, a drift coefficient, a drift rate coefficient and a clock error parameter reference time.

[0131] According to the first clock difference parameter, a first clock difference at a 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] Among them, Δt r is the relativistic correction term, e is the satellite orbit eccentricity, is the square root of the semi-major axis of the satellite orbit, E k is the satellite orbit anomaly angle, a0 is the satellite clock bias coefficient in the first message, a1 is the drift coefficient in the first message, a2 is the drift rate coefficient in the first message, t oc The reference time of the clock error parameter in the first message; Where μ is the gravitational constant and C is the speed of light.

[0135] The second clock error is calculated according to the second clock error parameter in the second message injected at the predetermined time; wherein the second clock error parameter includes the satellite clock deviation coefficient, the drift coefficient, the drift rate coefficient and the clock error parameter reference time; that is, the real-time clock error Δt is calculated using the corresponding satellite clock deviation coefficient a0′, the drift coefficient a1′ and the drift rate coefficient a2′ in the newly injected message. At this time, since t=T oc , so substituting into the above formula we get ΔT=a0′+ΔT r ; That is, the second clock difference is calculated based on the following formula:

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

[0137] Wherein, a0′ is the satellite clock deviation coefficient in the second message.

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

[0139] Furthermore, for the 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 the global messages (B1C and B2a) are respectively recorded using separate information types. After receiving the signals, the satellites broadcast the corresponding signal components; the regional signal messages and the global signal messages are calculated and generated separately using different signal monitoring quantities, and are recorded using different information types during the recording process, and can be considered independent of each other. Step S102 may specifically include:

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

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

[0142] Then, the second difference value can be obtained by subtracting Δt′ from Δt.

[0143] The first difference and the second difference calculated as above are compared with the corresponding first alarm threshold and second alarm threshold respectively 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, a clock difference alarm prompt is issued.

[0144] Optionally, before step S103, the method further includes:

[0145] Acquire multiple groups of the first difference and the second difference as sample values; specifically, the aforementioned method can be used to calculate multiple groups of first differences between extrapolated ephemeris and clock errors at the same frequency point and newly injected ephemeris and clock errors, as well as second differences between real-time ephemeris and clock errors at different frequencies, thereby obtaining multiple groups of sample data.

[0146] Statistics are performed on multiple groups of the first difference values ​​and the second difference values ​​to obtain a first standard deviation and a mean value.

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

[0148] The probability distribution and the second standard deviation of the normalized parameter are calculated.

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

[0150] Based on the standard deviation and integrity risk requirement, the first alarm threshold and the second alarm threshold are calculated. Since the parameter value obeys X~N(μ, f 2 σ 2 ) distribution, the probability of exceeding the threshold value is Pr(|x-μ|>fσ); for multi-source data, since the data from each source conforms to the Gaussian distribution. For multiple variables that conform to the Gaussian distribution, the probability distribution of data from multiple sources is calculated using the multivariate Gaussian distribution:

[0151]

[0152] For the two data sources of the same-frequency comparison and the different-frequency comparison, the probability of exceeding the threshold at the same time is Pr(x1-μ1|>f1σ1&x2-μ2|>f2σ2), which is the integrity risk corresponding to the threshold Th. According to the integrity risk requirements, the corresponding thresholds T1=μ1+σ1 and T2=μ2+σ2 are obtained, and two thresholds are obtained for the distribution of ephemeris and clock error parameters, namely and

[0153] In this embodiment, after the satellite is in orbit, the calculated new received ephemeris and historical ephemeris are used to calculate the difference between the new received ephemeris at different frequencies. When these two differences exceed the threshold at the same time, Then an alarm is issued, and the corresponding frequency ephemeris incompleteness flag can be set.

[0154] The calculated new receiving clock error and historical clock error can also be used to calculate the difference between the new receiving clock errors at different frequencies. The two differences exceed the threshold at the same time. An alarm is triggered, and a corresponding frequency clock error incompleteness flag can be set.

[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 monitoring other data in the message without being limited by the type of data.

[0156] In summary, the autonomous integrity monitoring method for navigation satellite telegrams described in the present invention can monitor the orbit and clock error in the telegrams at the same time, solving the problem that the telegrams cannot be autonomously monitored on the satellite. The present invention combines the same-frequency point and different-frequency point data, and performs joint comparison and judgment through multi-source data, thereby improving the accuracy of abnormal judgment; avoiding the situation where there is no abnormality in the ephemeris itself, but an abnormality occurs during the data transmission process, and an alarm (false alarm) is mistakenly issued. It also makes full use of the ground-based ephemeris data used in the satellite's on-orbit operation process, does not add new data, and has the characteristics of low implementation cost and can be upgraded on-orbit.

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

Claims

1. A method for autonomous integrity monitoring of navigation satellite messages, characterized in that: Includes steps: Acquire a first navigation message parameter in a first message, and extrapolate and calculate first state information at a predetermined time based on the first navigation message parameter, and calculate a first difference between the first state information and second state information; wherein the second state information is calculated based on the second navigation message parameter in the second message annotated at the predetermined time; the first message is a historically annotated message relative to the second message, and the first state information and the second state information are satellite positions or clock errors; Acquire a regional system signal and a global system signal registered at the same time, and calculate corresponding third state information and fourth state information based on message parameters in the regional system signal and the global system signal, and calculate a second difference between the third state information and the fourth state information; wherein the third state information and the fourth state information are satellite positions or clock errors; Determining whether the first difference exceeds a first alarm threshold, and whether the second difference exceeds a second alarm threshold; If the first difference exceeds a first alarm threshold and the second difference exceeds a second alarm threshold, an incompleteness alarm message is generated.

2. The method for autonomous integrity monitoring of navigation satellite messages according to claim 1, characterized in that: The first status information and the second status information are satellite positions; The acquiring of the first navigation message parameter in the first message, and extrapolating and calculating the first state information at a predetermined time based on the first navigation message parameter, and calculating the first difference between the first state information and the second state information comprises: Obtaining first ephemeris parameters in a first message; Calculate the extrapolated duration between the ephemeris reference time and the predetermined time in the first ephemeris parameter; Performing extrapolation calculation on each orbital element of the satellite based on the extrapolation duration; Calculate the three-dimensional coordinates of the first satellite at the predetermined time according to the orbital elements obtained by extrapolation; Calculating the three-dimensional coordinates of the second satellite according to the second ephemeris parameters in the second message annotated at the predetermined time; The position difference between the first satellite three-dimensional coordinates and the second satellite three-dimensional coordinates is calculated.

3. The method for autonomous integrity monitoring of navigation satellite messages according to claim 1, characterized in that: The extrapolated duration is calculated based on the following formula: t k =tt oe ; Wherein, t is the predetermined time, t oe is the ephemeris reference time; The extrapolating calculation of each orbital element of the satellite based on the extrapolation duration includes: According to the formula A0=A ref +ΔA, calculate the major semi-axis; where ΔA is the deviation of the major semi-axis from the reference value, A red In medium earth orbit, it is equal to 27906100m, A ref In inclined geosynchronous orbit or geostationary orbit equal to 42162200m; According to the formula Calculate the average angular velocity of the satellite; where μ is the 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 motion rate and the calculated value, The rate of change of the difference between the satellite's average motion velocity and the calculated value; According to the formula n A =n0+Δn A , calculate the average angular velocity of the satellite after bias correction; According to the 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 the formula M k =E k -esinE k , iteratively calculate the eccentric anomaly angle; where e is the eccentricity; According to the formula Calculate the true anomaly angle; According to the formula φ k =v k +ω, calculate the latitude argument; ω is the perigee declination; According to the formula Calculate the latitude argument correction, radial distance correction and orbital inclination correction respectively; where δu k is the dimension argument correction term, δr k is the radial distance correction term, δi k is the orbit inclination correction term, C us is the amplitude of the sine harmonic correction term of the latitude argument, C uc is the amplitude of the cosine harmonic correction term of the latitude argument, C rs is the amplitude of the sinusoidal harmonic correction term of the orbital radius, C rc is the amplitude of the cosine harmonic correction term of the orbital radius, C is is the amplitude of the sinusoidal harmonic correction term of the orbital inclination, C ic is the amplitude of the cosine harmonic correction term of the orbital inclination; According to the formula r k =A k (1-ecosE k )+δr k , calculate the corrected radial distance; where A k is the major semi-axis; According to the formula Calculate the coordinates of the satellite in the orbital plane; According to the formula Calculate the corrected longitude of the ascending node, where is the Earth's rotation rate, Ω0 is the reference ascending node longitude, is the rate of change of longitude of ascending node; According to the formula Calculate the orbital inclination of the ephemeris reference time; wherein i0 is the orbital inclination of the reference time, is the orbital inclination change rate.

4. The method for autonomous integrity monitoring of navigation satellite messages according to claim 3, characterized in that: Calculating the three-dimensional coordinates of the first satellite at the predetermined time according to the orbital elements obtained by extrapolation includes: According to the formula Calculate the three-dimensional coordinates of the first satellite at the predetermined time; The step of calculating the second satellite three-dimensional coordinates according to the second ephemeris parameter in the second message annotated at the predetermined time includes: According to the second ephemeris parameter and formula in the second message annotated at the predetermined time And let t k =0 to calculate the three-dimensional coordinates of the second satellite.

5. The method for autonomous integrity monitoring of navigation satellite messages according to claim 1, characterized in that: The first status information and the second status information are clock differences; The acquiring of the first navigation message parameter in the first message, and extrapolating and calculating the first state information at a predetermined time based on the first navigation message parameter, and calculating the first difference between the first state information and the second state information comprises: Acquire a first clock error parameter in the first message; wherein the first clock error parameter includes a satellite clock deviation coefficient, a drift coefficient, a drift rate coefficient, and a clock error parameter reference time; extrapolating and calculating a first clock difference at a predetermined time according to the first clock difference parameter; Calculate the second clock error according to the second clock error parameter in the second message annotated at the predetermined time; wherein the second clock error parameter includes a satellite clock deviation coefficient, a drift coefficient, a drift rate coefficient and a clock error parameter reference time; The time difference between the first clock difference and the second clock difference is calculated.

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

7. The method for autonomous integrity monitoring of navigation satellite messages according to claim 6, characterized in that: The second clock difference is calculated based on the following formula: Δt=a0 ′ +Δt r ; Among them, a0 ′ is the satellite clock deviation 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 a first alarm threshold and whether the second difference exceeds a second alarm threshold, the method further includes: Acquire multiple groups of the first difference values ​​and the second difference values ​​as sample values; Performing statistics on multiple groups of the first differences and the second differences to obtain a first standard deviation and a mean; Subtract the mean from the sample value and divide the result by the first standard deviation to obtain a normalized parameter value; Calculate the probability distribution and the second standard deviation of the normalized parameter; Multiplying the second standard deviation by the standard deviation expansion factor to obtain a standard deviation; The first alarm threshold and the second alarm threshold are calculated based on the standard deviation and the integrity risk requirement.

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

Images