python - 实现一维卡尔曼滤波器/平滑 Python

Question

我想测试卡尔曼滤波器以平滑我拥有的一组数据。请注意，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-');

每个点都是一个观察值。很明显，x=50 处的点是噪声。因此我希望卡尔曼滤波器的结果是这样的：

我不是数学专家，所以我不确定它是否重要，但我的数据不是速度或位置（我发现的所有卡尔曼示例都指的是那个案例）。问题是我不知道如何将这个相当简单的问题实现到 Python 中的卡尔曼滤波器。我看到很多人使用这个pykalman包

我的第一个问题是——卡尔曼滤波器可以处理不相等的时间间隔吗？如果答案是否定的，那么假设我的数据中的时间间隔相等，我仍然希望得到答案。我还在示例中看到数据应该是一种特定的方式，而不是像我的示例中那样“简单”的两个列表。所以我的第二个问题是，如何在 Python 中应用卡尔曼滤波器/平滑，盯着我的“简单”两个列表（如果出现问题，您可以将 x 间隔更改为相等）。

标签： pythonkalman-filterpykalman

Answer 1

卡尔曼滤波器可以处理不相等的时间间隔吗？

是的。您需要注意两件事 - 间隔之间的时间步长不同，您需要考虑这将对转换矩阵（描述系统动力学 - 这些通常具有 delta-t 依赖性）和协方差矩阵产生的影响 -特别是过渡协方差（观察之间的时间越长，系统如何演变的不确定性就越大。

我不确定这是否重要，但我的数据不是速度或位置（我发现的所有卡尔曼示例均指该案例）

您可以根据需要应用卡尔曼滤波器。但是请记住，卡尔曼滤波器实际上是一个状态估计器。特别是它是具有线性动力学和高斯噪声的系统的最佳状态估计器。术语“过滤器”可能有点误导。如果您没有要表示其动态的系统，则需要“弥补”一些动态以捕捉您对生成数据的物理过程的直觉/理解。

很明显，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()

这会产生以下情节：