Адаптивная обработка данных авиационной гравиметрии тема автореферата и диссертации по механике, 01.02.01 ВАК РФ

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ ИМЕНИ М.В. ЛОМОНОСОВА Механико-математический факультет

На правах рукописи

Дорошин Данила Рубенович

АДАПТИВНАЯ ОБРАБОТКА ДАННЫХ АВИАЦИОННОЙ ГРАВИМЕТРИИ

01.02.01 - іеоретическаи механика

АВТОРЕФЕРАТ

диссертации па соискание ученой степени кандидата физико-математических наук

Москва - 2012

Работа выполнена на кафедре прикладной механики и управления механико-математического факультета Московского государственного университета имени М.В. Ломоносова.

Научный руководитель: Болотин Юрий Владимирович.

доктор физико-математических наук, профессор

Официальные оппоненты: Семенихип Константин Владимирович, доктор физико-математических наук, профессор

Каршаков Евгений Владимирович, кандидат физико-математических наук, старший научный сотрудник

Ведущая организация: Институт физики Земли им. О.Ю. Шмидта Российской академии наук

Защита состоится 14 декабря 2012 г. в 16 часов 30 минут на заседании совета Д.501.001.22 при Московском государственном университете, имени М.В. Ломоносова по адресу: 119991, Москва, Ленинские горы, Главное здание МГУ, механико-математическнй факультет, ауд. 16-10.

С диссертацией можно ознакомится в читальном зале отдела диссертаций (Ломоносовский просп., 27. Фундаментальная библиотека, сектор А - 8 этаж, к.812)

Автореферат разослан 14 ноября 2012 г.

Ученый секретарь

диссертационного совета

кандидат физико-математических наук,

доцент Прошкин В.А.

РОССИЙСКАЯ

ГОСУДАРСТВЕННАЯ БИБЛИОТЕКА 2013

Общая характеристика работы

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

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

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

В лаборатории управления и навигации механико-математического факультета МГУ под руководством H.A. Парусникова, A.A. Голована и Ю.В. Болотина начиная с 1995 года ведутся работы но разработке математического и программного обеспечения для решения задачи авиационной гравиметрии. Диссертация является пртдолжением этих работ.

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

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

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

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

Теоретическая и практическая ценность.

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

Обсуждение работы и публикации. Результаты диссертации докладывались на международной конференции "Современные проблемы математики, механики и их приложе-ний"(Москвн, 2009 г.), на симпозиуме Международной ассоциации по геодезии (IAG, С.Петербург, 20X0 г.), на конгрессе международной федерации автоматического управления (The 18th World Congress of the International Federation of Automatic Control (IFAC), Италия, 2011 г.), на научном семинаре им. А.Ю. Ишлинского (МГУ, 2009, 2010, 2011, 2012 г.), и др.

Основные результаты диссертации опубликованы в пяти работах. Список публикаций приведен в конце автореферата. Работа над диссертацией выполнялась при поддержке РФФИ (проект 10-01-00703-а).

Структура и объем диссертации. Работа состоит из списка обозначений, трех глав, заключения и списка литературы (88 наименований). В работе приведено 24 рисунка, 5 таблиц. Общий объем диссертации составляет 112 страницы.

2 Содержание работы

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

В гравиметрии силу тяжести д принято представлять в виде суммы нормальной силы тяжести д0 и аномалии силы тяжести 6д:

д = до + 5д.

Нормальное поле силы тяжести до представляет иоле модельного эллипсоида вращения, поверхность которого близка к поверхности Земли. В гравиметрии в качестве референц-эллиисоида часто используется эллипсоид WGS-84. Эллипсоид вращения определяет систему координат A. h, где ц> — географическая широта, А — географическая долгота, h — высота над эллипсоидом. Нормальная сила тяжести ga{f< h) является известной функцией географической широты уэ м высоты h, и вычисляется по формуле Сомильяны1. Аномалия силы тяжести (ACT) определяет отклонение реального поля от нормального. Основное влияние на структуру ACT района съемок оказывает локальное распределение притягивающих масс. ACT является неизвестной пепараметрической функцией широты, долготы, и высоты над референц-эллипсоидом Sg(f,X.li).

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

Далее в первой главе рассматривается стандартная методика решения задачи авиационной

JTurge W. Geodesy. 2nd Edition Walter tie Gruyter., Berlin. 1991.

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

1. Определение координат и скоростей антенн бортовых приемников СНС.

2. Определение угловых ошибок горизонтирования платформы ИНС

3. 01гределение координат и скоростей чувствительного элемента гравиметра.

4. Решение основного гравиметрического уравнения — ЗАГТ.

5. Согласование галсов (устранение трендов гравиметра).

С. Построение карты ACT на высоте движения летательного аппарата.

7. Редукция карты на Землю, учет топографических поправок и т.д.

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

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

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

Далее в диссертации рассматривается ЗАГТ. Сила тяжести на траектории рассматривается кок функция полетного времени (:

<7(0 = лМО. Л(0) + àg{4>{t), W). Ht)).

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

V = f + д0 +ög + Адегу. (1)

Здесь / — вертикальная проекция удельной силы, действующей на чувствительный элемент со стороны подвеса; V — вертикальная компонента скорости чувствительного элемента; до — нормальная сила тяжести; ig — ACT; Agsrv -- поправка Этвеша. Измерениями в (1) являются оценка V' вертикальной скорости чувствительно элемента, вычисленная по показаниям СНС, и величина вертикальной проекции удельной силы /', действующей на чувствительный элемент. Уравнения измерений вертикальной скорости и удельной силы в частоте СНС можно записать в виде:

f'(th) = f(tk) + p(tk)..

V'(tk) = V(tk) + r(ifc), (2)

где tk — дискретные моменты времени измерения СНС. tk — 1.- At - интервал дискретизации СНС, риг- ошибки измерений.

Члены до и Ад&]-у в (1) с высокой степенью точности определяются но показаниям СНС. ЗАГТ можно сформулировать следующим образом: определить Sg(t) в (1) по измерениям (2). Далее формулируется ЗАГТ в дискретном времени. Используя (2) и известные go(tk), 9ETv{tk): вычисляются значения z[tк):

VVCitl

z(tk) = —дР - f'(tk) - g0(tk) - gETv(h),

где V — оператор конечной разности: W'(tk) = V'(tk) — V"(iJ._1). Заменив в (1) V на аппроксимацию в дискретном времени v д^', (1), (2) можно переписать в виде:

+ (3)

В дискретном времени ЗАГТ формулируется следующим образом: по измерениям z(tk), с учетом (3). определить Sg(tk)- Из-за наличия шума эта задача математически некорректна без дополнительных предположений о свойствах Sg(t,k) и шумов системы.

Далее в диссертации рассматриваются стандартные методы, позволяющие свести ЗАГТ к задаче оптимального стохастического оценивания, основанные на введении для ACT модели стационарного случайного процесса. Рассматриваются различные модели аномалии, применяемые в АГ. Основное внимание уделяется моделям, которые представляются во временной

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

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

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

В конце этой главы рассматривается вопрос о возможности применения различных адаптивных фильтров для решения ЗАГТ.

Во второй главе представлена разработанная методика адаптивного оценивания силы тяжести с учетом ее пространственной неоднородности. Задача оценивания решалась в стохастической постановке. Путем введения формирующих уравнений для аномалии силы тяжести и погрешностей измерений АГК задача оценивания сводится к задаче оценивания в пространстве состояний. В работе используется следующая система уравнений, связывающая измерения i(ijt) и искомое g(tk):

fv'^fo) ^ q(tk)

= ig{h) + ¿Vr(tt)

Здесь уравнение (3) записано без учета шума гравиметра р(1к), который на практике много меньше шума ^дН' ® качестве модели аномалии взято формирующее уравпсние в дискретном времени вида т-ого интеграла от белого шума, где V"1 - оператор m-ой конечной разности, q(lk) ~ гауссовый белый шум в дискретном времени:

Е [,(**)] = 0. Е [?(«*)?(*,)I = Q(h)Sk.,..

где f>k,i — символ Кронекера. Дисперсия Q(tk) полагается переменной для моделирования пространственной неоднородности. В качестве модели погрешностей СНС в определении скорости используется белый гауссовый шум:

Е [r(tt)] = 0, Е [r(it)r(i,)] = R(tk)6kJ.

Рассматривается модель шума с переменной дисперсией /?(£*)■

эБолотпн Ю. В., Попеленский М. Ю. Алализ точности |>ешспия задачи авиа1равимотрии при идентификации шьраметров гравиметра в полете /,/ Фу идам, и прикл. матем., 2005. 11. вьш. 7. 167-180.

Система (4) описывает процесс, авторегрессии, измеряемый с шумом. Путем введешш вектора состояния = ..., 6д(1/с_Г11+1), т{1к), система (4) сводится к системе в пространстве состояний:

(уМ^АЖь-О + Вф-к)

\г(1к) = Су(Ь) ^

Здесь А, В, С — матрицы системы, е(<л-) — центрированный двумерный гауссовый шум, составленный из формирующих шумов системы: е(1к) = г)^ .

ЗАГТ можно рассматривать как задачу оценивания вектора состояния системы (5). В случае, когда С^^ь), Я(^) — известные функции времени, задача оптимального оценивания вектора состояния может быть поставлена как задача оптимального сглаживания -

задача минимизации дисперсии ошибки оценки.

В авиационной гравиметрии, как правило, отсутствует точная априорная информация об интенсивности силы тяжести и об уровне ошибок СПС — (¡{Ьк)- неизвестны. Для поста-

новки задачи адаптивного оценивания вводится модель дисперсии шумов системы. Дисперсии Л(^) предполагаются кусочно-постоянными функциями, принимающими значения из конечных наборов:

<Э(Ь) € {<?!,<2*,}. Л(и) е {й,,й2,..,ЛЛ1}.

В каждый момент времени t.k дисперсии шумов могут принимать значения из нары {(Яп}-Общее число возможных пар - N = . Для описания дисперсий как функций времени

вводится марковская цепь, состоящая из N состояний. Пусть Е {1,2,..., ¿V} - состояние марковской цепи в момент времени 5 = [.ч(<о)> • • •, - траектория цепи, где Т -

количество отсчетов. Марковская цепь определяется набором переходных а,} и начальных вероятностей:

ац = Р.ч^Ы = = ¿}: г,] € {1,2,... ,ЛГ},

7Г, = Р.ч{и(Ь)=1}, »'€ {1,2,...,ЛГ}.

Все возможные пары дисперсий {<2^in}i нумеруются индексами г е {1.2..... N}. Каждое состояние марковской цепи ассоциируется с соответствующей парой дисперсий. Таким образом, состояние марковской цепи я(Ц) задает значения дисперсий шумов в момент траектория цепи 5 задает дисперсии С}{1к), Щк-) как функции времени. Представленная комбинация марковской цепи и модели п пространстве состояний является частным случаем модели линейной системы с марковскими скачками. Процесс генерации данных представленной модели выглядит следующим образом:

• реализация в момент состояния марковской цепи s(ijt) = г;

• состоянию марковской цепи с номером г соответствует наперед заданная пара {Qi, Ял},, шумы g(tk), r{tk) генерируются с дисперсиями Qi, R„;

• используя q(ti,), r(£fc) и y(tk-1), формируется вектор состояния y(ijt), далее формируется наблюдение z(ifc).

Задача оценивания аномалии на траектории ставится в диссертации как задача оценивания вектора состояния модели по критерию MAB:

y(tk) = arg max. fy(y(tk)\Z), (6)

где Z - [z(to), z(ti),.... ¿(tr-i)] -- набор измерений на галсе. При этом неизвестными являются как некоторые параметры модели Э = {{Qi, Ii,,},- "¡j> так и I1R наблюдаемая последовательность состояний марковской цепи S. Плотность распределения /y(y(f/C>\Z) представляет собой модель гауссовой смеси, что делает 'задачу оптимизации (ü) нелинейной. В литературе существует ряд субоптимальных методов, позволяющих решать (6). Большинство методов в том или ином виде сводятся к редукции (G) на задачу идентификации - оценивания неизвестных парамет1юв системы, и фильтрации - оценивания вектора состояния системы3. Данная методика использовалась в диссертации. Она позволяет оценить аномалию в несколько этапов, каждый из которых является оптимизацией по определенной группе параметров.

1. Идентификация.

• Обучение: оценивание всех возможных дисперсий шумов системы, переходных и начальных вероятностей марковской цепи. Критерий оптимизации - ММП:

в = шк1пах/2(2|в).

• Распознавание: оценивание пути марковской цепи. Оп]>еделеиие дисперс ий шумов системы как функций времени. Критерий оптимизации - MAB:

S = arg max fs(S\Z. 0).

2. Фильтрация: оценивание аномалии при известных дисперсиях шумов системы. Критерий оптимизации - MAB:

y(tk) = aigiiiaxfv(y(tk)\Z,S,e). »(ft)

3Ghahrnmani Z., Hinton G. Variational Learning for Switching Statß-Space Modul« Neural Computation, Vul 12, No. I. 831-864, April 2000.

Задача идентификации параметров рассматривается в стандартной постановке, принятой для классической СММ.

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

Существенной особенностью данных авиационной гравиметрии является низкое отношение сигнал-шум (ОСШ) — энергии аномалии Sg(t-k) на несколько порядков меньше энергии шума У д'*^ ■ В диссертации разработан алгоритм регуляризации данных авиационной гравиметрии, необходимый для решения задачи идентификации и являющийся ее частью. Алгоритм позволяет' сохранить структуру смеси скользящих средних (СС) для регуляризованных данных. Задачи обучения и распознавания решаются для регуляризованных данных.

Задача обучения для СММ-СС ставится как задача ММП. Метод оптимизации - ЕМ-алгоритм (expectation maximization), сходящийся к локальному максимуму фупкции правдоподобия4. Составной частью ЕМ-алгоритма является алгоритм прямого-обратного хода, разработанный в данной работе на случай процесса вида СММ-СС. Также разработан алгоритм обучения, позволяющий использовать данные со всех галсов съемки.

Задача распознавания для СММ-СС решиегся путем МАВ. Конечный радиус корреляции СММ-СС позволил свести задачу к алгоритму типа Витерби, являющемуся частным случаем алгоритма динамического программирования.

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

Идентификация. Измерения градиента аномалии

Задача идентификации решается, используя измерения ш-ой конечной разности аномалии силы тяжести x(tit) = Vm2(ijt), которые представляется в виде смеси СС:

фк)=я(*к) + ~Ут+1г(1к). (7)

Здесь q(tk) — тривиальное СС, отвечающее полезному сигналу - аномалии, ¿¡Vm+1r(i*) — СС. отвечающее шуму СНС. Отметим, что процесс x(t*) = x(t.k)/(VaAt)m, где Va — модуль вектора скорости летательного аппарата, можно рассматривать как измерения тп-ого пространственного градиента аномалии вдоль траектории. Размерность x(tk) — мт_1./с2. В ходе

4Dempster A., Laird N., Rubin D. Maximum likelihood from incomplete data via the EM algorithm. Royal Statistical Society, Series B, Vol. 39, No 1., 1-38, 1977

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

Регуляризация данных

Алгоритм регуляризации заключается в нормализации ОСШ путем сведения задачи в область низких частот, где энергия аномалии силы тяжести и шума соизмеримы. В процессе регуляризации данные авиационной гравиметрии сглаживаются оконным фильтром. В результате сглаживания коэффициенты СС (7) сворачиваются с коэффициентами окна, что приводит к существенному увеличению порядков СС. Для устранения избыточности и уменьшения порядков СС данные прореживаются. После прореживания структура смеси СС теряется. Алгоритм регуляризации позволяет сохранить структуру СММ-СС, для чего решается задача аппроксимации модели сглаженных данных на прореженной временной сетке:

к' V

£=ü /=0

Здесь Tfc — прореженное время: т* — т^-1 = nAt. х'(тк) - сглаженные и прореженные измерения (7), rj, d\, — коэффициенты СС, К', V — порядки СС, г/, г' — аппроксимационные гауссовые шумы в прореженном времени. Коэффициенты и порядки СС получаются в результате решения задачи ^-аппроксимации.Отметим, что алгоритм аппроксимации позволяет сохранить исходный набор дисперсий {Qi, Я,,}, для анпроксимационных шумов </', г'. При дальнейшем рассмотрении (8) штрихи для краткости опускаются.

□бучение

Модель (8) является частным случаем СММ-СС обтцего вида, которая представлена далее в диссертации. СММ-СС общего вида состоит из произвольного конечного числа СС. Шумы, соответствующие различным СС, полагаются попарно независимыми. Динамика параметров модели описывается марковской цепью с конечным числом состояний. В диссертации без потери общности рассмотрены решения задач обучения и распознавания на примере процесса (8), являющегося смесью двух СС. Для дальнейшего изложения значение частоты дискретизации в (8) не важно, далее будем пользоваться индексом измерения t б Z:

К L

*(t)=$>i(i-i) + £<M*-i)- (9)

1=0 1=0

Задача обучения, то есть оценивания параметров СММ-СС, поставлена как задача ММП. Неизвестными являются набор параметров 9 = {a1J; щ, Qi, Я,}, траектория марковской цепи (последовательность состояний) S — [я(-р),..., s(0), s(l)..... я(Т - 1)], где р = max{Ä\ L}, а

s(—p) — начальное состояние марковской цепи. Оценивание происходит при условии известного набора наблюдений X = [х(0). х(1).... ,х(Т — 1)]:

0 = argmax /a-(A"|0).

Обучение проводится путем итерационной оптимизации целевой функции ЕМ-алгоритма, доставляющего локальный максимум функции правдоподобия. На каждой n-ой итерации максимизируется функционал EM-алгоритма по 9 и проводится замена старых параметров 0„ на новые 8„+1:

вп+, = argmaxt/(0,0„).

Функционал EM-алгоритма имеет вид:

1/(0.0„) = 5>6(/xs(A-,S|e)) Ps{S\X, 0„), (10)

s

где fx.s(X, 5|0) — совместная плотность вероятности набора наблюдений X и последовательности состояний S при заданных параметрах 0; Ps{S|A, 0„} — апостериорная вероятность реализации последовательности состояний S при условии известных параметров, вычисленных на предыдущей итерации. Суммирование в (10) проводится по всем возможным последовательностям состояний.

Введем следующие обозначения для апостериорных вероятностей:

7t(0 = PsWi)=*l*!ön}, (11)

1t(i,j) = Ps{s(t-i) = z,s(t)=j\X,en).. (12)

7 ,.(SlT) = Ps{S<_r|X,0„}. (13)

Ведем обозначения для блоков данных: Xft = [x(ti),..., x(i2)]- St' = [s(ti)..... s(t2)]. Используя (11)-(13), целевая функция (10) представлена в следующем виде:

U( в:е„) = 53logr.(_p,7-pW-p)) +

<(-р) Т-1

+ Е Z l°gas(t_i)is(()7((s(i — l).s(t)) + (14)

t=-p+l <(t-l),a(l)

T-1

+ E E bg/.v(x(t)|A-(<:,;.5(l_ 2p,0)7((5it.2p).

Из (14) видно, что оптимизация параметров марковской цепи и параметров СС разделяют-

ся. Получены аналитические выражения для переопенивания параметров марковской цепи. Данные выражения по форме совпадают с аналогичными выражениями для переоцепивания параметров классической СММ:

„+, _ ЕГ=-р+1 УМ.})

: 7-Р(0-

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

Алгоритм прямого-обратного хода

Вычисление вероятностей (11)-(13) является одной из главных сложностей при решении задачи обучения модели линейной системы с марковскими скачкамп. В первую очередь это связано с бесконечным радиусом корреляции данных модели. СММ-СС обладает конечным радиусом корреляции, это позволяет реализовать алгоритм обучения в оптимальной постановке. Вероятности 7( могут быть представлены в следующем виде:

(15)

Основная идея, используемая в алгоритме прямого-обратного хода для подсчета 7, - представить числитель (15) в виде комбинации частей, зависящих от наблюдений, гцюведеиных до момента времени £. и от наблюдений, проведенных после этого момента. Для этого используются так называемые прямые и обратные вероятности. Прямые вероятности, зависящие от наблюдений до момента /, определяются следующим образом:

= 1хМК. |0П). (16)

Обратные вероятности, зависящие от наблюдений после момента I + р, определяются следующим выражением:

¿(¿г*-1) - /А-,.9(^1,51(:г1|.'(о-в„). (17)

Используя прямые и обратные вероятности, получено выражение для вычисления (15):

7/(3-*) = (18)

_ 511-2р1-в„,)о4_1(5'(%)а8((_1),а(1)Д(5Г2'>-1)

= /*(Х|в„) '

Здесь •^1+р'' '- — условная гауссовая плотность, последователь-

ность состояний бГ^Г1 полностью определяет параметры распределения Для пря-

мых и обратных вероятностей получены итерационные соотношения:

а(1-2р) «((+2р)

Начальные значения прямых и обратных вероятностей находятся следующим образом:

р-1

= /х(Хо'^З-р^вп) Ц а,(к-1).,(к)>г.( Р), к=-р-1

Т-1 к=Т-2р+1

С помощью значений прямых вероятностей, вычисляется значение функции правдоподобия, стоящее в знаменателе (18):

сХ-1 Т—2р

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

Распознавание

Задача распознавания траектории СММ-СС решается методом МАВ:

5 - а^шах РЯ(19|Л'. В). (19)

Здесь 5 — распознанная последовательность состояний, 0 — параметры ОММ-СС. В качестве В в (19) выбираются параметры, полученные в результате работы алгоритма обучения. Задача сведена к алгоритму типа алгоритма Витерби, который, и свою очередь, является

алгоритмом динамического программирования. Оптимизация (19) эквивалентна следующей задаче:

S = arg max log fs,x (S. X|9). (20)

Целевая функция (20) представлена в виде:

Т-1

loe/jwW.SÇ'le) + Н> + lüe«.(t-iW*)]-

к=Р

Алгоритм Витерби — итерационный алгоритм. В данной работе алгоритм реализован в обратном времени. Пусть последовательность состояний i^'lop фиксирована. Введем обозначения:

Т-1

MSlzD = max^[log/.v(3'(fc)|А"^,e) -ь loga,№_„^(4]. (21)

■S' k=t

Траектория Sf-\ максимизирующая выражение (21), является оптимальной траекторией, приводящей систему в фиксированную последовательность состояний ¿'/Гзр ■ Для выражения Jt(SltZ2р) справедливо итерационное соотношение:

MSlZ\р) = max [log fx (x(t.)\XlZp, S[_2p., в) + log a.(1_,),,(t) + JM )].

e(t)

Из данного соотношения следует, что оптимальная т раектория прихода системы в последонаг тельность состояний 5,'zip может быть найдена по соответствующей траектории, найденной для момента 1 + 1. Финальная максимизация (19), (20) описана ниже:

rnaxi./,^;1) + log/.v.sW-'.sVie)].

Фильтрация

Проведение идентификации позволяет оценить дисперсии Q(tk), R(.tk) как функции времени и свести задачу оценивания аномалии силы тяжести к оптимальному сглаживанию. Для сглаживания используется алгоритм BF-сглаживания5, не требующий вспомогательных обращений матриц. Также отметим, что в диссертации разработана методика оценивания градиента аномалии. Задача ставится как задача оценивания полезной компоненты процесса (9) (MAB) и сводится к фильтру Випера.

5Bryson А.Е., Frazier M. Smoothing for Linear anil Nonlinear Dynamic Systems / ' Prat:. Optimum Sys. Synthesis Con/., U. S. Air Foirc Tech. Rcpt. ASD-TDB, Feb., pp. 63-119.

Третья глава посвящена проверке разработанной методики адаптивного оценивания. Методика была опробована на данных авиационных съемок района реки Омолон6. Данные обладак>т сильно выраженной пространственной неоднородностью, обусловленной перемежат ющимися горными хребтами и равнинами в районе съемок. Аномалия силы тяжести моделировалась процессом первого порядка (т = 1). Марковская цепь состояла из двух состояний. Погрешности СНС моделировались стационарным шумом.

Частота съема СНС 10 Гц, скорость самолета (Ан-26) 100 м/с. В процессе регуляризации данные сглаживались фильтром Кайзера с частотой среза ш„ = 0.0237 Гц и частотой подавления Ljal„p = 0.05 Гц, что эквивалентно длине волны 4.22 км и 2 км соответственно. Далее данные прореживались в п = 100 раз, что эквивалентно пространственному разрешению в 1 км. В результате регуляризации аномалия и шум моделировались СС седьмого порядка.

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

Qi = 0.079 мГал2, Q2 = 0.443 мГал2, R = 1.36 см2/с*

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

а\ 1 = 0.986, ai2 = 0.014, a2i = 0.027, a22 = 0.973,'

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

На рис. 1 представлены регуляризованные измерения аномалии, градиента аномалии и соответствующие результаты распознавания для одного галса съемки. На графике выделяются два участки со сравнительно интенсивной аномалией (0—20 км п 120—180 км). Участки интенсивной аномалии обусловлены полетом летательного аппарата над горной местностью. На графике результатов распознавания виден участок с частыми переключениями, которые интерпретировались как ошибки работы алгоритма распознавания.

"Полетные данные были предосташюны компанией ЗАО "ГНПП Аэ1>огеофизнка"

Рис. 1: Регуляризованные измерения аномалии (верхний рисунок), градиента аномалии (средний рисунок) и соответствующая оценка траектории марковской цепи (нижний рисунок).

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

164 г

1571-1-i-1-i-i-i-i-i

65.2 65.4 65.6 65.8 66 66.2 66.4 66.6 66.8

latitude

Рис. 2: Карта распознанных состояний (верхний рисунок), сглаженная карта (нижний рисунок). Серым выделено первое состояние (Q1, R). черным - второе состояние (Q2. R).

На рис. 3 приведены результаты оценивания аномалии д,тя всех галсов съемки. При построении оценок использовались сглаженные результаты распознавания. Также в диссер-

Рис. 3: Результаты адаптивного оценивания аномалии для всех полетных галсов.

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

На примере исследуемых полетных данных проведено сравнение адаптивного и неадаптивного оценивания аномалии. Ошибка оценивания определялась отдельно для участков с различной интенсивностью аномалии (рис. 2). Под ошибкой понималось ее среднеквадрати-ческое значение (СКО). Для участка с менее интенсивной аномалией и для участка с более интенсивной аномалией проводилась обработка параметрическим семейством неадаптивных фильтров с заданными постоянными Л. Сравнение проводилось путем варьирования дисперсии формирующего шума аномалии ф или. что эквивалентно, частотой среза фильтра. За значение дисперсии шума СНС принималось значение, полученное на этапе идентификации. Перед определением точности проводилось согласование галсов.

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

164

1а!Ж1с1е

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

8 111! 1 | 1 ■ ■ Г----1-

7 0, = 0.0796 шва!2 а2 = 0.4433 тва!2

1 Х-12.77 кт >. = 8.31 кт

6 V ЭТО = 1.2827 шва! ЭТО = 2.83 тЗа!

И г О 5й Е \ -

О ь 4 ю :

3

2 I * в - -"" \ * ..--'■** 1 1 г 11111

0, т&!12

Ь___I___I____и-___-I.............. I__ „..., I ........ I .............'_-л

12,1 10,1 9,2 8,5 8,1 7,7 7,4 7,2 6,9

X, кт

Рис. 4: Оценка ошибки определения аномалии в зависимости от для участка с менее интенсивной аномалией (пунктирная линия) и более интенсивной аномалией (линия с точками).

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

1. Алгоритм позволил получить оценки параметров, близкие к оптимальным.

2. Использование неоптимальных параметров приводит к существенному ухудшению точности. Например, использование оптимальных параметров, соответствующих одному участку, для оценивания аномалии на другом участке, приведет к ухудшению точности оценивания в два раза.

3. На участке с менее интенсивной аномалией получена оценка точности 1.28 мГал. Длина волны аномалии, амплитуда которой подавлена вдвое — 12.77 км. Ожидаемая точность карты аномалии, с учетом осреднения по соседним галсам, не хуже 0.6 мГал. На участке с более интенсивной аномалией получена оценка точности 2.83 мГал. Длина волны аномалии, амплитуда которой подавлена вдвое — 8.31 км. Ожидаемая точность карты аномалии, с учетом осреднения по соседним галсам, не хуже 1.4 мГал. Данные результаты подтверждаются вычислением невязок па пересекающихся галсах.

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

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

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

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

Для решения задачи идентификации стандартный для классических скрытых марковских моделей подход, основанный на проведении обучения и распознавания, впервые распространен па случай СММ-СС общего вида. Алгоритм дает максимумы соответствующих функционалов правдоподобия и апостериорной вероятности.

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

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

Публикации по теме диссертации

1. Дорошин Д.Р. Применение скрытых марковских моделей для адаптивной фильтрат ции данных аэрогравиметрии // Международная конференция Современные проблемы математики, механики и их приложений. Москва, 2009, С. 277.

2. Болотин Ю.В., Дорошин Д Р. Адаптивная фильтрация данных авиационной гравиметрии /, Симпозиум Международной ассоциации по геодезии (IAG) "Наземная, морская и аэрогравиметрия: измерения на неподвижных и подвижных основаниях "(TG-SMM2010). Санкт-Петербург, 2010, С. 35.

3. Bolotin Yu.V.. Doroshm D.R. Adaptive filtering in airborne gravimetry with hidden Markov chains // The 18th World Congi-ess of the International Federation of Automatic Control (IFAC). Milan 2011, pp. 9996-10001.

4. Болотин Ю.В., Дорошин Д.Р. Адаптивная фильтрация данных авиагравиметрии с использованием скрытых марковских моделей ,'/ Вестник МГУ, Математика, механика. Москва, 2011, е 3, С. 36-42.

5. Дорошин Д. Р. Методика обработки данных авиационной гравиметрии с учетом пространственной неоднородности гравитационной аномалии /7 Электронный журнал Труды МАИ, 50, 2012.

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

13 -3650

2012252249

Подписано в печать 9 ноября 2012 г. Формат 60x90/16. Объем 1,25 п.л Тираж 100 экз. Заказ № 111239

Оттиражировано на ризографе в ООО «УннверПринт» ИНН/КПП 7728572912У772801001 Адрес: 105066, г. Москва, Лефортовский пер., дом 8. корпус 2. Тел. 728-97-17. +7(499)261-78-22. http://www.onlinecopy.ru