Las fórmulas de Vincenty

Methods in geodesy

Las fórmulas de Vincenty son dos métodos iterativos relacionados que se utilizan en geodesia para calcular la distancia entre dos puntos de la superficie de un esferoide , desarrollados por Thaddeus Vincenty (1975a). Se basan en el supuesto de que la figura de la Tierra es un esferoide achatado y, por lo tanto, son más precisos que los métodos que suponen una Tierra esférica , como la distancia del círculo máximo .

El primer método (directo) calcula la ubicación de un punto que se encuentra a una distancia y un acimut (dirección) determinados de otro punto. El segundo método (inverso) calcula la distancia geográfica y el acimut entre dos puntos determinados. Se han utilizado ampliamente en geodesia porque tienen una precisión de 0,5 mm (0,020  pulgadas) en el elipsoide terrestre .

Fondo

El objetivo de Vincenty era expresar los algoritmos existentes para las geodésicas en un elipsoide de una forma que minimizara la longitud del programa (Vincenty 1975a). Su informe inédito (1975b) menciona el uso de una calculadora de escritorio Wang 720, que tenía sólo unos pocos kilobytes de memoria. Para obtener una buena precisión para líneas largas, la solución utiliza la solución clásica de Legendre (1806), Bessel (1825) y Helmert (1880) basada en la esfera auxiliar. Vincenty se basó en la formulación de este método dada por Rainsford, 1955. Legendre demostró que una geodésica elipsoidal puede ser mapeada exactamente a un círculo máximo en la esfera auxiliar mapeando la latitud geográfica a la latitud reducida y estableciendo el acimut del círculo máximo igual al de la geodésica. La longitud en el elipsoide y la distancia a lo largo de la geodésica se expresan entonces en términos de la longitud en la esfera y la longitud del arco a lo largo del círculo máximo mediante integrales simples. Bessel y Helmert proporcionaron series de convergencia rápida para estas integrales, que permiten calcular la geodésica con una precisión arbitraria.

Para minimizar el tamaño del programa, Vincenty tomó estas series, las volvió a expandir usando el primer término de cada serie como el parámetro pequeño, ^{[ aclaración necesaria ]} y las truncó a . Esto dio como resultado expresiones compactas para las integrales de longitud y distancia. Las expresiones se pusieron en forma de Horner (o anidada ), ya que esto permite que los polinomios se evalúen usando solo un único registro temporal. Finalmente, se utilizaron técnicas iterativas simples para resolver las ecuaciones implícitas en los métodos directo e inverso; aunque estas son lentas (y en el caso del método inverso a veces no converge), dan como resultado el menor aumento en el tamaño del código. ${\displaystyle O(f^{3})}$

Notación

Defina la siguiente notación:

Problema inverso

Dadas las coordenadas de los dos puntos ( Φ ₁ ,  L ₁ ) y ( Φ ₂ ,  L ₂ ), el problema inverso encuentra los acimutes α ₁ , α ₂ y la distancia elipsoidal s .

Calcular U ₁ , U ₂ y L , y establecer el valor inicial de λ = L . Luego evaluar iterativamente las siguientes ecuaciones hasta que λ converja:

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

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

${\displaystyle \sigma =\operatorname {arctan2} \left(\sin \sigma ,\cos \sigma \right)}$^[1]

${\displaystyle \sin \alpha ={\frac {\cos U_{1}\cos U_{2}\sin \lambda }{\sin \sigma }}}$^[2]

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

${\displaystyle \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 }}}$^[3]

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

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

Cuando λ haya convergido al grado de precisión deseado (10 ⁻¹² corresponde aproximadamente a 0,06  mm), evalúe lo siguiente:

${\displaystyle {\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 {1}{6}}B\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 dos puntos casi antípodas, la fórmula iterativa puede no converger; esto ocurrirá cuando la primera estimación de λ calculada mediante la ecuación anterior sea mayor que π en valor absoluto.

Problema directo

Dado un punto inicial ( Φ ₁ , L ₁ ) y un acimut inicial, α ₁ , y una distancia, s , a lo largo de la geodésica, el problema es encontrar el punto final ( Φ ₂ , L ₂ ) y el acimut, α ₂ .

Comience calculando lo siguiente:

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

Luego, utilizando un valor inicial , itere las siguientes ecuaciones hasta que no haya ningún cambio significativo en σ : ${\displaystyle \sigma ={\tfrac {s}{bA}}}$

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

Una vez obtenido σ con suficiente precisión se evalúa:

${\displaystyle {\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 el punto inicial está en el polo norte o sur, entonces la primera ecuación es indeterminada. Si el acimut inicial es hacia el este o el oeste, entonces la segunda ecuación es indeterminada. Si se utiliza la función arcotangente estándar de 2 argumentos atan2 , estos valores se manejan normalmente correctamente. ^{[ aclaración necesaria ]}

Modificación de Vincenty

En su carta a Survey Review en 1976, Vincenty sugirió reemplazar sus expresiones de series para A y B con fórmulas más simples utilizando el parámetro de expansión de Helmert k ₁ :

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

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

dónde

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

Puntos casi antípodas

Como se señaló anteriormente, la solución iterativa al problema inverso no converge o converge lentamente para puntos casi antípodas. Un ejemplo de convergencia lenta es ( Φ ₁ ,  L ₁ ) = (0°, 0°) y ( Φ ₂ ,  L ₂ ) = (0,5°, 179,5°) para el elipsoide WGS84. Esto requiere alrededor de 130 iteraciones para dar un resultado preciso a 1 mm. Dependiendo de cómo se implemente el método inverso, el algoritmo puede devolver el resultado correcto (19936288,579 m), un resultado incorrecto o un indicador de error. Un ejemplo de un resultado incorrecto lo proporciona la utilidad en línea NGS, que devuelve una distancia que es aproximadamente 5 km demasiado larga. Vincenty sugirió un método para acelerar la convergencia en tales casos (Rapp, 1993).

Un ejemplo de un fallo del método inverso para converger es ( Φ ₁ ,  L ₁ ) = (0°, 0°) y ( Φ ₂ ,  L ₂ ) = (0,5°, 179,7°) para el elipsoide WGS84. En un informe no publicado, Vincenty (1975b) dio un esquema iterativo alternativo para manejar tales casos. Este converge al resultado correcto 19944127,421 m después de aproximadamente 60 iteraciones; sin embargo, en otros casos se requieren muchos miles de iteraciones.

Karney (2013) reformuló el problema inverso como un problema de búsqueda de raíces unidimensional ; esto puede resolverse rápidamente con el método de Newton para todos los pares de puntos de entrada.

Véase también

• Distancia geográfica
• Great-circle distance
• Meridian arc
• Geodesics on an ellipsoid
• Thaddeus Vincenty
• Geodesy

Notes

1. σ is not evaluated directly from sin σ or cos σ to preserve numerical accuracy near the poles and equator
2. If sin σ = 0 the value of sin α is indeterminate. It represents an end point coincident with, or diametrically opposed to, the start point.
3. Where the start and end point are on the equator, C = 0 and the value of ${\displaystyle \cos \left(2\sigma _{\text{m}}\right)}$ is not used. The limiting value is ${\displaystyle \cos \left(2\sigma _{\text{m}}\right)=-1}$.

References

• Bessel, Friedrich Wilhelm (2010). "The calculation of longitude and latitude from geodesic measurements (1825)". Astron. Nachr. 331 (8): 852–861. arXiv:0908.1824. Bibcode:2010AN....331..852K. doi:10.1002/asna.201011352. S2CID 118760590. English translation of Astron. Nachr. 4, 241–254 (1825).
• Helmert, Friedrich R. (1964). Mathematical and Physical Theories of Higher Geodesy, Part 1 (1880). St. Louis: Aeronautical Chart and Information Center. Retrieved 2011-07-30. English translation of Die Mathematischen und Physikalischen Theorieen der Höheren Geodäsie, Vol. 1 (Teubner, Leipzig, 1880).
• Karney, Charles F. F. (January 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. Addenda.
• 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 (1st sem): 130–161. Retrieved 2011-07-30.
• Rainsford, H. F. (1955). "Long geodesics on the ellipsoid". Bulletin Géodésique. 37: 12–22. Bibcode:1955BGeod..29...12R. doi:10.1007/BF02527187. S2CID 122111614.
• Rapp, Ricahrd H. (March 1993). Geometric Geodesy, Part II (Technical report). Ohio State University. Retrieved 2011-08-01.
• Vincenty, Thaddeus (April 1975a). "Direct and Inverse Solutions of Geodesics on the Ellipsoid with application of nested equations" (PDF). Survey Review. XXIII (176): 88–93. Bibcode:1975SurRv..23...88V. doi:10.1179/sre.1975.23.176.88. Retrieved 2009-07-11. In selecting a formula for the solution of geodesics it is of primary importance to consider the length of the program, that is the amount of core which it will occupy in the computer along with trigonometric and other required functions.
• Vincenty, Thaddeus (August 1975b). Geodetic inverse solution between antipodal points (PDF) (Technical report). DMAAC Geodetic Survey Squadron. doi:10.5281/zenodo.32999.
• Vincenty, Thaddeus (April 1976). "Correspondence". Survey Review. XXIII (180): 294.
• Geocentric Datum of Australia (GDA) Reference Manual. Intergovernmental committee on survey and mapping (ICSM). February 2006. ISBN 0-9579951-0-5. Archived from the original (PDF) on 2009-06-26. Retrieved 2009-07-11.

External links

• Online calculators from Geoscience Australia:
  • Vincenty Direct (destination point)
  • Vincenty Inverse (distance between points)
• Calculators from 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).
• Online calculators with JavaScript source code by 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