# A Solution for Real-Time Ionospheric Delay Using an Adaptive Kalman Filter Based on Estimating the Variance Component

Authors: Xu Yang; Qianxin Wang; Guobin Chang
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Real-time solution of Global Navigation Satellite System (GNSS) epoch-differenced ionospheric delay (DID) is of great significance for real-time cycle slip detection and repair of multi-GNSS dual-frequency or trifrequency undifferenced measurements under high ionospheric activity. We construct a dynamic model of DID and perform a real-time estimate of the noise level of DID based on estimating the variance component. The estimated and predicted values of DID are obtained by designing a new adaptive Kalman filter algorithm with colored noise. Combining the predicted value and the detection method for cycle slips for Melbourne-Wübbena (MW) and Geometry-Free (GF) combination and taking into account the correlation between the predicted value and the carrier signal, we estimate the cycle slip, N2, on the second frequency of the carrier signal. The prediction and estimate of DID and detection and repair of dual-frequency cycle slip of multisystem undifferenced phase observations are measured with the GNSS multisystem observational data at different sampling rates (30 s, 15 s, 10 s, and 1 s). The results show that the DID model constructed in this paper is correct. The predicted value of DID has a high accuracy, which can effectively assist in dual-frequency cycle slip detection and repair. (1) The obtained predicted values, the estimated value, and the difference value between the two values of DID are less than 1.2 cm (STD), 1.2 cm (STD), and 0.6 cm (STD), respectively; (2) the precisions of the detection of cycle slip for MW, GF, and N2 are less than 0.083 cycles (STD), 0.4 cm (STD), and 0.071 cycles (STD), respectively; (3) with the obtained predicted value of DID to aid the detection and repair of cycle slip in GNSS double-frequency signals, a success rate of 100% can be reached.


---

## Body

## 1. Introduction

Real-time multifrequency and multisystem positioning have been widely used because there are more and more GNSS satellite navigation systems and system frequencies, stricter real-time performance of data and products, and increasingly perfect mathematical models of fusion technology. To achieve high-precision real-time positioning performance, real-time detection and repair of cycle slip are two of the difficult problems that must be solved. In view of this, our paper starts with the calculation of the real-time DID to assist in real-time detection and repair of cycle slip in GNSS double-frequency undifferenced phase observation. This provides a useful reference for real-time cycle detection and repair of multifrequency and multisystem positioning under high ionospheric activity [1–6].In GNSS data processing, ionospheric delay is one of the most difficult errors to deal with [7]. The research on global or local high-precision simulation and correction methods is its main research content [8]. This paper will study the DID solution based on single station, slant path, and real-time. The following three aspects are considered.The first aspect is the modeling, estimation, and prediction of DID. Based on the short-term stability of the ionospheric delay, especially for the high frequency GNSS observations, the time difference of carrier phase observations can eliminate part of the ionospheric delay error [9–11]. However, due to the complexity of the ionosphere, a remaining ionospheric delay error is still inevitable after the time difference is handled [12–14]. By considering the time-varying characteristics of ionospheric delay [15–17] and the lack of an actual theoretical model and adopting the least squares variance component estimation (LS-VCE) method to adjust the noise level of DID [18, 19], we conduct online, adaptive modeling to obtain the estimated and predicted values of DID and analyze the magnitude and variation characteristics of DID.The second aspect is that the stochastic model of GNSS carrier phase time difference observation needs further description. In this paper, carrier phase time difference observation equations are used to calculate the DID. However, the existing stochastic model of the observations has not been correctly described and therefore has not been fully considered [20–22]. On the one hand, the time difference process obviously introduces a temporal correlation of the observation sequence [23–25]. In other words, the measurement noise of observations is colored noise. Further analysis finds that the correlation involved in this paper only exists in the observations of neighboring epochs. To be more specific, the noise sequence can be defined as {wt}. Then, we have cov⁡wt,wt-1≠0, cov⁡wt,wt-k=0, k≥2. This colored noise is a type of noise with a special nature. It is a first-order moving average sequence. That is, it is a finite impulse response (FIR). Obviously this finite impulse response cannot be strictly described by the commonly used first-order autoregressive model because the latter is an infinite impulse response [26–28]. On the other hand, at the statistical level, the real-time determination of DID is actually a recursive least squares estimate. If we focus on a single epoch, the single epoch actually degenerates into a least squares estimation (weighted) problem [29]. When considering the temporal correlation of observables, it is theoretically unreasonable to perform data processing on only a single epoch. The theoretical method, which is applicable to this problem, should be “recursive least squares method under colored noise.” However, the existing literature has no recursive least squares method that has been used to deal with the colored noise of finite impulse response [30]. In view of this, we will research and structure a practical random model that can correctly reflect the colored noise characteristics of the carrier phase time differential observations. The strict recursive least squares method with this model will be used to solve the DID. At the same time, considering the Kalman filter algorithm being a fast algorithm based on the recursive least squares [31] and online adjustment of DID time-varying characteristics, a new adaptive Kalman filter algorithm is needed to be designed.The third aspect is the correlation between predicted values of DID and epoch-differenced carrier phases. This correlation may not be ignored. The reason why this correlation exists is that both variables in the current epoch are affected by the phase error in the previous epoch. This correlation should be considered when we construct a stochastic model to estimate the narrow cycle slip, N2, by weighted least squares estimators.At present there are TurboEdit, Kalman filtering, MW, GF, and polynomial fitting methods for cycle detection of GNSS double-frequency undifferenced phase observations. In addition, there are many improvement and combination methods based on the traditional methods [32–34]. There are two main problems with these methods: (1) most of the methods are based on postprocessing or need to be further improved for real-time processing (such as polynomial fitting). (2) DID changes are large when the ionosphere delay varies drastically or when a magnetic storm occurs [3, 4]. This condition will lead to the fact that some methods do not work as well because of the lack of sufficient consideration of DID (such as the GF method). In view of the above shortcomings, after the high-precision solution of the DID by considering the above three aspects problems, high-precision cycle slip detection processes of MW, GF, and N2 are constructed by using the predicted values of DID and are all applied to the detection and repair of cycle slip. Corresponding experiments with the measured data at various sampling intervals are carried out to verify the high precision of DID values and the superiority of the cycle slip detection and repair method proposed in this paper.


## 2. Ionospheric Delay Prediction and Cycle Slip Repair with an Adaptive Kalman Filter Based on Estimating the Variance Component

### 2.1. Prediction of DID

The dynamic model of DID is constructed according to [35, 36](1)X˙=Δu˙Δu¨Δu⃛=010001000ΔuΔu˙Δu¨+00Δu⃛=AX+wThe covariance matrix of the state vectorX is given by (2)Q=cov⁡w=σΔu⃛2000000001In (1), the vector X includes the DID and its first and second derivative with respect to time. Its third derivative with respect to time is considered to be completely random noise. Assume the sampling interval is t; (1) is made discrete as (3), (4), and (5):(3)Xt=ΦXt-1+w-t-1with(4)Φ=exp⁡AT=I3+AT+12A2T2=1TT2201T001(5)Q-t-1=cov⁡w-t-1=cov⁡∫0Texp⁡AtwT-tdt=σΔu⃛2T520T48T36T48T33T22T36T22T=σΔu⃛2ΓIt is noted that, in (4), An=0, whenn ≥ 3. Based on (3), the prediction will be performed before the repair of the cycle slip of each epoch by using (6)X-t=ΦX^t-1(7)P-t=ΦP^t-1AT+σΔu⃛2ΓThe first element in prediction equation (6), the predicted of DID, is used to correct the pseudo-range and carrier observations of the current epoch, and this correction is used to repair the second and third cycle slip detection which will be discussed later. Meanwhile, the predicted variance at the top left of (7), the first element of P-t, will be used in the constructed stochastic model for the estimate of the narrow lane cycle slip, N2. After the detection and repair of the cycle slip the repaired carrier observations can be used to estimate the DID.


### 2.2. Estimation of DID

By combining together all measurement errors that are not related to ionospheric delay, the pseudo-range and phase measurement equations between the epochs are given by(8)Δpjt=Δρt+kjΔut+Δεjt(9)Δφjt=λjΔϕjt=Δρt-kjΔut+λjΔNjt+ΔξjtIn the above two equations, Δ denotes the epoch-differential operator.j = 1 or 2 denotes thej-th frequency.t denotes the epoch.f denotes the frequency. kj=f12/fj2, denotes the ionospheric delay coefficient relative to the first frequency. Δη denotes the variation of ionospheric delay on the first frequency. ΔN denotes the cycle slip (assuming the repair of the cycle slip, so it is a known value).λdenotes the wavelength. Δε and Δξrepresent the noise of the pseudo-range and carrier between the epochs, and its values are simply considered as std[Δεj]=2σp and std[Δξj]=2σφ, respectively, here assuming all these errors are independent of each other. Although the pseudo-range in (8) can be used to repair the MW detection, it is not used to estimate the DID due to its low accuracy. An observation equation is constructed to estimate the DID as (10)Δφ=φ1t-φ1t-1-λ1N∨1φ2t-φ2t-1-λ2N∨2=1-k11-k2ΔρΔu+Δξ1Δξ2=BΔρΔu+ΔξIn (10), N∨j denotes the rounded cycle slip. A least squares estimate is given by (11)Δρ~Δu~=BTB-1BTΔφSince the DID is only item we need to consider, its estimation equations are given by(12)Δu~=bTΔφ(13)σΔu~2=2bTbσφ2In the above two equations, vectorb is a known quantity. The estimate in (12) happens to be a measurement update of the subsequent update of the Kalman filter. It is worth noting that here the observation sequence is polluted by FIR colored noise, which will give full consideration to the design of the subsequent Kalman filter algorithm.The observations in (3) are redefined as (14)Δu~t=100Xt+dΔu~=mTXt+dΔu~In the above equation the operator “d” denotes “error in” or “uncertainty of”. Here the update to the Kalman filter measurement is converted to a least squares problem. A new pseudo-measurement is obtained by combining the predicted value in (6) and the observation in (14). Then the new observations are given by (15) and (16):(15)lt=X-tΔu~t=I3mTXt+dX-tdΔu~t=I3mTXt+ΦdX^t-1+w-t-1dΔu~t=MXt+Jwith(16)QJJ=P-tΦgt-1gt-1TΦTσΔu~2=ΦP^t-1ΦTΦgt-1gt-1TΦTσΔu~2+σΔu⃛2Γ000=Q-+σΔu⃛2ΨIn (16), gt-1=cov⁡dX^t-1,dΔu~t is the updated covariance during the measurement update phase of the previous epoch. Starting from least squares theory, there is an estimate given by (17)X^t=MTQJJ-1M-1MTQJJ-1lt=Mlt(18)P^t=MTQJJ-1M-1At the same time there is a covariance given by (19), which will be used in the measurement update phase of the next epoch.(19)
g
t
=
cov
⁡
d
X
^
t
,
d
Δ
u
~
t
+
1
=
cov
⁡
A
d
l
t
,
d
Δ
u
~
t
+
1
=
cov
⁡
A
0
3
×
1
d
Δ
u
~
t
,
d
Δ
u
~
t
+
1
=
cov
⁡
A
0
3
×
1
b
T
d
Δ
φ
t
,
b
T
d
Δ
ξ
t
+
1
=
cov
⁡
A
0
3
×
1
b
T
ξ
t
,
-
b
T
ξ
t
=
-
A
0
3
×
1
b
T
b
σ
φ
2For the first epoch (t = 1), it is assumed that there is no cycle slip or the cycle slip is repaired by other methods. At the same time there is no available prior information of DID, and there is only an estimated value in (14). The other state components, the first and the second derivatives of the DID, are unknown for the current epoch and hence are completely irrelevant for the DID of the next epoch. The equation for determining the covariance of the first epoch should be replaced by (20)g1=cov⁡dX^1,dΔu~2=cov⁡bTΔξ1,bTΔξ200T=cov⁡bTξ1,-bTξ100T=σφ2bTb00TAt this point, the correlation considerations in (16) and the correlation updates in (19) and (20) represent the special processing method of colored measurement noise in this paper.


### 2.3. Solution of the Noise Level of DID Based on Estimating the Variance Component

As shown above, there is another problem that has not been solved, namely, the adjustment of uncertain parameters (the noise level of DID) of the dynamic model. The adjustment can be solved using a new statistical method [37], and in this paper it is treated as a problem that estimates the variance components. Although there are many methods of estimating variance components without losing the generality, here we use the convenient and widely used least square method [38].Equation (21) is the singular value decomposition of the design matrix A in (15).(21)M=U∗ΣV=U∗Σ-3×30VThe result of the linear transformation of observation equation (15) is given by (22)st=Ult=Σ-0VXt+UJThe last element on the left of the above equation represents the information in addition to the estimated parameters. This additional information can be used to estimate the required solution of the variance component. Therefore, we haves4t=u4TJ and (23)s42≈var⁡s4t=u4TQ-u4+σΔu⃛2u4TΨu4=σ-2+σΔu⃛2r2whereu4T denotes the fourth row of theU vector. Finally, this variance component is simply calculated by (24)σ^Δu⃛2=s42-σ-2r2This is actually a degenerate version of the least squares estimate of the variance component [38], since weights are not necessary for a pure observation and parameters as shown in (23). It is noted here that there is only one degree of freedom for estimating a component, so the estimate is simple and straightforward. If r is not zero, the variance component can be estimated [39]. During the estimation process the variance component that is determined may be a negative value. The reasons for the negative variance are as follows: (1) there are not enough observations in the function model. (2) The stochastic model is improperly selected. (3) When the variance (covariance) component is solved, ana priori value is incorrect. To solve this problem when negative numbers appear, the component is simply set to zero [40]. Then the corresponding formula is defined as (25)σ^Δu⃛2=max⁡s42-σ-2r2,0Now there are two sets of information related to the level of process noise, that is, knowna prioriinformation and the currently calculated information in (24). We naturally want to combine them in a statistical way. However, this is difficult. Although the uncertainty or the estimate of variance in (24) is easily provided in least squares theory [38], the uncertainty of prior information is not so easy to determine. Instead of combining the noise in statistical significance (26) is used to combine them.(26)σΔu⃛2←1-hσΔu⃛2+hσ^Δu⃛2In the above equation,h is called the learning rate; for example, it can be 0.1. The learning rate balances the contributions of historical information and current information. In other words, what we are doing here is using forgotten memory [41] to reflect the active extent of ionospheric delay of the current epoch. For some of the epochs that have just started, such as the first ten epochs, with the relative weight of historical information equal to 0.9, which is too large, the update strategy in (26) is not used, instead we use (27)σΔu⃛2←t-1tσΔu⃛2+1tσ^Δu⃛2Ifr=0 in an epoch, it means that the variance component is not estimated, and the update in (26) or (27) is not performed. The variance value of the previous epoch is simply passed to this epoch. The update variance in (26) and (27) will be used in the time and measurement updates and in the update variances of the ionospheric delay of the next epoch. Finally, it is noted that (1) estimating the variance component and updating the noise variance will not be performed for the first epoch (t= 1) and (2) if the variance component has an anomaly (such as a serious deviation from the empirical value), the variance component is directly set to an empirical value.


### 2.4. Cycle Slip Repair

The cycle slip repair method proposed in this paper is performed in a step-by-step manner to handle real-time applications. Before the repair of cycle slip, the MW combined observation [32–34] is selected as the detection of the first cycle slip, the GF combined observation [32–34] repaired by the predicted value of DID is selected as the detection of the second cycle slip, and the cycle slip in the second frequency, N2, is selected as the detection of the third cycle. The first and third detection are used to detect and repair cycle slips. The second detection is used to check the cycle slip completeness (see later) on the first and third detection and to detect cycle slip as well. In a specific operation of the cycle slip repair we generally need to perform a step-by-step solution. First, we determine the first detection which is a wide-lane detection and is calculated by direct rounding rather than using a complex decorrelation or search algorithm. Second, we determine the second detection which is a narrow lane detection. Then we determine the third detection. Finally, the detection values of dual-frequency cycle slips can be solved immediately by combining the first and the third detection. Then, the cycle slips of carrier phase measurements can be repaired with the detection values.The following is a brief description of the repair of the third detection and the cycle slip completeness checking.(  1) Repair of the Third Detection. The fixed solution of the detection of the first cycle slip can be obtained by simply rounding the detection. We believe that the success rate of wide-lane combination rounding can reach 99.99% after using GF detection to check for cycle slip completeness for the first detection. The relationship among the cycle slip in the first frequency, the first detection, and the third detection is given by [3, 4, 32–34](28)ΔN1=ΔNw+ΔN2The third detection will be obtained by using the carrier signals on L1 and L2 in a statistically optimal manner. Since the pseudo-ranges have low precision and their contribution is very small, they are omitted here to save calculations and make the algorithm more suitable for real-time operations. The observation equation for this estimation problem is given by(29)Δφ=Δφ1-λ1ΔNwΔφ2=1λ11λ2ΔρΔN2-k1k2dΔu-+Δξ1Δξ2=CΔρΔN2-kΔu-+ΔξCombining the prediction variance in the upper left in (7), namely, σdΔu-2=var⁡dΔu-, we have the covariance given by (30)f=cov⁡Δξ,dΔu-=cov⁡Δξt,mTΦdX^t-1=cov⁡Δξt,mTΦMςt-1=cov⁡Δξt,mTΦM03×1dΔu~t-1=cov⁡Δξt,mTΦM03×1bTdΔφt-1=cov⁡-ξt-1,mTΦM03×1bTξt-1=-02×3σφ2bMTΦTmFinally, we have the covariance matrix for observation equation (29) given by (31)QΔφ=cov⁡-kdΔu-+Δξ=σdΔu-2kkT+2σϕ2I2-kfT-fkTSo we have the least squares estimate as(32)Δρ^ΔN^2=CTQΔφ-1C-1CTQΔφ-1ΔφWe hope that the fixed solution can be obtained by simply rounding the float solution after using GF detection to check for cycle slip completeness for the third detection, that is, obtaining the second element of (29) with a relatively satisfactory success rate. Finally, after plugging the first and the third rounded cycle slip detection into (28), the cycle slip on the first frequency can be solved.(  2) Cycle Slip Completeness Checking. The predicted value of DID can reach a high precision by using the method proposed above. The predicted value is used to repair the second detection. At the same time the second detection is also used to check cycle slip completeness on the first and third detection, which is shown in (33). For example, when the float solution of MW detection is between 0.5 and 1.5, the detection value is set to one after the direct rounding. If the GF detection is less than 0.05, then the MW detection rounding error may be determined first. And the detection rounding is zero again. Then the N2 detection is estimated, and a judgment similar to MW detection for N2 is performed. The cycle slip in the two frequencies is repaired. The estimated value of DID is calculated and the variance of DID is updated. If the ratio of variance at the current epoch to the previous epoch is small (generally less than ten; if the cycle slip is incorrectly repaired, the ratio would become a large value), then the detection rounded to zero is correct; otherwise the detection is still rounded to one. (33)GF<0.05MW=MW-0.5,MW∈0.5,1.0MW=MW-1,MW∈1.0,1.5MW=MW+0.5,MW∈-1.0,-0.5MW=MW+1,MW∈-1.5,-1,GF<0.05N2=N2-0.5,N2∈0.5,1.0N2=N2-1,N2∈1.0,1.5N2=N2+0.5,N2∈-1.0,-1.5N2=N2+1,N2∈-1.5,-1The algorithm design flow of the above-mentioned cycle slip detection and repair method is shown in Figure1, and the execution process is as follows: (1) initialize the filter parameters at the first epoch; (2) perform Kalman filter time update to obtain the predicted value and variance of DID; (3) use the pseudo-range and carrier observation data between the epochs, the predicted value, and its variance to construct the GF, MW, and N2 detection, then perform the cycle slip completeness checking, and perform the detection and repair of cycle slip; (4) Kalman filter measurement update: estimate the DID with the carrier signals repaired by the cycle slips, perform the Kalman filter measurement update, and update the variance of DID; (5) judge the variation of the DID variance: if the ratio of the DID variance of the current epoch to the DID variance of the previous epoch is larger than the threshold, namely, ΔσΔu⃛2≥σm2, the repair is considered to be unsuccessful; then perform the detection and repair from step (3) again, and do not check for cycle slip completeness; (6) return to step (2) to process the next epoch; otherwise terminate the process.Figure 1
Flow chart of the algorithm design for real-time detection and repair of cycle slip in GNSS double-frequency undifferenced phase observation with the predicted value of DID.


## 2.1. Prediction of DID

The dynamic model of DID is constructed according to [35, 36](1)X˙=Δu˙Δu¨Δu⃛=010001000ΔuΔu˙Δu¨+00Δu⃛=AX+wThe covariance matrix of the state vectorX is given by (2)Q=cov⁡w=σΔu⃛2000000001In (1), the vector X includes the DID and its first and second derivative with respect to time. Its third derivative with respect to time is considered to be completely random noise. Assume the sampling interval is t; (1) is made discrete as (3), (4), and (5):(3)Xt=ΦXt-1+w-t-1with(4)Φ=exp⁡AT=I3+AT+12A2T2=1TT2201T001(5)Q-t-1=cov⁡w-t-1=cov⁡∫0Texp⁡AtwT-tdt=σΔu⃛2T520T48T36T48T33T22T36T22T=σΔu⃛2ΓIt is noted that, in (4), An=0, whenn ≥ 3. Based on (3), the prediction will be performed before the repair of the cycle slip of each epoch by using (6)X-t=ΦX^t-1(7)P-t=ΦP^t-1AT+σΔu⃛2ΓThe first element in prediction equation (6), the predicted of DID, is used to correct the pseudo-range and carrier observations of the current epoch, and this correction is used to repair the second and third cycle slip detection which will be discussed later. Meanwhile, the predicted variance at the top left of (7), the first element of P-t, will be used in the constructed stochastic model for the estimate of the narrow lane cycle slip, N2. After the detection and repair of the cycle slip the repaired carrier observations can be used to estimate the DID.


## 2.2. Estimation of DID

By combining together all measurement errors that are not related to ionospheric delay, the pseudo-range and phase measurement equations between the epochs are given by(8)Δpjt=Δρt+kjΔut+Δεjt(9)Δφjt=λjΔϕjt=Δρt-kjΔut+λjΔNjt+ΔξjtIn the above two equations, Δ denotes the epoch-differential operator.j = 1 or 2 denotes thej-th frequency.t denotes the epoch.f denotes the frequency. kj=f12/fj2, denotes the ionospheric delay coefficient relative to the first frequency. Δη denotes the variation of ionospheric delay on the first frequency. ΔN denotes the cycle slip (assuming the repair of the cycle slip, so it is a known value).λdenotes the wavelength. Δε and Δξrepresent the noise of the pseudo-range and carrier between the epochs, and its values are simply considered as std[Δεj]=2σp and std[Δξj]=2σφ, respectively, here assuming all these errors are independent of each other. Although the pseudo-range in (8) can be used to repair the MW detection, it is not used to estimate the DID due to its low accuracy. An observation equation is constructed to estimate the DID as (10)Δφ=φ1t-φ1t-1-λ1N∨1φ2t-φ2t-1-λ2N∨2=1-k11-k2ΔρΔu+Δξ1Δξ2=BΔρΔu+ΔξIn (10), N∨j denotes the rounded cycle slip. A least squares estimate is given by (11)Δρ~Δu~=BTB-1BTΔφSince the DID is only item we need to consider, its estimation equations are given by(12)Δu~=bTΔφ(13)σΔu~2=2bTbσφ2In the above two equations, vectorb is a known quantity. The estimate in (12) happens to be a measurement update of the subsequent update of the Kalman filter. It is worth noting that here the observation sequence is polluted by FIR colored noise, which will give full consideration to the design of the subsequent Kalman filter algorithm.The observations in (3) are redefined as (14)Δu~t=100Xt+dΔu~=mTXt+dΔu~In the above equation the operator “d” denotes “error in” or “uncertainty of”. Here the update to the Kalman filter measurement is converted to a least squares problem. A new pseudo-measurement is obtained by combining the predicted value in (6) and the observation in (14). Then the new observations are given by (15) and (16):(15)lt=X-tΔu~t=I3mTXt+dX-tdΔu~t=I3mTXt+ΦdX^t-1+w-t-1dΔu~t=MXt+Jwith(16)QJJ=P-tΦgt-1gt-1TΦTσΔu~2=ΦP^t-1ΦTΦgt-1gt-1TΦTσΔu~2+σΔu⃛2Γ000=Q-+σΔu⃛2ΨIn (16), gt-1=cov⁡dX^t-1,dΔu~t is the updated covariance during the measurement update phase of the previous epoch. Starting from least squares theory, there is an estimate given by (17)X^t=MTQJJ-1M-1MTQJJ-1lt=Mlt(18)P^t=MTQJJ-1M-1At the same time there is a covariance given by (19), which will be used in the measurement update phase of the next epoch.(19)
g
t
=
cov
⁡
d
X
^
t
,
d
Δ
u
~
t
+
1
=
cov
⁡
A
d
l
t
,
d
Δ
u
~
t
+
1
=
cov
⁡
A
0
3
×
1
d
Δ
u
~
t
,
d
Δ
u
~
t
+
1
=
cov
⁡
A
0
3
×
1
b
T
d
Δ
φ
t
,
b
T
d
Δ
ξ
t
+
1
=
cov
⁡
A
0
3
×
1
b
T
ξ
t
,
-
b
T
ξ
t
=
-
A
0
3
×
1
b
T
b
σ
φ
2For the first epoch (t = 1), it is assumed that there is no cycle slip or the cycle slip is repaired by other methods. At the same time there is no available prior information of DID, and there is only an estimated value in (14). The other state components, the first and the second derivatives of the DID, are unknown for the current epoch and hence are completely irrelevant for the DID of the next epoch. The equation for determining the covariance of the first epoch should be replaced by (20)g1=cov⁡dX^1,dΔu~2=cov⁡bTΔξ1,bTΔξ200T=cov⁡bTξ1,-bTξ100T=σφ2bTb00TAt this point, the correlation considerations in (16) and the correlation updates in (19) and (20) represent the special processing method of colored measurement noise in this paper.


## 2.3. Solution of the Noise Level of DID Based on Estimating the Variance Component

As shown above, there is another problem that has not been solved, namely, the adjustment of uncertain parameters (the noise level of DID) of the dynamic model. The adjustment can be solved using a new statistical method [37], and in this paper it is treated as a problem that estimates the variance components. Although there are many methods of estimating variance components without losing the generality, here we use the convenient and widely used least square method [38].Equation (21) is the singular value decomposition of the design matrix A in (15).(21)M=U∗ΣV=U∗Σ-3×30VThe result of the linear transformation of observation equation (15) is given by (22)st=Ult=Σ-0VXt+UJThe last element on the left of the above equation represents the information in addition to the estimated parameters. This additional information can be used to estimate the required solution of the variance component. Therefore, we haves4t=u4TJ and (23)s42≈var⁡s4t=u4TQ-u4+σΔu⃛2u4TΨu4=σ-2+σΔu⃛2r2whereu4T denotes the fourth row of theU vector. Finally, this variance component is simply calculated by (24)σ^Δu⃛2=s42-σ-2r2This is actually a degenerate version of the least squares estimate of the variance component [38], since weights are not necessary for a pure observation and parameters as shown in (23). It is noted here that there is only one degree of freedom for estimating a component, so the estimate is simple and straightforward. If r is not zero, the variance component can be estimated [39]. During the estimation process the variance component that is determined may be a negative value. The reasons for the negative variance are as follows: (1) there are not enough observations in the function model. (2) The stochastic model is improperly selected. (3) When the variance (covariance) component is solved, ana priori value is incorrect. To solve this problem when negative numbers appear, the component is simply set to zero [40]. Then the corresponding formula is defined as (25)σ^Δu⃛2=max⁡s42-σ-2r2,0Now there are two sets of information related to the level of process noise, that is, knowna prioriinformation and the currently calculated information in (24). We naturally want to combine them in a statistical way. However, this is difficult. Although the uncertainty or the estimate of variance in (24) is easily provided in least squares theory [38], the uncertainty of prior information is not so easy to determine. Instead of combining the noise in statistical significance (26) is used to combine them.(26)σΔu⃛2←1-hσΔu⃛2+hσ^Δu⃛2In the above equation,h is called the learning rate; for example, it can be 0.1. The learning rate balances the contributions of historical information and current information. In other words, what we are doing here is using forgotten memory [41] to reflect the active extent of ionospheric delay of the current epoch. For some of the epochs that have just started, such as the first ten epochs, with the relative weight of historical information equal to 0.9, which is too large, the update strategy in (26) is not used, instead we use (27)σΔu⃛2←t-1tσΔu⃛2+1tσ^Δu⃛2Ifr=0 in an epoch, it means that the variance component is not estimated, and the update in (26) or (27) is not performed. The variance value of the previous epoch is simply passed to this epoch. The update variance in (26) and (27) will be used in the time and measurement updates and in the update variances of the ionospheric delay of the next epoch. Finally, it is noted that (1) estimating the variance component and updating the noise variance will not be performed for the first epoch (t= 1) and (2) if the variance component has an anomaly (such as a serious deviation from the empirical value), the variance component is directly set to an empirical value.


## 2.4. Cycle Slip Repair

The cycle slip repair method proposed in this paper is performed in a step-by-step manner to handle real-time applications. Before the repair of cycle slip, the MW combined observation [32–34] is selected as the detection of the first cycle slip, the GF combined observation [32–34] repaired by the predicted value of DID is selected as the detection of the second cycle slip, and the cycle slip in the second frequency, N2, is selected as the detection of the third cycle. The first and third detection are used to detect and repair cycle slips. The second detection is used to check the cycle slip completeness (see later) on the first and third detection and to detect cycle slip as well. In a specific operation of the cycle slip repair we generally need to perform a step-by-step solution. First, we determine the first detection which is a wide-lane detection and is calculated by direct rounding rather than using a complex decorrelation or search algorithm. Second, we determine the second detection which is a narrow lane detection. Then we determine the third detection. Finally, the detection values of dual-frequency cycle slips can be solved immediately by combining the first and the third detection. Then, the cycle slips of carrier phase measurements can be repaired with the detection values.The following is a brief description of the repair of the third detection and the cycle slip completeness checking.(  1) Repair of the Third Detection. The fixed solution of the detection of the first cycle slip can be obtained by simply rounding the detection. We believe that the success rate of wide-lane combination rounding can reach 99.99% after using GF detection to check for cycle slip completeness for the first detection. The relationship among the cycle slip in the first frequency, the first detection, and the third detection is given by [3, 4, 32–34](28)ΔN1=ΔNw+ΔN2The third detection will be obtained by using the carrier signals on L1 and L2 in a statistically optimal manner. Since the pseudo-ranges have low precision and their contribution is very small, they are omitted here to save calculations and make the algorithm more suitable for real-time operations. The observation equation for this estimation problem is given by(29)Δφ=Δφ1-λ1ΔNwΔφ2=1λ11λ2ΔρΔN2-k1k2dΔu-+Δξ1Δξ2=CΔρΔN2-kΔu-+ΔξCombining the prediction variance in the upper left in (7), namely, σdΔu-2=var⁡dΔu-, we have the covariance given by (30)f=cov⁡Δξ,dΔu-=cov⁡Δξt,mTΦdX^t-1=cov⁡Δξt,mTΦMςt-1=cov⁡Δξt,mTΦM03×1dΔu~t-1=cov⁡Δξt,mTΦM03×1bTdΔφt-1=cov⁡-ξt-1,mTΦM03×1bTξt-1=-02×3σφ2bMTΦTmFinally, we have the covariance matrix for observation equation (29) given by (31)QΔφ=cov⁡-kdΔu-+Δξ=σdΔu-2kkT+2σϕ2I2-kfT-fkTSo we have the least squares estimate as(32)Δρ^ΔN^2=CTQΔφ-1C-1CTQΔφ-1ΔφWe hope that the fixed solution can be obtained by simply rounding the float solution after using GF detection to check for cycle slip completeness for the third detection, that is, obtaining the second element of (29) with a relatively satisfactory success rate. Finally, after plugging the first and the third rounded cycle slip detection into (28), the cycle slip on the first frequency can be solved.(  2) Cycle Slip Completeness Checking. The predicted value of DID can reach a high precision by using the method proposed above. The predicted value is used to repair the second detection. At the same time the second detection is also used to check cycle slip completeness on the first and third detection, which is shown in (33). For example, when the float solution of MW detection is between 0.5 and 1.5, the detection value is set to one after the direct rounding. If the GF detection is less than 0.05, then the MW detection rounding error may be determined first. And the detection rounding is zero again. Then the N2 detection is estimated, and a judgment similar to MW detection for N2 is performed. The cycle slip in the two frequencies is repaired. The estimated value of DID is calculated and the variance of DID is updated. If the ratio of variance at the current epoch to the previous epoch is small (generally less than ten; if the cycle slip is incorrectly repaired, the ratio would become a large value), then the detection rounded to zero is correct; otherwise the detection is still rounded to one. (33)GF<0.05MW=MW-0.5,MW∈0.5,1.0MW=MW-1,MW∈1.0,1.5MW=MW+0.5,MW∈-1.0,-0.5MW=MW+1,MW∈-1.5,-1,GF<0.05N2=N2-0.5,N2∈0.5,1.0N2=N2-1,N2∈1.0,1.5N2=N2+0.5,N2∈-1.0,-1.5N2=N2+1,N2∈-1.5,-1The algorithm design flow of the above-mentioned cycle slip detection and repair method is shown in Figure1, and the execution process is as follows: (1) initialize the filter parameters at the first epoch; (2) perform Kalman filter time update to obtain the predicted value and variance of DID; (3) use the pseudo-range and carrier observation data between the epochs, the predicted value, and its variance to construct the GF, MW, and N2 detection, then perform the cycle slip completeness checking, and perform the detection and repair of cycle slip; (4) Kalman filter measurement update: estimate the DID with the carrier signals repaired by the cycle slips, perform the Kalman filter measurement update, and update the variance of DID; (5) judge the variation of the DID variance: if the ratio of the DID variance of the current epoch to the DID variance of the previous epoch is larger than the threshold, namely, ΔσΔu⃛2≥σm2, the repair is considered to be unsuccessful; then perform the detection and repair from step (3) again, and do not check for cycle slip completeness; (6) return to step (2) to process the next epoch; otherwise terminate the process.Figure 1
Flow chart of the algorithm design for real-time detection and repair of cycle slip in GNSS double-frequency undifferenced phase observation with the predicted value of DID.


## 3. Experimental Analysis

The example data for this paper is the dual-frequency observation data of the HKKT station from the Hong Kong CORS (Continuously Operating Reference Stations) network. The sampling interval is 1 sec, and the sampling duration is 4 h (0:00 to 4:00 AM, February 1, 2017). The observational data with sampling intervals of 10, 15, and 30 sec are extracted from the original observations to analyze the experimental effect under different sampling rates. The carriers phase (B1 and B2) BDS observational data of C01, C02, C03, C04, C05, C06, C07, C08, C09, and C13 and the carriers phase (L1 and L2) GPS observational data of G2, G5, G13, G15, G20, G21, G24, and G29 are used to analyze the experimental data. At the same time the dual-frequency observational data (0:00 to 5:00 AM, March 19, 2017, sampling interval being 30 s) of IGS (International GPS Service) station JFNG are used. Before the experiment the cycle slip detection with TEQC software was performed on the original observation data, and the detection results showed that there was no cycle slip. In the course of the experiment a given amount of cycle slips was artificially added to the original observational data. The three cycle slip detection processes (MW, GF, and N2, with cycle slip determination thresholds taken as 3 cycles, 0.05 m, 1 cycle, respectively) were used to detect and repair cycle slip for observational data at different sampling rates and different types of satellites to test the effectiveness of the new method.


### 3.1. Detection and Repair Results without Cycle Slips

(  1) Variation Results of DID. We performed detection and repair of the cycles for the raw observation data of the HKKT station. The statistical results of average value, maximum value, and STD were obtained for the following: the estimated values, predicted values, and difference between estimated and predicted values of DID of C01, G02, and the other 16 satellites under four different sampling intervals. The statistics of each satellite are shown in Figure 2. The statistics of estimated and predicted values for the 30 s sampling interval are shown in Figure 3. The statistics of all the satellites are shown in Table 1.Table 1
Statistical results of estimated and predicted values and difference between estimated and predicted values of DID for all satellites with four sampling intervals.
Sampling interval (s) Mean Maximum STD Difference values (m) Estimated values (m) Predicted values (m) Difference values (m) Estimated values (m) Predicted values (m) Difference values (m) Estimated values (m) Predicted values (m) 30 0.005 0.015 0.015 0.026 0.054 0.049 0.006 0.012 0.012  15 0.004 0.008 0.008 0.026 0.035 0.027 0.005 0.006 0.006 10 0.003 0.006 0.005 0.022 0.027 0.019 0.004 0.005 0.004 1 0.002 0.002 0.001 0.015 0.015 0.004 0.002 0.002 0.000 Note: the mean, maximum, and STD values in this table are the average values of the corresponding statistics for all of the satellites in Figure2.Figure 2
Statistical results of estimated values (a), predicted values (b), and difference between estimated and predicted values (c) of DID for each satellite under four sampling intervals.
(a) (b) (c)Figure 3
Variations of estimated and predicted values of DID for GPS and BDS satellites under 30-s sampling intervals.(a) As can be seen from Table1, the STDs of the predicted and the estimated values of DID do not exceed 0.012 m, and the STD of the difference between estimated and predicted values does not exceed 0.006 m. This indicates that the solutions of the predicted and the estimated values have a high accuracy, and the two values fit very well. At the same time the smaller the sampling interval, the higher the precision of the two values.(b) As can be seen from Figure3, the predicted values can well reflect the trend in the change of ionospheric delay, and they have a certain degree of smoothing effect, while the estimated values show a certain degree of fluctuation. The main reason is that the predicted values are calculated from a stable prediction model, and the estimated values are calculated from the interfered noise measurement model.(c) As can be seen from Figures2 and 3, the degree of variation for the DID varies from satellite to satellite. The variation of most GPS and BDS/MEO satellites is larger than BDS/IGSO and BDS/GEO satellites, and some GEO satellites have the smallest change. The maximum values of the estimated and predicted values of the DID can reach 0.17 m. If this predicted value is not used to assist the detection of other cycle slips, the accuracy of the detection of the cycle slip for N2 and GF will decrease dramatically, leading to the failure of cycle slip repair.(  2) Statistical Results of Cycle Slip Detection.The statistical results are obtained for average value, maximum value, and STD of the MW, GF, and N2 detection values of C01, G02, and the other 16 satellites with four different sampling intervals. The statistics of each satellite are shown in Figure 4. The statistics of all the satellites are shown in Table 2.Table 2
Statistical results of MW, GF, and N2 detections for all of the satellites with four sampling intervals.
Sampling interval Mean Maximum STD Number of error Success rate GF N2 MW GF N2 MW GF N2 MW (s) (m) (cycle) (cycle) (m) (cycle) (cycle) (m) (cycle) (cycle)  ( % ) 30 0.003 0.054 0.065 0.017 0.288 0.295 0.004 0.071 0.083 0 100% 15 0.002 0.042 0.038 0.017 0.271 0.209 0.003 0.057 0.050 0 100% 10 0.002 0.038 0.028 0.015 0.253 0.171 0.003 0.051 0.037 0 100% 1 0.001 0.023 0.009 0.010 0.175 0.063 0.002 0.029 0.012 0 100% Note: the mean, maximum, and STD values in this table are the average values of the corresponding statistics for all of the satellites in Figure4.Figure 4
Statistical results of MW (a), GF (b), and N2 (c) detection for each satellite with four sampling intervals.
(a) (b) (c)(a) As can be seen from Table2 the error detection number is zero at various sampling intervals; that is, the cycle slips of all satellites are correctly detected and repaired.(b) As can be seen from Figure4, the maximum values of MW, GF, and N2 detection are 0.648, 0.045, and 0.459, respectively, which are less than their set thresholds of 3, 0.05, and 1, indicating that the thresholds are reasonable. From Table 2, The STDs of the three detection are 0.083, 0.004, and 0.071, respectively, indicating that the three detection processes of cycle slip have a high accuracy.(c) As can be seen from Table2, the smaller the sampling interval, the higher the accuracy of the three detection processes of cycle slip. The main reason for this is that the smaller the sampling interval, the smaller the ionospheric delay variation and the more accurate the predicted ionospheric value.(d) As can be seen from Figure4, the accuracy of detection of cycle slips for GEO satellites is the highest. The main reason is that the GEO satellites have the smallest speeds, and the values of DID of the satellites are small. The reason why the accuracy of detection of cycle slips for BDS satellites is slightly better than that of GPS satellites is that the DID of the most BDS satellites is less than the GPS satellites during the observation period.


### 3.2. Effective Validation of Cycle Slip Detection and Repair

Figure5 shows the detection and repair results of G10 of the GFNG station without using the predicted values of DID to correct the carrier observation measurement in the estimation process of N2 detection. Without using the predicted values to correct the carrier measurement, especially under the conditions that the ionospheric delay variation is relatively large, it can be seen that the repair failure of the N2 detection appears first. Then the estimated and predicted values of DID are incorrectly calculated, which then cause the GF detection to deviate from the real situation. Finally, these effects lead to the detection and repair failure of subsequent epochs.Figure 5
Results of cycle slip detection and repair without using predicted values to correct carrier observations in the estimation process of N2 detection.Figure6 shows the detection and repair results of G10 of the GFNG station without using the predicted values of DID to repair the GF detection. It can be seen that the GF detection gradually increases with the constant increase in the ionospheric delay. When the value of GF exceeds its threshold, the cycle slips are not correctly detected, and the N2 detection is also affected. Then the estimated and predicted values of DID are incorrectly calculated. Finally, these effects lead to the error results of detection and repair of some epochs.Figure 6
Results of cycle slip detection and repair without using predicted values to repair GF detection.Figure7 shows the detection and repair results of G10 of the GFNG station without checking for completeness of the cycle slips. For some epochs with large pseudo-range noise, the MW detection values are more than 0.5 cycles; the one cycle deviation of MW detection is introduced because of simple rounding of the detection thus causing an incorrect solution of the N2 detection. Then, the estimated and predicted values of DID are incorrectly calculated, which in turn causes the GF detection to deviate from the real situation. Finally, these effects lead to the error results of detection and repair of some epochs.Figure 7
Results of cycle slip detection and repair without using GF detection to check for cycle slip completeness.The smoothed process in Figure8(a) refers to the multiple epochs smoothing of the detection of the MW cycle slip without checking for cycle slip completeness (here assuming no cycle slip). Figure 8(c) illustrates the real variations of the three detection processes after the unsmoothed process with the detection and repair method proposed in this paper (here the condition of cycle slip is unknown). The related results of ionospheric delay in the unsmoothed process are exactly the same as the results of Figures 8(b) and 8(d), and therefore they are not repeated. From Figure 8(a) it can be seen that the magnitude of MW detection gradually approaches a zero cycle after the smooth processing, and therefore the effect of the MW detection on the accurate determination of the DID and other detection of cycle slip can be neglected. From Figure 8(c) it can be seen that the magnitude of MW detection fluctuates within 1.5 cycles. From Figure 8 it can be seen that the results of cycle slip detection and repair using the smooth and unsmoothed MW cycle detection are exactly the same except for the variation of MW detection. This fully shows the validity and efficiency of checking for cycle slip completeness. It can also be clearly seen from Figure 8 that the variance of DID obtained by the method in this paper can accurately reflect the degree of change of the DID, that is, when the change is fast, the variance is large. Meanwhile, the predicted values of DID have a high accuracy, which ensures the detection of the GF and N2 cycle slip has a high accuracy and the cycle slip repairs are valid.Figure 8
Results of cycle slip detection and repair with the smooth and unsmoothed MW cycle detection. (a) and (c) show the variations of the GF, MW, and N2 detection with the smoothed process and without, respectively. (b) and (d) show the variations of estimated and predicted values of DID and the variance of DID with the smoothed process, respectively.
(a) (b) (c) (d)From the above analyses and figures, it can be determined that (1) the calculation of the predicted value of DID is interleaved with repair of the cycle slip, and they complement each other; that is, the predicted values can be used to repair the carrier phase measurement and GF and N2 detection, while the repaired cycle slips can be used to accurately calculate the predicted values; (2) checking for the cycle slip completeness can effectively avoid the one cycle deviation caused by the direct rounding of the detection of the MW and N2 cycle slips; (3) the predicted value is very important for estimating the N2 detection and checking for the cycle slip completeness; (4) the predicted value has better stability than the estimated value. It will not be easily affected by the failure of the cycle slip detection and repair and seriously deviates from the real situation. Meanwhile, the failure of the current epoch detection and repair does not have a significant impact on the detection and repair of subsequent epochs, since the method proposed in this paper is for real-time single epoch processing. That is, although the current epoch detection and repair is wrong, the next epoch may be right. The above analysis fully demonstrates that the constructed model of the DID has high accuracy and is particularly suitable for real-time processing.


### 3.3. Adding Cycle Slips Experiment

At first the known cycle slips are artificially added to the carrier phase measurement of G05 and C13 satellites of the HKKT station with four sampling intervals. Then the detection and repair of cycle slips with the method mentioned in this paper are carried out. The results of cycle slip detection and repair are shown in Tables3 and 4. The variations of MW, GF, and N2 detection are shown in Figures 9 and 12. The variations of estimated and predicted values of DID are shown in Figures 10 and 13. The variations of noise level of DID are shown in Figures 11 and 14.Table 3
Results of slip detection and repair for G05 satellite with the 30-s sampling interval.
L1 cycle ship L2 cycle ship Epoch GF detection N2 detection MW detection L1 repair L2 repair GF float value N2 float value MW float value (cycle) (cycle) (cycle) (cycle) (m) (cycle) (cycle) 1 1 59 Y Y N 1 1 -0.0543 1.0073 -0.0843 1 1 60 Y Y N 1 1 -0.0556 1.0307 -0.0507 -1 -1 109 Y Y N -1 -1 0.0591 -1.0969 -0.2524 -1 -1 110 Y Y N -1 -1 0.0568 -1.0530 0.0503 0 1 159 Y Y N 0 1 -0.2505 1.1175 -1.3559 2 0 160 Y N N 2 0 0.3789 0.0312 1.9861 4 3 180 N Y N 4 3 0.0406 2.7756 0.8113 -1 0 209 Y N N -1 0 -0.1826 -0.1423 -1.2365 0 -2 210 Y Y N 0 -2 0.4949 -2.1195 1.8622 -5 -4 220 N Y N -5 -4 0.0228 -3.9528 -1.0002 763 590 230 Y Y Y 763 590 1.1051 590.0917 173.1338 -763 -590 240 Y Y Y -763 -590 -1.1085 -590.0282 -172.9307 9 7 260 N Y N 9 7 0.0010 7.0397 2.3608 2 1 285 Y Y N 2 1 0.1280 1.1560 0.9846 -9 -7 310 N Y N -9 -7 -0.0086 -6.8993 -1.9787 3 3 319 Y Y N 3 3 -0.1634 3.0315 0.0866 3 -3 320 Y Y Y 3 -3 1.3074 -3.0721 5.9982 1 1 342 N Y N 1 1 -0.0615 1.1408 0.1181 -1 -2 352 Y Y N -1 -2 0.2899 -1.8476 0.5916 2 2 360 Y Y N 2 2 -0.1113 2.0646 -0.2162Table 4
Results of slip detection and repair for C13 satellite with the 30-s sampling interval.
B1 cycle ship B2 cycle ship Epoch GF detection N2 detection MW detection B1 repair B2 repair GF float value N2 float value MW float value (cycle) (cycle) (cycle) (cycle) (m) (cycle) (cycle) 1 1 59 Y Y N 1 1 -0.0564 1.0014 -0.0047 1 1 60 Y Y N 1 1 -0.0529 0.9388 -0.0398 -1 -1 109 Y Y N -1 -1 0.0565 -1.0034 0.0196 -1 -1 110 Y Y N -1 -1 0.0556 -0.9875 -0.0647 0 1 159 Y Y N 0 1 -0.2550 1.1177 -1.1477 2 0 160 Y N N 2 0 0.3753 0.1553 1.8488 4 3 180 N Y N 4 3 0.0207 3.0434 0.8433 -1 0 209 Y N N -1 0 -0.1887 -0.0601 -0.8586 0 -2 210 Y Y N 0 -2 0.4916 -1.9088 2.1022 -5 -4 220 N Y N -5 -4 0.0319 -3.9774 -0.9753 763 590 230 N Y Y 763 590 0.0038 589.9320 173.0592 -763 -590 240 N Y Y -763 -590 0.0010 -590.0178 -173.0423 9 7 260 Y Y N 9 7 -0.0093 6.9855 2.0975 2 1 285 Y Y N 2 1 0.1346 1.0204 0.8239 -9 -7 310 N Y N -9 -7 0.0131 -7.0537 -1.8256 3 3 319 Y Y N 3 3 -0.1683 2.9889 -0.1218 3 -3 320 Y Y Y 3 -3 1.3225 -3.0239 5.9734 1 1 342 Y Y N 1 1 -0.0594 1.0550 -0.0164 -1 -2 352 Y Y N -1 -2 0.3040 -1.9888 0.8781 2 2 360 Y Y N 2 2 -0.1089 1.9333 0.0064Figure 9
Variations of MW, GF, and N2 detection for G05 satellite after the same cycle slips being artificially added to the carrier phase measurement at same time point with four sampling intervals (30s, 15s, 10s, and 1s). (Note: the conditions of big cycle slips detection and repair are not shown in Figure9.)Figure 10
Variations of estimated and predicted values of DID for the G05 satellite with four sampling intervals.Figure 11
Variations of noise level of DID for the G05 satellite with four sampling intervals.Figure 12
Variations of MW, GF, and N2 detection for C13 satellite after the same cycle slips being artificially added to the carrier phase measurement at same time point with four sampling intervals (30s, 15s, 10s, and 1s). (Note: the conditions of big cycle slips detection and repair are not shown in Figure12.)Figure 13
Variations of estimated and predicted values of DID for the C13 satellite with four sampling intervals.Figure 14
Variations of noise level of DID for the C13 satellite with four sampling intervals.(1) From Figures9 and 12, it can be seen that since the same cycle slip is artificially added to the carrier phase measurement at same time point with different sampling intervals, the variations of the three detection processes are the same after performing the correct detection and repair of the cycle slip. That is, the cycle slip can be correctly detected and repaired by using the method proposed in this paper at different sampling intervals. From Tables 3 and 4 and Figures 9 and 12, it can be seen that most of the cycle slip combinations can be detected with the GF detection, but the combination of wavelength ratio reciprocal (such as 763:590 or 9:7) cannot be detected. The combinations of large cycle slip can be preferably detected with the MW detection, but the small cycle slips (less than one cycle) and the same cycle slips cannot be detected. All the cycle slips except the zero cycle in the second frequency can be detected with the N2 detection. All the small, large, continuous, and insensitive cycle slips can be correctly detected and repaired.(2) From Figures10 and 13 and Table 1, it can be seen that the ionospheric delay under a relatively large sampling interval varies greatly, and the estimated and predicted values of DID exceed 0.15 m. The smaller the sampling interval is, the smaller the estimated and predicted values of the ionospheric delay are. For the 1 s sampling interval the predicted values are close to zero and the STD of estimated values is not more than 0.002, which indicate that the influence of ionospheric delay can be ignored in cycle slip detection and repair with a small sampling interval. The estimated and the predicted values of the ionospheric delay of a GPS satellite (G05) are larger than the BDS satellite (C13).(3) From Figures11 and 14 it can be seen that the variance of the DID changes, which preferably reflects the variation degree of the ionospheric delay. That is, where the ionosphere delay varies greatly, the noise level is larger. Meanwhile, the smaller the sampling interval is, the greater the variance of the DID is. The variance is also different for different satellites. Meanwhile, the variance is also different for different sampling intervals and satellites.


## 3.1. Detection and Repair Results without Cycle Slips

(  1) Variation Results of DID. We performed detection and repair of the cycles for the raw observation data of the HKKT station. The statistical results of average value, maximum value, and STD were obtained for the following: the estimated values, predicted values, and difference between estimated and predicted values of DID of C01, G02, and the other 16 satellites under four different sampling intervals. The statistics of each satellite are shown in Figure 2. The statistics of estimated and predicted values for the 30 s sampling interval are shown in Figure 3. The statistics of all the satellites are shown in Table 1.Table 1
Statistical results of estimated and predicted values and difference between estimated and predicted values of DID for all satellites with four sampling intervals.
Sampling interval (s) Mean Maximum STD Difference values (m) Estimated values (m) Predicted values (m) Difference values (m) Estimated values (m) Predicted values (m) Difference values (m) Estimated values (m) Predicted values (m) 30 0.005 0.015 0.015 0.026 0.054 0.049 0.006 0.012 0.012  15 0.004 0.008 0.008 0.026 0.035 0.027 0.005 0.006 0.006 10 0.003 0.006 0.005 0.022 0.027 0.019 0.004 0.005 0.004 1 0.002 0.002 0.001 0.015 0.015 0.004 0.002 0.002 0.000 Note: the mean, maximum, and STD values in this table are the average values of the corresponding statistics for all of the satellites in Figure2.Figure 2
Statistical results of estimated values (a), predicted values (b), and difference between estimated and predicted values (c) of DID for each satellite under four sampling intervals.
(a) (b) (c)Figure 3
Variations of estimated and predicted values of DID for GPS and BDS satellites under 30-s sampling intervals.(a) As can be seen from Table1, the STDs of the predicted and the estimated values of DID do not exceed 0.012 m, and the STD of the difference between estimated and predicted values does not exceed 0.006 m. This indicates that the solutions of the predicted and the estimated values have a high accuracy, and the two values fit very well. At the same time the smaller the sampling interval, the higher the precision of the two values.(b) As can be seen from Figure3, the predicted values can well reflect the trend in the change of ionospheric delay, and they have a certain degree of smoothing effect, while the estimated values show a certain degree of fluctuation. The main reason is that the predicted values are calculated from a stable prediction model, and the estimated values are calculated from the interfered noise measurement model.(c) As can be seen from Figures2 and 3, the degree of variation for the DID varies from satellite to satellite. The variation of most GPS and BDS/MEO satellites is larger than BDS/IGSO and BDS/GEO satellites, and some GEO satellites have the smallest change. The maximum values of the estimated and predicted values of the DID can reach 0.17 m. If this predicted value is not used to assist the detection of other cycle slips, the accuracy of the detection of the cycle slip for N2 and GF will decrease dramatically, leading to the failure of cycle slip repair.(  2) Statistical Results of Cycle Slip Detection.The statistical results are obtained for average value, maximum value, and STD of the MW, GF, and N2 detection values of C01, G02, and the other 16 satellites with four different sampling intervals. The statistics of each satellite are shown in Figure 4. The statistics of all the satellites are shown in Table 2.Table 2
Statistical results of MW, GF, and N2 detections for all of the satellites with four sampling intervals.
Sampling interval Mean Maximum STD Number of error Success rate GF N2 MW GF N2 MW GF N2 MW (s) (m) (cycle) (cycle) (m) (cycle) (cycle) (m) (cycle) (cycle)  ( % ) 30 0.003 0.054 0.065 0.017 0.288 0.295 0.004 0.071 0.083 0 100% 15 0.002 0.042 0.038 0.017 0.271 0.209 0.003 0.057 0.050 0 100% 10 0.002 0.038 0.028 0.015 0.253 0.171 0.003 0.051 0.037 0 100% 1 0.001 0.023 0.009 0.010 0.175 0.063 0.002 0.029 0.012 0 100% Note: the mean, maximum, and STD values in this table are the average values of the corresponding statistics for all of the satellites in Figure4.Figure 4
Statistical results of MW (a), GF (b), and N2 (c) detection for each satellite with four sampling intervals.
(a) (b) (c)(a) As can be seen from Table2 the error detection number is zero at various sampling intervals; that is, the cycle slips of all satellites are correctly detected and repaired.(b) As can be seen from Figure4, the maximum values of MW, GF, and N2 detection are 0.648, 0.045, and 0.459, respectively, which are less than their set thresholds of 3, 0.05, and 1, indicating that the thresholds are reasonable. From Table 2, The STDs of the three detection are 0.083, 0.004, and 0.071, respectively, indicating that the three detection processes of cycle slip have a high accuracy.(c) As can be seen from Table2, the smaller the sampling interval, the higher the accuracy of the three detection processes of cycle slip. The main reason for this is that the smaller the sampling interval, the smaller the ionospheric delay variation and the more accurate the predicted ionospheric value.(d) As can be seen from Figure4, the accuracy of detection of cycle slips for GEO satellites is the highest. The main reason is that the GEO satellites have the smallest speeds, and the values of DID of the satellites are small. The reason why the accuracy of detection of cycle slips for BDS satellites is slightly better than that of GPS satellites is that the DID of the most BDS satellites is less than the GPS satellites during the observation period.


## 3.2. Effective Validation of Cycle Slip Detection and Repair

Figure5 shows the detection and repair results of G10 of the GFNG station without using the predicted values of DID to correct the carrier observation measurement in the estimation process of N2 detection. Without using the predicted values to correct the carrier measurement, especially under the conditions that the ionospheric delay variation is relatively large, it can be seen that the repair failure of the N2 detection appears first. Then the estimated and predicted values of DID are incorrectly calculated, which then cause the GF detection to deviate from the real situation. Finally, these effects lead to the detection and repair failure of subsequent epochs.Figure 5
Results of cycle slip detection and repair without using predicted values to correct carrier observations in the estimation process of N2 detection.Figure6 shows the detection and repair results of G10 of the GFNG station without using the predicted values of DID to repair the GF detection. It can be seen that the GF detection gradually increases with the constant increase in the ionospheric delay. When the value of GF exceeds its threshold, the cycle slips are not correctly detected, and the N2 detection is also affected. Then the estimated and predicted values of DID are incorrectly calculated. Finally, these effects lead to the error results of detection and repair of some epochs.Figure 6
Results of cycle slip detection and repair without using predicted values to repair GF detection.Figure7 shows the detection and repair results of G10 of the GFNG station without checking for completeness of the cycle slips. For some epochs with large pseudo-range noise, the MW detection values are more than 0.5 cycles; the one cycle deviation of MW detection is introduced because of simple rounding of the detection thus causing an incorrect solution of the N2 detection. Then, the estimated and predicted values of DID are incorrectly calculated, which in turn causes the GF detection to deviate from the real situation. Finally, these effects lead to the error results of detection and repair of some epochs.Figure 7
Results of cycle slip detection and repair without using GF detection to check for cycle slip completeness.The smoothed process in Figure8(a) refers to the multiple epochs smoothing of the detection of the MW cycle slip without checking for cycle slip completeness (here assuming no cycle slip). Figure 8(c) illustrates the real variations of the three detection processes after the unsmoothed process with the detection and repair method proposed in this paper (here the condition of cycle slip is unknown). The related results of ionospheric delay in the unsmoothed process are exactly the same as the results of Figures 8(b) and 8(d), and therefore they are not repeated. From Figure 8(a) it can be seen that the magnitude of MW detection gradually approaches a zero cycle after the smooth processing, and therefore the effect of the MW detection on the accurate determination of the DID and other detection of cycle slip can be neglected. From Figure 8(c) it can be seen that the magnitude of MW detection fluctuates within 1.5 cycles. From Figure 8 it can be seen that the results of cycle slip detection and repair using the smooth and unsmoothed MW cycle detection are exactly the same except for the variation of MW detection. This fully shows the validity and efficiency of checking for cycle slip completeness. It can also be clearly seen from Figure 8 that the variance of DID obtained by the method in this paper can accurately reflect the degree of change of the DID, that is, when the change is fast, the variance is large. Meanwhile, the predicted values of DID have a high accuracy, which ensures the detection of the GF and N2 cycle slip has a high accuracy and the cycle slip repairs are valid.Figure 8
Results of cycle slip detection and repair with the smooth and unsmoothed MW cycle detection. (a) and (c) show the variations of the GF, MW, and N2 detection with the smoothed process and without, respectively. (b) and (d) show the variations of estimated and predicted values of DID and the variance of DID with the smoothed process, respectively.
(a) (b) (c) (d)From the above analyses and figures, it can be determined that (1) the calculation of the predicted value of DID is interleaved with repair of the cycle slip, and they complement each other; that is, the predicted values can be used to repair the carrier phase measurement and GF and N2 detection, while the repaired cycle slips can be used to accurately calculate the predicted values; (2) checking for the cycle slip completeness can effectively avoid the one cycle deviation caused by the direct rounding of the detection of the MW and N2 cycle slips; (3) the predicted value is very important for estimating the N2 detection and checking for the cycle slip completeness; (4) the predicted value has better stability than the estimated value. It will not be easily affected by the failure of the cycle slip detection and repair and seriously deviates from the real situation. Meanwhile, the failure of the current epoch detection and repair does not have a significant impact on the detection and repair of subsequent epochs, since the method proposed in this paper is for real-time single epoch processing. That is, although the current epoch detection and repair is wrong, the next epoch may be right. The above analysis fully demonstrates that the constructed model of the DID has high accuracy and is particularly suitable for real-time processing.


## 3.3. Adding Cycle Slips Experiment

At first the known cycle slips are artificially added to the carrier phase measurement of G05 and C13 satellites of the HKKT station with four sampling intervals. Then the detection and repair of cycle slips with the method mentioned in this paper are carried out. The results of cycle slip detection and repair are shown in Tables3 and 4. The variations of MW, GF, and N2 detection are shown in Figures 9 and 12. The variations of estimated and predicted values of DID are shown in Figures 10 and 13. The variations of noise level of DID are shown in Figures 11 and 14.Table 3
Results of slip detection and repair for G05 satellite with the 30-s sampling interval.
L1 cycle ship L2 cycle ship Epoch GF detection N2 detection MW detection L1 repair L2 repair GF float value N2 float value MW float value (cycle) (cycle) (cycle) (cycle) (m) (cycle) (cycle) 1 1 59 Y Y N 1 1 -0.0543 1.0073 -0.0843 1 1 60 Y Y N 1 1 -0.0556 1.0307 -0.0507 -1 -1 109 Y Y N -1 -1 0.0591 -1.0969 -0.2524 -1 -1 110 Y Y N -1 -1 0.0568 -1.0530 0.0503 0 1 159 Y Y N 0 1 -0.2505 1.1175 -1.3559 2 0 160 Y N N 2 0 0.3789 0.0312 1.9861 4 3 180 N Y N 4 3 0.0406 2.7756 0.8113 -1 0 209 Y N N -1 0 -0.1826 -0.1423 -1.2365 0 -2 210 Y Y N 0 -2 0.4949 -2.1195 1.8622 -5 -4 220 N Y N -5 -4 0.0228 -3.9528 -1.0002 763 590 230 Y Y Y 763 590 1.1051 590.0917 173.1338 -763 -590 240 Y Y Y -763 -590 -1.1085 -590.0282 -172.9307 9 7 260 N Y N 9 7 0.0010 7.0397 2.3608 2 1 285 Y Y N 2 1 0.1280 1.1560 0.9846 -9 -7 310 N Y N -9 -7 -0.0086 -6.8993 -1.9787 3 3 319 Y Y N 3 3 -0.1634 3.0315 0.0866 3 -3 320 Y Y Y 3 -3 1.3074 -3.0721 5.9982 1 1 342 N Y N 1 1 -0.0615 1.1408 0.1181 -1 -2 352 Y Y N -1 -2 0.2899 -1.8476 0.5916 2 2 360 Y Y N 2 2 -0.1113 2.0646 -0.2162Table 4
Results of slip detection and repair for C13 satellite with the 30-s sampling interval.
B1 cycle ship B2 cycle ship Epoch GF detection N2 detection MW detection B1 repair B2 repair GF float value N2 float value MW float value (cycle) (cycle) (cycle) (cycle) (m) (cycle) (cycle) 1 1 59 Y Y N 1 1 -0.0564 1.0014 -0.0047 1 1 60 Y Y N 1 1 -0.0529 0.9388 -0.0398 -1 -1 109 Y Y N -1 -1 0.0565 -1.0034 0.0196 -1 -1 110 Y Y N -1 -1 0.0556 -0.9875 -0.0647 0 1 159 Y Y N 0 1 -0.2550 1.1177 -1.1477 2 0 160 Y N N 2 0 0.3753 0.1553 1.8488 4 3 180 N Y N 4 3 0.0207 3.0434 0.8433 -1 0 209 Y N N -1 0 -0.1887 -0.0601 -0.8586 0 -2 210 Y Y N 0 -2 0.4916 -1.9088 2.1022 -5 -4 220 N Y N -5 -4 0.0319 -3.9774 -0.9753 763 590 230 N Y Y 763 590 0.0038 589.9320 173.0592 -763 -590 240 N Y Y -763 -590 0.0010 -590.0178 -173.0423 9 7 260 Y Y N 9 7 -0.0093 6.9855 2.0975 2 1 285 Y Y N 2 1 0.1346 1.0204 0.8239 -9 -7 310 N Y N -9 -7 0.0131 -7.0537 -1.8256 3 3 319 Y Y N 3 3 -0.1683 2.9889 -0.1218 3 -3 320 Y Y Y 3 -3 1.3225 -3.0239 5.9734 1 1 342 Y Y N 1 1 -0.0594 1.0550 -0.0164 -1 -2 352 Y Y N -1 -2 0.3040 -1.9888 0.8781 2 2 360 Y Y N 2 2 -0.1089 1.9333 0.0064Figure 9
Variations of MW, GF, and N2 detection for G05 satellite after the same cycle slips being artificially added to the carrier phase measurement at same time point with four sampling intervals (30s, 15s, 10s, and 1s). (Note: the conditions of big cycle slips detection and repair are not shown in Figure9.)Figure 10
Variations of estimated and predicted values of DID for the G05 satellite with four sampling intervals.Figure 11
Variations of noise level of DID for the G05 satellite with four sampling intervals.Figure 12
Variations of MW, GF, and N2 detection for C13 satellite after the same cycle slips being artificially added to the carrier phase measurement at same time point with four sampling intervals (30s, 15s, 10s, and 1s). (Note: the conditions of big cycle slips detection and repair are not shown in Figure12.)Figure 13
Variations of estimated and predicted values of DID for the C13 satellite with four sampling intervals.Figure 14
Variations of noise level of DID for the C13 satellite with four sampling intervals.(1) From Figures9 and 12, it can be seen that since the same cycle slip is artificially added to the carrier phase measurement at same time point with different sampling intervals, the variations of the three detection processes are the same after performing the correct detection and repair of the cycle slip. That is, the cycle slip can be correctly detected and repaired by using the method proposed in this paper at different sampling intervals. From Tables 3 and 4 and Figures 9 and 12, it can be seen that most of the cycle slip combinations can be detected with the GF detection, but the combination of wavelength ratio reciprocal (such as 763:590 or 9:7) cannot be detected. The combinations of large cycle slip can be preferably detected with the MW detection, but the small cycle slips (less than one cycle) and the same cycle slips cannot be detected. All the cycle slips except the zero cycle in the second frequency can be detected with the N2 detection. All the small, large, continuous, and insensitive cycle slips can be correctly detected and repaired.(2) From Figures10 and 13 and Table 1, it can be seen that the ionospheric delay under a relatively large sampling interval varies greatly, and the estimated and predicted values of DID exceed 0.15 m. The smaller the sampling interval is, the smaller the estimated and predicted values of the ionospheric delay are. For the 1 s sampling interval the predicted values are close to zero and the STD of estimated values is not more than 0.002, which indicate that the influence of ionospheric delay can be ignored in cycle slip detection and repair with a small sampling interval. The estimated and the predicted values of the ionospheric delay of a GPS satellite (G05) are larger than the BDS satellite (C13).(3) From Figures11 and 14 it can be seen that the variance of the DID changes, which preferably reflects the variation degree of the ionospheric delay. That is, where the ionosphere delay varies greatly, the noise level is larger. Meanwhile, the smaller the sampling interval is, the greater the variance of the DID is. The variance is also different for different satellites. Meanwhile, the variance is also different for different sampling intervals and satellites.


## 4. Conclusion

In this paper, the dynamic observation model of DID is correctly constructed. The estimated and predicted values of DID are accurately obtained. Firstly, considering the lack of theoretical model of DID solution, a second-order dynamic observation model of DID is constructed. Then, considering the ionospheric delay itself having strong time-varying characteristics, the LS-VCE method is adopted to adjust the noise level of the DID online to make the level better to approach the actual degrees of change in DID. At the same time, after analysis of the colored noise characteristics of GNSS carrier phase time differential observations, the real-time estimated and predicted values of DID are obtained by designing an adaptive Kalman filter algorithm with a FIR colored noise.Combining the predicted value of DID, real-time detection, and repair of cycle slip in GNSS double-frequency undifferenced phase observations is presented in a stepwise manner. First, GF detection is accurately constructed by using the predicted value. MW detection is calculated by directly rounding the initial MW combined observation. At the same time the GF detection is used to check cycle slip completeness on the MW detection. Second, by considering the correlation between the predicted value of DID and the carrier phase measurement and using the predicted value, a rigorous stochastic model is used to determine the N2 detection. At the same time the GF detection is used to check cycle slip completeness on the N2 detection. Finally, a circulated process of real-time detection and repair of dual-frequency cycle slip is formed, combining the MW, GF, and N2 detection. That is, the carrier phase measurement repaired by the cycle slip is used to calculate the DID, and the predicted values of DID are also used to calculate the cycle slip.The results of experiment with the BDS and GPS observational data at four sampling intervals show the following: (1) the precision of the calculation of the predicted and estimated values of DID with four sampling intervals is higher (both STD less than 1.2 cm), and the two values have a good degree of coincidence (the STD of the difference values is less than 0.6 cm). At the same time, the smaller the interval is, the more accurate the estimated and predicted values of DID are. The larger the sampling interval, the greater the DID and the greater the impact on the detection and repair of the cycle slip. (2) The degree of variation for the DID varies from satellite to satellite. The variation of most GPS and BDS/MEO satellites is larger than BDS/IGSO and BDS/GEO satellites. The maximum value of DID estimation and prediction can reach 0.17m. (3) Compared with the estimated value, the predicted value has better stability and will not be easily deviated from the real situation due to the failure of the cycle slip detection and repair. At the same time, it can better reflect the change trend of DID and has a certain smoothing effect. (4) The online estimated DID variance with the LS-VCE method better reflects the degree of change of DID. (5) The predicted value is very important for N2 detection estimation and the cycle slip completeness checking. (6) The MW, GF, and N2 detection have a high accuracy, and their STDs are 0.083 cycles, 0.004m, and 0.071 cycles, respectively. (7) The proposed method can effectively detect small, large, continuous, and insensitive cycle slips (100% correct rate).If the related ambiguity search algorithm can be used to process the MW and N2 detection and the learning rate in (26) can be reasonably determined, the success rate of cycle slip detection and repair can be further improved. At the same time, we can further study the influence of DID on carrier smoothing pseudo-range and precision single-point positioning.



---

*Source: 1721973-2018-08-29.xml*