基于地質分區的伊犁河谷黃土滑坡滑移特征研究

摘" "要：滑坡給人類生命、財產和基礎設施帶來了極高風險。在阿拉爾村附近發現了一處潛在不穩定滑坡。本文通過ArcGIS軟件對潛在滑坡進行高程、坡度和坡向3方面分析，依據這3個地質因子對邊坡進行了地質分區；再運用顆粒流軟件PFC2D建立潛在滑坡數值模型，模擬潛在不穩定滑坡位移過程，分析其顆粒應變規律，通過雙軸試驗模型對滑坡土體的細觀參數進行標定，使其和實際相符；最后將數值模型的位移結果與地質分區相結合，得到該潛在滑坡的位移分析結果。結果表明：該潛在滑坡的位移主要集中在高高程、中等坡度區域，而該滑坡在此區域面積占比最大，研究為后續邊坡治理提供了一定理論依據。

關鍵詞：PFC2D顆粒流；ArcGIS；潛在滑坡；地質分區；位移特征

滑坡的發生是在重力作用下沿某一軟弱滑移面（帶）整體或分散下滑而形成，其往往受河流沖刷、地下水位變化、雨水浸潤、地震活動及人工切坡等因素的影響。位于天山西段的伊犁河谷鞏留縣發育大量的潛在黃土邊坡，滑坡、崩塌等，地質災害時有發生，如1960年6月發生的鐵列克薩溝頭滑坡導致100畝農田被淹沒，尤其在2003年前后，縣域東部山區發生的諸多地質災害造成大量人員傷亡和直接經濟損失。為保護當地居民人身和經濟財產安全，對當地潛在黃土滑坡進行土體位移分析刻不容緩，了解分析潛在滑坡運動規律及運動趨勢對于災害的預防起到很大作用。

滑坡運動過程模擬避免了傳統試驗規模和檢測方法的限制，能夠較為詳細的觀察滑坡變形破壞過程，是滑坡災變風險定量評估的重要研究方法之一[1]。Zhou Jian等嘗試在離散單元法中引入強度折減法和重力增加法來評價土坡安全系數[2]。趙洲等將顆粒流離散元方法應用于堆積體滑坡破壞模式和運動過程的數值模擬研究[3]。李子隆等運用經驗法和PFC2D顆粒流方法對滑坡失穩造成的堵江高度進行了預測[4]。Yang Bo等在室內外勘查和試驗基礎上采用PFC2D建立數值模型[5]，模擬了邊坡地震響應全過程。Zhou Jiawen等采用離散元法（黏結-接觸模型）研究了地震引起的沉積顆粒剛度、粘結力、摩擦系數和圍壓應力對其力學特性的影響及數值參數與實驗參數之間的關系[6]。Zhang Qingzhao等通過室內滑槽試驗與PFC數值模擬相結合的方法再現了滑坡-泥石流在山谷中的運動過程及其堆積模式[7]。了解滑坡堆積物的物理力學性質對研究滑坡的破壞機制和泥石流啟動機理具有重要意義，對于使用離散元法進行滑坡的穩定性分析，國內外已有眾多學者進行了大量研究[8-11]，此方法對于滑坡分析具有很大意義，但大多數研究僅立足于對滑坡土體數值模擬分析，未與滑坡實際地質因子相結合，而實際地質因子對于潛在滑坡運動滑移有重要的影響。因此，本文在野外實地調查基礎上，利用ArcGIS軟件對潛在滑坡進行地質分區，再結合顆粒流軟件PFC2D對潛在滑坡進行數值模擬分析，通過檢測顆粒位移和顆粒應力變化與地質分區結合，將滑坡分為強變形區、中變形區和輕變形區，研究結果為該地區潛在滑坡治理提供了理論依據。

1" 滑坡概況

阿拉爾潛在滑坡位于伊犁河谷鞏留縣阿拉爾村古滑坡后壁中部，該斜坡高152 m，坡度25°～35°，坡向34°，整體呈上陡下緩形態。該斜坡山頂處分布有長535 m的近NS向寬弧形張拉裂縫，裂縫寬5～50 cm，垂直方向位移10～60 cm，坡體中部變形較大，現場調查穩定性較差，嚴重威脅坡腳村民生命財產安全。該潛在滑坡屬降雨型黃土淺層滑坡（圖1），位于陰坡坡面，給水滲透時間較長，裂縫發育，斜坡穩定性差。

據伊犁地區地質災害易發分區結果，潛在滑坡易發區內大部分斜坡被上更新統風成黃土覆蓋，厚度從幾米至幾十米不等，顆粒成分以粉砂為主，礦物成分以石英和長石為主，黏土礦物較少。由于該類黃土孔隙大，礦物顆粒吸水性和持水性較強，在春季融雪和夏季降雨作用下土體易飽和。據勘查資料顯示，該滑坡黃土厚度40～60 m以上，推斷該滑坡為單一均質黃土滑坡。由于淺層黃土為區內主要易滑層，黃土下伏巖層以第三系泥巖為主，局部夾砂礫石，在融雪或降雨作用下，接觸面附近的巖土體物理力學性能惡化，常隨黃土沿軟弱面產生滑動。潛在滑坡典型剖面見圖2。

自2009年起阿拉爾潛在滑坡實施地質災害監測預警示范建設，先后累計建設雨量計、土層含水率監測儀、裂縫伸縮儀、裂縫位移計、深部變形監測儀等監測設備67臺套。據鞏留縣示范站監測數據顯示，該滑坡年均絕對位移累積沉降量為4 mm，單次最大沉降變幅22 mm，X方向最大位移量變幅173 mm，單次位移量最大變幅50 mm，裂縫累計最大位移為滑坡西北側位移計476 mm，滑坡深部變形主要集中發生于9—10月，土壤含水量監測數據顯示深部距地表40 cm處變化最明顯，其他深度變化幅度較小。阿拉爾滑坡示范站監測滑坡位移和土壤含水率變化情況見圖3。

2" 潛在滑坡的地質分區

2.1" 分區因子

據現場地質調查及地表監測點分布統計，結合相關單點滑坡地質分區文獻，對于單一滑坡地質地層構造、降雨及人類工程活動等未能體現影響范圍的情況下，參考選取變化較為明顯的高程、坡度、坡向等要素作為地質分區影響因子[12]，利用ArcGIS空間3D分析功能對滑坡進行幾何特征分析，選取高程、坡度、坡向作為邊坡地形因子，分析其空間分布特征。高程是影響滑坡災害發生的重要因素，坡體內的應力值往往受高程影響而發生變化。據統計數據，依托工程潛在滑坡高程范圍為1 187～1 517 m，基于Jecks將滑坡高程分成（1 187，1 228）m、（1 228，1 257）m、（1 257，1 286）m、（1 286，1 316）m、（1 316，1 346）m、（1 346，1 380）m、（1 380，1 415）m、（1 415，1 446）m、（1 446，1 479）m、（1 479，1 517）m 10個區間。該斜坡坡度范圍為0°～57°，同理將坡度劃分為0°～8°、8°～20°、20°～29°、29°～37°、37°～45°、45°～57° 6個區間。由于坡向決定了光照時間、水熱比、降水分布等對滑坡的影響，將坡向劃分為平面、N向、NE向、E向、SE向、S向、SW向、W向、NW向、N向10個區間（圖4），再通過柵格計算器計算出每個地質因子區間面積大小（表1）。

2.2" 地形因子疊加特征

據表2中給出的標準，對邊坡高程、坡度、坡向柵格圖層進行重分類，通過圖層疊加得到邊坡地質幾何特征分區圖。最終將邊坡分為132個特征區，總面積158 446.5 m2（圖5）。統計圖5中面積大于2 000 m2的特征區占岸坡中總面積的56％。從圖6可看出，面積最大的特征區為532號，面積為12 043 m2，即中等坡度，高高程，北坡向區。

3" 邊坡運動過程數值模擬

3.1" 邊坡模型細觀參數標定

顆粒流軟件的參數采用顆粒間接觸模型的細觀參數，在實際的巖土體中并不存在該參數，但兩者間存在一定聯系，目前尚無成熟理論和方法確定宏觀參數和細觀參數間的關系，目前最常用的方法是試錯法，本次模擬即采用試錯法確定細觀參數。降雨對滑坡的影響主要表現為降雨入滲導致斜坡土體飽和，甚至在下伏隔水層積水，增加滑體的重量，降低土體抗剪強度而發生滑坡。所以本模型是在飽和土基礎上考慮土的宏觀參數建立數值模型。顆粒接觸模型為PFC2D內置的接觸模型并聯粘結模型。基于Potyondy和Cundall等人提出的線性模型[13]，Pb模型（平行鍵）由無活性的粘合點和參考間隙為零的平行結組成。線性平行粘結模型提供了2個界面的行為：①一個無窮小的線性彈性（無張力）和承載一個力的摩擦界面；②一個有限大小的線性彈性和承載一個力和力矩的粘結界面（圖7，8）。

形成平行粘結后，力[Fi]和彎矩[M3]的初始值為零。每個計算時間步在接觸點處產生相應的位移和轉動增量，可通過方程（1）轉換為力和力矩增量[14]：

[ΔFni=（-knAΔUn）ni;ΔFsi=-ksAΔUsi;ΔM3=-knIΔθ3]（1）

式（1）中：[kn]、[ks]分別為并聯連接的法向剛度和切向剛度；[ni]為接觸點的法向量；[A]為并聯連接的截面面積；[I]為截面沿接觸點相對于轉動方向的轉動慣量。將力和力矩的增量分別加入到它們的初始值中，得到新的力和力矩：

[Fni←Fnni+ΔFni;Fsi←Fsi+Fsi;M3←M3+ΔM3]（2）

本模型使用雙軸試驗進行宏細觀參數標定（圖9），首先根據已有巖土體宏細觀參數相關性研究經驗給出潛在滑坡巖土體的初始細觀參數[15-17]，然后根據滑坡模型的顆粒接觸模型修改調整雙軸試驗模型，并對初始細觀參數進行不斷試錯調整，以獲取對應數值的宏觀力學參數，使得試驗所得宏觀力學參數與室內土工試驗所得參數相近，得到潛在滑坡數值模擬所需細觀參數（表3）。

3.2" 邊坡模型構建

顆粒流離散元是離散單元法的一種，離散元法的思想源自分子動力學[18]。Cundall于20世紀90年代末期首次提出了適用于解決巖石力學問題的顆粒流離散元PFC。基于PFC建立滑坡模型主要有Ball-Ball法和Ball-Wall法兩種。其中Ball-Ball法是用顆粒構建滑體、滑床及邊界，然后依據不同巖性給對應位置的顆粒賦予不同屬性，Ball-Wall法是僅用顆粒組建滑體，用wall來構建滑床及邊界。本文采用Ball-Ball法建立滑坡模型。模型建立步驟：①將滑坡剖面圖導入CAD生成dxf格式文件，通過geometry import命令導入PFC2D；②在特定幾何形狀中生成均勻顆粒后再刪除邊坡表面約束顆粒的墻體；③使用set gravity命令賦予顆粒y方向上的重力加速度，讓顆粒在自重作用下產生滑動；④使用measure create命令在邊坡表面生成5個測量圓，通過測量圓檢測顆粒應力變化規律和顆粒位移反映邊坡位移特征。

根據上述建模步驟生成長16 m、高13.4 m的滑坡模型（圖10）。模型顆粒總數13 905個，其中黃土顆粒約5 550個，粒徑0.02～0.05 m；砂巖顆粒約8 355個，粒徑0.05～0.08 m；并在潛在變形區特征位置設置5個一定半徑的測量圓。

3.3" 邊坡運動過程模擬及分析

等初始顆粒達到平衡后賦予滑坡顆粒粘聚力等細觀參數，并刪除坡面與巖土體接觸的墻體，同時用set gravity 9.81命令給模型施加重力作用，使滑坡顆粒在重力作用下產生滑移。滑坡在重力的作用下，運行到10 000步時坡面開始發生變形（圖11-a），位移大小隨滑坡高程增加而增加。運行至40 000步時位移主要集中在滑坡表層中高高程且高坡度處（圖11-b），坡腳處顆粒由于重力作用向下位移，最大位移達0.19 m。模型運行至55 000步時（圖11-c），坡底位移基本為0，進入穩定狀態。位移基本集中于中高高程且高坡度處。模型運行140 000步時（圖11-d），基本達到穩定狀態，最大位移達4.8 m。

運用設置好的測量圓對坡頂、坡腰和坡腳的x和y方向的應力進行監測（圖12）。從圖12-a可看出，1號測量圓位置初始階段顆粒向下發生位移，使x和y方向應力均大幅增加，且 y方向應力增加更大，約為x方向1.5倍，在0～0.5 s內達到應力峰值。2號測量圓也處于靠近坡頂位置（圖12-b），所以應力曲線波動和1號相似。由于坡度較緩，分別位于坡中和坡腳的3、4和5號測量圓內（圖12-c-e），應力均于約0.2 s達到應力峰值后趨于基本平衡狀態。

運用5個測量圓檢測坐標位置為A（-6，6），B（-3.8，4.3），C（-0.8，2），D（-1.3，1.5）和E（-4.8，1.2）的5個顆粒位移隨時間的變化，其中B顆粒位移最大且達速度峰值所用時間最少（圖13）。

選取阿拉爾滑坡前中后部3個監測位置點的監測預警數據顯示（圖14）：位于滑坡后緣高高程區地表位移監測點15-GP01最大水平位移達58.4 mm，垂直位移最大7.8 mm，位于滑坡中高程區的地表位移監測點17-GP01最大水平位移達45.1 mm，垂直位移最大-22 mm，位于滑坡低高程區的監測點60-GP01最大水平位移14.9 mm，垂直位移-55 mm。從監測數據對比分析發現，該滑坡變形主要集中于滑坡后緣高高程區，但在滑坡前部低高程區垂直位移量較高高程區大，說明滑坡為推移式滑坡，且垂直位移隨降雨量的累計發生呈明顯增加趨勢，結果與模擬數據一致。

綜上所述，阿拉爾潛在滑坡破壞模式為：在降雨條件下土體顆粒達到飽和狀態，強度降低后產生蠕滑變形，隨著土體顆粒強度不斷變小，蠕滑變形量隨時間增加且坡體中等坡度、高高程區域位移明顯。結合地質因子分區結果可知，該滑坡的中等坡度、高高程區域占比面積最大，所以對該潛在滑坡區域應多加關注并及時防護。

4" 結論

（1） 根據現場勘查數據，利用阿拉爾潛在滑坡的高程、坡度和坡向等地質因子對邊坡進行地質分區，通過柵格計算器計算出占比較高的幾個分區面積。

（2） 使用試錯法對邊坡的巖土體細觀參數進行標定，通過雙軸試驗模型參數計算巖土體宏觀參數模擬值，并與室內試驗實測值對比，最終計算數據誤差值均在合理范圍內，符合滑坡實際情況。

（3） 通過顆粒流PFC2D軟件對邊坡位移過程進行模擬分析，認為該潛在滑坡在初始階段位移隨高程增大而增大；在中后期位移主要集中在滑坡中等坡度、高高程區域。

參考文獻

[1] 曹文，李維朝，唐斌，等.PFC滑坡模擬二、三維建模方法研究[J].工程地質學報，2017，25（2）：455-462.

[2] 周健，王家全，曾遠，等.顆粒流強度折減法和重力增加法的邊坡安全系數研究[J].巖土力學，2009，30（6）：1549-1554.

[3] 趙洲，魏江波.基于顆粒流方法的堆積層滑坡運動過程模擬[J].煤田地質與勘探，2017，45（6）：111-116+122.

[4] 李子隆，王浪.基于PFC2D的山區中小型水庫滑坡堵江模擬研究[J].人民珠江，2019，45（6）：111-116+122.

[5] 楊博，田文通，孫軍杰，等.地震作用下高路塹黃土邊坡細觀力學研究[J].地震工程學報，2021，43（1）：187-194.

[6] Jiawen Z，Peng C，Xingguo Y.Effects of material composition and water content on the mechanical properties of landslide deposits triggered by the wenchuan earthquake[J].Acta Geologica Sinica，2016，90（1）：242-257.

[7] Zhang Q-Z，Pan Q，Chen Y，et al.Characteristics of landslide-debris flow accumulation in mountainous areas[J].Heliyon，2019，5（9）：1-10.

[8] 張志東，樊曉一，姜元俊. 巖土體顆粒級配對滑坡碎屑流沖擊力力鏈特征的影響[J].山地學報，2020，38（3）：402-415.

[9] 張志東，樊曉一，姜元俊.滑源區粒序分布及顆粒粒徑對碎屑流沖擊作用的影響研究[J].水文地質工程地質，2021，48（1）：49-59.

[10] Zhuang J，Jia K，Zhan J，et al.Scenario simulation of the geohazard dynamic process of large-scale landslides：a case study of the Xiaomojiu landslide along the Jinsha River[J].Natural Hazards，2022，112（2）：1337-1357.

[11] Zou Z，Tang H，Xiong C，et al.Kinetic characteristics of debris flows as exemplified by field investigations and discrete element simulation of the catastrophic Jiweishan rockslide，China[J].Geomorphology，2017，295：1-15.

[12] 謝良甫.基于地質幾何分區的反傾邊坡傾倒變形特征[J].科學技術與工程，2018，18（12）：32-37.

[13] Potyondy D O，Cundall P A.A bonded-particle model for rock[J].International Journal of Rock Mechanics and Mining Sciences，2004，41（8）：1329-1364.

[14] 吳順川，張曉平，劉洋.基于顆粒元模擬的含軟弱夾層類土質邊坡變形破壞過程分析[J].巖土力學，2008（11）：2899-2904.

[15] 陳建峰，李輝利，周健.黏性土宏細觀參數相關性研究[J].力學季刊，2010，31（2）：304-309.

[16] Li B，Gong W，Tang H，et al.Probabilistic analysis of a discrete element modelling of the runout behavior of the Jiweishan landslide[J].International Journal for Numerical and Analytical Methods in Geomechanics，2021，45（8）：1120-1138.

[17] 許江波，曹寶花，余洋林.基于PFC3D的黃土三軸試驗細觀參數敏感性分析[J].工程地質學報，2021，29（5）：1342-1353.

[18] 孟云偉，柴賀軍.顆粒流離散元在滑坡運動過程模擬中的應用[J].巖土力學，2006（3）：559-561.

Study on Characteristics of Loess Landslide and Slip in Yili River

Valley Based on Geological Zoning

Mao Wei1，2， Wang Yihui2， Ruheiyan·Muhemaier1，2， Xie Liangfu2， Liu Xuejun3， Mayila·Aishan4

（1.Xinjiang Vocational and Technical College of Communication，Urumqi，Xinjiang，831401，China;

2.Architectural Engineering Institute，Xinjiang University，Urumqi，Xinjiang，830000，China;3.Xinjiang

Academy of Building Research Co.，Ltd.，Urumqi，Xinjiang，830000，China;4.School of Geology and Mining

Engineering，Xinjiang University，Urumqi，Xinjiang，834000，China）

Abstract： Landslides pose serious risks to human life， property and infrastructure. Near the village of Alaer， a potentially unstable landslide was discovered. In this paper， ArcGIS software is used to analyze the elevation， slope and slope aspect of the potential landslide， and the slope is partitioned according to these three geological factors; then， the particle flow software PFC2D is used to establish the numerical model of potential landslide， simulate the displacement process of potential unstable landslide and analyze its particle stress law. The biaxial test model is used to calibrate the Meso-structure parameters of landslide soil to make it consistent with the actual situation. Finally， the displacement analysis results of the potential landslide are obtained by combining the displacement results of the numerical model with the geological zoning. The results show that the displacement of the potential landslide is mainly concentrated in the area of high elevation and moderate slope， and the area of this landslide occupies the largest proportion. This study provides a certain theoretical basis for the subsequent slope treatment.

Key words： PFC2D; ArcGIS; Potential landslide; Geological zoning; Characteristic of displacement