Подходы к выбору значимой плотности вероятности в фильтрах частиц (particle filters) Текст научной статьи по специальности «Физика»

Доклады БГУИР

2012 № 4 (66)

УДК 621.396.96

ПОДХОДЫ К ВЫБОРУ ЗНАЧИМОЙ ПЛОТНОСТИ ВЕРОЯТНОСТИ В ФИЛЬТРАХ ЧАСТИЦ (PARTICLE FILTERS)

А.В. ПАРАХНЕВИЧ, А С. СОЛОНАР, С.А. ГОРШКОВ

Военная академия Республики Беларусь Минск-57, 220057, Беларусь

Поступила в редакцию 20 апреля 2012

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

Ключевые слова: значимая плотность вероятности (importance density), метод Монте-Карло, обобщенный фильтр частиц (particle filter), фильтр Калмана.

Введение

При изучении алгоритмов работы фильтров частиц [1-6] одним из важных вопросов является выбор значимой плотности вероятности (ПВ) [2-4, 7, 8], которая аппроксимирует требуемую апостериорную плотность вероятности при помощи метода численного интегрирования Монте-Карло [1, 9]. Правильный выбор значимой плотности ведет к снижению ошибок нахождения моментов аппроксимируемой плотности вероятности и повышает точность использования метода Монте-Карло [1] в задачах дискретной байесовской фильтрации.

Наиболее точным является случай, когда дисперсия значимых весов оказывается минимальной (оптимальный подход) [1, 9]. Однако такой подход является достаточно сложным в вычислительном плане. Более простым в описании является подход, когда все плотности вероятности аппроксимируются гауссовыми (квазиоптимальный гауссов подход) [10]. Сложность возникает при расчете ковариационных матриц ошибок. И сложность тем выше, чем больше мерность вектора фильтруемых параметров. Тем не менее, наиболее часто используемым является случай, когда в качестве значимой плотности используется гауссова переходная плотность вероятности (модифицированный квазиоптимальный гауссов подход) [4].

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

Задачей дискретной нелинейной фильтрации является вычисление апостериорной плотности вероятности иа-мерного вектора состояния цели а при наличии и8-мерного вектора наблюдаемых параметров 8.

Введем следующие обозначения. Под Ак = {ак, ак_р..., а1} будем понимать обобщенный вектор состояния объекта в момент времени tk с учетом последовательности всех предыдущих состояний объекта а;- (у = 1, к^ . Совокупность векторов наблюдаемых параметров за к шагов наблюдения будем обозначать как 0к = { 8к,8к_1,...,81} .

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

Для решения поставленной задачи рассмотрим следующие варианты выбора значимой ПВ: оптимальный выбор значимой ПВ, квазиоптимальный гауссов и модифицированный квазиоптимальный гауссов методы. Сопоставительный анализ работы обобщенного фильтра частиц с перевыборкой весовых коэффициентов (SIR PF) и фильтра Калмана [11] был проведен для задающего воздействия для случая прямолинейно движущегося летательного аппарата, выполняющего разворот по окружности и затем продолжающего движение по прямолинейной траектории.

Значимая плотность вероятности

Одной из основных проблем в разработке фильтров частиц является выбор значимой плотности вероятности q(ak | ак8k), где i - номер частицы = 1, N^ . Как показано в [4, 5], оптимальная значимая плотность вероятности, минимизирующая дисперсию значимых весов, должна быть полностью равна аппроксимации апостериорной ПВ p(ak | ак_1,8k) [1, 4, 9]. В таком случае она будет определяться выражением:

q(ak |ak_1,8к)„ = p(aк |ak_,,8к) = p(8kK«■k|ak^. ш

P(8k|ak-1)

Как уже было рассмотрено в [1], нормированные веса частиц можно получить при помощи выражения:

< х wk[ p(8k|aak)/>(«k|a k _,) . (2)

q(ak| a k8k)

Подстановка (1) в (2) и упрощение p(8k | a'к,ak_j) ~ p(8k | ak) ([5], стр. 61) приведут к следующему выражению для расчета весов на текущем шаге:

wk х wk-ЛК_x). (3)

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

1) провести выборку случайных отсчетов из апостериорной ПВ p(ak | a'k-р8k) [4, 5];

2) вычислить для каждой i-й частицы значение условной функции правдоподобия вектора наблюдения:

P(8k | ak-1) = j P(8k | ak)P(ak | ak-1)dak. (4)

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

плотности вероятности p(ak | a'k-1,8k) становится возможной. Примером использования является случай фильтрации вектора состояния линейной системы со скачкообразным переключением модели движения при сопровождении маневрирующих целей. Второй случай: работа с

классом моделей, для которых плотность p(ak | ak-1,8k) является гауссовой ([5], стр. 62).

Для случая квазиоптимальной гауссовой аппроксимации ПВ примем: модель задающего воздействия с нелинейной регулярной частью fk_1(ak_1); модель входного воздействия (модель наблюдения) - линейной, а независимые аддитивные случайные элементы задающего и возму-

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

« * = Ак_1(а к-1) + V к _1, (5)

вк = hk(а к ) + ^ к = Н к а* + w к, (1)

где ук-1 и wk - взаимно независимые последовательности белых гауссовых шумов с нулевым математическим ожиданием и с ковариационными матрицами Мк-1 и Rв соответственно; ^(ак) -

функция пересчета вектора состояния ак в вектор наблюдения вк.

Нк - статическая матрица пересчета вектора состояния в вектор наблюдаемых параметров. Элементы матрицы определяются через частные производные:

Н [■, 7] = ^; (■ = 1,., Пв; 7 = 1,..., Па ), (7)

к [ 7 ]

где а к [ 7] - 7-й элемент вектора состояния; h к [7] - 7-я строка вектор функции h к (а к).

Для этого случая в [3, 4] было показано, что и оптимальная плотность вероятности р(ак | ак_1, вк), и априорная функция правдоподобия р(вк | ак_1) являются гауссовыми, и представляются как

р(ак | ак_1, вк) = N(аk ^ _1), (8)

р(вк|ак_1) = N (вк ;Ьк ,L к), (9)

где ак = fkч(ак_1) + Я_1Н 1 Яв1(вк -Ьк) - математическое ожидание ПВ р(ак | ак_1, вк);

R _1 = м к _1 _ Мк _1Н 1ь_;и к м к_1 ковариационная матрица ПВ р(а к 1 а к _р вк); ь к = Н к fk _1(а к _1)

- математическое ожидание ПВ р(вк | ак_1); Lк = НкМк1Н1 + Явк - ковариационная матрица

пв р(вкК 1).

В этом случае переходная ПВ и функция правдоподобия будут определяться следующими параметрами (с учетом допущений, относящихся к (5) и (6)):

р(ак1 ак_1) = N(ак 'А_1(ак_1),мк_l), р(вк I а к) = N (вк Н к а к Двк),

На практике такое аналитическое описание плотностей р(ак | ак_1, вк) и р(вк | ак_1), а также операции над матрицами, являются достаточно сложными. Это, в свою очередь, приводит к использованию модифицированного квазиоптимального гауссова подхода для аппроксимации оптимальных плотностей. Наиболее распространенной является аппроксимация значимой плотности переходной плотностью вероятности типа:

Ч(ак 1 ак_Р вк) = р(а к 1 а к _1). (10)

С учетом введенных обозначений переходная плотность запишется в виде:

р(ак I ак_1) = N(аk ^(а'^М^). (11)

В свою очередь, выражение (3) для расчета ненормированных весов частиц w'k на текущем временном шаге будет определяться значением функции правдоподобия в координатах частиц, экстраполированных с предыдущего шага а^ = Ак_1 (а7к 1) + ук_1. Пересчет осуществляется с учетом случайной составляющей ук_1 и значением веса частицы на предыдущем шаге

wkœ wk-1p(8k I ak )■

(12)

Если переходная плотность p(ak | ak-1) используется как значимая плотность вероятности с более широким распределением, чем функция правдоподобия p(8k | ak ), то в результате перемножения плотностей лишь нескольким частицам будут присвоены большие веса (рис. 1). Это приводит к быстрому вырождению частиц и к неработоспособности фильтра. Аналогичная проблема возникает вследствие работы фильтра с генерацией выборки весовых коэффициентов (SIS), рассмотренная в [2-4, 7].

В 1999-2003 гг. были разработаны методы для коррекции координат частиц с установкой их в «правильные» места (в зоны высокой вероятности) посредством объединения текущего наблюдения и результатов экстраполяции частиц. Наиболее широкое применение получил метод, реализуемый во вспомогательном фильтре частиц с перевыборкой (AUX PF) [4, 12].

В ряде случаев наиболее эффективными для квазиоптимальных аппроксимаций оптимальных плотностей вероятности могут быть приемы локальной линеаризации [7, 13]. Такая линеаризация может применяться для гауссовой значимой ПВ и реализуется в алгоритмах работы расширенного фильтра Калмана (extended Kalman filter) [6] и нецентральном фильтре Калмана (unscented Kalman filter) [6, 13].

Рис.1. Пример попадания одной частицы при аппроксимации функции правдоподобия для общего числа частиц Ж=10

Из известных публикаций выяснилось, что реализация эффективного фильтра частиц в большей степени искусство, чем наука. В настоящее время не существует четких теоретических правил для синтеза фильтра частиц, в то время как их использование требует выполнения процедуры перевыборки на каждой итерации или при периодическом выявлении вырождения. Кроме того, вопрос аналитического определения минимального числа частиц способных осуществлять статистически эффективную реализацию процедуры фильтрации, остается открытым. Ответить на него можно пока лишь опытным путем при проведении экспериментальных исследований. Немногочисленные исследования для случая неманеврирующих целей (раздел 6.5 в [6]) показывают, что увеличение числа частиц более, чем до ^=5000, приводят лишь к незначительному возрастанию эффективности алгоритмов фильтрации.

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

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

Пример работы фильтров частиц

Для проверки работоспособности фильтров частиц было проведено сопоставительное моделирование работы фильтра Калмана второго порядка и обобщенного фильтра частиц с оптимальной значимой ПВ. Входным воздействием являлись разовые оценки параметров траектории летательного аппарата, состоящей из трех участков (рис. 2). Первый и третий участки представляют собой прямолинейные траектории движения с постоянной скоростью. Второй участок движения - разворот по окружности радиусом 5000 м и нормальной перегрузкой Скорость полета на всем участке составляла 200 м/с, начальная точка (первая разовая оценка)

имеет дальность 120 км и азимут 330 градусов, начальный угол курса 90 градусов. Интервал обновления данных (период обзора) Т=5 с. Число разовых оценок М=120. Оценкам с 40-й по 68-ю соответствует участок движения по окружности. Вектор наблюдаемых параметров включал разовые оценки дальности г и азимута в: =|\гк вк| |Т. Затем осуществлялся пересчет разовых оценок в прямоугольную систему координат и вектор наблюдаемых параметров приобретал вид 02к =|\хк zk\|Т (ось Ох направлена на север, Oz - на восток). Вектор состояния фильтра Калмана включал в себя прямоугольные координаты х, z, скорости их изменения Vx,

Vz и ускорения ax, az: ak = \xk zk Vx Vz

a

. Вектор состояния фильтра частиц

SIR включал прямоугольные координаты x и z, полный вектор скорости V и угол курса ф:

Vk Ф

Координата Z, м

Рис. 2. Последовательность разовых оценок, образующих задающее воздействие Модель движения цели была линеаризована:

« к = ^к-1) + V к _1 = вка к-1 + V к _l, (13)

где Вк - динамическая матрица пересчета вектора состояния с (к-1) на к шаг измерения [14]; V к-1 - случайная составляющая модели задающего воздействия с ковариационной матрицей случайного маневра Мк-1 = SкSтам , распределенного по нормальному закону с нулевым МО и СКО ам =0,01 м/с ([14], стр. 346).

В модели наблюдения ковариационная матрица ошибок наблюдения определялась СКО ошибок разового оценивая дальности аг =70 м, азимута ар =0,2 градусов, значениями текущей дальности г до объекта и азимута в на объект:

R„

Пересчет ошибок измерения в прямоугольную систему координат осуществлялся при помощи матриц перехода, подробно описанных в ([14], стр. 354). Ковариационная матрица ошибок наблюдения для прямоугольных координат с учетом текущих значений дальности и азимута вектора наблюдаемых параметров 0Х имеет вид:

2 a r 0

0 2 aP

R

(arcosPk)2 + (aerk sinPk)

а;: sin рkcosPk " аrk sin PkcosP a2 sin PkcosPk - а2r2 sin PkcosPk (ark cosPk)2+(ал sinPk)2

Для сопоставительного сравнения анализировались следующие зависимости:

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

- ошибки фильтрации скорости (рис. 4) в фильтре Калмана определялись как корень суммы квадратов разности задающего значения скорости и отфильтрованного по каждой координате. Невязка скорости для фильтра частиц SIR определялась как модуль разности между полным вектором скорости задающего воздействия и отфильтрованного значения скорости;

T

a

x

k

T

a k = x

z

k

k

k

- ошибки фильтрации угла курса (рис. 5) определялась как разность между истинным и отфильтрованным значением угла курса. Значение угла курса в фильтре Калмана рассчитывалось на основании отфильтрованных оценок составляющих вектора скорости Ух и У2.

Данные ошибок усреднялись по 2000 экспериментам по каждой разовой оценке.

450-

0 10 20 30 40 50 60 70 80 90 100 110 120

Номер шага наблюдения, \

Рис. 3. Усредненная ошибка измерения положения объекта на плоскости разовых оценок фильтра Калмана и обобщенного фильтра частиц

Рис. 4. Усредненная ошибка измерения полной скорости объекта фильтра Калмана и обобщенного

фильтра частиц

I 0. 6

га и

& 0.5 >

> 0.4

* 0.2

X

"V

га 0.1

х

п

о ►

3 10 20 30 40 50 60 70 80 90 100 110 120

Номер шага наблюдения, i

Рис. 5. Усредненная ошибка измерения угла курса объекта фильтра Калмана и обобщенного фильтра

частиц

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

Так как фильтры Калмана охватывают линейную часть возможных задающих воздействий, то при фильтрации вектора состояния с нелинейным пересчетом параметров движения и при нелинейной взаимосвязи вектора состояния и вектора наблюдаемых параметров, фильтры частиц, как видно, имеют несомненные преимущества. Это подтверждается опубликованными в зарубежных источниках результатами сопоставительного анализа [7, 10, 13].

Заключение

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

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

APPROACES TO THE IMPORTANCE DENSITY CHOICE IN PARTICLE FILTERS

A.V. PARAKHNEVICH, A S. SOLONAR, S.A. GORSHKOV

Abstract

There are spent three methods of importance density choice (Gaussian, kvasi-Gaussian and modified kvasi-Gaussian) for posterior density approximation using Monte Carlo numerical integration method and considered comparing Kalman filter and Particle Filter implementation in the report.

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

1. Парахневич А.В., СолонарА.С., Горшков С.А. // Докл. БГУИР. 2012. №1(63). С. 22-28.

2. Парахневич А.В., СолонарА.С., Горшков С.А. // Докл. БГУИР. 2012. №з(б5). С. 49-55.

3. Chen Z. // IEEE A&E Systems Magazine. 2011. №4. P. 69.

4. Ristic B., Arulampalam S., Gordon N. Beyond the Kalman Filter. Particle filters for tracking applications. London, 2004.

5. Gustafsson F. // IEEE A&E SYSMEMS MAGAZINE. 2010. Vol. 25, №7. P. 53-81.

6. Daum F. // IEEE A&E Systems Magazine. 2005. Vol. 20, №8. P. 57-69.

7. Doucet A. // Signal Processing Group, Department of Engineering University of Cambridge. Technical report CUED/F-INFENG/TR.310. 1988.

8. Doucet A., Godsill S., Andrieu C. // Statistics and Computing. 2000. Vol. 10, №3. P. 197-208.

9. Соболь И.М. Численные методы Монте-Карло. М., 1973.

10. Bolic M. Architectures for Efficient Implementation of Particle Filters. Dissertation of Ph.D, Stony Brook University, 2004.

11. Kalman R.E. // Transactions of the ASME. 1960. P. 35-45.

12. PittM., ShephardN. // Journal of the American Statistical association. 1999. P. 590-599.

13. Merwe R., Doucet A., Freitas N. et al. // Tech. Rep. CUED/F-1NFENG/TR 380. Cambridge University Engineering Department, 2000.

14. Ширман Я.Д., Багдасарян С.Т., Горшков С.А. и др. Радиоэлектронные системы: Основы построения и теория. М., 2007.