简介这套MATLAB仿真程序包面向光纤光栅领域的科研人员、光通信工程师以及相关专业学生覆盖均匀光纤光栅、啁啾光纤光栅和长周期光纤光栅三类常见器件重点解决光栅结构参数设计、反射与透射光谱计算、色散补偿以及传感特性分析等实际问题。压缩包共包含30个文件以.m格式的脚本文件为主另有.fig图形文件与.asv格式的自动备份文件其中.m脚本对应传输矩阵法、常微分方程数值求解、符号计算、等效镜像等多种建模实现.fig文件用于展示界面或仿真结果整体仅181KB轻量且便于直接运行和二次修改。程序支持输入光栅周期、长度、啁啾度、折射率调制深度等参数直观输出波长响应、反射谱、透射谱与时延曲线。同时考虑了温度、应变等环境因素对光栅性能的影响覆盖从基础理论验证到工程优化设计的完整环节。已有1291人学习下载适合光通信、光纤传感方向的研究者与工程师借鉴使用。 很多人拿到光纤光栅相关的项目第一反应是找现成的仿真程序。做光纤传感解调的人想预估反射谱做光纤激光器的想算透射窗口做切趾工艺的想对比旁瓣抑制效果——需求各不相同但落到工具层面都需要一套能快速输出光谱的仿真代码。我最早开始写光纤光栅仿真程序也是被实验逼的紫外直写一次光栅成本不算高但如果周期、切趾函数、长度全靠试来回耗费的时间和材料足够让人头疼。后来我干脆把均匀、啁啾、相移、切趾几种常用光栅的仿真统一成一套框架改参数就能出结果省下来的时间不是一点点。这篇文章就把这套思路完整展开从算法选型到代码实现再到实际跑谱时容易踩的坑都讲清楚。适合刚接触光纤光栅仿真、或者想把手头零散脚本整理成规范程序的读者就算你之前只听说过传输矩阵法照着下面的框架也能在半天内跑出自己的反射谱。1. 多种光纤光栅仿真到底在仿什么1.1 四种常见结构和它们的光谱特征光纤光栅本质上是纤芯折射率沿轴向的周期调制光在其中传播时满足相位匹配条件的波长会发生耦合反射或者透射。按折射率调制方式的不同最常见的几种结构如下光栅类型结构变化典型光谱特征典型应用均匀布拉格光栅UFBG周期和折射率调制幅度恒定单个窄反射峰旁瓣随耦合强度增大而明显光纤传感、窄带滤波啁啾光栅CFBG周期沿轴向逐渐变化反射带宽显著展宽不同波长对应不同反射位置色散补偿、宽带滤波相移光栅PS-FBG在光栅中间插入一个或多个相移点反射带内打开一个窄透射窗口DFB光纤激光器、窄带滤波切趾光栅Apodized FBG折射率调制幅度沿轴向按函数包络变化旁瓣被抑制主峰更干净高精度传感、波分复用这四种结构看似各不相同但仔细看会发现它们改变的其实只是三个东西每段的周期值、每段的折射率调制幅度、段与段之间是否插入相位跳变。换句话说只要底层算法支持这三个自由度的变化一套代码就能覆盖绝大多数光栅类型。这也是我坚持用传输矩阵法而不是其他方法的核心原因。1.2 为什么传输矩阵法是主力方法光纤光栅的严格分析要从耦合模方程出发。对于均匀光栅耦合模方程有解析解可以直接写出反射率公式但对于啁啾光栅、相移光栅、切趾光栅这类非均匀结构解析解基本不存在只能靠数值手段。常见的数值方案无非三种直接数值积分耦合模方程、有限元/有限时域仿真、传输矩阵法。直接积分耦合模方程概念上清晰但遇到相移这种离散突变时要单独处理代码写起来麻烦有限元仿真比如用 COMSOL精度高可一次计算要跑很久根本不适合做参数扫描传输矩阵法把光栅切成很多小段每段近似当作均匀光栅处理用 2x2 矩阵描述输入输出关系再把所有段的矩阵乘起来。它既能处理非均匀结构计算速度又快在个人电脑上跑几千个波长点也就几秒的事。用生活类比来说传输矩阵法就像把一段复杂的道路拆成很多个直线路段每一段用简单的运动学公式描述最后把结果串起来。只要每段足够短误差就可以压得很低。2. 仿真程序的整体设计从参数到矩阵2.1 你真正需要输入的几个参数写仿真程序前先明确哪些参数是必须的。我常用的参数表如下参数符号典型值说明有效折射率neff1.45单模光纤1550nm附近通常取1.447~1.45布拉格波长λB1550 nm目标反射峰位置光栅周期Λ约534 nm由 λB 2 × neff × Λ 反算折射率调制深度Δn1e-4紫外直写常见量级飞秒写制可能更高光栅长度L10 mm常见几毫米到几厘米切趾函数-高斯/余弦可选参数均匀光栅不用啁啾系数C0.1~1 nm/cm啁啾光栅需要相移量φπ相移光栅需要有个细节值得注意λB 和 neff 给定之后周期 Λ 就不是独立变量了必须用 λB 2 neff Λ 去算而不是随便填一个。很多人第一次跑仿真反射峰不在预期位置基本都是这个原因。2.2 传输矩阵的构造与物理意义把长度为 L 的光栅均匀分成 N 段每段长度 h L/N。对第 i 段定义两个关键量交流耦合系数 κ π × Δn × v / λ其中 v 是条纹可见度通常取 1。直流自耦合系数 σ̂ 2 π neff / λ - π / ΛiΛi 是这一段的周期。如果存在平均折射率变化σ̂ 里还要加一项 2π × δn_dc / λ但很多仿真场景可以忽略。再定义 γ² κ² - σ̂²。第 i 段的传输矩阵是S_i [ cosh(γh) - i(σ̂/γ)sinh(γh) -i(κ/γ)sinh(γh) i(κ/γ)sinh(γh) cosh(γh) i(σ̂/γ)sinh(γh) ]光从一端进入经过所有小段后总的传输矩阵 S_total S_N × ... × S_2 × S_1。假设光从左边入射右边输出端没有后向入射波那么反射系数 r S_total[1,0] / S_total[0,0]反射率 R |r|²。这里要注意不同教材对矩阵元素和符号的约定可能略有差异有的把虚部符号写成相反方向有的把反射系数定义成别的分量。你只要固定使用同一套约定最后算出来的反射率 R 是一致的。我第一次对照文献核结果时就被符号问题绕了一下后来干脆只信自己的最终光谱输出。2.3 波长扫描范围和数值精度仿真程序本质上是在每个波长点上算一次传输矩阵波长必须离散化。扫描范围取决于光谱特征均匀布拉格光栅的3dB带宽通常在0.1~0.3 nm量级所以扫描范围至少要有 ±5 nm波长点数至少1000个以上否则谱线会变得很难看。啁啾光栅的情况不一样反射带宽可能达到几纳米甚至十几纳米扫描范围必须覆盖整个反射带。经验做法是先按理论公式估算带宽然后把扫描范围设成估算值的3~5倍跑完看谱图再调整。计算量方面N 取 500、波长点取 2000 时总共要算 100 万次小矩阵乘法纯 Python 循环大概需要几十秒。如果觉得慢可以用 NumPy 一次性向量化批量波长点性能能提升一个数量级。3. 一套代码实现多种光纤光栅3.1 均匀布拉格光栅先跑通最基础的反射谱我习惯用 Python 写因为 NumPy 和 Matplotlib 生态太方便了。下面这个函数是均匀布拉格光栅的核心实现去掉注释后不到30行import numpy as np import matplotlib.pyplot as plt def simulate_ufbg(neff1.45, L0.01, dN1e-4, lamB1550e-9, num_seg500, num_wavelength2000): # 由布拉格波长反算周期 period lamB / (2 * neff) dz L / num_seg # 波长扫描范围中心波长附近 ±5 nm lambdas np.linspace(lamB - 5e-9, lamB 5e-9, num_wavelength) reflectivity [] for lam in lambdas: # 每个波长点重新计算传播常数、失谐量、耦合系数 k0 2 * np.pi / lam sigma k0 * neff - np.pi / period # 直流自耦合项 kappa np.pi * dN / lam # 交流耦合系数 gamma np.sqrt(kappa**2 - sigma**2 0j) # 逐段乘传输矩阵 T np.eye(2, dtypecomplex) for _ in range(num_seg): S np.array([ [np.cosh(gamma * dz) - 1j * sigma / gamma * np.sinh(gamma * dz), -1j * kappa / gamma * np.sinh(gamma * dz)], [1j * kappa / gamma * np.sinh(gamma * dz), np.cosh(gamma * dz) 1j * sigma / gamma * np.sinh(gamma * dz)] ]) T T S r T[1, 0] / T[0, 0] reflectivity.append(abs(r)**2) return lambdas, np.array(reflectivity) lambdas, R simulate_ufbg() plt.plot((lambdas - 1550e-9) * 1e9, R) # 横轴转换成偏离1550nm的纳米数 plt.xlabel(Wavelength offset (nm)) plt.ylabel(Reflectivity) plt.show()用表格里的典型参数跑出来峰值反射率大约在 90% 以上3dB 带宽在 0.1 nm 量级旁瓣出现在主峰两侧。这个反射率对紫外直写光栅来说是合理的因为 κL 已经接近 2属于强耦合状态。如果你想要更低的反射率可以把 Δn 降到 5e-5或者把 L 缩短到 5 mm。跑通这个函数之后你就有了一个基准。后续所有变体都只需要在这个框架上做局部修改。3.2 啁啾与相移光栅只改周期和相位矩阵啁啾光栅的改动非常小每一段的周期不再恒定而是沿长度方向线性变化。假设中心布拉格波长为 λB0啁啾系数为 C单位用 nm/m那么第 i 段对应的布拉格波长是lamB_i lamB0 C * (i * dz - L / 2)然后每段用 lamB_i 重新计算周期 period_i lamB_i / (2 * neff)。代码里只需要把原来固定 period 的地方改成每段计算C 2e-9 # 2 nm/m即0.02 nm/cm for i in range(num_seg): z_i i * dz lamB_i lamB0 C * (z_i - L / 2) period_i lamB_i / (2 * neff) sigma k0 * neff - np.pi / period_i # ... 组装 S 矩阵注意啁啾系数的正负号决定长波长端在光栅的哪一侧实验上对应光从哪端入射。如果发现反射谱的群时延斜率方向不对多半就是符号选反了。相移光栅的改动更特殊不是在每段参数上做文章而是在段与段之间插入一个相移矩阵。比如在光栅正中间插入 π/2 相移phase np.pi / 2 P np.array([ [np.exp(-1j * phase / 2), 0], [0, np.exp(1j * phase / 2)] ]) T np.eye(2, dtypecomplex) mid num_seg // 2 for i in range(num_seg): # ... 计算当前段的 S 矩阵 T T S if i mid: T T P # 在中间插入相移相移矩阵的作用等价于在光栅中间人为引入一个光程突变。π 相移会在反射带中心打开一个非常窄的透射峰这个透射峰就是 DFB 光纤激光器的选模窗口。仿真时把透射率 T 1 - R 画出来就能看到这个窗口的宽度和深度。3.3 切趾光栅与长周期光栅的扩展思考切趾光栅处理起来也很简单让每段的 κ 乘以一个包络函数。最常用的是高斯切趾def gaussian_apod(z, L, strength8.0): # strength 越大边缘压制越狠 return np.exp(-strength * ((z - L / 2) / L) ** 2) for i in range(num_seg): z_i i * dz kappa_i (np.pi * dN / lam) * gaussian_apod(z_i, L) # ... 组装 S 矩阵时用 kappa_i切趾的效果非常直观均匀光栅反射谱两侧的旁瓣会被大幅压制主峰变得干净。代价是有效耦合长度变短峰值反射率会下降一点主瓣宽度略微变宽。工程上需要在这三个指标之间做权衡。长周期光栅LPG严格来说需要处理纤芯模和包层模之间的同向耦合传输矩阵的形式会变成 4x4 或者需要用两组耦合常数分别计算不能直接套用上面的 2x2 矩阵。不过如果你只需要粗略估算透射谱的包络也可以把包层模等效成一个有效折射率再用类似的矩阵框架去推。我自己在项目里做粗设计时用过这个近似趋势是对的但要精确定标还是有差距。4. 实测中那些坑排查与技巧4.1 分段数 N 和波长点数选多少这两个参数直接决定结果可信度。分段数 N 太小每段长度 h 相对周期来说太大“分段均匀”的近似就不成立光谱上会出现类似法布里-珀罗腔的纹波看起来像真实旁瓣其实是数值误差。经验值是每个周期至少要切 5~10 段。周期约 534 nmL10 mm 时光栅内有约 18700 个周期所以 N 取 500 时每段包含约 37 个周期是足够的N 取 100 时每段约 187 个周期也还能接受但边缘的啁啾细节会丢失。波长点数同样关键。带宽只有 0.1 nm 时如果你只扫 200 个点分辨率就是 0.05 nm几乎画不出谱线形状。我一般保证带宽范围内至少有 200 个波长点也就是总点数至少 2000。4.2 谱线不对称、多出毛刺的排查写程序时最容易出问题的几个点我按踩坑频率排个序第一σ̂ 的符号搞反。σ̂ 2π neff / λ - π / Λ 这个顺序不能乱反了之后谱线会左右反转看起来像镜像峰值波长还会偏移。第二啁啾光栅的周期更新位置不对。每段的周期必须基于该段的实际位置计算不能在循环外只算一次。否则你写的“啁啾”实际还是均匀光栅。第三相移矩阵的插入顺序。矩阵乘法不满足交换律T T P S 和 T T S P 结果完全不同。相移矩阵必须插在对应的两个段矩阵之间不能图省事统一在最后乘。第四强光栅下 sinh 函数溢出。κL 很大时γ h 的值可能很大sinh 函数会膨胀到天文数字矩阵元素数值溢出。解决方法是提高数值精度用 float128或者对强光栅改用超导近似公式。我自己的经验是反射率超过 99% 的强光栅用 float64 就可能出问题。4.3 从仿真到实验的衔接心得仿真跑得再漂亮最终都要和实验对得上。我的体会是仿真程序更适合做趋势预测和参数敏感性分析而不是指望它精确复现某一次实验光谱。实际光纤光栅和理想模型之间总有偏差比如紫外写制时折射率调制沿长度方向不均匀、切趾函数实际形状和理论有出入、光纤本身的双折射会导致偏振相关峰分裂、neff 随波长变化会产生二阶色散等等。把这些因素全考虑进去模型就太复杂了没那个必要。比较务实的做法是用仿真确定关键参数范围比如周期调多少能让反射峰落在目标波长、啁啾系数取多少能满足带宽要求、切趾强度选多大能把旁瓣压到 -30 dB 以下。确定之后再去写光栅实验数据反过来再校准仿真参数形成一个迭代闭环。我手头很多光栅的 Δn 实际值和标称值差 20% 以上都是这么校准出来的。5. 一点个人建议如果你准备自己写光纤光栅仿真程序我建议按照“均匀光栅 → 切趾 → 啁啾 → 相移”的顺序逐步扩展每加一种功能就画一次光谱和理论值对比。不要一上来就想写一个万能程序步子迈大了后面排查问题会非常痛苦。另外代码里尽量把参数和算法分开参数放配置文件或者函数入参里算法只负责算矩阵。这样以后换项目、改参数不需要动核心代码。我现在做新项目时基本就是复制这套框架改几个参数就能跑效率高很多。希望这篇文章能帮你少走一些弯路。本文还有配套的精品资源点击获取