基于变维滤波算法的Kalman最佳平滑器在滑坡监测数据处理中的应用

其中 S (k ) 是 v (k ) 的协方差矩阵，由于 v (k ) 是零均值高斯随机变量，所以 δ (k ) 服从 nz (nz 是测 量的维数) 个自由度的 χ2 分布，µ(k ) 亦服从 χ2 分布，且
k→∞
lim E [µ(k )] = nz /(1 − α).
588
工
程
数
学
学
报
第26卷
ˆ x ¨m (k − ∆/k − ∆) = 2[z (k − ∆) − z ˆ(k − ∆/k − ∆ − 1)]/T 2 , ˆ ˆ x ˙ m (k − ∆/k − ∆) = x ˙ (k − ∆ − 1/k − ∆ − 1) + T · x ¨m (k − ∆/k − ∆), ˆ x ¨m (k − ∆/k − ∆) = z (k − ∆). 协方差矩阵修正为 pm 11 (k − ∆/k − ∆) = R11 , pm 12 (k − ∆/k − ∆) = 2R11 /T,
第26卷 第4期 2009年08月
工
程
数
学
学
报
CHINESE JOURNAL OF ENGINEERING MATHEMATICS
Vol. 26 No. 4 Aug. 2009
文章编号:1005-3085(2009)04-0586-07
基于变维滤波算法的 Kalman 最佳平滑器 在滑坡监测数据处理中的应用
(16)
+8p12 (k − ∆ − 1/k − ∆ − 1)/T 3 +4p22 (k − ∆ − 1/k − ∆ − 1)/T 2 . 这样，由式 (2)-(17) 便构成了对滑坡自动跟踪的自适应变维滤波算法。 (17)
第4期
唐春艳，彭继兵：基于变维滤波算法的 Kalman 最佳平滑器在滑坡监测数据处理中的应用
关键词: Kalman 最佳平滑器；变维滤波；滑坡；数据处理 分类号: AMS(2000) 62M20
1
引言
滑坡是一个复杂的地质现象，且易受到各种诸如外界环境、人为误差、监测方法的不妥等 因素影响，所以必须对所获得的滑坡监测数据进行预处理，降低数据中的噪声和伪信息的影 响，以提高滑坡预报的精度。 滑坡的发展变化过程一般经历缓慢变形—匀速变形—加速变形—急剧变形四个阶段， 将 Kalman 滤波理论应用于滑坡变形监测的数据处理，传统方法是建立一个固定的状态方 程和观测方程来描述滑坡变形，并用标准 Kalman 滤波对滑坡位移监测数据进行跟踪滤波 去噪[1-3] 。该方法针对滑坡在匀速变形阶段是有效的，但当滑坡进入加速变形阶段时，标 准 Kalman 滤波可能会出现发散，所以需要采用自适应方法。 本文作者曾提出了基于交互多模算法的 Kalman 平滑器处理监测数据[4] ，虽然该方法能较好 处理滑坡监测数据，但是该方法使用的 Markov 链的转移概率矩阵不容易确定。 为此，经过进一步的研究，本文首次将变维滤波 (VD) 算法引入到滑坡监测数据处理中来， 为了提高滑坡预报精度并在此基础上给出了一种 Kalman 滤波平滑算法，该方法简洁明了，而 且也为寻找滑坡形变的变点提供了新的思路。
k
µa (k ) =
j =k−∆−1
δa (j ),
−1
(7) a ˆ(j/j ). (8)
m δa (j ) = a ˆT (j/j ) Pa (j/j )
m a ˆ(j/j ) 为加速度的估计值，Pa (j/j ) 为协方差矩阵的对应块，如果 µa (k ) < Ta ，则加速度估计 无显著性意义，滤波器退出机动模型，其中 Ta 的意义与 Th 相同。 3.2 滤波器初始化 当机动检测器在第 k 时刻检测到机动时，则滤波器假定在 k − ∆ − 1 开始有一恒定的加速， 在窗口 ∆ 内的状态估计应修正如下。
Z (k ) 为滑坡在匀速变形阶段的位置，H = [ 1 0 ]，V (k ) 为零均值、方差为 σ 2 的 Gauss 白噪声 序列，且与 W (k ) 不相关。 当滑坡进入加速和急剧加速变形的机动阶段，将加速度作为状态变量的一个分量，并视为 恒定的，将滑坡的地质本身和外界环境等因素引起的直线轨迹的偏离视为加速度的扰动，便可 建立滑坡在加速变形阶段的状态方程 X m (k + 1) = Φm X m (k ) + Γm W m (k ). 其中， m Xm = x ˙ x ¨m xm , Φm = 0 0 1 T 1 0 T 2 /2 T 1 , Γ = T /2 1 T 2 /4 , (4)
式 (23)-(25) 中 j = N − 1, N − 2, · · · , 0。这三式便构成了一套最佳固定区间平滑的逆向递推算 ˆ (N | N ) 和 P (N | N )，可由式 (18)-(22) 给出。 法，初值为 X
5
基 于 变 维 滤波 算 法 的 Kalman 最佳 平 滑 器
由 3.1 介绍的机动检测可知，通过机动监测器能够检测到滑坡进入机动的开始点和结束点， 也能够检测到非机动的开始点和结束点。这样便可在相应的机动区间和非机动区间分别进 行 Kalman 最佳固定区间平滑，最终获得滑坡在整个观测区间上的平滑值。具体算法步骤为： 1) 假定滑坡开始处于匀速变形阶段，故令非机动第一个开始点为 1，滤波器采用非机动模 型跟踪滑坡； 2) 在非机动阶段利用机动监测器检测机动，若在 k 时刻检测到机动，则立即对机动模型 进行初始化，并采用机动模型跟踪滑坡，相应地记下非机动第一个结束点 k ，机动第一个开始 点 k + 1，下转到 4)； 3) 若机动监测器在 k 时刻没有检测到机动，则继续利用非机动模型跟踪滑坡，上转 到 2)； 4) 在机动阶段利用机动监测器检测非机动，若在 k 时刻检测到非机动，则立即对非机动 模型进行初始化，并采用非机动模型跟踪滑坡，相应地记下机动第一个结束点 k ，非机动第二 个开始点 k + 1，上转到 2)； 5) 若机动监测器在 k 时刻没有检测到非机动，则继续利用机动模型跟踪滑坡，上转 到 4)； 6) 重复 2)-5)，直到 k = M ，M 为观测数据长度； 7) 分别在相应的非机动和机动阶段利用 Kalman 最佳固定区间平滑器对其进行平滑处理。
589
3.3 标准 Kalman 滤波算法 在相应的机动和非机动阶段利用标准 Kalman 滤波算法进行滤波计算，具体算法公式如下 ˆ (k/k − 1) = ΦX ˆ (k − 1/k − 1), X P (k/k − 1) = ΦP (k − 1/k − 1)ΦT + ΓQΓT , K (k ) = P (k/k − 1)H T HP (k/k − 1)H T + R
2
目标模型
滑坡在开始非机动阶段经历的缓慢变形和匀速变形运动可用下面的差分方程描述。 1 x(k + 1) = x(k ) + T x ˙ (k ) + a(k )T 2 , 2 x ˙ (k + 1) = x ˙ (k ) + T a ˙ (k ).
收稿日期: 2007-09-19. 作者简介: 唐春艳 (1979年4月生)，女，硕士，讲师. 研究方向：数学与计算机模拟.
唐春艳1 , 彭继兵2
(1- 吉林大学珠海学院公共基础教研中心数学教研室，珠海 519041; 2- 成都理工大学信息管理学院，成都 610059) 摘 要: 本文研究如何对滑坡监测数据进行去噪。利用变维滤波算法能够很好描述滑坡变形的特性，建立 了滑坡运动的数学模型，并引入变维滤波对滑坡监测数据进行降噪。为了提高滑坡预报精度， 提出了一种基于变维滤波的 Kalman 最佳平滑算法。实例仿真结果表明利用该方法不仅预报精度 高，而且也为寻找滑坡监测数据中的变点提供了新的思路。 中图分类号: P642.22 文献标识码: A
(1)
第4期
唐春艳，彭继兵：基于变维滤波算法的 Kalman 最佳平滑器在滑坡监测数据处理中的应用
587
其中，x(k ) 和 x ˙ (k ) 分别表示滑坡在 k 时刻的位置和速度，a(k ) 表示滑坡在 k 到 k + 1 时刻由于 降雨等外界因素引起的加速度，它考虑了滑坡地质本身和外界环境扰动所造成的不可预测行 2 为。由于滑坡处于匀速变形的非机动阶段，故可设 a(k ) 是服从零均值、方差为 σa 的正态分布 2 平稳随机序列，且 a(k ) 和 a(l)(k = l) 互不相关，即 E {a(k )} = 0、E {a(k )a(l)} = σa δ (k − l)， 其中， 1, k = 0 δ (k ) = . 0, k = 0 由以上分析可建立滑坡在匀速变形非机动阶段的状态方程和观测方程 X (k + 1) = ΦX (k ) + ΓW (k ), X (k ) = x(k ) x ˙ (k ) , Z (k ) = HX (k ) + V (k ), 1 T T 2 /2 , Γ = , Φ= 0 T T (2) (3) W (k ) = a(k ),
取 ∆ = (1 − α)−1 作为检测机动的有效窗口长度，则机动检测方法为：如果 µ(k ) Th ，则认为 目标在 k − ∆ − 1 时开始有一恒定的加速度加入，这时目标模型应由低阶非机动模型转向高阶 机动模型，其式中 Th 为 nz 个自由度的 χ2 分布对应于给定置信区间的置信度。 由高阶机动模型退回到低阶非机动模型的检测方法是检验加速度估计值是否有统计显著性 意义。令
−1
(18) (19) , (20) (21) (22)
ˆ (k/k ) = X ˆ (k/k − 1) + K (k ) Z (k ) − H X ˆ (k/k − 1) , X P (k/k ) = I − K (k )H P (k/k − 1).
4
Kalman 最 佳 平 滑器
最佳平滑问题分为三种情况：最佳固定区间平滑、最佳固定点平滑、最佳固定滞后平 滑[5-6] ，本文主要采用最佳固定区间平滑，它对应的算法公式如下 ˆ (j | N ) = X ˆ (j | j ) + K (j ) X ˆ (j + 1 | N ) − ΦX ˆ (j | j ) , X K (j ) = P (j | j )ΦT P T (j + 1 | j ), P (j/N ) = P (j/j ) − K (j ) P (j + 1 | j ) − P (j + 1 | N ) K T (j ). (23) (24) (25)