Obliczanie przecięcia dwóch kół?

29

Próbuję dowiedzieć się, jak matematycznie wyprowadzić wspólne punkty dwóch przecinających się kół na powierzchni Ziemi, biorąc pod uwagę środek Lat / Lon i promień dla każdego punktu.

Na przykład biorąc pod uwagę:

  • Lat / Lon (37,673442, -90,234036) Promień 107,5 NM
  • Lat / Lon (36.109997, -90.953669) Promień 145 NM

Powinienem znaleźć dwa punkty przecięcia z jednym z nich (36,948, -088.158).

Rozwiązanie problemu na płaskiej płaszczyźnie byłoby banalnie proste, ale nie mam doświadczenia w rozwiązywaniu równań na niedoskonałej kuli, takiej jak powierzchnia ziemi.

Wola
źródło
1
Jeśli wszystkie twoje promienie będą tak małe (mniej niż kilka kilometrów), wówczas ziemia jest zasadniczo płaska w tej skali i równie dobrze możesz wybrać dokładną, prostą projekcję i wykonać zwykłe obliczenia euklidesowe. Upewnij się, że obliczasz przecięcie z więcej niż trzema miejscami po przecinku - niedokładność w trzecim miejscu po przecinku jest tak duża, jak którykolwiek z twoich promieni!
whuber
1
Powinienem był dodać jednostki, te promienie są w NM, więc nadal jest niewielka odległość w stosunku do powierzchni ziemi, ale większa niż kilka kilometrów. Jak ta skala wpływa na zniekształcenie? Próbuję znaleźć rozwiązanie o dokładności mniejszej niż <1 nm, więc nie musi to być super precyzja. Dzięki!
Czy
To wszystko dobrze wiedzieć, ponieważ pokazuje, że możesz użyć sferycznego modelu ziemi - bardziej skomplikowane modele elipsoidalne są niepotrzebne.
whuber
@ whuber Czy to sugeruje, że problem można rozwiązać ponownie: znaleźć przecięcie 3 sfer, gdzie jedna z nich jest ziemią, a pozostałe dwie są wyśrodkowane na punktach o odpowiednich promieniach?
Kirk Kuykendall
@Kirk Tak, to jest sposób, aby to zrobić, zakładając sferyczny model powierzchni ziemi. Po wstępnych obliczeniach, które sprowadzają to do specjalnego przypadku problemu trilateracji w 3D. (Obliczenia są potrzebne do przeliczenia odległości wzdłuż łuków sferycznych na odległości wzdłuż cięciw sferycznych, które stają się promieniami dwóch mniejszych sfer.)
whuber

Odpowiedzi:

21

Gdy tylko to rozpoznasz, nie będzie to trudniejsze na kuli niż w samolocie

  1. Chodzi o wzajemne przecięcia trzech sfer: kuli wyśrodkowanej poniżej położenia x1 (na powierzchni ziemi) o danym promieniu, kuli wyśrodkowanej poniżej położenia x2 (na powierzchni ziemi) o danym promieniu i samej ziemi , która jest kulą wyśrodkowaną na O = (0,0,0) o danym promieniu.

  2. Przecięcie każdej z dwóch pierwszych sfer z powierzchnią ziemi to okrąg, który definiuje dwie płaszczyzny. Wzajemne przecięcia wszystkich trzech sfer leżą zatem na przecięciu tych dwóch płaszczyzn: linii .

W rezultacie problem sprowadza się do przecięcia linii ze sferą, co jest łatwe.


Oto szczegóły. Wejściami są punkty P1 = (lat1, lon1) i P2 = (lat2, lon2) na powierzchni ziemi, uważane za kulę, i dwa odpowiadające im promienie r1 i r2.

  1. Konwertuj (lat, lon) na (x, y, z) współrzędne geocentryczne. Jak zwykle, ponieważ możemy wybrać jednostki miary, w których Ziemia ma promień jednostki,

    x = cos(lon) cos(lat)
    y = sin(lon) cos(lat)
    z = sin(lat).
    

    W przykładzie P1 = (-90,234036 stopnia, 37,673442 stopień) ma współrzędne geocentryczne x1 = (-0,00323306, -0,7915, 0,61116), a P2 = (-90.953669 stopień, 36,109997 stopnia) ma współrzędne geocentryczne x2 = (-0,0134464, -0,8077 , 0,589337).

  2. Przelicz promienie r1 i r2 (mierzone wzdłuż kuli) na kąty wzdłuż kuli. Z definicji jedna mila morska (NM) to 1/60 stopnia łuku (co oznacza pi / 180 * 1/60 = 0,0002908888 radianów). Dlatego jako kąty

    r1 = 107.5 / 60 Degree = 0.0312705 radian
    r2 = 145 / 60 Degree = 0.0421788 radian
    
  3. Geodezyjnej okrąg o promieniu r1 około X1 przecięcie powierzchni Ziemi z euklidesowej kuli o promieniu sin (R1) wyśrodkowany cos (R1) * X1.

  4. Płaszczyzna wyznaczona przez przecięcie kuli o promieniu sin (r1) wokół cos (r1) * x1, a powierzchnia ziemi jest prostopadła do x1 i przechodzi przez punkt cos (r1) x1, skąd jej równanie wynosi x.x1 = cos (r1) („.” oznacza zwykły iloczyn skalarny ); podobnie w przypadku drugiej płaszczyzny. Na przecięciu tych dwóch płaszczyzn będzie unikalny punkt x0, który jest liniową kombinacją x1 i x2. Zapisując x0 = a x1 + b * x2 są dwa równania płaskie

    cos(r1) = x.x1 = (a*x1 + b*x2).x1 = a + b*(x2.x1)
    cos(r2) = x.x2 = (a*x1 + b*x2).x2 = a*(x1.x2) + b
    

    Biorąc pod uwagę fakt, że x2.x1 = x1.x2, które napiszę jako q, rozwiązanie (jeśli istnieje) jest podane przez

    a = (cos(r1) - cos(r2)*q) / (1 - q^2),
    b = (cos(r2) - cos(r1)*q) / (1 - q^2).
    

    W bieżącym przykładzie obliczam a = 0,973503 ib = 0,0260194.

    Najwyraźniej potrzebujemy q ^ 2! = 1. Oznacza to, że x1 i x2 nie mogą być ani tym samym punktem, ani antypodalnymi punktami.

  5. Teraz wszystkie inne punkty na linii przecięcia dwóch płaszczyzn różnią się od x0 pewną wielokrotnością wektora n, który jest wzajemnie prostopadły do ​​obu płaszczyzn. Produkt krzyżowy

    n = x1~Cross~x2
    

    wykonuje zadanie pod warunkiem, że n jest niezerowe: ponownie oznacza to, że x1 i x2 nie są ani zbieżne ani diametralnie przeciwne. (Musimy zadbać o obliczenie iloczynu krzyżowego z wysoką precyzją, ponieważ obejmuje on odejmowanie z dużą ilością anulowania, gdy x1 i x2 są blisko siebie.) W przykładzie n = (0,0272194, -0,00631254, -0,00803124) .

  6. Dlatego szukamy do dwóch punktów formy x0 + t * n, które leżą na powierzchni ziemi: to znaczy, że ich długość wynosi 1. Odpowiednio ich kwadratowa długość wynosi 1:

    1 = squared length = (x0 + t*n).(x0 + t*n) = x0.x0 + 2t*x0.n + t^2*n.n = x0.x0 + t^2*n.n
    

    Termin z x0.n znika, ponieważ x0 (będący liniową kombinacją x1 i x2) jest prostopadły do ​​n. Oba rozwiązania łatwo są

    t = sqrt((1 - x0.x0)/n.n)
    

    i jego negatywne. Ponownie wymagana jest wysoka precyzja, ponieważ gdy x1 i x2 są blisko, x0.x0 jest bardzo bliskie 1, co prowadzi do pewnej utraty precyzji zmiennoprzecinkowej. W przykładzie t = 1,07509 lub t = -1,07509. Dwa punkty przecięcia są zatem równe

    x0 + t*n = (0.0257661, -0.798332, 0.601666)
    x0 - t*n = (-0.0327606, -0.784759, 0.618935)
    
  7. Na koniec możemy przekonwertować te rozwiązania z powrotem na (lat, lon), przekształcając dane geocentryczne (x, y, z) na współrzędne geograficzne:

    lon = ArcTan(x,y)
    lat = ArcTan(Sqrt[x^2+y^2], z)
    

    Dla długości użyć uogólnionego tangens powrocie wartości w zakresie od -180 do 180 stopni (w zastosowaniach komputerowych, funkcja wykonuje zarówno X i Y, argumentów, a nie tylko w stosunku y / x, jest czasami nazywana „ATAN2”).

    Otrzymuję dwa rozwiązania (-88,151426, 36,989311) i (-92,390485, 38,238380), pokazane na rysunku jako żółte kropki.

Rysunek 3D

Osie wyświetlają współrzędne geocentryczne (x, y, z). Szara łata to część powierzchni ziemi od -95 do -87 stopni długości geograficznej, od 33 do 40 stopni szerokości geograficznej (zaznaczona siatką o jeden stopień). Powierzchnia ziemi została częściowo przezroczysta, aby pokazać wszystkie trzy sfery. Poprawność obliczonych rozwiązań jest widoczna dzięki temu, jak żółte punkty znajdują się na przecięciach sfer.

Whuber
źródło
Bill, to jest niesamowite. Jedno wyjaśnienie, które możesz dodać, oparte na osobie, która próbowała go wdrożyć. W kroku 2 nie podajesz jawnie przeliczenia stopni na radiany.
Jersey Andy,
@Jersey Dziękujemy za sugerowaną edycję. Zmieniłem to trochę, aby uniknąć nadmiarowości i zachować formuły tak jasne, jak to możliwe. Po przeczytaniu wątku, do którego się odwołujesz, dodałem również link wyjaśniający iloczyn skalarny.
whuber
8

Elipsoidalny przypadek:

Problem ten jest uogólnieniem problemu znajdowania granic morskich określonych jako „linie środkowe” i istnieje obszerna literatura na ten temat. Moim rozwiązaniem tego problemu jest wykorzystanie równoodległej projekcji azymutalnej:

  1. Zgadnij w punkcie przecięcia
  2. Rzutuj dwa punkty bazowe, używając tego odgadniętego punktu przecięcia jako środka równoodległego rzutu azymutalnego,
  3. Rozwiąż problem przecięcia w rzutowanej przestrzeni 2d.
  4. Jeśli nowy punkt przecięcia znajduje się za daleko od starego, wróć do kroku 2.

Algorytm ten zbiega się kwadratowo i daje dokładne rozwiązanie dla elipsoidy. (Dokładność jest wymagana w przypadku granic morskich, ponieważ określa prawa do połowów, ropy naftowej i minerałów).

Wzory podano w rozdziale 14 Geodezji na elipsoidzie obrotu . Elipsoidalna równoodległa projekcja azymutalna jest dostarczana przez GeographicLib . Wersja MATLAB jest dostępna w rzutach geodezyjnych dla elipsoidy .

cffk
źródło
+1 To niesamowity artykuł: twój skromny opis tutaj nie oddaje sprawiedliwości.
whuber
Zobacz także mój krótszy artykuł na temat geodezji „Algorytmy dla geodetyki” dx.doi.org/10.1007/s00190-012-0578-z (do pobrania bezpłatnie!) Oraz erratę i dodatki do tych artykułów geographiclib.sf.net/geod-addenda.html
cffk,
1

Oto trochę kodu R, aby to zrobić:

p1 <- cbind(-90.234036, 37.673442) 
p2 <- cbind(-90.953669, 36.109997 )

library(geosphere)
steps <- seq(0, 360, 0.1)
c1 <- destPoint(p1, steps, 107.5 * 1852)
c2 <- destPoint(p2, steps, 145 * 1852)

library(raster)
s1 <- spLines(c1)
s2 <- spLines(c2)

i <- intersect(s1, s2)
coordinates(i)

#        x        y
# -92.38241 38.24267
# -88.15830 36.98740

s <- bind(s1, s2)
crs(s) <- "+proj=longlat +datum=WGS84"
plot(s)
points(i, col='red', pch=20, cex=2)
Robert Hijmans
źródło
1

Po odpowiedzi na @ whuber , oto kod Java, który jest użyteczny z dwóch powodów:

  • podkreśla gotcha dotyczącą ArcTan (dla Javy i może innych języków?)
  • obsługuje możliwe przypadki krawędzi, w tym jeden nie wymieniony w odpowiedzi @ whuber.

Nie jest zoptymalizowany ani ukończony (pominąłem oczywiste klasy takie jak Point), ale powinien załatwić sprawę.

public static List<Point> intersection(EarthSurfaceCircle c1, EarthSurfaceCircle c2) {

    List<Point> intersections = new ArrayList<Point>();

    // project to (x,y,z) with unit radius
    UnitVector x1 = UnitVector.toPlanar(c1.lat, c1.lon);
    UnitVector x2 = UnitVector.toPlanar(c2.lat, c2.lon);

    // convert radii to radians:
    double r1 = c1.radius / RadiusEarth;
    double r2 = c2.radius / RadiusEarth;

    // compute the unique point x0
    double q = UnitVector.dot(x1, x2);
    double q2 = q * q;
    if (q2 == 1) {
        // no solution: circle centers are either the same or antipodal
        return intersections;
    }
    double a = (Math.cos(r1) - q * Math.cos(r2)) / (1 - q2);
    double b = (Math.cos(r2) - q * Math.cos(r1)) / (1 - q2);
    UnitVector x0 = UnitVector.add(UnitVector.scale(x1, a), UnitVector.scale(x2, b));

    // we only have a solution if x0 is within the sphere - if not,
    // the circles are not touching.
    double x02 = UnitVector.dot(x0, x0);
    if (x02 > 1) {
        // no solution: circles not touching
        return intersections;
    }

    // get the normal vector:
    UnitVector n = UnitVector.cross(x1, x2);
    double n2 = UnitVector.dot(n, n);
    if (n2 == 0) {
        // no solution: circle centers are either the same or antipodal
        return intersections;
    }

    // find intersections:
    double t = Math.sqrt((1 - UnitVector.dot(x0, x0)) / n2);
    intersections.add(UnitVector.toPolar(UnitVector.add(x0, UnitVector.scale(n, t))));
    if (t > 0) {
        // there's only multiple solutions if t > 0
        intersections.add(UnitVector.toPolar(UnitVector.add(x0, UnitVector.scale(n, -t))));
    }
    return intersections;
}

Co ważne, zwróć uwagę na użycie atan2- jest to odwrotność tego, czego można oczekiwać od odpowiedzi @ whubera (nie wiem dlaczego, ale działa):

    public static Point toPolar(UnitVector a) {
        return new Point(
                Math.toDegrees(Math.atan2(a.z, Math.sqrt(a.x * a.x + a.y * a.y))),
                Math.toDegrees(Math.atan2(a.y, a.x)));          
    }
ianhoolihan
źródło
0

Działający kod „R” dla odpowiedzi @wuhber.

P1 <- c(37.673442, -90.234036)
P2 <- c(36.109997, -90.953669) 

#1 NM nautical-mile is 1852 meters
R1 <- 107.5
R2 <- 145

x1 <- c(
  cos(deg2rad(P1[2])) * cos(deg2rad(P1[1])),  
  sin(deg2rad(P1[2])) * cos(deg2rad(P1[1])),
  sin(deg2rad(P1[1]))
);

x2 <- c(
  cos(deg2rad(P2[2])) * cos(deg2rad(P2[1])),
  sin(deg2rad(P2[2])) * cos(deg2rad(P2[1])),
  sin(deg2rad(P2[1]))
);

r1 = R1 * (pi/180) * (1/60)
r2 = R2 * (pi/180) * (1/60)

q = dot(x1,x2)
a = (cos(r1) - cos(r2) * q) / (1 - q^2)
b = (cos(r2) - cos(r1) * q)/ (1 - q^2)

n <- cross(x1,x2)

x0 = a*x1 + b*x2


t = sqrt((1 - dot(x0, x0))/dot(n,n))

point1 = x0 + (t * n)
point2 = x0 - (t * n)

lat1 = rad2deg(atan2(point1[2] ,point1[1]))
lon1= rad2deg(asin(point1[3]))
paste(lat1, lon1, sep=",")

lat2 = rad2deg(atan2(point2[2] ,point2[1]))
lon2 = rad2deg(asin(point2[3]))
paste(lat2, lon2, sep=",")
Sri
źródło
-1

Jeśli jednym z kół jest Nortstar, to najłatwiejszym sposobem jest użycie kuli jednostkowej.

Możesz zmierzyć swoją szerokość geograficzną za pomocą Nortstar. Następnie masz względną pozycję na tej sferze. v1 (0, sin (la), cos (la)) Znasz pozycję (kąt) innej gwiazdy (star2), z almanach. v2 (sin (lo2) * cos (la2), sin (la2), cos (lo2) * cos (la2)) Jego wektory. Z równania sfery.

lo2 jest względną długością geograficzną. Nie wiadomo .

Kąt między tobą a gwiazdą 2 również możesz zmierzyć, (m) I wiesz, iloczyn wewnętrzny wektora dwóch jednostek to cos (kąt) pomiędzy. cos (m) = kropka (v1, v2) możesz teraz obliczyć względną długość geograficzną (lo2). lo2 = acos ((cos (m) -sin (la) * sin (la2)) / (cos (la) * cos (la2)))

W końcu dodajesz prawdziwą długość geograficzną star2 do lo2. (lub okręt podwodny, zależą od jego zachodniej strony od ciebie lub wschodniej.) lo2 teraz jest twoją długością geograficzną.

Przepraszam za mój angielski, nigdy nie uczę się tego języka.


2 rzeczy: Northstar oznacza gwiazdę polarną.

Inne. Ponieważ kąt mierzony względem poziomu jest konieczny, zawsze potrzebna jest korekcja 90-kątowa. Obowiązuje również kąt m.

ps: średni kąt rzeczywisty: pozycja gwiazdy - korekta czasowa.

docwho
źródło
Nie wiadomo, jak to odpowiada na pytanie.
whuber