数值光学仿真前，为什么需要先做 Sanity Check

可运行代码

完整项目位于 GitHub 仓库，代码按职责组织如下：

一次性 sanity check 报告：propagation_sanity_check.py

纯辅助函数（不打印输出）：grid_sanity_check.py、field_support_info.py、output_grid_info.py 和 propagation_regime.py

完整可运行实验与绘图流程：three-propagation-models-clean-version.ipynb

即使物理方程和代码都没有错，数值波传播依然可能失败。

在比较 Fresnel、Fraunhofer 和角谱传播时，我遇到过这样一个结果：乍看合理，但一做定量比较，行为就不对了。问题并不出在衍射理论本身，而是出在用于表示光场的数值网格上。

这次实验带来了一个很简单的结论：

$$ \text{正确的传播方程} \neq \text{可靠的数值仿真} $$

在解释一个模拟光场之前，应该先检查它的数值采样条件。

1. 最初的实验

实验使用一个方形孔径，并分别通过三种方法传播：

角谱法（Angular Spectrum Method，ASM）
Fresnel 衍射
Fraunhofer 衍射

初始参数约为

$$ N=512, \qquad \Delta x=2\,\mu m, \qquad \lambda=532\,nm, $$

方形孔径宽度为

$$ a=100\,\mu m. $$

我将同一个光场传播到几个不同距离，并比较归一化后的中心强度剖面，如图 1 所示。

N=512、z=1 mm 时 ASM、Fresnel 与 Fraunhofer 的传播强度

N=512、z=10 mm 时 ASM、Fresnel 与 Fraunhofer 的传播强度

N=512、z=100 mm 时 ASM、Fresnel 与 Fraunhofer 的传播强度

图 1 比较了三个代表性距离下的三种传播模型。每一列分别对应 ASM、Fresnel 和 Fraunhofer；三行依次表示近场（$z=1,\mathrm{mm}$）、过渡区（$z=10,\mathrm{mm}$）和远场（$z=100,\mathrm{mm}$）。

在 $z=1,\mathrm{mm}$ 时，ASM 和 Fresnel 给出了几乎相同的紧凑衍射图样，而 Fraunhofer 模型本来就不应该在近场成立。在 $z=10,\mathrm{mm}$ 时，ASM 和 Fresnel 光场开始扩展，并出现更清楚的衍射结构，说明系统正在向远场过渡。

到了 $z=100,\mathrm{mm}$，ASM 和 Fresnel 光场已经扩展到几乎填满固定计算窗口。图样在边界附近出现截断，整个光场中还出现了额外的网格状数值结构。预期衍射尺度已经与数值窗口相当，原来的网格无法再可靠表示传播后的光场。

Fraunhofer 面板绘制在它自己的原生输出坐标上：

$$ x_{\mathrm{out}}=\lambda z f_x. $$

因此，不能把它在图上的空间范围，直接与 ASM 和传递函数 Fresnel 所使用的固定输出窗口比较。

这些现象说明，远场中的异常行为不只来自模型近似，也来自数值表示本身。因此，需要进一步检查计算窗口、空间频率采样和可能的 FFT 边界伪影。

N=512 时 ASM 与 Fresnel、Fresnel 与 Fraunhofer 的 MSE 随传播距离变化

图 2 展示了原始 $N=512$ 网格下，ASM–Fresnel 与 Fresnel–Fraunhofer 归一化中心强度剖面的 MSE。Fresnel–Fraunhofer 误差在 $z=100,\mathrm{mm}$ 时反常回升，说明数值限制已经开始污染远场比较。

图 2 给出了三种传播模型产生的归一化中心强度剖面之间的 MSE。ASM 与 Fresnel 在三个距离上都很接近，不过它们的差异在 $z=100,\mathrm{mm}$ 时略有增大。

Fresnel 与 Fraunhofer 的比较呈现了一个更可疑的趋势：误差在 $z=1,\mathrm{mm}$ 时很大，在 $z=10,\mathrm{mm}$ 时降到接近零，却在 $z=100,\mathrm{mm}$ 时再次升高。近场误差大是合理的，因为 Fraunhofer 近似在那里并不成立；但到了 $100,\mathrm{mm}$，误差反而上升就不合理了。随着传播进入远场，Fraunhofer 衍射应该逐渐接近 Fresnel，而不是偏离它。

结合图 1 中被截断的光场，这次回升说明大传播距离处的误差主要由数值表示主导，而不是物理近似本身。于是，下一步应该检查计算窗口和空间频率采样。

2. 诊断问题

图 1 和图 2 说明，传播方程本身不是唯一误差来源。接下来要问的就不再是“哪一种传播模型错了”，而是：当前数值网格是否有能力表示传播后的光场？

2.1 比较预期物理尺度与数值窗口

对于这个实验中的方形孔径，沿一个方向的第一个 Fraunhofer 零点大约位于

$$ x_1=\frac{\lambda z}{a}. $$

当

$$ z=100\,\mathrm{mm}, $$

且

$$ a=100\,\mu m, \qquad \lambda=532\,nm, $$

可得

$$ x_1\approx0.532\,\mathrm{mm}. $$

原始仿真使用

$$ N=512, \qquad \Delta x=2\,\mu m, $$

所以总计算窗口为

$$ L=N\Delta x=1.024\,\mathrm{mm}, $$

对应范围约为

$$ [-0.512,\,+0.512],\mathrm{mm}. $$

也就是说，预期的第一个衍射零点已经非常接近数值窗口边界。

这与图 1 中的视觉表现一致：在 $z=100,\mathrm{mm}$ 时，ASM 和 Fresnel 光场几乎填满了整个计算域。

2.2 数值参数彼此耦合

数值网格由几个相互关联的量决定：

$$ L=N\Delta x, $$

$$ \Delta f=\frac{1}{N\Delta x}, $$

以及

$$ f_{\mathrm{Nyquist}} = \frac{1}{2\Delta x}. $$

这些量不能彼此独立地选择。

当保持 $\Delta x$ 不变并增大 $N$ 时，

$$ N\uparrow \quad\Rightarrow\quad L\uparrow, $$

同时

$$ N\uparrow \quad\Rightarrow\quad \Delta f\downarrow. $$

因此，仿真会同时获得更大的物理窗口和更细的空间频率网格。

最大可表示空间频率

$$ f_{\mathrm{Nyquist}} $$

不会变化，因为它只由 $\Delta x$ 决定。

所以，增大 $N$ 不只是“增加像素数量”，而是在改变同一个物理系统的数值表示方式。

2.3 还必须正确理解输出坐标

三种传播方法不一定会在相同的原生输出坐标上产生光场。

对于 ASM 和传递函数 Fresnel 传播，

$$ \Delta x_{\mathrm{out}} = \Delta x_{\mathrm{in}}. $$

Fraunhofer 衍射则通过下面的关系把空间频率映射到物理位置：

$$ x_{\mathrm{out}} = \lambda z f_x. $$

因此，即使两个数组形状相同，也可能对应完全不同的物理视场。

在上面的比较中，我先把 Fraunhofer 光场转换到了它的物理输出坐标，再比较中心剖面。这样可以避免直接比较那些代表不同物理位置的像素索引。

但坐标一致，并不代表底层光场一定得到了充分采样。

3. 一个简单的收敛性测试

我没有立刻修改传播算法，而是在只改变采样数的情况下重复了同一个实验：

$$ N=512 \rightarrow N=1024. $$

以下物理参数保持不变：

$$ \Delta x, \qquad \lambda, \qquad z, \qquad a. $$

这一点非常重要：改变 $N$ 并不会改变物理系统处于 Fresnel 区还是 Fraunhofer 区。传播区域由孔径尺寸、波长和传播距离等物理量决定。

改变的是数值表示。

对于 $N=512$，

$$ L=1.024\,\mathrm{mm}, $$

而对于 $N=1024$，

$$ L=2.048\,\mathrm{mm}. $$

与此同时，

$$ \Delta f $$

减小了一半。

N=1024、z=1 mm 时 ASM、Fresnel 与 Fraunhofer 的传播强度

N=1024、z=10 mm 时 ASM、Fresnel 与 Fraunhofer 的传播强度

N=1024、z=100 mm 时 ASM、Fresnel 与 Fraunhofer 的传播强度

图 3 在保持 $\Delta x$、波长、孔径尺寸和传播距离不变的情况下，把网格从 $N=512$ 增大到 $N=1024$。在 $z=100,\mathrm{mm}$ 时，ASM 和 Fresnel 光场已经能够完整落在计算窗口内，两者的一致性也明显好于图 1。

图 3 展示了 $N=1024$ 时的同一组传播实验。与图 1 相比，最明显的改进出现在 $z=100,\mathrm{mm}$：ASM 和 Fresnel 光场现在可以舒适地落在扩大的计算窗口内，而不是一直延伸到边界。它们的整体衍射包络也非常接近，只剩下较弱的网格状结构。

在 $z=1,\mathrm{mm}$ 和 $z=10,\mathrm{mm}$ 时，ASM 与 Fresnel 仍像原始仿真一样几乎完全一致。主要变化发生在最大传播距离处，也就是原始 $N=512$ 网格无法可靠表示扩展光场的位置。

把 $N$ 从 512 增大到 1024 并没有改变物理传播区域，而是把计算窗口扩大了一倍：

$$ L:1.024\,\mathrm{mm}\rightarrow2.048\,\mathrm{mm}, $$

同时也让空间频率间隔

$$ \Delta f=\frac{1}{N\Delta x} $$

减小了一半。结果得到改善，说明图 1 中异常的远场行为强烈依赖数值网格。

4. 修正后的比较

接下来，我使用更大的网格重新计算中心强度剖面。

N=1024 时三个传播距离下 ASM、Fresnel 与 Fraunhofer 的中心强度剖面

图 4 是 $N=1024$ 时 ASM、Fresnel 和 Fraunhofer 传播的归一化中心强度剖面。从左到右，传播距离从 $1,\mathrm{mm}$ 增大到 $100,\mathrm{mm}$。Fraunhofer 衍射在近场差异很大，但随着系统进入远场，它会逐渐接近 Fresnel 结果。

图 4 比较了网格增大到 $N=1024$ 后的归一化中心强度剖面。三个面板分别对应 $z=1,\mathrm{mm}$、$10,\mathrm{mm}$ 和 $100,\mathrm{mm}$。

在 $z=1,\mathrm{mm}$ 时，ASM 和 Fresnel 仍然很接近，而 Fraunhofer 剖面无论宽度还是结构都明显不同。这是预期结果，因为远场近似在这么短的传播距离下并不成立。

在 $z=10,\mathrm{mm}$ 时，三条剖面开始相互接近。ASM 与 Fresnel 依然高度一致；Fraunhofer 已经抓住主要的中心特征，但还无法复现所有 Fresnel 衍射结构。

在 $z=100,\mathrm{mm}$ 时，三种方法的整体包络已经非常相似。尤其是 Fraunhofer 剖面已经紧跟 Fresnel 包络，这与“Fraunhofer 衍射是 Fresnel 衍射的远场极限”一致。ASM 和 Fresnel 曲线在平滑包络周围仍有少量高频波动，说明残余的数值离散效应依然存在。

从左到右的变化呈现了预期的物理关系：

$$ \text{ASM}\approx\text{Fresnel} \quad\text{并且随着 }z\text{ 增大，}\quad \text{Fraunhofer}\rightarrow\text{Fresnel}. $$

与原始 $N=512$ 仿真相比，修正后的剖面进一步说明：$z=100,\mathrm{mm}$ 处的大偏差在很大程度上来自数值网格，而不是传播模型本身。

N=512 时三个传播距离下 ASM、Fresnel 与 Fraunhofer 的中心强度剖面

图 5 是 $N=512$ 时 ASM、Fresnel 和 Fraunhofer 传播的归一化中心强度剖面。从左到右，传播距离从 $1,\mathrm{mm}$ 增大到 $100,\mathrm{mm}$。

对应的 MSE 比较也变得更符合物理预期。

N=1024 时 ASM 与 Fresnel、Fresnel 与 Fraunhofer 的 MSE 随传播距离变化

图 6 展示了 $N=1024$ 时归一化中心强度剖面的 MSE。ASM 与 Fresnel 在测试距离上始终高度一致；随着光场接近远场，Fresnel–Fraunhofer 的差异显著减小。$z=100,\mathrm{mm}$ 处仍有轻微回升，但远低于原始 $N=512$ 网格的结果。

图 6 展示了网格增大到 $N=1024$ 后，归一化中心强度剖面之间的 MSE。

在三个传播距离上，ASM–Fresnel 误差始终很小，不过在 $z=100,\mathrm{mm}$ 时略有增加。这与图 4 中两条剖面的高度一致相符，最长传播距离处只剩少量数值波动。

Fresnel–Fraunhofer 的差异对传播距离更敏感。在 $z=1,\mathrm{mm}$ 时，MSE 约为 $0.125$，反映了远场近似在近场失效；在 $z=10,\mathrm{mm}$ 时，随着 Fraunhofer 剖面开始接近 Fresnel，误差下降到接近零。

在 $z=100,\mathrm{mm}$ 时，误差再次轻微升高，但仍远小于原始 $N=512$ 仿真。结合图 4 的剖面比较，这次较小的回升更像是残余数值离散误差，而不是 Fraunhofer 近似失效。

最重要的是，增大 $N$ 显著降低了原始网格中的远场大误差。这确认了之前的误差回升强烈依赖数值表示，而不是衍射模型固有的特征。

5. 这个测试究竟证明了什么

人们很容易把原始问题归因于某一种具体数值效应。但在保持 $\Delta x$ 不变的情况下增大 $N$，会同时改变两个量：

$$ L\uparrow $$

和

$$ \Delta f\downarrow. $$

因此，仅凭这次实验还无法区分主要误差究竟来自有限计算窗口、空间频率网格分辨率不足，还是两者共同作用。

但它确立了一个更基础的事实：

$$ \text{原始仿真尚未达到数值收敛。} $$

如果合理改变采样网格后，物理结果发生显著变化，那么这个数值结果就不应该被相信。

这也自然引出了 sanity check 的必要性。

6. 仿真前应该检查什么

数值传播实验应该先从几个物理尺度与数值尺度检查开始。$N$、$\Delta x$、孔径尺寸、波长和传播距离彼此耦合，不能分别孤立地选择。

6.1 空间采样

计算窗口的物理尺寸为

$$ L=N\Delta x. $$

第一个问题是：输入光场本身是否在这个网格上得到了充分表示？

例如，一个 $100,\mu m$ 的孔径使用

$$ \Delta x=2\,\mu m $$

采样时，宽度方向大约包含 50 个采样点，足以表达它的基本几何形状。

对于更复杂的光场，快速变化的振幅或相位可能需要更细的采样。

最基础的几个问题是：

最小空间特征是否由足够多的像素表示？
输入光场是否只占据计算窗口中合理的一部分？
是否为光场传播后的扩展保留了足够的空白区域？

6.2 计算窗口

总视场为

$$ L=N\Delta x. $$

应该把它与传播后光场的预期物理尺寸进行比较。

对于宽度为 $a$ 的方形孔径，第一个 Fraunhofer 零点约为

$$ x_1\approx\frac{\lambda z}{a}. $$

如果 $x_1$ 已经与 $L/2$ 相当，那么即使是中心衍射结构也正在逼近数值边界。

这正是原始 $N=512$ 仿真在 $z=100,\mathrm{mm}$ 时发生的情况。

6.3 空间频率采样

傅里叶频率间隔为

$$ \Delta f=\frac{1}{N\Delta x}, $$

Nyquist 极限为

$$ f_{\mathrm{Nyquist}} = \frac{1}{2\Delta x}. $$

这两个量描述了频率网格的不同属性：

$$ \Delta f $$

决定空间频率域的采样精细程度，而

$$ f_{\mathrm{Nyquist}} $$

决定可以表示的最大空间频率。

因此，固定 $\Delta x$ 时，

$$ N\uparrow $$

会带来更细的频率分辨率；而

$$ \Delta x\downarrow $$

会扩展可表示的频率带宽。

两者都可能影响波传播仿真。

6.4 传播区域

传播距离也应该依据物理尺度选择，而不是任意指定。

对于横向特征宽度为 $D$ 的孔径，一个有用的衍射距离为

$$ z_{\mathrm{diff}} \sim \frac{D^2}{\lambda}. $$

它给出了区分不同传播区域的粗略尺度：

$$ z\ll z_{\mathrm{diff}} $$

对应相对近场传播；

$$ z\sim z_{\mathrm{diff}} $$

描述过渡区域；

$$ z\gg z_{\mathrm{diff}} $$

则逐渐接近 Fraunhofer 远场。

对于当前实验，

$$ D=100\,\mu m, \qquad \lambda=532\,nm, $$

因此

$$ z_{\mathrm{diff}} \approx \frac{(100\times10^{-6})^2}{532\times10^{-9}} \approx 18.8\,\mathrm{mm}. $$

所以，三个传播距离是有意选择的：

$$ 1\,\mathrm{mm} \ll 18.8\,\mathrm{mm}, $$

$$ 10\,\mathrm{mm} \sim 18.8\,\mathrm{mm}, $$

以及

$$ 100\,\mathrm{mm} \gg 18.8\,\mathrm{mm}. $$

它们大致代表近场、过渡区和远场传播。

这让比较具有了物理含义：实验并不是随意测试三个 $z$ 值，而是在追踪衍射图样如何跨越不同传播区域。

6.5 输出坐标

不同的数值传播方法可能使用不同的原生输出网格。

对于 ASM 和传递函数 Fresnel 传播，

$$ \Delta x_{\mathrm{out}} = \Delta x_{\mathrm{in}}, $$

而 Fraunhofer 衍射通过下面的关系把空间频率映射为物理位置：

$$ x_{\mathrm{out}} = \lambda z f_x. $$

因此，在比较两个传播光场之前，应该始终检查：

它们是否覆盖同一个物理区域？
它们是否具有相同的物理采样间隔？
如果没有，是否已经把其中一个光场插值到共同的物理坐标网格？

数组形状相同，并不足以说明它们可以直接比较。

6.6 数值收敛

即使所有解析估算看起来都合理，最终结果仍然应该经过数值测试。

一种简单的收敛性测试，是在改变一个数值参数后重复仿真，例如

$$ N=512,\quad1024,\quad2048, $$

然后在同一个物理区域内比较光场。

可靠的仿真最终应该近似满足

$$ I_{N=1024}(x) \approx I_{N=2048}(x). $$

如果改变数值网格后，传播光场仍然发生显著变化，那么仿真尚未收敛。

理想情况下，还应该分别改变 $N$ 和 $\Delta x$，因为它们控制不同的数值属性。

6.7 基础物理一致性检查

还有一些简单检查可以帮助发现实现错误：

当传播距离 $z=0$ 时，是否能复现输入光场？
对称输入是否产生预期的对称衍射图样？
在存在已知解析解时，结果是否与解析解一致？
光场是否正在接近数值边界？
强度归一化和物理单位是否处理一致？
在有效范围重叠的区域，定性上不同的传播方法是否相互一致？

与解释一个错误仿真相比，这些检查的成本很低。

因此，一套实用的 sanity check 顺序是：

$$ \text{输入采样} \rightarrow \text{窗口尺寸} \rightarrow \text{频率网格} \rightarrow \text{传播区域} \rightarrow \text{输出坐标} \rightarrow \text{收敛性测试} $$

目的并不是找到一套普遍正确的数值参数，而是确认当前网格确实有能力表示眼前这个具体的光学问题。

本文采用署名-非商业性使用-相同方式共享 4.0 国际许可协议，转载请注明出处。