基于MatDEM的煙家溝滑坡演化過(guò)程數(shù)值模擬分析

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

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

0 引言

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

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

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

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

1.1 地形地貌

滑坡現(xiàn)場(chǎng)遙感照片及斜坡剖面如圖1、2.滑坡區(qū)位于商洛市山陽(yáng)縣中村鎮(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ū)巖性復(fù)雜，地層出露較多，主要地層為震旦系燈影組、寒武系下統(tǒng)水溝口組、寒武系中統(tǒng)岳家坪組和在溝谷地帶大量分布的第四系［11］.詳見(jiàn)表1.

圖1 山陽(yáng)煙家溝滑坡遙感圖像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ā)育，且整個(gè)區(qū)域地層層序倒轉(zhuǎn)，使得水溝口組地層處于震旦系燈影組的地層之下，呈平行不整合接觸.坡體結(jié)構(gòu)上部堅(jiān)硬，下部軟弱.在區(qū)域構(gòu)造背景下，斷層和褶皺呈東西向分布，地應(yīng)力較為集中，巖石風(fēng)化程度高，斜坡破壞現(xiàn)象經(jīng)常發(fā)生.

1.4 滑坡運(yùn)動(dòng)過(guò)程特征

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

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

圖3 滑坡運(yùn)動(dòng)過(guò)程分析Fig.3 Movement process of landslide

2 數(shù)值模型建立

2.1 顆粒流離散元程序

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

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

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

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

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

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

Xd為斷裂變形.當(dāng)單元在法向上外力逐漸增加至超過(guò)Fnmax時(shí)，單元間膠結(jié)斷開(kāi)，拉力消失.單元膠結(jié)所能承受最大切向力滿(mǎn)足：

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

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

2.2 滑坡模型建立

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

3 模擬參數(shù)選取

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

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

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

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

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

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

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

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

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

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

4.1 滑坡滑移距離變化特征

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

4.2 滑坡速度變化特征

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

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

圖7 滑坡運(yùn)動(dòng)過(guò)程模擬Fig.7 Simulation of landsliding process

圖8 滑坡模型單元速度－時(shí)間變化曲線(xiàn)Fig.8 Speed vs.time curves of landslide model unit

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

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

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

4.3 滑坡堆積厚度分布特征

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

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

5 討論

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

6 結(jié)論

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

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

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

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

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