Simple Orbit Determination Using GPS Based on a Least-Squares Algorithm Employing Sequential Givens Rotations

A low-cost computer procedure to determine the orbit of an artificial satellite by using short arc data from an onboard GPS receiver is proposed. Pseudoranges are used as measurements to estimate the orbit via recursive least squares method. The algorithm applies orthogonal Givens rotations for solving recursive and sequential orbit determination problems. To assess the procedure, it was applied to the TOPEX/POSEIDON satellite for data batches of one orbital period (approximately two hours), and force modelling, due to the full JGM-2 gravity field model, was considered. When compared with the reference Precision Orbit Ephemeris (POE) of JPL/NASA, the results have indicated that precision better than 9 m is easily obtained, even when short batches of data are used.


Introduction
The Global Positioning System (GPS) provides a powerful and quick means to compute orbits for low Earth artificial satellites.Theoretically, four GPS satellites, simultaneously tracked, are enough for geometrical positioning of an artificial satellite carrying an onboard GPS receiver.Dynamical orbit determination is a nonlinear problem where the perturbing factors are not easily modelled.The GPS satellites transmit signals such that accurate measurements of distances are performed based on the comparison between received signals and template signals generated by the receiver [1].Through a GPS receiver onboard of an artificial satellite it is possible to obtain such measurements (pseudoranges) that can be processed to estimate a state vector (e.g., position and velocity vectors and parameters referred to the receiver clock bias).Using our knowledge about the dynamics of the motion and, at the same time, assuming statistics for the measurement errors, the state vector can be computed based on a set of observations.To this, the differences between the observed and modelled observations are minimised according to the least squares criterion [2][3][4].
The aim of this work is to determine the orbit of an artificial satellite carrying a GPS receiver, by a recursive least squares method based on L1-code GPS satellite-to-satellite tracking observations, and to make use of Givens rotations [5] in order to solve the problem in a recursive way.The Givens rotations algorithm is as stable as other orthogonalization algorithms (such as Cholesky decomposition, Gram-Schmidt, or Householder) and allows the processing for each observation epoch, avoiding the storage of large matrices.Advantages of the use of Givens rotations were already point out in the seventies of the last century in an improvement of a computational form of the discrete Kalman filter [6].Givens rotations have been successively used at Center for Space Research, University of Texas at Austin, Austin, Tex, USA [7].
Thus, this procedure for near real time positioning has been considered to be used in the next two Chinese-Brazilian remote sensing satellites CBERS-3 and CBERS-4 scheduled to fly in 2008-2010, according to a Brazil-China protocol agreement.These missions carry a main CCD camera payload with 20 m resolution and, therefore, a simple orbit determination scheme with 10 m accuracy is a desirable requirement for the image processing users.Due to onboard memory limitations it is likely that only a short batch of GPS data (one revolution) around the image scene will be available.
Initially, a simple dynamical model (geopotential perturbations only) was used to preliminary assess the procedure, by using the TOPEX/POSEIDON satellite as a test bed.However, the final purpose is to evolve to applications of the procedure considering real and suitable models for the perturbations well suited for the orbital geometry and physical characteristics of the CBERS mission satellites [8].

Dynamical model
The problem of dynamical orbit determination is essentially nonlinear.The orbital motion is described in an inertial frame by a set of ordinary differential equations: where r = (x, y,z) is the position vector, μ is the gravity parameter, and P stands for modelled perturbations.Formally, we include an additional equation ḃ = 0 with b = [b 0 ,b 1 ,b 2 ] solve-for parameters referred to the receiver clock bias, (bias, drift, and drift rate).Thus, the state vector to be estimated is defined by x ≡ r v b t with v being the velocity vector.The state transition matrix which relates the state between times t k and t k+1 can be computed by Φ(t,t k ) = F(x,t)Φ(t,t k ) with initial condition Φ(t k ,t k ) = I [3,9,10].The Jacobian matrix F(x,t) contains partial derivatives of the differential equations of motion.In general, the transition matrix is numerically integrated together with the orbit.In our case, where only the geopotential perturbation is modelled, the matrix F can be written as where f is the acceleration acting on the satellite, 0 3×3 is a 3 × 3 matrix of zeros, I 3×3 is a 3 × 3 identity matrix and, A 3×3 is the 3 × 3 matrix of gravity gradient 3 × 3 given by A 3×3 = ∂ v/∂r.In this work, the modelled force includes the geopotential, taking into account the spherical harmonic coefficients up to 50th order and degree of the JGM-2 model [11].The geopotential acceleration and its gradient matrix A 3×3 were computed using basically Pine's universal recursive formulation [12], however we use the improved numerical recursions as stated in Lundberg and Schutz [13].For assessing the proposed approach, we decided to consider only the geopotential perturbation, applied to TOPEX/Poseidon satellite for which results are broadly available for comparison purposes.Considering that the aimed application would be the CBERS satellite, the final procedure, as developed for the CBERS mission, will include also other relevant perturbations such as radiation pressure forces, lunisolar attraction, and atmospheric drag.

Measurement model
The general nonlinear equation representing the scalar model of observations at epoch. is given by y k = h k (x k ,t) + v k , where y k is the vector of m observations, h k (x k ) is the mdimensional nonlinear vector function of the state vector x k , v k is the m-dimensional vector of observation errors modelled as white noise.Pseudoranges are measurements of the distance between the GPS satellites and the receiver's antenna, referring to the epochs of emission and reception of the signals.Pseudorange is a typical type of measurement in an orbit determination process using GPS, and can be written as [1] , where P i is the pseudorange measured by the user with respect to the ith GPS satellite, ρ i is the geometric distance, c is the speed of light, dt is the user clock offset, dT i is the ith GPS satellite clock offset, D ion are ionospheric delays (not considered in this work), D trop are tropospheric delays (not applicable), and v denotes random measurement noise.Tests computing the first-order ionosphere effect using P 1 and P 2 measurements, were performed.They show that, given the aimed for precision in this work (order of 10 m), ionospheric delays can be neglected.The user clock offset is modelled as is a vector of parameters to be estimated representing the clock bias, bias rate, and rate of bias rate.To estimate the state vector parameters using the least squares algorithm (e.g., [14]), the measurement equation is linearized and the sensitivity matrix of partial derivatives with respect to the state vector is derived.

Recursive least squares method based on givens rotations
The numerical stability of a state vector estimator based on the least square method should be assured.To validate the method for orbit determination using a recursive approach, namely the recursive least squares method, one is tempted to use the Kalman form [2] in straightforward way.Extensive analyses of such algorithms have shown that they are unstable numerically and sensitive to error accumulation [2,3].On the other hand, the nonrecursive solution of normal equations requires a matrix inversion, another source of numerical errors, which should be avoided.In order to overcome such problems, several methods can be found in the literature using orthogonal transformations.They yield best numerical performance with respect to problems due to error propagation or uncertainty in the information.The main goal of applying orthogonal transformations to matrices and vectors in the least squares algorithms is the substitution of the conventional algebra based on matrix inversions by a method that is more robust and less prone to numerical errors.Among the several existing orthogonal transformations (e.g., Givens, Gram-Schmidt, Householder), the Givens rotations algorithm is the most adequate transformation to selectively annihilate elements of a matrix, making it easy to implement the least squares method in a recursive way.For a batch least squares method, Gram-Schmidt or Householder transformations are in general computationally less costly than Givens rotations.However, from a recursive formulation point of view, Givens rotations algorithm allows us to implement the processing of one measurement at a time, avoiding the need of large matrices for storage (see, e.g., [6,9,15]).Indeed, the direct solution of the normal equations can be quite sensitive to small errors in a sensitivity matrix H that are inevitable when forming the product (H T H), with a limited numerical machine accuracy.In order to avoid the normal equations of the form (H T H)x = H T y, the minimization of the least squares cost function using Q-R factorisation [6,16] is recommended: where the matrix H is factored into an orthogonal matrix Q and an upper triangular matrix R. Methods for performing the Q-R factorisation have been proposed [5,6,16,17], involving orthogonal transformations that successively annihilate the sub-diagonal elements of H. Nevertheless we concentrated our approach on the Givens rotations method for the reasons stated above.
The complete transformation can be written as the product of the sequence of orthogonal rotations required to transform Matrix H into the form of (4.1): denoting the sequence of orthogonal rotations required to put matrix H into the form of (4.1).Orthogonal transformations of matrices play a considerable role in the numerical computations of the least squares problems [18].In fact, the Euclidian norm of a vector does not change, the same accuracy is obtained with single computer floating point arithmetic that otherwise would require double precision and the problem is solved in a numerically robust fashion.
Several programming tests were carried out before choosing the Givens rotation algorithm for implementation in the recursive least squares method.The Kalman type algorithm led to truncation errors for long batches of data.The Gram-Schmidt and Householder method were not as well suited for recursive implementation.Thus the Givens algorithm was finally selected.The final combination of a simple orbit model and robust numerical properties of the Givens orthogonal transformation resulted in a very compact computer program with small requirements of storage, paving the way for further improvements in terms of more complex models for orbit and GPS measurements.

Results
Actual TOPEX/POSEIDON (T/P) satellite flight data are used to validate the proposed method.The T/P satellite was launched on August 10, 1992, and it orbits the Earth at an altitude of 1336 km in an orbit with an inclination of 66 • , with near-zero eccentricity, and an orbital period of about 1.87 hours.It has a GPS receiver onboard as experimental equipment to verify several precision methods for orbit determination.The receiver can track up to 6 GPS satellites simultaneously on two frequencies if antispoofing is inactive.RINEX format for T/P observation data, GPS group data, and navigation messages are provided by International GNSS Service (IGS), NASA, DC, USA [19], Goddard Space Flight Center (GSFC), NASA, DC, USA, and Jet Propulsion Laboratory (JPL), NASA, DC, USA.Here, the following files were used [20].
(1) T/P observation files: Pseudorange codes measured at two frequencies in GPS time steps of 10 seconds, and made available by the GPS data processing facility of JPL in RINEX format.(2) Files with the Precise Orbit Ephemeris (POE): generated by JPL in one minute UTC time steps in inertial true of date coordinates.(3) GPS navigation message files in RINEX format made available by the Crustal Dynamics Data Information System (CDDIS) of the GSFC.In this work, the estimated position and velocities were compared with the T/P Precise Orbit Ephemeris (POE) generated by JPL which provides position estimates with an estimated precision of 15 cm or better.The test conditions for the investigated problem are the following: (1) actual topex/poseidon pseudorange measurement dataset collected by the onboard GPS receiver on November 18 1993; (2) geopotential perturbations taking into account the spherical harmonic coefficients up to 50th order and degree of JGM-2 model; (3) GPS antenna offset neglected, as its effects is about 5 m; (4) pseudorange measurements at frequency L 1 (Code); (5) recursive least squares method with Givens rotations algorithm; (6) short dataset of 2 hours (about 1 orbital period).The estimates obtained by our method were compared with the POE reference and produced position error of better than 9 m, starting with the JPL/POE values at first epoch, but clock parameters at the 100 m error level.residuals are clearly unbiased with respect to the apriori residuals.These are the plots of pseudorange data residuals (observed minus computed) for different satellites, seen at different times around the orbit.The apriori covariance matrix was set to 1000 m and 10 m/s for each component of the position and velocity vectors, respectively, and 5 m standard deviation was adopted for the measurement error.

Conclusions
The aim of this work was to assess a recursive least squares method based on Givens rotations to determine the orbit of an artificial satellite using an onboard GPS receiver.A short batch of a two hours pseudorange measurement dataset (about one orbital period), collected by the GPS receiver onboard of the satellite, was used.The modelled forces included the geopotential taking into account the spherical harmonic coefficients up to 50th order and degree of the JGM-2 model.Results of this work, using T/P dataset from November 18, 1993, compared against the postprocessed GPS ephemeris POE/JPL, has demonstrated an accuracy better than 9 m, without requiring much computational effort.In accordance with the obtained results, we can conclude that the recursive least squares procedure implemented with Givens rotations algorithm is a simple, reliable, and numerically stable method for orbit determination by using GPS measurements.The attained precision can, of course, be improved, and the authors are working on it, considering other relevant perturbing effects [21].Finally, it is important to point out that for the target application, the CBERS-3 and 4 satellites, the perturbations due to atmospheric drag, solar radiation pressure, lunisolar gravity field and tides, as well as measurement errors (e.g., ionosphere) must be taken into account for precision applications.