Формулы Винсенти

Математические основы этой проблемы см. в разделе «Геодезика на эллипсоиде» .

Формулы Винсенти — это два связанных итеративных метода, используемых в геодезии для расчета расстояния между двумя точками на поверхности сфероида , разработанных Таддеусом Винсенти (1975a). Они основаны на предположении, что форма Земли представляет собой сплюснутый сфероид, и, следовательно, более точны, чем методы, предполагающие сферическую форму Земли, такие как расстояние по большому кругу .

Первый (прямой) метод вычисляет местоположение точки, которая находится на заданном расстоянии и азимуте (направлении) от другой точки. Второй (обратный) метод вычисляет географическое расстояние и азимут между двумя заданными точками. Они широко используются в геодезии, поскольку их точность   составляет не более 0,5 мм (0,020 дюйма) на земном эллипсоиде .

Фон

Целью Винсенти было выразить существующие алгоритмы геодезических на эллипсоиде в форме, минимизирующей длину программы (Винсенти 1975а). В его неопубликованном отчете (1975b) упоминается использование настольного калькулятора Wang 720, у которого было всего несколько килобайт памяти. Чтобы получить хорошую точность для длинных линий, в решении используется классическое решение Лежандра (1806 г.), Бесселя (1825 г.) и Гельмерта (1880 г.), основанное на вспомогательной сфере. Винсенти полагался на формулировку этого метода, данную Рейнсфордом в 1955 году. Лежандр показал, что эллипсоидальную геодезическую можно точно отобразить на большом круге на вспомогательной сфере, сопоставив географическую широту с уменьшенной широтой и установив азимут большого круга равным этому. геодезического. Долгота на эллипсоиде и расстояние вдоль геодезической затем выражаются через долготу на сфере и длину дуги вдоль большого круга с помощью простых интегралов. Бессель и Гельмерт дали для этих интегралов быстро сходящиеся ряды, которые позволяют вычислять геодезическую с произвольной точностью.

Чтобы минимизировать размер программы, Винсенти взял эти ряды и повторно расширил их, используя первый член каждого ряда в качестве маленького параметра: и сократил их до $O(f^{3})$. В результате были получены компактные выражения для интегралов долготы и расстояния. Выражения были помещены в форму Хорнера (или вложенную ), поскольку это позволяет оценивать полиномы, используя только один временный регистр. Наконец, для решения неявных уравнений прямым и обратным методами были использованы простые итерационные методы; хотя они и медленны (а в случае обратного метода он иногда не сходится), они приводят к наименьшему увеличению размера кода.

Обозначения

Определите следующие обозначения:

Обозначения	Определение	Ценить
а	длина большой полуоси эллипсоида (радиус на экваторе);	(6378137,0 метров в WGS-84 )
ƒ	уплощение эллипсоида;	(1/298,257223563 в WGS-84 )
б знак равно (1 - ƒ ) а	длина малой полуоси эллипсоида (радиус в полюсах);	(6356752,314245 метров в WGS-84 )
Ф 1 , Ф 2	широта точек;
U 1 = arctan( (1 - ƒ ) tan Φ 1 ), U 2 = arctan( (1 - ƒ ) tan Φ 2 )	уменьшенная широта (широта на вспомогательной сфере)
Л1 , LЛ2	долгота точек;
Л = Л 2 - Л 1	разница в долготе двух точек;
л	Разность долгот точек на вспомогательной сфере;
а 1 , а 2	вперед азимуты в точках;
а	прямой азимут геодезической на экваторе, если бы она была продлена так далеко;
с	эллипсоидное расстояние между двумя точками;
п	угловое расстояние между точками
п 1	угловое расстояние между точкой и экватором
σ м	угловое расстояние между серединой линии и экватором

Обратная задача

Учитывая координаты двух точек ( Φ ₁ , L ₁ ) и ( Φ ₂ , L ₂ ), обратная задача находит азимуты α ₁ , α ₂ и эллипсоидальное расстояние s .

Рассчитайте U ₁ , U ₂ и L и установите начальное значение λ = L . Затем итеративно вычисляйте следующие уравнения, пока λ не сходится:

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

Когда λ сошлась до желаемой степени точности (10 ⁻¹² соответствует примерно 0,06   мм), оцените следующее:

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

Между двумя почти противоположными точками итерационная формула может не сойтись; это произойдет, когда первое предположение о λ , вычисленное по приведенному выше уравнению, больше π по абсолютному значению.

Прямая проблема

Учитывая начальную точку ( Φ ₁ , L ₁ ) и начальный азимут α ₁ , а также расстояние s вдоль геодезической, задача состоит в том, чтобы найти конечную точку ( Φ ₂ , L ₂ ) и азимут α ₂ .

Начните с расчета следующего:

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

Затем, используя начальное значение $\sigma ={\tfrac {s}{bA}}$, повторяйте следующие уравнения до тех пор, пока не прекратится существенное изменение σ :

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

Как только σ будет получено с достаточной точностью, оцените:

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

Если начальная точка находится на Северном или Южном полюсе, то первое уравнение является неопределенным. Если начальный азимут направлен на восток или запад, то второе уравнение является неопределенным. стандартная функция арктангенса atan2 с двумя аргументами, то эти значения обычно обрабатываются правильно. Если используется

модификация Винсенти

В своем письме в Survey Review в 1976 году Винсенти предложил заменить свои выражения рядов для A и B более простыми формулами, использующими параметр разложения Гельмерта 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)$

где

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

Почти противоположные точки

Как отмечалось выше, итерационное решение обратной задачи не сходится или сходится медленно для почти противоположных точек. Примером медленной сходимости является ( Φ ₁ , L ₁ ) = (0°, 0°) и ( Φ ₂ , L ₂ ) = (0,5°, 179,5°) для эллипсоида WGS84. Для получения результата с точностью до 1 мм требуется около 130 итераций. В зависимости от того, как реализован обратный метод, алгоритм может возвращать правильный результат (19936288,579 м), неправильный результат или индикатор ошибки. Пример неверного результата дает онлайн-утилита NGS , которая возвращает расстояние, превышающее примерно 5 км. Винсенти предложил метод ускорения сходимости в таких случаях (Rapp, 1993).

Примером неудачи обратного метода является ( Φ ₁ , L ₁ ) = (0°, 0°) и ( Φ ₂ , L ₂ ) = (0,5°, 179,7°) для эллипсоида WGS84. В неопубликованном отчете Винсенти (1975b) предложил альтернативную итерационную схему для решения таких случаев. Это сходится к правильному результату 19944127,421 м примерно после 60 итераций; однако в других случаях требуются многие тысячи итераций.

Карни (2013) переформулировал обратную задачу как одномерную задачу поиска корня ; это можно быстро решить с помощью метода Ньютона для всех пар входных точек.

См. также

• Географическое расстояние
• Расстояние большого круга
• Меридианная дуга
• Геодезические на эллипсоиде
• Таддеус Винсенти
• Геодезия

Ссылки

• Бессель, Фридрих Вильгельм (2010). «Расчет долготы и широты по геодезическим измерениям (1825 г.)». Астрон. Нахр . 331 (8): 852–861. arXiv : 0908.1824 . Бибкод : 2010AN....331..852K . дои : 10.1002/asna.201011352 . S2CID   118760590 . Английский перевод Астрона. Нахр. 4 , 241–254 (1825).
• Хельмерт, Фридрих Р. (1964). Математические и физические теории высшей геодезии, Часть 1 (1880 г.) . Сент-Луис: Центр аэронавигационных карт и информации . Проверено 30 июля 2011 г. Английский перевод «Математических и физических теорий высшей геодезии» , Том 1 (Тойбнер, Лейпциг, 1880 г.).
• Карни, Чарльз Ф.Ф. (январь 2013 г.). «Алгоритмы геодезических работ» . Журнал геодезии . 87 (1): 43–55. arXiv : 1109.4448 . Бибкод : 2013JGeod..87...43K . дои : 10.1007/s00190-012-0578-z . Дополнения .
• Лежандр, Адриен-Мари (1806). «Анализ треугольников, нарисованных на поверхности сфероида» . Мемуары класса математических и физических наук Национального института Франции (1-й семестр): 130–161 . Проверено 30 июля 2011 г.
• Рейнсфорд, Х.Ф. (1955). «Длинная геодезическая на эллипсоиде». Бюллетень геодезии . 37 : 12–22. Бибкод : 1955BGeod..29...12R . дои : 10.1007/BF02527187 . S2CID   122111614 .
• Рэпп, Ричард Х. (март 1993 г.). Геометрическая геодезия. Часть II (Технический отчет). Университет штата Огайо . Проверено 1 августа 2011 г.
• Винсенти, Таддеус (апрель 1975 г.). «Прямые и обратные решения геодезических на эллипсоиде с применением вложенных уравнений» (PDF) . Обзор опроса . XXIII (176): 88–93. Бибкод : 1975SurRv..23...88V . дои : 10.1179/sre.1975.23.176.88 . Проверено 11 июля 2009 г. При выборе формулы решения геодезических линий первостепенное значение имеет длина программы, то есть объем ядра, который она будет занимать в компьютере вместе с тригонометрическими и другими необходимыми функциями.
• Винсенти, Таддеус (август 1975b). Геодезическое обратное решение между противоположными точками (PDF) (Технический отчет). Геодезическая исследовательская эскадрилья DMAAC. дои : 10.5281/zenodo.32999 .
• Винсенти, Таддеус (апрель 1976 г.). «Переписка». Обзор опроса . XXIII (180): 294.
• Справочное руководство по геоцентрическим данным Австралии (GDA) . Межправительственный комитет по съемке и картографии (ICSM). Февраль 2006 г. ISBN.  0-9579951-0-5 . Архивировано из оригинала (PDF) 26 июня 2009 г. Проверено 11 июля 2009 г.

Внешние ссылки

• Онлайн-калькуляторы от Geoscience Australia :
  • Винсенти Директ (пункт назначения)
  • Винсенти Инверс (расстояние между точками)
• Калькуляторы Национальной геодезической службы США :
  • Онлайновые и загружаемые на ПК утилиты для вычислений , включая прямые (прямые) и обратные задачи, как в двух, так и в трех измерениях (по состоянию на 1 августа 2011 г.).
• Онлайн-калькуляторы с исходным кодом JavaScript от Криса Венесса (лицензия Creative Commons Attribution):
  • Винсенти Директ (пункт назначения)
  • Винсенти Инверс (расстояние между точками)
• GeographicLib предоставляет утилиту GeodSolve (с исходным кодом, лицензированным MIT/X11) для решения прямых и обратных геодезических задач. По сравнению с Винсенти, это примерно в 1000 раз точнее (погрешность = 15 нм) и обратное решение является полным. Вот онлайн-версия GeodSolve .
• Полная реализация прямых и обратных формул Винсенти с исходным кодом, реализация Excel VBA Томаша Ястржебского.