：原理、Python实现与参数调优)
简介：面向雷达信号处理与逆合成孔径雷达成像研究者的自聚焦算法源码包，针对成像过程中因平台运动误差、相位失真等因素造成的图像散焦问题，提供相位梯度自聚焦算法的完整实现。压缩包体积约1KB，仅包含1个m格式的MATLAB脚本文件，代码围绕参数设定、相位误差估计、迭代校正、收敛判断等核心流程展开，结构清晰，既可直接运行验证效果，也可嵌入现有成像处理链路使用。已有1025人学习参考。通过研读这份脚本，读者可以快速理解算法从梯度计算到相位补偿的每一处实现细节，掌握在不依赖精确运动参数的前提下提升图像清晰度的经典自聚焦手段；同时可将其作为基线版本，用于算法改进、参数对比或课程实验。对于雷达信号处理、自动聚焦算法研究方向的初学者和工程师来说，是一份轻量但完整的参考实现，适合原理学习和二次开发。
1. ISAR自聚焦为什么绕不开PGA：相位误差把图像拖糊的硬现实
做ISAR成像的人都碰过这种情况：距离压缩做完了，包络对齐也做了，甚至运动补偿都按教科书走了一遍，成的图方位向还是糊的，像隔着毛玻璃看目标。你去查补偿参数，该校准的都校准了，可图像就是锐不起来。这时候真正的问题往往不在运动测量，而在回波里残留的相位误差。相位误差不像包络延迟那么好测，它对距离向几乎没有影响，却在方位向把能量抹开，让点目标变成一条条的拖尾。
自适应（autofocus）这类方法就是专门处理相位误差的，其中相位梯度自聚焦（Phase Gradient Autofocus，即PGA）是工程落地最稳的一支。它不依赖外部运动测量，不需要信标或参考点，直接从图像数据本身出发，把相位误差从强散射点的响应里估出来再补偿回去。这套方法的优点是：鲁棒、迭代次数少、适合ISAR这种目标散射点分布复杂、又没有人工参考点的场景。适合谁呢？做雷达信号处理、ISAR/SAR成像算法验证、以及想把实测数据聚焦质量拉起来的一线工程师和学生。下面按我自己的落地经验，把模型、实现、参数和坑一步步拆开。
2. PGA数学模型与ISAR数据适配：误差从哪来，输入怎么整理
2.1 相位误差的成像域表现：散焦、平移与方位模糊
ISAR成像通常先对距离压缩后的回波做包络对齐，再对每个距离单元的方位序列做FFT，得到距离-多普勒图像。理论上点目标在该图像里是一个尖锐的峰，峰值位置对应目标散射点的距离和横向位置。问题在于，回波的相位里如果不是理想的线性相位，而是叠加了一个随方位慢时间变化的误差项，FFT之后峰值就不再是一个点，而是被该误差的傅里叶变换撑开。
一个值得记住的结论是：恒定相位误差会让图像整体产生复常数偏移，不改变形状；线性相位误差带来整体平移；二次及更高阶项才造成散焦。机动目标或平台振动引入的往往是高阶项，所以散焦表现为点目标旁瓣抬升、主瓣变宽，严重时目标轮廓完全无法辨认。这也解释了为什么PGA只做相位校正，却能让图像清晰度发生质变——它等价于对每个方位时间的信号乘以一个估计出的复数校正因子，把高阶相位误差抹平。
2.2 从距离压缩数据到PGA输入：对齐、截断与去线性相位
PGA不是对原始回波跑的，它吃的是距离压缩之后、包络对齐之后的“准基带”数据。常见处理链是：脉压 → 包络对齐 → PGA → 方位FFT成图。包络对齐这一步非常关键，如果包络没对齐，PGA估计出的相位误差里会混入距离单元走动引起的线性相位项，导致校正方向错误。
我一般这样准备输入数据：
- 数据组织成二维复数矩阵 S[n, m]，n 为距离单元索引，m 为方位脉冲索引。
- 包络对齐算法用相邻相关法或全局最小熵法，先把整体时延搬平。
- 对每个距离单元做去均值，去掉直流分量，避免强静止杂波干扰后续选点。
- 可选做法是先把数据变换到图像域看一眼散焦程度，如果散焦实在太严重——比如主瓣宽度超过几十个单元——建议先做一次粗补偿（比如用特显点回波做相位梯度粗估），再进PGA迭代，否则第一次迭代就可能选到噪声上。
2.3 PGA迭代里的四个核心算子：选点、加窗、相位梯度估计、补偿
PGA的每次迭代由四个步骤组成，理解了这四个算子就能自己写实现。
第一步选点（dominant scatterer selection）。把 S[n, m] 沿 m 做FFT得到图像域 I[n, k]，在图像域找出每个距离单元的最大幅度位置，再选出幅度最大的若干个距离单元。这一步的本质是筛出信噪比高、相位展布真实的散射点，它们对相位误差的贡献是最诚实的。
第二步加窗（windowing）。对每个选中的距离单元，把其方位响应的峰值循环移位到图像中心，然后在峰值周围加一个窄窗，窗外的数据清零。窄窗把旁瓣和杂波截掉，让后续相位梯度估计不会被不相关的强散射干扰。窗口宽度是迭代的关键参数，后面第4章专门讲。
第三步相位梯度估计（phase gradient estimation）。对加窗后的数据做FFT回到方位时间域，然后相邻两个方位采样点共轭相乘，取辐角得到相位梯度。对整个窗内所有距离单元求平均，得到该次迭代的相位梯度估计值。
第四步补偿（correction）。把梯度积分得到相位误差函数，构造校正因子 exp(-j·φ(m))，乘到原始的 S[n, m] 上。然后重复迭代，直到梯度残差足够小或达到设定的最大迭代次数。
从这里可以看出，PGA的本质是把“估计相位误差”转化成“估计相位梯度”，因为梯度的统计估计比绝对相位的直接估计要稳定得多，这就是它比简单特显点测相法更耐噪声的根本原因。
3. 用Python跑通PGA自聚焦最小实现：从一维相位误差到聚焦图像
3.1 构造带相位误差的ISAR仿真数据
没有实测数据时，我习惯先用仿真数据验证算法流程。下面这段仿真构造了8个散射点，每个点有不同的距离和方位位置，并叠加了一个随方位时间变化的二次相位误差。注意这里故意只加了二次项，因为这是PGA最典型的应用场景。
import numpy as np def simulate_isar_data(n_range=128, n_pulse=256, n_scatter=8, phase_error_coef=0.8): """ 生成距离压缩后的ISAR回波数据（基带形式）。 参数说明: n_range: 距离单元数 n_pulse: 方位脉冲数（慢时间采样） n_scatter: 散射点数量 phase_error_coef: 二次相位误差系数，单位 rad/pulse^2 返回: data: (n_range, n_pulse) 的复数矩阵 true_phase: 真实的相位误差序列 """ t = np.arange(n_pulse) / n_pulse # 二次相位误差：系数越大散焦越严重 true_phase = phase_error_coef * (t - 0.5) ** 2 * n_pulse data = np.zeros((n_range, n_pulse), dtype=complex) rng = np.random.default_rng(42) for i in range(n_scatter): # 随机放置散射点，保证不同距离单元和方位位置 n_idx = rng.integers(5, n_range - 5) k_idx = rng.integers(0, n_pulse - 1) amp = 0.5 + 0.5 * rng.random() # 散射点贡献：距离向sinc旁瓣 + 方位向相位误差 sinc_r = np.sinc(np.arange(n_range) - n_idx) pulse_phase = 2 * np.pi * k_idx * t + true_phase data += amp * np.outer(sinc_r, np.exp(1j * pulse_phase)) # 加一点复高斯噪声 noise_std = 0.05 data += noise_std * (rng.normal(size=data.shape) + 1j * rng.normal(size=data.shape)) return data, true_phase这段代码里，散射点方位位置用k_idx决定，对应图像域里的横向位置；sinc_r模拟距离向点扩展函数；每个散射点的相位里都带同一个true_phase，所以整幅图像会被同一个相位误差函数散焦。噪声设为0.05量级，保证PGA还能正常工作。参数phase_error_coef决定误差强度，调试时可以从0.4逐级加到1.2，观察算法在散焦严重程度不同时的表现。
3.2 PGA核心迭代循环的实现
PGA的核心循环按前面说的四步展开。下面的实现保持了算法的模块化，方便你在自己的数据上替换输入。
def pga_autofocus(data, win_frac=0.25, max_iter=8, tol=1e-4): """ PGA自聚焦主循环。 参数说明: data: (n_range, n_pulse) 复数矩阵，已完成距离压缩和包络对齐 win_frac: 加窗宽度占方位FFT点数的比例，迭代中可自适应调整 max_iter: 最大迭代次数 tol: 相位梯度残差的收敛阈值，单位 rad 返回: data_corrected: 相位校正后的数据 phase_est: 估计出的相位误差序列 """ n_range, n_pulse = data.shape phase_est = np.zeros(n_pulse, dtype=float) # 每次迭代对当前相位校正后的数据重新选点加窗 data_work = data.copy() for it in range(max_iter): # 1. 变换到图像域并找最强散射点 img = np.fft.fft(data_work, axis=1) img_shifted = np.fft.fftshift(img, axes=1) amp = np.abs(img_shifted) # 每个距离单元的最大幅度位置及幅度值 peak_pos = np.argmax(amp, axis=1) peak_val = np.max(amp, axis=1) # 选幅度最大的前 25% 距离单元，且至少 8 个 n_select = max(8, int(n_range * 0.25)) select_idx = np.argsort(peak_val)[-n_select:] # 2. 对选中的距离单元做循环移位并加窗 # 窗宽随迭代进行可以由宽到窄，这里用固定比例 win_len = max(16, int(n_pulse * win_frac)) half = win_len // 2 phase_grad_sum = np.zeros(n_pulse, dtype=complex) for n_idx in select_idx: # 取出该距离单元的方位序列 row = data_work[n_idx, :] # 循环移位，使峰值位于中心 shift = peak_pos[n_idx] - n_pulse // 2 row_shifted = np.roll(row, -shift) # 加矩形窗，窗外置零 mask = np.zeros_like(row_shifted, dtype=bool) mask[max(0, n_pulse // 2 - half): n_pulse // 2 + half] = True row_windowed = row_shifted.copy() row_windowed[~mask] = 0 # 3. 相位梯度估计：FFT从图像域回到方位时域 # 相邻采样共轭相乘，取辐角 row_time = np.fft.ifft(row_windowed) grad_angle = np.angle(row_time[1:] * np.conj(row_time[:-1])) grad_complex = np.exp(1j * grad_angle) phase_grad_sum[1:] += grad_complex phase_grad_sum[0] += 1.0 # 平均梯度 grad_est = np.angle(phase_grad_sum) # 相位梯度积分得到相位误差 phase_new = np.cumsum(grad_est) phase_new -= phase_new.mean() # 4. 校正并更新相位估计 data_work *= np.exp(-1j * phase_new) phase_est += phase_new # 判断收敛：本次相位梯度修正量的RMS grad_rms = np.sqrt(np.mean(grad_est ** 2)) if grad_rms < tol: break return data_work, phase_est这段代码有几个实现细节值得注意。phase_grad_sum用复数累加而非直接对角度平均，是因为相邻共轭乘积本身是单位复指数，复数累加可以在平均时天然按幅度加权，避免角度平均在接近 ±π 跳变时产生错误结果。phase_est += phase_new是累积校正，每一轮迭代估计的是残余相位误差。循环移位用np.roll实现，配合对称窗可以避免截断造成的不连续。最大迭代次数我设置在8次，实际场景里3-6次基本就收敛，设多了也不会发散，只是浪费时间。
3.3 用图像熵和对比度量化聚焦效果
PGA跑完之后不能只看图，定量指标更重要。我最常用的两个指标是图像熵和对比度。图像熵越小代表能量越集中；对比度则对散焦更敏感，锐利的图像对比度显著更高。
def image_focus_metrics(img_2d): """ 计算二维复数图像的熵和对比度。 参数说明: img_2d: 复数图像域数据 返回: (entropy, contrast) 两个标量 """ amp2 = np.abs(img_2d) ** 2 total = amp2.sum() if total <= 0: return np.inf, 0.0 p = amp2 / total entropy = -np.sum(p * np.log(p + 1e-12)) / np.log(amp2.size) # 对比度：强度方差 / 强度均值 intensity = amp2 contrast = intensity.std() / (intensity.mean() + 1e-12) return entropy, contrast # 使用示例 data, true_phase = simulate_isar_data() data_fixed, est_phase = pga_autofocus(data) img_before = np.fft.fftshift(np.fft.fft(data, axis=1), axes=1) img_after = np.fft.fftshift(np.fft.fft(data_fixed, axis=1), axes=1) ent_before, con_before = image_focus_metrics(img_before) ent_after, con_after = image_focus_metrics(img_after) print(f"校正前: 熵={ent_before:.4f}, 对比度={con_before:.4f}") print(f"校正后: 熵={ent_after:.4f}, 对比度={con_after:.4f}")image_focus_metrics里的熵做了log(amp2.size)归一化，这样不同尺寸的图像之间可以比较。对比度用的是强度图的变异系数，散焦时旁瓣能量分散，对比度会明显下降。判断结果是否有效，一是看熵是否降低，二看对比度是否升高，两个指标方向一致才说明PGA起作用了。如果校正后熵反而变大，基本可以判定参数选得不对，往下看第4章。
4. PGA关键参数与选点细节：窗口宽度、阈值、迭代次数怎么定
4.1 加窗宽度：从宽到窄的迭代策略更稳
窗口宽度是整个PGA里最敏感的参数。窗开大了，窗内可能混入相邻散射点的旁瓣，相位梯度估计被污染；窗开小了，把有价值的散射响应截掉太多，估计噪声变大。我自己的经验是：第一轮迭代用宽窗，先捕捉大尺度的相位误差；随后每轮收窄，把精细结构修出来。
具体宽度一般用方位FFT点数的分数来表达，常见取值如下表：
| 迭代阶段 | 窗口宽度（占方位FFT点数） | 适用情况 |
|---|---|---|
| 第1轮 | 1/4 到 1/2 | 散焦严重，误差大 |
| 第2~3轮 | 1/8 | 误差已明显减小 |
| 第4轮及之后 | 1/16 到 1/32 | 精细校正 |
一种自适应做法是每轮迭代结束后，根据当前图像主瓣宽度重新估算窗口宽度。但如果图省事，固定使用win_frac=0.25也能稳定收敛，只是可能多迭代两轮。需要注意的是窗口宽度不能小于FFT点数除以目标横向尺寸对应的点数，否则截掉的不是杂波而是信号本体。
4.2 强散射点选择阈值与距离单元数目的平衡
选点时最常见的错误是只挑幅度最强的两三个距离单元，这样遇到目标上某个散射点特别强、其余散射点较弱的情况，PGA估计被强点“绑架”，弱目标区域的聚焦效果很差。正确做法是选足够多的距离单元，让相位梯度估计成为统计平均。
选择阈值两条经验：
- 幅度门限取最大峰值的 0.25~0.5 倍。低于0.25会把纯噪声距离单元纳入，高于0.5则选点太少。
- 距离单元数量至少15~30个。对128×256的小数据，选25%即32个距离单元是合理起点；对实测数据，距离单元数量多，可放宽到10%。
如果想进一步提升估计质量，可以在选点之后对每个距离单元的贡献按幅度加权。这样强点仍然主导，但弱散射点不至于完全没有投票权。加权方式和归一化幅度直接相乘即可，不需要复杂的自适应权重。
4.3 迭代次数与收敛判据：不要迷信固定轮数
很多工程实现把迭代次数定成固定值，比如10次，这是个隐蔽的坑。PGA在低信噪比条件下第1轮估计出的梯度可能含噪，第2轮往往最准；但如果信噪比极低，后续迭代可能在噪声上“过拟合”，越修越差。
我一般设两个停止条件，哪个先满足就停：
- 梯度残差RMS小于0.02~0.05 rad，说明相位误差基本被抹平。
- 连续两轮迭代的图像熵差小于0.1%，继续迭代收益很低。
# 修改 pga_autofocus 内部循环的停止判断示例 for it in range(max_iter): # ... 前面步骤不变 ... # 计算图像熵 ent_current = image_focus_metrics( np.fft.fftshift(np.fft.fft(data_work, axis=1), axes=1))[0] if it > 0 and abs(ent_current - ent_prev) < 0.001: break ent_prev = ent_current实测数据比仿真更复杂，第3轮之后经常出现熵曲线震荡，这时候不必强求完全收敛，取熵最小的那一轮结果重跑一次补偿即可。顺着这个思路，把每一轮的data_work存一份快照，最后选最优的来用，是一招很实用的“后悔药”。
4.4 分块处理与滑窗重叠：针对机动目标的参数调整
ISAR目标在成像积累时间内如果姿态变化明显，整段数据的相位误差不再能用同一个多项式描述。这时候常见的做法是分块——把方位向切成若干子块，每个子块独立跑PGA，再把校正后的数据拼起来。分块参数三个：子块长度、重叠率、块间相位对齐。
子块长度经验值取方位总点数的1/4到1/8，重叠率设50%左右。重叠的目的是保证块与块之间相位曲线的连续性，避免拼接处出现相位跳变。块间相位对齐则是把后一块的相位曲线整体平移，使其在重叠区域的相位均值与前一板块一致。实现上就是每块估计结束后记录相位曲线在重叠区的平均值，后续块差异对齐后再拼接。
这块挺玄学，不同雷达数据形态差异很大，我一般是先用1/4分块长度跑一版，看拼接处是否出现条带；有条带就减少子块长度加多重叠，直到边界消失。没有一组参数能通吃，这也是PGA调试里最耗时间但必须过的坎。
5. PGA常见问题避坑：相位解缠、强散射体与低信噪比场景
5.1 相位跳变导致图像出现横条纹
现象：校正后的图像方位向出现规律性横条纹，目标像被梳子梳过一样；或在图像边缘出现规则亮线。
原因：相位梯度估计值接近 ±π 时，np.angle的输出会从 π 跳变到 -π，积分后相位误差曲线出现一个2π的阶跃，补偿后等效于没补偿，反而引入周期性调制。
解决：主要靠两条。第一，对phase_grad_sum做幅度加权平均，强散射点贡献的对数幅度大，其梯度更可信，可以有效抑制噪声引起的相位跳变。第二，当检测到相邻梯度差值的绝对值大于π时，对该差值做加减2π修正，强行解缠。我在第3章代码里没有显式解缠，是因为仿真数据信噪比高；实测数据里这一步基本必须加。
# 对单距离单元的梯度做解缠修正 grad = np.angle(row_time[1:] * np.conj(row_time[:-1])) grad_unwrap = np.unwrap(grad)这里用np.unwrap比手动判断可靠，它默认的对相邻角度差值的π阈值正好适配相位梯度的特点。
5.2 强散射点主导，弱散射目标反而更糊
现象：PGA校正后，图像中最亮的点聚焦成尖锐峰，但周围的弱散射点区域比校正前更模糊，甚至出现虚假目标。
原因：选点时只选了幅度最大的少数距离单元，算法估计出的相位误差完全由强点决定。如果强点的相位误差特性和弱散射点不一致——比如强点是镜面反射，弱点是多次散射——强点的补偿量加到弱点上就是过补偿。
解决：把选点数量提高，同时引入幅度归一化后再累计梯度。另一种做法是先把最强点从数据里减除掉——用它的窗内估计重构该点响应并从原数据中扣除——再做一轮PGA，把次强点的信息挖出来。这种方法在ISAR动目标检测里很有效，但要注意减法后的残差不能引入额外偏差，重构响应要用加窗后的数据而非原始数据。
5.3 低信噪比下PGA不收敛，迭代反而恶化图像
现象：仿真数据信噪比降到5 dB以下时，PGA输出的图像熵不降反升，点目标周围出现大量碎斑。
原因：相位梯度的估计误差和信噪比成反比。信噪比低时，强的噪声单元被误选为散射点，相位梯度估计基本是随机的，补偿方向错了之后误差越积越大。
解决：第一个有效手段是限制迭代次数，低信噪比场景下2~3轮就停，不要追求梯度残差收敛到极小值。第二个手段是先用多视平滑预滤波，将邻近距离单元的幅度取平均后再选点，降低单单元噪声方差。第三是做特显点筛选的置信度判断——只保留在连续多轮迭代中峰值位置稳定的距离单元，抖动过大的直接丢弃。
5.4 多子块相位不一致，拼接后散焦
现象：分块PGA校正后，各块单独看都聚焦良好，但拼接起来目标边缘断裂，块与块之间有相位跳变的感觉。
原因：每一块独立估计的相位误差有一个整体常数偏移和可能的线性斜率差异。常数偏移不影响单块图像幅度，但拼接时块间数据同时在距离向上比，相位不一致就会在拼缝处产生相干抵消。
解决：分块处理时保留相邻块重叠区域，在重叠区内计算两块校正后数据的相位差均值，把后一块的数据整体乘一个相位校准因子。校准因子为exp(-j * mean(angle(block1 * conj(block2))))，重叠区越长，估计越稳，一般取50%重叠就能把拼接误差压到可忽略程度。
6. 把PGA嵌进ISAR完整处理链：验证方法、运动补偿衔接与进阶用法
6.1 用仿真的相位误差模型做定量验证
验证PGA实现是否正确，不能只靠看图。我在每次调试新数据前会先跑一个已知误差注入的仿真，量化PGA估计相位和真实相位的误差：
data, true_phase = simulate_isar_data() data_fixed, est_phase = pga_autofocus(data) # 估计相位和真实相位可能有整体常数差，去均值后比较 true_phase_c = true_phase - true_phase.mean() est_phase_c = est_phase - est_phase.mean() phase_mse = np.mean((true_phase_c - est_phase_c) ** 2) print(f"相位估计均方误差: {phase_mse:.6f} rad^2")相位误差的RMS小于0.05 rad基本可以认为实现正确。配合图像熵的前后对比，就能区分是算法问题还是参数问题。
6.2 PGA与包络对齐、初相校正的顺序关系
常见工程顺序是：距离压缩 → 包络对齐 → 初相校正（含PGA） → 方位FFT → 成像显示。PGA放在包络对齐后面是有道理的——包络对齐解决距离单元走动，PGA解决相位误差，两个问题基本解耦。如果先PGA后包络对齐，PGA会把包络错位引入的线性相位误差也当成待校正量，容易把距离单元搬乱。
对于机动目标，中间还要插一步：先粗包络对齐，PGA粗校正两轮，再做精包络对齐，最后PGA精校正。这样比一步到位的效果稳定很多。
6.3 多特显点加权PGA与子孔径PGA的取舍
实际工程的进阶方向基本是两个分支。一个是多特显点加权PGA（Weighted PGA），把每条距离单元的相位梯度估计按信噪比加权，适合实测数据中特显点密度不均的场景；另一个是子孔径PGA，把方位数据切短做多次PGA接力，适合目标大转角、相位误差随时间变化的场合。
从我自己的落地经验看，先固定用标准PGA把数据跑通，再做加权改造，收益比直接上复杂算法更高。我早期翻车最多的就是总想一步到位用加权PGA，结果参数太多，出了问题不知道怪谁。现在习惯是先标准PGA后看熵曲线，需要再往加权方向走。希望这个顺序能帮你在自己的数据上少走弯路，也希望这套参数和代码，能让你在下次被ISAR散焦折磨时，多一个顺手可用的工具。
本文还有配套的精品资源，点击获取
