Trilatération utilisant 3 points de latitude / longitude et 3 distances?

34

Je souhaite trouver une cible inconnue (coordonnées de latitude et de longitude). Il y a 3 points connus (paires de coordonnées de latitude et de longitude) et pour chaque point, une distance en kilomètres de l'emplacement cible. Comment puis-je calculer les coordonnées de l'emplacement cible?

Par exemple, disons que j'ai les points de données suivants

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

Ce que j'aimerais, c'est quel est le calcul pour une fonction qui prend cela en entrée et renvoie 37.417959, -121.961954 en sortie.

Je comprends comment calculer la distance entre deux points, à partir de http://www.movable-type.co.uk/scripts/latlong.html Je comprends le principe général voulant qu’avec trois cercles comme celui-ci, vous obtenez exactement un point de chevauchement. Ce qui me manque, c'est le calcul nécessaire pour calculer ce point avec cette entrée.

nohat
la source
Voici une page qui vous explique comment calculer le centre des trois coordonnées. Peut-être que cela pourrait aider d'une certaine manière. < mathforum.org/library/drmath/view/68373.html >
Jon Bringhurst
1
Est-ce que cela doit être sur la sphère / sphéroïde, ou un algorithme planaire est-il correct?
Fmark
1
Je ne peux pas vous donner la réponse, mais je pense pouvoir vous orienter dans la bonne direction. Trois coordonnées = trois points centraux. Trois distances = trois cercles. Deux cercles qui se croisent peuvent avoir une possibilité de solution nulle / une / deux. Trois cercles peuvent avoir aucun / un / ou un domaine comme solution. Obtenez la formule du cercle pour les trois cercles et résolvez-la avec Systèmes d’équations / Algèbre.
CrazyEnigma
En fait, vous n'avez même pas besoin de systèmes pour résoudre celui-ci. Il y a une ou deux possibilités, mais comme vous avez une valeur de distance, vous pouvez séparer la bonne réponse.
George Silva
1
+1 C'est une bonne question. Au début, j'ai pensé qu'une solution pourrait être facilement trouvée avec Google, mais apparemment pas. Le problème pourrait peut-être être posé de manière plus générale: étant donné que N points avec chaque point n’ayant pas seulement une distance mais aussi une marge d’erreur, trouve l’ellipse de confiance.
Kirk Kuykendall

Réponses:

34

Après quelques recherches sur Wikipedia et la même question / réponse sur StackOverflow , je me suis dit que j'essaierais de tenter ma chance et tenterais de combler les lacunes.

Tout d’abord, je ne sais pas où vous avez obtenu la sortie, mais cela semble être faux. J'ai tracé les points dans ArcMap, les ai mis en mémoire tampon aux distances spécifiées, puis je me suis rendu intersection sur les tampons, puis j'ai capturé le sommet de l'intersection pour obtenir les solutions. Votre sortie proposée est le point en vert. J'ai calculé la valeur dans la zone de légende, qui correspond à environ 3 mètres de ce qu'ArcMap a donné pour la solution dérivée de l'intersection.

texte alternatif

Les maths sur la page Wikipédia ne sont pas si mauvaises, il suffit de convertir vos coordonnées géodésiques en ECEF cartésien, que vous pouvez trouver ici . les termes a / x + h peuvent être remplacés par le rayon de la sphère automatique si vous n'utilisez pas d'ellipsoïde.

Le plus simple est probablement de vous donner un code bien documenté (?), Donc le voici 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
la source
1
J'allais préparer une réponse similaire, mais maintenant ce n'est plus nécessaire! Obtient mon upvote.
Wrass
Numpy à la rescousse! Il est compilé lorsque «triPt» est remplacé par «triLatPt», mais renvoie sinon 37.4191023738-121.960579208. Bon travail
WolfOdrade
Bon travail! Si je remplace le système de coordonnées géographiques par un système de coordonnées [cartésien] local, cela fonctionnera-t-il encore?
zengr
Pour ceux qui travaillent dans le domaine c ++, nous avons réuni un vrai rapide. pastebin.com/9Dur6RAP
raaj
2
Merci @wwnick! Je l'ai porté en JavaScript (destiné à Node, mais peut être facilement converti pour fonctionner dans le navigateur). gist.github.com/dav-/bb7103008cdf9359887f
DC_
6

Je ne suis pas sûr d'être naïf, mais si vous tamponnez chaque point en fonction de sa taille, puis que vous croisez les trois cercles, vous obtiendrez le bon emplacement?

Vous pouvez calculer l'intersection à l'aide d'API spatiales. Exemples:

  • GeoScript
  • Suite de topologie Java
  • Suite de topologie NET
  • GEOS
George Silva
la source
1
Exactement, il s'intéresse aux formules pour obtenir ce point d'intersection.
Vinko Vrsalovic
En utilisant une API spatiale, vous pouvez le faire sans utiliser des mathématiques pures.
George Silva
1
@ George pouvez-vous donner un exemple d'une telle API?
nohat
Poste édité pour refléter la demande de nohat.
George Silva
+1, bonne pensée latérale, même si ce n’est peut-être pas la plus efficace en calcul!
Fmark
2

Les remarques suivantes utilisent la géométrie planarithmique (c’est-à-dire que vous devez projeter vos coordonnées dans un système de coordonnées local approprié).

Mon raisonnement, avec un exemple concret en Python, est le suivant:

Prenez 2 des points de données (appelez-les aet b). Appelez notre point cible x. Nous connaissons déjà les distances axet bx. Nous pouvons calculer la distance en abutilisant le théorème de Pythagore.

>>> 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

Maintenant, vous pouvez calculer les angles de ces lignes:

>>> 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

Malheureusement, je suis à court de temps pour compléter la réponse pour vous. Cependant, maintenant que vous connaissez les angles, vous pouvez calculer deux emplacements possibles pour x. Ensuite, en utilisant le troisième point c, vous pouvez calculer quel emplacement est correct.

fmark
la source
2

Cela pourrait fonctionner. Rapidement à nouveau en python, vous pourriez le mettre dans le corps d'une fonction xN, yN = coordonnées des points, r1 & r2 = valeurs du rayon

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))

Les valeurs rx & ry sont les valeurs de retour (qui doivent figurer dans un tableau) des deux points d'intersection d'un cercle, si cela permet de clarifier les choses.

Faites cela pour les 2 premiers cercles, puis à nouveau pour le premier et le dernier. Si l'un des résultats de la première itération est comparable à celui de la seconde (dans une certaine tolérance, probablement), alors vous avez le point d'intersection. Ce n'est pas une bonne solution, surtout quand vous commencez à ajouter plus que des points dans le processus, mais c'est la plus simple que je puisse voir sans aller à la résolution d'un système d'équations.

WolfOdrade
la source
Que sont 'e' et 'k' dans votre code?
ReinierDG
Je ne me souviens pas :-) La réponse de wwnick va plus loin dans le sens de quelque chose que vous voudriez mettre en œuvre si vous n'avez que trois cercles.
WolfOdrade
1

Vous pouvez utiliser l’API spatiale de postgis (fonctions St_Intersection, St_buffer). Comme vous l'avez remarqué, vous devez également vous rappeler que Postgis utilise des algorithmes planaires, mais que pour les petites zones, l'utilisation d'une projection équi-distante n'introduit pas beaucoup d'erreur.

Stachu
la source
PostGIS peut effectuer des calculs sphéroïdaux en utilisant le GEOGRAPHYtype plutôt que le GEOMETRYtype.
Mark
1

Faites-le en langage PHP:

// en supposant que l'altitude = 0
$ earthR = 6371; // en km (= 3959 en 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;

/ *
# utilisant la sphère authalic
#si vous utilisez un ellipsoïde, cette étape est légèrement différente
#Convertir géodésique Lat / Long en ECEF xyz
# 1. Convertir Lat / Long en radians
# 2. Convertir Lat / Long (radians) en 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)));

/ *
INSTALLER:
sudo pear installer Math_Vector-0.7.0
sudo pear installer Math_Matrix-0.8.7
* /
// Inclure 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 (tableau ($ xA, $ yA, $ zA));
$ P2vector = new Math_Vector3 (tableau ($ xB, $ yB, $ zB));
$ P3vector = new Math_Vector3 (tableau ($ xC, $ yC, $ zC));

#de wikipedia: http://en.wikipedia.org/wiki/Trilateration
#transform pour obtenir le cercle 1 à l'origine
#transform pour obtenir le cercle 2 sur l'axe des x

// CALC EX
$ P2minusP1 = Math_VectorOp :: substract ($ P2vector, $ P1vector);
$ l = new Math_Vector ($ P2minusP1);
$ P2minusP1_length = $ l-> length ();
$ norm = new Math_Vector3 (tableau ($ P2minusP1_length, $ P2minusP1_length, $ P2minusP1_length));
$ d = $ norm; // sauvegarde 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 (tableau ($ 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 (tableau ($ 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 (tableau ($ 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
// le faire avant
$ d = floatval ($ d -> _ tuple-> getData () [0]);
// echo "d = $ d \ n";

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

# de wikipedia
#plug and chug en utilisant les valeurs ci-dessus
$ 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);

# un seul cas montré ici
$ z = sqrt (pow ($ DistA, 2) - pow ($ x, 2) - pow ($ y, 2));

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

#triPt est un tableau avec ECEF x, y, z du point de trilatération
$ 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 retour à lat / long de l'ECEF
#convertir en degrés
$ lat = rad2deg (asin ($ triPt_z / $ earthR));
$ lon = rad2deg (atan2 ($ triPt_y, $ triPt_x));

echo $ lat. ','. $ lon;
fquinto
la source