以速率為參量的GM(1,1)滑坡時間預報模型研究

，謝金

(廈門大學 a.建筑與土木工程學院, 福建 廈門 361005； b.深圳研究院， 廣東 深圳 518057)

1 研究背景

滑坡是世界上比較活躍的自然災害之一。一旦達到一定規(guī)模的滑體沿著軟弱結(jié)構(gòu)面下滑時，就有可能對其影響范圍內(nèi)的建筑物(或構(gòu)筑物)、人員或財產(chǎn)造成嚴重的損害或損失。因此，為了科學地防治滑坡，減少傷亡和損失，滑坡的預測預報顯得更為重要。

自20世紀60年代以來，國內(nèi)外學者先后提出了30多種滑坡預測預報的方法和模型。這些方法和模型大致可分為5類：①確定性預報模型,如齋藤模型[1]、福囿斜坡時間預報法[2]等；②統(tǒng)計預報模型,其代表性的模型有灰色GM(1,1)模型[3]、Verhulst模型[4-5]、等；③非線性預報模型,其代表性的模型有突變理論模型[6]、協(xié)同預報模型[7]、BP神經(jīng)網(wǎng)絡(luò)模型[8]等；④宏觀預報模型,該模型是根據(jù)滑坡演化過程中表現(xiàn)出的各種征兆、跡象，如裂縫寬度，滑帶土的孔隙率等建立各類滑坡特征相適應的預報模型，模型因滑坡個體差異而不同。⑤空間預報模型,如邏輯回歸模型[9-11]、隨機森林模型[12]。

從摩擦學的角度看，滑坡演化的本質(zhì)是構(gòu)成滑坡的顆粒或者顆粒團的內(nèi)摩擦行為。當用電子顯微鏡觀察滑帶土的特征時，可以發(fā)現(xiàn)滑帶土中有許多酥裂狀結(jié)構(gòu)，切礫微斷裂，以及放射狀或不規(guī)則狀裂紋；同時也以觀察到碎礫定向排列和顆粒壓扁拉長等變形現(xiàn)象，宏觀上表現(xiàn)為位移監(jiān)測曲線呈現(xiàn)波動上升的趨勢，這正是滑坡做黏滑運動和蠕滑運動的標志[13]。目前已有的預報模型大都沒有考慮滑坡的運動速率對摩擦力的影響。與此同時，現(xiàn)有的預報模型大多數(shù)以滑坡的位移為參量建立預報模型。

基于上述分析，有必要考慮滑坡運動速率對摩擦力的影響。該文以黏滑運動和灰色系統(tǒng)為理論基礎(chǔ)，建立以速率為參量的GM(1,1)預報模型，并結(jié)合黃茨滑坡、新灘滑坡進行實例分析，取得較好的預報效果。

2 滑坡孕育特征和變形運動特性

2.1 滑坡孕育特征

圖1 不同類型的滑坡位移隨時間的變化Fig.1 Displacement of different types of landslides against time

許強[14]根據(jù)斜坡的變形-時間曲線特征將滑坡分為：穩(wěn)定型滑坡、漸進型滑坡、突發(fā)型滑坡(圖1)。該文用殘余抗剪強度[15]σr、長期抗剪強度σl、峰值抗剪強度σp、流變下限強度σcl以及滑坡的下滑應力強度σ，這5個強度指標來區(qū)分滑坡的類型。當滑坡的下滑應力強度σ小于滑坡的流變下限強度σcl時，滑坡僅發(fā)生彈性變形，不發(fā)生宏觀的位移；當σcl<σ<σl時，滑坡存在流變變形，有宏觀的位移，但經(jīng)過一段時間后，流變停止。這主要是因為在沒有其他干擾(地震、突降暴雨等)的情況下，滑坡的殘余強度σr足以抵抗滑坡的下滑應力強度σ，一般把這種類型的滑坡稱為穩(wěn)定型滑坡(圖1)；當σl<σ<σp時，滑坡存在流變變形，經(jīng)過一段時間的發(fā)展，滑坡的殘余強度σr不能抵抗滑坡的下滑力σ，變形逐漸增大，漸漸進入加速破壞階段，直至整體失穩(wěn)破壞，這種具有滑動時效的滑坡稱為漸變型滑坡(圖1)；當σp<σ時，滑坡不發(fā)生流變變形，滑動面上的各點的強度均達到或超過滑坡的破壞強度，滑坡會突然失穩(wěn)，發(fā)生突發(fā)型滑坡(圖1)。由上述分析可見，不同的受力特征，滑坡表現(xiàn)出的變形特征也不同。因此，從分析滑坡的受力狀態(tài)出發(fā)，分析滑坡的變形以及運動特性，或許可以找到預報滑坡的方法。

圖2 巖土應力-應變關(guān)系的2種基本模式Fig.2 Two basic modes of stress-strain relationship of rock and soil

大量研究表明，滑坡不是瞬間發(fā)生的，通常具有一定的變形時效。滑坡主要發(fā)在坡度為30°～50°的斜坡上[16]。突變型一般具有較深的滑動面，往往發(fā)展為液化流動，初始蠕動階段、等速蠕變階段經(jīng)歷的時間很短，發(fā)生時往往沒有預兆；漸進型滑坡規(guī)模較大，演化時間較長，具有明顯的預兆，比如地表裂縫、隆起、下陷等。已有的研究表明，大多數(shù)的自然斜坡和人工邊坡，隨著時空的演進，多發(fā)生漸進型滑坡。漸進型滑坡的滑帶土一般夾雜著具有剪脹性的超固結(jié)黏土或密實的砂土，對于剪脹性的土而言，隨著下滑力和抗滑力的變化，土的密實性和黏聚力逐漸下降，表現(xiàn)出應變軟化的現(xiàn)象(圖2)。隨著變形的發(fā)展，摩擦系數(shù)逐漸降低，巖層間的錯動和彎曲越來越明顯。從時間上看，漸進式滑坡一般會經(jīng)歷3個階段：初始變形階段、等速變形階段、加速變形階段。從空間上也可以將滑坡的變形分為3個階段[14]：第1階段，當斜坡的下滑應力強度σ超過滑坡的流變下限σcl時，巖層向坡下彎曲，后緣產(chǎn)生拉力，坡體表層產(chǎn)生蠕滑。后緣的拉力使得坡體的潛在滑移面產(chǎn)生裂縫，并不斷延伸發(fā)展，逐漸形成主裂縫；與此同時，拉力的豎直分量使得坡體下錯，造成面向坡下的臺階；第2階段，隨著蠕變的不斷發(fā)展，中部剪應力集中的部位被擾動，使得斜坡的下半部分出現(xiàn)隆起或者鼓丘，裂縫繼續(xù)發(fā)展。第3階段，隨著裂縫的發(fā)展，潛在的滑移面逐漸貫通，斜坡破壞，形成滑坡。

2.2 滑坡變形運動特性

無論是穩(wěn)定型滑坡、漸進型滑坡還是突發(fā)型滑坡，滑坡的運動方式可以歸納為2種：一種是穩(wěn)定的滑動，即滑動過程中沒有應力降，如滑坡的等速滑動階段，這種滑動形式被稱為蠕滑；另一種是不穩(wěn)定的滑動，滑動過程中伴隨有周期性的應力降。前者滑動速度較慢，后者滑動速度較快?；碌脑杏桶l(fā)生過程實際上是這兩種運動交替發(fā)生的過程，顯微鏡下更為常見的情況是黏滑和蠕滑標志往往共存一處[17]。

在滑動體發(fā)生蠕動之前，潛在滑動面附近的顆?；蛘哳w粒團在地質(zhì)作用下隨機分布，具有高度的對稱性，處于無序狀態(tài)。當斜坡的下滑應力強度σ超過滑坡的流變下限σcl時，滑坡體開始蠕動，滑帶土內(nèi)的顆粒開始發(fā)生明顯的有序運動。通過顯微鏡觀察可以發(fā)現(xiàn)，滑帶土內(nèi)的軟弱結(jié)構(gòu)部分有從內(nèi)向外擴展的微裂縫。隨著微裂縫的不斷擴展，平行于滑動面的顆?；蝾w粒團在摩擦力或者剪應力的作用下不斷被擠壓或者拖拽。少量的硬度較低的顆?；蛘哳w粒團被壓碎、拉長或者形態(tài)逐漸向扁平化方向發(fā)展[18]，而硬度較大的顆?；蛘哳w粒團則產(chǎn)生滑移；垂直于滑動面的硬度比較低的顆?；蛘哳w粒團在正壓力的作用下長軸方向開始向滑動面傾斜，而硬度較高的顆?；蛘哳w粒團開始被剪斷；隨著蠕變變形的繼續(xù)發(fā)展，平行于滑動面的大部分顆?；蛘哳w粒團在剪應力作用下沿著滑動面發(fā)生定向運動[18]。

大量硬度比較低的顆?；蛘哳w粒團被壓碎或者呈現(xiàn)扁平狀，而硬度比較高的顆粒或者顆粒團逐漸填充到孔隙率比較高的區(qū)域。垂直于滑動面的大部分硬度比較低的顆粒或者顆粒團在正壓力的作用下長軸方向不向滑動面傾斜的程度不斷加大，而硬度較高的顆?；蛘哳w粒團由于被剪斷成三角楔，在摩擦力的作用下發(fā)生轉(zhuǎn)動[18]。通過顯微鏡觀察可以發(fā)現(xiàn)，滑帶土內(nèi)的軟弱結(jié)構(gòu)部分有放射狀或者不規(guī)則的裂縫形成。隨著剪切位移的不斷增大，裂縫不斷發(fā)展，顆?；蛘哳w粒團被壓碎、啃斷的現(xiàn)象越來越嚴重，潛在的滑移面逐漸貫通斜坡破壞，形成滑坡。

通過以上的分析發(fā)現(xiàn)，滑坡的孕育、發(fā)展和發(fā)生過程實際上是蠕滑和黏滑運動交替發(fā)生的過程。

3 滑坡演化的預測預報模型

3.1 傳統(tǒng)灰色預測預報模型及以速率為參量的GM(1,1)模型

3.1.1 以位移為參量的GM(1,1)預報模型

1986年，我國學者陳明東等[3]首次以位移切線角為參量建立灰色系統(tǒng)模型，對新灘滑坡進行了預報驗證，并取得成功。由于位移切線角獲取較困難，此小節(jié)將以位移為參量進行建模預報。

灰色系統(tǒng)模型的灰微分方程為

(1)

式中a1，b1為系數(shù)。對于一組等時距位移觀測序列

x(0)=[x(0)(1),x(0)(2),…,x(0)(n-1),x(0)(n)],

x(0)一次累加后得到

用最小二乘法求解系數(shù)a1，b1，計算表達式為

(2)

其中：

(3)

Z(1)(tk)=0.5(x(1)(tk-1)+x(1)(tk)),

k=2,3,…,n;

(4)

Y=[x(0)(t2),x(0)(t3),…,x(0)(tn)]T。

(5)

求解灰微分方程得到位移的時間響應公式為

(6)

式(6)兩邊對時間求導，得到速度隨時間的變化率為

3.1.2 以位移為參量的Verhulst預報模型

Verhulst是1987年德國物理學家發(fā)現(xiàn)的一種生物生長模型，由于滑坡的位移特征與生物的生長規(guī)律類似，很多學者用Verhulst模型對滑坡進行預報研究[4-5]。

Verhulst模型的灰微分方程為

(8)

式中a2，b2為系數(shù)。

Verhulst模型系數(shù)的求解方法和過程同GM(1,1)模型，唯一不同的是需要將式(3)中的矩陣B替換成如下形式，即

模型位移隨時間的演化方程為

3.1.3 以速率為參量的的GM(1,1)預報模型

絕大多數(shù)機械運動中存在摩擦振動現(xiàn)象(黏滑運動)，滑坡運動也不例外。為了深入認識這種運動，學者們對這種現(xiàn)象進行大量的研究。我國學者李春波[19]進行了一系列試驗，并提出用下述方程來描述摩擦力與滑動速度之間的依賴關(guān)系，即

(11)

(12)

取式 (12)前2項代入式(1)得

(13)

F=uNcosα-kv。

(14)

式中：u為滑動摩擦系數(shù)；N為垂直于坡面的支持力。

根據(jù)摩爾-庫倫準則，可以將抗滑力計算式表達為

F抗滑=uNcosα-kv+c。

(15)

式中c為黏聚力。

滑坡受自重力及內(nèi)力作用，且內(nèi)力平行于坡面(圖3)。

圖3 滑坡的簡化滑塊模型Fig.3 Simplified sliding model of landslide

則由牛頓第二定律可知：

(17)

式(17)即為描述邊坡演化的非線性動力學模型。

該模型是以速率為參量的灰色GM(1,1)預測模型。模型的求解方法同以位移為參量的灰色GM(1,1)預測模型相比，只需將位移參量替換成速率參量。

以速率為參量的GM(1,1)預測模型的速率計算式為

式 (18)對時間積分，可以得到以速率為參量的灰色GM(1,1)預測模型的位移時間公式為

(19)

3.2 模型預測預報判據(jù)和時間

許強等[20]認為改進的切線角>85°，斜坡變形進入加速變形階段。此方法的切線角85°判據(jù)是在對一個滑坡位移曲線進行分析的基礎(chǔ)上得出的，缺乏一般性。王家鼎[21]通過對大量滑坡實例的統(tǒng)計分析，提出將位移切線角89°～89.5°作為滑坡的預警判據(jù)。此后，李聰?shù)萚22]建立了滑坡數(shù)據(jù)庫，對不同類型的滑坡的臨界速率進行了統(tǒng)計，發(fā)現(xiàn)滑坡進入急劇變形階段位移臨界速率一般≤50 mm/d。根據(jù)切線角與速率之間的轉(zhuǎn)換關(guān)系v=tanθ，計算得到滑坡劇烈滑動時的位移切線角θ一般≤88.855°。由統(tǒng)計規(guī)律可知，滑坡實例數(shù)據(jù)樣本量越大，越能反映滑坡的變形規(guī)律。故到目前為止，根據(jù)切線角變量來確定預警判據(jù)的研究中，李聰?shù)呐R界速率判據(jù)最為科學和嚴密；賀小黑等[5]認為Verhulst模型以速率最大值為預報判據(jù)缺乏合理的依據(jù)，以加速度最大值作為預報判據(jù)更為合理。綜合上述分析，對于GM(1,1)模型，本文取臨界速率vcr=50 mm/d為滑坡發(fā)生的判據(jù)；對于Verhulst模型，本文以加速度最大值為判據(jù)。將vcr分別代入式(7)、式(18)，反解出滑坡發(fā)生時間。

以位移為參量的GM(1,1)模型預報時間的公式為

以速度為參量的GM(1,1)模型預報時間的公式為

以位移為參量的Verhulst模型預報時間的公式[6]為

(22)

式(20)、式(21)、式(22)中的Δt為位移監(jiān)測的等時間間隔天數(shù)。

4 實例說明

為了驗證以速率為參量的GM(1,1)模型的預測預報效果，分別對黃茨滑坡、新灘滑坡進行理論計算和分析。值得注意的是，實際監(jiān)測過程中，滑坡位移監(jiān)測數(shù)據(jù)往往會出現(xiàn)波動，速率會出現(xiàn)負值。為了減弱觀測誤差對數(shù)據(jù)分析的影響，可先用Vondrak平滑法[23]對原始數(shù)據(jù)進行預處理，處理后，如果速率沒有負值，則可直接進行建模計算。若速率依然有負值，可選用滑坡進入加速變形階段后的數(shù)據(jù)進行計算和分析。 為了方便進行對比分析，將以位移為參量的GM(1,1)模型編號為模型①，以位移為參量的Verhulst模型編號為模型②，以速率為參量的GM(1，1)模型編號為模型③。

4.1 實例1——黃茨滑坡

黃茨滑坡位于甘肅省永靖縣,該滑坡起始監(jiān)測時刻為1994年8月1日，監(jiān)測時間間隔為15 d，劇滑時間為1995年1月30日凌晨2點30分，監(jiān)測資料[5]如表1。

表1 黃茨滑坡B2監(jiān)測點的監(jiān)測數(shù)據(jù)Table 1 Monitored displacements at monitoring point B2

建模前，先對原始數(shù)據(jù)進行平滑處理，平滑方法采用Vondrak平滑法[23]，該法對等時距觀測數(shù)據(jù)和不等時間距離都能進行平滑處理。平滑效果如圖4。

圖4 黃茨滑坡原位移和平滑位移與測次曲線Fig.4 Original and fitted cumulative displacement vs.time curves of Huangci landslide

對于以位移為參量的GM(1,1)和Verhulst模型，用于建模的原始序列為

x(0)=[1.05，3.26，6.07，8.2，10.7，14.99，18.05，21.65，27.4，33.14，38.56，49.43，63.54]。

設(shè)時間序列為t=[1,2,3,4,5,6,7，8，9，10，11，12，13]，對x(0)進行一次累加后有：x(1)=[1，4.31，10.38，18.58，29.28，44.27，62.32，83.97，111.37，144.51，183.07，232.5，296.04]。

對于以速率為參量的GM(1,1)預報模型，由于滑坡監(jiān)測的物理量是位移，故需要將位移轉(zhuǎn)化為速率，才能用速率模型參量進行建模預報。位移轉(zhuǎn)化為速率的計算公式如下,即

計算得到建模的速率序列為

v(0)=[2.21，2.81，2.13，2.5，4.29，3.06，3.6，5.75，5.42，10.87，14.11]。

以1994年8月16日為速率模型建模的初始時間，v(0)對應的時間序列為

vt=[1,2,3,4,5,6,7,8,9,10,11,12]。

對v(0)進行一次累加后有

v(1)=[2.3，5.1，7.3，9.5，14.3，17.1，20.6，26.5，32.6,38,47.7,63.3],用最小二乘法求解，得到模型參數(shù)a1=-0.217 61，b1=5.155；a2=-0.223 98，b2=0.424 55；a3=0.308 06，b3=0.000 29。將求得的系數(shù)以及v(1)(t1)，vcr分別代入式(20)、式(21)、式(22)得模型的預報時間值tr1=108.779 26，tr2=115.715 41，tr3=271.671 06。需要說明的是，位移監(jiān)測單位是cm,vcr取5 cm/d；由于以速率為參量的GM(1,1)模型的計算時間從1994年8月16日開始，故代入式(20)時需要將式(20)的t替換成vt計算。把預報時間值轉(zhuǎn)換成日期，得到滑坡預報結(jié)果見表2。表中，Δ是預報絕對誤差(d)，是預報時間與滑坡實際發(fā)生時間的差值，即Δ=tcr-ts；δ為相對誤差(%)，相對誤差是絕對誤差Δ換算成天數(shù)后除以滑坡實際發(fā)生時間ts與時間序列t1之間的間隔天數(shù)再乘以100%。當滑坡預報時間與實際發(fā)生時間相同時，Δ=0，δ=0；當滑坡預報時間早于實際發(fā)生時間時，Δ>0，δ>0；當滑坡預報時間晚于實際發(fā)生時間時，Δ<0，δ<0。下文中的計算規(guī)則與此處相同。將預報結(jié)果列于表2中。

表2 黃茨滑坡模型預報結(jié)果Table 2 Forecasted occurrence time of HuangciLandslide by different models

圖5 黃茨滑坡模型位移預測值-時間曲線Fig.5 Comparison of predicted displacement vs.time curve among three models

從表2中可以看出，3個預測預報模型的預報誤差均在50%以內(nèi)，GM(1,1)模型的預報時間比Verhulst模型的預報時間早，能夠提前預報滑坡的發(fā)生。以速率為參量的GM(1,1)模型的預報時間更加接近滑坡實際發(fā)生時間，精度較高。

將求得的模型參數(shù)分別代入式(6)、式(10)、式(19)可以求得模型位移預測值。圖5是3個預測模型的位移時間曲線與原位移時間曲線和平滑位移曲線的對比圖。從圖5中可以看出，以速率為參量的GM(1,1)預測模型的位移預測效果明顯優(yōu)于Verhulst模型，基本上能夠反映這段時間內(nèi)位移隨時間的變化趨勢。

4.2 實例2——新灘滑坡

新灘滑坡位于湖北省秭歸縣境內(nèi)的新灘鎮(zhèn)長江北岸,其下游27 km 處是三峽大壩. 滑坡發(fā)生于1985年6月12日。將新灘滑坡B3監(jiān)測點的位移數(shù)據(jù)[7]列于表3中。為了分析以速率為參量的GM(1,1)模型的預測預報效果,首先對監(jiān)測點的位移數(shù)據(jù)進行平滑處理。平滑效果見圖6。

圖6 新灘滑坡原位移和平滑位移隨時間的變化曲線Fig.6 Original and fitted cumulative displacement vs.time curves of Xintan landslide

表3 新灘滑坡B3監(jiān)測點位移Table 3 Monitored displacements of point B3 of Xintan landslide

用平滑后的位移全部數(shù)據(jù)建立GM(1,1)模型。從圖中可以看出，第36個月以后，位移速率不斷增大，故認為從第37個月開始，滑坡進入加速變形階段。

對于Verhulst模型，選取第36個月以后的平滑數(shù)據(jù)進行預測預報。分別用最小二乘法求解得模型系數(shù)a1=-0.048 46,b1=110.605 22，a2=-0.046 45,b2=4.413 26，a3=0.066 64,b3=1.204 5×10-7。將求得的系數(shù)以及v(1)(t1)，vcr分別代入式(20)、式(21)、式(22)得模型的預報時間值tr1=47.423 99，tr2=54.002 07，tr3=79.789 39。將預報結(jié)果列于表4中。

表4 新灘滑坡模型預報結(jié)果Table 4 Forecasted occurrence time of XintanLandslide by different models

注:amax表示滑坡運動的加速度達到最大值

從表4中可知，3個預測預報模型的誤差均在50%以內(nèi)，GM(1,1)模型的預報時間比Verhulst模型的預報時間早，能夠提前預報滑坡的發(fā)生。以速率為參量的GM(1,1)模型的預報時間能做到提前預警。

將求得的模型參數(shù)分別代入式(6)、式(10)、式(19)可以求得模型位移預測值。圖7是3個預測模型的位移時間曲線與原位移時間曲線和平滑位移曲線的對比圖。

圖7 新灘滑坡模型位移預測值-時間曲線對比Fig.7 Comparison of predicted displacement vs. timecurve among three models

從圖7中可以看出，以速率為參量的GM(1,1)預測模型的位移預測效果明顯優(yōu)于Verhulst模型，基本上能夠反映位移隨時間的變化趨勢。

5 結(jié) 論

該文以摩擦學中的黏滑運動以及灰色系統(tǒng)為理論基礎(chǔ)，建立了以速率為參量的GM(1,1)預報模型，并結(jié)合黃茨滑坡、新灘頭滑坡實例進行了預測預報分析，主要得到以下結(jié)論：

(1)基于監(jiān)測位移建立模型對進行滑坡預測預報時，所建立的模型應該考慮滑坡的變形速率對摩擦力的影響。

(2)該文建立的以速率為參量的GM(1，1)預報模型的預報效果比傳統(tǒng)的以位移為參量的GM(1，1)預報模型和Verhulst預報模型更好。

(3)以速率為參量的GM(1,1)滑坡預測預報模型比以位移為參量的Verhulst預測預報模型更能夠反映位移隨時間的變化規(guī)律。