简介基于Gerchberg-Saxton迭代算法的光学相位恢复与波前重建仿真平台以Matlab为运行环境面向光学工程、物理及相关专业的学生与研究人员用于演示光强数据反演相位分布、重构波前传播形态的核心过程。资源包共30个文件主体为M脚本与FIG图像包含GS算法、干涉与衍射仿真、目标相位重建等模块程序另有JPG结果图、BIN数据文件及TXT/MD说明文档压缩包大小53.55MB结构清晰便于按需调用。已有63人浏览学习。通过修改输入参数或算法结构可观察不同条件下的恢复效果与误差表现适用于光学课程实验、全息成像与光学校正等场景的模拟验证。配套说明文件提供操作引导与理论注释有助于快速上手并深入理解GS算法的空间-频域迭代机制。1. 从一张强度图里“猜”出相位GS算法为什么值得做成仿真平台在无透镜成像、自适应光学和定量相位显微里你大概率绕不开一个名字Gerchberg-Saxton迭代算法。它解决的是光学相位恢复中最朴素也最磨人的那个问题——探测器只能记录光强相位信息在记录瞬间就丢了怎么把波前重建出来GS算法用“空间域振幅约束 频率域振幅约束”交替投影的方式把这个反问题变成一个可迭代的数值求解过程。但算法是一回事能稳定复现又是另一回事。多数人第一次跑GS都翻过车迭代几十轮不收敛、重建相位自带条纹、换张图就全面失效。与其在实验台上反复撞墙不如先把波前重建仿真平台搭起来——在电脑上把物理引擎、迭代求解器、收敛判据和评价指标一次配齐再用同一套代码去对接实验数据。这篇文章照着这个思路给你一条能从零跑到批量对比实验的完整路径。2. 相位恢复为什么绕不开GS单帧强度反演相位的数学地基2.1 相位丢失光学测量里被扔掉的那一半信息一个标量光场传播到探测器平面时完整的复振幅可以写成$$U(x,y) A(x,y)\cdot e^{j\varphi(x,y)}$$其中振幅 $A(x,y)$ 直接决定了光强 $I |A|^2$探测器能记录到的只有这个强度相位 $\varphi(x,y)$ 携带了物体厚度、表面形貌、波前畸变等信息却在光电转换的一瞬间彻底丢失。于是我们面对的是一个典型的相位恢复问题已知或能测到输入面光强分布以及经过一段距离衍射或经过傅里叶变换后的频谱面光强分布求原始波前的相位。GS算法的核心洞察在于虽然相位本身不可直接测量但光场在空间域和频率域之间存在确定的变换关系。对同一个光场施加傅里叶变换或衍射传播前后两个平面上的振幅都必须与各自的实验测量值一致。这个“双约束”构成了交替投影的数学基础也是GS算法在相位恢复领域屹立几十年的根本原因。2.2 GS算法的交替投影骨架前向传播、振幅替换、逆向回传用一句话概括GS反复在空间域约束和频率域约束之间来回投影每一次投影都保留当前估计的相位只把振幅替换成测量值。设输入面光强为 $I_1(x)$频谱面光强为 $I_2(u)$定义振幅约束 $A_1 \sqrt{I_1}$$A_2 \sqrt{I_2}$。第 $k$ 次迭代的估计光场为 $g_k(x) A_1 \cdot e^{j\varphi_k(x)}$一次完整迭代可以写成前向传播$G_k(u) \mathcal{F}{g_k(x)}$频谱域约束$\tilde{G}_k(u) A_2(u) \cdot e^{j\arg[G_k(u)]}$逆向回传$\tilde{g}_k(x) \mathcal{F}^{-1}{\tilde{G}_k(u)}$空间域约束$g_{k1}(x) A_1(x) \cdot e^{j\arg[\tilde{g}_k(x)]}$如果传播模型不是傅里叶变换而是菲涅尔衍射那就把 $\mathcal{F}$ 换成对应的衍射算子。伪代码写出来就是下面这个骨架import numpy as np def gs_skeleton(amp_input, amp_target, propagator, inv_propagator, iters100): GS迭代骨架propagator 和 inv_propagator 由具体衍射模型给出 # 初始相位随机生成振幅用输入面测量值 phase np.exp(1j * np.random.rand(*amp_input.shape) * 2 * np.pi) field amp_input * phase for _ in range(iters): # 前向传播到目标平面 field_fwd propagator(field) # 用目标振幅替换保留传播后的相位 field_fwd amp_target * np.exp(1j * np.angle(field_fwd)) # 逆向回传到输入平面 field_back inv_propagator(field_fwd) # 用输入振幅替换保留回传后的相位 field amp_input * np.exp(1j * np.angle(field_back)) return field逻辑说明前向传播把当前估计的光场送到目标平面在目标平面用实测振幅强修相位逆向回传把这个修正后的光场带回输入平面再被输入振幅约束一次。如此反复相位信息在两次约束之间被“挤”出来。参数说明amp_input和amp_target是两次测量的光强开根号propagator和inv_propagator必须是互逆的传播算子否则误差会越迭代越大iters的取值直接决定了收敛质量后文会专门讨论到底该迭代多少次。2.3 为什么不用单次傅里叶变换直接反演初学者最容易踩的坑是既然频谱面振幅已知直接对 $\sqrt{I_2} \cdot e^{j\psi}$ 做逆变换不就能拿到物面相位了吗问题是 $\psi$ 是未知的——频谱面的相位和物面相位一样丢得干干净净。单帧强度测量本质上是欠定问题不加约束根本无解。GS的做法是用两个平面的振幅约束互相“逼”出相位靠反复迭代逼近可行解。明白这一点你就知道为什么迭代次数、初始相位选择、约束窗口这些参数会影响最终结果而不是简单跑一遍就能拿到正确答案。3. 用Python搭建最小仿真平台物理引擎、GS求解器与第一个重建结果3.1 物理引擎选型为什么用角谱法而不是单次傅里叶变换仿真平台的地基是衍射传播模型。GS算法只提供了迭代框架真正决定重建精度的是前向传播算得准不准。常见的传播模型有三种傅里叶变换、菲涅尔衍射积分和角谱法。对于仿真平台来说我推荐角谱法作为默认物理引擎理由有三条第一角谱法基于瑞利-索末菲衍射积分在近场和远场都适用不需要满足菲涅尔近似条件。菲涅尔近似在传播距离极短或孔径极小的情况下会引入额外相位误差而波前重建恰恰经常面对这类场景。第二角谱法的实现只需要两次傅里叶变换和一次频域滤波计算效率高且与GS迭代天然兼容。第三它的逆传播实现简单——把传递函数共轭或者把传播距离取负即可保证前向与逆向严格互逆。菲涅尔衍射可以用一次傅里叶变换实现速度快内存省但适用条件限制在 $z^3 \gg \frac{\pi}{4\lambda}(x^2y^2)^2$ 附近一旦偏离条件就需要靠零填充凑近似。角谱法没有这个限制代价是必须保证采样间隔和传播距离避开频域混叠。在仿真平台里默认用角谱法遇到极大传播距离的场景再切换到菲涅尔。3.2 角谱法实现核心函数与参数说明角谱法的理论表达式是$$U_2(x,y) \mathcal{F}^{-1}\left{\mathcal{F}{U_1(x,y)} \cdot H(f_x, f_y)\right}$$其中传递函数为$$H(f_x, f_y) \exp\left(j\frac{2\pi z}{\lambda}\sqrt{1-\lambda^2 f_x^2-\lambda^2 f_y^2}\right)$$对应的Python实现如下这段代码就是整个仿真平台物理引擎的核心import numpy as np from numpy.fft import fft2, ifft2, fftfreq def angular_spectrum_propagate(u0, wavelength, pixel_size, z): 角谱法衍射传播 u0: 输入复振幅形状 (ny, nx) wavelength: 波长单位与 pixel_size 一致 pixel_size: 采样间隔像素尺寸 z: 传播距离正值表示前向传播 ny, nx u0.shape fx fftfreq(nx, dpixel_size) fy fftfreq(ny, dpixel_size) FX, FY np.meshgrid(fx, fy) # 频域平方根里的参数避免复数溢出 tmp 1.0 - (wavelength * FX) ** 2 - (wavelength * FY) ** 2 tmp[tmp 0] 0 # 倏逝波分量直接置零 H np.exp(1j * 2 * np.pi * z / wavelength * np.sqrt(tmp)) # 限制最大频率分量防止数值溢出 H[np.sqrt((wavelength * FX) ** 2 (wavelength * FY) ** 2) 1] 0 U0 fft2(u0) U1 U0 * H u1 ifft2(U1) return u1逻辑说明fftfreq生成与数组尺寸一致的频域坐标meshgrid构造成二维频域网格tmp中负数对应倏逝波成分物理上在传播中快速衰减置零处理是工程惯例传递函数乘以输入频谱后做逆变换得到传播后的复振幅。参数说明pixel_size是采样间隔单位建议全程统一用微米或纳米不要混用z的正负号直接控制传播方向负值即为逆向传播。注意当(wavelength * FX)^2 (wavelength * FY)^2接近 1 时sqrt(tmp)接近零传递函数进入高度振荡区此时需要检查采样是否满足奈奎斯特条件否则会出现明显的栅格伪影。3.3 跑通第一个GS相位重建从随机相位板到恢复结果物理引擎就位后我们来做一个完整的仿真实验生成一个随机相位板作为真值让光场传播一段距离在输出平面记录强度然后丢弃相位信息用GS迭代从强度中重建原始相位。import numpy as np from numpy.fft import fft2, ifft2, fftfreq # 仿真参数 np.random.seed(42) wavelength 0.5e-6 # 500nm 波长 pixel_size 3.45e-6 # 像元尺寸 3.45um size 256 # 网格 256x256 z 0.05 # 传播距离 5cm # 生成随机相位板作为真值振幅设为均匀分布 x np.linspace(-1, 1, size) X, Y np.meshgrid(x, x) true_phase 2.0 * np.exp(-(X ** 2 Y ** 2) / 0.3) 0.5 * np.random.randn(size, size) input_field np.exp(1j * true_phase) # 前向传播得到输出平面强度 output_field angular_spectrum_propagate(input_field, wavelength, pixel_size, z) target_amplitude np.abs(output_field) # 用GS迭代重建相位 iterations 200 est_phase np.exp(1j * np.random.rand(size, size) * 2 * np.pi) for i in range(iterations): # 正向传播当前估计 fwd angular_spectrum_propagate(np.exp(1j * np.angle(est_phase)), wavelength, pixel_size, z) # 替换振幅 fwd target_amplitude * np.exp(1j * np.angle(fwd)) # 逆向传播 back angular_spectrum_propagate(fwd, wavelength, pixel_size, -z) est_phase np.exp(1j * np.angle(back)) recovered_phase np.angle(est_phase)逻辑说明这个仿真实验完全模拟了真实场景——输入是已知的相位物体经过衍射后探测器只能看到强度GS迭代不断修正相位估计使传播后的强度逼近目标强度。这里振幅始终设为1聚焦纯相位重建。参数说明np.linspace(-1, 1, size)生成归一化坐标方便控制相位板的尺度true_phase里的高斯项模拟平滑相位变化随机项模拟细节这样的组合能同时考验算法的低频和高频重建能力iterations 200是经验值实际运行中建议记录每轮的误差曲线来判断是否收满。跑完这段代码你可以用np.corrcoef(true_phase, recovered_phase)计算相关系数通常能达到0.99以上。相关系数低于0.95时优先检查传播距离是否超出角谱法的适用范围或者初始相位随机种子是否导致落入局部极小值。3.4 平台目录结构与配置入口当代码量超过一个脚本能承载的范围时建议把它整理成标准的仿真平台目录结构。我的习惯是分四层物理引擎、求解器、度量函数、实验脚本。phase_retrieval_sim/ ├── engines/ │ └── propagation.py # 角谱法、菲涅尔、傅里叶传播 ├── solvers/ │ └── gs.py # GS迭代核心 ├── metrics/ │ └── quality.py # RMSE、相关系数、能量守恒 ├── experiments/ │ ├── run_basic_gs.py # 基本GS重建实验 │ └── run_param_scan.py # 参数扫描实验 └── config.py # 全局参数配置逻辑说明把物理引擎、求解器和实验脚本分离是为了让你在换算法比如把GS换成梯度下降类相位恢复时不需要动物理引擎换物理模型时不需要动求解器。config.py统一管理波长、像素尺寸、网格大小等参数避免散落在各个脚本里。参数说明这个结构中engines和solvers属于复用的核心层experiments是每次研究新建的脚本入口。一个实用技巧是把每组实验的参数和结果自动保存为JSON文件名带时间戳这样批量实验后复盘时能快速定位是哪组参数产出了哪张图。4. 仿真平台的参数设计与收敛判据从玄学到可复现4.1 四个必调参数迭代次数、过采样率、初始相位、约束窗口GS算法看起来只有两个振幅约束和一个传播算子实际调起来至少有四个参数会显著影响重建质量。第一个是迭代次数。常见做法是先用100次迭代快速看趋势再用300到500次收满但真正稳妥的方式是看收敛曲线——当误差在连续20次迭代内变化小于千分之一时可以认为迭代已经到平台期。盲目增大迭代次数不仅浪费算力还可能让算法在局部极小值附近反复震荡。第二个是过采样率这是仿真平台和实验场景差别最大的一个参数。当目标平面的采样不足以容纳衍射光场的完整带宽时频谱会混叠导致重建相位出现周期性条纹。工程上一般要求采样满足 $pixel_size \le \frac{\lambda z}{N \cdot pixel_size}$换算成过采样率通常取2到4。仿真时如果发现高频条纹优先提高网格分辨率或增大零填充而不是去调算法。第三个是初始相位。GS对初始相位敏感是出了名的随机相位初始化会在前几十次迭代产生较大的相位跳变。如果做多次重复实验我习惯固定随机种子或者在已知问题结构时用二次相位作为初始猜测这个习惯能显著减少结果方差。第四个是约束窗口。理论上的振幅约束作用于整个平面实际计算中目标平面往往只在一个有限区域内能量集中区域外是噪声。这时需要在振幅替换步骤加一个掩膜只在信号区域施加约束区域外允许自由传播。掩膜的半径通常取目标强度的能量占比阈值——我一般取总能量99%对应的半径效果比固定值稳妥得多。下面这段代码展示了如何把上述四个参数统一管理起来class GSConfig: GS仿真参数容器统一管理迭代次数、过采样率、随机种子和掩膜 def __init__(self, iterations300, oversample2.0, seed42, mask_energy0.99): self.iterations iterations self.oversample oversample self.seed seed self.mask_energy mask_energy def apply(self, amplitude): 根据能量阈值生成约束掩膜 threshold self.mask_energy * np.sum(amplitude ** 2) cumulative np.cumsum(np.sort(amplitude ** 2)[::-1]) cutoff_idx np.searchsorted(cumulative, threshold) # 简化实现实际可按半径做圆形掩膜 return amplitude np.sort(amplitude ** 2)[::-1][cutoff_idx]参数说明iterations建议根据收敛曲线动态判断初始给300是安全实践oversample用于校验采样是否满足带宽要求不直接参与计算seed保证实验可复现mask_energy控制约束区域大小。掩膜的逻辑在实验调试中经常被忽略但恰恰是影响真实数据重建质量的第一因素。4.2 收敛判据与停止条件如何确定需要迭代多少次收敛判据直接决定了“迭代多少次”这个问题。常见的做法是看输出平面振幅误差也叫振幅RMSEdef amplitude_rmse(field, target_amplitude): 输出平面振幅均方根误差用于GS收敛监控 estimated np.abs(field) return np.sqrt(np.mean((estimated - target_amplitude) ** 2))逻辑说明每一轮迭代结束后把当前估计光场传播到输出平面计算振幅与目标振幅的RMSE。这个指标理论上单调下降实际运行中会在小幅波动后进入平台期。平台期的判断标准是连续 $K$ 次迭代的RMSE相对变化小于 $10^{-3}$此时可以停止迭代。参数说明$K$ 一般取10到20。如果RMSE平台值偏高说明算法陷入局部极小此时继续增加迭代次数无济于事需要修改初始相位或放宽约束窗口。我见过太多人把迭代次数从300加到3000误差曲线纹丝不动——那不是迭代不够是参数配错了。相关系数比RMSE更适合评价相位重建质量def phase_correlation(true_phase, recovered_phase): 相位重建的相关系数接近1表示形态一致 p1 true_phase - np.mean(true_phase) p2 recovered_phase - np.mean(recovered_phase) covariance np.sum(p1 * p2) variance np.sqrt(np.sum(p1 ** 2) * np.sum(p2 ** 2)) return covariance / (variance 1e-12)逻辑说明相关系数对整体相位梯度不敏感非常适合评价“形状对不对”RMSE则对全局偏移敏感。实际评估时我两个都看——RMSE管收敛相关系数管最终质量。4.3 从纯仿真走向仿真-实验一致性噪声与像元尺寸建模仿真平台做到这一步还差一块拼图真实探测器有噪声、有像元量化、有填充因子损失。如果仿真全程用无噪声数据参数调到最优一换真实数据就失效这种“仿真-实验断崖”在相位恢复里非常普遍。给平台加噪声模型最简单有效的方式是在目标强度上叠加泊松噪声模拟光子的量子特性def add_poisson_noise(intensity, photon_count1e6): 给强度图叠加泊松噪声模拟真实探测器光子统计 norm intensity / np.sum(intensity) * photon_count noisy np.random.poisson(norm).astype(np.float64) return noisy / photon_count * np.sum(intensity)逻辑说明photon_count控制信噪比数量越小噪声越大对应低光照条件。加入噪声后再开根号作为目标振幅你会发现GS的收敛速度明显变慢RMSE平台值抬高相关系数下降——这是正常的。真实数据的相位恢复误差主要来自噪声而不是算法本身。参数说明photon_count从 $10^7$ 往下降每降一个数量级观察重建质量的退化趋势可以快速确定该算法能容忍的最低光子数。这个指标对你判断实验系统是否可行很有参考价值。除了噪声像元尺寸和填充因子影响的是采样带宽。仿真平台里只需要在传播后做一个降采样操作模拟探测器像元对连续光场的积分平均。加上这一步之后平台就能更真实地复现实验结果——从这个版本开始你在仿真里调的参数到实验台上大概率还能用。5. 波前重建仿真的四个翻车现场与排查清单5.1 迭代几十轮RMSE纹丝不动初始相位掉进了局部极小现象运行GS迭代绘图显示RMSE从前几轮开始就稳定在0.9附近没有任何下降趋势。把迭代次数从100加到1000曲线依然平坦。原因这是GS算法最经典的困境——交替投影在非凸问题上可能收敛到局部极小值也就是满足振幅约束但相位完全错误的状态。随机初始相位仍然是最大嫌疑其次是约束窗口和真实相位结构之间的不匹配。解决第一步固定随机种子连续跑三次确认不是单次随机性导致第二步改变初始相位策略改用二次相位或高斯相位作为起点第三步是缩小约束区域限制在信号支撑域内替换振幅能有效减少解空间的大小。如果三种方法都无效考虑改用多起始点策略从5个不同初始相位并行跑取相关系数最高的结果。5.2 重建相位自带周期性条纹采样频率不满足带宽约束现象重建出的相位图里出现规则排列的明暗条纹间距和像素尺寸呈固定比例改变迭代次数条纹不消失。原因输出平面的采样间隔过大导致衍射场的边缘频率折叠回来形成混叠伪影。这个坑在角谱法里尤其隐蔽因为传递函数的幅值仍然正确但不满足奈奎斯特条件的频谱被错误映射到了低频位置。解决检查采样条件 $\Delta x \le \frac{\lambda z}{N \Delta x}$如果左边大于右边就要减少物理尺寸或增大网格数也可以在被测物体周围补零扩大网格。补零是最省事的办法但会成倍增加内存占用网格从256补到512矩阵乘法耗时大约翻4倍。5.3 同一组参数换一张输入图就全面失效相位动态范围过载现象第一个测试相位板重建完美相关系数0.98换成相位动态范围更大的物体同样的参数跑出来直接是噪声。原因相位动态范围增大等价于频谱面高频分量变强原采样和约束窗口参数不再满足带宽条件。这不是算法不稳定而是仿真条件已经超出了当前配置的有效范围。解决在平台的参数配置里增加一个自动检查逻辑——根据输入的相位梯度估算最大空间频率再换算成最小需要的网格数和采样间隔。跑新图之前先过一遍检查而不是等结果出来再拍脑袋调参。这也是仿真平台相比单脚本的最大优势把经验固化成自动检查规则。5.4 仿真结果和实验数据对不上噪声模型与传播距离失配现象同一套代码仿真重建效果很好换实验强度图后相位一片模糊细节完全丢失。原因三个常见失配点。第一实验的噪声不是纯泊松还有读出噪声和暗电流低光强区域信噪比远低于仿真第二实验的传播距离标定有误差角谱法对距离参数很敏感差1毫米都可能在重建相位里引入二次相位项第三真实的物面振幅不是均匀的而仿真假设了纯相位物体。解决先在平台里用仿真强度图加不同等级的泊松和均匀噪声测试找出算法失效的噪声下限判断实验数据是否落在算法能力范围内。距离失配可以用自动对焦算法估计传播距离——在GS迭代外层包一个距离搜索每次迭代后计算重建相位的锐度指标如拉普拉斯能量取锐度最大时对应的距离。最后把物面振幅约束从均匀改为加权避免振幅假设和实验事实冲突。6. 让平台真正顶用批量参数扫描、收敛曲线横轴对齐与基线对比单次重建跑通只是开始。仿真平台真正的价值在于批量实验——同一组数据横向对比GS、梯度下降类算法和不同的参数组合。批量扫描最常见的问题是不同算法收敛速度差异巨大GS可能50轮就收敛另一个算法要200轮直接把横轴画成迭代次数曲线长度不一样没法直观对比。我的做法是先把横轴统一成“物理意义等价的迭代单位”比如都用“前向传播次数”作为单位。GS一轮迭代只做一次前向传播那它的横轴就是迭代次数本身其他算法一轮迭代可能要做两次前向传播计算梯度横轴就按两倍折算。折算之后所有算法的收敛曲线横轴数量级一致对比谁快谁慢才有意义。import json, os from datetime import datetime def run_param_scan(scan_space, base_config): 对参数网格做批量扫描保存收敛曲线与最终指标 results {} for name, params in scan_space.items(): config base_config.copy() config.update(params) # 实际在这里调用 GS 求解器 # 返回每次迭代的 RMSE、相位相关系数和耗时 rmse_curve, corr, elapsed run_single_gs(config) results[name] { params: config, rmse_curve: rmse_curve, correlation: corr, iterations: len(rmse_curve), elapsed_sec: elapsed, timestamp: datetime.now().isoformat() } os.makedirs(scan_results, exist_okTrue) with open(scan_results/scan_%s.json % datetime.now().strftime(%Y%m%d_%H%M%S), w) as f: json.dump(results, f, indent2) return results逻辑说明函数的职责很纯粹——把每次实验的迭代次数、RMSE曲线、相关系数和耗时原样存成JSON不预设好坏。跑完批量实验后你用脚本提取所有曲线的公共横轴范围对超过范围的曲线做截断或插值就能在同一个图上画出数量级一致的对比曲线。参数说明scan_space是一个嵌套字典外层是实验名内层是该实验的参数覆盖项run_single_gs返回的rmse_curve是完整列表而不是最终值这是为了方便后续画图。这个扫描框架不限于GS任何求解器只要遵循同样的输入输出约定就能在同一个平台里做公平对比。在算法对比这件事上我个人的习惯是永远保留一个“原生GS”作为基线因为论文里所有改进算法的优势都是相对传统GS而言的。批量扫描结果出来之后先看基线在哪个参数点表现最好再看新算法是否在同样的参数点上更好——如果新算法需要在不同参数下才能胜出这个对比的说服力就会打折扣。仿真平台做到这一步基本就把“从算法到结论”的链条闭环了。我自己做这类实验时的习惯是参数先小范围跑通再上批量扫描批量结果出来先看异常再看均值最终所有曲线、配置和原始结果一并归档方便几周后复盘。希望这个搭建思路能帮你少走一些弯路。本文还有配套的精品资源点击获取