GPS Coordination

Question

I have N GPS coordinates with N distances given in an unknown position that I want to determine.

My first approach was to use only three points and trilateration, as described here . This approach was already quite accurate (best error ~ 5 km), but I would like to improve this and increase reliability. Since these distances are not very accurate, I thought about using multiple measurements and multilateration. However, it turned out that this approach is much less accurate (the best error is ~ 100 km), although I provide more than 3 points / distances (verified to 6), and now I ask if anyone has an idea what I could have done wrong.

In short, my approach to multilateration is as follows:

• Convert all coordinates to ECEF
• Build the matrix as described in Eq.7 on wikipedia
• Use SVD to find a minimizer
• Since the solution is only to scale, I use the root search approach to determine normalization, so that the coordinates converted back to LLA lead to a height of 0 (my initial assumption is that all coordinates are at zero height)
• Convert back to LLA

The LLA / ECEF conversion is checked twice and correctly. Step 2 and 3 I checked with the Euclidean coordinates (and exact distances) and look right. I came up with step 4 myself, I don’t know how good this approach is, so suggestions are welcome.

+++ UPDATE

I put together a sample code in python to illustrate the problem with some truth. Trilateration reaches 400 m, and multilateration fluctuates here by 10-130 km. Because of the length, I put it on ideone

Answer 1

In the end, I figured it out myself - or at least greatly improved accuracy.

The approach described in wikipedia (Eq.7) does not seem to be very suitable for this application, but in this case it is much simpler.

Given Wikipedia equation 6, we can simplify it: R_0 can be guessed as the radius of the Earth, since the origin of the ECEF coordinates lies in the center of the Earth. Therefore, there is no need to shift everything to one point of origin, and we can use all the equations of N

In python with a P array of ECEF coordinates and dists distance to these points all boils down to

 R = 6378137 # Earth radius in meters A = [] for m in range(0,len(P)): x = P[m][0] y = P[m][1] z = P[m][2] Am = -2*x Bm = -2*y Cm = -2*z Dm = R*R + (pow(x,2)+pow(y,2)+pow(z,2)) - pow(dists[m],2) A += [[Am,Bm,Cm,Dm]] # Solve using SVD A = numpy.array(A) (_,_,v) = numpy.linalg.svd(A) # Get the minimizer w = v[3,:] w /= w[3] # Resulting position in ECEF 

With this approach, what I described as step 4 is no longer required. In fact, this even makes the solution worse.

Now the accuracy ranges from 2 to 275 m - in most cases better than the “optimal” trilateration with an error of 464 m.

Answer 2

Some comments:

1) You have already checked a few steps against the exact answers. I suggest you create toy problems with a known amount of random noise added to the observations. Since you know the correct answer in this case, you can see what happens with the propagation of errors. If your method works well, but poorly on real data, you might think of terrible behavior in real life, for example, that one or more of the distances is seriously wrong.

2) I don’t know why your solution only scales, since the basic data scales correctly - if I went there with taut ropes and tied them to fixed points, there would be no ambiguity. When you use SVD to solve equations (7), do you do something like www.cse.unr.edu/~bebis/MathMethods/SVD/lecture.pdf to get the least squares solution? That should give you x, y and z without ambiguity.

3) I am not at all sure how observational errors work (7). Firstly, I do not like all units. It might be worth writing the equation for the sum of the squared differences between the measured distances and the calculated distances given by x, y, z for an unknown position, and then minimizing this for x, y, z. The Wikipedia article discards this approach because of its cost, but may give you a more accurate answer, and calculating and comparing this answer may tell you something, even if you cannot use this method in practice.

Answer 3

I followed the @zerm code shown above and it worked quite well (yellow is a calculated point from 3 towers). The results are shown in a fragment of Folium. Multilateration using linalg.SVD

However, when I followed the same algorithm with the changes suggested by @mcdowella (No. 2) using the least squares solutions of the MxN system, the results are much better. Least Square Multiple Literature MxN

Here is the corrected code:

 A = [] b = [] for m in range(0,len(P)): x = P[m][0] y = P[m][1] z = P[m][2] Am = 2*x Bm = 2*y Cm = 2*z Dm = R*R + (pow(x,2)+pow(y,2)+pow(z,2)) - pow(dists[m],2) A += [[Am,Bm,Cm]] b += [[Dm]] # Solve using Least Squares of an MxN System # A*x = b --> x = (ATA)_inv.AT.b = A+.b A = np.array(A) b = np.array(b) AT = AT ATA = np.matmul(AT,A) ATA_inv = np.linalg.inv(ATA) Aplus = np.matmul(ATA_inv,AT) x = np.matmul(Aplus,b) # convert back to lat/long from ECEF # convert to degrees lat = math.degrees(math.asin(x[2] / R)) lon = math.degrees(math.atan2(x[1],x[0])) 

I am still studying other multilateration methods, but this post really allowed me to understand the basics of N-Points MLAT. Thanks!

Answer 4

To improve the accepted answer , one way to improve the SVD solution is to take into account latitude-based changes in the radius of the earth; This especially affects elevation estimates, but also has some knock effects on latitude and longitude. A “simple” solution would be to use the average value for R , which according to Wikipedia is 6371008.8 m, and not 6378137 m.

A more accurate estimate would be to adjust R for latitude:

 def EarthRadiusAtLatitude(lat): rlat = np.deg2rad(lat) a = np.float64(6378137.0) b = np.float64(6356752.3) rad = np.sqrt(((a*a*np.cos(rlat))**2 + (b*b*np.sin(rlat))**2) / ((a*np.cos(rlat))**2 + (b*np.sin(rlat))**2)) return rad 

Then set R depending on the latitude of one of the starting points. Or, if you have a large variation in latitude, you can calculate the SVD based on the R score and use the latitude of the solution to solve using the closer R estimate.

After completing this setup, in my experiments with both constructed examples and “real world” data based on eNodeB LTE latency values, the SVD solution is usually within one second of latitude and longitude, except in some degenerate cases, which is fairly comparable to the solution based on iterative optimization (i.e. minimization of residual distances).