基于無人機調查的露天礦滑坡運動特征數值研究

何旭乾，盧才武，李 萌，陳彥鈺，閆雪頌

(西安建筑科技大學資源工程學院，陜西 西安 710055)

露天礦山在開采或自然條件等因素影響下，一些邊坡會發生失穩，從而引發滑坡，對礦山開采和人員生命安全造成嚴重威脅。因此，對滑坡風險定量分析進行研究，對于滑坡災害的預防治理具有重要意義。隨著計算機性能的不斷進步，各種數值模擬軟件由于自身低成本、兼容強等優點，在滑坡災害模擬與致因分析過程中應用廣泛。顆粒流(particle flow code，PFC)程序是基于CUNDALL等[1]提出的離散元理論，通過Itasca開發的二維模擬程序或三維模擬程序。TANG等[2-3]應用PCD2D分別對1941年和1999年的滑坡的運動特性、裂紋擴展、斷裂特性進行分析；LO等[4]使用PCD3D對導致滑坡和破壞的運動學過程進行模擬，評估了滑坡機理、泥石流運動、沉積程度。由此可知，PFC軟件能夠模擬剛性粒子的運動和相互作用，允許其旋轉和移動，適用于滑坡形變破壞、滑移運動、堆積特征的分析。為了反映真實的滑坡運動，曹文等[5]指出了wall-ball法具備顆粒少、節省運算時間、結果準確度高的優點，但是由于PFC軟件在滑坡建模這一前期處理方面不具備優勢，因而受到一定的局限。

當今無人機技術日趨成熟，其具備的高精度、高效率、低成本、低風險等優點能夠很好地解決PFC軟件的不足。MA等[6]使用無人機劃分滑坡源區、重疊區和堆積區，結合DEM估算滑坡體積并分析DEM運動過程，定量描述了滑坡特征。因此，通過無人機建立高精度數字高程模型(digital elevation model，DEM)，能為PFC軟件的滑坡運動模擬提供準確的三維離散元模型。

綜上所述，本文對于滑坡運動過程離散元數值模擬主要流程為：首先通過無人機航攝影像建立滑坡區DEM并復現原始地形，構建出滑坡地形(wall)；然后以三軸壓縮試驗與BP神經網絡實現細觀參數的選取，對顆粒(ball)進行參數賦值；最后構建完整的ball-wall離散元三維滑坡模型，并從整體和局部的角度對滑坡區域的運動過程、堆積特征進行反演與分析，為滑坡的定量評估提供一套快捷、有效、完整的解決思路。

1 研究區概述

1.1 露天礦區基本工程概況

本文選取位于陜西省境內某露天煤礦作為研究對象(圖1)，該露天煤礦地處黃土高原，大部分地區為風沙灘地地貌，總體地勢平坦，海拔為1 060～1 230 m。當前采場非工作幫整體邊坡角為24°～32°，南端幫為15°，北端幫為17°。坡體巖性由上至下分別為松散沙層、粉質黏土、砂巖，粉質黏土顏色呈淺紅色或棕紅色，各種中砂、細砂、粉砂為淺黃色或褐黃色。雖然整體邊坡處于穩定狀態，但是根據現場勘查和對地質資料分析，部分位置在風化剝蝕作用下仍出現了局部滑坡現象。

圖1 露天礦現場情況與巖性分布Fig.1 Open-pit mine site situation and lithology distribution

1.2 滑坡現場地質情況

根據對采場現場的調查與分析，造成滑坡的主要影響因素包括節理裂隙、邊坡角設計、邊坡暴露時長等。該區域內黃土層和紅土層垂直節理發育，自然風化會引起局部淺層滑坡，邊坡巖體強度隨著時間增長呈下降趨勢，易產生滑坡。

該滑坡位于采場西北方向非工作幫的地表臺階上，臺階高度為20 m，整體地形主要由粉質黏土構成，采場區內土層垂直節理發育，經過長期暴露，邊坡巖體隨時間增長出現應力松弛、強度下降的現象，部分坡體向下滑落產生滑坡并形成堆積，根據滑坡后堆積DEM可知(圖2)，主要滑動幾乎是朝向東南方，滑移區域高約19 m，平均寬度約17 m，滑坡由頂端開始滑落，最終堆積層位于底端。

圖2 滑坡區域Fig.2 Landslide area

2 數值模型建立及參數確定

2.1 露天礦邊坡三維數值模型構建

本文使用大疆精靈4 Pro無人機對露天礦區進行航測，基于無人機影像數據，使用Agisoft Photoscan生成分辨率礦區的DEM數據。通過Global Mapper對進行影像裁剪處理，選出不穩定邊坡區域，并導入Rhinoceros軟件中，構建STL網格滑坡模型。在PFC3D顆粒流程序中將模型作為Geometry進行導入，Wall邊界設置為模型邊界。

在對滑坡區域進行反演的過程中需要對原始地形進行復現。根據DUMAN[7]的研究表明，通過對滑坡體積的估算來建立原始地形的方法具有可行性。在Rhinoceros軟件中導入滑坡模型后，生成滑動后地形等高線。通過參考ALOS PALSAR 12.5 m的DEM數據來生成滑坡前地形等高線，在Rhinoceros軟件根據等高線劃分滑移區域。由圖3可知，滑坡區域地形發生了明顯改變，滑坡區域呈現下窄上寬，上部寬約23.3 m，下部寬約7.0 m，相對高差19.3 m，滑坡面積約360 m2，體積約為207 m3，1 150 m高程以下的等高線發生明顯的彎曲，根據實際調查情況可以看出其下方形成了堆積。

圖3 滑坡區域滑移前后的高程線Fig.3 Elevation lines before and after slippage in the landslide area

由于滑移區和堆積區存在重疊，因此在完成地形重建后，可以看出滑坡的頂端較陡，滑移區呈現倒梯形，邊坡滑動面較淺，滑移區與堆積區有部分重疊，滑坡頂端與堆積體高差為14 m(圖4)。

圖4 滑坡區域三維地形Fig.4 3D topography of landslide area

2.2 巖土體的細觀力學參數確定

在PFC3D顆粒流軟件建立滑坡模型過程中，需要對模型賦予相應的參數，由于材料的細觀參數與實際巖土體宏觀力學參數沒有直接的量化關系式，為了使得模擬結果接近實際情況，需要通過三軸壓縮試驗，確定合適的細觀參數來準確反映實際的巖石力學特征。一般來說，這些細觀參數包括顆粒彈性模量Ec、顆粒法向剛度kn、顆粒切向剛度ks、摩擦系數μ、法向黏結強度σ、切向黏結強度c，這些參數對彈性模量E、泊松比?、黏聚力等宏觀力學參數具有較大影響。

三軸壓縮試驗的目的是對巖土體的宏細觀參數進行標定，以此獲得實驗數據。在數值模擬實驗中，首先根據以往對巖土體宏細觀參數的研究經驗，設置初始的細觀參數，然后在PFC3D軟件內建立三軸試驗模型，設置符合滑坡模型的顆粒大小和接觸模型，并對初始細觀參數進行模擬。通過對應變曲線的分析，可以獲取到對應數值的宏觀力學參數，將與實際試驗結果相近的參數值作為輸出樣本。

該滑坡體主要為粉質黏土，破壞誘發因素為自然風化與長期暴露。首先，按照巖樣尺寸生成邊界墻體，依照墻體生成緊密顆粒，指定最大顆粒與最小顆粒的半徑比，完成設定后PFC3D系統會按照預設值，將最大值和最小值隨機均勻附加在模型的每個顆粒上，最終壓縮巖樣如圖5所示。

圖5 三軸試驗數值模擬模型Fig.5 Numerical simulation model of triaxial test

為保證生成模型的顆粒尺寸和分布滿足實驗要求，使得顆粒之間具備較高的接觸精度，需要對模型進行均勻化處理。均勻化處理的目的是要保證生成的各相鄰顆粒之間的重疊度足夠小，讓顆粒與邊界接觸緊密，達到完全耦合的效果。根據圖6的均勻化處理，利用自編程序不斷調整模型邊界，使得顆粒之間接觸處于足夠理想的狀態。

圖6 三軸試驗均勻化處理Fig.6 Homogenization of triaxial test

接觸模型選擇為平行黏結模型，該模型可以同時傳遞力和力矩，受到破壞后剛度下降，與滑坡破壞特性相近。完成接觸模型設定后，將一組細觀參數值進行賦值，進行軸向加載，完成數據測定后繪制應力-應變曲線，根據測定數據計算宏觀力學參數。

為了確定細觀力學參數對巖土體宏觀力學參數與力學特性的影響，需要對細觀力學參數做敏感性分析，以此定性分析宏細觀參數間的關系。根據以往學者[9-11]的研究，主要包括黏結強度、摩擦系數、剛度比、顆粒形狀的影響。黏結強度決定了細觀層面上顆粒間膠結強弱度，從而影響宏觀的峰值強度和應力大小，若增大法向平行黏結強度與切向平行黏結強度，會引起彈性模量和峰值強度的增大。當摩擦系數增大時，彈性模量和峰值強度也在增長但漲幅較小。而顆粒的法向剛度與切向剛度之比將對泊松比產生影響，隨著剛度比的增大，泊松比也會增大。模擬實驗中顆粒形狀采用球形，不規則顆粒形狀雖然能夠更好地模擬現實中巖石的幾何形態，但考慮到該類塊體形態難以考慮到土石間凹凸不平引起的咬合作用，使得材料模擬結構與實際差異較大引起實驗結論不準確，因此后續模擬將不考慮這一因素。

一般獲取細觀力學參數的方法需要假定細觀力學參數，通過PFC3D軟件內三軸實驗獲取宏觀力學參數并與實際室內試驗數據進行對比來獲取符合參數。該方法較為耗時，因此本文采用機器學習方法中的BP神經網絡，對上述三軸試驗的數值模擬過程獲取數據進行訓練。具體流程是讓隨機組合的細觀參數，以及模擬相對應宏觀力學參數作為輸出樣本和輸入樣本，通過機器學習模型進行訓練，根據訓練樣本構建出宏觀參數與細觀參數之間的非線性映射關系，對細觀力學參數進行反演。將10組宏觀力學參數選做測試樣本來驗證模型，最終將實際的宏觀力學參數帶入模型中對細觀力學參數進行反演(表1)，分別獲取到最合適細觀標定參數。

應用反演精度公式[11]、平均殘差公式[12]、殘差均方差公式[13]來驗證BP神經網絡模型的性能，經過驗證其精度滿足建模要求。將表1獲得的細觀力學參數導入PFC3D軟件的三軸壓縮試驗中，獲得應力-應變曲線(圖7)。

表1 模型細觀力學參數Table 1 Mesomechanical parameters of model

由圖7可知，依照摩爾-庫倫強度理論計算滑坡黏聚力和內摩擦角[14-17]，計算同室內三軸壓縮實驗數據對比，由表2所得的對比值可知，其宏觀力學參數值非常相近，證明該方法可行，可以用于構建露天礦滑坡的顆粒流模型。

圖7 應力-應變曲線示意圖Fig.7 Schematic diagram of stress-strain curve

表2 模擬值的宏觀力學參數對比Table 2 Comparison of macro-mechanical parameters of simulated values

3 邊坡滑移過程模擬與分析

3.1 滑移運動過程模擬與分析

建立的滑坡三維模型(圖8)，滑移區域共生成11 002顆顆粒，顆粒高程差約為19 m。在對顆粒附加細觀參數后，施加重力即可進行滑坡運動特征的模擬?；履M過程中分別從整體區域滑移和區域特征點兩個方面分析，在滑坡區域分別選取9個特征點作為監測點(圖9)，在模擬中對其速度和位移進行監測，監測點的坐標分別為A(-8，13，10)、B(-3，14，10)、C(3，14，10)、D(-7.5，10，7)、E(-2.5，11，7)、F(2.5，11，7)、G(-5，5，3.5)、H(-1.5，5.5，3.5)、I(2，5.5，3.5)。

圖8 原始區域PFC滑坡模型Fig.8 Original regional PFC landslide model

圖9 監測點位置與編號Fig.9 Location and number of monitoring points

在將滑坡體的速度和位移清零后，完成顆粒與顆粒、顆粒與墻體間的細觀參數賦值，并刪除滑移坡面墻體，最后對滑坡模型施加重力，讓顆粒在重力作用下進行滑移，其整體運動過程如圖10所示。由圖10可知，滑坡從失穩到運動停止一共經歷20 s。在邊坡滑移前期中，1 s時坡面發生明顯形變，坡面頂端首先出現形變，坡表面顆粒速度達到了約10 m/s；到2 s時坡體迅速下滑并開始形成堆積，整體坡面繼續進行滑移，坡底端顆粒速度達到15 m/s，滑坡整體運動速度達到峰值；在3 s左右隨著堆積形成，滑移速度逐漸減慢，底端顆粒在堆積后繼續移動，底端局部速度為17.5 m/s；到達5 s時，頂端沉降速度開始減慢，隨著滑移、摩擦對能量損耗，滑坡運動速度顯著減慢，底端局部顆粒仍舊向前移動。隨著大部分顆粒的堆積，在10 s時滑坡體后端逐漸趨于穩定，前中端局部顆粒保持2.0 m/s速度繼續移動；在15 s時滑坡體前端整體開始穩定，直至20 s時大部分顆?；眷o止，最終形成完整的堆積體，只有小部分土體以較小速度繼續運動。

圖10 邊坡滑移過程Fig.10 Slope slip process

應用PFC程序自編FISH函數進行滑坡整體的位移與速度變化監測(圖11)。由圖11可知，滑坡平均速度在快速提升后又急速減緩，逐漸趨于平緩，在1.8 s時達到峰值約為5.6 m/s；平均位移不斷增加，但是增長速度逐漸減緩，最終位移距離約為13 m。

圖11 邊坡滑移平均速度和位移Fig.11 Average speed and displacement of slope slip

完成整體滑移運動分析后，分別對滑坡各區域速度和位移變化進行分析，監測點的顆粒運動速度變化見圖12。其中，圖12(a)中點C峰值速度最高，約為6.9 m/s，點A、點B的速度波動較為接近，峰值速度略低于整體平均速度，三個顆粒經過10 s都停止運動；圖12(b)中三個點峰值速度都與整體滑破平均速度接近，期間點F在3.5 s時有強烈波動，同頂端點運動停止時間相近；圖12(c)的底端點運動時間相比于中頂端的點，運動時間和整體滑破運動時間相同，經過20 s后點H、點I靜止，而點G還在緩慢移動。

圖12 滑坡各部分監測點速度變化曲線Fig.12 Velocity change curve of monitoring points in each part of the landslide

對比各點速度變化圖，對圖13中各部分監測點位移曲線變化進行分析。其中，圖13(a)中點B位移量最多為12.0 m；圖13(b)中點E位移量略高于點B，在10 s時達到12.2 m后靜止；圖13(c)中點H、點I經過20 s后停止移動，位移量分別為13.4 m和14.4 m，而點G經過22 s后位移量還在增加。

圖13 滑坡各部分監測點位移變化曲線Fig.13 Displacement change curve of monitoring points in each part of the landslide

通過對各部分顆粒速度位移變化分析可以看出，每個位置上監測點的速度變化中，靠近滑坡右側的點C、點F、點I的速度波動頻率最明顯，而點C、點F兩點對比相同位置上其他點位移量都最少，說明滑坡右側更為陡峭，而滑坡左側地勢較為平坦，所以底端點G位移量最多，并在滑坡整體基本靜止后仍有微量位移?；马敹撕椭卸私涍^10 s靜止，底端顆粒在形成堆積后繼續移動，監測點速度也在快速減緩后再次升高。

完成模擬后，對滑坡整體滑移峰值平均速度進行計算[15]，通過計算得到峰值為6.2 m/s，較高于模擬結果值。考慮到該理論忽略了巖體顆粒間的摩擦、碰撞的動能消耗，因此誤差可以近似忽略，可以看出滑坡的反演結果和實際滑坡堆積結果基本相似，證明了實驗的準確性。

3.2 滑移堆積過程分析

根據滑坡的滑移運動過程模擬分析，滑坡前邊坡最大寬度約為19 m，滑坡的滑移區域和堆積區域有部分重疊。為了更直觀展示滑坡的堆積特征，按照高程將這片區域的顆粒劃分為不同的組，對分組顆粒進行染色。圖14對比了滑坡前后的分組顆粒，能夠觀察到頂端顆粒、中端顆粒排序并未因滑坡發生而產生明顯變化，只有低端顆粒在運動時發生明顯混合。其中，頂端的兩個堆積區域和滑移區域發生重疊，左側部分顆粒被掩埋并隆起，中后端形成堆積平面。通過剖面圖也可看出堆積主要產生在后端，這與滑坡后實地調查堆積特征基本一致。

圖14 分組顆粒堆積分布Fig.14 Grouped particle packing distribution

通過滑坡運動軌跡和堆積特征分析，證明了離散元能夠有效地模擬礦山滑坡特征與運動趨勢，其數值模擬預測結能為滑坡的防治提供依據。

4 結 論

1) 本文將無人機調查同PFC3D三維數值建模相結合，通過無人機航攝構建的滑坡DEM與復現原始地形，以Geometry Import指令在PFC3D內構建出離散元三維數值模型，為后續滑坡的數值模擬提供了數據支撐。

2) 通過在三軸壓縮模擬實驗中加入BP神經網絡的方法，完成對宏細觀參數的反演與標定。經過驗證，將獲取的細觀參數帶入PFC3D內獲取數據與室內試驗對比，結果能夠滿足要求。

3) 滑坡整體歷時20 s，在1.8 s時平均速度達到峰值5.6 m/s。滑坡中頂端顆粒下落過程中呈現聚集狀，底端趨向于向左側運動，對比監測點速度位移曲線運動，滑坡體右端較為陡峭。堆積體主體主要為中頂端顆粒，滑坡與堆積體有部分重疊，與實際堆積特征基本相符，證明本研究的模擬方法可以應用于滑坡的運動特征研究。