多光束干涉Matlab仿真:从原理到参数扫描的完整实践
2026/9/10 0:01:06 网站建设 项目流程

从大学光学课第一次看到多光束干涉的公式开始,我就一直觉得这东西特别奇妙。明明都是同一束光分出来的,为什么双光束干涉条纹是正弦变化,而多光束干涉能出现又细又亮的锐利条纹?当时公式推导能看懂,但总觉得少了点直觉。后来工作里做光学仿真,才意识到靠手算是永远不可能“看见”多光束干涉的真正面貌的,必须借助数值工具。于是我用Matlab搭了一套多光束干涉仿真程序,今天把整个思路和踩过的坑都写出来。这篇文章适合正在学光学、做光电设计或者对Matlab光学仿真感兴趣的读者,核心是从原理到代码一步步复现多光束干涉,并学会用参数扫描看清物理本质。

我在实际仿真的过程中,最深的感受是:多光束干涉和双光束干涉最大的区别,不在于“多个光束”这个字面意思,而在于能量重新分布的方式发生了质变。只有把仿真做出来、把图样渲染出来、把参数扫起来,你才能真正理解那些教科书上写的“条纹锐利度随反射率增大而急剧提升”到底意味着什么。

1. 内容整体设计与思路拆解

1.1 多光束干涉:从物理本质到仿真价值

多光束干涉,顾名思义是两束以上的相干光在空间某一点叠加后,因为相位差形成明暗交替的强度分布。最常见的物理场景是法布里-珀罗干涉仪和衍射光栅:一束光在两面高反射镜之间来回反射,每一次透射出去的光都和前面透射出去的光相干叠加,这就是典型的多光束干涉。衍射光栅则是成千上万个狭缝的次级波相干叠加,也是多光束干涉的体现。

我在做仿真前花了很长时间想一个问题:为什么要在Matlab里做这件事?直接用公式画图不行吗?答案是:教科书只给你最终表达式,但仿真能让你看到每一项如何贡献到最终结果。比如多光束干涉的强度公式里有反射率R这个参数,理论上你知道R越大条纹越尖锐,但你不亲手扫一遍R从0.04到0.95的图样变化,很难建立起“条纹锐度”对这个参数的敏感度到底有多高。这种通过参数扫描获得物理直觉的过程,就是仿真的核心价值。

另外,Matlab做这类仿真还有一个天然优势:矩阵运算效率高,代码表达和数学公式几乎一一对应,排错和修改都很方便。相比用C++或者Python,Matlab在交互式探索光学图样这件事上确实更顺手。

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(氦氖激光),因为这样出来坐标轴带单位,后续如果要对接工程参数更方便。
  • 光束数N:N=2时是双光束干涉,N=5时条纹开始明显变锐,N=100时接近光栅的行为。建议从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、d=2μm、N=12条光束的情况,角度范围从-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是一个经过权衡的值。点太少,峰值位置和宽度会失真;点太多,比如50000,Matlab绘图响应会明显变慢。如果后续要做参数扫描循环,建议进一步降低到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/2,Matlab默认的极坐标图会在角度方向直接映射,可能导致图形看起来只占了半圈。我习惯把theta扩展到-pi到pi,同时在另一半补零,让极坐标图看起来更完整对称。另外polarplot中坐标轴字体属性调整和普通plot不太一样,它要用rticks、thetaticks这类专属命令。

3.5 参数扫描与动画:让物理“动”起来

静态图看多了,参数扫描才能带来突破性认知。我写过一个循环,扫描光束数N从2到50,每个N绘制一维强度曲线并保存为帧,最后合成动画。你会发现一个非常震撼的过程:N=2时条纹是宽宽的余弦峰,N=5时峰开始变窄,N=20时主极大周围出现了明显的次级峰,N=50时主极大锐利得像一根针。

核心循环代码如下:

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 不要忽视基础数学推导

做仿真最容易犯的错误是拿到公式就写代码,跳过物理推导。我建议不管代码多简单,先用手推一遍N=2和N=3的情况,明确每一项的物理意义。有了这一步,代码里的矩阵维度、相位符号、归一化系数才不容易出错。

6.2 参数归一化是一个双刃剑

很多教材喜欢把所有长度量归一化到波长,省略单位,结果图形确实简洁了,但物理直觉容易丢。我个人的习惯是:调试阶段全部用SI单位,让每个中间变量都有物理含义;等到生成论文插图或者做教学演示时,再归一化处理坐标轴。两种模式切换成本很低,但收益是双向的。

6.3 绘图技巧:让结果自己“讲故事”

Matlab的绘图功能很强,但默认配色和样式确实比较朴素。做多光束干涉仿真时,我推荐把颜色图设成hot或者turbo,这样干涉亮斑的层次感明显得多。曲线图默认的蓝色实线也建议加粗到1.5以上,否则发表到文档里看不太清。还有一点,所有图的坐标轴字体大小尽量统一,比如设为12号或14号,这会让整套仿真结果看起来非常专业。

6.4 从仿真到理解的最后一公里

必须承认,仿真的意义不在于“把图画出来”,而在于“从图里看出物理”。我建议你在跑完参数扫描后,尝试用自己的话解释以下现象。

  • 为什么N增大时主极大变窄?
  • 为什么主极大和次级极大之间有N-2个暗纹?
  • 为什么增加反射率等效于增加光束数?

如果你能流利回答这些问题,说明这个仿真真正起到了作用。否则,你只是按了运行按钮,并没有把它变成自己的知识。

我一直想做一个交互式的多光束干涉教学演示面板,把光束数、波长、间距、相位差都变成滑块,让操作者在屏幕上拖一拖就能看见干涉图样变化。这个想法后来用Matlab的App Designer实现了,过程不算复杂,但效果比我预想中好很多,很多朋友反馈说“玩着玩着就理解了”。如果你也想做,我建议从本文的核心函数出发,先做一个单参数滑块的版本,再逐步叠加。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询