做相干光成像计算的人十有八九都遇到过这种诡异时刻明明菲涅耳衍射积分公式背得滚瓜烂熟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^2y^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^2y^2)\right]\cdot \mathcal{F}\left[U_0(x_0,y_0)\exp\left(\frac{ik}{2z}(x_0^2y_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^2y^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次数输出网格适用距离实现难点推荐场景单次FFT1次随距离缩放$\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}$$N512$。则输入窗口$L_1 N\Delta x_1 5.12\text{mm}$。当传播距离$z1\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大部分光强被挤出窗口模拟结果当然是一片亮斑加生硬截断。而当$z50\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的高斯光束传播到$z0.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_00.5\text{mm}$的高斯光束传播到10cm处再和解析光斑半径对比。import numpy as np def gaussian_beam(dx, N, w0): 束腰在z0的高斯光束复振幅无初始二次相位。 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, ddx) fy np.fft.fftfreq(ny, ddx) FX, FY np.meshgrid(fx, fy) with np.errstate(invalidignore): 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_iter100): 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二值。这样迭代得到的相位不仅收敛更快整形效果也更好。我试过$n2$改成$n6$之后边缘波纹从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$窗口大得离谱。此时应该换用角谱法或卷积法而不是继续硬扛。还要确认输入面上的相位条件。高斯光束的公式要求束腰在输入面$z0$处。如果输入场带有一个曲率半径$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顺序出错。这套流程折腾下来相干光成像的精确计算会从“玄学”变成“工程”每个参数量化到位每条算法路线边界清晰每个异常都有对应的排查方向。我个人最大的体会是这类计算不怕算法“不够高级”就怕尺度关系没理顺就急着上迭代后者的时间成本往往是前者的好几倍。