卡尔曼滤波器进阶：二维线性动态系统与眼动追踪应用

摘要

本文深入探讨了二维卡尔曼滤波器的数学原理及其在眼动追踪数据分析中的应用。我们将从线性动态系统(LDS)的基础理论出发，逐步构建卡尔曼滤波器的实现，并通过实际眼动追踪数据展示其平滑效果。

1. 线性动态系统(LDS)基础

1.1 系统定义

线性动态系统由以下方程描述：

状态方程： sₜ = Dsₜ₋₁ + wₜ

观测方程： mₜ = Hsₜ + ηₜ

其中：

• sₜ ∈ ℝⁿ表示t时刻的隐状态
• mₜ ∈ ℝᵐ表示t时刻的观测值
• D ∈ ℝⁿˣⁿ是状态转移矩阵
• H ∈ ℝᵐˣⁿ是观测矩阵
• wₜ ~ N(0,Q)是过程噪声
• ηₜ ~ N(0,R)是观测噪声

1.2 系统采样实现

我们可以通过以下步骤生成LDS的样本轨迹：

def sample_lds(n_timesteps, params, seed=0):
    # 初始化随机数生成器
    np.random.seed(seed)
    
    # 预生成噪声样本
    mi = stats.multivariate_normal(cov=params['Q']).rvs(n_timesteps)
    eta = stats.multivariate_normal(cov=params['R']).rvs(n_timesteps)
    
    # 初始化状态和观测数组
    state = np.zeros((n_timesteps, params['D'].shape[0]))
    obs = np.zeros((n_timesteps, params['H'].shape[0]))
    
    # 模拟系统动态
    for t in range(n_timesteps):
        if t == 0:
            # 初始状态采样
            state[t] = stats.multivariate_normal(
                mean=params['mu_0'], 
                cov=params['sigma_0']).rvs()
        else:
            # 状态转移
            state[t] = params['D'] @ state[t-1] + mi[t]
        
        # 生成观测值
        obs[t] = params['H'] @ state[t] + eta[t]
    
    return state, obs

2. 卡尔曼滤波器原理

卡尔曼滤波器通过递归方式估计系统状态，包含两个主要步骤：

1. 预测步骤：

   • 预测下一状态：ŝₜ|ₜ₋₁ = Dŝₜ₋₁|ₜ₋₁
   • 预测误差协方差：Pₜ|ₜ₋₁ = DPₜ₋₁|ₜ₋₁Dᵀ + Q

2. 更新步骤：

   • 计算卡尔曼增益：Kₜ = Pₜ|ₜ₋₁Hᵀ(HPₜ|ₜ₋₁Hᵀ + R)⁻¹
   • 更新状态估计：ŝₜ|ₜ = ŝₜ|ₜ₋₁ + Kₜ(mₜ - Hŝₜ|ₜ₋₁)
   • 更新协方差估计：Pₜ|ₜ = (I - KₜH)Pₜ|ₜ₋₁

3. 眼动追踪数据应用

3.1 数据特点

眼动追踪数据通常包含：

• 二维坐标(x,y)表示注视点位置
• 高频采样(通常100-1000Hz)
• 存在测量噪声和微小抖动

3.2 卡尔曼滤波实现

def kalman_filter(observations, params):
    n_timesteps = observations.shape[0]
    n_dim_state = params['D'].shape[0]
    
    # 初始化估计值和协方差
    state_est = np.zeros((n_timesteps, n_dim_state))
    cov_est = np.zeros((n_timesteps, n_dim_state, n_dim_state))
    
    # 初始状态
    state_est[0] = params['mu_0']
    cov_est[0] = params['sigma_0']
    
    for t in range(1, n_timesteps):
        # 预测步骤
        state_pred = params['D'] @ state_est[t-1]
        cov_pred = params['D'] @ cov_est[t-1] @ params['D'].T + params['Q']
        
        # 更新步骤
        residual = observations[t] - params['H'] @ state_pred
        S = params['H'] @ cov_pred @ params['H'].T + params['R']
        K = cov_pred @ params['H'].T @ np.linalg.inv(S)
        
        state_est[t] = state_pred + K @ residual
        cov_est[t] = (np.eye(n_dim_state) - K @ params['H']) @ cov_pred
    
    return state_est, cov_est

3.3 效果评估

我们使用真实眼动数据比较原始观测和滤波后轨迹：

# 加载眼动数据
subjects, images = load_eyetracking_data()
gaze_data = subjects[0]['fixations'][0]  # 第一个受试者的第一次注视数据

# 设置卡尔曼滤波器参数
kf_params = {
    'D': 0.95 * np.eye(2),
    'Q': 0.1 * np.eye(2),
    'H': np.eye(2),
    'R': 1.0 * np.eye(2),
    'mu_0': gaze_data[0],
    'sigma_0': 0.1 * np.eye(2)
}

# 应用卡尔曼滤波
smoothed_trajectory, _ = kalman_filter(gaze_data, kf_params)

# 可视化结果
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
plot_gaze_data(gaze_data, img=images[0], ax=ax1)
ax1.set_title("原始眼动数据")
plot_gaze_data(smoothed_trajectory, img=images[0], ax=ax2)
ax2.set_title("卡尔曼滤波后")
plt.show()

4. 参数调优建议

1. 状态转移矩阵D：

   • 接近1的值表示状态变化缓慢
   • 对角线值可以不同以适应x,y方向的不同动态特性

2. 过程噪声Q：

   • 控制状态变化的灵活性
   • 值越大，滤波器对观测值变化越敏感

3. 观测噪声R：

   • 反映测量设备的精度
   • 值越大，滤波器对观测值的信任度越低

5. 总结

卡尔曼滤波器为处理噪声环境下的动态系统提供了强大的工具，特别适用于：

• 眼动追踪数据分析
• 目标跟踪系统
• 传感器融合应用
• 金融时间序列分析

通过适当调整参数，可以在平滑噪声和保持信号快速变化之间取得良好平衡。二维扩展使卡尔曼滤波器能够处理更复杂的空间动态问题。