粒子滤波实现MATLAB目标追踪

111 阅读1分钟

基于粒子滤波的运动目标跟踪MATLAB实现

粒子滤波(Particle Filter)是一种基于蒙特卡罗方法的非线性滤波技术,适用于目标跟踪中的状态估计问题。以下是完整的MATLAB实现方案:

1. 核心算法原理

设目标状态为 xt=[px,py,vx,vy]Tx_t = [p_x, p_y, v_x, v_y]^T,包含位置和速度分量。粒子滤波通过以下步骤迭代:

  1. 状态预测xt(i)=f(xt1(i))+wt(i),wt(i)N(0,Q)x_t^{(i)} = f(x_{t-1}^{(i)}) + w_t^{(i)},\quad w_t^{(i)} \sim \mathcal{N}(0,Q) 其中 f()f(\cdot) 是状态转移函数

  2. 权重更新wt(i)p(ztxt(i))w_t^{(i)} \propto p(z_t|x_t^{(i)}) 基于观测似然函数

  3. 重采样: 根据权重 {wt(i)}i=1N\{w_t^{(i)}\}_{i=1}^N 重新生成粒子集

2. MATLAB实现代码

% 粒子滤波目标跟踪主函数
function [estimated_pos] = particle_filter_tracking(observations, N_particles)
    % 初始化参数
    dt = 1;                      % 时间步长
    Q = diag([0.1, 0.1, 0.01, 0.01]); % 过程噪声协方差
    R = diag([1, 1]);             % 观测噪声协方差
    
    % 初始化粒子集
    particles = zeros(4, N_particles);
    weights = ones(1, N_particles)/N_particles;
    estimated_pos = zeros(2, size(observations,2));
    
    % 主循环
    for t = 1:size(observations,2)
        % 1. 状态预测
        particles = state_transition(particles, dt) + mvnrnd(zeros(4,1), Q, N_particles)';
        
        % 2. 权重更新
        weights = update_weights(particles, observations(:,t), R);
        
        % 3. 状态估计
        estimated_pos(:,t) = particles(1:2,:) * weights';
        
        % 4. 重采样
        [particles, weights] = systematic_resampling(particles, weights);
    end
end

% 状态转移函数 (匀速模型)
function new_particles = state_transition(particles, dt)
    A = [1 0 dt 0; 
         0 1 0 dt; 
         0 0 1 0; 
         0 0 0 1];
    new_particles = A * particles;
end

% 权重更新函数
function weights = update_weights(particles, observation, R)
    predicted_obs = particles(1:2,:);
    innovation = observation - predicted_obs;
    mahalanobis = sum(innovation .* (R \ innovation), 1);
    weights = exp(-0.5 * mahalanobis);
    weights = weights / sum(weights); % 归一化
end

% 系统重采样函数
function [new_particles, new_weights] = systematic_resampling(particles, weights)
    N = size(particles, 2);
    cum_weights = cumsum(weights);
    new_indices = zeros(1, N);
    u = ([0:N-1] + rand(1)) / N;
    
    j = 1;
    for i = 1:N
        while u(i) > cum_weights(j)
            j = j + 1;
        end
        new_indices(i) = j;
    end
    
    new_particles = particles(:, new_indices);
    new_weights = ones(1, N) / N;
end

3. 使用示例

% 生成测试数据 (真实轨迹 + 噪声观测)
T = 100;
true_trajectory = [linspace(0,10,T); linspace(0,5,T)];
observations = true_trajectory + randn(2,T)*2;

% 运行粒子滤波器
N_particles = 1000;
estimated_pos = particle_filter_tracking(observations, N_particles);

% 可视化结果
figure;
plot(true_trajectory(1,:), true_trajectory(2,:), 'b-', 'LineWidth', 2); hold on;
plot(observations(1,:), observations(2,:), 'ro', 'MarkerSize', 3);
plot(estimated_pos(1,:), estimated_pos(2,:), 'g--', 'LineWidth', 2);
legend('真实轨迹', '观测数据', '粒子滤波估计');
title('运动目标跟踪结果');
xlabel('X坐标'); ylabel('Y坐标');
grid on;

4. 关键参数说明

  • 粒子数量N_particles 建议在 500-2000 之间,平衡精度与计算效率
  • 过程噪声Q 矩阵影响状态预测的不确定性
  • 观测噪声R 矩阵决定观测数据的可信度
  • 重采样策略:采用系统重采样避免粒子退化

5. 性能优化建议

  1. 自适应粒子数:根据有效样本大小 Neff=1/(wi2)N_{eff} = 1/\sum(w_i^2) 动态调整粒子数
  2. 混合建议分布:结合最新观测信息改进状态预测
  3. 并行计算:利用 parfor 加速粒子预测步骤

注意:实际应用中需根据目标运动特性(如加速度模型)修改状态转移函数,并根据传感器特性(如雷达/视觉数据)调整观测似然函数。