吳川 廖坤陽 孫大齊

摘要：為了從時間尺度和空間變異性上研究石質(zhì)文物危巖體的結(jié)構(gòu)穩(wěn)定性，建立了巖體強度參數(shù)的空間分布和時間退化模型。選取麥積山石窟13號窟危巖體為研究對象，利用強度折減法對石窟巖體穩(wěn)定性進行模擬計算。模擬結(jié)果表明：① 麥積山13號窟現(xiàn)狀條件下最大位移為1.8 mm，中部位移方向略微傾向石窟外側(cè)，石窟頂部呈現(xiàn)局部塑性變形，但整體穩(wěn)定性較好，模擬結(jié)果與現(xiàn)狀調(diào)查較為一致。② 加入時間衰變模型后發(fā)現(xiàn)，13號窟現(xiàn)狀條件下頂部位移為1.2 mm;預(yù)測50 a后洞窟最大位移位于洞頂，為1.3 mm;100 a后洞窟最大位移從洞頂轉(zhuǎn)移至坡腳，為16.7 mm。③ 強度折減分析表明，100 a后模型產(chǎn)生大面積塑性區(qū)域貫通現(xiàn)象，洞頂發(fā)生拉張破壞，洞底剪切破壞。時空變異模型有效提高了傳統(tǒng)靜力學計算精度，通過時間退化模型真實地再現(xiàn)了石窟危巖體的動態(tài)破壞過程，在精確評價危巖體長期穩(wěn)定性方面具有廣泛應(yīng)用價值，可為危巖治理提供參考。

關(guān) 鍵 詞：

麥積山石窟; 時空變異性; 非平穩(wěn)隨機過程譜表示法; 危巖體穩(wěn)定性; 石質(zhì)文物

中圖法分類號： TU47

文獻標志碼： A

DOI：10.16232/j.cnki.1001-4179.2021.06.013

危巖體的穩(wěn)定性評價和分析一直是巖土工程、地質(zhì)工程領(lǐng)域的重要研究方向[1]。目前危巖體穩(wěn)定性分析方法主要分為兩類：第一類以剛體極限平衡法為基礎(chǔ)，如陳洪凱等[2]將三峽庫區(qū)危巖體失穩(wěn)破壞分為滑塌式、傾倒式和墜落式，并在考慮自重、地震和裂隙水壓力的情況下給出了3種危巖體的靜力學計算方法。第二類評價方法以數(shù)值方法（如限元法、有限差分和離散元）為基礎(chǔ)，如鄭穎人等[3]將有限元強度折減法應(yīng)用于巖土邊坡穩(wěn)定性評價中，開創(chuàng)了求解巖質(zhì)邊坡滑動面與穩(wěn)定安全系數(shù)的先河。前一類方法過于簡便。第二類評價法憑借計算機強大的運算效率劃分單元體進行穩(wěn)定性評價，卻對于巖體性質(zhì)的變異性考慮不足，如巖體參數(shù)（凝聚力、內(nèi)摩擦角、重度等）在空間分布的不均勻性和時間上是不斷變化的。Symbol`@@石質(zhì)文物由于遭受長年累月的風化侵蝕作用，整體穩(wěn)定性問題是文物保護工作的重點，同時由于其文物屬性無法采取大量破壞性的工程防護措施，因此其在時間和空間尺度上的穩(wěn)定性研究具有現(xiàn)實意義。

目前，考慮時空變異性的石質(zhì)文物危巖穩(wěn)定性評價的成果較少，大部分研究成果仍處于現(xiàn)象描述階段[4-7]。本文以麥積山石窟13號窟為研究對象，綜合以上研究成果，利用譜表示法所生成的隨機場來模擬巖體抗剪強度參數(shù)空間變化特征。作為一種嘗試，建立基于強度衰減規(guī)律的時間變化模型。借助有限元強度折減法動態(tài)分析了石窟危巖體穩(wěn)定性，并針對性地提出了適用于石質(zhì)文物保護的防治方案。

1 空間參數(shù)分布模型

各種成因的環(huán)境地質(zhì)作用會導致巖體在空間分布上出現(xiàn)差異性[8]，如發(fā)育于巖體中結(jié)構(gòu)面的幾何分布特征，對于巖體強度、變形和滲透性的影響不容忽視。由于巖體強度參數(shù)既具有變異性又具有相關(guān)性，無法用簡單線性函數(shù)表示，所以選擇非平穩(wěn)隨機過程譜表示法進行模擬。譜表示法于1990年被Shinozuka等用于模擬生成隨機場。本文采用了梁建文[9]推導的非平穩(wěn)隨機過程蒙特卡洛模擬的一個譜表示法，該方法產(chǎn)生的樣本函數(shù)能夠較準確反映巖體參數(shù)在水平方向的延續(xù)性與垂直方向上的變異性。

巖體強度參數(shù)是一個切實存在的變量（不能小于零）。依據(jù)摩爾-庫倫剪切強度準則可知，對于巖體穩(wěn)定性影響最為顯著的參數(shù)是凝聚力和內(nèi)摩擦角，而其他參數(shù)如容重、彈性模量以及摩擦角等空間變異性不明顯，對于巖體穩(wěn)定性分析貢獻很小?？紤]到數(shù)值計算效率，本文假定其他參數(shù)恒為常量，選用二維正態(tài)分布的非平穩(wěn)隨機場來表示基本符合模擬要求。二維對數(shù)正態(tài)分布表達式如下：

fx，y=exp{ln1+V2×

M-1i=0N-1j=0

σij[Vijθcosω1ix+ω2jy+

Wijθsinω1ix+ω2jy]+

lnμ-ln1+V22}（1）

式中：x和y為單元體空間坐標;V，μ為別為凝聚力和內(nèi)摩擦角的標準差和平均值;σij為（iM+j+1）項的標準差;Vij（θ），Wij（θ）分別為相對獨立的標準正態(tài)分布隨機變量;ω1i，ω2j則分別為頻率坐標。

其中，σij的計算公式如下：

σij=4Sω1i，ω2jΔω1Δω2（2）

式中：S為單邊功率譜密度函數(shù);Δω1，Δω2分別為頻率坐標軸的離散區(qū)間。

由于相關(guān)函數(shù)對于巖體結(jié)構(gòu)穩(wěn)定性指標影響不明顯[9]，對石質(zhì)文物穩(wěn)定性模擬中選擇計算量較小的二維指數(shù)函數(shù)作為相關(guān)函數(shù)：

ρx，y=exp-2xδx+2yδy（3）

式中：ρ代表自相關(guān)系數(shù);x，y代表水平和垂直方向距離，m;δx和δy分別代表水平方向和垂直方向的波動范圍，m。

現(xiàn)場測定樣本越多，就越能模擬出與實際情況相近的空間分布情況，理論上此種模擬方法可以獲得無限接近于真實情況的解，但在實際中是不可能的。通過譜表示法，根據(jù)有限個關(guān)鍵點出測得的強度參數(shù)（c，φ）統(tǒng)計值進行隨機場模擬，從而反演出可能的巖土體強度參數(shù)空間分布情況是可行的。由于文章篇幅限制，具體過程可參見文獻[8]。

2 時間衰變模型

描述巖石力學特性隨時間變化的規(guī)律是危巖體穩(wěn)定性分析中的難題，從時間尺度上評價巖體的穩(wěn)定性研究尤為重要。許宏發(fā)等[10]通過實驗定量研究了軟巖的長期強度特征，認為巖體存在一組初始強度和最終強度，隨時間推移初始強度逐漸降低并無限接近于最終強度，導致這一變化過程是巖體內(nèi)壁微觀損傷積累的結(jié)果。本文根據(jù)李連崇[11]關(guān)于巖體細微損傷隨時間變化試驗研究，嘗試建立了以指數(shù)衰減的力學退化模型（見圖1），用以描述危巖體的長期強度。

從圖1的退化模型中可知，巖體的強度（σt）隨時間逐漸衰退，初始強度（σ0）和長期強度（σl）固定不變;當t=0時，σt=σ0，當t→∞時，σt=σl。

由此將巖體微觀表征單元強度的時間退化模型用式（4）表示：

σt=σl+σ0-σle-kt（4）

式中：系數(shù)k和長期強度σl與巖體賦存環(huán)境有關(guān)。k反映強度衰減的快慢程度，若巖體所處環(huán)境有風化剝蝕、地震等劣化巖體結(jié)構(gòu)的因素存在，那么系數(shù)k勢必會較大，同時σl的大小受環(huán)境條件與巖性條件共同控制。對于以上兩種量值的選取需根據(jù)現(xiàn)場和室內(nèi)相關(guān)標定試驗來進行。

3 基于時空變異模型的評價方法

在將時間和空間變異模型加載之前，用強度折減法自動搜索滑動面，且應(yīng)用數(shù)值模擬軟件對其穩(wěn)定性進行分析[12-13]。本文將新鮮巖體和強風化巖體的強度指標作為初始強度σ0和長期強度σl。對于空間分布參數(shù)隨機場的模擬，選取麥積山石窟關(guān)鍵點的強度指標輸入譜表示計算法中，根據(jù)關(guān)鍵點強度參數(shù)分布規(guī)律，譜表示法反演出擬合最好的巖體整體空間參數(shù)分布情況。

（1） 基礎(chǔ)模型建立。首先在ANSYS中對模型進行剖分。通過程序轉(zhuǎn)換將建立好的三維模型導入FLAC3D。再在FLAC3D中對模型區(qū)域、邊界條件、材料屬性和模型等進行詳細設(shè)置。最后將完整的模型文件保存為后綴為xml的接口文件。

（2） 空間參數(shù)設(shè)定。根據(jù)所推演出的譜表示法模擬出不同方向（水平和垂直）波動范圍的隨機場，得到c和φ的n組隨機組合。然后將模型文件中的強度參數(shù)用以上n組數(shù)據(jù)替代，從而得到基于空間變異的新模型。

（3） 時間參數(shù)設(shè)定。根據(jù)強度隨時間衰變模型，將上一步所推導出的n組強度參數(shù)根據(jù)衰變公式折減后得n組新數(shù)據(jù)，基于這n組新數(shù)據(jù)的模型即為t時刻的危巖體三維模型，其具有時間和空間上的變異性。

（4） 穩(wěn)定性計算。利用WinbatchTM（批處理軟件）賦值，在FLAC3D中對新模型文件的穩(wěn)定性進行動態(tài)分析，得出計算結(jié)果。

4 算例分析

4.1 模型尺寸與網(wǎng)格劃分

采用三維激光掃描建立麥積山石窟的三維地形曲面。為了確保建模過程中模型不受網(wǎng)格離散化和建模范圍的影響，本文選取的模型范圍超過研究洞窟尺寸的2.5～3.0倍。依據(jù)三維激光掃描確定的模型曲面較為復(fù)雜和貼合實際情況，本文通過ANASYS采用四面體單元進行建模處理（見圖2），模型建好后導入到FLAC3D中進行數(shù)值分析。

由于麥積山石窟13號窟為南北向布置，為了便于計算，本文的模型坐標系設(shè)定豎直方向為Z軸，水平方向為X和Y軸，其中X軸正方向朝南，Y軸正方向朝東，石窟邊坡巖體走向與Y方向平行。最終確定模型尺寸為垂直方向（Z）高50 m，南北方向（X）長75 m，東西向（Y）長35 m。另外，考慮到石窟巖體表面性狀復(fù)雜，可能對后續(xù)的模型計算效率產(chǎn)生影響，因此本文對石窟中造像進行了一定程度簡化，并不影響石窟整體穩(wěn)定性計算結(jié)果，簡化后的模型見圖2。

4.2 邊界條件與模型參數(shù)

模型計算過程中設(shè)定頂面和底面為自由邊界，模型南側(cè)（前側(cè)）同樣為自由邊界，側(cè)面和后面為固定約束。破壞準則參考摩爾-庫侖破壞準則?，F(xiàn)場選取石窟基巖砂巖進行室內(nèi)物理力學試驗。由于研究對象為石質(zhì)文物，不宜進行大規(guī)模破壞性勘探，本文利用輕型便攜式鉆機沿石窟頂面進行鉆探，共布設(shè)鉆探剖面6條，鉆孔數(shù)量24個，鉆孔孔深2.0～5.0 m，通過對巖芯進行力學試驗發(fā)現(xiàn)弱風化巖體飽和單軸抗壓強度一般為30.0～40.0 MPa，強風化巖體一般為12.6～18.7 MPa，相較而言強風化巖體強度較弱風化強度衰減約38.1%～47.5%?，F(xiàn)場調(diào)查發(fā)現(xiàn)石窟巖體以弱風化和強風化巖體為主。模型建模所用巖體力學試驗結(jié)果見表1。

石窟模型在空間上的力學參數(shù)通過譜表示法來表征，根據(jù)有限個關(guān)鍵點測出的抗剪強度參數(shù)統(tǒng)計值進行隨機場模擬。由于巖體內(nèi)變形模量、泊松比和密度空間變化不大，因此將其簡化為常量。本文通過現(xiàn)場鉆探取樣過程中確定的關(guān)鍵點的物理力學特性指標，將其換算成有關(guān)二維正態(tài)分布的函數(shù)值，其計算參數(shù)見表2。

4.3 模擬成果分析

4.3.1 現(xiàn)狀位移分布規(guī)律

將空間分布模型導入麥積山石窟模型力學參數(shù)中，用FLAC3D進行麥積山石窟的位移和應(yīng)力分布規(guī)律分析，其模型的總位移矢量和等值線云圖見圖3。

由圖3可知，在自重應(yīng)力作用下，石窟總體的位移呈現(xiàn)從頂部到底部均勻降低的規(guī)律：在山頂部最大位移為1.8 mm，中部位移方向略微傾向石窟外側(cè)，底部位移基本為0，說明該模型自重應(yīng)力分布規(guī)律基本合理。

4.3.2 時間尺度預(yù)測

在石窟整體位移分析的基礎(chǔ)上，為了更全面分析和研究洞窟的時間衰變規(guī)律，對洞窟中間位置切割一處典型剖面，依據(jù)時間衰變模型設(shè)定強度衰減系數(shù)k為0.02，預(yù)測了13號窟隨時間的應(yīng)力及位移變化規(guī)律。

圖4（a）為t=0現(xiàn)狀條件下的石窟中軸線剖面位移云圖，可以很明顯發(fā)現(xiàn)剖面的總位移方向向模型下后方傾斜，并且最大位移位于石窟頂部，其值為1.8 mm，平臺底部位移為0。13號洞窟頂部位移為1.2 mm，底部位移為0.6 mm，位移方向均為下后方。

圖4（b）代表t=50 a即50 a后的石窟剖面位移云圖。對比發(fā)現(xiàn)模型總位移呈增大趨勢，最大值0.3 mm，同時平臺底部位移依舊為0，說明該段時間石窟穩(wěn)定性依舊較好。

當t=100 a時，模型頂部位移幾乎無變化，最大位移從石窟頂部轉(zhuǎn)移到底部平臺，石窟山腳的最大位移達到16.7 mm，13號洞窟頂板的沉降同樣達到13.7 mm，說明石窟已經(jīng)發(fā)生了不同程度的局部失穩(wěn)破壞。現(xiàn)場調(diào)查情況發(fā)現(xiàn)，13號洞窟砂巖巖體受風化和裂隙滲水侵蝕較為嚴重，已經(jīng)存在不同大小的風化凹槽，若不采取措施，可以預(yù)見100 a后石窟必然會發(fā)生一定程度的穩(wěn)定性破壞問題。

模擬發(fā)現(xiàn)50 a內(nèi)石窟山頂和山腳處應(yīng)力集中尚不存在，但是隨著時間推移應(yīng)力逐漸在兩個位置集中。在100 a后模型中山體向下應(yīng)力逐漸增大，洞窟頂板和前壁洞口均發(fā)生了應(yīng)力集中現(xiàn)象。

4.3.3 穩(wěn)定性時空評價

加入空間參數(shù)分布模型和時間衰變模型后，強度拆減理論計算量大大超過傳統(tǒng)計算方式。因此本文將折減系數(shù)限定在1.0～3.0進行試運算。研究發(fā)現(xiàn)，當t=0，R=1.42時，模型頂部的局部區(qū)域開始發(fā)生塑性破壞，說明現(xiàn)狀條件下13號窟整體穩(wěn)定性尚可，但是石窟頂部可能發(fā)生部分破壞。當t=100 a，R=1.10時，模型開始發(fā)生大面積塑性區(qū)域貫通現(xiàn)象。如圖5所示，在洞口上部黃色區(qū)域為大面積的拉張破壞，在洞窟底部紅色區(qū)域為剪切破壞，因此可以預(yù)測100 a石窟整體處于臨界破壞，隨時可能發(fā)生大規(guī)模的頂部巖體崩塌現(xiàn)象。

對于石質(zhì)文物巖體的穩(wěn)定安全性系數(shù)，目前并沒有統(tǒng)一的認識。參照現(xiàn)階段工程巖體邊坡相關(guān)規(guī)范，考慮到文物的價值屬性，石窟屬于重要保護文物，防治工程等級按一級考慮，危巖體防治安全系數(shù)應(yīng)為1.40～1.60。目前來看13號窟整體處于穩(wěn)定狀態(tài)。

5 結(jié) 論

（1） 巖體由于各種成因地質(zhì)作用，其強度參數(shù)（c，φ）在空間分布上具有不確定性，而空間分布參數(shù)模型的引入克服了傳統(tǒng)靜力學計算中這一缺點，使得計算更加精確和符合實際情況。石質(zhì)文物巖體由于巖體流變特性，其強度參數(shù)隨時間會逐漸降低，導致強度降低的主要因素是軟弱結(jié)構(gòu)面受風化剝蝕等內(nèi)外地質(zhì)作用。時間衰變模型基于實際勘測成果從整體上把握了巖體強度衰變規(guī)律，預(yù)測了麥積山石窟13號窟長時間的穩(wěn)定性，取得良好效果。

（2） 麥積山石窟現(xiàn)狀條件下最大位移為1.8 mm，石窟頂部呈現(xiàn)局部塑性變形，但整體穩(wěn)定性較好，模擬結(jié)果與現(xiàn)狀調(diào)查較為一致;加入時間衰變模型后發(fā)現(xiàn)洞窟頂部位移為1.2 mm，預(yù)測50 a后洞窟最大位移位于洞頂，為1.3 mm，100 a后洞窟最大位移從洞頂轉(zhuǎn)移至坡腳，為16.7 mm;強度折減分析表明100 a后模型產(chǎn)生大面積塑性區(qū)域貫通現(xiàn)象，洞頂發(fā)生拉張破壞，洞底發(fā)生剪切破壞。

（3） 目前，將強度參數(shù)的時間和空間變異性結(jié)合起來評價巖體穩(wěn)定性的研究尚處于初始階段。本文只是考慮了凝聚力和摩擦角的不確定性，對于重度、變形模量、邊界條件以及理論模型等假定其恒定不變，因此對于重要程度較高的復(fù)雜巖體，今后的穩(wěn)定性評價應(yīng)朝更精細、更復(fù)雜的方向發(fā)展。

參考文獻：

[1] 唐輝明.工程地質(zhì)學基礎(chǔ)[M]北京：化學工業(yè)出版社，2008：15-21.

[2] 陳洪凱，唐紅梅，王蓉.三峽庫區(qū)危巖穩(wěn)定性計算方法及應(yīng)用[J].巖石力學與工程學報，2004，23（4）：614-619.

[3] 鄭穎人，趙尚毅.有限元強度折減法在土坡與巖坡中的應(yīng)用[J].巖石力學與工程學報，2004，23（19）：3381-3388.

[4] 姜效璽，吳瓊，鄧小龍，等.復(fù)雜形態(tài)立柱式石質(zhì)文物的病害調(diào)查及防治對策：以浙江省柯巖云骨巖柱為例[J].工程勘察，2016，44（10）：1-5.

[5] 甄剛，秦立科，馬宏林.乾陵無字碑裂隙分布情況調(diào)查及其穩(wěn)定性分析[J].文博，2015（1）：94-98.

[6] 唐智亮.固結(jié)灌漿在保護云崗石窟石質(zhì)文物遺址的應(yīng)用研究[D].長春：吉林大學，2013.

[7] 田小甫.太原晉陽大佛邊坡巖體風化分級及地震穩(wěn)定性評價研究[D].北京：中國地質(zhì)大學（北京），2009.

[8] 潘別桐.巖體結(jié)構(gòu)面網(wǎng)絡(luò)模擬及應(yīng)用[M].武漢：中國地質(zhì)大學出版社，1987：1-2.

[9] 梁建文，肖笛.非平穩(wěn)隨機過程模擬的一個譜表示方法[J].應(yīng)用概率統(tǒng)計，2005，21（4）：375-386.

[10] 許宏發(fā).軟巖強度和彈模的時間效應(yīng)研究[J].巖石力學與工程學報，1997，16（3）：246-246.

[11] 李連崇，徐濤，唐春安，等.單軸壓縮下巖石蠕變失穩(wěn)破壞過程數(shù)值模擬[J].巖土力學，2007，28（9）：1978-1982，1986.

[12] 趙尚毅，鄭穎人，時衛(wèi)民，王敬林.用有限元強度折減法求邊坡穩(wěn)定安全系數(shù)[J].巖土工程學報，2002，24（3）：343-346.

[13] 向娟，陳占鋒，范文臣.考慮巖體節(jié)理密度影響的高陡邊坡穩(wěn)定性研究[J].中國科技論文，2018，13（24）：2790-2795.

（編輯：鄭 毅）

Stability analysis on dangerous rockmass based on spectral representation

method and time degradation model

WU Chuan1，LIAO Kunyang2，SUN Daqi3

（1.Chongqing Industry Polytechnic College，Chongqing 401120，China; 2.Fujian Forestry Vocational & Technical College，Longyan 353000，China; 3.China Academy of Railway Sciences，Beijing 100081，China;）

Abstract：

In order to study the structural stability of stone cultural relic dangerous rock masses from time scale and spatial variability，a model of spatial distribution and time degradation of rock mass strength parameters was established.The dangerous rock mass of Grotto 13 in the Maiji Mountain Grottoes was selected as the research object.The strength reduction method was used to simulate and calculate the rock mass stability.The simulation results show that the maximum displacement under the current conditions of the Maiji Mountain Grotto 13 is 1.8 mm，displacement direction of the middle part is slightly inclined to the outside of the cave，and the top of the cave shows local plastic deformation，however the overall stability is good，and the simulation results are consistent with the current survey.By using the time degradation model，it is found that the top displacement of Grotto 13 is 1.2 mm under current conditions，and the maximum displacement of the cave is predicted to be 1.3 mm after 50 years，still locating at the cave top，however the maximum displacement of the cave will be transferred from the cave top to the slope foot after 100 years，reaching 16.7mm.It shows that the model will have a large plastic area penetration phenomenon after 100 years，top of the cave shows tensile failure，and the bottom of the cave shows shear failure.The spatio-temporal variation model effectively improves the accuracy of traditional static mechanicalcal culations，and the time degradation model truly reproduces the dynamic failure process of dangerous rock masses in caves.It has extensive application values in accurately evaluating the long-term structural stability of dangerous rock masses，and can be used as references for dangerous rock treatment project.

Key words：

Maiji Mountain Grottoes;spatio-temporal variability;spectral representation method of nonstationary random process;stability of dangerous rock masses;stone cultural relics