0 引　言

主动声呐系统具有海洋环境中目标信号检测、目标估计和定位、目标分类、识别和目标跟踪等功能^[1]，影响主动声呐探测性能的因素包括船舶噪声、海洋环境噪声，混响是主要因素。混响来自海水不均匀性引起的声散射，以及海面和海床的散射。目前，主动声呐信号处理中混响抑制是一个重要问题。

学者们提出了很多抗混响方法，主要包括基于AR预白化的去混响方法^[2]、基于归一化匹配滤波的去混响方法^[3]、基于主成分求逆算法（principal component inversion，PCI）的抗混响方法^[4-5]以及基于低秩矩阵分解的抗混响方法^[6-7]等。舒象兰等^[2]在AR预白化的基础上，对接收信号进行二分法奇异值分解，以更好地抑制混响。但有一定的局限性，使用AR预白化方法时，确定AR系数的阶数非常困难。陈益明^[3]提出了混合块归一化匹配滤波算法，减少了由于混响非平稳性所导致的标准块归一化匹配滤波输出虚峰。但该方法只能一定程度上减弱相关峰，并不能完全抑制。许策、刘贯凌等^[4-5]利用PCI算法，根据混响和信号加背景噪声的功率差，将接收信号分为信号加背景噪声子空间和混响子空间，然后实现混响分离。但PCI方法很难估计秩，并且需要特定的先验知识。Wright等^[6-7]利用低秩矩阵分解法将接收信号分为低秩矩阵和稀疏矩阵。低秩矩阵是混响子空间，稀疏矩阵是将接收信号与混响分离的信号子空间。但低秩矩阵分解法计算量较大，在某些实用性要求较高的场合并不适用。

稀疏表示的关键问题之一是过完备字典的设计，这在很大程度上决定了信号的特征能否被有效表示。字典有固定完整字典和学习字典2种类型^[8]。固定完整字典主要通过一些数学函数变换，利用参数选择与信号相对应的原子，且要求原子数远大于信号长度。这种类型的常见形式包括离散小波变换（discrete wavelet transform，DWT）字典、离散余弦变换（discrete cosine transform , DCT）字典和Gabor字典。学习字典在固定词典的基础上添加了字典学习算法，并获得了第1类字典的最佳形式^[9]。字典学习是稀疏表示的一个典型示例，其目的是识别从信号本身学习字典中一组信号的稀疏表示^[10-11]。在迭代过程中，字典会根据给定的误差不断训练和更新原子，使其更接近信号的特征。与固定字典相比，使用训练字典可以降低算法的复杂度，更好地匹配信号的特征。

非负矩阵分解(non-negative matrix factorization, NMF) 是一种多变量子空间分析方法，自提出以来便得到深入研究和广泛应用^[12]。其核心思想是在一定约束条件下将原始非负矩阵分解为2个低维的非负矩阵^[13]。这种分解矩阵的非负性给这种描述数据的方式带来了一定的稀疏性。稀疏表示通过将主动声呐回波信号在所训练的字典上进行稀疏表示^[14]，使得信号具有一定稀疏性。因此，考虑将NMF应用到稀疏表示中，这样可以进一步提高主动声呐抗混响的性能。

本文介绍稀疏表示方法，提出一种基于稀疏表示的非负矩阵分解(non-negative matrix factorization, NMF)抗混响方法，并进行仿真模拟分析。

1 主动声呐回波信号模型

用于水下主动声呐目标检测的声呐接收人造目标反射回波、海洋环境噪声和周围环境中目标随机散射形成的混响。将携带目标参数信息的目标反射回波从混响和噪声中分离出来是主动声呐检测的主要任务。

亮点模型是工程应用中水下目标散射回波形成机理的高频近似简化模型。该模型假设复杂目标由多个散射亮点组成，总目标回波是每个亮点产生的散射回波分量的线性叠加。亮点分为几何亮点和弹性亮点。几何亮点回波主要包括光滑表面上的镜面反射回波和角反射回波，这与物体大小、入射角有关。这种回波具有很强的能量，在实际目标检测中很容易观测到。弹性亮点回波是声波与介质界面相互作用的响应，取决于物体的材料特性和内部结构。考虑到这些回波的传播衰减较大，且回波能量通常较小，在混响信号（SRR）较低时难以观测。因此本文仅考虑几何亮点。

高光模型将目标视为一个线性系统，目标散射回波信号的表达式用传递函数表示。对于图1所示的目标结构，对应于目标的每个亮点显示的传递函数可以定义为：

$ H(\omega ) = \sum\limits_{n = 1}^N {{A_n}(r,\omega ){e^{j\omega {\tau _n}}}{e^{j{\varphi _n}}}}。$

(1)

其中： $ {A_n} $ 为振幅反射系数； $ \tau $ 为延迟； $ \varphi $ 为目标边缘不连续处或亮点类别变化期间存在的相位跳变； $ m $ 是亮点数量。对于发射信号 $ s(t) $ ，目标回波 $ x(t) $ 表示为：

$ x(t) = \sum\limits_{n = 1}^N {{x_n}(t)} = \sum\limits_{n = 1}^N {{A_n}s(t - {\tau _n}){e^{j{\varphi _n}}}} 。$

(2)

采用单元散射模型描述混响的形成机理。混响被模拟为入射声波散射回波的线性叠加，这是由于许多水下散射而产生的：

$ R(t) = \sum\limits_{i = 1}^{N(t)} {K(t)s} (t - {\tau _i}){e^{j{\varphi _i}}}。$

(3)

其中： $ K(t) $ 为 $ t $ 时刻混响的随机振幅因子； $ N(t) $ 为 $ t $ 时刻产生混响的散射体数量； $ {\varphi _i} $ 是随机相位。主动声呐接受的回波信号用矩阵可以表示为：

$ {\boldsymbol{x}}(t) = {\boldsymbol{As}}(t) + {\boldsymbol{n}}(t) = {{\boldsymbol{x}}_T}(t) + {\boldsymbol{r}}(t) + {\boldsymbol{n}}(t) 。$

(4)

其中， $ {\boldsymbol{A}} $ 为从散射单元到接受系统的水下传输信道确定的混合矩阵。

2 稀疏表示 2.1 OMP算法

正交匹配追踪（OMP）算法是一种基于稀疏分解理论的稀疏信号表示算法^[15]，由匹配追踪（MP）算法改进而来。匹配追踪算法是贪婪算法（greedy algorithm）的一种，其通过每次在过完备字典中选择一个与原始信号最匹配的最优原子^[16]，将这些最优原子进行线性组合，从而表示原始信号。

OMP算法过程如下：

输入： $ {\text{N}} \times {\text{d}} $ 测量矩阵 ${\boldsymbol{\varPhi}}$ ， ${\boldsymbol{N}}$ 维数据向量 $ v $ ，理想信号稀疏度m；

输出：对 $ s $ 的稀疏表示 $ \mathop s\limits^ \wedge $ 。

具体步骤如下：

步骤1　初始化残差 $ {r_0} = v $ ，索引集 $ {\Lambda _0} = \Phi $ 和迭代计数器 $ t = 1 $ 。

步骤2　通过计算残差 $ r $ 与传感矩阵的列 $ {\varphi _j} $ 的积，选择的原子满足：

$ {\lambda _t} = \arg {\max _{j = 1,...,d}}\left| {\left\langle {{r_{t - 1,}}\left. {{\varphi _j}} \right\rangle } \right.} \right| 。$

(5)

步骤3　更新所选原子的索引集 $ {\Lambda _t} = {\Lambda _{t - 1}} \cup \left\{ {{\lambda _t}} \right\} $ ，记录找到的传感矩阵的重建原子集合：

$ {{\boldsymbol{\varPhi}} _t} = \left[ {{{\boldsymbol{\varPhi}} _{t - 1}}\left. {{\varphi _{{\lambda _t}}}} \right]}。\right. $

(6)

步骤4　更新稀疏矢量，由最小二乘法得：

$ {\widehat s_t} = \arg {\min _x}{\left\| {v - {{\boldsymbol{\varPhi }}_t}\hat s} \right\|_2} 。$

(7)

步骤5　计算新残差：

$ {r_t} = v - {{\boldsymbol{\varPhi}} _t}{\hat s_t} 。$

(8)

步骤6　判断是否满足 $ t < m $ 。若不满足，则停止迭代；若满足，则返回步骤2。

2.2 字典学习

信号的稀疏度很大程度上由字典决定，因此字典的选择至关重要^[17]。在进行稀疏表示时往往会先初始化一个字典，但它一般不是最优的。如固定的字典（DCT字典、Gabor字典等）对某些信号进行稀疏表示时效果很好，但不适用于其他信号。因此，有必要为特定的信号找到最合适的字典^[18]。为此，将使用字典学习算法获取主动声呐回波信号在字典下的稀疏表示。

K-SVD算法2006年由Aharon^[19]等提出。该算法推广了K-means聚类过程，通过在稀疏编码和更新字典原子之间交替执行训练超完备字典，以更好地表示数据。K-SVD算法过程为：

当给出训练集 $ Y = \left\{ {{y_i} \in {\text{ }}{R^m}} \right\}_{i = 1}^N $ ，通过解决以下优化问题，可以找到具有 $ K $ 个原子的最优字典 ${{D}}$ 稀疏表示 $ Y $ 。

$ \mathop {\min }\limits_{D,X} \left\{ {\left\| {Y - DX} \right\|_F^2} \right\}{\text{ s}}{\text{.t }}\forall i,{\left\| {{x_i}} \right\|_0} \leqslant {T_0} 。$

(9)

其中 ${{D}} \in {{{R}}^{m \times K}}(K \gg m)$ ， $ X \in {R^{K \times N}} $ 为 $ Y $ 的稀疏编码， $ {T_0} $ 为 $ X $ 的稀疏性阈值。

首先，设置具有 $ {l^2} $ 个归一化列的初始字典矩阵 ${{{D}}_0} \in {{{R}}^{m \times K}}(K \gg m)$ ，然后，重复2个阶段直到收敛。

1）稀疏编码阶段

使用追踪算法计算每个示例 $ {y_i} $ 对应表示向量 $ {x_i} $ ，如匹配追踪(MP) 或正交匹配追踪(OMP)。与最小二乘法相比，OMP具有更快的恢复速度，易于实现且速度更快。因此，本文选择OMP算法。最佳稀疏编码 $ X $ 可以通过下式获得:

$ X = \arg \mathop {\min }\limits_{{x_i}} \left\| {{y_i} - D{x_i}} \right\|_F^2{\text{ s}}{\text{.t }}\forall i,{\left\| {{x_i}} \right\|_0} \leqslant {T_0} 。$

(10)

2）字典更新阶段

通过求解一个小的奇异值分解（SVD）问题，一次只更新 ${{D}}$ 的一列。

$ {{\boldsymbol{E}}_k} = Y - \sum\limits_{j \ne k} {{d_j}x_T^j} ，$

(11)

$ {\boldsymbol{E}}_k^R = U\Delta {V^{\rm{T}}} 。$

(12)

其中： ${{\boldsymbol{E}}_k}$ 为整体表示误差矩阵； ${\boldsymbol{E}}_k^R$ 为 $ {E_k} $ 的限制矩阵； $ x_T^j $ 为 $ X $ 的第 $ j $ 行（对应于 ${{D}}$ 第 $ j $ 列的系数）。 $ U $ 的第一列作为D的更新列， $ V $ 的第一列乘以 $ \Delta (1,1) $ 更新相应的系数向量。

由SVD算法步骤可以看出，该算法计算量小，可以刻画出数据的重要特征，因此实际应用范围很广。

3 基于NMF的协方差估计

通过稀疏表示将主动声呐回波信号变换到稀疏域，改善信号的稀疏性。增加输出信号的稀疏性产生接近目标回波的信号，并根据这一特性，将其应用于去混响。

用NMF对稀疏表示信号 $ \widehat {\text{s}} $ 进行处理，首先设当前信号矩阵为 $ M $ ，则其低维形式表示为：

$ {\boldsymbol{M}} \approx {\boldsymbol{WH}},{{ {\boldsymbol{W}}}} \in {\boldsymbol{R}}_{0 + }^{i \times k},H \in {\boldsymbol{R}}_{0 + }^{k \times j}，$

(13)

协方差估计矩阵为：

$ V = WH 。$

(14)

式中： $ W $ 为基矩阵； $ H $ 为系数矩阵；且两者都是非负矩阵。

在NMF分解过程中， $ W $ 和 $ H $ 被限制为非负，实际对 $ M $ 的分解时，分解前后一定会存在误差。因此，非负矩阵分解过程是一个优化过程，需要选择目标函数结合该函数下的迭代规则进行矩阵分解。简单来说就是优化非负矩阵 $ {\boldsymbol{M}} $ ，以尽可能实现式（13）。

以KL为代价函数优化NMF的问题可以总结如下：

$ \min {D_{KL}}(M||WH) = \sum\limits_{i = 1}^I {\sum\limits_{j = 1}^J {\left( {{M_{ij}}\ln \frac{{{M_{ij}}}}{{{{\left[ {WH} \right]}_{ij}}}} - {M_{ij}} + {{\left[ {WH} \right]}_{ij}}} \right)} } ，$

(15)

梯度下降算法为：

$ \begin{gathered} {W_{ik}} = {W_{ik}} - {\mu _{ik}}\frac{{\partial {D_{KL}}(M||WH)}}{{\partial {W_{ik}}}}，\\ {H_{kj}} = {H_{kj}} - {\eta _{kj}}\frac{{\partial {D_{KL}}(M||WH)}}{{\partial {H_{kj}}}} 。\\ \end{gathered} $

(16)

其中：

$ \begin{gathered} \frac{{\partial {D_{KL}}(M||WH)}}{{\partial {W_{ik}}}} = - \sum\limits_{j = 1}^J {({H_{kj}}{M_{ij}}/{{[WH]}_{ij}} + {H_{kj}})}，\\ \frac{{\partial {D_{KL}}(M||WH)}}{{\partial {H_{kj}}}} = - \sum\limits_{j = 1}^J {({W_{ik}}{M_{ij}}/{{[WH]}_{ij}} + {W_{ik}})} 。\\ \end{gathered} $

(17)

根据梯度下降算法，转化为乘法算法，令

$ {\mu _{ik}} = \frac{1}{{\sum\limits_{j = 1}^J {{H_{kj}}} }},{\text{ }}{\eta _{kj}} = \frac{1}{{\sum\limits_{i = 1}^I {{W_{ik}}} }} ，$

(18)

将式（17）和式（18）分别代入式（16）中，则

$ \begin{gathered} {W_{ik}} = {W_{ik}}\frac{{\sum\limits_{j = 1}^J {{H_{kj}}{M_{ik}}/{{[WH]}_{ik}}} }}{{\sum\limits_{j = 1}^J {{H_{kj}}} }}，\\ {H_{kj}} = {H_{kj}}\frac{{\sum\limits_{i = 1}^I {{W_{ik}}{M_{ik}}/{{[WH]}_{ik}}} }}{{\sum\limits_{i = 1}^I {{W_{ik}}} }}。\\ \end{gathered} $

(19)

最后，将式（19）代入式（14），得到协方差估计矩阵 $ {\boldsymbol{V}} $ 。

去混响算法流程如图2所示。首先对主动声呐回波信号进行稀疏表示得到稀疏信号 $\widehat {{s}}$ ，然后对 $\widehat {{s}}$ 所构成的 $ {\boldsymbol{M}} $ 进行非负矩阵分解，最后利用式（14）求出协方差矩阵 $ {\boldsymbol{V}} $ 。

4 仿真验证

对基于稀疏表示的NMF去混响方法有效性进行仿真验证，发射信号为线性调频脉冲信号，频率范围为25~50 kHz，采样频率为500 kHz，发射信号脉宽为2 ms。

主动声呐目标回波仿真如图3所示。可以看出，对于主动声呐目标回波信号而言，对其进行稀疏表示是可行的，稀疏表示结果已经包含回波信号的大部分能量。因此，可以证明主动声呐回波可以进行稀疏表示，在混响背景下使用NMF抑制混响也具有可行性。

定义稀疏表示的相对均方误差为：

$ {{error}} = \frac{{\left\| {s - \mathop s\limits^ \wedge } \right\|_2^2}}{{\left\| s \right\|_2^2}}。$

(20)

其中： $ s $ 为主动声呐实际接收到的主动声呐回波信号； $ \mathop s\limits^ \wedge $ 为变换到稀疏域的稀疏信号。当稀疏向量中非零元素个数分别为0，100，200，300，400，500时，误差如图4所示。令非零元素个数为 $ N $ ， $ N $ 代表了稀疏表示信号中主要成分数量多少。当 $ N $ 比较小时，因为主要成分数量较少，所以误差比较大；随着非零元素数量的增大，主要成分数量也在增加，误差逐渐减小。但是当 $ N $ 进一步增大时，主要成分数量增加的同时冗余信息也在增加，因此信号的误差增大。

图5为原始信号以及采用预白化、常规匹配滤波和本文新方法对含混响的主动声呐回波信号进行处理后的时频图。

由图5（a）可以看出，亮斑1、亮斑2和亮斑3分别为目标亮点以及混响。由图5（b）可以清楚看到，经本文方法处理后，亮斑1即目标回波实现了增强，而混响2和混响3均被抑制。图5（c）和图5（d）分别为AR预白化、常规匹配滤波去混响后的信号时频分布，对比图5（b）可以明显看出，这2种方法虽去除了部分混响，但去混响效果远没有本文方法处理效果好。结果表明，该算法能有效抑制混响，提取目标回波成分。

进一步分析3种去混响方法，通过1000次的蒙特卡罗实验，在虚警概率不大于0.1的前提条件下，得到的性能比较结果如图6所示。

可知，本文方法优于其他2种方法。

5 结　语

针对主动声呐抗混响问题，本文提出基于稀疏表示的NMF抗混响方法。充分考虑稀疏信号的稀疏性以及NMF的特性，计算声呐接收的目标回波信号的协方差，即以KL为代价函数，构造一个NMF优化问题，利用梯度下降法算法估计协方差，得到了协方差矩阵的估计值，进一步提高主动声呐目标回波去混响性能。实验结果表明，相对于其他去混响算法，本文提出的基于稀疏表示的NMF去混响方法，在一定程度上提高了主动声呐去混响性能。