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

УДК 519.676

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

К.А. РЫБАКОВ

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

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

ВВЕДЕНИЕ

Ранее в работах [5], [6] были предложены новые методы и алгоритмы решения задачи оптимальной фильтрации сигналов в нелинейных стохастических системах диффузионного типа, т.е. в системах, модели объекта наблюдения и измерителя для которых описываются с помощью стохастических дифференциальных уравнений. Предложенные методы и алгоритмы базируются на общности структуры уравнений оптимальной нелинейной фильтрации для ненормированной апостериорной плотности вероятности вектора состояния объекта наблюдения, а именно уравнений типа Дункана-Мортенсена-Закаи при фиксированных измерениях, и обобщенного уравнения Фоккера-Планка-Колмогорова, включающего дополнительные слагаемые - функции поглощения и восстановления [1], [3], [4].

Один из основных результатов работ [5], [6] - это обоснование того, что задачу оптимальной фильтрации сигналов можно свести к задаче статистического анализа вспомогательной стохастической системы со случайными обрывами и ветвлениями траекторий, решение которой можно получить методом статистических испытаний или методом Монте-Карло. Есть два варианта подобных методов получения оптимальной оценки вектора состояния по результатам измерений: первый основан на классическом, а второй - на робастном уравнении Дункана-Мортенсена-Закаи. Оба варианта методов апробированы на решении модельных задач, они имеют свои достоинства и недостатки, в некотором смысле дополняя друг друга.

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

1 Работа выполнена при финансовой поддержке РФФИ (проект № 13-08-00323-а).

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

1. ПОСТАНОВКА ЗАДАЧИ

Как и в [8], будем рассматривать модель объекта наблюдения, описываемую стохастическим дифференциальным уравнением Ито с пуассоновской составляющей:

dX

где Xе Я" - вектор состояния; t е Т, Т = t1] - отрезок времени функционирования системы; /^, х): Т х Я" ^ Я" - вектор-функция " х1, о^, х) : Т х Я" ^ Я"Х!1 - матричная функция " х ^ ; Ж^) - ^ -мерный стандартный винеровский процесс, не зависящий от начального состояния Х0, которое определяется плотностью вероятности ф0(х). Далее, Q(t) - общий пуассоновский процесс, заданный в форме

р

Q(t) = Е А *.

к=1

В правой части последнего равенства Р^) - пуассоновский процесс, Лк - независимые случайные векторы из Я", распределение которых задано плотностью вероятности у(тк,Л), т.е. вектор состояния X получает случайные приращения в моменты времени т1,Т2,...е Т, образующие пуассоновский поток событий:

XО*) = X

Если величина приращения зависит от вектора состояния, то используется условная плотность вероятности у(тк, Л | £), характеризующая распределение Лк при условии X (тк - 0) = £. В частном случае у(тк, Л | £) = ^(хк, Л). Наряду с у(тк, Л | £) введем плотность вероятности п(тк, х | £), характеризующую распределение X(тк) при условии X(Тк - 0) = Е,.

Пуассоновский поток событий и, следовательно, моменты времени Т1, Т2,..., а также пуассоновский процесс Р^) определяются интенсивностью , х), т.е. условная вероятность события (2) при X ^) = х на промежутке ^, t + Лt ] задается равенством

Р^, t + Лt) = Pr(P(t + Л0 - P(t) = 11 X(t) = х) = Щ х^ + о(Л).

Модель измерительной системы проще, чем в [8], а именно здесь, как и в [6], будем рассматривать стационарную модель:

dY ^) = с

7

где У, 7 е Ит - векторы измерений; с(х): И" ^ К" - вектор-функция т х1, V(1) - т -мерный

стандартный винеровский процесс, не зависящий от Ж(1) и от начального состояния Х0, N(1) - т -

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

Задача оптимальной фильтрации состоит в нахождении оценки X(1) траектории случайного процесса X(1) по результатам измерений У0' = {У(т), те [10,1)} или 70 = {7(т), те [10,1)} в соответствии с заданным критерием оптимальности. Далее при записи соотношений будут использованы измерения У(1) или 7(1) в зависимости от удобства. Задачу оптимальной фильтрации можно рассматривать шире: как нахождение апостериорной плотности вероятности р(1, х | Уд ) вектора состояния X.

При использовании критерия минимума среднеквадратической ошибки оценивания имеем [10]:

Х

2. УРАВНЕНИЯ ДЛЯ НЕНОРМИРОВАННЫХ АПОСТЕРИОРНЫХ ПЛОТНОСТЕЙ ВЕРОЯТНОСТИ

Приведем уравнение Дункана-Мортенсена-Закаи [3], [4], [10] для ненормированной апостериорной плотности вероятности ф(1, х | У0\), которая связана с функцией р(1, х | У0) соотношением

р«ж) =, ф

.И ф( и х | У0) -х

и воспользуемся некоторыми результатами, полученными в [8].

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

Эф(\х|У°) = £ф(', х | У0) -х)ф(1, х | П) + К х | Е, | У)-% +

т -У (Г) + Е С а (х) -й-1 ф(1, х | У0), ф(^0, х) = ф0( х),

а=1 Ш

где

1 т

х | У) = А<&, х | У0) - - £ с2( х)ф(1, х | У0),

2 а=1

Лф(г, х | У0) = -£ ^ [ / (1, х)ф(1, х|У0) ] +1£ £ -Д— [ % (1, х)ф(1, х|У0 ) ],

1=1 дхг 2 1=1 ,=1 дхг

% (1, х) = о( х, х)от( х, х),

ИЛИ

где

Эф0= ^,х 170) - Ч*, х)ф(*, * | )+^ Ч(*, , г | $)ф(/, $ | 70 ж -

-д-

ш

- (х г) = x, ^ x, 2) < 0, + (х г) = Иx, ^ x, 2) > 0,

д (х 7) [0, д(х, 7) > 0, д (x, [0, д(х, 7) < 0;

д(х, 7) = £ С„ (х) ( га- ^).

Уравнение (8) по структуре соответствует обобщенному уравнению Фоккера-Планка-Колмогорова, при этом функции д- (X(*), 7(*)) и д+ (X(*), 7 (*)) - интенсивности обрывов и ветвлений траекторий случайного процесса X(*), а произведения д- (х, 7 (*))ф(*, х 170) и д+ (х, 7(*))ф(*, х | 70) - функции поглощения и восстановления соответственно [1]. Следовательно, условные вероятности обрывов и ветвлений при X(*) = х и 7 (*) = г на промежутке [*, * + А*] определяются равенствами:

Р - (*, * + А*) = д- (х, £)Аг + о( А), Р + (*, * + А*) = д+ (х, г)А* + о(А*).

Таким образом, функции , х | 70) и ф(*, х | 70 ) характеризуют распределение вектора

X - состояния объекта наблюдения, описываемого уравнением (1), - с учетом того, что траектории случайного процесса X(*) получают случайные приращения, ветвятся или обрываются. Напомним [8], что перечисленные события образуют неоднородные пуассоновские потоки с известными интенсивностями, при ветвлении в фиксированный момент времени может появиться только одна новая ветвь, при обрыве прекращается моделирование только одной ветви. Для удобства моделирования, как и в случае систем диффузионного типа [5; 6], каждая из новых ветвей должна рассматриваться как самостоятельная траектория.

Однако, как и для систем диффузионного типа, уравнение Дункана-Мортенсена-Закаи содержит множитель типа белого шума - случайный процесс N(*), наличие которого усложняет его решение с помощью метода статистических испытаний. Основная сложность состоит в том, что N(t) входит в интенсивности обрывов и ветвлений траекторий случайного процесса X (*), делая их, вообще говоря, неограниченными (см. (4) и (8)). В результате, как показано на примере в работе [6], для интенсивностей

д-

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

Чтобы минимизировать влияние указанного недостатка, проведем замену ненормированной апостериорной плотности вероятности [9]:

ф(*, х|70) = е7 )с(х) р(*, х|70).

Тогда[8]

ЭР

т т

27«(04Р(*,х | 70) + ^^7„ОД(0^,х | 70*).

а=1 в=1

а=1

В уравнении (10) Са = [Са, С], = 1[Са , £]], где [Са, С] и [Са, ^ ] -коммутаторы указанных в скобках операторов, Са - операторы умножения на функции са (х), а, в = 1,2,...,т.

Начальное условие для уравнения (10) совпадает с начальным условием для уравнения (7) или (8), поскольку е7Т(*)с(х) = 1 при * = *0, т.е. р(*0, х) = ф0(х).

Уравнение (8) отличается от робастного уравнения Дункана-Мортенсена-Закаи для систем диффузионного типа наличием слагаемых

, х)р(*, х|70) + | , $)п(, х|$)е

7 Т

р(*, $ | 70)Ш$,

поэтому можно воспользоваться результатами работы [6], где подробно показано, как получить другую форму записи робастного уравнения Дункана-Мортенсена-Закаи на основе тождественных преобразований выражения

т т т

Ср(*, х | 70 ) - 2 7„ (*)4р(*, х | 70 ) + 7« (*)7р (*)/>(*, х 170*),

а=1 в=1

которое можно представить как

Ар(*, х 1700)+х, 7(*))р(*, х 170),

где

Др(*, х 170*) = А[I(*, х, 7(*))р(*, х | 70)] +12£Ч^Т[х)р(*, х | 70)],

1=1 дх1 2 I=1 у=1 дх1 ах1

!

Эс( х) дх

у =1

Т

7 (*),

, х, 7

1 сТ

дх 2

2

дс( х) дх

g

дс( х) дх

7 (*).

Таким образом, имеем

Эр

+ v(t, x, Y(t))p(t, x | Y0), p(t0, x) = ф0(x).

3. СВЕДЕНИЕ К ЗАДАЧЕ АНАЛИЗА СИСТЕМ С ОБРЫВАМИ И ВЕТВЛЕНИЯМИ ТРАЕКТОРИЙ

Применим подход, использованный ранее для систем диффузионного типа [6], а именно сведем задачу оптимальной фильтрации к задаче анализа вспомогательной стохастической системы с обрывами и ветвлениями траекторий. Чтобы правильно интерпретировать входящие в уравнение (11) слагаемые при формировании алгоритма моделирования траекторий вспомогательного случайного процесса с обрывами, разрывами и ветвлениями, во-первых, представим функцию v(t, x, y) в виде:

v(t, x, y) = -v- (t, x, y) + v+(t, x, y),

где

, f-v(t,x, y), v(t,x, y) < 0, . fv(t,x, y), v(t,x, y) > 0,

v x y) _ \n „ ... v (t, x, y) = \n „

10, v(t, x, y) > 0, 10, v(t, x, y) < 0,

а во-вторых, определим функции

X

T, n(t V | ?)eyT((5)-C(x)) yT((5)-(x)) dx nt x I 5 __n(t, x | 5)e_

j" n(t, x | 5)eyT ((5)-c(x)) dx Тогда

dp

_ Дp(t, x | Y0) - X(t, x)p(t, x | Y0) + J X(t, 5, Y(t))n(t, x | 5, Y(t))p(t, 5 | Y0)d5 -

dt ' ' 07 v' /rv ' ' 07 У I" W/r-V.-SI

- v- (t, x, Y(t))p(t, x | Y0) + v+ (t, x, Y(t))p(t, x | Y0), p(t0, x) _ Ф0 (x),

где

£(*,х,у) > 0, П|(Х,х | у) > 0 и Гт^,х | у)йх = 1, Хе Т, х,^е Я", уе Ят,

•»Я

т.е. после преобразований робастное уравнение Дункана-Мортенсена-Закаи представлено так, что его структура совпадает со структурой уравнения (8).

Это уравнение описывает закон изменения ненормированной апостериорной плотности вероятности р(Х, х | У0Х) вектора состояния X = X(X) в момент времени X для вспомогательной

стохастической системы при условии, что вектор сноса - /(X, х, У(X)), а матрица диффузии совпадает с матрицей диффузии исходного объекта наблюдения (1) - g(X, х). Начальное состояние вспомогательной стохастической системы совпадает с начальным состоянием Х0

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

1) обрыв траектории с суммарной интенсивностью , X(*)) + у- (*, Х(*), 7(*)) , т.е. вероятность обрыва траектории при X(*) = х и 7(*) = у задается выражением

Р- (*, * + А*) = (Ч(*, х) + у- (*, х, у)) А* + о(А*);

2) ветвление траектории с интенсивностью (*,X(*),7(*)) , т.е. вероятность ветвления траектории при X(*) = х и 7(*) = у задается выражением

Р + (*, * + А) = (*, х, у) А + о(А*);

3) ветвление траектории со скачком с интенсивностью Ч(*, X(*),7(*)), т.е. вероятность такого события при X(*) = х и 7(*) = у задается выражением

Р*(*, * + А*) = Ч (*, х, у) А + о(А*),

а величина XX(тк) при ветвлении со скачком для новой ветви определяется условной плотностью вероятности П(тк, х | $, у) , если X(тк - 0) = $, 7(тк ) = у, Тк - момент ветвления. Выражение

(,х) + у-(*,х,7(*))),х | 70) определяет функцию поглощения, а выражение

Ч

- функцию восстановления траекторий случайного процесса X(*) [1].

Для этой системы распределением моментов времени появления обрывов и ветвлений управляет процесс 7(*), а не 7 (*), как это происходит со вспомогательной стохастической системой, построенной на основе классического, а не робастного уравнения Дункана-Мортенсена-Закаи. Этот процесс влияет и на поведение траекторий между моментами ветвлений, и на интенсивность ветвлений со скачком, и на распределение величины этого скачка. Основываясь на результатах моделирования, можно утверждать, что интенсивности

V-

меняются значительно медленнее по сравнению с интенсивностями

д- (О = д- 7 (*)) и д+ (*) = д+ 7(*)),

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

«баланса» между поглощением и восстановлением траекторий при разрыве траектории. Если для исходной стохастической системы интенсивности, входящие в функции поглощения и восстановления, одинаковы - они совпадают с интенсивностью X(t, х) разрыва траекторий случайного процесса X(t), то для вспомогательной стохастической системы они разные. Даже если в исходной постановке задачи X(t, х) = X = const, то для случайного процесса X(t) при восстановлении интенсивность будет переменной по времени и будет определяться текущим состоянием, поскольку она должна вычисляться на траекториях X(t) и Y(t). Эта интенсивность выражается через коэффициент сноса с(х) в уравнении (3) или (4) модели измерительной системы и явным образом зависит от Y(t).

Ненормированную апостериорную плотность вероятности p(t, х | Y0t) можно найти приближенно, моделируя траектории вспомогательного случайного процесса X(t) с учетом обрывов, разрывов и ветвлений. Как и в аналогичных алгоритмах, разработанных для более простых моделей систем наблюдения [5], [6], по ансамблю траекторий, полученному в результате применения методов численного решения стохастических дифференциальных уравнений и моделирования неоднородных пуассоновских потоков событий [2], [3], можно оценить функцию p(t, х | Y0t), например, с помощью построения гистограммы, и по этой оценке

найти апостериорную плотность вероятности p(t, х | Y0t) согласно формулам (6) и (9). По известной функции p(t, х | Y0t) может быть найдена и оптимальная оценка X(t), причем можно

использовать не только выражение (5), если модифицировать постановку задачи оптимальной фильтрации, заменив критерий минимума среднеквадратической ошибки оценивания каким-либо другим критерием [3].

Как отмечалось при формулировании постановки задачи, стационарность модели измерительной системы не имеет принципиального значения и предлагаемая методика решения задачи оптимальной фильтрации может применяться и для нестационарных моделей. При зависимости от времени коэффициента сноса в уравнении (3) или (4) модели измерительной системы в робастном уравнении Дункана-Мортенсена-Закаи появится дополнительное слагаемое, которое будет относиться к интенсивностям обрывов и ветвлений траекторий вспомогательного случайного процесса [7]. Кроме того, можно рассматривать вариант задачи оптимальной фильтрации, когда размерности вектора измерений и вектора шума в уравнении измерительной системы не совпадают, а в случае совпадения коэффициент при шуме необязательно представляет собой единичную матрицу и может быть матричной функцией времени, но не функцией вектора состояния [3], [4].

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

ЛИТЕРАТУРА

1. Казаков И.Е., Артемьев В.М., Бухалев В.А. Анализ систем случайной структуры. - М.: Физматлит, 1993.

2. Михайлов Г.А., Аверина Т.А. Алгоритм «максимального сечения» в методе Монте-Карло // Доклады АН. 2009. Т. 428. № 2. С. 163-165.

3. Пантелеев A.B., Руденко Е.А., Бортаковский A.C. Нелинейные системы управления: описание, анализ и синтез. - М.: Вузовская книга, 2008.

4. Параев Ю.И. Введение в статистическую динамику процессов управления и фильтрации. - М.: Советское радио, 1976.

5. Рыбаков К.А. Сведение задачи нелинейной фильтрации к задаче анализа стохастических систем с обрывами и ветвлениями траекторий // Дифференциальные уравнения и процессы управления. 2012. № 3. С. 91-110. [Электронный ресурс]. URL: http://www.math.spbu.ru/diffjournal.

6. Рыбаков К.А. Приближенное решение задачи оптимальной нелинейной фильтрации для стохастических дифференциальных систем методом статистических испытаний // Сибирский журнал вычислительной математики. 2013. Т. 16. № 4. С. 377-391.

7. Рыбаков К.А. О решении робастного уравнения Дункана-Мортенсена-Закаи для нестационарных систем // Информационные и телекоммуникационные технологии. 2014. № 22. С. 9-15.

8. Рыбаков К.А. Приближенный метод фильтрации сигналов в стохастических системах диффузионно-скачкообразного типа // Научный вестник МГТУ ГА. 2014. № 207. С. 54-60.

9. Рыбаков К.А. Решение робастного уравнения Дункана-Мортенсена-Закаи для систем диффузионно-скачкообразного типа на основе спектрального метода // Системи обробки шформацп. 2014. Вып. 7 (123). С. 143-147.

10. Сшитым И.Н. Фильтры Калмана и Пугачева. - М.: Логос, 2007.

FILTERING FOR JUMP-DIFFUSION MODELS BY STATISTICAL MODELING METHOD

Rybakov K.A.

This article develops new methods that reduce the optimal filtering problem for jump-diffusion models to the analysis problem for the special stochastic system with jumps, branching and terminating trajectories. Earlier appropriate methods and algorithms have been proposed and tested for diffusion models.

Keywords: branching processes, conditional density, Duncan-Mortensen-Zakai equation, jump-diffusion, Monte Carlo method, optimal filtering problem, statistical modeling, stochastic system.

REFERENCES

1. Kazakov I.Ye., Artem'ev V.M., Bukhalev V.A. Analiz sistem sluchaynoy struktury (Analysis of Systems with Random Structure), Moscow, Fizmatlit Publishing Company, Nauka Publishers, 1993.

2. Mikhaylov G.A., Averina T.A. The Maximal Section Algorithm in the Monte Carlo Method, Doklady Mathematics, 2009, vol. 80, no. 2, pp. 671-673.

3. Panteleev A.V., Rudenko Ye.A., Bortakovskiy A.S. Nelineynye sistemy upravleniya: opisanie, analiz i sintez (Nonlinear Control Systems: Description, Analysis, and Synthesis), Moscow, University Book, 2008.

4. Paraev Yu.I. Vvedenie v statisticheskuyu dinamiku protsessov upravleniya i fil'tratsii (Introduction to Statistical Dynamics of Control and Filtering Processes), Moscow, Soviet Radio, 1976.

5. Rybakov K.A. Differentsial'nye uravneniya i protsessy upravleniya, 2012, no. 3, pp. 91-110, available at: http://www.math.spbu.ru/diffjournal.

6. Rybakov K.A. Solving approximately an optimal nonlinear filtering problem for stochastic differential systems by statistical modeling, Numerical Analysis and Applications, 2013, vol. 6, no. 4, pp. 324-336.

7. Rybakov K.A. Informatsionnye i telekommunikatsionnye tekhnologii, 2014, no. 22, pp. 9-15.

8. Rybakov K.A. Nauchnyi vestnikMGTU GA, 2014, no. 207, pp. 54-60.

9. Rybakov K.A. Sistemy obrabotki informatsii, 2014, no. 7 (123), pp. 143-147.

10. Sinitsyn I.N. Fil'try Kalmana i Pugacheva (Kalman and Pugachev Filters), Moscow, Logos, 2007.

Сведения об авторе

Рыбаков Константин Александрович, 1979 г.р., окончил МАИ (2002), кандидат физико-математических наук, доцент кафедры математической кибернетики факультета «Прикладная математика и физика» МАИ, автор 110 научных работ, область научных интересов - анализ и синтез стохастических систем управления, спектральная форма математического описания систем управления, методы моделирования систем управления.