復雜滑坡三維離散建模方法及其在茂縣滑坡中的應(yīng)用*

陳 骎 范 剛② 周家文②

（①四川大學水利水電學院，成都610065，中國）

（②四川大學水力學與山區(qū)河流開發(fā)保護國家重點實驗室，成都610065，中國）

0 引 言

數(shù)值模擬可以重現(xiàn)滑坡破壞和失效后的復雜運動過程，是研究滑坡災(zāi)害的重要手段之一（杜鵑等，2015；李萌等，2018；馬鵬輝等，2018；施家杰等，2019）。3DEC（Three-dimensional Distinct Element Code）是一款由美國Itasca公司開發(fā)的三維塊體離散元計算軟件，采用顯示差分法進行求解運動方程，能夠準確地模擬外力作用下塊體的運動過程（Hart et al.，1988）。相對于基于連續(xù)介質(zhì)力學的數(shù)值方法（如FEM、FLAC等）和基于顆粒離散元的數(shù)值方法（如PFC、EDEM等），3DEC中的塊體接觸模型能夠更加準確地反映塊體間的相互作用，因此，3DEC更適用于滑坡的運動過程模擬（顧金等，2016；Wu et al.，2017，2018；Shen et al.，2019）。

建立幾何模型是進行數(shù)值模擬的第一步，3DEC軟件對于復雜幾何模型的建模提供了以下兩種方法：（1）在Rhinoceros軟件中建立邊坡solid實體，導出為VRML格式文件，再由內(nèi)置的VTPT插件將其轉(zhuǎn)為3DEC命令流；（2）通過Gridle插件將AutoCAD創(chuàng)建的閉合幾何圖形轉(zhuǎn)換為3DEC塊體（Itasca Consulting Group，Inc.，2016）。3DEC提供的兩種建模方法難以滿足復雜滑坡數(shù)值模擬快速建模的要求。為了克服3DEC軟件在建模方面的不足，一些學者進行了一些探索。韓浩亮等（2012）利用Surfer軟件和FISH語言編程，建立了復雜地表及簡單地層的地質(zhì)體模型。趙玉玲等（2018）采用MATLAB建立了復雜礦區(qū)的地表模型。石崇等（2016）提出一種自下向上的建模方法，利用AutoCAD、ANSYS軟件及FORTRAN語言實現(xiàn)了具有復雜坡面形態(tài)邊坡的模型構(gòu)建。

滑坡數(shù)值分析模型通?？煞譃榛病⒒w兩個部分，但是，對于高速滑坡而言，運動過程中常常伴隨有強烈的鏟刮現(xiàn)象，建模時還應(yīng)考慮滑動路徑上的基底物質(zhì)。一方面，滑動路徑上的基底物質(zhì)將影響滑坡物源的運動軌跡及運動距離；另一方面，由于基底物質(zhì)的存在，滑坡災(zāi)害往往體現(xiàn)出放大效應(yīng)，滑坡堆積區(qū)的物質(zhì)量遠大于物源區(qū)的物質(zhì)量（周家文等，2019）。因此，對滑坡運動路徑上的基底物質(zhì)進行精確建模并準確模擬滑坡失穩(wěn)過程中的鏟刮效應(yīng)是模擬滑坡災(zāi)害的關(guān)鍵。趙玉玲等（2018）所提出的方法只適用于簡單的滑坡建模，難以實現(xiàn)滑坡體和松散基底物質(zhì)的建模。韓浩亮等（2012）和石崇等（2016）提出的方法難以對基底物質(zhì)進行建模，且要使用多個軟件耦合建模與編程，過程較復雜。本文在總結(jié)前人經(jīng)驗的基礎(chǔ)上，提出了基于MATLAB平臺的3DEC復雜滑坡模型的前處理方法，并以四川省茂縣新磨村滑坡為例詳細敘述了滑坡模型建模步驟，對類似的滑坡運動過程模擬有一定的參考借鑒意義。

1 復雜滑坡建模方法

1.1 塊體模型生成

本方法采用三棱柱塊體構(gòu)建滑床與滑體，利用Voronoi多面體模擬滑動路徑上的基底物質(zhì)。下邊對上述2類塊體模型生成方法分別介紹。

1.1.1 三棱柱體生成

Distmesh是一款開源的MATLAB網(wǎng)格劃分工具箱（Persson et al.，2004），可對任意形狀的多邊形進行三角網(wǎng)格剖分。通過其中的distmesh2D函數(shù)生成二維三角網(wǎng)格，輸出結(jié)果為包含所有節(jié)點坐標和各個三角單元節(jié)點序號的數(shù)組。

進一步以滑坡前、后的表面點云三維坐標為基礎(chǔ)，通過線性插值得到網(wǎng)格節(jié)點的高程坐標，將二維網(wǎng)格擴展到三維。對于構(gòu)成滑床的三棱柱塊體，下表面節(jié)點的高程指定為一常數(shù)，上表面節(jié)點的高程由滑坡后地表點云進行插值得到。對于構(gòu)成滑體的三棱柱塊體，下表面即為滑動面，其高程由滑坡后的地表點云進行插值得到，上表面節(jié)點的高程由滑坡后的地表點云進行插值得到。需要注意的是，由于插值引入誤差以及坐標數(shù)據(jù)誤差的影響，可能出現(xiàn)上表面點高程低于下表面對應(yīng)點高程的情況，因此在插值后需對高程進行檢查，若上表面點高程小于下表面點對應(yīng)高程，則需對高程進行修正，本文方法取上表面點高程等于下表面點高程加上0.1im。生成的三棱柱體示意圖（圖1）。

1.1.2 Voronoi多面體生成

圖1 插值得到的三棱柱體Fig.1 Triangular prism generated by interpolation

MPT（Multi-Parametric Toolbox 3.0）是 在MATLAB環(huán)境下開發(fā)的用于約束最優(yōu)控制器建模、控制、分析和部署的算法集合（Herceg et al.，2013）。它有一個強大的幾何庫，其將工具箱的應(yīng)用擴展至優(yōu)化控制之外，以解決計算幾何中出現(xiàn)的各種問題。由幾何庫中的mpt＿voronoi函數(shù)可生成一系列Voronoi多面體。函數(shù)通過在指定凸多面體區(qū)域內(nèi)部進行Voronoi劃分以達到生成多個Voronoi多面體的目的，區(qū)域邊界由凸多面體的頂點控制，生成多面體的尺寸由外邊界與輸入凸多面體中心坐標決定，生成的凸多面體數(shù)量與輸入多面體中心坐標數(shù)相同。mpt＿voronoi函數(shù)的輸出結(jié)果中包含有各多面體的頂點坐標。

1.2 塊體拓撲信息獲取

描述一個三維塊體結(jié)構(gòu)，必須明確其包含的點、線、面間的拓撲關(guān)系。以3DEC中的塊體生成命令poly face為例，多面體由帶有face關(guān)鍵字的平面構(gòu)成，其中頂點坐標（x1，y1，z1），（x2，y2，z2）等定義了面的邊界，并且頂點必須按順時針順序輸入（從塊外部觀察塊面）。為此，需首先獲取塊體的所有面及面上點的信息，再將各面上點的坐標進行順時針排序。對于三棱柱體，上、下表面的節(jié)點序在之前的網(wǎng)格劃分與插值后已經(jīng)存在，對上、下表面的節(jié)點進行兩兩組合，得到3個側(cè)面的節(jié)點序號；而由MPT工具箱生成的Voronoi多面體P，由不等式組（1）定義：

式中：x為多面體內(nèi)部及面上的點坐標；H為m×3矩陣；m為多面體面的數(shù)量，矩陣的每一行為對應(yīng)面的法向量，H×x＝K即包含m個面的平面方程。計算所有頂點到每個面的距離，若距離為0，則此點位于該面上，通過循環(huán)計算即可獲取每個面所包含點的坐標。使用MATLAB中的二維cell數(shù)組來存放2類塊體的拓撲信息，每個塊體的信息按行進行存放，每個面的信息按列存放，每個非空的cell單元為n行3列的數(shù)組，用來存放每個面上n個節(jié)點的三維坐標。

接著對各面上的點進行順時針排序，對以上2類塊體均采用同一算法進行頂點排序。上述2類塊體皆為凸多面體，均由凸多邊形表面構(gòu)成，故只需對表面局部坐標系下的頂點按極角大小排序即可。如圖2所示，以Voronoi多面體為例，在多面體表面內(nèi)建立局部坐標系，以凸多邊形表面內(nèi)一點O為局部坐標系原點（本方法取各頂點坐標平均值），極軸方向為原點O指向該cell單元中第一個點P1方向，并記第一點的極角為0°，圖2中點C可為多面體內(nèi)部任意一點（本方法C點坐標取多面體所有頂點坐標的平均值）。該cell單元中第n（n＞1）個點Pn的極角α按式（2）計算：

最后按極角由大到小對各點重新排序，即實現(xiàn)面上點的順時針排序。

圖2 面上點順時針排序示意圖Fig.2 Schematic diagram of clockwise sorting of points on the surface

1.3 3DEC命令輸出

將存放于cell數(shù)組中的塊體坐標信息寫成滿足3DEC語法要求的命令流，并輸出為txt格式文本。poly prism是3DEC中專用于生成棱柱形塊體的命令，由于只需要輸入棱柱上下兩個面的頂點坐標，相比poly face命令更加簡便，讀取速度更快，韓浩亮等（2012）和趙玉玲等（2018）提出的方法均按poly prism命令格式輸出。但在實際運用中發(fā)現(xiàn)，當棱柱上、下表面距離很近時，3DEC中無法生成棱柱形的滑體。為避免這一缺陷，三棱柱體與Voronoi多面體均按poly face命令要求的格式輸出。

1.4 復雜滑坡建模實現(xiàn)方法

綜上所述，本研究提出的復雜滑坡的建模實現(xiàn)方法主要包括以下幾個步驟：

（1）進行塊體生成。采用三棱柱體構(gòu)建滑床與滑體，通過Distmesh工具將滑坡區(qū)域劃分為二維三角網(wǎng)格并通過高程插值將二維三角面拓展為三棱柱體。采用Voronoi多面體構(gòu)建基底物質(zhì)，通過MPT工具箱生成Voronoi多面體。

（2）進行塊體拓撲信息的獲取。通過節(jié)點組合的方式獲取三棱柱面的信息，通過判斷點與面的間距獲取Voronoi多面體面的信息，隨后對2類多面體面上點進行順時針排序。

（3）接著按poly face命令格式將塊體的拓撲信息轉(zhuǎn)化輸出為3DEC命令流文本文件。

（4）最后在3DEC中調(diào)用塊體生成命令流文本文件，即實現(xiàn)滑坡模型的建立。

復雜滑坡的建模方法步驟流程圖（圖3）。

圖3 滑坡建模流程圖Fig.3 Landslide modeling flow diagram

2 建模函數(shù)

為實現(xiàn)3DEC復雜滑坡建模，依據(jù)上節(jié)提出的復雜滑坡建模方法編寫了相應(yīng)的MATLAB函數(shù)。本節(jié)對建模過程中涉及的Distmesh、MPT工具箱中的部分函數(shù)及自編函數(shù)進行介紹。

2.1 Distmesh中的函數(shù)

2.1.1 function d＝drectangle（p，x1，x2，y1，y2）

本函數(shù)功能是計算平面點到矩形邊界的最短距離。輸入p為n×2的數(shù)組，存有n個點的x，y坐標；函數(shù)輸入?yún)?shù)x1，x2，y1，y2決定了矩形區(qū)域的4個頂點坐標（x1，y1）（x1，y2）（x2，y1）（x2，y2）。函數(shù)輸出d為包含n個點到矩形邊界距離的列向量，若點在矩形內(nèi)部，則距離為負，反之則為正。

2.1.2 function d＝dpoly（p，pv）

本函數(shù)功能是計算平面點到多邊形邊界的最短距離。輸入p為n×2的數(shù)組，存有n個點的x，y坐標；輸入pv為m×2的數(shù)組，為多邊形的m個頂點坐標。函數(shù)輸出d為包含n個點到多邊形邊界距離的列向量，若點在多邊形內(nèi)部，則距離為負，反之則為正。

2.1.3 function d＝ddiff（d1，d2）

本函數(shù)功能是取d為d1與－d2中各行元素的較大值。函數(shù)輸入?yún)?shù)d1與d2為兩個行數(shù)相同的列向量。函數(shù)輸出為與d1、d2行數(shù)相同的列向量。

2.1.4 function［p，t］＝distmesh2id（fd，fh，h0，bbox，pfix，varargin）

本函數(shù)功能為在指定區(qū)域內(nèi)生成三角形網(wǎng)格。函數(shù)輸入?yún)?shù)fd為點到網(wǎng)格邊界的距離函數(shù)，用以指定生成網(wǎng)格的區(qū)域邊界；輸入?yún)?shù)fh為確定網(wǎng)格節(jié)點間距的函數(shù)，用以決定網(wǎng)格的邊長；h0為網(wǎng)格節(jié)點間距的初始值；bbox為2×2數(shù)組［x min，y min；x max，y max］，表示生成初始網(wǎng)格節(jié)點的矩形區(qū)域，（x min，y min）、（x max，y max）為矩形的兩個頂點；pfix為n×2數(shù)組，表示固定的節(jié)點坐標；varargin為fd、fh函數(shù)的額外參數(shù)。函數(shù)輸出p為N×2數(shù)組，包含N個節(jié)點的x，y坐標；t為M×3數(shù)組，每行表示M個單元的3個節(jié)點序號。

2.2 MPT中的函數(shù)

2.2.1 function P＝polytope（V）

本函數(shù)功能為生成一個凸多面體。函數(shù)輸入?yún)?shù)V為n×3數(shù)組，存有凸多面體n個頂點的x，y，z坐標。函數(shù)輸出P為polytope結(jié)構(gòu)體，成員P.H與P.K分別為m×3和m×1數(shù)組，m為多面體面的個數(shù)，P＝{x｜H×x≤K}定義了P的范圍；成員P.V為n×3數(shù)組，存有凸多面體n個頂點的空間坐標。

2.2.2 function Pn＝mpt＿voronoi（points，Options）

本函數(shù)功能為在指定的邊界內(nèi)生成Voronoi多面體。函數(shù)輸入?yún)?shù)points為n×3數(shù)組，作為生成n個Voronoi多面體的中心坐標；函數(shù)輸入?yún)?shù)Options為結(jié)構(gòu)體，成員Options.bound為polytope結(jié)構(gòu)，作為多面體生成的邊界，成員Options.plot為1則繪制生成的Voronoi圖，為0則不繪圖。函數(shù)輸出Pn為n×1數(shù)組，n為生成多面體的個數(shù)，數(shù)組中每個元素為一個polytope結(jié)構(gòu)體。

2.3 自編函數(shù)

2.3.1 function pu＝interpo（pb，po）

本函數(shù)功能是利用下表面網(wǎng)格節(jié)點的平面坐標插值得到上表面網(wǎng)格節(jié)點的高程。函數(shù)輸入?yún)?shù)pb為下表面網(wǎng)格的節(jié)點坐標值，為n×3數(shù)組，po為插值需要的上表面節(jié)點坐標，為m×3數(shù)組。函數(shù)輸出pu為插值得到的平面坐標與pb對應(yīng)上表面節(jié)點坐標，為n×3數(shù)組。n為網(wǎng)格節(jié)點數(shù)目，m為原地形空間坐標點數(shù)目。

2.3.2 function C＝triprism（t，pu，pl）

本函數(shù)功能是根據(jù)上、下表面的點序號與點坐標值生成包含5個面的cell數(shù)組。函數(shù)輸入?yún)?shù)t為M×3數(shù)組，每行表示各個單元的3個節(jié)點序號，輸入?yún)?shù)pu為n×3數(shù)組，包含上表面網(wǎng)格節(jié)點的三維坐標，輸入?yún)?shù)pb為pu對應(yīng)的下表面網(wǎng)格節(jié)點坐標數(shù)組。函數(shù)輸出C為M×5的cell數(shù)組，包含各三棱柱體的面與點信息。n為網(wǎng)格節(jié)點數(shù)目，M為三棱柱體數(shù)目。

2.3.3 function C＝voronoiface（Pn）

本函數(shù)功能為生成Voronoi多面體cell數(shù)組。函數(shù)輸入?yún)?shù)Pn為M×1數(shù)組，其數(shù)組單元為polytope結(jié)構(gòu)體。函數(shù)輸出C為M×N的cell數(shù)組，M為Voronoi多面體的個數(shù)，N為各個多面體所包含面的數(shù)量的最大值。

2.3.4 function C＝clockwise（C）

本函數(shù)功能為對數(shù)組C中每一單元（對應(yīng)多面體的面）中的點進行順時針排序。

2.3.5 function C＝transfer（C，k，t，r）

本函數(shù)功能為對多面體進行縮放、平移和旋轉(zhuǎn)。函數(shù)輸入?yún)?shù)C為表示多面體的cell數(shù)組，k為縮放比例，t為1×3平移向量，包含沿x，y，z方向的移動距離，r為1×3向量，包含繞x，y，z軸的旋轉(zhuǎn)角度。函數(shù)輸出C為多面體進行縮放、平移和旋轉(zhuǎn)后的cell數(shù)組。

2.3.6 function polyface（C，filename）

本函數(shù)功能為將多面體拓撲信息輸出為3DEC命令流。輸入C為示多面體的cell數(shù)組，函數(shù)輸入?yún)?shù)filename為輸出路徑與文件名。函數(shù)輸出結(jié)果為指定路徑與名稱的txt格式文件。

3 實例應(yīng)用

新磨村滑坡位于四川省茂縣疊溪鎮(zhèn)新磨村（坐標東經(jīng)103°39′46″，北緯32°4′47″），為一典型高位順層巖質(zhì)滑坡，滑坡于2017年6月24日失穩(wěn)，共計造成83人死亡或失蹤，滑體堵塞松坪溝，形成堰塞湖（許強等，2017）。如圖4所示，滑源區(qū)位于松坪溝左岸富貴山西坡山脊處，滑體頂部高程3460im，剪出口高程約為3100im，滑動面平均坡度約47°?；磪^(qū)長和寬約分別為200im和300im，平均厚度約46im，體積約為450×104im3。滑源區(qū)出露基巖為三疊系中統(tǒng)雜谷腦組（T2z）變質(zhì)砂巖夾板巖，巖層面產(chǎn)狀為N70°W／SW∠51°，構(gòu)成滑體的底部邊界?；磪^(qū)巖體內(nèi)發(fā)育有2組結(jié)構(gòu)面，產(chǎn)狀分別為N7°E／NW∠71°和N40°E／NW∠29°，前者構(gòu)成滑體的兩側(cè)邊界（Su et al.，2017）?；磪^(qū)以下存在明顯的流通鏟刮區(qū)，流通鏟刮區(qū)主要為崩坡積碎石土，由巖塊、礫石和1933年疊溪地震誘發(fā)的古滑坡堆積體等組成，分布高程為2650～3100im，長度和寬度約分別為700im和410im，平均鏟刮深度約30im。滑體在向下運動中強烈刮鏟坡面原有松散堆積物，滑體體積不斷增大，最終形成總體積約為1300×104m3的滑坡堆積體。堆積體位于2650im高程以下至2300 m高程區(qū)域，呈扇形堆積，順滑向長度約1600im，順河向最大寬度為1080im。通過比較滑坡前后地形高程變化，滑坡的平均堆積厚度大于10 m，最大堆積厚度約為31 m（圖5）。

本文以新磨村滑坡為例，詳細敘述復雜滑坡的建模方法。

（1）數(shù)據(jù)準備。準備滑坡發(fā)生前、后的地表點坐標數(shù)據(jù)（可從CAD文件等高線數(shù)據(jù)中提?。┥W(wǎng)格節(jié)點高程插值使用，保存為txt格式文件，每行存放一個點的三維坐標數(shù)據(jù)。選取滑源區(qū)、流通鏟刮區(qū)邊界上的點，作為滑源區(qū)、流通鏟刮區(qū)網(wǎng)格劃分的邊界。

圖4 新磨村滑坡滑后數(shù)字高程模型Fig.4 Digital elevation model of Xinmo landslide after sliding

圖5 新磨滑坡高程變化圖Fig.5 Elevation variation of Xinmo landslide

（2）三角形網(wǎng)格劃分。讀取滑源區(qū)、流通鏟刮區(qū)邊界點文件并命名為bound1和bound2，滑源區(qū)和流通鏟刮區(qū)邊界點的命名為bound3。bound1為滑源區(qū)的外邊界，bound2為流通鏟刮區(qū)的外邊界，bound3為滑源區(qū)與流通鏟刮區(qū)的共同邊界，整個模型的外邊界為一矩形，用ddiff函數(shù)計算內(nèi)、外區(qū)域的差集。采用distmesh2d函數(shù)進行3個區(qū)域的網(wǎng)格劃分，網(wǎng)格邊長統(tǒng)一為40im。可在MATLAB中繪制網(wǎng)格并進行網(wǎng)格質(zhì)量檢查，生成的網(wǎng)格如圖6所示。進一步設(shè)置模型的底部高程為2000im，將二維三角網(wǎng)格轉(zhuǎn)化為三維網(wǎng)格。MATLAB命令如下：

＞＞bound1＝load（′E：＼滑源區(qū)邊界.txt′）；

＞＞bound2＝load（′E：＼流通鏟刮區(qū)邊界.txt′）；

＞＞bound3＝load（′E：＼內(nèi)邊界.txt′）；

＞＞［pb1，t1］＝distmesh2id（＠dpoly，＠huniform，40，［270，1600；750，2900］，bound1，bound1）；

＞＞［pb2，t2］＝distmesh2id（＠dpoly，＠huniform，40，［270，1600；750，2900］，bound2，bound2）；

＞＞fd＝inline（′ddiff（drectangle（pb3，－40，1370，10，2950），dpoly（pb3，bound3））′，′pb3′，′bound3′）；

＞＞pfix＝［bound3；1370，2950；1370，10；－40，2950；－40，10］；

＞＞［pb3，t3］＝distmesh2id（fd，＠huniform，40，［-40，10；1370，2950］，pfix，bound3）；

＞＞len＝［size（pb1，1），size（pb2，1），size（pb3，1）］；

圖6 Distmesh生成的三角網(wǎng)格Fig.6 Triangle meshes generated with Distmesh

＞＞pb1＝［pb1,2000*ones（len（1），1）］；

＞＞pb2＝［pb2,2000*ones（len（2），1）］；

＞＞pb3＝［pb3,2000*ones（len（3），1）］；

（3）底部地形建模。對于滑源區(qū)和流通區(qū)，上表面網(wǎng)格高程用滑后點云插值得到，對于外部區(qū)域，上表面網(wǎng)格用滑前點云插值得到。面上點順時針排序，輸出為3DEC命令文本。MATLAB命令如下：

＞＞pre＿surface＝load（′E：＼滑前.txt′）；

＞＞post＿surface＝load（′E：＼滑后.txt′）；

＞＞pu1＝interpo（pb1，post＿surface）；

＞＞pu2＝interpo（pb2，post＿surface）；

＞＞pu3＝interpo（pb3，pre＿surface）；

＞＞C1＝triprism（t1，pu1，pb1）；C2＝triprism（t2，pu2，pb2）；

＞＞C3＝triprism（t3，pu3，pb3）；

＞＞base＝［C1；C2；C3］；base＝clockwise（base）；

＞＞polyface（base，′E：＼base＿3DEC.txt′）；

（4）滑體建模?；磪^(qū)下表面網(wǎng)格采用上一步插值得到的t1與pu1，上表面網(wǎng)格利用滑前點云插值得到。對C4重新排序，設(shè)置輸出路徑與文件名，輸出為3DEC命令文本。MATLAB命令如下：

＞＞t4＝t1；pb4＝pu1；

＞＞pu4＝interpo（pb4，pre＿surface）；

＞＞C4＝triprism（t4，pu4，pb4）；

＞＞C4＝clockwise（C4）；

＞＞polyface（C4，′E：＼source＿3DEC.txt′）；

（5）流通鏟刮區(qū)堆積體建模。在長、寬、高分別為500im、300im、50im的長方體區(qū)域內(nèi)生成7500個Voronoi多面體。接著對所有多面體的面及面上的點進行計算、順時針排序。為產(chǎn)生堆積體松散堆積狀態(tài)，使塊體在重力作用下進行堆積，在輸出為3DEC命令文本前對塊體進行平移與旋轉(zhuǎn)。MATLAB命令如下：

＞＞P＝polytope（［0,0i0；300,0i0；300,500,0；0,500,0；0,0i50；300,0i50；300,500,50；0,500,50］）；

＞＞Options.plot＝1；Options.pbound＝P；

＞＞points＝［300*rand（1，7500）；500*rand（1，7500）；50*rand（1，7500）］；

＞＞Pn＝mpt＿voronoi（points′，Options）；

＞＞polytopes＝voronoi＿face（Pn）；

＞＞polytopes＝clockwise（polytopes）；

＞＞polytopes＝transfer（polytopes，1，［450，1960，2800］，［-30，0，0］）；

＞＞polyface（polytopes，′E：＼polytope3DEC.txt′）；

（6）在3DEC中使用call命令讀取MATLAB輸出的文本文件。塊體生成后，根據(jù)實際結(jié)構(gòu)面產(chǎn)狀重新劃分滑源區(qū)塊體，使用join on、jset命令生成結(jié)構(gòu)面。需要注意的是，若3DEC工程坐標系與大地坐標系不同，使用jset命令時需將結(jié)構(gòu)面產(chǎn)狀參數(shù)進行相應(yīng)變換。另外，考慮到計算效率，結(jié)構(gòu)面間距取20im。設(shè)置邊界條件、塊體與接觸本構(gòu)參數(shù)及重力作用方向，開始力學計算。在3DEC中，滑床作為固定邊界，塊體采用剛體本構(gòu)，接觸本構(gòu)為莫爾－庫侖滑動模型（Itasca Consulting Group，Inc.，2016），初始平衡階段計算參數(shù)取值見表1。

（7）達到平衡后使用geometry import命令導入dxf或stl格式的滑坡前地形，對流通鏟刮區(qū)的塊體進行分組，并刪除位于滑坡前地形高程以上的塊體，即完成流通鏟刮區(qū)堆積體的建模，建立的滑坡模型如圖7所示。受原始地形數(shù)據(jù)精度影響，建立的三維滑坡模型與實際原型存在一定差距。模型生成后，在3DEC中改變重力方向及接觸參數(shù)，開始滑坡運動過程模擬?；逻\動過程模擬階段參數(shù)取值見表1。

滑體對基底物質(zhì)的鏟刮包括滑體前部對基底物質(zhì)的鏟作用與滑體對底部基底物質(zhì)刮作用兩種形式（Barbolini et al.，2005；陸鵬源等，2016；Sovilla et al.，2006）?；w前部先與老滑坡堆積物接觸并將其鏟起（圖8a），強烈的碰撞導致塊體飛濺，隨著運動繼續(xù)發(fā)展，滑體翻越至老滑坡堆積物上部以刮的形式將其裹挾運移（圖8b），類似現(xiàn)象同樣在物理實驗中觀察到（陸鵬源等，2016）。圖9為模擬滑坡堆積結(jié)果，模擬堆積形態(tài)呈扇形，與實際堆積形狀接近。模擬出的堆積范圍順滑向水平長度約2450im，順河向水平寬度約1240im，與滑坡實際堆積范圍（順滑向長2500im，順河向?qū)?080im）較為吻合?；w在下滑過程中受流通鏟刮區(qū)西高東低地形影響，主滑方向略向東偏移，堆積厚度分布整體呈現(xiàn)出西低東高態(tài)勢（圖5）?？梢钥闯?，模擬堆積體厚度分布（圖10）與實際堆積趨勢一致（圖2），表明采用三棱柱體能夠有效地模擬出滑床的地形特征。圖11為堆積區(qū)1－1′剖面堆積情況，模擬堆積厚度在中部與實際情況吻合較好，而在東、西兩側(cè)大于實際堆積厚度。該差異是由于數(shù)值模擬中采用了較小的接觸摩擦角，只有極少的滑體停積滑源區(qū)與流通鏟刮區(qū)，導致模擬的整體堆積厚度比實際堆積厚度略大?？傮w來看，模擬結(jié)果與實際具有較好的一致性。

表1 數(shù)值模擬參數(shù)取值Table 1 Parameters for the numerical analysis

圖7 3DEC滑坡模型Fig.7 Landslide model in 3DEC

4 討論與結(jié)論

（1）本文基于MATLAB編寫了復雜滑坡建模相關(guān)程序，借助滑坡前、后的地表高程數(shù)據(jù)，采用三棱柱體及Voronoi多面體對滑體、滑床和松散堆積物進行構(gòu)建。此方法無須借助第三方軟件進行網(wǎng)格劃分，同時在MATLAB中完成塊體生成及3DEC命令流輸出，簡化了建模流程。

圖8 滑坡運動過程中的鏟刮現(xiàn)象Fig.8 Entrainment process of the landslide

圖9 模擬堆積區(qū)域Fig.9 Simulated depositional area

圖10 模擬堆積厚度Fig.10 Simulated deposit thickness

圖11 “1－1′”截面模擬與實際堆積厚度Fig.11 Simulated and actual deposit thickness along“1－1′”cross section

（2）本方法是基于分區(qū)建模思路對滑坡進行建模，更適用于滑動面位置為已知的邊坡，對滑動面未知的邊坡而言，可考慮利用其他程序得到最危險滑面位置再進行模型的分區(qū)構(gòu)建。另外，建立含結(jié)構(gòu)面的滑體模型還需在3DEC中對塊體進行重新分割，如何在MATLAB中實現(xiàn)這一功能還需進一步研究。

（3）本研究以新磨滑坡為例，詳細敘述了該方法的建模步驟，并對滑坡運動過程進行反演。模擬很好地展現(xiàn)了滑坡過程中的滑體對沿程堆積物的鏟、刮現(xiàn)象。同時，模擬得到的滑坡堆積形態(tài)，滑坡運動距離和滑坡堆積厚度均與實際情況表現(xiàn)出較好的一致性，驗證了該建模方法的可行性。