孔祥曌，李 濱，賀 凱，羅 浩，常文斌，邢愛國

（1.上海交通大學海洋工程國家重點實驗室，上海 200240；2.中國地質科學院地質力學研究所，北京 100081；3.中國地質調查局，北京 100037）

0 引言

柱狀巖體崩塌在我國西南地區(qū)頻繁發(fā)生，如2004年8月12日重慶甑子巖崩塌、2013年2月18日貴州凱里市老山新村崩塌等，其具有分布范圍廣、破壞能力強、影響范圍大等特點。在巖體崩塌過程中，危巖體不斷破碎解體，形成由不同粒徑巖塊構成的碎屑物質集合體，即崩塌碎屑流，其致災范圍常超過數(shù)百米。因此，對此類崩塌動力特征與破碎規(guī)律的研究十分重要。

當前國內外學者對塔柱狀巖體崩塌的動力特征與破碎規(guī)律研究較少，更多的是將塔柱狀巖體簡化為柱體崩塌研究其運動與堆積規(guī)律。Crosta 等[1]基于PFC2D研究了不同高寬比及連接準則的柱體崩塌，結合試驗資料驗證了離散元模擬柱體崩塌動力過程的可行性。Utili 等[2]基于離散元研究了不同高寬比和摩擦條件下顆粒柱體表面隨時間的演變過程，驗證了離散元模擬柱體崩塌運動過程的有效性。Zhou 等[3]基于離散元研究了顆粒柱的不同幾何尺寸與顆粒摩擦系數(shù)的對混合粒徑干顆粒流運動及堆積特征的影響。Huang 等[4]將大型柱狀巖體的崩塌概化為顆粒柱的崩塌，通過室內試驗研究了柱體中不同粒徑組分的運動特征。在實際崩塌中，節(jié)理、裂隙的存在對巖體的崩塌有著至關重要的作用，而地形又對堆積區(qū)的形態(tài)、顆粒破碎程度有較大影響，因此對真實崩塌的數(shù)值模擬研究尤為重要。

一些學者對甑子巖的演化機制、力學性能與失穩(wěn)模式進行研究。陳智強與李渝生[5]對甑子巖的形成演化機制與防治措施進行了研究。賀凱等[6-8]基于現(xiàn)場調查、室內試驗及數(shù)值計算等手段對甑子巖W12 危巖體崩塌進行了詳細研究，提出了塔柱狀巖體壓裂潰屈崩塌破壞模式，并從損傷力學角度揭示了塔柱狀底部巖體的強度劣化機制。馮振等[9]基于UDEC 對甑子巖初始變形破壞進行了研究，發(fā)現(xiàn)甑子巖崩塌的觸發(fā)條件為巖體軟弱基座的破壞。孫敬輝等[10]基于Rockfall 軟件對甑子巖不同尺寸崩塌落石的運動軌跡、速率、等進行研究，將崩塌落石區(qū)進行了風險評估與危險性分區(qū)。但當前尚未有學者對甑子巖崩塌破碎全過程進行模擬研究。

基于MatDEM 離散元軟件對甑子巖W12 危巖體崩塌動力特征與破碎規(guī)律進行了研究，依據(jù)實際節(jié)理分布建立了崩塌模型，實現(xiàn)了對崩塌運動全過程的模擬，通過與現(xiàn)場影像對比驗證了模型的有效性，在此基礎上，通過對MatDEM 二次開發(fā)，統(tǒng)計分析了運動過程中巖塊粒徑分布演化規(guī)律，并引入分形維數(shù)與雙參數(shù)Weibull 分布模型對顆粒破碎程度進行了研究。

1 崩塌特征與地質氣候條件

1.1 甑子巖崩塌特征

重慶市金佛山景區(qū)甑子巖危巖帶的W12 危巖體于2004年8月12日12 時53 分發(fā)生了大規(guī)模山體崩塌（文中簡稱甑子巖崩塌）。甑子巖崩塌體崩塌體高度為250 m，寬度為50 m，體積約為50×104m3，運動距離距崖腳約600 m，崩塌前后最大高差超過500 m （圖1、圖2）。崩塌體結構破碎，節(jié)理裂隙十分發(fā)育，在墜落撞擊基底層時與空氣相互作用劇烈，形成顯著的超前空氣沖擊效應，激起浮塵高度近150 m。由于及時預警，此次崩塌未造成人員傷亡[11]。

圖1 甑子巖地貌特征與崩塌源區(qū)Fig.1 Image of the landform and source area of the Zengziyan rockfall

圖2 甑子巖崩塌衛(wèi)星影像Fig.2 Aerial image of the Zengziyan rockfall

1.2 地質氣候條件

甑子巖崩塌源區(qū)巖體所在斜坡結構為上部堅硬，中下部軟弱的“上硬下軟”型巖體結構（圖3）。危巖帶由兩級陡崖組成，一級陡崖由棲霞組和茅口組一段石灰?guī)r組成；二級陡崖由茅口組三、四、五段石灰?guī)r組成；危巖帶中均有巖溶發(fā)育，二級陡崖多發(fā)育有溶洞，一級陡崖底部石灰?guī)r與炭質頁巖接觸面上則以巖溶大泉為主要特征，局部仍見有溶洞。

由圖3可知，崩塌體位于甑子巖二級陡崖，呈三棱柱狀。其高約250 m，寬約50 m。崩塌源區(qū)巖層產(chǎn)狀為280°～300°∠4°～5°，呈緩內傾巖體結構，崩塌體主要受三組裂隙切割控制（圖4），① 310°∠89°（J1），延伸長約200 m，貫通至頂，張開寬度2 m 左右，充填碎塊石土，下部張開20～30 cm，無充填，裂面見鈣泥質物，構成危巖體北側邊界；② 25°∠89°（J2），延伸長度200 m，張開約10 cm，無充填，局部充填碎石，裂面見鈣質物，貫通至頂，構成南側邊界；③ 組裂隙（J3）為層間裂隙在危巖體中部、下部形成凹巖腔，高0.50 m。J1、J2兩組節(jié)理將山體分割成三棱柱形危巖體，危巖體西南面臨空，加之危巖體底部基座為軟弱頁巖，導致巖體底部在上覆巖體的自重作用下，容易發(fā)生壓縮流變及剪切流變。

圖3 甑子巖崩塌斷面圖（沿圖2中A-B 斷面）Fig.3 Geological profile of the Zengziyan rockfall along line A-B in fig.2

圖4 甑子巖崩塌節(jié)理分布Fig.4 Stereographic projection of joints of the Zengziyan rockfall

金佛山地區(qū)屬亞熱帶溫濕氣候，具有“氣候溫和、雨量充沛、綿雨久、濕度大”等特征。危巖體崩塌發(fā)生于該地區(qū)雨季（5—9月），在節(jié)理裂隙、溶蝕管道等通道的優(yōu)勢入滲作用下，降水會顯著降低硬巖底部軟弱層強度，巖溶作用進一步加大了這種趨勢，引起差異沉降或局部崩塌，反過來又導致巖體節(jié)理裂隙擴大，形成惡性循環(huán)，最終使巖體穩(wěn)定性降低，發(fā)生崩塌。

2 離散元數(shù)值模擬

2.1 MatDEM 模型建立

甑子巖崩塌MatDEM 模型如圖5所示。模型寬1 200 m，高800 m，其中源區(qū)寬38～42 m，高約250 m，高寬比a為5.95～6.58。模型主要由崩塌源區(qū)與剛性基底層組成，崩塌源區(qū)由活動顆粒膠結而成，可遵循牛頓第二定律計算產(chǎn)生速率及位移，剛性基底層由邊界剛性單元組成，僅參與受力計算，并不產(chǎn)生速率與位移。

圖5 甑子巖崩塌MatDEM 模型Fig.5 MatDEM model of the Zengziyan rockfall

根據(jù)巖體中存在著的三組節(jié)理面將模型源區(qū)分為三部分，源區(qū)下部高度為100 m 的區(qū)域由強度較低的膠結單元組成，源區(qū)上部高度為150 m 的區(qū)域由強度較高的膠結單元組成。因J1 與J2 兩組節(jié)理在二維平面模型中表現(xiàn)一致，因此節(jié)理僅分為兩組。橫向節(jié)理J3 與水平方向夾角為0°～5°，間距約20 m；豎向節(jié)理在橫向上位于源區(qū)中部，豎直延伸，上下各留約50 m。

模型以破壞源區(qū)下部20%高度內的單元連接作為模型的觸發(fā)條件，模擬基底壓碎的情況。

2.2 MatDEM 模型參數(shù)確定

賀凱等[6-8]曾對甑子巖地區(qū)采集的巖樣開展了巖體物理力學性質試驗。甑子巖崩塌源區(qū)中上部主要為二疊系下統(tǒng)茅口組四、五段，因二者強度相近，故采用相同強度參數(shù)。源區(qū)下部主要為二疊系下統(tǒng)茅口組三段，受巖溶發(fā)育、人類工程活動及氣候因素影響，巖體強度較低。因源區(qū)巖體在失穩(wěn)前已經(jīng)歷嚴重的風化或擾動，其強度較完整巖體有較大折減，所以文章所采用的材料參數(shù)在賀凱所得試驗結果上有所折減。

離散元法通過堆積具有特定力學參數(shù)的顆粒以更好地實現(xiàn)對自然巖土體散粒特性的模擬。其整體力學性質與其單元間的接觸關系、力學參數(shù)和堆積方式等有關[12]。因此采用MatDEM 內置的材料參數(shù)標定程序BoxMatTraining 進行離散元宏微觀參數(shù)轉換，其原理為對基于隨機堆積而成的試樣做模擬三軸試驗，并根據(jù)試驗結果不斷調整輸入?yún)?shù)，直至試驗結果與真實情況基本一致，以完成材料參數(shù)的訓練。材料力學參數(shù)如表1所示。

表1 模型材料力學參數(shù)Table 1 Mechanical parameters of the Zengziyan rockfall model

2.3 MatDEM 二次開發(fā)

基于MatDEM 計算數(shù)據(jù)進行二次開發(fā)以研究巖體破碎過程，文章參考了奚悅[13]的研究作為巖塊的確定準則。在巖塊的運動過程中，由于碰撞等因素，部分顆粒單元間的膠結會遭到破壞，另一部分顆粒單元間仍保持膠結狀態(tài)。文章認為存在連接的兩個單元屬于同一個巖塊?？紤]天然巖體中節(jié)理、裂隙分布的多樣性，由顆粒膠結組成的條狀結構、多個顆粒膠結形成閉環(huán)與在閉環(huán)上向外延伸單元均被判定為巖塊。

許多學者使用等效粒徑的概念研究巖體破碎的粒徑分布，經(jīng)過對比試算，選用等效粒徑計算公式如下：

式中：d——巖塊的等效粒徑；

Vf——巖塊體積；

V0——模型總體積。

二次開發(fā)原理：在MatDEM 每一計算時間步中均記錄有單元接觸矩陣d.mo.nBall 與單元連接分布矩陣d.mo.bFilter，單元接觸包括壓縮、拉伸、剪切三種接觸模式，壓縮接觸的單元間并不一定存在連接，因此單元接觸矩陣并不能表示與一個單元膠結的所有單元，將兩矩陣結合起來可獲得與一個單元接觸膠結的所有單元，將單元連接矩陣中所有具有相同元素的行歸納在一起，去除重復元素并與單元接觸矩陣相對照，即可得到巖塊分布矩陣。

2.4 分維特征與Weibull 分布模型

巖體崩塌過程中巖塊將不斷發(fā)生摩擦、撞擊，導致顆粒破碎，因此引入分形幾何理論中的分形維數(shù)D描述巖塊破碎后的粒度屬性，當粒徑分布越分散，巖塊粒徑越小，則D值越大。采用Gates-Gaudin-Schuhmann（GGS）分布計算分形維數(shù)[14]。

式中：M(d＜di)——等效粒徑小于di巖塊的累積質量；

MT——巖塊總質量；

dmax——巖塊的最大等效粒徑。

對式（2）等式兩邊取自然對數(shù)，則可變?yōu)槭剑?）：

式中：n——因變量M(d＜di)/MT對自變量l n(di/dmax)直線的斜率。

Turcotte[15]推導出在GGS 分布中分形維數(shù)D與式（3）中的n關系如下：

引入雙參數(shù)Weibull 分布模型對崩塌前后巖塊的粒徑分布進行描述，公式為：

式中：dc——尺度參數(shù)，可以用于衡量粒徑分布中細粒巖塊占比，dc值越小，模型中的細粒巖塊占比越大；

β——形狀參數(shù)，可用于衡量粒徑分布的廣度，β值越大，粒徑分布越窄。

3 模擬結果分析

3.1 崩塌速率演化與有效性驗證

甑子巖崩塌過程速率演化如圖6所示。整個崩塌過程持續(xù)26.8 s。當運動開始t=1.6 s 時，崩塌巖體頂部有明顯位移差產(chǎn)生，巖體中下部有水平位移趨勢，受地形作用，巖體底部形成低速三角區(qū)域，其單元速率基本為0，巖體其他部位單元速率基本相同，為15 m/s。當t=3.2 s 時，巖體頂部位移差繼續(xù)增大，最大位移達50 m，受單元速率與地形影響，中下部崩塌巖體水平位移明顯，部分單元已滑出一級陡崖，崩塌巖體平均速率增加至20 m/s。t=4.8 s 時，上部巖體有明顯沿節(jié)理裂隙面破裂跡象，破裂后的兩部分運動速率有明顯差異，其平均速率分別為26.8 m/s 與29.2 m/s，崩塌體底部已開始凌空飛行，最大速率為35 m/s，在低速三角區(qū)域與下落運動部分間有明顯速率過渡區(qū)出現(xiàn)。

圖6 甑子巖崩塌速度演化Fig.6 Velocity evolution of the Zengziyan rockfall

t=6.4 s 時，崩塌巖體已沿節(jié)理裂隙面明顯分離，凌空飛行速率已超過40 m/s，崩塌體前緣少量顆粒已開始撞擊底面剛性模型。t=8.0 s 時，斜坡地面已形成堆積，凌空飛行的巖體沿內部橫向節(jié)理面有解體跡象，豎向節(jié)理右側巖體墜落速率明顯增加，已超過35 m/s。t=9.6 s時，崩塌體處低速三角區(qū)域外均以進入凌空飛行狀態(tài)，豎向節(jié)理左側巖體速率為35.6 m/s。當t=12.8 s 時，甑子巖上部巖體已基本結束凌空飛行墜落至斜坡地面，此時破碎巖體仍具有較高的速率，破碎巖體將繼續(xù)發(fā)生摩擦、移動、撞擊、破碎解體等過程，直至動能全部消散。t=26.8 s 時，崩塌巖體各部分動能基本已消散完成，堆積區(qū)已經(jīng)形成，通過測量可知堆積區(qū)長度為450 m，最大堆積厚度為40 m。

通過與甑子巖崩塌影像進行對比，可證實MatDEM崩塌模擬的有效性。圖7為甑子巖（W12）崩塌過程影像。當T=6 s 時，源區(qū)中上部巖體間有碎屑濺射，說明此時按已經(jīng)發(fā)生破碎，與MatDEM 模擬中t=6.4 s 時情況一致。當T=9.3 s 時，崩塌體整體已開始離開一級陡崖進入凌空墜落狀態(tài)，與MatDEM 模擬中t=9.6 s 時的巖體運動一致。在T=21.2 s 之后，甑子巖崩塌已基本完成，在影像中僅剩崩塌形成的沖擊氣浪仍在擴散，與模擬中的t=26.8 s 時的情況相同。通過對比影像資料與MatDEM 模型可證明離散元模擬的有效性，該模型可以反映甑子巖崩塌的動力特征。

圖7 甑子巖崩塌現(xiàn)場影像資料Fig.7 Video data of the Zengziyan rockfall

3.2 崩塌動力破碎演化

甑子巖崩塌過程中碎屑等效粒徑分布演化如圖8所示。當t=1.6 s 時，崩塌體底部內側相對較完整的巖體破碎成為基本單元，中上部巖體沿預設節(jié)理面破裂明顯，完整的巖體被分割為不同等效粒徑的巖體，等效粒徑范圍為8～21 m。t=3.2 s 時，崩塌體開始運動至一級陡崖外，中上部巖體自下而上開始破裂解體，等效粒徑降低至0～8 m。t=4.8 s 時，中上部巖體沿節(jié)理裂隙面破裂、墜落，因與低速三角區(qū)域摩擦、撞擊導致顆粒破碎明顯，等效粒徑降至4 m 左右，凌空飛行的巖體大部分已破碎為基本顆粒。

圖8 甑子巖崩塌巖塊粒徑演化Fig.8 Size evolution of the Zengziyan rockfall's fragments

t=6.4 s 時，巖體頂部基本未發(fā)生破碎，巖體沿節(jié)理面已明顯分離，基本顆粒單元包裹大粒徑巖體凌空飛行。t=8.0 s 時，豎向節(jié)理右側的巖塊已滑離一級陡崖，等效粒徑維持在8 m 左右，凌空飛行的中部巖體等效粒徑降至2 m。當t=9.6 s 時，除低速三角區(qū)域外崩塌源區(qū)已滑出一級陡崖，頂部巖體等效粒徑仍維持在20 m 左右，凌空飛行的巖塊與斜坡堆積區(qū)產(chǎn)生的碰撞對巖塊造成了再次破碎，大部分等效粒徑基本已降低至基本顆粒單元，少部分巖塊等效粒徑為2 m 左右。t=12.8 s 時，崩塌源區(qū)大部分巖體已形成堆積，仍有等效粒徑為4 m、8 m 的巖塊存在，頂部巖體發(fā)生破碎粒徑降低至0～13.5 m，根據(jù)圖6，此時顆粒仍具有較大動能，顆粒將在堆積區(qū)進一步撞擊破碎。t=26.8 s 時，崩塌運動結束，堆積區(qū)形成，大部分巖體已破碎為基本顆粒單元，但在減速堆積過程中巖體破碎不明顯，大粒徑巖塊基本集中在堆積體表面，呈現(xiàn)反粒序堆積分布特征。

巖體崩塌基本單元增長率與最大平均速率如圖9所示，基本單元增長率峰值點與甑子巖崩塌各階段密切相關，由此可確定甑子巖崩塌過程中存在四個顯著的破碎時刻，即基本單元增長率中的4 個明顯峰值點，其分別是崩塌源區(qū)底部巖體受壓破碎、中上部巖體撞擊低速三角區(qū)、中部巖體撞擊斜坡地面與上部巖體撞擊斜坡地面。因顆粒破碎原因不一致，因此基本單元增長率曲線與最大平均速率曲線變化規(guī)律并不完全一致。圖中最大平均速率曲線僅有一個峰值，是因為在崩塌過程中巖體首先進入受底部巖體約束的墜落狀態(tài)，之后滑出一級陡崖進入凌空飛行狀態(tài)，單元速率持續(xù)增加，直至巖體撞擊斜坡地面開始減速，因此最大平均速率最值點發(fā)生于大部分巖體在凌空飛行時的階段。

圖9 崩塌基本單元增長率與最大平均速率曲線Fig.9 Curve of growth rate of basic elements and maximum average velocity of the Zengziyan rockfall

3.3 分維特征與Weibull 曲線擬合

甑子巖崩塌前后巖塊質量分布如圖10所示。崩塌前后擬合直線的斜率分別為0.95 與0.88，R2值均為0.96，巖塊質量分布基本符合公式3 獲得的直線。通過計算得到崩塌前后巖塊的分維值分別為2.05 與2.12，分維值越高表示巖塊的破碎情況越顯著，由此可知甑子巖崩塌后巖塊破碎明顯。

圖10 甑子巖崩塌巖塊質量分布Fig.10 Mass distribution of rock block of the Zengziyan rockfall

巖體崩塌前后等效粒徑級配曲線如圖11所示，采用了雙參數(shù)Weilbull 分布模型進行曲線擬合，崩塌前后尺度參數(shù)dc由15.839 3 降低至3.767 5，說明巖塊尺寸明顯降低，細粒巖塊在粒徑分布占比中明顯增長。崩塌前后形狀參數(shù) β由1.153 升高至1.436，說明巖塊破碎后等效粒徑范圍明顯減小，崩塌后等效粒徑大于5 m 的巖塊明顯減少，巖塊粒徑主要集中在0～2.7 m。

圖11 崩塌堆積巖塊級配曲線Fig.11 Grain size distribution of rock blocks of the rockfall

4 討論

關于甑子巖崩塌的研究之前多集中于其啟動機理與穩(wěn)定性分析[5-6,8-9]，其運動過程的研究多為基于影像資料分析其運動特征[7]，文章實現(xiàn)了對甑子巖危巖體崩塌運動全過程的離散元模擬，模擬的崩塌過程與影像資料相吻合，并引入分形維數(shù)與雙參數(shù)Weibull 分布模型對堆積顆粒進行了分析。模擬崩塌過程中危巖體底部低速三角區(qū)域的出現(xiàn)及其形態(tài)特征與許多學者在室內試驗中觀測到的基本一致[16-17]，崩塌堆積區(qū)最終呈現(xiàn)反粒序堆積分布特征[18]。

許多學者將柱狀巖體崩塌問題簡化為顆粒柱崩塌進行研究，基于室內試驗與數(shù)值模擬研究柱體高寬比、顆粒間粘結狀態(tài)對運動過程與堆積分布的影響[2,16,19]。但如甑子巖危巖體、重慶箭穿洞危巖體、新疆蓋孜河危巖群體等一般位于高處[20-21]，其崩塌過程常伴隨凌空飛行階段與落地后撞擊斜坡發(fā)生二次破碎階段，破碎顆粒能量與粒徑均會發(fā)生較大改變，因此，之后可進一步研究顆粒柱高程、斜坡角度及坡面摩擦系數(shù)等對柱體崩塌破碎影響。一些學者對實際柱狀危巖體崩塌進行了數(shù)值模擬研究，但其關注點多在于崩塌運動過程的分析，較少引入顆粒單元破碎、分維特征等內容對崩塌堆積體進行分析，文章補充了這方面的空白，在未來也可通過現(xiàn)場調查甑子巖堆積體分布特征進行進一步優(yōu)化數(shù)值模擬模型[22-23]。

5 結論

（1）文中基于MatDEM 離散元軟件，實現(xiàn)了對按照真實節(jié)理分布的甑子巖W12 危巖體崩塌全過程模擬，并通過與影像資料對比驗證了數(shù)值模擬的有效性。

（2）通過對MatDEM 二次開發(fā)，統(tǒng)計分析了甑子巖崩塌全過程等效粒徑演化與破碎規(guī)律，確定了崩塌過程中巖體四個顯著的顆粒破碎時刻，分別是崩塌源區(qū)底部巖體受壓破碎、中上部巖體撞擊低速三角區(qū)、中部巖體撞擊斜坡地面與上部巖體撞擊斜坡地面。

（3）引入了分形維數(shù)與雙參數(shù)Weibull 分布模型對甑子巖崩塌前后進行分析，崩塌后分維值D增加至2.12，尺度參數(shù)dc降低至3.767 5，形狀參數(shù) β增加至1.436，說明甑子巖崩塌后巖塊破碎明顯，細粒巖塊占比明顯增長，粒徑分布范圍明顯減小。