摘" "要：滑坡給人類生命、財產(chǎn)和基礎(chǔ)設(shè)施帶來了極高風(fēng)險。在阿拉爾村附近發(fā)現(xiàn)了一處潛在不穩(wěn)定滑坡。本文通過ArcGIS軟件對潛在滑坡進行高程、坡度和坡向3方面分析，依據(jù)這3個地質(zhì)因子對邊坡進行了地質(zhì)分區(qū)；再運用顆粒流軟件PFC2D建立潛在滑坡數(shù)值模型，模擬潛在不穩(wěn)定滑坡位移過程，分析其顆粒應(yīng)變規(guī)律，通過雙軸試驗?zāi)Ｐ蛯峦馏w的細觀參數(shù)進行標定，使其和實際相符；最后將數(shù)值模型的位移結(jié)果與地質(zhì)分區(qū)相結(jié)合，得到該潛在滑坡的位移分析結(jié)果。結(jié)果表明：該潛在滑坡的位移主要集中在高高程、中等坡度區(qū)域，而該滑坡在此區(qū)域面積占比最大，研究為后續(xù)邊坡治理提供了一定理論依據(jù)。

關(guān)鍵詞：PFC2D顆粒流；ArcGIS；潛在滑坡；地質(zhì)分區(qū)；位移特征

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

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

1" 滑坡概況

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

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

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

2" 潛在滑坡的地質(zhì)分區(qū)

2.1" 分區(qū)因子

據(jù)現(xiàn)場地質(zhì)調(diào)查及地表監(jiān)測點分布統(tǒng)計，結(jié)合相關(guān)單點滑坡地質(zhì)分區(qū)文獻，對于單一滑坡地質(zhì)地層構(gòu)造、降雨及人類工程活動等未能體現(xiàn)影響范圍的情況下，參考選取變化較為明顯的高程、坡度、坡向等要素作為地質(zhì)分區(qū)影響因子[12]，利用ArcGIS空間3D分析功能對滑坡進行幾何特征分析，選取高程、坡度、坡向作為邊坡地形因子，分析其空間分布特征。高程是影響滑坡災(zāi)害發(fā)生的重要因素，坡體內(nèi)的應(yīng)力值往往受高程影響而發(fā)生變化。據(jù)統(tǒng)計數(shù)據(jù)，依托工程潛在滑坡高程范圍為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個區(qū)間。該斜坡坡度范圍為0°～57°，同理將坡度劃分為0°～8°、8°～20°、20°～29°、29°～37°、37°～45°、45°～57° 6個區(qū)間。由于坡向決定了光照時間、水熱比、降水分布等對滑坡的影響，將坡向劃分為平面、N向、NE向、E向、SE向、S向、SW向、W向、NW向、N向10個區(qū)間（圖4），再通過柵格計算器計算出每個地質(zhì)因子區(qū)間面積大小（表1）。

2.2" 地形因子疊加特征

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

3" 邊坡運動過程數(shù)值模擬

3.1" 邊坡模型細觀參數(shù)標定

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

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

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

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

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

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

3.2" 邊坡模型構(gòu)建

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

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

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

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

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

運用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個監(jiān)測位置點的監(jiān)測預(yù)警數(shù)據(jù)顯示（圖14）：位于滑坡后緣高高程區(qū)地表位移監(jiān)測點15-GP01最大水平位移達58.4 mm，垂直位移最大7.8 mm，位于滑坡中高程區(qū)的地表位移監(jiān)測點17-GP01最大水平位移達45.1 mm，垂直位移最大-22 mm，位于滑坡低高程區(qū)的監(jiān)測點60-GP01最大水平位移14.9 mm，垂直位移-55 mm。從監(jiān)測數(shù)據(jù)對比分析發(fā)現(xiàn)，該滑坡變形主要集中于滑坡后緣高高程區(qū)，但在滑坡前部低高程區(qū)垂直位移量較高高程區(qū)大，說明滑坡為推移式滑坡，且垂直位移隨降雨量的累計發(fā)生呈明顯增加趨勢，結(jié)果與模擬數(shù)據(jù)一致。

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

4" 結(jié)論

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

（2） 使用試錯法對邊坡的巖土體細觀參數(shù)進行標定，通過雙軸試驗?zāi)Ｐ蛥?shù)計算巖土體宏觀參數(shù)模擬值，并與室內(nèi)試驗實測值對比，最終計算數(shù)據(jù)誤差值均在合理范圍內(nèi)，符合滑坡實際情況。

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

參考文獻

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

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

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

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

[5] 楊博，田文通，孫軍杰，等.地震作用下高路塹黃土邊坡細觀力學(xué)研究[J].地震工程學(xué)報，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].山地學(xué)報，2020，38（3）：402-415.

[9] 張志東，樊曉一，姜元俊.滑源區(qū)粒序分布及顆粒粒徑對碎屑流沖擊作用的影響研究[J].水文地質(zhì)工程地質(zhì)，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] 謝良甫.基于地質(zhì)幾何分區(qū)的反傾邊坡傾倒變形特征[J].科學(xué)技術(shù)與工程，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] 吳順川，張曉平，劉洋.基于顆粒元模擬的含軟弱夾層類土質(zhì)邊坡變形破壞過程分析[J].巖土力學(xué)，2008（11）：2899-2904.

[15] 陳建峰，李輝利，周健.黏性土宏細觀參數(shù)相關(guān)性研究[J].力學(xué)季刊，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的黃土三軸試驗細觀參數(shù)敏感性分析[J].工程地質(zhì)學(xué)報，2021，29（5）：1342-1353.

[18] 孟云偉，柴賀軍.顆粒流離散元在滑坡運動過程模擬中的應(yīng)用[J].巖土力學(xué)，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