Performance Improvement of a Floating Solution Using a Recursive Filter

Abstract

In CDGPS, ambiguity resolution is determined by the performance of a floating solution, and thus, the performance needs to be improved. In the case of precise positioning at a stationary position, the batch method using multiple measurements is used for the accuracy improvement of a position. The position accuracy performance of a floating solution is outstanding, but it has a problem of high computation cost because all measurements are used. In this study, to improve the floating solution performance of the initial static user in CDGPS, a floating solution method using a recursive filter was implemented. A recursive filter estimates the position solution of the current epoch using the position solution of up to the previous epoch and the pseudorange measurement of the current epoch. The computation cost of the floating solution method using a recursive filter was found to be similar to that of the epoch-by-epoch method. Also, based on actual GPS signals, the floating solution performance was found to be similar to that of the batch method. The floating solution using a recursive filter could significantly improve the performance of the prompt initial position and ambiguity resolution of the initial static user.

1. INTRODUCTION

A precise positioning method using carrier phase measurements (CDGPS: Carrier phase Differential GPS) shows outstanding navigation performance (several cm ~ dozens of mm level) (Kaplan & Hegarty 2006), but has a problem of ambiguity resolution (Kaplan & Hegarty 2006). CDGPS is organized in the order of floating solution, ambiguity searching, and integer solution (Leick 2004). The floating solution is a stage in which a position solution and an ambiguity in the floating domain are determined ignoring the integer condition of an observation equation ambiguity. The ambiguity searching determines an ambiguity that satisfies the integer condition of the observation equation, from the ambiguity in the floating domain and covariance information, through searching. The integer solution is a stage in which a position solution in the integer domain is redetermined using the ambiguity determined by the searching. In the navigation field, a floating solution is determined by the epoch-by-epoch method (de Jonge et al. 2000), which has low computation cost. As it uses only the measurement of the current epoch, the computation cost is low, but the position accuracy decreases. However, in the geodetic survey field, to improve position accuracy performance, the batch method (Park 2001), which uses measurements of all epochs, is utilized. It uses all the measurements from the previous epoch to the current epoch, and thus has a problem of exponentially increasing computation cost.

Least SQuare (LSQ) including a recursive filter has been generally studied in many fields. In this study, a recursive filter was used for the determination of a floating solution that includes both a position solution and an ambiguity in the floating domain of CDGPS. For the floating solution using a recursive filter, a floating solution, which uses the position solution of the previous epoch and the pseudorange measurement of the current epoch, was defined based on a stationary user. The results indicated that the position accuracy was improved compared to the epoch-by-epoch method, which uses only the pseudorange measurement of the current epoch; and that the computation cost was similar to that of the case using a single method, resolving the problem of exponentially increasing computation cost in the batch method, which uses all measurements.

2. FLOATING SOLUTION IN CDGPS

For the floating solution of CDGPS, a position solution and an ambiguity in the floating domain are estimated using double differenced code measurement and carrier phase measurement. For the double differenced measurements in 0 baseline environment, common errors included in the measurements (satellite orbit, satellite clock error, ionospheric delay, tropospheric delay, and receiver clock error) are eliminated. Therefore, the double differenced code and carrier phase measurement observation equation of the \(k\)-th epoch can be expressed as Eq. (1).

\(\Psi(k)=r(k)+v(k)\)                                                                                                            

\(\Phi(k)=r(k)+\lambda a+w(k)\)                                                                                               (1)

where \(r(k)\) is the double differenced value of the calculated geometric distance between the satellite and the user, satellite and the user, \(\lambda\) is the wavelength, \(a\) is the double differenced ambiguity, and \(v(k)\) and \(w(k)\) are the double differenced code and carrier phase measurement errors. Eq. (2) can be obtained by the linearization of Eq. (1).

\(\rho(k)\equiv\Psi(k)-r_0(k)=H(k)b(k)+v(k)\)                                                                                               

\(l(k)\equiv\Phi(k)-r_0(k)=H(k)b(k)+\lambda a+w(k)\)                                                                                   (2)

where \(\rho(k)\) and \(l(k)\) are the double differenced code and carrier phase pseudorange measurements at the linearization reference point, \(r_0(k)\) is the calculated double differenced distance between the satellite and the linearization reference point, \(H(k)\) is the single differenced matrix of the line-of-sight vector from the linearization reference point to the satellite, and \(b(k)\) is the baseline vector between the reference station and the user. Eq. (2) can be expressed in a vector-matrix format, as shown in Eq. (3) (Han et al. 2009).

\(\left[\begin{matrix}l(k)\\\rho(k)\\\end{matrix}\right]=\left[\begin{matrix}H(k)\\H(k)\\\end{matrix}\right]b(k)+\left[\begin{matrix}\lambda I\\0\\\end{matrix}\right]a+\left[\begin{matrix}w(k)\\v(k)\\\end{matrix}\right]\)                                                                                     (3)

In this study, CDGPS consists of floating solution, ambiguity resolution, and integer solution stages. To improve the performance of a floating solution among the three stages, an initial position solution and an ambiguity in the floating domain of a stationary user in the existing epoch-by-epoch method that uses a single measurement were implemented using a recursive filter.

3. FLOATING SOLUTION USING RECURSIVE FILTER

For a floating solution using the batch method, a position solution and an ambiguity in the floating domain are obtained using all the measurements from the previous epoch to the current epoch. For the measurement of the N-epoch, N position vectors and one ambiguity vector are obtained if a user is in motion. However, if a user is stationary, it becomes a problem of obtaining one position vector and one ambiguity vector. In this study, ambiguity resolution is mostly performed at the early stage of navigation; and if an ambiguity is once determined, a fixed value is used or only a value that changes in combination with the epoch-by-epoch method can be obtained. Therefore, a floating solution method using a recursive filter for improving the initial floating solution performance of a stationary user is proposed. In the case of a stationary user, a measurement model for the N-epoch measurement can be expressed as Eq. (4). The subscript \(\underline{b}\) represents a solution obtained by the batch method (Park 2001).

\(\left[\begin{matrix}l(1)\\\rho(1)\\\vdots\\l(N)\\\rho(N)\\\end{matrix}\right]=\left[\begin{matrix}H(1)&\lambda I\\H(1)&0\\\vdots&\vdots\\H(N)&\lambda I\\H(N)&0\\\end{matrix}\right]\left[\begin{matrix}b\\a\\\end{matrix}\right]_{\underline{b}}+\left[\begin{matrix}w(1)\\v(1)\\\vdots\\w(N)\\v(N)\\\end{matrix}\right]\)                                                                           (4)

A floating solution can be calculated by applying LSQ to Eq. (4). However, if there are M double differenced pseudorange measurements (M+1 satellites), 2MN measurements need to be processed, and thus, the computation cost abruptly increases depending on the increase in the measurements.

In this study, to use a recursive filter for improving position accuracy and reducing computation cost, the observation equation of the N-th measurement using the floating solution of the previous epoch (k) and the measurement of the current epoch (k+1) was obtained as shown in Eq. (5). The subscript represents a solution obtained using a recursive filter.

\(\left[\begin{matrix}l(k+1)\\\rho(k+1)\\{\hat{b}}_r(k)\\{\hat{a}}_r(k)\\\end{matrix}\right]=\left[\begin{matrix}H(k+1)&\lambda I\\H(k+1)&0\\I&0\\0&I\\\end{matrix}\right]\left[\begin{matrix}b(k+1)\\a(k+1)\\\end{matrix}\right]_{\underline{r}}+\left[\begin{matrix}w(k+1)\\v(k+1)\\\eta_b(k)\\\eta_a(k)\\\end{matrix}\right]\)                                                                 (5)

where \(\left[\begin{matrix}\hat{b}(k)\\\hat{a}(k)\\\end{matrix}\right]_{\underline{r}}\) is the floating solution of up to the previous epoch (k), and \(\left[\begin{matrix}l(k+1)\\\rho(k+1)\\\end{matrix}\right]\) is the double differenced pseudorange measurement of the current epoch (k+1). The covariance of the double differenced pseudorange measurement is \(\left[\begin{matrix}w(k+1)\\v(k+1)\\\end{matrix}\right] \sim N \left ( 0,\left[\begin{matrix}\sigma_l^2Q&0\\0&\sigma_\rho^2Q\\\end{matrix}\right] \right ) \equiv N(0,Q_y)\), and the covariance , and the covariance of the floating solution is \(\left[\begin{matrix}\eta_b(k)\\\eta_a(k)\\\end{matrix}\right] \sim N \left ( 0,\left[\begin{matrix}Q_{\hat{b}}(k)&Q_{\hat{b}\hat{a}}(k)\\Q_{\hat{a}\hat{b}}(k)&Q_{\hat{a}}(k)\\\end{matrix}\right] \right )\equiv N(0,Q_x(k))\). Eq. (7) can be obtained by expressing Eq. (5) in a simplified format as shown in Eq. (6) and by applying LSQ.

\(\left[\begin{matrix}l(k+1)\\\rho(k+1)\\{\hat{b}}_r(k)\\{\hat{a}}_r(k)\\\end{matrix}\right]=\left[\begin{matrix}G(k+1)\\I\\\end{matrix}\right]\left[\begin{matrix}b(k+1)\\a(k+1)\\\end{matrix}\right]_{\underline{r}}+\left[\begin{matrix}w(k+1)\\v(k+1)\\\eta_b(k)\\\eta_a(k)\\\end{matrix}\right]\)                                                                     (6)

\(\left[\begin{matrix}\hat{b}(k+1)\\\hat{a}(k+1)\\\end{matrix}\right]_{\underline{r}}=(G(k+1)^TQ_y^{-1}G(k+1)+Q_x^{-1}(k))^{-1}G(k+1)^TQ_y^{-1}\left[\begin{matrix}l(k+1)\\\rho(k+1)\\\end{matrix}\right] + (G(k+1)^TQ_y^{-1}G(k+1)+Q_x^{-1}(k))^{-1}Q_x^{-1}(k)\left[\begin{matrix}\hat{b}(k)\\\hat{a}(k)\\\end{matrix}\right]_{\underline{r}}\)                

\(=(G(k+1)^TQ_y^{-1}G(k+1)+Q_x^{-1}(k))^{-1}(G(k+1)^TQ_y^{-1}G(k))\left[\begin{matrix}\hat{b}(k+1)\\\hat{a}(k+1)\\\end{matrix}\right]_{\underline{e}} (+(G(k+1)^TQ_y^{-1}G(k+1)+Q_x^{-1}(k))^{-1}Q_x^{-1}(k)\left[\begin{matrix}\hat{b}(k)\\\hat{a}(k)\\\end{matrix}\right]_{\underline{r}}\)   (7)

Eq. (7) is represented by the weighted sum of the estimated value obtained using only the measurement of the current epoch (k+1), \(\left[\begin{matrix}\hat{b}(k+1)\\\hat{a}(k+1)\\\end{matrix}\right]_{\underline{e}}\), and the estimated value obtained using the measurement of up to the previous epoch (k), \(\left[\begin{matrix}\hat{b}(k)\\\hat{a}(k)\\\end{matrix}\right]_{\underline{r}}\). Also, in Eq. (8), the sum of the weights becomes an identity matrix.

\((G(k+1)^TQ_y^{-1}G(k+1)+Q_x^{-1}(k))^{-1}((G(k+1)^TQ_y^{-1}G(k+1))+Q_x^{-1}(k))=I\)                                                                 (8)

In Eq. (9), the covariance of the estimated value is \(Q_x(k+1)=(G(k+1)^TQ_y^{-1}G(k+1)+Q_x^{-1}(k))^{-1}\).

\({\rm cov} \left ( \left[\begin{matrix}\hat{b}(k+1)\\\hat{a}(k+1)\\\end{matrix}\right]_{\underline{r}}\right ) =Q_x(k+1)G(k+1)^TQ_y^{-1}Q_yQ_y^{-1}G(k+1)Q_x(k+1)+Q_x(k+1)Q_x^{-1}(k)Q_x(k)Q_x^{-1}(k)Q_x(k+1)\)                               

\(=Q_x(k+1)[G(k+1)^TQ_y^{-1}G(k+1)+Q_x^{-1}(k)]Qx(k+1)\)                                                                                          

\(=Q_x(k+1)Q_x^{-1}(k+1)Q_x(k+1)=Q_x(k+1)\)                                                                                                      (9)

Based on the above results, a floating solution using a recursive filter can be organized as follows.

Summary of Recursive Algorithm for Stationary User

1) F or the measurement of the first epoch, \(\left[\begin{matrix}\hat{b}(1)\\\hat{a}(1)\\\end{matrix}\right]_{\underline{e}},Q_x(1)\) is calculated by applying LSQ to Eq. (3). 2) F or the measurement of the following epoch, a covariance and an estimated value are calculated using Eq. (10) and (11).

\(Q_x^{-1}(k+1)=G(k+1)^TQ_y^{-1}G(k+1)+Q_x^{-1}(k)\)                                                                                   (10)

\(\left[\begin{matrix}\hat{b}(k+1)\\\hat{a}(k+1)\\\end{matrix}\right]_{\underline{r}}=Q_x(k+1)G(k+1)^TQ_y^{-1}\left[\begin{matrix}l(k+1)\\\rho(k+1)\\\end{matrix}\right]+Q_x(k+1)Q_x^{-1}(k)\left[\begin{matrix}\hat{b}(k)\\\hat{a}(k)\\\end{matrix}\right]_{\underline{r}}\)                                                         (11)

In this study, a floating solution using a recursive filter utilizes the position solution of up to the previous epoch and the pseudorange measurement of the current epoch. It has the advantage of constant computation cost, and an accuracy improvement is expected similar to the case using all measurements from the position solution of the previous epoch. Therefore, for the performance improvement of a floating solution using a recursive filter, improvements in the ambiguity searching result at an initial stationary position and the integer solution performance of CDGPS could be expected.

4. PERFORMANCE ANALYSIS

In this study, the numerical example, navigation performance, and computation cost of the floating solution method was examined using actual GPS signals. For the collected signals, actual GPS signals were utilized using two NovAtel Propak-V3 GNSS receivers (NovAtel, Inc. 2013) at Chungnam National University. The distance between the two antennas was a short-distance baseline (5.57 m), and 10 visible satellites were observed. The navigation performance was examined using an L1 C/A code error standard deviation of 20 cm \((\sigma_\rho^2=(0.2)^2)\) and an L1 carrier of 2 mm \((\sigma_l^2=(0.002)^2)\).

4.1 Numerical example

The linearized double differenced code and carrier measurements and measurement matrix are as follows.

\(l(1)=\left[\begin{matrix}\mathrm{\ \ }0.663\\\mathrm{\ \ }0.572\\-0.190\\-0.662\\\mathrm{\ \ }1.041\\\mathrm{\ \ }0.192\\-5.426\\\mathrm{\ \ }4.090\\-0.952\\\end{matrix}\right],\)     \(l(2)=\left[\begin{matrix}\mathrm{\ \ }0.663\\\mathrm{\ \ }0.570\\-0.189\\-0.663\\\mathrm{\ \ }1.043\\\mathrm{\ \ }0.191\\-5.418\\\mathrm{\ \ }4.086\\-0.951\\\end{matrix}\right],\)     \(\rho(1)=\left[\begin{matrix}\mathrm{\ \ }0.568\\-0.205\\-0.069\\\mathrm{\ \ }0.254\\-0.181\\-0.147\\\mathrm{\ \ }0.036\\\mathrm{\ \ }0.165\\\mathrm{\ \ }0.039\\\end{matrix}\right],\)     \(\rho(2)=\left[\begin{matrix}\mathrm{\ \ }0.567\\-0.212\\-0.063\\\mathrm{\ \ }0.267\\-0.188\\-0.148\\\mathrm{\ \ }0.041\\\mathrm{\ \ }0.153\\\mathrm{\ \ }0.032\\\end{matrix}\right]\)

\(H(1)=\left[\begin{matrix}-0.672&-1.426&\mathrm{\ \ }0.507\\-0.323&-0.046&-0.577\\\mathrm{\ \ }0.634&\mathrm{\ \ }0.399&-0.445\\-1.145&-0.126&\mathrm{\ \ }0.369\\\mathrm{\ \ }1.419&\mathrm{\ \ }0.997&-0.276\\\mathrm{\ \ }0.039&-0.734&\mathrm{\ \ }0.994\\-0.545&-0.311&\mathrm{\ \ }0.202\\-0.488&\mathrm{\ \ }0.228&-1.450\\\mathrm{\ \ }1.219&\mathrm{\ \ }0.390&\mathrm{\ \ }0.945\\\end{matrix}\right],\)     \(H(2)=\left[\begin{matrix}-0.672&-1.426&\mathrm{\ \ }0.507\\-0.323&-0.046&-0.577\\\mathrm{\ \ }0.634&\mathrm{\ \ }0.399&-0.445\\-1.145&-0.126&\mathrm{\ \ }0.369\\\mathrm{\ \ }1.419&\mathrm{\ \ }0.997&-0.276\\\mathrm{\ \ }0.039&-0.734&\mathrm{\ \ }0.994\\-0.545&-0.311&\mathrm{\ \ }0.202\\-0.488&\mathrm{\ \ }0.228&-1.450\\\mathrm{\ \ }1.218&\mathrm{\ \ }0.390&\mathrm{\ \ }0.945\\\end{matrix}\right]\)

 

Based on the floating solution method, the position solution and ambiguity in the floating domain were compared using actual measurements. For the position solution in the floating domain, initial epoch-by-epoch method is performed in the first measurement, and thus, the same position solution error and floating ambiguity result were observed.

\({\hat{b}}_e(1)=\left[\begin{matrix}\mathrm{\ \ }0.0014\\-0.0016\\-0.0008\\\end{matrix}\right], \)  \({\hat{b}}_b(1)=\left[\begin{matrix}\mathrm{\ \ }0.0014\\-0.0016\\-0.0008\\\end{matrix}\right],\)  \({\hat{b}}_r(1)=\left[\begin{matrix}\mathrm{\ \ }0.0014\\-0.0016\\-0.0008\\\end{matrix}\right]\)

\({\hat{a}}_e(1)=\left[\begin{matrix}\mathrm{\ \ \ }0.7631\\\mathrm{\ \ \ \ }2.7168\\\ -0.6008\\\ -2.3422\\\mathrm{\ \ \ \ \ }5.1627\\\mathrm{\ \ \ \ \ }0.6790\\-27.3751\\\mathrm{\ \ \ }19.6656\\\ -4.0662\\\end{matrix}\right],\)  \({\hat{a}}_b(1)=\left[\begin{matrix}\mathrm{\ \ \ }0.7631\\\mathrm{\ \ \ \ }2.7168\\\ -0.6008\\\ -2.3422\\\mathrm{\ \ \ \ \ }5.1627\\\mathrm{\ \ \ \ \ }0.6790\\-27.3751\\\mathrm{\ \ \ }19.6656\\\ -4.0662\\\end{matrix}\right],\)  \({\hat{a}}_r(1)=\left[\begin{matrix}\mathrm{\ \ \ }0.7631\\\mathrm{\ \ \ \ }2.7168\\\ -0.6008\\\ -2.3422\\\mathrm{\ \ \ \ \ }5.1627\\\mathrm{\ \ \ \ \ }0.6790\\-27.3751\\\mathrm{\ \ \ }19.6656\\\ -4.0662\\\end{matrix}\right]\)

In the second measurement, the batch method and the floating solution method using a recursive filter had the position solution error and floating ambiguity result that are identical down to four decimal places, unlike the epoch-byepoch method.

\({\hat{b}}_e(2)=\left[\begin{matrix}-0.0106\\\mathrm{\ \ }0.0090\\\mathrm{\ \ }0.0075\\\end{matrix}\right], \)  \({\hat{b}}_b(2)=\left[\begin{matrix}-0.0049\\\mathrm{\ \ }0.0041\\\mathrm{\ \ }0.0037\\\end{matrix}\right],\)  \({\hat{b}}_r(2)=\left[\begin{matrix}-0.0049\\\mathrm{\ \ }0.0041\\\mathrm{\ \ }0.0037\\\end{matrix}\right]\)

\({\hat{a}}_e(2)=\left[\begin{matrix}\mathrm{\ \ \ \ \ }0.7745\\\mathrm{\ \ \ \ \ \ }2.7119\\\mathrm{\ \ }-0.5586\\\mathrm{\ \ }-2.4283\\\mathrm{\ \ \ \ \ }5.2181\\\mathrm{\ \ \ \ \ }0.6723\\-27.3831\\\mathrm{\ \ \ }19.6666\\\ -4.0511\\\end{matrix}\right],\)  \({\hat{a}}_b(2)=\left[\begin{matrix}\mathrm{\ \ \ \ \ }0.7697\\\mathrm{\ \ \ \ \ \ }2.7149\\\mathrm{\ \ }-0.5758\\\mathrm{\ \ }-2.3879\\\mathrm{\ \ \ \ \ }5.1915\\\mathrm{\ \ \ \ \ }0.6755\\-27.3799\\\mathrm{\ \ \ }19.6673\\\ -4.0589\\\end{matrix}\right],\)  \({\hat{a}}_r(2)=\left[\begin{matrix}\mathrm{\ \ \ \ \ }0.7697\\\mathrm{\ \ \ \ \ \ }2.7149\\\mathrm{\ \ }-0.5758\\\mathrm{\ \ }-2.3879\\\mathrm{\ \ \ \ \ }5.1915\\\mathrm{\ \ \ \ \ }0.6755\\-27.3799\\\mathrm{\ \ \ }19.6673\\\ -4.0589\\\end{matrix}\right]\)

It was numerically shown that the error of the floating solution method using a recursive filter proposed in this study was the same as that of the batch method, and that the floating solution method was superior to the epoch-byepoch method.

4.2 Floating position result

To examine the floating solution performance of CDGPS, it was compared with a true position using the same 300 epoch measurements. Figs. 1, 2, and 3 show the comparison of the X-, Y-, and Z-axis navigation errors for the epoch-byepoch method, batch method, and floating solution method using a recursive filter, respectively.

Fig. 1. X-Axis floating position error result.

Fig. 2 Y-Axis Floating Position Error Result.

Fig. 3. Z-Axis floating position error result.

The floating solution method using a recursive filter proposed in this study had superior floating solution performance compared to the existing epoch-by-epoch method, and had nearly identical navigation performance compared to the batch method.

4.3 Computation cost result

The floating solution method using a recursive filter proposed in this study is an algorithm for reducing the computation cost of the batch method. To compare the computation cost, the computation time of the floating solution was examined using the function of Matlab. The epoch-by-epoch & recursive filter methods had the same computation cost in every epoch. However, for the batch method, in the case of the same satellite in the N-th epoch, there were N times more measurements, and thus, the computation cost increased exponentially. Fig. 4 shows the computation time for every measurement depending on the floating solution method.

Fig. 4. Computation cost result.

As shown in the figure, for the batch method, the measurement size (matrix computation) exponentially increased as the epoch increased. The epoch-by-epoch method and the floating solution method using a recursive filter had the same computation time in every epoch. The average computation times of the epoch-by-epoch and recursive method were similar: 0.005747 seconds and 0.005968 seconds per epoch, respectively. However, for the batch method, the computation time of the 300th epoch was about 33 seconds, which increased substantially. The floating solution method using a recursive filter used in this study had similar computation time to that of the epoch-byepoch method, and it is thought to be appropriate for initial ambiguity resolution or real-time use.

5. CONCLUSIONS

To determine accurate initial position and ambiguity, a stationary user utilizes multi-epoch measurements. For the batch method using the existing multi epoch, as the number of multi-epoch measurements increases, the computation cost increases exponentially due to the increase in the matrix size. Thus, it is difficult to use the method for a dynamic user, and it has been mostly used for post-processing precise positioning. The floating solution method using a recursive filter used in this study had performance effects using multi-epoch measurements based on the computation cost reduction and floating solution accuracy improvement of the existing epoch-by-epoch and batch method. The results of the performance analysis indicated that the performance of the floating solution method was identical to that of the batch method, and that the computation cost was similar to that of the existing epoch-by-epoch method although the multi-epoch measurements increased. In the future, the floating solution method using a recursive filter could significantly contribute to prompt initial position and ambiguity resolution of a dynamic user, and further studies on floating solution design using a recursive filter that can be used in a dynamic scenario after the initial position of a dynamic user are required.