简介压缩感知CS理论为机载合成孔径雷达SAR成像提供了低速率采样与稀疏重建的新方案。MATLAB算法源码面向雷达信号处理、遥感图像反演方向的工程师与研究生针对传统SAR高采样率带来的存储、传输与实时处理难题适合作为教学演示与科研验证基础。压缩包采用zip格式内含1个.m文件大小仅2KB代码精炼规范模块划分清晰便于嵌入和修改。已有846人学习使用。文件内容将CS成像流程细化为回波信号构建、测量矩阵设计、稀疏基选择、基于L1范数的凸优化恢复、以及重构图像质量评估等完整环节并配有参数注释使用者可以据此掌握CS-SAR成像的算法骨架快速对比不同稀疏基与重构参数对图像质量的影响进而迁移至面目标、地物场景等更贴近实际的仿真分析中。这种实现方式还有助于降低数据采集成本、提升系统实时性为后续研究提供了可复现的基线参考。1. 机载SAR成像为什么把CS算法当作高分辨率标配机载SAR一旦把分辨率推到米级以内RD距离-多普勒算法里逐距离单元做插值那一步就会变成整条处理链最沉的部分。CS算法Chirp Scaling Algorithm线性调频变标换了一条路它不插值而是在距离-多普勒域用一个chirp信号和回波相乘把不同距离单元各自不同的距离徙动轨迹统一成同一条参考轨迹再在二维频域用一次相位相乘完成一致距离徙动校正RCMC。整个成像过程只靠FFT和相位乘法距离向、方位向的相位关系保持得很干净所以后面接干涉、极化处理都不吃亏。下面从机载SAR回波的信号模型开始把CS变标因子的来龙去脉讲清楚再给出一套可以直接跑的点目标仿真Python代码最后落到机载平台上最影响成像质量的几个参数上。2. 机载SAR回波模型与CS变标的数学原理2.1 从点目标回波看距离徙动量RD算法为什么卡在插值上正侧视条带机载SAR的几何可以用两个量概括平台速度V和点目标的最近斜距R0。慢时间η时刻目标斜距为$$R(\eta)\sqrt{R_0^2(V\eta-X_0)^2}$$发射信号是线性调频脉冲接收并解调到基带后的回波写成$$s(\tau,\eta)\sigma,\text{rect}!\left(\frac{\tau-\tau_d(\eta)}{T_p}\right)\exp!\left(j\pi K_r(\tau-\tau_d(\eta))^2\right)\exp!\left(-j\frac{4\pi}{\lambda}R(\eta)\right)$$其中τ是快时间K_r是发射chirp的调频率τ_d2R(η)/c最后一项是方位向相位历程。对慢时间做傅里叶变换后利用驻定相位法可以得到距离-多普勒域中的近似表达式其距离包络峰值位置落在$$\tau(f_\eta)\frac{2R_0}{c,D(f_\eta)},\qquad D(f_\eta)\sqrt{1-\left(\frac{\lambda f_\eta}{2V}\right)^2}$$从这个式子能读出三个重要信息第一目标回波在不同方位频率上落在不同的距离单元这就是距离单元徙动第二徙动量与R0成正比越远的距离单元弯路越大第三所有距离单元的徙动曲线形状相同只是幅度随R0缩放。把包络位置换算成距离域的徙动量就是$$\Delta R(f_\eta)R_0\left(\frac{1}{D(f_\eta)}-1\right)$$RD算法的做法是沿着这条双曲线轨迹逐距离单元插值把能量拉回直线。插值核长度、方位频率点数、数据量三者相乘就是可观的计算量而且插值本身会引入相位误差对后续干涉处理不友好。这是CS算法要解决的问题。2.2 Chirp Scaling的数学本质变标因子Cs是怎么来的CS算法利用的是chirp信号的一个代数性质两个chirp相乘会得到一个新的chirp其调频率和包络中心都发生偏移。设在距离-多普勒域某信号为$\exp(j\pi K_r(\tau-\tau_1)^2)$乘上变标函数$\exp(j\pi K_r C_s(\tau-\tau_{ref})^2)$展开后得到$$\exp!\left(j\pi K_r(1C_s)(\tau-\tau_{new})^2\right)\cdot e^{j\phi_{res}}$$其中新包络中心为$$\tau_{new}\frac{\tau_1C_s\tau_{ref}}{1C_s}$$现在把$\tau_12R_0/(cD)$和$\tau_{ref}2R_{ref}/(cD)$代入。我们希望变标后的包络中心中与R0有关的部分不再含1/D因子这样后续才能做一致性校正。整理后得到变标因子$$C_s(f_\eta)\frac{1}{D(f_\eta)}-1$$这个结果很干净在各种频谱里看到CSA的“变标因子”就是多普勒域中的$1/D-1$。把C_s代回去新的包络中心变成$$\tau_{new}\frac{2R_0}{c}\frac{2R_{ref}}{c}\left(\frac{1}{D}-1\right)$$第一项只与目标最近斜距有关第二项对场景内所有目标都相同是一致项。于是每个目标的徙动曲线都被“掰成”同一条参考曲线的形状剩下的统一在二维频域一次校掉。2.3 CS算法、RD算法与ωK算法的选型对比算法核心操作距离徙动校正方式主要计算量相位保持能力典型适用场景RD距离压缩、RCMC、方位压缩逐距离单元sinc插值插值运算占大头插值误差会损失相位分辨率要求不太高的常规条带CS频域变标、频域一致RCMC无插值纯相位相乘FFT与复乘好适合干涉正侧视、中小斜视角、高分辨率ωK二维频域Stolt映射频域重采样插值Stolt插值较重高斜视适应性强大斜视、超高分辨率机载平台和星载平台放在一起比的时候CS常被优先选择还有一个原因它不需要逐行插值数据流适合流水线实现实时成像处理器里容易把时延做到确定值。但机载平台惯导误差、运动误差大CS对平台速度误差比RD更敏感这是后续参数设计要重点盯的地方。3. CSA成像流程三步走变标、一致RCMC、方位压缩3.1 第一步在距离-多普勒域做CS变标原始回波先做方位向FFT进入距离-多普勒域。这里方位频率轴f_eta要经过fftshift使频率从负到正排列后面的相位因子才能和矩阵行对应上。变标相位因子写为$$H_1(\tau,f_\eta)\exp!\left(j\pi K_r C_s(f_\eta)\left(\tau-\tau_{ref}(f_\eta)\right)^2\right)$$其中$\tau_{ref}(f_\eta)2R_{ref}/(cD(f_\eta))$是参考斜距R_ref对应的徙动轨迹。这一步做完所有距离单元的徙动曲线都变成了参考距离的曲线形状因此同时为下一步的“一致距离徙动校正”创造了条件。实际实现里K_r应该用距离压缩前测得的实际调频率而不是发射参数表里的理论值。发射机增益曲线、功放非线性都会让理想chirp与真实chirp有偏差用实测值可以让变标后的残余相位更小。3.2 第二步在二维频域做距离压缩与一致RCMC变标后的信号做距离向FFT进入二维频域。此时距离chirp的调频率已经从K_r变成了$K_{eff}K_r(1C_s)$所以距离压缩匹配滤波器的调频率也要跟着变这是CSA中容易写错的地方。距离压缩项与一致距离徙动校正项合并写为$$H_2(f_\tau,f_\eta)\exp!\left(j\frac{\pi f_\tau^2}{K_r(1C_s)}\right)\cdot\exp!\left(-j\frac{4\pi R_{ref}}{c}\left(\frac{1}{D}-1\right)f_\tau\right)$$第一项是chirp匹配滤波第二个相位项完成参考距离处的一次性RCMC。注意R_ref要取场景中心斜距取偏了会让整个场景图像整体偏移但不影响聚焦取对了图像位置才是绝对定位正确的。这个域里还有一个经常被忽略的细节距离压缩的匹配滤波器带宽要和发射信号带宽匹配。采样率通常取带宽的1.2~1.4倍匹配滤波后在频带边缘会有幅度衰减如果后续要加窗压低旁瓣这部分衰减也要计入。3.3 第三步距离IFFT后做方位压缩与残余相位补偿二维频域处理完距离向IFFT回到距离-多普勒域。信号此时已经在正确的距离单元上剩下的方位向相位仍是一个关于f_eta的二次函数。方位压缩匹配滤波器为$$H_3(f_\eta,R_0)\exp!\left(j\frac{4\pi R_0}{\lambda}D(f_\eta)\right)$$其中R_0取当前距离单元对应的斜距$R_0c\tau/2$。H1引入的变标过程还会带来一个与τ无关的残余相位需要一并补偿$$\phi_{res}\pi K_r(1-D)\cdot\frac{4(R_0-R_{ref})^2}{c^2D^2}$$这一步只影响相位不改变幅度聚焦。做幅度成像时可以偷懒不乘但后续要做干涉测量或者极化定标残余相位必须补偿掉。下面的代码片段给出三个相位因子的构造方式频域轴与矩阵轴的对应关系以注释说明# 距离-多普勒域变标因子与H1 Cs 1.0 / D - 1.0 tau_ref 2.0 * R_ref / (c * D) H1 np.exp(1j * np.pi * Kr * Cs[:, None] * (tau[None, :] - tau_ref[:, None])**2) # 二维频域距离压缩与一致RCMC K_eff Kr * (1.0 Cs)[:, None] H2 np.exp(1j * np.pi * f_tau[None, :]**2 / K_eff) H2 * np.exp(-1j * 4.0 * np.pi * R_ref * (1.0 / D[:, None] - 1.0) * f_tau[None, :] / c) # 距离-多普勒域方位压缩与残余相位补偿 R_img c * tau / 2.0 H3 np.exp(1j * 4.0 * np.pi * R_img[None, :] * D[:, None] / lam) phi_res np.pi * Kr * (1.0 - D[:, None]) * ( 4.0 * (R_img[None, :] - R_ref)**2 / (c**2 * D[:, None]**2) ) H3 * np.exp(-1j * phi_res)代码里所有乘法和指数都是逐元素操作广播时距离频率、方位频率的轴关系必须和前面FFT后的矩阵保持一致。这三个相位因子分别作用在各自的域中顺序不能换H1在距离-多普勒域H2在二维频域H3在距离-多普勒域。后面第4章会把完整流程串起来。4. 用Python复现机载SAR的CS算法点目标仿真全流程4.1 机载SAR成像仿真参数怎么定点目标仿真不是为了模拟真实雷达回波的复杂电磁环境而是为了验证算法本身。参数取值要让距离、方位分辨率在同一量级并且保证方位向采样满足要求。下面这组参数接近小型无人机或通航飞机平台搭载的X波段SAR工作状态参数符号数值说明载频fc9.6 GHzX波段波长约3.1 cm距离带宽B180 MHz理论距离分辨率约0.83 m脉冲宽度Tp5 us线性调频脉宽距离采样率Fs220 MHz带宽的1.22倍留过采样平台速度V150 m/s典型小型平台速度平台高度H8000 m中低空机载场景中心斜距R_ref12000 m决定下视角约41.8度方位天线口径L_a2 m方位分辨率约1 m脉冲重复频率PRF320 Hz方位多普勒带宽约150 Hz余量2.1倍方位采样点数N_az1024覆盖完整合成孔径距离采样点数N_rg2048距离窗约1395 m多普勒带宽按$B_a2V/L_a$估算为150 HzPRF取320 Hz有足够余量但比理论最小值只高了1.1倍左右数据率仍在机载记录设备可承受范围内。快时间轴以场景中心双程延时为中心距离窗宽要覆盖目标之间的斜距差。4.2 回波生成代码与参数初始化import numpy as np c 299792458.0 fc 9.6e9 # 载频 lam c / fc B 180e6 # 距离向带宽 Tp 5e-6 # 脉冲宽度 Kr B / Tp # 调频率 Fs 220e6 # 距离向采样率 V 150.0 # 平台速度 m/s R_ref 12000.0 # 场景中心斜距 PRF 320.0 # 点目标列表(方位向位置, 最近斜距) targets [(0.0, R_ref), (-50.0, R_ref 60.0), (50.0, R_ref - 50.0)] # 慢时间轴以合成孔径中心为0 N_az 1024 eta (np.arange(N_az) - N_az // 2) / PRF # 快时间轴以场景中心双程延时为中心 N_rg 2048 tau np.arange(N_rg) / Fs tau tau - N_rg / (2.0 * Fs) 2.0 * R_ref / c # 逐脉冲生成基带回波 echo np.zeros((N_az, N_rg), dtypecomplex) for x_t, R0 in targets: R np.sqrt(R0**2 (V * eta - x_t)**2) # 每个慢时刻的斜距 tau_d 2.0 * R / c t_diff tau[None, :] - tau_d[:, None] in_pulse (t_diff 0) (t_diff Tp) phase (np.pi * Kr * t_diff**2 - 4.0 * np.pi * R[:, None] / lam) echo np.where(in_pulse, np.exp(1j * phase), 0.0) # 方位FFT进入距离-多普勒域 S_rd np.fft.fftshift(np.fft.fft(echo, axis0), axes0) f_eta np.fft.fftshift(np.fft.fftfreq(N_az, 1.0 / PRF))回波仿真直接在快时间采样点上离散计算chirp相位省去了插值步骤速度也快。注意相位里第一个指数项是基带chirp第二个指数项是载频对应的方位相位历程两者缺一不可。目标坐标里既有方位偏移又有斜距差是为了检查CS算法在二维方向上的聚焦能力而不是只测中心点。4.3 CS成像核心代码与参数说明# 距离-多普勒域的多普勒因子 D np.sqrt(1.0 - (lam * f_eta / (2.0 * V))**2) # 第一步CS变标 Cs 1.0 / D - 1.0 tau_ref 2.0 * R_ref / (c * D) H1 np.exp(1j * np.pi * Kr * Cs[:, None] * (tau[None, :] - tau_ref[:, None])**2) S_rd S_rd * H1 # 距离FFT进入二维频域 S_fd np.fft.fftshift(np.fft.fft(S_rd, axis1), axes1) f_tau np.fft.fftshift(np.fft.fftfreq(N_rg, 1.0 / Fs)) # 第二步距离压缩 一致RCMC K_eff Kr * (1.0 Cs)[:, None] H2 np.exp(1j * np.pi * f_tau[None, :]**2 / K_eff) H2 * np.exp(-1j * 4.0 * np.pi * R_ref * (1.0 / D[:, None] - 1.0) * f_tau[None, :] / c) S_fd S_fd * H2 # 距离IFFT回到距离-多普勒域 S_rd np.fft.ifft(np.fft.ifftshift(S_fd, axes1), axis1) # 第三步方位压缩 残余相位补偿 R_img c * tau / 2.0 H3 np.exp(1j * 4.0 * np.pi * R_img[None, :] * D[:, None] / lam) phi_res np.pi * Kr * (1.0 - D[:, None]) * ( 4.0 * (R_img[None, :] - R_ref)**2 / (c**2 * D[:, None]**2) ) H3 * np.exp(-1j * phi_res) S_rd S_rd * H3 # 方位IFFT输出复图像 img np.fft.ifft(np.fft.ifftshift(S_rd, axes0), axis0)这里做的是忽略二次距离压缩SRC的简化CS版本在正侧视、窄带条件下误差很小。决定是否要补上SRC项的判据是$\lambda^2 R_0 B_a^2/(16 c^2 V^2)$乘以适当系数后是否超过四分之一相位周期对上面的参数这顶多带来零点几弧度的相位误差肉眼分辨不出来但做干涉时要按完整CSA算一遍。代码里三个H因子的符号非常重要任何一个变号都会导致方位向散焦或者距离向偏移排查时先看距离压缩是不是把chirp压缩成单个尖峰。5. 机载PRF、斜视角与多普勒中心估计CSA的三个调参要点5.1 PRF与方位模糊机载SAR的采样率选择机载平台速度通常只有几十到两百米每秒方位多普勒带宽并不大但PRF不能只按带宽的奈奎斯特下限取。方位信号的多普勒带宽为$B_a2V/L_a$对应上面的参数是150 HzPRF取320 Hz看似没问题但要注意方位天线方向图的主瓣外还有旁瓣多普勒频率超过PRF/2时会折叠进主瓣区形成方位模糊。工程上常见做法是取PRF为多普勒带宽的1.2到1.5倍同时留出多普勒中心估计误差的余量。机载平台的GPS/惯导精度有限实际航线速度、偏流角都在缓慢变化多普勒中心不可能是零。如果PRF余量留得太小估计出的多普勒中心一旦超过PRF/2就会模糊到另一边图像方位向整体倒置这个问题比散焦更难查。在实时处理链里一般会先用惯性导航数据计算预期多普勒中心再用数据估计偏差两个值相差超过10%就说明数据或惯导有一方出了问题。5.2 多普勒中心不准会怎样先估计再成像CS算法推导时默认回波在方位谱上的中心是零频。实际机载SAR由于存在偏流角或天线安装角多普勒中心不为零直接送入CSA会让点目标在方位向偏移且天线方向图加权后的信噪比下降。先用回波数据估计多普勒中心是在成像前必做的一步。最简单的是能量重心法取一个较强距离单元做方位向FFT频谱重心就是多普勒中心的粗估计from numpy.fft import fft, fftshift, fftfreq row echo[500, :].copy() # 取某个距离单元 spec np.abs(fftshift(fft(row, axis0)))**2 f_eta_tmp fftshift(fftfreq(N_az, 1.0 / PRF)) fdc_est np.sum(f_eta_tmp * spec) / np.sum(spec) # 在方位时域去斜 eta2 (np.arange(N_az) - N_az // 2) / PRF echo echo * np.exp(-1j * 2.0 * np.pi * fdc_est * eta2[:, None])估计完多普勒中心后回波在方位向乘以线性相位exp(-j2πfd_est·η)把频谱中心平移回零频。之后再做CS成像图像里各目标方位位置才是真实位置。注意能量重心法对强点目标效果好一片均匀杂波时反而会偏更稳的是用相邻方位频谱的互相关法但粗估阶段重心法够用。5.3 大斜视角下CS变标失效的边界CS变标因子$C_s1/D-1$是在正侧视几何下推导的。当斜视角增大距离-多普勒域的精确表达式里会多出斜视角相关的耦合项简化CSA的点目标响应会从两侧旁瓣不对称逐渐演化为主瓣展宽。经验边界是10到15度以内用简化CSA问题不大超过这个范围要使用包含二次距离压缩的完整CSA或者直接换ωK算法。机载平台的斜视角通常不是定值。飞机受侧风影响实际航迹和天线指向的夹角在一条航线上可能波动好几度。惯导数据更新的速率如果低于PRF这种波动会直接变成方位向相位误差CSA这类频域算法对相位误差的敏感度又比RD高。我一般会在处理链里把回波按航迹段分块每块单独估算等效速度V和斜视角再用分块后的参数做CSA成像而不是全航线用一组参数。6. 点目标响应的定量评估用IRW、PSLR、ISLR看CS成像质量6.1 为什么评估点目标响应要插值到更高采样率成像结果栅格通常只按原采样率排列主瓣峰值两边只有几个采样点直接数主瓣宽度会引入很大的量化误差。评估点目标响应的常规做法是把峰值附近的距离剖面和方位剖面各取一段用零填充IFFT插值到8倍采样率再做3dB宽度和旁瓣统计。插值后主瓣不会被栅格“切歪”PSLR和ISLR也才稳定。6.2 用Python计算IRW、PSLR和ISLRdef measure_metrics(profile, dr, zoom8): # profile: 一维复数剖面dr为采样间隔 N len(profile) p np.abs(np.fft.irfft(np.fft.rfft(profile), N * zoom))**2 peak np.argmax(p) power p / p[peak] # 3dB主瓣宽度用峰值左右第一个低于0.5的点做线性插值 left np.where(power[:peak] 0.5)[0][-1] right np.where(power[peak:] 0.5)[0][0] peak irw (right - left) / zoom * dr # 主瓣边界取第一个零点 zleft np.where(power[:peak] 1e-3)[0][-1] zright np.where(power[peak:] 1e-3)[0][0] peak sidelobe np.concatenate([power[:zleft], power[zright:]]) pslr 10 * np.log10(np.max(sidelobe)) islr 10 * np.log10( (np.sum(power) - np.sum(power[zleft:zright])) / np.sum(power[zleft:zright]) ) return irw, pslr, islr主瓣边界用第一个零点近似划分实际处理中也可以固定取峰值两侧各若干个3dB宽度作为主瓣范围。PSLR反映强目标旁边的弱目标可检测能力ISLR反映分布式目标背景的泄漏水平。对未加窗的矩形谱响应PSLR接近-13.3 dBISLR接近-10 dB加Hamming窗后PSLR会压到-42 dB左右但IRW会扩展到约1.3倍。6.3 验证CS算法实现正确性的三个信号特征检查距离剖面和方位剖面的IRW是否接近理论值。距离向理论IRW为$0.886\cdot c/(2B)$方位向理论IRW为$0.886\cdot \lambda R_0/(2V T_{obs})$差距超过5%就说明某个相位因子有误差。第二三个点目标的峰值应出现在各自设定的斜距和方位坐标上距离向整体偏移说明R_ref取偏方位向出现随距离变化的位移说明多普勒中心估计遗漏。第三残余相位补偿H3被关掉时幅度图像看不出明显变化但点目标相位沿方位向会出现抛物线状残余用干涉或相位图检查能立刻暴露问题。这三个特征都正常CS成像链基本可以放心交给后面的自动聚焦和辐射定标模块处理。本文还有配套的精品资源点击获取 机载SAR的CS成像算法：原理、仿真与参数调优
简介压缩感知CS理论为机载合成孔径雷达SAR成像提供了低速率采样与稀疏重建的新方案。MATLAB算法源码面向雷达信号处理、遥感图像反演方向的工程师与研究生针对传统SAR高采样率带来的存储、传输与实时处理难题适合作为教学演示与科研验证基础。压缩包采用zip格式内含1个.m文件大小仅2KB代码精炼规范模块划分清晰便于嵌入和修改。已有846人学习使用。文件内容将CS成像流程细化为回波信号构建、测量矩阵设计、稀疏基选择、基于L1范数的凸优化恢复、以及重构图像质量评估等完整环节并配有参数注释使用者可以据此掌握CS-SAR成像的算法骨架快速对比不同稀疏基与重构参数对图像质量的影响进而迁移至面目标、地物场景等更贴近实际的仿真分析中。这种实现方式还有助于降低数据采集成本、提升系统实时性为后续研究提供了可复现的基线参考。1. 机载SAR成像为什么把CS算法当作高分辨率标配机载SAR一旦把分辨率推到米级以内RD距离-多普勒算法里逐距离单元做插值那一步就会变成整条处理链最沉的部分。CS算法Chirp Scaling Algorithm线性调频变标换了一条路它不插值而是在距离-多普勒域用一个chirp信号和回波相乘把不同距离单元各自不同的距离徙动轨迹统一成同一条参考轨迹再在二维频域用一次相位相乘完成一致距离徙动校正RCMC。整个成像过程只靠FFT和相位乘法距离向、方位向的相位关系保持得很干净所以后面接干涉、极化处理都不吃亏。下面从机载SAR回波的信号模型开始把CS变标因子的来龙去脉讲清楚再给出一套可以直接跑的点目标仿真Python代码最后落到机载平台上最影响成像质量的几个参数上。2. 机载SAR回波模型与CS变标的数学原理2.1 从点目标回波看距离徙动量RD算法为什么卡在插值上正侧视条带机载SAR的几何可以用两个量概括平台速度V和点目标的最近斜距R0。慢时间η时刻目标斜距为$$R(\eta)\sqrt{R_0^2(V\eta-X_0)^2}$$发射信号是线性调频脉冲接收并解调到基带后的回波写成$$s(\tau,\eta)\sigma,\text{rect}!\left(\frac{\tau-\tau_d(\eta)}{T_p}\right)\exp!\left(j\pi K_r(\tau-\tau_d(\eta))^2\right)\exp!\left(-j\frac{4\pi}{\lambda}R(\eta)\right)$$其中τ是快时间K_r是发射chirp的调频率τ_d2R(η)/c最后一项是方位向相位历程。对慢时间做傅里叶变换后利用驻定相位法可以得到距离-多普勒域中的近似表达式其距离包络峰值位置落在$$\tau(f_\eta)\frac{2R_0}{c,D(f_\eta)},\qquad D(f_\eta)\sqrt{1-\left(\frac{\lambda f_\eta}{2V}\right)^2}$$从这个式子能读出三个重要信息第一目标回波在不同方位频率上落在不同的距离单元这就是距离单元徙动第二徙动量与R0成正比越远的距离单元弯路越大第三所有距离单元的徙动曲线形状相同只是幅度随R0缩放。把包络位置换算成距离域的徙动量就是$$\Delta R(f_\eta)R_0\left(\frac{1}{D(f_\eta)}-1\right)$$RD算法的做法是沿着这条双曲线轨迹逐距离单元插值把能量拉回直线。插值核长度、方位频率点数、数据量三者相乘就是可观的计算量而且插值本身会引入相位误差对后续干涉处理不友好。这是CS算法要解决的问题。2.2 Chirp Scaling的数学本质变标因子Cs是怎么来的CS算法利用的是chirp信号的一个代数性质两个chirp相乘会得到一个新的chirp其调频率和包络中心都发生偏移。设在距离-多普勒域某信号为$\exp(j\pi K_r(\tau-\tau_1)^2)$乘上变标函数$\exp(j\pi K_r C_s(\tau-\tau_{ref})^2)$展开后得到$$\exp!\left(j\pi K_r(1C_s)(\tau-\tau_{new})^2\right)\cdot e^{j\phi_{res}}$$其中新包络中心为$$\tau_{new}\frac{\tau_1C_s\tau_{ref}}{1C_s}$$现在把$\tau_12R_0/(cD)$和$\tau_{ref}2R_{ref}/(cD)$代入。我们希望变标后的包络中心中与R0有关的部分不再含1/D因子这样后续才能做一致性校正。整理后得到变标因子$$C_s(f_\eta)\frac{1}{D(f_\eta)}-1$$这个结果很干净在各种频谱里看到CSA的“变标因子”就是多普勒域中的$1/D-1$。把C_s代回去新的包络中心变成$$\tau_{new}\frac{2R_0}{c}\frac{2R_{ref}}{c}\left(\frac{1}{D}-1\right)$$第一项只与目标最近斜距有关第二项对场景内所有目标都相同是一致项。于是每个目标的徙动曲线都被“掰成”同一条参考曲线的形状剩下的统一在二维频域一次校掉。2.3 CS算法、RD算法与ωK算法的选型对比算法核心操作距离徙动校正方式主要计算量相位保持能力典型适用场景RD距离压缩、RCMC、方位压缩逐距离单元sinc插值插值运算占大头插值误差会损失相位分辨率要求不太高的常规条带CS频域变标、频域一致RCMC无插值纯相位相乘FFT与复乘好适合干涉正侧视、中小斜视角、高分辨率ωK二维频域Stolt映射频域重采样插值Stolt插值较重高斜视适应性强大斜视、超高分辨率机载平台和星载平台放在一起比的时候CS常被优先选择还有一个原因它不需要逐行插值数据流适合流水线实现实时成像处理器里容易把时延做到确定值。但机载平台惯导误差、运动误差大CS对平台速度误差比RD更敏感这是后续参数设计要重点盯的地方。3. CSA成像流程三步走变标、一致RCMC、方位压缩3.1 第一步在距离-多普勒域做CS变标原始回波先做方位向FFT进入距离-多普勒域。这里方位频率轴f_eta要经过fftshift使频率从负到正排列后面的相位因子才能和矩阵行对应上。变标相位因子写为$$H_1(\tau,f_\eta)\exp!\left(j\pi K_r C_s(f_\eta)\left(\tau-\tau_{ref}(f_\eta)\right)^2\right)$$其中$\tau_{ref}(f_\eta)2R_{ref}/(cD(f_\eta))$是参考斜距R_ref对应的徙动轨迹。这一步做完所有距离单元的徙动曲线都变成了参考距离的曲线形状因此同时为下一步的“一致距离徙动校正”创造了条件。实际实现里K_r应该用距离压缩前测得的实际调频率而不是发射参数表里的理论值。发射机增益曲线、功放非线性都会让理想chirp与真实chirp有偏差用实测值可以让变标后的残余相位更小。3.2 第二步在二维频域做距离压缩与一致RCMC变标后的信号做距离向FFT进入二维频域。此时距离chirp的调频率已经从K_r变成了$K_{eff}K_r(1C_s)$所以距离压缩匹配滤波器的调频率也要跟着变这是CSA中容易写错的地方。距离压缩项与一致距离徙动校正项合并写为$$H_2(f_\tau,f_\eta)\exp!\left(j\frac{\pi f_\tau^2}{K_r(1C_s)}\right)\cdot\exp!\left(-j\frac{4\pi R_{ref}}{c}\left(\frac{1}{D}-1\right)f_\tau\right)$$第一项是chirp匹配滤波第二个相位项完成参考距离处的一次性RCMC。注意R_ref要取场景中心斜距取偏了会让整个场景图像整体偏移但不影响聚焦取对了图像位置才是绝对定位正确的。这个域里还有一个经常被忽略的细节距离压缩的匹配滤波器带宽要和发射信号带宽匹配。采样率通常取带宽的1.2~1.4倍匹配滤波后在频带边缘会有幅度衰减如果后续要加窗压低旁瓣这部分衰减也要计入。3.3 第三步距离IFFT后做方位压缩与残余相位补偿二维频域处理完距离向IFFT回到距离-多普勒域。信号此时已经在正确的距离单元上剩下的方位向相位仍是一个关于f_eta的二次函数。方位压缩匹配滤波器为$$H_3(f_\eta,R_0)\exp!\left(j\frac{4\pi R_0}{\lambda}D(f_\eta)\right)$$其中R_0取当前距离单元对应的斜距$R_0c\tau/2$。H1引入的变标过程还会带来一个与τ无关的残余相位需要一并补偿$$\phi_{res}\pi K_r(1-D)\cdot\frac{4(R_0-R_{ref})^2}{c^2D^2}$$这一步只影响相位不改变幅度聚焦。做幅度成像时可以偷懒不乘但后续要做干涉测量或者极化定标残余相位必须补偿掉。下面的代码片段给出三个相位因子的构造方式频域轴与矩阵轴的对应关系以注释说明# 距离-多普勒域变标因子与H1 Cs 1.0 / D - 1.0 tau_ref 2.0 * R_ref / (c * D) H1 np.exp(1j * np.pi * Kr * Cs[:, None] * (tau[None, :] - tau_ref[:, None])**2) # 二维频域距离压缩与一致RCMC K_eff Kr * (1.0 Cs)[:, None] H2 np.exp(1j * np.pi * f_tau[None, :]**2 / K_eff) H2 * np.exp(-1j * 4.0 * np.pi * R_ref * (1.0 / D[:, None] - 1.0) * f_tau[None, :] / c) # 距离-多普勒域方位压缩与残余相位补偿 R_img c * tau / 2.0 H3 np.exp(1j * 4.0 * np.pi * R_img[None, :] * D[:, None] / lam) phi_res np.pi * Kr * (1.0 - D[:, None]) * ( 4.0 * (R_img[None, :] - R_ref)**2 / (c**2 * D[:, None]**2) ) H3 * np.exp(-1j * phi_res)代码里所有乘法和指数都是逐元素操作广播时距离频率、方位频率的轴关系必须和前面FFT后的矩阵保持一致。这三个相位因子分别作用在各自的域中顺序不能换H1在距离-多普勒域H2在二维频域H3在距离-多普勒域。后面第4章会把完整流程串起来。4. 用Python复现机载SAR的CS算法点目标仿真全流程4.1 机载SAR成像仿真参数怎么定点目标仿真不是为了模拟真实雷达回波的复杂电磁环境而是为了验证算法本身。参数取值要让距离、方位分辨率在同一量级并且保证方位向采样满足要求。下面这组参数接近小型无人机或通航飞机平台搭载的X波段SAR工作状态参数符号数值说明载频fc9.6 GHzX波段波长约3.1 cm距离带宽B180 MHz理论距离分辨率约0.83 m脉冲宽度Tp5 us线性调频脉宽距离采样率Fs220 MHz带宽的1.22倍留过采样平台速度V150 m/s典型小型平台速度平台高度H8000 m中低空机载场景中心斜距R_ref12000 m决定下视角约41.8度方位天线口径L_a2 m方位分辨率约1 m脉冲重复频率PRF320 Hz方位多普勒带宽约150 Hz余量2.1倍方位采样点数N_az1024覆盖完整合成孔径距离采样点数N_rg2048距离窗约1395 m多普勒带宽按$B_a2V/L_a$估算为150 HzPRF取320 Hz有足够余量但比理论最小值只高了1.1倍左右数据率仍在机载记录设备可承受范围内。快时间轴以场景中心双程延时为中心距离窗宽要覆盖目标之间的斜距差。4.2 回波生成代码与参数初始化import numpy as np c 299792458.0 fc 9.6e9 # 载频 lam c / fc B 180e6 # 距离向带宽 Tp 5e-6 # 脉冲宽度 Kr B / Tp # 调频率 Fs 220e6 # 距离向采样率 V 150.0 # 平台速度 m/s R_ref 12000.0 # 场景中心斜距 PRF 320.0 # 点目标列表(方位向位置, 最近斜距) targets [(0.0, R_ref), (-50.0, R_ref 60.0), (50.0, R_ref - 50.0)] # 慢时间轴以合成孔径中心为0 N_az 1024 eta (np.arange(N_az) - N_az // 2) / PRF # 快时间轴以场景中心双程延时为中心 N_rg 2048 tau np.arange(N_rg) / Fs tau tau - N_rg / (2.0 * Fs) 2.0 * R_ref / c # 逐脉冲生成基带回波 echo np.zeros((N_az, N_rg), dtypecomplex) for x_t, R0 in targets: R np.sqrt(R0**2 (V * eta - x_t)**2) # 每个慢时刻的斜距 tau_d 2.0 * R / c t_diff tau[None, :] - tau_d[:, None] in_pulse (t_diff 0) (t_diff Tp) phase (np.pi * Kr * t_diff**2 - 4.0 * np.pi * R[:, None] / lam) echo np.where(in_pulse, np.exp(1j * phase), 0.0) # 方位FFT进入距离-多普勒域 S_rd np.fft.fftshift(np.fft.fft(echo, axis0), axes0) f_eta np.fft.fftshift(np.fft.fftfreq(N_az, 1.0 / PRF))回波仿真直接在快时间采样点上离散计算chirp相位省去了插值步骤速度也快。注意相位里第一个指数项是基带chirp第二个指数项是载频对应的方位相位历程两者缺一不可。目标坐标里既有方位偏移又有斜距差是为了检查CS算法在二维方向上的聚焦能力而不是只测中心点。4.3 CS成像核心代码与参数说明# 距离-多普勒域的多普勒因子 D np.sqrt(1.0 - (lam * f_eta / (2.0 * V))**2) # 第一步CS变标 Cs 1.0 / D - 1.0 tau_ref 2.0 * R_ref / (c * D) H1 np.exp(1j * np.pi * Kr * Cs[:, None] * (tau[None, :] - tau_ref[:, None])**2) S_rd S_rd * H1 # 距离FFT进入二维频域 S_fd np.fft.fftshift(np.fft.fft(S_rd, axis1), axes1) f_tau np.fft.fftshift(np.fft.fftfreq(N_rg, 1.0 / Fs)) # 第二步距离压缩 一致RCMC K_eff Kr * (1.0 Cs)[:, None] H2 np.exp(1j * np.pi * f_tau[None, :]**2 / K_eff) H2 * np.exp(-1j * 4.0 * np.pi * R_ref * (1.0 / D[:, None] - 1.0) * f_tau[None, :] / c) S_fd S_fd * H2 # 距离IFFT回到距离-多普勒域 S_rd np.fft.ifft(np.fft.ifftshift(S_fd, axes1), axis1) # 第三步方位压缩 残余相位补偿 R_img c * tau / 2.0 H3 np.exp(1j * 4.0 * np.pi * R_img[None, :] * D[:, None] / lam) phi_res np.pi * Kr * (1.0 - D[:, None]) * ( 4.0 * (R_img[None, :] - R_ref)**2 / (c**2 * D[:, None]**2) ) H3 * np.exp(-1j * phi_res) S_rd S_rd * H3 # 方位IFFT输出复图像 img np.fft.ifft(np.fft.ifftshift(S_rd, axes0), axis0)这里做的是忽略二次距离压缩SRC的简化CS版本在正侧视、窄带条件下误差很小。决定是否要补上SRC项的判据是$\lambda^2 R_0 B_a^2/(16 c^2 V^2)$乘以适当系数后是否超过四分之一相位周期对上面的参数这顶多带来零点几弧度的相位误差肉眼分辨不出来但做干涉时要按完整CSA算一遍。代码里三个H因子的符号非常重要任何一个变号都会导致方位向散焦或者距离向偏移排查时先看距离压缩是不是把chirp压缩成单个尖峰。5. 机载PRF、斜视角与多普勒中心估计CSA的三个调参要点5.1 PRF与方位模糊机载SAR的采样率选择机载平台速度通常只有几十到两百米每秒方位多普勒带宽并不大但PRF不能只按带宽的奈奎斯特下限取。方位信号的多普勒带宽为$B_a2V/L_a$对应上面的参数是150 HzPRF取320 Hz看似没问题但要注意方位天线方向图的主瓣外还有旁瓣多普勒频率超过PRF/2时会折叠进主瓣区形成方位模糊。工程上常见做法是取PRF为多普勒带宽的1.2到1.5倍同时留出多普勒中心估计误差的余量。机载平台的GPS/惯导精度有限实际航线速度、偏流角都在缓慢变化多普勒中心不可能是零。如果PRF余量留得太小估计出的多普勒中心一旦超过PRF/2就会模糊到另一边图像方位向整体倒置这个问题比散焦更难查。在实时处理链里一般会先用惯性导航数据计算预期多普勒中心再用数据估计偏差两个值相差超过10%就说明数据或惯导有一方出了问题。5.2 多普勒中心不准会怎样先估计再成像CS算法推导时默认回波在方位谱上的中心是零频。实际机载SAR由于存在偏流角或天线安装角多普勒中心不为零直接送入CSA会让点目标在方位向偏移且天线方向图加权后的信噪比下降。先用回波数据估计多普勒中心是在成像前必做的一步。最简单的是能量重心法取一个较强距离单元做方位向FFT频谱重心就是多普勒中心的粗估计from numpy.fft import fft, fftshift, fftfreq row echo[500, :].copy() # 取某个距离单元 spec np.abs(fftshift(fft(row, axis0)))**2 f_eta_tmp fftshift(fftfreq(N_az, 1.0 / PRF)) fdc_est np.sum(f_eta_tmp * spec) / np.sum(spec) # 在方位时域去斜 eta2 (np.arange(N_az) - N_az // 2) / PRF echo echo * np.exp(-1j * 2.0 * np.pi * fdc_est * eta2[:, None])估计完多普勒中心后回波在方位向乘以线性相位exp(-j2πfd_est·η)把频谱中心平移回零频。之后再做CS成像图像里各目标方位位置才是真实位置。注意能量重心法对强点目标效果好一片均匀杂波时反而会偏更稳的是用相邻方位频谱的互相关法但粗估阶段重心法够用。5.3 大斜视角下CS变标失效的边界CS变标因子$C_s1/D-1$是在正侧视几何下推导的。当斜视角增大距离-多普勒域的精确表达式里会多出斜视角相关的耦合项简化CSA的点目标响应会从两侧旁瓣不对称逐渐演化为主瓣展宽。经验边界是10到15度以内用简化CSA问题不大超过这个范围要使用包含二次距离压缩的完整CSA或者直接换ωK算法。机载平台的斜视角通常不是定值。飞机受侧风影响实际航迹和天线指向的夹角在一条航线上可能波动好几度。惯导数据更新的速率如果低于PRF这种波动会直接变成方位向相位误差CSA这类频域算法对相位误差的敏感度又比RD高。我一般会在处理链里把回波按航迹段分块每块单独估算等效速度V和斜视角再用分块后的参数做CSA成像而不是全航线用一组参数。6. 点目标响应的定量评估用IRW、PSLR、ISLR看CS成像质量6.1 为什么评估点目标响应要插值到更高采样率成像结果栅格通常只按原采样率排列主瓣峰值两边只有几个采样点直接数主瓣宽度会引入很大的量化误差。评估点目标响应的常规做法是把峰值附近的距离剖面和方位剖面各取一段用零填充IFFT插值到8倍采样率再做3dB宽度和旁瓣统计。插值后主瓣不会被栅格“切歪”PSLR和ISLR也才稳定。6.2 用Python计算IRW、PSLR和ISLRdef measure_metrics(profile, dr, zoom8): # profile: 一维复数剖面dr为采样间隔 N len(profile) p np.abs(np.fft.irfft(np.fft.rfft(profile), N * zoom))**2 peak np.argmax(p) power p / p[peak] # 3dB主瓣宽度用峰值左右第一个低于0.5的点做线性插值 left np.where(power[:peak] 0.5)[0][-1] right np.where(power[peak:] 0.5)[0][0] peak irw (right - left) / zoom * dr # 主瓣边界取第一个零点 zleft np.where(power[:peak] 1e-3)[0][-1] zright np.where(power[peak:] 1e-3)[0][0] peak sidelobe np.concatenate([power[:zleft], power[zright:]]) pslr 10 * np.log10(np.max(sidelobe)) islr 10 * np.log10( (np.sum(power) - np.sum(power[zleft:zright])) / np.sum(power[zleft:zright]) ) return irw, pslr, islr主瓣边界用第一个零点近似划分实际处理中也可以固定取峰值两侧各若干个3dB宽度作为主瓣范围。PSLR反映强目标旁边的弱目标可检测能力ISLR反映分布式目标背景的泄漏水平。对未加窗的矩形谱响应PSLR接近-13.3 dBISLR接近-10 dB加Hamming窗后PSLR会压到-42 dB左右但IRW会扩展到约1.3倍。6.3 验证CS算法实现正确性的三个信号特征检查距离剖面和方位剖面的IRW是否接近理论值。距离向理论IRW为$0.886\cdot c/(2B)$方位向理论IRW为$0.886\cdot \lambda R_0/(2V T_{obs})$差距超过5%就说明某个相位因子有误差。第二三个点目标的峰值应出现在各自设定的斜距和方位坐标上距离向整体偏移说明R_ref取偏方位向出现随距离变化的位移说明多普勒中心估计遗漏。第三残余相位补偿H3被关掉时幅度图像看不出明显变化但点目标相位沿方位向会出现抛物线状残余用干涉或相位图检查能立刻暴露问题。这三个特征都正常CS成像链基本可以放心交给后面的自动聚焦和辐射定标模块处理。本文还有配套的精品资源点击获取