Численно устойчивые реализации фильтра Калмана для оценивания линейных парных марковских моделей с гауссовым шумом Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 22, № 3, 2017

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

М.В. Куликова1, Ю.В. Цыганова2'*

1 Высший технический институт Лиссабонского университета, Португалия 2Ульяновский государственный университет, Россия *Контактный e-mail: TsyganovaJV@gmail.com

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

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

Ключевые слова: парная марковская модель, линейные стохастические системы, парный фильтр Калмана, устойчивость к ошибкам округления, U^-реализации.

Введение

Одной из важнейших задач теории статистических решений, имеющих большое практическое значение, является задача оптимальной фильтрации или оценивания неизвестного процесса Х^ = {xi}f=1 на основе доступного для наблюдений процесса Y-^ = {у}=1. Теория оптимальной нелинейной фильтрации для скрытой марковской модели (СММ) была развита в работах Стратоновича, она базируется на теории условных марковских процессов [1]. В СММ предполагается, что неизвестный процесс Х^ = {xi}f=1 марковский. В случае линейной дискретной СММ, возмущаемой гауссовым белым шумом, решением задачи оптимальной фильтрации является фильтр Калмана [2], а также фильтр Калмана — Бьюси, разработанный для непрерывных моделей [3].

Парные марковские модели (ПММ), изучаемые в данной работе, подразумевают, что пара (Х^,Y-^) — это марковский процесс. Из этого следует их основное отличие от СММ: в ПММ неизвестный процесс Х^ = значение которого необходимо оценить на основе доступного для наблюдений процесса Y-^ = {yi}j=1, не обязательно марковский; доказательство, примеры и детальное обсуждение можно найти

© ИВТ СО РАН, 2017

в [4, гл. 14], а также в [5, с. 165]. Класс ПММ, а также задача оптимальной фильтрации, поставленная для моделей данного типа, впервые описывались именно в отечественной литературе [4, 6]. Этот факт отмечается и в более поздних зарубежных источниках [7, разд. 3]. В отечественной литературе с конца 1960 г. изучению задачи фильтрации в классе ПММ или похожего типа моделей посвящены работы [8- 10] и др. За рубежом особый интерес к ПММ возник в последние годы, и связан он со свойствами моделей данного класса. В частности, показано, что использование ПММ вместо СММ в задачах обработки изображений позволяет сократить ошибку в два раза [11]. Таким образом, особое распространение класс ПММ получил в задачах обработки изображений [12, 13]. Общие свойства ПММ и других подобных классов моделей, а также их применение для решения прикладных задач рассмотрены в [14, 15].

Для линейной гауссовой ПММ можно сформулировать метод оценивания неизвестного процесса аналогично тому, как это сделано для СММ. Полученная техника оценивания получила название парного фильтра Калмана (ПКФ) [7, 12]. В ряде работ можно найти подтверждение того, что использование ПКФ вместо классического фильтра Калмана позволяет более точно решать те или иные практические задачи. В частности, в [16] решается задача слежения за движением зрачка глаза человека, где применение ПКФ позволяет получить более точные оценки, чем при использовании классического фильтра Калмана. С точки зрения алгоритмической реализации фильтров следует отметить, что парный фильтр (так же, как и классический фильтр Калмана [17]) чувствителен к ошибкам машинного округления. Чтобы сократить влияние ошибок округления на точность оценивания, ранее был развит устойчивый квадратно-корневой ПКФ [18, 19]. В указанных работах также предложен адаптивный вариант парного фильтра для одновременного оценивания неизвестных параметров модели и неизвестного вектора состояния динамической системы, построенный на основе ЕМ-алгоритма. Более подробная информация об оценивании неизвестных системных параметров методом максимума правдоподобия с помощью ЕМ-алгоритма приведена в работе [20]. Адаптивный ПКФ на основе градиентных методов оптимизации и квадратно-корневой формы фильтра предложен в [21].

Целью данной работы является изучение устойчивых реализаций парного фильтра. Рассмотрены особенности применения квадратно-корневых методов для реализации ПКФ, а также построен новый [/Д-алгоритм парной фильтрации. Для классического фильтра Калмана показано, что квадратно-корневые алгоритмы позволяют обрабатывать данные с двойной точностью1 [22, с. 729], однако их вычислительная сложность чуть выше, чем у стандартной реализации фильтра. Из класса квадратно-корневых методов выделяется класс [/Д-фильтров, который является их модификацией [23]. Он предполагает разложение матрицы ковариации на нижнюю треугольную, диагональную и верхнюю треугольную матрицы. Класс иД-реализаций обеспечивает такую же точность вычислений, как и все квадратно-корневые методы, но при этом позволяет сократить затраты машинного времени и избавляет от необходимости использовать опе-

хЗдесь имеется в виду, что число обусловленности положительно определенной матрицы ^(А) =

) = [^(й*)]2, где А = Б8Т и Б — нижняя треугольная матрица с положительными элементами на диагонали, т.е. Б — квадратный корень матрицы А, полученный в результате разложения Холецкого. Это означает, что если вычислительный алгоритм сформулирован в терминах матрицы А и влияние ошибок округления становится значительным при р(А) = 10р, то алгоритм, сформулированный в терминах матричного квадратного корня будет подвержен ощутимому влиянию ошибок округления только при р(А) = 102р, что эквивалентно вычислениям, проводимым с двойной точностью.

рацию извлечения квадратного корня. Несмотря на все их преимущества, [/Д-методы для ПКФ ранее не были построены. Здесь эта задача решена впервые.

Таким образом, в настоящей работе описан разработанный авторами [/Д-алгоритм для парного фильтра, а также показано несомненное преимущество квадратно-корневых и [/^-реализаций ПКФ по сравнению с его стандартной версией, которая работает с полными матрицами ковариации ошибки оценивания на каждом шаге алгоритма фильтрации.

1. Постановка задачи

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

хк+1 1 х,х 1 х,у хк

Уй Ру,Х ^у,у Ук-1

/

^к+1 Т

+

Щ

у

WI

™ к

Щ

У

Qx,

х,у

Qx,y Qy,y

а

Здесь € ), пх > 0, пу > 0, — составной вектор, в котором хк+1 € Е""* —

неизвестный вектор стохастической системы (1) и ук € — доступный вектор наблюдений; } — последовательность нормально распределенных случайных векторов размера пх + пу с нулевыми математическими ожиданиями и ковариационной матрицей О, € )х(пх+пу), где О, > 0; {щд.} не зависит от начального вектора состояния 1с = хо, хо - М (Хо, Ро), Ро € , где Р0 > 0 [7, 12].

Поставлена задача — оценить неизвестный процесс {х^}^=о на основе доступного для наблюдений процесса {уг}^=о. Парная марковская модель (1) — это более общий случай скрытой марковской модели — классической линейной дискретной стохастической системы в пространстве состояний, для оценивания неизвестного вектора состояния которой разработан фильтр Калмана [2]. Однако для парной линейной модели (1) также возможно сформулировать аналогичный метод оценивания (см. доказательство и вывод формул в [7, утверждение 1]). Такая техника получила название парного фильтра Калмана. Она позволяет находить среднюю квадратическую оценку х^ неизвестного вектора Хд. линейной парной марковской системы (1), возмущаемой гауссовым шумом, по доступным результатам наблюдений {уо,... , у д.}. Приведенный ниже алгоритм 1 представляет собой стандартную реализацию ПКФ (см. уравнения (20)-(26) в [7]). Краткое изложение алгоритма содержится в [18, с. 5-6].

Парный фильтр (как и классический фильтр Калмана [17]) чувствителен к ошибкам машинного округления, т. е. его стандартная реализация (алгоритм 1) численно неустойчива по отношению к ошибкам округления. Квадратно-корневые реализации представляют собой класс методов, которые позволяют обрабатывать данные с двойной точностью, однако их вычислительная сложность чуть выше, чем у стандартной реализации фильтра. Для классического фильтра Калмана первые квадратно-корневые алгоритмы были разработаны в [22, 24, 25]. Для парного фильтра Калмана, изучаемого в данной работе, квадратно-корневая реализация впервые предложена в [18], где авторы использовали идеи работы [22], развитые ранее для классического фильтра Калмана.

Алгоритм 1: Парный фильтр Калмана: стандартная реализация Входные данные : Начальные данные фильтра хо, Ро > 0 и результаты измерений

{yo,..., yn}. Предполагается, что Qy,y > 0. Выходные данные: Средние квадратические оценки {£j\j}, i = 1,..., N, по

доступным результатам наблюдений {y0,..., yn }. Инициализация фильтра: Присвоить хо\о = xo, Ро\о = Ро и определить

РХ'Х - РХ'Х Qx,y Q-У'У Ру'Х , Рх'У - ^Х'У Qx,y Qy,y Ру'У,

QX'X - QX'X Qx,yQy,yQy,X.

for k = 0 to N — 1 do

Этап экстраполяции. Используя полученную в момент времени к оценку Хк\к, предсказать значение xfc+1|fc для следующего момента времени к + 1:

оценка вектора состояния xfc+i\fc = Fx>xZk\k + Qx>yQ-- yk + Fx>yyk-i, ковариация ошибки предсказания Pk+l\k = Qx,x + Fx^xPk\kPl,x. (2)

Этап обработки измерений и фильтрации. Используя поступившие в момент времени к + 1 данные yk+i, обновить оценку вектора состояния:

невязка измерений фильтра ek+- = yk+i — Fy 'X xfc+l\fc Py,yyk,

ковариация невязки измерений Re,k+l = Qy,y + Py,xPk+l\kPy,x, (3)

коэффициент обратной связи Kk+l = Pk+l\kF^^iXRyyl.+l, (4)

ковариация ошибки оценивания Pk+l\k+l = Pk+l\k — Kk+lRe,k+lK]:+l, (5)

оценка вектора состояния xfc+i\fc+i = xfc+i\fc + Кк+-ек+-. (6)

end

В следующем разделе подробно рассмотрен класс квадратно-корневых методов для парного фильтра Калмана. Здесь А1/2 — верхняя треугольная матрица с положительными элементами на диагонали, полученная при разложении Холецкого вида А = Ат/2А1/2. Для удобства введены обозначения: АТ/2 = (А1/2)т, А-1/2 = (А1/2)-1, А-Т/2 = (А-1/2)т.

2. Класс квадратно-корневых алгоритмов парной фильтрации

Рассмотрим базовые принципы построения квадратно-корневых алгоритмов для классического фильтра Калмана. Они основаны на обновлении разностного уравнения Рик-кати в терминах треугольных матриц, полученных из матрицы ковариации Р^ с помощью разложения Холецкого. Отметим, что разложение Холецкого осуществляется только один раз — для начальных значений2, т. е. Р0 = Р$/2Р12, и далее обновление уравнений фильтра происходит в терминах треугольных матриц Р^/2. На последнем этапе

фильтрации, т. е. после обработки последнего измерения, вновь имеем Р^= Р^/мРЦ\м. Такая организация вычислений не свободна от влияния ошибок округления, однако она гарантирует симметричность ковариационных матриц и положительные диагональные элементы.

2 Для начальных условий разложение Холецкого матрицы Ро существует, так как Ро — вещественная симметрическая матрица с положительными ведущими минорами, поскольку Ро > 0, т. е. строго положительно определенная [26].

Алгоритм 2: Квадратно-корневой парный фильтр

Входные данные : Начальные данные фильтра хо, Ро > 0 и результаты измерений

(уо,..., ум}• Предполагается, что Qy,y > 0. Выходные данные: Средние квадратические оценки (х^}, г = 1,..., N, по

доступным результатам наблюдений (у0,..., ум }• Инициализация фильтра: Выполнить разложение Холецкого для Ро = РеТ^Ро^2 •

Присвоить начальные значения хо|о = хо, Рс

1/2 0|0

Р

1/2

Рх,х - Рх,х Qx,y^^ у,уFy,X, Рх,у - Рх,у Qx,y^^ у,уРг

-1

у,у

(ё^х,х - Qx,x Qx,y^у,уQy,X.

Выполнить разложение Холецкого для матриц ^ ^Т/2^1/2 ^ р.Т/2 р.1/2

ковариации Цуу — Чу,у Чу,у и Цх,х — Чх,х Чх,х •

Гсг к = 0 1с N - 1 ас

Этап экстраполяции. Используя полученную в момент времени к оценку х^, предсказать значение Хк+^к для следующего момента времени к + 1:

хк+1|к = Рх,ххк1к + Ях,у Я-^у к + Рх,уу k-1, 1/2

оценка вектора состояния квадратный корень матрицы

Р,

к+1|к 0

оТи

р1/2 рТ ■

к| к х, х

Vх,х

(7)

где — ортогональное преобразование, приводящее к верхнему треугольному виду матрицу, стоящую в левой части (7).

Этап обработки измерений и фильтрации. Используя поступившие в момент времени к + 1 данные ук+ъ обновить оценку вектора состояния:

найти значения

К

1/2 е,к+1

0

*I

1/2

к+1|к+1_

'/,к+1 р 1/2

2

ми к+1

о1/2

0

р 1/2 рТ р 1/2 _ к+1|к^у,х ^к+1|к_

(8

где 0,1М++и1 — матрица ортогонального преобразования к верхнему треугольному виду левой части (8); К/,к+1 = Рк+1|кР-у,х^—к/2 — нормализованный коэффициент

■Т/2,

обратной связи фильтра, т.е. К/,к+1 = Кк+1Ке'к

невязка измерений оценка состояния

хх

е к+1 = у к+1 — Ру,ххк+1|к — Ру,у у к,

— ~ I ^ Т)—Т /2

к+1 |к+1 = хк+1|к + Л/,к+1Ле,к+1ек+1.

(9)

епа

Современные квадратно-корневые методы, построенные для классического фильтра Калмана, используют ^Р-разложение для обновления квадратных корней матриц ковариации ошибки оценивания и предсказания. Каждая итерация таких реализаций фильтра имеет следующий вид: формируется блочная матрица Л, которая содержит всю известную на данном этапе информацию об элементах фильтра. Затем ортогональное преобразование Q применяется к Л таким образом, что матрица В, полученная в результате ортогонального преобразования, имеет специальную треугольную форму (блочная верхняя или нижняя треугольная). Имеем 0,Л = В. Ортогональное преобразование может быть осуществлено различными методами. Полученная блочная треугольная матрица В содержит обновленные величины фильтра, которые напрямую

считываются из В [27, 28]. Подобные методы не только более устойчивы по отношению к ошибкам округления (как и все квадратно-корневые алгоритмы), но и удобны для использования на практике. Каждая итерация фильтра может быть представлена в виде одного блочно-матричного преобразования, а все необходимые величины для следующей итерации фильтра напрямую доступны из В в виде отдельных блоков. Такая компактная запись дает возможность применять параллельные вычисления при практической реализации квадратно-корневых фильтров [29].

Для парного фильтра Калмана первая подобная квадратно-корневая реализация разработана в [18], где использованы идеи работы [22], развитые ранее для классического фильтра Калмана. Алгоритм 2 представляет собой краткое изложение квадратно-корневой реализации [18, с. 6]. Каждая итерация фильтра имеет вид 0,Л = В, где В — блочная верхняя треугольная матрица.

Утверждение 1. Алгоритм 2 алгебраически эквивалентен алгоритму 1.

Доказательство приведено в Приложении.

Отметим некоторые свойства предложенного в [18] квадратно-корневого парного фильтра Калмана (алгоритм 2). Заметим, что стандартная реализация фильтра (алгоритм 1) для вычисления оценки вектора состояния требует обращения матрицы ковари-ации невязки измерений Ке,к+\ (см. уравнения (4), (6)). Квадратно-корневая реализация (алгоритм 2) алгебраически эквивалентна стандартной реализации фильтра, но вместо вычисления В-\+1 требует обращения ее квадратного корня Д1/2+1. Действительно, в

формуле (9) величины К/,к+1 и Р1/к+1 доступны в виде отдельных блоков из матрицы, полученной в результате ортогонального преобразования и стоящей в левой части (8). Поскольку матрица Я1/2+1 — верхняя треугольная, то можно воспользоваться методом обратного хода для решения систем линейных уравнений [30, с. 258] или [31, с. 51].

3. Новая устойчивая реализация парного фильтра Калмана

Класс квадратно-корневых методов содержит иД-реализации дискретного фильтра Калмана, которые являются модификациями квадратно-корневых фильтров. Идея иД-фильтра принадлежит Дж. Бирману [23]. Она основана на представлении матрицы ковариации в виде Рк\к = иРщкВРщк Ц~рк]к, где иРщк — верхняя треугольная, а ВРщк — диагональная матрицы. Аналогичные обозначения введем для иД-разложения матриц ковариации По, Qx,x и Qyy. Класс [/^-реализаций обеспечивает такую же точность вычислений, как и все квадратно-корневые методы, но при этом позволяет сократить затраты машинного времени и избавляет от необходимости использовать операцию извлечения квадратного корня.

Для преобразования иД-факторов соответствующих матриц в алгоритме фильтрации используем алгоритм модифицированной взвешенной ортогонализации Грама — Шмидта (МВГШО) [32]. Доказано, что МВГШО в вычислительном плане эффективнее, чем классические алгоритмы ортогонализации Грама — Шмидта, при этом точность вычислений по данному алгоритму сравнима с точностью вычислений по алгоритму триангуляризации Хаусхолдера. Что касается реализации МВГШО, то известны два основных вида разложения прямоугольной матрицы А Е Ктхга с помощью МВГШО — это и ^-разложение (обратный порядок вычислений) и ¿^-разложение (прямой порядок вычислений) (детальное описание и доказательство алгоритма можно найти в [23, теорема У1.4.1, с. 127]). В настоящей работе для построения и ^-реализации

парного фильтра использован обратный порядок при реализации МВГШО. Обозначим такой вариант алгоритма как МВГШО-ИБ.

Алгоритм МВГШО-ИВ заключается в ортогонализации столбцов матрицы Лк относительно весовой матрицы и нахождении верхней треугольной матрицы Л\. и диагональной матрицы таких, что выполняется равенство

Лтк VкЛк = Л1VI (4у.

Разработаем новую иД-реализацию для парного фильтра Калмана, опираясь на идеи работ, посвященных и ^-реализациям для классического фильтра Калмана [33]. Новый парный [/Д-фильтр представлен ниже в виде алгоритма 3.

Утверждение 2. Алгоритм 3 алгебраически эквивалентен алгоритму 1.

Доказательство приведено в Приложении.

На последнем этапе фильтрации, т. е. после обработки последнего измерения, вновь имеем Рм= иРщмВРщм1Трм . Величина Кк+111Яек+1, стоящая в круглых скобках

в (14), получена непосредственно из Л\.. В (13) для вычисления оценки вектора состояния требуется обращение матрицы и—1к+1. Поскольку она верхняя треугольная, то можно воспользоваться методом обратного хода для решения систем линейных уравнений [30, с. 258] или [31, с. 51].

Отметим некоторые свойства нового парного иД-фильтра. В алгоритме 3 модифицированное разложение Холецкого осуществляется только один раз для начального значения матрицы Р0. Оно существует, так как вещественная симметрическая матрица Ро > 0 [26], т.е. Ро = иРоПРои]к0, и далее обновление уравнений фильтра происходит в терминах верхней треугольной и и диагональной И матриц.

4. Вычислительные эксперименты

Подтвердим справедливость утверждений 1 и 2 об алгебраической эквивалентности новых алгоритмов парного фильтра Калмана и его стандартной реализации (алгоритм 1) на практике. Проиллюстрируем работоспособность реализаций фильтра на примере решения задачи парной фильтрации из [12].

Пример 1. Рассмотрим парную линейную стохастическую систему вида (1), где

0.18 0.15 0.16" 0.15 0.18 0.14 0.16 0.14 0.18

с начальными условиями х0 = [0.5, 0.5]т и Ро = 2.5 12. Здесь 12 — единичная матрица размера 2.

Проведем следующую серию экспериментов. Прежде всего сгенерируем "истинное" решение парной линейной модели вида (1) с матрицами коэффициентов из примера 1 на отрезке наблюдения [0, 50]. Затем по найденной траектории хк сформируем реализации вектора измерений. После того как смоделированы измерения ук на всем отрезке наблюдения, решаем обратную задачу, т. е. задачу оценивания хщк.

Р Р

1 х,х 1 х,у

0.12 0.10 0.11 0.11 0.10 0.12 0.10 0.11 0.12

Qx,x

От

Qy,y

Алгоритм 3: Парный [/Д-фильтр

Входные данные : Начальные данные фильтра хо, Ро > 0 и результаты измерений

(уо,..., у N}. Предполагается, что Qy,y > 0. Выходные данные: Средние квадратические оценки (х^^}, г = 1,..., N, по

доступным результатам наблюдений (уо,..., у N}. Инициализация фильтра: Выполнить модифицированное разложение Холецкого для

Ро = иРоРр0и'р . Присвоить начальные значения: хо|о = хо,

^0|0 = иРо, ^Р0|0 = °Ро.

Рх,х — Рх,х ^х,у^^—,уРу, х, Рх,у - Рх,у ^х,у^^—,уРу,у,

х,х = Qx,x ^х,у ^—у,у Qy,X.

Выполнить модифицированное разложение Холецкого

ттТ

У,У Гу,у

для матриц ковариации Qy,y = UqViуDQyyу UQy и

0х,х - ир. РГ) .

РХ,Х РХ,Х Гх X

Гсг к = 0 1с N - 1 ас

Этап экстраполяции. Используя полученную в момент времени к оценку Хк|к, предсказать значение хк+1|к для следующего момента времени к + 1:

оценка вектора состояния хк+1|к = Рх^к |к + Ях,у Я—}у у к + Рх,у у к—1.

Вычислить РР-факторы матрицы ковариации Рк+1\к. Сформировать блочные матрицы

А

Т

Fx,xUPklk I U^х , Vk = Bvag{DPklk,DQx^x}

Выполнить МВГШО-преобразование

, МВГШО-UD г

{А^ Vk} _ {UPk+llk ,Dpk+llk}. (10)

Этап обработки измерений и фильтрации. Используя поступившие в момент

времени к + 1 данные у к+1, обновить оценку вектора состояния. Сформировать блочные матрицы

А

Т

UPk+llk 0

Fy,xUpk+iik Щу, Выполнить МВГШО-преобразование

Vk = Diag{Ppk+i|k ,DQy: у }.

г„ _ -.МВГШО-UD + +

{Ak, Vk} ^ {Ak, Vk}, (11)

xk> ^k

где

Vk = Diag{^Pk+i|k+i ,DRe> k+i}, Ak =

Вычислить нормированную невязку измерений

1

UPk+n k+i ^k+\URek+i

0 uRf

e,k + i

(12)

e k+i = U— [yk+i - Fy,xXk+ijk - Fy,yyk] (13

>-e,k+i

и оценку вектора состояния Xk+i|k+i = Xk+i|k + (Kk+iURek+i)ek+i. (14) end

В сравнительном анализе участвуют следующие три реализации фильтра:

• стандартный ПКФ (алгоритм 1);

• квадратно-корневой вариант ПКФ-ЯИ, (алгоритм 2);

• новый [/^-алгоритм ПКФ-ИБ (алгоритм 3).

Подчеркнем, что все реализации фильтра решают задачу из примера 1 в одних и тех же условиях, т. е. их начальные значения и используемые наблюдения одни и те же. Программы для всех численных методов, участвующих в исследовании, написаны на языке матричных вычислений ЫАТЬАБ. Результаты данной серии экспериментов представлены в виде графиков (см. рисунок). Анализируя полученные данные, легко видеть, что все три реализации парного фильтра дают одну и ту же оценку неизвестного вектора состояния парной линейной марковской модели. Это на практическом примере доказывает эквивалентность всех трех реализаций. В частности, показана работоспособность нового [/^-алгоритма. В заключение можно отметить, что оценка вектора состояния близка к "истинному" решению системы.

Пример 2. Рассмотрим задачу оценивания неизвестного вектора состояния парной линейной стохастической системы вида (1), где

F F

х,у

F F

У,У

"0.12 0.10 0.11 0.12'

0.11 0.10 0.12 0.10

1.10 1.10 0.10 0.11

1.10 1.10 + S 0.12 0.10

Q

х,у

Qx,y Qy,y

"0.18 0.15 0.0 0.0

0.15 0.18 0.0 0.0

0.0 0.0 ¿2 0.0

0.0 0.0 0.0 82

и начальные условия х0 = [0.5,0.5]т, Р0 = 2.512. Параметр б2 < еокр, но 6 > е

Величина

;окр-1.

- окр — параметр машинного округления, т. е. еокр + 1 = 1, но еокр/2 + 1

Выбор данной тестовой задачи не случаен. Здесь задачу из примера 1 мы дополнили плохо обусловленной схемой измерений, моделируя тем самым множество плохо обусловленных задач. Другими словами, при малых значениях 5 матрица является плохо обусловленной, а при 5 ^ еокр в результате округления становится вырожденной. Подобная ситуация для классического фильтра Калмана детально описана в [28].

1-1-1-1-1- _i-1-1-1-1-

0 10 20 30 40 к о 10 20 30 40 к

"Истинное" решение парной линейной модели (сплошная линия) и полученная оценка вектора состояния системы: стандартным ПКФ (штриховая линия с маркером о), квадратно-корневым ПКФ-8И (штриховая линия с маркером *) и ^^-фильтром ПКФ-ИБ

(сплошная линия с маркером □)

Накопленная средняя квадратическая ошибка оценивания ARMSE

S ПКФ ПКФ-SR ПКФ-UD S ПКФ ПКФ-SR ПКФ-UD

10-2 0.1797 0.1797 0.1797 10- 10 NaN 0.1792 0.1792

10-3 0.1735 0.1735 0.1735 10- 11 NaN 0.1748 0.1748

10-4 0.1651 0.1651 0.1651 10- 12 NaN 0.1778 0.1778

10-5 0.1793 0.1793 0.1793 10- 13 NaN 0.1711 0.1711

10-6 0.1769 0.1769 0.1769 10- 14 NaN 0.1747 0.1747

10-7 0.1766 0.1766 0.1766 10- 15 NaN 0.1751 0.1751

10-8 NaN 0.1751 0.1751 10- 16 NaN 0.1705 0.1692

10-9 NaN 0.1784 0.1784 10- 17 NaN 0.1787 0.1772

Цель следующей серии экспериментов — практическая оценка точности рассмотренных выше алгоритмов парной фильтрации и их сравнительный анализ на плохо обусловленной задаче примера 2. Итак, для каждого фиксированного значения параметра 8 генерируем "истинное" решение парной линейной модели вида (1) с матрицами коэффициентов из примера 2 на отрезке наблюдения [0,1000]. Затем по найденной траектории Хк сформируем реализации вектора измерений. После того как смоделированы измерения у к на всем отрезке наблюдения, решаем обратную задачу, т. е. задачу оценивания неизвестного вектора Хк|к (к = 1,..., N, N = 1000) парной стохастической модели. Сопоставляя полученные оценки с "истинным" значением вектора состояния Хк, к = 1,... , N, в одних и тех же временных точках, определяем величину ошибки оценивания. Подчеркнем, что все три фильтра используют одни и те же значения "истинного" решения, а также одинаково сгенерированные измерения в целях корректности сравнительного анализа. Поскольку система (1) стохастическая, проводим серии из 100 случайных экспериментов. Для каждой реализации парного фильтра на основе полученного статистического материала рассчитываем накопленную среднюю квадратическую ошибку (ЛИМБЕ) оценивания на всем интервале наблюдения по следующей формуле:

ARMSE

\

1 L N п

LN

xk *-fc|fc)

1=1 k=1 i=1

где верхний индекс г означает i-ю компоненту вектора состояния системы; пх = 2; I — номер эксперимента из L =100 случайных симуляций; к — дискретный момент времени и N — общее количество точек, в которых доступны измерения, т. е. N = 1000.

Расчеты приведенного выше примера выполнялись на персональном компьютере, обеспечивающем представление чисел с относительной погрешностью 10-16. Коды программ написаны на языке MATLAB, где параметр машинного округления еокр хранится в переменной eps со значением еокр = eps/2 = 1.1102 • 10-16. При 8 ^ еокр в результате округления матрица Re, к+1 становится вырожденной, что приводит к остановке алгоритма фильтрации. В системе матричных вычислений MATLAB, на которой написаны все коды программ, это приводит к появлению величины NaN3, которая означает "не число" и является результатом операций œ/œ, 0/0, 0 • œ. Полученные результаты эксперимента занесены в таблицу.

3Сокращение от англ. "Not a Number".

Проанализируем данные, представленные в таблице. Легко видеть, что при 8 > 10-7 все три алгоритма парной фильтрации обеспечивают одинаковую точность вычислений. Ошибка оценивания всех трех реализаций, участвующих в эксперименте, одна и та же. В частности, это еще раз доказывает эквивалентность всех трех методов. Однако при 8 = 10-8 и далее задача становится плохо обусловленной, что сказывается на точности вычислений. Так, при 8 ^ еокр стандартный ПКФ быстрее всех остальных реализаций теряет точность вычислений. Присутствие величины ЫаЫ в строках таблицы — признак полной неспособности стандартного метода справиться с поставленной задачей. Это происходит уже при 8 = 10-8. В то же время квадратно-корневой алгоритм ПКФ-БИ, и построенный в данной работе алгоритм ПКФ-ИВ парной фильтрации легко справляются с поставленной задачей, т. е. это более устойчивые реализации, чем стандартный ПКФ, по отношению к ошибкам округления.

Как отмечалось ранее, реализации ПКФ-БИ, и ПКФ-ИВ позволяют проводить вычисления с двойной точностью. Для данного примера это означает следующее. Легко видеть, что влияние ошибок округления на стандартную реализацию фильтра становится значительным при 8 = 10-8. Следовательно, алгоритмы, сформулированные в терминах квадратного корня или [/^-разложения, могут быть подвержены ощутимому влиянию ошибок округления только при 8 = 10-16, что эквивалентно вычислениям, проводимым с двойной точностью. Это и подтверждается результатами проведенных экспериментов, которые представлены в таблице. Легко видеть, что при 8 > 10-15 оба метода (ПКФ-БИ, и ПКФ-ИВ) справляются с поставленной задачей и их точность вычислений одинакова. Однако при 8 = 10-16 влияние ошибок округления проявляется и в устойчивых реализациях ПКФ. В частности, данная серия экспериментов подтверждает на практике, что построенный в работе новый ПКФ-ИВ-алгоритм несколько превосходит в точности вычислений ранее предложенную квадратно-корневую реализацию, поскольку АИМБЕ для ПКФ-Ш ниже, чем для ПКФ-БИ, при 8 < 10-16.

Заключение

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

В заключение отметим несколько вопросов в области парной калмановской фильтрации, которые до настоящего времени остаются открытыми. Прежде всего, в [17] представлен строгий теоретический анализ влияния ошибок округления на некоторые реализации классического фильтра Калмана. Насколько известно авторам данной статьи, подобный анализ для парного фильтра Калмана не проводился. Кроме того, для классического фильтра Калмана существуют так называемые быстрые алгоритмы фильтрации, основанные на 3-ортогональных преобразованиях. Для парного фильтра Калмана подобные реализации фильтра пока не разработаны. Решение указанных задач может быть продолжением исследований в области парной калмановской фильтрации.

Приложение

Доказательство утверждения 1.

Необходимо показать, что стандартная реализация ПКФ (алгоритм 1) алгебраически эквивалентна квадратно-корневой реализации фильтра (алгоритм 2). Отметим, что каждый этап алгоритма 2 представляет собой уравнение вида 0,А = В. Из свойств ортогональных матриц следует, что ВтВ = Ат 0Т 0.А = АтА. Сопоставляя левую и правую части, легко показать эквивалентность алгоритмов.

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

[Вт В]

Т/2 р1/2 _

р1 /2 р

[л- А]

(1,1) " к+ЦкР к+11к (1,1) ^х,х + Рх,хРк1кР:1.,х.

Аналогично докажем эквивалентность (3)-(5) и (8). Сопоставляя блочные элементы с индексами (1,1) и (1,2) матричных произведений, имеем

[ВТ В] (1,1) = ^!к+1&)!,к+1 = № А] (1,1) = Яу,у + Ру,хрк+11кГг

т

у,х.

Аналогично

[Вт В]

(1,2)

Ет/2 Кт

[лт л]

(1,2)

где

[лт л]

(1,2)

у,у 1/2 в

(0У,у)

Ру,х рТ/2 Рв рк+11к

0 р

1/2

к+Цк

Т

ру,х рТ/2 р1/2 = ру,х р Рв Рк+11кРк+11к = Рв рк+11к.

Отсюда следует, что

КТ

Ке^+лРв'* Рк+Цк,

т. е.

К/,к+1 = Рк+11к (Ре'Х)Т Ре,к+1 = Кк+1Р,

?Т/2 е,к .

Таким образом, эквивалентность (6) и (9) также доказана.

С учетом того, что К/,к+1 = , а также сопоставляя блочные элементы с ин-

дексом (2,2) в матричных произведениях ВтВ и АтА, легко видеть эквивалентность (5) и элемента (2,2) правой части уравнения (8). Утверждение 1 доказано.

Доказательство утверждения 2.

Алгоритм МВГШО заключается в ортогонализации столбцов матрицы Ак относительно весовой матрицы и нахождении верхней треугольной матрицы А\ и диагональной матрицы , таких, что выполняется равенство

4 VI (л

т

'15)

Для доказательства эквивалентности двух алгоритмов необходимо и достаточно доказать эквивалентность выражений (2) и (10), а также (3)-(5) и (11), (12). Эквивалентность выражений (6) и (14) очевидна.

Докажем эквивалентность выражений (2) и (10). Учитывая равенство (15) и выражение (10), запишем

А^ Vk Лк

Fx,xUPkik Uqx

Dp,

k| k

0 D

?т

TTl pl

U Рщк^ X,X

UÇ

QX,X

UPk+i,D

'k+i|kDpk+i|kUpk+i|k.

т

:i6)

Выполнив в (16) умножение матриц, приходим к (2).

Теперь докажем эквивалентность выражений (3)-(5) и (11), (12). Представим произведение Лк Vк Лк в виде

ЛкТ 'Dк Лк

A11 Ai2

A21 A22

:i7)

Определим элементы блочной матрицы (17):

An

A12 A22

uPk+ikDP

k+i|kDpk+i|kUPk+i|k = Рк+1|к

D иТ рТ

Pk+i|kDpk+i|kUPk+i|fc У,х

У

Т

Рк

к+1|к^уТх,

Py,xUPk+i|kDPk+i|kUPk+i|kFy,x + UQy,yDQy,yUQy,y = ^У,хРк+1|кРУ,х + Qy,y.

Аналогично представим произведение AlVl (Дк)т в виде

4 V\ (лк ) =

A A A11 A12

A A A21 A

22

:is)

Определим элементы блочной матрицы (18

A11 = Upk+i|k+iDpk+i|k+iuPk+i|k+i

= Рк+1|к+1 + Кк+1Ре,к+1КТ+1

^e,k + i Te,k + iU Te,k+i k+1

A A1

A

A0

Кk+1UTe,k+lDTe,k+lUTe,k+1 = Кк+1Ле,к+Ь

Т

Kk 11 Re

UTe,k+iDTe,k+iUTe,k+i = ^е,к+1.

Т

Re

Приравнивая левую (17) и правую (18) части равенства АкТ Т>кАк = АкТ>к(Дк)Т, приходим к соотношениям

A11 = AÎ

11

Рк

к+1|к

Рк+1|к+1 + Кк+1Ре,к+1Кк+1 ^

= Рк+1|к — Кк+1Ре,к+1Кк+1

^ Рк+1|к+1 : A12 = A12 ^ Рк+1|кРУ,х = Кk+1Pe,k+1, К,

A22 = A22 ^ Ру,хРк+1|кРУ,х + Qy,y = Ре,к+1 ^ (3).

(5),

k+1|k^y,x = Кк+1Ле,к+Ь Кк+1 = Рk+1|kFУТ,xRe,l+l

(4)

Утверждение 2 доказано.

Благодарности. Первый автор благодарит португальский фонд науки и технологии (Fundaçâo para a Ciencia e a Tecnologia) за финансовую поддержку в рамках гранта № UID/Multi/04621/2013. Второй автор благодарит РФФИ за финансовую поддержку в рамках научного проекта № 16-41-730784.

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

[1] Стратонович Р.Л. Условные марковские процессы и их применение к теории оптимального управления. М.: Изд-во МГУ, 1966. 319 с.

Stratonovic, R.L. Conditional Markov processes and their application to the theory of optimal control. Moscow: Izdat. Moskov. Univ., 1966. 319 p. (In Russ.)

[2] Kalman, R.E. A new approach to linear filtering and prediction problems // J. of Fluids Engineering. 1960. Vol. 82. P. 35-45.

[3] Kalman, R.E., Bucy, R.S. New results in linear filtering and prediction theory // J. of Basic Engineering. 1961. Vol. 83, No. 1. P. 95-108.

[4] Липцер Р.ШШ., ШШиряев А.Н. Статистика случайных процессов (нелинейная фильтрация и смежные вопросы). М.: Наука, 1974. 696 c.

Liptser, R.Sh., Shiryaev, A.N. Statistics of random processes (Nonlinear filtering and related problems) Moscow: Nauka, 1974. 696 p. (In Russ.)

[5] Lanchantin, P., Lapuyade-Lahorgue, J., Pieczynski, W. Unsupervised segmentation of randomly switching data hidden with non-Gaussian correlated noise // Signal Processing. 2011. Vol. 91, No. 2. P. 163-175.

[6] Lipster, R.S., Shiryaev, A.N. Statistics of conditionally gaussian random sequences // Sixth Berkeley Symp. on Mathematics, Statistics and Probability. California: University of California Press, 1972. Vol. 2. P. 389-422.

[7] Pieczynski, W., Desbouvries, F. Kalman filtering using pairwise Gaussian models // Proc. of the IEEE Intern. Conf. on Acoustics, Speech, and Signal Processing, Hong-Kong, 2003. Vol. VI. P. 57-60.

[8] Липцер Р.ШШ. Уравнения почти оптимального фильтра Калмана при особенной матрице ковариаций шума в наблюдениях // Автоматика и телемеханика. 1974. № 1. С. 35-41. Liptser, R.Sh. Equations of near-optimal Kalman filter with a singular matrix of noise covariations in observations // Automation and Remote Control. 1974. Vol. 35, No. 1. P. 29-35.

[9] Ершов А.А., Липцер Р.ШШ. Робастный фильтр Калмана в дискретном времени // Автоматика и телемеханика. 1978. № 3. С. 60-69.

Ershov, A.A., Liptser, R.Sh. A robust Kalman filter in discrete time // Automation and Remote Control. 1978. Vol. 39, No. 3. P. 359-367.

[10] Мильштейн Г.Н., Пьянзин С.А. Регуляризация и цифровое моделирование фильтра Калмана — Бьюси для систем с вырожденными шумами в наблюдениях // Автоматика и телемеханика. 1987. № 11. С. 80-92.

Mil'shtein, G.N., P'yanzin, S.A. Regularization and digital simulation of the Kalman — Bucy filter for systems where the observations are contaminated with singular noise // Automation and Remote Control. 1985. Vol. 46. P. 50-58.

[11] Derrode, S., Pieczynski, W. Signal and image segmentation using pairwise Markov chains // IEEE Transactions on Signal Processing. 2004. Vol. 52, No. 9. P. 2477-2489.

[12] Ait-el-Fquih, B., Desbouvries, F. Unsupervised signal restoration in partially observed Markov chains // Proc. of the IEEE Intern. Conf. on Acoustics, Speech and Signal Processing, Toulouse, France: IEEE, 2006. Vol. 3. P. 13-16.

[13] Derrode, S., Pieczynski, W. Exact fast computation of optimal filter in Gaussian switching linear systems // IEEE Signal Processing Letters. 2013. Vol. 20, No. 7. P. 701-704.

[14] Boudaren, M.E.Y., An, L., Pieczynski, W. Dempster-Shafer fusion of evidential pairwise Markov fields // Intern. J. of Approximate Reasoning. 2016. Vol. 74. P. 13-29.

[15] Boudaren, M.E.Y., Pieczynski, W. Dempster-Shafer fusion of evidential pairwise Markov chains // IEEE Transactions on Fuzzy Systems. 2016. Vol. 24, No 6. P. 1598-1610.

[16] Nemesin, V., Derrode, S. Inferring segmental pairwise Kalman filter with application to pupil tracking // Proc. of the Intern. Work. "Traitement et Analyse de l'Information Methodes et Applications". Hammamet, Tunisia, May 13-18, 2013. P. 1-6.

[17] Verhaegen, M., Van Dooren, P. Numerical aspects of different Kalman filter implementations // IEEE Trans. Automat. Control. 1986. Vol. AC-31, No. 10. P. 907-917.

[18] Nemesin, V., Derrode, S. Robust blind pairwise Kalman algorithms using QR decompositions // IEEE Trans. Signal Process. 2013. Vol. 61, No. 1. P. 5-9.

[19] Nemesin, V., Derrode, S. Robust partial-learning in linear Gaussian systems // IEEE Trans. Automat. Control. 2015. Vol. 60, No. 9. P. 2518-2523.

[20] Dempster, A., Laird, N., Rubin, D. Maximum likelihood from incomplete data via the ЕМ algorithm // J. of the Royal Statistical Society. Series B (Methodological). 1977. Vol. 39, No. 1. P. 1-38.

[21] Kulikova, M.V. Gradient-based parameter estimation in pairwise linear Gaussian system // IEEE Trans. on Automatic Control. 2017. Vol. 62, No. 3. P. 1511-1517.

[22] Kaminski, P.G., Bryson, A.E., Schmidt, S.F. Discrete square-root filtering: a survey of current techniques // IEEE Trans. Automatic Control. 1971. Vol. AC-16, No. 6. P. 727-735.

[23] Bierman, G.J. Factorization methods for discrete sequential estimation. N.Y: Acad. Press, 1977. 241 p.

[24] Potter, J.E., Stern, R.G. Statistical filtering of space navigation measurements // Proc. 1963 AIAA Guidance and Control Conf. N. Y.: AIAA, 1963. 13 p.

[25] Dyer, P., McReynolds, S. Extension of square-root filtering to include process noise // J. Optim. Theory Appl. 1969. Vol. 3, No. 3. P. 444-459.

[26] Голуб Дж., Ван Лоун Ч. Матричные вычисления. М.: Мир, 1999. 548 c.

Golub, G.H., Van Loan, C.F. Matrix computations. Baltimore: Johns Hopkins Univ. Press, 1996. 694 p.

[27] Kailath, T., Sayed, A.H., Hassibi, B. Linear estimation. New Jersey: Prentice Hall, 2000.

854 p.

[28] Grewal, M.S., Andrews, A.P. Kalman filtering: theory and practice. New Jersey: Prentice-Hall, 2001. 401 p.

[29] Park, P., Kailath, T. New square-root algorithms for Kalman filtering // IEEE Trans. Automatic Control. 1995. Vol. 40, No. 5. P. 895-899.

[30] Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.: Наука, 1987. 600 c.

Bakhvalov, N.S., Zhidkov, N.P., Kobel'kov, G.M. Numerical methods. Moscow: Nauka, 1987. 600 p.(In Russ.)

[31] Самарский А.А., Гулин А.В. Численные методы. М.: Наука, 1989. 432 c. Samarskii, A.A., Gulin, A.V. Numerical methods. Moscow: Nauka, 1989. 432 p. (In Russ.)

[32] Bjorck, A. Solving least squares problems by Gram —Schmidt orthogonalization // BIT. 1967. Vol. 7. P. 1-21.

[33] Jover, J.M., Kailath, T. A parallel architecture for kalman filter measurement update and parameter estimation // Automatica. 1986. Vol. 22, No. 1. P. 43-57.

Поступила в 'редакцию 28 октября 2016 г., с доработки — 6 февраля 2017 г.

60

М. В. KyuHKOBa, K. B. ^iraHOBa

Numerically stable Kalman filter implementations for estimating linear pairwise Markov models in the presence of Gaussian noise

Kulikova, Maria V.1, Tsyganova, Julia V.2'*

1CEMAT, Instituto Superior Tecnico, Universidade de Lisboa, 1049-001 Lisboa, Portugal 2Ulyanovsk State University, Ulyanovsk, 432970, Russian

* Corresponding author: Tsyganova, Julia V., e-mail: TsyganovaJV@gmail.com

This paper studies numerical methods of Kalman filtering for vector state estimation of linear Gaussian pairwise models. The pairwise Markov model generalizes the hidden Markov model and it attracted an increasing attention in recent years. For instance, the use of the pairwise Markov models instead of hidden Markov models in segmentation problems allows for dividing the error ratio by two. This paper explores such effective state estimation methods as the square-root pairwise Kalman filtering algorithms, including their array implementations. These filtering methods and their performance are studied in detail. The new UD factorization-based pairwise Kalman filtering approach has been developed. Other filtering algorithms applicable to the linear pairwise Markov model are discussed and numerically compared to the newly developed method on two examples.

Keywords : pairwise Markov model, linear stochastic system, optimal estimation of the state vector, e pairwise Kalman filter, numerically stable (with respect to round off errors) UD-based implementations of the pairwise Kalman filter.

Acknowledgements. The first author acknowledges the support from Portuguese National Funds through the Fundacao para a Ciencia e a Tecnologia (FCT) within project UID/Multi/04621/2013. The second author is grateful to RFBR for the financial support within project No. 16-41-730784.

Received 28 October 2016 Received in revised form 6 February 2017

© ICT SB RAS, 2016