简介光纤光栅是光纤通信与传感领域的核心器件其反射谱和群时延特性直接决定了色散补偿、脉冲整形等应用的性能。对于周期渐变的啁啾光栅传统解析方法难以处理而传输矩阵法通过将光栅分段近似为均匀光栅并连乘传输矩阵能高效仿真宽带反射谱与线性群时延。Matlab凭借强大的矩阵运算和信号处理函数成为实现该算法的理想工具。本文从耦合模理论出发介绍传输矩阵法原理、关键参数啁啾量、光栅长度、折射率调制深度的影响规律并结合可运行代码展示如何提取反射谱相位、计算群时延以及规避相位解包裹、单位换算等常见陷阱。掌握这一仿真框架可进一步扩展至切趾、相移光栅等复杂结构为光栅器件的自主设计与工程调试提供扎实基础。 前几天有位朋友在交流群里丢了个压缩包文件名写着“zhoujiu.zip_matlab光栅时延_matlab啁啾光栅的反射谱、时延”。群友都心知肚明这类资源多半是从某论坛或淘宝店转了几手的东西代码能跑但注释乱、公式缺、参数写死真正想拿来用在课题里反而要花大量时间反推。与其研究别人的古董代码不如把啁啾光纤光栅仿真这件事从原理到落地完整捋一遍。这篇内容我尽量讲透为什么啁啾光栅的反射谱是宽带、群时延为什么是线性以及如何用matlab自己搭一套传输矩阵法仿真。适合光纤通信、光纤传感方向的学生以及刚接触光栅器件仿真、想快速上手又不想只做“调参侠”的工程师。1. 啁啾光纤光栅入门周期在变的“彩虹镜子”1.1 从均匀光栅到啁啾光栅周期在变先回到光纤光栅本身。所谓光纤光栅就是利用紫外光在光纤纤芯里写下周期性的折射率调制形成一层“波长选择反射镜”。均匀光栅的周期是恒定的根据布拉格条件 λ_B 2·n_eff·Λ只有满足这个波长条件的光才会被强烈反射其他波长直接透过。所以均匀光栅的反射谱是一个极窄的尖峰典型半高全宽也就0.1~0.3 nm量级非常适合做窄带滤波和单波长激光器选模。啁啾光栅的区别在于周期不再恒定而是沿光栅长度方向线性变化最常见的是线性啁啾。也就是说光栅前端的局部布拉格波长是λ_B1后端是λ_B2中间连续过渡。这就像把一面单色镜子换成了一排按彩虹顺序排列的小镜子每个位置只反射对应的一小段波长。于是整体反射谱被展宽成一个宽带而且不同波长的光在不同位置反射天然带来了额外的时延。用生活化一点的方式来理解均匀光栅像百米赛跑的直线跑道所有运动员同时出发、同时冲线啁啾光栅则像一条从山脚蜿蜒到山顶的盘山道不同速度的运动员被分到不同路段跑完所用的时间自然不一样。这就是啁啾光栅时延特性最直观的来源。1.2 啁啾光栅的核心应用色散补偿与脉冲整形反射谱只是回答了“哪些波长被反射”时延则回答“这些波长分别延迟了多少”。两者合在一起才构成啁啾光栅最核心的两个指标。为什么这么看重时延因为啁啾光栅最重要的应用之一就是色散补偿。光纤通信里光脉冲在普通单模光纤中传播时不同波长分量因为群速度不同到达接收端的时间有差异脉冲就被“拉宽”了。这在高码率长距离传输里是致命的。啁啾光栅的做法是把不同波长的分量安排到光栅的不同深度反射短波长分量先返回、长波长分量后返回或者反过来人为制造一个与光纤相反的光程差从而把被拉开的脉冲重新“压”回去。除了色散补偿啁啾光栅还广泛用于超短脉冲整形、宽带滤波、多波长传感等领域。我见过不少做光纤传感的朋友把啁啾光栅当宽带反射器用配合可调谐激光器做解调也是常规操作。所以无论从通信还是传感角度反射谱和时延这两个仿真结果都是后续一切工作的基础。1.3 为什么选择matlab作为仿真工具matlab在这个领域的地位还是挺稳的。矩阵运算天然适合传输矩阵法内置的unwrap、gradient、plot函数让数据处理和可视化非常顺手而且网上现成资料多。但问题也出在“现成”上很多流传的代码要么参数单位混乱要么绕过了关键的物理项跑出来谱型怪怪的新手又不知道问题出在哪。我个人的建议是不要把“能跑”当作目标而是把底层原理弄明白自己搭一套干净的框架。这套框架以后加切趾、加相移、换成取样光栅都只是改几行参数的事。下面我从耦合模理论讲起把传输矩阵法的来龙去脉理清楚。2. 仿真方案选型为什么用传输矩阵法2.1 耦合模理论与传输矩阵法的关系光栅内部的电磁场分析经典路径是用耦合模理论。前向传播模式和反向传播模式之间因为折射率微扰而相互耦合最终得到一组耦合模方程。对于折射率调制均匀的理想光栅这组方程有严格的解析解解出来就是反射谱的表达式。但加上啁啾之后周期沿长度变化局部布拉格波长不再是常数严格解析解就变得很麻烦。工程上最通用的替代方案就是传输矩阵法Transfer Matrix MethodTMM。传输矩阵法的思路可以用一句话概括把整个光栅切成很多小段每一段足够短近似看成均匀光栅然后每段用解析解写成一个2×2矩阵最后把所有矩阵连乘得到整个光栅的传输特性。这就像用很多条直线段去逼近一条曲线切得越细、越接近真实情况只是计算量会上升。对20 mm长的光栅切成1000段每段只有20 µm局部的周期变化完全可以忽略不计每一段当作均匀光栅处理是完全合理的。具体到每一段的2×2矩阵核心物理量是三个失谐量σ表示入射波长偏离该段局部布拉格波长的程度交流耦合系数κ表示折射率调制幅度带来的前后向模式耦合强度段长Lseg即每段的物理长度。有了这三个量就可以构造出标准的均匀光栅传输矩阵再在波长循环里逐段连乘。虽然代码写起来只是几行矩阵乘法但背后是耦合模方程在分段区间上的解析解表达这就是整个方法的数学根基。2.2 时延计算的数学基础反射谱很容易得到反射系数r是一个复数反射率R |r|²。时延则要从r的相位里提取。设反射系数的相位为φ群时延的定义是τ dφ/dω。因为仿真里用的自变量是波长而不是角频率需要做一个变量替换dω -2πc/λ² · dλ换算过来就是τ -λ²/(2πc) · dφ/dλ这个负号经常把人绕晕。我自己的记法是波长增大时角频率减小所以用λ求导时要加负号把方向转回来。数值上dφ/dλ可以用一行gradient(phi, lambda)得到非常方便。这里有一个新手必踩的坑angle函数返回的相位是包裹在[-π, π]区间的直接求导会在跳变处产生巨大的尖峰所以必须先做unwrap。即使如此低反射区域的相位本身噪声很大unwrap之后也可能出现不合理的跳变后面第5章我会讲怎么处理。2.3 关键参数设定与单位处理仿真参数的单位问题是最容易翻车的地方。matlab没有量纲检查nm、m、ps混着写也不会报错但结果出来往往面目全非。我做这套代码时固定了一套单位约定建议新上手的朋友直接照抄。参数符号数值示例单位中心波长λ_B1550e-9m有效折射率n_eff1.45无光栅长度L20e-3m总啁啾量C2e-9m即2 nm折射率调制深度δn2e-4无分段数N1000无波长扫描范围span6e-9m波长采样点数Nλ2001无这里重点说一下总啁啾量C。它的含义是光栅首端到末端的局部布拉格波长一共变化了多少。C 2 nm意味着从一端到另一端布拉格波长连续从1550 nm渐变到1552 nm或反过来。光栅长度L和总啁啾量C的比值C/L才是真正决定反射带宽和时延斜率的物理量。为什么中心波长选1550 nm、折射率选1.45因为这是常规单模光纤在C波段的典型参数对应光纤通信最常用的窗口。折射率调制深度δn选1e-4到3e-4之间也是实际紫外写入光栅能达到的典型范围。如果你在做的是其他波段或特殊光纤只需要改这些基础值代码结构完全不动。3. Mathcad代码实现与结果解读3.1 完整可运行的仿真代码我直接给出一套完整代码基于传输矩阵法注释写在关键行。这段代码我在R2022b和R2024a上都跑过没有依赖额外工具箱纯基础函数就能运行。% chirped_fbg_sim.m % 啁啾光纤光栅反射谱与群时延仿真传输矩阵法 clear; close all; clc; %% 1. 基本参数 lambda_B 1550e-9; % 中心波长单位 m n_eff 1.45; % 有效折射率 L 20e-3; % 光栅总长度单位 m N 1000; % 分段数 C 2e-9; % 总啁啾量单位 m即 2 nm dn 2e-4; % 折射率调制幅度交流耦合相关 v 1; % 条纹可见度通常取1 %% 2. 波长扫描范围 span 6e-9; % 单侧扫描范围 ±6 nm Nlambda 2001; % 波长采样点数 lambda linspace(lambda_B - span, lambda_B span, Nlambda); %% 3. 传输矩阵求解 r zeros(1, Nlambda); % 反射系数 for k 1:Nlambda lam lambda(k); % 全局失谐量 delta以中心波长为参考 delta 2*pi*n_eff * (1/lam - 1/lambda_B); % 单位矩阵开始 M eye(2); for m 1:N % 第 m 段中心位置每段长度 Lseg Lseg L / N; z (m - 0.5) * Lseg; % 局部布拉格波长线性啁啾 lambda_local lambda_B - C/2 C * z / L; % 直流自耦合系数参考局部布拉格波长 sigma_hat 2*pi*n_eff * (1/lam - 1/lambda_local); % 交流耦合系数 kappa pi * v * dn / lambda_local; % 均匀光栅解析解 gamma sqrt(kappa^2 - sigma_hat^2); S sinh(gamma * Lseg); Cg cosh(gamma * Lseg); % 均匀光栅段传输矩阵 Mseg [Cg - 1i*sigma_hat/gamma*S, -1i*kappa/gamma*S; 1i*kappa/gamma*S, Cg 1i*sigma_hat/gamma*S]; M M * Mseg; end % 从总传输矩阵提取反射系数 r(k) -M(2,1) / M(2,2); end %% 4. 反射率、相位与群时延 R abs(r).^2; phi unwrap(angle(r)); c 3e8; dphi_dlambda gradient(phi, lambda); tau -lambda.^2 / (2*pi*c) .* dphi_dlambda * 1e12; % 单位ps %% 5. 绘图 figure(Position,[100 100 800 500]); subplot(2,1,1); plot((lambda-lambda_B)*1e9, R, b-, LineWidth, 1.5); xlabel(波长失谐 (nm)); ylabel(反射率); title(啁啾光纤光栅反射谱); grid on; xlim([-6 6]); subplot(2,1,2); plot((lambda-lambda_B)*1e9, tau, r-, LineWidth, 1.5); xlabel(波长失谐 (nm)); ylabel(群时延 (ps)); title(啁啾光纤光栅群时延); grid on; xlim([-6 6]);3.2 代码逐段拆解每一步在算什么这段代码核心在双重循环里。外层循环是波长内层循环是空间分段。很多人第一次看会问为什么不直接向量化确实可以向量化但双重循环的写法每一步物理意义直观调试起来也方便作为教学框架更合适。性能优化在第5章再谈。看内层循环几个关键量局部布拉格波长lambda_local的计算。lambda_local lambda_B - C/2 C*z/L当z0时是最短波长端zL时是最长波长端。光栅首端z0对应短波长末端对应长波长这是我自己习惯用的方向。如果你希望长波长端在入射端把C改成负值即可。这个方向不影响反射谱形状但会影响时延斜率的正负。sigma_hat是直流自耦合系数我这里用“入射波长相对于该段局部布拉格波长的失谐”来表示。它的物理含义是这一段光栅的“谐振中心”在哪里入射波长离它有多远。失谐越大耦合越弱反射越低。kappa是交流耦合系数直接正比于折射率调制深度dn。dn越大前后向模式的耦合越强也就是反射越强。这里用了一个近似κ π·v·dn/λ_local。严格公式里还需要考虑模式重叠因子但常规仿真中取1即可。gamma是传播常数差由 kappa 和 sigma_hat 共同决定。注意 gamma 可能是虚数或实数取决于 kappa² - sigma_hat² 的符号。当入射波长远离布拉格条件时sigma_hat 增大Gamma 变成纯虚数sinh/cosh 就自动转化为 sin/cos 的振荡行为。这个是传输矩阵法能同时描述通带和阻带的精髓。矩阵相乘完成后反射系数从M(2,1)和M(2,2)中提取。为什么要用 -M(2,1)/M(2,2)这是传输矩阵的边界条件决定的假设光从光栅前端入射后端没有反向输入只有前向输出。经过代数化简反射系数就落在这个表达式上。这个公式我每次都现推但写代码时直接背下来也行关键是记着不是M(1,2)/M(1,1)那是另一个方向的定义搞反了反射谱是对称的但初相位会差个符号时延结果就错了。3.3 仿真结果怎么解读跑完这段代码你应该看到两个图。反射谱在主带宽内反射率接近1带宽大约对应总啁啾量C的量级这里2 nm啁啾主反射带宽就接近2 nm并叠加一些带内振荡。时延曲线在反射带内是一条近似直线这就是啁啾光栅最标志性的特征。怎么确认结果是对的我有一个很好用的检查方法把C改成0这段代码就退化成均匀光栅。此时反射谱应该变成一个窄尖峰时延曲线在带内应该非常平缓而不是斜线。如果C0时跑出来还是宽的或斜的那说明啁啾逻辑写错了。先把均匀光栅调正确再加啁啾这种“最小可复现”的调试思路能帮你快速定位问题。再把C的方向反过来改成-2e-9反射谱形状不变但时延斜率会反向。这两个检查都通过基本可以确认代码逻辑没有问题。4. 关键参数对反射谱和时延的影响4.1 啁啾系数带宽与色散之间的权衡总啁啾量C是啁啾光栅设计的核心参数之一。我建议你做一个简单的参数扫描固定L20mm、dn2e-4分别让C等于0.5 nm、1 nm、2 nm、4 nm看反射谱和时延的变化。C增大时反射谱带宽几乎线性变宽但峰值反射率会略有下降。原因不难理解总调制深度dn不变啁啾量越大每个局部波长的“有效分段”越短相当于把同样的折射率调制能量摊到了更宽的波长范围内单点的反射强度自然会下降。同时时延曲线的斜率会变缓也就是单位波长间隔对应的物理反射位置差变小了。这个权衡在色散补偿设计中很关键。需要补偿的带宽越宽就要用更大的C但带来的时延斜率即色散量会变小。如果同时要大带宽和大色散量就只能加长光栅长度L这就是为什么工程用的色散补偿光栅动辄10 cm以上。总啁啾量C反射带宽峰值反射率时延斜率0.5 nm窄约0.5 nm高陡1 nm中等较高中等2 nm较宽略降平缓4 nm很宽明显下降很平缓4.2 光栅长度反射增强与纹波加剧固定C2nm、dn2e-4把L从5mm加到40mm你能看到两个明显趋势。第一反射率整体上升。光栅越长前后向模式的相互作用距离越长同样折射率调制下的反射就越高。这也解释了为什么弱光栅dn很小需要靠加长光栅长度来保证反射率。第二反射谱带内的振荡纹波会变得更加细密。原因是长光栅内部的多反射干涉效应更明显等效于一个更高质量的腔谐振峰更尖锐。对应地时延曲线上的纹波也会增大。这个现象在真实器件制造中是不可忽视的光谱上的纹波会导致色散补偿后的残余色散波动影响系统性能。如果既要长光栅的强反射又想抑制纹波常规做法就是切趾apodization也就是让折射率调制幅度从中心向两端平滑过渡到0削弱两端边界引起的干涉效应。切趾的代价是峰值反射率下降因为两端本来贡献反射的区域被削弱了。4.3 折射率调制深度强反射与旁瓣折射率调制深度dn对反射谱的影响非常直接。dn从5e-5增大到5e-4反射峰会从不足50%一路冲到接近100%。但与此同时带宽两侧的旁瓣也会迅速抬高反射谱的矩形度变差。这里要区分两种“旁瓣”一种是因为啁啾光栅本身是弱反射布拉格光栅的叠加带外残余反射天然存在另一种是光栅两端骤然的折射率突变引起的F-P干涉旁瓣这个更尖锐也是切趾主要要抑制的对象。实际写入光栅时dn并不是越大越好。过高的紫外曝光剂量会导致光纤损耗增加、折射率调制非线性饱和甚至损伤纤芯。仿真时可以大胆调dn但心里要有数物理上可实现的dn一般在10⁻⁵到10⁻³量级。4.4 切趾函数实操一行代码改变谱形切趾操作在代码里非常简单只需要给每段的κ乘一个窗函数。比如高斯切趾% 在循环外预计算切趾窗 z_all (0.5:N) * L / N; apod exp(-8 * ((z_all - L/2) / L).^2); % 在循环内给 kappa 乘以 apod(m) kappa pi * v * dn / lambda_local * apod(m);高斯切趾会显著压低旁瓣但也会让反射谱边缘变得圆滑带宽略有收缩时延曲线边缘的振荡也会变小。我个人做实际方案预研时通常先把无切趾的结果跑一遍了解器件“极限性能”再加切趾看“实用性能”两者对比更能说明问题。5. 实操中踩过的坑与排查技巧5.1 时延曲线边缘噪声严重怎么办时延曲线的边缘反射带外经常出现剧烈抖动看着像噪声其实是物理上必然的结果反射率极低时反射系数的相位处于“随机游走”状态数值求导自然放大了噪声。这在所有光栅仿真里都会遇到不必怀疑代码错了。处理办法有两种。第一种是显示时限制在反射率较高的区间比如反射率大于0.1的部分才显示时延曲线。第二种是对时延曲线做平滑滤波比如用movmean。但要注意平滑会抹掉真实纹波如果你关心的是带内时延纹波不要用太宽的窗。我一般这样处理% 只绘制反射率足够高的区域 idx R 0.1; plot((lambda(idx)-lambda_B)*1e9, tau(idx), r-);5.2 反射谱不对称是什么原因如果你跑出来的反射谱明显不对称主峰偏向一侧通常有三个原因。第一个原因最隐蔽直流折射率调制平均折射率变化被忽略了。实际光栅的折射率调制可以写成δn_eff(z)·[1 v·cos(...)]其中包含一个直流项。这个直流项会让整个光栅的有效光程发生变化在强光栅或高啁啾时会导致谱不对称。我这套教学代码里把它忽略了做定量分析时需要在sigma_hat里补上直流项。第二个原因是啁啾方向与局部布拉格波长的配合关系。如果你把lambda_local的顺序写反了反射谱主体仍然存在但相位变化不同时延斜率会反向叠加直流项后会出现轻微不对称。第三个原因是波长扫描范围不够。如果扫描范围边缘刚好覆盖了反射带边缘由于采样点不足边缘形态会失真。适当加大span或者增加Nlambda可以缓解。5.3 计算效率低下的优化方案双重循环写起来直观但Nlambda2001、N1000时循环次数是200万次在普通笔记本上要跑一段时间。如果你做了参数扫描比如遍历10组C时间会成倍上升这时候优化就很有必要。最有效的优化是消除内层循环把Mseg的构造向量化。具体做法是把所有段的sigma_hat和kappa预计算成数组然后利用matlab的数组运算一次性生成所有段的矩阵元素再用循环或cellfun连乘。不过连乘本身没法完全避免循环实际提升有限。更实用的方案是减少Nlambda比如从2001降到1001时间几乎减半而谱形变化不大。如果机器有多核外层波长循环可以直接用parfor替代for并行计算能带来几倍提速。注意parfor里不要依赖循环迭代顺序传输矩阵法每段之间没有历史依赖所以很适合并行。5.4 快速自查清单代码跑出异常结果时我建议按表格逐项排查效率最高检查项常见错误正确状态单位nm与m混写所有长度统一为m啁啾方向z起始方向写反时延斜率与预期一致分段数NN太小谱型粗糙N≥500推荐1000相位解包裹忘用unwrap时延无巨大跳变尖峰反射系数提取用错矩阵元素反射率峰值≤1时延单位忘记乘1e12时延值在几十ps量级这几项检查大概能覆盖9成的新手问题。尤其是单位问题我见过太多人把1550 nm直接写成1550结果反射谱跑到可见光波段去了找了半天bug才发现是单位没换算。最后再分享一个小技巧这套传输矩阵框架最大的价值不是“算啁啾光栅”而是它的可扩展性。我在实际项目里经常只改几行代码就把它变成相移光栅仿真在矩阵中插入一个纯相位矩阵、取样光栅仿真把折射率调制乘上周期采样函数、甚至是光纤法布里-珀罗腔仿真两段光栅中间夹一段无光栅区。每次接到新课题我都会先回到这个最基础的框架再往上叠加新的物理效应。你在跑通上面代码之后不妨也试着加一个切趾函数、插入一个相移区看看谱形怎么变这才是真正把传输矩阵法玩熟了。本文还有配套的精品资源点击获取
啁啾光纤光栅仿真入门：基于传输矩阵法的反射谱与群时延分析
简介光纤光栅是光纤通信与传感领域的核心器件其反射谱和群时延特性直接决定了色散补偿、脉冲整形等应用的性能。对于周期渐变的啁啾光栅传统解析方法难以处理而传输矩阵法通过将光栅分段近似为均匀光栅并连乘传输矩阵能高效仿真宽带反射谱与线性群时延。Matlab凭借强大的矩阵运算和信号处理函数成为实现该算法的理想工具。本文从耦合模理论出发介绍传输矩阵法原理、关键参数啁啾量、光栅长度、折射率调制深度的影响规律并结合可运行代码展示如何提取反射谱相位、计算群时延以及规避相位解包裹、单位换算等常见陷阱。掌握这一仿真框架可进一步扩展至切趾、相移光栅等复杂结构为光栅器件的自主设计与工程调试提供扎实基础。 前几天有位朋友在交流群里丢了个压缩包文件名写着“zhoujiu.zip_matlab光栅时延_matlab啁啾光栅的反射谱、时延”。群友都心知肚明这类资源多半是从某论坛或淘宝店转了几手的东西代码能跑但注释乱、公式缺、参数写死真正想拿来用在课题里反而要花大量时间反推。与其研究别人的古董代码不如把啁啾光纤光栅仿真这件事从原理到落地完整捋一遍。这篇内容我尽量讲透为什么啁啾光栅的反射谱是宽带、群时延为什么是线性以及如何用matlab自己搭一套传输矩阵法仿真。适合光纤通信、光纤传感方向的学生以及刚接触光栅器件仿真、想快速上手又不想只做“调参侠”的工程师。1. 啁啾光纤光栅入门周期在变的“彩虹镜子”1.1 从均匀光栅到啁啾光栅周期在变先回到光纤光栅本身。所谓光纤光栅就是利用紫外光在光纤纤芯里写下周期性的折射率调制形成一层“波长选择反射镜”。均匀光栅的周期是恒定的根据布拉格条件 λ_B 2·n_eff·Λ只有满足这个波长条件的光才会被强烈反射其他波长直接透过。所以均匀光栅的反射谱是一个极窄的尖峰典型半高全宽也就0.1~0.3 nm量级非常适合做窄带滤波和单波长激光器选模。啁啾光栅的区别在于周期不再恒定而是沿光栅长度方向线性变化最常见的是线性啁啾。也就是说光栅前端的局部布拉格波长是λ_B1后端是λ_B2中间连续过渡。这就像把一面单色镜子换成了一排按彩虹顺序排列的小镜子每个位置只反射对应的一小段波长。于是整体反射谱被展宽成一个宽带而且不同波长的光在不同位置反射天然带来了额外的时延。用生活化一点的方式来理解均匀光栅像百米赛跑的直线跑道所有运动员同时出发、同时冲线啁啾光栅则像一条从山脚蜿蜒到山顶的盘山道不同速度的运动员被分到不同路段跑完所用的时间自然不一样。这就是啁啾光栅时延特性最直观的来源。1.2 啁啾光栅的核心应用色散补偿与脉冲整形反射谱只是回答了“哪些波长被反射”时延则回答“这些波长分别延迟了多少”。两者合在一起才构成啁啾光栅最核心的两个指标。为什么这么看重时延因为啁啾光栅最重要的应用之一就是色散补偿。光纤通信里光脉冲在普通单模光纤中传播时不同波长分量因为群速度不同到达接收端的时间有差异脉冲就被“拉宽”了。这在高码率长距离传输里是致命的。啁啾光栅的做法是把不同波长的分量安排到光栅的不同深度反射短波长分量先返回、长波长分量后返回或者反过来人为制造一个与光纤相反的光程差从而把被拉开的脉冲重新“压”回去。除了色散补偿啁啾光栅还广泛用于超短脉冲整形、宽带滤波、多波长传感等领域。我见过不少做光纤传感的朋友把啁啾光栅当宽带反射器用配合可调谐激光器做解调也是常规操作。所以无论从通信还是传感角度反射谱和时延这两个仿真结果都是后续一切工作的基础。1.3 为什么选择matlab作为仿真工具matlab在这个领域的地位还是挺稳的。矩阵运算天然适合传输矩阵法内置的unwrap、gradient、plot函数让数据处理和可视化非常顺手而且网上现成资料多。但问题也出在“现成”上很多流传的代码要么参数单位混乱要么绕过了关键的物理项跑出来谱型怪怪的新手又不知道问题出在哪。我个人的建议是不要把“能跑”当作目标而是把底层原理弄明白自己搭一套干净的框架。这套框架以后加切趾、加相移、换成取样光栅都只是改几行参数的事。下面我从耦合模理论讲起把传输矩阵法的来龙去脉理清楚。2. 仿真方案选型为什么用传输矩阵法2.1 耦合模理论与传输矩阵法的关系光栅内部的电磁场分析经典路径是用耦合模理论。前向传播模式和反向传播模式之间因为折射率微扰而相互耦合最终得到一组耦合模方程。对于折射率调制均匀的理想光栅这组方程有严格的解析解解出来就是反射谱的表达式。但加上啁啾之后周期沿长度变化局部布拉格波长不再是常数严格解析解就变得很麻烦。工程上最通用的替代方案就是传输矩阵法Transfer Matrix MethodTMM。传输矩阵法的思路可以用一句话概括把整个光栅切成很多小段每一段足够短近似看成均匀光栅然后每段用解析解写成一个2×2矩阵最后把所有矩阵连乘得到整个光栅的传输特性。这就像用很多条直线段去逼近一条曲线切得越细、越接近真实情况只是计算量会上升。对20 mm长的光栅切成1000段每段只有20 µm局部的周期变化完全可以忽略不计每一段当作均匀光栅处理是完全合理的。具体到每一段的2×2矩阵核心物理量是三个失谐量σ表示入射波长偏离该段局部布拉格波长的程度交流耦合系数κ表示折射率调制幅度带来的前后向模式耦合强度段长Lseg即每段的物理长度。有了这三个量就可以构造出标准的均匀光栅传输矩阵再在波长循环里逐段连乘。虽然代码写起来只是几行矩阵乘法但背后是耦合模方程在分段区间上的解析解表达这就是整个方法的数学根基。2.2 时延计算的数学基础反射谱很容易得到反射系数r是一个复数反射率R |r|²。时延则要从r的相位里提取。设反射系数的相位为φ群时延的定义是τ dφ/dω。因为仿真里用的自变量是波长而不是角频率需要做一个变量替换dω -2πc/λ² · dλ换算过来就是τ -λ²/(2πc) · dφ/dλ这个负号经常把人绕晕。我自己的记法是波长增大时角频率减小所以用λ求导时要加负号把方向转回来。数值上dφ/dλ可以用一行gradient(phi, lambda)得到非常方便。这里有一个新手必踩的坑angle函数返回的相位是包裹在[-π, π]区间的直接求导会在跳变处产生巨大的尖峰所以必须先做unwrap。即使如此低反射区域的相位本身噪声很大unwrap之后也可能出现不合理的跳变后面第5章我会讲怎么处理。2.3 关键参数设定与单位处理仿真参数的单位问题是最容易翻车的地方。matlab没有量纲检查nm、m、ps混着写也不会报错但结果出来往往面目全非。我做这套代码时固定了一套单位约定建议新上手的朋友直接照抄。参数符号数值示例单位中心波长λ_B1550e-9m有效折射率n_eff1.45无光栅长度L20e-3m总啁啾量C2e-9m即2 nm折射率调制深度δn2e-4无分段数N1000无波长扫描范围span6e-9m波长采样点数Nλ2001无这里重点说一下总啁啾量C。它的含义是光栅首端到末端的局部布拉格波长一共变化了多少。C 2 nm意味着从一端到另一端布拉格波长连续从1550 nm渐变到1552 nm或反过来。光栅长度L和总啁啾量C的比值C/L才是真正决定反射带宽和时延斜率的物理量。为什么中心波长选1550 nm、折射率选1.45因为这是常规单模光纤在C波段的典型参数对应光纤通信最常用的窗口。折射率调制深度δn选1e-4到3e-4之间也是实际紫外写入光栅能达到的典型范围。如果你在做的是其他波段或特殊光纤只需要改这些基础值代码结构完全不动。3. Mathcad代码实现与结果解读3.1 完整可运行的仿真代码我直接给出一套完整代码基于传输矩阵法注释写在关键行。这段代码我在R2022b和R2024a上都跑过没有依赖额外工具箱纯基础函数就能运行。% chirped_fbg_sim.m % 啁啾光纤光栅反射谱与群时延仿真传输矩阵法 clear; close all; clc; %% 1. 基本参数 lambda_B 1550e-9; % 中心波长单位 m n_eff 1.45; % 有效折射率 L 20e-3; % 光栅总长度单位 m N 1000; % 分段数 C 2e-9; % 总啁啾量单位 m即 2 nm dn 2e-4; % 折射率调制幅度交流耦合相关 v 1; % 条纹可见度通常取1 %% 2. 波长扫描范围 span 6e-9; % 单侧扫描范围 ±6 nm Nlambda 2001; % 波长采样点数 lambda linspace(lambda_B - span, lambda_B span, Nlambda); %% 3. 传输矩阵求解 r zeros(1, Nlambda); % 反射系数 for k 1:Nlambda lam lambda(k); % 全局失谐量 delta以中心波长为参考 delta 2*pi*n_eff * (1/lam - 1/lambda_B); % 单位矩阵开始 M eye(2); for m 1:N % 第 m 段中心位置每段长度 Lseg Lseg L / N; z (m - 0.5) * Lseg; % 局部布拉格波长线性啁啾 lambda_local lambda_B - C/2 C * z / L; % 直流自耦合系数参考局部布拉格波长 sigma_hat 2*pi*n_eff * (1/lam - 1/lambda_local); % 交流耦合系数 kappa pi * v * dn / lambda_local; % 均匀光栅解析解 gamma sqrt(kappa^2 - sigma_hat^2); S sinh(gamma * Lseg); Cg cosh(gamma * Lseg); % 均匀光栅段传输矩阵 Mseg [Cg - 1i*sigma_hat/gamma*S, -1i*kappa/gamma*S; 1i*kappa/gamma*S, Cg 1i*sigma_hat/gamma*S]; M M * Mseg; end % 从总传输矩阵提取反射系数 r(k) -M(2,1) / M(2,2); end %% 4. 反射率、相位与群时延 R abs(r).^2; phi unwrap(angle(r)); c 3e8; dphi_dlambda gradient(phi, lambda); tau -lambda.^2 / (2*pi*c) .* dphi_dlambda * 1e12; % 单位ps %% 5. 绘图 figure(Position,[100 100 800 500]); subplot(2,1,1); plot((lambda-lambda_B)*1e9, R, b-, LineWidth, 1.5); xlabel(波长失谐 (nm)); ylabel(反射率); title(啁啾光纤光栅反射谱); grid on; xlim([-6 6]); subplot(2,1,2); plot((lambda-lambda_B)*1e9, tau, r-, LineWidth, 1.5); xlabel(波长失谐 (nm)); ylabel(群时延 (ps)); title(啁啾光纤光栅群时延); grid on; xlim([-6 6]);3.2 代码逐段拆解每一步在算什么这段代码核心在双重循环里。外层循环是波长内层循环是空间分段。很多人第一次看会问为什么不直接向量化确实可以向量化但双重循环的写法每一步物理意义直观调试起来也方便作为教学框架更合适。性能优化在第5章再谈。看内层循环几个关键量局部布拉格波长lambda_local的计算。lambda_local lambda_B - C/2 C*z/L当z0时是最短波长端zL时是最长波长端。光栅首端z0对应短波长末端对应长波长这是我自己习惯用的方向。如果你希望长波长端在入射端把C改成负值即可。这个方向不影响反射谱形状但会影响时延斜率的正负。sigma_hat是直流自耦合系数我这里用“入射波长相对于该段局部布拉格波长的失谐”来表示。它的物理含义是这一段光栅的“谐振中心”在哪里入射波长离它有多远。失谐越大耦合越弱反射越低。kappa是交流耦合系数直接正比于折射率调制深度dn。dn越大前后向模式的耦合越强也就是反射越强。这里用了一个近似κ π·v·dn/λ_local。严格公式里还需要考虑模式重叠因子但常规仿真中取1即可。gamma是传播常数差由 kappa 和 sigma_hat 共同决定。注意 gamma 可能是虚数或实数取决于 kappa² - sigma_hat² 的符号。当入射波长远离布拉格条件时sigma_hat 增大Gamma 变成纯虚数sinh/cosh 就自动转化为 sin/cos 的振荡行为。这个是传输矩阵法能同时描述通带和阻带的精髓。矩阵相乘完成后反射系数从M(2,1)和M(2,2)中提取。为什么要用 -M(2,1)/M(2,2)这是传输矩阵的边界条件决定的假设光从光栅前端入射后端没有反向输入只有前向输出。经过代数化简反射系数就落在这个表达式上。这个公式我每次都现推但写代码时直接背下来也行关键是记着不是M(1,2)/M(1,1)那是另一个方向的定义搞反了反射谱是对称的但初相位会差个符号时延结果就错了。3.3 仿真结果怎么解读跑完这段代码你应该看到两个图。反射谱在主带宽内反射率接近1带宽大约对应总啁啾量C的量级这里2 nm啁啾主反射带宽就接近2 nm并叠加一些带内振荡。时延曲线在反射带内是一条近似直线这就是啁啾光栅最标志性的特征。怎么确认结果是对的我有一个很好用的检查方法把C改成0这段代码就退化成均匀光栅。此时反射谱应该变成一个窄尖峰时延曲线在带内应该非常平缓而不是斜线。如果C0时跑出来还是宽的或斜的那说明啁啾逻辑写错了。先把均匀光栅调正确再加啁啾这种“最小可复现”的调试思路能帮你快速定位问题。再把C的方向反过来改成-2e-9反射谱形状不变但时延斜率会反向。这两个检查都通过基本可以确认代码逻辑没有问题。4. 关键参数对反射谱和时延的影响4.1 啁啾系数带宽与色散之间的权衡总啁啾量C是啁啾光栅设计的核心参数之一。我建议你做一个简单的参数扫描固定L20mm、dn2e-4分别让C等于0.5 nm、1 nm、2 nm、4 nm看反射谱和时延的变化。C增大时反射谱带宽几乎线性变宽但峰值反射率会略有下降。原因不难理解总调制深度dn不变啁啾量越大每个局部波长的“有效分段”越短相当于把同样的折射率调制能量摊到了更宽的波长范围内单点的反射强度自然会下降。同时时延曲线的斜率会变缓也就是单位波长间隔对应的物理反射位置差变小了。这个权衡在色散补偿设计中很关键。需要补偿的带宽越宽就要用更大的C但带来的时延斜率即色散量会变小。如果同时要大带宽和大色散量就只能加长光栅长度L这就是为什么工程用的色散补偿光栅动辄10 cm以上。总啁啾量C反射带宽峰值反射率时延斜率0.5 nm窄约0.5 nm高陡1 nm中等较高中等2 nm较宽略降平缓4 nm很宽明显下降很平缓4.2 光栅长度反射增强与纹波加剧固定C2nm、dn2e-4把L从5mm加到40mm你能看到两个明显趋势。第一反射率整体上升。光栅越长前后向模式的相互作用距离越长同样折射率调制下的反射就越高。这也解释了为什么弱光栅dn很小需要靠加长光栅长度来保证反射率。第二反射谱带内的振荡纹波会变得更加细密。原因是长光栅内部的多反射干涉效应更明显等效于一个更高质量的腔谐振峰更尖锐。对应地时延曲线上的纹波也会增大。这个现象在真实器件制造中是不可忽视的光谱上的纹波会导致色散补偿后的残余色散波动影响系统性能。如果既要长光栅的强反射又想抑制纹波常规做法就是切趾apodization也就是让折射率调制幅度从中心向两端平滑过渡到0削弱两端边界引起的干涉效应。切趾的代价是峰值反射率下降因为两端本来贡献反射的区域被削弱了。4.3 折射率调制深度强反射与旁瓣折射率调制深度dn对反射谱的影响非常直接。dn从5e-5增大到5e-4反射峰会从不足50%一路冲到接近100%。但与此同时带宽两侧的旁瓣也会迅速抬高反射谱的矩形度变差。这里要区分两种“旁瓣”一种是因为啁啾光栅本身是弱反射布拉格光栅的叠加带外残余反射天然存在另一种是光栅两端骤然的折射率突变引起的F-P干涉旁瓣这个更尖锐也是切趾主要要抑制的对象。实际写入光栅时dn并不是越大越好。过高的紫外曝光剂量会导致光纤损耗增加、折射率调制非线性饱和甚至损伤纤芯。仿真时可以大胆调dn但心里要有数物理上可实现的dn一般在10⁻⁵到10⁻³量级。4.4 切趾函数实操一行代码改变谱形切趾操作在代码里非常简单只需要给每段的κ乘一个窗函数。比如高斯切趾% 在循环外预计算切趾窗 z_all (0.5:N) * L / N; apod exp(-8 * ((z_all - L/2) / L).^2); % 在循环内给 kappa 乘以 apod(m) kappa pi * v * dn / lambda_local * apod(m);高斯切趾会显著压低旁瓣但也会让反射谱边缘变得圆滑带宽略有收缩时延曲线边缘的振荡也会变小。我个人做实际方案预研时通常先把无切趾的结果跑一遍了解器件“极限性能”再加切趾看“实用性能”两者对比更能说明问题。5. 实操中踩过的坑与排查技巧5.1 时延曲线边缘噪声严重怎么办时延曲线的边缘反射带外经常出现剧烈抖动看着像噪声其实是物理上必然的结果反射率极低时反射系数的相位处于“随机游走”状态数值求导自然放大了噪声。这在所有光栅仿真里都会遇到不必怀疑代码错了。处理办法有两种。第一种是显示时限制在反射率较高的区间比如反射率大于0.1的部分才显示时延曲线。第二种是对时延曲线做平滑滤波比如用movmean。但要注意平滑会抹掉真实纹波如果你关心的是带内时延纹波不要用太宽的窗。我一般这样处理% 只绘制反射率足够高的区域 idx R 0.1; plot((lambda(idx)-lambda_B)*1e9, tau(idx), r-);5.2 反射谱不对称是什么原因如果你跑出来的反射谱明显不对称主峰偏向一侧通常有三个原因。第一个原因最隐蔽直流折射率调制平均折射率变化被忽略了。实际光栅的折射率调制可以写成δn_eff(z)·[1 v·cos(...)]其中包含一个直流项。这个直流项会让整个光栅的有效光程发生变化在强光栅或高啁啾时会导致谱不对称。我这套教学代码里把它忽略了做定量分析时需要在sigma_hat里补上直流项。第二个原因是啁啾方向与局部布拉格波长的配合关系。如果你把lambda_local的顺序写反了反射谱主体仍然存在但相位变化不同时延斜率会反向叠加直流项后会出现轻微不对称。第三个原因是波长扫描范围不够。如果扫描范围边缘刚好覆盖了反射带边缘由于采样点不足边缘形态会失真。适当加大span或者增加Nlambda可以缓解。5.3 计算效率低下的优化方案双重循环写起来直观但Nlambda2001、N1000时循环次数是200万次在普通笔记本上要跑一段时间。如果你做了参数扫描比如遍历10组C时间会成倍上升这时候优化就很有必要。最有效的优化是消除内层循环把Mseg的构造向量化。具体做法是把所有段的sigma_hat和kappa预计算成数组然后利用matlab的数组运算一次性生成所有段的矩阵元素再用循环或cellfun连乘。不过连乘本身没法完全避免循环实际提升有限。更实用的方案是减少Nlambda比如从2001降到1001时间几乎减半而谱形变化不大。如果机器有多核外层波长循环可以直接用parfor替代for并行计算能带来几倍提速。注意parfor里不要依赖循环迭代顺序传输矩阵法每段之间没有历史依赖所以很适合并行。5.4 快速自查清单代码跑出异常结果时我建议按表格逐项排查效率最高检查项常见错误正确状态单位nm与m混写所有长度统一为m啁啾方向z起始方向写反时延斜率与预期一致分段数NN太小谱型粗糙N≥500推荐1000相位解包裹忘用unwrap时延无巨大跳变尖峰反射系数提取用错矩阵元素反射率峰值≤1时延单位忘记乘1e12时延值在几十ps量级这几项检查大概能覆盖9成的新手问题。尤其是单位问题我见过太多人把1550 nm直接写成1550结果反射谱跑到可见光波段去了找了半天bug才发现是单位没换算。最后再分享一个小技巧这套传输矩阵框架最大的价值不是“算啁啾光栅”而是它的可扩展性。我在实际项目里经常只改几行代码就把它变成相移光栅仿真在矩阵中插入一个纯相位矩阵、取样光栅仿真把折射率调制乘上周期采样函数、甚至是光纤法布里-珀罗腔仿真两段光栅中间夹一段无光栅区。每次接到新课题我都会先回到这个最基础的框架再往上叠加新的物理效应。你在跑通上面代码之后不妨也试着加一个切趾函数、插入一个相移区看看谱形怎么变这才是真正把传输矩阵法玩熟了。本文还有配套的精品资源点击获取