Многочастичный фильтр

Многочасти́чный фильтр^[1] (МЧФ, англ. particle filter — «фильтр частиц», «частичный фильтр», «корпускулярный фильтр») — последовательный метод Монте-Карло — рекурсивный алгоритм для численного решения проблем оценивания (фильтрации, сглаживания), особенно для нелинейных и не-гауссовских случаев. Со времени описания в 1993 году^[2] Н. Гордоном, Д. Салмондом и А. Смитом используется в различных областях — навигации, робототехнике, компьютерном зрении.

В сравнении с обычно применяемыми для подобных задач методами — расширенными фильтрами Кальмана (EKF) — многочастичные фильтры не зависят от методов линеаризации или апроксимации. Обычный EKF плохо справляется с существенно нелинейными моделями, а также в случае шумов системы и измерений, сильно отличающихся от гауссовых, поэтому были разработаны различные модификации, такие как UKF (англ. unscented KF), QKF (англ. Quadrature KF) и т. п.^[3]. Следует отметить, что в свою очередь многочастичные фильтры более требовательны к вычислительным ресурсам.

Термин «particle filter» был дан Дел Моралом в 1996 году^[4], а «sequential Monte Carlo» — Лю (Liu) и Ченом (Chen) в 1998.

Многие используемые на практике многочастичные фильтры выводятся применением последовательного метода Монте-Карло к последовательности целевых распределений^[5].

МЧФ предназначен для оценки последовательности скрытых переменных ${\displaystyle x_{n}}$ для ${\displaystyle n=1,2,\dots }$ на основании наблюдений ${\displaystyle y_{n}}$ при ${\displaystyle n=1,2,\dots }$. Для простоты изложения будем считать, что рассматривается динамическая система, и ${\displaystyle x_{n}}$ и ${\displaystyle y_{n}}$ — действительные вектора состояния и измерений соответственно^[1].

Стохастическое уравнение состояния системы имеет вид:

${\displaystyle x_{k}=f_{k}(x_{k-1},v_{k})}$,

где ${\displaystyle f_{k}}$ функция изменения состояния системы, ${\displaystyle v_{k}}$ — случайная величина, возмущающее воздействие.

Уравнение измерений:

${\displaystyle y_{k}=h_{k}(x_{k},w_{k})}$,

где ${\displaystyle h_{k}}$ функция измерения, ${\displaystyle w_{k}}$ — случайная величина, шум измерений.

Функции ${\displaystyle f_{k}}$ и ${\displaystyle h_{k}}$ в общем случае нелинейные, а статистические характеристики шума системы (${\displaystyle v_{k}}$) и измерений (${\displaystyle w_{k}}$) предполагаются известными.

Задачей фильтрации является получение оценки ${\displaystyle {\hat {x}}_{k}}$ на основе известных к моменту ${\displaystyle k}$ результатов измерений ${\displaystyle y_{1:k}}$.

Рассмотрим дискретный марковский процесс ${\displaystyle \{X_{n}\}_{n\geqslant 1}}$ со следующими распределениями вероятностей:

${\displaystyle X_{1}\sim \mu (x_{1})\quad }$ и ${\displaystyle X_{n}\mid (X_{n-1}=x_{n-1})\sim f(x_{n}\mid x_{n-1})}$,

(1)

где ${\displaystyle \mu (x)}$ — плотность вероятности, ${\displaystyle f(x_{n}\mid x_{n-1})}$ — условная плотность вероятности (переходная плотность вероятности) при переходе от ${\displaystyle x_{n-1}}$ к ${\displaystyle x_{n}}$.

Здесь нотация ${\displaystyle X\mid Y\sim f(\dots )}$ означает, что ${\displaystyle X}$ при условии ${\displaystyle Y}$ распределено как ${\displaystyle f(\dots )}$.

Реализации процесса ${\displaystyle \{X_{n}\}}$ (скрытые переменные ${\displaystyle x_{n}}$) наблюдаются посредством другого случайного процесса ${\displaystyle \{Y_{n}\}_{n\geqslant 1}}$ — процесса измерений — с маргинальными плотностями:

${\displaystyle Y_{n}\mid (X_{n}=x_{n})\sim h(y_{n}\mid x_{n})}$,

(2)

где ${\displaystyle h(y_{n}\mid x_{n})}$ — условная плотность вероятности (плотность измерений), измерения считаются статистически независимыми.

Модель может проиллюстрирована следующей диаграммой переходов:

${\displaystyle {\begin{array}{cccccccccc}X_{1}&\rightarrow &X_{2}&\rightarrow &X_{3}&\rightarrow &X_{4}&\rightarrow &\ldots &\\\downarrow &&\downarrow &&\downarrow &&\downarrow &&\ldots &\\Y_{1}&&Y_{2}&&Y_{3}&&Y_{4}&&\ldots &\end{array}}}$

Для простоты считаем, что переходная плотность и плотность измерений не зависят от ${\displaystyle n}$. Параметры модели считаются заданными.

Определённая таким образом модель системы и измерений известна как скрытая марковская модель^[6].

Уравнение (1) определяет априорное распределение для процесса ${\displaystyle \{X_{n}\}}$:

${\displaystyle p(x_{1:n})=\mu (x_{1})\prod _{k=2}^{n}f(x_{k}\mid x_{k-1})}$

(3)

Аналогично (2) задаёт функцию правдоподобия:

${\displaystyle p(y_{1:n}\mid x_{1:n})=\prod _{k=1}^{n}h(y_{k}\mid x_{k})}$,

(4)

Здесь и далее нотация ${\displaystyle x_{k:l}}$ для ${\displaystyle k\leqslant l}$ обозначает ${\displaystyle (x_{k},\dots ,x_{l})}$.

Таким образом, байесовский вывод для ${\displaystyle \{X_{1:n}\}}$ при известных реализациях измерений ${\displaystyle \{Y_{1:n}\}}$, обозначенных соответственно как ${\displaystyle \{x_{1:n}\}}$ и ${\displaystyle \{y_{1:n}\}}$, будет опираться на апостериорное распределение

${\displaystyle p(x_{1:n}\mid y_{1:n})={\frac {p(x_{1:n})p(y_{1:n}\mid x_{1:n})}{p(y_{1:n})}}}$,

(5)

где (здесь ${\displaystyle dx_{1:n}}$ — доминирующая мера):

${\displaystyle p(y_{1:n})=\int p(x_{1:n})p(y_{1:n}\mid x_{1:n})\,dx_{1:n}}$.

См. также Выборка по значимости.

Метод Монте-Карло позволяет оценивать свойства довольно сложных распределений вероятностей, например, путём вычисления средних и дисперсии в виде интеграла^[3]:

${\displaystyle {\bar {\theta }}=\int \theta (x)p(x)\,dx}$,

где ${\displaystyle \theta (x)}$ — функция для оценивания. Например, для среднего можно положить: ${\displaystyle \theta (x)=x}$.

В случае невозможности аналитического решения, задача может быть решена численно генерированием случайных выборок с плотностью ${\displaystyle p(x)}$, обозначим их как ${\displaystyle {x^{(i)}}_{1\leqslant i\leqslant N}}$, и получением среднего арифметического по точкам выборки^[3]:

${\displaystyle {\bar {\theta }}\approx {\frac {1}{N}}\sum _{i=1}^{N}\theta (x^{(i)})}$

В более общем случае, когда выборка из ${\displaystyle p}$ затруднена, применяется другое распределение ${\displaystyle q}$ (так называемое англ. instrumental or importance distribution), а для сохранения несмещённости оценки вводятся весовые коэффициенты ${\displaystyle w_{i}}$ на основе отношения ${\displaystyle r(x^{(i)})=p(x^{(i)})/q(x^{(i)})}$^[3]:

${\displaystyle w_{i}={\frac {r(x^{(i)})}{\sum _{j=1}^{N}r(x^{(j)})}}}$

после чего вычисляет взвешенное среднее:

${\displaystyle {\bar {\theta }}=\int \theta (x)r(x)q(x)\,dx\approx \sum _{i=1}^{N}w_{i}\theta (x^{(i)})}$,

Хотя вспомогательное распределение используется в основном для упрощения выборки из основного распределения ${\displaystyle p}$, часто применяется процедура «выборки и перевыборки по значимости» (англ. sampling importance resampling, SIR). Эта процедура состоит из двух этапов: собственно выборки по значимости с вычислением весов ${\displaystyle w_{i}}$, и дополнительной выборки точек, учитывающих эти веса^[3].

Перевыборка особенно необходима для последовательных фильтров^[3].

Методы многочастичной фильтрации и сглаживания являются наиболее известными примерами алгоритмов последовательного метода Монте-Карло (англ. sequential Monte Carlo, SMC). До такой степени, что в литературе часто не делают между ними различия. Тем не менее, SMC включает в себя более широкий класс алгоритмов, применимых для описания более сложных приблизительных методов фильтрации и сглаживания^[7].

Последовательные методы Монте-Карло являются классом методов Монте-Карло, которые производят последовательную выборку из последовательности целевых плотностей вероятностей ${\displaystyle \{f_{n}(x_{1:n})\}}$ увеличивающейся размерности, где каждое ${\displaystyle f_{n}(x_{1:n})}$ определено на декартовой степени ${\displaystyle {\mathcal {X}}^{n}}$^[5].

Если записать плотность как:^[5]

${\displaystyle f_{n}(x_{1:n})={\frac {\phi _{n}(x_{1:n})}{Z_{n}}}}$, где

${\displaystyle \phi _{n}\colon {\mathcal {X}}^{n}\to \mathbb {R} ^{+}}$ известна поточечно, а

${\displaystyle Z_{n}=\int \phi _{n}(x_{1:n})\,dx_{1:n}}$ — нормализующая, возможно неизвестная, постоянная, то

SMC-алгоритм будет находить приближения ${\displaystyle f_{k}(x_{1:k})}$ и оценки ${\displaystyle Z_{k}}$ для ${\displaystyle k=1,2,\dots }$.

Например, для случая фильтрации можно положить (см. (5)):

${\displaystyle \phi _{n}(x_{1:n})=p(x_{1:n})p(y_{1:n}\mid x_{1:n})}$ и

${\displaystyle Z_{n}=p(y_{1:n})}$,

из чего будем иметь:

${\displaystyle f_{n}(x_{1:n})={\frac {p(x_{1:n})p(y_{1:n}\mid x_{1:n})}{p(y_{1:n})}}=p(x_{1:n}|y_{1:n})}$.

Опуская вывод, схему предиктор-корректор можно представить в следующем виде^[3]:

${\displaystyle p(x_{1:n}\mid y_{1:n-1})=p(x_{1:n-1}\mid y_{1:n-1})f(x_{n}\mid x_{n-1})}$ — предиктор,

${\displaystyle p(x_{1:n}\mid y_{1:n})={\frac {h(y_{n}\mid x_{n})p(x_{1:n}\mid y_{1:n-1})}{p(y_{n}\mid y_{1:n-1})}}}$ — корректор.

Множитель ${\displaystyle (p(y_{n}\mid y_{1:n-1}))^{-1}}$ — нормализующая постоянная, которая не требуется для обычного SMC-алгоритма.

Типичный алгоритм многочастичного фильтра можно представить в следующем виде^[3]:

   Алгоритм МЧФ
   -- инициализация
   для i = 1...N:
     выборка ${\displaystyle \xi _{0}^{(i)}}$ из ${\displaystyle q_{0}(x_{0}\mid y_{0})}$
     -- начальные веса
     ${\displaystyle \omega _{0}^{(i)}:=h(y_{0}\mid \xi _{0}^{(i)})\mu (\xi _{0}^{(i)})\ /\ q_{0}(\xi _{0}^{(i)}\mid y_{0})}$ 
   кц
   для n = 1...T:
     если ПЕРЕВЫБОРКА то
       -- выбор индексов ${\displaystyle j_{i}\in \{1,\dots ,N\}}$ N частиц в соответствии с весами
       ${\displaystyle j_{1:N}}$ = SelectByWeight(${\displaystyle \{w_{n-1}^{(j)}\}}$)
       для i = 1...N:
         ${\displaystyle x_{n-1}^{(i)}:=\xi _{n-1}^{(j_{i})}}$
         ${\displaystyle w_{n-1}^{(i)}:=1/N}$
     иначе
       для i = 1...N:
         ${\displaystyle x_{n-1}^{(i)}:=\xi _{n-1}^{(i)}}$
     для i = 1...N:
       -- шаг распространения частицы
       ${\displaystyle \xi _{n}^{(i)}\sim q_{n}(\xi _{n}^{(i)}\mid \xi _{n-1}^{(i)},y_{n})}$
       -- обновление весов
       ${\displaystyle \omega _{n}^{(i)}:=w_{n-1}^{(i)}h(y_{n}\mid \xi _{n}^{(i)})f(\xi _{n}^{(i)}\mid x_{n-1}^{(i)})\ /\ q_{n}(\xi _{n}^{(i)}\mid x_{n-1}^{(i)},y_{n})}$ 
     кц
     -- нормализация весов
     ${\displaystyle s:=\sum _{j=1}^{N}\omega _{n}^{(j)}}$
     для i = 1...N:
       ${\displaystyle w_{n}^{(i)}:=\omega _{n}^{(i)}/s}$
   кц

• Фильтр Кальмана#UKF

• Doucet Arnaud, Johansen Adam M. A Tutorial on Particle Filtering and Smoothing: Fifteen Years Later // The Oxford Handbook of Nonlinear Filtering / D. Crisan, B. Rozovsky. — Oxford : Oxford University Press, 2011. — P. 656—704. — ISBN 978-0-19-953290-2.
• Cappé, Olivier and Godsill, Simon J. and Moulines, Eric. An Overview of Existing Methods and Recent Advances in Sequential Monte Carlo // Proceedings of the IEEE. — IEEE, 2007. — Т. 95, № 5. — P. 899—924. — ISSN 0018-9219. Архивировано 10 марта 2016 года.

• Doucet, Arnaud and de Freitas, Nando and Gordon, Neil. An Introduction to Sequential Monte Carlo Methods // Sequential Monte Carlo Methods in Practice / Doucet, Arnaud and de Freitas, Nando and Gordon, Neil. — Springer New York. — 3-14 p. — ISBN 978-1-4419-2887-0.
• Arulampalam, M.S. and Maskell, S. and Gordon, N. and Clapp, T. A Tutorial on Particle Filters for Online Nonlinear/non-Gaussian Bayesian Tracking (англ.) // Trans. Sig. Proc.. — IEEE Press, 2002. — Vol. 50, no. 2. — P. 174—188. — ISSN 1053-587X. См. также более раннюю версию (англ.)
• Gordon, N.J.; Salmond, D.J.; Smith, A.F.M. Novel approach to nonlinear/non-Gaussian Bayesian state estimation (англ.) // IEEE Proceedings F, Radar and Signal Processing. — IET, 1993. — Vol. 140, no. 2. — P. 107—113. — doi:10.1049/ip-f-2.1993.0015.
• Микаэльян С. В. Методы фильтрации на основе многоточечной аппроксимации плотности вероятности оценки в задаче определения параметров движения цели при помощи измерителя с нелинейной характеристикой // Наука и образование : электронное издание. — МГТУ им. Н. Э. Баумана, 2011. — ISSN 1994-0408. Архивировано 4 марта 2016 года.
• Ristic, B., Arulampalam, S., Gordon, N. Beyond the Kalman Filter — Particle Filters for Tracking Applications. — Artech House, 2004. — 299 p. — ISBN 9781580536318.

• Simon, Dan. 15 The particle filter // Optimal State Estimation: Kalman, H_∞, and Nonlinear Approaches. — Wiley-Interscience, 2006. — P. 461—480. — ISBN 0471708585.

• Particle Filter, SciPy Cookbook