Nonlinear analysis of position estimate in global navigation satellite systems

The paper presents a nonlinear analysis of position estimation based on a global navigation satellite system. A classical problem formulation and iterative solution that results in the weighted least squares estimate of the receiver state are assumed. The analysis employs the Taylor's theorem to express the nonlinear measurement model using the first order Taylor polynomial at the state estimate and the Lagrange form of the remainder. A sensitivity analysis of the Jacobian matrix pseudoinverse is performed, and an upper bound on the size of the Lagrange remainder is derived using eigenstructure of the Hessian matrix. The results obtained show that both the sensitivity of the pseudoinverse and the size of the quadratic term are not significant, and thus the linear approximation commonly used to derive stochastic properties of the state estimate is reasonable. Although this result has been experimentally confirmed by numerous successful applications, this analysis can serve as a more rigorous basis when the design procedures for a safety critical system have to be satisfied.


Introduction
Global navigation satellite systems (GNSS) have been continuously used in new kinds of applications that require positioning, navigation, or timing services [1,2]. Besides noncritical (car navigation, geocaching, etc.), there are safety critical applications (aviation navigation, train control systems, etc.) and liability critical applications (road tolling, traffic law enforcement, etc.) that impose increased requirements on accuracy, integrity, continuity, and availability of position, velocity and time information [3,4]. To make GNSS available for such applications, development has been carried out in several directions. The space and control segments of the fully operational Global Positioning System (GPS) are being gradually modernized, the new European GNSS Galileo is being developed with such applications in mind from the very beginning, satellite-based augmentation systems are being put in use, and GNSS receivers with advanced capabilities such as Receiver Autonomous Integrity Monitoring (RAIM) are being designed.
The development of a GNSS receiver with advanced capabilities to be used in a safety critical application has to comply with international and sector specific standards on the functional safety of systems [5,6]. Assuming that raw data like pseudoranges and ephemeris are available at sampling steps, the state, i.e., position and clock bias, of the GNSS receiver can be estimated using two different approaches. In the first approach, a static estimator that uses only data from the current sampling step is designed based on a measurement model. In the second approach, global or local nonlinear filtering techniques [7] are employed to design an estimator that uses a dynamical model in addition to the measurement model and data from the preceding sampling steps as well. Although the latter approach usually results into an estimator of a higher quality, it is more challenging to comply with safety standards as the dynamical model is application dependent and undetected faults may propagate into state estimates from the preceding sampling steps. Therefore, a static estimator will be considered in this paper.
The evaluation of state estimate quality and the design of a RAIM algorithm are almost exclusively performed using a linearized measurement model [8,9]. Nevertheless, a nonlinear analysis of the state estimation error based on two second order Taylor polynomials carried out in [10] showed that unlike the conclusions derived from the linearized model, the state estimate is biased and has a larger variance. It was also noted that the impact on a RAIM algorithm is difficult to evaluate in general because of combined effects. Although approximate expressions for the mean value and covariance matrix of the state estimation error were derived, the nonlinear analysis was not straightforward as it required several approximations to be made. Moreover, the analysis did not provide qualitative evaluation of nonlinear effects.
The aim of this paper is to present a more straightforward nonlinear analysis that utilizes Taylor's theorem with the first order Taylor polynomial and the Lagrange form of the remainder. A sensitivity of the Jacobian matrix pseudoinverse on the linearization point is examined and the eigenvalue analysis is employed to derive the lower and upper bounds on the linearization error. These results are subsequently used to quantify the nonlinear effects on the state estimate and the estimation error.
The paper is structured as follows. The problem of state estimation is presented in Section 2. Section 3 deals with the design of a state estimator that uses a nonlinear model and available measurements to compute the state of the receiver. Properties of the position estimator are analyzed in Section 4. First, the usual results based on linearized model are mentioned and then analysis utilizing the nonlinear model is presented. In Section 5, the results of the nonlinear analysis are demonstrated through a numerical example. Finally, Section 6 provides concluding remarks.

Problem formulation
The problem of positioning is considered in a common coordinate system (e.g. Earth-Centered-Earth-Fixed coordinate system) and only data from one time step are used. In this section, a nonlinear measurement model and a static state estimator are presented.

Nonlinear model of measurement
It is assumed that the measurement of pseudoranges in a GNSS can be described by the following static nonlinear measurement model where y ∈ R ny is a vector of measured pseudoranges to visible satellites, h : R nx → R ny is a known nonlinear vector-valued function of several variables, x 0 ∈ R nx is an unknown state of the receiver, f ∈ R ny represents unknown additive faults in the pseudoranges, and v ∈ R ny is the measurement noise that has the Gaussian distribution with the zero mean and known positive definite covariance matrix Σ. The state x 0 = [r T 0 , b 0 ] T consists of the unknown receiver position r 0 ∈ R 3 and unknown offset b 0 = c · dt 0 ∈ R caused by the receiver clock offset dt 0 ∈ R with respect to the global GNSS time. Note that c denotes the speed of light as defined in the GNSS. The individual functions of the nonlinear vector function h have the following form where · 2 denotes the L 2 -norm of a vector, x T = [r T , b] = [r 1 , r 2 , r 3 , b] is a point in R nx , s j = [s j,1 , s j,2 , s j,3 ] T ∈ R 3 denotes the known position of the j-th satellite, and dtsv j ∈ R denotes the known j-th satellite clock offset correction term. The measurement noise v represents a common level of uncertainty in the measurement of the pseudoranges whereas the faults f represent significant deviations that can be caused by unusual conditions in the ionosphere or troposphere, multipaths, errors in the navigation message, and others. In contrast to the measurement noise, the fault is not a random variable and no detailed information about its behavior is assumed.

State estimator
The objective is to design a state estimator that has the following form wherex ∈ R nx is an estimate of the unknown state x 0 and σ : R ny → R nx is a function that represents the state estimator.

Design criteria
The state estimator (3) is designed using the following criterion. The quality of the state estimator is evaluated through the weighted least squares criterion of the true and predicted measurements. Thus, the state estimate is computed by solving the partial optimization problem where the weighted least squares criterion V : R nx+ny → R is a function defined as and the equation error is e(x, y) = y − h(x).

Design of state estimator
The state estimator is obtained by solving the partial optimization problem defined in (4). For a given measurement y, the local minimum of the nonlinear least squares criterion V can occur at a stationary point where any of the first order partial derivatives of V does not exist or the condition is satisfied [11]. By taking the first order partial derivatives of the criterion V , the necessary condition (6) shows that the state estimatex has to satisfy the following system of nonlinear equations where H : R nx → R ny×nx is the Jacobian matrix of the nonlinear function h and its j-th row has the following form A geometric interpretation of the necessary condition is that the transformed equation error Σ −1 e(x, y) is perpendicular to the tangent hyperplane defined at the stationary pointx. Note that the domain of the Jacobian matrix is the whole R nx except for the set X = ∪ ny j=1 X j , where the set X j is defined as Since the least squares criterion V is a nonlinear function of x and the necessary condition for the stationary point is represented by the system of nonlinear equations (7), the state estimator cannot be expressed in the closed form. Instead, the state estimatex is computed using an iterative optimization method for a given measurement y. A frequently used iterative optimization method in GNSS is the Gauss-Newton (GN) method [1]. In its basic version, the GN method starts with an initial estimatex (0) and computes subsequent estimatesx (i) using the following update relation where i is the iteration index and the matrix valued function H † : R nx → R nx×ny defined as is the weighted pseudoinverse of the Jacobian matrix H (x). A state estimate from the previous time step or a fixed nonzero vector can be used as the initial estimatex (0) . The usual stopping conditions for the GN method are the maximum number of iterations i max and a limit GN on the L 2 -norm of the difference between two consecutive state estimates [12] i ≤ i max , In general, the basic version of the GN method may fail to converge. Even if the method converges, the attained minimum can be a local one instead of the global minimum because the nonlinear least squares criterion (5) is not necessarily a convex function. As the iterations of the GN method fail to converge to the point where the nonlinear least squares criterion V is not differentiable, it can be assumed that if the GN method converges, the resulting state estimatex is a stationary point, particularly a local minimum, of the criterion V and the necessary condition (7) is fulfilled up to an additive error that is proportional to GN . For the GN method to work, it is assumed that the Jacobian matrix has the full column rank at each estimatex (i) . Otherwise the GN method is terminated and the state estimate is not computed.

Analysis of state estimator properties
The previous section presented the state estimator that was designed by utilizing the nonlinear model. As pointed in [10], the properties of state estimate and a RAIM algorithm are almost exclusively derived from the linearized model despite that the system is nonlinear. To perform a nonlinear analysis, two second order Taylor polynomials of h, one defined at the true receiver state and the other one defined at its estimate, were used in [10]. It was shown that the linearization error can be regarded as an additive noise that makes the state estimate biased and its covariance matrix larger. However, the analysis did not include quantitative evaluation of the nonlinearity. In this section, a nonlinear analysis that uses a slightly different approach and allows the impact of the nonlinearity to be quantified is presented. The Taylor's theorem is employed to obtain the first order Taylor series with the Lagrange form of the remainder suitable for deriving bounds on the quadratic term of the Taylor series.
where Q j : R nx → R nx×nx is the Hessian matrix of the real-valued function h j , ξ j =x+η j (x 0 −x) is a point on the line segment betweenx and x 0 , and η j ∈ (0, 1) is an unknown parameter. By computing all second order partial derivatives of function h j , it can be derived that the Hessian matrix has the following form where I n is the identity matrix of order n and 0 m×n is the zero matrix of size m × n. The neighborhood O(x) can be defined for example as where R = min j r − s j 2 . For earth-based positioning applications, the condition x 0 ∈ O(x) required for Taylor's theorem to hold, can be considered satisfied because a simple detection technique is usually used to remove implausible position estimates that are too far from the Earth surface. Note that in the case of space-based positioning, a more elaborate analysis would be required. The nonlinear model (1) can be rewritten using the Taylor's theorem as where η = [η 1 , . . . , η ny ] T are unknown parameters and q : R 2nx+ny → R ny is a function whose j-th element is defined as Applying a simple algebraic manipulation, the model (16) can be written as As the state estimatex is the stationary point of the criterion V and thus meets the condition (7), the left hand side of the relation (18) is zero and the state estimation errorx = x 0 −x satisfies the relationx where

Linearized model
This subsection presents a common approach for deriving properties of the state estimate. Like the state estimatex the Jacobian matrix and the Hessian matrix are random variables. Since the state estimation errorx appears implicitly also on the right hand side of the relation (19), its exact distribution cannot be derived easily. To circumvent the issue, approximate properties of the state estimation error are usually derived using two following simplifying assumptions.
(i) The Jacobian matrix is evaluated at the computed realization of the state estimatex 1 and assumed to be constant when evaluating the mean value and the covariance matrix of the state estimation errorx. (ii) The function q that represents nonlinear behavior of the measurement model (1) is assumed to be the identically zero function.
Then the state estimation error has the Gaussian distribution with the following mean value and covariance matrix

Nonlinear model
A nonlinear analysis is performed to examine the validity of two above mentioned assumptions used in the previous subsection dealing with the linearized model. First, the sensitivity of the pseudoinverse H(x) † on x is examined, and then the worst case size of the quadratic term q is analyzed. To simplify the analysis, it is assumed that the covariance matrix Σ is diagonal where σ 2 is a variance of the measurement noise. Then the pseudoinverse of the Jacobian matrix is H( T be two state estimates whose position estimates satisfy the inequality where M is a given constant that is assumed to be significantly smaller than any distance between the receiver and the satellites.

Sensitivity of Jacobian matrix pseudoinverse
In the first part of the analysis, the sensitivity of the Jacobian matrix pseudoinverse to a change in the linearization point is evaluated. Although there are several ways of measuring the sensitivity of a pseudoinverse matrix, the norm H(x 2 ) † − H(x 1 ) † 2 is used because the change in the pseudoinverse matrix H(x) † induced by a change in the linearization pointx is of interest here. Another point of view is obtained when the expression (19) is used to derive an upper bound on the norm of the difference between two state estimation errors Thus the term H(x 2 ) † − H(x 1 ) † 2 determines how q(x 2 , x 0 , η 2 ) + f + v 2 is amplified when only the difference in the pseudoinverses is considered. Regardless of the other terms on the right hand side of (24), a general result on pseudoinverse sensitivity [13,14] leads to the following bound where E = H(x 2 )−H(x 1 ) ∈ R ny×nx is an additive perturbation matrix and µ = √ 2 is a constant that corresponds to the L 2 -norm. As the structure of the pseudoinverse is the same for both state estimatesx 1 andx 2 , upper bounds on H(x 1 ) † 2 2 and E 2 are derived. The squared norm of the pseudoinverse H(x 1 ) † satisfies where the third equality results from A T A and AA T having the same nonzero eigenvalues for any matrix A [15], σ max (A) denotes the maximum singular value of a matrix A, and λ max (A) and λ min (A) denote the maximum and minimum eigenvalues of a square matrix A, respectively. For certain satellite configurations, the Jacobian matrix H(x 1 ) can be almost a rank deficient matrix. In such a case, the minimum eigenvalue of H(x 1 ) T H(x 1 ) is arbitrarily close to zero and no meaningful upper bound on the squared norm of the pseudoinverse can be given without any additional assumption. Nevertheless, it is a common practice in GNSS applications to employ so called dilution of precision (DOP) metrics for measuring the quality of the satellite constellation and deciding whether or not the state estimate should be trusted. The geometric dilution of precision (GDOP) metric is defined as [16] GDOP = tr (H( and its relation to the eigenvalues of H(x 1 ) T H(x 1 ) −1 is given by where λ i (A) is the i-th eigenvalue of a square matrix A sorted in the descending order. If the maximum admissible geometric dilution of precision for a particular application denoted as GDOP max is chosen and considering that all eigenvalues of (H(x 1 ) T H(x 1 )) −1 are nonnegative, the following upper bound can be written Next, an upper bound on the norm of the perturbation matrix E is derived. When it is assumed that M is negligible compared to the distance between the receiver and a satellite, the j-th row of the perturbation matrix E can be approximated as follows The right eigenvector w 4 and zero eigenvalue λ 4 result from the function h j being linear in its fourth argument. The right eigenvector w 3 and zero eigenvalue λ 3 show that for the position estimation error r 0 −r lying entirely in the direction r − s j , the function h j behaves like it is linear. The first two eigenvectors w 1 and w 2 have the same nonzero eigenvalue and they are orthogonal to w 3 and w 4 . Thus, only the position estimation error that lies in the plane defined by these two eigenvectors makes the function h j nonlinear. Besides the eigenpair λ 4 , w 4 that is common to all functions h j , the other eigenpairs are different for each function q j and thus the value of the function q cannot be zero for nonzero state estimation errorx. Finally, as all eigenvalues are nonnegative, the Hessian matrix Q j is positive semidefinite and thus the range of the vector-valued function q is the positive orthant.
To continue in the analysis of the linearization error, the inequality for a quadratic norm with a symmetric matrix A is used [15]. It holds that where A ∈ R n×n is a square symmetric matrix and z ∈ R n is a vector. Using this inequality and substitutingr + η j (r 0 −r) − s j for a general point r in the expression for the eigenvalue λ 1 , the function q j is bounded as follows To eliminate the dependence on the unknown η j ∈ (0, 1), a lower bound on the denominator is derived. In general, it holds that and for η j ∈ (0, 1), the following strict inequality holds By combining (37) and (38), it is possible to derive the inequality which allows the right hand side of (36) to be further bounded as In the case of the GPS, the upper bound q j,max can be quantified approximately as follows. For near Earth's surface applications, a typical distance between a receiver and a satellite r 0 − s j is approximately 20.2 × 10 6 m. The dependence of the upper bound q j,max on the L 2 -norm of the position estimation error r 0 −r 2 is depicted in Fig. 1 and a detail for smaller position error estimates is given in Fig. 2. The figures show that the value of the function q j is lower than 2.5 m even if the L 2 -norm of the position estimation error reaches 10 km. It can be concluded that due to the great distance between the Earth's surface and satellites, the linear approximations of the functions h j are quite accurate even for a large neighborhood of the state estimate. It is also connected with a fast convergence of the GN method which is usually observed in practice.

Numerical example
The numerical example illustrates the bounds derived in the previous section. One particular state x 0 of the GPS receiver and two different satellite constellations are considered. The first constellation contains 10 visible satellites and is characterized by a good DOP with GDOP = 1.5269. The second constellation consists of 6 visible satellites only and has a worse DOP with GDOP = 5.9836 which is close to the chosen maximum allowable geometric dilution of precision GDOP max = 6. Since it is observed in practice that the norm of the position estimation error caused by the measurement noise rarely exceeds 100 m for low GDOP values, it is apriori assumed that M = 100 m. For both satellite constellations, the following results were obtained using 10000 Monte Carlo simulations where the pseudorange measurements are generated according to the model (1) with f = 0 and σ = 2.7. The histogram of the norm of difference between a chosen state estimate realization and other simulated realizations r l −r 1 2 for the good GDOP is depicted in Figure 3 and for the worse GDOP in Figure 4. It can be seen from these figures, that M is chosen such that the apriori assumption r l −r 1 2 ≤ M is satisfied. The computed bounds are summarized in Table 1. Note that the worse constellation has a smaller value of H max mainly because it contains fewer satellites compared to the good constellation. The realization of H(x l ) † − H(x 1 ) † 2 for the good constellation and the worse constellation are depicted in Figure 5 and Figure 6, respectively. By comparing the maximum value of H(x l ) † − H(x 1 ) † 2 with the upper bound H max in Table 1, it can be observed that the derived upper bound is not particularly tight and true values are much smaller. The realizations of the values of function q j are depicted in Figures 7 and 8. These histograms show that the nonlinear term is indeed always positive and like in the case of the bound on H(x 2 ) † − H(x 1 ) † 2 , this bound is not tight either.  Figure 3. Histogram of r l −r 1 2 for the good satellite constellation.

Conclusion
The paper dealt with a nonlinear analysis of the state estimate properties in GNSS positioning when a static estimator is used. Basic techniques such as Taylor's theorem and eigenvalue analysis were applied to analyze the state estimate and its error. Also the conditions under which the analysis holds were highlighted. The evaluation of the linearization error and sensitivity of the Jacobian matrix pseudoinverse showed that the linear approximation commonly used when deriving the properties of the state estimate is indeed reasonable. The analysis could be further extended with minor modifications to the case when a dynamical model and a nonlinear filter are used. Another potential extension consists in deriving a robust bound on the estimation error covariance matrix that can be used in a RAIM algorithm. Although the linear approximations have proven to be adequate by their successful use in many noncritical applications so far, the presented nonlinear analysis can serve as a more rigorous basis when GNSS positioning algorithms should be approved for use in safety critical applications.