顾及未知系统误差的变形监测滤波算法

左廷英，曾 磊

(中南大学 地球科学与信息物理学院，湖南 长沙，410083)

滑坡的变形过程受多方面因素的影响[1-2]。人们对边坡的物理状态信息的认识和假设等，都会影响各测点的位移速率，使滑坡的位移时序曲线表现出波动性，从而使滑坡动态分析的难度增加，建立的模型无疑会包含一定程度的误差[3-8]，进而导致形变分析结果不精确。计算边坡形变的一种合理的方法是：将边坡物理模型估计的形变量与几何观测量相结合建立混合模型。这种混合模型最早由 Schwintzer[9]提出；Bock等[10-11]对此进行了进一步的研究，他们根据地球物理模型和实测拟合模型之间的差异, 应用抗差估计来调整先验参数对计算结果的影响。Segall等[12]提出“模型调整”法，即通过几何模型计算的位移和物理模型预测位移之间的差异达到最小的原则来求解形变量。在实际跟踪滑坡的动态变化过程中，总是存在一些未知边坡物理信息，它们通常只有微小的变化，在固定的观测历元间可以看成是常量，或围绕某常量随机变化，因此，可以采用固定窗口内的观测残差和状态预测残差对它们进行拟合[13]，并给出相应协方差矩阵的近似估计。这种估计方法的优点是：在 Kalman滤波过程中，不仅能够减小因物理信息不充分而导致的系统偏差的影响，而且可以求出这些系统偏差，即从函数模型和随机模型2方面同时提高滤波结果的可靠性。

1 带有模型误差的滤波算法

在动态方程中，带有未知模型误差的滤波模型为：

其中：Xk为 tk时刻的状态向量；Фk,k-1为状态转移矩阵；sk为未知的模型误差；Wk为动力模型噪声向量；Lk为观测向量；Hk为设计矩阵；ek为观测误差向量。

sk=0时的滤波模型为：

它的预测向量为：

实践中，如何求解sk的估计值仍需进一步研究。在固定的观测历元间视模型系统误差为常量，或围绕某常量随机变化，则可以采用固定窗口内的观测残差和状态预测残差进行拟合，最简单的方法是取平均值。根据前面的定义，sk为动力模型的系统误差，而动力模型系统误差应主要反映状态预报值的偏差。根据式(5)，在tk-i，应有：

若假设动力学模型系统偏差在短时间(tk-N～tk)内维持微小变化，即满足 E ( Wk-i) =sk，则预测状态向量残差的期望应该为0，即) = 0。显然，将(6)式两端取和，再除以N，并考虑，有：

可以证明：由式(7)求得的sˆk为sk的无偏估计。因为已经假定由表示的含有系统误差，故理论上应有，式(7)可改写为

式(9)表明：由式(7)求得的动力模型系统误差sˆk是sk的无偏估计。

2 边坡滑坡的滤波模型

根据滑动面的类型，边坡滑坡可分为平面型、契型、曲面型和倾倒型等多种形式[14]。若把边坡滑坡视为平面问题，则它们都可用一个块体系统来描述(见图1), 即滑体可以视为由许多小的块体组成。本文中只考虑图1所示的块体系统。为了使问题简化，假设块体都是刚体。若块体的几何形状和力学性质都是已知的，则可以导出块体系统的运动方程。由牛顿第二定律可知，块体系统中的任何一块刚体的受力状态见图2，其运动可表示为：

其中：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)可构成如下方程组：

M 和A1由Mi和A1i组成。

方程(16)表示的是 1个刚体系统的运动方程，若系统中只有1个块体，则必有：

为了提高计算精度，对参数 Y进行如下变换：Y=Y0+ΔY(其中，Y0是 Y在边坡滑体处于极限平衡状态下的取值)。由于在极限平衡状态下，任意块体的加速度都为0，即axi=azi=0，故由式(16)可得：

G0和C0是G在极限平衡状态下的取值。由式(16)和(18)可得：

其中：Δ Y = Y -Y0; X1=(x1, z1, … ,xi, zi…)T，为位移向量； Δ G = G - G0= ( 0, Δ m1g , … , 0, Δ mig ,…)；V = X˙1=… ,vxi, vzi,…)T，为速度向量；a=V˙= ( ax1, az1,… ,axi, azi,…)T。显然，ΔY 表示边坡滑体的当前状态与稳定状态之间的差别。边坡滑体的状态可用下列状态向量来描述：

按照刚体的运动方程，当刚体的运动从状态k 转移到k+1时，其位移和速度按下式变化：

除了雨后引起的地下水位变化和地震引起的震动外，在实际工作中，作用在边坡滑体上的外力一般是保持不变的；因此，可以假定ΔG=0及ΔYk+1=ΔYk。若外力产生变化，则只要这种变化很小，可以视为状态转移误差(系统噪声)，根据这一假设和方程(19)可得：

由式(20)和(21)可得状态转移方程：

模型(12)考虑了边坡运动的加速度，但在一般情况下，边坡在发生滑坡突变之前的滑移过程中，移动的加速度非常小，在绝大多数情况下，加速度如果能达1 cm/d2将预示滑坡将发生，而加速度1 cm/d2在运动学上则非常小，因此，可在模型(22)中删除加速度项，得到下列模型：

由于边块的物理信息并不充分，不能得到sk的具体值，因此，在滤波过程中，把它看作未知输入信息。用Xk表示状态变量代替式(23)中的Xk′，由式(23)可以得到带有未知输入的变形监测滤波模型(见式(1))。

3 实例解算与分析

以湖南某高速公路边坡监测为例。该边坡已采用抗滑桩进行处治，边坡质量未知，在抗滑桩上布设观测点，观测其三维位移。采用GPS连续静态模式观测，采样时间为15 s，基线每3 h计算1次结果，共观测6月，约1 500个结果。 X1,k+1=(Xk, Yk, Zk)T；V1,k+1=(其中：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.