简介面向合成孔径雷达SAR成像研究与工程应用的MATLAB实现专注于聚束式SAR中的极坐标格式算法PFA适用于遥感测绘、地形侦察、目标识别等成像处理场景也适合雷达信号处理方向的高校学生作为课程设计或毕业设计参考。压缩包内仅包含一个PFA.m脚本文件包体约2KB代码量不大却覆盖PFA核心处理链路需要读者具备一定的雷达原理与数字信号处理基础。目前已有1635人学习或下载具有一定的实践参考热度。脚本围绕数据处理与图像重建展开涉及原始回波校准去噪、距离压缩、多普勒参数估计、方位向FFT聚焦、极坐标转换以及最终图像重构等关键步骤可以直观展示从回波到SAR图像的主要处理流程同时方便结合MATLAB信号处理与图像处理工具箱进行二次开发和算法调试为深入理解聚束式SAR成像机制提供简洁的代码级参考。1. PFA 到底在解决 SAR 里的哪一个矛盾传统 SAR 成像最常用的 Range-Doppler 算法沿距离向做脉冲压缩、沿方位向做匹配滤波但它在聚束 SAR 或者大斜视角条带 SAR 下会碰到一个根本矛盾目标回波的方位向相位历史不是时间的线性函数而是随目标与平台之间的斜距变化呈双曲线甚至更复杂的形态。若直接在原始回波域做二维傅里叶变换散焦会大到主瓣完全分裂根本无法形成图像。PFA 的思路是把回波数据从极坐标采样网格转换到直角坐标网格让变换后的数据满足二维傅里叶变换的条件。它并不是对原始回波的直接压缩而是先做几何重排再聚焦这也是 PFA 区别于 Range-Doppler 和 Chirp Scaling 的根本特征。适合做机载或星载聚束 SAR 的成像处理也适用于条带 SAR 的子孔径成像对工程中已经做过去调频、采集到基带信号的雷达数据尤其合适。2. 距离向处理Dechirp、快时间 FFT 与斜距标定2.1 从线性调频回波到差频相位先去调频再谈成像SAR 发射的线性调频信号可以写成 $s_t(\tau) \exp[j2\pi(f_c\tau \frac{1}{2}K_r\tau^2)]$其中 $\tau$ 是快时间$K_r$ 是调频率$f_c$ 是载频。接收回波相对发射信号有一段延时 $\tau_d 2R/c$这里的 $R$ 是目标在某个方位时刻的斜距。如果不做任何处理直接采样保存那每个脉冲都是完整的一段 chirp数据率极高而且后续距离脉冲压缩需要对全带宽做匹配滤波计算量非常大。所以实际系统在射频端就用本地参考 chirp 和接收信号做混频这就是 Dechirp也叫去调频、去斜、Stretch 处理。去调频以后距离向信号变成了单频信号。目标距离越远差频频率越高。这个差频信号再做一次 FFT得到的频率轴就对应距离轴不需要再做传统意义上的匹配滤波一次 FFT 就完成了距离压缩。这里的关键是差频相位中还包含一个二次相位项称为残余视频相位RVP它来源于发射 chirp 的相位项与参考 chirp 相位项没有完全对消的部分。RVP 在距离向 FFT 后是一个与距离频率相关的线性相位如果不补偿点目标在距离向上会出现位置偏移方位向聚焦也会受影响。工程上常用的做法是用一个称为 Deskew 的过程消除 RVP把 Dechirp 数据换算成理想脉压后的数据。2.2 快时间维 FFT 和距离轴标定PFA 输入数据的第一步假设去调频后的基带信号为 $s_{if}(f_\tau, t)$其中 $f_\tau$ 是快时间频率$t$ 是慢时间方位维。对每一列慢时间数据做 FFT得到以快时间频率表示的距离压缩数据import numpy as np from scipy.fft import fft, fftshift def range_compress_dechirp(data_if, n_fft_rNone): data_if: 2D array, shape (n_azimuth, n_fast_time) 每行是一个脉冲去调频后的 I/Q 数据 n_fft_r : 距离向 FFT 点数, 默认等于快时间采样点数 n_az, n_fast data_if.shape if n_fft_r is None: n_fft_r n_fast # 沿第二维做 FFT, 对应距离向脉冲压缩 range_comp fftshift(fft(data_if, nn_fft_r, axis1), axes1) return range_comp这段代码的核心在于fft(data_if, nn_fft_r, axis1)指定了沿快时间维做 FFTfftshift把零频挪到中间便于后续按距离频率索引。注意 Dechirp 数据的 FFT 不需要做加窗因为距离维的窗函数效应已经由发射信号带宽决定了但实际系统中为了让旁瓣可控往往会在此处额外乘一个距离窗。这里最容易被忽略的是距离轴的标定快时间频率 $f_\tau$ 到斜距 $R$ 的换算关系是 $f_\tau -\frac{2K_r}{c}(R - R_{ref})$其中 $R_{ref}$ 是 Dechirp 所用的参考斜距。做极坐标格式转换时需要根据这个关系把频率轴映射到空间频率域不能直接用 FFT 的 bin 号当距离。2.3 Dechirp 参数选择参考距离、调频率与采样点数的约束Dechirp 的参考信号通常取场景中心处的回波延时参考距离一旦选偏距离维的中心频率就会偏离零频导致距离压缩后的数据频谱中心不在带内。参考距离偏差 $dR$ 引起的差频中心偏移是 $\Delta f -2K_r dR / c$所以 Dechirp 的参考斜距误差必须控制在一个距离单元内。调频率 $K_r$ 的标定更是重要如果系统给出的 $K_r$ 与实际有偏差Dechirp 后信号残余二次相位会增大直接表现为距离向主瓣展宽。PFA 对距离向调频率误差比 Range-Doppler 更敏感因为在极坐标重采样时残余二次相位会随方位角变化被映射到二维相位误差中形成空变散焦。实际系统中调频率的标定误差通常控制在 0.1% 以内。验证办法是检查距离压缩后点目标响应的峰值相位跨方位向的平坦度。理论上一个理想点目标在 Dechirp 后距离压缩输出的相位沿方位向只包含由斜距变化引起的多普勒相位不应该存在快时间频率维的二次项。做 PFA 之前最好先用点目标仿真数据把 Dechirp 链路调干净否则后面极坐标重采样和二维聚焦的问题会被误判为插值精度或运动补偿的问题。3. 极坐标格式转换距离向重采样与方位向插值的配合3.1 为什么原始数据是极坐标而不是直角坐标聚束 SAR 在成像期间天线始终指向场景中心所以每个方位时刻 $t$ 对应的波数矢量 $\mathbf{k}$ 的方向都在变化。目标到雷达的斜距 $R(t)$ 随平台位置变化回波在空间频率域的位置可以写成 $k_x \frac{4\pi}{\lambda}\cos\theta(t)$$k_y \frac{4\pi}{\lambda}\sin\theta(t)$ 的形式其中 $\theta(t)$ 是方位角。波长 $\lambda$ 对每个距离频率单元不同所以不同距离频率的数据对应不同的 $k_x, k_y$ 值。将全部脉冲的回波数据画在 $(k_x, k_y)$ 平面上会呈现扇形极坐标网格而不是标准的均匀直角栅格。二维 IFFT 要求数据在直角网格上均匀分布因此 PFA 的核心工作就是把扇形网格重建成矩形网格这包括沿距离向的频率重采样消除距离弯曲导致的距离单元徙动和沿方位向的插值。3.2 距离向重排将极坐标网格映射到直角坐标网格的编程实现PFA 常见的实现方式是先做距离向插值再做方位向插值。距离向插值的对象是快时间频率轴。假设成像场景中心位于 $R_c$目标点相对于场景中心的斜距为 $r_x$。通过在极坐标域用场景中心作为参考对距离向频率轴做坐标变换把每个方位角度下的频谱重排到直角坐标的某一个高度这样目标在距离向上的位置与方位角度不再耦合也就是把极坐标展开“拉直”。这里给出 PFA 中最关键的距离向重采样代码。数据已经完成距离向 FFT每一行对应一个方位脉冲行内是距离频率轴上的复数值。def polar_to_rect_resample(data_rg, freq_r, freq_x_new): 距离向极坐标到直角坐标的重采样。 data_rg : 2D array, shape (n_az, n_range), 距离压缩后数据 freq_r : 原始距离频率轴, 1D array, 长度 n_range freq_x_new: 新的直角坐标频率轴, 1D array, 长度 n_range_x n_az, n_range data_rg.shape n_rx freq_x_new.shape[0] out np.zeros((n_az, n_rx), dtypecomplex) for i in range(n_az): row data_rg[i, :] # 实际工程中可按方位角度实时计算映射关系 # 这里用线性插值做示意, 精确实现应使用 sinc 插值 out[i, :] np.interp(freq_x_new, freq_r, row.real) \ 1j * np.interp(freq_x_new, freq_r, row.imag) return out线性插值只适合快速验证流程真正成像时插值核长度至少要 8 点以上。原因是距离向重采样会引入插值误差这个误差表现为相位误差而 SAR 成像对相位误差极为敏感。通常是构造一个 8 点或 16 点的 sinc 核并加上 Kaiser 窗抑制截断旁瓣。freq_x_new的生成不是等间距就完事需要按照成像中心点的波数坐标来定义$k_x \frac{4\pi}{c}(f_c f_\tau)\cos\theta$每个方位角下的距离频率采样点都映射到直角坐标的 $k_x$ 位置新网格的起始和终止要覆盖全部数据的 $k_x$ 范围。3.3 方位向插值把极坐标的非均匀方位采样变成均匀栅格距离向重采样完成后数据在距离维上已经近似为直角网格但方位维上的频谱位置仍然随方位角变化数据点不在均匀的方位频率栅格上。需要把每一列对应一个距离频率的数据从极坐标方位角映射到直角坐标方位频率 $k_y$。这一步的本质是沿方位维做一次插值把每个距离频率单元在方位向上“搬到”正确的位置。方位向插值的方式有两种。第一种是先做方位向 FFT 到多普勒域再在频域做插值第二种是直接在时域慢时间域插值然后做 FFT。工程中常用后者因为 PFA 的方位向处理必须精确控制每个距离单元的相位时域插值更容易和运动补偿结合。插值的过程可以用下面的代码表达from scipy.interpolate import interp1d def azimuth_reposition(data_rect, k_y_old, k_y_new): 方位向重采样。 data_rect : 距离向重排后的数据, shape (n_az, n_rx) k_y_old : 每个方位时间对应的原方位频率 k_y_new : 均匀的直角坐标系方位频率 n_az, n_rx data_rect.shape out np.zeros((n_az, n_rx), dtypecomplex) for j in range(n_rx): col data_rect[:, j] real_interp interp1d(k_y_old, col.real, kindlinear, bounds_errorFalse, fill_value0) imag_interp interp1d(k_y_old, col.imag, kindlinear, bounds_errorFalse, fill_value0) out[:, j] real_interp(k_y_new) 1j * imag_interp(k_y_new) return outk_y_old在聚束 SAR 中的计算是 $k_y -k_x \tan\theta$也就是把每个方位时间点的斜视角映射到波数域的角度。k_y_new是等间距网格其间距要满足奈奎斯特条件通常是 $\Delta k_y \le \pi / X_{scene}$这里的 $X_{scene}$ 是方位向成像场景尺寸。如果网格间距选得太大图像方位向会混叠选得太小计算量增大但分辨率不会提升因为分辨率由合成孔径长度决定与网格密度无关。3.4 两次插值的顺序不能换讲清楚为什么距离向重采样必须在方位向插值之前完成。如果先做方位向插值数据在距离维仍然是极坐标分布方位向插值后每个距离单元内部的距离频率仍然与方位角耦合后续方位向 IFFT 会因为距离单元徙动没有被校正而散焦。距离向重采样本质上完成了距离走动与距离弯曲的校正它把弯曲的轨迹“拉直”成直线这样方位向处理才是真正的一维信号处理。但也有一种特殊情况如果场景尺寸很小距离弯曲量不足一个距离分辨单元可以省去距离向重采样只做方位向插值。这种简化 PFA 在小场景成像中很常见。判断准则很简单计算场景边缘目标的最大距离弯曲量 $\Delta R_{max} \frac{L_s^2}{8R_c}$其中 $L_s$ 是场景方位向尺寸$R_c$ 是场景中心斜距。若 $\Delta R_{max}$ 大于四分之一距离分辨率就要做距离向重采样。实际工程中这个判断非常重要能省掉一半计算量。4. 二维聚焦实现RVP 补偿、频域滤波与成像参数判定4.1 二维 IFFT 前的相位校正残余视频相位与自动聚焦入口极坐标重采样完成后数据已经位于均匀的直角网格上理论上做二维 IFFT 就能得到聚焦图像。但工程中还有两个环节不能跳过。第一个是 RVP 补偿。Dechirp 处理残留的 RVP 项在距离压缩后体现为一个关于快时间频率的线性相位这个相位在极坐标重采样后会映射到二维波数域成为沿 $k_x$ 方向的相位斜坡。补偿方法是把距离压缩后的数据乘以一个相位因子 $\exp(-j\pi f_\tau^2 / K_r)$。注意这个因子需要在极坐标重采样之前乘上否则无法和各方位角完全对齐。第二个是自动聚焦Autofocus的入口。PFA 对运动误差的敏感度非常高即使系统标定完美平台运动误差和大气扰动仍然会在波数域留下二维相位误差。因为极坐标格式本身把数据转换成了二维频域信号所以非常适合做相位梯度自聚焦PGA。通常在二维 IFFT 得到粗聚焦图像后选取若干个强散射点做 PGA 迭代估计残余相位然后回到波数域补偿。4.2 二维频域匹配滤波的编程实现与滤波参数说明PFA 在距离向已经做过脉冲压缩方位向匹配滤波实际包含在极坐标重采样的操作中但工程实现里往往还需要一个残余的频域滤波器来修正波数谱的幅相特性。下面给出一段完整的 PFA 聚焦函数包含 RVP 补偿和二维 IFFT。def pfa_focus(data_rect, K_r, freq_r, fc, c299792458.0): data_rect : 极坐标重采样后的二维频谱数据 freq_r : 距离频率轴 K_r : 调频率 (Hz/s) fc : 载频 (Hz) # RVP 补偿 rvp_phase np.exp(-1j * np.pi * freq_r**2 / K_r) data_rvp data_rect * rvp_phase[np.newaxis, :] # 二维 IFFT, 注意必须是 ifft2 而不是 fft2 img np.fft.ifft2(np.fft.ifftshift(data_rvp, axes1)) # 幅值图像取模, 保留复数数据用于 PGA return imgifftshift是要害。极坐标重采样之后零频位于数组中心如果直接用ifft2图像会发生整体偏移。ifft2和fft2在频域表示的零点位置不一样必须先ifftshift把零频挪到角落再做逆变换。rvp_phase的计算中freq_r必须是真实的频率值单位 Hz不能直接用 bin 索引换算否则补偿相位是错的。补偿完之后数据在距离维已经等效于理想脉冲压缩后的谱距离向分辨率完全由发射带宽决定。4.3 PFA 的聚焦深度与空变性参数表与散焦判据PFA 是空变的场景中心聚焦最好越往边缘散焦越严重。这个“有效成像场景尺寸”受限于极坐标重采样的近似条件。判断 PFA 是否适用可以用下面这张表来快速决策参数判定条件说明场景方位向尺寸 $L_s$$\frac{L_s^2}{8R_c} \le \frac{\rho_r}{4}$超过此值需要做距离向重采样距离向相对带宽 $B_r / f_c$$\frac{B_r}{f_c} \le 0.1$超过此值需要考虑高阶项改用子孔径或后向投影合成孔径累积角 $\Delta\theta$$\frac{\Delta\theta^2}{4} \le \frac{\rho_r}{4R_c}$累积角过大会引入显著的波前弯曲误差方位向分辨率 $\rho_a$$\rho_a \ge \frac{\lambda_c}{4\Delta\theta}$分辨率要求超过此界限时需增大合成孔径角但会加剧空变这张表的核心思想是 PFA 把回波信号的波前近似为平面波。当场景尺寸或累积角过大波前弯曲误差超过相位误差容限通常取 $\pi/4$PFA 就会失效。补救手段是将大场景划分为多个子块每个子块中心重新定义参考斜距和参考角度也就是 Subaperture PFA或者直接切换到后向投影算法BPA它在任何几何下都适用代价是计算量大幅上升。4.4 PFA 成像质量验证方法与失败特征对照拿到图像后最直接的验证是用一个角反射器目标布设在场景角落看它的响应函数。理想情况下点目标响应的峰值旁瓣比应该在 -13 dB 左右矩形窗加窗后更低。如果观察到的点目标响应沿距离向或方位向出现双峰、非对称旁瓣大概率是插值核不足或 RVP 补偿不对。我常用的验证流程是先跑点目标仿真验证 PFA 的核心链路正确性再跑真实数据用场景中的强散射点做 PGA 精聚焦。PGA 的输入是距离压缩后的复图像要求目标响应孤立且信噪比足够高。如果图像整体聚焦但存在区域性的模糊往往是运动补偿残余的空变相位这时要回到原始回波域检查平台轨迹拟合的阶数。轨迹拟合至少做到二次项对应加速度补偿拟合阶数不足会直接表现为图像方位向散焦且散焦程度随方位位置变化。5. PFA 的真实用法用自检验证聚焦、控制插值核长度和处理稀疏孔径的折中处理真实数据时我一般先在数据里找一个相位和幅度都稳定的强散射点比如铁路桥的金属护栏或者刻意布设的角反射器。这个点用于三个目的一是验证距离压缩后的峰值相位是否在方位向上平滑变化二是做 PGA 的种子点三是检查插值核长度设置是否合理。插值核长度是 PFA 里最难调的参数。单精度浮点下8 点 sinc 核能达到约 60 dB 的旁瓣水平而 4 点核只能到 30 dB 出头。对大多数成像场景8 点核够用但如果你做的是高分辨率星载 SAR距离向和方位向各做一次 8 点插值后两级误差累加可能让峰值旁瓣比恶化到 -20 dB 以上这时需要上 16 点核。核长了计算量线性上涨一个 16384 × 16384 的复图像用 16 点核插值单次插值的复数乘法次数在 40 亿次量级GPU 是必需品。还有一个常用技巧数据本身做过去调频之后距离向带宽通常不需要完整 FFT 点数。极坐标格式转换前先把距离频谱截断到实际信号带宽对应的通道数能显著降低插值次数。截断的位置要留 10% 的保护带防止频谱搬移时产生混叠。PFA 输出的是复图像相位信息保留完整。干涉应用中两幅 PFA 图像的相位差可以直接用于高程反演但前提是两幅图像的极坐标重采样网格完全一致否则相干性会因插值核不同而下降。做干涉测量时我通常固定插值核参数只改场景中心位置这样系统误差至少是相关的后续可以用干涉相位定标统一扣除。PFA 的复杂度集中在插值和几何标定但它的效率远高于后向投影是工程上做星载、机载聚束 SAR 成像时最常见的算法之一。本文还有配套的精品资源点击获取