引 言

隧道工程在施工過程中經常受到地質條件和地（dì）形地貌等因素的製約，尤（yóu）其是複雜的地質條件可能誘發的（de）多種地質災害，嚴（yán）重（chóng）地（dì）影（yǐng）響工程施工進度，增加建設成本（běn），甚至造成傷亡事故.常（cháng）見的地質災害有岩溶、斷層、岩爆、突水（shuǐ）、突泥、大變形、瓦斯等.隧道（dào）工程施工過程中（zhōng）進行超前地質預報，探明前方工程地質狀況，進而指導隧道施工的順利進行有著重要（yào）的作用(薛國強和李貅，2008).

目前我國用（yòng）於隧道地質超（chāo）前預報有多種（zhǒng）方法（fǎ），主要（yào）有隧道地震超（chāo）前預報係統(TSP)、探地雷達法(GPR)、陸地（dì）聲呐法、水（shuǐ）平聲波剖麵（miàn）法(HSP)、超前鑽孔法等(宋先海等，2006；薛國（guó）強和李貅，2008).不同的方法有自身的（de）適應性和缺陷.探（tàn）地雷達法具有對施工影響、預報效率高等方麵優勢，但是在（zài）預報的解譯和探測的深度方麵存在不足.本文分析了基於有限差分法GPRMAX正演模擬的理論基礎，討論了二維FDTD法基本參數的選取，對隧道前方的不良（liáng）地（dì）質體進行正（zhèng）演模擬得到反射圖譜，並（bìng）結合工（gōng）程實例，說（shuō）明該正演模（mó）擬能夠指導解譯工作.

1 GPRMAX正演模擬的理論基礎

宏觀來看，所有的電磁現象都可通過Maxwell方程組來描述，它由兩個旋度方程和兩個散度方程組成（chéng）的，電（diàn）磁波分（fèn）為TEM、TM、TE三種模式.探地雷達中采用（yòng）橫磁模（mó）(TM型電磁波)，TM波方程為(葛德彪和閆玉波，2011；郭立（lì）等，2012)：

 

式中：ε為介質的介電常數（shù）；μ為磁導率；σ為電（diàn）導率；σm為等效磁阻率，引入等效磁阻（zǔ）率其目的主要（yào）在於使方程具有對稱（chēng）性，在大（dà）多數電磁場問（wèn）題中計算（suàn）空間內不包含磁性介（jiè）質，則有σm=0.

 

從上式可以看出，TM波隻有Z方向電場強度（dù）Ez、X、Y方向的磁場強度Hx、Hy三個分量（liàng）.運用Yee氏網（wǎng）格模型，利用中心差分代替對時間、空間坐標的微分，把連（lián）續變量離散（sàn）化，推導出探地雷達正（zhèng）演模擬方（fāng）程(李靜等，2011；曾昭發等，2010)：

 

式中：

Δt為時間步長；Δx、Δy為空間步長（zhǎng），ε(i，j)、σ(i，j)、μ(i，j)為(i，j)節點處的介電常數、電導率和磁導率.2 二維（wéi）FDTD法基本參數確定

2.1 空間步長的確定

用FDTD法對Maxwell方程（chéng）進行數值計算（suàn）時（shí），在時域（yù）有限差分網格中（zhōng），數值模擬的傳播速（sù）度將隨（suí）頻（pín）率改變，導致非物理（lǐ）因素（sù）引起的脈衝波形（xíng）畸變、人為的各向異性及虛假的折射（shè）現象，即出現數值色散現象，從（cóng）而造成數值不穩定性（xìng）(郭立等，2012).

二維空間中TM波的數值色散方程為（wéi）

式中：kx、ky分別（bié）為波（bō）矢量沿x、y方向的分量，ω為角頻率，υ為被模（mó）擬的均勻介質中的光速.

當Δt、Δx、Δy均趨於（yú）零時，色散可以減小到（dào）任意程度，但是由於（yú）計算機空間及速度等因素，決定時間步長和空間步長不可能無限的小，故選擇合適的時間步長和（hé）空間步長不可避免.

根據文獻(曾昭發等，2010)，當時，相速度的最大（dà）誤差為1.3%，當時，為0.31%，網格減小一倍，相位誤差減小為四（sì）分之一.本文Δl取略（luè）小於下的合適數據，選取原則是盡（jìn）量降低Δl的有效小數數字長（zhǎng）度，例如本文模型1中（zhōng）=0.036514837 m，Δl取0.03 m.

2.2 時間步長的確定

在FDTD中（zhōng），通過按時間步推進計算空間內的變化規律，時間步長與空間步（bù）長之間必須滿足一定（dìng）條件，解的（de）結果才能穩定，否則差分網格不穩定，隨著時（shí）間步數的增加，計算結果（guǒ）也將無限製地增加.

二維空間中TM波的穩定（dìng）性條件為

 

式中：υ為計算空間中的最大速度.

考慮到計算（suàn）的簡便性，探地雷達中取：

 

則：

2.3 吸收邊界厚度確定

時域差分法是在電磁場全部（bù）空間建立Yee矢網格計算空間的，但（dàn）是由於計算機（jī）的存儲空間都是有限（xiàn）的，在無限大網格空間中計算電（diàn）磁場根本不（bú）可能，因此在計算中（zhōng）總是在某處把網格空間（jiān）截斷，使之成為有限空間(張鴻飛等，2009).網格截斷處不引起波的明顯反射，對向外（wài）傳播的波而言（yán）就像在大空間傳播一樣.因此就需要在截斷處設置某種邊界條件（jiàn），使傳（chuán）輸到截斷（duàn）處的波被邊界條件吸收而不產生發射，從而起到模擬無（wú）限（xiàn）空間的目的.

邊界吸收條件發展經曆從開始簡單（dān）的（de）插值（zhí）邊界條件，到後來廣泛采用的Mur邊（biān）界條件階段，以及近（jìn）二十年來發展的完全（quán）匹配層(PML)吸收邊界階段，隨邊界條件發展其吸收效果越來越（yuè）好（hǎo）(葛德彪等，2011).在PML邊界條件（jiàn）中，層的厚度通（tōng）常有3~9個網格，本文邊界層選（xuǎn）取20個網（wǎng）格單（dān）元（yuán）.

3 隧道超（chāo）前地質預報數值模擬實例（lì）

根（gēn）據隧道異常地質的特點，建立異常探地雷達探測正演模型，得到模擬的結果，為探地（dì）雷達進行超前地質（zhì）預（yù）報的解譯工作提供參考（kǎo）.

3.1 正前方含水溶洞模型（xíng）

隧道掌子麵正前方30 m處有一圓形含（hán）水溶洞（dòng），直徑為0.5 m，如圖 1a所示.取雷達頻率為100 MHz，隧道介質的介電常數ε1=7.5，電導率σ1=0.00001 S/m；圓形溶洞中含水介（jiè）質的（de）介電常數ε2=81，電導率σ2=0.005 S/m；取網格的空間步長Δl=0.03 m，時間步長Δt=7.0711×10-11 s，時間窗口（kǒu）深度894 ns，迭代次（cì）數12635次，共模擬87道雷達波，模擬區域大小9 m×50 m，邊界層厚（hòu）度取0.6 m.

通過模擬，得（dé）到雷達波形圖如（rú）圖 1b和（hé）c所示.其（qí）中b圖為未去除（chú）直達波的反射剖麵（miàn）波形圖，可（kě）以看出，由於直達波的影響，直達波信號幾乎完全淹沒（méi）了溶洞的回波信號（hào）.圖 1c為經過零點平均、直達波置零方法去除直達波後的反射剖麵（miàn）波形圖，可以看出，從（cóng）30 m處（chù）開始，雷達波多（duō）次出現強（qiáng）烈反射，其後能量衰減明顯.可推測該處地質異常，圖形類似雙曲線狀，應為塊狀(或點狀)異常，雙曲線最高點應為異常點位（wèi）置.

另外本文將模擬的數據利用Matlab轉化成Mala格式，利用REFLEXW軟件進行處理，實（shí）現數（shù）據處理工具之間的互補，如圖 2所示（shì）為（wéi）模擬（nǐ）數據利用REFLEXW軟件（jiàn）顯示的波形圖，與圖 1b比較（jiào）可（kě）知，Matlab中顯示的（de）波形圖與REFLEXW下顯示的（de）結果一致，但是REFLEXW顯示的效果受到直達波的（de）影響較小（xiǎo），其REFLEXW下探地雷達常規處理方法也（yě）較成熟和專業化.

 

3.2 前（qián）方邊界溶（róng）洞一（yī）角模型

隧道掌（zhǎng）子麵前方30 m處有圓形含水溶洞一角(半圓形)，直徑為0.5 m，如圖 3a所示.共設置模擬88道雷達波，其（qí）餘參數選取同3.1中模型.

 

通過模擬，得到雷達波形圖，如圖 3b和c所示.可以看出，雷達（dá）波反射波形圖（tú）與3.1中模型（xíng）類似（sì），但由於溶洞是（shì）在邊界處，故（gù）其反射波形是雙曲線的一部分，凸點處就是溶洞位置.

3.3 斷層模型

隧道掌子麵前方30 m處（chù）有一（yī）含水斷層，斷層水平方向厚度為10 cm(即（jí）兩（liǎng）斷層垂直厚度約為7.1 cm)，傾角45°，如圖 4a所示.共設置模擬87道雷達波，其餘參（cān）數選取同3.1中模型.

 

去除直達波影響後得到如圖 4b所（suǒ）示的波形圖（tú），觀（guān）察可知，45°斷層（céng）反射波形為弧線形，反射能量主要集中在斷層位置處（chù），振蕩較少，而對於溶洞模型來說，其後為多次反射.

3.4 淤泥模型

隧道掌子（zǐ）麵（miàn）前方30 m處有一淤泥地帶（dài），設為一方形構造（zào），傾角45°，邊（biān）長1 m，如圖 5a所示.共設置模擬87道雷（léi）達波，其餘參數選取同3.1中模型.

 

去除直達（dá）波影響後得到（dào）如圖 5b所示的波形圖，可以看出，在（zài）淤泥處也形成多次連續明顯的強反射，波形為雙曲線，但是相（xiàng）對（duì）於3.1中模（mó）型，該模型多次反射間（jiān）距離明顯減小，信號衰減迅速.通（tōng）過（guò）兩者幾何模型比較，在電磁（cí）參數相同的情況下，與異常體的幾何形式及傾角有關係，結合3.3中（zhōng）斷（duàn）層模型，推斷出水平方向傾角越小（xiǎo），異常體處能量衰減速度越快.

4 工程應用實例

如（rú）圖 6為某地使用瑞（ruì）典Mala探地雷（léi）達在隧道掌子麵進行超（chāo）前地質預報時經數據處理後得（dé）到的GPR圖像，可以看出圖中圓（yuán）圈處的波形帶有明顯的雙曲線特征，且（qiě）隨深度(向右)增加，信號振蕩並逐漸減弱，判定為溶洞，後經工程施工證實該處判斷正確，如圖 7所示為該處現場圖像.

如圖 8所示（shì），為另一處隧道掌子麵超前地質預報GPR圖像（xiàng），可以看出（chū）圖中圓圈處呈現雙曲線波形，且迅速衰減，該處電磁信號頻率（lǜ）較低，振幅強，存在大量的多次振蕩，初步斷定為含水量較大的岩溶或淤泥，因此給出施工在該處需注意（yì）安全（quán）的建議，後經（jīng）證實該處為岩溶.

如圖 9所示，亦為隧道掌子（zǐ）麵超前地質（zhì）預報GPR圖像，圖中（zhōng）兩條黑線（xiàn）為近（jìn）似弧線線形，震蕩較（jiào）少，預估為裂隙發育，圖中裂隙右側區（qū）域信號雜亂，多處（chù）有雙曲線狀波形，推測為破碎區域，局部富水（shuǐ），後（hòu）亦已得到（dào）驗證.

上述通過對現場應（yīng）用實例（lì）的分（fèn）析，印證了（le）本文的模擬（nǐ）結果.反之（zhī），可以通過對探地雷達隧道超前預（yù）報的正演（yǎn）模（mó）擬的研究，指導（dǎo）工程的解譯工作.但（dàn）是需要注意的是由於工（gōng）程實際情況複雜，介質為非均勻物質，其物性參數存在或多或少的差異等原因，因此造成工程中探地雷達（dá）圖像顯得雜亂，這就需要工程（chéng）人員在解譯過程中仔細觀察圖像，從中找出圖（tú）像中的特征進行解譯（yì）.

5 結 語

5.1     通過分析GPRMAX正演模擬的理論基礎，利用中心差分（fèn）代替微分，推導出探地雷達正演模擬方程.

 

5.2     從降低數值色散現象，避免出現（xiàn）數值（zhí）不穩定性出發，依據（jù）數值色散方程確定出FDTD法空間步長的基本參數；依據解的穩定性條件，確定基本參（cān）數時間步長；以及吸收邊界層厚度（dù）的確定.

 

5.3     采用GPRMAX程（chéng）序對正前方含水溶洞、隧道邊界（jiè）處溶洞（dòng）、斷層模型、淤泥模（mó）型進行了正（zhèng）演模擬，模擬結果表明：圓形含水溶洞產生雙曲線弧線波形，雙曲線最（zuì）高（gāo）點即為溶洞位置；斷層射（shè）波形為弧線形，反射能量（liàng）主要集中在（zài）斷層（céng）位（wèi）置處，振（zhèn）蕩較少；不規則的淤泥模型也形成雙（shuāng）曲線波形，但是信號（hào）衰減迅速.通過模擬以及工程驗證，便於工作人員依據圖（tú）形特征（zhēng）準確進行（háng）地質解譯.

 

5.4     利用Matlab將模擬的數據轉（zhuǎn）化成Mala格式，利用REFLEXW軟件進行數據處理，結果（guǒ）一致，從而實現數據處理工具之（zhī）間的（de）互補.