Преобразование Геодезических Систем Координат

Применение различных систем координат при решении практических задач геодезии, картографии, навигации и в геоинформационных системах неизбежно. Различают несколько типов преобразования координат: переход между различными форматами координат, переход между различными системами координат и картографическими проекциями, а также преобразование датумов. Все перечисленные виды преобразования будут рассмотрены в данной статье.

Обозначить географическое место обычно значит передать широту и долготу места. Числовые значения для широты и долготы могут быть представлены в нескольких различных видах единиц и форматах:

шестидесятеричная система: градусы, минуты и секунды: 40° 26′ 46″ N 79° 58′ 56″ W

градусы и десятичные минуты: 40° 26.767′ N 79° 58.933′ W

десятичные градусы: 40,446° N 79,982° W

В градусе 60 минут, а в минуте 60 секунд. Поэтому для перевода из формата градусы/минуты/секунды в формат десятичных градусов можно использовать формулу:

десятичные градусы=градусы+минуты/60+секунды/3600.

Для обратного перевода из формата десятичных градусов в формат градусы/минуты/секунды можно использовать формулы:

градусы=[десятичные градусы]

минуты=[60*(десятичные градусы-градусы)]

секунды=3600*(десятичные градусы-градусы)-60*минуты

где обозначение [x] значит, что надо взять целую часть x, и обратиться к «полочной функции».

Преобразование системы координат — это переход из одной системы координат в другую, причём обе системы координат основаны на одном и том же геодезическом датуме. Часто задача преобразования заключается в переходе от геодезической системы координат к прямоугольным координатам либо переход от одной картографической проекции к другой.

Из геодезической системы координат в прямоугольную

Прямоугольные координаты точек в пространстве можно вычислить по известным геодезическим координатам этих точек (широта B, долгота L, высота H) по формулам:

${\displaystyle {\begin{aligned}X&=\left(N(B)+H\right)\cos {B}\cos {L}\\Y&=\left(N(B)+H\right)\cos {B}\sin {L}\\Z&=\left({\frac {b^{2}}{a^{2}}}N(B)+H\right)\sin {B}\end{aligned}}}$ 

где

${\displaystyle N(B)={\frac {a^{2}}{\sqrt {a^{2}\cos ^{2}B+b^{2}\sin ^{2}B}}}={\frac {a}{\sqrt {1-e^{2}\sin ^{2}B}}},}$ 

где ${\displaystyle a}$  и ${\displaystyle b}$  — экваториальный (большая полуось) и полярный радиусы (малая полуось), соответственно. ${\displaystyle e^{2}={\frac {a^{2}-b^{2}}{a^{2}}}}$  — квадрат первого эксцентриситета эллипсоида. ${\displaystyle \,N(B)}$  радиус кривизны первого вертикала — расстояние по нормали к эллипсоиду от точки пересечения поверхности эллипсоида нормалью до оси oZ (Рис. 1).

Из прямоугольной системы координат в геодезическую

При переходе от прямоугольных пространственных координат к геодезической системе координат (такой как WGS84) геодезические широты B и высоты H зачастую требуется вычислять итеративными методами, то есть выполнением последовательных приближений. Что касается долгот L, то они вычисляются обычным путём.

Существует несколько методов для вычислений геодезических широт и высот, рассмотрим два из них.

Метод Ньютона-Рафсона

Следующее иррациональное уравнение Боуринга для геодезической широты решается итерационным методом Ньютона — Рафсона:

${\displaystyle \kappa -1-{\frac {e^{2}a\kappa }{\sqrt {p^{2}+(1-e^{2})Z^{2}\kappa ^{2}}}}=0,}$ 

где ${\displaystyle p={\sqrt {x^{2}+y^{2}}}}$ ,

Широту B можно найти из уравнения ${\displaystyle \kappa ={\frac {p}{Z}}\tan(B)}$ .

Высота H рассчитывается как:

${\displaystyle H=e^{-2}(\kappa ^{-1}-{\kappa _{0}}^{-1}){\sqrt {p^{2}+Z^{2}\kappa ^{2}}}.}$ 

Итерация может быть преобразована к следующему виду:

${\displaystyle \kappa _{i+1}={\frac {c_{i}+\left(1-e^{2}\right)Z^{2}\kappa _{i}^{3}}{c_{i}-p^{2}}}=1+{\frac {p^{2}+\left(1-e^{2}\right)Z^{2}\kappa _{i}^{3}}{c_{i}-p^{2}}},}$ 

где ${\displaystyle c_{i}={\frac {\left(p^{2}+\left(1-e^{2}\right)Z^{2}\kappa _{i}^{2}\right)^{3/2}}{ae^{2}}},}$  ${\displaystyle \kappa _{0}=\left(1-e^{2}\right)^{-1}.}$ 

Постоянная ${\displaystyle \,\kappa _{0}}$  является хорошим начальным значением для итерации, когда ${\displaystyle H\approx 0}$ . Боуринг показал, что в таком случаи уже первая итерация дает достаточно точное решение. Он использовал дополнительные тригонометрические функции в своей первоначальной формулировке.

Решение Ferrari

Вышеизложенное уравнение для ${\displaystyle \kappa }$  можно решить методом Феррари:

${\displaystyle {\begin{aligned}\zeta &=(1-e^{2})z^{2}/a^{2},\\[6pt]\rho &=(p^{2}/a^{2}+\zeta -e^{4})/6,\\[6pt]s&=e^{4}\zeta p^{2}/(4\rho ^{3}a^{2}),\\[6pt]t&={\sqrt[{3}]{1+s+{\sqrt {s(s+2)}}}},\\[6pt]u&=\rho (1+t+1/t),\\[6pt]v&={\sqrt {u^{2}+e^{4}\zeta }},\\[6pt]w&=e^{2}(u+v-\zeta )/(2v),\\[6pt]\kappa &=1+e^{2}({\sqrt {u+v+w^{2}}}+w)/(u+v).\end{aligned}}}$ 

Применение решения Феррари

Существует ряд методов и алгоритмов, но наиболее точным, согласно Чжу, является следующая последовательность, установленная Хейккиненом. Предполагается, что геодезические параметры известны.

${\displaystyle {\begin{matrix}r&=&{\sqrt {X^{2}+Y^{2}}}\\e'^{2}&=&(a^{2}-b^{2})/b^{2}\\F&=&54b^{2}Z^{2}\\G&=&r^{2}+(1-e^{2})Z^{2}-e^{2}(a^{2}-b^{2})\\c&=&{\frac {e^{4}Fr^{2}}{G^{3}}}\\s&=&{\sqrt[{3}]{1+c+{\sqrt {c^{2}+2c}}}}\\P&=&{\frac {F}{3\left(s+{\frac {1}{s}}+1\right)^{2}G^{2}}}\\Q&=&{\sqrt {1+2e^{4}P}}\\r_{0}&=&{\frac {-(Pe^{2}r)}{1+Q}}+{\sqrt {{\frac {1}{2}}a^{2}\left(1+1/Q\right)-{\frac {P(1-e^{2})Z^{2}}{Q(1+Q)}}-{\frac {1}{2}}Pr^{2}}}\\U&=&{\sqrt {(r-e^{2}r_{0})^{2}+Z^{2}}}\\V&=&{\sqrt {(r-e^{2}r_{0})^{2}+(1-e^{2})Z^{2}}}\\z_{0}&=&{\frac {b^{2}Z}{aV}}\\H&=&U\left(1-{\frac {b^{2}}{aV}}\right)\\B&=&\arctan \left[{\frac {Z+e'^{2}z_{0}}{r}}\right]\\L&=&\arctan 2[Y,X]\end{matrix}}}$ 

Примечание: arctan2 [Y, X] — обратная касательная к четырем квадрантам.

Степенной ряд

Для малых е² степенной ряд ${\displaystyle \kappa =\sum _{i\geq 0}\alpha _{i}e^{2i}}$  начинается с

${\displaystyle \alpha _{0}=1;}$ 
${\displaystyle \alpha _{1}={\frac {a}{\sqrt {Z^{2}+p^{2}}}};}$ 
${\displaystyle \alpha _{2}={\frac {aZ^{2}{\sqrt {Z^{2}+p^{2}}}+2a^{2}p^{2}}{2(Z^{2}+p^{2})^{2}}}.}$ 

Переход из геодезической системы координат в ENU и обратно

Преобразование из геодезических координат в топоцентрические координаты ENU состоит их двух этапов:

1. Преобразование координат из геодезической системы в прямоугольную.
2. Преобразование координат из прямоугольной в топоцентрическую систему координат ENU.

Преобразование координат из прямоугольной в топоцентрическую систему координат ENU

Чтобы преобразовать прямоугольные координаты в топоцентрические необходимо знать начальную точку топоцентрической системы координат, обычно она располагается в некоторой точке наблюдений. Если наблюдение производится в точке ${\displaystyle \{X_{r},Y_{r},Z_{r}\}}$  , а наблюдаемый объект в ${\displaystyle \{X_{p},Y_{p},Z_{p}\}}$  тогда радиус-вектор этого направления в системе координат ENU имеет вид:

${\displaystyle {\begin{bmatrix}x\\y\\z\\\end{bmatrix}}={\begin{bmatrix}-\sin L_{r}&\cos L_{r}&0\\-\sin B_{r}\cos L_{r}&-\sin B_{r}\sin L_{r}&\cos B_{r}\\\cos B_{r}\cos L_{r}&\cos B_{r}\sin L_{r}&\sin B_{r}\end{bmatrix}}{\begin{bmatrix}X_{p}-X_{r}\\Y_{p}-Y_{r}\\Z_{p}-Z_{r}\end{bmatrix}}}$ 

Преобразование координат из топоцентрической системы координат ENU в прямоугольную.

Путём обратного преобразования координат из прямоугольной системы получим топоцентрическую систему координат:

${\displaystyle {\begin{bmatrix}X\\Y\\Z\\\end{bmatrix}}={\begin{bmatrix}-\sin \lambda &-\sin \phi \cos \lambda &\cos \phi \cos \lambda \\\cos \lambda &-\sin \phi \sin \lambda &\cos \phi \sin \lambda \\0&\cos \phi &\sin \phi \end{bmatrix}}{\begin{bmatrix}x\\y\\z\end{bmatrix}}+{\begin{bmatrix}X_{r}\\Y_{r}\\Z_{r}\end{bmatrix}}}$ 

Переход к другой картографической проекции

Преобразование координат и положений на карте между различными картографическими проекциями карты, привязанные к одной и тоже геодезической поверхности, может быть выполнено либо с помощью формул прямого перехода из одной проекции в другую, либо сначала проекция ${\displaystyle A}$  преобразуется в промежуточную систему координат, такую ​​как прямоугольная, а уже из неё в проекцию ${\displaystyle B}$ . Используемые формулы могут быть сложными, в некоторых случаях преобразование не имеет решения в замкнутой форме, и необходимо использовать приближенные методы. Обычно для выполнения задач преобразования координат используются компьютерные программы, например, с программой GEOTRANS, поддерживаемой DoD и NGA.

 

Преобразования между датумами могут быть выполнены различными способами. Существуют преобразования, которые позволяют совершить прямой переход от геодезических координат одного датума к геодезическим координатам другого датума. Существуют менее прямые переходы, которые преобразуют геодезические координаты в геоцентрические (ECEF), преобразуют геоцентрические координаты из одного датума в другой, затем преобразуют геоцентрические координаты другого датума обратно в геодезические. Также существуют проекционные преобразования, которые позволяют совершить прямой переход из одной (датум, проекция) пары к другой (датум, проекция) паре.

Проекционные преобразования

 

Проекционные преобразования позволяют совершить прямой переход от координат на карте для одной (картографическая проекция, датум) пары к координатам на карте для другой (картографическая проекция, датум) пары. Как пример можно привести метод NADCON для преобразования от Северо — Американского датума (North American Datum) (NAD) 1927 года к датуму NAD^[en] 1983 года. Высокоточная эталонная сеть (The High Accuracy Reference Network) (HARN), высокоточная версия преобразований NADCON, имеет точность приблизительно 5 сантиметров. Национальное преобразование версия 2(The National Transformation version 2) (NTv2) является канадской версией NADCON для перехода между NAD 1927 и NAD 1983^[en]. Методы HARN также известны, как NAD 83/91 и Высокоточные проекционные сети (High Precision Grid Networks) (HPGN). Впоследствии, Австралия и Новая Зеландия адаптировали для себя формат NTv2 для того, чтобы создать методы проекционных преобразований для переходов между их собственными местными датумами.

Как и преобразования с помощью уравнений множественной регрессии, проекционные методы используют интерполяцию низких порядков для преобразования координат карты, но в двух пространствах вместо трех. NOAA предоставляет программное обеспечение (как часть NGS Geodetic Toolkit) для производства преобразований NADCON.

Преобразование Молоденского

Преобразование Молоденского позволяет совершить прямой переход между геодезическими координатами разных датумов без необходимости промежуточного перехода к геоцентрическим координатам. Для него требуются три смещения между центрами систем координат и разницы между большими полуосями и параметрами сжатия референц-эллипсоидов.

Преобразование Молоденского используется Национальным агентством геопространственной разведки (National Geospatial-Intelligence Agency) (NGA) в их техническом сообщении TR8350.2, а также в программе GEOTRANS, поддерживаемой агентством. Преобразование Молоденского было популярно до прихода эры современных компьютеров, и метод является частью многих геодезических программ.

Уравнения множественной регрессии

Преобразования датумов с использованием эмпирических методов множественной регрессии были созданы для достижения большей точности для маленьких географических регионов чем стандартные преобразования Молоденского. Данные преобразования используются для преобразования местных датумов, которые создаются для континентов или меньших регионов, в глобальные датумы, такие как WGS 84. В стандарте NIMA TM 8350.2, Приложение D перечислены преобразования с помощью уравнений множественной регрессии от нескольких местных датумов к WGS 84, с точностью около 2 метров.

Метод уравнений множественной регрессии позволяет совершить прямое преобразование геодезических координат без промежуточного перевода в геоцентрические координаты. Геодезические координаты ${\displaystyle \phi _{B},\lambda _{B},h_{B}}$  в новом датуме B моделируются как полиномы вплоть до девятой степени в геодезических координатах ${\displaystyle \phi _{A},\lambda _{A},h_{A}}$ изначального датума A. Например, приращение ${\displaystyle \phi _{B}}$ может быть разложено как (показано разложение только до квадратичных членов):

${\displaystyle \Delta \phi =a_{0}+a_{1}U+a_{2}V+a_{3}U^{2}+a_{4}UV+a_{5}V^{2}+\cdots }$ 

где

${\displaystyle {\begin{aligned}a_{i}&={\mbox{параметры подобранные множественной регрессией}}\\U&=K(\phi _{A}-\phi _{m})\\V&=K(\lambda _{A}-\lambda _{m})\\K&={\mbox{фактор масштаба}}\\\phi _{m},\lambda _{m}&={\mbox{начало датума }}A\\\end{aligned}}}$ 

для ${\displaystyle \Delta \lambda }$  и ${\displaystyle \Delta h}$  выстраиваются схожие уравнения. При достаточном количестве пар координат (A, B) для точек в обоих датумах для хорошей статистики используются множественные методы регрессии для подбора параметров этих полиномов. Полиномы, вместе с подобранными коэффициентами, формируют уравнения множественной регрессии.

Преобразование Гельмерта

Использование преобразования Гельмерта^[en] при переходе от геодезических координат датума ${\displaystyle A}$  к геодезическим координатам датума ${\displaystyle B}$  происходит в три шага:

1 Преобразование геодезических координат датума ${\displaystyle A}$  в геоцентрические;

2 Применение преобразования Гельмерта, с соответствующими параметрами преобразования для ${\displaystyle A\to B}$ , для перехода от геоцентрических координат датума ${\displaystyle A}$  к геоцентрическим координатам датума ${\displaystyle B}$ ;

3 Преобразование геоцентрических координат в геодезические координаты для датума ${\displaystyle B}$ .

Для геоцентрических координат XYZ преобразование Гельмерта^[en] имеет форму:

${\displaystyle {\begin{bmatrix}X_{B}\\Y_{B}\\Z_{B}\end{bmatrix}}={\begin{bmatrix}c_{x}\\c_{y}\\c_{z}\end{bmatrix}}+(1+s\times 10^{-6})\cdot {\begin{bmatrix}1&-r_{z}&r_{y}\\r_{z}&1&-r_{x}\\-r_{y}&r_{x}&1\end{bmatrix}}\cdot {\begin{bmatrix}X_{A}\\Y_{A}\\Z_{A}\end{bmatrix}}.}$ 

Преобразование Гельмерта это преобразование с семью элементами, с тремя параметрами смещения ${\displaystyle c_{x},c_{y},c_{z}}$ , тремя параметрами разворота ${\displaystyle r_{x},r_{y},r_{z}}$  и одним масштабным параметром ${\displaystyle s}$ . Преобразование Гельмерта это приближенный метод, который можно считать точным только, когда параметры преобразования малы по сравнению с величинами векторов геоцентрической системы координат. С такими условиями преобразование можно считать обратимым.

Преобразование Гельмерта с четырнадцатью параметрами, с линейной зависимостью от времени для каждого параметра, может быть использовано для наблюдений за изменением во времени географических координат из-за геоморфологических процессов, таких как континентальный дрейф и землетрясения. Оно было преобразовано в программное обеспечение, такое как инструмент Horizontal Time Dependent Positioning (HTDP) в программном обеспечении U.S. NGS.

Преобразование Молоденского — Бадекаса

Для исключения связи между смещениями и разворотами преобразования Гельмерта могут быть использованы три дополнительных параметра, чтобы получить новый XYZ центр разворота ближе к преобразуемым координатам. Это преобразование с десятью параметрами называется преобразованием Молоденского — Бадекаса, его не стоит путать с более простым преобразованием Молоденского.

Как при использовании преобразования Гельмерта, использование преобразования Молоденского — Бадекаса состоит из трех шагов:

1. Преобразование геодезических координат датума ${\displaystyle A}$  в геоцентрические.
2. Применение преобразования Молоденского — Бадекаса, с соответствующими параметрами преобразования для ${\displaystyle A\to B}$ , для перехода от геоцентрических координат датума ${\displaystyle A}$  к геоцентрическим координатам датума ${\displaystyle B}$ .
3. Преобразование геоцентрических координат в геодезические координаты для датума ${\displaystyle B}$ .

Преобразование имеет форму:

${\displaystyle {\begin{bmatrix}X_{B}\\Y_{B}\\Z_{B}\end{bmatrix}}={\begin{bmatrix}X_{A}\\Y_{A}\\Z_{A}\end{bmatrix}}+{\begin{bmatrix}\Delta X_{A}\\\Delta Y_{A}\\\Delta Z_{A}\end{bmatrix}}+{\begin{bmatrix}1&-r_{z}&r_{y}\\r_{z}&1&-r_{x}\\-r_{y}&r_{x}&1\end{bmatrix}}\cdot {\begin{bmatrix}X_{A}-X_{A}^{0}\\Y_{A}-Y_{A}^{0}\\Z_{A}-Z_{A}^{0}\end{bmatrix}}+\Delta S{\begin{bmatrix}X_{A}-X_{A}^{0}\\Y_{A}-Y_{A}^{0}\\Z_{A}-Z_{A}^{0}\end{bmatrix}}.}$ 

где ${\displaystyle (X_{A}^{0},~Y_{A}^{0},~Z_{A}^{0})}$  — начало отсчёта для разворота и масштабного преобразования, а ${\displaystyle \Delta S}$  — масштабный фактор.

Преобразование Молоденского — Бадекаса используется для преобразования местных геодезических датумов в глобальные датумы, такие как WGS 84. В отличие от преобразования Гельмерта, преобразование Молоденского — Бадекаса необратимо из-за того, что начало отсчёта для разворота относится к изначальному датуму.

• Геодезическая система координат
• Прямоугольная система координат
• Список картографических проекций
• Проекция Гаусса — Крюгера
• Универсальная поперечная проекция Меркатора