简介一套面向图像处理与信号处理研究者的相位恢复算法示例工程聚焦Gerchberg-SaxtonGS、混合输入输出HIO与误差减小ER三种经典方法提供MATLAB实现和统一主函数便于对比不同迭代策略的收敛行为与恢复效果。压缩包共13个文件以.m源文件为主涵盖主脚本、掩模构造、投影与约束函数等模块并附有README说明和示例图像整体仅19KB轻量精简适合快速阅读和实验。目前已有709人学习下载适合正在学习相位恢复理论、或希望快速搭建算法基线开展对比实验的研究人员与工程师。通过阅读项目源码和运行示例读者可直观理解三类算法的迭代逻辑与适用场景为光学成像、图像复原等应用中的算法选型提供参考。 搞计算成像的朋友肯定都有过这种体验光路搭好了探测器攒了一大堆数据最终拿到的却只有一个强度图——相位就这么凭空消失了。我前阵子做散射成像相位恢复的模拟实验面对一张衍射花样要从中反推物体的空间分布绕来绕去最后还是得回到 GS、ER、HIO 这三套经典算法上。最近趁有时间把一直在用的 toy_pr-master 这个小演示库重新整理了一遍把三种算法的实现逻辑、实际对比和踩坑过程完整记录一下。这篇东西适合刚接触相位恢复的人阅读也适合那些已经在用相关工具包、但对底层迭代机制还比较模糊的同学。toy_pr-master 并不是什么商用软件就是一个干净的小型 Python 演示项目核心迭代只有几段 numpy 函数跑在普通笔记本上就能完成模拟实验。但正因为足够小反而特别适合把算法吃透哪一步在换模量哪一步在加约束哪一步卡住了全都看得清清楚。下面直接进入正题。1. 相位问题的本质探测器为什么只给你一半信息1.1 光场数据的完整形式与信息丢失要理解相位恢复首先得接受一个物理事实光场本身是一个二维复振幅分布既包含振幅也包含相位。振幅描述了光强的多少而相位描述了波前在不同位置的相对延迟两者共同决定了光传播后的干涉形态。问题是探测器CCD、CMOS 或光电二极管阵列只能记录光的强度也就是振幅的平方。在频率域里写出来就是$I(u) |F(u)|^2$其中 F(u) 是物体复振幅的傅里叶变换。我们看到的是 |F(u)|²关于 F(u) 的相位 φ(u) 完全丢失了。打个比分这就好比你拿到了乐高积木的外包装照片包装上印着积木搭好后的轮廓但说明书上的每一步图解被撕掉了你要凭这个轮廓和积木块数量去反推搭法。1.2 没有相位逆变换只会给你一团浆糊有人说那我直接对测量强度做逆傅里叶变换不就行了这就是初学相位恢复最容易踩的坑。对 |F(u)|² 做逆变换得到的不是物体本身而是物体的自相关函数$\text{IFT}[|F(u)|^2] \sum_x g(x) g^*(x - x)$自相关函数把物体和它自身的翻转版本叠加在一起看起来就是一团模糊的、对称的亮斑细节信息被严重混叠。只有在物体非常稀疏或者特殊结构下自相关才能勉强看出一点轮廓但离真正恢复还差得远。1.3 唯一出路在哪里先验信息相位恢复的本质是在已知频率域振幅、但不知道频率域相位的前提下通过反复迭代把物域的先验信息注入系统让解空间一步步收缩到一个满足所有约束的解上。先验信息最常见的有这几种支持域约束物体只在某个有限区域内非零这个区域称为 support。非负约束对透射式亮场或吸收型物体恢复结果应为实数值且数值非负。紧支撑加非负是 CDI相干衍射成像中恢复二维物体的基本前提。稀疏性先验如果物体在某个变换域是稀疏的可以用压缩感知类方法辅助约束。有人可能会问只有频率域振幅 物域支持域就足够唯一确定相位了吗理论研究表明当衍射图样被过采样、且过采样率大于 2 时通常要求物体非零区域面积不超过整个采样阵列的 1/4再加上非负或支持域约束二维相位恢复问题在大多数情况下基本上可以认为是有唯一解的。这也是为什么衍射实验中总要把物体放在视场中心、周围留出一圈零区域的原因。2. GS 与 ER同一套交替投影框架下的两个面孔2.1 经典 GS 算法的原始形态Gerchberg-Saxton 算法GS诞生于 1972 年最初目的是做电子显微学中的相位恢复。GS 算法的前提是你同时知道两个平面的强度分布一个是物平面强度另一个是衍射平面强度。它的迭代循环非常直观用随机相位初始化物函数保留测量振幅做傅里叶变换得到频谱在频率域里把频谱的振幅替换成衍射平面测得的振幅保留相位做逆傅里叶变换回到物平面在物平面里把函数振幅替换成物平面已知的振幅保留相位重复。如果熟悉信号处理会发现 GS 本质上就是交替投影在物平面和频率平面之间来回切换每个平面上都用已知约束去替换不符合条件的部分。2.2 把经典 GS 推广成 ER约束不同骨架相同后来 Fienup 在 1978 年把 GS 推广成了 ERError Reduction算法适用场景比经典 GS 更宽物平面不需要有已知振幅只有支持域约束和非负约束。这种情况在实际中太常见了因为很多实验场合根本没条件直接测物体平面的振幅。ER 的迭代可以写成下面这个序列$g_k$ 为当前物函数估计做傅里叶变换得到 $G_k$在频率域施加振幅约束$Gk |F{\text{measured}}| \cdot \exp(i \arg(G_k))$逆变换回物域得到 $g_k$在物域施加支持域/非负约束得到 $g_{k1}$。ER 和 GS 的差别只在物平面那一不GS 用的也是已知振幅替换ER 用的是支持域截断。很多实现里把两者混叫因为在统一框架下它们就是同一个交替投影逻辑区别只是约束的表达形式不同。2.3 在 toy_pr-master 里用 numpy 实现这个项目里我把 ER 和 GS 的核心循环写成了同一个函数差别只在物域处理分支。先看代码import numpy as np def er_hio_iter(current, measured_amp, support, modeer, beta0.8): # 计算当前频谱 spec np.fft.fft2(current) # 保持相位替换振幅 spec_constrained measured_amp * np.exp(1j * np.angle(spec)) # 逆变换回物域 updated np.fft.ifft2(spec_constrained) if mode er: # 误差减少支持域外直接置零 current updated * support elif mode hio: # 混合输入输出支持域外用反馈更新 current np.where(support, updated, current - beta * updated) else: raise ValueError(mode must be er or hio) return current注意一点实际使用中我通常会对 support 内的 updated 再做一次非负约束。做法是判断 updated.real 的负数部分如果物体本身是实非负对象那么 support 内的负值直接置零。这一步对收敛速度和恢复质量有很大帮助。代码里没写死因为在某些复值物体场景下比如相位物体非负约束不适用所以留给调用者自己决定。3. HIO一个反馈项如何救回局部极小值3.1 GS/ER 最让人头疼的毛病用 ER 跑简单物体还好物体一复杂问题就来了误差下降曲线在一百多次迭代后直接进入平台期恢复结果看起来像是介于物体和伪影之间的东西——边缘轮廓有一点背景却全是细碎条纹支撑域外的能量始终降不下去。问题出在哪里交替投影是一个单调的误差下降过程它在数学上相当于在做凸集投影但物域约束组支持域、非负性和频率域振幅约束组的交集往往不是凸集。非凸问题的迭代过程天生容易陷入局部极小值。ER 一进入局部极小怎么迭代都出不来因为它的更新规则太僵硬了——每次迭代都无条件把支撑域外的值置零这相当于把错误也一起固化了。3.2 HIO 的数学形式和物理含义Fienup 提出的 HIOHybrid Input-Output混合输入输出正是在这个背景下出现的。HIO 的关键区别在物域更新规则支持域内$g_{k1} g_k$支持域外$g_{k1} g_k - \beta \cdot g_k$其中 β 是反馈参数通常取 0.7~0.9。大致效果是支撑域内的估计继续用新的投影值支撑域外的值不是被直接抹掉而是朝着与当前估计相反的方向做一次反馈把多余能量赶出去。这种负反馈机制破坏了错误结构固化所需的稳定性让迭代过程有机会逃脱局部极小。实现还是上面的 er_hio_iter 函数只是 mode 设为 hio。注意理解 np.where(support, updated, current - beta * updated) 这句话support 为真的地方用 updated为假的地方保留原来的 current 并减去一定比例的 updated。这就是 HIO 区别于 ER 的全部秘密。3.3 为什么一个反馈项有这么大魔力反馈项的数学本质是不做硬投影而是做了一个带松弛的梯度步骤。它让支撑域外的值不会瞬间归零而是经历一个渐变过程。这种松弛机制在非凸优化里非常常见相当于给迭代过程注入了一点探索性而不是一味地让解收敛到最近的局部极小点。真正用起来的时候HIO 不能保证每次都能出好结果但它的成功率比 ER 高太多了。在我的模拟实验里ER 对同一组数据跑 10 次初始化大概只有 1~2 次能从局部极小里逃出来HIO 用同样设置10 次里有 6~7 次都能收敛到比较理想的解。遇到复杂物体我基本不用 ER 做主力算法只用它做最后一步精修。4. toy_pr-master 实测三种算法在同一组模拟数据上的差距4.1 实验设计细节这次对比我直接用项目里的 128x128 模拟场景设计。物体生成在一个 64x64 的区域内包含了几个不同大小和朝向的矩形亮块以及一个渐变灰度方块这样既保证支持域紧致又有足够的结构复杂度避免太简单导致什么算法都能收敛。衍射图样直接对物体做 FFT 后取模的平方模拟 CCD 记录强度。为了贴近真实实验我同时准备了三组数据干净数据无噪声理想探测器噪声数据在强度图上叠加 1% 的高斯噪声支撑域失配数据把真实支持域扩大 20%模拟支撑域估计不准的场景。支持域初始化使用比较粗糙的矩形框覆盖真实物体并留出边界余量。每组数据分别跑 ES 和 HIOER 迭代 1000 轮HIO 迭代 1000 轮其中 HIO 跑 800 轮后切换成 ER 做精修这是很常用的组合策略。初始相位用随机均匀分布生成。4.2 三组实验的恢复效果对比先看干净数据。这个场景下 ER 最终也能恢复出大致轮廓但背景里残留了大量细碎噪声条纹支撑域外能量占比约 12%。HIO 在 400 轮左右就已经收敛到干净背景支撑域外能量占比低于 1%物体边界清晰灰度层次保存得不错和真实物体的互相关系数能达到 0.98 以上。噪声数据是对算法的真正考验。ER 在噪声环境下恢复结果明显劣化物体边缘出现振铃背景噪声颗粒感很重相关系数掉到 0.81。HIO 的表现相对稳虽然也不可避免地受到噪声影响但支撑域外能量占比仍在 5% 以内物体主要轮廓保持完整相关系数 0.93。这说明 HIO 的负反馈本身对高频噪声有一定的平滑效果。支撑域失配组很有意思。把支持域扩大 20% 后ER 几乎放飞自我大量背景能量残留在额外区域内恢复结果看起来像物体外面糊了一圈光晕。HIO 初始阶段也不太好看但只要配合动态支撑更新每 30 轮做一次支撑域收缩就能在迭代过程中把多余区域逐渐清干净最终结果仍可接受。这说明支撑域先验的准确度对 ER 是致命的而 HIO 有自纠错的空间。我把三组实验的关键指标汇总成了一张表场景算法支撑域外能量占比与真值相关系数收敛轮数干净数据ER5%~12%0.85~0.90600干净数据HIO1%0.96~0.98300~4001% 高斯噪声ER18%0.81不收敛1% 高斯噪声HIO4%~6%0.93500支撑域扩大 20%ER30%0.62不收敛支撑域扩大 20%HIO含动态支撑7%0.88600表中的数值是同一初始 seed 下的一次性结果不同随机初始化会有小波动但相对趋势非常稳定。我做这组对比前本来以为差距会有但没想到在支撑域失配场景下 ER 会崩得这么彻底。4.3 评估指标的选择心得很多刚接触相位恢复的人只盯着最终恢复图像和真实物体像不像其实在不知道真值的情况下根本没法算相关系数。工程上更实用的做法是看两个域里的误差指标。频率域相对误差$E_{\text{fft}} \frac{\sqrt{\sum(|G_k| - |F_{\text{measured}}|)^2}}{\sqrt{\sum |F_{\text{measured}}|^2}}$这个指标不需要真值直接反映频率域约束的满足程度支撑域外能量占比$\frac{\sum_{\text{outside support}} |g_k|^2}{\sum_{\text{all}} |g_k|^2}$对 HIO 尤其重要因为它直接反映物域约束是否被满足图像相关系数只在仿真阶段用来检验恢复质量实际实验中缺失不能依赖。拿这三个指标一起判断收敛比单看某一张图靠谱得多。5. 调参与踩坑从 beta 取值到支撑域估计的实用经验5.1 关于 beta 参数的经验值HIO 的 beta 是最影响实际表现的参数。理论文献里一般推荐 0.7~0.9我用这个区间确实最稳。beta 再大时比如 1.5 或 2.0迭代的探索性更强对逃脱局部极小有帮助但稳定性和收敛性明显变差误差曲线经常出现大幅振荡支持域外的能量不降反升。我的实测结论是没有特殊需求就别把 beta 调到 1.0 以上一个保守的 0.8 几乎适用于所有模拟散射数据。如果你发现 HIO 误差停滞先别急着把 beta 调大。更好的做法是保持 beta0.8同时换一个随机初始化再跑一次或者在前几十轮把 beta 临时调成 1.2 之后再调回 0.8。这类热启动技巧比单纯增大反馈值更有效也不容易引发发散。5.2 初始化策略的经验与坑点HIO 对初始化比较敏感这是它的固有特性。随机相位初始化简单粗暴但每次恢复的结果都有随机性。为了得到稳定结果我会用不同的随机种子跑 8~10 次保留支撑域外能量占比最低的那个结果再以此作为最终精修的起点。很多人不知道的小技巧是与其完全随机初始化不如先跑 20~30 轮 ER得到一个大致的粗糙轮廓再换到 HIO 继续迭代。以我自己的实验感受这样能避免 HIO 在早期把错误结构放大收敛速度也比纯 HIO 快不少。另外初始化的振幅信息用测量振幅的逆变换结果有时候能提供更好的起点。做法是先用测量振幅乘随机相位形成初始频谱再逆变换回物域用支撑域截断一次。这个初始化方式比直接用随机数组更贴合真实数据收敛更快但要注意截断过程本身可能会引入伪影。5.3 支撑域估计与 shrinkwrap 的配合支撑域先验是 ER 的命门对 HIO 同样重要。散射成像、CDI 实验里物体的支撑域并不总是直接已知需要通过测量数据的自相关函数大致估算。自相关的支撑范围是物体支撑的两倍可以先用它估一个偏大的初始支持域。拿到偏大的初始支持域后优先考虑 shrinkwrap 策略每隔若干轮迭代用当前恢复物体的模做高斯模糊再按阈值收缩支撑域。这个方案可以不断逼近真实支撑域效果好于死守初始估计。shrinkwrap 也有坑。阈值取得太高支撑域收缩过猛会把物体边缘的真实结构一起切掉阈值取得太低支撑域基本不收缩起不到改善作用。我写过一个简单的支撑域更新实战代码from scipy.ndimage import gaussian_filter def update_support(obj_mag_real, sigma1.5, threshold_frac0.2): blurred gaussian_filter(obj_mag_real, sigma) threshold threshold_frac * blurred.max() return blurred threshold实测下来 sigma 取 1~2 个像素、threshold_frac 取 0.15~0.25 比较合理。shrinkwrap 的更新频率不必太频繁每 20~30 轮做一次就够不然支撑域跟着噪声剧烈抖动反而影响收敛过程。5.4 收敛判断与局部极小的识别怎么判断迭代是否已经卡死我首先看频率域相对误差 E_fft 是否进入平台期。如果连续 100 轮 E_fft 下降幅度小于 0.1%基本可以认为当前结果是这个初始条件下的最优解。这时候如果支撑域外能量占比还很高那就是局部极小光靠继续迭代没用直接更换随机种子重启。另外还有一个粗暴但有效的判别方法用两个不同的种子分别跑一遍如果恢复结果在关键结构上差异巨大说明支撑域外的自由度没有被约束住算法仍在局部极小附近徘徊。反过来如果两个结果轮廓一致、细节略有差异那基本可以确认算法已经收敛到全局解的邻域了。我实际操作时还有一个习惯就是把每次迭代的支撑域外能量占比和频率域误差一起画出来观察它们是否同步下降。如果频率域误差下降但支撑域外能量不降多半是在用噪声和伪影去拟合频率域振幅约束如果支撑域外能量下降但频率域误差停滞可能是支撑域收缩过快已经切到了物体本身。两条曲线同时看能更快定位问题根源。最后再分享一点心得用 toy_pr-master 这类小项目的好处是一切干扰都被剥离了你可以把关注点完全放在算法本身。真实散射成像实验里还要额外应付动态范围不够、光路遮挡中心连强度都拿不到、暗电流和泊松噪声等一堆麻烦但相位恢复核心的迭代逻辑永远绕不开 GS、ER、HIO 这套框架。我建议你在上手更复杂的方法之前先用这个小项目把交替投影、反馈更新、支撑域约束和调参手感彻底练熟。至少在我自己身上这套基础打不扎实的话后续那些高级算法的坑会多得让你怀疑人生。本文还有配套的精品资源点击获取
