Трилатерация с использованием 3 точек широты / долготы и 3 расстояний?

34

Я хочу узнать неизвестное местоположение цели (координаты широты и долготы). Существует 3 известных точки (пары координат широты и долготы), и для каждой точки расстояние в километрах до местоположения цели. Как я могу рассчитать координаты целевого местоположения?

Например, скажем, у меня есть следующие данные

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

То, что я хотел бы, это то, что математика для функции, которая принимает это в качестве входных данных и возвращает 37,417959, -121,961954 в качестве выходных данных.

Я понимаю, как рассчитать расстояние между двумя точками, из http://www.movable-type.co.uk/scripts/latlong.html. Я понимаю общий принцип, согласно которому с тремя такими кружками вы получаете ровно одну точку перекрытия. Что мне не нравится, так это математика, необходимая для вычисления этой точки с помощью этого ввода.

nohat
источник
Вот страница, которая проведет вас по математике нахождения центра трех координат. Возможно, это может как-то помочь. < mathforum.org/library/drmath/view/68373.html >
Джон Брингхерст,
1
Это должно быть на сфере / сфероиде, или плоский план в порядке?
Fmark
1
Я не могу дать вам ответ, но я думаю, что могу указать вам правильное направление. Три координаты = три центральные точки. Три расстояния = три круга. Два круга, которые пересекаются, могут иметь возможность ни одного / одного / двух решений. Три Круга не могут иметь ни одного / одного / или области в качестве решения. Получите формулу окружности для трех окружностей и решите ее с помощью системы уравнений / алгебры.
CrazyEnigma
На самом деле, вам даже не нужны системы, чтобы решить эту проблему. Есть одна или две возможности, но так как у вас есть значение расстояния, вы можете отделить правильный ответ.
Джордж Сильва
1
+1 Это хороший вопрос. Сначала я подумал, что решение можно легко найти с помощью Google, но, видимо, нет. Возможно, проблема может быть сформулирована в более общем виде: учитывая N точек, в которых каждая точка имеет не только расстояние, но и предел погрешности, найдите эллипс доверия.
Кирк Куйкендалл

Ответы:

34

После некоторого осмотра Википедии и того же вопроса / ответа в StackOverflow , я подумал, что сделаю удар и попробую заполнить пробелы.

Во-первых, не уверен, где вы получили вывод, но, похоже, это неправильно. Я построил точки в ArcMap, разместил их в буфере на заданных расстояниях, пересек их с буферами и затем захватил вершину пересечения, чтобы получить решения. Ваш предложенный результат - точка зеленого цвета. Я рассчитал значение в окне выноски, которое составляет около 3 метров от того, что ArcMap дал для решения, полученного из пересечения.

альтернативный текст

Математика на странице википедии не так уж и плоха, просто нужно преобразовать ваши геодезические координаты в декартовую ECEF, которую можно найти здесь . члены a / x + h могут быть заменены автономным радиусом сферы, если вы не используете эллипсоид.

Вероятно, проще всего просто дать вам хорошо (?) Документированный код, так что вот он на 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
источник
1
Я собирался собрать аналогичный ответ, но теперь нет необходимости! Получает мой upvote.
Wrass
Numpy на помощь! Он компилируется, когда «triPt» заменяется «triLatPt», но в противном случае возвращает 37.4191023738 -121.960579208. Хорошая работа
WolfOdrade
Хорошая работа! Если я заменим географическую систему координат на локальную [декартову] систему координат, это все еще будет работать?
Зенг
для тех, кто в домене c ++ .. взломал вместе очень быстро pastebin.com/9Dur6RAP
raaj
2
Спасибо @wwnick! Я перенес это в JavaScript (предназначен для Node, но его можно легко преобразовать для работы в браузере). gist.github.com/dav-/bb7103008cdf9359887f
DC_
6

Я не уверен, что я наивен, но, если вы буферизуете каждую точку по размеру, а затем пересекаете все три круга, чтобы получить правильное местоположение?

Вы можете вычислить пересечение, используя пространственные API. Примеры:

  • GeoScript
  • Java Topology Suite
  • NET Topology Suite
  • GEOS
Джордж Сильва
источник
1
Точно, он интересуется формулами, чтобы получить эту точку пересечения.
Винко Врсалович
Используя пространственный API, вы можете сделать это без использования чистой математики.
Джордж Сильва
1
@ Джордж, можешь ли ты привести пример такого API?
нох
Отредактированный пост, чтобы отразить просьбу ничего.
Джордж Сильва
+1, хорошее боковое мышление, даже если не самое эффективное в вычислительном отношении!
Fmark
2

В следующих примечаниях используется планарифмическая геометрия (т.е. вам придется проецировать ваши координаты в соответствующую локальную систему координат).

Мои рассуждения с работающим примером на Python следующие:

Возьмите 2 точки данных (назовите их aи b). Позвоните в нашу целевую точку x. Мы уже знаем расстояния axи bx. Мы можем рассчитать расстояние, abиспользуя теорему Пифагора.

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

Теперь вы можете отработать углы этих линий:

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

К сожалению, у меня не хватает времени, чтобы завершить ответ за вас, однако теперь, когда вы знаете углы, вы можете рассчитать два возможных местоположения x. Затем, используя третью точку c, вы можете рассчитать, какое местоположение является правильным.

fmark
источник
2

Это может сработать. Снова быстро в Python, вы можете поместить это в тело функции xN, yN = координаты точек, r1 & r2 = значения радиуса

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 - это возвращаемые значения (должны быть в массиве) двух точек пересечения на окружности, если это помогает прояснить ситуацию.

Сделайте это для первых 2 кругов, затем снова для первого и последнего. Если какой-либо из результатов первой итерации сравнивается с результатами второй (возможно, в пределах некоторого допуска), то у вас есть точка пересечения. Это не очень хорошее решение, особенно когда вы начинаете добавлять больше, чем просто точки в процесс, но это самое простое, что я вижу, не решая систему уравнений.

WolfOdrade
источник
Что такое «е» и «к» в вашем коде?
ReinierDG
Я не помню :-) Ответ wwnick больше похож на то, что вы хотели бы реализовать, если у вас всего три круга.
WolfOdrade
1

Вы можете использовать пространственный API из postgis (функции St_Intersection, St_buffer). Как заметил fmark, вы также должны помнить, что Postgis использует плоские алгоритмы, но для небольших областей использование равноудаленной проекции не приводит к большой ошибке.

stachu
источник
PostGIS может выполнять сфероидальные вычисления, используя GEOGRAPHYтип, а не GEOMETRYтип.
fmark
1

Сделайте это на языке PHP:

// предполагается, что высота = 0
$ earthR = 6371; // в км (= 3959 в милях)

Лата = 37,418436;
$ LonA = -121,963477;
$ DistA = 0,265710701754;

Лат = 37,417243;
$ LonB = -121,961889;
$ DistB = 0,234592423446;

LatC = 37,418692;
$ LonC = -121,960194;
$ DistC = 0,0548954278262;

/ *
#using автономная сфера
# если использовать эллипсоид, этот шаг немного отличается
# Преобразовать геодезическую широту / долготу в ECEF xyz
# 1. Преобразовать широту / долготу в радианы
# 2. Перевести лат / лонг (радианы) в 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)));

/ *
УСТАНОВИТЬ:
sudo pear установить Math_Vector-0.7.0
sudo pear установить Math_Matrix-0.8.7
* /
// Включить 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 (массив ($ xA, $ yA, $ zA));
$ P2vector = new Math_Vector3 (массив ($ xB, $ yB, $ zB));
$ P3vector = new Math_Vector3 (массив ($ xC, $ yC, $ zC));

# из википедии: http://en.wikipedia.org/wiki/Trlateration
# преобразовать, чтобы получить круг 1 в начале координат
# преобразовать, чтобы получить круг 2 на оси х

// CALC EX
$ P2minusP1 = Math_VectorOp :: substract ($ P2vector, $ P1vector);
$ l = новый Math_Vector ($ P2minusP1);
$ P2minusP1_length = $ l-> length ();
$ norm = new Math_Vector3 (массив ($ P2minusP1_length, $ P2minusP1_length, $ P2minusP1_length));
$ d = $ норма; // сохранить калькулятор D
$ ex = Math_VectorOp :: div ($ 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 (массив ($ 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 = новый Math_Vector3 (массив ($ 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 = новый Math_Vector ($ P3P1iex);
$ P3P1iex_length = $ l-> length ();
$ norm = new Math_Vector3 (массив ($ P3P1iex_length, $ P3P1iex_length, $ P3P1iex_length));
// echo "norm:". $ norm-> toString (). "\ n";
$ ey = Math_VectorOp :: split ($ P3P1iex, $ норма);
// 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 (массив ($ ey_x, $ ey_y, $ ey_z));

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

// CALC D
// сделать это раньше
$ d = floatval ($ d -> _ tuple-> getData () [0]);
// echo "d = $ d \ n";

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

# из википедии
#plug и chug, используя вышеуказанные значения
$ 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 / $ к) * $ х);

# здесь показан только один случай
$ z = sqrt (pow ($ DistA, 2) - pow ($ x, 2) - pow ($ y, 2));

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

#triPt - это массив с точками трилатерации x, y, z ECEF
$ 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]);


# конвертировать обратно в широту / долготу из ECEF
# конвертировать в градусы
$ lat = rad2deg (asin ($ triPt_z / $ earthR));
$ lon = rad2deg (atan2 ($ triPt_y, $ triPt_x));

echo $ lat. ','. $ LON;
fquinto
источник