PGA自聚焦原理与MATLAB实现：ISAR相位误差校正实战

简介：在雷达成像领域，ISAR与SAR成像常因目标非合作运动而引入高阶相位误差，导致图像方位向散焦、轮廓模糊甚至出现重影。相位梯度自聚焦（PGA）作为经典的运动补偿算法，无需依赖孤立强散射点，通过循环移位、加窗和相位梯度估计的迭代流程，可有效校正平动与转动引起的相位污染。其技术价值在于兼具鲁棒性和工程可行性，广泛用于高分辨率雷达目标识别、舰船与飞机成像质量提升等场景。本文从相位误差来源讲起，深入PGA核心原理，给出可直接运行的MATLAB最小实现、参数选取经验与常见避坑指南，并结合图像熵与组合自聚焦技巧，帮助SAR/ISAR成像工程师快速掌握从散焦数据到清晰图像的完整处理路径。

1. PGA自聚焦：为什么ISAR图像模糊这件事绕不开它

当一组ISAR回波做完距离压缩和包络对齐，图像仍然在方位向糊成一片时，PGA自聚焦是我第一个想到的算法。PGA（Phase Gradient Autofocus，相位梯度自聚焦）是ISAR运动补偿里处理高阶相位误差最常用的工程方法。它不依赖预先挑选所谓“孤立强散射点”，也不假设目标转台模型精确成立，而是在方位向傅里叶变换后，通过相邻慢时间采样之间的相位差分，逐步逼近真实相位误差并补偿回去。我处理过不少X波段实测数据，最终让图像从“拖尾散焦”恢复到清晰轮廓，靠的就是这套迭代流程。这篇笔记写给做SAR/ISAR成像算法、雷达目标识别以及被散焦图像折磨的同学，从原理、最小实现到常见翻车点一次讲完。

2. ISAR成像里的相位误差从哪来：PGA要对付的是哪些项

2.1 一阶相位误差与方位散焦

ISAR成像的基本处理是把回波变换到距离-慢时间域，做完距离压缩后，每个距离单元内信号的相位可以写成一个沿慢时间变化的复指数序列。理想转台模型下，目标上某个散射点的回波相位近似为：

φ(t_m) ≈ (4π/λ)[ R0 + x_r·Ω·t_m + (y_r·Ω²/2)·t_m² + … ]

其中 t_m 是慢时间，λ 是载波波长，R0 是目标中心到雷达的参考距离，x_r、y_r 是散射点在目标坐标系里的位置，Ω 是等效转动角速度。

这里最关键的是 t_m 的高次项。一次项 (4π/λ)·x_r·Ω·t_m 对应散射点的多普勒中心频率，只会让图像整体在方位向平移，通常在做平动补偿时已经处理掉。真正让图像散焦的是二次项和更高次项。二次项相当于方位向存在一个线性调频率，如果成像处理时没有精确匹配这个调频率，点目标的方位响应就从一根谱线展宽成一个主瓣变胖的sinc函数。信噪比稍差时，整个目标看起来像是被横向抹了一层雾。

PGA 的特别之处在于，它并不试图单独估计这个调频率，而是把一次、二次、更高次相位误差当作一个整体，用相位梯度的方式逐轮迭代消除。这意味着它对模型失配的容忍度比“估计一个二次函数再补偿”的做法高得多。

2.2 二阶及高阶误差与图像重影/畸变

如果相位误差只有二次项，图像的主要问题是方位向散焦，至少目标的主体结构还在，只是轮廓变模糊。但ISAR目标通常是飞行器、舰船这类非合作目标，等效转动角速度 Ω 可能随姿态变化而改变，甚至目标自身存在振动部件。这种非均匀转动会在相位上引入三次、四次项，甚至周期性抖动。

高阶相位误差带来的图像效应不再是对称展宽，而是主瓣位置偏移、旁瓣不对称、强散射点周围出现“彗尾”状伪影。实际数据里最常见的表现是：做完常规运动补偿后，图像中心区域看似有一团亮斑，但目标边缘轮廓撕裂，几个强散射点之间出现一连串假目标。很多人把这种现象误判为“分辨率不够”，其实问题出在相位误差没有补偿干净。

PGA 选择的相位梯度估计方式，天然能适配上述误差中的慢变部分。对周期性抖动这类高阶快变误差，PGA 也能通过多轮迭代逐步压缩，只是收敛速度和初始窗宽的选择会更苛刻。

2.3 包络对齐和相位补偿的分工

ISAR 运动补偿通常拆成两段：包络对齐和相位补偿。包络对齐解决的是“目标回波在哪个距离单元”的问题。由于目标运动使每次脉冲的回波延迟不断变化，必须先把所有脉冲的包络对齐到同一距离位置，否则每个距离单元里的信号来自目标的不同部位，后续相位处理无从谈起。包络对齐最常用的是相邻互相关法或最小熵全局对齐法，精度通常能做到子距离单元量级。

但包络对齐只能把回波“搬”到同一个距离门，无法修正每个脉冲内部的载波相位。相位补偿必须达到波长量级精度，PGA 做的就是这后半段。它要求输入的复数数据矩阵已经完成包络对齐，且每个距离单元内的信噪比足够支撑相位估计。如果包络对齐残差太大，比如超出半个距离单元，PGA 会把包络错位误当成相位误差来补偿，结果图像距离向出现撕裂状条纹。这不是 PGA 的算法缺陷，而是输入数据没达到它的工作前提。

3. PGA自聚焦核心原理：循环移位、加窗与相位梯度估计

3.1 为什么PGA不挑特显点

传统自聚焦方法里，特显点法（也叫强散射点法）依赖在图像里挑选一两个幅度明显高于背景的孤立点，用它们的相位差来估计误差。这个方法在特显点稳定时很快，但一旦特显点被旁边的散射体干扰，或信噪比下降，估计结果会出现跳变。实测 ISAR 目标往往没有那么多理想孤立强点，尤其是舰船目标和隐身外形目标，散射分布很碎。

PGA 的鲁棒性来源在于：它同时使用几十个距离单元的信号，对每个距离单元先做循环移位，把最强散射点搬到孔径中心，再用幅度加权平均来估计相位梯度。幅度加权意味着信号强的距离单元贡献大，信号弱、噪声主导的距离单元贡献小。即使个别距离单元内没有强散射点，它们也不会把估计结果拉偏太多。这种“多个单元联合、按信噪比加权”的策略，让 PGA 在低信噪比和散射点密集的场景下依然可用。

3.2 三步核心操作的物理含义

PGA 每次迭代做三件核心操作：

首先是循环移位。把当前图像里每个距离单元的方位向幅度峰值移到孔径中心。为什么要移？因为散射点在目标上的位置不同，它们本身就有各自的多普勒频率，这个一阶项会污染相位梯度估计。循环移位相当于把每个距离单元的“局部最强散射点”先对齐到零频位置，消掉它们各自的线性相位项，剩下的就是真正需要补偿的相位误差。

其次是加窗。移位后的图像会在孔径中心附近形成一个能量集中的亮区。在这个区域外，主要剩下噪声和旁瓣。加窗就是只保留中心附近一个范围内的数据，把窗外数据置零。窗宽的选择很讲究，太宽会把噪声放进来，太窄又会把真实信号截断。PGA 的做法是让窗宽随迭代递减，第一轮用比较宽的窗，把大尺度散焦能量都包进来；后面几轮逐步收缩，聚焦越来越好，散射能量越来越集中，窗宽也可以越来越窄。

最后是相位梯度估计。把加窗后的频域数据做逆傅里叶变换回到慢时间域，然后对相邻两个慢时间采样做共轭相乘。假设真实相位误差是 φ(t_m)，共轭相乘后的相位就是 φ(t_m) - φ(t_{m-1})，也就是相位误差在相邻脉冲之间的变化量，即梯度。对所有参与估计的距离单元按幅度加权平均，得到当前迭代的相位梯度估计值。

3.3 迭代与窗宽收缩：PGA的收敛逻辑

单靠一次循环移位和加窗很难把相位误差全部估计出来。因为第一轮窗宽大，估计出的梯度是宽频带内误差的平均效果，只能补偿掉主要的大尺度散焦。补偿之后残余误差变小，散射能量更集中，这时再把窗宽收缩一轮，就能估计出更高频的残余相位分量。如此循环迭代，就是一个由粗到细的逼近过程。

我习惯把 PGA 的迭代流程写成六步：

1. 方位向 FFT，得到当前图像。
2. 对每个选中的距离单元做循环移位，把幅度最大值搬到中心。
3. 按当前窗宽加窗，滤除窗外噪声。
4. 逆 FFT 回慢时间域，相邻采样共轭相乘，幅度加权求相位梯度。
5. 对梯度做积分得到相位误差，取共轭乘回数据做补偿。
6. 更新图像，减小窗宽，回到第 1 步。

需要提醒的是，这个循环里“积分得到相位误差”并不是单纯累加梯度就完事。每次估计出的梯度要去掉一个常数偏置再积分，否则循环移位带来的整体相位偏移会被错误地当作真实误差补偿进去，导致图像在方位向整体搬移。这个细节在新手实现里非常容易漏。

4. 用MATLAB跑通PGA自聚焦：最小实现和参数选择

4.1 输入数据约定和主循环代码

先约定输入格式。这是一个大小为 [N_rng, N_az] 的复数矩阵，N_rng 是距离单元数，N_az 是方位向脉冲数。数据必须已经完成距离压缩和包络对齐，每个距离单元内是沿慢时间排列的复基带信号。

下面是一版可以直接用到自家数据上的 MATLAB 实现，核心循环保留，注释写在实际动手时容易踩的位置。

function [s_out, phi_est] = pga_autofocus(s, rbins, max_iter, w0) % PGA 相位梯度自聚焦 % s: 距离压缩+包络对齐后的复数矩阵, [N_rng, N_az] % rbins: 参与估计的距离单元序号, 传空数组则自动挑能量最大的32个 % max_iter: 最大迭代次数 % w0: 初始窗宽(方位向点数), 默认取全孔径 % 返回: % s_out: 相位误差校正后的复数矩阵 % phi_est: 估计出的相位误差向量, 长度 N_az [N_rng, N_az] = size(s); if nargin < 4 || isempty(w0), w0 = N_az; end if nargin < 3 || isempty(max_iter), max_iter = 8; end if nargin < 2 || isempty(rbins) energy = sum(abs(s).^2, 2); [~, order] = sort(energy, 'descend'); rbins = order(1:min(32, N_rng)); end s_corr = s; % 工作副本 phi_est = zeros(1, N_az); % 累积相位误差 ctr = floor(N_az/2) + 1; % 方位向中心索引 for it = 1:max_iter % 1. 方位向成像 img = fft(s_corr, [], 2); img_rb = img(rbins, :); % 2. 循环移位: 每个距离单元的最大值移到中心 [~, max_pos] = max(abs(img_rb), [], 2); shift = ctr - max_pos; % circshift 向右移 shift 个点 img_fix = zeros(numel(rbins), N_az); for rr = 1:numel(rbins) img_fix(rr, :) = circshift(img_rb(rr, :), [0, shift(rr)]); end % 3. 加窗: 只保留中心附近 w0 个点 win = zeros(1, N_az); half = floor(w0 / 2); left = max(1, ctr - half); right = min(N_az, ctr + half - 1); win(left:right) = 1; img_win = img_fix .* win; % 4. 逆FFT回慢时间域, 相邻采样共轭相乘估计相位梯度 data = ifft(img_win, [], 2); prod = data(:, 2:end) .* conj(data(:, 1:end-1)); ang = angle(prod); % 幅度加权: 幅度大的采样点更可信 wgt = abs(data(:, 1:end-1)) .* abs(data(:, 2:end)); grad = sum(ang .* wgt, 1) ./ (sum(wgt, 1) + eps); grad = grad - mean(grad); % 去常数偏置 % 5. 积分得到相位误差并补偿 phi = cumsum(grad); phase_mat = repmat(phi, N_rng, 1); s_corr = s_corr .* exp(-1j * phase_mat); phi_est = phi_est + phi; % 6. 收缩窗宽 w0 = round(w0 * 0.5); if w0 < 2, w0 = 2; end end s_out = s_corr; end

这段代码的逻辑对应前面讲的六步流程。第 4 步里grad = grad - mean(grad)是必须的，它去掉的是循环移位带来的整体线性相位偏置；如果不做这一步，最终图像会在方位向整体平移，看起来好像聚焦了，但目标位置对不上。第 5 步直接把相位误差乘回全矩阵，不需要对每个距离单元单独做 FFT，因为相位误差在慢时间维是公共的，一次相乘即可完成校正。

4.2 特显点距离单元的挑选与加权方式

函数里默认挑了能量最大的 32 个距离单元参与估计。这个数量的选择是有讲究的：太少，统计平均不够，估计出的相位梯度噪声大；太多，会把大量以噪声为主的低能量距离单元也包括进来，反而拖低信噪比。实测中如果目标回波只占十几个距离单元，那就全部参与，不必强行凑 32 个。

挑选时只按能量排序不够稳。我一般会在排序前先看一眼距离压缩后的数据，把明显受干扰、整行都有连续强杂波的距离单元剔除，然后再交给 PGA 处理。有些实现还会检查每个距离单元方位向幅度分布的尖锐度，只选那些“主瓣越窄越好”的距离单元，因为这类单元内往往存在单个主导散射点，相位梯度估计更稳定。

幅度加权那里，wgt 用的是相邻两个慢时间采样幅度的乘积。这个权重设计对应最大似然意义下的最优加权，前提是噪声近似高斯白噪声。如果你的数据里存在射频干扰或脉冲性干扰，加权会被少数异常采样点带偏，这时候可以把 wgt 做一次截断，比如超过三倍中位数的权重全部压平，能有效提升鲁棒性。

4.3 迭代次数与窗宽：怎么定，怎么收

PGA 的性能高度依赖两个参数：初始窗宽和窗宽收缩节奏。初始窗宽最好由散焦程度估计，而不是盲目取全孔径。如果第一轮图像里目标方位向铺开的宽度约占整个孔径的 40%，那初始窗宽取 0.5 倍孔径是合理的。窗宽收缩用 0.5 衰减系数是经典做法，但实际数据里这个值常常偏激进，后文避坑部分会展开。

迭代次数上，我见过有人固定跑 10 轮，结果后面几轮窗宽只剩两个点，估计的相位梯度基本是噪声，把已经聚焦的图像又搞出新的随机相位。更稳的做法是每轮计算当前图像的熵，熵不再明显下降就停止迭代，这个方法在最后一章给代码。

5. PGA自聚焦避坑：5个能让你白跑一夜的常见问题

5.1 包络对齐残差把PGA带偏

现象：PGA 跑完以后，图像距离向出现斜纹或撕裂状结构，目标轮廓像被“搓”了一下，而不是整体变清晰。

原因：包络对齐残差超过半个距离单元时，相位误差不再只是慢时间的公共函数，而是随距离单元变化。PGA 用同一个相位向量去补偿全图像，自然会把一部分距离单元补对、另一部分补错。我在一次实测数据处理里，包络对齐残差约 0.3 个距离单元，PGA 第一轮迭代后图像熵反而上升，就是这个原因。

解决：回到包络对齐阶段，先用相邻互相关法做一个粗对齐，再用最小熵全局法做精对齐，把残差压到 0.1 个距离单元以内再进 PGA。检查方法也很简单，对齐完成后把每个距离单元的包络峰值位置画出来，若曲线有明显跳变，就说明对齐没做完。

5.2 窗宽收敛太快导致误判强散射体

现象：迭代两三轮后，相位梯度估计值不再变化，但图像还没完全聚焦，甚至出现新的模糊。

原因：窗宽按 0.5 系数每轮减半，第二轮窗宽只剩四分之一孔径，可能把真实散射体的一部分截在窗外。窗外数据被置零后，参与相位梯度估计的信号能量大幅下降，噪声相对占比变高，估计结果自然不可信。

解决：把衰减系数从 0.5 改成 0.75 左右，或者用自适应判据：只有当窗内能量占整段数据总能量的比例超过 70% 时才允许继续收窗，否则保持当前窗宽再迭代一轮。

5.3 循环移位把强散射体中心移错

现象：某一轮迭代估计出的相位梯度出现大幅跳变，相位误差曲线不连续。

原因：循环移位直接用max(abs(img))找最大值位置，如果这个最大值是旁瓣、噪声尖峰或两个散射体干涉的结果，移位后中心位置并不对应真实散射点，后面所有相位梯度估计都建立在错误的位置上。

解决：改用局部质心估计。先以最大值位置为中心，取左右各几个点，计算这些点的幅度平方加权重心，再把重心位置作为移位目标。我在代码里通常这样处理：

hw = 3; idx = max(1, mi-hw):min(N_az, mi+hw); weight = abs(img_rb(rr, idx)).^2; cen = sum((idx - mi) .* weight) / sum(weight); shift(rr) = ctr - round(mi + cen);

注意：质心估计只取局部点，不要取整个孔径，否则会把旁瓣能量也拉进来。实测中这一步能把移位精度从“一个采样点”提升到“亚采样点”，相位梯度估计的跳变概率明显下降。

5.4 加权平均被少数大目标支配

现象：图像里最强的那一两个散射点聚焦得不错，周围弱散射点的清晰度反而变差，整体图像熵没有同步下降。

原因：幅度加权平均本质上是用幅度平方做权重，强散射点贡献占比过高，相当于整个相位梯度估计被少数距离单元垄断。如果这些强距离单元内部本身的相位误差与全局不一致，误差就被放大。

解决：对参与估计的距离单元先做幅度归一化，每个距离单元的采样信号除以该单元的最大幅度，再按归一化后的幅度加权。这样每个距离单元的贡献被限制在相似量级，弱散射点也能参与“投票”。更进一步可以用截断权重，把超过中位数三倍的权重直接压到三倍中位数，鲁棒性更好。

5.5 高阶相位误差出现“假聚焦”

现象：迭代后图像熵在下降，但目标形状明显失真，强点周围出现不对称的“彗尾”，而不是干净的旁瓣结构。

原因：PGA 的相位梯度估计对慢变相位误差效果好，但对快速抖动、周期性振动这类高阶误差响应不足。窗宽收缩到很窄以后，窗内剩余能量主要来自高频相位分量，此时梯度估计逼近的是这部分误差，但积分后无法完全还原真实的非线性相位轨迹。

解决：改用子孔径 PGA。把方位向脉冲序列分成四段，每段独立做 PGA，得到四段相位误差曲线，再按重叠区域相位差最小原则拼接成一条完整曲线。分段长度要保证每段内目标转动变化仍近似均匀，一般取总孔径的 1/4 左右。这个方案对非均匀转动目标效果改善非常明显，也是我在处理舰船目标时的默认升级项。

6. 把PGA用到准：收敛判据与组合自聚焦技巧

6.1 用图像熵判断什么时候该停

固定迭代次数不够稳，我习惯在迭代循环里记录每轮图像的幅度熵。定义归一化幅度图像的熵为：

E = -Σ p(m,n) · ln p(m,n)，其中 p(m,n) = |I(m,n)|² / Σ|I(m,n)|²

每一轮迭代完成后计算一次熵，连续两轮的熵变化率小于 0.1% 时就停止迭代。原因是：图像越聚焦，能量越集中，熵越低；熵不再下降，说明残余相位误差已经小到对成像质量没有实质影响，继续迭代只会增加计算量甚至引入噪声。

6.2 最小熵法粗聚焦 + PGA精修的组合打法

纯 PGA 对大的二次相位误差收敛较慢，尤其当初始调频率失配严重时，前两轮窗宽大、估计精度低，浪费了很多迭代。工程上更稳的做法是：先用最小熵法或对比度最大化法估计方位向调频率，做一个粗聚焦，把散焦宽度压缩到主瓣的几倍以内，然后把粗聚焦后的数据交给 PGA 做精修。

这组组合的好处是分工明确：粗聚焦解决“大散焦”问题，PGA 解决“残余高阶误差”问题。实测数据里，粗聚焦加上 3~4 轮 PGA 通常就能达到纯 PGA 跑 8 轮的效果。

6.3 用注入误差验证PGA的性能边界

在把 PGA 结果交给下游识别方才放心，我最后都会做一次注入验证。对同一份数据人为注入一个已知的二次相位误差，跑完 PGA 后比较估计相位与注入相位的残差。

t_m = (0:N_az-1) / prf; % 慢时间轴线 k_qpe = 8; % 人为二次相位系数 s_test = s .* exp(1j * pi * k_qpe * t_m.^2); [~, phi_out] = pga_autofocus(s_test, [], 6, N_az); residual = wrapToPi(phi_out - pi * k_qpe * t_m.^2); rms_err = sqrt(mean(residual.^2));

如果 RMS 残差小于 5 度，说明这条数据链路的信噪比和 PGA 参数都没问题。如果残差明显偏大，我会回头检查包络对齐质量、参与估计的距离单元数，以及初始窗宽是否太激进。这套验证流程帮我挡掉过不少“看似聚焦实则虚焦”的坏结果，也让我对 PGA 的实际鲁棒性心里有底。希望帮到你。

本文还有配套的精品资源，点击获取