基于FLAC3D的貴州馬達嶺煤炭采動型滑坡機理研究

潘網生，傅良同，李 鑫，張曉周，韓國棟

(1.黔南民族師范學院 旅游與資源環境學院，貴州 都勻 558000；2.長安大學 環境科學與工程學院，陜西 西安 710054)

貴州省喀斯特出露面積約11×104km2，約占全省國土面積61.9%[1-3]，其喀斯特地貌形態以錐狀為主，類型多樣、形成發育過程相對復雜，其形態結構獨特、錐形明顯，單體高度為60～250m，坡度為45°～ 47°[4-6]。近幾年，貴州喀斯特區域滑坡地質災害頻發，生命與財產損失巨大，如關嶺縣山體滑坡(2010年6月28日)，死亡99人；福泉市道坪鎮英坪村大型山體滑坡(2014年8月27日)，死亡23人；大方縣理化鄉偏坡村山體滑坡(2016年7月1日)，死亡21人；玉屏縣田坪鎮巖屋口村山體滑坡(2017年7月7日)，死亡2人；納雍縣張家灣鎮大型山體滑坡(2017年8月28日)，死亡15人；三荔高速荔波段滑坡(2017年9月23日)，死亡3人，六盤水六枝特區關寨鎮大型山體滑坡(2018年7月9日)，因預警及時以致無人員傷亡。上述滑坡地質災害已經成為制約貴州社會經濟可持續發展的重要因素之一。由于喀斯特滑坡災害機理研究可以揭示滑坡災害成因、發展及演變規律，也可以為早期預警和地質災害防治提供科學依據和理論支撐。因此，研究煤炭采動型滑坡機理具有重要的實踐意義。

國內學術界關于貴州喀斯特滑坡機理研究，總體概括為巖性和巖土結構等(內因)，降雨和不科學的人類工程活動等(外因)。分以下兩種類型：其一是由自然地質沉積所形成的軟弱結構層在地表滲流和地下水溶蝕作用下，其粘聚力和內摩擦角減小，成為天然滑動面，導致軟弱結構層上方巖土體靜力平衡被打破，進而失穩致滑。如貴州關嶺“6.28” 喀斯特山體滑坡地質災害等[7-13]；其二是人為斜坡開挖、礦產資源開采等人類工程活動導致邊坡應力平衡被打破，在降雨或震動條件下誘發山體滑坡。如貴州印江巖口大型喀斯特山體滑坡地質災害等[14-18]。

貴州都勻馬達嶺滑坡屬于第二種類型—人為開挖煤層導致采空區上方巖體滑坡，也稱煤炭采動型滑坡[19]。趙建軍[20]基于物理模型研究其地質力學模式；李禹霏[21]通過自動化監測研究其穩定性；亓星[22]研究其滑坡碎屑流形成泥石流的特征；王玉川[23]對其開展巖石力學試驗。本文在地質野外調查基礎上，擬采用有限差分方法研究采動滑坡的變形過程及特征，以揭示其滑坡機理和驅動模式，相關研究成果可以應用于類似喀斯特礦山開采的安全管理和地質災害防治。

1 數據來源及研究方法

研究區域地質資料來源于前人研究成果及貴州省區域地質志(1987版、2017版)[24，25]，地形地貌數據來源于國家地理空間數據云服務平臺及野外現場勘察、測繪。基于CAD-ANSYS構建原始三維邊坡模型，導入FLAC3D2600軟件后模擬煤層開采對邊坡穩定性的影響。

FLAC3D由美國ITASCA公司研究開發，是一款基于拉格朗日算法有限差分的顯式數值分析模擬軟件。該軟件包含11種彈塑性材料本構模型，非常適合模擬巖土體的屈服、塑性流動、軟化及大變形問題，可以快速求解場的力學方程，準確分析巖土邊坡的穩定性[26]。拉格朗日差分計算原理是基于顯式差分求解偏微分方程[27，28]，即通過運動方程計算節點運動速率，根據運動時間Δt，計算節點(ni)時步位移，由節點位移推算單元(ei)之間的相對位移(單元ei由邊長lj，lj+1，…lj+k構成，邊長lj由節點nk，nk+1和構成)，進而求解單元應變增量(△εi，j)。由單元材料本構方程獲取應變條件下的應力增量(△σi，j)，在此基礎上求解各節點的不平衡力(即合力)，根據中心差分原理[29]，計算節點的速度與加速度，進而再一次獲得各節點的時步位移。按照上述步驟往復迭代，擴展到整個區域，直至各節點的不平衡力足夠小或者各節點的位移趨于平衡[30]。

2 研究區域概況及地質模型

馬達嶺滑坡位于貴州省都勻市江洲鎮富溪村干壩馬達嶺溝，溝口地理坐標為26°10′18.9″N，107°17′23.4″E[22](圖1)。馬達嶺溝流域面積0.86km2，流域全長1.62km，最大高程1570m(北部山脊)，最小高程1150m(南部溝口)。溝道走向大致由北向南，中游呈“V”型，上、下游呈“U”型，溝道平均縱向比降約28.2%，上陡下緩。災害點地處亞熱帶季風濕潤氣候，雨熱同期，4—9月為雨季，年平均降雨量約1400mm，且主要集中在6—8月[22]。災害點地層屬濱海沉積相，地質構造相對穩定，溝口上方有一條東西走向的正斷層。地層地質時代劃歸早石炭世大塘早期，為祥擺組(C1x)。地層結構由灰黃色的石英砂巖、黃褐色砂質頁巖、灰黑色炭質頁巖及黑色煤層組成，其中夾少量泥灰巖。地層主要特征為砂頁巖夾煤層，且被認為是貴州石炭系出露最廣、唯一含煤的巖組。由Ⅰ-Ⅰ′地質剖面(圖2)可知，巖層產狀280°∠18°，緩傾坡內。煤層主要有3層，自下而上分別為M1、M2和M3，其中M1和M2煤層厚度1.9～2.4m，M3煤層厚度0.2～0.9m不等(圖3)。因煤炭開采，坡頂產生縱橫交錯、寬度不等的裂隙，發育形成兩組優勢結構面[20]，大致產狀分別為34°∠87°、118°∠86°，受其影響滑坡體與山體拉裂，形成“圈椅狀”滑坡后壁。滑坡位于馬達嶺坡頂，滑坡體后緣肩部高程約1573.6m，滑坡前緣坡腳高程約1481.5m，相對高差約92.1m(圖3)。滑坡滑動主方向約204.6 °，滑坡體沿該主方向長度約88m，平均寬度102m，平均厚度92m，滑坡體體積達8.26×105m3。此外，滑坡物堆積于馬達嶺溝中部，堆積體平均厚度約6m[23]。

圖1 滑坡區域平面圖(滑坡前)

圖2 Ⅰ-Ⅰ′剖面地層結構

圖3 滑坡全貌圖

由于FLAC3D構建不規則三維模型具有一定的難度，因此本文以Ⅰ-Ⅰ′剖面為基礎，在ANSYS軟件平臺上恢復原始三維邊坡概化模型，之后再導入FLAC3D(圖4)。該邊坡模型x軸方向長0～426.49m，y軸方向寬0～100m，z軸方向高1400.00～1577.09m，模型共產生381786個節點和365950個單元，分19個地層，分別為砂巖、炭質頁巖和煤巖互層。巖性參考文獻[23][31]，本文數值模擬計算時假設各層砂巖的巖性一致，各層炭質頁巖的巖性一致，各層煤巖的巖性也一致(表1)。由于FLAC3D關于巖體變形模量僅接受體積模量K和剪切模量G，因此本文設定砂巖、炭質頁巖及煤巖均彈塑性材料，服從莫爾-庫倫屈服準則，并依據文獻[32]確定體積模量K、剪切模量G與彈性模量E及泊松系數μ的關系。

圖4 初始概化三維邊坡模型

3 煤層開挖過程的邊坡變形分析

為深入分析煤層開挖對邊坡穩定性的影響，本文設計5種工況，分別探討開挖進深、開挖作業面寬度及開挖煤層的先后次序等對馬達嶺邊坡變形的影響。

表1 巖石物理力學參數

3.1 工況設計

1)工況一：開挖M1、M2煤層至x=300m，作業面寬度8m。通過FLAC3D模擬M1、M2煤層開挖過程，并監測邊坡應力應變狀態。M1、M2煤層開挖最終各形成5個工作面，每個工作面寬8m，工作面之間以煤柱作為支撐，其中，M1煤層開挖至x=300m，開挖進深為245m，M2煤層也開挖至x=300m，開挖進深為187m，工作面y軸方向的坐標分別為：4～12m；24～32m；44～52m；64～72m；84～92m，開挖位置及尺寸如圖5(a)所示。模擬過程中，在邊坡體中選取有代表性的點位A(256，0，95)、B(280，0，130)和C(360，0，140)(見圖2)作為監測點，以記錄節點豎向位移、節點合速率等參數隨運行步的變化，其在格網中的節點編號分別為A(GP4104)、B(GP1815)、C(GP1414)。

2)工況二：開挖M1、M2煤層至x=400m，作業面寬度8m。M1煤層開挖至x=400m，開挖進深345m，M2煤層也開挖至x=400m，開挖進深287m，其他條件同開挖進深至x=300m一致，開挖位置及尺寸如圖5(a)所示。

3)工況三：開挖M1、M2煤層至x=400m，作業面寬度增加至12～16m。M1煤層開挖至x=400m，開挖進深為345m，M2煤層也開挖至x=400m，開挖進深為287m，工作面12～16m不等，y軸方向坐標分別為：4～20，24～40，44～60，64～80，84～92，開挖位置及尺寸見圖5(b)各工作面之間以煤柱作為支撐。其他條件同工況二。

4)工況四：在工況二基礎上開挖M3煤層，開挖至x=300m。在工況一基礎上開挖M3煤層，開挖至x=300m，M3煤層工作面在y軸方向分別為：4～12m；24～32m；44～52m；64～72m；84～92m，開挖位置及尺寸如圖5(c)所示。

5)工況五：同時開挖M1、M2、M3煤層，其他開挖條件同工況四，開挖位置及尺寸如圖5(c)所示。

圖5 各工況煤層開挖位置及工作面寬度(m)

3.2 邊坡變形結果分析

3.2.1 開挖進深對邊坡變形的影響

根據工況一和工況二的邊坡變形計算結果(圖6)，當M1、M2煤層開挖進深由300m增加至400m，邊坡整體變形增加了約2倍。開挖之初，坡腳前緣部分最先產生塑性變形，但由于坡腳前緣煤層上方的巖體厚度較薄，較小的勢能和載荷使得坡腳變形相對較小。隨著開挖進深增加，邊坡變形越來越大，且變形最大區域逐步由坡腳前緣向坡頂和邊坡內部轉移。

圖6 不同開挖進深下的邊坡合位移

根據邊坡內部A、B、C節點監測數據，節點產生位移的先后次序分別為A、B、C，且由于節點在y軸方向的水平位移最小，在x軸方向的水平位移次之，而在z軸豎向位移最大，因此本文選取豎向位移作為節點位移監測指標，另選取合速率作為節點速率監測指標。由工況一的邊坡變形和節點豎向位移、合速率計算結果(圖6、圖7)可知，當M1、M2煤層開挖進深由300m增加至400m，節點合速率增加了約2倍。節點豎向位移與節點合速率與開挖進深均呈正比例關系(圖7)，即隨開挖進深的增加而增大。坡頂部位合位移將沿Z軸負向推進，其位移大小在任一Y剖面上沿X軸正向呈遞減分布，即形成不同的合位移區間。因此，邊坡巖體存在合位移差異，導致邊坡垂向節理和裂隙形成并發育，且隨著開挖進深的增加而增大，逐步向邊坡內部推進，直至下切到M1、M2和M3煤層，形成潛在的優勢滑動結構面。

圖7 不同開挖進深下的節點豎向位移與合速率監測

3.2.2 開挖作業面寬度對邊坡變形的影響

根據工況三的邊坡合位移計算可知，當M1、M2煤層的作業面開挖寬度由8m增加至12～16m時(2倍)，邊坡合位移增加了約10倍(圖8)。

圖8 增加作業面開挖寬度下的邊坡合位移

根據邊坡內部A、B、C節點豎向位移監測數據記錄，增加作業面寬度后，A、B節點的豎向位移增大了約11倍，C節點的豎向位移增加了約12.3倍。此外，增加作業面寬度后(2倍)，峰值合速率增大約4倍(圖9)。上述數據說明，邊坡內部節點豎向位移與作業面寬度呈正比關系，即作業面寬度越大，節點豎向位移就越大，即作業面寬度對邊坡變形的影響要遠大于開挖進深。

圖9 增加開挖寬度下的節點豎向位移與合速率監測

由此表明，隨著M1、M2煤層開挖作業面寬度的增加，邊坡合位移、節點豎向位移和節點合速率均呈不同程度增加，且節點豎向位移增幅遠大于節點合速率增幅。

3.2.3 開挖煤層先后次序對邊坡變形的影響

由工況四的邊坡變形計算(圖10a，b)可知，邊坡最終合位移由兩次開挖疊加導致，且開挖M3煤層所產生的合位移主要集中于M3煤層上方巖體，合位移約為單獨開挖M1、M2煤層合位移的10 %～20%，說明深部煤層開挖對邊坡變形的影響要遠大于淺層煤層開采。此外，根據工況五的邊坡合位移計算結果(圖10c)，同時開挖煤層M1、M2和M3(工況五)產生的最大變形區域主要分布于邊坡后緣中下部，其整體變形較分步驟開挖煤層產生的變形(工況四)稍大，節點豎向位移和節點合速率記錄也驗證此現象(圖11)。總體而言，開挖煤層先后次序對邊坡變形存在一定的影響，且同時開挖各煤層對邊坡變形的影響程度較分步驟開挖各煤層要大，但相對于開挖進深和開挖作業面寬度而言，其對邊坡變形的影響程度又相對較小。

圖10 不同開挖次序情況下的邊坡合位移

圖11 不同開挖次序下的節點豎向位移與合速率監測

4 基于數值分析的滑坡機理推定

本文在現場調查和定期的邊坡變形監測(GPS)基礎上，開展基于FLAC3D的數值模擬分析，得到與現場調查、監測相一致的結論。說明本文選取邊坡合位移、節點豎向位移和節點合速率等指標監測邊坡變形過程具有可操作性和實際意義，為揭示邊坡滑動機理和本質特征提供了重要的觀測和量化工具。由此，借助該工具推定馬達嶺滑坡機理和過程如下：

單獨開挖M1、M2煤層時，隨煤層開挖進深的推進，邊坡內部因變形位移差異產生沿x軸方向的若干條垂向優勢結構面(e、f、g、h)(圖12)。該類優勢結構面受煤層開挖作業面寬度控制，隨作業面寬度的增加而拉裂擴張。其向上延伸至坡頂，形成若干條地裂縫，成為坡面雨水入滲邊坡的優勢通道，向下則切割邊坡，直至M1、M2煤層。

圖12 滑坡機理示意圖

在此基礎上若進一步開挖M3煤層，邊坡上方前緣再次遭受采掘擾動。此時位于M3開挖進深范圍內、前期所形成的垂向優勢結構面，將隨M3煤層作業面寬度的增加而進一步拉裂擴張。受M3煤層采空影響，M3煤層上方巖體沿z軸負向產生位移，同時M3下方巖體因采空卸荷而反彈，產生若干卸荷拉張裂隙面，即橫向優勢結構面(d的下部)。在滿足一定的降雨條件時，M3煤層采空區匯集大量裂隙水，一部分裂隙水沿M3煤層開挖工作面入口溢出，另一部分則沿橫向優勢結構面繼續下滲。當M3煤層上方巖體發生失穩垮塌，在重力作用下沖擊M3煤層下方巖體時，橫向優勢結構面鎖固段被打通，優勢結構面(d)在水、巖等綜合作用下實現連續和貫通。此時，上方巖體沿優勢結構面在剖面D處臨空剪出，剪出巖體受重力作用拋向馬達嶺溝谷，遭受撞擊后形成大小不等塊體散落于溝谷中段，滑坡地質災害由此發生。

馬達嶺滑坡地質災害發生后，該邊坡災害隱患依然存在，甚至更為嚴峻。一方面，拋落在溝谷中的大小巖塊堆積物在極端天氣條件下可能會再次發生滑動，威脅溝谷下游安全。據2018年6月22日強降雨后的現場觀測記錄數據，該原滑坡堆積區(D下方)邊緣產生5～10cm寬度不等的蠕滑裂縫。另一方面，開挖M1、M2過程中所產生的一系列垂向優勢結構面e、f、g、h等，已經將邊坡切割為若干部分，這些被切割巖體之所以能處于相對穩定狀態，是因為邊坡巖層反傾向于坡內，并缺少足夠水文地質條件。然而，該邊坡及后緣山體巖層產狀為NW282°∠18°，在水、巖長期相互作用下，山體垂向優勢結構面和巖層水平節理、裂隙會持續擴張。滿足一定降雨條件后，地表水將沿后緣山體垂向優勢結構面滲入巖體層狀節理、裂縫，引起巖層順層錯動，直接擠壓和推動前述邊坡的被切割巖體沿優勢結構面(如：g-M1)剪出，到時可能會形成規模更大的滑坡地質災害。

5 結論與討論

5.1 結論

1)將基于快速拉格朗日差分計算的FLAC3D軟件應用于貴州省都勻市馬達嶺滑坡機理分析，通過記錄邊坡變形最大不平衡力、單元和節點的豎向位移、合速率等參數變化來揭示邊坡滑動機理和過程，具有明顯的技術比較優勢，可以在西南喀斯特崩塌和滑坡等地質災害研究中進一步推廣和應用。

2)煤層開挖進深控制邊坡垂向優勢結構面的空間分布，煤層開挖作業面寬度控制垂向優勢結構面擴張程度，后者對邊坡穩定性的影響遠大于前者。

3)深層煤層開挖對邊坡變形的影響遠大于淺層煤層開挖。此外，煤層開挖的先后次序對邊坡變形存在一定的影響，但相對于開挖進深和開挖作業面寬度而言，其影響較小。

4)馬達嶺滑坡的形成與發育可以分為兩個階段，第一階段培育了垂向優勢結構面(開挖M1、M2煤層)和橫向優勢結構面(開挖M3煤層)；第二個階段(M3上層巖體重力失穩垮塌，水、巖相互作用)實現橫向、垂向優勢結構面連續和貫通，在水、巖等綜合作用下，巖體沿優勢結構面臨空剪出。

5.2 討論

1)貴州省煤炭資源廣泛分布，有些礦區就在旅游景區周邊，甚至位于景區必經之路旁，對旅游景區和游客安全構成極大威脅。正如馬達嶺坡頂西側山體正在開發農業觀光項目，建議相關部門應該針對此項風險進行專門評估。

2)開挖M1、M2煤層所培育的大量垂向優勢結構面與M1、M2煤層采空區交會，實現了對后緣邊坡的完全切割。不利的巖層產狀(NW282°∠18°)及大量節理裂隙水的存在，可能正在孕育一個規模更大的滑坡地質災害。為此，將在后期通過持續的調查和監測來修正和完善上述結論。

3)基于FLAC3D的邊坡變形過程分析是否適用于含溶洞的喀斯特邊坡穩定性分析，以及是否適用于喀斯特溶洞結構安全評價，尚無相關研究報道，值得進一步深入研究。此外，FLAC3D在這一領域的應用，其可靠性依賴于現場試驗技術和手段的先進性，特別是巖體微細觀節理、裂隙等空間探測技術的重大突破。