陳 骎 范 剛② 周家文②

（①四川大學(xué)水利水電學(xué)院，成都610065，中國(guó)）

（②四川大學(xué)水力學(xué)與山區(qū)河流開(kāi)發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，成都610065，中國(guó)）

0 引 言

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

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

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

1 復(fù)雜滑坡建模方法

1.1 塊體模型生成

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

1.1.1 三棱柱體生成

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

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

1.1.2 Voronoi多面體生成

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

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

1.2 塊體拓?fù)湫畔@取

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

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

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

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

圖2 面上點(diǎn)順時(shí)針排序示意圖Fig.2 Schematic diagram of clockwise sorting of points on the surface

1.3 3DEC命令輸出

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

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

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

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

（2）進(jìn)行塊體拓?fù)湫畔⒌墨@取。通過(guò)節(jié)點(diǎn)組合的方式獲取三棱柱面的信息，通過(guò)判斷點(diǎn)與面的間距獲取Voronoi多面體面的信息，隨后對(duì)2類多面體面上點(diǎn)進(jìn)行順時(shí)針排序。

（3）接著按poly face命令格式將塊體的拓?fù)湫畔⑥D(zhuǎn)化輸出為3DEC命令流文本文件。

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

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

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

2 建模函數(shù)

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

2.1 Distmesh中的函數(shù)

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

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

2.1.2 function d＝dpoly（p，pv）

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

2.1.3 function d＝ddiff（d1，d2）

本函數(shù)功能是取d為d1與－d2中各行元素的較大值。函數(shù)輸入?yún)?shù)d1與d2為兩個(gè)行數(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為點(diǎn)到網(wǎng)格邊界的距離函數(shù)，用以指定生成網(wǎng)格的區(qū)域邊界；輸入?yún)?shù)fh為確定網(wǎng)格節(jié)點(diǎn)間距的函數(shù)，用以決定網(wǎng)格的邊長(zhǎng)；h0為網(wǎng)格節(jié)點(diǎn)間距的初始值；bbox為2×2數(shù)組［x min，y min；x max，y max］，表示生成初始網(wǎng)格節(jié)點(diǎn)的矩形區(qū)域，（x min，y min）、（x max，y max）為矩形的兩個(gè)頂點(diǎn)；pfix為n×2數(shù)組，表示固定的節(jié)點(diǎn)坐標(biāo)；varargin為fd、fh函數(shù)的額外參數(shù)。函數(shù)輸出p為N×2數(shù)組，包含N個(gè)節(jié)點(diǎn)的x，y坐標(biāo)；t為M×3數(shù)組，每行表示M個(gè)單元的3個(gè)節(jié)點(diǎn)序號(hào)。

2.2 MPT中的函數(shù)

2.2.1 function P＝polytope（V）

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

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

本函數(shù)功能為在指定的邊界內(nèi)生成Voronoi多面體。函數(shù)輸入?yún)?shù)points為n×3數(shù)組，作為生成n個(gè)Voronoi多面體的中心坐標(biāo)；函數(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為生成多面體的個(gè)數(shù)，數(shù)組中每個(gè)元素為一個(gè)polytope結(jié)構(gòu)體。

2.3 自編函數(shù)

2.3.1 function pu＝interpo（pb，po）

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

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

本函數(shù)功能是根據(jù)上、下表面的點(diǎn)序號(hào)與點(diǎn)坐標(biāo)值生成包含5個(gè)面的cell數(shù)組。函數(shù)輸入?yún)?shù)t為M×3數(shù)組，每行表示各個(gè)單元的3個(gè)節(jié)點(diǎn)序號(hào)，輸入?yún)?shù)pu為n×3數(shù)組，包含上表面網(wǎng)格節(jié)點(diǎn)的三維坐標(biāo)，輸入?yún)?shù)pb為pu對(duì)應(yīng)的下表面網(wǎng)格節(jié)點(diǎn)坐標(biāo)數(shù)組。函數(shù)輸出C為M×5的cell數(shù)組，包含各三棱柱體的面與點(diǎn)信息。n為網(wǎng)格節(jié)點(diǎn)數(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多面體的個(gè)數(shù)，N為各個(gè)多面體所包含面的數(shù)量的最大值。

2.3.4 function C＝clockwise（C）

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

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

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

2.3.6 function polyface（C，filename）

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

3 實(shí)例應(yīng)用

新磨村滑坡位于四川省茂縣疊溪鎮(zhèn)新磨村（坐標(biāo)東經(jīng)103°39′46″，北緯32°4′47″），為一典型高位順層巖質(zhì)滑坡，滑坡于2017年6月24日失穩(wěn)，共計(jì)造成83人死亡或失蹤，滑體堵塞松坪溝，形成堰塞湖（許強(qiáng)等，2017）。如圖4所示，滑源區(qū)位于松坪溝左岸富貴山西坡山脊處，滑體頂部高程3460im，剪出口高程約為3100im，滑動(dòng)面平均坡度約47°?；磪^(qū)長(zhǎng)和寬約分別為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，長(zhǎng)度和寬度約分別為700im和410im，平均鏟刮深度約30im?；w在向下運(yùn)動(dòng)中強(qiáng)烈刮鏟坡面原有松散堆積物，滑體體積不斷增大，最終形成總體積約為1300×104m3的滑坡堆積體。堆積體位于2650im高程以下至2300 m高程區(qū)域，呈扇形堆積，順滑向長(zhǎng)度約1600im，順河向最大寬度為1080im。通過(guò)比較滑坡前后地形高程變化，滑坡的平均堆積厚度大于10 m，最大堆積厚度約為31 m（圖5）。

本文以新磨村滑坡為例，詳細(xì)敘述復(fù)雜滑坡的建模方法。

（1）數(shù)據(jù)準(zhǔn)備。準(zhǔn)備滑坡發(fā)生前、后的地表點(diǎn)坐標(biāo)數(shù)據(jù)（可從CAD文件等高線數(shù)據(jù)中提?。┥W(wǎng)格節(jié)點(diǎn)高程插值使用，保存為txt格式文件，每行存放一個(gè)點(diǎn)的三維坐標(biāo)數(shù)據(jù)。選取滑源區(qū)、流通鏟刮區(qū)邊界上的點(diǎn)，作為滑源區(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ū)邊界點(diǎn)文件并命名為bound1和bound2，滑源區(qū)和流通鏟刮區(qū)邊界點(diǎn)的命名為bound3。bound1為滑源區(qū)的外邊界，bound2為流通鏟刮區(qū)的外邊界，bound3為滑源區(qū)與流通鏟刮區(qū)的共同邊界，整個(gè)模型的外邊界為一矩形，用ddiff函數(shù)計(jì)算內(nèi)、外區(qū)域的差集。采用distmesh2d函數(shù)進(jìn)行3個(gè)區(qū)域的網(wǎng)格劃分，網(wǎng)格邊長(zhǎng)統(tǒng)一為40im?？稍贛ATLAB中繪制網(wǎng)格并進(jìn)行網(wǎng)格質(zhì)量檢查，生成的網(wǎng)格如圖6所示。進(jìn)一步設(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）底部地形建模。對(duì)于滑源區(qū)和流通區(qū)，上表面網(wǎng)格高程用滑后點(diǎn)云插值得到，對(duì)于外部區(qū)域，上表面網(wǎng)格用滑前點(diǎn)云插值得到。面上點(diǎn)順時(shí)針排序，輸出為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)格利用滑前點(diǎn)云插值得到。對(duì)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ū)堆積體建模。在長(zhǎng)、寬、高分別為500im、300im、50im的長(zhǎng)方體區(qū)域內(nèi)生成7500個(gè)Voronoi多面體。接著對(duì)所有多面體的面及面上的點(diǎn)進(jìn)行計(jì)算、順時(shí)針排序。為產(chǎn)生堆積體松散堆積狀態(tài)，使塊體在重力作用下進(jìn)行堆積，在輸出為3DEC命令文本前對(duì)塊體進(jìn)行平移與旋轉(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ù)實(shí)際結(jié)構(gòu)面產(chǎn)狀重新劃分滑源區(qū)塊體，使用join on、jset命令生成結(jié)構(gòu)面。需要注意的是，若3DEC工程坐標(biāo)系與大地坐標(biāo)系不同，使用jset命令時(shí)需將結(jié)構(gòu)面產(chǎn)狀參數(shù)進(jìn)行相應(yīng)變換。另外，考慮到計(jì)算效率，結(jié)構(gòu)面間距取20im。設(shè)置邊界條件、塊體與接觸本構(gòu)參數(shù)及重力作用方向，開(kāi)始力學(xué)計(jì)算。在3DEC中，滑床作為固定邊界，塊體采用剛體本構(gòu)，接觸本構(gòu)為莫爾－庫(kù)侖滑動(dòng)模型（Itasca Consulting Group，Inc.，2016），初始平衡階段計(jì)算參數(shù)取值見(jiàn)表1。

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

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

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

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

4 討論與結(jié)論

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

圖8 滑坡運(yùn)動(dòng)過(guò)程中的鏟刮現(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′”截面模擬與實(shí)際堆積厚度Fig.11 Simulated and actual deposit thickness along“1－1′”cross section

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

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