The total optimal search criterion in solving the mixed integer linear model with GNSS carrier phase observations

Abstract

Existing algorithms for GPS ambiguity determination can be classified into three categories, i.e. ambiguity resolution in the measurement domain, the coordinate domain and the ambiguity domain. There are many techniques available for searching the ambiguity domain, such as FARA (Frei and Beutler in Manuscr Geod 15(4):325–356, 1990), LSAST (Hatch in Proceedings of KIS’90, Banff, Canada, pp 299–308, 1990), the modified Cholesky decomposition method (Euler and Landau in Proceedings of the sixth international geodetic symposium on satellite positioning, Columbus, Ohio, pp 650–659, 1992), LAMBDA (Teunissen in Invited lecture, section IV theory and methodology, IAG general meeting, Beijing, China, 1993), FASF (Chen and Lachapelle in J Inst Navig 42(2):371–390, 1995) and modified LLL Algorithm (Grafarend in GPS Solut 4(2):31–44, 2000; Lou and Grafarend in Zeitschrift für Vermessungswesen 3:203–210, 2003). The widely applied LAMBDA method is based on the Least Squares Ambiguity Search (LSAS) criterion and employs an effective decorrelation technique in addition. G. Xu (J Glob Position Syst 1(2):121–131, 2002) proposed also a new general criterion together with its equivalent objective function for ambiguity searching that can be carried out in the coordinate domain, the ambiguity domain or both. Xu’s objective function differs from the LSAS function, leading to different numerical results. The cause of this difference is identified in this contribution and corrected. After correction, the Xu’s approach and the one implied in LAMBDA are identical. We have developed a total optimal search criterion for the mixed integer linear model resolving integer ambiguities in both coordinate and ambiguity domain, and derived the orthogonal decomposition of the objective function and the related minimum expressions algebraically and geometrically. This criterion is verified with real GPS phase data. The theoretical and numerical results show that (1) the LSAS objective function can be derived from the total optimal search criterion with the constraint on the fixed integer ambiguity parameters, and (2) Xu’s derivation of the equivalent objective function was incorrect, leading to an incorrect search procedure. The effects of the total optimal criterion on GPS carrier phase data processing are discussed and its practical implementation is also proposed.

Appendix: Proof of orthogonal decomposition of objective function

In this appendix we will derive the orthogonal decomposition in linear model algebraically and geometrically. We begin with the linear Gauss–Markov model

$$ \begin{aligned}& {\mathbf{y}} = {\mathbf{A}}{\varvec{\upxi}} + {\mathbf{B}}{\varvec{\upeta}} + {\mathbf{e}} = {\mathbf{G}}{\varvec{\upgamma}} + {\mathbf{e}} \\ &{\text{where}}\;{\mathbf{G}} = [{\mathbf{A\;B}}],{\varvec{\upgamma}} = [{\varvec{\upxi^{\prime}\;{\varvec{\upeta^{\prime}}}}}]^{\prime},{\mathbf{y}} \in {\mathbb{R}}^{n},{\varvec{\upxi}} \in {\mathbb{R}}^{{m_{1}}},{\varvec{\upeta}} \in {\mathbb{R}}^{{m_{2}}} \\ & {\it{E}}\{{\mathbf{y}}\} = {\varvec{\rm{A}\upxi}} + {\varvec{\rm{B}\upeta}},{\it{D}}\{{\mathbf{y}}\} = {\varvec{\Upsigma}}_{{\mathbf{y}}}. \\ \end{aligned} $$

(30)

The best linear unbiased estimator (BLUE) is

$$ \begin{aligned}{\hat{\varvec{\upgamma}}}& =\left[{\begin{array}{l}{{\hat{\varvec{\upxi}}}} \\ {{\hat{\varvec{\upeta}}}} \\\end{array}} \right] =\left\{{[{\mathbf{A}}\;{\mathbf{B}}]^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} [{\mathbf{A}}\;{\mathbf{B}}]} \right\}^{- 1}[{\mathbf{A}}\;{\mathbf{B}}]^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} {\mathbf{y}} \\ {\hat{\varvec{\upxi}}} &\in{\mathbb{R}}^{{m_{1}}},\quad{\hat{\varvec{\upeta}}} \in{\mathbb{R}}^{{m_{2}}} \\ \end{aligned} $$

(31)

together with the residual vector

$$ {\hat{\mathbf{e}}} = {\mathbf{y}} - {\mathbf{A}}\hat{\varvec{\upxi}} - {\mathbf{B}}{\hat{\varvec{\upeta}}} = {\mathbf{y}} - {\mathbf{G}}{\hat{\varvec{\upgamma}}} $$

(32)

First we address the orthogonality of error, residual and parameter vectors and then study the orthogonality of the two parameter spaces when one set is conditioned on other set.

Applying the Pythagorean Theorem in a linear model, we have

$$ {\mathbf{y^{\prime}\varvec{\Upsigma}}}_{{\mathbf{y}}}^{- 1} {\mathbf{y}} = {\hat{\mathbf{y}}}^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} {\hat{\mathbf{y}}} + {\hat{\mathbf{e}}}^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} {\hat{\mathbf{e}}} $$

(33)

or in the format of the weighted squared norm

$$ ||{\mathbf{y}}||_{{{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1}}}^{2} =||{\hat{\mathbf{y}}}||_{{{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1}}}^{2} +||{\hat{\mathbf{e}}}||_{{{\varvec{\Upsigma}}_{{{\mathbf{y}}}}^{- 1}}}^{2} $$

(34)

i.e. the total sum of squares equals the sum of squares of adjusted observation vector and squares of the residual vector. Further, we can express this with the adjusted parameter vector

$$ \begin{aligned} ({\mathbf{y}} - {\mathbf{G}} {\varvec{\upgamma}})^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} ({\mathbf{y}} - {\mathbf{G}}{\varvec{\upgamma}})&= [{\mathbf{y}} - {\hat{\mathbf{G}}}{\varvec{\upgamma}} + {\mathbf{G}}({\hat{\varvec{\upgamma}}} - {\varvec{\upgamma}})]^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} [{\mathbf{y}} - {\hat{\mathbf{G}}}{\varvec{{\upgamma}}} + {\mathbf{G}}({\hat{\varvec{\upgamma}}} - {\varvec{\upgamma}})] \\ & = ({\mathbf{y}} - {\mathbf{G}}{\hat{\varvec{\upgamma}}})^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} ({\mathbf{y}} - {\mathbf{G}}{\hat{\varvec{\upgamma}}}) + ({\hat{\varvec{\upgamma}}} - {\varvec{\upgamma}})^{\prime}{\mathbf{G}}^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} {\mathbf{G}}({\hat{\varvec{\upgamma}}} - {\varvec{\upgamma}}) \\ &\quad + ({\mathbf{y}} - {\mathbf{G}}{\hat{\varvec{\upgamma}}})^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} {\mathbf{G}}({\mathbf{y}} - {\mathbf{G}}{\varvec{\upgamma}}) + ({\hat{\varvec{\upgamma}}} - {\varvec{\upgamma}})^{\prime}{\mathbf{G}}^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} ({\mathbf{y}} - {\mathbf{G}}{\hat{\varvec{\upgamma}}}) \\ & = ({\mathbf{y}} - {\mathbf{G}}{\hat{\varvec{\upgamma}}})^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} ({\mathbf{y}} - {\mathbf{G}}{\hat{\varvec{\upgamma}}}) + ({\hat{\varvec{\upgamma}}} - {\varvec{\upgamma}})^{\prime}{\mathbf{G}}^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} {\mathbf{G}}({\hat{\varvec{\upgamma}}} - {\varvec{\upgamma}}) \\ & ={\hat{\mathbf{e}}}^{\prime}{{\varvec{\Upsigma}}}_{{\mathbf{y}}}^{- 1} {\hat{\mathbf{e}}} + ({\hat{\varvec{\upgamma}}} - {\varvec{\upgamma}})^{\prime}{\mathbf{G}}^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} {\mathbf{G}}({\hat{\varvec{\upgamma}}} - {\varvec{\upgamma}}) \\ & = {\hat{\mathbf{e}}}^{\prime}\varvec{\Upsigma}_{{\mathbf{y}}}^{- 1} {\hat{\mathbf{e}}} + {\varvec{\Updelta}} \\ \end{aligned} $$

(35)

since ${\mathbf{G}}^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} ({\mathbf{y}} - {\mathbf{G}}{\hat{\varvec{\upgamma}}}) = ({\mathbf{y}} - {\mathbf{G}}{\hat{\varvec{\upgamma}}})^{\prime}{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1} {\mathbf{G}} = 0.$ This can be also represented in the weighted squared norm

$$ ||{\mathbf{y}} - {\mathbf{G}}{\varvec{\upgamma}}||_{{{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1}}}^{2} =||{\mathbf{y}} - {\mathbf{G}}{\hat{\varvec{\upgamma}}}||_{{{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1}}}^{2} +||{\hat{\varvec{\upgamma}}} -{\varvec{\upgamma}}||_{{{\varvec{\Upsigma}}_{{{\hat{\varvec{\upgamma}}}}}^{- 1}}}^{2} $$

(36)

This is the first step of orthogonal decomposition of the objective function. The geometric illustration of both orthogonal decompositions can be found in Fig. 1a.

Fig. 1

A geometry of the orthogonal decomposition of $||{\mathbf{y}}||_{{{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1}}}^{2} \,{\text{and}}\,||{\mathbf{y}} - {\mathbf{G}}{\varvec{\upgamma}}||_{{{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1}}}^{2}$

Let us consider the second term of the sum of squares of the error vector in Eq. 36 and assume that the second parameter set is given as ${\varvec{\upeta}} = {\bar{\varvec{\upeta }}},$ which has to be conditioned by the first parameter set, then we have

$$ \begin{aligned}||{\hat{\varvec{\upgamma}}} -{\varvec{\upgamma}}||_{{{\varvec{\Upsigma}}_{{{\hat{\varvec{\upgamma}}}}}^{- 1}}}^{2} &= ({\hat{\varvec{\upgamma}}} -{\varvec{\upgamma}})^{\prime}{\mathbf{G^{\prime}\varvec{\Upsigma}}}_{{\mathbf{y}}}^{- 1} {\mathbf{G}}({\hat{\varvec{\upgamma}}} -{\varvec{\upgamma}}) \\ & = [({\hat{\varvec{\upxi}}} -{\varvec{\upxi}})^{\prime}({\hat{\varvec{\upeta}}} -{\bar{\varvec{\upeta}}})^{\prime}]\left[{\begin{array}{ll}{{\mathbf{A^{\prime}\varvec{\Upsigma}}}_{{\mathbf{y}}}^{- 1}{\mathbf{A}}} & {{\mathbf{A^{\prime}\varvec{\Upsigma}}}_{{\mathbf{y}}}^{-1} {\mathbf{B}}} \\ {{\mathbf{B^{\prime}\varvec{\Upsigma}}}_{{\mathbf{y}}}^{- 1} {\mathbf{A}}} &{{\mathbf{B^{\prime}\varvec{\Upsigma}}}_{{\mathbf{y}}}^{- 1}{\mathbf{B}}} \\ \end{array}} \right]\left[{\begin{array}{ll} {{\hat{\varvec{\upxi}}} -{\varvec{\upxi}}} \\ {{\hat{\varvec{\upeta}}} -{\bar{\varvec{\upeta}}}} \\ \end{array}} \right] \\ & =[({\hat{\varvec{\upxi}}} -{\varvec{\upxi}})^{\prime}({\hat{\varvec{\upeta}}} -{\bar{\varvec{\upeta}}})^{\prime}]\left[{\begin{array}{ll}{{\mathbf{N}}_{11}} & {{\mathbf{N}}_{12}} \\ {{\mathbf{N}}_{21}} & {{\mathbf{N}}_{22}} \\ \end{array}} \right]\left[{\begin{array}{ll} {{\hat{\varvec{\upxi}}} -{\varvec{\upxi}}} \\ {{\hat{\varvec{\upeta}}} -{\bar{\varvec{\upeta}}}} \\ \end{array}} \right] \\ & =({\hat{\varvec{\upxi}}} -{\varvec{\upxi}})^{\prime}{\mathbf{N}}_{11} ({\hat{\varvec{\upxi}}} - {\varvec{\upxi}}) + ({\hat{\varvec{\upeta}}} -{\bar{\varvec{\upeta}}})^{\prime}{\mathbf{N}}_{21}({\hat{\varvec{\upxi}}} - {\varvec{\upxi}}) \\ &\quad+({\hat{\varvec{\upxi}}} -{\varvec{\upxi}})^{\prime}{\mathbf{N}}_{12} ({\hat{\varvec{\upeta}}} - {\bar{\varvec{\upeta}}}) + ({\hat{\varvec{\upeta}}} -{\bar{\varvec{\upeta}}})^{\prime}{\mathbf{N}}_{22}({\hat{\varvec{\upeta}}} - {\bar{\varvec{\upeta}}}) \\ & =({\hat{\varvec{\upxi}}} -{\varvec{\upxi}})^{\prime}{\mathbf{N}}_{11} ({\hat{\varvec{\upxi}}} - {\varvec{\upxi}}) + 2({\hat{\varvec{\upxi}}} -{\varvec{\upxi}})^{\prime}{\mathbf{N}}_{12} ({\hat{\varvec{\upeta}}} - {\bar{\varvec{\upeta}}}) + ({\hat{\varvec{\upeta}}} -{\bar{\varvec{\upeta}}})^{\prime}{\mathbf{N}}_{22}({\hat{\varvec{\upeta}}} - {\bar{\varvec{\upeta}}}) \\\end{aligned} $$

(37)

With the relationship 27 we have $({\hat{\varvec{\upxi}}} - {\bar{\varvec{\upxi}}}_{{|{\varvec{\upeta}} = {\bar{\varvec{\upeta}}}}}) = -{\mathbf{N}}_{11}^{- 1} {\mathbf{N}}_{12} ({\hat{\varvec{\upeta}}} - {\bar{\varvec{\upeta}}})$ and note that the variance-covariance of the estimators with the Schur complement are

$$ \begin{aligned}{\varvec{\Upsigma}}_{{\hat{\varvec{\upeta}}}} &=({\mathbf{N}}_{22} - {\mathbf{N}}_{21} {\mathbf{N}}_{11}^{- 1}{\mathbf{N}}_{12})\\{\varvec{\Upsigma}}_{{\bar{\varvec{\upxi}}}_{{{\varvec{|\upeta}}} ={{{{{\bar{\varvec{\upeta}}}}}}}}} &= {\varvec{\Upsigma}}_{{\hat{\varvec{\upxi}}}} - {\mathbf{N}}_{11}^{- 1} {\mathbf{N}}_{12}\Upsigma_{{\hat{\varvec{\upeta}}}} {\mathbf{N}}_{21} {\mathbf{N}}_{11}^{-1} = {\mathbf{N}}_{11}^{- 1}\end{aligned} $$

(38)

Then

$$\begin{aligned}||{\hat{\varvec{\upgamma}}} -{\varvec{\upgamma}}||_{{{\varvec{\Upsigma}}_{{\hat{\varvec{\upgamma}}}}^{- 1}}}^{2} & = ({\hat{\varvec{\upxi}}} -{\varvec{\upxi}})^{\prime}{\mathbf{N}}_{11} ({\hat{\varvec{\upxi}}} - {\varvec{\upxi}}) - 2({\hat{\varvec{\upxi}}} -{\varvec{\upxi}})^{\prime}{\mathbf{N}}_{11} ({\hat{\varvec{\upxi}}} - {\varvec{\bar{\upxi}}}_{{|{\varvec{\upeta}} = {{\bar{\varvec{\upeta}}}}}}) +({\hat{\varvec{\upeta}}} - {\bar{\varvec{\upeta}}})^{\prime}{\mathbf{N}}_{22} ({\hat{\varvec{\upeta}}} -{\bar{\varvec{\upeta}}}) \\ & = ({\hat{\varvec{\upxi}}} -{\varvec{\upxi}})^{\prime}{\mathbf{N}}_{11} ({\hat{\varvec{\upxi}}} - {\varvec{\upxi}}) - ({\hat{\varvec{\upxi}}} -{\varvec{\upxi}})^{\prime}{\mathbf{N}}_{11} ({\hat{\varvec{\upxi}}} - {\bar{\varvec{\upxi}}}_{{|{\varvec{\upeta}} = {{\bar{\varvec{\upeta}}}}}}) -({\hat{\varvec{\upxi}}} -{\varvec{\upxi}})^{\prime}{\mathbf{N}}_{11} ({\hat{\varvec{\upxi}}} - {\bar{\varvec{\upxi}}}_{{|{\varvec{\upeta}} = {{\bar{\varvec{\upeta}}}}}}) +({\hat{\varvec{\upeta}}} - {\bar{\varvec{\upeta}}})^{\prime}{\mathbf{N}}_{22} ({\hat{\varvec{\upeta}}} -{\bar{\varvec{\upeta}}}) \\ & = ({\hat{\varvec{\upxi}}} -{\varvec{\upxi}})^{\prime}{\mathbf{N}}_{11} [({\hat{\varvec{\upxi}}} - {\varvec{\upxi}}) - ({\hat{\varvec{\upxi}}} -{\bar{\varvec{\upxi}}}_{{|{\varvec{\upeta}} = {{\bar{\varvec{\upeta}}}}}})] -({\hat{\varvec{\upxi}}} -{\varvec{\upxi}})^{\prime}{\mathbf{N}}_{11} ({\hat{\varvec{\upxi}}} - {\bar{\varvec{\upxi}}}_{{|{\varvec{\upeta}} = {{\bar{\varvec{\upeta}}}}}}) +({\hat{\varvec{\upeta}}} - {\bar{\varvec{\upeta}}})^{\prime}{\mathbf{N}}_{22} ({\hat{\varvec{\upeta}}} -{\bar{\varvec{\upeta}}}) \\ & = ({\bar{\varvec{\upxi}}}_{{|{\varvec{\upeta}} = {{\bar{\varvec{\upeta}}}}}} -{\varvec{\upxi}})^{\prime}{\mathbf{N}}_{11} ({\bar{\varvec{\upxi}}}_{{|{\varvec{\upeta}} = {{\bar{\varvec{\upeta}}}}}} - {\varvec{\upxi}}) +({\hat{\varvec{\upeta}}} - {\bar{\varvec{\upeta}}})^{\prime}({\mathbf{N}}_{22} - {\mathbf{N}}_{21}{\mathbf{N}}_{11}^{- 1} {\mathbf{N}}_{12})({\hat{\varvec{\upeta}}} - {\bar{\varvec{\upeta}}}) \\ & = ({\bar{\varvec{\upxi}}}_{{|{\varvec{\upeta}} = {{\bar{\varvec{\upeta}}}}}} -{\varvec{\upxi}})^{\prime}{\varvec{\Upsigma}}_{{{\bar{\varvec{\upxi}}}_{{|{\varvec{\upeta}} = {{\bar{\varvec{\upeta}}}}}}}}^{- 1}({\bar{\varvec{\upxi}}}_{{|{\varvec{\upeta}} = {{\bar{\varvec{\upeta}}}}}} -{\varvec{\upxi}}) + ({\hat{\varvec{\upeta}}} - {\bar{\varvec{\upeta}}})^{\prime}{\varvec{\Upsigma}}_{{{\hat{\varvec{\upeta}}}}}^{-1} ({\hat{\varvec{\upeta}}} - {\bar{\varvec{\upeta}}}) \\\end{aligned}$$

(39)

In terms of the weighted squared norm, it can be expressed as

$$ \begin{gathered} ||{\mathbf{y}} - {\mathbf{G}}\varvec{\upgamma }||_{{{\varvec{\Upsigma}}_{{\mathbf{y}}}^{ - 1} }}^{2} = ||{\mathbf{y}} - {\mathbf{G}}{\hat{\varvec{\upgamma }}}||_{{{\varvec{\Upsigma}}_{{\mathbf{y}}}^{ - 1} }}^{2} + ||\hat{{\varvec{\upgamma}} } - \varvec{\upgamma }||_{{{\varvec{\Upsigma}}_{{\mathbf{y}}}^{ - 1} }}^{2} \hfill \\ {\text{i}} . {\text{e}} .\hfill \\ ||{\mathbf{y}} - {\mathbf{A}}\varvec{\upxi } - {\mathbf{B}}\varvec{\upeta}||_{{{\varvec{\Upsigma}}_{{\mathbf{y}}}^{ - 1} }}^{2} = ||{\hat{\mathbf{e}}}||_{{{\varvec{\Upsigma}}_{{\mathbf{y}}}^{ - 1} }}^{2} + ||\bar{\varvec{\upxi }}_{{{|\varvec{{\upeta}= {\bar{\varvec{\upeta}}}}}}} - \user2{\upxi }||_{{{\varvec{\Upsigma}}_{\bar{{\varvec{\upxi}}}|_{\varvec{\upeta} = {\bar{\varvec{\upeta}}}}}^{ - 1} }}^{2} + ||{\hat{\varvec{\upeta }}} - {\bar{\varvec{\upeta} }}||_{{{\varvec{\Upsigma}}_{\hat{\varvec{\upeta}}}^{ - 1} }}^{2} \hfill \\ \end{gathered} $$

(40)

Equation (40) explains the so-called orthogonal decomposition (Eq. 6) generally used in mixed integer least squares. Fig. 2a illustrates the geometry of the above.

Fig. 2

A geometry of the orthogonal decomposition of $||{\mathbf{y}} - {\mathbf{A}}{\varvec{\upxi}} - {\mathbf{B}}{\varvec{\upeta}}||_{{{\varvec{\Upsigma}}_{{\mathbf{y}}}^{- 1}}}^{2}$when the coordinate parameter vector is conditioned on the given ambiguity vector