2.1 概述 卡尔曼滤波是有卡尔曼博士为满足计算机用于人造地球卫星定轨和导航等计算的要求， 在 1960 年提出的一种新的线性滤波模型，通常称他为卡尔曼滤波的方法。卡尔曼滤波是一 种具有无偏性的递推线性最小方差估计，及估计误差的均值或数学期望为零。在计算方法上， 卡尔曼滤波采用递推形式，即在 t-1 时刻估值的基础上，利用 t 时刻的观测值，递推得到 t 时刻的状态估值。由于一次仅处理一个时刻的观测值，无需存储先前的观测数据，因而计算 量大大减少。卡尔曼滤波能有效地消除干扰噪声，获得逼近真实情况的有用信息。这种滤波 方法很适合处理动态系统的观测数据，因而广泛用于动态测量、飞行器导航、导弹制导和工 业自动化等。 卡尔曼滤波作为一种比较成熟的熟悉方法，近几年来也应用于监测网的形变分析。利用 新的观测值，通过不断地预测和修正即可估计出系统新的状态值，这正使用预处理多期重复 观测的 GPS 监测网的观测数据，从而对监测网进行形变分析。在系统噪声和观测噪声均为 白噪声序列的情况下，采用卡尔曼滤波可以估计和预测边坡的三维变形量及形变速率。由于 边坡上监测点的型变量是随着时间不断变化的状态参数，因而利用卡尔曼滤波可以估计这些 状态参数，使它们能更全面地反映边坡的运动状态。 2.2 卡尔曼滤波方法 卡尔曼滤波借助于系统建模的状态矩阵和观测数据，利用数据来估计随时间不断变化的 状态向量，实时、最优的估计系统的状态量，并对未来时刻的系统状态量进行预报。利用卡 尔曼滤波进行状态估计，首先要建立系统的状态方程和观测方程。 假设连续线性系统有如下状态方程和观测方程，即 式中，X(t)为 t 时刻系统状态向量，A(t)为 t 时刻系统状态矩阵，F(t)为 t 时刻系统状态噪声输 入矩阵，Z(t)为 t 时刻系统观测向量，H(t)为 t 时刻系统观测矩阵，W(t)为 t 时刻系统状态噪 声，V(t)为 t 时刻系统状态观测噪声。 式中，Q(t)为系统状态噪声 W(t)的方差阵；R(t)为系统观测噪声 V(t)的方差阵；δ(t-t’)为 Dirac 函数，且有（在 0 点的值为无穷大，其余点的值为 0） X t =AtXt +FtWt Zt =HtXt +Vt 假设状态噪声和观测噪声满足约束条件EWt =0 EVt =0 EWtVT(t') =0 EWtWT(t') =Qtδ(t−t') EVtVT(t') =Rtδ(t−t') δt−t' = ∞,t=t' 0,t≠t' t)dt=1 ∞δ( −∞ Xk=Φk,k−1Xk−1+Γk,k−1Wk−1 Zk=HkXk+Vk Φk,k−1=I+A(tk)Δtk+o(Δtk) 将(2.1)式离散化可得 且有 (2.1) (2.2) (2.3) (2.4) (2.5) (2.6) 

Γk,k−1= I+AtkΔtk⋯FtkΔtk=FtkΔtk+o(Δtk) 式(2.5)中，Xk为 k 时刻系统状态向量，Φk,k−1为 k-1 到 k 时刻的系统状态转移矩阵，Γk,k−1 为 k-1 到 k 时刻系统状态噪声输入矩阵，Zk为 k 时刻系统观测向量，Hk为 k 时刻系统的观测 矩阵，Vk为 k 时刻系统观测噪声。 转移矩阵：转移概率矩阵是俄国数学家马尔科夫提出的，他在 20 世纪初发现：一个系统的 某些因素在转移中，第 n 次结果只受第 n-1 的结果影响，即只与当前所处状态有关，而与过 去状态无关。 在马尔科夫分析中，引入状态转移这个概念。所谓状态是指客观事物可能出 (2.7) 现或存在的状态；状态转移是指客观事物由一种状态转移到另一种状态的概率。 对应于你提的天气预报的问题，若天气状态转移概率表如下： 今天 明天 晴 阴 雨 晴 3/4 1/2 1/4 阴 1/8 1/4 1/2 雨 1/8 1/4 1/4 其中转移矩阵 A 的每一个元素都表示从今天的一种状态到明天的一种状态的概率，例如， 第 2 行第 3 列(最后一列)的值为 1/2,，这给出了今天下雨而明天转阴的概率。所以称 A 为转 移矩阵。 状态噪声和观测噪声满足如下约束条件，即E(Wk)=0 E(Vk)=0 E(WkVjT)=0 E(WkWjT)=Qkδkj E(VkVjT)=Rkδkj 式中，δkj为 Kronecker 函数 δkj= 1,k=j 0,k≠j Qk为系统状态噪声Wk的方差阵，即 Qk=Q(t)∆tk Rk为系统观测噪声Vk的方差阵，即 Rk=R(t)∆tk Xk，k−1=Φk,k−1Xk−1 Xk=Xk,k−1+Kk(Zk−HkXk,k−1) Kk=Pk,k−1HkT(HkPk,k−1HkT+Rk)−1 状态一步预测误差协方差阵Pk,k−1=Φk,k−1Pk−1Φk，k−1 T 因此，卡尔曼滤波由如下递推方程组成。 状态一步预测 状态量估计 滤波增益矩阵 (2.8) (2.9) (2.10) (2.11) (2.12) （2.13） （2.14） （2.15） (2.16） +Γk,k−1Qk−1Γk，k−1 T 

状态估计误差协方差阵 Pk=(I−KkHk)Pk,k−1(I−KkHk)T+KkHkKkT 由式(2.13)~式(2.17)可以看到，卡尔曼滤波方程式一组递推计算公式，整个计算过程是 一个不断的预测、修正的过程。在估计状态量是不需要存储大量的观测数据，并且当得到的 新的观测数据时，可以随时得到状态量新的滤波值，从而便于实时处理观测结果。图显示了 卡尔曼滤波算法的框图。 (2.17） 卡尔曼滤波方法具有以下特点： 1）卡尔曼滤波递推公式是基于广义最小二乘原理得到的，因而它具有无偏和方差最小 的特点，所得到的估值比最小二乘估值具有更高的精度。 不依赖于状态初始值的选取，因而卡尔曼滤波稳态值与初始值选取无关。 2）滤波递推充分多步后，滤波误差协方差阵Pk将不依赖与初始值P0，状态滤波结果也 3）增益矩阵Kk与观测值无关，因而可以预先计算，从而可减少实时处理的计算工作量。 的Qk和Rk矩阵。 4）应用卡尔曼滤波算法，除了需要精确地建立动态方程和观测方程，还需要选取合适 给定初值 X0，P0，k=1 Φ，Г，Q H，R 状态一步预测误差方程 +Γk,k−1Qk−1Γk，k−1 Pk,k−1=Φk,k−1Pk−1Φk，k−1 T T Kk=Pk,k−1HkT(HkPk,k−1HkT+Rk)−1 Pk=(I−KkHk)Pk,k−1(I−KkHk)T 最优滤波增益矩阵 滤波误差方程阵 进一步状态预测 K=k+1 Zk Xk，k−1=Φk,k−1Xk−1 最优滤波值 Xk=Xk,k−1+Kk(Zk−HkXk,k−1) 图 2.1 卡尔曼滤波算法框图 2.3 基于卡尔曼滤波的边坡变形预报 2.3.1 状态方程的建立 

假设边坡上某监测点 i 在 t 时刻的三维形变量为si(t)，瞬时形变速度为vi(t)，并将瞬时 速度看做一阶马尔科夫的随机过程，则有如下方程，即 式中，τ是相关的时间常数。将边坡上监测点的三维形变量和三维形变速度考虑为状态向量 X(t)，即 sit =vi(t) vit =−1τvit +W(t) siz(t) vix(t) viy(t) viz(t)]T siy(t) Xt =[six(t) Wt =[ωx ωy ωz]T I3×3 Xt = O3×3 Xt + O3×3 I3×3 W(t) O3×3 −1τI3×3 Xk=Φk,k−1Xk−1+Γk,k−1Wk−1 Φk,k−1= I3×3 ΔtkI3×3 (1−1τΔtk)I3×3 O3×3 Γk,k−1= O3×3 ΔtkI3×3 sixtsiytsizt = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 Zt = sixtsixtsizt Zk+1=HXk+1+Vk+1 0 0 1 0 0 0 H= 1 0 0 0 1 0 0 0 0 0 0 0 sixtsiytsiztvixtviytvizt + VxVyVz (2.18) (2.19) (2.20) (2.21) (2.22) (2.23) (2.24) (2.25) 且 则式(2.18)可以写成 将(2.20)离散化得到 且有 令 式中 则离散化的观测方程为 2.3.2 观测方程的建立 考虑边坡上某监测点 i 在 t 时刻的观测方程为 式(2.24)和式(2.25)组成了边坡形变量鼓励及的卡尔曼滤波状态方程和观测方程。 2.3.3 卡尔曼滤波初值确定 从卡尔曼滤波递推公式和图 2.1 可以看出，要确定系统在tk时刻的状态，必须要先知道 系统在t0时刻的初始值。在实际应用中，系统的初始值是难以精确先知的，一般只能近似的 

给定。但如果给定初始值偏差较大，那么卡尔曼滤波的估计值收敛速度就较慢。因此，合理 的给定系统初值是很重要的。 卡尔曼滤波的初始值包括初始状态向量X0和初始协方差阵P0，对于便携监测系统分两种 情况来确定滤波初始值。 1）将监测点的初始值作为状态初始值，适用于变形观测时间间隔较短的连续观测情况。 2）将监测点的三维形变量和在相应方向上的形变速率作为状态初值，适用于便携观测 时间间隔相对较长的周期性观测的情况。