CC BY
262. Точка Р: vïrl} = U+«)p, = )р, = С) + ^^ ^
Работа выполнена при поддержке РФФИ, проекты № 08-01-00231, 08-01-00251 и 08-01-00353.
СПИСОК ЛИТЕРАТУРЫ
1. Александров В.М., Сметанин Б.Н., Соболь Б.В. Тонкие концентраторы напряжений в упругих телах. М.: Наука, 1993.
2. Назаров С.А. Асимптотическая теория тонких пластин и стержней. Т. 1: Понижение размерности и интегральные оценки. Новосибирск: Научная книга, 2002.
3. Никабадзе М.У. Некоторые вопросы варианта теории тонких тел с применением разложения по системе многочленов Чебышева второго рода // Изв. РАН. Механ. твердого тела. 2007. № 3. 73-106.
4. Гольденвейзер А.Л. Построение приближенной теории оболочек при помощи асимптотического интегрирования уравнений теории упругости // Прикл. матем. и механ. 1963. 27, вып. 4. 593-608.
5. Победря Б.Е. Численные методы в теории упругости и пластичности. М.: Изд-во МГУ, 1995.
6. Кривцов А.М., Морозов Н.Ф. Аномалии механических характеристик наноразмерных объектов // Докл. РАН. 2001. 381, № 3. 345-347.
7. Иванова Е.А., Индейцев Д.А., Морозов Н.Ф. Об определении параметров жесткости нанообъектов // Докл. РАН. 2006. 410, № 6. 754-758.
8. Александров В.М. Асимптотическое решение контактной задачи для тонкого упругого слоя // Прикл. матем. и механ. 1969. 33, вып. 1. 61-73.
9. Назаров С.А. Введение в асимптотические методы теории упругости. Л.: Изд-во ЛГУ, 1993.
10. Аргатов И.И. Об улучшении асимптотического решения, получаемого по методу сращиваемых разложений в контактной задаче теории упругости // Журн. вычисл. матем. и матем. физ. 2000. 40, № 4. 623-632.
11. Georgievskii D.V. Asymptotics with respect to a small geometric parameter for solutions of three-dimensional Lame equations // Russ. J. Math. Phys. 2009. 16, N 1. 74-80.
12. Georgievskii D.V. General asymptotic expansions by low geometric parameter in problems of thin solid mechanics / / Proc. 7th EUROMECH Solid Mechanics Conf. Lisbon. 7-11 September 2009. Lisbon, 2009. 35-36.
Поступила в редакцию 16.11.2009
УДК 531.38
АДАПТИВНАЯ ФИЛЬТРАЦИЯ ДАННЫХ АВИАГРАВИМЕТРИИ С ИСПОЛЬЗОВАНИЕМ СКРЫТЫХ МАРКОВСКИХ МОДЕЛЕЙ
Ю.В. Болотин1, Д. Р. Дорошин2
Рассматривается задача адаптивной фильтрации данных авиационной гравиметрии на траектории полета летательного аппарата с учетом неоднородности поля сил тяжести Земли, вызванной неоднородным строением земной коры, рельефом местности и другими факторами. Неоднородность описывается на языке скрытых марковских моделей с конечным числом состояний, где каждому состоянию отвечает определенный характер профиля аномалии силы тяжести. При решении задачи учитывается характер погрешностей спутниковой навигационной системы.
Ключевые слова: гравиметрия, аэрогравиметрия, адаптивная фильтрация, скрытые марковские модели, скользящее среднее.
The paper studies the adaptive filtering of airborne gravimetry data on a survey line, where the inhomogeneous gravity in time is induced by the unknown structure of the Earth's
1 Болотин Юрий Владимирович — доктор физ.-мат. наук, проф. каф. прикладной механики и управления мех.-мат. ф-та МГУ; e-mail: [email protected].
2 Дорошин Данила Рубенович — асп. каф. прикладной механики и управления мех.-мат. ф-та МГУ; e-mail: [email protected].
core, topography, etc. Gravity is modeled by a hidden Markov chain with a finite number of states, each state corresponds to a certain type of the gravity profile. The problem under study is solved with consideration of the type of noises of the global positioning system.
Key words: gravimetry, airborne gravimetry, adaptive filtering, hidden Markov models, moving average.
Задача аэрогравиметрии состоит в определении аномалии силы тяжести по измерениям гравиметра на траектории полета летательного аппарата. В аэрогравиметрии часто используется оптимальное стохастическое оценивание, основанное на стохастической априорной информации об аномалии, которая обычно предполагается стационарным случайным процессом [1]. В настоящей работе предлагается метод адаптивного стохастического оценивания аномалии, для описания которой используется нестационарное формирующее уравнение первого порядка, причем параметры последнего переменны в пространстве, неизвестны и определяются в процессе адаптации. Модель служит для приближенного учета неоднородного распределения притягивающих масс в земной коре. Нестационарная модель аномалии дополняется нестационарной во времени моделью сбоев в измерениях спутниковой навигационной системы. Данные аэрогравиметрии представляются в виде комбинации стационарных участков, где параметры распределений постоянны во времени, и участков переходов между стационарными состояниями, где параметры распределений меняются. Этот процесс описывается при помощи скрытой марковской модели (СММ) смеси скользящих средних (СС). Для построения траектории СММ и последующего определения аномалии разработан алгоритм, основанный на сочетании метода максимума правдоподобия (EM-алгоритм), алгоритма Витерби и сглаживателя Винера.
Для применения предложенного подхода к реальным данным авиагравиметрии в силу слабой обусловленности задачи необходимо ее регуляризовать. Алгоритм регуляризации основан на сглаживании данных и тем самым сведении задачи в область низких частот, где энергии аномалии и шумов измерений сопоставимы.
1. Сведение задачи авиагравиметрии к оцениванию смеси СС. В состав типичной аэрогравиметрической системы входит спутниковая навигационная система (СНС) и набор высокоточных гравиметров (акселерометров), измеряющих так называемую удельную силу — проекцию силы, действующей на чувствительный элемент (ЧЭ) единичной массы со стороны подвеса, на оси чувствительности прибора. Методами инерциальной навигации определяется проекция удельной силы на географическую вертикаль [1]. Запишем основное уравнение скалярной аэрогравиметрии [2] в проекции на географическую вертикаль:
V = /з + go + 5g + A^etv • (1)
Здесь /з — вертикальная проекция удельной силы, действующей на ЧЭ; V — вертикальная компонента скорости ЧЭ; go — нормальное поле силы тяжести; 5g — гравитационная аномалия (в свободном воздухе); AgETV — поправка Этвеша [2]. Члены go и AgETV с высокой степенью точности определяются по показаниям СНС, /з и V измеряются, а аномалия 5g неизвестна и является искомой величиной.
Измерениями в уравнениях (1) служат оценка V' вертикальной скорости ЧЭ, вычисленная по показаниям СНС, и средняя величина удельной силы за такт гравиметра. Частота измерений СНС обычно существенно ниже частоты измерений гравиметра [2], так что измерения гравиметра можно путем интерполяции с высокой точностью привести к частоте СНС, считая, что измеряется интеграл Vf от удельной силы за такт времени СНС. Уравнения измерений в частоте СНС можно записать в виде
t
V;(t)= j /з(т)dT + 5Vf(t), V'(t) = V(t) + 5V(t)• (2)
t-At
Здесь t = kAt, k = 1, — дискретные моменты времени измерения СНС; At — интервал между
измерениями; 5Vf (t) и 5V(t) — ошибки измерений удельной силы и вертикальной скорости. Отметим, что величину Vf' можно интерпретировать как приращение вертикальной скорости ЧЭ, вызванное удельной силой.
Проинтегрировав (1) на интервале At, подставив в него уравнения измерений (2) и перенеся все известные величины в левую часть, получим уравнение в вариациях
t
y(t) = V5V(t)+ j 5g(T) dT• (3)
t- At
Здесь у(Ь) обозначает известную левую часть:
г
у(Ь) = УУ'(1) — У^(Ь) - I [до(т)+АдШу(т)] йт,
г-Аг
где V — взятие конечной разности, соответствующей рассматриваемому интервалу дискретизации: ^У'(Ь) = У'(Ь) — У'(Ь — АЬ). При выводе (3) мы пренебрегли влиянием ошибок гравиметра ЬУ$, так как они, как правило, значительно меньше погрешностей СНС [2].
Задачу аэрогравиметрии можно рассматривать как задачу определения по измерениям (3) величины Ьд. Для дальнейшего следует ввести статистические предположения о составляющих (3). Предполагается, что аномалия Ьд(Ь) является процессом типа т-го интеграла от белого гауссовского шума:
,Я(Ь), Ь = з; VmЬg(t) = д(Ь), М [д(Ь)д(з)] = '
0, Ь = з.
Правомерность данного подхода показана в [3]. Выбор параметра т позволяет строить адекватные модели для аномалий в геологически разных регионах.
В качестве модели шума СНС рассматривается нестационарный центрированный гауссовский белый шум [3]
¡Е(Ь), Ь = з; М [ЬУ(Ь)ЬУ(з)] = [
[0, Ь = з.
Аппроксимируя интеграл в (3) выражением АЬЬд(Ь), далее применяя к (3) оператор Vm и обозначая х(Ь) = Vmу(Ь), т(Ь) = ЬУ(Ь), получаем
х(Ь) = АЬд(Ь) + ^т+1т(Ь), Ь = кАЬ, к = 1,2,.... (4)
Левая часть (3) является функцией измеряемых величин. В правой части стоит сумма (смесь) гауссовских скользящих средних. Таким образом задача аэрогравиметрии свелась к задаче оценивания смеси СС.
С методической точки зрения удобно вместо (4) рассматривать общий случай представления измеряемого сигнала в виде смеси двух скользящих средних: х(Ь) = хг(Ь) + хд(Ь), где
ь к
хг(Ь) = ^ йтт(Ь — тАЬ), хд(Ь) = ^ стд(Ь — тАЬ). (5)
т=0 т=0
Коэффициенты ст, йт называются весовыми коэффициентами СС.
2. Регуляризация данных. Особенность данных аэрогравиметрии состоит в том, что шум СНС на несколько порядков больше полезного сигнала — аномалии. Поэтому задача выделения аномалии плохо обусловлена в том смысле, что второй член в уравнениях (4), (5) на несколько порядков больше первого члена. С целью повышения обусловленности задача сводится в область низких частот, где энергии двух сигналов сопоставимы. Эта операция проводится при помощи сглаживания и прореживания данных. Для сглаживания используется окно Кайзера [4] с частотной характеристикой Н(и). Характеристики окна Кайзера выбираются из априорных соображений о структуре сигнала. Частота подавления фильтра и3 выбирается так, чтобы в области |и| < и3 энергия компоненты аномалии и энергия компоненты шума СНС (4), (5) после сглаживания были соизмеримы. В результате коэффициенты скользящих средних в (5) сворачиваются с коэффициентами фильтра.
Далее наблюдения прореживаются в п раз до интервала дискретизации АЬ' = пАЬ так, чтобы частота Найквиста = п/АЬ' прореженного сигнала х'(Ь), Ь = кАЬ', к = 1, 2,..., совпадала с частотой подавления фильтра и3 [4]. В предположении отсутствия маскировки частот [4] в интервале |и| < и^ спектральная плотность прореженного сигнала х'(Ь) может быть записана в виде Б'(и) = Б'д(и) + Б'г(и), где
Б>(Ш) = |ям|2|я,м|2, = ^ |ян|2|ягн|2, (6)
а Hq(и), Hr(u) — частотные характеристики СС, связанные с коэффициентами cm, dm из (5) формулами разложения в комплексный ряд Фурье.
Фильтрация и прореживание данных в общем случае приводят к тому, что данные не представляются в виде смеси СС. Чтобы вернуть структуру смеси СС, в прореженной временной сетке строится новая аппроксимирующая модель x'(t) = x'q(t) + x'r(t) + e'(t), аналогичная (5), где
K' L'
x'q(t) = 52 Cmq'(t — mAt'), x'r(t) = ^ d'mr'(t — mAt'), (7)
m=0 m=0
q'(t), r'(t) — формирующие шумы в прореженном времени, e'(t) — некоррелированная с q'(t) и r'(t) ошибка аппроксимации.
Аппроксимация строится спектральными методами. Спектральную плотность сигнала x'(t), представленного в виде смеси скользящих средних и ошибки аппроксимации (7), в интервале |и| < w'N можно записать в виде суммы спектральных плотностей компонент: S'x(и) = S"(u) + S"(u) + S"(и), где S'f!(u) — спектральная плотность погрешности e'(t),
= ^ |Я»|2, S'M =
а H'q(и), H'r(и) — частотные характеристики СС, связанные с коэффициентами c'm, d'm формулами разложения в комплексный ряд Фурье.
Порядки скользящих средних K', L' и коэффициенты c'm, d'm выбираются так, чтобы спектральные плотности S'q(и), S'r(и) компонент смеси (7) наилучшим образом аппроксимировали компоненты (6) спектральной плотности прореженного сигнала x'(t):
j \Hq(u)H(и) — H'q(и)|2dw min, J \Hr(u)H(и) — H'r(и)\2(1и dimL min .
Дальнейшие преобразования ведутся с моделью (7), при этом влиянием ошибки аппроксимации пренебрегаем, а штрихи для краткости опускаем.
3. Адаптивное оценивание. Задача оценивания траекторий смеси СС (7) решается в предположении, что интенсивности Q(t), R(t) процессов q(t), r(t) — случайные функции времени. В этих условиях задача относится к классу задач адаптивного оценивания [5]. Для ее решения в настоящей работе была разработана методика адаптивного оценивания, основанная на СММ. Существенным фактором при построении данной методики является конечный радиус корреляции процесса (7).
Введем дополнительные предположения о переменных во времени дисперсиях Q(t), R(t). Первое допущение заключается в том, что дисперсия формирующего шума аномалии Q(t) и шума СНС R(t) может принимать значения соответственно из конечных наборов {Qs, s G [1,... , Nq]}, {Rs, s G [1,... , Nr]}. В конкретный момент времени система может находиться в N = Nq ■ Nr состояниях (s(t) — состояние в момент времени t). Переходы между состояниями определяются вероятностями начальных состояний п и переходными вероятностями aij [5], где i и j — целочисленные индексы состояния, принимающие значения в интервале [1,..., N]. Комбинация уравнения измерений (7) и соотношений марковской цепи для дисперсий составляет скрытую марковскую модель для смеси скользящих средних. Скрытыми в данной постановке являются состояния марковской цепи, которые недоступны наблюдению.
Задача оценивания аномалии по вектору измерений X = (x(0), x(At),...)T может быть поставлена как задача максимизации соответствующей функции правдоподобия. Неизвестными являются набор © = {aij,ni,Qs,Rs} параметров задачи и траектория (последовательность состояний) S = {s(t)} марковской цепи. В данной задаче максимизация функции правдоподобия сводится к алгоритмам с большой вычислительной сложностью. В статье предлагается субоптимальный метод оценивания аномалии, построенный в три этапа, каждый из которых представляет собой оптимизацию по определенной группе параметров:
1) задача обучения. По набору наблюдений X определить неизвестные параметры модели © — матрицу переходных вероятностей марковской цепи, вектор начальных вероятностей, параметры распределений
СС;
2) задача распознавания. По набору наблюдений X при известных параметрах модели оценить статистические характеристики шумов Q(t), R(t) как функции времени, т.е. оценить путь марковской цепи S;
3) задача фильтрации. При известном наборе наблюдений X и известных статистических характеристиках шумов Б, определенных как функции времени, оценить значения компоненты смеси, отвечающей полезному сигналу.
Решение задачи обучения проводится путем итерационной оптимизации функции правдоподобия при помощи ЕМ-алгоритма [2], доставляющего локальный максимум функции правдоподобия. Функционал ЕМ-алгоритма имеет вид
да, @т) = ^ Ъg(fx,S(X, Б1@))Рз{Я|Х, вт}. (8)
в
Здесь вт — оценка параметров на т-й итерации; !(щ(/х,в(X, Б|в)) — совместная плотность вероятности набора наблюдений X и последовательности состояний Б при заданных параметрах в; Рв{БX, вт} — апостериорная вероятность реализации последовательности состояний Б при условии известных параметров, вычисленных на предыдущей итерации. Суммирование в (8) проводится по всем возможным последовательностям состояний. Задача оптимизации заключается в максимизации данного функционала по в и замене старых параметров на новые: вт+\ = aтgmax ©Q(в, вт).
В данной работе ЕМ-алгоритм используется для оценки матрицы переходов и вектора начальных вероятностей марковской цепи. Пусть ) = Рв(з(Ь — 1) = г,з(Ь) = ] X, вт) — вероятность того, что
система находилась в момент Ь — 1 в состоянии г, а в момент Ь — в состоянии ] при условии наблюдаемой последовательности X и известных параметрах вт. Обозначим через ^г(г) = Рв(з(Ь) = IX, вт) аналогичную вероятность пребывания системы в состоянии г в момент Ь. Данные вероятности вычисляются при помощи алгоритма прямого-обратного хода [5], модифицированного на случай СММ, построенной на смеси СС. Приведем выражение для переоценки переходных вероятностей:
т+1 =
ЕЫО '
В качестве переоценки начальных вероятностей пг берутся условные вероятности 7г(г), вычисленные для начального момента времени.
Использование ЕМ-алгоритма для оценки дисперсий формирующих шумов приводит к решению нелинейных уравнений. Для оценки неизвестных дисперсий предлагается субоптимальный линейный алгоритм оценивания, который основан на оценке корреляционной функции процесса для каждого из стационарных состояний. Вид корреляционной функции процесса х(Ь) зависит от набора состояний, при которых реализуются формирующие шумы модели (7). Обозначим множество таких последовательностей состояний через Б1, и пусть Чг(Бг) — условная вероятность реализации последовательности Б1 для момента времени Ь, аналогичная введенным выше условным вероятностям для одного и двух состояний. Оценка
корреляционной функции строится путем взвешивания данных с условными вероятностями 7г(Б1):
= (9)
Ъг Ъ-г (Б\)
Значения корреляционной функции линейно зависят от неизвестных дисперсий. Оценив корреляционную функцию (9) и разрешив систему переопределенных уравнений, можно оценить неизвестные дисперсии системы.
Задача распознавания решается методом максимизации апостериорной вероятности. В данном случае она сводится к задаче динамического программирования (алгоритм Витерби) [5, 6]:
Б = ащтахвР(БX, в). (10)
Здесь Б — распознанная последовательность состояний. В качестве в в (10) выбираются параметры, полученные в результате работы алгоритма обучения.
Задача фильтрации решается как задача минимизации дисперсии ошибки оценки полезной компоненты смеси (7). Оптимальный алгоритм решения данной задачи при известной последовательности состояний марковской цепи — фильтр Винера [5]. Оценка хд(Ь) компоненты смеси хд(Ь) в момент времени Ь (см. (7)) вычисляется путем взвешивания наблюдений X с вектором весовых коэффициентов wг согласно формуле [5]
хд (Ь) = X, т, = Я-£ТХЧ (г)х.
Здесь Кхх = М [XXт] — ковариационная матрица процесса X при известной последовательности состояний, тХд(г)х = М[хд(t)X] — корреляция X с компонентой хд(Ь) в момент Ь.
Рис. 1. Аномалия силы тяжести 6д (кривая 1), регуляризованные измерения уген (кривая 2)
Приведем результаты работы алгоритма с данными, полученными путем использования гравиметра GT1A на прямолинейном участке полета ЛА со скоростью 300 км/ч. На этом участке аномалия имеет три участка с малой и два участка с большой интенсивностью. На рис. 1 приведены графики регуляризо-ванных измерений гравиметра уге8 и истинной аномалии 5д. Регуляризованные измерения уге8 получены путем фильтрации исходных измерений гравиметра у (3) фильтром Н(и) (6) (на рис. 1, а и Ь представлены увеличенные фрагменты графиков для участков с различной интенсивностью аномалии).
Для регуляризации данных проводилось осреднение с разрешением в 2,3 км. Аномалия силы тяжести описывалась интегралом от белого шума, интенсивность которого задавалась марковской моделью с двумя состояниями, шум СНС предполагался стационарным. Для описания регуляризованных аномалии и шума СНС использовались СС седьмого порядка. Алгоритм обучения позволил оценить обе дисперсии формирующего шума аномалии и дисперсию шума СНС: Q\ = 0,37 м2/с4, Q2 = 45,09 м2/с4, К = 7,31 м2/с2.
Сравним полученную оценку аномалии с Рис. 2. Аномалия 5д (кривая 1), адаптивная оценка ано-оценкой, вычисленной с помощью стационар- малии 5да (кривая 2), стационарная оценка аномалии 6д8 ного алгоритма фильтрации [2]. В качестве па- (кривая 3)
раметров стационарного алгоритма выступали параметры аномалии для равнинного участка. На рис. 2 представлены фрагменты графиков оценок (5да — адаптивная оценка, 5д3 — оценка, полученная стационарным алгоритмом) на фоне истинной аномалии 5д, соответствующие первому гористому участку, показанному на рис. 1. Из рис. 2 видно, что неточность значений параметров для стационарного алгоритма привела к эффекту пересглаженности оценки [2]. В свою очередь адаптивный алгоритм определил моменты смены стационарных состояний и подстроился под изменение структуры аномалии, что позволило получить более точную оценку.
СПИСОК ЛИТЕРАТУРЫ
1. Болотин Ю.В., Голован А.А., Кручинин П.А., Парусников Н.А., Тихомиров В.В., Трубников С.А. Задача авиационной гравиметрии. Некоторые результаты испытаний // Вестн. Моск. ун-та. Матем. Механ. 1999. № 2. 36-41.
2. Голован А.А., Болотин Ю.В., Парусников H.A. Уравнения аэрогравиметрии. Алгоритмы и результаты испытаний. М.: Изд-во ЦПИ при механико-математическом факультете МГУ, 2002.
3. Болотин Ю.В., Попеленскии М.Ю. Анализ точности решения задачи авиагравиметрии при идентификации параметров гравиметра в полете // Фунд. и прикл. матем. 2005. 11, вып. 7. 167-180.
4. Хемминг Р.В. Цифровые фильтры. М.: Недра, 1987.
5. Vaseghi S. Advanced Digital Signal Processing and Noise Reduction. 3d ed. Hoboken: John Wiley & Sons Ltd., 2006.
6. Беллман Р. Динамическое программирование. М.: ИЛ, 1960.
Поступила в редакцию 18.06.2010
