Способ решения навигационной задачи

В морской практике навык решения навигационных задач очень нужен и важен.
Конечно, сейчас существует целая прорва различных навигаторов, картплоттеров и прочего оборудования, НО, НИЧЕГО НАДЁЖНЕЕ БУМАЖНОЙ КАРТЫ ещё не придумано! Ведь любая техника может выйти из строя по множеству причин (а в море их тьма-тьмущая) и тогда спасти Вас, ваше судно и экипаж сможет только Ваше умение по старинке пользоваться картой и компасом.

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

1. ПРЯМЫЕ задачи на вычисление истинного курса (ИК- направления куда мы реально идём), когда мы знаем компáсный курс (КК — на практике он снимается непосредственно с компаса) и склонение (d)

2. ОБРАТНЫЕ задачи — когда нам известны ИК (т.е., истинное направление, чтобы прийти туда, куда надо) и склонение (d), нам нужно вычислить КК (ведь истинные и компасные курсы различаются на величину поправки магнитного компаса — ΔΚ и если это не учесть, то можно прийти туда, куда вовсе не собирались).

3. Задачи на вычисление:
ИК — истинного курса;
МК — магнитного курса;
КК — компасного курса;
КУ — курсового угла;
ПУ — путевого угла;
ИП — истинного пеленга;
ОИП — обратного истинного пеленга;
МП — магнитного пеленга;
ОМП — обратного магнитного пеленга;
КП — компасного пеленга;
ОКП — обратного компасного пеленга, когда известны некоторые из этих данных, а также, вычисление:
• Cклонения (d) для различных районов плавания;
• Девиации (δ) для разных КК;
• ΔΚ — поправки магнитного компаса.

4. Задачи на вычисление дальности видимости объектов наблюдения.

И для самых дотошных (или продвинутых)

5. Задачи по учёту дрейфа (под действием ветра) и сноса (под действием водных течений).

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

А сейчас приступим, собственно, к решениям задач из «Дельты»:

1. Решение ПРЯМОЙ задачи, вычисление девиации и ИК, когда известны КК и склонение.

В сущности, это довольно просто и быстро, если вы помните школьный курс математики 🆒 и обладаете задатками пространственного мышления. Тем более, что самое мучительное — вычисление значения девиации для заданного КК в этом случае нужно произвести всего один раз 😁

Итак, для примера решим следующую задачу:
КК=003°, d=-2,7° (т.е.W), рассчитать δ по таблице и вычислить ИК. Значение округлить до целого.

Девиация берётся из таблицы для КК, но, в ней даны значения только для целых десятков градусов, а у нас КК=3°, поэтому для начала нам нужно вычислить δ для 3°.
Это делается методом интерполяции.
Берём два ближайших значения, у нас это δ0°и δ10°.
Теперь вычисляем среднее значение (Δδ) для одного градуса на данном отрезке:
Δδ=(δ0°- δ10°)/10= (-2,5 -(-0,5))/10
НЕ ЗАБЫВАЕМ ПРО ЗНАКИ (кстати, формулы алгебраические поэтому «-»на«-» даёт«+», а «+»на«-» даёт «-»)!! И, надеюсь, не надо обьяснять почему делим на десять?
Итак, получилось Δδ= -0,2°

Затем, Δδ умножаем на количество градусов, которое нужно добавить (или отнять) к(от) целому(го) десятку(а), т.е., в нашем случае получается 3Δδ и 7Δδ.
Теперь узнаём δ3° и заодно проверяем:
δ0° — 3Δδ = δ10° + 7Δδ
-2,5 — 3(-0,2) = -0,5 + 7(-0,2)
-1,9 = -1,9°
Уравнение верное, УРА. Мы нашли девиацию для нашего компасного курса.
δКК= -1,9° т.е. Вестовое (W, западное)

Да, кстати, при построении формулы
δ0° — 3Δδ = δ10° + 7Δδ
нужно учесть в какую сторону изменяется значение девиации и соответственно расставить знаки.
У нас от 0° до 10° девиация уменьшается, поэтому δ0° — 3Δδ, а в обратную сторону увеличивается и, мы прибавляем δ10° + 7Δδ.

Закипели? Но это ещё не всё! Нам же ИК нужно узнать!

А тут уже всё элементарно, для этого у нас есть формулы:

ΔΚ=d+δ,
где ΔΚ — поправка магнитного компаса,
d — склонение (известно из условия задачи),
δ — девиация (которою мы только что вычислили)

ИК = КК+ΔΚ = ΚΚ+(d+δ) =
3°+(-2,7°- 1,9°)= 358,4°
(на шарике-то у нас где 0°, там и 360°, так что вычисляем правильно)

По условиям задачи округляем до целого, получается ИК=358°
Вуаля!

Теперь вроде всё объяснил.

Р.S. Таблица девиации для предлагаемых задач везде одна и та же, поэтому в дальнейшем рассмотрении дублировать её не буду. Пусть здесь повисит:

Источник

Через всю географию: навигационные и геодезические задачи на разных языках

Приветствую вас, глубокоуважаемые!

«… истинное место судна хотя и неизвестно, но оно не случайно, оно есть, но неизвестно в какой точке» Алексишин В. Г. и др. Практическое судовождение, 2006. стр. 71

«С двух краев галактики вышли пешеходы. » (С) Сергей Попов (Астрофизик)

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

расстояние между двумя географическими точками
определение точки по известной, расстоянию до нее и азимутальному углу
определение положения точки по измеренным дальностям до известных точек (TOA, TOF)
определение положения точки по измеренным временам прихода сигнала (TDOA)

Все это на C#, Rust и Matlab, на сфере и эллипсоидах, с картинками, графиками, исходным кодом — под катом.

А это, релевантная КДПВ:

Для тех, кто спешит (я и сам такой), вот репозиторий на GitHub, где лежат все исходники с тестами и примерами.

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

C#
Rust
Matlab

Наиболее полная реализация на C#: в отличие от остальных в ней присутствуют методы т.н. виртуальной длинной базы — это когда объект, положение которого необходимо определить неподвижен, и есть измеренные дальности до него из разных точек, с известным положением.

Чтобы посмотреть, как все работает, с какими параметрами вызывает и что возвращает, и провести разведку боем, есть разные демки и тесты:

Тестовое консольное приложение на C#
Тест всей библиотеки на Matlab
Демонстрационный скрипт по TOA/TDOA с красивыми картинками на Matlab
Скрипт на Matlab для сравнения точности решений геодезически задач на сфере (Haversine equations) и на эллипсоиде (Vincenty Equations)
Для реализации на Rust в коде библиотеки присутствуют тесты. И можно посмотреть как все работает просто запустив команду «Cargo -test»

Я постарался сделать библиотеку как можно более независимой и самодостаточной. Чтобы при желании можно было просто взять нужный кусок (сославшись на источник, конечно), не таская за собой все остальное.

Почти всегда углы — в радианах, расстояния в метрах, время в секундах.

Теперь, начнем, пожалуй, с начала:

Геодезические задачи

Есть две типовые геодезические задачи: прямая и обратная.

Если например, я знаю свои текущие координаты (широту и долготу), а потом прошагал 1000 ~~кило~~метров строго на северо-восток, ну или на север. Какие теперь у меня будут координаты? — Узнать, какие у меня будут координаты — значит решить прямую геодезическую задачу.

То есть: Прямая геодезическая задача — это нахождение координат точки по известной, дистанции и дирекционному углу.

С обратной задачей все совсем понятно — например, я определил свои координаты, а потом прошагал сколько-то по прямой и снова определил свои координаты. Найти, сколько я прошел — значит решить обратную геодезическую задачу.

То есть: Обратная геодезическая задача — это нахождение расстояния между двумя точками с известными географическими координатами.

Решать эти задачи можно несколькими способами, в зависимости от необходимой точности и времени, которое вы готовы на это потратить.

Самый простой способ — представить что земля ~~плоская~~ — это сфера. Давайте попробуем.
Вот формула для решения прямой задачи (источник):

Здесь , — широта и долгота исходной точки, — дирекционный угол, отсчитывающийся по часовой стрелке от направления на север (если смотреть сверху), — угловое расстояние d/R. d — измеренное (пройденное) расстояние, а R — радиус земли. , — широта и долгота искомой точки (ту, в которую мы пришли).

Для решения обратной задачи есть другая (не менее простая формула):

Где , и , -координаты точек, R — земной радиус.

Описанные формулы называются Haversine Equations.

В реализации на C# соответствующие функции называются HaversineDirect и HaversineInverse и живут в Algorithms.cs.
В реализации на Rust это функции haversine_direct и haversine_inverse.
И наконец, на Matlab функции хранятся в отдельных файлах и вот обе функции:
HaversineDirect и HaversineInverse

Для C# я буду приводить названия функций и ссылку на файл, где они находятся. Для Rust — только названия функций (коль скоро вся библиотека лежит в одном файле), а для Matlab — ссылку на соответствующий файл скрипта, потому что в Matlab одна функция — один скрипт.

Очевидно, что здесь есть какой-то подвох: земля не сфера, ~~а плоскость~~ и это как-то должно отражаться на применимости этих формул и/или на точности решения.

И действительно. Но для того, чтобы определиться с этим, нужно с чем-то сравнивать.
Еще в 1975 году Тадеуш Винценти (Thaddeus Vincenty) опубликовал вычислительно эффективное решение прямой и обратной геодезической задач на поверхности сфероида (известного более под ником ~~Эллипсоид Революции, товарищ!~~ Эллисоид Вращения), ставшее почти стандартом.

Описание устройства метода тянет на отдельную статью, поэтому я ограничусь лишь отсылкой на оригинальную работу Винценти и на онлайн-калькулятор с описанием алгоритма.

В библиотеке UCNLNav решение прямой и обратной геодезической задач по формулам Винценти лежит в следующих функциях:

C#: VincentyDirect и VincentyInverse в Algorithms.cs.
Rust: vincenty_inverse и vincenty_direct
Matlab: Nav_vincenty_direct и Nav_vincenty_inverse.

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

Самое важное отличие — данные формулы выполняют решение на сфероиде, и его параметры нужно передавать в функции. Для этого есть простая стуктура, которая его описывает.

Во всех реализациях можно в одну строчку получить один из стандартных эллипсоидов. (Сплошь и рядом применяется WGS84 [https://en.wikipedia.org/wiki/World_Geodetic_System] и его приведем в качестве примера):

На C#: В Algorithms.cs есть статическое поле Algorithms.WGS84Ellipsoid — его можно передавать в методы.
На Rust:
На Matlab:

Наименование остальных параметров вполне очевидное и не должно вызвать неясностей.

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

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

Логика его работы скрипта:

1. Берем точку с произвольными координатами

и произвольное направление (я выбрал примерно на запад):

2. Шагаем от нее на все увеличивающуюся дистанцию. Для чего сразу задаемся числом шагов и размером шага:

3. Для каждого шага решаем прямую геодезическую задачу на сфере и на эллипсоиде, получая искомую точку:

4. Для каждого шага решаем обратные геодезические задачи — вычисляем расстояния между результатами, полученными на сфере и эллипсоиде:

5. Проверяем прямые решения обратными для обоих методов:

В скрипте эта последовательность выполняется сначала для шага = 1000 м, а потом для шага = 1 метр.

Сначала посмотрим, насколько отличаются результаты прямых решений по координатам (широте и долготе), для чего вычислим векторы «дельт», благо на Matlab все пишется в одну строчку:

По оси абцисс будем отображать в логарифмическом масштабе, т.к. у нас расстояния меняются от 1 до 10000 км:

В результате получаем такие графики для широты:

Я плохо понимаю в градусах, всегда руководствуюсь методом для прикидки «на глазок»:
1° чего-нибудь это в среднем 100-110 км. И если ошибка больше миллионной или хотя бы стотысячной части градуса — это плохие новости.

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

В результате получаем такую картину:

Получается, что на дальностях 10000 км методы расходятся на 10 км.

Если теперь все повторить для шага в 1000 раз меньше, т.е. когда весь диапазон по оси Х будет не 10000 км а всего 10 км, то картина выходит следующая:

То есть, на дальности 10 км набегает всего 20 метров, а на 1-2 метра формулы расходятся только на дистанциях порядка 1000 метров.

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

Ну, а для тех, кому миллиметровой точности недостаточно, в 2013 году была опубликована работа с описанием решения геодезических задач с нанометровой (!) точностью. Не уверен, что могу сходу придумать, где такое может понадобится — разве что при геодезических изысканиях при постройке гравитационно-волновых детекторов или чего-то совершенно фантастического ).

Теперь перейдем к самому вкусному:

Решение навигационных задач

На данный момент библиотека умеет определять:

Местоположение объекта по дальностям до точек, с известными координатами в 2D и 3D. Такое мы называет TOA — Time Of Arrival (или что более правильно TOF — Time Of Flight)
Местоположение объекта по разностям времен прихода в 2D и 3D. Такое мы называем TDOA (Time Difference Of Arrival).

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

То, что нужно минимизировать, называется функцией невязки.

Для задач TOA она выглядит так:

Где — значение функции невязки для некоей точки с координатами ; N — число опорных точек, имеющих координаты , — измеренные расстояния от опорных точек до позиционируемого объекта.

А для задач TDOA вот так:

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

А вот так эти функции выглядят в коде:

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

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

Остановимся чуть более подробно на внутреннем устройстве библиотеки

На данном этапе библиотека работает с 2D и 3D задачами и сам решатель не знает и не хочет знать как выглядит минимизируемый функционал. Это достигается следующим способом.

У решателя есть две ипостаси: 2D и 3D решатели, основанные на методе Нелдера-Мида или, как еще его называют, метода Симплекса.

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

В C# и Rust 2D и 3D Решатели — Generic-методы:

Пример вызова самого решателя:

Все идентично, с точностью до синтаксиса языка.

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

И соответственно, вызов решателя выглядит так:

Для демонстрации решения TOA и TDOA задач есть специальный скрипт на Matlab.

Демонстрация в 2D выбрана не случайно — я не уверен что могу придумать, как просто и информативно отобразить трехмерную функцию невязки =)

Итак. В начале скрипта есть параметры, которые можно менять:

Положение искомой точки задается случайным образом в указанной области:

Далее, случайно располагаем опорные точки, вычисляем дистанцию от искомой до них и отображаем все:

В итоге получаем такую картинку:

Добавляем к измерениям дистанций случайные ошибки:

Строим функцию невязки для выбранной области с некоей децимацией — иначе расчеты могут занять ощутимое время. Я выбрал размер области 1000 х 1000 метров и считаю функцию невязки по всей области через 10 метров:

Вот так выглядит функция невязки:

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

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

Теперь отобразим все вместе. На поверхности это сделать сложно — разные величины по вертикальной оси. Поэтому удобно все нарисовать на двумерном срезе:

В результате получается примерно так:

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

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

Вот как это выглядит:

Вспомним, что у нас амплитуда случайной ошибки — 1% от дальности, в среднем дальность

200-400 метров, т.е. амплитуда ошибки составляет порядка 2-4 метров. При поиске решения мы ошиблись всего на 70 сантиметров.

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

Строим и рисуем поверхность ошибок:

Получается что-то такое:

И вид «сверху» с опорными точками, реальным и вычисленным положениями искомой точки:

И более детально, расхождение реального и вычисленного местоположения:

В этом конкретном случае решение по TDOA оказалось даже лучше, чем по TOA — абсолютная ошибка составляет 0.3 метра.

Хорошо в модели — всегда точно знаешь, где фактически расположена искомая точка. На воздухе хуже — может быть несколько точек зрения, под водой ты просто что-то вычислил и все — в 99% случаев, чтобы вычислить отклонение от фактического местоположения, его (это местоположение) тоже сначала надо вычислить.

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

Финальный аккорд

Максимально приблизим ситуацию к реальной жизни:

пусть у нас опорные точки имеют встроенные GNSS-приемники и мы знаем только их географические координаты
вертикальная координата нам неизвестна (3D Задача)
мы измеряем только времена прихода сигнала от опорных точек на искомой или наоборот

Такая ситуация описана в самом последнем тесте во всех трех реализациях. Я как-то обделил Rust, и финальный пример разберу на нем.

Итак, самый последний тест в библиотеке. В качестве координат искомой точки я выбрал место в парке, где часто гуляю с собакой.

В результате имеем:
Реальное местоположение (Lat, Lon, Z): 48.513724 44.553248 25
Вычисленное положение (Lat, Lon, Z): 48.513726 44.553252 45.6
Расстояние между точками по поверхности (м): 0.389
Разность по координате Z (м): 20.6

Совпадение «в плане» — очень хорошее, ошибка составляет всего 40 сантиметров, а по вертикальной координате — 20 метров. Почему так происходит предлагаю подумать читателям =)

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

На этом разрешите откланяться, спасибо за внимание. Буду бесконечно рад любому feedback.
Надеюсь, статья и библиотека будут полезны.
Про любые ошибки (грамматические и логические) сообщайте — я исправлю.

На всякий случай приведу здесь ссылку на онлайн (и не только) интерпретаторы Matlab/Octave, которыми пользуюсь сам:

Источник

Способ решения навигационной задачи

Навигационные задачи и методы их решения