CC BY
112
правление подвижными объектами и навигация
УДК 629.7.05
АЭРОМАГНИТНАЯ ГРАДИЕНТОМЕТРИЯ И ЕЕ ПРИМЕНЕНИЕ В НАВИГАЦИИ
Е.В. Каршаков, М.Ю. Тхоренко, Б.В. Павлов
Описаны существующие методы измерения магнитного поля на борту летательного аппарата. Рассмотрен стохастический алгоритм компенсации девиации показаний аэромагнитометра и аэромагнитного градиентометра. Кратко описан алгоритм комплексирования инерциальных и корреляционно--экстремальных навигационных систем. Обосновано преимущество использования поля магнитного градиента в качестве навигационной информации. Приведены результаты численного моделирования работы алгоритма комплек-сирования.
Ключевые слова: корреляционно-экстремальная навигационная система, аэромагнитная градиен-тометрия, магнитная компенсация, интегрированная навигационная система.
ВВЕДЕНИЕ
В настоящее время рядом научных коллективов решаются задачи создания автономных средств навигации по физическим полям Земли — магнитному, гравитационному, рельефа поверхности, оптическому, тепловому и др. [1—4]. Решение задачи навигации в соответствующих системах основано на сопоставлении информации, получаемой с помощью бортовой системы измерения параметров поля, с информацией о поле, хранящейся в памяти бортового компьютера — картой поля или эталоном. Сопоставление осуществляется обычно посредством вычисления некоторого функционала (корреляционной функции) и определения экстремума этой характеристики. Поэтому системы навигации по геофизическим полям называют корреляционно-экстремальными навигационными системами (КЭНС) [5, 6]. Признанными основоположниками этого направления в нашей стране считаются академик А.А. Красовский [7] и профессор В.П. Тарасенко [8], предложившие корреляционно-экстремальные методы навигации по геофизическим полям Земли и радиолокационному изображению земной поверхности в начале 1960-х гг. В дальнейшем при разработке алгоритмов активное применение получили методы ли-
нейной и нелинейной фильтрации [6, 9—12]. Современное развитие теории связано с методами многоальтернативной фильтрации [13] и последовательными методами Монте-Карло [9, 14].
Непрерывное развитие методов решения навигационных задач позволяет разрабатывать новые подходы, обеспечивающие повышение точности и эффективности уже существующих систем, а также создание новых бортовых комплексов.
Для целей навигации особый интерес представляют системы измерения градиента того или иного физического поля. Подобные измерения обладают повышенной чувствительностью и помехоустойчивостью; в отличие от полей поверхностного излучения — видимого, инфракрасного, гамма-излучения — рассматриваемые поля не зависят от инсоляции, метеоусловий или времени года, что дает определенное преимущество при решении навигационных задач.
На современном этапе среди бортовых гради-ентометрических систем можно и необходимо рассматривать векторные и тензорные магнитогра-диентометры, а также тензорные гравиградиенто-метры. Заметим также, что из всех перечисленных систем только магнитные при установке на борту ЛА дают результаты измерений, сопоставимые с наземными по точности и пространственному разрешению [15]. По этой причине в настоящей
работе именно магнитоградиентометрические системы рассматриваются как наиболее перспективные [16]. Рассмотрение градиентометрических систем другого типа возможно с применением тех же подходов.
1. ПРИБОРЫ И МЕТОДЫ БОРТОВЫХ ИЗМЕРЕНИЙ МАГНИТНОГО ПОЛЯ
Количественными характеристиками магнитного поля служат вектор магнитной индукции В (здесь и в дальнейшем будем записывать трехмерные векторы как матрицы-столбцы, состоящие из трех элементов), которая измеряется в теслах (Тл) и вектор напряженности магнитного поля Н, которая измеряется в амперах на метр (А/м). Данные величины связаны между собой материальными уравнениями [17]: В = цц0Н (в системе СИ). Здесь ц — относительная магнитная проницае-
_7
мость среды, ц0 = 4п -10 Гн/м — магнитная постоянная. Вообще говоря, для анизотропных сред параметр ц — это тензор. Однако, поскольку в настоящей работе речь идет об измерениях магнитного поля на борту, этот параметр можно считать постоянной скалярной величиной, равной относительной магнитной проницаемости воздуха ц = 1,00000037.
Таким образом, векторы напряженности и индукции магнитного поля в рассматриваемых системах отличаются лишь масштабным коэффициентом и единицами измерений. В дальнейшем, в том случае, когда это различие не принципиально, будет употребляться словосочетание «вектор магнитного поля».
Для измерения компонент вектора магнитного поля на борту применяются, как правило, ферро-зондовые магнитометры [18]. Разработанные для бортовых применений феррозондовые магнитометры имеют чувствительность порядка 0,1 нТл и погрешность измерений до 10 нТл.
Для измерения модуля вектора индукции магнитного поля, определяемого соотношением
|В | = л/ ВТ В, применяются, как правило, квантовые датчики разных модификаций [18]. Разработанные для бортовых применений квантовые датчики имеют чувствительность порядка 0,001 нТл и погрешность измерений до 0,1 нТл.
Векторным магнитным градиентометром называют конструкцию, включающую в себя несколько квантовых магнитометров, работающих в дифференциальном режиме при фиксированном известном разносе датчиков, который обычно составляет
1—10 м. Чувствительность измерений определяется жесткостью базовой линии между датчиками и ее длиной и для существующих систем имеет порядок 1 пТл/м в полосе 1 Гц [19]. Погрешность же измерения компонент градиента зависит от условий калибровки. Учитывая конструктивные особенности, такие системы трудно калибровать в лабораторных условиях, поскольку даже после калибровки каждого датчика в отдельности останется неучтенным влияние конструкции градиентометра. Оценка поля конструктивных элементов градиентометра выполняется в специальном калибровочном полете на большой высоте. Как правило, градиент внешнего поля в таком эксперименте имеет значение порядка 10 пТл/м, что сопоставимо с погрешностью квантового магнитометра, отнесенного к длине базовой линии.
Чтобы получить компоненты тензора градиента
магнитного поля VB , необходимо продифференцировать все компоненты вектора поля. Очевидно, что при использовании феррозондовых магнитометров тензорный градиентометр получится на несколько порядков грубее векторного, построенного на квантовых датчиках. Однако в последнее время д елаются попытки построить тензорный градиентометр на основе компонентных квантовых магнитометров — СКВИДов (от англ. SQUID — Superconducting Quantum Interference Device — сверхпроводящий квантовый интерферометр) [20].
Заметим, что компоненты вектора индукции магнитного поля, также как компоненты вектора и тензора градиента магнитного поля, измеряются в осях, связанных с измерительным устройством. Для того чтобы перепроектировать значения компонент данных величин в любую другую систему координат, необходимо учитывать ориентацию измерителя во время полета.
2. АНАЛИЗ ИЗМЕНЧИВОСТИ ПАРАМЕТРОВ МАГНИТНОГО ПОЛЯ ЗЕМЛИ НА ОСНОВАНИИ РЕАЛЬНЫХ ДАННЫХ
Для анализа устойчивости параметров магнитного поля во времени использовались результаты аэрогеофизических исследований одного и того же участка в 1999 и в 2011 гг. [21]. Размер участка примерно 7 х 5 км. На рис. 1 (см. третью стр. обложки) приведены карты модуля индукции магнитного поля, соответствующие этим двум съемкам (высота полета составляла 70—80 м). При этом, по горизонтальным осям (на рисунке не показаны) отложены координаты точек измерения в проекции Гаусса — Крюгера, значения измеренной величи-
Рисунки к статье Е.В. Каршакова, М.Ю. Тхоренко, Б.В. Павлова «АЭРОМАГНИТНАЯ ГРАДНЕНТОМЕТРНЯ И ЕЕ ПРИМЕНЕНИЕ
В НАВИГАЦИИ» (см. с. 72 - 80)
54513 64590 54667 54744 54821 54897 54974 54513 64590 54667 54744 54821 54897 54974
Рис. 1.
а -2011 г.; б- 1999 г.
Модуль вектора индукции магнитного поля ¿9 , нТл.:
1
0,012 0,033 0,046 0,059 0,073 0,087 0,103 0,123 0,01 а 0,034 0,048 0,061 0,075 0,089 0,105 0,125 6
0,2 0,9 1,6 2,2 2,9 3,6 4,3 4,9 5,6 6,3 6,9 7,6 8,3
Рис. 2. Разница значений модуля индукции магнитного по ля 2011 и 1999 гг. Ы -|5| , нТл
I 12011 I 11999 '
0,00128 0,01026 0,01923 0,02821 0,03718 0,04615
Рис. 3. Горизонтальный градиент модуля индукции магнитного поля Уя|.£>|, иТл/м:
а — 2011 г.; б — 1999 г.
Рис. 4. Разница градиентов модуля индукции магнитного поля 2011 и 1999 гг.
УяИ20„-УяК999 'НТЛ/М
ны обозначаются цветом в соответствии с выбранной шкалой; для получения непрерывного изображения по данным измерений, проведенных в отдельных точках, применяется билинейная интерполяция. Диапазон изменения индукции магнитного поля на данном участке составляет 54,5—55 мкТл.
Разность измерений приведена на рис. 2 (см. третью стр. обложки). Учитывая точность современных квантовых датчиков магнитного поля, видно, что поле значительно поменялось на всей территории, что ограничивает возможности использования значения модуля индукции магнитного поля для задач высокоточной корреляционно-экстремальной навигации, так как изменение со временем модуля индукции магнитного поля значительно превышает погрешность измерения поля квантовым датчиком.
В то же время, использование градиента магнитного поля имеет ряд преимуществ. Прежде всего, аномалии градиента обусловлены, главным образом, локальными особенностями среды. Они в меньшей степени подвержены глобальным изменениям магнитного поля. Далее, измерения градиента на относительно короткой базе (~ 1 м) позволяют пренебречь вариационной составляющей магнитного поля, чего нельзя сделать для измерений самого поля. Наконец, аномальное поле градиента имеет меньший радиус корреляции, поскольку влияние оказывают в первую очередь приповерхностные объекты, которые дают более контрастные аномалии.
На рис. 3 (см. третью стр. обложки) приведены карты значений горизонтального градиента модуля магнитной индукции (горизонтальный градиент модуля магнитной индукции — величина, определяемая соотношением
^ = /И?)2+(I?) ■
причем оси Оу1 и Оу2 горизонтальны), вычисленных по данным, отображенным на рис. 1. Можно видеть, что на рис. 2 и 3 есть повторяющиеся участки (контуры аномалии в центральной части рис. 2 совпадают с контурами аномалии в центральной части рис. 3, а и б), которые связаны с геологическими особенностями данного региона. Это означает, что временным изменениям подвержено не только нормальное магнитное поле, но и аномальное, которое предполагается использовать для решения задач навигации. Поэтому, при отсутствии актуальных магнитных данных для КЭНС необходимо использовать только аномалии значи-
тельного размера и амплитуды, изменение со временем которых мало. Как следствие, в КЭНС на основе магнитного поля можно применять современные феррозондовые датчики.
На рис. 4 (см. третью стр. обложки) показана разность градиентов, полученных по данным 2011 и 1999 гг. Если исключить из рассмотрения разности, обусловленные аномалиями техногенного характера и сбоями измерительной аппаратуры, максимальное расхождение составляет порядка 0,01 нТл/м. Это соответствует уровню погрешности градиентометров, построенных на основе квантовых датчиков с оптической накачкой. Таким образом, с одной стороны, квантовые векторные магнитоградиентометрические системы могут применяться для решения задач навигации уже при современном уровне развития инструментальной базы, а с другой стороны, поле градиента в достаточной мере устойчиво, чтобы полагаться на карты градиента магнитного поля десятилетиями.
3. КОМПЕНСАЦИЯ ВЛИЯНИЯ МАГНИТНОГО ПОЛЯ ЛЕТАТЕЛЬНОГО АППАРАТА НА ПОКАЗАНИЯ МАГНИТОГРАДИЕНТОМЕТРА
На бортовую измерительную систему, кроме магнитного поля Земли, действует также магнитное поле самого ЛА, связанное с его намагничением и вихревыми токами в его корпусе. Понятно, что поле ЛА создает помехи и снижает точность бортовых измерений. Поэтому возникает задача компенсации влияния поля ЛА на бортовой измеритель магнитного поля (называемая также задачей компенсации девиации).
Методика компенсации влияния магнитного поля ЛА требует проведения калибровочного полета, необходимого для определения параметров модели этих помех. Калибровочный полет предполагает подъем ЛА на максимальную высоту к ~ 1000 м для минимизации влияния аномального магнитного поля Земли. На этой высоте совершается ряд эволюций с углами порядка 5° по рысканию, крену и тангажу, при полете вдоль прямой на четырех существенно различных курсах, необходимых для изменения ориентации вектора магнитного поля относительно носителя.
Для получения измерений в задаче компенсации девиации показаний магнитоградиентометра необходим векторный магнитометр — феррозон-довый магнитометр, чтобы получить направление вектора ВР ^ — поле по данным векторного датчика) магнитного поля. В случае тензорного магнитоградиентометра источником информации
служит любой из датчиков, входящих в градиентометр.
При решении задачи компенсации совместно с результатами измерений градиента и вектора магнитного поля необходима информация о координатах и скорости носителя на этапе калибровки, которая получается по результатам работы приемника спутниковой навигационной системы. Предполагается, что собственное поле ЛА описывается линейной моделью, зависящей от ряда постоянных параметров. Тогда решение задачи компенсации заключается в определении параметров данной модели.
Предлагается описанную задачу компенсации свести к стандартной постановке стохастической задачи оценивания. Таким образом, требуется поставить замкнутую задачу оценивания, т. е. получить формирующие уравнения для аномального поля, уравнения для параметров девиации и уравнения измерений.
Математическое выражение для девиации ДВт (получено Пуассоном в 1824 г.), создаваемой магнитными массами, имеет вид [22—24]:
хастической моделью, параметры которой выбираются на основании анализа реальных данных.
После нормализации параметров с учетом характерных значений параметров девиации, можно переписать выражение В = В0 + АВт + ДВ. для измеряемого вектора магнитного поля с учетом магнитных помех в матричных обозначениях, учитывая порядок малости входящих в них величин, в виде [26]:
Ь = Ь0 + е3К + е3ЬтЬ0 + Е%МТ<е , 0 0 0 ат
(1)
где Ь — безразмерный вектор измерений, связанный с вектором В соотношением В = В*Ь (где В * = 50000 нТл), Ь0 — безразмерный вектор магнитного поля Земли, все элементы матриц К, Ь и М — величины порядка 1, е = 0,1, р0 = Ь0/|Ь0| — вектор направляющих косинусов, т — безразмерное время [27]. Точность выражения (1) изменяется
от 0(е6) до 0(е8) в зависимости от высоты съемки [26]. Переписывая соотношение (1) в тензорных обозначениях, получаем:
ДВт = К + ЬВ0.
Здесь В0 — вектор магнитного поля Земли; К — вектор постоянной или «жесткой» компоненты девиации, обусловленной полем магнитожестких материалов ЛА; Ь — матрица размерности 3 х 3 индуктивной или «мягкой» компоненты, обусловленной полем магнитомягких материалов ЛА.
Девиация ДВ., создаваемая индукционными токами, возникает при вариациях во времени магнитного поля, при неравномерных движениях носителя или при его движении в поле с большими горизонтальными или вертикальными градиентами [25]:
ДВ. = М^,
г (Н
где М — матрица влияния токов размерности 3 х 3.
Модель измерений в задаче оценивания получается, если совместно использовать модели полной девиации ДВыт = ДВт + ДВ. для магнитогра-диентометра и феррозондового датчика.
Для получения замкнутой модели оценивания модель измерений дополняется разностными уравнениями, описывающими компоненты аномального магнитного поля и параметров девиации. В нашей постановке задачи параметры девиации полагаются постоянными. Поле описывается сто-
ар1
Ь = Ьы + е3К. + е3Ь/ + е^мО?.
(2)
При этом считается, что здесь и далее все индексы изменяются от 1 до 3, по повторяющимся индексам производится суммирование и вводится
обозначение Ьь = ^Ь0 ¡Ь10 .
Теперь, используя выражение (2) и вводя обозначение Г.. = У (а Г0 .. = У ¡Ь0) для тензора градиента, модель девиации показаний тензорного маг-нитоградиентометра можно записать как
Г, =
йг.
е 0. + К. + + еМ "ОТ
к 0
ат
(3)
где величины К ., Ь к. и М к. имеют порядок 1. При использовании высокоточных СКВИД-датчиков на этапе начальной калибровки нужно вводить модель градиента.
Из свойств тензора градиента магнитного поля вытекают следующие свойства тензоров намагничения:
КН = КЦ, К.к = Ь
. к
м к = .
К. = 0, Ь\к = 0, М.к = 0.
к
Следовательно, остается 35 параметров, требующих определения.
В рамках стохастического подхода градиент магнитного поля в точках измерения считается дискретным процессом Гаусса — Маркова второго порядка (см. также работу [28]). Данный процесс удовлетворяет разностному уравнению:
(4)
X + 1 = X + иу + 1 = и + ду м,
2 2 Х0, и0 = 0, ду е N(0, ст ), Е[дд] = ст 5у
где Ху — у-е измерение любой из компонент тензора (или вектора) градиента магнитного поля, иу — производная компоненты градиента вдоль направления полета в момент времени у, дj — вектор шума в момент времени у, у. — скорость ЛА в момент времени у, Дt — интервал времени между измерениями, Щ-) — нормальное распределение, Е(-) —
математическое ожидание. Параметр ст2 подбирается, исходя из анализа статистических характеристик магнитного поля по реальным данным.
Заметим, что в модель (4) входит производная компоненты градиента вдоль направления полета. Поэтому для того, чтобы использовать описанную модель при постановке задачи оптимального оценивания, надо предварительно перепроектировать значения магнитного поля, полученные в путевой системе координат, одна из осей которой совпадает с направлением полета, в систему координат, связанную с измерителем. Пусть ось Оу1 совпадает с направлением полета; обозначим через О матрицу перехода от путевой системы к системе, связанной с измерителем.
Для того чтобы математически корректно поставить задачу оптимального оценивания параметров магнитных помех и при этом избежать громоздких выкладок, введем ряд вспомогательных обозначений. Пусть для произвольных матриц X,
V т"> 3 х 3 / т
лк е к (где к — множество действительнозначных матриц размерности к х т):
Х{у} = (Х11, Х12, Х13, Х22, Х23),
Х{у}к = (Х11к, Х12к, Х13к, Х22к, Х23к)
>кх т
и
Х0 = (Мр VlVlBl, ^^ ^^ VlB2, ^^ ^2В2, ^2^, V2Bз, ^Д,).
Введем матрицы:
Л = | 1 УД' 0 1
О =
( \
0 0
ч 0 у!м2ст2 ,
4 = {¿г 4 / / 4 ^Ь
О/ = {О, О, О, О, Оу, ®35х35} — блочно-диаго-нальные матрицы, где индексом у обозначен номер текущего измерения, у — скорость ЛА в момент измерения, 1п и ©пХ к — единичная и нулевая матрицы соответствующей размерности. Обозначим для -го измерения:
Н = Н1, ©5Х1, Я,, ©5x1, Н3, ®5х1, Н4, ®5х1, Н5, ©5x1,
^5, ^5е1, ^5е2, Ье3, ^5 е1, ^5 е2, ^5 е3,
где индекс у всюду опущен, чтобы избежать громоздкости, через е' обозначена конечно-разно -стная производная:
(еу-у-), у >
у Д t ' J '
- _ (е2 - е1)
е1 =
Дt
у = 1,
где е у — вектор направляющих косинусов е1 = Ь/\Ъ'\
в момент времени у, матрицы Н е Я5 имеют вид
(где Як — множество действительнозначных векторов размерности к):
(H1, H2, H3, H4, Н5) =
О2 О11 - О2 О13 О2 О21 О223 Оц О21 - О13 О23 О11О31 - О13О33 О21О31 - О23 О33
О212 О213 О222 О223 О12 О22 - О13 О23 О12О32 - О13О33 О22 О32 - О23 О33
2 Оц О12 2 О21 О22 О12О21 + Оц О22 О12 О31 + О11О32 О22 О31 + О21 О32
2 Оц О13 2 О21 О23 О13 О21 + Оц О23 О13 О31 + О11О33 О23 О31 + О21 О33
2 О12 О13 2 О22 О23 О13 О22 + О12 О23 О13 О32 + О12 О33 О23 О32 + О22 О33 )
Введем вектор состояния х е Я45, х = (х0, К^,
^^л/}!' L^m, А/^ Мш» Мм-» МШ3)Т и вектоР изме-
5Т
рений ^ е Я ^ = Г^, удовлетворяющие, с учетом выражений (3) и (4), соотношениям:
х. + ! = А. + д., д. е N(0, <), Е[д^] = <3.5., Р = Е [х.хТ ь ^ = Н/ х + г., Г е N(0, Д), Е [г^ ] = Д 5., где г. е Я5 — вектор шумов тензорного измерителя
с д иагональной ковариационной матрицей В. Как известно, оптимальная оценка вектора состояния в такой задаче может быть получена с помощью алгоритма дискретного фильтра Калмана, при этом начальные значения х0 = ©45 х1, а Р0 подбирается по результатам анализа экспериментальных данных (см. также работу [27]). Для анализа сходимости оценок возможно использование так называемых стохастических м ер оцениваемости [29].
При измерениях вектора градиента м одуля м аг-нитной индукции задача компенсации принимает иной вид. В этом случае каждая компонента градиента представляет собой разницу измерений скалярных магнитометров. Для каждого из них модель девиации выглядит следующим образом [27]:
Ь = Ь0 + е3 вТрК+ 0,5е3 eTLSbp + е4Ь/М^Т + °(е6).
метрии, что и для тензорного градиентометра. Поэтому после декомпозиции задача (5) распадается на три задачи оценивания, в каждой из которых требуется определить в сумме 16 параметров тензоров К, L и М. Общее число параметров — 48.
Кроме того, заметим, что чувствительность скалярных квантовых магнитометров вполне достаточна, чтобы измерить градиент на высоте калибровочного полета, поэтому для полной модели задачи оценивания необходима модель градиента поля.
Для того чтобы поставить задачу оценивания параметров магнитной девиации, введем, как и в предыдущем случае, ряд вспомогательных обозначений. Пусть:
4 = А А .
< = < 3р ©48Х481 и для -го измерения:
3x1
02, ©3x1,
°3, ©3x1, Т3е1, Т3е2, Т3е3,
Т3е2е3, Т3е1 е1, 13е2е2, Т3е1 е2, Т3е1 е3,
Н = (01, ©
Т 2 Т 2
Т3 е1 , Т3 е2 ,
Т3 е2 е 1, Т3 е2 е3, Т3 е3 е1, Т3 е2 е ),
где индекс . опущен, чтобы не загромождать обозначения, а через 0{ обозначен /-й столбец матрицы поворота О, введенной выше. Обозначим, кроме того:
х0 = (У^, У^В!, У^, У^В У^, У^ДО),
Здесь LS = L + L — симметричная матрица, ер = Ьр/|Ьр | — вектор направляющих косинусов по данным векторного магнитометра. Очевидно, что измерений скалярных д атчиков для решения задачи компенсации не достаточно. При анализе наблюдаемости полученной модели можно показать, что она содержит 16 наблюдаемых параметров
К1, К2, К3, AL11, AL22, L12, L13, L23, АМ11, М12, M13, М21, AM22, M23, M31, М3
АА/ = А/ - ^
32
АМ/ = Ми - м33.
По аналогии с измерением тензора градиента магнитного поля может быть получено выражение для девиации компонент вектора градиента, для которого введено обозначение g. = У .Ь:
= е3£п. + Ке. + Д,, е{ Ь^ +
" у^р ир
йе{
+ еМук ^Ьр + о(е6),
(5)
где тензоры намагничения (К.., Lуk, М.к) не совпадают, вообще говоря, с аналогичными тензорами, фигурирующими в выражении (3). Значит, эти величины уже не обладают теми же свойствами сим-
х/ (К1/, К2 /, К3/, L11 /, L22 /, ^2/, ^^ 13/, L23/, ^ М22 М12Р М13Р М21Р М23 Р M31/, М32г).
Введем вектор состояния х е Я54 х = (х0, х 1,
~ ~ Т 3 Т
х2, х3) и вектор измерений г е Я , г = (^1, £2, £3) ,
удовлетворяющие, с учетом выражений (4) и (5), соотношениям:
х. + 1 = Ах + д, е N(0, <3.), Е[д$] = <2,5
Р = Е [х^ ^
г.. = н.х + г., г.. е N(0, Д), Е [г. = Д 5
3
где г.. е к — вектор шумов векторного измерителя
с диагональной ковариационной матрицей Д. Оптимальная оценка вектора состояния может быть получена, как и в предыдущем случае, с помощью алгоритма дискретного фильтра Калмана, при этом начальные значения х0 = ©54 х 1, а Р0 подбирается по результатам анализа экспериментальных данных.
Отметим, прежде всего, что магнитотометри-ческая система, как уже упоминалось выше, требу-
ет наличия средств измерения углов ориентации — это необходимо, поскольку м од ель поля задается в географической системе координат, а измерения производятся в связанной с ЛА системе координат. Далее, для увязки измерений на различных курсах необходимо ввести условие совпадения векторов градиента, которое даст дополнительное корректирующее измерение.
После завершения процедуры калибровки девиация может быть компенсирована в соответствии с соотношениями для тензорного градиентометра:
Г0И = % - Н(©10хр K{j■}, L{iЦV K{ij}2, ^У^ МЦЛ1, М{У}2, МШ3)Т для векторного градиентометра:
§01 = § - Н (©1 x6, ^ ~X3)T,
где все величины относятся к одному и тому же измерению.
4. АЛГОРИТМ КОМПЛЕКТОВАНИЯ МАГНИТОМЕТРИЧЕСКОЙ И ИНЕРЦИАЛЬНОЙ НАВИГАЦИОННЫХ СИСТЕМ
Алгоритм комплексирования инерциальной и корреляционно-экстремальной навигационных систем (ИНС-КЭНС) рассматривается в следующей постановке. В качестве ИНС выбрана двух-компонентная платформенная ИНС с горизон-тируемой платформой и относительно свободной ориентацией в азимуте. Для определения высоты применяется баровысотомер. В качестве корректирующей информации используется информация о значениях вектора градиента магнитного поля § = (§1, §2, §3)т. Для реализации алгоритмов ком-плексирования необходимо, чтобы в бортовом вычислителе осуществлялось не только счисление навигационной информации по показаниям инер-циальных датчиков, но и решение разностных уравнений для вектора ошибок и его ковариационной матрицы, которое будет необходимо на этапе коррекции. Ковариационная матрица должна вычисляться с учетом вида уравнений ошибок.
Предварительно измеренные значения градиента магнитного поля задаются в виде набора измерений, выполненных в ряде точек пространства, лежащих в некоторой трехмерной ограниченной области. С помощью какого-либо алгоритма интерполяции (в данной работе применялась кубическая интерполяция) по этим измерениям могут быть построены непрерывные функции, связывающие вектор градиента магнитного поля и его производные по у1 и у2 с координатами точки из-
мерения М(у1, у2, у3). Обозначим эти функции как § = Ф(М) и У1§ = У1Ф(М), У2§ = У2Ф(М), при этом будем считать, что компоненты данных векторно-значных функций записаны в гринвичской системе координат.
Алгоритм комплексирования состоит из двух этапов [21]. Первый этап предназначен для грубой оценки ошибок позиционирования по измерениям параметров физического поля. Данный этап наступает, когда начинается движение по эталон -ному участку согласно показаниям непрерывно работающей ИНС. Для этого этапа необходимо проанализировать статистические параметры эталонной карты и определить радиус корреляции для текущей высоты полета. Далее зона возможных значений вектора ошибок ИНС разбивается на некоторое число равных участков, каждый из которых по размеру меньше радиуса корреляции. Для всех участков инициируется траектория в пространстве вектора ошибок, которая рассчитывается согласно уравнениям ошибок с начальными условиями, отвечающими центру зоны. Для каждого момента времени можно определить значение поля по эталону, отвечающее координатам, равным текущему показанию инерциальной навигационной системы. Для всех траекторий необходимо вычислить интеграл, определяющий степень правдоподобия гипотезы. Сравнивая коэффициент правдоподобия с заданным порогом, можно отбросить гипотезы, не удовлетворяющие условиям правдоподобия. Постепенно число гипотез будет уменьшаться. После того, как останется единственная гипотеза, можно переходить ко второму этапу.
На втором этапе задача коррекции ИНС по измерениям вектора градиента м агнитного поля сводится к стохастической задаче оптимального оценивания. Чтобы показать, как это сделать, введем ряд обозначений. Пусть для у-го измерения:
Л" =
¡2 ¡2 А1 ©2 х 3 ©2х2 У
-Ю()/2А1 ¡2 + 2и3 ТА1 ©2х 3 ©2х2 М?
©3x2 ©3 Х2 13 + СО Д1 ©3х 2 /3Д/
©2x2 ©2 х2 ©2 х 3 ¡2 ©2x3
©3x2 ©3 х2 ©3 х 3 ©3 х 2 ¡3
т = I 0 1 -1 0
Н" = (—Ар ©3Х2, -ОтФу (М) + Д2У, ©3Х5),
Фу(М) = ОдФ(М>,
Д1 = От (^Фу(М) + У2Фу(М>), Д2 = (Д1 ©3 Х1),
где М — модельная точка с модельными координатами у е Я3, м? е Я3 — модельная относительная
скорость, и е Я3 — угловая скорость вращения Земли в проекции на оси модельного трехгранника, ю е Я — абсолютная угловая скорость модельного трехгранника, ю0 — частота Шулера, 0Я — матрица поворота от гринвичского трехгранника к модельному трехграннику, 0Т — матрица поворота от приборного трехгранника к трехграннику, связанному с измерителем (ошибкой определения данной матрицы пренебрегаем), индекс измерения . всюду опущен. Кроме того, для произвольного вектора / е Я введены обозначения:
Рис. 5. Результаты численного моделирования алгоритма комп-лексирования
на каждом шаге алгоритма производится коррекция ИНС.
/ =
и
/ =
0 /3 -/2
-/3 0 /1 /2 -/1 0
0 /3 -/2
-/3 0 /1
5. ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ
Запишем для .-го измерения вектор состояния
12 Т Т Т Т Т Т
х е Я : х = (5у , 5м , в , А/ , V ) и вектор из-
32
мерений г е Я : г = 8 — 07Фу(М), где 5у е Я — динамические ошибки модельных координат, 5м е Я — динамические ошибки модельных скоростей, в е Я — угол малого поворота от модельного трехгранника к приборному, А/ е Я — вектор инструментальных погрешностей ньютонометров,
33
V е Я — дрейф гироплатформы [30], 8 е Я — измеренный вектор градиента магнитного поля. Учтем, что вектор 8 может быть записан как
8 = 0Т (Т3 + в )Ф(М), где М' — точка с координатами у' = у — 5у — в У, 5у3 = 0, тогда выражение для г можно переписать в виде (учитывая только линейные члены в разложении в ряд Тейлора)
г = Нтх + г, где г е Я3 — вектор шума датчиков. Значит, введенные выше векторы удовлетворяют соотношениям:
Для того чтобы продемонстрировать эффективность применения описанного алгоритма интеграции данных платформенной ИНС и измерений градиента магнитного поля, было проведено численное моделирование работы данного алгоритма с использованием реальных данных о магнитном поле, полученных в ходе тестового полета на одном из участков. При моделировании предполагалось, что ЛА движется с постоянной скоростью вдоль прямой, при этом его ориентация медленно меняется. Модель ошибок инерциальной системы была подобрана так, чтобы правильно описывать результаты работы датчиков средней точности. В результате была получена зависимость ошибки определения координат ЛА от пройденного им расстояния ^
Результаты моделирования представлены на рис. 5. Сплошной л инией изображена зависимость проекции градиента магнитного поля на направление полета 1; штриховой линией — значение горизонтальной составляющей ошибки определе-
ния положения 5у
22 (у 1 - у1 ) + ( у2 - у2 ) по дан-
ным, полученным в результате комплексирования. Видно, что ошибка быстро убывает до уровня меньше 20 м. Это подтверждает высокую эффективность предложенного решения по интеграции данных измерений градиента магнитного поля и ИНС.
7 + 1
А/пх] + д,, е N(0, 0Ш), Е[д.дТ ]
'Г V
01П5.,
Zj = Н'пх] + г,, . е N(0, Дт)
г )
Е [г. =
0 т п1п
и Я определяются характеристиками инерциальных и магнитных датчиков. Оптимальная оценка компонент вектора состояния дается алгоритмом дискретного фильтра Калмана; с помощью полученной оценки
ЗАКЛЮЧЕНИЕ
Представлен аналитический обзор существующих решений, предназначенных для измерения градиента магнитного поля на борту летательного аппарата. Показано, что ключевые особенности градиента магнитного поля, такие, как временная стабильность и малый радиус корреляции, обусловливают высокую перспективность примене-
ния корреляционно-экстремальных навигационных систем, основанных на измерении градиента магнитного поля. Значительное место в работе уделено задаче компенсации девиации показаний магнитоградиентометра. Для ее решения развит подход, основанный на стохастической модели аномального поля, предложен калмановский алгоритм для оценивания коэффициентов компенсации. Численное моделирование на реальных данных показало эффективное снижение ошибок инерциальной навигационной системы при комп-лексировании ее с магнитоградиентометричес-кой корреляционно-экстремальной навигационной системой.
ЛИТЕРАТУРА
1. Джанджгава Г.И. Навигация по аномальному гравитационному полю Земли. Выбор структуры и обоснование требований к системе навигации с учетом возможностей существующего картографического и аппаратного обеспечения / Г.И. Джанджгава, Л.И. Августов, А.И. Сорока // Авиакосмическое приборостроение. — 2002. — № 6. — С. 63—68.
2. Пешехонов В.Г. Навигационные системы // Вестник Российской академии наук. — 1997. — Т. 67, № 1. — С. 43—52.
3. Щербинин В.В. Алгоритмы фрагментации цветных фотоснимков для формирования разносезонных эталонных изображений оптических корреляционно-экстремальных систем навигации ЛА / В.В. Щербинин, Е.В. Шевцова // Известия ЮФУ. Технические науки. — 2010. — № 3. — С. 87—92.
4. Xiaoming Z. and Yan Z. Analysis of Key Technologies in Geomagnetic Navigation // Proc. SPIE 7128, Seventh International Symposium on Instrumentation and Control Technology: Measurement Theory and Systems and Aeronautical Equipment, 71282J (October 13, 2008); doi:10.1117/12.807129.
5. Баклицкий В.К. Корреляционно-экстремальные методы навигации и наведения. — Тверь: ТО «Книжный клуб», 2009. — 360 с.
6. Белоглазов И.Н. Основы навигации по геофизическим полям / И.Н. Белоглазов, Г.И. Джанджгава, Г.П. Чигин. — М.: Наука, 1985. — 328 с.
7. Красовский А.А. Теория корреляционно-экстремальных навигационных систем / А.А. Красовский, И.Н. Белоглазов, Г.П. Чигин. — М.: Наука, 1979. — 448 с.
8. Белоглазов И.Н., Тарасенко В.П. Корреляционно-экстремальные системы. — М.: Сов. радио, 1974. — 392 с.
9. Bergman N. Recursive Bayesian estimation. Navigation and Tracking Applications. — Sweden: Linkoping: Linkoping University, 1999. — 219 p.
10. Дмитриев С.П. Шимелевич Л.И. Нелинейные задачи обработки навигационной информации: обзор. — Л.: ЦНИИ «Румб», 1977. — 84 с.
11. Степанов О.А. Применение теории нелинейной фильтрации в задачах обработки навигационной информации. — СПб.: ЦНИИ «Электроприбор», 1998. — 369 с.
12. Степанов О.А., Торопов А.Б. Методы нелинейной фильтрации в задаче навигации по геофизическим полям. Ч. 1. Обзор алгоритмов // Гироскопия и навигация. — 2015. — № 3 (90). — С. 102—125.
13. Дмитриев С.П., Степанов О.А. Многоальтернативная фильтрация в задачах обработки навигационной информации // Радиотехника. — 2004. — № 7. — С. 11—17.
14. Степанов О.А., Торопов А.Б. Использование последовательных методов Монте-Карло в задаче корреляционно-экстремальной навигации // Изв. вузов. Приборостроение. — 2010. — Т. 53, № 10. — С. 49—54.
15. Контарович Р. С., Бабаянц П.С. Аэрогеофизика — эффективный инструмент решения геолого-поисковых задач // Разведка и охрана недр. — 2011. — № 7. — С. 3—10.
16. Волковицкий А.К., Каршаков Е.В, Павлов Б.В., Тхоренко М.Ю. Измерения физических полей на борту летательного аппарата при решении навигационных задач // 7-я Российская мультиконференция по проблемам управления // Материалы XXIX Конференции памяти выдающегося конструктора гироскопических приборов Н.Н. Острякова. — СПб., 2014. — С. 232—241.
17. Парселл Э. Электричество и магнетизм: учебное руководство. — М.: Наука, 1983. — 415 с.
18. Telford W.M., Geldart L.R., Sheriff R.E. Applied Geophysics. — Cambridge: Cambridge University Press, 2004. — 744 p.
19. Noriega G., Aeromagnetic Compensation in Gradiometry — Performance, Model Stability, and Robustness // IEEE Geo-science and Remote Sensing Letters. — 2014. Vol. PP, iss. 99 (early publication). — P. 1—5.
20. Foley C.P., Tilbrook D.L., Leslie K.E., et al. Geophysical Exploration Using Magnetic Gradiometry Based on HTS SQUIDs // IEEE Trans. on Applied Superconductivity. — 2001. — Vol. 11, iss. 1. — P. 1375—1378.
21. Волковицкий А.К., Каршаков Е.В., Мойланен Е.В., Павлов Б.В. Integration Magnetic Gradiometer Correlation-Extremal and Inertial Navigation Systems // Proc. of the 19th Saint Petersburg International Conference on Integrated Navigation Systems. — СПб., 2012. — С. 182—184.
22. Каршаков Е.В., Павлов Б.В., Тхоренко М.Ю. Модели и структура бортовых измерений физических полей // Тр. XII Всеросс. совещания по проблемам управления (ВСПУ 2014) / Москва, ИПУ РАН. — М., 2014. — С. 7032—7043.
23. Лысенко А.П. Теория и методы компенсации магнитных помех // Геофизическое приборостроение — 1960. — Вып. 7.
24. Волковицкий А.К., Каршаков Е.В, Харичкин М.В. Система аэромагнитной съемки аномалий магнитного поля // Датчики и системы. — 2007. — № 8. — С. 17—21.
25. Leliak P. Identification and Evaluation of Magnetic Field Sources of Magnetic Airborne Detector Eqquipped Aircraft // IRE Trans. on Aerospace and Navigational Electronics. — 1961. — P. 95—105.
26. Каршаков Е.В, Харичкин М.В. Стохастическая задача оценивания при компенсации девиации аэромагнитометра // Автоматика и телемеханика. — 2008. — № 7. — С. 68—77.
27. Новожилов И.В. Фракционный анализ. — М.: Изд-во мех.-мат. ф-та МГУ, 1995. — 224 с.
28. Болотин Ю.В., Голован А.А., Парусников Н.А. Уравнения аэрогравиметрии. Алгоритмы и результаты испытаний. — М: Из-во Центра прикладных исследований при механико-математическом ф-те МГУ, 2002.
29. Варавва В.Г., Голован А.А., Парусников Н.А. О стохастической мере оцениваемости / В сб.: Коррекция в навигационных системах и системах ориентации ИСЗ. — М.: Изд-во мех.-мат. ф-та МГУ, 1987.
30. Голован А.А., Парусников Н.А. Математические основы навигационных систем. Часть 1. Математические модели инерциальной навигации. — 2-е изд. — М.: Изд-во МГУ, 2010. — 126 с.
Статья представлена к публикации членом редколлегии В.Н. Афанасьевым.
Каршаков Евгений Владимирович — канд. физ.-мат. наук, ст. науч. сотрудник, Н [email protected],
Павлов Борис Викторович — д-р техн. наук, гл. науч. сотрудник, И [email protected],
Тхоренко Максим Юрьевич — ст. инженер, И [email protected],
Институт проблем управления им. В.А. Трапезникова РАН, г. Москва.
