基于探地雷達和位移反分析的灰巖邊坡穩定性數值模擬

溫世儒，吳 霞

(1.江西理工大學 土木與測繪工程學院，江西 贛州 341000；2.江西應用技術職業學院 建筑工程學院，江西 贛州 341000)

0 引言

在我國廣西、貴州、云南等西南省份以及廣東、江西等省份的部分地區存在較為廣泛的灰巖地層。灰巖是一類可溶解的碳酸鹽類巖石，在流水的作用下易產生溶解、分崩和破碎并逐漸形成溶洞(穴)、巖溶管道和破碎帶等溶蝕現象。水流是灰巖產生溶蝕作用的關鍵因素，與季節性、周期性的地表徑流相比，豐富發達且持久的天然地下水動力場是灰巖溶蝕作用的主要貢獻者[1]。因此，灰巖地層的地質條件與風化程度在時間和空間上存在顯著的變異性，地表的地質條件往往與內部地層之間存在較大的差異，且這種差異性隨著時間的推移將變得更加復雜多變。理論研究和工程實踐表明[2]，正是由于上述顯著的時空變異性導致灰巖邊坡與花崗巖、大理巖、白云巖等難溶解的巖體邊坡相比，在施工擾動、降雨沖刷等因素的影響下更易失穩下滑，危及施工安全，阻礙施工進度。對此，在施工階段對其穩定性進行預測評估成為灰巖邊坡工程設計與施工中的重要內容。

目前，理論估算法、現場監測法、土工模型試驗以及數值模擬分析法是邊坡工程穩定性評估預測的常規方法。在工程實踐中，通常會采用多種方法進行預測以實現相互驗證及提高預測準確性的目的，如現場監測法常與數值模擬分析法相互配合應用。數值模擬分析法涵蓋多種二級方法，如：有限元法、有限差分法、離散元法、邊界元法和混合元法，其基本思路均為需將巖土體進行網格單元化處理并通過力學參數賦值的方式進行抗力設置[3]。因此，在進行數值模擬分析時，需準確掌握巖土體的物理力學參數，這是保證模擬分析結果準確性的關鍵步驟。

為了保證設計、施工的時效性和力學參數獲取的便捷性，通常會依據相關規范或者根據工程類比對巖土體的力學參數進行賦值[4]。不可否認的是，這種做法雖然簡單、快速，但有一個難以避免的缺陷，即對技術人員本身的個體主觀性和經驗性具有極大的依賴性，即便面對相同的巖土體，不同的技術人員也極有可能會選擇完全不同的力學參數，導致分析結果離散化甚至完全失效，這顯然難以乃至無法為施工過程提供有效的參考信息。此外，天然狀態下的真實巖土體，其富水性、破碎性在三維空間里都是不均勻的，地表特征與深部特征之間存在差異,且規模越大差異越明顯，這在灰巖溶蝕地區尤為典型。在這種條件下，簡單籠統地根據地表特征難以準確地對巖土體進行級別歸類和參數賦值，也很難保證計算結果的有效性及為工程實踐提供指導信息。

為了應對這一不足，一種以巖土體實際位移、應力為原始依據并通過一定的數學模型對力學參數進行反算的方法被提出，也就是反分析法。巖土體的實際位移、應力和數學模型是反分析法的兩個關鍵客觀要素，在分析過程中可較好地避免個體主觀性和經驗性的影響，近年來得到了廣泛應用，如：凌同華等[5]以里巖壟隧道YK11+150斷面為實際依托，采用監控量測獲取了圍巖的位移取值,利用一種改進的粒子群算法和BP神經網絡建立了待反演參數與實測位移值之間的非線性關系，并對圍巖的力學參數進行了反演，正交試驗的結果表明其相對誤差較小；徐中華等[6]依托上海長峰大酒店基坑工程，依據基坑的實測變形，綜合采用UCODE反分析軟件和ABAQUS有限元分析在二維平面內對基坑土層的豎向抗力單參數m值進行了反演，并以上海地區20個基坑工程為實例進行驗證，獲得了上海地區典型土層的m值及其取值范圍；Caudal等[7]以加拿大溫石棉礦東墻的尾礦邊坡為實際依托，綜合采用遙感與數值模擬技術進行了滑坡預測，并采用改進權重系數與傳遞函數的方法對集群應力反分析模型實施了優化，同時對尾礦大型巖塊的深部應力和三維位移分別實施了反求與模擬預測；張繼勛等[8]依托某水電站引水隧洞工程，基于遺傳算法基本模型編寫了改進程序,對ABAQUS商業軟件進行了二次開發并將其應用于地下洞室的圍巖力學參數反分析，發現圍巖的反演參數大于設計值，利用反演參數進行計算時位移值小于初始參數的計算值；Mahmoudi等[9]為了提高卡拉季(位于伊朗德黑蘭以西36 km處)地鐵隧道圍巖參數的估算精度及為隧道維護提供指導，利用FLAC3D軟件中的備選單變量反分析程序，基于隧道前期監測獲得的位移值對圍巖的彈性模量進行了反分析，實踐表明基于該參數可將初期噴射混凝土的厚度減少43%；為了根據有限的現場測量數據對初始應力進行反分析，Pu等[10]利用多輸出決策樹回歸器(DTR)建立了初始地應力場與其影響因子之間的關系模型，使用系數優化糾正子模型策略建立并計算獲得了400個DTR全尺寸有限元模型訓練樣本，結果表明實測值與反分析值的相關系數r達到了0.92，證明了采用DTR及神經網絡存儲全局初始地應力場的初步可行性；納啟財等[11]以姜路嶺隧道為依托，采用室內直剪壓縮試驗獲取了圍巖的天然密度和含水率，并借助數值反分析的手段獲得了圍巖的彈性模量及其增量取值。毫無疑問，上述研究對于豐富、完善反分析方法與技術具有重要的理論參考價值和實踐意義。

然而，諸如上述反分析研究，多數主要集中于算法模型，如函數改進、迭代簡化、算法編制等模型優化，忽視了巖土體本身物理性質的時空變異性，且多以單參數或兩參數反求為主，忽視了多參數之間的相互影響。這顯然不符合天然巖土體的真實條件，也就限制了研究成果的進一步推廣與應用。

為了對上述不足進行完善改進，基于理論分析，以廣西宜州-河池高速公路No.2標段某灰巖邊坡為依托，應用探地雷達實測和BMP90位移反分析模型對灰巖邊坡巖體的力學參數進行反求，并據此采用MIDAS-GTS有限元數值模擬軟件對邊坡的穩定性進行分析，以期為溶蝕區灰巖邊坡的力學參數求取與數值模擬計算提供參考。

1 工程概況與研究方案

1.1 工程概況

宜州-河池高速公路是廣西壯族自治區高速公路網絡規劃“四縱六橫三支線”的重要組成部分，是連接桂中-桂西北的重要通道。其設計時速為120 km/h，荷載等級為公路-I級，處于黔中高原向廣西盆地過渡地帶。巖溶特征顯著，以丘陵、峰叢洼地為主，植被發育。谷地狹長且邊坡較陡，最大標高和切割深度分別達到530 m和200 m。沿線基巖出露，地表多見巖溶管道和落水洞，亦可見石芽和出露鐘乳石，局部被地表殘坡積黏土覆蓋。第四系巖性為殘坡積層(Q4el+dl)粉質黏土、黏土和碎石；基巖為石炭系上統馬平組白云質灰巖、石炭系中統可溶碳酸鹽巖、石炭系下統馬平組強風化白云質灰巖和泥盆系上統白云質灰巖。

所在地四季分明，自然降水雨量充沛，雨季集中在5—8月份。根據宜州市水文站提供的水文統計資料顯示，年平均降水量為1 470 mm。天然水動力場發達，沿線河流屬珠江流域西江水系，含龍江、大環江、小環江；地下水類型主要為第四系松散層孔隙水潛水、基巖裂隙水，對混凝土無腐蝕性。

依托邊坡位于宜州市德勝鎮境內，緊鄰測垌隧道，樁號里程范圍為K49+105～136，沿路線走向的長度約31 m，屬于路塹巖質邊坡。邊坡整體稍陡峭，坡外角約118°～130°，坡頂至坡底的高程為11.2 m，坡面風化程度不均，坡頂和左側1/3坡面以強-全風化為主，表層多為全風化土塊夾雜灰巖碎石，完整性、穩定性差；坡面中部-右側整體以中風化為主，但巖體不完整，有明顯的節理、裂隙，呈塊-塊碎狀結構；可見該邊坡的整體完整性一般且穩定性較差。

2020年5月12—20日當地連續降雨，地表濕潤。5月22日上午，邊坡中部與右側部分表層坡體產生滑塌，導致應急停車道旁的邊溝被掩埋及部分隔離護欄變形損壞。為了防止落石危險，應急部門于當天下午對邊坡實施了簡單掛網處理。圖1所示為滑坡現場，可見落石已侵入應急車道，所幸未對過往車輛造成不利影響。因此，盡快對整個坡體的穩定性進行預測評估是接下來需要完成的重要工作，以便為后續處治提供相關參考。

圖1 滑坡現場

1.2 研究方案

通過坡體表層地質素描和現場踏勘可知，該邊坡的淺層完整性和風化程度顯著不均。對于整個邊坡淺層坡體而言，全風化、強風化和中風化巖體的比例大致相當，無論將其視為何種級別巖體均不合理。此外，該邊坡屬于溶蝕區灰巖邊坡，當地巖溶作用特征顯著，天然的地下水動力場彼此輻射相連，坡體內部易發育溶洞(穴)、破碎帶等不良地質體，僅依靠表層地質調查難以發現此類不良地質體,從而影響邊坡的整體性和穩定性預測結果。

為此，摒棄以往對整體坡體進行一次性定級歸類的做法，首先進行坡面地質素描并采用探地雷達對邊坡內部的富水性、破碎性進行探測，然后再依據《工程巖體分級標準》(GB/T 50218—2014)的建議對邊坡進行分塊與單元級別劃分。

一般而言，當巖體的級別得到確定時，可以便捷地根據規范進行力學參數的賦值。但這忽視了巖體物理特征的時間-空間變異性，沒有考慮力學參數隨時間-空間的動態變化性，更多的在于依據工程類比和地區經驗實施靜態賦值，這顯然不符合工程實際。

根據巖體本構關系可知，塊體的形變與位移是彈性抗力與效應二者之間的平衡結果，而巖體的抗力又與其自身的抗剪強度指標、泊松比等參數具有直接關系[12]。此外，由于接觸與擠壓，塊體之間存在位移和應力的傳遞，從而產生位移和應力的疊加效應。工程實踐表明，這種位移和應力是可以通過在巖塊上安設具有一定錨固長度的可產生協調變形的柔性構件進行測量估算的，一旦獲得了各個級別巖塊的位移或者應力取值，就可以通過巖體本構關系對巖塊的力學參數進行反求，最終得到不同巖塊的力學參數。

為此，基于前述探地雷達探測及巖塊單元劃分，引入位移反分析技術對完整邊坡不同單元巖塊的力學參數進行反演。將得到的力學參數應用于有限元數值模擬計算，從而對邊坡的穩定性進行預測。最后，采用現場監控量測對模擬分析結果進行對比驗證。技術路線如圖2所示。

圖2 技術路線

2 探地雷達探測

2.1 探測原理

探地雷達是近年來在土建領域被廣泛采用的一種時效性強、設備簡易、操作便捷的電磁勘探技術，其有效探測深度可達20～35 m。探地雷達系統由3部分組成，即：控制電腦、連接電纜和發射天線。探測時，首先在控制主機中設定好發射參數，然后發射天線將根據該參數向目標地層發射具有初始頻率、相位、能量和振幅的入射電磁波。

與地震波等機械波不同，電磁波本質上屬于交替變換的電場-磁場所形成的統一場在空間中的傳播[13]。天然狀態下的真實巖土體是電介質(既有導電性又有絕緣性)，土體內部不同組成物之間的接觸面屬于電性質突變帶。當入射電磁波到達該突變帶時將產生反射和折射，且這種反射和折射在土層內部將不斷循環直至能量衰減至無法繼續傳播。理論與試驗研究表明[14]，土層的物理特征會影響其內部組成物的電性質，據此可通過分析反射波的波形特征對內部土層的地質條件實施預測。

在有效探測深度范圍內，探地雷達不僅對地層中的大斷層、大空洞等大目標體具有較高的分辨率，還對節理、裂隙、小型溶洞(穴)、軟弱夾層等具有良好的識別度，對自由水也具有靈敏度高、可判性好的優勢，是短距離無損探測的代表[15]。收-發分離式探地雷達的探測原理如圖3所示。

圖3 探地雷達探測原理示意圖

2.2 現場探測與單元分塊

在前期地質勘察資料的基礎上，采用青島電波所生產的LTD-2500型探地雷達實施現場探測，相關探測參數見表1。需要說明的是，探地雷達的探測深度隨天線中心頻率(75 MHz～2.5 GHz)的提高而降低，常規地質預報選用的中心頻率通常為100 MHz，此時其探測深度可達20～35 m且同時能滿足分辨率的要求，并能完整地對邊坡內部的地質情況進行探測，故天線頻率選用100 MHz。

表1 探測參數

由于邊坡陡峭且高度較高，實測時無法用人力完成探測，故需要事先在鏟車的鏟斗上焊接好鋼架操作平臺，然后選配2名探測人員進入平臺并手持雷達天線且要求緊貼坡面。鏟車司機需嚴格控制沿道路縱向的移動速度從而使天線在坡面上緩慢移動以形成水平向測線、測點以保證探測質量。由于不易在豎向控制鏟斗的下放速度且容易導致天線脫離坡面，因而僅沿坡面水平向設置測線和測點。

測線的起點、終點和測點個數視現場操作難易而定(如坡面是否凹凸不平)，間距為1～2 m。為了完整地對邊坡進行探測，除坡面外還需在邊坡周圍3～5 m范圍內布置測點以擴大探測范圍。

探測完畢后，采用自帶的IDSP6.0探地雷達分析軟件對探測文件進行背景去除、反褶積、信號增益、0點校正等處理。剔除無效文件后，獲得線測圖像12副(如圖4(a)所示)和點測圖像37副(如圖4(b)所示)。綜合分析前期地質勘察資料和探測圖像的波形特征，對邊坡內部的地質條件進行預測，再結合淺層坡體的地質條件綜合對邊坡進行單元劃分。

圖4 探測圖像

統計表明，該邊坡共劃分為15個單元塊，其中包括6個V級單元巖塊，位于坡頂和坡體左側；7個IV級單元巖塊，位于邊坡中部和右側坡體上部；2個III級單元巖塊，位于右側坡體的中部和下部。至此，邊坡單元分塊完畢，后續還需在不同的單元巖塊上安設位移計以繼續完成單元巖塊的力學參數反分析。此外，通過圖4可知，回波反射強烈，線測圖像波形雜亂，同相軸斷裂不連續，且從圖像中段開始存在“雪花狀”特征，表明內部巖體多數較破碎，完整性差且潮濕富水[16]。

3 位移反分析

天然真實巖體由微小單元巖塊和結構面共同構成，單個巖塊自身與結構面共同形成相互約束與協調變形，最終在宏觀上表現出位移與應力特征，且可以通過安設一定的測試元件進行測量。表面應變計、表面壓力盒等面接觸式元件雖然也可測量位移和應力，但通常更適合測量接觸面之間的接觸應力與相互位移，不適用于對巖體內部全尺度范圍內的位移進行測量。因此，為了獲取不同單元塊的全尺度位移值和即時力學參數，首先需通過鉆孔在不同單元塊上安設可與巖塊產生協調變形的位移計(型號：WYJ-100，量測范圍0～100 mm)，然后基于位移本構關系和位移反分析模型反求不同單元巖塊的力學參數。

事實上，采用人工神經網絡開展巖土體位移反分析的研究與實踐相繼出現且具有一定的可行性，但對有效數據樣本的容量、樣本訓練時間具有較高的要求，時效性與便捷性較差，更適用于非緊急情況[17]。對此，采用一種對樣本數據容量依賴性更低、效率更高的“試湊式”位移反分析模型——BMP90位移反分析模型對不同巖塊的力學參數進行反求。

3.1 反求模型與原理

BMP90位移反分析模型是李世輝等人根據邊界元法基本理論研制開發的一款適用于各類型工程巖土體的位移反分析算法，早期多應用于隧道工程巖體，后來逐漸在各類型巖土體中得到了應用和完善。其基本原理為：通過不斷地調整參數取值使實測值與擬定值之間的差值滿足某容許范圍，此時擬定值視為可用，如式(1)所示[18]：

(1)

式中，R為實測值與擬定值之間的差值的容許范圍；K1，K2分別為擬定值和實測值。

因此，基于上述不同單元巖塊的位移實測數據，通過該反分析程序來獲取單元巖塊的彈性模量E、泊松比μ等力學參數，并將其應用于后續的有限元數值模擬計算中。

3.2 反求結果

根據前述方案，完成了位移計的安設并開始實施位移測量，由于當邊坡內部1 m 范圍內發育有溶洞(穴)時無法安設位移計，因此避開內部溶洞(穴)位置后共安設13根位移計。位移計的安設示意如圖5所示，鉆孔直徑為76 mm；反分析結果見表2。

圖5 位移計安設示意圖

表2 位移反分析結果

由此可見，不同級別單元巖塊的力學參數變化規律基本上與其理論本構關系一致，且即便是相同級別的巖塊，其力學參數亦不相同甚至差別較大，這較為真實地反映了不同巖塊本身所具有的即時力學特征，與根據規范對整個邊坡進行統一性賦值相比具有顯著的優勢。

4 有限元數值模擬

為了對邊坡的下滑穩定性進行預測，采用MIDAS-GTS有限元數值模擬軟件建立邊坡的網格模型并對位移進行計算，同時采用現場位移監測對計算結果進行驗證對比。

4.1 網格模型

以邊坡的起止里程樁號為參照建立網格模型，橫向寬度為31 m，豎向高度為11.2 m，坡頂往外側的范圍取5 m。根據邊界約束條件設置準則，在模型左右兩側設置X方向的約束，在模型底部設置Y方向的約束，默認邊坡表面為自由面且考慮基巖裂隙水滲流。與三邊形網格相比，四邊形網格單元具有更好的力傳遞效率，因而采用四邊形單元生成網格[19]。巖體內部力和位移的傳遞本質上依賴于塊體接觸與相互運動，因而考慮到邊坡在三維空間里的地質條件差異性，模型中的網格單元不宜采用等距離設置，而應按照1∶1.25的比例對V，IV級單元巖塊的網格進行加密處理，從而提高此類巖體內部節點對位移和力的傳遞效率以及提高運算精度。

為了節省運算存儲空間與加快運算效率，每一次迭代都設置應力釋放系數[20]。假定邊坡巖體為彈塑性體，遵守M-C應力準則。為了充分實現巖體內部應力的重新分布與平衡，視應力分配達到平衡為基本原則由軟件自行運算，因而不設置固定的迭代次數[21]。網格模型見圖6，圖6(a)為添加基巖裂隙水滲流壓力后顯示的網格模型，為了直觀顯示已隱藏了網格線,圖6(b)為未隱藏網格線的模型圖。

圖6 網格模型

4.2 結果分析

計算時，各單元巖塊分別采用兩種力學參數，一種是由前述反分析所得參數，一種是依據《工程巖體分級標準》(GB/T 50218—2014)的建議進行賦值。重度γ取22 kN/m3。具體而言，首先利用BLOCK-SET命令，將每一個單元塊設置成隔離體，由此保證獨立的隔離體之間具有互不相同的力學參數；然后再對每一個單元塊的力學參數進行依次逐個輸入；當全部彼此相互獨立的單元塊的參數都輸入完畢時，整體邊坡不同級別單元便獲得了各自獨立的力學參數，最終避免規范建議賦值所存在的盲目性和籠統性。

網格模型的位移云圖顯示的是整個邊坡的位移狀態，無法顯示不同巖塊的位移變化趨勢，對于分析單元巖塊的位移特征而言沒有意義。對此，采用GTS內置的節點提取技術對不同單元巖塊的位移值進行提取。為了加快提取速度，每個單元巖塊提取一個節點，取不同力學參數時的提取結果見表3。

表3 單元巖塊位移

由此可見，基于位移反分析所得力學參數的計算值與實測值更為接近，這說明與規范建議賦值相比，通過雷達實測、單元分塊及反分析獲取的力學參數更加可靠。此外，表中多數數值模擬計算的位移值小于實測值，這是由于網格邊界條件限定了巖體的變形范圍和方向，而真實巖體是三向變形的。

5 結論

本研究在理論分析的基礎上，以廣西壯族自治區典型溶蝕區灰巖邊坡為實際依托，采用探地雷達實測和位移反分析獲取邊坡不同單元巖塊的力學參數，并利用有限元數值模擬計算和現場量測對邊坡的位移進行了預測對比，結論如下：

(1)灰巖邊坡的溶蝕作用強烈，求取力學參數時需考慮其地質條件的典型時空變異性，宜先對邊坡進行分塊處理再實施單元塊定級歸類。

(2)基于探地雷達實測的單元劃分可較好地避免定級時的主觀性和經驗性影響，在此基礎上采用BMP90位移反分析模型可反演得到更加真實的力學參數。

(3)獲取有效的力學參數是保證數值模擬計算結果可靠性的重要前提，今后可嘗試綜合利用探地雷達和其他人工智能反分析模型對巖體的力學參數進行求取，這是一種新思路、新方法，應繼續加以重視。