一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法与流程

本发明属于水声工程、海洋工程和声纳技术等领域，涉及一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法，是基于循环调制谱和多重信号分类算法相结合的改进噪声包络信号识别方法用来提取舰船螺旋桨的轴频及各阶的谐波线谱，该方法适用于舰船辐射噪声的目标识别。

背景技术

理论研究及实测表明，实船的螺旋桨噪声是非平稳的宽带连续信号，其频谱信号由连续谱噪声与低频离散噪声两部分组成。连续谱的主要声源是螺旋桨的空化噪声，低频离散噪声又被称为线谱噪声，其频率通常在0-300hz之间，各条线谱相互之间不连续且呈周期性出现，其幅值高于邻近的连续谱信号，而且频率越低幅值越高，一般频率最低幅值最高的线谱就是螺旋桨的轴频，其余那些频率周期性成倍递增的线谱就是螺旋桨轴频的各阶次谐波，各阶次谐波的频率大致上等于轴频与其阶次相乘。螺旋桨噪声中包含着丰富的目标信息，其中螺旋桨的轴频物理特征最为明显，再加上离散线谱噪声易检测，呈周期性出现等特点，螺旋桨轴频始终被视为目标识别的重要特征之一。因此改进螺旋桨线谱的搜索方法，提高轴频检测的精度，对于海上目标的分类和识别有着重要意义。目前常用的螺旋桨线谱搜索算法大都是基于lofar(low-frequencyanalysisandrecording)谱和demon(detectionofenvelopemodulationonnoise)谱的各类改进算法，最新提出的线谱搜索方法是先基于线谱噪声的二阶循环平稳特性进行循环调制谱的计算，然后利用压缩感知方法从循环调制谱中分离线谱，该方法参见“compressivesensingfordetectingshipswithsecond-ordercyclostationarysignatures”，该文2017年9月发表于《ieeejournalofoceanicengineering》第99期，页码为1-13。该方法虽然能够有效的提取线谱，但是循环迭代的收敛速度过慢导致计算时间过长，多次频谱计算会导致频率偏移和频谱泄露，而且利用压缩感知方法会有部分数据在数据压缩的过程中丢失，最终降低轴频检测的精度。

技术实现要素：

要解决的技术问题

为了避免现有技术的不足之处，本发明提出一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法，克服现有的噪声信号识别方法在估计目标轴频时存在的计算量过大、迭代收敛速度慢以及估计精度较差等缺点。

技术方案

一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法，其特征在于步骤如下：

步骤1：对水听器接收到的时域原始信号x(t)做i层的小波包分解，得到最后一层的各段小波包信号，设为第i层小波包分解得到的第j个小波包信号，而原始信号其分解算法为：

其中，参数k为位置坐标，gk和hk为小波包分解滤波器的系数，且gk与hk相互正交，满足gk＝(-1)^kh1-k；最后一层小波包分解得到的每一段小波信号为

步骤2：对最后一层小波包分解得到的每一段小波信号做希尔伯特-黄变换，首先对每一段小波信号进行经验模态分解；然后进行希尔伯特变换求解包络信号，得到连续时域信号为：

式中：为卷积运算符；为小波信号经过m次模态分解后得到信号序列；

步骤3：对所有经过希尔伯特-黄变换的小波信号进行重构得到新的包络信号h[r(t)]，其第i层的第j个小波包的重构计算公式如下：

其中：i＝i,i-1,...2,1，j＝2ⁱ,2ⁱ-1,...2,1，则重构后最终得到的噪声包络信号为h[r(t)]；

步骤4、对重构的包络信号h[r(t)]做加权傅立叶变换计算，求其时间相关谱：

先对连续时间信号h[r(t)]采样得到n点离散信号h[r(n)]，之后采用长度为l2的汉宁窗函数对h[r(n)]进行分段加权计算，这里h[r(n)]被分成l1段，每段的长度含有个数据点，窗函数的重叠率为50％，窗系数为w(n)；因此，对于第l段信号h[r(n)]，进行短时离散傅立叶计算得到时间相关谱定义如下：

w(n)＝0.5-0.5cos(2πn/n)

其中：zl(n,f)是h[r(n)]的时间相关谱，f是时间相关谱信号的频率，fs是整段信号的采样频率；

步骤5、对时间相关谱信号zl(n,f)采用多目标分类算法处理：先计算zl(n,f)的相关矩阵r，并对r做特征值分解，其计算过程如下：

其中：λ1,λ2,...λ为矩阵r中p个大的特征矢量，σ²为噪声的方差；用p个大的特征矢量构成信号子空间us，用n-p个小的特征值对应的特征矢量来构成噪声子空间un；

利用噪声子空间un与信号子空间us构造空间谱：

将pmusic在2～15hz以内的峰值对应的频点记录下来，这些点的频率组成集合f1＝{fq|fq＝fmusic,2＜fq＜15,q＝1,2,3…q}，fmusic为空间谱pmusic的频率，空间谱公式构造为：

其中：a为搜索矢量，a^h为a的共轭转置，为un的共轭转置矩阵；

其中的pmusic为2～15hz频段内的极大值即为疑似轴频

将f1中的每个元素fq视为疑似轴频，分别以fq为中心频率，以0.01＜α＜0.1大小的步长生成q个集合，每个集合为这里g为中心频率加减步长的个数；

步骤6：对步骤4中计算得到的zl(n,f)的平方再一次做离散傅立叶计算，将时间变量n转换为双循环频率变量f双，推导得到二阶循环相关信号频谱为：

其中：f双为集合中的每个疑似轴频；f既是每段时间相关谱zl(n,f)的频率也是循环调制谱中载频f载的中心频率，即上述公式计算得到的循环调制谱被最终表示为p(fq,g,f载)，是关于载频f载和双循环频率变量的频谱；

步骤7：对于只包含单个螺旋桨噪声的原始信号，定义每个集合的频率统计量为pg＝argmax(p(fq,g,f载))，fq,g是步骤6得到的二阶循环频率分量，定义轴频是max(pg)对应的频率；

步骤8：以步骤6检测得到的循环调制谱p(fq,g,f载)在1～100hz频段内的载频f载分别与比较搜索，搜索公式如下：

其中：r是对取整后四舍五入的值；是选定的搜索各阶次谐波的偏移精度阈值，取值为小于等于频谱分辨率；

步骤9：将满足公式所得到每个有效频点f载按从小到大排序，得到频率成倍数递增的一系列线谱，为搜索得到的目标螺旋桨的轴频与各阶次的谐波，线谱序列中的第一条谱线就是轴频剩余的线谱f载就是各阶次的谐波。

所述步骤5的p个大的特征矢量的p≥10。

所述步骤5的中心频率加减步长的个数

有益效果

本发明提出的一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法，先对原始水听器信号求包络，也就是对该信号进行多尺度小波包分解，然后对每个小波包信号做希尔伯特-黄变换得到小波包络信号，之后将所有的小波包信号进行小波重构，得到改进的噪声包络信号。然后对包络信号进行分段加权计算其时间相关谱，接下来结合多目标分类算法处理时间相关谱信号，取其在低频段峰值突出的线谱作为疑似轴频，根据需求设定小于频率分辨率的步长，以疑似轴频为中心组成疑似轴频集合，再对各疑似轴频的时间相关谱信号做循环相关计算，得到循环调制谱的频率统计量，最终根据该频率统计量来搜索目标的轴频。本发明方法校正了传统的循环调制谱方法在频谱计算过程中产生的轴频偏移误差，提高了轴频检测的精度，极大地减小了循环迭代的计算量，而且在工程应用方面便于实现。

本发明有益效果是：步骤一、二、三改进了传统噪声包络信号求解方法，步骤一中该方法利用权利所述的尺度函数与小波函数构造小波包函数集合对时域滤波信号做i层的小波包分解，最后得到共2ⁱ个小波包分解信号，这些分解信号将频带进行了多层次划分，能够提高信号的时频分辨率；步骤二对这2ⁱ个小波包信号进行希尔伯特-黄变换来求解信号的包络，通过多次经验模态分解去除信号中包含的连续谱噪声分量；步骤三中对所有的小波包络信号进行小波包重构，得到改进后的时域噪声包络信号。接下来步骤四、五、六对时域噪声的包络信号进行改进循环调制谱分析，步骤五中该方法利用多目标分类算法处理通过步骤四的计算得到的时间相关谱，将信号在2-15hz频段的一系列极大值视作疑似轴频，之后以各个疑似轴频为中心频率以适宜步长衍生出多个轴频集合；步骤六中仅仅需要计算轴频集合中的各频率分量的循环调制谱，不需要计算所有频率分量，不仅达到了抑制非二次相位耦合的谐波分量，增强周期性的线谱的目的，还极大减轻了传统循环调制谱的计算量。步骤七与八中统计了每个轴频集合的频率统计量，分别针对单目标信号与多目标信号进行轴频搜索，从中选出符合要求的目标轴频与各阶次谐波。该方法避免了人工挑选疑似轴频，衍生的轴频集合不仅修正了传统的循环调制谱两次频谱计算时带来的误差，还减少了循环调制谱迭代计算所花费的时间。该方法通过实验验证适用于工程上检测目标轴频，即使在低信噪比情况下也具有较高的精度。

附图说明

图1改进噪声包络信号谱的轴频检测流程示意图

图2实测噪声信号的时域波形图及频域波形图，图2(a)为噪声信号的时域波形图，图2(b)为噪声信号的频谱波形图。

图3传统噪声包络信号的时频图(即lofar图)。

图4改进噪声包络信号的时频图(即lofar图)。

图5多目标分类识别算法得到的疑似轴频的频率空间谱图

图6改进后的噪声信号的循环调制谱图

图7目标螺旋桨的轴频及各阶谐波的频谱图

图8根据轴频计算得到的船速和实际船速的对比

具体实施方式

现结合实施例、附图对本发明作进一步描述：

一种基于改进的循环调制谱方法的螺旋桨轴频搜索方法，其特征在于：首先是通过小波包变换和希尔伯特-黄变换相结合的方式，对经典的包络调制谱中求解包络的步骤加以改进。之后在循环调制谱的第一次频谱计算之后添加多目标分类算法这一步骤，来降低第二次调制谱的计算量，提高轴频检测的精度。具体来说是先对原始信号进行多尺度小波包分解，再对每个小波包信号做希尔伯特-黄变换求包络，然后对所有的小波包络信号重构得到改进的噪声包络信号；之后，改进求解包络信号的循环调制谱，在循环调制谱的两次频谱计算之间添加了一个步骤，即先计算包络信号的时间相关谱，再利用多目标分类算法计算该时间相关谱信号在低频段的局部最大值点，接着以局部最大值点为中心和小于频率分辨率的步长组成疑似轴频集合，之后的循环调制谱仅需要计算集合中的每一个疑似轴频而不是全频段所有频率分量，计算结果作为频率统计量，最终根据该频率统计量来搜索目标的轴频。其改进过程分为以下步骤：

步骤一：先对水听器接收到的时域原始信号x(t)做i层的小波包分解，得到最后一层的各段小波包信号，设为第i层小波包分解得到的第j个小波包信号，而原始信号其分解算法为：

其中，参数k为位置坐标，gk和hk为小波包分解滤波器的系数，且gk与hk相互正交，满足gk＝(-1)^kh1-k。最后一层小波包分解得到的每一段小波信号为

假定小波包变换的尺度函数是小波函数是ψ(x)，小波包函数的定义如下：

这里，定义函数集合{ψn(t)}为正交尺度函数所确定的小波包。

步骤二：对最后一层小波包分解得到的每一段小波信号做希尔伯特-黄变换。希尔伯特-黄变换操作包含两部分，一是经验模态分解；二是希尔伯特变换求解包络信号。经验模态分解是对输入信号求取极大值点和极小值点，之后对极大值点和极小值点采用三次样条函数插值来构造的上下包络，再计算上下包络的均值函数；减去包络平均后的均值函数得到本征模态分量sv(t)，若sv(t)还存在负的局部极大值和正的局部极小值，则需要继续进行经验模态分解。定义信号减去均值函数得到的本征模态分量sv(t)，然后以sv(t)作为下一次经验模态分解的输入信号不断的进行模态分解，直到满足条件假设条件才会结束，即：

其中，sv(t)为每一次模态分解得到的本征模态分量，为m次模态分解后减去所有本征模态分量最终得到的信号序列。当sv(t)同时满足条件：(1)sv(t)的均值函数的平均值趋近于0；(2)极值点个数(包括极大值点个数+极小值点个数)和其同y轴的交点个数之差不能大于1(小于等于1)时，模态分解结束。

小波信号经过m次模态分解，最终得到信号序列将作为输入信号进行希尔伯特变换得到连续时间信号的希尔伯特变换等价于，将信号与冲激响应为的线性系统做卷积运算，对做希尔伯特，其本质是令在频域的频率分量幅度不变，产生90度相移，对于连续时域信号其公式计算如下：

式中为卷积运算符。

步骤三：对所有经过希尔伯特-黄变换的小波信号进行重构得到新的包络信号h[r(t)]，其第i层的第j个小波包的重构计算公式如下：

其中i＝i,i-1,...2,1，j＝2ⁱ,2ⁱ-1,...2,1，gk和hk分别为小波包重构的滤波器系数，与公式(1)中相同，则重构后最终得到的包络信号为h[r(t)]。

步骤四：对重构的包络信号h[r(t)]做加权傅立叶变换计算，求其时间相关谱。先对连续时间信号h[r(t)]采样得到n点离散信号h[r(n)]，之后采用长度为l2的汉宁窗函数对h[r(n)]进行分段加权计算，这里h[r(n)]被分成l1段，每段的长度含有个数据点，窗函数的重叠率为50％，窗系数为w(n)。因此，对于第l段信号h[r(n)]，进行短时离散傅立叶计算得到时间相关谱定义如下：

w(n)＝0.5-0.5cos(2πn/n)(8)

其中，zl(n,f)是h[r(n)]的时间相关谱，f是时间相关谱信号的频率，fs是整段信号的采样频率。

步骤五：对时间相关谱信号zl(n,f)采用多目标分类算法处理。即先计算zl(n,f)的相关矩阵r，并对r做特征值分解，其计算过程如下：

其中，λ1,λ2,...λ为矩阵r中p个大的特征矢量，σ²为噪声的方差。用p个大的特征矢量构成信号子空间us，用n-p个小的特征值对应的特征矢量来构成噪声子空间un。根据以往对螺旋桨噪声数据处理的经验：理想状况下轴频是幅值最大且频率最小的线谱，而随着阶次的递增各阶次谐波的幅值会逐渐降低，通常谐波阶次大于10的线谱往往会淹没在连续谱中，没有特殊需要一般也不会去检测阶次大于10的谐波。因此，采用多目标分类算法时对特征子空间的选取一般要满足p≥10。

利用噪声子空间un与信号子空间us构造空间谱pmusic，并求取谱函数的极大值点。由于2hz以下频点易产生强干扰，将pmusic在2～15hz以内的峰值对应的频点记录下来，这些点的频率组成集合f1＝{fq|fq＝fmusic,2＜fq＜15,q＝1,2,3…q}，fmusic为空间谱pmusic的频率。空间谱公式构造为：

其中a为搜索矢量，a^h为a的转置，un为噪声子空间。pmusic在2-15hz频段内的极大值就是疑似轴频。

将f1中的每个元素fq视为疑似轴频，分别以fq为中心频率，以0.01＜α＜0.1大小的步长生成q个集合，每个集合为这里g为中心频率加减步长的个数，取值要满足

步骤六：对步骤四中计算得到的zl(n,f)的平方再一次做离散傅立叶计算，将时间变量n转换为双循环频率变量f双，最终推导得到二阶循环相关信号频谱如下：

其中，f双的取值仅选择集合中的每个疑似轴频，fs是zl(n,f)的采样频率，f既是每段时间相关谱zl(n,f)的频率也是循环调制谱中载频f载的中心频率，即公式(11)计算得到的循环调制谱被最终表示为p(fq,g,f载)，是关于载频f载和双循环频率变量的频谱，该频谱不需要计算每段时间相关谱zl(n,f)在全频域内的所有频率分量，仅需计算疑似轴频，极大减少了传统循环调制谱的计算量。

步骤七：对于只包含单个螺旋桨噪声的原始信号，定义每个集合的频率统计量为pg＝argmax(p(fq,g,f载))，fq,g是步骤六中计算得到的二阶循环频率分量，定义轴频是max(pg)对应的频率。

步骤八：理想状态下，轴频对应的各阶次谐波的频率应该是其本身的倍数，即第i阶谐波的频率值为但实际上由于螺旋桨叶片的不均匀导致各阶次谐波的频率并不全是轴频的倍数，存在一定的偏移差。因此当需要识别各阶次谐波的频率时应选择合适的阈值作为偏移精度范围。

此时，用步骤六检测得到的循环调制谱p(fq,g,f载)在1-100hz频段内的载频f载分别与比较搜索，搜索公式如下：

其中：r是对取整后四舍五入的值；是选定的搜索各阶次谐波的偏移精度阈值，取值一般要小于等于频谱分辨率。

步骤九：记录下满足公式(12)搜索得到每个有效频点f载，将其按从小到大排序，可以得到频率成倍数递增的一系列线谱，这就是搜索得到的目标螺旋桨的轴频与各阶次的谐波，该线谱序列中的第一条谱线就是轴频剩余的线谱f载就是各阶次的谐波。

图1是本发明方法提出的基于循环调制谱和多目标分类算法相结合的改进噪声包络信号识别方法的流程图。如图所示，本发明方法的具体流程是对水听器接收到的原始信号先做小波包分解，每组的小波包信号都经过希尔伯特-黄变换后再重构成一个包络信号，然后对包络信号分段加权计算得到时间相关谱，再利用多目标分类算法处理该时间相关谱，筛选出该时间相关频谱在0-15hz频段的极大值，再以这些极大值的频点作为频率中心以适宜步长扩展，得到疑似轴频集合，最后计算疑似轴频的循环调制谱，得到频率统计量后搜索目标的轴频。

图2是本发明方法在海上试验中实际采集到的舰船噪声，图2(a)为实测舰船噪声的时域波形图，图2(b)为实测舰船噪声的频谱图。如图2(b)所示，舰船噪声的线谱淹没在连续谱中，无法有效识别。

图3是本发明方法采用短时傅立叶变换处理传统的噪声包络信号，得到的低频搜索和测距谱(简称lofar)图，图3中的频谱信号是先对实测舰船信号带通滤波后做希尔伯特变换得到包络信号，之后再做短时傅立叶变换得到的。带通滤波采用有限冲激响应数字滤波器，参数设置通带一般在1-300hz左右，阻带衰减不小于6db，包络信号是采用希尔伯特变换计算带通滤波器的输出信号得到的。根据图3，尽管能够确定信号所包含的离散线谱的大体位置，但无法做到准确识别。

图4是本发明方法采用短时傅立叶变换处理改进过的噪声包络信号，得到的低频段的低频搜索和测距谱(简称lofar)图。图中的改进噪声包络信号是利用步骤一到三处理实测舰船噪声得到的，该实测舰船噪声与图2所示噪声相同。对比图3，经过小波包变换和希尔伯特-黄变换计算后，图4的改进噪声包络信号可以从频谱中识别出一系列短时间内频率保持不变的离散线谱，而这些离散的线谱中就包含着螺旋桨的轴频与各阶次谐波。

图5是本发明方法采用多目标分类算法处理包络信号的时间相关谱，得到的0-15hz频段内的频谱。如图5所示，0-15hz频段内的极大值频点3.052hz、6.104hz、7.935hz、10.99hz、14.04hz，设定信号子空间维数为p＝12，分别以3.052hz、6.104hz、7.935hz、10.99hz、14.04hz为中心频率，以0.01大小的步长扩展生成了5个疑似轴频集合。

图6是本发明方法采用循环调制谱处理疑似轴频集合中每个频率分量，得到的低频搜索和测距谱(简称lofar)图。如图所示，图中不随载频f变化的离散谱线就是检测得到的频率分量具有二阶循环平稳的特性，且保留了轴频与各阶次谐波的周期性。对比图3和图4，图6中对离散线谱的检测精度有了很大提高，且计算量远小于传统的循环调制谱。

图7是本发明方法利用步骤七、八计算频率统计量，得到的目标螺旋桨轴频与各阶次谐波的频谱。如图7所示，图中标注出的10.38hz是利用本发明方法最终检测到的实测舰船辐射噪声信号中的轴频；其余幅值随频率递增而逐渐降低的线谱为搜索到的各阶次谐波，设定搜索阈值为各阶次谐波根据公式(12)从载频f载中搜索得到。图7中检测到的目标轴频10.38hz相比图5采用多目标分类算法中得到的极大值频率10.99hz有0.61hz的偏移校正。

图8是采用全球定位系统(gps)测定的实验船航速与采用本发明方法检测所得轴频计算出的实验船航速对比图，实际上是利用实测数据来验证本发明提出的轴频检测方法的有效性。具体实施过程是先通过全球定位系统(gps)得到不同航速下实验船的辐射噪声数据，记录下当前实验船的航速，再利用本方法计算实验船的辐射噪声信号得到轴频估计值，然后通过公式(13)计算得到当前实验船的瞬时估计航速，最后与全球定位系统记录下的航速对比。

式中：f轴为轴频，v为航速，d为实验船的螺旋桨直径，j为螺旋桨的进速系数，每条船的d、j都是固定不变的。

图中离散的点“*”代表本发明方法从108个辐射噪声数据中搜索得到轴频后，依据公式(13)计算出对应时刻的实验船航速，“o”代表全球定位系统测定的整条航线上的实验船航速变化。两者之间吻合度非常高，证明了本发明方法对于实测舰船噪声数据的有效性。

本发明在典型实施例中取得了明显的实施效果，与现有技术相比其优越性在于无需频率误差校正，轴频搜索精确且计算量小，便于应用。