Carrier Phase-Based Ionospheric Gradient Monitor Under the Mixed Gaussian Distribution

Abstract

Anomalous ionospheric gradient is a critical risk to be monitored by ground-based augmentation systems (GBASs) utilized for safety-of-life navigation applications. A dual-frequency carrier phase-based ionospheric gradient monitoring method is proposed under the mixed Gaussian distribution. The minimum detection error of the proposed method can be greatly reduced by allowing acceptable ambiguity resolution failure modes, given the required averaging length. The real BeiDou navigation satellite system data were utilized to test the proposed method. The experimental results showed that the minimum detection error (MDE) of the proposed dual-frequency ionospheric gradient monitoring method can be reduced by at least 30% in comparison with the maximum acceptable anomalous ionospheric gradient of category III GBAS. This study demonstrated that the proposed method can be used to protect against the ionospheric gradient for a ground-based augmentation system.

1. Introduction

Global navigation satellite systems (GNSS) play an important role in civil aviation applications [1,2]. The ground-based augmentation system (GBAS), based on the differential architecture of GNSS, is a new generation of civil aviation navigation systems to replace the traditional instrument landing system [3,4,5]. The local ground facility of GBAS consists of multiple spatially separated GNSS antennas and receivers. The primary reason for the multiple reference receiver configuration at the local ground facility is to detect and isolate the faulty receivers. The multiple reference receiver configuration provides a hardware foundation for space atmosphere gradient monitoring based on double-difference carrier phase observations [6,7]. Anomalous ionospheric gradient is one of the most challenging risk sources during the precision approach [8,9,10,11].

Since the variation of ionosphere is highly influenced by solar activity intensity and magnetic latitude, it is difficult to eliminate the positioning errors caused by ionospheric delay through a precise ionosphere model [12,13,14]. Since the baseline between reference stations is relatively short, when the ionosphere is quiet, the differential positioning errors caused by ionospheric delay is usually less than 10cm, which is acceptable for the aircraft on category III precision approach supported by GBAS [15]. However, during an extreme ionospheric storm, the anomalous fluctuation of electron concentration may produce a spatial decorrelation ionospheric gradient that can hardly be eliminated by differential processing. It can be anticipated that the undetected ionospheric gradient will pose a threat to positioning. According to Luo et al., Pullen et al., Kim et al., Lee et al., and Wang et al. [9,16,17,18,19], anomalous ionospheric gradients have been detected in continuously operating reference stations (CORS) in at least four regions, i.e., the United States, South Korea, Brazil, and China. Specifically, the minimum ionospheric gradient observed in South Korea also reaches 160mm/km, and when the geometric dilution of precision is poor, it can cause up to 10 m of vertical positioning error in 5 km baseline. Obviously, such a large error will endanger the precise approach and landing guidance services supported by GBAS. Therefore, the anomalous ionospheric gradient must be monitored for safety-of-life navigation applications to avoid unacceptable positioning errors or even availability loss of GNSS-based positioning.

Compared with code observations, carrier phase observations can detect ionospheric gradient with higher sensitivity when the ambiguity is correctly resolved. However, in the presence of the ionospheric storms and scintillations, it is a challenging task to correctly resolve the ambiguity to achieve high sensitivity test statistics. At present, according to the different processing methods of ambiguity in the construction of test statistics, there are mainly two types of differential carrier phase-based ionospheric gradient monitoring (IGM) methods.

The first type is defined as the instantaneous IGM method without the accumulation of historical observations [6,7,20]. The instantaneous IGM method is implemented under the Gaussian distribution assumption to control false alarm and missed detection errors, simultaneously. Since the instantaneous IGM method pays attention to the detection initialization time, the performance of this method is limited by the coupling between the ambiguity and the ionospheric gradient, which yields the dead zone of IGM, i.e., the amplitude of anomalous ionospheric gradient is integer multiplied by wavelength. Therefore, it can be found that the instantaneous IGM method is insufficient to satisfy the requirements of integrity and continuity for precise approach.

To improve the sensitivity to anomalous ionospheric gradients, Patel et al. [21] proposed a decoupling IGM method to avoid the dead zone of IGM. The ionosphere-free model is utilized for ambiguity resolution to ensure the anomalous ionospheric gradient is independent of the ambiguity. However, the reliability of ambiguity resolution is affected by the ionosphere-free combination noise. Therefore, the moving average method is applied in the decoupling IGM method to improve the reliability of ambiguity resolution. When the ambiguity is incorrectly resolved, the integer cycle bias is coupled with the test statistics, which yields the failure of continuity or integrity. Therefore, the reliability of ambiguity resolution directly affects the performance of the decoupling IGM.

It is worth noting that the over-long moving average length will impact the real-time performance of IGM during the precision approach. Moreover, the effect of ambiguity resolution failure caused by insufficient averaging length on the minimum detection error (MDE) is still ambiguous. The aim of this work was to reduce the MDE for ionospheric anomaly detection by allocating the probabilities of false alarm and missed detection among different ambiguity resolution modes. Furthermore, with the forthcoming multi-GNSS modernization, since GBAS approach service type F (GAST-F) has been designed for the provision of category III services using multi-constellation (MC) and multi-frequency (MF) corrections [22,23], we propose a dual-frequency carrier phase-based IGM method in order to improve the sensitivity of the IGM method and reduce the averaging length, simultaneously.

In this contribution, we developed an IGM method under the mixed Gaussian distribution. The detection threshold and MDE constraint method were analyzed by allowing acceptable ambiguity resolution failure modes. Then, the dual-frequency carrier phase was used to construct highly sensitive test statistics for IGM, and the MDE is compared with single- and dual-frequency IGM methods. Finally, we present the effectiveness of the proposed method in IGM by using BeiDou navigation satellite system (BDS) dual-frequency data and summarize the research findings.

This paper is organized as follows. The physical basis for the test statistics of the proposed IGM method is presented in Section 2. The modified MDE algorithm and the simulation result used in the study are presented in Section 3 and Section 4. The experiment results and discussion are given in Section 5.

2. Test Statistic for the IGM

The basic configuration for the proposed IGM method is the short baseline reference receivers involving code and phase observations. Both the code and phase observations can be used to construct the test statistics. Due to the relatively larger observation noise, the test statistics based on code observation have limitations in obtaining higher sensitivity of IGM. Therefore, the carrier phase-based IGM method has attracted great attention on monitoring the anomalous ionospheric gradients with high sensitivity. However, it is still a challenge to correctly resolve the ambiguity in the case of anomalous ionospheric gradient. We assume only one satellite m will be affected by anomalous ionospheric gradient at any given time. For a faulty satellite m and a fault-free satellite n, the test statistic q for the carrier phase-based IGM method can be computed by using double-difference phase observations as

$$q_{a,b,f}^{m,n}=\mathsf{\Delta }\nabla {\varphi }_{a,b,f}^{m,n}-\nabla {\mathit{e}}^T{\mathit{x}}_{a,b}^{}-{\lambda }_f\mathsf{\Delta }\nabla N_{a,b,f}^{m,n}={\mu }_f\mathsf{\Delta }\nabla I_{a,b}^{m,n}+\mathsf{\Delta }\nabla T_{a,b}^{m,n}+{\sigma }_{\mathsf{\Delta }\nabla \varphi }$$

(1)

where $\mathsf{\Delta }\nabla$ denotes the double-difference operator between the receivers a-b and satellites m-n, $\varphi$ is the phase observations in unit of meters, f is the frequency of phase observation, ${\mathit{e}}^T$ is the differential user satellite line of sight unit vector, $x$ is the baseline vector between reference receiver antennae, N_f is the unknown integer ambiguity with the corresponding wavelength ${\lambda }_f$, ${\mu }_f$ is the ratio of the ionospheric delay I at the f frequency to the delay at the L1 frequency, T is the tropospheric delay, and ${\sigma }_{\mathsf{\Delta }\nabla \varphi }$ is the random noise of normal distribution with zero mean [24].

It can be found from Equation (1) that the reliability of ambiguity resolution is of significance to obtain a high sensitivity in the carrier phase-based IGM method. The single-frequency ionosphere-free combination model is selected for ambiguity resolution to avoid the influence of the potential ionospheric anomaly on the ambiguity resolution [21,25]. In addition, in order to accurately estimate the success rate of the ambiguity resolution, we selected the rounding estimator for ambiguity resolution in this contribution [26,27,28]. Since the accuracy of ionosphere-free ambiguity float resolution is degraded by the combination observation noise, we utilized the moving averaging method to reduce the noise and achieve reliable ambiguity resolution.

When the ambiguity is correctly fixed, the test statistics can be utilized to monitor either the ionospheric or the tropospheric gradient anomalies. It has been widely proven that the ionospheric gradient and the tropospheric gradient have strong temporal and spatial correlation under the nominal atmospheric condition [29,30]. In category III of the precision approach, the horizontal distance between the aircraft decision height point and the glide slope intercept point (GPIP) is defined as 5km, which belongs to the case of short baseline [31]. Therefore, we ignored the influence of double-difference tropospheric residuals in the following analysis. When the observation noise is effectively suppressed, the test statistics can be highly sensitive to the anomalous ionospheric gradients.

3. Minimum Detectable Error

GBAS must monitor the quality of the observations to guarantee the integrity and continuity of the precision approach guidance service. Murphy et al. [15] confirmed that a hazardous situation occurs when the differential ranging errors between aircraft decision height point and the ground reference receiver antenna exceed 1.5 m. According to the ionospheric gradient model, when the horizontal baseline length between the aircraft-determined altitude point and the ground reference receiving antenna is 5 km, the local ground facility of GBAS needs to successfully monitor anomalous ionospheric gradients greater than 300 mm/km [11,16]. When the MDE of IGM is less than the limit, the local ground facility can guarantee the normal operation of landing guidance service with the acceptable false alarm and missed detection probability. Therefore, the MDE can be used to guarantee the availability of IGM when the GBAS is utilized for the safety-of-life navigation services.

The MDE of IGM is associated with the constructed test statistics. The reliability of ambiguity resolution directly affects the statistical distribution of test statistics. We define $\mathsf{\Delta }\nabla \overline{N}$ and $\mathsf{\Delta }\nabla N_i$ as the true double-difference ambiguity and the estimated double-difference ambiguity, respectively, where the subscripts i denote the integer estimation biases between $\mathsf{\Delta }\nabla \overline{N}$ and $\mathsf{\Delta }\nabla N_i$, i.e., $i=\mathsf{\Delta }\nabla \overline{N}-\mathsf{\Delta }\nabla N_i$. It can be found from Equation (1) that the test statistics follow the Gaussian distribution when the ambiguity is correctly resolved, i.e., i = 0. Otherwise, the test statistics follow the mixed Gaussian distribution when the ambiguity resolution is a failure, i.e., $i\ne 0$. We analyzed the MDE under the assumptions of Gaussian distribution and the mixed Gaussian distribution, respectively.

3.1. Gaussian Distribution

The hypotheses for the ionospheric gradients were set as fault-free hypothesis H₀ and the alternative faulty hypothesis H₁. Under the H₀ hypothesis, the false alarm error will present when the test statistics are beyond the detection threshold. Therefore, given the required false alarm probability P_FA, the detection threshold T can be determined as

$$T=k_{fa}{\sigma }_{\mathsf{\Delta }\nabla \varphi }$$

(2)

where k_fa is the false alarm factor, which is derived from a probability of false alarm requirement [7,31,32]. If the test statistic with respect to the nth satellite is beyond the detection threshold, the monitor must flag the nth satellite ionospheric as faulty. Otherwise, we evaluate the missed detection errors under the H₁ hypothesis.

The constructed test statistics follow non-zero mean Gaussian distribution in the presence of anomalous ionospheric gradients. When the detection threshold T is determined by Equation (2), the MDE can be computed as

$$MDE=(k_{fa}+k_{md}){\sigma }_{\mathsf{\Delta }\nabla \varphi }$$

(3)

where k_md is the missed detection factor derived from a probability of missed detection requirement [7,31,32]. It can be found from Equations (2) and (3) that the MDE is designed to satisfy the requirements of false alarm and missed detection probabilities, simultaneously. Therefore, the sensitivity of the IGM can be characterized by the MDE under the Gaussian distribution.

The single dimensional simulation was carried out to verify the MDE with different averaging length under the Gaussian distribution by varying the noise of double-difference observations, as shown in Figure 1. In addition, we varied the prior probability of the anomalous ionospheric gradient as 10⁻¹, 10⁻², 10⁻³, 10⁻⁴, and 10⁻⁵. Note that the GBAS continuity risk of 10⁻⁸ was equally allocated to each satellite as the predefined false alarm probability. The integrity risk for IGM was set to 10⁻⁹ for category III precision approach. It can be found that the MDE decreased for a smaller prior faulty probability. Since the determination of the prior probability of ionosphere anomaly was not the focus in this contribution, we selected 10⁻³ as the prior probability in the subsequent analysis [21]. It can be observed that both the numerical results of MDE under different observation noise were within the limit of 300mm/km, which meant that the IGM under the Gaussian distribution could be applied to protect against the ionospheric gradient. It should be noted that the extraordinary MDE was based on the assumption that the ambiguity was correctly resolved. However, the successful ambiguity resolution assumption was overoptimistic because the ratio of noise to wavelength was limited by the amplified observation noise from ionosphere-free combination. Furthermore, the rounding estimator for ambiguity resolution was undermined because the correlation between observations was ignored.

In order to reveal the mixed Gaussian distribution caused by the incorrect ambiguity resolution, we display the distribution of test statistics based on 2 days of real BDS data in Figure 2. It can be observed that the test statistics had a multi-peak feature that was inconsistent with the Gaussian distribution. Specifically, the peaks were concentrated at integers, such as −1, 0, and 1, which was mainly caused by the incorrect ambiguity resolution. Although the longer averaging length can make the Gaussian distribution assumption more reasonable, the over-long averaging length cannot satisfy the real-time requirements of IGM. Therefore, the MDE under the Gaussian distribution is not rigorous enough to express the sensitivity of the IGM, particularly in the presence of ambiguity resolution failure. It is necessary to discuss the MDE under the mixed Gaussian distribution.

3.2. Mixed Gaussian Distribution

There are theoretically infinite ambiguity resolution failure modes to obtain the distribution of test statistics. We denote the ambiguity resolution failure mode as F_i ($i\in Z$). It always holds true that ${\displaystyle \sum _{i\in Z}^{}P(F_i)}=1$, in which the corresponding prior probability of failure mode F_i can be expressed as [33,34]

$$P(F_i)={\displaystyle {\int }_{i-\frac{1}{2}}^{i+\frac{1}{2}}\frac{1}{\sqrt{2\pi }{\sigma }_{amb}}}\exp \left(-\frac{v^2}{2{\sigma }_{amb}^2}\right)dv=\mathsf{\Phi }\left(\frac{2i+1}{2{\sigma }_{amb}}\right)-\mathsf{\Phi }\left(\frac{2i-1}{2{\sigma }_{amb}}\right)$$

(4)

where $\mathsf{\Phi }\left(x\right)={\displaystyle {\int }_{-\infty }^x\frac{1}{\sqrt{2\pi }}\exp \left(-\frac{1}{2}{v}^2\right)}dv$ is the Gaussian cumulative probability function, and ${\sigma }_{amb}$ is the standard deviation of float ambiguity in cycles. Moreover, for each failure mode F_i, the test statistic follows a Gaussian distribution with the bias of $i\lambda$ [35]. Therefore, the distribution of the test statistic is defined as a linear combination of different independent Gaussian distributions, as shown in Figure 3, and the probability density function of the mixed Gaussian distribution can be expressed as

$$F\left(x\right)={\displaystyle \sum _{i\in Z}^{}f\left(x|F_i\right)}P(F_i)$$

(5)

where $f\left(x|F_i\right)=\frac{1}{\sqrt{2\pi }{\sigma }_{\mathsf{\Delta }\nabla \varphi }}\exp \left(-\frac{{\left(x-i\lambda \right)}^2}{2{\sigma }_{\mathsf{\Delta }\nabla \varphi }^2}\right)$ is the probability density function for failure mode F_i.

Unlike the Gaussian distribution, each ambiguity resolution failure mode F_i will cause a peak in the probability density function of the mixed Gaussian distribution. The multimodal feature brings challenges to the IGM because it is difficult to control the false alarm and missed detection error simultaneously when there are infinite ambiguity resolution failure modes to be considered. Since most failure modes have little impact on the overall distribution, we only therefore show the failure modes from F₋₄ to F₄ in Figure 3. Under the ionospheric fault-free hypothesis, when the test statistics are beyond the detection threshold, the probability of false alarm can be computed as

$$P_{FA}=P\left(\left|q\right|>T|H_0\right)=2{\displaystyle {\int }_{-\infty }^{-T}F\left(x|H_0\right)dx}={\displaystyle \sum _{i\in Z}^{}{P}_{FA|F_i}}P(F_i)$$

(6)

where Z represents the set of integers, and $P_{FA|F_i}=\mathsf{\Phi }\left(\frac{-T-i\lambda }{{\sigma }_{\mathsf{\Delta }\nabla \varphi }}\right)+1-\mathsf{\Phi }\left(\frac{T-i\lambda }{{\sigma }_{\mathsf{\Delta }\nabla \varphi }}\right)$ represents the conditional false alarm probability in the presence of the ambiguity resolution failure mode $F_i$. Compared with the Gaussian distribution, it can be noted that the probability of false alarm is the result of each conditional false alarm probability. It can be induced that it is practically not possible to integrate all the series terms in Equation (6) to control the false alarm errors under the mixed Gaussian distribution. Therefore, the effect of each ambiguity resolution failure mode $F_i$ on the false alarm event FA needs be clarified to construct a tight boundary on the probability of false alarm.

When the false alarm error is flagged, the probability of the false alarm error caused by failure mode F_i can be computed as

$$P\left(F_i|FA\right)=\frac{P_{FA|F_i}P(F_i)}{{\displaystyle \sum _{i\in Z}^{}{P}_{FA|F_i}}P(F_i)}=\frac{P_{FA|F_i}P(F_i)}{P_{FA}}$$

(7)

where P_FA is the required false alarm probability. It can be found that the $P\left(F_i|FA\right)$ is affected by both the conditional false alarm probability $P_{FA|F_i}$ and the a priori ambiguity resolution failure probability $P(F_i)$. It is obvious that the observation noise plays an important role in computing the $P\left(F_i|FA\right)$. Therefore, in order to reveal the contribution of each mode F_i to the false alarm probability, we varied the averaging length from 200 to 1600 epochs to indicate different combined observation noise levels. The standard deviations of double-difference observations noises for carrier phase and code were conservatively set as 1cm and 1m, respectively. Due to the symmetry feature of the zero-mean Gaussian distribution under H₀ hypothesis, we found that $P_{FA|F_i}=P_{FA|F_{-i}}$ and $P(F_i)=P(F_{-i})$. The numerical probability of the false alarm error caused by mode $F_{\pm i}$ is shown in

Table 1. Contribution of mode F_i to the overall false alarm event.

Averaging Length	200	400	600	800	1000	1200	1400	1600
$P\left(F_0|FA\right)$	0	0	0	0	0	0	0	0
$P\left(F_{\pm 1}|FA\right)$	0	0	0	0	0	0	0	91.83%
$P\left(F_{\pm 2}|FA\right)$	0	0	96.68%	100%	100%	100%	100%	8.17%
$P\left(F_{\pm 3}|FA\right)$	0	100%	3.32%	0	0	0	0	0
$P\left(F_{\pm 4}|FA\right)$	99.45%	0	0	0	0	0	0	0
$P\left(F_{\pm 5}|FA\right)$	0.55%	0	0	0	0	0	0	0
$\cdots$	0	0	0	0	0	0	0	0

. We denoted an ambiguity failure mode $F_{i_{FA}}$ such that $P\left(F_{\pm i_{FA}}|FA\right)$ was nearly 100% as the dominate mode to trigger the false alarm errors. For example, when the averaging length was greater than 600 and less than 1400 epochs, $i_{FA}$ was determined as 2. Furthermore, it can be found that for the failure modes $i<i_{FA}$, the probability of false alarm can barely be affected by the mode $F_{\pm i}$ because of $\frac{P_{FA|F_{\pm i}}P(F_{\pm i})}{P_{FA}}\ll {10}^{-2}$, and thus we can ignore the impact of failure mode $F_{\pm i}$ on the overall false alarm error. Therefore, the probability of false alarm can be simplified as

$$P_{FA}={\displaystyle \sum _{i\in N}^{}{P}_{FA|F_{\pm i}}}P(F_{\pm i}),\left(i\ge i_{FA}\right)$$

(8)

where N represents the set of natural numbers. On the other hand, for the ambiguity resolution failure mode $i>i_{FA}$, we conservatively assumed that the failure modes needed be accounted for in terms of the false alarm errors, i.e., $P_{FA|F_{\pm i}}=1$. Therefore, the false alarm probability can be upper bounded as

$$\begin{array}{l}{P}_{FA}=P_{FA|F_{\pm i_{FA}}}P(F_{\pm i_{FA}})+{\displaystyle \sum _{i\in N,i>i_{FA}}^{}{P}_{FA|F_{\pm i}}P(F_{\pm i})}\le P_{FA|F_{\pm i_{FA}}}P(F_{\pm i_{FA}})+{\displaystyle \sum _{i\in N,i>i_{FA}}^{}P(F_{\pm i})}\\ =P_{FA|F_{\pm i_{FA}}}P(F_{\pm i_{FA}})+1-{\displaystyle \sum _{i=0}^{i_{FA}}P(F_{\pm i})}\end{array}$$

(9)

Then, the conservative detection threshold $T_{mix}$ under the mixed Gaussian distribution can be computed by bringing Equation (9) into Equation (2):

$$T_{mix}=i_{FA}\lambda +k_{fa\_mix}{\sigma }_{\mathsf{\Delta }\nabla \varphi }$$

(10)

where $i_{FA}\lambda$ is the bias caused by the ambiguity resolution failure mode $i_{FA}$, and $k_{fa\_mix}$ is the false alarm quantile under the mixed Gaussian distribution,

$$k_{fa\_mix}={\mathsf{\Phi }}^{-}\left(1-\frac{P_{FA}+{\displaystyle \sum _{i=0}^{i_{FA}}P(F_{\pm i})}-1}{P(F_{\pm i_{FA}})}\right){\sigma }_{\mathsf{\Delta }\nabla \varphi }$$

(11)

Since the dominant mode depends on the performance of ambiguity resolution, we found that the detection threshold $T_{mix}$ was sensitive to averaging length. The ionospheric gradient of this satellite pair can be flagged as anomalous when the test statistic is beyond the detection threshold. The specific ionospheric anomaly satellite can be further indicated by fault detection of all test statistics. Otherwise, the probability of missed detection can be computed when the test statistics are within the protection of detection threshold. With the detection threshold $T_{mix}$ computed from the false alarm error constraint, the resulting probability of missed detection P_MD can be computed as

$$P_{MD}=P\left(\left|q\right|<T|H_1\right)={\displaystyle {\int }_{-T}^{T}F\left(x|H_1\right)dx}={\displaystyle \sum _{i\in Z}^{}{P}_{MD|F_i}}P(F_i)$$

(12)

where $P_{MD|F_i}=\mathsf{\Phi }\left(\frac{T-\mu }{{\sigma }_{\mathsf{\Delta }\nabla \varphi }}\right)-\mathsf{\Phi }\left(\frac{-T-\mu }{{\sigma }_{\mathsf{\Delta }\nabla \varphi }}\right)$ represents the conditional missed detection probability under the ambiguity resolution failure mode F_i, and $\mu$ represents the unknown biases caused by the coupling of the incorrect ambiguity resolved and anomalous ionospheric gradient. The missed detection error can be constrained by comparing the probability of the missed detection computed by Equation (12), with the required probability of missed detection. Furthermore, the probability of the missed detection error caused by the failure mode F_i can be computed as

$$P\left(F_i|MD\right)=\frac{P_{MD|F_i}P(F_i)}{{\displaystyle \sum _{i\in Z}^{}{P}_{MD|F_i}}P(F_i)}=\frac{P_{MD|F_i}P(F_i)}{P_{MD}}$$

(13)

The numerical probability of $P\left(F_i|MD\right)$ by varying averaging length is shown in

Table 2. Contribution of mode F_i to the overall missed detection event.

Averaging Length	200	400	600	800	1000	1200	1400	1600
$P\left(F_0|MD\right)$	0	0	0	0	0	0	0	0
$P\left(F_1|MD\right)$	0	0	0	0	86.11%	98.24%	99.88%	99.99%
$P\left(F_2|MD\right)$	0	93.88	100%	100%	13.89%	1.76%	0.12%	0.01%
$P\left(F_3|MD\right)$	93.04%	3.06%	0	0	0	0	0	0
$P\left(F_4|MD\right)$	6.96%	3.06%	0	0	0	0	0	0
$\cdots$	0	0	0	0	0	0	0	0

. When the contribution of failure mode F_i to the overall missed detection event MD is nearly 100%, it can be also denoted as the dominate mode $F_{i_{MD}}$. When the failure mode $i<i_{MD}$, we found that the missed detection probability is barely be affected by the mode $F_i$ because of $\frac{P_{MD|F_i}P(F_i)}{P_{MD}}\ll {10}^{-2}$, and thus we can also ignore the impact of mode F_i on the missed detection probability. Therefore, the missed detection probability P_MD can be simplified as

$$P_{MD}={\displaystyle \sum _{i\in N}^{}{P}_{MD|F_i}}P(F_i),\left(i\ge i_{MD}\right)$$

(14)

Similarly, for the failure mode $i>i_{MD}$, we conservatively assumed that incorrect ambiguity resolution will inevitably lead to missed detection errors. Therefore, the resulting missed detection probability can be upper bounded as

$$P_{MD}\le P_{MD|F_{i_{MD}}}P(F_{i_{MD}})+1-{\displaystyle \sum _{i=0}^{i_{MD}}P(F_i)}$$

(15)

Then, the conservative MDE under the mixed Gaussian distribution can be computed as

$${MDE}_{mix}=i_{MD}\lambda +\left(k_{fa\_mix}+k_{md\_mix}\right){\sigma }_{\mathsf{\Delta }\nabla \varphi }$$

(16)

where $i_{MD}\lambda$ is the biases caused by ambiguity resolution failure mode $i_{MD}$, and $k_{md\_mix}$ is the missed detection quantile under mixed Gaussian distribution, which can be computed as

$$k_{md\_mix}={\mathsf{\Phi }}^{-1}\left(1-\frac{P_{MD}+{\displaystyle \sum _{i=0}^{i_{MD}}P(F_i)}-1}{P(F_{i_{MD}})}\right)$$

(17)

We found that when compared with Gaussian distribution, part of the sensitivity performance under the mixed Gaussian distribution was lost by allocating the probabilities of false alarm and missed detection among different ambiguity resolution modes. In return, the requirement for the success rate of ambiguity resolution was reduced. Moreover, it should be noted that when a long enough averaging length is used for ambiguity resolution, the MDE under the mixed Gaussian distribution can achieve comparable performance to the Gaussian distribution. Therefore, the counterbalance between the averaging length and the MDE performance is critical to improve the IGM performance.

4. Benefits of Dual Frequency

We proved that the MDE performance under the mixed Gaussian distribution depends on the reliability of ambiguity resolution. With dual-frequency signals available for civil aviation, we proposed the dual-frequency IGM method to release the limitation of the averaging length in the MDE. L2 signal is more frequently loss of locks than L1 signal, particularly at high latitudes [36,37]. The signals acquisition of L1 and L2 frequencies simultaneously is essential to obtain the test statistics for the proposed IGM method. However, the signal acquisition is not the focus of this contribution. We assumed the availability of the dual-frequency observations. To be immune from the influence of potential ionospheric anomalies on ambiguity resolution, we proposed a two-step ionosphere-free method for ambiguity resolution. Firstly, the wide-lane (WL) ambiguity was estimated by using the double-difference Melbourne–Wübbena (MW) combination,

$$\mathsf{\Delta }\nabla {\widehat{N}}_{WL}=\frac{\left(\frac{f_1\mathsf{\Delta }\nabla {\varphi }_1-f_2\mathsf{\Delta }\nabla {\varphi }_2}{f_1-f_2}\right)-\left(\frac{f_1\mathsf{\Delta }\nabla p_1+f_2\mathsf{\Delta }\nabla p_2}{f_1+f_2}\right)}{{\lambda }_{wl}}$$

(18)

where $\mathsf{\Delta }\nabla {\widehat{N}}_{WL}$ is the float solution of WL ambiguity, and the wavelength of WL observation ${\lambda }_{wl}=\frac{c}{f_1-f_2}$, with the uncertainty of the $\mathsf{\Delta }\nabla {\widehat{N}}_{WL}$ being computed as

$${\sigma }_{\mathsf{\Delta }\nabla {\widehat{N}}_{WL}}=\frac{1}{{\lambda }_{wl}^{}}\sqrt{\left(\frac{f_1^2+f_2^2}{{\left(f_1-f_2\right)}^2}\right){\sigma }_{\mathsf{\Delta }\nabla \varphi }^2+\left(\frac{f_1^2+f_2^2}{{\left(f_1+f_2\right)}^2}\right){\sigma }_{\mathsf{\Delta }\nabla p}^2}$$

(19)

where ${\sigma }_{\mathsf{\Delta }\nabla \varphi }^2$ and ${\sigma }_{\mathsf{\Delta }\nabla p}^2$ are the variance of double-difference carrier phase and code observations, respectively. Secondly, the ionosphere-free combination observations can be computed as

$$\mathsf{\Delta }\nabla {\varphi }_{IF}=\frac{f_1^2\mathsf{\Delta }\nabla {\varphi }_1-f_2^2\mathsf{\Delta }\nabla {\varphi }_2}{f_1^2-f_2^2}$$

(20)

where the wavelength of ionosphere-free observation ${\lambda }_{IF}=\frac{c}{f_1+f_2}$. The narrow-lane (NL) ambiguity can therefore be estimated as

$$\mathsf{\Delta }\nabla {\widehat{N}}_1=\frac{\mathsf{\Delta }\nabla {\varphi }_{IF}-\frac{c{f}_2}{f_1^2-f_2^2}\mathsf{\Delta }\nabla N_{WL}-\nabla {\mathit{e}}^T\mathit{x}}{{\lambda }_{IF}}$$

(21)

where the noise of the $\mathsf{\Delta }\nabla {\widehat{N}}_1$ can be computed as

$${\sigma }_{\mathsf{\Delta }\nabla {\widehat{N}}_1}=\frac{1}{{\lambda }_{IF}^{}}\sqrt{\frac{f_1^4+f_2^4}{{\left(f_1^2-f_2^2\right)}^2}}{\sigma }_{\mathsf{\Delta }\nabla \varphi }^{}$$

(22)

Since the precise coordinates of reference station antennae are known accurately in the local ground facility of GBAS, we could resolve the integer ambiguity by the rounding estimator when the noise of $\mathsf{\Delta }\nabla {\widehat{N}}_{WL}$ and $\mathsf{\Delta }\nabla {\widehat{N}}_1$ were effectively suppressed by the moving averaging method. Then, the test statistics could be obtained with the fixed NL ambiguity.

In order to verify the sensitivity improvement of dual-frequency IGM method, we conducted a single-dimensional simulation to compare the MDE between the single- and the dual-frequency IGM method, as shown in Figure 4. We conservatively assumed that the double-difference carrier phase observation error follows the zero-mean Gaussian distribution with a standard deviation of 1 cm. The double-difference code observation error followed the zero-mean Gaussian distribution with a standard deviation of 1 m. The baseline length was selected as 1 km. The required false alarm probability was set to 10⁻⁸ to satisfy the continuity requirement of GBAS category III. The prior probability of ionospheric anomaly was set up as 10⁻³. The resulted MDE was evaluated by varying the averaging length.

It can be seen from Figure 4 that both the MDE of single- and dual-frequency IGM method cannot satisfy the requirement of GBAS approach service type D category III when the averaging length was relatively short. In addition, we observed that the MDE increased with the averaging length. Compared with the single-frequency IGM method, the comparable MDE can be obtained with a shorter averaging length for the dual-frequency IGM method. The simulation result demonstrates the effectiveness of the dual-frequency IGM on achieving a better counterbalance between the averaging length and the MDE. A better integrity monitoring performance of ionospheric gradient can therefore be anticipated.

5. Experiment and Discussion

Since the MDE indicates the sensitivity of test statistics to anomalous ionospheric gradients, we tested the proposed IGM method under the fault-free and the ionospheric anomaly conditions. The dual-frequency BDS observations from 19 March 2017 to 20 March 2017 were collected from https://www.geodetic.gov.hk/sc/satref/SatRef.htm. The sampling period was 30 s, with a total of 5760 epochs. Two international GNSS service stations HKWS and HKSS were selected to form a local ground facility. The cutoff elevation angle was set as 15 degrees. The observation noise model adopted the elevation-dependent weighting, i.e., $\sigma (\theta )={\sigma }_0\left(1+1/\sin \theta \right)$. The required false alarm and missed detection probabilities were set up to 10⁻⁸ and 10⁻⁶, respectively. The cycle slip detection was implemented with the multiple-epochs ambiguity consistency check [38].

In order to verify that the anomalous atmospheric gradient was absent for the selected data, we calculated the double-difference carrier phase residuals by real-time kinematic (RTK) algorithms, as shown in Figure 5. Since the double-difference carrier phase residuals of each satellite-pairs had zero mean characteristics, we induced that the carrier phase residuals could be assumed to follow the zero mean Gaussian distribution as the observation noise. Therefore, we can prove that the anomalous atmospheric gradient was absent for the selected short-baseline dataset.

The averaging length was related to the IGM performance because the ambiguity resolution success rate was affected by averaging length. The reference for the proposed ambiguity resolution came from the output of the LAMBDA algorithm. The ambiguity resolution success rates based on the ionosphere-free combination is shown in Figure 6. Since the multiple epochs averaging method was beneficial in reducing the observation noise, we found that when the averaging length was less than 300 epochs, both the single- and dual-frequency ambiguity resolution success rates were increased with averaging length. Specifically, the best ambiguity resolution performance of the dual-frequency ionosphere-free model was obtained when the averaging length was 300 epochs. Moreover, since the cycle slip detection with the multiple-epoch ambiguity consistency check increased the probability of averaging reinitialization, we have shown that the ambiguity resolution performance deteriorated when the averaging length exceeded 300 epochs. Therefore, we chose the averaging length as 300 epochs in the following experiments.

We first evaluated the performance improvement of the proposed dual-frequency IGM method in the fault-free case. Two groups of test statistics constructed by single- and dual-frequency IGM methods are given in Figure 7. It can be observed that there were many cases wherein the test statistics of the single-frequency IGM method were beyond the protection of the detection threshold. This can be attributed to the incorrect ambiguity fixing. In addition, we that the detection thresholds of dual-frequency IGM method were stricter than the single-frequency IGM because there were less ambiguity resolution failure modes to be accounted for when the proposed dual-frequency ionosphere-free ambiguity resolution method was utilized. Moreover, we observed that the false alarm error of the dual-frequency IGM method was absent. It can be induced that the proposed dual-frequency IGM method can improve the false alarm performance under the fault-free case.

There are two groups of MDE computed by single- and dual-frequency IGM method displayed in Figure 8. It can be observed that the single-frequency MDE was larger than 300 mm/km; therefore, the single-frequency IGM method cannot satisfy the minimum acceptable ionospheric gradient requirement of the GBAS category III precision approach. In contrast, it can be found that the dual-frequency IGM method can monitor the anomalous ionospheric gradient larger than 200 mm/km. It has been proven that the proposed dual-frequency IGM method has higher sensitivity to the ionospheric gradient.

Secondly, we tested the IGM in the case of ionospheric anomaly. It should be noted that the instantaneous IGM method from Khanafseh and Pullen (2012) was chosen to reveal the sensitivity improvement of the proposed IGM method, while the single-frequency IGM method was used to compare the improvement of the proposed IGM method in the required averaging length. Since we verified that the ionosphere gradient was absent for the selected dataset, we simulated a virtual ionospheric gradient from 60 mm/km to 300 mm/km to verify the MDE of the three IGM methods to the ionospheric gradients. The virtual ionospheric gradient was injected into the reference satellite to affect all the test statistics. Moreover, we selected the geostationary earth orbit satellite C01 as the reference satellite to ensure the consistency of experimental conditions. The state of simulated virtual ionosphere gradient injection and the detection process of three IGM methods are shown in Figure 9.

As seen from Figure 9, when the virtual ionosphere gradient was injected, both the three IGM methods can respond immediately. It can be observed that there were more observation noises in the test statistics of instantaneous IGM methods, which limited the detection power of the instantaneous IGM method. Therefore, we found that some of the test statistics were within the protection of the detection threshold when the virtual ionosphere gradient was injected. In addition, we found that although the virtual ionosphere injection caused the test statistics of single-frequency IGM method to fluctuate, most of the test statistics were still within the protection of detection threshold. The main reason was that the detection threshold was too loose due to accounting for too many ambiguity resolution failure modes. In contrast, the test statistics of dual-frequency IGM method were obviously beyond the protection of the detection threshold when the virtual ionospheric gradient was injected. It can be demonstrated that neither the instantaneous nor the single-frequency IGM method had limitation in monitoring the small-scale ionospheric gradient. Moreover, we proved that the dual-frequency IGM method had higher sensitivity to the small-scale ionospheric gradients with a shorter averaging length.

The numerical results of false alarm and missed detection probability shown in

Table 3. P_FA performance comparison without ionospheric gradients.

	Instantaneous IGM	SF-IGM	DF-IGM
PFA (%)	0	0.19	0

and

Table 4. Comparison of detection sensitivity in the presence of ionospheric gradients.

Ionospheric Gradients (mm/km)	PMD (%)
	Instantaneous IGM	SF-IGM	DF-IGM
60	94.55	1	15.64
90	55.27	98.54	1.82
120	11.27	89.45	1.09
150	70.55	78.55	0.72
180	89.82	44.36	0.01
210	76.36	35.27	0
240	48.00	26.91	0
270	35.27	14.18	0
300	9.82	4.36	0

were used to quantitatively verify the MDE of the proposed IGM methods. When the numerical results of false alarm and missed detection probability satisfied the requirements of integrity and continuity, the MDE of the IGM methods could simultaneously be verified. It should be noted that the P_FA was determined by the ratio of the epochs, showing that the test statistics were beyond the detection threshold and the total epochs. On the other hand, the P_MD was determined by the ratio of the epochs, showing that the resulted missed detection probability was beyond the required missed detection probability and the total epochs of virtual ionospheric gradient being injected.

For the instantaneous IGM method, we observed that although the false alarm performance could be successfully constrained, the probability of missed detection could hardly be reduced with the increase of the virtual ionospheric gradient. This was because the insensitivity originated from the dead zone of instantaneous IGM method. It can be induced that more baselines are required to improve the MDE of the instantaneous IGM method. In addition, we found that the single-frequency IGM method can monitor the anomalous ionospheric gradients with unacceptable false alarm and missed detection case, even if the virtual ionospheric gradient reached 300 mm/km, which is consistent with the numerical MDE results. It should be noted that longer averaging length is needed to improve MDE of single-frequency IGM method. Since the large anomalous ionospheric gradient can easily be monitored by the IGM, we found that the missed detection resulting from dual-frequency IGM method decreased gradually with the increased virtual ionospheric gradients. Specially, when the virtual ionospheric gradient reached 210 mm/km, the dual-frequency IGM method could successfully monitor the anomalous ionospheric gradients with the required false alarm and missed detection probability, which indicates that the MDE resulting from the dual-frequency IGM method can be reduced to 210 mm/km given the averaging length. This is 30% less than the maximum acceptable anomalous ionospheric gradients for category III GBAS. The comparison results indicated that the proposed dual-frequency IGM method can provide more of a monitoring margin for category III GBAS. Furthermore, in order to compare the sensitivity between the dual-frequency and the single-frequency IGM methods, we further increased the injected virtual ionospheric gradient for the single-frequency IGM method. When the virtual ionospheric gradient reached 540mm/km, we found that the single-frequency IGM method could satisfy the required probability of missed detection. It revealed that the 61.1% sensitivity improvement of dual-frequency IGM method can be achieved, given the required averaging length when compared with the single-frequency IGM method. The results can also indicate that the proposed dual-frequency IGM method achieved a better counterbalance between the averaging length and MDE.

6. Concluding Remarks

We proposed a high sensitivity carrier phase-based IGM method to detect the anomalous ionospheric gradients. On the basis of allocating the probabilities of false alarm and missed detection among different ambiguity resolution modes, we characterized the sensitivity of carrier phase-based IGM method by the computation of MDE under the mixed Gaussian distribution. The sensitivity of anomalous ionospheric gradient monitoring is therefore improved by allowing acceptable ambiguity resolution failure modes. Moreover, the averaging length of the ambiguity resolution for the dual-frequency carrier phase-based IGM method is greatly reduced when compared with single-frequency observations.

On the basis of the dual-frequency BDS real data, we assessed the reliability of ambiguity resolution for single- and dual-frequency IGM methods with different averaging lengths. The false alarm performance of the proposed carrier phase-based IGM method was verified under fault-free conditions. The experimental result indicates that the proposed dual-frequency IGM method can satisfy the false alarm requirement with a shorter averaging length. The sensitivity improvement of the proposed IGM method was verified by injecting the virtual ionospheric gradient. The numerical results showed that the proposed dual-frequency IGM method can be sensitive to the virtual ionospheric gradient of 210 mm/km given the required probabilities of false alarm and missed detection, which is 30% less than the limit of category III GBAS. In contrast, the single-frequency IGM method could achieve comparable performance when the injected ionospheric gradient was 540 mm/km. It was demonstrated that the proposed dual-frequency IGM method achieved a better counterbalance performance between the MDE and the averaging length. The proposed dual-frequency IGM methods can therefore provide real-time protection against the anomalous ionospheric gradient for GBAS category III services.