CN108921014B - 一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法 - Google Patents

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

Info

Publication number
CN108921014B
CN108921014B CN201810487439.XA CN201810487439A CN108921014B CN 108921014 B CN108921014 B CN 108921014B CN 201810487439 A CN201810487439 A CN 201810487439A CN 108921014 B CN108921014 B CN 108921014B
Authority
CN
China
Prior art keywords
frequency
spectrum
signal
carrier
noise
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810487439.XA
Other languages
English (en)
Other versions
CN108921014A (zh
Inventor
杨坤德
周星月
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Northwestern Polytechnical University
Original Assignee
Northwestern Polytechnical University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Northwestern Polytechnical University filed Critical Northwestern Polytechnical University
Priority to CN201810487439.XA priority Critical patent/CN108921014B/zh
Publication of CN108921014A publication Critical patent/CN108921014A/zh
Application granted granted Critical
Publication of CN108921014B publication Critical patent/CN108921014B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/02—Preprocessing
    • G06F2218/04—Denoising
    • G06F2218/06—Denoising by applying a scale-space analysis, e.g. using wavelet analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/12—Classification; Matching

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Signal Processing (AREA)
  • Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

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

Description

一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法
技术领域
本发明属于水声工程、海洋工程和声纳技术等领域，涉及一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法，是基于循环调制谱和多重信号分类算法相结合的改进噪声包络信号识别方法用来提取舰船螺旋桨的轴频及各阶的谐波线谱，该方法适用于舰船辐射噪声的目标识别。
背景技术
理论研究及实测表明，实船的螺旋桨噪声是非平稳的宽带连续信号，其频谱信号由连续谱噪声与低频离散噪声两部分组成。连续谱的主要声源是螺旋桨的空化噪声，低频离散噪声又被称为线谱噪声，其频率通常在0-300Hz之间，各条线谱相互之间不连续且呈周期性出现，其幅值高于邻近的连续谱信号，而且频率越低幅值越高，一般频率最低幅值最高的线谱就是螺旋桨的轴频，其余那些频率周期性成倍递增的线谱就是螺旋桨轴频的各阶次谐波，各阶次谐波的频率大致上等于轴频与其阶次相乘。螺旋桨噪声中包含着丰富的目标信息，其中螺旋桨的轴频物理特征最为明显，再加上离散线谱噪声易检测，呈周期性出现等特点，螺旋桨轴频始终被视为目标识别的重要特征之一。因此改进螺旋桨线谱的搜索方法，提高轴频检测的精度，对于海上目标的分类和识别有着重要意义。目前常用的螺旋桨线谱搜索算法大都是基于LOFAR(Low-frequency Analysis and Recording)谱和DEMON(Detection of Envelope Modulation on Noise)谱的各类改进算法，最新提出的线谱搜索方法是先基于线谱噪声的二阶循环平稳特性进行循环调制谱的计算，然后利用压缩感知方法从循环调制谱中分离线谱，该方法参见“Compressive Sensing for Detecting Shipswith Second-Order Cyclostationary Signatures”，该文2017年9月发表于《IEEEJournal of Oceanic Engineering》第99期，页码为1-13。该方法虽然能够有效的提取线谱，但是循环迭代的收敛速度过慢导致计算时间过长，多次频谱计算会导致频率偏移和频谱泄露，而且利用压缩感知方法会有部分数据在数据压缩的过程中丢失，最终降低轴频检测的精度。
发明内容
要解决的技术问题
为了避免现有技术的不足之处，本发明提出一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法，克服现有的噪声信号识别方法在估计目标轴频时存在的计算量过大、迭代收敛速度慢以及估计精度较差等缺点。
技术方案
一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法，其特征在于步骤如下：
步骤1：对水听器接收到的时域原始信号x(t)做I层的小波包分解，得到最后一层的各段小波包信号，设
Figure GDA0002984751130000021
为第i层小波包分解得到的第j个小波包信号，而原始信号
Figure GDA0002984751130000022
其分解算法为：
Figure GDA0002984751130000023
其中，参数k为位置坐标，gk和hk为小波包分解滤波器的系数，且gk与hk相互正交，满足gk＝(-1)kh1-k；最后一层小波包分解得到的每一段小波信号为
Figure GDA0002984751130000024
步骤2：对最后一层小波包分解得到的每一段小波信号
Figure GDA0002984751130000025
做希尔伯特-黄变换，首先对每一段小波信号
Figure GDA0002984751130000026
进行经验模态分解；然后进行希尔伯特变换求解包络信号，得到连续时域信号为：
Figure GDA0002984751130000027
式中：
Figure GDA0002984751130000028
为卷积运算符；
Figure GDA0002984751130000029
为小波信号
Figure GDA00029847511300000210
经过m次模态分解后得到信号序列；
步骤3：对所有经过希尔伯特-黄变换的小波信号
Figure GDA0002984751130000031
进行重构得到新的包络信号H[r(t)]，其第i层的第j个小波包的重构计算公式如下：
Figure GDA0002984751130000032
其中：i＝I,I-1,...2,，j＝2i,2i-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段，每段的长度含有
Figure GDA0002984751130000033
个数据点，窗函数的重叠率为50％，窗系数为w(n)；因此，对于第l段信号H[r(n)]，进行短时离散傅立叶计算得到时间相关谱定义如下：
Figure GDA0002984751130000034
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做特征值分解，其计算过程如下：
Figure GDA0002984751130000035
其中：λ1,λ2,...λ为矩阵R中p个大的特征矢量，σ2为噪声的方差；用p个大的特征矢量构成信号子空间Us，用N-p个小的特征值对应的特征矢量来构成噪声子空间UN；
利用噪声子空间UN与信号子空间Us构造空间谱：
将PMUSIC在2～15Hz以内的峰值对应的频点记录下来，这些点的频率组成集合F1＝{fq|fq＝fMUSIC,2＜fq＜15,q＝1,2,3…Q}，fMUSIC为空间谱PMUSIC的频率，空间谱公式构造为：
Figure GDA0002984751130000041
其中：a为搜索矢量，aH为a的共轭转置，
Figure GDA0002984751130000042
为UN的共轭转置矩阵；
其中的PMUSIC为2～15Hz频段内的极大值即为疑似轴频
将F1中的每个元素fq视为疑似轴频，分别以fq为中心频率，以0.01＜α＜0.1大小的步长生成Q个集合，每个集合为
Figure GDA0002984751130000043
这里g为中心频率加减步长的个数；
步骤6：对步骤4中计算得到的zl(n,f)的平方再一次做离散傅立叶计算，将时间变量n转换为双循环频率变量f双，推导得到二阶循环相关信号频谱为：
Figure GDA0002984751130000044
其中：f双为集合
Figure GDA0002984751130000045
中的每个疑似轴频；f既是每段时间相关谱zl(n,f)的频率也是循环调制谱中载频f载的中心频率，即上述公式计算得到的循环调制谱被最终表示为P(fq,g,f载)，是关于载频f载和双循环频率变量
Figure GDA0002984751130000046
的频谱；
步骤7：对于只包含单个螺旋桨噪声的原始信号，定义每个集合
Figure GDA0002984751130000047
的频率统计量为Pg＝argmax(P(fq,g,f载))，fq,g是步骤6得到的二阶循环频率分量，定义轴频
Figure GDA0002984751130000048
是max(Pg)对应的频率；
步骤8：以步骤6检测得到的循环调制谱P(fq,g,f载)在1～100Hz频段内的载频f载分别与
Figure GDA0002984751130000049
比较搜索，搜索公式如下：
Figure GDA00029847511300000410
其中：R是对
Figure GDA0002984751130000051
取整后四舍五入的值；
Figure GDA0002984751130000052
是选定的搜索各阶次谐波的偏移精度阈值，取值为小于等于频谱分辨率；
步骤9：将满足公式所得到每个有效频点f载按从小到大排序，得到频率成倍数递增的一系列线谱，为搜索得到的目标螺旋桨的轴频与各阶次的谐波，线谱序列中的第一条谱线就是轴频
Figure GDA0002984751130000053
剩余的线谱f载就是各阶次的谐波。
所述步骤5的p个大的特征矢量的p≥10。
所述步骤5的中心频率加减步长的个数
Figure GDA0002984751130000054
有益效果
本发明提出的一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法，先对原始水听器信号求包络，也就是对该信号进行多尺度小波包分解，然后对每个小波包信号做希尔伯特-黄变换得到小波包络信号，之后将所有的小波包信号进行小波重构，得到改进的噪声包络信号。然后对包络信号进行分段加权计算其时间相关谱，接下来结合多目标分类算法处理时间相关谱信号，取其在低频段峰值突出的线谱作为疑似轴频，根据需求设定小于频率分辨率的步长，以疑似轴频为中心组成疑似轴频集合，再对各疑似轴频的时间相关谱信号做循环相关计算，得到循环调制谱的频率统计量，最终根据该频率统计量来搜索目标的轴频。本发明方法校正了传统的循环调制谱方法在频谱计算过程中产生的轴频偏移误差，提高了轴频检测的精度，极大地减小了循环迭代的计算量，而且在工程应用方面便于实现。
本发明有益效果是：步骤一、二、三改进了传统噪声包络信号求解方法，步骤一中该方法利用尺度函数与小波函数构造小波包函数集合对时域滤波信号做I层的小波包分解，最后得到共2I个小波包分解信号，这些分解信号将频带进行了多层次划分，能够提高信号的时频分辨率；步骤二对这2I个小波包信号进行希尔伯特-黄变换来求解信号的包络，通过多次经验模态分解去除信号中包含的连续谱噪声分量；步骤三中对所有的小波包络信号进行小波包重构，得到改进后的时域噪声包络信号。接下来步骤四、五、六对时域噪声的包络信号进行改进循环调制谱分析，步骤五中该方法利用多目标分类算法处理通过步骤四的计算得到的时间相关谱，将信号在2-15Hz频段的一系列极大值视作疑似轴频，之后以各个疑似轴频为中心频率以适宜步长衍生出多个轴频集合；步骤六中仅仅需要计算轴频集合中的各频率分量的循环调制谱，不需要计算所有频率分量，不仅达到了抑制非二次相位耦合的谐波分量，增强周期性的线谱的目的，还极大减轻了传统循环调制谱的计算量。步骤七与八中统计了每个轴频集合的频率统计量，分别针对单目标信号与多目标信号进行轴频搜索，从中选出符合要求的目标轴频与各阶次谐波。该方法避免了人工挑选疑似轴频，衍生的轴频集合不仅修正了传统的循环调制谱两次频谱计算时带来的误差，还减少了循环调制谱迭代计算所花费的时间。该方法通过实验验证适用于工程上检测目标轴频，即使在低信噪比情况下也具有较高的精度。
附图说明
图1改进噪声包络信号谱的轴频检测流程示意图
图2实测噪声信号的时域波形图及频域波形图，图2(a)为噪声信号的时域波形图，图2(b)为噪声信号的频谱波形图。
图3传统噪声包络信号的时频图(即LOFAR图)。
图4改进噪声包络信号的时频图(即LOFAR图)。
图5多目标分类识别算法得到的疑似轴频的频率空间谱图
图6改进后的噪声信号的循环调制谱图
图7目标螺旋桨的轴频及各阶谐波的频谱图
图8根据轴频计算得到的船速和实际船速的对比
具体实施方式
现结合实施例、附图对本发明作进一步描述：
一种基于改进的循环调制谱方法的螺旋桨轴频搜索方法，其特征在于：首先是通过小波包变换和希尔伯特-黄变换相结合的方式，对经典的包络调制谱中求解包络的步骤加以改进。之后在循环调制谱的第一次频谱计算之后添加多目标分类算法这一步骤，来降低第二次调制谱的计算量，提高轴频检测的精度。具体来说是先对原始信号进行多尺度小波包分解，再对每个小波包信号做希尔伯特-黄变换求包络，然后对所有的小波包络信号重构得到改进的噪声包络信号；之后，改进求解包络信号的循环调制谱，在循环调制谱的两次频谱计算之间添加了一个步骤，即先计算包络信号的时间相关谱，再利用多目标分类算法计算该时间相关谱信号在低频段的局部最大值点，接着以局部最大值点为中心和小于频率分辨率的步长组成疑似轴频集合，之后的循环调制谱仅需要计算集合中的每一个疑似轴频而不是全频段所有频率分量，计算结果作为频率统计量，最终根据该频率统计量来搜索目标的轴频。其改进过程分为以下步骤：
步骤一：先对水听器接收到的时域原始信号x(t)做I层的小波包分解，得到最后一层的各段小波包信号，设
Figure GDA0002984751130000071
为第i层小波包分解得到的第j个小波包信号，而原始信号
Figure GDA0002984751130000072
其分解算法为：
Figure GDA0002984751130000073
其中，参数k为位置坐标，gk和hk为小波包分解滤波器的系数，且gk与hk相互正交，满足gk＝(-1)kh1-k。最后一层小波包分解得到的每一段小波信号为
Figure GDA0002984751130000074
假定小波包变换的尺度函数是
Figure GDA0002984751130000075
小波函数是ψ(x)，小波包函数的定义如下：
Figure GDA0002984751130000076
Figure GDA0002984751130000081
这里，定义函数集合{ψn(t)}为正交尺度函数
Figure GDA0002984751130000082
所确定的小波包。
步骤二：对最后一层小波包分解得到的每一段小波信号
Figure GDA0002984751130000083
做希尔伯特-黄变换。希尔伯特-黄变换操作包含两部分，一是经验模态分解；二是希尔伯特变换求解包络信号。经验模态分解是对输入信号
Figure GDA0002984751130000084
求取极大值点和极小值点，之后对极大值点和极小值点采用三次样条函数插值来构造
Figure GDA0002984751130000085
的上下包络，再计算上下包络的均值函数；
Figure GDA0002984751130000086
减去包络平均后的均值函数得到本征模态分量sv(t)，若sv(t)还存在负的局部极大值和正的局部极小值，则需要继续进行经验模态分解。定义
Figure GDA0002984751130000087
信号减去均值函数得到的本征模态分量sv(t)，然后以sv(t)作为下一次经验模态分解的输入信号不断的进行模态分解，直到满足条件假设条件才会结束，即：
Figure GDA0002984751130000088
其中，sv(t)为每一次模态分解得到的本征模态分量，
Figure GDA0002984751130000089
为m次模态分解后
Figure GDA00029847511300000810
减去所有本征模态分量最终得到的信号序列。当sv(t)同时满足条件：(1)sv(t)的均值函数的平均值趋近于0；(2)极值点个数(包括极大值点个数+极小值点个数)和其同y轴的交点个数之差不能大于1(小于等于1)时，模态分解结束。
小波信号
Figure GDA00029847511300000811
经过m次模态分解，最终得到信号序列
Figure GDA00029847511300000812
将
Figure GDA00029847511300000813
作为输入信号进行希尔伯特变换得到
Figure GDA00029847511300000814
连续时间信号的希尔伯特变换等价于，将信号与冲激响应为
Figure GDA00029847511300000815
的线性系统做卷积运算，对
Figure GDA00029847511300000816
做希尔伯特，其本质是令
Figure GDA00029847511300000817
在频域的频率分量幅度不变，产生90度相移，对于连续时域信号其公式计算如下：
Figure GDA00029847511300000818
式中
Figure GDA00029847511300000819
为卷积运算符。
步骤三：对所有经过希尔伯特-黄变换的小波信号
Figure GDA0002984751130000091
进行重构得到新的包络信号H[r(t)]，其第i层的第j个小波包的重构计算公式如下：
Figure GDA0002984751130000092
其中i＝I,I-1,...2,1，j＝2i,2i-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段，每段的长度含有
Figure GDA0002984751130000093
个数据点，窗函数的重叠率为50％，窗系数为w(n)。因此，对于第l段信号H[r(n)]，进行短时离散傅立叶计算得到时间相关谱定义如下：
Figure GDA0002984751130000094
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做特征值分解，其计算过程如下：
Figure GDA0002984751130000095
其中，λ1,λ2,...λ为矩阵R中p个大的特征矢量，σ2为噪声的方差。用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的频率。空间谱公式构造为：
Figure GDA0002984751130000101
其中a为搜索矢量，aH为a的转置，UN为噪声子空间。PMUSIC在2-15Hz频段内的极大值就是疑似轴频。
将F1中的每个元素fq视为疑似轴频，分别以fq为中心频率，以0.01＜α＜0.1大小的步长生成Q个集合，每个集合为
Figure GDA0002984751130000102
这里g为中心频率加减步长的个数，取值要满足
Figure GDA0002984751130000103
步骤六：对步骤四中计算得到的zl(n,f)的平方再一次做离散傅立叶计算，将时间变量n转换为双循环频率变量f双，最终推导得到二阶循环相关信号频谱如下：
Figure GDA0002984751130000104
其中，f双的取值仅选择集合
Figure GDA0002984751130000105
中的每个疑似轴频，fs是zl(n,f)的采样频率，f既是每段时间相关谱zl(n,f)的频率也是循环调制谱中载频f载的中心频率，即公式(11)计算得到的循环调制谱被最终表示为P(fq,g,f载)，是关于载频f载和双循环频率变量
Figure GDA0002984751130000106
的频谱，该频谱不需要计算每段时间相关谱zl(n,f)在全频域内的所有频率分量，仅需计算疑似轴频，极大减少了传统循环调制谱的计算量。
步骤七：对于只包含单个螺旋桨噪声的原始信号，定义每个集合
Figure GDA0002984751130000107
的频率统计量为Pg＝argmax(P(fq,g,f载))，fq,g是步骤六中计算得到的二阶循环频率分量，定义轴频
Figure GDA0002984751130000111
是max(Pg)对应的频率。
步骤八：理想状态下，轴频对应的各阶次谐波的频率应该是其本身的倍数，即第i阶谐波的频率值为
Figure GDA0002984751130000112
但实际上由于螺旋桨叶片的不均匀导致各阶次谐波的频率并不全是轴频的倍数，存在一定的偏移差。因此当需要识别各阶次谐波的频率时应选择合适的阈值
Figure GDA0002984751130000113
作为偏移精度范围。
此时，用步骤六检测得到的循环调制谱P(fq,g,f载)在1-100Hz频段内的载频f载分别与
Figure GDA0002984751130000114
比较搜索，搜索公式如下：
Figure GDA0002984751130000115
其中：R是对
Figure GDA0002984751130000116
取整后四舍五入的值；
Figure GDA0002984751130000117
是选定的搜索各阶次谐波的偏移精度阈值，取值一般要小于等于频谱分辨率。
步骤九：记录下满足公式(12)搜索得到每个有效频点f载，将其按从小到大排序，可以得到频率成倍数递增的一系列线谱，这就是搜索得到的目标螺旋桨的轴频与各阶次的谐波，该线谱序列中的第一条谱线就是轴频
Figure GDA0002984751130000118
剩余的线谱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变化的离散谱线就是检测得到的频率分量
Figure GDA0002984751130000121
Figure GDA0002984751130000122
具有二阶循环平稳的特性，且保留了轴频与各阶次谐波的周期性。对比图3和图4，图6中对离散线谱的检测精度有了很大提高，且计算量远小于传统的循环调制谱。
图7是本发明方法利用步骤七、八计算频率统计量，得到的目标螺旋桨轴频与各阶次谐波的频谱。如图7所示，图中标注出的10.38Hz是利用本发明方法最终检测到的实测舰船辐射噪声信号中的轴频；其余幅值随频率递增而逐渐降低的线谱为搜索到的各阶次谐波，设定搜索阈值为
Figure GDA0002984751130000131
各阶次谐波根据公式(12)从载频f载中搜索得到。图7中检测到的目标轴频10.38Hz相比图5采用多目标分类算法中得到的极大值频率10.99Hz有0.61Hz的偏移校正。
图8是采用全球定位系统(GPS)测定的实验船航速与采用本发明方法检测所得轴频计算出的实验船航速对比图，实际上是利用实测数据来验证本发明提出的轴频检测方法的有效性。具体实施过程是先通过全球定位系统(GPS)得到不同航速下实验船的辐射噪声数据，记录下当前实验船的航速，再利用本方法计算实验船的辐射噪声信号得到轴频估计值，然后通过公式(13)计算得到当前实验船的瞬时估计航速，最后与全球定位系统记录下的航速对比。
Figure GDA0002984751130000132
式中：f轴为轴频，V为航速，D为实验船的螺旋桨直径，J为螺旋桨的进速系数，每条船的D、J都是固定不变的。
图中离散的点“*”代表本发明方法从108个辐射噪声数据中搜索得到轴频后，依据公式(13)计算出对应时刻的实验船航速，“o”代表全球定位系统测定的整条航线上的实验船航速变化。两者之间吻合度非常高，证明了本发明方法对于实测舰船噪声数据的有效性。
本发明在典型实施例中取得了明显的实施效果，与现有技术相比其优越性在于无需频率误差校正，轴频搜索精确且计算量小，便于应用。

Claims (3)

1.一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法，其特征在于步骤如下：
步骤1：对水听器接收到的时域原始信号x(t)做I层的小波包分解，得到最后一层的各段小波包信号，设
Figure FDA0002984751120000011
为第i层小波包分解得到的第j个小波包信号，而原始信号
Figure FDA0002984751120000012
其分解算法为：
Figure FDA0002984751120000013
其中，参数k为位置坐标，gk和hk为小波包分解滤波器的系数，且gk与hk相互正交，满足gk＝(-1)kh1-k；最后一层小波包分解得到的每一段小波信号为
Figure FDA0002984751120000014
步骤2：对最后一层小波包分解得到的每一段小波信号
Figure FDA0002984751120000015
做希尔伯特-黄变换，首先对每一段小波信号
Figure FDA0002984751120000016
进行经验模态分解；然后进行希尔伯特变换求解包络信号，得到连续时域信号为：
Figure FDA0002984751120000017
式中：
Figure FDA0002984751120000018
为卷积运算符；
Figure FDA0002984751120000019
为小波信号
Figure FDA00029847511200000110
经过m次模态分解后得到信号序列；
步骤3：对所有经过希尔伯特-黄变换的小波信号
Figure FDA00029847511200000111
进行重构得到新的包络信号H[r(t)]，其第i层的第j个小波包的重构计算公式如下：
Figure FDA00029847511200000112
其中：i＝I,I-1,...2,，j＝2i,2i-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段，每段的长度含有
Figure FDA0002984751120000021
个数据点，窗函数的重叠率为50％，窗系数为w(n)；因此，对于第l段信号H[r(n)]，进行短时离散傅立叶计算得到时间相关谱定义如下：
Figure FDA0002984751120000022
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做特征值分解，其计算过程如下：
Figure FDA0002984751120000023
其中：λ1,λ2,...λ为矩阵R中p个大的特征矢量，σ2为噪声的方差；用p个大的特征矢量构成信号子空间Us，用N-p个小的特征值对应的特征矢量来构成噪声子空间UN；
利用噪声子空间UN与信号子空间Us构造空间谱：
将PMUSIC在2～15Hz以内的峰值对应的频点记录下来，这些点的频率组成集合F1＝{fq|fq＝fMUSIC,2＜fq＜15,q＝1,2,3…Q}，fMUSIC为空间谱PMUSIC的频率，空间谱公式构造为：
Figure FDA0002984751120000024
其中：a为搜索矢量，aH为a的共轭转置，
Figure FDA0002984751120000025
为UN的共轭转置矩阵；
其中的PMUSIC为2～15Hz频段内的极大值即为疑似轴频
将F1中的每个元素fq视为疑似轴频，分别以fq为中心频率，以0.01＜α＜0.1大小的步长生成Q个集合，每个集合为
Figure FDA0002984751120000026
这里g为中心频率加减步长的个数；
步骤6：对步骤4中计算得到的zl(n,f)的平方再一次做离散傅立叶计算，将时间变量n转换为双循环频率变量f双，推导得到二阶循环相关信号频谱为：
Figure FDA0002984751120000031
其中：fq,g为集合
Figure FDA0002984751120000032
中的每个疑似轴频；f既是每段时间相关谱zl(n,f)的频率也是循环调制谱中载频f载的中心频率，即上述公式计算得到的循环调制谱被最终表示为P(fq,g,f载)，是关于载频f载和双循环频率变量fq,g的频谱；
步骤7：对于只包含单个螺旋桨噪声的原始信号，定义每个集合
Figure FDA0002984751120000033
的频率统计量为Pg＝arg max(P(fq,g,f载))，fq,g是步骤6得到的二阶循环频率分量，定义轴频
Figure FDA0002984751120000034
是max(Pg)对应的频率；
步骤8：以步骤6检测得到的循环调制谱P(fq,g,f载)在1～100Hz频段内的载频f载分别与
Figure FDA0002984751120000035
比较搜索，搜索公式如下：
Figure FDA0002984751120000036
其中：R是对
Figure FDA0002984751120000037
取整后四舍五入的值；
Figure FDA0002984751120000038
是选定的搜索各阶次谐波的偏移精度阈值，取值为小于等于频谱分辨率；
步骤9：将满足公式所得到每个有效频点f载按从小到大排序，得到频率成倍数递增的一系列线谱，为搜索得到的目标螺旋桨的轴频与各阶次的谐波，线谱序列中的第一条谱线就是轴频
Figure FDA0002984751120000039
剩余的线谱f载就是各阶次的谐波。
2.根据权利要求1所述基于改进噪声包络信号识别的螺旋桨轴频搜索方法，其特征在于：所述步骤5的p个大的特征矢量的p≥10。
3.根据权利要求1所述基于改进噪声包络信号识别的螺旋桨轴频搜索方法，其特征在于：所述步骤5的中心频率加减步长的个数
Figure FDA0002984751120000041
CN201810487439.XA 2018-05-21 2018-05-21 一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法 Active CN108921014B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810487439.XA CN108921014B (zh) 2018-05-21 2018-05-21 一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810487439.XA CN108921014B (zh) 2018-05-21 2018-05-21 一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法

Publications (2)

Publication Number Publication Date
CN108921014A CN108921014A (zh) 2018-11-30
CN108921014B true CN108921014B (zh) 2021-05-14

Family

ID=64403926

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810487439.XA Active CN108921014B (zh) 2018-05-21 2018-05-21 一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法

Country Status (1)

Country Link
CN (1) CN108921014B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109840495B (zh) * 2019-01-29 2020-09-18 浙江大学 一种低信噪比条件下的轴频线谱增强方法
CN110596458B (zh) * 2019-07-16 2021-02-02 西北工业大学 Demon谱谐波线谱和基频自动估计方法
CN110798419A (zh) * 2019-10-28 2020-02-14 北京邮电大学 一种调制方式识别方法及装置
CN113095113B (zh) * 2019-12-23 2024-04-09 中国科学院声学研究所 一种用于水下目标识别的小波线谱特征提取方法及系统
CN111160207B (zh) * 2019-12-24 2023-08-15 浙江大学 一种基于辐射噪声调制的桨叶数特征提取方法
CN111368679B (zh) * 2020-02-26 2022-03-22 西北工业大学 一种低秩矩阵分解的谱线检测方法
CN111735525B (zh) * 2020-05-28 2023-03-31 哈尔滨工程大学 一种适用于无人声纳的demon谱特征提取方法
CN114757242B (zh) * 2022-06-16 2022-09-23 中国空气动力研究与发展中心低速空气动力研究所 基于循环维纳滤波的直升机噪声增强方法以及检测方法
CN117008863B (zh) * 2023-09-28 2024-04-16 之江实验室 一种lofar长数据处理及显示方法和装置
CN118471264A (zh) * 2024-07-09 2024-08-09 威海晟南船舶技术服务有限公司 船舶目标识别方法及装置

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101539597A (zh) * 2009-04-29 2009-09-23 哈尔滨工程大学 一种分辨具有相同频带辐射噪声的多目标方法
CN102183366A (zh) * 2011-03-08 2011-09-14 上海大学 滚动轴承振动测量和故障分析装置及方法
CN104359674A (zh) * 2014-10-20 2015-02-18 广东电网有限责任公司电力科学研究院 基于时域与频域状态监测的高速滚动轴承故障诊断方法
CN106483520A (zh) * 2016-09-27 2017-03-08 哈尔滨工程大学 一种船舶辐射噪声调制系数估计方法
CN107024718A (zh) * 2017-05-31 2017-08-08 西南石油大学 基于ceemd‑spwvd时频谱分析的叠后地震流体预测方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7872948B2 (en) * 2008-04-14 2011-01-18 The Boeing Company Acoustic wide area air surveillance system
US8612182B2 (en) * 2010-04-13 2013-12-17 General Electric Company Methods and systems for isolating a frequency in a rotating machine

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101539597A (zh) * 2009-04-29 2009-09-23 哈尔滨工程大学 一种分辨具有相同频带辐射噪声的多目标方法
CN102183366A (zh) * 2011-03-08 2011-09-14 上海大学 滚动轴承振动测量和故障分析装置及方法
CN104359674A (zh) * 2014-10-20 2015-02-18 广东电网有限责任公司电力科学研究院 基于时域与频域状态监测的高速滚动轴承故障诊断方法
CN106483520A (zh) * 2016-09-27 2017-03-08 哈尔滨工程大学 一种船舶辐射噪声调制系数估计方法
CN107024718A (zh) * 2017-05-31 2017-08-08 西南石油大学 基于ceemd‑spwvd时频谱分析的叠后地震流体预测方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Sea trial researches on extraction of propeller shaft frequency;Niu F.等;《Proceedings 2013 International Conference on Mechatronic Sciences, Electric Engineering and Computer (MEC)》;20140828;第1306-1309页 *
低信噪比下的螺旋桨轴频估计方法;陈敬军 等;《信号处理》;20041031;第20卷(第5期);第461-464页 *

Also Published As

Publication number Publication date
CN108921014A (zh) 2018-11-30

Similar Documents

Publication Publication Date Title
CN108921014B (zh) 一种基于改进噪声包络信号识别的螺旋桨轴频搜索方法
CN108548957B (zh) 基于循环调制频谱和分段互相关相结合的双谱分析方法
CN111024209B (zh) 一种适用于矢量水听器的线谱检测方法
CN111222088B (zh) 一种改进的平顶自卷积窗加权电力谐波幅值估计方法
CN106199532B (zh) 基于混合傅立叶-小波分析的探地雷达信号降噪方法
CN104142425B (zh) 一种正弦信号频率估计的相位匹配方法
CN102393423A (zh) 基于自适应阈值正交小波变换的兰姆波去噪方法
CN112183225B (zh) 一种基于概率潜在语义分析的水下目标信号特征提取方法
CN114785379A (zh) 一种水声janus信号参数估计方法及系统
CN111856401A (zh) 一种基于互谱相位拟合的时延估计方法
CN113899444A (zh) 一种基于汉宁双窗的振弦传感器共振频率测量方法
Lu et al. Fundamental frequency detection of underwater acoustic target using DEMON spectrum and CNN network
CN111175727B (zh) 一种基于条件波数谱密度的宽带信号方位估计的方法
Niu et al. Mode separation with one hydrophone in shallow water: A sparse Bayesian learning approach based on phase speed
CN109101462A (zh) 可降低时间损耗的傅里叶变换算法
RU2464588C1 (ru) Устройство обнаружения шумовых гидроакустических сигналов в виде звукоряда на основе вычисления интегрального вейвлет-спектра
Leiber et al. Differentiable short-time Fourier transform with respect to the hop length
CN112505640B (zh) 基于参数自适应的扩展b分布脉冲信号时频分析方法
CN110865375A (zh) 一种水中目标检测方法
CN116612735A (zh) 一种水声音频去噪方法
CN115755106A (zh) 一种线性调频干扰的扫频周期估计方法
CN112162235B (zh) 一种光滑分段随机共振增强的声矢量信号定向方法
CN115015942A (zh) 一种自适应水下目标声学测速装置及方法
CN107957571A (zh) 水听器测向方法、装置、计算机可读存储介质及计算机设备
CN115166648B (zh) 一种低信噪比雷达信号处理方法及装置

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant