0 引言

    目前，水聲監(jiān)聽是海洋調(diào)查研究的重要手段之一，監(jiān)聽信號頻帶范圍較廣，主要包括寬帶中低頻信號和窄帶高頻信號，如5～10 Hz的海面波浪噪聲；覆蓋在10～200 kHz在以上頻率的海洋哺乳動物發(fā)聲；中心頻率集中在1～200 kHz或更高的窄帶內(nèi)的商用聲吶；1 kHz以下的艦船噪聲等^[1]。在實際水聲監(jiān)聽過程中，記錄以上各類有效特征信息的同時會摻雜大量復雜噪聲信息。因此，如何最大程度實現(xiàn)未知監(jiān)聽信息分類提取是研究的關鍵。

    小波包算法因具有將信號分解到各頻段范圍的優(yōu)勢，被廣泛應用于分頻段閾值降噪、能量特征提取和水聲信號處理等方面，如鐘孟春等改進的小波包能量分段閾值降噪方法，按最優(yōu)小波包能量法分頻段進行降噪^[2]。史秋亮等基于小波包分解與能量特征提取的相關分析法，按子帶能量提取特征值^[3]。楊亞菁等最佳能量小波包技術在海洋水聲信號處理中的應用，提出基于分類距離標準小波包基能量法實現(xiàn)水聲信號的分類和標識^[4]。劉深等基于IMF能量譜的水聲信號特征提取與分類研究，通過經(jīng)驗模態(tài)分解，按本征模式分量的能量占比提取和分類^[5]。以上研究方法均實現(xiàn)信號的精細化分頻，但利用能量法進行特征提取，具有很大的局限性，容易將弱能量有效信號忽略，高能量噪聲依然保留。趙超等基于EMD和小波包能量法的信號去噪，利用EMD和小波包相結(jié)合的方式具有很高的參考價值，但僅實現(xiàn)了高斯白噪聲的去噪，未涉及其他成分噪聲和特征信號提取^[6]。李茂等基于EMD及主成分分析的缺陷超聲信號特征提取研究，利用EMD分量的能量占比選取主成分進行降維處理，實現(xiàn)缺陷超聲信號檢測，該方法只針對缺陷超聲信號，且超聲發(fā)射探頭和接收探頭的頻率已知^[7]。陳功等提出希爾伯特-黃變換在微弱被動瞬態(tài)魚聲信號中的檢測，利用經(jīng)驗模態(tài)分解實現(xiàn)魚類瞬態(tài)檢測，該方法僅針對魚類這一特定目標信號^[8]。蘇祖強等基于小波包分解與主流形識別的非線性降噪，針對軸承的振動信號進行小波包分解系數(shù)相空間主流行識別重構實現(xiàn)信號分離^[9]。以上方法針對性較強，目標較為明確。但水聲監(jiān)聽探測的聲頻帶范圍較寬，信號記錄過程中包含各個頻段的信息，有效特征信息可能分布于高、中、低不同頻段且能量占有比大小不一，其頻段分布和能量大小未知，因此，不能單純以能量高低的方法來判斷有效特征信息，更不能以單一方法提取不同頻段特征信息，否則，易造成有效頻段特征信息丟失，噪聲作為有效特征信息保留的情況。

    小波包能量法特征提取、希爾伯特-黃能量法特征提取等，容易將微弱的高頻特征段遺漏，能量高的低頻噪聲引入，也容易將噪聲頻段作為特征頻段提取，特征頻段作為噪聲忽略。采用單一方法對于特定頻段目標提取是可靠的，但并不適應于寬頻帶水聲監(jiān)聽信號多個目標特征頻段的提取。針對以上問題，采用小波包、相關系數(shù)、經(jīng)驗模態(tài)分解相結(jié)合的方式實現(xiàn)優(yōu)勢互補提取有效特征頻段。小波包能量法按最優(yōu)分解層對未知水聲監(jiān)聽信號從低頻到高頻無差別分段處理。相關系數(shù)法將小波包分解后的所有子帶進行相關性分析，提取含有效特征信息的子帶，舍棄噪聲子帶，其優(yōu)勢在于可判別信號子帶的有效性，不依靠能量占比高低來識別有效信號。經(jīng)驗模態(tài)分解將相關系數(shù)法提取的有效子帶進一步剔除噪聲，提高有效子帶的純度，其優(yōu)勢在于可將自相關系數(shù)判別的有效子帶內(nèi)含有的噪聲量進一步剔除，最終將純度較高的子帶通過小波包重構獲得特征信息頻段。本文方法相當于在小波包分解和重構的過程中引入了相關系數(shù)分析法和經(jīng)驗模態(tài)分解，頻帶信號處理范圍較寬，特定頻段聚焦分析能力較強，比單一小波包閾值處理法可靠性更高。經(jīng)仿真和實測數(shù)據(jù)處理結(jié)果表明該方法在寬帶水聲監(jiān)聽信號頻段聚焦分析與提取方面優(yōu)勢明顯。

1 信號提取的基本原理與方法

    水聲監(jiān)聽信號時域信息用以下公式表示：

式中，R(t)是水聲監(jiān)聽原始信號，s(t)代表有效信號，n(t)代表環(huán)境噪聲信息。本文主要對有效信號s(t)保護和提取，盡最大努力忽略噪聲信號n(t)，流程圖如圖1所示。

1.1 小波包基本方法

    小波包可實現(xiàn)水聲監(jiān)聽信號頻段多層次劃分，既可分解低頻，也可分解高頻，與小波變換相比，具有更高的時頻分辨率，可對寬頻信號進行更好的時頻分析^[10]。小波與小波包的子帶分解對比如圖2所示，A代表低頻部分，D代表高頻部分。  

    設尺度函數(shù)和小波函數(shù)分別為v(t)和p(t)，令U₀=v(t)，U₁=p(t)，定義函數(shù){U_n}為尺度函數(shù)v(t)的小波包，如式(2)所示，式中h_k和g_k為多分辨率分析中濾波器系數(shù)。

    所謂小波包就是一個函數(shù)族構造規(guī)范的正交基庫，此庫包含許多規(guī)范正交基，小波正交基只是其中一組，因此，小波包也是小波概念的延伸。

    (1)小波包分解。選擇最優(yōu)小波包基并確定最優(yōu)分解層，對信號進行小波包分解。小波包分解實際就是對小波包系數(shù)的分解。U(i，j)為第j層第i個小波包系數(shù)，設R(t)的小波包系數(shù)為U(0，0)，分解1層后系數(shù)為U(0，1)，U(1，1)，其分解后數(shù)據(jù)長度為U(0，0)的1/2。分解層數(shù)為N，則第N層每個節(jié)點的系數(shù)長度為U(0，0)的1/2^N。

    (2)閾值與閾值函數(shù)選取。對每個小波包分解系數(shù)，需選擇一個恰當?shù)拈撝岛烷撝岛瘮?shù)進行閾值量化處理。常用閾值方法有：固定閾值法、無偏似然閾值法、啟發(fā)式閾值、極大極小閾值法。閾值函數(shù)通常有硬閾值函數(shù)和軟閾值函數(shù)，以及另外幾種改進的閾值函數(shù)，如加權平均值函數(shù)、半軟閾值函數(shù)等^[11]。

    (3)小波包重構。對處理后的低頻和高頻系數(shù)進行數(shù)據(jù)整合后重構。

1.2 希爾伯特-黃算法

    HHT處理非平穩(wěn)信號的基本過程是：首先，利用EMD方法將給定的信號分解為若干IMF，為把復雜信號分解為簡單的單分量信號的組合，在進行EMD時，所獲得的IMF必須滿足下列兩個條件：(1)在整個信號長度上，所有極值點和過零點數(shù)目必須相等或者最多差一點；(2)任意時刻，局部極大值點的上包絡和極小值點的下包絡平均值為零，也就是對稱于時間軸。然后，對每一個IMF進行Hilbert變換，得到相應的Hilbert譜，即將每個IMF表示在聯(lián)合的時頻域中。

    EMD是使用由局部最大值和最小值分別定義的，一旦確定了極值，所有的局部最大值就會被三次樣條函數(shù)擬合曲線連接起來，局部最小值生成最低層的過程是不斷重復的^[12]。如給定一個信號x(t)，得到上下包絡的平均值m₁，判斷是否滿足IMF的兩個條件；若不滿足，重復循環(huán)n次直到滿足IMF的兩個條件。記為信號經(jīng)EMD得到的第1個IMF1分量，將IMF1從信號x(t)中分離出來得到R₁=x(t)-IMF1，這一處理過程不斷被重復提取n次，R_n=R_n+1-IMF1。最終合成公式為：

1.3 相關系數(shù)

    互相關函數(shù)是衡量兩個時間序列x(t)和y(t)在兩個不同時刻t₁，t₂之間的相似程度，通?？梢杂糜谠陂L序列中尋找一個特定的短序列。在數(shù)理統(tǒng)計中，互相關表示兩個隨機序列的相關性，設采集樣本為n，互相關記為R_n(x(t_i)，y(t_i))，如式(4)所示。R_n(x(t_i)，y(t_i))的取值在-1與+1之間，變量值為正數(shù)，表明是正相關，若為負數(shù)，表明是負相關。R_n(x(t_i)，y(t_i))絕對值越大表明相關性越強^[13]。

    自相關函數(shù)是互相關的一種特殊情況，反映信號n(t)在兩個不同時刻的關聯(lián)程度，如式(5)所示。歸一化自相關函數(shù)的公式如式(6)所示，式中R_n(0)為信號自身在0點時刻的函數(shù)值。噪聲信號在0點處存在最大自相關，其余點自相關系數(shù)迅速衰減為0^[14]。一般信號在0點處取值最大，但其他點隨著時間緩慢減小為0。

2 模擬信號仿真

    通過模擬信號y=0.5sin(2πf₁t)+sin(2πf₂t)對上述方法仿真計算，其中f₁=1 kHz，f₂=2 kHz，添加noise=randn(size(y))·0.4隨機噪聲信號，采樣頻率100 kHz。原始及加噪后信號如圖3中(a)和(b)所示。小波包分解層數(shù)過低去噪不完全，層數(shù)過高存在過分消噪和運算成本問題，因此需找出最優(yōu)分解層。原始信號已知情況下，可通過信噪比和均方根差的方式來判斷最優(yōu)分解層^[15]。利用小波對加噪信號分別采取1～5層分解與重構，每層重構信號的信噪比和均方根差如表1所示，通過數(shù)據(jù)對比得出第3層信號的信噪比最高，均方根差最低，因此，小波包最優(yōu)分解層為3。

    對模擬信號最優(yōu)3層分解的7～14節(jié)點做歸一化自相關處理，根據(jù)文獻[6]、[16]論證可知：噪聲具有隨機和非周期性質(zhì)，其自相關系數(shù)曲線在0點以外急劇衰減為零；而一般信號具有非隨機性質(zhì)，因此一般信號自相關系數(shù)曲線并沒有迅速衰減到很小值，而是隨著時間差變化慢慢衰減到0，變化規(guī)律明顯有別于噪聲信號的自相關系數(shù)變化規(guī)律。僅節(jié)點7符合一般信號自相關特性，同時節(jié)點7歸一化自相關系數(shù)變化與原始信號y成一定比例，如圖3(d)所示。其余節(jié)點均與節(jié)點8的系數(shù)類似符合噪聲信號特征如圖3(e)所示。

    節(jié)點7進行經(jīng)驗模態(tài)分解得到9個頻率從高到低的模態(tài)分量IMF₁～IMF₉和一個殘余分量Residual。按照式(3)分別計算IMF₁～IMF₉與節(jié)點7的互相關系數(shù)，IMF₁～IMF₉系數(shù)記為x₁～x₉，節(jié)點7系數(shù)記為u，互相關系數(shù)記為R(x_i，u)，如表2所示，最大互相關系數(shù)為0.971 8，按照參考文獻[6]小于最大互相關系數(shù)的1/10判斷分量的有效性，主要分量分布于x₁、x₂，將2個分量重構后數(shù)據(jù)與原節(jié)點7的互相關系數(shù)為0.996 7，通過經(jīng)驗模態(tài)分解重構處理后的互相關系數(shù)可以看出與原節(jié)點信號一致性較強。本文方法處理后數(shù)據(jù)如圖3(f)所示，通過與單純的小波包分解相比，本文方法的信噪比和均方根誤差更優(yōu)，如表3所示。

3 監(jiān)聽水聲信號去噪分析

    實測水聲監(jiān)聽信息來自布放于阿拉斯加州冰灣海域，水深位置在90 m被動水聽器所記錄的被動水聲信號，其主要貢獻者為阿拉斯加費班克大學和華盛頓大學研究人員^[17]。監(jiān)聽采樣頻率為100 kHz，從0.1 kHz到50 kHz分為64個頻段。讀取2011年某段原始音頻文件，其時頻特性標度如圖4所示。

    由圖4可以看出原始水聲監(jiān)聽信號頻率0.1 kHz～1 kHz范圍內(nèi)存在較強信號分布，其余信號相對較弱，但特征頻率信息提取過程中不能按照信號強弱判斷信號的有效性，信號未知的情況下，強信號可能為噪聲信號，弱信號也可能含有效信號。本文核心工作是在強弱不同的頻段信號中，提取可靠的監(jiān)聽信號特征頻率信號，特別是高頻弱信號。按照小波包最優(yōu)分解4層分解，節(jié)點頻率排序及強弱分布如圖5所示。

    小波包分解第四層節(jié)點系數(shù)15～30按照頻率順序由下往上進行排序，其中橫軸為時間，左邊為頻率順序，右邊為節(jié)點按頻率范圍排序。將0～50 kHz的頻率信號共分成16個頻段，由圖5可知每個節(jié)點的頻率的強弱分布，但信號的強弱并不能區(qū)分有效信號和噪聲信號的頻段分布，本文通過歸一化自相關系數(shù)法對4層所有節(jié)點頻率信息進行分析，區(qū)分監(jiān)聽信號有效特征節(jié)點和噪聲節(jié)點頻段分布情況。

    實測的水聲監(jiān)聽信號為未知信號，低頻信息和高頻信息當中都有可能存在有效信號，利用小波包對水聲監(jiān)聽信號在最優(yōu)分解層分解后，節(jié)點系數(shù)按頻率從低頻到高頻進行排序，在未知條件下，高頻節(jié)點信號可能存在有效信號，低頻節(jié)點信號也可能存在低頻噪聲，信號降噪處理過程中不能單純將高頻或低頻節(jié)點作為噪聲直接舍棄。

    根據(jù)文獻[6]、[16]的驗證分析得到的噪聲與一般信號的自相關系數(shù)變化規(guī)律。對圖4原始水聲監(jiān)聽信號和15～30的節(jié)點做歸一化自相關處理，其中節(jié)點15、19、23、25、27共5個節(jié)點均符合一般信號的自相關特性^[18]，如圖6所示。節(jié)點15與原始信號的自相關系數(shù)一致性較強，該頻段信號有效成分較多。節(jié)點19和節(jié)點27在零點處的最大值偏小，與原始信號自相關系數(shù)相類似，含部分有效成分。節(jié)點25與原始信號自相關系數(shù)一致性較差，包含較少有效信號，但其符合一般信號特征。節(jié)點23是高頻部分，能量較小，通過自相關系數(shù)看出，該節(jié)點在零點處最大值較大，且與原始信號自相關系數(shù)相類似，有效成分含量高，該部分是在常規(guī)數(shù)據(jù)提取中最容易忽略的。

    由圖5節(jié)點系數(shù)頻率分布圖可看出，節(jié)點16、18能量相對較高，但能量高并不代表信號的有效性，也可能只是噪聲能量高。通過歸一化自相關系數(shù)對能量占有較高的16、18節(jié)點進行分析，由圖7歸一化自相關系數(shù)可看出其符合噪聲信號在t=0時刻迅速衰減的特點。剩余節(jié)點的歸一化自相關系數(shù)均與節(jié)點16、18節(jié)點的情況基本相似，可按噪聲節(jié)點舍棄，由于篇幅限制不再將其一一列出。

    利用小波包算法對監(jiān)聽數(shù)據(jù)進行節(jié)點主成分初判以后，通過希爾伯特算法分別對含有效成分的15、19、23、25、27節(jié)點做進一步消噪提取。以節(jié)點15為例，對其進行經(jīng)驗模態(tài)分解得到頻率從高到低的模態(tài)分量IMF1～IMF10和一個殘余分量Residual。按照式(3)分別計算IMF1～IMF10與節(jié)點15的互相關系數(shù)，IMF1～IMF10系數(shù)記為c₁～c₁₀，節(jié)點15系數(shù)記為s，互相關系數(shù)記為R(c_i，s)，如表4所示，最大互相關系數(shù)為0.771 9，按照參考文獻[6]小于最大互相關系數(shù)的1/10判斷分量的有效性，其中分量c₅、c₆、c₇、c₈、c₉、c₁₀與原始節(jié)點15的互相關系數(shù)較小，主要分量分布于c₁、c₂、c₃、c₄，將4個分量重構后數(shù)據(jù)與原節(jié)點15的互相關系數(shù)為0.994 5，通過相關系數(shù)可看出與原節(jié)點信號一致性較強。利用上述方法分別對剩余19、23、25、27節(jié)點系數(shù)實現(xiàn)經(jīng)驗模態(tài)的分解與重構，最終通過小波包將含有效成分的5個節(jié)點系數(shù)15、19、23、25、27重構，獲得監(jiān)聽的特征頻率信息信號，通過本文方法實現(xiàn)寬頻帶信號分頻段提取，即提取到了低頻強信號頻段，也提取到了高頻弱信號特征頻段，如圖8所示。

4 驗證分析

    小波包分解與重構過程中引入自相關系數(shù)和經(jīng)驗模態(tài)分解的分析，旨在聚焦寬頻帶水聲監(jiān)聽信號的各特征頻段信息提取，盡最大可能保留有效聲頻信息，排除噪聲干擾，進一步提高分辨率，精細和優(yōu)化提取結(jié)果，凸顯高頻弱信號提取能力。由圖8可以看出本文方法對寬頻帶水聲監(jiān)聽信號既提取到能量較高、頻率相對低的信號，也提取到了能量較低、頻率較高信號。特征頻段提取的有效性就是盡最大可能保留和凸顯有效信號，弱化噪聲信號，圖9(a)為沒有經(jīng)過數(shù)據(jù)處理的原始水聲監(jiān)聽信號成分信息時頻分析結(jié)果，圖9(b)為經(jīng)過小波包閾值降噪處理后時頻分析結(jié)果，圖9(c)為采用本文方法處理后的時頻分析結(jié)果。圖9(b)與圖9(a)相比，信號頻段提取較為明顯，提取頻段大致為0～10 kHz、15～25 kHz、26～36 kHz、38～48 kHz共4段，特征頻段有效信號進一步凸顯。圖9(c)與圖9(b)相比看出：圖9(c)有效特征頻段層為0～5 kHz、7～10 kHz、15～18 kHz、20～30 kHz、33～43 kHz、45～50 kHz共6段，比圖9(b)多分化2個頻段，分辨率進一步提高，高頻段的提取相對突出，45～50 kHz范圍高頻段弱信號提取較為明顯。圖9(b)的小波包閾值方法按照平均10 kHz的頻段間隔進行消噪提取，并受節(jié)點頻段能量大小約束，因此，存在一定局限性。而本文方法引入了自相關系數(shù)提高了有效子帶信號判斷準確率，并通過經(jīng)驗模特分解對有效子帶內(nèi)的噪聲進一步分離，使信號的提取更精細化。由對比分析圖可看出本文方法優(yōu)勢更加明顯。

5 結(jié)論

    本文方法集中小波包優(yōu)化寬帶信號頻段劃分、自相關系數(shù)判別子帶信號有效性、經(jīng)驗模態(tài)分解進一步分離子帶中有效信號和噪聲這三方面的優(yōu)勢，針對含有多頻段特征信息的寬帶水聲監(jiān)聽信號實現(xiàn)有效頻段特征信息提取，解決依靠能量占有比判斷信號有效性的弊端，以及單一方法提取寬頻帶信號的局限性，進一步提高寬帶水聲信號特征頻段提取的分辨率，優(yōu)化了高頻弱信號的分析分辨能力。經(jīng)仿真分析和實測水聲監(jiān)聽數(shù)據(jù)處理，本文方法優(yōu)勢明顯，既可剔除摻雜在高能量有效信號頻段內(nèi)的低頻噪聲，也可提取高頻弱信號頻段的有效信息。該方法可在寬頻帶水聲監(jiān)聽信號處理方面推廣和應用，并為后續(xù)水下目標識別處理提供可靠樣本數(shù)據(jù)。

參考文獻

[1] 張國勝，顧曉曉，邢彬彬，等.海洋環(huán)境噪聲的分類及其對海洋動物的影響[J].大連海洋大學學報，2012，27(1)：89-94.

[2] 鐘孟春，張春林，李華，等.改進的小波包能量分段閾值降噪方法[J].計算機工程與應用，2015，51(5)：204-207.

[3] 史秋亮，林江.基于小波包分解與能量特征提取的相關分析法[J].聲學與電子工程，2010(4)：18-20，24.

[4] 楊亞菁，彭宏.最佳能量小波包技術在海洋水聲信號處理中的應用[J].計算機應用與軟件，2005(4)：21-22，55.

[5] 劉深，張小薊，牛奕龍，等.基于IMF能量譜的水聲信號特征提取與分類[J].計算機工程與應用，2014，50(3)：203-206，226.

[6] 趙超，楊慶東.基于EMD和小波包能量法的信號去噪[J].北京信息科技大學學報(自然科學版)，2019，34(1)：75-79.

[7] 李茂，楊錄，張艷花.基于EMD及主成分分析的缺陷超聲信號特征提取研究[J].中國測試，2018，44(2)：118-121，133.

[8] 陳功，王平波，鮑玉軍，等.希爾伯特-黃變換在微弱被動瞬態(tài)魚聲信號中的檢測[J].海洋科學，2016，40(10)：91-96.

[9] 蘇祖強，蕭紅，張毅，等.基于小波包分解與主流形識別的非線性降噪[J].儀器儀表學報，2016，37(9)：1954-1961.

[10] 宋保業(yè)，徐繼偉，許琳.基于小波包變換-主元分析-神經(jīng)網(wǎng)絡算法的多電平逆變器故障診斷[J].山東科技大學學報(自然科學版)，2019，38(1)：111-120.

[11] 閻妍，行鴻彥.基于小波包多閾值處理的海雜波去噪方法[J].電子測量與儀器學報，2018，32(8)：172-178.

[12] 譚秋衡，吳量，李波.基于EMD及非平穩(wěn)性度量的趨勢噪聲分解方法[J].數(shù)學物理學報，2016，36(4)：783-794.

[13] 周進，卜英勇，羅柏文.互相關函數(shù)法在海底微地形測量中的應用研究[J].測繪科學，2009，34(1)：178-179.

[14] 楊濤，樂友喜，曾賢德，等.基于自相關函數(shù)特性的CEEMD全局閾值去噪方法研究[J].地球物理學進展，2018，33(4)：1622-1628.

[15] 王維，張英堂，任國全.小波閾值降噪算法中最優(yōu)分解層數(shù)的自適應確定及仿真[J].儀器儀表學報，2009，30(3)：526-530.

[16] 余發(fā)軍，周鳳星.基于EEMD和自相關函數(shù)特性的自適應降噪方法[J].計算機應用研究，2015，32(1)：206-209.

[17] Passive acoustic recording data from the Alaskan Beaufort Sea[DB/OL].https：/data.eol.ucar.edu/dataset/106.437.

[18] 王佳飛，關添，姜宇程，等.主動噪聲控制平臺的FPGA實現(xiàn)[J].電子技術應用，2018，44(2)：59-61，65.

作者信息：

趙  杰1，2，3，楊俊賢1，2，3，惠  力1，2，3，王  志1，2，3，初士博1，2，3

(1.齊魯工業(yè)大學(山東省科學院)，山東省科學院海洋儀器儀表研究所，山東 青島266061；

2.山東省海洋環(huán)境監(jiān)測技術重點實驗室，山東 青島266061；3.國家海洋監(jiān)測設備工程技術研究中心，山東 青島266061)