做相干光成像计算的人，十有八九都遇到过这种诡异时刻：明明菲涅耳衍射积分公式背得滚瓜烂熟，FFT也老老实实调了，出来的图却带一圈圈明暗环纹，或者光斑半径越传播越小，甚至模拟高斯光束自由传播都会算出“空心光束”。我早先第一次上手相干光成像的精确计算时，就因为这些现象整整排查了一晚上，最后发现根本不是代码写错，而是对离散衍射的物理尺度理解偏了。
相干光成像的精确计算是计算全息、相位恢复、光束整形、数字全息这些方向的共同地基。它不像几何光学那样画几条光线就能完事，每一步离散化都藏着尺度换算的坑：采样间隔、衍射距离、像素数量三者是互相咬死的，一个对不上，算出来的“图像”就只是数学伪影。这篇文章把我重新梳理整套相干光成像精确计算流程时的完整思路写出来，包含三种可落地的衍射计算路线、一组动手前必须算清的采样约束公式、几个能直接抄的代码实例，以及我一路踩过来的排错经验。无论你是刚进光学仿真方向的研究生，还是在做全息显示、激光加工的工程师，这份流程都能帮你把“算不准”变成“算得可控”。
1. 一动手就翻车：相干光成像计算里的三类隐藏问题
1.1 连续积分变成离散求和后，世界就不一样了
相干光成像的传播计算，核心是标量衍射理论。最常见的菲涅耳衍射积分长这样：
$$U(x,y,z)=\frac{e^{ikz}}{i\lambda z}\iint U_0(x_0,y_0)\exp\left[\frac{ik}{2z}\left((x-x_0)^2+(y-y_0)^2\right)\right]dx_0dy_0$$
看着很规整，但计算机没法处理连续积分，所有东西都要离散成网格阵列。这一离散，三件事同时发生：采样间隔、观察窗口、频域范围被绑在了一起。
用生活化的方式理解：连续积分就像用眼睛直接看一幅画的全部细节，离散化则是拿一张网孔固定的纱窗挡住眼睛看画。纱窗网孔太粗，你会漏掉画里的细条纹；纱窗尺寸太小，你只能看到画的一个角落。很多计算翻车，不是因为“画”错了，而是纱窗本身不合适。菲涅耳积分里那个二次相位项$\exp[\frac{ik}{2z}(x^2+y^2)]$，在网格边缘变化非常快，如果采样间隔太大，相邻像素之间的相位差就会超过$\pi$，直接造成混叠。这一点在近场计算里尤其致命。
1.2 最容易踩的三个误区
误区一：以为FFT之后取模就是衍射场。很多人写代码时直接把$U_0$丢进fft2，再乘个大常数就完事了。这是夫琅禾费远场的做法，在菲涅耳区是错的。菲涅耳衍射在FFT之前要先乘一个输入面上的二次相位因子，FFT之后还要再乘一个输出面上的二次相位因子，少了任何一个chirp项，相位信息都是错的。
误区二：以为输入输出网格必须一致。FFT本身不改变数据尺寸，但菲涅耳衍射的物理坐标在离散后会悄悄变化。单次FFT算法的输出采样间隔是$\Delta x_2 = \lambda z / (N \Delta x_1)$，和距离$z$强相关。距离一变，输出平面的像素尺寸就变，观察窗口也变。我那次“光斑越传越小”的诡异现象，根源就在这里：$z$越小，输出窗口被压缩得越厉害，光束明明没聚焦，却因为窗口变小而看起来像被“掐”成一个小亮斑。
误区三：以为网格越大越精确。把$N$从512加到2048确实能提高分辨率，但如果物理窗口大小不变，像素尺寸变小，频域范围变大，计算量暴涨的同时，欠采样的区域可能还是一样欠采样。近场问题尤其如此，光靠加$N$解决不了chirp项振荡过快的问题，需要配合换算法或零填充策略。
2. 三条可落地的衍射计算路线：单次FFT、角谱与卷积
2.1 单次FFT：快，但输出网格会“变”
单次FFT算法是将菲涅耳积分直接改写成傅里叶变换形式：
$$U(x,y,z)=\frac{e^{ikz}}{i\lambda z}\exp\left[\frac{ik}{2z}(x^2+y^2)\right]\cdot \mathcal{F}\left[U_0(x_0,y_0)\exp\left(\frac{ik}{2z}(x_0^2+y_0^2)\right)\right]$$
实现起来非常快，整个计算只有一次FFT。核心步骤是：
- 构造输入面坐标网格，乘以输入chirp因子；
- 做二维FFT；
- 在频域坐标$f_x = x/(\lambda z)$处读取结果，乘以输出chirp因子和常数$\frac{e^{ikz}}{i\lambda z}$。
这里的坑在于：FFT输出数组的顺序对应频率从0到N-1，而我们习惯把坐标原点放在数组中心，所以必须用ifftshift再fft2，得到频谱后再fftshift回到中心。很多人写代码时省掉这一步，结果相位的空间分布整个颠倒，光斑自然对不上。
这个方法适合中等以上距离。它的最大限制就是输出采样间隔被$\lambda z/(N\Delta x_1)$锁定，无法和输入网格保持一致。
2.2 角谱法：近场与网格不变的首选
角谱法不经过菲涅耳近似，属于严格的瑞利-索末菲衍射的频域实现：
$$U(x,y,z)=\mathcal{F}^{-1}\left[\mathcal{F}(U_0)\cdot H(f_x,f_y)\right]$$
其中传递函数为：
$$H(f_x,f_y)=\exp\left[ikz\sqrt{1-(\lambda f_x)^2-(\lambda f_y)^2}\right]$$
这个方法的优点是输出网格和输入网格完全相同，观察窗口不会随距离变化。因此在近场模拟、需要网格严格对齐的迭代算法里非常好用。
注意两个细节。第一，当$(\lambda f_x)^2+(\lambda f_y)^2 > 1$时，对应的是倏逝波分量，在自由空间传播会迅速衰减，数值上直接将该处传递函数置零即可。第二，它的适用距离也不是无限的。距离$z$很大时，传递函数的相位在频域高频端振荡非常剧烈，如果采样点数不足，同样会有欠采样问题。粗略的约束条件是$z$不要超过约$2N\Delta x_1^2/\lambda$这个量级，否则需要增大$N$。
2.3 卷积型方法：通用解法与循环卷积警告
菲涅耳衍射本质上是一个卷积：
$$U(x,y,z)=U_0(x,y)\otimes h(x,y),\quad h(x,y)=\frac{e^{ikz}}{i\lambda z}\exp\left[\frac{ik}{2z}(x^2+y^2)\right]$$
用离散FFT实现卷积，是线性卷积还是循环卷积，这是个必须当面问清楚的问题。直接对$U_0$和$h$分别做FFT、相乘、逆FFT，得到的是循环卷积，边缘的部分会“卷”到对面去，产生周期条纹。解决办法是零填充：把$U_0$和$h$都补零到2N×2N再卷积，最后裁出中间N×N区域，这样才等价于线性卷积。
卷积法的优点是输出网格不变，且不像角谱法那样存在$z$过大时传递函数剧烈振荡的局限，因为脉冲响应$h$可以预先解析构造。代价是计算量更大，一次传播需要三次FFT（$U_0$一次、$h$一次、结果逆变换一次），内存占用也更多。好在工程中常常可以把$h$的FFT提前算好缓存，实际运行就退化到两次FFT的量级。
2.4 三条路线怎么选：一张对比表
| 方法 | FFT次数 | 输出网格 | 适用距离 | 实现难点 | 推荐场景 |
|---|---|---|---|---|---|
| 单次FFT | 1次 | 随距离缩放，$\Delta x_2=\lambda z/(N\Delta x_1)$ | 中远距离 | 坐标映射、shift顺序 | 快速模拟、实时计算 |
| 角谱法 | 2次 | 与输入完全一致 | 近场、中距离 | 传递函数相位处理 | GS迭代、相位恢复、网格对齐要求高的场合 |
| 卷积法 | 3次 | 与输入完全一致 | 适用范围广 | 循环卷积需补零 | 通用验证、中等距离精确计算 |
我的习惯是：近场用角谱法，远场用单次FFT，拿不准的时候用卷积法做交叉验证。三条路线在参数合理的情况下应该给出几乎一致的结果，如果差异明显，那大概率是采样或边界出了问题，而不是算法本身不对。
3. 动手前必须算清的采样账：网格、距离与像素的三角关系
3.1 采样间隔变形的根源：一对耦合公式
单次FFT算法里的输出采样间隔：
$$\Delta x_2 = \frac{\lambda z}{N \Delta x_1}$$
这可能是整个相干光成像数值计算里最容易被忽略、却影响最大的一个公式。它告诉我们：输入像素尺寸、像素数量、波长、传播距离四个量，共同决定了输出平面的观察窗口。
举一个我常用来和学生算的实例。设$\lambda = 632.8\text{nm}$，$\Delta x_1 = 10\mu\text{m}$，$N=512$。则输入窗口$L_1 = N\Delta x_1 = 5.12\text{mm}$。当传播距离$z=1\text{cm}$时，输出采样间隔：
$$\Delta x_2 = \frac{632.8\times 10^{-9}\times 0.01}{512\times 10\times 10^{-6}} \approx 1.236\mu\text{m}$$
输出窗口$L_2 = 512 \times 1.236\mu\text{m} \approx 0.633\text{mm}$。原本5.12mm的光束，在输出平面看到的“视野”只剩下0.633mm，大部分光强被挤出窗口，模拟结果当然是一片亮斑加生硬截断。而当$z=50\text{cm}$时，$\Delta x_2 \approx 61.8\mu\text{m}$，$L_2$膨胀到约31.6mm，光束又显得非常小。
这就是很多人模拟高斯光束传播时“光斑尺寸忽大忽小”的真相：不是光斑变了，是观察窗口在变。
3.2 临界距离：让输入输出窗口一致的临界面
如果希望输出窗口$L_2$和输入窗口$L_1$相等，即$N\Delta x_2 = N\Delta x_1$，代入上面的耦合公式，得到：
$$z_c = \frac{N \Delta x_1^2}{\lambda}$$
这个$z_c$我习惯叫它临界距离。用上面的参数计算：
$$z_c = \frac{512\times (10\times 10^{-6})^2}{632.8\times 10^{-9}} \approx 0.0809\text{m} \approx 8.09\text{cm}$$
意思是，只有传播距离在8cm附近时，单次FFT算法的输入输出窗口才恰好一致。低于这个距离，输出窗口收缩；高于这个距离，输出窗口扩张。
明白了这个关系，很多“bug”都能解释。比如要模拟一个束腰半径1mm的高斯光束传播到$z=0.3\text{m}$处，用单次FFT且保持输入窗口5.12mm，输出窗口会膨胀到约19mm，如果只盯着中心区域看，光斑会显得比实际的小。正确的做法是：要么主动用临界距离反推合适的输入窗口，要么干脆换角谱法，让网格尺寸不随距离漂移。
动手写代码前，我强烈建议先画一张参数表：波长、像素尺寸、像素数、输入窗口、输出窗口、传播距离、临界距离、菲涅耳数。把这些数字填满，再谈写代码的事。
3.3 判断“算得对不对”的三种方法
第一种方法，和解析解对比。高斯光束的自由传播有解析公式，束腰半径$w_0$、波长$\lambda$的高斯光束传播到$z$处，光斑半径满足：
$$w(z)=w_0\sqrt{1+\left(\frac{z\lambda}{\pi w_0^2}\right)^2}$$
用数值算法算出二维强度分布，再用二阶矩法提取光斑半径，和这个公式对比。误差在1%以内，说明你的传播算法和参数设置是可信的。这个验证手段是整套流程的“体检项目”，建议每次换参数后都跑一遍。
第二种方法，能量守恒检查。不考虑边界截断的理想情况下，衍射传播前后总能量应该守恒。算出$\sum |U|^2 \Delta x^2$，传播前后的值应该几乎相等。如果能量掉了20%，说明窗口截断严重；如果能量多了，说明出现了数值混叠。
第三种方法，多算法互验。同一个物理场景，分别用角谱法和卷积法计算，对比结果。两种算法实现思路差别很大，如果它们对不上，几乎可以断定是采样参数或边界条件的问题，而不是巧合。
4. 从模拟到全息图：三个高频场景的完整实操
4.1 高斯光束传播：跑通任何算法的第一关
我习惯用高斯光束验证每个传播函数。下面这段代码用角谱法模拟束腰半径$w_0=0.5\text{mm}$的高斯光束传播到10cm处，再和解析光斑半径对比。
import numpy as np def gaussian_beam(dx, N, w0): """束腰在z=0的高斯光束复振幅，无初始二次相位。""" x = (np.arange(N) - N // 2) * dx X, Y = np.meshgrid(x, x) return np.exp(-(X**2 + Y**2) / w0**2) def prop_asm(u0, lam, dz, dx): """角谱法衍射传播，输入输出网格一致。""" ny, nx = u0.shape k = 2 * np.pi / lam fx = np.fft.fftfreq(nx, d=dx) fy = np.fft.fftfreq(ny, d=dx) FX, FY = np.meshgrid(fx, fy) with np.errstate(invalid='ignore'): H = np.exp(1j * k * dz * np.sqrt(1 - (lam * FX)**2 - (lam * FY)**2)) H[(lam * FX)**2 + (lam * FY)**2 > 1] = 0 # 倏逝波 U0 = np.fft.ifftshift(u0) U1 = np.fft.fft2(U0) * H u1 = np.fft.fftshift(np.fft.ifft2(U1)) return u1 # 参数 lam = 632.8e-9 dx = 10e-6 N = 512 w0 = 0.5e-3 z = 0.1 u0 = gaussian_beam(dx, N, w0) u1 = prop_asm(u0, lam, z, dx) # 解析光斑半径 w_analytic = w0 * np.sqrt(1 + (z * lam / (np.pi * w0**2))**2) # 数值光斑半径（二阶矩法） I1 = np.abs(u1)**2 x = (np.arange(N) - N // 2) * dx X, Y = np.meshgrid(x, x) r2 = X**2 + Y**2 w_numeric = np.sqrt(2 * np.sum(r2 * I1) / np.sum(I1)) print(f"解析光斑半径: {w_analytic*1e3:.4f} mm") print(f"数值光斑半径: {w_numeric*1e3:.4f} mm")用这个参数跑，两条结果误差通常不到百分之一。如果这一步就对不上，不要急着做全息，先把采样参数表拿出来一项项核对。
注意，测试时把目标图像的尺寸、位置和能量都先归一化，避免后续迭代时被极端值带偏。
4.2 计算全息图：GS算法的工程化实现
做计算全息，最常见的任务是给定目标强度图像，反演出纯相位全息图。盖师贝格-撒克斯通（GS）算法是入门的首选：物面振幅固定为激光照明振幅，像面振幅替换为目标振幅，在物面和像面之间反复正逆传播迭代。
我这里直接复用上面的prop_asm，写出一个最小可运行的GS循环：
def gs_phase(target_amp, lam, z, dx, n_iter=100): """GS算法生成纯相位全息图，物面振幅为均匀平面波。""" ny, nx = target_amp.shape phase = np.random.rand(ny, nx) * 2 * np.pi field = np.ones((ny, nx)) * np.exp(1j * phase) # 物面振幅取1 for _ in range(n_iter): img = prop_asm(field, lam, z, dx) # 正传 amp = np.abs(img) phi_img = np.angle(img) img_corrected = target_amp * np.exp(1j * phi_img) # 像面振幅约束 field_back = prop_asm(img_corrected, lam, -z, dx) # 反传 phi_obj = np.angle(field_back) field = np.ones((ny, nx)) * np.exp(1j * phi_obj) # 物面振幅约束 return np.angle(field)几个工程细节：
第一，目标振幅target_amp一定要和光源总能量匹配。GS算法把物面振幅固定为1，目标振幅的总能量应该等于$N^2$这个量级。如果目标图像本身特别亮或特别暗，迭代会震荡，收敛很慢。最简单的做法是把目标振幅归一化为$\text{target_amp} = \sqrt{I / \bar{I}}$，让平均振幅接近1。
第二，初始相位用随机相位，不要用全零或常数相位，否则算法几乎必然陷入一个很差局部极小值。多跑几次取重建误差最小的一次，是实用中的常见策略。
第三，如果全息图直接重建时中心出现一个特别亮的零级光斑，通常是因为相位分布含有大块的低频分量或直流偏置。可以在目标图像后面乘一个离轴载频相位，把有效信息从零级衍射中分离出来。这个技巧在做显示类全息图时几乎是必须的。
4.3 光束整形与相位恢复：迭代算法的实战细节
光束整形是另一个高频场景，典型任务是把高斯光束变成平顶光束、环形光束或任意自定义形状。实现思路和GS基本一致，只是目标振幅从“图像”变成了“想要的远场或焦面光强分布”。
这里最容易踩的坑是目标轮廓设计得太“硬”。比如要一个矩形平顶光斑，直接在目标平面上给一个理想的矩形函数，迭代结果边缘会出现明显的振铃振荡。原因是硬边界的傅里叶频谱衰减很慢，算法很难在有限迭代次数内兼顾边缘陡峭和带内均匀。
解决办法是给目标轮廓加过渡带，比如用超高斯函数：
$$A_{\text{target}}(r)=\exp\left[-\left(\frac{r}{R}\right)^{n}\right]$$
其中$n$取4到8，而不是直接设0/1二值。这样迭代得到的相位不仅收敛更快，整形效果也更好。我试过$n=2$改成$n=6$之后，边缘波纹从15%降到了3%以内。
相位恢复方面，经典的GS算法有个老大难问题：容易陷入局部极小。对同一个目标图像，稍微改一下初始随机种子，迭代结果可能差出一大截。工程上的改良思路有两个方向。一个是多初始值：跑20次GS，每次用不同随机种子，取重建误差最小的结果。另一个是改用混合输入输出（HIO）算法，在像面约束里引入一个负反馈项，破坏局部极小的吸引域。对于大多数简单目标，多跑几次GS已经够用。
5. 现场排错实录：条纹、光斑失真与近场发散
5.1 条纹状伪影与网格纹“查不完”
最经典的故障画面是：重建图像背景里出现一排排规则明暗条纹，或者边缘有鳞片状网格纹，移动几张图像后条纹方向还会变。
这类问题的根源几乎都是周期延拓。FFT假定输入信号是无限周期重复的，如果计算窗口边缘的光场幅度没有衰减到0，边缘的不连续就会在频谱里产生高频扩散，表现成周期性条纹。
处理办法按优先级排：
第一，扩大计算窗口。把物场零填充到2N甚至4N大小，传播完成后只取中心N×N区域。这个方法能同时改善边缘不连续和高频采样密度，是我最常用的手段。
第二，如果零填充后内存压力大，也可以用边缘窗函数，比如二维汉宁窗，让窗口边缘平滑衰减到0。但要注意，加窗会把窗口内的有效信息也削弱，只适合边缘区域本来就接近0的场景。
第三，检查频谱里有没有特别亮的直流峰。如果有，说明光场里叠加了一个常数偏置。处理方式是在输入面减掉均值，或者改用离轴载频把信号从频域中心移开。
5.2 光斑尺寸总对不上理论值
这个问题的排查路径非常固定，按顺序查：
先打印输出平面的实际采样间隔。如果用单次FFT，确认输出坐标是用$\Delta x_2 = \lambda z/(N\Delta x_1)$重新生成的，而不是沿用了输入平面的$\Delta x_1$。很多人栽在这里。
再查传播距离和临界距离的关系。如果$z$远小于$z_c$，单次FFT的输出窗口会缩得很小；如果$z$远大于$z_c$，窗口大得离谱。此时应该换用角谱法或卷积法，而不是继续硬扛。
还要确认输入面上的相位条件。高斯光束的公式要求束腰在输入面$z=0$处。如果输入场带有一个曲率半径$R$的二次相位，传播到后面的光斑半径表达式会多一项，不能直接套用$w(z)=w_0\sqrt{1+(z\lambda/\pi w_0^2)^2}$。处理方法是先用数字方式“反推”束腰位置，或者干脆用完整的$q$参数法解析式。
5.3 近场与极近场发散
算近场衍射时经常遇到这个问题：距离$z$只有几个毫米甚至微米，角谱法的结果出现高频噪声，单次FFT的结果则完全无法辨认。
原因是菲涅耳近似在极近场条件下误差变大，同时chirp项的空间频率太密，网格根本采不住。这种情况下：
第一选择是用角谱法，因为它是严格亥姆霍兹方程的解，不引入菲涅耳近似。
第二，如果角谱法也有振荡，问题多半是频域采样不足。可以把传递函数在频域补零插值后再乘，等效于提高频域分辨率。也可以用卷积法加零填充，让脉冲响应的相位在相邻像素之间变化平缓。
第三，引入菲涅耳数做预判：
$$F = \frac{a^2}{\lambda z}$$
其中$a$是孔径半径。当$F$远大于1时属于近场，此时优先用角谱法；当$F$接近1或更小时，单次FFT开始变得可靠。这个无量纲数是我做参数表时必填的一栏，它能帮你在一开始就避开近场的地雷。
5.4 我调试衍射计算时固定先做的三件事
经验积累到现在，我每次调试相干光成像计算，都会强迫自己按顺序做三件事。
第一件事，写参数表。波长、像素尺寸、像素数、输入窗口、输出窗口、传播距离、临界距离、菲涅耳数，全部列清楚。参数表写完了，超过一半的“疑难杂症”自己就暴露了。
第二件事，用高斯光束预热。任何新写的传播函数或迭代算法，先拿高斯光束传播和解析解对比。这一步通过了，再碰全息图和相位恢复。不要图省事直接拿目标图像开跑，否则算法错了都不知道错在哪个环节。
第三件事，全程保持复数思维。衍射场是复数，强度只是幅值平方。调试的时候不要只盯着abs()看，要把相位图也打出来看看。相位分布里一旦出现密集的螺旋位错或规则的梳状条纹，那基本就是网格映射或shift顺序出错。
这套流程折腾下来，相干光成像的精确计算会从“玄学”变成“工程”：每个参数量化到位，每条算法路线边界清晰，每个异常都有对应的排查方向。我个人最大的体会是，这类计算不怕算法“不够高级”，就怕尺度关系没理顺就急着上迭代，后者的时间成本往往是前者的好几倍。
