劉康琦 劉紅巖 祁小博

(①中國地質(zhì)大學(xué)(北京)工程技術(shù)學(xué)院， 北京 100083， 中國) (②自然資源部深部地質(zhì)鉆探技術(shù)重點實驗室， 北京 100083， 中國) (③中國地質(zhì)環(huán)境監(jiān)測院， 北京 100081， 中國)

0 引 言

土石混合體邊坡是第四紀(jì)以來形成的、由土和塊石組成的一種極度不均勻的巖土體邊坡。不同于一般的土質(zhì)邊坡或巖質(zhì)邊坡，由于塊石與土的共同存在，使土石混合體邊坡具有極大的不均勻性和不連續(xù)性。目前針對土石混合體邊坡的穩(wěn)定性已有較多研究，油新華等(2002)通過對土石混合體進行野外水平推剪試驗，得出土石混合體的變形特點和相關(guān)的抗剪強度參數(shù)，并發(fā)現(xiàn)土石混合體在剪切破壞時會繞過大的石塊。徐文杰等(2009)將土石混合體定義為第四紀(jì)以來形成的，由具有一定工程尺度、強度較高的塊石、細(xì)粒土體及孔隙構(gòu)成且具有一定含石量的極端不均勻松散巖土介質(zhì)系統(tǒng)，并通過可視粒徑、土石閾值、“土”與“石”的強度特征及含石量4個參數(shù)對土石混合體的概念進行了補充。Zhou et al.(2019)利用室內(nèi)循環(huán)凍融試驗和PFC3D數(shù)值模擬研究了凍融循環(huán)對土石混合體強度的影響。徐文杰(2009)基于對周家灣滑坡的大量地質(zhì)調(diào)查研究，對土石混合體滑坡體進行了工程地質(zhì)分區(qū)，并建立了三維地質(zhì)模型，運用有限元強度折減法分析了其穩(wěn)定性。楊忠平等(2017)基于PFC2D接觸黏結(jié)模型的離散元數(shù)值模型，研究了含石量變化對土石混合體剪切特性的影響。龔健等(2017)通過對不同含石量的土石混合體進行靜力超載試驗，得出土石混合體的力學(xué)行為與含石量的關(guān)系，另基于PIV分析方法，研究了剪切破壞時剪切帶的發(fā)展規(guī)律。張森等(2016)提出基于顆粒流區(qū)域填充技術(shù)的土石混合體有限元隨機生成技術(shù)，并利用強度折減法研究了含石量對邊坡穩(wěn)定性的影響規(guī)律。Sun et al.(2014)通過對冷谷水電站馬河堆積體進行工程類比、顆粒分析試驗、直剪試驗、現(xiàn)場大型推剪試驗、數(shù)值模擬試驗、參數(shù)反演等研究，對比分析了土石混合體的抗剪強度，并對馬河堆積體的穩(wěn)定性進行了分析。Fan et al.(2015)使用離散元數(shù)值模擬方法，研究了土石混合體在循環(huán)荷載作用下的力學(xué)特性，結(jié)果表明峰值應(yīng)力會隨著動應(yīng)力幅值和頻率的增大而增大，但會隨著塊石所占比例的增大而減小，同時彈性模量會隨著動應(yīng)力幅值和頻率的增大而減小。鄭博寧等(2019)提出了基于CT掃描技術(shù)和PFC三維顆粒流程序的土石混合體三維建模技術(shù)，并通過與三軸剪切試驗過程中礫石的空間運動規(guī)律進行對比分析，驗證了其建模方法的可行性。Wang et al.(2019)利用離心模型試驗通過考慮不同體積塊石含量，研究了土石混合體邊坡在地表荷載作用下的破壞行為。Lu et al.(2018)使用離散元法研究了含石量、坡高比和塊石大小對土石混合體邊坡穩(wěn)定性的影響。Yang et al.(2019)利用修正數(shù)值流形法研究了連續(xù)開挖條件下土石混合體邊坡的穩(wěn)定性。

然而目前對土石混合體邊坡穩(wěn)定性的研究，基本都忽略了土體的蠕變效應(yīng)。研究表明(陳衛(wèi)兵等， 2008； 蔣海飛等， 2013； 陳廷君等， 2019)，土體蠕變會降低坡體的穩(wěn)定性系數(shù)，對邊坡穩(wěn)定性造成不利影響。本文首先利用數(shù)字圖像處理技術(shù)對我國某水電站庫區(qū)的一個土石混合體邊坡進行處理，最大化的獲得土石混合體邊坡的細(xì)觀模型，然后利用FLAC3D程序中強度折減法對其進行穩(wěn)定性分析，并考慮了土體的蠕變性質(zhì)，對比分析了考慮土體蠕變的土石混合體邊坡與不考慮土體蠕變的土石混合體邊坡的變形特征。然后通過對采用不同土體蠕變參數(shù)的土石混合體邊坡進行穩(wěn)定性分析，討論了土體蠕變參數(shù)對土石混合體邊坡穩(wěn)定性的影響。

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

1.1 數(shù)字圖像處理

數(shù)字圖像處理(Digital Image Processing)被最早用來改善圖像的質(zhì)量，而隨著計算機科學(xué)的不斷發(fā)展，其開始廣泛應(yīng)用于航空航天、醫(yī)學(xué)、軍事、工程等領(lǐng)域。近年來，數(shù)字圖像處理技術(shù)也越來越多地被應(yīng)用到巖土工程界(Chen et al.， 2004； Xu et al.， 2008; 徐文杰等， 2008； 廖秋林等， 2010； 王學(xué)濱， 2018)，它將圖像經(jīng)過數(shù)學(xué)變換后轉(zhuǎn)換為數(shù)字信號并利用計算機對其進行分析處理，可將圖像定義為一個非零的和有限的二維函數(shù)f(x，y)，其中(x，y)表示空間或平面坐標(biāo)，f表示圖像在該點處的灰度值。在對土石混合體邊坡進行處理時，可根據(jù)土與塊石灰度值的不同，通過設(shè)定閾值，將土與石塊分別開來，這樣利于提取土石混合體邊坡中土與塊石的分布信息。

在進行圖像處理時依托于OpenCV強大的圖像處理功能及Python語言的簡潔高效，OpenCV是由C++語言編寫的一種計算機視覺庫(Bradski， 2000)，具有輕量、高效且開源的優(yōu)點，并提供了多種語言的接口。

圖 1 實際的土石混合體邊坡圖像Fig. 1 Image of soil-rock mixture slope

圖 1為我國某水電站庫區(qū)的一個土石混合體邊坡的斷面圖像。對其具體處理時過程如下：

(1)新建.py文件。與需要處理的土石混合體邊坡圖像放置在同一文件夾下。

(2)邊坡圖像二值化處理。依次對圖像進行均值濾波器去噪處理→圖像灰度化處理，將原來的RGB彩色圖像轉(zhuǎn)化為灰度圖像→灰度圖像高斯濾波器去噪處理，以消除灰度圖內(nèi)的高斯噪聲→自適應(yīng)閾值二值化處理，可在圖像局部范圍內(nèi)自適應(yīng)一個閾值對圖像進行二值化處理，將圖像上石塊的灰度值設(shè)置為0，呈黑色，土體的灰度值設(shè)置為255，呈白色，使用自適應(yīng)閾值二值化可最大化消除光照等因素對圖像造成的影響→腐蝕算法去噪處理。

(3)細(xì)節(jié)調(diào)整。使用Photoshop軟件對二值化后的邊坡圖像內(nèi)不理想的區(qū)域進行手工處理，使二值化圖像最大化接近真實土石分布情況。

(4)光滑及矢量化處理。由于生成的二值化圖像內(nèi)石塊邊界粗糙會對數(shù)值模擬造成困難甚至?xí)箶?shù)值模擬失敗，因此有必要對其進行邊界光滑處理。另數(shù)字圖像并不能直接用于數(shù)值模擬，因此需要對其進行矢量化處理，生成AutoCAD可接受的DXF文件。圖 2即為對圖 1進行矢量化處理后的結(jié)果，土石比約為14︰1，最大粒徑為1.37im。

圖 2 矢量化土石混合體邊坡模型Fig. 2 Vector model of soil-rock mixture slope

1.2 建立模型

在進行數(shù)值計算時塊石距離邊坡邊界過近會產(chǎn)生邊界效應(yīng)(徐文杰等， 2008)，為減少由于邊界效應(yīng)產(chǎn)生的影響，在邊坡模型的右側(cè)及左側(cè)坡腳處使用均質(zhì)土體進行了延拓，在AutoCAD中面域化后生成DXF文件導(dǎo)入ANSYS軟件進行網(wǎng)格劃分，然后轉(zhuǎn)化為可供FLAC3D軟件接受的文件格式。最終在FLAC3D中生成模型(圖 3)，共生成節(jié)點10i906個，劃分網(wǎng)格單元5455個，其中土體單元4736個，塊石單元719個。

圖 3 FLAC3D中生成的土石混合體邊坡模型圖Fig. 3 Model of soil-rock mixture slope in FLAC3D

2 考慮蠕變性質(zhì)的強度折減法

2.1 強度折減法

由于傳統(tǒng)的靜力平衡法需要預(yù)先指定滑動面，不適于內(nèi)部塊石分布復(fù)雜的土石混合體邊坡，因此采用強度折減法對邊坡進行穩(wěn)定性分析(Roth et al.， 1999)。強度折減法中將安全系數(shù)定義為邊坡剛好達(dá)到臨界失穩(wěn)狀態(tài)時，對巖土體的抗剪強度的折減程度。以莫爾-庫侖強度屈服準(zhǔn)則為例，采用強度折減法對邊坡進行穩(wěn)定性分析時，通過不斷的對邊坡的黏聚力c和內(nèi)摩擦角φ除以折減系數(shù)進行試算，直到邊坡達(dá)到失穩(wěn)狀態(tài)，此時的折減系數(shù)即是該邊坡的安全系數(shù)，其折減過程如下：

(1)

(2)

式中，ctrial、φtrial為土體折減后的黏聚力、內(nèi)摩擦角；Ftrial為折減系數(shù)；c、φ為巖土體折減前的黏聚力、內(nèi)摩擦角。

2.2 蠕變本構(gòu)模型的選取

在對蠕變邊坡進行穩(wěn)定性分析時，選取合適的黏彈塑性蠕變本構(gòu)模型及參數(shù)可充分反映蠕變對邊坡變形及穩(wěn)定性的影響，直接關(guān)系到分析結(jié)果的合理性與正確性。一般來說，在選取蠕變模型時，應(yīng)先進行室內(nèi)試驗，在試驗數(shù)據(jù)的基礎(chǔ)上利用蠕變理論選取合適的本構(gòu)模型，然后利用回歸分析和最小二乘法得出模型的參數(shù)。本文選取FLAC3D中Burgers黏彈性模型與Mohr-coulomb模型串聯(lián)而成的Cvisc模型，其構(gòu)成如下(圖4)：

圖 4 FLAC3D中Cvisc模型示意圖Fig. 4 Illustrations of Cvisc rheological model in FLAC3D

在FLAC3D中，Cvisc模型由Kelvin體、Maxwell體和一個塑性體構(gòu)成，其中，σ表示作用在巖土體上的應(yīng)力；σt表示巖土體的屈服強度；EM、EK分別表示彈性模量和黏彈性模量；ηM、ηK分別表示Maxwell黏性參數(shù)和Kelvin黏性參數(shù)。在巖土體的蠕變變形中，若其具有穩(wěn)定蠕變或加速蠕變性質(zhì)，則即使不對巖土體的黏聚力和內(nèi)摩擦角進行折減，邊坡也會一直變形直至破壞，這樣則無法體現(xiàn)邊坡的變形失穩(wěn)與巖土體強度之間的關(guān)系。因此需要將巖土體設(shè)定為只具有衰減蠕變性質(zhì)，在FLAC3D中，只需將Maxwell黏性系數(shù)ηM設(shè)置為無限大即可，此時，Cvisc模型則退化為廣義Kelvin黏彈塑性模型。

三維形式下模型應(yīng)變增量偏量由Maxwell體、Kelvin體和塑性體3部分組成。分別為：

(3)

(4)

(5)

(6)

2.3 邊坡失穩(wěn)判據(jù)的選取

在使用強度折減法對邊坡進行穩(wěn)定性分析時，需要選取合適的失穩(wěn)判據(jù)來判斷邊坡是否達(dá)到臨界失穩(wěn)狀態(tài)。傳統(tǒng)的強度折減法中，一般有3種方法用來作為判據(jù)：

(1)數(shù)值計算不收斂。當(dāng)巖土體的強度參數(shù)折減到一定程度使邊坡達(dá)到失穩(wěn)破壞狀態(tài)時，F(xiàn)LAC3D程序無法得到一個確定的解來滿足靜力平衡和應(yīng)力-應(yīng)變關(guān)系。

(2)關(guān)鍵點位移突變。可以通過監(jiān)測若干關(guān)鍵點的位移來判定邊坡是否破壞，當(dāng)強度折減到一定程度使邊坡達(dá)到臨界破壞狀態(tài)時，位移會突然增大。

(3)具有貫通的塑性區(qū)。當(dāng)塑性區(qū)從坡頂一直貫通到坡腳時，可認(rèn)為邊坡已達(dá)到失穩(wěn)狀態(tài)。在FLAC3D中，可通過查看剪應(yīng)變速率云圖來觀察塑性區(qū)是否貫通。

以上3種方法即為傳統(tǒng)的強度折減法中邊坡是否失穩(wěn)的判據(jù)，而在考慮蠕變的土石混合體邊坡穩(wěn)定性分析中并不適用。一方面，由于塊石的存在，很可能使邊坡在破壞時不會出現(xiàn)明顯貫通的塑性區(qū)； 另一方面，在考慮巖土體的蠕變性質(zhì)時，很難以數(shù)值收斂作為邊坡是否失穩(wěn)破壞的判據(jù)。因此，可將關(guān)鍵點位移在經(jīng)過一定時間后是否穩(wěn)定及位移在強度折減到某個程度時是否突變作為土石混合體邊坡是否失穩(wěn)的判據(jù)(蔣海飛等， 2013)。

2.4 模型驗證

為驗證本文方法的正確性，使用自編蠕變強度折減法對陳衛(wèi)兵等(2008)的算例進行對比。在對算例進行試算時，取土體參數(shù)如下：重度γ=20ikN·m-3，黏聚力c=35 kPa，內(nèi)摩擦角φ=18°，剪脹角Ψ=18°，彈性模量和黏彈性模量分別為E=50 MPa、Ek=300 MPa，黏性系數(shù)ηk=5 GPa·d， 泊松比v=0.￣35，蠕變時間為1ia。通過不斷折減土體強度參數(shù)，最后求得在考慮蠕變效應(yīng)時該算例的安全系數(shù)為1.48，這與陳衛(wèi)兵等(2008)的結(jié)果一致，因此可以認(rèn)定本文所使用自編強度折減法是可靠的。圖 5為該算例與使用自編強度折減法取折減系數(shù)為1.48時的剪應(yīng)變速率云圖。

圖 5 折減系數(shù)為1.48時剪應(yīng)變速率云圖Fig. 5 Contours of shear strain rate at reduction factor=1.48a. 本文計算結(jié)果； b. 文獻(xiàn)計算結(jié)果(陳衛(wèi)兵等，2008)

3 土石混合體邊坡穩(wěn)定性分析

計算模型的邊界條件為：模型右側(cè)及左側(cè)邊界采用水平約束，邊坡底面采用水平、垂直全部約束。

3.1 考慮蠕變的土石混合體邊坡穩(wěn)定性分析

為作對比，首先研究不考慮蠕變效應(yīng)時的土石混合體邊坡穩(wěn)定性，土體和塊石均采用Mohr-Coulomb模型?；谑Y海飛等(2013)的研究，本文所選巖土體模型參數(shù)如表 1所示。

表 1 不考慮蠕變時模型參數(shù)表Table 1 Physico-mechanical parameters of S-RMS without considering rheology

圖 6 不考慮蠕變效應(yīng)時邊坡剪應(yīng)變速率云圖Fig. 6 Contours of shear strain rate without considering rheology

表 2 土體蠕變模型參數(shù)表Table 2 Physico-mechanical parameters of soil with considering rheology

圖 6為不考慮蠕變效應(yīng)時土石混合體邊坡強度折減法的分析結(jié)果，可以看到，由于塊石的存在，在邊坡內(nèi)部形成了多條滑帶，而且滑帶在靠近塊石的部位有明顯的“繞石”現(xiàn)象，由強度折減法得出的折減系數(shù)為1.85，即不考慮蠕變效應(yīng)時此土石混合體邊坡的安全系數(shù)為1.85。

在對考慮蠕變效應(yīng)的土石混合體邊坡進行穩(wěn)定性分析時，在可能的滑帶上及坡頂坡腳處選擇關(guān)鍵點來監(jiān)測其水平位移情況(陳衛(wèi)兵等， 2008)。石塊采用Mohr-Coulomb模型，土體采用Cvisc黏彈塑性蠕變模型，蠕變時間為1ia。其Cvisc模型參數(shù)如表 2所示。

圖 7 考慮蠕變效應(yīng)時邊坡剪應(yīng)變速率云圖Fig. 7 Contours of shear strain rate with considering rheology

圖 8 折減系數(shù)為1.46時監(jiān)測點位移曲線Fig. 8 Horizontal displacement vs time of monitoring points at reduction factor=1.46

圖 9 折減系數(shù)為1.47時監(jiān)測點位移曲線Fig. 9 Horizontal displacement vs time of monitoring points at reduction factor=1.47

圖 10 關(guān)鍵點水平位移與折減系數(shù)的關(guān)系曲線Fig. 10 Horizontal displacements of monitoring points at different reduction factors with considering rheology

通過自編強度折減法來對考慮蠕變效應(yīng)的土石混合體邊坡進行穩(wěn)定性分析，其中蠕變參數(shù)不變，折減系數(shù)從1.30增加到1.47。圖 7為折減系數(shù)為1.47時的剪應(yīng)變速率云圖，圖 8和9為折減系數(shù)分別是1.46和1.47時監(jiān)測點水平位移曲線，可以看到，當(dāng)折減系數(shù)為1.46時，各監(jiān)測點水平位移在經(jīng)過短暫的瞬時變形后，位移增加量逐漸減小，最終趨于穩(wěn)定，位移增加量近似為0。當(dāng)折減系數(shù)增加到1.47時，雖然監(jiān)測點3和監(jiān)測點4的位移最終趨于穩(wěn)定，但監(jiān)測點1和監(jiān)測點2在經(jīng)歷1ia的蠕變時間之后，位移仍未收斂。

各關(guān)鍵點水平位移與強度折減系數(shù)的關(guān)系曲線如圖 10所示，可以看出，折減系數(shù)由1.30增加到1.46時，關(guān)鍵點總位移并未發(fā)生太大變化，而當(dāng)折減系數(shù)增加到1.47后，可以明顯地看到位移曲線有一個突變，位移增大較為明顯，而此時其蠕變曲線也未到達(dá)穩(wěn)定狀態(tài)，因此可以認(rèn)定，考慮蠕變時該土石混合體邊坡的安全系數(shù)為1.47，比不考慮蠕變時該土石混合體邊坡的安全系數(shù)下降了20.5%。以監(jiān)測點2為例，取折減系數(shù)為1.46，不考慮土體蠕變效應(yīng)時，其水平位移僅為0.002im，而考慮土體蠕變效應(yīng)時，水平位移為0.32im，是不考慮蠕變效應(yīng)時的160倍，可見土體蠕變對邊坡變形影響極大； 在不考慮蠕變效應(yīng)的情況下取折減系數(shù)為1.84，監(jiān)測點2的水平位移為0.15im，因此若將0.15im作為預(yù)測邊坡失穩(wěn)的位移量則顯得過于保守。

圖 11 不同黏性系數(shù)時邊坡安全系數(shù)Fig. 11 FOS with different viscosity coefficient

3.2 蠕變參數(shù)對安全系數(shù)的影響

巖土體的黏性系數(shù)是反映其蠕變程度的一個重要參數(shù)，為探索黏性系數(shù)對土石混合體邊坡安全系數(shù)的影響關(guān)系，在其他力學(xué)參數(shù)不變的情況下，分別取Cvisc模型中的Kelvin黏性系數(shù)為1.19e3iPa·d、1.19e4iPa·d、1.19e5iPa·d和1.19e6iPa·d，并對邊坡進行穩(wěn)定性分析。圖 11為取不同黏性系數(shù)時的邊坡安全系數(shù)關(guān)系圖，由圖可以看出，黏性系數(shù)對于邊坡的變形及穩(wěn)定影響較大，隨著黏性系數(shù)的增大，邊坡安全系數(shù)隨之減小，說明黏性系數(shù)的增大對于邊坡的變形及穩(wěn)定具有不利影響。

4 結(jié) 論

土石混合體邊坡是自然界中廣泛存在的一種巖土體邊坡，以我國某水電站庫區(qū)的一個土石混合體邊坡為例，首先運用數(shù)字圖像處理技術(shù)對該邊坡進行處理，得到土石混合體邊坡的細(xì)觀物理模型，然后在FLAC3D中運用強度折減法對該邊坡進行了穩(wěn)定性分析，并考慮了土體的蠕變性質(zhì)，與未考慮蠕變性質(zhì)的土石混合體邊坡穩(wěn)定性分析結(jié)果進行比較，得出以下結(jié)論：

(1)土石混合體邊坡不同于一般均質(zhì)土體或巖質(zhì)邊坡，在達(dá)到臨界失穩(wěn)狀態(tài)時，在邊坡內(nèi)部會出現(xiàn)多條滑帶，且滑帶具有明顯的“繞石”現(xiàn)象。

(2)不考慮土體蠕變性質(zhì)的土石混合體邊坡的安全系數(shù)為1.85，而考慮了土體蠕變性質(zhì)的土石混合體邊坡的安全系數(shù)為1.47，相比下降了20.5%； 在取折減系數(shù)同為1.46的情況下，考慮土體蠕變效應(yīng)的邊坡監(jiān)測點水平位移是不考慮蠕變效應(yīng)時的160倍之多。因此巖土體的蠕變性質(zhì)不利于邊坡的穩(wěn)定，在對巖土體邊坡進行穩(wěn)定性分析時不應(yīng)忽略其蠕變性質(zhì)。

(3)研究了Kelvin黏性系數(shù)對邊坡進行長期穩(wěn)定性的影響，認(rèn)為邊坡安全系數(shù)隨著Kelvin黏性系數(shù)的增大而減小，說明黏性系數(shù)對于巖土體邊坡的變形及穩(wěn)定具有不利影響。