简介面向光学工程与MATLAB仿真的光栅反射特性资源包内容覆盖均匀布拉格光栅、切趾光栅、相移布拉格光栅、啁啾光栅与取样光栅的反射/透射谱模拟适用于课程设计、课题预研以及光栅基础理论学习。压缩包内共12个文件包含6个可直接运行的.m脚本和6张对应结果的.jpg图片整体仅88KB结构清晰便于快速定位不同光栅类型对应的仿真程序。脚本中运用傅里叶变换、光栅方程等原理构造光栅模型可修改参数观察谱线变化图片直观展示均匀、切趾、取样等光栅的反射谱形态方便对照理解算法输出。目前已有297人浏览学习特别适合刚接触光学仿真、希望通过MATLAB快速上手光栅计算的学习者。借助这些脚本读者不仅能复现多种光栅反射谱还能掌握二维傅里叶变换、光栅掩模生成等关键代码思路为后续自建光栅模型或扩展至更复杂光学系统模拟打下基础。1. 为什么光栅反射仿真答案藏在取样结构里做光纤通信和激光器设计的人几乎都会撞上同一个问题普通均匀布拉格光栅反射谱只有单一峰值做单波长选择没问题但一到波分复用、多波长激光器这类场景单峰就成了瓶颈。取样光栅Sampled Grating用周期性通断调制把反射谱打散成一组离散梳状峰相当于一个光域梳状滤波器。问题是这种结构的折射率调制沿轴向剧烈变化解析解基本失效手算只能给出包络。拿 MATLAB 做光栅反射仿真最务实的路径不是去求解复杂的亥姆霍兹方程而是用传输矩阵法把非均匀光栅切成几百上千段均匀小段逐段累乘几分钟内得到完整反射谱。本文从耦合模方程讲起落在取样光栅的 MATLAB 实现、参数扫描与数值收敛性检查上。2. 光栅反射的 MATLAB 物理建模从耦合模方程到传输矩阵2.1 为什么传输矩阵法能算取样光栅均匀光栅的反射谱可以用耦合模方程解析求解但取样光栅的折射率调制函数是一个矩形脉冲序列——光栅段与空白段交替出现这就是一个强非均匀结构解析解只对均匀段成立。更麻烦的是实际设计的取样光栅可能还要叠加切趾、相移或啁啾任何一项都会破坏解析公式的前提。传输矩阵法Transfer Matrix Method, TMM的思路比较直接把光栅沿传播方向切成 N 段每一段足够短以至于段内折射率可视为常量。对每一段写出一个 2×2 矩阵描述前向波和后向波在该段两端的传播关系最后把所有矩阵按传播顺序相乘。这么做的好处是无论折射率分布多复杂只要切得足够细传输矩阵的乘积就能逼近真实响应。对 MATLAB 而言矩阵乘法本身就是原生运算几百段连乘只是几毫秒的事所以 TMM 是取样光栅仿真的最优选型没有之一。2.2 MATLAB 里构造折射率切片与光栅矩阵2.2.1 全局仿真参数与切片初始化动手写代码之前先把仿真参数定清楚。常见做法是把布拉格波长、有效折射率、调制深度、光栅长度、取样周期、占空比全部集中在脚本头部方便后续扫描。以下是最小参数集% 全局仿真参数 lambda_B 1550e-9; % 布拉格波长单位米 n_eff 1.45; % 有效折射率光纤基模典型值 delta_n 1e-4; % 折射率调制深度均匀段内 L_g 1e-3; % 光栅总长度单位米 N_seg 500; % 沿光栅方向的分段数 dz L_g / N_seg; % 每段长度单位米 lambda linspace(1548e-9, 1552e-9, 2001); % 波长扫描范围与点数这里的核心参数是N_seg它决定仿真的空间分辨率。N_seg太小相位积累误差会让反射峰位置偏移N_seg太大矩阵相乘次数增加数值噪声可能被放大。经验法则是每段长度要小于光栅周期的 1/10也就是dz lambda_B / (2 * n_eff) / 10算下来 1mm 光栅最少需要 200 段实际取 500 到 1000 比较稳妥。2.2.2 单段传输矩阵的推导与实现先看均匀段的情况。耦合模理论给出该段的传输矩阵形式% 单段均匀光栅的传输矩阵 function M uniform_seg_matrix(kappa, delta_beta, dz) % kappa: 交流耦合系数单位 1/m % delta_beta: 失谐量波长偏离布拉格条件的传播常数差单位 1/m gamma sqrt(kappa^2 - delta_beta^2 1i * eps); % 加 eps 防止复数域奇异 S sinh(gamma * dz); C cosh(gamma * dz); M [C - 1i * delta_beta .* S ./ gamma, -1i * kappa .* S ./ gamma; 1i * kappa .* S ./ gamma, C 1i * delta_beta .* S ./ gamma]; end这个矩阵的四个元素分别对应前向波透射、前向波反射、后向波透射和后向波反射。gamma的开方操作决定了矩阵元素是振荡|delta_beta| kappa还是指数增长|delta_beta| kappa这是布拉格反射“带内反射、带外透射”物理本质的数学来源。delta_beta的计算公式是2 * pi * n_eff * (1 ./ lambda - 1 / lambda_B)由波长偏离布拉格条件产生。2.3 最小可运行样例与自检规则把所有段的矩阵连乘起来得到整个光栅的传输矩阵。反射率的定义是后向波系数与前向波系数之比在入射端的投影。完整的最小脚本如下% 均匀光栅反射谱计算最小可运行 lambda_B 1550e-9; n_eff 1.45; delta_n 1e-4; L_g 1e-3; N_seg 500; dz L_g / N_seg; lambda linspace(1549e-9, 1551e-9, 1001); kappa pi * delta_n / lambda_B; % 交流耦合系数 R zeros(size(lambda)); for i 1:length(lambda) delta_beta 2 * pi * n_eff * (1 / lambda(i) - 1 / lambda_B); M_total eye(2); % 单位矩阵作为累乘起点 for j 1:N_seg M_total M_total * uniform_seg_matrix(kappa, delta_beta, dz); end % 反射系数 r -M(2,1)/M(2,2)注意矩阵元素顺序 r -M_total(2,1) / M_total(2,2); R(i) abs(r)^2; end plot((lambda - lambda_B) * 1e9, R); xlabel(波长偏移 (pm)); ylabel(反射率);运行后会看到典型的均匀光栅反射峰中心波长处反射率接近 1两侧出现旁瓣。自检规则有两条第一M_total的行列式必须等于 1无损耗介质中能量守恒如果偏离超过1e-6说明分段数不够或数值溢出第二中心波长反射率不应超过 1一旦超过就是矩阵元素计算有误回头检查gamma的复数分支。3. 取样光栅的 MATLAB 实现占空比与采样周期如何决定通道数3.1 取样光栅的物理结构决定梳状谱取样光栅的本质是在均匀光栅上叠加一个周期性的通断调制。设采样周期为P_s占空比为D则每个周期内有光栅的部分长度为D * P_s空白段长度为(1 - D) * P_s。从傅里叶角度看这种周期性调制会在波矢空间产生一系列离散的傅里叶分量每个分量对应一个反射峰峰间隔由采样周期决定参数符号典型值对反射谱的影响采样周期P_s50~200 μm通道间隔周期越短间隔越大占空比D0.2~0.8通道包络宽度与通道数量总光栅长度L_g0.5~5 mm每个通道的带宽调制深度delta_n0.5e-4~5e-4整体反射强度通道间隔的理论公式是delta_lambda lambda_B^2 / (2 * n_eff * P_s)。以 1550nm 波长、1.45 有效折射率、100μm 采样周期为例算出来通道间隔约 8.3nm这个值在 MATLAB 里跑通后可以作为第一道校验。3.2 在 TMM 框架中加入取样调制取样光栅的折射率调制函数可以写成矩形波形式有光栅段delta_n取正常值空白段取 0。实现方式是在分段循环里判断当前位置是否落在有光栅段内function delta_n_profile sampled_grating_profile(N_seg, P_s, D, L_g, delta_n) % 生成取样光栅的折射率调制空间分布 dz L_g / N_seg; delta_n_profile zeros(1, N_seg); % 当前段中心距离光栅起点的位置 for j 1:N_seg z (j - 0.5) * dz; % 判断 z 落在采样周期内的哪个位置 z_mod mod(z, P_s); if z_mod D * P_s delta_n_profile(j) delta_n; else delta_n_profile(j) 0; end end end这里z_mod D * P_s的判断方式决定了调制是“先开后关”。P_s必须能被dz整除否则采样周期的边界会产生额外相位误差表现为谱线底部出现细碎毛刺。如果P_s不是dz的整数倍可以在生成分布前做一次round(P_s / dz)取整处理。注意取样光栅不能用均匀段的单一kappa值走全程因为空白段的kappa为 0该段的传输矩阵会退化为纯传播矩阵。这恰恰是 TMM 的优势每段的矩阵根据该段的局部参数实时计算天然适配非均匀结构。3.3 一个完整采样仿真脚本与结果判读把上面两块拼起来就可以计算取样光栅反射谱。脚本结构是外层波长扫描、内层分段累乘但取样光栅还要在中层做一次折射率分布读取% 取样光栅反射谱完整计算 lambda_B 1550e-9; n_eff 1.45; delta_n 1e-4; L_g 2e-3; N_seg 1000; dz L_g / N_seg; P_s 100e-6; D 0.5; lambda linspace(1540e-9, 1560e-9, 4001); prof sampled_grating_profile(N_seg, P_s, D, L_g, delta_n); R zeros(size(lambda)); for i 1:length(lambda) delta_beta 2 * pi * n_eff * (1 / lambda(i) - 1 / lambda_B); M_total eye(2); for j 1:N_seg if prof(j) 0 kappa pi * prof(j) / lambda_B; else kappa 0; end M_total M_total * uniform_seg_matrix(kappa, delta_beta, dz); end r -M_total(2,1) / M_total(2,2); R(i) abs(r)^2; end plot((lambda - lambda_B) * 1e9, R);运行后能看到一排反射峰峰包络的形状接近 sinc 函数中心通道反射率最高向两侧衰减。判读结果时留意三件事第一通道间隔是否与理论值一致偏差超过 2% 就要检查P_s和n_eff的设置第二中心波长是否仍在 1550nm偏移了说明delta_beta的参考波长有误第三最外侧通道是否出现不对称这通常是取样段的边界处产生了非预期的相移可以尝试把取样边界对齐到网格点上。4. 光栅反射谱的特征参数扫描中心波长、旁瓣抑制与通道均匀性4.1 三个必须调好的参数及其灵敏度行为取样光栅的设计本质上是调节三个参数采样周期P_s、占空比D、调制深度delta_n。它们分别控制通道间隔、包络宽度、反射强度彼此独立适合用 MATLAB 的for循环逐参数扫描。以占空比为例扫描脚本的核心结构如下% 占空比扫描观察通道包络变化 D_list [0.3, 0.5, 0.7]; colors [r, b, k]; figure; for k 1:length(D_list) prof sampled_grating_profile(N_seg, P_s, D_list(k), L_g, delta_n); R zeros(size(lambda)); for i 1:length(lambda) delta_beta 2 * pi * n_eff * (1 / lambda(i) - 1 / lambda_B); M_total eye(2); for j 1:N_seg if prof(j) 0 kappa pi * prof(j) / lambda_B; else kappa 0; end M_total M_total * uniform_seg_matrix(kappa, delta_beta, dz); end r -M_total(2,1) / M_total(2,2); R(i) abs(r)^2; end plot((lambda - lambda_B) * 1e9, R, colors(k)); hold on; end扫描D从 0.3 到 0.7会发现两件事其一占空比越大中心通道反射峰越窄旁瓣越多其二包络零点的位置随占空比改变第一零点的偏移量近似等于lambda_B^2 / (2 * n_eff * D * P_s)。这对应一个经典 trade-off占空比小的时候通道数多但平坦度差占空比大的时候中心通道突出但边通道衰减快。设计多波长激光器时通常希望边通道尽量平坦这时占空比取 0.5 附近是安全起点。4.2 反射率大于 1 的仿真陷阱取样光栅仿真里最常出现的“假信号”是反射率超过 1。出现这个结果不是因为物理上实现了增益而是传输矩阵数值出错。常见原因有两个一是gamma在delta_beta^2 kappa^2时变成纯虚数如果sinh和cosh用了复数输入而没有正确处理符号矩阵元素会出现非物理的指数增长二是分段数太少导致每一段的相位积累不连续矩阵相乘后数值噪声被放大。排查方法很直接先跑一遍均匀光栅把解析反射率公式tanh^2(kappa * L)的变体与 TMM 结果对比。两者在中心波长处的差应小于 1%否则优先修正gamma的计算分支。其次是检查总传输矩阵的行列式det(M_total)偏离 1 超过1e-8时就该考虑增加N_seg。4.3 用窗函数压制旁瓣取样光栅的反射谱天然带 sinc 包络旁瓣第一个旁瓣只比主峰低约 9.5dB这在滤波应用中不够干净。工程上的标准做法是引入切趾apodization也就是让调制深度沿光栅长度方向按窗函数变化% 高斯切趾调制深度沿长度方向呈高斯分布 z (0.5:N_seg) * dz; w exp(-((z - L_g/2).^2) / (2 * (L_g/5)^2)); % 高斯窗 prof_apod prof .* w; % 对调制深度加权高斯窗的宽度取L_g/5时主瓣展宽约 20%但旁瓣可以压到 25dB 以下。注意切趾函数作用在delta_n上而不是kappa上两者是线性关系所以直接加权没问题。切趾后的取样光栅通道间隔不变各通道带宽展宽取舍标准是滤波器类应用追求旁瓣抑制优先切趾激光器类应用追求窄线宽则减小切趾强度。5. 仿真发散与收敛性把分段数和扫描步长调到什么程度最稳5.1 分段数、波长步长与谱线宽度的三角关系取样光栅仿真出现发散的表现有两种一是反射谱出现高频振荡相邻波长点的反射率跳变幅度远超物理预期二是增加分段数后反射谱不收敛而是持续变化。前者通常是波长扫描步长太大谱线太窄采样点不足后者才是真正的数值发散。这里有一个三角约束关系。设每个通道的带宽约为delta_lambda_ch lambda_B^2 / (2 * n_eff * L_g)波长扫描步长要小于这个带宽的 1/5而分段数要保证每段相位积累小于 π/10即dz lambda_B / (20 * n_eff)。以 2mm 光栅、1550nm 波长为例通道带宽约 0.4nm波长步长要小于 0.08nm分段数要超过2e-3 / 5.34e-8约 3745 段。工程上取 5000 段起步比较稳妥。% 收敛性对比矩阵行列式偏差随分段数的变化 N_list [500, 1000, 2000, 5000]; det_dev zeros(size(N_list)); for k 1:length(N_list) N_seg N_list(k); dz L_g / N_seg; M_total eye(2); for j 1:N_seg M_total M_total * uniform_seg_matrix(kappa, delta_beta, dz); end det_dev(k) abs(det(M_total) - 1); end运行这段代码会发现det_dev随N_seg增大会逐渐下降到一个平台但如果继续增大到几万段det_dev反而可能反弹。这是因为浮点数的舍入误差在几千次矩阵相乘后会累积放大。因此收敛性检查不能只看分段数还要设一个行列式偏差阈值。5.2 三招自检仿真发散第一招零反射极限检查。把delta_n设成 0整个光栅变成纯波导反射率必须在所有波长上恒等于 0。如果出现非零反射说明空白段的矩阵赋值出了问题。第二招是能量守恒检查所有波长点的R T都必须等于 1无损耗模型任何一处偏离都表明该波长点的矩阵元素计算有误。第三招是对照解析极限当占空比 D1 时取样光栅退化为均匀光栅此时反射谱峰值应与均匀光栅解析公式tanh^2(kappa * L)一致。% 自动发散检查脚本 assert(max(R) 1 1e-10, 反射率超过1检查gamma分支); M_total_det det(M_total); assert(abs(M_total_det - 1) 1e-6, 传输矩阵行列式偏离1); % 零调制检查 prof_zero zeros(1, N_seg); R_zero calc_reflectance(prof_zero, lambda, lambda_B, n_eff, dz); % 自定义函数 assert(max(R_zero) 1e-12, 零调制时反射率非零);这三个断言建议直接写进仿真脚本的尾部每次参数调整后自动跑一遍。取样光栅设计流程中参数改动频繁人工目测反射谱很容易漏掉数值问题自动检查能在第一时间暴露矩阵计算的异常。6. 通道数翻倍与平坦化取样光栅的进阶仿真验证均匀取样光栅的通道包络受 sinc 函数限制通道数基本由占空比决定。想在不增大P_s的前提下增加可用通道数业界标准做法是引入相位采样Phase Sampling即在每个采样周期内给光栅段附加不同的相位偏移让傅里叶分量的幅度重新分配。MATLAB 里可以用优化工具箱配合 TMM 计算来设计相位序列每个采样周期一个相位值作为优化变量目标函数是通道间反射率的均方差最小化约束为每个相位的取值范围在 0 到 2π 之间。ga函数在这个规模下通常几百代就能收敛比手工调整占空比高效得多。相位采样仿真落地时把每个采样周期的均匀段矩阵左乘一个对角相位矩阵即可% 给第 k 个采样周期附加相位 exp(i*phi_k) M_phase [exp(1i * phi_k), 0; 0, exp(-1i * phi_k)]; M_total M_total * M_phase;签名phi_k的优化结束后通道平坦度最高与最低通道反射率之差能从 5dB 压到 0.8dB 以内代价是中心波长反射率整体下降 1~2dB。最后一步用谱线位置做交叉验证把仿真得到的通道间隔与理论值lambda_B^2 / (2 * n_eff * P_s)逐通道比对偏差超过 0.1% 就重新检查n_eff与温度相关的折射率修正量。确认通道位置符合设计预期后再把相位分布连同结构参数一起保存为.mat文件供后续掩膜版制作或加工公差分析使用。本文还有配套的精品资源点击获取