В интернете, в том числе и на хабре, можно найти много информации про фильтр Калмана. Но тяжело найти легкоперевариваемый вывод самих формул. Без вывода вся эта наука воспринимается как некое шаманство, формулы выглядят как безликий набор символов, а главное, многие простые утверждения, лежащие на поверхности теории, оказываются за пределами понимания. Целью этой статьи будет рассказать об этом фильтре на как можно более доступном языке.
Фильтр Калмана — это мощнейший инструмент фильтрации данных. Основной его принцип состоит в том, что при фильтрации используется информация о физике самого явления. Скажем, если вы фильтруете данные со спидометра машины, то инерционность машины дает вам право воспринимать слишком быстрые скачки скорости как ошибку измерения. Фильтр Калмана интересен тем, что в каком-то смысле, это самый лучший фильтр. Подробнее обсудим ниже, что конкретно означают слова «самый лучший». В конце статьи я покажу, что во многих случаях формулы можно до такой степени упростить, что от них почти ничего и не останется.

Ликбез

Перед знакомством с фильтром Калмана я предлагаю вспомнить некоторые простые определения и факты из теории вероятности.

Случайная величина

Когда говорят, что дана случайная величина $\xi$ , то имеют ввиду, что эта величина, может принимать случайные значения. Разные значения она принимает с разной вероятностью. Когда вы кидаете, скажем, кость, то выпадет дискретное множество значений: $\{1,2,3,4,5,6\}$ . Когда речь идет, например, о скорости блуждающей частички, то, очевидно, приходится иметь дело с непрерывным множеством значений. «Выпавшие» значения случайной величины $\xi$ мы будем обозначать через $x_1,x_2...$ , но иногда, будем использовать ту же букву, которой обозначаем случайную величину: $\xi_1, \xi_2,...$ .
В случае с непрерывным множеством значений случайную величину характеризует плотность вероятности $\rho(x)$ , которая нам диктует, что вероятность того, что случайная величина «выпадет» в маленькой окрестности точки $x$ длиной $dx$ равна $\rho(x)dx$ . Как мы видим из картинки, эта вероятность равна площади заштрихованного прямоугольника под графиком:

Довольно часто в жизни случайные величины распределены по Гауссу, когда плотность вероятности равна $\varphi(x)\sim e^{ -\frac{(x-\mu)^2}{2\sigma^2}}$ .

Мы видим, что функция $\rho(x)$ имеет форму колокола с центром в точке $\mu$ и с характерной шириной порядка $\sigma$ .
Раз мы заговорили о Гауссовом распределении, то грешно будет не упомянуть, откуда оно возникло. Также как и числа $e$ и $\pi$ прочно обосновались в математике и встречаются в самых неожиданных местах, так и распределение Гаусса пустило глубокие корни в теорию вероятности. Одно замечательное утверждение, частично объясняющее Гауссово всеприсутствие, состоит в следующем:
Пусть есть случайная величина $\xi$ имеющая произвольное распределение (на самом деле существуют некие ограничения на эту произвольность, но они совершенно не жесткие). Проведем $n$ экспериментов и посчитаем сумму $\xi_1+...+\xi_n$ «выпавших» значений случайной величины. Сделаем много таких экспериментов. Понятно, что каждый раз мы будем получать разное значение суммы. Иными словами, эта сумма является сама по себе случайной величиной со своим каким-то определенным законом распределения. Оказывается, что при достаточно больших $n$ закон распределения этой суммы стремится к распределению Гаусса (к слову, характерная ширина «колокола» растет как $\sqrt{n}$ ). Более подробно читаем в википедии: центральная предельная теорема. В жизни очень часто встречаются величины, которые складываются из большого количества одинаково распределенных независимых случайных величин, поэтому и распределены по Гауссу.

Среднее значение

Среднее значение случайной величины — это то, что мы получим в пределе, если проведем очень много экспериментов, и посчитаем среднее арифметическое выпавших значений. Среднее значение обозначают по-разному: математики любят обозначать через $M\xi$ (математическое ожидание), а заграничные математики через $E\xi$ (expectation). Физики же через $\overline{\xi}$ или $\textless\xi\textgreater$ . Мы будем обозначать на заграничный лад: $E\xi$ .
Например, для Гауссова распределения $\varphi(x)\sim e^{ -\frac{(x-\mu)^2}{2\sigma^2}}$ , среднее значение равно $\mu$ .

Дисперсия

В случае с распределением Гаусса $\varphi(x)\sim e^{ -\frac{(x-\mu)^2}{2\sigma^2}}$ мы совершенно четко видим, что случайная величина предпочитает выпадать в некоторой окрестности своего среднего значения $\mu$ . Как видно из графика, характерный разброс значений порядка $\sigma$ . Как же оценить этот разброс значений для произвольной случайной величины, если мы знаем ее распределение. Можно нарисовать график ее плотности вероятности и оценить характерную ширину на глаз. Но мы предпочитаем идти алгебраическим путем. Можно найти среднюю длину отклонения (модуль) от среднего значения: $E|\xi-E\xi|$ . Эта величина будет хорошей оценкой характерного разброса значений $\xi$ . Но мы с вами очень хорошо знаем, что использовать модули в формулах — одна головная боль, поэтому эту формулу редко используют для оценок характерного разброса.
Более простой способ (простой в смысле расчетов) — найти $E(\xi-E\xi)^2$ . Эту величину называют дисперсией, и часто обозначают как $\sigma_\xi^2$ . Корень из дисперсии называют среднеквадратичным отклонением. Среднеквадратичное отклонение — хорошая оценка разброса случайной величины.
Например, для распределение Гаусса $\varphi(x)\sim e^{ -\frac{(x-\mu)^2}{2\sigma^2}}$ можно посчитать, что определенная выше дисперсия $\sigma_\xi^2$ в точности равна $\sigma^2$ , а значит среднеквадратичное отклонение равно $\sigma$ , что очень хорошо согласуется с нашей геометрической интуицией.
На самом деле тут скрыто маленькое мошенничество. Дело в том, что в определении распределения Гаусса под экспонентой стоит выражение $-\frac{(x-\mu)^2}{2\sigma^2}$ . Эта двойка в знаменателе стоит именно для того, чтобы среднеквадратичное отклонение $\sigma_\xi$ равнялось бы коэффициенту $\sigma$ . То есть сама формула распределения Гаусса написана в виде, специально заточенном для того, что мы будем считать ее среднеквадратичное отклонение.

Независимые случайные величины

Случайные величины бывают зависимыми и нет. Представьте, что вы бросаете иголку на плоскость и записываете координаты ее обоих концов. Эти две координаты зависимы, они связаны условием, что расстояние между ними всегда равно длине иголки, хотя и являются случайными величинами.
Случайные величины независимы, если результат выпадения первой из них совершенно не зависит от результата выпадения второй из них. Если случайные величины $\xi_1$ и $\xi_2$ независимы, то среднее значение их произведения равно произведению их средних значений:
$E(\xi_1\cdot\xi_2)=E\xi_1\cdot E\xi_2$

Доказательство

Например, иметь голубые глаза и окончить школу с золотой медалью — независимые случайные величины. Если голубоглазых, скажем $20\%=0.2$ а золотых медалистов $5\%=0.05$ , то голубоглазых медалистов $0.2\cdot 0.05=0.01=1\%.$ Этот пример подсказывает нам, что если случайные величины $\xi_1$ и $\xi_1$ заданы своими плотностями вероятности $\rho_1(x)$ и $\rho_2(y)$ , то независимость этих величин выражается в том, что плотность вероятности $\rho(x,y)$ (первая величина выпала $x$ , а вторая $y$ ) находится по формуле:
$\rho(x,y)=\rho_1(x)\cdot\rho_2(y)$
Из этого сразу же следует, что:
$\begin{array}{l} \displaystyle E(\xi_1\cdot\xi_2)=\int xy\rho(x,y)dxdy=\int xy\rho_1(x)\rho_2(y)dxdy=\\ \displaystyle \int x\rho_1(x)dx\int y\rho_2(y)dy=E\xi_1\cdot E\xi_2 \end{array}$
Как вы видите, доказательство проведено для случайных величин, которые имеют непрерывный спектр значений и заданы своей плотностью вероятности. В других случаях идея доказательтсва аналогичная.

Фильтр Калмана

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

Обозначим за $x_k$ величину, которую мы будем измерять, а потом фильтровать. Это может быть координата, скорость, ускорение, влажность, степень вони, температура, давление, и т.д.
Начнем с простого примера, который и приведет нас к формулировке общей задачи. Представьте себе, что у нас есть радиоуправляемая машинка, которая может ехать только вперед и назад. Мы, зная вес машины, форму, покрытие дороги и т.д., расcчитали как контролирующий джойстик влияет на скорость движения $v_k$ .

Тогда координата машины будет изменяться по закону:

$x_{k+1}=x_k+v_kdt$

В реальной же жизни мы не можем учесть в наших расчетах маленькие возмущения, действующие на машину (ветер, ухабы, камушки на дороге), поэтому настоящая скорость машины, будет отличаться от расчетной. К правой части написанного уравнения добавится случайная величина $\xi_k$ :

$x_{k+1}=x_k+v_k dt+\xi_k$

У нас есть установленный на машинке GPS сенсор, который пытается мерить истинную координату $x_k$ машинки, и, конечно же, не может ее померить точно, а мерит с ошибкой $\eta_k$ , которая является тоже случайной величиной. В итоге с сенсора мы получаем ошибочные данные:

$z_k=x_k+\eta_k$

Задача состоит в том, что, зная неверные показания сенсора, $z_k$ найти хорошее приближение $x^{opt}_k$ для истинной координаты машины $x_k$ .
В формулировке же общей задачи, за координату $x_k$ может отвечать все что угодно (температура, влажность…), а член, отвечающий за контроль системы извне мы обозначим за $u_k$ (в примере c машиной $u_k=v_k\cdot dt$ ). Уравнения для координаты и показания сенсора будут выглядеть так:

(1)

Давайте подробно обсудим, что нам известно:

$u_k$ — это известная величина, которая контролирует эволюцию системы. Мы ее знаем из построенной нами физической модели.
Ошибка модели $\xi_k$ и ошибка сенсора $\eta_k$ — случайные величины. И их законы распределения не зависят от времени (от номера итерации $k$ ).
Средние значения ошибок равны нулю: $E\xi_k=E\eta_k=0$ .
Сам закон распределения случайных величин может быть нам и не известен, но известны их дисперсии $\sigma_\xi^2$ и $\sigma_\eta^2$ . Заметим, что дисперсии не зависят от $k$ , потому что законы распределения не зависят от него.
Предполагается, что все случайные ошибки независимы друг от друга: какая ошибка будет в момент времени $k$ совершенно не зависит от ошибки в другой момент времени $k'$ .

Нелишним будет отметить, что задача фильтрации — это не задача сглаживания. Мы не стремимся сглаживать данные с сенсора, мы стремимся получить наиболее близкое значение к реальной координате $x_k$ .

Алгоритм Калмана

Мы будем рассуждать по индукции. Представьте себе, что на $k$ -ом шаге мы уже нашли отфильтрованное значение с сенсора $x^{opt}_k$ , которое хорошо приближает истинную координату системы $x_k$ . Не забываем, что мы знаем уравнение, контролирующее изменение нам неизвестной координаты:

$x_{k+1}=x_k+u_k+\xi_k$

поэтому, еще не получая значение с сенсора, мы можем предположить, что на шаге $k+1$ система эволюционирует согласно этому закону и сенсор покажет что-то близкое к $x^{opt}_{k}+u_k$ . К сожалению, пока мы не можем сказать ничего более точного. С другой стороны, на шаге $k+1$ у нас на руках будет неточное показание сенсора $z_{k+1}$ .
Идея Калмана состоит в следующем. Чтобы получить наилучшее приближение к истинной координате $x_{k+1}$ , мы должны выбрать золотую середину между показанием $z_{k+1}$ неточного сенсора и нашим предсказанием $x^{opt}_{k}+u_k$ того, что мы ожидали от него увидеть. Показанию сенсора мы дадим вес $K,$ а на предсказанное значение останется вес $(1-K)$ :

$x^{opt}_{k+1}=K\cdot z_{k+1}+(1-K)\cdot(x^{opt}_k+u_k)$

Коэффициент $K$ называют коэффициентом Калмана. Он зависит от шага итерации, поэтому правильнее было бы писать $K_{k+1}$ , но пока, чтобы не загромождать формулы расчетах, мы будем опускать его индекс.
Мы должны выбрать коэффициент Калмана $K$ таким, чтобы получившееся оптимальное значение координаты $x^{opt}_{k+1}$ было бы наиболее близко к истинной $x_{k+1}$ . К примеру, если мы знаем, что наш сенсор очень точный, то мы будем больше доверять его показанию и дадим значению $z_{k+1}$ больше весу ( $K$ близко единице). Eсли же сенсор, наоборот, совсем не точный, тогда больше будем ориентироваться на теоретически предсказанное значение $x^{opt}_{k}+u_k$ .
В общем случае, чтобы найти точное значение коэффициента Калмана $K$ , нужно просто минимизировать ошибку:

$e_{k+1}=x_{k+1}-x^{opt}_{k+1}$

Используем уравнения (1) (те которые в на голубом фоне в рамочке), чтобы переписать выражение для ошибки:

$\begin{array}{l} e_{k+1}=(1-K)(e_k+\xi_k)-K\eta_{k+1} \end{array}$

Доказательство

$\begin{array}{l} e_{k+1}=x_{k+1}-x^{opt}_{k+1}=x_{k+1}-Kz_{k+1}-(1-K)(x^{opt}_k+u_k)=\\ =x_k+u_k+\xi_k-K(x_k+u_k+\xi_k+\eta_{k+1})-(1-K)(x^{opt}_k+u_k)=\\ =(1-K)(x_k-x_k^{opt}+\xi_k)-K\eta_{k+1}=(1-K)(e_k+\xi_k)-K\eta_{k+1} \end{array}$

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

$E(e^2_{k+1})\rightarrow \mathrm{min}$

Распишем последнее выражение:

$E(e_{k+1}^2)=(1-K)^2(Ee_k^2+\sigma_\xi^2)+K^2\sigma_\eta^2$

Доказательство

Из того что все случайные величины, входящие в выражение для $e_{k+1}$ , независимы, следует, что все «перекрестные» члены равны нулю:
$E(\xi_{k}\eta_{k+1})=E(e_k\xi_k)=E(e_k\eta_{k+1})=0$
Мы использовали тот факт, что $E\eta_{k+1}=E\xi_{k}=0$ , тогда формула для дисперсии выглядит намного проще: $\sigma^2_\xi=E\xi^2$ .

Это выражение принимает минимальное значение, когда(приравниваем производную к нулю):

$K_{k+1}=\frac{Ee_k^2+\sigma_\xi^2}{Ee_k^2+\sigma_\xi^2+\sigma_\eta^2}$

Здесь мы уже пишем выражение для коэффициента Калмана с индексом шага $k+1$ , тем самым мы подчеркиваем, что он зависит от шага итерации.
Подставляем полученное оптимальное значение $K_{k+1}$ в выражение для $E(e^2_{k+1})$ , которую мы минимизировали. Получаем;

$E(e^2_{k+1})=\frac{\sigma_\eta^2(Ee_k^2+\sigma_\xi^2)}{Ee_k^2+\sigma_\xi^2+\sigma_\eta^2}$ .

Наша задача решена. Мы получили итерационную формулу, для вычисления коэффициента Калмана.
Давайте сведем, наши полученные знания в одну рамочку:

Пример

На рекламной картинке вначале статьи отфильтрованы данные с вымышленного GPS сенсора, установленного на вымышленной машине, которая едет равноускоренно c известным вымышленным ускорением $a$ .

$x_{t+1}=x_t+at\cdot dt+\xi_t$

Код на матлабе

clear all; N=100  % number of samples a=0.1 % acceleration sigmaPsi=1 sigmaEta=50; k=1:N x=k x(1)=0 z(1)=x(1)+normrnd(0,sigmaEta); for t=1:(N-1)   x(t+1)=x(t)+a*t+normrnd(0,sigmaPsi);     z(t+1)=x(t+1)+normrnd(0,sigmaEta); end; %kalman filter xOpt(1)=z(1); eOpt(1)=sigmaEta; for t=1:(N-1)   eOpt(t+1)=sqrt((sigmaEta^2)*(eOpt(t)^2+sigmaPsi^2)/(sigmaEta^2+eOpt(t)^2+sigmaPsi^2))   K(t+1)=(eOpt(t+1))^2/sigmaEta^2  xOpt(t+1)=(xOpt(t)+a*t)*(1-K(t+1))+K(t+1)*z(t+1) end; plot(k,xOpt,k,z,k,x)

Анализ

Если проследить, как с шагом итерации $k$ изменяется коэффициент Калмана $K_k$ , то можно показать, что он всегда стабилизируется к определенному значению $K_{stab}$ . К примеру, когда среднеквадратичные ошибки сенсора и модели относятся друг к другу как десять к одному, то график коэффициента Калмана в зависимости от шага итерации выглядит так:

В следующем примере мы обсудим как это поможет существенно облегчить нашу жизнь.

Второй пример

На практике очень часто бывает, что нам вообще ничего не известно о физической модели того, что мы фильтруем. К примеру, вы захотели отфильтровать показания с вашего любимого акселерометра. Вам же заранее неизвестно по какому закону вы намереваетесь крутить акселерометр. Максимум информации, которую вы можете выцепить — это дисперсия ошибки сенсора $\sigma^2_{\eta}$ . В такой непростой ситуации все незнание модели движения можно загнать в случайную величину $\xi_k$ :

$x_{k+1}=x_k+\xi_k$

Но, откровенно говоря, такая система уже совершенно не удовлетворяет тем условиям, которые мы налагали на случайную величину $\xi_k$ , ведь теперь туда запрятана вся неизвестная нам физика движения, и поэтому мы не можем говорить, что в разные моменты времени ошибки модели независимы друг от друга и что их средние значения равны нулю. В этом случае, по большому счету, теория фильтра Калмана не применима. Но, мы не будем обращать внимания на этот факт, а, тупо применим все махину формул, подобрав коэффициенты $\sigma^2_\xi$ и $\sigma^2_\eta$ на глаз, так чтобы отфильтрованные данные миленько смотрелась.
Но можно пойти по другому, намного более простому пути. Как мы видели выше, коэффициент Калмана $K_k$ с увеличением $k$ всегда стабилизируется к значению $K_{stab}$ . Поэтому вместо того, чтобы подбирать коэффициенты $\sigma^2_\xi$ и $\sigma^2_\eta$ и находить по сложным формулам коэффициент Калмана $K_k$ , мы можем считать этот коэффициент всегда константой, и подбирать только эту константу. Это допущение почти ничего не испортит. Во-первых, мы уже и так незаконно пользуемся теорией Калмана, а во-вторых коэффициент Калмана быстро стабилизируется к константе. В итоге все очень упростится. Нам вообще никакие формулы из теории Калмана не нужны, нам просто нужно подобрать приемлемое значение $K_{stab}$ и вставить в итерационную формулу:

$x^{opt}_{k+1}=K_{stab}\cdot z_{k+1}+(1-K_{stab})\cdot x^{opt}_k$

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

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

Обсуждение

Как мы увидели, основная идея фильтра Калмана состоит в том, чтобы найти такой коэффициент $K$ , чтобы отфильтрованное значение

$x^{opt}_{k+1}=Kz_{k+1}+(1-K)(x^{opt}_k+u_k)$

в среднем меньше всего отличалось бы от реального значения координаты $x_{k+1}$ . Мы видим, что отфильтрованное значение $x^{opt}_{k+1}$ есть линейная функция от показания сенсора $z_{k+1}$ и предыдущего отфильтрованного значения $x^{opt}_k$ . А предыдущее отфильтрованное значение $x^{opt}_k$ является, в свою очередь, линейной функцией от показания сенсора $z_k$ и предпредыдущего отфильтрованного значения $x^{opt}_{k-1}$ . И так далее, пока цепь полностью не развернется. То есть отфильтрованное значение зависит от всех предыдущих показаний сенсора линейно:

$x^{opt}_{k+1}=\lambda+\lambda_0z_0+...+\lambda_{k+1}z_{k+1}$

Поэтому фильтр Калмана называют линейным фильтром.
Можно доказать, что из всех линейных фильтров Калмановский фильтр самый лучший. Самый лучший в том смысле, что средний квадрат ошибки фильтра минимален.

Многомерный случай

Всю теорию фильтра Калмана можно обобщить на многомерный случай. Формулы там выглядят чуть страшнее, но сама идея их вывода такая же, как и в одномерном случае. В этой прекрасной статье вы можете увидеть их: http://habrahabr.ru/post/140274/.
А в этом замечательном видео разобран пример, как их использовать.

Литература

Оригинальную статью Калмана можно скачать вот тут: http://www.cs.unc.edu/~welch/kalman/media/pdf/Kalman1960.pdf.

ссылка на оригинал статьи http://habrahabr.ru/post/166693/

Фильтр Калмана

Ликбез

Случайная величина

Среднее значение

Дисперсия

Независимые случайные величины

Фильтр Калмана

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

Алгоритм Калмана

Пример

Анализ

Второй пример

Обсуждение

Многомерный случай

Литература

Комментарии

Добавить комментарий Отменить ответ