¿Trilateración usando 3 puntos de latitud / longitud y 3 distancias?

34

Quiero encontrar una ubicación objetivo desconocida (coordenadas de latitud y longitud). Hay 3 puntos conocidos (pares de coordenadas de latitud y longitud) y para cada punto una distancia en kilómetros a la ubicación del objetivo. ¿Cómo puedo calcular las coordenadas de la ubicación de destino?

Por ejemplo, supongamos que tengo los siguientes puntos de datos

37.418436,-121.963477   0.265710701754km
37.417243,-121.961889   0.234592423446km
37.418692,-121.960194   0.0548954278262km

Lo que me gustaría es cuál es la matemática para una función que toma eso como entrada y devuelve 37.417959, -121.961954 como salida.

Entiendo cómo calcular la distancia entre dos puntos, de http://www.movable-type.co.uk/scripts/latlong.html Entiendo el principio general de que con tres círculos como estos obtienes exactamente un punto de superposición. De lo que estoy confuso es de las matemáticas necesarias para calcular ese punto con esta entrada.

Sin sombrero
fuente
Aquí hay una página que lo guía a través de las matemáticas para encontrar el centro de tres coordenadas. Quizás podría ayudar de alguna manera. < mathforum.org/library/drmath/view/68373.html >
Jon Bringhurst el
1
¿Esto necesita estar en la esfera / esferoide, o está bien un algoritmo plano?
fmark
1
No puedo darte la respuesta, pero creo que puedo orientarte en la dirección correcta. Tres coordenadas = tres puntos centrales. Tres distancias = tres círculos. Dos círculos que se cruzan pueden tener la posibilidad de ninguna / una / dos soluciones. Tres círculos pueden tener ninguno / uno / o un área como su solución. Obtenga la fórmula del círculo para los tres círculos y resuélvala con Sistemas de ecuaciones / álgebra.
CrazyEnigma
En realidad, ni siquiera necesitas sistemas para resolver este. Hay una o dos posibilidades, pero como tiene un valor de distancia, puede separar la respuesta correcta.
George Silva
1
+1 Esta es una buena pregunta. Al principio pensé que se podría encontrar fácilmente una solución con google, pero aparentemente no. Quizás el problema podría ser más general: dado que N puntos con no solo una distancia sino también un margen de error, encuentre la elipse de confianza.
Kirk Kuykendall

Respuestas:

34

Después de mirar un poco en Wikipedia y la misma pregunta / respuesta en StackOverflow , pensé que podría darle una puñalada e intentar llenar los vacíos.

En primer lugar, no estoy seguro de dónde obtuviste la salida, pero parece estar mal. Tracé los puntos en ArcMap, los amortigué a las distancias especificadas, corrí intersectando en los búferes y luego capturé el vértice de la intersección para obtener las soluciones. Su salida propuesta es el punto en verde. Calculé el valor en el cuadro de llamada, que es aproximadamente 3 metros de lo que ArcMap dio para la solución derivada de la intersección.

texto alternativo

La matemática en la página de Wikipedia no es tan mala, solo necesita convertir sus coordenadas geodésicas a la ECEF cartesiana, que se puede encontrar aquí . los términos a / x + h pueden reemplazarse por el radio de la esfera autálica, si no está utilizando un elipsoide.

Probablemente lo más fácil es darle un código bien documentado (?), Así que aquí está en Python

import math
import numpy

#assuming elevation = 0
earthR = 6371
LatA = 37.418436
LonA = -121.963477
DistA = 0.265710701754
LatB = 37.417243
LonB = -121.961889
DistB = 0.234592423446
LatC = 37.418692
LonC = -121.960194
DistC = 0.0548954278262

#using authalic sphere
#if using an ellipsoid this step is slightly different
#Convert geodetic Lat/Long to ECEF xyz
#   1. Convert Lat/Long to radians
#   2. Convert Lat/Long(radians) to ECEF
xA = earthR *(math.cos(math.radians(LatA)) * math.cos(math.radians(LonA)))
yA = earthR *(math.cos(math.radians(LatA)) * math.sin(math.radians(LonA)))
zA = earthR *(math.sin(math.radians(LatA)))

xB = earthR *(math.cos(math.radians(LatB)) * math.cos(math.radians(LonB)))
yB = earthR *(math.cos(math.radians(LatB)) * math.sin(math.radians(LonB)))
zB = earthR *(math.sin(math.radians(LatB)))

xC = earthR *(math.cos(math.radians(LatC)) * math.cos(math.radians(LonC)))
yC = earthR *(math.cos(math.radians(LatC)) * math.sin(math.radians(LonC)))
zC = earthR *(math.sin(math.radians(LatC)))

P1 = numpy.array([xA, yA, zA])
P2 = numpy.array([xB, yB, zB])
P3 = numpy.array([xC, yC, zC])

#from wikipedia
#transform to get circle 1 at origin
#transform to get circle 2 on x axis
ex = (P2 - P1)/(numpy.linalg.norm(P2 - P1))
i = numpy.dot(ex, P3 - P1)
ey = (P3 - P1 - i*ex)/(numpy.linalg.norm(P3 - P1 - i*ex))
ez = numpy.cross(ex,ey)
d = numpy.linalg.norm(P2 - P1)
j = numpy.dot(ey, P3 - P1)

#from wikipedia
#plug and chug using above values
x = (pow(DistA,2) - pow(DistB,2) + pow(d,2))/(2*d)
y = ((pow(DistA,2) - pow(DistC,2) + pow(i,2) + pow(j,2))/(2*j)) - ((i/j)*x)

# only one case shown here
z = numpy.sqrt(pow(DistA,2) - pow(x,2) - pow(y,2))

#triPt is an array with ECEF x,y,z of trilateration point
triPt = P1 + x*ex + y*ey + z*ez

#convert back to lat/long from ECEF
#convert to degrees
lat = math.degrees(math.asin(triPt[2] / earthR))
lon = math.degrees(math.atan2(triPt[1],triPt[0]))

print lat, lon
wwnick
fuente
1
Iba a armar una respuesta similar, ¡pero ahora no hay necesidad! Obtiene mi voto a favor.
Wrass
numpy al rescate! Se compila cuando 'triPt' se reemplaza por 'triLatPt', pero de lo contrario devuelve 37.4191023738 -121.960579208. Buen trabajo
WolfOdrade
¡Buen trabajo! Si reemplazo el sistema de coordenadas geográficas por un sistema de coordenadas local [cartesiano], ¿seguirá funcionando?
zengr
para aquellos en el dominio c ++ ... hackearon uno realmente rápido pastebin.com/9Dur6RAP
raaj
2
Gracias @wwnick! Lo porté a JavaScript (destinado para Node pero se puede convertir fácilmente para que funcione en el navegador). gist.github.com/dav-/bb7103008cdf9359887f
DC_
6

No estoy seguro de si soy ingenuo, pero, ¿si amortigua cada punto por tamaño y luego se cruza con los tres círculos que le darían la ubicación correcta?

Puede calcular la intersección utilizando API espaciales. Ejemplos:

  • GeoScript
  • Conjunto de topología de Java
  • Conjunto de topología NET
  • GEOS
George Silva
fuente
1
Exactamente, está interesado en las fórmulas para obtener ese punto de intersección.
Vinko Vrsalovic
Usando una API espacial puede hacerlo sin usar matemática pura.
George Silva
1
@ George ¿Puedes dar un ejemplo de tal API?
nohat
Publicación editada para reflejar la solicitud de nohat.
George Silva
+1, buen pensamiento lateral, ¡aunque quizás no sea el más eficiente desde el punto de vista informático!
fmark
2

Las siguientes notas utilizan geometría planarítmica (es decir, tendría que proyectar sus coordenadas en un sistema de coordenadas local apropiado).

Mi razonamiento, con un ejemplo trabajado en Python, sigue:

Tome 2 de los puntos de datos (llámelos ay b). Llama a nuestro punto objetivo x. Ya sabemos las distancias axy bx. Podemos calcular la distancia abusando el teorema de Pitágoras.

>>> import math
>>> a = (1, 4)
>>> b = (3, 6)
>>> dist_ax = 3
>>> dist_bx = 5.385
# Pythagoras's theorem
>>> dist_ab = math.sqrt(abs(a[0]-b[0])**2 + abs(a[1]-b[1])**2)
>>> dist_ab
2.8284271247461903

Ahora, puedes calcular los ángulos de estas líneas:

>>> angle_abx = math.acos((dist_bx * dist_bx + dist_ab * dist_ab - dist_ax * dist_ax)/(2 * dist_bx * dist_ab))
>>> math.degrees(angle_abx)
23.202973815040256
>>> angle_bax = math.acos((dist_ax * dist_ax + dist_ab * dist_ab - dist_bx * dist_bx)/(2 * dist_ax * dist_ab))
>>> math.degrees(angle_bax)
134.9915256259537
>>> angle_axb = math.acos((dist_ax * dist_ax + dist_bx * dist_bx - dist_ab * dist_ab)/(2 * dist_ax * dist_bx))
>>> math.degrees(angle_axb)
21.805500559006095

Desafortunadamente, tengo poco tiempo para completar la respuesta, sin embargo, ahora que conoce los ángulos, puede calcular dos ubicaciones posibles para x. Luego, utilizando el tercer punto c, puede calcular qué ubicación es la correcta.

fmark
fuente
2

Esto podría funcionar. Rápidamente de nuevo en python, podría poner esto en el cuerpo de una función xN, yN = coordenadas de puntos, r1 y r2 = valores de radio

dX = x2 - x1
dY = y2 - y1

centroidDistance = math.sqrt(math.pow(e,2) + math.pow(dY,2)) #distance from centroids
distancePL = (math.pow(centroidDistance,2) + (math.pow(r1,2) - math.pow(r2,2))) / (2 * centroidDistance) #distance from point to a line splitting the two centroids

rx1 = x1 + (dX *k)/centroidDistance + (dY/centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))
ry1 = y1 + (dY*k)/centroidDistance - (dX /centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))

rx2 = x1 + (dX *k)/centroidDistance - (dY/centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))
ry2 = y1 + (dY*k)/centroidDistance + (dX /centroidDistance) * math.sqrt(math.pow(r1,2) - math.pow(distancePL,2))

Los valores de rx & ry son los valores de retorno (deben estar en una matriz) de los dos puntos de intersección en un círculo, si eso ayuda a aclarar las cosas.

Haga esto para los primeros 2 círculos, luego nuevamente para el primero y el último. Si alguno de los resultados de la primera iteración se compara con los resultados de la segunda (dentro de cierta tolerancia, probablemente), entonces tiene el punto de intersección. No es una gran solución, especialmente cuando comienzas a agregar más de puntos en el proceso, pero es lo más simple que puedo ver sin resolver un sistema de ecuaciones.

WolfOdrade
fuente
¿Qué son 'e' y 'k' en su código?
ReinierDG
No recuerdo :-) la respuesta de wwnick es más similar a algo que te gustaría implementar si solo tienes tres círculos.
WolfOdrade
1

Puede usar API espacial de postgis (St_Intersection, St_buffer funciones). Como notó fmark, también debe recordar que Postgis usa algoritmos planos, pero para áreas pequeñas, el uso de la proyección equidistante no introduce muchos errores.

stachu
fuente
PostGIS puede hacer cálculos esferoidales utilizando el GEOGRAPHYtipo en lugar del GEOMETRYtipo.
fmark
1

Hazlo en lenguaje PHP:

// suponiendo elevación = 0
$ tierraR = 6371; // en km (= 3959 en millas)

$ LatA = 37.418436;
$ LonA = -121.963477;
$ DistA = 0.265710701754;

$ LatB = 37.417243;
$ LonB = -121.961889;
$ DistB = 0.234592423446;

$ LatC = 37.418692;
$ LonC = -121.960194;
$ DistC = 0.0548954278262;

/ *
#uso de esfera authalic
# si usa un elipsoide, este paso es ligeramente diferente
#Convertir geodésico Lat / Long a ECEF xyz
# 1. Convertir Lat / Long a radianes
# 2. Convertir Lat / Long (radianes) a ECEF
* /
$ xA = $ earthR * (cos (deg2rad ($ LatA)) * cos (deg2rad ($ LonA)));
$ yA = $ earthR * (cos (deg2rad ($ LatA)) * sin (deg2rad ($ LonA)));
$ zA = $ earthR * (sin (deg2rad ($ LatA)));

$ xB = $ earthR * (cos (deg2rad ($ LatB)) * cos (deg2rad ($ LonB)));
$ yB = $ earthR * (cos (deg2rad ($ LatB)) * sin (deg2rad ($ LonB)));
$ zB = $ earthR * (sin (deg2rad ($ LatB)));

$ xC = $ earthR * (cos (deg2rad ($ LatC)) * cos (deg2rad ($ LonC)));
$ yC = $ earthR * (cos (deg2rad ($ LatC)) * sin (deg2rad ($ LonC)));
$ zC = $ earthR * (sin (deg2rad ($ LatC)));

/ *
INSTALAR:
sudo pear install Math_Vector-0.7.0
sudo pear install Math_Matrix-0.8.7
* /
// Incluir PEAR :: Math_Matrix
// /usr/share/php/Math/Matrix.php
// include_path = ".: / usr / local / php / pear /"
require_once 'Math / Matrix.php';
require_once 'Math / Vector.php';
require_once 'Math / Vector3.php';


$ P1vector = new Math_Vector3 (array ($ xA, $ yA, $ zA));
$ P2vector = new Math_Vector3 (array ($ xB, $ yB, $ zB));
$ P3vector = new Math_Vector3 (array ($ xC, $ yC, $ zC));

# de wikipedia: http://en.wikipedia.org/wiki/Trilateration
#transformar para obtener el círculo 1 en origen
#transformar para obtener el círculo 2 en el eje x

// CALC EX
$ P2minusP1 = Math_VectorOp :: substract ($ P2vector, $ P1vector);
$ l = nuevo Math_Vector ($ P2minusP1);
$ P2minusP1_length = $ l-> length ();
$ norma = nuevo Math_Vector3 (matriz ($ P2minusP1_length, $ P2minusP1_length, $ P2minusP1_length));
$ d = $ norma; // guardar calc D
$ ex = Math_VectorOp :: divide ($ P2minusP1, $ norma);
// echo "ex:". $ ex-> toString (). "\ n";
$ ex_x = floatval ($ ex -> _ tuple-> getData () [0]);
$ ex_y = floatval ($ ex -> _ tuple-> getData () [1]);
$ ex_z = floatval ($ ex -> _ tuple-> getData () [2]);
$ ex = new Math_Vector3 (matriz ($ ex_x, $ ex_y, $ ex_z));

// CALC i
$ P3minusP1 = Math_VectorOp :: substract ($ P3vector, $ P1vector);
$ P3minusP1_x = floatval ($ P3minusP1 -> _ tuple-> getData () [0]);
$ P3minusP1_y = floatval ($ P3minusP1 -> _ tuple-> getData () [1]);
$ P3minusP1_z = floatval ($ P3minusP1 -> _ tuple-> getData () [2]);
$ P3minusP1 = new Math_Vector3 (array ($ P3minusP1_x, $ P3minusP1_y, $ P3minusP1_z));
$ i = Math_VectorOp :: dotProduct ($ ex, $ P3minusP1);
// echo "i = $ i \ n";

// CALC EY
$ iex = Math_VectorOp :: scale ($ i, $ ex);
// echo "iex =". $ iex-> toString (). "\ n";
$ P3P1iex = Math_VectorOp :: substract ($ P3minusP1, $ iex);
// echo "P3P1iex =". $ P3P1iex-> toString (). "\ n";
$ l = nuevo Math_Vector ($ P3P1iex);
$ P3P1iex_length = $ l-> length ();
$ norma = nuevo Math_Vector3 (matriz ($ P3P1iex_length, $ P3P1iex_length, $ P3P1iex_length));
// echo "norma:". $ norma-> toString (). "\ n";
$ ey = Math_VectorOp :: divide ($ P3P1iex, $ norma);
// echo "ey =". $ ey-> toString (). "\ n";
$ ey_x = floatval ($ ey -> _ tuple-> getData () [0]);
$ ey_y = floatval ($ ey -> _ tuple-> getData () [1]);
$ ey_z = floatval ($ ey -> _ tuple-> getData () [2]);
$ ey = new Math_Vector3 (array ($ ey_x, $ ey_y, $ ey_z));

// CALC EZ
$ ez = Math_VectorOp :: crossProduct ($ ex, $ ey);
// echo "ez =". $ ez-> toString (). "\ n";

// CALC D
// hazlo antes
$ d = floatval ($ d -> _ tupla-> getData () [0]);
// echo "d = $ d \ n";

// CALC J
$ j = Math_VectorOp :: dotProduct ($ ey, $ P3minusP1);
// echo "j = $ j \ n";

# de wikipedia
#plug and chug usando los valores anteriores
$ x = (pow ($ DistA, 2) - pow ($ DistB, 2) + pow ($ d, 2)) / (2 * $ d);
$ y = ((pow ($ DistA, 2) - pow ($ DistC, 2) + pow ($ i, 2) + pow ($ j, 2)) / (2 * $ j)) - (($ i / $ j) * $ x);

# solo se muestra un caso aquí
$ z = sqrt (pow ($ DistA, 2) - pow ($ x, 2) - pow ($ y, 2));

// echo "x = $ x - y = $ y - z = $ z \ n";

#triPt es una matriz con ECEF x, y, z del punto de trilateración
$ xex = Math_VectorOp :: scale ($ x, $ ex);
$ yey = Math_VectorOp :: scale ($ y, $ ey);
$ zez = Math_VectorOp :: scale ($ z, $ ez);

// CALC $ triPt = $ P1vector + $ xex + $ yey + $ zez;
$ triPt = Math_VectorOp :: add ($ P1vector, $ xex);
$ triPt = Math_VectorOp :: add ($ triPt, $ yey);
$ triPt = Math_VectorOp :: add ($ triPt, $ zez);
// echo "triPt =". $ triPt-> toString (). "\ n";
$ triPt_x = floatval ($ triPt -> _ tuple-> getData () [0]);
$ triPt_y = floatval ($ triPt -> _ tuple-> getData () [1]);
$ triPt_z = floatval ($ triPt -> _ tuple-> getData () [2]);


#convertir de nuevo a lat / long desde ECEF
#convertir a grados
$ lat = rad2deg (asin ($ triPt_z / $ earthR));
$ lon = rad2deg (atan2 ($ triPt_y, $ triPt_x));

echo $ lat. ','. $ lon;
fquinto
fuente