# Research on GNSS Receiver Autonomous Integrity Monitoring Method Based on M-Estimation

Authors: Pengfei Zhang; Pengyun Chen; Dan Song; Guochao Fan
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Receiver Autonomous Integrity Monitoring (RAIM) method is an effective means to provide integrity monitoring for users in time. In order to solve the misjudgment caused by the interference of gross error to the least squares algorithm, this paper proposes a RAIM method based on M-estimation for multiconstellation GNSS. Based on five programs, BDS, GPS/BDS, and GPS/BDS/GLONASS at the current stage, the future Beidou Global Navigation Satellite System, and the future GPS/BDS/GLONASS/Galileo system, the new RAIM method is compared with the traditional least squares method by simulation. The simulation results show that, with the increase of constellations, RAIM availability, fault detection probability, and fault identification probability will be improved. Under the same simulation conditions, the fault detection and identification probabilities based on M-estimation are higher than those based on least squares estimation, and M-estimation is more sensitive to minor deviation than least squares estimation.


---

## Body

## 1. Introduction

Receiver Autonomous Integrity Monitoring (RAIM) is an effective method of integrity monitoring [1]. With the rapid development of Chinese Beidou navigation system and European Galileo system and recovering of Russia’s GLONASS constellation, interoperability among the four global navigation satellite systems (GNSS) has become an inevitable trend. In addition, research on multiconstellation RAIM algorithm has also been promoted [2–5]. RAIM can respond quickly and completely to the satellite fault and aerial abnormality, without any external intervention; it can also provide alarm information for users timely and effectively. Moreover, the user-level integrity monitoring can be realized simply and the input cost is low [6]. Therefore, research on RAIM algorithm is very necessary, and reliable RAIM algorithm can guarantee the integrity of navigation and positioning results. In 1987, Kalafus first introduced the concept of RAIM, and the consistency check for the current time redundancy variable was carried out, and at least five visible satellites are required for fault detection; at least six visible satellites are required for fault identification and exclusion. Afterward, many domestic and foreign experts proposed different snapshot algorithm and filter algorithm [6–9], but most of the above studies are based on GPS, GLONASS, and Galileo systems. The RAIM method research on the Beidou satellite navigation system and its combination with other GNSS is still in its infancy and development stage. This paper mainly applies the snapshot algorithm to multiconstellation fault detection and identification and proposes a RAIM method based on M-estimation for multiconstellation. Five programs, BDS, GPS/BDS, GPS/BDS/GLONASS at the current stage, the future Beidou Global Navigation Satellite System, and the future GPS/BDS/GLONASS/Galileo system, are compared with the traditional RAIM method based on least squares in simulation analysis.


## 2. RAIM Method Based on Least Squares for Multiconstellation

Suppose that the receiver pseudorange observation equation is(1)Y=HX+ε,where ε is a n×1 vector, representing the pseudorange measurement error for each satellite.And the least square positioning solution is(2)X^=HTH-1HTY.Make(3)A=HTH-1HTand then(4)X^=AY,Y^=HX^.Pseudorange residual vector:(5)w=Y-Y^=Y-HHTH-1HTY=In-HHTH-1HTY=In-HHTH-1HTHX+ε=In-HHTH-1HTε.Make(6)S=In-HHTH-1HT.S is called residual sensitivity matrix; then the residual sum of squares can be expressed as(7)SSE=wTw=εTS2ε=εTSε.S is symmetric, idempotent, and its sum of squares of each row and column is equal to the corresponding diagonal element, and the sum of each row or column is equal to 0.MakeG=HHTH-1HT, and then(8)rankG=m+3,rankSSE=rankS=rankI-G=n-m-3,where n is the number of visible satellites and m is the number of navigation satellites involved in the calculation.The error in the a posteriori unit weight of the pseudorange residual vector is(9)σ^=SSEn-m-3.Therefore, the unit weight error of the pseudorange residual vectorσ^ is calculated by the sum of squares of the pseudorange residuals. Under the normal circumstance of the system, the residual of the pseudorange is small, and the a posteriori unit weight error σ^ is also small. When the deviation of the measurement pseudorange is large, σ^ will become larger and it needs to be detected. Assuming that there is no fault, each component of the distance residual vector w is independent of the normal distribution random error with the mean zero and variance σ02. Because the residual sensitivity matrix S is a real symmetric matrix whose rank is equal to n-m-3, according to the statistical distribution theory, SSE/σ02 obeys chi-square distribution with freedom degree n-m-3; if there is a fault and the mean value of the distance residual vector w is not zero, then SSE/σ02 obeys the noncentral chi-square distribution with freedom degree n-m-3 [9]. Therefore, σ^ can be used as a test statistic. Let the test statistic be(10)TS=σ^=SSEn-m-3.The detection thresholdTD can be calculated by the allowed maximum false alarm probability Pfa. The false alarm is an indication that the user is notified of a positioning fault when no positioning fault occurs. The detection threshold is obtained by detecting the probability density function from the detection threshold to the infinite integral. This problem belongs to hypothesis test problem in mathematical statistics. Suppose the following:H0: no fault occurred.H1: a fault occurred.Then the false alarm probability(11)Pfa=PSSEσ02≥TD2∣H0.According to the above equation, the detection threshold can be obtained. IfTD>TS, it indicates that the fault is detected, and the alarm will be sent to the user.Fault detection is based on the test of pseudorange residuals sum of squares, and fault identification is based on the test of pseudorange residual element. And the basic idea is based on Baarda’s data snooping method [10], constructing a statistic based on least squares residual vector. The statistic obeys a certain distribution, and given a significant level, the residuals can be judged by the statistical test whether a gross error exists. From the relationship between residual and observation error, we can make the fault identification test statistic(12)TSi=wiσ0QVii,where QVii represents the i row and i column element of the pseudorange residual vector cofactor matrix QV.(13)QV=W-1-HHTWH-1HT.From (12), we can see that the statistical distribution of TSi is consistent with wi.When there is no fault,TSi~N(0,1); when there is a fault, TSi~N(δi,1).Whereδi is the statistic offset parameter, if satellite i has a pseudorange error bi, then(14)δi=QViiWiibiσ0.n visible satellites can obtain n test statistics. Given the total false alarm probability pfa, then the false alarm probability of each test statistic is pfa/n, and the detection threshold TDi of fault identification can be calculated from the false alarm probability of each test statistic.(15)PTSi>TDi∣H0=2∫Td∞hxdx=pfan,where(16)hx=12πe-x2/2x∈R.From the above formula we can calculate the identification thresholdTDi corresponding to each test statistic TSi. Compared to the test statistic with the identification threshold, if TSi>TDi, then the ith satellite is faulty and should be excluded.


## 3. RAIM Method Based on M-Estimation for Multiconstellation

M-estimation used in this paper is an iterative weighted least squares estimator. Different weights are applied to different points according to the pseudorange residual vector; that is, the points with small residuals are given a larger weight, while those with larger residuals are given a smaller weight. And weighted least squares estimation is then established, repeatedly iterating to improve the weight coefficient.Different from the least squares making pseudorange residual sum as the extreme function, the extreme function of M-estimation is(17)PX=∑i=1npiρwi=∑i=1npiρhiX-Yi,where ρ(x)=-ln⁡f(x), f(x) being the probability density function, take the derivative of w and let it equal zero, make ψ(wi)=∂ρ/∂wi, and then (18)∑i=1npiψwihi=0.Makeψ(wi)/wi=Wi (the weighting factor), and introduce equivalent weight element p-ii=piWi; then the above formula is rewritten as(19)HTP¯w=0,where P¯ is the equivalent weight matrix (diagonal matrix), and the element is p-ii; bring it into the error equation; then(20)HTP-HX-HTP-Y=0.Thus the M estimated value of the robustness of the parameter vector is(21)X=HTP-H-1HTP-Y,where P--1=p-1120⋯00p-222⋯0⋯⋯⋯⋯00⋯p-nn2.There are many methods to construct the equivalent weight matrix, but the robust estimates are much the same, and a “normalized” residual indexui is used in each method, which is defined as(22)ui=wis=0.6745×wimedwi-medwi,where med is median, and s is residual scale. The Huber method is chosen; that is,(23)p-ii=1ui≤chchuiui>ch,where ch generally takes 1.345 [11].The following fault detection and identification methods are similar to the least squares RAIM method, and construct test statistics:(24)TS=σ^=wTP-wn-m-3.No fault assumptionH0: σ^2/σ02~χ2(n-m-3).A fault assumptionH1: σ^2/σ02~χ2(n-m-3,λ).λ is the noncentralization parameter, and σ02 is the a priori variance of the pseudorange residual; calculate the corresponding detection threshold from the corresponding false alarm probability. Then the fault detection and fault identification are carried out.


## 4. Simulation Analysis

This paper designs the following five programs, using the self-compiled software for simulation analysis.Program 1. It is the currently operating Beidou regional navigation system BD2 (5GEO + 5IGSO + 4MEO), with a total of 14 satellites.Program 2. It is the Beidou system, BDS, for the future global navigation satellite system (5GEO + 3IGSO + 27MEO), with a total of 35 satellites.Program 3. It is the currently operating GPS (32 satellites) + BD2 (5GEO + 5IGSO + 4MEO), with a total of 46 satellites.Program 4. It is the currently operating GPS (32 satellites) + BD2 (5GEO + 5IGSO + 4MEO) + GLONASS (24 satellites), with a total of 70 satellites.Program 5. It is a complete combination of four major satellite navigation systems GPS (32 satellites) + BDS (5GEO + 3IGSO + 27MEO) + GLONASS (24 satellites) + Galileo system (27 satellites) in the future, with a total of 118 satellites.The Beidou regional constellation BD2, GPS constellation, and GLONASS constellation all adopt the broadcast ephemeris of 2015-01-12. The Beidou system, BDS, for the future global navigation satellite system will be simulated with 35 satellites (5GEO, 3IGSO, and 27MEO); the Galileo system is simulated with 27 satellites; the specific parameters are shown in Table1.Table 1
Constellation parameters for BDS and Galileo system.
Navigation system BDS Galileo system Orbit type GEO IGSO MEO MEO  Orbital plane 1 3 3 3 Satellite number 5 3 27 27 Semimajor axis/km 42164 42164 27906 29978 Orbit eccentricity 0 0 0 0 Orbit inclination 0° 55° 55° 56° Ascending node Fixed in the longitude of 58.75°, 80°, 110.5°, 140°, 160° The intersection longitude is east longitude 118° 70°, 190°, 310° 0°, 120°, 240° Mean anomaly The initial time near point angle of first satellite in each orbit is, respectively, 0°, 15°, 30°, followed by an increase of 40° Argument of perigee 0° 0° 0° 0°First, select Program 2 (GPS + BD2) double constellations, and two kinds of RAIM methods are simulated and compared; the basic simulation conditions are shown in Table2.Table 2
Basic simulation conditions.
Program Condition Ephemeris reference time 2015.1.12 00:00:00 Simulation constellation GPS + BD2 Simulation area (40°N, 116°E) Standard deviation of pseudorange noise 5 m Obstacle angle 10° Simulation step 5 s Simulation cycle 24 h False alarm probability 1/3000000During the specified simulation time interval, 60 m deviation is injected in a satellite. Test statistics and detection threshold statistics of all sampling points of two methods are shown in Figure1.Figure 1
Test statistic and threshold comparison of the two methods separately.It can be seen from Figure1 that both of two methods can be used for multiconstellation fault detection, but under the same circumstances, fault detection probability based on M-estimation is significantly higher than that based on least squares. Since the M-estimation is able to amplify the pseudorange residuals of the failed satellite in the test statistic, it is more sensitive than the least squares to the minor deviations; that is, the correct warning probability of the least squares to the minor deviations is lower than that of M-estimation.In addition, the M-estimation robustness performance is better; when the fault is not eliminated, an iteratively weighted method based on the deviation is selected for M-estimation, and the impact of the deviation on positioning performance can be reduced eventually, as shown in Figure2.Figure 2
Positioning errors of the two methods on the condition that the fault is not excluded.
(a)
Least squares (b)
M-estimationAs shown in Figure2, 10 m deviation is added in a specified period. If the fault is not excluded, it directly participates in the positioning calculation. And the positioning results of two methods are shown. Obviously, M-estimation has better robustness performance.Therefore, this paper selects the multiconstellation RAIM method based on M-estimation to simulate and analyze five programs separately. Add the deviation from 5 m to 120 m in the faulty satellite, and the step is 5 m. Simulation with the Monte Carlo method, fault detection probability, and fault identification probability corresponding to five programs are shown in Figures3 and 4, respectively.Figure 3
Fault detection probability of M-estimation.Figure 4
Fault identification probability of M-estimation.From Figures3 and 4, it can be seen that, with the increase of fault deviation, the fault detection probability and fault identification probability of each program are improved. When the deviation is greater than a certain value, the fault detection probability and fault identification probability can reach 100%. In addition, multiconstellation GNSS fault detection capability and fault identification capability are better than those of single-constellation. With the increase of constellation numbers, fault detection capability and fault identification capability are gradually increased. Multiconstellation GNSS is more sensitive to smaller fault deviations. This is because the multiconstellation GNSS increases the redundant information of RAIM by increasing the number of visible satellites, thus improving the fault detection probability and fault identification probability. At the same time, under the same fault deviation, the fault identification probability of the same program is lower than the fault detection probability, which is due to the fact that the fault identification is more strict than the fault detection. At current stage, for BD2, in the service area, when the fault deviation is greater than 100 m, the fault detection probability can reach 100%; when the fault deviation is greater than 110 m, the fault identification probability can reach 100%. For BD2/GPS, in the global area, when the fault deviation is greater than 75 m, the fault detection probability can reach 100%; when the fault deviation is greater than 85 m, the fault identification probability is 100%. For BD2/GPS/GLONASS, in the global area, when the fault deviation is greater than 55 m, the fault detection probability can reach 100%; when the fault deviation is greater than 65 m, the fault identification probability can reach 100%. For BDS, in the future, the fault detection probability can reach 100% when the fault deviation is greater than 85 m, and the fault identification probability can reach 100% when the fault deviation is greater than 95 m. For BDS/GPS/GLONASS/Galileo system, in the future, when the fault deviation is more than 50 m, the fault identification probability can reach 100%, and the fault detection probability can reach 100%.


## 5. Conclusion

This paper firstly introduces the RAIM method of multiconstellation based on traditional least squares and deduces the test statistic and threshold calculation process of fault detection and fault identification. Then, a multiconstellation RAIM method based on M-estimation is proposed in this paper; at the same time, the fault detection and identification process are deduced. Finally, for the five programs, including BD2, GPS/BD2, and GPS/BD2/GLONASS at the current stage and Beidou Global Navigation Satellite System and GPS/BDS/GLONASS/Galileo System in the future, the RAIM method based on M-estimation is compared with the traditional RAIM method based on least squares by simulation. The simulation results show that the availability of RAIM method increases with the number of constellations, and the fault detection probability and fault identification probability also increase. Under the same condition, the fault detection probability and fault identification probability based on M-estimation method are higher than those based on least squares method, and the M-estimation is more sensitive to the minor deviation than the least squares estimation.



---

*Source: 2563202-2018-03-28.xml*