企业官网建设流程全解析

简介：基于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), sigma=16) def forward(field): return np.fft.fftshift(np.fft.fft2(np.fft.fftshift(field))) def gs_minimal(I_obj, I_far, phase_init, iters=200): 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)), iters=300) print(err[0], err[-1])

这段脚本里有几个值得注意的点。初值我用了全零，在物面本身相位起伏不大的情况下，GS 很容易收敛；但如果真值相位 RMS 很大，全零初值会很快陷入局部极小，这是第 4 章要展开讨论的问题。误差记录的是焦面振幅残差平方，它在一开始下降非常快，几十次迭代后进入平台期，这时继续加迭代次数，相位残差也不会有明显改善。

平台期的成因要从交替投影的几何意义解释：它等价于在两个约束集合之间做最速下降，而振幅匹配约束是一个非凸集合，投影点不唯一，迭代序列很容易被吸引到局部极小点。误差平台期是 GS 最常见的正常行为，不代表算法坏了，只说明该换策略了，比如换初值、加支撑域，或者改用带反馈的输入输出算法。

3. 从光强图到波前图：仿真平台的数据流与GS工程实现

3.1 仿真平台的数据流设计：四个模块串起来

一个能反复使用、能定量评估算法改动效果的 GS 仿真平台，我一般按四个模块来组织：物理参数配置、正向数据生成、相位恢复算法、误差评估。各模块的输入输出关系可以用一张表说清楚：

这样组织的平台能做的事包括：批量生成随机像差，验证 GS 在不同 RMS 和噪声水平下的表现；对比不同初值策略；对比 GS 与 HIO 的收敛曲线。我在做这类仿真时会把真实相位保留住作为真值，这跟真实实验不同。真实实验里永远拿不到真值，只能靠焦面光强残差间接判断是否收敛。仿真平台的价值恰恰在这里：有真值，才能回答“算法恢复得好不好”这个最根本的问题。

3.2 生成带真值的模拟数据：相位屏、衍射光强与噪声

先写一个可复用的数据生成函数，它负责产出一组“模拟实测”数据：给定波前畸变 RMS，生成平滑相位屏，正向传播得到焦面光强，再叠加噪声。

import numpy as np from scipy.ndimage import gaussian_filter def simulate_measurement(N=256, rms_wave=0.3, noise_level=1e-3): """ 生成模拟波前、物面光强、焦面光强 rms_wave : 波前畸变RMS，单位是波长（lambda） noise_level : 焦面光强的相对高斯噪声方差 """ # 1. 平滑随机相位屏，RMS换算到弧度 phase_true = gaussian_filter(np.random.randn(N, N), sigma=N / 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 或 512，2 的幂在 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, iters=300, tol=1e-6, verbose=False): """ 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(f"iter {k}, far-field amp err {err:.3e}") if err < tol: break return best_phase, errs

这里有两个工程点值得展开。第一个是 fftshift：FFT 得到的结果零频在四个角，不 shift 的话等效于把图像平移了半个周期，前后两次 shift 保证坐标和物理坐标一致。第二个是终止条件：我记录每次迭代后的焦面振幅残差，同时保留历史最优相位，因为 GS 在接近收敛时可能在最优解附近来回抖动，最后一次迭代的结果反而不一定是全过程中最好的。

迭代次数一般给 300 到 500，tol 给 1e-6 到 1e-8。实际跑的时候你会发现前 50 次迭代误差下降极快，之后基本进入水平线。这时 tol 设得太小只会空转浪费时间，更合理的做法是把 iters 设为硬上限，同时用 tol 做提前退出。

3.4 初值策略：全零、随机还是泽尼克拟合

初值选择直接决定 GS 是收敛到全局解还是局部极小。三种常见初值各有适用场景：

我一般先用全零初值跑一遍，如果误差平台期过后相位残差仍然明显，就换随机初值跑 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, alpha=0.5, iters=300): 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, beta=0.9, iters=500): """ 单强度相干衍射成像的混合输入输出算法 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_around=False)

两个使用要点。第一，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 表示波前质量接近理想衍射极限。

仿真平台上计算 Strehl 比可以借助一版很短的代码。核心是用残余相位在瞳孔区域上做积分，再归一化到理论最大值：

def strehl_ratio(phase_residual, aperture): """ phase_residual : 残余相位，恢复相位-真值相位，解包裹后 aperture : 瞳孔mask，1为通光区域 """ 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 的结果做工程判断。希望你在一开始就把评估流程搭好，省掉这一步的弯路。希望帮到你。

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