0 引　言

回波是发射信号与目标相互作用的结果，必然包含目标特征，因此通过分析回波信号的特征可以完成对水下目标的检测^[1]。而混响是由声波在不规则散射体、海面海底上散射产生的^{[2 − 3]}，是影响主动声呐目标检测的重要干扰之一。在接收信号的混响抑制问题上，传统的混响抑制方法都是基于单帧的混响特性进行处理，而忽视了前后帧数据间混响的关联特性。事实上，混响在波束空间中的特性和各帧之间的关联性具有很大的利用价值，如果充分利用多帧波束空间混响数据内在属性，即充分利用多帧波束空间混响数据统计相关信息，则可以有效提高主动声呐探测时对混响的抑制能力，近年来兴起的低秩稀疏矩阵分解算法是一个很好的研究方向。Toh K-C等^[4]采用具有较快收敛速率的加速近端梯度法（Accelerated Proximal Gradient, APG）解决了以凸优化问题为基础建立的矩阵低秩稀疏分解问题。陈明等^[5]进一步考虑了同一凸优化问题的对偶问题，设计了基于对偶问题实现数据分离的算法，并考虑以增广拉格朗日乘子法解决类似的问题^[6]。周涛等^[7]以快速双边随机投影法迭代解决具有低秩稀疏非凸约束的优化分解问题。

近年来，低秩稀疏矩阵分解在水声领域的尝试和突破越来越多，主要体现在不同场景混响背景下的应用以及各种分解算法的实现。李伟昌等^[8]提出采用联合子空间和稀疏过滤算法，解决传统子空间技术不适用于强混响背景下高速移动目标的检测问题，该算法将原始数据分解成低秩、稀疏和噪声3个部分，提高了低信噪比情况下的移动目标跟踪能力。葛凤翔等^[9]利用混响背景的潜在相干结构，将波束空间图像数据表示为混响低秩矩阵和目标信号稀疏矩阵的和，将检测和跟踪表述为一个结构化的凸优化问题，并通过加速近端梯度（APG）算法求解，恢复了相干背景混响和运动目标信号对应的低秩和稀疏分量。刘建设等^[10]重点研究了基于低秩稀疏矩阵分解算法的冰下运动目标主动声呐探测方法，将混响与运动目标分离，并对比了传统的背景差分法，结果表明低秩稀疏矩阵分解算法有更好的适应性。朱广平等^[11]根据目标回波信号与混响信号在时频域上能量分布的相关性不同，利用自适应核时频分析方法将接收信号转换到时频域上进行分析，通过低秩稀疏矩阵分解算法将目标回波与混响分别分解到稀疏矩阵和低秩矩阵中，实现目标回波与混响的分离，降低了混响对回波信号的干扰。吴耀文等^[12]针对水声信号的非线性、非高斯、非平稳和低信噪比等特点，将基于RPCA算法的Godec算法、非负矩阵分解算法和维纳滤波应用于低频段的水声场景，建立分解模型，将接收信号表示为低秩、稀疏、噪声3个部分，对低频段信号降噪难题提出了有效的解决方案。陈勇^[13]通过分析数据协方差矩阵的低秩稀疏矩阵分解问题，提出了一种基于最小方差无失真响应和低秩稀疏矩阵分解的波达方向估计方法，提高了波束形成器波达方向估计的算法性能。贺玉梁^[14]提出了低秩稀疏矩阵分解的各类离线算法和在线算法，并且对比总结了PCA、APG和FDPM算法的性能优劣性。

以上研究都取得了较好的处理效果，但没有完全抑制混响以及噪声的影响，处理结果依然存在不同程度的随机干扰。本文在低秩稀疏矩阵分解的基础上利用多帧之间目标运动的连续特性，提出一种定位窗滤波方法，滤除了稀疏杂点，净化主动声呐显示图像，提高对动目标的自主检测性能。

1 低秩稀疏矩阵分解算法 1.1 PCA算法

PCA（Principal Component Analysis）算法是一种常见的数据处理方法。PCA通过线性变换将原始数据变换为一组各维度线性无关的表示，可用于提取数据的主要特征分量，常用于高维数据到低维空间的降维运算。原始数据矩阵$ \boldsymbol{X} $的前r个最大的奇异值对应的左特征向量展成的特征子空间近似为混响的主要成分确定的子空间。矩阵$ \boldsymbol{X} $先进行奇异值分解，即$ \boldsymbol{X}={U}{\Sigma }{{V}}^{{{\mathrm{T}}}} $，其中$ {\boldsymbol{U}} \in {{\boldsymbol{R}}^{m\times m}} $和$ {\boldsymbol{V}} \in {{\boldsymbol{R}}^{m \times n}} $均为正交矩阵，$ {{\boldsymbol \Sigma} _r} = diag ({\sigma _1},{\sigma _2},...,{\sigma _r}) \in {{\boldsymbol{R}}^{r\times r}} $为对角矩阵，其对角线元素为正数且随下标递减，则矩阵$ {\boldsymbol{U}} $近似确定了混响主成分的子空间，分离出的各成分可简单表示为：

$ {\boldsymbol{A}} = {{\boldsymbol{U}}_R}{\boldsymbol{U}}_R^{\mathrm{T}}{\boldsymbol{X}}，$

(1)

$ {\boldsymbol{E}} = {{\boldsymbol{P}}_\Omega }({\boldsymbol{X}} - {\boldsymbol{L}})，$

(2)

$ {\boldsymbol{G}} = {\boldsymbol{X}} - {\boldsymbol{A}} - {\boldsymbol{E}}。$

(3)

式中：矩阵A为混响背景，是一个低秩矩阵；矩阵E为目标信号，是一个稀疏矩阵；矩阵G为噪声矩阵；$ {{\boldsymbol{P}}_\Omega }(z) $表示一个已知的投影算子；$ \Omega $是列向量z的前k个元素下标组成的集合，该投影保持$ \Omega $索引的元素不变，其余均为零元素。输入为矩阵时，其代表对该矩阵的每一列进行相应操作。

1.2 鲁棒PCA算法

考虑到实际子空间成分相互泄露，以及存在其余的成分，会降低分解算法的效果，采用鲁棒PCA（RPCA）算法。定义混响背景矩阵为A，回波信号强度矩阵为E，则在混响下，A是低秩矩阵，E是稀疏矩阵，它们的叠加即为接收机接收信号，定义为D=A+E，从分解算法设计的角度出发，矩阵D也称为原始矩阵。则优化问题如下：

$ \mathop {{\mathrm{min}}}\limits_{A,E} rank({\boldsymbol{A}}) + \lambda ||{\boldsymbol{E}}|{|_0} 。$

(4)

其中，$ {\boldsymbol{A}} + {\boldsymbol{E}} = {\boldsymbol{D}} $。

式(4)中目标函数为矩阵A的秩与矩阵E的零范数，$ \lambda $表示回波信号强度所占的权重，而式中秩函数和零范数均为非凸函数，需要对问题进行转松弛。

经过对实际应用场景的分析，利用矩阵的核函数来近似矩阵的秩，用矩阵的1-范数来近似0-范数，由范数相关理论知，两者皆为凸函数，则上述优化问题转化为：

$ \mathop {{\mathrm{min}}}\limits_{A,E} ||A|{|_*} + \lambda ||E|{|_1}。$

(5)

其中，$ A + E = D $。

对该问题的解读如下：矩阵A的秩取最小时，认为A是低秩的；矩阵E是稀疏的，则用E的0-范数，即E中非零元素的个数表征，当矩阵中非0元素的个数越少，E就越稀疏，因此得到这个优化问题。其中λ指的是回波信号强度所占的权重，是一个已知量，但是在这个式子中的2个函数都是非凸的，所以无法求出最优解。因此用矩阵A的核范数（矩阵A的奇异值之和）来近似A的秩，此时核范数越小时，可近似认为矩阵A的秩也越低；而用矩阵E的1范数来近似矩阵E的0-范数时，因为1-范数是指矩阵中某一列元素的绝对值之和取最大时，这个求和值记为1-范数，那么假如这个和已经很小了，趋向于0，则说明该矩阵其他列元素的绝对值之和也很小，因此可以近似认为非0元素的个数很少。由此得到最终的优化问题式(5)，两者均为凸函数，求和后存在最优解。

对优化问题式(5)，利用增广拉格朗日乘子法来求解。构造拉格朗日函数如下：

$ \begin{split} {\boldsymbol{L}}({\boldsymbol{A}},{\boldsymbol{E}},{\boldsymbol{Y}},\mu ) = & ||{\boldsymbol{A}}|{|_*} + \lambda ||{\boldsymbol{E}}|{|_1}+ < {\boldsymbol{Y}},{\boldsymbol{D}} - {\boldsymbol{A}} - {\boldsymbol{E}} >+ \\ & \frac{\mu }{2}||{\boldsymbol{D}} - {\boldsymbol{A}} - {\boldsymbol{E}}||_2^{\boldsymbol{F}} 。\end{split} $

(6)

式中：Y为拉格朗日乘子（足够大才能使得A和E逼近D）；<X，Y>代表矩阵X和Y的内积运算；$ \mu $为将有约束问题（D=A+E）转化为无约束问题时产生的代价因子，是一个正标量（足够大），该项称为惩罚项；$ \left|\right|{\boldsymbol{X}}|{|}_{{F}} $代表矩阵X的Frobenius范数。

对式(6)固定其他量，分别化简变量为A和E对应的最小值函数表达式：

1）固定A，Y，$ \mu $对目标变量为E的函数进行化简：

$ \begin{gathered} \mathop {{\mathrm{min}}}\limits_{{\boldsymbol{A}},{\boldsymbol{E}}} ||{\boldsymbol{A}}|{|_*} + \lambda ||{\boldsymbol{E}}|{|_1} + < {\boldsymbol{Y}},{\boldsymbol{D}} - {\boldsymbol{A}} - {\boldsymbol{E}} >+ \frac{\mu }{2}||{\boldsymbol{D}} - \\ {\boldsymbol{A }}- {\boldsymbol{E}}||_2^F \geqslant \mathop {{\mathrm{min}}}\limits_{\boldsymbol{E}} \frac{\lambda }{\mu }||{\boldsymbol{E}}|{|_1} + \frac{1}{2}||{\boldsymbol{E}} - \left({\boldsymbol{D}} - {\boldsymbol{A}} + \frac{{\boldsymbol{Y}}}{\mu }\right)||_F^2 。\\ \end{gathered} $

(7)

2）固定E，Y，$ \mathrm{\mu } $对目标变量为A的函数进行化简：

$ \begin{gathered} \mathop {{\mathrm{min}}}\limits_{{\boldsymbol{A}},{\boldsymbol{E}}} ||{\boldsymbol{A}}|{|_*} + \lambda ||{\boldsymbol{E}}|{|_1} + < {\boldsymbol{Y}},{\boldsymbol{D}} - {\boldsymbol{A}} - {\boldsymbol{E }}>+ \frac{\mu }{2}||{\boldsymbol{D}} - \\{\boldsymbol{A}} - {\boldsymbol{E}}||_2^F{\geqslant} \mathop {{\mathrm{min}}}\limits_{\boldsymbol{A}} \frac{1}{\mu }||{\boldsymbol{A}}|{|_*} + \frac{1}{2}||{\boldsymbol{A}} - \left({\boldsymbol{D}} - {\boldsymbol{E}} + \frac{{\boldsymbol{Y}}}{\mu }\right)||_F^2 。\\ \end{gathered} $

(8)

则式(7)和式(8)是原函数表达式分别对A和E求最小值后的化简式子，这2个式子取到最小值时对应的A和E即为相应分解结果。关于式(7)和式(8)的最小值求解，分别考虑如下两类优化问题：

$ {\rm{arg}}\mathop {\rm{min}}\limits_{{X}} \varepsilon ||{\boldsymbol{X}}|{|_1} + \frac{1}{2}||{\boldsymbol{X}} - {\boldsymbol{M}}||_F^2 ，$

(9)

$ {\rm{arg}}\mathop {\rm{min}}\limits_{{X}} \varepsilon ||{\boldsymbol{X}}|{|_*} + \frac{1}{2}||{\boldsymbol{X}} - {\boldsymbol{M}}||_F^2，$

(10)

通过求解，利用软阈值函数：

$ {T}_{\varepsilon}(x)= \left\{\begin{array}{l} x-\varepsilon ，x > \varepsilon ，\\ x+\varepsilon ，x < \varepsilon ，\\ 0，其他。\end{array}\right. $

(11)

在化简式(9)的过程中，不仅要X的1范数取最小，还需要第2项尽可能等于0 ，这样关于X的这个函数才会达到最小，所以近似认为X = M，这样后面才会等于0。因为近似两者相等，所以对M进行操作的时候也就是对于X进行操作。$ \varepsilon $足够小，则当M中的元素大于$ \varepsilon $时就会减去$ \varepsilon $直到最后收敛到0，小于$ - \varepsilon $同理，那么就这样对$ {\boldsymbol M} $中的元素进行迭代更新，也就是近似对M进行更新，使得最终得到的X为一个稀疏矩阵。因此利用式(11)的函数，可得式(9)优化问题的解为：

$ {\boldsymbol{X}} = {{{T}}_\varepsilon }({\boldsymbol{M}}) 。$

(12)

同理，对式(10)中进行类似的运算，但如果对式(10)完全相同的操作，则会使得其结果也稀疏。因此对矩阵M先进行奇异值分解，即$ {\boldsymbol{M}} = {\boldsymbol{U}}\Sigma {{\boldsymbol{V}}^{\rm T}} $，其中$ {{\boldsymbol{U}}}\in{{\boldsymbol{R}}}^{{m}{\times}{m}} $和$ {{\boldsymbol{V}}}\in{{{\boldsymbol{R}}}}^{{m}{\times}{n}} $均为正交矩阵，$ {\Sigma _r} = diag ({\sigma _1},{\sigma _2},..., {\sigma _r}) \in {R^{r*r}} $为对角矩阵 ，其对角线元素为正数且随下标递减。同样利用式(11)的函数，可得式(10)优化问题的解为：

$ X = U{T_\varepsilon }(\Sigma ){V^{\rm T}} 。$

(13)

注意，上述过程中将式(11)的函数$ {T_\varepsilon }(x) $用到矩阵上时，指的是对该矩阵的每一个元素进行函数映射运算。对比式(9)和式(10)表达式以及解式(12)和解式(13)，可得式(7)和式(8)的解分别为：

$ {\boldsymbol{E}} = {{\boldsymbol{T}}_{\frac{\lambda }{\mu }}}({\boldsymbol{D}} - {\boldsymbol{A}} - \frac{{\boldsymbol{Y}}}{\mu })，$

(14)

$ {\boldsymbol{A}} = {\boldsymbol{U}}{{\boldsymbol{T}}_{\frac{\lambda }{\mu }}}(\Sigma ){{\boldsymbol{V}}^{\rm T}}。$

(15)

其中，$ {\boldsymbol{D}} - {\boldsymbol{A}} - \displaystyle\frac{{\boldsymbol{Y}}}{\mu } = {\boldsymbol{U}}\Sigma {{\boldsymbol{V}}^{\rm T}} $。

在分别更新完矩阵E和A之后，再对拉格朗日算子Y和代价因子$ \mu $进行更新：

$ {\boldsymbol{Y}} = {\boldsymbol{Y }}+ \mu ({\boldsymbol{D}} - {\boldsymbol{A}} - {\boldsymbol{E}})，$

(16)

$ \mu = \rho \mu 。$

(17)

当矩阵A和E收敛时，分解算法完成。

2 定位窗滤波方法

由于实际数据存在噪声、随机起伏等影响，分解的稀疏矩阵结果必然会存在很多稀疏杂点，但是这些稀疏杂点和目标点是存在本质特征区别的。显然，目标的运动物理连续，不会出现间断跳跃的情况，但是稀疏杂点的出现位置随机变化。针对此特征差异，本文提出一个定位窗方法，在低秩稀疏矩阵分解处理后进一步处理提高主动声呐图像显示。

实现流程如下：根据目标回波长度、屏幕刷新速率以及目标运动速度极限值，设置M×N的定位窗，每一个窗都聚焦前一帧稀疏矩阵的稀疏点，并且对下一帧矩阵同样位置的定位窗进行扫描。若存在稀疏点，即定位窗为非零矩阵，则认为存在目标点，将其保留；若不存在稀疏点，即定位窗为全零矩阵，则将其认定为随机杂点。多帧处理并重复上述过程，最终得到目标轨迹。

定位窗的大小设置很关键，如果过大则会包含更多的稀疏杂点，从而把这些点认定为目标点，如果过小则会使得下一帧的目标点超出定位窗框定范围，从而认定为稀疏杂点，被滤除。假设预检测的目标速度最大为$ {v_{\max }} $，单帧周期为$ T $，稀疏点中心坐标$ (X,Y) $，点大小为$ (m,n) $，其中行代表距离信息，列代表角度信息，最小探测距离为$ {d_{\min }} $，运动态势示意图见图1。

分析可知，当目标在最小探测距离圆的圆边时对应目标单帧理论最大移动角度，通过几何关系计算即可得到目标单帧的最大移动角度为：

$ {\theta _{\max }} = \frac{{{v_{\max }}T}}{{{d_{\min }}}}\frac{{{{180}^ \circ }}}{\text π} 。$

(18)

然而式(18)仅为理论上的极限边界值，实际若采用这个值来设定定位窗大小的话，势必过大。因此在实际场景中使用时，必须将最小探测距离$ {d_{\min }} $转变为期望探测距离$ {d_{\rm{exp ect}}} $：

$ {\theta _{\max }} = \frac{{{v_{\max }}T}}{{{d_{\rm{exp ect}}}}}\frac{{{{180}^ \circ }}}{\text π} 。$

(19)

即事先预估能够探测到动目标的最小距离，由此来提高定位窗的滤波能力。同样的，期望距离越近，定位窗就越大，其滤波能力也就越差。

综上所述，定位窗矩阵设定为：

$ \left[ {\begin{array}{*{20}{c}} {\left( {X - m - {v_{\max }}T,Y - n - {\theta _{\max }}} \right)}& \cdots &{}&{}&{}& \cdots &{\left( {X - m - {v_{\max }}T,Y + n + {\theta _{\max }}} \right)} \\ \vdots &{}&{}&{}&{}&{}& \vdots \\ {}&{}&{\left( {X - m,Y - n} \right)}& \cdots &{\left( {X - m,Y + n} \right)}&{}&{} \\ {}&{}& \vdots &{}& \vdots &{}&{} \\ {}&{}&{\left( {X + m,Y - n} \right)}& \cdots &{\left( {X + m,Y + n} \right)}&{}&{} \\ \vdots &{}&{}&{}&{}&{}& \vdots \\ {\left( {X + m + {v_{\max }}T,Y - n - {\theta _{\max }}} \right)}& \cdots &{}&{}&{}& \cdots &{\left( {X - m + {v_{\max }}T,Y + n + {\theta _{\max }}} \right)} \end{array}} \right]。$

(20)

3 仿真及试验研究 3.1 仿真分析

多帧混响仿真条件如下：单基元全向发射，16元标准线列阵接收，在水深50 m，分别考虑单频脉冲（CW）信号和调频线性（LFM）信号，中心频率$ {{f}}_{{c}}=30\ \mathrm{k}\mathrm{H}\mathrm{z} $，信号带宽$ {B}=5\ \mathrm{k}\mathrm{H}\mathrm{z} $，脉宽$ {T}=5\ \mathrm{m}\mathrm{s} $，声速$ {c}=1\ 500\ \mathrm{m}/\mathrm{s} $，采样频率$ {{f}}_{{s}}=200\ \mathrm{k}\mathrm{H}\mathrm{z} $，单帧接收信号截止时间$t_ {\mathrm{e}\mathrm{n}\mathrm{d}}=0.3\ \mathrm{s} $。发射信号的$ 80\ \mathrm{m}\mathrm{s} $内不接收回波，即单帧接收信号起始时间$ {t}_{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{r}\mathrm{t}}=0.08\ \mathrm{s} $。取11帧混响信号，并对多帧混响方位距离图信息进行联合。

运动目标仿真条件如下：

目标角度：[120,118,116,114,112,110,108,106,103,100,96]（单位：°）；

目标距离：[200,196,192,188,184,180,176,172,168,164,160]（单位：m）。

两者一一对应，构成了11帧目标运动信息，角度和距离向量的元素序号即为对应帧号。将基元端信混比设置为−5 dB，通过波束形成，分别得到发射信号为CW信号和LFM信号的目标方位距离图。原始数据以及PCA分解结果、RPCA分解结果如图2～图4所示。

由于单帧接收起始时间为$ 0.08\ {\mathrm{s}} $，因此距离轴起始距离显示为$ 60\ {\mathrm{m}} $。图2为第一帧原始方位距离图，对比可知LFM信号产生的混响更平稳，方位距离图中目标更清晰；对比图3和图4可知，PCA分解结果的稀疏杂点更多，目标轨迹周围存在更多的随机影响，不利于后续的目标检测，所以RPCA算法比PCA算法稀疏分离的效果更好。

在RPCA分解的基础上，取$ {d_{\exp {\mathrm{ect}}}} = 100\ {\mathrm{m}} $，$ {v_{\max }} = 10\ {\mathrm{m/s}} $，代入式(19)计算得到$ {\theta _{\max }} \approx {2^ \circ } $，则定位窗的大小设置为7×5，聚焦前一帧稀疏矩阵的稀疏点，扫描下一帧相同位置的定位窗，认为是目标的将其显示输出，否则就滤除，定位窗处理结果如图5所示。

对比图4和图5可知，定位窗处理后，稀疏矩阵的目标轨迹更加明显，目标周围区域的稀疏杂点被完全滤除，说明定位窗滤波效果很好，提高了动目标检测能力。

然而，从主动声呐显示界面的角度来说，实际检测过程判决门限依然不容易设定，因此进一步地从时频域的角度考量。将每个角度的时域信号分段做短时傅里叶变换，并按时间分段叠加，转换到时频域，则不同帧的同一角度所对应的时频信息具备低秩性。处理流程和结果将在后续试验数据处理中体现。

3.2 试验数据处理

试验在千岛湖进行，试验水域水深50 m，利用线列阵进行相控发射CW信号，信号脉宽200 ms，接收阵为线列阵，配合目标为一艘旅游船，最大吃水深度1.2 m，航速约10 kn，船长38 m。

图6是对第一帧数据作多普勒频移检测结果，中心是强混响的时频结果，其左侧可以隐约看到目标回波的位置。

对连续4帧处在77°方向的数据进行低秩稀疏矩阵分解，并将分解后的稀疏矩阵进行叠加输出。第一帧分解稀疏结果如图7所示，将所有帧稀疏矩阵叠加后结果如图8所示，对比图6可以看到，低秩稀疏矩阵分解后强混响（低秩矩阵）被滤除，但是除了目标回波以外，还有许多由噪声、随机起伏带来的杂点，影响了对目标的判断。

将前述定位窗目标单帧最大移动角度$ {\theta _{\max }} $等价替换为前后帧速度变化带来的多普勒频移变化$ \Delta {f_d} $，极限变化情况对应最大多普勒偏移：

$ \begin{aligned} \max (\Delta {f_d}) = &{f_d}{{\text{|}}_{v = {v_{\max }}}} - {f_d}{{\text{|}}_{v = - {v_{\max }}}}= \\ &\left(\frac{{c + {v_{{\mathrm{max}}}}}}{{c - {v_{\max }}}} - \frac{{c - {v_{{\mathrm{max}}}}}}{{c + {v_{\max }}}}\right){f_0} 。\end{aligned} $

(21)

当目标靠近时$ {v_{\max }} $取正值，远离时$ {v_{\max }} $取负值。经过计算分析可得，此处选用10×3的定位窗，聚焦第一帧矩阵的稀疏点，对后一帧稀疏矩阵中的稀疏点定位，重复前述操作直至最后一帧。定位窗处理结果如图9所示，对比图8可知，经过定位窗的筛选，声呐图像仅有连续的目标轨迹，而随机位置的杂点被滤除，显著改善了主动声呐的检测性能。

4 结　语

本文提出了基于低秩稀疏矩阵分解的目标定位窗滤波方法。先利用低秩稀疏矩阵分解算法分离具有低秩特性的混响和具有稀疏特性的目标，得到包含目标和稀疏杂点的稀疏矩阵，再根据目标运动具有的连续性特征和噪声、随机起伏等干扰具有的随机性特征差异，利用本文提出的定位窗方法对稀疏杂点进行滤除。仿真和试验数据处理结果表明，本算法对方位历程图和多普勒频移图均有着不错的滤波效果，其中在方位历程图中起到了突显目标轨迹的效果，在多普勒频移图中发挥了净化图像输出的作用，可显著提高主动声呐的自主检测性能，且算法原理简单，计算量小，易于工程实现。