Clutter Suppression Technology Based Space-time Adaptive ANM-ADMM-Net for Bistatic SAR
-
摘要: 双基合成孔径雷达(BiSAR)在实现对地面运动目标检测和成像时,需要抑制地面背景杂波。然而由于双基SAR收发分置的空间构型,会导致主瓣杂波出现严重的空时非平稳问题,从而恶化杂波抑制性能。基于稀疏恢复空时自适应处理方法(SR-STAP)虽然可以通过降低样本数量减少非平稳的影响,但是在处理过程中会出现字典离网问题,从而导致空时谱估计效果下降。并且大部分现有的典型SR-STAP方法虽然具有明确的数学关系和可解释性,但在针对复杂、多变场景时,也存在参数设置不恰当、运算复杂等问题。为解决上述一系列问题,该文提出了一种适用于双基SAR空时自适应杂波抑制处理的基于交替方向乘子法(ADMM)的复值神经网络ANM-ADMM-Net。首先,基于原子范数最小化(ANM)构建双基SAR连续空时域下杂波谱的稀疏恢复模型,克服传统离散字典模型下的离网问题;其次,采取ADMM对该双基SAR杂波谱稀疏恢复模型进行快速迭代求解;然后,根据迭代流程和数据流图进行网络化处理,将人工超参数迭代过程转换为网络可学习的ANM-ADMM-Net;再次,设置归一化均方根误差网络损失函数,并利用获取的数据集对网络模型进行训练;最后,利用训练后的ANM-ADMM-Net网络架构对双基SAR回波数据进行快速迭代处理,从而完成双基SAR杂波空时谱的精确估计和高效抑制。该文通过仿真试验和实测数据处理,表明该方法具有更好的杂波抑制性能和更加高效的运算效率。
-
关键词:
- 双基合成孔径雷达(BiSAR) /
- 稀疏恢复 /
- 空时处理 /
- 杂波抑制 /
- 复值神经网络
Abstract: Bistatic Synthetic Aperture Radar (BiSAR) needs to suppress ground background clutter when detecting and imaging ground moving targets. However, due to the spatial configuration of BiSAR, the clutter poses a serious space-time nonstationary problem, which deteriorates the clutter suppression performance. Although Space-Time Adaptive Processing based on Sparse Recovery (SR-STAP) can reduce the nonstationary problem by reducing the number of samples, the off-grid dictionary problem will occur during processing, resulting in a decrease in the space-time spectrum estimation effect. Although most of the typical SR-STAP methods have clear mathematical relations and interpretability, they also have some problems, such as improper parameter setting and complicated operation in complex and changeable scenes. To solve the aforementioned problems, a complex neural network based on the Alternating Direction Multiplier Method (ADMM), is proposed for BiSAR space-time adaptive clutter suppression. First, a sparse recovery model of the continuous clutter space-time domain of BiSAR is constructed based on the Atomic Norm Minimization (ANM) to overcome the off-grid problem associated with the traditional discrete dictionary model. Second, ADMM is used to rapidly and iteratively solve the BiSAR clutter spectral sparse recovery model. Third according to the iterative and data flow diagrams, the artificial hyperparameter iterative process is transformed into ANM-ADMM-Net. Then, the normalized root-mean-square-error network loss function is set up and the network model is trained with the obtained data set. Finally, the trained ANM-ADMM-Net architecture is used to quickly process BiSAR echo data, and the space-time spectrum of BiSAR clutter is accurately estimated and efficiently restrained. The effectiveness of this approach is validated through simulations and airborne BiSAR clutter suppression experiments. -
1. 引言
合成孔径雷达(Synthetic Aperture Radar, SAR)作为一种全天时、全天候的微波成像雷达,被广泛应用于环境探测和军用侦查等领域[1]。双基SAR (Bistatic Synthetic Aperture Radar, BiSAR)是收发装置分载于不同平台的雷达系统,相比单基SAR,具有抗干扰能力强,可实现前视成像等优势,是雷达领域研究热点之一[2]。随着需求的不断提升,双基SAR除了需要完成对地高分辨成像外,还需具备地面运动目标指示(Ground Moving Target Indication, GMTI)的能力,即实现双基SAR-GMTI功能。但是在双基SAR-GMTI中,强地面杂波常常会干扰运动目标检测,因此,杂波抑制是双基SAR-GMTI处理过程中的关键环节[3−7]。如何利用有效的杂波抑制方法去除背景杂波也成为目前学者关注的重点之一,并为此进行了大量的研究。
杂波抑制技术领域目前根据系统通道数量的差异,可以分为单通道抑制方法和多通道抑制方法,单通道杂波抑制方法只利用一个接收通道的回波进行处理,主要有基于多普勒滤波方法,基于调频率滤波方法等[8−10],其优点在于系统结构相对简单,但问题在于此类方法通常难以检测主瓣杂波内的慢速运动目标。多通道方法结合通道之间的空间信息差异,可以进一步提高杂波抑制效果。偏置相位中心天线技术(Displaced Phase Center Antenna, DPCA)是其中主流算法之一[11],其方法通过联合处理多个脉冲和多个天线相位中心的数据实现杂波对消,从而有效实现杂波抑制。而在双基SAR中,由于收发平台的分置,难以满足DPCA所要求的空时等效条件。文献[12]为此提出一种多脉冲移相中心天线方法,去除空时域之间的耦合关系,并克服DPCA的严格应用条件,解决双基SAR杂波的非平稳问题。除此之外,空时自适应处理方法(Space-Time Adaptive Processing, STAP)也是一种广泛应用的多通道杂波抑制技术,它利用杂波和运动目标在空时平面上的差异,实现杂波与动目标的分离,进而达到杂波有效抑制的目的[13]。根据RMB (Reed-Mallet-Brennan)准则[14],传统STAP方法在估计检测单元的杂波协方差矩阵时,需要的训练样本数至少满足系统自由度的两倍,并且样本之间的数据要求独立同分布。但是由于双基SAR收发平台的灵活构型配置和真实复杂环境的影响,杂波出现非平稳/非均匀特性[15],这将导致难以获得足够的有效训练样本,从而影响STAP方法对权矢量的估计,导致杂波抑制性能下降。为此,文献[16]提出一种基于级联抵消的非平稳杂波抑制方法,该方法首先在空时域利用STAP技术进行杂波抑制后,在图像域利用DPCA技术进一步抑制残余杂波。通过两步级联的方式缓解传统STAP单一方法的不足。
随着压缩感知与稀疏恢复理论的发展,基于稀疏恢复的空时自适应处理(STAP based on Sparse Recovery, SR-STAP)受到众多学者关注,该方法是利用杂波的空时稀疏特性,通过稀疏恢复算法估计杂波的空时功率谱,利用功率谱和杂波协方差矩阵的相关性重建杂波协方差矩阵,然后设计对应滤波器进行杂波抑制。其优点在于只需要少量训练样本便可以实现杂波协方差矩阵的估计,能够缓解杂波非平稳与非均匀问题[17−19]。传统SR-STAP方法中,通过对空时平面进行均匀离散化构建相应的超完备字典,结合稀疏恢复方法实现对杂波空时功率谱的高精度估计,进而设计相应的滤波器进行杂波抑制。但是不同于单基SAR系统中杂波的空时特性呈线性分布的特性,双基SAR系统杂波特性呈现出非线性高阶耦合的特征,这将导致字典原子无法与杂波脊进行精确匹配,即离网问题[20−22]。这将导致空时谱估计存在误差,造成杂波抑制效果下降。
针对双基SAR空时非线性造成的离网问题,也有学者做出相应研究,文献[23]提出基于降维局部搜索的OMP算法(Reduced-Dimension Local Search OMP, RD-LS-OMP),该方法是在传统正交匹配追踪(Orthogonal Matching Pursuit, OMP)算法的基础上[24],对匹配到的原子再次进行局部网格化,在局部网格字典中再次进行原子匹配。此方法进一步提高了原子匹配精度,但是局部网格字典仍是均匀离散的,在面对双基SAR的非线性杂波特性上,依然会出现原子匹配不准确的问题,并且局部字典的增加,也会出现较多冗余的计算。文献[25]提出了一种局部网格分割算法,通过子空间投影选择残差最相关的空时向量原子,最后围绕选定的原子进行局部网格分割,从而获得较高的匹配精度。文献[26]则提出一种基于杂波脊匹配的非平稳杂波抑制方法,利用先验知识构造尽可能与杂波空时特性匹配的空时字典,从而缓解离网问题,但是该方法需要精确的先验信息构造匹配空时字典,在先验信息有误的情况下,仍会导致杂波抑制效果下降。总体而言,此类基于字典网格的SR-STAP算法不可避免会出现字典原子失配问题,从而造成空时谱出现估计误差。为此文献[27]提出基于原子范数最小化(Atomic Norm Minimization, ANM)的稀疏恢复方法,该方法在连续空时平面内,即无限精度内的原子集合中寻求匹配原子,消除了离网问题,具有非常高的精确度,但是相应的问题在于运算量过大,导致效率低下,难以进行实际应用。文献[28]基于交替方向乘子法(Alternating Direction Method of Multipliers, ADMM)提出相应的ANM-ADMM方法提高ANM-STAP的计算效率,但是存在ADMM超参数需要人工设置的问题,当设置不恰当的超参数时,会严重影响ADMM算法的收敛速度和精度,从而恶化后续杂波抑制效果。
针对上述这类需要设置参数的模型驱动方法所面临的问题,深度学习相关的研究给予了专家学者启发。传统迭代算法通常通过分析物理过程进行模型设计,而深度学习方法则试图从数据中进行学习,优化网络模型参数。而展开技术(Unrolling)为“模型驱动方法”与“数据”结合提供了有效桥梁。它可以将一个算法模型的迭代优化转换为一个动态网络系统,将迭代过程离散化为若干可学习模块,进而通过数据驱动进行参数演化,最终形成性能更加优异、可解释性更强的深度网络模型[29,30]。利用该技术的特点,本文针对ANM-ADMM算法中存在的超参数设置问题,结合深度展开技术,构建了一种基于ADMM的复值神经网络ANM-ADMM-Net,该方法在利用ADMM算法提高计算效率的同时,通过网络方法避免了参数设置不恰当的问题,从而降低人为设置参数导致的性能恶化,可以有效提高算法性能,并利用较少网络层数,进一步提高算法运算效率。首先构建对应双基SAR回波模型,利用ADMM算法对ANM-STAP算法进行迭代估计;其次,根据迭代流程和数据流图对每一步迭代过程进行网络化处理,从而构建将人工超参数转换为网络可学习参数的ANM-ADMM-Net;然后,设置对应的网络损失函数,并利用相应数据集对网络模型进行训练,获取最优网络架构;最后,利用训练后网络架构对训练距离单元数据进行处理,获取对应杂波空时谱的准确估计,设计对应的空时滤波器进行杂波抑制。仿真结果表明,相比于其他SR-STAP方法,本文所提算法利用数据驱动的方式避免了传统参数的设置不当问题,在具有更高估计效果的同时,进一步降低了运算复杂度。
2. 双基SAR信号模型及SR-STAP方法
2.1 信号模型
双基SAR几何构型如图1所示,设发射机初始坐标为(xT,yT,zT),接收机的初始坐标为(xR,yR,zR)。接收站沿航迹方向设置N元均匀线阵,通道间隔为d。P(x,y)是地面目标点,VT为发射机速度,VR为接收机速度;δT与δR分别为发射机飞行方向和接收机飞行方向,θT与φT分别是发射机相对于目标点P的方位角和俯仰角,θR与φR分别是接收机相对于目标点P的方位角和俯仰角,vx与vy分别是运动目标点P在x轴和y轴的速度分量。L是发射机与接收机之间的基线在地面的投影距离。
其中,RsT与RsR分别是目标点距离发射机与接收机的距离,对应的接收机的回波双基距离和可以表示为
Rsn(x,y;t)=RsT(x,y;t)+RsRn(x,y;t) (1) 其中,t表示慢时间,RsT(x,y;t)是发射机到目标点P的距离历史,RsRn(x,y;t)是接收机第n个通道到目标点P的距离历史,具体可以表示为
RsT(x,y;t)=√(xT+(VTx−vx)t−x)2+(yT+(VTy−vy)t−y)2+z2T (2) RsRn(x,y;t)=√(xR+(VRx−vx)t−x)2+(yR+(n−1)d+(VRy−vy)t−y)2+z2R (3) 其中,VTx与VTy是发射机在x, y轴的速度分量,VRx, VRy是接收机在x, y轴的速度分量,d是接收通道之间的间距,vx与vy分别表示目标点对应的x轴、y轴的速度分量。
在已知目标点距离双基距离和后,可得对应回波信号模型,具体公式如下:
Sn(t,τ)=ωr(τ−Δτ)ωa(t)exp{jπBT(τ−Δτn)2}⋅exp(−j2πfcΔτn) (4) 其中,t与τ分别是方位向时间变量和距离向时间变量,Δτn=Rsn/Rsncc是双基距离和引起的时延,n表示第n个阵元通道,ωr(⋅)表示距离向窗函数,ωa(⋅)表示方位向窗函数,B/BTT是LFM信号的调频率,其中B是信号带宽,T是信号时宽。
对回波数据进行距离向脉冲压缩及快速傅里叶变换(Fast Fourier Transform, FFT)处理,可得对应回波为
Sn(t,fτ) = ωr(fτ)ωa(t)exp{−j2πfτ+fccRsn} (5) 其中,fτ是FFT处理后的频域变量,fc是载频频率。
可将式(1)进行泰勒公式展开,可得
Rsn(x,y;t)=Rsn(x,y;0)+R′sn(x,y;0)t+12R″sn(x,y;0)t2+⋯ (6) 将式(6)代入式(5)可得
Sn(t,fτ)=ωr(fτ)ωa(t)⋅exp{−j2πfτc(Rsn(x,y;0)+R′sn(x,y;0)t+12R″sn(x,y;0)t2+⋯)}⋅exp{−j2πfcc(Rsn(x,y;0)+R′sn(x,y;0)t+12R″sn(x,y;0)t2+⋯)} (7) 由式(7)可以发现,在第1个指数项中,回波信号的频率变量fτ与时间t耦合,该耦合导致回波数据出现距离徙动问题,为此需要对回波信号进行距离徙动校正[31−32],最终得到对应基带回波信号:
Sn(t,τ)=sinc{B(τ−Rsnc)}ωa(t)exp{−j2πfccRsn} (8) 2.2 双基SAR杂波非平稳特性分析
双基SAR回波信号的多普勒频率和空间频率具有二维耦合特性,其中,双基机载雷达的多普勒频率和空域频率分别可以表示为
fd=VRλfrcos(θR−δR)cosφR+VTλfrcos(θT−δT)cosφT (9) fs=dλcos(θR−δR)cosφR (10) 其中,λ与fr分别为信号波长与脉冲重复频率(Pulse Repetition Frequency, PRF)。双基SAR由于收发分置,多普勒频率具有收发平台的两个分量,空时频率不再是线性分布,并且不同距离单元下具有不同的分布关系。通过仿真验证进一步分析双基SAR杂波空时特性,设置仿真参数:发射机方位与速度分别为(−4,−3.5,4)km和(0,130,0)m/s,接收机方位与速度分别为(0,−4,4)km和(0,130,0)m/s,脉冲重复频率为fr=1000Hz,信号波长为λ=0.03m。
图2中不同颜色代表不同距离单元,在单基正侧视构型下,不同距离单元具有相同的空时分布曲线,杂波特性不依赖于距离分布。在双基平飞构型下,不同距离单元下双基SAR的杂波脊空时特性具有不同的分布,并且呈曲线形状。这种情况下,待检测距离单元样本与邻近单元样本空时特性不再满足独立同分布特性,利用2NM个临近样本估计得到的空时谱会出现严重展宽。如图3所示,图3(a)为理想的待检测单元空时谱,图3(b)为传统估计所得待检测单元空时谱,对应杂波空时谱展宽从而淹没动目标,影响后续目标检测与估计。
2.3 SR-STAP原理及离网问题
不同于传统STAP方法需要利用2NM个单元样本去估计待检测单元的协方差矩阵,基于稀疏恢复的空时自适应处理方法结合回波信号在空时谱中的稀疏特性,利用少量样本进行空时稀疏恢复求解,估计得到杂波的高精度空时谱。在传统SR-STAP[21−26]过程中,将空时平面进行归一化离散处理,杂波信号可以表示为
y=Φα+n (11) 其中,y∈CNM×1是某一距离单元内的杂波数据,Φ∈CNM×NsMd是将空时平面均匀离散处理后形成的空时过完备字典,α∈CNsMd×1为对应字典原子的幅值,n∈CNM×1为噪声。其中Φ划分大小为Ns=ρsN,Md=ρdM,ρs, ρd (ρs>1, ρd>1)是离散量化因子,字典Φ中每个原子都代表特定的空间频率和多普勒频率,即
Φ=[s(f1d,f1s),s(f1d,f2s),⋯,s(f1d,fMds),⋯,s(fid,fjs),⋯,s(fNsd,fMds)] (12) 其中,s(fid,fjs)=s(fid)⊗s(fjs)为对应空时二维导向矢量,⊗表示Kronecker积,s(fid)与s(fjs)分别表示杂波的时间导向矢量和空间导向矢量,具体表示如下:
{s(fid)=[1,ej2πfid,ej2π⋅2fid,⋯,ej2π⋅(M−1)fid]Ts(fjs)=[1,ej2πfjs,ej2π⋅2fjs,⋯,ej2π⋅(N−1)fjs]T (13) 式(11)完成杂波特性的稀疏表示后,通过稀疏求解得到对应的稀疏解α,从而恢复出高精度空时谱。
一般稀疏恢复会采用待检测单元周围若干样本,从而构造如下的约束优化问题进行求解:
ˆa=argmin (14) 其中, {\boldsymbol{Y}} = [{{\boldsymbol{y}}_1},{{\boldsymbol{y}}_2}, \cdots ,{{\boldsymbol{y}}_L}] \in {\mathbb{C}^{NM \times L}} 为L个距离单元的杂波样本,|| \cdot |{|_{\text{F}}}表示Frobenius范数,{\boldsymbol{A}} = [{{\boldsymbol{\alpha}} _1}, {{\boldsymbol{\alpha}} _2}, \cdots ,{{\boldsymbol{\alpha}} _L}] \in {\mathbb{C}^{NM \times L}}分别为每个距离单元对应的稀疏解, || \cdot |{|_{2,0}} 表示先对稀疏矩阵解集各行计算{l_2}范数,然后对列求{l_0}范数。
通过相应稀疏恢复算法,可以通过求解对应稀疏解,得到空时谱中对应空时频率的幅值, {a_i} 表示每个空时频率点上的杂波幅值:
{\boldsymbol{\alpha}} = {[{a_1},{a_2}, \cdots, {a_q}, \cdots ,{a_{{N_{\text{s}}}{M_{\text{d}}}}}]^{\text{T}}} (15) 根据字典与稀疏解,进行杂波协方差估计:
{{\boldsymbol{R}}_{\text{c}}} = E \left[{{\boldsymbol{y}}_l}{\boldsymbol{y}}_l^{\text{H}}\right] = {{\boldsymbol{\varPhi}} ^{\text{H}}}{\text{diag}}({\boldsymbol{P}}){\boldsymbol{\varPhi}} (16) 其中,{\boldsymbol{P}} = [a_1^2,a_2^2,\cdots,a_{{N_{\text{s}}}{M_{\text{d}}}}^2]^{\mathrm{T}}表示对应空时平面的杂波功率谱,{\text{diag}}( \cdot )表示矩阵对角化。
最终得到对应的空时滤波器权矢量:
{{\boldsymbol{W}}_{{\text{opt}}}} = \mu {\boldsymbol{R}}_{\text{c}}^{ - 1}{\boldsymbol{s}} (17) 其中,s是目标对应的空时导向矢量,\mu = {1/ {({{\boldsymbol{s}}^{\text{H}}}{\boldsymbol{R}}_{\text{c}}^{ - 1}{\boldsymbol{s}})}}为常数项。
在单基正侧视场景下,空时谱通常表现为一条倾斜直线,在稀疏恢复过程中,可以较为准确的定位在空时字典中,如图4(a)所示。
然而,在双基SAR构型下,空时杂波脊不再是一个均匀的直线分布。如图4(b)所示,杂波脊整体呈现出曲线形状。在稀疏恢复过程中,杂波能量会在求解过程中被分散到临近平面内的若干网格点中,从而导致能量泄露。最终表现为空时谱展宽,降低SR-STAP的杂波抑制性能。虽然通过调节网格密度,可以一定程度缓解逸散情况,但是会增加网络原子之间的相关性,造成恢复效果的下降。
3. ANM-ADMM-Net方法
针对离网问题,文献[27]提出一种基于原子最小范数的SR-STAP方法,但是该方法存在计算复杂度高的缺点,针对该类优化问题,可以利用ADMM算法[28]加快运算效率,但是ANM-ADMM仍存在超参数需人工设置的缺陷,本文基于展开技术进一步将对应ANM-ADMM进行网络化处理,利用数据驱动的方式解决人工参数设置的问题,并基于网络的简洁架构,减少迭代次数,进一步加快运算效率,我们首先导出ANM-ADMM算法[27,28],然后进一步在此基础上建立ANM-ADMM-Net网络架构。
3.1 ANM-ADMM算法
为避免离网问题,可以将式(14)中对应的优化问题转换为原子范数最小问题。为此首先需要将离散字典转换为空时导向矢量连续的连续原子集 \mathcal{A} ,具体表示为
\begin{split} \mathcal{A} \triangleq \,&\left\{ {\boldsymbol{s}}\left( {{f_{\text{s}}},{f_{\text{d}}}} \right) = {\boldsymbol{s}}\left( {{f_{\text{s}}}} \right) \otimes {\boldsymbol{s}}\left( {{f_{\text{d}}}} \right)\left| {f_{\text{s}}} \in \left( { - 1,1} \right), \right.\right.\\ & \left. {f_{\text{d}}} \in \left( { - 1,1} \right) \right\} \end{split} (18) 利用杂波在空时谱中的稀疏特性,可以定义杂波信号{{\boldsymbol{y}}_l}的原子范数如下:
{\left\| {{{\boldsymbol{y}}_l}} \right\|_{\mathcal{A},0}} \triangleq \mathop {\inf }\limits_{\substack{ {\alpha _i} \in \mathbb{C} \\ f_{\text{s}}^i \in \left( { - 1,1} \right),f_{\text{d}}^i \in \left( { - 1,1} \right) }} \sum\limits_{i = 1}^O {{\alpha _i}{\boldsymbol{s}}\left( {f_{\text{s}}^i,f_{\text{d}}^i} \right)} (19) 其中, {\boldsymbol{s}}\left( {f_{\text{s}}^i,f_{\text{d}}^i} \right) 为来自字典 \mathcal{A} 的空时导向矢量, {\alpha _i} 为其对应的幅度。通过求解如下的原子范数最小化问题得到对应的稀疏恢复解:
\min {\left\| {\boldsymbol{y}} \right\|_{\mathcal{A},0}}\;{\mathrm{s.t}}.\left\| {{{\boldsymbol{y}}_l} - {\boldsymbol{y}}} \right\|_2^2 < {\varepsilon _n} (20) 利用杂波原子范数与其协方差矩阵的秩大小相同,即 {\left\| {{{\boldsymbol{y}}_l}} \right\|_{\mathcal{A},0}} = {\text{rank}}({{\boldsymbol{R}}_l}) = {M_{\text{r}}} ,以及杂波协方差矩阵{{\boldsymbol{R}}_l}的块Toeplitz矩阵特性,通过遵循ANM问题求解思路,可以将式(20)优化问题最终转换如下形式:
\begin{split} & \mathrm{min}\;m+\text{tr}\left({\boldsymbol{R}}\right)\\ & {\mathrm{s.t}}.\;{\boldsymbol{T}}=\left[\begin{array}{cc}m& {{\boldsymbol{y}}}^{\text{H}}\\ {\boldsymbol{y}}& {\boldsymbol{R}}\end{array}\right]\ge 0,{\boldsymbol{R}}={S}\left({\boldsymbol{U}}\right),\;{\Vert {{\boldsymbol{y}}}_{l}-{\boldsymbol{y}}\Vert }_{2}^{2} < {\varepsilon }_{n} \end{split} (21) 其中,m为引入的未知参数,表示未知的杂波协方差矩阵的秩。{\text{tr}}( \cdot )表示矩阵的迹。y为输入的杂波数据,R表示对应的杂波协方差矩阵,并且是一个半正定的块Toeplitz矩阵,可以表示为 {\boldsymbol{R}} = {{S}}\left( {\boldsymbol{U}} \right) ,是一个由大小为 \left( {2N - 1} \right)\left( {2M - 1} \right) 的复数矩阵U生成的N \times N块Toeplitz矩阵,其形式为
{{S}}\left( {\boldsymbol{U}} \right) = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{U}}_0}}&{{{\boldsymbol{U}}_1}}& \cdots &{{{\boldsymbol{U}}_M}} \\ {{{\boldsymbol{U}}_{ - 1}}}&{{{\boldsymbol{U}}_0}}& \cdots &{{{\boldsymbol{U}}_{M - 1}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{{\boldsymbol{U}}_{ - M}}}&{{{\boldsymbol{U}}_{1 - M}}}& \cdots &{{{\boldsymbol{U}}_0}} \end{array}} \right] (22) 其中,{{\boldsymbol{U}}_n}是为M \times M的Toeplitz矩阵,由矩阵U的第n行生成,具体定义为
{{\boldsymbol{U}}_n} = \left[ {\begin{array}{*{20}{c}} {{u_0}}&{{u_1}}& \cdots &{{u_N}} \\ {{u_{ - 1}}}&{{u_0}}& \cdots &{{u_{N - 1}}} \\ \vdots & \vdots & \ddots & \vdots \\ {{u_{ - N}}}&{{u_{1 - N}}}& \cdots &{{u_0}} \end{array}} \right] (23) 当采用多个距离单元样本进行联合稀疏恢复时,式(21)对应的优化问题可以扩展如下形式:
\begin{split} & \mathrm{min}\,\text{tr}({\boldsymbol{Z}})+\text{tr}\left({\boldsymbol{R}}\right)\\ & {\mathrm{s.t}}.\;{\boldsymbol{T}}=\left[\begin{array}{cc}{\boldsymbol{Z}}& {{\boldsymbol{Y}}}^{{\mathrm{H}}}\\ {\boldsymbol{Y}}& L{\boldsymbol{R}}\end{array}\right]\ge 0,\;{\boldsymbol{Z}}={{\boldsymbol{Z}}}^{\text{H}},\\ & {\boldsymbol{R}}={S}\left({\boldsymbol{U}}\right),\;\;{\Vert {{\boldsymbol{Y}}}_{l}-{\boldsymbol{Y}}\Vert }_{2}^{2} < {\varepsilon }_{n} \end{split} (24) 其中, {\boldsymbol{Y}} = [{{\boldsymbol{y}}_1},{{\boldsymbol{y}}_2}, \cdots ,{{\boldsymbol{y}}_L}] \in {\mathbb{C}^{NM \times L}} ,{\boldsymbol{ Z}} \in {\mathbb{C}^{L \times L}} 是Hermitian矩阵,L为距离单元样本数。通过求解上述优化问题后,所得R不是对杂波协方差矩阵的准确估计,而是对应的杂波子空间估计,在求得对应R后,最终杂波协方差矩阵可以被估计为
{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{R}}} _{{\text{ANM}}}} = \frac{1}{L}{\boldsymbol{\varGamma}} {\text{diag}}\left({\left| {{{\boldsymbol{\varGamma}} ^{ - 1}}{\boldsymbol{y}}} \right|^2}\right){{\boldsymbol{\varGamma}} ^{\text{H}}} + \sigma _n^2{{\boldsymbol{I}}_{NM}} (25) 其中,{\boldsymbol{\varGamma}} 是估计杂波协方差矩阵进行特征值分解 {\boldsymbol{R}} = {\boldsymbol{\varGamma}} {\boldsymbol{\varSigma}} {{\boldsymbol{\varGamma}} ^{ - 1}} 后的特征向量矩阵,{\boldsymbol{\varSigma}} 为特征值组成的对角矩阵, \sigma _n^2 为噪声方差, {{\boldsymbol{I}}_{NM}} 表示为NM阶单位矩阵。
在此基础上,利用ADMM算法对式(24)的优化问题进行求解,可改写为如下所示的增广拉格朗日形式:
\begin{split} & {\mathcal{L}}({\boldsymbol{Z}},{{\boldsymbol{Y}}}_{\text{C}},{\boldsymbol{U}},{\boldsymbol{T}},{\boldsymbol{\varLambda}} )\\ & \quad =\frac{\lambda }{2}\Vert {\boldsymbol{Y}}-{{\boldsymbol{Y}}}_{\text{C}}{\Vert }_{\text{F}}^{2}+\text{Tr}(\mathcal{S}({\boldsymbol{U}}))+\text{Tr}({\boldsymbol{Z}})\\ &\qquad + \left\langle {\boldsymbol{\varLambda}} ,{\boldsymbol{T}}-\left[\begin{array}{cc}{\boldsymbol{Z}}& {{\boldsymbol{Y}}}_{\text{C}}^{\text{H}}\\ {{\boldsymbol{Y}}}_{\text{C}}& \mathcal{S}(U)\end{array}\right] \right\rangle\\ & \qquad+\frac{\rho }{2}{\left\| {\boldsymbol{T}}-\left[\begin{array}{cc}{\boldsymbol{Z}}& {{\boldsymbol{Y}}}_{\text{C}}^{\text{H}}\\ {{\boldsymbol{Y}}}_{\text{C}}& \mathcal{S}({\boldsymbol{U}})\end{array}\right]\right\| }_{\text{F}}^{2} \end{split} (26) 其中,\rho > 0是惩罚因子,\lambda > 0是归一化参数,使用 \mathcal{S}({\boldsymbol{U}}) 的形式代表协方差矩阵R。T与 {\boldsymbol{\varLambda}} 均为半正定Hermitian矩阵, {\boldsymbol{\varLambda}} 拉格朗日乘子矩阵,均可以表示为
\left. {{{\boldsymbol{T}}^t} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{T}}_{\boldsymbol{Z}}^t}&{{{({\boldsymbol{T}}_{{{\boldsymbol{Y}}_{\text{C}}}}^t)}^{\text{H}}}} \\ {{\boldsymbol{T}}_{{{\boldsymbol{Y}}_{\text{C}}}}^t}&{{\boldsymbol{T}}_{\mathcal{S}({\boldsymbol{U}})}^t} \end{array}} \right.} \right] (27) \left. {{{\boldsymbol{\varLambda}} ^t} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{\varLambda}} _Z^t}&{{{({\boldsymbol{\varLambda}} _{{{\boldsymbol{Y}}_{\text{C}}}}^t)}^{\text{H}}}} \\ {{\boldsymbol{\varLambda}} _{{{\boldsymbol{Y}}_{\text{C}}}}^t}&{{\boldsymbol{\varLambda}} _{\mathcal{S}(U)}^t} \end{array}} \right.} \right] (28) 通过ADMM求解上述优化问题,则得到对应第t+1次迭代时对应的更新步骤为
\left\{ \begin{gathered} {\boldsymbol{Y}}_{\mathrm{C}}^{t + 1} = \frac{1}{{2\rho + \lambda }}\left(\lambda {\boldsymbol{Y}} + 2{\boldsymbol{\varLambda}} _{{{\boldsymbol{Y}}_{\text{C}}}}^t + 2\rho {\boldsymbol{T}}_{{{\boldsymbol{Y}}_{\text{C}}}}^t \right) \\ {{\boldsymbol{Z}}^{t + 1}} = {\rho ^{ - 1}}{\boldsymbol{\varLambda}} _{\boldsymbol{Z}}^t + {\boldsymbol{T}}_{\boldsymbol{Z}}^t - {\rho ^{ - 1}}{{\boldsymbol{I}}_L} \\ {{\boldsymbol{U}}^{t + 1}} = {\mathcal{S}^*}\left({\rho ^{ - 1}}{\boldsymbol{\varLambda}} _{\mathcal{S}({\boldsymbol{U}})}^t + {\boldsymbol{T}}_{\mathcal{S}({\boldsymbol{U}})}^t \right) - {\rho ^{ - 1}}{{\boldsymbol{E}}_{M,N}} \\ {\boldsymbol{T}}_{{\text{un}}}^{t + 1} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{Z}}^{t + 1}}}&{{{({\boldsymbol{Y}}_{\text{C}}^{t + 1})}^{\text{H}}}} \\ {{\boldsymbol{Y}}_{\text{C}}^{t + 1}}&{\mathcal{S}({{\boldsymbol{U}}^{t + 1}})} \end{array}} \right] - {\rho ^{ - 1}} {{\boldsymbol{\varLambda}} ^t} \\ {{\boldsymbol{T}}^{t + 1}} = {\boldsymbol{G}}{\text{diag}}({\{ {\delta _g}\} _ + }){{\boldsymbol{G}}^{ - 1}},\;g = 1,2, \cdots ,MN + L \\ {{\boldsymbol{\varLambda}} ^{t + 1}} = {{\boldsymbol{\varLambda}} ^t} + \rho \left( {{{\boldsymbol{T}}^{t + 1}} - \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{Z}}^{t + 1}}}&{{{({\boldsymbol{Y}}_{\text{C}}^{t + 1})}^{\text{H}}}} \\ {{\boldsymbol{Y}}_{\text{C}}^{t + 1}}&{\mathcal{S}({{\boldsymbol{U}}^{t + 1}})} \end{array}} \right]} \right) \\ \end{gathered} \right. (29) 其中,上标t表示迭代次数,\mathcal{S}({\boldsymbol{U}})表示将参数矩阵U转化为对应的块Toeplitz矩阵。 {\mathcal{S}^*}({\rho ^{ - 1}}{\boldsymbol{\varLambda}} _{\mathcal{S}({\boldsymbol{U}})}^t + {\boldsymbol{T}}_{\mathcal{S}({\boldsymbol{U}})}^t) 表示对 {\rho ^{ - 1}}{\boldsymbol{\varLambda}} _{\mathcal{S}({\boldsymbol{U}})}^t + {\boldsymbol{T}}_{\mathcal{S}({\boldsymbol{U}})}^t 进行逆Toeplitz化变换[28],其中对{{\boldsymbol{I}}_{MN}}进行逆Toeplitz化可得 {{\boldsymbol{E}}_{M,N}} , {{\boldsymbol{E}}_{M,N}} 是一个大小为(2M - 1)(2N - 1)的矩阵,第(M,N)个元素为1,其他元素均为0。G表示{\boldsymbol{T}}_{{\text{un}}}^{t + 1}的特征向量矩阵,{\delta _g}是{\boldsymbol{T}}_{{\text{un}}}^{t + 1}的特征值,{\{ {\delta _g}\} _ + }表示保留非负特征值,即将负特征值置零。在上述求解过程中,需要设置恰当的超参数\rho 和\lambda ,通过迭代求解获取最终的协方差矩阵。
在实际应用中,首先需要不断人为调整合适的超参数,并且上述算法仍需要高维的迭代次数,每次迭代过程中涉及矩阵的特征分解和块Toeplitz变换,使得估计杂波协方差矩阵耗时严重,难以满足STAP实时处理要求,并且难以设置最为合适的超参数,导致算法收敛速度下降和估计精度不足。为解决此问题,本文针对上述方法进行网络化处理,利用数据驱动的方式获取最优参数,并提高运行效率。
3.2 ANM-ADMM-Net
针对式(26)所表示的优化问题,ANM-ADMM方法在求解过程中需要设置超参数,当参数设置不恰当时,会恶化算法的求解效率和求解结果,为此将ANM-ADMM进行网络化展开,构建ANM-ADMM-Net,利用目标场景对应数据集进行网络参数训练,从而获取对应最优参数。并且基于网络的简洁架构加快运算效率,利于工程化实时处理。
根据式(29)对应的算法步骤,可以将对应过程转换为数据流图的形式,如图5所示每一块代表一层网络结构,对应ADMM的一次完整迭代过程,数据通过箭头指示传递,每层网络都有4个子结构,依次对应ADMM算法每次迭代过程中的更新步骤,分别为数据重构层 \mathcal{X} ,非线性变换层 \mathcal{S} ,特征分解层\mathcal{R},乘子更新层\mathcal{M}。其中作为网络与外部数据进行交互的初始输入和最终输出结构,与其他中间结构有所不同,这里表示为 {\text{net\_org}} 阶段与 {\text{net\_end}} 输出阶段。
网络的输入输出架构和各层子结构具体如下:
(1) 初始输入 {\text{net\_org}} :
如图5所示,将ANM-ADMM算法进行数据流图等效后,数据将以矩阵的形式进行传递,初始输入的作用是将ADMM中的矩阵和拉格朗日乘子矩阵进行初始化设置,即将式(29)中的T矩阵与拉格朗日乘子矩阵{\boldsymbol{\varLambda}} 初始化为全0矩阵,也分别对应初始阶段的特征分解层{\mathcal{R}^0}和拉格朗日乘子更新层{\mathcal{M}^0},并将杂波数据{\boldsymbol{Y}} \in {\mathbb{C}^{NM \times L}}传递至数据重构层。
(2) 数据重构层\mathcal{X}:
数据重构层包含了ADMM中的3个步骤,作用是将杂波数据及对应乘子参数等整合为一整体矩阵并进行更新,即对应ADMM中的 {\boldsymbol{Y}}_{\text{C}}^{t + 1} , {{\boldsymbol{Z}}^{t + 1}} 以及 {{\boldsymbol{U}}^{t + 1}} 中的非线性部分更新步骤。在该层主要表现为对杂波数据Y,秩矩阵Z的迭代更新,其中输入原始杂波数据{\boldsymbol{Y}} \in {\mathbb{C}^{NM \times L}}, {{\boldsymbol{T}}^t}, {{\boldsymbol{\varLambda}} ^t} 以及可学习参数\rho 和\lambda 。将其中每个矩阵对应的数据部分进行提取,根据式(30)—式(32)进行更新:
\quad {\boldsymbol{Y}}_{\text{C}}^{t + 1} = \frac{1}{{2\rho + \lambda }}\left(\lambda {\boldsymbol{Y}} + 2{\boldsymbol{\varLambda}} _{{{\boldsymbol{Y}}_{\text{C}}}}^t + 2\rho {\boldsymbol{T}}_{{Y_{\text{C}}}}^t \right) (30) \quad {{\boldsymbol{Z}}^{t + 1}} = {\rho ^{ - 1}}{\boldsymbol{\varLambda}} _{\boldsymbol{Z}}^t + {\boldsymbol{T}}_{\boldsymbol{Z}}^t - {\rho ^{ - 1}}{{\boldsymbol{I}}_L} (31) \quad {{\boldsymbol{U}}^{*t + 1}} = \left({\rho ^{ - 1}}{\boldsymbol{\varLambda}} _{\mathcal{S}({\boldsymbol{U}})}^t + {\boldsymbol{T}}_{\mathcal{S}({\boldsymbol{U}})}^t \right) (32) 更新后仍保持在矩阵对应位置中,可以通过设置对应掩模矩阵{{\boldsymbol{Q}}_{{\mathrm{mask}}}}进行完成:
\left.\begin{aligned} & {{\boldsymbol{Q}}_{{{\boldsymbol{Y}}_{\mathrm{C}}}}} = \left[ {\begin{array}{*{20}{c}} {{{\mathbf{0}}_L}}&{{{\mathbf{1}}_{L \times MN}}} \\ {{{\mathbf{1}}_{MN \times L}}}&{{{\mathbf{0}}_{MN}}} \end{array}} \right] \\ & {{\boldsymbol{Q}}_{\boldsymbol{Z}}}^{} = \left[ {\begin{array}{*{20}{c}} {{{\mathbf{1}}_L}}&{{{\mathbf{0}}_{L \times MN}}} \\ {{{\mathbf{0}}_{MN \times L}}}&{{{\mathbf{0}}_{MN}}} \end{array}} \right] \\ & {\boldsymbol{{Q}}_{{{\boldsymbol{U}}^*}}} = \left[ {\begin{array}{*{20}{c}} {{{\mathbf{0}}_L}}&{{{\mathbf{0}}_{L \times MN}}} \\ {{{\mathbf{0}}_{MN \times L}}}&{{{\mathbf{1}}_{MN}}} \end{array}} \right] \end{aligned}\right\} (33) 设Yimp为
{\boldsymbol{Y}}_{\mathrm{imp}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{0}}_L}}&{{{({\boldsymbol{Y}})}^{\text{H}}}} \\ {\boldsymbol{Y}}&{{{\boldsymbol{0}}_{MN}}} \end{array}} \right] (34) 则该层输出为
\begin{split} {\boldsymbol{H}}_{{\text{un}}}^{t + 1} = \,& \frac{1}{{2\rho + \lambda }}{{\boldsymbol{Q}}_{{{\boldsymbol{Y}}_{\text{C}}}}}\left(\lambda {\boldsymbol{Y}}_M + 2{{\boldsymbol{\varLambda}} ^t} + 2\rho {{\boldsymbol{T}}^t}\right) \\ \,& + {{\boldsymbol{Q}}_{\boldsymbol{Z}}}\left({\rho ^{ - 1}}{{\boldsymbol{\varLambda}} ^t} + {{\boldsymbol{T}}^t} - {\rho ^{ - 1}}{{\boldsymbol{I}}_{(NM + L)}}\right) \\ & + {{\boldsymbol{Q}}_{{{\boldsymbol{U}}^*}}}({\rho ^{ - 1}}{{\boldsymbol{\varLambda}} ^t} + {{\boldsymbol{T}}^t}) \end{split} (35) 具体表现形式如图6所示。
(3) 非线性变换层 \mathcal{S} :
该层对应ADMM迭代过程中的 {{\boldsymbol{U}}^{t + 1}} 更新步骤中的非线性变换过程,将重构层输出矩阵 {\boldsymbol{H}}_{{\text{un}}}^{t + 1} 以及可学习参数\rho 和 \lambda 作为输入,将其对应的协方差矩阵部分进行逆Toeplitz化,该过程为一个非线性变换过程,对应过程如下:
{{\boldsymbol{U}}^{t + 1}} = {S^*}({\rho ^{ - 1}}{\boldsymbol{\varLambda }}_{\mathcal{S}({\boldsymbol{U}})}^t + {\boldsymbol{T}}_{\mathcal{S}({\boldsymbol{U}})}^t) - {\rho ^{ - 1}}{{\boldsymbol{E}}_{M,N}} (36) 完成非线性变换后,将对应 {{\boldsymbol{U}}^{t + 1}} 进行块Toeplitz化处理,即式(22)的操作处理,最后输出整体整合矩阵{{\boldsymbol{H}}^{t + 1}},具体过程定义为
{{\boldsymbol{H}}^{t + 1}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{Z}}^{t + 1}}}&{{{({\boldsymbol{Y}}_{\text{C}}^{t + 1})}^{\text{H}}}} \\ {{\boldsymbol{Y}}_{\text{C}}^{t + 1}}&{\mathcal{S}({{\boldsymbol{U}}^{t + 1}})} \end{array}} \right] (37) 该层输出与数据重构层输出的不同在于 {{\boldsymbol{U}}^{t + 1}} 进行了一次非线性变换更新。具体过程如图7所示。
(4) 特征分解层\mathcal{R}:
该层结构处理过程对应式(29)中 {\boldsymbol{T}}_{{\text{un}}}^{t + 1} 求解并更新为 {{\boldsymbol{T}}^{t + 1}} 的过程,该过程主要为对非线性变换层 \mathcal{S} 的输出进行求解得到等式约束项 {\boldsymbol{T}}_{{\text{un}}}^{t + 1} 并进行特征值分解,并将其中的负特征值置零,即保证整体矩阵 {{\boldsymbol{T}}^{t + 1}} 的半正定特性,以确保后续迭代过程中的半正定约束条件成立,处理过程表示如下:
\left\{ \begin{aligned} & {\boldsymbol{T}}_{{\text{un}}}^{t + 1} = {{\boldsymbol{H}}^{t + 1}} - {\rho ^{ - 1}} {{\boldsymbol{\varLambda}} ^t} \\ & {\boldsymbol{T}}_{{\text{un}}}^{t + 1} = {\boldsymbol{G\varSigma}} {{\boldsymbol{G}}^{ - 1}} \\ & {{{\boldsymbol{T}}^{t + 1}} = {\boldsymbol{G}}{{\boldsymbol{\varSigma}} _ + }{{\boldsymbol{G}}^{ - 1}}} \end{aligned} \right. (38) 其中, {{\boldsymbol{H}}^{t + 1}} 来自同一阶段的非线性变换层的输出, {{\boldsymbol{\varLambda}} ^t} 来自前一阶段的拉格朗日乘子更新层的输出,\rho 为第t层的可学习参数。{\boldsymbol{T}}_{{\text{un}}}^{t + 1} = {\boldsymbol{G\varSigma}} {{\boldsymbol{G}}^{ - 1}}表示特征值分解,G为对应的特征向量矩阵,{\boldsymbol{\varSigma}} 是特征值对角矩阵,{{\boldsymbol{\varSigma}} _ + }为将负特征值置零后的特征值矩阵, {{\boldsymbol{T}}^{t + 1}} 为特征分解层的输出矩阵。
(5) 拉格朗日乘子更新层\mathcal{M}:
该层对应式(29)中的拉格朗日乘子矩阵 {{\boldsymbol{\varLambda}} ^{t + 1}} 的更新,具体表示为
{{\boldsymbol{\varLambda}} ^{t + 1}} = {{\boldsymbol{\varLambda}} ^t} + \rho \left( {{{\boldsymbol{T}}^{t + 1}} - {{\boldsymbol{H}}^{t + 1}}} \right) (39) 其中,该层输入包含可学习参数\rho 、前一阶段的乘子更新层的输出结果 {{\boldsymbol{\varLambda}} ^t} 、当前阶段的特征分解层输出 {{\boldsymbol{T}}^{t + 1}} 以及当前阶段的非线性变换层输出结果{{\boldsymbol{H}}^{t + 1}},最终输出该阶段更新后的拉格朗日乘子矩阵 {{\boldsymbol{\varLambda }}^{t + 1}} ,作为下一阶段各个层的输入。
(6) 最终输出 {\text{net\_end}} :
网络最终需要输出的结果是杂波协方差矩阵的估计,对应为非线性变换层输出结果中的 \mathcal{S}({{\boldsymbol{U}}^{{\text{end}}}}) 部分。因为 \mathcal{S}({{\boldsymbol{U}}^{{\text{end}}}}) 只对应杂波子空间R,最终结果还需要通过式(25)进行重构,从而输出最后训练结果{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{R}}} _{{\text{net}}}}。
3.3 数据集构建
本文所提网络采用监督学习方式,基于数据驱动进行网络参数学习,为确保网络具有优异的性能与泛化能力,需要构建合理的训练数据集。训练数据集主要需要表征杂波的空时谱特性,在双基构型下,杂波的空时频率为非线性关系,其影响因素主要有收发站位置、载机速度以及距离单元位置等,其能量分布主要受场景中杂波幅值和噪声的影响。在双基SAR成像过程中,通常具有较长的合成孔径时间,在同一距离单元内,相干脉冲数相对较大。在进行数据集构建前,需要进行相应的多普勒预滤波和距离徙动校正。在数据集构建过程中,为获取场景内足够的有效训练样本,在距离向分别选取不同距离单元的回波采样数据;针对同一距离单元内样本,进行子孔径划分,获取若干子孔径相干训练样本数据。具体处理过程如图8所示,即回波数据集中既包含不同距离单元的空时特性,也包含了同一距离单元中不同相干脉冲下对应的回波数据的空时特性。
通过上述处理获取K个空时回波数据,随机划分其中T个作为训练数据集,剩余I = K - T个作为测试样本集。为获取有效的测试标签集,通常需要对训练样本进行预求解,进而获取对应精确的空时谱。因为该数据集是对仿真场景所有点进行回波信号生成,无法获取每个点精确的空时特性和幅度特性,而采用ANM-ADMM算法需要进行人为参数设置,难以获取最为理想的超参数。并且采用该算法进行标签集生成,会造成网络训练结果也趋于人工设置结果,无法凸显网络效果。为此本文采用最为基础的ANM-STAP算法进行回波信号空时谱的求解,该算法虽然计算复杂度高,但是无需设置人工参数,可以估计得到较为精确的空时谱信息,从而作为训练和测试标签集。
3.4 网络训练
基于3.3节构建的训练数据集进行网络训练,首先进行网络初始化设置,其中网络学习参数为\rho 与\lambda 。在初始化设置前,首先进行一定理论分析确定参数的设置范围。在此基础上,采用ANM-ADMM算法进行参数验证,选取当前效果最佳的参数,从而在一定程度上加快网络收敛速度,尽可能缩短训练时间。初始化完成后,给定网络层数Q,并定义归一化均方根误差(Normalized Root Mean Square Error, NRMSE)作为损失函数,NRMSE损失可以衡量网络输出和标签矩阵之间的相似性,NRMSE损失的表达式可以表示为
{\text{los}}{{\text{s}}_{{\text{NRMSE}} }} = \frac{{{{\left\| {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{R}}} }_{{\text{net}}}} - {{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{R}}} }_{{\text{label}}}}} \right\|}_{\text{F}}}}}{{{{\left\| {{{\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{R}}} }_{{\text{label}}}}} \right\|}_{\text{F}}}}} (40) 其中, {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{R}}} _{{\text{net}}}} 表示网络的输出结果, {\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{{\boldsymbol{R}}} _{{\text{label}}}} 是标签矩阵,分子部分表示网络输出与标签之间的差值大小,分母部分为归一化操作。当NRMSE损失较小时,说明网络的输出与目标结果的相似性较高。通过降低网络损失值更新网络中对应参数,可以得到输出理想的网络结果。ANM-ADMM-Net中的超参数更新通过后向传播完成,从而使最终的估计结果逼近理想的情况。
经过训练获取最优参数后,即得到最终的最优网络结构,可以进行回波信号的空时估计求解。首先基于测试集进行网络效果测试,其中输入数据为测试集回波信号,通过网络求解输出最终杂波协方差矩阵,利用式(17)获取对应滤波器权矢量,进行滤波操作。
3.5 计算复杂度说明
本节通过对比分析ANM-ADMM-Net的计算复杂度,在分析网络计算复杂度过程中,不包括网络训练阶段中的运算过程。其中ANM-ADMM算法的计算复杂度与网络方法之间只有迭代次数的差异,单次迭代过程中,两者具有相同的计算复杂度。表1列举了ANM-CVX-STAP算法, FOCUSS-STAP算法, SBL-STAP算法以及ANM-ADMM单次迭代求解过程中对应计算复杂度。
表 1 不同算法的计算复杂度Table 1. Computational complexity of different algorithms算法 计算复杂度 ANM-CVX-STAP O({({L^2} + (2M - 1)(2N - 1) + MNL)^2}{(L + MN)^{2.5}}) FOUCSS-STAP O(NM{N_{\text{s}}}{M_{\text{d}}} + {(NM)^3} + 2{(NM)^2}{N_{\text{s}}}{M_{\text{d}}} + NM{({N_{\text{s}}}{M_{\text{d}}})^2}) SBL-STAP O(NM{N_{\text{s}}}{M_{\text{d}}} + {(NM)^3} + 3{(NM)^2}{N_{\text{s}}}{M_{\text{d}}} + 2NM{({N_{\text{s}}}{M_{\text{d}}})^2}) ANM-ADMM O({(MN + L)^3} + {(MN)^2} + 6MN + {L^2} + L) 图9展示了上述4种算法在不同阵元数和不同字典大小下的计算复杂度增长情况,其中在图9(a)中,固定字典大小为{N_{\text{s}}} = {M_{\text{d}}} = 100,训练样本数L = 8,相干脉冲数M = 10时,各类算法随着阵元数的增长情况对比可得,ANM-CVX-STAP计算复杂度最高,而ANM-ADMM的计算复杂度远小于其他3类算法。图9(b)对比了不同字典大小下4类算法的计算复杂度,其中固定阵元数与相干脉冲数N = M = 10,训练样本数L = 8,字典多普勒维度{M_{\text{d}}} = 100。相比可以看出,ANM-CVX-STAP与ANM-ADMM均不需要字典进行求解,对应计算复杂度恒定,但是ANM-CVX-STAP计算复杂度仍为最高,ANM-ADMM计算复杂度最低。在单次迭代中,ANM-ADMM算法具有最小的运算复杂度,而ANM-ADMM-Net则在此基础上,进一步降低迭代次数,因此进一步降低了运算复杂度。
为进一步验证算法的运算性能,在运行内存为32 GB,Intel(R) Core(TM) i7-13700H CPU计算平台进行100次蒙特卡罗实验仿真,对应平均运行时间如表2所示。
表 2 算法运行时间对比(s)Table 2. Run time of different algorithms (s)算法 平均运行时间 ANM-CVX-STAP 32.7150 FOUCSS-STAP 5.3151 SBL-STAP 10.9429 ANM-ADMM 1.9462 通过表2可得,ANM-ADMM算法具有最快的运行效率,而传统凸优化算法耗时最久,相比而言,ANM-ADMM-Net在ANM-ADMM算法的基础上,可以进一步减少迭代次数,更有利于实际工程应用。
4. 仿真验证
本节通过双基SAR杂波抑制仿真实验,验证ANM-ADMM-Net方法的有效性,并通过与典型SR-STAP算法以及ANM-ADMM算法进行对比分析,具体仿真参数如表3所示。
表 3 双基SAR仿真参数Table 3. Simulation parameters of BiSAR参数 数值 载频 10 GHz 信号带宽 150 MHz 脉冲重复频率 1000 Hz 天线通道数 5 相干脉冲数 10 发射机初始位置 (–5000, –3000, 4000) m 接收机初始位置 (0, –5000, 3000) m 发射机速度矢量 (0, 150, 0) m/s 接收机速度矢量 (0, 150, 0) m/s 运动目标初始位置 (0, 0, 0) m 运动目标速度矢量 (–4, 4, 0) m/s 基于表3对应的载机构型,针对某一具体SAR图像场景进行回波生成,杂噪比(Clutter signal to Noise Ratio, CNR)设置为40 dB。设置保护单元个数为8,针对某一待检测单元获取该单元左右共100个临近距离单元数,进行对应图8的空时矢量化处理,对数据进行长度为5 \times 10的滑窗取样处理,为避免SAR过长合成孔径时间的影响,针对某一滑窗获取样本,选取周围临近150个样本作为数据集,其中随机选择T = 130个作为训练数据集,其余为测试数据集。ANM-ADMM-Net网络层数设置为15层,相当于ADMM算法迭代15次,设置损失函数NRMSE小于 - 20\;{\text{dB}}时停止训练,得到最终的训练后网络。
4.1 杂波抑制效果验证
本节验证ANM-ADMM-Net空时谱估计效果,并与传统STAP,基于CVX, FOCUSS, SBL有字典的SR-STAP方法,以及无字典的ANM-ADMM进行对比验证。对比效果如图10所示。
如图10所示,当传统STAP使用充足的样本时,会导致空时谱严重展宽,从而淹没附近的微小动目标。当6种算法的估计样本数为7时,通过空时谱结果对比,传统STAP因为样本数过少,估计效果非常不理想,导致能量损失严重。基于空时字典的SR-STAP算法中,CVX与FOCUSS算法均有严重的展宽效应,SBL算法虽然展宽效应不明显,但是具有一定的精度和能量损失。相比于典型的SR-STAP算法,基于ANM-ADMM的STAP算法与ANM-ADMM-Net在空时谱估计中具有更加理想的效果。
使用SCNR损失函数进一步衡量不同方法的杂波抑制效果,结果如图11所示,通过损失曲线对比可得,ANM-ADMM-Net估计所得损失曲线凹口最窄,ANM-ADMM算法因人为设置参数的局限性,导致其凹口仍有一定展宽。在归一化频率为0时,ANM-ADMM-Net相比于传统STAP, CVX, FOCUSS, SBL以及ANM-ADMM算法其SCNR损失分别改善了15.02 dB, 11.53 dB, 10.03 dB, 6.54 dB以及4.69 dB。
为进一步验证ANM-ADMM-Net的双基SAR杂波抑制性能,下面进行了双基SAR杂波抑制效果仿真验证,仿真参数如表3所示。对回波数据进行多普勒预滤波和距离徙动校正后的回波域与成像结果如图12所示。其中图12(a)、图12(b)分别是经过距离徙动校正后的原始回波和图像域,可以看到在未进行杂波抑制处理之前,运动目标回波淹没在静止杂波之中,无法被有效检测。
基于SCNR损失函数对比结果,本节只使用传统STAP, SBL-STAP, ANM-ADMM以及ANM-ADMM-Net进行杂波抑制处理效果对比,结果如图13所示,对应成像结果如图14所示。
图13(a)使用传统STAP算法进行杂波抑制处理,大量残余杂波被保留,目标回波能量无法凸显。图13(b)使用SBL-STAP进行抑制处理后,相比于传统STAP,SBL-STAP进一步提高了杂波抑制效果,但是目标回波相对过于微弱,不利于后续动目标检测。图13(c)与图13(d)分别为ANM-ADMM方法与ANM-ADMM-Net处理结果,两种方法均过滤大量杂波能量,目标信号得到了完整的保留,但是ANM-ADMM仍残留了少量杂波,并且相比于网络,ANM-ADMM迭代次数在2000~5000次,耗时严重。图14展示了对应杂波抑制处理后的成像结果,传统STAP与SBL-STAP均留有大量杂波能量,使得后续无法实现有效的动目标检测。ANM-ADMM在保留动目标的同时,残余部分静止杂波能量,也会对后续动目标检测造成一定干扰。
为进一步比较不同方法的检测性能,图15给出了上述4种方法进行抑制后沿Y轴的图像域剖面图,对应结果如图15所示。
图15中,灰色虚线表示抑制前杂波幅度,使用传统STAP和SBL-STAP方法后目标能量并未有效保留,使用ANM-ADMM以及本文所提网络算法进行杂波抑制后,目标能量被有效保留,其中ANM-ADMM输出SCNR为21.75 dB,本文所提算法输出SCNR为22.71 dB,其他典型算法中,传统STAP未检测出运动目标,SBL-STAP输出SCNR低于7 dB,抑制效果不明显。以上结果均有效证明本文所提网络的有效性。
4.2 实测数据处理
本节中,将3.2节提出的ANM-ADMM-Net应用于实测飞行实验中。该飞行实验是在2020年10月于中国银川进行,该实验采用塞斯纳208B飞机,接收机与发射机初始坐标分别为(0, –2860, 2400) m和(1840, –1140, 2000) m,采用双基前视构型。收发载机飞行速度均为 70\;{\text{m/s}} 。在本实验中,机载双基SAR系统工作带宽为 300\;{\text{MHz}} ,脉冲重复频率为 3000\;{\text{Hz}} 。接收机采取前视模式,具有3个接收通道,通道间隔为0.8 m,发射机采取斜视模式。对应几何构型如图16所示,移动目标是一个速度为 20\;{\text{km/h}} 的皮卡车。
飞行实验得到对应回波数据,经过距离徙动校正等预处理后,其对应回波域结果如图17(a)所示。分别对实测数据采用传统STAP, SBL-STAP和所提ANM-ADMM-Net进行对比验证,其中ANM-ADMM-Net采用由仿真数据训练所得网络结构,由于地面仿真数据与实测数据具有类似的杂波分布特征,该网络结构在实测数据中也具体较好的抑制性能。结果如图17(b)、图17(c)、图17(d)所示,经传统STAP方法、SBL-STAP与本文提出的方法进行处理后,可以看到经过传统STAP处理后,仍有明显杂波保留,动目标不明显,经计算目标输出SCNR在3.28 dB左右,经过SBL-STAP处理后,目标输出SCNR为4.68 dB,并且保留了部分残余杂波。而经过ANM-ADMM-Net处理后杂波被明显抑制,而动目标仍保留在回波域中,目标输出SCNR大致在14.36 dB。对应图17(e)展示了杂波抑制前后的信号幅度,可以看出,目标SCNR被明显提升。
对应的成像结果如图18所示,图18(a)为抑制前成像结果,目标淹没在场景杂波中,传统STAP和SBL-STAP经过处理后,仍残余大量背景杂波,而本文所提方法进行处理后,有效保留了动目标信号,过滤了场景杂波。
5. 结语
本文针对在双基SAR构型下传统SR-STAP方法出现的离网问题、参数不能自适应以及运算复杂问题,提出了一种基于ANM-ADMM算法的网络化处理方法。该方法基于原子范数最小化思想,将双基SAR离散空时域求解模型转换为连续空时域下的稀疏恢复模型,有效解决了传统方法中的网格失配问题;其次,通过ADMM算法对模型进行快速迭代求解,然后针对ADMM算法中存在的参数设置困难的问题,根据迭代步骤和数据流图构建ANM-ADMM-Net网络模型,利用数据驱动的方式有效解决参数设置问题。通过仿真实验对本文所提网络方法进行了验证,结果表明ANM-ADMM-Net可以有效利用数据获取最优参数,与现有算法相比,该网络方法具有更优异的空时谱估计性能,其SCNR性能提高了3.49~10.33 dB。并在兼顾准确性的前提下,有效降低了运算复杂度,通过与现有SR-STAP定量对比,所提方法的计算复杂度降低了至少一个数量级,证明了所提方法具有更高效的运行效率。最后利用实测数据验证了本文所提方法的有效性。
-
表 1 不同算法的计算复杂度
Table 1. Computational complexity of different algorithms
算法 计算复杂度 ANM-CVX-STAP O({({L^2} + (2M - 1)(2N - 1) + MNL)^2}{(L + MN)^{2.5}}) FOUCSS-STAP O(NM{N_{\text{s}}}{M_{\text{d}}} + {(NM)^3} + 2{(NM)^2}{N_{\text{s}}}{M_{\text{d}}} + NM{({N_{\text{s}}}{M_{\text{d}}})^2}) SBL-STAP O(NM{N_{\text{s}}}{M_{\text{d}}} + {(NM)^3} + 3{(NM)^2}{N_{\text{s}}}{M_{\text{d}}} + 2NM{({N_{\text{s}}}{M_{\text{d}}})^2}) ANM-ADMM O({(MN + L)^3} + {(MN)^2} + 6MN + {L^2} + L) 表 2 算法运行时间对比(s)
Table 2. Run time of different algorithms (s)
算法 平均运行时间 ANM-CVX-STAP 32.7150 FOUCSS-STAP 5.3151 SBL-STAP 10.9429 ANM-ADMM 1.9462 表 3 双基SAR仿真参数
Table 3. Simulation parameters of BiSAR
参数 数值 载频 10 GHz 信号带宽 150 MHz 脉冲重复频率 1000 Hz 天线通道数 5 相干脉冲数 10 发射机初始位置 (–5000, –3000, 4000) m 接收机初始位置 (0, –5000, 3000) m 发射机速度矢量 (0, 150, 0) m/s 接收机速度矢量 (0, 150, 0) m/s 运动目标初始位置 (0, 0, 0) m 运动目标速度矢量 (–4, 4, 0) m/s -
[1] ROBINSON P N. Synthetic array radar[J]. IEEE Potentials, 1997, 16(1): 8–11. doi: 10.1109/45.565604. [2] 杨建宇. 双基地合成孔径雷达技术[J]. 电子科技大学学报, 2016, 45(4): 482–501. doi: 10.3969/j.issn.1001-0548.2016.04.001.YANG Jianyu. Bistatic synthetic aperture radar technology[J]. Journal of University of Electronic Science and Technology of China, 2016, 45(4): 482–501. doi: 10.3969/j.issn.1001-0548.2016.04.001. [3] WILDEN H and BRENNER A R. The SAR/GMTI airborne radar PAMIR: Technology and performance[C]. IEEE MTT-S International Microwave Symposium, Anaheim, USA, 2010: 534–537. doi: 10.1109/MWSYM.2010.5518080. [4] LI Zhongyu, WU Junjie, HUANG Yulin, et al. Ground-moving target imaging and velocity estimation based on mismatched compression for bistatic forward-looking SAR[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(6): 3277–3291. doi: 10.1109/TGRS.2016.2514494. [5] LI Junao, LI Zhongyu, YANG Qing, et al. Joint clutter suppression and moving target indication in 2-D azimuth rotated time domain for single-channel bistatic SAR[J]. IEEE Transactions on Geoscience and Remote Sensing, 2023, 61: 5202516. doi: 10.1109/TGRS.2023.3237553. [6] LI Zhongyu, WU Junjie, YANG Jianyu, et al. Bistatic SAR Clutter Suppression[M]. Singapore: Springer, 2022: 8–19. doi: 10.1007/978-981-19-0159-1. [7] 李中余. 双基地合成孔径雷达动目标检测与成像技术研究[D]. [博士论文], 电子科技大学, 2017: 11–59.LI Zhongyu. Research on bistatic SAR moving target detection and imaging technology[D]. [Ph.D. dissertation], University of Electronic Science and Technology of China, 2017: 11–59. [8] CHEN H C and MCGILLEM C D. Target motion compensation by spectrum shifting in synthetic aperture radar[J]. IEEE Transactions on Aerospace and Electronic Systems, 1992, 28(3): 895–901. doi: 10.1109/7.256313. [9] FIENUP J R. Detecting moving targets in SAR imagery by focusing[J]. IEEE Transactions on Aerospace and Electronic Systems, 2001, 37(3): 794–809. doi: 10.1109/7.953237. [10] MOREIRA J R and KEYDEL W. A new MTI-SAR approach using the reflectivity displacement method[J]. IEEE Transactions on Geoscience and Remote Sensing, 1995, 33(5): 1238–1244. doi: 10.1109/36.469488. [11] LIGHTSTONE L, FAUBERT D, and REMPEL G. Multiple phase centre DPCA for airborne radar[C]. IEEE National Radar Conference, Los Angeles, USA, 1991: 36–40. doi: 10.1109/NRC.1991.114720. [12] LI Zhongyu, LI Shanchuan, LIU Zhutian, et al. Bistatic forward-looking SAR MP-DPCA method for space–time extension clutter suppression[J]. IEEE Transactions on Geoscience and Remote Sensing, 2020, 58(9): 6565–6579. doi: 10.1109/TGRS.2020.2977982. [13] 谢文冲, 段克清, 王永良. 机载雷达空时自适应处理技术研究综述[J]. 雷达学报, 2017, 6(6): 575–586. doi: 10.12000/JR17073.XIE Wenchong, DUAN Keqing, and WANG Yongliang. Space time adaptive processing technique for airborne radar: An overview of its development and prospects[J]. Journal of Radars, 2017, 6(6): 575–586. doi: 10.12000/JR17073. [14] REED I S, MALLETT J D, and BRENNAN L E. Rapid convergence rate in adaptive arrays[J]. IEEE Transactions on Aerospace and Electronic Systems, 1974, AES-10(6): 853–863. doi: 10.1109/TAES.1974.307893. [15] KLEMM R. Comparison between monostatic and bistatic antenna configurations for STAP[J]. IEEE Transactions on Aerospace and Electronic Systems, 2000, 36(2): 596–608. doi: 10.1109/7.845248. [16] LIU Zhutian, YU Huaiqin, LI Zhongyu, et al. Non-stationary clutter suppression approach based on cascading cancellation for bistatic forward-looking SAR[C]. 2019 IEEE Radar Conference, Boston, USA, 2019: 1–5. doi: 10.1109/RADAR.2019.8835707. [17] LI Junao, LI Zhongyu, YANG Qing, et al. Efficient matrix sparse recovery STAP method based on Kronecker transform for BiSAR sea clutter suppression[J]. IEEE Transactions on Geoscience and Remote Sensing, 2024, 62: 5103218. doi: 10.1109/TGRS.2024.3362844. [18] 马泽强, 王希勤, 刘一民, 等. 基于稀疏恢复的空时二维自适应处理技术研究现状[J]. 雷达学报, 2014, 3(2): 217–228. doi: 10.3724/SP.J.1300.2014.14002.MA Zeqiang, WANG Xiqin, LIU Yimin, et al. An overview on sparse recovery-based STAP[J]. Journal of Radars, 2014, 3(2): 217–228. doi: 10.3724/SP.J.1300.2014.14002. [19] SUN Ke, ZHANG Hao, LI Gang, et al. A novel STAP algorithm using sparse recovery technique[C]. 2009 IEEE International Geoscience and Remote Sensing Symposium, Cape Town, South Africa, 2009: V-336–V-339. doi: 10.1109/IGARSS.2009.5417664. [20] TANG Gongguo, BHASKAR B N, SHAH P, et al. Compressed sensing off the grid[J]. IEEE Transactions on Information Theory, 2013, 59(11): 7465–7490. doi: 10.1109/TIT.2013.2277451. [21] YE Hongda, LI Zhongyu, LIU Zhutian, et al. Clutter-ridge matched SR-STAP technique for non-stationary clutter suppression[C]. 2020 IEEE Radar Conference, Florence, Italy, 2020: 1–4. doi: 10.1109/RadarConf2043947.2020.9266628. [22] DUAN Keqing, LIU Weijian, DUAN Guangqing, et al. Off-grid effects mitigation exploiting knowledge of the clutter ridge for sparse recovery STAP[J]. IET Radar, Sonar & Navigation, 2018, 12(5): 557–564. doi: 10.1049/iet-rsn.2017.0425. [23] LI Zhihui, ZHANG Yongshun, GE Qichao, et al. Off-grid STAP algorithm based on reduced-dimension local search orthogonal matching pursuit[C]. 2019 IEEE 4th International Conference on Signal and Image Processing, Wuxi, China, 2019: 187–191. doi: 10.1109/SIPROCESS.2019.8868509. [24] 段克清, 王泽涛, 谢文冲, 等. 一种基于联合稀疏恢复的空时自适应处理方法[J]. 雷达学报, 2014, 3(2): 229–234. doi: 10.3724/SP.J.1300.2014.13149.DUAN Keqing, WANG Zetao, XIE Wenchong, et al. A space-time adaptive processing algorithm based on joint sparse recovery[J]. Journal of Radars, 2014, 3(2): 229–234. doi: 10.3724/SP.J.1300.2014.13149. [25] HE Pengyuan, HE Shun, YANG Zhiwei, et al. An off-grid STAP algorithm based on local mesh splitting with bistatic radar system[J]. IEEE Signal Processing Letters, 2020, 27: 1355–1359. doi: 10.1109/LSP.2020.3010161. [26] LI Zhongyu, YE Hongda, LIU Zhutian, et al. Bistatic SAR clutter-ridge matched STAP method for nonstationary clutter suppression[J]. IEEE Transactions on Geoscience and Remote Sensing, 2022, 60: 5216914. doi: 10.1109/TGRS.2021.3125043. [27] FENG Weike, GUO Yiduo, ZHANG Yongshun, et al. Airborne radar space time adaptive processing based on atomic norm minimization[J]. Signal Processing, 2018, 148: 31–40. doi: 10.1016/j.sigpro.2018.02.008. [28] LI Zhongyue and WANG Tong. ADMM-based low-complexity off-grid space-time adaptive processing methods[J]. IEEE Access, 2020, 8: 206646–206658. doi: 10.1109/ACCESS.2020.3037652. [29] ZOU Bo, WANG Xin, FENG Weike, et al. DU-CG-STAP method based on sparse recovery and unsupervised learning for airborne radar clutter suppression[J]. Remote Sensing, 2022, 14(14): 3472. doi: 10.3390/rs14143472. [30] SU Hanning, BAO Qinglong, and CHEN Zengping. ADMM–net: A deep learning approach for parameter estimation of chirp signals under sub-nyquist sampling[J]. IEEE Access, 2020, 8: 75714–75727. doi: 10.1109/ACCESS.2020.2989507. [31] RICHARDS M A. The keystone transformation for correcting range migration in range-doppler processing[J]. Pulse, 2014, 1000(1). [32] LIU Zhutian, LI Zhongyu, YU Huaiqin, et al. Bistatic forward-looking SAR moving target detection method based on joint clutter cancellation in echo-image domain with three receiving channels[J]. Sensors, 2018, 18(11): 3835. doi: 10.3390/s18113835. -