水體的探地雷達成像數(shù)值模擬與實驗比較研究

孫 偉 徐 青 藍朝楨 李建勝

(1.信息工程大學測繪學院，河南 鄭州 450052；2.蘇州市數(shù)字城市工程研究中心，江蘇 蘇州 215021)

引 言

探地雷達(GPR)是一種新型的地下目標探測手段[1]。它利用電磁波對地表的穿透能力，從地表上向地下發(fā)射脈沖電磁波，電磁波在地下阻抗特性發(fā)生變化的界面上發(fā)生反射，通過接收反射回波信號，根據(jù)其時延、形狀及頻譜特性等參數(shù)，可以解譯出目標深度和介質特性。在地下淺層目標的探測和識別方面，除了進行探測設備的研制和改進外，理論上的模擬計算也是這一領域的重要內容之一。GPR的數(shù)值模擬是分析探測問題，研究電磁波在介質中傳播規(guī)律的有效手段，對于提高探測的效果和解釋的準確性具有重要意義。此外，利用正演結果進行反演，可驗證反演算法的正確性，提高解釋精度。目前，國內外學者開展了大量的探地雷達數(shù)值模擬工作，取得了很大的進展[2-8]。時域有限差分方法(FDTD)已經成為主流算法被廣泛采用[9]，目前正朝著大規(guī)模并行計算方向以及快速算法發(fā)展以解決三維GPR數(shù)值模擬計算量大的問題[10]。數(shù)值模擬的結果需要與實測結果進行對比才能判斷數(shù)值模擬的正確性和可靠性，兩者的對比研究不僅在理論上而且在工程應用上都具有重要意義。然而令人遺憾的是，已經公開發(fā)表的文獻較少涉及兩者的對比研究。如何利用數(shù)值模擬結果更好地為提高GPR探測效果以及解釋的準確性服務變得越來越重要。下面從數(shù)值模擬和實驗室實測兩個方面入手來進行對比研究。

1. 算法原理

GPR技術的應用領域逐漸由低損耗介質向高損耗介質拓展，如堤壩、水體等，而介質的導電性會造成電磁波能量衰減，這是高頻GPR只能進行淺層探測的原因。在早期探地雷達的資料處理和解釋中,曾不考慮導電率的作用，隨著探測精度和要求的提高，介質的導電性越來越重要。含衰減的GPR電磁波傳播基本規(guī)律可用經典的Maxwell方程來描述

(1)

其中：E為電場強度，單位為V/m；H為磁場強度，單位為A/m；ε為介電常數(shù)，單位為F/m；μ為磁導率，單位為H/m；J是電流密度，單位為A/m2；Jm為磁流密度，單位為V/m2；σ為電導率，單位為S/m；σm為等效磁阻率，單位為Ω/m.

在二維剖面情況下，電磁波有TE和TM兩種模式，分別對應電場矢量平行和垂直于電磁波入射平面。其中TM模的Maxwell方程如下

(2)

采用FDTD來數(shù)值解上述方程組，在TM模式下，把電場安排在整數(shù)網格上，把磁場放置在半網格上，離散化方程為

(3)

上式中：Δx，Δy，Δt分別是x，y方向的網格步長以及時間步長，需滿足如下Courant穩(wěn)定性條件

(4)

其中,Vmax是計算域的最大波速。有限的計算空間需要增加吸收層才能模擬電磁波在無限空間傳播，采用8層完全匹配層(PML)[11]作為邊界條件來吸收外向電磁波。

電磁波激勵源采用超寬帶Ricker子波

A(t)=A0[1-2(πfmt)]e-(πfmt)2

(5)

式中：fm是超寬帶激勵源中心頻率；A0是脈沖輸入最大值。收發(fā)天線間距為0.

2. 實驗情況

圖1給出了GPR探測的實驗裝置，在一個57 cm×43 cm×34 cm的箱體內充滿水，箱體壁厚約5 mm，金屬或聚氯乙烯(PVC)管線插入水中。使用儀器為意大利IDS公司生產的RIS型探地雷達，天線中心頻率1.6 GHz.測線布置在箱體的外側面中心水平方向，GPR天線緊貼著箱壁沿著測線從右向左水平掃描。之所以采取側面掃描，而不采用水面朝下掃描的方式，主要是為了實驗的簡便以及GPR天線的定位。管線采用較有代表性的金屬管線和PVC管線兩種。由于天線頻率較高，水中電磁波的衰減較大，因此，探測深度有限，管線放置的位置距離前壁300 mm以內。

圖1 GPR探測實驗裝置

3. 結果分析與討論

實驗裝置中PVC管線的壁厚僅僅 1.5 mm，介電常數(shù)約為3，內外都充滿介電常數(shù)為81的水。為了建立薄壁PVC管線模型，必須用很小的網格尺寸。設置網格步長為0.2 mm，計算空間總網格數(shù)為2 900×2 100.由于單次計算量較大，為了減少總的計算量，計算過程中僅均勻設置44個測點，每個測點間距1 cm.另外，計算模型包含了箱體前壁，但不包含側壁和后壁，箱體材料的介電常數(shù)約為3.

3.1 淡水中的Ф50 mm金屬和Ф50 mm PVC管線

淡水中介質損耗相對較小，電導率設為0.5 mS/m，磁導率為0.金屬管線和PVC管線的直徑均為50 mm，兩種管線距離前壁為100 mm，兩者間距170 mm，他們的實驗布局和數(shù)值模擬模型如圖2所示。

圖2 實驗與數(shù)值模擬的管線相對位置

(a) 實測 (b) 數(shù)值模擬圖3 Ф50 mm金屬和Ф50 mm PVC管線成像比較

圖3給出了實驗探測的雷達剖面和數(shù)值模擬成像的對比。從圖中可以清晰辨別右邊較亮的雙曲線屬于金屬管線成像，而左邊較暗的雙曲線屬于PVC管線成像。左右兩張圖的雙曲線位置、形狀參數(shù)基本一致。因為水箱前壁會部分反射電磁波，反射下行的電磁波又被管線散射后上行進入雷達天線，所以可以看到對于一個管線的雷達剖面是由許多個上下排列、強弱分布的雙曲線共同組成。實測成像和數(shù)值模擬結果很好地證實了這個現(xiàn)象。在實際應用當中，類似的情況同樣會出現(xiàn)，下行的電磁波在被探測物體與地表間來回反射，形成一系列明暗的雙曲線。然而，該現(xiàn)象會對圖像的解釋產生干擾。由于GPR本身的信噪比以及環(huán)境因素，GPR剖面經過一系列后處理以后，多次反射產生的后續(xù)雙曲線可能要強于第一個雙曲線，或者金屬的后續(xù)雙曲線可能強于介質的第一個雙曲線，正如圖3的PVC管線成像所示。這對圖像的解釋能力提出了很高的要求，如何從雷達剖面里區(qū)分出金屬管線和非金屬管線以及它們的相對位置需要有豐富的經驗。但要從圖3左圖中識別出管線類型以及每種管線的數(shù)量是困難的，數(shù)值模擬可以在一定程度上幫助解釋人員完成這項任務。

3.2 淡水中的Ф75 mm PVC和Ф50 mm PVC管線

相同材料、不同直徑的管線的識別在實際應用中很重要，這需要探地雷達有足夠的橫向分辨率。圖4給出了直徑分別為75 mm和50 mm PVC管線的實測和數(shù)值模擬雷達剖面。它們與前壁距離以及兩者間距與3.1節(jié)相同。從圖4中兩圖對比可以看到：對于不同的管徑，雙曲線的形狀參數(shù)基本能區(qū)分出來。右圖的雙曲線比左圖復雜得多，這主要是邊界條件殘留的反射造成的。和圖3比較，可以看出相同材料的管線比不同材料的管線更容易從雷達剖面中識別出來，金屬比非金屬容易識別。

圖5是將兩個PVC管線到前壁的距離增大到200 mm后的雷達剖面。顯然，距離越大，管線越不容易被探測到。由于水對電磁波的衰減，實際的GPR信噪比較低，雙曲線變得模糊且不完整，已經較難從左圖的實測圖像中根據(jù)公式計算出管徑。但是，右圖管線的雙曲線依然非常清晰，這是因為數(shù)值模擬沒有噪音以及其他各種影響，盡管信號很微弱，但信噪比足夠使得雙曲線清晰、完整。從這也可以看出：盡管數(shù)值模擬能模擬介質衰減等各種情況，但缺乏GPR實際工作條件下的背景噪音，使得在有耗介質存在時的數(shù)值模擬結果與實測結果差別較大。

(a)實測 (b)數(shù)值模擬圖4 Ф50 mm和Ф75 mm PVC管線成像比較

(a)實測 (b)數(shù)值模擬圖5 Ф50 mm和Ф75 mm PVC管線成像比較

3.3 鹽水中的Ф50 mm金屬+Ф50 mm PVC 管線

GPR的探測深度與介質的導電率直接相關，3.2節(jié)的結果是在淡水條件下得到的，電導率很小?，F(xiàn)在水中加入160 g鹽，按照1 μS/cm=0.55～0.75 mg/l計算，電導率σ在0.2～0.3 S/m之間，取中間值σ=0.25 S/m作為電導率計算參數(shù)。

圖6給出了實測所得的雷達剖面，金屬和PVC管線直徑均為50 mm，對應管線到前壁的距離(探測深度)分別從40 mm增大到150 mm.從圖可以看出：金屬管線的最大探測深度明顯大于PVC管線，金屬管線在120 mm處雷達剖面依然有較強的雙曲線，而PVC管線在80 mm以上就無法成像。在較大的電導率的介質中，只有當管線在淺層(如圖6(a))時才能在雷達剖面中觀察到如前面所述的多個雙曲線上下排列的現(xiàn)象；在當管線埋藏較深時，多次反射的信號基本上被介質吸收，這就是圖6(b)～(f)一個管線在雷達剖面上對應一個雙曲線的原因。

(a) h=40 mm (b) 60 mm

(c) h=80 mm (d) 100 mm

(e) h=120 mm (f) 150 mm圖6 鹽水中不同探測深度對應的實測雷達剖面圖

圖7是與圖6一一對應的數(shù)值模擬雷達剖面圖。根據(jù)前面的討論，由于沒有信噪比的限制，在管線與前壁任何間距下，數(shù)值模擬的雷達剖面都可以清晰辨別兩個管線對應的兩條雙曲線。只是隨著管線-前壁的距離增加，接收的電場強度減弱，雷達剖面上的雙曲線變得暗淡。通過以下兩種方法可以使數(shù)值模擬與實際相符：第一種方法是在數(shù)值模擬中對接收的電場強度設置某一個與實際噪音相關的臨界值，當電場強度低于該臨界值時表示電場強度太弱被探測物無法成像； 第二種方法是在激勵Ricker子波中加入噪音。但受到FDTD數(shù)值色散的影響，第二種方法可能會受到限制。

(a) h=40 mm (b) 60 mm

(c) h=80 mm (d) 100 mm

(e) h=120 mm (f) 150 mm圖7 鹽水中不同探測深度對應的數(shù)值模擬圖

4. 結 論

開展了中心頻率為1.6 GHz的GPR水箱管線成像研究，并建立了二維的數(shù)值模擬程序和實驗進行比較。針對有代表性的金屬和PVC管線，在不同深度、不同管徑組合以及介質損耗大小等情況下獲得了實測的雷達剖面。相同模型下的數(shù)值模擬結果與實驗符合得較好。討論了一些以前較少關注的問題，如電磁波在被探物體和水箱前壁間多次反射形成的多條雙曲線問題。對于在較強介質損耗情況下的實測與數(shù)值模擬不符現(xiàn)象，也給出了可能的解決辦法。

[1] 蔚建斌, 陳自力, 江 濤. 探地雷達在地下介質中傳播路徑的計算[J]. 電波科學學報, 2009, 24(5): 925-933.

WEI Jianbin, CHEN Zili, JIANG Tao. Calculation of propagation path of GPR in media subsurface[J]. Chinese Journal of Radio Science, 2009, 24(5): 925-933. (in Chinese).

[2] BOURGEOIS J M, SIMITH G S. A fully three-dimensional simulation of a ground-penetrating radar: FDTD theory compared with experiment[J]. IEEE Trans Geosci Remote Sensing, 1996, 34(1): 36-44.

[3] CHEN Howwei, HUANG Taimin. Finite-difference time-domain simulationof GPR data[J]. Journal of Applied Geophysics, 1998, 40(1/3): 139-163.

[4] TEIXEIRA F L, WENG C C, STRAKA M, et al. Finite-difference time-domain simulation of ground penetrating radar on dispersive, inhomogeneous, and conductive soils[J]. IEEE Trans Geosci Remote Sensing, 1998, 36(6): 1928-1937.

[5] 周 峰, 潘和平, 文國軍, 等. 基于有限差分的井中激發(fā)極化三維正演[J]. 電波科學學報, 2010, 25(4): 785-792.

ZHOU Feng, PAN Heping, WEN Guojun, et al. Three-dimensional forward modeling of induced of polarization borehole using finite difference method [J]. Chinese Journal of Radio Science, 2010, 25(4): 785-792. (in Chinese)

[6] ROBERTS R L, DANIELS J. Finite difference time domain (FDTD) forward modeling of GPR data[C]// 5th Int Conf Ground Pentrating Radar. Waterloo, 1994: 185-204.

[7] 鄭奎松, 葛德彪, 葛 寧. 三維電磁散射的網絡并行FDTD計算和加速比分析[J]. 電波科學學報, 2004, 19(6): 767-771.

ZHENG Kuisong, GE Debiao, GE Ning. Parallel FDTD computing for 3D EM scattering problem and its speedup factor analysis[J]. Chinese Journal of Radio Science, 2004, 19(6): 767-771. (in Chinese).

[8] 馮德山, 戴前偉, 何繼善. 探地雷達的正演模擬及有限差分波動方程偏移處理[J]. 中南大學學報:自然科學版, 2006, 37(2): 361-365.

FENG Deshan, DAI Qianwei, HE Jishan. Forward smulation of ground penetrating radar and its finite difference method wave equation migration processing[J]. Journal of Central South University: Science and Technology, 2006, 37(2): 361-365. (in Chinese).

[9] YEE K S. Numerical solution of initial boundary value problems involving Maxwell’s equation in isotropic media [J]. IEEE Trans Antennas Propagat, 1966, 14(3): 302-307.

[10] 馮德山. 基于小波多分辨探地雷達正演及偏移處理研究[D]. 長沙: 中南大學信息物理工程學院, 2006.

FENG Deshan. Research Forward Smulation and Migration Processing of GPR Based on the Wavelet Multi-Resolution[D]. Changsha: Central South University, 2006. (in Chinese).

[11] BERENGER J P. A perfectly matched layer for the absorption of electromagnetic waves [J]. Journal of Computational Physics, 1994, 114 (2): 185-200.