水下目标识别技术是水声信息处理与现代声呐系统的关键技术。最早提出和使用的水下目标特征是时域波形特征^[1]，如波长差分布特征和过零点分布特征。快速傅里叶变换算法提出后，基于功率谱估计的特征提取方法^[2]得到了快速发展，出现了基于现代谱估计的特征提取方法，以及能够表征水声信号非高斯分布特点的高阶谱特征^[3]。为了反映出非平稳信号的时变信息，出现了基于短时傅里叶变换、小波变换和希尔伯特黄变换的时频域特征提取方法^[4-6]。基于短时傅里叶变换的LOFAR谱分析是目前人工分析时广泛采用的技术手段。在非线性分析领域，研究者利用水下目标噪声信号中存在混沌现象，提出了混沌特征^[7]。随着人耳听觉模型研究的深入，模拟人耳听觉特性的听觉特征^[8]被应用于水下目标识别，其中梅尔频率倒谱系数(Mel-frequency cepstral coefficients，MFCC)特征^[9]在水下目标识别中有较多的应用。

水下目标辐射噪声包括机械噪声、螺旋桨噪声和水动力噪声^[10]。在频谱上看，主要包括宽带连续谱和窄带线谱。水下目标辐射噪声中的强能量窄带线谱成分往往包含水下目标的重要特征信息。相比于连续谱，目标信号的线谱成分往往更加稳定，而提取的特征中稳定的信息在识别中起到的作用更大^[11]，若能更好地反映线谱成分，可以有效地提高特征的噪声鲁棒性。而强能量窄带线谱成分在频谱中的分布是具有稀疏性的。

近年来，结构化稀疏信号的恢复重构问题一直是热门研究方向，其中一个具有代表性的工作是Zhang等^[12]提出的块结构化稀疏贝叶斯学习框架。通过对水下目标辐射噪声的短时傅里叶分析，时间连续的若干帧信号中强能量窄带线谱在各自频谱中出现的位置往往相同或相近。而稀疏的强能量窄带线谱成分的这种时间相关性正是一种时间相关的结构化稀疏特性。本研究受到这一特点的启发，将信号基于傅里叶字典分解，借助稀疏贝叶斯学习框架和结构化稀疏恢复算法，对分解系数中稀疏分布的大元素进行一定程度地增强和恢复，并将分解系数用于特征提取，以提高特征的噪声鲁棒性。

1 水下目标信号的稀疏性分析

由于水下目标辐射噪声本身包含宽频带连续谱成分，导致水下目标信号能量在频谱上的分布范围很广，因此频谱不具有稀疏性。但结合水下目标辐射噪声的特点可知，水下目标信号中的频谱中往往存在能量较强的窄带线谱成分，这些成分是区分水下目标的重要信息，而且它们在频谱中的分布具有稀疏性。

1.1 信号的短时傅里叶时频分析

水声信号可以被视为一种具有短时平稳特点的随机信号，短时傅里叶分析是一种基础的时频分析手段，能够获取信号的动态特征信息。通过对信号进行短时傅里叶变换，将获得的频谱映射到以时间为横轴、频率为纵轴，像素点颜色代表频谱幅值相对大小(图 1中由蓝到黄对应幅值由小到大)的二维图被称为短时傅里叶时频分析图。

图 1是数据集中A类目标信号的短时傅里叶谱图。从图中可以发现，在50~200 Hz范围内，存在稀疏分布的呈块状聚集的强能量窄带线谱成分，且其出现位置具有时间连续性和相关性。数据集的另外三类目标也各自在一定的频段内存在类似的呈横向条纹状分布的强能量窄带线谱成分。

1.2 信号基于离散傅里叶字典的分解

离散傅里叶变换是信号的一种线性变换，因而可以写成信号基于离散傅里叶字典的分解形式：

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{Dx}} $

(1)

式中：y是时域信号，分解系数x是相应的频域变换向量。D是离散傅里叶字典，它的每一列都是一个傅里叶变换基，其形式为

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{D}} = }\\ {\frac{1}{{\sqrt N }}\left[ {\begin{array}{*{20}{c}} 1&1&1&1& \cdots &1\\ 1&{W_N^1}&{W_N^2}&{W_N^3}& \cdots &{W_N^{N - 1}}\\ 1&{W_N^2}&{W_N^4}&{W_N^6}& \cdots &{W_N^{2\left( {N - 1} \right)}}\\ 1&{W_N^3}&{W_N^6}&{W_N^9}& \cdots &{W_N^{3\left( {N - 1} \right)}}\\ \vdots&\vdots&\vdots&\vdots&\ddots&\vdots \\ 1&{W_N^{N - 1}}&{W_N^{2\left( {N - 1} \right)}}&{W_N^{3\left( {N - 1} \right)}}& \cdots &{W_N^{\left( {N - 1} \right)\left( {N - 1} \right)}} \end{array}} \right]} \end{array} $

(2)

式中: $W_N^k = {\rm{exp}}({\rm{j2 \mathsf{ π} }}k/N) $，N为信号长度。

从信号基于傅里叶字典分解的角度来看，包含强能量窄带线谱的信号分解系数x中的大值元素是稀疏分布的，而且往往呈块状聚集出现，这被称为块稀疏结构。除此之外，通过短时傅里叶变换，相邻时间帧信号的分解系数x中大值元素出现的位置具有时间相关性，这是一种时间相关的结构化稀疏特性。

考虑含有噪声v的情形，信号基于傅里叶字典的分解问题可以表示为

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{Dx}} + \mathit{\boldsymbol{v}} $

(3)

为了利用分解系数x中大值元素具有时间相关的块结构化稀疏特点，本文采用层次贝叶斯模型对式(3)的信号分解问题进行建模。

2 多帧联合结构化稀疏特征的提取

在稀疏贝叶斯模型中，式(3)的分解系数x是随机变量。引入先验分布后，利用贝叶斯变分算法对模型中的各个隐变量进行推断，最终得到分解系数x的后验估计。通过合理设置超参数的共享方式，通过超参数的不断迭代更新过程，就能够在估计稀疏分解系数x的同时，考虑x中大元素的结构化稀疏特点，增强结构信息在分解系数x估计结果中的表达。

2.1 稀疏贝叶斯学习模型

借助稀疏贝叶斯学习模型^[13]，模型假设N维时域信号y基于离散傅里叶字典D_N×L的L维分解系数x是服从多元高斯分布的随机变量。x的各个元素x_i的先验分布是均值为0，精确度(方差倒数)为α_i的高斯分布，且各个元素x_i之间独立，则x的概率分布密度为

$ p\left( {\mathit{\boldsymbol{x}}\left| \mathit{\boldsymbol{\alpha }} \right.} \right) = \sum\limits_{i = 1}^L {N\left( {{\mathit{\boldsymbol{x}}_i}\left| 0 \right.,\mathit{\boldsymbol{\alpha }}_i^{ - 1}} \right)} $

(4)

其中α={α₁, α₂, …, α_N}的各个分量独立，且服从Gamma分布。其概率密度函数为

$ p\left( {\mathit{\boldsymbol{\alpha }}\left| {\mathit{\boldsymbol{c}},\mathit{\boldsymbol{d}}} \right.} \right) = \prod\limits_{i = 1}^L {{\rm{Gamma}}\left( {{\mathit{\boldsymbol{\alpha }}_i}\left| {{\mathit{\boldsymbol{c}}_i},{\mathit{\boldsymbol{d}}_i}} \right.} \right)} $

(5)

噪声v各分量服从均值为0，精确度为α₀的高斯分布。并假设α₀是服从Gamma分布的随机变量。其概率密度函数为

$ p\left( {{\alpha _0}\left| {a,b} \right.} \right) \sim {\rm{Gamma}}\left( {{\alpha _0}\left| {a,b} \right.} \right) $

(6)

在此基础上，式(3)所表示的含噪结构化稀疏信号的分解问题可以表示为以下条件概率模型：

$ p\left( {\mathit{\boldsymbol{y}}\left| {\mathit{\boldsymbol{x}};{\alpha _0}} \right.} \right) \sim {N_{y\left| x \right.}}\left( {\mathit{\boldsymbol{Dx}},\alpha _0^{ - 1}\mathit{\boldsymbol{I}}} \right) $

(7)

根据高斯分布的条件分布性质可得x的后验分布仍然是高斯分布，其均值μ和协方差Σ分别为

$ \mathit{\boldsymbol{\mu }} = {\alpha _0}\mathit{\Sigma }{\mathit{\boldsymbol{D}}^{\rm{T}}}\mathit{\boldsymbol{y}} $

(8)

$ \mathit{\Sigma } = {\left( {{\alpha _0}{\mathit{\boldsymbol{D}}^{\rm{T}}}\mathit{\boldsymbol{D}} + \mathit{\boldsymbol{ \boldsymbol{\varLambda} }}} \right)^{ - 1}} $

(9)

式中：Λ=diag(α)。

通过积分将隐随机变量α消去，x的边缘分布p(x|c, d)是服从t分布的。t分布在0值处有一个尖峰，因而可以用来增加稀疏性。

为了可以有效利用各次观测向量间存在的时间相关性信息来提高混有噪声的结构化稀疏信号的恢复性能，借助稀疏贝叶斯多任务学习模型^[14]，对相邻M帧信号{y⁽¹⁾, y⁽²⁾, …, y^(M)}的分解问题进行联合求解。{y⁽¹⁾, y⁽²⁾, …, y^(M)}基于同一离散傅里叶字典D的分解系数为{x⁽¹⁾, x⁽²⁾, …, x^(M)}。模型假设噪声n⁽ⁱ⁾=y⁽ⁱ⁾-Dx⁽ⁱ⁾, i=1, 2, …, M都服从均值为0，精确度为α₀的高斯分布。将α⁽ⁱ⁾, i=1, 2, …, M按列排列组成矩阵A，即

$ \mathit{\boldsymbol{A}} = \frac{1}{M}\left( {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{\alpha }}^{{x_1}}}}&{{\mathit{\boldsymbol{\alpha }}^{{x_2}}}}& \cdots &{{\mathit{\boldsymbol{\alpha }}^{{x_M}}}} \end{array}} \right) $

(10)

这样x⁽ⁱ⁾的高斯后验分布的均值μ_xⁱ和方差Σ_xⁱ分别为

$ {\mathit{\Sigma }_{{x^i}}} = {\mathit{\Sigma }_0} - {\mathit{\Sigma }_0}{\mathit{\boldsymbol{D}}^{\rm{T}}}{\left( {{\alpha _0}I + \mathit{\boldsymbol{D}}{\mathit{\Sigma }_0}{\mathit{\boldsymbol{D}}^{\rm{T}}}} \right)^{ - 1}}\mathit{\boldsymbol{D}}{\mathit{\Sigma }_0} $

(11)

$ {\mathit{\boldsymbol{\mu }}_{{x^i}}} = {\mathit{\Sigma }_0}{\mathit{\boldsymbol{D}}^{\rm{T}}}{\left( {{\alpha _0}I + \mathit{\boldsymbol{D}}{\mathit{\Sigma }_0}{\mathit{\boldsymbol{D}}^{\rm{T}}}} \right)^{ - 1}}{\mathit{\boldsymbol{y}}_i} $

(12)

其中Σ₀满足：

$ \mathit{\Sigma }_0^{ - 1} = {\rm{diag}}\left( {\begin{array}{*{20}{c}} {\frac{1}{M}\sum\limits_{j = 1}^M {{\mathit{\boldsymbol{A}}_{1j}}} }&{\frac{1}{M}\sum\limits_{j = 1}^M {{\mathit{\boldsymbol{A}}_{2j}}} }& \cdots &{\frac{1}{M}\sum\limits_{j = 1}^M {{\mathit{\boldsymbol{A}}_{Lj}}} } \end{array}} \right) $

(13)

相邻多帧样本的分解模型中都使用同一组a、b、c、d超参数，这样相邻多帧样本对该组参数的更新过程都有影响，且这些参数也会反过来影响到相邻的多帧样本，这样相邻帧间的相关性可以在超参数的迭代更新过程中得以体现。

2.2 贝叶斯变分推断算法

在推断模型隐随机变量时，需要计算该隐随机变量后验分布，这往往涉及到复杂积分的求解。贝叶斯变分算法通过一个较易处理的概率分布来近似所需要求解的后验分布。在推断模型中的某个隐随机变量时，其他隐随机变量使用上一次的估计值，通过EM迭代算法使需推断的隐随机变量后验概率最大。经过多次迭代可以依次推断出模型的各个隐随机变量。

贝叶斯推断和EM算法涉及到多次迭代过程，计算量非常大。在Zhao等提出的MO-SBMTL快速算法^[15]的基础上，通过对超参数共享机制的合理设置，经过多次迭代过程，推断出各个隐变量的后验分布，最终用均值μ作为分解系数x的估计。

记相邻M帧信号{y⁽¹⁾, y⁽²⁾, …, y^(M)}基于同一离散傅里叶字典D的分解系数为{x⁽¹⁾, x⁽²⁾, …, x^(M)}。给隐随机变量a、b、c、d设置初值a₀、b₀、c₀、d₀。a₀=1/std(Y)²，其中列向量Y= (y^(1)T y^(2)T … y^(M)T)^T，std表示求标准差。b₀=e^-3, c₀和d₀中的各个元素可设为e^-4。为了降低模型复杂度，a和c在迭代过程中保持不变，a=a₀+MN，c=c₀+1/M(1 1 … 1)_N×1。

x⁽ⁱ⁾的均值μ_xⁱ和Σ_xⁱ按照式(10)和(11)迭代更新。b的迭代更新公式为

$ b = \sum\limits_{i = 1}^M {\left( {{b_0} + \left\| {{\mathit{\boldsymbol{y}}^{\left( i \right)}} - \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{\mu }}_{{x^i}}}} \right\|_2^2 + {\rm{trace}}\left( {\mathit{\boldsymbol{D}}{\mathit{\Sigma }_{{x^i}}}{\mathit{\boldsymbol{D}}^{\rm{T}}}} \right)} \right)} $

(14)

噪声精确度α₀的迭代更新公式为

$ {a_0} = a/b/k $

(15)

式中：k值可自行选取，研究中取k=1.2。d的迭代公式为

$ \mathit{\boldsymbol{d}} = {\mathit{\boldsymbol{d}}_0} + {\left( {\begin{array}{*{20}{c}} {\sum\limits_{i = 1}^M {\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{P}}_{1i}}} }&{\sum\limits_{i = 1}^M {\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{P}}_{2i}}} }& \cdots &{\sum\limits_{i = 1}^M {\mathit{\boldsymbol{D}}{\mathit{\boldsymbol{P}}_{Li}}} } \end{array}} \right)^{\rm{T}}} $

(16)

其中矩阵DP的第i列为

$ \mathit{\boldsymbol{D}}{\mathit{\boldsymbol{P}}_i} = \frac{1}{M}{\left[ {\begin{array}{*{20}{c}} {{\rm{abs}}{{\left( {{\mathit{\boldsymbol{\mu }}_{{x^1}}}^{\rm{T}}{\mathit{\boldsymbol{\mu }}_{{x^1}}} + {\rm{diag}}\left( {{\mathit{\Sigma }_{{x^1}}}} \right)} \right)}^{\rm{T}}}}\\ {{\rm{abs}}{{\left( {{\mathit{\boldsymbol{\mu }}_{{x^2}}}^{\rm{T}}{\mathit{\boldsymbol{\mu }}_{{x^2}}} + {\rm{diag}}\left( {{\mathit{\Sigma }_{{x^2}}}} \right)} \right)}^{\rm{T}}}}\\ \vdots \\ {{\rm{abs}}{{\left( {{\mathit{\boldsymbol{\mu }}_{{x^M}}}^{\rm{T}}{\mathit{\boldsymbol{\mu }}_{{x^M}}} + {\rm{diag}}\left( {{\mathit{\Sigma }_{{x^M}}}} \right)} \right)}^{\rm{T}}}} \end{array}} \right]^{\rm{T}}} $

(17)

α的迭代公式为

$ {\mathit{\boldsymbol{\alpha }}_{ij}} = \frac{{{\mathit{\boldsymbol{c}}_{ij}}}}{{{\mathit{\boldsymbol{d}}_{ij}}}},1 < i < L,1 < j < M $

(18)

若迭代次数达到上限，停止迭代。用均值μ_xⁱ作为分解系数x⁽ⁱ⁾的估计。

2.3 多帧联合稀疏特征的提取方法

将一段连续时间的信号分帧并作去均值和能量归一化预处理后，共T个帧信号样本。以连续相邻的M帧信号进行分组，假设M帧连续帧信号的下标为d₁, d₂, …, d_M，初始时设置d₁=1, d₂=2, …, d_M=M。采用贝叶斯变分算法推断模型中的隐随机变量，求解该组M帧信号基于离散傅里叶字典的M个分解系数。然后更新d₁=2, d₂=3, …, d_M=M+1，重复以上过程，直到d_M>T时停止。最后一共有T-M+1组。每一组包含M帧信号及其对应的M个分解系数，取最中间的一组分解系数(若M为偶数，取x^(M/2)；若M为奇数，取x^(M+1)/2)，对分解系数取幅值，并作能量归一化处理，作为各组M帧信号的多帧联合结构化稀疏特征，这样一共有T-M+1个特征样本。在本次研究中，为兼顾运算量和分类性能取M=3，分解系数能量归一化到0~1。在构造离散傅里叶字典D时，各傅里叶变换基的选取需对应于包含有线谱的频段范围。

多帧联合稀疏特征能够在一定程度上恢复和加强强能量窄带线谱频率成分，同时抑制能量较弱的时变性强的连续谱噪声。在测试样本加有高斯白噪声的失匹配条件下，该特征仍然能较好的恢复出强能量窄带线谱频率成分。

图 2为加入高斯白噪声后信噪比为-5 dB时，A类目标在50~200 Hz间短时傅里叶时频谱图。在相同条件下，选取50~200 Hz频段，分辨率为1 Hz的151个傅里叶基构成离散傅里叶字典D。图 3是将加噪信号基于字典D分解得到的多帧联合稀疏特征的示意图。通过对比可以发现在-5 dB的低信噪比条件下，多帧联合稀疏特征依然能够较好地恢复出强能量窄带线谱频率成分。

3 仿真实验和结果分析

实验中使用的数据集包含4类目标。数据首先经过分帧，去均值和能量归一化等预处理过程。共有2 400个样本，每帧时长1 s。训练样本数占总样本数的1/3。提取特征后，将特征样本和其相应的类标送入SVM分类器进行分类。训练期间通过交叉验证找到最佳模型参数用于测试阶段。

由于稀疏特征是通过基于离散傅里叶字典的分解得到的，为了进行对比，选择与字典中傅里叶变换基相对应频率点处的幅度谱系数，将其能量归一化后作为对比基准特征。同时选取了常用的MFCC特征作为对比。在送入SVM分类器前，所有待测试特征都会做归一化处理。图 4是多帧联合结构化稀疏特征与相应频段内的幅度谱特征以及MFCC特征的识别分类正确率对比图。

图 4中clean标注表示测试样本不添加额外的白噪声，和训练样本保持同一噪声环境的情况。从图 4中可以看出，多帧联合稀疏特征在所有测试条件下的识别正确率均高于幅度谱特征，在高信噪比失匹配条件下识别正确率与MFCC特征持平。在低信噪比条件下，MFCC特征识别率大幅下降，对于加性白噪声的鲁棒性较差，而多帧联合稀疏特征仍然能够保持较高的分类识别正确率，是对噪声具有鲁棒性的特征。

4 结论

1) 本文提出的多帧联合结构化稀疏特征，在训练样本和测试样本噪声环境不匹配(失匹配)条件下，能够有效恢复和增强窄带线谱成分，对噪声具有鲁棒性。

2) 通过一组实测数据进行分类的实验表明，在训练样本和测试样本噪声环境相同时，幅度谱特征、基于倒谱的MFCC特征和多帧联合结构化稀疏特征都有较高的识别正确率。但是在测试样本加入白噪声的失匹配条件下，幅度谱特征和MFCC特征识别率均出现明显下降，而多帧联合结构化稀疏特征仍能够保持较高的识别率，表明该特征在测试样本和训练样本失匹配条件下对噪声具有鲁棒性。

提取具有噪声鲁棒性的特征一直是水下目标识别的重要内容，本研究从目前热门的稀疏分解理论出发，利用结构化稀疏算法提取基于幅度谱上的特征，是一种新的尝试。本文利用幅度谱上的窄带线谱结构提取稀疏特征，需事先选定合适的频率范围，此外，大量迭代计算需要较高的计算资源，这给实际应用带来了困难，改进和优化快速算法是下一步要解决的问题。