简介基于GS算法的光学相位恢复与波前重建仿真平台是一套面向光学工程、信息光学与计算成像方向学习者和研究者的Matlab仿真资源利用GS算法在空间域与频域间反复迭代逼近真实波前适用于全息成像、显微成像、光学校正和精密测量等场景。压缩包共30个文件总大小约53.55MB包括10个Matlab脚本、5个fig图形、5张jpg结果图、4个zbak备份、2个bin数据文件以及txt和md文档分别用于算法仿真、结果可视化、调试备份与操作说明。通过调整迭代次数、初始相位或光强分布可直观观察GS算法收敛过程及不同参数对重建精度的影响并借助说明文档快速掌握运行方式同时可对照图形与图像分析重建前后的相位分布差异。这套资源既适合光学课程实践与课程设计也能为相关科研与算法二次开发提供验证基础帮助深入理解相位恢复原理并提升仿真分析能力。目前已有63人学习资源来源于网络分享仅限学习交流使用。1. 相位恢复为什么绕不开Gerchberg-Saxton一个只能测强度的问题现在的CCD和CMOS相机只能记录光强的空间分布相位信息在探测过程中直接丢失。可光学设计、像差测量、无透镜成像和自适应光学偏偏都要用到波前数据。Gerchberg-Saxton迭代算法以下按工程习惯简称GS算法用两个不同平面的强度图交替迭代把丢失的相位一步一步“猜”回来是光学相位恢复与波前重建领域最基础也最常用的一类算法。围绕它搭一个仿真平台好处是可以不花钱搭光路就能批量验证算法在不同像差、噪声和初值条件下的表现把恢复误差量化出来。这篇笔记适合正在做计算成像、光束整形、波前传感或者自适应光学预研的工程师和学生也适合手里有实测光强图、想知道GS到底能不能用的场景。2. 交替投影的数学直觉两个振幅约束如何逼近丢失的相位2.1 相位是光学里的“黑匣子”为什么两幅强度图就能约束它用复振幅描述一束光场可以写成 U(x,y)A(x,y)exp(iφ(x,y))。探测器记录的是光强 I|U|²A²相位 φ(x,y) 在平方操作里被抹掉了。单拿一幅强度图只能得到振幅 A对相位几乎没有任何约束信息量上天然不足。但如果让这束光继续传播到另一个平面比如经过透镜后在焦面上再采集一幅光强图两幅强度图就构成了一组约束物面振幅已知、焦面振幅已知、相位未知。在夫琅禾费近似下焦面的复振幅正好是物面复振幅的二维傅里叶变换。于是问题被改写成已知傅里叶变换两端的振幅分布求两端的相位。这个约束关系是非线性的没有办法直接求逆只能靠迭代逼近。GS算法的核心思路是交替投影在两个约束集合之间反复横跳用“替换振幅、保留相位”这个单一操作让结果同时靠近两端的振幅约束。每完成一轮替换相位就在两个平面之间被“传递”一次最终收敛到一个同时满足两端振幅约束的解那个解对应的相位就是恢复出的波前。这套思想最早由 Gerchberg 和 Saxton 在 1972 年提出后来 Fienup 等人发展出带支撑域反馈的输入输出算法成为 GS 家族里最重要的变体。理解交替投影是理解后面所有改进算法的起点。2.2 一次GS循环的四步正变换、替换、逆变换、再替换GS 的每一次循环可以拆成四个动作用伪代码写出来非常紧凑E0(x,y) sqrt(I_obj) * exp(i * phase_init) # 物面初始复振幅 循环直到收敛: E0_far FFT(E0) # 正向传播傅里叶变换 E0_far sqrt(I_far) * exp(i * angle(E0_far)) # 焦面振幅替换为测量值 E1 IFFT(E0_far) # 逆向传播逆傅里叶变换 E1 sqrt(I_obj) * exp(i * angle(E1)) # 物面振幅替换为测量值 E0 E1这里唯一的自由度是相位两端的振幅被强行钉在实测值上。算法能做的只有不断调整相位让物面和焦面的振幅约束同时成立。当两端约束都被满足到足够好的程度时恢复出的相位就与真实波前一致。这个过程中 FFT 和 IFFT 是成对出现的实际工程里还会涉及 fftshift 和像素坐标与物理坐标的换算这些放到第 3 章的仿真平台实现里细说。2.3 收敛的边界误差曲线为什么先陡降后平台为了验证 GS 在理想条件下的行为我写了一段极简脚本先随机生成一个平滑相位正向模拟得到两组光强再跑 GS 恢复并记录焦面振幅误差。import numpy as np from scipy.ndimage import gaussian_filter N 256 phase_true gaussian_filter(np.random.randn(N, N), sigma16) def forward(field): return np.fft.fftshift(np.fft.fft2(np.fft.fftshift(field))) def gs_minimal(I_obj, I_far, phase_init, iters200): amp_obj np.sqrt(I_obj) amp_far np.sqrt(I_far) field amp_obj * np.exp(1j * phase_init) err [] for _ in range(iters): far forward(field) far amp_far * np.exp(1j * np.angle(far)) field np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(far))) field amp_obj * np.exp(1j * np.angle(field)) err.append(np.mean((np.abs(far) - amp_far) ** 2)) return np.angle(field), err obj_field np.exp(1j * phase_true) I_obj np.abs(obj_field) ** 2 far_field forward(obj_field) I_far np.abs(far_field) ** 2 rec_phase, err gs_minimal(I_obj, I_far, np.zeros((N, N)), iters300) print(err[0], err[-1])这段脚本里有几个值得注意的点。初值我用了全零在物面本身相位起伏不大的情况下GS 很容易收敛但如果真值相位 RMS 很大全零初值会很快陷入局部极小这是第 4 章要展开讨论的问题。误差记录的是焦面振幅残差平方它在一开始下降非常快几十次迭代后进入平台期这时继续加迭代次数相位残差也不会有明显改善。平台期的成因要从交替投影的几何意义解释它等价于在两个约束集合之间做最速下降而振幅匹配约束是一个非凸集合投影点不唯一迭代序列很容易被吸引到局部极小点。误差平台期是 GS 最常见的正常行为不代表算法坏了只说明该换策略了比如换初值、加支撑域或者改用带反馈的输入输出算法。3. 从光强图到波前图仿真平台的数据流与GS工程实现3.1 仿真平台的数据流设计四个模块串起来一个能反复使用、能定量评估算法改动效果的 GS 仿真平台我一般按四个模块来组织物理参数配置、正向数据生成、相位恢复算法、误差评估。各模块的输入输出关系可以用一张表说清楚模块输入输出关键依赖物理参数配置波长、像素数、像素尺寸、传播距离坐标网格、采样参数numpy正向数据生成真值相位屏、物体复振幅物面光强、焦面光强、噪声版本scipy、numpy相位恢复算法两组光强、初值、迭代参数恢复相位、迭代误差曲线numpy误差评估恢复相位、真值相位相位RMSE、相关系数、Strehl比scikit-image、numpy这样组织的平台能做的事包括批量生成随机像差验证 GS 在不同 RMS 和噪声水平下的表现对比不同初值策略对比 GS 与 HIO 的收敛曲线。我在做这类仿真时会把真实相位保留住作为真值这跟真实实验不同。真实实验里永远拿不到真值只能靠焦面光强残差间接判断是否收敛。仿真平台的价值恰恰在这里有真值才能回答“算法恢复得好不好”这个最根本的问题。3.2 生成带真值的模拟数据相位屏、衍射光强与噪声先写一个可复用的数据生成函数它负责产出一组“模拟实测”数据给定波前畸变 RMS生成平滑相位屏正向传播得到焦面光强再叠加噪声。import numpy as np from scipy.ndimage import gaussian_filter def simulate_measurement(N256, rms_wave0.3, noise_level1e-3): 生成模拟波前、物面光强、焦面光强 rms_wave : 波前畸变RMS单位是波长lambda noise_level : 焦面光强的相对高斯噪声方差 # 1. 平滑随机相位屏RMS换算到弧度 phase_true gaussian_filter(np.random.randn(N, N), sigmaN / 8.0) phase_true - phase_true.mean() phase_true * (rms_wave * 2 * np.pi) / np.sqrt(np.mean(phase_true ** 2)) # 2. 物面复振幅假设均匀振幅只看相位 field_obj np.exp(1j * phase_true) # 3. 夫琅禾费衍射到焦面 field_far np.fft.fftshift(np.fft.fft2(np.fft.fftshift(field_obj))) # 4. 光强计算与归一化 I_obj np.abs(field_obj) ** 2 I_far np.abs(field_far) ** 2 I_far I_far / I_far.max() if noise_level 0: I_far I_far noise_level * np.random.randn(N, N) I_far np.maximum(0, I_far) return phase_true, I_obj, I_far参数设置上有几个经验值。N 取 256 或 5122 的幂在 FFT 上计算效率更高。rms_wave 控制波前起伏强度实际光学加工误差通常在 1/20 到 1 个波长之间大于 1 个波长的波前畸变对 GS 来说已经属于高难度恢复。z 和像素尺寸没有直接体现在这段代码里因为夫琅禾费近似的 FFT 已经把常数相位因子略掉了如果要做菲涅尔衍射就得把传播距离 z 和像素尺寸放到传递函数里fftshift 的位置也会变化。为什么要做光强归一化GS 每次迭代都执行振幅替换如果两个平面的光强峰值差了好几个数量级振幅替换会引入数值上的不均衡。把焦面光强归一化到峰值 1再配合接下来的振幅归一化算法会稳定很多。这个细节经常被忽略但它直接影响收敛质量。3.3 GS核心迭代的工程实现终止条件与fftshift的意义把极简版扩展成带终止条件的工程函数def gs_recover(I_obj_measured, I_far_measured, phase_init, iters300, tol1e-6, verboseFalse): GS双平面相位恢复 I_obj_measured : 物面光强 I_far_measured : 焦面光强 phase_init : 物面相位初值弧度 返回: 恢复相位, 误差序列 amp_obj np.sqrt(np.maximum(I_obj_measured, 0)) amp_far np.sqrt(np.maximum(I_far_measured, 0)) # 振幅归一化让两个平面能量量级一致 amp_obj amp_obj / np.sqrt(np.mean(amp_obj ** 2)) amp_far amp_far / np.sqrt(np.mean(amp_far ** 2)) field amp_obj * np.exp(1j * phase_init) errs [] best_phase None best_err np.inf for k in range(iters): # 正向传播fftshift让零频移动到数组中心 far np.fft.fftshift(np.fft.fft2(np.fft.fftshift(field))) far amp_far * np.exp(1j * np.angle(far)) # 逆向传播 field np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(far))) field amp_obj * np.exp(1j * np.angle(field)) # 重新正变换检查焦面振幅残差 far_check np.fft.fftshift(np.fft.fft2(np.fft.fftshift(field))) err np.mean((np.abs(far_check) - amp_far) ** 2) errs.append(err) # 保存历史最优结果因为GS在接近收敛时会抖动 if err best_err: best_err err best_phase np.angle(field) if verbose and k % 50 0: print(fiter {k}, far-field amp err {err:.3e}) if err tol: break return best_phase, errs这里有两个工程点值得展开。第一个是 fftshiftFFT 得到的结果零频在四个角不 shift 的话等效于把图像平移了半个周期前后两次 shift 保证坐标和物理坐标一致。第二个是终止条件我记录每次迭代后的焦面振幅残差同时保留历史最优相位因为 GS 在接近收敛时可能在最优解附近来回抖动最后一次迭代的结果反而不一定是全过程中最好的。迭代次数一般给 300 到 500tol 给 1e-6 到 1e-8。实际跑的时候你会发现前 50 次迭代误差下降极快之后基本进入水平线。这时 tol 设得太小只会空转浪费时间更合理的做法是把 iters 设为硬上限同时用 tol 做提前退出。3.4 初值策略全零、随机还是泽尼克拟合初值选择直接决定 GS 是收敛到全局解还是局部极小。三种常见初值各有适用场景初值类型做法优点风险常数相位全零或全π收敛快低频像差效果好对RMS较大的像差容易陷入局部极小随机相位均匀分布[-π,π]有概率跳出局部极小每次结果不重复需要多次运行取最优泽尼克拟合初值先估计低频像差系数再生成初值接近真值收敛稳定需要额外的波前估计手段我一般先用全零初值跑一遍如果误差平台期过后相位残差仍然明显就换随机初值跑 5 到 10 次取焦面残差最小的一次结果。这个“多初值取最优”的做法不需要改算法只是批量跑在仿真平台里实现成本很低。平台设计时初值策略应该做成可配置参数而不是写死在代码里否则后面所有算法对比都会带上前置偏差。4. 让GS真正收敛的三个参数和两类改进算法4.1 三个必调参数迭代次数、支撑域、光强归一化迭代次数要根据误差曲线来调而不是拍脑袋给个固定值。如果跑了 500 次误差还在平台期那问题基本不在次数上而在初值或者约束上。支撑域约束是加速收敛最有效的手段当物面是已知孔径或者有限尺寸物体时可以在每次振幅替换之后把支撑域外的振幅直接置零。这会强行压缩解空间让迭代更快进入稳定区。支撑域的代码实现很简短support amp_obj 0.3 * amp_obj.max() # 或者根据真实孔径边界生成 field np.where(support, amp_obj * np.exp(1j * np.angle(field)), 0)这个操作必须放在物面振幅替换之后并且要注意支撑域给错的情况给大了等于没约束给小了会把真实分量切掉。对于平滑波前支撑域通常是整个光瞳边缘外直接置零即可。光强归一化是第三件必做事项。把物面和焦面的方根振幅都归一化到平均能量一致防止两个平面量级差异带来的数值病态。很多初次接触 GS 的人只做光强归一化而不做振幅归一化其实 GS 迭代操作的对象是振幅对振幅做归一化才更直接。4.2 迭代停滞加权GS与混合输入输出HIO当 GS 陷入停滞时一个简单的改进是不完全替换振幅。把焦面替换操作从“直接换成测量值”改成“新旧值按比例混合”def weighted_gs_recover(I_obj, I_far, phase_init, alpha0.5, iters300): amp_obj np.sqrt(np.maximum(I_obj, 0)) amp_far np.sqrt(np.maximum(I_far, 0)) amp_obj amp_obj / np.sqrt(np.mean(amp_obj ** 2)) amp_far amp_far / np.sqrt(np.mean(amp_far ** 2)) field amp_obj * np.exp(1j * phase_init) for _ in range(iters): far np.fft.fftshift(np.fft.fft2(np.fft.fftshift(field))) # 不完全替换新旧振幅按比例混合alpha越小越稳定 far_amp (1 - alpha) * np.abs(far) alpha * amp_far far far_amp * np.exp(1j * np.angle(far)) field np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(far))) field amp_obj * np.exp(1j * np.angle(field)) return np.angle(field)alpha 取值在 0.2 到 0.8 之间alpha 越小越稳定但收敛越慢。加权版本牺牲单次下降速度换取跳出局部极小的能力在噪声较强的场景下比纯 GS 健壮得多。更系统的改进是混合输入输出算法也就是工程里常说的 HIO。它的思路是在物面施加支撑域约束支撑域内接收投影结果支撑域外保留旧值并叠加一个负反馈项让错误分量被逐步推出去。def hio_recover(I_far, support, phase_init, beta0.9, iters500): 单强度相干衍射成像的混合输入输出算法 I_far : 远场衍射光强 support : 物面支撑域布尔mask beta : 反馈系数一般0.7~0.9 amp_far np.sqrt(np.maximum(I_far, 0)) amp_far amp_far / np.sqrt(np.mean(amp_far ** 2)) obj support * np.exp(1j * phase_init) for _ in range(iters): far np.fft.fftshift(np.fft.fft2(np.fft.fftshift(obj))) far amp_far * np.exp(1j * np.angle(far)) obj_proj np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(far))) # 支撑域内接受投影支撑域外保留旧值并加负反馈 obj np.where(support, obj_proj, obj - beta * obj_proj) return np.angle(obj)注意这里的 HIO 只用了一幅远场强度图物面没有振幅测量这是它与双平面 GS 的本质区别。如果物面色一个已知孔径HIO 的收敛半径通常比 GS 更大对噪声的健壮性也更好。beta 一般取 0.9反馈太强会导致振荡太弱又退化成普通投影。我习惯在仿真平台上把 GS、加权 GS 和 HIO 的误差曲线画在同一张图里对比你会看到 HIO 前期可能更慢但中期经常能走出 GS 的平台期。这不是绝对的数据集过差时 HIO 也有自己的停滞模式但作为默认对照算法它值得加进任何一套 GS 仿真平台。4.3 相位解包裹恢复结果里绕不开的2π跳变np.angle 输出的相位范围是[-π, π]。真实波前 RMS 只要超过 0.5 个波长恢复相位里就会出现大量 2π 跳变。这些跳变是反正切函数的数学特性造成的不是算法错误。直接用跳变后的相位去算 RMSE误差会大得离谱甚至得出结论“算法完全不能用”。正确的做法是先做相位解包裹把跳变后的相位还原成连续曲面from skimage.restoration import unwrap_phase phase_unwrapped unwrap_phase(recovered_phase, wrap_aroundFalse)两个使用要点。第一unwrap_phase 求解的是相对值如果真实相位里带整体倾斜解包裹结果会差一个整体倾斜项必须在评估阶段做去倾斜处理。第二wrap_around 参数用于环形数据普通波前重建不要打开。解包裹应该放在评估之前并且尽量在低光强区域生成掩膜之后再解否则噪声区域的相位会被解包裹放大成大片伪条纹这个问题在第 5 章会再展开。5. GS相位恢复的五个高频翻车场景现象、原因与解法5.1 迭代误差在几十次后进入平台期恢复的相位依然模糊现象误差曲线先陡降随后几乎不再下降恢复出来的相位和真值相比有明显差距边缘细节模糊。原因这是 GS 陷入局部极小的典型表现。振幅约束集合是非凸的交替投影无法保证全局收敛。初值质量差、波前 RMS 过大、缺少支撑域约束都会加剧这一现象。解决换随机初值跑 5 到 10 次取焦面误差最小的一次给物面加上支撑域约束或者换用加权 GS 和 HIO。在仿真平台上把关这批方案固化成自动化流程比每次手动调参要高效得多。5.2 重建波前出现棋盘格或条纹状伪影现象恢复相位里有规则的棋盘格或平行条纹看起来像高频噪声又不像真实像差。原因最常见的元凶是 fftshift 使用不一致。FFT 和 IFFT 之间如果 shift 的次数不匹配等效于图像做了亚像素平移高频分量会产生周期性伪影。另一个来源是图像边界不连续。FFT 默认把图像左右、上下边界当作周期性边界处理如果光强图边界有陡峭跳变就会出现网格状振铃。解决正向传播和逆向传播保持成对的 fftshift 调用输入光强先做零填充或者加窗函数再接 FFT。零填充直接在数组外围扩展一圈零值可以显著缓解边界效应。5.3 两个平面的光强中心对不齐重建波前整体漂移现象恢复相位和真值相比有整体平移同一块像差出现在不同位置。原因物面和焦面坐标系基准不一致。比如物面用了 fftshift 而焦面没用的或者实验采集时两个探测器像素起点不同。仿真平台里这个错位不容易出现但换真实数据时几乎必踩。解决跑 GS 之前先做质心对齐把两幅光强图的质心移动到数组中心。对仿真平台从物理参数生成坐标时就把原点统一到中心对真实数据这个对齐步骤不能省直接决定重建波前是否与实测坐标系吻合。5.4 光强带噪声时算法发散或误差在高位振荡现象噪声稍微大一点GS 误差曲线不再下降而是在某个数值附近持续跳动恢复结果里有大量椒盐状噪声。原因振幅替换会把光强噪声直接转成振幅扰动整个迭代没有任何正则化机制去抑制它。泊松噪声、高斯噪声都会被逐次放大尤其是在焦面光强接近零的区域噪声相对幅度会非常高。解决光强数据先滤波中值滤波或高斯滤波都可以改用加权 GS把替换强度从 1 降到 0.5 左右或者在物面加 TV 正则化约束。仿真时如果模拟光子散粒噪声按泊松分布建模比高斯噪声更贴近实际探测器行为。5.5 相位解包裹后出现大片错误条纹现象unwrap_phase 之后原本只是局部噪声的区域被放大成连片条纹整个波前图看起来像被污染过。原因解包裹算法依赖相邻像素的相位差低光强区域的相位是噪声主导那里的相位梯度会被误判为真实跳变解包裹结果被这些错误跳变污染并通过路径传播扩散开。解决先根据光强分布生成一个有效区域 mask只对 mask 内像素做解包裹。mask 之外不参与相位评估计算 RMSE 也放在 mask 内进行。这相当于把低信噪比区域的干扰彻底隔离掉是最稳妥的评估方式。6. 用相位RMSE和Strehl比给GS重建结果打分验证一套 GS 算法实现是否真的可用我一般同时看三个指标相位 RMSE、复振幅相关系数和 Strehl 比。相位 RMSE 要在解包裹和去倾斜之后计算否则数值没有意义。相关系数衡量整体分布相似度对局部细节误差不敏感。Strehl 比是光学系统点扩散函数峰值强度的归一化指标对残余像差非常敏感接近 1 表示波前质量接近理想衍射极限。指标计算方式适用场景相位RMSE解包裹后逐像素求差去均值后算均方根总体误差量化相关系数恢复相位与真值相位归一化互相关形状相似度判断Strehl比残余相位在瞳孔上的衍射峰值与理想值之比光学质量评估仿真平台上计算 Strehl 比可以借助一版很短的代码。核心是用残余相位在瞳孔区域上做积分再归一化到理论最大值def strehl_ratio(phase_residual, aperture): phase_residual : 残余相位恢复相位-真值相位解包裹后 aperture : 瞳孔mask1为通光区域 if aperture is None: aperture np.ones_like(phase_residual) n_pixels np.sum(aperture) E np.exp(1j * phase_residual) * aperture return np.abs(np.sum(E)) ** 2 / n_pixels ** 2需要提醒的是Strehl 比只有在残余相位 RMS 小于约 1 弧度时有明确物理意义残余像差太大时它迅速逼近 0区分度反而下降。这时候要同时看衍射光强分布图而不是只盯一个数字。固定评估流程是我吃过亏之后养成的习惯解包裹去倾斜生成有效区域 mask最后三个指标一起算。我最早就因为没解包裹、没去倾斜就直接拿恢复相位和真值算 RMSE得到一份很难看的数字一度差点否定了整条算法路线。后来把评估流程固定下来才真正敢用 GS 的结果做工程判断。希望你在一开始就把评估流程搭好省掉这一步的弯路。希望帮到你。本文还有配套的精品资源点击获取
