CC BY
302
ГЕОДЕЗИЯ И МАРКШЕЙДЕРИЯ
УДК 528.2
АЛГОРИТМЫ НЕПОСРЕДСТВЕННОГО ВЫЧИСЛЕНИЯ ГЕОДЕЗИЧЕСКОЙ ШИРОТЫ И ГЕОДЕЗИЧЕСКОЙ ВЫСОТЫ ПО ПРЯМОУГОЛЬНЫМ КООРДИНАТАМ
Павел Александрович Медведев
Омский государственный аграрный университет, 644008, Россия, г. Омск, Институтская пл., 1, доктор технических наук, профессор кафедры геодезии и дистанционного зондирования, тел. (381)265-26-72, e-mail: omgau-math@rambler.ru
Борис Тимофеевич Мазуров
Сибирский государственный университет геосистем и технологий, 630108, Россия, г. Новосибирск, ул. Плахотного, 10, доктор технических наук, профессор кафедры физической геодезии и дистанционного зондирования, тел. (383)343-29-11, e-mail: btmazurov@mail.ru
В настоящее время активно осуществляется практическая реализация государственной геодезической системы координат 2011 г. на территории Российской Федерации. В связи с этим объективно существует необходимость уточнения математической и методологической основы сравнения параметров земного эллипсоида в государственных системах координат, схем преобразования координат, решения проблем, возникающих при преобразовании координат из местных систем координат в единую государственную. При этом должны учитываться результаты выполнения программы по построению современной спутниковой государственной геодезической сети России трех уровней (ФАГС, ВГС и СГС-1), а также точность ее связи с геодезическими сетями триангуляции и полигонометрии 1-4-го классов. При построении глобальных геодезических сетей, связанных с обработкой спутниковых наблюдений, и решении геодезических задач по определению положения точек земной поверхности применяются пространственные прямоугольные координаты X, Y, Z и геодезические координаты B, L, H. В связи с этим возникает задача по преобразованию этих систем координат. В статье математически обоснованы высокоточные неитеративные алгоритмы для прямого вычисления геодезической широты и высоты по исходным данным - пространственным прямоугольным координатам.
Ключевые слова: математика, астрономия, картография, геодезия, геодинамика.
ALGORITHMS FOR DIRECT COMPUTATION GEODETIC LATITUDE AND GEODETIC HEIGHT IN RECTANGULAR COORDINATES
Pavel A. Medvedev
Omsk State Agrarian University, 644008, Russia, Omsk, 2 Institutskaya Sq., D. Sc., Professor, Department of Geodesy and Remote Sensing, tel. (381)265-26-72, e-mail: omgau-math@rambler.ru
Boris T. Mazurov
Siberian State University of Geosystems and Technologies, 630108, Russia, Novosibirsk, 10 Plakhotnogo St., D. Sc., Professor, Department of Physical Geodesy and Remote Sensing, tel. (383)343-29-11, e-mail: btmazurov@mail.ru
Currently active is a practical implementation of the state geodetic coordinate system, 2011 on the territory of the Russian Federation. In this regard, entify a need clarification the mathematical and methodological foundation for the comparison of the parameters of the earth ellipsoid under transformations of coordinate systems, coordinate transformations, and solution of problems arising from the coordinate transformation from the local coordinate systems in national. Addition should be taken into consideration account the outcomes of the program by the construction of modern satellite State Geodetic Network Russian of three levels (TAGS, GHS and GHS-1), as well as the accuracy of its connection with geodetic networks of triangulation and polygonometry 1-4-th class. When construction of global geodetic networks associated with processing of satellite observations, and the solution of geodetic problems to determine the position of points on the earth surface are applied the spatial rectangular coordinates and geodetic coordinates . In this context arise the problem of transforming these systems of coordinates. The article proved mathematically sound high-precision non-iterative method of procedure for the direct calculation of geodetic latitude and altitude from the initial data - three-dimensional rectangular coordinates.
Key words: mathematics, astronomy, cartography, geodesy, geodynamics.
В работах [1, 2] показаны перспективы дальнейшего развития системы координат Российской Федерации 2011 г. на период до 2020 г. По принципам ориентировки в теле Земли ГСК-2011 идентична Международной земной опорной системе координат ITRS, установленной в соответствии с рекомендациями Международной службы вращения Земли (International Earth Rotation and Reference Systems Service - IERS) и Международной ассоциации геодезии (International Association of Geodesy - IAG).
Успешное решение этой государственной задачи предполагает соответствующее математическое и методологическое обоснование с учетом схемы построения государственной геодезической сети России трех уровней (ФАГС, ВГС и СГС-1) [3]. Остается также актуальной задача по преобразованию пространственных прямоугольных координат X, Y, Z и геодезических координат B, L, H.
Математическое обоснование алгоритмов преобразований имеет большую историю. Упомянем, в первую очередь работы B. R. Bowring, например, [4]. В настоящее время продолжаются научные исследования по данному вопросу многими другими учеными. Сравниваются вычислительные методы [5], предлагаются новые алгоритмы [6-8]. Традиционно, вычисление пространственных прямоугольных координат по заданным геодезическим координатам выполняют по соотношениям, как, например, в работе [9]:
' X = (N + H) cos B ■ cos L; Н = (N + H)cosB ■ sinL;
Z = ( N (1 - e2) + H) sin B,
где N = а л/1 - е2 зт2 В - радиус кривизны первого вертикала;
а - большая полуось;
е - первый эксцентриситет земного эллипсоида.
Выражения (1) относительно простые и позволяют легко определить искомые величины X, У, Z. На поверхности конкретного эллипсоида равенства (1) представляют собой функции, выраженные только через исходные данные, т. е. X = X (В, Ь, Н), У = У (В, Ь, Н), 2 = 2 (В, Ь, Н). Естественно, в такой же форме: В = В(Х, У, 2), Ь = Ь(Х, У, 2), Н = Н(X, У, 2) - желательно получить и алгоритмы для обратного перехода. Однако эта задача, решаемая обращением выражений (1), приводит к громоздким и сложным закономерностям. Из трех определяемых величин наиболее простой алгоритм предложен для вычисления геодезической долготы Ь в работе [10]:
1. L0 = 2 arctg (УЦХ + R)), R = VX 2 + Y 2 ф 0.
2. L =
IL0, если Y > 0, L0 + 2л, если Y < 0;
3. L =
f0, при X > 0, л, при X < 0, если Y = 0.
Приведенный выше алгоритм по сравнению с системой формул (1) содержит больший объем вычислительных операций и имеет более сложную структуру. Для вычисления геодезической высоты в публикациях [9, 11] и даже в Госстандарте РФ [12] рекомендуются формулы:
H = R/cos B - N; (2)
H = Z¡sin B - N(1 - e2), (3)
которые, как установлено в работе [13], при малой погрешности широты AB образуют большую погрешность AH при вычислениях высоты.
Более стабильные по точности результаты получаются при использовании формулы
H = RcosB + ZsinB - ayl 1 - e2sin2 B . (4)
В этом случае погрешность формулы (4), обусловленная неточной величиной широты, связана представленной в работе [10] зависимостью
AH = - 2 (a + H )-AB 2. (5)
Основной недостаток приведенных алгоритмов (3), (4) состоит в том, что они не образуют прямого пути вычисления высоты по исходным данным, так как для их применения нужно знать определяемую в процессе решения задачи геодезическую широту В.
Вместе с тем, в работах [14-16] предлагается совместное вычисление геодезической широты В и геодезической высоты Н путем решения системы уравнений способом последовательных приближений. Недостатки этого пути решения задачи приведены в [10, 17] и заключаются в исключении при исследованиях условий сходимости итерационных процессов, отсутствии анализа скоростей приближений к определяемым величинам B и H, оценки точности по итерациям и конечных результатов.
Наиболее часто геодезическую широту В определяют как, например, в работе [4] путем решения трансцендентного уравнения
tg B = (z + e 2N sin B)/r . (6)
Такие решения без должного анализа результатов выполняют методом последовательных приближений [9, 18], методом разложений в ряды и дифференциальных поправок и другими приближенными способами.
Точные алгоритмы решения уравнения (6) получают путем его преобразования с помощью тригонометрических тождеств к алгебраическому виду, как в работе [19] и др. Недостатком этих решений является громоздкость и сложность полученных алгоритмов и практическая непригодность их для теоретических исследований и решения смежных задач.
Исследования разных подходов и методов решения уравнения (6) выполнены в работах [10, 13] и основаны на разработанной теории отделения корня этого уравнения. Определение промежутка изоляции тангенса искомой широты позволило выполнить сравнительный анализ решений, предложенных как отечественными, так и зарубежными учеными, применить для вычисления широты методы хорд и касательных, как в [13].
Полученные результаты исследований в [13] были использованы для решения геодезических задач в пространстве. При геодезических высотах H = 0 решения задач выполняются на поверхности эллипсоида, теоретические основы которых были заложены еще Л. Эйлером [20].
Из итерационных методов наиболее быструю квадратичную сходимость обеспечивает решение уравнения методом касательных с приведенной широтой и [13]:
tg u = W1 - e 2 + ae 2sin3 u 0 = П1 (7)
R - ae cos u0 П
2
Такая закономерность позволяет на основе уравнения (7) построить неитеративный высокоточный алгоритм, выраженный только через исходные данные.
Для этого воспользуемся начальным приближением [13]:
zVi
tg uo ; Л0 4R2 + Z7(i - e2), (8)
R11 -ae /a0l
где a o - большая полуось вспомогательного эллипсоида, проходящего через точку M((, Y, Z) и имеющего тот же эксцентриситет и центр, что и земной эллипсоид.
С помощью зависимости (8) исключим содержащиеся в (7) величины sin u0, cos u0. Для удобства преобразований введем обозначения:
P = 1 - ae2/a0 ; Di = ^(RP)2 +(l - e2 )• Z2 ; Ri = R/Di; Zi = Z¡DX. (9)
Тогда tg u0 = zVi - e2 /(HP); cos u0 = PRi; sin u0 = Z^i - e2 и выражения формулы (7) преобразуются к виду
П2 = R - ae2 cos3 u0 = Ri (di - ae2P3R2);
Пi = Zx4i - e2 + ae2 sin3 u0 = Z^i - e2 (di + ae2 (i - e2)Z2). Полученное соотношение ni выразим через Ri на основе зависимости (9)
D2 = (PR)2 +(i - e 2) Z2.
2
Разделив обе части этого равенства на Di , получим
i = P 2 Ri2 +(i - e 2) Zi2.
В этом случае
ni = Z^i - e2 (d1 + ae2 (i - P2R2)).
С полученными результатами П1 и П 2 алгоритм (7) принимает вид
А zJí - e2 D1 + ae2 c0 - ae2 P 2 R12
tgu = ---1-0 2 2 3—L.
R1 D1 - ae 2 Ri2 P3
После деления числителя на знаменатель, с учетом зависимости
tg u ^ 1 - e2 tg B (10)
получаем формулу для вычисления тангенса геодезической широты
( л (л тЛ \2 Л
tg B_R
1 + 1 -(1 - P ).(PR )2
, Dj(ae2)- P • (PR )2
(11)
Погрешность формулы (11) определяется по выражению, предложенному в работе [8]
3 a3 H2 e10 7 3 AB" = 3р" • a H e 5 sin7 B cos3 B. (12)
2 (a + H )5
Исследованиями на экстремум функции (12) устанавливаем, что наибольшая
2
погрешность широты достигается при B = 56°47' и H = 3a « 4 250 км и составляет AB " = 6,8"-10 "9. Чтобы получить результат с такой точностью, нужно по формуле (11) вести вычисления на 16-разрядной сетке. Для точек земной поверхности при |h| < 10 км и B « 57°, максимальная погрешность AB" = 4,8"-10"13
будет определяться при вычислениях на 20-разрядной вычислительной сетке.
Выведем теперь формулу для прямого вычисления высоты Н путем исключения в равенстве (4) функций широты с помощью начального приближения (8) и введенных величин (9). В этом случае, с учетом соотношения (10), получаем:
, B Z RP RP RP
tg B =-; cos B = - - 1 •
VP2 R2 + Z2 VD2 + e2 Z2 VI+TZf
inB = zJ^1 + e2Z2 ; aV1 - e2 sin2 B _ a/+ e2Z2
и формула (4) принимает вид
_ D| (PR2 + Zj2)- a
H= ^ 1 1 . (13)
V1 + e 2 Z12
На основе зависимости (5) с учетом погрешности начального приближения (8), описанного в [13]
AB _ aHe42 sin3 BcosB, (a + H )2
для оценки точности формулы (13) получаем выражение
АН = —•
1 а 2Н 2 е8
2 (а + Н )3
Бт3 В соб В
Г
(14)
Методами дифференциального исчисления устанавливаем, что наибольшая погрешность формулы (13) получается при Н = 2а и В = 60° и составляет |АН| = 0,1 мм при вычислениях с 13-значащими цифрами.
Для точек земной поверхности при |Н < 10 км ее погрешность не превосходит |АН| < 2 • 10 6 мм, что обеспечивает высокоточный результат по сравнению с полученными в [10, 16] алгоритмами. Даже при Н = 1000 км АН = 0,01 мм.
Для предвычисления точности результатов преобразования координат и определения количества значащих цифр, с которыми нужно выполнять решения, выразим погрешности (12), (14) через исходные данные, полагая приближенно:
Н = а 0 - а; а/ (а + Н )= а/а 0 ; Н/ (а + Н ) = 1 -
а а,
0 >
Бт В = 2/а0 ; соб В = Я/
а
Тогда
АВ' = 3 (• е10 )•
í \3 í а
V а0 У
1 -
а
а
í ^ \
0
2
V а0 у
V а0 у
(15)
АН
ае
а0
4^0 у
а
\2
а0
V
г 2 V
V а0 у
Я а0
2
м
(16)
При оценке погрешностей по выражениям (12), (14)-(16) достаточно вести вычисления с тремя верными значащими цифрами. Исключение составляет определение высоты точек земной поверхности по формуле Н = а0 - а. В этом случае ввиду близости значений а0 и а происходит потеря точности. Поэтому вычисление а0 нужно выполнить по формуле (8) с большим числом десятичных разрядов.
Итак, по результатам исследований нами предлагаются следующие высокоточные неитеративные алгоритмы для прямого вычисления геодезической широты В и высоты Н по исходным данным X, У, 2.
1. Определяем постоянные величины:
¿0 = 1 - е2 = (( -1)//)2; кх = ае 2 = а (2 / -1V/ 2 ,
где / - знаменатель сжатия а = ( а -Ь)/а = 1//.
На эллипсоиде Красовского: а = 6 378 245 м; / = 298,3.
2
7
3
Для системы координат ГСК-2011 [1] эти параметры: а = 6 378136,5 м; / = 298,256 4151.
2. Вычисляем вспомогательные величины:
а) R = Vx2 + 7 2 ; б) ^ = W 1 + (Z/RT/k0 ;
в) P = 1 - kila0 ; г) D1 =Ryjp2 + k0(Z/R)
д) Zx = ; R1 = R/A .
3. Определяем геодезические широту В и высоту Н:
( Л и тЛ \2 Л
B = arctg
Z R
1 + 1 -(1 - P )-(PRj )2
Д/kx -P-(PR!)2
= Dj (PR2 + Zf)- a
H= ,-
V1 +e 2 Z12
Таким образом, здесь математически обоснованы высокоточные неитеративные алгоритмы для прямого вычисления геодезической широты и высоты по исходным данным - пространственным прямоугольным координатам.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Горобец В. П., Ефимов Г. Н., Столяров И. А. Опыт Российской Федерации по установлению государственной системы координат 2011 года // Вестник СГУГиТ. - 2015. -Вып. 2 (30). - С. 24-37.
2. Голякова Ю. Е., Касаткин Ю. В., Щукина В. Н. Анализ установления единых государственных систем координат // Вестник СГУГиТ. - 2015. - Вып. 2 (30). - С. 55-61.
3. Анализ состояния государственной геодезической сети России с учетом существующих и перспективных требований / Е. М. Мазурова, К. М. Антонович, Е. К. Лагутина, Л. А. Липатников // Вестник СГУГиТ. - 2014. - Вып. 3 (27). - С. 84-89.
4. Bowring B. R. The accuracy of geodetic latitude and height equations. Surv. Rev., 28, 202-206, 1985.
5. Burtch R. A comparison of methods used in rectangular to geodetic coordinate transformations. Proceedings of the American Congress on Surveying and Mapping (ACSM) Annual Conference and Technology Exhibition, Orlando, April 21-26, 2006.
6. Pick M. Closed formulas for the transformation of the Cartesian coordinate system into a system of geodetic coordinates. Stud. Geophys. Geod., 29, 112-119, 1985.
7. Sofair I. Improved method for calculating exact geodetic latitude and altitude revisited. J. Guid. Control Dyn., 23, 369-369, 2000.
8. Vermeille H. Computing geodetic coordinates from geocentric coordinates. J. Geodesy, 78, 94-95, doi:10.1007/s00190-004-0375-4, 2004.
9. Морозов В. П. Курс сфероидической геодезии. - М. : Недра, 1979. - 296 с.
10. Медведев П. А. Анализ преобразований пространственных координат точек земной поверхности // Геодезия и картография. - 2014. - № 4. - С. 2-8.
11. Ганьшин В. Н. Переход от пространственных прямоугольных координат к геодезическим // Геодезия и картография. - 1970. - № 1. - С. 17-19.
12. ГОСТ Р 52572-2006. Географические информационные системы. Координатная основа. Общие требования. - М. : Стандартинформ, 2006. - 19 с.
13. Медведев П. А. Анализ преобразований пространственных прямоугольных координат в геодезические. - Омск : Изд-во ОмГАУ, 2000. - 104 с.
14. Огородова Л. В. Совместное вычисление геодезической широты и высоты точек поверхности Земли // Геодезия и картография. - 2011. - № 9. - С. 11-15.
15. Огородова Л. В. Простой и надежный способ вычисления геодезической широты и высоты точек поверхности Земли по прямоугольным координатам // Изв. вузов. Геодезия и аэрофотосъемка. - 2014. - № 2. - С. 63-66.
16. Алгоритм вычисления геодезической высоты по пространственным прямоугольным координатам / В. Н. Баландин, М. Я. Брынь, С. П. Имшенецкий, А. Ю. Матвеев, А. В. Юськевич // Геодезия и картография. - 2006. - № 6. - С. 15-16.
17. Медведев П. А. Определение погрешностей геодезической высоты, широты и долготы аналитическими методами // Геодезия и картография. - 2009. - № 1. - С. 25-27.
18. ГОСТ Р 51794-2008. Глобальные навигационные спутниковые системы. Системы координат. Методы преобразований координат определяемых точек. - М. : Изд-во стандартов, 2008. - 16 с.
19. Пенев П., Пенева Е. Преобразование прямоугольных геоцентрических координат в геодезические без применения итераций // Изв. вузов. Геодезия и аэрофотосъемка. - 2012. -№ 3. - С. 34-38.
20. Мазуров Б. Т., Медведев П. А. Леонард Эйлер - вклад для астрономии, небесной механики, геодезии, картографии, геодинамики // Интерэкспо ГЕО-Сибирь-2014. Х Между-нар. науч. конгр. : Междунар. науч. конф. «Геодезия, геоинформатика, картография, маркшейдерия» : сб. материалов в 2 т. (Новосибирск, 8-18 апреля 2014 г.). - Новосибирск : СГГА, 2014. Т. 1. - С. 186-191.
Получено 10.04.2016
© П. А. Медведев, Б. Т. Мазуров, 2016
