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

во- первых: это не фильтр, это оцениватель-корректор

во-вторых: это рекурсивный метод избавления от шума и инструментальных ошибок

предположим есть некая величина x изменяющаяся во времени и у вас есть некий прибор, измеряющий величину x с некоторой погрешностью и выдающий результат измерения y. вы хотите оценить, чему НА САМОМ ДЕЛЕ равна величина x - "отфильтровать погрешность измерений и шум среды"

алгоритм Калмана позволяет оценить истинное значение x_t на основе ИСТОРИИ ИЗМЕНЕНИИЙ измеряемой величины y. при этом делается ряд предположений:

предположение #1

величина x в момент времени t представляет собой линейную комбинацию величины x, управляющего сигнала u и шума w в прошлый момент времени t-1:

    xt = A * xt-1 +  B * ut-1 + Q * wt-1

где A - матрица системы , B - матрица управления , Q - матрица дисперсий шума

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

шум должен иметь нормальное распределение. самое важное - дисперсия шума dw должна быть ОДИНАКОВА ПО ВСЕЙ ВЫБОРКЕ (это не так в финансах - кластеризация волатильности - хорошо известный факт)

предположение #2

измеряемая величина y - это линейная комбинация реальной величины x (с учетом параметров измерителя H) и погрешности измерения v (с учетом R - матрицы дисперсий измерений):

    yt = H * xt  +  R * vt

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

очень редко на практике H = I (так как разные орты значения могут иметь различные инструментальные погрешности - именно так ВСЕГДА обстоят дела в финансах)

что мы только что сделали? мы построили модель, описывающую нашу систему

что мы хотим сделать дальше? мы хотим получить некое уравнение:

      xt = f (xt-1 , xt-2 , .. , x0 , yt , yt-1 , .. , t0 , dw , dv)

где dw и dv - дисперсии шума и погрешностей измерения, позволяющее нам оценивать истинное значение x исходя из исторических оценок и текущих измеряемых величин

• A — это известная величина, которая контролирует эволюцию системы. мы ее знаем из построенной нами модели
• w - ошибка модели и v - ошибка сенсоров — это случайные величины. их законы распределения не зависят от времени (от номера итерации). СРЕДНИЕ ЗНАЧЕНИЯ ОШИБОК РАВНЫ НУЛЮ. закон распределения случайных величин МОЖЕТ БЫТЬ нам и НЕ ИЗВЕСТЕН, но известны их дисперсии Q и R

будем рассуждать по индукции

представьте себе, что на t-ом шаге мы уже нашли отфильтрованное значение x'(t+1), которое хорошо ПРИБЛИЖАЕТ истинную величину. поэтому, еще не получая ФАКТИЧЕСКОЕ значение y(t+1) с сенсоров, мы можем предположить, что на шаге t+1 система эволюционирует согласно этому закону и сенсоры покажут что-то близкое к x'(t+1)

с другой стороны, на шаге t+1 у нас на руках будет НЕТОЧНОЕ показание сенсоров y(t+1)

итак, что же получается? на шаге (t+1) мы имеем

• во-первых, неточное показание сенсоров y(t+1)
• во-вторых, неточно РАССЧИТАННОЕ значение x'(t+1)

идея Калмана состоит в том, что чтобы получить наилучшее приближение к истинному значению - мы должны выбрать золотую середину между НЕТОЧНЫМ показанием y(t+1) сенсоров и нашим, тоже НЕТОЧНЫМ, предсказанием x'(t+1) того, что мы ожидали от них получить

показанию сенсоров мы дадим вес K, меняющийся где-то в районе от 0 до 1 - то больше, то меньше

сам алгоритм (псевдокод):

  шаг 0 - иницализация:
          имеем матрицы A , B , H
          имеем матрицы дисперсий случайных величин Q , R
          задаем н.у. - параметры x₀ , u₀ , y₀ , K₀ = 0 , P₀ != 0
    
  шаг 1 - считаем априорное состояние из параметров
             x' = A * x   +   B * u  

  шаг 2 - считаем априорное с датчиков
             y' = H * x'
  
  шаг 3 - проецируем из предиктора
             P' = A  *  P   *  A⁺   +   Q
             P" = H  *  P'  *  H⁺   +   R
  
  шаг 4 - вычисляем корректор по проекциям:
             K =  P'  *  H⁺  *  P"

  шаг 5 - вычисляем предиктор для следующей итерации:
             P = (I - K * H) * P'

  шаг 6 - берем показания y с датчиков
          и вычисляем отфильтованное от шума и инструментальных ошибок
          состояние:
              x = x' +  K * (y - y') 
          это - (вместе с шагом 5) результат очередной итерации
  
  переходим к шагу 1, передав полученные значения предиктора P и величины x
  (а также новые u и y извне) как параметры для следующей итерации 

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

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

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

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

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

а теперь немного теории из статьи Калмана за 1960 год:

---

множество 𝕐(t) = y(t), .. , y(0) образует линейное многообразие (векторное пространство) и ЗНАЯ ЕГО, мы должны оценить x с минимальной ошибкой. если время лежит в интервале [0,t), то это - интерполяционная задача. если же вермя равно t - то это задача фильтрации. если время >t, то это - задача экстраполяции. и все это вместе называется "оценкой"

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

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

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

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

линейная динамическая система может быть представлена дифференциальным уравнением:

    dx/dt = F(t) * x + D(t) * u
    
    y(t)  = M(t) * x 

где x - это вектор "состояния" системы, u - вектор "входов" системы. вектор y - "выходы" системы

если все коэффициенты матриц F,D,M - константы, то такая динамическая система называется ЛИВ-системой (Линейной Инвариантой во Времени)

                   ∞
    Ф(m) = exp F = Σ Fi/m! 
                  i=0
  

где F⁰ = I

если система не стационарна, то мы не можем выразить x точно, но лишь как:

    x(t+1) = Ф(t+1, t) * x(t) + Δ(t) * u(t)  

матрица Ф называется матрицей перехода системы из состояния t в состояние (t+1). естественно, что Ф(t,t) = I

в ЛИВ-системе

    Ф(t,τ) = exp F(t-τ) 

........... TODO ..............

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

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

вывод программы:

λ> main
input  0.39	result  0.39
input  0.50	result  0.39
input  0.48	result  0.40
input  0.29	result  0.40
input  0.25	result  0.39
input  0.32	result  0.38
input  0.34	result  0.37
input  0.48	result  0.37
input  0.41	result  0.38
input  0.45	result  0.39

а выведя все в файл "kalman.dat" можем и график построить:

gnuplot> plot "kalman.dat" using 2 with lines, "kalman.dat" using 4 with lines