简介：面向本科、硕士阶段雷达信号处理教研学习的一套MATLAB实现，聚焦机载雷达时空自适应处理（STAP）核心算法，兼顾理论演示、仿真复现与工程实践。资源包含完整运行结果与可视化脚本，可直接在MATLAB 2019a中运行复现，也可作为基础教程配合教材，完成对照推导、仿真验证与结果分析。压缩包共204个文件，大小约4.25MB，其中54个m脚本为源码主体，106个png图片展示各环节输出结果，44个html文件提供图形化中间过程预览，便于快速定位关键步骤、检查运行是否正常；按fig编号组织多组实验，脚本注释清晰，可修改参数重新运行。已有167人浏览学习，适合雷达课程设计、毕业设计、教研备课及自学入门。针对运行中遇到的问题，作者支持私信沟通，能有效降低上手门槛，也可支撑课堂演示、结课报告与论文复现。
1. 为什么机载雷达 STAP 值得用一套 MATLAB 代码反复拆
拿到这份机载雷达时空自适应处理附 matlab 代码时，我建议先别急着打开 fig65.html。STAP 的入门难点其实不在公式推导，而在于回波数据怎么组织、训练样本怎么取、权矢量算出来之后怎么验证。这套基于 MATLAB 2019a 的代码包把阵元维、脉冲维和距离门维都做成了可运行的脚本，运行后直接输出 fig48、fig51、fig65 到 fig69 等 HTML 结果图。本科、硕士阶段做教研项目时，最有效的用法不是一键跑完，而是把 N、M、β 这几个关键量逐个改一遍，观察特征谱和多普勒凹口的变化。它会让你直观理解：空时自适应处理之所以比纯空域滤波强，是因为它同时利用空间、时间两个自由度去拟合杂波。
2. 空时自适应处理的数学模型：先把手头数据写成可计算的张量
机载雷达在正侧视工作时，地面静止杂波会同时出现在空间频率和多普勒频率上，二者以直线形式耦合在一起；运动目标的多普勒则偏离这根杂波脊。这是 STAP 能工作的物理前提，也是仿真中检验数据组织是否正确的最直观标尺。因此，动手写代码之前，必须先把“空时快拍”这个结构定义清楚，否则后面算协方差矩阵、求权矢量都会对不上号。
2.1 空时快拍：把回波组织成 N×M 矩阵再按列展开
一个相干处理间隔（CPI）内，阵列有 N 个阵元，每个阵元在一个 CPI 内采样 M 个脉冲。对某一个距离门，收到的数据可以排成一个 N×M 矩阵，行是阵元维，列是脉冲维。MATLAB 对这种矩阵做x(:)展开时，按列优先读取，也就是第一列前 N 个元素对应第 1 个脉冲在所有阵元上的采样，接着是第 2 个脉冲的信号。这个顺序直接影响导向矢量的构造，也是最容易写反的地方。
% 假设 sensor_data 是 N×M×L 的三维回波数据，L 为距离门数量 N = 8; % 阵元数 M = 16; % 脉冲数 L = 128; % 距离门数 l_gate = 64; % 取中间距离门观察 x_nm = sensor_data(:, :, l_gate); % 单个距离门的 N×M 回波 x = x_nm(:); % N*M×1 空时快拍向量这里的关键是x(:)后，索引前 N 个数据是第 1 个脉冲，第 N+1 到 2N 个数据是第 2 个脉冲。如果后面用kron生成导向矢量，就必须保持同样的顺序。常见做法是定义目标空间导向矢量s_s、时间导向矢量s_t，再构造kron(s_t, s_s)，因为s_t的每个元素对应一个脉冲块，块内顺序由s_s决定。
2.2 最优 STAP 权矢量：从最小方差响应说起
STAP 的优化目标可以写成：在保证目标方向增益不变的约束下，最小化输出功率。用公式表达就是最小化w^H R w，同时满足w^H s = 1。其中R是杂波加噪声协方差矩阵，s是目标空时导向矢量。解出来的最优权是：
w = α R^(-1) s，其中α = 1 / (s^H R^(-1) s)。
在 MATLAB 中，s^H对应s'，不是s.'。很多初学者在这里把转置和共轭转置混用，导致算出的权矢量相位错误，最终看到的输出根本不是零陷而是随机起伏。理论中的R是统计期望，实际系统只能拿到有限样本估计出的R_hat，所以后面所有工程技巧都围绕着如何让R_hat更接近真实干扰环境。
2.3 杂波秩：为什么不是 N×M 而是远小于 N×M
如果只看自由度，STAP 似乎需要在一个 N×M 维空间里求逆。但机载正侧视雷达的地杂波是有结构的，杂波协方差矩阵的秩大约只有N + (M-1) * β，其中β = 2v / (d * PRF)是杂波脊的斜率系数。这个结论非常重要，它告诉我们全维 STAP 的理论自由度并没有被全部用完，也为后面降秩处理提供了依据。
| 参数 | 符号 | 典型设置 | 对 STAP 的影响 |
|---|---|---|---|
| 阵元数 | N | 8 | 决定空间自由度和波束宽度 |
| 脉冲数 | M | 16 | 决定多普勒分辨率和时间自由度 |
| 载机速度 | v | 100 m/s | 改变杂波多普勒展宽程度 |
| 阵元间距 | d | λ/2 | 影响空间采样和无模糊范围 |
| PRF | fr | 2000 Hz | 决定多普勒无模糊区间 |
| 杂波秩系数 | β | 约 0.67 | 决定可检测的自由度数量 |
按照这份参数，β约等于 0.67，杂波秩大约为8 + 15 * 0.67 ≈ 18。但全维 STAP 的自由度是8 * 16 = 128，如果直接估计 128 维协方差矩阵，训练样本需求会非常夸张。这就是为什么仿真代码里通常会先输出特征谱图，让你直观看到只有少数大特征值占主导。
3. MATLAB 仿真实现：从回波生成到最优 STAP 权矢量
仿真代码的核心不是把公式敲进 MATLAB，而是先把几何模型和坐标系统一。很多人在这一步栽跟头：空间频率、多普勒频率、波达方向三者的定义互相矛盾，画出的图永远不对。下面这段实现按照正侧视模型展开，杂波脊会穿过原点，方便后续验证。
3.1 从正侧视几何确定仿真参数
正侧视意味着载机飞行方向与天线法向垂直，地面散射点相对阵列的运动速度是v * cos(θ)，这里θ是来波方向与载机速度方向的夹角，正侧视时θ = 90°。参数设置直接影响杂波秩和凹口位置，所以把它们集中放在脚本开头，方便反复修改。
N = 8; % 阵元数，空间自由度 M = 16; % 一个 CPI 内脉冲数，时间自由度 lambda = 0.3; % 载波波长，单位：米 d = lambda / 2; % 阵元间距，半波长设置 v = 100; % 载机速度，单位：m/s fr = 2000; % 脉冲重复频率 PRF，单位：Hz beta = 2 * v / (d * fr); % 杂波秩斜率系数 beta fs_dop = 2 * v / lambda; % 最大多普勒频率，用于后续坐标这段代码里最值得注意的是beta。它不是一个可以随便给的数，而是由速度、阵元间距、PRF 共同决定的。改变速度或 PRF 后，杂波凹口的斜率会跟着变，这也直接决定局域 STAP 需要选取多少个邻域通道。lambda取 0.3 米对应对应 1 GHz 左右频段，教学用足够了。
3.2 生成杂波、目标与噪声，构造训练样本集
生成杂波时，常用做法是把连续地物离散成大量独立散射点，每个散射点具有随机的幅度和位置，位于同一个距离门内的散射点相干叠加。这种建模方式能保留杂波的统计特性，又能用 MATLAB 向量化快速计算。
num_scat = 500; % 散射点数量 theta_scat = linspace(1, 179, num_scat); % 从机头到机尾的方向角 clutter = zeros(N, M); for k = 1:num_scat f_s = d / lambda * cosd(theta_scat(k)); % 空间频率 f_d = 2 * v / lambda * cosd(theta_scat(k)); % 多普勒频率 s_s = exp(1j * 2 * pi * f_s * (0:N-1)') / sqrt(N); s_t = exp(1j * 2 * pi * f_d * (0:M-1)') / sqrt(M); clutter = clutter + s_s * s_t.'; % 外积构成 N×M 回波 end x_c = clutter(:); % 杂波快拍 noise_power = 1e-6; x_n = sqrt(noise_power/2) * (randn(N*M,1) + 1j*randn(N*M,1)); x = x_c + x_n; % 当前距离门的观测快拍这段代码用s_s * s_t.'生成单个散射点的空时矩阵，再把所有散射点累加。注意这里s_t是列向量，转置后变成行向量，乘出来的矩阵行对应阵元、列对应脉冲。最后用x(:)展开时，向量中的排列顺序和前面讲的kron(s_t, s_s)一致。如果后面构建目标导向矢量，也需要保持同样顺序。
如果要加入运动目标，需要单独构造目标导向矢量，目标速度不只是沿视线方向，所以多普勒频率和空间频率没有固定关系。
t_theta = 60; % 目标方向角 t_v = 30; % 目标径向速度，单位：m/s f_st = d / lambda * cosd(t_theta); f_dt = 2 * t_v / lambda; s_st = exp(1j * 2 * pi * f_st * (0:N-1)') / sqrt(N); s_tt = exp(1j * 2 * pi * f_dt * (0:M-1)') / sqrt(M); s_target = kron(s_tt, s_st); % 与 x(:) 的排列一致 amp_t = 0.01; % 目标相对幅度 x = x + amp_t * s_target; % 将目标注入观测快拍训练样本通常从相邻距离门获取，但必须避开目标所在距离门及其附近 2 到 3 个保护距离门，否则目标信号会进入协方差矩阵估计，导致自适应权把目标当作干扰抑制掉。
3.3 采样协方差矩阵与最优权：不要直接用 inv 求逆
得到训练样本后，协方差矩阵用R_hat = X * X' / K估计。X 是 N*M 行、K 列的矩阵，K 是训练样本数。理论上样本数要至少大于自由度，实际上最好取自由度的 2 到 5 倍。对角加载是为了保证矩阵可逆，同时牺牲极小一部分自适应性能，换取数值稳定。
K = 4 * N * M; % 训练样本数，手动设定 X_train = zeros(N*M, K); for k = 1:K X_train(:, k) = generate_snapshot_no_target(k); end R_hat = (X_train * X_train') / K; diag_ratio = 1e-4; % 对角加载比例 R_ld = R_hat + diag_ratio * trace(R_hat) / (N*M) * eye(N*M); w = R_ld \ s_target; w = w / (s_target' * w); % 归一化，满足约束这段代码里我用\而不是inv。inv会显式计算逆矩阵，数值稳定性更差；反除法的求解速度更快，且对接近奇异的情况更鲁棒。trace(R_hat) / (N*M)相当于每个维度的平均功率，用它作为对角加载的基准，不会因为目标信号功率过大而把加载量算偏。通常diag_ratio从 1e-6 开始尝试，看到输出曲线失稳再逐步增大。
3.4 从 fig48 到 fig69：一套结果图按什么顺序看
打包里的 fig48.html、fig51.html、fig52a/52b、fig64 到 fig69.html 都是 MATLAB publish 导出的 HTML 图形，浏览器双击即可打开。按我的习惯，第一看 fig65 或 fig66，检查杂波凹口是否出现在预期位置；第二看 fig68 和 fig69，确认目标多普勒附近的输出是否保留；最后回看 fig48 和 fig51，判断数据维度和幅度是否正常。如果自己修改参数后，凹口消失或者出现多个伪零陷，几乎可以锁定是角度定义、kron顺序或训练样本混入了目标信号。
4. 降秩处理与参数调优：让权矢量在杂波秩面前真正可用
全维 STAP 在理论上是干净的，但在实际参数下很难直接用。自由度越大，需要的独立同分布训练样本就越多；机载雷达飞过一个均匀杂波区的时间窗口有限，可用距离门不多，因此必须把自由度压下来，让权矢量在有限样本下仍然稳定。
4.1 全维 STAP 的样本需求在机载环境中不现实
以 N=8、M=16 为例，全维 STAP 要做 128 维协方差矩阵求逆。杂波秩大约 18，但矩阵求逆仍然希望样本数达到 256 个以上。地面反射特性跨越几十公里后变化很大，事实上很难找到足够多“同分布”的距离门做样本。降秩 STAP 的思路就是只保留与杂波强相关的局部空时区域，比如目标多普勒附近的三五个多普勒通道，以及目标方向附近的三五个空间波束，构成小型化矩阵，再在这个小矩阵里做自适应滤波。
4.2 局域 STAP 实现：用 3×3 空时模块替换全维矩阵
局域 STAP 是工程上最常见的一种降维方案。下面代码把 128 维降到 9 维，生成一个选取矩阵 T，用 T 对快拍、导向矢量和协方差矩阵同时做投影。
doppler_ref = 8; % 目标所在多普勒通道索引 b_dopp = doppler_ref + (-1:1); % 多普勒邻域，共3个通道 b_elem = round(N/2) + (-1:1); % 空间波束邻域，共3个阵元 T = zeros(N*M, 9); idx = 0; for dd = 1:3 for ee = 1:3 idx = idx + 1; % x(:) 中阵元维变化最快，所以线性索引如下 lin_idx = (b_dopp(dd) - 1) * N + b_elem(ee); T(lin_idx, idx) = 1; end end s_red = T' * s_target; % 降维后的导向矢量 R_red = T' * R_hat * T; % 降维后的协方差矩阵 w_red = R_red \ s_red; w_red = w_red / (s_red' * w_red); % 归一化权这段代码里lin_idx的构造必须和x(:)的排列一致，也就是先变阵元索引，再变多普勒通道索引。如果写成(b_elem-1)*N+b_dopp，选出来的是另一个子空间，结果表现会完全失真。这也是很多教材代码跑不出理想凹口的原因之一。降维后只有 9 个自由度，训练样本需求从 256 降到 18 左右，在单帧数据里基本可以满足。
全维和局域 STAP 的差异对比如下：
| 指标 | 全维 STAP | 局域 3×3 STAP |
|---|---|---|
| 自由度 | N*M=128 | 9 |
| 训练样本需求 | ≥256 | ≥18 |
| 计算量 | 高，适合离线分析 | 低，适合实时处理 |
| 杂波抑制能力 | 理论最优 | 局部最优，凹口略宽 |
4.3 参数调优顺序：从 rcond 到对角加载，再到多普勒凹口
拿到代码后不要急着跑完整流程，先按下面顺序调参数。第一步检查rcond(R_hat)，如果小于 1e-12，说明训练样本折叠或者矩阵病态，必须增加对角加载。第二步画 Capon 谱，看看杂波脊是否和beta计算结果吻合。
capon = 1 ./ abs(s_target' * (R_ld \ s_target)); % Capon 谱输出Capon 谱是功率的无偏估计，杂波所在位置会出现明显峰值。若峰值位置和理论杂波脊不一致，先检查坐标定义，而不是急着调权矢量。第三步才是调整对角加载系数。对角加载太大，凹口变宽，目标也会被压掉；太小，矩阵奇异，输出噪声反而增大。我一般以 1e-6 为起点，每次乘 10，观察目标输出功率达到稳定就不再增加。
5. 用改善因子曲线和 HTML 结果图做一次可复现的验收
只看权矢量系数无法判断 STAP 算法是否有效，更可靠的方式是计算改善因子曲线。改善因子是输出信干噪比和输入信干噪比的比值，反映了自适应处理对不同多普勒频率的保留和抑制程度。对机载雷达教研项目来说，这条曲线比某个点的方向图更有解释力。
5.1 改善因子曲线：比驻波图更接近 STAP 的本质
改善因子曲线的横轴是多普勒频率，纵轴是分贝形式的改善能力。杂波所在频率处会出现深凹口，目标频率处保持较高增益。计算时固定权矢量不变，扫描不同多普勒通道的导向矢量，观察输出响应。
fd_axis = linspace(-fs_dop, fs_dop, 256); IF_curve = zeros(size(fd_axis)); for k = 1:numel(fd_axis) s_dop = exp(1j * 2 * pi * fd_axis(k) * (0:M-1)' / fs_dop); s_cell = kron(s_dop, s_st); % 保持与 x(:) 的结构一致 IF_curve(k) = abs(w_red' * s_cell)^2 / (w_red' * R_red * w_red); end plot(fd_axis, 10*log10(IF_curve));这里s_st是固定的空间导向矢量，扫描时只改变多普勒维。w_red和R_red都来自局域 STAP，因此曲线凹口位置对应的是降维后的杂波范围。如果凹口不够深，返回上一步调整对角加载；如果凹口位置偏移，则需要检查多普勒通道索引doppler_ref是否正确。
5.2 用 publish 一键导出 HTML 图并与 fig68/fig69 对照
每次修改参数后，重新跑主脚本并手动截图效率太低。代码包里已有的 fig48、fig69 等 HTML 页面，是 MATLAB publish 生成的。可以在主脚本末尾加上 publish 命令，自动把运行结果输出成同名格式，方便和原始包对比。
% 一键导出 HTML 结果到 html 目录 opts = struct('format', 'html', 'outputDir', './html'); publish('stap_main.m', opts);运行之后检查新生成的 fig68.html 和 fig69.html。如果凹口深度与原始包差异超过 2 dB，优先怀疑训练样本集合、保护距离门数量或者对角加载系数。最后把rcond(R_hat)和改善因子凹口深度两行记录追加到一个日志文件里，保存为一列基线数据，下次调参时先看这两列，再决定要不要动对角加载。
本文还有配套的精品资源，点击获取
