Analysis of BDS GEO satellite multipath effect for GNSS integrity monitoring in civil aviation

With the great development of Global Navigation Satellite System (GNSS), multi-GNSS constellations (GPS, BDS, GLONASS, and Galileo) are able to provide users with more accurate positioning result. For civil aviation, to guarantee user’s safety, multi-constellation GNSS needs to meet the integrity requirement. Using conservative error models, Multiple Hypothesis Solution Separation (MHSS) Advanced Receiver Autonomous Integrity Monitoring (ARAIM) is proposed to evaluate GNSS integrity. Current multipath error model in ARAIM algorithm is based on data of GPS Medium Earth Orbit (MEO) satellites. But BDS is a hybrid constellation. For BDS II, there has 5 Geosynchronous Earth Orbit (GEO) satellites. Previous studies have shown that the multipath effect of GEO satellites has statistical characteristics different from MEO satellites. Meanwhile, the multipath magnitude of GEO satellites is larger than that of MEO satellites. This paper mainly focuses on validating whether the multipath error model in ARAIM algorithm is conservative enough for GEO satellites. In this paper, Code Minus Carrier (CMC) residuals are calculated for BDS GEO dual-frequency signals. Then the Standard Deviation (STD) of CMC residuals can be conservatively estimated by bounded Cumulative Distribution Function (CDF). After eliminating interference from receiver noise, STD of GEO multipath can be obtained. Comparing the STD of GEO multipath effect with ARAIM multipath error model, a conclusion could be drawn that current multipath error model in ARAIM algorithm is no longer able to conservatively estimate the statistical characteristics of GEO satellites.


Introduction
With the ability of providing positioning result worldwide, satellite-based navigation has been widely used in civil aviation. For civil aircraft, ensuring user's safety is of crucial significance. To guarantee the safety during the flight, International Civil Aviation Organization (ICAO) put forward four criteria to evaluate the performance of satellite-based navigation system: integrity, continuity, accuracy, availability [1]. Among them, two of the most challenging requirements for civil aviation are integrity and continuity [2]. Integrity is a measure of trust, which is used to determine whether the positioning results provided by navigation system is correct [3]. Continuity measures the ability of navigation system to operate without unplanned interruptions [1,4].
Before 2000s, Receiver Autonomous Integrity Monitoring (RAIM) had been developed to measure the integrity for civil aviation [5][6][7]. By checking the consistency of positioning residuals, RAIM becomes a supplemental navigation means for en-route through non-precision approach (NPA) phases of flight [8]. However, RAIM has a lot of drawbacks. One of the major drawbacks is occasional lack of availability [9]. With the navigation satellites increasing, using dual-frequency signals and multi-constellation positioning, MHSS ARAIM theory is proposed. MHSS ARAIM has better ability to provide integrity service. Many researchers have made great contributions to ARAIM user algorithm. Juan Blanch and Boris Pervan proposed Fault Detection (FD) scheme and Fault Exclusion (FE) scheme for ARAIM user algorithm [10,11]. To determine the error magnitude of GPS satellite, Todd Walter analyzed GPS clock error as well as ephemeris error in detail and established the corresponding error model for these two errors [12]. In [13], Todd Walter specified two critical parameters for ARAIM user algorithm: the probability of satellite fault and the probability of constellation fault. As for the multipath error model, Tim Murphy used GPS data from Boeing 777-300ER and 737-NG airplanes to build the STD-elevation error model [14].
BDS constellation is composed of MEO, Inclined Geosynchronous Orbit (IGSO) and GEO satellites. In Wide Area Augmentation System (WAAS), the ranging error of GEO satellites stands in a bias over a long time [15,16] analyzed the fading characteristics of satellites in different obits to explain this phenomenon. Gao et al. [16] shows that the orbit type of navigation satellites influences the multipath fading characteristic. As for GEO satellite in BDS II, [17] demonstrates that the multipath magnitude of BDS GEO satellites is larger than those of IGSO and MEO satellites. Moreover, the multipath effect is difficult to mitigate by increasing the observation period or averaging over multiple epochs [18]. Therefore, multipath error is one of error sources that could affect GNSS integrity. To determine the airframe multipath error model for GPS constellation, the initial model is proposed by [19] and [20] using the airframe multipath data. However, this model is shown to be insufficient [21]. Thus, the revised model is proposed by [22], which is also the current multipath error model for ARAIM algorithm. Current multipath error model in ARAIM user algorithm is based on GPS MEO satellite data. But this error model does not take BDS GEO satellites into consideration. When implementing GNSS integrity monitoring involving BDS constellation, the multipath error model of GEO needs to be determined.
To determine whether current multipath error model in ARAIM algorithm is conservative for BDS GEO satellites, this paper uses 24-h data with 1 Hz sample rate to validate the error model. The data is collected from chock ring antenna in Shanghai Jiao Tong University, which could mitigate the multipath effect. Meanwhile, comparing to the reflection of aircraft, since there are no obstacles around the antenna, the multipath magnitude of antenna is lower than that of aircraft. Code Minus Carrier (CMC) method is firstly used to calculate the margin of single frequency multipath residual. In this process, two methods are used to avoid cycle split effect. Due to ARAIM algorithm needs to conservatively guarantee user's safety [23], the worst-case scenario of multipath effect needs to be taken into account. Bounded CDF is used to determine the worst-case STD of multipath residual for BDS GEO satellites. However, the worst-case STD of multipath residual is contaminated by receiver noise. In the following step, the STD of receiver noise is analyzed using GNSS zero baseline data. After eliminating the interference of receiver noise, the multipath statistic characteristic of BDS GEO satellites can be obtained.
The paper is structured as follows. In the first part, the basic outline of MHSS based ARAIM user algorithm is given. The second part shows the method to calculate the multipath residual. This part also introduces how to determine the conservative STD for multipath effect of BDS GEO satellites. After eliminating the interference of receiver noise, this part calculates statistic characteristic of multipath error. The last part is conclusion. where H is a n × (3 + m) matrix. n is the number of pseudorange measurements and m is the number of constellations.

ARAIM user algorithm
x is a (3 + m) × 1 vector, which includes the user's position and clock bias between different GNSS constellations. f is a n × 1 vector, which indicates the fault magnitude in navigation system. is the noise vector, which includes the noise from satellite clock, ephemeris, troposphere, multipath and user's receiver. In order to calculate x , Weighted Least Square (WLS) method is used to solve Eq. (1): In (2), W is the weighted matrix, which is a diagonal matrix and can be derived from the error modes of ARAIM user algorithm. For SS ARAIM, we build N fault subsets to detect the possible faults which can severely affect GNSS accuracy. For detection function, positioning result of the i -th fault subset can be defined as where W i is calculated by weight matrix W after removing the measurements corresponding to the i th fault subset. Here gives an example of calculating W i . Assuming the algorithm only monitors single satellite (the i th fault subset corresponds to the i th measurement) fault, W i can be expressed as To detect satellite fault in navigation system, detection statistic needs to be calculated. We define x 0 as all in view positioning result, then the detection statistic can be calculated according to Based on the probability of false alert, the threshold T i for different detection statistics can be determined. If Δ i > T i , it means the fault has occurred in navigation system. To meet the continuity requirement of navigation system, fault exclusion function needs to be operated.
MHSS ARAIM builds N 1 sub-subsets to determine whether the exclusion function is implemented correctly. The positioning result of sub-subset can be expressed as Here gives an example of calculating W i,j . Assuming only single satellite fault is monitored by ARAIM algorithm, then both subset and sub-subset contain one satellite. If algorithm has detected GNSS fault. Fault exclusion function tries to exclude the faulty satellite to improve navigation performance. Firstly, FE algorithm chooses the satellite in the i th subset as the exclusion option. For a specific sub-subset j , W i,j can be denoted as To measure the correctness of exclusion option, exclusion statistic needs to be calculated. Based on subset and sub-subset positioning result, the exclusion statistic can be calculated according to Using the probability of fault detection but not exclusion, the threshold T ex,ij for different exclusion statistics can be obtained. If exclusion option i is determined and ∀{j ∈ [1, N]|j ≠ i} meet Δ i,j < T ex,ij , it represents the i th subset is the one that need to be excluded. Otherwise, the algorithm will change the exclusion option.

Error model of ARAIM
The error sources of GNSS can be divided into three types: (a) the error comes from satellite, which includes satellite clock and ephemeris error. (b) the error comes from signal propagation, which includes tropospheric and ionospheric error. (c) the error comes from user's environment, which includes multipath and receiver noise error.
To determine the statistic characteristics of these errors, [23] proposed corresponding error model to estimate the STD of them. For errors in (a), their STD can be measured by URE .
URE can be obtained from broadcast ephemeris. The STD of tropospheric error tropo in (b) can be measured by mode listed in (9) [24]. For MHSS ARAIM algorithm, it does not consider the ionospheric error, because the ionospheric error can be canceled using dual frequency signals. For airborne GPS users, error models to determine the STD of error in (c) is well developed. The error models can be measured by (10), (11) and (12). In (10), f L1 and f L5 represent the signal frequency of GPS L1 and L5 signal. user is the STD for airborne user error. MP and nosie represent the STD of multipath and receiver noise, respectively: For every GNSS satellite, the Eq. (9), (10) and URE can derive the weight matrix that used for WLS method to solve (1). The weight matrix can be expressed as where 2 i = 2 URE,i + 2 tropo,i + 2 user,i , the subscript i represents the i-th diagonal element of W.

Code minus carrier method
For GNSS receiver, multipath delay results in tracking error of code correlator. The tracking error causes GNSS receiver to generate faulty code and carrier phase measurement. Multipath effect can severely affect the magnitude of code measurement, because the wavelength of carrier phase is much smaller than chip interval [25]. According to this property, multipath residuals for different GNSS signal frequencies can be calculated through the combinations of code measurement and carrier phase measurement in (14) and (15) [18,18,18]: In (14) and (15), MP f 1 and MP f 2 are multipath residuals of frequency f 1 and f 2 , respectively.N f 1 and N f 2 are integer ambiguities which can be seemed as constant bias when there is no cycle slip occurs [26]. 1 and 2 are the wavelength of GNSS signal at f 1 and f 2 frequency, respectively. is GNSS receiver noise. Actually, the residuals in (14) and (15) contain the integer ambiguities which needs to be eliminated.

Cycle slip detection and repairation
Before eliminating the constant term in (14) and (15), the effect of cycle slip needs to be avoided. Based on [28], we combine Total Electron Contents Rate (TECR) and Melbourne-Wübbena wide-lane (MWWL) method to detect and repair the cycle slip.
Assume that the ionospheric delay changes very slowly with time. Test statistic can be derived from ionospheric residual from epoch k and epoch k − 1 . (16) is the equation to calculate TECR test statistic. In (16), 1,k is carrier phase observation of f 1 frequency signal at epoch k . 2,k is the carrier phase observation of f 2 frequency signal at epoch k − 1 . According to the law of error propagation, the threshold T TECR,k of Δ TECR,k is set to 0.07 [29]. (17) is the equation that can be used to calculate the cycle split magnitude ΔN 1,k and ΔN 2,k at frequency f 1 and f 2 at epoch k: TECR has its own drawback. When 1 ΔN 1,k = 2 ΔN 2,k , TECR cannot detect the cycle split. To solve this problem and repair the cycle split. MWWL test statistic Δ MWWL,k is introduced [29]: In (18), PR 1,k is pseudorange measurement of frequency f 1 at epoch k . PR 2,k is pseudorange measurement of frequency f 2 at epochk . The threshold T MWWL,k of test statistic Δ MWWL,k is set to 1 based on [29]. (20) can be used to calculate the cycle split magnitude ΔN 1,k and ΔN 2,k at frequency f 1 and f 2 at epochk.
If |Δ MWWL,k | > T MWWL,k or |Δ TECR,k | > T TECR,k , the cycle split is detected. Then, combining (17) and (20), the magnitude of cycle split can be obtained. To validate the effectiveness of this method, we first select 200 s GNSS observation data from BDS GEO 01 satellite which is not affected by cycle slip. Figure 1 shows the magnitude of TECR and MWWL test statistic.
Assuming that B1 frequency carrier phase observations are affected by cycle slip at 100 s. The magnitude of cycle slip is set to 100 cycle. Figure 2 shows the result of TECR and MWWL test statistic.
In Fig. 2, both TECR and MWWL method can detect the cycle slip. Based on (17) and (20), the magnitude of cycle slip in B1 and B2 frequencies can be expressed as follows: where [] represents the rounding process. Because the magnitude of cycle slip can only be integer, the repair results of cycle slip need to be rounded. In this case, the repair results are 99 cycle and 1 cycle for B1 and B2 frequency, respectively. The results show that the method can effectively repair the cycle slip.

Multipath residual calculation
After removing the cycle split effect from carrier phase measurement, the constant term in (14) and (15)

Receiver noise estimation using zero baseline test
In (23), the multipath residual still contains the receiver noise.
To remove the interference of receiver noise, GNSS zero baseline data is used to estimate the receiver noise. Zero baseline test is a method to assess the intrinsic noise characteristics of GNSS receiver [30]. The test requires two receivers of the same type. Two receivers need to connect to the same antenna. The experiment environment is shown in Fig. 3. The pseudorange of receiver a and b can be expressed in (24) and (25):  SD is receiver noise of SD result. In SD result, external errors such as atmospheric errors, satellite-dependent errors and multipath can be canceled. Receiver noise and clock bias are the only two errors in SD result.
Since SD result is affected by clock bias, DD residual ∇ΔPR ba i which only contains receiver noise can be represented by (27). The subscript j is the reference satellite with high Carrier to Noise Ratio (CNR). DD is receiver noise of DD result.
Time series of receiver noise for different satellite can be obtained from DD residual. However, the DD residual of the reference satellites cannot be obtained. Therefore, based on [31], this paper uses DD result to recover SD sequence which is not affected by receiver clock bias. Assuming that the first satellite is reference satellite, DD sequences can be expressed as where the subscript 1i represents the DD residual between the first and the i th satellite. The subscript of i represents the i th satellite. Since GNSS intrinsic receiver noise is Gaussian white noise, the sum of SD sequence obeys Gaussian DD,1n = SD,1 − SD,n distribution with zero bias. Using this character and solving (29), SD sequence can be obtained: After bounding SD sequence and multipath residual, the statistic characteristic of BDS II GEO multipath can be calculated. By comparing with STD calculated from ARAIM model, this paper examines whether ARAIM multipath model could be used for BDS GEO satellite.

Multipath residual for BDS GEO satellite
24-h data is used to calculate the multipath residual of BDS GEO satellites. The STD of multipath residual can be conservatively estimated by bounding process. Due to the low elevation angle, the signal from BDS 05 satellite is highly susceptible to the interference. So, in this paper, we only analyze four BDS II GEO satellites. Table 1 shows the bounded STD value MR of BDS 01 to 04 satellite at different frequencies (B1 and B2 frequency). Figure 4 shows the multipath residual and bounding results. Take BDS 01 satellite as an example, Figure 4a, b show the variation of multipath residual and elevation over time. The orange line and blue line represent the elevation angle and multipath residual. Figure 4c,d shows the bounding results, the blue dot represents the folded CDF of multipath residual data. The green line is the bounded CDF, which can conservatively estimate the STD of multipath residual data.

SD sequence for BDS GEO satellite
Even the STD of multipath residual can be obtained, the result still affected by receiver noise. To eliminate the interference of receiver noise, receiver noise evaluation is implemented. For receiver noise evaluation, this paper uses 8-h data to obtain SD sequence for different GEO satellites. This paper only analyzes the receiver noise for BDS 01 to 04. Table 2 shows the bounded STD value SD of SD sequence at different frequencies (B1 and B2 frequency). Figure 5 shows the time series of GNSS receiver noise. At the same time, Fig. 5 also shows the bounding result of receiver noise for different BDS GEO satellite at difference frequencies. Take GEO 01 satellite as an example, Fig. 5a, b shows the SD sequence result of receiver noise. Figure 5c, d shows the bounding results which can conservatively estimate the STD of SD sequence.

Multipath magnitude evaluation for BDS GEO satellite
Based on previous analysis in this paper, multipath residual is affected by receiver noise. Therefore, the STD of multipath residual MR can be expressed as where noise is receiver noise. Based on (30), the STD of multipath can be obtained by following: The relationship between the STD of SD result and noise can be described as where SD can be obtained from Table 2. Using (30), (31) and (32), MP for BDS GEO 01 to 04 is shown in Table 3. Figure 6 shows the difference of STD between ARAIM model and BDS GEO satellites. The circles in Fig. 6 represents different frequencies of BDS. The dots with different color in the circle represent different GEO satellites.
As shown in Fig. 6, whether at B1 or B2 frequency, the STD values of BDS GEO satellites are far larger than ARAIM error model.

Conclusions
This work investigated whether the multipath error model in ARAIM algorithm is suitable for BDS II GEO satellites. This paper firstly calculates the multipath residual using CMC method. In this step, to avoid the interference of cycle slip, TECR and MWWL method is used to detect and repair the cycle slip. Since the results of multipath residual are contaminated by receiver noise, GNSS zero baseline data is used to estimate the receiver noise. After eliminating the interference of receiver noise, we conservatively estimate the STD of multipath residual for BDS II GEO satellites. We find that the error model in ARAIM user algorithm are not suitable for BDS II GEO satellites. In the future, for GNSS integrity monitoring system using BDS II GEO satellites, the multipath error model of these satellites should be refined.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.