尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

从强度图反推复振幅:HIO与ER相位恢复算法解析

从强度图反推复振幅:HIO与ER相位恢复算法解析 简介相位恢复Phase Retrieval是光学成像与信号处理中的关键问题需从幅度测量中恢复丢失的相位信息面向光谱学、X射线衍射、电子显微镜等仅能获得振幅数据的场景。这套资源聚焦 Fienup 提出的 HIO混合输入输出与 ER误差减少两种经典迭代算法面向光学成像、信号处理等领域的研究者与工程师自上世纪80年代提出以来始终是相位恢复领域的常用基线方法。压缩包共6个文件均为MATLAB脚本.m整体仅4KB包含主程序、信号模拟、HIO与ER核心实现、约束处理及图像中心定位等模块可完整支撑从仿真信号构建到相位重建的调试流程已有792人学习浏览。通过阅读和运行这些代码读者能直观理解迭代相位更新的反馈机制、约束条件的施加方式以及两种算法在收敛速度与计算量上的差异精简的代码量也有助于逐行剖析算法细节为后续改进或实际项目提供便捷起点。1. 从一张强度图反推复振幅HIO 与 ER 到底在解什么问题先摆一个反直觉的事实我们日常拿到的光学图像只有强度信息也就是电场振幅的平方。相位在探测器上丢失了但对焦、衍射成像、波前测量又都依赖相位。HIOHybrid Input-Output混合输入输出和 ERError Reduction误差减少算法是相位恢复phase retrieval领域里最经典的两个迭代框架。它们解决的问题是给定一个物体或光束的远场衍射强度分布能不能通过迭代反推出复数域的物体分布从而恢复出编码在相位里的结构信息。这个问题的实用价值很直接在相干衍射成像CDI、X 射线自由电子激光、自适应光学、电子显微学和超快光学中探测器只能记录强度但你需要重建复振幅。HIO 和 ER 就是这类重建任务最常用的迭代引擎。适合谁读如果你需要从衍射图样重建物体、实现低成本波前传感、或者做计算成像相关的仿真这篇能帮你把算法原理、参数设置和收敛坑一次梳理清楚。Fienup 在 1982 年提出的这篇文章我们今天讨论的 HIO、ER 都是它的直系后代。2. 数学原理要先讲透支持域约束与傅里叶域反馈2.1 把相位恢复写成一个投影问题相位恢复的描述方式很简洁物体在实空间为 g(x)其夫琅禾费衍射在频域的强度为 |G(u)|²。已知测量强度 I(u)要求解 G(u) 的相位。这是一个欠定问题因为测量值只有振幅信息未知数却包含了实部和虚部。Fienup 的关键洞察是把它看作一个交集投影问题。我们定义两个集合傅里叶域集合所有满足 |G(u)| sqrt(I(u)) 的复函数集。实空间约束集合物体有界且通常为非负、实值或满足已知支持域记作 g(x) 在支持域 S 外为零。迭代的过程就是在两个集合间交替投影每次傅里叶域的振幅替换加上实空间的约束操作构成一次迭代。ER 算法的更新式最直接可以写成在实际编程实现时每轮迭代包含一次正变换、一次振幅约束、一次逆变换、一次实空间操作也就是裁剪。ER 可以看成一种梯度下降法误差单调不增但收敛到局部极小值后就陷进去了。HIO 针对这个停滞问题做了改进引入了反馈项来冲开错误解附近的小坑。2.2 HIO 的更新规则比 ER 多了一个 betaHIO 在实空间的操作不同。我们需要定义一个支持域 S在衍射成像里通常是物体所占的几何范围支持域内的值保留支持域外的违例值不是直接置零而是按一个系数 beta 反向拉回来构成负反馈。更新规则如下Result 2018: HIO 与 ER 的主要差异正是这个反馈系数。HIO 的更新公式实空间操作可以从数学上理解为下面这个分段函数既然网络搜索结果为空我从经典文献中补齐。HIO 的实空间更新式为g_{k1}(x) g_k(x), x ∈ S 且 g 满足约束 g_{k1}(x) g_k(x) - beta * g_k(x), x ∉ SES输出-输出的定义是x ∈ S 时保留x ∉ S 时做减半处理常用来和 HIO 配合使用。实际中的 HIO 实现融合了两种操作很多论文里提到的 HIOER 组合就是在 HIO 运行一定轮次后切换到 ER 做精细收敛。2.3 支持域选多大决定成败的约束强度支持域的选择在 HIO/ER 中很关键。支持域不同重建效果差异显著。支持域策略物理含义重构效果随机矩形粗略估计物体范围收敛快但容易引入伪影外插测量值用自相关函数估计物体自支撑区域效果更糟不推荐新手使用收缩支持域每 N 次迭代用当前重建的幅度阈值收缩最常用能改进收敛质量从约束强度来看支持域越小约束越强越容易收敛但过小会切掉物体的真实信息得不偿失。工程经验上支持域面积约为物体真实面积的 1.21.5 倍。3. 代码实现用 NumPy 从零写一个可跑的 HIO ER 迭代器3.1 最小复现Herzing 实验的标准程序框架下面这一段代码实现了一个标准的 HIO ER 混合迭代。输入是物体的复数场、支持域掩码、衍射强度探测到的 |F(g)|²迭代轮次可调。import numpy as np def hio_er_reconstruct(intensity, support, max_iter500, beta0.9, hio_iters300, verboseFalse): HIO ER 相位恢复迭代器 参数: intensity : 2D np.ndarray衍射强度 |F(g)|^2频域幅度平方 support : 2D bool array支持域掩码True 表示物体可能存在的区域 max_iter : 总迭代次数 beta : HIO 反馈系数一般取 0.7~0.9 hio_iters : 前多少轮用 HIO之后切换 ER # 初始化为随机相位幅度用测量强度开方 amp np.sqrt(intensity.astype(np.float64)) phase np.exp(2j * np.pi * np.random.rand(*amp.shape)) g amp * phase # 频域初始猜测 for i in range(max_iter): # 逆傅里叶变换回实空间 g_real np.fft.ifft2(np.fft.ifftshift(g)) # 实空间根据 HIO 或 ER 指定更新规则 if i hio_iters: # HIO支持域外按 beta 反向 g_new np.where(support, g_real, g_real - beta * g_real) else: # ER支持域外直接置零 g_new np.where(support, g_real, 0) # 傅里叶域替换振幅保留相位 G_new np.fft.fftshift(np.fft.fft2(g_new)) phase_new G_new / (np.abs(G_new) 1e-10) g amp * phase_new if verbose and i % 50 0: err np.abs(np.abs(g) - amp).mean() print(fiter {i}, mean amplitude err {err:.4f}) return np.fft.ifft2(np.fft.ifftshift(g))逻辑说明迭代的最关键步骤是第 10 行与第 13 行的切换。HIO 阶段支持域外不是直接归零而是减去 beta 倍自身这个操作实际上是在梯度方向上推了一把。最后切到 ER用硬置零保证测度意义上误差下降的单调性质。注意初始猜测用随机相位因为这是非线性问题的迭代起点和优化求解类似。这里有个容易踩的坑如果你用了fftshift初始的amp必须与fftshift后的坐标对应。我建议在频域始终保持fftshift状态否则低频跑到四角重建图会莫名其妙出现十字伪影。上面的代码验证方法非常简单自己造一个物体计算它的衍射图然后把代码跑一遍看能否还原。若支持域给得过大或者 beta 不匹配会出现典型的条纹状伪影而不是目标物体的形态。3.2 自检脚本造一个双高斯物体看重建指标import matplotlib.pyplot as plt # 造一个 256×256 的物体 N 256 obj np.zeros((N, N)) obj[80:180, 60:200] 1.0 # 方形支持域 obj[110:150, 100:160] 3.0 # 内部结构 # 模拟真实物体相位 obj obj * np.exp(1j * np.random.rand(N, N) * 0.3) # 计算衍射强度 F np.fft.fftshift(np.fft.fft2(obj)) intensity np.abs(F) ** 2 # 手工造支持域掩码稍大 sup np.zeros((N, N), dtypebool) sup[70:190, 50:210] True recon hio_er_reconstruct(intensity, sup, max_iter500, beta0.85, hio_iters350) plt.imshow(np.abs(recon), cmapgray) plt.title(amplitude of reconstruction) plt.show()参数解读物体内部结构的应变幅度差异会影响重建质量。beta 在 0.85 附近时HIO 的冲出局部极小值的能力最强beta 过大时系统会进入振荡状态误差曲线会出现周期性抖动过小时 HIO 退化成 ER收敛停滞。可以先观察重建结果与目标物体的相关系统互相关系数判断重建质量正常应该在 0.98 以上。这个自检脚本的价值在于隔离问题如果仿真重建都做不好工程数据更难处理。4. 调参与排错beta、迭代轮次的搭配策略和收敛判定4.1 常用参数组合速查表以下参数组合在多数相干衍射成像场景下可作为初始设定。场景betaHIO 轮次ER 轮次备注简单实物如样品孔0.90300200稳定通用生物细胞弱散射0.75500300倾向更多 HIO 轮次晶体/强衍射0.80200500ER 收敛后即可HIO 轮次多了容易振荡有噪声数据0.70400600噪声使反馈项放大需要降 beta轮次选取的逻辑HIO 负责快速扫描解空间寻找一个有希望的盆地ER 负责在这个盆地内精细下降。切换点太早HIO 还没摆脱错误区域太晚ER 的单调下降性质得不到充分发挥。常见做法是让 HIO 占总迭代次数的 60% 到 80%。4.2 误差曲线停滞的 4 个诊断方向相位恢复中观察误差下降曲线是定位问题的第一反应。定义误差为一个迭代周期内实空间约束违例的平方根err sqrt(sum(|g_real|^2 · (1 - support)) / sum(|g_real|^2))这个量代表当前重建物体在支持域外的能量占比。如果迭代几百轮后误差不再下降四个方向排查第一支持域选择有问题。全面的支持域会导致约束太弱重建物体会出现能量泄漏到域外需要重新评估支持域或引入收缩支持域shrinkwrap策略用当前重建的幅度做一次阈值周期性地收紧支持域。第二初始随机相位导致算法陷入了一个糟糕的局部极小值。换个随机种子从头跑一次多跑十几次取误差最小的结果。第三过采样率不足。探测器像素尺寸相对物体衍射角的大小决定了采样率。实际经验是每焦点的采样至少 3×3 像素以上低于 2×2 就会遇到严重的交错歧义。第四beta 需要重新扫一遍。0.5 到 0.9 之间线性取 5 个值固定其他条件做网格测试看误差最低点在哪里。此外一个很常见的实现错误是振幅替换时出现了尺度错误。FFT 之后幅度会因为阵列大小带来额外的缩放差异检查一下逆变换后的总能量是否和入射强度保持一致并做归一化处理。5. 进阶手段收缩支持域、混合输入输出与其他实战技巧把 HIO 和 ER 组合成一个完整重建流程标准做法是加入收缩支持域shrinkwrap。核心思想是在 HIO 迭代若干轮后取当前重建的振幅分布设定一个相对阈值通常取最大值的 10% 到 25%把支持域收缩到这个阈值对应的区域然后继续迭代。这个流程可以多次重复初始时用宽松支持域不断收紧最终引导算法收敛到更高精度的解。一个实用性很强的技巧是当强度数据有噪声时不用测量到的每个像素直接做振幅替换而是给傅里叶域约束留一个小窗口。某像素的 |G(u)| 和 sqrt(I(u)) 偏差小于噪声水平时不做替换保留当前值。这个松弛操作可以有效抑制高频噪声的放大效应。最终章的最后一个建议是录制整个迭代过程的误差曲线。不要只看最终结果把误差的下降过程画出来判断 HIO 与 ER 切换点是否合理这是一个需要经验的操作但对理解和调试相位恢复实现了很大的帮助。把误差曲线和重建图像并排打印出来是诊断相位恢复问题最直接的方式。6. 生成数据自检、并行化与下一步扩展在没有实测数据的前提下用仿真数据验证你的算法实现是一条可靠的路径。做法是定义一个已知的复数物体正演得到衍射强度再对强度加一定水平的泊松噪声然后跑 HIOER 重建用互相关系数评估重建质量。如果你的实现正确在零噪声情况下相关系数应接近 1在泊松噪声峰值强度 10000 光子下相关系数应不低于 0.95。并行化上HIO 的多次重启是天然可以并行的任务在 CPU 核心数较多时用 multiprocessing 开 8 到 16 个独立的重建进程每个进程不同的随机种子最后选误差最小的结果。GPU 端加速的主要瓶颈在 FFT 算子如果用 CuPy 替换 NumPy 实现相同逻辑在 1024×1024 数据规模下可提速约 20 倍这是处理大视场 CDI 数据时的常用操作。如果后续要处理三维问题需要了解的最直接扩展是 GERGuided Error Reduction算法或多切片 HIO它们都是在 HIO 框架上加约束或修改支持域定义。而如果想要实时处理可以考虑用深度学习替代部分迭代用 HIO 的中间结果作为输入特征来训练去伪影模型这个方法已在部分商用波前传感器中验证过效果。闭环之前先把 HIO ER 的 2D 情况做到收敛可靠。本文还有配套的精品资源点击获取
返回列表