Las fórmulas de Vincenty

Métodos en geodesia

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 utilizando 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
• Distancia del círculo máximo
• Arco meridiano
• Geodésicas sobre un elipsoide
• Tadeo Vincenty
• Geodesia

Notas

1. σ no se evalúa directamente a partir de sen  σ o cos  σ para preservar la precisión numérica cerca de los polos y el ecuador
2. Si sen σ = 0 el valor de sen α es indeterminado. Representa un punto final coincidente con el punto inicial o diametralmente opuesto a él.
3. Cuando el punto inicial y el final están en el ecuador, C = 0 y no se utiliza el valor de . El valor límite es . ${\displaystyle \cos \left(2\sigma _{\text{m}}\right)}$ ${\displaystyle \cos \left(2\sigma _{\text{m}}\right)=-1}$

Referencias

• Bessel, Friedrich Wilhelm (2010). "El cálculo de longitud y latitud a partir de mediciones geodésicas (1825)". Astron. Nachr . 331 (8): 852–861. arXiv : 0908.1824 . Código Bibliográfico :2010AN....331..852K. doi :10.1002/asna.201011352. S2CID  118760590.Traducción al inglés de Astron. Nachr. 4 , 241–254 (1825).
• Helmert, Friedrich R. (1964). Teorías matemáticas y físicas de la geodesia superior, parte 1 (1880). St. Louis: Aeronautical Chart and Information Center . Consultado el 30 de julio de 2011 . Traducción al inglés de Die Mathematischen und Physikalischen Theorieen der Höheren Geodäsie , vol. 1 (Teubner, Leipzig, 1880).
• Karney, Charles FF (enero de 2013). "Algoritmos para geodésicas". Journal of Geodesy . 87 (1): 43–55. arXiv : 1109.4448 . Código Bibliográfico :2013JGeod..87...43K. doi : 10.1007/s00190-012-0578-z .Adendas.
• Legendre, Adrien-Marie (1806). "Análisis de las huellas de los triángulos sobre la superficie de una esphėroïde". Mémoires de la classe des sciences mathématiques et physiques de l'Institut National de France (primer semestre): 130–161 . Consultado el 30 de julio de 2011 .
• Rainsford, HF (1955). "Geodésicas largas en el elipsoide". Bulletin Géodésique . 37 : 12–22. Bibcode :1955BGeod..29...12R. doi :10.1007/BF02527187. S2CID  122111614.
• Rapp, Ricahrd H. (marzo de 1993). Geodesia geométrica, parte II (informe técnico). Universidad Estatal de Ohio . Consultado el 1 de agosto de 2011 .
• Vincenty, Thaddeus (abril de 1975a). "Soluciones directas e inversas de geodésicas en el elipsoide con aplicación de ecuaciones anidadas" (PDF) . Survey Review . XXIII (176): 88–93. Bibcode :1975SurRv..23...88V. doi :10.1179/sre.1975.23.176.88 . Consultado el 11 de julio de 2009 . Al seleccionar una fórmula para la solución de geodésicas, es de importancia primordial considerar la longitud del programa, es decir, la cantidad de núcleo que ocupará en la computadora junto con las funciones trigonométricas y otras funciones requeridas.
• Vincenty, Thaddeus (agosto de 1975b). Solución inversa geodésica entre puntos antípodas (PDF) (Informe técnico). Escuadrón de reconocimiento geodésico DMAAC. doi :10.5281/zenodo.32999.
• Vincenty, Thaddeus (abril de 1976). "Correspondencia". Survey Review . XXIII (180): 294.
• Manual de referencia del sistema de referencia geocéntrica de Australia (GDA). Comité intergubernamental de topografía y cartografía (ICSM). Febrero de 2006. ISBN 0-9579951-0-5Archivado desde el original (PDF) el 26 de junio de 2009. Consultado el 11 de julio de 2009 .

Enlaces externos

• Calculadoras en línea de Geoscience Australia :
  • Vincenty Direct (punto de destino)
  • Inversa de Vincenty (distancia entre puntos)
• Calculadoras del Servicio Geodético Nacional de Estados Unidos :
  • Utilidades de cálculo ejecutables para PC en línea y descargables, incluyendo problemas directos e inversos, tanto en dos como en tres dimensiones (consultado el 1 de agosto de 2011).
• Calculadoras en línea con código fuente JavaScript de Chris Veness (licencia Creative Commons Attribution):
  • Vincenty Direct (punto de destino)
  • Inversa de Vincenty (distancia entre puntos)
• GeographicLib ofrece una utilidad, GeodSolve (con código fuente con licencia MIT/X11), para resolver problemas geodésicos directos e inversos. En comparación con Vincenty, esta es aproximadamente 1000 veces más precisa (error = 15 nm) y la solución inversa es completa. Aquí se incluye una versión en línea de GeodSolve.
• Implementación completa de fórmulas directas e inversas de Vincenty con código fuente, implementación de Excel VBA por Tomasz Jastrzębski