0 引　言

舰船辐射噪声（SR-N）信号由船舶的各种振动和声源引起，包括机械噪声、螺旋桨噪声和水动力噪声^[1]。SR-N包含了关于船舶特性的多种信息，如船型、航速、吨位等。SR-N去噪是水声信号处理的基础和前提，对SR-N信号进行降噪处理，有利于进一步的特征提取、分类和检测。然而，SR-N不仅具有非线性、非平稳特性和非高斯特性，还具有混沌特性。因此，有必要寻找一种有效的SR-N降噪方法。

传统的SR-N降噪方法主要包括傅里叶变换、小波变换、局部投影和压缩感知等。上述方法虽然在一定程度上可以降低水声信号的噪声，但都有不同程度的局限性。例如，傅里叶变换方法不适用于非线性和非平稳信号，小波变换的基函数无法自适应选择，局部投影的邻域半径设置困难和压缩感知方法对于非高斯背景噪声的舰船辐射信号存在去噪效果不理想的问题。

1998年，李紫鹏等^[2]提出了经验模态分解（EMD）算法，是通过一种由数据本身驱动的自适应分解方法，解决了小波变换需事先设置基函数的问题，但是这种方法存在着严重的模态混叠问题。

针对该问题，文献[3 − 4]提出了在原始信号上加入不同高斯白噪声的集成经验模态分解方法，这是一种噪声辅助信号处理方法，可以有效衰减EMD方法中存在的模态混叠现象。但是由于集成经验模态分解在目标信号中多次添加了白噪声成分，存在引入的额外噪声部分不能够完全抵消的问题。为了避免集合经验模态分解时引入额外白噪声成分，抑制白噪声对于原始信号重构时产生的不良影响，文献[5 − 7]提出了一种改进的集合经验模态分解方法，称为互补集合经验模态分解（CEEMD），该方法很好地解决了经验模态分解和集合经验模态分解的不足，消除因添加白噪声而导致的重构结果出现误差的问题，并且能够在确保分解质量的前提下减少算法迭代所需次数，有效降低了计算成本。

熵是对系统有序程度进行衡量的指标，Shannon在1948年将熵的概念引入通信领域，从而形成了香农熵的基础理论。在香农熵的基础上，Bandt等提出了香农熵的改进算法排列熵（PE），排列熵的计算基于时间序列的排列模式，算法简单且计算复杂度低，是一种应用广泛的熵，可以通过排列熵来衡量模态的复杂度进而分离出噪声模态。而对于这些噪声模态，一些学者采用小波阈值去噪法来去除噪声并重构残差分量，然而低频IMF分量可能仍含有低频间谐波，其成分复杂并会随结构运行条件的变化而变化，导致小波阈值去噪法无法彻底消除噪声。对此，通过奇异谱分析（SSA）从噪声信号中提取特征分量^[8]，并通过这些分量对信号进行重构，实现进一步提取信号有效成分，从而提升信号降噪的效果。

为此，提出一种基于CEEMD算法，联合PE、WST和SSA的SR-N降噪新方法。

1 基本原理 1.1 CEEMD算法

该算法的基本原理是向目标信号中添加多组正负成对的白噪声，然后对每组噪声与信号进行单独的经验模态分解。接着，将得到的同阶模态进行求和并取平均，最终得到本征模态函数（IMF）。此方法不仅能够保持与集合经验模态分解相似的分解效果，还有效地解决了由额外白噪声引起的重构误差问题。互补集合经验模态分解算法的实现步骤为：

步骤1　给定原始信号为$ X $，通过向$ X $中分别加入$ K $组正负白噪声$ Y $后，得到$ 2K $个需要分解的信号，对这$ 2K $个信号进行EMD，可以得到对应的本征模态函数组合和余量。

步骤2　 将第$ i $个需要分解的信号中的第$ j $个本征模态函数记为$ {c}_{ij} $，最终CEEMD的本征模态函数可表示为：

$ {c}_{j}=\frac{1}{2K}\sum \limits_{i=1}^{2K}{c}_{ij} 。$

(1)

式中：$ {c}_{j} $为最终的本征模态函数。

步骤3　将原始信号表示为：

$ X=\sum \limits_{j=1}^{L}{c}_{j}(t)+r(t) 。$

(2)

式中：$ L $为本征模态函数的数量；$ r(t) $为余量。

1.2 排列熵

排列熵算法可以分析一维时间序列的复杂程度，与Lyapunov指数等计算复杂度的参数类似，但该算法耗时短且不易受噪声干扰。该算法步骤通过以下计算为：

步骤1　相空间重构。设原始时间序列为长度为$ N $的$ X=\left\{{x}_{1},{x}_{2},\cdots ,{x}_{N}\right\} $，为揭示序列的局部动态特征，构造重构向量为：

$ X_{i}^{}=\left[{x}_{i},{x}_{i+1},\ldots ,{x}_{i+(m+1)}\right], 1\leqslant i\leqslant N-(m-1)\tau。$

(3)

式中：$ m $为嵌入维数；$ \tau $为时间延迟。

步骤2　排列模式构建和排序。对每个向量$ X_{i}^{} $，依据其元素的大小关系定义一个排列模式$ \text{π}_{i} $，使得向量$ X_{i}^{} $里的元素从小到大排序：

$ {x}_{i+\sigma (0)\tau }\leqslant {x}_{i+\sigma (1)\tau }\leqslant \cdots \leqslant {x}_{i+\sigma (m-1)\tau }。$

(4)

式中：$ \sigma $为索引集合$ \{0,1,\ldots ,m-1\} $的一个排列，使得该排列对应分量满足升序条件。

相应的，排列模式$ \text{π}_{i} $的定义为：

$ {\text π }_{i}=(\sigma (0),\sigma (1),\ldots ,\sigma (m-1))。$

(5)

排列模式$ \text{π}_{i} $描述了每个嵌入向量中数值的相对大小关系。在嵌入维数为$ m $的条件下，总共存在最多$ m! $种可能的排列模式，记为$ \text{π}_{1},\text{π}_{2},\ldots ,\text{π}_{m!} $。

步骤3　排列熵的计算。计算每种排列出现的次数，可以得到其出现的概率为$ {p}_{j}={a}_{j}/N-m+1 $，其中$ j=1, 2, \ldots ,m! $。在获得所有排列模式的概率分布$ \left\{{p}_{1},{p}_{2}, \ldots , {p}_{m!}\right\} $之后，排列熵进行归一化后可以定义为：

$ {H}_{\mathrm{{norm}}}(m)=\frac{-\displaystyle\sum \limits_{j=1}^{m!}{p}_{j}\log \left({p}_{j}\right)}{\log (m!)},0\leqslant {H}_{\mathrm{{norm}}}(m)\leqslant 1 。$

(6)

归一化排列熵越大，表明时间序列的随机性越强；反之，越小则表示序列越具有规律性或确定性。

1.3 小波阈值降噪

设$ f(t) $表示待处理的观测信号，其表达式为：

$ f(t)=s(t)+n(t)。$

(7)

式中：$ s(t) $为待提取的目标信号；$ n(t) $为背景噪声信号。

将$ f(t) $进行离散小波变换，其表达式为：

$ {W}_{j,k}={2}^{\frac{j}{2}}\sum \limits_{t=0}^{N}f(t)\psi \left({2}^{-j}t-k\right) 。$

(8)

式中：参数$ N $为信号序列长度；$ \psi (t) $为所选用的母小波函数；系数$ {W}_{j,k} $为对应的小波变换系数。

在进行小波阈值去噪时，首先依据信号特征确定适当的小波函数类型及其分解层数，然后经过离散小波变换得到小波系数，小波系数可以分为包含有效信号的系数和包含噪声的系数。在降噪过程中，有效信号和噪声在数学特性上存在差异，差异主要表现为有效信号在小波域中往往集中在少数几个大系数上，这表明其具有较高的能量且幅值较大，这部分系数对信号重构起主要作用，噪声通常会散布在许多小系数中，其幅值较小。因此，通过小波变换后，不同类型的系数会表现出不同的幅值大小，接下来选择合适的阈值，并利用硬阈值或软阈值函数对每一层的小波系数进行量化。最后，通过重构处理后的系数就可以获得降噪后的信号波形。

1.4 奇异谱分析

奇异谱分析（SSA）是一种基于时间序列数据的信号处理方法，适用于分析具有噪声成分的信号。与小波变换不同，SSA无需预先设定基函数，因此在复杂信号分析中具有较强的适应性，具体步骤如下：

步骤1　嵌入过程。对于一个原始的序列数据，序列长度为$ N $，在选择合适的窗口大小$ L $后构建轨迹矩阵$ {\boldsymbol{X}} $，该轨迹矩阵是一个${\boldsymbol{ L}}\times {\boldsymbol{K}} $维的矩阵，具体形式为：

$ {\boldsymbol{X}}=\left[\begin{matrix}\boldsymbol{X}_{1} & \boldsymbol{X}_{2} & \cdots & \boldsymbol{X}_{K}\\ \boldsymbol{X}_{2} & \boldsymbol{X}_{3} & \cdots & \boldsymbol{X}_{K+1}\\ \boldsymbol{X}_{3} & \boldsymbol{X}_{4} & \cdots & \boldsymbol{X}_{K+2}\\ \vdots & \vdots & \ddots & \vdots \\ \boldsymbol{X}_{L} & \boldsymbol{X}_{L+1} & \cdots & \boldsymbol{X}_{N}\\ \end{matrix}\right]。$

(9)

式中：$ K=N-L+1 $。

2）奇异值分解。对轨迹矩阵$ X $进行奇异值分解，通过求解矩阵$ S=X{X}^{\mathrm{{T}}} $的特征方程得到特征值$ {\lambda }_{1}\geqslant {\lambda }_{2} \geqslant \ldots \geqslant {\lambda }_{L}\geqslant 0 $和对应的特征向量$ U=[{U}_{1}, {U}_{2}, \ldots,{U}_{L}] $，进而计算右奇异向量$ {V}_{i} $，最终将轨迹矩阵$ {\boldsymbol{X}} $表示为：

$ \boldsymbol{\boldsymbol{X}}=\sum_{i=1}^d\limits\sqrt{\lambda_i}U_iV_i^{\mathrm{T}}=\boldsymbol{\boldsymbol{X}}_1+\boldsymbol{\boldsymbol{X}}_2+...+\boldsymbol{\boldsymbol{X}}_d。$

(10)

式中：$ V_i^{\mathrm{T\mathrm{ }}} $为右奇异向量$ V_{i}$的转置；$ d $为非零特征值的个数。

步骤3　分量分组。在获得奇异值分解结果后，需要根据奇异值的大小和分布特征，结合奇异向量所呈现的周期性特征以及实际物理意义来确定分组方案。将指标集$ \left\{1,2,...,d\right\} $划分为$ m $个不相交的子集$ {I}_{1},{I}_{2},...,{I}_{m} $，则对应的分组矩阵可表示为：

$ \boldsymbol{X}_{{{I}_{k}}}=\sum \limits_{i\in {I}_{k}}\boldsymbol{X}_{i}，k=1,2,\ldots,m 。$

(11)

通常较大的奇异值及其对应的初等矩阵反映了时间序列的主要特征，最大奇异值通常对应着序列的总体趋势，而较小的奇异值对应的初等矩阵则可能代表噪声。

步骤4　重构过程。对分组后的矩阵进行对角平均操作。设$ x_{i,j,k}^{} $为第$ k $个分组矩阵$ \boldsymbol{X}_{{{I}_{k}}} $中位置$ \left(i,j\right) $的元素，则重构后的第$ k $个分量的第$ n $个值可表示为：

$ {{\tilde{X}}_{n,k}=\left\{\begin{aligned} &\frac{1}{n}\sum \limits_{i=1}^{n}x_{i,n-i+1,k}，1\leqslant n\leqslant L-1，\\& \frac{1}{L}\sum \limits_{i=1}^{L}x_{i,n-i+1,k}，L\leqslant n\leqslant K-1，\\& \frac{1}{N-n+1}\sum \limits_{i=n-K+1}^{N}x_{i,n-i+1,k}，K\leqslant n\leqslant N。\\ \end{aligned} \right.}$

(12)

原始时间序列可以表示为：

$ {X}_{n}=\sum \limits_{k=1}^{m}{\tilde{X}}_{n,k}。$

(13)

重构能有效去除原始时间序列中的噪声成分，提取关键特征、降低噪声、识别重要模式，从而提高序列预测的准确性和可靠性。

2 降噪方法和评价指标 2.1 基于CEEMD的联合降噪方法

利用CEEMD的分解特点，联合排列熵、小波阈值去噪及奇异谱分析算法，提出一种基于CEEMD的联合降噪方法。该降噪方法的主要步骤如下：

步骤1　将输入的原始信号采用CEEMD方法自适应分解为一组按频率从高到低排列的IMF分量。

步骤2　根据排列熵的计算方法计算各个IMF分量的排列熵值，通过排列熵的值来划分出含噪IMF分量有效IMF分量。

步骤3　采用WST方法对含噪IMF分量进行降噪，可得到一组降噪后的含噪IMF分量。

步骤4　将降噪后的含噪IMF分量和有效IMF分量进行重构。

步骤5　对重构信号进行奇异谱分析得到最终的降噪信号。

所提算法流程图如图1所示。

2.2 评价指标

为了定量评价降噪结果，通过信噪比（SNR），均方根误差（RMSE），最大李雅普诺夫指数（MLE）的方法来对比降噪前后的效果。

信噪比（SNR）表示信号平均功率与噪声平均功率的比值，信号的信噪比越大，表示噪声含量越少。信噪比定义如下：

$ SNR=10{\log }_{10}\left(\frac{\parallel x(n){\parallel }^{2}}{\parallel \hat{x}(n)-x(n){\parallel }^{2}}\right) 。$

(14)

式中：$ x(n) $为无噪声信号；$ \hat{x}(n) $为去噪信号。

均方根误差（RMSE）表示去噪信号和无噪声信号之间的差值。RMSE越小，去噪信号和无噪声信号的时域波形越接近，去噪效果越好。均方根误差定义如下：

$ RMSE=\sqrt{\frac{\parallel \hat{x}(n)-x(n){\parallel }^{2}}{N}}。$

(15)

式中：$ N $为时间序列的范围。

李雅普诺夫指数是衡量动态系统复杂性的重要指标，指在相空间中初始时无限接近的2个轨道，随着时间的不断推移按指数收敛或发散的平均变化率，可以定量描述混沌系统在局部范围里系统轨道间的分离程度。在实际应用中，最大李雅普诺夫指数意义重大，其计算公式为：

$ \lambda =\frac{1}{{t}_{m}-{t}_{0}}\sum \limits_{k=1}^{M}\ln \frac{D\left({t}_{k}\right)}{D\left({t}_{k-1}\right)} 。$

(16)

式中：$ \lambda $为最大李雅普诺夫指数；$ D\left({t}_{k}\right) $为在$ {t}_{k} $时刻相空间中相邻轨道间的距离；$ {t}_{m}-{t}_{0} $为总观测时间区间；$ M $为迭代计算过程中的累积步数。该公式通过计算相空间中相邻轨道距离比值的对数平均值来量化系统的动力学特性。具体而言，系统的动力学行为可通过$ \lambda $值进行判定：当$ \lambda \gt 0 $时，相空间中的轨道呈指数分离，系统表现为混沌动力学特征；当$ \lambda =0 $时，轨道既不收敛也不发散，系统处于临界稳定状态；当$ \lambda \lt 0 $时，相空间中的轨道趋于收敛，系统呈现稳定的动力学行为。

3 实　验 3.1 仿真信号降噪

通过Blocks、Bumps、Doppler和Heavysine等常用的仿真信号来研究不同信噪比条件下的降噪效果，从而验证提出降噪算法的有效性，4类仿真信号如图2所示，采样频率为1000 Hz，数据长度为1024点。

以Blocks信号为例，数据长度为1024点，通过叠加高斯白噪声得到信噪比为6 dB的含噪信号，如图3所示。可知，Blocks信号已经完全淹没在噪声中，分别用EMD和CEEMD对含噪信号进行分解，分别得到一组IMF，如图4所示。

分别计算CEEMD分解得到的各模态的排列熵，如表1所示。随着IMF阶数的增加，其排列熵值呈现出递减的趋势。根据排列熵对模态进行划分，前4个IMF的排列熵值大于0.5，可视其为含噪模态，将前四阶含噪模态进行小波阈值降噪，其余排列熵小于0.5的高阶模态为有效模态。

针对前四阶含噪模态采用小波软阈值进行降噪，选取小波基函数为4 dB，分解层数为4。将降噪以后的含噪模态与有效模态重构得到降噪后的Blocks信号。图5为4种降噪方法的结果，其中直接采用EMD，CEEMD分解，剔除噪声模态并重构的方法可分别表示为EMD-PE，CEEMD-PE，采用CEEMD分解将噪声模态进行小波软阈值降噪后并重构的方法表示为CEEMD-PE-WST，所提降噪方法表示为CEEMD-PE-WST-SSA。

为了评判几种降噪方法的效果，采用不同降噪方法分别进行降噪试验，利用信噪比和均方根误差来衡量降噪的效果。表2～表5为4种降噪方法的信噪比和均方根误差均值，可知所提出的降噪方法具有较高的信噪比和较低的均方根误差，降噪效果优于其他方法。

3.2 混沌信号降噪

由于水声信号具有混沌特性，因此利用提出的方法对Lorenz混沌系统进行降噪处理。Lorenz混沌系统Lorenz在1963年发现，从此揭开了混沌研究的序幕。混沌系统是指在一个确定性系统中，存在着貌似随机的不规则运动，其表现为不确定性，不可重复和不可预测。随着混沌理论的发展，许多新的混沌系统被提出并且广泛应用于各个领域。采用典型洛伦兹混沌系统来检验降噪算法的有效性。

洛伦兹系统可以表示为：

$ \left\{\begin{aligned}&x=a\left(y-x\right)，\\& y=cx-y-xz，\\ &z=xy-bz。\\ \end{aligned} \right.$

(17)

在$ a $=10，$ b $=8/3，$ c $=28时呈现混沌态，采用龙格-库塔迭代法计算步长为0.01的$ x $分量，选取2048个点的$ x $分量信号作为洛伦兹信号，在此基础上添加−5 dB的高斯白噪声，采用提出的CEEMD-PE-WST-SSA降噪方法对不同信噪比的含噪洛伦兹信号进行降噪，通过比较降噪前后的波形和吸引子轨迹来验证降噪方法的有效性。

由图6～图8可知，无噪声Lorenz吸引子轨迹光滑且有一定规律，而含噪Lorenz信号的吸引子轨迹无明显规律。采用提出的降噪方法对−5 dB和0 dB含噪Lorenz信号降噪后，其时域波形更接近原信号，降噪后的Lorenz吸引子轨迹也较光滑且有规律。

CEEMD-PE-WST-SSA对不同SNR的Lorenz信号降噪结果如表6所示，信噪比和均方根均值得到明显改善。综上所述，针对含噪混沌信号，提出的降噪方法可以有效去除噪声，实现混沌信号降噪。

3.3 船舰辐射噪声信号降噪

将提出的CEEMD-PE-WST-SSA方法应用于4种实测船舰辐射噪声信号，这四类实测的船舰辐射噪声信号来源于DeepShip水声数据库，对其进行降噪处理，降噪前后的4类舰船信号和吸引子轨迹图及其对比图如图9～图12所示。

可知，4类水声信号存在明显的噪声且吸引子轨迹杂乱无章，采用提出的方法降噪处理后4类水声信号的是波形图更加接近原信号且吸引子轨迹具有一定的规律性。因此，提出的方法可以有效去除噪声，实现水声信号的降噪。为了使结果更具有说服力，采用改进排列熵和最大李雅普诺夫指数2种指标对去噪效果进行评价。熵可以代表时间序列的复杂性，熵减小意味着降低了复杂度；最大李雅普诺夫指数用来衡量动态系统复杂性，最大李雅普诺夫指数越小，表明系统越稳定。如表7所示，降噪后的PE值和MLe值均小于降噪前说明水声信号得到了有效的降噪。

4 结　语

在互补集合经验模态分解的基础上，结合排列熵、小波阈值降噪和奇异谱分析，提出了CEEMD-PE-WST-SSA降噪方法，该方法将CEEMD分解得到的IMF通过PE识别出含噪IMF，将识别出含噪IMF进行小波阈值降噪，将降噪后的含噪IMF和有效IMF重构，重构后的信号经过SSA后实现信号降噪。实验结果表明，不管是Blocks、Bumps、Doppler和Heavysine 4类信号，还是混沌和实测舰船信号，所提出的CEEMD-PE-WST-SSA降噪方法均可以较好地实现降噪。