馮德山，戴前偉，余凱,

(1. 中南大學(xué) 地球科學(xué)與信息物理學(xué)院，湖南 長沙，410083；2. 水利部黃委會(huì)勘測(cè)規(guī)劃設(shè)計(jì)研究院，河南 鄭州，450003)

在探地雷達(dá)(Ground penetrating radar，GPR)的探測(cè)過程中，由于地下介質(zhì)結(jié)構(gòu)復(fù)雜、物性參數(shù)迥異和各種噪聲、雜波干擾的存在，導(dǎo)致雷達(dá)波在地層中傳播時(shí)變得十分復(fù)雜[1]。因此，如何有效地對(duì)GPR這種非平穩(wěn)、非線性信號(hào)[2]進(jìn)行處理，突破以處理穩(wěn)態(tài)數(shù)據(jù)為主的Fourier變換制約，成為GPR實(shí)際應(yīng)用中亟待解決的問題。柳剛等[3]應(yīng)用小波變換對(duì)低信噪比雷達(dá)信號(hào)進(jìn)行處理，論述了小波變換在處理低信噪比的瞬變信號(hào)的優(yōu)勢(shì)；張志禹等[4]對(duì)雷達(dá)記錄進(jìn)行小波分解并結(jié)合KL變換實(shí)現(xiàn)了目標(biāo)回波信號(hào)的串?dāng)_抑制；吳寶杰等[5]應(yīng)用S變換進(jìn)行雷達(dá)數(shù)據(jù)處理，去除了部分噪聲干擾，突出有效信號(hào)；吳健生等[6]應(yīng)用 Radon變換對(duì)GPR剖面中存在的“X”形同相軸等線段形進(jìn)行處理，消除了旁側(cè)干擾影響等。以上時(shí)頻方法在一定程度上提高了雷達(dá)的數(shù)據(jù)解譯精度，但大都屬于小波變換或以傅里葉變換為其最終理論依據(jù)，仍存在各自難以克服的缺點(diǎn)。以小波變換為例，雖然它具備了多分辨的性質(zhì)，但仍存在預(yù)選小波基函數(shù)、軟硬閾值參數(shù)設(shè)置等制約因素，限制了小波變換在GPR信號(hào)處理中的應(yīng)用和推廣。經(jīng)驗(yàn)?zāi)B(tài)分解 (Empirical mode decomposition，EMD)是一種分析非線性非平穩(wěn)信號(hào)的新方法[7?8]。它首先利用EMD方法將信號(hào)分解為若干固有模態(tài)函數(shù)(Intrinsic mode function，IMF)，然后，將 Hilbert變換作用在每一個(gè) IMF上，得到相應(yīng)的Hilbert瞬時(shí)譜，通過分析各個(gè)分量及其Hilbert譜，揭示原信號(hào)的多尺度振蕩變化特征。它既能使信號(hào)分解具有唯一性，又能在時(shí)域和頻域同時(shí)具有良好的局部化性質(zhì)。利用 EMD分解對(duì)信號(hào)分解完畢之后，可以根據(jù)工程問題的需要靈活地對(duì)信號(hào)實(shí)現(xiàn)重構(gòu)。它沒有先驗(yàn)的變換函數(shù)，但是，它具有小波變換的多分辨率性質(zhì)，又不需要像小波那樣考慮小波基的選取問題。經(jīng)過十幾年的發(fā)展，EMD分解在海洋氣象過程[9]、光譜數(shù)據(jù)預(yù)處理[10]、地球物理[11]等領(lǐng)域得到了廣泛研究和應(yīng)用。Flandrin等[12]通過對(duì)高斯白噪聲的分析，發(fā)現(xiàn)EMD方法表現(xiàn)為時(shí)域的二進(jìn)濾波；Wu等[13]把EMD應(yīng)用于基于噪聲統(tǒng)計(jì)規(guī)律的去噪；Damerval等[14]提出了二維EMD快速算法；湯井田等[15]應(yīng)用EMD分解對(duì)大地電磁信號(hào)進(jìn)行工頻干擾抑制和基線飄移矯正；段生全等[16]以地震資料為應(yīng)用實(shí)例，在物理意義、精確性和自適應(yīng)性方面全面比較了EMD分解(HHH變換)、Fourier 變換和小波分析法異同及優(yōu)劣；余志雄等[17?18]運(yùn)用Hilbert變換將探地雷達(dá)實(shí)信號(hào)轉(zhuǎn)換成復(fù)信號(hào)，分別提取瞬時(shí)振幅、瞬時(shí)相位、瞬時(shí)頻率波形圖，形成3個(gè)參數(shù)獨(dú)立剖面，提高了雷達(dá)解釋精度。在此，本文作者首先對(duì)低信噪比GPR信號(hào)進(jìn)行EMD分解，得到從高頻至低頻的GPR信號(hào)的IMF分量，通過消去屬于噪聲部分的本征模態(tài)函數(shù)分量，達(dá)到提高GPR數(shù)據(jù)信噪比的目的；然后，對(duì)去噪后的雷達(dá)剖面進(jìn)行重構(gòu)，再利用Hilbert變換求取GPR復(fù)信號(hào)，并提取瞬時(shí)振幅、瞬時(shí)相位、瞬時(shí)頻率3個(gè)參數(shù)，形成3個(gè)相互獨(dú)立的剖面，通過多參數(shù)波形剖面相互參照、綜合分析，避免了使用單一時(shí)距剖面分析所造成的解釋偏差。它與傳統(tǒng)復(fù)信號(hào)分析方法不同之處在于輸入的信號(hào)是經(jīng)過EMD分解之后的GPR信號(hào)，能根據(jù)信號(hào)處理的目的，針對(duì)性地去除低信噪比中各種噪聲，在賦于瞬時(shí)參數(shù)物理意義的同時(shí)，提高了GPR數(shù)據(jù)的解析精度。

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

經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)是Hilbert-Huang變換的關(guān)鍵組成部分，它有3個(gè)假定條件[19]：

(1) 待處理的信號(hào)中至少存在1個(gè)極大值和1個(gè)極小值。

(2) 由極值點(diǎn)間的間隔決定特征時(shí)間尺度。

(3) 若數(shù)據(jù)序列僅僅包含有拐點(diǎn)，可通過求 1階或者多階導(dǎo)數(shù)來確定極值點(diǎn)，并且最終結(jié)果也可由求積分來獲得。

在滿足了3個(gè)假定條件后，經(jīng)驗(yàn)?zāi)B(tài)分解方法認(rèn)為所有的信號(hào)都可由不同的本征模態(tài)函數(shù)IMF組成，其中，任意1個(gè)IMF都可以是線性的或者是非線性的。IMF都必須滿足以下2個(gè)條件：

(1) 對(duì)于 1列數(shù)據(jù)，極值點(diǎn)和過零點(diǎn)的個(gè)數(shù)都必須相等或者至多相差1點(diǎn)。

(2) 在任意點(diǎn)，由局部極大點(diǎn)和極小點(diǎn)構(gòu)成的 2條包絡(luò)線的平均值為0。每個(gè)IMF可以認(rèn)為是信號(hào)中固有的1個(gè)模態(tài)函數(shù)。

EMD的具體實(shí)現(xiàn)步驟可描述如下：設(shè)信號(hào)序列為f(t)，首先找出 f(t)中所有的極大值點(diǎn)和極小值點(diǎn)，通過三次樣條擬合，獲得f(t)的上包絡(luò)線u1(t)和下包絡(luò)線v1(t)，計(jì)算上、下包絡(luò)線在每點(diǎn)的平均值，從而獲得1條平均值曲線 m1(t)。用原始信號(hào) f(t)減去所得到的均值，得到1個(gè)新的數(shù)據(jù)序列h1(t)：

根據(jù)上面給出的IMF的判定條件對(duì)新的數(shù)據(jù)序列h1(t)進(jìn)行判定。若 h1(t)不滿足判定條件，則它不是 1個(gè)IMF分量序列，為此，重復(fù)進(jìn)行上述處理過程n次，使所得的hk(t)滿足IMF的判定條件，此時(shí)，hk(t)就是第 1階(IMF)c1(t)，這個(gè)分量代表信號(hào) f(t)中的最高頻率分量。然后，原始信號(hào)f(t)減去c1(t)，則可以得到去掉了高頻部分的差值信號(hào)序列r1(t)：

此時(shí)，將r1(t)作為待處理的原始數(shù)據(jù)序列，重復(fù)上述步驟，這樣，就可以得到第2個(gè)IMFc2(t)。重復(fù)操作n次，則可以得到n個(gè)本征模態(tài)函數(shù)分量，即：

當(dāng)rn(t)變成了1個(gè)常量或者成為1個(gè)單調(diào)函數(shù)時(shí)就停止運(yùn)算。原始信號(hào)f(t)可以表示為所有的IMF及剩余量之和：

式中：rn(t)為最后得到的殘余函數(shù)(單調(diào)函數(shù))，它代表了信號(hào)的平均趨勢(shì)。所有IMF分量都反映了信號(hào)的特征尺度，其尺度依次由小到大。所以，各個(gè)IMF分量相應(yīng)地包含了由高到低不同頻率段的成分，代表非線性形號(hào)的內(nèi)在模態(tài)特征，它隨著信號(hào)本身的變化而變化。

2 復(fù)信號(hào)分析理論

復(fù)信號(hào)分析又稱解析信號(hào)分析，就是把與記錄道相關(guān)的信息在時(shí)間域上直接分解為瞬時(shí)振幅、瞬時(shí)相位、瞬時(shí)頻率的一種處理和解釋技術(shù)[20]。在進(jìn)行復(fù)信號(hào)分析之前，首先要進(jìn)行Hilbert變換。Hilbert變換能有效、真實(shí)地獲取信號(hào)中所含的有效信息，它本質(zhì)上是一個(gè)全通濾波器[20]。設(shè)輸入的GPR信號(hào)為x(t)，它是經(jīng)過EMD分解并去噪處理的，通過濾波器H(ω)濾波后輸出信號(hào)為x?(t)。若H(ω)具有的幅頻特征是全通型的且相頻具有－90°相移，則

這時(shí)濾波器的輸出 x?(t)稱為x(t)的Hilbert變換。顯然，x?(t)與 x(t)正交。濾波器 H(ω)稱為 Hilbert濾波器。x(t)的 Hilbert變換可以記作 x?(t)或者 H·x(t)。而 x?(t)的定義為：

將 x(t)和它的 Hilbert變換x?(t)結(jié)合起來，組成 1個(gè)復(fù)信號(hào)，即[21]：

式中：u(t)為x(t)的復(fù)信號(hào)，又稱為解析信號(hào)。由于x(t)可以分解為三角函數(shù)形式，設(shè)x(t)=A(t)cos[ω0t+φ(t)]，(t)亦可表示為(t) = A(t) sin [ ω0t+ ? (t )](其中 ω0=2πf)，因此，x(t)的復(fù)信號(hào)又可表示為：

顯然，A(t)和 q(t)都隨時(shí)間而變化。A(t)稱為 u(t)的瞬時(shí)振幅；θ(t)=ω0t+φ(t)稱為u(t)的瞬時(shí)相位，相位的時(shí)間變化率為

S(t)即為所謂u(t)的瞬時(shí)頻率。當(dāng)φ(t)不變或變化不大時(shí)，φ′(t)可視為0或常數(shù)C，即S(t)=ω0+C只與頻率有關(guān)。

對(duì)于瞬時(shí)振幅、瞬時(shí)相位和瞬時(shí)頻率，可以用以下方法計(jì)算：首先由EMD分解后的GPR記錄x(t)經(jīng)Hilbert變換求得?(t)，然后計(jì)算瞬時(shí)振幅：

瞬時(shí)振幅是時(shí)間變量t的函數(shù)，與相位θ(t)無關(guān)。瞬時(shí)相位為

或

瞬時(shí)頻率 S(t)是瞬時(shí)相位函數(shù)對(duì)時(shí)間的變化率，即對(duì)θ(t)求導(dǎo)得：

復(fù)信號(hào)分析技術(shù)可以將GPR記錄中的瞬時(shí)振幅、瞬時(shí)相位和瞬時(shí)頻率分離出來，同一GPR的3個(gè)參數(shù)瞬時(shí)譜從不同的角度反映地下介質(zhì)的特性。瞬時(shí)振幅是對(duì)電磁波反射強(qiáng)度的量度，與該時(shí)刻GPR信號(hào)總能量的平方根成正比，利用這種特征可以對(duì)特殊巖層的變化作出判定。當(dāng)?shù)貙又写嬖谟斜容^明顯的分層、空隙或者是地下水的分界面時(shí)，在瞬時(shí)振幅剖面圖中會(huì)看到非常明顯的異常反應(yīng)。瞬時(shí)相位是GPR剖面上同相軸連續(xù)性的量度。雷達(dá)波在物性參數(shù)相同或相近的地下介質(zhì)中傳播，在相位圖中表現(xiàn)出來其相位是連續(xù)變化的；而雷達(dá)波在物性參數(shù)不同或相差較大的地下介質(zhì)中傳播，在異常體處相位會(huì)表現(xiàn)出相位錯(cuò)斷的強(qiáng)烈反應(yīng)。瞬時(shí)頻率所表示的是相位的時(shí)間變化率即相位的導(dǎo)數(shù)，反映出其相速度。在地層中介質(zhì)巖性的變化會(huì)引起頻率的變化，因而利用瞬時(shí)頻率也作為對(duì)地層分界面進(jìn)行判斷的依據(jù)。在這3個(gè)參數(shù)中，瞬時(shí)相位譜的分辨率最高，而瞬時(shí)頻率譜和瞬時(shí)振幅譜的變化也較直觀，所以，通常根據(jù)瞬時(shí)頻率譜和瞬時(shí)振幅譜來確定地下異?；蚍謱拥拇蟾盼恢茫缓?，利用瞬時(shí)相位譜精確確定異常位置和分層輪廓線。有時(shí)也可以直接利用瞬時(shí)相位譜來確定地下異常的位置。

3 基于 EMD分解的雷達(dá)數(shù)據(jù)處理實(shí)例

3.1 單道雷達(dá)信號(hào)的EMD分解

經(jīng)驗(yàn)?zāi)B(tài)分解是自適應(yīng)的，其分解快速而且有效，同時(shí)它又是基于信號(hào)的局部變化特性的，所以，特別適合于 GPR這種非線性、非平穩(wěn)信號(hào)的分析。圖1所示為1個(gè)單道雷達(dá)數(shù)據(jù)EMD分解結(jié)果。分析圖1中的各個(gè)IMF曲線的波形可知：由IMF1到IMF5，頻率成分隨著分量階數(shù)的增加移向低頻；第1階分量(IMF1)頻率最高，首先被提取出來，其中可能包含了高頻噪聲，余下分量的噪聲依次降低，說明 EMD分解具有去除高頻噪聲的功能；第 2和第 3階分量(IMF2，IMF3)代表了原雷達(dá)道的主要成分，IMF4及其之后的分量相對(duì)分量幅值較前3個(gè)分量小，不是雷達(dá)波的主要成分，余量反映了信號(hào)的平均趨勢(shì)。此外，原始圖、IMF1和IMF2中600到800號(hào)采樣點(diǎn)間反映并明顯的異常信息，在IMF3中反映非常清晰、突出，由此可見，經(jīng)驗(yàn)?zāi)B(tài)分解還具有分離高、低頻異常的作用。

3.2 GPR模擬數(shù)據(jù)的EMD分解去噪及復(fù)信號(hào)分析

圖1 單道GPR信號(hào)EMD分解圖Fig.1 Single trace GPR signal EMD map

通常使用GPR進(jìn)行勘探時(shí)，干擾普遍存在，如外部噪聲的干擾、復(fù)雜的地下圍巖的干擾等等。因而，在GPR信號(hào)處理過程中，如何消除這些干擾并從信號(hào)中得到有效的信息，這對(duì)于GPR資料的處理及解釋有著重大意義。以GPR正演模擬為基礎(chǔ)，并在GPR正演數(shù)據(jù)中加入隨機(jī)噪聲，以分析 EMD分解的去噪效果。設(shè)定雷達(dá)模型如圖2所示，其中：最上層是厚度為0.5 m的空氣介質(zhì)，空氣層以下為素填土層，素填土的介電常數(shù)為10.0，電導(dǎo)率為0.002 S/m；在素填土層中設(shè)置5個(gè)金屬球狀異常體，金屬球體的上頂埋度為6.0 m，其電導(dǎo)率為1.0×108S/m，半徑為1.0 m，金屬球體間距為1.5 m；模擬區(qū)域長×寬為20.0 m×10.0 m。脈沖波形為Ricker子波，頻率為60 MHz，脈沖位置距上邊界0.48 m處，計(jì)算的網(wǎng)格長與寬均為0.1 m；時(shí)窗為400 ns；采樣點(diǎn)從左邊0.2 m開始，每隔0.2 m采1個(gè)點(diǎn)，采樣過程中發(fā)送點(diǎn)與接收點(diǎn)是同步移動(dòng)的，收發(fā)距為0.8 m，每道波形有1 024個(gè)采樣點(diǎn)，共采集了94道雷達(dá)數(shù)據(jù)。采用有限差分法對(duì)該模型進(jìn)行正演，其正演合成剖面如圖3所示。圖4所示為雷達(dá)正演數(shù)據(jù)加入隨機(jī)噪聲后的剖面。由圖4可見：由于繞射現(xiàn)象及隨機(jī)噪聲的存在，導(dǎo)致信號(hào)異常混亂，GPR剖面的信噪比與分辨力大大降低。

為此，利用EMD對(duì)加噪后的GPR正演剖面進(jìn)行分解，得到圖5所示的分解后的本征模態(tài)函數(shù)圖IMF1～ IMF7。由于EMD分解實(shí)際上是一個(gè)高通濾波的過程，所以，在對(duì)信號(hào)的分解過程中高頻成分被慢慢濾出，圖5(a)～(b)中前2個(gè)分離出來的IMF完全是后期加入的高頻噪聲；圖5(c)～(d)中2個(gè)IMF是雷達(dá)剖面信號(hào)的主要成分，其中圖5(c)中的 IMF4為金屬球體的繞射異常；圖5(e)所示為分解得到的IMF7，由于GPR信號(hào)的有效部分基本被分離出來，故屬于金屬球體的異常在圖中得不到反映，后面的IMF分量圖不再列出。

圖2 雷達(dá)模型示意圖Fig.2 Sketch map of GPR model

圖3 雷達(dá)正演剖面圖Fig.3 GPR forward simulation map

圖4 加噪后的雷達(dá)正演剖面圖Fig.4 GPR forward simulation map of addition noise

通過分析EMD分解結(jié)果可知：加入的噪聲主要存在于IMF1～I(xiàn)MF2剖面中，為此去除IMF1～I(xiàn)MF2并對(duì)其它IMF分量中的有效信號(hào)進(jìn)行重構(gòu)，重構(gòu)后的新分量如圖5(i)所示。分析圖5(i)可知：通過EMD分解后去除主要由噪聲組成的 IMF1～I(xiàn)MF2后，大大降低了信號(hào)中的噪聲干擾，起到了類似濾波的效果。盡管對(duì)于GPR剖面的信噪比有較大提高，然而，剖面中異常的顯示并不很清晰。為此，在EMD分解的基礎(chǔ)上，對(duì)重構(gòu)后的圖5(i)信號(hào)進(jìn)行Hilbert變換，分別求得圖6所示的GPR瞬時(shí)振幅剖面、圖7所示的瞬時(shí)相位剖面和圖8所示的瞬時(shí)頻率剖面。由于這3個(gè)瞬時(shí)特征參數(shù)所對(duì)應(yīng)地下介質(zhì)的物性參數(shù)不同，其剖面圖表像各異，故根據(jù)不同的地質(zhì)狀況選擇合適的剖面圖對(duì)信號(hào)進(jìn)行分析比直接利用原始剖面圖要更加清晰、有效。對(duì)比原始剖面圖，瞬時(shí)相位圖與瞬時(shí)頻率圖更好地反映了異常信息。

3.3 對(duì)實(shí)測(cè)GPR數(shù)據(jù)進(jìn)行基于EMD分解的復(fù)信號(hào)分析

湖南長沙黑麋峰抽水蓄能電站進(jìn)廠交通隧洞未襯砌的J0+200.00 m～J0+984.00 m段圍巖地質(zhì)條件復(fù)雜，需詳細(xì)查明隧洞中存在的破碎帶、裂隙等不良地質(zhì)體，以確保進(jìn)出廠人員的作業(yè)安全。為此，在該隧洞段開展了探地雷達(dá)隱患探測(cè)。在探測(cè)過程中，采用 GSSI公司生產(chǎn)的 SIR?3000型雷達(dá)數(shù)據(jù)采集單元，選取頻率為900 MHz天線，點(diǎn)采樣方式，點(diǎn)距為0.20 m，每道掃描512個(gè)采樣點(diǎn)。

圖9(a)所示為該段中比較典型的雷達(dá)原始剖面圖，可見：由于隧洞中存在的鋼拱架、鐵絲網(wǎng)、電纜干擾等，導(dǎo)致該原始剖面異常復(fù)雜，直觀感覺剖面中存在異常體，但又不能準(zhǔn)確定位。顯然，僅根據(jù)該剖面圖無法對(duì)探測(cè)段進(jìn)行準(zhǔn)確的判斷和地質(zhì)解釋。為此，先對(duì)實(shí)測(cè)信號(hào)進(jìn)行EMD分解，得到了圖9(b)～(e)中分解后的本征模態(tài)函數(shù)圖IMF1～ IMF4(后面分量為無效信號(hào)，沒有列出)。由于EMD分解實(shí)際上是一個(gè)高通濾波的過程，信號(hào)中的高頻噪聲成分主要集中在IMF1剖面中，故去除IMF1可達(dá)到濾波消噪的目的。圖9(f)所示為去除IMF1后的重構(gòu)圖，該剖面中的異常信息比其在原始剖面圖中要更加清晰，但仍未達(dá)到準(zhǔn)確指導(dǎo)異常解釋的分辨率。對(duì) EMD分解后的信號(hào)進(jìn)行Hilbert變換，在此基礎(chǔ)上，進(jìn)行復(fù)信號(hào)分析以提取瞬時(shí)特征參數(shù)。圖10所示為重構(gòu)后的GPR瞬時(shí)振幅剖面，圖11所示為重構(gòu)后的GPR瞬時(shí)相位剖面，圖12所示為重構(gòu)后的GPR瞬時(shí)頻率剖面。對(duì)比3種瞬時(shí)剖面可知：經(jīng)過復(fù)信號(hào)分析處理的GPR剖面更加清晰直觀，分辨率更高，由于它們分別對(duì)應(yīng)地下介質(zhì)不同的物性參數(shù)，故其剖面圖表象各異，在這3個(gè)瞬時(shí)參數(shù)譜之中，瞬時(shí)相位譜效果最好，瞬時(shí)頻率譜次之。從3個(gè)瞬時(shí)剖面中，可以清楚地分辨位于道號(hào)50～125之間采樣點(diǎn)數(shù)150～200之間的破碎帶位置，后經(jīng)鉆探驗(yàn)證。因此，根據(jù)不同的地質(zhì)狀況選擇合適的剖面圖，對(duì)信號(hào)進(jìn)行分析會(huì)比直接利用原始剖面圖要更加清晰、有效，它能使實(shí)測(cè)信號(hào)中異常信息的識(shí)別變得更加簡單、準(zhǔn)確。若綜合3個(gè)參數(shù)的瞬時(shí)譜，利用多個(gè)參數(shù)聯(lián)合評(píng)估、解釋，則能達(dá)到相互驗(yàn)證，提高解釋精度的目的。

圖5 EMD分解得到的IMF1～I(xiàn)MF8及其重構(gòu)圖Fig.5 Maps of IMF1～I(xiàn)MF8 after EMD decomposition and reconstruction

圖6 正演資料瞬時(shí)振幅剖面Fig.6 Instantaneous amplitude profile of simulated data

圖7 正演資料雷達(dá)瞬時(shí)相位剖面Fig.7 Instantaneous phase profile of simulated data

圖8 正演資料雷達(dá)瞬時(shí)頻率剖面Fig.8 Instantaneous frequency profile of simulated data

圖9 雷達(dá)原始剖面與EMD分解得到的IMF1～I(xiàn)MF4及其重構(gòu)圖Fig.9 Original profiles of GPR and IMF1～I(xiàn)MF4 after EMD decomposition and reconstruction profile

圖10 實(shí)測(cè)資料瞬時(shí)振幅剖面Fig.10 Instantaneous amplitude profile of field data

圖11 實(shí)測(cè)資料雷達(dá)瞬時(shí)相位剖面Fig.11 Instantaneous phase profile of field data

圖12 實(shí)測(cè)資料雷達(dá)瞬時(shí)頻率剖面Fig.12 Instantaneous frequency profile of field data

4 結(jié)論

(1) EMD分解的本質(zhì)是從GPR剖面中從高至低逐步分離出不同頻率成分，并通過分屏的方式把各個(gè)頻率段所對(duì)應(yīng)的異常信息更有效地顯示出來，根據(jù)不同的數(shù)據(jù)處理目的，可靈活地對(duì)這些IMF雷達(dá)剖面圖進(jìn)行處理。EMD分解對(duì)于低信噪比雷達(dá)數(shù)據(jù)具有較強(qiáng)的適應(yīng)性，把 EMD分解應(yīng)用于含噪的雷達(dá)信號(hào)，去除含噪的IMF分量再重構(gòu)GPR剖面，可以較好地實(shí)現(xiàn)對(duì)低信噪比GPR數(shù)據(jù)的噪聲去除，達(dá)到改善GPR信號(hào)分析效果、提高信噪比和數(shù)據(jù)解釋精度的目的。

(2) 對(duì) GPR數(shù)據(jù)先進(jìn)行 EMD分解，然后結(jié)合Hilbert變換和復(fù)信號(hào)分析求取雷達(dá)信號(hào)的瞬時(shí)特征參數(shù)，并用所求得的瞬時(shí)振幅、瞬時(shí)相位、瞬時(shí)頻率的三瞬信息特征圖與原始信號(hào)的剖面圖相結(jié)合，應(yīng)用多參數(shù)分析方法來表達(dá)雷達(dá)剖面中蘊(yùn)含的信息，能夠使雷達(dá)剖面中的反射界面更加清晰可見，有助于突出雷達(dá)剖面中異常體特征，提高GPR數(shù)據(jù)處理及解釋的精度。

[1]曾昭發(fā), 劉四新, 王者江, 等. 探地雷達(dá)方法原理及應(yīng)用[M].北京: 科學(xué)出版社, 2006: 159?161.ZENG Zhao-fa, LIU Si-xin, WANG Zhe-jiang, et al. Theory and application of ground penetrating radar[M]. Beijing: Science Press, 2006: 159?161.

[2]翟波, 楊峰, 孫水明, 等. 基于二維濾波的探地雷達(dá)數(shù)據(jù)去噪研究[J]. 南京師范大學(xué)學(xué)報(bào): 工程技術(shù)版, 2007, 7(3): 79?83.ZHAI Bo, YANG Feng, SUN Shui-ming, et al. Research of GPR signal de-noising based on 2-D in ensional filtering[J]. Journal of Nanjing Normal University:Engineering and Technology Edition,2007, 7(3): 79?83.

[3]柳剛, 李術(shù)才, 薛翊國, 等. 基于小波變換的雷達(dá)低信噪比信號(hào)處理技術(shù)及應(yīng)用研究[J]. 工程勘察, 2009(9): 85?90.LIU Gang, LI Shu-cai, XUE Yi-guo, et al. GPR signal processing approach under low signal to noise ratio based on wavelet transforms and its application[J]. Geotechnical Investigation & Surveying, 2009(9): 85?90.

[4]張志禹, 劉亞麗. 基于小波域KL變換的探地雷達(dá)信號(hào)處理[J].煤田地質(zhì)與勘探, 2007, 35(2): 70?72.ZHANG Zhi-yu, LIU Ya-li. Signal processing of ground penetrating radar based on KL transforms in wavelet domain[J].Coal Geology & Explorations, 2007, 35(2): 70?72.

[5]吳寶杰, 楊樺, 張偉光. 探地雷達(dá)數(shù)據(jù)的 S變換時(shí)頻分析[J].上海地質(zhì), 2008, 107(3): 13?15, 19.WU Bao-jie, YANG Hua, ZHANG Wei-Guang. S-transform time-frequency analysis of GPR data[J]. Shanghai Geology,2008, 107(3): 13?15, 19.

[6]吳健生, 張昊. 拉東變換在探地雷達(dá)資料處理中的應(yīng)用[J].同濟(jì)大學(xué)學(xué)報(bào): 自然科學(xué)版, 2005, 33(9): 1270?1273.WU Jian-sheng, ZHANG Hao. Application of radon transform in ground penetrating radar data processing[J]. Journal of Tongji University :Natural Science, 2005, 33(9): 1270?1273.

[7]Huang N E, Shen Z, Long S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationarity time series analysis[J]. Proceedings of the Royal Society London: Series A, 1998, 454: 903?995.

[8]Huang N E, Wu M C, Long S R, et al. A confidence limit for the empirical mode decomposition and the Hilbert spectral analysis[J]. Proceedings of the Royal Society London: Series A,2003, 31: 417?457.

[9]熊學(xué)軍, 郭炳火, 胡筱敏, 等. EMD方法和Hilbert譜分析法的應(yīng)用與探討[J]. 黃渤海海洋, 2002, 20(2): 12?21.XIONG Xue-jun, GUO Bing-huo, HU Xiao-min, et al. The application and research of EMD and Hilbert[J]. Jounal of Oceanography of HuangHai & BoHai Seas, 2002, 20(2): 12?21.

[10]蔡劍華, 王先春. 基于經(jīng)驗(yàn)?zāi)B(tài)分解的近紅外光譜預(yù)處理方法[J]. 光學(xué)學(xué)報(bào), 2010, 30(1): 267?271.CAI Jian-hua, WANG Xian-chun. Near-Infrared spectrum pretreatment based on empirical mode decomposition[J]. Acta Optica Sinica, 2010, 30(1): 267?271.

[11]皮紅梅, 劉財(cái), 王典. 利用 Hilbert-Huang變換提取地震信號(hào)瞬時(shí)參數(shù)[J]. 石油地球物理勘探, 2007, 42(4): 418?424.PI Hong-mei, LIU Cai, WANG Dian. Using Hilbert-Huang transform to pick up instantaneous parameters of seismic signal[J]. Oil Geophysical Prospecting, 2007, 42(4): 418?424.

[12]Flandrin P, Rilling G, Goncalves P. Empirical mode decomposition as a filter bank[J]. Singnal Processing Letters,IEEE, 2004, 11(2): 112?114.

[13]Wu Z, Huang N E. A study of the characteristics of white noise using the empirical mode decomposition method[J]. Proceedings of the Royal Society London: Series A, 2004, 460: 1597?1611.

[14]Damerval C, Meignen S, Perrier V. A fast algorithm for bidimensional EMD[J]. Signal Processing Letters, IEEE, 2005,12(10): 701?704.

[15]湯井田, 化希瑞, 曹哲民, 等. Hilbert-Huang變換與大地電磁噪聲壓制[J]. 地球物理學(xué)報(bào), 2008, 51(2): 603?610.TANG Jing-tian, HUA Xi-rui, CAO Zhe-min, et al.Hilbert-Huang transformation and noise suppression of magnetotelluric sounding data[J]. Chinese Journal of Geophysics,2008, 51(2): 603?610.

[16]段生全, 賀振華, 黃德濟(jì). HHT方法及其在地震信號(hào)處理中的應(yīng)用[J]. 成都理工大學(xué)學(xué)報(bào): 自然科學(xué)版, 2005, 32(4):396?400.DUAN Sheng-quan, HE Zhen-hua , HUANG De-ji. Application of the Hilbert-Huang transform to the analysis of seismic signal[J]. Journal of Chengdu University of Technology: Science& Technology Edition, 2005, 32(4): 396?400.

[17]余志雄, 薛桂玉, 周創(chuàng)兵. 復(fù)信號(hào)分析技術(shù)及其在地質(zhì)雷達(dá)數(shù)字處理中的應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報(bào), 2005, 24(5):798?802.YU Zhi-Xiong, XUE Gui-yu, ZHOU Chuang-bin. The application of complex signals in GPR[J]. The Jurnal of Rock Mechanic and Engineering, 2005, 24(5): 798?802.

[18]謝雄耀, 萬明浩. 復(fù)信號(hào)分析技術(shù)在地質(zhì)雷達(dá)信號(hào)處理中的應(yīng)用[J]. 物探化探計(jì)算技術(shù), 2000, 22(2): 108?112.XIE Xiong-yao, WAN Ming-Hao. The application of complex signal in treatment of signal of GPR[J]. Geophysical Exploration and Chemical Exploration, 2000, 22(2): 108?112.

[19]湯井田, 蔡劍華, 任政勇, 等. Hilbert-Huang變換與大地電磁信號(hào)的時(shí)頻分析[J]. 中南大學(xué)學(xué)報(bào): 自然科學(xué)版, 2009, 40(5):1439?1405.TANG Jing-tian, CAI Jian-hua, REN Zheng-yong.Hilbert-Huang transform and time-frequency analysis of magnetotelluric signal[J]. Journal of Central South University of Technology: Science & Technology, 2009, 40(5): 1439?1405.

[20]肖兵, 鮑光淑, 趙秋梅, 等. 探地雷達(dá)復(fù)信號(hào)分析及改進(jìn)[J].中南工業(yè)大學(xué)學(xué)報(bào): 自然科學(xué)版, 1997, 28(1): 1?3.XIAO Bin, BAO Guang-shu, ZHAO Qiu-mei, et al. The improvement in complex signal of GPR[J]. Journal of Central South University Technology: Natural Science, 1997, 28(1): 1?3.

[21]楊秋芬. 基于EMD分解的探地雷達(dá)信號(hào)瞬時(shí)頻率分析[J]. 煤田地質(zhì)與勘探, 2009, 37(4): 64?67.YANG Qiu-fen. The instantaneous frequency analysis of GPR data using empirical mode decomposition[J]. Coal Geology &Exploration, 2009, 37(4): 64?67.