尹光輝， 馮雨寧， 張懷凱， 馮興樂

(1.長安大學(xué)　信息工程學(xué)院,西安　710064；2.北京大學(xué)　地球與空間科學(xué)學(xué)院,北京　100871)

?

基于GprMax軟件的道路路基空洞探地雷達正演模擬

尹光輝1， 馮雨寧2， 張懷凱1， 馮興樂1

(1.長安大學(xué)信息工程學(xué)院,西安710064；2.北京大學(xué)地球與空間科學(xué)學(xué)院,北京100871)

探地雷達具有高效、無損、分辨率高、使用靈活的特點，在道路檢測方面的應(yīng)用越來越廣泛。根據(jù)時域有限差分(FDTD)算法，建立理想條件下道路路基充氣和充水空洞模型，用GprMax及Matlab編程分別對不同形狀、不同激勵源頻率、不同大小條件下空洞模型進行正演模擬，并對空洞模型正演模擬雷達圖像進行分析。模擬結(jié)果有助于對實際探測空洞識別和雷達圖像的解釋，為探地雷達在路基空洞檢測應(yīng)用提供參考。

探地雷達； 空洞； 正演模擬； FDTD算法； GprMax

0　引言

近年來，我國的道路空洞形成的塌陷事故愈演愈烈，給人民的生命財產(chǎn)安全造成了極大危害。與傳統(tǒng)鉆芯取樣道路檢測方法相比，由于探地雷達具有高效、無損、分辨率高的特點，且道路結(jié)構(gòu)層及土基和下方潛伏的充氣空洞、充水空洞之間的介電常數(shù)差異大，因此探地雷達技術(shù)適合應(yīng)用在道路空洞的檢測。在實際應(yīng)用過程中，缺乏探地雷達空洞異常圖像評判客觀標(biāo)準，道路空洞與探地雷達異常圖像的對應(yīng)關(guān)系缺乏系統(tǒng)地研究，且雷達資料解釋具有多解性和主觀經(jīng)驗性，導(dǎo)致同一個異常圖像不同的人得出不同的解釋結(jié)果，導(dǎo)致探地雷達技術(shù)在道路地下空洞探測中的應(yīng)用難以獲得理想效果[1]。因此，利用計算機開展道路路基空洞正演模擬，研究各種空洞模型的探地雷達探測圖像特征和規(guī)律，不必花費大量人力、物力和時間進行實地探測，就可以得到接近實地探測的效果。這對于認識探地雷達實測結(jié)果，識別空洞具有實際意義[2]。作者通過建立道路路基下充氣與充水空洞模型，用GprMax及Matlab編程實現(xiàn)道路路基空洞正演模擬，對空洞模型正演模擬圖像形態(tài)和特征進行分析，并通過實例驗證正演模擬的有效性。這里模擬結(jié)果有助于對實際探測空洞識別和探地雷達圖像地解釋，提高檢測結(jié)果的可靠性和檢測評判標(biāo)準的統(tǒng)一性，為探地雷達在路基空洞檢測應(yīng)用積累理論經(jīng)驗。

1　時域有限差分法正演的基本理論

探地雷達探測路基空洞中高頻電磁波以寬頻帶短脈沖形式，通過發(fā)射天線被定向送入地下，經(jīng)存在電性差異的地下地層或空洞反射后返回地面，由接收天線所接收。高頻電磁波在介質(zhì)中傳播時，其傳播路徑、電磁場強度與波形隨待測介質(zhì)的電性特征與幾何形態(tài)而變化[4]。因此，通過采集時域波形、并進行處理和分析，可確定地下分界面或地質(zhì)體的空間位置及結(jié)構(gòu)。

為了模擬探地雷達探測時的電磁波在地下介質(zhì)層及空洞中傳播，K.S.Yee提出一種計算電磁波數(shù)值的時域有限差分法，基于Yee[3]網(wǎng)格空間離散方式，將地下待探測區(qū)域劃分為由一定數(shù)量的空間網(wǎng)格組成的仿真空間。在FDTD離散中，Yee網(wǎng)格空間中電場和磁場各節(jié)點的空間排布見圖1。每個磁場分量由四個電場分量環(huán)繞，每個電場分量由四個磁場分量環(huán)繞。然后用有限差分式代替Maxwell方程時域場旋度方程中的微分式，在設(shè)置初始的場值和邊界條件后，通過逐步遞進的方法求得各個時刻的空間電磁波電磁場的分布數(shù)值，就獲得了電磁波在整個模型仿真空間的傳播模擬結(jié)果。

圖1　FDTD算法空間Yee網(wǎng)格Fig.1　Yee space grid of FDTD algorithm

由于FDTD算法是用有限差分式代替Maxwell方程組，離散后的差分方程組的解必須收斂穩(wěn)定。在FDTD中，時間步長Δt和空間網(wǎng)格步長Δx、Δy、Δz不是相互獨立的，它們的取值滿足一定的關(guān)系，才能避免數(shù)值結(jié)果的不穩(wěn)定[5-8]，解得穩(wěn)定性條件如式(1)所示。

(1)

用差分法對麥克斯韋方程進行數(shù)值計算將引起波的色散，即在FDTD網(wǎng)格中波的傳播速度將隨波長而改變。這種色散將導(dǎo)致非物理原因引起的脈沖波形畸變，使得時域數(shù)值計算產(chǎn)生誤差。因此，空間網(wǎng)格的幾何長度要滿足數(shù)值色散條件：

(2)

2　GprMax正演模擬

GprMax是基于FDTD算法的探地雷達正演模擬軟件，可以用于各向同性均勻介質(zhì)中電磁波的傳播規(guī)律以及電磁波與目標(biāo)體的相互作用，通過編寫正演模擬模型的輸入文件程序，保存后綴為*.in的文件，輸入程序文件路徑進行仿真。計算機進行數(shù)值計算后得到關(guān)于目標(biāo)體的探地雷達正演模擬后綴為*.out和*.geo的二進制文件。通過編寫Matlab程序讀入二進制文件數(shù)據(jù)，處理后綴*.geo文件獲得模型結(jié)構(gòu)圖，處理后綴*.out文件生成探地雷達正演模擬圖像(圖2)。

圖2　GprMax正演模擬流程圖Fig.2　Flow of GprMax forward simulation

2.1道路結(jié)構(gòu)模型參數(shù)

為了建立地下待探測區(qū)域的物理模型，需要設(shè)置相對介電常數(shù)、電導(dǎo)率參數(shù)。根據(jù)實際道路結(jié)構(gòu)，把道路模型定為三層：第一層是道路面層，由厚度為10 cm的瀝青材料構(gòu)成，瀝青材料相對介電常數(shù)=4，電導(dǎo)率σ=0.005 S/m；第二層為道路基層，由厚度為15 cm混凝土材料構(gòu)成，相對介電常數(shù)εr=9，電導(dǎo)率σ=0.05 S/m；第三層為土基層，由厚度為20 cm壓實處理后的土壤構(gòu)成并向下延伸到2 m，相對介電常數(shù)=12，電導(dǎo)率σ=0.1 S/m。水的相對介電常數(shù)εr=81，電導(dǎo)率σ=0.03 S/m。

3　道路路基空洞模型正演模擬

3.1不同形狀充氣及充水空洞正演模擬

根據(jù)道路結(jié)構(gòu)模型，分別構(gòu)造位于土基層的矩形、圓形空洞(圖3)，空洞中心點深度距離地表為0.6 m。矩形空洞左下頂點坐標(biāo)(1,1.2)，右上頂點坐標(biāo)(1.5,1.6); 圓形空洞圓心坐標(biāo)(1.25,1.4)，直徑為0.4 m;結(jié)構(gòu)激勵源為300 MHz的ricker子波，其他參數(shù)如表1所示。正演模擬結(jié)果如圖4所示。

圖3　不同形狀類型空洞結(jié)構(gòu)圖Fig.3　Structure of cavity in different shapes(a)矩形空洞; (b)圓形空洞

表1　FDTD正演模擬參數(shù)

對比分析不同形狀空洞正演模擬圖像，矩形、圓形充氣空洞正演模擬結(jié)果均呈現(xiàn)雙曲線特征，曲線開口向下。雙曲線頂部反射振幅最強，兩端振幅較弱。地下空洞在地面的水平中心投影位置可以由雙曲線頂點確定。矩形空洞由于頂面水平方向具有較長延伸，頂點反射較為明顯，雙曲線有水平弧度?？偟膩碚f圓形、矩形正演模擬圖像形成的雙曲線較為近似，很難通過雙曲線形狀估計地下空洞的大致形狀。

提取矩形充氣及充水空洞中心測線電磁波幅值數(shù)據(jù)得到電磁波波形圖(圖5、圖6)。發(fā)射天線電磁波從空氣射入面層、基層、路基時，介電常數(shù)依次增大，在每層分界面電磁波反射系數(shù)都為負值，因此在0～1 500測點數(shù)電磁波發(fā)生四次反向，充水型矩形空洞介電常數(shù)為81，遠大于路基層介電常數(shù)，電磁波反射系數(shù)為負值，中心測線處電磁波進入充水型空洞發(fā)生反向，充氣型空洞反射系數(shù)為正值，電磁波不發(fā)生反向。對比相同形狀但是內(nèi)部填充介質(zhì)不同的充水空洞，充水類型空洞正演模擬結(jié)果形成兩條雙曲線波形。以圓形充水空洞為例，電磁波在水中的傳播速度0.033 m/ns，傳播時長約為25 ns，兩者相乘除以2即為空洞的直徑0.412 5 m與模型參數(shù)設(shè)置的0.4 m接近。因此，實際探地雷達探測過程中，可以根據(jù)雙曲線的條數(shù)確定空洞內(nèi)部填充介質(zhì)是空氣還是水，充水型空洞可以根據(jù)兩條雙曲線的間隔時間，計算空洞的垂直距離。

3.2不同激勵源中心頻率正演模擬

在不同中心頻率正演模擬中，道路及空洞模型與圖2的圓形空洞模型一致，改變ricker子激勵源中心頻率，依次為200 MHz、300 MHz、400 MHz、500 MHz、600 MHz，研究不同中心頻率激勵源條件下雷達波成像規(guī)律，不同頻率下正演模擬圖像如圖7所示。

對比不同頻率正演模擬圖像，200 MHz圖像分辨率最低，隨著激勵源中心頻率的增大，雷達波正演模擬圖像分辨率逐漸提升。其中以400 MHz和600 MHz的雷達波正演模擬圖像分辨率最好，但不是頻率越高正演模擬圖像分辨率越好，500 MHz激勵源正演模擬圖像中出現(xiàn)許多雜波。

3.3不同大小圓形空洞正演模擬

不同大小圓形空洞頂面埋深一致，半徑分別為0.2 m、0.4 m、0.6 m，ricker激勵源頻率為400 MHz，其他參數(shù)與前面一致。

圖8中的正演模擬結(jié)果表明，隨著圓形空洞半徑的增大，正演模擬圖像生成的雙曲線雷達波圖像長度增大，曲率半徑增大。但是雙曲線延長線水平位置坐標(biāo)與實際空洞位置并不符合，不能通過雙曲線長度獲得空洞水平距離。

圖4　不同填充介質(zhì)正演模擬圖Fig.4　Simulation of cavity in different medium(a)矩形充氣空洞; (b)矩形充水空洞; (c)圓形充氣空洞; (d)圓形充水空洞

圖5　矩形充氣空洞中心測線波形圖Fig.5　Center line waveform of rectangle aeration cavity

圖6　矩形充水空洞中心測線波形圖Fig.6　Center line waveform of rectangle watery hole

圖7　不同頻率ricker激勵源正演模擬圖Fig.7　Simulation of different frequencies ricker excitation(a)頻率200 MHz；(b)頻率300 MHz；(c)頻率400 MHz；(d)頻率500 MHz；(e)頻率600 MHz

圖8　不同半徑圓形充氣空洞正演模擬圖Fig.8　Simulation of aeration circular cavity aeration in different radius(a)半徑0.2 m；(b)半徑0.4 m；(c)半徑0.6 m

4　探測實例正演模擬驗證

利用天毅達無線探地雷達對某處瀝青路面進行檢測，探測結(jié)果如圖9所示。在雷達增益為110倍時，實時偽彩圖顯示清晰雙曲線，與正演模擬結(jié)果呈現(xiàn)一致性，初步推斷為空洞。經(jīng)過探地雷達數(shù)據(jù)處理，顯示地下埋深約為0.5 m，實際驗證為一豎直排水口。根據(jù)探地雷達實測資料以及參數(shù)設(shè)置對此排水口進行正演模擬，構(gòu)造0.6 m×1.3 m×0.85 m三維空間，空間網(wǎng)格大小為0.01 m×0.01 m×0.01 m，路面介質(zhì)材料為瀝青，介電常數(shù)為5，電導(dǎo)率為0.01 S/m。構(gòu)造埋深約為0.5 m的矩形空洞，矩形空洞三維坐標(biāo)(0.2, 0.55, 0.1)(0.4, 0.75, 0.3)，收發(fā)天線坐標(biāo)(0.3, 0.125, 0.55)(0.3, 0.375, 0.55)，天線沿Y軸步進距離為0.04 m，測線道數(shù)21條，ricker激勵源發(fā)射頻率為900 MHz，時窗15 ns。正演模擬圖像如圖10所示。對比正演模擬圖像與實測圖像，路面結(jié)構(gòu)層分層時間都在3 ns，雙曲線頂點時間都在6 ns，圖像吻合度較好，驗證了基于時域有限差分法的GprMax結(jié)合Matlab進行探地雷達探測空洞正演模擬結(jié)果的有效性。

圖9　探地雷達實測圖像Fig.9　Image of GPR detection

圖10　實例正演模擬圖像Fig.10　Simulation of instance cavity

5　結(jié)論

基于時域有限差分法，通過GprMax及Matlab編程實現(xiàn)對道路路基空洞模型的正演模擬，模擬結(jié)果表明：

1)不同形狀空洞正演模擬結(jié)果雙曲線類似，難以通過雙曲線反推空洞形狀。可根據(jù)雙曲線頂點計算空洞頂面埋深，根據(jù)兩條雙曲線相差走時計算充水型空洞垂直距離，根據(jù)雙曲線條數(shù)推斷空洞充水還是充氣類型。

2)不同激勵源中心頻率，影響探地雷達圖像分辨率，一定頻率范圍內(nèi)中心頻率越高，探地雷達圖像分辨率越高。

3)不同大小規(guī)?？斩葱纬傻碾p曲線曲率不同、延長度不同，空洞規(guī)模越大，雙曲線延長度越大，曲率越大。

通過探地雷達實測圖像與正演模擬結(jié)果對比，驗證了正演模擬結(jié)果的有效性，可用于指導(dǎo)探地雷達空洞探測實例識別，為建立探地雷達空洞模型數(shù)據(jù)庫提供數(shù)據(jù)資料。

[1]葛廣志.北京市道路地下空洞探地雷達探測正演模擬研究[D].北京:中國地質(zhì)大學(xué),2010.

GE G Z.Forward simulation of underground cavity in Beijing by GPR[D].Beijing:China University of Geosciences,2010.(In Chinese)

[2]郭成超,王復(fù)明.探地雷達電磁波正演模擬研究[J].公路交通科技,2008,25(8):37-41.

GUO C C,WANG F M.Study on forward modeling of GPR electromagnetic wave propagation[J].Journal of Highway and Transportation Research and Development,2008,25(8):37-41.(In Chinese)

[3]彭湘佳.探地雷達在北京市道路下方病害檢測中的研究[D].北京:中國地質(zhì)大學(xué),2011.

PENG X J.Research on the disease below road in Beijing by GPR[D].Beijing:China University of Geosciences,2010.(In Chinese)

[4]馮彥謙,王銀,劉四新.基于FDTD的探地雷達數(shù)值仿真與成像研究[J].鐵道工程學(xué)報,2009,26(6):17-20.

FENG Y Q,WANG Y,LIU S X.Research on the numerical simulation and imaging of GPR based on the FDTD[J].Journal of Railway Engineering Society,2009,26(6):17-20.(In Chinese)

[5]周奇才,李炳杰,鄭宇軒,等.基于GPRMax2D的探地雷達圖像正演模擬[J].工程地球物理學(xué)報,2008,5(4):396-399.

ZHOU Q C,LI Y J,ZHENG Y X,et al.Forward simulation of GPR image based on GPRMax2D[J].Chinese Journal of Engineering Geophysics,2008,5(4):396-399.(In Chinese)

[6]謝雄耀,于超,趙永輝,等.山嶺隧道空洞災(zāi)害雷達探測圖像影響因素的正演模擬研究[J].工程地質(zhì)學(xué)報,2010,18(1):81-87.

XIE X Y,YU C, ZHAO Y H,et al.The forward simulation research of cavity radar image influence factors in mountain tunnel[J].Journal of Engineering Geology,2010,18(1):81-87.(In Chinese)

[7]徐波,蒙愛軍,吳鳳昌.利用數(shù)值模擬仿真探地雷達地下障礙物探測[J].路基工程,2010,27(5):170-171.

XU B,MENG A J,WU F C.Underground obstacle detection of GPR by numerical simulation[J].Subgrade Engineering,2010,27(5):170-171.(In Chinese)

[8]A.GIANNOPOULOS.Modeling ground penetrating radar by GprMax[J].Construction and Building Materials,2005,19(10):755-762.

Forward simulation of ground penetration radar based on the GprMax for the roadbed cavity

YIN Guang-hui1, FENG Yu-ning2, ZHANG Huai-kai, FENG Xing-le

(1.School of Information Engineering, Chang 'an University, Xi'an710064, China;(2.School of Earth and Space Sciences, Peking University, Beijing100871,China)

Ground penetrating radar has been widely used in road detection due to the features of efficient, non-destructive, high resolution and flexible. Based on the finite different time domain algorithm (FDTD), the ideal models of roadbed cavity filled with air and water are built, GprMax and Matlab programming are used to realize forward simulation which aimed at different shapes, different center frequencies and different sizes, and then the cavity modeling simulation radar images are analyzed. Simulation results contribute to the detection of cavity and the interpretation of the radar image and references are provided for the ground penetrating radar application in cavity detection.

ground penetration radar; cavity; forward simulation; FDTD algorithm; GprMax

2015-06-15改回日期：2015-09-05

國家自然科學(xué)基金(41404095)

尹光輝(1990-)，男，碩士，主要研究方向為探地雷達檢測，E-mail:412245505@qq.com。

1001-1749(2016)04-0480-07

P 631.3

A

10.3969/j.issn.1001-1749.2016.04.07