做光纤传感或者光通信相关的仿真多半绕不开光纤布拉格光栅FBG。均匀光栅作为最基础的一类它的反射谱和透射谱几乎可以看作整个光栅理论的地基。今天这篇文章我就用MATLAB把均匀光纤布拉格光栅的传输矩阵法完整实现一遍——从矩阵怎么推导、参数怎么选到代码怎么写、结果怎么读一次性讲清楚。无论你是刚接触光栅仿真的学生还是想快速验证方案的工程师这套代码和思路都能直接拿走用。顺便说一句我最早做FBG仿真也是用这个思路到现在做光纤传感标定时还经常翻出这套代码改一改算是压箱底的工具了。1. 均匀光栅为什么要用传输矩阵法1.1 先搞明白均匀FBG到底在算什么光纤布拉格光栅简单说就是在光纤芯层写入一段周期性折射率调制周期通常在几百纳米量级。光在光纤里传播时遇到这种周期结构会发生布拉格反射满足相位匹配条件的波长被反射回来其他波长继续透射。中心反射波长由布拉格条件决定[ \lambda_B 2 n_{eff} \Lambda ]这里的 (n_{eff}) 是光纤有效折射率(\Lambda) 是光栅周期。举个例子普通单模光纤 (n_{eff}1.45)想做1550 nm的中心波长周期就是 (1550/(2 \times 1.45) \approx 534.5) nm。这个量级比可见光波长还小工艺上并不容易写但仿真里就是一行代码的事。均匀FBG的定义是全段光栅周期和折射率调制幅度都不变。它的反射谱最大值就在 (\lambda_B) 处谱形两侧会出现一串对称的旁瓣像 sinc 函数那样的起伏。为什么会出现旁瓣因为有限长光栅本身相当于一个有限长的周期性扰动傅里叶变换后必然带旁瓣。这个特征既是理解FBG光谱的钥匙也是后面做切趾优化的动机。从数值仿真角度看真正需要求解的是前向模和后向模沿着光栅长度的耦合过程。光栅内部每个位置前向波一边前进一边把部分能量转换到后向波后向波同样也在转换前向波。对这个过程建模最经典的路径有两个直接解耦合模方程或者用传输矩阵法Transfer Matrix MethodTMM把问题离散化。1.2 传输矩阵法比直接解耦合模方程好在哪里耦合模方程对均匀光栅其实是有解析解的所以有人会问既然有解析解为什么还要传输矩阵法这个问题的答案很直接解析解虽然漂亮但换到啁啾光栅、切趾光栅、相移光栅这些非均匀结构时立刻失效。而传输矩阵法的思路是“分段解析”——把光栅沿轴向切成很多小段每一段都近似看成均匀光栅用解析矩阵描述它的输入输出关系再把所有小段矩阵按顺序乘起来就得到整段光栅的响应。这个思路最大的优势是通用性强。均匀光栅只是每段参数相同的特例把代码里的周期、折射率调制改成向量立刻就能仿真啁啾光栅把某些段的参数做特殊处理就能仿真相移光栅。也就是说你花半小时把传输矩阵法写明白后面各种花式光栅仿真都是在这套代码上改参数不用推倒重来。另外传输矩阵法对分段数的要求很宽松。每段的物理长度不需要小到和光栅周期一个量级因为段内用的是均匀光栅的解析解不是简单的相位累积。只要段内参数近似均匀分段数取几百到几千就足够精确。对均匀光栅来说分段矩阵连乘的结果和整段解析解完全一致这也方便我们用理论值校验代码正确性。2. 传输矩阵法的核心公式与参数换算2.1 2×2矩阵每一段光栅用四个元素描述传输矩阵法的基本单元是一个2×2复矩阵描述一段长度为 (dz) 的均匀光栅前后两端的前向波 (A(z)) 和后向波 (B(z)) 之间的关系。我习惯写成下面的形式[ \begin{bmatrix} A(zdz) \ B(zdz) \end{bmatrix}\begin{bmatrix} T_{11} T_{12} \ T_{21} T_{22} \end{bmatrix} \begin{bmatrix} A(z) \ B(z) \end{bmatrix} ]对于均匀FBG一段矩阵元素是[ T_{11} \cosh(\gamma dz) - i\frac{\hat{\sigma}}{\gamma}\sinh(\gamma dz) ][ T_{12} -i\frac{\kappa}{\gamma}\sinh(\gamma dz) ][ T_{21} i\frac{\kappa}{\gamma}\sinh(\gamma dz) ][ T_{22} \cosh(\gamma dz) i\frac{\hat{\sigma}}{\gamma}\sinh(\gamma dz) ]其中 (\kappa) 是交流耦合系数描述折射率周期性调制造成的前后向波耦合强度(\hat{\sigma}) 是直流自耦合系数描述波长偏离布拉格条件时的相位失配(\gamma \sqrt{\kappa^2 - \hat{\sigma}^2})。这两个系数的表达式是[ \kappa \frac{\pi \Delta n v}{\lambda_B} ][ \hat{\sigma} 2\pi n_{eff}\left(\frac{1}{\lambda} - \frac{1}{\lambda_B}\right) ]这里 (\Delta n) 是折射率调制幅度(v) 是条纹可见度通常取1。注意 (\hat{\sigma}) 在波长等于布拉格波长时为零此时前后向波完全相位匹配反射最强波长偏离越远(\hat{\sigma}) 越大失配越严重反射率急剧下降。这就是FBG反射谱为什么会有一个窄带主峰的根本原因。有些教材在直流自耦合系数里还会加一项平均折射率增量导致的失配项。如果写光栅时折射率调制里带有直流分量 (\delta n_{mean})公式变成 (\hat{\sigma} 2\pi n_{eff}(1/\lambda - 1/\lambda_B) 2\pi \delta n_{mean}/\lambda)。均匀光栅仿真里我通常先设 (\delta n_{mean}0)跑通了再考虑这个细节。2.2 参数换算从设计波长到光栅周期仿真开始前一定要把所有物理量统一到国际单位。这个坑我踩过太多次把波长写成纳米直接代入公式结果是反射谱主峰位置完全对不上查半天代码才发现是单位问题。建议的参数初值如下参数符号建议初值说明有效折射率(n_{eff})1.45普通单模光纤典型值中心波长(\lambda_B)1550 nmC波段常用光栅周期(\Lambda)534.48 nm由 (\lambda_B/(2n_{eff})) 计算光栅长度(L)10 mm影响带宽和反射率折射率调制(\Delta n)5e-5典型范围1e-5到1e-4条纹可见度(v)1理想正弦调制分段数(M)500数值仿真离散段数周期必须由中心波长反算不能随便拍脑袋填。如果你想要仿真温度或者应变传感直接把 (\lambda_B) 或 (\Lambda) 改成移动后的值就行比如温度升高使光栅周期变大中心波长大约按10 pm/℃的量级往长波方向漂移。2.3 边界条件与反射率公式怎么来的整段光栅的总传输矩阵 (T_{total}) 是各段矩阵按顺序相乘的结果。有了总矩阵还需要正确的边界条件才能算出反射率和透射率。物理上光从光栅左侧 (z0) 入射所以前向波初始值 (A(0)1)光栅右侧 (zL) 之后没有反向波返回所以 (B(L)0)。这是TMM最容易出错的地方。写成矩阵方程就是[ \begin{bmatrix} A(L) \ 0 \end{bmatrix}\begin{bmatrix} T_{11} T_{12} \ T_{21} T_{22} \end{bmatrix} \begin{bmatrix} 1 \ B(0) \end{bmatrix} ]第二个等式展开(B(L) T_{21} \cdot 1 T_{22} \cdot B(0) 0)所以反射系数[ r \frac{B(0)}{A(0)} -\frac{T_{21}}{T_{22}} ]反射率 (R |r|^2)无损耗情况下透射率 (T 1 - R)。这个公式推导虽然简单但如果你在代码里把 (-T_{21}/T_{22}) 写成了 (T_{12}/T_{22})光谱形态会完全错乱。我第一次写就是这里搞反了反射谱看起来像一条诡异的曲线排查了两个小时才发现。3. MATLAB代码实现从零搭一个均匀FBG光谱仿真器3.1 主程序参数声明、波长扫描、绘图下面这段MATLAB代码可以直接复制运行。主程序负责设置参数、产生波长扫描点、调用子函数计算反射率和透射率最后把光谱画出来。clear; clc; close all; % 光栅参数 neff 1.45; % 有效折射率 lambda_B 1550e-9; % 设计中心波长 1550 nm Lambda lambda_B / (2 * neff); % 光栅周期约 534.48 nm L 10e-3; % 光栅长度 10 mm dn 5e-5; % 折射率调制幅度 v 1; % 条纹可见度 M 500; % 分段数 % 波长扫描范围 lambda_list linspace(1545e-9, 1555e-9, 2001); % 预分配数组 R_list zeros(size(lambda_list)); T_list zeros(size(lambda_list)); % 逐波长计算 for i 1:length(lambda_list) [R_list(i), T_list(i)] uniformFBG_TMM( ... lambda_list(i), neff, Lambda, L, dn, v, M); end % 绘图 figure; plot(lambda_list * 1e9, R_list, LineWidth, 1.5); hold on; plot(lambda_list * 1e9, T_list, LineWidth, 1.5); xlabel(波长 (nm)); ylabel(反射率 / 透射率); legend(反射率 R, 透射率 T); title(均匀光纤布拉格光栅光谱); grid on; xlim([1545, 1555]);这一段主程序非常直白。2001个波长点每个点做500次矩阵乘法整个计算在普通电脑上只需要一两秒。如果你只需要看主峰附近的细节可以把扫描范围缩小到1549到1551 nm点数保持2001分辨率会更高。3.2 核心子函数单波长传输矩阵实现真正干活的是这个子函数。输入一个波长和光栅参数输出该波长下的反射率和透射率。为了以后扩展我把所有参数都显式传进去没有使用全局变量。function [R, T] uniformFBG_TMM(lambda, neff, Lambda, L, dn, v, M) % 传输矩阵法计算均匀FBG反射率与透射率 % 输入 % lambda : 当前波长单位米 % neff : 有效折射率 % Lambda : 光栅周期单位米 % L : 光栅长度单位米 % dn : 折射率调制幅度 % v : 条纹可见度 % M : 分段数 % 输出 % R : 反射率 % T : 透射率 dz L / M; % 每段长度 lambda_D 2 * neff * Lambda; % 布拉格波长 % 交流耦合系数和直流自耦合系数 kappa pi * dn * v / lambda_D; sigma_dc 2 * pi * neff * (1 / lambda - 1 / lambda_D); gamma sqrt(kappa^2 - sigma_dc^2); % 单段传输矩阵 Tseg zeros(2, 2); Tseg(1,1) cosh(gamma * dz) - 1i * sigma_dc / gamma * sinh(gamma * dz); Tseg(1,2) -1i * kappa / gamma * sinh(gamma * dz); Tseg(2,1) 1i * kappa / gamma * sinh(gamma * dz); Tseg(2,2) cosh(gamma * dz) 1i * sigma_dc / gamma * sinh(gamma * dz); % 所有段矩阵连乘 Ttotal eye(2); for m 1:M Ttotal Ttotal * Tseg; end % 边界条件 A(0)1, B(L)0 r -Ttotal(2,1) / Ttotal(2,2); R abs(r)^2; T 1 - R; end这里有一个容易困惑的点(\gamma) 在波长远离布拉格条件时是纯虚数(\sinh(\gamma dz)) 变成 (\sin(|\gamma| dz)) 的形式。MATLAB的内置复数双曲函数会自动处理这种情况所以不用手动分情况写。第一次写代码时我为了“优化性能”用简单判断分开算实数虚数结果反而引入了很多边界问题后来直接交给MATLAB的复数运算代码简洁还稳定。另外虽然均匀光栅每段矩阵完全相同按说可以写成Ttotal Tseg^M但我还是建议用循环连乘。原因有两个一是以后改成啁啾光栅时每段矩阵本来就不同循环结构可以直接复用二是循环连乘的代码更贴近“传输矩阵”的物理含义别人看你代码时更容易理解。3.3 批量扫描折射率调制dn对光谱的影响单项参数扫描是仿真的日常操作。比如我想看不同折射率调制幅度下反射谱怎么变化只需要在外层再加一个循环。下面这段代码可以作为主程序的一部分追加运行用子图形式对比四条曲线。dn_list [1e-5, 2e-5, 5e-5, 1e-4]; figure; for j 1:length(dn_list) dn dn_list(j); for i 1:length(lambda_list) [R_list(i), T_list(i)] uniformFBG_TMM( ... lambda_list(i), neff, Lambda, L, dn, v, M); end subplot(2, 2, j); plot(lambda_list * 1e9, R_list, LineWidth, 1.5); xlabel(波长 (nm)); ylabel(反射率); title([dn , num2str(dn)]); grid on; xlim([1545, 1555]); end这个批量扫描跑起来也很快四个子图十来秒就出来了。你会看到随着 (\Delta n) 增大主峰反射率升高、宽度变宽旁瓣幅度也跟着抬升。这个趋势对光栅设计和传感方案选型非常关键。4. 仿真结果解读与光谱特征4.1 主反射峰与透射凹陷默认参数 (L10) mm、(\Delta n5 \times 10^{-5}) 跑出来的结果主峰应该准确落在1550 nm峰值反射率大约在0.58附近。为什么是0.58左右均匀FBG在布拉格波长处的理论最大反射率可以写成[ R_{max} \tanh^2(\kappa L) ]把 (\kappa \pi \times 5 \times 10^{-5} / 1550 \times 10^{-9}) 算出来大概是 (101.4) 每米(\kappa L \approx 1.014)所以 (R_{max} \approx \tanh^2(1.014) \approx 0.58)。把这个理论值和仿真输出对照一下如果对得上说明代码基本没问题。这个校验方法非常实用我在写任何新脚本时都会先用中心波长反射率做一次自检。透射谱在1550 nm处出现一个凹陷凹陷深度刚好是反射峰的高度。因为无损耗光栅满足能量守恒(R T 1)透射凹陷越深说明反射越强。如果你的透射率不是 (1-R)而是也通过传输矩阵直接推出来的那结果应该一致但实际项目里用 (T1-R) 就够了省时省力。4.2 旁瓣与带宽均匀光栅的“签名”均匀FBG反射谱最显著的特征除了主峰之外就是两侧的旁瓣。这些旁瓣是有限长周期性结构固有的来源可以粗浅理解为光栅长度有限两端像“截断窗口”对折射率调制的空间傅里叶变换产生了旁瓣。主峰两侧第一对旁瓣的幅度大约在主峰的十分之一量级按dB算大概相差10 dB左右。线性坐标下这些旁瓣看起来像是主峰两侧的小鼓包很多人会忽略。但如果把纵轴换成dB坐标旁瓣结构会变得非常清楚。我在实际分析中几乎都是看 (10\log_{10}(R)) 曲线的。下面这段代码可以追加在主程序后快速切换dB显示RdB 10 * log10(R_list 1e-12); % 加一个小量防止 log10(0) figure; plot(lambda_list * 1e9, RdB, LineWidth, 1.5); xlabel(波长 (nm)); ylabel(反射率 (dB)); title(均匀FBG反射谱dB坐标); grid on; xlim([1545, 1555]); ylim([-60, 0]);从dB曲线能直观看到主瓣宽度、旁瓣位置和深度也方便和实验测到的光谱做对比。做WDM滤波器的人最讨厌这些旁瓣因为它们会串扰邻近信道做传感器的人反而对主瓣和旁瓣的相对位置不敏感更关心峰值波长漂移。同样是这套光谱不同应用场景的关注点完全不同。如果要定量提取主瓣带宽可以用一条简单的逻辑找半高位置lambda_nm lambda_list * 1e9; half_idx R_list 0.5 * max(R_list); fwhm_nm lambda_nm(half_idx(end)) - lambda_nm(half_idx(1)); fprintf(3dB带宽约 %.4f nm\n, fwhm_nm);这个FWHM值会随 (\Delta n) 和 (L) 变化。我做参数优化时经常把这个提取逻辑和批量扫描结合起来自动找出满足带宽要求的光栅参数组合。4.3 dB坐标下旁瓣的真实面孔前面提过用dB坐标看旁瓣这里我再多说一句。均匀FBG线性光谱里旁瓣高度看起来很低好像只比基底高一点点但一旦换成dB你会发现第一旁瓣其实只比主峰低10 dB左右并没有想象中那么“弱”。这个认知对实际滤波器设计很重要——想要把旁瓣压到-30 dB以下均匀光栅根本做不到必须做切趾。切趾的思路很简单把折射率调制幅度 (\Delta n) 在光栅两端逐渐减小到零中间维持最大最常用的是高斯窗、余弦窗或者升余弦窗。从傅里叶变换的角度理解均匀光栅相当于矩形窗频域旁瓣高加窗之后等效窗函数变平滑旁瓣自然就被压低了。代价是主瓣变宽、峰值反射率略微下降。这套取舍在通信滤波器设计里尤其关键。4.4 中心波长附近的规律与传感启发把扫描范围缩小到1549到1551 nm你会看到反射谱在中心波长附近几乎是关于 (\lambda_B) 对称的。这个对称性来自均匀FBG结构的对称性——折射率调制在轴向是对称的反射谱自然对称。如果哪天你仿真出来的均匀FBG光谱不对称先不要怀疑物理大概率是代码里某个符号写错了。从传感角度看FBG最迷人的一点是环境温度、应变会改变 (n_{eff}) 和 (\Lambda)直接移动中心波长。温度升高10℃中心波长大约往长波方向漂移0.1 nm左右应变量变1个微应变更接近1 pm的漂移。所以哪怕只看中心波长的移动就能反推外界物理量。仿真里模拟这个现象非常方便比如把 (\Lambda) 变大0.1%重新跑一遍你会发现整个光谱平移形状几乎不变。我用这套代码做温度标定前置仿真时就是靠这个办法预估探测器分辨率够不够用。5. 常见问题与排查技巧实录5.1 反射率全零或全是NaN先查单位这是新手最容易栽的坑。如果你的反射率全部都是0先看一眼是不是把波长以纳米为单位直接代入公式了。比如lambda_D 1550而不是1550e-9算出来的 (\hat{\sigma}) 量级完全错误失配量可能大到 (10^{14}) 量级反射率自然趋近于零。如果是NaN或者Inf常见原因有两个一是扫描范围设置太大(\gamma dz) 的实部很大时(\cosh) 和 (\sinh) 数值溢出二是某个参数意外变成了0比如gamma0导致除法分母为零。其实第一种情况在正常参数范围内不太容易发生除非你从1300 nm扫到1800 nm还不收敛。解决办法是把扫描范围缩小到感兴趣的区域或者对矩阵元素做数值归一化处理。5.2 主峰位置不对或谱形不对称检查矩阵方向和边界符号主峰不在设计波长多半是布拉格波长算错或者Lambda没有根据lambda_B反算。检查一下lambda_D 2 * neff * Lambda是否等于你期望的中心波长。谱形不对称特别是整体看起来“歪了”十有八九是矩阵连乘方向写反或者边界条件里的负号丢了。可以回到2.3节的推导在纸上把矩阵方程写一遍和代码逐行对照。这里有一个快速自检方法把仿真得到的中心波长反射率和理论值 (\tanh^2(\kappa L)) 对比如果一致说明矩阵逻辑基本没问题如果差很多肯定是矩阵公式的问题。5.3 分段数M是不是越大越好很多初学者觉得数值仿真的分段数越大越精确这其实是个误区。传输矩阵法每一段内的传播和耦合是用解析公式描述的分段剖分只是为了逼近非均匀结构。对均匀FBG来说理论上M取1都能得到精确解M取500纯粹是为了后续扩展。实际操作中M取200到2000之间都是合理范围。M取太大的问题不是精度而是速度和数值累计误差。矩阵连乘次数多了舍入误差可能慢慢累积。如果做啁啾光栅仿真分段数应该根据光栅参数变化的剧烈程度来定变化快的地方可以加密变化慢的地方可以稀疏一点。均匀光栅仿真不用纠结这个问题。5.4 仿真跑得慢怎么办默认2001点、M500的情况下跑完只需要一两秒。如果你觉得慢大概率是扫描点设置太多或者把批量扫描嵌在了多重循环里。我的经验是先粗扫定位主峰位置再缩小范围细扫效率最高。比如先在1540到1560 nm用501个点粗扫找到主峰再在1549到1551 nm用2001个点细扫既快又准。如果要做大批量参数扫描MATLAB的parfor也可以省时间。把主程序里的for i 1:length(lambda_list)改成parfor i 1:length(lambda_list)并行计算立即生效。要注意子函数uniformFBG_TMM里不要依赖外部变量传参我写的这个函数是纯函数直接支持并行很方便。下面是一个常见的排查速查表你可以贴在工位旁边。症状可能原因排查方向反射率全为0波长单位没转成米检查lambda和Lambda出现NaN或Infgamma0或者扫描范围太宽打印gamma值缩小范围主峰位置偏了Lambda没按lambda_B反算核对lambda_D光谱左右不对称矩阵连乘顺序或符号问题对照2.3节手推一遍反射率恒等于1且极宽dn设得太大检查调制幅度通常小于1e-3RT不等于1代码里引入了损耗模型无损耗时直接用T1-R6. 实操扩展思路与个人经验这套代码我现在还经常用。前阵子做温度传感标定客户要求波长漂移分辨率到1 pm我先在MATLAB里把不同温度下反射谱扫了一遍确认峰值跟踪算法能稳定找到中心波长再上实验台用实际光源和光谱仪验证整个过程省了不少调试时间。仿真和实测当然会有差异比如实际光栅的写入误差、光纤双折射、光源谱宽等都会展宽光谱但趋势判断和参数预估MATLAB仿真完全足够。如果你想把代码能力再往前推一步我建议把子函数里的Lambda和dn从标量改成向量这样就实现了啁啾光栅和切趾光栅的仿真框架。均匀FBG只是所有分段参数恒定的特例改起来只需要把每段的矩阵用不同参数计算外层逻辑几乎不动。我自己当初就是从均匀光栅起步把代码改造成啁啾、切趾、相移通用版之后做光纤传感和通信滤波器的项目都顺手多了。最后再分享一个小技巧仿真脚本里最好保留一个“自校验开关”每跑完一组参数就输出中心波长反射率和理论值对比一旦偏差超过1e-6就报警提示。这个习惯帮我挡掉了好几次因为参数单位写错导致的低级失误。光纤世界里的奥秘往往就藏在那些看起来不起眼的数值细节里。