Schnittpunkt zweier Kreise berechnen?

Ich versuche herauszufinden, wie man mathematisch die gemeinsamen Punkte von zwei sich schneidenden Kreisen auf der Erdoberfläche ableitet, wenn man einen Mittelpunkt Lat / Lon und einen Radius für jeden Punkt erhält.

Zum Beispiel gegeben:

• Lat / Lon (37,673442, -90,234036) Radius 107,5 NM
• Lat / Lon (36,109997, -90,953669) Radius 145 NM

Ich sollte zwei Schnittpunkte finden, von denen einer (36.948, -088.158) ist.

Es wäre trivial einfach, dies auf einer flachen Ebene zu lösen, aber ich habe keine Erfahrung damit, Gleichungen auf einer unvollkommenen Kugel wie der Erdoberfläche zu lösen.

Auf der Kugel ist es nicht viel schwieriger als im Flugzeug, wenn man das einmal erkennt

1. Die fraglichen Punkte sind die gegenseitigen Schnittpunkte von drei Kugeln: einer Kugel, die unter dem Ort x1 (auf der Erdoberfläche) eines gegebenen Radius zentriert ist, einer Kugel, die unter dem Ort x2 (auf der Erdoberfläche) eines gegebenen Radius zentriert ist, und der Erde selbst Dies ist eine Kugel, die bei O = (0,0,0) eines gegebenen Radius zentriert ist.

2. Der Schnittpunkt jeder der ersten beiden Kugeln mit der Erdoberfläche ist ein Kreis, der zwei Ebenen definiert. Die gegenseitigen Schnittpunkte aller drei Sphären liegen also auf dem Schnittpunkt dieser beiden Ebenen: einer Linie .

Infolgedessen reduziert sich das Problem darauf, eine Linie mit einer Kugel zu schneiden, was einfach ist.

Hier sind die Details. Die Eingaben sind Punkte P1 = (lat1, lon1) und P2 = (lat2, lon2) auf der Erdoberfläche, die als Kugel betrachtet werden, und zwei entsprechende Radien r1 und r2.

1. Konvertiert (lat, lon) in (x, y, z) geozentrische Koordinaten. Wie üblich, weil wir Maßeinheiten wählen können, in denen die Erde einen Einheitsradius hat,

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

   Im Beispiel hat P1 = (-90,234036 Grad, 37,673442 Grad) geozentrische Koordinaten x1 = (-0,00323306, -0,7915, 0,61116) und P2 = (-90,953669 Grad, 36,109997 Grad) geozentrische Koordinaten x2 = (-0,0134464, -0,807775 0,589337).

2. Wandeln Sie die Radien r1 und r2 (die entlang der Kugel gemessen werden) in Winkel entlang der Kugel um. Per Definition ist eine Seemeile (NM) 1/60 Bogengrad (was pi / 180 * 1/60 = 0,0002908888 Bogenmaß ist). Daher als Winkel,

   r1 = 107.5 / 60 Degree = 0.0312705 radian
   r2 = 145 / 60 Degree = 0.0421788 radian
   

3. Der geodätische Kreis mit dem Radius r1 um x1 ist der Schnittpunkt der Erdoberfläche mit einer bei cos (r1) * x1 zentrierten euklidischen Kugel mit dem Radius sin (r1).

4. Die Ebene, die durch den Schnittpunkt der Kugel mit dem Radius sin (r1) um cos (r1) * x1 und der Erdoberfläche bestimmt wird, verläuft senkrecht zu x1 und durch den Punkt cos (r1) x1, dessen Gleichung x.x1 = cos lautet (r1) (das "." stellt das übliche Skalarprodukt dar ); ebenfalls für das andere Flugzeug. Es wird einen eindeutigen Punkt x0 auf dem Schnittpunkt dieser beiden Ebenen geben, der eine lineare Kombination von x1 und x2 ist. Wenn Sie x0 = a x1 + b * x2 schreiben, sind die beiden ebenen Gleichungen

   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
   

   Unter Verwendung der Tatsache, dass x2.x1 = x1.x2, die ich als q schreiben werde, wird die Lösung (falls vorhanden) durch gegeben

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

   Im laufenden Beispiel berechne ich a = 0,973503 und b = 0,0260194.

   Offensichtlich brauchen wir q ^ 2! = 1. Dies bedeutet, dass x1 und x2 weder der gleiche Punkt noch Antipodenpunkte sein können.

5. Nun unterscheiden sich alle anderen Punkte auf der Schnittlinie der beiden Ebenen von x0 um ein Vielfaches eines Vektors n, der senkrecht zu beiden Ebenen steht. Das Kreuzprodukt

   n = x1~Cross~x2
   

   erledigt die Aufgabe, vorausgesetzt, n ist ungleich Null: Dies bedeutet wiederum, dass x1 und x2 weder zusammenfallen noch sich diametral gegenüberliegen. (Wir müssen darauf achten, das Kreuzprodukt mit hoher Genauigkeit zu berechnen, da es Subtraktionen mit viel Löschung beinhaltet, wenn x1 und x2 nahe beieinander liegen.) Im Beispiel ist n = (0.0272194, -0.00631254, -0.00803124) .

6. Deshalb suchen wir bis zu zwei Punkte der Form x0 + t * n, die auf der Erdoberfläche liegen: Das heißt, ihre Länge ist gleich 1. Entsprechend ist ihre quadratische Länge 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
   

   Der Term mit x0.n verschwindet, weil x0 (eine lineare Kombination von x1 und x2) senkrecht zu n ist. Die beiden Lösungen sind leicht

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

   und es ist negativ. Auch hier ist hohe Präzision gefragt, denn wenn x1 und x2 nahe beieinander liegen, liegt x0.x0 sehr nahe bei 1, was zu einem gewissen Verlust an Gleitkommapräzision führt. Im Beispiel ist t = 1,07509 oder t = -1,07509. Die beiden Schnittpunkte sind also gleich

   x0 + t*n = (0.0257661, -0.798332, 0.601666)
   x0 - t*n = (-0.0327606, -0.784759, 0.618935)
   

7. Schließlich können wir diese Lösungen zurück in (lat, lon) konvertieren, indem wir geozentrische (x, y, z) in geografische Koordinaten konvertieren:

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

   Verwenden Sie für den Längengrad den generalisierten Arkustangens-Rückgabewert im Bereich von -180 bis 180 Grad (in Computeranwendungen verwendet diese Funktion sowohl x als auch y als Argumente und nicht nur das Verhältnis y / x; manchmal wird sie als "ATan2" bezeichnet).

   Ich erhalte die beiden Lösungen (-88.151426, 36.989311) und (-92.390485, 38.238380), die in der Abbildung als gelbe Punkte dargestellt sind.

Die Achsen zeigen die geozentrischen Koordinaten (x, y, z) an. Der graue Fleck ist der Teil der Erdoberfläche von -95 bis -87 Grad Länge, 33 bis 40 Grad Breite (markiert mit einem Grad-Gitter). Die Erdoberfläche wurde teilweise transparent gemacht, um alle drei Kugeln zu zeigen. Die Richtigkeit der berechneten Lösungen wird dadurch deutlich, wie die gelben Punkte an den Schnittpunkten der Kugeln liegen.

Der ellipsoidale Fall:

Dieses Problem ist eine Verallgemeinerung des Auffindens von Seegrenzen, die als "Mittellinien" definiert sind, und es gibt eine umfangreiche Literatur zu diesem Thema. Meine Lösung für dieses Problem besteht darin, die äquidistante azimutale Projektion zu nutzen:

1. Errate den Schnittpunkt
2. Projizieren Sie die beiden Basispunkte mit diesem erratenen Schnittpunkt als Mittelpunkt einer äquidistanten Azimutalprojektion.
3. Lösen Sie das Überschneidungsproblem im 2d projizierten Raum.
4. Wenn der neue Schnittpunkt zu weit vom alten entfernt ist, kehren Sie zu Schritt 2 zurück.

Dieser Algorithmus konvergiert quadratisch und liefert eine genaue Lösung auf einem Ellipsoid. (Bei Seegrenzen ist Genauigkeit erforderlich, da hierdurch die Fischerei-, Öl- und Mineralienrechte festgelegt werden.)

Die Formeln für ein Rotationsellipsoid sind in Abschnitt 14 der Geodäten angegeben . Die ellipsoidale azimutale Projektion mit gleichem Abstand wird von GeographicLib bereitgestellt . Eine MATLAB-Version ist unter Geodätische Projektionen für ein Ellipsoid verfügbar .

Hier ist ein R-Code, um dies zu tun:

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)

Nach der Antwort von @ whuber folgt ein Java-Code, der aus zwei Gründen nützlich ist:

• Es wird eine Zusammenfassung zu ArcTan angezeigt (für Java und möglicherweise für andere Sprachen?)
• es behandelt die möglichen Randfälle, einschließlich eines in @ whubers Antwort nicht erwähnten.

Es ist nicht optimiert oder vollständig (ich habe offensichtliche Klassen wie ausgelassen Point), sollte aber den Trick machen.

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

Beachten Sie auch die Verwendung von atan2- es ist das Gegenteil von dem, was Sie von der Antwort von @ whuber erwarten würden (ich weiß nicht warum, aber es funktioniert):

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

Working 'R' Code für @wuhber Antwort.

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=",")

Wenn einer der Kreise der Nortstar ist, gibt es den einfachsten Weg mit der Einheitenkugel.

Sie können Ihren Breitengrad mit Nortstar messen. Dann haben Sie eine relative Position auf dieser Kugel. v1 (0, sin (la), cos (la)) Sie kennen die Position (Winkel) eines anderen Sterns (star2) aus Almanach. v2 (sin (lo2) * cos (la2), sin (la2), cos (lo2) * cos (la2)) Seine Vektoren. Aus Kugelgleichung.

lo2 ist die relative Länge. Es ist unbekannt .

Der Winkel zwischen dir und star2, den du auch messen kannst, (m) Und du weißt, das innere Produkt von zwei Einheitsvektoren ist cos (Winkel) von dazwischen. cos (m) = Punkt (v1, v2) u kann nun die relative Länge (lo2) berechnen. lo2 = acos ((cos (m) -sin (la) · sin (la2)) / (cos (la) · cos (la2)))

Schließlich addieren Sie die wahre Länge von star2 zu lo2. (oder U-Boot, hängt von der Westseite von Ihnen ab oder von der Ostseite.) lo2 ist jetzt Ihre Länge.

Entschuldigung für mein Englisch, ich lerne diese Sprache nie.

2 Dinge: Northstar bedeutet Polstar.

Ein weiterer. Da der Winkel relativ zum Horizont gemessen wird, muss immer eine 90-Winkel-Korrektur durchgeführt werden. Es gilt auch für m-Winkel.

ps: realer Winkel Mittelwert: Sternposition - Zeitkorrektur.