0 引　言

在海洋环境中，由于传播信道中海洋介质随机散射及信道的多途、频散等水声环境效应，使得水中目标辐射噪声信号波形发生畸变，从而改变接收端信号的频谱特征。为研究水中目标辐射噪声通过水下声信道传播后噪声特征的变化规律，需要分析表征目标辐射噪声典型特征参数的变化。描述目标辐射噪声的特征参数包括目标辐射噪声的基频、谐波振幅及各阶谐波的信噪比^[1]。水下目标辐射噪声中的谐波分量包含了反映目标自身本质特性的信息，能否有效提取目标谐波特征关系到目标类型识别的效果。

为此，本文提出了一种水声目标辐射噪声谐波特征提取算法，结合仿真信号与实测数据进行对比，验证了该算法的可行性。

1 水声目标辐射噪声信号模型

目标辐射噪声由宽带噪声与多个正弦信号组成，可表示为^[2]：

$\begin{split}r(t) = & s(t) + n(t)=\\& \sum\limits_{h = 1}^H {{A_h}} \cos (2\pi h\gamma t + {\varphi _h}) + n(t)\end{split}\text{。} $

(1)

其中，H为谐波阶数，A_h和φ_h为第h阶谐波分量的振幅和相位，γ为基频，与目标的引擎转速等参数有关，因此可将r（t）表示为 $r(t,\theta )$ ，其中θ表示待估计参数集合，包括γ及A_h，即 $\theta = \left\{ {\gamma ,{A_h}} \right\}$ 。

由于引擎速度的变化，多普勒频移以及传播环境的影响，被动声呐接收到的信号与目标实际辐射的噪声信号不同。由于基频时变，将γ表示为 $\gamma (t) = {f_o} + \Delta f(t)$ ，其中f_o为基频，∆f（t）为基频随时间的变化量，因此被动声呐的接收信号 $r(t,\theta )$ 可表示为：

$\begin{split}r(t,\theta ) = & s(t,\theta ) + n(t)=\\& \sum\limits_{h = 1}^H {{A_h}} \cos (2\pi h\gamma (t)t + {\varphi _h}) + n(t)\end{split}\text{。}$

(2)

被动声呐接收信号是频率时变的连续波信号，在信号处理前须对其进行采样。利用短时傅里叶变换（STFT）按一定窗口长度滑动，并计算每个窗内数据的傅里叶变换，得到信号频率随时间变化关系^[3]。假设时间窗的长度足够小，则窗内数据对应的频率变化可近似忽略，即 $\Delta f(t) \cong \Delta f({t_k})$ ，则 ${\gamma _k} = {f_o} + \Delta f({t_k})$ ，其中t_k为第k段数据快拍的中心时刻。每个窗内数据对应的STFT是信号数据傅里叶变换与时域窗傅里叶变换的卷积，即 $S(f,\theta ) = F\left\{ {s(t,\theta )} \right\}*F\left\{ {{\rm{rectwin}}(t)} \right\}$ $S(f,\theta )$ 也可表示为辛格函数的加权和，即

$S(f,{\theta _k}) = \sum\limits_{h = 1}^H {\frac{{{A_h}}}{2}} \sin {\mathop c\nolimits} [\pi (f - h{\gamma _k})]\text{，}$

(3)

假设 $r(t,\theta )$ 时长为T_r，采样频率为f_s，采样周期 $\Delta t = 1/{f_s}$ ，因而 $r(t,\theta )$ 可表示为 $r(j\Delta t,\theta )$ ，j = {0, 1, …, N_r-1}，且 ${N_r} = {T_r}/\Delta t$ 。信号采样可分为K个重叠数据快拍，每个快拍长度为T，或 ${N_s} = T/\Delta t$ 个采样点。利用 ${r_k}(n\Delta t,{\theta _k})$ 表示第K个快拍的第n个采样点，而 $r(n\Delta t,\theta )$ 的维数为 ${N_S} \times K$ 。对每个快拍计算N_S点FFT，即

${R_{ko}}(m\Delta f,{\theta _k}) \!=\! \frac{2}{{{N_S}}}\sum\limits_{n = - {N_S}/2}^{{N_S}/2 - 1} {{r_k}} ({t_k} \!\!+\! n\Delta t,{\theta _k}){e^{ - jm2\pi \Delta fn\Delta t}}\text{。}\!\!$

(4)

其中 ${R_{ko}}(m\Delta f,{\theta _k})$ 是频率分量 $m\Delta f$ 的傅里叶变换系数，频率分辨率为 $\Delta f = 1/T$ 。记 ${f_m} = m\Delta f$ ，其中f_m= ∆f {−N_s/2, ···, 0, ···, N_s/2−1}，即 ${R_{{\rm{ko}}}}(m\Delta f,{\theta _k}) = {R_{{\rm{ko}}}}({f_m},{\theta _k})$ ，维数为 ${N_S} \times K$ 。

2 水声目标辐射噪声谐波特征提取方法

水下目标谐波特征提取算法的总体思路。主要分为 3 个步骤：

步骤 1　计算信号频谱与梳状滤波器的频率响应间的互相关，完成对基频的估计；

步骤 2　利用卡尔曼滤波器跟踪基频量测随时间的变化，得到基频的最优估计；

步骤 3　利用基频的最优估计，提取相应谐波的振幅信息。

2.1 基频估计方法

由式（3）可知，目标辐射噪声的频谱可表示为一组在t_k时刻基频为γ_k的sinc函数的加权和，利用该模型构建梳状滤波器的频率响应为：

$\begin{split}{C_{mp}} =& S({f_m},{\theta _k}|{\gamma _k} = {\zeta _p},{A_h} = 2)=\\& \sum\limits_{h = 1}^H {\sin {\mathop c\nolimits} \left[ {\pi ({f_m} - h{\zeta _p})} \right]}\text{。}\end{split}$

(5)

在此假设所有谐波的振幅相等，而C_mp是矩阵 ${\bar { C}}$ 的一行， ${\bar { C}}$ 的维数为 ${N_\gamma } \times {N_S}$ 。在此基频量化为间隔为∆γ的N_γ个频点，ζ表示每个梳状滤波器频率响应的基频，ζ_p为 ${\zeta _p} = {\gamma _{\min }} + p\Delta \gamma $ 。

对于均值分别为 ${\mu _x}\left( t \right)$ 和 ${\mu _y}\left( t \right)$ 的随机过程X（t）和Y（t），其互相关函数可以表示为：

$\rho \left( t \right) = \frac{{E\left[ {\left( {X\left( t \right) - {\mu _x}\left( t \right)} \right)\left( {Y\left( t \right) - {\mu _y}\left( t \right)} \right)} \right]}}{{\sqrt {E\left[ {{{\left( {X\left( t \right) - {\mu _x}\left( t \right)} \right)}^2}} \right]E\left[ {{{\left( {Y\left( t \right) - {\mu _y}\left( t \right)} \right)}^2}} \right]} }}\text{，}$

(6)

ρ（t）取值在[–1，1]间。要估计t_k时刻接收信号 ${R_k}({f_m},{\theta _k})$ 的基频，则需计算样本信号模型与接收信号间的互相关函数^[4]：

${\hat \rho _{pk}}({\zeta _p},{t_k}) = \frac{{\sum\nolimits_{fm} {({C_{mp}} - {\mu _C}){R_k}({f_m},{\theta _k})} }}{{\sqrt {\sum\nolimits_{fm} {{{({C_{mp}} - {\mu _C})}^2}} } \sqrt {\sum\nolimits_{fm} {{R_k}{{({f_m},{\theta _k})}^2}} } }}\text{。}$

(7)

其中 ${\hat \rho _{pk}}({\zeta _p},{t_k})$ 表示第k个快拍 ${R_k}({f_m},{\theta _k})$ 与基频为ζ_p的梳状滤波器响应间的相关系数，简化记为 ${\hat \rho _{pk}} \equiv {\hat \rho _{pk}}({\zeta _p},{t_k})$ 。如此计算得到每个快拍与相应梳状滤波器间的相关系数，则得到维数为 ${N_\gamma } \times K$ 的矩阵 ${ {\bar \rho}} $ 。当样本信号模型C_mp的基频ζ_p与接收信号 ${R_k}({f_m},{\theta _k})$ 的基频γ_k相等时，样本信号模型与接收信号高度相关，则在 $\bar \rho $ 中将出现峰值，而该峰值对应频率为 ${\hat \gamma _k}$ 。

2.2 基频最优估计 2.2.1 基于卡尔曼滤波理论的基频跟踪

利用式（7）计算 ${R_k}({f_m},{\theta _k})$ 与所有可能样本信号C_mp的相关系数。基频估计就是找k时刻对应的ζ_p，即 ${\hat \rho _{pk}}$ 最大时，得 ${\hat \gamma _k} = \arg {\max _p}{\hat \rho _{pk}}$ 。由于目标辐射噪声可能存在多个不同的频率，如引擎的齿轮齿数比不是整数时，引擎谐波对应的基频不同于轴和螺旋桨产生的谐波对应的基频，所以直接取相关系数最大值对应的频率作为基频的估计不够准确。考虑到基频的时变性，借鉴卡尔曼滤波理论，对每个峰值对应的时变基频进行估计。 $\bar \rho $ 中每个峰值及其对应的频率 ${\hat \gamma _k}$ 将组合成量测集Z_k，作为卡尔曼滤波的量测。得到k–1时刻的状态估计 ${{\hat{ \varGamma }}_{k - 1}}$ 及其误差协方差估计 ${P_{k - 1}}$ 后，通过式（8）得到k时刻的状态预测值 ${{\hat{ \varGamma }}_k}$ 及其协方差矩阵 ${ P}_k^ - $ ，进一步预测^[5]：

${\hat{ \varGamma }}_k^ - = {{F}}{{\hat{ \varGamma }}_{k - 1}} \text{，}\;\;\;{{P}}_k^ - = {{F}}{{{P}}_{k - 1}}{{F}^{\rm T}} + {{Q}} \text{。}$

(8)

其中F是状态转移矩阵，Q是过程噪声协方差。利用k时刻量测z_k，通过式（9）进行状态更新得到k时刻的状态估计 ${{\hat{ \varGamma }}_k}$ 及其误差协方差 ${{P}}_k^{}$ ，即状态更新：

$\left\{ \begin{aligned}& {K_k} = {{P}}_k^ - {{{H}}^{\rm T}}{({{HP}}_k^ - {{{H}}^{\rm T}} + {{R}})^{ - 1}}\\& {{{\hat{ \varGamma }}}_k} = {\hat{ \varGamma }}_k^ - + {K_k}({{\bf{z}}_k} - {{{\rm H}\hat \Gamma }}_k^ - )\\& {{{P}}_k} = ({{I}} - {K_k}{{H}}){{P}}_k^ - \end{aligned} \right.$

(9)

其中K_t是k时刻的卡尔曼滤波增益，R是量测噪声协方差矩阵，H为量测矩阵。在下一时刻，将当前时刻的状态估计值反馈到预测方程中，并且重复式（8）和（9）的预测更新过程，最终得到所有时刻的均值和误差协方差的滤波估计，形成频率跟踪航迹。

2.2.2 基频跟踪算法

借鉴近似的匀速运动模型，基频频率变化跟踪，即

$\left\{ {\begin{array}{*{20}{l}}{{{{\varGamma }}_k} = {{F}}{{\bf{\varGamma }}_{k - 1}} + {{{w}}_k}}\text{，}\\{{{{\varGamma }}_k} = {{[{\gamma _{KF,k}},{{\dot \gamma }_{KF,k}},{\rho _{KF,pk}}]}^{\rm T}}}\text{，}\\{{{F}} = \left[ {\begin{array}{*{20}{c}}1 & 1 & 0\\0 & 1 & 0\\0 & 0 & 1\end{array}} \right]}\text{，}\\{{{{w}}_k} \sim N({{0}},{{Q}})}\text{。}\end{array}} \right.$

(10)

其中 ${{{\varGamma }}_k}$ 为单个频率航迹的状态向量， ${\gamma _{KF,k}}$ 和 ${\dot \gamma _{KF,k}}$ 分别为频率和k–1时刻到k时刻的频率变化量， ${\rho _{KF,pk}}$ 为频率 ${\gamma _{KF}}_{,k}$ 的相关系数，F为状态转移矩阵， ${{{w}}_k}$ 为均值为0，协方差矩阵为Q的高斯白噪声^[6]。目标位置的预测可通过前一时刻的位置加上速度与时间步长的乘积来得到，由于跟踪算法采用匀速运动，则k–1时刻到k时刻的频率变化可以忽略，那么k时刻的频率预测为 ${\gamma _{KF,k - 1}} + {\dot \gamma _{KF,k}}$ 。而实际频率的变化可用过程噪声项来描述，将k时刻的状态及量测z_k的关系表示为：

$\left\{ {\begin{array}{*{20}{l}}{{{{z}}_k} = {{H}}{{{\varGamma }}_k} + {v_k}}\text{，}\\{{{{z}}_k} = {{[{{\hat \gamma }_k},{{\hat \rho }_{pk}}]}^{\rm T}}}\text{，}\\{{{H}} = \left[ {\begin{array}{*{20}{c}}1 & 0 & 0\\0 & 0 & 1\end{array}} \right]}\text{，}\\{{v_k} \sim N(0,R)}\text{。}\end{array}} \right.$

(11)

其中z_k为频率 ${\hat \gamma _k}$ 和相关值 ${\hat \rho _{pk}}$ 的测量值，H为描述状态与量测关系的测量矩阵。假设v_k是均值为零，方差为R的高斯白噪声，过程噪声协方差矩阵Q和量测噪声协方差矩阵R为：

${{Q}} = \left[ {\begin{array}{*{20}{c}}{{Q_\gamma }} & 0 & 0\\0 & {{Q_{\dot \gamma }}} & 0\\0 & 0 & {{Q_\rho }}\end{array}} \right]\text{，}\;\;\;{{R}} = \left[ {\begin{array}{*{20}{c}}{{R_\gamma }} & 0\\0 & {{R_\rho }}\end{array}} \right]\text{。}$

(12)

实际上，过程噪声或量测噪声的大小蕴含了利用该模型进行状态估计时的可信度。

2.2.3 基频跟踪逻辑

卡尔曼滤波需要一个初始的状态估计，跟踪过程需要航迹起始和结束的逻辑判断^[7]。首先设置相关系数门限λ_γ，对每个快拍对应的相关系数 $\bar \rho $ 进行判断，可能得到若干个基频和数据吻合。集合Z_k中的基频就是k时刻跟踪器的量测，利用这些检测得到的频率用于航迹起始。航迹可能会出现3种状态：1）起始：如果在连续N个快拍内检测到基频M次，则认为航迹已建立；2）标记：如果在量测Z_k中，没有与航迹关联的观测，则对该航迹做结束标记；3）结束：如果连续N_F个快拍都做航迹结束标记，则结束该条航迹。每一步都需要将观测集与航迹进行关联，若观测量与某一航迹满足：

$({z_k} - H\hat \Gamma _k^ - ){(HP_k^ - {H^{\rm T}} + R)^{ - 1}}{({z_k} - H\hat \Gamma _k^ - )^{\rm T}} < {\chi ^2}\text{。}$

(13)

则该观测量与航迹相关，χ²度量的是每个测量频率和卡尔曼滤波预测得到频率的差，χ²值越小，表明测量频率和预测频率越接近。当有多个观测频率满足上述关系式时，最小χ²值对应的频率为预测频率。

2.3 谐波特征提取

实际中，卡尔曼滤波输出可能有多条航迹，每条航迹是 ${R_k}({f_m},{\theta _k})$ 中谐波分量对应基频的时变估计。为得到与实际噪声拟合最好的航迹，利用相关系数计算参数 ${\bf{\psi }}$ ，得到：

${\bf{\psi }}\left[ k \right] = \sqrt {\frac{1}{L}\sum\limits_{k = 1}^L {{{\left| {{\rho _{KF,pk}}\left[ k \right]} \right|}^2}} } \text{。}$

(14)

其中k为航迹号， ${\rho _{KF,pk}}\left[ k \right]$ 为k航迹的相关系数经卡尔曼滤波后得到的估计值，L为快拍中航迹的长度。 ${\bf{\psi }}$ 值最大的航迹与实际噪声数据拟合得最好，对应基频的最优估计。 ${\bf{\psi }}$ 值为接收信号数据和样本模型间拟合好坏的度量。 ${\bf{\psi }}$ 值等于1表示接收信号数据与样本模型完全匹配。 ${\bf{\psi }}$ 值同样可度量信号中谐波分量的多少，因为参数 ${\bf{\psi }}\left[ k \right]$ 是相关系数估计值 ${\rho _{KF,pk}}$ 在整条航迹持续时间内的均值， ${\bf{\psi }}\left[ k \right]$ 越大，在谱图上可见的谐波分量越多。相反， ${\bf{\psi }}\left[ k \right]$ 越低，在谱图上可见的谐波分量就越少。

谐波特征是目标的“声音指纹”。由于航迹 ${\gamma _{KF,k}}$ 与数据拟合的性能可用 ${\bf{\psi }}$ 来衡量，为获得谐波特征，在所有谐波分量中，首先把拟合最好的频率航迹投影到声谱图上。对于每个快拍，航迹在声谱图上的投影就是找出最接近 ${\gamma _{KF,k}}$ 整数倍的频率单元f_m。在投影频率附近的小范围内搜索 ${R_{ko}}({f_m},{\theta _k})$ 的峰值，得到各阶谐波的振幅 ${\hat A_h}$ 。通过平均峰值两侧窗内的谱图估计得到各阶谐波的背景噪声，而窗的大小与目标运动的快慢和给定快拍内发生变化的谐波分量数密切相关。对整条航迹的振幅和背景噪声进行平均，得到谐波谱图。实际中，目标很难保持匀速运动，因此对于任意给定的快拍，谐波频率都可能扩展到多个频率单元。为保证估计性能，需要加大峰值搜索范围。

3 仿真结果与分析

按照式（1）模拟生成被动声呐的接收信号，谐波数为5，宽带噪声的功率为5。假设所有谐波振幅相等均为1，且信号时长有限，而噪声与信号时长相等。信号采样频率为8 kHz，采样点数为65 536，信号基频在[5，65] Hz间服从均匀分布。仿真参数见表1。

图1给出的是一个段时长8.192 s的目标噪声对应的声谱图，为便于显示，X轴截选显示500~900 Hz范围内频率。图2给出了各阶谐波与背景噪声的相对振幅特性，每根谱线给出的是各阶谐波的相对振幅，每个峰值周围的背景噪声用实线表示，图中所有谐波分量全部可见。图1中基频在整个航迹跟踪范围内有一些微小变化。在频率[500，550] Hz范围内，基频跟踪在谱图上的投影明显更亮一些，说明在这个范围内 ${\bf{\psi }}$ 值更高，对应的仿真信号与样本信号模型拟合得更好；而在频率范围[550，900] Hz内，基频跟踪在谱图上的投影亮度稍微暗一些，说明在这一范围内的 ${\bf{\psi }}$ 值低，拟合度稍差。从图2提取的各阶谐波振幅也可以说明这一点，图中有5个谐波幅度明显高于后面谐波的幅度，且谐波振幅在[500，550] Hz范围内附近有抖动，在[550，900] Hz范围内的变化较小，这与模拟仿真生成的谐波数对应。仿真结果表明，论文提出的算法能够提取得到目标辐射噪声的基频和各阶谐波的特征。

4 结　语

谐波特征是目标的“声音指纹”。每个目标自身辐射噪声对应的“声音指纹”都是由其引擎转速、汽缸数、齿轮齿数比、叶片数量等参数决定。本文提出了一种基于最大似然估计和卡尔曼滤波理论的目标辐射噪声谐波特征的提取与分析算法，经验证该算法的稳健性较好，且提取的谐波信息包含了反映目标自身本质特性的信息，可进一步用于不同类型的水下目标识别。