Расширенный фильтр калмана пример. Фильтр Калмана — Введение

28.04.2019

Привет всем

Вот краткое руководство для реализации фильтра Калмана. Первоначально я написал это для статьи общества Robot несколько лет назад. Я пересмотрел это немного более четко и исправил некоторые ошибки в исходном сообщении.

Наслаждайтесь!

Введение

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

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

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

Фильтр

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

(У меня вопрос о том, почему я перечисляю положение и скорость. Ответ прост, если вы подумаете о квадрокоптере, его можно отметить в одном направлении во время полета / движения в другом направлении.)

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

УРАВНЕНИЯ ПРОГНОЗИРОВАНИЯ

УРАВНЕНИЯ ОБНОВЛЕНИЯ

Эти уравнения могут выглядеть страшно, так как есть много переменных, так что мы теперь уточним, что из себя представляет каждая переменная.
x , X = Эти переменные представляют свое состояние (то есть вещи, о которых Вы заботитесь и / или пытаетесь фильтровать). Например, если вы управляете механическим манипулятором с тремя шарнирами ваше состояние может быть:

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

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

u = какое бы ни было действие для трех соединений механической руки, действия могут быть:

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

где отклонение - variance - представляет собой квадрат стандартного отклонения ошибки.

Q = Является ковариацией шума процесса (т.е. действие). Он формируется по аналогии с выше за исключением того, что она является матрицей, которую вы определяете, и не обновляется с помощью фильтра. Эта матрица рассказывает фильтру Калмана, сколько ошибок в каждом действии с момента выдачи вами заданного напряжения, пока он на самом деле не происходит.
K = Эта матрица обновляется как часть измерения стадии обновления.
H = Является моделью датчиков, но ее трудно определить. Простой подход- инициализировать ее как диагональную матрицу идентичности и настроить, чтобы улучшить конечные результаты фильтрации.
R = Аналогично Q за исключением того, это определяет уверенность в каждом из датчиков. Это ключевая матрица для проведения синтеза датчика
z = Это измерения, которые возвращаются из датчиков в каждом месте. Фильтр Калмана обычно используется для очистки шума от этих сигналов или для оценки этих параметров, когда нет датчика.
I = единичная матрица (также диагональная)
Следующие переменные, которые мы должны определить, являются A и B. Эти матрицы - это модель, которая показывает, как ваша система переходит из одного состояния в другое.
Так как фильтр Калмана предназначен для линейных систем, мы можем предположить, что, если начать с самого начала (t = 0) и запустить вашу систему к следующему состоянию (t = 1), и запустите ее снова (t = 2), количество изменения будет таким же, так что мы не можем использовать это изменение, независимо от того, где (t = везде) в системе мы на самом деле находимся.

Самый простой способ найти матрицы А и В, если вы в состоянии собирать данные из вашей системы, это выполнить эксперименты перед созданием фильтра с системой без принуждения (входы должны быть 0). newState и lastState являются матрицы, чьи столбцы являются входные и выходные состояния, измеренные в ходе экспериментов. Для того, чтобы это работало, вы должны быть в состоянии измерить полное состояние, что часто не представляется возможным, если вы используете фильтр Калмана в качестве оценки состояния (который является частым в использовании). Это позволяет решить для А и В в цифровой форме. Вам нужно повторить эксперимент n-ное количество раз, где n есть размерность матрицы A.

тогда мы можем записать его в виде:

Следуя этой логике, как только у вас есть А вы можете найти B

Зная текущее состояние, какое последнее состояние было, A, и действие, которое вызвало изменение состояния, вы можете решить для B.

Примечание: В системах, где вы только фильтруете или дублируете входы (т.е. датчики), так что вы не можете быть способны иметь какие-либо действия, то термин B * действия отменяет и все, что вам нужно знать, это ваш термин А.

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

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

Таким образом, в псевдокоде это будет:

X 0 = состояние, которое вы линеаризуете;
U 0 = Напряжения, необходимые для производства (обратная динамика);

delta = малое число;
// Нахождение матрицы A

для ii =1 через # ofStates {
X = X 0;
X (ii )= X (ii )+ delta ; //беспокойное состояние x (i ) посредством delta
X 1 = f (X , U 0); // f () является моделью системы
A (:, ii ) = (X 1- X 0)/ delta ;
}

// Нахождение матрицы B

для ii =1 посредством # ofInputs {
U = U 0;
U (ii )= U (ii )+ delta ; // беспокойное состояние U (i ) посредством delta
X 1 = f (X 0, U ); // f () является моделью системы
B (:, ii ) = (X 1- X 0)/ delta ;
}

f(X0,U) это ваша модель. Эта модель может занять состояние и действие и определить, какое следующее состояние будет. Этот шаг включает в себя кинематику и динамику ровера. Помните, что если ваша модель не прогнозирует того, что происходит в реальной жизни, то ваша модель ошибается, и вы должны это исправить!

Теперь, когда мы знаем, что это все за переменные, я хочу повторить алгоритм со словами. Начнем с нашего текущего состояния x k -1 , текущей ковариационной матрицы P k -1 , а также текущего входа u k -1 чтобы получить ваше прогнозируемое состояние X и предсказанную ковариационную матрицу p. После этого вы можете получить измерения y k и использовать уравнения обновления, чтобы исправить свои предсказания для того, чтобы получить новое состояние матрицы x k и новый ковариационный P k . После того, как вы сделаете все это вы просто держите на итерации те шаги, где вы сделали прогноз, а затем обновите предсказание на основе некоторой известной информации.

Современные методы фильтров

Для нелинейной системы существует два основных подхода. Первый заключается в разработке расширенного фильтра Калмана (EKF - Extended Kalman Filter). Для EKF вам необходимо линеаризовать модель, а затем сформировать ваши матрицы А и В. Этот подход включает в себя немного математики и нечто, называемое Jacobean, который позволяет масштабировать различные значения по-разному. Второй и более простой подход состоит в использовании кусочного приближения. Для этого вы делаете разбивку данных в регионах, близких к линейным и образуете различные матрицы А и В для каждого региона. Это позволяет проверить данные и использовать соответствующие матрицы А и В в фильтре, чтобы точно предсказать, что переходное состояние будет.

Образец кода

Вот c++ код для фильтра Калмана предназначенный для PUMA 3 DOF роботизированного манипулятора. Этот код используется для оценки скорости, так как он является гораздо более точным, чем просто дифференциация позиции.

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

Я сделал этот код уже давно. И теперь немного смущен тем, как код выглядит, но им я все равно буду делиться.

/*******************************************************
* роботизированный манипулятор PUMA 3 DOF фильтра Калмана
********************************************/
# include < fstream >
# include < iostream >
# include < unistd . h >
# include < stdlib . h >
# include < string . h >
# include < stdio . h >
// моя матричная библиотека, вы можете использовать свою любимую матричную библиотеку
# include “ matrixClass . h ”
# define TIMESTEP 0.05

используя пространство имен std ;

индекс резерва времени, m _ a 1, m _ a 2, m _ a 3, c _ a 1, c _ a 2, c _ a 3, tau 1, tau 2, tau 3;
резерв a 1, a 2, a 3, a 1 d , a 2 d , a 3 d , a 1 dd , a 2 dd , a 3 dd ;
резерв new _ a 1 d , new _ a 2 d , new _ a 3 d ;
резерв TIME =0;
матрица A (6,6);
матрица B (6,3);
матрица C = matrix :: eye (6); // инициализировать их как 6×6 матриц идентичности
матрица Q = matrix :: eye (6);
матрица R = matrix :: eye (6);
матрица y (6,1);
матрица K (6,6);
матрица x (6,1);
матричное состояние(6,1);
матричное действие(3,1);
матрица lastState (6,1);
матрица P = matrix :: eye (6);
матрица p = matrix :: eye (6);
матричный замер(6,1);

пустой initKalman (){

резерв a ={
{1.004, 0.0001, 0.001, 0.0014, 0.0000, -0.0003 },
{0.000, 1.000, -0.00, 0.0000, 0.0019, 0 },
{0.0004, 0.0002, 1.002, 0.0003, 0.0001, 0.0015 },
{0.2028, 0.0481, 0.0433, 0.7114, -0.0166, -0.1458 },
{0.0080, 0.0021, -0.0020, -0.0224, 0.9289, 0.005 },
{0.1895, 0.1009, 0.1101, -0.1602, 0.0621, 0.7404 }
};

резерв b = {
{0.0000, 0.0000 , 0.0000 },
{0.0000, 0.0000, -0.0000 },
{0.0000, -0.0000, 0.0000 },
{0.007, -0.0000, 0.0005 },
{0.0001, 0.0000, -0.0000 },
{0.0003, -0.0000, 0.0008 }
};

/* загружает матрицы А и В из выше */
для (int i = 0; i < 6; i ++){
для (int j = 0; j < 6; j ++){
A [ i ][ j ]= a [ i ][ j ];
}
}
для (int i = 0; i < 6; i ++){
для (int j = 0; j < 3; j ++){
B [ i ][ j ]= b [ i ][ j ];
}
}

/* инициализирует состояние*/
state=0.1;
state=0.1;
state=0.1;
state=0.1;
state=0.1;
state=0.1;

lastState= состояние ;

пустой kalman(){
lastState= состояние ;
state=c_a1;
state=c_a2;
state=c_a3;
state=a1d;
state=a2d;
state=a3d;

measurement=m_a1;
measurement=m_a2;
measurement=m_a3;
measurement=a1d;
measurement=a2d;
measurement=a3d;

action=tau1;
action=tau2;
action=tau3;

matrix temp1(6,6);
matrix temp2(6,6);
matrix temp3(6,6);
matrix temp4(6,1);
/************ Уравнения п рогнозировани я *****************/
x = A*lastState + B*action;
p = A*P*A’ + Q;
/************ Уравнения обновления **********/
K = p*C*pinv(C*p*C’+R);

y=C*state;

состояние = x + K*(y-C*lastState);

P = (eye(6) – K*C)*p;

a1=state;
a2=state;
a3=state;
a1d=state;
a2d=state;
a3d=state;
}

/* Эта функция не используется, так как я использую позиции, чтобы получить скорости (т.е. дифференциацию).
* Тем не менее я считаю, что это полезно включить, если вы хотите скорости и позицию
* от ускорения вы будете использовать его */
пустой integrate(){
new_a1d = a1d + a1dd*TIMESTEP;
a1 += (new_a1d + a1d)*TIMESTEP/2;
a1d = new_a1d;
new_a2d = a2d + a2dd*TIMESTEP;
a2 += (new_a2d + a2d)*TIMESTEP/2;
a2d = new_a2d;
new_a3d = a3d + a3dd*TIMESTEP;
a3 += (new_a3d + a3d)*TIMESTEP/2;
a3d = new_a3d;
TIME+=TIMESTEP;
}

/*Это дает мне скорость от позиции*/
пустой differentiation(){
a1d=(state-lastState)/TIMESTEP;
a2d=(state-lastState)/TIMESTEP;
a3d=(state-lastState)/TIMESTEP;
TIME+=TIMESTEP;
}

int main () {
initKalman();
char buffer;
ifstream readFile (“DATA.txt”); // вот тут я читал мои данные, так как я обрабатывал все это оффлайн

while (!readFile.eof()){
readFile.getline (buffer,500);
sscanf(buffer, “%f %f %f %f %f %f %f %f %f %f “,&timeIndex,&m_a1,&m_a2,&m_a3,&c_a1,&c_a2,&c_a3,&tau1,&tau2,&tau3);

kalman();
differentiation();
//integrate();

/* вот тут я вношу в журнал результаты и / или выношу их на экран */
FILE *file=fopen(“filterOutput.txt”, “a”);
fprintf(file,”%f %f %f %f %f %f %f %f %f %f \n”,TIME,a1,a2,a3,a1d,a2d,a3d,tau1,tau2,tau3);
fprintf(stderr,”%f %f %f %f %f %f %f %f %f %f \n”,TIME,a1,a2,a3,a1d,a2d,a3d,tau1,tau2,tau3);
fclose(file);
}
return 1;
}

Заключение

Этот пост был сосредоточен на реализации фильтра Калмана. Надеюсь, вы сможете принять эту информацию, чтобы улучшить и усовершенствовать ваши робототехнические проекты. Для получения дополнительной информации я рекомендую Введение в фильтрацию Калмана Грега Уэлча и Гэри Бишопа http://www.cs.unc.edu/~welch/kalman/ и пост из TKJ Electronics .

Фильтр Калмана - это, наверное, самый популярный алгоритм фильтрации, используемый во многих областях науки и техники. Благодаря своей простоте и эффективности его можно встретить в 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)*P0; } } // Применение... 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):

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

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

UPD: Список источников.

Tutorial

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

Ликбез

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

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

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

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

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

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

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

Дисперсия

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

Еще раз полюбоваться распределением Гаусса

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

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

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

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

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

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

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

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

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

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

ключ к доказательству

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

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

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

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

Все формулы в одном месте

Пример

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

Еще раз посмотреть на результат фильтрования

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

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; % eOpt(t) is a square root of the error dispersion (variance). It"s not a random variable. 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)

Анализ

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

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

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

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

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

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

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

Обсуждение

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

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

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

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

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

Ликбез

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

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

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

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

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

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

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

Дисперсия

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

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

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

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

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

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

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

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

(1)

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

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

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

Коэффициент называют коэффициентом Калмана. Он зависит от шага итерации, поэтому правильнее было бы писать , но пока, чтобы не загромождать формулы расчетах, мы будем опускать его индекс.
Мы должны выбрать коэффициент Калмана таким, чтобы получившееся оптимальное значение координаты было бы наиболее близко к истинной . К примеру, если мы знаем, что наш сенсор очень точный, то мы будем больше доверять его показанию и дадим значению больше весу ( близко единице). 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)

Анализ

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

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

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

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

Обсуждение

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

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

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

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

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

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

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

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

Введение

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

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

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

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

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

Расширенный фильтр Калмана (EKF, Extended Kalman filter). Сведение нелинейных моделей наблюдений и формирующего процесса с помощью линеаризации посредством разложения в ряд Тейлора .
Unscented Kalman filter (UKF). Используется в задачах, в которых простая линеаризация приводит к уничтожению полезных связей между компонентами вектора состояния. В этом случае «линеаризация» основана на unscented -преобразовании.
Ensemble Kalman filter (EnKF). Используется для уменьшения размерности задачи.
Возможны варианты с нелинейным дополнительным фильтром, позволяющим привести негауссовские наблюдения к нормальным.
Возможны варианты с «обеляющим» фильтром, позволяющим работать с «цветными» шумами
и т. д.

Используемая модель динамической системы

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

При использовании фильтра Калмана для получения оценок вектора состояния процесса по серии зашумленных измерений необходимо представить модель данного процесса в соответствии со структурой фильтра - в виде матричного уравнения определенного типа. Для каждого такта k работы фильтра необходимо в соответствии с приведенным ниже описанием определить матрицы: эволюции процесса F k ; матрицу наблюдений H k ; ковариационную матрицу процесса Q k ; ковариационную матрицу шума измерений R k ; при наличии управляющих воздействий - матрицу их коэффициентов B k .

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

Модель системы/процесса подразумевает, что истинное состояние в момент k получается из истинного состояния в момент k −1 в соответствии с уравнением:

F k - матрица эволюции процесса/системы, которая воздействует на вектор x k −1 (вектор состояния в момент k −1 );
B k - матрица управления, которая прикладывается к вектору управляющих воздействий u k ;
w k - нормальный случайный процесс с нулевым математическим ожиданием и ковариационной матрицей Q k , который описывает случайный характер эволюции системы/процесса:

В момент k производится наблюдение (измерение) z k истинного вектора состояния x k , которые связаны между собой уравнением:

где H k - матрица измерений, связывающая истинный вектор состояния и вектор произведенных измерений, v k - белый гауссовский шум измерений с нулевым математическим ожиданием и ковариационной матрицей R k :

Начальное состояние и векторы случайных процессов на каждом такте {x 0 , w 1 , …, w k , v 1 , …, v k } считаются независимыми .

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

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

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

Состояние фильтра задается двумя переменными:

Итерации фильтра Калмана делятся на две фазы: экстраполяция и коррекция. Во время экстраполяции фильтр получает предварительную оценку состояния системы (в русскоязычной литературе часто обозначается , где означает «экстраполяция», а k - номер такта, на котором она получена) на текущий шаг по итоговой оценке состояния с предыдущего шага (либо предварительную оценку на следующий такт по итоговой оценке текущего шага, в зависимости от интерпретации). Эту предварительную оценку также называют априорной оценкой состояния, так как для её получения не используются наблюдения соответствующего шага. В фазе коррекции априорная экстраполяция дополняется соответствующими текущими измерениями для коррекции оценки. Скорректированная оценка также называется апостериорной оценкой состояния, либо просто оценкой вектора состояния . Обычно эти две фазы чередуются: экстраполяция производится по результатам коррекции до следующего наблюдения, а коррекция производится совместно с доступными на следующем шаге наблюдениями, и т. д. Однако возможно и другое развитие событий, если по некоторой причине наблюдение оказалось недоступным, то этап коррекции может быть пропущен и выполнена экстраполяция по нескорректированной оценке (априорной экстраполяции). Аналогично, если независимые измерения доступны только в отдельные такты работы, всё равно возможны коррекции (обычно с использованием другой матрицы наблюдений H k ).

Этап экстраполяции

Этап коррекции

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

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

Инварианты

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

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

где - математическое ожидание .

Расчетные матрицы ковариаций экстраполяций, оценок состояния системы и вектора ошибок совпадают с истинными матрицами ковариаций:

Пример построения фильтра

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

В данной задаче матрицы F , H , R и Q не зависят от времени, опустим их индексы. Кроме того, мы не управляем вагонеткой, поэтому матрица управления B отсутствует.

Координата и скорость вагонетки описывается вектором в линейном пространстве состояний

где - скорость (первая производная координаты по времени).

Будем считать, что между (k −1 )-ым и k -ым тактами вагонетка движется с постоянным ускорением a k , распределенным по нормальному закону с нулевым математическим ожиданием и среднеквадратическим отклонением σ a . В соответствии с механикой Ньютона можно записать

Ковариационная матрица случайных воздействий

(σ a - скаляр).

На каждом такте работы производится измерение положения вагонетки. Предположим, что погрешность измерений v k имеет нормальное распределение с нулевым математическим ожиданием и среднеквадратическим отклонением σ z . Тогда

и ковариационная матрица шума наблюдений имеет вид

Начальное положение вагонетки известно точно

, .

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

, .

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

Вывод формул

Ковариационная матрица оценки вектора состояния

По определению ковариационной матрицы P k |k

подставляем выражение для оценки вектора состояния

и расписываем выражение для вектора ошибок

и вектора измерений

выносим вектор погрешности измерений v k

так как вектор погрешности измерений v k не коррелирован с другими аргументами, получаем выражение

в соответствии со свойствами ковариации векторов данное выражение преобразуется к виду

заменяя выражение для ковариационной матрицы экстраполяции вектора состояния на P k |k −1 и определение ковариационной матрицы шумов наблюдений на R k , получаем

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

Оптимальная матрица коэффициентов усиления

Фильтр Калмана минимизирует сумму квадратов математических ожиданий ошибок оценки вектора состояния.

Вектор ошибки оценки вектора состояния

Стоит задача минимизировать сумму математических ожиданий квадратов компонент данного вектора

что эквивалентно минимизации следа ковариационной матрицы оценки вектора состояния P k |k . Подставим в выражение для ковариационной матрицы оценки вектора состояния имеющиеся выражения и дополним до полного квадрата:

Заметим что, последнее слагаемое является ковариационной матрицей некоторой случайной величины, поэтому его след неотрицателен. Минимум следа достигнется при обнулении последнего слагаемого:

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

Ковариационная матрица оценки вектора состояния при использовании оптимальной матрицы коэффициентов

Выражение для ковариационной матрицы оценки вектора состояния P k |k при использовании оптимальной матрицы коэффициентов примет вид:

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