許勁峰 鄭 威

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

1 引言

艦船輻射噪聲主要包括機(jī)械噪聲、螺旋槳噪聲和水動(dòng)力噪聲，其頻譜往往表現(xiàn)為寬帶連續(xù)譜和線譜的疊加，不同類(lèi)型的艦船螺旋槳在不同條件下航行，使得它們?cè)诼菪龢肼暦矫嫘纬闪瞬町悾捎糜谂灤瑱z測(cè)和分類(lèi)識(shí)別，得到了廣泛的研究。在艦船輻射噪聲中，螺旋槳節(jié)拍對(duì)寬帶輻射噪聲存在著明顯的振幅調(diào)制。這種幅度調(diào)制信號(hào)即包絡(luò)信號(hào)，攜帶有許多重要的信息，其調(diào)制頻率和調(diào)制深度與螺旋槳轉(zhuǎn)動(dòng)的軸頻、葉片頻及航速有關(guān)。對(duì)調(diào)制譜的準(zhǔn)確提取，可有效地運(yùn)用到推算艦船的速度和類(lèi)型，具有重要意義［1～4］。

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

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

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

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

2 原理分析

2.1 EMD原理

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

經(jīng)驗(yàn)?zāi)B(tài)分解也稱篩選過(guò)程［21］（The sifting processing），分解過(guò)程具體如下：

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

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

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

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

h1k(t)是不是一個(gè)IMF，需要由篩選過(guò)程終止準(zhǔn)則來(lái)確定，終止準(zhǔn)則是利用兩個(gè)連續(xù)處理結(jié)果之間的標(biāo)準(zhǔn)差Sd的值來(lái)作為判斷是否終止的依據(jù)，Sd的定義如下：

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

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

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

將r1(t)看作是一個(gè)新的信號(hào)重復(fù)上述經(jīng)驗(yàn)?zāi)B(tài)分解的過(guò)程，經(jīng)過(guò)多次運(yùn)算得到全部剩余分量的ri(t)：

當(dāng)滿足條件：小于預(yù)定的誤差或成為一單函數(shù)，即不可能再?gòu)闹刑崛MF時(shí)，就終止分解過(guò)程。

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

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

圖1 EMD算法流程圖

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

2.2.1 （1/2）維譜的定義

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

上式也可以通過(guò)頻域計(jì)算，如下：

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

離散信號(hào)的計(jì)算方法如下：

設(shè)原始信號(hào)為{x1,x2,…,xN=KM}共K段，每段的長(zhǎng)度是M，計(jì)算1（1/2）維譜的算法如下：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖2 原始信號(hào)和疊加噪聲后的信號(hào)

圖3 信號(hào)功率譜圖

圖4 信號(hào)的1（1/2）維譜圖

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

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

圖5 信號(hào)功率譜圖

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

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

3 特征線譜提取

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

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

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

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

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

圖8 艦船輻射噪聲原始信號(hào)與功率譜圖

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

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

圖9 EMD分解圖

圖10 EMD分解功率譜圖

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

5 結(jié)語(yǔ)

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