Iterative Least Squares Method for Global Positioning System

. The efﬁcient implementation of positioning algorithms is investigated for Global Positioning System (GPS). In order to do the positioning, the pseudoranges between the receiver and the satellites are required. The most commonly used algorithm for position computation from pseudoranges is non-linear Least Squares (LS) method. Linearization is done to convert the non-linear system of equations into an iterative procedure, which requires the solution of a linear system of equations in each iteration, i.e. linear LS method is applied iteratively. CORDIC-based approximate rotations are used while computing the QR decomposition for solving the LS problem in each iteration. By choosing accuracy of the approximation, e.g. with a chosen number of optimal CORDIC angles per rotation, the LS computation can be simpliﬁed. The accuracy of the positioning results is compared for various numbers of required iterations and various approximation accuracies using real GPS data. The results show that very coarse approximations are sufﬁcient for reasonable positioning accuracy. Therefore, the presented method reduces the computational complexity signiﬁcantly and is highly suited for hardware implementation.

The result of the free availability of satellite positioning parameters has led to wide adoption of the GPS systems.Building GPS devices in commercially available cell phones has been achieved by mobile device providers, such that the number of cell phones equipped with GPS functionality has rapidly grown in the last few years.
In order to do the positioning, an initial set of pseudoranges between the receiver and the satellites is needed.Nonlinear LS is the most common method to determine the receiver's position from the pseudoranges.Usually, linearization is done to convert the non-linear problem into an iterative algorithm, which requires the solution of an overdetermined system of linear equations in each iteration step itr, i.e. linear LS method is applied in each iteration step itr.For solving the linear LS problems in each iteration step an iterative version of the QR decomposition (QRD) [Gtze (1994)] is applied in this paper.Instead of annihilating the lower diagonal elements during the QRD, CORDIC-based approximate rotations are used.By choosing the accuracy of the approximation, e.g. by choosing itg optimal CORDIC angles per rotation, the LS computation can be simplified.However, we only obtain an approximate solution to the LS problem, whose accuracy depends on itg.The accuracy of the positioning results of GPS method is compared for varying numbers of iterations itr of the positioning algorithms and varying numbers of iterations itg of the iterative QRD using real GPS data.The results show that very coarse approximations (small itg) are sufficient for obtaining a reasonable position estimate.Therefore the presented methods reduce the computational complexity and the required power consumption significantly.
GPS positioning is introduced in Sec. 2 resulting in an algorithm, which requires the solution of LS problems in each iteration itr (itr = 1, 2, • • • , itr max ).For solving these LS problems an iteration version of the QRD is presented in Sec. 3 using CORDIC-based approximate rotations, where itg denotes the approximation accuracy (itg = 1, 2, • • • , itg max ).The trade-off between itr (iteration of the positioning method) and itg (iteration of the iterative QRD) is investigated in Sec. 4, where experimental results are given using real GPS data.The paper finishes with a conclusion and an outlook to the future work in Sec. 5.

GPS Positioning
The whole GPS positioning procedure includes three tasks: acquisition, tracking and positioning [Borre (2007)].The acquisition tries to find satellites and to get their positions.It gives rough estimates of signal parameters.Tracking keeps track of these parameters as the signal properties change over time.After tracking, the navigation data can be extracted and pseudoranges (measured distance from satellites to GPS receiver) can be computed.The final task of the receiver is to compute the user position.
The GPS satellites' arrangement ensures that every point on our planet is in contact with at least six satellites at all times.Each satellite k continuously broadcasts a digital radio signal that includes its position (X k , Y k , Z k ) and its time t k .On board atomic clocks ensure an accurate time to a billionth of a second.The radio signal of satellite spreads with c = 3 * 10 8 in universe, the velocity of light in vacuum.GPS receivers measure the time delay τ k of the signal from each satellite k to the receiver, so τ k = t − t k , where t is time of receiver.The measurement of t in the receiver is not very accurate (as compared to the satellite time t k ).Furthermore the speed of radio signal from the satellites cuts down because of ionosphere and troposphere.Signal propagation duration from satellites to the receiver is longer than expected.Therefore the measured distance from the satellite to the receiver P k , measured by P k = τ k • c, is a rough distance estimate called 'Pseudorange' (see Fig. 1).The receiver simultaneously collects these measurements from at least four satellites and processes them to solve for position and time measurement error.

Observation Equation
The most commonly used algorithm for position computation from pseudoranges is based on the LS method.
This method is used to find the receiver position from pseudoranges to four or more satellites.
The basic observation equation for the pseudorange P k is ρ k is the geometrical range between satellite k and receiver, which can be computed as: where (X, Y, Z) is the position of receiver (see Fig. 2).dt denotes the receiver clock offset and dt k is the satellite clock offset.From the ephemerids, which also include information on the satellite clock offset dt k , the position of the satellite (X k , Y k , Z k ) can be computed.T k is the tropospheric error and k is the ionospheric error.These two errors are computed from a priori models, whose coefficients are part of the broadcast ephemerids.e k is the observation error of the pseudorange.Therefore, Eq. 1 contains four unknowns X, Y, Z and dt.The error terms are minimized by using the LS method.Eq. 2 is nonlinear with respect to the receiver position (X, Y, Z), so the equation is linearized before using the LS method.The nonlinear term in Eq. 2: (3) is linearized.Starting from an initial position for the receiver (X 1 , Y 1 , Z 1 ), the position estimate is improved iteratively.The center of the Earth (0, 0, 0) can be chosen as the initial point, if no a-priori-information (as e.g. from AGPS) is available.Let itr be the number of the iterations (itr = 1, 2, • • • , itr max ).The increments ∆X itr , ∆Y itr , ∆Z itr update the receiver coordinates as follows: Eq. 5 includes only first-order terms, and the partial derivatives are ) 2 be the range computed from the satellite (X k , Y k , Z k ) to the approximate receiver position (X itr , Y itr , Z itr ), so the first-order linearized observation equation becomes where dt itr is the estimated clock error at the receiver.

Applying Least-Squares Method
Rearranging Eq. 6, it can be rewritten as A unique solution can not be found from a single equation.Therefore m ≥ 4 satellites are required to form a system of linear equations (usually m ≥ 6 satellites are available).Let Then we obtain the LS problem: where and x itr = [∆X itr ∆Y itr ∆Z itr cdt itr ] T .The solution ∆X itr , ∆Y itr , ∆Z itr has to be added to the approximate receiver position to get the next approximate position as in Eq. 4. These iterations continue until the solution ∆X itr , ∆Y itr , ∆Z itr is at meter level.
For solving Eq. 8 we need to find xitr which minimizes the length of the error vector êitr = b itr − A itr xitr with ||e itr || 2 = (b itr − A itr x itr ) T (b itr − A itr x itr ) the sum of squares of the m separate errors.Minimizing this quadratic gives the normal equations 9) and the error vector is There are various algorithms for solving LS problems (Eq.8) [Golub (1996)].In the subsequent section, we use the QRD by Givens rotations.

Iterative Solution of LS
In each iteration of the positioning algorithm a LS problem must be computed.Since the pseudoranges are subject to measurement errors and the convergence (number of required iterations) of the algorithm depends on the accuracy of these LS solutions, it is worthwhile to investigate the use of an iterative LS solver and the trade-off between the number of iterations of the positioning method and the number of iterations of the iterative LS solver.
The iterative version of the QR decomposition (QRD) presented in [Gtze (1994)] is used for the iterative solution of Eq. 8, since it is well suited for hardware implementation, suitable for the adaption to measurement errors (pseudoranges), and yields a solution vector which converges linearly to the exact solution.Here we will briefly review this iterative QRD.
The QR decomposition of a m × n matrix A = Q • R can be computed by applying a sequence of Givens rotations G(φ ij ) (φ ij = arctan(a ij /a jj )) to the matrix such that the matrix entries below the diagonal of A are annihilated, i.e. generate a ij = 0 for

Yuheng He et al.: Iterative Least Squares Method for Global Positioning System
The resulting matrix is upper triangular and denoted by R. The product of all required Givens rotations forms the or- and the solution w can be computed by back substitution.
One iteration of the iterative version of the QRD works exactly like the QRD but instead of using exact rotations G(φ ij ) that annihilate a ij (a ij = 0) CORDIC-based approximate rotations are used resulting in The CORDIC-based approximate Givens rotations are computed by determining the shift value , ∈ {0, 1, 2, . . ., b} (b wordlength) corresponding to the CORDIC angle φ ij ( ) = arctan 2 − which is closest to the exact rotation angle φ ij .Instead of annihilating the matrix elements during the course of the QRD an approximate rotation will only reduce the matrix elements by the maximal factor possible with itg specific CORDIC angles φ ij ( ).This CORDIC-based approximate rotation is applied to the QRD for solving LS problems.Since the matrix elements below the diagonal are no longer annihilated the QRD procedure using CORDIC-based approximate rotations must iteratively be applied until the matrix is ultimately upper triangular.One obtains iterative versions of the QRD distinguished by itg, which determines the accuracy of the approximation of the rotations.
Defining the lower diagonal quantity for iteration itg: the above algorithm guarantees When we apply this method to the LS problems, which must be solved for positioning in Eq. 8, very few steps, i.e. itg b, of the CORDIC-based approximate rotation are sufficient to obtain similar results as using exact rotations.This yields a significant reduction in computational complexity by itg/b (we will show in the following section for real GPS data that even itg = 1 gives reasonable results).Furthermore, this method only requiring shift and add operations is very well suited for hardware implementation.

GPS raw data collection
For GPS positioning, a SiGe GN3S Sampler v2 is used to capture the raw GPS data, which are low level signal data (raw intermediate frequency samples) being delivered by the GPS satellite network and processed by the SiGe radio front end [ GN3S (2009)].Each GPS satellite continuously broadcasts a navigation message at a rate of 50 bits per second.
The obtained raw GPS data at each measurement point includes information of m = 7 satellites with matched pseudoranges.78 raw GPS data records are gained.Afterwards position calculation is done for GPS method for various numbers of required iterations itr (see Sec. 2) and various approximation accuracies itg (see Sec. 3).

Experimental results
The mean values of the positioning results for the 78 measurements are presented in Tab. 1.The calculated positions are compared with the exact positions (known coordinates at the measure points from land surveying office in Bochum Germany) and the accuracy of the positioning results of GPS method for varying itr and itg is compared.The number of required iterations itr is at least two.The second column itg in Tab. 1 is the approximation accuracy in QRD.Fig. 3 and Fig. 4 show the position estimates of GPS positioning.In Tab. 1 the 3rd column GP S − itr is the accuracy of GPS positioning with exact LS method and the 4th column GP S − itg is the accuracy of GPS positioning using QRD with CORDIC-based approximate rotations for solving LS problems.The results show that with the increasing number of iterations in QRD, the accuracy of GP S − itg also increases.If itg ≥ 8, the accuracy of the exact LS method is achieved.Especially, if itr is big enough, e.g.itr ≥ 3,  only one iteration in QRD itg = 1 is required to compute the position with reasonable positioning accuracy, which means very coarse approximations are sufficient.The position estimates of 78 measurements are shown in Fig. 3. GPS position with exact LS (black circles), its mean value (black cross), GPS position using iterative QRD with itg = 1 (blue stars), its mean value (blue cross), GPS position using iterative QRD with itg = 3 (green squares), its mean value (green cross) and the exact position (red cross) are shown in the figure.All the position results are considered as accuracy of position, i.e. the position estimates subtracted by the exact position (so the exact position is set to (0, 0)).
Fig. 4 is an enlarged part of Fig. 3.It is obvious to notice that GPS position using QRD with itg = 3 (the green squares) are more closer to GPS position with exact LS (black circles) than GPS position using QRD with itg = 1 (blue stars).If the number of iterations in QRD increases, i.e itg increases, the position results will become more and more similar to the results of exact LS.The mean value of itg = 3 is almost the same as the mean value of exact LS (green cross is almost at the same position as black cross in Fig. 4).

Conclusion
An iterative LS approach and an iterative version of the QRD using CORDIC-based approximate rotation are applied to position computation.The accuracy of the positioning results of GPS method is compared for various numbers of required iterations itr and various approximation accuracies itg of CORDIC-based approximate rotations by using real GPS data.It is shown that a significant reduction concerning computational complexity and hardware requirements can be obtained.Furthermore, this method only requiring shift and add operations is very well suited for hardware implementation.
The presented method is very efficient for the implementation of the standard triangulation method based on nonlinear LS.In future work we apply this idea to recursive computations of the position estimates using Orthogonal DGPS [Chang (2003)] and Kalman Filter based recursive GPS algorithms [Grewall (2001)].In the first case no iterative LS method is required for positioning which leads to further reduction of the computational effort and the power consumption.In the second case using the square root version of the Kalman filter [Merwe (2001)] QRD can be applied and the required approximation accuracy (number of itg) can be investigated to obtain desired position accuracy.

Table 1 .
Mean value of position-accuracy results (in meter) by GPS with exact LS method and CORDIC-based approximate rotations