CN106595836B - 水下双声源频谱混叠情况下的频率提取方法 - Google Patents

Abstract

水下双声源频谱混叠情况下的频率提取方法，涉及水下声源的频率提取方法。本发明为了解决目前还没有一种能够针对水下多声源在频谱相互混叠的情况进行频率分离和提取的有效方法的问题。本发明首先采集含有水下双声源频率信息的水面振动信息，得到含有水下双声源频率信息的水面振动信息电信号；然后将所述的水面振动信息电信号经过去噪滤波处理后进行希尔伯特变换，将水面振动信息电信号转换为解析信号，然后取解析信号的模值，再取模值的平方；最后进行平滑处理后再进行傅里叶变换，得到频谱分析图，根据频谱分析图上的频移值得出水下两个声源的频率信息。本发明适用于水下声源的频率提取。

Description

水下双声源频谱混叠情况下的频率提取方法

技术领域

本发明涉及水下声源的频率提取方法。

背景技术

水下声信号的特征提取一直是海洋工程和国防的热点研究问题之一。不管是民用方面的水下工程建设，海洋资源开发还是军事方面的对潜通信、作战目标定位等都以水下声信号特征的提取为基础。发展水下声信号的特征提取技术具有非常重要的战略意义。

水下声信号的众多特征中，声信号的频率特征提取一直深受学者们的关注。在以往公开报道的文献中，水下声信号频率提取主要针对两类声信号，一是单一频率的水下声信号频率提取，二是频带完全分离的多频率水下声信号频率提取。目前已经有很多有效的手段实现频带完全分离的多频率水下声信号频率提取，但是目前还没有一种能够针对水下多声源在频谱相互混叠的情况进行频率分离和提取的有效方法。

发明内容

本发明为了解决目前还没有一种能够针对水下多声源在频谱相互混叠的情况进行频率分离和提取的有效方法的问题。

水下双声源频谱混叠情况下的频率提取方法，是基于激光干涉法探测水下声源频率的探测装置实现的，所述的激光干涉法探测水下声源频率的探测装置光路如下：

激光器射出一束激光，激光与立方体半反射镜的镜片呈45度角射到立方体半反射镜上，激光发生反射和透过，分别得到反射激光和透射激光；所述的反射激光透过一片与反射激光垂直的衰减片垂直射到反射镜上，然后经过反射再次垂直穿过衰减片达到立方体半反射镜上，记为参考光；所述的透射激光射到与透射激光呈45度角的反射镜上，然后透射激光射到水面上，透射激光经过水面反射按照原光路返回到立方体半反射镜上，记为测量光；返回到立方体半反射镜上的参考光和测量光发生干涉后被光电接收器接收，并通过数据采集模块到达上位机；

本发明所述的频率提取方法包括以下步骤：

步骤1：利用激光干涉法探测水下声源频率的探测装置采集含有水下双声源频率信息的水面振动信息；将含有水下双声源频率信息的水面振动信息经过光电接收器接收，然后经过数据采集卡转换，得到含有水下双声源频率信息的水面振动信息电信号；

步骤2、将所述的水面振动信息电信号经过去噪滤波处理后进行希尔伯特变换，将水面振动信息电信号转换为解析信号，然后取解析信号的模值，再取模值的平方；

步骤3、将步骤2取模值平方后的信号平滑处理后再进行傅里叶变换，得到频谱分析图，根据频谱分析图上的频移值得出水下两个声源的频率信息。

优选地，步骤1中所述的含有水下双声源频率信息的水面振动信息经过光电接收器接收后的电信号U表示为：

其中，U₀为电信号幅值；k为波数，k＝2π/λ，λ为激光波长；ΔL为水面静止时参考光和测量光的初始光程差；t为时间；A_s、f_s、分别为自然水表面波的幅值、频率和初始相位；s为自然水表面波的个数；A₁、A₂为水下声源引起的水表面声波的振幅，f₁、f₂为水下声源频率，为水下声源初始相位。

优选地，所述步骤2的具体过程包括以下步骤：

针对

为了方便表述，设

第一类的N阶贝塞尔函数J_N(x)的展开式如下：

其中，x表示自变量；ω为角频率，i为正整数；J₀(x)表示0阶贝塞尔函数，J_2i-1(x)表示(2i-1)阶贝塞尔函数；

利用三角函数公式和第一类的N阶贝塞尔函数J_N(x)的展开式对步骤一中的水面振动信息电信号进行整理，得到：

其中，ω₁＝2πf₁、ω₂＝2πf₂分别表示水下声源的角频率；m、n为正整数；

结合希尔伯特变换以及希尔伯特变换调制性质，对上式进行希尔伯特变换：

得到信号的解析形式：

其中，H(U)表示希尔伯特变换的信号，Z(U)表示解析信号；

取解析信号模值的平方：|Z(U)|²＝U²+H(U)²；并利用第一类的n阶贝塞尔函数J_n(x)的展开式展开cosφ和sinφ；为了方便表述，设水表面波仅存在一种频率的波动，即φ＝2k(△L+A_n sinω_nt)；其中，ω_n为自然水表面波的角频率，ω_n＝2πf_n；

多个频率的波动只是在此基础上的简单叠加；Z(U)为U的函数，U为t的函数，将信号表示为以下形式：

其中，第一项表示以ω₁为中心，左右频移iω_n的频带；第二项表示以ω₂为中心，左右频移jω_n的频带；第三项是ω₁和ω₂的差频项和倍频项组成的频带；M_i、N_j′、Q_p分别为相应项的幅值调制函数，i、j′、p均表示正整数；M_i、N_j′的值受自然水表面波幅值的影响，且随着i、j′值的增大，急速减小，使得以ω₁和ω₂为中心，左右频移的项很少，从而使得以ω₁和ω₂为中心的频带很窄，两个频带得意区分；Q_p相对M_i、N_j′很小，忽略不计；根据傅里叶变化的性质，以ω₁和ω₂为中心的频带分别关于ω₁和ω₂左右对称，因此通过频移值计算便可得到ω₁和ω₂的值。

优选地，步骤3中所述的根据频谱分析图上的频移值得出水下两个声源的频率信息的过程包括以下步骤：

在频谱分析图中，以信号频率为中心，很窄的频带范围内存在较大幅值，通过频移值计算便能够得到ω₁和ω₂的值。

优选地，步骤2中结合希尔伯特变换以及希尔伯特变换调制性质对上式进行希尔伯特变换所述的希尔伯特变换通过以下希尔伯特变换对实现：

H(cosωt)＝sinωt和H(sinωt)＝-cosωt。

优选地，步骤2中结合希尔伯特变换以及希尔伯特变换调制性质对上式进行希尔伯特变换所述的希尔伯特变换调制性质如下：

两个信号为f(t)、g(t)，分别为f(t)、g(t)的Fourier变换；若使得当|ω|>a时，当|ω|<a时，则两个信号乘积的Hilbert变换为：H[f(x)g(x)]＝f(x)H[g(x)]；f(x)相对于g(x)变化缓慢时，两个信号乘积的Hilbert变换取决于高频信号。

本发明具有以下有益效果：

本发明利用激光干涉法探测水下声源频率的探测装置采集到干涉信号，对干涉信号进行处理后，可以分离水下两个声源频率接近时导致频谱混叠在一起的频谱，使得两个频带得以区分。仿真实验结果表明：当水下两个声源的频率间隔大于最大自然水表面波频率时，本发明的频率提取方法能够准确的分离区分两个频带，准确提取出两个声信号频率。

附图说明

图1为激光干涉法探测水下声源频率的探测装置光路图；

图2为水面振动信息电信号时域波形图；

图3为水面振动信息时域波形图经过傅里叶变换得到的信号频谱图；

图4为希尔伯特变换后再取解析信号模值平方后的信号时域波形图；

图5为经过平滑处理后再进行傅里叶变换得到的信号频谱图；

图6为图5的局部放大图；

图7为信号发生器设置频率为4kHz和4.1kHz时，上位机采集到的信号时域波形图；

图8为信号频域图；

图9为经过本方法处理后得到的信号频域分析图；

图10为图9的局部放大图。

具体实施方式

具体实施方式一：

水下双声源频谱混叠情况下的频率提取方法，是基于激光干涉法探测水下声源频率的探测装置实现的，如图1所示，所述的激光干涉法探测水下声源频率的探测装置光路如下：

激光器1射出一束激光，激光与立方体半反射镜2的镜片呈45度角射到立方体半反射镜2上，激光发生反射和透过，分别得到反射激光和透射激光；所述的反射激光透过一片与反射激光垂直的衰减片3垂直射到反射镜4上，然后经过反射再次垂直穿过衰减片达到立方体半反射镜上，记为参考光；所述的透射激光射到与透射激光呈45度角的反射镜5上，然后透射激光射到水面上，透射激光经过水面反射按照原光路返回到立方体半反射镜上，记为测量光；返回到立方体半反射镜上的参考光和测量光发生干涉后被光电接收器6接收，并通过数据采集模块7到达上位机8；

当水中存在声源A和声源B时，声源A和声源B的波动信号会与水面波动信号叠加，利用以上光路装置和本实施方式的方法进行声源A和声源B的频率提取。

本实施方式所述的频率提取方法包括以下步骤：

步骤1：利用激光干涉法探测水下声源频率的探测装置采集含有水下双声源频率信息的水面振动信息；将含有水下双声源频率信息的水面振动信息经过光电接收器接收，然后经过数据采集卡转换，得到含有水下双声源频率信息的水面振动信息电信号；

步骤2、将所述的水面振动信息电信号经过去噪滤波处理后进行希尔伯特变换，将水面振动信息电信号转换为解析信号，然后取解析信号的模值，再取模值的平方；

步骤3、将步骤2取模值平方后的信号平滑处理后再进行傅里叶变换，得到频谱分析图，根据频谱分析图上的频移值得出水下两个声源的频率信息。

具体实施方式二：

本实施方式步骤1中所述的含有水下双声源频率信息的水面振动信息经过光电接收器接收后的电信号U表示为：

其中，U₀为电信号幅值；k为波数，k＝2π/λ，λ为激光波长；ΔL为水面静止时参考光和测量光的初始光程差；t为时间；A_s、f_s、分别为自然水表面波的幅值、频率和初始相位；s为自然水表面波的个数；A₁、A₂为水下声源引起的水表面声波的振幅，f₁、f₂为水下声源频率，为水下声源初始相位。

其他步骤和参数与具体实施方式一相同。

具体实施方式三：

本实施方式所述步骤2的具体过程包括以下步骤：

针对

为了方便表述，设

第一类的N阶贝塞尔函数J_N(x)的展开式如下：

其中，x表示自变量；ω为角频率，i为正整数；J₀(x)表示0阶贝塞尔函数，J_2i-1(x)表示(2i-1)阶贝塞尔函数；

利用三角函数公式和第一类的N阶贝塞尔函数J_N(x)的展开式对步骤一中的水面振动信息电信号进行整理，得到：

其中，ω₁＝2πf₁、ω₂＝2πf₂分别表示水下声源的角频率；m、n为正整数；

从上式可以看出，信号的频谱结构如下：在低频段存在由自然水表面波引起的频带，在高频段存在多个频带，分别以nω₁、mω₂、nω₁±mω₂为中心左右频移对称，其频移宽度受自然水表面波影响，幅值不均匀，由贝塞尔函数值决定；当ω₁和ω₂比较接近时，则频谱混叠在一起，频谱不再遵循以信号频率为中心左右对称，且各频移分量幅值无明显特征，因而无法识别区分两个频带，从而无法提取出两个声源频率信息；

结合希尔伯特变换以及希尔伯特变换调制性质，对上式进行希尔伯特变换：

得到信号的解析形式：

其中，H(U)表示希尔伯特变换的信号，Z(U)表示解析信号；

取解析信号模值的平方：|Z(U)|²＝U²+H(U)²；并利用第一类的n阶贝塞尔函数J_n(x)的展开式展开cosφ和sinφ；为了方便表述，设水表面波仅存在一种频率的波动，即φ＝2k(△L+A_n sinω_nt)；其中，ω_n为自然水表面波的角频率，ω_n＝2πf_n；

多个频率的波动只是在此基础上的简单叠加；Z(U)为U的函数，U为t的函数，将信号表示为以下形式：

其中，第一项表示以ω₁为中心，左右频移iω_n的频带；第二项表示以ω₂为中心，左右频移jω_n的频带；第三项是ω₁和ω₂的差频项和倍频项组成的频带；M_i、N_j′、Q_p分别为相应项的幅值调制函数，i、j′、p均表示正整数；M_i、N_j′的值受自然水表面波幅值的影响，且随着i、j′值的增大，急速减小，使得以ω₁和ω₂为中心，左右频移的项很少，从而使得以ω₁和ω₂为中心的频带很窄，两个频带得意区分；Q_p相对M_i、N_j′很小，忽略不计；根据傅里叶变化的性质，以ω₁和ω₂为中心的频带分别关于ω₁和ω₂左右对称，因此通过频移值计算便可得到ω₁和ω₂的值。

其他步骤和参数与具体实施方式一或二相同

具体实施方式四：

本实施方式步骤3中所述的根据频谱分析图上的频移值得出水下两个声源的频率信息的过程包括以下步骤：

在频谱分析图中，以信号频率为中心，很窄的频带范围内存在较大幅值，通过频移值计算便能够得到ω₁和ω₂的值。

其他步骤和参数与具体实施方式一至三之一相同

具体实施方式五：

本实施方式步骤2中结合希尔伯特变换以及希尔伯特变换调制性质对上式进行希尔伯特变换所述的希尔伯特变换通过以下希尔伯特变换对实现：

H(cosωt)＝sinωt和H(sinωt)＝-cosωt。

其他步骤和参数与具体实施方式二至三之一相同

具体实施方式五：

本实施方式步骤2中结合希尔伯特变换以及希尔伯特变换调制性质对上式进行希尔伯特变换所述的希尔伯特变换调制性质如下：

设两个信号为f(t)、g(t)，分别为f(t)、g(t)的Fourier变换；若(a为常数)，使得当|ω|>a时，当|ω|<a时，则两个信号乘积的Hilbert变换为：H[f(x)g(x)]＝f(x)H[g(x)]；f(x)相对于g(x)变化缓慢时，两个信号乘积的Hilbert变换取决于高频信号。

其他步骤和参数与具体实施方式二至四之一相同

实施例

利用MATLAB软件进行仿真计算，目的在于验本发明的正确性及效果。参数设置如表1所示。

表1

图2为水面振动信息电信号时域波形图，由于自然水表面波幅值(微米量级)远远大于水下声源引起的水表面声波幅值(纳米量级)，因此波形图中整体明暗相间变化的干涉条纹是由自然水表面波引起的，加载在干涉条纹上相位高速变化的条纹是由水下声源引起的。图3为水面振动信息时域波形图经过傅里叶变换得到的信号频谱图，从图中可以看出，两个频带混叠在一起，频带不再遵循关于信号频率为中心左右对称，因而无法区分两个频带，提取出水下两个声源信号的频率。

将水面振动信息时域波形图经过希尔伯特变换之后得到信号的解析形式，取解析信号模值的平方得到如图4所示的变换后的信号时域波形图。与原始信号相比，其波形去掉了趋势项。图5为将图4所示的波形经过平滑处理后再进行傅里叶变换得到的信号频谱图，从图中可以看出在信号频率处存在两个明显很窄的频带。图6为图5的局部放大图，从图中可以清楚的看出两个很窄其明显分离的频带，分别以设置的两个信号频率为中心左右对称。利用MATLAB软件编程提取出两个频带的中心，得到信号频率为4.001kHz和4.102kHz。与设置的信号频率基本一致，验证了理论分析的正确性及本方法信号提取的可行性。

为了证明探测到的水下声源频率与信号发生器设置的信号频率的一致性，从而验证本方法的提取水下双声源频谱混叠情况下的频率的可行性，进行了多组对比实验。图7为信号发生器设置频率为4kHz和4.1kHz时，上位机采集到的信号时域波形图，图8为信号频域图，图9为经过本方法处理后得到的信号频域分析图，图10为图9的局部放大图。各个信号波形图与理论分析和仿真结果相一致。各组实验数据、程序解调结果及误差分析如表2所示。

表2

从多组实验数据可以看出，最大标准偏差为2.3866Hz，选取3组数据的平均值作为测量值，与频率发生器设置的信号频率进行对比，最大相对误差为0.052％。实验证明本方法能有效提取水下双声源频谱混叠情况下的两个声源频率。

Claims (6)

1.水下双声源频谱混叠情况下的频率提取方法，是基于激光干涉法探测水下声源频率的探测装置实现的，所述的激光干涉法探测水下声源频率的探测装置光路如下：

激光器射出一束激光，激光与立方体半反射镜的镜片呈45度角射到立方体半反射镜上，激光发生反射和透过，分别得到反射激光和透射激光；所述的反射激光透过一片与反射激光垂直的衰减片垂直射到反射镜上，然后经过反射再次垂直穿过衰减片达到立方体半反射镜上，记为参考光；所述的透射激光射到与透射激光呈45度角的反射镜上，然后透射激光射到水面上，透射激光经过水面反射按照原光路返回到立方体半反射镜上，记为测量光；返回到立方体半反射镜上的参考光和测量光发生干涉后被光电接收器接收，并通过数据采集模块到达上位机；

其特征在于，所述的频率提取方法包括以下步骤：

步骤1：利用激光干涉法探测水下声源频率的探测装置采集含有水下双声源频率信息的水面振动信息；将含有水下双声源频率信息的水面振动信息经过光电接收器接收，然后经过数据采集卡转换，得到含有水下双声源频率信息的水面振动信息电信号；

步骤2、将所述的水面振动信息电信号经过去噪滤波处理后进行希尔伯特变换，将水面振动信息电信号转换为解析信号，然后取解析信号的模值，再取模值的平方；

步骤3、将步骤2取模值平方后的信号平滑处理后再进行傅里叶变换，得到频谱分析图，根据频谱分析图上的频移值得出水下两个声源的频率信息。

2.根据权利要求1所述的水下双声源频谱混叠情况下的频率提取方法，其特征在于，步骤1中所述的含有水下双声源频率信息的水面振动信息经过光电接收器接收后的电信号U表示为：

其中，U₀为电信号幅值；k为波数，k＝2π/λ，λ为激光波长；ΔL为水面静止时参考光和测量光的初始光程差；t为时间；A_s、f_s、分别为自然水表面波的幅值、频率和初始相位；s为自然水表面波的个数；A₁、A₂为水下声源引起的水表面声波的振幅，f₁、f₂为水下声源频率，为水下声源初始相位。

3.根据权利要求2所述的水下双声源频谱混叠情况下的频率提取方法，其特征在于，所述步骤2的具体过程包括以下步骤：

针对

为了方便表述，设

第一类的N阶贝塞尔函数J_N(x)的展开式如下：

其中，x表示自变量；ω为角频率，i为正整数；J₀(x)表示0阶贝塞尔函数，J_2i-1(x)表示(2i-1)阶贝塞尔函数；

利用三角函数公式和第一类的N阶贝塞尔函数J_N(x)的展开式对步骤一中的水面振动信息电信号进行整理，得到：

其中，ω₁＝2πf₁、ω₂＝2πf₂分别表示水下声源的角频率；m、n为正整数；

结合希尔伯特变换以及希尔伯特变换调制性质，对上式进行希尔伯特变换：

得到信号的解析形式：

其中，H(U)表示希尔伯特变换的信号，Z(U)表示解析信号；

取解析信号模值的平方：|Z(U)|²＝U²+H(U)²；并利用第一类的n阶贝塞尔函数J_n(x)的展开式展开cosφ和sinφ；为了方便表述，设水表面波仅存在一种频率的波动，φ＝2k(ΔL+A_n sinω_nt)；其中，ω_n为自然水表面波的角频率，ω_n＝2πf_n；

多个频率的波动只是在此基础上的简单叠加；Z(U)为U的函数，U为t的函数，将信号表示为以下形式：

其中，第一项表示以ω₁为中心，左右频移iω_n的频带；第二项表示以ω₂为中心，左右频移jω_n的频带；第三项是ω₁和ω₂的差频项和倍频项组成的频带；M_i、N_j′、Q_p分别为相应项的幅值调制函数，i、j′、p均表示正整数；M_i、N_j′的值受自然水表面波幅值的影响，且随着i、j′值的增大，急速减小，使得以ω₁和ω₂为中心，左右频移的项很少，从而使得以ω₁和ω₂为中心的频带很窄，两个频带得意区分；Q_p相对M_i、N_j′很小，忽略不计；根据傅里叶变化的性质，以ω₁和ω₂为中心的频带分别关于ω₁和ω₂左右对称，因此通过频移值计算便可得到ω₁和ω₂的值。

4.根据权利要求3所述的水下双声源频谱混叠情况下的频率提取方法，其特征在于，步骤3中所述的根据频谱分析图上的频移值得出水下两个声源的频率信息的过程包括以下步骤：

在频谱分析图中，以信号频率为中心，很窄的频带范围内存在较大幅值，通过频移值计算便能够得到ω₁和ω₂的值。

5.根据权利要求3或4所述的水下双声源频谱混叠情况下的频率提取方法，其特征在于，步骤2中结合希尔伯特变换以及希尔伯特变换调制性质对式子(4)进行希尔伯特变换，所述的希尔伯特变换通过以下希尔伯特变换对实现：

H(cosωt)＝sinωt和H(sinωt)＝-cosωt。

6.根据权利要求5所述的水下双声源频谱混叠情况下的频率提取方法，其特征在于，所述的希尔伯特变换调制性质如下：

两个信号为f(t)、g(t)，分别为f(t)、g(t)的Fourier变换；若使得当|ω|>a时，当|ω|<a时，则两个信号乘积的Hilbert变换为：H[f(x)g(x)]＝f(x)H[g(x)]；f(x)相对于g(x)变化缓慢时，两个信号乘积的Hilbert变换取决于高频信号。