Конспект по фильтрации

Мальцева Н. В.

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

51

2.Перейд¼м к последовательной обработке измерений, при которой любая группа из мерений yk(i) обрабатывается отдельно. Коэффициент усиления, оценку и ковариаци онную матрицу, полученные на основании данных измерений yk(i), будем обозначать индексом i. Тогда для i = 1 можно записать:

P (1)

(k|k)

−1

= (P (k|k − 1))−1 +

Hk(1)

T

Rk(1)

−1

Hk(1)

(15.3)

= P (1)(k|k)

Hk(1)

T

(15.4)

Kk(1)

Rk(1)

−

1

x(k|k − 1) ,

x(1)(k|k) = x(k|k − 1) + Kk

yk

− Hk

(15.5)

(1)

(1)

(1)

где значения x(k|k−1) и P (k|k−1) получаются с помощью основных формул фильтра

Калмана. Теперь, используя x(1)		è	P	(1)	обработаем измерение yk(2)		. Т. к. измерение
	(k|k)			(k|k)
yk(2)	относится к тому же моменту времени tk, ÷òî è yk(1)					, то экстраполяции не требу

ется. Тогда, аналогично предыдущим соотношениям, можно записать

P (2)(k|k) −

1

= P (1)(k|k − 1) −T

+ Hk(2)1

T

Rk(2)

−

1

1

K(2)

= P (2)(k

k)

H

(2)

R(2)

−

(2)

k

|

k (2)

(2)

k

x(2)(k|k) = x1(k|k − 1) + Kk

yk

− Hk

x(1)(k|k

Hk(1) (15.6)

(15.7)

− 1) , (15.8)

Подставляя в выражение для		P (2)(k k)		−1
получим		|			(15.6) выражение для
P (2)(k|k) −1	= P	(1)(k|k − 1) −1				2
						+ i=1 Hk(i) T Rk(i)
						X

P (1)(k|k) −1 (15.3),

−1 Hk(i) .

После обработки q систем данных измерений имеем:

P (q)(k|k) −

1

=

P (1)(k|k − 1) −

1

q

Hk

T

Rk

−

1

Hk

(15.9)

+ i=1

def

X

(i)

(i)

(i)

Сопоставляя (15.9) и (15.2), убеждаемся в их идентичности. Осталось получить выражение для x(k|k).

Подставляя (15.5) в (15.8), имеем:

x

(k|k) = x(k|k − 1) +

I − Kk

Hk

Kk

Kk

×

"y(2)

−H

(2)x(1)

(k|

k−

1)#

(2)

h

(2)

(2)

(1)

(2)

i

yk(1)

Hk(1)x(k k

1)

k

−

k

|

−

Сопоставляя аналогичным образом все q систем измерений, получим:

(q) | def

x (k k) =

(1)

(1)

x(k|k) = x(k|k − 1) + Kk ·

yk

− Hk

...x(k|k

− 1)

,

(15.10)

e

yk(2)

−

Hk(2)x(1)(k

k

−

1)

|

Мальцева Н. В.

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

52

ãäå

"

Kek = I − Kk(q)Hk(q) I − Kk(q−1)Hk(q−1) · . . . · I − Kk(2)Hk(2) Kk(1) . . .

#

. . . I − Kk(q)Hk(q) Kk(q−1) Kk(q)

Сравнивая (15.1) и (10.10) видим, что надо доказать равенство коэффициентов уси

m

ления Kk è Kek. Покажем, что они имеют одинаковые блоки. n Kk

Начн¼м с последнего блока .

В коэффициенте Kek:

îäíà èç

(q) ôîðì Kkq

Kk =

ïî

условию

теоремы

Pdef

(k|k)

Hk

T

Rk

−

1

= P (k|k)

Hk

T

Rk

−

1

(q)

(q)

(q)

(q)

(q)

|

= P{z(k|k)}

áëîê â Kk

|

последний

{z

}

= последние блоки совпадают.

Рассмотрим предпоследний блок .

|

условию

|

−

|

P (k

k) = P

(q)(k

k) =

I

K

(q)H(q)

P (q−1)

(k

k)

(15.11)

ïî

k

k

теоремы

Домножим это равенство справа на Hk(q−1) T Rk(q−1) −1

:

P (k|k) Hk(q−1) T Rk(q−1) −1 = I − Kk(q)Hk(q) Kk(q−1)

|

{z

K

}

|

{z

Kk

}

предпоследний блок в

предпоследний блок в

e

Далее, используя P (q−1)(k|k) =

I − Kk(q−1)Hk(q−1)

P (q−2)(k|k) и, подставляя это ра

венство в (15.11), имеем:

P (k|k) = I − Kk(q)Hk(q) I − Kk(q−1)Hk(q−1) P (q−2)(k|k)

Умножив последнее равенство справа на

H

(q−2)

T

R(q−2)

−1

, получим равенство

третьих блоков от конца матриц

Kk

è

k

k

Kk. Повторяя заданную процедуру для номеров

i = q − 3, . . . , 1, находим, что Kk = Kk,e÷. ò. ä.

последовательная обработка измерений да¼т одина

Таким образом, одновременная и

e

ковые результаты.

Мальцева Н. В.

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

53

16Методы фильтрации на основе использования квадратных корней из матриц

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

при потере матрицей положительной определ¼нности;

при плохой обусловленности обращаемых матриц.

Потеря положительной определ¼нности матрицами P (k|k) или P (k|k − 1) при исполь зовании ЦВМ может произойти в двух случаях:

при формировании апостериорной корреляционной матрицы с помощью соотношения P (k|k) = (I −KkHk)P (k|k −1) при условии близости к 1 какого-либо из диагональных элементов матрицы KkHk.

при условии округления матрицы Φ(k+1|k)P (k|k)ΦT(k+1|k)+ (k+1|k)Qk T(k+1|k) до матрицы Φ(k + 1|k)P (k|k)ΦT(k + 1|k) и близости к 0 какого-либо из диагональных элементов этой матрицы.

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

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

Мы рассмотрим несколько вариантов метода извлечения корней квадратных из матрицы

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

Основная èäåÿ метода фильтрации с использованием корней квадратных из матриц состоит в следующем:

1)экстраполированную и апостериорную корреляционную матрицы оцениваемых па раметров представляют в виде произведения некоторых невырожденных матриц на транспонированные по отношению к ним матрицы;

2)соотношения для расч¼та корреляционных матриц и коэффициента усиления заменя ют эквивалентными соотношениями, позволяющими:

неявно гарантировать положительную определ¼нность корреляционных матриц;

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

Впервые эта идея была реализована Поттером при отсутствии ошибок модели систе

мы для скалярных измерений . Корреляционная матрица ошибок оценивания представля лась в виде P = SST, где S некоторая невырожденная матрица.

Итак, рассмотрим подробнее алгоритм Поттера.

Алгоритм (Поттера). Äàíî:

Мальцева Н. В. Рекуррентная и адаптивная фильтрация								54
1. Модель системы				xk = Φ(k|k − 1)xk−1
2. Модель измерений
	yk = Hkxk + vk ;					yk скаляр R1 .
3. Априорные данные
vk N(0, Rk) ;	cov(vk, vj) = Rkδkj ;					x0 N(x0, P0) ;	cov(vk, x0) = 0 ;
Доказать:
1. Алгоритм экстраполяции значений
x(k|k − 1) = E(xk|Y1k−1)					è	P (k|k − 1) = cov(xk, xk|Y1k−1) :
	x(k|k − 1) = Φ(k|k − 1)x(k − 1|k − 1)							(16.1)
	S(k|k − 1) = Φ(k|k − 1)S(k − 1|k − 1)							(16.2)
	P (k|k − 1) = S(k|k − 1)ST(k|k − 1)							(16.3)
2. Алгоритм фильтрации значений
x(k|k) = E(xk|Y1k)					è	P (k|k) = cov(xk, xk|Y1k) :
x(k|k) = x(k|k − 1) + Kk					yk − Hkx(k|k − 1)			(16.4)
Kk = αkS(k|k − 1) · Fk								(16.5)
S(k k) = S(k k − 1)			I − αkγkFkFkT					(16.6)
P (k|k) = S(k|k)ST(k k) ,								(16.7)
|	|		|
ãäå
Fk = ST(k|k − 1)HkT								(16.8)
αk = (FkTFk + Rk)−1						скаляр		(16.9)
γk = (1 + p				)−1 ;		γk > 0
		Rkαk						(16.10)

Доказательство. Используя соотношения классических формул фильтра Калмана приме нительно к рассматриваемым моделям, имеем:

x(k|k − 1)		= Φ(k|k − 1)x(k − 1|k − 1)								(16.11)
P (k|k − 1)		= Φ(k|k − 1)P (k − 1|k − 1)ΦT(k − 1|k − 1)								(16.12)
x(k|k)		= x(k|k − 1) + Kk(yk − Hkx(k|k − 1))							)−1	(16.13)
K	k	= P (k k	1)HT(H	P (k	k	−	1)HT + R			(16.14)
		| −	k k	|			k	k
P (k|k)		= (I − KkHk)P (k|k − 1)								(16.15)

Отсюда следует (16.2).

Мальцева Н. В. Рекуррентная и адаптивная фильтрация 55

Согласно теореме о существовании квадратного корня из симметричной положительно определ¼нной матрицы, для корреляционных матриц P (k|k) и P (k|k − 1) соответственно

существуют неособенные матрицы S(k|k) и S(k|k − 1) такие, что возможно представление

P (k|k)	= S(k|k) · ST(k|k)	(16.16)
P (k|k − 1)	= S(k|k − 1) · ST(k|k − 1)	(16.17)

Подставим (16.17) в (16.12) (вид прогнозируемой матрицы):

h ih iT

P (k|k − 1) = S(k|k − 1) · ST(k|k − 1) = Φ(k|k − 1)S(k − 1|k − 1) Φ(k|k − 1)S(k − 1|k − 1)

(S(k|k − 1) = Φ(k|k − 1)S(k − 1|k − 1))

Подставим (16.17) в (16.14) (вид Kk):

Kk = S(k|k − 1)ST(k|k − 1)HkThHkS(k|k − 1)ST(k|k − 1)HkT + Rki−	(16.18)
	1

Введ¼м обозначения:

1

Fk = ST(k|k − 1)HkT	1	n Fk
def
αk = HkS(k|k − 1)ST(k|k − 1)HkT + Rk −		= (FkTFk + Rk)−1 ,
def

где α число.

Тогда соотношение (16.18) принимает вид требуемый вид (16.5)

− 1)Fk)

вектор (16.19)

(16.20)

(Kk = αkS(k|k −

Наконец, подставим (16.17) и (16.5) в (16.15). Учитывая (16.19), имеем:

P (k|k) = S(k|k)ST(k|k) =

(I−KkHk)P (k|k−1)=P (k|k−1)−KkHkP (k|k−1)

Kk èç (16.5)

HkS(k|k−1)

(16.19)

FkT

=

z

}|

{ z

}|

{

P (

= S(k|k − 1)ST(k|k − 1) − αkS(k|k − 1)Fk Hk S(k|k − 1)ST(k|k − 1) =

= S(k k 1) |

I αkFk

{z

}

k

|

−

1)

k

|

−

−

F T

ST(k

k

−

1) (16.21)

k

|

Пусть γk неопредел¼нный множитель, который обеспечивает представление

I − αkFkFKT = I − αkγkFkFkT I − αkγkFkFkT T

(16.22)

В этом случае из (12.10) и (16.22) непосредственно следует (16.6). Оста¼тся только найти положительный неопредел¼нный множитель γk.

Раскроем соотношение (16.22):

I − αkFkFKT = I − 2αkγkFkFkT + αk2γk2 FkFkT FkTFk ,
èëè
αk FkTFk γk2 − 2γk + 1 = 0	(16.23)

Мальцева Н. В. Рекуррентная и адаптивная фильтрация							56
Равенство (16.23) квадратное уравнение относительно неизвестного множителя γk:
γk = 1 ± q						αk FkTFk
			1 − αk FkTFk				(16.24)
Из (16.20) следует:
FkTFk = 1/αk − Rk αkFkTFk = 1 − αkRk							(16.25)
Подставляя (16.25) в (16.24):
	γk = 1 ± p				. (1 − αkRk)
				αkRk			(16.26)
Т. к. по предположению γk > 0, то в (16.26) необходимо выбрать один знак
γk = 1 − p		. h(1 − αkRk)(1 + αkRk)i = (1 + αkRk)−1
	αkRk

и тем самым, равенство (16.10) доказано.

Преимущества метода Поттера :

гарантированное обеспечение неотрицательной определ¼нности матриц P (k|k − 1) и P (k|k) (вид матриц);

обусловленность матриц S(k|k − 1) и S(k|k) лучше, чем P (k|k − 1) и P (k|k)

Число обусловленности	def	λmax
матрицы A	= N(A) =		,
		λmin

ãäå λmax, λmin максимальное и минимальное собственное число матрицы A.

Тогда, можно показать, что для матрицы B:

r

A = BTB N(B) = λmax , ò. å. N(A) = [N(B)]2

λmin

Как это влияет?

При вычислении в десятичной системе исчисления с l значащими разрядами чис ленные трудности из-за округления результата возникают тогда, когда N(A) прибли жается к 10l. Фильтр же Поттера может работать со значениями N(B) = 10l èëè

N(A) = 102l.

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

Недостатки метода Поттера :

обработка только скалярных измерений;

модель системы не содержит шумов.

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

Мальцева Н. В. Рекуррентная и адаптивная фильтрация

57

Алгоритм (Эндрюса). Äàíî:

1.

Модель системы

xk = Φ(k|k − 1)xk−1 + wk−1

2.

Модель измерений

yk = Hkxk + vk

3.

Априорные данные

vk N(0, Rk) ;

wk N(0, Rk) ;

x0 N(x0, P0) ;

k 6= j

cov(vk, vj) = cov(vk, wj) = cov(vk, x0) = cov(wj, x0) = 0 ,

Доказать:

1.

Алгоритм экстраполяции значений

x(k|k − 1) = E(xk|Y1k−1)

è

P (k|k − 1) = cov(xk, xk|Y1k−1) :

x(k|k − 1) = Φ(k|k − 1)x(k − 1|k − 1)

(16.27)

ϕk = Φ(k|k − 1)S(k − 1|k − 1)

(16.28)

S(k|k − 1) =

ϕkϕkT + Qk

1/2

(16.29)

P (k k

1) =

k

k

1) ,

(16.30)

|

−

S(k

−

1)ST(k

−

|

|

ãäå S(k|k − 1) нижняя треугольная матрица. 2. Алгоритм фильтрации значений

x(k|k) = E(xk|Y1k)

è

P (k|k) = cov(xk, xk|Y1k) :

K = S(k k 1)F G

−

G−

x(k|k) = x(k|k −

1) + Kk

yk − Hkx(k|k

− 1)

(16.31)

S(k k) = S(k|k −

I

T

1

1

1)

F

GT k−1 (G

+ B )

1F T

(16.33)

k

h

k

k

k

k

k i

(16.32)

|

| −T

−

k

k −

(16.34)

P (k|k) = S(k|k)S

(k|k) ,

ãäå S(k|k) нижняя треугольная матрица

Fk = ST(k|k − 1)HkT

(16.35)

Gk = FkTFk + Rk

1/2

(16.36)

Bk

= k

R1/2 .

(16.37)

Доказательство. Соотношения фильтра Калмана в данном случае имеют вид (16.11), (16.13) (16.15) и

P (k|k − 1) = Φ(k|k − 1)P (k − 1|k − 1)ΦT(k|k − 1) + Qk

(16.38)

Мальцева Н. В.

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

58

Выразим корреляционные матрицы P (k|k) и P (k|k−1) (симметричные, положительно опре дел¼нные), соответственно, через неособенные нижнетреугольные матрицы 4 Sk è Sk−1:

P (k|k) = S(k|k)ST(k|k) ;

P (k|k − 1) = S(k|k − 1)ST(k|k − 1) ;

Тогда подставим (16.39) в (16.38):

S(k|k − 1)ST(k|k − 1) = Φ(k|k − 1)S(k − 1|k − 1) Φ(k|k − 1)S(k − 1|k − 1) T +

Обозначим

def

ϕk = Φ(k|k − 1)S(k − 1|k − 1) ,

и придадим равенству (16.40) иной вид:

(16.39)

Qk

(16.40)

(16.41)

S(k|k − 1)ST(k|k − 1) = ϕkϕkT + Qk = ϕkϕkT + Qk 1/2 ϕkϕkT + Qk 1/2 T

(16.42)

| {z }

симметричная

положительно

определ¼нная

матрица

Èç (16.42) = (16.29)

= справедливость алгоритма экстраполяции доказана.

Получим новое соотношение для коэффициента усиления. В соотношение (16.14) (Kk) подставим (16.39).

Fk

|

−

z

|

}|

{

= S(k k 1)F

GT

−1 G−1

, (16.43)

Kk = S(k k

1) ST(k k − 1)HkT

HkS(k|k − 1)ST(k|k − 1)HkT

+ Rk

−1

= подставим (16.35), (16.36)

|

−

k

k

k

а это есть соотношение (16.32) (вид Kk), что и хотелось получить.

Преобразуем соотношение для апостериорной корреляционной матрицы (16.15), подста вив в него (16.39) и (16.43):

P (k|k) = S(k|k − 1)ST(k|k − 1) = S(k|k − 1) hI − Fk GTk −1 G−k 1FkTiT ST(k|k − 1) . (16.44)

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

(Gk + Bk)(Gk + Bk)T = Gk(Gk + Bk)T + (Gk + Bk)GkT + (BkBkT − GkGkT) .

(16.45)

4

Определение. Квадратная матрица A размера n × n, удовлетворяющая условию R(A) = n, называется неособенной (невырожденной).

Мальцева Н. В. Рекуррентная и адаптивная фильтрация					59
ïî (16.37)	1/2			1/2
Åñëè Bk =	Rk , òî, ñ ó÷¼òîì (16.36)		Gk = FkTFk + Rk
	1/2	1/2		T
	BkBkT − GkGkT = Rk	Rk		− FkTFk − Rk = −FkTFk	(16.46)

Единичную матрицу с уч¼том (16.46) запишем в виде:

I= (Gk + Bk)−1(Gk + Bk)(Gk + Bk)T (Gk + Bk)T −1 =

=(Gk + Bk)−1 Gk(Gk + Bk)T + (Gk + Bk)GTk − FkTFk (Gk + Bk)T −1 = раскрыли скобки

=(Gk + Bk)−1Gk + GTk (Gk + Bk)T −1 − (Gk + Bk)−1FkTFk (Gk + Bk)T −1 . (16.47)

Преобразуем выражение в квадратных скобках в (16.44), учтя (16.47):

I− Fk GTk −1 G−k 1FkT = вставим I èç (16.47)

=I − Fk GTk −1 h(Gk + Bk)−1Gk + GTk (Gk + Bk)T −1 −

−(Gk + Bk)−1FkTFk (Gk + Bk)T −1 iG−k 1FkT =

=hI − Fk GTk −1 (Gk + Bk)−1FkTi hI − Fk GTk −1 (Gk + Bk)−1FkTiT . (16.48)

Подставив (16.48) в (16.44), получим требуемое соотношение (16.33), ч. т. д.

Недостатки метода Эндрюса

1.Потеря матрицей S(k|k) треугольного вида при использовании соотношения (16.33).

На любом шаге вычислений необходимо эту матрицу возводить в квадрат, а затем извлекать квадратный корень из матрицы S(k|k)ST(k|k). Такая процедура позволяет восстановить треугольный вид матрицы S(k|k).

2.Увеличение необходимого объ¼ма памяти ЦВМ и количества вычислений по сравне нию с калмановским фильтром.

Замечание. При доказательстве использовался тот факт, что если матрица A симметрич на и положительно определена, то существует нижнетреугольная не особенная матрица B такая, что A = BBT. Матрица B имеет следующие элементы.

b

=

aii1/2 ,

i

1

1/2

ïðè

i = 1

ii

aii

−

2

,

ïðè

i = 2, 3, . . . , n ;

−

bik

!

X

0 ,

i < j

ïðè

k=1

b

=

aijb−1

,

ïðè

j = 1 ; i = 2, 3, . . . , n

ij

jj

j

1

−

1

ïðè

aii −

bikbjk! bjj−

,

i > j > 1

X

k=1

Этот метод называется методом Чолеского.

Мальцева Н. В. Рекуррентная и адаптивная фильтрация 60

17 Интерполяция компонент гауссовских процессов

Интерполяция случайного процесса xt состоит в построении оценки x(t|s) первой ком поненты xt гауссовского процесса {xt, yt} по наблюдениям Y0s = {yi, i [0, s]}, s > t è

представляет собой более общую задачу оценивания, чем задача фильтрации, при которой строится оценка x(t|t) по наблюдениям Y0t.

Полезность интерполяции состоит в том, что она позволяет улучшить оценку x(t|t),

являющуюся результатом решения задачи фильтрации по наблюдениям Y0t, çà ñ÷¼ò ïðè влечения избыточной измерительной информации Yts, полученной после момента времени

t, на который строится оценка процесса xt.

Рассмотрим решение тр¼х классов задач дискретной интерполяции:

1.Интерполяция в закрепл¼нной точке, или интерполяция в прямом времени, при кото рой строится оценка x(N|k), k = N + 1, N + 2, . . . на фиксированный момент времени tn по нарастающей выборке измерений Y1k.

2.Интерполяция на закрепл¼нном интервале , или интерполяция в обратном времени, при которой определяется оценка x(k|N), k = N − 1, N − 2, . . . , 1, 0, на любой из моментов времени tk [t0, tN ] после получения выборки измерений Y0N , относящейся ко всему интервалу времени [t0, tN ]. Отправной точкой для таких оценок является x(N|N).

3.Интерполяция с постоянным запаздыванием , при которой формируется оценка x(k|k+

+ N) на последовательно возрастающие моменты времени tk, k = 1, 2, . . . , запазды

вающие на постоянную величину (tk+N − tk) относительно времени tk+N последнего измерения, принадлежащего выборке Y0k+N .

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

Теорема (о нормальной корреляции). Пусть

x N(mx, Pxx) ,				x Rn
y N(my, Pyy) ,		Pyy		y Rm
{x, y} N my	, Pyx		,	(x, y) Rn+m
mx	Pxx	Pxy
Тогда условное математическое ожидание			mx|y = E(x|y) и условная корреляционная

матрица Pxx|y = cov(x, x|y) удовлетворяют соотношениям:

mx|y = mx + K(y − my) ;

K = Pxy · Pyy−1 ;

Pxx|y = Pxx − KPxyT

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

Следствие. 1. Пусть {(x1, . . . , xn)T, (y1, . . . , ym)T, (z1, . . . , zl)T} условно-гауссовский вектор, т. е. π(x, y|z) гауссовская условная плотность распределения вероятно

ñòåé.