A method for calculating the protection level of GNSS positioning state domain

Abstract

A method for calculating the protection level of a GNSS positioning state domain, which performs fault detection and elimination by constructing state domain fault detection statistics and detection thresholds, and calculates the protection level caused by positioning deviation based on the false alarm rate and the protection level caused by positioning noise based on the missed detection rate, and finally determines the final protection level of the current epoch, and issues an alarm to the user when the protection level exceeds the alarm limit of a specific application field; when the protection level is less than or equal to the alarm limit of a specific application field, the protection level calculation for the next epoch is performed. The present invention calculates the protection level from the state domain, avoids the indirectness of the measurement residual relative to the state domain protection level, makes full use of the state estimation historical data, obtains the optimal protection level on the premise of meeting the integrity risk and continuity risk requirements, reflects the real navigation positioning error, and improves the reliability and calculation efficiency of the protection level.

Description

Technical Field

[0001] The invention relates to a method for calculating a protection level of a GNSS positioning state domain, and belongs to the technical field of satellite positioning protection level calculation. Background Art

[0002] With the development of global satellite navigation system (GNSS) positioning technology, the integration of multiple constellations and the increase in the number of satellites have significantly improved the accuracy and robustness of GNSS positioning, and also made it possible for receiver autonomous integrity monitoring (RAIM) to meet the navigation performance requirements of various fields on a global scale. RAIM uses consistency detection theory to detect and eliminate faults at the user end, and quantifies the integrity risks caused by undetectable faults and incorrect elimination of faulty satellites in the form of protection levels. In practical applications, RAIM is based on the integrity risk and continuity risk requirements of professional fields (horizontal / vertical protection level, H / VPL). When H / VPL is less than the horizontal / vertical alarm limit (horizontal / vertical alarm limit, H / VAL), it indicates that the current RAIM result is available, otherwise it is unavailable. Therefore, H / VPL is an important output parameter of RAIM.

[0003] H / VPL is not only a function of the random error model, fault model and satellite geometry, but also the construction method of the test statistic will affect the calculation of H / VPL. There are two main ways to construct fault detection statistics in RAIM, one is the detection statistic of the multiple hypothesis solution separation (MHSS) method, and the other is the residual chi-square test method (RCTM) statistic based on the measurement residual. Among them, MHSS divides the Kalman filter into three levels, and performs multiple fault detection and elimination through the difference of the state estimation value of the main filter, sub-filter and secondary sub-filter. Since the integrity risk is defined in the position domain, MHSS can better quantify the integrity risk and then calculate H / VPL; but MHSS calculation is complex, the amount of calculation will increase exponentially with the increase of visible satellites, and it cannot cope with slowly varying slope faults. RCTM detects and eliminates multiple faults through new information. The threshold of new information detection is relatively large and the protection level of the relative position domain is indirect, so H / VPL is relatively large. Autonomous integrity monitoring with an extrapolation method (AIME) takes into account the possible impact of slow-varying ramp faults on positioning results, reconstructs the calculation formula of H / VPL, and obtains a protection level smaller than MHSS in simulation data. However, AIME assumes that the number of satellites in each epoch is the same when it is designed, which greatly limits its application in practice. Therefore, under the premise of meeting the requirements of integrity risk and continuity risk, how to make full use of historical status information to obtain the optimal protection level is still a problem that needs to be solved urgently. Summary of the invention

[0004] The present invention provides a method for calculating the protection level of a GNSS positioning state domain. The method can carry out calculations from the state domain and fully utilize historical state information to obtain an optimal protection level while meeting the requirements of integrity risk and continuity risk.

[0005] In order to achieve the above object, the present invention provides a method for calculating the protection level of a GNSS positioning state domain, comprising the following steps:

[0006] S1: Determine the recursive period m based on the GNSS positioning random model. If m>k, where k is the current epoch number, obtain the posterior state estimate of the current epoch through Kalman filtering. and its covariance matrix P xx,k , calculate the prior state estimate and the posterior state estimate The difference between the two k and its covariance matrix Pdd,k ; If m≤k, use the state recursor execution time update to obtain the m recursive state estimates within the recursive cycle and its covariance matrix At the same time, the extended Kalman filter is used to obtain m a posteriori state estimates within the recursive cycle. and its covariance matrix P xx,k-m+i , calculate the difference between the recursive state estimate and the posterior state estimate within the recursive cycle and its covariance matrix The state recursor obtains the recursive state estimate and its covariance matrix, which means updating the posterior state estimate and its covariance matrix of km epoch to k epoch in time;

[0007] S2: When m>k, use the d of the current epoch k and its covariance matrix P dd,k Perform fault detection and elimination; when m≤k, calculate the difference between the recursive state estimate and the posterior state estimate within the recursive cycle The weighted average and its covariance matrix Construct state domain fault detection statistics and detection thresholds to perform fault detection and elimination;

[0008] S3: Based on the false alarm rate P FA Calculate the protection level PL1 caused by positioning deviation;

[0009] S4: Based on the missed detection rate P MD Calculate the protection level PL2 caused by positioning noise;

[0010] S5: Determine the final protection level PL of the current epoch, and issue an alarm to the user when the protection level PL exceeds the alarm limit AL of the specific application field; when the protection level PL is less than or equal to the alarm limit AL of the specific application field, calculate the protection level of the next epoch.

[0011] Furthermore, in step S1, when m>k, P dd,k =P x ' x,k -P xx,k , where P x ' x,k is the prior state estimate covariance matrix of the current epoch; when m≤k, the recursive state estimate Recursive state estimation covariance matrix

[0012] Among them, the initial recursive state estimate Its covariance matrix f(·) is the state transition function, Indicates that the state vector x obeys the mean Covariance matrix cov[x] = P x ' x Gaussian distribution, Q is the system noise covariance matrix, the difference between the recursive state estimate and the posterior state estimate The covariance matrix of

[0013] Furthermore, when m>k in step S2, the fault detection statistic is The detection threshold is in, P dd,k The eigenvector matrix of dd,k P dd,k The eigenvalue matrix of For dd,k The maximum eigenvalue of The i in the equation is the component corresponding to the maximum eigenvalue, Φ -1 (·) is the inverse function of the normal cumulative distribution function, P FA is the false alarm rate based on continuous risk allocation; when m≤k, within the recursive cycle The weighted average The covariance matrix of The fault detection statistic is The detection threshold is in for The eigenvector matrix of for The eigenvalue matrix of for The maximum eigenvalue of k and β k If any one of the above is greater than the threshold, it means that there is a fault. After eliminating the fault, that is, when the detection statistic α k and β k If both are less than or equal to the threshold, re-enter step S1. k and β k If both are less than or equal to the threshold, it means there is no fault, and the process goes to step S3;

[0014] Furthermore, the protection level PL1 caused by the positioning deviation in step S3 is divided into a horizontal protection level HPL1 and a vertical protection level VPL1. When only the positioning deviation of the current epoch is considered, in, and and They are Λdd,k The eigenvalue components in the east, north and sky directions; when considering the positioning deviation within the recursive period, we have in, and and They are Eigenvalue components in the east, north, and sky directions.

[0015] Furthermore, in step S4, the protection level PL2 caused by the positioning noise is divided into a horizontal protection level HPL2 and a vertical protection level VPL2. in, and and P xx,k In the three directions of east, north and sky, P MD is the missed detection rate based on integrity risk allocation.

[0016] Furthermore, in step S5, the final protection level PL of the current epoch is divided into a horizontal protection level HPL and a vertical protection level VPL. Considering only the positioning deviation of the current epoch, HPL = HPL1 + HPL2, VPL = VPL1 + VPL2; considering the positioning deviation within the recursive period, The alarm limit AL is divided into the horizontal alarm limit HAL and the vertical alarm limit VAL; when the horizontal protection level HPL is greater than the horizontal alarm limit HAL or the vertical protection level VPL is greater than the vertical alarm limit VAL, an alarm is issued to the user.

[0017] Furthermore, in step S2, the process of fault detection and troubleshooting is as follows:

[0018] (1) Rearrange the GNSS observations of the current faulty epoch from small to large according to the normalized new information, and take the first n satellite observations to calculate the posterior state estimate and the covariance matrix P 0,xx,k , and based on the posterior state estimate and the covariance matrix P 0,xx,k Construct a fault detection statistic and compare it with the detection threshold. When the fault detection statistic is greater than the detection threshold, it means that the quality of the observation data of the current epoch is poor and the fault cannot be eliminated. An alarm is issued to the user and the fault elimination of the current epoch is terminated. When the fault detection statistic is less than or equal to the detection threshold, it means that the fault can be eliminated and proceed to the next step.

[0019] (2) Delete the GNSS observation corresponding to the maximum innovation and use the remaining satellite observations to calculate the posterior state estimate and the covariance matrix P test,xx,k , the posterior state estimate and the covariance matrix P test,xx,k Construct a fault detection statistic and compare it with the detection threshold. When the fault detection statistic is greater than the detection threshold, it means that if a fault exists, it means that the fault has not been completely eliminated. Delete the GNSS observation value corresponding to the maximum new information in the remaining observations, and repeat this step until the fault is eliminated. Perform fault detection and elimination for the next epoch. When the fault detection statistic is less than or equal to the detection threshold, it means that the fault has been eliminated, and directly perform fault detection and elimination for the next epoch.

[0020] The present invention calculates the protection level directly from the state domain of the GNSS system, constructs state domain fault detection statistics and detection thresholds for fault detection and elimination, calculates the protection level caused by positioning deviation based on the false alarm rate, and calculates the protection level caused by positioning noise based on the missed detection rate, and finally determines the final protection level of the current epoch. Compared with the residual chi-square test method, the indirectness of the measured residual relative to the state domain protection level is avoided. Compared with the multi-hypothesis solution separation method that can only use the current epoch data to calculate the protection level, the method of the present invention makes full use of the state estimation historical data, and obtains the optimal protection level on the premise of meeting the integrity risk and continuity risk requirements, reflecting the real navigation positioning error, and improving the reliability and calculation efficiency of the protection level. In addition, the present invention can not only be used for single systems such as GPS, GLONASS, BDS and GALILEO, but also can be used for combined positioning of different satellite systems and inertial systems, and has a wide range of applications. BRIEF DESCRIPTION OF THE DRAWINGS

[0021] Figure 1 is a flow chart of the present invention;

[0022] Figure 2 It is the positioning error diagram before and after the SPP mutation fault simulation data experiment fault elimination;

[0023] Figure 3 It is a comparison chart of protection level of SPP sudden fault simulation data experimental level;

[0024] Figure 4 This is a comparison chart of vertical protection levels in the SPP sudden fault simulation data experiment;

[0025] Figure 5 This is a comparison chart of the protection level of the SPP slow-changing fault simulation data experiment level;

[0026] Figure 6 This is a comparison chart of vertical protection levels in the SPP slow-changing fault simulation data experiment. DETAILED DESCRIPTION

[0027] The present invention will be further described below in conjunction with the accompanying drawings.

[0028] like Figure 1 As shown, a method for calculating the protection level of a GNSS positioning state domain comprises the following steps:

[0029] S1: Determine the recursive period m based on the GNSS positioning random model. If m>k, where k is the current epoch number, obtain the posterior state estimate of the current epoch through Kalman filtering. and its covariance matrix P xx,k , calculate the prior state estimate and the posterior state estimate The difference between the two k and its covariance matrix P dd,k ; If m≤k, use the state recursor execution time update to obtain the m recursive state estimates within the recursive cycle and its covariance matrix At the same time, the extended Kalman filter is used to obtain m a posteriori state estimates within the recursive cycle. and its covariance matrix P xx,k-m+i , calculate the difference between the recursive state estimate and the posterior state estimate within the recursive cycle and its covariance matrix

[0030] S2: When m>k, use the d of the current epoch k and its covariance matrix P dd,k Perform fault detection and elimination; when m≤k, calculate the difference between the recursive state estimate and the posterior state estimate within the recursive cycle The weighted average and its covariance matrix Construct state domain fault detection statistics and detection thresholds to perform fault detection and elimination;

[0031] S3: Based on the false alarm rate P FA Calculate the protection level PL1 caused by positioning deviation;

[0032] S4: Based on the missed detection rate P MD Calculate the protection level PL2 caused by positioning noise;

[0033] S5: Determine the final protection level PL of the current epoch, and issue an alarm to the user when the protection level PL exceeds the alarm limit AL of the specific application field; when the protection level PL is less than or equal to the alarm limit AL of the specific application field, calculate the protection level of the next epoch.

[0034] The random model in GNSS positioning is affected by multiple factors such as the positioning model, observation environment, and data quality, and has uncertainty. The recursive period m depends on the specific situation of the noise size in the random model, and its value range is recommended to be 5 to 10 seconds. If it is set too small, the state information in a short period of time cannot detect slow-changing faults. If it is set too large, the noise accumulated over a long period of time in the random model will mask the impact of slow-changing faults, because the state recurser only updates time without measurement updates, and the state noise will be superimposed over time.

[0035] Embodiment: The GNSS positioning state domain multi-fault elimination method of the present invention can be used not only for single systems such as GPS, GLONASS, BDS and GALILEO, but also for combined positioning of different satellite systems and inertial systems. The following embodiment is the application of the present invention in BDS pseudorange single point positioning, which specifically includes the following steps:

[0036] S1: Determine the recursive period m. If m>k, obtain the posterior state estimate of the current epoch through the BDS pseudo-range single point positioning extended Kalman filter. and its covariance matrix P xx,k , calculate the prior state estimate and the posterior state estimate The difference between the two k and its covariance matrix P dd,k ; If m≤k, use the state recursion device to obtain the m recursive state estimates within the recursive cycle and its covariance matrix At the same time, the extended Kalman filter is used to obtain m a posteriori state estimates within the recursive cycle. and its covariance matrix P xx,k-m+i , calculate the difference between the recursive state estimate and the posterior state estimate within the recursive cycle and its covariance matrix

[0037] S1.1: Define the receiver initial state value and false alarm rate

[0038] Define the a priori state estimate of the receiver at epoch k and its covariance matrix is ​​P x ' x,k , And determine the false alarm rate as P FA ,in Indicates that the state vector x obeys the mean Covariance matrix cov[x] = P x ' x Gaussian distribution, let k = 0, then the system is in the initial state;

[0039] S1.2: Get the pseudo-range observation value of the BDS system:

[0040] P=ρ-cδt+I+T+ε (1)

[0041] Where P is the pseudorange, ρ is the geometric distance from the satellite to the receiver, cδt is the speed of light multiplied by the receiver clock error, I is the ionospheric delay error, T is the tropospheric delay error, and ε is the pseudorange observation noise. Formula (1) is linearized at the approximate coordinates of the station:

[0042] P-ρ0-D=ldX+mdY+ndZ+cδt (2)

[0043] Where ρ0 is the geometric distance from the satellite to the approximate coordinates of the receiver, D = I + T + ε, (l, m, n) is the direction cosine from the receiver to the satellite, and (dX, dY, dZ) is the coordinate correction number;

[0044] S1.3: Posterior state estimate for k epochs via BDS extended Kalman filter and its covariance matrix P xx,k :

[0045] S1.3.1: Define the state equation and measurement equation of the BDS system:

[0046] x k+1 =f k (x k )+w k (3)

[0047] z k =h k (x k )+v k (4)

[0048] In the formula, x k is the state vector, including the receiver coordinate correction and clock correction; f k is the state transfer function; w k is the system process noise, and its covariance matrix is ​​Q k ; z k is the BDS pseudorange measurement value in equation (2), h k is the measurement function, v k is the measurement noise, and its covariance matrix is ​​R k ;

[0049] S1.3.2: Extended Kalman filtering is performed on the BDS to obtain the posterior state estimate of k epochs and its covariance matrix P xx,k :

[0050]

[0051] P xx,k =P x ' x,k -K k P z ' z,k (K k ) T (6)

[0052] in:

[0053] K k = P x ' z,k (P z ' z,k ) -1 (7)

[0054] P xx,k =P x ' x,k -K k P z ' z,k (K k ) T (8)

[0055]

[0056]

[0057] In the formula, K k is the Kalman filter gain, is the prior measurement value, P x ' z,k is the state measurement cross covariance matrix, P z ' z,k is the measurement covariance matrix;

[0058] S1.4: When m>k, construct the fault detection residual d of the current epoch k and its covariance matrix P dd,k :

[0059]

[0060] P dd,k =P x ' x,k -P xx,k (13)

[0061] When m≤k, construct the fault detection residuals for each epoch in the recursive period and its covariance matrix

[0062]

[0063] The recursive state estimate and its covariance matrix for:

[0064]

[0065] Initial recursive state estimate Its covariance matrix

[0066] S2: If m≤k, calculate the weighted average of the difference between the state recursive estimate and the posterior state estimate within the recursive cycle and its covariance matrix Construct state domain fault detection statistics and detection thresholds to perform fault detection and troubleshooting:

[0067] S2.1: When m≤k, the fault detection residual within the recursive cycle The weighted average and its covariance matrix for:

[0068]

[0069] S2.2: When m>k, the fault detection statistic α k With detection threshold They are:

[0070]

[0071]

[0072] in, P dd,k The eigenvector matrix of dd,k P dd,k The eigenvalue matrix of For dd,k The maximum eigenvalue of The i in the equation is the component corresponding to the maximum eigenvalue, Φ -1 (·) is the inverse function of the normal cumulative distribution function, P FA is the false alarm rate based on the continuous risk allocation;

[0073] When m≤k, the fault detection statistic β k With detection threshold They are:

[0074]

[0075] in, for The eigenvector matrix of for The eigenvalue matrix of for The maximum eigenvalue of

[0076] S2.3: If the detection statistic is greater than the threshold, it means there is a fault. After eliminating the fault, go back to step S1.4. If the detection statistic is less than the threshold, it means there is no fault, and go to step S3;

[0077] S3: Based on the false alarm rate P FA Calculate the protection level PL1 caused by the positioning deviation. PL1 includes the horizontal protection level HPL1 and the vertical protection level VPL1. If only the positioning deviation of the current epoch is considered, then

[0078]

[0079] in and and They are Λ dd,k The eigenvalue components in the east, north, and sky (ENU) directions;

[0080] If the positioning deviation within the recursive cycle is considered, then

[0081]

[0082] in and and They are The eigenvalue components in the east, north, and sky (ENU) directions;

[0083] S4: Based on the missed detection rate P MD Calculate the protection level PL2 caused by positioning noise. PL2 includes the horizontal protection level HPL2 and the vertical protection level VPL2. Since the noise comparison deviation will not be accumulated over multiple epochs,

[0084]

[0085] in and and P xx,k In the three directional components of the northeast sky (ENU), P MD is the missed detection rate based on integrity risk allocation;

[0086] S5: Determine the final protection level PL of the current epoch, and issue an alarm to the user when the protection level PL exceeds the alarm limit AL of the specific application field; when the protection level PL is less than or equal to the alarm limit AL of the specific application field, calculate the protection level of the next epoch:

[0087] S5.1: The final protection level PL of the current epoch is divided into the horizontal protection level HPL and the vertical protection level VPL. If only the positioning deviation of the current epoch is considered, then the comprehensive equations (24), (25), (28) and (29), HPL and VPL are:

[0088] HPL=HPL1+HPL2 (30)

[0089] VPL=VPL1+VPL2 (31)

[0090] If the positioning error within the recursive cycle is considered, then the HPL and VPL are respectively:

[0091]

[0092] S5.2: If the protection level exceeds the alarm limit AL for a specific application area, an alarm is issued to the user. Otherwise, the protection level calculation for the next epoch is continued, and the prior state estimate for the next epoch is predicted by the following formula: and its covariance matrix P x ' x,k+1 :

[0093]

[0094] Let k=k+1 and return to step S1.2.

[0095] To further illustrate the characteristics and advantages of the present invention, SPP simulation data experiments are used to compare and verify the protection levels given by different methods after fault detection and elimination, including MHSS (multi-hypothesis solution separation method), RCTM (residual chi-square test method), SRCTM (state domain robust chi-square test method), AIME (autonomous integrity extrapolation method) and SRAIME (state domain robust autonomous integrity extrapolation method). Among them, SRCTM is the method of calculating the protection level using the current epoch data of the present invention, and SRAIME is the method of calculating the protection level using the data within the recursive cycle of the present invention. For the convenience of comparison, according to the integrity index required by LPV-200, the integrity risk in the experiment is 10 -7 , the continuity risk is 4*10 -6 , and evenly distributed in both vertical and horizontal directions and in each fault hypothesis, with the horizontal alarm limit being 40m and the vertical alarm limit being 35m;

[0096] This experiment uses a modified version of the GPSoft toolbox to generate GPS single-frequency pseudorange observations with a standard deviation of 1 m and pseudorange increments with a standard deviation of 0.1 m / s for simulation experiments. The reference trajectory is directly derived from the pre-designed mathematical model. Figure 2 The reference trajectory and estimated trajectory of the SPP simulation data experiment are as follows: the sampling frequency of the receiver is 1 Hz, the dimension of the measurement domain is related to the number of visible satellites and is constant at 10, and the recursion period of AIME and SRAIME is 5 s;

[0097] Since there is no fault in the simulated observations in the simulation data, the positioning result has high accuracy. In order to compare the differences between different protection level calculation methods under the influence of different faults, sudden fault and slow fault schemes are set:

[0098] ① Artificially add a 20m mutation fault to a pseudorange observation value at 100s, 300s and 500s respectively;

[0099] ② Slow-varying faults with different growth rates (v1 = 0.0005*(k-200), v2 = 0.00075*(k-200), v3 = 0.001*(k-200)) are artificially added to the receiver clock drift in 200-400s;

[0100] from Figure 2 It can be seen that after adopting the fault detection and elimination algorithm, the added faults are successfully eliminated and the positioning accuracy is improved;

[0101] Figure 3 and Figure 4 They are H / VAL calculated from the SPP sudden fault simulation data experiments of different methods, where H / VPE represents the horizontal / vertical positioning true error. It can be seen that in a short period of time after the fault is eliminated, the protection levels of the five methods fluctuate to varying degrees and then quickly converge to the normal level. This is because the reduction in the fault epoch observation value causes the calculated posterior state covariance matrix to become larger, and then the posterior state covariance matrix converges to the normal level after the Kalman filter. In addition, it can be found through observation that RCTM has the largest protection level, because the detection threshold of RCTM is relatively large and the protection level in the relative position domain is indirect, while the protection level of SRCTM is greater than MHSS, because the residual error of SRCTM fault detection is greater than MHSS. After reconstructing the protection levels of RCTM and SRCTM taking into account the influence of slow-varying slope faults, it can be seen that the protection level of AIME is less than MHSS, while the protection level of SRAIME is the smallest;

[0102] Figure 5 and Figure 6The H / VPL under slow-changing faults with different growth rates are given respectively. It can be seen that the added slow-changing fault has little impact on the horizontal direction, and the positioning error almost does not exceed the protection level. However, the vertical direction is greatly affected by the slow-changing fault, and it exceeds the protection level after a certain period of growth. The faster the growth rate, the earlier the positioning error exceeds the protection level. Since SRAIME has the smallest vertical protection level, it can detect the slow-changing slope fault first, which reflects the superiority of SRAIME.

Claims

1. A method for calculating the protection level of a GNSS positioning state domain, characterized in that: The steps include: S1: Determine the recursive period m based on the GNSS positioning random model. If m>k, where k is the current epoch number, obtain the posterior state estimate of the current epoch through Kalman filtering. and its covariance matrix P xx,k , calculate the prior state estimate and the posterior state estimate The difference between the two k and its covariance matrix P dd,k ; If m≤k, use the state recursor execution time update to obtain the m recursive state estimates within the recursive cycle and its covariance matrix At the same time, the extended Kalman filter is used to obtain m a posteriori state estimates within the recursive cycle. and its covariance matrix P xx,k-m+i , calculate the difference between the recursive state estimate and the posterior state estimate within the recursive cycle and its covariance matrix S2: When m>k, use the d of the current epoch k and its covariance matrix P dd,k Perform fault detection and elimination; when m≤k, calculate the difference between the recursive state estimate and the posterior state estimate within the recursive cycle The weighted average and its covariance matrix Construct state domain fault detection statistics and detection thresholds to perform fault detection and elimination. The specific process is as follows: When m>k, the fault detection statistic is The detection threshold is in, P dd,k The eigenvector matrix of dd,k P dd,k The eigenvalue matrix of For dd,k The maximum eigenvalue of The i in the equation is the component corresponding to the maximum eigenvalue, Φ -1 (·) is the inverse function of the normal cumulative distribution function, P FA is the false alarm rate based on the continuous risk allocation; When m≤k, within the recursive cycle The weighted average The covariance matrix of The fault detection statistic is The detection threshold is in for The eigenvector matrix of for The eigenvalue matrix of for The maximum eigenvalue of k and β k If any one of the above is greater than the threshold, it means that there is a fault. After eliminating the fault, that is, when the detection statistic α k and β k are less than or equal to the threshold, re-enter step S1. If the detection statistic α k and β k If both are less than or equal to the threshold, it means there is no fault, and the process goes to step S3; S3: Based on the false alarm rate P FA Calculate the protection level PL1 caused by positioning deviation; S4: Based on the missed detection rate P MD Calculate the protection level PL2 caused by positioning noise; S5: Determine the final protection level PL of the current epoch, and issue an alarm to the user when the protection level PL exceeds the alarm limit AL of the specific application field; when the protection level PL is less than or equal to the alarm limit AL of the specific application field, calculate the protection level of the next epoch.

2. The method for calculating the protection level of the GNSS positioning state domain according to claim 1, characterized in that: In the step S1, when m>k, P dd,k =P′ xx,k -P xx,k , where P′ xx,k is the prior state estimate covariance matrix of the current epoch; when m≤k, the recursive state estimate Recursive state estimation covariance matrix Among them, the initial recursive state estimate Its covariance matrix f(·) is the state transition function, Indicates that the state vector x obeys the mean Covariance matrix cov[x] = P′ xx Gaussian distribution, Q is the system noise covariance matrix, the difference between the recursive state estimate and the posterior state estimate The covariance matrix of 3. The method for calculating the protection level of the GNSS positioning state domain according to claim 1, characterized in that: The protection level PL1 caused by the positioning deviation in step S3 is divided into a horizontal protection level HPL1 and a vertical protection level VPL1. When only the positioning deviation of the current epoch is considered, there is in, and and They are Λ dd,k The eigenvalue components in the east, north and sky directions; when considering the positioning deviation within the recursive period, we have in, and and They are Eigenvalue components in the east, north, and sky directions.

4. The method for calculating the protection level of the GNSS positioning state domain according to claim 3, characterized in that: In step S4, the protection level PL2 caused by the positioning noise is divided into a horizontal protection level HPL2 and a vertical protection level VPL2. in, and and P xx,k In the three directions of east, north and sky, P MD is the missed detection rate based on integrity risk allocation.

5. The method for calculating the protection level of the GNSS positioning state domain according to claim 4, characterized in that: In step S5, the final protection level PL of the current epoch is divided into a horizontal protection level HPL and a vertical protection level VPL. Considering only the positioning deviation of the current epoch, HPL = HPL1 + HPL2, VPL = VPL1 + VPL2; considering the positioning deviation within the recursive period, The alarm limit AL is divided into the horizontal alarm limit HAL and the vertical alarm limit VAL; when the horizontal protection level HPL is greater than the horizontal alarm limit HAL or the vertical protection level VPL is greater than the vertical alarm limit VAL, an alarm is issued to the user.

6. The method for calculating the protection level of the GNSS positioning state domain according to claim 1, characterized in that: In step S2, the process of fault detection and elimination is as follows: (1) Rearrange the GNSS observations of the current faulty epoch from small to large according to the normalized new information, and take the first n satellite observations to calculate the posterior state estimate and the covariance matrix P 0,xx,k , and based on the posterior state estimate and the covariance matrix P 0,xx,k Construct a fault detection statistic and compare it with the detection threshold. When the fault detection statistic is greater than the detection threshold, it means that the quality of the observation data of the current epoch is poor and the fault cannot be eliminated. An alarm is issued to the user and the fault elimination of the current epoch is terminated. When the fault detection statistic is less than or equal to the detection threshold, it means that the fault can be eliminated and proceed to the next step. (2) Delete the GNSS observation corresponding to the maximum innovation and use the remaining satellite observations to calculate the posterior state estimate and the covariance matrix P test,xx,k , the posterior state estimate and the covariance matrix P test,xx,k Construct a fault detection statistic and compare it with the detection threshold. If the fault detection statistic is greater than the detection threshold, it means that a fault exists and the fault has not been completely eliminated. Delete the GNSS observation value corresponding to the maximum innovation in the remaining observation values. Repeat this step until the fault is completely eliminated, and then perform fault detection and elimination for the next epoch. When the fault detection statistic is less than or equal to the detection threshold, it means that the fault has been eliminated and the next epoch of fault detection and elimination can be directly carried out.

Images