简介该资源是一套基于 MATLAB 的合成孔径雷达压缩感知CS成像算法实现代码面向具备一定雷达信号处理基础的研究生、科研人员与工程实践者用于解决低采样率条件下 SAR 高分辨率图像重构的学习与仿真验证问题。压缩包共 58 个文件约 351KB其中 11 个 .m 脚本与函数文件承担核心计算涵盖 soumekh 主流程、ftx 与 fty 方位向和距离向处理、iftx 与 ifty 逆变换、range 与 crange 距离压缩、stripmap 与 spotlight 两种成像模式及雷达方向图等模块46 个 .ps 文件为相关图表与流程示意另有 1 个 readme.txt 提供使用说明。已有 309 人学习关注说明该方向具备一定参考需求。读者可据此完整走通从信号采集建模、稀疏表示选择、优化问题设定到恢复算法执行的链路并结合条带式与聚束式仿真脚本调试算法深入理解压缩感知在 SAR 成像中的落地方式适合作为课程设计、课题预研与算法复现的实践素材。1. 解压 daima.zip 后先别找 cs.mftx/fty 才是合成孔径雷达 MATLAB 仿真的真正入口daima.zip 解压后你会看到一堆ftx.m、fty.m、iftx.m、ifty.m和soumekh开头的脚本。很多人想找合成孔径雷达 CS 算法 MATLAB 实现结果先被文件名劝退。这里的关键不是某个cs.m而是fty_ftx这一对方向变换函数它们决定了雷达回波的距离-方位频谱怎么搬、怎么逆也直接决定后续压缩感知观测矩阵的形状。在这个包里ftx.m与fty.m分别处理距离向快时间和方位向慢时间的一维 FFTiftx.m与ifty.m是对应的逆变换其余range.m、crange.m、stripmap.m、spotlight.m都建立在它们之上。这个组合适合两类人一类是正在复现经典 SAR 成像链路、想从匹配滤波切到稀疏重构的研究生另一类是已经在 MATLAB 里写过雷达成像但需要一份能快速替换观测算子的工程底稿。2. 压缩感知前传距离向方位向双域稀疏性与观测矩阵的设计空间2.1 SAR 回波在哪个域稀疏决定了你这套 CS 算法能不能收敛在写代码之前先明确一个问题你要恢复的“图像”到底在哪个变换域里稀疏。常规 SAR 匹配滤波输出的是地形后向散射系数在方位-距离网格上的分布这个分布本身通常不是严格稀疏的但点目标、角反射器和强散射体场景满足近似稀疏性。更关键的是回波信号被建模成二维复正弦的叠加方位向由平台运动产生多普勒历史距离向由线性调频脉冲产生距离频率。因此对条带和聚束雷达而言经过ftx/fty变换后的二维频域天然接近稀疏支撑这就是为什么很多 CS-SAR 实现直接拿傅里叶基当稀疏字典而不是像自然图像处理那样先试小波或 DCT。这里要区分两个概念稀疏基和观测矩阵。常见做法是用傅里叶字典Psi kron(F_az, F_range)表示场景到回波的映射再设计一个随机采样掩码作为观测矩阵Phi。实际 MATLAB 实现里Psi被拆成两步——先对矩阵的每一列做ftx再对每一行做fty。这个拆分不是实现细节而是物理含义距离向快时间 FFT 把脉冲时延变成距离频率方位向慢时间 FFT 把多普勒历程变成方位频率。如果两个方向的处理顺序颠倒得到的频谱仍然是正确的但与你后续要做的距离徙动校正crange.m就对不上号。在设计观测矩阵时我一般会避免显式构造一个M*N x M*N的大矩阵那会让内存先爆掉。更好的方式是写函数句柄A (x) fty(ftx(x)); % 二维频域观测算子 AH (y) iftx(ifty(y)); % 对应的伴随算子用于梯度回代这里的A接收一个Nr x Na的回波矩阵fty(ftx(x))先做距离向 FFT 再做方位向 FFT物理上等价于二维频域采样。AH是它的共轭转置效应用于迭代重构中的反向投影。注意AH必须与A严格匹配否则迭代算法会在前几步就发散。常见错误是只写iftx(ifty(y))忘了调换顺序或者把fft写成ifft少了一个尺度因子。2.2 欠采样方式、恢复算法与 MATLAB 里的矩阵表达CS 的收益来自欠采样。对 SAR 来说最实用的欠采样方式有三种距离向随机脉冲抽取、方位向随机孔径缺失、二维频域随机掩码。三者对应不同的物理约束也对应不同的观测矩阵写法。欠采样方式观测矩阵写法MATLAB 风格适用场景恢复难点距离向随机脉冲mask rand(size(s)) 0.3后与A(x)逐点相乘发射脉冲受限距离向旁瓣抬高方位向随机孔径对矩阵的列随机置零航线不连续多普勒模糊二维频域掩码mask zeros(N,M); mask(randperm(N*M, K))1;数据存储后处理观测矩阵自洽性下面这段代码是二维频域掩码下的迭代软阈值骨架适合先用小矩阵验证链路x zeros(128, 128); x(64, 64) 1 1i; % 稀疏点目标 x(32, 96) 0.5; % 第二个弱散射点 mask zeros(size(x)); Pn randperm(numel(x), round(0.3 * numel(x))); mask(Pn) 1; % 30% 频域采样点 y mask .* fty(ftx(x)); % 欠采样观测 xs zeros(size(x)); % 初始解 mu 0.1; % 步长过大会震荡 for k 1:200 grad iftx(ifty(mask .* (mask .* fty(ftx(xs)) - y))); xs xs - mu * grad; xs max(abs(xs) - 0.01, 0) .* exp(1i * angle(xs)); % 复软阈值 end imagesc(abs(xs)); title(IHT 重构结果);这里mask就是观测矩阵的载体fty(ftx(xs))每次都走一次完整二维 FFT。iftx(ifty(...))是反向投影把频域残差拉回图像域。max(abs(xs)-0.01,0)是软阈值收缩0.01是稀疏度阈值场景越稀疏可以设得越大。这段代码不是最终优化版本但能验证你的ftx/fty方向约定是否正确如果重构出来是一条对角线而不是两个点说明伴随算子与正变换不一致先回头检查 FFT 的fftshift位置。3. stripmap/spotlight 里的 ftx/fty 实现函数约定、参数边界与常见错位3.1 四个方向变换函数的实现约定这个压缩包里的ftx.m、fty.m、iftx.m、ifty.m看起来像四行封装却是整个 SAR 仿真的地基。不同作者的代码里ftx到底沿行还是沿列 FFT 并不统一在 Soumekh 的体系里通常把矩阵的每一行当成一个距离向快时间采样序列每一列当成一个方位向慢时间采样序列。对应的典型实现是function Sx ftx(s) % 距离向快时间 FFT沿矩阵列方向对每个方位单元做 FFT % s : Nr x NaNr 为距离采样点数Na 为方位采样点数 Sx fftshift(fft(fftshift(s, 1), [], 1), 1); endfunction Sy fty(s) % 方位向慢时间 FFT沿矩阵行方向对每个距离单元做 FFT Sy fftshift(fft(fftshift(s, 2), [], 2), 2); end这里fftshift(s,1)先把距离向零频移到序列中央fft(...,[],1)按列做 FFT之后fftshift(...,1)再搬回来得到零频居中的距离向频谱。fty同理只是维度从 1 换成了 2。iftx与ifty应当把fft换成ifft同时fftshift要成对出现。我见过的最隐蔽错误是fftshift只做了一次正变换时搬一次逆变换时又搬一次结果图像中心翻转 180 度。验证方向约定的最快方式不是读注释而是构造一个只在中心列有值的矩阵s zeros(64, 64); s(32, 32) 1; % 点目标 Sx ftx(s); Sy fty(s); figure; subplot(121); imagesc(abs(Sx)); axis image; title(ftx 结果); subplot(122); imagesc(abs(Sy)); axis image; title(fty 结果);如果ftx是距离向 FFT点目标在距离向会铺开成一条竖线如果fty是方位向 FFT点目标会铺开成横线。看到结果再下结论比查谁调用谁更靠谱。3.2 strip map 与 spotlight 主脚本里的调用关系stripmap.m和spotlight.m是这个压缩包里最值得读的两个入口。stripmap.m负责连续条带场景的回波仿真与成像spotlight.m负责聚束模式下对同一区域持续照射的处理。两者的数据矩阵形状不同条带模式的方位向采样次数与合成孔径长度成正比聚束模式则在数据录取过程中始终保持波束指向同一中心。因此stripmap.m里调用fty做方位向压缩时需要配合平台速度V与脉冲重复频率PRF计算多普勒频移而spotlight.m会先调用crange.m做距离徙动校正再进入二维频域聚焦。压缩包里的rad_pat.m通常承担天线方向图加权sig_sub_a.m和sig_sub_b.m则常用于分段子孔径或子带回波生成它们都会在某个环节回调ftx/fty。常见的参数维度关系如下表这个表在排错时比代码注释更管用参数符号量纲在 MATLAB 里的常见命中位置距离采样率FsHzrange.m的时延轴t (0:N-1)/Fs载频fcHzfty的方位向频率轴生成平台速度Vm/sstripmap.m的多普勒相位项脉冲宽度Tpsrange.m的矩形包络判断我在处理这类脚本时习惯先把stripmap.m里的主循环找出来而不是看文件头注释。soumekh作为主脚本通常会调用stripmap生成回波再调用ftx/fty做匹配滤波成像。如果你想改成 CS 流程正确做法是保留stripmap里生成回波的部分删掉它后面对fty/ftx的匹配滤波调用把回波矩阵交给自己的稀疏恢复函数。注意不要直接替换stripmap内部某一行因为你可能把它内部的参考函数也一起换掉了。4. 把 CS 恢复接到 soumekh 数据链路上一套可复现的 MATLAB 改造流程4.1 解压后第一步建立方向变换的回归基线拿到daima.zip后我一般不会立刻跑soumekh而是先建一条回归基线确保ftx/fty/iftx/ifty这四个函数在你机器上行为一致。先解压到D:\radar_lab\soumekh_cs然后在 MATLAB 里执行addpath(D:\radar_lab\soumekh_cs); s0 randn(64, 64) 1i * randn(64, 64); % 复随机矩阵 s1 iftx(ifty(fty(ftx(s0)))); fprintf(重建误差: %.2e\n, norm(s1 - s0, fro) / norm(s0, fro));理想情况下误差在1e-14量级。如果误差接近1e-1大概率是这四个函数里混入了真正的fft归一化因子或者fty/ftx的方向反了。这一步必须通过后面所有 CS 实验才有意义。接下来用stripmap.m生成一组小回波把它的输出记为s_echo并把参考成像结果压缩包里对应的P*.ps图保存成.mat文件作为回归基准。之后每改一次代码都重跑一次对比。4.2 重建回波域的观测矩阵与恢复脚本CS 重构的目标是从少量回波采样y中恢复完整图像x。这里的关键是让观测算子A与你手里的回波数据维度对齐。stripmap.m输出的回波矩阵一般是Nr x Na即每一列是一个方位向采样脉冲每一行是一个距离向快时间窗。那么观测算子就定义为Nr size(s_echo, 1); Na size(s_echo, 2); A (z) fty(ftx(z)); % 正变换图像 - 二维频域 AH (r) iftx(ifty(r)); % 伴随频域残差 - 图像 mask rand(Nr, Na) 0.25; % 保留 75% 频域采样 y mask .* A(s_echo); % 欠采样后的观测数据这里mask在rand 0.25条件下保留了 75% 的数据对应 25% 的欠采样率。你可以把0.25改成0.5体验更极端的压缩感知场景。注意mask必须在每次实验中固定否则观测矩阵不满足约束等距性RIP重构结果会出现随机闪烁。一般做法是在实验最前面用rng(2024)固定随机种子。有了观测数据下一步是求解下面的稀疏优化问题lambda 0.05; % 稀疏正则权重 xC zeros(Nr, Na); % 初始图像 for iter 1:100 res mask .* A(xC) - y; % 频域残差 grad AH(res); % 反投影到图像域 xC xC - 0.02 * grad; % 梯度下降 % 软阈值每个像素做复数软阈值收缩 xC sign(xC) .* max(abs(xC) - lambda, 0); end这个循环就是最简单的迭代软阈值算法ISTA。res是当前估计的频域投影与真实欠采样数据之差AH(res)把频域误差搬回图像域0.02是步长过大会变成发散过小则收敛太慢。sign(xC) .* max(abs(xC) - lambda, 0)是对复数矩阵做逐点软阈值其中lambda控制稀疏性场景越复杂越应该调大但调得过大时弱散射点会被直接置零在最终图像里表现为目标缺失。4.3 从条带模式切到聚束模式的差异点spotlight.m的接入套路与条带模式几乎一样但有一个必须先处理的额外模块距离徙动校正RMC。聚束模式下目标回波的距离弯曲量随方位角度变化同一目标的回波会跨越多个距离门。如果不先调用crange.m校正二维频域里的目标能量会沿着距离方向拖尾CS 重构时无法用有限个稀疏系数描述它。所以在聚束模式下我会把重构链改成s_rmc crange(s_echo, rng_axis, params); % 先做距离徙动校正 y mask .* fty(ftx(s_rmc)); % 再进入频域观测crange.m输出的矩阵s_rmc与s_echo尺寸一致。如果你在聚束模式下跳过这一步得到的频谱会呈现出著名的“弯脊椎”形状后续只能靠增大稀疏度阈值来掩盖错误获得的图像分辨率会明显下降。另外spotlight.m的方位向有效孔径长度与目标位置有关在生成观测掩码时不要用条带模式的整个矩阵维度而应根据目标方位区间把掩码裁剪到有效孔径区间内。5. 用 P*.ps 比对图反推重构质量CS-SAR 排错与参数校调三招5.1 把 .ps 图转成 MATLAB 可读格式建基准压缩包里那些P1.1a.ps、P3.7a.ps文件是 Soumekh 书里的仿真结果图不能直接被 MATLAB 读取。我一般用ps2pdf转 PDF 后截图或者用pstoedit转 PNG再存进ref_imgs/目录。这样做的目的是在改代码后对比图像的主瓣宽度CS 重构如果正确点目标主瓣宽度应当接近匹配滤波的参考图如果旁瓣出现非对称优先怀疑观测矩阵的方向约定而不是去调正则参数。5.2 从二维频谱中心线的能量残差判断方向约定当重构结果出现整幅图像旋转 90 度或上下翻转时用imagesc看频谱就能定位。ftx和fty如果各少一次fftshift重构图像会翻转如果两个函数方向互换重构图像会转置。我的排查方法是找一条只在某一行有值的原始回波做一次正反变换后看峰值是否回到原行。下面这段自洽性断言是每个 main 脚本里都应该保留的assert(norm(ifty(iftx(fty(ftx(x0)))) - x0, fro) ... 1e-10, 正逆变换不自洽检查 fftshift 或用错维度);norm计算 Frobenius 范数1e-10是一个经验阈值如果失败先把assert后面的报错信息改成你正在查的函数名再运行能直接定位是哪个文件出了问题。5.3 欠采样率与稀疏阈值的配合最后一条技巧是关于参数配合的欠采样率越低lambda就要越接近0.5 * max(abs(xC(:)))否则弱散射点会被整体吞掉。常用做法是先用 full 数据跑一次匹配滤波统计图像幅度的 90% 分位点把lambda设成该值的十分之一。这个经验值在Nr/Na从 128 到 1024 的范围内都稳定但当目标数量接近网格点数 10% 时稀疏假设已经失效此时应该改用分块稀疏重构不要再指望全局软阈值。在main_cs.m顶部写上rng(2024)并保存mask.mat后续无论怎么调参观测算子都不会变.ps参考图对比也才有意义。本文还有配套的精品资源点击获取 
长河君 · 2026/9/16 22:24:15 · 阅读 1.2k SAR压缩感知重构的MATLAB实现：从ftx/fty到稀疏重建
简介该资源是一套基于 MATLAB 的合成孔径雷达压缩感知CS成像算法实现代码面向具备一定雷达信号处理基础的研究生、科研人员与工程实践者用于解决低采样率条件下 SAR 高分辨率图像重构的学习与仿真验证问题。压缩包共 58 个文件约 351KB其中 11 个 .m 脚本与函数文件承担核心计算涵盖 soumekh 主流程、ftx 与 fty 方位向和距离向处理、iftx 与 ifty 逆变换、range 与 crange 距离压缩、stripmap 与 spotlight 两种成像模式及雷达方向图等模块46 个 .ps 文件为相关图表与流程示意另有 1 个 readme.txt 提供使用说明。已有 309 人学习关注说明该方向具备一定参考需求。读者可据此完整走通从信号采集建模、稀疏表示选择、优化问题设定到恢复算法执行的链路并结合条带式与聚束式仿真脚本调试算法深入理解压缩感知在 SAR 成像中的落地方式适合作为课程设计、课题预研与算法复现的实践素材。1. 解压 daima.zip 后先别找 cs.mftx/fty 才是合成孔径雷达 MATLAB 仿真的真正入口daima.zip 解压后你会看到一堆ftx.m、fty.m、iftx.m、ifty.m和soumekh开头的脚本。很多人想找合成孔径雷达 CS 算法 MATLAB 实现结果先被文件名劝退。这里的关键不是某个cs.m而是fty_ftx这一对方向变换函数它们决定了雷达回波的距离-方位频谱怎么搬、怎么逆也直接决定后续压缩感知观测矩阵的形状。在这个包里ftx.m与fty.m分别处理距离向快时间和方位向慢时间的一维 FFTiftx.m与ifty.m是对应的逆变换其余range.m、crange.m、stripmap.m、spotlight.m都建立在它们之上。这个组合适合两类人一类是正在复现经典 SAR 成像链路、想从匹配滤波切到稀疏重构的研究生另一类是已经在 MATLAB 里写过雷达成像但需要一份能快速替换观测算子的工程底稿。2. 压缩感知前传距离向方位向双域稀疏性与观测矩阵的设计空间2.1 SAR 回波在哪个域稀疏决定了你这套 CS 算法能不能收敛在写代码之前先明确一个问题你要恢复的“图像”到底在哪个变换域里稀疏。常规 SAR 匹配滤波输出的是地形后向散射系数在方位-距离网格上的分布这个分布本身通常不是严格稀疏的但点目标、角反射器和强散射体场景满足近似稀疏性。更关键的是回波信号被建模成二维复正弦的叠加方位向由平台运动产生多普勒历史距离向由线性调频脉冲产生距离频率。因此对条带和聚束雷达而言经过ftx/fty变换后的二维频域天然接近稀疏支撑这就是为什么很多 CS-SAR 实现直接拿傅里叶基当稀疏字典而不是像自然图像处理那样先试小波或 DCT。这里要区分两个概念稀疏基和观测矩阵。常见做法是用傅里叶字典Psi kron(F_az, F_range)表示场景到回波的映射再设计一个随机采样掩码作为观测矩阵Phi。实际 MATLAB 实现里Psi被拆成两步——先对矩阵的每一列做ftx再对每一行做fty。这个拆分不是实现细节而是物理含义距离向快时间 FFT 把脉冲时延变成距离频率方位向慢时间 FFT 把多普勒历程变成方位频率。如果两个方向的处理顺序颠倒得到的频谱仍然是正确的但与你后续要做的距离徙动校正crange.m就对不上号。在设计观测矩阵时我一般会避免显式构造一个M*N x M*N的大矩阵那会让内存先爆掉。更好的方式是写函数句柄A (x) fty(ftx(x)); % 二维频域观测算子 AH (y) iftx(ifty(y)); % 对应的伴随算子用于梯度回代这里的A接收一个Nr x Na的回波矩阵fty(ftx(x))先做距离向 FFT 再做方位向 FFT物理上等价于二维频域采样。AH是它的共轭转置效应用于迭代重构中的反向投影。注意AH必须与A严格匹配否则迭代算法会在前几步就发散。常见错误是只写iftx(ifty(y))忘了调换顺序或者把fft写成ifft少了一个尺度因子。2.2 欠采样方式、恢复算法与 MATLAB 里的矩阵表达CS 的收益来自欠采样。对 SAR 来说最实用的欠采样方式有三种距离向随机脉冲抽取、方位向随机孔径缺失、二维频域随机掩码。三者对应不同的物理约束也对应不同的观测矩阵写法。欠采样方式观测矩阵写法MATLAB 风格适用场景恢复难点距离向随机脉冲mask rand(size(s)) 0.3后与A(x)逐点相乘发射脉冲受限距离向旁瓣抬高方位向随机孔径对矩阵的列随机置零航线不连续多普勒模糊二维频域掩码mask zeros(N,M); mask(randperm(N*M, K))1;数据存储后处理观测矩阵自洽性下面这段代码是二维频域掩码下的迭代软阈值骨架适合先用小矩阵验证链路x zeros(128, 128); x(64, 64) 1 1i; % 稀疏点目标 x(32, 96) 0.5; % 第二个弱散射点 mask zeros(size(x)); Pn randperm(numel(x), round(0.3 * numel(x))); mask(Pn) 1; % 30% 频域采样点 y mask .* fty(ftx(x)); % 欠采样观测 xs zeros(size(x)); % 初始解 mu 0.1; % 步长过大会震荡 for k 1:200 grad iftx(ifty(mask .* (mask .* fty(ftx(xs)) - y))); xs xs - mu * grad; xs max(abs(xs) - 0.01, 0) .* exp(1i * angle(xs)); % 复软阈值 end imagesc(abs(xs)); title(IHT 重构结果);这里mask就是观测矩阵的载体fty(ftx(xs))每次都走一次完整二维 FFT。iftx(ifty(...))是反向投影把频域残差拉回图像域。max(abs(xs)-0.01,0)是软阈值收缩0.01是稀疏度阈值场景越稀疏可以设得越大。这段代码不是最终优化版本但能验证你的ftx/fty方向约定是否正确如果重构出来是一条对角线而不是两个点说明伴随算子与正变换不一致先回头检查 FFT 的fftshift位置。3. stripmap/spotlight 里的 ftx/fty 实现函数约定、参数边界与常见错位3.1 四个方向变换函数的实现约定这个压缩包里的ftx.m、fty.m、iftx.m、ifty.m看起来像四行封装却是整个 SAR 仿真的地基。不同作者的代码里ftx到底沿行还是沿列 FFT 并不统一在 Soumekh 的体系里通常把矩阵的每一行当成一个距离向快时间采样序列每一列当成一个方位向慢时间采样序列。对应的典型实现是function Sx ftx(s) % 距离向快时间 FFT沿矩阵列方向对每个方位单元做 FFT % s : Nr x NaNr 为距离采样点数Na 为方位采样点数 Sx fftshift(fft(fftshift(s, 1), [], 1), 1); endfunction Sy fty(s) % 方位向慢时间 FFT沿矩阵行方向对每个距离单元做 FFT Sy fftshift(fft(fftshift(s, 2), [], 2), 2); end这里fftshift(s,1)先把距离向零频移到序列中央fft(...,[],1)按列做 FFT之后fftshift(...,1)再搬回来得到零频居中的距离向频谱。fty同理只是维度从 1 换成了 2。iftx与ifty应当把fft换成ifft同时fftshift要成对出现。我见过的最隐蔽错误是fftshift只做了一次正变换时搬一次逆变换时又搬一次结果图像中心翻转 180 度。验证方向约定的最快方式不是读注释而是构造一个只在中心列有值的矩阵s zeros(64, 64); s(32, 32) 1; % 点目标 Sx ftx(s); Sy fty(s); figure; subplot(121); imagesc(abs(Sx)); axis image; title(ftx 结果); subplot(122); imagesc(abs(Sy)); axis image; title(fty 结果);如果ftx是距离向 FFT点目标在距离向会铺开成一条竖线如果fty是方位向 FFT点目标会铺开成横线。看到结果再下结论比查谁调用谁更靠谱。3.2 strip map 与 spotlight 主脚本里的调用关系stripmap.m和spotlight.m是这个压缩包里最值得读的两个入口。stripmap.m负责连续条带场景的回波仿真与成像spotlight.m负责聚束模式下对同一区域持续照射的处理。两者的数据矩阵形状不同条带模式的方位向采样次数与合成孔径长度成正比聚束模式则在数据录取过程中始终保持波束指向同一中心。因此stripmap.m里调用fty做方位向压缩时需要配合平台速度V与脉冲重复频率PRF计算多普勒频移而spotlight.m会先调用crange.m做距离徙动校正再进入二维频域聚焦。压缩包里的rad_pat.m通常承担天线方向图加权sig_sub_a.m和sig_sub_b.m则常用于分段子孔径或子带回波生成它们都会在某个环节回调ftx/fty。常见的参数维度关系如下表这个表在排错时比代码注释更管用参数符号量纲在 MATLAB 里的常见命中位置距离采样率FsHzrange.m的时延轴t (0:N-1)/Fs载频fcHzfty的方位向频率轴生成平台速度Vm/sstripmap.m的多普勒相位项脉冲宽度Tpsrange.m的矩形包络判断我在处理这类脚本时习惯先把stripmap.m里的主循环找出来而不是看文件头注释。soumekh作为主脚本通常会调用stripmap生成回波再调用ftx/fty做匹配滤波成像。如果你想改成 CS 流程正确做法是保留stripmap里生成回波的部分删掉它后面对fty/ftx的匹配滤波调用把回波矩阵交给自己的稀疏恢复函数。注意不要直接替换stripmap内部某一行因为你可能把它内部的参考函数也一起换掉了。4. 把 CS 恢复接到 soumekh 数据链路上一套可复现的 MATLAB 改造流程4.1 解压后第一步建立方向变换的回归基线拿到daima.zip后我一般不会立刻跑soumekh而是先建一条回归基线确保ftx/fty/iftx/ifty这四个函数在你机器上行为一致。先解压到D:\radar_lab\soumekh_cs然后在 MATLAB 里执行addpath(D:\radar_lab\soumekh_cs); s0 randn(64, 64) 1i * randn(64, 64); % 复随机矩阵 s1 iftx(ifty(fty(ftx(s0)))); fprintf(重建误差: %.2e\n, norm(s1 - s0, fro) / norm(s0, fro));理想情况下误差在1e-14量级。如果误差接近1e-1大概率是这四个函数里混入了真正的fft归一化因子或者fty/ftx的方向反了。这一步必须通过后面所有 CS 实验才有意义。接下来用stripmap.m生成一组小回波把它的输出记为s_echo并把参考成像结果压缩包里对应的P*.ps图保存成.mat文件作为回归基准。之后每改一次代码都重跑一次对比。4.2 重建回波域的观测矩阵与恢复脚本CS 重构的目标是从少量回波采样y中恢复完整图像x。这里的关键是让观测算子A与你手里的回波数据维度对齐。stripmap.m输出的回波矩阵一般是Nr x Na即每一列是一个方位向采样脉冲每一行是一个距离向快时间窗。那么观测算子就定义为Nr size(s_echo, 1); Na size(s_echo, 2); A (z) fty(ftx(z)); % 正变换图像 - 二维频域 AH (r) iftx(ifty(r)); % 伴随频域残差 - 图像 mask rand(Nr, Na) 0.25; % 保留 75% 频域采样 y mask .* A(s_echo); % 欠采样后的观测数据这里mask在rand 0.25条件下保留了 75% 的数据对应 25% 的欠采样率。你可以把0.25改成0.5体验更极端的压缩感知场景。注意mask必须在每次实验中固定否则观测矩阵不满足约束等距性RIP重构结果会出现随机闪烁。一般做法是在实验最前面用rng(2024)固定随机种子。有了观测数据下一步是求解下面的稀疏优化问题lambda 0.05; % 稀疏正则权重 xC zeros(Nr, Na); % 初始图像 for iter 1:100 res mask .* A(xC) - y; % 频域残差 grad AH(res); % 反投影到图像域 xC xC - 0.02 * grad; % 梯度下降 % 软阈值每个像素做复数软阈值收缩 xC sign(xC) .* max(abs(xC) - lambda, 0); end这个循环就是最简单的迭代软阈值算法ISTA。res是当前估计的频域投影与真实欠采样数据之差AH(res)把频域误差搬回图像域0.02是步长过大会变成发散过小则收敛太慢。sign(xC) .* max(abs(xC) - lambda, 0)是对复数矩阵做逐点软阈值其中lambda控制稀疏性场景越复杂越应该调大但调得过大时弱散射点会被直接置零在最终图像里表现为目标缺失。4.3 从条带模式切到聚束模式的差异点spotlight.m的接入套路与条带模式几乎一样但有一个必须先处理的额外模块距离徙动校正RMC。聚束模式下目标回波的距离弯曲量随方位角度变化同一目标的回波会跨越多个距离门。如果不先调用crange.m校正二维频域里的目标能量会沿着距离方向拖尾CS 重构时无法用有限个稀疏系数描述它。所以在聚束模式下我会把重构链改成s_rmc crange(s_echo, rng_axis, params); % 先做距离徙动校正 y mask .* fty(ftx(s_rmc)); % 再进入频域观测crange.m输出的矩阵s_rmc与s_echo尺寸一致。如果你在聚束模式下跳过这一步得到的频谱会呈现出著名的“弯脊椎”形状后续只能靠增大稀疏度阈值来掩盖错误获得的图像分辨率会明显下降。另外spotlight.m的方位向有效孔径长度与目标位置有关在生成观测掩码时不要用条带模式的整个矩阵维度而应根据目标方位区间把掩码裁剪到有效孔径区间内。5. 用 P*.ps 比对图反推重构质量CS-SAR 排错与参数校调三招5.1 把 .ps 图转成 MATLAB 可读格式建基准压缩包里那些P1.1a.ps、P3.7a.ps文件是 Soumekh 书里的仿真结果图不能直接被 MATLAB 读取。我一般用ps2pdf转 PDF 后截图或者用pstoedit转 PNG再存进ref_imgs/目录。这样做的目的是在改代码后对比图像的主瓣宽度CS 重构如果正确点目标主瓣宽度应当接近匹配滤波的参考图如果旁瓣出现非对称优先怀疑观测矩阵的方向约定而不是去调正则参数。5.2 从二维频谱中心线的能量残差判断方向约定当重构结果出现整幅图像旋转 90 度或上下翻转时用imagesc看频谱就能定位。ftx和fty如果各少一次fftshift重构图像会翻转如果两个函数方向互换重构图像会转置。我的排查方法是找一条只在某一行有值的原始回波做一次正反变换后看峰值是否回到原行。下面这段自洽性断言是每个 main 脚本里都应该保留的assert(norm(ifty(iftx(fty(ftx(x0)))) - x0, fro) ... 1e-10, 正逆变换不自洽检查 fftshift 或用错维度);norm计算 Frobenius 范数1e-10是一个经验阈值如果失败先把assert后面的报错信息改成你正在查的函数名再运行能直接定位是哪个文件出了问题。5.3 欠采样率与稀疏阈值的配合最后一条技巧是关于参数配合的欠采样率越低lambda就要越接近0.5 * max(abs(xC(:)))否则弱散射点会被整体吞掉。常用做法是先用 full 数据跑一次匹配滤波统计图像幅度的 90% 分位点把lambda设成该值的十分之一。这个经验值在Nr/Na从 128 到 1024 的范围内都稳定但当目标数量接近网格点数 10% 时稀疏假设已经失效此时应该改用分块稀疏重构不要再指望全局软阈值。在main_cs.m顶部写上rng(2024)并保存mask.mat后续无论怎么调参观测算子都不会变.ps参考图对比也才有意义。本文还有配套的精品资源点击获取