Abstract
Mixed integer–real least squares (MIRLS) estimation still has two open scientific problems, i.e., the validation of results and computational efficiency for a large number of satellites. This paper presents and discusses a non-conventional approach to MIRLS estimation, which belongs to the ambiguity function method (AFM) class. Because the solution is searched for in the constant three-dimensional coordinate domain instead of the n-dimensional ambiguity domain, the computational efficiency does not depend as much on the number of satellites as it does in conventional MIRLS estimation. Simple numerical pretests have shown that the reliability and precision of results from the presented approach and the conventional MIRLS estimation are exactly the same. Hence, the presented approach, contrary to AFM, may be treated as MIRLS estimation. Furthermore, the presented approach is a few hundred times faster than AFM and may be considered in (near) real-time GNSS positioning. In light of the above, the new field of research on MIRLS estimation may be opened.
Similar content being viewed by others
Explore related subjects
Discover the latest articles, news and stories from top researchers in related subjects.Avoid common mistakes on your manuscript.
Introduction
Precise GNSS positioning requires resolving the so-called mixed integer–real problem to determine carrier phase integer ambiguities. Teunissen developed and ordered the entire family of ambiguity estimators and not only the integer ones. Three classes may be distinguished in this family. Integer (I) estimators constitute the first and most important class (Teunissen 1999a). Well-known examples are integer least squares (ILS), integer bootstrapping (IB), and integer rounding (IR). These estimators recognize to some degree the correlation between estimated ambiguities. Integer aperture (IA) estimators belong to the second class (Teunissen 2003). These estimators unify I-estimation with validation and can adopt both integer and real values. Estimators in this class differ in applied methods of I-estimation and validation. Integer equivariant (IE) estimators constitute the third and most general class (Teunissen 2002). These estimators are always real-valued. Partial (P) estimation is a significant expansion of the above ambiguity estimation theory. This approach was first introduced in Teunissen et al. (1999). In P-estimation, only a subset of ambiguities can be an integer, while the other ambiguities must be real-valued. Examples of P-estimators are the PI-estimator (Verhagen et al. 2011; Brack 2016) and PIA-estimator (Brack 2015). A common feature of the above family of estimation is that searching for the final solution is performed in the ambiguity domain.
The second category of algorithms includes the very first ambiguity estimation technique developed, namely the ambiguity function method (AFM). This method was first introduced by Counselman and Gourevitch (1981). Remondi (1984, 1990) used this method extensively for GPS static positioning and also for pseudokinematic positioning. These methods utilize certain properties of the chosen trigonometric functions, which have known values for the integer arguments, without determining ambiguities. The solution search proceeds in the coordinate domain.
The last category includes the simplest techniques that use C/A or P-code pseudoranges directly to estimate the ambiguities of corresponding carrier phase measurements. The precision of the pseudoranges is not sufficient to determine the integer ambiguities, and linear inter-frequency combinations are usually used to estimate the ambiguities. Computations proceed in the measurement domain.
In the Gaussian case, the ILS estimator is the optimal solution (Teunissen 1999b). The optimality criterion used is that of maximizing the probability of correct integer ambiguity estimation, the so-called success rate (SR). However, the search for the ILS estimator is time-consuming and hence ineffective in (near) real-time applications. This results from the occurrence of strong correlations between estimated ambiguities. To resolve the problem, numerous methods were developed, increasing the computational efficiency of the ILS estimation (Frei and Beutler 1990; Teunissen 1995; Chen and Lachapelle 1995; Kim and Langley 2000a, b; Xu 2001, 2012; Chang et al. 2005; Zhou 2011; Jazaeri et al. 2014). The computational efficiency of the ILS estimation methods existing now is satisfactory for (near) real-time positioning. Hence, these methods, in particular, the least squares ambiguity decorrelation adjustment (LAMBDA), are today widely used for the mixed integer–real least squares (MIRLS) estimation. However, this applies to only up to a dozen or so ambiguities. In the case of a larger number of ambiguities, the computation time may be too long for (near) real-time applications (Jazaeri et al. 2014). It is still a challenge for GNSS investigators. The validation of estimated ambiguities is a second open problem (Verhagen and Teunissen 2013).
We discuss and pretest a completely different approach which may open the new field of research on the MIRLS estimation. Contrary to the conventional MIRLS estimation, the search for a solution occurs in a three-dimensional coordinate domain and is much faster than in AFM. Because the solution is searched for in the constant domain instead of the n-dimensional domain, the computational efficiency does not depend as much on the number of satellites as it does in the conventional MIRLS estimation. The presented approach is, in fact, a variant of the modified ambiguity function approach (Cellmer et al. 2010; Cellmer 2012). Hence, it has been conventionally named the modified ambiguity function approach-integer least squares (MAFA-ILS).
Mixed integer–real least squares estimation
The equation of double difference (DD) carrier phase measurements may be represented as:
where \( {\varvec{\uprho}} \) is the vector of DD geometrical distances, \( \lambda \) is the wavelength of the signal, a is the vector of DD carrier phase ambiguities, and e is the vector of DD carrier phase measurement errors. It is assumed that biases such as satellite ephemeris errors, tropospheric and ionospheric delays and ranging errors caused by multipath are absent. Linearization of the observation equations, with respect to the unknown parameters, gives the well-known mixed integer–real linear model:
where y is the vector of ‘observed minus computed’ DD carrier phase measurements, A and B are the design matrices, b is the vector that contains the increments of the unknown baseline components, Zn is the space of integers, R3 is the space of reals, \( \sigma_{0}^{2} \) is the variance factor, and Q y and C y are the cofactor and covariance matrix of DD carrier phase measurements, respectively. In order to solve for this system of equations, the constrained LS principle is applied:
where \( {\hat{\mathbf{e}}} \) and \( {\hat{\mathbf{a}}} \) are the LS estimate of e and a, respectively, and \( {\mathbf{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{a} }} \) is the ILS estimate of a. The final constrained baseline solution has the following form:
where \( {\hat{\mathbf{b}}} \) is the LS estimate of b. The ILS estimation of ambiguities is the basic task. Unfortunately, the real-valued LS estimator of ambiguities is usually highly correlated, and its confidence region is usually extremely elongated. Hence, finding the ILS solution may be extremely time-consuming. As it has already been mentioned, the ambiguity search space is most often reduced to increase the computational efficiency. The integer ‘decorrelation’ plays a significant role. The newly transformed ambiguities are almost independent, and the search for the ILS solution is then relatively quick.
Modified Ambiguity Function Approach
Typically, the noise of carrier phase measurement is about 0.01 cycle. Thus, the errors of DD carrier phase measurements should be much less than half a cycle. Hence, on the basis of (1), and taking into account the integer nature of ambiguities, the vector of ambiguities shall then have the form: \( {\mathbf{a}} = round({\varvec{\Phi}} - {\varvec{\uprho}}/\lambda ) \), where the round is a function of rounding to the nearest integer value. For m epochs, the vector of ambiguities shall obviously have the form of \( \left[ {\begin{array}{*{20}c} \ldots & {{\mathbf{a}}_{j}^{\text{T}} } & \ldots \\ \end{array} } \right]^{\text{T}} \), \( j = 1, \ldots ,m \). After substituting the above vector of ambiguities to (1), one obtains the following equation:
For the purpose of linearization, a differentiable function in the place of the term on the right side of (5) is proposed (Cellmer 2012). A differentiable single component has the following form:
where \( r_{i} = \varPhi_{i} - \rho_{i} /\lambda \), \( \rho_{i} = \rho \left( {\mathbf{x}} \right) \), \( {\mathbf{x}} = \left[ {\begin{array}{*{20}c} x & y & z \\ \end{array} } \right]^{T} \) is the vector of the receiver coordinates and \( i = 1, \ldots ,n \). The derivative of (6) is: \( \partial e_{i} /\partial {\mathbf{x}} \) = (\( \partial e_{i} /\partial r_{i} \))·(\( \partial r_{i} /\partial {\mathbf{x}} \)) = 1· \( ( - \partial \rho_{i} /\partial {\mathbf{x}})/\lambda \). Hence, after a Taylor series expansion of (6):
where \( x^{0} ,y^{0} ,z^{0} \) are the approximate coordinates of the receiver and \( dx,dy,dz \) are the increments of the unknown baseline components, one obtains a linear mathematical model:
where \( {\varvec{\updelta}} = \left( {{\varvec{\Phi}} - {\varvec{\uprho}}^{0} /\lambda } \right) - round\left( {{\varvec{\Phi}} - {\varvec{\uprho}}^{0} /\lambda } \right) \) and \( {\varvec{\uprho}}^{0} = \rho ({\mathbf{x}}^{0} ) \) is the vector of computed DD geometrical distances. In order to solve for this equation, the LS principle is applied:
where \( {\mathbf{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{e} }} \) is the LS estimate of e, but conditioned on \( {\mathbf{a}} \in Z^{n} \). The constrained baseline solution is the following estimator:
The ambiguity parameters are not present in model (8). Nevertheless, the above formulation yields results that fulfill the condition of integer ambiguities.
Problem of the approximate coordinates
After linearization, the vector of ambiguities has the form:
The following condition must be met for (11) to be correct:
where \( {\mathbf{e}}_{{\rho^{0} }} \) is the vector of the computed DD geometrical distance errors. A so-called good Voronoi cell is a graphical interpretation of Eq. (11) in the coordinate domain (Cellmer 2012). The entire coordinate domain is filled—without gaps or overlaps—with Voronoi cells, V, i.e.,
and all coordinates, x, from the given cell pull in to the center of this cell, x(v). The coordinates of the center of each Voronoi cell realize some integer ambiguities, v, and the coordinates of the center of good Voronoi cell realize true integer ambiguities, a. In the GNSS literature, the pull-in regions from the Teunissen theory are sometimes referred to also as Voronoi cells, where the distance is measured in the metric of \( {\mathbf{Q}}_{{\hat{a}}} \) (Xu 2006). However, the Voronoi cells in MAFA are different. Figure 1 presents an example of a few Voronoi cells and some cloud of float positions in 2D space. The good Voronoi cell is located in the center. These cells were empirically generated based on (11) from a dense set of points, for some model (8). The shape of these cells depends on the geometrical configuration of satellites and their size on the wavelength of the phase signal.
Condition (12) is a necessary and sufficient condition to obtain a correct solution in MAFA. The determination of coordinates that satisfy this condition, i.e., which are located in the good Voronoi cell as shown in Fig. 1, is the basic problem. Unfortunately, approximate coordinates obtained, for example, from the float solution frequently are not located in the good Voronoi cell. Therefore, a cascade adjustment was suggested in MAFA, using a few linear combinations of phase signals. The performed tests have shown that such an approach significantly increases the probability of determining the coordinates located in the good Voronoi cell and hence the probability of good solution. Some improvement to the above approach was also suggested in Cellmer (2013). Unfortunately, it is still not possible to always determine the coordinates which are accurate enough, i.e., which are located in the good Voronoi cell. There is always a small risk that these coordinates may be located in the bad Voronoi cell. However, another solution is possible here. Namely, it is always possible to find such approximate coordinates, which shall give the ILS solution!
Modified ambiguity function approach-integer least squares
The optimization problem in MAFA has the form of (8) and (9). However, despite the application of the LS principle and enforcing an integer ambiguity nature, the MAFA method is not the same as the ILS ambiguity estimation. The MAFA method is the LS estimation of the baseline components, using certain, most often true, integer ambiguities. Instead, a variant of MAFA has been suggested in this section, which ‘estimates’ ambiguities and specifically implements the ILS ambiguities.
Keep in mind that in the MIRLS estimation one searches for such a vector of integer ambiguities which is the solution of optimization (2) and (3). Determination of coordinates, which are located in a good Voronoi cell, is the basic problem of MAFA. As it has already been mentioned, determination of such coordinates is not always possible. However, finding such coordinates, which shall be the global solution of optimization (8) and (9), is always possible! Such solution is named MAFA-ILS. Figure 2 shows the graphic interpretation of the MAFA-ILS minimization form (9), against a background of the orthogonal decomposition of the constrained LS minimization forms (3).
The first step in MAFA-ILS is intended only to provide initial approximate coordinates, to determine the search region center. However, step 2 is of basic importance. This step is intended for finding such approximate coordinates, which give the global solution of optimization (8) and (9). It should be emphasized that such coordinates are not unique. These could be any coordinates located in the ILS Voronoi cell, where the coordinates at the center of the ILS Voronoi cell realize the ILS ambiguities. It is important to note that the ILS Voronoi cell shall not always be a good Voronoi cell. This obviously results from the stochastic nature of model (8) and depends on its strength. It is a general problem that for weak models the correct solution and the ILS solution are sometimes in different places.
Trivial search procedure
A trivial search procedure may consist of generating an appropriate cloud of points and select such point, which is located in the ILS Voronoi cell. Such a cloud must obviously contain at least one point located in the ILS Voronoi cell. The most trivial solution may consist in generating an appropriate cube of grid points around the initial position, e.g., float position, \( {\hat{\mathbf{x}}} \). Such a set of points may be easily generated by means of the following recurrent formula:
where
and \( {\mathbf{G}}_{1} (1 \times g) = [\begin{array}{*{20}c} \ldots & { - 2} & { - 1} & 0 & 1 & 2 & \ldots \\ \end{array} ] \) is the matrix of grid points around the initial position, g is the number of grid layers, d is the distance between neighboring points, \( {\mathbf{X}}_{f} (3 \times g^{3} ) = {\hat{\mathbf{x}}}(3 \times 1) \cdot {\mathbf{1}}\left( {1 \times g^{3} } \right) \) is the matrix of initial coordinates (the same in each column), \( i = 1, \ldots ,c \), \( c = g^{3} \) means the number of all points, and ‘\( \otimes \)’ is the Kronecker product symbol. Now it is necessary to find the point located in the ILS Voronoi cell. A most trivial way may consist in calculating the square form (9) for each candidate with respect to model (8). The candidate, which will give the smallest value of this form, will be the point located in the ILS Voronoi cell. For this point, the ambiguity values (11) shall be the same for all m epochs and equal to the ILS ambiguity estimator (3). This hypothesis will be studied in the numerical section.
Density and size of the appropriate search region
The crucial question is: What could be the minimum size of the search region and the minimum density of grid candidates so that the set would contain at least one point located in the ILS Voronoi cell? It may be noticed that to satisfy the above condition at least one candidate must fall into each Voronoi cell. So, the density of grid candidates depends on the size and shape of Voronoi cells. Figure 3 shows an example of a Voronoi cell. This cell was generated empirically based on (11) from a dense set of points, for some model (8). Voronoi cells are symmetric and convex; therefore, an upper bound for the minimum density, dmax, can be calculated empirically based on given model (8) as follows:
where \( max\left( {\Delta x} \right) \), \( max\left( {\Delta y} \right) \), \( max\left( {\Delta z} \right) \) are maximum differences of the coordinates of grid points, making an empirical Voronoi cell, in directions passing through the center of the cell along the axes of the coordinate system as shown in Fig. 3. The upper bound (17) can be helpful in determining the minimum density of appropriate search cube.
Instead, the minimum size of an appropriate search region, e.g., the minimum length of the cube side, s, depends on the distance between the ILS Voronoi cell and the center of the search region. This problem may also be resolved empirically, for example, based on the Monte Carlo (MC) simulations. Specifically, it is possible to generate a cloud of float positions, which concern a unique ILS solution. Based on such a cloud, it is possible to determine the coordinate differences between this solution and the most distant float position. On the basis of such differences, it is now possible to determine, for example, the minimum length of the search cube side. It is easily noticeable that this length should not be smaller than a double value of the largest difference of coordinates between the unique ILS solution and the most distant float position, i.e.,
where \( {\hat{\mathbf{X}}}\left( {3 \times N} \right) \) is the matrix of float positions, which concern the unique ILS solution, and \( {\mathbf{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{X} }}\left( {3 \times N} \right) \) is the matrix of the unique ILS positions (the same in each column), where N is the number of simulations. Then, the point located in the ILS Voronoi cell shall be in the search cube with probability \( P = 1 - 1/N \). Figure 4 shows some unique ILS Voronoi cell, its corresponding cloud of float positions for N = 5000, and the appropriate search cube. In this case, the maximum coordinate difference was \( max\left( {\Delta Y} \right) = 0.68{\text{ m}} \); hence, \( s_{\min} = 1.36{\text{ m}} \).
It should be noted that the above solution may be applied to formulate a certain universal solution. Namely, for many configurations of satellites the values of (18) can be calculated, and then, the relationships between those values and values characterizing the satellite configurations, e.g., the PDOP values, may be determined. In this way, a universal formula for calculation of the minimum size of appropriate search cube for any satellite configuration can be empirically determined.
Necessary condition for the correct solution
A cloud of float positions for a unique ILS Voronoi cell, Fig. 4, creates a certain region to be called the region of the ILS float positions, \( S_{{{\mathbf{x}}({\mathbf{v}})}} :R^{3} \). In practice, this region can be approximated, for example, by means of a convex hull, which is the smallest surface comprising all points of the given cloud as shown in Fig. 5.
The coordinates of this region center are also the coordinates of the center of the ILS Voronoi cell and realize integer ambiguities. Hence, in the context of the Teunissen theory, it is possible to say that the region of the ILS float positions represents the pull-in region of the ILS ambiguities, but in the coordinate domain. Moreover, the Teunissen theory states that if some float ambiguities are located in a pull-in region of the ILS ambiguities, at which center there are correct integer ambiguities, then the ILS ambiguity estimator will be correct. This is a definition of necessary and sufficient condition to obtain a correct solution in the MIRLS estimation. Also in MAFA-ILS, we have a certain analogy. That is, if the float position is located outside the region of the ILS float positions, corresponding to correct ambiguities, then the MAFA-ILS ambiguity estimator will not be correct, i.e.,
In the context of the Teunissen theory, this may be interpreted as a definition of a necessary condition for a correct solution in MAFA-ILS. For example, for the model in Fig. 4, it is possible to state that if a true error of float position would be larger than 0.68 m, then MAFA-ILS would give an incorrect solution. Instead, the problem of sufficient condition for a correct solution in MAFA-ILS has no analogy to the MIRLS estimation. Namely, if the float position is located in the region of ILS float positions, corresponding to correct ambiguities, then despite that, the MAFA-ILS ambiguity estimator may not be correct, i.e.,
This results from the fact that individual regions of the ILS float positions have gaps although the space filled with these regions has no gaps, and these regions overlap, i.e.,
Thus, the sufficient condition for the correct solution in MAFA-ILS cannot be defined.
Implicit integer ambiguities
Each candidate pulls in to the center of its Voronoi cell, and the coordinates of the center of each Voronoi cell realize some integer ambiguities, in accordance with (11). So, one may say that in MAFA-ILS, as in the case of MIRLS estimation, integer ambiguities are also considered, but in this case in an implicit form. Moreover, as in MIRLS estimation, the number of all possible integer ambiguities in the appropriate search region is reduced here. However, this reduction is different. In the MIRLS estimation, an appropriate set of all possible integer ambiguities is reduced most often by means of integer ‘decorrelation.’ Instead, in MAFA-ILS this set is reduced by filtering using function (11) instead of by decreasing. More precisely, this function takes into account only such vectors of integer ambiguities for which all observations fit model (8) well. ‘Well’ means that the misclosures are relatively small, i.e., a priori \( {\varvec{\updelta}} \in \left\langle { - 0.5,0.5} \right\rangle \) cycle and a posteriori:
where \( {\mathbf{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{e} }} \) is the residual vector of optimization (8) and (9). Such an implicit set of integer ambiguities, expressed in the three-dimensional coordinate domain, is significantly reduced. This is why, the number of points in MAFA-ILS, even with a trivial search procedure, can be relatively small.
MAFA-ILS in the AFM context
Broadly and without loss of generality, in the case of AFM one searches for station coordinates mostly that maximize/minimize the cosine/sine of the residuals. If one further assumes that all observations are equally weighted, then by means of AFM one can look for an approximate ILS solution for one epoch. For a single baseline, this task may be presented in the following form (Leick et al. 2015, p. 317):
A relatively dense grid cube, e.g., \( d = 0.01{\text{ m}} \), is a set of candidates. Moreover, this value of density is a possible error in the determination of the minimum of function (23).
Let it be noted now that, with respect to the symbols of model (8), the optimization condition (23) may be—with a good approximation—written as:
because \( \left| {2\delta_{k} } \right| \approx \left| {\sin \left\{ {\pi \left( {\varPhi_{k} - \rho_{k}^{0} /\lambda } \right)} \right\}} \right| \). It is seen that in the case of the approximate ILS solution, AFM minimizes squares of a priori misclosures of model (8). Whereas if one assumes that all observations are equally weighted, MAFA-ILS minimizes squares of a posteriori misclosures of model (8), i.e.,
Generally, this is a theoretical relationship between MAFA-ILS and an approximate ILS solution by means of AFM.
Also, certain practical differences result from the above. Specifically, MAFA-ILS has two advantages in comparison with an approximate ILS solution by means of AFM. First, MAFA-ILS is exact. That is, the minimum of function (25) is determined precisely in the estimation process. Second, and the most important, MAFA-ILS is more computationally efficient. Figure 6 presents a fragment of the set of candidates for AFM (left) and MAFA-ILS (right). The coordinates \( {\mathbf{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x}}} \) are the sought LS solution. It is seen that in AFM one must approach this LS solution manually. Therefore, the number of candidates is very large. The point, which is closest to the minimum of function (23), is the final LS solution. The denser the grid region is, the smaller a possible error in the determination of the minimum of function (23) can be. Instead, the situation in MAFA-ILS is different. Namely, one must put at least one candidate into the ILS Voronoi cell. Afterward, in the estimation process, this point is pulled precisely to the sought LS solution (25), as shown in Fig. 6. This is why the number of points in MAFA-ILS may be much smaller than in AFM.
Numerical examples
The first numerical test is based on a three-epoch simulated observational session with an interval of 90 s from six satellites. A baseline is short, about 230 m, located in Olsztyn in Poland as shown in Fig. 7.
GPS single-frequency L1 signals (\( \lambda = 0.1903{\text{ m}} \)) are used. The theoretical coordinates of satellites and receivers are accessible at: www.gnss.5v.pl/appendix.txt. These coordinates were used for computation of theoretical DD geometric distances. Atmospheric delays are assumed to be sufficiently suppressed because of a very short baseline. The noise of simulated DD carrier phase measurements is assumed as zero-mean Gaussian noise, \( {\mathbf{e}}\left( {m \cdot n \times 1} \right) \sim N\left( {{\mathbf{0}},\sigma^{2} \cdot diag\left( {\begin{array}{*{20}c} \ldots & {{\mathbf{Q}}_{y,j} } & \ldots \\ \end{array} } \right)} \right) \), where
is the cofactor matrix for the subset of DD carrier phase measurements from the jth epoch, and \( \sigma^{2} \) is the variance of simulated DD carrier phase measurements. The Box–Muller (1958) method has been used for this purpose. The vectors of independent normal random variables were thus obtained, \( {\mathbf{e}}' \). Next, these vectors were changed into the vectors of dependent normal random variables by means of the following linear transformation
where \( chol( \cdot ) \) is an upper triangular matrix from the Cholesky decomposition. The above vectors of errors (26) were generated for three different values of DD carrier phase measurement noise: \( \sigma = 0.02 \), \( \sigma = 0.03 \), and \( \sigma = 0.04 \) cycle. The elevation-dependent weighting factor and the biases such as satellite ephemeris errors, ranging errors caused by multipath or cycle slips, have been ignored. The simulated errors (26) and some integer values of ambiguities have been added to the theoretical DD geometric distances. In this way, simulated DD carrier phase measurements have been generated. The data processing was performed with the application of the MIRLS estimation and MAFA-ILS. The MIRLS estimation was implemented by means of LAMBDA. The float position was the initial approximate position in MAFA-ILS.
Reliability
The first test was intended only for comparing success rates (SRs) of the ambiguity estimation for the MIRLS estimation and MAFA-ILS. Integer ambiguities in MAFA-ILS were calculated based on formula (11). The minimum size of the appropriate search cube in MAFA-ILS was calculated empirically based on formula (18). Additionally, the upper bound for the density of grid candidates was calculated empirically based on formula (17). These calculations gave: \( s_{\min} = 1.359{\text{ m}} \) and \( d_{\max} = 0.087{\text{ m}} \), for all three cases. However, for cognitive purposes, the SRs were calculated for different, appropriate, and inappropriate parameters (Table 1). For each variant, 10,000 simulated DD carrier phase measurements were generated, as described previously.
The results of above calculations show that the suggested solutions (17) and (18) are valid. It turns out that MAFA-ILS obtains exactly the same SR values as MIRLS. It is seen that in case of MAFA-ILS, the minimum parameters of appropriate search cube were equal: \( d = 0.6\lambda { (}0.114{\text{ m)}} \) and \( s = 12d \, (1.370{\text{ m}}) \). For these parameters, the ambiguity values (11) were always the same for all epochs and equal to the ILS ambiguity estimator (3). Furthermore, the above results also confirm the fact that the minimum values of appropriate grid parameters do not depend on a measurement noise. That is, for all three cases, the same values proved sufficient. Thus, these values depend only on the functional part of model (8). In addition, Fig. 8 presents the relationship between the SR value and the size of the search cube in MAFA-ILS. SRs obtained in the MIRLS estimation are presented in gray lines.
Precision
The results of calculations thus far have shown that the reliability of both approaches is exactly the same. However, this does not necessarily mean that both approaches give the same precision of estimated position. Therefore, the second test was carried out in order to compare the empirical precision of estimated position. Empirical covariance matrices were calculated based on the MC simulations. In addition, for cognitive purposes, analytical covariance matrices were calculated from the propagation law for formulas (4) and (10). The same computations were carried out for the float and the AFM solution. In the case of MAFA-ILS and AFM, the same size of appropriate search cube was taken, \( s = 12d \) \( (d = 0.6\lambda = 0.114{\text{ m)}} \) and \( s = 137d \) \( (d = 0.01{\text{ m)}} \), respectively. For AFM, due to a relatively long computations time, only 1000 simulations were performed, and 10,000 for the other solutions. Results are presented in Table 2.
One can see that MAFA-ILS always gives exactly the same standard deviations as the MIRLS estimation, and these values are significantly lower than in AFM. Moreover, the values of analytical standard deviations are always smaller than those of empirical ones. This obviously results from the fact that the covariance matrix calculated from the propagation law assumes that \( {\mathbf{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{a} }} = {\mathbf{a}} \) (\( SR = 100\% \)). Below, to evaluate more completely the empirical precision of both approaches, Fig. 9 presents the empirical probability density function (PDF) for the L2-norm vector of the true errors of estimated position, \( \left\| {{\mathbf{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\smile}$}}{x} }} - {\mathbf{x}}} \right\| \), separately for a success and for a failure. Despite the fact that the empirical standard deviations are exactly the same for both approaches (Table 2), the plot of empirical PDFs differs as shown in Fig. 9. For example, small differences are visible in the enlarged PDF fragments for success. This could be explained by the fact that both approaches are based on different values of approximate coordinates. Namely, MAFA-ILS is based on approximate coordinates from the ILS Voronoi cell, while the MIRLS estimation is based on different approximate coordinates, for example, from the single-point positioning. However, the above differences are not significant, and it is possible to suggest the thesis that the estimated position has the same probabilistic properties in both approaches.
Table 2 also presents the information on the computational efficiency of MAFA-ILS, AFM, and MIRLS estimation. In fact, the computation time in MAFA-ILS is a few hundred times shorter than that for AFM, but it is ten or so times longer than in the MIRLS estimation.
The second numerical test is based on three-epoch sessions of the GPS real measurements with an interval of 90 s. The baseline WRLG–WT21 is very short, about 60 m, located in Wettzell in Germany. Data were obtained from BKG (Bundesamt für Kartographie und Geodäsie) GNSS Data Center for the time interval 12:00:00–13:00:00 UTC on July 19, 2017. The coordinates of the satellites for the epochs of signal emission were interpolated based on precise orbits SP3, taking into account all required time corrections. An elevation mask angle was equal to 10°. The observations did not contain cycle slips. Atmospheric delays were assumed to be sufficiently suppressed because of a very short baseline. The measurements on the L1 frequency were processed with the application of the MIRLS estimation and MAFA-ILS based on ten, three-epoch sessions. The first session concerns the epochs 12:00:00, 12:01:30 and 12:03:00, the second the epochs 12:03:30, 12:05:00 and 12:06:30, etc. The parameters of the search cube in MAFA-ILS were calculated empirically based on formulas (17) and (18). The float position was calculated based on L1 carrier phase measurements, and contrary to the previous simulation test, additionally C/A pseudoranges. As a consequence, the search cube was smaller than previously. The computation time was equal to and about 0.02 and 0.05 s in the MIRLS estimation and MAFA-ILS, respectively. Figure 10 presents the estimated coordinates of the WT21 point. The ‘true’ position was obtained from the float solution of a 120-epoch session (12:00:00–13:00:00; interval 30 s). One can see that the results of above field experiment are consistent with the results of the previous simulated experiment. Both approaches give exactly the same estimated positions, which lead us to the previous conclusion that MAFA-ILS can be treated as MIRLS estimation.
Summary and conclusions
We discussed and pretested a non-conventional approach to the MIRLS estimation which may open the new field of research. Contrary to the conventional MIRLS estimation, the search for a solution occurs in a constant three-dimensional coordinate domain and is much faster than in AFM. The performed experiments have shown that the suggested approach may be considered in (near) real-time GNSS positioning. Despite that, the search for a solution is sometimes about ten times slower than in the conventional MIRLS estimation. However, in this research, the center of the search region was based on the float position, which in terms of computational efficiency is suboptimal. Second, the problem of the optimization of the search procedure was not the objective of this research, and a very trivial procedure has been used. The optimization of the search procedure in the coordinate domain is a separate, currently open problem. In Cellmer et al. (2018), we propose to use a grid ellipsoid instead of a grid cube and a cascade processing with a wide-lane combination in a first step. These solutions can significantly improve the computational efficiency.
References
Box GEP, Muller ME (1958) A note on the generation of random normal deviates. Ann Math Stat 29(2):610–611
Brack A (2015) On reliable data-driven partial GNSS ambiguity resolution. GPS Solut 19(3):411–422
Brack A (2016) Optimal estimation of a subset of integers with application to GNSS. Artif Satell 51(4):123–134
Cellmer S (2012) A graphic representation of the necessary condition for the MAFA method. IEEE Trans Geosci Remote 50(2):482–488
Cellmer S (2013) Search procedure for improving modified ambiguity function approach. Surv Rev 45(332):380–385
Cellmer S, Wielgosz P, Rzepecka Z (2010) Modified ambiguity function approach for GPS carrier phase positioning. J Geod 84(4):264–275
Cellmer S, Nowel K, Kwaśniak D (2018) The new search method in precise GNSS positioning. IEEE Trans Aerosp Electron Syst. https://doi.org/10.1109/TAES.2017.2760578
Chang X, Yang X, Zhou T (2005) MLAMBDA: a modified LAMBDA algorithm for integer least-squares estimation. J Geod 79(9):552–565
Chen D, Lachapelle G (1995) A comparison of the FASF and least-squares search algorithms for on-the-fly ambiguity resolution. Navigation 42(2):371–390
Counselman C, Gourevitch S (1981) Miniature interferometer terminals for earth surveying: ambiguity and multipath with global positioning system. IEEE Trans Geosci Remote 19(4):244–252
Frei E, Beutler G (1990) Rapid static positioning based on the fast ambiguity resolution approach FARA: theory and first results. Manuscr Geodaet 15(6):325–356
Jazaeri S, Amiri-Simkooei AR, Sharifi M (2014) On lattice reduction algorithms for solving weighted integer least squares problems: comparative study. GPS Solut 18(1):105–114
Kim D, Langley RB (2000a) GPS ambiguity resolution and validation: methodologies, trends and issues. In: Proceedings of the GNSS workshop international symposium on GPS/GNSS, Seoul, Korea, November–December, vol 30(2.12), pp 1–9
Kim D, Langley RB (2000b) A search space optimization technique for improving ambiguity resolution and computational efficiency. Earth Planets Space 52(10):807–812
Leick A, Rapoport L, Tatarnikov D (2015) GPS satellite surveying, 4th edn. Wiley, New York
Remondi B (1984) Using the global positioning system (GPS) phase observable for relative geodesy: modelling, processing and results. Dissertation, University of Texas at Austin
Remondi B (1990) Pseudo-kinematic GPS results using the ambiguity function method. Navigation 38(1):17–36
Teunissen PJG (1995) The least-squares ambiguity decorrelation adjustment: a method for fast GPS integer ambiguity estimation. J Geod 70(1):65–82
Teunissen PJG (1999a) The probability distribution of the GPS baseline for a class of integer ambiguity estimators. J Geod 73(5):275–284
Teunissen PJG (1999b) An optimality property of the integer least-squares estimator. J Geod 73(11):587–593
Teunissen PJG (2002) A new class of GNSS ambiguity estimators. Artif Satell 37(4):111–120
Teunissen PJG (2003) Integer aperture GNSS ambiguity resolution. Artif Satell 38(3):79–88
Teunissen PJG, Joosten P, Tiberius CCJM (1999) Geometry-free ambiguity success rates in case of partial fixing. In: Proceedings of the ION GPS 1999, Institute of Navigation, San Diego, USA, September 14–17, pp 201–207
Verhagen S, Teunissen PJG (2013) The ratio test for future GNSS ambiguity resolution. GPS Solut 17(4):535–548
Verhagen S, Teunissen PJG, van der Marel H, Li B (2011) GNSS ambiguity resolution: which subset to fix. In: Proceedings of IGNSS symposium, Sydney, Australia, November 15–17
Xu P (2001) Random simulation and GPS decorrelation. J Geod 75(7):408–423
Xu P (2006) Voronoi cells, probabilistic bounds, and hypothesis testing in mixed integer linear models. IEEE Trans Inform Theory 52(7):3122–3138
Xu P (2012) Parallel Cholesky-based reduction for the weighted integer least squares problem. J Geod 86(1):35–52
Zhou Y (2011) A new practical approach to GNSS high-dimensional ambiguity decorrelation. GPS Solut 15(4):325–331
Acknowledgements
We would like to thank anonymous reviewers for their time and constructive comments, which improved the quality of this paper. This research was supported by a grant with Agreement Number UMO-2014/13/B/ST10/02547 from Polish National Center of Science.
Author information
Authors and Affiliations
Corresponding author
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
About this article
Cite this article
Nowel, K., Cellmer, S. & Kwaśniak, D. Mixed integer–real least squares estimation for precise GNSS positioning using a modified ambiguity function approach. GPS Solut 22, 31 (2018). https://doi.org/10.1007/s10291-017-0694-6
Received:
Accepted:
Published:
DOI: https://doi.org/10.1007/s10291-017-0694-6