Сплайн-интерполяция в минимаксной фильтрации для решения задач навигации по геофизическим полям мирового океана1 Текст научной статьи по специальности «Математика»

УДК 681.514

А. Н. Розенбаум, А. И. Дешнер

СПЛАЙН-ИНТЕРПОЛЯЦИЯ В МИНИМАКСНОЙ ФИЛЬТРАЦИИ ДЛЯ РЕШЕНИЯ ЗАДАЧ НАВИГАЦИИ ПО ГЕОФИЗИЧЕСКИМ ПОЛЯМ МИРОВОГО ОКЕАНА1

Введение

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

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

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

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

При отсутствии знания стохастических характеристик самого оцениваемого процесса и связанных с ним наблюдений использование среднестатистических алгоритмов, в том числе и фильтров калмановского типа, для обработки данных измерений сопряжено с принятием ряда труднопроверяемых гипотез (нормальность возмущений, независимость измерений, принадлежность обрабатываемой выборки генеральной совокупности и т. д.). При этом оказывается невозможным оценить достоверность результатов фильтрации. Гарантия достоверности в рассматриваемой ситуации при определении величин коррекции ИНС может быть обеспечена применением неклассической фильтрации калмановского типа. Теоретической основой для такой фильтрации может быть расчет на «наихудший» случай, т. е. принцип минимакса. В настоящей работе предлагается метод использования геофизических полей для построения морских навигационных систем, в котором получение оценки текущих координат объекта производится с помощью теории гарантированного оценивания, что позволяет корректно обрабатывать неслучайные ошибки модели состояния и движения объекта, а также учитывать неполноту и неопределенность исходной информации.

Постановка задачи

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

х1+1 = Агхг + ег, 1 = 1, N, г1 = е0, (1)

при линейных измерениях

У; = Нгхг + вг, 1 = 1, N, (2)

1 Работа поддержана грантом Отделения энергетики, механики и процессов управления РАН 09-1-ОЭМПУ-01, 2009-2010.

где Л*, - известные матрицы соответствующей размерности; е*, в*, - неопределенные векторы

размерности т и п соответственно. Априорная неопределенность в системе и условиях ее функционирования заключена в наборе ю = со1(е^ 8^), где е^ = со1(е0, ..., е^, = со1(8ь ..., влт).

Про вектор ю неизвестных параметров известно лишь, что он стеснен ограничением

N N

21Iе* - т* II В +2118* - Ч* 11с, £ г2, (3)

1=0 ' *=1 '

где В*, С* - положительно определенные матрицы; г - заданная постоянная; т*, ч* - известные

векторы. Множество

Г N N Г

Ег = 1®: 21 Iе*- т*1В. +21К-- ч*11С.£ г2 г (4)

I *=0 ' *=1 '

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

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

Соотношения минимаксного фильтра

Итак, по измерениям ул = со1(уь ., уЛ) требуется оценить в момент времени N + 1 вектор состояния гл+1 системы. Стратегию оценивания * N+1 = ил (ул) будем искать в классе линейных

N

и N (ул )= 2 (^у* + V), (5)

*=1

где V* - весовые матрицы фильтра; V* - векторы настройки фильтра.

Требуется выбрать такие наборы (V, V;) весовых матриц V и векторов V*, / = 1, N, чтобы была минимальной наибольшая погрешность оценивания

Ф^, V) = 8ир |\х N+1 (ю) - и N+1 (у(®))||2 , (6)

юе Ег

т. е.

(V*, V*) = а^ шт Ф^, V) = а^ тт шах||хN+1 (ю) - иN+1 (у (ю))||2 . (7)

В,Ь В,Ь юеЕг

Таким образом, критерий качества алгоритма оценивания является минимаксным.

Синтез минимаксного фильтра [2] приводит к рекуррентным соотношениям:

* * * * тг*(тт 1 \ 1

v = Vлг+1, где Vк+1 = Vк -Ук{Нк ■ Ьк + qk)-Ьк;

V* = 0, к = 1, N, bj+1 = Ibj + т j, Ь0 = 0, j = 0, N;

А А А

V* = Ук • К, где Ук =(Ак+1 - ^к+1 • Нк+1 )Ук+1;

^N+1 = I, к = N,1, Кк = ладГ(с*-1 +НкРкН)-1;

Рк+1 = АкРкЛгк + 5;1 -ЛкРкНгкКгк , р = £0-1.

Ф(У , V ) = г1, где 1 - максимальное собственное значение матрицы Р+;

А А ( А

г к+1( у *) = Лк г к(ук-1) + Кк • I У*- Чк- Нк г к(ук-1)

Приведенные соотношения по форме аналогичны соотношениям, определяющим классический фильтр Калмана. Матрицы Рк в данных соотношениях не являются в общем случае ковариационными матрицами погрешностей оценивания, т. к. в рассматриваемой модели (1), (2) помехи ег-, 8г- являются неопределенными. Оценки у(0ор‘ являются однозначными, несмещенными и удовлетворяют условию сходимости [2].

Использование геофизических полей

Использование фильтрации Калмана в задачах высокоточной навигации требует определения матрицы Н (в уравнениях (8)), которая может быть получена путем аппроксимации областей цифровых карт аналитическими поверхностями. При этом аппроксимация может касаться только областей неопределенности (областей возможных вариаций модельных ошибок и ошибок наблюдения) измерения у(7), лежащих на траектории движения объекта.

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

В данной работе предлагается использование интервальных сплайн-поверхностей, с помощью которых можно не только восстанавливать исходную поверхность, но и учитывать неопределенность, возникающую при наличии ошибок измерений непосредственно геофизических полей с помощью датчиков, ошибок позиционирования при подготовке эталонных карт, а также ошибок, возникающих при преобразовании данных в определенный формат данных ЭВМ. Для представления данных предлагается использование квадратичных однородных интервальных сплайн-поверхностей на интервале zi е [г“, ], где и, I - соответственно верхняя

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

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

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

+ т

к, У0 = т0.

(9)

(8)

Интервал X определяется как множество X = {х | X1 < х < Xй} , где индексы и и I -

соответственно верхняя и нижняя границы интервала. Далее используется X = ^1,Xй ] для обозначения границ интервала X .

XI + xu

Среднее интервала есть X =----------, ширина интервала - w( X) = Xй — X1, внутренность

интервала /(X) = еопу{х | X1 < х < Xй }, где еопу - выпуклая оболочка множества.

Интервальные арифметические операции определяются как обычные алгебраические операции X о у = еопу{х о у | х е X, у е У} , при этом каждая операция дает новый интервал.

Интервальные функции определяются как функции действительного параметра с интервальными коэффициентами. Например, функция / (V) = [/и (V), / (О] при V = а вычисляется как

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

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

V- < V-+1, 0 < , < к + 2(М — 1),

где к - число сплайн-коэффициентов. Тогда действительная функция /(V) в области [о, Vк+2(м —1) ] называется сплайном порядка М или степени М — 1, если многочлен степени М — 1 на каждом подынтервале [^, V,+1] и его первые М — р — 1 производных, где р - количество контрольных узлов, непрерывны на всем интервале [^, Vк+2(М — 1) ]. Более того, производные высших порядков /(V) непрерывны везде, кроме V-, 0 < 1 < к + 2(М — 1). Вектор

Т = {о, ^, к, Vк+2(М —1)} называется вектором контрольных узлов, и значения ^0, ..., Vк+2(М —1)

называются внутренними узлами 5-сплайна.

Имеется рекуррентное соотношение для базисных функций

и}

В-,1 = Ю IВ',м (V) = / ' 1 , Вш —1(0 + I1 +1,М , Х В,+1,М —1(0, М > 1.

I0 +1]. Ч +м —1 — Ч Ч +м — +1

5-сплайн-поверхность определяется прямоугольным множеством контрольных точек , 0 < 1 < т — 1, 0 < ] < п — 1, являющихся вершинами контрольного полиэдра Z, и двумя

векторами Т и S, связанными с параметрами и и V . Для интервального 5-сплайна каждая контрольная точка заменяется интервалом

II т—1п—1

%М,N (u, v),%М,N (u, v)]= 22 %-,}В1,М (и)В],N (v) .

,=0 }=0

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

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

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

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

Алгоритмы триангуляции используют процедуру сортировки для предварительной обработки массива данных, при этом время выполнения сортировки оценивается как 0(п ^2 п). Индекс, образующийся в результате сортировки, хранится в массиве размером Ь < 6п — 9, еще один массив требуется для хранения пар соседних вершин, рассчитанных в результате выполнения триангуляции.

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

ха уа 1

хь уь 1 < 0,

хс ус 1

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

и / \ ^ Ах + Ву + В

Линейная интерполяция точки (х, у) рассчитывается по формуле-------------—-------, где

А = (у — у0 )(г2 — г0)- (г1— г0 )(у 2 — у0),

В = (х2 — х0 )(г1 — г0 )— (х1 — х0 )(г2 — г0 ),

С = (х1 — х0Xу2 — у0 )—(у1 — у0 Хх2 — х0 ),

В = —(Ах1 + Ву1 + Сг1).

Для построения 5-сплайн-поверхности требуется задание внутренних узлов. При этом их количество, в принципе, может быть произвольным, ограничение состоит лишь в том, чтобы количество внутренних узлов превышало порядок поверхности. При этом целесообразно выбрать количество внутренних узлов сплайна таким образом, чтобы обеспечить определенную точность восстановления исходной поверхности. Для выделения в исходном массиве данных сегментов используем деление двумерного массива на прямоугольные области с помощью тернарного дерева. Каждая вершина тернарного дерева имеет не менее двух и не более четырех потомков, вершина дерева представляет всю исходную прямоугольную область, а ветви дерева получаются делением прямоугольных областей на 4 подобласти. Каждому узлу дерева присваивается значение «истина» или «ложь» либо производится дальнейшее деление до получения какого-либо значения (классификация) либо по достижении заданной глубины дерева.

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

Общая задача заключается в построении поверхности X = X(х,у), минимизирующей ошибку аппроксимации и имеющей форму согласно заданному вектору контрольных узлов. При этом интервальная 5-сплайн-поверхность строится таким образом, чтобы объем, ограниченный 5-сплайн-поверхностью и билинейной поверхностью, был минимальным. Определим

2(и, V) ї 2и (и, V), 21 (и, V)]

й %и (и,V), X' (и,V)

как интервальную 5-сплайн-поверхность, где и = Ах + В, V = Су + В, А, В, С, В выбираются из условия и, V е [0,1] . Требуется минимизировать

11 ки к"

11 ки г

р=112 2 К В, (и) в} (^и^, с=—{{22 X, В (и)В}

0 0,=0 у=0 0 0 ,=0 у=0

по Х,иу , X1 у , где ки, ^ - число внутренних узлов (количество контрольных узлов минус

единица), необходимых для каждого сегмента в параметрических направлениях и и V .

Интеграл от 5-сплайн-поверхности вычисляется с помощью выражения

11 ки к"

— к к 1 ки к"

{{22 Р, у В, (и) Ву ^и^ = 22 р, у 0г+N — V' )(5у + М — 5у )

0 0,=0 у=0 М7У ,=0 у=0

где М = N = 3 для двунаправленных квадратичных сплайнов, откуда получаем

ки Г Л ки г

1 Л Л-

^=шII2,У('..N -^)«

і=0 У=0

У +М

1 Л Л-

)•а=-ТЖII2,у«,+» -)(*/+м - *,) ■

і=0 У=0

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

В качестве примера выполнен расчет 5-сплайн-поверхности двух участков морского дна, изображенных на рис. 1, 2. Оба участка представлены цифровой картой на сетке размерами 200 х 200.

а 100

131.7 131.8 131.9 132 132.1 132.2 132.3 132.4 132.5

Рис. 1

а 4-

131.3 131.32 131.34 131.36 131.38 131.4 131 42

Рис. 2

42.5

42.3

42.2

42.1

п

2

п

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

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

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

1000 0

це Р1 было задано значение Ю00 . На рис. 3, а приведен пример расчета оценок коорди-

нат при использовании некоторой модели движения, а на рис. 3, б - графики изменения величины максимального собственного значения матрицы Рк (8). Масштаб изображений - 1 : 100.

Рис. 3

б

а

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

Заключение

В работе рассмотрен метод использования квадратичный 5-сплайнов для представления геофизических полей в задачах навигации по полям Мирового океана. Использование моделей полей в виде сплайнов позволяет получать достаточно приемлемую точность при одновременно небольших требованиях к объему памяти бортовой ЭВМ и в комбинации с минимаксным фильтром калмановского типа позволяет получать оценки необходимой точности при ограниченной информации о стохастических характеристиках наблюдаемых процессов и возмущений модели движения.

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

1. Генике А. А., Малорацкий Л. Г., Фрумович В. Л. Высокоточные системы ближней навигации // Зарубежная радиоэлектроника. - 1980. - № 10. - С. 87-93.

2. Малышев В. В., Кибзун А. И. Анализ и синтез высокоточного управления летательными аппаратами. -М.: Машиностроение, 1987. - 304 с.

3. Пузырев В. А., Гостонина М. А. Алгоритмы оценивания движения летательных аппаратов // Зарубежная радиоэлектроника. - 1981. - № 4. - С. 3-25.

4. Ржевкин В. А. Автономная навигация по картам местности // Зарубежная радиоэлектроника. - 1981. -№ 10. - С. 3-28.

Статья поступила в редакцию 8.05.2009

SPLINE INTERPOLATION IN MINIMAX FILTRATION FOR SOLUTION NAVIGATION TASKS

OF GEOPHYSICAL FIELDS OF THE WORLD OCEAN

A. N. Rozenbaum, A. I. Deshner

Problems of navigation, tracking moving objects and control of dynamic systems are solved using information about coordinates and motion parameters obtained through means of tracking systems and various sensors including gravity, radio and optical sensors. As the speed of motion is increased, the requirements to dynamic accuracy of control systems increase as well. Evolution of methods of the theory of optimal control as well as computing technology allows the development of algorithms and technical means of implementing high accuracy navigation systems and high accuracy tracking systems. The method of quadratic fi-spline interpolation use in the geophysical fields of marine navigation is described. Using the spline models gives rather acceptable accuracy under concurrently low requirements to the volume of board computer memory and combined with Kalman-type minimax filtering techniques it allows estimation of required accuracy in case when stochastic characteristics of the observed processes and the disturbances of the moving model are uncertain.

Key words: geophysical navigation, minimax filtration, spline-interpolation.