Методы трансформации геодезических и пространственных прямоугольных координат, их алгоритмы, параметры, точность Текст научной статьи по специальности «Физика»

УДК 528.236.4

Вестник СПбГУ. Сер. 7. 2012. Вып. 3

Г. Д. Курошев, А. А. Харунжий

МЕТОДЫ ТРАНСФОРМАЦИИ ГЕОДЕЗИЧЕСКИХ И ПРОСТРАНСТВЕННЫХ ПРЯМОУГОЛЬНЫХ КООРДИНАТ, ИХ АЛГОРИТМЫ, ПАРАМЕТРЫ, ТОЧНОСТЬ

1. Введение

В современных инструментальных ГИС и прикладных геодезических программах реализовано большое количество методов трансформации систем координат. И несмотря на то что используются одни и те же алгоритмы, на практике в большинстве случаев такие программы представляют собою «черный ящик», так как код изначально закрыт, а руководство пользователя подчас не содержит исчерпывающей информации о параметрах и точности алгоритма. Таким образом, ставится задача исследования результатов работы той или иной программы на предмет корректности и точности. Цель данной статьи — произвести обзор имеющихся методов трансформации координат, дать априорную оценочную характеристику результатов их применения в зависимости от значений входных параметров и особенности реализации алгоритмов.

2. Методы трансформации координат

Рассмотрим геодезическую систему координат на эллипсоиде В, L, начало которой может быть как в центре земли (общеземной эллипсоид), так и смещено от него на некоторый вектор (референц-эллипсоид), а также прямоугольную геоцентрическую систему координат X, Y, Z. Связь между геодезическими и прямоугольными координатами с началом в общей точке и одинаковой ориентировкой осей легко находится методами дифференциальной геометрии и выглядит следующим образом [1, 2]:

X = (N + H) cos B cos L,

Y = (N + H) cos B sin L,

2 (1)

Z = ((1 -e )N + H)sinB,

a

N =

V1 - e2 sin2 B

где N — радиус кривизны эллипсоида первого вертикала, H — геодезическая высота.

Данные формулы не обратимы в элементарных функциях, обратный переход может быть осуществлен методом последовательных приближений с любой наперед заданной точностью. Вывод и описание этого алгоритма подробно изложены в [1, 2], для удобства дальнейшего изложения обозначим его (2).

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

© Г. Д. Курошев, А. А. Харунжий, 2012

образом, насчитывает семь возможных параметров — три угла поворота вокруг каждой из осей (ых, ыу, шг) три компонента вектора сдвига (ДХ, ДY, ДZ) и один масштабный коэффициент т. При различных способах выбора углов (например, углы Эйлера) и порядка в последовательности поворотов будут получаться разные матрицы поворота. За стандартный подход обычно принимают последовательное вращение вокруг осей X, Y, Z на углы (ых, шг), при этом угол считается положительным при вращении по часовой стрелке, если смотреть из начала координат вдоль положительного направления соответствующей оси. В случае малых значений углов значение синуса можно положить равным его аргументу, значение косинуса — единице, и тогда вид матрицы поворота существенно упрощается, вопрос о корректности подобного допущения будет обсуждаться в главе 4. Подобным образом получается наиболее распространенное в геодезии семипараметрическое преобразование Хельмерта [2]:

(3)

где угловые элементы (ых, ыу, шг) выражены в радианах.

Матрица поворота в формуле (3) легко обратима и обратный переход выглядит следующим образом:

f X ^ ( 1 Л 1 ®г -®y f X ^ ГДХ ^

Y = (1 + m) -®z 1 ®x Y + ДY

Z ® 1 Z ДZ

V У B V y x У А V У

f X ^

Y = (1 - m)

Z

V /А

1

v-®y

-®z 1

®y

-®x 1

f x ^ ГДХ ^

Y - ДY

Z V У B ДZ V У

(4)

Трансформацию систем координат с использованием формулы (2) в зарубежной литературе принято называть coordinate frame rotation (код EPSG: 9607 [3], ArcGIS Macro: PE_MTH_Coordinate_Frame [4]). Существует другой ее вариант, отличающийся от приведенных противоположным направлением вращения осей и, следовательно, инвертированными значениями углов поворота, он носит название формулы Бурсы-Вольфа, в зарубежной литературе — position vector (код EPSG: 9606 [3], ArcGIS Macro: PE_ MTH_Position_Vector [4]). В случае, когда угловые элементы и масштабный коэффициент равны нулю, матрица в (3) становится единичной и переход заключается лишь в пространственном сдвиге. В зарубежной литературе трансформация подобного рода носит название geocentric translation (код EPSG: 9601 [3], ArcGIS Macro: PE_MTH_Geo-centric_Translation [4]).

Существует расширенная модификация формулы (3), в которой добавляется возможность выбора точки вращения системы координат. Подобный подход потенциально может увеличить точность трансформации, но требует определения трех дополнительных параметров: Ха Yc, Zc — координат точки вращения в исходной системе. Эта трансформация, называемая десятипараметрической, задается следующей формулой:

f X ^ ' 1 ®z -®У ' f X - Xc

Y = (1 + m) -®z 1 ®x Y - Yc

Z V У B V® У -®x 1 У Z - V Zc

У А

' Л

Yc

V Zc У

АДХ л ДY Д Z

В зарубежных источниках десятипараметрическую трансформацию принято называть методом Молоденского-Бадекиса (код EPSG: 9636 [3], ArcGIS Macro: PE_MTH_ Molodensky_Badekas [4]).

Еще одной модификацией метода Хельмерта является трансформация, учитывающая скорости изменения семи параметров перехода во времени. Для их вычисления необходимо иметь сведения об изменениях в координатах опорных пунктов в результате переуравнивания на основе новых измерений. Поэтому данный метод используется в высокоточных сетях, на которых производятся регулярные наблюдения. Он задается четырнадцатью параметрами — семью параметрами метода Хельмерта и семью значениями скоростей их изменения, его полная формула представлена, например в [5].

Десяти и четырнадцатипараметрический методы используются нечасто, так как не всегда имеется достаточное количество данных для вычисления необходимых параметров трансформации. Однако математическое строение этих методов родственно методу Хельмерта и, следовательно, точностные характеристики их численной реализации также будут схожи (см. гл. 4).

Переход между системами геодезических координат может быть произведен последовательным применением формул (1), (3), (2); соответственно обратный переход — применением формул (1), (4), (2). Вместо формулы (3) может быть использовано десятипараметрическое (5) и четырнадцатипараметрическое преобразование, а их обращение производится аналогичным методу Хельмерта образом.

Существуют также прямые методы, не использующие прямоугольные системы координат как вспомогательные. Наиболее общий из них задается так называемыми дифференциальными формулами преобразования [2, 5, 6], которые получаются с использованием усеченных разложений функций в ряды. В качестве его входных параметров используются те же семь параметров перехода (ДХ, Д Y, AZ, wx, wy, wz, m) плюс геометрические параметры эллипсоидов исходной и конечной систем геодезических координат (значения эксцентриситетов — eA, eB и больших полуосей — aA, aB):

Бв = Ba + AB = Ba +

P

f ^2 a 2 л e2 sin B cos BAa + (— + 1)N sin B cos B--

a a~

(AX cos L + AY sin L)sin B + AZ cos B

„2___О D\ „ ___т /! , „2___о D\ _____2 .

(M + H)

- юх sin L(1 + e2 cos 2B) + юу cos L(1 + e2 cos 2B) - pme2 sin B cos B;

lb = la +al = la +

P

(N + H)cos B

(-AX sin L + AY cos L) +

(6)

+ tgB(1 -e2)(юх cosL + юу sinL)-az;

HB = HA +AH = HA - a Aa + N sin2 B — + (AX cos L + AY sin L)cos B + AZ sin B -

N

(

-Ne sin B cos B

• r T

sin L--— cos L

Л f 2 a

— + H

N

m,

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

2 + &B

e = , Да = aB - a

Де2 = е2в - e2,

B

2

М, N — радиусы кривизны меридианного сечения и первого вертикала соответственно, р — число угловых секунд в одном радиане (р = 206264,806").

Вычисления по формулам (6) производятся методом итераций. Для вычисления первого приближения значениям В, Ь, Н присваивают их значения в исходной системе, то есть ВА, LA, НА соответственно; в последующих итерациях используют арифметические средние между значениями В, L, Н на последнем и предпоследнем шаге:

Согласно [2] погрешность в вычислениях при выполнении одной итерации не превышает 0,3 м в линейной мере, при двух итерациях — 0,001 м. Данный метод рекомендован к использованию государственными стандартами РФ [2].

Упрощенная версия формул (6), в которой не учитываются члены, содержащие углы поворота осей (wx, wy, wz) и масштабный коэффициент m, а также с несколько видоизмененным слагаемым, содержащим sinBcosB, называется методом Молоденско-го (код EPSG: 9604 [3], ArcGIS Macro: PE_MTH_Molodensky [4]). Эта трансформация реализует только пространственный сдвиг аналогично преобразованию geocentric translation. Усеченный метод Молоденского, используемый для упрощения вычислений, не учитывает еще и сдвиг по высоте (код EPSG: 9605 [3], ArcGIS Macro: PE_MTH_ Molodensky_Abridged [4]). Полные формулы этих методов приведены, например в [5], а вычисления по ним ведутся в один прием, без итераций.

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

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

3. Наборы параметров трансформации

параметров трансформаций требует подробного изучения методики и качества произведенных геодезических работ и выходит за рамки данного исследования. Здесь мы попробуем дать оценку ожидаемой разницы при использовании разных наборов параметров трансформации между двумя заданными системами координат.

В качестве иллюстрации рассмотрим наборы параметров трансформации российской референцной системы координат СК-42 в общеземную систему WGS-84 [2-4, 7, 8]. Этот переход часто выполняется в процессе решения типовых картографических задач, связанных с использованием данных спутниковой навигации, и в различных программах представлен более чем пятью возможными наборами параметров (табл. 1). Согласно государственному стандарту на территории Российской Федерации для трансформации координат следует производить последовательный переход от СК-42 к ПЗ-90 или к ПЗ-90.02 и уже потом к WGS-84; численные значения параметров этих трансформаций даны в табл. 2 [9, 10, 11].

Таблица 1. Параметры трансформации координат из СК-42 в WGS-84 (обозначение типа трансформации: g.t. — geocentric translation, c.f.r. — coordinate frame

rotation)

№ п/п ГИС и название трансформации Тип AX, м AY, м AZ, м шу, " ^ " m-106

1 NIMA-2000; EPSG:1254; ArcGIS: Pulkovo_1942_ To_WGS_1984; ERDAS Img.: Pulkovo 1942-1 (Russia) g.t. 28 -130 -95 0 0 0 0

2 MapInfo 1001 c.f.r. 24 -123 -94 -0,02 0,25 0,13 0,11

3 ERDAS Img.: Pulkovo 1942 c.f.r. 27 -135 -84,5 0 0 -0,554 0,2263

4 EPSG:15865 c.f.r 25 -141 -78,5 0 -0,35 -0,736 0

Таблица 2. Параметры трансформации координат согласно ГОСТ Р 51794-2008 (обозначение типа трансформации: g.t. — geocentric translation, c.f.r. — coordinate frame

rotation)

№ п/п ГИС и название трансформации Тип AX, м AY м AZ, м " шу, " ^ " 6 m-10

1 Из СК-42 в ПЗ-90.02 c.f.r. 23,93 -141,03 -79,98 0 -0,35 -0,79 -0,22

2 Из ПЗ-90.02 в WGS-84 g.t. -0,36 0,08 0,18 0 0 0 0

3 Из СК-42 в ПЗ-90 c.f.r. 25 -141 -80 0 -0,35 -0,66 0

4 Из ПЗ-90 в WGS-84 c.f.r -1,1 -0,3 -0,9 0 0 -0,2 -0,12

Все расчеты, результаты которых приведены ниже, велись в разработанной А. А. Харунжим программе. При этом для каждого набора параметров производилось сравнение с результатами аналогичных преобразований в 9.2, Мар!п6э 10.0,

PHOTOMOD GeoCalculator, Credo_ТРАНСКОР. Во всех случаях совпадение происходило не менее, чем до седьмого знака после запятой в значениях В, Ь, выраженных в десятичных градусах. Такую точность можно считать удовлетворительной (см. гл. 4). Величина расстояния между образами преобразований вычислялась как длина геодезической линии на эллипсоиде, графическая интерпретация полученных данных реализована в среде ArcGIS 9.2.

Первое замечание, которое необходимо сделать, состоит в том, что последовательное использование трансформаций, заданных параметрами первой и второй строк табл. 2, не дает удовлетворительного схождения с трансформацией по третьей и четвертой строкам той же таблицы, хотя параметры заданы одним источником и теоретически должны быть согласованы. Распределение разницы в полученных координатах в графическом виде дано на рис. 1; расстояние на эллипсоиде может достигать двух и более метров. Полученный результат дополнительно указывает на некорректность назначения какого-то определенного набора параметров трансформации оптимальным в глобальном смысле.

Рис. 1. Распределение разницы (м) в трансформированных координатах из СК-42 в WGS-84 по параметрам первой, второй строк табл. 2 и третьей, четвертой строк табл. 2

Для дальнейшего исследования будем считать трансформацию через ПЗ-90 (по третьей и четвертой строкам табл. 2) эталонной. На рис. 2-5 в графическом виде представлены расхождения в координатах при использовании параметров эталонного перехода и параметров, данных в табл. 1. Как видно, величина расхождения и ее распределение различны, максимум варьируется от 6 до 58 м, что в общем случае является недопустимой величиной для картографирования в масштабе вплоть до 1:300 000.

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

0,000-9,000 | | 9,001-18,000 I I 18,001-27.000 И 27.001-36,000 И 36.001-46.000

Рис. 2. Распределение разницы (м) в трансформированных координатах из СК-42 в ШОБ-84 по параметрам первой строки табл. 1 и третьей, четвертой строк табл. 2

0,000-12,000 [ | 12,001-23,000 | | 23,001-35,000 Щ 35,001-47,000 Щ 47,001-58,000

Рис. 3. Распределение разницы (м) в трансформированных координатах из СК-42 в ШОБ-84 по параметрам второй строки табл. 1 и третьей, четвертой строк табл. 2

Рис. 4. Распределение разницы (м) в трансформированных координатах из СК-42 в WGS-84 по параметрам третьей строки табл. 1 и третьей, четвертой строк табл. 2

Рис. 5. Распределение разницы (м) в трансформированных координатах из СК-42 в WGS-84 по параметрам четвертой строки табл. 1 и третьей, четвертой строк табл. 2

В работе не ставится цель выяснить, какой именно набор параметров следует считать оптимальным, но на основе полученных данных можно дать оценку ожидаемой разницы при использовании различных наборов параметров и методов трансформации. В любом случае для конкретных задач следует выбирать тот тип и параметры трансформации, которые дают минимальную разницу в координатах на ближайших к картографируемой местности опорных пунктах. Например, набор параметров в строке первой табл. 1 рекомендован к использованию исключительно на территории европейской части России. При наличии достаточного количества достоверных данных имеет смысл вычислять параметры трансформации для локальных областей самостоятельно методом наименьших квадратов. Данная функция реализована в программном продукте А. А. Харунжего, а также во многих прикладных геодезических программах, например Сге^_ТРАНСКОР.

4. Точность реализации алгоритма трансформации

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

Измерения и расчеты в ходе решения геодезических задач обычно ведутся с учетом величин порядка 10-3. В угловых значениях широты и долготы эта величина соответствует 10-8 десятичных долей градуса. Точность определения координат для задач картографирования ограничена графической читаемостью и соответствует величине 0,2-0,4 мм в масштабе карты, т. е. имеет порядок 2-4 • М • 10-4 м, где М — знаменатель масштаба. Соответственно для крупномасштабного картографирования в масштабе 1:500 эта величина равна 0,1-0,2 м, для масштаба 1:50 000 — 10 м. Поскольку основная часть задач, решаемых при помощи ГИС, имеет картографический характер, можно считать погрешность, имеющую порядок меньше, чем 10-2 м, несущественной.

На точность работы любого численного алгоритма влияет выбор количества разрядов в представлении числа с плавающей точкой. Наиболее распространены на данный момент следующие типы: real (6 байт, 11-12 значащих цифр), double (8 байт, 15-16 значащих цифр) и extended (10 байт, 19-20 значащих цифр). Относительный порядок желаемой точности и максимальных величин в задачах трансформации (107 м — порядок полуосей земного эллипсоида и, соответственно, максимума абсолютных значений пространственных координат) указывает на то, что для качественных расчетов необходимо по крайней мере 9 значащих цифр. Следовательно, даже с использованием типа данных real должен получаться приемлемый результат. Это предположение частично подтвердилось в ходе практических расчетов: при однократном применении любой из операций, описанных в гл. 2, разница в первых девяти цифрах результата расчета не зависела от выбора между real, double и extended. При многократном циклическом трансформировании (30 и более прямых и обратных переходов) ошибка накапливается и начинает превышать допустимую. Поэтому рекомендуется выполнять расчеты, используя как основной тип данных double или extended.

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

можно сказать и о методе Хельмерта (3) вместе со всеми его модификациями. Единственный нюанс состоит в моменте замены синуса значением его аргумента, а косинуса — единицей. При подробном рассмотрении заметим, что величина углов поворота для большинства трансформаций не превышает десятых долей секунды, что соответствует 10-7 рад. При значениях такого порядка разность между синусом и его аргументом имеет порядок 10-22. После умножения на координату порядок увеличивается максимум до 1015, что, согласно ранее произведенным расчетам, является совершенно несущественной величиной. Вычислительные эксперименты с методом Хельмерта подтвердили факт появления разницы не более, чем в пятнадцатой значащей цифре.

Особое внимание следует уделить исследованию точности итерационных алгоритмов, коим является, например, переход от прямоугольных координат к геодезическим [1, 2]. Существует несколько вариантов такого алгоритма, все они используют некоторую малую константу в качестве допуска, останавливающего проведение итераций [1, 2, 6, 8]. В ходе вычислительных экспериментов выяснено, что заявленная точность иногда обеспечивается лишь при однократном выполнении алгоритма, поэтому имеет смысл уменьшать рекомендуемый допуск на один-два порядка даже в ущерб скорости работы алгоритма.

Дополнительно следует отметить, что все вышеупомянутые алгоритмы не показали принципиальной чувствительности к диапазону значений аргументов В, Ь, Н и X, Y, X. Этот факт не распространяется на трансформацию с использованием дифференциальных формул (6) и метода Молоденского. На рис. 6 в графическом виде представлено распределение разницы в координатах между результатами семипараметрической трансформации по формуле (6) и по формулам (1), (3), (2); использовался набор параметров четвертой строки табл. 1. На рис. 7 дан аналогичный пример распределения

(1 1 3 Г Й^3 а 4-,

... > ■V ш « у

т о 1ГМ Ь У £

\ г Л £

ч IV А. Ж > > .

■{х J

Ч. I П г® у. их Л ■гО

N л. и. л,

V ■т >

Г г У г

я уу * —' "л > V,

¥ С г

-

___ - ■ =: Й

0,000-0,003 | | 0,003-0,009 | | 0,009-0,019 Щ 0,019-0,040 Щ 0,040-0,052

Рис. 6. Распределение разницы в координатах между результатами семипараметриче-ской трансформации с помощью дифференциальных формул и методом Хельмерта для набора параметров из четвертой строки табл. 1

Рис. 7. Распределение разницы в координатах между результатами трехпараметриче-ской трансформации с помощью дифференциальных формул и методом Хельмерта для набора параметров из первой строки табл. 1

Рис. 8. Распределение разницы в координатах между результатами трехпараметриче-ской трансформации методом Молоденского и методом Хельмерта для набора параметров из первой строки табл.1

разницы для трехпараметрического преобразования, заданного параметрами первой строки табл. 1.

Как видно, разница между результатами применения дифференциальных формул и метода Хельмерта почти всюду мала. Согласно выдвинутому ранее заключению о нечувствительности переходов (1), (2), (3) к диапазону значений аргументов, можно сделать вывод, что формулы (6) несколько теряют точность в приполярных районах, а в широтах больших 89° не могут применяться в принципе. Этот факт легко объясним самим видом формул: при стремлении значения B к ±90° второй член формулы (6) для L стремится к бесконечности. Аналогично ситуация обстоит с методом Молоденского и его модификациями, что наглядно проиллюстрировано на рис. 8.

Можно сделать вывод, что метод Хельмерта и методы, основанные на дифференциальных формулах, везде кроме полярных районов дают идентичные результаты, поскольку порядок разницы в трансформированных координатах не превышает 10-3 м.

Литература

1. Морозов В. П. Курс сфероидической геодезии. М.: Недра, 1979. 296 с.

2. Глобальные навигационные спутниковые системы. Системы координат. Методы преобразований координат определяемых точек. ГОСТ Р 51794-2008. М.: Стандартинформ, 2009. 16 с.

3. EPSG Geodetic Parameter Dataset, Version 7 — June 2011 [Электронный ресурс]. URL: http:// www.epsg.org/guides/docs/G7-1.pdf (дата обращения: 15.09.2011)

4. ArcGIS: Supported geographic (datum) transformations [Электронный ресурс]. URL: http:// downloads.esri.com/support/whitepapers/ao_/geographic_transformations9.pdf (дата обращения: 15.09.2011)

5. Комаровский Ю. А. Использование различных референц-эллипсоидов в судовождении: учебное пособие. Владивосток: ИПК МГУ им. адм. Г. И. Невельского, 2005. 342 с.

6. Герасимов А. П., Назаров В. Г. Местные системы координат. М.: Издательство Проспект, 2010. 62 с.

7. MapInfo Professional 10.0, руководство пользователя [Электронный ресурс]. М.: ООО "ЭСТИ МАП", 2009. 500с. URL: http://74.54.97.210/~estimap/download/docs/mi100ug.pdf (дата обращения: 15.09.2011)

8. National imagery and mapping agency technical report 8352.2 third edition [Электронный ресурс]. MD USA, 2000. URL: http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin. pdf (дата обращения: 15.09.2011)

9. Галазин Л. А., Каплан Б.Л., Лебедев М. Г. и др. Система геодезических параметров Земли «Параметры Земли 1990 года» (ПЗ-90): справочный документ. М.: Координационный научно-информационный центр, 1998. 37 с.

10. Введение новой государственной референцной системы геодезических координат 1995 года (СК-95). Технико-экономический доклад. М.: ЦНИИГАИК, 1998. 72 с.

11. Руководство пользователя по выполнению работ в системе координат 1995 года (СК-95). М.: ЦНИИГАИК, 2004. 84 с.

Статья поступила в редакцию 23 марта 2012 г.