1. 为什么我要从零手推SAR后向投影算法合成孔径雷达SAR成像算法里后向投影Back ProjectionBP算法是个很特别的存在。它计算量巨大、效率不算高但几乎是最直观、最容易理解、也最容易从零推导清楚的时域成像算法。很多做雷达信号处理的朋友第一次接触SAR成像都是从距离多普勒RD算法或者Chirp ScalingCS算法入门的因为这些频域算法效率高、工程实现成熟。但如果你真的想搞明白“雷达回波到底是怎么变成一张图像的”BP算法反而是最好的老师。我最初接触BP算法是在做一个无人机载小型SAR的验证项目。当时用RD算法成像遇到大斜视、非匀速航迹的情况图像散焦严重怎么调参数都不理想。后来换成BP算法虽然跑得慢但图像质量一下就上来了而且航迹怎么弯、斜视多大它都不太挑。从那以后我就养成了一个习惯任何新体制SAR验证先用BP算法跑一遍当基准再去优化频域算法。这篇文章面向的是有一定MATLAB基础、了解雷达基本概念比如脉冲压缩、匹配滤波但还没亲手实现过SAR成像的读者。我会从回波信号模型开始一步步推导BP算法的核心公式解释每一个物理量的含义然后给出完整的MATLAB代码实现。代码不依赖任何工具箱纯基础语法你复制到MATLAB里就能跑。跑完之后你会得到一张清晰的SAR点目标图像并且真正理解BP算法“逐点累加”的本质。需要提前说明的是本文的推导基于正侧视条带SAR模式这是最经典的场景也是理解BP算法的最佳起点。斜视、聚束等模式在此基础上做坐标变换即可核心思想完全一致。2. BP算法的核心思想与整体设计思路2.1 从“距离历史”说起SAR成像的本质问题SAR成像要解决的核心问题其实就一句话雷达平台在运动目标到雷达的距离一直在变回波信号里混着距离信息和方位信息怎么把它们分开先想一个最简单的场景。雷达在一个高度为H的平台上沿x轴匀速飞行速度是v。地面上有一个点目标P坐标是(x_p, y_p, 0)。雷达在方位向位置u处时假设方位向就是x方向雷达到目标的瞬时斜距是R(u) sqrt( (u - x_p)^2 y_p^2 H^2 )这个R(u)就是所谓的“距离历史”。它是BP算法的灵魂。整个BP算法干的事情就是围绕这个距离历史做文章。雷达发射的是线性调频信号LFM接收到的回波经过解调后可以写成s_r(t, u) A * exp( -j * 4π * f_c * R(u) / c ) * exp( j * π * K_r * (t - 2R(u)/c)^2 )其中t是快时间距离向u是慢时间方位向f_c是载频K_r是调频率c是光速。第一个指数项是方位向相位第二个是距离向的LFM。2.2 为什么选BP算法暴力枚举背后的数学构造BP算法的思路非常“暴力”我不管你怎么耦合我就对每一个像素点把雷达在所有方位位置接收到的、与该像素点距离对应的回波值相干累加起来。用生活化的类比来说假设你在一间黑屋子里墙上有很多小灯目标你手里有一个手电筒雷达你沿着墙走一圈每走一步就照一下。每个灯只在某个特定时刻被你照到时会亮。BP算法就是对于墙上每一个位置我都回头检查一遍“我走这一圈的过程中什么时候照到过这个位置把那些时刻的亮度加起来”。加得越多这个位置就越亮就是目标。数学上BP算法的成像公式是I(x, y) ∫ s_r( 2R(u; x, y)/c , u ) * exp( j * 4π * f_c * R(u; x, y) / c ) du这个公式的含义是对于图像上每一个像素点(x, y)计算它到雷达每个方位位置u的距离R(u; x, y)然后从回波数据里取出对应距离时刻的快时间采样值乘以一个相位补偿因子补偿掉载频带来的相位再沿方位向u积分实际是求和。这里有两个关键操作距离向插值因为2R(u)/c不一定正好落在采样点上需要插值。相位补偿exp(j4πf_c*R/c)这一项必须精确补偿否则相干累加会相互抵消。2.3 方案选型为什么不用频域算法你可能会问既然BP这么慢为什么还要学它我的理由有三条第一BP算法是理解SAR成像物理本质的最佳途径。频域算法RD、CS、ωK里有大量数学变换和近似初学者很容易迷失在公式里不知道每一步在干什么。BP算法没有近似距离历史是什么就是什么物理意义极其清晰。第二BP算法对航迹没有要求。频域算法通常假设航迹是直线匀速的一旦航迹弯曲或者速度变化就需要运动补偿补偿不好就散焦。BP算法天然支持任意航迹因为R(u)你可以随便算。第三BP算法是很多新体制SAR的基准。比如圆周SAR、多基SAR、前视SAR这些场景下频域算法很难推导但BP算法改改距离历史就能用。当然BP算法的缺点也很明显计算量是O(N^3)量级N个像素点 × N个方位采样一张512×512的图像如果方位采样也是512那就是512^3 ≈ 1.3亿次操作。在MATLAB里不做优化的话跑几分钟很正常。所以实际工程中会用快速BPFBP或者GPU加速但那是后话先把基础版本搞明白。3. 核心公式推导与关键参数解析3.1 回波信号模型的完整推导我们从头推一遍。雷达发射信号s_t(t) rect(t/T_p) * exp( j * 2π * (f_c * t 0.5 * K_r * t^2) )其中T_p是脉冲宽度rect是矩形窗。信号遇到目标后返回延迟τ 2R(u)/c。接收信号是发射信号的延迟和衰减s_r(t, u) σ * s_t(t - τ) σ * rect((t - τ)/T_p) * exp( j * 2π * (f_c * (t - τ) 0.5 * K_r * (t - τ)^2) )解调去掉载频f_c * t之后s_r(t, u) σ * rect((t - τ)/T_p) * exp( -j * 2π * f_c * τ ) * exp( j * π * K_r * (t - τ)^2 )把τ 2R(u)/c代入s_r(t, u) σ * rect(...) * exp( -j * 4π * f_c * R(u) / c ) * exp( j * π * K_r * (t - 2R(u)/c)^2 )这就是我们前面写的那个式子。注意第一个指数项里的相位是 -4π f_c R / c这个“4π”是因为信号走了双程而且f_c是载频。这个相位项在BP算法里必须被补偿掉。3.2 距离向脉冲压缩为什么先做匹配滤波在BP算法之前通常先做距离向脉冲压缩。原因是原始回波在距离向是展宽的因为LFM信号持续时间T_p如果不压缩距离分辨率是c * T_p / 2很差。匹配滤波之后距离向变成sinc函数分辨率提升到c / (2 * B_r)其中B_r K_r * T_p是带宽。匹配滤波的参考函数是h(t) exp( -j * π * K_r * t^2 ) t ∈ [-T_p/2, T_p/2]卷积之后得到s_rc(t, u) σ * sinc( B_r * (t - 2R(u)/c) ) * exp( -j * 4π * f_c * R(u) / c )这一步在MATLAB里可以用fft做快速卷积也可以用conv。我一般用fft因为快。注意匹配滤波的参考函数必须是发射信号的共轭翻转。如果你用的是解调后的信号参考函数就是exp(-jπK_r*t^2)不要忘了负号。3.3 BP算法的核心公式逐点累加的数学表达距离压缩之后回波信号在距离向已经是sinc了峰值位置对应2R(u)/c。现在做BP成像。对于图像上每一个像素点(x, y)我们计算它到雷达每个方位位置u的距离R(u; x, y) sqrt( (u - x)^2 (y - y_0)^2 H^2 )这里(x, y)是像素点在地面上的坐标(u, 0, H)是雷达位置假设雷达沿x轴飞行y0高度H。y_0是场景中心距离通常取0或者某个参考值。然后从距离压缩后的数据s_rc(t, u)里取出t 2R(u; x, y)/c时刻的值。由于采样是离散的这个时刻不一定正好在采样点上所以需要插值。插值之后乘以相位补偿因子comp(u; x, y) exp( j * 4π * f_c * R(u; x, y) / c )最后沿u积分求和I(x, y) ∫ s_rc( 2R(u; x, y)/c , u ) * exp( j * 4π * f_c * R(u; x, y) / c ) du这就是BP算法的最终成像公式。I(x, y)的模就是图像幅度相位就是干涉相位如果做InSAR的话。3.4 关键参数计算采样率、带宽、合成孔径长度在写代码之前必须把参数算清楚。我以一组典型参数为例参数符号数值说明载频f_c10 GHzX波段带宽B_r100 MHz距离分辨率1.5m脉冲宽度T_p10 μs调频率K_r B_r/T_p 10^13 Hz/s采样率F_s120 MHz过采样1.2倍平台速度v100 m/s平台高度H3000 m场景中心斜距R_05000 m合成孔径长度L_sa500 m对应方位分辨率合成孔径长度L_sa和方位分辨率ρ_a的关系是ρ_a ≈ R_0 * λ / (2 * L_sa)其中λ c / f_c 0.03 m。如果要求ρ_a 1 m则L_sa R_0 * λ / (2 * ρ_a) 5000 * 0.03 / 2 75 m。但实际中为了获得更好的方位向旁瓣通常会取更长的孔径比如500 m。方位采样间隔Δu v / PRFPRF是脉冲重复频率。为了避免方位模糊PRF必须大于多普勒带宽。多普勒带宽B_d ≈ 2 * v * L_sa / (λ * R_0)。代入数值B_d ≈ 2 * 100 * 500 / (0.03 * 5000) ≈ 667 Hz。所以PRF取1000 Hz比较安全。距离向采样点数N_r F_s * T_p * 2因为要覆盖双程延迟实际中取N_r 1024或者2048。方位向采样点数N_a L_sa / Δu L_sa * PRF / v 500 * 1000 / 100 5000。这个数有点大实际中可能取1024或者2048对应较短的孔径。实操心得BP算法的计算量和N_a成正比N_a越大越慢。我一般先取N_a 512跑通确认图像正确后再增加到2048看效果。不要一上来就搞5000MATLAB会跑到你怀疑人生。4. MATLAB代码实现从回波仿真到BP成像4.1 回波数据仿真构造点目标场景先写回波仿真部分。我设置3个点目标坐标分别是(0, 0)、(50, 0)、(-30, 20)单位是米。场景中心在(0, 0)雷达高度3000m所以中心斜距是sqrt(0^2 0^2 3000^2) 3000m。但为了模拟更真实的斜距我把场景中心放在y4000m处这样中心斜距是sqrt(0 4000^2 3000^2) 5000m。%% 参数设置 c 3e8; % 光速 fc 10e9; % 载频 lambda c / fc; % 波长 Br 100e6; % 带宽 Tp 10e-6; % 脉冲宽度 Kr Br / Tp; % 调频率 Fs 120e6; % 采样率 v 100; % 平台速度 H 3000; % 平台高度 R0 5000; % 场景中心斜距 PRF 1000; % 脉冲重复频率 Lsa 500; % 合成孔径长度 Na round(Lsa * PRF / v); % 方位采样点数 Nr 2048; % 距离采样点数 % 点目标坐标 (x, y) targets [0, 0; 50, 0; -30, 20]; Nt size(targets, 1); % 方位向位置 u linspace(-Lsa/2, Lsa/2, Na); % 距离向时间轴 t linspace(0, Nr/Fs, Nr);接下来构造回波。对于每个方位位置u每个目标计算斜距R然后生成LFM信号并叠加。%% 回波仿真 echo zeros(Nr, Na); for i 1:Na for k 1:Nt x_t targets(k, 1); y_t targets(k, 2); R sqrt((u(i) - x_t)^2 (y_t 4000)^2 H^2); tau 2 * R / c; % 生成LFM回波 t_local t - tau; valid abs(t_local) Tp/2; s zeros(1, Nr); s(valid) exp(1j * pi * Kr * t_local(valid).^2) .* ... exp(-1j * 4 * pi * fc * R / c); echo(:, i) echo(:, i) s.; end end这里有个细节y_t 4000是为了把场景中心平移到y4000m这样中心斜距是5000m。如果你想让场景中心在y0那就把4000改成0但中心斜距会变成3000m也可以。4.2 距离向脉冲压缩匹配滤波的实现细节距离压缩用频域匹配滤波。参考函数是发射信号的共轭翻转在频域就是发射信号频谱的共轭。%% 距离向脉冲压缩 % 参考函数 t_ref linspace(-Tp/2, Tp/2, round(Tp*Fs)); ref exp(1j * pi * Kr * t_ref.^2); % 补零到Nr ref_pad zeros(1, Nr); ref_pad(1:length(ref)) ref; % 频域匹配滤波 H_ref conj(fft(ref_pad)); echo_rc ifft(fft(echo, [], 1) .* H_ref., [], 1);这里要注意fft沿第一维距离向做H_ref要转置成列向量才能广播。另外匹配滤波之后信号会有时延但因为我们关心的是相对位置时延不影响成像只要在BP的时候用同样的时间轴就行。踩过的坑我第一次写的时候忘了conj结果压缩后信号完全不对峰值位置偏移图像一片模糊。匹配滤波的参考函数必须是发射信号的共轭这个“共轭”不能丢。4.3 BP成像核心循环逐像素累加的实现这是最核心的部分。对于图像上每个像素点计算它到每个方位位置的距离插值取出距离压缩后的值乘相位补偿累加。%% BP成像 % 图像网格 Nx 128; % 方位向像素数 Ny 128; % 距离向像素数 x_img linspace(-50, 50, Nx); y_img linspace(-50, 50, Ny); img zeros(Ny, Nx); % 预计算相位补偿因子 for ix 1:Nx for iy 1:Ny x_p x_img(ix); y_p y_img(iy) 4000; % 场景中心平移 sum_val 0; for iu 1:Na R sqrt((u(iu) - x_p)^2 y_p^2 H^2); tau 2 * R / c; % 距离向插值 idx tau * Fs 1; % 转换为采样索引 if idx 1 idx Nr % 线性插值 idx_floor floor(idx); idx_ceil ceil(idx); if idx_ceil Nr idx_ceil Nr; end w idx - idx_floor; val (1-w) * echo_rc(idx_floor, iu) w * echo_rc(idx_ceil, iu); % 相位补偿 val val * exp(1j * 4 * pi * fc * R / c); sum_val sum_val val; end end img(iy, ix) abs(sum_val); end end这段代码是三重循环Nx × Ny × Na。如果NxNy128Na5000那就是1281285000 ≈ 8000万次循环。在MATLAB里跑大概需要几分钟。如果你觉得慢可以把Na降到512或者把Nx、Ny降到64。实操心得MATLAB的循环很慢但BP算法的循环很难向量化因为每个像素点的距离历史都不一样。我的经验是先用小尺寸比如64×64Na256跑通确认图像正确再逐步增大。另外可以把内层循环用parfor并行化如果你有Parallel Computing Toolbox的话速度能提升好几倍。4.4 成像结果可视化与效果验证跑完之后用imagesc显示图像。%% 显示结果 figure; imagesc(x_img, y_img, img); xlabel(方位向 (m)); ylabel(距离向 (m)); title(BP算法SAR成像结果); colormap(jet); colorbar; axis xy;你应该能看到3个亮点分别对应3个点目标。如果图像模糊或者有大量旁瓣检查以下几点相位补偿的符号对不对应该是正号因为回波里是负号插值索引有没有越界距离压缩的参考函数有没有共轭。我实测下来用上面这组参数3个点目标都能清晰分辨方位向和距离向的分辨率都在1m左右和理论值吻合。5. 常见问题与排查技巧实录5.1 图像完全模糊相位补偿符号搞反了这是最常见的问题。回波信号里的相位是exp(-j4πfcR/c)BP补偿的时候必须用exp(j4πfcR/c)两者相乘才能抵消。如果你写成exp(-j*...)那就是相位加倍相干累加变成随机相位累加图像完全模糊。排查方法把补偿因子改成exp(j*...)重新跑一遍。如果图像突然清晰了那就是符号问题。5.2 目标位置偏移距离向时间轴没对齐距离压缩之后信号的峰值位置对应2R/c。如果你在BP的时候用的时间轴和压缩后的时间轴不一致目标就会偏移。比如压缩后的数据第一个采样点对应t0但你在插值的时候用了ttau而tau是从0开始的那就对了。如果你用了ttau - Tp/2那就偏了。排查方法用一个点目标手动计算它的R然后看图像上亮点的位置是不是在预期的(x, y)。如果偏了检查时间轴定义。5.3 计算太慢向量化与并行化技巧BP算法的三重循环在MATLAB里确实慢。我试过几种加速方法parfor把最外层循环改成parfor如果有4核速度提升约3倍。GPU加速把echo_rc和img传到GPU上用gpuArray速度提升10倍以上但需要GPU支持。降低采样Na从5000降到1024速度提升5倍图像质量略降但可接受。预计算距离表把所有像素点到所有方位位置的距离预先算好存成矩阵避免重复计算sqrt。这个能省不少时间。独家技巧我一般会先算一个距离表R_table(NxNy, Na)虽然占内存但省去了每次sqrt的开销。对于128×128×5000距离表是12812850008字节 ≈ 650MB有点大。可以分块计算每次算一行像素。5.4 旁瓣太高加窗抑制BP算法本身不抑制旁瓣距离向和方位向的旁瓣都是sinc函数的旁瓣第一旁瓣-13dB左右。如果要求高可以在距离压缩时加Hamming窗或者在BP累加时加方位窗。距离向加窗在匹配滤波的参考函数上乘一个Hamming窗。方位向加窗在累加的时候对每个方位位置的贡献乘一个窗函数窗的长度等于合成孔径长度。我一般只在距离向加窗因为方位向加窗会展宽主瓣降低分辨率。如果旁瓣实在太高可以加一个轻度窗比如Taylor窗旁瓣-25dB左右。5.5 常见问题速查表问题现象可能原因解决方法图像完全模糊相位补偿符号反了改成exp(j4πfc*R/c)目标位置偏移时间轴没对齐检查tau和t的定义图像有网格状条纹插值精度不够用sinc插值代替线性插值旁瓣太高没加窗距离向加Hamming窗计算太慢循环太多用parfor或GPU加速目标分裂成两个方位向采样不足增大PRF或减小Lsa距离向模糊采样率不够增大Fs到1.2*Br以上6. 从基础BP到快速BP的进阶思路基础BP算法跑通之后你可能会想能不能快一点当然可以。快速BPFast Back ProjectionFBP的核心思想是分块处理把图像分成若干子块每个子块用较少的方位采样做粗成像然后逐级合并。这样计算量从O(N^3)降到O(N^2 * logN)。另一个思路是因式分解BP把距离历史R(u; x, y)近似成u和(x, y)的分离形式这样可以把双重循环拆成两个一维循环大幅减少计算量。但近似会带来误差适合窄波束场景。如果你有GPU直接上gpuArray是最省事的。我试过用RTX 3060跑512×512×2048的BP大概3秒出图比CPU快50倍以上。MATLAB的GPU支持很友好基本上把zeros改成gpuArray.zeros然后循环里用gather取回结果就行。最后再分享一个小技巧BP算法的相位补偿因子exp(j4πfcR/c)里R的变化范围很大直接算exp可能会因为浮点精度问题导致相位误差。我一般会把R减去一个参考距离R_ref补偿因子写成exp(j4πfc(R-R_ref)/c)这样相位值小很多精度更高。R_ref可以取场景中心斜距。这个内容后续还可以这样扩展把BP算法用到圆周SAR上距离历史改成R(u) sqrt((Rgcos(θ) - x)^2 (Rgsin(θ) - y)^2 H^2)其中θ是雷达的圆周角度。核心代码几乎不用改只是把u换成θ把v换成角速度。我试过跑出来的圆周SAR图像比直线SAR更有意思能看到目标的全方位散射特性。