光纤光栅传感器现在用得越来越普遍尤其是在结构健康监测、航空航天、土木工程和油气管道领域。它的核心优势是抗电磁干扰、耐腐蚀、适合分布式测量而应变测量是最经典、最成熟的应用方向。做研究和工程应用时第一步往往不是直接去拉光纤而是先在MATLAB里把反射谱仿真做出来再决定光栅参数怎么设计、解调算法怎么匹配。这篇内容我围绕“MATLAB实现均匀应变和非均匀应变光纤光栅仿真”展开会把两种应变场景下的建模思路、公式推导、代码实现和踩坑记录都写清楚全程基于我自己调试过、验证过的方案。无论你是刚接触光纤光栅的硕士生还是在做传感器设计的工程师这套仿真流程都能直接拿过去用。1. 为什么用MATLAB做光纤光栅应变仿真1.1 这个仿真到底在解决什么问题光纤光栅FBG本质上是一段纤芯折射率周期性调制的结构。当入射光进入光栅后满足布拉格条件的波长会被反射回来。用公式表达就是布拉格条件λ_B 2 × n_eff × Λ其中n_eff是有效折射率Λ是光栅周期。任何改变这两个物理量的外部因素都会导致布拉格波长移动。应变的作用就是通过弹光效应改变折射率同时直接拉伸或压缩光栅周期从而实现波长编码的传感。所以仿真首先要回答的问题是在一个给定的应变量下反射谱会怎么移动、怎么变形。如果光栅受到的是均匀轴向应变那整个光栅的周期和折射率变化是一致的反射谱只发生整体平移形状基本不变。这个场景用解析公式就能算。但实际情况往往没那么理想。比如应变施加在结构表面时由于胶水厚度不均匀、粘贴区域局部受力、结构本身存在应力集中区光栅沿线感受到的应变很可能是不均匀的。此时光栅被分段拉伸不同段对应的布拉格波长不同反射谱不再是一个规整的尖峰而会出现展宽、分裂、多峰等现象。这个仿真的价值就在于第一可以在不花钱买光纤器件的情况下先把各种光栅参数和应变分布对光谱的影响摸清楚第二给解调算法提供理论依据知道在非均匀应变下光谱会呈现什么特征才能设计合适的中心波长提取方法。1.2 均匀与非均匀应变的核心差异均匀应变情况下整个光栅的等效布拉格波长是标量整个光栅如同一个新周期、新折射率的均匀光栅仿真结果是一条平滑的高斯状反射峰应变大小直接对应中心波长偏移量。这种情况下用公式Δλ_B λ_B × (1 - p_e) × ε就能快速算出应变灵敏度系数。对于1550nm波段的光纤光栅p_e约为0.22因此理论上灵敏度系数约为1.2 pm/με。这意味着1微应变产生大约1.2皮米的波长偏移。非均匀应变情况下光栅不同位置的有效布拉格波长不同整个结构变成了一个啁啾光栅或分段啁啾光栅。光谱的形状就是这些不同局部反射叠加干涉的结果。此时解析法虽然也能勉强处理某些特殊分布比如线性啁啾但通用性很差。一旦应变分布是任意的、分段跳变的就必须用数值方法最常用的是传递矩阵法。2. 均匀应变FBG解析法的MATLAB实现2.1 理论铺垫耦合模方程与反射率公式均匀光纤光栅的反射率可以从耦合模方程推导出来最终得到的是包含双曲函数的解析表达式。设光栅长度L耦合系数κ失谐量Δβ则定义直流自耦合系数σ的实部为失谐量交流耦合系数κ表示折射率调制的反射强度。在忽略损耗的理想情况下反射率公式为r -iκ × sinh(γL) / (γ × cosh(γL) i × Δβ × sinh(γL))其中γ² κ² - Δβ²失谐量Δβ 2π × n_eff × (1/λ - 1/λ_B)。反射率功率为R |r|²。耦合系数κ由折射率调制深度δn决定。实际光栅写入时如果折射率调制为δn_eff则交流耦合系数近似为κ π × δn_eff / λ需要注意的是不同文献中δn_eff的定义可能略有差异。有的把折射率调制幅度直接叫δn有的用Δn指的是折射率变化峰值的一半。我在仿真中使用的是“有效折射率调制深度”这个量含义是折射率余弦调制的振幅值。如果要在MATLAB里定义这个参数建议先确认你参考的文献的定义否则算出来的峰值反射率会有偏差。均匀应变建模的核心就是把应变映射到λ_B上。应变ε作用下光栅周期变为Λ(1ε)有效折射率由于弹光效应变为n_eff(1 - p_e × ε)。因此应变量ε对应的布拉格波长为λ_B(ε) 2 × n_eff × (1 - p_e × ε) × Λ(1ε)忽略高阶小量后就是Δλ_B λ_B × (1 - p_e) × ε。在仿真代码里直接更新λ_B代入反射率公式即可。2.2 MATLAB代码实战从参数定义到反射谱绘制这套代码我在实际项目里跑过很多次步骤很简单关键是参数的单位和物理量定义要一致。下面给出一段完整的MATLAB脚本示例实现不同均匀应变下反射谱的计算和绘制。% 均匀应变FBG反射谱仿真 % 思路更新应变后的布拉格波长代入耦合模解析反射率公式 clear; close all; clc; % 基础参数 n_eff 1.468; % 有效折射率 Lambda 528.5e-9; % 光栅周期单位米 L 10e-3; % 光栅长度单位米10mm d_n 2.5e-4; % 折射率调制幅度交流耦合用 p_e 0.22; % 有效弹光系数 lambda_B0 2 * n_eff * Lambda; % 初始布拉格波长 % 波长扫描范围以布拉格波长为中心 lambda_range linspace(lambda_B0 - 1.5e-9, lambda_B0 1.5e-9, 5000); % 应变列表单位微应变 strain_list [0, 250, 500, 1000]; % 0、250、500、1000 με figure; hold on; grid on; box on; for i 1:length(strain_list) eps_applied strain_list(i) * 1e-6; % 应变转为无量纲 % 应变后的布拉格波长 lambda_B lambda_B0 * (1 (1 - p_e) * eps_applied); % 计算各波长下的耦合系数与失谐量 kappa pi * d_n ./ lambda_range; delta_beta 2 * pi * n_eff .* (1 ./ lambda_range - 1 / lambda_B); gamma sqrt(kappa.^2 - delta_beta.^2); % 反射系数计算 r_num -1i .* kappa .* sinh(gamma .* L); r_den gamma .* cosh(gamma .* L) 1i .* delta_beta .* sinh(gamma .* L); r r_num ./ r_den; R abs(r).^2; % 归一化波长偏移便于同一坐标下比较 plot((lambda_range - lambda_B0) * 1e9, R, LineWidth, 1.5, ... DisplayName, sprintf(%d με, strain_list(i))); end xlabel(相对波长偏移 (nm)); ylabel(反射率); title(均匀应变FBG反射谱); legend(show);这段脚本的核心逻辑就三件事定义静态参数、把应变折合成布拉格波长、再挨个波长算反射率。运行速度很快5000个波长点加4条光谱在老一点的笔记本上也就是一两秒的事。我在调试时遇到过一个问题当kappa大于delta_beta时gamma是实数反射谱呈现主峰加多个旁瓣结构这是正常的。旁瓣的高度和数量与kappa×L的乘积有关这个乘积越大主峰反射率越高但旁瓣也越明显。如果写光栅时过曝折射率调制过大旁瓣会变得非常高传感解调时候会造成干扰。2.3 均匀应变下的波长漂移量检验跑完上面的代码后你会发现一个直观的规律均匀应变下反射谱的波形几乎不变只是整体往长波长方向平移。以1550nm左右的光栅为例1000微应变大约对应1.2nm的偏移。这个偏移量非常大在光谱仪上非常明显。为了精确读取中心波长通常用峰值搜索法直接找反射率最大值所在的波长。更精细的做法是高斯拟合或质心法。在均匀应变下峰值搜索法和质心法结果一致因为光谱是对称的。这个一致性本身就是判断应变是否均匀的一个重要指标。所以我建议你仿真时顺手把中心波长提取写进去比如用findpeaks函数[peak_val, peak_idx] max(R); lambda_peak lambda_range(peak_idx); fprintf(应变 %d με 时峰值波长偏移: %.4f nm\n, ... strain_list(i), (lambda_peak - lambda_B0) * 1e9);实测下来计算得到的灵敏度大约在1.203 pm/με附近与理论值1.2 pm/με对得上。这种一致性说明代码的物理建模是可靠的。3. 非均匀应变FBG传递矩阵法实现3.1 为什么均匀模型在非均匀应变下失效非均匀应变下光栅周期Λ和有效折射率n_eff沿长度方向都是位置z的函数。布拉格波长λ_B(z) 2×n_eff(z)×Λ(z)也跟着变化每个位置的局部反射波长不同。如果强行用均匀反射率公式只能算出一条整体移位的尖峰无法捕捉光谱的展宽和多峰结构。从物理图像上看非均匀应变下的光栅相当于很多个子光栅串联起来每个子光栅的反射波长略有差异。光在传播过程中每个界面都会产生反射这些反射波会相互干涉最终叠加出一个复杂的反射谱。这个“分段叠加”的思想正是传递矩阵法的物理基础。传递矩阵法把光栅分成M段每一段足够短可以近似看成均匀FBG。对每一小段用耦合模解析解构造一个2×2的传输矩阵然后把所有段的矩阵按顺序乘起来得到整个光栅的传输矩阵。最后从总矩阵中提取反射系数。注意这里假设每段内的应变是常数这要求分段数要足够多。分段数太少每段内的应变梯度被粗暴平均掉光谱细节会失真分段数太多计算矩阵连乘的时间会线性增加。实践下来10mm光栅分100段就是比较稳的平衡点如果追求精细结构可以分200段。3.2 传递矩阵法的基本思想与分段建模传递矩阵法的公式推导很多教材都有我这里只写关键步骤和MATLAB实现中容易出错的地方。每段均匀FBG的传输矩阵形式为M_i [ cosh(γ_i dz) - i(Δβ_i/γ_i)sinh(γ_i dz), -i(κ_i/γ_i)sinh(γ_i dz) ; i(κ_i/γ_i)sinh(γ_i dz), cosh(γ_i dz) i(Δβ_i/γ_i)sinh(γ_i dz) ]其中dz是分段长度Δβ_i是第i段对应的失谐量κ_i是第i段的耦合系数。失谐量的计算需要用到该段的局部布拉格波长λ_B,i。应变量ε(z)作用在第i段时λ_B,i λ_B0 × (1 (1 - p_e) × ε_i)其中ε_i是第i段中点位置的应变值。这样就把应变分布转换成了每段的布拉格波长序列。整个光栅的总传输矩阵为M_total M_M × M_{M-1} × ... × M_2 × M_1反射系数通过总矩阵的第二行第一列与第一行第一列的比值得到r -M_total(2,1) / M_total(1,1)这里需要特别注意符号约定。我开始自己写代码时反射谱没有画出来反射率一直是随着波长平缓变化而没有峰折腾了很久发现是符号问题。不同的文献对传输矩阵内场的正方向定义不同有的定义前向波在第二行有的在第一行导致M21/M11的符号和正负号差异。如果你用这个公式画出来反射率1那大概率是边界条件取反了。功率谱R|r|²不涉及符号但如果后续你要做相位分析这个约定必须统一。3.3 核心MATLAB实现分段FBG矩阵连乘下面这段代码是三段式可以自由修改应变分布函数。我给了三个预置场景均匀、线性梯度、分段跳变。你自己用的时候只需要改strain_profile这个函数的返回值就可以了。% 非均匀应变FBG反射谱仿真传递矩阵法 clear; close all; clc; % 基础参数 n_eff 1.468; Lambda 528.5e-9; L 10e-3; M 200; % 分段数 d_n 2.2e-4; p_e 0.22; lambda_B0 2 * n_eff * Lambda; % 波长扫描范围 lambda_range linspace(lambda_B0 - 2e-9, lambda_B0 2e-9, 3000); % 定义应变分布微应变随位置z变化z从0到L z_pos linspace(0, L, M); strain_profile 400 * (z_pos / L); % 线性非均匀0到400 με % 将应变转化为各段的局部布拉格波长 eps_val strain_profile * 1e-6; lambda_B_local lambda_B0 * (1 (1 - p_e) * eps_val); % 逐波长计算反射率 R_spectrum zeros(size(lambda_range)); dz L / M; for idx 1:length(lambda_range) lambda lambda_range(idx); % 总传输矩阵初始化为单位矩阵 M_total eye(2); for i 1:M % 本段失谐量 delta_beta 2 * pi * n_eff * (1/lambda - 1/lambda_B_local(i)); kappa pi * d_n / lambda; gamma sqrt(kappa^2 - delta_beta^2); % 构造单段传输矩阵 cg cosh(gamma * dz); sg sinh(gamma * dz); Mi [cg - 1i*delta_beta/gamma*sg, -1i*kappa/gamma*sg; 1i*kappa/gamma*sg, cg 1i*delta_beta/gamma*sg]; % 矩阵连乘注意顺序 M_total Mi * M_total; end % 反射系数 r -M_total(2,1) / M_total(1,1); R_spectrum(idx) abs(r)^2; end % 绘制反射谱 figure; plot((lambda_range - lambda_B0)*1e9, R_spectrum, LineWidth, 1.5); xlabel(相对波长偏移 (nm)); ylabel(反射率); title(非均匀应变FBG反射谱线性梯度应变); grid on;这段代码的结构是“外层扫波长内层扫分段”。对每个波长点都要把200个段的矩阵乘起来所以计算量是均匀模型的好几百倍。3000个波长点叠加200段在老机器上可能要跑十几秒。如果你想加速可以先用粗波长间隔扫一遍把光谱的大致轮廓摸清楚再在兴趣峰附近用细间隔加密。另外矩阵元素如果有接近无穷大的情况可以检查是不是gamma*dz过大一般不会出现因为dz很小。4. 仿真结果解读与典型场景分析4.1 线性梯度应变光谱展宽与啁啾现象线性梯度应变是最常见的非均匀应变场景之一比如悬臂梁表面粘贴光栅时应变沿光栅长度方向呈线性变化。仿真的结果中反射谱会明显展宽主峰变钝有时候甚至出现多个相邻的小峰。我们可以从光学原理上理解线性梯度应变相当于把光栅变成了线性啁啾光栅。光栅各段的布拉格波长从低到高依次排列整个结构像一个波长选择反射器不同波长的光在不同位置被反射回来。因此反射谱的带宽由应变梯度决定梯度越大展宽越明显。在做传感器设计时如果你发现实测光谱展宽了说明光栅粘贴区域的应变确实不均匀。这时候中心波长的定义就要谨慎了。用3dB带宽中点定义中心波长会比峰值搜索更稳定但物理含义有所变化。我会建议使用质心法这个指标在光谱展宽情况下变化相对平滑。4.2 局部应变区旁瓣特征与定位分析另一个典型的非均匀场景是局部应变例如光纤光栅中间某一段受到了额外应力而其他区域应变为零。此时反射谱会在主峰旁边出现一个或多个额外的峰这些峰的位置对应局部应变区域的布拉格波长。这段仿真代码只需要改动应变分布函数即可strain_profile zeros(M, 1); strain_profile(80:120) 600; % 中间约2mm区域受到600με应变跑出来的反射谱如果出现一个明显的次级峰就可以根据该峰与主峰的波长间隔估算局部应变的大小。这个特性被用在分布式光纤传感中通过分析光谱形状来判断应变发生的位置和大小。不过要提醒一句局部应变区域非常短比如小于光栅长度的1/10时次级峰幅度会很低可能淹没在旁瓣或噪声里。要提高检测灵敏度可以增加光栅折射率调制深度或者用更短的光栅。这是传感方案设计时要权衡的点。4.3 应变传感设计的参数选择建议综合以上仿真经验我给出几个参数选择建议这些都是在实际项目里调试过程中总结出来的比理论推导更贴近工程需求。第一光栅长度。光栅越长反射谱越窄峰值反射率越高波长解调精度越好。但光栅越长对非均匀应变的敏感性也越高光谱容易展宽。所以测量大梯度应变时建议用短光栅比如2mm到5mm测量微小的均匀应变时可以用长光栅提高精度。第二折射率调制深度。增大d_n可以提高反射率但旁瓣也会明显增大。如果解调算法比较弱建议用切趾技术这在光栅写入时就要设计。仿真里可以加一个切趾函数比如高斯窗对κ进行调制代码如下apod exp(-0.5 * ((z_pos - L/2) / (L/6)) .^ 2); kappa_i pi * d_n ./ lambda_range(i) * apod(i);加了切趾之后旁瓣能压掉很多代价是峰值反射率略微下降带宽稍微变宽。第三波长扫描范围与分辨率。反射谱宽度和应变范围成正比建议仿真时扫描范围设为预期应变对应的波长跨度再加20%余量。波长点数至少2000否则展宽谱的形状不够平滑后续解调算法的精度验证也没有意义。5. 常见缺陷与排错记录5.1 反射率大于1转移矩阵方向与符号约定这个问题我在网上帮别人看代码时遇到过很多次。反射率大于1显然不符合物理事实根因几乎都出在传输矩阵的边界条件上。一种情况是总矩阵的连乘顺序反了我代码里写的是M_total Mi * M_total如果写成M_total M_total * Mi矩阵乘法不满足交换律结果就乱了。另一种情况是反射系数公式的符号约定不对乘以负号或者取倒数位置。排查方法很简单在均匀应变极限下也就是所有ε_i相同时用传递矩阵法算出来的反射谱应该和均匀解析法的结果完全重合。如果两条曲线对不上就是矩阵实现有问题。我建议你写代码时始终保留均匀解析法作为基准很多逻辑错误都能通过对比暴露出来。5.2 光谱锯齿状分段数与波长分辨率的关系如果把反射谱放大看出现密密麻麻的小锯齿多数是分段数不够或者波长点数不够。分段数不够时光栅被离散成粗格子每个格子就是一个均匀段段与段之间的边界会引入周期性的反射干涉形成锯齿。波长点数不够时扫出来的曲线本身就不平滑。经验值是10mm光栅分100到200段波长点数不低于2000反射谱基本是光滑的。如果仍然有锯齿检查应变分布是否连续分段跳变本身就会造成干涉条纹这是物理现象而不是数值问题。5.3 波长漂移不准单位换算与有效折射率修正这是一个特别容易踩的坑。应变值在工程中习惯用微应变με表示但在物理公式里要换算成无量纲的ε除以10^6。如果忘记换算算出来的波长漂移会大1000倍明显偏离理论值。另外光弹系数p_e是有效值它综合了光纤材料的弹光张量、泊松比等因素。不同文献给的值略有不同从0.204到0.24都有这和光纤掺杂类型有关。如果你仿真结果与标称灵敏度1.2 pm/με对不上先检查p_e取值再检查单位换算。这两个点都排除后还是不匹配那就要看n_eff和Λ的取值是否一致了因为λ_B0 2×n_eff×Λ必须落在你实际使用的波段。比如1550nm波段的光栅如果n_eff取1.468Λ应该在528.5nm左右这是同一个物理事实的两种表达。最后再分享一个小技巧转移矩阵法的代码结构非常适合做参数扫描比如你要分析不同应变梯度下的光谱演变只需要把外层加一个for循环把应变梯度作为循环变量最后把所有反射谱叠在一张图上就能看到光谱展宽的连续变化过程。这种可视化对写论文插图非常有帮助也能帮你快速判断设计参数的敏感区间。我实际调试时还发现一个事非均匀应变仿真如果加了下拉菜单或者交互式滑块用来动态调整应变分布函数整个分析过程会高效得多。MATLAB的App Designer可以实现但那是另一个更大的话题了。先用这一套静态脚本把物理机制吃透再去折腾交互界面路径更稳。