Trilaterasi menggunakan 3 titik lintang / bujur, dan 3 jarak?

34

Saya ingin mengetahui lokasi target yang tidak diketahui (koordinat lintang dan bujur). Ada 3 titik yang diketahui (pasangan koordinat lintang dan bujur) dan untuk setiap titik jarak dalam kilometer ke lokasi target. Bagaimana saya bisa menghitung koordinat lokasi target?

Sebagai contoh, katakan saya memiliki poin data berikut

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

Apa yang saya suka adalah apa matematika untuk fungsi yang menganggapnya sebagai input dan mengembalikan 37.417959, -121.961954 sebagai output.

Saya mengerti bagaimana menghitung jarak antara dua titik, dari http://www.movable-type.co.uk/scripts/latlong.html Saya mengerti prinsip umum bahwa dengan tiga lingkaran seperti ini Anda mendapatkan tepat satu titik tumpang tindih. Yang saya kabur adalah matematika yang diperlukan untuk menghitung titik itu dengan input ini.

nohat
sumber
Berikut adalah halaman yang menuntun Anda melalui matematika untuk menemukan pusat tiga koordinat. Mungkin ini bisa membantu. < mathforum.org/library/drmath/view/68373.html >
Jon Bringhurst
1
Apakah ini perlu di sphere / spheroid, atau apakah algoritma planar baik-baik saja?
fmark
1
Saya tidak bisa memberi Anda jawabannya, tetapi saya pikir saya bisa mengarahkan Anda ke arah yang benar. Tiga koordinat = tiga titik pusat. Tiga jarak = tiga lingkaran. Dua lingkaran yang berpotongan, dapat memiliki kemungkinan tidak ada / satu / dua solusi. Tiga Lingkaran tidak dapat memiliki / satu / atau suatu area sebagai solusinya. Dapatkan rumus lingkaran untuk tiga lingkaran, dan pecahkan dengan Sistem Persamaan / Aljabar.
CrazyEnigma
Sebenarnya, Anda bahkan tidak memerlukan sistem untuk menyelesaikannya. Ada satu atau dua kemungkinan, tetapi karena Anda memiliki nilai jarak, Anda dapat memisahkan jawaban yang benar.
George Silva
1
+1 Ini adalah pertanyaan yang bagus. Pada awalnya saya menemukan solusi yang dapat dengan mudah ditemukan dengan google, tetapi ternyata tidak. Mungkin masalahnya bisa lebih umum dinyatakan: diberi N poin dengan setiap titik tidak hanya memiliki jarak tetapi juga margin kesalahan, menemukan elips kepercayaan.
Kirk Kuykendall

Jawaban:

34

Setelah beberapa melihat-lihat Wikipedia dan pertanyaan / jawaban yang sama di StackOverflow , saya pikir saya akan menikamnya, dan mencoba mengisi kekosongan.

Pertama, Tidak yakin dari mana Anda mendapatkan output, tetapi tampaknya salah. Saya memplot titik-titik di ArcMap, buffer mereka ke jarak yang ditentukan, berlari berpotongan ke pada buffer, dan kemudian menangkap titik persimpangan untuk mendapatkan solusi. Output yang Anda usulkan adalah titik berwarna hijau. Saya menghitung nilai di kotak info, yaitu sekitar 3 meter dari apa yang ArcMap berikan untuk solusi yang berasal dari intersect.

teks alternatif

Matematika di halaman wikipedia tidak terlalu buruk, hanya perlu menyelaraskan koordinat geodetik Anda dengan kartesian ECEF, yang dapat ditemukan sini . istilah a / x + h dapat digantikan oleh jari-jari bola authalic, jika Anda tidak menggunakan ellipsoid.

Mungkin yang paling mudah hanya memberi Anda beberapa kode yang didokumentasikan dengan baik (?), Jadi ini dia dengan 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
Wickick
sumber
1
Saya akan mengumpulkan jawaban yang sama, tetapi sekarang tidak perlu! Mendapat upvote saya.
Wrass
numpy untuk menyelamatkan! Ia mengkompilasi ketika 'triPt' diganti dengan 'triLatPt', tetapi sebaliknya mengembalikan 37.4191023738 -121.960579208.
Kerja
Kerja bagus! Jika saya mengganti sistem koordinat geografis ke sistem koordinat [Cartesian] lokal, apakah ini masih berfungsi?
zengr
bagi mereka yang ada di domain c ++ ..hacked bersama yang sangat cepat pastebin.com/9Dur6RAP
raaj
2
Terima kasih @wickick! Saya telah mem-porting ini ke JavaScript (dimaksudkan untuk Node tetapi dapat dengan mudah dikonversi agar berfungsi di browser). gist.github.com/dav-/bb7103008cdf9359887f
DC_
6

Saya tidak yakin apakah saya naif, tetapi, jika Anda buffer setiap titik berdasarkan ukuran, dan kemudian memotong ketiga lingkaran yang akan membuat Anda mendapatkan lokasi yang benar?

Anda dapat menghitung persimpangan menggunakan API spasial. Contoh:

  • GeoScript
  • Suite Java Topology
  • Suite Topologi NET
  • GEOS
George Silva
sumber
1
Tepatnya, dia tertarik pada formula untuk mendapatkan titik persimpangan itu.
Vinko Vrsalovic
Menggunakan API spasial Anda dapat melakukannya tanpa menggunakan matematika murni.
George Silva
1
@ George dapatkah Anda memberikan contoh API seperti itu?
nohat
Posting yang diedit untuk mencerminkan permintaan nohat.
George Silva
+1, pemikiran lateral yang baik, bahkan jika bukan yang paling efisien secara komputasi!
fmark
2

Catatan berikut menggunakan geometri planaritmik (yaitu Anda harus memproyeksikan koordinat Anda ke sistem koordinat lokal yang sesuai).

Alasan saya, dengan contoh yang berhasil di Python, mengikuti:

Ambil 2 titik data (panggil mereka adan b). Hubungi titik target kami x. Kami sudah tahu jarak axdan bx. Kita dapat menghitung jarak abmenggunakan teorema Pythagoras.

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

Sekarang, Anda bisa menentukan sudut garis-garis ini:

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

Sayangnya saya kekurangan waktu untuk menyelesaikan jawaban untuk Anda, Namun, sekarang Anda tahu sudutnya, Anda dapat menghitung dua lokasi yang mungkin untuk x. Kemudian, menggunakan titik ketiga c Anda dapat menghitung lokasi mana yang benar.

fmark
sumber
2

Ini mungkin berhasil. Dengan cepat lagi dalam python, Anda bisa meletakkan ini di tubuh fungsi xN, yN = koordinat titik, r1 & r2 = nilai radius

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 nilai adalah nilai kembali (harus dalam array) dari dua titik persimpangan pada lingkaran, jika itu membantu memperjelas hal-hal.

Lakukan ini untuk 2 lingkaran pertama, lalu lagi untuk yang pertama dan terakhir. Jika salah satu hasil dari iterasi pertama dibandingkan dengan hasil dari iterasi kedua (dalam beberapa toleransi, mungkin), maka Anda memiliki titik persimpangan. Ini bukan solusi yang bagus terutama ketika Anda mulai menambahkan lebih dari beberapa poin ke dalam proses, tetapi merupakan yang paling sederhana yang dapat saya lihat tanpa harus menyelesaikan sistem persamaan.

WolfOdrade
sumber
Apa 'e' dan 'k' dalam kode Anda?
ReinierDG
Saya tidak ingat :-) Jawaban wwnick lebih sesuai dengan sesuatu yang ingin Anda terapkan jika Anda hanya memiliki tiga lingkaran.
WolfOdrade
1

Anda dapat menggunakan API spasial dari postgis (fungsi St_Intersection, St_buffer). Seperti yang diperhatikan oleh fmark, Anda juga harus ingat bahwa Postgis menggunakan algoritma planar, tetapi untuk area kecil, menggunakan proyeksi equi-jauh tidak menimbulkan banyak kesalahan.

stachu
sumber
PostGIS dapat melakukan perhitungan spheroid menggunakan GEOGRAPHYtipe daripada GEOMETRYtipe.
fmark
1

Lakukan dalam bahasa PHP:

// mengasumsikan ketinggian = 0
$ earthR = 6371; // dalam km (= 3959 mil)

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

/ *
#menggunakan bola authalic
#Jika menggunakan ellipsoid, langkah ini sedikit berbeda
# Konversi geodetik Lat / Panjang ke ECEF xyz
# 1. Konversi Lat / Panjang ke radian
# 2. Konversi Lat / Panjang (radian) ke 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)));

/ *
MEMASANG:
sudo pear instal Math_Vector-0.7.0
sudo pear instal Math_Matrix-0.8.7
* /
// Sertakan 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));

#dari wikipedia: http://en.wikipedia.org/wiki/Trilateration
#transformasi untuk mendapatkan lingkaran 1 di tempat asalnya
#transformasi untuk mendapatkan lingkaran 2 pada sumbu x

// CALC EX
$ P2minusP1 = Math_VectorOp :: subtract ($ P2vector, $ P1vector);
$ l = new Math_Vector ($ P2minusP1);
$ P2minusP1_length = $ l-> length ();
$ norm = new Math_Vector3 (array ($ P2minusP1_length, $ P2minusP1_length, $ P2minusP1_length));
$ d = $ norm; // simpan 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 :: subtract ($ 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 :: subtract ($ 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
// lakukan sebelumnya
$ d = floatval ($ d -> _ tuple-> getData () [0]);
// echo "d = $ d \ n";

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

#dari wikipedia
# colok dan cabut menggunakan nilai di atas
$ 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);

# hanya satu kasing yang ditunjukkan di sini
$ z = sqrt (pow ($ DistA, 2) - pow ($ x, 2) - pow ($ y, 2));

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

#triPt adalah array dengan ECEF x, y, z dari titik trilateration
$ xex = Math_VectorOp :: scale ($ x, $ ex);
$ yey = Math_VectorOp :: scale ($ y, $ ey);
$ zez = Math_VectorOp :: scale ($ z, $ ez);

// CALC $ triPt = $ P1vektor + $ 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 kembali ke lat / long dari ECEF
#konversi ke derajat
$ lat = rad2deg (asin ($ triPt_z / $ earthR));
$ lon = rad2deg (atan2 ($ triPt_y, $ triPt_x));

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