两三年前我第一次在MATLAB里跑通Chirp Scaling算法时盯着屏幕上那九个清晰的点目标亮斑愣了好一会儿。当时距离多普勒RD算法也能成像但每次做距离徙动校正都绕不开插值慢不说插值核选不好图像上还会冒出细小伪影。而Chirp Scaling算法把整条处理链化简成几次FFT和相位相乘思路干净实现也优雅。这篇就是围绕“用MATLAB把合成孔径雷达CS算法从回波仿真到成像完整做一遍”写的里面包括了参数设计、回波构造、CS算法每一步的代码以及我调试时踩过的坑。适合正在啃SAR成像教材、准备课程设计或者想快速验证自己理解的读者。1. 从RD到CS为什么我觉得CS值得吃透1.1 RD算法那片“插值乌云”先稍微回顾一下距离多普勒算法。RD算法把原始回波做方位向FFT之后得到的是距离-多普勒域。在这个域里不同斜距处的目标其距离徙动曲线的弯曲程度不一样也就是所谓“距离单元徙动是随距离变化的”。RCMC这一步要在距离-多普勒域把弯曲的轨迹一路平移回去常规做法是插值要么sinc插值要么8点立方卷积要么频域补零处理。插值在仿真里不是不能跑但放在真实回波上会让人头疼。插值核一长运算量成倍上涨核函数截断又会带来近似误差直接表现为旁瓣抬升、相位残差甚至图像里出现细小的“鬼影”。我在学习阶段一直有个疑问为什么不能把不同距离上的徙动曲线统一成同一条再校正Chirp Scaling算法回答的就是这个问题。1.2 CS算法把“曲线校直”变成了两次缩放CS算法最核心的一点是利用线性调频信号在时延轴上的尺度不变性。雷达发射的LFM信号乘上某个二次相位因子之后只是包络位置被拉伸或压缩信号本身仍是线性调频的。这样我们就能在距离-多普勒域中构造一个与多普勒频率相关的“尺度因子”把不同斜距目标各自的徙动曲线斜率统一拉到一个参考值上。这个过程听起来玄用生活里的比喻就很简单大家排队进站本来每个人走的通道宽窄不同CS算法相当于在入口处放了一台可调的“宽度处理器”先把每个人的步伐调成一致再统一带队通过。后面距离向FFT在距离频域里乘一个匹配滤波器一次同时完成距离压缩、二次距离压缩和一致徙动校正。全程没有插值全是相位相乘。所以在MATLAB里看CS算法它给人的直观感受就是那种特别规整的矩阵运算二维回波矩阵在距离向和方位向反复变换每次变换都乘上一个大小为“距离×方位”的相位矩阵。你不需要像RD算法那样去挪动量数据只需要把相位构造对结果自然就聚焦了。2. 仿真前先把这几个参数定死很多人在CS算法上栽跟头不是算法本身没看懂而是仿真参数从一开始就自相矛盾。比如距离向采样率小于带宽方位向脉冲重复频率不够或者场景布得太大导致回波矩阵根本覆盖不了目标。参数定死后面代码才有意义。2.1 距离向、方位向分辨率怎么算SAR两点分辨率公式是绕不开的先摆出来距离向分辨率δr c / (2·Bw)c是光速Bw是发射信号带宽。方位向分辨率δa La / 2La是天线实际孔径长度。注意这里和平台高度、距离没有直接关系这是合成孔径雷达最反直觉的地方。举个例子带宽400 MHz距离向分辨率约0.375 m。如果天线孔径长度La4 m方位向分辨率约2 m。点目标仿真的场景尺寸和间隔都应该比分辨率大一个数量级才看得出效果。我习惯让点目标之间距离至少25 m到100 m这样旁瓣不会互相干扰成像后每个峰都能独立评估。2.2 PRF、采样率和点目标布阵的选择方位向PRF的选择主要受多普勒带宽限制。正侧视条带模式下多普勒带宽近似为2·v/Lav是平台速度。比如平台速度7062 m/s、天线孔径4 m多普勒带宽约3531 Hz那PRF至少取4000 Hz以上才能躲开方位混叠。当然如果为了跑通点目标合成孔径时间短一点PRF取2000 Hz也能看到聚焦结果但严格来说要按多普勒带宽来定。距离向采样率更直接必须大于等于信号带宽。带宽400 MHz采样率建议取500 MHz或600 MHz。我一开始图省事用250 MHz结果距离向频谱混叠成像剖面上出现大量虚线折腾了半小时才反应过来是采样率问题。点目标布阵方面3×3是最经典的配置中心一个近距、远距、方位向不同位置各布几个。这样既能验证算法在场景边缘的表现又能看距离徙动校正是否彻底。间隔我通常取50 m到100 m不要挤在一起。2.3 我用的仿真参数表直接给出一组能跑通的参数方便照抄。参数数值说明载频 fc5.3 GHzC波段经典星载SAR常用带宽 Bw400 MHz决定距离分辨率脉冲宽度 Tp10 μs决定发射能量和距离覆盖距离向采样率 fs600 MHz高于带宽留有余量平台速度 v7062 m/s近似低轨卫星运行速度PRF2000 Hz满足本次点目标多普勒带宽场景中心斜距 R0_ref850 km参考距离距离向采样点数4096覆盖约1 km的距离窗方位向脉冲数1024对应约0.5 s合成孔径时间点目标布阵3×3距离间隔100 m方位间隔100 m这个表里的方位向脉冲数其实非常“教学化”真实星载SAR一景可能有几万个脉冲但仿真只是为了验证算法没必要把矩阵堆到几百兆。4096×1024的复数矩阵在普通PC上大概占用67 MB左右MATLAB处理起来也没压力。3. 构造回波3×3点目标场景的MATLAB实现CS算法要跑起来第一步不是写成像代码而是先给自己“造一个雷达回波”。这一步不能偷懒因为只有你自己知道回波里每个目标的真实距离和强度后面才能校验成像结果到底准不准。3.1 回波模拟的几何关系正侧视条带模式下雷达平台沿方位向飞行。假设第k个点目标位于最近斜距R0_k那么在方位慢时间ta时刻雷达与目标的瞬时斜距是R(ta) sqrt(R0_k^2 (v·ta - X_k)^2)其中X_k表示目标在方位向相对场景中心的位置。这个公式本质是勾股定理SAR的大部分几何关系都能简化成这个模型。回波的快时间信号在距离向是一段LFMs(tr, ta) exp(-jπKr(tr - 2R(ta)/c)^2) · exp(-j4πR(ta)/λ)第一项是距离向延迟第二项是方位向历史相位。Kr Bw/Tp也就是调频率。点目标回波模拟就是把每个目标的贡献叠加进二维矩阵。3.2 生成回波矩阵的代码与关键细节下面是生成3×3点目标回波的示例代码。我通常把距离向作为行、方位向作为列这样后续方位FFT沿dim2距离FFT沿dim1不容易错。clear; clc; c 3e8; fc 5.3e9; lambda c/fc; Bw 400e6; Tp 10e-6; Kr Bw/Tp; fs 600e6; PRF 2000; v 7062; R0_ref 850e3; Nr 4096; Na 1024; tr linspace(2*R0_ref/c - (Nr/2)/fs, 2*R0_ref/c (Nr/2)/fs, Nr).; ta (-Na/2 : Na/2-1) / PRF; [Ta, Tr] meshgrid(ta, tr); targets []; for k -1:1 for m -1:1 targets(end1,:) [R0_ref k*100, m*100]; end end S_raw zeros(Nr, Na); for k 1:size(targets,1) R0_k targets(k,1); X_k targets(k,2); R sqrt(R0_k^2 (v.*Ta - X_k).^2); tau_delay 2*R/c; S_k exp(-1j*pi*Kr*(Tr - tau_delay).^2) .* exp(-1j*4*pi*R/lambda); S_raw S_raw S_k; end这里有一个细节经常被忽略距离向快时间tr的原点应该对应什么位置。如果统一以参考斜距R0_ref的往返时延作为原点后成像得到的距离向坐标很容易和真实目标位置对应起来。我的习惯是让tr从“2R0_ref/c之前半个距离窗”开始这样场景中心的回波正好落在矩阵中间。网格矩阵Ta和Tr都会各占内存。实际工程里更推荐用单精度或直接用矩阵广播但教学演示用meshgrid看得最清楚。跑完点目标后记得逐点核对各目标的理论时延和构造位置这一步错了后面全错。4. CS成像主流程四次FFT和三次相位乘法的代码这部分是整个CS算法的核心。流程写出来其实很简短方位向FFT进入距离-多普勒域。第一次相位相乘Chirp Scaling因子。距离向FFT进入距离频域-方位多普勒域。第二次相位相乘距离压缩、二次距离压缩、一致RCMC一次完成。距离向IFFT回到距离-多普勒域。第三次相位相乘方位匹配滤波和残余相位补偿。方位向IFFT得到图像。我在这里把每一步单独拆开代码和原理绑在一起讲。4.1 第一次相位相乘Chirp Scaling因子到底在干嘛先做方位向FFT注意我在FFT前都加了ifftshiftFFT后再用fftshift把零频挪回中间。MATLAB自带的FFT基于数字频率0到2π不处理原点就会让图像劈成两半。S_azfft fftshift(fft(ifftshift(S_raw, 2), [], 2), 2);方位频率轴和多普勒域的距离徙动表达式有关fa (-Na/2 : Na/2-1) * (PRF / Na); D sqrt(1 - (lambda .* fa / (2*v)).^2); C 1 ./ D - 1; % Chirp Scaling 因子 tau_ref 2*R0_ref./(c*D); % 参考距离在距离-多普勒域的时延其中D是距离徙动因数正侧视情况下它只和多普勒频率有关。C这个量就是用来做尺度变换的。第一次相位相乘的公式写作φ_CS exp(-jπKr·C·(tr - tau_ref).^2)在代码里tr是Nr×1列向量tau_ref是1×Na行向量自动广播成Nr×Na矩阵。phase_cs exp(-1i*pi*Kr .* C .* (Tr - tau_ref).^2); S_cs S_azfft .* phase_cs;这一步做的是把不同距离目标的距离徙动曲线“斜率”调整到一致。C和D的关系是C 1/D - 1所以当D1也就是多普勒频率为0的时候C0该处相位因子退化为1不做任何改动。多普勒频率越大的目标被“缩放”得越明显。4.2 距离FFT后的距离压缩与一致徙动校正CS相位乘完后信号在距离向依然是LFM只是调频率变了。现在距离向FFT把信号变换到距离频域。ft (-Nr/2 : Nr/2-1) * (fs / Nr); [Ft, Fa] meshgrid(ft, fa); % 注意这里Fa是NaN? fa为1×Naft为1×Nr矩阵维度Na×Nr如果你跟我一样习惯距离向是行、方位向是列那么这里要小心S_cs是Nr×Na距离FFT沿dim1。用meshgrid生成网格时ft应为Nr×1fa为1×Na。下面用矩阵广播来构造相位。S_distfft fftshift(fft(ifftshift(S_cs, 1), [], 1), 1); Keff Kr ./ (1 C); % 参考距离处等效调频率1×Na phase_rcmc exp(1i*pi .* (Ft.^2) ./ Keff) ... .* exp(1i*4*pi .* Ft .* R0_ref .* C / (c.*D)); S_rcmc S_distfft .* phase_rcmc;第一个指数项完成距离压缩和二次距离压缩第二个指数项是一致距离徙动校正。很多教材里会把这几个相位拆得很细但在MATLAB实现里把它们合在一起乘是完全等价的。这个归并逻辑我在调试时验证过很多次如果你以后看到其他版本的代码所有相位函数分开写本质也是这两件事。距离压缩完成后做距离向IFFT回到距离时间域。S_t ifft(ifftshift(S_rcmc, 1), [], 1);到了这里理论上每个目标的能量都应该落在它真实斜距对应的距离单元里距离徙动已经被校正干净了。4.3 方位向FFT与匹配滤波最后一步聚焦等等方位向FFT还在最前面没动过没错CS算法全程只在方位多普勒域做完距离处理最后方位压缩仍然在方位多普勒域做。此时S_t是距离时域-方位多普勒域我们直接乘方位匹配滤波器。方位匹配滤波器要补偿两部分一是方位向历史相位二是前面Chirp Scaling操作引入的残余相位。我在初学阶段只乘第一个指数项结果点目标虽然能聚焦但旁瓣明显不对称就是因为忘了补偿残余相位。完整形式如下R0_t c * Tr / 2; % 距离时域对应的斜距 phase_az exp(1i*4*pi*R0_t.*D./lambda) ... .* exp(-1i*4*pi*Kr.*C.*(1C).*(Tr - tau_ref).^2); S_az S_t .* phase_az; img ifft(ifftshift(S_az, 2), [], 2);第一个指数项是方位压缩的匹配滤波器第二个是CS残余相位补偿。最后方位向IFFT得到二维复数图像。取模abs(img)后9个点目标应该出现在对应的距离和方位坐标上。需要特别提醒R0_t用的是距离时域对应的斜距它是二维矩阵。tau_ref是1×Na会自动广播到每个距离单元上。如果你的矩阵方向和我不同这里维度对不上MATLAB会直接报错报错信息会指向广播尺寸不一致不用慌把Tr和tau_ref转置对齐即可。5. 成像质量怎么评估峰值旁瓣比、隔离度和图像倾斜跑出图像只是第一步更重要的是知道图像到底好不好。很多教程把“图像有亮斑”就当成算法正确但真实工程里还要看旁瓣结构、相位特性、定位精度。5.1 肉眼观测与量化指标先看整幅图9个点应该在各自理论位置成像中心点强度最高边缘8个点强度一致。理想点目标响应是一个二维sinc函数距离向和方位向都有第一旁瓣。量化指标上最重要的是峰值旁瓣比。理想LFM匹配滤波后峰值旁瓣比大约为-13.26 dB。如果算法实现正确距离向和方位向剖面的PSLR都应该非常接近这个值。计算方式很简单取目标附近的最大值再找最大旁瓣值换算成分贝。[~, idx] max(abs(img(:))); [pr, pa] ind2sub(size(img), idx); profile_r abs(img(pr, :)); profile_a abs(img(:, pa)); % 峰值旁瓣比粗略计算 mainlobe max(profile_r); [~, locs] findpeaks(profile_r, SortStr, descend);在实际点目标上旁瓣会成对出现在主瓣两侧。找到除主瓣外的第一个峰值用20*log10(旁瓣/主瓣)就能得到PSLR。如果算出来在-10 dB以上要么是CS残余相位没补偿要么是距离向采样率不够。另一个容易看漏的是隔离度。3×3布阵中相邻目标之间的响应不能互相污染。假如中心目标的距离向剖面里相邻目标的旁瓣泄漏进了主瓣窗口就说明布阵间隔太小或者压缩后旁瓣过高。隔离度一般要求低于-20 dB仿真时只要间隔足够大通常不是问题。5.2 旁瓣为什么高出-13dB相位参考的锅我仿真过程中遇到过一种情况距离向剖面的PSLR看着正常方位向剖面却比-13 dB高出不少图像还有轻微“涂抹”。后来发现是方位向相位函数的参考时延和距离向网格原点没对齐。CS算法残余相位项里用到的tau_ref本质是“参考距离上点目标随多普勒变化的位置”。如果距离向快时间的零点定义和R0_ref不一致这个残余相位就算错了。修法很简单保证tr以2R0_ref/c为中心构造同时在方位压缩前用R0_t c*Tr/2重新构建立体坐标而不是直接用场景中心R0_ref代全部距离单元。这类问题最麻烦的地方在于不会是“完全不成像”而是“像差一点”。你能看到亮点但亮点的旁瓣结构是畸形的。所以我的经验是点目标仿真不要只看图像一定要把距离向和方位向剖面画出来和-13.26 dB这条理论线做对比。5.3 为什么我的图是斜的点目标图像理论上应该在行列方向上呈“十”字形。如果倾斜45度通常说明距离徙动校正没有做干净目标能量沿倾斜轨迹扩散了。出现这种情况时优先检查C因子正侧视模式下C 1/D - 1如果斜视模式下用了正侧视的D表达式距离徙动轨迹会有偏差图上就斜了。另外方位向FFT前的ifftshift如果记错顺序图像会沿方位向有半个矩阵的偏移看起来像是所有目标整体平移。这种情况不是斜而是错位。代码里统一用fftshift(fft(ifftshift(...)))保持对称性就不会出大错。6. 我踩过的最典型的三个坑以及排查链路最后这部分对我来说价值最大因为每一个坑都花了一个晚上才走出来。如果你正卡在某一步可以直接对照检查。6.1 坑一距离向采样率不够频谱混叠直接把焦点打散症状图像能出亮斑但距离向剖面主瓣外出现额外的高旁瓣看似随机实则周期。排查链路先画回波的距离向频谱观察频谱边缘是否突然隆起。如果频谱在±Bw/2处被硬切说明采样率低于奈奎斯特率。此时相位因子在离散频点上有混叠距离匹配滤波后能量无法集中。解决办法最粗暴也最有效把fs从400 MHz调到600 MHz甚至800 MHz。点目标仿真里距离向采样率多花一点内存换来的是省下半天排查时间。这个坑最大的教训是别因为“点目标带宽窄”就低估采样率。6.2 坑二fftshift用错位置图像横在边界上症状图像能聚焦但9个目标分别散落在图像边缘中心反而是黑的一片或者整个图像是四块对拼。排查链路检查距离向FFT和方位向FFT是否都用了fftshift(fft(ifftshift(...)))组合。MATLAB的FFT默认0频率在第一个点图像显示时想看到零频在中间就需要shift。如果只在FFT前shift不在FFT后shift或者反过来频谱原点会偏半个栅格。我的习惯是写一个简单函数包起来function Y fft_mid(X, dim) Y fftshift(fft(ifftshift(X, dim), [], dim), dim); end之后所有FFT都调用这个函数少踩一半的坑。IFFT则用对应的ifft_mid保证正逆变换自洽。6.3 坑三残余相位补偿项方向搞反症状图像聚焦良好旁瓣却不对称左边和右边旁瓣一个高一个低。排查链路重点看方位压缩公式中残余相位项的正负号。我在第4.3节给出的exp(-1i*...)是配合前面CS相位推导出来的。如果你用的是自己推导的公式正负号必须和CS相位保持“时域对偶”关系。简单验证方式正反两个符号都试一次看旁瓣对称的那个就是正确的。这个“试错法”虽然不够优雅但对快速验证很有用。理论上正负号可以从公式推出来但调试时先用对称性判断会快很多。真正的工程实现要严格跟踪每一处相位不能依赖试错但在学习阶段图像旁瓣结构就是我们最好的“示波器”。CS算法在MATLAB里的仿真实现说到底就是“构造参数-仿真回波-三次相位相乘-评估图像”这四件事。很多人一上来就抄版本代码跑通了也不知道每一步为什么要乘那个相位换个斜视角参数就崩。我自己更建议用点目标仿真扎扎实实看清每一个阶段数据流的变化。你只要在调代码时盯着距离-多普勒域的信号相位看到不同距离目标的徙动轨迹被CS因子“掰直”的那一刻这个算法就算真的学明白了。
