0 引言

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

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

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

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

    水聲監聽信號時域信息用以下公式表示：

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

1.1 小波包基本方法

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

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

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

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

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

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

1.2 希爾伯特-黃算法

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

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

1.3 相關系數

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

    自相關函數是互相關的一種特殊情況，反映信號n(t)在兩個不同時刻的關聯程度，如式(5)所示。歸一化自相關函數的公式如式(6)所示，式中R_n(0)為信號自身在0點時刻的函數值。噪聲信號在0點處存在最大自相關，其余點自相關系數迅速衰減為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)所示。小波包分解層數過低去噪不完全，層數過高存在過分消噪和運算成本問題，因此需找出最優分解層。原始信號已知情況下，可通過信噪比和均方根差的方式來判斷最優分解層^[15]。利用小波對加噪信號分別采取1～5層分解與重構，每層重構信號的信噪比和均方根差如表1所示，通過數據對比得出第3層信號的信噪比最高，均方根差最低，因此，小波包最優分解層為3。

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

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

3 監聽水聲信號去噪分析

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

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

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

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

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

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

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

4 驗證分析

    小波包分解與重構過程中引入自相關系數和經驗模態分解的分析，旨在聚焦寬頻帶水聲監聽信號的各特征頻段信息提取，盡最大可能保留有效聲頻信息，排除噪聲干擾，進一步提高分辨率，精細和優化提取結果，凸顯高頻弱信號提取能力。由圖8可以看出本文方法對寬頻帶水聲監聽信號既提取到能量較高、頻率相對低的信號，也提取到了能量較低、頻率較高信號。特征頻段提取的有效性就是盡最大可能保留和凸顯有效信號，弱化噪聲信號，圖9(a)為沒有經過數據處理的原始水聲監聽信號成分信息時頻分析結果，圖9(b)為經過小波包閾值降噪處理后時頻分析結果，圖9(c)為采用本文方法處理后的時頻分析結果。圖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的頻段間隔進行消噪提取，并受節點頻段能量大小約束，因此，存在一定局限性。而本文方法引入了自相關系數提高了有效子帶信號判斷準確率，并通過經驗模特分解對有效子帶內的噪聲進一步分離，使信號的提取更精細化。由對比分析圖可看出本文方法優勢更加明顯。

5 結論

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

參考文獻

[1] 張國勝，顧曉曉，邢彬彬，等.海洋環境噪聲的分類及其對海洋動物的影響[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] 陳功，王平波，鮑玉軍，等.希爾伯特-黃變換在微弱被動瞬態魚聲信號中的檢測[J].海洋科學，2016，40(10)：91-96.

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

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

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

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

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

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

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

[16] 余發軍，周鳳星.基于EEMD和自相關函數特性的自適應降噪方法[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實現[J].電子技術應用，2018，44(2)：59-61，65.

作者信息：

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

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

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