Введение

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

Классическое решение для автономной навигации – применение инерциальных навигационных систем (ИНС), дающих «грубые» координаты движущегося аппарата, и использование для корректировки этих координат данных о параметрах различных геофизических полей [1; 2, с. 8–35; 3].

Подходящим, в смысле автономной навигации, является аномальное (связанное с поверхностным слоем Земли) магнитное поле. Очевидно, что оно доступно в любой точке, его практически невозможно намеренно заблокировать или исказить, особенно на расстоянии, оно не зависит от времени суток, сезона и т. д. Проблемой, однако, являются различные естественные возмущения магнитного поля, к которым можно отнести временные изменения, суточные колебания, магнитные бури. Так, если естественные вариации магнитного поля в спокойные периоды не превосходят 5–20 нТл, то во время магнитных бурь они достигают 100–500 нТл [4]. Существуют также проблемы неточности и неполноты измерений его параметров. Для решения проблемы нестабильности значений навигационного поля могут быть использованы два подхода.

Первый – фильтровать навигационные данные. Так поступила, например, компания SandboxAQ (USA) с навигационной системой AQNav (Quantum – AI Navigation System, 2024 год) [5].

Второй – рассматривать навигационные данные как нечеткие величины с учетом ошибок измерения, картографирования и возможных флуктуаций. Такой подход был предложен автором этой статьи в работах [6; 7].

Цель исследования – тестирование нечеткого алгоритма и программного кода для определения местоположения движущегося аппарата по аномальному магнитному полю Земли в заданной области коррекции.

Материалы и методы исследования

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

Пусть имеем некоторую траекторию в заданной области коррекции, вдоль которой проводятся измерения геофизического поля fk* в моменты времени tk, и пусть имеется цифровая карта поля fn(φn,λn), соответствующая высоте/глубине траектории движения аппарата. Здесь φn – географическая широта, λn – долгота n-й точки.

Пусть Sk* = (φk*,λk*) – выдаваемые ИНС координаты местоположения в моменты проведения измерений, а rk – максимальная погрешность ИНС в данной точке маршрута. Для каждой точки Sk*, начиная с первой, рассмотрим поисковую окрестность Ok c центром в точке Sk* и радиусом rk. Пусть Sk – множество точек карты, попадающих в поисковую окрестность Ok, а fk – множество значений поля в этих точках (далее – эталонные данные). Преобразуем четкие значения поля из множества fk и измеренное значение fk* в нечеткие и выполним их попарные нечеткие сравнения, результатами которых являются меры сходства сравниваемых значений.

Для фазификации (введения нечеткости) четких значений поля здесь использовались треугольные функции принадлежности. Однако интервалы неопределенности, определяющие ширину носителей функций принадлежности нечетких значений карты Δэ и измеренных значений Δи рассматривались отдельно, а не объединялись в один, как в [6; 7]. Считалось, что величина Δэ определяется ошибкой картографирования, а Δи задается ошибкой измерения и оценкой возможных флуктуаций магнитных аномалий во время измерений. Такой подход к фазификации представляется более простым и не менее эффективным, чем рассмотренный ранее в [7] способ через введение лингвистической переменной [9].

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

Δи = 2·(δизм + δф). (1)

Следует отметить, что абсолютная погрешность измерений современных магнитометров не превосходит нескольких нТл. При корректной установке датчика и стабильной калибровке влияние магнитных свойств конструкции аппарата также обычно мало [4]. Основной вклад в неопределенность измеряемого значения поля вносит величина δф, обусловленная, например, магнитной бурей. Сложность оценки δф состоит в том, что магнитные бури имеют различные интенсивности (в том числе в зависимости от географической широты), различные частоты и фазы, то есть для различных областей оценки Δи будут различны. В рассмотренной задаче величина Δи менялась в диапазоне 100–400 нТл. Предполагалось, что интервал неопределенности Δи = 100 нТл с учетом ошибки измерения соответствует слабым флуктуациям, например, суточным колебаниям, а интервал неопределенности Δи = 400 нТл –возможным флуктуациям во время магнитных бурь в средних широтах.

Погрешность картографирования, определяющая неопределенность эталонных данных, тоже может изменяться в широком диапазоне значений. Пусть Δэ также принадлежит интервалу 100–400 нТл.

Треугольные функции принадлежности нечетких эталонных Вэ и измеренных Ви значений поля даются формулой

(2)

где b = Вэ, ω = Δэ / 2 или b = Вu, ω = Δu / 2.

Индекс Жаккара для этих функций принадлежности определяется по формуле

. (3)

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

Тогда меру сходства нечетких измеренного и эталонного значений поля предлагается вычислять по формуле

(4)

где k – нормировочный коэффициент, т. ч.

Для описанных выше треугольных функций принадлежности с различными носителями Δu и Δэ:

k = 1 / Jmax , (5)

где Jmax – отношение ширины меньшего интервала неопределенности к ширине большего.

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

Далее переходим к построению эталонных траекторий. Обозначим через Ak подмножества множеств Sk у которых вычисленная мера сходства выше некоторого заданного уровня. Рассмотрим все возможные отрезки, соединяющие точки соседних множеств Ak и Ak+1. Если отрезок удовлетворяет характеристикам реальной траектории по скорости аппарата (длине отрезка) и углу от направления на север с учетом соответствующих погрешностей, то он рассматривается как участок эталонной траектории.

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

В качестве эталонных данных могут быть рассмотрены цифровые карты для определенных регионов [10]. Общедоступными являются данные всемирной карты магнитных аномалий WDMAM [11; 12], которые загружаются в текстовом формате для выбранной пользователем территории. WDMAM имеет разрешение 3 угловые минуты для высоты 5 км над эллипсоидом WGS84. Для перерасчета карты на заданную высоту/глубину могут быть использованы специальные программы, например Geomag-IMP [13].

Результаты исследования и их обсуждение

Пусть аппарат движется прямолинейно в области коррекции между 55°35’ и 55°56’ с. ш., 37° и 37°20’ в. д. Считаются заданными скорость аппарата, азимут и погрешности их определения. Заданы координаты ИНС точек вдоль траектории аппарата, в которых считается измеренным модуль аномальной магнитной индукции. Значения параметров модели приведены в табл. 1, их следует рассматривать как условные, поскольку не ставилась задача привязки к аппарату конкретного вида. При моделировании использовалась готовая цифровая карта аномалий магнитного поля для указанной области коррекции, при уменьшении шага сетки проводилась интерполяция кубическими сплайн-функциями.

Координаты аппарата на выходе из области коррекции, полученные по результатам расчетов, приведены в табл. 2. Координаты аппарата, моделируемые по заданным измерениям поля: (55°54’, 37°10’).

По результатам табл. 1 были вычислены средняя абсолютная ошибка (6), среднеквадратичная ошибка (7), выборочная дисперсия оценок (8). Результаты расчетов показаны в табл. 3.

(6)

(7)

(8)

где – оценки истинного значения xист

Таблица 1

Исходные данные

Примечание: составлена автором на основе полученных данных в ходе исследования.

Таблица 2

Результаты расчетов координат аппарата

Примечание: составлена автором на основе полученных данных в ходе исследования.

Таблица 3

Оценки ошибок

Примечание: составлена автором на основе полученных данных в ходе исследования.

Таблица 4

Сравнение разброса оценок координат с размером ячейки

Примечание: составлена автором на основе полученных данных в ходе исследования

Сравнение полученной выборочной дисперсии с размером ячейки в терминах евклидовой метрики () дано в табл. 4.

По результатам расчетов, показанных в табл. 2–4, можно сделать следующие выводы.

1. Как и следовало ожидать, наиболее точные оценки местоположения аппарата после прохождения заданной области коррекции могут быть получены на сетке с шагом 0,25’. Это связано не с уточнением данных карты, а с более точным учетом геометрических погрешностей участков эталонных траекторий, что также подтверждается характером изменения оценок координат: при шаге 1 минута они «связаны» с одним узлом сетки, при уменьшении шага сетки оценки координат начинают изменяться между соседними узлами.

2. Среднеквадратичная ошибка оценок координат для сетки с шагом 0,25’ в евклидовой метрике составляет ~555 м.

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

4. Сравнение разброса оценок координат (по статистике s2) с размером ячейки карты свидетельствует об устойчивости нечеткого метода к изменениям ширин интервалов неопределенности.

5. В рассмотренном вычислительном эксперименте полученные оценки координат характеризуются положительным смещением (в большей мере по долготе) относительно моделируемого положения аппарата. При использовании нечеткого метода на практике следует либо проводить предварительный расчет для заданных области коррекции и характеристик движения аппарата с целью вычисления систематического смещения, учет которого значительно повысит точность определения координат, либо комбинировать нечеткий алгоритм с другим методом. Например, сначала применить нечеткий метод, а после уменьшения ошибки ИНС до ошибки нечеткого метода, использовать расширенный фильтр Калмана [14; 15].

Для наглядности некоторые полученные расчетные результаты показаны на рис. 1–3. На всех рисунках нижняя прямая изображает траекторию аппарата по данным ИНС, лежащие выше отрезки – участки возможных эталонных траекторий. Точками внутри поисковых окрестностей отмечены те узлы сетки, в которых мера сходства соответствующих нечеткого измеренного и эталонного значений модуля аномальной магнитной индукции выше порогового значения, равного 0,4.

Для оценки вычислительной сложности рассмотренного нечеткого алгоритма введем обозначения: N – полное число точек на карте области коррекции, M – число точек внутри одной поисковой окрестности (зависит от шага сетки карты и радиуса),C – число поисковых окрестностей.

Рис. 1. Построение эталонных траекторий на карте магнитных аномалий области коррекции с шагом сетки 1’ (Δэ = 200 нТл, Δи = 300 нТл) Примечание: составлен автором по результатам данного исследования

Рис. 2. Построение эталонных траекторий на карте магнитных аномалий области коррекции с шагом сетки 0,5’ (Δэ = 400 нТл, Δи = 400 нТл) Примечание: составлен автором по результатам данного исследования

Рис. 3. Построение эталонных траекторий на карте магнитных аномалий области коррекции с шагом сетки 0,25’ (Δэ = 300 нТл, Δи = 400 нТл) Примечание: составлен автором по результатам данного исследования

Учитывая, что сложность одной операции нечеткого сравнения, вычисления углового отклонения и длины участка эталонной траектории оценивается как O(1), можно дать следующие оценки сложности основных функций алгоритма:

− загрузка данных – O(N);

− отбор точек карты, попадающих в поисковые окрестности, – O(CN);

− отбор «лучших» (по результатам нечетких сравнений) точек внутри всех поисковых окрестностей – O(CM);

− построение участков эталонных траекторий по отобранным точкам соседних окрестностей – O(CM 2);

− ранжирование эталонных траекторий – O(1)

Тогда итоговая оценка по доминирующим членам – O(C(N+M 2)).

Заключение

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

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