卡尔曼滤波在 GNSS 时间序列上的应用

GNSS 时间序列的特点

• 每日一个解，含位置 + 不确定度
• 短期噪声：多路径、电离层残差，±3~10 mm
• 中期漂移：天线墩温度变形、参考框架调整
• 长期趋势：构造运动、边坡位移

直接看日解时间序列，趋势被噪声掩盖。做卡尔曼滤波可以同时估计位置、速度、加速度，输出平滑趋势。

状态向量

位移监测的状态向量：

• x: 位置
• v: 速度
• a: 加速度

状态方程（匀加速运动）：

• x[k] = x[k-1] + v[k-1] * dt + 0.5 * a[k-1] * dt²
• v[k] = v[k-1] + a[k-1] * dt
• a[k] = a[k-1]

写成矩阵形式： `` [x] [1 dt 0.5dt²] [x] [noise_x] [v] = [0 1 dt ] [v] + [noise_v] [a] [0 0 1 ] [a] [noise_a] ``

观测方程

GNSS 日解只观测位置 x： `` [z] = [1 0 0] [x v a]^T + noise_z ``

z 是 GNSS 测得的位置，x 是真实位置，观测噪声是 GNSS 解的不确定度。

卡尔曼滤波步骤

预测

`` x_pred = A * x_est P_pred = A * P * A^T + Q `` A 是状态转移矩阵，P 是状态协方差，Q 是过程噪声。

更新

`` K = P_pred * H^T * (H * P_pred * H^T + R)^-1 x_est = x_pred + K * (z - H * x_pred) P = (I - K * H) * P_pred `` K 是卡尔曼增益，H 是观测矩阵，R 是观测噪声协方差。

参数调优

Q（过程噪声）

Q 描述"系统状态自身变化的不确定度"：

• Q 大：滤波跟踪快，噪声压不下去
• Q 小：滤波平滑，但响应慢

GNSS 位移监测的合理 Q：

• 位置噪声 0.1 mm²/天（边坡稳定，位移慢）
• 速度噪声 0.01 mm²/天²
• 加速度噪声 0.001 mm²/天³

R（观测噪声）

R 是 GNSS 日解的不确定度（标准差平方）：

• 日解 RMS 3 mm → R = 9 mm²
• 如果接收机开放了不确定度估计，用每个观测的实际 R

实战代码（Python）

```python import numpy as np

def kalman_gnss(positions, uncertainties, dt=1.0): """GNSS 位移时间序列的卡尔曼滤波

positions: 日解位置数组 (mm) uncertainties: 每日解的不确定度 (mm) dt: 时间步长（天） """ n = len(positions) # 状态向量 [x, v, a] x = np.array([positions[0], 0, 0]) P = np.diag([uncertainties[0]**2, 1, 1])

A = np.array([ [1, dt, 0.5 * dt**2], [0, 1, dt], [0, 0, 1] ]) H = np.array([[1, 0, 0]]) Q = np.diag([0.1, 0.01, 0.001])

filtered = np.zeros(n) velocities = np.zeros(n)

for k in range(n): # 预测 x_pred = A @ x P_pred = A @ P @ A.T + Q

# 更新 R = uncertainties[k]**2 K = P_pred @ H.T / (H @ P_pred @ H.T + R) z = positions[k] x = x_pred + K.flatten() * (z - H @ x_pred) P = (np.eye(3) - K @ H) @ P_pred

filtered[k] = x[0] velocities[k] = x[1]

return filtered, velocities

模拟：真实位移 0.01 mm/天，GNSS 噪声 ±3 mm

np.random.seed(42) true_pos = np.cumsum(np.full(100, 0.01)) # 1 mm 总位移 noisy = true_pos + np.random.randn(100) * 3 unc = np.full(100, 3.0)

filtered, vels = kalman_gnss(noisy, unc) print(f"噪声 RMS: {np.std(noisy - true_pos):.2f} mm") print(f"滤波后 RMS: {np.std(filtered - true_pos):.2f} mm") print(f"速度估计（最后10天平均）: {np.mean(vels[-10:]):.4f} mm/天") print(f"真实速度: 0.01 mm/天") ```

调参经验

Q 太大

滤波"跟踪噪声"，输出仍抖。降噪效果差。

Q 太小

滤波过度平滑，对真实变形响应慢。加速变形阶段可能漏掉。

折中

边坡监测推荐 Q = (0.1, 0.01, 0.001)，对应位置噪声 0.3 mm/天，速度噪声 0.1 mm/天，加速度噪声 0.03 mm/天²。

自适应 Q

高变形期 Q 自适应增大，低变形期 Q 减小。基于"残差大小"自适应调整：

• 残差大（实际位移超预测）→ Q 增大，跟踪快
• 残差小 → Q 减小，平滑好

常见错

错一：用固定 R

GNSS 解的不确定度日变化大（多路径强的天 R 大）。用固定 R 会"信任坏数据"或"不信任好数据"。正确做法：用每个观测的实际不确定度。

错二：忽略粗差

GNSS 日解偶尔有粗差（多路径 + 维护操作），卡尔曼会被"拉偏"。要在滤波前做粗差剔除（如 Hampel 滤波）。

错三：状态向量过简

只用 [x, v] 不带加速度，对加速变形响应慢。边坡监测推荐 [x, v, a] 三状态。

结语

卡尔曼滤波能把 GNSS 日解噪声从 ±3 mm 降到 ±1 mm，并实时估计位移速率。关键是Q 与 R 的合理设置——Q 太大跟踪噪声，太小漏掉变形。边坡监测推荐 Q = (0.1, 0.01, 0.001)，并配合粗差剔除预处理。