發(fā)射脈沖波形對探地雷達(dá)響應(yīng)的影響

余小東，陸從德，杜興忠，武 瑩

（1.成都理工大學(xué) 地球探測與信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，成都 610059；2.貴陽勘測設(shè)計(jì)研究院 工程物探測試分院，貴陽 550081）

0 前言

探地雷達(dá)以其高分辨率、高效率、無損、結(jié)果直觀等優(yōu)點(diǎn)，已成為淺層地球物理勘探中一種常規(guī)的探測技術(shù)，其應(yīng)用領(lǐng)域不斷擴(kuò)展，并受到人們越來越多的重視［1－5］。因?yàn)榉直媛屎吞綔y深度是探地雷達(dá)研究的重要內(nèi)容，而影響分辨率和探測深度的因素有很多，為此許多學(xué)者在這方面展開了深入地研究。目前探地雷達(dá)分辨率和探測深度的影響因素可以分為電性參數(shù)、系統(tǒng)參數(shù)和處理方法三大類別。Smith－Rose和Scott等人［6－7］通過實(shí)驗(yàn)分析了電導(dǎo)率和介電常數(shù)等電性參數(shù)與頻率之間的關(guān)系，指出電性參數(shù)與使用的頻率有著緊密的聯(lián)系，而頻率又對探地雷達(dá)分辨率和探測深度有著很大的影響。Jol、王水強(qiáng)和戴前偉等［8－10］對探地雷達(dá)的發(fā)射天線中心頻率、激發(fā)能量、時(shí)窗、采樣率、測點(diǎn)間距和收發(fā)天線間距等系統(tǒng)參數(shù)，與探測深度、分辨率和反射連續(xù)性的關(guān)系進(jìn)行實(shí)際案例分析。這些研究不僅對工程應(yīng)用有重要的指導(dǎo)作用，而且還促進(jìn)了探地雷達(dá)理論的發(fā)展。此外，探地雷達(dá)數(shù)據(jù)處理方法對探測深度和分辨率也有影響，例如曾昭發(fā)等［1］認(rèn)為，有效識別并消除各種干擾噪聲是地質(zhì)雷達(dá)數(shù)據(jù)處理與資料解釋的關(guān)鍵環(huán)節(jié)，也是改善數(shù)據(jù)質(zhì)量、提高探地雷達(dá)探測分辨率的重要手段。

實(shí)際上在電磁勘探中，激發(fā)波形對電磁響應(yīng)的影響是比較顯著的，例如在航空瞬變電磁研究中，Liu［11］分析了激發(fā)波形對航空瞬變電磁響應(yīng)的影響。相似地，在探地雷達(dá)探測中，其激發(fā)波形也對電磁響應(yīng)有重要的影響，因此對激發(fā)波形的影響進(jìn)行分析和研究成為了本文的主要研究內(nèi)容。作者首先對漢寧、高斯、柯西和泊松四種窗函數(shù)及相應(yīng)的激發(fā)波形進(jìn)行分析；然后通過二維時(shí)域有限差分正演模擬獲得上述不同激發(fā)波形在不同地電模型中的響應(yīng)；最后分析這些響應(yīng)結(jié)果并總結(jié)出不同發(fā)射脈沖波形對探地雷達(dá)響應(yīng)的影響。

1 發(fā)射脈沖波形分析

目前，大多商用探地雷達(dá)系統(tǒng)采用高斯脈沖形式的調(diào)幅脈沖源，其子波的表達(dá)式為［1］：

其中 ω0為發(fā)射天線中心頻率；α為電磁脈沖的衰減系數(shù)。為了提高探地雷達(dá)的分辨率，其天線的頻率響應(yīng)要保證頻帶寬度等于其中心頻率，即當(dāng)α＝0.93ω0時(shí)，頻帶寬度可滿足要求［12］。圖1為中心頻率分別為100MHz、300MHz、500MHz的脈沖子波歸一化后的波形圖。

圖1 調(diào)幅脈沖源子波波形Fig.1 The wavelet waveform of amplitude pulse

從圖1可以看出，脈沖子波的寬度與天線中心頻率之間有著直接地聯(lián)系，隨著中心頻率的增大，脈沖子波的寬度逐漸變窄，從而使得探測分辨率隨之提高，但是這同時(shí)也會減小探測深度。這說明發(fā)射脈沖波形也是影響探地雷達(dá)分辨率和探測深度的因素之一。

Irving［13］在進(jìn)行探地雷達(dá)正演模擬時(shí)，采用Blackman－Harris窗函數(shù)的一階導(dǎo)數(shù)，來產(chǎn)生一個雙極性的脈沖子波作為發(fā)射波。Harris［14］給出了漢寧窗、高斯窗等多種窗函數(shù)的離散表達(dá)式，并分析了各種窗函數(shù)的特點(diǎn)。在窗函數(shù)的窗寬相同的條件下，不同類型窗函數(shù)的形狀不同，主要表現(xiàn)在能量的集中度不同以及上升沿或下降沿時(shí)間不同，因此各種窗函數(shù)的一階導(dǎo)數(shù)所形成的雙極性波形也存在著差異。也就是說當(dāng)發(fā)射天線中心頻率相同時(shí)，不同的窗函數(shù)可以形成不同的發(fā)射脈沖波形。實(shí)際上，由于激發(fā)波形的能量集中度不一樣，所以不同激發(fā)波形對探地雷達(dá)響應(yīng)的影響存在差異，這正是本文分析的理論基礎(chǔ)。下面主要討論漢寧、高斯、柯西、泊松這四種窗函數(shù)形成的發(fā)射脈沖波形對探地雷達(dá)響應(yīng)的影響。

（2）高斯窗函數(shù)的離散表達(dá)式為式（3）。

（3）柯西窗函數(shù)的離散表達(dá)式如式（4）表示。

（4）泊松窗函數(shù)的離散表達(dá)式如式（5）表示。

其中 N為窗寬離散采樣點(diǎn)數(shù)；α為窗的控制參數(shù)。

這四種窗函數(shù)的波形和其一階導(dǎo)數(shù)的波形分別如圖2、圖3所示，其中高斯、柯西和泊松窗函數(shù)的控制參數(shù)α分別為3.5、4.0和3.0。由圖2和圖3可知，從漢寧窗依次到泊松窗的能量集中度逐漸變小，通過一階求導(dǎo)后所得到的雙極性脈沖波的寬度也逐漸變小。因此不同窗函數(shù)在窗寬相同的條件下，可以產(chǎn)生不同寬度的雙極性脈沖波，使得在探地雷達(dá)中，在發(fā)射天線中心頻率保持不變的情況下，可以產(chǎn)生能量集中度不同的發(fā)射脈沖波形。

圖2 窗函數(shù)波形Fig.2 The waveform of window function

圖3 窗函數(shù)一階導(dǎo)數(shù)波形Fig.3 The waveform of the first derivative of window function

在探地雷達(dá)二維時(shí)域有限差分正演模擬中，發(fā)射脈沖波是關(guān)于時(shí)間的函數(shù)，因此需要將上面幾種窗函數(shù)的表達(dá)式表示成關(guān)于時(shí)間的離散函數(shù)。設(shè)T為發(fā)射脈沖持續(xù)的時(shí)間，fc為發(fā)射天線的中心頻率，Δt為離散時(shí)間間隔，則有

將式（6）代入式（2）可得漢寧窗函數(shù)的時(shí)間離散表達(dá)式為［15］式（7）：

設(shè)置w ij和w jk的初始連接權(quán)值，初始連接權(quán)值是在(-1，1)區(qū)間隨機(jī)選取的非零值，同時(shí)給定計(jì)算精度值ε(ε＞0)。

同理，將式（6）分別代入式（3）、式（4）、式（5）可得下面三種離散表達(dá)式。

（1）高斯窗函數(shù)的時(shí)間離散表達(dá)式。

（2）柯西窗函數(shù)的時(shí)間離散表達(dá)式。

（3）泊松窗函數(shù)的時(shí)間離散表達(dá)式。

2 二維時(shí)域有限差分正演模擬原理

時(shí)域有限差分（Finite Difference Time Domain Method，簡稱 FDTD）方法，是1966年 Yee［16］首次提出的一種以物理光學(xué)原理為基礎(chǔ)的高頻電磁場模擬方法，該方法具有存儲量小、運(yùn)算快、程序通用以及適于并行計(jì)算等優(yōu)點(diǎn)，是一種比較理想的探地雷達(dá)正演模擬工具［17］。

Yee氏網(wǎng)格模型（即Yee元胞）描述了時(shí)域有限差分法中電磁場空間分布情況，即每個電場E的分量周圍有四個磁場H分量環(huán)繞，并且每個磁場H的分量周圍有四個電場E分量環(huán)繞。時(shí)域有限差分法是對電磁場E和H分量在時(shí)間和空間上采用半步長交替網(wǎng)格的離散形式，并應(yīng)用這種離散方式將含時(shí)間變量的麥克斯韋旋度方程轉(zhuǎn)化為一組差分方程，然后在時(shí)間軸上逐步推進(jìn)地求解空間的電磁場［18］。

在二維正演模擬中，電磁場可分為TM波和TE波兩種。這里主要討論TM波的二維時(shí)域有限差分正演，這時(shí)有Ex＝Ez＝Hy＝0，從而獲得Ey、Hx、Hz空間網(wǎng)格分布如圖4所示。

圖4 二維時(shí)域有限差分網(wǎng)格Fig.4 Spatial arrangement of Ey，Hz，Hxfield components for TM－model modeling

利用四階中心差分格式將連續(xù)的電磁場變量進(jìn)行離散化，從而得到TM模式下探地雷達(dá)二維時(shí)域有限差分正演模擬方程［13］：

3 正演結(jié)果分析

目前所用的雙天線探地雷達(dá)測量方式主要有剖面法和寬角法兩種。作者采用剖面法的測量方式進(jìn)行正演模擬，發(fā)射天線和接收天線置于地表，按固定間距沿測線同步移動發(fā)射和接收天線。為了更真實(shí)地模擬雷達(dá)波在地下的傳播特征，創(chuàng)建模型時(shí)，在模型地表上添加一層空氣層，空氣層厚度設(shè)為模型厚度的1／10，空氣層的相對介電常數(shù)設(shè)為“1.0”，電導(dǎo)率設(shè)為“0”；邊界條件采用卷積完全匹配層（CPML）吸收邊界條件［19］，邊界設(shè)為10層網(wǎng)格。

3.1 鋼筋模型（橫向分辨率分析）

鋼筋模型的模型區(qū)域?yàn)?.0m×0.6m，背景區(qū)域參數(shù)εr＝7，σ＝0.005S／m；模型區(qū)域內(nèi)有四根圓形鋼筋異常體，其電性參數(shù)相同，即εr＝1.0×106，σ＝500S／m。左邊第一根鋼筋模型中心位置為（0.25m，0.25m），其余三根鋼筋中心位置和左邊第一根深度相同，測線方向相距0.15m，模型圖形詳見圖5。此外，網(wǎng)格步長Δx＝Δy＝0.01m，發(fā)射天線中心頻率為800MHz，發(fā)射和接收天線間距為0.05m，發(fā)射和接收天線每次同步移動的距離為0.02m，采樣時(shí)間點(diǎn)數(shù)取768點(diǎn)。采樣時(shí)間步長由程序根據(jù)最大時(shí)間步長確定且最大時(shí)間步長計(jì)算式為式（12）［13］。

圖5 鋼筋模型Fig.5 Rebar model

對于漢寧、高斯、柯西、泊松四種窗函數(shù)所產(chǎn)生的發(fā)射脈沖波形（分別簡稱：漢寧波形、高斯波形、柯西波形和泊松波形），其正演結(jié)果的wiggle圖和掃描圖分別如圖6～圖9所示。

圖6 漢寧波形的正演結(jié)果Fig.6 The forward results of Hanning waveform（a）wiggle圖；（b）掃描圖

圖7 高斯波形的正演結(jié)果Fig.7 The forward results of Gaussian waveform（a）wiggle圖；（b）掃描圖

圖8 柯西波形的正演結(jié)果Fig.8 The forward results of Cauchy waveform（a）wiggle圖；（b）掃描圖

圖9 泊松波形的正演結(jié)果Fig.9 The forward results of Poisson waveformm（a）wiggle圖；（b）掃描圖

從圖6－圖9可以看出，隨著發(fā)射脈沖波形的脈沖能量集中度變化，其正演結(jié)果圖形各不相同。①漢寧波形的正演結(jié)果效果最差，不能很好地將鋼筋分辨出來；②高斯波形的正演結(jié)果比漢寧波形的正演結(jié)果要好，可以將鋼筋模型分辨出來，并且可以判斷出鋼筋模型在水平方向的位置，但和柯西波形的正演結(jié)果相比，其分辨率相對較低；③泊松波形和柯西波形的正演結(jié)果都能很好地將鋼筋模型分辨出來，同時(shí)還可以計(jì)算出鋼筋模型的位置，且與設(shè)置的模型基本一致。但是在較深部的介質(zhì)中，泊松波形的正演結(jié)果存在的干擾信息比較嚴(yán)重（見圖9中大于7ns的響應(yīng)），使得其結(jié)果較柯西波形的正演結(jié)果差，限制了其探測深度。

在實(shí)驗(yàn)中我們還發(fā)現(xiàn)，如果繼續(xù)減小鋼筋模型之間的間距，漢寧波形和高斯波形的正演結(jié)果就很難分辨，而柯西波形和泊松波形可以獲得更高分辨率的正演結(jié)果。因此在保證發(fā)射天線中心頻率不變的情況下，使用寬度較窄的發(fā)射脈沖波形，可以改善探地雷達(dá)的橫向分辨率，但是激發(fā)波形的寬度并不是越小越好，因?yàn)楫?dāng)寬度越小，其干擾就會越大。

3.2 層狀模型（縱向分辨率分析）

層狀模型的模型區(qū)域?yàn)?.0m×1.0m，背景區(qū)域參數(shù)εr＝7，σ＝0.005S／m，模型區(qū)域內(nèi)有一薄層，薄層范圍為0.5m～0.58m，其電性參數(shù)設(shè)置為εr＝20，σ＝0.005S／m，模型圖形詳見圖10。網(wǎng)格步長Δx＝Δy＝0.01m，采樣時(shí)間點(diǎn)數(shù)取1 024點(diǎn)，發(fā)射天線中心頻率為800MHz，發(fā)射和接收天線間距以及每次同步移動的距離同鋼筋模型相同。對于漢寧波形、高斯波形、柯西波形和泊松波形，其正演結(jié)果wiggle圖和掃描圖分別如圖11～圖14所示。

圖10 薄層模型Fig.10 Thin layer model

從圖11－圖14可以看出：①漢寧波形和高斯波形的正演結(jié)果中可以很清楚地將薄層分辨出來，通過計(jì)算可以得出薄層的深度和厚度，不過高斯波形的正演結(jié)果比漢寧波形要好；②柯西波形和泊松波形的正演結(jié)果中，干擾信息較多，其縱向分辨率比漢寧波形和高斯波形都低，并且很難從其正演結(jié)果中正確地分辨出地下薄層模型深度和厚度。

在實(shí)驗(yàn)中，當(dāng)繼續(xù)減小薄層模型的厚度時(shí)，漢寧波形要先于高斯波形不能分辨出地下的薄層模型，因?yàn)闈h寧波形的脈沖寬較比高斯波形要寬，其縱向分辨率較低。例如，從圖15可看出，當(dāng)薄層模型厚度減小到0.05m時(shí)，就只有高斯波形的正演結(jié)果還能清晰地分辨出薄層模型。

圖11 漢寧波形的正演結(jié)果Fig.11 The forward results of Hanning waveform（a）wiggle圖；（b）掃描圖

圖12 高斯波形的正演結(jié)果Fig.12 The forward results of Gaussian waveform（a）wiggle圖；（b）掃描圖

圖13 柯西波形的正演結(jié)果Fig.13 The forward results of Cauchy waveform（a）wiggle圖；（b）掃描圖

圖14 泊松波形的正演結(jié)果Fig.14 The forward results of Poisson waveform（a）wiggle圖；（b）掃描圖

圖15 薄層厚度為0.05m的正演結(jié)果Wiggle圖Fig.15 The wiggle figure of forward results（a）漢寧波形；（b）高斯波形

4 結(jié)論

通過上述實(shí)驗(yàn)結(jié)果及分析，我們得出以下五點(diǎn)結(jié)論：

（1）在保持發(fā)射天線中心頻率不變的前提下，采用不同的窗函數(shù)可以產(chǎn)生不同脈沖寬度的激發(fā)波形。

（2）隨著窗函數(shù)能量集中度越高，一階求導(dǎo)后所產(chǎn)生的雙極性脈沖波形的脈沖寬度越小，探地雷達(dá)正演模擬結(jié)果分辨率越高，但其探測深度越淺。

（3）發(fā)射脈沖波形的寬度越寬，相應(yīng)波形具有的能量就越大，其穿透深度也就越大。因此在保持發(fā)射天線中心頻率不變的情況下，采用脈沖較寬的發(fā)射波形可以提高探測深度，但是探測分辨率會有所下降。

（4）當(dāng)發(fā)射脈沖寬度逐漸減小時(shí)，到一定勘探深度后，就會產(chǎn)生干擾，這也就限制了探測深度。所以當(dāng)探測淺部的目標(biāo)體時(shí)，應(yīng)采用寬度較窄的發(fā)射波形；而當(dāng)探測深部的目標(biāo)體時(shí)，則應(yīng)采用較寬的發(fā)射波形。

（5）對本文中的四種窗函數(shù)產(chǎn)生的發(fā)射波形而言，柯西脈沖波具有較好的橫向分辨率，而高斯脈沖波在橫向和縱向上都具有較好的分辨率。

［1］曾昭發(fā)，劉四新，馮暄.探地雷達(dá)原理與應(yīng)用［M］.北京：電子工業(yè)出版社，2010.

［2］陳義群，肖柏勛.論探地雷達(dá)現(xiàn)在與發(fā)展［J］.工程地球物理學(xué)報(bào)，2005，2（2）：149－155.

［3］李嘉，郭成超，王復(fù)明，等.探地雷達(dá)應(yīng)用概述［J］.地球物理學(xué)進(jìn)展，2007，22（2）：629－637.

［4］李華，魯光銀，何現(xiàn)啟，等.探地雷達(dá)的發(fā)展歷程及其前景探討［J］.地球物理學(xué)進(jìn)展，2010，25（4）：1492－1502.

［5］SLOB E，SATO M，OLHOEFT G.Surface and borehole ground－penetrating－radar development［J］.Geophysics，2010，75（5）：A103－A120.

［6］SMITH－ROSE R L.The electrical properties of soil［J］.Proceedings of the royal society，1933，841（140）：359－377.

［7］SCOTT J H，CARROLL R D，CUNNINGHAM D R.Dielectric constant and electrical conductivity measurements of moist rock：a new laboratory method ［J］.Journal of geophysical research，1967，72：5101－5115.

［8］JOL H M.Ground penetrating radar antennae frequencies and transmitter powers compared for penetration depth，resolution and reflection continuity［J］.Geophysical prospecting，1995，43（5）：693－709.

［9］王水強(qiáng)，萬明浩，謝雄耀，等.不同地質(zhì)雷達(dá)測量參數(shù)對數(shù)據(jù)采集效果的影響［J］.物探與化探，1999，23（3）：211－217.

［10］戴前偉，呂紹林，肖彬.地質(zhì)雷達(dá)的應(yīng)用條件探討［J］.物探與化探，2000，24（2）：157－160.

［11］LIU G M.Effect of transmitter current waveform on airborne TEM response［J］.Exploration Geophysics，1998，29（1）：35－38.

［12］馮德山.地質(zhì)雷達(dá)二維時(shí)域有限差分正演［D］.長沙：中南大學(xué)，2003.

［13］IRVING J，KNIGHT R.Numerical modeling of ground－penetrating－radar in 2－D using MATLAB［J］.Computers＆ Geosciences，2006，32：1247－1258.

［14］HAMIS F J.On the Use of Windows for Harmonic A－nalysis with the Discrete Fourier Transform ［J］.Proceedings of The IEEE，1978，66（1）：51－81.

［15］CHEN Y H，CHEW W C，ORISTAGLIO M L.Application of perfectly matched layers to the transient modeling of subsurface EM problems［J］.Geophysics，1997，62（6）：1730－1736.

［16］YEE K.S.Numerical Solution of Initial Boundary Val－ue Problems Involving Maxwell's Equation in Isotropic Media［J］.IEEE Transactions on Antennas and Propagation，1966，AP－14（3）：302－307.

［17］詹應(yīng)林，昌彥君，曹中林.基于UPML吸收邊界的探地雷達(dá)數(shù)值模擬研究［J］.資源環(huán)境與工程，2008，22（2）：35－38.

［18］葛德彪，閆玉波.電磁波時(shí)域有限差分方法（第二版）［M］.西安：西安電子科技大學(xué)出版社，2005.

［19］RODEN J A，GEDNEY S D.Convolution PML（CPML）：an efficient FDTD implementation of the CFS－PML for arbitrary media ［J］.IEEE Microwave and Optical Technology Letters，2000，27（5）：334－339.