马赫-曾德干涉仪这个东西光学实验室里见得最多但真正把它吃透的人其实不多。我最早接触它是在做光纤振动传感的时候那会儿被条纹漂移搞得头大后来花了一整个月把波动光学仿真做通才算彻底搞明白两条光路里每一处相位和振幅变化到底是怎么回事。这篇博文就把我从建模思路到仿真实操、再到精度控制和排坑的经验一次性写出来适合正在做干涉仪仿真、光学系统设计或光学测量相关工作的朋友。如果你刚接触这个仪器也没关系这里面的每一步我尽量把“为什么要这么做”也讲清楚。1. 马赫-曾德干涉仪从实验台到仿真屏1.1 这台“光学听力计”到底在测什么马赫-曾德干涉仪Mach-Zehnder interferometer简称MZI是光学干涉家族里非常经典的一员。它的结构一句话就能说清一束光先被分束器一分为二两路光分别走两条独立的臂然后再被第二个分束器合到一块最后形成干涉。这个“先分开、后合拢”的结构让两条光路的物理空间彼此独立所以你可以在一条臂上随意放样品、加相位调制器、做时间延迟而不用怕干扰另一条臂。正是这种灵活性让它成了折射率测量、温度传感、声波检测、气体浓度分析以及量子光学实验中出镜率极高的装置。我举个具体场景你就理解了假设信号臂上有一小段光纤暴露在外外部气体浓度的变化会导致这段光纤的有效折射率发生微小改变光走过的光程就会变最后干涉条纹会跟着移动。只要你把条纹移动量测准反推回去就是折射率变化量。这个过程本质上就是利用“波动的相位”做一把极其灵敏的尺子。之所以强调“波动”是因为干涉现象本身就是相位的直接体现几何光学那套“一条线”的光线模型在这里根本不够用你必须用复振幅和相位来描述整个系统。1.2 为什么必须做波动光学仿真很多刚入门的同学会问我这光路看着这么简单直接套双光束干涉公式不就行了为什么还要仿真答案是公式只能给你理想平面波的结果而真实的光是高斯光束是有束腰半径、有发散角的分束器有尺寸镜面有孔径传播一段距离后还会出现衍射效应。你要是用纯几何方式去算最后得到的干涉条纹可能只在“中心区域”勉强对得上边缘部分几乎全是错误的。波动光学仿真要做的事情就是在计算机里把“复振幅”这个物理量完整地模拟出来一束高斯光出发经过分束器分成两束各自走完几十毫米甚至几百毫米的路径期间经历了衍射、波前畸变、相位延迟最后再次相遇叠加。这个过程尤其适合研究下面这几类问题探测器平面上干涉条纹的对比度受什么影响两条臂的光程差略有偏差时条纹状态如何变化光束尺寸、传播距离对输出端两路光功率分配的影响干涉仪对微小相位扰动的响应灵敏度和线性范围。如果你后续要做实验仿真还能帮你在搭光路之前就把参数定下来知道用多大焦距的透镜、探测器放在什么位置避免一到实验室就来回搬镜架。1.3 仿真能帮你回答的四个典型问题我把自己做仿真时最常遇到的四个问题列在这里后文会逐一展开问题类型仿真作用关键指标条纹质量判断干涉条纹是否清晰有没有变形可见度对比度相位提取从干涉图中反演相位变化换算物理量相位误差灵敏度极限确定能探测到的最小相位扰动相位噪声、探测下限输出串扰两输出端口之间的能量分配波动分光比变化这些问题如果只靠手推公式要么算不了非理想条件要么算出来的结果让你心里没底。而一版靠谱的波动光仿真代码能在几秒钟内给你一组可视化结果还能顺便帮你把数据导出成表格或者图片用来写报告和论文非常方便。2. 波动光学仿真背后的数学工程2.1 不是“画光路图”复振幅才是主角在波动光学框架下光场不再是一堆几何射线而是用一个复数场 (E(x,y,z)) 来描述。它的模表示振幅辐角表示相位。光强则是 (I |E|^2)。干涉的本质就是两个或多个复场叠加后模平方出现的交叉项也就是 (2|E_1||E_2|\cos(\Delta \phi)) 这个来源。拿马赫-曾德干涉仪来说假设分束器分光比是 50:50合束前的两束光复振幅分别是 (E_1) 和 (E_2)那么输出端的光强是[ I_{out} \frac{1}{4}|E_1 E_2|^2 ]展开后就是直流项各自的强度加上交流项干涉项。当 (E_1) 和 (E_2) 的相位差为 0 时亮纹为 (\pi) 时暗纹。这个表达式看似简单但真正的功夫在于如何精确算出 (E_1) 和 (E_2) 在合束时的复振幅分布——因为它们在各自臂上传播时已经发生了衍射和相位累积。这个“算传播”的过程就是波动光学仿真需要解决的核心问题。2.2 角谱法在频域里做衍射自由空间传播的仿真方法有挺多直接对菲涅尔衍射积分做数值计算、用SASD单步角谱衍射方法、或者用光束传播法BPM逐步推进。在我的MZI仿真里最稳也最常用的是角谱法Angular Spectrum Method。它的思想并不复杂把光场从空间域变换到空间频率域每一个频率分量看作一个传播方向确定的平面波在频域里给每个分量乘上一个相位因子再做逆傅里叶变换回到空间域。用公式表示[ E(x,y,z\Delta z) \mathcal{F}^{-1}\left{ \mathcal{F}{E(x,y,z)} \cdot \exp\left[i k \Delta z \sqrt{1 - (\lambda f_x)^2 - (\lambda f_y)^2}\right] \right} ]其中 (f_x, f_y) 是空间频率(k2\pi/\lambda)。这个方法的优点是它对传播距离没有像近轴近似的强限制只要采样满足条件在短距离和中等距离下精度都相当高。而且实现起来就是对数组做两三次FFT计算速度非常快。我做仿真时通常会封装一个propagate_ASM(field, wavelength, delta_x, delta_y, distance)函数只要传入当前光场和传播距离就能输出传播后的光场。这样两臂就可以统一走同一套传播代码唯一区别是臂长和相位延迟参数不同。2.3 采样条件与网格设计角谱法看似只需两行代码但坑非常深最大的坑就是采样条件。因为角谱法相当于把光场分解成一系列平面波而离散傅里叶变换对应的空间频率范围是有限的最高频率由采样间隔决定。如果你选的网格太粗那些大角度的平面波分量会被混叠到低频区域仿真结果会出现明显的周期性假条纹。我自己的经验是先估算光场中的最大空间频率 (f_{max} \frac{\sin\theta_{max}}{\lambda})其中 (\theta_{max}) 是你要模拟的最大衍射角。对于高斯光束发散角约为 (\theta \approx \frac{\lambda}{\pi w_0})(w_0) 是束腰半径。采样间隔必须满足[ \Delta x \le \frac{1}{2 f_{max}} \frac{\lambda}{2 \sin\theta_{max}} ]否则就会“欠采样”。比如波长 632.8nm、束腰半径 1mm 的光束发散角大约 0.2 mrad那么网格尺寸取几十个微米量级就非常安全。为了留出余量我一般取 (\Delta x \lambda / 20 \sim \lambda / 10)既保证精度又不会让矩阵大到内存爆炸。另外网格总大小要足够覆盖光斑扩散后的范围。传播距离比较长时光斑会扩展如果你的窗口太小边缘的光场会被硬生生截断产生边缘衍射伪影。这里我的建议是让活动网格边长为光束衍射后尺寸的约2~3倍这样既不会浪费计算资源也不会引入明显的边缘效应。3. 用 Python 从零搭一台“数值马赫-曾德”3.1 搭建仿真框架和分束器模型我习惯用 Python NumPy 做仿真偶尔用 SciPy 的信处理工具。原因在于NumPy的数组运算和FFT特别顺手代码量少调试也方便。下面给出大体框架。首先定义一个二维网格import numpy as np wavelength 632.8e-9 w0 0.5e-3 # 束腰半径 N 1024 # 网格尺寸 L 8e-3 # 网格边长 x y np.linspace(-L/2, L/2, N) X, Y np.meshgrid(x, y) delta_x L / N分束器在仿真里不需要真去渲染一块玻璃板只需要把它当作一个空间位置上的算符入射复振幅乘以反射/透射系数。def beam_splitter(field_in, t1, r1, t2, r2): # 返回透射臂和反射臂的光场 return field_in * t1 * r2, field_in * t2 * r1对于理想50:50分束器相位关系特别重要反射和透射光之间通常有 (\pi/2) 的相位差。在仿真中我一般这样处理反射分量乘 (i)也就是相位增加 (\pi/2)透射分量保持不变这样两束光在合束后的干涉条纹位置才是物理正确的。3.2 参考臂与信号臂的光束传播仿真两臂里的光束传播直接调用角谱法函数就可以def propagate_ASM(field, wavelength, delta_x, delta_y, distance): k 2*np.pi/wavelength fx np.fft.fftfreq(N, ddelta_x).reshape(-1,1) fy np.fft.fftfreq(N, ddelta_y).reshape(1,-1) # 空间频率网格 FX, FY np.meshgrid(fx.squeeze(), fy.squeeze()) # 频域传播相位 H np.exp(1j*k*distance*np.sqrt(1 - (wavelength*FX)**2 - (wavelength*FY)**2)) # 注意np.fft.fftfreq返回顺序与fft2一致 F np.fft.fft2(field) out np.fft.ifft2(F * H) return out信号臂里的相位调制可以简单地在光场上乘一个相位屏比如温度变化造成的光程变化我习惯建模为def add_phase(field, phase_map): return field * np.exp(1j * phase_map)例如要在中心区域做一个折射率变化相位分布约等于[ \phi(x,y) \frac{2\pi}{\lambda} \Delta n \cdot L_{cell} ](\Delta n) 是折射率改变量(L_{cell}) 是光与样品相互作用的长度。这个模型在仿真里非常常用比实打实算多层薄膜透过率快得多而且便于参数扫参。3.3 合束、计算干涉条纹与相位提取两臂光场传播到合束位置后合束器同样是一个复振幅叠加的过程。以理想合束器为例输出端光场为E_out1 (E_ref E_sig) / np.sqrt(2) E_out2 (E_ref - E_sig) / np.sqrt(2)然后光强可以直接算I1 np.abs(E_out1)**2 I2 np.abs(E_out2)**2如果你只是要看条纹长什么样画I1的热力图就能获得清晰的干涉图样。但很多时候我们还需要提取相位用得较多的是傅里叶变换法对条纹图做FFT在频谱里提取 ±1 级旁瓣逆变换后再取辐角。def extract_phase(I_all): F np.fft.fft2(I_all) F np.fft.fftshift(F) # 假定条纹空间频率已知, 将中心移到旁瓣 F_shift np.roll(F, (cx_f, cy_f), axis(0,1)) # 只保留旁瓣周围区域 mask np.zeros_like(F_shift, dtypebool) mask[cx_f-size:cx_fsize, cy_f-size:cy_fsize] True F_shift[~mask] 0 phase np.angle(np.fft.ifft2(np.fft.ifftshift(F_shift))) return phase这里有个很关键的操作必须先把F的零频移到中心再用roll把旁瓣搬到零频附近最后逆变换得到的相位才是缓变的相位分布而不是高频载波相位。3.4 仿真参数的选择与调试思路参数的选取直接决定仿真的是否“像真的”。我常用的初始参数如下参数数值说明波长632.8 nmHe-Ne激光典型波长束腰半径0.2~0.5 mm普通氦氖激光直接输出的光斑尺寸两臂长度差0~20 mm决定初始条纹间距和方向分束器分光比50:50保证高对比度干涉网格尺寸1024×1024兼顾精度和速度网格边长4~10 mm要足够覆盖光斑扩展后的范围调试顺序上我建议先跑一个最简单的“零光程差”场景确认两个输出端中有一个是全亮、另一个是全暗——这是MZI合束的最基本特征。如果发现两端光强几乎一样大概率是分束器的相位关系没写对回去检查是不是少乘了 (i)。4. 典型案例微小折射率变化的灵敏度分析4.1 案例设定我用马赫-曾德干涉仪做过一个可追溯复现的仿真案例信号臂里放一个长度 (L_{cell}10) mm 的液体样品池里面液体的折射率相对空气略大当浓度变化导致 (\Delta n 1 \times 10^{-5}) 时我要观察探测器平面上的条纹变化并计算相位变化量。这个场景非常实际比如测血糖浓度、测溶液浓度梯度都是这个思路。4.2 仿真步骤与代码实现核心操作是在信号臂的样品池区域引入一个均匀相位变化[ \Delta\phi \frac{2\pi}{\lambda} \cdot \Delta n \cdot L_{cell} ]代入数值(\Delta\phi \frac{2\pi}{632.8\times10^{-9}} \times 10^{-5} \times 10^{-2} \approx 1.0) rad。也就是大约 0.16 个条纹周期。这么小的相位变化用肉眼直接看条纹图几乎分辨不出来但用相位提取算法可以非常清晰地看到整体相移。整个流程写的伪代码如下# 1. 生成高斯光束 r np.sqrt(X**2 Y**2) E0 np.exp(-(r**2 / w0**2)) # 2. 分束器1 E_ref, E_sig beam_splitter(E0, t11/np.sqrt(2), r11j/np.sqrt(2), t21/np.sqrt(2), r21j/np.sqrt(2)) # 3. 两臂传播 E_ref propagate_ASM(E_ref, wavelength, delta_x, delta_y, L_ref) E_sig propagate_ASM(E_sig, wavelength, delta_x, delta_y, L_sig) # 4. 给信号臂加相位 phase_screen np.ones_like(X) * (2*np.pi/wavelength * delta_n * L_cell) E_sig add_phase(E_sig, phase_screen) # 5. 合束 E_out1 (E_ref E_sig) / np.sqrt(2) E_out2 (E_ref - E_sig) / np.sqrt(2)跑完之后你会发现虽然光强图上的条纹几乎没有可察觉的移动但相位图上亮暗边界的位置整体偏了一个量这个偏移就是灵敏度分析的基础。4.3 结果分析与灵敏度推导灵敏度定义为相位变化量与折射率变化量的比值。对这个模型相位变化近似为[ \frac{d\phi}{dn} \frac{2\pi L_{cell}}{\lambda} ]代入 (L_{cell}10) mm、(\lambda632.8) nm灵敏度约为 (9.9 \times 10^4) rad/RIURIU 折射率单位。什么意思呢折射率每变化 1相位走将近十万弧度如果探测器能分辨 0.01 rad 的相位理论上能探测的折射率变化就是 (1 \times 10^{-7}) 量级。这也是为什么马赫-曾德干涉仪可以作为高灵敏度传感器的原因。不过要注意这只是理想情况。实际中光源功率波动、探测器噪声、温度漂移都会拉低真实探测能力。仿真能帮你确认“理想情况下系统能做到什么水平”进而判断瓶颈到底在光学部分还是电路部分。5. 仿真结果的验证与精度控制5.1 用解析解检验仿真精度仿真写得再漂亮也得有办法验证对不对。最直接的验证方式是对照理想双光束干涉解析解。对于完全均匀的平面波理想MZI输出光强应该是[ I I_0 \left[1 \cos(\phi_1 - \phi_2)\right] ]我可以让仿真里两臂均为均匀平面波也就是忽略高斯振幅记录不同相位差下的中心点光强和解析式画在同一张图里看是否重合。这个测试只需要几分钟就能完成但能迅速暴露分束器相位写错、网格尺寸过小等低级错误。此外还有一个检验方法把两臂长度设成相同再做一组小角度倾斜第二分束器的仿真。理论预测条纹间距应正比于两束光夹角的正弦倒数。如果仿真出来的条纹间距符合经验公式就说明角谱法的传播模型工作正常。5.2 影响精度的三个“隐形敌人”实操中影响仿真精度的因素不少我踩过的坑主要集中在下面三个第一个是空间混叠。网格尺寸太大会让高频衍射分量混叠成低频表现为干涉条纹周围出现歪歪扭扭的额外纹路。解决办法就是满足采样条件最直接的检测方法把网格尺寸减半跑一次看结果是否明显变化。如果不变化说明此前设置已经收敛如果变了那就继续加密网格。第二个是边缘衍射。计算窗口面积有限光场边缘被突然截断会产生边缘衍射波。这种效应在光束尺寸和窗口尺寸相当的时候非常严重看起来像是条纹被扰动。解决思路是扩大窗口或者给入射光场加一个超高斯孔径函数让光场在边缘平滑过渡到零。第三个是相位卷绕。用反正切函数提相位时相位限制在 ((-\pi, \pi]) 范围内。如果真实相位变化超过 (2\pi)提取结果就会出现台阶状跳变也就是所谓“包裹相位”。用np.unwrap可以解开一维相位但二维相位展开需要更细致的路径算法。最好的办法是控制相位变化在每相邻像素之间小于 (\pi)也就是保证相位在空间上有足够的过采样。5.3 提高精度的实用手段如果你发现数值结果始终有微弱的周期误差优先做下面三件事网格零填充在光场周围补一圈零相当于扩大计算窗口能有效降低边缘衍射高斯窗滤波在提取相位前对光强图乘一个高斯窗压制旁瓣泄漏草坪式参数扫描固定两臂光程差微小改变网格尺寸做收敛性验证找到“结果不再变化”的最低网格需求。我通常会写一个简短参数扫描脚本自动改变 (N) 从256到2048记录中心相位结果当数值变化小于千分之一时就用这组参数。这个习惯在我后来做多组仿真时省了大量时间因为它把“靠经验调参”变成“靠数据说话”。6. 常见问题与调试技巧实录6.1 仿真条纹模糊或不清晰如果你明明设置了两臂的相位差却发现干涉条纹对比度非常低先检查两个分束器输出的两束光在合束点的振幅是否接近。对比度 (V) 的公式很简单[ V \frac{I_{max}-I_{min}}{I_{max}I_{min}} \frac{2\sqrt{I_1 I_2}}{I_1I_2} ]当两束光的振幅差别大时条纹就被“稀释”了。常见原因光束经过长距离传播后由于衍射振幅分布发生改变或者分束器的分光比设成了 90:10。如果你要的是高对比度就得让两臂光强接近相等。这个经验在实验中同样适用。6.2 相位提取结果跳变用傅里叶变换法提取相位时如果载频选择不准提取的相位会带有周期性条纹残留。解决办法是把条纹频率选在频谱旁瓣的中心而不是随便用理论计算值。实际操作中我会先对干涉图做FFT找到幅度最大的非零频峰再把旁瓣中心对准那个峰值。这个步骤看似麻烦但比手动输入频率值要可靠得多。6.3 不同波长和光束尺寸下仿真的变化波长变化会同时影响衍射传播相位和干涉条纹灵敏度这是新手最容易漏掉的地方。仿真时如果从 632.8 nm 换成 1550 nm不仅条纹间距会变衍射效果也会增强因为衍射角正比于 (\lambda/d)。所以一套参数不能通吃所有场景至少要做一次“波长-网格采样”的自检。光束尺寸也一样。束腰越小发散越快同样的传播距离下光斑扩展越大。如果你的网格只按初始光斑尺寸算跑到一半光就顶到边界了。我的习惯是把网格边长设在光斑最大尺寸的2~3倍宁可多算一点像素也不要边界反射干扰真实信号。6.4 个人心得从仿真到实验的几点建议仿真跑通之后最好顺手把结果和实际实验对比一轮。我每次搭MZI之前都会先做一个“预演仿真”用实验室真实的光束参数跑一遍看看探测器各个位置的光强分布大概是什么量级提前判断该用多大增益、会不会饱和。这套流程看起来不起眼但真的帮我避免过多次烧光电探测器或者白调一晚上光路的尴尬。另外提醒一句仿真里的“分束器相位关系”是很多人容易搞错的地方。物理上反射和透射存在相位差但不同教材的约定可能不同。我的建议是不管约定是反射加 (\pi/2) 还是透射加 (\pi)只要你在仿真里保持自洽并且保证合束后两路光的相对相位关系正确最后的物理结论就不会出错。可如果符号不一致你会在条纹亮暗位置和相位正负上反复栽跟头。最后再分享一个实用小技巧在做参数扫描时把每次的两臂相位差记录下来画成“相位-响应”曲线往往能直接看出干涉仪的线性范围和周期特性。这比盯着单幅条纹图反复看要高效得多而且写报告时这张图非常能说明问题。接下来如果再往下扩展你还可以加入透镜聚焦、光纤耦合、非理想分束器波长依赖等模块让仿真一步步逼近真实系统。