简介面向本科、硕士阶段雷达信号处理教研学习的一套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 个脉冲第 N1 到 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 的影响阵元数N8决定空间自由度和波束宽度脉冲数M16决定多普勒分辨率和时间自由度载机速度v100 m/s改变杂波多普勒展宽程度阵元间距dλ/2影响空间采样和无模糊范围PRFfr2000 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 的样本需求在机载环境中不现实以 N8、M16 为例全维 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)*Nb_dopp选出来的是另一个子空间结果表现会完全失真。这也是很多教材代码跑不出理想凹口的原因之一。降维后只有 9 个自由度训练样本需求从 256 降到 18 左右在单帧数据里基本可以满足。全维和局域 STAP 的差异对比如下指标全维 STAP局域 3×3 STAP自由度N*M1289训练样本需求≥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)和改善因子凹口深度两行记录追加到一个日志文件里保存为一列基线数据下次调参时先看这两列再决定要不要动对角加载。本文还有配套的精品资源点击获取