李秀坤，李婷婷，馬 濤

(哈爾濱工程大學 水聲技術(shù)實驗室，哈爾濱 150001)

水聲傳播理論發(fā)展的近百年過程中，主要關(guān)注的是聲波在海水介質(zhì)中的傳播。根據(jù)淺海波導聲傳播的理論，簡正波在波導中傳播的頻率低于截止頻率時不能在水中傳播［1］。尤其在夏季淺海中普遍存在負梯度水文現(xiàn)象，水中傳播的聲波將會急劇地折射向海底，此時水中的聲波更易于通過海底傳播，而不是通過水中傳播。近幾十年來人們對海底地質(zhì)聲學特性的認識不斷深入，推動了對低頻和極低頻聲傳播研究的不斷發(fā)展［2－4］。基于陸地的探測手段可以不受海洋水文環(huán)境的限制，全天候的對潮汐引起的地殼運動、海洋環(huán)境參數(shù)進行測量、探測水中目標。

理想情況下，測量系統(tǒng)需要一個絕熱、隔振的環(huán)境，以避免環(huán)境因素如:振動、溫度、濕度、氣壓等對測量的影響。但在實際的測量環(huán)境中，接收系統(tǒng)很難達到絕熱真空處理，并且對空氣擾動和低頻振動極其敏感，這些隨機干擾將直接影響接收系統(tǒng)的接收效果，所以在接收信號中會存在隨機的強低頻干擾，這類干擾稱為趨勢項干擾。

趨勢項的存在，會使時域中的相關(guān)分析和頻率中的功率譜產(chǎn)生大的誤差，甚至使低頻譜完全失去真實性［5］?，F(xiàn)有的趨勢項消除方法有幾種［6－8］，如平均斜率法、差分法、低通濾波法及最小二乘擬合方法等。這些方法通常需要預(yù)先假定信號中趨勢項的類型，如線性趨勢或指數(shù)趨勢等等，不適用于具有復雜變化趨勢或隨即變化趨勢的信號，因此不具有普遍的適用性。經(jīng)驗?zāi)B(tài)分解(EMD)方法是一種全新的處理方法，該方法依據(jù)信號本身固有的特征自然的分解信號，無需設(shè)置先驗的分解基函數(shù)，因而適用于處理各種不同類型的數(shù)據(jù)。

1 經(jīng)驗?zāi)B(tài)分解方法

希爾伯特－黃變換(HHT:Hilbert-Huang Transform)是美國工程院院士Huang及其合作者［9］在20世紀末提出的一種新的非線性、非平穩(wěn)時間序列分析方法。它基于經(jīng)驗?zāi)B(tài)分解，將信號分解為一系列固有模態(tài)函數(shù)(IMF:Intrinsic Mode Function)之和，對每一個固有模態(tài)函數(shù)進行Hilbert變換，得到其瞬時頻率，并將其以時間—頻率為坐標的平面表示出來，得到時間－頻率－能量的分布，稱為Hilbert譜。從理論上分析，HHT不僅解決了多輻射聲源時小波分析的分辨率問題和對不同聲源信號的適應(yīng)性問題，而且解決了Winger-Ville分布的多分量信號交叉項問題。本文采用EMD方法將實測信號分解成一系列IMF，其結(jié)果是將信號中不同尺度的波動或趨勢逐級分解，產(chǎn)生一系列具有不同特征尺度的數(shù)據(jù)序列，達到分離趨勢項的目的。

1.1 經(jīng)驗?zāi)B(tài)分解

EMD分解方法基于以下幾點假設(shè):(1)信號至少有兩個極值點，一個極大值和一個極小值，或者極大值或極小值數(shù)目比零點的數(shù)目多2個(或2個以上);(2)信號的特征時間尺度是由極值點之間的時間間隔確定的;(3)如果數(shù)據(jù)中缺乏極值點，但存在奇異點，可以通過一次或多次差分來求出極值點。

EMD分解的具體處理方法是，找到數(shù)據(jù)x(t)的極大值點集和極小值點集后，用插值法擬合x(t)的上、下兩條包絡(luò)線。計算兩條包絡(luò)線的平均值，記為m1(t)。計算得到x(t)與m1(t)的差，記為h1(t)，即:

一般來講，h1(t)仍然不是一個IMF分量，為此需要對它重復上述處理過程。即將h1(t)視為新的數(shù)據(jù)序列，擬合其上、下包絡(luò)，得到兩條包絡(luò)線的平均值m11(t)，并計算h1(t)與m11(t)的差值，記為h11(t)，即:

重復以上操作i次，直到h1i(t)滿足IMF的條件為止，至此得到了信號x(t)的第一個IMF分量，記為:

從數(shù)據(jù)中分離出imf1，即:

若剩余量r1(t)仍含有較長的周期成分，則將其看作新的數(shù)據(jù)，并重復以上步驟，分解出新的IMF:

當剩余量變?yōu)橐粋€常量，或者一個單調(diào)函數(shù)，再或者一個有且僅有一個極點的函數(shù)時，數(shù)據(jù)篩分結(jié)束。此時的rn(t)稱為余項，每一個篩分出的IMF分量與前一個篩分出的IMF分量相比，含有較低的頻率特性。最后，原始的數(shù)據(jù)序列即可由這些IMF分量以及一個均值或趨勢項表示:

1.2 EMD方法的終止條件

文中采用簡單收斂準則，即只要信號的極值點和過零點的數(shù)目相等時，篩選過程就終止準則。這種停止準則很簡單，同時考慮了IMF的定義，更加趨于合理化。通過仿真研究，在實際中很難達到信號的極值點和過零點的數(shù)目相等。因此，本文在使用此準則時，將終止條件設(shè)置為信號的極值點與過零點的數(shù)目相等或至多相差1，即可停止篩選。并設(shè)置最多循環(huán)1200次如不滿足簡單收斂準則，則強制退出。

1.3 趨勢項的確定

緩慢趨勢項的估計可能是分解后的余項和幾個低頻模式分量的和，即認為從第D個基本模式分量到最后一個基本模式分量和余項的和。由此可以得到信號的緩慢趨勢項為:

其中，n和rn(t)分別是經(jīng)EMD分解得到的基本模式分量的個數(shù)和余項。

為了確定D值的大小，首先定義前d個基本模式分量的和為:

根據(jù)EMD分解原理可知，每個模式分量的均值都應(yīng)該為零。但是由于一個IMF可能包含非整數(shù)個周期，所以在實際應(yīng)用中Rimf(d)可能是非零的。并且隨著分解階數(shù)i的增大，IMF的周期也將增大，隨之由于非整數(shù)周期引起的前i個IMF分量絕對值的和也單調(diào)增加。因此，可以觀察Rimf(d)隨分解階數(shù)i的變化情況，在d≥D后明顯地偏離零點便可以認為此時的IMF進入了緩慢變化的狀態(tài)。取使Rimf(d)明顯偏離零點時的d值作為D，然后用式(7)進行重構(gòu)，就可以獲得所求的緩變趨勢項。

2 實測信號的數(shù)據(jù)處理

測量系統(tǒng)在進行接收信號，地聲測量系統(tǒng)采集的數(shù)據(jù)的時域和頻域如圖1。觀察實測數(shù)據(jù)的時域和頻域上的特性不難看出，數(shù)據(jù)中存在極強的趨勢項干擾，并且趨勢項的形式非常復雜，而感興趣的信號部分則相對較弱。為了實現(xiàn)強干擾的抑制和提高弱信號的檢測與識別能力，便需要進行強干擾的分離與去除。

圖1 實測信號的時域波形和功率譜Fig.1 Time-domain waveforms and power spectrum of achual signal

圖2為采用最小二乘法擬合指數(shù)型的趨勢項所構(gòu)造得到的對于趨勢項的近似估計和去除趨勢項后信號的時域結(jié)果。采用最小二乘法只能近似擬合數(shù)據(jù)中趨勢項，對于復雜形式的趨勢項不能夠很好的擬合出來。并且，這種趨勢項的構(gòu)造方法并不是按照信號組成成分的意義進行構(gòu)造的，容易破壞有用信號成分或是沒有完全去除趨勢項部分。因此，利用這種方法去除趨勢項具有很大的局限性。

現(xiàn)采用EMD方法將實測信號分解成一系列的IMF，其結(jié)果是將信號中不同尺度的波動或趨勢逐級分解，產(chǎn)生一系列具有不同特征尺度的數(shù)據(jù)序列，從而達到分離強干擾的目的。對實測數(shù)據(jù)經(jīng)過EMD分解得到的12個IMF分量和一個余項，如圖3。Rimf(d)隨分解階數(shù)i的變化情況如圖4所示。

圖2 最小二乘法去趨勢項Fig.2 Tendency removing by least square method

按照文中所指出的趨勢項判斷方法分析，當d≥4時，Rimf(d)的值明顯地偏離了零點。因此取D=4，即利用第5到第12個基本模式分量和余項根據(jù)式(7)重構(gòu)就可以得到信號的緩變趨勢項。處理結(jié)果如圖5所示，其中上圖為消除趨勢項后的信號的時域圖，圖5為所去除趨勢項的時域圖。

圖3 EMD分解的IMF分量和余項Fig.3 IMF and residue of the EMD

圖4 前i個IMF分量的和Fig.4 Summation of the first i IMFs

圖5 期望信號和所剔除的趨勢項的時域圖Fig.5 Time-domain waveforms of the expect signal and the tendency

處理結(jié)果表明，利用EMD分解可以有效的提取出實測信號中的趨勢項。與基于最小二乘的趨勢項擬合方法相比，這種方法能有效地提取并分離出復雜的趨勢項成分。

3 結(jié)論

實測信號是多源混合信號，存在形式復雜的強趨勢項干擾，而真正所感興趣的信號部分則相對較弱。EMD算法簡單，無需任何的先驗知識就可以自適應(yīng)的把一個多分量信號進行分解，并且不損失信號成分。針對實測信號的特點，文中給出了根據(jù)IMF分量進行趨勢項判別與重構(gòu)的方法。實際數(shù)據(jù)的處理結(jié)果表明，文中方法可以有效地提取出復雜的趨勢項。

［1］何祚鏞，趙玉芳.聲學理論基礎(chǔ)［M］.北京:國防工業(yè)出版社，1981.

［2］Chapman N R，Chin-Bing S，King D，et al.Special issue on geoacoustic inversion in range-dependent shallow-water environments［J］.IEEE J.Oceanic Eng，2003，28(4):317－319.

［3］Chapman N R，Chin-Bing S，King D，et al.Inversion techniques workshop:Range dependent test cases［J］.IEEE J.Oceanic Eng，2003，28:320 －330.

［4］Fulford J K，King D，Chin-Bing S，et al.Benchmarking geoacoustic inversion methods using range-dependent field data［J］.IEEE J.Oceanic Eng，2004，29:3 －12.

［5］高品賢.趨勢項對時域參數(shù)識別的影響及消除［J］.震動、測試與診斷.1996，14(2):20 －26.

［6］Bendat JS， PiersolA G. Random dataanalysisand measurement procedures［J］.Second Edition.John Wiley ＆Sons，1986，12:396 －398.

［7］王廣斌，劉義倫，金曉宏，等.基于最小二乘原理的趨勢項處理及其 MATLAB的實現(xiàn)［J］.有色設(shè)備，2005，1003－8884，05 －04 －05.

［8］王宏禹.非平穩(wěn)隨即信號分析與處理［J］.北京:國防工業(yè)出版社，1999.

［9］Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition and Hilbert spectrum for nonlinear and nonstationary time series analysis［J］.Proc.Roy.Soc.London A，1998，454:903 －995.