Solving the Multilateration Problem without Iteration

The process of positioning, using only distances from control stations, is called trilateration (or multilateration if the problem is over-determined). The observation equation is Pythagoras’s formula, in terms of the summed squares of coordinate differences and, thus, is nonlinear. There is one observation equation for each control station, at a minimum, which produces a system of simultaneous equations to solve. Over-determined nonlinear systems of simultaneous equations are typically solved using iterative least squares after forming the system as a truncated Taylor’s series, omitting the nonlinear terms. This paper provides a linearization of the observation equation that is not a truncated infinite series—it is exact—and, thus, is solved exactly, with full rigor, without iteration and, thus, without the need of first providing approximate coordinates to seed the iteration. However, there is a cost of requiring an additional observation beyond that required by the non-linear approach. The examples and terminology come from terrestrial land surveying, but the method is fully general: it works for, say, radio beacon positioning, as well. The approach can use slope distances directly, which avoids the possible errors introduced by atmospheric refraction into the zenith-angle observations needed to provide horizontal distances. The formulas are derived for two- and three-dimensional cases and illustrated with an example using total-station and global navigation satellite system (GNSS) data.


Introduction
A station's coordinates can be determined from distances-with no angles-using a technique known as trilateration, or multilateration if there are redundant observations. The U.S. Bureau of Land Management [1] defines trilateration as, "A method of determining horizontal positions by measuring the lengths of triangle sides, usually with the use of electronic instruments." This definition predates global navigation satellite systems (GNSSs) and is a little dated. The method is used extensively nowadays because it is how GNSSs are used to determine the coordinates of a receiver [2,3], and such positions are inherently three dimensional. For terrestrial surveying, slope distances observed with an electronic distance meter (EDM) can be used, likewise producing three-dimensional positions. The trilateration problem can be stated: given the coordinates of at least three/four control marks and the distances from the control marks to a target mark, determine the horizontal/3-d coordinates of a target mark. This is similar to the resection problem, but resection uses directions instead of distances [4].
Pythagoras's formula is the observation equation, written with the lengths of the shorter sides are given as differences in Cartesian coordinates, Equation (3). There is one equation for each control station, so there is a set of simultaneous equations solve. The system has an analytical solution-there are three equations and three unknowns in a quadratic-however, the exact solution is somewhat unwieldy and, worse, cannot handle an over-determined system, having more observations than the minimum of three. The typical solution strategy is to linearize the system by expanding the equations in a Taylor's series and retain only the first-order terms, then solve by iterative least squares. However, as will be shown, taking the difference of the observation equations linearizes the system without any approximations, which allows the solution to be found without iteration.
Multilateration uses mark-to-mark slope distances. If the height of the prism target above the unknown station, hr, equals the height of the optical center of the telescope in a total station, hi, then the observed slope distance s equals mark-to-mark slope distance s. Otherwise, let ζ denote the observed zenith angle. Then the mark-to-mark zenith angle is [5] cot ζ = csc ζ hi − hr s + cot ζ (1) and the mark-to-mark slope distance is Hereafter, all slope distance will be assumed to be mark-to-mark. Suppose we are to determine the coordinates of station i = (e i , n i , u i ) given the coordinates of a set of control stations C = e j , n j , u j , . . . (e m , n m , u m ) , and |C| = n, i.e., there are n control stations. The slope distance from station i to station j, s ij , is related to the stations' Cartesian coordinates by [6].
Similarly, the slope distance from station i to station k, s ik , is related to the stations' coordinates by Expanding the quadratic terms, showing only the easting terms for brevity, gives Taking s 2 ij − s 2 ik eliminates the e 2 i , n 2 i , and u 2 i terms, leaving s 2 ij − s 2 ik = e 2 j + 2e j e i + e 2 i + · · · − e 2 k + 2e k e i + e 2 i + · · · = e 2 j − e 2 k + 2e i e j − e k + · · · (6) Putting constants on the right and factoring on the left, gives The only unknown terms are e i , n i , and u i . If station i has been observed from stations j . . . m, then these equations can be written as a matrix expression as Or So, if M is 3 × 3, If M has more than three rows, the system is over-determined, and it is solved with a pseudo-inverse in the usual way: Matrix M has the differences of pairs of control coordinates for its elements. There are n(n − 1)/2 ways to pick two items from a set of n things, and each of these could be a row in M. That would result in M having quite a few rows producing, supposedly, a lot of built-in redundancy. However, the GNSS-positioning concept of "independent baselines" applies here: many of the pairs would be trivial inverses of others, and no real additional information is added. Therefore, like with GNSS, pick one station to be a common station and form all differences with the common station. Station j is the common station in Equation (8). The constant term k in Equation (9) can be written as k is the square of distance from station k to the origin of the coordinate system; similar for station j. This suggests two things. First, picking the control station closest to the centroid of all the control stations to be the common differencing station minimizes |K| 2 − |J| 2 . Second, it is sensible to use localized coordinates to keep the coordinates themselves as small as possible so that the quadratic terms do not suffer undue truncation and round-off errors. For example, if station i is set to be a local origin, then all other stations' coordinates become e j − e i , n j − n i , u j − u i .
Our method is similar to that developed by Navidi et al. [7]. Their approach depends on a special control station that they call the reference point that is required to be at the centroid of the other control stations. They then add and subtract the coordinates of the reference point from those of the other control stations to form a third type of distance that, when added to the distances from the other control stations to the unknown point, causes most of the nonlinear terms to cancel. Placing the reference point at the centroid is highly auspicious because it causes the observation matrix to become orthogonal to the parameter vector, which eliminates the remaining nonlinear term. Ignoring numerical stability issues, our method and theirs must be equivalent: they are both least squares estimators, and, since least squares is a best linear unbiased estimator, any two different approaches to the same problem using least squares must be equivalent. However, the problem formulations can differ and that can lead to implementation issues. As just stated above, we recommend that the common differencing station be near the centroid, but it is not a requirement. The method by Navidi et al. requires that the reference point be at the centroid, which could be onerous depending on how difficult it is to stake out that position, assuming that it is possible at all.
The positioning equation for GNSS uses multilateration but it includes additional terms [8], p. 410: where x, y, and z are geocentric Cartesian coordinates in the reference frame of the broadcast ephemeris, c is the speed of light in a vacuum, b u is the unknown time offset of the receiver's clock to GPS time (called a time bias), and pi is a "catch-all" error term for ionospheric and tropospheric delays, multipath, and c. Unfortunately, our scheme does not apply here because squaring the right-hand side of the equation mixes the terms together so that the quadratic terms do not cancel after differencing. Therefore, the equation is linearized by a truncated Taylor's series and solved with iteration.

Stochastic Model
One can reasonably assume that the slope-distance observations are independent and normally distributed as s j ∼ N(µ j , σ 2 j ) where µ j is the true distance, and σ j is the EDM standard deviation given by σ j = c + s j /1000 × ppm, where c is the constant term, and ppm is the "parts per million" millimeters per kilometer term. This model ignores setup uncertainty at the instrument and at the target. As such, the n × n covariance matrix of the observations Q = diag σ 2 1 , σ 2 2 , . . . , σ 2 n is a diagonal matrix. The inverse of the covariance matrix W = Q −1 = diag 1/σ 2 1 , 1/σ 2 2 , . . . , 1/σ 2 n can be used as a weight matrix, and is applied as which is solved by The matrix M t W −1 M −1 is called a cofactor matrix, and it is proportional to the estimated covariance of the estimations. The constant of proportionality is called the variance of unit weight, or just "unit variance," and the estimated covariance of the The unit variance can itself be estimated as follows, [9,10].
where, n = |k| is the number of observations and u = |x| is the number of coordinates to estimate (i.e., three in this case).

Coordinate Systems
The formulas require Cartesian coordinates for the control stations. There are three choices: geocentric Cartesian (called XYZ), topocentric Cartesian (called ENU), or grid coordinates. Grid coordinates require additional reductions be applied to the slope distances, which is thoroughly covered in standard surveying texts and an unnecessary distraction here. In fact, XYZ coordinates are well suited for this application, at least in principle.

Example: More Control than the Minimum
Hereafter, all linear quantities will be given in meters. In this example, we pretend to confirm the coordinates for the CTMA CORS ARP NAD 83(2011) coordinates (https://www.ngs. noaa.gov/cgi-bin/ds_cors.prl?CorsSelected=\T1\textbar{}CTMA&CorsTypeSelected=Arp, accessed on 28 June 2021), which are x = 1, 456, 379.711, y = −4, 539, 030.822, z = 4, 223, 420.343. Suppose there are five control stations available whose XYZ coordinates were determined using GNSS positioning as tabulated in Table 1. In fact, these stations' coordinates were generated in Mathematica™ using random-number generators, such that the horizontal coordinates follow uniform distributions ±300 m from CTMA and the vertical coordinates follow a uniform distribution ±5 m from CTMA. Inversing the control coordinates with CTMA produces exact slope distances that a perfect EDM would observe. These slope distances from the five control stations to the "unknown station" CTMA are 44.947, 22.106, 183.812, 82.023, and 355.566, to three significant digits. First, use the coordinates as given. The system of equations is Using Mathematica, the solution works out to x = 1, 456, 379.711, y = −4, 539, 030.822, z = 4, 223, 420.343. These coordinates differ from CTMA by −9.8 × 10 −6 , 0.00003, −0.00003, which is good but it is not even machine precision. The issue is the large values in k, which involve the squares of the coordinate values, the right side of Equation (8). The solution can be improved by localizing the coordinates, which makes them much smaller numbers and, thus, improves numerical stability. For example, subtract station A's coordinates from the other control stations, resulting in the coordinates in Table 2.

ENU Coordinates
If the control stations' coordinates had been determined by radial surveying from station A using a perfect total station, we would have the control coordinates shown in Table 3. ENU coordinates are XYZ coordinates that have been translated and rotated-no scaling or stretching-which means the distances and angles between the stations are preserved exactly in the two systems. Therefore, the ENU solution must, in principle, be identical to the XYZ solution, and would be exactly identical if exact arithmetic were available. The implication here is that, in principle, there is no need to reduce slope distances to horizontal distances.
Good strength-of-figure results from surrounding the unknown station with observations from all directions. It is unlikely, however, that the unknown station can be observed from far below or above. There might be a situation where an observation could be, say, made from a rooftop, but it will usually happen that the control stations are generally around the same height. The concepts of dilutions of precision pertain here. Borrowing from GNSS [11], the dilutions of precisions are various combinations of the diagonal elements of the estimated position's covariance matrix, being the variances of the estimated coordinates. As with GNSS positioning, and fundamentally for the same reason, the horizontal dilution of precision will be better (lower) than the vertical dilution of precision because, for GNSS, it is impossible to observe the satellites (a.k.a. space vehicles [SVs]) below the horizon so the positioning cannot be controlled from below. If all the Up-coordinates are the same so that their differences are zero, then matrix M =   a b 0 c d 0 e f 0   has a row echelon form of   1 0 0 0 1 0 0 0 0   , which shows the matrix cannot be inverted. More generally, co-planar control stations produce an M matrix that cannot be inverted. However, M can be inverted, in principle, if the stations are not exactly coplanar, but M can become ill-conditioned [10]. For example, replacing the Up-coordinates in the previous example with normally distributed random numbers u ∼ N(0, 0.01) produced a M T W M matrix with elements differing by ten orders of magnitude, resulting in a position estimate of (0.024, 0.023, −503.181) instead of (0,0,0). The accuracies of the horizontal coordinates are compromised by the ill-conditioning brought about by the vertical coordinates. Not including the vertical coordinates mitigates this, so we now re-work the problem for horizontal positioning.

Horizontal Positioning
The horizontal-only equations come from dropping the Up-coordinates in Equation (8) and replacing the slope distances with horizontal distances, s = s cos(ζ), where ζ is the mark-to-mark zenith angle. Equation (8) Solving Equation (17) with horizontal distances from the previous examples produces a position estimate of −1.5 × 10 −14 , 1.0 × 10 −14 , which is around machine precision for double-precision floating point arithmetic.

Fieldwork Example
Observations were collected on the campus of New Mexico State University in May 2020. Five stations were positioned using a static GNSS survey lasting for at least 25 min, with a 1-s sampling interval, using Topcon HiPer Lite+ GNSS receivers. The specifications of the receivers, as reported by the manufacturer, are horizontal accuracies of ±3 mm ± 0.5 ppm and vertical accuracies of ±5 mm ± 0.5 ppm for GNSS static surveys. Data were then processed using the Online Positioning User Service (OPUS) tool provided by NGS (https://geodesy.noaa.gov/OPUS/, accessed on 28 June 2021), [12]. Seven CORSs provided control, being an average of around 100 km away (DK7580, DM3997, DJ8945, DG7417, DG6517, DG5392, DG5481). The average RMS values for all coordinates was less than 3 cm, as reported by OPUS. We also measured the mark-to-mark distances using a Leica TS02 total station and accompanied Leica reflector, both erected atop a slip-leg tripod and a tribrach with an optical plummet. Leica reports the accuracy of the EDM as 1.5 mm + 2 ppm. Each distance was measured three times in both the left and right faces. The height of instrument and the height of the reflector were measured and input to the total station. We then documented the observed horizontal and slope distances ( Table 4).
The International Terrestrial Reference Frame (ITRF) 2014 coordinates are tabulated in Table 5. Stations A, B, C, D, and U were positioned with OPUS, and station CG is the average of stations A, B, C, D. Station U is the check station, whose coordinates are to be determined using stations A, B, C, D. The XYZ coordinates were converted to geodetic longitude, latitude, and height (λ, φ, h), and then to local topocentric ENU using CG as the origin. The direct-inverse observations where averaged to mitigate collimation error, [13] and then all corrected pairs averaged to a single value. This is not necessary: M and k can be formed from the individual collimation-corrected averages but using the single ensemble average values yielded nearly identical estimates and makes the examples below easier to read. In practice, it might be preferable to use the individual collimation-corrected averages for blunder detection. • These coordinates differ from the check coordinates by (1.773, 5.405, −3.606 ), which is 6.735 m in error (3-d).

2.
Estimate station U's coordinates using slope distances and XYZ coordinates localized to their centroid, columns x , y , z in Table 5. • The XYZ coordinates are recovered by adding the coordinates of CG to the estimated coordinates.

•
The estimated coordinates differ from the check coordinates by (1.757, 5.356, −3.573), which is 6.674 m in error (3-d), slightly worse than using the XYZ coordinates directly.

•
The matrix M is identical for both "absolute" and these relative coordinates because M is formed from the differences in coordinates.

•
The k vector's values are much smaller for the relative coordinates because the distances from the stations to the local origin (CG) are far shorter than the distances of the XYZ stations to their origin (Earth's center of mass). These smaller values will generally lead to better numerical precision, which should have a good influence on improving the estimates, but that was not the case here.

•
The XYZ coordinates yield terrible results because, although these are geometrically 3-d coordinates, the stations are nearly coplanar, topologically 2-d. As discussed above, this produces very poor strength-of-figure for the vertical coordinate, but there is no notion of the vertical in an XYZ coordinate system so the inherent weakness in the vertical direction bleeds into all three coordinates.

3.
Estimate station U's coordinates using topocentric Cartesian ENU coordinates using station CG as the ENU coordinate system's origin. See the e, n, u columns in Table 5. The numerical instability due to the ill-conditioned vertical coordinate has been removed so the estimates are very good. Notice that these EN coordinates have been estimated better than their ENU counterparts, because the deleterious effect of the Up coordinate has no effect here. • Notice that M is almost an identity matrix, which means that k has almost solved the problem itself.

Discussion and Conclusions
A multilateration method for terrestrial surveying was presented using differenced observation equations, which eliminated the quadratic terms in the unknown coordinates and, thus, linearized the system to be solved. Consequently, the position can be determined with a pseudo-inverse in the usual way, an operation that is well within reach of any spreadsheet's computational capabilities, and still provides full statistical rigor. The linearization does have a cost: using differenced observations means one additional observation set than is required by the non-linear equivalent (four instead of three).
Understanding the geometrical layout of the survey is difficult when looking at the XYZ coordinates because they impart no notion of cardinal directionality or the vertical. The ENU coordinates readily reveal that the stations are nearly at the same height and occupy the four quadrants around the check station U. Slope distances can be used with XYZ and ENU coordinates alike, and slope distances are largely unaffected by atmospheric refraction and vertical-angle pointing errors that affect horizontal distances. However, as the examples show, poor strength-of-figure can completely overwhelm these benefits when the control stations are nearly coplanar.
Standard practice is to observe all distances at least twice in independent observations. A typical protocol would be to observe in different setups, not merely observing in the direct and inverted positions of the telescope.
M T M is invertible, in theory, so long as the control stations are not coplanar, but M T M becomes ill conditioned if the stations are bunched together or if they are nearly coplanar. If M T M is ill conditioned, all three coordinates will have less accurate estimations because the estimates are correlated. However, as demonstrated in the example, the Up-coordinate can be especially erroneous, but perhaps sometimes not so much as to be blatantly incorrect, which is a very insidious situation. In contrast, there is no numerical stability concerns if solving for just the horizontal coordinates, provided that the control stations surround the unknown station. The best evidence for whether the Up-coordinate is trustworthy is, unsurprisingly, its standard deviation (and common sense). The covariance matrix for the above ENU example works out to be 10 −6   2.5 −1.4 6.8 −1.4 1.1 −17 6.8 −17 769   so the standard deviations of the E-, N-, and U-coordinates are 1.6 mm, 1.1 mm, and 27.7 mm. As with GNSS positioning, the horizontal coordinates should have similar standard deviations and the Up-coordinate's standard deviation should be larger. However, for GNSS, the Up-coordinate is expected to be around, say, three times larger, and here, the Up-coordinate is around 20 times larger. This is because, for GNSS, one expects there to be SVs spread across the whole sky at a broad span of zenith angles, whereas for terrestrial surveying one expects the control stations to be nearly at the same height, creating a scenario equivalent for GNSS positioning in which all the SVs are within, say, 5 • in elevation angle from one another. Thus, in general, the Up-coordinate cannot be expected to have as good precision as the horizontal coordinates in most cases. There is no limit, per se, to the number of stations that can be positioned simultaneously with this method, but every unknown station must tie to enough control stations to be positioned individually. In contrast, a typical traverse will run through "interior" stations that have no, or insufficient, direct ties to control stations, such as stations 2 and 4 in Figure 1. That is not allowed here. The reason is in the M matrix: if any of M's elements (coordinate differences) are not unknown a priori, then M will contain variable elements, and the system becomes nonlinear. This is true for k, too: k will not be a vector of numbers only. This is not the case when solving the multilateration problem with a nonlinear system. The nonlinear solution requires a priori estimates for all the unknown coordinates, so the Jacobian matrix's terms contain only numbers and, thus, the nonlinearity vanishes in the calculation. Nevertheless, the traverse in Figure 1 can be solved in a bootstrapping fashion by positioning stations 1 and 5, and then using those to control 2 and 4. However, if station g were unknown, nothing could be done.
Geomatics 2021, 1, FOR PEER REVIEW 10 and the system becomes nonlinear. This is true for k, too: k will not be a vector of numbers only. This is not the case when solving the multilateration problem with a nonlinear system. The nonlinear solution requires a priori estimates for all the unknown coordinates, so the Jacobian matrix's terms contain only numbers and, thus, the nonlinearity vanishes in the calculation. Nevertheless, the traverse in Figure 1 can be solved in a bootstrapping fashion by positioning stations 1 and 5, and then using those to control 2 and 4. However, if station g were unknown, nothing could be done.