基于EMD-1（1/2）维谱的舰船辐射噪声调制特征提取

许劲峰 郑 威

（江苏科技大学 镇江 212003）

1 引言

舰船辐射噪声主要包括机械噪声、螺旋桨噪声和水动力噪声，其频谱往往表现为宽带连续谱和线谱的叠加，不同类型的舰船螺旋桨在不同条件下航行，使得它们在螺旋桨噪声方面形成了差异，可用于舰船检测和分类识别，得到了广泛的研究。在舰船辐射噪声中，螺旋桨节拍对宽带辐射噪声存在着明显的振幅调制。这种幅度调制信号即包络信号，携带有许多重要的信息，其调制频率和调制深度与螺旋桨转动的轴频、叶片频及航速有关。对调制谱的准确提取，可有效地运用到推算舰船的速度和类型，具有重要意义［1～4］。

经过专家学者们几十年的研究，对舰船辐射噪声研究取得了一系列的进展，获得了许多显著成果。通过分析舰船噪声产生的原理，对各部分机械推进器和其他装置进行了深入研究，对各噪声的纵向分布和噪声的非高斯性做了大量工作。陶笃纯等分析了舰船辐射噪声中调制包络的物理特性，并通过研究指出了舰船辐射噪声中几种常见的调制类型。将舰船辐射噪声中包络调制作为具有相同形状和重复周期、随机幅度的脉冲序列来处理，并在这种理论下研究计算了各种类型调制包络的功率谱密度函数［5］。史广智，胡均川研究了大量的舰船辐射噪声样本数据，通过舰船辐射噪声调制谱谐波族的结构特点建立了螺旋桨空化噪声信号模型，并根据空化噪声模型推出了谐波族结构特性表达式。然而，通过该模型建立的谐波族结构特性容易受到舰船噪声的信噪比、平稳性等因素的影响［6］。文献［7～8］在对线谱进行研究时，讨论了基于功率谱的方法提取线谱的主要内容以及优缺点。文献［9］利用离散傅里叶变换和Welch算法对舰船辐射噪声信号进行了线谱提取，主要是用Welch算法对信号比较密集频段进行提取，但这样提取的线谱信息不能代表全局的舰船线谱信息，对后续研究的结果影响较大。文献［10］利用辐射噪声中线谱的低起伏特性，使用了短时傅里叶分析法（STFT）对辐射噪声特征线谱进行了提取。文献［11～12］通过对舰船噪声的摸底试验，较为详尽地研究了它的组成形式，并且根据线谱的相关文献资料及其参数数据的研究，建立了舰船辐射噪声线谱的提取模型。紧接着对线谱的平稳性以及特有性质进行了论证，在提取舰船线谱特征之后，通过神经网络等方法对舰船特征线谱的类型进行划分。

文献［13～15］通过在线谱图的分析中，提出了一种较为优越的优化算法，其主要是利用代价函数对线谱图进行分析，并通过这种函数直接对线谱进行提取。在研究过程中，为了避免检测效率低下、虚警概率较高的情况，文中提出了一种新的智能的方式，即利用对双门限检测流程的模仿，把线谱判别、检测以及追踪相结合。这样就能够在比较低的信噪比中完成对舰船特征线谱的自动检测。

文献［16～17］根据DEMON谱法研究了舰船噪声的频谱以及频率特征，因为调制线谱的方位和轴频与叶片频率相对应，所以对解调出的调制线谱进行谱分析，得到轴频和叶片频调制谱，进一步分析舰船线谱中的叶片数量以及螺旋桨的速度，为被动声纳目标检测和分类提供了有力的依据。文献［18］通过对噪声特性的研究，结合经典功率谱和DEMON谱分析法对辐射噪声线谱与连续谱分离、去除虚警及归并线谱，有效地提取了特征线谱。文献［19］综合运用改进的最大公约数算法和余数门限算法提取DEMON谱中的轴频和叶频，解决了传统最大公约数算法提取轴频叶频误差较大的问题。

传统的DEMON谱分析在提取调制谱时，需要人为将信号分成几个频段，然后进行宽带或窄带解调，最后滤波进行傅里叶计算，结果稳定性差、误差较大。EMD方法依据数据自身的时间尺度特征来进行信号分解，无须预先设定任何基函数，这比小波分解，DEMON谱分析的宽带、窄带解调更具灵活性和适用性。1（1/2）维谱具有抑制高斯噪声，分析谐波信号时，信号的基频分量可得到加强，非谐波项可被清除等性质。本文就是结合EMD分解和1（1/2）维谱分析的优点，提取舰船噪声中的有效信息。

2 原理分析

2.1 EMD原理

EMD方法是美籍华人N.E.Huang等于1998年提出的，目的是通过对非线性非平稳信号的分解获得一系列表征信号特征时间尺度的IMF，使得各个IMF是窄带信号，可以进行HS分析。IMF要满足两个条件：1）整个数据集的极大极小值数目与过零点数目相等或最多相差一个；2）数据集的任意点上，由极大值确定的包络与由极小值确定的包络的均值始终为零。这两个条件实际上使得分解得到的IMF是窄带信号。而且EMD分解基于下面的假设：1）信号至少有两个极值，一个极大值和一个极小值；2）信号特征时间尺度是由极值间的时间间隔确定的；3）如果数据中缺乏极值点，但存在缺陷点，可以通过微分、分解、再积分的方法获得IMF［20］。

经验模态分解也称筛选过程［21］（The sifting processing），分解过程具体如下：

1）找出原始信号x(t)上的所有极大值和极小值点，然后分别对所有的极大值点和极小值点用三次样条函数进行连接，拟合出原数据序列上的包络线，上包络xmax(t)和下包络线xmin(t)。此时，这两条包络线就包含了原信号的所有数据，对上下包络线求均值，得到一条均值线m1(t)。

2）再用原信号x(t)减去上下包络均值m1(t)得到一个新的数据序列：

3）对于不同的信号x(t)，可能存在着负的局部极大值和正的极小值点，这时它并不满足IMF的条件，需要把h1(t)作为原信号继续进行“筛选”，重复上述步骤，得到：

上式中，m11(t)是h1(t)的上下包络线的均值，如果h11(t)不是固有模态函数，则继续进行筛选，重复上述处理过程k次，直到得到的h1k(t)符合IMF的要求：

h1k(t)是不是一个IMF，需要由筛选过程终止准则来确定，终止准则是利用两个连续处理结果之间的标准差Sd的值来作为判断是否终止的依据，Sd的定义如下：

上式中，h1(k-1)(t)和h1k(t)是在筛选IMF过程中，两个连续处理结果的时间序列。Huang建议Sd的典型取值范围在0.2～0.3之间。同时，Terradas指出将Sd的取值范围定为0.05～0.3之间会产生更好的结果，通过对太阳日冕数据的研究，发现其第一个IMF（高频噪声项）和趋势项基本保持不变，具有一定的鲁棒性［24］。

4）当h1k(t)满足Sd的要求时，则h1k(t)为第一阶固有模态函数，此时的h1k(t)就是最高频率分量c1(t)，即

从x(t)中分离出c1(t)分量，得到一个差值信号r1(t)，即

将r1(t)看作是一个新的信号重复上述经验模态分解的过程，经过多次运算得到全部剩余分量的ri(t)：

当满足条件：小于预定的误差或成为一单函数，即不可能再从中提取IMF时，就终止分解过程。

信号经EMD分解后，原信号x(t)就可以由n阶分量和残余分量函数rn(t)构成，即

图1给出了EMD算法的流程图。

图1 EMD算法流程图

2.2 1（1/2）维谱的性质与算法

2.2.1 （1/2）维谱的定义

对于随机变量x(t)，它的三阶累积量C3x(τ1,τ2)的对角切片为C3x(τ , τ)(τ1= τ2=τ)，定义该对角切片的傅立叶变换为随机变量 x(t)的1（1/2）维谱C(w)。

上式也可以通过频域计算，如下：

其中，X(w)为x(t)的傅里叶变换，X*(w)为X(w)的复共轭。

离散信号的计算方法如下：

设原始信号为{x1,x2,…,xN=KM}共K段，每段的长度是M，计算1（1/2）维谱的算法如下：

1）对每段数据去直流；

2）分别计算每段数据的三阶累积量；

式中i=1,2,…,K，s1=max(0 ,τ)，s2=min(M-1,M-1-τ)

3）对每段数据的c(i)(τ)取平均，得

这种变化的根本原因是在恒温、恒定表面面积条件下，对表面层施加压力，液体体积会发生变化，导致液体表面层体积变化(见图6)，压力升高会使液体的表面张力数值下降。随着体系压力的升高，烃与CO2的接触面积逐渐减小，界面张力逐渐降低。烃组分界面张力降低的幅度差，表现为CO2与原油的不同烃组分动态混相过程的难易程度的差别。烃组分最外层不断被CO2高密度流体所饱和并不断溶解到CO2中，轻烃不断地补充到混相界面层，混相界面层为促进混相的过渡带，烃组分碳数越小，混相层中的组分溶解速度越快，最终越容易达到完全混相。

4）对 c˄(τ)做一维傅里叶变换，得到信号的1（1/2）维谱。

在对非线性相位耦合使用1（1/2）维谱分析时，需要适当的给离散信号序列进行分段，M和K的值需要折中选取，可以获得比较小的方差。

2.2.2 1（1/2）维谱的性质与仿真

1）1（1/2）维谱具有以下性质：

性质1 设x(t)为零均值、基频为w0的n次实谐波信号，在幅值a相等，相位为零的情况下，当|wm|＜ | wl|时，有

其 中 wm=mw0,m=±1,±2,…±n,wl=lw0,l=±1,±2,…,±n 。该性质说明，当选用1（1/2）维谱分析谐波信号时，信号的基频分量可以得到加强，这有利于抽取信号中较弱的基频分量。要说明的是该性质中的零相位假设，显然，在实际中各谐波分量的相位不可能为零，但这不是问题，可以在计算1（1/2）维谱前对信号序列进行自相关处理，正弦信号的自相关函数仍为同频率的正弦信号，但其相位是零，幅值为原幅值平方的二分之一。

性质2 设n(t)是零均值的随机噪声，任何两个不同时刻都互不相关，且概率密度函数 f(n)为对称分布，则有

该性质表明当信号中混有对称分布的随机噪声时，理论上也可被1（1/2）维谱完全抑制掉。

性质3 设n(t)为零均值的高斯噪声，则有

该性质表明1（1/2）维谱同样具有抑制高斯随机噪声的特点。

性质4 设x(t)是谐波信号，wm、wp、wq为其中三个谐波分量，若wm≠wp+wq，则有

从该性质可知，当信号中含有非相位耦合的谐波项时，通过1（1/2）维谱的处理，这些谐波项可以被清除掉。

2）1（1/2）维谱特性仿真分析

（1）仿真研究1（1/2）维谱的去噪能力和加强基频分量特性。仿真信号为频率分别为50Hz、100Hz、150Hz、200Hz、250Hz的正弦信号和高斯白噪声的叠加，信噪比为-10dB。仿真结果如图2、图3、图4所示，图3为信号的功率谱，图4为1（1/2）维谱，可以看出采用1（1/2）维谱分析后，信噪比提高了很多，信号的基频分量也得到了加强。

图2 原始信号和叠加噪声后的信号

图3 信号功率谱图

图4 信号的1（1/2）维谱图

（2）1（1/2）维谱剔除信号中非耦合谐波分量特性。仿真信号为正弦信号叠加高斯白噪声，信噪比-10dB，其中正弦信号为

它的频率分别为 30Hz、50Hz、100Hz、150Hz、200Hz。图5为它的功率谱，图6是1（1/2）维谱。从图中的结果可以看出，采用功率谱分析时，信号中所有频率成分均清晰可见，而采用1（1/2）维谱分析后，只有那些为整数倍的频率成分清晰可见，而非耦合谐波部分被剔除。

图5 信号功率谱图

图6 1（1/2）维谱剔除信号中非耦合谐波分量

1（1/2）维谱的这一性质有利于分析提取舰船辐射噪声信号，本文采用对提取包络信号求1（1/2）维谱来求其基频和谐波，从而判断分析舰船航速，螺旋桨等相关信息。

3 特征线谱提取

传统的DEMON谱分析法，首先对信号进行宽带解调，需要事先确定好频带范围，然后进行解调，这样得出的结果误差比较大，EMD分解是依据数据自身的时间尺度特征来进行信号分解，不需要事先设定任何基函数，这比小波分解，DEMON谱分析的宽带、窄带解调更具灵活性和适用性，而1（1/2）维谱分析对比功率谱分析具有抑制高斯噪声；分析谐波信号时，信号的基频分量可得到加强，非谐波项可被清除等优势，因此结合EMD分解和1（1/2）维谱分析对DEMON谱分析法进行改进，分解过程如下图7所示。

图7 基于EMD-1（1/2）改进的DEMON谱分析法流程图

4 数据处理与分析

本次实测数据为某大型货船顺风高速航行实测数据，为了去除不确定随机误差，选取了中间较为稳定测量信号，其长度为10s，采样频率为10000Hz，然后对噪声信号进行EMD分解。观察发现，信号分解出的IMF包含了原始信号的绝大多数信息，而最后的残余分量是信号的趋势项没有保留，分解出了六阶IMF分量作进一步特征提取对象，并对此阶分量作了谱分析。前三阶IMF分量可以看出明显的调制现象，且前三节IMF分量包含了信号的主要组成部分，所以对这三阶IMF进行平方低通解调，然后进行谱分析，观察看到了得到了辐射噪声信号线谱信息。

首先，由于实测数据量信息庞大，对数据进行截取，得到中间比较稳定的部分，然后对数据进行谱分析，观察辐射噪声信号的频率大致分布情况，如图8所示，信号能量主要集中在高频段，包含了螺旋桨对宽带噪声的调制。

图8 舰船辐射噪声原始信号与功率谱图

对截取的数据进行EMD分解，得到了前六阶IMF分量，如图9所示。然后对这6个分量进行谱分析，如图10所示，观察发现前三阶分量中，有明显的调制包络现象，然后进一步谱分析，观察线谱信息。

对EMD分解的前三阶IMF分量进行谱分析，图11是功率谱分析，图12是1（1/2）维谱分析，对比发现都能看到线谱信息，但是1（1/2）维谱分析抑制噪声效果明显，得到的线谱比较清晰。图13是1（1/2）维谱分析过后的局部放大图，可以看到辐射噪声的轴频基频为32Hz，且二次谐波分量和三次谐波分量信息明显。

图9 EMD分解图

图10 EMD分解功率谱图

图11 前三阶IMF分量平方检波功率谱图

5 结语

针对传统的DEMON谱分析舰船辐射噪声调制特征线谱的缺点，即在提取调制谱时，效果不太理想，存在对噪声抑制能力差、提取后的谐波特征不明显等缺点。通过研究EMD分解和1（1/2）维谱分析的特性，对传统DEMON谱分析作进一步改进，得到了基于EMD-1（1/2）维谱的舰船辐射特征提取方法。通过仿真分析和实测数据的处理，可以看出本方法在舰船辐射噪声特征提取方面比传统的DEMON谱分析更具优势。