王 超， 沈斐敏

(1.福建船政交通職業(yè)學(xué)院， 福建 福州 350007； 2.福州大學(xué) 土木工程學(xué)院， 福建 福州 350002)

?

基于EMD的探地雷達(dá)信號(hào)直達(dá)波去除方法研究

王 超1，2， 沈斐敏2

(1.福建船政交通職業(yè)學(xué)院， 福建 福州 350007； 2.福州大學(xué) 土木工程學(xué)院， 福建 福州 350002)

結(jié)合探地雷達(dá)直達(dá)波的特點(diǎn)，提出了基于EMD的探地雷達(dá)直達(dá)波去除方法的基本思想，在去除直達(dá)波方法中引入了差異擴(kuò)大和能量均衡，并給出了計(jì)算流程。通過(guò)對(duì)數(shù)值模擬信號(hào)和現(xiàn)場(chǎng)實(shí)測(cè)信號(hào)的處理表明，基于EMD的直達(dá)波去除方法可以有效去除直達(dá)波信號(hào)，使用差異擴(kuò)大和能量均衡后可以較好地將有效信號(hào)顯現(xiàn)。

探地雷達(dá)信號(hào)； 直達(dá)波去除； EMD； 差異擴(kuò)大； 能量均衡

1 概述

探地雷達(dá)主要采用高頻電磁波(主頻為數(shù)十兆赫至數(shù)百兆赫以至千兆赫)以寬頻帶短脈沖形式，由地面通過(guò)發(fā)射天線送入被測(cè)對(duì)象，經(jīng)被測(cè)對(duì)象或目標(biāo)體反射后返回地面，為地面接收天線所接收，通過(guò)對(duì)接收波場(chǎng)的成像分析，獲取地下目標(biāo)的探測(cè)圖像。當(dāng)前探地雷達(dá)在隧道超前地質(zhì)預(yù)報(bào)、隧道襯砌質(zhì)量和公路路面及結(jié)構(gòu)質(zhì)量檢測(cè)等方面得到較大的應(yīng)用[1-4]。但是對(duì)于近距離瞬態(tài)測(cè)量，接收信號(hào)不可避免的含有發(fā)射信號(hào)的直接耦合部分，另外收發(fā)天線與被測(cè)物體并非緊密貼合，因此接收到的信號(hào)中含有信號(hào)經(jīng)被測(cè)對(duì)象表面反射回來(lái)的強(qiáng)反射信號(hào)，這些信號(hào)的能量遠(yuǎn)遠(yuǎn)強(qiáng)于目標(biāo)發(fā)射體的回波信號(hào)(本文將這些信號(hào)統(tǒng)稱為直達(dá)波信號(hào))。當(dāng)進(jìn)行淺層目標(biāo)體探測(cè)時(shí)，目標(biāo)體的發(fā)射信號(hào)與直達(dá)波信號(hào)相重疊，此時(shí)直達(dá)波信號(hào)對(duì)目標(biāo)體的識(shí)別產(chǎn)生嚴(yán)重的干擾，因此去除探地雷達(dá)直達(dá)波信號(hào)對(duì)淺層地下目標(biāo)體的識(shí)別、定位及解釋有重要的意義[5-7]。

去除直達(dá)波的方法有時(shí)間門限法、抵消法[8]、相干濾波法、模型法、小波濾波法[3，9]等。其中最簡(jiǎn)單的方法為抵消法，但是此方法不適用于目標(biāo)反射信號(hào)較弱的情況[9]；模型法、時(shí)間門限法都有其局限性，限制了方法的使用[8]。由于經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)在信號(hào)分解中可將高低頻分解的功能，因此本文嘗試使用EMD的方法對(duì)直達(dá)波進(jìn)行去除。

2 基于EMD去除探地雷達(dá)信號(hào)直達(dá)波原理

2.1 基本思想

基于EMD法去除直達(dá)波的基本思想為：在探地雷達(dá)淺表層探測(cè)中，反射信號(hào)只有在點(diǎn)目標(biāo)(如空洞、垂直放置的鋼筋、地雷等點(diǎn)狀目標(biāo))上方附近處的若干個(gè)掃描點(diǎn)上較強(qiáng)，在離開(kāi)目標(biāo)較遠(yuǎn)的掃描點(diǎn)上強(qiáng)度快速下降，而直達(dá)波在整個(gè)掃描區(qū)域都較強(qiáng)，而且變化幅度不大。如果將水平方向看做是信號(hào)的話，則直達(dá)波信號(hào)成分在水平方向表現(xiàn)為平穩(wěn)的穩(wěn)態(tài)信號(hào)，而點(diǎn)狀目標(biāo)反射信號(hào)成分在水平方向上則表現(xiàn)為高強(qiáng)度的快速變化的非穩(wěn)態(tài)信號(hào)。因此，本文依據(jù)此思想，將各成分在同一水平方向上的表現(xiàn)看做信號(hào)(本文稱為“似信號(hào)”)進(jìn)行EMD處理，從而達(dá)到將直達(dá)波信號(hào)去除的目的。

2.2 經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)原理

經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)基于3個(gè)假設(shè)： ①待處理信號(hào)中至少有一個(gè)極大值和一個(gè)極小值； ②時(shí)域特性由極值間隔決定； ③如果數(shù)據(jù)序列完全缺乏極值但是僅包含拐點(diǎn)，那么它也可通過(guò)一階求導(dǎo)或多階求導(dǎo)來(lái)表示極值點(diǎn)，而最終結(jié)果可以由這些成分求積分來(lái)獲得[10，11]。EMD的大體思路是利用時(shí)間序列上下包絡(luò)的平均值確定“瞬時(shí)平衡位置”，進(jìn)而提取固有模態(tài)函數(shù)。具體方法是由一個(gè)“篩選”過(guò)程完成的：

a. 找出原信號(hào)s(t)的局部極大值和極小值，利用三次樣條函數(shù)擬合成上、下包絡(luò)線。

b. 計(jì)算上、下包絡(luò)線的均值，記為m1(t)，把原數(shù)據(jù)序列s(t)減去該均值即可得到一個(gè)去掉低頻的新數(shù)據(jù)序列h1：

h1(t)=s(t)-m1(t)

(1)

c. 因?yàn)閔1(t)一般仍不是一個(gè)IMF分量序列，為此需要對(duì)它重復(fù)進(jìn)行上述處理過(guò)程。重復(fù)進(jìn)行上述處理過(guò)程k次，直到h1(t)符合IMF的2個(gè)條件： ①極值點(diǎn)數(shù)目和過(guò)零點(diǎn)數(shù)目相等或最多相差1個(gè)； ②上下包絡(luò)關(guān)于時(shí)間軸局部對(duì)稱。這樣就得到了第1階IMF分量c1(t)，它代表信號(hào)s(t)中最高頻率的分量：

h1k(t)=h1(k-1)(t)-m1k(t)

(2)

c1(t)=h1k(t)

(3)

d. 將c1(t)從s(t)中分離出來(lái)，即得到一個(gè)去掉高頻分量的差值信號(hào)c1(t)，即有：

r1(t)=s(t)-c1(t)

(4)

將r1(t)作為原始數(shù)據(jù)，重復(fù)步驟(a)，(b)和(c)，得到第二階IMF分量c2(t)，重復(fù)n次，得到n階IMF分量。這樣就有：

(5)

當(dāng)cn(t)或rn(t)滿足給定的終止條件(通常使rn(t)成為一個(gè)單調(diào)函數(shù))時(shí)，循環(huán)結(jié)束，由上面2個(gè)式子可以得到：

(6)

式中：rn(t)為殘余函數(shù)，代表信號(hào)的平均趨勢(shì)。而各個(gè)IMF分量c1(t)，c2(t)…cn(t)分別包含了信號(hào)不同時(shí)間特征尺度大小的成分，其尺度依次由小到大。因此，各分量也就相應(yīng)地包含了從高到底不同頻率段的成分，每一個(gè)頻率段所包含的頻率成分都是不同的，且隨信號(hào)本身的變化而變化。

2.3 差異擴(kuò)大

在探地雷達(dá)信號(hào)中，將受到直達(dá)波影響的區(qū)域內(nèi)選定，按照水平方向?qū)υ搮^(qū)域進(jìn)行EMD處理。由于EMD處理后，往往低階的固有模態(tài)函數(shù)(Intrinsic Mode Function，IMF)高頻成分較多，更多的表示為探測(cè)目標(biāo)體的特征；而高階的IMF信號(hào)低頻成分較多，更多的表現(xiàn)為受到直達(dá)波影響的特征，該部分往往變化比較平緩。因此可以將高階的IMF去除后，將剩余的各階IMF合成后表示原信號(hào)，由于反射目標(biāo)體信號(hào)成分在水平方向上的能量相對(duì)較弱，因此合成后的信號(hào)區(qū)分度不大，所以需要對(duì)該水平信號(hào)進(jìn)行差異擴(kuò)大。

本文采用反正弦函數(shù)表達(dá)式，如公式(7)所示，對(duì)數(shù)據(jù)歸一化后通過(guò)該函數(shù)處理，將使經(jīng)EMD處理后的似信號(hào)幅值表現(xiàn)為“強(qiáng)者更強(qiáng)、弱者更弱”的特征，數(shù)據(jù)差異的擴(kuò)大增加了數(shù)據(jù)的區(qū)分程度。如圖1所示為反正弦函數(shù)在[-1，1]區(qū)間曲線，可以看出，信號(hào)越接近于兩端極值，縱坐標(biāo)變換速度越快。

y=arcsin(t)

(7)

圖1 反正弦函數(shù)曲線Figure 1 Arc sine function curve

2.4 能量均衡

由于直達(dá)波信號(hào)主要表現(xiàn)在垂直方向上局部區(qū)域內(nèi)，因此只要對(duì)該部分似信號(hào)進(jìn)行處理即可。但是由于經(jīng)過(guò)處理后的似信號(hào)能量產(chǎn)生了變化，此時(shí)如果不進(jìn)行處理，該部分的似信號(hào)圖像與未選定區(qū)域沒(méi)有經(jīng)過(guò)處理的信號(hào)圖像不能夠很好的銜接，同樣無(wú)法顯示出直達(dá)波信號(hào)下淹沒(méi)的有用信號(hào)。因此需要對(duì)上述處理后的信號(hào)再進(jìn)行能量均衡。

信號(hào)能量公式：

(8)

若直達(dá)波信號(hào)在第a行至第b行區(qū)間內(nèi)，取a行似信號(hào)之前的前5行似信號(hào)的平均能量Ea賦給a行，b行似信號(hào)之后的后15行似信號(hào)的平均能量Eb賦值給b行，則線性均衡后直達(dá)波信號(hào)內(nèi)的第i行(i∈[a，b])似信號(hào)的平均能量為：

(9)

第i行似信號(hào)的經(jīng)過(guò)能量擬合后的信號(hào)為：

(10)

式中：Ei0為原信號(hào)x的能量。

通過(guò)上述各式，可以有效的將第a至第b行范圍內(nèi)的有效似信號(hào)與原處理信號(hào)很好的融合在一起，有效反映信號(hào)特征。

2.5 計(jì)算流程

在信號(hào)處理過(guò)程中，主要采用如圖2所示的計(jì)算流程。

圖2 計(jì)算流程Figure 2 Calculation process

3 探地雷達(dá)信號(hào)直達(dá)波去除舉例

3.1 模擬信號(hào)的直達(dá)波EMD處理

設(shè)計(jì)如圖3(a)所示探地雷達(dá)應(yīng)用幾何模型，模型長(zhǎng)1.0 m，高0.65 m，在坐標(biāo)(0.5，0.58)處有一圓形孔洞，半徑為2.25 cm。利用GPRMAX軟件進(jìn)行探地雷達(dá)探測(cè)正演數(shù)值模擬，天線發(fā)射頻率為800 M，從而獲得雷達(dá)波反射圖如圖3(b)所示，可以看出：在圖3(b)中直達(dá)波信號(hào)已經(jīng)將空洞處的反射信號(hào)完全淹沒(méi)。如圖4(a)所示為通過(guò)孔洞中心處的第40道雷達(dá)波波形圖，從中可以看出：在第40道雷達(dá)波信號(hào)中已經(jīng)完全淹沒(méi)了孔洞的反射信號(hào)。

圖3 幾何模型和雷達(dá)反射剖面Figure 3 Geometric modal and radar reflection profile

圖4 第40道雷達(dá)波波形和經(jīng)EMD處理后的第40道雷達(dá)波波形圖Figure 4 No.40 waveform and No.40 Waveform after direct wave removing

通過(guò)對(duì)第50行至170行的似信號(hào)進(jìn)行EMD分解，去除最后一階IMF和殘余函數(shù)后進(jìn)行合成，并進(jìn)行差異擴(kuò)大和能量均衡化處理后得到如圖5(a)反射圖像，從中可以看出：直達(dá)波信號(hào)幾乎完全被去除掉，孔洞的反射信號(hào)也有效的顯現(xiàn)出來(lái)。如圖5(b)為第50行至170行信號(hào)經(jīng)過(guò)EMD分解被去除的直達(dá)波信號(hào)，可以看出被處理區(qū)域的直達(dá)波信號(hào)已經(jīng)明顯發(fā)生畸變，但與圖3(b)比較，觀察出直達(dá)波信號(hào)振幅強(qiáng)弱區(qū)域仍然很明顯。如圖4(b)為經(jīng)過(guò)EMD處理后的第40道雷達(dá)波波形圖，比較圖4(a)可以看出圖4(b)的信號(hào)主要應(yīng)為孔洞的反射信號(hào)，直達(dá)波信號(hào)已經(jīng)被基本去除，但是去除后的信號(hào)產(chǎn)生了弱噪聲信號(hào)。

圖5 去除直達(dá)波后剖面和經(jīng)EMD分解被去除的直達(dá)波信號(hào)Figure 5 Profile after direct wave removing and direct wave signal after direct wave removing by EMD

3.2 工程信號(hào)的EMD分解與分析

某隧道進(jìn)行襯砌質(zhì)量檢測(cè)，如圖6(a)為未經(jīng)處理的探地雷達(dá)信號(hào)。從中可以看出：有效信號(hào)被直達(dá)波信號(hào)掩蓋的非常明顯。

根據(jù)上述處理程序，本次處理選取第20行到第80行之間的數(shù)據(jù)進(jìn)行處理，得到如圖6(b)的經(jīng)去除直達(dá)波處理后的圖像。從中可以看出在第20～80行數(shù)據(jù)內(nèi)的直達(dá)波信號(hào)幾乎被完全消除掉，隱藏在直達(dá)波內(nèi)部的信號(hào)顯露出來(lái)。如圖7(a)為被去除的直達(dá)波信號(hào)，如圖7(b)為未經(jīng)擴(kuò)大和能量均衡的去除直達(dá)波后的圖像。將圖6(b)與圖7(b)比較，可以看出經(jīng)過(guò)擴(kuò)大和能量均衡后，去除直達(dá)波之后淺部信號(hào)可以較好地顯示出來(lái)。

(a) 去除直達(dá)波前(b) 去除直達(dá)波后

Figure 6 Profile before direct wave removing and after direct wave removing

(a) 直達(dá)波信號(hào)(b) 雷達(dá)剖面

圖7 經(jīng)EMD分解被去除的直達(dá)波信號(hào)和未經(jīng)擴(kuò)大和能量均衡的去除直達(dá)波后的雷達(dá)剖面

圖7 Direct wave signal after direct wave removing by EMD & profile after direct wave removing without expanding and energy balance

4 結(jié)論

① 本文嘗試采用EMD的原理在水平方向?qū)μ降乩走_(dá)信號(hào)中的直達(dá)波進(jìn)行去除，效果表明使用該方法直達(dá)波信號(hào)可以得到了有效的去除，被強(qiáng)振幅的直達(dá)波信號(hào)淹沒(méi)的信號(hào)則顯露出來(lái)；

② 為增加信號(hào)的區(qū)分度，在經(jīng)EMD處理的信號(hào)歸一化后采用反正弦函數(shù)形式進(jìn)行處理；

③ 為使處理后圖像具有連續(xù)性的特點(diǎn)，文中采用能量均衡的方法將未經(jīng)處理信號(hào)與已處理的后能量產(chǎn)生改變的信號(hào)有效拼接起來(lái)；

④ 基于EMD的方法并非替代探地雷達(dá)直達(dá)波去除的傳統(tǒng)的方法，僅僅是對(duì)傳統(tǒng)方法的補(bǔ)充，當(dāng)使用傳統(tǒng)方法無(wú)法有效獲得有效信號(hào)時(shí)，可以嘗試此方法。

[1] 趙云威，朱自強(qiáng)，魯光銀，等.地質(zhì)雷達(dá)在玄武巖地區(qū)溶洞勘察中的應(yīng)用[J].公路工程，2012，37(3)：204-208.

[2] 陳培德.地質(zhì)雷達(dá)檢測(cè)技術(shù)在梧村隧道襯砌質(zhì)量檢測(cè)中的應(yīng)用[J].公路工程，2010，35(1)：134-137.

[3] 陳文超，師振盛，汪文秉，等.小波變換在去除探地雷達(dá)信號(hào)直達(dá)波的應(yīng)用[J].電波科學(xué)學(xué)報(bào)，2000，15(3)：352-357.

[4] 褚偉.GPR在張花高速公路結(jié)構(gòu)層厚度控制中的應(yīng)用[J].公路工程，2013，38(5)：34-37.

[5] 李昂，蔣延生，張安學(xué)，等.自適應(yīng)對(duì)消在去除探地雷達(dá)信號(hào)直達(dá)波的應(yīng)用[J].電波科學(xué)學(xué)報(bào)，2004，19(2)：223-227.

[6] 楊威，朱自強(qiáng)，魯光銀，等.Kirchhoff成像技術(shù)在路面無(wú)損檢測(cè)中的應(yīng)用[J].公路工程，2012，37(5)：79-82.

[7] 肖虎.淺層埋藏地雷GPR目標(biāo)回波信號(hào)處理研究[D].青島：中國(guó)海洋大學(xué)，2013.7：27.

[8] 高守傳，粟毅，黃春琳，等.用平均法實(shí)現(xiàn)瞬態(tài)信號(hào)接收中的直達(dá)波抑制[J].系統(tǒng)工程與電子技術(shù)，2004.1，26(1)：21-25.

[9] 楊秋芬，石顯新.探地雷達(dá)信號(hào)直達(dá)波的去除[J].煤炭學(xué)報(bào)，2008，33(7)：770-774.

[10] 孫曉云，周文佳，程久龍，等.基于二維HHT的隧道超前探測(cè)圖像識(shí)別與檢測(cè)[J].廣西師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2014，32(1)：26-31.

[11] 周寶峰，任葉飛，溫瑞智，等.EMD方法在強(qiáng)震記錄處理中的應(yīng)用研究[J].地震工程與工程振動(dòng)，2013，33(6)：9-15.

Study of the Direct Wave Removing in GPR Signal Based on EMD

WANG Chao1，2， SHEN Feimin2

(1.Fujian Chuanzheng Communications College， Fuzhou， Fujian 350007， China； 2.College of Civil Engineering， Fuzhou University， Fuzhou， Fujian 350002， China)

Combined with the GPR direct wave characteristics，it proposed the basic idea of the wave removing based on EMD，and introduced the expansion of difference and the energy balance，and provided the calculation process.By means of numerical simulation signal and the signal processing field measurement it showed that the direct wave removing method based on EMD can effectively remove direct wave signal，and effective signal display very well by using the expansion of difference and the energy balance.

the GPR signal； the direct wave removing； empirical mode decomposition； the expansion of difference； the energy balance

2014 — 11 — 13

福建省交通廳科技項(xiàng)目(201338)

王 超(1981 — )，男，江蘇徐州人，博士研究生，講師，研究方向：工程災(zāi)害防范技術(shù)。

U 456.3

A

1674 — 0610(2016)05 — 0208 — 04