HHT方法在探地雷達回波信號特征提取上的應用

甘 露 ，周 龍 ，尤新革 ，2

（1.武漢工業(yè)學院 湖北 武漢 430023；2.華中科技大學 湖北 武漢 430074）

探地雷達的信號處理結(jié)果決定著其在應用領域中的有效性，其處理過程主要是通過分析探地雷達接收到的有效反射電磁波的特征來推斷地下介質(zhì)的空間分布狀態(tài)。其中關(guān)于回波信號特征提取是最重要、最關(guān)鍵也最困難的問題之一，它直接關(guān)系到對地下目標探測的準確性和可靠性。傳統(tǒng)探地雷達資料常用的分析方法是Fourier分析[1-2]，然而，在傳統(tǒng)的探地雷達回波信號的表達中，一般都認為回波信號是平穩(wěn)的。但是，在實際應用中，大地的電磁特性呈現(xiàn)隨機性，導致接收的回波信號具有不同的幅度和相位，因此探地雷達的回波信號具有非平穩(wěn)的特征。在回波數(shù)據(jù)中，收發(fā)天線直接耦合波、地面反射波與它們之間多次波疊加，引起的反射波稱統(tǒng)稱為直達波。直達波的能量很大，某些情況下甚至能夠掩蓋目標。

希爾伯特一黃變換（Hilbert-Huang Transform，HHT）是一種新的非平穩(wěn)信號的處理技術(shù)，該方法由經(jīng)驗模態(tài)分解（Empirical Mode Decomposition，EMD）與Hilbert譜分析（Hilbert Spectral Analysis，HSA）兩部分組成：任意的非平穩(wěn)信號首先經(jīng)過EMD方法處理后被分解為若干個本征模態(tài)函數(shù)（Intrinsic Mode Function，IMF），然后對每個IMF分量進行Hilbert譜分析得到相應分量的Hilbert譜，最后匯總所有IMF分量Hilbert譜就得到了原始非平穩(wěn)信號的Hilbert譜[3-5]。按照這種方法分解所得到的IMF分量也具備明確的物理意義。這些IMF刻畫了信號在每一個局部的振蕩結(jié)構(gòu)或頻率結(jié)構(gòu)。為此，本文以探地雷達實測數(shù)據(jù)位研究對象，引入HHT將各種頻率成分以IMF的形式從原始信號中分離出來，對探地雷達回波信號進行特征提取，并給出探地雷達信號EMD分解算法。

1 HHT方法

Hilbert-Huang變換按照如下的2個步驟來分析數(shù)據(jù)。1）用經(jīng)驗模式分解對數(shù)據(jù)作預處理，通過分解我們可以得到一組“本征模式函數(shù)”，也即我們可以用得到的“基”展開數(shù)據(jù)。2）將分解得到的“本征模式函數(shù)”作希爾伯特變換并構(gòu)造能量-時間-頻率分布，即希爾伯特譜，它將保存事件的時間局部性。換句話說，我們需要的是瞬時的頻率和能量，而不是傅立葉譜分析中的全局頻率和能量。經(jīng)驗模式分解（EMD）的本質(zhì)是根據(jù)數(shù)據(jù)的特征時間尺度來經(jīng)驗地識別固有振蕩模式，然后據(jù)此來分解數(shù)據(jù)。不同于信號的Fourier分解和小波分解.EMD沒有確定的基，它的“基”是根據(jù)信號而自適應產(chǎn)生的，這使得它不僅具有很高的分解效率，同時也具有良好的時頻局部性。EMD分解的具體步驟如下：

設 x（t）為原始數(shù)據(jù)。

步驟 1 初始化:令 r0（t）=x（t），i=1

步驟 2 提取第i個IMF：

1）初始化:h0（t）=xi-1（t），j=1

2）求 hj-1（t）的局部極大值和局部極小值；

3）以局部極大值作為節(jié)點作三次樣條插值形成hj-1（t）的上包絡 u pper（t），類似的形成下包絡 l ower（t）；

5）令 hj（t）=hj-1（t）-mj-1（t）；

6）若hj（t）的穿越零點的個數(shù)與極值點的個數(shù)最多相差為 1 ，且

小于某個預先給定的常數(shù)，則 I MFi（t）=hj（t），否則轉(zhuǎn)至2），并令 j =j+1；

步驟 3 令 ri（t）=ri-1（t）-IMFi（t）；

步驟 4 若ri（t）至少有2個極值點，轉(zhuǎn)至步驟2并令i=i+1；否則，ri（t）分解結(jié)束，為余量。這樣，原信號 X（t）可表示為

2 系統(tǒng)硬件設計

雷達系統(tǒng)硬件主要由控制顯示單元、雷達主機和雷達天線及發(fā)射接收系統(tǒng)3部分組成，系統(tǒng)硬件結(jié)構(gòu)圖如圖1所示[6]。系統(tǒng)采用“PC＋DSP＋FPGA”三層控制方案。底層的控制器由以Altera公司的EP1K30TC144-3為核心的高速控制邏輯電路構(gòu)成，該芯片具有1728個邏輯單元和171個I/O引腳，可以滿足雷達時序系統(tǒng)的存儲容量需要和接口需要。FPGA不僅可以實現(xiàn)對時變增益放大器以及模數(shù)轉(zhuǎn)換器的同步控制，還同時生成發(fā)射和取樣觸發(fā)脈沖源。該脈沖源經(jīng)放大整形后送往天線系統(tǒng)；中層控制由AD公司的ADSP21065L芯片完成，該芯片采用超級哈佛（Super Harvard）結(jié)構(gòu)，其峰值運算速度可以達到198MFLOPS，它還提供了多種外部接口，具有10個DMA通道，可以為雷達數(shù)據(jù)的采集傳送提供快速通道。DSP不僅可實現(xiàn)與PC機的通信，接收并執(zhí)行由PC機發(fā)送的控制命令，實施對底層控制器的控制，同時還可將采集數(shù)據(jù)讀出并迅速上傳至PC機進行處理；上層控制由PC機完成對探地雷達功能和參數(shù)的設置，并向DSP發(fā)送控制命令，同時還實時地接收DSP上傳的回波信號數(shù)據(jù)，進行基于HHT的時頻分析并顯示圖象及波形。

3 基于HHT的探地雷達信號特征提取

根據(jù)上述分析，使用HHT方法進行探地雷達信號特征提取[7]的過程如下：

圖1 雷達硬件結(jié)構(gòu)圖Fig.1 Structure diagram of the radar hardware system

1）移動雷達對地下埋藏物的狀態(tài)進行數(shù)據(jù)采集，得到探地雷達回波信號的時域波形。

2）1）中所得的時域波形為一維信號，信號中包含直達波及目標信號。設接收到信號表示為d（x，t），其中x表示測點，t表示時間，則任意測點x處探地雷達信號的EMD分解可表示為：

式中，k表示第級IMF分量。在對實際探地雷達信號進行EMD分解時，當殘留分量剖面的能量遠遠小于原始剖面能量時，EMD分解就可以停止，則確定了其分解級數(shù)。

3）通過EMD分解，得到具有窄帶特性的內(nèi)蘊模式函數(shù)剖面imfk（x，t），進一步求取其瞬時頻率作為探地雷達回波數(shù)據(jù)特征。

4 實驗結(jié)果與分析

文中選取探地雷達實測數(shù)據(jù)中的一組數(shù)據(jù)進行測試。實驗場地選在戶外的一段水泥路面上，已知該路段下方有一水管道。將探地雷達工作參數(shù)的設置為：發(fā)射脈沖重復頻率為300 kHz，掃描速率為 256 scans/Sec，采樣率為 512 samp/Scan。采集到的回波信號如圖2所示。

圖2 探地雷達回波信號Fig.2 GPR echo signal

由圖2中可以看出在100點左右與400點左右均存在著明顯的幅值越變，而且幅值變化較劇烈。而位于前端的即為直達波信號，后端的為埋藏物信號。顯然，當直達波信號較強烈時，容易引起誤判。對圖2所示信號進行EMD分解，即可得到其IMF分量，如圖3所示。

圖3 回波信號的IMF分量Fig.3 IMF component of echo signal

由圖3可見，對圖2中的回波信號用EMD方法分解可得到3個IMF（imf1，imf2，imf3）分量和一個殘余分量（res），分解較為徹底，imf1～imf3頻率由高到低。然而，無法分辨出直達波與目標信號。進一步求取圖3中IMF分量的瞬時頻率，如圖4所示。

圖4 回波信號IMF分量的瞬時頻率Fig.4 EMD result of echo signal

由圖4可見，imf1的瞬時頻率表征了直達波與埋藏物信號處存在明顯變化。而imf2的瞬時頻率僅表征了埋藏物所處位置的頻率變化，說明在雷達經(jīng)過埋藏物上方時的回波信號頻率為明顯的高頻分量，體現(xiàn)了沖擊頻率。imf3和imf4的瞬時頻率分量則不包含過多信息。實驗結(jié)果表明：IMF的瞬時頻率能有效的提取出信號的主要特征信息，較好消除直達波的影響。

5 結(jié) 論

HHT作為一種新的信號處理方法，在非平穩(wěn)非線性信號的分析上有著獨特的優(yōu)勢，本文將HHT方法應用于探地雷達信號特征提取，并通過對實測數(shù)據(jù)的分析，得到各階固有模態(tài)函數(shù)信號。分析結(jié)果表明，采用HHT方法得到的IMF分量的瞬時頻率可較好的消除直達波影響，提取出地下埋藏物的特征，可供信號的進一步資料解釋。

[1]Maida A，Pennock S，Shepherd P.Improving ground penetrating radar signal analysis through FFT superimposition[J].Antennas and Propagation Society International Symposium，2005 IEEE，2005（4）:118.

[2]SONG Jia-yu，LIU Qing-Huo，Torrione P，et al.Two-dimensional and three-dimensional NUFFT migration method for landmine detection using ground-penetrating Radar[J].Geoscience and Remote Sensing，IEEE Transactions on，2006（6）:1462-1466.

[3]Huang N E，Shen Z，Long S R.A new view of nonlinear water waves：the Hilbert spectrum[J].Annual Review of Fluid Mechanics，1999（31）:417-457.

[4]Montesinos M E，Munoz-Cobo J L，Perez C.Hilbert-Huang analysis of BWR neutron detector signals：application to DR calculation and to corrupted signal analysis[J].Annals of Nuclear Energy，2003（30）：715-727.

[5]Phillips S C，Gledhill R J，Essex J W，et al.Application of the Hilbert-Huang transform to the analysis of molecular dynamics simulations[J].Journal of Physical Chemistry Ser A，2003（24）:4869-4876.

[6]GAN Lu，GAN Liang-cai，TIAN Mao.Step-system of high resolution ground penetrating radar[J].Chinese Journal of Radio Science，2008（3）:555-559.

[7]孫百紅，宋少偉.基于小波分析的發(fā)動機轉(zhuǎn)動慣量測量信號特征提取[J].火箭推進，2010（06）:52-55.SUN Bai-hong，SONG Shao-wei.Feature extraction for measurement signal of engine rotary inertia moment based on wavelet analysis[J].Journal of Rocket Propulsion，2010（06）:52-55.