简介：面向光学成像与相位恢复研究的MATLAB算法包，聚焦ER（误差减小）与Fienup HIO（混合输入输出）两种经典迭代算法，适合需要从幅度数据恢复相位信息的研究人员、学生及工程师。包内共6个文件，均为.m脚本，整体仅4KB，结构紧凑，涵盖主流程、相机信号模拟、两种核心算法实现、物理约束处理以及图像中心定位等关键模块，便于直接运行与二次开发。已有794人学习下载。通过研读这些代码，可深入理解迭代相位恢复的更新策略、约束条件与收敛机制，掌握Fienup算法在X射线衍射、电子显微镜等场景中的实际用法；同时，脚本模块化拆分清晰，也方便替换约束函数或扩展其他相位恢复变体，对学习与实验验证均有实用价值。

1. 衍射强度摆在面前，相位却丢了：HIO 和 ER 是来干什么的

只有一张衍射强度图样摆在手里，想把样品原本的复振幅分布反演出来——相位信息在记录时已经丢了，直接逆傅里叶变换只会得到一堆鬼影。HIO 和 ER 这两类 Fienup 迭代算法，正是应对这个问题最经典的起点。所谓 HIO + ER，指的是在一次重建里先用混合输入输出（HIO）跳出局部极小，再切回误差减少（ER）把解稳定下来，这是相干衍射成像、电子显微和光学显微重建里最常见的组合玩法。拿到一个命名成“HIO + ER.zip”的代码包，你要解决的大概率就是这类重建管线。这篇按 ER 原理、HIO 实现、混合策略、踩坑和验证的顺序展开，适合手上有衍射数据要出图、又不想把重建过程当黑匣子的工程师照着复现。

2. 把 ER 算法拆开看：傅里叶域换模、实空间砍支撑，为什么它容易卡死

2.1 相位恢复的问题模型与两个约束域

相位恢复用一句话描述：已知探测器记录下来的衍射强度 I(u)，希望恢复物体复振幅 f(x)，使得 |F(f)(u)|² = I(u)，其中 F 代表离散傅里叶变换。之所以不能直接做逆变换，是因为 F(f)(u) 的相位在记录那一刻就丢了，而物体的结构信息恰恰同时编码在振幅和相位里。缺少相位，|f| 的空间分布就还原不出来。

要缩小解空间，光靠强度图样不够，必须引入先验。使用频率最高的先验是有限支撑：物体只占据视场内某个有限区域，支撑域 S 之外必须为零。这样一来，整个迭代就被夹在两个约束之间：

• 傅里叶域约束：计算得到的傅里叶振幅必须替换成测量振幅 A(u)=√I(u)，相位保持自由。
• 实空间约束：估计值在支撑域 S 内可以非零，在 S 外必须清零。

Fienup 算法这个叫法，通常就是指在这两个约束之间交替投影的整套框架。它把 Gerchberg-Saxton 那种“两个平面强度互相投影”的思路，改成了“一个强度测量 + 一个支撑先验”，这让单帧衍射数据也能重建。ER 和 HIO 都是这个框架下的具体实例，区别只在实空间约束这一步怎么执行。

这里还有个容易忽视的边界条件：物体支撑的最大直径不能超过视场边长的一半。说得直白些，物体在视场里必须“小”，衍射图样才能被充分过采样，支撑约束才能把解卡住。如果物体铺满整个视场，ER 和 HIO 几乎没有恢复可能，这是实验端就要解决的问题，不是算法端能救的。

2.2 ER 的四步迭代与误差单调性

ER 算法的核心迭代体很短，单步更新不到十行：

def er_update(g, support, measured_amplitude): # g: 当前实空间估计，复数数组 # support: 布尔数组，True 表示支撑域内 # measured_amplitude: 衍射振幅，由 sqrt(强度) 得到 G = np.fft.fft2(g) # 1. 正变换到傅里叶域 amp = np.abs(G) G_new = np.zeros_like(G) valid = amp > 1e-12 # 保护除零 G_new[valid] = measured_amplitude[valid] * G[valid] / amp[valid] # 2. 换模保相：用测量振幅替换计算振幅，保留原相位 g_prop = np.fft.ifft2(G_new) # 3. 逆变换回实空间 return g_prop * support # 4. 支撑域外硬清零

四步的逻辑很清晰：先把当前估计正变换到傅里叶域，做振幅替换；逆变换回实空间后，把支撑域外的一切乘 0。第 2 步是整条管线里唯一接触测量数据的动作，第 4 步是把支撑先验压进解的动作。两步交替执行，直到误差不再明显下降。

第 4 步的硬清零，正是“误差减少”（Error Reduction）这个名字的来源：支撑域外的残余能量每轮都被直接砍掉，因此实空间误差只会下降或持平，不会上升。这个单调性是 ER 最迷惑人的地方——它的误差曲线永远是平滑下降的，看起来一切正常，实际上可能只是在某个局部极小附近越走越深，图像细节一直糊着。记住这一点，后面调试时能少走很多弯路。

代码里有个细节值得注意：对 amp 接近 0 的位置，代码直接把 G_new 置零，不做任何替换。这在大部分实验数据里够用，但如果你的衍射图样有大量暗像素，可以改成对零振幅位置注入随机相位，能略微改善暗区重建质量。

2.3 局部极小根源：傅里叶模约束集不凸

交替投影要稳定收敛到全局最优，需要两个约束集都是凸集。但相位恢复里的傅里叶模约束集 M={g: |FFT(g)|=A} 不是凸集——两个都满足“振幅等于 A”的场的线性叠加，振幅通常不再等于 A。这一条性质直接决定了 ER 的收敛结果与初始猜测强相关，不同随机种子跑出不同解是常态。

另一个容易被忽略的问题是：当物体是实值、非负的时候，ER 的支撑投影并不限制负值或虚部，解空间比真实物体大得多，更容易落进错误配置。换句话说，ER 的每一步都在做同一件事——硬掰。支撑域外有多少能量，直接清零；傅里叶振幅差多少，直接替换。这种硬投影在远离真解时效率很高，一旦靠近局部极小，就彻底失灵。

实际表现就是：误差曲线前几十次迭代快速下降，然后进入一条几乎水平的平台，恢复图像模模糊糊“像那么回事”，但细节全是糊的。这是相位恢复里最常见的翻车现象。遇到这种情况，不要盲目加大迭代次数，那只会把时间耗在同一个坑里。正确做法是换用 HIO，或者让 HIO 和 ER 交替跑，把支撑域外的能量从简单清零改成反向反馈。

3. 复现 HIO 算法：反馈系数 β 与支撑域更新的最小实现

3.1 从 ER 的硬清零到 HIO 的反向反馈

HIO 对 ER 的改动只在支撑域外这一个分支，但效果完全不同。它的更新公式是：

• 支撑域内：g_{k+1}(x) = g''_k(x)
• 支撑域外：g_{k+1}(x) = g_k(x) - β·g''_k(x)

其中 g''_k 是经过傅里叶域换模并逆变换回来的估计。支撑域内照常接受新的“输出”，支撑域外不再强制为 0，而是把当前输入减去 β 倍的输出，形成一个纠错信号。β=1 时支撑域外的新估计是 -g''_k，相当于把漏到外面的能量反向抛回迭代，这是最激进的反馈；β 取 0.7 到 0.9 时，反馈强度被压低，稳定性更好。

这个机制的直观理解是：HIO 允许解在支撑域外“临时犯错”，再用反馈把错误能量逐步泵回支撑域内。相比 ER 的硬清零，HIO 不会因为一次错误的支撑投影就把解卡死在某个状态，这也是它能跳出一部分局部极小的原因。代价是 HIO 不再保证误差单调，误差曲线会带起伏，迭代后期甚至小幅振荡，所以它通常不单独用，而是和 ER 搭配。

3.2 Python 实现 HIO 核心迭代代码

HIO 的单步更新代码和 ER 非常像，只改最后一行：

def hio_update(g, support, measured_amplitude, beta=0.9): # beta: 反馈系数，典型取值 0.7 ~ 1.0 G = np.fft.fft2(g) amp = np.abs(G) G_new = np.zeros_like(G) valid = amp > 1e-12 G_new[valid] = measured_amplitude[valid] * G[valid] / amp[valid] g_prop = np.fft.ifft2(G_new) g_next = np.empty_like(g) g_next[support] = g_prop[support] # 支撑域内沿用 ER 行为 g_next[~support] = g[~support] - beta * g_prop[~support] # 支撑域外：输入减去 beta 倍的输出，形成反馈 return g_next

换模部分与 ER 完全一致，真正不同的是最后两行。支撑域内直接取 g_prop，支撑域外用 g - β·g_prop 做反馈。注意这里用的是布尔索引，support 和 ~support 分别是支撑域内外的两个掩膜，这也要求支撑域必须是一个明确的二值数组，而不是模糊的权重。

初始化方式对结果影响很大。我常用的初始化是：在支撑域内生成随机复数场，支撑域外置零。

rng = np.random.default_rng(0) g = (rng.standard_normal(intensity.shape) + 1j * rng.standard_normal(intensity.shape)) g *= support

固定随机种子很重要。同一份数据，换个种子可能得到不同结果，这是非凸问题的固有属性。调试时固定种子，确认算法本身没问题后，再放开种子做多轮启动、挑误差最小的一次输出。

3.3 β 怎么调：0.7、0.9 和 1.0 的差异

β 是 HIO 里唯一值得认真调的超参数，其他参数都相对皮实。我按使用频率排序给出经验值：

• β=0.9：默认起点。反馈够强，能较快跳出局部极小，振荡幅度在可接受范围。
• β=0.7：更保守。迭代平稳，噪声容忍度好，适合信噪比低的衍射数据，但收敛偏慢。
• β=1.0：反馈拉满。前期收敛速度最快，但误差曲线锯齿明显，后期容易在真解附近来回振荡不收敛。

有一种做法是给 β 做退火：前期用 0.9 快速探索，跑两三百轮后降到 0.5 到 0.6 稳定细节。这个策略在低信噪比数据上比固定 β 更稳，代价是多一组超参数要调。

支撑域本身的更新也很关键。如果初始支撑给得不准，可以每隔 20 到 50 轮用当前重建振幅重新估计支撑域：

def refresh_support(g, threshold=0.15, dilate_pixels=1): amp = np.abs(g) amp_max = amp.max() new_support = amp > threshold * amp_max from scipy.ndimage import binary_dilation return binary_dilation(new_support, iterations=dilate_pixels)

这个手法叫 shrinkwrap，本质是用当前估计的轮廓反推支撑，再外扩几个像素防止过紧。它能让支撑域从“宽松但安全”逐步收紧到“紧凑且准确”，是处理未知支撑时最可靠的手段。注意：外扩像素不能省，支撑域一旦过紧，等于在每一轮都砍掉真实信号的一部分，后面怎么迭代都补不回来。

4. HIO + ER 混合策略：先 HIO 跳出局部极小，再 ER 把解焊死

4.1 纯 HIO 为什么在真解附近振荡

HIO 不保证误差单调，这个特性既是优点也是缺点。优点在前探索阶段明显——它能从局部极小里爬出来；缺点在后期暴露——接近真解时，支撑域外那些微小残差会被反馈项反复放大和修正，误差曲线在一个低水平上锯齿状波动，解在真解附近来回抖动，始终达不到高精度。

ER 恰好是另一个极端：它单调收敛，但容易在局部极小处卡死。不过有一个关键经验：如果 ER 的起点已经离真解足够近，它的单调性会把解稳稳“焊”在真解上。这个“足够近”不需要很精确，HIO 跑出来的解通常就满足条件。

所以 HIO + ER 的混合策略本质上是分工：HIO 负责在全局范围内探索，把解带到一个靠谱的盆地；ER 负责在盆地里做精细收敛，消除 HIO 的锯齿尾巴。这个组合比单独跑任何一个算法都稳定，也是大量开源实现里默认的流程。

4.2 切换判据：实空间误差进入平缓期

HIO 跑多少轮再切 ER，最稳妥的判断方式不是凭感觉定死一个轮数，而是监控实空间误差的变化率。实空间误差的定义是支撑域外能量占比：

def real_space_error(g, support): outside_energy = np.sum(np.abs(g[~support]) ** 2) total_energy = np.sum(np.abs(g) ** 2) + 1e-12 return np.sqrt(outside_energy / total_energy)

当这个误差连续一段迭代不再明显下降，就说明 HIO 的探索阶段已经结束，该切 ER 了。这里有个细节：由于 HIO 误差会锯齿波动，不能看单次值，要看滑动平均。我通常每 10 轮记录一个点，然后比较最近一段与上一段的均值变化。

def is_plateau(errors, window=3, tol=1e-3): # errors 是每 10 轮记录一次的误差序列 if len(errors) < 2 * window: return False recent = errors[-window:] prev = errors[-2 * window:-window] return abs(np.mean(recent) - np.mean(prev)) / (abs(np.mean(prev)) + 1e-12) < tol

这里 window=3 表示比较最近 3 个记录点（约 30 轮）与之前 3 个点（约 30 轮）的平均误差，变化率小于 0.1% 就认定进入平台期。阈值不用太严，HIO 本身就会振荡，抓大趋势即可。

4.3 混合流程的完整函数与调用方法

把上述逻辑拼成一个可直接调用的函数：

def hio_plus_er(intensity, support, beta=0.9, hio_iters=400, er_iters=100, plateau_switch=True, seed=0): measured_amplitude = np.sqrt(intensity) rng = np.random.default_rng(seed) g = (rng.standard_normal(intensity.shape) + 1j * rng.standard_normal(intensity.shape)) g *= support errors = [] for i in range(hio_iters): g = hio_update(g, support, measured_amplitude, beta) if i % 10 == 0: errors.append(real_space_error(g, support)) if plateau_switch and i > 60 and is_plateau(errors): break # 平台期提早结束 HIO 阶段 for i in range(er_iters): g = er_update(g, support, measured_amplitude) if i % 10 == 0: errors.append(real_space_error(g, support)) return g, errors

参数说明：hio_iters 建议 200 到 500 轮，模拟数据 300 轮足够，实验数据可能要上千轮。er_iters 50 到 100 轮就够，ER 在局部收敛很快，跑太多轮除了耗时没有明显收益。plateau_switch 默认开启，会在误差进入平台期时提前切到 ER。返回的 errors 列表可以直接画图，用于判断这次重建是否正常。

还有一个我常用的变体：如果 HIO + ER 跑完一轮后傅里叶误差仍然偏高，不必从头随机重启，而是把当前解当作新初始值，再跑一轮短 HIO（100 轮左右），通常能把误差再压下去一截。这种“引导重启”比完全重新随机初始化有效得多，建议在正式重建流程里加上。

5. 相位恢复避坑记录：从解压代码包到看懂误差曲线的 5 个翻车点

5.1 压缩包解压报错 could not find eocd：先验货再动手

现象：下载的“HIO + ER.zip”双击解压就直接失败，解压工具报 invalid zip archive，错误信息里带着 could not find eocd。

原因：eocd 是 zip 格式的尾部中央目录记录，找不到它基本等于文件被截断。八成是下载过程断流、网盘中转时文件不完整，或者浏览器缓存了部分内容，跟算法本身没有任何关系。

解决：先不要急着换解压工具，用命令行验证文件完整性。Linux 和 macOS 下直接跑unzip -t HIO.zip，Windows 下用tar -tf HIO.zip或 7-Zip 的测试功能。如果测试失败，重新下载，并且对比文件字节数是否与发布方标注一致。文件完整但依然报错，再考虑是不是压缩包本身加密或损坏。

5.2 支撑域给太大，误差“假收敛”图像一团糊

现象：HIO 阶段误差曲线下降很快，切入 ER 后也一路平滑，最终误差值挺好看，但恢复图像细节模糊、边界发散，看起来像蒙了一层雾。

原因：支撑域给得比真实物体大太多，实空间约束形同虚设，解可以在大半个视场内自由分布。误差低是因为能量分散开了，不是因为结构收敛了。很多初次跑实验数据的人，习惯把视场中央一个大矩形当成支撑，结果就是这种“假收敛”。

解决：用强度图的自相关来估计初始支撑。自相关 c(x)=F⁻¹{I}(x) 的显著区域宽度，约等于物体支撑直径的两倍，取一半再做一圈膨胀，就是个很稳的初始值。配合 shrinkwrap 在迭代中逐步收紧。另外检查过采样条件：物体支撑的线性尺寸尽量不要超过视场边长的一半，这是唯一性能否成立的前提。

5.3 把衍射强度 I 直接当振幅 |F| 用，低频直接失真

现象：重建结果出现巨大的直流分量，物体中心一片亮斑，周围结构对比度极低，换什么算法都救不回来。

原因：代码里把探测器读出的强度 I 直接当成了振幅 A 代进换模步骤。物理量搞混了，强度是振幅的平方，低频部分被平方后会被严重放大，换模时低频权重失衡，迭代等于在错误的约束上跑。

解决：在任何代码入口处明确区分两个变量，命名上强制体现物理量。我习惯把探测器数据叫 intensity，然后用measured_amplitude = np.sqrt(intensity)单独算一次，之后所有换模操作只允许碰 amplitude，不让 intensity 进入迭代函数。这一步几乎是所有“看起来哪都对就是结果不对”案例的头号原因。

5.4 中心亮斑压过细节，迭代被低频牵着走

现象：误差曲线在持续下降，但重建图像整体偏灰，边缘和精细结构出不来，好像迭代全在“照顾”中心最强的那几个像素。

原因：衍射图样中心是零频和低频分量，强度比高频细节大几个数量级。换模时振幅直接替换，低频误差在总误差里占绝对主导，HIO 的反馈几乎全花在修正低频上，高频细节得不到足够的梯度。

解决：两个方向。一是实验端，在探测器中心放一个小遮挡块（beamstop）挡住零级光，很多 CDI 实验就是这么做的；二是算法端，对误差函数里的不同频率做加权，或者在做换模前对振幅谱做一次白化处理，压制低频主导。最省事的做法是先把强度图的低频部分做局部归一化，再进入迭代。注意这一条和 5.3 不同，5.3 是用错了物理量，这里是物理量对了但权重失衡。

5.5 HIO 切 ER 过早，误差曲线当场反弹

现象：HIO 跑了不到 100 轮就看曲线“差不多平了”，切到 ER 后前几十轮还行，后面误差突然反弹，甚至比切换前还高。

原因：HIO 的误差曲线在探索阶段也会出现短暂平台——它可能卡在一个局部极小里还没跳出来。这时切到 ER，等于用 ER 的单调性把这个局部极小焊死了。ER 的反弹其实不是误差上升，而是它把 HIO 从错误盆地拉出来之前的挣扎，但结果一样是坏的。

解决：切换前加一个硬性的最小迭代数，比如至少 200 轮 HIO 才允许切。同时把平台期判断窗口放宽，别用 10 轮 20 轮的短窗口。我一般用 30 轮以上的窗口做滑动平均比较。切到 ER 之后如果误差在前 20 轮明显上升，马上中断，换一个随机种子重启 HIO，不要硬等 ER 自己缓过来。

6. 自检脚本：用模拟数据验证你的 HIO + ER 实现能不能真的收敛

6.1 生成模拟物体与衍射强度

调试任何重建算法，第一件事不是碰实验数据，而是先在模拟数据上跑通。模拟物体的生成和衍射强度计算只需要几行：

size = 128 obj = np.zeros((size, size), dtype=complex) obj[44:84, 44:84] = np.random.default_rng(1).standard_normal((40, 40)) intensity = np.abs(np.fft.fft2(obj)) ** 2 support = np.zeros((size, size), dtype=bool) support[44:84, 44:84] = True

物体是 40×40 的随机区块，放在 128×128 视场中央，线性过采样比约 3.2，面积比约 10，满足唯一性条件。用这个设置跑 hio_plus_er，正确的实现在几百轮内就能把支撑域外误差压到接近零。

6.2 三个定量指标判断收敛

单看误差曲线不够，我习惯同时盯三个指标，做完一次重建后把它们全部打印出来：

相关系数需要回归原物体对比，注意相位恢复存在平移、共轭、整体相位旋转移位等平凡歧义。对比前先让恢复结果在图像上对齐，再取振幅算相关，否则相关系数会被这些歧义拖低，让人误以为算法有问题。

6.3 一个收尾习惯

我现在的习惯是：所有重建脚本都内置一个--self-test开关，跑实验数据前强制先生成一组模拟物体，用 HIO + ER 完整跑一遍并打印三个指标。全部达标才碰真实数据。这个习惯帮我挡掉了至少一半的“算法不收敛”问题——大多数时候不是算法的问题，是支撑域给错、振幅和强度搞混、或者切 ER 时机不对。把自检流程做进管线里，比每次手动查那些坑省力得多。希望帮到你。

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