李秀坤，李婷婷，馬 濤

(哈爾濱工程大學(xué) 水聲技術(shù)實(shí)驗(yàn)室，哈爾濱 150001)

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

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

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

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

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

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

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

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

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

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

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

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

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

1.2 EMD方法的終止條件

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

1.3 趨勢(shì)項(xiàng)的確定

緩慢趨勢(shì)項(xiàng)的估計(jì)可能是分解后的余項(xiàng)和幾個(gè)低頻模式分量的和，即認(rèn)為從第D個(gè)基本模式分量到最后一個(gè)基本模式分量和余項(xiàng)的和。由此可以得到信號(hào)的緩慢趨勢(shì)項(xiàng)為:

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

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

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

2 實(shí)測(cè)信號(hào)的數(shù)據(jù)處理

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

圖1 實(shí)測(cè)信號(hào)的時(shí)域波形和功率譜Fig.1 Time-domain waveforms and power spectrum of achual signal

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

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

圖2 最小二乘法去趨勢(shì)項(xiàng)Fig.2 Tendency removing by least square method

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

圖3 EMD分解的IMF分量和余項(xiàng)Fig.3 IMF and residue of the EMD

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

圖5 期望信號(hào)和所剔除的趨勢(shì)項(xiàng)的時(shí)域圖Fig.5 Time-domain waveforms of the expect signal and the tendency

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

3 結(jié)論

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

［1］何祚鏞，趙玉芳.聲學(xué)理論基礎(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ì)項(xiàng)對(duì)時(shí)域參數(shù)識(shí)別的影響及消除［J］.震動(dòng)、測(cè)試與診斷.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］王廣斌，劉義倫，金曉宏，等.基于最小二乘原理的趨勢(shì)項(xiàng)處理及其 MATLAB的實(shí)現(xiàn)［J］.有色設(shè)備，2005，1003－8884，05 －04 －05.

［8］王宏禹.非平穩(wěn)隨即信號(hào)分析與處理［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.