简介:面向合成孔径雷达(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_r=None): """ 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, n=n_fft_r, axis=1), axes=1) return range_comp这段代码的核心在于fft(data_if, n=n_fft_r, axis=1)指定了沿快时间维做 FFT,fftshift把零频挪到中间,便于后续按距离频率索引。注意 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), dtype=complex) 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), dtype=complex) for j in range(n_rx): col = data_rect[:, j] real_interp = interp1d(k_y_old, col.real, kind='linear', bounds_error=False, fill_value=0) imag_interp = interp1d(k_y_old, col.imag, kind='linear', bounds_error=False, fill_value=0) 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, c=299792458.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, axes=1)) # 幅值图像取模, 保留复数数据用于 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 成像时最常见的算法之一。
本文还有配套的精品资源,点击获取