FAST INTEGER AMBIGUITY RESOLUTION IN GPS KINEMATIC POSITIONING USING LEFT NULL SPACE AND MULTI-TIME ( INVERSE ) PAIRED CHOLESKY

Aiming at the problems that huge amount of computation in ambiguity resolution with multiple epochs and high-order matrix inversion occurred in the GPS kinematic relative positioning, a modified algorithm for fast integer ambiguity resolution is proposed. Firstly, Singular Value Decomposition (SVD) is applied to construct the left null space matrix in order to eliminate the baselines components, which is able to separate ambiguity parameters from the position parameters efficiently. Kalman filter is applied only to estimate the ambiguity parameters so that the real-time ambiguity float solution is obtained. Then, sorting and multi-time (inverse) paired Cholesky decomposition are adopted for decorrelation of ambiguity. After diagonal elements preprocessing and diagonal elements sorting according to the results of Cholesky decomposition, the efficiency of decomposition and decorrelation is improved. Lastly, the integer search algorithm implemented in LAMBDA method is used for searching the integer ambiguity. To verify the validity and efficacy of the proposed algorithm, static and kinematic tests are carried out. Experimental results show that this algorithm has good performance of decorrelation and precision of float solution, with computation speed also increased effectively. The final positioning accuracy result with static baseline error less than 1 cm and kinematic error less than 2 cm, which indicates that it can be used for fast kinematic positioning of high precision carrier. 833 Fast integer ambiguity... Bol. Ciênc. Geod., sec. Artigos, Curitiba, v. 21, n 4, p.832-847, out-dez, 2015.


Introduction
Fast ambiguity resolution is critical to GPS carrier phase measurements for high precision kimematic positioning (Hofmann-wellenhoff et al., 2001;Leick, 2004), which has been studied by many researchers during the last two decades.Examples of proposed methods were dual-frequency P code pseudo range algorithm (Blewitt, 1989), least squares (LS) searching algorithm (Hatch, 1990), ambiguity function algorithm (Remondi, 1991), fast ambiguity resolution approach (FARA) ( Frei and Beulter, 1990), Cholesky decomposition (Euler and Landau, 1992), OMEGA (Kim and Langley, 1999), LAMBDA (Teunissen, 1994(Teunissen, , 1995;;Teunissen and Verhagen, 2008), etc.Recently, the fast integer least-squares estimation for high-dimensional ambiguity resolution using lattice theory was proposed (Jazaeri et al., 2012(Jazaeri et al., , 2014)).Among them, some were based on dual frequency measurements.When applied for single frequency kinematic positioning, those methods require more epochs, resulting in large amount of calculation caused by higher-order matrix inversion operation and a long time to fix the ambiguity.Thus, those methods could not meet the requirement of real time kinematic applications.
To solve this problem, a new algorithm is proposed to implement fast integer ambiguity resolution of kinematic application.Firstly, SVD decomposition and transformation are applied to construct a left null space matrix to remove baseline coordinate parameters and separate ambiguity parameters from position parameters.Kalman filter is used to estimate only the ambiguity parameters that can be used to acquire real-time float solution of integer ambiguity.Then, diagonal entries of covariance matrix are sorted, and multi-time (inverse) paired Cholesky decomposition is applied for the decorrelation of ambiguity.After diagonal elements preprocessing and diagonal elements sorting according to the results of Cholesky decomposition, the efficiency of decomposition and decorrelation is improved.Finally, integer ambiguity is estimated by the integer search algorithm of LAMBDA method.Static and kinematic experiments prove the correctness and feasibility of the new algorithm.1 n  NZ is the unknown DD ambiguity parameters vector with n-1 dimensions, which is independent of the epoch.

Fast Calculation of Ambiguity Float Solution
is the coefficient matrix with n-1 dimensions, where 1 L  is the L1 carrier wavelength.i ε is the measurement noise vector.
In fast positioning, it is preferable that fewer epochs or even single epoch can implement positioning.
According to Equations 1, it has 3 rank deficiencies when calculating in single epoch, thus the LS method is unavailable.What is adopted in common is to increase the number of observation epochs, i.e., to increase equations.For m epochs, corresponding equations are
High-order matrix inversion is a problem (Liu et al., 2005;Liu et al., 2013) of huge computation in directly solving Equations 3. In order to obtain high precision ambiguity float solution, at least 200 to 400 epochs is needed, which results in very high-order matrix with about 600 to 1200 dimensions.Thus, it can not meet the requirement of real-time kinematic positioning.In fact, only the ambiguity float solution and its covariance are concerned.Therefore this paper applies SVD decomposition to coefficient matrix (Timothy, 2010).Baseline component correction vector is eliminated by constructing left null space matrix A L according to U matrix features.In this way, ambiguity parameters can be successfully separated from position parameters, with reduced matrix dimension.
Based on SVD decomposition, transformation steps are as follows (1) Solve the left null space matrix of coefficient matrix and A and i  denotes all non-zero singular values of matrix i A .
(2) Divide the matrix U as . What can be easily proved is the following formula:  Where V is the measurement white noise with mean zero and variance cov( , ) I .Since the state vector of integer ambiguity N is constant, time update process (prediction process) during Kalman filtering can be simplified as follows Kalman filter Measurement update process (calibration process) is given as Among Equations 8, k K is filtering gain; I is unit matrix.According to the initial value of state vector 0 N and the estimation error covariance ˆ0 N Q , the optimal estimation of ambiguity state vector and estimation error covariance at any time can be obtained.As indicated from the Filtering equations above, state parameters, such as position parameters and velocity parameters, are eliminated during the process of calculation.Therefore, the computation load is greatly reduced by avoiding of high-order matrix inversion, which contributes to ambiguity estimation in real-time kinematic applications.

Sorting and Multi-Time (Inverse) Paired Cholesky Decomposition for Decorrelation
Real-time float solution of ambiguity can be obtained based on SVD decomposition and Kalman filtering.Nevertheless, during a short period of observing time in actual kinematic positioning environment, double difference observation may lead to inferior geometric structure between the ground station and satellites, which in turn leads to strong correlation between DD measurements.This strong correlation extends multi-dimensional ellipsoid search space and integer ambiguity search results are far from expectation.To solve this problem, the covariance matrix of float solution needs to be decomposed and diagonalized as much as possible.In this way, the correlation between ambiguities of DD phase measurements can be reduced, which makes the search ellipsoid space closer to sphere, minimizing length of search interval and improving the efficiency in the search of integer ambiguity fixed solution.Concerning that Gauss transformation for decorrelation is not stable, and its calculation amount is double of Cholesky decomposition, this paper proposes a method of continuously implementing modified upper triangular Cholesky decomposition and lower triangular Cholesky decomposition to Bol.Ciênc.Geod., sec.Artigos, Curitiba, v. 21, n o 4, p.832-847, out-dez, 2015.realize continuous decorrelating (Xu, 2001;Zhou, 2011;Zhou and He, 2014).Before Cholesky decomposition, sorting the diagonal elements in ambiguity covariance matrix (ascending or descending) can reduce condition number of matrix and effectively improve the performance of matrix decomposition (Liu et al., 1999;Chen and Wang, 2002;Huang and Chen, 2010).Compared to this traditional method, this paper applies a method of diagonal elements preprocessing after Cholesky decomposition.New algorithm orders diagonal elements of the matrix according to values after decomposition, in contrast to existing method according to values before decomposition.Consequently, our new method is closer to the goal of arranging the larger diagonal element to the smaller row before decomposition for UDUT, and LDLT decomposition is on the contrary.
The sorting depends on the relative size of the diagonal candidate values   . Where ( , )   ii iM

S
is adjusting matrix, obtained by exchange of .
3) Let entries in upper triangular be integers and calculate the transformation matrix   , where ( , )   ii iM G is adjusting matrix, obtained by exchange of ith row and i M th row in unit matrix nn I  .In particular, ( , ) Modified LDL T decomposition is obtained by the previous step, which is Let entries in lower triangular be integers, calculate the transformation matrix   indicates that the correlation between the undetermined ambiguity is quite strong, and the process is back to step I) and go on until   1 1  L becomes unit matrix, where the decomposition stops (Details of the proof of the correctness of this decorrelation method are discussed in Appendix).
Conduct the modified integer UDU T decomposition and LDL T decomposition repeatedly.Assuming that iteration is executed for m times, and   1 1  L is transformed to be a unit matrix, the final transformation matrix is going to be Correspondingly, the covariance matrix of ambiguity after the transformation is Integer ambiguity float solution is Finally, integer ambiguity is estimated by the integer search algorithm of LAMBDA method.z is searched using Equations 16 to minimize the object function as the estimated ẑ

Then perform inverse transformation
The original space of ambiguity is therefore obtained.

Experimental Analysis
In order to verify the accuracy, correctness and effectiveness of proposed algorithm, static and kinematic experiments are carried out respectively.Actual collected receiver data is processed to calculate the navigation result.Then, the collected data is compared with the standard value to prove the reliability of this algorithm.

Static Positioning and Result Analysis
In order to validate the correctness and feasibility, static test was carried out.The data was acquired at ten thirty-one on May 16, 2013, on the top Air Force Engineering University Taoyuan campus laboratory building.Two NovAtel receivers were deployed with type OEM628 as one base station and another mobile station.All stations are installed GPS-702-GG GPS dual-frequency antenna.The recording sampling rate is 1 s.In this case, the baseline length is 7.812 m.GPS data were collected in half an hour tracking 7 satellites with 10° elevation mask: PRN6, PRN8, PRN11, PRN15, PRN17, PRN24, and PRN28.In order to reduce the errors' effect, PRN24 with maximal elevation was chosen as the reference satellite.Actual data was processed by proposed algorithm in simulation.After 120 epochs, integer ambiguity was estimated.Then, integer ambiguity was substituted back to the algorithm to compute the fixed solution, which is used to make a comparison with the real value.Figure 1 shows the length of baseline and its error.As can be seen from the figure baseline length error of proposed new algorithm is less than 1 cm, which indicates some good performance of accuracy.

Kinematic Positioning and result analysis
To further validate the applicability of the algorithm in a kinematic environment and a relative longer baseline situation, a vehicle test was carried out.

Verification scheme
Test conditions being the same with that of static test, dual mobile station are applied, due to the fact that dynamic changes of baseline with a much longer distance and accuracy of the algorithm can not be validated directly, as shown in Figure 2. Base station is placed in open area under the Lab building and two mobile stations are placed on the top of the vehicle.One is at the head and another one is at the rear, 2.115 m apart.After a short static period, the vehicle began to go around in circles, 200 meters away from the building.Original measurements data and epoch data (L1 frequency) were collected separately and processed by the proposed fast new algorithm one epoch to another.As is shown in Figure 2, the baseline between antenna 2 and antenna1 is computed and named as

Experiments and Results Analysis
Using the proposed algorithm, the integer ambiguity between mobile station 2 and the base station 1 is resolved, so is the integer ambiguity between mobile station 3 and the base station 1.Two receivers work simultaneously with 7 visible satellites (satellite cut-off angle sets to 10°).In order to reduce the errors' effect, the satellite with maximal elevation is chosen as the reference satellite.Finally, it can be combined to obtain six double differential ambiguities [N 21 , N 31 , N 41 , N 51 , N 61 , N 71 ].
The real double difference ambiguity (i.e., N=[3, 8, -1, 6 ,4, -6]) between station 2 and station 1 is used as the reference value.Then the effect of ambiguity float solution solved by this new algorithm using Kalman filter is analyzed.All the float solution of double difference ambiguities are shown in Figure 3.As can be seen from Figure 3, all the float ambiguities present a trend that approximately equal to the true value at the time of epoch 168, and the filter starts to become stabilize state.At the same time, different epochs of data are resolved by using the proposed approach and the conventional Kalman filter respectively.The results comparison of two methods in which parts of continuous epochs around epoch 168 (i.e., epoch 168 as the point) are selected is shown in Table 1.

Table 1: Results comparison of two methods
Table 1 shows that precision of ambiguity float solution calculated by the proposed method is better than conventional method.In addition, the average calculation time of both methods for Table 1 are respectively counted.Average calculation time of this new method is 1067 ms, while the conventional Kalman filtering method is 1703 ms.The reason is that the conventional Kalman filter needs to simultaneously estimate all the parameters such as position parameters, velocity parameters and ambiguity parameters, while in the new algorithm, SVD decomposition transform of the coordinate coefficient matrix is applied to construct the left null space matrix in order to eliminate the baseline coordinate vector parameters, thereby the Kalman filter equations can be established to estimate only the ambiguity parameters.This will not only improve the computing speed, but also improve the accuracy of float solution, which contributes to the fast resolution of ambiguity.
Similarly, the double difference ambiguity between station 3 and station 1 can be resolved.Therefore, integer ambiguities are fast estimated when epoch number is 168 and 177 respectively, the Kalman filter is becoming stable gradually and the integer ambiguity can be computed rapidly.The relative position between mobile station 2 and station 3 is computed by the validation scheme described above.The length and error of baseline are shown in Figure 4 As is shown in Figure 4, the computed position relationship between mobile station and base station is consistent with the reality.The errors of elevation and azimuth are stable in a small angle because of the circle motion of vehicle on the ground.During the test, with good signal quality, baseline error is limited within 2 cm, which indicates that ambiguity is correctly estimated.As the length of baseline, elevation and azimuth are all computed based on the relative position between mobile stations, so making the length of baseline as criterion is practical and scientific, which further verifies the correctness of the algorithm.Additionally, the performance of decorrelation is analyzed with the decorrelation coefficient r and spectral condition number e as standard to evaluate the performance of decorrelation (Huang and Chen, 2010).As is indicated in Figure 5(a), decorrelation coefficient rbefore is near 0 before decorrelation, which indicates strong correlation between ambiguities.The spectral condition number ebefore is large in Figure 5(b), which indicates a flat search space, also reflecting the strong correlation.After decorrelation, the decorrelation coefficient rafter becomes larger and approaches 1.At the same time, its spectral condition number eafter is greatly reduced and covariance matrix is closer to diagonal one.Correlation between different ambiguities is reduced successfully, which indicates the feasibility and good performance of decorrelation of the proposed algorithm.

Conclusions
In this paper, the static and kinematic tests and analysis show that the proposed algorithm for ambiguity fast resolution is feasible, the static baseline error is less than 1cm and kinematic baseline error is less than 2cm.These results verify the correctness and effectiveness of the algorithm.SVD decomposition is applied to construct the left null space of matrix to eliminate the baseline coordinate parameters which can separate the ambiguity parameters from the position parameters.Thus Kalman filter is used to estimate only the ambiguity parameters in the new algorithm, which greatly reduces the amount of computation.Computation speed are increased, which means its real-time capacity.Sorting and multiple (inverse) Cholesky decompositions are performed for ambiguity decorrelation, adopting method of diagonal entries pre-processing and adjusting the order of diagonal entries according to values by Cholesky decompositions.The effectiveness of matrix decomposition is well ensured and much smaller conditional number is obtained, thereby performance of decorrelation is improved, which contributes to ambiguity search efficiency and correctness.
As the new algorithm in this paper can implement fast solution of integer ambiguity, acquiring high accurate position quickly, it could be used for 'BDS-2' system in the near future.The algorithm may Bol.Ciênc.Geod., sec.Artigos, Curitiba, v. 21, n o 4, p.832-847, out-dez, 2015.
have a broad application prospects for airborne platform fast positioning and attitude determination and for precision approach and landing system.The proposed method also needs further improvement in some specific engineering implementation.
Based on SVD DecompositionSuppose that base station and mobile station observe n satellites synchronously, and each epoch can construct n-1 double difference carrier phase measurement equations.For the ith epoch, GPS double difference carrier phase linear observation equations are generally expressed as observation vector of double difference (DD) carrier phase at the ith epoch, which is the difference between observed value and the calculated one.ith epoch.i  X is the unknown parameter vector of 3D baseline.
formula above is rewritten as

L
on both sides; and the transformed equation is Considering the characteristics that integer ambiguity of each epoch is equal when there is no cycle slip, Kalman filtering state equation and measurement equation(Tomoji and Akio, 2010) are respectively written as Bol.Ciênc.Geod., sec.Artigos, Curitiba, v. 21, n o 4, p.832-847, out-dez, 2015.
computed by pre-compute formula, instead of relative size of the diagonal elements of the original covariance matrix.Adjusting matrix could only be determined after pre-computing the candidate elements vp of diagonal elements di, which ensures that the adjusted covariance matrix after Cholesky decomposition acquire the best decorrelation performance.After sorting based on the result of Cholesky decomposition, bigger diagonal entries of UDU T result is adjusted to a smaller I positions in diagonal line, and LDL T vice versa.The value of each diagonal element after adjustment is much closer to one another, and the validity of matrix decomposition is improved, so is the decorrelation.Then, the efficiency and quality of ambiguity discrete search is improved.The steps of sorting and multiple (inverse) paired Cholesky decomposition are as follows I) Modified integer UDU T decomposition candidate entry in unit upper triangular matrix U and p v denotes the candidate in diagonal matrix D. elements in modified Cholesky decomposition U and D respectively.c) Adjust variance-covariance matrix according to the index number constructed from the last step: 1 UDU T decomposition is obtained by the previous step, where [] .Geod., sec.Artigos, Curitiba, v. 21, n o 4, p.832-847, out-dez, 2015. 1 the following steps to 0  J row by row.a) Pre-compute candidate elements In Equations 11, pj f is the candidate element in unit lower triangular matrix L and p v denotes the candidate element in diagonal matrix D. elements in modified Cholesky decomposition L and D respectively.c) Adjust variance-covariance matrix according to the index number from the last step: 1 the covariance matrix after modified LDL T transformation is III) Check whether   1 1  L is an unit matrix.If so, then calculation completes; otherwise, it Bol.Ciênc.Geod., sec.Artigos, Curitiba, v. 21, n o 4, p.832-847, out-dez, 2015.

Figure 1 :
Figure 1: Baseline Length and its error.

1d
vector, so was the baseline between antenna 3 and antenna 1, named as 2 d vector.The baseline between antenna 2 and antenna 3 is fixed and named as d .It is easy to know .Then it is compared with fixed length of baseline d to verify the accuracy of relative position.

Figure 3 :
Figure 3: All the float solution of double difference ambiguities.
(a) and the elevation and azimuth of baseline are shown in Figure 4(b).

Figure 4 :
Figure 4: Relative position of the mobile station 2 and station 3.
Figure 5(a) and Figure 5(b) shows the curves of decorrelation and spectral condition number before and after decorrelation respectively.

Figure 5 :
Figure 5: The r and e of ambiguity covariance matrix before and after decorrelation.