王 超， 沈斐敏

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

?

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

王 超1，2， 沈斐敏2

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

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

探地雷達信號； 直達波去除； EMD； 差異擴大； 能量均衡

1 概述

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

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

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

2.1 基本思想

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

2.2 經驗模態(tài)分解(EMD)原理

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

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

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

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

(1)

c. 因為h1(t)一般仍不是一個IMF分量序列，為此需要對它重復進行上述處理過程。重復進行上述處理過程k次，直到h1(t)符合IMF的2個條件： ①極值點數(shù)目和過零點數(shù)目相等或最多相差1個； ②上下包絡關于時間軸局部對稱。這樣就得到了第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ù)據，重復步驟(a)，(b)和(c)，得到第二階IMF分量c2(t)，重復n次，得到n階IMF分量。這樣就有：

(5)

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

(6)

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

2.3 差異擴大

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

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

y=arcsin(t)

(7)

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

2.4 能量均衡

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

信號能量公式：

(8)

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

(9)

第i行似信號的經過能量擬合后的信號為：

(10)

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

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

2.5 計算流程

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

圖2 計算流程Figure 2 Calculation process

3 探地雷達信號直達波去除舉例

3.1 模擬信號的直達波EMD處理

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

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

圖4 第40道雷達波波形和經EMD處理后的第40道雷達波波形圖Figure 4 No.40 waveform and No.40 Waveform after direct wave removing

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

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

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

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

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

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

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

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

圖7 經EMD分解被去除的直達波信號和未經擴大和能量均衡的去除直達波后的雷達剖面

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

4 結論

① 本文嘗試采用EMD的原理在水平方向對探地雷達信號中的直達波進行去除，效果表明使用該方法直達波信號可以得到了有效的去除，被強振幅的直達波信號淹沒的信號則顯露出來；

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

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

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

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

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

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

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

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

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

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

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

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

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

[11] 周寶峰，任葉飛，溫瑞智，等.EMD方法在強震記錄處理中的應用研究[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 — )，男，江蘇徐州人，博士研究生，講師，研究方向：工程災害防范技術。

U 456.3

A

1674 — 0610(2016)05 — 0208 — 04