Gerchberg-Saxton相位恢复算法：从交替投影到Python仿真实践

简介面向光学成像与波前测量方向的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, supportNone, iterations100): 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/mtemp对应角谱理论里的传播条件当空间频率的平方和超过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这个算法本身并不神秘它的价值全在工程细节里——域约束怎么设、支撑集怎么估、初始相位怎么初始化、传播算子怎么换。把这些细节写进平台并逐条验证比追求更花哨的新算法更管用。我第一次把这套流程跑通时花了三个晚上几乎全耗在支撑集微调引发的跳变上后来验证集成了标准动作再没让这类问题过夜。希望这份排查思路能帮你少走这段路。本文还有配套的精品资源点击获取