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

АДАПТИВНАЯ ПРЕДЕЛЬНО ОПТИМАЛЬНАЯ ФИЛЬТРАЦИЯ ПРИ НЕИЗВЕСТНОЙ КОВАРИАЦИИ ВОЗМУЩЕНИЙ

А. Е. Барабанов1, Д. В. Ромаев2

1. С.-Петербургский государственный университет,

д-р физ.-мат. наук, профессор, andrey.barabanov@gmail.com

2. С.-Петербургский государственный университет, аспирант, romaev@yandex.ru

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

[1], как правило, оказывается лучшим способом оценивания фазового вектора системы даже в случае нарушения некоторых предположений, при которых он оптимален

[2]. Точность оценивания фильтра Калмана—Бьюси существенно зависит от матриц ковариаций возмущений в объекте и шумов измерения. В частности, известно, что любой устойчивый фильтр является оптимальным при подходящем выборе этих матриц. Основной трудностью в синтезе оптимальных фильтров является отсутствие необходимых данных о свойствах полезного сигнала и помехи. Характеристики возмущений обычно неизвестны даже приблизительно. Фильтры, в которых собственные параметры настраиваются в ходе активного процесса оценивания, относятся к адаптивным [3-5]. Если параметры восстанавливаются достаточно эффективно, то со временем фильтр приобретает оптимальные свойства — является предельно оптимальным [2]. В данной статье представлен адаптивный фильтр Калмана—Бьюси, содержащий состоятельные оценки неизвестной матрицы ковариаций возмущений. Оценивание матрицы ковариаций возмущений производится методом стохастической аппроксимации. Результат является обобщением решения задачи предельно оптимальной адаптивной фильтрации [6] на случай векторного возмущения в объекте наблюдения. Методы исследования и доказательства существенно изменены.

1. Постановка задачи. В стандартной задаче Калмана—Бьюси объект наблюдения описывается уравнениями

а(д-1)хг = Ь(д-1)щ + с(д-1)->мг,

Уг = X + V,

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

оценить вектор хг.

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

© А. Е. Барабанов, Д. В. Ромаев, 2011

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

2. Оптимальный стационарный фильтр. Решение задачи оптимального оценивания вектора xг по измерениям увплоть до момента к = Ь включает предварительные преобразования. Матричный многочлен ¿(г) определяется из уравнения факторизации

а(г )Еу а* (г) + с(г)Кю с* (г) = ¿(г)3*(г)

с дополнительным условием: все нули С расположены вне единичного круга, что означает, что С устойчив.

Лемма 1. Оптимальный фильтр в поставленной задаче оценивания по критерию минимума дисперсии ошибки определяется уравнением

Хг = Уг - 7$г,

где последовательность (6г) определяется из уравнения фильтра

¿(д-1)5г = а(д-1 )уг — Ь(д-1)щ и 7 определяется через свободные члены многочленов а(г) и ¿(г),

7 = Еу а° (С° )-1.

Минимальная матрица ковариаций ошибки оценивания ег = хг — хг равна

^е,ш[и а0 [¿0С0 ] а0^.

Доказательство. Искомую передаточную функцию фильтра от уг к оценке

хг обозначим 0(г). Передаточная функция от иг к хг подбирается так, чтобы она

сократилась в разности ег = хг — хг. Передаточные функции от (гшг,уг) к ошибке оценивания ег суть

= (I — С(г ))а(г)-1с(г),

Ше/у = —0(г),

где I — единичная матрица. Если объект наблюдения неустойчив и матрица а(г) имеет нули в единичном круге, то дополнительным ограничением на О является отсутствие полюсов у функции (I — О(г))а(г)-1 в замкнутом единичном круге. Ковариация ошибки оценивания есть

Т = Ееге* = У (Ше/ш Ше*/ш + Ше/УКу) ¿т(г) =

= !(Оа-1С — У) (Оа-1С — У)* — УУ* + а-1сКшс*(а-1)* ¿т(г),

У = а-1сЕт с* (сТ1)*, где Ст(г) обозначает нормированную меру Лебега на единичной окружности.

В соответствии с теорией Винера—Колмогорова разделим функцию У на две части: устойчивую У+ и полностью неустойчивую У—, причём У—(г) ^ 0 при г ^то:

У+(г) = а-1 С — Еу а(0)т (¿(0)т )-1,

УТ(г) = —Еу(а*(г) — а(0)т(¿(0)т)-1 С*(г)) С*(г)-1.

Оптимальная функция фильтра имеет вид

О(г) = У+ (г)С(г)-1а(г) = I — Еу а(0)т (¿(0)т )-1 ¿(г )-1а(г).

Оптимальный фильтр определяется уравнением

Хг = Уг — чС(д-1)-1(а(д-1)уг — Ь(д-1)щ), 7 = Еу а(0)т (С(0)т )-1.

При этом ошибка оценивания не зависит от управлений:

ег = Хг — хг = Уг — ^¿(д-1)-1 с^-1^.

Минимальная матрица ковариаций ошибки оценивания равна

УУ* —УУ*+а-1сЕшс* (а-1)*

Ст(г) = Еу—Еу а(0)т (¿(0)т) 1 ¿(0) 1а(0)Еу.

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

Лемма 2. Для каждой симметричной матрицы Еы определим устойчивый матричный многочлен ¿(г) равенством

С(г)С* (г) = а(г)Еу а* (г) + с(г)Еыс* (г),

если это выражение неотрицательно на единичной окружности.

Определим последовательности (6г), (уг) и (рг) уравнениями

с(д-1)6г = £г, С(д-1)Мг = £г, ^г = Мг — ¿-1со5г, £г = а(д-1 )уг — Ь(д-1)щ,

где Со и со — свободные члены матричных многочленов ¿(г) и с(г) соответственно.

Тогда среднее значение матричной квадратичной формы Ег = ]1гь'т' равно

Е^г = ¡Цг)^ — Еы)Ь*(г) Ст(г),

Ь(г) = ¿(г)-1 с(г) — ¿-1с о.

Доказательство. Из уравнения объекта следует, что

£г = а(д-1)уг + с(д-1)мг,

и это есть случайный стационарный процесс со спектральной плотностью

ЯЕ(г) = а(г)Еу а* (г) + с(г )Еы с* (г).

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

Е^г = J ¿(г)-1 (а(г)Еуа* (г) +с(г)Еыс* (г))(е* (г)-1 — с* (г)-1 ¿Г(К)-1) ¿т(г) =

= J ¿(г)-1(¿(г)¿* (г) + с(г)(Еы — Еы)с* (г))(3.'* (г)-1 — с* (г)-1с°(3)-1) ¿т(г) =

= 1 — J ¿* (г)с* (г)-1с°(3)-1 ¿т(г) — J ¿(г)-1с(г)(Еы — Еы)с°(3)-1) ¿т(г) +

+1- Еы)с*()Я*(гГ' ¿т(г).

По теореме о вычетах

J ¿* (г)с* (г)-1ст(3)-1 ¿т(г) = J ¿(г)т(с(г)-1)т(3)-1 ¿т(г) = 1,

j ¿(г)-1с(г)(Еы — Еы)ст(3)-1) ¿т(г) = 31 со(Еш — Еы)ст(3)-1) =

= I ¿-1со(Еы — Еы)с*(г)3*(г)-1) ¿т(г).

Отсюда

Е^г = / (¿(г)-1с(г) — 5-^1сЛ (Ет — Еы)(¿(г)-1с(г) — с-1^ ¿т(г).

4. Сходимость алгоритма оценивания. Пусть /г(Е) —некоторая функция от измеряемых величин в момент £ и от матрицы Е, которая заменяет в расчётах матрицу Еы. Пусть для каждого фиксированного значения Е матрицы /г(Е) в среднем сходятся к некоторой матрице ](К):

1т

-5>(Д)^ЕЛ(Д) = /(Д).

г=1

Рассмотрим рекуррентный алгоритм типа стохастической аппроксимации:

Д(+1 — Д* + £1/2+е/‘№)>

где к > 0, £ € (0,1/2). Функция f определяет следующее ассоциированное дифференциальное уравнение к этому алгоритму:

Ё = / (К).

Если это уравнение имеет стационарную точку К0, которая асимптотически устойчива, и если последовательность (Кг) попадает в инвариантную относительно этого уравнения область притяжения этой точки бесконечное количество раз почти наверное, то по теореме Л. Льюнга [7] последовательность (Кг) сходится к К0 с вероятностью 1.

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

R — р (Rw — R), р

d(z) 1c(z) — d- 1gq

2

dm(z) > 0,

очевидно, устойчиво, и все решения сходятся к значению К0 = Кш. Поэтому рекуррентный алгоритм сходится к истинному значению ковариации возмущения Кш, а фильтр сходится к оптимальному. Таким образом, в скалярном случае фильтр Калмана—Бьюси с дополнительной подстройкой неизвестного параметра Кш является адаптивным и предельно оптимальным [6].

В общем матричном случае рассмотрим отображение С множества всех симметричных матриц размера N х N в себя:

С : X ^ Y — J L(z)XL*(z) dm(z).

Если ядро этого оператора тривиальное, то С — биекция, и существует отображение С-1.

Теорема 1. Пусть R — некоторое выпуклое замкнутое множество симметричных матриц, содержащее Rw внутри себя. Пусть для любой матрицы R &R выполнены перечисленные ниже свойства.

1. Для любого числа z на единичной окружности матрица

D (z) — a(z) Rv a* (z) + c(z) Rc* (z)

положительно определена.

2. Пусть d(z) —результат факторизации D(z),

D(z) — d(z)d*(z),

причём матрица d(z) невырождена при \z\ < 1. Определим квадратную матричную функцию

L(z) — d(z )-1c(z) — d-1 co.

Тогда уравнение

CX — j L(z )XL(z )* dm(z)—0

имеет только нулевое решение в классе симметричных матриц X. При этом условии линейное отображение С обратимо.

За,висимость введённых величин от R будем обо,зна,чать С — Cr и d — d,R.

Rt+i = P

где к > 0 —параметр алгоритма, £ € (0, 1/2), последовательности (рг) и (иг) определяются уравнениями

с(ч-1')5ь = £г,

¿Яь (ч-1)^г = £г,

Vг = Цг — (Лпь )0 со5г^

где (¿Яь )о и со — свободные члены матричных многочленов ¿я (г) и с(г) соответственно, и

£г = а,(д-1)уг - Ь(д-1)щ.

Отображение V есть ортогональный проектор на множество матриц ^ в следующей метрике:

N

IRII2 = £ \Rij \2, R = (Rij)Nj=i ■

i,j=l

Тогда последовательность матриц Rt сходится к матрице Rw с вероятностью 1.

Доказательство. Для каждой матрицы R определим

F = % = С-1 Ft,

где последовательности (pt) и (vt) определены, как в лемме 2. Тогда по лемме 2 и из определения оператора L следует, что при фиксированной матрице R

Eft = L-l(EFt) = Rw - R■

Поэтому ассоциированное дифференциальное уравнение имеет вид

R = Rw — R,

и оно асимптотически устойчиво. Функцией Ляпунова является

V = IIR — Rw II2.

Ортогональная проекция P на множество R уменьшает эту функцию.

Выполнены все условия теоремы Л. Льюнга о сходимости рекурсивных алгоритмов [7]. По этой теореме последовательность Rt сходится к Rw почти наверное. □

Следствие 1. Пусть выполнены условия теоремы 1. Пусть (Rt) — последовательность матриц, определённая в этой теореме и dRt — соответствующие устойчивые матричные многочлены. Тогда фильтр

xt = yt — Yt^t,

является предельно оптимальным в смысле минимума дисперсии ошибки оценивания, где последовательность St определяется уравнением

dRt (q-í)St = a(q-1)yt — b(q-1)ut, а Yt —через свободные члены многочленов a(z) и d,Rt (z),

Jt = Rva0 ((dRt )0 ) ■

Вычисление значения оператора L-1 можно выполнить, например, следующим образом. Найдём такие матрицы p, q и D, что

L(z) = d(z)-1c(z) — d—1co = zp(I — Dz)-1q,

причём все собственные числа D лежат внутри единичного круга. Пусть

Y = LX = í L(z)XL*(z) dm(z)■

J\z\ = 1

Требуется найти X при известном Y. Рассмотрим уравнение

Q = í (I — Dz)-1qXqT (I — DT z-1)-1 dm(z )■

J\z\ = 1

Нетрудно проверить, что

Q — DQDT = qXqT, pQpT = Y.

В этой системе количество линейных уравнений совпадает с количеством неизвестных в матрицах Q и X. Отсюда они определяются.

В частном случае, когда уравнение объекта записано в пространстве состояний, степени многочленов a и d равны 1. Тогда матрицы p и q квадратные. Их невырожденность равносильна биективности оператора L, и матрица X легко находится через Y .

5. Пример. В качестве примера рассмотрим систему

Xt+1 = AXt + wt,

yt = Xt + vt,

где

. = (0.95 0.2\

A = ^—0.3 0.8y ,

wt и vt — белые шумы, с матрицами ковариаций

R = (и о\ R = (l 0\

Rw ^0 А) , Rv = ^0 l) ■

Данные моделировались с помощью датчиков случайных чисел в системе Matlab. На рис. 1 приведены данные по подстройке матрицы ковариации возмущений. По

Рис. 1.

7000 7200 7400 7600 7800 8000

Рис. 2.

7 6 5 4 3 2 1 О

7000 7200 7400 7600 7800 8000

Рис. 3.

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

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

Литература

1. Bucy R. S., Kalman R. E. New results in linear filtering and prediction theory // J. Basic. Eng. ASME. 1961. Vol. 83. N 1. P. 95-108.

2. Фомин В. Н. Рекуррентное оценивание и адаптивная фильтрация. М., 1977.

3. Mehra R. K. On the identification of variances and adaptive Kalman filtering. IEEE Trans. on AC. 1970. Vol. AC-15. P. 175-184.

4. Isaksson A. Identification of time varying systems through adaptive Kalman filtering. Proc. of the 10th IFAC Congress. Munich. 1987. Vol. 10. P. 305-310.

5. Ljung L., Soderstrom T. Theory and practice of recursive identification. Cambridge: MIT Press. 1983.

6. Барабанов А. Е., Лукомский Ю. А., Мирошников А. Н. Адаптивная фильтрация при неизвестной интенсивности возмущений и шумов измерений // Автоматика и телемеханика. 1992. №11. С. 93-101.

7. Ljung L. Analysis of recursive stochastic algorithms. IEEE Trans. on AC. 1977. Vol. AC-22. P. 551-571.

Статья поступила в редакцию 11 марта 2011 г.