Formules de Vincenty

Les formules de Vincenty sont deux méthodes itératives apparentées utilisées en géodésie pour calculer la distance entre deux points à la surface d'un sphéroïde, développées par Thaddeus Vincenty (en) en 1975.

Ces formules utilisent l'hypothèse que la figure de la Terre est un sphéroïde aplati aux pôles (en), ce qui permet d'obtenir des résultats plus précis qu'avec la distance du grand cercle (ou distance orthodromique) qui suppose que la Terre est sphérique.

La première méthode (directe) calcule l'emplacement d'un point situé à une distance et un azimut (direction) donnés d'un autre point. La seconde méthode (inverse) calcule la distance géographique et l'azimut entre deux points donnés. Elles ont été largement utilisées en géodésie car elles sont précises à 0,5 mm (0,020 in) sur l'ellipsoïde terrestre.

Historique

L'objectif de Vincenty était d'exprimer les algorithmes existants pour les géodésiques sur un ellipsoïde sous une forme qui minimise la longueur du programme (Vincenty 1975a). Son rapport non publié (1975b) mentionne l'utilisation d'une Wang 720, qui n'avait que quelques kilo-octets de mémoire. Pour obtenir une bonne précision pour les longues lignes, la solution utilise la solution classique de Legendre (1806), Bessel (1825) et Helmert (1880) basée sur la sphère auxiliaire. Vincenty s'est appuyé sur la formulation de cette méthode donnée par Rainsford, 1955. Legendre a montré qu'une géodésique ellipsoïdale peut être exactement cartographiée sur un grand cercle de la sphère auxiliaire en cartographiant la latitude géographique en latitude réduite et en fixant l'azimut du grand cercle à celui de la géodésique. La longitude sur l'ellipsoïde et la distance le long de la géodésique sont alors données en termes de longitude sur la sphère et de longueur d'arc le long du grand cercle par de simples intégrales. Bessel et Helmert ont donné des séries rapidement convergentes pour ces intégrales, ce qui permet de calculer la géodésique avec une précision arbitraire.

Afin de minimiser la taille du programme, Vincenty a pris ces séries, les a réexpandues en utilisant le premier terme de chaque série comme petit paramètre,^[pas clair] et les a tronquées en $O(f^{3})$. Cela a permis d'obtenir des expressions compactes pour les intégrales de longitude et de distance. Les expressions ont été mises sous forme de Horner (ou "imbriquées"), car cela permet d'évaluer les polynômes en n'utilisant qu'un seul registre temporaire. Enfin, des techniques itératives simples ont été utilisées pour résoudre les équations implicites dans les méthodes directe et inverse ; bien qu'elles soient lentes (et dans le cas de la méthode inverse, il arrive qu'elle ne converge pas), elles entraînent la plus faible augmentation de la taille du code.

Définitions et notations

Notation	Définition	Valeur dans WGS 84
a	Longueur du demi-grand axe de l'ellipsoïde (rayon à l'équateur)	6 378 137,0 m
ƒ	Aplatissement de l'ellipsoïde	1/298,257223563
b = (1 − ƒ) a	Longueur du demi-petit axe de l'ellipsoïde (rayon aux pôles)	6 356 752,314 245 m
Φ1, Φ2	Latitude des points
U1 = arctan( (1 − ƒ) tan Φ1 ), U2 = arctan( (1 − ƒ) tan Φ2 )	Latitude réduite (latitude sur la sphère auxiliaire)
L1, L2	Longitude des points
L = L2 − L1	Différence en longitude des deux points
λ	Différence de longitude des points de la sphère auxiliaire
α1, α2	Azimuts de la géodésique aux points
α	Azimut à l'équateur de la géodésique, si elle est prolongée jusque-là
s	Distance ellipsoïdale entre les deux points
σ	Distance angulaire entre les points
σ1	Distance angulaire entre le point et l'équateur
σm	Séparation angulaire entre le point médian de la ligne et l'équateur

Problème inverse

Étant donné les coordonnées des deux points (Φ₁, L₁) et (Φ₂, L₂), le problème inverse permet de trouver les azimuts α₁, α₂ et la distance ellipsoïdale s.

Calculer U₁, U₂ and L, et fixer la valeur initiale de λ = L. Puis évaluer itérativement les équations suivantes jusqu'à ce que λ converge:

$\sin \sigma ={\sqrt {\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 =\operatorname {arctan2} \left(\sin \sigma ,\cos \sigma \right)$

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

$\cos ^{2}\alpha =1-\sin ^{2}\alpha$

$\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\}$

Lorsque λ a convergé vers le degré de précision souhaité (10⁻¹² correspond à environ 0,06mm), évaluez ce qui suit:

${\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}&=\operatorname {arctan2} \left(\cos U_{2}\sin \lambda ,\cos U_{1}\sin U_{2}-\sin U_{1}\cos U_{2}\cos \lambda \right)\\\alpha _{2}&=\operatorname {arctan2} \left(\cos U_{1}\sin \lambda ,-\sin U_{1}\cos U_{2}+\cos U_{1}\sin U_{2}\cos \lambda \right)\end{aligned}}$

Entre deux points presque antipodaux, la formule itérative peut ne pas converger ; cela se produit lorsque la première estimation de λ, telle que calculée par l'équation ci-dessus, est supérieure à π en valeur absolue.

Problème direct

Étant donné un point initial (Φ₁, L₁) et l'azimut initial, α₁, et une distance, s, le long de la géodésique, le problème est de trouver le point d'arrivée (Φ₂, L₂) et l'azimut, α₂.

Commencez par calculer les éléments suivants :

${\begin{aligned}U_{1}&=\arctan \left[(1-f)\tan \phi _{1}\right]\\\sigma _{1}&=\operatorname {arctan2} \left(\tan U_{1},\cos \alpha _{1}\right)\\\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}}}\right)\\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}}$

Ensuite, en utilisant une valeur initiale $\sigma ={\tfrac {s}{bA}}$, itérer les équations suivantes jusqu'à ce qu'il n'y ait pas de changement significatif dans σ :

${\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}}$

Une fois que σ est obtenu avec une précision suffisante, évaluez :

${\begin{aligned}\phi _{2}&=\operatorname {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 &=\operatorname {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 \left[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}&=\operatorname {arctan2} \left(\sin \alpha ,-\sin U_{1}\sin \sigma +\cos U_{1}\cos \sigma \cos \alpha _{1}\right)\end{aligned}}$

Si le point initial est au pôle Nord ou Sud, la première équation est indéterminée. Si l'azimut initial est orienté vers l'est ou l'ouest, la deuxième équation est indéterminée. Si la fonction standard d'arctangente à 2 arguments atan2 est utilisée, ces valeurs sont généralement traitées correctement.^[pas clair]

Modification de Vincenty

Dans sa lettre à Survey Review en 1976, Vincenty a suggéré de remplacer ses expressions de séries pour A et B par des formules plus simples utilisant le paramètre d'expansion k₁ de Helmert :

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

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

où

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

Points presque antipodaux

Comme indiqué ci-dessus, la solution itérative du problème inverse ne converge pas, ou converge lentement, pour des points presque antipodaux. Un exemple de convergence lente est (Φ1, L1) = (0°, 0°) et (Φ2, L2) = (0,5°, 179,5°) pour l’ellipsoïde WGS 84. Cela nécessite environ 130 itérations pour obtenir un résultat précis à 1 mm près. Selon la manière dont la méthode inverse est implémentée, l’algorithme peut renvoyer le résultat correct (19 936 288,579 m), un résultat incorrect ou un indicateur d’erreur. Un exemple de résultat incorrect est donné par l’utilitaire en ligne de la NGS, qui renvoie une distance trop longue d’environ 5 km. Vincenty a proposé une méthode pour accélérer la convergence dans de tels cas (Rapp, 1993).

Un exemple d’échec de convergence de la méthode inverse est (Φ1, L1) = (0°, 0°) et (Φ2, L2) = (0,5°, 179,7°) pour l’ellipsoïde WGS 84. Dans un rapport non publié, Vincenty (1975b) a présenté un schéma itératif alternatif pour traiter de tels cas. Celui-ci converge vers le résultat correct (19 944 127,421 m) après environ 60 itérations ; toutefois, dans d’autres cas, plusieurs milliers d’itérations sont nécessaires.

Karney (2013) a reformulé le problème inverse comme un problème de recherche de racine unidimensionnel ; celui-ci peut être résolu rapidement à l’aide de la méthode de Newton pour toutes les paires de points en entrée.

Voir aussi

• Modèle ellipsoïdal de la Terre

Références

(en) Cet article est partiellement ou en totalité issu de l’article de Wikipédia en anglais intitulé « Vincenty's formulae » (voir la liste des auteurs).

Lien externe

• Calculateurs de Geoscience Australia :
  • solution directe (emplacement d'un point en fonction d'un autre point, d'une distance et d'un azimut) ;
  • solution inverse (distance et azimut entre deux points).
• Calculateurs du U.S. National Geodetic Survey :
  • outils de calcul en ligne et téléchargeable des solutions directes et inverses, en 2 et 3 dimensions (consulté le 1^er août 2011).
• Calculateurs en ligne avec code source JavaScript par Chris Veness (licence Creative Commons Attribution) :
  • solution directe ;
  • solution inverse.
• GeographicLib fournit l'outil GeodSolve permettant de calculer les solutions directe et inverse. Comparé aux formules de Vincenty, il est environ 1000 fois plus précis (avec une erreur de 15 nm), et la solution inverse est complète. Une version en ligne est aussi disponible.
• Implémentation complète des formules de Vincenty en Excel VBA, par Tomasz Jastrzębski