公路隧道巖溶病害探測(cè)數(shù)值模擬研究*

鐘安然，王應(yīng)偉，陳海濤

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

0 引言

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

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

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

假定在空間中一定區(qū)域內(nèi)沒有電磁場(chǎng)場(chǎng)源，且區(qū)域內(nèi)介質(zhì)是各向同性的，那么Maxwell方程組中的2個(gè)旋度方程在該特定區(qū)域內(nèi)可寫成[5-6]：

(1)

(2)

式中：E為電場(chǎng)強(qiáng)度(V/m)；H為磁場(chǎng)強(qiáng)度(A/m)；ε為介電常數(shù)(F/m)，σ為介質(zhì)電導(dǎo)率(S/m)；μ為磁導(dǎo)率(H/m)；ρ為磁損耗磁阻率(Ω/m)；t為時(shí)間(s)。

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

(3)

(4)

(5)

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

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

圖1 二維Yee氏差分網(wǎng)格示意

2 二維時(shí)域有限差分基本參數(shù)選取原則

2.1 空間步長的確定

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

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

(6)

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

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

2.2 時(shí)間步長的確定

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

(7)

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

考慮到探地雷達(dá)正演計(jì)算的簡(jiǎn)便性，可令Δl= Δx=Δy，繼而有：

(8)

2.3 吸收邊界厚度的確定

時(shí)域有限差分法需在電磁場(chǎng)全部空間建立Yee氏網(wǎng)格計(jì)算空間，但由于計(jì)算機(jī)內(nèi)存空間有限，不可能利用計(jì)算機(jī)直接在無限大的網(wǎng)格空間中計(jì)算電磁場(chǎng)，因此在計(jì)算過程中須在某處將網(wǎng)格空間截?cái)?，使之成為有限空間[11-12]。為使電磁波不在差分網(wǎng)格截?cái)嗵幇l(fā)生明顯反射，須設(shè)置吸收邊界條件，將傳播至邊界處的電磁波吸收，進(jìn)而確保計(jì)算機(jī)模擬的有限空間與自然界無窮大空間之間的差異達(dá)到最小[13]。

在PML邊界條件中，層厚通常為3～9個(gè)網(wǎng)格單元厚度，本文在進(jìn)行探地雷達(dá)時(shí)域有限差分正演模擬時(shí)，采用PML吸收邊界條件，邊界層選取10個(gè)網(wǎng)格單元。

3 隧道巖溶數(shù)值模擬分析

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

3.1 單個(gè)溶洞模型的建立

為驗(yàn)證數(shù)值模擬方法的正確性，并了解灰?guī)r中典型巖溶病害波場(chǎng)特征，對(duì)隧道底板中單個(gè)圓形溶洞內(nèi)全部充填空氣、全部充填土、全部充填水、上半部分充填空氣下半部分充填水、上半部分充填空氣下半部分充填土、上半部分充填水下半部分充填土的模型進(jìn)行了正演模擬。各模型參數(shù)如表1所示，相對(duì)介電常數(shù)參考實(shí)際圍巖情況而定。模型1～6模擬的地質(zhì)模型如圖2所示。

表1 單個(gè)溶洞模型參數(shù)

圖2 模型1～6模擬的地質(zhì)模型

3.2 單個(gè)溶洞正演模擬與波場(chǎng)特征分析

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

圖3 模型1～6正演雷達(dá)剖面

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

模型3雷達(dá)響應(yīng)特征明顯，正演雷達(dá)剖面主要表現(xiàn)為：洞頂界面(剖面長度7m、雙程走時(shí)49ns處)雷達(dá)波反射同相軸能量較空洞和土洞強(qiáng)，呈弧形強(qiáng)負(fù)反射同相軸，且洞底界面(剖面長度7m、雙程走時(shí)91ns處)反射波更清晰，并與頂界面完全分離。這是因?yàn)槟Ｐ椭性O(shè)定的介質(zhì)電導(dǎo)率小，雷達(dá)波在純水中的衰減較慢，而水與圍巖的介電常數(shù)差異大，洞頂界面處雷達(dá)波由圍巖進(jìn)入水中反射系數(shù)為負(fù)，極性發(fā)生反轉(zhuǎn)，導(dǎo)致頂界面處呈強(qiáng)負(fù)反射。

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

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

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

表2 模型1～6正演雷達(dá)波形特征

3.3 巖溶裂隙模型正演模擬與波場(chǎng)特征分析

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

圖4 模型7

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

圖5 模型7正演雷達(dá)剖面

4 隧道巖溶探測(cè)實(shí)例與分析

某在建高速公路隧道隧址區(qū)巖溶發(fā)育強(qiáng)烈，施工過程中偶有底板坍塌導(dǎo)致施工機(jī)械被困的安全事故發(fā)生，為探明隧道底板巖溶發(fā)育情況，采用PULSE EKKO PRO型探地雷達(dá)，配置100MHz天線進(jìn)行探地雷達(dá)巖溶探測(cè)，測(cè)線縱向布置于隧道底板中線，長60m，時(shí)窗設(shè)定為350ns，點(diǎn)測(cè)模式，點(diǎn)距0.3m。

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

圖6 底板巖溶探測(cè)探地雷達(dá)剖面

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

圖7 模型8正演雷達(dá)剖面

對(duì)比實(shí)測(cè)雷達(dá)剖面與模擬雷達(dá)剖面，發(fā)現(xiàn)如不考慮水平干擾波，模擬雷達(dá)剖面與實(shí)測(cè)雷達(dá)剖面特征基本相似。為驗(yàn)證探測(cè)結(jié)果，在實(shí)測(cè)雷達(dá)剖面長度15.5m處進(jìn)行了鉆孔及井中電視驗(yàn)證。井中電視結(jié)果如圖8所示，由圖8可知，從深度2.1m起，有明顯溶洞出現(xiàn)，說明探地雷達(dá)探測(cè)結(jié)果準(zhǔn)確，正演解釋推斷合理。

圖8 井中電視結(jié)果

5 結(jié)語

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

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

3)探地雷達(dá)隧道巖溶探測(cè)工程實(shí)例解釋結(jié)果得到了鉆孔電視的驗(yàn)證，表明探地雷達(dá)正演模擬有助于探測(cè)成果圖像中異常的識(shí)別，尤其是在灰?guī)r地區(qū)，利用探地雷達(dá)進(jìn)行隧道巖溶探測(cè)是有效的。