Trilateration mit 3 Breiten- / Längenpunkten und 3 Entfernungen?

34

Ich möchte einen unbekannten Zielort (Längen- und Breitengradkoordinaten) herausfinden. Es gibt 3 bekannte Punkte (Koordinatenpaare von Breite und Länge) und für jeden Punkt eine Entfernung in Kilometern zum Zielort. Wie kann ich die Koordinaten des Zielorts berechnen?

Angenommen, ich habe die folgenden Datenpunkte

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

Was ich möchte, ist, was die Mathematik für eine Funktion ist, die das als Eingabe annimmt und 37.417959, -121.961954 als Ausgabe zurückgibt.

Ich verstehe, wie man den Abstand zwischen zwei Punkten berechnet, von http://www.movable-type.co.uk/scripts/latlong.html Ich verstehe das allgemeine Prinzip, dass man mit drei Kreisen wie diesen genau einen Überlappungspunkt erhält. Worauf ich mich einlasse, ist die Mathematik, die benötigt wird, um diesen Punkt mit dieser Eingabe zu berechnen.

nohat
quelle
Hier ist eine Seite, die Sie durch die Mathematik des Findens des Zentrums von drei Koordinaten führt. Vielleicht könnte es irgendwie helfen. < mathforum.org/library/drmath/view/68373.html >
Jon Bringhurst
1
Muss dies auf der Kugel sein, oder ist ein planarer Algorithmus in Ordnung?
Markieren Sie den
1
Ich kann Ihnen keine Antwort geben, aber ich glaube, ich kann Sie in die richtige Richtung weisen. Drei Koordinaten = drei Mittelpunkte. Drei Entfernungen = drei Kreise. Zwei Kreise, die sich kreuzen, können die Möglichkeit haben, keine / eine / zwei Lösungen zu finden. Drei Kreise können keine / eine / oder eine Fläche als Lösung haben. Sie erhalten die Kreisformel für die drei Kreise und lösen sie mit Gleichungssystemen / Algebra.
CrazyEnigma
Eigentlich brauchen Sie nicht einmal Systeme, um dieses Problem zu lösen. Es gibt ein oder zwei Möglichkeiten, aber da Sie einen Entfernungswert haben, können Sie die richtige Antwort trennen.
George Silva
1
+1 Das ist eine gute Frage. Zuerst dachte ich, dass eine Lösung mit Google leicht gefunden werden könnte, aber anscheinend nicht. Vielleicht könnte das Problem allgemeiner ausgedrückt werden: Wenn N Punkte gegeben sind, wobei jeder Punkt nicht nur eine Entfernung, sondern auch eine Fehlertoleranz hat, finde die Vertrauensellipse.
Kirk Kuykendall

Antworten:

34

Nach einigem Hin und Her bei Wikipedia und der gleichen Frage / Antwort bei StackOverflow dachte ich, ich würde es versuchen, die Lücken zu füllen.

Zunächst einmal: Ich bin mir nicht sicher, woher die Ausgabe stammt, aber sie scheint falsch zu sein. Ich habe die Punkte in ArcMap geplottet, sie auf die angegebenen Entfernungen gepuffert, Schnittpunkte auf den Puffern erstellt und dann den Schnittpunkt erfasst, um die Lösungen zu erhalten. Ihre vorgeschlagene Ausgabe ist der Punkt in Grün. Ich habe den Wert in der Callout-Box berechnet, der ungefähr 3 Meter des Werts entspricht, den ArcMap für die vom Schnittpunkt abgeleitete Lösung ergab.

Alt-Text

Die Mathematik auf der Wikipedia-Seite ist nicht allzu schlecht. Sie müssen nur Ihre geodätischen Koordinaten in den kartesischen ECEF umwandeln, den Sie hier finden . Die a / x + h-Terme können durch den authalen Kugelradius ersetzt werden, wenn Sie kein Ellipsoid verwenden.

Am einfachsten ist es wahrscheinlich, wenn Sie einen gut (?) Dokumentierten Code angeben. Hier also in 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
quelle
1
Ich wollte eine ähnliche Antwort zusammenstellen, aber jetzt gibt es keine Notwendigkeit! Ruft meine Zustimmung.
Wrass
zur Rettung betäubt! Es wird kompiliert, wenn 'triPt' durch 'triLatPt' ersetzt wird, gibt aber ansonsten 37.4191023738 -121.960579208 zurück. Gute Arbeit
WolfOdrade
Gute Arbeit! Funktioniert das immer noch, wenn ich das geografische Koordinatensystem durch ein lokales [kartesisches] Koordinatensystem ersetze?
Zengr
für die in der c ++ domain..hacked zusammen eine ganz schnell ein pastebin.com/9Dur6RAP
Raaj
2
Danke @wwnick! Ich habe dies nach JavaScript portiert (für Node gedacht, kann aber problemlos für die Arbeit im Browser konvertiert werden). gist.github.com/dav-/bb7103008cdf9359887f
DC_
6

Ich bin mir nicht sicher, ob ich naiv bin, aber wenn Sie jeden Punkt nach Größe puffern und dann alle drei Kreise schneiden, die Ihnen die richtige Position bringen würden?

Sie können die Schnittmenge mit räumlichen APIs berechnen. Beispiele:

  • GeoScript
  • Java Topology Suite
  • NET Topology Suite
  • GEOS
George Silva
quelle
1
Genau, er interessiert sich für die Formeln, um diesen Schnittpunkt zu erhalten.
Vinko Vrsalovic
Mit einer räumlichen API können Sie dies ohne reine Mathematik tun.
George Silva
1
@ George können Sie ein Beispiel für eine solche API geben?
Mitternacht,
Beitrag bearbeitet, um die Anfrage von nohat widerzuspiegeln.
George Silva
+1, gutes Querdenken, auch wenn es vielleicht nicht das recheneffizienteste ist!
Markieren Sie den
2

Die folgenden Hinweise verwenden planarithmische Geometrie (dh Sie müssten Ihre Koordinaten in ein geeignetes lokales Koordinatensystem projizieren).

Meine Argumentation mit einem Beispiel in Python lautet wie folgt:

Nimm 2 der Datenpunkte (nenne sie aund b). Rufen Sie unseren Zielpunkt an x. Wir kennen bereits die Entfernungen axund bx. Wir können die Entfernung abmit dem Satz von Pythagoras berechnen .

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

Nun können Sie die Winkel dieser Linien berechnen:

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

Leider fehlt mir die Zeit, um die Antwort für Sie zu vervollständigen. Nachdem Sie jedoch die Winkel kennen, können Sie zwei mögliche Positionen für berechnen x. Mit dem dritten Punkt c können Sie dann berechnen, welcher Ort korrekt ist.

fmark
quelle
2

Das könnte funktionieren. In Python können Sie dies schnell wieder in den Rumpf einer Funktion einfügen: xN, yN = Koordinaten von Punkten, r1 & r2 = Radiuswerte

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 & ry-Werte sind die Rückgabewerte (sollten in einem Array liegen) der beiden Schnittpunkte auf einem Kreis, wenn dies zur Klärung beiträgt.

Tun Sie dies für die ersten 2 Kreise, dann erneut für den ersten und den letzten. Wenn eines der Ergebnisse der ersten Iteration mit den Ergebnissen der zweiten Iteration verglichen wird (wahrscheinlich innerhalb einer gewissen Toleranz), haben Sie den Schnittpunkt. Es ist keine großartige Lösung, besonders wenn Sie anfangen, mehr als Punkte in den Prozess einzufügen, aber es ist die einfachste, die ich sehen kann, ohne ein Gleichungssystem zu lösen.

WolfOdrade
quelle
Was sind "e" und "k" in Ihrem Code?
ReinierDG
Ich erinnere mich nicht :-) Die Antwort von wwnick ist eher so, wie man es mit etwas umsetzen möchte, wenn man nur drei Kreise hat.
WolfOdrade
1

Sie können die räumliche API von postgis verwenden (St_Intersection, St_buffer-Funktionen). Wie fmark bemerkt hat, müssen Sie sich auch daran erinnern, dass Postgis planare Algorithmen verwendet, aber für kleine Bereiche führt die Verwendung von äquidistanten Projektionen nicht zu großen Fehlern.

Stachu
quelle
PostGIS kann sphärische Berechnungen anhand des GEOGRAPHYTyps und nicht anhand des GEOMETRYTyps durchführen.
Markieren Sie den
1

Mach es in PHP Sprache:

// Unter der Annahme von Höhe = 0
$ earthR = 6371; // in km (= 3959 in Meilen)

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

/ *
# unter Verwendung der authalen Kugel
#wenn ein Ellipsoid verwendet wird, unterscheidet sich dieser Schritt geringfügig
#Geodätisches Lat / Long in ECEF xyz konvertieren
# 1. Lat / Long in Radiant umrechnen
# 2. Lat / Long (Bogenmaß) in ECEF umrechnen
* /
$ 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)));

/ *
INSTALLIEREN:
sudo pear installiert Math_Vector-0.7.0
sudo pear install Math_Matrix-0.8.7
* /
// PEAR :: Math_Matrix einbinden
// /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));

#aus Wikipedia: http://en.wikipedia.org/wiki/Trilateration
#transformiere, um Kreis 1 am Ursprung zu erhalten
#transformiere, um Kreis 2 auf der x-Achse zu erhalten

// CALC EX
$ P2minusP1 = Math_VectorOp :: substract ($ P2vector, $ P1vector);
$ l = neuer Math_Vector ($ P2minusP1);
$ P2minusP1_length = $ l-> length ();
$ norm = new Math_Vector3 (Array ($ P2minusP1_length, $ P2minusP1_length, $ P2minusP1_length));
$ d = $ norm; // Berechnung speichern 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 = neuer 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
// mach es vorher
$ d = floatval ($ d -> _tuple-> getData () [0]);
// echo "d = $ d \ n";

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

#aus Wikipedia
#plug und chug mit den obigen Werten
$ 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);

# Hier wird nur ein Fall angezeigt
$ z = sqrt (pow ($ DistA, 2) - pow ($ x, 2) - pow ($ y, 2));

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

#triPt ist ein Array mit ECEF x, y, z des Trilaterationspunkts
$ 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]);


# Von ECEF zurück nach Lat / Long konvertieren
#in Grad umrechnen
$ lat = rad2deg (asin ($ triPt_z / $ earthR));
$ lon = rad2deg (atan2 ($ triPt_y, $ triPt_x));

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