基于小波脊的水下声信号瞬时频率解调方法

摘要

基于小波脊的水下声信号瞬时频率解调方法，涉及水下声信号瞬时频率的解调，属于水下声信号探测技术领域。本发明解决了现有的常规的信号解调方法是频谱分析方法，频谱分析方法不能确定水下声信号的瞬时频率；计算水下声信号的频率时，对频率的探测精度依赖于频谱阈值的设定，常常出现粗大误差的问题。本发明的技术方案为：对信号小波变换得到小波变换系数；根据局部模极大值方法得到各个时刻小波脊点的分布位置；根据小波脊点的分布位置计算出各个时刻的小波脊中心尺度；根据小波中心尺度反求各个时刻的瞬时频率。本发明的适用领域：水下发声目标发声频率的探测、液体表面声波频率探测、微小振幅机械振动频率的测量。

说明书

技术领域

本发明涉及水下声信号瞬时频率的解调，尤其涉及一种基于小波脊的水下声信号瞬时 频率解调方法，属于水下声信号探测技术领域。

背景技术

水下声源向外辐射声波，能在水表面形成横向传播的水表面声波，水表面声波的振动 频率与水下声源的发声频率一致，因此水表面声波携带了水下声信号的相关信息。因此， 由于水下目标探测的需要，越来越重视对水表面声波的探测和利用。激光干涉方法能够从 动态水表面中获取水下声信号的频率信息，利用激光干涉来探测水表面声波以获得水下声 信号的信息逐渐成为研究热点。目前利用激光多普勒效应探测水下声信号已成为获取水下 声信号的重要手段之一。

然而当前利用激光干涉方法的水下声信号探测仅限于对稳频信号的探测上，实际应用 中，水下声信号的发声频率往往是时变的，常规的信号解调方法是频谱分析方法，该方法 不能确定水下声信号的瞬时频率；不仅如此，通过频谱分析方法计算水下声信号的频率时， 这种方法对频率的探测精度依赖于频谱阈值的设定，常常出现粗大误差。

发明内容

本发明的目的是提出一种基于小波脊的水下声信号瞬时频率解调方法，以解决针对现 有的常规的信号解调方法是频谱分析方法，频谱分析方法不能确定水下声信号的瞬时频 率；计算水下声信号的频率时，对频率的探测精度依赖于频谱阈值的设定，常常出现粗大 误差的问题。

本发明为解决上述技术问题所采用的技术方案是：

本发明所述的一种基于小波脊的水下声信号瞬时频率解调方法，是按照以下步骤实现 的：

步骤一、利用迈克尔逊激光干涉系统探测水表面；

步骤二、采集激光干涉信号，所述采集过程为：干涉信号被光电探测器接收，并将其 转化为电信号，该电信号通过数据采集卡采集后，传送给上位机；

步骤三、对采集到的干涉信号进行频谱分析；

步骤四、确定小波变换的观测尺度；

步骤五、对干涉信号进行小波变换；

步骤六、根据小波变换系数的模提取小波脊；

步骤七、滤除噪声点；

步骤八、确定各个时刻的小波脊中心尺度；

步骤九、根据小波脊中心尺度计算出各个时刻水下声信号瞬时频率。

本发明的有益效果是：

1、本发明能够实现水下声信号瞬时频率的提取，即能够实现变频水下声信号的探测， 解决了常规不能解调瞬时频率的缺点。

2、实验表明，本发明瞬时频率的解调精度优于100Hz。

3、对于稳频的水下声信号频率的提取，本发明采用的方法比常规的频谱分析方法具 有更好的抗干扰性，本发明对稳频水下声信号的频率测量，精度可优于1Hz。

附图说明

图1为典型的迈克尔逊干涉系统基本光路结构图；

图2为具体实施方式一中，小波脊方法确定水下声信号瞬时频率的一般步骤示意图；

图3为本发明实验验证中，不同时刻的小波变换的脊特征示意图，其中(a)为t＝t₃₅₀₀时刻的小波变换的脊特征示意图，(b)为t＝t₃₆₀₀时刻的小波变换的脊特征示意图；

图4为本发明实验验证中，小波系数模值分布图；

图5为本发明实验验证中，小波脊分布图；

图6为本发明实验验证中，频率变化直线拟合结果示意图；

具体实施方式

结合图1、图2进一步详细说明本发明的具体实施方式。

具体实施方式一：本实施方式所述的一种基于小波脊的水下声信号瞬时频率解调方 法，包括以下步骤：

步骤一、利用迈克尔逊激光干涉系统探测水表面；

步骤二、采集激光干涉信号，所述采集过程为：干涉信号被光电探测器接收，并将其 转化为电信号，该电信号通过数据采集卡采集后，传送给上位机；

步骤三、对采集到的干涉信号进行频谱分析，用于确定水下声信号的频率分布范围；

步骤四、确定小波变换的观测尺度；

经过步骤三后，就可以知道信号谱的频率范围，根据这一范围结合母小波中心频率、 采样频率即可确定小波变换的观测尺度。

步骤五、对干涉信号进行小波变换；

步骤三所述的频谱分析是为了确定尺度范围，不是对信号本身处理的必要程序。这里 的干涉信号就是指数据采集卡采集到的信号。

步骤六、根据小波变换系数的模提取小波脊；

步骤七、滤除噪声点；

就是通过设定阈值的方法来滤除，就是对小波脊对应的小波系数的模的最小值有一个设 定，只有当小波脊对应的小波系数的模大于该值，该脊点才被认为是有效的。

步骤八、确定各个时刻的小波脊中心尺度；

步骤九、根据小波脊中心尺度计算出各个时刻水下声信号瞬时频率。

激光干涉法探测水下声信号介绍如下：

水下声源发出声信号会引起水表面横向传播的表面声波，通常采用零差干涉光路原理 对水表面进行探测，典型的零差干涉法光路如图1所示。

激光干涉探测水下声信号都通过探测水表面声波的方式来实现，对于零差干涉法探 测系统光路通过分光镜将激光束分为探测光和参考光，探测光入射到水表面，受到水表面 声波调制后返回与参考光束相遇发生干涉，干涉信号被光电探测器接收，由光电探测器接 收干涉信号，并将其转化为电信号，该电信号通过数据采集卡和采集程序可被上位机读取， 被上位机读取的干涉信号做进一步处理可解调出水下声源的频率，例如做FFT变换。

根据前述测量原理，参考光的振幅分布可由下式表示：

E_b(t)＝A_b sin(ω₀t+kz_b+φ)   (1)

式中，ω₀——参考光角频率，

k——波数，$k=\frac{2\pi }{\lambda };$

Ф——激光束的初始相位。

测量光的振幅分布可有下式表述：

E_m(t)＝A_msin(ω₀t+kz_m+φ)   (2)

两束光汇合后的合振幅由下式表述：

E(t)＝E_b(t)+E_m(t)＝A_m sin(ω₀t+kz_m+φ)+A_b sin(ω₀t+kz_b+φ)   (3)

因此干涉光强可由下式表述：

光强由光电探测器接收，探测器不对前3个高频项响应(即转变为直流信号)，因此， 去除直流分量接收到的光强信号可由下式表述：

I_d＝KA_mA_b cos[k(z_m-z_b)]＝A_d cos[k(z_m-z_b)]   (5)

式中，K——光电探测器放大倍率。

由上式可知，光电探测器接收的信号是关于参考光和测量光光程差的调制信号，两束 激光的光程差由两个臂的初始长度和水表面振动幅度有关，光程差z_m-z_b可由下式表示：

$z_m-z_b=2(L_1-L_2)+2\underset{n}{\Sigma }\left(A_n\sin (2\pi f_{n}t+{\phi }_n)\right)+{2A}_s\sin (2\pi f_{s}t+{\phi }_s)---\left(6\right)$

式中，为自然水表面的振动，A_ssin(2πf_st+φ_s)为水下声源引起的水 表面振动，(L₁-L₂)为水面静止时两路光的初始光程差，记为ΔL。因此光电探测器接收到 的光强信号可由式(7)表述，对该信号进行频谱分析即可得到水表面声波的频率，即水 下声信号的发声频率。

$I_d=A_d\cos \left\{k\right[2\mathrm{\Delta L}+2\underset{n}{\Sigma }\left(A_{\text{n}}\sin (2\pi f_{n}t+{\phi }_n)\right)+2A_s\sin (2\pi f_{s}t+{\phi }_s)\left]\right\}---\left(7\right)$

基于小波脊的水下声信号激光干涉探测方法基本原理如下：

针对频谱分析方法不能实时地确定水下声信号频率这一缺点，研究利用小波变换的方 法处理干涉信号。不失一般性地假定水表面的自然波动为A_nsin(2πf_nt+φ_n)，那么根据三角函 数的和差化积公式，以及贝塞尔恒等式可将式(7)进一 步分解为f_n和f_s的整数倍谐波项以及它们的和频项和差频项之和，如式(8)所示，式中， x₁＝2kA_n，x₂＝2kA_s。

$\left(\begin{array}{c}{I}_d=A_d\cos \left(2\mathrm{k\Delta L}\right)J_0\left(x_1\right)J_0\left(x_2\right)\\ -2A_d\sin \left(2\mathrm{k\Delta L}\right)J_0\left(x_2\right)\underset{m=0}{\overset{\infty }{\Sigma }}{J}_{2m+1}\left(x_1\right)\sin \left((2m+1)(2\pi f_{n}t+{\phi }_n)\right)\\ +2A_d\cos \left(2\mathrm{k\Delta L}\right)J_0\left(x_2\right)\underset{m=0}{\overset{\infty }{\Sigma }}{J}_{2m}\left(x_1\right)\cos \left(2m({2\mathrm{\pi f}}_{n}t+{\phi }_n)\right)\\ -2A_d\sin \left(2\mathrm{k\Delta L}\right)J_0\left(x_1\right)\underset{m=0}{\overset{\infty }{\Sigma }}{J}_{2m+1}\left(x_2\right)\sin \left((2m+1)({2\mathrm{\pi f}}_{s}t+{\phi }_s)\right)\\ +2A_d\cos \left(2\mathrm{k\Delta L}\right)J_0\left(x_1\right)\underset{m=1}{\overset{\infty }{\Sigma }}{J}_{2m}\left(x_2\right)\cos \left(2m(2\pi f_{s}t+{\phi }_s)\right)\\ -4A_d\cos \left(2\mathrm{k\Delta L}\right)\left[\underset{m=0}{\overset{\infty }{\Sigma }}{J}_{2m+1}\left(x_1\right)\sin \left((2m+1)({2\mathrm{\pi f}}_{n}t+{\phi }_n)\right)\right]\left[\underset{m=0}{\overset{\infty }{\Sigma }}{J}_{2m+1}\left(x_2\right)\sin \left((2m+1)(2\pi f_{s}t+{\phi }_s)\right)\right]\\ -4A_d\sin \left(2\mathrm{k\Delta L}\right)\left[\underset{m=1}{\overset{\infty }{\Sigma }}{J}_{2m}\left(x_1\right)\cos \left(2m(2\pi f_{n}t+{\phi }_n)\right)\right]\left[\underset{m=0}{\overset{\infty }{\Sigma }}{J}_{2m+1}\left(x_2\right)\sin \left((2m+1)(2\pi f_{s}t+{\phi }_s)\right)\right]\\ +4A_d\cos \left(2\mathrm{k\Delta L}\right)\left[\underset{m=1}{\overset{\infty }{\Sigma }}{J}_{2m}\left(x_1\right)\cos \left(2m(2\pi f_{n}t+{\phi }_n)\right)\right]\left[\underset{m=1}{\overset{\infty }{\Sigma }}{J}_{2m}\left(x_2\right)\cos \left(2m(2\pi f_{s}t+{\phi }_s)\right)\right]\\ -4A_d\sin \left(2\mathrm{k\Delta L}\right)\left[\underset{m=0}{\overset{\infty }{\Sigma }}{J}_{2m+1}\left(x_1\right)\sin \left((2m+1)(2\pi f_{n}t+{\phi }_n)\right)\right]\left[\underset{m=1}{\overset{\infty }{\Sigma }}{J}_{2m}\left(x_2\right)\cos \left(2m(2\pi f_{s}t+{\phi }_s)\right)\right]\end{array}\right)---\left(8\right)$

观察分解后的干涉信号表达式可知：由于自然水表面的波动为低频波动，干涉信号的 频谱分布出现了一个非常重要的特点，即在以f_s为中心的一个频带内出现较密集的谱线， 且谱线的分布呈现严格的对称性。各谱线对应的频率为f_s±2nf_n或f_s±(2n-1)f_n，n∈N⁺。 它们的幅值应分别正比于2A_dsin(2k△L)J_2n(x₁)J₁(x₂)和2A_dcos(2k△L)J_2n-1(x₁)J₁(x₂)，而f_s处谱线的幅值 则正比于2A_d sin(2k△L)J₀(x₁)J₁(x₂)。

根据上述特点可知，在该频带内频谱峰值的位置可能有两种情况，视贝塞尔函数值 J₀(x₁)以及初始相位差2kΔL的正余弦值而定，或位于f＝f_s处，或对称地分布在f＝f_s的左 右两侧，即使自然水表面的波动有一定的随机性，这样的频谱分布特点仍然贯穿整个时域 范围。因此，当我们对接收到的干涉信号进行小波变换，在特定的观察尺度内，小波脊点 总是沿着f_s所在尺度分布，当f_s是随时间变化的量，那么小波脊点则沿着f_s所对应的尺度 时间曲线来分布。因此根据这一结论我们可以通过计算小波脊点的位置来对水表面声波的 瞬时频率(即水下声信号的瞬时频率)进行判定。

具体实施方式二：本实施方式与具体实施方式一不同的是：步骤三所述的频谱分析， 是指对采集到的干涉信号采用快速傅里叶变换算法，即FFT算法。其它步骤与具体实施 方式一相同。

具体实施方式三：本实施方式与具体实施方式一或二不同的是：步骤五所述的小波变 换的具体过程为：

选择Morlet小波为母小波，其近似数学表达式为

$\psi \left(t\right)=\exp (-\frac{t^2}{2\sigma })\exp \left(j{\omega }_{0}t\right)---\left(\text{10}\right)$

式中，ω₀为母小波的中心角频率，σ为高斯函数参数，j是虚数符号，t为时间，其相应 的傅里叶变换为：

$\psi \left(\omega \right)=\sqrt{2\mathrm{\pi \sigma }}\exp \left[\frac{-\sigma {(\omega -{\omega }_0)}^2}{2}\right]---\left(11\right)$

ω是角频率，则对干涉信号的小波变换系数计算公式如下：

$W_f(a,b)=\frac{1}{\sqrt{\alpha }}\int f\left(t\right)\psi \left(\frac{t-b}{a}\right)\mathrm{dt}---\left(12\right)$

式中，是以Morlet小波为母小波的一个函数族，a为尺度因子，b为平移因 子，即是通过对小波母函数的平移和伸缩得到的一个函数族，f(t)为采集到的干涉 信号，W_f(a,b)为小波变换系数。

其它步骤与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一至三之一不同的是：步骤六所述的根 据小波变换系数的模提取小波脊的具体过程为：

根据小波变换系数计算公式得到小波变换系数W_f(a,b)，对小波变换系数取模值， 然后利用局部模极大值算法来计算局部模极大值，局部模极大值点即为小波脊点，得到小 波脊点的分布位置(t₀,a₀)，这里的t₀就是小波脊点所在时刻，a₀就是小波脊点所在的尺 度。

实际小波脊点分布位置的计算过程如下：完成小波变换后得到的是一系列小波系数， 如果时间长度为m，计算的尺度长度为n，那么小波系数就是m*n的矩阵，我们想知道 t₀时刻的小波脊点分布状态，就把t＝t₀这一时刻的小波系数拿出来，这时候拿出来的小波 系数就是一个(1*n)的序列，计算这一序列的极大值出现哪个位置，也就是说计算序列 的第几个值是极大值，这个极大值对应的尺度值就是小波脊点分布位置的尺度值a₀。其 它步骤与具体实施方式一至三之一相同。

具体实施方式五：本实施方式与具体实施方式一至四之一不同的是：步骤八所述的确 定各个时刻的小波脊中心尺度的具体过程为：

小波脊点的位置或位于水下声信号频率对应的尺度线上，或对称分列在尺度线的上下。 因此，对于只有一个脊点的时刻，该脊点所在尺度即为小波脊中心尺度；

对于有两个脊点的时刻，则这两个脊点的中点所在尺度为小波脊中心尺度，小波脊中 心尺度对应的频率为该时刻的瞬时频率。

小波脊中心尺度就是指各小波脊点中心所对应的尺度，下文有解释，不同的时刻有两 种小波脊状态，一种是只有一个脊点，一种是两个脊点，对于只有一个脊点的情况，小波 脊中心尺度就是该脊点对应的尺度，对于两个脊点的情况，小波脊中心尺度就是这两个小 波脊点的中点所对应的尺度。其它步骤与具体实施方式一至四之一相同。

具体实施方式六：本实施方式与具体实施方式一至五之一不同的是：步骤九所述的根 据小波脊中心尺度计算出各个时刻水下声信号瞬时频率的具体过程为：

假定t时刻干涉信号的小波脊中心尺度为a_c，那么水下声信号t时刻的瞬时频率f计 算公式如式(13)所示，式中f_s为干涉信号采集模块的采样率

f＝f_a·f_s/a_c   (13)

其中，f_a是Morlet小波的中心频率，是一个设定值。其它步骤与具体实施方式一至五之 一相同。

本发明的实验验证如下：

为观察小波脊点的分布特点，利用前述零差干涉法测量系统对频率为2kHz的水表面 声波进行激光干涉探测，得到了受水表面声波调制的干涉信号。以时频域内都有良好聚焦 性的解析小波作为小波基对该干涉信号进行小波变换，得到小波系数矩阵，根据设定的小 波中心频率和信号采样的时间间隔计算出水下声信号频率2kHz对应的尺度a_s＝10.40。图 3所示为两个不同时刻的小波系数模值随尺度变化的分布情况，由图可知，水表面声波激 光干涉探测信号的小波变换存在两种截然不同的脊特征：小波脊点或位于水下声信号频率 对应的尺度上，或分列在该尺度的两侧，与前文理论分析一致。

为验证小波脊方法提取水下声信号瞬时频率的准确性，对频率呈线性变化的水下声信 号进行了探测，设定水下声信号的频率变化率为8kHz/s，对该水表面声波进行激光干涉 探测，得到干涉信号进行小波变换后系数模值的分布如图4所示，由图可知小波系数模值 的局部极大值所在尺度随着时间的变化而变化，也即小波脊的分布位置随着水表面声波频 率的变化而变化。小波脊点分布如图5所示，小波脊点的中心位置由尺度为7逐渐降低到 尺度为5，这表明小波脊能够实时跟踪水表面声波频率的变化，小波脊所在尺度能够反映 水下声信号的瞬时频率。

根据瞬时频率的确定方法，计算出观测时间内各个时刻的瞬时频率，然后拟合频率时 间直线，计算出频率变化率是否与实际值相吻合。拟合结果如图6所示，图中直线即为频 率时间变化曲线，计算出来的频率变化率(即斜率)为8.12kHz/s，与理论值基本相吻合， 其误差是由于频率分辨率及小波脊提取时的噪声点等因素造成的，这证明了小波脊方法提 取水下声信号瞬时频率的准确性。

本发明的技术方案主线为：1)对信号小波变换得到小波变换系数；2)根据局部模极 大值方法得到各个时刻小波脊点的分布位置；3)根据小波脊点的分布位置计算出各个时 刻的小波脊中心尺度；4)根据小波中心尺度反求各个时刻的瞬时频率。本发明适用于水 下发声目标发声频率的探测、液体表面声波频率探测、微小振幅机械振动频率的测量。