左廷英，曾 磊

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

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

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

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

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

sk=0時(shí)的濾波模型為：

它的預(yù)測向量為：

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

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

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

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

2 邊坡滑坡的濾波模型

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

其中：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為相鄰塊對(duì)它的作用力；t為時(shí)間。

方程 (14) 也可表示為：

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

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

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

M 和A1由Mi和A1i組成。

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

為了提高計(jì)算精度，對(duì)參數(shù) Y進(jìn)行如下變換：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 表示邊坡滑體的當(dāng)前狀態(tài)與穩(wěn)定狀態(tài)之間的差別。邊坡滑體的狀態(tài)可用下列狀態(tài)向量來描述：

按照剛體的運(yùn)動(dòng)方程，當(dāng)剛體的運(yùn)動(dòng)從狀態(tài)k 轉(zhuǎn)移到k+1時(shí)，其位移和速度按下式變化：

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

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

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

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

3 實(shí)例解算與分析

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

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

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

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

圖3和圖4分別給出了方案1和方案2的計(jì)算結(jié)果。對(duì)比圖3和圖4可以看出：

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

(2) 在方案 2中，因?yàn)樵黾恿它c(diǎn)位的物理信息，采用未知的物理信息進(jìn)行擬合，部分地平衡了地球物理模型信息和觀測信息對(duì)滑坡預(yù)測的貢獻(xiàn)，精度明顯提高(見圖4)。

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

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

圖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 結(jié)論

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

(2) 本文提出的算法其應(yīng)用前提是：在一個(gè)固定的時(shí)間窗口內(nèi)，未知信息是1個(gè)常量。然而，在實(shí)際變形觀測中，未知的信息會(huì)不斷變化，因此，該解算方法還有待于進(jìn)一步研究。

[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ì)、力學(xué)信息和監(jiān)測數(shù)據(jù)的滑坡動(dòng)態(tài)模型[J]. 測繪學(xué)報(bào), 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] 楊元喜, 張雙成. 導(dǎo)航解算中的系統(tǒng)誤差及其協(xié)方差矩陣擬合[J]. 測繪學(xué)報(bào), 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.