基于精細有限元模型的云岡石窟穩定分析與安全評估

中圖分類號：TU457 文獻標志碼：A 文章編號：2096-6717（2025）03-0066-

竇超1，丁會方1，閆丁2，任建光2，王娟1，3，楊娜1，3（1.北京交通大學土木建筑工程學院，北京100044；2.云岡研究院，山西大同037007；3.北京古遺址保護研究中心，北京100044）

Stability analysis and safety evaluation of Yungang Grottoes based on fine finite elementmodel

DOU Chao1，DING Huifang1，YAN Ding2，REN Jianguang2，WANG Juan ^1，3 YANG Na1.3

（1.School ofCivilEngineering，BeijingJiaotong University，Beijing lOoo44，P.R.China;2.Yungang Grottoes Research Institute，Datong O37Oo7，Shanxi，P.R.China；3.Beijing Ancient SiteProtection Research Center， Beijing 100044，P.R.China）

Abstract： Subjected to various natural forces，the disease of groto temple is serious.The damage degree of groto temple is increasing with time，and the preventive protection for groto temple is becoming more essential. The preventive protection measures mainly consists of the numerical evaluation and on-site monitoring. In order to ensure the accuracy of numerical calculation，this paper proposes the following methodology： First，the three-dimensional laser scanning cloud image of the grottoes temple is used to generate a fine geometric model of the grottoes temple to accurately reflect its complex geometric shape; Secondly，H-B criterion is used to fully consider the deterioration of rock mass mechanical properties caused by weathering， cracks，joints and other diseases，and to obtain reliable rock physical and mechanical parameter，and then establish a refined numerical calculation model that more truly reflects the boundary conditions，stress displacement and possible failure path.Finally，the current stability of the groto temple is analyzed by strength reduction method，and the danger zone of the groto after the deterioration of the rock massis predicted by the longitudinal wave velocity of the rock mass，and an early warning system is established to realize the preventive protection of the grotto temple. Taking Grottoes No.9 and 1O of Yungang Grottoes as an example，systematic modeling，force analysis and stability discussion are caried out by using the above research ideas.The results show that the grottoes maintain a stable state as a whole at present，and stress concentration exists in some areas.By simulating rock strength deterioration，the potential instability risk area of the cave is predicted，and an early warning system based on rock P-wave velocity is established.

Keywords： grotto temple；stability analysis；safety evaluation；H-B criterion； refined numerical model;precautionary analysis

作為一種重要的歷史文化載體，石窟寺具有較高的文化、歷史和觀賞等價值，但在風化、雨蝕以及人類活動等各種營力作用下，各石窟寺受到了不同程度的破壞，其預防性保護研究亟待開展[1-2]。

對于石窟寺的預防性保護研究工作，常見方法有數值分析和現場監測等[3-4]。數值分析作為一種行之有效的研究方法，有著很廣泛的應用。為保證數值分析的準確性，構建的模型需盡量與實際相符，包括幾何信息、力學參數以及邊界條件等。已有不少專家學者對云岡石窟[5-7]、龍游石窟[8-9]、敦煌莫高窟[10-11]、圓覺洞石窟[12-13]等眾多石窟寺開展了數值分析研究，取得了一定成果。但目前對石窟寺的穩定性數值分析研究工作或多或少存在以下問題：1）在建立數值模型時，一般對石窟寺幾何形狀進行了簡化處理，常將復雜的石窟幾何形狀簡化為規則幾何形狀（如文獻[14-15]等）；2）對巖體力學參數的選取過于簡化，未考慮節理分層、裂隙、風化等缺陷等對巖石性能的綜合性影響（如文獻[16-17]等）；3）研究結果針對當時研究的情形，未能建立石窟寺的預防性保護體系（如文獻[18-19]等）。針對上述問題，筆者在盡量符合石窟真實幾何尺寸、現存條件以及應力狀態的基礎上，建立更為精細的幾何和力學數值分析模型。首先，利用石窟寺的掃描云圖生成精細的石窟幾何模型，準確反映石窟內部的復雜幾何形狀；其次，利用H-B準則充分考慮風化、裂隙、節理等病害造成的巖體力學性能劣化，得到接近實際情況的巖石物理力學參數值，進而建立能較為真實反映石窟邊界條件、應力位移及可能破壞情況的精細化數值計算模型；最后，結合強度折減法對石窟的穩定性進行分析，依據巖體的縱波波速對因巖體劣化造成的石窟危險區域進行預測，并建立預警系統，從而實現對石窟寺的預防性保護。

1 云岡石窟概況

作為中國首批國家重點文物保護單位和世界文化遺產，云岡石窟具有極高的歷史文化價值。但由于遭受近1500a自然風化營力的作用和人為損毀，洞窟與石雕受到不同程度的損壞，局部窟頂巖石有崩塌、掉塊的危險。在構造裂隙和風化卸荷裂隙切割的影響下，云岡石窟的第9、10號洞窟巖體的穩定性有所降低，且外立柱風化嚴重，洞窟失穩風險大[20]。為驗證研究思路的可行性，選取具有典型石窟寺病害特點的云岡石窟第9、10號洞窟作為研究對象進行分析。

1.1 洞窟分布

云岡石窟位于山西省大同市西郊云岡鎮，距大同市中心約 16km 。窟區南側、西側毗鄰十里河，北依武周山。云岡石窟現存洞窟53個，大小佛龕1100多個，東西綿延約 1km 。窟區分為東部（1～4窟）中部（5＼～20窟）西部（ 21～53 窟）3部分，為中國規模最大的古代石窟群之一。石窟平面分布見圖1。

云岡石窟第9、10號洞窟位于云岡石窟中部的五華洞，二者同期開鑿，都分前后兩室。窟前都有立柱支撐，其平面布置圖見圖2。

圖2中標號1～7分別為：1為第9、10號洞窟前室隔墻；2為第9、10號洞窟前、后室隔墻；3為第9、10號洞窟后室隔墻；4為第9號洞窟后室佛像；5為第10號洞窟后室佛像；6為第9號洞窟拱門所處位置；7為第10號洞窟拱門所處位置。

第9、10號洞窟的寬度（東西長）、進深（南北長）高度和底板標高見表1。

1.2 石窟地質地貌

云岡石窟開鑿于十里河三級基座階地砂巖中。石窟南臨十里河，窟頂北部為低山丘陵，地勢平緩，呈波狀起伏，地形北高南低，最大高差 55m 左右；石窟東側為沖溝，西側為十里河谷。窟區被南北向大沖溝（東谷、西谷）分為東、中、西3部分。

窟區地層結構較簡單，屬中生代中侏羅統上部云岡組（ 和第四系中上部的中更新統（殘積-沖積）、上更新統（沖積-洪積）全新統（沖積-坡積）20]可將云岡組對石窟有影響的地層劃分為6個巖性段，石窟開鑿于云岡組第 2～4 巖性段。第9、10號洞窟主要開鑿于第4段 ，該段巖性為黃褐色、灰黃色、黃色及肉紅色含礫粗中粒長石巖屑砂巖，粗粒鐵質、鈣質長石巖屑砂巖，中粗粒長石巖屑砂巖及中粗粒鈣質長石巖屑砂巖。巖石質地較疏松，易風化，具有交錯層理（地層分布圖見圖3）。

窟區內巖層產狀平緩，一般傾角為 3^°～4^° ，傾向多為北北東向。斷層和節理均不發育，區內構造以節理為主。區內裂隙以構造裂隙為主，亦發育有風化裂隙。風化裂隙的發育很不規則，多呈不規則的網狀，裂面不平整。第9、10號洞窟所在區段平均節理0.14條/ ^'m² ，該區段節理優勢方向以 81^°～90^° 為主， 41^°～60^° 次之， 351^°～360^° 方向節理也有一定發育，即近東西向張節理（迭加卸荷）在本區段最發育，其次為北東向剪切節理。

2 精細幾何模型的建立

石窟寺數字化保護的發展以及各種建模、計算軟件的交互式應用給石窟寺數值模型的建立帶來了新的思路[13，21]。使用三維激光掃描技術獲取石窟寺的點云圖，利用Geomagicstudio軟件的點云處理功能和Hypermesh的網格劃分功能，建立石窟寺的精細化幾何模型，準確反映洞窟的內部尺寸和形狀細節[17]，精細化模型建立流程圖見圖4。

2.1 點云數據處理

利用三維激光掃描技術掃描石窟內部構造及石窟入口處的外立柱，獲得云岡石窟第9、10號洞窟的點云數據，將三維掃描點云圖導人Geomagic件hm文件，劃分建穩定性分析 立job并提交計算 文件studio軟件（圖5）。初始的點云圖包括543008個點云數據，去除無用的點云數據，并生成三角面片（圖6）。

在圖6的三角面片圖形中，存在對石窟穩定性分析無用的區域，如立柱的護欄，洞窟內部的雜物等，將上述內容去除，并把模型中的孔洞區域進行修補，經過降噪、松弛、構造曲面片等操作生成NURBS曲面，隨后導出為iges文件。

2.2 模型網格處理

將iges文件導人Hypermesh中保存為 hm 文件，并將曲面模型實體化，利用Hypermesh的 3D tetra $$ volumetetra模塊進行網格劃分。由于重點關注對象為洞窟內部窟壁的穩定性，為了減少不必要的計算，在進行網格劃分時，將實體模型分為兩個部分（內部區域尺寸為 31.5m×27.7m×19.8m 總尺寸為 130m×120m×70m） ，內部區域的網格劃分較為精細，非重點關注的外部區域網格劃分較為粗糙[19]，具體為：內部區域平均網格尺寸為 250mm ，在模型網格化后仍能保持真實的輪廓形狀。內部區域單元數量為864497，節點數量為163671，外部區域單元數量為371392，節點數量為69406，采用的單元類型為C3D4，如圖7所示。

為方便對洞窟內部各部位的應力-應變狀態進行描述，對第9、10號洞窟進行剖切，以說明洞窟內部各位置分布情況及各部位名稱，如圖8所示。

2.3 計算范圍的確定

文獻[16]提到，石窟開挖對周圍巖體的應力及位移有明顯影響的范圍大約是開挖輪廓尺寸的2.5～3.0 倍。為了得到更加符合實際的模型尺寸，以第9、10號洞窟開挖輪廓尺寸的倍數建立相應的幾何模型，同時模型地面以上的豎向尺寸（地面至第9、10號洞窟所在邊坡頂部的距離）保持為實際高度。

建立幾何模型之后，在Hypermesh中賦予模型邊界條件。為使模型能在豎向自由沉降變形而不發生側向移動，模型的邊界條件設置為：約束石窟模型底部節點 x，y，z 三個方向的位移， y 方向邊界約束其y方向的位移， x 方向邊界約束其 x 方向的位移。最后，將模型導出為inp文件，再導人ABAQUS中，賦予材性和荷載參數進行計算。

圖9所示為在自重作用下第9、10號洞窟的右3柱上部區域一點的最小主應力值隨不同邊界尺寸的變化曲線。可以看到，該點應力值隨著模型尺寸的增大而增加，但增長速度逐漸減緩。當模型尺寸由洞窟開挖輪廓尺寸的4.0倍增加到7.0倍時，該點應力值的增幅只有 2.8% ，因此，可認為模型外圍尺寸為第9、10號洞窟開挖輪廓尺寸的4.0倍時，便能較好地考慮洞窟邊界范圍的影響。為了減少計算量，取第9、10號洞窟開挖輪廓尺寸的4.0倍建立最終模型。

3 實際破損巖體的力學參數取值

第9、10號洞窟巖體的物理力學參數從《云岡石窟防水保護工程巖土工程勘察報告》20]（以下簡稱《勘察報告》）獲取。出于文物保護的需要，《勘察報告》中獲取洞窟巖體物理參數的巖體樣本并非取自洞窟內部本體，而是取自洞窟所在山體后方的巖層。對后方巖層鉆孔取芯，并對巖體樣本進行室內試驗，得到第9、10號洞窟巖石的物理參數如表2所示。

由于巖體的尺寸效應[22以及洞窟內部及鄰近部注： H 為第9、10號石窟所在邊坡的高度 _;ρ 為巖體密度； σ_c 為巖體單軸抗壓強度； σ_cw 為巖體飽水單軸抗壓強度；u為巖體泊松比： _;E 為巖體彈性模量； φ 為巖體內摩擦角；c為巖體內聚力。

位的巖體存在風化、裂隙、節理等病害，相比新鮮完整的巖石，其力學性能存在一定程度的劣化。而《勘察報告》中的巖體樣本取自較遠區域的深層巖石，不能反映洞窟巖體的真實力學性能，因此，數值分析中一項重要工作就是對表2中的力學參數進行折減，以獲取存在病害和缺陷的石窟巖體的實際力學性能參數[12]。

利用H-B強度準則來考慮巖體病害和缺陷對巖體強度的折減。H-B強度準則是基于對大量巖體強度曲線的系統研究結果提出的，其經驗公式的表達式為[23-24]

σ₁=σ₃+σ_c（m_bσ₃/σ_c+s）^a

式中： σ₁，σ₃ 分別為實際巖體的最大、最小主應力值；σ_c 為完整未破損巖石的單軸抗壓強度； m_b，s 和 α_a 為半經驗參數， m_b 與巖石軟硬程度有關； s 反映巖體破損程度，取值范圍為 0～1 ，完整巖石 s=1;a 與巖體類型有關。

從式（1）看到，H-B強度準則中采用 m_b，s 和 a 三個關鍵參數來反映巖體病害對其強度的折減效果。m_b，s 和 a 可由式（2）～式（4）確定。

式中： D 為巖石擾動參數，取決于外界因素對原位巖體的擾動程度，如爆破、巖體開挖、巖體卸荷等行為，取值范圍為 0～1 ，外界擾動越大，其值越大，可按照式（5）取值[25]。

D=1-K_v

式（5）中 K_v 為巖體完整性系數，可根據巖體中的縱波波速進行取值[26]。

式中： V_ml 為實際巖體中的縱波波速； V_cl 為完整未破損巖石中的縱波波速。 K_v 值越小，表明巖體的完整性越差，節理化程度越嚴重。

根據《勘察報告》可知，第9、10號洞窟的 V_ml=

由式（5）和式（6）可得，第9、10號洞窟的完整性系數 K_v=0.689 ，擾動參數 D=0.311 。

G（GSI） 為巖石的地質強度指標，其值在 0～ 100，取決于巖體的結構、節理等發育程度。實際巖體的地質狀況越好， G 值越大，相應的計算式為[26]

G=1.418B^0.6241-5

B=90+3σ_cw+250K_v

當 σ_cw?90K_v+30 時，以 σ_cw=90K_v+30 代人求 B 值；當 K_v≥0.04σ_cw+0. 4時，以 K_v=0.04σ_cw+ 0.4代人求 B 值。 σ_cw 為巖體飽水單軸抗壓強度； K_v為前述巖體完整性系數。

對于第9、10號洞窟，巖石飽水單軸抗壓強度σ_cw=32.53MPa （表2），則： ：B=359.76，G=50.87 0

m_i 為巖石類型參數，與組成巖石的礦物成分有關，巖體質量越好其值越大。

可查詢表格25獲得 m_i 值，部分取值見表3。根據《勘察報告》，第9、10號洞窟的巖體以多層黃灰色薄層狀細粒～中細粒泥質長石巖屑砂巖為主，因此，取 m_i=17 。

表3 m_i 近似值取值表（部分）

由于巖體的劣化，實際巖體的彈性模量[25也與《勘察報告》中的試驗結果有所差別。可由式（9）、式（10）確定折減后的巖體彈性模量。

當 σ_c?100 時：

當 σ_cgt;100 時：

由表2可知，未破損巖石的單軸抗壓強度 σ_c= 59.42?100 ，則采用式（9），計算得到 E_m=6.844 GPa。

由式（2）～式（4）可得H-B強度準則的半經驗參數 m_b=2.124，s=2.25×10^-3，a=0.505

4有限元分析的M-C準則參數

盡管Hoek-Brown強度準則應用廣泛，能夠有效考慮結構面等缺陷的影響以及巖體結構的非線性破壞特征[2]，但當前H-B本構模型僅集成在少數專用有限元軟件，無法直接在ABAQUS等通用有限元軟件中使用。相比之下，M-C準則因參數簡單、物理意義明確而被廣泛應用于各數值分析軟件中[28]。并且已有研究表明，利用H-B強度準則的等效M-C參數在實際工程應用中效果不錯，如利用等效的M-C參數對某地下工程[29]、公路隧道工程[30]等實際工程問題進行研究，驗證了使用等效M-C參數進行數值分析的可行性。

M-C準則表達式為

式中： σ₁，σ₃ 分別為巖體的最大、最小主應力值； φ 為巖體內摩擦角； c 為巖體內聚力。

將H-B準則的有關參數結果應用于基于ABAQUS軟件的石窟穩定性分析的關鍵在于，如何將H-B強度準則參數等效為M-C強度準則參數φ 和 ∣c∣ ，進而實現在ABAQUS軟件中考慮病害和缺陷影響的石窟巖體的實際力學性能。

為得到等效強度參數，可在如圖10所示的特定最小主應力區間內，基于H-B準則和M-C準則的覆蓋面積相等的原則[28.31]，計算等效M-C參數 φ 和 ∣c∣ 。

等效M-C參數計算式為

式中： σ_3n=σ_3max/σ_c 。

對于邊坡工程，側限應力上限值 σ_3max 由式（14）確定。

式中： γ 為巖體重度； H 為埋深（坡高）； σ_cm 為巖體的整體強度，可由式（15）確定。

由式（12）～式（15）可計算得到等效的M-C強度參數 c=1.015MPa，φ=30.6^°

將計算得到的巖體抗剪強度指標匯總，見表4。

巖土材料破壞的主要形式是剪切破壞和拉伸破壞，對于石窟巖體可能發生破壞的判定，可選擇M-C破壞準則[32]，其表達式為

f_t=σ₃-σ_t

式中： σ₁ 為最大主應力； σ₃ 為最小主應力； c 為巖體內聚力； φ 為巖體內摩擦角； σ_t 為巖體抗拉強度。

當巖體中某一點應力滿足 f_slt;0 ，則發生剪切 破壞。當巖體中某一點應力滿足 f_tgt;0 ，則發生拉伸 破壞。

在ABAQUS軟件中以拉應力為正、壓應力為負，因此，最大主應力云圖一般反映拉應力集中現象，拉應力集中區易發生拉伸破壞；最小主應力云圖一般反映壓應力集中現象，壓應力集中區易發生剪切破壞。可根據石窟應力集中區分布和石窟的塑性區發展情況定性地判斷洞窟可能發生破壞的位置，然后通過判定公式 f_s 和f定量地判斷巖體某點是否發生剪切破壞和拉伸破壞。

對于洞窟的整體穩定性可根據洞窟塑性區的發展情況進行判斷，若洞窟在自重作用下形成了連續的塑性貫通區或滑移面[33-34]，則認為洞窟整體發生了失穩破壞。

5 數值計算結果及穩定性分析

基于得到的M-C巖土力學參數，利用ABAQUS進行數值計算，得到第9、10號洞窟在重力載荷作用下的應力-變形分布云圖，并對其可能的破壞情況進行分析。

5.1石窟塑性區發展

圖11為第9、10號洞窟在重力荷載作用下的塑性區發展情況。可以看到，第9、10號洞窟除了立柱上部和底部、前室中間隔墻下部存在塑性發展區，其他部位沒有塑性區形成。其中，第9號洞窟左側立柱（右2柱）底部等效塑性應變最大。

塑性區呈現零散狀態，未形成塑性貫通區，說明石窟雖然局部出現破損，但在重力作用下整體尚處于穩定狀態。隨著時間推移、巖體強度進一步劣化，上述區域發展為塑性貫通區而發生破壞的風險性較大，在日常巡護中應格外關注。

5.2 石窟最小主應力分布

圖12為石窟在重力荷載作用下的最小主應力分布情況及對應的剖面圖。可以看到，第9、10號洞窟整體應力分布較為均勻（圖12（b）），但在第9號洞窟最右側入口坡腳處（圖12（c））及第10號洞窟最左側入口坡腳處（圖12（d））壓應力集中現象較為突出，此區域易發生剪切破壞，其他部位無明顯應力集中現象。提取壓應力集中區最危險點的應力值，即第9號洞窟最右側入口坡腳處某點的應力值為σ₁=-4.34MPa ， σ₃=-0.58MPa 。由此可得 f_s= ，說明壓應力集中區未發生剪切破壞，目前處于穩定狀態。

5.3石窟最大主應力分布

圖13為石窟在重力荷載作用下的最大主應力分布情況及相應剖面圖。可以看到，大部分區域應力分布較為均勻，但第9、10號洞窟內部以下位置拉應力集中現象較為嚴重：第9號洞窟和第10號洞窟的拱門底部（圖13（c））、9號洞窟前室頂部、9、10號洞窟立柱門梁處及兩窟明窗頂部位置（圖13（d））。由于拉應力集中，這些部位易發生拉伸破壞，造成原有裂隙進一步發展，出現洞窟底部開裂，頂部掉塊等現象，需要引起重視。

5.4 立柱安全系數

對于有立柱承重的石窟寺，立柱的穩定性對整個石窟寺安全十分重要，作為主要承力部位，云岡石窟第9、10號洞窟的外立柱風化嚴重，亟待對其開展安全性評價，立柱的安全系數可通過強度折減法得到。強度折減法的思路為[35]：對巖體的力學參數內聚力 和內摩擦角 φ 的初始值進行人為折減（式（18）），將不同折減系數k值對應的 c_k 和 φ_k 賦予有限元模型進行計算，當石窟某處剛好達到失穩臨界破壞狀態時，對應的折減系數即為安全系數 k_s 。

式中： c_k 為折減后的內聚力； φ_k 為折減后的內摩擦角；k為折減系數。

臨界失穩破壞狀態的判斷依據有3種：數值計算不收斂；特征部位位移曲線出現拐點；形成連續的塑性貫通區[36]。一般而言，以數值計算不收斂為判斷依據會導致結構安全系數偏大，利用塑性貫通區作為判斷依據增加了不確定的人為因素[3]，因此，選擇洞窟立柱特征部位（立柱上部區域，圖11）位移拐點作為立柱安全性評價標準，繪制折減系數-位移（k-u）曲線，見圖14。

根據圖14所示的折減系數-位移曲線可知，石窟立柱的安全系數均在2.5左右，處于15.315（未考慮巖體劣化）和 1.43^[16] （把巖體視為強風化巖體）之間，可認為結果可靠。目前石窟立柱尚處于安全狀態[38]，但隨著風化的加劇，其安全系數會不斷減小，立柱在未來仍存在失穩的可能性。

6考慮巖石強度劣化的石窟危險性預測

對于石窟寺，隨著時間的推移，風化等病害發展會導致巖體整體性和強度降低，可能引起石窟產生穩定性問題。因此，預測石窟發生破壞的位置和時刻，對于石窟安全預警和預防性保護具有重要意義。

由式 （1）～ 式（8）和式（12）～式（15）可以發現，只要確定了石窟的完整性系數 K_v （即巖體的縱波波速 V_ml） ，就可以確定參數 D 和 G （式（5）式（7）式（8））的值，同時 m_i 可由巖石類型確定。從而確定參數 m_b，s 和 ?α （式（2）～式（4）），最后由式 （12）～ 式（15）可以確定M-C等效參數 φ 和 ∣c∣ 。即巖體的縱波波速 V_ml 與巖石力學參數值 φ 和 c 存在一一對應關系，見圖15。而參與數值計算的巖體力學參數 c_k 和φ_k 與折減系數k之間存在對應關系（式（17）），由此可見，折減系數 k 與巖體的縱波波速 V_ml 存在一一對應關系。因此，可基于ABAQUS軟件，通過強度折減法計算不同 c_k 和 φ_k 下石窟的塑性區發展情況，預測未來石窟發生破壞的位置，確定薄弱部位位置，并根據折減系數 k 與巖體的縱波波速 V_ml 的關系建立關于巖體縱波波速的危險性預測方法。

由圖11可知，在當前巖體強度參數數值及重力荷載作用下，第9、10號洞窟的立柱上部和底部、前室中間隔墻下部存在塑性發展區。現利用強度折減法進行數值計算，得到不同折減系數下的洞窟塑性區發展情況，見圖16。

由圖16可以發現，隨著巖體強度的降低，即隨著折減系數的增大，不斷有新的塑性區產生。當折減系數為1.20時，第9、10號洞窟的前室隔墻中部區域及后室隔墻下部區域開始進入塑性（圖16（b））；當折減系數為1.72時，第9、10號洞窟拱門及第9、10號洞窟前、后室隔墻底部區域開始進入塑性（圖16（c））；當折減系數為2.4時，第9號洞窟和第10號洞窟后室佛像部分區域開始進人塑性（16（d））。

由分析可知，隨著巖體的劣化，第9、10號洞窟最先可能破壞的位置為石窟立柱，然后是前室隔墻，后室隔墻，第9、10號洞窟的拱門，第9、10號洞窟前、后室的隔墻、佛像等區域。這些部位應在日常巡檢中引起重視。

根據不同折減系數對應洞窟塑性區的發展情況，并結合折減系數與巖體縱波波速的對應關系（圖15），確定對應的巖體縱波波速的取值范圍，從而建立石窟的縱波波速危險性預測方法，見表5。在石窟的日常維護中可通過對測量潛在失穩區域巖體的縱波波速，來判斷其是否處于危險狀態，提前預測巖體的穩定狀態。

1）構建接近真實的石窟寺幾何模型，綜合考慮各病害對石窟寺巖體物理力學參數的折減，構建的計算模型能較好地應用在石窟寺的預防性保護工作中。

2）通過對云岡石窟第9、10號洞窟分析發現，目前石窟在局部和整體上均處于穩定狀態；對未來因巖體劣化導致的失穩區域進行了預測，建立了以縱波波速為參考的預警系統。

（b）折減系數為1.2（剖切位置1，視圖方向1）

（c）折減系數為1.72（剖切位置2，視圖方向2）

7 結論與討論

7.1 結論

對具有典型病害特征的石窟寺，提出了一種新的預防性保護研究思路，通過建立精細化幾何模型、力學參數模型以及預警系統，以實現對石窟寺的現存穩定性及未來危險區域預測的分析研究，并以云岡石窟第9、10號洞窟為例進行了分析，得到如下主要結論：

7.2 討論

1）研究的優勢與不足

以云岡石窟為例，利用石窟寺的三維激光掃描云圖，并運用H-B準則充分考慮石窟寺的風化、裂隙、節理等病害對巖體力學性能劣化的影響，建立了較為接近實際情況的石窟寺精細化數值計算模型；利用強度折減法對石窟寺的現有穩定性進行分析，結合巖體的縱波波速對因巖體劣化的石窟危險區域進行預測，并建立預警系統，實現對石窟寺的預防性保護，證明石窟寺預防性保護研究思路具備可行性。與之前的研究相比，建立的幾何模型更為精細，巖體物理力學參數取值更符合實際情況。

但是，由于現場勘探數據的不足和文物保護的現實需求，云岡石窟巖體物理力學參數及相關數據獲取有限。基于現有數據，將石窟巖體視為均一整體，采用力學參數折減的思路綜合考慮病害的影響。實際上，石窟巖體在風化營力、構造運動、人類活動等多因素的影響下產生了異質性劣化，不同部位的完整性、風化程度存在差異。因此，算例模型存在一定不足，離反映實際石窟巖體的真實受力狀態還存在差距。

2）未來展望

在后續的工作中需開展詳細的現場勘察工作，并可以結合巖體的物理力學參數試驗和數值反分析方法，減輕勘探作業對石窟寺的損傷，同時獲取盡可能多的真實的巖體數據，在此基礎上對具有不同病害程度的區域分區建模，并賦予不同參數進行分析計算，以提升模型計算的準確性和科學性。

此外，構建的預測模型在后續實際應用中還可以與現場監測相結合，基于現場監測數據以及利用機器學習等算法進行相關病害的綜合性評估預測，建立考慮多因素的綜合性預警系統。

參考文獻

[1］蘭恒星，呂洪濤，包含，等.石窟寺巖體劣化機制與失 穩機理研究進展[J].地球科學，2023，48（4）：1603-1633. LAN H X，LU H T，BAO H，et al.Advances in degradation and instability mechanism of grotto temple rock mass[J].Earth Science，2023，48（4）：1603-1633. （in Chinese）

[2］王金華，陳嘉琦.我國石窟寺保護現狀及發展探析[J]. 東南文化，2018（1）：6-14. WANG J H，CHEN JQ. Current status and future development of cave temples protection in China [J]. Southeast Culture，2018（1）： 6-14.（in Chinese）

[3］安程，呂寧，張榮，等.預防性保護理念對我國石窟寺 保護的影響與實踐[J].東南文化，2020（5）：13-19. AN C，LYU N， ZHANG R，et al. Preventive conservation and its influence to the conservation practice on Chinese grotto temples [J]. Southeast Culture，2O2O（5）：13- 19. （in Chinese）

[4］龔曉南.對巖土工程數值分析的幾點思考[J].巖土力 學，2011，32（2）：321-325. GONG X N. Reflections on numerical analysisof geotechnical engineering [J]. Rock and Soil Mechanics， 2011，32（2）： 321-325.（in Chinese）

[5］晏鄂川，方云.云崗石窟立柱巖體安全性定量評價[J]. 巖石力學與工程學報，2004，23（增刊2）：5046-5049. YAN E C，FANG Y.Quantitative evaluation on the safety of the column rock masses in Yungang Groto [J]. Chinese Journal of Rock Mechanics and Engineering， 2004，23（Sup 2）： 5046-5049.（in Chinese）

[6］李天祺，董松源，周太全.云岡石窟在風化作用下巖體 穩定性的數值模擬分析[J].自然災害學報，2012，21 （6）： 124-130. LITQ，DONGSY，ZHOUTQ.Numerical simulation analysis of stability of Yungang Cavern's rock mass underweathering action [J]. Journal of Natural Disas ters，2012，21（6）：124-130.（in Chinese）

「7］劉仝子彭亞雄旦等云雙窟窟宏圍巖巒形分

析[J].科學技術與工程，2015，15（31）：226-229. LIU JM，PENG Y X，LV Q，et al. Research on stress and deformation of cave chambers in Yungang doublecave [J]. Science Technology and Engineering，2015，15 （31）：226-229. （in Chinese）

［8］高千，王靖，楊志法，等.龍游石窟1～3號洞室穩定性 分析及安全通道路線的選擇[J].巖土力學，2009，30 （9）： 2713-2721. GAO Q，WANG J，YANG Z F，et al.Stability analysis of the large ancient underground rock Caverns in Longyou and the selection of maximum-security routes inthe Caverns [J].Rock and Soil Mechanics，2O09，30 （9）： 2713-2721. （in Chinese）

［9］廖小輝，王雅南，劉浩，等.關于龍游石窟23號古地下 洞室23-1號斜坡柱穩定問題的討論[J].工程地質學報， 2020，28（6）： 1406-1414. LIAO X H，WANG Y N，LIU H，et al. The engineering scientific of the No.23-1 inclined rock pillar in No.23 cavern in the ancient underground rock Caverns inLongyou[J].Journal of Engineering Geology，2020，28（6）： 1406-1414. （in Chinese）

[10]李利莎，杜建國，劉洪麗，等.罕遇地震下敦煌莫高窟 彩塑動力響應數值模擬[J].防災減災工程學報，2017， 37（2）： 201-207. LI L S，DU JG，LIU HL，et al. Numerical simulation of dynamic responses of painted sculptures of Dunuang Mogao Grottoes under rare earthquakes [J]. Journal of Disaster Prevention and Mitigation Engineering，2017， 37（2）：201-207.（in Chinese）

[11]王旭東，張景科，孫滿利，等.莫高窟南區洞窟古人開 鑿智慧初探[J].巖石力學與工程學報，2018，37（增刊 2）： 4207-4220. WANG X D， ZHANG JK， SUN M L， et al. Study on theancient wisdom in the excavation of caves in the south area of Mogao Grottoes [J].Chinese Journal of Rock Mechanics and Engineering，2018，37（Sup 2）： 4207-4220.（in Chinese）

[12]孟志剛.圓覺洞石窟頂板變形機理及控制對策研究 [D].北京：中國礦業大學（北京），2021. MENG Z G. Study on deformation mechanism and control countermeasures of roof in Yuanjue Cave [D]. Bejing：China Universityof Mining amp;.Technology， Beijing，2021. （in Chinese）

[13]劉長青，包含，蘭恒星，等.石窟寺多尺度巖體結構發 育特征與三維精細化建模方法研究——以安岳圓覺洞 為例[J].工程地質學報，2024，32（6）：1904-1915. LIUCQ，BAOH，LANHX，etal.Multi-scalerock structures development characteristicsand 3D refinement modeling method of grottoes： A case study of Anyue Yuanjue Cave [J]. Journal of Engineering Geology， 2024，32（6）：1904-1915. （in Chinese）

[14］郭志謙.敦煌莫高窟南區密集洞窟群穩定性及危巖體 風險評估[D].蘭州：蘭州大學，2018. GUOZ Q.Stability of dense caves in the south area of Dunhuang mogao grottoes and risk assessment of dangerous rock mass [D]. Lanzhou： Lanzhou University， 2018. （in Chinese）

[15]盧帥明.云岡石窟9、10號窟柱體穩定性分析[D].沈陽： 沈陽建筑大學，2018. LU S M. Stability analysis of Pillars in caves 9 and 10 of Yungang Grottoes [D]. Shenyang： Shenyang Jianzhu University，2018.（in Chinese）

[16]劉積魁.云岡石窟9、10窟三維穩定性分析與地震動力 響應模擬[D].武漢：中國地質大學，2011. LIUJ K. Three-dimensional stability analysisand seismic dynamic response simulation of caves 9 and 10 in Yungang Grottoes [D]. Wuhan： China University of Geosciences，2011.（in Chinese）

[17]黃皇.圓覺洞紋層狀砂巖圍巖穩定性分析[D].北京：中 國地質大學（北京），2021. HUANG H. Stability analysis of surrounding rock of layered sandstone with circular cave pattern [D]. Beijing： ChinaUniversity of Geosciences，2O2l.（in Chinese）

[18］馬雨君.基于精細化勘察的臥佛佛手及目犍連造像穩 定性研究[D].蘭州：蘭州大學，2022. MA Y J. Study on the stability of the statues of Buddha's finger grown and sleeping Buddha’s finger grown based on fine investigation [D].Lanzhou：Lanzhou University， 2022. （in Chinese）

[19］黃蕾.大足石刻臥佛造像區局部巖體穩定性評價[D]. 蘭州：蘭州大學，2023. HUANG L. Stability evaluation of local rock mass in the reclining Buddha statue area of Dazu stone carving [D]. Lanzhou：Lanzhou University，2O23.（in Chinese）

[20]云岡石窟防水保護工程巖土工程勘察報告[R].建設綜 合勘察研究設計院，2003. Geotechnical investigation report of Yungang Grottoes waterproofing protection project [R]. Construction of Comprehensive Survey，Research and Design Institute， 2003.（in Chinese）

[21]李敏，刁常宇，葛云飛，等.石窟寺文物的數字化保護 與利用[J].遙感學報，2021，25（12）：2351-2364. LIM，DIAOCY，GEYF，et al.Digital protection and utilization of Grotto cultural relics [J].National Remote Sensing Bulletin，2021，25（12）： 2351-2364. （in Chinese）

[22]蘭恒星，包含，孫巍鋒，等.巖體多尺度異質性及其力 學行為[J].工程地質學報，2022，30（1）：37-52. LAN H X，BAO H，SUN W F，et al.Multi-scale heterogeneity of rock mass and its mechanical behavior [J]. Journal of Engineering Geology，2022，30（1）： 37-52. （in Chinese）

[23]HOEK E，BROWN E T. Practical estimates of rock mass strength [J]. International Journal of Rock Mechanics and Mining Sciences，1997，34（8）： 1165-1186.

[24]HOEK E，BROWN E T.The Hoek-Brown failure criterion and GSI-2O18 edition [J]. Journal of Rock Mechanics and Geotechnical Engineering，2Ol9，11（3）： 445-463.

[25]陳鵬.基于Hoek-Brown準則的邊坡等效MohrCoulomb參數估算[D].西安：長安大學，2021. CHEN P. Estimation of equivalent Mohr-Coulomb parameters of slope based on Hoek-Brown criterion [D]. Xi'an：Chang'an University，2O21. （in Chinese）

[26]工程巖體分級標準：GB/T 50218—2014[S].北京：中 國計劃出版社，2015. Standard for engineering classification of rock mass： GB/ T 50218—2O14 [S]. Beijing：China Planning Press， 2015.（in Chinese）

[27］金俊超，余成學，尚朋陽.基于Hoek-Brown準則的應 變軟化模型有限元數值實現研究[J].工程力學，2020， 37（1）： 43-52. JINJC，SHEC X，SHANG PY.A finite element implementation of the strain-softening model based on thehoek-brown criterion [J]. Engineering Mechanics， 2020，37（1）： 43-52.（in Chinese）

[28]吳順川，張敏，張詩淮，等.修正Hoek-Brown準則的等 效Mohr-Coulomb強度參數確定方法研究[J].巖土力 學，2019，40（11）： 4165-4177. WU S C， ZHANG M， ZHANG S H，et al. Study on determination method of equivalent Mohr-Coulomb strength parameters of a modified Hoek-Brown failure criterion [J]. Rock and Soil Mechanics，2019，40（11）： 4165-4177.（in Chinese）

[29］田茂霖.HOEK-BROWN強度準則的參數取值研究及 其工程應用[D].山東 青島：山東科技大學，2017. TIAN ML. Study on parameter selection of HOEKBROWN strength criterion and its engineering application [D].Qingdao，Shandong： Shandong University of Science and Technology，2Ol7. （in Chinese）

[30]胡其志，劉倩，丁志剛，等.基于Hoek-Brown準則的強 度等效方法對比及應用研究[J].隧道建設（中英文）， 2021，41（10）： 1662-1671. HU Q Z，LIU Q，DING ZG，et al. Comparison and application of strength equivalent methods based on Hoek-Brown criterion [J]. Tunnel Construction，2021， 41（10）： 1662-1671. （in Chinese）

[31] SONG Y H，FENG M，CHEN P.Modified minimum principal stress estimation formula based on Hoek-Brown criterion and equivalent Mohr-Coulomb strength parame ters[J]. Scientific Reports，2023，13： 6409.

[32]劉煜，崔德山，楊赫楠，等.須彌山石窟5號窟窟檐的 地震動力響應模擬[J].安全與環境工程，2022，29（4）： 101-109. LIUY，CUI DS，YANGHN，et al.Seismic dynamic response simulation ofthe eaves of cave No.5of Xumishan Grottoes [J]. Safety and Environmental Engineering，2022，29（4）：101-109.（inChinese）

[33]候玲，薛海斌，周澤華，等.基于ABAQUS二次開發 平臺的邊坡有限元強度折減法研究[J].西安理工大學 學報，2016，32（4）：449-454，499. HOUL，XUEHB，ZHOUZH，etal.Researchon the strength reduction FEM of the slope by applying the secondary development platform of ABAQUS[J]. Journal ofXi'anUniversity of Technology，2Ol6，32（4）： 449-454，499.（inChinese）

[34]尹小濤，薛海斌，湯華，等.邊坡局部和整體穩定性評 價方法的辯證統一[J].巖土力學，2018，39（Supl）： 98-104. YIN X T，XUE HB，TANG H，et al. Dialectical unity ofslope local and global stability analysis methods[J]. Rockand Soil Mechanics，2018，39（Supl）：98-104.（in Chinese）

[35]HUACY，YAOLH，SONGCG，etal.Variational method for determining slope instability based on the strength reduction method [J]. Bulletin of Engineering Geology and the Environment，2022，8l（1O）： 395.

[36]張鵬，張紹和，劉磊磊.基于物質點強度折減法的二維 均質 c-φ 邊坡穩定性分析圖表[J].工程地質學報， 2023，31（2）：596-606. ZHANGP，ZHANGSH，LIULL.Two-dimensional homogeneous c-φ slope stability analysis charts based on material point strength reduction method [J]. Journal of EngineeringGeology，2023， 31（2）： 596-606．（in Chinese）

[37]商兆濤，姚家李，夏琴，等.基于強度折減理論的深基 坑穩定性分析[J]合肥工業大學學報（自然科學版）， 2022，45（9）：1209-1215. SHANG Z T，YAO JL，XIA Q，et al.Stability analysis of deep foundation pit based on strength reduction method [J]. Journal of Hefei Universityof Technology （Natural Science），2022，45（9）：1209-1215. （in Chinese）

[38]建筑地基基礎設計規范：GB50007—2011[S].北京：中 國計劃出版社，2012. Code for design of building foundation：GB 5ooo7—2011 [S].Beijing：China Planning Press，2Ol2.（in Chinese）

（編輯胡英奎）