做光纤光栅仿真的朋友大概率都绕不开Matlab。我第一次接触FBG仿真时也试过直接用商业光学软件点几个参数出反射谱确实快可一旦涉及多级光栅、非均匀切趾、双光栅腔型结构或者要在设计空间里反复扫描优化商业软件那套交互就会拖慢节奏。后来我自己用Matlab实现了FBG和DFBG双光纤光栅的仿真并基于传递矩阵法搭了一套可扩展的模型实测下来效果很稳调参和批量计算都灵活得多。这篇博文就把整套思路摊开讲从物理模型、仿真环境准备到代码实现、结果校验再到我实际踩过的坑和调试经验尽量让还没动手写过的朋友也能顺着走一遍。1. 仿真前先搞清FBG和DFBG在算什么1.1 FBG的物理模型与布拉格条件FBG本质上是在光纤纤芯内写入周期性折射率调制相当于一段天然滤波器。当宽带光进入光栅区域只有满足布拉格条件的光会被强烈反射其余波长继续传输。布拉格条件写出来很简单λ_B 2 * n_eff * Λ其中n_eff是纤芯基模的有效折射率Λ是光栅周期。这个公式是FBG所有特性的起点。仿真里我们其实在算不同波长下光在光栅里往返传输后的反射率R(λ)和透射率T(λ)所以核心任务就是建立折射率扰动δn(z)沿光栅长度方向的分布模型然后求电磁场在结构中的响应。对于均匀光栅δn(z)沿轴向是常数只在光栅区段内存在。若是切趾光栅或者啁啾光栅δn(z)或Λ(z)会随位置变化此时要分段甚至逐周期建模。绝大多数FBG仿真基于耦合模理论把前向模与后向模的振幅耦合方程写出来求解出反射系数。耦合系数κ的计算式为κ π * Δn_eff / λ这里的Δn_eff是有效折射率调制深度。在弱导近似下Δn_eff可取折射率调制的平均值实际计算时还要乘以条纹可见度系数一般为1或小于1的常数。这个系数直接影响反射率峰值大小仿真中调参时可以把它当作一个独立变量。1.2 DFBG的结构类型与仿真目标双光纤光栅DFBG这个缩写在不同文献里含义不完全一致。我在实际工作中遇到过两类常见指代一是在同一根光纤上先后写入两个不同中心波长的布拉格光栅用于双参数或多波长测量二是将两个相同参数的光栅间隔一定距离排列形成光纤法布里-珀罗F-P腔型结构这类结构在传感中常用于同时测量静态应变和动态振动。本文仿真重点放在第二类也就是两个均匀FBG级联、中间隔一段普通光纤的结构因为它的谱线调制特性更能体现“双光栅”的价值。级联双FBG的反射谱并不是两个独立FBG反射谱的简单叠加。布拉格波长处的反射光会在两个光栅之间的腔内多次往返产生干涉调制具体表现为反射谱在中心波长附近出现一系列等间隔的细小波纹或陷波。相邻腔模的波长间隔与腔长L有关近似为Δλ ≈ λ_B² / (2 * n_eff * L_cav)这里的L_cav对应两个光栅之间的物理间距加上光栅内部的等效穿透深度。实际仿真时如果我们把中间间隔光纤也建模成一段传播矩阵这个腔模间隔会自然出现在结果中不需要额外修正。若用两个不同中心波长的FBG串联反射谱会呈现两个独立峰每个峰的位置由各自布拉格条件决定互相干扰很小仿真模型与单FBG情况几乎一样只需把不同光栅的参数分别设置。1.3 为什么选Matlab做这套仿真市面上有OptiGrating、OptiWave、RSoft等专业光器件仿真软件它们内置了完整的耦合模求解器做单次仿真很方便。但选择Matlab构建自己的仿真模型有几点非常实际的优势代码透明每一步物理量都能打印出来核对结果数组可以直接和后处理脚本对接便于大规模参数扫描多光栅结构、非均匀切趾、温度应变梯度等自定义场景在Matlab里只是改参数和矩阵连接关系而在商业软件里往往受固定模板限制。另外做课题或写论文时Matlab代码通常比商业软件的截图更容易被复现和验证这也是很多研究组愿意把算法核心放在Matlab里的原因。仿真精度上传递矩阵法配合合理的分段数完全可以达到nm级甚至pm级反射谱精度足够绝大多数工程分析使用。2. Matlab仿真环境的搭建与常见小坑2.1 安装与配置要点Matlab的安装本身不复杂最稳妥的方式是到MathWorks官网下载对应版本安装包用学校或单位的教育License激活。如果你是在国内环境建议直接选用R2021b之后的版本这些版本对中文编码和系统字体支持更好。安装完成后我会第一时间做两件事一是确认“C/C编译器”可用因为后续如果涉及编译MEX代码加速仿真需要用到二是预设一个固定的工作目录用“cd”命令把默认路径切换到项目文件夹避免后续代码里出现相对路径混乱。对于光纤光栅仿真这种计算密度较高的任务Matlab版本之间的性能差异并不明显真的瓶颈在于循环写法。我电脑上从R2018a到R2023a都用过同一份传递矩阵代码在R2023a上并没有快很多反而是改成向量化和预分配数组后提速了好几倍。因此别把性能希望寄托在版本升级上关键是代码本身。2.2 中文注释乱码与编码问题很多朋友在Matlab里写中文注释换一个机器打开就变成乱码这大概率是文件编码不一致导致的。早期Matlab默认以系统本地编码保存.m文件Windows下通常是GBK而Linux或Mac下常是UTF-8。跨平台打开时如果识别错编码中文注释就会乱码。这个问题从R2020b左右开始有所改善但仍不是百分百可靠。我的建议是一律采用英文注释或者将Matlab编辑器预设的编码格式统一为UTF-8。如果你的项目里已经存在大量GBK编码的旧文件可以通过“主页-预设-编辑器/调试器-语言”调整编码优先级再把旧代码另存为UTF-8格式。实操中如果想保存为UTF-8可以在编辑器中执行“文件-另存为-编码为UTF-8”这样代码里带中文注释就不会在跨平台时出乱码了。更省心的做法是只用英文变量和英文注释这个词在社区里能少一大半。2.3 代码组织与工程习惯仿真代码刚开始写的时候可以是一长条脚本但一旦进入参数扫描阶段脚本复用性会很差。我习惯把模型拆成三层参数定义脚本struct或单独的.m文件核心计算函数以及绘图与数据导出脚本。例如建立一个fbg_sim_params.m存放所有物理参数一个calc_fbg_spectrum.m负责矩阵计算一个main_run.m负责调用并出图。这样做的好处是换一组参数时只需要改参数脚本计算函数完全不动避免在脚本里来回改数值导致记录混乱。另外强烈建议在代码里加入“单位说明”。FBG仿真中长度单位常用nm表示波长mm表示光栅长度折射率调制深度Δn_eff通常写作1e-4到1e-3量级。每一次参数赋值都写上注释比如% L_cm in mm能够避免日常写作时把1e-3误写成1e-4这种低级错误。3. 传递矩阵模型搭建和代码讲解3.1 传递矩阵法的数学原理与分段思路传递矩阵法把光栅沿轴向划分成M段每一段近似为均匀周期光栅或均匀波导再用一个2x2矩阵描述该段的输入输出关系最终将M段矩阵依次相乘得到整个结构的传输矩阵。相比直接求解耦合微分方程传递矩阵法更直观也更适合处理级联光栅和中间插入普通光纤段的情况。对于一段长度为dz、耦合系数为κ、相位失配为Δβ的均匀光栅其传递矩阵T可表示为T [ cosh(γdz) - j(Δβ/γ)sinh(γdz), -j*(κ/γ)sinh(γdz); j*(κ/γ)sinh(γdz), cosh(γdz) j(Δβ/γ)sinh(γdz) ]其中γ² κ² - Δβ²Δβ β - π/Λβ 2π * n_eff / λ。要注意矩阵元素里j的正负号顺序不同教材的符号约定可能不同最终结果只要自洽就行。我习惯采用上述定义计算反射系数时使用r -T21/T22透射系数t 1/T22反射率R |r|²透射率T |t|²。如果计算的是两段均匀光栅之间的普通光纤段没有调制κ0γ j*Δβ矩阵退化为对角传播矩阵P [ exp(-jβL_fiber), 0; 0, exp(jβL_fiber) ]这里的β是普通光纤中基模的传播常数L_fiber是两个光栅之间的光纤长度。这个传播矩阵是DFBG腔型结构仿真的关键。分段数M的选择直接影响精度。M取得太小每个子段的近似误差会被放大尤其是高反射率光栅反射谱可能出现振荡或峰值偏移M取得太大计算时间成倍上升。经过实测均匀FBG取M在200到500之间就能在实验室常见参数下获得平滑光谱。如果做啁啾光栅M需要增加到1000以上因为周期在变化每段内的均匀近似依赖细小的划分。3.2 单FBG的仿真代码详解下面这段代码是我常用的均匀FBG反射谱计算函数包含注释可以直接复制到.m文件中运行。function [lambda, R, T] calc_uniform_fbg(neff, L, Lambda, dn, n_seg, lambda_range) % 均匀FBG传递矩阵仿真 % 输入: % neff - 有效折射率 % L - 光栅长度 (m) % Lambda - 光栅周期 (m) % dn - 折射率调制深度 % n_seg- 分段数 % lambda_range - [lambda_min, lambda_max] (m) % 输出: % lambda, R, T - 波长向量、反射率、透射率 lambda_min lambda_range(1); lambda_max lambda_range(2); lambda linspace(lambda_min, lambda_max, 5000); % 波长扫描点数 dz L / n_seg; % 每段长度 kappa pi * dn / lambda; % 耦合系数随波长变化 % 预分配 R zeros(size(lambda)); T zeros(size(lambda)); for idx 1:length(lambda) lam lambda(idx); beta 2 * pi * neff / lam; delta_beta beta - pi / Lambda; gamma sqrt(kappa(idx)^2 - delta_beta.^2); % 对于|gamma|很小的情况用数值近似避免sinh/cosh参数虚部过大 M eye(2); for seg 1:n_seg if abs(gamma) 1e-9 T_seg [1 - 1j*delta_beta*dz, -1j*kappa(idx)*dz; -1j*kappa(idx)*dz, 1 1j*delta_beta*dz]; else cosh_g cosh(gamma*dz); sinh_g sinh(gamma*dz); T_seg [cosh_g - 1j*delta_beta/gamma*sinh_g, ... -1j*kappa(idx)/gamma*sinh_g; 1j*kappa(idx)/gamma*sinh_g, ... cosh_g 1j*delta_beta/gamma*sinh_g]; end M M * T_seg; end r -M(2,1)/M(2,2); t 1/M(2,2); R(idx) abs(r)^2; T(idx) abs(t)^2; end end这段代码的逻辑很直观外层循环扫描波长内层循环把整段光栅串起来。实际运行中如果只算一条反射谱5000个波长点配合400个分段在普通电脑上大约需要几秒到十几秒。如果要频繁扫描参数建议先降低波长点数到2000左右快速观察趋势再提高精度做最终确认。一个容易忽视的细节是当波长刚好在布拉格中心附近时delta_beta接近0此时γ接近κ矩阵中的cosh和sinh参数不再是纯虚数反射率接近峰值。若delta_beta比较大γ会变为虚数cosh(gamma*dz)等于cos(|gamma|*dz)此时光谱会出现旁瓣结构。这些数值行为正是反射谱纹理的来源理解后能帮你判断代码哪里写错了。3.3 级联DFBG的建模与仿真级联DFBG就是在单FBG仿真基础上把两个光栅的传递矩阵与中间光纤段的传播矩阵乘起来。整体结构可以表示为M_total M_FBG2 * P_fiber * M_FBG1这里的矩阵乘法顺序表示光从左向右穿过结构光栅1在左光栅2在右中间是一段普通光纤。注意我在代码里让光栅1作为入射端最终反射系数仍用r -M_total(2,1)/M_total(2,2)。如果两个光栅参数不同需要分别计算两个FBG段的矩阵。DFBG仿真代码片段如下% 级联双FBG (DFBG) 示例 % 两个光栅参数相同长度分别为15mm和10mm中间间隔20mm neff 1.46; Lambda 0.535e-6; % 对应中心波长约1562nm dn 1.5e-4; L1 15e-3; % 第一个光栅长度 L2 10e-3; % 第二个光栅长度 L_fiber 20e-3; % 中间光纤长度 n_seg1 400; n_seg2 300; n_seg_fiber 20; lambda linspace(1560e-9, 1564e-9, 10000); R_total zeros(size(lambda)); T_total zeros(size(lambda)); for idx 1:length(lambda) lam lambda(idx); kappa pi * dn / lam; beta 2 * pi * neff / lam; M1 fbgn_tmatrix(neff, L1, Lambda, kappa, beta, n_seg1); M2 fbgn_tmatrix(neff, L2, Lambda, kappa, beta, n_seg2); % 中间光纤段的传播矩阵 P [exp(-1j*beta*L_fiber), 0; 0, exp(1j*beta*L_fiber)]; M_total M2 * P * M1; r -M_total(2,1) / M_total(2,2); t 1 / M_total(2,2); R_total(idx) abs(r)^2; T_total(idx) abs(t)^2; end这里的fbgn_tmatrix是一个辅助函数内部就是根据给定参数计算一段均匀FBG的2x2矩阵。如果你只把两个光栅的矩阵简单相乘而不插入中间的传播矩阵得到的结果是两个独立FBG反射率的叠加完全没有腔模信息。加上了P后反射谱中会自然出现周期性调制条纹。细心的读者会发现我把波长扫描点数提高到10000因为腔模的条纹间隔通常很窄如果扫描点太稀条纹会被直接跳过结果看起来只是一条平滑曲线从而误认为双光栅没起作用。3.4 仿真结果的物理校验仿真跑完后第一件要做的不是调参美化而是校验结果是否物理合理性。对均匀FBG可以对比理论峰值反射率。均匀光栅在布拉格波长处的峰值反射率可用公式计算R_peak tanh²(κL)用前面的参数计算一下如果dn1e-4λ1550nmκ约为π*1e-4/1.55e-6 ≈ 202.8 m⁻¹若L10mm则κL≈2.028tanh²(2.028)≈0.937。仿真结果应该非常接近这个值。如果仿真得到的反射率明显偏离大概率是分段数太少或矩阵元素符号写错了。对于级联DFBG可以验证腔模间隔与理论公式是否一致。假设两个光栅之间物理间隔20mm有效折射率1.46中心波长1562nm那么考虑到光栅内部有效穿透深度大约为1/κ量级实际等效腔长要比物理间隔略大一些。我用传递矩阵算出的条纹间隔通常在几十pm量级而用简化公式计算的Δλ ≈ λ²/(2nL)也在这个范围内两者可以互相印证。如果条纹间隔完全对不上就要检查有没有把光栅内部带来的穿透深度考虑进去。4. 参数选择、常见错误和提速技巧4.1 参数选择的“度”仿真参数不是越多越好关键是找到“平衡点”。波长扫描范围要覆盖光栅反射谱的主要部分。普通切趾FBG的半峰宽可能只有0.1~0.5nm如果扫描范围设成全波段谱线根本看不出细节反之如果只扫中心附近0.1nm又必然错过旁瓣和DFBG的条纹调制。我的经验是先用宽范围快速粗扫比如中心波长±5nm步长10pm看一下整体形态和峰值位置然后再加密到步长1pm只扫中心波长±1nm的区域。这样既不会浪费算力也不会丢失关键特征。分段数M需要和折射率调制深度、光栅长度协同考虑。对于弱调制光栅dn~1e-5量级400段已经足够对于强调制光栅dn~1e-3量级反射率很高能量在光栅内部快速交换分段不足会导致矩阵元素累计误差反射率曲线在带内出现不合理的振荡。遇到这种情况把n_seg提高到1000观察结果是否稳定如果谱线变化不大说明精度已收敛。判断标准很简单同一参数跑两次一次500段、一次1000段两条曲线完全重合就认定结果可靠。4.2 常见报错与不收敛问题仿真过程中最常见的报错是“Matrix dimensions must agree”或“Subscript indices must either be real positive integers”。前者通常是beta、kappa和lambda三个向量的尺寸对不上比如计算kappa时用了整个lambda向量但在内层循环里lambda是标量忘了重新赋值。后者常见于不小心用idx索引到了非整数数组尤其是当你把波长扫描点数写成linspace(1560e-9,1564e-9,10000)后又试图用得到的小数以数组索引访问数据。这个问题我早期几乎每周都会遇到解决方法是尽量避免在循环体内使用复数计算得到的浮点数作为索引所有索引都用整数变量。另一个更隐蔽的问题是反射率不为正或在某些波段超过1。反射率理论范围是0到1如果算出来有负数或超过1原因多半是矩阵构造时符号写反。常见错误出现在delta_beta的定义上有些教材写成β - π/Λ有些写成π/Λ - β两者会改变传递矩阵的共轭对称性。不同写法最后算出的反射率其实一致但中途矩阵的虚部符号会有差别。如果不小心把sinh项前面的j放错位置就会导致能量不守恒。此时检查准则透射率加反射率应恒等于1忽略损耗时如果发现RT明显偏离1优先级最高的排查方向就是矩阵元素符号。4.3 提速与批量扫描技巧如果只是做单条光谱仿真上述代码足够快。但如果要扫描切趾参数、光栅长度、折射率调制深度这三种变量每条曲线10000个波长点叠加1000段分段循环次数就会达到千万甚至亿级速度会变得很难看。这时需要三步优化。第一步向量化波长扫描。矩阵计算里的很多操作其实是逐元素运算可以写成向量形式但传递矩阵每段相乘时必须按波长方向迭代很难完全向量化。不过对于线性问题可以采用差分方法或状态空间近似不过实现复杂度较高。第二步使用parfor并行。把外面一层的for idx 1:length(lambda)改成parfor idx 1:length(lambda)就能利用多核CPU同时计算不同波长点的传递矩阵。要注意的是parfor循环体内不能写全局变量或随机数否则结果可能不一致。第三步预分配数组。不要用动态增长的方式添加R和T在循环前用zeros预先分配好否则内存重新分配的开销会拖慢整体速度。我实测过一组参数波长点数5000分段数400单次仿真约8秒。改成parfor并行并同时分配到4个workder后耗时降到3秒左右。再对中间结果做缓存如果很多次仿真的光栅参数完全相同只是波长范围不同那就把对应光栅的传递矩阵缓存下来避免重复计算。4.4 和实验对照时的关键点仿真和实验对照时最容易出问题的地方不是反射率峰值而是旁瓣结构。实际光纤光栅写入过程中折射率调制的上升沿和下降沿不是理想矩形光栅两端会有渐变区域。如果仿真里硬用均匀折射率突变旁瓣幅度往往看起来比实验高谱线也显得更“脏”。这时候需要在模型中引入两端渐变段通常取上升沿长度约为总长度的5%~10%用高斯或升余弦函数过渡。仿真结果会和实验更接近。DFBG仿真中另一个常被忽略的细节是光纤的双折射。若用于传感的光栅是写在保偏光纤上的两个正交偏振方向的等效折射率不同布拉格波长也会分裂为两个峰。严格来说这种情况要用各向异性模型但很多情况下可以先当作两个独立FBG分别仿真再把结果做矢量叠加。这样处理虽然精度有限但能快速估计双峰间距和偏振串扰水平。若想更精细需要在耦合模方程中加入偏振项代码复杂度会上一个台阶。5. 我在实际仿真中积累的几点体会说了这么多最后聊点实操层面踩过的坑。第一点永远不要先调参数再检查物理合理性。我无数次先把光栅长度设成15mm又改了8mm结果反射峰位置整体偏移最后才发现是周期和波长没对应上。仿真第一步永远是拿简单的已知条件去验证代码比如弱光栅的峰值反射率对照tanh²(κL)公式偏差小于1%再开始正式调参。第二点保存仿真脚本时把关键参数和结果图像名同步写进文件名。比如“DFBG_L115mm_L210mm_Lcav20mm_dn1p5e-4.png”。我吃过太多次亏同一批仿真跑了几十条曲线隔几天再回来看文件夹根本分不清哪张图对应哪组参数。没有规范化命名后续整理数据时整理到怀疑人生。第三点DFBG的仿真结果和实验对照时不能只看峰值位置和带宽。因为实验光栅写入过程中的紫外光强波动、纤芯折射率不均匀都会让谱线变得不再平滑。我的经验是先把两个光纤光栅分别做实验测试拿到每个光栅的实际参数再在仿真模型里使用这些“实测参数”而不是理论标称值。这样才能把DFBG干涉条纹的细节对得起来。做FBG和DFBG仿真说到底是把物理概念变成数学模型再让Matlab替我们把模型算出来。只要传递矩阵思路理解透后面扩展到啁啾光栅、相移光栅、取样光栅都只是换一下矩阵构造方式而已。希望这篇分享能帮你减少一些试错时间动手写出一版稳定好用的仿真代码然后在这个基础上加上你自己的应用场景。
