基于雙邊濾波石質文物逆向建模及其穩定性分析

關鍵詞：大足石刻；石質文物；雙邊濾波；強度折減法；Hoek-Brown準則；穩定性評價

中國現存石窟寺及摩崖造像共5986處，因其蘊含的歷史、藝術、科學價值受到廣泛關注。不可移動石質文物在千百年自然環境與人為因素影響下，其本體結構不斷劣化，出現變形、開裂，甚至傾倒等病害，對文物本體穩定性產生不利影響，嚴重損害了文物價值保存。因此，針對石窟石刻的穩定性研究，一直是中國文化遺產保護工作的重點之一。

目前，圍繞石質文物穩定性的研究尚未形成完整的學科體系，缺少規范概念；眾多學者借鑒工程地質學、巖石力學等學科知識，將巖土體穩定性分析方法及理論應用于文物研究。隨著計算機的發展，數值模擬分析方法成為石質文物穩定性評價的主要方法。Ambraseys等[1]通過離散元法（discreteelementmethod，DEM）對阿波羅雕像穩定性做出評估，并將分析結果與地震發生后雕像旋轉情況進行對比，以驗證分析方法的合理性。Guo等[2?3]和Peng等[4]結合現場調查信息，利用FLAC3D分別研究了敦煌莫高窟南區洞窟和云岡石窟在靜力及動力狀態下的穩定性分析。Chen等[5]為克服大地測量法在文物研究的局限性，基于差分層析SAR（DTomoSAR）和有限元法（finiteelementmethod，FEM），通過有效監測吳哥窟的變形，實現對潛在不穩定結構發生變形破壞的預測。已有研究在對石質文物進行穩定性分析的過程中，以契合研究區域幾何特點為主。隨著三維激光掃描技術在文物保護領域的應用，大量文物的點云信息被采集，鄧小龍[6]和高相波[7]提出了一種基于非接觸式大空間三維激光掃描儀（3DLaserScanning，3DLS）來獲取復雜地質體表面信息并將其轉化為三維數值計算模型，使用有限差分法（finitedifferencemethod，FDM）探討了危巖體的穩定性。Meng等[8]通過3DSL方法獲取的點云信息建立了圓覺洞復雜裂隙下的高精度圍巖分析模型，研究其穩定性。然而，三維掃描獲取的點云數據，動輒百萬個甚至上億個，含有較多噪聲點，為兼顧模型的精確性和數據處理的快速性，點云數據處理至關重要。Wei等[9]提出了一種基于激光掃描獲取的點云分割方法并建立云岡石窟18窟的曲面模型。趙夫群等[10]基于層次化，依據點云曲率實現對文物點云數據的去噪濾波處理。Fleishman等[11]將雙邊濾波函數應用于3D不規則網格去噪，在平滑光順區域時較好保存了特征區。曹爽等[12]提出了利用鄰域點信息來判斷任意區域是否為特征區域，并賦予其不同雙邊濾波因子，實現基于特征選擇的雙邊濾波點云去噪。王曉輝等[13]基于法向量距離識別分類方法，采用加權局部最優投影算法和雙邊濾波分別對平滑區域和尖銳區域去噪。其中，雙邊濾波具有保持邊緣與降噪平滑的功能，可處理具有復雜特征的石質文物點云信息。但國內外石質文物數值模擬穩定性研究多聚焦于大型洞窟數值模擬結果中穩定性評價，更注重整體穩定性，而對于具有復雜特征的石質文物模型建立精細程度不夠，目前未形成一套有效的包含精細化模型建立在內的穩定性分析方法。

鑒于此，研究以大足石刻目犍聯菩薩像為例，基于雙邊濾波法，結合PCL平臺及GeomagicStudio逆向工程軟件實現石質文物的精細化構建，總結出提高石質文物逆向建模模型精度的方法，并對建立的精細化目犍聯菩薩像模型運用強度折減法定量分析石質文物穩定性，形成一套適應于具有復雜幾何特征的精美石質文物數值模擬建模及分析方法。期望對具有復雜幾何特征的石質文物穩定性研究提出科學的指導意義。

1 大足石刻目犍連菩薩像概況

大足石刻最早雕刻于唐永輝至乾封年間（公元650-668年），先后雕鑿了北山、石篆山、南山、石門山、寶頂山等眾多摩崖造像和銘文，現存造像共計5×104余尊，其中北山、寶頂山2處規模最大、保存完整、造像精美[14]。

目犍連菩薩像（下稱“菩薩像”）位于寶頂山大佛灣東側釋迦牟尼涅槃圣跡圖（即臥佛）處，菩薩像整體高度2.03m，是臥佛前眾多半身菩薩造像之一（如圖1所示）。2018年，菩薩像背部已產生大面積水平裂隙，此時頭部左側裂隙深度小，僅存在于表面；2019年，菩薩像頭部裂隙向內延伸，與背部向上延伸裂隙相貫通，穩定性受到威脅，故采取臨時支護進行保護；2021年，為增強支護文物觀賞性，更換臨時支撐，采用木支護結構，將菩薩像面部形象再次展露于游客。現階段菩薩像頭部左側及后側出現大范圍開裂，左側花冠處的裂隙L1向下延伸至左耳耳垂，與頭部后裂隙貫通，并于1.46m高度裂隙張開處發現內部有一鋼筋，外露鋼筋表面局部已發生銹蝕；后側裂隙向下張開至背部高度0.84m處，由此向2側各產生一水平裂隙（如圖2所示），菩薩像裂隙統計見表1。

2 石質文物逆向建模

基于改進后的雙邊濾波法，在PCL軟件開發平臺下，結合Geomagic逆向建模手段，對獲取的目犍連造像點云信息進行處理，在提高建模效率的同時，增加穩定性分析結果的可靠性，具體流程圖如圖3所示。

2.1 法向量估計

點云的法向量估計常采用主成分分析法[15]，為求任意點法向量，需要選取鄰域范圍的近鄰點擬合平面。基于局部點云信息上1組坐標建立點集，擬合平面方程為：ax+by+cz=0（去中心化，且a2+b2+c2=0），式中，（a，b，c）為該平面的單位法向量，為使平面擬合效果最佳，使得近鄰點更貼近平面，目標方程如下

其中：di為鄰域點到擬合平面的距離，第i個近鄰點坐標為（xi，yi，zi）。

根據條件極值存在的必要條件

將式（2）用矩陣方程AX=0表示

令式（3）為AX=0，ATA=0的最小特征值對應的特征向量即為AX=0的最小二乘解，即為該點點云的法向量。此時確定的法向量具有二義性，需要為點云集中的1點設定方向，以臨近點法向方向為參考（如圖4所示），定義為

2.2 雙邊濾波法

三維點云的雙邊濾波（bilateralfilter，BF）是二維圖像雙邊濾波的拓展[16]。采用加權平均的方式，依靠局部鄰域特征信息，聯合色度差權重因子和空間距離權重因子，在實現對噪聲光順的同時，較好保留模型的細部特征。雙邊濾波可定義為

式中：pi為點云數據P中任意一點；p′i為濾波處理后的點云；α為雙邊濾波權重因子；n為法向量。一般情況下，雙邊濾波采用高斯分布的加權平均，分別為中心點與其鄰域半徑r內的鄰域點空間域權重因子σd和頻率域權重因子σs，σd，σs分別反應鄰域范圍內的點對中心點p沿切線和法線方向的影響。其對應的高斯核函數可表達為

在雙邊濾波的運用過程中，研究人員發現依賴鄰域特征信息[17]。因此Digne[18]提出改進后的雙邊濾波因子，以便于處理大尺度噪聲。從式（9）中可以得出：當點云信息確定后，鄰域點與中心點p距離差值越大，影響越小；σd一定時，σs越大，不同頻率對應的高斯函數值越大，雙邊濾波近似退化為高斯濾波，可以看出雙邊濾波在保留特征信息時是跳躍變化的過程。

2.3 點云去噪實施

在VS2017和PCL1.11.1開發平臺下，對大足目犍連菩薩像點云數據進行去噪實驗，雙邊濾波輸入參數如表2所示。

原始數據與經過雙邊濾波處理的菩薩像如圖5所示。與原始點云數據相比，濾波后菩薩像的主要特征得到了保留，不同參數下的模型表面噪聲均有一定程度的降低，例如，服飾表面褶皺留存明顯，平滑區域內小凹陷、小突起的數量減少。參數K越大，σd越大，表面光順效果越強；當σs較大時，表面微小信息保留越多；最鄰近點個數K的變化對點云去噪效果影響明顯，臨近點個數量選取越小，特征信息保留越多。其中Model_4較其他結果的表面信息處理效果更好，構建的菩薩像模型忽略了“無意義的點”，同時較為全面地描述了特征信息。選擇Model_4進行模型建立。

2.4 點云重采樣

掃描后石質文物的點云數據密度各不相同，且數據拼接過程中部分點云信息重復疊加，為減少建模工作量，對其進行重采樣是必要的[19?20]。不同方法下的采樣效果顯示（如圖6所示）：直接進行70%曲率采樣的菩薩像表面出現多處孔洞，對其進行大量修補及修補效果的檢驗是不可行的；而采用70%隨機采樣表面噪聲過濾效果較差。采用雙邊濾波去噪后的模型在70%重采樣下，與單獨的曲率采樣/隨機采樣相比，信息保留更完整，在減少點云數量的同時保留菩薩像表面特征信息。

2.5 目犍連菩薩像建立

將上述預處理后的模型信息，在GeomagicStudio平臺下實現三維精細化曲面模型建立。曲面模型菩薩像偏差分析結果表明其均偏差為-5.0～4.0mm，標準值為1.1mm，均方根估計誤差1.1mm（如圖7所示）。在三維數據掃描期間（2018年），菩薩像背部裂隙張開度小，頭部裂隙未貫通至耳部，裂隙結構面平整性差，掃描儀難以獲取內部裂隙信息。因此，在曲面模型導入到Rhnio中完成結構面添加時，主要依據裂隙現場勘察數據（如表1所示），結合三維激光獲取的菩薩像表面裂隙分布形態，重新構建裂隙于對應位置，建立裂隙時將其進行必要的簡化處理，刪除過于尖銳的裂隙幾何特征。利用Hypermesh軟件劃分高質量2D表面網格，網格尺寸由0.15～0.005m過渡變化，在此基礎上完成三維有限元模型建立，共106464個單元，模型如圖8所示。

2.6 力學參數確定

由于文物本體取樣困難，且造像經歷多次修復，內部存在數量眾多的結構面或破碎巖體，獲取原位樣本進行試驗存在諸多困難，難以獲得可靠的力學參數。近年來，研究人員將Hoek-Brown準則引入力學參數分析過程中，以室內力學試驗為基準，通過巖體的超聲波測量結果對力學參數進行修正[21?22]。

以大佛灣巖樣的物理力學試驗結果（如表3所示）[23]為基準，通過對菩薩像原位進行表面硬度測試（Equotip3（D型））及超聲波探測儀探測（ProceqPunditLab+），儀器名稱內均勻性，并對力學試驗結果進行修正。由于菩薩像頭部表面起伏變化大，超聲波換能器無法緊貼造像表面，故僅測量其身體部分，波速測點及對應波速結果如圖9所示。

由菩薩像波速及硬度測試結果可知（如圖9～10所示），由上至下，菩薩像表層硬度、波速均呈現減小趨勢，且背部測量結果小于正面測量結果。其中，硬度云圖局部出現的幾處較大“孤島”，與波速測試結果表現出一致性。因此，參考菩薩像頭部硬度結果，對無法進行波速測量的部分以硬度變化規律作為補充。

根據室內試驗測得臥佛高度內新鮮砂巖縱波波速為3541m/s，但菩薩像下部測點波速主要集中于1600～2900m/s之間，菩薩像表面巖石完整性較差，且不同位置的力學參數差異性顯著。基于Hoek-Brown準則[24]，以波速預測巖體力學參數式（11）（12），選用波速算術平均值1953m/s作為計算依據，此時對應物理力學參數為c=25kPa，φ=26°。

局部鋼筋外露導致其工作性能劣化，材料間的黏結性能降低，本文以丁萬濤[25]對錨固支護鋼筋銹蝕的參數折減為依據，選用9%的銹蝕度，其中，E降低10%，fv或ft降低21%。

3 目犍連菩薩像穩定性分析

3.1 初始狀態穩定性分析

固定菩薩像底部邊界所有節點的X，Y，Z方向，以菩薩像頭部裂隙開裂之前的平衡狀態作為初始狀態，利用FLAC3D進行自重加載，并以最大不平衡力1×10-5N作為結束標志。為直觀觀察菩薩像內部受力情況，沿菩薩像中部做切面分析，圖11（a）（b）。自然狀態下，菩薩像呈現低眉頷首姿態，頭部重心外移，在頭部與身體連接處出現了2處應力集中區域。其中，菩薩像頸部受壓，項部受拉，由應力集中區域向內，最大主應力值以拋物線形式逐漸減小。菩薩像頭部位移層帶狀分布明顯，沿面部方向和位移增加。

3.2 現存狀態的穩定性分析

菩薩像中部切面云圖顯示頭部裂隙的產生使項部應力得到釋放（如圖12所示），菩薩像內部應力發生重分布。此時頭部應力最大值增量較小，拉應力集中區由項部向裂隙端移動，但總體仍表現為前側受壓，后側受拉的狀態。

以靠近背部一側的鋼筋上部端為起點向下選取一路徑（如圖13所示），鋼筋表面應力沿路徑方向增大，當到達頭部裂隙端時，拉應力迅速增大至最大值（95kPa）。與菩薩像巖體最大主應力相比，內部鋼筋受力明顯大于巖體所受應力，由此可見鋼筋的存在改善了菩薩像項部及裂隙端的受力，對其穩定性起到了一定的承載作用。

3.3 強度折減分析

3.3.1 塑性區分析

不同折減系數下菩薩像頭部的塑性區分布表明，圖14表示在折減初期，塑性區主要分布頸部、項部的局部區域，僅外側少量巖體達到剪切破壞，范圍較小未出現貫通情況。當折減系數增大，裂隙端開始出現塑性區，此時塑性區以頸部、裂隙端為起點，向內部逐漸延伸，且延伸速度逐漸加快；當折減系數Fs=1.20時，在頸部范圍內均出現較大范圍的剪切塑性區，此時頭部塑性區已基本貫通，可能出現剪切破壞。

3.3.2 最大剪應變分析

菩薩像的破壞過程受結構面抗剪強度的影響較大，頭部裂隙的產生促進了向頸部前側剪切帶的形成，同時造成兩端剪切應變增量集中，如圖14（b）所示。隨折減系數增大，剪切應變分布由裂隙端、頸部同時向內側擴大。當Fs=1.10時，逐漸呈現出拋物線分布狀態；Fs繼續增大，剪應變增量潛在滑動面增大，由7.50×10-5增大至6.00×10-4；當Fs=1.28時，最大剪應變增量增大至1.13×10-3，在菩薩像頭部形成剪切滑移帶。

3.3.3 位移分析

菩薩像水平位移與豎向位移隨折減系數變化基本趨于一致如圖14（c）所示：當Fs在1.0～1.2時，監測點各向位移量隨Fs增大緩慢變化，位移量恒定于較小數值，基本為μm級；當Fsgt;1.28后，菩薩像頭部監測點位移突增，各測點位移量級由μm級增加至mm級；在同一折減系數下，頭部前側監測點位移量大于同高度后側監測點的位移量，頂部位移變化大于底部監測點的位移變化。尤其測點1、測點3對折減系數變化最為敏感，說明頭部前側監測點位移變化敏感度均大于后側。最大位移預警值位于1.28～1.30之間，即菩薩像頭部水平（豎向）位移值達到-73.1μm（-8.99μm）后，頭部極有可能發生破壞。相較于塑性區貫通現象出現的Fs，監測點位移發生突變時的Fs略大，即塑性區先發生貫通，繼而監測點出現位移突變。在菩薩像逐漸發生破壞的過程中，位移變化主要發生于頭部前側，其半身部分及裂隙后頭部整體位移量較小（如圖15所示）。表現為同一高度處，頭部裂隙結構面前側巖體位移大于后側位移，臨空面位移大于內側位移，裂隙頂端位移大于內部端點位移；沿結構面向外，各點位移矢量方向由垂直裂隙面逐漸向下發展，到達菩薩像表面時各點傾向近似垂直向下，位移矢量方向變化快。由位移矢量變化可知其破壞過程表現為沿頭部裂隙下端點向外轉動傾倒的趨勢。結合上述位移、塑性區分析結果可知，菩薩像頭部前側穩定性最低。

4 菩薩像穩定性分析討論

現存石刻佛像大多呈現低眉頷首、頭部前傾的特點，且頭部尺寸較大，重心外移。頭部裂隙出現后，裂隙兩側巖體相對位移與折減系數變化呈正相關。基于強度折減法，以塑性區貫通為判斷依據時的安全系數為1.20，以監測點位移出現無限位移為判斷依據時的安全系數為1.28。不同判斷依據下的安全系數相差0.08，但菩薩像的破壞形式均表現為沿頭部支撐點發生傾倒破壞。

與菩薩像形態類似的倒懸危巖體，也表現出因局部凌空而重心外移的特點。自由端的巖體在重力作用下所產生的復合彎矩，使上部巖體隨力學參數降低產生張拉裂隙，降低危巖體穩定性[26?27]。當危巖體內部出現裂隙后，受結構面切割的危巖體在自身重力的持續作用下，使得裂隙端拉應力增大，沿內部損傷面產生向支撐角點方向的破壞（如圖16所示）。

依據現階段對菩薩像的臨時支護情況（如圖1所示），在頭部底端施加法向約束。與未支護效果相比，監測點水平位移、豎向位移均得到了大幅度減小；破壞巖體由頭部結構面前側轉為后側部分（如圖17所示）。當Fs=3.88時，頭部后側巖體產生向臨空面的破壞。因此，對于菩薩像來說，雖然臨時支護增加了頭部的穩定性，但從文物的長期穩定性來說，臨時支撐僅對頭部裂隙前側巖體進行加固，使得破壞巖體轉移至后部，并未增強頭部穩定的一致性。內部鋼筋為頭部的長期穩定性提供了積極作用，但在潮濕多雨的環境下將加速鋼筋的銹蝕，會導致連接巖體的開裂及銹脹裂隙的發展。這些對菩薩像的長期保存產生不利影響，因此有必要采取合理的加固手段以增強文物價值的延續。

5 結論

1）將基于法向量的雙邊濾波方法應用于建模前的去噪處理，經過目犍連菩薩像應用及偏差分析，結果表明良好的去噪方式可以減少建模困難度，完整保留文物特征信息。去噪后的點云模型分布均勻程度得到改善，增強平滑區光順效果，同時較好地保留有尖銳特征信息。

2）點云處理后得到高精度曲面模型，基于此建立包含裂隙結構面在內的菩薩像三維精細化實體模型，并采用基于巖體波速的Hoek-Brown準則預測得到石質文物巖體力學參數，形成了一套“模型前處理-模型建立-巖體參數獲取-穩定性模擬分析”高精度石質文物逆向建模及數值模擬流程。

3）重力荷載下，菩薩像頸部、項部兩側應力差異大。隨著裂隙產生，應力發生重分布，最大拉應力集中區向裂隙端移動。以塑性區貫通、發生無限位移為破壞依據，對菩薩像力學參數進行強度折減。結果表明：自重荷載下，Fs=1.20時，菩薩像頭部塑性區貫通；當Fs=1.28時，頭部前側位移達到臨界值，此時增加安全系數，位移敏感度劇增。即菩薩像以項部為支撐點，沿裂隙端發生傾倒破壞。

4）將數值模擬結果與含豎向裂隙的危巖體失穩破壞進行對比，危巖體在自重荷載下裂隙端拉應力集中，沿內部損傷面向支撐角點方向產生破壞，驗證了本文理論方法在菩薩像穩定性分析中應用效果的合理性。