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)