公路隧道巖溶病害探測數值模擬研究*

鐘安然，王應偉，陳海濤

(1.昭通市交通建設工程質量安全監督局，云南 昭通 657000； 2.昭通市宜昭高速公路項目指揮部，云南 昭通 657000；3.中鐵大橋科學研究院有限公司，湖北 武漢 430034)

0 引言

巖溶災害是目前國內外隧道建設面臨的重難點問題。在巖溶區施工的隧道常發生突泥、突水、塌方等現象，造成人員傷亡、施工設備損壞，嚴重影響工期。因此，在隧道施工過程中及時進行掌子面前方及底板巖溶探測，可起到指導隧道安全施工的重要作用。目前用于隧道巖溶探測的方法主要有隧道超前地質預報法、探地雷達法、水平聲波剖面法、超前鉆孔法等[1]，不同方法具有各自適應性和缺點，探地雷達法因具有對施工影響小、預報效率高等特點，在隧道巖溶探測中廣泛應用[2]，但在探地雷達工程檢測和資料解釋中，受現場檢測條件、環境噪聲干擾及技術人員經驗制約，常導致在具體異常解釋推斷方面出現分歧，且缺乏足夠的數值模擬和實測典型數據支持[3]。因此，筆者將數值模擬與工程實踐相結合，為探地雷達技術在隧道巖溶病害探測中的應用進行研究。

1 基于GprMax軟件的正演模擬基本原理

宏觀上，所有的電磁現象均可通過Maxwell方程組進行描述，早在1966年，Yee[4]通過將Maxwell旋度方程引入空間離散方式，將其由微分方程轉化為差分方程，并模擬出理想導體在時間域的電磁響應，經完善后發展形成時域有限差分(finite difference time domain，FDTD)模擬電磁場傳播的數值計算方法。

假定在空間中一定區域內沒有電磁場場源，且區域內介質是各向同性的，那么Maxwell方程組中的2個旋度方程在該特定區域內可寫成[5-6]：

(1)

(2)

式中：E為電場強度(V/m)；H為磁場強度(A/m)；ε為介電常數(F/m)，σ為介質電導率(S/m)；μ為磁導率(H/m)；ρ為磁損耗磁阻率(Ω/m)；t為時間(s)。

在隧道巖溶探測中，模擬的對象是二維剖面，且探地雷達探測中主要利用TM(橫磁模式)電磁波[3,6]，則Maxwell方程組中的2個旋度方程在直角坐標系中的分量形式可簡化為[5]：

(3)

(4)

(5)

式中：Hx，Hy分別為磁場x，y向分量；Ez為電場z向分量。

由式(3)～式(5)可知，TM電磁波僅有Hx，Hy，Ez分量。運用Yee氏網格模型，利用中心差分代替對時間、空間坐標的微分，將連續變量離散化，即可推導出探地雷達正演模擬方程[7]。在二維情況下，Yee氏差分網格如圖1所示。

圖1 二維Yee氏差分網格示意

2 二維時域有限差分基本參數選取原則

2.1 空間步長的確定

在時域有限差分網格中，數值模擬的傳播速度將隨頻率改變，導致非物理因素引起的脈沖波形畸變、人為的各向異性及虛假的折射現象，即出現數值色散現象，造成數值不穩定[7]。

二維空間中TM電磁波數值色散方程為[8]：

(6)

式中：kx，ky分別為波矢量沿x，y向的分量；ω為角頻率；v為被模擬均勻介質中的光速；Δt為時間步長；Δx，Δy分別為x，y向空間步長。

當Δt，Δx，Δy均趨于0時，色散可減小至任意程度，但由于計算機內存空間及運算速度等因素限制，時間步長和空間步長不可能無限小，因此需選擇合適的時間步長和空間步長。通常在時域有限差分中，網格空間步長Δl最大值不超過電磁波波長λ的1/10。因數值色散引起的誤差是可接受的，通常取Δl=0.1λ，有效數字保留至小數點后三位即可，本研究案例據此取Δl=0.030m。

2.2 時間步長的確定

早在1975年，Taflove[9]探討了Yee氏差分算法的穩定性問題，并給出了時間步長Δt的限定條件。對于二維TM電磁波正演問題，其對應的穩定條件為[10]：

(7)

式中：c為真空中的光速。

考慮到探地雷達正演計算的簡便性，可令Δl= Δx=Δy，繼而有：

(8)

2.3 吸收邊界厚度的確定

時域有限差分法需在電磁場全部空間建立Yee氏網格計算空間，但由于計算機內存空間有限，不可能利用計算機直接在無限大的網格空間中計算電磁場，因此在計算過程中須在某處將網格空間截斷，使之成為有限空間[11-12]。為使電磁波不在差分網格截斷處發生明顯反射，須設置吸收邊界條件，將傳播至邊界處的電磁波吸收，進而確保計算機模擬的有限空間與自然界無窮大空間之間的差異達到最小[13]。

在PML邊界條件中，層厚通常為3～9個網格單元厚度，本文在進行探地雷達時域有限差分正演模擬時，采用PML吸收邊界條件，邊界層選取10個網格單元。

3 隧道巖溶數值模擬分析

巖溶發育種類繁多，形態各異，成因不一，影響巖溶發育的主要因素包括巖性、巖層產狀與構造、流水及地下水作用。巖溶發育過程可簡化為：溶蝕裂隙→溶孔→溶洞，溶洞內充填的物質主要為空氣、水、黏土(夾雜碎石)等。

3.1 單個溶洞模型的建立

為驗證數值模擬方法的正確性，并了解灰巖中典型巖溶病害波場特征，對隧道底板中單個圓形溶洞內全部充填空氣、全部充填土、全部充填水、上半部分充填空氣下半部分充填水、上半部分充填空氣下半部分充填土、上半部分充填水下半部分充填土的模型進行了正演模擬。各模型參數如表1所示，相對介電常數參考實際圍巖情況而定。模型1～6模擬的地質模型如圖2所示。

表1 單個溶洞模型參數

圖2 模型1～6模擬的地質模型

3.2 單個溶洞正演模擬與波場特征分析

根據實際檢測使用的參數，模擬時采用的中心頻率f=100MHz，采樣時間長度為150ns，道間距為0.3m。模型尺寸為7m×15m(深度×寬度)，模擬網格Δx=Δy=0.02m，總道數為45。PML吸收邊界采用10個網格單元厚度，激勵源選用雷克子波。模型1～6正演雷達剖面如圖3所示。

圖3 模型1～6正演雷達剖面

模型1雷達響應特征明顯，正演雷達剖面主要表現為：洞頂(剖面長度7m、雙程走時49ns處)附近呈明顯的弧形正反射同相軸；洞底(剖面長度7m，雙程走時52.5ns處)為較弱的負反射同相軸，且反射能量低于洞頂；洞頂、底界面波形未完全分離。模型2雷達響應特征明顯，正演雷達剖面主要表現為：洞頂(剖面長度7m、雙程走時49ns處)附近表現為能量較弱的弧形負反射同相軸；洞底(剖面長度7m、雙程走時70ns處)為較弱的正反射同相軸；洞頂、底界面波形完全分離。

模型3雷達響應特征明顯，正演雷達剖面主要表現為：洞頂界面(剖面長度7m、雙程走時49ns處)雷達波反射同相軸能量較空洞和土洞強，呈弧形強負反射同相軸，且洞底界面(剖面長度7m、雙程走時91ns處)反射波更清晰，并與頂界面完全分離。這是因為模型中設定的介質電導率小，雷達波在純水中的衰減較慢，而水與圍巖的介電常數差異大，洞頂界面處雷達波由圍巖進入水中反射系數為負，極性發生反轉，導致頂界面處呈強負反射。

由于雷達波在水中的波速最小，雷達波穿越溶洞的雙程走時最大，故充水溶洞頂、底界面分離程度最大，然后為充土溶洞，最后為充氣溶洞。這說明當溶洞內充填單一介質時，根據雷達剖面反射特征可識別出溶洞頂部位置，根據充填介質的不同，可較好地識別出溶洞底部位置，以便對溶洞尺寸進行分析。

由圖3可知，洞頂(剖面長度7m、雙程走時49ns處)界面反射均呈弧形同相軸形態，洞內2種充填介質的分界面(模型4剖面長度7m、雙程走時52.5ns處，模型5剖面長度7m、雙程走時52.5ns處，模型6剖面長度7m、雙程走時91ns處)較易識別，但洞底反射界面較難識別。由此可推斷，當洞內充填介質數量≥2時，可確定洞頂界面，但不易確定洞底界面，故無法準確對溶洞實際尺寸進行評估。

綜上所述，模型1～6正演雷達波形特征如表2所示。

表2 模型1～6正演雷達波形特征

3.3 巖溶裂隙模型正演模擬與波場特征分析

根據常見的巖溶裂隙發育特征，在模型0的基礎上添加1組全部充填土的巖溶裂隙模型，即模型7(見圖4)。模型內設定1條斜向巖溶裂隙，裂隙寬度為1.0～1.5m。模型尺寸為7m×15m(深度×寬度)，模擬網格Δx=Δy=0.02m。模擬時采用的中心頻率f=100MHz，采樣時間長度為150ns，道間距為0.3m，總道數為45。PML吸收邊界采用10個網格單元厚度，激勵源選用雷克子波。

圖4 模型7

模型7正演雷達剖面如圖5所示。由圖5可知，剖面內有1條明顯的傾斜同相軸，同相軸呈強反射特征，整體連續，表明模擬結果與設定模型基本吻合，但正演雷達剖面中的裂隙傾角小于模型中設定的傾角，這是由于受裂隙傾角較大的影響，雷達記錄的異常反射信息并非來自該點的正下方，從而產生偏移現象。為得到正確的解釋結果，須進行偏移歸位處理，常見的偏移處理方法包括Kirchhoff偏移、時域有限差分偏移、F-K偏移等[14]。

圖5 模型7正演雷達剖面

4 隧道巖溶探測實例與分析

某在建高速公路隧道隧址區巖溶發育強烈，施工過程中偶有底板坍塌導致施工機械被困的安全事故發生，為探明隧道底板巖溶發育情況，采用PULSE EKKO PRO型探地雷達，配置100MHz天線進行探地雷達巖溶探測，測線縱向布置于隧道底板中線，長60m，時窗設定為350ns，點測模式，點距0.3m。

該段底板巖溶探測探地雷達剖面如圖6所示。由圖6可知，剖面內有1條明顯的斜向同相軸，同相軸呈強反射特征，整體連續，同相軸下方伴有多個具有明顯強反射特征的弧形繞射，初步推斷本次探測范圍內同時發育溶洞與裂隙，且巖溶裂隙向隧道底板延伸，在雷達剖面長度15.5m處，巖溶裂隙距隧道底板最小距離約為2m。

圖6 底板巖溶探測探地雷達剖面

為驗證正演模擬在實際探測中的有效性，根據實測雷達剖面初步推斷結果，進行地質模型正演模擬，建立模型8，在其內設定1條斜向巖溶裂隙及14處埋深不一、半徑各異的溶洞，模型總長60m，深16m，選取探地雷達天線的中心頻率為100MHz，天線距為0.6m，道間距取0.3m，時窗設定為350ns，共模擬采集了200道數據，每道數據各4 947個采樣點?？紤]到實際探測過程中，天線不可能完全緊貼隧道底板，因此，模型中將天線置于遠離地面0.02m高的位置處。由于該模型較大，為盡量節省計算時間，并保證模擬效果，選取的空間步長Δx=Δy=0.03m，模型8正演雷達剖面如圖7所示。

圖7 模型8正演雷達剖面

對比實測雷達剖面與模擬雷達剖面，發現如不考慮水平干擾波，模擬雷達剖面與實測雷達剖面特征基本相似。為驗證探測結果，在實測雷達剖面長度15.5m處進行了鉆孔及井中電視驗證。井中電視結果如圖8所示，由圖8可知，從深度2.1m起，有明顯溶洞出現，說明探地雷達探測結果準確，正演解釋推斷合理。

圖8 井中電視結果

5 結語

1)通過分析GprMax正演模擬理論基礎，依據數值色散方程給出了時域有限差分法空間步長確定原則。依據解的穩定性條件，確定了基本參數時間步長。根據PML邊界條件，給出了吸收邊界層厚選取的建議值。

2)采用GprMax軟件對隧道底板下方不同類型溶洞模型進行了正演模擬。模擬結果表明，當隧道底板下方基巖中存在溶洞時，因圍巖與溶洞充填物交界面兩側存在較大的電性差異，從而易形成強烈的反射波和繞射波，這種回波在時間剖面上表現為不同極性的弧形同相軸，弧形同相軸頂部為溶洞距隧道底板的最小距離；當溶洞內部為單一介質充填時，可通過識別洞頂、底界面位置評估溶洞尺寸；當溶洞內部充填的介質數量≥2時，可確定洞頂界面，但不易確定洞底界面，無法準確對溶洞實際尺寸進行評估；當基巖中存在巖溶裂隙發育帶時，會產生相應的回波，常呈現為非水平的較連續同相軸特征；當巖溶裂隙傾角較大時，雷達剖面中顯示的異常信息與實際情況存在偏移現象，需進行偏移歸位處理。

3)探地雷達隧道巖溶探測工程實例解釋結果得到了鉆孔電視的驗證，表明探地雷達正演模擬有助于探測成果圖像中異常的識別，尤其是在灰巖地區，利用探地雷達進行隧道巖溶探測是有效的。