Performance improvement for GPS single frequency kinematic relative positioning under poor satellite visibility

Reliable ambiguity resolution in difficult environments such as during setting/rising events of satellites or during limited satellite visibility is a significant challenge for GPS single frequency kinematic relative positioning. Here, a recursive estimation method combining both code and carrier phase measurements was developed that can tolerate recurrent satellite setting/rising and accelerate initialization in motion. We propose an ambiguity dimension expansion method by utilizing the partial ambiguity relevance of previous and current observations. In essence, this method attempts to integrate all useful information into the recursive estimation equation and performs a better least squares adjustment. Using this method, the success rate of the extended ambiguity estimation is independent of the satellite setting and shows robust performance despite poor satellite visibility. Our model allows integration of other useful information into the recursive process. Actual experiments in urban environments demonstrate that the proposed algorithm can improve the reliability and availability of relative positioning.

assumption for most applications given good satellite visibility. However, limited satellite visibility due to the loss of lock may frequently occur due to natural or man-made obstructions such as trees, buildings, bridges, or mountains and this impedes accurate real-time kinematic relative navigation. Additionally, the number of measurement equations may be different for different epochs, leading to different dimensions of the double difference (DD) ambiguities vector. As a consequence, the continuity of the recursive algorithm decreases in schemes that rely on measurements from common satellites of consecutive epochs. Additionally, the information may not be utilized sufficiently, which may prolong its initialization time. To avoid information loss and accelerate initialization, we have developed a new recursive estimator which integrates almost all the available observations for ambiguity resolution. The proposed algorithm can tolerate satellite setting and rising in difficult environments and supports initialization in motion. This is beneficial for some special dynamic applications, which are unable to stop and switch to standby for ambiguity initialization or re-initialization. To achieve reliable ambiguity resolution under poor satellite visibility, we exploited two schemes to strengthen the mathematical model. Firstly, we extended the ambiguity vector by utilizing the partial ambiguity relevance of previous and current observations, and performed a better least squares adjustment. When satellite visibility is limited, it is more reliable to estimate the extended ambiguity vector, since previous observations of invisible satellites may still strengthen the current estimation. Secondly, we designed the algorithm to be efficient, allowing integration of other information to further strengthen the model, such as the integration of the predicted baseline of velocity measurements. We tested the proposed algorithm with actual experiments in urban environments and found improved reliability and availability for kinematic relative positioning.

The recursive mathematical model
For GPS single frequency relative positioning applications, the LAMBDA method is efficient and optimal (Teunissen 1995(Teunissen , 1999Verhagen 2004). In this method, continuous carrier phases through enough epochs must be observed to obtain an accurate float solution. The code measurement gives an approximation of the true range between the receiver and the satellite and it would be beneficial to improve the accuracy of the float solution using both the code and the carrier phase measurements. When processing these code measurements with a least squares method, both functional and stochastic models need to be carefully defined (Rizos 1997;Teunissen and Kleusberg 1998). Many realistic modeling approaches are well documented in the literatures (Langley 1997;Tiberius and Kenselaar 2000;Wang et al. 2002;Li et al. 2008;Luo 2013). For simplicity, we developed our algorithm based on the standard GNSS baseline model.
To distinguish the observations of different time points, we utilized the subscript k to the kth time step, which is referred to as the kth epoch. We assume that the proposed method starts at epoch k, when the number of common GPS visible satellites m k exceeds or is equal to 4 and all the lines of sight (LOS) of two antennas to visible satellites have been calculated successfully for both the base station and rover station. The double-difference (DD) baseline model that combines both the code and carrier phase measurements can be written as: where y is the given GNSS data vector, µ is the noise vector with its variance-covariance (v-c) matrix given by the positive definite matrix Q y , a and b are the integer ambiguity vector and the real-valued baseline vector of order m k − 1 and 3 respectively, and A and B are the given design matrices that link the data vector to the unknown parameters. By applying Cholesky factorization, we decompose the variance-covariance matrix Q y k into Then we calculate the inverse matrix of lower triangular matrix L k :

Multiplying (1) by (3) from the left, we obtain
We then define the following notations: Equation (4) can then be written as Note that the transformed noise vector follows a standard normal distribution. Then let the QR factorization of B k be where Q k is an orthogonal matrix and usually is the product of Householder transformations or Givens rotation matrices (Chang and Paige 2003b), and R k is a 3 × 3 nonsingular upper triangular matrix. Since we assume that m k ≥ 4 so that B k almost always has full column rank, we can partition Note that U k has the same number of rows as R k . After the multiplication on the left of (6) by (8), the following equations can be obtained: Now Eq. (9) can be separated into two parts, and the first one is related to the ambiguity vector and the baseline vector: The second part is only related to the ambiguity vector: (1) The transformed noise vector follows the same distribution because orthogonal transformation does not change the statistical properties of white noise (Chang and Paige 2003b;Chen and Qin 2012): For epoch k we obtain V kȳk ∈ R (2m k −5)×1 and V kĀk ∈ R (2m k −5)×(m k −1) . Note that the columns of V kĀk are equal to the order of the ambiguity vector, and less than or equal to the rows of V kĀk , with the assumption that m k ≥ 4. Next we perform the following orthogonal transformation: where T k is orthogonal and S k is nonsingular upper triangular. The orthogonal transformation can be implemented by a sequence of Householder transformations and the zero matrix vanishes for m k = 4. Next we divide the orthogonal matrix T k into two parts as follows: Multiplying (11) by (14) from the left gives Thus, we can rewrite the top part of (15) as where ŵ k ≡ J k V kȳk ,μ k ≡ J k (V kμk ). Due to the statistical properties of (12), the transformed noise vector still follows the same distribution: By solving the upper triangular system (16), the float solution of ambiguity vector at epoch k can be estimated as Obviously we have Thus, the variance-covariance matrix Qâ k can be written as: (18) and (20), the double difference integer ambiguities can be estimated with the well-known LAMBDA method. Once this has been done successfully, the carrier phase data will act as very precise pseudorange data, allowing accurate baseline solution. The conditional least-squares solution can be written as In general, the reliability of ambiguity resolution is described by the probability of correct integer ambiguity estimation or the so-called ambiguity success rate. For the ILS estimator of the unconstrained GNSS baseline model, Teunissen (1999) showed that P ADOP can be used as an approximation to the ILS success rate, i.e., with n is the dimension of ambiguity vector and the cumulative normal distribution is The ADOP, or Ambiguity Dilution of Precision, is given by Qâ 1 2n , which is a diagnostic that captures the main characteristics of ambiguity precision (Teunissen and Odijk 1997). When the ambiguities are completely decorrelated, the ADOP equals the geometric mean of the standard deviations of the ambiguities, hence it can be considered as a measure of the average ambiguity precision (Verhagen et al. 2013).The success rate of ambiguity resolution is determined by the strength of the underlying GNSS model; the stronger the model, the higher the success rate (Verhagen 2005;Teunissen 2010;Buist 2013;Chen and Li 2014). In order to improve the strength of the single frequency model, the float solution should be estimated with data from more epochs.
Next, we will deduce the recursive least squares approach based on the baseline model (1) of epoch k + 1. The problems of satellite setting and rising are treated with the permutation matrix. This treatment has been well documented in Tiberius (1998) and Chang and Paige (2003b). Here, we employ similar processing in the recursive model, which works very well in our single-frequency case.
At epoch k + 1, despite the differences between the baseline vector b k+1 and b k , we can obtain the following equation using the same processing approach: If the tracking is continued without loss of lock for each satellite, the ambiguity vector will remain at current estimation. However, satellite rising/setting often occurs when the satellite is in motion due to obstructions. If we assume that tracking is continued for the reference satellite, we can compare the ambiguities of epoch k and partition the ambiguity vector a k+1 into two parts by reordering the elements of the ambiguity vector: where the elements of a rem k+1 correspond to the non-reference satellites which are visible at epoch k and remain at epoch k + 1, and those of a new k+1 correspond to the non-reference satellites which appear between epochs k and k + 1. To perform the rearrangement of ambiguities, we construct a permutation matrix G k+1 = G rem k+1 G new k+1 such that where G rem k+1 reorders the ambiguity elements of a k+1 to obtain a rem k+1 and G new k+1 reorders the ambiguity elements of a k+1 to obtain a new k+1 . Note that the permutation matrix is orthogonal as Using (26) and (27), the following equations can be obtained: Applying (28) to (24) gives At this point, we have an equivalent form of (24), which shares common ambiguities with (16), for example, a rem k+1 . To integrate (29) into (16) based on the common ambiguities, we use another permutation matrix k+1 = his k+1 rem k+1 such that where his k+1 reorders the ambiguity elements of a k to obtain a his k+1 and rem k+1 reorders the ambiguity elements of a k to obtain a rem k+1 . Here the elements of a his k+1 correspond to the DD integer ambiguity of non-reference satellites that are not visible before epoch k + 1. In other words, it includes all the 'historical' DD integer ambiguities. Essentially, although some satellites may be invisible during the current epoch, their mathematical equations at previous epoch can be integrated into the current model, since the partial ambiguities remain. By using similar transformations as performed in (27) and (28), the upper triangular system (16) can be written as Then stacking (31) on top of (29) gives This is similar to (13), and we can transform the coefficient matrix of the ambiguity vector to upper triangular form with the following orthogonal transformation: where P k+1 is orthogonal and S k+1 is nonsingular upper triangular with the same number of rows as the ambiguity vector. Note that this 'new' ambiguity vector is dimension-extended and incorporates both the historical and the current DD ambiguities as follows: Multiplying (32) by the orthogonal matrix P T k+1 from the left, we obtain To perform the recursive estimation, we can partition T so that the number of rows of C k+1 is equal to the order of a E k+1 . Thus, we obtain the following upper triangular system: . Note that Eq. (36) is in similar form as (16) and recursive estimation can be achieved by solving the upper triangular system as In this recursive estimator, the ambiguity vector is extended by allowing satellite rising and setting. Once the extended ambiguity vector a E k+1 is resolved, we can extract the current ambiguity vector a k+1 and calculate the current conditional baseline vector in the same way as (21): The recursive process can always be achieved if m k+1 ≥ 4 and if a rem k+1 is not empty. To more clearly demonstrate the rearrangement of the ambiguities, we have presented a case in Table 1 which shows how the extended ambiguity vector changes as satellite setting/rising occurs. The notation PRN epoch denotes the DD integer ambiguity of the nonreference satellite PRN to the reference satellite with its 'birth moment' denoted with the superscript 'epoch' . Table 1 shows that satellite setting cannot change the number of ambiguities.
The dimension of the extended ambiguity vector cannot be reduced by satellite setting. However, the dimension increases as the number of visible satellites increases. Thus, the progressive increase of ambiguities may increases the burden of computation. To deal with this problem during real-time application, a fast integer least-squares estimation is required for high-dimensional ambiguity resolution, which has been described previously (Chang et al. 2005;Zhou 2010;Jazaeri et al. 2012). Another simple approach is to limit the growth of the dimension of the historical DD ambiguity vector a his k+1 . For instance, when it  exceeds the preset threshold, we can remove the first ambiguity element of the upper triangular system (36), thus reducing the total ambiguity dimension by 1. Here, the threshold was set as 10 and little effect on the global success rate was found.

Ambiguity validation for ambiguity initialization
The final aim of initialization is to perform a rapid and successful ambiguity resolution for real-time application. Thus, the validation procedure is required for initialization. We used the ratio test, a commonly used validation procedure described in Verhagen (2004) and Teunissen and Verhagen (2009). In this method, the ratio of the squared norm of the ambiguity residuals in the metric of the covariance for the best and second best integer solution, as obtained from LAMBDA, is examined to perform a validation calculation: where ⌣ a E k+1 and ã E k+1 denote the best and second best integer solution, respectively. The critical value is often accomplished using empirical results. Many software packages set this value as 3 (Leick 2004;Kroes 2006;Chen et al. 2015). Here, we employed the same value for the ambiguity validation, which worked adequately in our single-frequency case. Note that when the correct ambiguity candidate is resolved, the norm of the noise term of upper triangular system (36) is equal to the squared norm of the ambiguity residuals from LAMBDA, according to the following derivation: Because the noise error term obeys a standard normal distribution with a mean of zero, the residual term of the correct candidate will maintain normality. The Kolmogorov-Smirnov residual test (K-S test) can be used to compare a sample with a reference normal probability distribution. The K-S statistic quantifies the distance between the empirical distribution function of the sample and the cumulative distribution function of the reference distribution (Marsaglia et al. 2003). Here the null hypothesis for the K-S test is that μ k+1 has a normal distribution N (0, I n ) where n is the dimension of μ k+1 . The alternative hypothesis is that μ k+1 does not have that distribution. If one cannot reject the null hypothesis, the integer ambiguity vector and the baseline vector may be correct (Chen and Qin 2012). The squared norm of the optimal ambiguity candidate residuals should be checked by a user-defined Chi square threshold as follows: where α is the probability of rejecting the null hypothesis given that it is true. Once the test above is rejected, it indicates that the measurements may not fit the model of the null hypothesis. Possible explanations such as cycle slips and abnormal data or systematic error should be identified. A good overview of the commonly used quality control procedures and testing theory is presented in Teunissen and Kleusberg (1998) and Teunissen (2000).

Using the structure to improve the model strength with the baseline state prediction
Once the initialization is achieved by validation, all subsequent relative position estimations can be done with the continuous tracking. However, the resolved ambiguity cannot be used for the current baseline estimation if the corresponding satellite is invisible. We designed the algorithm to be efficient and to use the mathematical structure of the recursive estimation to allow the possible integration of other available measurements. For example, in order to improve the success rate under poor satellite visibility, we can integrate the resolved baseline of last epoch into the recursive model. When initialization is achieved, we obtain the first precise estimation of baseline and then predict its value for the next epoch based on the velocity information. For most GPS receivers, accurate position information is provided as well as 3D velocity. As is shown in Fig. 2, the baseline vector b is defined as a vector from reference antenna A to another antenna B, and the measurements of the velocity vectors of both antennas are given as v A k and v B k at epoch k. Assume that the time span is very short and the accuracy of velocity measurements is sufficiently high, the baseline vector at epoch k + 1 can be predicted as follows: where ⌣ b k ⌣ a k is the fixed ambiguity vector at epoch k and w k,k+1 is the measurement noise vector. The variance-covariance matrix can be obtained as Fig. 2 Motion of the baseline So far, if we define p k+1 ≡ y k+1 − B k+1bk+1 , the following equation can be derived from (1) The v-c matrix is written as By applying the Cholesky factorization, we decompose the v-c matrix Q p k+1 into Multiplying (48) by inverse of L p k+1 from the left and combining it with (36) gives We are still able to obtain the upper triangular system by applying the orthogonal transformation on both side of (51) in a similar method as used in (35). The approach above demonstrates the flexibility of the recursive model. There are many approaches that use velocity to predict the baseline, such as the state-space models and Kalman filter (Chang and Huang 2005;Wang and Lai 2015). In these approaches, the user motion can be modeled using one of many common dynamic models (random walk, constant velocity, 1st order GM, etc.) while the ambiguities can be modelled as random constants.

Results of kinematic experiment
The advantage of this proposed method is that the initialization phase can be achieved in motion, since only (11) is utilized to construct the recursive model, which is independent of the baseline status. The success rate of integer ambiguity estimation is also insensitive to satellite setting. To test the actual performance of this extended ambiguity resolution technique with limited satellite visibility, two experiments were performed in urban environments. Both experiments used the receiver onboard the remote control car (telecar) as the rover, shuttling in the street by receiving remote control commands. The reference station was located at a crossroad. A JAVAD Alpha 2 GPS receiver was used, as shown in Fig. 3. Since the advanced multipath reduction techniques are utilized −1 u k Fig. 3 The JAVAD Alpha2 GPS receiver for both carrier phase and code measurements, almost all anomalies and satellite multipath were removed with Alpha 2. GPS antenna design can further reduce the effect of multipath. Here we used two Trimble ® Zephyr ™ 2 antennas that show excellent performance in multipath mitigation, as shown in Fig. 4. For the first experiment, we controlled the telecar moving towards the north. A bird's eye view of the study area is given in Fig. 5. There are many trees and buildings on both sides of the street, especially in the middle segment. Figure 6 shows the resolved baseline north/east/up components and the baseline length. The ambiguity validation was achieved at epoch 67, marked with the vertical line in Fig. 6. Before this epoch, we estimated the baseline with the float ambiguity vector. After that, the recursive process continued and the optimal ambiguity vector of LAMBDA method was utilized to calculate the baseline. Although the proposed strengthening scheme was not exploited in this experiment, the resolved east, north and up baseline components and baseline length are consistent with the movement of the vehicle. Small fluctuations were found in the east and up components due to rough pavement and inaccurate control. Figure 7 shows the global/local success rate and the number of visible satellites during the first experiment. As is shown, the global success rate was higher than the local success rate and was not sensitive to the drop of the number of visible satellites, even when the number of satellites drops to five. However, unlike the global success rate, the local  of the extended ambiguity vector is more reliable. The lower panel of Fig. 7 shows that the dimension of the extended ambiguity vector increases progressively as satellite rising occurs. Recurrent setting and rising events of satellites will accelerate the growth of the dimension of the extended ambiguity vector. To reduce the necessary computation, the threshold is set as 10 which showed little effect on the global success rate.
In our second experiment, we investigated the actual performance of the proposed method in a more challenging environment. A bird's eye view of the study area is provided in Fig. 8. The telecar moved towards the east, gradually approaching the reference station. Compared with the first experiment, there were more buildings and trees along the road. Once initialization was achieved, the proposed model-strengthening scheme was exploited during the subsequent epochs. The global/local success rate and the number of visible satellites are shown in Fig. 9, and the vertical line indicates the end of the initialization. The top panel of Fig. 9 shows that the local success rate increased dramatically when the predicted baseline equation was integrated into the recursive algorithm, see (51). This modification made the model more 'robust' to satellite setting and the recovery from decline required only a few seconds. The middle panel of Fig. 9 shows the number of visible satellites. Compared with the first experiment, satellite setting and rising was more frequent and the number of visible satellites dropped to 4 many times. However, the global success rate remained close to 1, and was insensitive to the poor satellite visibility. Figure 10 shows the resolved baseline north/east/up components and the baseline length of the second experiment. There was no obvious failure of the baseline estimation, and the inclusion of this estimation in the model allows increased reliability of ambiguity resolution and improved continuity and availability.

Results of static experiment
Validation of corrected position data is not done or possible for the kinematic tests. Thus, we instead used a static test with two JAVAD Alpha 2 receivers to validate the ambiguity-fixed positions by comparison to a known ground-truth position. Data were collected at 1 Hz with a zero cut-off elevation angle for a total of 1408 epochs logged. To obtain a known ground-truth relative position, the distance between the centers of both antennas was measured and the baseline length was approximately 30.08 m. Both antennas were placed on one side of a wide road near the south-north direction. The up component of baseline was close to zero and the north component of baseline was close to  includes epoch 1 to 60, including the initialization phase (Figs. 11,12), and the second one includes epochs 61 to 1408 (Figs. 13, 14). Figures 11 and 13 show the global/local success rate and the number of visible satellites (the vertical line denotes that the initialization is achieved), and Figs. 12 and 14 show the baseline components and resolved baseline length. With the fixed ambiguities, the average value of north, east, and up measurements of the baseline vector were 30.0857, 0.2497, −0.0466 meters, with standard deviations of 0.0052, 0.0024, 0.0095 meters, respectively, indicating a sub-centimeter level relative position accuracy.
The average of the vertical baseline component was close to zero, but the precision was poorer than the horizontal baseline components because only satellites above the horizon are tracked. The resolved baseline length was very close to the measured value of 30.08 m and the resolved north component was much larger than the east component, consistent with the baseline placement.Thus, the ambiguity-fixed relative position was validated with sub-centimeter level precision.

Conclusions
A recursive least squares estimator was presented for reliable GPS single frequency kinematic relative positioning in difficult environments. The initialization can be achieved in motion, even if satellite setting and rising events occur frequently during this phase. In order to improve the success rate of ambiguity resolution, an ambiguity dimension expansion method was developed by integrating previous and available information into the current ambiguity estimation equation. The recursive model can be further strengthened using the predicted baseline model or the state-space model. The success rates and relative positioning performance were measured in actual experiments in urban environments. The results reveal that reliable ambiguity estimation can be achieved in motion during initialization. The proposed algorithm also increases the continuity and availability level of relative positioning, even if the number of visible satellites frequently decreases to 4 during the observation.

Author information
Wantong Chen is a lecturer in Tianjin Key Lab for Advanced Signal Processing at Civil Aviation University of China. He received his Ph.D. from Beihang University in 2013, for work on GNSS high-precision attitude determination applications. His current area of research interest is the use of GNSS relative positioning technique for spacecraft rendezvous/docking and aircraft precision landing.