CC BY
78
УДК 528.235
С. Н. Голубков, О. А. Павлова, Е. А. Паниди, В. М. Щербаков
автоматизированная система для анализа основных метрических свойств картографического изображения
Рассмотренный в данной работе программный комплекс (ПК) является следующей версией ранее созданного пакета утилит генерации номенклатуры российских топографических карт, вычисления координат углов трапеций карт и перевычисления координат точек из плоских прямоугольных систем координат в геодезические и обратно [1].
Обновленный ПК, который носит название «Топографический калькулятор» (“Topographic Calculator”), расширен новыми функциями: вычисление площади сфероидической трапеции, расчет длины дуги локсодромии и расчет длины дуги ортодромии.
Актуальность исследования в данном направлении выражается в том, что разработчики большинства используемых в настоящее время геоинформационных систем и другого специализированного программного обеспечения, как правило, не публикуют алгоритмы и формулы, использованные ими в своих программах, вследствие чего, оценить корректность и точность вычислений становится весьма проблематично.
Кроме того, следует отметить, что реализация приведенных вычислений во многих программных продуктах зачастую неудобна для пользователя. В качестве примеров можно привести широко известные в нашей стране ГИС MapInfo и ArcGIS. Так в MapInfo все вычисления производятся на плоскости или сфере, но не на эллипсоиде, что не применимо в случае, когда пользователь располагает достаточно точными данными. В ArcGIS вычисление расстояний между двумя точками, протяженности линейных объектов и площадей полигонов производится в единицах измерения используемой системы координат, то есть если данные представлены в геодезической системе координат и координаты измеряются в градусах, то расстояния, длины и площади, как ни странно, тоже будут вычислены в градусах, в случае, когда для вычисления используются стандартные приемы. Таким образом, с точки зрения корректного использования метрики картографического изображения геометрические характеристики картографируемых объектов можно вычислить только по данным, представленным в какой-либо проекции, учитывая, что в этом случае они будут искажены за счет свойств проекции.
Алгоритмы использованные в «Топографическом калькуляторе» разработаны для производства вычислений на поверхности эллипсоида. Так как эллипсоид вращения по своей форме более близок к форме геоида чем сфера, то в общем случае автоматизированные вычисления расстояний, длин и площадей удобнее производить именно на поверхности эллипсоида. С одной стороны это позволяет достичь более высокой точности самих вычислений, с другой — дает возможность сделать алгоритм расчетов более универсальным. Несмотря на то, что вычисления на поверхности эллипсоида более трудоемки по сравнению с вычислениями на поверхности сферы, этот недостаток успешно компенсируется скоростью расчетов на современной вычислительной технике.
© С. Н. Голубков, О. А. Павлова, Е. А. Паниди, В. М. Щербаков, 2008
Алгоритм вычисления площади сфероидической трапеции
Площадь бесконечно малой сфероидической трапеции вычисляется следующим образом [3]:
dS = MN cos ф dф d 1,
a(1 - e2) a 2 a2 - b2
гдеM = , 2 ; N = , 2 ; e2 =-------- —; 1 — долгота; ф — широта;M — радиус
V(1 - e2sm^)J V1 - e2sm^ a2
кривизны меридиана; N — радиус кривизны первого вертикала; e — первый эксцентриситет эллипсоида; a — большая полуось эллипсоида; b — малая полуось эллипсоида
Соответственно формула для расчета площади трапеции на поверхности эллипсоида, ограниченной широтами ф; и ф2 и долготами 11 и 12, будет представлять собой интеграл следующего вида:
12 ф-2 ф_2 sin ф2
S = j jMNcosф dipdX = (12 - 11)¿2 j С(7.^ )2 = (12 - ^!)b2 I (1 ds2iricp2
і І І (1 - e2sin2cp)2 J (1 - e2sin2p)2
Xi фі ф
Алгебраическое решение данного уравнения, выглядит следующим образом.
sin ф2
Примем: Z = sin ф. Тогда S = (1 - l^b2 |
dZ
2у2\2 '
sin фі
(1 - e2z2)
г dz _ г г dz г г dz г г dz г г dz
1 (1 - e2z2)2 = 4 1 (1 - ez)2 + 4 1 1 - ez + 4 1 1 + ez + 4 J (1 + ez)
sin фі sin ф1 sin фі sin фі sin фі
sin ф2 /1 11 + (— ln(1 + ez) - ■
sin ф1
4e 4e 1 + ez
1 1 ln(1 - e sin ф2) ln(1 - e sin ф1)
4e(1 - e sin ф2) 4e(1 - e sin ф1) 4e 4e
ln(1 + e sin ф2) ln(1 + e sin ф1) 11
4e 4e 4e(1 + e sin ф2) 4e(1 + e sin ф1)
1 ( (1 + e sin ф2)(1 - e sin ф1^ 2e sin ф2 2e sin ф1
4e V (1 + e sin ф!)(1 - e sin ф2) (1 - e sin ф2)(1 + e sin ф2) (1 - e sin ф j)(1 + e sin ф^
Произведем следующие замены:
q = e sin фь p = e sin ф2, u = 1 - q, v = 1 + q, i = 1 -p, j = 1 + p.
sin 2
Ы11 cp 1
В результате, с учетом введенных обозначений, формула для расчета площади трапеции на поверхности эллипсоида преобразуется в следующее выражение:
^ = (12 - ЮЪ2 Л ] • и х 2р 2д
4е ^ 1 • у 1 • ] и • V
Более подробно решение представлено в виде:
Г аС = Г ос + 1 Г ос + 1 Г ос + _1. Г ос
Г (1 - е2С 2)2 4 Г (1 - еС)2 4 Г 1 - еС 4 Г1 + еС 4 Г (1 + еС)2’
и в соответствии с одним из стандартных приемов интегрирования, преобразуем подин-тегральную функцию
ос = -1 г а(-еС) = -1 |-а(1 - еС)=-1 а=-^, =
(1 - е2С 2)2 е Г (1 - еС)2 е Г (1 - еС)2 е Г ?2 е Г
1 ? 1 1 1
е -1 е? е(1 - еС)
где ? = 1 - еС.
Принимая во внимание равенства О(-ах) = -аОх и О(х + 1) = Ох, получаем:
Г
Х'Ох = ----------;
П + 1
г ОС 1 |-О(1 - еС) 1п(1 - еС)
Г 1 - еС еГ 1 - еС е
Используя равенство |йХ/х = 1п х, приходим к следующим выражениям:
г ОС _ 1г О(еС) _ 1 г О(1 + еС) _ 1п(1 + еС)
Г 1 + еС е Г 1 + еС еГ 1 + еС е
г ОС 1 Г О(еС) 1 Г О(1 + еС) 1 1
Г (1 + еС)2 е Г (1 + еС)2 е Г
(1 + еС)2 ^ (1 + еС)2 ^ (1 + еС)2 е 1 + еС
Алгоритм вычисления длины дуги локсодромии
Формула расчета длины дуги локсодромии для эллипсоида имеет вид [2]:
¿2 - ^
cos а
где: S1 — длина дуги меридиана от экватора до начальной точки дуги локсодромии, S2 — длина дуги меридиана от экватора до конечной точки дуги локсодромии, a — угол (азимут) под которым локсодромия пересекает меридианы
Расчет величины азимута локсодромии производится по следующей формуле [2]:
, ( х2 - X,'
a ' 310,8 bwT
где Xj, Y, — координаты начальной точки дуги локсодромии в проекции Меркатора; X2, Y2 — координаты конечной точки дуги локсодромии в проекции Меркатора.
Формулы нормальной равноугольной цилиндрической проекции Меркатора представляются следующими равенствами [2]:
r Io
X = , Y = r0 • ln(U):
a • oos cp 0 180 tg(45° + ф)
где r0 = , , p '-------------, U =-г-r-, y ' arosin (e • sin p).
V1 - e • sinp p ,ge ( 45° + W
Здесь I — долгота точки; ф — широта точки; r0 — радиус стандартной параллели проекции; e — первый эксцентриситет эллипсоида; a — большая полуось эллипсоида; b — малая полуось эллипсоида; p ' 3,1415926536; po = 57,2957795131 (радиан выраженный в градусах); U — изометрическая широта точки.
Формула расчета длины дуги меридиана от экватора до заданной параллели представляется рядом [4]:
S = a • ф • ( 1 - n + — (n2 - n3) + —(n4 - n5) + ...
3 ( 7 , л 55
- — a • sin (2p) • I n - n + — (n3 - n4) + — n + ..
15 3
+ — a • sin (4p) • I n2 - n3 + — (n4 - n5) + ...
35 ( , „ 1K > 315 ( ,
- —- a • sin (6p) • n3 - n + — n + ... + —— a • sin (8p) • n4 - n + ...
48 I 16 j 512 I
где n '--------; a — большая полуось эллипсоида; b — малая полуось эллипсоида;
a + b
p — широта данной параллели.
Алгоритм вычисления длины дуги ортодромии
В основу алгоритма вычисления длины дуги ортодромии был положен способ геометрического построения ортодромии на поверхности эллипсоида. Этот способ построения заключается в следующем. Точки A и B на поверхности эллипсоида, между которыми необходимо построить дугу ортодромии соединяются хордой. Из центра хорды
к поверхности эллипсоида опускается нормаль. Точка пересечения построенной нормали и поверхности эллипсоида (С), также будет лежать и на дуге ортодромии соединяющей А и В. Далее строятся хорды АС и ВС. Дальнейшие построения продолжаются таким же образом. При достижении достаточно малых величин выстраиваемых хорд, разницей между суммой длин хорд, соединяющих исходные точки А и В поверхности эллипсоида и дугой ортодромии, соединяющей эти же точки, можно пренебречь. Тогда, длину дуги ортодромии между двумя исходными точками можно представить как сумму длин всех построенных малых хорд с заданной вычислительной точностью.
В «Топографическом калькуляторе» описанное построение производится аналитически в прямоугольной системе координат с началом координат в центре эллипсоида — (X, У, 7). В этой же системе вычисляются и длины хорд.
Для задания исходных пунктов используются геодезические координаты (долгота, широта, высота над поверхностью эллипсоида). Переход от геодезических координат к прямоугольным и обратно осуществляется по нижеприведенным формулам [5].
Переход от геодезических координат к прямоугольным:
Широта вычисляется методом итераций. Здесь N — радиус кривизны первого вертикала в данной точке; е — первый эксцентриситет эллипсоида; а — большая полуось эллипсоида; Ь — малая полуось эллипсоида.
Вычисление длины отрезка (хорды) производится по формуле:
где: Х1, У1, 71, Х2, У2, 72 — координаты, соответственно первой и второй точек.
Для вычисления координат точки пересечения поверхности эллипсоида, с нормалью опущенной к ней из центра хорды, сначала необходимо вычислить прямоугольные координаты центра хорды:
X \ ((Ы + И) • cos (ф) • cos (1) ^
У = (Ы + И) • cos (ф) • sin (1)
Ч 2 ) \ ([1 - е2] • Ы + Ъ) • ^п (ф) у
Переход от прямоугольных к геодезическим координатам:
Ъ = X • sec (1) • sec (ф) - N.
ь = -X - X,)2 + (У - У,)2 + 2 - г,)2,
(X ^ ( (Х2 - X,) / 2 Л
У = (У - У,) / 2 ^ г J ^ г - г,) / 2 J
Затем от прямоугольных координат требуется перейти к геодезическим. В результате получаем геодезические координаты центра хорды, в которых высота будет отрицательна, так как центр хорды всегда будет лежать в теле эллипсоида, то есть ниже его поверхности. Если же в полученных геодезических координатах центра хорды высоту приравнять нулю, то получим координаты точки пересечения нормали к поверхности эллипсоида (так как высота в данном случае откладывается по нормали), проходящей через центр хорды, с поверхностью эллипсоида. После этого полученные геодезические координаты нового узла дуги ортодромии можно пересчитать в прямоугольные.
Последовательно вычисляя координаты новых узлов дуги ортодромии, и сумму длин соединяющих их хорд, можно с необходимой точностью вычислить длину дуги ортодромии. Для этого, требуется постоянно контролировать разность текущего результата вычисления длины с предыдущим результатом. И если разность эта не будет превышать заданного значения, то процесс вычисления можно остановить.
Подводя итог, еще раз отметим, что картометрические вычисления на поверхности эллипсоида позволяют минимизировать искажения результата вычислений, особенно в сравнении с вычислениями в плоскостях различных картографических проекций. Подобный подход может быть успешно реализован в среде современных ГИС. Однако, на сегодняшний день, ни в одной из широко распространенных ГИС не существует аппарата предназначенного для производства широкого спектра картометрических работ на поверхности Земного эллипсоида. Таким образом, очевидна актуальность развития вычислительных процедур и методик вычисления других картометрических характеристик на поверхности эллипсоида.
Summary
Golubkov S. N., Pavlova O. A., Panidi E. A., Shcherbakov V. M. Automated system for analysis of cartographic imagery basic metrical characteristics.
An updated version of the computer program “Topographic calculator” created at the Cartography Department of St.-Petersburg state university is considered. The chartometric work technique on a terrestrial ellipsoid surface is considered. Algorithms for calculation of a great circle and rhumb line lengths on ellipsoid surface and an algorithm for calculation of a ellipsoidal trapeze area are analyzed. All mentioned algorithms are realized in “Topographic calculator”.
Key words: cartometry, spheroid trapezium, geodesic line, loxodrome.
Литература
1. Афанасьев В. А., Павлова О. А., Паниди Е. А, Щербаков В. М. Автоматизированная система генерации номенклатуры российских топографических карт и пересчёта координат между различными системами // Сб. статей. Теория и практика эколого-географических исследований. СПб., 2005.
2. Бугаевский Л. М. Математическая картография. М., 1998.
3. Закатов П. С. Курс высшей геодезии. М., 1976.
4. The universal grids: Universal Transverse Mercator (UTM) and Universal Polar Stereographic (UPS) — Fairfax: Defense Mapping Agency, 1989.
5. World Geodetic System — 1984 (WGS-84) Manual — Montreal: International Civil Aviation Organization, 2002.
