3

我目前正在尝试使用Free C++ Extended Kalman Filter Library。我了解卡尔曼滤波器的基础知识，但是我遇到了使用此库生成 NaN 值的问题。SO上有没有人有使用卡尔曼滤波器算法来发现我的错误的经验？

这是我的过滤器：

class PointEKF : public Kalman::EKFilter<double,1,false,true,false> {
public:
        PointEKF() : Period(0.0) {
            setDim(3, 1, 3, 1, 1);
        }

        void SetPeriod(double p) {
            Period = p;
        }
protected:
        void makeBaseA() {
            A(1, 1) = 1.0;
            //A(1, 2) = Period;
            //A(1, 3) = Period*Period / 2;
            A(2, 1) = 0.0;
            A(2, 2) = 1.0;
            //A(2, 3) = Period;
            A(3, 1) = 0.0;
            A(3, 2) = 0.0;
            A(3, 3) = 1.0;
        }
        void makeBaseH() {
            H(1, 1) = 1.0;
            H(1, 2) = 0.0;
            H(1, 3) = 0.0;
        }
        void makeBaseV() { 
            V(1, 1) = 1.0;
        }
        void makeBaseW() {
            W(1, 1) = 1.0;
            W(1, 2) = 0.0;
            W(1, 3) = 0.0;
            W(2, 1) = 0.0;
            W(2, 2) = 1.0;
            W(2, 3) = 0.0;
            W(3, 1) = 0.0;
            W(3, 2) = 0.0;
            W(3, 3) = 1.0;
        }

        void makeA() {
            double T = Period;
            A(1, 1) = 1.0;
            A(1, 2) = T;
            A(1, 3) = (T*T) / 2;
            A(2, 1) = 0.0;
            A(2, 2) = 1.0;
            A(3, 3) = T;
            A(3, 1) = 0.0;
            A(3, 2) = 0.0;
            A(3, 3) = 1.0;
        }
        void makeH() {
            double T = Period;
            H(1, 1) = 1.0;
            H(1, 2) = T;
            H(1, 3) = T*T / 2;
        }
        void makeProcess() {
            double T = u(1);
            Vector x_(x.size());
            x_(1) = x(1) + x(2) * T + (x(3) * T*T / 2);
            x_(2) = x(2) + x(3) * T;
            x_(3) = x(3);
            x.swap(x_);
        }
        void makeMeasure() {
            z(1) = x(1);
        }

        double Period;
};

我使用它如下：

void init() {
    int n = 3;
    static const double _P0[] = {
                                 1.0, 0.0, 0.0,
                                 0.0, 1.0, 0.0,
                                 0.0, 0.0, 1.0
                                };
    Matrix P0(n, n, _P0);
    Vector x(3);
    x(1) = getPoint(0);
    x(2) = getVelocity(0);
    x(3) = getAccleration(0);
    filterX.init(x, P0);
}

和，

    Vector measurement(1), input(1), u(1);
    u(1) = 0.400;
    double start = data2->positionTimeCounter;
    double end = data->positionTimeCounter;
    double period = (end - start) / (1000*1000);
    filterX.SetPeriod(period);
    measurement(1) = getPoint(0);
    input(1) = period;
    filterX.step(input, measurement);
    auto x = filterX.predict(u);

注意： 我使用的数据是从单位圆生成的 x 点。

4

1 回答 1

4

如果您使用矩阵的基本版本：

A = [ 1 0 0;
      0 1 0;
      0 0 1 ];
H = [ 1 0 0 ];

您没有可观察的系统，因为您的测量仅捕获第一个状态（位置），并且在 A 矩阵中，位置及其导数（速度、加速度）之间没有耦合。可观察性矩阵如下：

O = [ H;
      H*A;
      H*A*A ];
O = [ 1 0 0;
      1 0 0;
      1 0 0 ];

这显然是奇异的，即您的系统是不可观察的。并且通过 EKF 算法输入应该会产生错误（这种情况应该由算法检测到），但如果没有检测到，它将导致估计中的 NaN 结果，与您所经历的完全一样。

现在，makeA() 函数中的 A 矩阵更合适：

A = [ 1 h h*h/2;
      0 1 h;
      0 0 1 ];
H = [ 1 0 0 ];       // use this H matrix (not [ 1 h h*h/2 ])

导致可观察性矩阵：

O = [ 1    0      0;
      1    h  h*h/2;
      1  2*h  2*h*h ];

这是全秩（不是单数），因此，你有一个可观察的系统。

卡尔曼滤波算法可能对矩阵的条件非常敏感，这意味着如果时间步长非常小（例如1e-6），则需要使用连续时间版本。此外，NaN 的问题可能来自 KF 算法所需的线性求解器（求解线性方程组）。如果库的作者使用了一种简单的方法（例如，高斯消元法、带或不带枢轴的 LU 分解、没有枢轴的 Cholesky 等），那么这将使数值调节问题变得更糟。

注意你应该从一个非常高的 P 矩阵开始你的 KF 过滤，因为初始 P 应该反映你的初始状态向量的不确定性，这通常非常高，所以 P 应该在1000 * identity.

于 2013-01-10T16:47:14.273 回答