从大学光学课第一次看到多光束干涉的公式开始我就一直觉得这东西特别奇妙。明明都是同一束光分出来的为什么双光束干涉条纹是正弦变化而多光束干涉能出现又细又亮的锐利条纹当时公式推导能看懂但总觉得少了点直觉。后来工作里做光学仿真才意识到靠手算是永远不可能“看见”多光束干涉的真正面貌的必须借助数值工具。于是我用Matlab搭了一套多光束干涉仿真程序今天把整个思路和踩过的坑都写出来。这篇文章适合正在学光学、做光电设计或者对Matlab光学仿真感兴趣的读者核心是从原理到代码一步步复现多光束干涉并学会用参数扫描看清物理本质。我在实际仿真的过程中最深的感受是多光束干涉和双光束干涉最大的区别不在于“多个光束”这个字面意思而在于能量重新分布的方式发生了质变。只有把仿真做出来、把图样渲染出来、把参数扫起来你才能真正理解那些教科书上写的“条纹锐利度随反射率增大而急剧提升”到底意味着什么。1. 内容整体设计与思路拆解1.1 多光束干涉从物理本质到仿真价值多光束干涉顾名思义是两束以上的相干光在空间某一点叠加后因为相位差形成明暗交替的强度分布。最常见的物理场景是法布里-珀罗干涉仪和衍射光栅一束光在两面高反射镜之间来回反射每一次透射出去的光都和前面透射出去的光相干叠加这就是典型的多光束干涉。衍射光栅则是成千上万个狭缝的次级波相干叠加也是多光束干涉的体现。我在做仿真前花了很长时间想一个问题为什么要在Matlab里做这件事直接用公式画图不行吗答案是教科书只给你最终表达式但仿真能让你看到每一项如何贡献到最终结果。比如多光束干涉的强度公式里有反射率R这个参数理论上你知道R越大条纹越尖锐但你不亲手扫一遍R从0.04到0.95的图样变化很难建立起“条纹锐度”对这个参数的敏感度到底有多高。这种通过参数扫描获得物理直觉的过程就是仿真的核心价值。另外Matlab做这类仿真还有一个天然优势矩阵运算效率高代码表达和数学公式几乎一一对应排错和修改都很方便。相比用C或者PythonMatlab在交互式探索光学图样这件事上确实更顺手。1.2 方案选型为什么用解析建模而非FDTD光学仿真领域有多种工具例如FDTD Solutions、COMSOL等它们做的是电磁场数值求解把空间离散成网格逐时间步推进。这类方法精度高能处理复杂结构但计算量大、学习成本高而且对多光束干涉这种物理图像非常清晰的问题来说属于“杀鸡用牛刀”。我的思路是采用解析建模直接用光的叠加原理把多束光在观察面上的复振幅叠加起来再模方得到强度分布。这样做的好处有几点。计算速度快参数改变后几乎实时更新结果。代码逻辑清晰每一步都能和物理公式对应。便于深挖参数影响因为每个变量都是显式的不存在数值噪声干扰判断。解析建模也有局限比如无法处理偏振态在界面的复杂变化、近场效应等。但对多光束干涉的教学演示和工程预研来说这个精度完全够用。1.3 Matlab光学工具箱的定位与本文的路径搜Matlab的光学工具箱时我发现其实Matlab没有专门的多光束干涉工具箱它有的更多是图像处理、信号处理方面的能力。有些人在做类似仿真时会把光学问题抽象成图像卷积另一些人用Simulink做光学通信系统仿真。但我们要做的是基础物理过程的建模所以最直接的路径是用Matlab的基础矩阵运算和绘图函数自己写。你也可以尝试用Symbolic Math Toolbox做符号推导或者用App Designer做一个交互界面但我建议第一步还是把核心物理模型写清楚、画出来。工具是辅助物理思路才是主线。2. 核心细节解析与实操要点2.1 光波叠加的数学模型从双光束到N光束光波是电磁波在空间某一点的电场可以用复振幅表示。假设有N束光在观察面上叠加每束光的振幅为E₀先假设一致相位为φᵢ则合成复振幅为E_total Σ E₀ · exp(i·φᵢ)光强I与|E_total|²成正比。双光束时这个求和能得到简单的余弦表达式。多光束时引入等比数列求和公式可以推导出解析表达式。但在Matlab里我选择直接用复数数组模拟求和过程代码写起来更直观也更容易扩展到任意N值。我需要小心相位差的来源。在多光束干涉里相位差来源于光程差。对于等间距排列的N个点光源或N个狭缝相邻光束的光程差为 Δ d·sinθ对应的相位差为 δ (2π/λ)·d·sinθ。于是第i个光束的相位为 φᵢ i·δ加上初始相位偏移。这个模型是光栅方程和多光束干涉统一的基础。我建议编程时把波长λ、间距d、光束数N、初始相位四个参数单独设置方便后续扫描。2.2 干涉图样的两种观察方式远场角度分布与空间分布仿真中有一个很容易混淆的地方干涉图样到底是“空间位置”的函数还是“角度”的函数远场观察方式比如光栅的夫琅禾费衍射横轴是sinθθ是观察方向与法线的夹角。强度分布是角度的函数。近场观察方式比如双缝后方某个平面横轴是空间坐标x强度分布是位置的函数。两种方式公式形式非常相似只要把x/D近似为sinθD是观察屏到光源的距离就可以互相转换。我在代码里统一用角度变量theta在区间[-π/2, π/2]上扫描既方便画一维曲线也方便生成二维图样时建立空间网格。2.3 参数选择对仿真结果的决定性影响做仿真不是随便填几个参数就完事。我总结出几个关键经验和大家分享。波长λ的选择可见光380~780nm。如果在仿真里把λ设成1归一化那么相位差δ 2π·d·sinθ就变成纯几何关系。归一化后图形更好看但物理直觉会减弱。我的习惯是保留SI单位制比如λ632.8nm氦氖激光因为这样出来坐标轴带单位后续如果要对接工程参数更方便。光束数NN2时是双光束干涉N5时条纹开始明显变锐N100时接近光栅的行为。建议从2开始逐步增加观察干涉图样是如何从正弦条纹演化为锐利主极大加次级极大。间距d与波长λ的比值这个比值决定了干涉极大的角度位置。如果d/λ太小只有零级附近有几个极大角度范围很窄如果d/λ很大角度方向上会出现非常多的极大图样会密集到难以分辨。刚开始仿真时建议取d/λ2~5。2.4 代码实现的三个模块划分为了保持代码清晰我把整个仿真拆成三个模块参数定义模块设置波长、间距、光束数、振幅、观察角度范围。核心计算模块计算各光束相位叠加复振幅得到强度。可视化模块绘制一维强度曲线、二维干涉图样、极坐标图等。这种划分方便后续扩展。比如你想改成研究随机相位扰动只需要在核心计算模块里加入随机数即可你想改成研究不同入射角只需要在相位表达式中加入入射角项。3. 实操过程与核心环节实现3.1 环境准备Matlab版本与依赖我使用的是Matlab R2021a实际上从R2016b开始这段代码都能直接跑不需要额外的工具箱。基础的矩阵运算和plot、imagesc、surf等绘图函数都属于Matlab核心功能。如果你用的是更老的版本只需要注意Implicit Expansion特性在R2016b引入是否可用否则要用bsxfun来扩展数组维度。有一点要提醒很多人在网上下载的Matlab安装包可能版本较旧如果遇到函数兼容性问题优先检查你的数组维度操作是否适合当前版本。比如下述代码中theta和delta数组的形状要匹配老版本可能需要加repmat处理。3.2 一维多光束干涉强度分布先看最核心的计算代码。我定义了一个函数输入是波长lambda、间距d、光束数N、观察角度theta数组输出是归一化强度I。这段代码是我后来重构过的版本最初一版用了两个嵌套循环虽然也能跑但速度慢很多改成向量化之后快了十倍不止。function I multi_beam_interference(lambda, d, N, theta) % 多光束干涉强度分布 % lambda: 波长 (m) % d: 相邻光束间距 (m) % N: 光束数目 % theta: 观察角度数组 (rad) k 2 * pi / lambda; % 波数 delta k * d * sin(theta); % 相邻光束相位差数组维度与theta相同 % 构建N行length(theta)列的相位矩阵 % 第i行第j列表示第i束光在角度theta(j)处的相位 i_idx (0:N-1); % 列向量 phase_matrix i_idx * delta; % 隐式扩展得到 N x M 矩阵 % 复振幅叠加 E sum(exp(1j * phase_matrix), 1); % 沿第一维求和得到1 x M行向量 % 强度并归一化 I abs(E).^2 / N^2; % 除以N^2使最大强度为1 end这段代码中最关键的一行是phase_matrix i_idx * delta。如果要兼容老版本Matlab可以改成phase_matrix repmat(i_idx, 1, length(delta)) .* repmat(delta, N, 1)。这里展开的原因是我需要让每一束光在每一个角度下都计算一次相位本质上是一个二维矩阵运算而向量化编码正好契合Matlab的设计哲学。调用方式很简单。比如我想算λ632.8nm、d2μm、N12条光束的情况角度范围从-30度到30度lambda 632.8e-9; d 2e-6; N 12; theta linspace(-pi/6, pi/6, 2000); I multi_beam_interference(lambda, d, N, theta); figure; plot(theta * 180/pi, I, b-, LineWidth, 1.5); xlabel(观察角度 (度)); ylabel(归一化强度); title(多光束干涉一维强度分布); grid on;这里角度采样点数2000是一个经过权衡的值。点太少峰值位置和宽度会失真点太多比如50000Matlab绘图响应会明显变慢。如果后续要做参数扫描循环建议进一步降低到1000点对图形趋势的影响几乎看不出来。3.3 二维干涉图样把强度映射成平面图案一维图只能看某一条线上的强度二维图才能直观展示干涉图样的空间分布。实际做法是把二维观察屏上的每个像素看成一个观察方向用坐标换算得到该像素对应的角度。假设观察屏在距离光源L远处像素坐标为(x, y)则角度约为sinθ_x ≈ x / sqrt(x² y² L²)sinθ_y ≈ y / sqrt(x² y² L²)如果只关心小角度区域还可以直接近似为sinθ_x ≈ x/L。我写了一个生成二维图样的脚本lambda 632.8e-9; d 2e-6; N 10; L 1; % 观察屏距离1m pixels 500; % 每边像素数 range 0.02; % 观察屏范围米 x linspace(-range, range, pixels); y linspace(-range, range, pixels); [X, Y] meshgrid(x, y); theta_x atan(X / L); theta_y atan(Y / L); % 二维情况下相位差是x和y方向相位差的矢量和 kx 2 * pi / lambda; ky 2 * pi / lambda; delta_x kx * d * sin(theta_x); delta_y ky * d * sin(theta_y); delta delta_x delta_y; % 严格来说要看你光栅刻线方向这里假设x方向刻线 i_idx (0:N-1); phase_matrix i_idx .* reshape(delta, 1, pixels, pixels); E sum(exp(1j * phase_matrix), 1); I squeeze(abs(E).^2 / N^2); imagesc(x * 1e3, y * 1e3, I); axis image; colormap(hot); colorbar; xlabel(x (mm)); ylabel(y (mm)); title(sprintf(%d光束干涉二维图样 (lambda%.1fnm, d%.1fum), N, lambda*1e9, d*1e6));这段代码需要注意squeeze的使用。三维数组经过sum后中间多了一个长度为1的维度必须压掉才能用imagesc。我第一次写的时候忘了加squeeze结果imagesc把三维数组当成RGB数据画出来全图一团黑排查了半天才找到原因。3.4 极坐标可视化MATLAB polarplot的妙用有时候笛卡尔坐标下的强度曲线不够直观尤其是在描述多光束干涉的方向性时极坐标图能更清楚地表达“哪些方向有光、哪些方向没光”。Matlab的polarplot函数可以帮助我们画极坐标下的强度分布。theta linspace(-pi/2, pi/2, 2000); I multi_beam_interference(632.8e-9, 2e-6, 8, theta); figure; polarplot(theta, I, b-, LineWidth, 1.5); title(多光束干涉极坐标强度分布);需要注意polarplot对数据范围比较敏感。如果theta从-pi/2到pi/2Matlab默认的极坐标图会在角度方向直接映射可能导致图形看起来只占了半圈。我习惯把theta扩展到-pi到pi同时在另一半补零让极坐标图看起来更完整对称。另外polarplot中坐标轴字体属性调整和普通plot不太一样它要用rticks、thetaticks这类专属命令。3.5 参数扫描与动画让物理“动”起来静态图看多了参数扫描才能带来突破性认知。我写过一个循环扫描光束数N从2到50每个N绘制一维强度曲线并保存为帧最后合成动画。你会发现一个非常震撼的过程N2时条纹是宽宽的余弦峰N5时峰开始变窄N20时主极大周围出现了明显的次级峰N50时主极大锐利得像一根针。核心循环代码如下lambda 632.8e-9; d 2e-6; theta linspace(-0.3, 0.3, 3000); N_list [2 3 4 5 8 10 15 20 30 50]; figure(Position, [100 100 800 500]); for idx 1:length(N_list) N N_list(idx); I multi_beam_interference(lambda, d, N, theta); plot(theta * 180/pi, I, b-, LineWidth, 1.5); xlabel(角度度); ylabel(归一化强度); title(sprintf(N %d 光束干涉, N)); ylim([0 1]); grid on; drawnow; frame getframe(gcf); [A, map] rgb2ind(frame.cdata, 256); if idx 1 imwrite(A, map, N_scan.gif, gif, LoopCount, Inf, DelayTime, 0.6); else imwrite(A, map, N_scan.gif, gif, WriteMode, append, DelayTime, 0.6); end end这个动画我发给过不少学生反馈都说“看了动画才真正理解多光束干涉和双光束干涉的区别”。我也建议你自己跑一遍观察主极大半宽的变化规律。4. 常见问题与排查技巧实录4.1 主极大位置偏移相位计算中的符号问题我做仿真时遇到的第一个诡异问题是主极大位置不对。理论上N个等间距光束的主极大应该出现在满足d·sinθ mλ的角度上但实测仿真得到的极大角总是往负方向偏移一点。排查半天发现相位矩阵构建时我把方向搞反了。第0个光束的相位如果是0第i个光束的相位如果是i·δ那么叠加结果对应的是正向传播。但如果我不小心用了负号phase_matrix i_idx * delta写成i_idx * (-delta)整个图样就会镜像翻转。这个问题在参数对称时会看不出来但只要把入射角设成非对称的问题立刻暴露。4.2 条纹过密无法分辨角度范围与像素数的平衡当d/λ较大比如10的时候sinθ每变化0.1就会有多个主极大出现如果角度范围设得太大比如-60度到60度整个图像会密密麻麻全是条纹根本看不出结构。解决办法有两个一是缩小观察角度范围聚焦在零级附近二是增加采样点。我实践中发现对于d/λ5的情况0.6弧度的角度范围内至少需要3000个采样点否则峰形会严重变形。如果你只是为了看趋势2000个点够用如果要做精确的半高宽分析建议5000点以上。4.3 内存溢出与运行过慢向量化与并行化的取舍我第一次写二维干涉图样时用了三层嵌套循环遍历每个像素、每束光结果500×500分辨率的图案跑了将近3分钟内存占用也高得吓人。换成向量化后同样的结果0.2秒就出来了。Matlab里向量化永远优先于显式循环。如果你的参数扫描需要同时遍历N、d、λ、相位等多个维度建议先写好单次计算的核心函数然后用parfor并行起来。我试过在一台6核机器上把4个参数的扫描任务并行化速度提升接近5倍。但要提醒并行循环里的绘图需要特别小心不能直接在每个worker里画图而要收集结果后统一绘图。4.4 归一化陷阱强度最大值未必在0级多光束干涉强度公式的归一化有一个陷阱理论上N束光同相位时强度可以达到N²E₀²归一化后是1。但如果观察角度范围不包含主极大对应的角度max(I)就小于1导致后续分析出现偏差。我的处理方式是在绘图前先检查max(I)是否接近1如果不是说明角度范围或参数设置有问题。4.5 常见问题速查表现象可能原因排查方法主极大位置偏移相位符号反了检查phase_matrix中i_idx和delta的乘法顺序与正负号图样过于密集角度范围太大或d/λ太大缩小theta范围或减少d/λ强度最大只有0.5观察角度不含主极大扩展theta范围二维图像全黑未用squeeze压缩维度对sum的结果加squeeze运行速度极慢使用了显式循环改成矩阵向量化运算条纹有锯齿采样点不足增加采样点数目极坐标图只显示半圈theta范围过窄扩展到-pi到pi对称区间5. 实战扩展从基础仿真到工程应用5.1 高斯光束与多光束干涉的结合实际激光工程里极少有理想平面波做多光束干涉的。激光器输出的是高斯光束振幅不是均匀的而是随位置呈高斯分布。如果你想仿真这个效果只需要在核心计算模块中为每束光乘上高斯包络因子。这样得到的干涉图样会出现很明显的“中心亮、边缘暗”的调制和实验照片更接近。具体实现是在相位矩阵计算完成后新建一个振幅矩阵每一列乘上对应的角向高斯权重然后叠加。你会发现高斯包络的作用是抑制旁瓣让图样看起来更干净。这在光学相控阵的设计中非常有用。5.2 随机相位扰动的影响仿真真实光路中不可能做到完全相干。环境的震动、温度引起的光程抖动都会给每束光引入随机相位。把这个扰动加入仿真只需要在相位矩阵上叠加一个随机矩阵。我做过蒙特卡洛模拟重复500次后统计平均强度可以看到随机相位会显著降低条纹对比度。这个结果对评估光学系统的稳定性很有参考价值。% 在核心计算中加入随机相位扰动 phase_noise sigma_noise * randn(N, length(theta)); phase_matrix_with_noise phase_matrix phase_noise; E sum(exp(1j * phase_matrix_with_noise), 1); I_noisy abs(E).^2 / N^2;sigma_noise从0逐渐增大到2π你会看到干涉条纹从清晰到完全消失的完整过程。光学里管这叫“退相干”。这个仿真花不了几分钟但对理解相干性的物理意义帮助极大。5.3 与MATLAB图像处理模块的联动分析当你把二维干涉图样生成后可以调用图像处理工具箱做进一步分析。比如用imregionalmax找出亮斑质心用bwdist做条纹间距测量或者用fft2做空间频率分析。这一步可以把“物理仿真”和“图像诊断”打通非常接近实际工程中光学测量的流程。我甚至试过把干涉图样保存成图片再用图像处理流程自动统计条纹间距反过来推算光源波长。仿真精度高的时候反推结果和真实波长误差能控制在0.1%以内。这说明无论仿真还是实际实验物理图像的一致性都是可靠的。5.4 适合进阶的扩展方向进阶可以考虑的问题还有不少。比如把一维光栅改成二维光栅阵列观察点阵状干涉图样把等间距改成啁啾间距看看聚焦效应把相位调制做成随机编码模拟波前整形甚至可以把核心函数打包成App Designer应用做一个交互式光学教学小工具。我在工程中实际用到的是把多光束干涉仿真拓展到相控阵天线的方向图计算上。微波和光学的数学基础完全一致把波长换成天线工作波长把光束间距换成阵元间距干涉公式就直接变成了阵列天线的方向图公式。这也是很多电扫阵列建模的底层原理。我当时拿到一个8×8阵列的波束扫描需求第一反应就是套用这套多光束干涉的代码只不过把一维扩展成二维把光频换成微波频率几个小时内就给出了初步的方向图分析效率比从零开始写快了太多。6. 实操心得我踩过的坑和想对你说的6.1 不要忽视基础数学推导做仿真最容易犯的错误是拿到公式就写代码跳过物理推导。我建议不管代码多简单先用手推一遍N2和N3的情况明确每一项的物理意义。有了这一步代码里的矩阵维度、相位符号、归一化系数才不容易出错。6.2 参数归一化是一个双刃剑很多教材喜欢把所有长度量归一化到波长省略单位结果图形确实简洁了但物理直觉容易丢。我个人的习惯是调试阶段全部用SI单位让每个中间变量都有物理含义等到生成论文插图或者做教学演示时再归一化处理坐标轴。两种模式切换成本很低但收益是双向的。6.3 绘图技巧让结果自己“讲故事”Matlab的绘图功能很强但默认配色和样式确实比较朴素。做多光束干涉仿真时我推荐把颜色图设成hot或者turbo这样干涉亮斑的层次感明显得多。曲线图默认的蓝色实线也建议加粗到1.5以上否则发表到文档里看不太清。还有一点所有图的坐标轴字体大小尽量统一比如设为12号或14号这会让整套仿真结果看起来非常专业。6.4 从仿真到理解的最后一公里必须承认仿真的意义不在于“把图画出来”而在于“从图里看出物理”。我建议你在跑完参数扫描后尝试用自己的话解释以下现象。为什么N增大时主极大变窄为什么主极大和次级极大之间有N-2个暗纹为什么增加反射率等效于增加光束数如果你能流利回答这些问题说明这个仿真真正起到了作用。否则你只是按了运行按钮并没有把它变成自己的知识。我一直想做一个交互式的多光束干涉教学演示面板把光束数、波长、间距、相位差都变成滑块让操作者在屏幕上拖一拖就能看见干涉图样变化。这个想法后来用Matlab的App Designer实现了过程不算复杂但效果比我预想中好很多很多朋友反馈说“玩着玩着就理解了”。如果你也想做我建议从本文的核心函数出发先做一个单参数滑块的版本再逐步叠加。