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



Математическое моделирование: методы, алгоритмы, технологии

УДК 519.6:527

В.А. Ботнев, С.М. Устинов

методы решения прямой и обратной геодезических задач

с высокой точностью

V.A. Botnev, S.M. Ustinov

methods for direct and inverse geodesic problems solving

with high precision

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

МЕТОДЫ МОДЕЛИРОВАНИЯ; ПОТЕРЯ ТОЧНОСТИ; МАШИННОЕ ЭПСИЛОН; АСИМПТОТИЧЕСКОЕ РАЗЛОЖЕНИЕ; ГЕОДЕЗИЧЕСКИЕ ЗАДАЧИ; ЛОКСОДРОМИЯ.

Proposed algorithms for navigational calculations provide an extremely low error over the entire range of input parameters and allow avoiding the loss of accuracy typical for these tasks. These methods demonstrate an acceptable performance for practical applications and we can recommend using them in commercial navigational software products meeting international navigational standards.

MODELLING METHODS; LOSS OF ACCURACY; MACHINE EPSILON; ASIMPTOTIC EXPANSIONS; GEODESIC PROBLEMS; LOXODROME.

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

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

Методы решения обоих вариантов задач хорошо известны и на их основе реализованы многочисленные программные продукты, успешно используемые на практике. Первый вариант задачи (ортодромия)

4

Локсодрома на эллипсоиде

глубже исследован, а соответствующее программное обеспечение хорошо апробировано и имеется в открытом доступе [3].

Методики решения задач во второй постановке (локсодромии) отличаются большим разнообразием [4—7] и среди них нет единой общепризнанной и всесторонне проверенной на практике.

Так, например, в работах [5, 7] указывается на то, что существует некоторое многообразие различных методов и формул, часть из которых появилась относительно недавно. В [5] отмечается также, что в попытке создать простые и эффективные алгоритмы некоторые авторы предлагают методы расчета с гигантскими ошибками. С другой стороны, даже подзадачи являются предметом отдельных статей [8, 9], и для них также не существует единого метода решения. Так, например, в работе [8] разбираются несколько приемов для нахождения меридионального расстояния (эта подзадача является одним из составных элементов прямой и обратной задач). Авторы статьи [8] не останавливают свой выбор на каком-то одном из рассматриваемых ими методов, а только лишь забраковывают некоторые из них.

В этих и в ряде других публикаций указывается на целый ряд проблем, возника-

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

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

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

используемые алгоритмы ориентированы на поиск различного рода экстремальных точек, требования к точности решений геодезических задач ужесточаются. Необходимая погрешность в них не должна превышать 10-15. Для арифметики с двойной точностью, где машинное эпсилон порядка 10-16, это достижимо, но труднореализуемо. Учетверенная же точность в настоящее время доступна на аппаратном уровне только на весьма ограниченном круге платформ, а надежные кроссплатформенные библиотеки для работы с числами с плавающей точкой учетверенной точности отсутствуют. К тому же из-за программной реализации скорость работы с числами учетверенной точности не очень велика, в то время как предлагаемый в работе алгоритм для чисел двойной точности в несколько раз опережает по скорости максимально упрощенный вариант алгоритма с числами учетверенной точности.

Кроме сказанного выше имеется целый ряд широко распространенных устройств (планшетные компьютеры, встраиваемые контроллеры и пр.) с величиной машинного эпсилон порядка 10-7. При этом существуют многочисленные приложения, в которых такая точность допустима. Предлагаемые ниже алгоритмы актуальны и в этом случае. Они обеспечивают погрешность, немногим превышающую машинное эпсилон, независимо от того, какая точность (е « 10-7или е « 10-16) обеспечивается на аппаратном уровне. В случае е « 10-7 учитывается меньшее число слагаемых соответствующих разложений.

Методика решения прямой геодезической задачи

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

Заданы точка с широтой ф1 и долготой А1? азимут а и расстояние Б до второй точки. Требуется найти широту ф2и долготу А2второй точки.

Эта геодезическая задача для локсодромии может быть разделена на следующие этапы:

1. Вычисление меридионального расстояния М1 для первой точки (см. рисунок).

2. Вычисление разности меридиональных расстояний ДМ и географической широты ф2 второй точки.

3. Вычисление разности изометрических широт.

4. Вычисление разности долгот второй и первой точек ДА,.

Вычисление меридионального расстояния для первой точки. Меридиональное расстояние М1 первой точки рассчитывается по хорошо известной формуле:

ф1 а(1 - е2)аф

М1 = |

(1 - е2 зт2 ф)3/2'

(1)

Известные методы вычисления интеграла (1) условно разделяются на три группы.

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

2) Разложение подынтегрального выражения в степенной ряд по различным параметрам с последующим интегрированием этого ряда. Простейший вариант — выбрать в качестве малого параметра е2 = /(2 - /) (где / — это коэффициент сжатия земного эллипсоида), которое действительно мало и для модели Земли, принятой для международного мореплавания ^0884), составляет около 0,0067. Но существует

и более удобная для этой цели величина:

/

т. н. третий эксцентриситет п = ^ - , который еще меньше и приблизительно равен 0,0017. Важно отметить, что выбор этого параметра сокращает число слагаемых ряда и не приводит к усложнению вычисления коэффициентов. Как следствие, это позволяет уменьшить время расчета и снизить накапливание ошибок. Используемое разложение в ряд имеет вид:

М(ф) = а(1 - п)(1 - п2)5(ф),

(2)

где

V

S(ф) « aoф + £ (-1) a2k з1л(2кф), p = 6 (3)

1 , 9 ,„2 . 4 . 1225 <6

а„ = 1 +--П + 1

225 4 1225

-п +-п6 + ...

64 256

3 45 3 525

п +

а = — п + 2 2 16

п + ...

15

а. = — 4 16

п2 +

128

105 4 4725 6

-п +-п + ...

64 2048

а =

35 3 315

п +

48

256 2079

п + ...

315 4 20/9 6

а„ =-п +--п + ...

8 512 2048

693 1280

п +...

1001 6

а12 =-п + ...

12 2048

Разложение этого вида было впервые предложено Гельмертом с точностью до а8 [10] и уточнено до a10 в работе [7]. Для достижения необходимой точности нужно было продолжить представленный ряд до a12 (р = 6 в формуле (3)) и вычислить необходимые коэффициенты, содержащие слагаемые с П .

Введение a12 влияет на точность результата не слишком сильно, в то время кфк учет слагаемых си6в коэффициентах a0, a4 и а8 крайне важен и позволяет ограничить максимальную ошибку сверху величиной 2е. В противном случае она достигает 4е.

Суммирование ряда (2) осуществляется по технике Кленшоу [11], что позволяет сократить объем вычислений и получить более точный результат. На практике эта процедура сводится к следующему. Зададим

Ср+ 1(ф) = СР+ 2 (ф) = Тогда, вычисляя последовательно по мере уменьшения индексов все ср (ф) по формуле

ск (ф) = а2к + 2ск+1(ф) СОЗ(2ф) к ск+2(ф), 1 < к < р, получим с1(ф) и

S (ф) = а0 ф + с1(ф) Б1л(2ф). (5)

Как показано в работе [11], формула

(4)

(5) дает для S(ф) тот же результат, что и

(3). Эта техника подобна схеме Горнера и позволят свести вычисление суммы к вычислению только двух тригонометрических функций (зт(2ф) и соэ(2ф)). Таким образом, расчеты значительно ускоряются. Кроме того, т. к. первое слагаемое в формуле

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

3) Представление интеграла (1) в виде комбинации легко вычислимых эллиптических интегралов. Исходный интеграл (1) представляется как комбинация интегралов в форме Карлсона [12, 13]:

Г

М (ф) = а 8т(ф)

Я (соз2 ф,1 - е2 81л2 ф, 1) -

е2 81л2 фЯв (сО82 ф,1 - е2 81л2 ф, 1) 3

2

СО8(ф)

+ (6)

д/1 - е2 81

81Л2 ф

где Яф и — базовые эллиптические интегралы в форме Карлсона, которые могут быть вычислены эффективно по специальным алгоритмам. Погрешность данного метода чуть больше 2е. Этот результат оказывается чуть менее точным, чем разложение в ряд с использованием техники Кленшоу.

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

Вычисление разности меридиональных

к

к=1

расстояний ДМ и географической широты ф2 второй точки. Нужно найти такое ф2, что выполнено равенство:

ф2 а(1 - е2)аф

М = J

(1 - в2 sin2 Ф)3/2'

(7)

где М 2 можно получить из разности изометрических расстояний ДМ = Б соБ(а), М2 = М1 + ДМ. Для определения ф2 = ф1 + Дф найдем сначала разность широт Дф. Использование в дальнейших расчетах величины Дф вместо ф2 вызвано тем, что при близких значениях ф1 и ф2 происходит катастрофическая потеря точности в операторе ф2 - ф1. С этой целью обратим ряд (3) и вместо М(ф) построим другой ряд ф(М), воспользовавшись техникой, изложенной в [7]. При этом, так же как и в предыдущем подразделе, необходимо в коэффициентах разложения учитывать слагаемые с п5 и п6, чтобы не использовать метод Ньютона, как это было предложено в [7].

Первоначально строится ряд:

ф(М) = 5 + ¿^1п25 +

+ b4 sin 4S + b6 sin 6S + b8 sin 8S + + b10 sin 10S + b12 sin 12S + ...

(8)

где

- 3 27 з

b = — n--n +

32

269 512

n5 +.

21

, _ 2 55 4 6759

b = — n--n +--

16 32 4096

151 3 417 5

b6 =-n--n + .

6 96 128

n6 + .

b =

1097 4 15543 6 n--n +

512

2560

8011 5

b10 = —— n - . 10 2560

293393 6

b12 = -----n - .

S =

61440 1

a(1 - n)(l - n2 ) M

1 + 9 n2 +

225 4 1225 n +

64

256

n6 +.

Если не учитывать коэффициент Ь12, то погрешность может достигать семи машин-

ных эпсилон, в то время как c его учетом мы получаем ошибку не более 2е. Далее,

Дф = ф(М + ДМ) - ф(М) = = AS + 2b sin ДS • cos| S + ^ +

AS 2

AS 2

+ 2b sin 4 AS • cos41 S + —

+ 2b sin 2AS • cos 21 S +

+ 2b sin 3AS • cos31 S + — I + (9)

+ 2b„ sin 5AS • cos 51 S + I +

+ 2b-, sin 6AS • cos 61 S +

Суммирование этого ряда также осуществляется по технике Кленшоу [11]. Следует при этом отметить, что относительная погрешность для Дф не превышает 2е.

После вычисления Дф легко определить ф2 = ф1 + Дф. В итоге ошибка по широте не превышает 15 нм в терминах расстояния вдоль меридиана.

Вычисление разности изометрических широт. Изометрическая широта 0 является функцией обычной широты:

Q = ln

^ 1 - в sin ф 4 2 Jl 1 + в sin ф

(10)

В дальнейших расчетах потребуется значение разности изометрических широт

Д 0 = 02 - 01.

Чтобы не потерять точность, необходимо избегать непосредственного вычитания выражений, полученных в результате промежуточных вычислений [14]. С этой целью вместо ф2 будем использовать выражение ф1 + Дф :

AQ = ln

tan | П + ф1 +Дф

tan | П + ф 1 4 2

(11)

X

+ ® 1п

(1 - е эт( ф1 + Лф))(1 + е эт ф:) (1 + е эт(ф1 + Лф))(1 - е эт ф1)

.(11)

При малых значениях Лф числители и знаменатели обоих логарифмов в (11) близки друг к другу, и здесь также может происходить потеря точности при вычитании близких чисел. Поэтому перспективным представляется разложение этих логарифмов в ряд как по самому параметру Лф, так и какому-либо другому малому параметру, зависящему от Лф.

Наибольшую проблему представляет первое слагаемое в (11), т. к. выражение под логарифмом может быть сколь угодно близким к нулю и сколь угодно большим в допустимом диапазоне изменения широт. Задача вычисления второго логарифма заметно проще, т. к. его аргументом является выражение, отстоящее от единицы на величину порядка 2е или 0,16 для земного эллипсоида. Рассмотрим четыре варианта решения возникших трудностей.

Вариант 1. Воспользуемся известным разложением для функции 1п(1 + х) при малом значении х:

X 2 X 3

1п(1 + X) = X--+--...

4 ' 2 3

(12)

Предварительно преобразуем первое

слагаемое в (11): Л(ф1, Лф) = 1п

1ап| П + ф1 +Лф ' 4 2

1ап| П + ф

= 1п

1ап1 П + ф11 + 1ап (Лф

1ап| П+ф

1 - 1ап| П + ьш (Лф

= 1п

1 +

1 - а • %

(13)

где а = Ьап | — + ^ 1 и % = Ьап ( — I. С уче-

.4 2 ) у 2

том диапазона изменения географической

п п

широты ф1 е | - — , — | параметр а не может

быть отрицательным.

Теперь, руководствуясь (12), получаем ряд по степеням %. Ему свойственна относительно медленная сходимость. Наиболее эффективный расчет будет получаться только в окрестности экватора, где а = 1, что приведет к сокращению всех четных степеней % в разложении:

'1 + % и ( z2 %4 % —■ 1 = 2% 1 + — + — + — 1 - ^ I 3 5 7 ^

Для упрощения дальнейшего анализа заметим, что

Д(ф1, Лф) = Л(-ф1, -Лф). (14)

Это позволяет ограничиться рассмотрением только положительных значений ф1? что гарантирует неравенство а > 1. Так как

1п

в этом случае

< \а%\, наибольшие про-

блемы вызывает разложение знаменателя дроби под логарифмом в правой части (13). Численные эксперименты показали, что использовать разложение в ряд выгоднее из соображений погрешности только при |а%| < 0,5. В противном случае без потери точности можно использовать непосредственное вычисление логарифма.

Этот критерий стал определяющим для стратегии вычисления Л(ф1, Лф). При граничном значении |а%| = 0,5 для вычисления Л(ф1, Лф) с максимальной точностью в 3е необходимо учесть 48 слагаемых разложения (13), что не очень эффективно. Реально эта погрешность будет несколько выше, т. к. имеется погрешность в вычислении Лф. Второе слагаемое в (11) во всех случаях разлагалось в ряд (12), т. к. его сходимость не представляла проблемы при любых значениях входных аргументов. В итоге погрешность вычисления только Л0 (считая аргументы абсолютно точными) достигала 6е.

Вариант 2. Второй метод заключается в использовании формулы:

х • 1п(1 + х)

1п(1 + х) =

(1 + х) -1

(15)

применимость которой оправдана при |х| < 0,5 [15]. При этом если сложение 1 и х дает 1 (х меньше машинного эпсилон), то достаточно в качестве ответа вернуть само

а

значение х (так как ln(l + х) « х).

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

Вариант 3. В широко распространенной библиотеке Boost Math [15] для вычисления 1п(1 + х) при малом значении х, если компилятор не позволяет реализовать подход, описанный в варианте 2, предлагается использовать Паде-аппроксимацию. Многочисленные эксперименты с этим алгоритмом показали, что погрешность достигает 6s, что даже хуже варианта 1.

Вариант 4. Наилучшие показатели обеспечил следующий прием:

'1 + y y

ln X = ln

1 - у

= 2 у

f у2 у4 у6 1 + — + — + — + ,

Л

(16)

где у =

X -1

и X > 0. Погрешность этого

X + 1

метода не превышает 2е, а проигрыш по времени варианту 2 (самому быстрому из всех методов) составляет около 14 %.

Существенная потеря точности обычного логарифма, как уже отмечалось, происходит в диапазоне х е [0,5; 1,5], что для переменной у отвечает диапазону 1 1" 3; 5 ^

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

у е

Таким образом, наблюдает-

Реализация алгоритма сводится к следующему. Чтобы для первого слагаемого в (11) обеспечить вид

ln

tan I П + ф1 +Аф 4 2

tan I П + ф 4 2

= ln

1 + у у

1 - у,

sin

Аф

следует выбрать у =

cos I ф1 +

Аф

Анало-

гично, второе слагаемое в (11) приобретает вид:

ln

(1 - е sin(ф1 + Аф))(1 + е sin ф11 (1 + е sin(ф1 + Аф))(1 - е sin ф11

1 + у

= ln

1-у

если у =

т . Аф f Аф -2е sin —-í- cos I ф, + —-

2 Г1 2

1 - е2 sin ф1 sin I ф1 +

Аф

Модуль второго слагаемого в (11) почти на два порядка меньше, чем модуль первого. Это позволяет складывать их непосредственно, т. к. потеря точности при этом возникнуть не может. Для суммирования знакопостоянного ряда (16) применялась схема Горнера, сокращающая объем вычислений и минимизирующая погрешность.

Вычисление разности долгот. Разность долгот вычисляется по формуле:

ДА = AQ • tan а. (17)

Отсюда легко получить А2 = А1 + ДА. Относительная погрешность ДА не превышает нескольких машинных эпсилон, что подтвердилось всесторонними экспериментами во всем допустимом диапазоне входных параметров. Вместе с тем выявился набор исходных данных, который отвечает движению к околополярной точке с многочисленным количеством витков вокруг полюса. В этом случае погрешность в определении второй точки могла составить до 360 нм. Справедливости ради следует отметить, что такая траектория лишена какого-либо практического смысла.

2

Обратная геодезическая задача

Даны: широта ф1 и долгота первой точки и широта ф2 и долгота Х2 второй точки. Нужно определить азимут а и дистанцию, которую нужно пройти, двигаясь постоянно под углом а к направлению на северный полюс, чтобы из первой точки попасть во вторую.

Основные этапы решения обратной геодезической задачи для локсодромии:

1) Вычисление разности изометрических широт Л0.

2) Вычисление азимута а.

3) Вычисление меридионального расстояния ЛМ.

4) Вычисление расстояния между точками Б.

Вычисление разности изометрических широт. Эта величина определяется по методике варианта 4, только вместо ф1 + Лф записываем ф2 :

AQ = ln

tanl П +

tan I П + ф,

+ ® ln

4 2

(1 - e sin( ф2))(1 + e sin ф,)

(18)

(1 + e sin( ф2))(1 - e sin ф1)

Первое слагаемое принимает требуемый

Ф2 -Ф1

а второе слага-

sin

вид для y =

2

cos

Ф2 +Ф1 2

Ф2 - ф1 ф2 + ф1

-2 sin ——— cos ———

емое — при y =

л 2 • • Ф2 + Ф, 1 - e sin ф1 sin ———

Вычисление азимута:

а = arctan

X 2 - X, AQ

(19)

С учетом знаков числителя и знаменателя а однозначно определяется в диапазоне [-п; п].

Погрешность при определении азимута, как показали многочисленные эксперименты, не превышает 3е.

Вычисление меридионального расстояния:

<^2

AM = J

a(l - e2)dф (1 - e2 sin2 ф)3/2'

(20)

Это выражение записывается как М(ф2) - М(ф1). Уменьшаемое и вычитаемое раскладываются в ряд (2) по степеням п, затем ряды вычитаются подобно тому, как это делалось в предыдущем разделе.

Возможен и другой подход. В этом случае применяется метод Карлсона, описанный выше.

Вычисление расстояния:

ЛМ

D =

(21)

cos а

Специальный случай возникает при п

а = ± — и cos а = 0, что отвечает движению вд оль параллели. Здесь можно легко определить D, зная радиус параллели

„ а cos ф,

R = . и радианную меру дуги

дД - e2 sin2 ф, по известным долготам X, и X2. Погрешность при определении расстояния не превышает 5е.

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

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

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

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

• Для вычисления разности изометрических широт предложены два метода (вариант 1 и вариант 4) и выполнено их сравнение с использующимися (вариант 2 и вариант 3). Вариант 3 (Паде-аппроксимация) оказался заметно хуже всех по точности. Вариант 4, хотя и уступил по объему вычислений варианту 2 на 14 %, однако, в отличие от последнего, гарантированно работает на всех платформах. Еще один предложенный метод (вариант 1) также может быть рекомендован для решения поставленной задачи. При этом, если одна из точек расположена на экваторе, объем вычислений и точность методов (варианты 1 и 4) совпадают, а в остальных случаях эти характеристики для варианта 4 несколько лучше.

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

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

СПИСОК ЛИТЕРАТУРЫ

1. Alexander J. Loxodromes: A Rhumb Way to Go // College Mathematics Journal. 2004. No. 77(5). Pp. 349-356.

2. Weintrit A., Kopacz P. A Novel Approach to Loxodrome (Rhumb-Line), Orthodrome (Great Circle) and Geodesic Line in ECDIS and Navigation in General // Methods and Algorithms in Navigation, Marine Navigation and Safety of Sea Transportation. Leiden: A Balkema Book, CRC Press, Taylor & Francis Grup, 2011. Pp. 123-132.

3. Karney C.F.F. Algorithms for geodesics // J. Geodesy. 2013. No. 87(1). Pp. 43-55.

4. Petrovic M. Differential Equation of a Loxodrome on the Spheroid // Nase more. 2007. No. 54(3-4). Pp. 87-89.

5. Pallikaris A., Paradissis D., Tsoulos L. New calculation algorithms for GIS navigational systems and receivers // European Navigational Conf. ENC-GNSS. Naples, Italy, 2009.

6. [Электронный ресурс]/ URL: http://koti. mbnet.fi/jukaukor/loxodrome.pdf

7. Deakin R.E., Hunter M.N. Geometric Geodesy. Part A // Lecture Notes. School of Mathematical & Geospatial Sciences, RMIT University, Melbourne, Australia, 2008. 140 p.

8. Kos S., Pogany T.K. On the mathematics

of navigational calculations for meridian sailing, Solstice // J. of Geography and Mathematics. 2012. No. 23(2).19 p.

9. Weintrit A., Kopacz P. On Computational Algorithms Implemented in Marine Navigational Software Used in Marine Navigation Electronic Devices and Systems // Annual of Navigation. 2012. No. 19. Part 2. Pp. 171-184.

10. Helmert F.R. Die mathematischen und physikalischen Theorem der höheren Geodäsie // Die mathematischen Theorem. Leipzig, 1880. Vol. 1.

11. Clenshaw C.W. A note on the summation of Chebyshev series // Math. Tables Aids Comput.1955. No. 9. Pp. 118-120.

12. Carlson B.C. Numerical computation of real or complex elliptic integrals // Numerical Algorithms. 1995. No. 10. Pp. 13-26.

13. Carlson B.C. A table of elliptic integrals of second kind // Math. Comp. 1987. No. 49. Pp. 595-606.

14. Устинов С.М., Зимницкий В.А. Вычислительная математика. СПб.: БХВ-Петербург, 2009. 336 с.

15. [Электронный ресурс]/ URL: http://boost.

org

REFERENCES

1. Alexander J. Loxodromes: A Rhumb Way to 2. Weintrit A., Kopacz P.A Novel Approach to

Go,College Mathematics Journal, 2004, No. 77(5), Loxodrome (Rhumb-Line), Orthodrome (Great Pp. 349—356. Circle) and Geodesic Line in ECDIS and Navigation

in General,Methods and Algorithms in Navigation, Marine Navigation and Safety of Sea Transportation. Leiden: A Balkema Book, CRC Press, Taylor & Francis Grup, 2011, Pp. 123-132.

3. Karney C.F.F. Algorithms for geodesies, Journal Geodesy, 2013, No. 87(1), Pp. 43-55.

4. Petrovic M. Differential Equation of a Loxodrome on the Spheroid, Nase more, 2007, No. 54(3-4), Pp. 87-89.

5. Pallikaris A., Paradissis D., Tsoulos L.New calculation algorithms for GIS navigational systems and receivers, European Navigational Conference ENC-GNSS. Naples, Italy, 2009.

6. Available: http://koti.mbnet.fi/jukaukor/ loxodrome.pdf

7. Deakin R.E., Hunter M.N. Geometric Geodesy. Part A, Lecture Notes. School of Mathematical & Geospatial Sciences, RMIT University, Melbourne, Australia, 2008, 140 p.

8. Kos S., Pogany T.K. On the mathematics of navigational calculations for meridian sailing, Solstice, Journal of Geography and Mathematics,

2012, No. 23(2), 19 p.

9. Weintrit A., Kopacz P. On Computational Algorithms Implemented in Marine Navigational Software Used in Marine Navigation Electronic Devices and Systems, Annual of Navigation, 2012, No. 19, Part 2, Pp. 171-184.

10. Helmert F.R. Die mathematischen und physikalischen Theorem der höheren Geodäsie, Die mathematischen Theorem, Leipzig, 1880, Vol. 1.

11. Clenshaw C.W. A note on the summation of Chebyshev series, Math. Tables Aids Comput., 1955, No. 9, Pp. 118-120.

12. Carlson B.C. Numerical computation of real or complex elliptic integrals, Numerical Algorithms, 1995, No. 10, Pp. 13-26.

13. Carlson B.C. A table of elliptic integrals of second kind, Math. Comp., 1987, No. 49, Pp. 595-606.

14. Ustinov S.M., Zimnitskiy V.A. Vychislitelnaya matematika, St. Petersburg: BKhV-Petersburg Publ., 2009, 336 p. (rus)

15. Available: http://boost.org

БОТНЕВ Виктор Александрович — старший инженер-программист ЗАО «Транзас». 199178, Россия, Санкт-Петербург, Малый пр. В.О., д. 54, корп. 4, лит. В. E-mail: [email protected]

BOTNEV, Victor A. Transas.

199178, Maly Ave. 54-4, V.O., St. Petersburg, Russia. E-mail: [email protected]

УСТИНОВ Сергей Михайлович — профессор Санкт-Петербургского государственного политехнического университета, доктор технических наук.

195251, Россия, Санкт-Петербург, ул. Политехническая, д. 29. E-mail: [email protected]

USTINOV, Sergey M. St. Petersburg State Polytechnical University. 195251, Politekhnicheskaya Str. 29, St. Petersburg, Russia. E-mail: [email protected]

© Санкт-Петербургский государственный политехнический университет, 2014