Trilateracja przy użyciu 3 punktów szerokości / długości geograficznej i 3 odległości?

34

Chcę znaleźć nieznaną lokalizację docelową (współrzędne szerokości i długości geograficznej). Istnieją 3 znane punkty (pary współrzędnych szerokości i długości geograficznej) i dla każdego punktu odległość w kilometrach od miejsca docelowego. Jak mogę obliczyć współrzędne docelowej lokalizacji?

Powiedzmy na przykład, że mam następujące punkty danych

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

To, co chciałbym, to matematyka dla funkcji, która przyjmuje to jako dane wejściowe i zwraca 37,417959, -121,961954 jako dane wyjściowe.

Rozumiem, jak obliczyć odległość między dwoma punktami, na stronie http://www.movable-type.co.uk/scripts/latlong.html Rozumiem ogólną zasadę, że przy takich trzech okręgach otrzymujesz dokładnie jeden punkt nakładania się. Mętnie rozumiem matematykę potrzebną do obliczenia tego punktu za pomocą tych danych wejściowych.

nohat
źródło
Oto strona, która przeprowadzi cię przez matematykę znajdowania środka trzech współrzędnych. Być może może to w jakiś sposób pomóc. < mathforum.org/library/drmath/view/68373.html >
Jon Bringhurst
1
Czy to musi być na kuli / sferoidzie, czy algorytm planarny jest w porządku?
fmark
1
Nie mogę dać ci odpowiedzi, ale myślę, że mogę skierować cię w dobrym kierunku. Trzy współrzędne = trzy punkty środkowe. Trzy odległości = trzy koła. Dwa koła, które się przecinają, mogą mieć możliwość braku jednego / jednego / dwóch rozwiązań. Trzy Kręgi mogą mieć brak / jeden / lub obszar jako rozwiązanie. Uzyskaj formułę okręgu dla trzech kół i rozwiąż ją za pomocą układów równań / algebry.
CrazyEnigma
W rzeczywistości nie potrzebujesz nawet systemów do rozwiązania tego. Istnieje jedna lub dwie możliwości, ale ponieważ masz wartość odległości, możesz oddzielić poprawną odpowiedź.
George Silva,
1
+1 To dobre pytanie. Na początku pomyślałem, że w Google można łatwo znaleźć rozwiązanie, ale najwyraźniej nie. Być może problem można by bardziej ogólnie powiedzieć: biorąc pod uwagę N punktów, przy czym każdy punkt ma nie tylko odległość, ale także margines błędu, znajdź elipsę zaufania.
Kirk Kuykendall

Odpowiedzi:

34

Po tym, jak rozejrzałem się po Wikipedii i tym samym pytaniu / odpowiedzi na StackOverflow , pomyślałem, że zrobię sobie dźgnięcie i spróbuję uzupełnić luki.

Po pierwsze, Nie jestem pewien, skąd masz wynik, ale wydaje się, że jest nieprawidłowy. Wykreśliłem punkty w ArcMap, zbuforowałem je do określonych odległości, pobiegłem przeciąć na buforach, a następnie uchwyciłem wierzchołek przecięcia, aby uzyskać rozwiązania. Proponowany wynik jest zielony. Obliczyłem wartość w oknie objaśnień, czyli około 3 metry tego, co ArcMap dał dla rozwiązania uzyskanego z przecięcia.

alternatywny tekst

Matematyka na stronie wikipedii nie jest taka zła, wystarczy zakryć swoje współrzędne geodezyjne kartezjańskim ECEF, który można znaleźć tutaj . warunki a / x + h można zastąpić promieniem kuli authalicznej, jeśli nie używasz elipsoidy.

Prawdopodobnie najłatwiej jest po prostu podać ci dobrze (?) Udokumentowany kod, więc tutaj jest w Pythonie

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
źródło
1
Zamierzałem ułożyć podobną odpowiedź, ale teraz nie ma takiej potrzeby! Dostaje moje poparcie.
Wrass
numpy na ratunek! Kompiluje się, gdy „triPt” jest zastępowany przez „triLatPt”, ale w przeciwnym razie zwraca 37,4191023738 -121,960579208. Dobra robota
WolfOdrade
Dobra robota! Jeśli zastąpię układ współrzędnych geograficznych lokalnym układem współrzędnych [kartezjańskich], czy to nadal będzie działać?
zengr
dla osób w domenie c ++..hacked naprawdę szybki pastebin.com/9Dur6RAP
raaj
2
Dzięki @wwnick! Przeniesiłem to na JavaScript (przeznaczony dla Node, ale można go łatwo przekonwertować do pracy w przeglądarce). gist.github.com/dav-/bb7103008cdf9359887f
DC_
6

Nie jestem pewien, czy jestem naiwny, ale czy buforujesz każdy punkt według rozmiaru, a następnie przecinasz wszystkie trzy koła, które zapewniłyby ci poprawną lokalizację?

Możesz obliczyć skrzyżowanie za pomocą przestrzennych interfejsów API. Przykłady:

  • GeoScript
  • Java Topology Suite
  • NET Topology Suite
  • GEOS
George Silva
źródło
1
Dokładnie, on jest zainteresowany formułami, aby uzyskać ten punkt przecięcia.
Vinko Vrsalovic
Korzystając z przestrzennego interfejsu API, możesz to zrobić bez użycia czystej matematyki.
George Silva,
1
@George czy możesz podać przykład takiego interfejsu API?
no
Edytowany post odzwierciedlający prośbę nohat.
George Silva,
+1, dobre myślenie boczne, nawet jeśli nie najbardziej wydajne obliczeniowo!
fmark
2

Poniższe uwagi wykorzystują geometrię planaryczną (tj. Musiałbyś rzutować swoje współrzędne na odpowiedni lokalny układ współrzędnych).

Moje rozumowanie z działającym przykładem w Pythonie wygląda następująco:

Weź 2 punkty danych (zadzwoń do nich ai b). Zadzwoń do naszego punktu docelowego x. Znamy już odległości axi bx. Odległość możemy obliczyć za abpomocą twierdzenia Pitagorasa.

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

Teraz możesz obliczyć kąty tych linii:

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

Niestety brakuje mi czasu na udzielenie odpowiedzi, jednak teraz znasz kąty, możesz obliczyć dwie możliwe lokalizacje x. Następnie, korzystając z trzeciego punktu c, możesz obliczyć, która lokalizacja jest poprawna.

fmark
źródło
2

To może zadziałać. Szybko ponownie w pythonie, możesz umieścić to w ciele funkcji xN, yN = współrzędne punktów, r1 i r2 = wartości promienia

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

Wartości rx i ry to wartości zwracane (powinny znajdować się w tablicy) dwóch punktów przecięcia na okręgu, jeśli to pomaga wyjaśnić.

Zrób to dla pierwszych 2 kręgów, a następnie dla pierwszego i ostatniego. Jeśli którykolwiek z wyników z pierwszej iteracji porówna się z wynikami z drugiego (prawdopodobnie w ramach pewnej tolerancji), to masz punkt przecięcia. Nie jest to świetne rozwiązanie, zwłaszcza gdy zaczynasz dodawać do procesu więcej niż punkty, ale jest to najprostsze, jakie widzę bez rozwiązywania układu równań.

WolfOdrade
źródło
Co to jest „e” i „k” w twoim kodzie?
ReinierDG
Nie pamiętam :-) odpowiedź wwnick jest bardziej podobna do czegoś, co chciałbyś wdrożyć, gdybyś miał tylko trzy kręgi.
WolfOdrade
1

Możesz używać przestrzennego API z postgis (funkcje St_Intersection, St_buffer). Jak zauważył fmark, należy również pamiętać, że Postgis stosuje algorytmy planarne, ale w przypadku małych obszarów stosowanie równomiernego przesuwania nie wprowadza dużego błędu.

stachu
źródło
PostGIS może wykonywać obliczenia sferoidalne przy użyciu GEOGRAPHYtypu, a nie GEOMETRYtypu.
fmark
1

Zrób to w języku PHP:

// przy założeniu rzędnej = 0
$ earthR = 6371; // w km (= 3959 w milach)

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

/ *
#z wykorzystaniem kuli autologicznej
#if przy użyciu elipsoidy ten krok jest nieco inny
# Przelicz geodezyjny Lat / Long na ECEF xyz
# 1. Konwertuj Lat / Long na radiany
# 2. Konwertuj Lat / Long (radiany) na 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)));

/ *
ZAINSTALOWAĆ:
sudo gruszka zainstalować Math_Vector-0.7.0
sudo gruszka zainstalować Math_Matrix-0.8.7
* /
// Dołącz PEAR :: Math_Matrix
// /usr/share/php/Math/Matrix.php
// include_path = ".: / usr / local / php / pear /"
requ_once 'Math / Matrix.php';
requ_once 'Math / Vector.php';
requ_once 'Math / Vector3.php';


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

# z wikipedii: http://en.wikipedia.org/wiki/Trilateration
#transformuj, aby uzyskać krąg 1 w punkcie początkowym
#transformuj, aby uzyskać okrąg 2 na osi x

// CALC EX
$ P2minusP1 = Math_VectorOp :: subtract ($ P2vector, $ P1vector);
$ l = nowy Math_Vector ($ P2minusP1);
$ P2minusP1_length = $ l-> length ();
$ norm = nowy Math_Vector3 (tablica ($ P2minusP1_length, $ P2minusP1_length, $ P2minusP1_length));
$ d = $ norma; // zapisz oblicz 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 = nowy Math_Vector3 (tablica ($ ex_x, $ ex_y, $ ex_z));

// CALC i
$ P3minusP1 = Math_VectorOp :: subtract ($ P3vector, $ P1vector);
$ P3minusP1_x = zmiennoprzecinkowy ($ P3minusP1 -> _ tuple-> getData () [0]);
$ P3minusP1_y = floatval ($ P3minusP1 -> _ tuple-> getData () [1]);
$ P3minusP1_z = floatval ($ P3minusP1 -> _ tuple-> getData () [2]);
$ P3minusP1 = nowy Math_Vector3 (tablica ($ 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 :: subtract ($ P3minusP1, $ iex);
// echo "P3P1iex =". $ P3P1iex-> toString (). "\ n";
$ l = nowy Math_Vector ($ P3P1iex);
$ P3P1iex_length = $ l-> length ();
$ norm = nowy Math_Vector3 (tablica ($ 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 = nowy Math_Vector3 (tablica ($ ey_x, $ ey_y, $ ey_z));

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

// CALC D.
// zrób to wcześniej
$ d = floatval ($ d -> _ tuple-> getData () [0]);
// echo "d = $ d \ n";

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

# z wikipedii
#plug and chug przy użyciu powyższych wartości
$ 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);

# pokazano tylko jeden przypadek
$ z = sqrt (pow ($ DistA, 2) - pow ($ x, 2) - pow ($ y, 2));

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

#triPt to tablica z ECEF x, y, z punktu trilateracji
$ 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 = zmiennoprzecinkowy ($ triPt -> _ tuple-> getData () [0]);
$ triPt_y = zmiennoprzecinkowe ($ triPt -> _ tuple-> getData () [1]);
$ triPt_z = floatval ($ triPt -> _ tuple-> getData () [2]);


#convert powrót do lat / long z ECEF
#konwertuj na stopnie
$ lat = rad2deg (asin ($ triPt_z / $ earthR));
$ lon = rad2deg (atan2 ($ triPt_y, $ triPt_x));

echo $ lat. ','. $ lon;
fintinto
źródło