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 qui suppose que la Terre est sphérique. (distance orthodromique)

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.

Notation

Définition et notation:

Notation	Definition	Valeur
a	longueur du demi-grand axe de l'ellipsoïde (rayon à l'équateur) ;	(6378137.0 mètre WGS-84)
ƒ	aplatissement de l'ellipsoïde;	(1/298.257223563 dans WGS-84)
b = (1 − ƒ) a	longueur du demi-petit axe de l'ellipsoïde (rayon aux pôles);	(6356752.314245 mètre WGS-84)
Φ1, Φ2	latitude des points;
U1 = arctan( (1 − ƒ) tan Φ1 ), U2 = arctan( (1 − ƒ) tan Φ2 )	reduced latitude (latitude sur la sphère auxillière)
L1, L2	longitude des points;
L = L2 − L1	difference en longitude des deux points;
λ	Différence de longitude des points de la sphère auxiliaire ;
α1, α2	forward azimuths at the points;
α	forward azimuth of the geodesic à l'équateur, si il est prolongé 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 les points presque antipodaux. Un exemple de convergence lente est (Φ₁, L₁) = (0°, 0°) and (Φ₂, L₂) = (0.5°, 179.5°) for the WGS84 ellipsoid. This requires about 130 iterations to give a result accurate to 1 mm. Depending on how the inverse method is implemented, the algorithm might return the correct result (19936288.579 m), an incorrect result, or an error indicator. An example of an incorrect result is provided by the NGS online utility, which returns a distance that is about 5 km too long. Vincenty suggested a method of accelerating the convergence in such cases (Rapp, 1993).

An example of a failure of the inverse method to converge is (Φ₁, L₁) = (0°, 0°) and (Φ₂, L₂) = (0.5°, 179.7°) for the WGS84 ellipsoid. In an unpublished report, Vincenty (1975b) gave an alternative iterative scheme to handle such cases. This converges to the correct result 19944127.421 m after about 60 iterations; however, in other cases many thousands of iterations are required.

Karney (2013) reformulated the inverse problem as a one-dimensional root-finding problem; this can be rapidly solved with Newton's method for all pairs of input points.

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

• Calculateur de Geoscience Australia:
  • Vincenty Direct (destination point)
  • Vincenty Inverse (distance between points)
• Calculateurs de the U.S. National Geodetic Survey:
  • Online and downloadable PC-executable calculation utilities, including forward (direct) and inverse problems, in both two and three dimensions (accessed 2011-08-01).
• Calculateur en ligne avec code source JavaScript par Chris Veness (Creative Commons Attribution license):
  • Vincenty Direct (destination point)
  • Vincenty Inverse (distance between points)
• GeographicLib provides a utility GeodSolve (with MIT/X11 licensed source code) for solving direct and inverse geodesic problems. Compared to Vincenty, this is about 1000 times more accurate (error = 15 nm) and the inverse solution is complete. Here is an online version of GeodSolve.
• Complete Vincenty's direct and inverse formulae implementation with source code, Excel VBA implementation by Tomasz Jastrzębski