PGA自聚焦算法源码解析：MATLAB实现ISAR图像相位误差校正

简介：基于惩罚函数梯度（PGA）的自聚焦算法程序，服务于逆合成孔径雷达图像处理，面向雷达信号处理方向的研究者与工程师，可解决因系统参数误差、目标运动不确定等因素造成的图像散焦问题。压缩包内仅含一个MATLAB脚本文件，大小约一千字节，代码结构完整，涵盖参数初始化、图像质量指标计算、梯度求解、相位校正因子迭代更新以及收敛判断等关键步骤，可直接运行或嵌入成像处理流程。目前已有一千零二十五人学习下载，属于小而精的算法示例。通过学习这份代码，读者能理解自聚焦算法如何通过梯度信息优化相位校正因子，提升逆合成孔径雷达图像的清晰度。在此基础上，可进行二次开发，服务于自身的科研或工程项目。

1. PGA自聚焦源码包：给ISAR图像散焦上一道保险

做ISAR成像的人多半都遇到过这种场面：距离-多普勒矩阵拉出来，目标轮廓有了，方位向却糊成一片，像没对准焦的照片。PGA自聚焦就是专门处理这个问题的。这个源码包的核心是一份PGA.m，用MATLAB实现相位梯度自聚焦算法，思路是估计回波信号里的相位误差，再做补偿，把散焦的ISAR图像重新拉清晰。它解决的是运动补偿之后残余相位误差导致的图像质量劣化，适合做雷达信号处理、ISAR成像算法验证的工程师和研究生。下面按照原理、代码拆解、参数整定、踩坑记录一路讲下来，最后给一套验证方法。

2. ISAR相位误差模型：为什么PGA的“梯度思路”能去散焦

把PGA跑起来之前，得先建立相位误差的数学模型。ISAR成像的基本流程是对回波做距离压缩、运动补偿、方位向FFT，理想情况下每个散射点应该是一个冲激响应。但实际数据里，目标不是匀速转动的，平台自身的抖动、信道衰落、系统参数不精确，都会让散射点的相位偏离理想值，结果就是方位向主瓣展宽、旁瓣抬高，图像看起来是糊的。

2.1 相位误差从哪来

ISAR本质上是在小转角内观测目标，方位向的分辨率靠多普勒频率差来实现。运动补偿通常只能去掉一阶项（平动），二阶以上的残余运动会导致相位误差随方位向慢变化。另一个来源是脉冲间系统的时钟抖动和本振相位噪声，这个误差在慢时间维度上叠加。

相位误差的数学表现形式是：每个距离单元的回波序列在方位向上被乘以一个复指数相位误差项exp(jφe(n))，n是方位向脉冲序号。自聚焦要做的就是从距离-多普勒域图像里估计φe(n)，然后乘以exp(-jφe(n))，把散焦搬回聚焦。

这里要提一句：有些资料把PGA全称写成Penalty Function Gradient，但ISAR工程圈的通用叫法是Phase Gradient Autofocus，也就是相位梯度自聚焦，重点在于利用相位误差在距离单元间的梯度结构来估计相位，而不是什么惩罚函数。这篇内容按工程惯例处理。

2.2 相位梯度自聚焦的核心思想

PGA的出发点是：目标里总有几个强散射点，它们散焦之后能量会展布到多个方位单元。如果能找到这些强点的响应，把它们的相位梯度提取出来，就能还原真实相位误差的差分序列。

实现上分四步：第一步，把强散射点所在的距离单元在方位向做循环移位，让峰值移到图像中心，等效于去掉常数相位项；第二步，加窗只保留主瓣和少量旁瓣，抑制噪声和邻近散射点干扰；第三步，用最大似然估计提取相邻方位单元的相位差，得到相位梯度；第四步，对梯度做积分还原相位校正曲线。这一步思路的关键是：散焦导致的相位误差在慢时间维度上是平滑变化的，梯度估计天然对噪声鲁棒。

实际工程中，PGA一般迭代运行。每一轮迭代先用当前相位校正曲线补偿回波，重新成像，再做循环移位、加窗、梯度估计，直到相位梯度足够小。迭代的好处是，估计精度会随补偿越来越准，强散射点也会被聚焦得越来越锐利。

2.3 不同自聚焦算法怎么选

ISAR自聚焦算法常见的有对比度最大法、最小熵法、子空间法和PGA。对比度最大法把自聚焦变成一个最优化问题，用图像对比度做代价函数，效果不错但计算量大，而且容易陷入局部极值；最小熵法类似，熵值作为代价函数；子空间法精度高但要在快拍维度上求协方差矩阵，数据量小的场景不稳定。

PGA的优势在于：它不需要对整幅图像做非线性搜索，每一轮迭代的开销主要在FFT和复数乘加运算上；通过加窗把污染点挡在外面，稳健性也够。缺点是强散射点不明显时估计方差会变大，所以很多工程实现会先用对比度最大法做一轮粗补偿，再用PGA精细补偿。常见做法是混合使用，这个后面讲代码时再展开。

3. PGA.m逐段拆解：从距离-多普勒数据到聚焦图像

拿到了源码包，先别急着跑，把PGA.m拆开看。它的输入一般是距离-多普勒数据矩阵，行对应距离单元，列对应方位向脉冲，输出是补偿后的聚焦图像和估计出来的相位校正序列。这一章按函数骨架、循环移位、加窗、梯度估计、迭代更新五个模块来拆。

3.1 主函数骨架与输入输出约定

主函数签名常见做法是这样：

function [img_focused, phi_est] = pga_autofocus(img, opts) % img: 二维复数矩阵，距离-多普勒域数据，尺寸 [Nr, Na] % opts: 结构体，包含迭代次数、窗口宽度系数、收敛阈值等 % img_focused: 相位补偿后的ISAR聚焦图像 % phi_est: 估计的相位校正序列，长度 Na，单位弧度 if nargin < 2 opts = struct('iter_max', 15, ... 'win_coef', 0.2, ... 'tol', 0.01, ... 'n_scatter', 12); end Na = size(img, 2); phi_est = zeros(1, Na); img_cur = img; for iter = 1:opts.iter_max % 每轮迭代的核心步骤 [grad, idx] = estimate_phase_gradient(img_cur, opts); if isempty(grad) break; end phi_est = phi_est - grad; % 负梯度方向修正 comp = exp(-1i * phi_est); img_cur = img .* repmat(comp.', size(img, 1), 1); % 收敛判断 if norm(grad) / sqrt(Na) < opts.tol break; end end img_focused = img_cur; end

这段代码的核心逻辑是把每次估计出来的相位梯度累加到总校正序列phi_est上，再用负梯度方向更新补偿。这里有几个细节：img_cur每一轮都基于原始img乘上累积补偿，而不是在上一轮结果上继续乘，这样避免误差逐轮累积放大。收敛判断用的是梯度的均方根值，比用最大绝对值更稳。

参数方面，iter_max控制最大迭代次数，一般取10到20就够；win_coef控制加窗宽度，默认0.2意味窗口占方位向总长度的20%；tol是梯度收敛阈值，单位是弧度；n_scatter是参与估计的强散射点数量。这些参数在第四章还会细说。

3.2 循环移位与加窗

循环移位是PGA里容易被忽略的步骤。强散射点散焦后分布在多个方位单元，但它的峰值位置不一定在图像中心。如果不把峰值移到零频位置，后面加窗的时候就会把主瓣削掉一半。循环移位的代码实现并不复杂：

function img_shifted = circular_shift(img, idx_peak) % 将每个距离单元沿方位向循环移位，使峰值位于方位中心 % idx_peak: 每个距离单元强散射点的峰值索引 [Nr, Na] = size(img); img_shifted = zeros(size(img)); for r = 1:Nr shift = idx_peak(r) - round(Na/2); img_shifted(r, :) = circshift(img(r, :), -shift); end end

注意是负号：如果峰值在Na/2的右边，shift为正，circshift负移把峰值拉回左边到中心。循环移位并不会改变信号的幅值谱，只改变线性相位项，这正是我们想要的，因为线性相位项对应目标在方位向上的位置，不影响聚焦质量。

加窗是为了把目标距离单元里其他散射点和噪声的影响排除掉。散焦后的主瓣宽度和相位误差大小有关，窗口宽度不能固定。我一般会取方位向多普勒谱峰附近主瓣3dB宽度的2到3倍，上限设为图像宽度的一半，避免把有效信号全截掉。实现上用汉明窗：

win_len = round(Na * opts.win_coef); win = hamming(win_len).'; pad_left = max(1, round(Na/2) - floor(win_len/2)); sel = pad_left : pad_left + win_len - 1; img_shifted(:, ~ismember(1:Na, sel)) = 0;

这里采用只保留窗口内数据的做法，窗口外的方位单元直接置零。如果不用窗函数而是用矩形截断，旁瓣会比较高，相位梯度估计受干扰明显。汉明窗在抑制旁瓣和保留主瓣能量之间比较均衡。

3.3 相位梯度估计

相位梯度估计是PGA的核心，思路是取加窗后的相邻方位单元做共轭相乘，提取相位差。最大似然估计的实现方式是对多个距离单元的相位差加权平均：

function grad = estimate_phase_gradient(img_w, idx_strong) [Nr, Na] = size(img_w); % 只选取能量最强的几个距离单元 power_r = sum(abs(img_w).^2, 2); [~, sort_idx] = sort(power_r, 'descend'); idx_use = sort_idx(1:min(length(idx_strong), Nr)); accum = zeros(1, Na-1); for r = 1:length(idx_use) s = img_w(idx_use(r), :); diff_prod = s(1:end-1) .* conj(s(2:end)); accum = accum + angle(diff_prod .* abs(diff_prod).^2); end grad = accum / (length(idx_use) + 1); grad = [0, grad]; % 首位补零，长度对齐Na end

为什么用angle(diff_prod .* abs(diff_prod).^2)而不是直接用angle(diff_prod)？把幅度平方加权进去，可以让强散射点对梯度估计有更大话语权。角度提取放在加权之后，这个顺序不能反，因为angle操作是非线性的，权重要作用在复数量上才有意义。相邻方位单元相位差的物理含义是相位误差的一阶差分，对差分做累加恢复相位曲线时，DC分量也就是常数相位，不影响图像聚焦，所以补零可以从第一个单元开始。

3.4 迭代更新与收敛检查

相位校正序列的更新公式是phi_est = phi_est - grad。为什么用负梯度方向？因为估计出来的相位差来自散焦信号，散焦的相位误差方向和校正方向恰好相反，减掉梯度等效于在对相位误差做最速下降。迭代过程中要注意相位卷绕问题，grad是从复数角度提取的，值域范围在[-pi, pi]，但累加后的phi_est可能超过这个范围，所以在乘exp(-1i*phi_est)时MATLAB会自动按复数模运算处理，不会出错。

收敛检查还有一层更稳健的做法，就是用连续两轮图像熵的变化量来判断：

entropy_prev = compute_entropy(img_cur); % 迭代体 entropy_new = compute_entropy(img_cur); if abs(entropy_new - entropy_prev) < 1e-3 break; end

图像熵比梯度范数更贴近“看起来清楚不清楚”这个直观感受，如果梯度在跳动但熵值已经稳定，说明图像质量不再有提升，继续迭代只会浪费时间。我一般会把这两种判断都写进代码里，两者满足其一就退出循环。

4. 参数整定：窗口宽度、迭代次数、收敛阈值怎么配

PGA的参数不是越激进越好。窗口宽度、迭代次数、强散射点个数、收敛阈值这四组参数，直接决定了算法的收敛速度和最终图像质量。这一章针对不同数据特征给一组调参经验。

4.1 加窗宽度：主瓣与估计偏置的平衡

加窗宽度是最敏感的参数。窗口太宽，会把干扰散射点和噪声放进估计区间，导致相位梯度估计方差变大；窗口太窄，主瓣被削掉一部分，强散射点的能量损失，估计结果出现偏置。经验上，窗口宽度对应散焦主瓣宽度的1.5到3倍比较合适。

判断窗口宽度是否合适的办法很简单：看第一轮迭代后图像熵是降还是升。如果熵值在第一轮就明显上升，多半是窗口太宽把旁瓣干扰放进来；如果是收敛很慢、十几轮还在波动，多半是窗口太窄。另一个技巧是先跑一轮对比度最大法得到粗聚焦图像，再从粗聚焦图像里量主瓣宽度，反推窗口宽度，这个做法在工作站上跑几秒钟就能出来。

4.2 迭代次数与收敛阈值

迭代次数一般设置在5到15之间。少于5轮，相位梯度估计还没有收敛到稳定值，图像提升不明显；多于20轮，边际收益递减，而且会把估计噪声慢慢积累进相位校正序列。阈值tol取值在0.005到0.02弧度之间。取0.01作为默认值，能在收敛速度和精度之间取得平衡。

注意阈值太小会有问题：梯度中总有残留噪声，阈值太小会让算法进入震荡状态，表现为图像熵在某个值附近来回跳动。如果观察到这种情况，把阈值调大一个数量级，或者增加一个早停条件，比如连续三轮梯度范数变化小于5%，就强制退出。

4.3 强散射点选取数量

参与梯度估计的距离单元数，影响最大的是估计方差。选取的数量太少，比如只有2到3个强点，遇到某个距离单元里散射点分布不佳时，估计会失控；选择太多，比如超过距离单元总数的三分之一，会把很多能量不高的弱散射点也卷进来，估计反而不准。

实际做法是按下式选择参与单元：先按距离单元功率降序排列，取前10个，再删掉那些峰值不在方位向中心附近的单元。PGA算法有效的前提是选中的距离单元里有明显的孤立强散射点，如果整个目标都是分布式散射，PGA的效果会很差，需要用其他自聚焦算法。这个在第五章里还会细说。

4.4 典型场景参数推荐表

这个表不是通用的，不同雷达参数下主瓣宽度、过采样率、信噪比不同，参数要按实测量做微调。跑参数扫描也不难：把窗口宽度系数做一组循环，比如0.1、0.15、0.2、0.25、0.3，每轮输出图像熵，选熵值最小的一档作为最终参数。

5. PGA自聚焦避坑指南：5个高频问题与排查流程

PGA在工程上不是一跑就灵的，以下列出的问题都是在实拍数据上翻过车之后总结出来的。每一条都按现象、原因、解决来写。

5.1 聚焦后图像反而变差

现象：运行PGA后图像熵上升，原本能看到的目标轮廓变得更模糊，甚至方位向出现双重成像。原因往往不在PGA本身，而是输入数据没有做过距离对齐。距离单元在慢时间上有走动时，同一个散射点的能量会分配到多个距离单元，PGA在逐个距离单元上做循环移位和加窗时，估计出来的相位梯度就不一致。

解决：在PGA之前先做距离对齐，常用做法是包络互相关法或最小熵距离对齐。某次处理一批舰船数据时，偷懒跳过了距离对齐，结果PGA跑完图像分裂成两个位置，补上对齐后再跑，一次就聚焦了。从那以后我每次处理ISAR数据都强制把距离对齐作为前置步骤，宁可多花几毫秒也绝不跳过。

5.2 相位估计出现卷绕锯齿

现象：估计出的相位校正序列phi_est在某个方位索引处突然跳变约2π，相位曲线出现锯齿状。原因是相邻方位单元的相位差超过π时，angle函数把差值卷绕到[-π,π]区间，这个跳变叠加到积分序列里就形成锯齿。

解决：在累加相位梯度之前，先对梯度做去卷绕处理，把相邻梯度差值限制在[-π,π]范围内：

grad_unwrap = grad; for n = 2:length(grad) delta = grad(n) - grad_unwrap(n-1); delta = delta - 2*pi*round(delta/(2*pi)); grad_unwrap(n) = grad_unwrap(n-1) + delta; end

这个操作在相位梯度和总校正序列上都要做一遍。如果你发现PGA.m里没有这段逻辑，建议自己补上，否则误差一大就会翻车。

5.3 迭代不收敛甚至发散

现象：运行前几轮迭代，图像熵在下降，但从第5轮或第6轮开始反向上升，梯度范数持续增大。原因一般有两个：一是加窗宽度太大，旁瓣干扰主导了梯度估计；二是更新步长太激进，直接减掉全量梯度导致过冲。

解决：先把窗口宽度系数往下调，比如从0.25调到0.15，再给梯度更新加一个松弛因子：

phi_est = phi_est - 0.6 * grad_unwrap;

松弛因子取0.6到0.8之间，付出一点收敛速度，换回稳定性。做完这两步，迭代发散的场景十有八九能救回来。

5.4 加窗宽度设太窄

现象：算法收敛了，图像也聚焦了，但强散射点的位置和真实位置对不齐，或者图像里有规则的波纹条纹。原因是窗口宽度低于主瓣3dB宽度，强散射点的部分主瓣能量被切掉，估计出的相位梯度出现恒定偏置。

解决：如果怀疑窗口太窄，最直接的检查方法是把加窗输出累加能量和未加窗时做对比，如果能量损失超过10%，就要放宽窗口。另外注意，循环移位和加窗的顺序不能反，先移位再开窗，保证窗口始终对准峰值中心。

5.5 整距离单元的循环移位失效

现象：图像中某个距离单元在多个方位位置都出现亮线，循环移位无法把能量集中到一个峰值。原因是目标转动非均匀或存在越分辨单元走动，同一个距离单元内可能包含多个多普勒频率成分，峰值位置不确定。

解决：这种情况不能直接对整距离单元做循环移位，需要把数据分成子孔径，分别估计相位梯度再做融合。子孔径长度通常取全孔径的四分之一到二分之一，重叠率50%左右。融合时对各个子孔径的phi_est做最小二乘拼接，务必在拼接前做去卷绕。

6. 验证技巧：用点目标仿真把PGA的效果“照”出来

自己改过PGA代码之后，第一件事不是拿实弹数据去试，而是用点目标仿真做冒烟测试。有误差和无误差的成像结果摆在一起，算法的好坏一目了然。仿真数据构造的关键是给回波注入一段已知的相位误差序列，比如正弦形式的相位误差，幅度0.1到0.5弧度、频率随机，这样PGA跑完和原始回波对比就知道估得准不准。

生成仿真的代码片段如下：

Na = 256; % 方位向脉冲数 t = (0:Na-1) / Na; phi_true = 0.3 * sin(2*pi*3*t) + 0.1 * randn(1, Na); sig = exp(1i * 2*pi * 20 * t) .* exp(1i * phi_true); img_raw = fftshift(ifft(sig), 2);

评估指标我用两个：归一化图像熵和对比度。图像熵越低说明聚焦越好，对比度则能反映主瓣的锐利程度。观察PGA补偿后这两项指标能否恢复到接近未加相位误差时的水平。部分实现里还会计算相位误差估计值与真实值的均方根误差作为算法精度的量化指标：

rmse_phase = sqrt(mean((angle(exp(1i*(phi_est - phi_true)))).^2));

用是因为相位误差存在周期性，直接拿phi_est减phi_true会算出很大的假误差，先转成单位复数再取角度才是物理上合理的差值。

验证在整条ISAR处理链路上的位置也要注意。PGA通常嵌在距离对齐完成之后、方位向成像FFT之前。实际处理流程建议：回波数据读入，距离压缩，距离对齐，PGA自聚焦，方位向FFT，幅度检波，最后做图像增强和显示。PGA输出的是相位校正序列，要乘在没有成像的慢时间回波上，而不是乘在已经做过方位FFT的图像上，这一点最容易搞错。跑仿真数据时，先让PGA对这种已知相位误差的估计误差在10%以内，再喂实测数据，这样每一步的翻车点都能定位到具体模块。

拿到这份源码包后的完整路线就是：先把PGA.m跑通，再用仿真数据验证，最后调参数处理实测数据。从那以后我每次拿到ISAR数据处理项目，都强制自己先把这个流程走一遍，确认每个环节到位了才敢继续往下做。希望帮到你。

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