每日算法 — 使用java实现卡尔曼滤波:最优状态估计与预测更新

卡尔曼滤波是一种高效的递归滤波器,能够从一系列包含噪声的测量中估计动态系统的状态。它广泛应用于导航定位、目标跟踪、自动驾驶等领域。本文用Java完整实现卡尔曼滤波算法,讲解预测-更新双阶段模型、协方差矩阵传播与最优增益计算,帮助读者掌握这一经典状态估计算法。

一、算法原理

卡尔曼滤波的核心思想是融合两种信息来源:

  1. 系统模型预测:根据上一时刻的状态和控制输入,预测当前时刻的状态;
  2. 传感器测量:通过观测设备获取当前时刻的带噪声测量值。

通过最优增益矩阵,卡尔曼滤波将预测值与测量值加权融合,得到统计意义上的最优估计。

数学模型

状态转移方程:

x_k = F * x_{k-1} + B * u_k + w_k

观测方程:

z_k = H * x_k + v_k

其中:
– x_k 为时刻k的状态向量
– F 为状态转移矩阵
– B 为控制输入矩阵
– u_k 为控制向量
– w_k 为过程噪声,协方差为Q
– z_k 为观测向量
– H 为观测矩阵
– v_k 为观测噪声,协方差为R

预测阶段

x̂_k|k-1 = F * x̂_{k-1|k-1} + B * u_k
P_k|k-1 = F * P_{k-1|k-1} * F^T + Q

更新阶段

K_k = P_k|k-1 * H^T * (H * P_k|k-1 * H^T + R)^(-1)
x̂_k|k = x̂_k|k-1 + K_k * (z_k - H * x̂_k|k-1)
P_k|k = (I - K_k * H) * P_k|k-1

二、Java实现

矩阵工具类

import java.util.Arrays;

/**
 * 简单矩阵运算工具类
 * 支持基本的矩阵乘法、转置、求逆(用于低维矩阵)
 */
public class Matrix {
    private final double[][] data;
    private final int rows;
    private final int cols;

    public Matrix(int rows, int cols) {
        this.rows = rows;
        this.cols = cols;
        this.data = new double[rows][cols];
    }

    public Matrix(double[][] data) {
        this.rows = data.length;
        this.cols = data[0].length;
        this.data = new double[rows][cols];
        for (int i = 0; i < rows; i++) {
            System.arraycopy(data[i], 0, this.data[i], 0, cols);
        }
    }

    public double get(int i, int j) { return data[i][j]; }
    public void set(int i, int j, double val) { data[i][j] = val; }
    public int getRows() { return rows; }
    public int getCols() { return cols; }

    /**
     * 矩阵乘法
     */
    public Matrix multiply(Matrix other) {
        if (this.cols != other.rows) {
            throw new IllegalArgumentException("矩阵维度不匹配");
        }
        Matrix result = new Matrix(this.rows, other.cols);
        for (int i = 0; i < this.rows; i++) {
            for (int j = 0; j < other.cols; j++) {
                double sum = 0;
                for (int k = 0; k < this.cols; k++) {
                    sum += this.data[i][k] * other.data[k][j];
                }
                result.set(i, j, sum);
            }
        }
        return result;
    }

    /**
     * 矩阵转置
     */
    public Matrix transpose() {
        Matrix result = new Matrix(cols, rows);
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < cols; j++) {
                result.set(j, i, data[i][j]);
            }
        }
        return result;
    }

    /**
     * 矩阵加法
     */
    public Matrix add(Matrix other) {
        if (this.rows != other.rows || this.cols != other.cols) {
            throw new IllegalArgumentException("矩阵维度不匹配");
        }
        Matrix result = new Matrix(rows, cols);
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < cols; j++) {
                result.set(i, j, this.data[i][j] + other.data[i][j]);
            }
        }
        return result;
    }

    /**
     * 矩阵减法
     */
    public Matrix subtract(Matrix other) {
        if (this.rows != other.rows || this.cols != other.cols) {
            throw new IllegalArgumentException("矩阵维度不匹配");
        }
        Matrix result = new Matrix(rows, cols);
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < cols; j++) {
                result.set(i, j, this.data[i][j] - other.data[i][j]);
            }
        }
        return result;
    }

    /**
     * 标量乘法
     */
    public Matrix scale(double scalar) {
        Matrix result = new Matrix(rows, cols);
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < cols; j++) {
                result.set(i, j, data[i][j] * scalar);
            }
        }
        return result;
    }

    /**
     * 2x2矩阵求逆(简化实现)
     */
    public Matrix inverse2x2() {
        if (rows != 2 || cols != 2) {
            throw new IllegalArgumentException("仅支持2x2矩阵求逆");
        }
        double a = data[0][0], b = data[0][1];
        double c = data[1][0], d = data[1][1];
        double det = a * d - b * c;
        if (Math.abs(det) < 1e-10) {
            throw new ArithmeticException("矩阵奇异,无法求逆");
        }
        Matrix result = new Matrix(2, 2);
        result.set(0, 0, d / det);
        result.set(0, 1, -b / det);
        result.set(1, 0, -c / det);
        result.set(1, 1, a / det);
        return result;
    }

    /**
     * 单位矩阵
     */
    public static Matrix identity(int n) {
        Matrix result = new Matrix(n, n);
        for (int i = 0; i < n; i++) {
            result.set(i, i, 1.0);
        }
        return result;
    }

    @Override
    public String toString() {
        StringBuilder sb = new StringBuilder();
        for (int i = 0; i < rows; i++) {
            sb.append("[");
            for (int j = 0; j < cols; j++) {
                sb.append(String.format("%.4f", data[i][j]));
                if (j < cols - 1) sb.append(", ");
            }
            sb.append("]\n");
        }
        return sb.toString();
    }
}

卡尔曼滤波器

/**
 * 卡尔曼滤波器实现
 * 支持任意维度的状态向量和观测向量
 */
public class KalmanFilter {

    // 状态向量维度
    private final int stateDim;
    // 观测向量维度
    private final int obsDim;

    // 状态转移矩阵 F
    private Matrix F;
    // 控制输入矩阵 B
    private Matrix B;
    // 观测矩阵 H
    private Matrix H;
    // 过程噪声协方差 Q
    private Matrix Q;
    // 观测噪声协方差 R
    private Matrix R;

    // 当前状态估计 x̂
    private Matrix x;
    // 当前估计误差协方差 P
    private Matrix P;

    public KalmanFilter(int stateDim, int obsDim) {
        this.stateDim = stateDim;
        this.obsDim = obsDim;

        // 初始化矩阵
        this.F = Matrix.identity(stateDim);
        this.B = new Matrix(stateDim, 1);
        this.H = new Matrix(obsDim, stateDim);
        this.Q = Matrix.identity(stateDim).scale(0.01);
        this.R = Matrix.identity(obsDim).scale(0.1);

        // 初始状态为0,协方差为单位矩阵
        this.x = new Matrix(stateDim, 1);
        this.P = Matrix.identity(stateDim);
    }

    /**
     * 设置状态转移矩阵
     */
    public void setF(Matrix F) {
        if (F.getRows() != stateDim || F.getCols() != stateDim) {
            throw new IllegalArgumentException("F矩阵维度错误");
        }
        this.F = F;
    }

    /**
     * 设置控制输入矩阵
     */
    public void setB(Matrix B) {
        if (B.getRows() != stateDim) {
            throw new IllegalArgumentException("B矩阵维度错误");
        }
        this.B = B;
    }

    /**
     * 设置观测矩阵
     */
    public void setH(Matrix H) {
        if (H.getRows() != obsDim || H.getCols() != stateDim) {
            throw new IllegalArgumentException("H矩阵维度错误");
        }
        this.H = H;
    }

    /**
     * 设置过程噪声协方差
     */
    public void setQ(Matrix Q) {
        if (Q.getRows() != stateDim || Q.getCols() != stateDim) {
            throw new IllegalArgumentException("Q矩阵维度错误");
        }
        this.Q = Q;
    }

    /**
     * 设置观测噪声协方差
     */
    public void setR(Matrix R) {
        if (R.getRows() != obsDim || R.getCols() != obsDim) {
            throw new IllegalArgumentException("R矩阵维度错误");
        }
        this.R = R;
    }

    /**
     * 设置初始状态
     */
    public void setInitialState(Matrix x0, Matrix P0) {
        this.x = x0;
        this.P = P0;
    }

    /**
     * 预测步骤
     * @param u 控制输入(可选,为null时忽略)
     */
    public void predict(Matrix u) {
        // x̂_k|k-1 = F * x̂_{k-1|k-1} + B * u
        x = F.multiply(x);
        if (u != null && B != null) {
            x = x.add(B.multiply(u));
        }

        // P_k|k-1 = F * P_{k-1|k-1} * F^T + Q
        P = F.multiply(P).multiply(F.transpose()).add(Q);
    }

    /**
     * 更新步骤
     * @param z 观测值
     */
    public void update(Matrix z) {
        // 计算卡尔曼增益
        // K = P * H^T * (H * P * H^T + R)^(-1)
        Matrix S = H.multiply(P).multiply(H.transpose()).add(R);
        Matrix K = P.multiply(H.transpose()).multiply(S.inverse2x2());

        // x̂_k|k = x̂_k|k-1 + K * (z - H * x̂_k|k-1)
        Matrix y = z.subtract(H.multiply(x)); // 新息(观测残差)
        x = x.add(K.multiply(y));

        // P_k|k = (I - K * H) * P_k|k-1
        Matrix I = Matrix.identity(stateDim);
        P = I.subtract(K.multiply(H)).multiply(P);
    }

    /**
     * 执行一次完整的预测-更新循环
     */
    public void step(Matrix u, Matrix z) {
        predict(u);
        update(z);
    }

    public Matrix getState() { return x; }
    public Matrix getCovariance() { return P; }

    public static void main(String[] args) {
        // 示例:一维匀速运动模型
        // 状态向量:[位置, 速度]^T
        // 观测:仅测量位置

        KalmanFilter kf = new KalmanFilter(2, 1);

        // 时间间隔
        double dt = 1.0;

        // 设置状态转移矩阵 F = [[1, dt], [0, 1]]
        Matrix F = new Matrix(new double[][]{{1, dt}, {0, 1}});
        kf.setF(F);

        // 设置观测矩阵 H = [[1, 0]]
        Matrix H = new Matrix(new double[][]{{1, 0}});
        kf.setH(H);

        // 设置过程噪声 Q
        Matrix Q = new Matrix(new double[][]{{0.1, 0}, {0, 0.01}});
        kf.setQ(Q);

        // 设置观测噪声 R
        Matrix R = new Matrix(new double[][]{{0.5}});
        kf.setR(R);

        // 初始状态
        Matrix x0 = new Matrix(new double[][]{{0}, {1}}); // 位置0,速度1
        Matrix P0 = Matrix.identity(2).scale(10);
        kf.setInitialState(x0, P0);

        // 模拟10步
        System.out.println("卡尔曼滤波跟踪结果:");
        System.out.println("步骤\t真实位置\t观测位置\t估计位置\t估计速度\t位置方差");

        double truePos = 0;
        double trueVel = 1;
        java.util.Random rand = new java.util.Random(42);

        for (int k = 1; k <= 10; k++) {
            // 真实状态演进
            truePos += trueVel * dt;

            // 带噪声的观测
            double measurement = truePos + rand.nextGaussian() * Math.sqrt(0.5);
            Matrix z = new Matrix(new double[][]{{measurement}});

            // 执行卡尔曼滤波
            kf.step(null, z);

            Matrix state = kf.getState();
            Matrix cov = kf.getCovariance();

            System.out.printf("%d\t%.3f\t\t%.3f\t\t%.3f\t\t%.3f\t\t%.4f%n",
                k, truePos, measurement, state.get(0, 0), state.get(1, 0), cov.get(0, 0));
        }
    }
}

三、应用场景扩展

卡尔曼滤波可以扩展到多维场景。以下是一个二维目标跟踪的示例配置:

/**
     * 二维目标跟踪配置示例
     * 状态:[x, y, vx, vy]^T
     * 观测:[x, y]^T
     */
    public static KalmanFilter create2DTracker(double dt) {
        KalmanFilter kf = new KalmanFilter(4, 2);

        // F = [[1, 0, dt, 0], [0, 1, 0, dt], [0, 0, 1, 0], [0, 0, 0, 1]]
        Matrix F = new Matrix(new double[][]{
            {1, 0, dt, 0},
            {0, 1, 0, dt},
            {0, 0, 1, 0},
            {0, 0, 0, 1}
        });
        kf.setF(F);

        // H = [[1, 0, 0, 0], [0, 1, 0, 0]]
        Matrix H = new Matrix(new double[][]{
            {1, 0, 0, 0},
            {0, 1, 0, 0}
        });
        kf.setH(H);

        // Q:过程噪声
        Matrix Q = Matrix.identity(4).scale(0.05);
        kf.setQ(Q);

        // R:观测噪声(GPS精度约5米)
        Matrix R = Matrix.identity(2).scale(5.0);
        kf.setR(R);

        return kf;
    }

四、复杂度分析

  • 预测阶段:主要开销为矩阵乘法F×P×F^T,时间复杂度为O(n³),n为状态维度;
  • 更新阶段:主要开销为求逆运算,对于m维观测,复杂度为O(m³);
  • 整体:单次迭代复杂度为O(n³ + m³)。对于常见应用(n≤10, m≤5),计算量极小,完全满足实时性要求。

五、总结

本文用Java完整实现了卡尔曼滤波算法,核心要点包括:

  1. 预测-更新框架:清晰分离系统模型预测与观测融合两个步骤;
  2. 最优增益推导:通过最小化估计误差协方差,得到理论最优的卡尔曼增益;
  3. 矩阵运算实现:构建基础的矩阵工具类,支撑滤波器的核心计算;
  4. 多维扩展:从一维匀速模型自然扩展到二维目标跟踪等复杂场景。

卡尔曼滤波是状态估计领域的基石算法,理解它对于后续学习扩展卡尔曼滤波(EKF)、无迹卡尔曼滤波(UKF)以及粒子滤波等高级方法具有重要的铺垫意义。