中交路橋科技是從事工程檢測監(jiān)測、城市安全監(jiān)測預(yù)警與評價、數(shù)字智能化研發(fā)為一體的復(fù)合型高新技術(shù)集團(tuán)企業(yè)。
隧道工程在施工過程中經(jīng)常受到地質(zhì)條件和地形地貌等因素的制約，尤其是復(fù)雜的地質(zhì)條件可能誘發(fā)的多種地質(zhì)災(zāi)害，嚴(yán)重地影響工程施工進(jìn)度，增加建設(shè)成本，甚至造成傷亡事故.常見的地質(zhì)災(zāi)害有巖溶、斷層、巖爆、突水、突泥、大變形、瓦斯等.隧道工程施工過程中進(jìn)行超前地質(zhì)預(yù)報，探明前方工程地質(zhì)狀況，進(jìn)而指導(dǎo)隧道施工的順利進(jìn)行有著重要的作用(薛國強(qiáng)和李貅，2008).
目前我國用于隧道地質(zhì)超前預(yù)報有多種方法，主要有隧道地震超前預(yù)報系統(tǒng)(TSP)、探地雷達(dá)法(GPR)、陸地聲吶法、水平聲波剖面法(HSP)、超前鉆孔法等(宋先海等，2006；薛國強(qiáng)和李貅，2008).不同的方法有自身的適應(yīng)性和缺陷.探地雷達(dá)法具有對施工影響、預(yù)報效率高等方面優(yōu)勢，但是在預(yù)報的解譯和探測的深度方面存在不足.本文分析了基于有限差分法GPRMAX正演模擬的理論基礎(chǔ)，討論了二維FDTD法基本參數(shù)的選取，對隧道前方的不良地質(zhì)體進(jìn)行正演模擬得到反射圖譜，并結(jié)合工程實例，說明該正演模擬能夠指導(dǎo)解譯工作.
1 GPRMAX正演模擬的理論基礎(chǔ)
宏觀來看，所有的電磁現(xiàn)象都可通過Maxwell方程組來描述，它由兩個旋度方程和兩個散度方程組成的，電磁波分為TEM、TM、TE三種模式.探地雷達(dá)中采用橫磁模(TM型電磁波)，TM波方程為(葛德彪和閆玉波，2011；郭立等，2012)：
式中：ε為介質(zhì)的介電常數(shù)；μ為磁導(dǎo)率；σ為電導(dǎo)率；σm為等效磁阻率，引入等效磁阻率其目的主要在于使方程具有對稱性，在大多數(shù)電磁場問題中計算空間內(nèi)不包含磁性介質(zhì)，則有σm=0.
從上式可以看出，TM波只有Z方向電場強(qiáng)度Ez、X、Y方向的磁場強(qiáng)度Hx、Hy三個分量.運用Yee氏網(wǎng)格模型，利用中心差分代替對時間、空間坐標(biāo)的微分，把連續(xù)變量離散化，推導(dǎo)出探地雷達(dá)正演模擬方程(李靜等，2011；曾昭發(fā)等，2010)：
2.1 空間步長的確定
用FDTD法對Maxwell方程進(jìn)行數(shù)值計算時，在時域有限差分網(wǎng)格中，數(shù)值模擬的傳播速度將隨頻率改變，導(dǎo)致非物理因素引起的脈沖波形畸變、人為的各向異性及虛假的折射現(xiàn)象，即出現(xiàn)數(shù)值色散現(xiàn)象，從而造成數(shù)值不穩(wěn)定性(郭立等，2012).
二維空間中TM波的數(shù)值色散方程為
當(dāng)Δt、Δx、Δy均趨于零時，色散可以減小到任意程度，但是由于計算機(jī)空間及速度等因素，決定時間步長和空間步長不可能無限的小，故選擇合適的時間步長和空間步長不可避免.
根據(jù)文獻(xiàn)(曾昭發(fā)等，2010)，當(dāng)時
，相速度的最大誤差為1.3%，當(dāng)時，為0.31%，網(wǎng)格減小一倍，相位誤差減小為四分之一.本文Δl取略小于
下的合適數(shù)據(jù)，選取原則是盡量降低Δl的有效小數(shù)數(shù)字長度，例如本文模型1中=0.036514837 m，Δl取0.03 m.
2.2 時間步長的確定
在FDTD中，通過按時間步推進(jìn)計算空間內(nèi)的變化規(guī)律，時間步長與空間步長之間必須滿足一定條件，解的結(jié)果才能穩(wěn)定，否則差分網(wǎng)格不穩(wěn)定，隨著時間步數(shù)的增加，計算結(jié)果也將無限制地增加.
二維空間中TM波的穩(wěn)定性條件為
考慮到計算的簡便性，探地雷達(dá)中?。?/p>
則：
時域差分法是在電磁場全部空間建立Yee矢網(wǎng)格計算空間的，但是由于計算機(jī)的存儲空間都是有限的，在無限大網(wǎng)格空間中計算電磁場根本不可能，因此在計算中總是在某處把網(wǎng)格空間截斷，使之成為有限空間(張鴻飛等，2009).網(wǎng)格截斷處不引起波的明顯反射，對向外傳播的波而言就像在大空間傳播一樣.因此就需要在截斷處設(shè)置某種邊界條件，使傳輸?shù)浇財嗵幍牟ū贿吔鐥l件吸收而不產(chǎn)生發(fā)射，從而起到模擬無限空間的目的.
邊界吸收條件發(fā)展經(jīng)歷從開始簡單的插值邊界條件，到后來廣泛采用的Mur邊界條件階段，以及近二十年來發(fā)展的完全匹配層(PML)吸收邊界階段，隨邊界條件發(fā)展其吸收效果越來越好(葛德彪等，2011).在PML邊界條件中，層的厚度通常有3~9個網(wǎng)格，本文邊界層選取20個網(wǎng)格單元.
3 隧道超前地質(zhì)預(yù)報數(shù)值模擬實例
根據(jù)隧道異常地質(zhì)的特點，建立異常探地雷達(dá)探測正演模型，得到模擬的結(jié)果，為探地雷達(dá)進(jìn)行超前地質(zhì)預(yù)報的解譯工作提供參考.
3.1 正前方含水溶洞模型
隧道掌子面正前方30 m處有一圓形含水溶洞，直徑為0.5 m，如圖 1a所示.取雷達(dá)頻率為100 MHz，隧道介質(zhì)的介電常數(shù)ε1=7.5，電導(dǎo)率σ1=0.00001 S/m；圓形溶洞中含水介質(zhì)的介電常數(shù)ε2=81，電導(dǎo)率σ2=0.005 S/m；取網(wǎng)格的空間步長Δl=0.03 m，時間步長Δt=7.0711×10-11 s，時間窗口深度894 ns，迭代次數(shù)12635次，共模擬87道雷達(dá)波，模擬區(qū)域大小9 m×50 m，邊界層厚度取0.6 m.
 |
圖 1 含水溶洞模型及波形圖Fig. 1 Water cave model and waveform figure |
通過模擬，得到雷達(dá)波形圖如圖 1b和c所示.其中b圖為未去除直達(dá)波的反射剖面波形圖，可以看出，由于直達(dá)波的影響，直達(dá)波信號幾乎完全淹沒了溶洞的回波信號.圖 1c為經(jīng)過零點平均、直達(dá)波置零方法去除直達(dá)波后的反射剖面波形圖，可以看出，從30 m處開始，雷達(dá)波多次出現(xiàn)強(qiáng)烈反射，其后能量衰減明顯.可推測該處地質(zhì)異常，圖形類似雙曲線狀，應(yīng)為塊狀(或點狀)異常，雙曲線最高點應(yīng)為異常點位置.
另外本文將模擬的數(shù)據(jù)利用Matlab轉(zhuǎn)化成Mala格式，利用REFLEXW軟件進(jìn)行處理，實現(xiàn)數(shù)據(jù)處理工具之間的互補(bǔ)，如圖 2所示為模擬數(shù)據(jù)利用REFLEXW軟件顯示的波形圖，與圖 1b比較可知，Matlab中顯示的波形圖與REFLEXW下顯示的結(jié)果一致，但是REFLEXW顯示的效果受到直達(dá)波的影響較小，其REFLEXW下探地雷達(dá)常規(guī)處理方法也較成熟和專業(yè)化.
 |
圖 2 REFLEXW中顯示的波形圖Fig. 2 Waveform figure in REFLEXW |
隧道掌子面前方30 m處有圓形含水溶洞一角(半圓形)，直徑為0.5 m，如圖 3a所示.共設(shè)置模擬88道雷達(dá)波，其余參數(shù)選取同3.1中模型.
 |
圖 3 邊界有含水溶洞模型Fig. 3 Border water cave model |
通過模擬，得到雷達(dá)波形圖，如圖 3b和c所示.可以看出，雷達(dá)波反射波形圖與3.1中模型類似，但由于溶洞是在邊界處，故其反射波形是雙曲線的一部分，凸點處就是溶洞位置.
3.3 斷層模型
隧道掌子面前方30 m處有一含水?dāng)鄬?，斷層水平方向厚度?0 cm(即兩斷層垂直厚度約為7.1 cm)，傾角45°，如圖 4a所示.共設(shè)置模擬87道雷達(dá)波，其余參數(shù)選取同3.1中模型.
 |
圖 4 斷層模型Fig. 4 Fault model |
去除直達(dá)波影響后得到如圖 4b所示的波形圖，觀察可知，45°斷層反射波形為弧線形，反射能量主要集中在斷層位置處，振蕩較少，而對于溶洞模型來說，其后為多次反射.
3.4 淤泥模型
隧道掌子面前方30 m處有一淤泥地帶，設(shè)為一方形構(gòu)造，傾角45°，邊長1 m，如圖 5a所示.共設(shè)置模擬87道雷達(dá)波，其余參數(shù)選取同3.1中模型.
 |
圖 5 淤泥模型Fig. 5 Silt model |
去除直達(dá)波影響后得到如圖 5b所示的波形圖，可以看出，在淤泥處也形成多次連續(xù)明顯的強(qiáng)反射，波形為雙曲線，但是相對于3.1中模型，該模型多次反射間距離明顯減小，信號衰減迅速.通過兩者幾何模型比較，在電磁參數(shù)相同的情況下，與異常體的幾何形式及傾角有關(guān)系，結(jié)合3.3中斷層模型，推斷出水平方向傾角越小，異常體處能量衰減速度越快.
4 工程應(yīng)用實例
如圖 6為某地使用瑞典Mala探地雷達(dá)在隧道掌子面進(jìn)行超前地質(zhì)預(yù)報時經(jīng)數(shù)據(jù)處理后得到的GPR圖像，可以看出圖中圓圈處的波形帶有明顯的雙曲線特征，且隨深度(向右)增加，信號振蕩并逐漸減弱，判定為溶洞，后經(jīng)工程施工證實該處判斷正確，如圖 7所示為該處現(xiàn)場圖像.
 |
圖 6 現(xiàn)場實測GPR圖像1Fig. 6 GPR image 1 in project |
 |
圖 7 工程現(xiàn)場Fig. 7 Engineering field |
如圖 8所示，為另一處隧道掌子面超前地質(zhì)預(yù)報GPR圖像，可以看出圖中圓圈處呈現(xiàn)雙曲線波形，且迅速衰減，該處電磁信號頻率較低，振幅強(qiáng)，存在大量的多次振蕩，初步斷定為含水量較大的巖溶或淤泥，因此給出施工在該處需注意安全的建議，后經(jīng)證實該處為巖溶.
 |
圖 8 現(xiàn)場實測GPR圖像2Fig. 8 GPR image 2 in project |
如圖 9所示，亦為隧道掌子面超前地質(zhì)預(yù)報GPR圖像，圖中兩條黑線為近似弧線線形，震蕩較少，預(yù)估為裂隙發(fā)育，圖中裂隙右側(cè)區(qū)域信號雜亂，多處有雙曲線狀波形，推測為破碎區(qū)域，局部富水，后亦已得到驗證.
 |
圖 9 現(xiàn)場實測GPR圖像3Fig. 9 GPR image 3 in project |
上述通過對現(xiàn)場應(yīng)用實例的分析，印證了本文的模擬結(jié)果.反之，可以通過對探地雷達(dá)隧道超前預(yù)報的正演模擬的研究，指導(dǎo)工程的解譯工作.但是需要注意的是由于工程實際情況復(fù)雜，介質(zhì)為非均勻物質(zhì)，其物性參數(shù)存在或多或少的差異等原因，因此造成工程中探地雷達(dá)圖像顯得雜亂，這就需要工程人員在解譯過程中仔細(xì)觀察圖像，從中找出圖像中的特征進(jìn)行解譯.
5 結(jié) 語
5.1 通過分析GPRMAX正演模擬的理論基礎(chǔ)，利用中心差分代替微分，推導(dǎo)出探地雷達(dá)正演模擬方程.5.2 從降低數(shù)值色散現(xiàn)象，避免出現(xiàn)數(shù)值不穩(wěn)定性出發(fā)，依據(jù)數(shù)值色散方程確定出FDTD法空間步長的基本參數(shù)；依據(jù)解的穩(wěn)定性條件，確定基本參數(shù)時間步長；以及吸收邊界層厚度的確定.
5.3 采用GPRMAX程序?qū)φ胺胶芏?、隧道邊界處溶洞、斷層模型、淤泥模型進(jìn)行了正演模擬，模擬結(jié)果表明：圓形含水溶洞產(chǎn)生雙曲線弧線波形，雙曲線最高點即為溶洞位置；斷層射波形為弧線形，反射能量主要集中在斷層位置處，振蕩較少；不規(guī)則的淤泥模型也形成雙曲線波形，但是信號衰減迅速.通過模擬以及工程驗證，便于工作人員依據(jù)圖形特征準(zhǔn)確進(jìn)行地質(zhì)解譯.
5.4 利用Matlab將模擬的數(shù)據(jù)轉(zhuǎn)化成Mala格式，利用REFLEXW軟件進(jìn)行數(shù)據(jù)處理，結(jié)果一致，從而實現(xiàn)數(shù)據(jù)處理工具之間的互補(bǔ).
