0 引　言

随着海洋科学技术的发展，对舰船辐射噪声特征提取与识别技术展开深入研究，因水下环境、地势等多种因素，特征提取与识别分类工作尤为困难。舰船辐射噪声产生的机理较为复杂，并且受到海洋环境的影响，具有非线性、非平稳和非高斯特性，传统的处理水声信号方法比较有局限性，已经不能有效的进行处理。

近几年来，经验模态分解(empirical mode decomposition, EMD)的出现，舰船辐射噪声“三非”的问题得到了解决。EMD^[1]可以将非平稳信号分解成多个相对平稳具有不同尺度的振动信号特性的固有模态函数(intrinsic mode function, IMF)，但EMD具有因极值点分布不均匀导致模态混叠等不足问题^[2]。WU等^[3]提出了集合经验模态分解(ensemble empirical mode decomposition, EEMD)方法解决了EMD模态混叠的问题，该方法利用白噪声的特性，多次加入不同幅值白噪声改变信号极值点特性进行分解。李余兴等^[4-5]提出利用EEMD对舰船辐射噪声分解，采用中心频率或能量差的方式进行特征提取。刘千里^[6]提出EEMD的自适应线谱及连续谱提取方法，将EEMD分解的IMF进行线谱的提取，使用余量和剩余的IMF进行连续谱的准确估计。但EEMD属于递归分解，受到端点效应影响，也会存在一定的分解误差。Dragomiretskiy 等^[7]提出了变分模态分解方法，VMD通过迭代搜寻变分模型最优解来确定所知的IMF及其对应的中心频率和带宽。文献[8]提出的VMD和中心频率相结合的特征提取方法，获取信号的特性，有效进行识别分类。在故障诊断中，文献[9]采用VMD与谱峭度法结合对齿轮箱故障信号进行特征提取，该方法成功地检测出齿轮故障。VMD与EMD，EEMD相比，具有良好的鲁棒性和坚实的数学理论基础。

一般来说，舰船辐射噪声可以看作是具有非线性的时间序列，可以从时间序列复杂度的角度进行分析。有许多表示时间序列复杂度特征的算法：如样本熵^[10]、近似熵^[11]、模糊熵^[12]等。其中由Bandt等提出的排列熵(permutation entropy, PE)^[13]是检测时间序列复杂度的新方法，PE算法具有计算简单、抗噪能力强等特点，故在舰船噪声特征提取有所应用^[14]。但排列熵没有考虑幅值问题，因此XIA等^[15]在排列熵的基础上提出加权排列熵（weighted permutation entropy, WPE），然而WPE只考虑单一尺度，忽略其他尺度的有用信息。因此YIN等^[16]提出了MWPE描述不同尺度下时间序列的复杂性，该算法已经广泛应用于故障诊断和管道侧漏领域^[17-19]。但MWPE在粗粒化过程中忽略其他粗粒化的特征信息，影响熵值的准确性。

本文提出一种基于VMD与IMWPE的舰船辐射噪声特征提取方法。IMWPE是对MWPE的一种改进，通过平移均值法弥补MWPE的不足。首先将舰船辐射噪声进行VMD分解成多个IMF，选取最优的IMF作为研究对象，并通过IMWPE算法进行特征提取。最后将所提取的特征输入PSO-SVM进行识别分类验证，实验结果表明了该方法的稳定性和优越性。

1 基础理论 1.1 VMD理论

VMD是继EMD和EEMD方法之后的新一种方法，具体步骤如下：

1）构造变分问题

变分模态分解算法将输入信号分解成不同中心频率和有限带宽 $\{ {u_k}(t)\} (k = 1,2, \cdots ,K)$ 的子序列。通过求解变分问题进行分解，约束变分模型如下式：

$\left\{\begin{array}{l}\underset{\left\{uk\right\},\left\{\omega k\right\}}{\mathrm{min}}\left\{\displaystyle\sum _{k=1}^{K}{\partial }_{t}\left[\left(\delta (t)+\frac{j}{\text{π} t}\right){u}_{k}(t)\right]{{\rm{e}}}^{-j{\omega }_{k}t}{}_{2}^{2}\right\}，\\ \text{s}\text{.t}\text{.}\displaystyle\sum _{k=1}^{K}{u}_{k}(t)=f(t)。\end{array}\right.$

(1)

式中： $\{ {\omega _k}\}$ 为模态函数的中心频率； $\{ {u_k}\}$ 表示各个模式； ${u_k}(t)$ 为第k阶调频信号； $\delta (t)$ 为单位脉冲函数； $f(t)$ 为输入信号。

2）求解变分问题

使用惩罚因子 $\alpha$ 和拉格朗日乘算子 $\theta (t)$ ，再使用二次惩罚因子重建信号，使用贝叶斯推导出惩罚权重，噪声越大惩罚权重越小。 $\theta (t)$ 用于负责条件的严格约束，定义增广拉格朗日 $L$ 公式如下：

$\begin{split} L\left( {\left\{ {{u_k}} \right\},\left\{ {{\omega _k}} \right\},\theta } \right) =& \alpha \sum\limits_{k = 1}^K {\left\| {{\partial _t}\left[ {\left( {\delta \left( t \right) + \frac{j}{{\text{π} t}}} \right){u_k}\left( t \right)} \right]{{\rm{e}}^{ - j{\omega _k}t}}} \right\|_2^2} + \\ & {\mkern 1mu} \left\| {f\left( t \right) - \sum\limits_{k = 1}^K {{u_k}\left( t \right)} } \right\|_2^2 + \\ & < \theta (t),f(t) - \sum\limits_{k = 1}^K {{u_k}(t)} > ，\end{split}$

(2)

采用一种称为交替方向乘子法(ADMM)的迭代优化方法求取式(2)中的“鞍点”，则式(1)最优解求解步骤为：

步骤1　初始化 $\left\{ {\hat u_k^1} \right\},\left\{ {\hat \omega _k^1} \right\},\left\{ {{{\hat \lambda }^1}} \right\},n$ ；

步骤2　循环n=n+1；

步骤3　对所有 $\omega > 0$ ，更新 $\hat u_k^{}$ ， $k \in \{ 1,2,\cdots,K\} ;$ 更新 ${\omega _k}$ ；更新 $\lambda$ ；

$\widehat u_k^{n + 1}(\omega ) \leftarrow \frac{{\hat x(\omega ) - \displaystyle\sum\limits_{i < k} {\hat u_i^{n + 1}(\omega )} - \sum\limits_{i < k} {\hat u_i^n(\omega )} + {{\hat \lambda }^n}(\omega )/2}}{{1 + 2\alpha {{(\omega - \omega _k^n)}^2}}} ，$

(3)

$\omega _k^{n + 1} \leftarrow \frac{{\displaystyle\int_0^\infty {\omega |\hat u_k^{n + 1}(\omega ){|^2}{\rm{d}}\omega } }}{{\displaystyle\int_0^\infty {|\hat u_k^{n + 1}(\omega ){|^2}{\rm{d}}\omega } }},k \in \{ 1,K\} ，$

(4)

${\hat \lambda ^{n + 1}}(\omega ) \leftarrow {\hat \lambda ^n}(\omega ) + \tau \left[\hat x(\omega ) - \sum\limits_k {\hat u_k^{n + 1}(\omega )} \right]。$

(5)

步骤4　重复步骤2和步骤3操作，直到满足式(6)停止迭代，即得到K个IMF分量

$\sum\limits_k {\left\| {\hat u_k^{n + 1} - \hat u_k^n} \right\|} _2^2/\left\| {\hat u_k^n} \right\|_2^2 < \varepsilon 。$

(6)

1.2 加权排列熵

排列熵算法描述了时间序列的复杂性，它只考虑时间序列的顺序结构而忽略序列的幅值信息，WPE解决了对此不足的问题，算法如下：

1）原始时间序列 $X = \{ x(i),i = 1,2,\cdots ,N\}$ 进行相空间重构

$X_i^{(m)} = \{ x_i^{(m)},x_{i + \tau }^{(m)},\cdots ,x_{i + (m - 1)\tau }^{(m)}\} ，$

(7)

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

2）计算每个子序列的权重值

${w_i} = \frac{1}{m}\sum\limits_{k = 1}^m {{{({x_{i + (k - 1)\tau }} - \bar X_k^{(m)})}^2}}，$

(8)

$\bar X_k^{(m)}{\text{ = }}\frac{1}{m}\sum\limits_{k = 1}^m {{x_{i + (k - 1)\tau }}} ,k = 1,2,\cdots m 。$

(9)

3）计算每一种排列出现的概率

任意每个子序列 $X_i^{(m)}$ 的特征信息用加权排列熵和 ${\pi _k}$ 表示，每一种排列出现的概率 ${P_m}({\pi _k})$ 为：

${P_w}({\pi _k}) = \frac{{\displaystyle\sum {\{ {w_i}} |1 \leqslant i \leqslant N - (m - 1)\tau ,i \in {Z^ + },N(X_i^{(m)}){\text{\} }}}}{{\displaystyle\sum {{w_i}} }}。$

(10)

其中， $N(X_i^{(m)})$ 为 $X_i^{(m)}$ 的排列模式 ${\pi _k}$ 。

4）计算加权排列熵的值

$WPE(X,m,\tau ) = - \frac{1}{{\ln m!}}\sum\limits_{k = 1}^K {{P_w}({\pi _k})\ln {P_w}({\pi _k})} 。$

(11)

1.3 多尺度加权排列熵

多尺度加权排列熵是多尺度熵与加权排列熵的结合，是用来描述多尺度时间序列复杂度的方法。MWPE在不同尺度下求取粗粒化时间序列的熵值结果，研究各尺度下熵值的特征，从而反映出舰船辐射噪声的变化情况，可以比单一尺度的加权排列熵获取更多状态量。计算步骤如下:

1）进行粗粒化处理

原始时间序列 $X = \{ x(i),i = 1,2,\cdots ,N\}$ 进行粗粒化处理，得到序列 $y_j^{(s)}$ ，即

$y_j^{(s)} = \frac{1}{s}\sum\limits_{i = (j - 1)s + 1}^{js} {{x_i}} ,j = 1,2,\cdots ,[N/S] 。$

(12)

其中：s为尺度因子；[N/S] 向下取整。

2）计算每个粗粒化序列的排列熵

$MWPE(X,m,\tau ,s) = WPE({y^{(s)}},m,\tau ) 。$

(13)

2 改进多尺度加权排列熵

MWPE通过粗粒化方式解决了WPE尺度单一的问题，但该算法在舰船辐射噪声特征提取上仍存在缺陷：当尺度因子较大时，粗粒化序列会变短，会使熵值产生误差；在粗粒化过程中，忽略了时间序列上的自相关性。自相关度是指信号在同一个时间序列中，某一采样点的瞬时值与另一个采样点瞬时值之间的依赖程度。自相关系数的绝对值接近1就代表相关性很强，绝对值接近0代表相关性很弱。自相关函数公式为：

$R(n) = \frac{{\displaystyle\sum\limits_{i = 1}^{n - h} {[(xi - \mu )(xi + h - \mu )]} }}{{\displaystyle\sum\limits_{t = 1}^n {{{(xi - \mu )}^2}} }}。$

(14)

式中：n为样本容量；h为阶数；μ为该时间序列的均值；x_i代表各采样点。

通过对时间序列的分析，任意一个采样点与其前2个滞后变量呈强自相关度，因此MWPE以s为尺度因子分段时，未考虑分段之间前后的自相关度，从而导致不能够有效地进行特征提取。

针对MWPE存在的不足，提出IMWPE算法解决MWPE单一粗粒化问题。为了考虑相邻时间序列之间的有用信息，本文采用均值移动法计算同一尺度下多个粗粒化时间序列的WPE值。在较大的尺度因子下，IMWPE可以最大程度保留原始时间序列的信息特征。由于在采集舰船辐射噪声信号时会在采集前后给出预留值，因此舰船噪声信号的中间部分最能反映噪声的真实特征，在时间序列上，IMWPE采用累加求均值法将计算资源集中在噪声的中间部分，多次计算时间序列中段来获取重要信息，具体操作步骤如图1所示。

IMWPE整体计算步骤如下：

将信号进行改进粗粒化过程处理，得到s组新序列：

$\left\{ \begin{gathered} z_i^{(s)} = \{ y_{i,1}^{(s)},y_{i,2}^{(s)},\cdots \} ，\\ y_{i,j}^{(s)} = \frac{{\displaystyle\sum\limits_{p = 0}^{s - 1} {{x_{p + i + s(j - 1)}}} }}{s} ，\\ \end{gathered} \right.\quad i = 1,2,\cdots ,s 。$

(15)

计算s个时间序列熵的平均值，如下式：

$IMWPE(x,m,\tau ,s) = \frac{1}{s}\sum\limits_{i = 1}^s {WPE(z_i^{(s)},m,\tau )}。$

(16)

3 特征提取识别分类方法

基于VMD与改进多尺度加权排列熵的特征提取和识别分类流程如图2所示。

1）对获取的三类舰船辐射噪声仿真信号进行EMD分解，确定VMD的分解阶数；

2）VMD分解后得到一定数量的IMF，计算各IMF的能量，确定能量最大的IMF为最优IMF；

3）提取最优IMF的IMWPE；

4）采用PSO-SVM方法对多类舰船噪声VMD分解后最优IMF的IMWPE特征参数进行分类识别并得出结果。

4 实验验证 4.1 舰船噪声模态分解

采用舰船辐射噪声仿真信号模拟真实舰船信号数据，按照线谱与连续谱幅值不同分为A、B、C三类舰船，分别进行EMD分解，确定分解阶数为8，因此设置参数K=8，采样频率为5120 Hz，数据长5120点， $\alpha$ =2000。根据设置的参数，将三类仿真信号再进行VMD分解，图3为3种舰船辐射噪声VMD分解图。

由图3看出，舰船辐射噪声信号分解出8个IMF和一个残差信号，在一定数量的IMF中，至少有一个IMF信号幅值明显高于其他分量，一般将幅值最高的IMF称为信号的主能量，可以反映舰船辐射噪声信号的主要特征。根据文献[4]提出的相关公式，计算每个IMF。设第k阶模态有N个采样点，第i个采样点瞬时频率为 ${f_{ki}}$ ，瞬时振幅为 ${b_{ki}}$ ，瞬时强度为 ${Q_{ki}} = b_{ki}^2$ ，IMF强度为：

${\hat B_k} = \frac{{\displaystyle\sum\limits_{i = 1}^N {{Q_{ki}}} }}{N} 。$

(17)

则能量强度最大的IMF为：

${\hat B_{\max }} = \max ({\hat B_1},{\hat B_2},{\hat B_3}\cdots ,{\hat B_K}) 。$

(18)

根据式(17)和式(18)求出能量最强的IMF，用EIMF表示，则A类与B类的EIMF阶数都为3，C类的EIMF阶数为1。

4.2 IMWPE特征提取

3种舰船辐射噪声信号通过VMD分解后，根据式(17)和式(18)提取能量最大的EIMF作为研究对象，然后提取A类、B类、C类舰船辐射噪声的加权排列熵，3种舰船噪声的样本数各为150个，样本总量为450个，实验结果如图4和表1所示。

由图4和表1可知，加权排列熵特征分布均匀，表明了加权排列熵算法具有较强的稳定性。由于舰船辐射噪声产生的机理不同，导致3种舰船噪声的时间复杂度存在一定的差异。根据平均加权排列熵值得出，B类的噪声信号复杂度最高，其次是A类与C类，从加权排列熵特征分布来看，虽然A类与B类的加权排列熵特征分布较为重叠，但加权排列熵特征分布还是可以很好地区分3种舰船辐射噪声。

加权排列熵是从单尺度方面考虑时间序列复杂度问题，在此引入多尺度加权排列熵从多个尺度了解舰船辐射噪声特性情况。由于多尺度加权排列熵粗粒化过程忽略其他序列信息问题，本文提出改进多尺度加权排列熵算法，提取出的EIMF分别使用MWPE和IMWPE计算后进行对比。一般来说，嵌入维数m通常取值为3～7，m取值太小，重构时间序列中的重要信息量会减少，m值太大，计算时间过长，不能揭示信号变化，本文设置m=5。信号数据长度也会影响熵值的精确度，信号越长，熵值精确度越高，太多的数据点会消耗大量的计算资源，信号过短，则不能准确反映信号的特征变化，文献[19]建议取值为N≥5m!，因此本文取N=5120。延迟时间对熵值计算影响较小，则设 $\tau$ =1，根据研究需要设s=10。图5为提取的EIMF中三类舰船噪声的随机一种熵值对比图。

可以看出，3种舰船在各尺度因子上分布均匀，但提取的MWPE熵值曲线比IMWPE波动较大。这是因为IMWPE提取熵值采用改进的同一尺度多个粗粒化方法，比MWPE算法更为稳定。表2给出2种熵值的标准差对比，说明IMWPE比MWPE更具有稳定性和优越性。

4.3 分类识别检测

为了进一步比较3种舰船辐射噪声特征提取方法，将随机提取的150个样本中70%作为训练样本，30%作为测试样本，并把提取的WPE，MWPE和IMWPE采用PSO-SVM进行分类识别验证。仿真实验中，设置粒子群算法的种群规模为20，最高迭代50次，学习因子C₁和C₂均为1.5^[20]，分类结果如表3所示。

可以看出，WPE特征提取算法总识别率较低，由图4可知，A类与B类的WPE值略有重叠，并且只从单一尺度方面考虑，A类与B类在识别分类中产生误差导致总识别率降低。从图5和表3看出，MWPE与IMWPE都是从多尺度角度分析，C类舰船基本上可以完全识别。因为IMWPE采用平移均值法解决了单一尺度的粗粒化问题，所以总识别率高于MWPE和WPE特征提取算法。

5 结　语

本文从时间复杂度的角度去分析舰船噪声，提出了一种基于VMD与IMWPE相结合的舰船辐射噪声特征提取方法，先将信号通过VMD进行分解，选取能量最大的EIMF作为研究对象，并计算IMWPE作为特征向量，采用PSO-SVM进行分类识别。经实验数据表明：由于不同种舰船辐射噪声产生的机理不同，它们的WPE也大不相同；MWPE忽略了单一粗粒化相邻序列之间的关系，IMWPE算法解决了MWPE不足的问题，具有良好的稳定性和一致性；在识别分类方面中，提取的IMWPE识别率效果最好。