Фильтрация сигналов со скачками, возникающими в дискретном времени и с конечным горизонтомфильтрация сигналов со скачками, возникающими в дискретном времени и с конечным горизонтом Текст научной статьи по специальности «Математика»

РАДИОФИЗИКА

УДК 519.21

Г.И. Белявский, И.В. Мисюра

Южный федеральный университет

фильтрация сигналов со скачками, возникающими в дискретном времени и с конечным горизонтом

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

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

Введение

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

с1X = а+ р^Щ, I е [0, Т], Х0 = 0,

где — стандартное броуновское движение; а,, Р, — случайные коэффициенты (предсказуемые либо детерминированные), которые удовлетворяют условию Липшица по временной переменной t почти всюду.

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

В современных исследованиях значительное внимание уделяется процессам со скачками (процессы Леви, аддитивные и другие), и для ознакомления с состоянием вопроса следует обратить внимание на монографии [3, 4].

Для описания процесса со скачками целесообразно рассмотреть уравнение вида

¿Х1 = а, Л + р , с1Ж1 + у ¿Ср,

в котором СР1 — составной процесс Пуассона, не зависящий от ^ (т. е. от стандартного броуновского движения). При этом на любом конечном интервале [0, Т] число скачков данного процесса конечно. Скачки происходят в случайные моменты времени т1, т2, ..., хк; приращения Дт 1 — это независимые и одинаково распределенные показательные случайные величины с интенсивностью X, поэтому дифференциал составного процесса Пуассона выражается как

где — независимые, одинаково распределенные случайные величины.

Фильтрация процессов со скачками упоминается также в монографии [5]. По аналогии с этим упоминанием рассмотрим дискретный аналог диффузии со скачками для равномерного разбиения интервала с шагом к:

Р( 1)

Хк

ДХ; = а-к + + у ,_1 £ (1)

I=1

где Pfh) — последовательность независимых пуассоновских случайных величин с общей интенсивностью Ah; ) — последовательность всех прочих независимых случайных величин с общим законом распределения; приращение AWj распределено по нормальному закону N (0, h).

Для 1) е N(0,1) уравнение (1) приобретает более простой вид:

AXj = a i_ih + (в j_iJh + Y i ^ Щ, (2)

где M. — независимые стандартные нормальные случайные величины.

Уравнение (2) является базовым для описания процессов, рассмотренных далее.

Под сигналом мы будем далее понимать случайный вектор X = (Xt)", и модели поведения сигнала будем описывать в конечномерном линейном пространстве Rn. С сигналом будем связывать естественный стохастический базис

(0,(F, )0, Fn),

где Fo =а{0,0}, F} = B{R1}.

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

Возможны две постановки задачи: априорная и апостериорная. в первом случае рассматривается семейство конечномерных условных законов распределения Law (X1 / Y1 ) и последовательность оценок X1 вектора XJ по наблюдаемым векторам Y1 (верхний индекс обозначает проекцию соответствующего вектора на подпространство R1). В качестве наилучшей считается, как правило, оценка, доставляющая минимум среднеквадратического отклонения E(X1 - X1 )2, при условии, что оценка XJ е ct(Y1 ). Известно, что наилучшая в среднеквадратическом смысле оценка — это условное математическое ожидание:

X1 = E(X1 / Y1). (3)

При существовании условной плотности f (х / YJ), кроме условного математического ожидания, используется также

оценка максимального правдоподобия.

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

При апостериорной постановке задачи вектор У доступен полностью и вычисляется одна апостериорная оценка сигнала, например

X = Е(Х / У).

Общая модель описания зашумленных процессов

Предположим, что в уравнении (2) коэффициенты являются константами и, не нарушая общности, будем считать, что а . = 0. Кроме того, будем считать, что интенсивность скачков или шаг разбиения настолько малы, что вероятностью Р(Р^) > 2) можно пренебречь. Благодаря этим предположениям, уравнение (2) приобретает более компактный вид:

ДХ. = (Р>/А + у8 j. )М.,

где 8. £ {0, 1}.

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

£ {ст2, ст2 + с1} (^ > 0).

Если = ст2, то в 1-й момент времени скачка нет; если же = ст2 + то в /-й момент времени скачок есть. в результате модель задается системой стохастических разностных уравнений:

У = X + , АХ = БМ; (4)

при этом векторы N М независимы и распределены по нормальному закону с нулевым математическим ожиданием и единичной ковариационной матрицей; А — матрица дискретного оператора дифференцирования.

Рис. 1. Траектория сигнала со скачками в моменты времени т1 = 12 о.е., т2 = 34 о.е., т3 = 67 о.е.; значения параметров: а1 = 0,2, а2 = 0,1, с1 = 1,0

На рис. 1 представлена траектория сигнала с тремя скачками в определенные моменты времени и соответствующие значения параметров о и ¿.

Модель 1

В данной модели предполагается, что число скачков ограничено небольшим значением К. Вид случайной матрицы Б полностью определяется поведением К-мерной случайной величины I = (11,12, ..., ¡к), для которой ¡1 е {0, 1, ..., п}. Если = (3,, то ¡1+1 = ¡1+2 = ... = 0; если ¡1 ф (3, то ¡1+1 = 0 или ¡1 < ¡1+1. Вычисление оценки связано с решением системы линейных алгебраических уравнений (СЛАУ), имеющей матрицу, которая отличается от матрицы с известным сингулярным разложением на матрицу ранга к. Поскольку значение к не оказывает существенного влияния на приведенные далее вычислительные процедуры, мы остановимся на случае к = 1 (чтобы не загромождать статью выкладками).

Рассмотрим случайную величину ¡ е {0, 1, ..., п} с законом распределения Р (! = I) = р1. Равенство ¡ = 0 означает, что все = ст2; ¡ = 1 ф 0 означает, что

d =

i ф j;

ст2 + d, i - j.

Для модели 1 фильтрация рассматривается при двух априорных предположениях.

Предположение первое. Закон распределения случайной величины ¡ известен. В этом случае формула (3) приобретает вид

1 -£ E(X / Y, l)p.

(5)

Условный закон распределения Law (X/Y, l) является нормальным, с математическим ожиданием

XеI -Ь

(

-1 I + ATD-2 A

VCTi

Y

Y.

(6)

В результате оценка вектора X имеет

вид

Xе -I X'pt .

(7)

Оценка (7) является линейной оценкой и вместе с удалением шума из сигнала сглаживает его резкие изменения, что не всегда оправдано, например, при выделении границ (edge detection). В отличие от (7), оценка

X - Xj, (8)

где

~(Y - X1 )2

j = arg min

2ст2

-ln p

относится к нелинейным оценкам.

Кроме того, оценка (8) позволяет получить ориентировочное значение момента разладки случайного процесса, если под

l-0

-1

разладкой понимать разрыв траектории. введем обозначения

И, = -1I + АТБ— А. СТ1

Нетрудно установить связь между матрицами Иг (/ ф 0) и матрицей И0; эта связь выражается следующим образом:

И, = И0 + ааТ, (9)

где а, — вектор,

^ст\ + (ст2 + с1 )2 ст21

Ае,

(е, — /-й орт);

матрица И0 = -1 I + -1 АТА.

СТ1 СТ2

Равенство (8) позволяет выразить оценки Х1 через оценку Х0 :

X' = Х0 - , {ИпХ) ) И-1а,. (10)

1 + (И- а,, а,)

Основная вычислительная трудность связана с обращением матрицы И0. в работе [6] представлено спектральное разложение матрицы

И0 = Ш ЛШТ,

,. = ЯП

518'п Iп

1 - С08

п + 1

где

Это разложение затем использовалось в работе [7]. Тогда обратная матрица имеет вид И0-1 = Ш Л-1ШТ.

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

тт| Щ, (11)

X еЬ 11 11

где Ь — выпуклая оболочка, натянутая на векторы: У - Х,.

оценка сигнала имеет вид

Х = У - г*, (12)

где 2* — решение задачи (11).

одновременно с этим находится априорное распределение случайной величины ,

а/ =

Рис. 2. Оценка сигнала, выраженного формулой (12); использован эмпирический байесовский подход

что позволяет решить задачу разладки случайного процесса. Задача (11) хорошо изучена [9], и для ее решения существуют эффективные алгоритмы.

На рис. 2 представлена оценка сигнала (12) по его зашумленной реализации.

Модель 2

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

Б( Е) = Ша§( ст2 + Цй).

Условная оценка имеет вид

- 1 Г1 -1

X(Е) = -1 -11 + ЛТБ(Е)-\ СТ1 1ст1

Y, (13)

(14)

а оценка сигнала —

X = ЕХ (Е).

Далее будем считать, что Е — однородная марковская последовательность с двумя состояниями и матрицей переходных вероятностей

Г ? 1 - ? " 1 0

Единица во второй строке означает, что скачки не следуют один за другим. Это ограничение несущественно для вычисления оценки. Для вычисления по формуле (14) используем метод Монте-Карло:

Q =

1 N _

x « — у X (L)),

(15)

где Е1) — сгенерированная марковская последовательность.

Этот метод фильтрации состоит из двух элементов: генерации бинарной однородной марковской последовательности и решения СЛАУ с трехдиагональной матрицей. Каждый из элементов не требует больших вычислительных затрат.

Как и в модели 1, рассмотрим оценку, которая получается в результате решения задачи:

min У

X,L 4-f

(Y - X' )2 + (X' - X'-1 )2

2 ст2

2( ст2 + Lid )2

+ 1п(ст2 + Ltd) + Li ln1—q +

+Li Li-1 ln q + L(-1 ln1

(16)

1 - q

причем X,, = Е0 = 0. Решение задачи (16) может быть выполнено методом динамического программирования.

Уравнение Беллмана

Для вывода указанного уравнения, связанного с фильтрацией сигнала, рассмотрим последовательность функций:

" ъ - X )2,

ф^ (Xk , Lk) =

+ (X,. - X,.-1)2 2( ст2 + L.d )2

min

..,Xn ;Lk+1

У

2ст2

+ ln(a2 + Ltd) + L, ln

1 - q

+LiLi-1+ Li-1ln1 1 - q q

Уравнение Беллмана имеет следующий вид:

Фк-1(xk-l, 4-1) = mL1

Y - Xk )2

2ст2

(Xk - Xk-1)2 2( ст2 + Lkd )2

+ ln(a2 + Lkd) + Lk ln

1 - q

+44-1 ^ + Lk-1 1 + Фk (Xk, Lk) 1 - q q

где Xk (Xk-1, Lk-1), Lk (Xk-1, Lk-1) - аргументы, на которых достигается минимум.

Краевое условие выражается как

Ф-Ф!, 4+1) = 0.

Непосредственное решение уравнения Беллмана, связанное с вычислением последовательности функций, требует больших вычислительных затрат. Поэтому в качестве альтернативного метода рассмотрим спуск по обобщенным координатам (X, L). Метод спуска - итерационный, на каждой итерации которого решаются две вспомогательные задачи:

min F (X, L));

min F (X(t+1), L), где F - целевая функция в задаче (16).

(17)

(18)

Вектор Х('+1) вычисляется по формуле (13) с использованием вектора Ь). Вектор Ь('+1) — решение задачи (18), которое находится методом динамического программирования. Последний применяется для оценки (16) при условии, что вектор Х фиксирован, т. е. Х = Х('+1). Итерации продолжаются до тех пор, пока не выполнится равенство Ь+1) = Ь(').

на рис. 3 представлен результат фильтрации, выполненной с помощью итерационного алгоритма.

некоторым обоснованием алгоритма является следующее утверждение.

Утверждение. Итерационный алгоритм останавливается за конечное число шагов.

Действительно, по переменной Х целевая функция является квадратичной, с положительно определенной матрицей, диагональ которой зависит от Ь. Отсюда следует, что на каждой итерации происходит уменьшение целевой функции. Поскольку множество значений переменной Ь — конечно, то алгоритм останавливается за конечное число шагов.

Двухкритериальная оптимизационная задача фильтрации

В связи с описанием и использованием модели 2 рассмотрим еще один подход к фильтрации сигналов со скачками. Опре-

делим два критерия: критерий близости и критерий гладкости.

критерий близости выразим в виде Х) = (У - Х )2, критерий гладкости —

Ш2(Х) = (АТАХ, Х).

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

Математическое ожидание первого критерия

ЕШ (Х) = пст2,

а второго —

ЕШ2(Х) = пст2 + (2ст2С + С2) х

х т - еп )(1 - е )-1 бр е2),

где е1 = (1, 0), е2 = (1, 0).

В результате для оценки необходимо найти

шт

Х

(У - Х)2 + (АТАХ, Х)

пст

\ = пст2 + (2ст2С + С2) х

х е( I - еп)(I - е) "Ч, е2],

что приводит к уравнению

(АТА + 91 )Х = У,

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

где

9 = п-1ст-2[пст2 + (2ст2й + й2) х

х (е (I - оп)(/ - е )-1 ¿1,«!)].

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

min(ЛTЛX, X), (19)

при ограничении (У - X)2 < а.

В отличие от предыдущего случая, решение задачи (19) сводится к решению двух уравнений:

(ЛГЛ + 9/)Х = У; (У - X)2 = а. (20)

Используем спектральное представление матрицы

ЛТЛ + 91 = иТ Л(9)и.

Ортогональная матрица и — та же, что и выше, а диагональные элементы матрицы Л(9) выражаются как

(

X (9) -9 + 2

1 - cos

j П и + 1

В результате уравнения (20) преобразуются к следующему виду:

X - UTA-1(9)Z,; £ ./-1

1-

X ((0)

Y2 - а, (21)

где У = иУ.

Данная оценка является нелинейной. При варьировании величины а получаются как различные оценки X (а), так и различные значения критерия ¥г( X (а)). Параметр а можно связать с математическим ожиданием Е¥1( X), например, с помощью равенства а = кпи^. Константу к можно выбрать из условия Р((У - X)2 < а) < п.

Заключение

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

Работа выполнена при финансовой поддержке РФФИ, проект № 140100579.

список литературы

1. Липцер Р.Ш., Ширяев А.Н. Статистика случайных процессов. М.: Наука, 1974. 696 с.

2. Липцер Р.Ш., Ширяев А.Н. Нелинейная фильтрация диффузионных марковских процессов. Исследования по математической статистике // Труды МИАН СССР. 1968. Т. 104. С. 135-180.

3. Applebaum D. Levy processes and stochastic calculus. Cambridge: Cambridge University Press, 2004. 384 p.

4. Bertoin J. Levy processes. Cambridge: Cambridge university Press, 1996. 266 p.

5. Oksendal B., Sulem A. Applied stochastic control of jump diffusion. New York: Springer Ver-

lag, 2004. 208 p.

6. Беллман Р., Энджел Э. Динамическое программирование и уравнения в частных производных. М.: Мир, 1974. 204 c.

7. Мисюра И.В. Один метод фильтрации случайного сигнала // Обозрение прикладной и промышленной математики. 2010. Т. 17. Вып. 6. С. 911-912.

8. Белицер Э., Еникеева Ф.Н. Адаптивная фильтрация случайного сигнала в гауссовском белом шуме // Проблемы передачи информации. 2008. Т. 44. Вып. 4. С. 39-51.

9. Демьянов В.Ф., Малоземов В.Ф. Введение в минимакс. М.: Наука, 1972. 368 с.

сведения об авторах

БЕЛЯВСКИй Григорий Исаакович — кандидат физико-математических наук, профессор кафедры высшей математики и исследования операций Южного федерального университета. 344090, Россия, г. Ростов-на-Дону, ул. Мильчакова, 8а. [email protected]

МИСЮРА Илья Владимирович — аспирант кафедры высшей математики и исследования операций Южного федерального университета.

344090, Россия, г. Ростов-на-Дону, ул. Мильчакова, 8а. [email protected]

Beliavsky G.I., Misyura I.V. SIGNAL FILTERING WITH JUMPS DURING DISCRETE TIME AND UNDER FINITE HORIZON.

The problem of filtering signals with jumps occurring at random times on a background of white noise in discrete time and with finite horizon has been considered. Two signal models being discrete analogs of diffusion with jumps were developed.

signal filtering, jump diffusion, markov chain, bellman equations, discrete time, finite horizon.

references

1. Liptser R.Sh., Shiryaev A.N. Statistika sluchaynykh protsessov. Moscow: Nauka, 1974, 696 p. (rus)

2. Liptser R.Sh., Shiryaev A.N. Nelineynaya fil'tratsiya diffuzionnykh markovskikh protsessov. Issledovaniya po matematicheskoy statistike. Trudy MIAN SSSR, 1968, Vol. 104, pp. 135-180. (rus)

3. Applebaum D. Levy processes and stochastic calculus. Cambridge: Cambridge University Press, 2004, 384 p.

4. Bertoin J. Levy processes. Cambridge: Cambridge University Press, 1996, 266 p.

5. Oksendal B., Sulem A. Applied stochastic control of jump diffusion. New York: Springer Verlag,

2004, 208 p.

6. Bellman R., Endzhel E. Dinamicheskoe programmirovanie i uravneniya v chastnykh proizvodnykh. Moscow: Mir, 1974, 204 p. (rus)

7. Misyura I.V. Odin metod fil'tratsii sluchaynogo signala. Obozrenie prikladnoy i promyshlennoy matematiki, 2010, Vol. 17, Iss. 6, pp. 911-912. (rus)

8. Belitser E., Enikeeva F.N. Adaptivnaya fil'tratsiya sluchaynogo signala v gaussovskom belom shume. Problemy peredachi informatsii, 2008, Vol. 44, Iss. 4, pp. 39-51. (rus)

9. Dem'yanov V.F., Malozemov V.F. Vvedenie v minimaks. Moscow: Nauka, 1972, 368 p. (rus)

the authors

MISYURA Il'ya V.

Southern Federal University

8a Milchakova St., Rostov-on-Don, 344090, Russia [email protected]

BELYAVSKIY Grigoriy I.

Southern Federal University

8a Milchakova St., Rostov-on-Don, 344090, Russia [email protected]

© Санкт-Петербургский государственный политехнический университет, 2014