基于传递矩阵法的FBG与DFBG仿真:Matlab实现与调试指南
2026/9/24 22:37:19 网站建设 项目流程

做光纤光栅仿真的朋友,大概率都绕不开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)

用前面的参数计算一下:如果dn=1e-4,λ=1550nm,κ约为π*1e-4/1.55e-6 ≈ 202.8 m⁻¹,若L=10mm,则κL≈2.028,tanh²(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(忽略损耗时),如果发现R+T明显偏离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_L1=15mm_L2=10mm_Lcav=20mm_dn=1p5e-4.png”。我吃过太多次亏,同一批仿真跑了几十条曲线,隔几天再回来看文件夹,根本分不清哪张图对应哪组参数。没有规范化命名,后续整理数据时整理到怀疑人生。

第三点,DFBG的仿真结果和实验对照时,不能只看峰值位置和带宽。因为实验光栅写入过程中的紫外光强波动、纤芯折射率不均匀,都会让谱线变得不再平滑。我的经验是先把两个光纤光栅分别做实验测试,拿到每个光栅的实际参数,再在仿真模型里使用这些“实测参数”,而不是理论标称值。这样才能把DFBG干涉条纹的细节对得起来。

做FBG和DFBG仿真,说到底是把物理概念变成数学模型,再让Matlab替我们把模型算出来。只要传递矩阵思路理解透,后面扩展到啁啾光栅、相移光栅、取样光栅,都只是换一下矩阵构造方式而已。希望这篇分享能帮你减少一些试错时间,动手写出一版稳定好用的仿真代码,然后在这个基础上加上你自己的应用场景。

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

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

立即咨询