GS 算法详解：从相位恢复到纯相位全息图

如果相机只能拍到光强，怎样知道光波的相位？反过来，如果一块空间光调制器只能改变相位，怎样让一束高斯光在远处形成指定的光强图案？

这两个问题分别属于相位恢复与计算全息设计。它们都可以从同一个简单的迭代思想出发：在两个由传播关系连接的平面之间往返，每到一个平面，就把振幅改成已知值，同时保留当前估计的相位。

这就是本文讨论的 Gerchberg–Saxton 算法，简称 GS 算法。这里的 GS 不是数值线性代数中的 Gauss–Seidel 迭代。Gerchberg 与 Saxton 在 1972 年提出的方法针对成像面和衍射面强度已知时的相位问题；此后，这种迭代结构也成为计算全息的重要起点。早期工作的背景可参见 Gerchberg 的同期实验论文。

GS 的两个已知振幅平面，以及需要设计的输入相位

图 1｜已知两端的光强，中间缺少的是相位。中图是本文九光斑算例计算出的输入相位，左右分别为入射光强和目标光强。两张光强图分别按自身范围显示，不能凭颜色比较两平面的峰值；相位色标单位为弧度，孔径外振幅为零、相位无物理意义，统一用浅灰色遮罩。本文六组图均由明确的公式、几何或数值计算生成，源码见文末。

本文先建立数学模型，再推导每一步为什么这样做，最后用可复现的九光斑算例解释算法的能力和局限。阅读时最需要分清三件事：振幅不等于强度、强制约束后的场不等于实际传播的场、误差下降不等于相位唯一或全局最优。

1. 相机究竟丢掉了什么？

1.1 复光场、振幅、相位与强度

在单色、标量模型下，一个横截面上的复光场可以写成

u(x,y)=a(x,y)e^{i\phi(x,y)},\qquad a(x,y)\geq0.

其中 $a$ 是振幅， $\phi$ 是相位。省略介质阻抗等固定系数后，相机测到的强度是

$$ I(x,y)=|u(x,y)|^2=a^2(x,y). $$

因此，从强度得到振幅要取平方根： $a=\sqrt I$ 。如果把灰度图的像素值定义为目标强度，送入 GS 振幅约束的应当是灰度值的平方根。直接使用灰度值作为振幅，会把期望的强度轮廓再次平方。

相位描述不同位置的波振动到哪一步。相机不能在一个像素的强度读数中直接给出这个信息，但相位会影响不同位置的光传播后如何相干叠加，所以会改变另一个平面的强度。GS 就是利用传播把相位信息间接编码到第二张强度图中的事实。

1.2 为什么不是一次逆傅里叶变换就够了？

设目标面复场为 $v=b e^{i\psi}$ 。如果振幅 $b$ 和相位 $\psi$ 都已知，那么 $u=\mathcal F^{-1}v$ 就给出输入复场。但实验或设计任务通常只给出 $b$ 。

直接计算 $\mathcal F^{-1}b$ ，实际上偷偷假设了目标面相位处处为零。得到的输入振幅通常既不等于测得的物面振幅，也不等于照射在 SLM 上的固定光束振幅。GS 的迭代，就是逐步调整这个未知目标相位，使反传结果更符合输入面的要求。

2. 一个数学框架，两类不同任务

本文将离散输入场记为 $u$ ，输出场记为 $v$ ，并采用单位化二维离散傅里叶变换：

v=\mathcal Fu,\qquad u=\mathcal F^{-1}v, \qquad \|\mathcal Fu\|_2=\|u\|_2.

最后一个等式是 Parseval 能量关系。这里的范数包含数组中全部像素。

问题	输入面已知什么	输出面已知什么	最终想得到什么
双平面相位恢复	测得的振幅 $a=\sqrt{I_{\mathrm{in}}}$	测得的振幅 $b=\sqrt{I_{\mathrm{out}}}$	与测量及先验相容的未知相位
纯相位全息设计	光束和孔径决定的振幅 $a$	人为设定的目标振幅 $b$	能产生接近目标强度的输入相位

相位恢复中，理想无噪声数据应来自同一个真实场；全息设计中，两个振幅往往是分别指定的，未必存在任何相位能同时精确满足它们。后者即使实现完全正确，也可能只能得到近似解。

2.1 先满足必要的功率条件

在本文单位化、完整采样且无损的离散模型中，精确解必须满足

\sum_{m,n}a_{mn}^2=\sum_{p,q}b_{pq}^2.

设计目标只指定形状时，可以先令 $b_0=\sqrt{I_{\mathrm{target}}}$ ，再统一缩放：

b=b_0\frac{\|a\|_2}{\|b_0\|_2}.

这是功率匹配，不代表两个约束已经兼容。若功率没有匹配，反三角不等式直接给出

\big\||\mathcal Fu|-b\big\|_2 \geq \left|\|a\|_2-\|b\|_2\right|.

也就是说，单是总功率不一致，就已经规定了一个无法靠迭代消除的误差下限。

真实测量要考虑像素面积、曝光时间、相机增益、光路损耗以及有限探测范围。不能将两幅未经标定、视场不同的相机照片简单归一化，就宣称它们构成一致的完整双平面数据。

3. 为什么“保留相位、替换振幅”有道理？

暂时只看一个像素。当前复数为 $z=r e^{i\theta}$ ，已知这个像素应有振幅 $b$ ，所以允许的复数必须落在半径为 $b$ 的圆周上。

在这些允许点中，寻找离当前值最近的点：

\min_{\alpha}\left|r e^{i\theta}-b e^{i\alpha}\right|^2.

将模平方展开：

\left|r e^{i\theta}-b e^{i\alpha}\right|^2 =r^2+b^2-2rb\cos(\theta-\alpha).

当 $r,b>0$ 时，余弦取最大值，距离就最小，因此 $\alpha=\theta\pmod{2\pi}$ 。最小距离是 $|r-b|$ 。于是逐像素投影可以写成

P_b(z)=b e^{i\arg z}.

复平面上的振幅投影：保持相角，将半径改为目标振幅

图 2｜振幅替换是欧氏距离意义下的最近点投影。圆周上的每一点都满足振幅要求，但集合不是凸集：两个允许点的平均值一般不在圆周上。图中半径归一为 1。

两个边界情况需要说清楚：当 $b=0$ 时，投影就是零；当 $z=0$ 而 $b>0$ 时，相位没有定义，圆周上所有点都同样近。代码中使用 np.angle(0)=0 是一种确定性的选择，也可以采用已有相位等约定。

不要不加区分地写成 b * z / (abs(z) + eps)：它避免了除零，却让非零点的输出振幅也偏离 $b$ ，在严格零点更无法产生目标非零振幅。若要实现数学上的精确振幅投影，显式处理相位或零点更清楚。

4. 从初值到完整的一次迭代

4.1 初始化

给输入面一个相位初值 $\phi_0$ ，构造

u_0=a e^{i\phi_0}.

常见选择是独立随机相位，也可以使用平滑相位、透镜相位或上一次优化结果。初值会影响迭代路径。零相位并非永远错误，但在对称目标上可能让迭代长时间停留在特殊的对称结构中。

4.2 正向传播，并观察当前误差

v_k=\mathcal Fu_k=|v_k|e^{i\psi_k}.

此时 $|v_k|^2$ 才是第 $k$ 个输入相位实际产生的模拟强度。它通常不等于目标 $b^2$ 。

4.3 在目标面施加振幅约束

\hat v_k=P_b(v_k)=b e^{i\psi_k}.

帽子表示“施加约束后的场”。这一步人为使目标振幅正确，但 $\hat v_k$ 不一定能由当前输入振幅和某个可加载相位实现。

4.4 反向传播，再施加输入振幅约束

\tilde u_k=\mathcal F^{-1}\hat v_k, \qquad u_{k+1}=P_a(\tilde u_k) =a e^{i\arg\tilde u_k}.

反传得到的振幅通常偏离 $a$ 。将它恢复为 $a$ 后，保留下来的新相位就是下一轮使用的相位。两端约束互相影响，必须反复协调。

GS 正向变换、目标振幅投影、逆向变换和输入振幅投影的完整流程

图 3｜每一轮由两次变换和两次振幅投影构成。每次计算质量指标，都要使用目标面投影前的传播场。

将上述步骤合成一个式子：

\boxed{u_{k+1}=P_a\mathcal F^{-1}P_b\mathcal Fu_k.}

结束迭代后，应重新从最后的 $u_K$ 正向传播一次，计算 $v_K=\mathcal Fu_K$ 。如果输出了更新后的相位，却报告更新前的强度或误差，就会出现相位图与验证结果相差一轮的问题。

5. GS 的误差为什么会下降，又为什么仍会失败？

5.1 将两个约束放进同一个空间

定义输入振幅约束集合和傅里叶振幅约束集合：

\mathcal A=\{u:|u|=a\},\qquad \mathcal B=\{u:|\mathcal Fu|=b\}.

因为 $\mathcal F$ 是酉变换，所以 $P_{\mathcal B}=\mathcal F^{-1}P_b\mathcal F$ 也是对应欧氏距离的最近点投影。GS 就是在 $\mathcal A$ 和 $\mathcal B$ 之间交替投影。

5.2 误差不增的简洁推导

令 $w_k=P_{\mathcal B}u_k$ ，定义未归一化振幅误差

d_k=\big\||\mathcal Fu_k|-b\big\|_2 =\|u_k-w_k\|_2.

由于 $u_{k+1}$ 是 $w_k$ 在 $\mathcal A$ 中的最近点，而旧点 $u_k$ 本来就属于 $\mathcal A$ ，所以

\|u_{k+1}-w_k\|_2\leq\|u_k-w_k\|_2.

下一轮中， $w_{k+1}$ 是 $u_{k+1}$ 在 $\mathcal B$ 中的最近点，而 $w_k$ 也属于 $\mathcal B$ ，所以

\begin{aligned} d_{k+1} &=\|u_{k+1}-w_{k+1}\|_2\\ &\leq\|u_{k+1}-w_k\|_2\\ &\leq d_k. \end{aligned}

因此，在固定振幅约束、精确投影、酉传播算子这些条件下，目标面振幅误差不增加。这里将经典误差下降论证用统一范数和投影记号展开；相关算法比较与原始讨论见 Fienup，1982。

5.3 这个证明没有保证什么？

它保证非负标量序列 $d_k$ 有极限，却没有证明这个极限为零，也没有证明复场序列一定收敛到唯一解。

原因在于：两个约束是非凸的，可能不相交；即使存在交点，局部迭代也可能停在误差非零的停滞状态，或者很长时间只发生极小改进。于是“迭代 300 次后曲线平了”不能直接翻译为“求出了正确相位”。

更不能把上述结论直接套到其他指标。振幅误差不增，并不保证强度误差、每个光斑的均匀性或某个区域内的效率逐轮改善。 对加权目标、量化相位、裁剪传播、多平面组合等修改，也必须重新检查相应的算子和目标函数。

5.4 相位本身可能不唯一

如果 $u$ 是一个解，那么 $u e^{i\theta_0}$ 具有相同的输入、输出强度。因此，全局相位通常不可观测。只有单幅傅里叶强度时，还常有平移、共轭反演等歧义；增加已知输入振幅或其他先验后，这些变换未必仍然满足全部约束，不能一概保留或一概排除。

比较恢复相位和真值时，需要先消除可接受的全局相位差，并用环绕相位差 $\arg(e^{i(\phi_{\mathrm{rec}}-\phi_{\mathrm{true}})})$ ，而不是直接相减。振幅为零的位置没有可比较的物理相位，应从相位误差统计中排除。

6. 怎样评价一个 GS 结果？

最基本的指标是实际传播场的相对振幅误差：

E_A(k)=\frac{\big\||v_k|-b\big\|_2}{\|b\|_2}.

全息图的视觉效果直接对应强度，因此还应报告

E_I(k)=\frac{\big\||v_k|^2-b^2\big\|_2}{\|b^2\|_2}.

这两个量都无量纲，但含义不同。由于 $I=a^2$ ，在较亮区域的小振幅偏差会变成更显著的强度偏差；不能将 $E_A=0.1$ 直接称为“强度误差 10%”。

如果关心目标区域 $\Omega$ 内收集到多少光，可以定义

\eta_\Omega=\frac{\sum_{(p,q)\in\Omega}|v_{pq}|^2} {\sum_{p,q}|v_{pq}|^2}.

对光斑阵列，还可先在每个光斑周围同样大小的窗口中积分，得到功率 $P_j$ ，再用 $\mathrm{CV}=\mathrm{std}(P_j)/\mathrm{mean}(P_j)$ 衡量均匀性。窗口尺寸必须一起报告，否则不同算例的效率或均匀性无法直接比较。

最容易制造“完美结果”的错误，是拿 $\hat v_k$ 算误差。因为它的振幅被强制设为 $b$ ，此时误差几乎必然是零。正确做法是用可实现输入场 $a e^{i\phi_k}$ 再传播，评估得到的 $v_k$ 。

停止条件可以结合最大迭代次数、目标误差阈值，以及若干轮内的相对改进量。只有“改进很小”而没有达到任务精度要求时，应称为停滞，不能称为任务已经达标。

7. 可复现实验：让高斯光变成九个光斑

7.1 数值条件

这里设计的是纯相位全息图，不假设存在一个等待恢复的真实相位。模型采用无损、单色、标量、理想傅里叶传播；以下像素坐标都是模拟网格坐标。

参数	取值
网格	$256\times256$
输入振幅	$a(x,y)\propto\exp[-(x^2+y^2)/65^2]$
输入圆孔径	半径 100 像素，外部振幅为零
目标	$3\times3$ 个等强二维高斯光斑
光斑中心	横、纵坐标各取 $-28,0,28$ 像素
单光斑强度标准差	$\sigma=3.5$ 像素
总功率	输入和目标均归一为 1
初相位	NumPy `default_rng(7)`，在 $[-\pi,\pi)$ 均匀分布
GS 更新次数	300

输入高斯表达式中的 65 是振幅降到 $1/e$ 的半径，也是强度降到 $1/e^2$ 的半径。目标中的 3.5 则是强度高斯的标准差，两种参数定义不能混淆。

九光斑目标以及第 0、1、10、50、300 次更新后的真实模拟强度

图 4｜目标及迭代结果的中心区域，所有子图共享同一个线性色标；超过目标峰值的部分按色标上限饱和显示，峰值差异另见图 5 的定量截面。第 0 次是随机相位初始结果；后续每幅图都从满足输入振幅约束的相位重新传播得到。显示裁剪到中心区域，但迭代及全局误差使用完整的 $256\times256$ 网格；初始随机场的大量功率分散在图框之外。

第一轮更新就把能量明显集中到目标区域，但继续迭代后仍可看到背景、局部起伏和光斑形状误差。GS 很快产生“像目标的图”，精确逼近则困难得多。

7.2 实际计算结果

完成的更新次数	相对振幅误差 $E_A$
0	1.277382
1	0.478160
10	0.353803
50	0.308642
300	0.288974

300 次更新后的相对强度误差为 0.209418。以九个光斑中心为圆心、半径 8 像素的九个不重叠圆盘作为信号区域，重建光的区域功率占比为 **84.02%**，九个圆盘积分功率的变异系数为 **4.21%**。

同样的圆盘只包含目标高斯总功率的 **92.24%**，因为高斯本身有尾部。因此，84.02% 是针对该区域定义的模拟功率占比，不能直接等同于 SLM 的实验衍射效率，也不能把区域外所有光都算成算法产生的杂散光。

图 5｜左：随机种子 1、7、11、23 的振幅误差都不增，但到达不同的最终值。右：种子 7 的中心横截面，直接比较实际强度与目标强度，未对两条曲线分别按峰值归一化。

四个种子的最终振幅误差约为 0.295570、0.288974、0.291237、0.296787。这说明多初值比较有实际价值，但这四次计算既不是全局最优性证明，也不能证明目标绝对不可实现。

7.3 数值实现做了哪些检查？

随文程序检查了奇数和偶数网格上的 FFT 往返、Parseval 关系、与显式 DFT 矩阵的独立比较、输入振幅约束、全局相位不改变强度，以及全部四条误差曲线的不增性质。

本次计算中，FFT 往返相对误差小于 $4.4\times10^{-16}$ ，显式 DFT 与 FFT 的相对差为 $2.03\times10^{-15}$ 以内，最终功率相对偏差小于 $2.3\times10^{-16}$ 。这些结果验证了离散实现内部的一致性，并不替代实验光路、像素标定或真实器件响应的验证。

8. Python 实现：核心并不长

下面代码给出与本文算例相同的核心算法。数组必须同形状、有限、非负且具有非零总功率；示例在进入迭代前完成统一归一化。

import numpy as np

def ft(u):
    return np.fft.fftshift(
        np.fft.fft2(np.fft.ifftshift(u), norm="ortho")
    )

def ift(v):
    return np.fft.fftshift(
        np.fft.ifft2(np.fft.ifftshift(v), norm="ortho")
    )

def gs(a, b, steps=300, seed=7):
    rng = np.random.default_rng(seed)
    phi = rng.uniform(-np.pi, np.pi, a.shape)
    u = a * np.exp(1j * phi)
    errors = []

    for k in range(steps + 1):
        v = ft(u)
        # Measure BEFORE imposing the target amplitude.
        errors.append(np.linalg.norm(np.abs(v) - b)
                      / np.linalg.norm(b))
        if k == steps:
            break
        v_projected = b * np.exp(1j * np.angle(v))
        u_back = ift(v_projected)
        u = a * np.exp(1j * np.angle(u_back))

    # v corresponds to exactly the returned phase.
    return np.angle(u), np.abs(v)**2, np.asarray(errors)

n = 256
q = np.arange(n) - n // 2
x, y = np.meshgrid(q, q)
a = np.exp(-(x*x + y*y) / 65**2)
a *= (x*x + y*y <= 100**2)
a /= np.linalg.norm(a)

target = sum(
    np.exp(-((x-cx)**2 + (y-cy)**2) / (2*3.5**2))
    for cy in (-28, 0, 28)
    for cx in (-28, 0, 28)
)
target /= target.sum()
b = np.sqrt(target)
phase, reconstructed, error = gs(a, b)
print(error[0], error[-1])

norm="ortho" 使变换酉归一化，输入、输出数组的模平方和可以直接比较。常规 FFT 默认归一化也能实现 GS，但必须一致处理功率尺度及误差定义。

ifftshift → fft2/ifft2 → fftshift 则保证以数组中央作为坐标原点。不要因为偶数网格下 fftshift 和 ifftshift 看起来可以互换，就在所有位置使用同一个移位操作；奇数网格会暴露这种错误。

如果总像素数为 $M=N_xN_y$ ，一轮需要两次二维 FFT，计算复杂度约为 $O(M\log M)$ ，主工作数组的内存复杂度为 $O(M)$ 。目标图、相位图、多个中间复场和保存的历史都会增加实际内存用量。

完整程序还会绘制本文六组图，比较四个初值，并导出误差指标与相位数组。它是可下载后直接运行的单文件程序，不依赖本博客代码。

9. FFT 数组怎样对应真实光学坐标？

9.1 傅里叶变换成立的光学条件

本文的传播关系适用于理想远场或适当的透镜傅里叶面模型。下面画出输入平面位于薄透镜前焦面、观察平面位于后焦面的标准结构，假定单色、近轴标量传播以及足够大的透镜孔径。

SLM 输入面、透镜与后焦面相机的傅里叶光学示意图

图 6｜理想前焦面—后焦面系统，图中距离均为 $f$ 。反射式 SLM 在实际光路中还需要折转、偏振设置和可能的中继系统；图示只用于说明传播模型。

物理后焦面场与输入面连续傅里叶变换对应，其频率坐标满足

\nu_x=\frac{X}{\lambda f},\qquad \nu_y=\frac{Y}{\lambda f}.

整体常数、坐标尺度和由具体参考平面产生的已知相位因子，要与光学布局一致。若 SLM 实际位于透镜面，而不是前焦面，后焦面表达式会包含相应的已知相位因子；即使它不改变单个观察面的强度，也应在反向传播或扩展到其他平面时正确处理。

9.2 像素间距决定输出采样

输入面计算像素间距为 $\Delta x$ 、横向点数为 $N_x$ 时，DFT 的空间频率间距为 $1/(N_x\Delta x)$ 。于是后焦面像素间距为

\Delta X=\frac{\lambda f}{N_x\Delta x},\qquad L_X=N_x\Delta X=\frac{\lambda f}{\Delta x}.

例如使用 $\lambda=532\thinspace\mathrm{nm}$ 、 $f=200\thinspace\mathrm{mm}$ 、 $\Delta x=8\thinspace\mu\mathrm m$ 和 $N_x=256$ ，则

\Delta X\approx51.95\,\mu\mathrm m, \qquad L_X\approx13.30\,\mathrm{mm}.

这是把当前 256 点网格解释为 8 微米采样间距时的示例，不是某款 SLM 的规格。若有中继放大率，要使用折算到所选输入平面的有效采样间距；相机坐标也要按成像倍率配准。

9.3 补零、分辨率与视场不要混为一谈

在固定输入采样间距和实际孔径下增加外围零填充，会减小输出采样间距，帮助更细地采样衍射图样，但不会增加真实孔径，也不会自动提升物理分辨率。物理细节尺度还受有效口径、波长和数值孔径约束。

对于全息设计，零填充区域应始终属于输入零振幅区域，不能把新增像素当成可控制的 SLM 像素。DFT 还隐含周期延拓；图样靠近计算边界、传播区域过小或裁剪过强时，需要检查重复像和混叠影响。

10. 从理想相位图到 SLM，还差哪些步骤？

算法输出的是弧度相位，SLM 接收的通常是数字灰度。最简单的线性映射只是理想化示例，真实液晶器件应使用工作波长、偏振和温度条件下标定的灰度—相位响应。

首先，用实际入射光束的振幅和有效孔径替换理想高斯模型。其次，确认器件相位范围、位深、像素填充因子和坏点，并把这些约束反映到前向计算中。若只在最后将连续相位量化，必须用量化后的相位重新计算结果。

常见加载相位还包含用于选级的载波和像差补偿：

\phi_{\mathrm{load}}=\operatorname{mod} \left(\phi_{\mathrm{GS}}+\phi_{\mathrm{carrier}} -\phi_{\mathrm{aberr}},2\pi\right).

这里 $\phi_{\mathrm{aberr}}$ 指与该输入平面等效的已测相位像差，符号约定应与标定一致。不同位置的像差不一定都能直接压缩成这张单平面相位图。线性载波会平移衍射图样，有助于通过空间滤波区分目标衍射级和未调制零级光，但要同时检查可用带宽、目标偏移和系统口径。

实验检验的顺序应当是：加载实际可实现相位，拍摄目标面，做暗场及响应校正、坐标配准，再比较强度和积分功率。若模型与实测存在系统偏差，可以引入相机反馈来更新目标权重或补偿像差。这样形成的闭环包含真实测量误差，不能用理想仿真结果代替。

11. 什么时候需要 GS 的改进版本？

11.1 多初值：先用最直接的办法诊断停滞

使用多个随机种子，在相同计算预算下比较最终误差、效率与均匀性，往往比只给单个初值无限增加迭代次数更有信息量。不要只挑最好的一张图；至少保留初值、迭代次数和评价指标，才能复现。

11.2 加权 GS：针对光斑功率不均

对离散光斑阵列，可以根据当前各光斑的功率偏差调整下一轮目标权重，让偏弱的光斑得到更大的振幅要求。一个示意性的阻尼反馈形式是

w_j^{(k+1)}=w_j^{(k)} \left(\frac{P_j^\star}{P_j^{(k)}+\epsilon}\right)^{\gamma/2}, \qquad 0<\gamma\leq1.

这里 $w_j$ 是目标光斑振幅权重， $P_j^\star$ 和 $P_j^{(k)}$ 分别是期望和当前积分功率。指数中的 $1/2$ 来自振幅与功率的平方关系； $\epsilon$ 与功率同单位，用于限制暗光斑引起的数值发散。

这只是便于理解的反馈构造，不是所有“加权 GS”实现的统一公式。实际还需定义测量窗口、限制权重、统一功率尺度，并重新验证稳定性。因为目标权重随迭代改变，前面对固定目标的单调性证明不能直接沿用。

11.3 信号区与噪声区：把自由度用在重要位置

若实验只关心输出面的有限区域，就没必要要求整个平面都精确匹配一幅图。混合区域振幅自由度方法（MRAF）在信号区控制目标振幅，同时允许其他区域有更大振幅自由度，以换取关注区域的质量；代价可能是更多功率进入区域外。具体更新及性能应以 Pasienski 与 DeMarco 的原论文为准。

这也提示一个设计原则：背景是否必须严格为零，是物理任务要求，不应只是“目标图片的黑色像素”无意中施加的苛刻约束。

11.4 更换约束或传播模型

只有傅里叶强度、物面只有支持域或非负性先验时，问题与双振幅约束 GS 已不同，可以考虑误差下降法的其他形式或 HIO 等策略。HIO 包含反馈更新，不能简单理解为“把 GS 多跑几次”。

对于近场、多个离焦平面或多层目标，需要将 $\mathcal F$ 换成正确的菲涅耳或角谱传播算子，并为各平面定义约束与损失。在多个目标平面上达到一个约束，可能破坏另一个约束，更新顺序及权重也需要设计。

尤其要区分逆算子和伴随算子：无损且离散为酉算子的模型中二者相同；带孔径裁剪、探测丢失、吸收或非酉重采样时通常不相同。不能在有信息丢失的模型中直接“反向传播一下”，就假设上述误差下降证明仍然成立。此时常需要显式优化目标及匹配的梯度、正则化或统计噪声模型。

12. 常见问题与排查顺序

现象	优先检查
重建亮暗层次不对	是否把目标强度误当成振幅，漏了平方根
图案跑到角落、镜像或偏一像素	移位约定、变换符号、坐标方向和奇偶网格
第一轮之后报告误差为零	是否评估了人工替换振幅后的场
模拟误差长期平台	多初值对比、功率匹配、约束兼容性、目标是否过苛刻
固定目标下振幅误差明显上升	正反变换是否配对，投影是否精确，有无额外量化或裁剪
仿真好，实验差	入射光束、像差、偏振、SLM 相位标定、零级光及相机配准
光斑看起来一样亮，但捕获效果不同	比较同一窗口的积分功率，而不仅是单像素峰值
换了网格尺寸，物理光斑位置变了	是否同步更新了傅里叶面采样和物理坐标

理解 GS 最关键的一点，是把它看成一个由传播模型与约束集合共同定义的求解过程。相位不是凭空猜出来的，而是在两端振幅约束的反复协调中更新；结果是否可信，则要由实际可实现的输入场重新传播，并结合误差、能量、采样和实验标定来判断。

参考文献与复现资料

R. W. Gerchberg and W. O. Saxton, A practical algorithm for the determination of phase from image and diffraction plane pictures, Optik 35, 237–246 (1972)。GS 原始论文。同期应用与该文引用见 Gerchberg, Nature 240, 404 (1972)。
J. R. Fienup, Phase retrieval algorithms: a comparison, Applied Optics 21, 2758–2769 (1982)。作者机构提供的全文，DOI。
M. Pasienski and B. DeMarco, A high-accuracy algorithm for designing arbitrary holographic atom traps, Optics Express 16, 2176–2190 (2008)。作者预印本，DOI。

本文的单像素投影展开、数值例子、作图和实现检查均在文中或源码中给出，图 1–6 为本文原创示意与模拟，不是论文实验图片。

安装 NumPy 和 Matplotlib 后，运行 python gs_demo.py 即可在当前目录的 gs-output/ 下生成全部配图、结果数组和验证报告。中文图注需要系统中有微软雅黑或黑体等字体；缺少中文字体只影响标注显示，不影响计算。本文数值结果使用 NumPy 1.26.4、Matplotlib 3.10.0 生成。