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


Этот фильтр применяют в разных областях – от радиотехники до экономики. Здесь мы обсудим основную идею, смысл, суть данного фильтра. Излагаться она будет максимально простым языком.
Предположим, что у нас есть необходимость в измерениях некоторых величин некоего объекта. В радиотехнике чаще всего имеют дело с измерениями напряжений на выходе некоего устройства (датчика, антенны и т.д.). В примере с электрокардиографом (см. ) мы имеем дело с измерениями биопотенциалов на теле человека. В экономике, например, измеряемой величиной могут быть курсы валют. Каждыё день курс валют разный, т.е. каждый день “его измерения” дают нам разную величину. А если обобщать, то можно сказать, что большая часть деятельности человека (если не вся) сводится именно к постоянным измерениям-сравнениям тех или иных величин (см. книгу).
Итак, предположим, что мы что-то постоянно измеряем. Так же предположим, что наши измерения всегда идут с некоторой ошибкой – оно и понятно, ведь нет идеальных измерительных приборов, и каждый выдаёт результат с ошибкой. В простейшем случае описанное можно свести к следующему выражению: z=x+y, где x – истинное значение, которое мы хотим измерить и которое измерили бы если бы у нас был идеальный измерительный прибор, y – ошибка измерения, вносимая измерительным прибором, а z – измеренная нами величина. Так вот задача фильтра Калмана состоит в том, чтобы по измеренной нами z всё-таки догадаться (определить), а какое же истинное значение x было, когда мы получали нашу z (в которой "сидит" истинное значение и ошибка измерения). Необходимо отфильтровать (отсеять) из z истинное значение x – убрать из z искажающий шум y. То есть, имея на руках только лишь сумму нам необходимо догадаться о том, какие слагаемые дали эту сумму.
В свете вышеописанного сформулируем теперь всё следующим образом. Пусть есть всего лишь два случайных числа. Нам даётся только их сумма и от нас требуется по этой сумме определить, какими являются слагаемые. Например, нам дали число 12 и говорят: 12 – это сумма чисел x и y, вопрос – чему равны x и y. Чтобы ответить на этот вопрос, составляем уравнение: x+y=12. Мы получили одно уравнение с двумя неизвестными, поэтому, строго говоря, найти два числа которые и дали эту сумму не возможно. Но кое-что об этих числах мы всё-таки можем сказать. Мы можем сказать, что это были либо числа 1 и 11, либо 2 и 10, либо 3 и 9, либо 4 и 8 и т.д., также это либо 13 и -1, либо 14 и -2, либо 15 и -3 и т.д. То есть мы можем по сумме (в нашем примере 12) определить множество возможных вариантов, которые дают в сумме именно 12. Один из этих вариантов – это искомая нами пара, которая на самом деле прямо сейчас и дала 12. Нелишне так же отметить, что все варианты пар чисел дающих в сумме 12 образуют прямую, изображённую на рис.1, которая и задаётся уравнением x+y=12 (y=-x+12).

Рис.1

Таким образом, искомая нами пара лежит где-то на этой прямой. Повторюсь, выбрать из всех этих вариантов ту пару, которая была на самом деле – которая дала число 12, не владея какими-либо дополнительными подсказками, невозможно. Однако, в ситуации, для которой изобретён фильтр Калмана, такие подсказки есть . Там заранее о случайных числах кое-что известно. В частности там известна так называемая гистограмма распределения для каждой пары чисел. Она обычно бывает получена после достаточно длительных наблюдений за выпадениями этих самых случайных чисел. То есть, например, из опыта известно, что в 5% случаев обычно выпадает пара x=1, y=8 (обозначим эту пару так: (1,8)), в 2% случаев пара x=2, y=3 (2,3), в 1% случаев пара (3,1), в 0.024% случаев пара (11,1) и т.д. Повторюсь, эта гистограмма задана для всех пар чисел, в том числе и для тех, что образуют в сумме 12. Таким образом, для каждой пары, что даёт в сумме 12, мы можем сказать, что, например, пара (1, 11) выпадает в 0.8% случаев, пара (2, 10) – в 1% случаев, пара (3, 9) – в 1.5% случаев и т.д. Таким образом, мы можем по гистограмме определить, в скольких процентах случаев сумма слагаемых пары равна 12. Пусть, например, в 30% случаев сумма даёт 12. А в остальных 70% выпадают остальные пары – это (1,8), (2,3), (3,1) и т.д. – те, что в сумме дают числа отличные от 12. Причём пусть, например, пара (7,5) выпадает в 27% случаев в то время, как все остальные пары, что дают в сумме 12, выпадают в 0.024%+0.8%+1%+1.5%+…=3% случаев. Итак, по гистограмме мы выяснили, что числа дающие в сумме 12 выпадают в 30% случаев. При этом мы знаем, что если выпало 12, то чаще всего (в 27% из 30%) причиной этого является пара (7,5). То есть если уже выпало 12, то мы можем сказать, что в 90% (27% из 30% – или, что то же самое 27 раз из каждых 30-ти) причиной выпадения 12 является пара (7,5). Зная, что чаще всего причиной получения суммы равной 12 является пара (7,5) логично предположить, что, скорее всего, она выпала и сейчас. Конечно, всё-таки не факт, что на самом деле сейчас число 12 образовано именно этой парой, однако, в следующие разы, если нам попадётся 12, и мы опять предположим пару (7,5), то где-то в 90% случаев из 100% окажемся правы. А вот если мы будем предполагать пару (2, 10), то окажемся правы лишь в 1% из 30% случаев, что равно 3.33% правильных догадок по сравнению с 90% при предположении пары (7,5). Вот и всё – в этом и состоит смысл алгоритма фильтра Калмана. То есть фильтр Калмана не гарантирует, что не ошибётся в определении слагаемого по сумме, однако он гарантирует, что ошибётся минимальное количество раз (вероятность ошибки будет минимальна), так как использует статистику – гистограмму выпадения пар чисел. Так же необходимо подчеркнуть, что часто в алгоритме фильтрации Калмана используется так называемая плотность распределения вероятности (ПРВ). Однако необходимо понимать, что смысл там тот же, что и у гистограммы. Более того, гистограмма – это функция, построенная на основе ПРВ и являющаяся её приближением (см., например, ).
В принципе мы эту гистограмму можем изобразить в виде функции двух переменных – то есть в виде некоей поверхности над плоскостью xy. Там, где поверхность выше, там выше и вероятность выпадения соответствующей пары. На рис.2 изображена такая поверхность.


рис.2

Как видно над прямой x+y=12 (которая есть варианты пар дающих в сумме 12) расположены точки поверхности на разной высоте и наибольшая высота у варианта с координатами (7,5). И когда нам встречается сумма равная 12, в 90% случаев причиной появления этой суммы является именно пара (7,5). Т.е. именно эта пара, дающая в сумме 12, имеет наибольшую вероятность появления при условии, что сумма равна 12.
Таким образом, здесь описана идея лежащая в основе фильтра Калмана. Именно на ней и построены всевозможные его модификации – одношаговые, многошаговые рекуррентные и т.д. Для более глубокого изучения фильтра Калмана рекомендую книгу: Ван Трис Г. Теория обнаружения, оценок и модуляции.

p.s. Для того, кто интересуется объяснениями понятий математики что называется "на пальцах" можно посоветовать вот эту книгу и в частности главы из её раздела "Математика" (саму книгу или отдельные главы из неё вы можете приобрести ).

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

Ликбез

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

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

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

Довольно часто в жизни случайные величины распределены по Гауссу, когда плотность вероятности равна .

Мы видим, что функция имеет форму колокола с центром в точке и с характерной шириной порядка .
Раз мы заговорили о Гауссовом распределении, то грешно будет не упомянуть, откуда оно возникло. Также как и числа и прочно обосновались в математике и встречаются в самых неожиданных местах, так и распределение Гаусса пустило глубокие корни в теорию вероятности. Одно замечательное утверждение, частично объясняющее Гауссово всеприсутствие, состоит в следующем:
Пусть есть случайная величина имеющая произвольное распределение (на самом деле существуют некие ограничения на эту произвольность, но они совершенно не жесткие). Проведем экспериментов и посчитаем сумму «выпавших» значений случайной величины. Сделаем много таких экспериментов. Понятно, что каждый раз мы будем получать разное значение суммы. Иными словами, эта сумма является сама по себе случайной величиной со своим каким-то определенным законом распределения. Оказывается, что при достаточно больших закон распределения этой суммы стремится к распределению Гаусса (к слову, характерная ширина «колокола» растет как ). Более подробно читаем в википедии: центральная предельная теорема . В жизни очень часто встречаются величины, которые складываются из большого количества одинаково распределенных независимых случайных величин, поэтому и распределены по Гауссу.

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

Среднее значение случайной величины - это то, что мы получим в пределе, если проведем очень много экспериментов, и посчитаем среднее арифметическое выпавших значений. Среднее значение обозначают по-разному: математики любят обозначать через (математическое ожидание), а заграничные математики через (expectation). Физики же через или . Мы будем обозначать на заграничный лад: .
Например, для Гауссова распределения , среднее значение равно .

Дисперсия

В случае с распределением Гаусса мы совершенно четко видим, что случайная величина предпочитает выпадать в некоторой окрестности своего среднего значения . Как видно из графика, характерный разброс значений порядка . Как же оценить этот разброс значений для произвольной случайной величины, если мы знаем ее распределение. Можно нарисовать график ее плотности вероятности и оценить характерную ширину на глаз. Но мы предпочитаем идти алгебраическим путем. Можно найти среднюю длину отклонения (модуль) от среднего значения: . Эта величина будет хорошей оценкой характерного разброса значений . Но мы с вами очень хорошо знаем, что использовать модули в формулах - одна головная боль, поэтому эту формулу редко используют для оценок характерного разброса.
Более простой способ (простой в смысле расчетов) - найти . Эту величину называют дисперсией, и часто обозначают как . Корень из дисперсии называют среднеквадратичным отклонением. Среднеквадратичное отклонение - хорошая оценка разброса случайной величины.
Например, для распределение Гаусса можно посчитать, что определенная выше дисперсия в точности равна , а значит среднеквадратичное отклонение равно , что очень хорошо согласуется с нашей геометрической интуицией.
На самом деле тут скрыто маленькое мошенничество. Дело в том, что в определении распределения Гаусса под экспонентой стоит выражение . Эта двойка в знаменателе стоит именно для того, чтобы среднеквадратичное отклонение равнялось бы коэффициенту . То есть сама формула распределения Гаусса написана в виде, специально заточенном для того, что мы будем считать ее среднеквадратичное отклонение.

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

Случайные величины бывают зависимыми и нет. Представьте, что вы бросаете иголку на плоскость и записываете координаты ее обоих концов. Эти две координаты зависимы, они связаны условием, что расстояние между ними всегда равно длине иголки, хотя и являются случайными величинами.
Случайные величины независимы, если результат выпадения первой из них совершенно не зависит от результата выпадения второй из них. Если случайные величины и независимы, то среднее значение их произведения равно произведению их средних значений:

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

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

Из этого сразу же следует, что:

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

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

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

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

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

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

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

Задача состоит в том, что, зная неверные показания сенсора, найти хорошее приближение для истинной координаты машины .
В формулировке же общей задачи, за координату может отвечать все что угодно (температура, влажность…), а член, отвечающий за контроль системы извне мы обозначим за (в примере c машиной ). Уравнения для координаты и показания сенсора будут выглядеть так:

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

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

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

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

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

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

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

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

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

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

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

Из того что все случайные величины, входящие в выражение для , независимы, следует, что все «перекрестные» члены равны нулю:

Мы использовали тот факт, что , тогда формула для дисперсии выглядит намного проще: .

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

Здесь мы уже пишем выражение для коэффициента Калмана с индексом шага , тем самым мы подчеркиваем, что он зависит от шага итерации.
Подставляем полученное оптимальное значение в выражение для , которую мы минимизировали. Получаем;

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

Пример

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

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)

Анализ

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

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

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

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

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

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

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

Обсуждение

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

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

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

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

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

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

Одним из эффективных алгоритмов фильтрации является фильтр Калмана. Фильтр Калмана - рекурсивный фильтр, оценивающий вектор состояния динамической системы, используя ряд неполных и зашумленных измерений . Фильтр Калмана использует динамическую модель системы (к примеру, физический закон движения), управляющие воздействия и множество последовательных измерений для формирования оптимальной оценки состояния. Алгоритм состоит из двух повторяющихся фаз: предсказание и корректировка. На первом этапе рассчитывается предсказание состояния в последующий момент времени (с учетом неточности их измерения). На втором, новая информация с датчика корректирует предсказанное значение (также с учетом неточности и зашумленности этой информации).

На этапе предсказания происходит:

  1. Предсказание состояния системы:

где – предсказание состояния системы в текущий момент времени; – матрица перехода между состояниями (динамическая модель системы); – предсказание состояния системы в предыдущий момент времени; – матрица применения управляющего воздействия; – управляющее воздействие в предыдущий момент времени.

  1. Предсказание ошибки ковариации:

где – предсказание ошибки; – ошибка в предыдущий момент времени; – ковариация шума процесса.

На этапе корректировки происходит:

  1. Вычисление усиления Калмана:

где – усиление Калмана; – матрица измерений, отображающая отношение измерений и состояний; – ковариация шума измерения.

где – измерение в текущий момент времени.

  1. Обновление ошибки ковариации:

где – матрица идентичности.

Если состояние системы описывается одной переменной, то = 1, а матрицы вырождаются в обычные уравнения.

Чтобы наглядно продемонстрировать эффективность фильтра Калмана, был проведён эксперимент с датчиком громкости AVR PIC KY-037, который подключен к микроконтроллеру Arduino Uno. На рисунке 1 представлен график показаний датчика без применения фильтра (линия 1). Хаотичные колебания значения на выходе датчика свидетельствуют о наличии шумов.

Рисунок 1. График показаний датчика без применения фильтра

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

Расчёт переменной произведём в программе Microsoft Excel 2010. Для этого необходимо рассчитать среднеквадратичное отклонение для выборки значений показаний датчика. = 0,62. подбирается в зависимости от требуемого уровня фильтрации, принимаем = 0,001. На рисунке 2 второй линией представлен график показаний датчика с применением фильтра.

Рисунок 2. График показаний датчика с применением фильтра Калмана

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

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

Рисунок 3. Реакция фильтра Калмана на скачок измеряемой величины

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

Из проведённого эксперимента можно сделать вывод о том, что фильтр Калмана целесообразно применять для фильтрации показаний датчиков в системах с низким быстродействием.

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

  1. ГОСТ 17657-79. Передача данных. Термины и определения. – Москва: Изд-во стандартов, 2005. – 2 с.
  2. Фильтр Калмана // Википедия. . Дата обновления: 26.04.2017. URL: http://ru.wikipedia.org/?oldid=85061599 (дата обращения: 21.05.2017).

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

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

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

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

Задача оценки параметров

Одной из задач теории статистических решений, имеющих большое практическое значение, является задача оценки векторов состояния и параметров систем, которая формулируется следующим образом. Предположим, необходимо оценить значение векторного параметра $X$, недоступного непосредственному измерению. Вместо этого измеряется другой параметр $Z$, зависящий от $X$. Задача оценивания состоит в ответе на вопрос: что можно сказать об $X$, зная $Z$. В общем случае, процедура оптимальной оценки вектора $X$ зависит от принятого критерия качества оценки.

Например, байесовский подход к задаче оценки параметров требует полной априорной информации о вероятностных свойствах оцениваемого параметра, что зачастую невозможно. В этих случаях прибегают к методу наименьших квадратов (МНК), который требует значительно меньше априорной информации.

Рассмотрим применения МНК для случая, когда вектор наблюдения $Z$ связан с вектором оценки параметров $X$ линейной моделью, и в наблюдении присутствует помеха $V$, некоррелированная с оцениваемым параметром:

$Z = HX + V$, (1)

где $H$ – матрица преобразования, описывающая связь наблюдаемых величин с оцениваемыми параметрами.

Оценка $X$, минимизирующая квадрат ошибки, записывается следующим образом:

$X_{оц}=(H^TR_V^{-1}H)^{-1}H^TR_V^{-1}Z$, (2)

Пусть помеха $V$ не коррелирована, в этом случае матрица $R_V$ есть просто единичная матрица, и уравнение для оценки становится проще:

$X_{оц}=(H^TH)^{-1}H^TZ$, (3)

Запись в матричной форме сильно экономит бумагу, но может быть для кого то непривычна. Следующий пример, взятый из монографии Коршунова Ю. М. "Математические основы кибернетики", все это иллюстрирует.
Имеется следующая электрическая цепь:

Наблюдаемые величины в данном случае – показания приборов $A_1 = 1 A, A_2 = 2 A, V = 20 B$.

Кроме того, известно сопротивление $R = 5$ Ом. Требуется оценить наилучшим образом, с точки зрения критерия минимума среднего квадрата ошибки значения токов $I_1$ и $I_2$. Самое важное здесь заключается в том, что между наблюдаемыми величинами (показаниями приборов) и оцениваемыми параметрами существует некоторая связь. И эта информация привносится извне.

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

Итак, знание законов Кирхгофа, никак не связанное с теорией статистических решений, позволяет установить связь между наблюдаемыми значениями и оцениваемыми параметрами (кто изучал электротехнику – могут проверить, остальным придется поверить на слово):

$$z_1 = A_1 = I_1 + \xi_1 = 1$$

$$z_2 = A_2 = I_1 + I_2 + \xi_2 = 2$$

$$z_2 = V/R = I_1 + 2 * I_2 + \xi_3 = 4$$

Это же в векторной форме:

$$\begin{vmatrix} z_1\\ z_2\\ z_3 \end{vmatrix} = \begin{vmatrix} 1 & 0\\ 1 & 1\\ 1 & 2 \end{vmatrix} \begin{vmatrix} I_1\\ I_2 \end{vmatrix} + \begin{vmatrix} \xi_1\\ \xi_2\\ \xi_3 \end{vmatrix}$$

Или $Z = HX + V$, где

$$Z= \begin{vmatrix} z_1\\ z_2\\ z_3 \end{vmatrix} = \begin{vmatrix} 1\\ 2\\ 4 \end{vmatrix} ; H= \begin{vmatrix} 1 & 0\\ 1 & 1\\ 1 & 2 \end{vmatrix} ; X= \begin{vmatrix} I_1\\ I_2 \end{vmatrix} ; V= \begin{vmatrix} \xi_1\\ \xi_2\\ \xi_3 \end{vmatrix}$$

Считая значения помехи некоррелированными между собой, найдем оценку I 1 и I 2 по методу наименьших квадратов в соответствии с формулой 3:

$H^TH= \begin{vmatrix} 1 & 1& 1\\ 0 & 1& 2 \end{vmatrix} \begin{vmatrix} 1 & 0\\ 1 & 1\\ 1 & 2 \end{vmatrix} = \begin{vmatrix} 3 & 3\\ 3 & 5 \end{vmatrix} ; (H^TH)^{-1}= \frac{1}{6} \begin{vmatrix} 5 & -3\\ -3 & 3 \end{vmatrix} $;

$H^TZ= \begin{vmatrix} 1 & 1& 1\\ 0 & 1& 2 \end{vmatrix} \begin{vmatrix} 1 \\ 2\\ 4 \end{vmatrix} = \begin{vmatrix} 7\\ 10 \end{vmatrix} ; X{оц}= \frac{1}{6} \begin{vmatrix} 5 & -3\\ -3 & 3 \end{vmatrix} \begin{vmatrix} 7\\ 10 \end{vmatrix} = \frac{1}{6} \begin{vmatrix} 5\\ 9 \end{vmatrix}$;

Итак $I_1 = 5/6 = 0,833 A$; $I_2 = 9/6 = 1,5 A$.

Задача фильтрации

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

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

Вначале получим оценку текущего значения $x_n$ по имеющимся $k$ последним значениям временного ряда $z_n, z_{n-1},z_{n-2}\dots z_{n-(k-1)}$. Модель наблюдения та же, что и в задаче оценки параметров:

Понятно, что $Z$ – это вектор–столбец, состоящий из наблюдаемых значений временного ряда $z_n, z_{n-1},z_{n-2}\dots z_{n-(k-1)}$, $V$ – вектор–столбец помехи $\xi _n, \xi _{n-1},\xi_{n-2}\dots \xi_{n-(k-1)}$, искажающий истинный сигнал. А что означают символы $H$ и $X$? О каком, например, векторе–столбце $X$ может идти речь, если все, что необходимо, – это дать оценку текущего значения временного ряда? А что понимать под матрицей преобразований $H$, вообще непонятно.

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

В данном случае оцениваются параметры именно этой модели. При выборе подходящей модели генерации сигнала вспомним о том, что любую аналитическую функцию можно разложить в ряд Тейлора. Поразительное свойство ряда Тейлора заключается в том, что форма функции на любом конечном расстоянии $t$ от некой точки $x=a$ однозначно определяется поведением функции в бесконечно малой окрестности точки $x=a$ (речь идет о ее производных первого и высшего порядков).

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

$x_{n-i} = F_{-i}x_n$, (4)

$$X_n= \begin{vmatrix} x_n\\ x"_n\\ x""_n \end{vmatrix} ; F_{-i}= \begin{vmatrix} 1 & -i & i^2/2\\ 0 & 1 & -i\\ 0 & 0 & 1 \end{vmatrix} $$

То есть формула 4, при заданном порядке полинома (в примере он равен 2) устанавливает связь между $n$-ым значением сигнала во временной последовательности и $(n-i)$–ым. Таким образом, оцениваемый вектор состояния в данном случае включает в себя, помимо собственно оцениваемого значения, первую и вторую производную сигнала.

В теории автоматического управления такой фильтр назвали бы фильтром с астатизмом 2-го порядка. Матрица преобразования $H$ для данного случая (оценка осуществляется по текущему и $k-1$ предшествующим выборкам) выглядит так:

$$H= \begin{vmatrix} 1 & -k & k^2/2\\ - & - & -\\ 1 & -2 & 2\\ 1 & -1 & 0.5\\ 1 & 0 & 0 \end{vmatrix}$$

Все эти числа получаются из ряда Тейлора в предположении, что временной интервал между соседними наблюдаемыми значениями постоянный и равен 1.

Итак, задача фильтрации при принятых нами предположениях свелась к задаче оценки параметров; в данном случае оцениваются параметры принятой нами модели генерации сигнала. И оценка значений вектора состояния $X$ осуществляется по той же формуле 3:

$$X_{оц}=(H^TH)^{-1}H^TZ$$

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

Формула 3 легко реализуется программно, для этого нужно заполнить матрицу $H$ и вектор столбец наблюдений $Z$. Такие фильтры называются фильтры с конечной памятью , так как для получения текущей оценки $X_{nоц}$ они используют последние $k$ наблюдений. На каждом новом такте наблюдения к текущей совокупности наблюдений прибавляется новое и отбрасывается старое. Такой процесс получения оценок получил название скользящего окна .

Фильтры с растущей памятью

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

Чтобы бороться с этим недостатком, необходимо перейти от фильтра с постоянной памятью к фильтру с растущей памятью . В таком фильтре число наблюдаемых значений, по которым производится оценка, должна совпадать с номером n текущего наблюдения. Это позволяет получать оценки, начиная с числа наблюдений, равного числу компонент оцениваемого вектора $X$. А это определяется порядком принятой модели, то есть сколько членов из ряда Тейлора используется в модели.

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

Дело в том, что к моменту n мы уже имеем оценку $X_{(n-1)оц}$, в которой содержится информация обо всех предыдущих наблюдениях $z_n, z_{n-1}, z_{n-2} \dots z_{n-(k-1)}$. Оценку $X_{nоц}$ получаем по очередному наблюдению $z_n$ с использованием информации, хранящейся в оценке $X_{(n-1)}{\mbox {оц}}$. Такая процедура получила название рекуррентной фильтрации и состоит в следующем:

  • по оценке $X_{(n-1)}{\mbox {оц}}$ прогнозируют оценку $X_n$ по формуле 4 при $i = 1$: $X_{\mbox {nоцаприори}} = F_1X_{(n-1)оц}$. Это априорная оценка;
  • по результатам текущего наблюдения $z_n$, эту априорную оценку превращают в истинную, то есть апостериорную;
  • эта процедура повторяется на каждом шаге, начиная с $r+1$, где $r$ – порядок фильтра.

Окончательная формула рекуррентной фильтрации выглядит так:

$X_{(n-1)оц} = X_{\mbox {nоцаприори}} + (H^T_nH_n)^{-1}h^T_0(z_n - h_0 X_{\mbox {nоцаприори}})$, (6)

где для нашего фильтра второго порядка:

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

При практической реализации этой формулы необходимо помнить, что входящая в него априорная оценка определяется по формуле 4, а величина $h_0 X_{\mbox {nоцаприори}}$ представляет собой первую компоненту вектора $X_{\mbox {nоцаприори}}$.

У фильтра с растущей памятью имеется одна важная особенность. Если посмотреть на формулу 6, то окончательная оценка есть сумма прогнозируемого вектора оценки и корректирующего члена. Эта поправка велика при малых $n$ и уменьшается при увеличении $n$, стремясь к нулю при $n \rightarrow \infty$. То есть с ростом n сглаживающие свойства фильтра растут и начинает доминировать модель, заложенная в нем. Но реальный сигнал может соответствовать модели лишь на отдельных участках, поэтому точность прогноза ухудшается.

Чтобы с этим бороться, начиная с некоторого $n$, накладывают запрет на дальнейшее уменьшение поправочного члена. Это эквивалентно изменению полосы фильтра, то есть при малых n фильтр более широкополосен (менее инерционен), при больших – он становится более инерционен.

Сравните рисунок 1 и рисунок 2. На первом рисунке фильтр имеет большую память, при этом он хорошо сглаживает, но в силу узкополосности оцениваемая траектория отстает от реальной. На втором рисунке память фильтра меньше, он хуже сглаживает, но лучше отслеживает реальную траекторию.

Литература

  1. Ю.М.Коршунов "Математические основы кибернетики"
  2. А.В.Балакришнан "Теория фильтрации Калмана"
  3. В.Н.Фомин "Рекуррентное оценивание и адаптивная фильтрация"
  4. К.Ф.Н.Коуэн, П.М. Грант "Адаптивные фильтры"

Фильтр Калмана - это, наверное, самый популярный алгоритм фильтрации, используемый во многих областях науки и техники. Благодаря своей простоте и эффективности его можно встретить в GPS-приемниках, обработчиках показаний датчиков, при реализации систем управления и т.д.

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

Для чего он нужен?

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

Фильтр - это алгоритм обработки данных, который убирает шумы и лишнюю информацию. В фильтре Калмана есть возможность задать априорную информацию о характере системе, связи переменных и на основании этого строить более точную оценку, но даже в простейшем случае (без ввода априорной информации) он дает отличные результаты.

Рассмотрим простейший пример - предположим нам необходимо контролировать уровень топлива в баке. Для этого в бак устанавливается емкостный датчик, он очень прост в обслуживании, но обладает некоторыми недостатками - например, зависимость от заправляемого топлива (диэлектрическая проницаемость топлива зависит от многих факторов, например, от температуры), большое влияние «болтанки» в баке. В итоге, информация с него представляет типичную «пилу» с приличной амплитудой. Такого рода датчики часто устанавливаются на тяжелой карьерной технике (не смущайтесь объемам бака):

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

Немного отвлечемся и познакомимся с самим алгоритмом. Фильтр Калмана использует динамическую модель системы (например, физический закон движения), известные управляющие воздействия и множество последовательных измерений для формирования оптимальной оценки состояния. Алгоритм состоит из двух повторяющихся фаз: предсказание и корректировка. На первом рассчитывается предсказание состояния в следующий момент времени (с учетом неточности их измерения). На втором, новая информация с датчика корректирует предсказанное значение (также с учетом неточности и зашумленности этой информации):

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

Разберемся сначала в обозначениях: подстрочный индекс обозначает момент времени: k - текущий, (k-1) - предыдущий, знак «минус» в верхнем индексе обозначает, что это предсказанное промежуточное значение.

Описание переменных представлены на следующих изображениях:

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

Опробуем в деле

Вернемся к примеру с датчиком уровня топлива, так как состояние системы представлено одной переменной (объем топлива в баке), то матрицы вырождаются в обычные уравнения:
Определение модели процесса
Для того, чтобы применить фильтр, необходимо определить матрицы/значения переменных определяющих динамику системы и измерений F, B и H:

F - переменная описывающая динамику системы, в случае с топливом - это может быть коэффициент определяющий расход топлива на холостых оборотах за время дискретизации (время между шагами алгоритма). Однако помимо расхода топлива, существуют ещё и заправки… поэтому для простоты примем эту переменную равную 1 (то есть мы указываем, что предсказываемое значение будет равно предыдущему состоянию).

B - переменная определяющая применение управляющего воздействия. Если бы у нас были дополнительная информация об оборотах двигателя или степени нажатия на педаль акселератора, то этот параметр бы определял как изменится расход топлива за время дискретизации. Так как управляющих воздействий в нашей модели нет (нет информации о них), то принимаем B = 0.

H - матрица определяющая отношение между измерениями и состоянием системы, пока без объяснений примем эту переменную также равную 1.

Определение сглаживающих свойств
R - ошибка измерения может быть определена испытанием измерительных приборов и определением погрешности их измерения.

Q - определение шума процесса является более сложной задачей, так как требуется определить дисперсию процесса, что не всегда возможно. В любом случае, можно подобрать этот параметр для обеспечения требуемого уровня фильтрации.

Реализуем в коде
Чтобы развеять оставшиеся непонятности реализуем упрощенный алгоритм на C# (без матриц и управляющего воздействия):

class KalmanFilterSimple1D
{
public double X0 {get; private set;} // predicted state
public double P0 { get; private set; } // predicted covariance

Public double F { get; private set; } // factor of real value to previous real value
public double Q { get; private set; } // measurement noise
public double H { get; private set; } // factor of measured value to real value
public double R { get; private set; } // environment noise

Public double State { get; private set; }
public double Covariance { get; private set; }

Public KalmanFilterSimple1D(double q, double r, double f = 1, double h = 1)
{
Q = q;
R = r;
F = f;
H = h;
}

Public void SetState(double state, double covariance)
{
State = state;
Covariance = covariance;
}

Public void Correct(double data)
{
//time update - prediction
X0 = F*State;
P0 = F*Covariance*F + Q;

//measurement update - correction
var K = H*P0/(H*P0*H + R);
State = X0 + K*(data - H*X0);
Covariance = (1 - K*H)*F;
}
}

// Применение...

Var fuelData = GetData();
var filtered = new List();

Var kalman = new KalmanFilterSimple1D(f: 1, h: 1, q: 2, r: 15); // задаем F, H, Q и R
kalman.SetState(fuelData, 0.1); // Задаем начальные значение State и Covariance
foreach(var d in fuelData)
{
kalman.Correct(d); // Применяем алгоритм

Filtered.Add(kalman.State); // Сохраняем текущее состояние
}

Результат фильтрации с данными параметрами представлен на рисунке (для настройки степени сглаживания - можно изменять параметры Q и R):

За рамками статьи осталось самое интересное - применение фильтра Калмана для нескольких переменных, задание взаимосвязи между ними и автоматический вывод значений для ненаблюдаемых переменных. Постараюсь продолжить тему как только появится время.

Надеюсь описание получилось не сильно утомительным и сложным, если остались вопросы и уточнения - добро пожаловать в комментарии)