顾及未知系统误差的变形监测滤波算法
左廷英，曾磊
(中南大学 地球科学与信息物理学院，湖南 长沙，410083)
摘要：通过建立一种带有未知系统误差的滤波模型，研究变形监测滤波解算中一些未知边坡物理信息的处理问题，给出一种基于移动窗口的系统误差自适应拟合法。研究结果表明：该算法充分融合了包含在观测中的统计信息、边坡滑坡的有关力学状态和地质条件所提供的信息，减弱了观测异常所带来的影响, 提高了形变参数解算精度。
关键词：系统误差；边坡；变形监测；卡尔曼滤波
中图分类号：P22 文献标志码：A 文章编号：1672-7207(2011)03-0738-06
Filter algorithm in consideration of unknown systematic errors for deformation monitoring
ZUO Ting-ying, ZENG Lei
(School of Earth Science and Info-physics, Central South University, Changsha 410083, China)
Abstract: By constructing a filter models with unknown systematic error, slope’s physical information disposal was studied in the course of deformation monitoring filter calculation, and an adaptive fittings algorithm of systematic errors was given using moving windows. The results show that the method is fully fused with statistic information in the course of observation, dynamical state of slope slippage and geophysical information, it can eliminate the influence of abnormal observation and improve the precision of deformation parameter estimation.
Key words: systematic errors; slope; deformation monitoring; Kalman filtering
滑坡的变形过程受多方面因素的影响[1-2]。人们对边坡的物理状态信息的认识和假设等，都会影响各测点的位移速率，使滑坡的位移时序曲线表现出波动性，从而使滑坡动态分析的难度增加，建立的模型无疑会包含一定程度的误差[3-8]，进而导致形变分析结果不精确。计算边坡形变的一种合理的方法是：将边坡物理模型估计的形变量与几何观测量相结合建立混合模型。这种混合模型最早由Schwintzer[9]提出；Bock 等[10-11]对此进行了进一步的研究，他们根据地球物理模型和实测拟合模型之间的差异, 应用抗差估计来调整先验参数对计算结果的影响。Segall等[12]提出“模型调整”法，即通过几何模型计算的位移和物理模型预测位移之间的差异达到最小的原则来求解形变量。在实际跟踪滑坡的动态变化过程中，总是存在一些未知边坡物理信息，它们通常只有微小的变化，在固定的观测历元间可以看成是常量，或围绕某常量随机变化，因此，可以采用固定窗口内的观测残差和状态预测残差对它们进行拟合[13]，并给出相应协方差矩阵的近似估计。这种估计方法的优点是：在Kalman滤波过程中，不仅能够减小因物理信息不充分而导致的系统偏差的影响，而且可以求出这些系统偏差，即从函数模型和随机模型2方面同时提高滤波结果的可靠性。
1 带有模型误差的滤波算法
在动态方程中，带有未知模型误差的滤波模型为：
(1)
其中：Xk为tk时刻的状态向量；Фk,k-1为状态转移矩阵；sk为未知的模型误差；Wk为动力模型噪声向量；Lk为观测向量；Hk为设计矩阵；ek为观测误差向量。
sk=0时的滤波模型为：
(2)
它的预测向量为：
的残差向量为：
式中：
为的残差向量。Lk的残差向量为：
。假设Wk和ek的协方差矩阵分别为
和
，且Wk，Wi，ek及ei互不相关。由于
(3)
在无误差情况下，
，故有：
(4)
显然，
，即
有偏。若已求得模型系统误差估值
，则修正后的动力模型预报值的误差方程为：
(5)
式中：
为在含有系统误差情形下的残差向量。它们与观测残差和预测残差及状态预测残差向量 不同。
实践中，如何求解sk的估计值仍需进一步研究。在固定的观测历元间视模型系统误差为常量，或围绕某常量随机变化，则可以采用固定窗口内的观测残差和状态预测残差进行拟合，最简单的方法是取平均值。根据前面的定义，sk为动力模型的系统误差，而动力模型系统误差应主要反映状态预报值的偏差。根据式(5)，在tk-i，应有：
i=1，…，N (6)
若假设动力学模型系统偏差在短时间
内维持微小变化，即满足
，则预测状态向量残差的期望
应该为0，即
。显然，将(6)式两端取和，再除以N，并考虑，有：
(7)
可以证明：由式(7)求得的为sk的无偏估计。因为已经假定由
表示的
含有系统误差，故理论上应有
，式(7)可改写为
(8)
式中：
和
分别为
和
的真误差，其数学期望为0。对式(8)取期望值得：
(9)
式(9)表明：由式(7)求得的动力模型系统误差是sk的无偏估计。
为求状态向量预报值的协方差矩阵
，可先求预报残差的近似协方差矩阵。因的期望为0，故近似地有
(10)
由于
(11)
取
，考虑式(10)有：
在一般情况下，希望直接求解
，则可取：
；i=1，…，N (12)
又由于
，考虑式(11)，有：
(13)
求得动力模型系统误差以及预测状态向量的协方差矩阵
后，可将代入观测方程和状态预测方程，并利用所求得的即可进行动态Kalman 滤波。
2 边坡滑坡的滤波模型
根据滑动面的类型，边坡滑坡可分为平面型、契型、曲面型和倾倒型等多种形式[14]。若把边坡滑坡视为平面问题，则它们都可用一个块体系统来描述(见图1), 即滑体可以视为由许多小的块体组成。本文中只考虑图1所示的块体系统。为了使问题简化，假设块体都是刚体。若块体的几何形状和力学性质都是已知的，则可以导出块体系统的运动方程。由牛顿第二定律可知，块体系统中的任何一块刚体的受力状态见图2，其运动可表示为：
(14a)
(14b)
其中：aix和aiz分别为第i块刚体在x方向和z方向的加速度；N为正压力；R为摩擦阻力；Nx和Nz分别为N在x和z方向的分量；Rx和Rz分别为R在x和z方向的分量；m为块体的质量；g为重力加速度；xi和zi为块体在x和z方向的位移；P和T为相邻块对它的作用力；t为时间。
方程 (14) 也可表示为：
图1 块体系统
Fig.1 Block system
图2 某一刚体的受力状态
Fig.2 Geometry and forces associated with a rigid block
(15)
其中：
对于整个块体系统，由式(15)可构成如下方程组：
(16)
M 和A1 由Mi和A1i 组成。
方程(16)表示的是1个刚体系统的运动方程，若系统中只有1个块体，则必有：
(17)
为了提高计算精度，对参数Y进行如下变换：Y=Y0+DY(其中，Y0是Y在边坡滑体处于极限平衡状态下的取值)。由于在极限平衡状态下，任意块体的加速度都为0，即axi=azi=0，故由式(16)可得：
(18)
G0和C0是G在极限平衡状态下的取值。由式(16)和(18)可得：
(19)
其中：
，为位移向量；
0, 
；
，为速度向量；a= 
。显然，DY 表示边坡滑体的当前状态与稳定状态之间的差别。边坡滑体的状态可用下列状态向量来描述：
按照刚体的运动方程，当刚体的运动从状态k 转移到k+1时，其位移和速度按下式变化：
(20)
除了雨后引起的地下水位变化和地震引起的震动外，在实际工作中，作用在边坡滑体上的外力一般是保持不变的；因此，可以假定DG=0及DYk+1=DYk。若外力产生变化，则只要这种变化很小，可以视为状态转移误差(系统噪声)，根据这一假设和方程(19)可得：
(21)
由式(20)和(21)可得状态转移方程：
(22)
模型(12)考虑了边坡运动的加速度，但在一般情况下，边坡在发生滑坡突变之前的滑移过程中，移动的加速度非常小，在绝大多数情况下，加速度如果能达1 cm/d2将预示滑坡将发生，而加速度1 cm/d2在运动学上则非常小，因此，可在模型(22)中删除加速度项，得到下列模型：
；
(23)
其中：
；
；
；Lk为观测向量；ek为观测误差。
由于边块的物理信息并不充分，不能得到sk的具体值，因此，在滤波过程中，把它看作未知输入信息。用Xk表示状态变量代替式(23)中的Xk′，由式(23)可以得到带有未知输入的变形监测滤波模型(见式(1))。
3 实例解算与分析
以湖南某高速公路边坡监测为例。该边坡已采用抗滑桩进行处治，边坡质量未知，在抗滑桩上布设观测点，观测其三维位移。采用GPS连续静态模式观测，采样时间为15 s，基线每3 h计算1次结果，共观测6月，约1 500个结果。
；
(其中：X和Y分别为x和y水平方向位移，Z为沉降)。由于监测点布置在抗滑桩上，抗滑桩主要向公路倾斜(即水平位移)，沉降较小。为了说明所提出的该方法的优越性，本文按照刚体的运动方程，当刚体的运动从状态k 转移到k+1时，其位移和速度按式(20)变化，设计了以下2种计算方案。
方案1 采用常见的卡尔曼滤波中最常见的动态方程，在这个动态方程里只考虑了刚体的运动，不考虑边坡的物理信息。
方案2 除了雨后引起的地下水位变化和地震引起的震动外，在实际工作中，认为作用在边坡滑体上的外力一般是保持不变的。因此，可以假定ΔG=0及ΔYk+1=ΔYk。若外力产生变化，则只要这种变化很小，可以视为状态转移误差(系统噪声)。
由于边块的物理信息并不充分，不能得到sk的具体值，因此，在滤波过程中，把它看作未知输入信息，采用本文给出的算法进行滤波解算。
图3和图4分别给出了方案1和方案2的计算结果。对比图3和图4可以看出：
(1) 受观测误差和基准点系统误差的双重影响，采用方案1求得的点位形变不仅不光滑, 精度也较低(见图3)。由于物理模型的误差具有系统性质, 由此求得的形变位移仍然存在较明显的系统误差。
(2) 在方案2中，因为增加了点位的物理信息，采用未知的物理信息进行拟合，部分地平衡了地球物理模型信息和观测信息对滑坡预测的贡献，精度明显提高(见图4)。
(3) 在边坡监测实践中，人们往往不能预知观测方程和动力学方程是否含有系统误差，于是，可采用移动窗口的系统误差进行拟合。
(4) 移动窗口的系统误差拟合在实践中也存在一定困难，因为预先很难确定窗口的宽度。在本次试算中，采用窗口宽度N=10。若系统误差变化区间有较大变化，则N也相应变化。如何选取窗口的宽度是实践中的难点，一般需通过多次试验经计算确定。
图3 方案1在X，Y和Z方向位移曲线
Fig.3 Displacement curve of Scheme 1 in X,Y and Z direction
图4 方案2在X，Y和Z方向位移曲线
Fig.4 Displacement curve of Scheme 2 in X,Y and Z direction
4 结论
(1) 在变形监测中，当边坡的物理信息不充分时，可以把它们看作一种未知的系统误差，利用滤波输出的观测残差和状态预报值残差对它们进行开窗拟合，这样就可以降低未知的系统偏差影响，从函数模型和随机模型2个方面同时改进滤波结果的可靠性。
(2) 本文提出的算法其应用前提是：在一个固定的时间窗口内，未知信息是1个常量。然而，在实际变形观测中，未知的信息会不断变化，因此，该解算方法还有待于进一步研究。
参考文献：
[1] Zangerl C, Eberhardt E, Perzlmaier S. Kinematic behaviour and velocity characteristics of a complex deep-seated crystalline rockslide system in relation to its interaction with a dam reservoir[J]. Engineering Geology, 2010, 112: 53-67.
[2] Bonzanigo L, Eberhardt E, Loew S. Long-term investigation of a deep-seated creeping landslide in crystalline rock-geological and hydromechanical factors controlling the campo vallemaggia landslide[J]. Canadian Geotechnical Journal, 2007, 44(10): 1157-1180.
[3] 朱建军, 丁晓利, 陈永奇. 集成地质、力学信息和监测数据的滑坡动态模型[J]. 测绘学报, 2003, 32(3): 261-266.
ZHU Jian-jun, DING Xiao-li, CHEN Yong-qi. Dynamic landsliding model with integration of monitoring information and mechanic information[J]. Acta Geodaetica et Cartographica Sinca, 2003, 32(3): 261-266.
[4] YANG Yuan-xi, ZENG An-min. Adaptive filtering for deformation parameter estimation in consideration of geometrical models[J]. Science in China Series D: Earth Science, 2009, 52(8): 1216-1222.
[5] Shi G H. Applications of discontinuous deformation analysis (DDA) to rock engineering[J]. Computtional Mechanics, 2009(2): 136-147.
[6] Bonaldi P. Displacement forecasting for concrete dams via deterministic mathematical models[J]. Water Power & Dam Construction, 1977, 29(9): 74-78.
[7] Purer E. Application of statistical methods in monitoring dam behavior[J]. Water Power & Dam Construction, 1986, 38(12): 16-19.
[8] DeSortis A, Paoliani P. Statistical analysis and structural identification in concrete dam monitoring[J]. Engineering Structures, 2007, 29: 110-120.
[9] Schwintzer P. Generalization for deformation vector with hybrid model[C]//Jo? I, Detrek?i A, eds. Deformation Measurements. Budapest: Akademiai kiadó, 1982: 453-463.
[10] Bock Y. Estimating crustal deformations from a combination of baseline measurements and geophysical models[J]. Bull Geod, 1983, 57: 294-311.
[11] Bock Y, Schaffrin B. Robust predication of the Earth’s crustal movements from precise geodetic data and a vague geophysical mode[C]//The First World Congress of Bernoulli Society on Mathematical Statistics. Taschkent (USSR), 1986 361-367.
[12] Segall P, Matthews M V. Displacement calculations from geodetic data and the testing of geophysical deformation model[J]. J Geophys Res, 1988, 93(B12): 14954-14966.
[13] 杨元喜, 张双成. 导航解算中的系统误差及其协方差矩阵拟合[J]. 测绘学报, 2004, 33(3): 189-194.
YANG Yuan-xi, ZHANG Shuang-cheng. Fittings of systematic errors and covariance matrices in navigation[J]. Acta Geodaetica et Cartographica Sinca, 2004, 33(3): 189-194.
[14] Lü W C, Xu S Q. Kalman filtering algorithm research for the deformation information series of the similar single difference model[J]. Journal of China University of Mining and Technology, 2004, 14(2): 189-194.
(编辑 陈灿华)
收稿日期：2010-05-10；修回日期：2010-07-25
基金项目：国家自然科学基金资助项目(40874005)
通信作者：左廷英(1964-)，女，山西大同人，博士研究生，副教授，从事测量数据处理等研究；电话：0731-88837860；E-mail：ZTY2003@163.com
