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

34

Quiero descubrir una ubicación de destino 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 de destino. ¿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, desde enlace Entiendo lo general principio de que con tres círculos como estos se obtiene exactamente un punto de superposición. Lo que me molesta es la matemática necesaria para calcular ese punto con esta entrada.

    
pregunta nohat 22.07.2010 - 22:19

6 respuestas

34

Después de observar a Wikipedia y la misma pregunta / respuesta en StackOverflow , pensé que lo intentaría y trataría de rellenar el brechas.

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

Lasmatemáticasenlapáginadewikipedianosontanmalas,solonecesitaconvertirsuscoordenadasgeodésicasalECEFcartesiano,quesepuedeencontrar here . los términos a / x + h pueden ser reemplazados por el radio de la esfera auténtica, si no estás usando un elipsoide.

Probablemente, lo más fácil es simplemente 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
    
respondido por el wwnick 28.07.2010 - 13:45
6

No estoy seguro de ser ingenuo, pero, ¿si amortiguas cada punto por tamaño y luego intersectas los tres círculos que te darían la ubicación correcta?

Puedes calcular la intersección usando API espaciales. Ejemplos:

  • GeoScript
  • Suite de topología de Java
  • NET Topology Suite
  • GEOS
respondido por el George Silva 22.07.2010 - 22:30
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:

Toma 2 de los puntos de datos (llámalos a y b ). Llame a nuestro punto de destino x . Ya conocemos las distancias ax y bx . Podemos calcular la distancia ab usando 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

Lamentablemente no tengo 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 correcta.

    
respondido por el fmark 26.07.2010 - 02:35
2

Esto podría funcionar. Rápidamente de nuevo en python, puedes poner esto en el cuerpo de una función xN, yN = coordenadas de puntos, r1 & 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))

rx & Los valores de 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.

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

    
respondido por el WolfOdrade 28.07.2010 - 03:57
1

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

    
respondido por el stachu 23.07.2010 - 10:09
1

Hazlo en lenguaje PHP:

//assuming elevation = 0
$earthR = 6371; // in km ( = 3959 in miles)

$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 *(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)));

/*
INSTALL:
sudo pear install Math_Vector-0.7.0
sudo pear install Math_Matrix-0.8.7
*/
// Include 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));

#from wikipedia: http://en.wikipedia.org/wiki/Trilateration
#transform to get circle 1 at origin
#transform to get circle 2 on x axis

// CALC EX
$P2minusP1 = Math_VectorOp::substract($P2vector, $P1vector);
$l = new Math_Vector($P2minusP1);
$P2minusP1_length = $l->length();
$norm = new Math_Vector3(array($P2minusP1_length,$P2minusP1_length,$P2minusP1_length));
$d = $norm; //save calc D
$ex = Math_VectorOp::divide($P2minusP1, $norm);
//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(array($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 = new Math_Vector($P3P1iex);
$P3P1iex_length = $l->length();
$norm = new Math_Vector3(array($P3P1iex_length,$P3P1iex_length,$P3P1iex_length));
//echo "norm: ".$norm->toString()."\n";
$ey = Math_VectorOp::divide($P3P1iex, $norm);
//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
// do it before
$d = floatval( $d->_tuple->getData()[0] ) ;
//echo "d = $d\n";

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

#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 = sqrt( pow($DistA,2) - pow($x,2) - pow($y,2) );

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

#triPt is an array with ECEF x,y,z of trilateration point
$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] ) ;


#convert back to lat/long from ECEF
#convert to degrees
$lat = rad2deg(asin($triPt_z / $earthR));
$lon = rad2deg(atan2($triPt_y,$triPt_x));

echo $lat .','. $lon;
    
respondido por el fquinto 25.08.2014 - 14:06

Lea otras preguntas en las etiquetas