基于张量结构的快速三维稀疏贝叶斯学习STAP方法崔 宁①② 行 坤① 段克清*③ 喻忠军*①②①(中国科学院空天信息创新研究院 北京 100094)②(中国科学院大学电子电气与通信工程学院 北京 100049)③(中山大学电子与通信工程学院 广州 510006)摘 要：当机载雷达处于非正侧视工作模式时，非平稳杂波会对运动目标检测造成严重干扰。传统三维空时自适应处理(3D-STAP)方法通过构造俯仰-方位-多普勒三维自适应滤波器，可有效抑制非平稳杂波，然而巨大的系统自由度导致其在非均匀杂波环境下训练样本严重不足。虽然稀疏恢复(SR)技术可有效改善样本需求，但庞大的运算开销又使得该技术难以应用于实际。针对上述问题，该文结合机载雷达回3阶张量结构提出一种新的快速三维稀疏贝叶斯学习STAP方法，通过采用运算开销更低的张量处理将大规模矩阵求解拆分为多个小规模矩阵计算，从而大幅降低运算复杂度。详尽的数值实验验证了所提张量基SR-STAP方法可在维持SR-STAP小样本处理性能不变的基础上，将运行时间直接降低数个量级，因此是一种更适用于实际工程的SR-STAP处理方式。关键词：三维空时自适应处理；稀疏恢复；机载雷达；非平稳杂波；张量结构中图分类号：TN957.51 文献标识码：A 文章编号：2095-283X(2021)06-0919-10DOI: 10.12000/JR21140引用格式：崔宁, 行坤, 段克清, 等. 基于张量结构的快速三维稀疏贝叶斯学习STAP方法[J]. 雷达学报, 2021,10(6): 919–928. doi: 10.12000/JR21140.Reference format: CUI Ning, XING Kun, DUAN Keqing, et al. Fast tensor-based three-dimensional sparseBayesian learning space-time adaptive processing method[J]. Journal of Radars, 2021, 10(6): 919–928. doi:10.12000/JR21140.Fast Tensor-based Three-dimensional Sparse Bayesian LearningSpace-Time Adaptive Processing MethodCUI Ning①② XING Kun① DUAN Keqing*③ YU Zhongjun*①②①(Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China)②(School of Electronic, Electrical and Communication Engineering, University ofChinese Academy of Sciences, Beijing 100049, China)③(School of Electronics and Communication Engineering, Sun Yat-sen University,Guangzhou 510006, China)Abstract: When airborne radar is applied to the non-side-looking mode, moving target detection performanceconsiderably degrades because of the nonstationary clutter. Conventional three-dimensional (3D) Space-TimeAdaptive Processing (STAP) can effectively eliminate the nonstationary clutter via adaptively constructing anelevation-azimuth-Doppler 3D filter. However, large system degrees of freedom lead to a shortage of trainingsamples in a heterogeneous environment. Although introducing the Sparse Recovery (SR) technologysubstantially reduces the sample requirement, the practical application of this technology is limited bycomputational complexities. To solve the above problems, this paper proposes a fast 3D sparse Bayesian 收稿日期：2021-09-26；改回日期：2021-12-03；网络出版：2021-12-23*通信作者： 喻忠军 yuzj@ucas.ac.cn； 段克清 duankeqing@aliyun.com*Corresponding Author: YU Zhongjun, yuzj@ucas.ac.cn; DUAN Keqing, duankeqing@aliyun.com基金项目：国家自然科学基金(61871397)Foundation Item: The National Natural Science Foundation of China (61871397)责任主编：冯大政 Corresponding Editor: FENG Dazheng第10卷第6期 雷 达 学 报 Vol. 10No. 62021年12月 Journal of Radars Dec. 2021learning STAP, based on the third-order tensor structure of echo data. In the proposed method, large-scalematrix calculation is decomposed into small-scale matrix calculation using a low-complexity tensor-basedoperation, thus considerably reducing the computational load. Exhaustive numerical experiments verify that theproposed method directly reduces the computational load by several orders of magnitude compared with that ofthe existing SR-STAP algorithms, while maintaining the SR-STAP performance. Therefore, the tensor-basedmethod is a superior processing method than the vector-based method in engineering.Key words: Three-dimensional Space-Time Adaptive Processing (3D-STAP); Sparse Recovery (SR); Airborneradar; Non-stationary clutter; Tensor structure 1 引言空时自适应处理(Space-Time AdaptiveProcessing, STAP)作为一种实用的杂波抑制技术，可有效地改善机载雷达在强杂波环境下的运动目标检测能力。该技术主要利用杂波方位和多普勒之间的耦合特性，在空时二维平面内构建最优自适应滤波器，因此具有比传统多普勒滤波器更好的杂波抑制能力。然而，待检测区域实际杂波分布特性通常是未知的，需要利用待检测距离单元(CellUnder Test, CUT)相邻数据作为训练样本经似然估计获取。为分析杂波协方差矩阵(ClutterCovariance Matrix, CCM)估计误差对STAP性能影响，Reed, Mallet和Brenna共同提出RMB准则，即输出信杂比相较于最优性能损失不超过3 dB所需最小独立同分布(Independent and IdenticallyDistributed, IID)训练样本数为2倍系统自由度(De-grees Of Freedom, DOF)[1]。当机载雷达处于非正侧工作时，杂波多普勒频率与空间方位关系随距离改变而变化，即存在杂波非平稳现象。同时，由于机载预警雷达常工作于中/高脉冲重复频率(Pulse Repetition Frequency,PRF)模式，因此雷达回波往往存在严重的距离模糊。上述非平稳杂波和距离模糊交织在一起，进而导致后续传统二维STAP (Two-dimensional STAP,2D-STAP)处理存在以下问题：(1)各距离单元杂波分布特性不同，使得用于估计CCM的IID训练样本数严重不足，引起STAP性能恶化；(2)CUT中存在多个具有不同多普勒频率的方位主瓣杂波，大大增加了盲速范围。为消除非平稳杂波影响，国外学者陆续提出各类杂波谱补偿方法，如多普勒补偿法、高阶多普勒补偿法、角度-多普勒补偿法、自适应角度-多普勒补偿法等[2]，通过各类数学变换，将不同距离单元杂波谱校正一致。然而，上述补偿类方法无法同时补偿由距离模糊导致的重叠杂波，因此仅适用于低PRF工作模式[3]。三维STAP(Three-dimensional STAP, 3D-STAP)在方位-多普勒二维处理的基础上利用平面阵列增加了俯仰维度自适应处理能力，可在距离模糊情况下取得很好的杂波抑制效果。然而，更多的系统DOF令3D-STAP处理所需IID样本数远大于传统2D-STAP，为解决上述问题，降维技术(Reduced-Dimensional, RD)由二维拓展至三维，形成3D-RD-STAP方法[4]。RD处理可将IID样本需求降至2倍降维后的系统DOF，但在极端非均匀杂波环境下依然难以获得充足的训练样本。近年来，与稀疏恢复(Sparse Recovery, SR)技术相结合的STAP处理方法受到了广泛关注，大量研究表明SR-STAP只需极少训练样本甚至单样本即可获得良好的杂波抑制能力[5,6]。最近有学者将SR应用于3D-STAP并提出3D-SR-STAP方法[7]，结果表明仅需数个样本3D-SR-STAP即可实现非平稳杂波的有效抑制。虽然SR-STAP在小样本条件下相比传统方法具有显著的优势，但巨大的计算开销严重制约了它的实际应用。目前，已有学者针对SR-STAP计算开销问题提出许多快速算法，例如，结合降维方法形成局域稀疏处理[8]，通过谱先验知识辅助对空时格点预筛选[9]以及设计SR-STAP快速处理算法[10]。然而，以上方法均是对矢量化后的回波数据进行计算，所需处理的字典维度将随各维度数据增加而急剧上升，对硬件处理平台的计算负载和存储能力提出了极高要求[11]。张量作为一种处理多维数据的高效方式，已在机器学习[12]和图像处理[13]等领域取得了广泛的研究和应用。不同于传统矢量基处理，张量处理借助数据Kronecker结构特性，将大规模矩阵运算拆解为多个小尺度矩阵运算，避免因重复计算而产生的额外开销。据此，本文结合平面阵机载雷达回波三阶结构特性提出一种基于张量的快速三维稀疏贝叶斯学习STAP (3D Sparse Bayesian Learning STAP,3D-SBL-STAP)方法。首先，建立基本的三维信号模型；然后在此模型基础上，本文将2D-SR-STAP问题拓展至3D-SR-STAP问题，并给出了多观测样本矢量(Multiple Measurement Vectors, MMV)下3D-SR-STAP优化模型；接着对所提张量MSBL方法的数学原理进行描述，并从运算复杂度对所提方法的性能优势进行分析；最后，数值实验验证了所提方法。 920 雷 达 学 报 第10卷2 三维信号模型Hvp MN dKφp如图1所示，载机平台在高度处以恒定速度沿y轴正向飞行，其上装配 个理想全向阵元构成的均匀平面机载相控阵雷达，任意相邻阵元间距均为半波长。假定待探测区域位于天线远场，则对于不同阵元，其来向回波可视为具有相同入射角。在每个相干处理间隔内，雷达通过收/发天线以恒定的PRF向外发射个窄带调频相参脉冲串并接收回波。和分别表示阵列与地面杂波块之间的方位角和俯仰角。为天线阵列摆放方向与载机飞行方向之间的夹角。NMK1经T/R组件和A/D采集处理后， 基x带回波信号矢量可表示为x=xc+n=Nr ∑i=0Nc ∑j=1i;js(j;φi) +n (1)xc Nr Nci;j n2s(;φ)2CNMK1其中， 代表杂波成分； 和 分别表示距离模糊数和单个距离环内所划分的杂波块数； 为第i个模糊距离的第j个杂波块所对应的复反射系数；通常假定为零均值、方差 的复高斯白噪声；为三维空时导向矢量，可进一步展开写为s(;φ) = [se(φ)sa(;φ)]sd(;φ) (2)其中se(φ) =[1exp(j2πdsinφ)...exp(j2π(M 1)dsinφ)]T;sa(;φ) =[1exp(j2πdcosφcos)...exp(j2π(N 1)dcosφcos)]T;sd(;φ) =[1exp(j4πvpfrcosφcos(+p))...exp(j(K 1)4πvpfrcosφcos(+p))]Tse(φ) sa(;φ) sd(;φ) fr ()T其中， , 和 分别为俯仰、方位和多普勒导向矢量；代表雷达PRF；为雷达工作波长； 为转置操作；表示Kronecker积。进而CCM可表示为Rc =E{xcxHc}=Nr ∑i=0Nc ∑j=1Nr ∑i′=0Nc ∑j′=1E{i;ji′;j′}s(j;φi)sH(j′;φi′)=Nr ∑i=0Nc ∑j=1E{ji;jj2}s(j;φi)sH(j;φi)(3) Efg ()jj ()HRc^R=1L∑Ll=1xlxHlL其中， , , 和 分别为期望、共轭、取模和共轭转置操作。实际中，由于先验杂波分布未知，常使用CUT相邻数据作为训练样本进行最大似然估计(Maximum Likelihood Estimation,MLE)来替代 ，即 , 为训练样本数。最优三维空时自适应处理优化问题定义为{minwwHRws:t: wHst= 1(4)st其中，代表目标方向的三维空时导向矢量。求解上述问题得最优3D-STAP权值为w=^R 1stsHt^R 1st(5)() 1^R其中， 为矩阵求逆操作。当MLE样本不足时，无法直接求逆，此时可通过对角加载改善性能[14]。 3 三维稀疏STAP基本原理NaMeKd Na =
aN Me =
eM Kd =
dK
a
e
d S首先，对方位-俯仰-多普勒三维连续空间做网格离散化处理，将其均匀划分为 个空间格点。 其中， , 和分别代表方位、俯仰和多普勒维划分网格数； ,和 为各维度的网格划分系数。所有格点对应的三维空时导向矢量可组成过完备字典，具体可写为S= (SeSa)Sd(6) Se =[s1e;s2e;...;sMee]2CMMeSa =[s1a;s2a;...;sNaa]2CNNa其中， 为俯仰维字典， 为方位维字典， 杂波块距离环模糊距离环模糊距离环Hvpθpθφzyx 图 1 平面阵机载雷达几何照射图Fig. 1 The geometry of planar phased-array airborne radar第6期 崔 宁等：基于张量结构的快速三维稀疏贝叶斯学习STAP方法 921Sd =[s1d;s2d;...;sKdd]2CKKd为多普勒维字典。若网格划分足够精细，包含实际杂波格点，则式(1)可重新表述为x=Sξ+n (7) ξ= [12...NaMeKd]Tξ其中， 为反射系数矢量。由于在方位-俯仰-多普勒三维空间内杂波呈现稀疏分布特性，可知中大部分元素为0[5]，满足稀疏假设，因此可通过SR来估计杂波区域。3D-SR-STAP可表述为如下优化问题minξ∥ξ∥1; s:t:∥Sξ x∥22" (8)∥∥1=∑jj ξl1∥∥22=∑jj2l2 εξ^R其中， 表示描述稀疏性的范数；表示范数，用来度量恢复结果与实际数据之间的相近程度；为误差容限，与噪声强度有关。恢复变量即为三维空时网格点对应强度，可用来构建：^R=Sdiag(jξj2)SH+I (9)diag()I其中， 表示将矢量排成对角矩阵；为对角加载值；为单位对角矩阵。在获得估计CCM后，可进一步通过式(5)计算STAP权值。上述模型基于单观测矢量(Single Measurement Vector, SMV)，但实际中常采用MMV模型以进一步减小噪声影响[15]。MMV模型下的3D-SR-STAP优化问题可描述为minΞ∥Ξ∥2;1; s:t:∥SΞ X∥2F" (10)∥∥2;1 l2/l1 l1l2∥∥2F=∑∑jj2Ξ= [ξ1;ξ2;...;ξL] ^R其中， 是 混合范数，具体为对空时维度(列)使用范数，对距离维度(行)使用范数；为斐波那契范数，用来度量恢复矩阵和原始矩阵之间的相近程度；代表恢复矩阵。相应地，MMV模型下的构造式为^R=Sdiag(1LL ∑l=1jξlj2)SH+I (11) 4 张量基处理结构分析X 2CNMKNMK1 xSa Se Sd SNMKNaMeKd SXSa Se Sd常规矢量基处理先将高维回波 排为 的向量，然后再进行处理。矢量处理实现简单，但在计算过程中，单维字典 , 和通过Kronecker乘积耦合，令完备字典尺度扩充至 ，由此产生的大量冗余项给后续操作带来严重的计算负担。张量基运算直接在多维数据域处理，无须构建高维字典，从而避免由此导致的重复计算。具体地，对3阶平面阵数据而言，在进行乘法操作时，沿方位、俯仰和多普勒维分别与单维字典 , 和 进行3次小规模矩阵乘法，每次待处理矩阵的维度仅与上一次处理有关，中间无重复项生成，因此计算效率更高。下面将对张量处理过程展开详细分析。Vi1fg i11 Mi1;i2fg i1i2Fi1;i2;...;iQfg i1i2...iQ为方便后续说明，先对操作符号进行定义：代表将数据排为 的向量；表示将数据排为 的矩阵； 将数据排为 的张量。由Kronecker结构特性知乘法操作的张量等效形式为X=Z1Sa2Se3Sd(12) Z=FNa;Me;Kd;Lfξg ()n ZZ=MN;MKfZg SHaZZ′NaNMK MeMNaK KdKNaMeNaK(NM+MeM+ KdMe)其中， 为待乘张量； 为张量的n阶模乘。以方位-俯仰-多普勒模乘顺序为例，张量模乘的具体实现过程如图2所示。首先在方位模乘中，将按方位模作切片展开，即 ，然后作矩阵相乘 得，所需乘法次数为 ；依次在俯仰模乘和多普勒模乘中采取相同操作，所需乘法次数分别为 和 ；最后可得3阶张量模乘总的乘法次数为。SHξ SHaZ SHeZ′SHeZ′′NMKNaMeKdNaK(NM+MeM+KdMe) O(NMKNaMeKd) O(NaKKdMe)通过上述分析，可知张量模乘将矢量基乘法拆解为3个小尺度矩阵相乘： , 和 ，从而令所需乘法次数 中部分乘积项变为求和项 ，进而使计算复杂度由 降至 。进一步如图3所示，对比两种处理方式可发现在矢量基乘法中存在大量由Kronecker乘积产生的冗余项，而张量运算由于分维度实现，不产生额外的冗余项，也就避免无用的运算开销。 5 三维张量稀疏贝叶斯学习STAPΞ假定中元素相互独立，且每列均服从零均值复高斯先验分布p(Ξ:;l;Γ) =1πNaMeKdjΓje ΞH:;lΓ 1Ξ:;l(13)Γ= diag(γ) γ2RNaMeKd+ 2X其中， , 为稀疏超参数。同时，令噪声满足均值为0，方差为的复高斯分布假设，则 的似然函数为p(X;Ξ;2)=1(π2)NaMeKdLe 12L ∑l=1∥X:;l SΞ:;l∥22(14) p(ΞjX;Γ;2)CN(Y;Σ)结合式(13)和式(14)，由贝叶斯准则可得后验概率 也服从复高斯分布[16]，其中{Σ=Γ ΓSH(2I+SΓSH) 1SΓY=ΓSH(2I+SΓSH) 1X(15)ΞY进而可由均值估计。未知参数采用期望最922 雷 达 学 报 第10卷大化方法计算，求解过程如下：首先，基于对数最大似然构造损耗函数L(γ;2)=LlogjCj +L ∑l=1XH:;lC 1X:;l(16)C=2I+SΓSH其中， 。利用矩阵恒等变换，式(16)可改写为L(γ;2)≜L[NaMeKd ∑i=1logi+NMKlog2+log|Σ 1|]+L ∑l=1(2∥X:;l SY:;l∥22+YH:;lΓ 1Y:;l)(17) γ2然后，分别令其对和 的1阶偏导数为0，得每次迭代的最优估计为^γi =1LL ∑l=1Tr(Y:;lYH:;l)+i;i^2=∥X SY∥2FNMKL+2NMKNaMeKd ∑i=1(1 Di) (18)Di =i;i/iTr() γi i;i其中， ; 代表矩阵的迹。以上为MSBL方法，下面引入张量形式。估计过程仅使用 ，因此将其放入循环内部求解i;i =i Qi;:T:;i(19) Q=ΓSHC 1T=SΓ CYCTT 2CNMKNaMeKd其中， , 。经简化处理后，复杂度最高步骤仅余矩阵和的计算。首先，考虑的求解，将排为4阶张量 ，则T:;:;:;i =FN;M;KfT:;ig=FN;M;K{i(Sd:;locdSa:;locaSe:;loce)}=i1Sd:;locd2Sa:;loca3Se:;loce(20)locdlocaloce iSd Sa Se其中， , 和 分别代表第个三维空时导向矢量对应字典 , 和 中的列。式(20)等效矩阵实现方式如下：T:;:;:;i =FN;M;K{VNK{iSd:;locd(Sa:;loca)T}(Se:;loce)T}(21)TNMKNaMeKd TSH由于采用矩阵存储，式(21)具有更快的求解速度。此外，建议按照维度由小到大顺序计算以获得最高效率。类似的，再次整理4阶张量为，则 的张量形式为TSH=MNMK;NaMeKdfTg(SHdSHaSHe)=MNMK;NMK{T2SHd3SHa4SHe}(22) 张量模乘张量模乘的处理结构方位字典俯仰字典多普勒字典方位维模乘俯仰维模乘多普勒维模乘KdKKdMeMeMMNKNKMNaNaKKNaMeMeMKNaNaN 图 2 张量基处理流程Fig. 2 The flow of tensor-based processing 矢量基乘法张量基乘法×××===× = 图 3 两种不同结构对比Fig. 3 The comparison of two different structures第6期 崔 宁等：基于张量结构的快速三维稀疏贝叶斯学习STAP方法 923CY至此矩阵张量重构完毕，下面将阐述的求解。ΓSHSHC 1ΓC 1C 12CNMKNMKY由于 不具备Kronecker结构无法直接分解，因此调整矩阵运算顺序为先计算 ，再计算与其相乘。将 表示为 ，进而计算可改写为Y=ΓMNaMeKd;NMK{C 11SHd2SHa3SHe}X=MNaMeKd;NMK{D⊙(C 11SHd2SHa3SHe)}X(23) ⊙D 2CNaMeKdNMK NMK NaMeKdNMK DDD()i 1其中，为Hamdamard乘积；为系数张量，定义如下：将复制 次生成系数矩阵 ，紧接着把 折叠为4阶系数张量。张量MSBL (Tensor MSBL, TMS-BL)整体流程如表1所示，其中 代表前一次迭代结果。进一步归纳可知张量处理前提为待相乘的两个矩阵至少一个具备Kronecker结构，满足使用条件即可采用张量处理。 6 计算复杂度分析rs =min(M;Nse)rsa rsa Nse本节将TMSBL算法与多种SR方法的计算复杂度进行比较，其中正交匹配追踪(OrthogonalMatching Pursuit, OMP)[17]为贪婪类算法，具有极快计算效率；多迭代自适应(Multiple IterativeAdaptive Approach, MIAA)[18]和多欠定系统聚焦式求解(Multiple Focal Underdetermined SystemSolver，MFOCUSS)[19]为SR常用方法；MSBL[16]为较新方法，其不依赖参数，因此具有较高的稳健性；多样本快速收敛SBL (Multiple Fast Conver-ging SBL, MFCSBL)[10]是在MSBL基础上改进，因此具有更快求解速度。以L个训练样本求解所需复乘法次数作为计算复杂度量化统计，同时由于后续构造CCM步骤相同，这里为简化分析不计入此步操作。表2为各算法计算复杂度，其中OMP算法稀疏度以平面阵杂波秩 来近似估计[4]， 和 分别为方位线阵杂波秩和俯仰维子阵KOMP KMIAA KMFOC KMSBL KMFCSBL KTMSBL KMSBL=KTMSBL数。 , , , , 和分别代表各算法单样本求解所需迭代次数，通常MFOCUSS和MSBL迭代次数要大于MI-AA。此外，由于张量处理不影响算法固有性能，可知 。一般而言SR-STAP仅需数个样本即可有效恢复，因此L要远小于其他参数。表 1 TMSBL算法Tab. 1 TMSBL algorithmSaSe Sd X 输入：方位字典 ，俯仰字典 ，多普勒字典 ，训练样本集合 。Ξ 输出：稀疏系数S=Se Sa Sdγ0=eNaMeKdY0= 1 20 imax 初始化：过完备字典 ，稀疏控制系数 ，均值 ，噪声方差 ，最大迭代 次数 ，收敛阈值。i = 1 :imax 1：对于 执行j= 1 :NaMeKd 2: 对于 执行localocelocd 3: 计算方位索引 ，俯仰索引 和多普勒索引T:;:;:;j=FN;M;K{VNK{jSd:;locd(Sa:;loca)T}(Se:;loce)T} 4: 5: 结束循环C=MNMK;NMK{T 2SHd3SHa4SHe}+2I 6: Y=MNaMeKd;NMK{D ⊙(C 11SHd2SHa3SHe)}X 7: j= 1 :NaMeKd 8: 对于 执行j;j =j Qj;:T:;j 9: j+1=1L∥Yj;:∥22+j;j 10: Dj =j;jj 11: 2=1NMKL∥X SY∥2F+2i 1NMKNaMeKd ∑j=1(1 Dj) 12: 结束循环 13: ∥Y Yi 1∥2F/∥Y∥2F 14: 如果 跳出循环 15：结束循环表 2 计算复杂度Tab. 2 Computational complexity方法 复乘法次数 计算复杂度OMP(NMKNaMeKd+r3s+NMKr2s+ 2NMKrs)LKOMP O(NMKNaMeKdLKOMP)MIAA(2NaMeKd(NMK)2+ (NMK)3+ (L+ 1)NaMeKdNMK+LNaMeKd)KMIAA O(NaMeKd(NMK)2KMIAA)MFOCUSS(NMKNaMeKdL+ (NMK)3+ 2(NMK)2NaMeKd+NMK(NaMeKd)2)KMFOC O(NMK(NaMeKd)2KMFOC)MSBL((NaMeKd)3+ 4NMK(NaMeKd)2+(3(NMK)2+L2+ 2NMKL)NaMeKd+4(NMK)3/3 + (NMK)2+NMKL)KMSBLO((NaMeKd)3KMSBL)MFCSBL(NaMeKd(5(NMK)2+ 2LNMK+ 4NMK+ 3)+ (NMK)3+LNMK+L)KMFCSBL O(NaMeKd(NMK)2KMFCSBL)TMSBL(NaMeKd(3NMK+NM+N+L2+ 2NMKL+NMK2+NKM2) + 2(NMK)3/3+(NMK)2(1 +Kd+Me) +NMK(NaMeKN+NMNaKd+L))KTMSBL O(NaMeKdNMK2KTMSBL)924 雷 达 学 报 第10卷NaMeKdO(NaMeKdNMK2KTMSBL)O(NMKNaMeKdLKOMP)从表2的统计结果来看，OMP由于只涉及低复杂度的矩阵乘法，具有最高的计算效率；MIAA,MFOCUSS和MSBL计算复杂度中分别包含的一次项、二次项和三次项，因此运算开销呈指数级依次增加。虽然OMP求解的复杂度要远低于后几种SR算法，但局部贪婪式搜索性能在非正侧视情况下会急剧恶化，无法满足应用要求；而后几种方法尽管估计精度较高，但高度复杂的求解过程对实际硬件处理提出了很高的要求，难以实现。所提TMSBL在运算效率上复杂度量级为 ，略高于OMP的，这表明TMSBL所需求解时间接近于OMP的量级；在性能上，TMSBL通过张量结构改进MSBL方法，之前分析已表明张量处理加速核心是利用Kronecker性质来避免冗余项计算，所以TMSBL具有与MSBL相同的处理性能。综上分析可知TMSBL借助张量计算完美地解决了速度和性能之间的矛盾，因此相比其他SR算法具有更显著优势。 7 数值实验实验基于仿真回波数据来验证所提3D-TMSBL-STAP方法的性能。在仿真过程中，天线阵列工作在前视模式，具体仿真参数如表3所示。代码中涉及的张量计算由MATLAB内部Tensor工具箱实现[20]。 7.1 三维谱估计结果比较
a= 4
e= 2
d= 4 = 10 5L= 6本节检验TMSBL的杂波谱估计效果，其中三维离散网格划分选取 , 和 ，迭代终止门限 ，训练样本数 。以CCM的最小方差无失真响应(Minimum Variance Distor-tionless Response, MVDR)谱作为对比，定义如下：P(fa;fe;fd) =1|s(fa;fe;fd)HR 1s(fa;fe;fd)|(24)s(fa;fe;fd) P(fa;fe;fd)其中， 为对应频点的三维空时导向矢量，为谱强度。
a= 8
e= 4
d= 8图4为确知CCM的MVDR谱和TMSBL稀疏恢复谱。由图4(a)可知待估计杂波主要包含两部分：近程杂波和交叠的远程杂波。从图4(b)TMSBL的估计结果来看，在小样本条件下不同杂波区域均能有效恢复。当继续增加格点划分密度至 , 和 时，如图4(c)所示，由网格失配导致的杂波展宽区域明显缩小，从而估计精度得到提升。 7.2 稀疏方法杂波抑制性能对比rs下面将3D-TMSBL-STAP的杂波抑制性能分别与现有3D-OMP-STAP, 3D-MIAA-STAP, 3D-MFOCUSS-STAP, 3D-MFCSBL-STAP和3D-MS-BL-STAP方法对比，其中网格划分、迭代条件和样本选取均与之前相同。特别地，对于3D-OMP-STAP，以杂波秩作为稀疏度估计，设置其略大于。杂波抑制性能使用信杂噪比(Signal-to-Clutter-and-Noise Ratio, SCNR)损失评价，定义如下[21]：SCNR Loss=2nNMK|wHst|2(wHRw)(25)图5为信杂噪比损失结果。在非平稳杂波条件下，3D-OMP-STAP处理性能最差，这是由于较大 (a) MVDR (c) TMSBL: 8×4×8 (b) TMSBL: 4×2×4归一化方位频率归一化多普勒频率归一化俯仰频率01 0 -1-0.5-1.0-1 0 10-10-20-30-40-50-60归一化方位频率归一化多普勒频率归一化俯仰频率01 0 -1-0.5-1.0-1 0 10-10-20-30-40-50-60归一化方位频率归一化多普勒频率归一化俯仰频率01 0 -1-0.5-1.0-1 0 10-10-20-30-40-50-60 图 4 三维杂波谱估计结果Fig. 4 The estimation of three-dimensional clutter spectrum表 3 雷达系统参数Tab. 3 Radar system parameters参数 符号 数值载机速度vp150 m/s载机高度H8000 m行/列阵元数M/N6/8脉冲数K8雷达波长0.1 m脉冲重复频率fr8100 Hz阵元间距d0.05 m主波束指向/φ–90°/0°杂噪比 – 60 dB第6期 崔 宁等：基于张量结构的快速三维稀疏贝叶斯学习STAP方法 925的网格失配导致贪婪类算法在搜索过程中受到严重干扰，从而性能损失严重。与之对比，3D-SR-STAP方法即使在该条件下依然能取得较好的抑制效果，其中SBL类方法性能最优，信杂噪比损失结果要好于其他SR方法。需要注意的是，所提3D-TMSBL-STAP方法与3D-MSBL-STAP方法性能一致，说明张量处理对算法固有杂波抑制性能无损失。 7.3 传统方法杂波抑制性能对比L= 6 rs MN本节对比传统3D-STAP方法与3D-TMSBL-STAP方法在小样本 条件下的杂波抑制性能，所选取的3D-STAP方法包括：样本矩阵求逆(Sample Matrix Inversion, SMI)[1]、局域联合处理(Joint Domain Localized, JDL)[22]、多级维纳滤波(Multistage Wiener Filter, MWF)[23]、直接数据域(Direct Data Domain, DDD)[24]和子阵合成(Sub-array Synthesis, SS)[4]。具体地，JDL在方位-俯仰-多普勒向选取5×6×5波束；MWF的分解次数设为；DDD的滑动窗口选取各维度的一半；SS合成行列阵元数为 和的十字子阵。从图6可以看到，传统3D-STAP方法在样本短缺情况下形成凹口远宽于实际主杂波区域，此时会对主杂波区域外目标强度造成严重损失。虽然3D-DDD-STAP方法仅靠单样本即可构建杂波滤波器，但由于存在孔径损失，可以看到非主杂波区域性能下降严重。而3D-TMSBL-STAP方法基于SR理论，对样本数依赖并不严重，在小样本条件下具有远好于其他传统方法的杂波抑制性能。如图7所示，以信杂噪比损失的平均值(平均SCNR损失)作为评价[4]，当可用训练样本数增加时，传统方法性能逐渐改善。JDL和MWF两种方法由于降维损失在L大于50后开始趋于稳定，此时平均信杂噪比损失相比于3D-TMSBL-STAP方法依然相差4～5 dB。而SS和SMI方法尽管最终性能要好于3D-TMSBL-STAP方法，但需要充足样本才能实现，不适应于严重非均匀杂波环境。以上结果进一步说明3D-TMSBL-STAP方法相比于传统3D-STAP方法更适用于非均匀杂波环境。 7.4 运行时间对比最后将对比不同N, M 和K大小下3D-TMSBL-STAP方法与现有3D-SR-STAP算法的实际运行时间。参数设置与之前相同，不再赘述。共选取7组数据，对应维度为俯仰×方位×多普勒，相邻数据仅在单维尺度上变化。各算法均执行于MATLAB R2019a，硬件配置为Intel Core i9-9900KF 3.6 GHz加64 GB内存，采用MATLAB内部计时器统计运行时间。NaMeKdNaMeKd如图8所示，从整体运行时间上看，MSBL要远高于其他方法1～3个数量级，原因是计算过程中涉及两个 的大规模矩阵乘法计算。紧接着是MIAA, MFCSBL和MFOCUSS，三 归一化多普勒频率-1.0 -0.5 0 0.5 1.0-60-50-40-30-20-100SCNR损失 (dB)3D-TMSBL-STAP3D-OMP-STAP3D-MFOCUSS-STAP3D-MIAA-STAP3D-MSBL-STAP3D-MFCSBL-STAPOPT 图 5 不同SR-STAP方法的SCNR损失结果Fig. 5 The SCNR loss results of different SR-STAP methodss 归一化多普勒频率-1.0 -0.5 0 0.5 1.0-60-50-40-30-20-100SCNR损失 (dB)3D-TMSBL-STAP3D-SMI-STAP OPT3D-SS-STAP 3D-JDL-STAP3D-MWF-STAP 3D-DDD-STAP 图 6 传统STAP方法的SCNR损失结果Fig. 6 The SCNR loss results of conventional STAP methods 训练样本数10 20 30 40 50 60-50-40-30-20-100平均SCNR损失 (dB)3D-TMSBL-STAP3D-SMI-STAP OPT3D-SS-STAP 3D-JDL-STAP3D-MWF-STAP 3D-DDD-STAP 图 7 平均SCNR损失结果Fig. 7 The result of average SCNR loss926 雷 达 学 报 第10卷rs者的运行时间处于中等水平。运算开销最小的是OMP和TMSBL，其中TMSBL运算时间要小于OMP，这是因为当工作于大尺度数据时，OMP受约束迭代次数较多。综上所述，可知TMSBL借助张量处理取得了显著的速度提升。 8 结论3D-SR-STAP仅需数个样本即可实现有效的非平稳杂波抑制，但巨大的计算开销限制了它的应用。本文从空时字典的Kronecker结构特性出发，提出基于张量处理的3D-MSBL-STAP方法，通过低复杂度张量操作加速求解过程，从而获得远快于矢量基方法的运算速度。仿真结果表明在本文条件下，所提方法可在维持原算法性能的同时，将计算速度提升2～3个数量级，因此更具实际工程价值。此外，还可将所提张量处理方法拓展应用于其他矢量基求解算法，以提升其运算速度。未来，本团队将改进张量的读写方式，以提升其处理速度，并尝试融入并行处理，以期获得进一步加速。参 考 文 献REED I S, MALLETT J D, and BRENNAN L E. Rapidconvergence rate in adaptive arrays[J]. IEEE Transactionson Aerospace and Electronic Systems, 1974, AES-10(6):853–863. doi: 10.1109/TAES.1974.307893.[1]谢文冲, 段克清, 王永良. 机载雷达空时自适应处理技术研究综述[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 ofRadars, 2017, 6(6): 575–586. doi: 10.12000/JR17073.[2]MENG Xiangdong, WANG Tong, WU Jianxin, et al. Short-range clutter suppression for airborne radar by utilizingprefiltering in elevation[J]. IEEE Geoscience and RemoteSensing Letters, 2009, 6(2): 268–272. doi: 10.1109/LGRS.2008.2012126.[3]DUAN Keqing, XU Hong, YUAN Huadong, et al. Reduced-DOF three-dimensional STAP via subarray synthesis fornonsidelooking planar array airborne radar[J]. IEEETransactions on Aerospace and Electronic Systems, 2020,56(4): 3311–3325. doi: 10.1109/TAES.2019.2958174.[4]SUN Ke, MENG Huadong, WANG Yongliang, et al. Directdata domain STAP using sparse representation of clutterspectrum[J]. Signal Processing, 2011, 91(9): 2222–2236. doi:10.1016/j.sigpro.2011.04.006.[5]YANG Zhaocheng, LI Xiang, WANG Hongqiang, et al. Onclutter sparsity analysis in space-time adaptive processingairborne radar[J]. IEEE Geoscience and Remote SensingLetters, 2013, 10(5): 1214–1218. doi: 10.1109/LGRS.2012.2236639.[6]DUAN Keqing, XU Hong, YUAN Huadong, et al. Three-dimensional sparse recovery space-time adaptive processingfor airborne radar[J]. The Journal of Engineering, 2019,2019(19): 5478–5482. doi: 10.1049/joe.2019.0343.[7]GUO Yiduo, LIAO Guisheng, and FENG Weike. Sparserepresentation based algorithm for airborne radar in beam-space post-Doppler reduced-dimension space-time adaptiveprocessing[J]. IEEE Access, 2017, 5: 5896–5903. doi:10.1109/ACCESS.2017.2689325.[8]HAN Sudan, FAN Chongyi, and HUANG Xiaotao. A novelSTAP based on spectrum-aided reduced-dimension cluttersparse recovery[J]. IEEE Geoscience and Remote SensingLetters, 2017, 14(2): 213–217. doi: 10.1109/LGRS.2016.2635104.[9]WANG Zetao, XIE Wenchong, DUAN Keqing, et al.Clutter suppression algorithm based on fast convergingsparse Bayesian learning for airborne radar[J]. SignalProcessing, 2017, 130: 159–168. doi: 10.1016/j.sigpro.2016.06.023.[10]QIU Wei, ZHOU Jianxiong, ZHAO Hongzhong, et al.Three-dimensional sparse turntable microwave imagingbased on compressive sensing[J]. IEEE Geoscience andRemote Sensing Letters, 2015, 12(4): 826–830. doi: 10.1109/LGRS.2014.2363238.[11]SIDIROPOULOS N D, DE LATHAUWER L, FU Xiao,et al. Tensor decomposition for signal processing and machinelearning[J]. IEEE Transactions on Signal Processing, 2017,65(13): 3551–3582. doi: 10.1109/TSP.2017.2690524.[12]ZHAO Rongqiang, WANG Qiang, FU Jun, et al. Exploitingblock-sparsity for hyperspectral Kronecker compressivesensing: A tensor-based Bayesian method[J]. IEEE[13] 数据大小4×8×85×8×86×8×86×10×86×12×86×12×106×12×12104102运行时间 (s)3D-TMSBL-STAP3D-OMP-STAP3D-MFOCUSS-STAP3D-MIAA-STAP3D-MSBL-STAP3D-MFCSBL-STAP 图 8 不同数据大小下实际运行时间Fig. 8 The running time of different input data size第6期 崔 宁等：基于张量结构的快速三维稀疏贝叶斯学习STAP方法 927Transactions on Image Processing, 2020, 29: 1654–1668.doi: 10.1109/TIP.2019.2944722.姜磊, 王彤. 机载雷达自适应对角加载参数估计方法[J]. 电子与信息学报, 2016, 38(7): 1752–1757. doi: 10.11999/JEIT151003.JIANG Lei and WANG Tong. An adaptive estimationmethod of diagonal loading parameter for airborne radar[J].Journal of Electronics & Information Technology, 2016,38(7): 1752–1757. doi: 10.11999/JEIT151003.[14]MA Zeqiang, LIU Yimin, MENG Huadong, et al. Jointlysparse recovery of multiple snapshots in STAP[C]. 2013IEEE Radar Conference, Ottawa, Canada, 2013: 1–4. doi:10.1109/RADAR.2013.6586083[15]DUAN Keqing, WANG Zetao, XIE Wenchong, et al.Sparsity-based STAP algorithm with multiple measurementvectors via sparse Bayesian learning strategy for airborneradar[J]. IET Signal Processing, 2017, 11(5): 544–553. doi:10.1049/iet-spr.2016.0183.[16]TROPP J A and GILBERT A C. Signal recovery fromrandom measurements via orthogonal matching pursuit[J].IEEE Transactions on Information Theory, 2007, 53(12):4655–4666. doi: 10.1109/TIT.2007.909108.[17]YANG Zhaocheng, LI Xiang, WANG Hongqiang, et al.Adaptive clutter suppression based on iterative adaptiveapproach for airborne radar[J]. Signal Processing, 2013,93(12): 3567–3577. doi: 10.1016/j.sigpro.2013.03.033.[18]YANG Shuyuan, LI Bin, WANG Min, et al. Compressivedirection-of-arrival estimation via regularized multiple[19]measurement FOCUSS algorithm[C]. 2014 InternationalJoint Conference on Neural Networks, New York, USA,2014: 2800–2803. doi: 10.1109/IJCNN.2014.6889967.doi:10.1109/IJCNN.2014.6889967.BADER B W and KOLDA T G. Matlab tensor toolboxversion 2.6[EB/OL]. http://www.sandia.gov/tgkolda/TensorToolbox. 2015.[20]GUERCI J R. Space-Time Adaptive Processing forRadar[M]. Boston: Artech House, 2014: 52–73.[21]HALE T B, TEMPLE M A, RAQUET J F, et al. Localizedthree-dimensional adaptive spatial-temporal processing forairborne radar[C]. 2002 International Radar Conference,Edinburgh, UK, 2002: 191–195. doi: 10.1049/cp:20020275.[22]洪玺, 王文杰, 殷勤业. 基于多级维纳滤波器的空时自适应信号处理及其在无线通信系统中的应用[J]. 信号处理, 2017,33(3): 430–436. doi: 10.16798/j.issn.1003-0530.2017.03.026.HONG Xi, WANG Wenjie, and YIN Qinye. Multistagewiener filter based space and time adaptive signalprocessing and its application in wireless communicationsystem[J]. Journal of Signal Processing, 2017, 33(3):430–436. doi: 10.16798/j.issn.1003-0530.2017.03.026.[23]王璐, 吴仁彪. 直接数据域空时自适应单脉冲方法[J]. 系统工程与电子技术, 2016, 38(12): 2738–2744. doi: 10.3969/j.issn.1001-506X.2016.12.09.WANG Lu and WU Renbiao. Direct data domain space-time adaptive monopulse method[J]. Systems Engineeringand Electronics, 2016, 38(12): 2738–2744. doi: 10.3969/j.issn.1001-506X.2016.12.09.[24]作 者 简 介 崔 宁(1995–)，男，辽宁抚顺人，中国科学院大学在读博士研究生，主要研究方向为空时自适应处理和雷达动目标检测等。行 坤(1980–)，男，陕西西安人，中国科学院空天信息创新研究院副研究员，中国科学院大学硕士生导师，主要研究方向为机载多功能雷达、地基监视雷达系统总体设计与性能仿真，地面、海面、空中动目标检测与跟踪技术等。段克清(1981–)，男，河北石家庄人，中山大学副教授，硕士生导师，主要研究方向为阵列信号处理、空时自适应处理、压缩感知、深度学习及其在雷达系统中的应用。喻忠军(1980–)，男，四川成都人，中国科学院空天信息创新研究院研究员，中国科学院大学博士生导师，主要研究方向为微系统微集成、微波毫米波电路模块、先进相控阵天馈技术等。928 雷 达 学 报 第10卷
