Vincentys formler

För den matematiska bakgrunden till detta problem, se Geodesik på en ellipsoid .

Vincentys formler är två relaterade iterativa metoder som används inom geodesi för att beräkna avståndet mellan två punkter på ytan av en sfäroid, utvecklade av Thaddeus Vincenty ( 1975a). De är baserade på antagandet att jordens figur är en oblate sfäroid , och därför är mer exakta än metoder som antar en sfärisk jord, såsom storcirkelavstånd .

  Den första (direkta) metoden beräknar platsen för en punkt som är ett givet avstånd och azimut (riktning) från en annan punkt. Den andra (omvända) metoden beräknar det geografiska avståndet och azimuten mellan två givna punkter. De har använts i stor utsträckning inom geodesi eftersom de är exakta till inom 0,5 mm (0,020 tum) på jordens ellipsoid .

Bakgrund

Vincentys mål var att uttrycka befintliga algoritmer för geodetik på en ellipsoid i en form som minimerar programlängden (Vincenty 1975a). Hans opublicerade rapport (1975b) nämner användningen av en Wang 720 skrivbordsräknare, som bara hade några få kilobyte minne. För att få bra noggrannhet för långa linor använder lösningen den klassiska lösningen av Legendre (1806), Bessel (1825) och Helmert (1880) baserad på hjälpsfären. Vincenty förlitade sig på formuleringen av denna metod som gavs av Rainsford, 1955. Legendre visade att en ellipsoidal geodetik exakt kan mappas till en storcirkel på hjälpsfären genom att kartlägga den geografiska latituden till reducerad latitud och sätta storcirkelns azimut lika med den av det geodetiska. Longituden på ellipsoiden och avståndet längs geodetiken ges sedan i termer av longituden på sfären och båglängden längs storcirkeln av enkla integraler. Bessel och Helmert gav snabbt konvergerande serier för dessa integraler, vilket gör att geodetiken kan beräknas med godtycklig noggrannhet.

För att minimera programstorleken tog Vincenty dessa serier, expanderade dem på nytt med den första termen i varje serie som den lilla parametern, [ och trunkerade dem till $O(f^{ 3})$ . Detta resulterade i kompakta uttryck för longitud och distansintegraler. Uttrycken sattes i Horner (eller kapslade ) form, eftersom detta gör att polynom kan utvärderas med endast ett enda temporärt register. Slutligen användes enkla iterativa tekniker för att lösa de implicita ekvationerna i de direkta och inversa metoderna; även om dessa är långsamma (och i fallet med den inversa metoden konvergerar den ibland inte), resulterar de i den minsta ökningen av kodstorlek.

Notation

Definiera följande notation:

Notation	Definition	Värde
a	längden på ellipsoidens halvstora axel (radie vid ekvatorn);	(6378137,0 meter i WGS-84 )
ƒ	tillplattning av ellipsoiden;	(1/298.257223563 i WGS-84 )
b = (1 − ƒ ) a	längden på ellipsoidens halva axel (radie vid polerna);	(6356752.314245 meter i WGS-84 )
Φ 1 , Φ 2	latitud för punkterna;
U 1 = arctan( (1 − ƒ ) tan Φ 1 ), U 2 = arctan( (1 − ƒ ) tan Φ 2 )	reducerad latitud (latitud på hjälpsfären)
L 1 , L 2	punkternas longitud ;
L = L 2 − L 1	skillnad i longitud för två punkter;
λ	Skillnad i longitud för punkterna på hjälpsfären;
α 1 , α 2	framåtriktade azimuter vid punkterna;
α	framåtriktad azimut för geodetiken vid ekvatorn, om den sträcktes så långt;
s	ellipsoidalt avstånd mellan de två punkterna;
σ	vinkelavstånd mellan punkter
σ 1	vinkelseparation mellan punkten och ekvatorn
σ m	vinkelseparation mellan linjens mittpunkt och ekvatorn

Omvänt problem

Givet koordinaterna för de två punkterna ( Φ ₁ , L ₁ ) och ( Φ ₂ , L ₂ ), finner det omvända problemet azimuterna α ₁ , α ₂ och det ellipsoidala avståndet s .

Beräkna U ₁ , U ₂ och L , och ställ in initialvärdet på λ = L . Utvärdera sedan iterativt följande ekvationer tills λ konvergerar:

$\sin \sigma = {\sqr \left(\cos U_{2}\sin \lambda \right)^{2}+\left(\cos U_{1}\sin U_{2}-\sin U_{1}\cos U_{2}\ cos \lambda \right)^{2}}}$

$\cos \sigma =\sin U_{1 }\sin U_{2}+\cos U_{1}\cos U_{2}\cos \lambda \,$

$\sigma =\operatörsnamn {arctan2 } \left(\sin \sigma ,\cos \sigma \right)$

$\sin \alpha ={\frac {\cos U_ {1}\cos U_{2}\sin \lambda }{\sin \sigma }}$

$\cos \left(2\sigma _{\text{m}}\right)=\cos \sigma -{\frac {2 \sin U_{1}\sin U_{2}}{\cos ^{2}\alpha }}=\cos \sigma -{\frac {2\sin U_{1}\sin U_{2}}{1 -\sin ^{2}\alpha }}$

$C={\frac {f}{16}}\ cos ^{2}\alpha \left[4+f\left(4-3\cos ^{2}\alpha \right)\right]$

$\lambda =L+(1-C)f\sin \alpha \left\{\sigma +C\sin \sigma \left[\cos \left(2\sigma _{\text{m}}\right)+C\cos \sigma \left(-1+2\cos ^ {2}\left(2\sigma _{\text{m}}\right)\right)\right]\right\}$

  När λ har konvergerat till önskad noggrannhetsgrad (10 ⁻¹² motsvarar ungefär 0,06 mm), utvärdera följande:

${\begin{aligned}u^{2}&=\cos ^{2}\alpha \left({ \frac {a^{2}-b^{2}}{b^{2}}}\right)\\A&=1+{\frac {u^{2}}{16384}}\left(4096 +u^{2}\left[-768+u^{2}\left(320-175u^{2}\right)\right]\right)\\B&={\frac {u^{2}} {1024}}\left(256+u^{2}\left[-128+u^{2}\left(74-47u^{2}\right)\right]\right)\\\Delta \sigma &=B\sin \sigma \left\{\cos(2\sigma _{\text{m}})+{\frac {1}{4}}B\left(\cos \sigma \left[-1 +2\cos ^{2}\left(2\sigma _{\text{m}}\right)\right]-{\frac {B}{6}}\cos \left[2\sigma _{\ text{m}}\right]\left[-3+4\sin ^{2}\sigma \right]\left[-3+4\cos ^{2}\left(2\sigma _{\text{ m}}\right)\right]\right)\right\}\\s&=bA(\sigma -\Delta \sigma )\,\\\alpha _{1}&=\operatörsnamn {arctan2} \left( \cos U_{2}\sin \lambda ,\cos U_{1}\sin U_{2}-\sin U_{1}\cos U_{2}\cos \lambda \right)\\\alpha _{2 }&=\operatörsnamn {arctan2} \left(\cos U_{1}\sin \lambda ,-\sin U_{1}\cos U_{2}+\cos U_{1}\sin U_{2}\cos \lambda \right)\end{aligned}}$

Mellan två nästan antipodala punkter kan den iterativa formeln misslyckas med att konvergera; detta inträffar när den första gissningen på λ , beräknad av ekvationen ovan, är större än π i absolut värde.

Direkt problem

Givet en initialpunkt ( Φ ₁ , L ₁ ) och initial azimut, α ₁ , och ett avstånd, s , längs geodetiken är problemet att hitta slutpunkten ( Φ ₂ , L ₂ ) och azimut, α ₂ .

Börja med att räkna ut följande:

${\begin{aligned}U_{1}&=\arctan \left[ (1-f)\tan \phi _{1}\right]\\\sigma _{1}&=\operatörsnamn {arctan2} \left(\tan U_{1},\cos \alpha _{1}\ höger)\\\sin \alpha &=\cos U_{1}\sin \alpha _{1}\\u^{2}&=\cos ^{2}\alpha \left({\frac {a^ {2}-b^{2}}{b^{2}}}\right)=\left(1-\sin ^{2}\alpha \right)\left({\frac {a^{2} -b^{2}}{b^{2}}}\höger)\\A&=1+{\frac {u^{2}}{16384}}\left(4096+u^{2}\left [-768+u^{2}(320-175u^{2})\right]\right)\\B&={\frac {u^{2}}{1024}}\left(256+u^{ 2}\left[-128+u^{2}\left(74-47u^{2}\right)\right]\right)\end{aligned}}$

Använd sedan ett initialt värde ${\displaystyle \sigma ={\tfrac {s}{bA}}} ,$ iterera följande ekvationer tills det inte finns någon signifikant förändring i σ :

$displaystyle {\begin{aligned}2 \sigma _{\text{m}}&=2\sigma _{1}+\sigma \\\Delta \sigma &=B\sin \sigma \left\{\cos \left(2\sigma _{\ text{m}}\right)+{\frac {1}{4}}B\left(\cos \sigma \left[-1+2\cos ^{2}\left(2\sigma _{\text {m}}\right)\right]-{\frac {B}{6}}\cos \left[2\sigma _{\text{m}}\right]\left[-3+4\sin ^ {2}\sigma \right]\left[-3+4\cos ^{2}\left(2\sigma _{\text{m}}\right)\right]\right)\right\}\\ \sigma &={\frac {s}{bA}}+\Delta \sigma \end{aligned}}}$

När σ har erhållits med tillräcklig noggrannhet, utvärdera:

${\begin{aligned}\phi _{2}&=\operatörsnamn {arctan2} \left(\sin U_{1}\cos \sigma +\cos U_{1}\sin \sigma \cos \ alpha _{1},(1-f){\sqrt {\sin ^{2}\alpha +\left(\sin U_{1}\sin \sigma -\cos U_{1}\cos \sigma \cos \alpha _{1}\right)^{2}}}\right)\\\lambda &=\operatörsnamn {arctan2} \left(\sin \sigma \sin \alpha _{1},\cos U_{1 }\cos \sigma -\sin U_{1}\sin \sigma \cos \alpha _{1}\right)\\C&={\frac {f}{16}}\cos ^{2}\alpha \ vänster[4+f\left(4-3\cos ^{2}\alpha \right)\right]\\L&=\lambda -(1-C)f\sin \alpha \left\{\sigma +C \sin \sigma \left(\cos \left[2\sigma _{\text{m}}\right]+C\cos \sigma \left[-1+2\cos ^{2}\left(2\ sigma _{\text{m}}\right)\right]\right)\right\}\\L_{2}&=L+L_{1}\\\alpha _{2}&=\operatörsnamn {arctan2 } \left(\sin \alpha ,-\sin U_{1}\sin \sigma +\cos U_{1}\cos \sigma \cos \alpha _{1}\right)\end{aligned}}$

Om den initiala punkten är vid nord- eller sydpolen är den första ekvationen obestämd. Om den initiala azimuten är rakt åt öster eller väster är den andra ekvationen obestämd. Om en atan2 med dubbelt värde används, hanteras dessa värden vanligtvis korrekt.

Vincentys modifikation

I sitt brev till Survey Review 1976 föreslog Vincenty att ersätta sina serieuttryck för A och B med enklare formler med hjälp av Helmerts expansionsparameter k ₁ :

$A={\frac {1+{\frac {1}{4}}{k_{1}}^{2}}{1-k_{ 1}}}$

$B=k_{1}\left(1-{\frac {3}{8}}{k_{1}}^{ 2}\right)$

var

$k_{1}={\frac {{\sqrt {1+u^{2}}}-1}{{\sqrt {1+ u^{2}}}+1}}$

Nästan antipoda punkter

Som noterats ovan misslyckas den iterativa lösningen av det omvända problemet att konvergera eller konvergerar långsamt för nästan antipoda punkter. Ett exempel på långsam konvergens är ( _Φ1 , _L1 ) = (0°, 0°) och ( _Φ2 , L2 ) = (0,5°, _179,5 °) för WGS84-ellipsoiden . Detta kräver cirka 130 iterationer för att ge ett resultat exakt 1 mm. Beroende på hur den omvända metoden implementeras kan algoritmen returnera det korrekta resultatet (19936288.579 m), ett felaktigt resultat eller en felindikator. Ett exempel på ett felaktigt resultat tillhandahålls av NGS onlineverktyg , som returnerar en sträcka som är cirka 5 km för lång. Vincenty föreslog en metod för att påskynda konvergensen i sådana fall (Rapp, 1993).

Ett exempel på ett misslyckande med den inversa metoden att konvergera är ( Φ ₁ , L ₁ ) = (0°, 0°) och ( Φ ₂ , L ₂ ) = (0,5°, 179,7°) för WGS84-ellipsoiden. I en opublicerad rapport gav Vincenty (1975b) ett alternativt iterativt schema för att hantera sådana fall. Detta konvergerar till det korrekta resultatet 19944127.421 m efter cirka 60 iterationer; i andra fall krävs dock många tusen iterationer.

Karney (2013) omformulerade det omvända problemet som ett endimensionellt problem för att hitta roten ; detta kan snabbt lösas med Newtons metod för alla par av ingångspunkter.

Se även

• Geografiskt avstånd
• Storcirkelavstånd
• Meridianbåge
• Geodesik på en ellipsoid
• Thaddeus Vincenty
• Geodesi

Anteckningar

•   Bessel, Friedrich Wilhelm (2010). "Beräkningen av longitud och latitud från geodetiska mätningar (1825)". Astron. Nachr . 331 (8): 852–861. arXiv : 0908.1824 . Bibcode : 2010AN....331..852K . doi : 10.1002/asna.201011352 . S2CID 118760590 . Engelsk översättning av Astron. Nachr. 4 , 241–254 (1825).
• Helmert, Friedrich R. (1964). Matematiska och fysiska teorier om högre geodesi, del 1 (1880) . St Louis: Aeronautical Chart and Information Center . Hämtad 2011-07-30 . Engelsk översättning av Die Mathematischen und Physikalischen Theorieen der Höheren Geodäsie , Vol. 1 (Teubner, Leipzig, 1880).
• Karney, Charles FF (januari 2013). "Algorithms for geodesics" . Journal of Geodesy . 87 (1): 43–55. arXiv : 1109.4448 . Bibcode : 2013JGeod..87...43K . doi : 10.1007/s00190-012-0578-z . Tillägg .
• Legendre, Adrien-Marie (1806). "Analyse des triangles tracės sur la surface d'un sphėroïde" . Mémoires de la classe des sciences mathématiques et physiques de l'Institut National de France (1:a sem): 130–161 . Hämtad 2011-07-30 .
•   Rainsford, HF (1955). "Lång geodetik på ellipsoiden". Bulletin Geodésique . 37 : 12–22. Bibcode : 1955BGeod..29...12R . doi : 10.1007/BF02527187 . S2CID 122111614 .
• Rapp, Ricahrd H. (mars 1993). Geometrisk geodesi, del II (Teknisk rapport). Ohio State University . Hämtad 2011-08-01 .
• Vincenty, Thaddeus (april 1975a). "Direkta och omvända lösningar av geodesik på ellipsoiden med tillämpning av kapslade ekvationer" ( PDF) . Enkätgranskning . XXIII (176): 88–93. doi : 10.1179/sre.1975.23.176.88 . Hämtad 2009-07-11 . När man väljer en formel för lösningen av geodetik är det av största vikt att beakta längden på programmet, det vill säga mängden kärna som det kommer att uppta i datorn tillsammans med trigonometriska och andra nödvändiga funktioner.
• Vincenty, Thaddeus (augusti 1975b). Geodetisk omvänd lösning mellan antipodalpunkter (PDF) (Teknisk rapport). DMAAC Geodetic Survey Squadron. doi : 10.5281/zenodo.32999 .
• Vincenty, Thaddeus (april 1976). "Korrespondens". Enkätgranskning . XXIII (180): 294.
•   Geocentric Datum of Australia (GDA) Referensmanual (PDF) . Mellanstatlig kommitté för kartläggning och kartläggning (ICSM). Februari 2006. ISBN 0-9579951-0-5 . Hämtad 2009-07-11 .

externa länkar

• Onlineräknare från Geoscience Australia :
  • Vincenty Direct (destination)
  • Vincenty Inverse (avstånd mellan punkter)
• Miniräknare från US National Geodetic Survey :
  • Online- och nedladdningsbara PC-körbara beräkningsverktyg , inklusive framåt (direkt) och inversa problem, i både två och tre dimensioner (tillgänglig 2011-08-01).
• Onlineräknare med JavaScript-källkod av Chris Veness (Creative Commons Attribution-licens):
  • Vincenty Direct (destination)
  • Vincenty Inverse (avstånd mellan punkter)
• GeographicLib tillhandahåller ett verktyg GeodSolve (med MIT/X11-licensierad källkod) för att lösa direkta och omvända geodetiska problem. Jämfört med Vincenty är detta cirka 1000 gånger mer exakt (fel = 15 nm) och den omvända lösningen är komplett. Här är en onlineversion av GeodSolve .
• Slutför Vincentys direkta och omvända formlerimplementering med källkod, Excel VBA-implementering av Tomasz Jastrzębski