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

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

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

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

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

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

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

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

$x_k=f_k(x_{k-<!--\; -\; -->1},v_k)$$x_{k}=f_{k}(x_{k-1},v_{k})$ ,

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

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

$y_k=h_k(x_k,w_k)$$y_{k}=h_{k}(x_{k},w_{k})$ ,

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

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

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

Скрытая марковская модель и байесовский выводПравить

Основная статья: Скрытая марковская модель

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

$X_1\sim <!--\; \sim \; -->\mathrm{\mu <!--\; \mu \; -->}(x_1)$$X_{1}\sim \mu (x_{1})\quad$  и $X_n\mid <!--\; \mid \; -->(X_{n-<!--\; -\; -->1}=x_{n-<!--\; -\; -->1})\sim <!--\; \sim \; -->f(x_n\mid <!--\; \mid \; -->x_{n-<!--\; -\; -->1})$$X_{n}\mid (X_{n-1}=x_{n-1})\sim f(x_{n}\mid x_{n-1})$ ,

(1)

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

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

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

$Y_n\mid <!--\; \mid \; -->(X_n=x_n)\sim <!--\; \sim \; -->h(y_n\mid <!--\; \mid \; -->x_n)$$Y_{n}\mid (X_{n}=x_{n})\sim h(y_{n}\mid x_{n})$ ,

(2)

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

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

${\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}}$ 

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

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

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

$p(x_{1:n})=\mathrm{\mu <!--\; \mu \; -->}(x_1)\underset{k=2}{\overset{n}{\prod <!--\; \prod \; -->}}f(x_k\mid <!--\; \mid \; -->x_{k-<!--\; -\; -->1})$$p(x_{1:n})=\mu (x_{1})\prod _{k=2}^{n}f(x_{k}\mid x_{k-1})$

(3)

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

$p(y_{1:n}\mid <!--\; \mid \; -->x_{1:n})=\underset{k=1}{\overset{n}{\prod <!--\; \prod \; -->}}h(y_k\mid <!--\; \mid \; -->x_k)$$p(y_{1:n}\mid x_{1:n})=\prod _{k=1}^{n}h(y_{k}\mid x_{k})$ ,

(4)

Здесь и далее нотация $x_{k:l}$$x_{k:l}$  для $k⩽<!--\; ⩽\; -->l$$k\leqslant l$  обозначает $(x_k,\dots <!--\; \dots \; -->,x_l)$$(x_{k},\dots ,x_{l})$ .

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

$p(x_{1:n}\mid <!--\; \mid \; -->y_{1:n})=\frac{p(x_{1:n})p(y_{1:n}\mid <!--\; \mid \; -->x_{1:n})}{p(y_{1:n})}$$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)

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

$p(y_{1:n})=\int <!--\; \int \; -->p(x_{1:n})p(y_{1:n}\mid <!--\; \mid \; -->x_{1:n})d{x}_{1:n}$$p(y_{1:n})=\int p(x_{1:n})p(y_{1:n}\mid x_{1:n})\,dx_{1:n}$ .

Выборка по значимостиПравить

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

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

$\stackrel{¯<!--\; ¯\; -->}{\mathrm{\theta <!--\; \theta \; -->}}=\int <!--\; \int \; -->\mathrm{\theta <!--\; \theta \; -->}(x)p(x)dx$${\bar {\theta }}=\int \theta (x)p(x)\,dx$ ,

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

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

$\stackrel{¯<!--\; ¯\; -->}{\mathrm{\theta <!--\; \theta \; -->}}\approx <!--\; \approx \; -->\frac{1}{N}\underset{i=1}{\overset{N}{\sum <!--\; \sum \; -->}}\mathrm{\theta <!--\; \theta \; -->}(x^{(i)})$${\bar {\theta }}\approx {\frac {1}{N}}\sum _{i=1}^{N}\theta (x^{(i)})$ 

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

$w_i=\frac{r(x^{(i)})}{\underset{j=1}{\overset{N}{\sum <!--\; \sum \; -->}}r(x^{(j)})}$$w_{i}={\frac {r(x^{(i)})}{\sum _{j=1}^{N}r(x^{(j)})}}$ 

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

$\stackrel{¯<!--\; ¯\; -->}{\mathrm{\theta <!--\; \theta \; -->}}=\int <!--\; \int \; -->\mathrm{\theta <!--\; \theta \; -->}(x)r(x)q(x)dx\approx <!--\; \approx \; -->\underset{i=1}{\overset{N}{\sum <!--\; \sum \; -->}}{w}_i\mathrm{\theta <!--\; \theta \; -->}(x^{(i)})$${\bar {\theta }}=\int \theta (x)r(x)q(x)\,dx\approx \sum _{i=1}^{N}w_{i}\theta (x^{(i)})$ ,

ПеревыборкаПравить

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

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

Последовательный метод Монте-КарлоПравить

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

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

Если записать плотность как:

$f_n(x_{1:n})=\frac{{\mathrm{\varphi <!--\; \varphi \; -->}}_n(x_{1:n})}{Z_n}$$f_{n}(x_{1:n})={\frac {\phi _{n}(x_{1:n})}{Z_{n}}}$ , где

${\mathrm{\varphi <!--\; \varphi \; -->}}_n:<!--\; :\; -->{\mathcal{X}}^n\to <!--\; \to \; -->{\mathbb{R}}^{+}$$\phi _{n}\colon {\mathcal {X}}^{n}\to \mathbb {R} ^{+}$  известна поточечно, а

$Z_n=\int <!--\; \int \; -->{\mathrm{\varphi <!--\; \varphi \; -->}}_n(x_{1:n})d{x}_{1:n}$$Z_{n}=\int \phi _{n}(x_{1:n})\,dx_{1:n}$  — нормализующая, возможно неизвестная, постоянная, то

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

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

${\mathrm{\varphi <!--\; \varphi \; -->}}_n(x_{1:n})=p(x_{1:n})p(y_{1:n}\mid <!--\; \mid \; -->x_{1:n})$$\phi _{n}(x_{1:n})=p(x_{1:n})p(y_{1:n}\mid x_{1:n})$  и

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

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

$f_n(x_{1:n})=\frac{p(x_{1:n})p(y_{1:n}\mid <!--\; \mid \; -->x_{1:n})}{p(y_{1:n})}=p(x_{1:n}|y_{1:n})$$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})$ .

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

$p(x_{1:n}\mid <!--\; \mid \; -->y_{1:n-<!--\; -\; -->1})=p(x_{1:n-<!--\; -\; -->1}\mid <!--\; \mid \; -->y_{1:n-<!--\; -\; -->1})f(x_n\mid <!--\; \mid \; -->x_{n-<!--\; -\; -->1})$$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})$  — предиктор,

$p(x_{1:n}\mid <!--\; \mid \; -->y_{1:n})=\frac{h(y_n\mid <!--\; \mid \; -->x_n)p(x_{1:n}\mid <!--\; \mid \; -->y_{1:n-<!--\; -\; -->1})}{p(y_n\mid <!--\; \mid \; -->y_{1:n-<!--\; -\; -->1})}$$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})}}$  — корректор.

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

АлгоритмПравить

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

   Алгоритм МЧФ
   -- инициализация
   для i = 1...N:
     выборка ${\mathrm{\xi <!--\; \xi \; -->}}_0^{(i)}$$\xi _{0}^{(i)}$  из $q_0(x_0\mid <!--\; \mid \; -->y_0)$$q_{0}(x_{0}\mid y_{0})$ 
     -- начальные веса
     ${\mathrm{\omega <!--\; \omega \; -->}}_0^{(i)}:=h(y_0\mid <!--\; \mid \; -->{\mathrm{\xi <!--\; \xi \; -->}}_0^{(i)})\mathrm{\mu <!--\; \mu \; -->}({\mathrm{\xi <!--\; \xi \; -->}}_0^{(i)})/q_0({\mathrm{\xi <!--\; \xi \; -->}}_0^{(i)}\mid <!--\; \mid \; -->y_0)$$\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:
     если ПЕРЕВЫБОРКА то
       -- выбор индексов $j_i\in <!--\; \in \; -->\{1,\dots <!--\; \dots \; -->,N\}$$j_{i}\in \{1,\dots ,N\}$  N частиц в соответствии с весами
       $j_{1:N}$$j_{1:N}$  = SelectByWeight($\{w_{n-<!--\; -\; -->1}^{(j)}\}$$\{w_{n-1}^{(j)}\}$ )
       для i = 1...N:
         $x_{n-<!--\; -\; -->1}^{(i)}:={\mathrm{\xi <!--\; \xi \; -->}}_{n-<!--\; -\; -->1}^{(j_i)}$$x_{n-1}^{(i)}:=\xi _{n-1}^{(j_{i})}$ 
         $w_{n-<!--\; -\; -->1}^{(i)}:=1/N$$w_{n-1}^{(i)}:=1/N$ 
     иначе
       для i = 1...N:
         $x_{n-<!--\; -\; -->1}^{(i)}:={\mathrm{\xi <!--\; \xi \; -->}}_{n-<!--\; -\; -->1}^{(i)}$$x_{n-1}^{(i)}:=\xi _{n-1}^{(i)}$ 
     для i = 1...N:
       -- шаг распространения частицы
       ${\mathrm{\xi <!--\; \xi \; -->}}_n^{(i)}\sim <!--\; \sim \; -->q_n({\mathrm{\xi <!--\; \xi \; -->}}_n^{(i)}\mid <!--\; \mid \; -->{\mathrm{\xi <!--\; \xi \; -->}}_{n-<!--\; -\; -->1}^{(i)},y_n)$$\xi _{n}^{(i)}\sim q_{n}(\xi _{n}^{(i)}\mid \xi _{n-1}^{(i)},y_{n})$ 
       -- обновление весов
       ${\mathrm{\omega <!--\; \omega \; -->}}_n^{(i)}:=w_{n-<!--\; -\; -->1}^{(i)}h(y_n\mid <!--\; \mid \; -->{\mathrm{\xi <!--\; \xi \; -->}}_n^{(i)})f({\mathrm{\xi <!--\; \xi \; -->}}_n^{(i)}\mid <!--\; \mid \; -->x_{n-<!--\; -\; -->1}^{(i)})/q_n({\mathrm{\xi <!--\; \xi \; -->}}_n^{(i)}\mid <!--\; \mid \; -->x_{n-<!--\; -\; -->1}^{(i)},y_n)$$\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})$  
     кц
     -- нормализация весов
     $s:=\underset{j=1}{\overset{N}{\sum <!--\; \sum \; -->}}{\mathrm{\omega <!--\; \omega \; -->}}_n^{(j)}$$s:=\sum _{j=1}^{N}\omega _{n}^{(j)}$ 
     для i = 1...N:
       $w_n^{(i)}:={\mathrm{\omega <!--\; \omega \; -->}}_n^{(i)}/s$$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