二、空时自适应处理简介
空时自适应处理是一个用来描述同时处理空域和时域的自适应阵列的术语，信号的空域分量由阵列传感器收集（与所有阵列工作相同），而信号的时域分量用每个阵列传感器后等间隔延时单元产生，为此目的，一个尺寸N阵列有N个子通道(每个传感器后面对应一个）在每个子通道中，来自第j个距离单元的信号组成了M个脉冲，这M个脉冲按照雷达脉冲间隔交错。然后来自所有M个延时响应的输出进行相干累加，以产生合成的阵列响应，假设阵列输入由目标回波、杂波回波及干扰信号回波组成。
1 运算量和误差问题
自1973年STAP技术提出至今，如何在保持STAP方法杂波抑制性能的基础上降低计算复杂度、提高算法稳健性，一直是学者们所关注的重点问题。
目前，降低STAP计算复杂度的途径多是采用优化迭代方法或降维方法。考虑到共轭梯度法在矩阵求逆中的计算优势，文献[2]提出了一种结合CLEAN方法[3]和共轭梯度法的STAP算法，可获得比迭代自适应方法（iterative adaptive approach,IAA）更低的计算成本，并且改善了IAA的失配问题，一定程度上减少传统STAP方法对辅助样本数的需求。文献[4]利用了杂波协方差矩阵的斜对称结构，与局域联合处理法（joint domain localized,JDL）这一传统降维STAP方法相结合，可在降低计算复杂度的同时充分利用辅助样本。文献[5]则将多普勒域的降维变换引入互质阵雷达STAP中，有效地降低了运算量，收敛速度也有所提升。文献[6]关注于ΣΔ-3DT降维STAP方法的硬件实现，利用高速数字信号处理器（digital signal processor,DSP）的多核特性，在多核上并行处理多个子段的雷达数据；并对样本协方差矩阵的更新进行优化处理。该方法可降低计算复杂度，缩短处理时间，满足工程所需的实时性要求。
通过修正数据模型或施加额外约束，可提高STAP对误差的稳健性。针对分布式天线阵列，文献[7]证明了幂律相位噪声功率谱密度模型的适用性，同时提出一个分布式STAP系统的广义运动模型，分析了各种误差因素对分布式STAP性能的影响，这些工作有助于为提高分布式STAP的鲁棒性提供理论依据和算法改进思路。文献[8]为了提高基于参数化协方差矩阵估计（parametric covariance matrix estimation,PCE）的STAP在系统参数和杂波分布存在误差时的稳健性，首先利用稀疏恢复与Radon变换得到杂波分布，再使用归一化广义内积统计量对杂波进行修正，最后通过PCE估计出杂波协方差矩阵（clutter covariance matrix,CCM）。文献[9]将幅度约束和相位约束分别施加到空域和时域中，通过迭代优化得到相应的鲁棒权重即设计空时可分滤波器。此方法可保证STAP处理在阵列校准误差、导向矢量失配等偏差下的鲁棒性。
2 非均匀/非平稳杂波问题
STAP的实际应用场景多为非均匀、非平稳杂波环境，难以获得大量独立同分布的辅助样本，以致于对杂波协方差矩阵的估计精度下降，严重影响STAP性能。近年来，稀疏恢复（sparse recovery,SR）方法在STAP中的应用受到了广泛关注，该方法可在辅助样本数严重不足时仍然保持良好的杂波抑制性能。
三维（three-dimensional,3D)STAP是抑制非平稳杂波的重要方法，文献[10]采用基于多测量矢量的稀疏贝叶斯学习（sparse bayesian learning,SBL）算法，针对迭代过程效率低、计算量大的问题，选择了一种扩展相关向量机的快速算法来估计超参数。改进后的3D-STAP方法可在非平稳、非均匀杂波环境中较为准确地恢复出杂波三维谱，并增强了近距离杂波和俯仰角模糊杂波的可分性。大量文献已经证明，仅使用一个辅助样本时，SR-STAP算法的杂波抑制性能将会大大降低[11,12,13,14,15,16,17,18]，文献[19,20,21]致力于提高基于单样本的SR-STAP算法性能。文献[19]利用杂波回波对称性这一先验知识（图2[19]），重构一些丢失的杂波分量，得到更为准确的CCM；文献[20]没有按照传统的SR-STAP直接计算自适应权值，而是利用SR方法估计出的CCM设计降维变换矩阵，再使用基于广义旁瓣相消器（generalized sidelobe canceller,GSC）的降维STAP方法。尽管只有一个辅助样本所获得CCM的精度不足以获得合适的自适应权值向量，但它可以为选择最合适的辅助通道提供很好的支持。文献[19,20]是基于杂波空时功率谱的稀疏性，而文献[21]则是基于STAP滤波器权值的稀疏性。该方法利用角度-多普勒域权值向量的稀疏性设计降维变换矩阵，仅选择权重值较大的角度多普勒通道进行杂波抑制。与利用辅助样本进行杂波抑制的上述研究不同，文献[22]利用杂波脊斜率的先验知识，预先从接收到的数据中构造斜投影来剔除潜在的目标成分，再利用CCM的结构特点重构出CCM，实现基于SR的直接数据域（direct data domain,D3)STAP算法，不存在网格失配和系统自由度损失的问题。
除上述SR-STAP算法以外，还有其他解决非均匀、非平稳杂波问题的算法研究。文献[23]提出了一种从训练数据中剔除潜在运动目标信号和强散射体的运动目标信号抑制模型，与针对多通道雷达STAP的4种训练数据选择方法相结合并进行了比较评估。所提出的模型在保持较低的虚警概率情况下，可提高目标检测概率。针对于端射阵列机载雷达（end-fire array airborne radar,EAAR），文献[24]引入时间协方差矩阵黎曼平均距离的概念量化杂波的非平稳性，通过自适应分割将杂波距离多普勒谱分为纯噪声区、平稳杂波区和非平稳杂波区，分别对应使用脉冲多普勒（pulse Doppler,PD）、二维（two-dimensional,2D)STAP、3D-STAP进行处理。与传统3D-STAP相比，该方法降低了计算复杂度和采样要求；与传统2D-STAP相比，该方法能有效地提高非平稳杂波区域的杂波抑制性能。
三、部分源代码
clear
tic
%杂波仿真参数
N = 12; % 阵元个数
M = 10; % 相干脉冲数
CNR = 30; % 杂噪比
beta = 1; % 杂波折叠系数(beta = 2*v*T/d)
sita_a = -90:.9:90.; % 杂波单元个数
sita = sita_a*pi/180;
[NN N_bin] = size(sita);
%目标参数
sita_t = -25; % 目标DOA
omiga_t = 0.4; % 目标Doppler
SNR = 0; % 信噪比
%空间导向矢量和时间导向矢量
%空间频率和Dopple频率满足 omiga_d = beta * omiga_s
omiga_s = pi*sin(sita);
omiga_d = beta*omiga_s;
aN = zeros(N,N_bin);
bN = zeros(M,N_bin);
aN = exp(-j*[0:N-1]'*omiga_s)./sqrt(N);
bN = exp(-j*[0:M-1]'*omiga_d)./sqrt(M);
%目标空时信号
aN_t = zeros(N,1);
bN_t = zeros(M,1);
aN_t = exp(-j*pi*[0:N-1]'*sin(sita_t*pi/180))/sqrt(N);
bN_t = exp(-j*pi*[0:M-1]'*omiga_t)/sqrt(M);
S_t = zeros(M*N,1);
S_t = kron(aN_t,bN_t);
%计算杂波协方差矩阵
R = zeros(M*N,M*N);
S = zeros(M*N,N_bin);
ksai = 10^(CNR/10)*(randn(1,N_bin)+j*randn(1,N_bin))/sqrt(2); %服从正态分布的随机幅值，方差为1
for ii = 1:N_bin
S(:,ii) = kron(aN(:,ii),bN(:,ii));
R = R + ksai(ii).*(S(:,ii)*S(:,ii)');
end
%干扰协方差矩阵，杂噪比为30dB
R = R +eye(M*N); %CNR = 30dB
inv_R = inv(R); %逆矩阵
%求特征值谱
[u s v] = svd(R); %特征值分解
figure(1);
plot(10*log10(diag(s)));
title('阵元数N=12, 相干脉冲数M=10');
axis([0 120 -10 50]);
xlabel('特征值数目');
ylabel('特征值(dB)');
grid on
P_f = zeros(N_bin,N_bin);
P_min_var = zeros(N_bin,N_bin);
%求杂波谱
for ii = 1:N_bin
for jj = 1:N_bin
SS = kron(aN(:,ii),bN(:,jj));
P_f(ii,jj) = SS'*R*SS; %傅氏谱
P_min_var(ii,jj) = 1./(SS'*inv_R*SS);
end
end
%最小方差功率谱
figure(2)
mesh(sin(sita),omiga_d/pi,20*log10(abs(P_min_var)));
title('阵元数N=12, 相干脉冲数M=10');
xlabel('方位余弦');
ylabel('归一化Dopple频率');
zlabel('功率(dB)');
grid on
%空时最优权向量
w_opt = inv(R)*S_t./(S_t'*inv_R*S_t);
%w_opt = inv(RR)*a_t;
%求最优空时响应
for ii = 1:N_bin
for jj = 1:N_bin
SSS = kron(aN(:,ii),bN(:,jj));
res_opt(ii,jj) = SSS'*w_opt;
四、运行结果
五、matlab版本及参考文献
1 matlab版本
2014a
2 参考文献
[1] 沈再阳.精通MATLAB信号处理[M].清华大学出版社，2015.
[2]高宝建,彭进业,王琳,潘建寿.信号与系统——使用MATLAB分析与实现[M].清华大学出版社，2020.
[3]王文光,魏少明,任欣.信号处理与系统分析的MATLAB实现[M].电子工业出版社，2018.
