企业官网建设流程全解析

光纤光栅传感器现在用得越来越普遍，尤其是在结构健康监测、航空航天、土木工程和油气管道领域。它的核心优势是抗电磁干扰、耐腐蚀、适合分布式测量，而应变测量是最经典、最成熟的应用方向。做研究和工程应用时，第一步往往不是直接去拉光纤，而是先在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可以实现，但那是另一个更大的话题了。先用这一套静态脚本把物理机制吃透，再去折腾交互界面，路径更稳。