基于EMD-1（1/2）維譜的艦船輻射噪聲調(diào)制特征提取

許勁峰 鄭 威

（江蘇科技大學 鎮(zhèn)江 212003）

1 引言

艦船輻射噪聲主要包括機械噪聲、螺旋槳噪聲和水動力噪聲，其頻譜往往表現(xiàn)為寬帶連續(xù)譜和線譜的疊加，不同類型的艦船螺旋槳在不同條件下航行，使得它們在螺旋槳噪聲方面形成了差異，可用于艦船檢測和分類識別，得到了廣泛的研究。在艦船輻射噪聲中，螺旋槳節(jié)拍對寬帶輻射噪聲存在著明顯的振幅調(diào)制。這種幅度調(diào)制信號即包絡信號，攜帶有許多重要的信息，其調(diào)制頻率和調(diào)制深度與螺旋槳轉動的軸頻、葉片頻及航速有關。對調(diào)制譜的準確提取，可有效地運用到推算艦船的速度和類型，具有重要意義［1～4］。

經(jīng)過專家學者們幾十年的研究，對艦船輻射噪聲研究取得了一系列的進展，獲得了許多顯著成果。通過分析艦船噪聲產(chǎn)生的原理，對各部分機械推進器和其他裝置進行了深入研究，對各噪聲的縱向分布和噪聲的非高斯性做了大量工作。陶篤純等分析了艦船輻射噪聲中調(diào)制包絡的物理特性，并通過研究指出了艦船輻射噪聲中幾種常見的調(diào)制類型。將艦船輻射噪聲中包絡調(diào)制作為具有相同形狀和重復周期、隨機幅度的脈沖序列來處理，并在這種理論下研究計算了各種類型調(diào)制包絡的功率譜密度函數(shù)［5］。史廣智，胡均川研究了大量的艦船輻射噪聲樣本數(shù)據(jù)，通過艦船輻射噪聲調(diào)制譜諧波族的結構特點建立了螺旋槳空化噪聲信號模型，并根據(jù)空化噪聲模型推出了諧波族結構特性表達式。然而，通過該模型建立的諧波族結構特性容易受到艦船噪聲的信噪比、平穩(wěn)性等因素的影響［6］。文獻［7～8］在對線譜進行研究時，討論了基于功率譜的方法提取線譜的主要內(nèi)容以及優(yōu)缺點。文獻［9］利用離散傅里葉變換和Welch算法對艦船輻射噪聲信號進行了線譜提取，主要是用Welch算法對信號比較密集頻段進行提取，但這樣提取的線譜信息不能代表全局的艦船線譜信息，對后續(xù)研究的結果影響較大。文獻［10］利用輻射噪聲中線譜的低起伏特性，使用了短時傅里葉分析法（STFT）對輻射噪聲特征線譜進行了提取。文獻［11～12］通過對艦船噪聲的摸底試驗，較為詳盡地研究了它的組成形式，并且根據(jù)線譜的相關文獻資料及其參數(shù)數(shù)據(jù)的研究，建立了艦船輻射噪聲線譜的提取模型。緊接著對線譜的平穩(wěn)性以及特有性質(zhì)進行了論證，在提取艦船線譜特征之后，通過神經(jīng)網(wǎng)絡等方法對艦船特征線譜的類型進行劃分。

文獻［13～15］通過在線譜圖的分析中，提出了一種較為優(yōu)越的優(yōu)化算法，其主要是利用代價函數(shù)對線譜圖進行分析，并通過這種函數(shù)直接對線譜進行提取。在研究過程中，為了避免檢測效率低下、虛警概率較高的情況，文中提出了一種新的智能的方式，即利用對雙門限檢測流程的模仿，把線譜判別、檢測以及追蹤相結合。這樣就能夠在比較低的信噪比中完成對艦船特征線譜的自動檢測。

文獻［16～17］根據(jù)DEMON譜法研究了艦船噪聲的頻譜以及頻率特征，因為調(diào)制線譜的方位和軸頻與葉片頻率相對應，所以對解調(diào)出的調(diào)制線譜進行譜分析，得到軸頻和葉片頻調(diào)制譜，進一步分析艦船線譜中的葉片數(shù)量以及螺旋槳的速度，為被動聲納目標檢測和分類提供了有力的依據(jù)。文獻［18］通過對噪聲特性的研究，結合經(jīng)典功率譜和DEMON譜分析法對輻射噪聲線譜與連續(xù)譜分離、去除虛警及歸并線譜，有效地提取了特征線譜。文獻［19］綜合運用改進的最大公約數(shù)算法和余數(shù)門限算法提取DEMON譜中的軸頻和葉頻，解決了傳統(tǒng)最大公約數(shù)算法提取軸頻葉頻誤差較大的問題。

傳統(tǒng)的DEMON譜分析在提取調(diào)制譜時，需要人為將信號分成幾個頻段，然后進行寬帶或窄帶解調(diào)，最后濾波進行傅里葉計算，結果穩(wěn)定性差、誤差較大。EMD方法依據(jù)數(shù)據(jù)自身的時間尺度特征來進行信號分解，無須預先設定任何基函數(shù)，這比小波分解，DEMON譜分析的寬帶、窄帶解調(diào)更具靈活性和適用性。1（1/2）維譜具有抑制高斯噪聲，分析諧波信號時，信號的基頻分量可得到加強，非諧波項可被清除等性質(zhì)。本文就是結合EMD分解和1（1/2）維譜分析的優(yōu)點，提取艦船噪聲中的有效信息。

2 原理分析

2.1 EMD原理

EMD方法是美籍華人N.E.Huang等于1998年提出的，目的是通過對非線性非平穩(wěn)信號的分解獲得一系列表征信號特征時間尺度的IMF，使得各個IMF是窄帶信號，可以進行HS分析。IMF要滿足兩個條件：1）整個數(shù)據(jù)集的極大極小值數(shù)目與過零點數(shù)目相等或最多相差一個；2）數(shù)據(jù)集的任意點上，由極大值確定的包絡與由極小值確定的包絡的均值始終為零。這兩個條件實際上使得分解得到的IMF是窄帶信號。而且EMD分解基于下面的假設：1）信號至少有兩個極值，一個極大值和一個極小值；2）信號特征時間尺度是由極值間的時間間隔確定的；3）如果數(shù)據(jù)中缺乏極值點，但存在缺陷點，可以通過微分、分解、再積分的方法獲得IMF［20］。

經(jīng)驗模態(tài)分解也稱篩選過程［21］（The sifting processing），分解過程具體如下：

1）找出原始信號x(t)上的所有極大值和極小值點，然后分別對所有的極大值點和極小值點用三次樣條函數(shù)進行連接，擬合出原數(shù)據(jù)序列上的包絡線，上包絡xmax(t)和下包絡線xmin(t)。此時，這兩條包絡線就包含了原信號的所有數(shù)據(jù)，對上下包絡線求均值，得到一條均值線m1(t)。

2）再用原信號x(t)減去上下包絡均值m1(t)得到一個新的數(shù)據(jù)序列：

3）對于不同的信號x(t)，可能存在著負的局部極大值和正的極小值點，這時它并不滿足IMF的條件，需要把h1(t)作為原信號繼續(xù)進行“篩選”，重復上述步驟，得到：

上式中，m11(t)是h1(t)的上下包絡線的均值，如果h11(t)不是固有模態(tài)函數(shù)，則繼續(xù)進行篩選，重復上述處理過程k次，直到得到的h1k(t)符合IMF的要求：

h1k(t)是不是一個IMF，需要由篩選過程終止準則來確定，終止準則是利用兩個連續(xù)處理結果之間的標準差Sd的值來作為判斷是否終止的依據(jù)，Sd的定義如下：

上式中，h1(k-1)(t)和h1k(t)是在篩選IMF過程中，兩個連續(xù)處理結果的時間序列。Huang建議Sd的典型取值范圍在0.2～0.3之間。同時，Terradas指出將Sd的取值范圍定為0.05～0.3之間會產(chǎn)生更好的結果，通過對太陽日冕數(shù)據(jù)的研究，發(fā)現(xiàn)其第一個IMF（高頻噪聲項）和趨勢項基本保持不變，具有一定的魯棒性［24］。

4）當h1k(t)滿足Sd的要求時，則h1k(t)為第一階固有模態(tài)函數(shù)，此時的h1k(t)就是最高頻率分量c1(t)，即

從x(t)中分離出c1(t)分量，得到一個差值信號r1(t)，即

將r1(t)看作是一個新的信號重復上述經(jīng)驗模態(tài)分解的過程，經(jīng)過多次運算得到全部剩余分量的ri(t)：

當滿足條件：小于預定的誤差或成為一單函數(shù)，即不可能再從中提取IMF時，就終止分解過程。

信號經(jīng)EMD分解后，原信號x(t)就可以由n階分量和殘余分量函數(shù)rn(t)構成，即

圖1給出了EMD算法的流程圖。

圖1 EMD算法流程圖

2.2 1（1/2）維譜的性質(zhì)與算法

2.2.1 （1/2）維譜的定義

對于隨機變量x(t)，它的三階累積量C3x(τ1,τ2)的對角切片為C3x(τ , τ)(τ1= τ2=τ)，定義該對角切片的傅立葉變換為隨機變量 x(t)的1（1/2）維譜C(w)。

上式也可以通過頻域計算，如下：

其中，X(w)為x(t)的傅里葉變換，X*(w)為X(w)的復共軛。

離散信號的計算方法如下：

設原始信號為{x1,x2,…,xN=KM}共K段，每段的長度是M，計算1（1/2）維譜的算法如下：

1）對每段數(shù)據(jù)去直流；

2）分別計算每段數(shù)據(jù)的三階累積量；

式中i=1,2,…,K，s1=max(0 ,τ)，s2=min(M-1,M-1-τ)

3）對每段數(shù)據(jù)的c(i)(τ)取平均，得

這種變化的根本原因是在恒溫、恒定表面面積條件下，對表面層施加壓力，液體體積會發(fā)生變化，導致液體表面層體積變化(見圖6)，壓力升高會使液體的表面張力數(shù)值下降。隨著體系壓力的升高，烴與CO2的接觸面積逐漸減小，界面張力逐漸降低。烴組分界面張力降低的幅度差，表現(xiàn)為CO2與原油的不同烴組分動態(tài)混相過程的難易程度的差別。烴組分最外層不斷被CO2高密度流體所飽和并不斷溶解到CO2中，輕烴不斷地補充到混相界面層，混相界面層為促進混相的過渡帶，烴組分碳數(shù)越小，混相層中的組分溶解速度越快，最終越容易達到完全混相。

4）對 c?(τ)做一維傅里葉變換，得到信號的1（1/2）維譜。

在對非線性相位耦合使用1（1/2）維譜分析時，需要適當?shù)慕o離散信號序列進行分段，M和K的值需要折中選取，可以獲得比較小的方差。

2.2.2 1（1/2）維譜的性質(zhì)與仿真

1）1（1/2）維譜具有以下性質(zhì)：

性質(zhì)1 設x(t)為零均值、基頻為w0的n次實諧波信號，在幅值a相等，相位為零的情況下，當|wm|＜ | wl|時，有

其 中 wm=mw0,m=±1,±2,…±n,wl=lw0,l=±1,±2,…,±n 。該性質(zhì)說明，當選用1（1/2）維譜分析諧波信號時，信號的基頻分量可以得到加強，這有利于抽取信號中較弱的基頻分量。要說明的是該性質(zhì)中的零相位假設，顯然，在實際中各諧波分量的相位不可能為零，但這不是問題，可以在計算1（1/2）維譜前對信號序列進行自相關處理，正弦信號的自相關函數(shù)仍為同頻率的正弦信號，但其相位是零，幅值為原幅值平方的二分之一。

性質(zhì)2 設n(t)是零均值的隨機噪聲，任何兩個不同時刻都互不相關，且概率密度函數(shù) f(n)為對稱分布，則有

該性質(zhì)表明當信號中混有對稱分布的隨機噪聲時，理論上也可被1（1/2）維譜完全抑制掉。

性質(zhì)3 設n(t)為零均值的高斯噪聲，則有

該性質(zhì)表明1（1/2）維譜同樣具有抑制高斯隨機噪聲的特點。

性質(zhì)4 設x(t)是諧波信號，wm、wp、wq為其中三個諧波分量，若wm≠wp+wq，則有

從該性質(zhì)可知，當信號中含有非相位耦合的諧波項時，通過1（1/2）維譜的處理，這些諧波項可以被清除掉。

2）1（1/2）維譜特性仿真分析

（1）仿真研究1（1/2）維譜的去噪能力和加強基頻分量特性。仿真信號為頻率分別為50Hz、100Hz、150Hz、200Hz、250Hz的正弦信號和高斯白噪聲的疊加，信噪比為-10dB。仿真結果如圖2、圖3、圖4所示，圖3為信號的功率譜，圖4為1（1/2）維譜，可以看出采用1（1/2）維譜分析后，信噪比提高了很多，信號的基頻分量也得到了加強。

圖2 原始信號和疊加噪聲后的信號

圖3 信號功率譜圖

圖4 信號的1（1/2）維譜圖

（2）1（1/2）維譜剔除信號中非耦合諧波分量特性。仿真信號為正弦信號疊加高斯白噪聲，信噪比-10dB，其中正弦信號為

它的頻率分別為 30Hz、50Hz、100Hz、150Hz、200Hz。圖5為它的功率譜，圖6是1（1/2）維譜。從圖中的結果可以看出，采用功率譜分析時，信號中所有頻率成分均清晰可見，而采用1（1/2）維譜分析后，只有那些為整數(shù)倍的頻率成分清晰可見，而非耦合諧波部分被剔除。

圖5 信號功率譜圖

圖6 1（1/2）維譜剔除信號中非耦合諧波分量

1（1/2）維譜的這一性質(zhì)有利于分析提取艦船輻射噪聲信號，本文采用對提取包絡信號求1（1/2）維譜來求其基頻和諧波，從而判斷分析艦船航速，螺旋槳等相關信息。

3 特征線譜提取

傳統(tǒng)的DEMON譜分析法，首先對信號進行寬帶解調(diào)，需要事先確定好頻帶范圍，然后進行解調(diào)，這樣得出的結果誤差比較大，EMD分解是依據(jù)數(shù)據(jù)自身的時間尺度特征來進行信號分解，不需要事先設定任何基函數(shù)，這比小波分解，DEMON譜分析的寬帶、窄帶解調(diào)更具靈活性和適用性，而1（1/2）維譜分析對比功率譜分析具有抑制高斯噪聲；分析諧波信號時，信號的基頻分量可得到加強，非諧波項可被清除等優(yōu)勢，因此結合EMD分解和1（1/2）維譜分析對DEMON譜分析法進行改進，分解過程如下圖7所示。

圖7 基于EMD-1（1/2）改進的DEMON譜分析法流程圖

4 數(shù)據(jù)處理與分析

本次實測數(shù)據(jù)為某大型貨船順風高速航行實測數(shù)據(jù)，為了去除不確定隨機誤差，選取了中間較為穩(wěn)定測量信號，其長度為10s，采樣頻率為10000Hz，然后對噪聲信號進行EMD分解。觀察發(fā)現(xiàn)，信號分解出的IMF包含了原始信號的絕大多數(shù)信息，而最后的殘余分量是信號的趨勢項沒有保留，分解出了六階IMF分量作進一步特征提取對象，并對此階分量作了譜分析。前三階IMF分量可以看出明顯的調(diào)制現(xiàn)象，且前三節(jié)IMF分量包含了信號的主要組成部分，所以對這三階IMF進行平方低通解調(diào)，然后進行譜分析，觀察看到了得到了輻射噪聲信號線譜信息。

首先，由于實測數(shù)據(jù)量信息龐大，對數(shù)據(jù)進行截取，得到中間比較穩(wěn)定的部分，然后對數(shù)據(jù)進行譜分析，觀察輻射噪聲信號的頻率大致分布情況，如圖8所示，信號能量主要集中在高頻段，包含了螺旋槳對寬帶噪聲的調(diào)制。

圖8 艦船輻射噪聲原始信號與功率譜圖

對截取的數(shù)據(jù)進行EMD分解，得到了前六階IMF分量，如圖9所示。然后對這6個分量進行譜分析，如圖10所示，觀察發(fā)現(xiàn)前三階分量中，有明顯的調(diào)制包絡現(xiàn)象，然后進一步譜分析，觀察線譜信息。

對EMD分解的前三階IMF分量進行譜分析，圖11是功率譜分析，圖12是1（1/2）維譜分析，對比發(fā)現(xiàn)都能看到線譜信息，但是1（1/2）維譜分析抑制噪聲效果明顯，得到的線譜比較清晰。圖13是1（1/2）維譜分析過后的局部放大圖，可以看到輻射噪聲的軸頻基頻為32Hz，且二次諧波分量和三次諧波分量信息明顯。

圖9 EMD分解圖

圖10 EMD分解功率譜圖

圖11 前三階IMF分量平方檢波功率譜圖

5 結語

針對傳統(tǒng)的DEMON譜分析艦船輻射噪聲調(diào)制特征線譜的缺點，即在提取調(diào)制譜時，效果不太理想，存在對噪聲抑制能力差、提取后的諧波特征不明顯等缺點。通過研究EMD分解和1（1/2）維譜分析的特性，對傳統(tǒng)DEMON譜分析作進一步改進，得到了基于EMD-1（1/2）維譜的艦船輻射特征提取方法。通過仿真分析和實測數(shù)據(jù)的處理，可以看出本方法在艦船輻射噪聲特征提取方面比傳統(tǒng)的DEMON譜分析更具優(yōu)勢。