Реализация 1D фильтра Калмана / сглаживания Python - PullRequest
Новая
       88

Реализация 1D фильтра Калмана / сглаживания Python

0 голосов
/ 08 мая 2020

Я хотел бы протестировать фильтр Калмнана, чтобы сгладить набор данных, которые у меня есть. Обратите внимание, что интервалы по оси x не равны. 

x = [1,10,22,35,40,51,59,72,85,90,100]
y = [0.2,0.23,0.3,0.4,0.5,0.2,0.65,0.67,0.62,0.5,0.4]
plt.plot(x,y, 'go-');

enter image description here

Где каждая точка является наблюдением. Очевидно, что точка при x = 50 - это шум. Следовательно, я ожидаю, что результаты фильтра Калмана будут примерно такими:

enter image description here

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

Мой первый вопрос: может ли фильтр Калмана обрабатывать интервалы времени, которые не равны? Если ответ отрицательный, то я все равно хотел бы получить ответ при условии, что интервалы времени в моих данных равны. Я также видел в примерах, что данные должны быть специфическим c способом, а не такими "простыми" двумя списками, как в моем примере. Итак, мой второй вопрос: как я могу применить фильтр Калмана / сглаживание в Python, глядя на мои «простые» два списка (вы можете изменить интервалы x, чтобы они были равны, если это проблема).

1 Ответ

0 голосов
/ 22 мая 2020

может ли фильтр Калмана обрабатывать временные интервалы, которые не равны?

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

Я не уверен, имеет ли это значение, но мои данные - это не скорость или местоположение (все примеры Калмана, которые я нашел, относятся к этому случаю)

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

Очевидно, что точка x = 50 - это шум.

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

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

from pykalman import KalmanFilter
import numpy as np
import matplotlib.pyplot as plt
import copy

outlier_thresh = 0.95

# Treat y as position, and that y-dot is
# an unobserved state - the velocity,
# which is modelled as changing slowly (inertia)

# state vector [y,
#               y_dot]

# transition_matrix =  [[1, dt],
#                       [0, 1]]

observation_matrix = np.asarray([[1, 0]])

# observations:
t = [1,10,22,35,40,51,59,72,85,90,100]

# dt betweeen observations:
dt = [np.mean(np.diff(t))] + list(np.diff(t))
transition_matrices = np.asarray([[[1, each_dt],[0, 1]]
                                    for each_dt in dt])

# observations
y = np.transpose(np.asarray([[0.2,0.23,0.3,0.4,0.5,0.2,
                              0.65,0.67,0.62,0.5,0.4]]))

y = np.ma.array(y)


leave_1_out_cov = []

for i in range(len(y)):
    y_masked = np.ma.array(copy.deepcopy(y))
    y_masked[i] = np.ma.masked

    kf1 = KalmanFilter(transition_matrices = transition_matrices,
                   observation_matrices = observation_matrix)

    kf1 = kf1.em(y_masked)

    leave_1_out_cov.append(kf1.observation_covariance[0,0])

# Find indexes that contributed excessively to observation covariance
outliers = (leave_1_out_cov / np.mean(leave_1_out_cov)) < outlier_thresh

for i in range(len(outliers)):
    if outliers[i]:
        y[i] = np.ma.masked


kf1 = KalmanFilter(transition_matrices = transition_matrices,
                   observation_matrices = observation_matrix)

kf1 = kf1.em(y)

(smoothed_state_means, smoothed_state_covariances) = kf1.smooth(y)


plt.figure()
plt.plot(t, y, 'go-', label="Observations")
plt.plot(t, smoothed_state_means[:,0], 'b--', label="Value Estimate" )
plt.legend(loc="upper left")
plt.xlabel("Time (s)")
plt.ylabel("Value (unit)")

plt.show()

В результате получается следующий график:

enter image description here

...