王述紅，朱承金，張紫杉，任藝鵬,王鵬宇，邱 偉

(東北大學(xué) 資源與土木工程學(xué)院，遼寧 沈陽 110819)

高邊坡由于具有規(guī)模大、地質(zhì)條件復(fù)雜、易受周邊環(huán)境干擾而破壞失穩(wěn)的特點(diǎn)，在初次失穩(wěn)后易發(fā)生二次失穩(wěn)，形成次級(jí)滑動(dòng)面。巖質(zhì)高邊坡受層理結(jié)構(gòu)的影響，使得邊坡的穩(wěn)定問題更加復(fù)雜和突出，故對(duì)潛在多滑面的巖質(zhì)高邊坡穩(wěn)定性進(jìn)行深入探究任重而道遠(yuǎn)[1-3]。

邊坡穩(wěn)定性分析中有2種常用的計(jì)算方法:極限平衡法(LEM)[4]和有限元法(FEM)[5]。與有限元相比，極限平衡法計(jì)算效率高，通過幾何假定并基于已知或假定的規(guī)則滑面求解安全系數(shù)，但其忽視了邊坡滑落過程中巖土體的本構(gòu)關(guān)系[6]。近年來有限元快速發(fā)展并將強(qiáng)度折減法結(jié)合，該技術(shù)不必預(yù)先假設(shè)滑動(dòng)面的位置，且易得安全系數(shù)及相應(yīng)潛在滑動(dòng)面，其局限性在于安全系數(shù)及滑動(dòng)面求解的單一性。CALA等[7]開創(chuàng)了新型強(qiáng)度折減法，在常規(guī)方法基礎(chǔ)上得出最危險(xiǎn)失穩(wěn)面，進(jìn)而確定首次滑落的步數(shù)Nr，繼續(xù)增大折減系數(shù)，隨即1.1Nr步將會(huì)被計(jì)算，相繼求得對(duì)應(yīng)次級(jí)潛在失穩(wěn)面及安全系數(shù)。李小春等[8]利用檢索邊坡體單元的方法，對(duì)不同折減系數(shù)下各范疇邊坡單元進(jìn)行統(tǒng)籌，以此得到多個(gè)滑落面及相應(yīng)安全系數(shù)。

非連續(xù)性是巖體固有的屬性，由SHI[9]首次提出的非連續(xù)變形分析(DDA)在模擬塊體大位移、大變形方面獨(dú)具優(yōu)勢(shì)，在邊坡工程領(lǐng)域也得到迅速發(fā)展。MACLAUGHLIN等[10]模擬了傾斜邊坡的平面和弧形破壞模式，發(fā)現(xiàn)其結(jié)果比常規(guī)分析方法具有更好的精度。FU等[11]將矢量和方法(VSM)應(yīng)用到DDA中，基于實(shí)際應(yīng)力狀態(tài)和矢量和算法計(jì)算邊坡穩(wěn)定安全系數(shù)。

發(fā)展邊坡穩(wěn)定性評(píng)價(jià)核心問題在于安全系數(shù)及相應(yīng)滑落面的求解。過去研究中的強(qiáng)度折減法大都是對(duì)巖體統(tǒng)一折減，在實(shí)際情況下，主要是結(jié)構(gòu)面的力學(xué)性質(zhì)決定了巖體的力學(xué)性質(zhì)，且?guī)r體結(jié)構(gòu)面破壞時(shí)并不是同時(shí)損傷同等程度;再之邊坡首次發(fā)生失穩(wěn)后塊體位置及應(yīng)力已重新分布，故應(yīng)使用首次失穩(wěn)結(jié)束后的邊坡模型對(duì)次級(jí)滑落面進(jìn)行分析，而以往對(duì)邊坡次級(jí)滑動(dòng)面的研究中鮮有考慮該實(shí)際狀況。為了解決該問題，筆者借鑒DDA計(jì)算位移的優(yōu)勢(shì)，考慮塊體間相對(duì)位移，提出基于DSR-DDA法的邊坡多滑面搜索并進(jìn)行邊坡穩(wěn)定性評(píng)價(jià)，利用傾斜平面滑塊經(jīng)典案例測(cè)試位移閾值，驗(yàn)證方法的可行性及計(jì)算精度，并嘗試將其應(yīng)用到撫順西露天礦巖質(zhì)高邊坡的穩(wěn)定性分析中。

1 DSR-DDA法

1.1 DDA基本原理

DDA利用一階位移函數(shù)來表述每個(gè)塊體的運(yùn)動(dòng)參量，并假設(shè)應(yīng)力和應(yīng)變是恒定的[12]。每個(gè)塊體的位移矢量包含6個(gè)變量:

D=(u0,v0,r0,εx,εy,γxy)T

(1)

其中，(u0,v0)為塊體內(nèi)特定點(diǎn)(x0,y0)的剛性位移;r0為塊體繞特定點(diǎn)(x0,y0)的旋轉(zhuǎn)角度;εx,εy,γxy分別為該塊體的正應(yīng)變與切應(yīng)變。塊體內(nèi)任意點(diǎn)(x,y)的位移(u,v)為

(2)

基于最小勢(shì)能原理構(gòu)建的系統(tǒng)整體平衡方程為

(3)

式中，系數(shù)矩陣中Kij為6×6子矩陣;Kii由塊體單元的材料屬性和幾何參數(shù)決定;Kij(i≠j)則由塊體i和塊體j間的接觸條件而決定;[Di]和[Fi]為6×1子矩陣;Di為塊體i的變形變量(d1i,d2i,d3i,d4i,d5i,d6i);Fi為塊體i上分配給6個(gè)變形變量的荷載。

1.2 DSR-DDA法位移閾值

為了考慮真實(shí)情況中塊體結(jié)構(gòu)面損傷的程度不一性，提出對(duì)位移變化不小于閾值的塊體結(jié)構(gòu)面剪切強(qiáng)度進(jìn)行折減變化，利用經(jīng)典斜坡滑塊案例測(cè)試剪切強(qiáng)度折減的位移閾值，模型如圖1所示，巖體參數(shù)見表1。在程序計(jì)算中，塊體結(jié)構(gòu)面抗剪強(qiáng)度參數(shù)每隔一定時(shí)步(10 000步，每一時(shí)步為0.001 s)按照塊體相對(duì)位移變化與閾值的關(guān)系只針對(duì)位移變化符合要求的結(jié)構(gòu)面進(jìn)行折減，折減系數(shù)從1.000開始，依次增加0.001，直至發(fā)生急劇變形，則此時(shí)其取值為該邊坡滑落面安全系數(shù)。該程序計(jì)算所模擬的折減過程在一定程度上與實(shí)際情況下巖體材料的損傷退化相符。

圖1 滑塊沿斜坡滑動(dòng)模型Fig.1 Model of sliding block along slope

該簡(jiǎn)易模型運(yùn)動(dòng)問題安全系數(shù)解析解為

(4)

該算例結(jié)構(gòu)面選取各內(nèi)摩擦角對(duì)應(yīng)解析解值見表2。圖2為結(jié)構(gòu)面內(nèi)摩擦角φ取不同值下，模型安全系數(shù)及安全系數(shù)與解析解的誤差百分比隨位移閾值的變化曲線，由圖2可得，不論內(nèi)摩擦角φ取何值，當(dāng)位移閾值接近1 mm時(shí)，安全系數(shù)取值接近最大值，且與解析解的誤差百分比有一個(gè)明顯的凹槽，說明這時(shí)的誤差與解析解相對(duì)最小且在0.5%以內(nèi)，所以位移閾值暫時(shí)取為1 mm，與解析解的結(jié)果保持一致。從而說明了DSR-DDA法是滿足計(jì)算要求的，并且能保證一定的計(jì)算準(zhǔn)確性。

表1 邊坡模型巖體參數(shù)Table 1 Mechanical parameters of the slope model

表2 安全系數(shù)解析解Table 2 Theoretical safety factor

圖2 安全系數(shù)及其誤差隨位移閾值變化曲線Fig.2 Curves of safety factor and its error with displacement threshold

塊體間接觸力的精確求解是解決總平衡方程的核心步驟。在每一時(shí)步內(nèi)，都要重新確定彈簧的施加與否及彈簧的位置，需要反復(fù)生成求解總剛度矩陣，剛性彈簧的施加與去除過程稱之為開-合迭代。接觸有3種狀態(tài):張開、滑動(dòng)和鎖定。模式變化的判定準(zhǔn)則見表3，其中，N為法向位移，N>0為張開;t為剪切位移矢量;f為摩擦力矢量；T為剪切位移;‖為兩個(gè)矢量方向相同。

基于以上分析，DSR-DDA法在程序中實(shí)現(xiàn)如圖3所示。

表3 接觸狀態(tài)Table 3 Contact status

圖3 DDA程序步驟流程Fig.3 Flow chart outlining the procedures of the DDA

圖4 邊坡全貌及地質(zhì)剖面Fig.4 Slope topography and geological section

2 撫順西露天礦巖質(zhì)邊坡模型及參數(shù)

2.1 撫順西露天礦巖質(zhì)邊坡概況

露天開采在煤炭資源開采中占據(jù)至關(guān)重要的地位[13-14]。撫順西露天礦地處撫順市中心，其邊坡穩(wěn)定性直接決定了礦坑與城市交界處的安全與否。由于受到渾河斷裂的牽引作用,向斜北翼地層遭受了強(qiáng)烈改造和破壞。向斜軸部呈圓弧型褶曲,軸面呈現(xiàn)曲面狀,向北傾斜,傾角30°左右軸跡方向?yàn)镹E60°,長(zhǎng)為20 km,兩翼間距為3 km。北翼產(chǎn)狀各異,傾角15°～60°,間或伴有小型褶曲。目前,采坑北幫邊坡滑坡、崩塌、地裂縫等地質(zhì)災(zāi)害仍持續(xù)不斷，不時(shí)產(chǎn)生大型滑坡。故選取北幫W區(qū)某標(biāo)段巖質(zhì)高邊坡作為案例。圖4為該標(biāo)段邊坡全貌及地質(zhì)剖面圖，為監(jiān)測(cè)邊坡實(shí)時(shí)位移，在鉆孔中布置深部位移測(cè)斜儀。

2.2 計(jì)算分析模型及參數(shù)

圖5 UAV技術(shù)Fig.5 UAV technology

節(jié)理、裂隙等不連續(xù)結(jié)構(gòu)面對(duì)巖體的變形起到控制作用，結(jié)構(gòu)面信息采集精度對(duì)數(shù)值模擬分析準(zhǔn)確性至關(guān)重要。運(yùn)用高精度、高效率的UAV(圖5)技術(shù)對(duì)確定性結(jié)構(gòu)面信息進(jìn)行非接觸精準(zhǔn)獲取，生成點(diǎn)云模型。假定結(jié)構(gòu)面平面方程為

D=AX+BY+C

(5)

其中，A，B，C，D為結(jié)構(gòu)面平面參數(shù)，提取結(jié)構(gòu)面所在平面的法向量n=(-A,-B,1)，并對(duì)結(jié)構(gòu)面進(jìn)行點(diǎn)的拾取，可得

(6)

借助最小二乘法解(A，B，C)，結(jié)構(gòu)面傾向及傾角分別為

(7)

式中，α，β分別為結(jié)構(gòu)面傾向和結(jié)構(gòu)面傾角。

對(duì)結(jié)構(gòu)面產(chǎn)狀進(jìn)行計(jì)算，并對(duì)提取出的結(jié)構(gòu)面產(chǎn)狀信息導(dǎo)出，部分計(jì)算結(jié)果見表4，結(jié)構(gòu)面產(chǎn)狀統(tǒng)計(jì)云圖如圖6所示。

研究表明結(jié)構(gòu)面產(chǎn)狀服從一定規(guī)律分布，對(duì)確定性結(jié)構(gòu)面分組處理，并依照雙正態(tài)密度分布生成隨機(jī)結(jié)構(gòu)面產(chǎn)狀，其標(biāo)準(zhǔn)變化形式為

(8)

式中，σα為確定性結(jié)構(gòu)面傾向的標(biāo)準(zhǔn)差;σβ為確定性結(jié)構(gòu)面傾角的標(biāo)準(zhǔn)差;ρ為結(jié)構(gòu)面傾向、結(jié)構(gòu)面傾角的相關(guān)系數(shù);μα為確定性結(jié)構(gòu)面傾向的均值;μβ為確定性結(jié)構(gòu)面傾角的均值。

表4 產(chǎn)狀計(jì)算結(jié)果Table 4 Yield calculation results

圖6 結(jié)構(gòu)面產(chǎn)狀統(tǒng)計(jì)云圖Fig.6 Statistical cloud chart of structure surface

將獲取的結(jié)構(gòu)面信息導(dǎo)入DDA程序中，結(jié)構(gòu)面網(wǎng)格切割巖體，形成DDA的塊體單元，即算例數(shù)值計(jì)算分析模型(圖7)。

圖7 撫順西露天礦北幫邊坡DDA數(shù)值計(jì)算分析模型Fig.7 DDA numerical calculation and analysis model of slope in West Fushun Open-Pit Mine

天然應(yīng)力場(chǎng)僅將重力作用納入范圍，不計(jì)算區(qū)域構(gòu)造應(yīng)力影響。根據(jù)地質(zhì)勘探和試驗(yàn)成果，得表5所示各項(xiàng)參數(shù)。程序計(jì)算中，塊體結(jié)構(gòu)面抗剪強(qiáng)度參數(shù)每隔一定時(shí)步(10 000步，每一時(shí)步約為0.001 s)針對(duì)塊體相對(duì)位移變化達(dá)到或超過位移閾值的結(jié)構(gòu)面進(jìn)行折減，折減系數(shù)從1.000開始，依次增加0.001，在一定程度上與實(shí)際情況下巖體材料的損傷退化相符。

表5 邊坡模型巖體參數(shù)Table 5 Mechanical parameters of the slope model

3 撫順西露天礦巖質(zhì)邊坡穩(wěn)定分析

3.1 現(xiàn)場(chǎng)位移監(jiān)測(cè)及數(shù)值計(jì)算分析

測(cè)斜孔布置參考圖4，測(cè)點(diǎn)間隔為1 m，近2個(gè)月內(nèi)各測(cè)斜孔不同深度測(cè)得的累計(jì)位移值如圖8所示。

圖8 測(cè)斜孔測(cè)點(diǎn)累積總位移值Fig.8 Cumulative displacements of monitoring point in surveying slant hole

分析圖8可得，測(cè)斜孔69002內(nèi)前100 m左右測(cè)點(diǎn)總位移值隨時(shí)間增長(zhǎng)變化極大，且在深度為100 m附近存在位移值劇增，說明測(cè)斜孔附近巖體已發(fā)生相對(duì)滑移，并有繼續(xù)滑動(dòng)的趨勢(shì)，且滑動(dòng)面在距坡面100 m附近。測(cè)斜孔55026內(nèi)前20 m左右測(cè)點(diǎn)總位移值增加相對(duì)迅速，且存在繼續(xù)增長(zhǎng)的趨勢(shì)，后續(xù)測(cè)點(diǎn)位移變化效果不明顯，且趨于穩(wěn)定，推測(cè)該測(cè)斜孔上部巖體存在個(gè)別不穩(wěn)定塊體。測(cè)斜孔74003內(nèi)測(cè)點(diǎn)位移隨時(shí)間增長(zhǎng)較慢，無明顯規(guī)律且趨于穩(wěn)定，說明暫時(shí)該區(qū)域巖體穩(wěn)定性較好。

為進(jìn)一步確定邊坡潛在滑落面位置，根據(jù)DSR-DDA法，動(dòng)態(tài)折減抗剪強(qiáng)度參數(shù)，對(duì)邊坡漸進(jìn)失穩(wěn)過程進(jìn)行表征，并給出滑落面與相應(yīng)的安全系數(shù)。選取該邊坡最危險(xiǎn)滑落面與次級(jí)滑面驗(yàn)證了方法的可行性，圖9～10為邊坡首次失穩(wěn)與二次失穩(wěn)過程。

第1.360×106時(shí)步前，各塊體監(jiān)測(cè)點(diǎn)位移及邊坡整體基本保持平穩(wěn)狀態(tài)，當(dāng)運(yùn)算至第1.360×106時(shí)步時(shí)，即折減系數(shù)達(dá)到1.136時(shí)，邊坡發(fā)生急劇變形，相應(yīng)塊體位移激增，且后續(xù)變形亦呈逐步增加趨勢(shì)，并且從邊坡整體變形中可以明顯看出下部滑坡體沿相應(yīng)滑落面發(fā)生滑移，說明邊坡發(fā)生首次失穩(wěn)，最危險(xiǎn)滑落面所對(duì)應(yīng)安全系數(shù)為1.136。

邊坡首次失穩(wěn)結(jié)束后，塊體位置及內(nèi)部應(yīng)力重新分布。運(yùn)用首次失穩(wěn)過程結(jié)束后的邊坡模型繼續(xù)運(yùn)行程序，旨在尋找次級(jí)滑落面。

圖9 撫順西露天礦北幫邊坡首次失穩(wěn)過程Fig.9 First failure process of the north slope in West Fushun Open-Pit Mine

圖10 撫順西露天礦北幫邊坡二次失穩(wěn)過程Fig.10 Second failure process of the north slope in West Fushun Open-Pit Mine

第1.890×106時(shí)步前，邊坡整體趨于穩(wěn)定，當(dāng)運(yùn)算至第1.890×106時(shí)步時(shí)，即折減系數(shù)達(dá)到1.189時(shí)，邊坡再次發(fā)生急劇變形，形成次級(jí)滑落面，其對(duì)應(yīng)安全系數(shù)為1.189，且較首次邊坡失穩(wěn)其滑落面范圍更大，而且是在首次失穩(wěn)發(fā)生后邊坡強(qiáng)度略微減小后，與常規(guī)分析得出的結(jié)果大相徑庭(常規(guī)分析中不能剝離首次失穩(wěn)滑坡體，導(dǎo)致其對(duì)二次滑坡體起到保護(hù)作用，從而計(jì)算得到的安全系數(shù)偏高)。究其主要原因是邊坡首次發(fā)生失穩(wěn)后塊體位置及應(yīng)力已重新分布，故應(yīng)使用首次失穩(wěn)結(jié)束后的邊坡模型對(duì)次級(jí)滑落面進(jìn)行分析，而以往研究中并未充分考慮這一因素。

撫順西露天礦邊坡案例破壞模式為牽引式滑動(dòng)破壞，中下部油母頁巖對(duì)邊坡穩(wěn)定性起關(guān)鍵作用，故不可繼續(xù)開采剩余油母頁巖，以免引起邊坡上部大面積滑坡及礦坑-城市邊界滑坡。

3.2 DSR-DDA法與GeoSMA-3D耦合

GeoSMA-3D(巖土工程結(jié)構(gòu)與模型分析系統(tǒng))是一款以Key Block Theory為基礎(chǔ)，在C++框架下編譯源程序，團(tuán)隊(duì)研發(fā)的三維關(guān)鍵塊體可視化分析軟件。該軟件能實(shí)現(xiàn)工程巖體空間結(jié)構(gòu)的建模、結(jié)構(gòu)面的空間表征、模型表面跡線顯示等多項(xiàng)功能，并且其精度已經(jīng)過大量實(shí)例驗(yàn)證[15-17]。

將結(jié)構(gòu)面產(chǎn)狀導(dǎo)入GeoSMA-3D軟件，完成關(guān)鍵塊體搜索過程，并實(shí)現(xiàn)其三維表征;DSR-DDA法中最先發(fā)生滑落的塊體即為關(guān)鍵塊體，關(guān)鍵塊體耦合如圖11所示，且耦合效果良好。

圖11 關(guān)鍵塊體耦合Fig.11 Key block coupling

4 結(jié) 論

(1)提出了受位移閾值控制的DSR-DDA法，給出了DSR-DDA法計(jì)算邊坡多滑面安全系數(shù)的基本過程，為邊坡多滑面穩(wěn)定性分析貢獻(xiàn)了一種新手段。

(2)利用經(jīng)典斜坡滑塊案例驗(yàn)證了DSR-DDA法的可行性與計(jì)算精度，并給出了位移閾值為1mm，對(duì)不同破壞程度下巖體結(jié)構(gòu)面動(dòng)態(tài)折減不同的強(qiáng)度，解決了數(shù)值分析中巖體結(jié)構(gòu)面損傷的程度不一性問題。

(3)基于DSR-DDA法分析了撫順西露天礦北幫邊坡多滑面穩(wěn)定性，得到了其最危險(xiǎn)滑落面及次級(jí)滑落面位置，與常規(guī)分析結(jié)果大相徑庭，佐證了邊坡多滑面穩(wěn)定性分析中模型需動(dòng)態(tài)變化的必要性。該邊坡破壞模式為牽引式滑動(dòng)破壞，中下部油母頁巖對(duì)邊坡穩(wěn)定性起關(guān)鍵作用，故不可繼續(xù)開采剩余油母頁巖，以免引起邊坡上部大面積滑坡及礦坑-城市邊界滑坡。

致謝衷心感謝石根華博士對(duì)筆者的熱忱幫助。