简介：面向光学成像与波前测量方向的Matlab仿真工具包，适合高校相关专业学生、科研人员及对相位恢复算法感兴趣的开发者。内容围绕Gerchberg-Saxton迭代算法展开，完整演示从强度数据反演相位、再到波前重建的核心流程，覆盖空间域与频域交替迭代、傅里叶变换及衍射模拟等关键环节，可用于全息成像、显微成像和光学校正等场景的实验验证与教学演示。压缩包共30个文件，以m脚本与fig图像为主，辅以jpg预览图、bin数据文件、说明文档及备份文件，整体大小约53.55MB，结构清晰，便于按流程运行和二次修改。已有63人学习使用。资料包含可直接运行的GS算法核心代码、多种光场仿真脚本、结果图像与理论说明文档，并提供了参数调整接口与性能评估思路，帮助使用者深入理解算法原理，同时提升在光学系统模拟中的实际操作与问题分析能力。
1. 光学相位恢复为什么绕不开Gerchberg-Saxton：从迭代算法到仿真平台
做波前重建的人迟早要被相位问题卡住：探测器只记录光强，相位在曝光瞬间就丢了，可波前重建、光束整形、相干成像全都指着相位过活。Gerchberg-Saxton迭代算法，业内习惯叫GS，是解决这个问题最朴素也活得最久的方案——在两个傅里叶相关的平面之间来回投影，轮流让振幅去凑测量值，相位在迭代里一点点逼近真实波前。跑过几轮你就会发现它不挑物理模型，这正是光学仿真平台普遍内置GS求解器的原因。这篇笔记从GS的数学骨架拆到Python实现、参数选择和收敛失败的排查，适合刚开始做光学相位恢复、以及被不收敛反复折腾的人。
2. Gerchberg-Saxton的收敛骨架：交替投影、约束条件与唯一性边界
2.1 物面和频域两个平面之间的四步循环在做什么
GS算法的出发点是两个平面之间的傅里叶变换关系。物平面上的复振幅分布 (g(x)=A(x)e^{i\phi(x)})，经过聚焦透镜或者远场衍射之后，在频域平面上变成 (G(u)=B(u)e^{i\psi(u)})，两者之间是严格的一对一傅里叶变换。探测器能拿到的是强度 (|G(u)|^2)，开根号之后就是频域振幅 (B(u))。物面振幅 (A(x)) 在多数实验里也知道——比如入射光束是高斯光或者平面波。问题是物面相位 (\phi(x)) 拿不到，频域相位 (\psi(u)) 也拿不到，要从四个量里已知的两个反推另外两个，这本身是个欠定问题。
GS给出的解法非常直观：先猜一个物面相位，做一次傅里叶变换到频域，把算出来的频域幅度强行替换成实测值，相位留着；然后做逆傅里叶变换回物面，再把物面幅度强行替换成已知入射值，相位继续留着。这两步替换加上两个傅里叶变换，构成一个循环，跑几十轮之后，物面相位就会稳定在一个解上。
从数学角度看，这个循环是在两个集合之间做交替投影。第一个集合是「所有经傅里叶变换后幅度等于实测频域幅度的复场分布」，第二个集合是「所有物面幅度等于已知入射幅度（或满足支撑约束）的复场分布」。GS每轮做的事情，就是往这两个集合上各投影一次。如果两个集合都是凸集，交替投影理论保证收敛到交集；可惜物面幅度约束这个集合不是凸的，所以GS不保证收敛到全局最优，却在实际问题里出奇地耐用。
做仿真平台时，我习惯把这个循环拆成四个算子而不是一个整体函数：正傅里叶变换、频域振幅替换、逆傅里叶变换、物域约束。拆开的好处是后面换近场衍射时，只替换第一和第三个算子，替换振幅和约束的逻辑完全不动。
2.2 支撑约束与幅度约束怎么选：GS变体的分水岭
上一节说的物面幅度替换，有一个前提假设——物面振幅分布已知。这在纯光束整形问题里成立：入射光是什么形状、透过率多少，都是事先标定过的。但相位恢复实验里这个前提经常崩掉。你想恢复的是一个样品透射后的波前，样品内部有吸收，出射振幅变成未知量，这时候还强制物面幅度等于入射值，相当于把样品吸收信息消灭了，恢复出来的相位一定是错的。
这种情况下能用的约束退化成支撑约束：在支撑集内部，复振幅不做任何限制；在支撑集外部，复振幅强制为零。支撑集就是样品在物面上占据的区域，可以从可见光图像或者衍射图样的自相关函数估出来。支撑约束比幅度约束弱得多，但它保留了对样品复振幅的全部自由度，这也是相干衍射成像领域敢用GS系列算法恢复复杂样品的原因。
频域一侧的约束则基本固定为幅度替换，因为探测器直接给强度。不过这里有一个工程细节：强度值开根号之前要做背景校正，否则CCD暗电流和读出噪声会被开根号操作扭曲，低频区域出现一层缓慢起伏的伪纹。仿真平台里如果直接拿原始强度图开根号，迭代结果经常是「看起来收敛了但波前形状不对」，这类问题几乎都出在强度预处理这一步。
2.3 为什么GS容易停滞：局部极小值与Fienup的HIO解法
GS在实际使用中最大的问题是停滞。所谓停滞，是指能量残差不再下降，但恢复出的相位和真实值还差着一截。原因是物面幅度约束构成的集合非凸，交替投影可能落进局部极小点，而这个局部极小点和真实解之间的距离没有上界，有时候差0.1 rad，有时候差1 rad，全看初始相位猜得怎么样。
Fienup在1982年提出的混合输入输出算法（HIO）正是冲着这个缺点来的。HIO把物面更新公式从「替换」改成混合更新：支撑集内部用新计算出的场值，支撑集外部用上一轮场值的反向分量，相当于给迭代加了一个「反弹力」。这个改动让迭代在陷入局部极小的时候获得一个扰动，有机会跳出来继续搜索。仿真平台上最常见的搭配是先HIO跑20到30轮跳出停滞区，再切回GS做精细收敛。
需要注意HIO不是直接替代GS，它在快速跳出局部极小之后，收敛行为反而没有GS稳定，所以工程上把两者串起来用是常态。切换时有一个关键动作：HIO最后一步的复场相位必须作为GS的初始相位传入，而不是重新随机初始化，否则前面的搜索全部白费。
3. 在Python里跑通GS仿真平台：从NumPy核心循环到角谱传播
3.1 二十行代码实现物面到频域的振幅替换
先给一个最干净的GS核心实现，它假设物面振幅已知，频域目标振幅已知，迭代100轮返回恢复相位。这段代码就是整个仿真平台的地基。
import numpy as np def gs_core(input_amp, target_amp, support=None, iterations=100): """ GS 核心迭代：交替满足物面与频域振幅约束 input_amp : 物面振幅分布, shape=(H, W) target_amp: 频域目标振幅, shape=(H, W)，由实测强度开根号得到 support : 物面支撑掩模, 1 为有效区域，0 为强制置零区域 iterations: 迭代轮数 """ h, w = input_amp.shape # 初始相位：均匀随机分布在 [0, 2pi) 区间 phase = 2.0 * np.pi * np.random.rand(h, w) field = input_amp * np.exp(1j * phase) for _ in range(iterations): # 正向傅里叶变换 spectrum = np.fft.fft2(field) # 频域振幅替换，相位保留 spectrum = target_amp * np.exp(1j * np.angle(spectrum)) # 逆向傅里叶变换回物面 field = np.fft.ifft2(spectrum) # 物域约束：先乘支撑掩模，再替换振幅 if support is not None: field = field * support field = input_amp * np.exp(1j * np.angle(field)) return np.angle(field)逻辑说明：第5行到第8行完成初始随机相位猜测；第10行到第20行是四步循环。第13行把频谱幅度强制替换成实测值，这是GS的核心动作——只改幅度、不碰相位；第16行逆变换回物面；第18行先乘支撑掩模，把支撑外残留值清零；第19行再把幅度替换成入射值，相位保留。第18行和第19行的顺序是有讲究的：如果先替换幅度再乘支撑，支撑外的值会被幅度归一化重新拉大，等于没约束干净。
参数说明：input_amp和target_amp必须归一化到同一能量尺度，否则替换振幅时会引入系统性能量跳变，迭代过程出现锯齿状波动。iterations取100对256×256阵列在普通PC上大约几秒，但这个函数没做频域归一化，实际平台里要用fftshift统一坐标原点，我后面会专门提这个坑。
3.2 把纯FFT换成角谱法：近场衍射与采样间隔的共同约束
上一版实现隐含一个光学假设：物面和频域面之间是严格的傅里叶变换关系。这个假设只在夫琅禾费远场或者透镜焦平面成立。仿真平台如果要对近场衍射建模，FFT这个传播子就需要替换成角谱法。角谱法和傅里叶变换的差异在于多了一个传递函数，这个传递函数描述了自由空间里各空间频率分量的相位积累。
def angular_spectrum_propagate(field, wavelength, distance, dx, dy): """ 角谱法近场传播算子 field : 输入复振幅, shape=(H, W) wavelength: 波长, 单位 m distance : 传播距离, 单位 m dx, dy : 物面采样间隔, 单位 m """ ny, nx = field.shape fx = np.fft.fftfreq(nx, dx) fy = np.fft.fftfreq(ny, dy) FX, FY = np.meshgrid(fx, fy) k = 2.0 * np.pi / wavelength # 角谱传递函数中的开根号项 temp = 1.0 - (wavelength * FX)**2 - (wavelength * FY)**2 temp[temp < 0] = 0 # 倏逝波分量在自由空间中不传播 H = np.exp(1j * k * distance * np.sqrt(temp)) spectrum = np.fft.fft2(field) return np.fft.ifft2(spectrum * H)逻辑说明：fftfreq生成频率坐标，单位是1/m；temp对应角谱理论里的传播条件，当空间频率的平方和超过1/λ²时变成负值，这些高频分量对应倏逝波，不能传播到远场，直接置零；传递函数H的相位就是平面波走过distance距离的积累。用这个函数替换gs_core里的fft2和ifft2，就得到近场版GS。
参数说明：distance必须精确到波长量级，误差超过一个波长，相位梯度会整体偏移。dx和dy的取值有一个硬限制：必须小于λ/(2·NA)，否则高频空间频率分量收敛不到真实值，恢复相位表面会出现细密的水波纹。这个限制和奈奎斯特采样定理是一回事，它在仿真平台里表现为——缩小像素尺寸不总是好事，像素太密导致阵列太大，频域里有效信息的占比萎缩，迭代效率反而下降。
3.3 三个输出指标：RMSE、能量残差与有效区域信噪比
仿真平台最终给用户的是一个相位矩阵，但矩阵本身不能说明好坏。实际项目里我同时挂三个指标，分别监控不同层面的质量。
第一个是物面均方根误差RMSE，只在做仿真验证（已知真实相位）时才有意义。第二个是频域能量残差，定义为每次迭代后模拟频域幅度与实测幅度之差的平方和，除以实测幅度平方和。这个指标不需要知道真实相位，能实时观察迭代是否在收敛。第三个是有效区域信噪比，把恢复相位在支撑集内部的均值与支撑集外噪声的标准差做比值，用来衡量相位图的对比度是否足够支撑下游处理。
能量残差有个典型的陷阱：它会在前20轮快速下降，然后在某轮突然反弹。反弹的原因是物面振幅约束把支撑集外残留的能量重新注入频域，频域替换又把这部分能量压回去，两个约束打架。如果反弹幅度超过前20轮下降幅度的20%，基本可以判定振幅约束条件与实验实际情况不匹配，需要把物面约束切换成纯支撑约束再看一轮。
4. 参数工程：像素尺寸、支撑集与初始相位怎么调到收敛
4.1 阵列尺寸和采样间隔：奈奎斯特限制下的网格设计
仿真平台里最常被忽略的是网格设计。许多人直接从实验装置拷贝像素尺寸，阵列大小随便填一个2的幂次，结果跑出来的相位要么边缘振铃、要么中心区域信息挤成一团。网格设计要同时看两个量：采样间隔决定了频域覆盖范围，阵列尺寸决定了频域分辨率。
具体来说，物面采样间隔dx对应的最大频域范围是1/(2dx)，这决定了能恢复的最高空间频率。样品特征尺寸a对应的频域宽度约为1/a，这个值必须落在1/(2dx)之内，最好只占中间三分之一区域。如果占满了整个频域，高频信息在频域阵列边缘被截断，恢复相位外圈会出现一圈环形伪影，而且无论如何增加迭代次数都消不掉。
阵列尺寸N的选择则看频域分辨率需求。频域采样间隔是1/(N·dx)，它决定了恢复相位在频域里的像素化程度。N太小，频域振幅约束变得粗糙，GS每轮替换幅度时引入的量化误差大，最终相位普遍偏平滑、细节丢失。经验值是让N·dx至少等于样品尺寸的4倍，这样频域中间还能留出足够的过采样倍数来抑制伪影。
4.2 支撑集宽度：改两个像素就能让结果翻车
支撑集是GS里最敏感的参数，我在这上面踩过真金白银的坑。支撑集设得比样品实际尺寸大，支撑外的弱散射信息被当作噪声保留下来，恢复相位边缘会出现一圈强度渐变的光晕；支撑集设得比样品实际小，样品的低频信息被截断，恢复相位变成模糊的椭圆形，两条路径都会让最终波前失真。
比较稳健的做法是从衍射图样的自相关函数估支撑集。强度图的自相关宽度大约是样品自相关宽度的两倍，所以把强度图自相关函数大于零的区域缩小到一半，可以作为支撑集初值。得到初值之后先跑20轮GS，观察恢复相位在支撑边缘是否有剧烈跳变。有跳变说明支撑集偏小，往外扩3到5个像素；没有跳变但边缘区域出现均匀模糊，说明支撑集偏大，往里收2到3个像素。这个「先粗跑、再微调」的流程，比任何理论公式都实用。
支撑集还有一个容易被忽略的性质：它必须是非对称的。如果样品本身是中心对称图形，GS恢复出的相位会陷入镜像简并，两个对称解轮流出现，迭代曲线呈周期性震荡。解决办法是把支撑集做一个明显的不对称偏移，或者干脆在物面一角加一个已知的标记物，打破对称性。
4.3 初始相位随机种子与迭代次数：统计策略比玄学可靠
初始相位是随机给的，这一度让GS被戏称为玄学算法。同一个强度图，不同随机种子跑出来的相位差异客观存在，差异幅度和约束条件的强度成反比。约束越强（比如物面幅度完全已知），随机种子影响越小；约束越弱（只有支撑约束），随机种子影响越大，极端情况下可以差出0.5 rad。
我的平台里对付这个问题的方法是批量初始化：一次跑8到16个随机种子，每个种子迭代200轮，然后取能量残差最低的那个作为输出。这个方法本质上是把随机初值的方差摊平到并行计算里，代价是计算量增加，但在现代机器上256×256阵列跑16个种子也就一两分钟，性价比远高于单个种子迭代2000轮。迭代轮数本身不需要贪多，200到300轮之后能量残差的下隆曲线已经平坦，继续跑只是在局部极小里原地踏步。
迭代过程的监控不能只看残差数值，还要看相位图在每50轮之间的变化。如果两轮快照之间相位差异仍然明显，说明还没收敛；如果相位图几乎不动但残差在下降，说明在做慢速精修，可以提前停止。
5. 避坑与排查：GS恢复不收敛的五个典型现场
5.1 相位图冒出一片椒盐噪点
现象：恢复相位在有效区域内部均匀分布着细密的高频噪点，像照片里的颗粒感，能量残差曲线仍然单调下降，看起来一切正常。
原因：频域目标幅度里残留了探测器读出噪声。开根号操作会让低强度区域的相对噪声变大，GS每轮做振幅替换时，这些噪点被反复注入相位，最终均匀扩散到物面。能量残差下降是因为噪声成分也参与了幅度匹配。
解决：在开根号之前对强度图做一次3×3中值滤波，再做背景减除；如果噪点仍然明显，把GS里频域替换的幅度从硬约束改成软约束——每次替换时按75%新值加25%旧值混合，给迭代一个平滑缓冲。这个方法会让收敛变慢，但能明显压低高频噪声。
5.2 能量残差停滞但波前还在缓慢漂移
现象：能量残差在某个数值附近震荡，不再下降；但是每隔100轮保存的相位图对比，波前边缘的等位线还在缓慢移动，说明状态没真停，而是在局部极小附近游荡。
原因：物面幅度约束太硬，把迭代压在一个非凸集合的局部极小里出不去。这种情况在只有支撑约束的GS里尤其常见。
解决：切到HIO跑20轮再切回GS。HIO的反弹项会给迭代一个持续的扰动，大概率能把状态从局部极小里推出来。如果切换一次不够，就在切换时机上做随机化——让HIO在支撑集里也保留一小部分上一轮值，制造更随机的跳跃路径。
5.3 恢复相位出现镜像孪生像
现象：恢复出的相位分布关于坐标原点呈中心对称，像是原相位绕原点转了180度。两个对称解在迭代过程中交替出现，导致能量残差曲线呈周期性起伏，相位图在两次快照之间整体翻转。
原因：频域振幅和物面振幅只能确定相位到「相位φ(x)和相位-φ*(-x)同为解」的程度。这是傅里叶相位恢复问题本身的简并性，不是算法缺陷。物体越接近中心对称，这个现象越严重。
解决：在物面支撑集里人为制造不对称性。最简单的方法是把支撑集沿一个方向平移几个像素，让样品偏离坐标原点；或者用倾斜照明，使频域强度图整体偏移，打破对称。平台里如果跑出来孪生像，优先检查样品是否被居中摆放了。
5.4 支撑集微调后结果剧烈跳变
现象：支撑集只往外扩了3个像素，恢复相位结构彻底变了，之前看起来合理的波前一下子变成另一幅面孔。
原因：支撑集的微小变化改变了非凸集合的边界，局部极小的位置发生跳变，迭代从原来的极小值点滑落到另一个极小值点。这不代表扩支撑集是错误的，只是当前迭代状态对约束边界敏感。
解决：支撑集变化之后不要延续旧的迭代状态，而是重新初始化相位再跑一轮。同时在微调支撑集时，把扩大的像素先用渐变权重纳入约束，第一轮乘0.3、第二轮乘0.6、第三轮才完全生效，给迭代一个平滑过渡期。我后来把这段逻辑做成了平台里的「支撑集退火」按钮，基本消除了这类跳变。
5.5 仿真平台输出达标但实验对不上
现象：仿真里恢复的相位RMSE很低、残差收敛漂亮，但同样的参数套到实验数据上，恢复波前和独立测量的结果明显不符。
原因：仿真里没有建模真实光学系统的误差源——畸变、像差、照明不均匀、CCD非线性响应。GS在仿真里收敛，恰恰利用了目标幅度是理想参数生成这一便利；实验数据里低频背景不均匀，相当于给强度图叠了一个未知增益分布。
解决：在仿真平台里预置一个失真模型：给目标幅度乘上一个缓变的空间增益函数，再叠加高斯噪声。跑通失真模型之后，再回到实验数据时就不会措手不及。此外，实验数据在进GS之前要做平场校正和暗场减除，这两个步骤省略掉，再好的迭代算法也救不回来。
6. 把GS平台从「能跑」变成「可信」：已知相位验证集与HIO混合求解
6.1 用已知相位分布做回归验证
平台搭好之后第一件事不是接实验数据，而是做回归验证。生成一组已知相位分布（比如泽尼克多项式的前几项），用它和设定好的入射振幅合成物面场，正向传播生成频域强度，然后把强度作为GS的输入，看恢复相位和标准答案差多少。这个验证能一次性暴露网格设计、支撑集宽度、振幅归一化三类问题，比对着实验数据盲调快得多。我的习惯是保存标准答案和恢复结果的相位差矩阵，画成伪彩图查看残差的空间分布——如果残差集中在边缘，支撑集偏小；如果残差呈条纹状，采样间隔不达标。
6.2 引入HIO分支处理停滞，GS只做最后精修
验证通过之后，把平台的求解器升级为GS/HIO混合模式。整体流程是：HIO跑30轮作为粗搜索，选取能量残差最低的快照，把它最后的物面复场作为GS的初始值，再跑200到300轮精修。这个组合把HIO的跳出能力和GS的稳定性都用了上，对付局部极小和孪生像我至今没遇到过失效的情况。迭代过程中保留每一轮的相位快照，方便事后追溯是哪个环节引入了异常。
说句实在话，GS这个算法本身并不神秘，它的价值全在工程细节里——域约束怎么设、支撑集怎么估、初始相位怎么初始化、传播算子怎么换。把这些细节写进平台并逐条验证，比追求更花哨的新算法更管用。我第一次把这套流程跑通时花了三个晚上，几乎全耗在支撑集微调引发的跳变上；后来验证集成了标准动作，再没让这类问题过夜。希望这份排查思路能帮你少走这段路。
本文还有配套的精品资源，点击获取
