Satellite integrity monitoring for satellite-based augmentation system: an improved covariance-based method – Satellite Navigation

ByShuaiyong Zheng, Mengzhi Gao, Zhigang Huang, Xiaoqin Jin and Kun Li

May 9, 2022

In this section, a scale model and a groove model are developed for SBAS to perform satellite integrity monitoring.

Scale model based on multiple error sources

The pseudorange correction error ({text{d}}rho), namely User Equivalent Range Error (UERE), for a specific satellite is computed by

$${text{d}}rho = Delta rho – Delta {varvec{R}} cdot {varvec{I}} + Delta B$$

(2)

where the variables (Delta rho ,Delta {varvec{R}},Delta B) represent the synchronized pseudorange residual from monitor stations, the total long-term corrections computed by long-term satellite error corrections, and the total fast corrections computed by fast corrections and range-rate corrections, respectively. The 4 × 1 vector ({varvec{I}}) consists of a unit vector and one element − 1. The first three terms of this vector are the components of the unit vector along the line of sight in Earth-Centered Earth-Fixed (ECEF) coordinates (Wu & Peck, 2002).

Specifically, the pseudorange residual in (2) due to the ephemeris and clock for a satellite-user pair with the line of sight ({varvec{I}}) is theoretically given by (Blanch et al., 2012)

$${text{d}}rho_{{{text{SCE}}}} = {varvec{I}}^{{text{T}}} left( {{varvec{x}}_{{{text{BD}}}} – {varvec{x}}} right)$$

(3)

where ({varvec{x}}) and ({varvec{x}}_{{{text{BD}}}}) represent the true clock-ephemeris, and the clock-ephemeris computed by broadcast ephemeris and SBAS corrections, respectively. UDRE along with MT28 is the upper bound on this pseudorange residual for each satellite-user pair. According to SC-159 (2016), the error bound is expressed in the following form (Blanch et al., 2014)

$$L = K_{7} sigma_{{{text{flt}}}} = K_{7} sigma_{{{text{UDRE}}}} sqrt {{varvec{I}}^{{text{T}}} {varvec{C}}_{{{text{MT2}}8}}^{{{text{ov}}}} {varvec{I}}}$$

(4)

where the matrix ({varvec{C}}_{{{text{MT2}}8}}^{{{text{ov}}}}) is a 4 by 4 matrix. The parameter (K_{7}) denotes the quantile of 99.99999%. The parameter (sigma_{{{text{flt}}}}) denotes the clock-ephemeris standard deviation determined by UDRE and MT28.

To deduce the error bound for clock-ephemeris, there are four cases to be considered (Blanch et al., 2012): (1) nominal errors from the monitoring network receivers; (2) nominal biases or antenna biases; (3) satellite correction errors; (4) possibly undetected errors in the monitoring network receivers where one station is assumed to return erroneous measurements.

As for case (1), the relationship between the bound on the estimation error and the probability (P_{{{text{HMI}}}}) of Hazardously Misleading Information (HMI) is described by

$$Pleft( {left| {{varvec{I}}^{{text{T}}} left( {hat{user2{x}} – {varvec{x}}} right)} right| ge K_{{{text{HMI}}}} sqrt {{varvec{I}}^{{text{T}}} {varvec{P}}_{{{text{SCE}}}} {varvec{I}}} } right) = P_{{{text{HMI}}}}$$

(5)

where (K_{{{text{HMI}}}}) represents the quantile related to the probability (P_{{{text{HMI}}}}) and ({varvec{P}}_{{{text{SCE}}}}) denotes the covariance for the state (hat{user2{x}}) of clock-ephemeris. The probability (P_{{{text{HMI}}}}) is determined by an integrity allocation strategy (Lu et al., 2021; SC-167, 1992; Schempp et al., 2001; Wu & Peck, 2002).

According to (5), the error bound on the estimation error under the nominal conditions is given by

$$L_{1} = K_{{{text{HMI}}}} sqrt {{varvec{I}}^{{text{T}}} {varvec{P}}_{{{text{SCE}}}} {varvec{I}}}$$

(6)

The bound (L_{1}) is one upper bound on clock-ephemeris errors before the clock-ephemeris covariance matrix is broadcast. The nominal biases ({varvec{b}}) such as antenna biases and Code Noise and Multi-Path (CNMP) termed as the case (2) is described as a Gaussian vector with expectation (overline{user2{b}}) and covariance ({varvec{W}}^{ – 1}). As for the line of sight, the contribution of these biases can be given by ({varvec{I}}^{{text{T}}} {varvec{Hb}}). An upper bound of this variable is deduced by (Blanch et al., 2012)

$$K_{{{text{bias}}}} = mathop {max }limits_{{varvec{b}}} sqrt {{varvec{b}}^{T} {varvec{Wb}}} = sqrt {{varvec{b}}_{max }^{{text{T}}} {varvec{Wb}}_{max } }$$

(7)

The parameter (K_{{{text{bias}}}}) associated with antenna biases is calculated in real time as a function of its covariance and the maximum biases (Shallberg & Sheng, 2008). Therefore, when case (2) is considered, the error bound on the estimation error is adjusted by

$$L_{2} = left( {K_{{{text{HMI}}}} + K_{{{text{bias}}}} } right)sqrt {{varvec{I}}^{{text{T}}} {varvec{P}}_{{{text{SCE}}}} {varvec{I}}}$$

(8)

where (L_{2}) denotes the bound for clock-ephemeris errors under cases (1) and (2).

When the error from the broadcast clock-ephemeris or the quantization error of satellite corrections termed as the case (3) is considered, the upper bound for this error can be deduced by the following inequality

begin{aligned} left| {I^{T} left( {x_{BD} – hat{x}} right)} right| & = left| {I^{T} P_{SCE}^{0.5} P_{SCE}^{ – 0.5} left( {x_{BD} – hat{x}} right)} right| \ & le sqrt {left( {x_{BD} – hat{x}} right)^{T} P_{SCE}^{ – 1} left( {x_{BD} – hat{x}} right)} sqrt {I^{T} P_{SCE} I} \ end{aligned}

(9)

Let

$$K_{{{text{pfa}}}} = mathop {max }limits_{i} sqrt {left( {{varvec{x}}_{{{text{BD}},i}} – hat{user2{x}}_{i} } right)^{{text{T}}} {varvec{P}}_{{{text{SCE}}}}^{{ – {1}}} left( {{varvec{x}}_{{{text{BD}},i}} – hat{user2{x}}_{i} } right)}$$

(10)

The designator i represents the epoch of each update interval. When case (3) is taken into consideration, the error bound on the estimation error is updated by

$$L_{3} = left( {K_{{{text{HMI}}}} + K_{{{text{bias}}}} + K_{{{text{pfa}}}} } right)sqrt {{varvec{I}}^{{text{T}}} {varvec{P}}_{{{text{SCE}}}} {varvec{I}}}$$

(11)

Obviously, the new error bound satisfies

$$Pleft( {left| {{varvec{I}}^{{text{T}}} left( {hat{user2{x}} – {varvec{x}}} right)} right| ge left( {K_{{{text{HMI}}}} + K_{{{text{bias}}}} + K_{{{text{pfa}}}} } right)sqrt {{varvec{I}}^{{text{T}}} {varvec{P}}_{{{text{SCE}}}} {varvec{I}}} } right) le P_{{{text{HMI}}}}$$

(12)

which is in accordance with the constraint (5).

As for case (4), the reliability of each monitor station can be guaranteed by checking the observations of multi-set receivers which refers to a problem of data quality monitoring or system reliability (Hamada, 2008; Yin & Chai, 2020). As for WAAS monitor stations, three sets of receivers are equipped to collect observations (Parkinson et al., 1996). The observations from three receivers are used to conduct the cross check among three threads to remove erroneous observations (Parkinson et al., 1996; Shallberg & Sheng, 2008). Moreover, dual frequency pseudoranges can be smoothed by dual frequency carriers, and the related algorithm like IFree filter is chosen to improve the quality of the observations simultaneously with ionospheric delay removed (Hwang et al., 1999; Konno et al., 2006).

Before message type 28 is broadcast, the covariance matrix of clock-ephemeris needs to be adjusted to meet the integrity requirement. According to (4) and (11), the adjusted covariance matrix can be obtained by

$${varvec{P}}_{1} = left( {frac{{L_{3} }}{L}} right)^{2} {varvec{P}}_{{{text{SCE}}}} = left( {frac{{K_{{{text{HMI}}}} + K_{{{text{bias}}}} + K_{{{text{pfa}}}} }}{{K_{7} }}} right)^{2} {varvec{P}}_{{{text{SCE}}}}$$

(13)

The matrix denotes the clock-ephemeris covariance considering four cases.

After the message type 28 is available, the quantization error related to this message type needs to be protected either by the term (varepsilon_{C}) or by increasing the broadcast UDRE (Walter et al., 2001). Computing the theoretically optimal matrix ({varvec{P}}_{{{text{brdc}}}}) refers to a complicated mathematical problem, and a near-optimal method is adopted in practice. The authors choose a scaling method where the matrix ({varvec{P}}_{{{text{brdc}}}}) is obtained by finding a suitable shape covariance matrix and then scaling it so that the integrity safety condition is met (Walter et al., 2001; Wu & Peck, 2002). The integrity condition to be satisfied is

$$sigma_{{{text{UDRE}}}} delta_{{{text{UDRE}}}} ge sqrt {{varvec{I}}^{{text{T}}} {varvec{P}}_{1} {varvec{I}}}$$

(14)

$$frac{{sigma_{{{text{UDRE}}}} }}{{sqrt {P_{max } } }} ge frac{{sqrt {{varvec{I}}^{{text{T}}} {varvec{P}}_{1} {varvec{I}}} }}{{delta_{{{text{UDRE}}}} sqrt {P_{max } } }} = frac{{sqrt {{varvec{I}}^{{text{T}}} {varvec{CI}}} }}{{delta_{{{text{UDRE}}}} }} = frac{{sqrt {{varvec{I}}^{{text{T}}} {varvec{CI}}} }}{{sqrt {{varvec{I}}^{{text{T}}} {varvec{C}}_{{{text{brdc}}}} {varvec{I}}} + varepsilon_{C} }}$$

(15)

Then, the covariance matrix can be adjusted again by

$${varvec{P}}_{2} = left( {frac{{sigma_{{{text{UDRE}}}}^{{2}} }}{{P_{max } }}} right)left( {frac{{K_{{{text{HMI}}}} + K_{{{text{bias}}}} + K_{{{text{pfa}}}} }}{{K_{7} }}} right)^{2} {varvec{P}}_{{{text{SCE}}}}$$

(16)

where the parameter Pmax is deduced by

$$P_{max } = mathop {max }limits_{{{varvec{I}} in {varvec{I}}_{{{text{sv}}}} }} left( {{varvec{I}}^{{text{T}}} {varvec{P}}_{{{text{SCE}}}} {varvec{I}}} right) = {varvec{I}}_{{{text{WUL}}}}^{{text{T}}} {varvec{P}}_{{{text{SCE}}}} {varvec{I}}_{{{text{WUL}}}}$$

(17)

The set ({varvec{I}}_{{{text{SV}}}}) represents all the 4 by 1 vectors located in the Service Volume (SV) of a monitored satellite. The Worst User Location (WUL) is determined by an analytic method (Shao et al., 2011a; Zhao et al., 2014).

In (16), the true value of (sigma_{{{text{UDRE}}}}) is not known and therefore an overbound must be determined. Under the assumption that the range error ({text{d}}rho) from the clock-ephemeris satisfies the condition ({text{d}}rho sim {text{N}}left( {mu ,sigma^{2} } right)), the following formula can be obtained

$$Pleft( {left| {{text{d}}rho – mu } right| ge K_{{{text{HMI}}}} sigma } right) le P_{{{text{HMI}}}}$$

(18)

Considering the inequality

$$Pleft( {left| {{text{d}}rho – mu } right| ge K_{{{text{HMI}}}} sigma } right) ge Pleft( {left| {{text{d}}rho } right| – left| mu right| ge K_{{{text{HMI}}}} sigma } right)$$

(19)

the following inequality can be obtained

$$Pleft( {left| {{text{d}}rho } right| ge K_{{{text{HMI}}}} left( {frac{left| mu right|}{{K_{{{text{HMI}}}} }} + sigma } right)} right) le P_{{{text{HMI}}}}$$

(20)

Hence, the Gaussian distribution ({text{N}}left( {mu ,left( {frac{left| mu right|}{{K_{{{text{HMI}}}} }} + sigma } right)^{2} } right)) is used to find the bound of ({text{d}}rho) and its standard deviation is computed by

$$s = frac{{left| {{text{d}}overline{rho }} right|}}{{K_{{{text{HMI}}}} }} + sigma left( {{text{d}}rho } right)$$

(21)

Several situations are taken to tackle the integrity threats. Firstly, to meet the strict Gaussian overbounding properties required by the WAAS integrity monitoring, the CNMP algorithm including mean filter and mean error function is adopted to reduce the effects of multipath (Decleene, 2000; Shallberg et al., 2001). Secondly, the right tail Cumulative Distribution Function (CDF) is used to bound the probability of HMI, and the thresholds for the error in the corrections and the noise in the measurements can be derived (Schempp et al., 2001). More importantly, the authors also develop a model to improve the monitoring capability of monitor stations, which will be introduced below.

Groove model based on satellite-station geometry

In this subsection, a model is constructed to improve the performance of the bound of clock-ephemeris range errors. The algorithm for s is adjusted under different situations. A groove model is proposed to provide a solution.

To improve the insufficient monitoring capability of monitor stations when a satellite is moving over the boundary of the monitoring network, the geometry between the satellite and monitor stations is introduced as a kind of prior information to overcome the shortcomings of the method to determine the parameters.

The geometry of these monitor stations tracking a specific satellite is evaluated by Monitoring Geometric Dilution Of Precision (MGDOP), or further (ln (x_{MGDOP} )) which is computed by the corresponding geometry matrix (Chen et al., 2017a). The relationship between MGDOP and UDREI is analyzed in (Chen et al., 2017a; Shao et al., 2009, 2011a). Based on this, the geometric information is used to adjust the shape of the clock-ephemeris covariance to ensure the satellite integrity. Taking satellite PRN 6 as an example, the number of the monitor stations tracking PRN6 and the WAAS-reported (sigma_{{{text{UDRE}}}}) of PRN6 for one day are shown in Fig. 2. The lateral axis denotes time (t) in unit of day (d).

As shown in Fig. 2, the number of the monitor stations varies from 0 to over 30 and changes fast. Accordingly, (ln (x_{MGDOP} )) varies conversely and synchronously. The trend of (sigma_{{{text{UDRE}}}}) is the same as that of (ln (x_{MGDOP} )), and (ln (x_{MGDOP} )) can be taken as supplementary information to deduce s. Apparently, the trend of (sigma_{{{text{UDRE}}}}) or (ln (x_{MGDOP} )) for all satellites likes a groove. Therefore, the geometry of a specific satellite and the monitor stations, namely, MGDOP, can be taken as the boundary condition and further used to compensate the algorithm for deducing s, denoted as (s_{{{text{DOP}}}}). Based on these, a groove model is developed to describe the trend of s as illustrated in Fig. 3.

According to this model, the algorithms to compute s are given by

$$s_{1} = left{ begin{gathered} max left( {s,s_{{{text{limit}}}} } right),{text{N}}left( {U_{90,15} } right) ge N_{0} hfill \ max left( {s,s_{{{text{limit}}}} ,s_{{{text{DOP}}}} } right),{text{N}}left( {U_{90,15} } right) < N_{0} & {text{N}}left( {U_{90,5} } right) ge 4 hfill \ {text{NaN}},{text{otherwise}} hfill \ end{gathered} right.$$

(22)

where the parameters s and (s_{{{text{DOP}}}}) are computed by user equivalent range errors and geometry information, respectively, and NaN denotes not a number. For simplicity, let (U_{{alpha_{{text{EL(2)}}} ,alpha_{{text{EL(1)}}} }}) denote the set of UEREs with their ELevation angle (EL) (alpha_{{{text{EL}}}}) under the condition (alpha_{{text{EL(1)}}} le alpha_{{{text{EL}}}} le alpha_{{{text{EL}}(2)}}) and ({text{N}}left( {U_{{alpha_{{{text{EL}}(2)}} ,alpha_{{text{EL(1)}}} }} } right)) represent the sample size of set (U_{{alpha_{{{text{EL}}(2)}} ,alpha_{{text{EL(1)}}} }}). As described in (22) and Fig. 3, the algorithm for s1 is divided into three parts by the sample size N0 of set (U_{{alpha_{{{text{EL}}(2)}} ,alpha_{{text{EL(1)}}} }}). After the analysis of the requirement, the value of N0 is set as 8.

There are three intermediate variables to obtain the final parameter s. The parameter (s_{{{text{limit}}}}) is used to limit or bound the noise of pseudorange residuals relative to many error sources and can be computed by the standard deviation of the pseudorange residuals (or UEREs) (Blanch et al., 2012; Chen et al., 2017b, 2018). The variable s is used to evaluate pseudoranges (or UEREs) in real time, and further monitor the state of a specific satellite. UEREs are processed by CNMP algorithm and right tail CDF, and then used to deduce the variable s. The variable (s_{{{text{DOP}}}}) takes the geometry between the satellite and monitor stations as another information source to compensate for the insufficient monitoring capability of monitor stations when the satellite is moving over the boundary of the monitoring network. The geometry information can be translated into (s_{{{text{DOP}}}}) according to the relationship between them (Chen et al., 2017b; Shao et al., 2009, 2011b).

Finally, the covariance matrix for clock-ephemeris after the above processes is updated by

$${varvec{P}}_{3} = frac{{s_{1}^{2} }}{{P_{max } }}left( {frac{{K_{{{text{HMI}}}} + K_{{{text{bias}}}} + K_{{{text{pfa}}}} }}{{K_{7} }}} right)^{2} {varvec{P}}_{{{text{SCE}}}}$$

(23)

The matrix is used as the final covariance for clock-ephemeris errors which will be formatted into UDRE and MT28.

To broadcast the clock-ephemeris covariance, one needs to compress the matrix ({varvec{P}}_{3}) into UDRE and MT28. Firstly, ({varvec{P}}_{3}) is projected along the vector from the monitor station to a satellite, and the maximum of the projection is obtained by searching the service volume (Shao et al., 2011a; Zhao et al., 2014). Secondly, UDRE can be obtained by searching the lookup table of UDREI to bound the projection (SC-159, 2016). Thirdly, the matrix ({varvec{C}}) to be broadcast can be determined by finding the optimal solution with the minimum of quantization error for the clock-ephemeris covariance (Chen et al., 2018; Kailath et al., 2000). Finally, UDRE and MT28 are determined and updated within their update interval.