混响抑制是水声信号处理的难题之一。声呐平台运动时，混响频谱会产生扩展，给滤波及目标检测增加了难度。运动平台声呐混响具有空时二维耦合特性，这决定了混响抑制属于二维滤波问题。空时自适应处理(space time adaptive processing, STAP)采用二维联合滤波自适应地感知干扰环境、抑制混响，是声呐抗混响检测的重要理论。

STAP是Brennan等^[1]首次提出的，并应用于机载相控阵雷达来补偿其中的运动效应。在水声信号处理领域，Klemm将其应用于运动声呐的低速目标检测，获得理想的抗混响性能^[2]。之后，人们在舰载声呐、鱼雷自导等水下目标检测中进行了较多的研究和实验^[3-6]。

在STAP中，待检测距离单元(rangecell under test, RUT)的干扰协方差矩阵估计是一个重要问题。Reed等在分析STAP收敛性时指出：有限样本的协方差矩阵估计会引起系统输出信干比损失，信干比损失不超过3 dB时，所需要的样本数应大于系统自由度的2倍^[7]。STAP系统自由度大，通常采用降维、降秩、模型参数化等次优处理以及知识辅助训练样本挑选等方法来减小系统对样本数的需求。

干扰样本获取也是一个重要问题。由于假设RUT位置还存在目标，传统STAP在独立同分布(independent identically distributed, IID)假设下，采用邻近RUT的辅助数据作为干扰样本来估计协方差矩阵。当海洋环境变化较大时，混响不满足IID，辅助数据的统计特性难以代表待检数据的统计特性，将降低STAP的性能。特别当RUT还存在目标时，STAP会将RUT中的目标当作干扰抑制掉。

稀疏重构(sparse recovery，SR)可利用极少观测样本高精度地恢复信号，为混响样本的高精度恢复提供了新的思路。

Mallat和Zhang首次提出在冗余字典进行稀疏分解(sparse decomposition)的思想，并提出匹配追踪(matching pursuit, MP)算法^[8]。之后，Pati提出正交匹配追踪(orthogonal matching pursuit, OMP)算法^[9]，改善了MP过匹配、收敛慢的问题，成为MP类算法的代表。1998年，Chen和Donoho等基于最优化理论，提出了基追踪(basis pursuit, BP)算法^[10]。通过研究l₁与l₀范数解的等效性，将一个多项式时间内无法求解的NP难问题转化为一个凸优化问题，开辟了稀疏重构的新途径。2006年，Candès及华裔科学家Tao等提出压缩感知(compressive sensing, CS)理论^[11-13]，从信号稀疏表示角度研究数据采集与恢复问题，发展了稀疏重构理论。

近年来，国内外研究者开始将稀疏分解理论用于STAP，以改善小样本条件下信号检测性能^[14-16]，研究工作主要集中在稀疏空时功率谱估计^[17]、利用辅助样本估计干扰协方差矩阵^[18-20]等方面。阳召成、段克清等对这部分工作进行了阶段性总结^[21-22]，分析了稀疏重构STAP(SR-STAP)对非平稳杂波的适应性、杂波样本获取、稀疏重构算法参数选择、多测量向量等应用中面临的关键问题，为深入开展研究提供了重要参考。

本文应用稀疏分解研究STAP中混响样本获取及协方差矩阵估计问题，所提出的方法可从待检测距离单元中重构混响样本，有效估计混响协方差矩阵，提高了STAP的抗混响性能。

1 混响的空时阵列信号模型

声呐阵列与待检距离环的几何关系如图 1所示。由图可以得到混响、回波多普勒频移与空域入射锥角余弦的耦合关系。入射锥角定义为阵列轴向与观测视线的夹角。

前视声呐的水声阵列排列在x轴上, 法线与声呐平台运动速度v一致。入射锥角余弦为cosα_i=sinθ_i·cosφ, 其中θ_i是混响单元或目标的方位角, φ是俯仰角。待检距离环上混响、回波的多普勒频移为

(1)

式中: λ是波长; v_i是混响单元或目标的速度。混响散射体的速度很小, 约等于0, f_di只与角度α_i有关。以不同锥角进入水声阵列的混响多普勒频移是不同的, 混响多普勒频移范围为[0, 2v/λ]。(1)式是圆方程, 在空时平面上, 混响分布在这个方程确定的半圆线上。运动目标回波的空时位置不在混响分布的半圆线上, 与混响是分离的。

侧视声呐的阵元排列在z轴上, 阵列法线与速度方向垂直。入射锥角余弦为cosα_i=cosθ_i·cosφ。混响、回波的多普勒频移为

(2)

多普勒频移与入射锥角余弦是线性关系, 混响分布在斜直线上, 多普勒频移范围为±2V/λ。同样, 运动目标回波的空时位置与混响是分离的。

按照混响单元散射模型，假设待检距离环上存在d个混响单元, 并考虑存在一个目标, 则当混响、回波进入到M个阵元、N次快拍的水声阵列时, 在一定假设条件下, 空时阵列信号模型表示为

(3)

式中, X是M×N维空时阵列采样矩阵, 其元素x_mn是第m个阵元、第n次快拍的采样数据, x是MN×1维空时采样向量, 为X的向量形式, vec(·)是向量化符号。A_i是幅度, 不妨记A_d+1是回波幅度, 其余为混响幅度。是相应的空时导向向量, 其中的是Kronecker直积运算符号, S_t(f_di)=[1, e^jω_ti, …, e^j(N-1)ω_ti]^T是频率导向向量, 其中ω_ti是时间角频率, 均匀采样时, 有ω_ti=2π(f₀+f_di)/f_s, f₀是发射信号频率, S_s(α_i)=[1, e^jω_si, …, e^j(M-1)ω_si]^T是方位导向向量, ω_si是空间角频率, 均匀线列阵的空间角频率与入射锥角有关, 为ω_si=2πΔ/λ·cosα_i, Δ是阵元间距。空时导向向量集合记为, 其中, 是多普勒频移、入射锥角组成的集合。

2 空时自适应处理分析 2.1 空时最优滤波器

STAP采用二维联合滤波自适应地感知干扰环境和抑制混响, 使输出信混比达到最大。Brennan根据似然比检测理论推导出空时最优滤波器的处理结构。

对于空时采样向量x, 设计MN×1维加权向量w作为空时滤波器, 它是多普勒频移、入射锥角的函数, 则空时滤波器输出为

(4)

若期望信号输出为y_s=A_sw^HS(f_d, α), 其中A_s是信号幅度, S(f_d, α)是其空时导向向量, 混响的输出功率为P=w^HR_cw, 其中R_c是混响协方差矩阵。在最小方差无失真响应(MVDR)准则下, 构造如下约束最优化问题

(5)

当R_c满秩时, 采用拉格朗日乘子法求解, 可得最优权向量为

(6)

式中, μ=1/[S^H(f_d, α)R_c^-1S(f_d, α)]是归一化空时导向向量。最优权向量与R_c有关, 依赖于输入混响的二阶统计性能。因此, w_opt是空时自适应权向量, 能够自适应地感知混响环境。

将最优权向量代入(4)式, 得空时最优滤波器的输出为

(7)

式中, μR_c^-1是对混响的白化处理, S(f_d, α)是对信号的广义匹配滤波。当混响是高斯白噪声n时, x=A_sS(f_d, α)+n。将其代入(7)式, 并由R_c=6n²I, 可得

(8)

容易看到, 最优权向量对期望信号进行了相位补偿, 实现了相干相加。噪声不能被补偿, 输出以功率相加。空时最优滤波器具有空间处理增益, 高斯噪声中的增益分贝数为10lg(MN)。

经过空时自适应处理, 就可直接在滤波器输出进行信号恒虚警检测。

2.2 混响样本获取

STAP的性能用改善因子(improvement factor, IF)来衡量。

IF定义为处理器输出信混比与输入信混比之比。设输入信混比I_SRR=A_s²/δ_c², 其中δ_c²是混响方差, 不难得到输出信混比为O_SRR=A_s²S^H(f_d, α)R_c^-1S(f_d, α)。因此, 空时最优滤波器的I_f为

(9)

I_f是多普勒频移和入射锥角的函数。通常使用其随多普勒频移变化的曲线表示处理器性能。

改善因子依赖于混响协方差矩阵R_c, 协方差矩阵估计对于I_f有很大影响。实际中, 混响协方差矩阵是由邻近RUT的辅助数据估计的, 如图 2所示。

雷达采用多脉冲体制时, 接收数据的组织形式是三维立方块, 如图 2a)所示。接收数据按照距离门方式排列, 相同距离门的数据堆砌成ML×1维观测样本向量, M是阵元数, L是脉冲数。

水下声速小, 传播慢, 因此声呐采用单脉冲体制, 1个探测周期只发射1个脉冲, 并完成数据采集和信号检测。采样数据的组织结构是二维形式, 如图 2b)所示。接收数据按照距离单元排列, M个阵元在在一个距离单元有L次快拍, 通常L等于发射脉宽τ。一个距离单元的数据堆砌成ML×1维观测样本向量。RUT中被假设为存在信号, 因此由邻近的辅助数据估计协方差矩阵。

不难看出, 单脉冲体制声呐能提供的样本数据少。不仅如此, 当海洋环境变化较大时, 辅助数据的统计特性难以代表RUT中混响的统计特性。当距离单元划分不合理导致辅助数据还存在目标时, STAP还会将RUT中的目标当作干扰抑制掉。因此, 声呐采用STAP时, 一般很难获得足够并满足统计性能要求的混响样本。

2.3 混响协方差矩阵估计的样本需求

混响协方差矩阵的最大似然(ML)估计为

(10)

式中: r_i是混响样本向量; K代表样本数。最大似然估计是大样本估计, K越大, 协方差矩阵估计越接近真值。以代替R_c, 则与之对应的最优权向量为

(11)

改善因子为

(12)

采用性能损失来衡量有限样本带来的影响, 将其定义为该改善因子与理想情况改善因子的比值, 可得

(13)

性能损失指标是的函数, 取值范围为。按照Reed等的研究结果^[7], 在K个样本服从IID条件下, 服从beta分布, 概率密度函数为

(14)

式中, L=MN是处理器自由度。性能损失指标的期望为

(15)

可见, 增大样本量可以减小性能损失。若使平均性能损失不超过3 dB, 即E(ρ)≥0.5, 则K≥2MN-3, 样本数约大于等于2倍系统自由度, 这是混响协方差矩阵估计所需要的最少样本数。不同自由度时的性能损失曲线如图 3所示。

降维、降秩等次最优处理就是通过减小系统自由度达到降低样本需求目的的。

3 混响样本的稀疏重构 3.1 空时导向字典设计

由模型(3)稀疏重构混响时, 需要密集划分空时平面来设计空时导向字典。

将多普勒频移、入射锥角划分成P×Q个密集网格, 则序号为(p, q)的网格对应的空时导向向量为。由空时导向向量构成冗余字典, 其中Ω= 是多普勒频移、入射锥角组成的集合, 则x的空时稀疏表示模型为

(16)

式中: c_{p, q}是稀疏分解系数, 表示网格上散射波的幅度; c是PQ×1维的稀疏分解系数向量。冗余字典是空时导向向量的扩展形式, 每一列S(f_dp, α_q)是一个原子, 共有PQ个原子。不考虑离网格(off-grid)的模型失配问题, 空时导向向量集合S(Θ)是冗余字典D(Ω)的子集。

按照模型(3), 待检距离环上存在d个混响单元和一个目标。当网格密集划分时, , 因此c中的元素大部分为0, 少量元素不为0, 混响是稀疏的。稀疏分解系数c_{p, q}与幅度A_i的关系为

(17)

不为0的元素等于混响或回波幅度。稀疏分解系数的平方就是空时二维谱估计。

注意到, 前视声呐的混响分布在一个半圆线上, 侧视声呐混响分布在一条斜直线上, 而运动目标回波的空时参数不在混响分布线上, 与混响位置是分离的。根据混响的空时分布特点, 若沿着混响分布线设计字典, 则该字典只包含混响单元的空时导向向量, 而不包含回波空时导向向量。将信号x在这个字典上进行稀疏分解, 则混响单元的分解系数就等于该混响幅度。因此, 由稀疏分解可重构RUT的混响样本。

为此, 将声呐工作扇面按照入射锥角等间隔划分为Q个角度网格。设第q个网格的方位锥角为α_q, 相应的多普勒频率为f_dq。在混响分布线上, f_dq是α_q的函数, 即f_dq=f_d(α_q), 因此, 沿着混响分布线设计的冗余字典为

(18)

式中, 是所有入射锥角组成的集合。冗余字典的每一列S(f_d(α_q), α_q)是一个原子, 共有Q个原子。将x在这个字典上做稀疏分解, 则空时稀疏表示模型为

(19)

式中, c_q是稀疏分解系数, c=[c₁, …, c_Q]^T是Q×1维向量。比较(19)式和(3)式可知, 在混响分布线上, 当离散网格上存在混响散射体时, 稀疏重构系数就等于混响幅度, 否则, 等于0。

由于待检距离环上只有d个混响单元, 而网格划分较密时, , 因此c中的元素大部分为0, 少量元素不为0, 混响是稀疏的。

数学上, 稀疏重构是一个约束最优化问题。对于空时稀疏表示模型(19), 考虑到噪声存在时, 稀疏表示的约束优化模型为

(20)

式中: ‖·‖₀是l₀范数, 表示c中非零元素的个数; ε是一个很小的量, 表示重构精度。因此, 稀疏重构就是利用极少观测样本高精度地恢复信号, 对于混响的精确表示与恢复具有优势。

由约束优化模型(20)求得系数向量c后, 稀疏重构的混响样本为

(21)

式中, r是MN×1维向量, 其中的元素r_i是混响瞬时幅度。

由于混响样本是由待检数据重构的, 而不依赖于临近的辅助数据, 基于这些样本可有效估计RUT的协方差矩阵。因此, 基于混响稀疏重构的STAP不仅适合于混响统计特性不变的情况, 也适合于统计特性变化的情况。不仅如此, 稀疏重构有降噪能力, 可实现混响样本的精确重构。

3.2 混响稀疏重构算法

求解约束优化模型(20), 可得混响幅度。但是l₀范数是非凸的, 采用遍历方法求是一个NP难问题。现在, 人们已提出许多优秀算法来求解这个问题。本文采用Pati的OMP算法, 这是一个简便、有效的算法。

OMP是迭代算法。每次迭代都是按最大内积准则在字典中寻找与上次迭代分解产生的残留误差最匹配的原子, 并对得到的原子进行正交化处理。

设R⁰x=x, 表示还未分解的残留误差, 则第m次迭代就是将残留误差R^mx对D的原子进行正交投影。假设第一次迭代时, 原子S₀是最匹配原子, 令u₀=S₀。第m次迭代时, R^mx选择的最匹配原子S_m为

(22)

式中: 0 < α≤1是优化因子; Γ_a是原子序号集合。对该原子做Gram-Schmidt正交化处理, 表示为

(23)

则有

(24)

式中, R^m+1x是第m次迭代后的残留误差。Mallat和Zhang证明了当每步中选择的原子满足最大内积条件时, 残留误差能量将按照指数规律衰减。

由正交化表示(23)式得

(25)

从而可得

(26)

此时, R^mx是在向量空间{S_k}_{0≤k < m}的补空间上投影, 投影得到的分解系数c_m=(R^mx)^TS_m/‖u_m‖²。由于R^mx与前面选择到的所有原子正交, {S_k}_{0≤k < m}之间是线性独立的, 因此可得

(27)

Pati已经证明, 采用OMP算法可使得残留误差在有限的迭代次数内减小到零^[9]。

4 混响协方差矩阵估计 4.1 混响的协方差矩阵估计

混响协方差矩阵估计是大样本估计, 按照性能损失分析, 所需样本数应大于2倍系统自由度。为此, 需要建立混响的概率分布模型, 产生足够的混响样本。

混响样本是网格上大量散射波的叠加。按照水声学原理, 可建立混响瞬时幅度的正态分布模型, 即r~N(μ, σ²), 其中μ和σ²分别是r的均值向量和方差向量, N(·, ·)表示正态分布。混响样本r是MN×1维随机向量, 当空时采样相互独立时, 联合概率密度为

(28)

式中, μ_i, σ_i²分别是r中元素的均值和方差。

产生混响样本时, 需要估计正态分布模型的参数。实际应用时, 可对阵列信号进行高速采样, 然后使用待检数据中的降采样数据来稀疏重构混响样本, 并采用这些重构样本估计正态模型参数。

假设由K组降采样数据重构的混响样本为r_i, 那么, 概率模型中均值和方差参数的最大似然估计为

(29)

(30)

有了混响的正态模型, 计算机仿真容易产生足够的高斯样本。

4.2 小样本时的对角加载

混响样本数较少时, 其协方差矩阵估计对于STAP性能有较大影响。采用对角加载可改善协方差矩阵估计的性能。

将的特征值记为, 相应特征向量为, 噪声的特征值记为, 则最优权向量可表示为

(31)

在样本数较多的理想情况下, 小特征值等于噪声功率σ_n², 但在有限样本时, 的小特征值在真实值σ_n²附近波动将使权向量发生变化, 影响处理性能。在实际情况中, 样本数受限时, 采用对角加载减小的波动。加载协方差矩阵估计为

(32)

式中, σ_L²为对角加载量。对角加载相当于增大了噪声, 对处理器的抗混响性能带来损失。一般情况下, 对角加载量范围为σ_n²≤σ_L²≤10σ_n²。

5 计算机仿真

分别采用声呐前视阵和侧视阵来验证本文方法的正确性, 并与采用辅助数据的传统空时最优滤波器进行性能比较。

仿真参数:设前视阵和侧视阵都是均匀线列阵, 阵元数为14, 半波长布阵。高频窄带脉冲宽度为5 ms。平台速度约为25.72 m/s。设待检数据包含有回波、混响和噪声, 而辅助数据只有混响和噪声。

按照单元散射模型对待检测距离单元和邻近距离单元的混响建模。将入射锥角等间隔划分, 间隔为5°, 共37个散射单元。混响幅度服从高斯分布, 信混比-10 dB, 信噪比20 dB。

目标在前视阵中方位角10°, 速度约为6.17 m/s, 目标在侧视阵中方位角40°, 速度约为6.17 m/s。

本文方法和传统空时最优滤波器的方位搜索步长1°, 多普勒搜索步长1 Hz。

仿真实验1:空时二维谱估计

图 4是空时二维谱估计比较。图 4a)至4b)分别是前视阵信号的傅氏空时谱估计和稀疏重构空时谱估计。图 4c)至4d)分别是侧视阵信号的傅氏空时谱估计和稀疏重构空时谱估计。冗余字典的网格间隔参数为1°×1 Hz。

在空时平面, 混响分布在多普勒频移与入射锥角余弦确定的曲线上, 前视阵的混响分布曲线是半圆, 侧视阵是斜直线。目标回波的空时位置与混响是分离的。傅氏空时谱估计的分辨率低, 不能分辨目标与混响分布线, 而稀疏重构空时谱估计的方位分辨率突破瑞利限, 频率分辨率突破傅立叶分辨率, 分辨性能优于傅氏空时谱, 清楚地分辨出了目标与混响分布线。

仿真实验2:空时自适应处理比较

图 5是前视阵列的空时自适应处理结果比较, 比较了本文方法与传统空时最优滤波器的性能。其中本文方法采用待检数据稀疏重构混响样本, 冗余字典的网格间隔参数为1°, 传统方法采用辅助数据作为混响样本。图 5a)是2类方法混响协方差矩阵估计的特征值比较, 协方差矩阵估计均采用了对角加载技术, 加载量取噪声功率。曲线明显下降之前的部分是混响特征值, 之后是噪声特征值。图 5b)是传统方法的空时响应。图 5c)是本文方法的空时响应。图 5d)是改善因子比较。

从特征值来看, 稀疏重构的降噪作用减小了混响样本中的噪声成分, 使得R_c中描述噪声的特征值很小, 这有助于减小噪声对STAP的影响。传统方法以辅助数据作为混响样本, 含有较大的噪声, R_c中的噪声特征值也较大。从空时响应看, 2类方法均在混响的半圆分布位置出现深凹槽, 表现出良好的抗混响性能。但本文方法的凹槽更深, 性能更优。

再从I_f来看, 改善因子在混响多普勒频移位置有一凹槽, 表示了方法的混响抑制性能, 凹槽越深, 抗混响性能越好。本文方法以较高的精度稀疏重构了混响样本, 其I_f凹槽平均值比传统方法降低约10 dB, 具有更强的抗混响能力。

图 6是侧视阵列的空时自适应处理结果比较，其中的混响样本获取方法、冗余字典网格间隔参数与前视阵列情况相同。图 6a)是混响协方差矩阵估计的特征值比较。图 6b)是传统方法的空时响应。图 6c)是本文方法的空时响应。图 6d)是改善因子比较。

同样，本文方法重构的混响样本中含有较小的噪声，描述噪声的特征值较小，空时响应在混响分布位置的凹槽更深，I_f凹槽平均值比传统方法降低约10 dB，表现出更强的抗混响能力。

6 结论

本文采用稀疏分解来重构STAP的混响样本，以较高的精度从待检数据中重构了混响样本，并基于混响模型的先验信息生成足够数量的样本来估计混响协方差矩阵，改善了STAP的抗混响性能。由于从待检数据来重构混响样本，而不是依赖于邻近RUT的辅助数据，在混响统计特性变化的海底环境中也具有较大应用前景。

声呐前视阵列和侧视阵列的计算机仿真结果验证了本文方法的性能。