在网格化数据上高效实现卡尔曼滤波器

Question

我写了一个非常简单的卡尔曼滤波器，它对时间序列进行操作（包括数据间隙）。它工作得很好，但我碰巧有一个数据立方体（例如形状Nt, Ny, Nx数组），我想对数据立方体中的每个像素应用时间卡尔曼滤波器。我已经完成了显而易见的（在最后两个维度上循环），但这需要相当长的时间。

最终，我总是不得不从单个"像素"中提取数据，并构建相关的矩阵/向量，因此该过程非常缓慢（请注意，每个单独时间序列中的间隙是不同的，通常，将状态与观察结果联系起来的 H 矩阵也是如此）。我不熟悉cython，它可能会有所帮助（只是我不熟悉它）。

我只是想知道对问题的巧妙改写或巧妙的数据结构是否可以更有效地进行时间过滤。我宁愿这只使用 numpy/scipy，而不是 OpenCV，否则取决于额外的包会很麻烦。

Answer 1

我做了一个简单的矢量化卡尔曼滤波器，像这样处理电影帧。它非常快，但目前仅限于 1D 输入和输出，并且它不对任何滤波器参数进行 EM 优化。

import numpy as np
def runkalman(y, RQratio=10., meanwindow=10):
    """
    A simple vectorised 1D Kalman smoother
    y       Input array. Smoothing is always applied over the first 
            dimension
    RQratio     An estimate of the ratio of the variance of the output 
            to the variance of the state. A higher RQ ratio will 
            result in more smoothing.
    meanwindow  The initial mean and variance of the output are 
            estimated over this number of timepoints from the start 
            of the array
    References:
        Ghahramani, Z., & Hinton, G. E. (1996). Parameter Estimation for
        Linear Dynamical Systems.
        http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.55.5997
        Yu, B., Shenoy, K., & Sahani, M. (2004). Derivation of Extented
        Kalman Filtering and Smoothing Equations.
        http://www-npl.stanford.edu/~byronyu/papers/derive_eks.pdf
    """
    x,vpre,vpost = forwardfilter(y, RQratio=RQratio, meanwindow=meanwindow)
    x,v = backpass(x, vpre, vpost)
    x[np.isnan(x)] = 0.
    return x
def forwardfilter(y, RQratio=10., meanwindow=10):
    """
    the Kalman forward (filter) pass:
        xpost,Vpre,Vpost = forwardfilter(y)
    """
    y = np.array(y,copy=False,subok=True,dtype=np.float32)
    xpre = np.empty_like(y)
    xpost = np.empty_like(y)
    Vpre = np.empty_like(y)
    Vpost = np.empty_like(y)
    K = np.empty_like(y)
    # initial conditions
    pi0 = y[:meanwindow].mean(0)
    ystd = np.std(y,0)
    R = ystd * ystd
    Q = R / RQratio
    V0 = Q
    xpre[0] = xpost[0] = pi0
    Vpre[0] = Vpost[0] = V0
    K[0] = 0
    # loop forwards through time
    for tt in xrange(1, y.shape[0]):
        xpre[tt] = xpost[tt-1]
        Vpre[tt] = Vpost[tt-1] + Q
        K[tt] = Vpre[tt] / (Vpre[tt] + R)
        xpost[tt] = xpre[tt] + K[tt] * (y[tt] - xpre[tt])
        Vpost[tt] = Vpre[tt] - K[tt] * (Vpre[tt])
    return xpost,Vpre,Vpost
def backpass(x, Vpre, V):
    """ 
    the Kalman backward (smoothing) pass:
        xpost,Vpost = backpass(x,Vpre,V)
    """
    xpost = np.empty_like(x)
    Vpost = np.empty_like(x)
    J = np.empty_like(x)
    xpost[-1] = x[-1]
    Vpost[-1] = V[-1]
    # loop backwards through time
    for tt in xrange(x.shape[0]-1, 0, -1):
        J[tt-1] = V[tt-1] / Vpre[tt]
        xpost[tt-1] = x[tt-1] + J[tt-1] * (xpost[tt] - x[tt-1])
        Vpost[tt-1] = V[tt-1] + J[tt-1] * (Vpost[tt] - Vpre[tt]) * J[tt-1]
    return xpost,Vpost 

如果有人知道 Python 中支持多维输入/输出和 EM 参数优化的矢量化卡尔曼更平滑的实现，我很想听听！我向 PyKalman 维护者提出了一个功能请求，但他们说清晰度比速度更重要。