基于MatDEM的煙家溝滑坡演化過程數值模擬分析

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

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

0 引言

2015年8月12日，陜西省山陽縣煙家溝村發生一起突發性滑坡災害，坡體的總體積達168×104m3，經過煙家溝支溝左、右兩側山體的先后阻擋，滑體的運動方向先后兩次發生偏轉，最終沿溝谷下游運動，給當地居民的生命財產造成了嚴重損失［1］.

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

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

1 滑坡區域地質條件

1.1 地形地貌

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

1.2 地層巖性

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

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

圖2 滑前斜坡地質結構剖面圖Fig.2 Geological profile of the slope before sliding

1.3 地質構造

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

1.4 滑坡運動過程特征

據野外調查發現，滑源區巖土體在重力作用下失穩破壞以66°方向開始滑動，在受到煙家溝左側山體碰撞后朝105°方向繼續滑動，進而受到煙家溝右側山體阻擋后轉向60°方向運動直至停止（圖3）.滑坡運動過程呈現典型的右旋特征.

表1 滑坡區地層巖性Table 1 Stratigraphic lithology of landslide area

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

2 數值模型建立

2.1 顆粒流離散元程序

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

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

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

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

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

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

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

其中，Fs0為初始抗剪力，μp為摩擦系數.在實際模擬過程中，切向力大于Fsmax時，單元膠結亦發生斷裂，單元產生錯動，此時單元間只剩滑動摩擦－μpFn.

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

2.2 滑坡模型建立

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

3 模擬參數選取

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

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

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

其中，d為顆粒直徑，I為比例系數，

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

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

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

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

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

4 模擬結果及分析

4.1 滑坡滑移距離變化特征

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

4.2 滑坡速度變化特征

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

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

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

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

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

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

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

4.3 滑坡堆積厚度分布特征

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

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

5 討論

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

6 結論

通過上述研究取得了以下結論：

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

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

3）經過MatDEM數值模擬后所呈現出的堆積特征和現場調查結果完全符合，證明矩陣離散元法模擬滑坡演化過程結果具有可靠性.

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