Wie berechnet man den Breiten- und Längengrad eines Punktes auf der Erde mit Hilfe von Azimut und Höhe, gesehen von einem geostationären Satelliten?

Hinweis: Ich glaube, die Erdkoordinaten sollten mit Großbuchstaben geschrieben werden$(X, Y, Z)$aber Zitate aus der anderen Antwort verwenden Kleinbuchstaben$(x, y, z)$. Da der radikale Begriff auch kaum zu Kleinbuchstaben passt, habe ich an den meisten Stellen die Kleinschreibung beibehalten. Bei dieser Antwort stehen Groß- und Kleinschreibung für dieselben Koordinaten.

tl; dr: Stecken Sie die Position Ihres Satelliten ein ($x, y, z$)

$$\mathbf{x_0} = x\mathbf{\hat{x}} + y\mathbf{\hat{y}} + z\mathbf{\hat{z}}$$ und die Normale des Vektors ($u, v, w$) $$\mathbf{\hat{n}} = u\mathbf{\hat{x}} + v\mathbf{\hat{y}} + w\mathbf{\hat{z}}$$ zu bekommen$t$die Länge des Vektors vom Satelliten zum Ellipsoid, dann verwenden $$\mathbf{x_{intersect}} = \mathbf{x_0} + t\mathbf{\hat{n}}$$ um die Koordinaten des Schnittpunkts mit dem Referenzellipsoid zu erhalten.

Verwenden Sie dann entweder die analytische Lösung im zweiten Abschnitt unten oder eine der iterativen Lösungen der Referenzen dort, um diese in Längen- und Breitengrad umzuwandeln.

---

Schnittpunkt-Mathematik, um (x, y, z) zu erhalten

Aus dieser praktischen Antwort in GIS SE:

Das WGS84-Referenzellipsoid ist ein zweiachsiges (und abgeflachtes) Ellipsoid. Er ist an den Polen kürzer als der Äquator, und der Äquator ist ein Kreis.

Die Gleichung dafür lautet:

$$\frac{x^2}{a^2} + \frac{y^2}{a^2} + \frac{z^2}{b^2} - 1 = 0$$

Wo$x, y, z$ein Punkt auf der Oberfläche des Ellipsoids ist, und

Wo$a$ist die große Halbachse (6.378.137 Meter) und$b$ist die kleine Halbachse des WGS84-Ellipsoids (6.356.752,3142 Meter).

Ich werde die Mathematik aus dieser Antwort weiterhin verwenden und in MathJax und später in Python neu formatieren:

Die Position Ihres Satelliten ist$x,y, z$und die Richtung normal ist$u, v, w$. Die Länge des Vektors vom Satelliten zum ersten und nächsten Schnittpunkt (vorausgesetzt natürlich, dass sich Ihr Satellit außerhalb der Erde befindet :

t = -(1/(b^2 (u^2 + v^2) +  a^2 w^2)) * (b^2 (u x + v y) + a^2 w z + 1/2 Sqrt[
     4 (b^2 (u x + v y) + a^2 w z)^2 - 
     4 (b^2 (u^2 + v^2) + a^2 w^2) (b^2 (-a^2 + x^2 + y^2) + a^2 z^2)])

In MathJax:

$$A = -\frac{1}{b^2(u^2+v^2) + a^2w^2}$$

$$B = b^2(ux + vy) + a^2wz$$

$$C = \sqrt{(b^2(ux + vy) + a^2wz)^2 - (b^2(u^2 + v^2) + a^2w^2) (b^2(-a^2 + x^2 + y^2) + a^2z^2)}$$

$$t = A(B+C)$$

was einfacher zu lesen ist als die Art und Weise, wie MathJax alles in einer Zeile anzeigt:

$$t = -\frac{1}{b^2(u^2+v^2) + a^2w^2} \left( b^2(ux + vy) + a^2wz + \sqrt{(b^2(ux + vy) + a^2wz)^2 - (b^2(u^2 + v^2) + a^2w^2) (b^2(-a^2 + x^2 + y^2) + a^2z^2)} \right)$$

Implementiert in Python:

import numpy as np
a, b = 0.9, 1.1
asq, bsq = a**2, b**2
x, y, z = 5.0, 0.0, 0.0
u, v, w = -np.sqrt(1 - 0.1**2 - 0.1**2), 0.1, 0.1
xsq, ysq, zsq = x**2, y**2, z**2
usq, vsq, wsq = u**2, v**2, w**2

A = -(1/(bsq*(usq + vsq) +  asq*wsq))
B = bsq*(u*x + v*y) + asq*w*z
C = 0.5*np.sqrt(4*(bsq*(u*x + v*y) + asq*w*z)**2 -
                      4*(bsq*(usq + vsq) + asq*wsq) *
                      (bsq*(-asq + xsq + ysq) + asq*zsq))
t = A * (B + C)
print "t: ", t
xyz, uvw   = np.array([x, y, z]), np.array([u, v, w])
xyzi       = xyz + t*uvw
xi, yi, zi = xyzi
print "point of intersection: ", xyzi
print "check, is it zero within roundoff? ", xi**2/asq + yi**2/asq + zi**2/bsq - 1

es ergibt:

t:  4.33961998769
point of intersection:  [ 0.70399539  0.433962    0.433962  ]
check, is it zero?  -8.54871728961e-15

Die Lösung fällt tatsächlich auf das Ellipsoid.

Wandeln Sie (x, y, z) in Längen- und Breitengrad um

Wir müssen die Gleichungen (in dieser Antwort gezeigt ) umkehren, um den Breiten- und Längengrad auf dem WSG84-Ellipsoid zu erhalten, der diesen Koordinaten entspricht.

Laut [Geographic_coordinate_conversion#From_geodetic_to_ECEF_coordinates][1] von Wikipedia

Die kartesischen 3D-Koordinaten$X, Y, Z$in erdzentrierten, erdfesten Koordinaten unter der Annahme einer Ellipsoidform ist gegeben durch:

$$X = \left(N(\phi) + h \right) \cos\phi \cos\lambda$$

$$Y = \left(N(\phi) + h \right) \cos\phi \sin\lambda$$

$$Z = \left(\frac{b^2}{a^2} N(\phi) + h \right) \sin\phi$$

Wo$\phi, \lambda, h$sind Breite, Länge und Höhe, und$a, b$sind die äquatorialen und polaren Radien des verwendeten Ellipsoids, und

$$N(\phi) = \frac{a^2}{\sqrt{a^2\cos^2\phi + b^2 \sin^2\phi}}.$$

Die Umkehrung wird in Datumstransformationen von GPS-Positionen behandelt; Application Note finden Sie in dieser GIS-Frage sowie in diesem Link in Kommentaren.

Aus der Application Note zeige ich die analytische Lösung .

Definitionen:

$$a = 6378137$$ $$b = a(1-f) = 6356752.31424518$$ $$f = \frac{1}{298.257223563}$$ $$e = \sqrt{\frac{a^2-b^2}{a^2}}$$ $$e' = \sqrt{\frac{a^2-b^2}{b^2}}.$$

Das Auflösen nach Längengrad ist trivial. Die Anwendungsnotiz gibt Folgendes an:

$$\lambda = \arctan\left(\frac{y}{x}\right),$$

In der Praxis können Sie Arctan jedoch nicht blind verwenden, da es nicht alle vier Quadranten unterscheiden kann. Also stattdessen verwenden

$$\lambda = \arctan2(y, x).$$

Analytisch nach Breite auflösen, gibt die Application Note an

$$\phi= \arctan\left( \frac{z + e'^2b \sin^3(\theta)}{p-e^2a\cos^3(\theta)} \right)$$

Wo

$$p = \sqrt{x^2+y^2}$$ $$\theta = \arctan\left(\frac{za}{pb}\right).$$

Wieder sollten Sie wahrscheinlich verwenden$arctan2()$anstatt sich um die Quadranten zu kümmern.

Iterative Lösungen finden Sie auch in dieser Referenz sowie den anderen (bereits oben erwähnten). Iterative Lösungen waren in der Vergangenheit wahrscheinlich deutlich schneller (und könnten es auch jetzt sein), was wichtig wäre, wenn Sie Raytracing für Millionen von Punkten durchführen, z. B. ein Bild der Erdoberfläche aus der Sicht der Satelliten mit einem genauen erzeugen Darstellung.