我们做光学的迟早会碰上多光束干涉。不管是法布里-珀罗腔、衍射光栅还是薄膜滤光片背后都是同一套物理多列相干光叠加形成尖锐或平缓的干涉条纹。我当年读书时第一次用Matlab跑出多光束干涉图样看到主极大旁边整齐排列的次级峰突然觉得以前死活记不住的光栅方程变得特别直观。这篇文章就把这个仿真项目完整拆开讲从物理模型到代码实现再到调试经验一次说透。1. 项目概述与整体设计思路1.1 这个仿真项目到底在做什么多光束干涉简单说就是N列频率相同、相位差固定的光波在空间某一点叠加。常见场景有两类一类是N个等间距的相干光源比如衍射光栅的N条缝另一类是光在平行板之间来回反射多次后透射或反射比如法布里-珀罗干涉仪。用Matlab做这个仿真核心就一句话把每束光的复振幅算出来逐点叠加再取模平方得到光强分布。听起来简单但真做起来有不少细节——相位差怎么算、边界怎么处理、N取多少合适、角度范围怎么选都直接影响结果好不好看、物理上对不对。我选择Matlab而不是Python或COMSOL原因很实在Matlab的数组运算天然适合这种逐点扫描的物理仿真一行exp(1i * phi)就能生成整个角度范围内的复振幅序列画图交互也方便改参数重跑一遍非常顺手。如果你手头只有Python思路完全一样把语法翻译过去就行。1.2 模型选择为什么从光栅衍射入手多光束干涉的入门模型我推荐用**N缝衍射光栅**而非F-P腔。原因有两个光栅模型的物理图像最清晰N列光从不同位置出发在远场某一点汇聚几何关系直接相位差好算。它和F-P腔的数学结构完全同构都是叠加N个等幅、相位等差递增的复振幅光栅给出的是透射方向分布F-P腔给出的是频率/相位分布。等你把光栅模型跑通了再把相位差公式换成F-P腔的往返相位代码改动不超过五行。提示如果你的目标是用Matlab做电扫阵列建模与仿真这套多光束干涉代码其实就是相控阵天线方向图仿真的光学版。阵列天线的阵因子、栅瓣条件、主瓣宽度与光栅的干涉项、主极大、谱线半宽在数学上完全一一对应。跑通这个光学仿真等于把阵列信号处理的物理直觉也顺带建立了。1.3 技术路线总览整个项目的技术链路如下建立坐标系与参数波长、缝间距、缝数、观察角范围。计算相邻光束的相位差delta k * d * sin(theta)。用等比数列求和公式计算总复振幅得到解析解。用逐束累加法验证解析解排除公式写错的可能。绘制归一化光强分布观察主极大、次级大、缺级现象。若继续做F-P腔替换相位差模型改变参数R观察精细度变化。这套流程从简到繁每步都有明确的物理意义和验证手段跑数不心虚。2. 多光束干涉的核心物理与数学基础2.1 相位的来源光程差决定一切多光束干涉的灵魂是相位差。N束光在观察点P叠加时各束光的实际传播路径长度不同导致到达P点时有不同的相位偏移。对光栅来说相邻缝到P点的光程差是 d·sinθ因此相邻光束的相位差是delta k * d * sinθ (2π / λ) * d * sinθ这个公式是整套仿真的地基。你要注意几个细节θ是观察方向与光栅法线的夹角不是与光栅面的夹角。d是相邻缝中心的间距不是缝宽。只考虑远场夫琅禾费衍射所以光程差近似用平行光模型这里的近似在仿真里其实是精确的因为我们直接把观察角θ作为自变量扫描。理解了这个相位差后面的叠加就水到渠成。2.2 等比数列求和多光束叠加的优雅解法N束等幅、相位依次相差delta的光总复振幅可以写成U 1 exp(i·delta) exp(i·2·delta) ... exp(i·(N-1)·delta)这是个标准的等比数列公比 q exp(i·delta)求和后取模平方得到I(θ) I0 * [sin(N·delta/2) / sin(delta/2)]^2这个公式在Matlab里一行代码就能实现但有几个坑sin(delta/2) 0 时分子分母同时为0需要特殊处理在编程时这会导致0/0的NaN。我见过很多人直接跑出满屏NaN后开始怀疑人生。解决办法很简单在分母上加一个极小量或者独立检测分母为零的点并单独赋值。还有一个细节总能量守恒的归一化。N束光叠加主极大峰值强度正比于N²所以绘图前要把光强除以N²否则N增大时曲线纵坐标会飞出坐标轴。2.3 从N2到N6看清干涉条纹怎么变尖锐如果你只仿真N2双缝干涉会得到经典的余弦条纹宽、平、对比度100%。但你加到N6时会看到主极大明显变锐主峰两侧出现N-24个次级小峰。这就是多光束干涉最迷人的地方——叠加的光束越多干涉条纹越尖锐。继续增大到N100主极大变得像刀锋一样次级峰虽然还在但已经很弱。这正是衍射光栅能精确测量波长的原因条纹越锐分辨能力越强。在Matlab里观察这个变化非常直观设一个for循环让N从2递增到20然后生成一个动态图肉眼可见地看到谱线长出锐利的峰。这一步做完你对光栅分辨本领正比于N这句话的感受会深很多。2.4 多光束干涉与法布里-珀罗腔的统一视角如果你继续深入会发现光栅的多光束干涉和F-P腔的多次反射干涉数学形式出奇地一致。F-P腔的透射率公式是T 1 / (1 F·sin²(δ/2))其中F 4R/(1-R)²是精细度系数δ是光在腔内往返一次的相位差。这里也是多光束叠加只不过每束光的振幅不是相等的而是经过反射镜后逐次衰减——振幅按R^(m-1)递减相位等差递增。所以做多光束干涉仿真时我建议你一次性把这两个模型都建了光栅模型等幅、等差相位对应sin(Nx)/sin(x)形式。F-P模型衰减振幅、等差相位对应Airy函数形式。两个模型对比着看,能清晰理解振幅均匀性如何影响条纹锐度。3. Matlab代码实现与参数选择3.1 光栅多光束干涉的完整代码先上代码再逐段解释。这是我的版本你直接复制就能跑% 多光束干涉仿真光栅模型 % 适合Matlab R2016b及以上版本 clear; close all; clc; %% 参数设置 lambda 632.8e-9; % 波长He-Ne激光 632.8nm d 50e-6; % 缝间距 50um N 6; % 光束数量缝数 k 2*pi/lambda; % 波数 % 观察角范围这里用弧度 % 主极大出现在 d*sinθ m*λ % 一阶极大角度 approx lambda/d 0.01266 rad theta linspace(-0.05, 0.05, 4001); theta_deg theta * 180/pi; %% 方法一等比数列求和公式 phi k * d * sin(theta); % 处理除零把分母极小化 epsilon 1e-15; U_formula exp(1i*(N-1)*phi/2) .* sin(N*phi/2) ./ (sin(phi/2) epsilon); I_formula abs(U_formula).^2 / N^2; %% 方法二逐束累加验证用 U_sum zeros(size(theta)); for n 0:N-1 U_sum U_sum exp(1i*n*phi); end I_sum abs(U_sum).^2 / N^2; %% 绘图 figure(Position, [100 100 900 400]); plot(theta_deg, I_formula, b-, LineWidth, 2); hold on; plot(theta_deg, I_sum, r--, LineWidth, 1); xlabel(观察角 θ (度), FontSize, 12); ylabel(归一化光强 I/I_0, FontSize, 12); title([N , num2str(N), 束光干涉, d , num2str(d*1e6), μm, λ , num2str(lambda*1e9), nm], FontSize, 12); legend({解析公式, 逐束累加}, Location, north); grid on; xlim([-3 3]); % 检查两种方法最大偏差 max_dev max(abs(I_formula - I_sum)); fprintf(解析法与累加法最大偏差: %.3e\n, max_dev);运行这段代码你会看到两条曲线几乎完全重合最大偏差在10的负15次方量级主要来自浮点误差和epsilon修正。这说明公式写对了代码逻辑也没问题。3.2 关键参数的选择逻辑与经验值参数怎么取这里有几个经验值供参考波长lambda用632.8nm氦氖激光最经典这是光学实验室最常见的激光波长。用别的波长也行但632.8nm的好处是很多教材的仿真图都用这个值方便对照。缝间距dd取50μm这样d/λ约为79一阶主极大在sinθ 0.0127附近观察角范围取±0.05弧度约±2.9度刚好能看到两三个主极大。如果你的电脑性能好theta阵列可以取4001个点曲线很光滑。如果跑不动2001个点也足够。光束数NN6是我推荐的教学参数。N2时看双缝干涉余弦条纹N6时已经能清晰看到4个次级峰N30以上时效应开始趋于光栅极限。N步进太大时主峰太尖锐肉眼反而不容易看细节。注意当N很大比如500以上逐束累加法的循环会变慢。这时请用解析公式或者把循环换成矩阵运算n_vec (0:N-1);% U_matrix exp(1in_vecphi); % I_matrix abs(sum(U_matrix, 1)).^2 / N^2;这种向量化写法在N很大时比for循环快好几个数量级。3.3 F-P腔仿真代码五分钟改出来如果你想继续做F-P腔代码改动非常小。核心是换成Airy分布% 法布里-珀罗干涉仪透射谱仿真 clear; close all; clc; R 0.9; % 镜面反射率 F 4*R / (1-R)^2; % 精细度系数 delta linspace(0, 4*pi, 4001); % 单程相位差 T 1 ./ (1 F * sin(delta/2).^2); figure(Position, [100 100 900 400]); plot(delta/pi, T, b-, LineWidth, 2); xlabel(单程相位差 δ/π, FontSize, 12); ylabel(透射率 T, FontSize, 12); title([F-P腔透射谱, R , num2str(R)], FontSize, 12); grid on; xlim([0 4]); ylim([0 1]);把R改成0.04、0.5、0.9、0.99对比透射峰的宽度变化你会直观看到反射率越高、精细度越大、峰越尖锐。这个峰越尖锐直接对应F-P腔在光谱分析中的高分辨率应用。3.4 扩展到二维干涉图样一维曲线看够了很多人会想画二维干涉图样。思路是把观察角从一维linspace扩展成二维meshgrid光强公式不变用imagesc或surf画出来。% 二维多光束干涉图样 theta_x linspace(-0.02, 0.02, 600); theta_y linspace(-0.02, 0.02, 600); [ThetaX, ThetaY] meshgrid(theta_x, theta_y); % 二维光栅x和y方向各有N列光源 phi_x k * d * sin(ThetaX); phi_y k * d * sin(ThetaY); U_x sin(N*phi_x/2) ./ sin(phi_x/2 1e-15); U_y sin(N*phi_y/2) ./ sin(phi_y/2 1e-15); I_2d (abs(U_x).^2 ./ N^2) .* (abs(U_y).^2 ./ N^2); figure; imagesc(theta_x*180/pi, theta_y*180/pi, I_2d); axis image; colormap(hot); colorbar; xlabel(θ_x (度)); ylabel(θ_y (度)); title(二维多光束干涉图样);跑出来你会看到整齐排列的亮点阵列——这就是二维光栅比如生物传感器的纳米孔阵列的远场衍射图样。如果你是做阵列天线方向图仿真这个二维版本也直接对应平面阵列的阵因子。4. 仿真结果分析与数据可视化4.1 主极大位置验证理论的铁证跑完仿真第一件事是验证主极大位置是否符合理论。从光栅方程d·sinθ m·λ一阶主极大m1处 sinθ λ/d 632.8e-9 / 50e-6 ≈ 0.01266对应角度0.725度。你在Matlab里找峰值坐标[~, locs] findpeaks(I_formula, MinPeakHeight, 0.5); theta_deg(locs)你会看到峰值出现在约±0.725度、±1.45度、±2.17度。这跟理论值完全吻合。这个验证很重要——它告诉你仿真代码没有犯方向性错误物理模型是自洽的。4.2 缝数N对条纹锐度的影响从余弦到刀锋把N从2扫到10、50、100光强分布会呈现如下特征缝数N主极大半宽近似次级峰数量条纹特征2较宽cos²形式0经典双缝条纹平缓6sin(Nx)/sin(x)型4主峰锐化次级峰可见30更锐28主峰细锐次级峰小而密100极锐98近似δ函数接近光栅极限主极大半宽大约等于2π/(N·d·cosθ)N越大峰越窄。这就是光栅分辨本领R mN的物理来源。4.3 次级峰的位置与缺级现象多光束干涉的光强分布中主极大之间出现了N-2个次级峰。这些次级峰的位置满足sin(N·delta/2) ±1 而 sin(delta/2) ≠ 0直观上这些峰是N个相位矢量部分抵消、未完全抵消的结果。次级峰的强度大约为主极大的1/(N·sin(3π/2N))²量级远小于主极大。更精彩的是缺级现象当你把单缝衍射的包络由缝宽a决定叠加上去后某些级次的主极大恰好落在单缝衍射的暗纹位置这一级就消失了。这时你需要把代码扩展为% 单缝衍射包络缝宽a a 20e-6; % 缝宽 single_slit (sin(k*a*sin(theta)/2) ./ (k*a*sin(theta)/2 1e-15)).^2; I_total single_slit .* I_formula; % 多光束干涉 × 单缝衍射叠加后绘制的曲线你会看到有的主极大被吃掉了。这个现象在光栅光谱仪设计中非常关键——选择适当的缝宽与缝间距比例可以抑制不需要的衍射级次提高信噪比。4.4 数据可视化经验让你的图能放进论文Matlab绘图有几个小技巧能让仿真结果更专业使用set(gca, FontName, Times New Roman)统一字体论文投稿时更美观。用yyaxis在一张图上同时显示线性坐标和对数坐标的光强分布主极大和次级峰都能看清。保存图像用exportgraphics(gcf, fig.png, Resolution, 300)避免低分辨率锯齿。% 对数坐标显示次级峰细节 figure; semilogy(theta_deg, I_formula, b-, LineWidth, 1.5); xlabel(观察角 θ (度)); ylabel(光强 I/I_0 (对数坐标)); title(多光束干涉光强分布对数坐标); grid on; xlim([0 3]); ylim([1e-6 1]);在普通线性坐标下次级峰几乎贴着零线看不见但在对数坐标下次级峰的精细结构一览无余。做光学仿真一定要善用对数坐标。5. 从多光束干涉到电扫阵列跨学科迁移5.1 阵列天线与光栅的同构关系如果你用Matlab做过电扫阵列建模与仿真一定会觉得光栅的多光束干涉公式眼熟。没错阵列天线的阵因子公式是AF(θ) Σ A_n · exp(i·n·k·d·sinθ)这就是多光束干涉的精确翻版。区别仅在于光学中n是缝的序号射频中是阵元的序号。光学中激励幅度通常均匀射频中可以做泰勒加权、切比雪夫加权。光学中关注光强分布射频中关注辐射方向图。基于多光束干涉代码改成阵列天线方向图仿真只需要% 电扫阵列因子仿真8元均匀线阵示例 N_elem 8; d_elem 0.5 * lambda_rf; % 半波长间距 theta_scan 30 * pi/180; % 波束指向30度 % 相邻阵元相位差补偿β -k*d*sin(θ_scan) beta -2*pi*d_elem/lambda_rf * sin(theta_scan); psi 2*pi*d_elem/lambda_rf * sin(theta) beta; AF sin(N_elem*psi/2) ./ (N_elem * sin(psi/2));跑出来就是电扫阵列的方向图主瓣在30度方向旁瓣、栅瓣条件与光栅完全一致。这就是为什么光学仿真能帮你理解射频阵列物理模型是相通的。5.2 主瓣宽度、栅瓣抑制与光栅色散的类比在阵列天线设计中有一个经典问题阵元间距超过半个波长时方向图会出现与主瓣同等幅度的栅瓣grating lobe。这在光学中就是多光束干涉的高阶主极大用你的光栅仿真代码把d从0.5λ逐步增大到1λ、1.5λ观察主极大间距变化d 0.5λ±1级主极大在sinθ ±2的位置不在可见空间-1到1之间无栅瓣。d 1λ±1级主极大恰好在sinθ ±1边缘有栅瓣危险。d 1.5λ±1级主极大出现在sinθ ±0.67明显的栅瓣进入可见空间。这个模拟直观得让人叹气——原来天线阵列中的栅瓣抑制和光栅衍射中高阶极大抑制是同一个物理问题。5.3 波束扫描与色散效应电扫阵列通过改变相邻阵元之间的激励相位差beta来偏转波束方向。这个相移扫描的光学类比就是改变入射光的角度或改变介质折射率两者都会导致主极大方向改变。在你的光栅仿真中加入一个初始相位偏移phi_0phi k * d * sin(theta) phi_0; % phi_0 是初始相位差观察phi_0从0扫到π时主极大位置的移动。这本质上就是相控阵的波束扫描动画。配合animatedline函数可以做出动态扫描效果figure; h animatedline(Color, b, LineWidth, 2); for phi0 linspace(0, pi, 50) clearpoints(h); phi k*d*sin(theta) phi0; I (sin(N*phi/2) ./ (N*sin(phi/2))).^2; addpoints(h, theta_deg, I); drawnow; pause(0.05); end你会看到主极大像探照灯一样扫过整个观察角范围。这个画面让我当年对相控阵的理解从公式变成了直觉。6. 常见问题与调试技巧实录6.1 结果出现NaN或Inf除零问题排查现象光强曲线在特定位置出现NaN或峰值无穷大。原因sin(phi/2)在phi 2πm处为零此时分子sin(N·phi/2)也为零0/0在浮点运算中得到NaN。解决% 代替直接除零使用epsilon修正 epsilon 1e-15; I_formula abs(sin(N*phi/2) ./ (sin(phi/2) epsilon)).^2 / N^2;这样在零点附近表达式会给出一个有限值约为N²除以N²后趋于1正好是主极大。6.2 曲线毛刺很多网格点数不足或过密现象曲线在次级峰处有锯齿状毛刺或用findpeaks找峰时出现虚假峰。原因theta网格点数不足导致峰值采样点太少形状失真。反过来点数过多时内存和计算量增加但一般不会出错。经验值角度范围±0.05弧度至少2001个点推荐4001。角度范围±0.5弧度大角度至少8001个点。如果使用了meshgrid二维扩展600×600网格是起步1200×1200需要较好内存。6.3 主极大强度不归一化N越大值越大现象N从6变成100光强峰值从1变成了100甚至上万。原因忘了除以N²。总复振幅U的幅度在主极大处为N光强为N²不归一化的话曲线纵坐标随N剧烈变化。解决在绘图前统一I I / max(I)或I I / N^2。推荐后者因为这样你看到的纵坐标直接对应物理上的功率占比。6.4 相位差公式用错角度单位与弧度陷阱现象结果看起来完全不对峰的位置离理论值差很远。原因Matlab的sin函数默认输入是弧度很多人直接输入角度值。解决始终用弧度计算显示时再转角度。或者在开头加一行theta_deg linspace(-3, 3, 4001); theta theta_deg * pi / 180;6.5 代码性能优化从for循环到向量化现象N很大比如1000时逐束累加的for循环要跑很久。解决优先用解析公式如果必须用逐束累加比如验证或自定义振幅用矩阵运算一次完成% 向量化累加N×M矩阵每行对应一束光 n_vec (0:N-1); U_matrix exp(1i * n_vec * phi); % phi 是 1×M 行向量 U_vec sum(U_matrix, 1); I_vec abs(U_vec).^2 / N^2;这个写法在N1000、M4001时依然流畅。6.6 F-P腔仿真中R1时公式爆炸现象F4R/(1-R)²在R1时分母为零透射率公式无穷大。解决实际镜面反射率永远小于1建议R的取值范围设为0到0.9999。如果你看到透射峰极度尖锐到出现数值溢出那就是R太接近1了适当调小即可。7. 实操心得与项目扩展建议这个仿真项目做完我的体会是多光束干涉不是一个孤立的知识点而是一张连接了光学、射频、信号处理的网络。你在Matlab里敲下的每一行代码都在把抽象的公式变成可视化的物理图像。如果你想继续扩展这里有几个方向供参考非均匀振幅把等幅叠加改成高斯振幅分布观察条纹副瓣抑制效果——这对应阵列天线中的泰勒加权。相位噪声给每束光的相位叠加随机扰动模拟真实光源的相干性下降观察条纹对比度退化。啁啾光栅让缝间距d随位置缓慢变化模拟啁啾光纤光栅的反射谱。动态动画把d或lambda做成滑动条用uicontrol做交互式参数调节变成一个小型虚拟光学实验室。我个人觉得做这类仿真最忌讳的就是跑通就跑通完全不看中间过程。多光束干涉的美妙之处在于它的每一步中间结果都有明确的物理含义——相位差对应光程叠加对应干涉平方对应能量。把这些含义和代码对应起来你对光的理解就不再是公式而是图像。这也是Matlab这类工具给我最大的帮助它不是替你思考而是帮你看清思考的对象。