A Real-Time Detection Method for BDS Signal in Space Anomalies

Signal In Space (SIS) anomalies in satellite navigation systems can degrade satellite-based navigation and positioning performance. The occurrence of SIS anomalies from the BeiDou navigation satellite System (BDS) may be more frequent than for the Global Positioning System (GPS). In order to guarantee the integrity of BDS users, detecting and excluding SIS anomalies is indispensable. The traditional method through the comparison between the final precision ephemeris and the broadcast ephemeris is limited by the issue of long latency of precision ephemeris release. Through the statistical characteristics analysis of Signal In Space User Range Error (SISURE), we propose a real-time Instantaneous SISURE (IURE) estimation method by using the Kalman filtering-based carrier-smoothed-code to detect and exclude BDS SIS anomalies, in which the threshold for BDS IURE anomaly detection are obtained from the integrity requirement. The experimental results based on 1 Hz data from ground observations show that the proposed method has an estimation accuracy of 1.1 m for BDS IURE. The test results show that the proposed method can effectively detect the SIS anomalies caused by either orbit faults or clock faults.


Introduction
The BeiDou navigation satellite System (BDS) will provide Positioning, Navigation and Timing (PNT) services for the Asia-Pacific region before 2020 [1]. For most of real-time PNT users, the orbit and clock parameters are derived from broadcast ephemeris. In order to ensure the reliability of positioning, the health flag contained in the broadcast ephemeris is widely used to indicate the availability of satellites [2]. However, the broadcast ephemeris belongs to the predicted ephemeris. Due to the possible inconsistency between the prediction model and the true kinetic model, the predicted ephemeris sometimes cannot accurately describe the true satellite orbit and clock, and cause largely biased position errors. Therefore, real-time detection of satellite orbit and clock anomaly is indispensable for the reliability of positioning.
Signal In Space (SIS) errors are mainly composed of satellite orbit and clock errors, which are the major error sources of ranging observations [3][4][5]. Signal In Space User Range Error (SISURE) is the projection of SIS errors on the direction of line-of-sight, which represents the impact of SIS errors on code observations. SISURE varies with the location of users, therefore, the Instantaneous SISURE (IURE) and global average SISURE (avgURE) are widely used to describe the statistics of SIS errors and detect the SIS anomaly [2].
Post processing and the real-time estimation are two widely-used methods for SIS anomaly monitoring [6]. The post processing method obtains the SISURE by comparing the broadcast ephemeris and the final precision ephemeris [7]. Cohenour and Kovach used the post-processing method to International Terrestrial Reference Frame 2008 (ITRF2008). The satellite coordinate obtained from final precise ephemeris is based on the Center of Mass (CoM) [3]. Fortunately, the difference between CGCS2000 and ITRF2008 is at centimeter level and can be ignored when calculating SISURE [3,4,21,22]. In order to eliminate the difference between APC and CoM, the antenna Phase Center Offset (PCO) is usually corrected before comparing the broadcast ephemeris with the final precision ephemeris [2,3,21,22]. In order to coordinate the time system, the group-delay and time-reference must be considered. The clock of broadcast ephemeris is measured by using B3I signal, whereas the clock of precision ephemeris is estimated by using the dual-frequency combination of B1I and B2I signals [3,21]. In order to eliminate the bias caused by different measuring methods, the group-delay of satellites should be eliminated:  represent original clock derived from the broadcast ephemeris and clock corrected by groupdelay, respectively. Different satellite clock products refer to different time-reference. There is an unignorable bias between two different clock products [21,22]. In order to eliminate the bias, we use the weighted average method to estimate the common bias of all satellites' clocks: where ˆ() k  is the common bias between two clock products at kth epoch, ˆj t  and j t  are the satellite clock of broadcast ephemeris and final precision ephemeris for satellite j . N is the total number of all valid satellites.

IURE Computation
Based on the final precision ephemeris, the SIS errors of broadcast ephemeris can be calculated. Firstly, the time and space references of two types of ephemeris must be coordinated. For the spatial reference, the coordinate system and the center of satellite must be considered [21]. The orbital position coordinate system of BDS broadcast ephemeris is China Geodetic Coordinate System 2000 (CGCS2000). The satellite coordinate obtained from broadcast ephemeris is with respect to Antenna Phase Center (APC) [1]. The orbital position coordinate system of the final precise ephemeris is International Terrestrial Reference Frame 2008 (ITRF2008). The satellite coordinate obtained from final precise ephemeris is based on the Center of Mass (CoM) [3]. Fortunately, the difference between CGCS2000 and ITRF2008 is at centimeter level and can be ignored when calculating SISURE [3,4,21,22]. In order to eliminate the difference between APC and CoM, the antenna Phase Center Offset (PCO) is usually corrected before comparing the broadcast ephemeris with the final precision ephemeris [2,3,21,22].
In order to coordinate the time system, the group-delay and time-reference must be considered. The clock of broadcast ephemeris is measured by using B3I signal, whereas the clock of precision ephemeris is estimated by using the dual-frequency combination of B1I and B2I signals [3,21]. In order to eliminate the bias caused by different measuring methods, the group-delay of satellites should be eliminated: where f 1 , and f 2 represent the frequency of BDS B1I and B2I signals separately.
δt TGD1 and δt TGD2 denote group-delay of B1I and B2I signals, respectively. δt brd and δt represent original clock derived from the broadcast ephemeris and clock corrected by group-delay, respectively. Different satellite clock products refer to different time-reference. There is an unignorable bias between two different clock products [21,22]. In order to eliminate the bias, we use the weighted average method to estimate the common bias of all satellites' clocks: whereμ(k) is the common bias between two clock products at kth epoch, δt j and δt j are the satellite clock of broadcast ephemeris and final precision ephemeris for satellite j. N is the total number of all valid satellites. w j represents the weight of satellite j when calculating the common bias. To avoid the potential impact of the satellites clock outliers, the δt j (k) − δt j (k) with a larger residual is weighted toward less contribution by w j = 1/ δt j (k) − δt j (k) −μ(k) . Iterations are applied at each epoch until µ(k) converges to an acceptable level. The common biasμ(k) is the time-reference bias between two clock products. The clock error of broadcast ephemeris can be obtained by subtracting the time-reference bias from the difference between broadcast ephemeris and final precision ephemeris clock product.
When the time and space reference of the precise and the broadcast ephemeris are consistent, the IURE of broadcast ephemeris can be calculated by taking the final precise ephemeris as true values: where IURE j represents the IURE of satellite j,X j and X j represent the satellite position calculated by using broadcast ephemeris and final precision ephemeris, respectively. c represents the speed of light. X r represents the position of the station x = √ x · x T . In order to analyze the distribution of IURE, the outliers of IURE should be eliminated. At present, BDS has not released the relevant indicators of integrity. We assume the URA definition of BDS is similar with GPS [2,3,25]. Therefore, the outlier is eliminated with a threshold of 4.42 × URAUB, where URAUB indicates the upper boundary of URA.

IURE Model
The accurate IURE model relies on the statistical characteristics of IURE. It can be seen from (3) that the IURE changes with the coordinates of stations due to the effect of line-of-sight. Due to the limitations of the content, we take "CUT0" station as shown in Figure 1 as an example to analyze the statistical characteristics of IURE. The broadcast ephemeris comes from the "BRDM" products provided by International GNSS Service (IGS) (ftp://cddis.gsfc.nasa.gov/pub/gps/data/ campaign/mgex/daily/rinex3/). The final precision ephemeris is provided by WuHan University (WHU) (ftp://igs.ign.fr/pub/igs/products/mgex/). The time range is from February 18, 2013 to June 2, 2018, and the sampling interval is 15 minutes. The broadcast ephemeris-based satellite orbit and clock are interpolated based on BDS interface control document [1]. A total of 185472 samples are collected for the last five years. Other stations have similar results. Only three typical satellites are plotted, and the IURE within ±10 m is shown in Figure 2.
where j IURE represents the IURE of satellite j , ˆj X and j X represent the satellite position calculated by using broadcast ephemeris and final precision ephemeris, respectively. c represents the speed of light. r X represents the position of the station In order to analyze the distribution of IURE, the outliers of IURE should be eliminated. At present, BDS has not released the relevant indicators of integrity. We assume the URA definition of BDS is similar with GPS [2, 3,25]. Therefore, the outlier is eliminated with a threshold of 4.42 × URAUB, where URAUB indicates the upper boundary of URA.

IURE Model
The accurate IURE model relies on the statistical characteristics of IURE. It can be seen from (3) that the IURE changes with the coordinates of stations due to the effect of line-of-sight. Due to the limitations of the content, we take "CUT0" station as shown in Figure 1 as an example to analyze the statistical characteristics of IURE. The broadcast ephemeris comes from the "BRDM" products provided by International GNSS Service (IGS) (ftp://cddis.gsfc.nasa.gov/pub/gps/data/campaign/ mgex/daily/rinex3/). The final precision ephemeris is provided by WuHan University (WHU) (ftp://igs.ign.fr/pub/igs/products/mgex/). The time range is from February 18, 2013 to June 2, 2018, and the sampling interval is 15 minutes. The broadcast ephemeris-based satellite orbit and clock are interpolated based on BDS interface control document [1]. A total of 185472 samples are collected for the last five years. Other stations have similar results. Only three typical satellites are plotted, and the IURE within ±10 meters is shown in Figure 2.  It can be seen from the Figure 2 that the unimodal characteristics of IURE distribution are obvious for MEO and IGSO satellites, but GEO satellites have multiple peaks. In addition, the IURE It can be seen from the Figure 2 that the unimodal characteristics of IURE distribution are obvious for MEO and IGSO satellites, but GEO satellites have multiple peaks. In addition, the IURE distributions have skewed bias of −0.78 m, −1.15 m, and 0.73 m, respectively. The skewed bias distribution is caused by a systematic bias originated from satellite clock errors [3]. Therefore, the distribution of IURE does not follow zero-mean normal distribution.
Based on the time series of IURE from Figure 3, we can find that the IURE presents a certain periodicity. Therefore, we establish the IURE model as: where IURE trend and IURE stoc represent the trend term and the random term of IURE. p(1), p(2), p(3) and p(4) represent the parameters of trend term. ∆t = t−t 0 T s represents the normalized time parameters. t, t 0 and T s represent the observation time, the initial time and the sampling period, respectively. T b is the normalized period, and pi represents the circumference ratio. We make T b equal to a sidereal day. distributions have skewed bias of −0.78 m, −1.15 m, and 0.73 m, respectively. The skewed bias distribution is caused by a systematic bias originated from satellite clock errors [3]. Therefore, the distribution of IURE does not follow zero-mean normal distribution. Based on the time series of IURE from Figure 3, we can find that the IURE presents a certain periodicity. Therefore, we establish the IURE model as: T is the normalized period, and pi represents the circumference ratio. We make b T equal to a sidereal day.  As shown in Figure 3, three typical satellites are plotted and the IURE within ±10 meters are shown. We can see that each satellite has a systematic bias, which is mainly affected by the clock error [3]. In addition, it can be seen from the zoomed-in portion that the IURE of GEO satellite has obvious periodic feature. The periodic fluctuation is mainly caused by the period of satellite motion [20]. Although the trend IURE depends on the line-of-sight as well as the user location, we can see the longterm stability of trend IURE for the past five years. IUREtrend is a function of time with the coefficients p(i) (i = 1, 2, 3, 4) which can be estimated by the least square method based on the historical IURE data. As shown in Figure 3, three typical satellites are plotted and the IURE within ±10 m are shown. We can see that each satellite has a systematic bias, which is mainly affected by the clock error [3]. In addition, it can be seen from the zoomed-in portion that the IURE of GEO satellite has obvious periodic feature. The periodic fluctuation is mainly caused by the period of satellite motion [20]. Although the IURE trend depends on the line-of-sight as well as the user location, we can see the long-term stability of IURE trend for the past five years. IURE trend is a function of time with the coefficients p(i) (i = 1, 2, 3, 4) which can be estimated by the least square method based on the historical IURE data. Figure 4 plots the histogram of IURE stoc for CUT0 station. The red fitting curve can fit the distribution of IURE. The mean of IURE stoc is very close to 0, especially for GEO satellites, we can see that the multi-peak is eliminated. The distribution of IURE stoc is unimodality and symmetry, which indicates that the IURE stoc can be assumed to follow zero-mean normal distribution. Overall, the IURE can be modelled as the combination of the zero-mean normal distribution with a trend bias. can be assumed to follow zero-mean normal distribution. Overall, the IURE can be modelled as the combination of the zero-mean normal distribution with a trend bias.

IURE Anomaly Detection
Based on the constructed IURE model, we can achieve real-time estimation of IURE. Highprecision IURE estimation can be beneficial for not only accurately characterizing the ranging error, but also quickly detecting SIS anomalies with a reasonable threshold. Therefore, we propose a realtime IURE estimation method by using the Kalman filtering-based carrier-smoothed-code to improve the accuracy of IURE estimation. Furthermore, we obtain the threshold to detect SIS anomaly based on the integrity requirement.

Real-time IURE Estimation
Since the phase observation suffers from the issue of ambiguity resolution, we therefore utilize the code observation for the IURE real-time estimation. The code observation can be modelled as: where r is the actual distance between the ground station and the satellite. In order to improve the estimation performance of IURE, it is necessary to reduce the effect of various atmospheric errors. For dual-frequency signals, the Ionosphere-Free (IF) combination is used:

IURE Anomaly Detection
Based on the constructed IURE model, we can achieve real-time estimation of IURE. High-precision IURE estimation can be beneficial for not only accurately characterizing the ranging error, but also quickly detecting SIS anomalies with a reasonable threshold. Therefore, we propose a real-time IURE estimation method by using the Kalman filtering-based carrier-smoothed-code to improve the accuracy of IURE estimation. Furthermore, we obtain the threshold to detect SIS anomaly based on the integrity requirement.

Real-Time IURE Estimation
Since the phase observation suffers from the issue of ambiguity resolution, we therefore utilize the code observation for the IURE real-time estimation. The code observation can be modelled as: where r is the actual distance between the ground station and the satellite. δt u , δt brd and δt r represent the clock bias caused by receiver, satellite and relativity, respectively. T denotes tropospheric delay error. ρ i represents the code observation of carrier f i . I is the ionospheric delay on carrier f 1 . ε ρ i is the observed noise of carrier f i . In addition, the influence of the earth's rotation is corrected in the calculation of r.
In order to improve the estimation performance of IURE, it is necessary to reduce the effect of various atmospheric errors. For dual-frequency signals, the Ionosphere-Free (IF) combination is used: where ρ IF denotes IF observation. ε ρ IF is the model error and observed noise of ρ IF . In general, we utilize the traditional carrier-smoothed-code to suppress the amplified IF observation noise: where ρ(k), φ(k) andρ(k) represent the code, carrier phase and carrier-smoothed-code of the satellite at kth epoch, respectively. M denotes the traditional smoothing length.
The standard implementation of Local Area Augmentation System (LAAS) recommends the smoothing length is fixed to 100 s for GPS [26]. At present, there is no official recommendation of smoothing length for BDS [1, 25,27]. Due to the different signal systems, BDS and GPS are different in terms of code observation multipath and noise [28,29]. Therefore, the smoothing length of GPS may not be suitable for BDS. In addition, the fixed smoothing length cannot adapt to the observation quality of different satellites because BDS has three types of satellites. In order to select the optimal smoothing length, we rearrange (7) to construct the state and measuring model of Kalman filter: where W(k) and V(k) are the system process noise and observation noise, respectively. The gain K of Kalman filter is the reciprocal of M. The carrier-smoothed-code based on Kalman filter can realize the optimal control of M by using the gain K to suppress the code noise. It should be noted that it has been widely proven that raw BDS code and phase observation follow Gaussian distribution [28,29]. Therefore, we select the traditional linear Kalman filter for the carrier-smoothed-code.
We define r =r − ε e , δt = δt brd − l 1 δt TGD1 + l 2 δt TGD2 = δt − ε c , in whichr represents the range calculated by the broadcast ephemeris. ε e and ε c denote the projection of orbit and clock error on the line-of-sight, respectively. Since the IURE can be modeled as IURE = IURE trend + IURE stoc = ε e − cε c , the IURE can therefore be calculated as: whereρ IF and ερ IF are smoothed code and the corresponding noise.r can be calculated by the satellite position derived from the broadcast ephemeris and the station position obtained from IGS.ρ IF , δt, and δt r can be obtained from the receiver observation file and the broadcast ephemeris. Most of T errors can be eliminated by Saastamoinen model [30]. Therefore, we define: The IURE can be computed by combining (9) and (10): It can be found that the receiver clock error δt u is unknown at the right side of (11). For any satellite j with the same observation time, we have IURE j = IURE j trend + IURE j stoc = ρ j cal + cδt u + ε jρ IF . Since the IURE stoc follows the zero-mean normal distribution: where IURE j trend represents the trend term of satellite j that can be derived from historical data. For zero-mean observation noise, it can be assumed that 1 N N ∑ j=1 ε jρ IF = 0, therefore, the receiver clock can be estimated as: and then, the IURE can be estimated as: Based on the statistical distribution of IURE and the instantaneous observation collected at the station with precise locations, the real-time IURE estimation method can be implemented. Particularly, besides the atmospheric error correction, the receiver clock is estimated and extracted based on the zero-mean normal distribution assumption of IURE stoc . Regarding to the code observation noise, the smoothing length of carrier-smoothed-code is adjusted in real time by the constructed Kalman filter from (8), which adaptively suppresses the observation noise. Therefore, the estimation performance of IURE can be improved with deliberate observation error correction.

Detection Threshold
Threshold plays a crucial role for the anomaly detection. GPS had provided the SIS integrity definition and requirement in 2008 [31]. The threshold for the anomaly detection can be obtained as 4.42 × URAUB: where 1 × 10 −5 is the integrity requirement given by GPS [31]. At present, although BDS characterizes the accuracy of SISURE by URA contained in the broadcast ephemeris, BDS has not released the definitions of integrity. If (15) is used to determine the threshold, the IURE must follow zero-mean normal distribution. According to the previous analysis, it is found that the statistical distribution of BDS IURE are quite different from zero-mean normal distribution. Therefore, the threshold of 4.42 × URAUB is inappropriate for the anomaly detection of BDS IURE. Fortunately, the distribution of IURE stoc can be approximated as zero-mean normal distribution, so the threshold for IURE anomaly detection can be obtained as: We can take IURE trend ± 4.42 × σ IURE stoc as the threshold to detect BDS IURE anomaly, where IURE trend and σ IURE stoc can be derived from the historical data analysis.
Based on above analysis, a potential BDS SIS anomaly will be claimed when an IURE exceeds the threshold and the broadcast ephemeris is flagged as healthy. In general, the real-time detection of SIS anomaly needs to cross check the IURE anomalies of worldwide multiple stations in order to ensure the reliability of SIS anomaly detection. In this contribution, although the detection process of single station is implemented, the detection method can be the foundation for the cross check at multiple stations.
Based on the zero-mean normal distribution of IURE stoc , the threshold for BDS IURE anomaly detection is derived from the integrity requirement. It can be anticipated that the performance of SIS anomaly detection will be enhanced.

Experiments and Discussion
To verify the performance of the proposed IURE estimation method, we analyze the estimation errors of IURE based on raw code observation, denoted as the RCO method, the traditional carrier-smoothed-code, denoted as the TCSC method, and the Kalman filtering-based carrier-smoothed-code, denoted as the KCSC method. Moreover, we use the Root Mean Square (RMS) and STandard Deviation (STD) of the estimation errors to verify the accuracy and reliability of real-time IURE estimation, respectively. In order to demonstrate the edge of threshold derived from the integrity requirement, we analyze the effect of two thresholds on the IURE anomaly detection, i.e., the traditional threshold ±4.42 × URAUB, denoted as TTH, and the proposed threshold IURE trend ± 4.42 × σ IURE stoc , denoted as PTH. In addition, we take CUT0 and "NNOR" stations provided by IGS as the examples to prove the efficiency of the proposed anomaly detection method by using SPP technology.

IURE Estimation and Analysis
In order to sufficiently evaluate the estimation performance of IURE, we collected raw data from six stations as shown in Figure 1. The data sampling rate is 1 Hz and the experiment time is from 20 May 2018 to 16 June 2018, a total of 2,419,200 samples in 4 weeks. The elevation cutoff is set to 20 • . The traditional smoothing length M is set to 100s. The observation of stations comes from IGS (ftp://cddis.gsfc.nasa.ov/highrate/2018/) except CUT0 whose observation comes from Curtin GNSS Research Centre (http://saegnss2.curtin.edu.au/ldc/rinex3/daily/2018/). The reference of IURE is calculated by the post processing method. Due to the limitation of real-time data distribution, we use the offline data to simulate the real-time operation. The IURE estimation result of CUT0 station are given in Figure 5.    As observed from Figure 5, there are three parts of data missing as shown in gray area of each subpanel, in which the first and third segments are the absence of the final precise ephemeris, and the second segment is caused by BDS signal interruption. In addition, C05 has been deleted because the frequent signal interruption. Figure 5c contains a zoomed-in portion of C07 satellite. As shown in the red circle of zoomed-in portion, some newly observed satellites have large estimation errors because either low elevation or poor observation quality, and the estimation errors gradually converged when the elevation increases. Moreover, we can find that there is a systematic bias for each satellite as shown in the green circle of zoomed-in portion, which may be caused by the multipath and the unmodelled bias. It can be seen from Figure 5 that the IURE estimation accuracy of KCSC method and TCSC method are better than the RCO method, which illustrate the necessity of code noise suppression. Furthermore, the estimation errors of KCSC method are relatively less than the TSCS method. It reveals the efficiency of utilizing the Kalman filter to determine the optimal smoothing length. The higher IURE estimation accuracy will provide better performance of SIS anomaly detection.
According to the three types of BDS satellites, Tables 1 and 2 give the numerical statistics of the estimation errors for six stations. Table 1 summarizes the RMS of IURE estimation errors. As seen from Table 1, a minimum RMS of 0.74 m can be obtained from the MEO satellite at JFNG station by using the KCSC method. Moreover, we find that the GEO satellite have the largest RMS of IURE estimation errors for three methods, which can be explained by the lowest precision of final precision ephemeris for GEO satellites compared with IGSO and MEO satellites [32]. Meanwhile, it can be seen that stations at low latitudes can achieve smaller RMS than the high latitudes, which may be because stations at low latitudes can observe more satellites and achieve higher estimation accuracy of receiver clock. Regardless of the satellite types and the location of stations, it can be found that the KCSC method can achieve better estimation accuracy than the RCO method and the TCSC method, which proves the estimation accuracy of real-time IURE is benefit from the Kalman filtering-based carrier-smoothed-code algorithm. Compared with the RCO and TCSC methods, the RMS of IURE estimation errors obtained by the KCSC method is reduced by 17.3% and 8.3% respectively. The STD of IURE estimation errors is shown in Table 2, from which we find that a minimum STD of 0.39 m can be obtained from the GEO satellite at NTUS station by using the KCSC method. The MEO satellite have the largest STD compared with GEO and IGSO satellites for all methods as shown in Table 2. This mainly because the observation of MEO satellites is discontinuous, and there are more convergence processes. In addition, there is no obvious relationship between STD and station location. Compared with he RCO and TCSC methods, the STD obtained by the KCSC method is the smallest for all satellite types and all sites, which indicates that the proposed algorithm can effectively improve the reliability of real-time IURE estimation. Overall, compared with RCO and TCSC methods, the STD of IURE estimation errors obtained by the KCSC method is reduced by 38.2% and 22.5%, respectively.
The IURE estimation accuracy of 1.1 m is sufficient for satisfying the integrity requirement of 1 × 10 −5 , which is significant for applying the proposed method for the BDS-based safety-critical services.

Detection Performance of SIS Anomaly Based on the IURE Estimation
We define the SIS anomaly is caused by orbit fault when the contribution of the orbit error to the IURE is greater than the clock error. If the contribution of the clock error to the IURE is greater than the orbit error, we define the anomaly is caused by a clock fault.
The SIS anomaly detection process of C09 is shown in Figure 6. This SIS anomaly is caused by orbit fault and is detected at CUT0 station. It can be seen from the Figure 6 that the real-time IURE can track the IURE reference. Regardless of detection thresholds, the SIS can be marked as anomaly from the 20th to the 21th hour.

Detection Performance of SIS Anomaly Based on the IURE Estimation
We define the SIS anomaly is caused by orbit fault when the contribution of the orbit error to the IURE is greater than the clock error. If the contribution of the clock error to the IURE is greater than the orbit error, we define the anomaly is caused by a clock fault.
The SIS anomaly detection process of C09 is shown in Figure 6. This SIS anomaly is caused by orbit fault and is detected at CUT0 station. It can be seen from the Figure 6 that the real-time IURE can track the IURE reference. Regardless of detection thresholds, the SIS can be marked as anomaly from the 20th to the 21th hour. To further verify the efficiency of anomaly detection, we analyze the impact of C09 SIS anomaly on the SPP positioning errors. Figure 7 shows the positioning errors whether or not C09 is engaged in SPP solution from the 20th to the 21th hour. In the normal state of Figure 7c,d, there are several blue dots that exceed the red threshold, but the positioning error still meets the positioning accuracy requirement which demands that the error does not exceed 10 m under the metric of 95% probability [25]. As shown in Figure 7a,c, both the horizontal and vertical errors exceed 10 m when C09 is engaged in the position computation. When the observation of C09 is removed, the positioning error To further verify the efficiency of anomaly detection, we analyze the impact of C09 SIS anomaly on the SPP positioning errors. Figure 7 shows the positioning errors whether or not C09 is engaged in SPP solution from the 20th to the 21th hour. In the normal state of Figure 7c,d, there are several blue dots that exceed the red threshold, but the positioning error still meets the positioning accuracy requirement which demands that the error does not exceed 10 m under the metric of 95% probability [25]. As shown in Figure 7a,c, both the horizontal and vertical errors exceed 10 m when C09 is engaged in the position computation. When the observation of C09 is removed, the positioning error meets the positioning accuracy requirement as shown in Figure 7b,d. Therefore, the SIS anomaly detecting is essential to improve the accuracy and the reliability of positioning.
The SIS anomaly detection process of C11 is shown in Figure 8. The anomaly of C11 is caused by clock fault and is detected at NNOR station. Note that NNOR is located in Australia. The dual frequency code and carrier observations of NNOR are used for SIS anomaly detection. As seen from the Figure 8, the real-time IURE of C11 surpassed PTH and then exceeded TTH over time. On 27 December 2016, the IURE reference is interrupted due to the lack of the final precise ephemeris. This also implies that the real-time estimation of the IURE can effectively compensate for the absence of the final precise ephemeris for the anomaly detection.
meets the positioning accuracy requirement as shown in Figure 7b,d. Therefore, the SIS anomaly detecting is essential to improve the accuracy and the reliability of positioning. The SIS anomaly detection process of C11 is shown in Figure 8. The anomaly of C11 is caused by clock fault and is detected at NNOR station. Note that NNOR is located in Australia. The dual frequency code and carrier observations of NNOR are used for SIS anomaly detection. As seen from the Figure 8, the real-time IURE of C11 surpassed PTH and then exceeded TTH over time. On December 27, 2016, the IURE reference is interrupted due to the lack of the final precise ephemeris. This also implies that the real-time estimation of the IURE can effectively compensate for the absence of the final precise ephemeris for the anomaly detection.
To further verify the efficiency of anomaly detection, we analyze the impact of C11 SIS anomaly on the SPP positioning errors. Figure 9 shows the positioning errors whether or not C11 is engaged in SPP solution.
As shown in Figure 9a,c, when the IURE is greater than PTH but less than TTH in the normal/anomaly state, the SIS anomaly of C11 causes vertical positioning errors exceeding the position accuracy requirement. If C11 is removed, the positioning accuracy will be constrained to an acceptable level, as shown in Figure 9b  The invisible state indicates that the receiver has not captured the satellite signal or the elevation is lower than 20°. The state of "lack of precision ephemeris" denotes that the final precise ephemeris is missing. The normal/anomaly state represents that the real-time IURE is greater than PTH but less than TTH when the broadcast ephemeris is marked as healthy. The anomaly state indicates that the real-time IURE is larger than PTH when the broadcast ephemeris is marked as healthy. RIURE and PIURE represent the real-time IURE estimation and the reference of IURE respectively. UTC denotes coordinated universal time. The invisible state indicates that the receiver has not captured the satellite signal or the elevation is lower than 20 • . The state of "lack of precision ephemeris" denotes that the final precise ephemeris is missing. The normal/anomaly state represents that the real-time IURE is greater than PTH but less than TTH when the broadcast ephemeris is marked as healthy. The anomaly state indicates that the real-time IURE is larger than PTH when the broadcast ephemeris is marked as healthy. RIURE and PIURE represent the real-time IURE estimation and the reference of IURE respectively. UTC denotes coordinated universal time.
To further verify the efficiency of anomaly detection, we analyze the impact of C11 SIS anomaly on the SPP positioning errors. Figure 9 shows the positioning errors whether or not C11 is engaged in SPP solution.
December, 2016. The invisible state indicates that the receiver has not captured the satellite signal or the elevation is lower than 20°. The state of "lack of precision ephemeris" denotes that the final precise ephemeris is missing. The normal/anomaly state represents that the real-time IURE is greater than PTH but less than TTH when the broadcast ephemeris is marked as healthy. The anomaly state indicates that the real-time IURE is larger than PTH when the broadcast ephemeris is marked as healthy. RIURE and PIURE represent the real-time IURE estimation and the reference of IURE respectively. UTC denotes coordinated universal time.
(a) (c) (b) (d) Figure 9. SPP errors of NNOR station. The blue dots and red lines indicate the positioning errors and the threshold of 10 m. The normal state indicates that the positioning error is consistent with the position accuracy requirement, whereas the anomaly state means the position accuracy requirement cannot be satisfied. The normal/anomaly state represents that the real-time IURE is greater than PTH but less than TTH when the broadcast ephemeris is marked as healthy. (a,c) are the horizontal and vertical positioning errors with C11 when the SIS is marked as anomaly, respectively. (b,d) denote the horizontal and vertical positioning errors without C11 when the SIS is marked as anomaly, respectively. UTC denotes coordinated universal time. Figure 9. SPP errors of NNOR station. The blue dots and red lines indicate the positioning errors and the threshold of 10 m. The normal state indicates that the positioning error is consistent with the position accuracy requirement, whereas the anomaly state means the position accuracy requirement cannot be satisfied. The normal/anomaly state represents that the real-time IURE is greater than PTH but less than TTH when the broadcast ephemeris is marked as healthy. (a,c) are the horizontal and vertical positioning errors with C11 when the SIS is marked as anomaly, respectively. (b,d) denote the horizontal and vertical positioning errors without C11 when the SIS is marked as anomaly, respectively. UTC denotes coordinated universal time.

Conclusions
As shown in Figure 9a,c, when the IURE is greater than PTH but less than TTH in the normal/anomaly state, the SIS anomaly of C11 causes vertical positioning errors exceeding the position accuracy requirement. If C11 is removed, the positioning accuracy will be constrained to an acceptable level, as shown in Figure 9b,d. It implies that the proposed threshold IURE trend ± 4.42 × σ IURE stoc can better guarantee the accuracy and reliability of positioning than the traditional threshold ±4.42 × URAUB, which proves the edge of the proposed threshold based on integrity requirement.

Conclusions
Through the post-processing of the BDS data collected from 2013 to 2018, we have explored the statistical distribution of IURE. It is found that the PCO should be corrected with different correction products at different stages. Furthermore, it has been shown the IURE does not follow zero-mean normal distribution, but the IURE can be modelled as the combination of the zero-mean normal distribution with a trend bias. With the atmospheric error correction and receiver clock estimation, we develop the IURE estimation method by using the Kalman filtering-based carrier-smoothed-code algorithm in order to suppress the observation noise adaptively. Meanwhile, we obtain the threshold derived from the integrity requirement of 1 × 10 −5 to detect IURE anomaly. Based on the experiment of six stations, the IURE estimation results has revealed that the RMS and the STD of estimation errors can achieve 1.1 m and 0.55 m, respectively. Compared with the RCO and TCSC methods, the proposed IURE estimation method can obtain the best IURE estimation accuracy regardless of satellite types and locations. Compared with the RCO and TCSC methods, the proposed method can improve the IURE estimate accuracy by 17.3% and 8.3% respectively, the STD of IURE estimation errors obtained by the proposed method is reduced by 38.2% and 22.5% respectively. The test results of CUT0 and NNOR stations show that the proposed detection threshold can provide better accuracy and reliability of user positioning than the traditional detection threshold. Meanwhile, it has been demonstrated that the proposed SIS anomaly detection method can effectively identify the SIS anomalies caused by either orbit faults or clock faults based on the SPP result analysis. Note that the SIS anomalies detected in this contribution are obtained based on the observation of single station. More reliable SIS anomaly detection requires cross-checking among worldwide multiple stations, which will be our future work topic.