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

王 超， 沈斐敏

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

?

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

王 超1，2， 沈斐敏2

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

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

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

1 概述

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

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

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

2.1 基本思想

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

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

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

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

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

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

(1)

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

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

(2)

c1(t)=h1k(t)

(3)

d. 將c1(t)從s(t)中分離出來，即得到一個去掉高頻分量的差值信號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)成為一個單調(diào)函數(shù))時，循環(huán)結(jié)束，由上面2個式子可以得到：

(6)

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

2.3 差異擴(kuò)大

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

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

y=arcsin(t)

(7)

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

2.4 能量均衡

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

信號能量公式：

(8)

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

(9)

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

(10)

式中：Ei0為原信號x的能量。

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

2.5 計算流程

在信號處理過程中，主要采用如圖2所示的計算流程。

圖2 計算流程Figure 2 Calculation process

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

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

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

圖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

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

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

3.2 工程信號的EMD分解與分析

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

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

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

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

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

圖7 經(jīng)EMD分解被去除的直達(dá)波信號和未經(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á)信號中的直達(dá)波進(jìn)行去除，效果表明使用該方法直達(dá)波信號可以得到了有效的去除，被強振幅的直達(dá)波信號淹沒的信號則顯露出來；

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

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

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

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

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

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

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

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

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

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

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

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

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

[11] 周寶峰，任葉飛，溫瑞智，等.EMD方法在強震記錄處理中的應(yīng)用研究[J].地震工程與工程振動，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

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

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

U 456.3

A

1674 — 0610(2016)05 — 0208 — 04