简介相位恢复是图像与信号处理中的关键问题。这份MATLAB源码包聚焦三种经典迭代算法Gerchberg-SaxtonGS、混合输入输出HIO与误差减小ER并集成振幅约束、支撑域更新及误差反馈等模块帮助研究者在统一框架下对比各算法的收敛速度与恢复精度。资源面向光学成像、X射线衍射、通信信号处理等领域的工程师和学生通过单一主程序即可调用不同算法调整掩膜与支撑参数来适应不同退化场景。压缩包共13个文件以9个.m脚本为主涵盖主函数、投影算子、掩码生成、误差统计及绘图工具每个模块职责清晰另有说明文档与示例图像整体仅19KB结构简洁轻量便于直接阅读与二次开发。目前已有709人学习适合作为相位恢复算法教学演示或预研验证的入门参考。借助这套代码使用者可以快速建立迭代重建流程直观比较HIO、ER与GS的差异为算法选型和后续优化提供实验依据。1. 相位恢复不是玄学toy_pr-master 里的 GS、ER、HIO 到底能做什么拿到toy_pr-master_HIO_ER_GS_相位恢复_这个项目名时第一反应是先拆词toy 表示这是一个轻量级、用来验证思路的 toy 级实现PR 是 phase retrieval相位恢复后面三个缩写 HIO、ER、GS 正是这个领域最经典的三种迭代算法。相位恢复要解决的问题非常具体探测器只记录光的强度相位信息在记录瞬间就丢了而散射成像相位恢复这类场景里你手里往往只有一张衍射光斑强度图却要逆向重建出样品的复振幅分布。这个项目能帮你把三种算法的行为差异、参数边界一次看明白适合刚接触相干衍射成像或散射成像重建的工程师和学生也适合要快速跑 baseline 的熟手。2. 从 GS 到 HIO三种迭代相位恢复算法的递进关系动手跑代码之前先把三种算法的关系理清楚。很多人在 toy_pr 上翻车不是因为代码写错而是把 GS、ER、HIO 当作三个独立算法去背没意识到它们其实是同一套傅里叶域替换幅度、空域做投影框架下的三个变体。理解了递进关系后面调参才不是玄学。2.1 只测强度也能找回相位反问题为什么能解探测器记录的是强度图物理上可以写成 I(u)|FT{o(x)}|²也就是样品复振幅 o(x) 做傅里叶变换后取模平方。取模这个动作把相位信息扔掉了。直觉上少了一半信息为什么还能逆回来原因有两条。第一傅里叶变换是全局算子空域里任何一点的变化都会牵动频域所有点的取值强度图里实际嵌入了大量关于物体结构的统计信息只是不以直观形式呈现。第二也是更关键的重建时永远有先验约束。比如你知道物体只存在于一个有限区域或者物体是实值、非负的这些约束会把解从无穷多个压到很小的范围。相干衍射成像和散射成像相位恢复做的就是这件事用支持域约束配合过采样把丢失的相位信息补回来。原理上这不是猜而是信息量充足前提下的反问题求解只是迭代过程不保证每次都收敛到你想要的那个解。2.2 GS 算法双平面强度之间的乒乓迭代GS 算法是 Gerchberg 和 Saxton 在 1972 年提出的有资料也追溯到 1960 年代的工作它最朴素也最好理解。适用场景是你知道两个平面上的强度分布比如样品面的已知振幅 A₁ 和衍射面的已知振幅 A₂求相位分布。迭代过程就像一个乒乓球在两个平面之间来回弹用已知振幅 A₁ 和随机初相位组合成初始场 u₀对 u₀ 做傅里叶变换得到 G用衍射面已知振幅 A₂ 替换 |G|保留相位反变换回空域再用 A₁ 替换空域振幅保留相位重复第 2、3 步直到相位收敛。这个流程里的每一步都是往一个约束集合上做投影频域约束是振幅必须等于测量值空域约束是振幅必须等于已知值。两个集合都是凸集的话交替投影保证收敛。实际中 GS 收敛很快几轮就能看到结构。它的最大限制是必须有第二个面的强度如果只有单张衍射强度图空域约束就不成立了。2.3 ER 算法把第二个强度约束换成支持域回到最常见的实验场景探测器只有一张衍射强度图空域里没有已知振幅可用。这时 Fienup 在 1978 年提出的 Error ReductionER算法把空域约束替换成支持域物体只可能存在于一个有限支撑区域内区域外完全是空的。每一轮迭代是对当前估计 g_k 做傅里叶变换在频域用实测振幅替换G_k sqrt(I_measured) * exp(i·arg G_k)相位保留反变换回空域得到 g_k空域更新支持域内直接取 g支持域外强制置零。ER 的名字很直白每次更新后误差单调不增一定会往下降。问题在于支持域约束的集合不是凸的交替投影容易卡进局部极小。典型现象是误差曲线前 30 轮降得很快之后进入平台期跑再多轮也不动。很多人在这一步以为迭代次数不够把 ER 跑到 5000 轮误差曲线还是平的——这不是程序员的问题是算法的几何特性决定的。2.4 HIO 算法用负反馈跳出停滞区HIOHybrid Input-Output是 Fienup 在同一篇工作里提出的修正方案和 ER 共享频域幅度替换步骤差别只在空域更新规则。HIO 的空域更新不再把支持域外的值硬置零而是保留上次估计并加入一个负反馈项支持域内g_{k1} g_k支持域外g_{k1} g_k - β·g_k这里的 β 是反馈系数典型取值范围 0.5 到 0.9。大白话解释这个更新规则支持域外的值不清零而是把溢出物放大 β 倍后反号叠加上一轮估计相当于往迭代里注入一个探测性扰动。这样下一次傅里叶变换时历史信息和反变量一起参与约束不容易在同一个局部极小点死磕。这也是 HIO 能成为相干衍射成像领域主力算法的原因。需要提醒一个容易混淆的点β1 时 HIO 并不等于 ER因为 ER 是支持域外直接置零而 HIO 在 β1 时支持域外是 g_k - g_k两者数学形式不一样。toy 级项目里经常把 ER 写成支持域外置零的掩膜操作把 HIO 写成一行带 feedback 的更新两者代码只差一行但收敛行为差很远。算法空域约束空域更新规则β 参数典型场景GS已知空域振幅替换振幅、保留相位无光束整形、全息相位计算ER支持域支持域外置零无CDI 基础验证易停滞HIO支持域 反馈支持域外 g-βg0.5~0.9相干衍射成像、散射重建主力3. 跑通最小流程模拟数据自检与 HIO 核心循环现在进入正题。toy_pr-master 这类项目最方便的地方是把三套算法骨架放在一个极小代码量里你一眼能看出 ER 和 HIO 的差别只在空域更新那一行。我建议拿到手先不要碰真实数据用模拟数据把接口跑通、把误差曲线看熟再上实验数据。否则真实数据里的噪声、坏点、饱和会让你以为是算法写错了。3.1 数据流强度图怎么变成算法可吃的输入标准数据流分四步。第一步读取衍射强度图 I它是一个二维实数数组第二步取振幅 Asqrt(max(I, 0))因为算法替换的是振幅而不是强度第三步对振幅做归一化比如除以最大值避免数值量级影响傅里叶变换的浮点精度第四步构造初始复数场 u₀支持域内是随机相位乘以常数振幅支持域外置零。有的 toy 实现会略过归一化模拟数据没问题实测数据不归一化会导致后期迭代误差计算溢出或精度丢失。另外多数 toy 项目默认数据已经做了暗场减除实测数据如果没减背景重建结果会带一圈低频晕染这点先记住第 5 章还会展开。3.2 生成模拟衍射数据先给算法一份标准答案自己造模拟数据相当于带着标准答案考试。下面这个函数生成一个圆形支持域内的复数物体并返回它的远场衍射强度足够用来验证后面 HIO 循环。import numpy as np def make_sample(n256, support_radius24, seed1): rng np.random.default_rng(seed) x np.linspace(-1, 1, n) X, Y np.meshgrid(x, x) # 圆形支持域半径按相对坐标计算 support (X**2 Y**2) (support_radius / n * 2) ** 2 obj np.zeros((n, n), dtypecomplex) # 支持域内随机振幅 随机相位模拟复值样品 obj[support] rng.normal(0, 1, support.sum()) * np.exp( 1j * rng.uniform(0, 2 * np.pi, support.sum()) ) far_field np.fft.fft2(obj) amp np.abs(far_field) intensity amp**2 / amp.size # 归一化强度 return obj, support, intensity逻辑说明support_radius / n * 2把像素半径换算成 -1 到 1 的坐标比例圆形 mask 比矩形更接近真实样品约束。物体是复值随机场模拟散射样品的振幅和相位都未知的情况。fft2后取平方再除以像素数是为了让强度量级不随图像尺寸膨胀。参数上需要注意support_radius和n的比例决定了过采样比示例里 n 是 256、半径 24每个维度上支撑占比约 0.19满足过采样要求。如果想模拟噪声可以对 intensity 做一次泊松采样。3.3 HIO 主循环与误差记账核心 HIO 迭代可以压缩成一个很短的函数这也是 toy 级项目最迷人的地方。def hio_recon(intensity, support, beta0.8, n_iter300): amp np.sqrt(intensity) # 从强度恢复振幅 rng np.random.default_rng(0) u np.zeros_like(support, dtypecomplex) u[support] np.exp(1j * rng.uniform(0, 2*np.pi, support.sum())) err_curve [] for i in range(n_iter): G np.fft.fft2(u) G_new amp * np.exp(1j * np.angle(G)) # 频域幅度替换 u_new np.fft.ifft2(G_new) # 反变换回空域 # HIO 更新支持域内接受新值支持域外做负反馈 u np.where(support, u_new, u - beta * u_new) # 计算频域振幅相对误差用来监控收敛 err np.linalg.norm(np.abs(G) - amp) / np.linalg.norm(amp) err_curve.append(err) return u, err_curve obj, support, intensity make_sample(seed42) recon, errs hio_recon(intensity, support, beta0.8, n_iter300)逻辑说明np.angle(G)保留当前相位amp * exp(i·angle)完成频域幅度替换这一行所有迭代算法共用。np.where(support, u_new, u - beta * u_new)是 HIO 和 ER 在实现上的唯一区别把后半段换成np.zeros_like(u_new)就是 ER。误差计算的分子是频域振幅残差分母是测量振幅的范数它能稳定反映迭代是否收敛。参数上 beta 取 0.8 是稳妥起点n_iter 取 300 对模拟数据足够如果误差在 100 轮后还在显著下降说明初始相位或支持域有问题而不是迭代次数不够。3.4 重建验收相关系数和视觉检查模拟数据有标准答案可以用相关系数验收重建质量。一个容易被忽略的细节是相位恢复对整体旋转、共轭、全局相位常数有天然的不敏感性直接比较复数场往往得到很低的分数所以更稳妥的做法是分别比较重建幅度和真值幅度且只在支持域内计算。def evaluate(recon, obj, support): rec_amp np.abs(recon) obj_amp np.abs(obj) # 只在支持域内计算皮尔逊相关系数 cc np.corrcoef(rec_amp[support], obj_amp[support])[0, 1] # 支持域外的能量泄漏比例越小说明约束执行得越干净 leakage np.sum(np.abs(recon[~support])**2) / np.sum(np.abs(recon)**2) return cc, leakage cc, leakage evaluate(recon, obj, support) print(fsupport_mask内相关系数: {cc:.4f}, 域外能量比例: {leakage:.3%})逻辑说明np.corrcoef对两组振幅分布做归一化内积0.95 以上说明重建幅度可信leakage统计支持域外的能量占比理想情况在 1% 以下。如果 CC 低于 0.8 但 leakage 很低问题大概率在初始化或局部极小如果 leakage 偏高支持域约束没有执行干净先检查np.where的分支顺序是否写反。4. 调参核心支持域、反馈系数 beta、过采样比跑通模拟数据只证明代码逻辑没错。真正决定重建好不好看的是三个参数支持域、反馈系数、过采样比。这三个参数互相牵制一个调坏整体翻车。4.1 支持域没它不行太大太小都翻车支持域是 HIO 和 ER 唯一的空域约束它的大小直接决定解空间的收缩程度。支持域太大空域约束失去意义解会被傅里叶域约束牵着走出现孪生像和重影支持域太小把真实物体的边缘截断约束与数据矛盾误差根本降不下去。手工画矩形往往不靠谱因为真实样品很少是规则的。常见的做法是先估计再收缩用强度图的反变换估计自相关支撑。from scipy.ndimage import binary_dilation, label def estimate_support(intensity, threshold0.05, pad3): # 强度做逆傅里叶变换得到自相关图自相关的支撑约等于物体支撑的自卷积 autocorr np.abs(np.fft.ifft2(intensity)) autocorr / autocorr.max() mask autocorr threshold # 取最大的连通域避免噪声孤立点干扰 labeled, n label(mask) size np.bincount(labeled.ravel()) mask labeled np.argmax(size[1:]) 1 mask binary_dilation(mask, iterationspad) # 外扩几像素留余量 return mask逻辑说明ifft2(intensity)在理想情况下等于物体的自相关函数自相关的支撑范围比物体本身大一倍左右但至少给出了一个靠谱的上界。取主连通域后外扩 3 个像素是为了给后续 shrinkwrap 留收缩空间。这里threshold0.05是经验起点噪声大时可以调到 0.1太小会把背景噪点并进支持域。支持域不是一成不变的迭代过程中可以每 30 到 50 轮做一次 shrinkwrap用当前重建振幅的阈值重新收紧 mask。def shrinkwrap(u, threshold0.1, pad3): amp np.abs(u) thr threshold * amp.max() mask amp thr return binary_dilation(mask, iterationspad)注意 shrinkwrap 不能太激进。稀疏弱散射物体如果振幅主要集中在一个亮点阈值一高就把弱信号砍没了。我的习惯是前 100 轮用固定的大支持域等大致结构出来后再 shrinkwrap阈值从 0.1 起步每轮只收缩一点点。4.2 beta 与迭代节奏HIO 热身ER 收尾反馈系数 beta 只在 HIO 里有。beta 取 0.5 时反馈弱、收敛慢但稳定取 0.9 时收敛快但误差曲线会带明显波动超过 1.2 基本必发散结果是满屏散斑噪声。新手最容易犯的错是希望数字大一点收敛快一点结果把 beta 调到 1.5然后回头怀疑算法写错了。beta 应该配合迭代阶段变化常见做法是两段式前 200 到 500 轮用 HIO 加 beta0.7~0.9把大致结构拉出来后 50 到 100 轮切成 ER 做抛光把支持域外的杂散能量干净地压掉。ER 抛光阶段误差可能轻微抬升这不要紧视觉上重建会更干净。语气上可以这样记HIO 负责探路ER 负责收尾迭代次数不是投资进入平台期就该调参而不是加跑几千轮。4.3 过采样比为什么 2 倍是个坎过采样比指采样密度相对物体支撑的比值。理论上的关键结论是强度图必须在每个维度上对自相关支撑做到 2 倍以上过采样否则相位恢复问题信息量不足支持域约束救不回来。很多 toy 项目里模拟数据天然满足这个条件所以跑得通一到实验数据探测器尺寸、物距、波长一算过采样比不够重建就崩溃。这里有一个常见的误解有人想通过 zero padding 把强度数组变大来提高过采样比。padding 只能提高频域插值的视觉细腻度不能增加实际信息量物体的真实采样密度没有变。如果实验上过采样不够正确解法是改变成像系统放大率或改用扫描式方法收集多帧数据而不是在数据处理阶段硬补。这一点特别值得写在实验记录本第一页。参数建议值越界后果支持域初值自相关支撑外扩 1.2~1.5 倍或 3~5 px太大出孪生像太小截断物体shrinkwrap 阈值0.1~0.2 每 30~50 轮阈值过高砍掉弱散射信号HIO beta0.5~0.90.3 收敛慢1.0 发散迭代节奏HIO 200~500 轮 ER 50~100 轮平台期继续跑无收益过采样比每个维度 ≥2无法唯一重建初始化随机相位支持域内振幅取常数 1随机振幅引入高频噪声4.4 初始化随机相位之外还能做什么初始化通常用随机相位就够但有一个细节很多人忽视不要在支持域内用随机复数值初始化而是用常数振幅比如 1.0乘以随机相位。随机振幅会给初始场引入大量高频成分傅里叶域替换后这些高频分量会撑大支持域外的误差延长前期收敛时间。u np.zeros((n, n), dtypecomplex) u[support] np.exp(1j * rng.uniform(0, 2 * np.pi, sizesupport.sum()))如果对样品有额外先验比如样品的吸收已知很弱、接近纯相位物体初始化时可以把支持域内振幅设成 1.0 并对相位做小扰动这样收敛更快且更容易避开孪生像。最后记住一条朴素的工程经验固定随机种子跑出来的单个结果永远不要全信换两个 seed 跑一下如果五官不同就不是迭代次数问题是支持域或过采样比还没有调到合理区间。5. 真实踩坑记录5 个让相位恢复重建翻车的现场模拟数据一路绿灯换到真实数据就翻车这是相位恢复最常见的剧本。下面五条都是从实际项目里反复踩出来的问题每一条都按现象 → 原因 → 解决的顺序写希望能帮你省掉几轮不必要的 debug。5.1 重建结果是孪生像镜像翻转还是共轭翻转现象重建的物体轮廓、相对位置都对但细节左右颠倒或者幅度分布和真值互为镜像相关系数测出来只有 0.6 左右但人眼粗看觉得有点对。原因傅里叶域强度对物体的共轭和镜像具有对称性如果支持域约束不够紧或过采样比刚好卡在 2 附近迭代收敛到孪生像是正常现象不是算法 bug。解决最直接的办法是收紧支持域并做 shrinkwrap让空域约束把这个对称解排除掉其次是跑多次随机初始化取误差最小的一次如果数据允许用两个不同方向的照明各拍一张用共同约束打破对称性。5.2 ER 误差曲线后期像心电图纹丝不动现象ER 算法前 20 轮误差快速下降之后变成一条平线跑到 2000 轮也不变换一个初始化 seed终点误差还更高。原因ER 的硬置零投影让迭代被限制在一个局部极小附近支持域集合的非凸性决定了这个问题无法靠加迭代次数绕开。这不是参数没调好是算法本身的选择。解决切到 HIO用 beta0.8 先跑 300 轮再回到 ER 抛光 50 轮。如果必须只用 ER尝试每 50 轮对支持域做一次轻微膨胀再收缩手动把迭代从局部极小里抖出来但这不如直接换 HIO 干净。5.3 换上真实衍射图重建结果全是条纹现象模拟数据重建很干净同一套代码换真实衍射强度图重建出现横向或纵向条纹支持域内结构模糊误差曲线下不去。原因实测强度图通常有中心低频饱和、坏像素点、暗场偏置这三个问题。中心饱和会把低频分量压到一个错误量级坏像素在衍射图里表现为一个常量背景的相干叠加暗场偏置直接给振幅域加了一个直流项。解决重建前先做预处理——暗场减除、坏像素邻域插值、把中心过曝区域用 mask 盖掉不参与幅度替换。衍射强度动态范围很大习惯先取 log 显示检查一遍看到中心一片死白就说明要处理饱和。5.4 支持域手画小了能量满屏泄漏现象重建幅度图里本该为空白的区域出现明显的残余强度误差曲线长期高位震荡把支持域 mask 和重建幅度叠在一起看发现 mask 边缘正好切在物体高亮度区上。原因支持域太小把真实物体的某一部分排除在约束之外。空域约束和频域幅度约束互相矛盾算法只能在两者之间反复横跳能量自然会往支持域外泄漏。解决支持域宁大勿小。先用自相关支撑估一个偏大的 mask迭代中再 shrinkwrap 收紧。不要凭肉眼在强度图上画矩形强度图是频域量看不清物体边界。检查方法很简单把当前 mask 和重建幅度叠加显示看物体的边缘是否被 mask 硬生生裁掉。5.5 beta 越高收敛越快结果越来越像噪声现象把 HIO 的 beta 从 0.8 调到 1.5期望加快收敛结果运行几十轮后重建变成随机散斑误差曲线不降反升。原因beta 超过稳定区间后支持域外的负反馈项不再起纠正作用而是变成正反馈振荡。还有一个常见原因是更新公式写反把u - beta * u_new错写成u beta * u_new反馈立刻变负反馈同样发散。解决把 beta 拉回 0.5~0.9先确认更新分支是np.where(support, u_new, u - beta * u_new)。判断稳定性的标准是误差曲线连续 20 轮上升就直接停不要等它自己回来。想用自适应策略的话可以每 50 轮把 beta 从 0.9 逐步衰减到 0.5而不是一开始就取大值。提示以上五条有一个共同规律——先查预处理再查约束参数最后才怀疑迭代算法本身。相位恢复的坑大多不在算法在数据进入迭代之前的状态。6. 验证重建可复现误差曲线、多初始化与 FSC最后给你三个比看图说话可靠一点的验收手段也是我跑实验数据的固定流程。第一个是读误差曲线。模拟数据阶段把频域振幅误差画出来正常形状是前 50 轮快速下降、之后趋于平稳如果出现锯齿状持续波动说明支持域或 beta 边界出了问题。我现在的习惯是连续 20 轮误差不降就重新配置参数而不是加跑几千轮等它自己解套。这个习惯帮我省下了大量无意义的算力。第二个手段是多初始化。相位恢复是随机初始化下的非凸优化单次结果不具备可复现性。跑 5 到 10 次初始化取误差最小的结果作为最终重建更进一步可以在归一化对齐后对复数场取平均把随机的孪生像分量互相抵消。代价很低这是整个流程里唯一稳赚不赔的操作。第三个手段是 FSCFourier Ring Correlation实验数据专用的验证尺子。把实测数据按奇偶帧分成两组分别重建在傅里叶域的环形带里计算两组重建的相关系数以 FSC 曲线降到 1/2-bit 阈值的位置作为分辨率判据。它不仅证明你的重建不是噪声拟合还能直接报出分辨率是说服自己和审稿人的硬指标。模拟数据阶段不强制做但上真实数据前建议提前写进脚本。可能的话每次换新数据先把第 3 章的模拟自检跑一遍确认当前的参数区间仍然成立再上实测强度图实测数据第一轮只求误差曲线平稳下降不求图像漂亮。这套流程多花十分钟比对着满屏散斑猜参数值得多。希望帮到你。本文还有配套的精品资源点击获取
