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

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

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

?

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

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

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

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

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

0　引言

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

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

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

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

圖1　FDTD算法空間Yee網格Fig.1　Yee space grid of FDTD algorithm

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

(1)

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

(2)

2　GprMax正演模擬

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

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

2.1道路結構模型參數

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

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

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

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

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

表1　FDTD正演模擬參數

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

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

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

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

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

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

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

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

圖4　不同填充介質正演模擬圖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　探測實例正演模擬驗證

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

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

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

5　結論

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

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

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

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

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

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

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

[2]郭成超,王復明.探地雷達電磁波正演模擬研究[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].北京:中國地質大學,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的探地雷達數值仿真與成像研究[J].鐵道工程學報,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].工程地球物理學報,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]謝雄耀,于超,趙永輝,等.山嶺隧道空洞災害雷達探測圖像影響因素的正演模擬研究[J].工程地質學報,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]徐波,蒙愛軍,吳鳳昌.利用數值模擬仿真探地雷達地下障礙物探測[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

國家自然科學基金(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