简介：这份基于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 步到第 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, cmap='jet') plt.title('true phase') plt.subplot(1, 3, 2) plt.imshow(recovered_phase, cmap='jet') 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 换成 532nm，dx 变小，同样的 N 下物体视场变窄。若测试波前尺寸不变，可能出现视场装不下的问题，这一点在对比不同激光器仿真结果时经常踩坑。

支撑宽度要区分两种情况。如果物体天然就在孔径内，支撑约束选真实孔径的 1.0～1.2 倍比较稳妥；太紧会切断真实信号的高频边缘，恢复结果出现振铃假象，太宽则约束形同虚设。如果仿真只针对纯相位波前，第三步的幅度归一化本身就是最强约束，不需要额外支撑。噪声水平的控制没有万能公式，可以先在无噪声情况下把基线跑通，再逐步加大噪声观察残差曲线的抬升幅度，确定你算法能接受的信噪比下限。

提示：常见做法是先在 64×64 小网格上调通行为，再放大到 256×256。小网格上每步迭代只要几毫秒，可以快速试探初始猜测、噪声和支撑三者的组合，省下的时间远比看着大图跑几百次迭代划算。

4. 波前重建仿真里的 5 个常见翻车现场：现象、原因与避坑排查

GS 的代码短短几十行，真正难的是调试。下面这 5 个问题是仿真平台上最高频的翻车点，每条都按“现象 → 原因 → 解决”来写，可以直接对照排查。

4.1 恢复相位出现大块错乱条纹，收敛曲线却已经压平

现象：迭代到 100 次以后，傅里叶域残差几乎不再下降，但恢复相位图里仍有大面积的条纹错乱，和真实相位对不上。

原因：GS 掉进了局部极小。相位恢复问题非凸，没有理论保证收敛到全局最优。傅里叶域残差平坦只能说明“在某个解附近稳定”，不能说明解是对的。

解决：多初始猜测是最直接的手段。用 20～50 组不同的随机相位初始化，分别跑同一套迭代，记录每组结束时的残差，选残差最小的一次当作最终结果。残差相同时再参考物面支撑区外的能量，支撑外残余能量越少通常越好。

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 是一个复数，包含幅度缩放和相位对齐两层补偿。真实波前是纯相位时幅度本来就是 1，c 的模接近 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 值 1～4 波长的量级。第二步，用同一套传播模型生成模拟强度，加上目标信噪比的噪声。第三步，运行 GS 主循环，记录两个残差曲线和最终波前 RMS。第四步，判定是否通过：波前 RMS 误差低于 λ/20 且两个残差都进入平台期，判通过；RMS 在 λ/10 到 λ/20 之间，需要记录初始猜测组数和噪声水平；RMS 高于 λ/10，回到参数整定。

这个阈值本身是个经验范围，不是行业标准。自适应光学系统里常要求 λ/10 或更严，仿真里如果连 λ/20 都达不到，基本要先怀疑支撑或初始化，不要急着加迭代次数。GS 到了平台期之后加再多迭代也无济于事，这是血泪教训。

6. 把 GS 仿真平台迁移到真实数据前要养成的三个习惯

仿真平台跑通了，接下来面对真实探测器数据时，有三个习惯比调算法参数更管用。

第一个习惯是先做强度图预处理再进 GS。真实 CCD 有暗电流、坏像素和响应不均匀，直接开方成振幅喂进迭代，那些坏点会在模值替换时被放大成星芒状噪声。我会先拍暗场取平均，用平场校正把响应归一化，再做中值滤波剔除坏点。这一步和光学无关，纯粹是给算法一个干净的输入。

第二个习惯是初始猜测不要只用纯随机相位。随机相位覆盖面广，但收敛路径不可控；更稳妥的做法是用低阶像差的先验来初始化。具体说，先用一个包含离焦和像散项的解析相位作为初始猜测，再在它上面叠加小幅随机扰动。这个方式在多组初始化时能更快找到低残差解，代价是可能漏掉某些高频结构——所以我会各跑几组，对比后再定。

第三个习惯是保留多份不同条件下的测量数据验证解的唯一性。GS 恢复的问题是可能有多组不同物面对应同一组强度。真实实验中换一个焦点位置、换一个孔径形状再测一次，得到的冗余测量能显著压缩解空间。仿真里模拟这种冗余很简单，就是跑两个不同噪声种子下的强度，看两组恢复结果是否一致；若残差接近但相位明显不同，说明解不唯一，需要增加约束。

最后说个实际事例。年前做一次镜片面形恢复，我图省事只用了单次测量强度，恢复结果在支撑边缘始终有条假条纹。加上辅助离焦测量之后，同样的 GS 代码，假条纹直接被抑制掉了。自那以后我做任何相位恢复，都坚持至少留一张冗余视角数据。GS 这个活了五十多年的算法，边界在哪我还在摸，但先把输入和初始化做扎实，往往比折腾迭代公式更快。希望帮到你。

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