栗曉松，范 文，2，曹琰波，全倬梁

1.長安大學 地質(zhì)工程與測繪學院，陜西 西安710054；2.工業(yè)和信息化部電子綜合勘察設計院，陜西 西安710054

0 引言

2015年8月12日，陜西省山陽縣煙家溝村發(fā)生一起突發(fā)性滑坡災害，坡體的總體積達168×104m3，經(jīng)過煙家溝支溝左、右兩側(cè)山體的先后阻擋，滑體的運動方向先后兩次發(fā)生偏轉(zhuǎn)，最終沿溝谷下游運動，給當?shù)鼐用竦纳敭a(chǎn)造成了嚴重損失［1］.

正確認識滑坡的運動、堆積特征是對滑坡災害進行合理的評價，提出合理的治理措施的前提.滑坡體滑動速度、滑動距離和堆積特征是對滑坡災害進行有效評估的重要參考指標.滑坡的發(fā)生一般較突然且危害程度大［2-3］，許多學者針對滑坡運動速度、堆積特征等方面開展了廣泛的研究.魯曉兵等從理論分析出發(fā)，主要研究滑坡坡角、土體的內(nèi)摩擦角對斜坡失穩(wěn)下滑的運動特征和最終堆積分布的影響［4］；樊曉一等從滑槽模型試驗入手，對滑坡體碎屑流前緣運動速度研究，重點考慮坡度、顆粒級配及坡角等因素［5］；郝明輝等在現(xiàn)場試驗的基礎上，結(jié)合調(diào)查結(jié)果，具體地研究了滑坡體在下滑過程中顆粒分選性、反序列堆積的問題［6］；Hungr等人采用數(shù)值模擬的方法，基于DAN模擬了滑坡發(fā)生過程，得出溝谷形態(tài)是影響滑坡運動特征和堆積分布的重要影響因素［7］.目前，雖大多數(shù)學者采用多種研究方法，如現(xiàn)場調(diào)查研究、理論分析計算、模型試驗和數(shù)值模擬等對滑坡運動過程進行研究，但是總體來說針對滑坡運動發(fā)生過程的直觀性和整體性呈現(xiàn)還相對較弱.

基于此，本文以山陽縣煙家溝滑坡作為研究對象，利用高性能離散元軟件MatDEM［8-9］，依據(jù)無人機現(xiàn)場實拍，獲得煙家溝滑坡高程信息，經(jīng)ArcGIS初步柵格化處理后，結(jié)合滑坡實際特征建立初始三維地質(zhì)模型，并對模型賦予真實的力學性質(zhì)參數(shù)，模擬真實情況下滑坡啟動發(fā)生的全過程.

1 滑坡區(qū)域地質(zhì)條件

1.1 地形地貌

滑坡現(xiàn)場遙感照片及斜坡剖面如圖1、2.滑坡區(qū)位于商洛市山陽縣中村鎮(zhèn)，為大西溝和煙家溝的交界處.滑坡區(qū)總體上屬基巖中山陡坡區(qū)，地質(zhì)構(gòu)造發(fā)育，處于倒轉(zhuǎn)褶皺的一側(cè)［10］.滑坡區(qū)原始斜坡地形陡峻，坡度約為35°，高程1010～1300 m，高差290 m.前緣大西溝溝谷深切，坡向與地層產(chǎn)狀一致，呈突出的山梁，構(gòu)成陡傾順層邊坡.

1.2 地層巖性

滑坡區(qū)巖性復雜，地層出露較多，主要地層為震旦系燈影組、寒武系下統(tǒng)水溝口組、寒武系中統(tǒng)岳家坪組和在溝谷地帶大量分布的第四系［11］.詳見表1.

圖1 山陽煙家溝滑坡遙感圖像Fig.1 Remote sensing image of Yanjiagou landslide in Shanyang County

圖2 滑前斜坡地質(zhì)結(jié)構(gòu)剖面圖Fig.2 Geological profile of the slope before sliding

1.3 地質(zhì)構(gòu)造

滑坡區(qū)地處秦嶺南部，屬薄皮逆沖推覆構(gòu)造帶，是舒家壩－劉嶺逆沖推覆構(gòu)造帶和鳳縣－鎮(zhèn)安逆沖推覆構(gòu)造帶的接合部位［12］，斷層和褶皺較發(fā)育，且整個區(qū)域地層層序倒轉(zhuǎn)，使得水溝口組地層處于震旦系燈影組的地層之下，呈平行不整合接觸.坡體結(jié)構(gòu)上部堅硬，下部軟弱.在區(qū)域構(gòu)造背景下，斷層和褶皺呈東西向分布，地應力較為集中，巖石風化程度高，斜坡破壞現(xiàn)象經(jīng)常發(fā)生.

1.4 滑坡運動過程特征

據(jù)野外調(diào)查發(fā)現(xiàn)，滑源區(qū)巖土體在重力作用下失穩(wěn)破壞以66°方向開始滑動，在受到煙家溝左側(cè)山體碰撞后朝105°方向繼續(xù)滑動，進而受到煙家溝右側(cè)山體阻擋后轉(zhuǎn)向60°方向運動直至停止（圖3）.滑坡運動過程呈現(xiàn)典型的右旋特征.

表1 滑坡區(qū)地層巖性Table 1 Stratigraphic lithology of landslide area

圖3 滑坡運動過程分析Fig.3 Movement process of landslide

2 數(shù)值模型建立

2.1 顆粒流離散元程序

離散元法（DEM）是由Cundall等［13］首次提出，主要是用以研究顆粒狀物質(zhì)的運動及相互作用.隨著近年來計算機技術(shù)進步和離散元理論的發(fā)展，離散元法在多種領域得到廣泛應用［14］，包括固體力學領域中的顆粒堆積體力學性質(zhì)研究［15］，地質(zhì)領域中各類構(gòu)造的形成和演化研究等.近年來很多研究人員［16-17］采用離散元數(shù)值模擬方法對滑坡運動過程進行反演分析，取得了較好的成果.因此，本文選擇矩陣離散元法進行煙家溝滑坡運動數(shù)值模擬.

MatDEM采用矩陣離散元計算法和三維接觸算法，在MATLAB語言運行環(huán)境下自主建立模型并完成離散元生態(tài)下的能量守恒計算.該算法運行速度更快，模擬效果更直觀，對模擬滑坡、巖爆、撞擊破壞、樁土作用等具有突出的優(yōu)點.MatDEM軟件可將一定數(shù)量的具有特定力學性質(zhì)的單元按照相應組合方式膠結(jié)在一起，建構(gòu)出多種符合實際結(jié)構(gòu)性的三維巖土體（圖4）.本文采用線彈性接觸方式建立滑坡模型，單元間膠結(jié)力可用法向彈簧力和切向彈簧力來表示（圖5）.彈簧接觸和切向彈簧接觸.

圖4 顆粒膠結(jié)三維球體模型Fig.4 3D sphere model of particle cementing

通常情況下，顆粒單元之間的法向力和法向變形、剪切力和剪切變形都維持在穩(wěn)定的限度內(nèi).滿足：

圖5 線彈性接觸方式示意圖（據(jù)文獻［18］）Fig.5 Schematic diagram of linear elastic contact mode（From Reference［18］）

其中，F(xiàn)n為法向力，Kn為法向剛度，Xn為法向變形；Fs為剪切力，Ks為切向剛度，Xs為剪切變形.顆粒單元間的最大法向力滿足：

Xd為斷裂變形.當單元在法向上外力逐漸增加至超過Fnmax時，單元間膠結(jié)斷開，拉力消失.單元膠結(jié)所能承受最大切向力滿足：

其中，F(xiàn)s0為初始抗剪力，μp為摩擦系數(shù).在實際模擬過程中，切向力大于Fsmax時，單元膠結(jié)亦發(fā)生斷裂，單元產(chǎn)生錯動，此時單元間只剩滑動摩擦－μpFn.

MatDEM數(shù)值模擬計算符合力－位移定律和牛頓第二定律，模型時間步（dT）很小，此間顆粒的加速度及速度保持不變，通過時間步迭代來計算顆粒的受力、加速度、速度和位移，通過反復迭代最終實現(xiàn)離散元法動態(tài)模擬［18］.

2.2 滑坡模型建立

建立準確的三維模型是模擬真實滑坡的基礎，通過無人機航拍獲取高精度高程坐標，經(jīng)ArcGIS軟件柵格化后，導入MatDEM軟件進行后續(xù)處理分析，最終建立起山陽煙家溝滑坡的三維地形模型（圖6）.通過滑前、滑后數(shù)字高程模型對比并結(jié)合現(xiàn)場實際勘察發(fā)現(xiàn)，滑坡體堆積分布分區(qū)較明顯，無明顯刮鏟區(qū)，且?guī)r體單元的半徑小，坡體前緣破碎嚴重，故依據(jù)坡體形態(tài)和物質(zhì)組成特點把滑坡區(qū)域概括為滑源區(qū)和堆積區(qū).

3 模擬參數(shù)選取

顆粒單元的力學性質(zhì)通常會在其重新堆積和膠結(jié)后發(fā)生較大改變.在本研究中，單元的力學參數(shù)主要通過宏微觀參數(shù)轉(zhuǎn)換公式計算以及數(shù)值模擬測試來確定.通過采用這樣的參數(shù)建立的模型，能夠使得巖土體在最大程度上符合天然狀態(tài)下的結(jié)構(gòu)性和復雜性.

圖6 滑坡前、后數(shù)字高程模型及分區(qū)Fig.6 Digital elevation model and zoning before and after landsliding

宏微觀參數(shù)轉(zhuǎn)換公式計算：線彈性接觸模型的5個微觀力學參數(shù)，分別是摩擦系數(shù)（μp）、初始抗剪力（Fs0）、正向勁度系數(shù)（Kn）、切向勁度系數(shù)（Ks）和斷裂位移（Xb），這些力學參數(shù)經(jīng)由所選材料的5個宏觀力學性質(zhì)（泊松比v、楊氏模量E、抗壓強度Cu、抗拉強度Tu和內(nèi)摩擦系數(shù)μi）結(jié)合宏微觀參數(shù)轉(zhuǎn)換公式計算得到.轉(zhuǎn)換公式如下所示［19］：

其中，d為顆粒直徑，I為比例系數(shù)，

獲得參數(shù)具體步驟：①經(jīng)過室內(nèi)物理實驗測得巖體力學參數(shù)值（表2），結(jié)合轉(zhuǎn)換公式求取初始數(shù)值模型的物理力學參數(shù)，將參數(shù)賦值給顆粒單元并建立滑坡堆積模型；②運行MatDEM內(nèi)置的抗拉強度、單軸壓縮和抗壓強度測試程序，求取模型測試的楊氏模量及抗拉強度等性質(zhì)；③將數(shù)值測試的模型力學性質(zhì)與實際測量的巖體力學性質(zhì)對比得到一個尺度值；④最后把實際測量的巖體力學參數(shù)與尺度值對比，再次放到上述轉(zhuǎn)換公式中計算，把此次的單元力學參數(shù)賦值給顆粒單元.

表2 實測滑坡體巖石力學性質(zhì)Table 2 Mechanical properties of measured landslide dolomite

經(jīng)過至少3次②③④過程可獲取誤差小于2%的單元力學參數(shù).多次迭代計算后得出用于滑坡模型單元的宏微觀力學參數(shù)值（如表3）.

表3 宏、微觀力學參數(shù)值Table 3 Macro and micro mechanical parameters

根據(jù)前人的試驗研究及其相關(guān)數(shù)值模擬結(jié)果［20］得知，巖土體的抗壓強度、抗拉強度及彈性模量值有“尺度效應”的現(xiàn)象.在試驗室當中實測的小尺寸巖體的物理力學性質(zhì)通常情況下要比現(xiàn)場的大尺度巖體的物理力學性質(zhì)偏大.本文數(shù)值模擬計算過程實際使用的強度參數(shù)值約等于室內(nèi)實測試驗值的40%.

4 模擬結(jié)果及分析

4.1 滑坡滑移距離變化特征

結(jié)合表3中所得模擬參數(shù)進行迭代計算，對煙家溝滑坡運動情況進行反演分析，模擬過程如圖7所示.圖7a、b、c、d分別截取了第6、10、16、28 s時的位移距離.經(jīng)讀取模型數(shù)據(jù)得，滑坡堆積體長度約為510 m，最大滑移距離約為570 m，與現(xiàn)場調(diào)查得到的最遠滑動距離600 m較接近［11］，進一步驗證了選取參數(shù)的合理性，有效地再現(xiàn)了煙家溝滑坡滑動和堆積過程.

4.2 滑坡速度變化特征

通過提取滑坡體活動單元的速度、滑坡表面跳躍單元的速度及對應時間，生成滑坡巖土體速度－時間變化曲線（圖8）.滑坡體內(nèi)部在應力條件下發(fā)生失穩(wěn)破壞，在5～12 s時間內(nèi)完成了啟動、加速過程，最大速度達到44 m/s；滑坡啟動約9 s時，滑體平均速度達到最大值，約23.5 m/s，滑坡啟動28 s后滑體平均速度幾乎收斂于0，滑體運動趨于停止，僅有零星的表面單元向下跳躍運動.

根據(jù)滑坡運動速度、滑移距離及堆積位置，將滑坡堆積過程分為加速碰撞堆積階段、整體滑動堆積階段和減速堆積階段.

圖7 滑坡運動過程模擬Fig.7 Simulation of landsliding process

圖8 滑坡模型單元速度－時間變化曲線Fig.8 Speed vs.time curves of landslide model unit

加速碰撞堆積階段（圖7b）：滑源區(qū)前緣顆粒單元碰撞解體滑落，在內(nèi)部應力釋放及重力作用下沿著坡向運動，到達坡腳下“V”型谷時平均速度達到約23 m/s.由于受到溝谷左側(cè)山體的阻擋，大部分滑體單元在此處堆積；其余顆粒單元在慣性力影響下以較高速度向前沖出，并順著溝谷方向繼續(xù)運動.

整體滑動堆積階段（圖7c）：剩余滑體單元整體以較高速度下滑并到達第二次碰撞點，在此階段溝谷寬度變大，部分滑坡單元能量得以充分釋放并沿途堆積，其余部分滑體單元會沿溝谷發(fā)生偏轉(zhuǎn)繼續(xù)運動.

減速堆積階段（圖7d）：在16～28 s時滑體單元速度逐漸降低，直至到達坡底，主體僅發(fā)生局部滑動和調(diào)整，堆積基本完成.

4.3 滑坡堆積厚度分布特征

根據(jù)滑坡前后高程數(shù)據(jù)對比，獲取滑坡堆積厚度分布（圖9）.煙家溝滑坡發(fā)生后，滑源區(qū)單元碰撞解體滑落，故滑源區(qū)堆積厚度為負值，從圖中可以看出滑源區(qū)最大深度約為43 m.隨著滑體下滑，在運動方向會受到溝谷左側(cè)山體的阻擋，顆粒之間相互擠壓碰撞，導致大部分的顆粒單元能量損耗［21］.顆粒大量在溝谷處堆積，使得平均堆積厚度大幅增加，最大堆積厚度為48 m.另外部分單元在慣性作用力下保持原有運動趨勢向下游方向繼續(xù)運動，此階段溝谷寬度變大，顆粒堆積厚度減小，且有小范圍“爬高現(xiàn)象”.剩余顆粒單元會沿下游方向繼續(xù)運動并不斷沿途堆積，堆積厚度逐漸減小.堆積區(qū)堆積厚度呈現(xiàn)出先增加后減少的總體變化趨勢.

圖9 沿滑動方向的堆積厚度分布Fig.9 Distribution of accumulation thickness along the sliding direction

5 討論

本文基于MatDEM將巨量的顆粒單元用于數(shù)值建模，借助無人機航拍等手段獲得準確高程信息，于三維模型中呈現(xiàn)出復雜的滑坡要素特征，改進了傳統(tǒng)二維模型在滑坡運動發(fā)生過程中直觀性和整體性呈現(xiàn)不足的缺陷，很好地再現(xiàn)了煙家溝滑坡的運動、堆積過程.盡管精確的滑坡高程信息和MatDEM軟件高效的矩陣離散元計算法在滑坡演化過程模擬中具有無可比擬的優(yōu)勢，但在地質(zhì)構(gòu)造、坡度、顆粒級配等因素如何具體影響滑坡的啟動發(fā)生問題上還有待進一步研究.實際的滑坡運動及堆積過程非常復雜，影響因素眾多，后續(xù)還需要進一步結(jié)合理論推導、模型實驗等手段來進一步分析.本文旨在通過以上的數(shù)值模擬研究，為以后類似滑坡的運動和堆積特征分析、成災范圍劃定及災害防治提供一定的技術(shù)參考.

6 結(jié)論

通過上述研究取得了以下結(jié)論：

1）山陽縣煙家溝滑坡從啟動到最終靜止整個過程持續(xù)時間約28 s，滑坡體單元最大速度可達44 m/s，平均速度最高可達23.5 m/s，故將煙家溝滑坡歸為高速滑坡；滑坡堆積體長度約為510 m，最大滑移距離約為570 m，最大堆積厚度約為48 m.該類滑坡對人員及建筑物具有較強的破壞性.

2）根據(jù)滑坡運動速度和距離將滑坡堆積過程分為加速碰撞堆積階段、整體滑動堆積階段和減速堆積階段.順著滑坡運動路徑的方向，堆積厚度呈現(xiàn)先增加后減小的總體變化趨勢，堆積厚度在堆積體的后部最大.

3）經(jīng)過MatDEM數(shù)值模擬后所呈現(xiàn)出的堆積特征和現(xiàn)場調(diào)查結(jié)果完全符合，證明矩陣離散元法模擬滑坡演化過程結(jié)果具有可靠性.

致謝：本文使用的高性能離散元軟件MatDEM由南京大學劉春團隊自主研發(fā)，模擬過程中使用的代碼為在其源代碼基礎上進行的改編.