简介基于Gerchberg-Saxton迭代算法的光学仿真平台面向光学工程、信息光学等方向的学子与科研人员用于解决光强数据中相位恢复和波前重建问题。资源包共30个文件以Matlab脚本.m、图形文件.fig、图像.jpg及数据文件.bin为主整体约53.55MB涵盖理论说明、操作指引与核心仿真代码。通过调整输入参数或算法结构可模拟不同条件下的波前恢复效果适用于全息成像、显微成像和光学校正等场景。内含的多个算法脚本与示例图形能帮助使用者深入理解GS算法的迭代原理并掌握从光强分布到相位重建的完整流程。目前已有62人学习下载适合用于光学课程实践、算法验证及相关课题预研。1. 从强度测量到波前重建为什么GS迭代算法适合做仿真平台的底层核心做光学干涉测量、衍射成像或者自适应光学的人经常会碰到同一个尴尬相机只能记录光场的强度分布而相位信息在探测时被丢掉了。可是波前重建偏偏又依赖相位比如要反推一个畸变镜的面形或者再现一个全息物体的相位结构。Gerchberg-Saxton迭代算法就是用来干这件事的它利用物面与像面两个域上已知的振幅约束反复做傅里叶变换和逆变换把丢失的相位一步步“猜”回来。这个算法只有几十行核心代码不依赖昂贵的干涉光路很适合作为光学相位恢复仿真平台的数学引擎。对于想理解相位恢复原理、快速验证重建效果、或者把算法嵌入到更大波前传感系统里的工程师来说GS算法是最容易跑通的第一块跳板。需要明确的是GS算法并不是全能的。它要求知道物面振幅和像面强度并且光场传播满足傅里叶变换关系比如夫琅禾费衍射。实际仿真平台里我们可以用采样光场近似这些条件让GS算法在二维矩阵上迭代从而直观看到每一次迭代如何修正相位。2. GS迭代算法的数学基础与收敛边界2.1 傅里叶光学里的相位恢复问题建模在标量衍射理论中一束单色相干光从物面传播到像面满足傅里叶变换关系。设物面复振幅为f(x, y) |f(x, y)| * exp(i * φ(x, y))像面复振幅为F(u, v) |F(u, v)| * exp(i * ψ(u, v))。探测器能测量到的是强度I(u, v) |F(u, v)|^2相位ψ完全丢失。如果我们还提前知道物面的振幅分布比如知道入射光均匀照明那么相位恢复问题就变成在“物面振幅已知”和“像面强度已知”两个约束下求出满足傅里叶变换关系的相位分布。GS算法的核心是把重建问题转化为交替投影问题。每次投影分两步在物面域用已知振幅替换当前重建场振幅保留相位在像面域用测量振幅替换变换后的振幅保留相位。反复执行直到相位收敛。下表是仿真平台里需要明确的两个域约束。定义域已知条件每次迭代施加的约束物面振幅 |f_known|通常为常量或已知光瞳保留计算相位强制振幅等于已知值像面强度 I_measured对应振幅 sqrt(I)保留计算相位强制振幅等于测量值的平方根这种交替投影可以看成是两个集合之间的交替投影所有可能物面的集合和所有傅里叶积分后满足测量强度的像面的集合。当两个集合都是凸集时交替投影可以保证收敛但相位恢复问题的集合不是凸的这也是GS算法容易陷入局部极小、收敛停滞的原因。2.2 标准GS迭代步骤与伪代码整个GS迭代过程的核心操作只有四个傅里叶变换、像面约束、逆傅里叶变换、物面约束。写成伪代码非常直观。# 初始化复振幅使用已知物面振幅与随机相位 field amplitude_known * np.exp(1j * np.random.rand(M, N) * 2 * np.pi) for k in range(max_iter): # 1. 正向傅里叶变换到像面 spectrum np.fft.fft2(field) # 2. 像面约束用测量振幅替换计算振幅 spectrum measured_amplitude * np.exp(1j * np.angle(spectrum)) # 3. 逆傅里叶变换回物面 field np.fft.ifft2(spectrum) # 4. 物面约束用已知物面振幅替换计算振幅 field amplitude_known * np.exp(1j * np.angle(field))这里的关键参数有三个max_iter控制迭代次数amplitude_known是物面振幅先验measured_amplitude是像面测量幅度的平方根。要注意fft2得到的结果默认零频在左上角实际物理仿真里需要用fftshift处理中心化但这对GS迭代本身没有影响因为傅里叶变换的可交换性保证了交替投影依然成立。2.3 收敛边界为什么GS算法在仿真平台上需要加入松弛因子GS算法虽然简单但收敛性质并不好。两个集合的非凸性导致迭代可能在一个环里震荡无法收敛到全局最优。具体表现是均方误差曲线在若干次迭代后不再下降但重建相位仍然有空间上的条纹误差。一个实用的缓解手段是引入松弛参数beta把物面约束改成混合更新f_new (1 - beta) * f_old beta * amplitude_known * exp(i * angle(f_old))当beta取1时退化为标准GS小于1时更平滑但收敛变慢。大于1可能加速也更容易发散。在我的仿真平台里通常把beta作为外部可调参数方便对比不同对象的重建效果。另一个容易被忽视的问题是采样率。物面和像面都需要满足奈奎斯特采样否则傅里叶变换的能量分布会失真。仿真平台里建立光场网格时必须把传播距离、波长、像素尺寸对应到正确的傅里叶变换标定上否则GS算法迭代到的是一个错误的变换尺度相位恢复自然失败。3. 在仿真平台里写出可复现的GS重建核心3.1 选择仿真工具和光场数据结构常见做法是使用Python做光学仿真因为NumPy和SciPy已经提供了完整的快速傅里叶变换流程Matplotlib可以方便画出波前分布。波前重建仿真平台的核心数据结构一般是一个二维复数数组表示物面上的复振幅。相位恢复的目标就是从强度测量数据中恢复这个数组的相位部分。仿真平台的光场参数建议单独定义在一个配置对象里至少包含参数说明典型值pixel_size物面像素尺寸10 μmwavelength光波波长633 nmprop_distance传播距离0.1 mgrid_size采样网格点数256×256amplitude_known物面振幅分布圆形光瞳3.2 基于NumPy的完整GS核心函数下面是一段可直接放进仿真平台运行的GS算法实现。它支持自定义物面振幅、像面测量振幅和迭代次数并且返回重建相位。import numpy as np def gerchberg_saxton(measured_amplitude, amplitude_known, max_iter200, beta1.0, verboseFalse): 使用 Gerchberg-Saxton 迭代算法恢复相位。 参数: measured_amplitude: 像面测量振幅二维数组 amplitude_known: 物面已知振幅二维数组 max_iter: 最大迭代次数 beta: 松弛因子(0, 1] 表示阻尼更新 返回: recovered_phase: 重建相位单位为弧度 # 初始化指定振幅随机相位 phase_init 2 * np.pi * np.random.rand(*amplitude_known.shape) field amplitude_known * np.exp(1j * phase_init) for k in range(max_iter): # 正向传播到像面 spectrum np.fft.fft2(field) # 像面约束保留相位替换振幅 spectrum measured_amplitude * np.exp(1j * np.angle(spectrum)) # 逆传播回物面 field_next np.fft.ifft2(spectrum) # 物面约束阻尼更新 constrained amplitude_known * np.exp(1j * np.angle(field_next)) field beta * constrained (1 - beta) * field if verbose and k % 20 0: error np.mean(np.abs(np.abs(field) - amplitude_known) ** 2) print(fIteration {k}: amplitude error {error:.6f}) return np.angle(field)需要特别说明的是beta的写法我给constrained和旧field做加权平均在非凸情况下能有效减少迭代震荡。verbose参数方便在平台上实时观察收敛情况。实际调用时measured_amplitude来自探测器可能带有噪声这时需要在仿真里加入高斯噪声或者泊松噪声来测试算法的鲁棒性。3.3 初始化策略与误差函数初始化对GS算法影响很大。随机相位是最常见的选择但如果已知相位大致是平滑波前用一个二次曲面初始化会收敛更快。例如在自适应光学仿真平台里可以先假设一个离焦波前作为初始相位再交给GS迭代去细修。为了在仿真平台上量化重建精度通常计算重建相位与真实相位之间的均方根误差RMSE。由于光场存在全局相位偏移和整体常数相位比较之前需要把重建相位减去均值并且把像素值转到[-π, π]区间。这一步在平台实现里非常关键否则RMSE会被2π跳变直接拉高。4. 用仿真平台做波前重建实战从生成像面到恢复相位4.1 构造模拟测量数据仿真平台的价值在于可以生成已知真实相位的波前然后模拟探测器记录强度再用GS算法恢复相位最后对比恢复结果与真实相位。整个过程不需要任何真实光学器件。下面这段代码生成一个带有离焦和像散的真实相位模拟像面强度并调用上文的GS函数。from scipy.fft import fft2, ifft2 import matplotlib.pyplot as plt # 生成坐标网格 N 256 x np.linspace(-1, 1, N) X, Y np.meshgrid(x, x) # 真实相位离焦 像散 true_phase 2.0 * (X**2 Y**2) 1.5 * (X**2 - Y**2) # 物面振幅圆形光瞳 amplitude_known ((X**2 Y**2) 0.9**2).astype(float) # 物面复振幅 field_obj amplitude_known * np.exp(1j * true_phase) # 模拟传播到像面得到测量振幅 spectrum fft2(field_obj) measured_amplitude np.abs(spectrum) # 调用GS算法 recovered_phase gerchberg_saxton( measured_amplitude, amplitude_known, max_iter300, beta0.8 ) # 对比真实相位与重建相位 fig, axes plt.subplots(1, 3, figsize(12, 4)) im0 axes[0].imshow(true_phase, cmapviridis); plt.colorbar(im0, axaxes[0]). # 以下省略绘图细节需要留意的是物面振幅如果是二值光瞳约束时会把光瞳外的振幅强制设为0光瞳内设为1。这个支撑约束非常强能帮助GS算法更快收敛。像面测量振幅如果不做低通滤波GS会把探测器噪声当成真实高频信息导致重建相位出现椒盐状伪影。因此在仿真平台里我通常会在测量振幅上先做一个适度的高斯平滑再交给算法迭代。4.2 评价重建质量的三个指标波前重建仿真平台不能只看肉眼效果需要量化指标。指标计算方法合格范围RMSEsqrt(mean((reconstructed - true)^2))小于0.1 radStrehl比接近1表示衍射极限0.8相位残差峰谷值max(reconstructed - true) - min(...)与原始PV相比下降一个数量级RMSE需要在去掉整体活塞相位后计算。Strehl比可以从残余波前方差近似推导S ≈ exp(-var(残差))非常方便。如果仿真结果始终达不到这些指标大概率不是GS迭代次数不够而是物面约束或像面振幅预处理有问题。4.3 平台化的参数扫描与批量实验一个仿真平台通常要支持多个波长、多种像差、不同噪声水平的测试。我的做法是把GS核心函数包一层实验驱动函数用字典传入参数批量跑结果并自动生成收敛曲线。python run_gs_experiment.py --pixel-size 10e-6 --wavelength 633e-9 --noise 0.05命令行参数可以直接映射到平台配置把GS迭代结果写入CSV文件便于统计。这样做的好处是可以快速验证一个想法比如想知道增加迭代次数能不能抵消噪声影响跑一批实验就能得出结论。5. 让GS算法在仿真平台里更稳的几个实践技巧5.1 遇到迭代震荡时先调松弛因子标准GS在相位分布剧烈变化时容易震荡表现为误差曲线不下降。这时不要盲目增大迭代次数而是把beta调到0.5~0.9之间可能会立刻恢复单调收敛。如果调小后收敛太慢可以配合动态增大beta例如前100次迭代用0.6后面逐渐提高到1.0。5.2 用混合输入-输出算法处理支撑约束如果物面光瞳约束不可靠比如只知道一个大致的支撑区域可以用HIO算法替换第4章里的物面约束。HIO的更新公式为field_new field - beta * (field - constrained) outside_support ~support_mask field_new[outside_support] field[outside_support] - beta * field[outside_support]这比标准GS更适合处理非凸问题尤其适合扩展物体。仿真平台里常常在标准GS不稳定时自动切到HIO模式。5.3 验证重建结果的唯一性最容易被忽略的验证是把重建相位代回傅里叶变换看计算强度是否与测量强度吻合。如果吻合而相位与真实相位差异仍大说明问题可能出现在物面振幅约束太弱或者存在孪生相位解。这时需要用多平面强度测量来打破歧义或者引入预先标定的波前先验。GS算法不是一个黑盒它的每一次迭代都是在傅里叶变换和约束之间找平衡。把上面这些技巧落进仿真平台你会看到收敛速度明显变快重建结果也从“看着像”变成“数值上合格”。本文还有配套的精品资源点击获取