简介这份基于Gerchberg-Saxton迭代算法的光学相位恢复与波前重建仿真平台面向光学工程、物理及相关专业的学习者与研究人员适用于课程实践、算法验证及成像系统模拟等场景。压缩包共30个文件约53.55MB以Matlab脚本.m、图形界面.fig、图像结果.jpg及二进制数据.bin为主另含README说明文档便于快速定位核心代码与实验数据。实现内容涵盖GS迭代算法、傅里叶变换FT2Dc/IFT2Dc、单缝衍射、光栅衍射及高斯光模拟等典型模块用户可通过调整输入参数或算法结构观察不同条件下的波前恢复效果。资源还包含备份文件与附加说明进一步支撑算法原理讲解和操作指引。已有63人学习使用适合希望掌握相位恢复基本原理、提升光学仿真与数据分析能力的学习者深入研读。1. 光学相位恢复为什么绕不开 Gerchberg-Saxton 迭代算法要做波前重建第一道坎就是光学相位恢复。CCD 和 CMOS 探测器只能记录光强光场的相位分布几乎全部丢失Gerchberg-Saxton 迭代算法常缩写为 GS用一种很朴素的手法把相位找回来在物体平面与焦平面之间反复做正反傅里叶变换每变换一个来回就用测量强度替换模值用先验约束校正另一个平面上的场。过去几十年它从相干衍射成像一路用到自适应光学到现在依然是相位恢复问题里成本最低、最好上手的基线方案。这篇内容从数学原理、Python 实现、参数整定到常见翻车点最后落到一套能改参数、能验收的仿真平台上。2. GS 迭代算法的数学骨架强度约束与物面约束的交替投影写代码之前先把 GS 算法到底在解什么说清楚。假设一束相干光场的复振幅是 f(x,y)经过一个 4f 系统或单透镜后在焦平面上的复振幅分布 U(u,v) 与 f(x,y) 构成二维傅里叶变换关系。探测器记录的强度图就是 I(u,v)|U(u,v)|²这里相位 arg[U] 完全丢掉而波前重建的目标恰恰是把这一项找回来。直接求逆是行不通的因为测量信息量少于未知量一个 N×N 的复场有 2N² 个实数自由度强度测量只给了 N² 个约束。GS 把已知信息切成两块一块是傅里叶域的强度约束另一块是物面域的先验约束然后交替投影让解在两个约束集之间反复逼近。比起后面发展出来的 Wirtinger Flow、半正定规划等算法GS 的单次迭代只是两次 FFT复杂度 O(N²logN)在大规模波前重建里仍然是最划算的起点。这里要泼一盆冷水相位恢复问题的可行域不是凸集GS 的交替投影没有严格的全局收敛保证。实际表现里它经常收敛到一个局部极小解也就是所谓“看起来每条曲线都平坦了但恢复出来的相位翻车了”。所以下面每一章里的参数整定和验收方式本质上都是在跟“局部极小”这个老对手周旋。2.1 傅里叶域的强度约束CCD 只记录模值相位信息丢在哪一步以最常用的傅里叶变换光路为例。探测器放在透镜后焦平面入射波前 f(x,y) 经过透镜后的衍射场在傍轴近似下就是它的二维傅里叶变换。这里有一个物理尺度换算必须记牢焦平面上的空间坐标 x_f 与空间频率 u 的关系是 x_f λ f u其中 λ 是波长f 是焦距。也就是说CCD 的每个像素不光是一个空间位置它同时对应一个空间频率采样点。有了这个对应关系仿真平台里的“强度约束”就非常直接给定实测强度数组 I_meas先把每个像素开方得到振幅图 A_meas(u,v)√I_meas这是唯一能从探测器里拿到的硬信息。GS 在傅里叶域做的事就是让当前猜测的复振幅 G_k(u,v) 的模值等于 A_meas同时保留 G_k 的相位。换句话说傅里叶域投影算子可以写成G_k(u,v) A_meas(u,v) · G_k(u,v) / |G_k(u,v)|当 |G_k| 接近零的时候这个除法会放大数值噪声这是后面第 4 章里要重点排查的坑之一。为什么要保留当前相位而不是随便给一个因为相位是自由度算法只能通过反复修正来让它逼近真实解如果每次迭代都把相位随机重抽信息就永远积累不起来。保留相位、替换模值本质上是把测量强度当作“锚点”让相位在迭代中慢慢长出来。2.2 物平面的支撑与实值约束一次完整迭代的四步与算法流程物面域的约束取决于具体应用。最常见的两类先验一类是支撑约束即知道物体或波前只在某个孔径范围内有值支撑区外的直接置零另一类是纯相位约束即已知物体幅度均匀或接近均匀把幅度归一化到 1只保留相位。仿真平台里这两类常组合使用也可以单独用。一次 GS 迭代拆成四步步骤操作使用的先验1当前物面场 g 做正傅里叶变换得到 G传播模型2保留 G 的相位用实测振幅替换模值测量强度3对修正后的 G 做逆傅里叶变换回到物面传播模型4对物面场施加支撑或纯相位约束物面先验第 1 步到第 3 步解决的是“测量一致性”第 4 步解决的是“物理合理性”。整个算法就是在这两个约束集合之间来回踢皮球直到两边都不再明显变化。初始猜测的选择是个容易被低估的变量。物面初始相位通常用均匀随机相位这一步没有太严格的理论依据更像是一场赌博。如果初始猜测离真实解太远GS 可能从第一步开始就走进了局部极小多跑几组不同的随机初始猜测选傅里叶域残差最小的那个结果是成本很低的后悔药。这一习惯在许多开源相位恢复代码里都有体现仿真平台上尤其值得坚持。3. 搭建 GS 光学相位恢复仿真平台Python 最小实现与参数整定这章直接给一个能跑的最小平台。代码不依赖光学仿真库只用 numpy 的 FFT 和 matplotlib 画图。先搭出“真实波前 → 模拟 CCD 强度 → GS 恢复 → 对比真值”的完整链路后续所有参数实验都在这套骨架上改。3.1 最小可复现代码numpy 的正反 FFT 循环与关键参数import numpy as np import matplotlib.pyplot as plt # ---------- 参数区 ---------- N 256 # 网格尺寸2 的幂最省事 pix 3.45e-6 # CCD 像元尺寸单位 m wl 632.8e-9 # 波长单位 m f_lens 50e-3 # 傅里叶透镜焦距单位 m iters 300 # GS 迭代次数 # 物体平面等效采样间隔由焦距和探测器共同决定 dx wl * f_lens / (N * pix) # 物体平面像素间距 x np.linspace(-N * dx / 2, N * dx / 2, N) XX, YY np.meshgrid(x, x) # ---------- 生成一个带像差的已知测试波前 ---------- r_normal N * dx / 2 # 归一化半径 true_phase 1.2 * (XX / r_normal) ** 2 0.8 * (YY / r_normal) ** 2 true_field np.exp(1j * true_phase) # 纯相位物体幅度为 1 # ---------- 模拟 CCD 测量正变换取模加噪声 ---------- F_meas np.fft.fftshift(np.fft.fft2(true_field)) intensity np.abs(F_meas) ** 2 noise np.random.normal(0, 0.01 * intensity.max(), intensity.shape) intensity_noisy np.abs(intensity noise) amplitude_meas np.sqrt(intensity_noisy) # 探测器只能拿到振幅 # ---------- GS 主循环 ---------- rng np.random.default_rng(42) guess np.exp(1j * rng.uniform(-np.pi, np.pi, (N, N))) # 随机初始相位 err_history [] for k in range(iters): # 1. 正傅里叶变换到焦平面 G_est np.fft.fftshift(np.fft.fft2(guess)) # 2. 保留相位替换振幅 phase G_est / (np.abs(G_est) 1e-12) G_corr amplitude_meas * phase # 3. 逆傅里叶变换回物面 field_est np.fft.ifft2(np.fft.ifftshift(G_corr)) # 4. 物面约束这里假设纯相位物体归一化幅度 field_est field_est / (np.abs(field_est) 1e-12) guess field_est # 记录傅里叶域残差 err np.mean((np.abs(G_est) - amplitude_meas) ** 2) / np.mean(amplitude_meas ** 2) err_history.append(err) # ---------- 输出对比 ---------- recovered_phase np.angle(guess) plt.figure(figsize(12, 4)) plt.subplot(1, 3, 1) plt.imshow(true_phase, cmapjet) plt.title(true phase) plt.subplot(1, 3, 2) plt.imshow(recovered_phase, cmapjet) plt.title(recovered phase) plt.subplot(1, 3, 3) plt.plot(err_history) plt.yscale(log) plt.title(fourier residual) plt.show()这段代码的逻辑分成三段。第一段是物理参数换算物体平面的采样间隔 dx 不是随便取的必须满足 dx λf/(Npix) 这个关系才能保证 FFT 结果和真实焦平面上的强度分布一一对应。很多刚接触 GS 的人直接在物体域随便给个网格最后恢复出的相位跟真实测量怎么都对不上多半是这里出了问题。第二段是模拟测量。因为仿真平台里“真实波前”是我们自己造出来的才能事后做定量验收真实实验里没有这一步直接把 CCD 强度图喂给循环。这里人为加了高斯噪声模拟相机本底实际散粒噪声占主导时建议改成 np.random.poisson 再开方。第三段是 GS 主循环四个步骤跟第 2 章完全对应。最容易被忽略的是 phase 那行里的 1e-12当某个频率上实测振幅几乎为零时G_est 的模也很小直接除会得到巨大的补偿值导致整个迭代发散。加一个极小量是常规操作。最后记住把收敛历史 err_history 打出来看曲线不要只看恢复图。3.2 网格规模、波长焦距、支撑宽度与噪声水平四个参数怎么选网格规模 N 是最直接的选择。128 以下跑得飞快但相位图的细节丢失严重1024 以上 FFT 的内存压力和耗时成倍增加且对恢复精度提升有限。我的经验是测试阶段用 256正式验收换 512除非你关心的是高频像差否则没有必要上 2048。波长和焦距在迭代里不参与数学计算它们只负责把物理坐标与频率坐标对应起来。改 λ 或 f 不改变 GS 的收敛行为但会改变物体平面采样范围。比如波长从 632.8nm 换成 532nmdx 变小同样的 N 下物体视场变窄。若测试波前尺寸不变可能出现视场装不下的问题这一点在对比不同激光器仿真结果时经常踩坑。支撑宽度要区分两种情况。如果物体天然就在孔径内支撑约束选真实孔径的 1.01.2 倍比较稳妥太紧会切断真实信号的高频边缘恢复结果出现振铃假象太宽则约束形同虚设。如果仿真只针对纯相位波前第三步的幅度归一化本身就是最强约束不需要额外支撑。噪声水平的控制没有万能公式可以先在无噪声情况下把基线跑通再逐步加大噪声观察残差曲线的抬升幅度确定你算法能接受的信噪比下限。提示常见做法是先在 64×64 小网格上调通行为再放大到 256×256。小网格上每步迭代只要几毫秒可以快速试探初始猜测、噪声和支撑三者的组合省下的时间远比看着大图跑几百次迭代划算。4. 波前重建仿真里的 5 个常见翻车现场现象、原因与避坑排查GS 的代码短短几十行真正难的是调试。下面这 5 个问题是仿真平台上最高频的翻车点每条都按“现象 → 原因 → 解决”来写可以直接对照排查。4.1 恢复相位出现大块错乱条纹收敛曲线却已经压平现象迭代到 100 次以后傅里叶域残差几乎不再下降但恢复相位图里仍有大面积的条纹错乱和真实相位对不上。原因GS 掉进了局部极小。相位恢复问题非凸没有理论保证收敛到全局最优。傅里叶域残差平坦只能说明“在某个解附近稳定”不能说明解是对的。解决多初始猜测是最直接的手段。用 2050 组不同的随机相位初始化分别跑同一套迭代记录每组结束时的残差选残差最小的一次当作最终结果。残差相同时再参考物面支撑区外的能量支撑外残余能量越少通常越好。4.2 恢复相位整体偏移明明像差形状对但数值整体不对现象形状和真实相位一致但整体差了一个常数相位或者沿横向有个平移导致定量比较时误差很大。原因GS 解的固有歧义。焦平面强度对整体常数相位不敏感也就是 exp(iθ₀) 乘在整个场上不会改变强度另外物体平面平移等价于傅里叶域乘线性相位因子如果不做配准恢复场会出现平移。解决做量化评估之前必须先做全局相位补偿。用复振幅对齐的方式计算一个全局缩放因子 c把恢复场乘以 c 的共轭后再和真值比较。具体代码在第 5 章给出这一步不做后面的波前 RMS 全是虚的。4.3 迭代次数增加后误差反而上升曲线出现锯齿现象前 50 次迭代残差下降正常50 次之后突然跳高之后反复震荡甚至彻底发散。原因模值替换步骤里除以接近零的振幅导致数值爆掉。尤其是激光照明边缘的暗区或探测器上的坏像素那里强度接近零开方后的振幅很小除出来的相位权重被噪声主导一次迭代就能污染整个物面场。解决给振幅加一个底噪上限比如把低于全局峰值 1% 的振幅强行抬升到 1% 峰值或者直接加一个 Tikhonov 形式的正则项。仿真里最省事的做法就是代码里的 1e-12但真实数据里 1e-12 往往不够要用相对底噪来设置。4.4 恢复结果边缘出现环状涟漪中间区域尚可现象视场中心恢复得很好边缘一圈却出现类似衍射环的同心波纹且波纹随迭代次数增强。原因周期性混叠加上支撑边界效应。FFT 默认把输入当作周期信号处理如果真实波前在视场边缘没有衰减到零周期性延拓就会在边界产生伪边缘支撑约束的外沿又把能量反推回中心形成环状假象。解决让真实测试波前在采样窗口内占半幅以上但边缘留出渐晕区让幅度自然衰减到零。更严格的方法是在物面加一个超高斯窗例如乘上 exp(-(r/r₀)⁴)让波前在边界处平滑过渡。恢复结果评估时只取中间 70% 区域避开边界效应。4.5 改了一个参数之后之前能收敛的配置全都不收敛现象波长、焦距或像素尺寸换了一组其他代码完全没动原来很好的恢复结果突然变差。原因物理尺度换算脱节。改 λ 或 f 之后物体平面的采样间隔 dx 跟着变如果真实波前是用像素坐标定义的那它对应的物理尺寸就变了重新生成测试波前时像差的归一化半径必须按新的 dx 重新计算。解决所有测试波前都用物理坐标生成不要直接写“第 100 个像素开始”。我在实际项目中踩过一次把波长从 633nm 换到 405nm 后忘记改归一化半径离焦项整体缩小了 1.5 倍排查了两天才发现是坐标尺度问题。物理单位贯穿到底能省掉这类隐性坑。5. 恢复结果怎么验收定量指标、补偿技巧与整套验收流程相位恢复仿真的验收不能只靠肉眼“看着像”。人眼对相位的判断很不敏感尤其当真实相位和恢复结果同时被显示成彩色图时颜色相近不代表数值接近。这一章给出三件套先补偿全局相位再看两组残差最后跑一套标准流程。5.1 全局相位补偿直接盯相位真值之前先做的一次对齐由于整体常数相位不影响强度GS 恢复出来的复场与真实复场之间天然差一个未知的常数因子 e^{iθ₀}。直接对相位图做差会得到一个遍布全图的误差根本没法用。正确做法是在复振幅层面找出最优对齐因子再比较相位。# recovered 来自上一章 GS 主循环true_field 是生成测试波前时保存的字段 # 1. 计算复振幅对齐因子 c np.vdot(recovered, true_field) / np.vdot(recovered, recovered) # 2. 恢复场乘上 c 的共轭让它和真实场对齐 aligned recovered * np.conj(c) # 3. 相位差注意相位可能被包裹到 (-pi, pi] phase_diff np.angle(aligned * np.conj(true_field)) # 4. 波前 RMS 误差以波长为单位 rms_rad np.sqrt(np.mean(phase_diff ** 2)) rms_lambda rms_rad / (2 * np.pi) print(f波前 RMS 误差 {rms_lambda:.4f} lambda)np.vdot 对两个复场做内积得到的 c 是一个复数包含幅度缩放和相位对齐两层补偿。真实波前是纯相位时幅度本来就是 1c 的模接近 1主要补偿就是相位。算完 phase_diff 后建议用 np.angle 而非直接求实部因为复数乘法的相位差天然包裹实部接近 1 时角度才有意义。5.2 两个独立的收敛判据傅里叶域残差与物面约束残差多数 GS 代码只记录傅里叶域残差也就是第 3 章代码里的 err_history。这个量能反映测量一致性但一个局部极小解同样可以做到残差很低。所以还要加一个物面约束残差用来衡量支撑区外能量占比support_outside ~support_region # 布尔掩膜支撑区外为 True object_err np.mean(np.abs(guess[support_outside]) ** 2) / np.mean(np.abs(guess) ** 2)如果支撑约束是可靠先验object_err 应该随迭代逐渐下降。如果它先降后升往往是模值替换步骤在边缘引入了高频能量支撑区外开始积累假信号。两个残差一起看一个管“测量符合度”一个管“先验符合度”同时平坦才说明迭代真的停了。5.3 可复制的验收流程从生成测试波前到判定通过的三段式一套能反复使用的验收流程我通常固定为四步每次都走完再下结论。第一步指定测试波前。用 Zernike 像差生成真实相位例如离焦加球差幅度设成 PV 值 14 波长的量级。第二步用同一套传播模型生成模拟强度加上目标信噪比的噪声。第三步运行 GS 主循环记录两个残差曲线和最终波前 RMS。第四步判定是否通过波前 RMS 误差低于 λ/20 且两个残差都进入平台期判通过RMS 在 λ/10 到 λ/20 之间需要记录初始猜测组数和噪声水平RMS 高于 λ/10回到参数整定。这个阈值本身是个经验范围不是行业标准。自适应光学系统里常要求 λ/10 或更严仿真里如果连 λ/20 都达不到基本要先怀疑支撑或初始化不要急着加迭代次数。GS 到了平台期之后加再多迭代也无济于事这是血泪教训。6. 把 GS 仿真平台迁移到真实数据前要养成的三个习惯仿真平台跑通了接下来面对真实探测器数据时有三个习惯比调算法参数更管用。第一个习惯是先做强度图预处理再进 GS。真实 CCD 有暗电流、坏像素和响应不均匀直接开方成振幅喂进迭代那些坏点会在模值替换时被放大成星芒状噪声。我会先拍暗场取平均用平场校正把响应归一化再做中值滤波剔除坏点。这一步和光学无关纯粹是给算法一个干净的输入。第二个习惯是初始猜测不要只用纯随机相位。随机相位覆盖面广但收敛路径不可控更稳妥的做法是用低阶像差的先验来初始化。具体说先用一个包含离焦和像散项的解析相位作为初始猜测再在它上面叠加小幅随机扰动。这个方式在多组初始化时能更快找到低残差解代价是可能漏掉某些高频结构——所以我会各跑几组对比后再定。第三个习惯是保留多份不同条件下的测量数据验证解的唯一性。GS 恢复的问题是可能有多组不同物面对应同一组强度。真实实验中换一个焦点位置、换一个孔径形状再测一次得到的冗余测量能显著压缩解空间。仿真里模拟这种冗余很简单就是跑两个不同噪声种子下的强度看两组恢复结果是否一致若残差接近但相位明显不同说明解不唯一需要增加约束。最后说个实际事例。年前做一次镜片面形恢复我图省事只用了单次测量强度恢复结果在支撑边缘始终有条假条纹。加上辅助离焦测量之后同样的 GS 代码假条纹直接被抑制掉了。自那以后我做任何相位恢复都坚持至少留一张冗余视角数据。GS 这个活了五十多年的算法边界在哪我还在摸但先把输入和初始化做扎实往往比折腾迭代公式更快。希望帮到你。本文还有配套的精品资源点击获取
