土石混合體邊坡模型隨機構建與穩定性分析

賀 勇, 張宗堂,楊期君,張巨峰,李文菁

（1. 長沙市規劃勘測設計研究院，長沙 410000； 2. 湖南科技大學巖土工程穩定控制與健康監測湖南省重點實驗室，湖南 湘潭 411201； 3. 湖南科技大學資源環境與安全工程學院，湖南 湘潭 411201）

土石混合體（soil rock mixture）是指由具有一定尺寸且強度較高的塊石、強度較低的土體及孔隙構成的極端不均勻的松散巖土體，是一種在自然界廣泛存在的典型非均質巖土材料[1-5]；Medley也稱之為bimrock（block-in-matric rock），意思是“由較細結構的粘結基質與具有巖土工程意義的塊體組成”，其中“具有巖土工程意義的塊體”表示塊體與弱基質之間存在足夠的力學對比，并且塊體的尺寸分布和體積比例會影響工程規模下的巖體性質[6,7]。在巖土工程領域，土石混合體材料十分常見，例如：川藏公路林芝至八宿段全長427 km范圍內，地質災害的主要物質是土石混合體[8]；浙江省上（虞）三（門）高速公路K92-K95段的邊坡為土石混合體，滑坡治理費用高達億元[9]；攀西地區發生的816個滑坡中，土石混合體邊坡就有500個，占61.3%[3,10]。

考慮到這些土石混合體的巨大空間與力學變異性，以往學者們為了簡化計算，常常選擇使用較弱強度和變形特征好模擬的土壤來開展工作[11-12]，忽略了塊體對整體土石混合體強度的貢獻。然而，正如許多文獻所報道的一樣，這種簡化可能導致較大的計算誤差。也有實踐表明，由于土石混合體本身的多向性及不連續性等特征，使用傳統的均質巖土材料假定應用于土石混合體巖土構筑物的分析存在明顯不足[13]。在目前的研究中缺乏與自然界真實情況符合程度較高的土石混合體模型。

在過去的幾十年里，人們進行了大量的研究，以建立正確描述土石混合體特征的系統方法。學者們大都致力于選擇適當的強度和變形參數，并進行恰當的數值模擬，以便在復雜地層中正確地開展工程活動。在土石混合體室內試驗結果的基礎上，一些學者提出了初步的強度準則：假定土石混合體是均勻各向同性的且具有等效的力學性能，則土石混合體的強度可根據其塊石含量和土體強度參數確定[14]。還有一些學者分析了在用確定性方法分析非均質地層中邊坡穩定性時塊石的影響[15]，他們強調了這種土石混合體建模方法產生的特殊結果，與分析均質土體時得到的結果有顯著不同。

基于前人的研究成果，本文提出了一種土石混合體的隨機生成算法。首先采用數字圖像處理來獲取自然界的真實塊石輪廓，然后通過MATLAB編程來生成具有特定級配和塊石含量土石混合體邊坡模型。最后分別基于FLAC2D與OPTUM G2對其進行穩定性分析；為了提高結果的統計有效性，對于每種含石量的邊坡（含石量0%的邊坡除外）進行了5次平行試驗，并對兩種方法的計算結果進行了對比分析。

1 土石混合體邊坡模型的構建

1.1 參數說明

本文旨在探討含石量對土石混合體邊坡穩定性的影響，為此，進行了大量的二維邊坡穩定性分析。對于簡單的邊坡模型，分別采用了有限差分法（FDM）和極限分析有限元法（FELA）對邊坡模型進行分析。其中，所有邊坡模型的坡度為45°，邊坡尺寸如圖1所示，級配均一致（圖2），但塊體含量不同。分別計算了塊石含量為0%（均值土體）、20%、40%和60%的土石混合體邊坡模型。由于它們是二維模型，塊體的面積即被認為等于塊石的體積。

圖1 邊坡尺寸示意圖（單位：m）Fig. 1 Schematic diagram of slope size (unit: m)

圖2 級配曲線Fig. 2 Grading curve

為了考慮塊石在空間和尺寸上的變化，作者編寫了一個MATLAB程序，用于在邊坡區域中生成隨機位置和大小的塊石顆粒。對于每個給定的級配和含石量，共隨機構建了5個不同的邊坡模型，即共進行了5次穩定性分析。設置5組平行試驗的目的是消除隨機誤差，以增加結果的可靠性。

此外，還分析了含石量為0%的邊坡，以評估地質勘探人員在僅根據土體的強度和變形特性進行設計時可能造成的潛在誤差。表1顯示了分析中使用的塊石和土體參數[1-2,16-18]，兩種材料的力學參數取值符合要求（Erock/Ematrix＞2和tanφrock/tanφmatrix＞2）[6,19]。假定它們遵循莫爾-庫侖破壞準則，并具有理想彈塑性行為；根據Lindquist和Kalender準則[20]，假設塊石為均勻各向同性材料，并通過不同方法對其進行了數值模擬。總共對32個土石混合體邊坡分別采用FLAC2D（有限差分法）和OPTUM G2（極限分析有限元）進行了穩定性分析，評估并比較了安全系數。

表1 參數取值Table 1 Parameter values

2 邊坡模型隨機構建

2.1 數字圖像處理獲取塊石顆粒

通過對真實塊石顆粒進行拍照并進行數字圖像處理，來提取塊石顆粒的輪廓坐標。采用數碼相機對顆粒進行拍照，相機型號為：佳能600D，高級攝影系統（APS）框架和單鏡頭反射（SLR）數碼相機，1800萬像素APS-C CMOS傳感器（尺寸：22.3*14.9 mm）和EXPEED 2圖像處理芯片，1.6倍鏡頭轉換系數；鏡頭型號為：EF-S 18-135 mm，67 mm鏡頭；得到的原始塊石照片如圖3a所示。首先，利用MATLAB中的rgb2gray函數對塊石顆粒的二維照片進行灰度處理（圖3b）；然后，利用MATLAB中的fspecial和imfilter函數對圖像進行高斯模糊平滑處理（圖3c），降低了圖像的噪聲；使用MATLAB中的im2bw函數對照片進行二值化處理（圖3d）后，可以看到二值化處理后顆粒區域為黑色，背景區域為白色；但是顆粒上仍然有一些部分是白色，這是因為顆粒上面的顏色分布不均勻，表面有少量白色區域造成的，不過這并不會影響輪廓的提取。最后，使用MATLAB中的bwboundaries函數搜索邊界像素來提取顆粒的輪廓（圖3e）。

需要注意的是，輪廓提取后（圖3e）得到的顆粒坐標是像素坐標，意味著一個顆粒的輪廓點有幾千上萬個，這樣的塊石顆粒直接用于數值計算是不現實的，因為輪廓點太多會造成計算十分耗時。因此，需要對塊石的輪廓點進行簡化處理。常用的簡化方法是等角度隨機采樣（圖3e到圖3f），詳細過程如圖4所示。當采樣角度為θ°時，粒子可簡化為N（N=360/θ）個點（在本研究中θ=10°，N=36），最終得到簡化后的顆粒輪廓如圖3f所示。等角度隨機采樣的詳細過程參見文獻[21]。

圖3 數字圖像處理獲取塊石輪廓Fig. 3 Block stone outline acquired by digital image processing

圖4 等角度采樣示意圖Fig. 4 Schematic diagram of equal angle sampling

2.2 MATLAB編程構建邊坡模型

在考慮塊石存在的情況下，利用MATLAB編寫程序進行蒙特卡羅模擬，根據給定的塊石含量和級配曲線，在邊坡區域內生成任意直徑和位置的真實塊石顆粒。對于每個不同的含石量（僅20%、40%和60%，因為0%代表均值土體），隨機構建了5個土石混合體邊坡模型，即為了實現結果的統計有效性，進行了5次穩定性分析。

當計算得到的含石量（計算為所有塊體的面積總和除以邊坡模型的面積）對應于輸入中定義的含石量時，使用“while”循環停止提取，允許誤差為±2.5%。

為了最大限度地提高代碼性能，加快計算速度，塊石的放置是根據塊石粒徑的大小從大到小進行的。MATLAB代碼滿足的其他要求是：（1）塊石不能相互穿透，否則就沒有物理意義；（2）塊石不能與它們所在的斜坡模型的外部邊界相交。實際上，任何交叉點都會導致塊體體積的部分損失，從而導致對塊石含量的低估，這將與規定的塊石含量不符。為了達到這些目的，將兩個塊石之間以及塊石與外部斜坡邊界之間的最小距離設置為5 cm。最終構建的土石混合體模型如圖5所示（含石量40%，級配如圖2所示）。

圖5 含石量為40%的土石混合體邊坡模型示例Fig. 5 Example of soil rock mixture slope model with 40% stone content

MATLAB代碼的主要輸出結果是一個文本文件，其中包含每個塊石的輪廓坐標。為了導入FLAC2D和OPTUM G2，所有MATLAB輸出文件都轉換成腳本文件，以便在AutoCAD軟件中可視化，然后以DXF格式保存。

3 土石混合體邊坡模型穩定性分析

分別采用有限差分法(FDM)與極限分析有限元法(FELA)，對含石量為0%（均值土體）、20%、40%和60%的土石混合體進行穩定性分析。除均值土體外，每個含石量隨機構建了5個邊坡進行分析，共16個邊坡，即采用兩種方法一共進行了32次穩定性分析。

3.1 有限差分法(FDM)分析

基于軟件FLAC2D對土石混合體邊坡進行有限差分法分析。首先將前文中得到的DXF格式的邊坡模型導入MIDAS GTS中進行網格劃分，然后編寫一個MIDAS2FLAC2D的接口程序將模型導入到FLAC2D中。分別對不同含石量的邊坡采用FLAC2D內置的有限差分強度折減法進行穩定性分析，其計算結果如圖6所示，計算所得的安全系數如表2和圖7所示。

由圖6可以看到，均質土體的塑性區為一條類圓弧滑動帶，而土石混合體的塑性區則較為復雜。具體來說，滑動帶沿著塊石周圍進行擴展，呈現出“繞石效應”。這是由于塊石的強度相對土體較大，滑動帶無法貫穿塊石，因此只能繞過塊石進行擴展。同時可以看到，隨著含石量的增大，邊坡的滑動帶變得更為破碎和曲折，這與徐文杰等[1]的報導一致。從安全系數計算結果上來看，雖然在相同含石量下，不同組邊坡計算所得的安全系數存在一定差異，但是總體上，邊坡的安全系數隨著含石量的增大而增大，這一點可以在安全系數平均值上（表2和圖7）體現出來。同時可以發現，隨著含石量的增大，5組邊坡的安全系數的標準差也隨之增大，說明含石量越大，塊石分布的隨機性對邊坡安全系數的影響越大。

圖6 有限差分法（其中一組平行試驗）分析結果Fig. 6 Analysis results of finite difference method (one group of parallel tests)

圖7 邊坡的安全系數平均值Fig. 7 Average safety factor of slope

表2 有限差分法安全系數計算結果Table 2 Calculation results of safety factor by finite difference method

3.2 極限分析有限元(FELA)分析

基于軟件OPTUM G2對土石混合體邊坡及逆行極限分析有限元分析（取上限解）。將前文中得到的DXF文件導入G2中進行網格劃分，然后賦予材料參數（表1）與邊界條件（G2中的“標準邊界條件”）對不同含石量的土石混合體邊坡進行穩定性分析，其計算結果如圖8所示，計算所得的安全系數如表3和圖7所示。

由圖8可以看到，FELA分析得到的邊坡的破壞模式與圖6中（FDM分析）邊坡的破壞模式基本一致，即均質土體的塑性區為一條類圓弧滑動帶，而土石混合體的塑性區則較為復雜，且隨著含石量的增大，土石混合體邊坡的滑動帶變得更加破碎。由表3和圖7可以看出，邊坡的安全系數隨含石量的增大而增大，且含石量越大，5組邊坡的安全系數的標準差也越大。

表3 極限分析有限元法安全系數計算結果Table 3 Calculation results of safety factor of limit analysis finite element method

圖8 極限分析有限元法（其中一組平行試驗）分析結果Fig. 8 Analysis results of limit analysis finite element method (one group of parallel tests)

3.3 FDM與FELA的比較

由上述兩種計算方式對比可以看到，采用不同方法進行分析的土石混合體邊坡的破壞模式相差不大，但是安全系數存在一定的差異。具體來說：FELA方法計算得到的邊坡安全系數的平均值會略大于FDM方法計算得到的安全系數的平均值，且差值隨著含石量的增大而增大（圖7）。同時，對于同一個含石量，FELA計算得到的邊坡安全系數的標準差要大于FDM；這是由于在采用FELA方法進行分析時，我們采用的是極限分析有限元的上限解。所以，FELA計算出的結果總是偏安全的，即計算出的邊坡安全系數偏大。

4 結論

本文通過數字圖像處理獲取了自然界真實顆粒的數字化模型，然后基于MATLAB進行特定級配與含石量的土石混合體邊坡模型的構建，最后分別采用FDM與FELA方法對不同含石量的土石混合體邊坡進行了穩定性分析。結果表明：

（1）采用本文提出的土石混合體隨機生成方法，將塊石和土體分開考慮進行計算，可有效評估土石混合體邊坡穩定性。

（2）含石量對邊坡穩定性影響較大，具體表現為：隨著含石量的增大，土石混合體邊坡的滑動帶變得更為破碎和曲折；邊坡的安全系數與5組邊坡安全系數的標準差隨著含石量的增大而增大。

（3）對于土石混合體邊坡，FELA方法計算得到的邊坡安全系數的平均值大于FDM方法計算的平均值。