左廷英，曾 磊

(中南大學 地球科學與信息物理學院，湖南 長沙，410083)

滑坡的變形過程受多方面因素的影響[1-2]。人們對邊坡的物理狀態(tài)信息的認識和假設等，都會影響各測點的位移速率，使滑坡的位移時序曲線表現(xiàn)出波動性，從而使滑坡動態(tài)分析的難度增加，建立的模型無疑會包含一定程度的誤差[3-8]，進而導致形變分析結果不精確。計算邊坡形變的一種合理的方法是：將邊坡物理模型估計的形變量與幾何觀測量相結合建立混合模型。這種混合模型最早由 Schwintzer[9]提出；Bock等[10-11]對此進行了進一步的研究，他們根據(jù)地球物理模型和實測擬合模型之間的差異, 應用抗差估計來調(diào)整先驗參數(shù)對計算結果的影響。Segall等[12]提出“模型調(diào)整”法，即通過幾何模型計算的位移和物理模型預測位移之間的差異達到最小的原則來求解形變量。在實際跟蹤滑坡的動態(tài)變化過程中，總是存在一些未知邊坡物理信息，它們通常只有微小的變化，在固定的觀測歷元間可以看成是常量，或圍繞某常量隨機變化，因此，可以采用固定窗口內(nèi)的觀測殘差和狀態(tài)預測殘差對它們進行擬合[13]，并給出相應協(xié)方差矩陣的近似估計。這種估計方法的優(yōu)點是：在 Kalman濾波過程中，不僅能夠減小因物理信息不充分而導致的系統(tǒng)偏差的影響，而且可以求出這些系統(tǒng)偏差，即從函數(shù)模型和隨機模型2方面同時提高濾波結果的可靠性。

1 帶有模型誤差的濾波算法

在動態(tài)方程中，帶有未知模型誤差的濾波模型為：

其中：Xk為 tk時刻的狀態(tài)向量；Фk,k-1為狀態(tài)轉移矩陣；sk為未知的模型誤差；Wk為動力模型噪聲向量；Lk為觀測向量；Hk為設計矩陣；ek為觀測誤差向量。

sk=0時的濾波模型為：

它的預測向量為：

實踐中，如何求解sk的估計值仍需進一步研究。在固定的觀測歷元間視模型系統(tǒng)誤差為常量，或圍繞某常量隨機變化，則可以采用固定窗口內(nèi)的觀測殘差和狀態(tài)預測殘差進行擬合，最簡單的方法是取平均值。根據(jù)前面的定義，sk為動力模型的系統(tǒng)誤差，而動力模型系統(tǒng)誤差應主要反映狀態(tài)預報值的偏差。根據(jù)式(5)，在tk-i，應有：

若假設動力學模型系統(tǒng)偏差在短時間(tk-N～tk)內(nèi)維持微小變化，即滿足 E ( Wk-i) =sk，則預測狀態(tài)向量殘差的期望應該為0，即) = 0。顯然，將(6)式兩端取和，再除以N，并考慮，有：

可以證明：由式(7)求得的s?k為sk的無偏估計。因為已經(jīng)假定由表示的含有系統(tǒng)誤差，故理論上應有，式(7)可改寫為

式(9)表明：由式(7)求得的動力模型系統(tǒng)誤差s?k是sk的無偏估計。

2 邊坡滑坡的濾波模型

根據(jù)滑動面的類型，邊坡滑坡可分為平面型、契型、曲面型和傾倒型等多種形式[14]。若把邊坡滑坡視為平面問題，則它們都可用一個塊體系統(tǒng)來描述(見圖1), 即滑體可以視為由許多小的塊體組成。本文中只考慮圖1所示的塊體系統(tǒng)。為了使問題簡化，假設塊體都是剛體。若塊體的幾何形狀和力學性質(zhì)都是已知的，則可以導出塊體系統(tǒng)的運動方程。由牛頓第二定律可知，塊體系統(tǒng)中的任何一塊剛體的受力狀態(tài)見圖2，其運動可表示為：

其中：aix和aiz分別為第i塊剛體在x方向和z方向的加速度；N為正壓力；R為摩擦阻力；Nx和Nz分別為N在x和z方向的分量；Rx和Rz分別為R在x和z方向的分量；m為塊體的質(zhì)量；g為重力加速度；xi和zi為塊體在x和z方向的位移；P和T為相鄰塊對它的作用力；t為時間。

方程 (14) 也可表示為：

圖1 塊體系統(tǒng)Fig.1 Block system

圖2 某一剛體的受力狀態(tài)Fig.2 Geometry and forces associated with a rigid block

對于整個塊體系統(tǒng)，由式(15)可構成如下方程組：

M 和A1由Mi和A1i組成。

方程(16)表示的是 1個剛體系統(tǒng)的運動方程，若系統(tǒng)中只有1個塊體，則必有：

為了提高計算精度，對參數(shù) Y進行如下變換：Y=Y0+ΔY(其中，Y0是 Y在邊坡滑體處于極限平衡狀態(tài)下的取值)。由于在極限平衡狀態(tài)下，任意塊體的加速度都為0，即axi=azi=0，故由式(16)可得：

G0和C0是G在極限平衡狀態(tài)下的取值。由式(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 表示邊坡滑體的當前狀態(tài)與穩(wěn)定狀態(tài)之間的差別。邊坡滑體的狀態(tài)可用下列狀態(tài)向量來描述：

按照剛體的運動方程，當剛體的運動從狀態(tài)k 轉移到k+1時，其位移和速度按下式變化：

除了雨后引起的地下水位變化和地震引起的震動外，在實際工作中，作用在邊坡滑體上的外力一般是保持不變的；因此，可以假定ΔG=0及ΔYk+1=ΔYk。若外力產(chǎn)生變化，則只要這種變化很小，可以視為狀態(tài)轉移誤差(系統(tǒng)噪聲)，根據(jù)這一假設和方程(19)可得：

由式(20)和(21)可得狀態(tài)轉移方程：

模型(12)考慮了邊坡運動的加速度，但在一般情況下，邊坡在發(fā)生滑坡突變之前的滑移過程中，移動的加速度非常小，在絕大多數(shù)情況下，加速度如果能達1 cm/d2將預示滑坡將發(fā)生，而加速度1 cm/d2在運動學上則非常小，因此，可在模型(22)中刪除加速度項，得到下列模型：

由于邊塊的物理信息并不充分，不能得到sk的具體值，因此，在濾波過程中，把它看作未知輸入信息。用Xk表示狀態(tài)變量代替式(23)中的Xk′，由式(23)可以得到帶有未知輸入的變形監(jiān)測濾波模型(見式(1))。

3 實例解算與分析

以湖南某高速公路邊坡監(jiān)測為例。該邊坡已采用抗滑樁進行處治，邊坡質(zhì)量未知，在抗滑樁上布設觀測點，觀測其三維位移。采用GPS連續(xù)靜態(tài)模式觀測，采樣時間為15 s，基線每3 h計算1次結果，共觀測6月，約1 500個結果。 X1,k+1=(Xk, Yk, Zk)T；V1,k+1=(其中：X和Y分別為x和y水平方向位移，Z為沉降)。由于監(jiān)測點布置在抗滑樁上，抗滑樁主要向公路傾斜(即水平位移)，沉降較小。為了說明所提出的該方法的優(yōu)越性，本文按照剛體的運動方程，當剛體的運動從狀態(tài)k 轉移到k+1時，其位移和速度按式(20)變化，設計了以下2種計算方案。

方案 1 采用常見的卡爾曼濾波中最常見的動態(tài)方程，在這個動態(tài)方程里只考慮了剛體的運動，不考慮邊坡的物理信息。

方案 2 除了雨后引起的地下水位變化和地震引起的震動外，在實際工作中，認為作用在邊坡滑體上的外力一般是保持不變的。因此，可以假定ΔG=0及ΔYk+1=ΔYk。若外力產(chǎn)生變化，則只要這種變化很小，可以視為狀態(tài)轉移誤差(系統(tǒng)噪聲)。

由于邊塊的物理信息并不充分，不能得到sk的具體值，因此，在濾波過程中，把它看作未知輸入信息，采用本文給出的算法進行濾波解算。

圖3和圖4分別給出了方案1和方案2的計算結果。對比圖3和圖4可以看出：

(1) 受觀測誤差和基準點系統(tǒng)誤差的雙重影響，采用方案1求得的點位形變不僅不光滑, 精度也較低(見圖3)。由于物理模型的誤差具有系統(tǒng)性質(zhì), 由此求得的形變位移仍然存在較明顯的系統(tǒng)誤差。

(2) 在方案 2中，因為增加了點位的物理信息，采用未知的物理信息進行擬合，部分地平衡了地球物理模型信息和觀測信息對滑坡預測的貢獻，精度明顯提高(見圖4)。

(3) 在邊坡監(jiān)測實踐中，人們往往不能預知觀測方程和動力學方程是否含有系統(tǒng)誤差，于是，可采用移動窗口的系統(tǒng)誤差進行擬合。

(4) 移動窗口的系統(tǒng)誤差擬合在實踐中也存在一定困難，因為預先很難確定窗口的寬度。在本次試算中，采用窗口寬度N=10。若系統(tǒng)誤差變化區(qū)間有較大變化，則N也相應變化。如何選取窗口的寬度是實踐中的難點，一般需通過多次試驗經(jīng)計算確定。

圖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) 在變形監(jiān)測中，當邊坡的物理信息不充分時，可以把它們看作一種未知的系統(tǒng)誤差，利用濾波輸出的觀測殘差和狀態(tài)預報值殘差對它們進行開窗擬合，這樣就可以降低未知的系統(tǒng)偏差影響，從函數(shù)模型和隨機模型2個方面同時改進濾波結果的可靠性。

(2) 本文提出的算法其應用前提是：在一個固定的時間窗口內(nèi)，未知信息是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] 朱建軍, 丁曉利, 陳永奇. 集成地質(zhì)、力學信息和監(jiān)測數(shù)據(jù)的滑坡動態(tài)模型[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] 楊元喜, 張雙成. 導航解算中的系統(tǒng)誤差及其協(xié)方差矩陣擬合[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.