基于隧道超前預(yù)報的GPR正演模擬及應(yīng)用分析

陳 松，周黎明，羅士新

(1.武漢地質(zhì)調(diào)查中心勘查技術(shù)室，武漢 430205;2.長江科學(xué)院 水利部巖土力學(xué)與工程重點實驗室，武漢 430010)

1 研究背景

我國地質(zhì)預(yù)報工作開展相對較晚，在20世紀(jì)70年代才開始進行隧道地質(zhì)預(yù)報的應(yīng)用和研究。20世紀(jì)80年代，在南嶺隧道、大瑤山隧道的掘進過程中開展了地質(zhì)預(yù)報的試驗研究，且大瑤山隧道較成功地引用了地震反射法等［1］。C.W.Zang(2004)等［2］把地質(zhì)雷達法、TSP法、超前鉆探法集合起來，對圍巖的完成性、強度以及水系發(fā)育的分級等進行了研究。齊傳生(2004)等［3］在圓梁山隧道中采用了地質(zhì)雷達法、TSP202法、HSP法、紅外探測法、超前鉆等方法綜合分析，總結(jié)出了各種地質(zhì)預(yù)報的特點、工序。曲海峰(2006)等［4］以TSP法、地質(zhì)雷達法、地質(zhì)調(diào)查法為手段，研究了綜合預(yù)報的方法。由上可以看出地質(zhì)雷達是隧道預(yù)報工作者最常選用的一種方法，因此，如何能利用好這一物探方法是預(yù)報的前提和基礎(chǔ)。

近年來也有不少學(xué)者做過探地雷達正演模擬相關(guān)的研究工作，如鄧世坤等(1993)［5］給出了用射線追蹤法合成的二維地電模型的探地雷達圖像，用衍射波公式合成了二維管狀模型的雷達圖像，并與物理模擬作了比較。吳寶杰(2007)［6］、趙宇(2009)［7］詳細介紹了地質(zhì)雷達的原理及二維正演模擬。趙峰(2012)等［8］分析了探地雷達在探測隧道襯砌空洞的有效性以及襯砌空洞的雷達波反射特征。胡俊(2012)等［9］對比分析了探地雷達對多年凍土區(qū)厚層地下冰的分布、埋藏深度、賦存情況以及對凍土類型識別的有效性。孫忠輝(2013)等［10］利用軟件GprMax模擬地質(zhì)雷達在初襯中的檢測是準(zhǔn)確可行的，能夠清晰地反映隧道襯砌中鋼拱架和鋼筋的成像特征。

在隧道環(huán)境下，多數(shù)正演模擬的對象主要集中在隧道的襯砌上，使用的天線頻率高，模擬的厚度小，消耗的時間短，最后的效果比較理想，而將隧道掌子面前方的圍巖作為正演模型研究對象的相對較少。因此，本文主要從大尺度的模型出發(fā)，通過設(shè)計預(yù)報的地質(zhì)模型，使用低頻天線來探討地質(zhì)雷達在不同模型下的波場響應(yīng)特征及其對應(yīng)的解釋方法，并通過實測數(shù)據(jù)的應(yīng)用效果對比分析，證明地質(zhì)雷達在隧道環(huán)境下進行超前預(yù)報的可靠性。

2 地質(zhì)雷達正演基本原理

地質(zhì)雷達GPR(Ground Penetrating Radar)是通過發(fā)射天線將高頻電磁波以脈沖形式從隧道的掌子面?zhèn)鬟f到隧道前方的圍巖中，經(jīng)過圍巖界面的反射返回隧道掌子面，并被雷達的接收天線接收。通過對回波信號進行分析解釋，就可以獲得隧道掌子面前方的圍巖信息。本文中正演的基本原理是基于1966年Kane S.Yee提出的 FDTD(Finite Difference Time Domain)差分算法和 PML(perfectly matched layer)吸收邊界條件，其中差分算法是把Maxwell方程組用中心差分近似替代，該方法運算速度快，可以解決各種復(fù)雜點的散射［11－12］。

假設(shè)介質(zhì)是各向同性的，在二維空間下，若Z方向上只有電場分量Ez，磁場分量Hz為0，那么麥克斯韋的旋度方程可表示為

式中:H為磁場強度(A/m);E為電場強度(V/m);μ為磁導(dǎo)率(H/m);σ為電導(dǎo)率(s/m);ε為介電常數(shù)(F/m);σm為導(dǎo)磁率(Ω/m)。

K.S.Yee根據(jù)上面的3個偏微分方程引入FDTD算法，通過在空間建立矩形的差分網(wǎng)格，一個網(wǎng)格節(jié)點對應(yīng)的時變參量為f(x，y，t)，則可以得到FDTD表達式為

式中:i，j為x，y方向的網(wǎng)格號;Δx，Δy為矩形網(wǎng)格分別沿著x，y方向的空間步長;Δt為時間步長。麥克斯韋方程要保持穩(wěn)定性和收斂性，則需要滿足公式

考慮到離散化后的頻散現(xiàn)象，如果通過減小網(wǎng)格而抑制頻散的話，就會增大相應(yīng)的計算量，所以，在考慮壓制頻散與控制計算量的情況下，可以使網(wǎng)格盡量縮?。?3－15］。

網(wǎng)格步長要小于波傳播時最小波長 λmin的1/10，則有

3 隧道模型下的正演模擬

3.1 吸收邊界

受到計算量的影響，用FDTD計算求解時，一般只求取空間有限區(qū)域的迭代解，在實際的計算中，通常用有限計算空間來代替無限的空間。實際應(yīng)用表明，PML在大的入射角上吸收的效果很好，在一般的計算中該算法會顯得很復(fù)雜，占用很大的內(nèi)存［16－19］。

3.2 模型參數(shù)設(shè)置

隧道模型參照隧道開挖面，設(shè)計了2個側(cè)壁和1個掌子面。下面將給定模擬的相關(guān)參數(shù)。我們選取的天線頻率f0為50 MHz。根據(jù)電磁波在介質(zhì)中傳播的最小波速Vmin=c/=0.6 × 108，則有模擬的最高頻率 fm=3f0=150 MHz，λmin=0.4 m。根據(jù)式(3)和式(4)，有穩(wěn)定性條件 Δx≤0.04m，Δy≤0.04m，Δt≤30 ps，給定的時間窗口為400 ns。根據(jù)模擬的效果，本節(jié)選取的網(wǎng)格步長為Δx=Δy=0.01 m。天線間距為0.22 m，移動步長為0.20 m。天線探測方向為從掌子面的左端測向右端。模型中的幾種介質(zhì)電性參數(shù)見表1。

表1 試算模型電性參數(shù)Table 1 Dielectric parameters for model calculation

本文采用的是GprMax2D軟件進行正演分析，給出了9種地質(zhì)預(yù)報的模型，具體參數(shù)見表2。

表2 試算模型屬性Table 2 Properties of model for the calculation

通過試算并根據(jù)隧道超前預(yù)報的場地條件，本文設(shè)計的隧道正演模型總長為30 m，總寬度為20 m，左右側(cè)壁寬度均為5 m，掌子面處于已開挖的5 m處，開挖空間長5 m，寬10 m。模型Ⅰ至模型Ⅵ為設(shè)計的簡單模型，分別為距掌子面前方10 m處發(fā)育有水平狀的破碎帶、距掌子面5 m處垂直掌子面的破碎帶、左傾30°的破碎帶、右傾30°的破碎帶、中心距掌子面10 m處發(fā)育有半徑為1 m的空溶洞及充水溶洞;模型Ⅶ至模型Ⅸ為花崗巖地層、頁巖地層、灰?guī)r地層與水系(水系帶中部為梯形石塊)、溶洞(分別充滿空氣、水、破碎塊)、透鏡體的組合模型。以上設(shè)計的模型在隧道開挖中會經(jīng)常遇見，具有代表性。隧道圍巖介質(zhì)中花崗巖的相對介電常數(shù)為5，電導(dǎo)率為0.001(s/m);頁巖的相對介電常數(shù)為7，電導(dǎo)率為0.01(s/m);灰?guī)r的相對介電常數(shù)為8，電導(dǎo)率為1×10－9(s/m);破碎帶的相對介電常數(shù)為25，電導(dǎo)率為0.003(s/m);透鏡體參照破碎帶的電性參數(shù)。隧道模型見圖1。

3.3 正演模擬計算及分析

把模型參數(shù)文件導(dǎo)入到GprMax2D軟件進行計算，可以得到正演模擬數(shù)據(jù);再把得到的數(shù)據(jù)文件導(dǎo)入Fortran程序，經(jīng)過增益恢復(fù)、能量均衡等初步處理生成*.grd格式數(shù)據(jù)，利用Surfer軟件成圖即完成整個正演的工作，成果圖如圖2所示。

從圖2正演成果分析可知:地質(zhì)雷達識別異常體的頂界面位置非常準(zhǔn)確，如模型Ⅰ頂界面的埋深10 m、模型Ⅱ的頂界面埋深5 m、模型Ⅴ及模型Ⅵ溶洞的頂界面埋深9 m、模型Ⅷ組合溶洞埋深7 m等位置均與正演結(jié)果對應(yīng)，但是由于電磁波受到破碎帶、水等介質(zhì)的干擾及吸收影響，地質(zhì)雷達對異常體下界面位置判斷相對有出入;同時地質(zhì)雷達對斷層破碎帶的傾角判斷也比較敏感，如模型Ⅰ、模型Ⅲ、模型Ⅳ正演結(jié)果分別對應(yīng)了水平層狀傾向、左傾30°、右傾30°的角度變化，模型 IX中傾斜斷層面(圖1中標(biāo)注為9)與正演的同相軸(圖2中標(biāo)注為9)位置也比較吻合，由于受到上層透鏡體的反射吸收影響，斷層面反射波同相軸的能量較弱;地質(zhì)雷達對圓形溶洞的形狀預(yù)測比較準(zhǔn)確，如模型Ⅴ、模型Ⅵ、模型Ⅷ正演剖面中反射波的弧形頂部對應(yīng)了溶洞的頂部形狀，但是對于模型Ⅷ中同時出現(xiàn)3個相鄰的溶洞識別效果一般;地質(zhì)雷達對于巖層的界面位置識別也較為準(zhǔn)確，如模型Ⅶ的頁巖與灰?guī)r界面3、模型Ⅷ的花崗巖與頁巖界面5及頁巖與灰?guī)r界面6對應(yīng)了正演剖面中的位置(圖2中分別標(biāo)注為3，5，6);由于模型吸收邊界沒有起到對隧道側(cè)壁邊界反射波完全吸收的作用，邊界拐點的繞射和散射干擾等對最后的正演結(jié)果有一定的影響。因此，地質(zhì)雷達在隧道環(huán)境下的正演結(jié)果與理論模型對應(yīng)得比較一致，通過以上分析電磁波的波場特征可以為解釋實測的雷達數(shù)據(jù)作理論依據(jù)，對提高解釋的精度有很大幫助。

圖1 隧道模型示意圖Fig.1 Schematic diagram of tunnel models

圖2 正演成果Fig.2 Results of forward simulation

圖3為從模型Ⅰ至模型Ⅸ正演結(jié)果中提取異常體中心位置的單道數(shù)據(jù)，通過對比分析總結(jié)有:電磁波遇到斷層破碎帶及溶洞等界面的反射，電磁波的相位會發(fā)生變化，如圖3中紅色標(biāo)注的為地質(zhì)體引起的相位變化，模型Ⅰ—Ⅳ、模型Ⅵ、模型Ⅶ電磁波由介電常數(shù)小的介質(zhì)進入介電常數(shù)大的介質(zhì)，波相位由正變負，而模型Ⅴ電磁波從圍巖進入溶洞的空氣介質(zhì)，介電常數(shù)由大變小，對應(yīng)的波相位由負變正，因此，根據(jù)波相位的正負相對變化也可以預(yù)測掌子面前方圍巖的變化情況。同時透過破碎帶、透鏡體區(qū)域后的電磁波衰減很快，波的振幅變小，隨著深度的增加能量已經(jīng)很弱。在后期資料解釋中，我們要結(jié)合單道數(shù)據(jù)綜合判斷異常體引起的波形特征，提高解釋的精度。

對比圖2和圖3可以發(fā)現(xiàn)，地質(zhì)雷達在理論的隧道模型下預(yù)測掌子面前方的地質(zhì)異常體比較準(zhǔn)確，波形圖上能夠定性地反映地質(zhì)體的形態(tài)，并能定量地給出埋深及發(fā)育的規(guī)模、傾向傾角等重要的信息，同時也可以根據(jù)單道數(shù)據(jù)的波相位變化分析反射層兩端巖體的介電常數(shù)變化，進而推斷掌子面前方巖性的變化。因此，通過正演模擬，從理論上驗證了地質(zhì)雷達在隧道環(huán)境下做超前預(yù)報的可靠性，同時也為實際應(yīng)用提供了分析波形的依據(jù)。

通過大尺度的模型計算，可以很明顯地看到計算量和計算時間在大大增加，一個剖面的計算時間大約是23 h，這也是不同于對襯砌模擬(計算1個剖面大約數(shù)分鐘的時間)的地方。同時模型尺度變大，電磁波會隨著深度的增加迅速衰減，50 MHz的雷達天線探測的有效深度在30 m以內(nèi)。

圖3 單道數(shù)據(jù)振幅圖Fig.3 Amplitude of single-channel data

4 GPR實測資料特征分析

隧道Ⅰ:此隧道區(qū)域地質(zhì)露頭多為牛屋組(Ptn)中厚層狀風(fēng)化巖層的巖體，其節(jié)理裂隙比較發(fā)育。布置的雷達測線在掌子面上的長度為13 m，雷達的天線中心頻率是50 MHz，天線間距1.5 m，天線移動步長0.25 m，時間窗口長度1 200 ns。

隧道Ⅱ:此隧道區(qū)地貌單元屬構(gòu)造剝蝕低山—丘陵地貌，隧道走向呈現(xiàn)近北東—南西向，與山脊線走向近于平行。該區(qū)大地構(gòu)造屬于昆侖—秦嶺緯向構(gòu)造體系東秦嶺南亞帶，該構(gòu)造體系是區(qū)內(nèi)成生發(fā)展最早的構(gòu)造體系，為區(qū)內(nèi)主要構(gòu)造格架，其他后期構(gòu)造在此構(gòu)造之上疊加、穿切、橫跨，構(gòu)造形態(tài)復(fù)雜，褶皺發(fā)育，斷裂交織，具有多期活動特征，是個相對復(fù)雜的構(gòu)造地區(qū)。布置的雷達測線長度為5.5 m，天線的中心頻率為100 MHz，天線間距1 m，天線移動步長0.25 m，時間窗口長360 ns。

從圖4(a)雷達剖面分析，可以發(fā)現(xiàn)在深度位置17～20 m處出現(xiàn)較為明顯的連續(xù)的強能量同相軸，且同相軸的展布范圍較大，與本文正演的模型I的波形特征比較接近，推斷此處圍巖的節(jié)理裂隙較發(fā)育，有破碎帶的存在，且破碎帶呈現(xiàn)平行掌子面的發(fā)育方向。后期隧道開挖過程中，掌子面前方17 m后巖層比較破碎，掉塊嚴(yán)重，說明巖體節(jié)理裂隙發(fā)育，跟探測的結(jié)果比較一致，達到了超前預(yù)報的效果。同時也可以看到，雷達剖面深度在30 m以后，同相軸的能量已經(jīng)很弱，無法提供有效的信息。

根據(jù)圖4(b)地質(zhì)雷達剖面圖，可以比較清楚的看到在掌子面右上方向前50～80 ns的區(qū)段內(nèi)(深度2.8～4.0 m，參考波速0.1 m/ns)，同相軸呈彎曲拋物線狀態(tài)，與本文模型Ⅵ正演的波形特征非常接近，而且信號振幅較強，與周圍波形存在明顯差異，該異常區(qū)由巖溶洞穴帶形成，結(jié)合隧道地質(zhì)調(diào)查，初步判斷此區(qū)段為巖溶洞穴。經(jīng)過后期隧道開挖，驗證了此位置預(yù)測的洞穴的存在，達到了預(yù)報的效果。

圖4 隧道雷達剖面Fig.4 GPR image for tunnels

5 結(jié)論

地質(zhì)雷達在短距離內(nèi)識別隧道掌子面前方的斷層破碎帶、溶洞等地質(zhì)體上高效、快捷、準(zhǔn)確，本文不僅從理論模型上進行了驗證，采集的實測數(shù)據(jù)也說明了其可靠性。根據(jù)正演及應(yīng)用實例分析可以總結(jié)以下幾點:

(1)利用GprMax2D軟件與Fortran語言可以很好地實現(xiàn)隧道超前預(yù)報模型下的數(shù)值模擬，效果理想。

(2)分析正演成果圖可知，電磁波在識別破碎帶、溶洞的頂界面位置、發(fā)育方向及形狀上準(zhǔn)確，波相位會隨著破碎帶兩端介質(zhì)電性差異產(chǎn)生變化，綜合這2點可以很好地獲取掌子面前方反射波帶回來的圍巖信息，為探明隧道前方圍巖情況提供很好的理論支撐。

(3)采集的2段隧道超前預(yù)報數(shù)據(jù)分別驗證了水平發(fā)育方向破碎帶、溶洞下的波場特征，反射波形特征給解釋人員提供了很重要的參考資料，可指導(dǎo)施工的順利進行。

［1］徐濟川，黃少霞.大瑤山隧道的突泥涌水機制［J］.鐵道工程學(xué)報，1996，(2):83－89.(XU Ji-chuan，HUANG Shao-xia.Mechanism of Mud and Water Burst in Dayao Mountain Tunnel［J］.Journal of Railway Engineering Society，1996，(2):83－89.(in Chinese))

［2］ZANG C W，HUANG H W，ZHANG Z X.Forecasting the Strata Condition of a Long Road Tunnel by Using Fuzzy Synthetic Judgment［J］.International Journal of Rock Mechanics Sciences，2004，41(3):406－407.

［3］齊傳生，王洪勇.圓梁山隧道綜合超前地質(zhì)預(yù)報技術(shù)［J］.鐵道勘察，2004，(5):52－56.(QI Chuan-sheng，WANG Hong-yong.Technique of Comprehensive Geological Prediction in Yuanliangshan Tunnel［J］.Railway Investigation and Surveying，2004，(5):52－56.(in Chinese))

［4］曲海峰，劉志剛，朱合華.隧道信息化施工中綜合超前地質(zhì)預(yù)報技術(shù)［J］.巖石力學(xué)與工程學(xué)報，2006，25(6):1246－ 1250.(QU Hai-feng，LIU Zhi-gang，ZHU He-hua.Technique of Synthetic Geologic Prediction Ahead in Tunnel Informational Construction［J］.Chinese Journal of Rock Mechanics and Engineering，2006，25(6):1246－1250.(in Chinese))

［5］鄧世坤，王惠濂.探地雷達圖像的正演合成與偏移處理［J］.地球物理學(xué)報，1993，36(4):528－535.(DENG Shi-kun，WANG Hui-lian.Forward Synthesizing and Migration Processing for GPR Image［J］.Chinese Journal of Geophysics，1993，36(4):528－535.(in Chinese))

［6］吳寶杰.探地雷達二維正演模擬及其工程實例［D］.杭州:浙江大學(xué)，2007.(WU Bao-jie.Two-dimensional Forward Simulation of Ground Penetrating Radar and Engineering Examples［D］.Hangzhou:Zhejiang University，2007.(in Chinese))

［7］趙 宇.GPR及TSP在隧道超前地質(zhì)預(yù)報中的解譯標(biāo)志研究［D］.成都:成都理工大學(xué)，2009.(ZHAO Yu.Study on the Interpretation of GPR and TSPin Tunnel Geological Forecast［D］.Chengdu:Chengdu University of Technology，2009.(in Chinese))

［8］趙 峰，周 斌，武永勝.探地雷達在隧道襯砌空洞檢測的正演模擬應(yīng)用研究［J］.鐵道建筑，2012，(8):99－103.(ZHAO Feng，ZHOU Bin，WU Yong-sheng.Study on Application of Forward modeling of Inspecting Void in Tunnel Lining by Ground Probing Radar［J］.Railway Engineering，2012，(8):99－ 103.(in Chinese))

［9］胡 俊，俞祁浩，游艷輝，等.探地雷達在多年凍土區(qū)正演模型研究及應(yīng)用［J］.物探與化探，2012，36(3):457－460.(HU Jun，YU Qi-hao，YOU Yan-hui，et al.Study and Application of Ground Penetration Radar to the Forward Modeling of Permafrost Area［J］.Geophysical＆Geochemical Exploration，2012，36(3):457－460.(in Chinese))

［10］孫忠輝，劉金坤，張新平，等.基于GprMax的隧道襯砌地質(zhì)雷達檢測正演模擬與實測數(shù)據(jù)分析［J］.工程地球物理學(xué)報，2013，10(5):730－735.(SUN Zhonghui，LIU Jin-kun，ZHANG Xin-ping，et al.The Tunnel Lining Detection Forward Numeral Simulation and Measured Data Analysis Based on GprMax［J］.Chinese Journal of Engineering Geophysics，2013，10(5):730－735.(in Chinese))

［11］周奇才，李炳杰，鄭宇軒，等.基于GPRMax2D的探地雷達圖像正演模擬［J］.工程地球物理學(xué)報，2008，5(4):396－ 399.(ZHOU Qi-cai，LI Bing-jie，ZHENG Yu-xuan，et al.Forward Simulation of GPR Image Based on GPRMax2D［J］.Chinese Journal of Engineering Geophysics，2008，5(4):396－399.(in Chinese))

［12］YEEK S.Numerical Solution of Initial Boundary Value Problem Involving Maxell’s Equations in Isotropic Media［J］.IEEE Transactions on Antennas and Propagation，1966，14(3):302－307.

［13］喻振華，馮德山，戴前偉，等.復(fù)雜地電模型的探地雷達時域有限差分正演［J］.物探化探計算技術(shù)，2005，27(4):279－283.(YU Zhen-hua，F(xiàn)ENG De-shan，DAI Qian-wei，et al.FDTD Method Forward Simulation of the Complex Geoelectricity GPR Model［J］.Computing Techniques for Geophysical and Geochemical Exploration，2005，27(4):279－283.(in Chinese))

［14］馮德山，戴前偉，何繼善，等.探地雷達GPR正演模擬的時域有限差分實現(xiàn)［J］.地球物理學(xué)進展，2006，21(2):630－636.(FENG De-shan，DAI Qian-wei，HE Ji-shan，et al.Finite Difference Time Domain Method of GPR Forward Simulation［J］.Progress in Geophysics，2006，21(2):630－636.(in Chinese))

［15］何兵壽，張會星.地質(zhì)雷達正演中的頻散壓制和吸收邊界改進方法［J］.地質(zhì)與勘探，2000，36(3):59－63.(HE Bing-shou，ZHANG Hui-xing.The Suppression of Numerical Dispersion and Improvement of Absorbing Boundary Conditions in Forward Modeling of GPR［J］.Geology and Prospecting，2000，36(3):59－ 63.(in Chinese))

［16］BERENGER J P.A Perfectly Match Layer for the Absorption of Electromagnetic Waves［J］.Journal of Computation Physics，1994，114(2):185－200.

［17］HASTINGS F D，SCHNEIDER J B，BROSCHAT S L.Application of the Perfectly Matched Layer(PML)Absorbing Boundary Condition to Elastic Wave Propagation［J］.The Journal of the Acoustical Society of America，1996，100(5):3061－3069.

［18］凌 飛，肖宏躍，朱夏樂，等.基于黏彈性介質(zhì)的隧道地震二維正演模擬［J］.長江科學(xué)院院報，2015，32(5):121－ 126.(LING Fei，XIAO Hong-yue，ZHU Xia-le，et al.Two-dimensional Forward Seismic Modeling for Tunnel Based on Viscoelastic Medium［J］.Journal of Yangtze River Scientific Research Institute， 2015，32(5):121－126.(in Chinese))

［19］呂喬森，徐 穎，胡志強.提高隧道地震波地質(zhì)預(yù)報(TSP)精度方法研究［J］.長江科學(xué)院院報，2012，29(12):41－ 45.(LV Qiao-sen，XU Ying，HU Zhiqiang.Methods of Improving Tunnel Seismic Prediction Accuracy［J］.Journal of Yangtze River Scientific Research Institute，2012，29(12):41－45.(in Chinese))