A GNSS outlier detection method and system based on residual poly-chi-square distribution used on satellite and storage medium

Abstract

The application provides a GNSS outlier detection method, system and storage medium based on residual multi-order chi-square distribution used on a satellite, comprising the following steps: step one, acquiring the position, velocity and satellite time in the earth-fixed system output by the GNSS; step two, converting the output data of the GNSS into orbital 14 plane roots; step three, calculating the residual and the residual rate of change to make a first heavy outlier judgment, then calculating one-dimensional chi-square test outliers, then calculating two-dimensional chi-square test outliers after the one-dimensional chi-square test is passed, calculating five-dimensional chi-square test outliers after the two-dimensional chi-square test is passed, if all the tests are passed, updating the parameters; if one test is not passed, determining that it is an outlier, and simultaneously counting the continuous invalid times, when the continuous invalid times are greater than a threshold value, the threshold value is adjusted along with the increase of the continuous invalid times, so that the system divergence is avoided. The application has the beneficial effects that the application can obviously improve the rejection rate of outliers, and provides correct data for the use on the satellite.

Description

Technical Field

[0001] This invention relates to the field of spacecraft applications, and in particular to a GNSS outlier detection method, system, and storage medium based on residual multi-order chi-square distribution for use on satellite. Background Technology

[0002] Satellites rely on onboard GNSS (Global Navigation Satellite System) units to provide real-time orbit and attitude reference information, and the accuracy of the output data directly affects control performance. The urgent need for low cost in commercial spaceflight has driven the application of low-cost GNSS components, but this has also significantly increased the risk of outliers in the output data. Therefore, how to effectively eliminate these outliers under the constraints of limited onboard computing power and real-time performance has become a critical problem that urgently needs to be solved. Existing GNSS outlier elimination methods mainly rely on judging the magnitude of position and velocity.

[0003] Existing GNSS outlier removal methods mainly use position and velocity thresholds for judgment, which can only roughly identify outliers with large errors. Furthermore, unreasonable threshold settings can prevent the correct GNSS data from being introduced. As satellites decay, the thresholds need to be updated regularly, resulting in low recognition rate and poor reliability, which cannot meet the current requirements for high-precision satellite control. Summary of the Invention

[0004] To address the issues of low identification rate and poor reliability in existing GNSS outlier removal methods, this invention provides a satellite-based GNSS outlier detection method based on residual multi-order chi-square distribution, comprising the following steps:

[0005] Step 1: Acquire GNSS output data, including position, velocity, and satellite time in the Earth-fixed system, where position and velocity are three-dimensional vectors;

[0006] Step 2: Convert the GNSS output data obtained in Step 1 into orbital square root counts;

[0007] Step 3: Calculate the residuals and residual change rates to determine the first field value. If the first field value passes, calculate the one-dimensional chi-square test field value. If the one-dimensional chi-square test passes, calculate the two-dimensional chi-square test field value. If the two-dimensional chi-square test passes, calculate the five-dimensional chi-square test field value. If all tests pass, update the parameters, including the state variables. X The coefficient of the first-order long-term term of the right ascension of the pseudo-elevation intersection. Argument of Latitude Long-term coefficient Argument of Latitude Second-order long-term coefficients If any test fails, the GNSS data captured at that time is judged as an outlier and marked as invalid. At the same time, the number of consecutive invalid data is counted. When the number of consecutive invalid data exceeds the threshold, the threshold will be adjusted as the number of consecutive invalid data increases to prevent the system from diverging.

[0008] As a further improvement of the present invention, in step two, the 14 square roots include: the relative amount of the semi-major axis. ; coefficients of the first-order long-term term of the quasi-flat semi-major axis ; Inclination of the proposed level track ; Proposed ascending intersection right ascension The coefficient of the first-order long-term term of the right ascension of the pseudo-elevation intersection. ; for Eccentricity component 1 at time. for The eccentricity component at time 2, where , This is the converted approximate eccentricity. The zero-order coefficient of the average latitude argument , , for The angle of near-perimeter at any given time, for The pseudo-mean perigee argument at time. The argument of the perigee; First-order long-term coefficients , The eccentricity component is 1; First-order long-term coefficients , The eccentricity component is 2; the coefficient of the first-order long-term term of the pseudo-planar perigee argument is... ; Long-term coefficient , Argument of latitude; Second-order long-term coefficients ; The relative star time cumulative seconds corresponding to the square root of 14 .

[0009] As a further improvement of the present invention, step three further includes:

[0010] Step 1: Select the state variables from the square root of track 14 ,in It is a relative quantity of the semi-major axis. To simulate the inclination angle of the track, To simulate the right ascension of the nodes, , , The zero-order coefficient of the argument of latitude. To simulate the argument of the perigee, The converted approximate eccentricity is T, where T represents the transpose.

[0011] Step 2: Calculate the state prediction, using the following formula:

[0012] ,

[0013] in Uh = 3.986004418e+14, where 3 is the Earth's gravitational constant. aE = 6378137, which is the average radius of the Earth. dt The difference in seconds between the predicted time and the current time. This indicates that the first element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the second element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the third element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the fourth element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the fifth element of the state variable X in step 1 is... k The predicted value at time +1; Represents the relative quantity of the semi-major axis at time k. The inclination of the pseudo-horizontal orbit at time k, The right ascension of the pseudo-ascending node at time k is denoted as . Represents the approximate eccentricity at time k. Indicates the mean latitude argument at time k;

[0014] Step 3: Calculate the covariance matrix prediction, the formula is as follows:

[0015] ,

[0016] in, express k+ The covariance matrix predicted at time 1. express k The state transition matrix at time t, express k The covariance matrix at time t, Represents the prediction noise matrix. T Indicates transpose;

[0017] Step 4: Use the corresponding parameters of the 14 square roots obtained in Step 2 as observations. The formula is as follows, which is the parameter extracted from the 14 roots obtained by converting the position, velocity, and satellite time from the GNSS output:

[0018] ;

[0019] Step 5: Calculate the residuals and the rate of change of the residuals, using the following formulas:

[0020] ,

[0021] in, express k The news at +1 hour, express k+ The observation at time 1, express k+ The predicted value of state variable X at time 1. express k+ 1-moment news and k The difference in information at any given moment, express k The new information of the moment;

[0022] Step 6: Determine whether the absolute values ​​of the residuals and residual change rates calculated in Step 5 are greater than the set threshold. If any element of the absolute value of the residuals and residual change rates is greater than the set threshold, it is determined to be an outlier, and Step 8 is executed.

[0023] Step 7: Calculate the chi-square value and make a judgment; the calculation matrix for the chi-square value is:

[0024] ,

[0025] in, Represents the new information covariance matrix. express k+ The covariance matrix predicted at time 1. Represents the observation noise matrix. Elements representing state variable X i The one-dimensional chi-square value, express k The news at +1 hour, The element identifier represents the state variable X. This represents the two-dimensional chi-square values ​​of elements 1 and 5 of the state variable X. This represents the five-dimensional chi-square value of the state variable X. Let element 2 represent the one-dimensional chi-square value of the state variable X. Let element 3 represent the one-dimensional chi-square value of the state variable X. The one-dimensional chi-square value of element 4 represents the state variable X;

[0026] If the formula

[0027] ,

[0028] If any one of the conditions is met, it is determined to be a wild value, and step 8 is executed, where, For threshold coefficient, The maximum threshold coefficient, For continuous invalid counting, select according to actual needs. To allow invalid counts, Elements of state variable X i The chi-square threshold, The coefficient of thermal expansion is 1 / 3. For the two-dimensional chi-square thresholds of elements 1 and 5, Let X be the five-dimensional chi-square coefficient of the state variable X; if all decisions pass, then perform a state update:

[0029] ,

[0030] use GNSS The calculation results are updated;

[0031]

[0032] in, This is the gain coefficient. It is a five-dimensional identity matrix. To predict the initial value of the noise matrix, This indicates that the state variable X in step 1 is in k The predicted value at time +1;

[0033] Step 8: Use the predicted values ​​of the state variable X and the covariance matrix P as the optimal estimates for the current time step, and increment the invalid count by 1. The formula is as follows:

[0034] ,

[0035] in, Represents the prediction noise matrix The coefficient of thermal expansion, This indicates that the state variable X in step 1 is in k The optimal estimate at time +1.

[0036] As a further improvement of the present invention, in step 7... chr_thr(i) The chi-square value corresponding to a significance level of α = 0.03 is 4.709. hello_thr The chi-square value with a significance level of α = 0.05 was 5.991. The chi-square value with a significance level of α=0.1 was 9.236; allowct was 100, chi2_factor was 1, and chi2_factor_Max was 5.

[0037] This invention also discloses a satellite-based GNSS outlier detection system based on residual multi-order chi-square distribution, comprising:

[0038] Parameter acquisition module: used to acquire GNSS output data, including position, velocity and satellite time in the Earth-fixed system, where position and velocity are three-dimensional vectors;

[0039] Parameter conversion module: used to convert the GNSS output data acquired by the parameter acquisition module into orbital square roots (14).

[0040] Outlier detection and judgment module: Used to calculate residuals and residual change rates for the first round of outlier judgment. If the first round passes, the module calculates the one-dimensional chi-square test outlier. If the one-dimensional chi-square test passes, the module calculates the two-dimensional chi-square test outlier. If the two-dimensional chi-square test passes, the module calculates the five-dimensional chi-square test outlier. If all tests pass, the module updates the parameters, including the state variables. X The coefficient of the first-order long-term term of the right ascension of the pseudo-elevation intersection. Argument of Latitude Long-term coefficient Argument of Latitude Second-order long-term coefficients If any test fails, the GNSS data captured at that time is judged as an outlier and marked as invalid. At the same time, the number of consecutive invalid data is counted. When the number of consecutive invalid data exceeds the threshold, the threshold will be adjusted as the number of consecutive invalid data increases to prevent the system from diverging.

[0041] As a further improvement of the present invention, in the parameter conversion module, the 14 square root count includes: the 14 square root count includes: the relative amount of the semi-major axis. ; coefficients of the first-order long-term term of the quasi-flat semi-major axis ; Inclination of the proposed level track ; Proposed ascending intersection right ascension The coefficient of the first-order long-term term of the right ascension of the pseudo-elevation intersection. ; for Eccentricity component 1 at time. for The eccentricity component at time 2, where , The converted approximate eccentricity, The zero-order coefficient of the average latitude argument , , for The angle of near-perimeter at any given time, for The pseudo-planar perigee argument at time. The argument of the perigee; First-order long-term coefficients , The eccentricity component is 1; First-order long-term coefficients , The eccentricity component is 2; the coefficient of the first-order long-term term of the pseudo-planar perigee argument is... ; Long-term coefficient , Argument of latitude; Second-order long-term coefficients ; The relative star time cumulative seconds corresponding to the square root of 14 .

[0042] As a further improvement of the present invention, the outlier detection and judgment module further includes:

[0043] Unit 1: Used to select state variables from the square root of orbit 14. ,in It is a relative quantity of the semi-major axis. To simulate the inclination angle of the track, To simulate the right ascension of the nodes, , , The zero-order coefficient of the argument of latitude. To simulate the argument of the perigee, The converted approximate eccentricity is T, where T represents the transpose.

[0044] The second unit is used to calculate state predictions, and its formula is as follows:

[0045] ,

[0046] in Uh = 3.986004418e+14, where 3 is the Earth's gravitational constant. aE = 6378137, which is the average radius of the Earth. dt The difference in seconds between the predicted time and the current time. This indicates that the first element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the second element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the third element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the fourth element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the fifth element of the state variable X in step 1 is... k The predicted value at time +1; Represents the relative quantity of the semi-major axis at time k. The inclination of the pseudo-horizontal orbit at time k, The right ascension of the pseudo-ascending node at time k is denoted as . Represents the approximate eccentricity at time k. Indicates the mean latitude argument at time k;

[0047] Unit 3: Used to calculate and predict the covariance matrix, the formula of which is as follows:

[0048] ,

[0049] in, express k+ The covariance matrix predicted at time 1. express k The state transition matrix at time t, express k The covariance matrix at time t, Represents the prediction noise matrix. T Indicates transpose;

[0050] The fourth unit is used to take the corresponding parameters of the 14 square roots obtained by the parameter conversion module as observations. The formula is as follows, which is the parameter extracted from the 14 roots obtained by converting the position, velocity, and satellite time from the GNSS output:

[0051] ;

[0052] Unit 5: Used to calculate residuals and the rate of change of residuals, the formulas are as follows:

[0053] ,

[0054] in, express k The news at +1 hour, express k+ The observation at time 1, express k+ The predicted value of state variable X at time 1. express k+ 1-moment news and k The difference in information at any given moment, express k The new information of the moment;

[0055] The sixth unit is used to determine whether the absolute values ​​of the residuals and residual change rates calculated by the fifth unit are greater than the set threshold. If any element of the absolute value of the residuals and residual change rates is greater than the set threshold, it is determined to be an outlier and the eighth unit is activated.

[0056] Unit 7: Used to calculate and evaluate chi-square values; the calculation matrix for chi-square values ​​is:

[0057] ,

[0058] in, Represents the new information covariance matrix. express k+ The covariance matrix predicted at time 1. Represents the observation noise matrix. Elements representing state variable X i The one-dimensional chi-square value, express k The news at +1 hour, The element identifier represents the state variable X. This represents the two-dimensional chi-square values ​​of elements 1 and 5 of the state variable X. This represents the five-dimensional chi-square value of the state variable X. Let element 2 represent the one-dimensional chi-square value of the state variable X. Let element 3 represent the one-dimensional chi-square value of the state variable X. The one-dimensional chi-square value of element 4 represents the state variable X;

[0059] If the formula

[0060] ,

[0061] If any one of the conditions is met, it is determined to be a wild value, and the eighth unit is activated. Here, chi2_factor_Max is the threshold coefficient, and chi2_factor_Max is the maximum threshold coefficient. ct For continuous invalid counting, select according to actual needs. allow To allow invalid counts, chr_thr(i) Elements of state variable X i The chi-square threshold, crrate The coefficient of thermal expansion is 1 / 3. hello_thr For the two-dimensional chi-square thresholds of elements 1 and 5, sense_thr Let X be the five-dimensional chi-square coefficient of the state variable X; if all decisions pass, then perform a state update:

[0062] ,

[0063] use GNSS The calculation results are updated;

[0064] ,

[0065] in, This is the gain coefficient. It is a five-dimensional identity matrix. To predict the initial value of the noise matrix, This indicates that the state variable X in step 1 is in k The predicted value at time +1;

[0066] Unit 8: Used to take the predicted values ​​of state variable X and covariance matrix P as the optimal estimate at the current time, incrementing the invalid count by 1. The formula is as follows:

[0067] ,

[0068] in, Represents the prediction noise matrix The coefficient of thermal expansion, This indicates that the state variable X in step 1 is in k The optimal estimate at time +1.

[0069] As a further improvement of the present invention, in the seventh unit... chr_thr(i) The chi-square value corresponding to a significance level of α = 0.03 is 4.709. hello_thr The chi-square value with a significance level of α = 0.05 was 5.991. The chi-square value with a significance level of α=0.1 was 9.236; allowct was 100, chi2_factor was 1, and chi2_factor_Max was 5.

[0070] The present invention also discloses a computer-readable storage medium storing a computer program configured to implement the steps of the method described in the present invention when invoked by a processor.

[0071] The beneficial effects of this invention are: in the application of low-cost GNSS components in commercial aerospace, the outlier removal method provided by this invention can significantly improve the outlier removal rate compared with existing methods, providing correct data for use on the satellite. This method has strong robustness and can significantly reduce satellite development costs while ensuring satellite control accuracy. Attached Figure Description

[0072] Figure 1 This is a flowchart of the GNSS outlier detection method of the present invention. Detailed Implementation

[0073] like Figure 1As shown, this invention discloses a GNSS outlier detection method based on residual multi-order chi-square distribution used on satellites, comprising the following steps:

[0074] Step 1: Acquire GNSS output data, including position, velocity, and satellite time in the Earth-fixed system, where position and velocity are three-dimensional vectors;

[0075] Step 2: Convert the GNSS output data obtained in Step 1 into orbital 14 square root numbers. The 14 square root numbers refer to a set of square root numbers containing 14 parameters, as shown in Table 1:

[0076] Table 1

[0077] ,

[0078] Step 3: Calculate the residuals and residual change rates to determine the first field value. If the first field value passes, calculate the one-dimensional chi-square test field value. If the one-dimensional chi-square test passes, calculate the two-dimensional chi-square test field value. If the two-dimensional chi-square test passes, calculate the five-dimensional chi-square test field value. If all tests pass, update the parameters, including the state variables. X The coefficient of the first-order long-term term of the right ascension of the pseudo-elevation intersection. Argument of Latitude Long-term coefficient Argument of Latitude Second-order long-term coefficients If any test fails, the GNSS data captured at that time is judged as an outlier and marked as invalid. At the same time, the number of consecutive invalid data is counted. When the number of consecutive invalid data exceeds the threshold, the threshold will be adjusted as the number of consecutive invalid data increases to prevent the system from diverging.

[0079] Step three also includes:

[0080] Step 1: Select the state variables from the square root of track 14 ,in It is a relative quantity of the semi-major axis. To simulate the inclination angle of the track, To simulate the right ascension of the nodes, The zero-order coefficient of the argument of latitude. To simulate the argument of the perigee, The converted approximate eccentricity is T, where T represents the transpose.

[0081] Step 2: Calculate the state prediction, using the following formula:

[0082] ,

[0083] in Uh = 3.986004418e+14, where 3 is the Earth's gravitational constant. aE = 6378137, which is the average radius of the Earth. dt The difference in seconds between the predicted time and the current time. This indicates that the first element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the second element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the third element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the fourth element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the fifth element of the state variable X in step 1 is... k The predicted value at time +1; Represents the relative quantity of the semi-major axis at time k. The inclination of the pseudo-horizontal orbit at time k, The right ascension of the pseudo-ascending node at time k is denoted as . Represents the approximate eccentricity at time k. Indicates the mean latitude argument at time k;

[0084] Step 3: Calculate the covariance matrix prediction, the formula is as follows:

[0085] ,

[0086] in, express k+ The covariance matrix predicted at time 1. express k The state transition matrix at time t, express k The covariance matrix at time t, Represents the prediction noise matrix. T Indicates transpose;

[0087] Step 4: Use the corresponding parameters of the 14 square roots obtained in Step 2 as observations. The formula is as follows, which is the parameter extracted from the 14 roots obtained by converting the position, velocity, and satellite time from the GNSS output:

[0088] ,

[0089] Step 5: Calculate the residuals and the rate of change of the residuals, using the following formulas:

[0090] ,

[0091] in, express kThe news at +1 hour, express k+ The observation at time 1, express k+ The predicted value of state variable X at time 1. express k+ 1-moment news and k The difference in information at any given moment, express k The new information of the moment;

[0092] Step 6: Determine whether the absolute values ​​of the residuals and residual change rates calculated in Step 5 are greater than the set threshold. If any element of the absolute value of the residuals and residual change rates is greater than the set threshold, it is determined to be an outlier, and Step 8 is executed.

[0093] Step 7: Calculate the chi-square value and make a judgment; the calculation matrix for the chi-square value is:

[0094] ,

[0095] in, Represents the new information covariance matrix. express k+ The covariance matrix predicted at time 1. Represents the observation noise matrix. Elements representing state variable X i The one-dimensional chi-square value, express k The news at +1 hour, The element identifier represents the state variable X. This represents the two-dimensional chi-square values ​​of elements 1 and 5 of the state variable X. This represents the five-dimensional chi-square value of the state variable X. Let element 2 represent the one-dimensional chi-square value of the state variable X. Let element 3 represent the one-dimensional chi-square value of the state variable X. The one-dimensional chi-square value of element 4 represents the state variable X;

[0096] If any condition in the formula is met, it is determined to be an outlier, and step 8 is executed, as follows:

[0097] First, calculate the new information covariance matrix. And one-dimensional chi-square value you see ( i ),

[0098] If ∃ exists you see ( i )> a ( i If the result is positive, proceed to step 8; otherwise, calculate the two-dimensional chi-square value. hello, If it exists hello > b If the result is positive, proceed to step 8; otherwise, calculate the five-dimensional chi-square value.

[0099] sense, If it exists sense > c If the condition is met, proceed to step 8; otherwise, if all checks pass, perform a state update. For threshold coefficient, The maximum threshold coefficient, For continuous invalid counting, select according to actual needs. To allow invalid counts, Elements of state variable X i The chi-square threshold, The coefficient of thermal expansion is 1 / 3. For the two-dimensional chi-square thresholds of elements 1 and 5, Let X be the five-dimensional chi-square coefficient of the state variable X; It is recommended to use a chi-square value of 4.709 corresponding to a significance level of α=0.03. It is recommended to use a chi-square value of 5.991 with a significance level of α=0.05. A chi-square value of 9.236 is recommended with a significance level of α = 0.1. `ct` represents continuous invalid counts, and `allowct` represents the allowed invalid counts; these should be selected according to actual needs, and are set to 100 in this invention. `chi2_factor` is the threshold coefficient, and it is recommended to set it to 1. The maximum coefficient is set to 5 in this invention, based on actual settings.

[0100] If all checks pass, the state is updated using the following formula:

[0101] ,

[0102] use GNSS The calculation results are updated;

[0103] ,

[0104] in, This is the gain coefficient. It is a five-dimensional identity matrix. To predict the initial value of the noise matrix, This indicates that the state variable X in step 1 is in k The predicted value at time +1;

[0105] Step 8: If it is determined to be an outlier, the predicted values ​​of the state variable X and the covariance matrix P are used as the optimal estimates at the current time, and the invalid count is incremented by 1. The formula is as follows:

[0106] ,

[0107] in, Represents the prediction noise matrix The coefficient of thermal expansion, This indicates that the state variable X in step 1 is... k The optimal estimate at time +1.

[0108] The working principle of the GNSS outlier detection method based on residual multi-order chi-square distribution disclosed in this invention is as follows: After acquiring GNSS data, it is converted into orbital square root (14). The residuals, residual change rate, one-dimensional chi-square distribution, two-dimensional chi-square distribution, and five-dimensional chi-square distribution are calculated. If any one of these conditions is not met, the GNSS data is considered invalid, and the predicted values ​​of X and P are taken as the optimal estimates for the current moment, with the invalidity count incremented by 1. Otherwise, the GNSS data is considered valid, the invalidity count is reset to zero, and the system is updated. .

[0109] This invention also discloses a satellite-based GNSS outlier detection system based on residual multi-order chi-square distribution, comprising:

[0110] Parameter acquisition module: used to acquire GNSS output data, including position, velocity and satellite time in the Earth-fixed system, where position and velocity are three-dimensional vectors;

[0111] Parameter conversion module: used to convert the GNSS output data acquired by the parameter acquisition module into orbital square roots (14).

[0112] Outlier detection and judgment module: Used to calculate residuals and residual change rates for the first round of outlier judgment. If the first round passes, the module calculates the one-dimensional chi-square test outlier. If the one-dimensional chi-square test passes, the module calculates the two-dimensional chi-square test outlier. If the two-dimensional chi-square test passes, the module calculates the five-dimensional chi-square test outlier. If all tests pass, the module updates the parameters, including the state variables. X The coefficient of the first-order long-term term of the right ascension of the pseudo-elevation intersection. Argument of Latitude Long-term coefficient Argument of Latitude Second-order long-term coefficients If any test fails, the GNSS data captured at that time is judged as an outlier and marked as invalid. At the same time, the number of consecutive invalid data is counted. When the number of consecutive invalid data exceeds the threshold, the threshold will be adjusted as the number of consecutive invalid data increases to prevent the system from diverging.

[0113] In the parameter conversion module, the number of square roots (14) includes:

[0114] The 14 square roots include: relative semi-major axis. ; coefficients of the first-order long-term term of the quasi-flat semi-major axis ; Inclination of the proposed level track ; Proposed ascending intersection right ascension The coefficient of the first-order long-term term of the right ascension of the pseudo-elevation intersection. ; for Eccentricity component 1 at time. for The eccentricity component at time 2, where , This is the converted approximate eccentricity. The zero-order coefficient of the average latitude argument , , for The angle of near-perimeter at any given time, for The pseudo-mean perigee argument at time. The argument of the perigee; First-order long-term coefficients , The eccentricity component is 1; First-order long-term coefficients , The eccentricity component is 2; the coefficient of the first-order long-term term of the pseudo-planar perigee argument is... ; Long-term coefficient , Argument of latitude; Second-order long-term coefficients ; The relative star time cumulative seconds corresponding to the square root of 14 .

[0115] The outlier detection and judgment module also includes:

[0116] Unit 1: Used to select state variables from the square root of orbit 14. ,in It is a relative quantity of the semi-major axis. To simulate the inclination angle of the track, To simulate the right ascension of the nodes, , , The zero-order coefficient of the argument of latitude. To simulate the argument of the perigee, The converted approximate eccentricity is T, where T represents the transpose.

[0117] The second unit is used to calculate state predictions, and its formula is as follows:

[0118] ,

[0119] in Uh = 3.986004418e+14, where 3 is the Earth's gravitational constant. aE = 6378137, which is the average radius of the Earth. dt The difference in seconds between the predicted time and the current time. This indicates that the first element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the second element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the third element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the fourth element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the fifth element of the state variable X in step 1 is... k The predicted value at time +1; Represents the relative quantity of the semi-major axis at time k. The inclination of the pseudo-horizontal orbit at time k, The right ascension of the pseudo-ascending node at time k is denoted as . Represents the approximate eccentricity at time k. Indicates the mean latitude argument at time k;

[0120] Unit 3: Used to calculate and predict the covariance matrix, the formula of which is as follows:

[0121] ,

[0122] in, express k+ The covariance matrix predicted at time 1. express k The state transition matrix at time t, express k The covariance matrix at time t, Represents the prediction noise matrix. T Indicates transpose;

[0123] The fourth unit is used to take the corresponding parameters of the 14 square roots obtained by the parameter conversion module as observations. The formula is as follows, which is the parameter extracted from the 14 roots obtained by converting the position, velocity, and satellite time from the GNSS output:

[0124] ;

[0125] Unit 5: Used to calculate residuals and the rate of change of residuals, the formulas are as follows:

[0126] ,

[0127] in, express k The news at +1 hour, express k+ The observation at time 1, express k+ The predicted value of state variable X at time 1. express k+ 1-moment news and k The difference in information at any given moment, express k The new information of the moment;

[0128] The sixth unit is used to determine whether the absolute values ​​of the residuals and residual change rates calculated by the fifth unit are greater than the set threshold. If any element of the absolute value of the residuals and residual change rates is greater than the set threshold, it is determined to be an outlier and the eighth unit is activated.

[0129] Unit 7: Used to calculate and evaluate chi-square values; the calculation matrix for chi-square values ​​is:

[0130] ,

[0131] in, Represents the new information covariance matrix. express k+ The covariance matrix predicted at time 1. Represents the observation noise matrix. Elements representing state variable X i The one-dimensional chi-square value, express k The news at +1 hour, The element identifier represents the state variable X. This represents the two-dimensional chi-square values ​​of elements 1 and 5 of the state variable X. This represents the five-dimensional chi-square value of the state variable X. Let element 2 represent the one-dimensional chi-square value of the state variable X. Let element 3 represent the one-dimensional chi-square value of the state variable X. The one-dimensional chi-square value of element 4 represents the state variable X;

[0132] If the formula

[0133] ,

[0134] If any one of the conditions is met, it is determined to be a wild value, and the eighth unit is activated. Here, chi2_factor_Max is the threshold coefficient, and chi2_factor_Max is the maximum threshold coefficient. ct For continuous invalid counting, select according to actual needs. allow To allow invalid counts, chr_thr(i) Elements of state variable X i The chi-square threshold, crrate The coefficient of thermal expansion is 1 / 3. hello_thr For the two-dimensional chi-square thresholds of elements 1 and 5, sense_thr Let X be the five-dimensional chi-square coefficient of the state variable X; if all decisions pass, then perform a state update:

[0135] ,

[0136] use GNSS The calculation results are updated;

[0137] ,

[0138] in, This is the gain coefficient. It is a five-dimensional identity matrix. To predict the initial value of the noise matrix, This indicates that the state variable X in step 1 is in k The predicted value at time +1;

[0139] Unit 8: Used to take the predicted values ​​of state variable X and covariance matrix P as the optimal estimate at the current time, incrementing the invalid count by 1. The formula is as follows:

[0140] ,

[0141] in, Represents the prediction noise matrix The coefficient of thermal expansion, This indicates that the state variable X in step 1 is in k The optimal estimate at time +1.

[0142] In Unit 7, chr_thr(i) The chi-square value corresponding to a significance level of α = 0.03 is 4.709. hello_thr The chi-square value with a significance level of α = 0.05 was 5.991. The chi-square value with a significance level of α=0.1 was 9.236; allowct was 100, chi2_factor was 1, and chi2_factor_Max was 5.

[0143] The present invention also discloses a computer-readable storage medium storing a computer program configured to implement the steps of the method described in the present invention when invoked by a processor.

[0144] The above description, in conjunction with specific preferred embodiments, provides a further detailed explanation of the present invention. It should not be construed that the specific implementation of the present invention is limited to these descriptions. For those skilled in the art, various simple deductions or substitutions can be made without departing from the concept of the present invention, and all such modifications and substitutions should be considered within the scope of protection of the present invention.

Claims

1. A method for GNSS outlier detection based on residual poly-chi-squared distribution used on-board satellites, characterized in that, Includes the following steps: Step 1: Acquire GNSS output data, including position, velocity, and satellite time in the Earth-fixed system, where position and velocity are three-dimensional vectors; Step 2: Convert the GNSS output data obtained in Step 1 into orbital square root counts; Step 3: Calculate the residuals and residual change rates to determine the first field value. If the first field value passes, calculate the one-dimensional chi-square test field value. If the one-dimensional chi-square test passes, calculate the two-dimensional chi-square test field value. If the two-dimensional chi-square test passes, calculate the five-dimensional chi-square test field value. If all tests pass, update the parameters, including the state variables. X The coefficient of the first-order long-term term of the right ascension of the pseudo-elevation intersection. Argument of Latitude Long-term coefficient Argument of Latitude Second-order long-term coefficients If any test fails, the GNSS data captured at that time is judged as an outlier and marked as invalid. At the same time, the number of consecutive invalidation is counted. When the number of consecutive invalidation exceeds a threshold, the threshold will be adjusted as the number of consecutive invalidation increases to prevent system divergence. In step two, the 14 square roots include: the relative amount of the semi-major axis. ; coefficients of the first-order long-term term of the quasi-flat semi-major axis ; Inclination of the proposed level track ; Proposed ascending intersection right ascension The coefficient of the first-order long-term term of the right ascension of the pseudo-elevation intersection. ; : The eccentricity component at time 1, where , for The pseudo-mean perigee argument at time. This is the converted approximate eccentricity. First-order long-term coefficients , The eccentricity component is 1; : Eccentricity component 2 at time. ; First-order long-term coefficients , The eccentricity component is 2; the coefficient of the first-order long-term term of the pseudo-planar perigee argument is... Argument of Latitude , , for The angle of near-perimeter at any given time, The argument of the perigee; Long-term coefficient , Argument of latitude; Second-order long-term coefficients ; The relative star time cumulative seconds corresponding to the square root of 14 .

2. The GNSS outlier detection method of claim 1, wherein, Step three also includes: Step 1: Select the state variables from the square root of track 14 ,in It is a relative quantity of the semi-major axis. To simulate the inclination angle of the track, To simulate the right ascension of the nodes, , , The zero-order coefficient of the argument of latitude. To simulate the argument of the perigee, The converted approximate eccentricity is T, where T represents the transpose. Step 2: Calculate the state prediction, using the following formula: , in uE = 3.986004418e+14, where 3 is the Earth's gravitational constant. aE = 6378137, which is the average radius of the Earth. dt The difference in seconds between the predicted time and the current time. This indicates that the first element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the second element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the third element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the fourth element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the fifth element of the state variable X in step 1 is... k The predicted value at time +1; Represents the relative quantity of the semi-major axis at time k. The inclination of the pseudo-horizontal orbit at time k, The right ascension of the pseudo-ascending node at time k is denoted as . Represents the approximate eccentricity at time k. Indicates the mean latitude argument at time k; Step 3: Calculate the covariance matrix prediction, the formula is as follows: , in, express k+ The covariance matrix predicted at time 1. express k The state transition matrix at time t, express k The covariance matrix at time t, Represents the prediction noise matrix. T Indicates transpose; Step 4: Use the corresponding parameters of the 14 square roots obtained in Step 2 as observations. The formula is as follows, which is the parameter extracted from the 14 roots obtained by converting the position, velocity, and satellite time from the GNSS output: , Step 5: Calculate the residuals and the rate of change of the residuals, using the following formulas: , in, express k The news at +1 hour, express k The observation at time +1, express k The predicted value of state variable X at time +1. express k +1 hour news and k The difference in information at any given moment, express k The new information of the moment; Step 6: Determine whether the absolute values ​​of the residuals and residual change rates calculated in Step 5 are greater than the set threshold. If any element of the absolute value of the residuals and residual change rates is greater than the set threshold, it is determined to be an outlier, and Step 8 is executed. Step 7: Calculate the chi-square value and make a judgment; the calculation matrix for the chi-square value is: , in, Represents the new information covariance matrix. express k The covariance matrix predicted at time +1 Represents the observation noise matrix. Elements representing state variable X i The one-dimensional chi-square value, express k The news at +1 hour, The element identifier represents the state variable X. This represents the two-dimensional chi-square values ​​of elements 1 and 5 of the state variable X. This represents the five-dimensional chi-square value of the state variable X. Let element 2 represent the one-dimensional chi-square value of the state variable X. Let element 3 represent the one-dimensional chi-square value of the state variable X. The one-dimensional chi-square value of element 4 represents the state variable X; If the formula , If any one of the conditions is met, it is determined to be a wild value, and step 8 is executed, where, For threshold coefficient, The maximum threshold coefficient, For continuous invalid counting, select according to actual needs. To allow invalid counts, Elements of state variable X i The chi-square threshold, The coefficient of thermal expansion is 1 / 3. For the two-dimensional chi-square thresholds of elements 1 and 5, Let X be the five-dimensional chi-square coefficient of the state variable X; if all decisions pass, then perform a state update: , use GNSS The calculation results are updated; ， in, This is the gain coefficient. It is a five-dimensional identity matrix. To predict the initial value of the noise matrix, This indicates that the state variable X in step 1 is in k The predicted value at time +1; Step 8: Use the predicted values ​​of the state variable X and the covariance matrix P as the optimal estimates for the current time step, and increment the invalid count by 1. The formula is as follows: , in, Represents the prediction noise matrix The coefficient of thermal expansion, This indicates that the state variable X in step 1 is in k The optimal estimate at time +1.

3. The GNSS outlier detection method according to claim 2, characterized in that, In step 7, chr_thr(i) The chi-square value corresponding to a significance level of α = 0.03 is 4.

709. chiau_thr The chi-square value with a significance level of α = 0.05 was 5.

991. The chi-square value with a significance level of α=0.1 was 9.236; allowct was 100, chi2_factor was 1, and chi2_factor_Max was 5.

4. A satellite-based GNSS outlier detection system based on residual multi-order chi-square distribution, characterized in that, include: Parameter acquisition module: used to acquire GNSS output data, including position, velocity and satellite time in the Earth-fixed system, where position and velocity are three-dimensional vectors; Parameter conversion module: used to convert the GNSS output data acquired by the parameter acquisition module into orbital square roots (14). Outlier detection and judgment module: Used to calculate residuals and residual change rates for the first round of outlier judgment. If the first round passes, the module calculates the one-dimensional chi-square test outlier. If the one-dimensional chi-square test passes, the module calculates the two-dimensional chi-square test outlier. If the two-dimensional chi-square test passes, the module calculates the five-dimensional chi-square test outlier. If all tests pass, the module updates the parameters, including the state variables. X The coefficient of the first-order long-term term of the right ascension of the pseudo-elevation intersection. Argument of Latitude Long-term coefficient Argument of Latitude Second-order long-term coefficients If any test fails, the GNSS data captured at that time is judged as an outlier and marked as invalid. At the same time, the number of consecutive invalidation is counted. When the number of consecutive invalidation exceeds a threshold, the threshold will be adjusted as the number of consecutive invalidation increases to prevent system divergence. In the parameter conversion module, the 14 square roots include: relative quantity of semi-major axis. ; coefficients of the first-order long-term term of the quasi-flat semi-major axis ; Inclination of the proposed level track ; Proposed ascending intersection right ascension The coefficient of the first-order long-term term of the right ascension of the pseudo-elevation intersection. ; : The eccentricity component at time 1, where , for The pseudo-mean perigee argument at time. This is the converted approximate eccentricity. First-order long-term coefficients , The eccentricity component is 1; : Eccentricity component 2 at time. ; First-order long-term coefficients , The eccentricity component is 2; the coefficient of the first-order long-term term of the pseudo-planar perigee argument is... Argument of Latitude , , for The angle of near-perimeter at any given time, The argument of the perigee; Long-term coefficient , Argument of latitude; Second-order long-term coefficients ; The relative star time cumulative seconds corresponding to the square root of 14 .

5. The GNSS outlier detection system according to claim 4, characterized in that, The outlier detection and judgment module also includes: Unit 1: Used to select state variables from the square root of orbit 14. ,in It is a relative quantity of the semi-major axis. To simulate the inclination angle of the track, To simulate the right ascension of the nodes, , , The zero-order coefficient of the argument of latitude. To simulate the argument of the perigee, The converted approximate eccentricity is T, where T represents the transpose. The second unit is used to calculate state predictions, and its formula is as follows: ， in uE = 3.986004418e+14, where 3 is the Earth's gravitational constant. aE = 6378137, which is the average radius of the Earth. dt The difference in seconds between the predicted time and the current time. This indicates that the first element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the second element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the third element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the fourth element of the state variable X in step 1 is... k The predicted value at time +1, This indicates that the fifth element of the state variable X in step 1 is... k The predicted value at time +1; Represents the relative quantity of the semi-major axis at time k. Indicates the pseudo-horizontal orbital inclination at time k. The right ascension of the pseudo-ascending node at time k is denoted as . Represents the approximate eccentricity at time k. Indicates the mean latitude argument at time k; Unit 3: Used to calculate and predict the covariance matrix, the formula of which is as follows: ， in, express k+ The covariance matrix predicted at time 1. express k The state transition matrix at time t, express k The covariance matrix at time t, Represents the prediction noise matrix. T Indicates transpose; The fourth unit is used to take the corresponding parameters of the 14 square roots obtained by the parameter conversion module as observations. The formula is as follows, which is the parameter extracted from the 14 roots obtained by converting the position, velocity, and satellite time from the GNSS output: , Unit 5: Used to calculate residuals and the rate of change of residuals, the formulas are as follows: ， in, express k The news at +1 hour, express k The observation at time +1, express k The predicted value of state variable X at time +1. express k +1 hour news and k The difference in information at any given moment, express k The new information of the moment; The sixth unit is used to determine whether the absolute values ​​of the residuals and residual change rates calculated by the fifth unit are greater than the set threshold. If any element of the absolute value of the residuals and residual change rates is greater than the set threshold, it is determined to be an outlier and the eighth unit is activated. Unit 7: Used to calculate and evaluate chi-square values; the calculation matrix for chi-square values ​​is: ， in, Represents the new information covariance matrix. express k+ The covariance matrix predicted at time 1. Represents the observation noise matrix. Elements representing state variable X i The one-dimensional chi-square value, express k The news at +1 hour, The element identifier represents the state variable X. This represents the two-dimensional chi-square values ​​of elements 1 and 5 of the state variable X. This represents the five-dimensional chi-square value of the state variable X. Let element 2 represent the one-dimensional chi-square value of the state variable X. Let element 3 represent the one-dimensional chi-square value of the state variable X. The one-dimensional chi-square value of element 4 represents the state variable X; If the formula ， If any one of the conditions is met, it is determined to be a wild value, and the eighth unit is activated. Here, chi2_factor_Max is the threshold coefficient, and chi2_factor_Max is the maximum threshold coefficient. ct For continuous invalid counting, select according to actual needs. allowct To allow invalid counts, chr_thr(i) Elements of state variable X i The chi-square threshold, chrrate The coefficient of thermal expansion is 1 / 3. chiau_thr For the two-dimensional chi-square thresholds of elements 1 and 5, chiall_thr Let X be the five-dimensional chi-square coefficient of the state variable X; if all decisions pass, then perform a state update: , use GNSS The calculation results are updated; ， in, This is the gain coefficient. It is a five-dimensional identity matrix. To predict the initial value of the noise matrix, This indicates that the state variable X in step 1 is in k The predicted value at time +1; Unit 8: Used to take the predicted values ​​of state variable X and covariance matrix P as the optimal estimate at the current time, incrementing the invalid count by 1. The formula is as follows: ， in, Represents the prediction noise matrix The coefficient of thermal expansion, This indicates that the state variable X in step 1 is in k The optimal estimate at time +1.

6. The GNSS outlier detection system according to claim 5, characterized in that, In the seventh unit, chr_thr (i) The chi-square value corresponding to a significance level of α = 0.03 is 4.

709. chiau_thr The chi-square value with a significance level of α = 0.05 was 5.

991. The chi-square value with a significance level of α=0.1 was 9.236; allowct was 100, chi2_factor was 1, and chi2_factor_Max was 5.

7. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores a computer program configured to implement the steps of the method according to any one of claims 1-3 when invoked by a processor.

Images