滑坡監測的自適應約束抗差濾波算法研究*

左廷英 吳蕓蕓 宋迎春

(中南大學信息物理工程學院，長沙 410083)

滑坡監測的自適應約束抗差濾波算法研究*

左廷英 吳蕓蕓 宋迎春

(中南大學信息物理工程學院，長沙 410083)

在滑坡的動態監測中通常存在先驗幾何信息和物理信息，現有的濾波算法并不能夠充分利用這些信息，為控制幾何觀測異常對形變參數估計的影響，提出一種利用滑坡力學狀態信息的約束濾波模型。同時，利用約束方程獲得有益信息，精確地確定測量的等價權，給出相應的自適應約束抗差濾波算法，并通過實例說明了帶有約束的濾波算法更能提高狀態參數估計的可靠性。

約束信息;變形監測;抗差濾波;約束濾波;自適應濾波

1 引言

邊坡巖體由初始變形發展到破壞性滑坡，是一個極其復雜的發展演變過程［1，2］，國內外學者已建立了許多動態的滑坡預測預報模型［3-8］。然而，由于影響滑坡的因素非常復雜，而人們的認識又十分有限，因此，滑坡動態分析的難度通常很大，建立的模型常常包含一定程度的誤差，從而導致不精確的形變分析結果。在變形觀測中，大地測量方面的學者只關注幾何觀測數據，對形變的物理背景知識了解甚少，當在觀測過程出現了異常數據容易造成計算結果的偏差。所以，在建立邊坡濾波模型計算邊坡形變時，一種合理的方法是將物理力學模型與幾何觀測模型相結合。Schwintzer［9］于1982年提出了這種建模思想，Bock等［10，11］對這一建模思想進行了進一步的拓展，他們根據地球物理模型和實測幾何數據建立變形監測模型，應用抗差估計理論來調整先驗信息對計算結果的影響。Segall和Matthews［12］建立了一個幾何觀測位移量和物理模型預測位移量間的差異達到最小的原則，用來解算變形體的形變量。為了精確地跟蹤滑坡的動態變化，我們提出一種新的方法，即把邊坡的幾何約束信息、力學狀態結合起來，建立一個混合的濾波模型，合理地利用先驗的幾何信息和物理信息建立關于狀態參數的約束方程，再利用這些約束信息來控制幾何觀測異常對形變參數估計的影響。

2 滑坡的濾波模型

依據滑動面的類型不同，滑坡類型可分為平面型、契型、曲面型和傾倒型等多種形式［13］。若把滑坡視為一個平面問題，那么可以用一個塊體系統來描述(圖1)［14］。為了讓問題變得簡單，假設塊體為剛體，并且塊體的幾何形狀和力學性質都是預先知道的，由牛頓第二定律可知塊體系統中的任何一塊的運動方程可表示為(圖2)。

圖1 塊體系統Fig.1 Block system

圖2 剛體受力狀態示意圖Fig.2 Geometry and forces associated with a rigid block

式中，aix、aiz表示第i塊滑塊的加速度，N表示滑塊的正壓力，R表示滑塊所受的摩擦阻力，Nx、Nz、Rx、Rz為N和R中的在x和z方向受力的分量，m為塊體的質量，g為重力加速度，P和T是相鄰塊對它的作用力，xi、zi是塊體的位移量。

方程(1)、(2)也可表示為

這里

對于整個塊體系統，由式(3)可構成

這里，M和A1由Mi和A1i組成，

方程(3)表示的是一個剛體系統的運動方程，如果系統中只有一個塊體，必有：

假設塊體為剛體，并且各塊體面保持相互接觸，這時滑體將沿滑動面向下運動，在滑動面的法方向上的運動顯然為0，即

如果系統中有n塊滑體，則必有n個類似式(6)的方程。

此外，由于各塊體面之間保持相互接觸，相鄰塊體的位移、速度和加速度在垂直于接觸面的法線方向上都是相同的，即對于加速度有：

如果有n個滑動塊體，則必有n-1個類似式(7)的方程。對整個滑體系統，方程(6)、(7)可以組合成：

A2按(6)、(7)組成。

摩擦力為：

其中，k表示摩擦系數，c表示內聚力，對n個塊體的塊體系統，就有2n-1個類似式(9)的方程。所有這些方程寫成

為了方便計算，我們對參數Y進行如下變換，即讓Y=Y0+ΔY，這里Y0是Y在滑體處于極限平衡狀況下的取值。因為，在極限平衡狀態下，任意塊體的加速度值應該為0，所以有axi=azi，由式(4)可得：

式中，G0、C0是G在極限平衡狀態下的取值。由式(4)、(8)和(11)可得：

這里，ΔY=Y-Y0，ΔG=G-G0=(0，Δm1g，…，0，Δmig，…)，顯然ΔY是表示滑體當前狀態與穩定狀態之間的差別。

X1是位移向量，V=是速度向量，滑體的狀態可用X=(，，aT，ΔYT)T狀態向量來描述。其中，X1=(x1，z1，…，xi，zi，…)T，V==(vx1，vz1，…，vxi，vzi，…)T，a==(ax1，az1，…，axi，azi，…)T。

按照剛體運動方程可以得到，剛體的運動從狀態k轉移到k+1時，其位移和速度有：

除一些特殊原因外，作用在滑體上的外力通常沒有大的變化。因此，可以假定ΔG=0及ΔYk+1= ΔYk。當外力發生微小變化時，則可以把它們看成中狀態轉移誤差(系統噪聲)。因此，由方程(12)可得：

由式(13)和(14)可得狀態轉移方程為

在邊坡監測中觀測方程可表示為

這里，Hk=(I 0)，Lk表示觀測向量，ek表示觀測誤差。

在剛體假設下我們可以建立某些約束條件，例如由式(12)，有：

式(17)可以用作為幾何約束。此外，相鄰兩塊體的位移和速度接觸面的法向上應該相同，因此，類似于式(10)，有：

式(17)和(18)可合并為DkXk=dk，這里

由式(15)～(17)得到的變形監測的濾波模型為：

其中，Xk為tk時刻的狀態向量，Φk為狀態轉移矩陣，Wk為動態噪聲向量。式中Lk為觀測向量，Hk為設計矩陣，ek為觀測噪聲向量，ek和Wk不相關。

3 自適應約束濾波算法

考慮測量誤差方程Vk=Hk-Lk，則帶約束的狀態估計為［17］：

容易發現，式(22)的第一項與自適應卡爾曼濾波估計類似。設

估計式(25)幾乎和帶等式約束的標準卡爾曼濾波遞推公式相同，自適應因子通過一步向前估計式(23)影響狀態估計，假若和被看成是先驗信息，那么自適應等式約束卡爾曼濾波是非常類似標準卡爾曼濾波。

4 自適應因子的確定

4.1 狀態不符值統計量

一個統一的自適應因子能夠方便地調節測量信息對預測狀態估計的貢獻，簡化實時處理的計算。文獻［18］給出的由狀態不符值所構建的自適應因子為：

4.2 預測殘差統計量

當tk歷元觀測信息Lk可靠時，其相應的預測殘差=Hk-Lk，預測殘差的理論協方差陣為=HkHk+Σk。反映了動力學模型的誤差大小，利用它可以構造自適應因子［18］：

5 實例解算與分析

以湖南某高速公路邊坡監測為例，該邊坡已采用抗滑樁進行處治，邊坡重量大約10噸，實測α1= 45°。在抗滑樁上布設觀測點，觀測三維位移。采用兩臺GPS接收機連續靜態模式觀測，采樣率為15 s，基線每3小時計算一次結果，共觀測了6個月。由于觀測期間有一臺GPS接收機信號缺失，無法獲得正常GPS觀測值，導致某些時間段的觀測數據出現了一些異常，因此位移的計算結果脫離實際，本例將通過先驗約束信息對其中一個邊坡觀測數據進行校正，濾波結果與正常接收機的計算進行比較。設計了3種計算方案：

方案1：采用常見的卡爾曼濾波計算。

方案2：采用無約束的自適應抗差濾波式(23)計算，自適應因子按照式(26)計算。

方案3：采用約束自適應抗差濾波式(25)計算，自適應因子按照式(26)計算。

圖3～5給出了3種方案的計算結果(每個歷元的結果是按照3小時的靜態計算為一個結果)，從分析計算結果可以看出：

1)由于接收機信號缺失導致觀測誤差異常，利用方案1的標準卡爾曼濾波進行狀態參數估計求得的形變位移參數，與用正常觀測數據求得的形變位移參數相比較，存在較明顯的偏差(圖3)，特別是在歷元800～1 000附近出現了較大的抖動，這說明用標準的卡爾曼濾波算法不能有效地校正因觀測異常導致的濾波結果的偏差。

2)方案2中，因為采用了無約束的自適應抗差濾波，其形變參數的計算精度均有很大的改善(表1)，說明自適應因子起到了調控作用，它能很好地平衡物理力學信息和觀測信息對形變參數估計的貢獻，然而由于觀測本身并不可靠，要完全修正觀測結果是比較困難的。從圖5可見，對于歷元800～1000附近出現的抖動仍然無法修復。

表13 個方向的RMS表(單位：mm)Tab.1 RMS of three components(unit：mm)

圖3 方案1的X、Y、Z方向誤差曲線圖Fig.3 Displacement curves in X，Y，Z of Scheme 1

3)方案3中，由于增加了約束信息，當觀測不足情況下，起到了補充觀測信息的作用，因此，能夠修補因為長時間觀測異常引起的結果偏差(圖5)，說明用帶有狀態約束的抗差自適應濾波算法可以有效地校正因觀測異常導致濾波結果的偏差。

圖4 方案2的X、Y、Z方向位移曲線圖Fig.4 Displacement curves in X，Y，Z of Scheme 2

圖5 方案3的X、Y、Z方向位移曲線圖Fig.5 Displacement curves in X，Y，Z of Scheme 3

6 小結

模型的不確定性或模型誤差會嚴重影響濾波的分析結果，甚至會導致濾波的發散。本文采用附加先驗約束信息來消除或削弱模型誤差的影響。先驗約束信息常存在于邊坡的動態監測模型中，但是現有的濾波算法并不能夠充分利用這一信息，針對這一點，提出了一種狀態變量帶有約束的抗差濾波算法，由于約束信息能夠改善量測殘差，使得濾波解算過程可以從約束方程中獲得有益信息，因此，帶有狀態約束的抗差濾波算法能夠改善狀態估計的精度。可靠的殘差可以導致精確地確定測量的等價權，因而可以改善濾波解的可靠性。約束濾波遞推解與無約束的卡爾曼濾波遞推解非常類似，無約束的狀態估計可以用來評價狀態估計，然后約束方程再來對狀態估計進行更新。在抗差估計中，無約束濾波過程首先提供一個初始的狀態估計解，然而再綜合利用約束信息進行更新。實例解算說明，抗差自適應濾波和帶約束的抗差自適應濾波都能夠改善狀態參數的精度。然而，考慮邊坡的幾何約束信息、力學狀態，將變形監測數據通過卡爾曼濾波就可以結合起來，更能夠加強濾波解算結果的可靠性，有效消除或削弱模型誤差的影響。

1 Zangerl C，Eberhardt E and 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 and Loew S.Long-term investigation of a deep-seated creeping landslide in crystalline rockgeological and hydromechanical factors controlling the Campo Vallemaggia landslide［J］.Canadian Geotechnical Journal，2007，44(10)：1 157-1 180.

3 朱建軍，丁曉利，陳永奇.集成地質、力學信息和監測數據的滑坡動態模型［J］.測繪學報，2003，32(3)：261-266.(Zhu Jianjun，Ding Xiaoli and Chen Yongqi.Dynamic landsliding model with integration of monitoring information and mechanic information［J］.Acta Geodaetica et Cartographica Sinca，2003，32(3)：261-266)

4 Yang Yuanxi and Zeng Anmin.Adaptive filtering for deformation parameter estimation in consideration of geometrical models［J］.Science in China Series D：Earth Science，2009，52(8)：1 216-1 222.

5 Shi G H.Block system modeling by discontinuous deformation analysis［M］.Computational Mechanics Publications，Southampton UK and Boston USA，1993.

6 Bonaldi D.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［A］.In：I Joó，A.Detrek?i，eds.Deformation Measurements［C］.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［R］.The first World Congress of BernoulliSociety on MathematicalStatistics. Taschkent (USSR)，1986.

12 Segall P and Matthews M V.Displacement calculations from geodetic data and the testing of geophysical deformation model［J］.J Geophys Res.，1988，93(B12)：14 954-14 966.

14 Lü W C and 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.

15 Hoek E and Bray J W.Rock slope engineering［R］.The Institution of Mining and Metallurgy，London，1981.

16 Ren D and Ding X.Dynamic deformation analysis of open pit slopes［R］.The 8th Fig International Symposium on Deformation Measurements，25-28 June 1996，Hongkong：157-163.

17 Yang Y，He H and Xu G.Adaptively robust filtering for kinematic geodetic positioning［J］.Journal of Geodesy，2001，75(2/3)：109-116.

18 Yang Y，Gao W and Zhang X.Robust Kalman filtering with constraints：a case study for integrated navigation［J］.Journal of Geodesy，2010，84：373-381.

19 Yang Y and Gao W.An optimal adaptive Kalman filter［J］.Journal of Geodesy，2006，80： 177-183.

ADAPTIVELY ROBUSTLY CONSTRAINED FILTERING ALGORITHM FOR SLOPE MONITORING

Zuo Tingying，Wu Yunyun and Song Yingchun
(School of Info.Physics and Geomatics Engineering，Central South University，Changsha 410083)

For the rational use of geometric and physical information of slope，and in order to control the influence of outliers on deformation parameter estimation，a constrained filtering model using mechanical state information of slope is given.At the same time，the algorithm utilizes the useful information in constrained equation to accurately determine the equivalent weights of measurements，thus the relevant adaptive robust constrained filtering algorithm is given.By a pratical example，it is proved that this algorithm can improve the reliability of state parameter estimation.

constrained information;deformation monitoring;robust filtering;constrained filtering;adaptive filtering

1671-5942(2011)06-0094-06

2011-04-23

國家自然科學基金(40874005);教育部博士點基金(200805331086)

左廷英，女，1964年生，副教授，研究方向：測量數據處理.E-mail：zty2003@163.com

P207

A