★ FEATURED ARTICLE 简介一份面向受激布里渊散射SBS光纤传感研究的MATLAB仿真压缩包聚焦泵浦光、信号光与布里渊散射光之间的三波耦合方程求解适合非线性光学、分布式光纤传感方向的初学者和研究人员。SBS源于光与声波的非线性相互作用在分布式温度/应变传感中应用广泛脚本正是对该过程的数值化模拟。包内仅含1个m脚本文件整体大小约2KB代码结构精简便于快速运行、修改参数与二次开发。已有609人学习下载说明该脚本在相关课题中有一定参考价值。通过调整泵浦光功率、光纤长度、光波频率等参数可直观观察布里渊散射强度变化并对照理论结果验证三波耦合模型脚本可作为光纤传感系统仿真、课程设计或论文复现的实用起点。1. 受激布里渊散射仿真SBSa.m 把三波耦合方程变成可跑的实验受激布里渊散射Stimulated Brillouin Scattering, SBS在光纤里是典型的低门槛高翻车率非线性效应阈值只有毫瓦级增益谱窄到几十兆赫还牵着一根看不见的声波耦合通道。做光纤传感的都知道布里渊频移对温度和应变敏感这正是全分布式传感的理论根基。SBSa.m 这个 MATLAB 脚本把三波耦合方程做成可改参数的数值求解器泵浦功率、光纤长度、增益线宽都在文件头部定义跑完直接输出泵浦耗尽曲线和布里渊增益谱。对想验证模型、标定传感系数、或者给 BOTDA 方案打底的人这个脚本能省掉大量重复造轮子的时间。2. 三波耦合方程的数值求解从声波场物理到 SBSa.m 迭代骨架2.1 泵浦光、斯托克斯光与声波场三波到底在耦合什么SBS 的物理图像可以这样理解强泵浦光在光纤中传播时通过电致伸缩效应在纤芯内激发一个声波场。这个声波场本质上是随声速传播的密度光栅折射率被周期性调制泵浦光在这个运动光栅上发生布拉格衍射产生频率下移的斯托克斯光。斯托克斯频率比泵浦低一个布里渊频移ν_B 2·n·V_a / λ_p对 1550 nm 窗口的标准单模光纤n≈1.45石英声速 V_a≈5960 m/s算出来 ν_B≈11.14 GHz。但实际 G.652 光纤因纤芯掺锗声速和折射率都有偏移实测通常在 10.811.0 GHz 附近所以仿真里别把理论值当死数这一点在避坑章节我会专门展开。三波耦合方程的三波指泵浦光、斯托克斯光和声波场。严格形式是三个复数场振幅的偏微分方程组∂E_p/∂z (iω_p/(2nc))·Q·E_s ∂E_s/∂z -(iω_s/(2nc))·Q*·E_p ∂Q/∂t Γ_B·Q i·g₂·E_p·E_s*这里 Q 是声波场复振幅Γ_B 是声波阻尼率与增益线宽直接相关Γ_B π·Δν_B。方程组的本质是三个场互相喂能量泵浦和斯托克斯通过拍频驱动声波声波反过来把泵浦光子散射成斯托克斯光子这个循环建立起来的稳态就是布里渊增益的源头。SBSa.m 里实际用了稳态近似把 ∂Q/∂t 置零Q 直接用 E_p·E_s* 的代数表达式替换再带回前两个方程得到简化后的泵浦-斯托克斯强度耦合方程dI_p/dz -g_B·I_p·I_s - α·I_p dI_s/dz -g_B·I_p·I_s α·I_s这里 g_B 是布里渊增益系数本身是频率失谐的函数呈洛伦兹线型g_B(ν) g_B0·(Δν_B/2)² / [(ν - ν_B)² (Δν_B/2)²]峰值增益 g_B0 在石英光纤中约 35×10⁻¹¹ m/W线宽 Δν_B 典型值 3050 MHz。斯托克斯方程里 α 项符号为正是因为它反向传播在 z 坐标看能量从泵浦转移到斯托克斯光纤损耗对斯托克斯来说表现为传播方向上的衰减。把这两个方程的方向约定搞清楚是读懂 SBSa.m 迭代代码的前提。2.2 SBSa.m 的离散化与固定点迭代这两个方程看起来简单但它是两点边值问题泵浦的边界条件在 z0斯托克斯的边界条件在 zL不能直接往前积分一步到位必须迭代。SBSa.m 用的常见做法是固定点迭代先给一组初值正向更新泵浦、反向更新斯托克斯反复刷直到两条功率曲线不再变化。下面是脚本核心循环的骨架实际 SBSa.m 的逻辑与此一致% 空间网格 Nz 2001; % 网格数决定空间步长 z linspace(0, L, Nz); dz z(2) - z(1); % 边界条件泵浦在 z0 入射斯托克斯在 zL 入射 Ip zeros(1, Nz); Is zeros(1, Nz); Ip(1) Pp_in / Aeff; % 泵浦入纤功率换算成光强 Is(Nz) Ps_in / Aeff; % 斯托克斯种子光强 max_iter 500; tol 1e-6; for iter 1:max_iter Ip_old Ip; Is_old Is; % 正向推进泵浦dIp/dz -gB*Ip*Is - alpha*Ip for k 2:Nz Ip(k) Ip(k-1) - dz*(gB*Ip(k-1)*Is(k-1) alpha*Ip(k-1)); end % 反向推进斯托克斯dIs/dz -gB*Ip*Is alpha*Is for k Nz-1:-1:1 Is(k) Is(k1) dz*(gB*Ip(k1)*Is(k1) - alpha*Is(k1)); end % 收敛判据两次迭代的最大绝对变化小于 tol if max(abs(Ip - Ip_old)) tol max(abs(Is - Is_old)) tol fprintf(第 %d 次迭代收敛\n, iter); break; end end正推泵浦时用 k-1 处的值这是显式欧拉格式简单但要求 dz 不能太大否则数值上会发散。反推斯托克斯时因为它是反向波更新方向从 zL 往 z0 走公式里 g_B 项为正、α 项为负对应它边传播边被放大、同时被损耗。注意收敛判据用的是绝对差值如果泵浦功率接近 SBS 阈值迭代可能在两个状态之间振荡。这种情况我一般会加亚松弛% 亚松弛每次只取一半新值抑制迭代振荡 relax 0.5; for k 2:Nz Ip_new Ip(k-1) - dz*(gB*Ip(k-1)*Is(k-1) alpha*Ip(k-1)); Ip(k) relax*Ip_new (1-relax)*Ip(k-1); end斯托克斯方向同样处理。亚松弛的代价是收敛变慢但稳定性明显改善尤其在 g_B·I_p·I_s 项接近阈值时。后面避坑章我会讲一个具体的振荡案例。3. 跑通 SBSa.m参数清单、扫频循环与布里渊频移的读数方法3.1 运行前的物理参数清单与单位换算SBSa.m 的所有物理量集中在脚本头部定义改参数不用动主循环。这张表是我每次跑脚本前必查的清单单位错一位结果就全错参数符号典型值单位说明泵浦波长λ_p1550nm通信波段常用有效折射率n1.45-掺锗光纤略高声速V_a5960m/s纯石英值掺锗后会降布里渊频移ν_B约 11GHz由 n、V_a、λ 算出增益线宽Δν_B3050MHz洛伦兹线型半高全宽峰值增益g_B04×10⁻¹¹m/W与纤芯掺杂有关有效模场面积Aeff80μm²G.652 典型值损耗系数α0.2dB/km换算线性损耗要除以 4.343泵浦功率P_p120mW超过阈值会出现耗尽种子功率P_s10⁻⁷10⁻⁶W等效自发散射噪声这里最容易翻车的换算是损耗0.2 dB/km 换算成线性吸收系数 α公式是 α 0.2/(4.343×1000)单位 1/m。直接拿 0.2 进方程会让损耗项大几千倍泵浦光在几米内就被算没了。提示0.2 dB/km 换算线性吸收系数时除以 4.343这是光纤损耗仿真里最常见的单位陷阱SBSa.m 脚本内已处理但你自己另写脚本时别漏。3.2 主循环之外的扫频模块三波耦合方程解出来的是某个固定失谐频率下的功率分布。要得到完整布里渊增益谱需要对失谐频率扫点在每个频点按洛伦兹线型算一个 g_B重跑一遍耦合方程记录 z0 处输出的斯托克斯功率。SBSa.m 里这个扫频循环大致长这样% 扫频范围布里渊频移附近 ±0.5 GHz nu_scan linspace(nu_B - 0.5e9, nu_B 0.5e9, 501); Ps_out zeros(size(nu_scan)); for m 1:length(nu_scan) % 按洛伦兹线型计算当前失谐频率下的增益系数 gB gB0 * (dnu_B/2)^2 / ((nu_scan(m) - nu_B)^2 (dnu_B/2)^2); % 重新初始化并迭代求解调用封装好的求解器 [Ip, Is] solve_sbs(gB, L, alpha, Pp_in, Ps_in, Aeff, Nz); % 记录 z0 处斯托克斯输出功率 Ps_out(m) Is(1) * Aeff; end % 画增益谱 figure; plot((nu_scan - nu_B)/1e9, 10*log10(Ps_out/max(Ps_out)), LineWidth, 1.5); xlabel(失谐频率 [GHz]); ylabel(相对斯托克斯功率 [dB]); title(SBS 布里渊增益谱); grid on;这里把求解器封装成函数形式输入当前频点 g_B 和固定参数输出两条功率分布。501 个频点每个都要迭代几百次MATLAB 跑起来几十秒到几分钟取决于 Nz 和收敛速度。如果只关心峰值频移201 个点足够要做高精度洛伦兹拟合再加密到 1001 个点。3.3 从增益谱读出布里渊频移洛伦兹拟合扫完频画出增益谱峰值位置就是布里渊频移 ν_B。但人眼读图精度有限增益谱顶部偏平直接用 max 找峰值会损失精度。工程上我一般用洛伦兹拟合% 洛伦兹拟合y A/((x-x0)^2 (gamma/2)^2) B fitfun (p, x) p(1)./((x - p(2)).^2 (p(3)/2)^2) p(4); p0 [max(Ps_out), nu_B, 35e6, min(Ps_out)]; p lsqcurvefit(fitfun, p0, nu_scan, Ps_out); fprintf(拟合布里渊频移: %.6f GHz\n, p(2)/1e9); fprintf(拟合线宽 FWHM: %.2f MHz\n, p(3)/1e6);拟合得到的 p(2) 是亚兆赫兹精度的布里渊频移p(3) 是线宽。做光纤传感标定时温度每变化 1℃ 频移移动约 1.18 MHz应变每变化 1 με 频移移动约 0.05 MHz所以拟合精度直接决定你能分辨多小的温度和应变变化。这也是我强调别用 max 找峰值的原因——max 的离散误差很容易到几个兆赫兹比温度系数还大。实际项目中我还会同时打印峰值和 FWHM 两个量线宽窄说明声波阻尼小、相干性好适合高精度传感线宽被泵浦功率展宽时说明泵浦已靠近阈值系统信噪比会恶化。4. 从仿真到分布式光纤传感温度应变标定与阈值自检4.1 布里渊频移的温度与应变系数石油钻探里的井筒温度剖面、大跨度桥梁的结构健康监测、地震台网的应变场重建靠的都是光纤沿线几千个点的频移信息。SBS 能做成全分布式传感器核心在于布里渊频移 ν_B 对温度和应变敏感。工程上用的线性模型是Δν_B C_T·ΔT C_ε·ΔεC_T 在 1550 nm 窗口大约 1.01.2 MHz/℃C_ε 大约 0.0480.052 MHz/με。这个线性关系在常规工作范围内足够准确但温度跨度超过 100℃ 或应变超过几千 με 时就要考虑二阶项。物理机制上温度升高时声速下降、折射率微变共同导致 ν_B 下移应变拉伸光纤时声速和折射率也变。严格说这两个系数不是独立常数工程标定时会分段拟合比如 -2060℃ 用一条直线60150℃ 用另一条交点处误差最小。仿真里做这件事很方便因为你可以精确控制每个温度点的声速值而实验里很难把光纤段均匀加热。4.2 在 SBSa.m 上做参数扫描标定系数怎么算出来我一般会在 SBSa.m 外面再套一层参数扫描脚本模拟不同温度下的增益谱从拟合频移里反推 C_T% 温度扫描模拟 -20℃ 到 80℃ 的频移响应 T_range -20:10:80; nuB_fit zeros(size(T_range)); CT 1.18e6; % 经验初值单位 Hz/℃ for t 1:length(T_range) % 温度改变量折算成布里渊频移偏移 nu_B_sim nu_B CT * (T_range(t) - 25); % 扫频并洛伦兹拟合复用第三章流程 nuB_fit(t) fit_brillouin_peak(nu_B_sim); end % 线性拟合得到仿真标定系数 p_lin polyfit(T_range, (nuB_fit - nu_B)/1e6, 1); fprintf(仿真标定系数 C_T %.3f MHz/℃\n, p_lin(1));这里的设计取舍值得说我用经验值 CT 偏移 ν_B再通过完整仿真流程把频移拟合回来实际上是在验证输入偏移→仿真→拟合输出这条链路有没有额外误差。如果拟合出来的 C_T 和输入经验值偏差超过 0.01 MHz/℃说明扫频步长或拟合初值有问题先修仿真再谈标定。这类自检逻辑在真实传感系统里同样适用——先注入已知扰动看系统能不能准确还原。应变扫描改法一模一样把 T_range 换成应变数组、CT 换成 C_ε 就行。仿真标定的意义在于实验里很难精确控制光纤的均匀应变仿真里每个参数都可控能干净地分离温度和应变两个变量的贡献这在真实实验里恰好是最难的部分。4.3 SBS 阈值功率的仿真界定传感系统设计里另一个关键参数是 SBS 阈值——泵浦功率超过它之后大量能量反向散射前向信号被严重损耗。仿真里看阈值很方便固定其他参数扫描泵浦功率记录 z0 处斯托克斯输出功率功率曲线由线性增长转为指数增长的拐点就是阈值。经验公式是 P_th ≈ 20·A_eff/(g_B·L_eff)其中 L_eff (1 - e^(-αL))/α 是有效相互作用长度。对 20 km 光纤、Aeff80 μm²、g_B4×10⁻¹¹ m/WL_eff 约 19.6 km算出来 P_th 大约 2 mW。这个数值和仿真扫描结果对得上说明模型设置没问题——我习惯把这事当成 SBSa.m 的出厂自检。扫描泵浦功率还有个附带收益不同光纤长度下阈值不同。比如你要在 40 km 光纤上做传感阈值会比 20 km 低一半左右泵浦功率选高一点增益大但超过阈值会引发串扰仿真能帮你找到折中值。5. SBS 仿真避坑指南五个实测翻车现象与排查方法5.1 增益谱不收敛迭代在振荡现象主循环迭代几百次不收敛泵浦功率曲线在两个形状之间来回跳tol 一直压不下来。原因泵浦功率接近或超过 SBS 阈值时耦合方程进入高增益区固定点迭代的收敛半径很小显式欧拉格式的数值稳定性也不够两者叠加就振荡。还有个隐蔽因素收敛判据用绝对容差种子光功率本身只有 10⁻⁷ W 量级绝对容差 10⁻⁶ 根本压不住相对误差。解决两步走。先把网格加密一倍看是否改善不行就加亚松弛relax 取 0.5 甚至 0.3。同时把收敛判据改成相对变化max(abs(Ip-Ip_old)./max(abs(Ip))) tol。如果还没改善检查泵浦功率是否设到了几十毫瓦那已经远离 SBS 仿真通常的工作区间。5.2 布里渊频移和实测对不上现象脚本里用 n1.45、V_a5960 算出 ν_B11.14 GHz实测光纤只有 10.85 GHz偏差近 300 MHz。原因ν_B 对纤芯材料极敏感。掺锗会降低声速也微调折射率G.652 光纤的有效声速实际在 53005600 m/s 范围不是纯石英的 5960。解决不要迷信理论公式的出厂值。跑仿真前查所用光纤型号的手册或者用窄线宽激光器实测一次布里渊谱反推等效声速填进脚本。Corning SMF-28 的有效声速约 5840 m/s算出来 ν_B≈10.92 GHz比纯石英的 11.14 接近实测得多。做一次校准后面所有频移分析才有意义。5.3 泵浦功率不低仿真里却没有增益现象P_p5 mW 时理论上早该有可观的布里渊增益但输出的斯托克斯功率几乎等于输入的种子功率增益谱是平的。原因常见有三个——种子功率设太低且迭代容差太大小信号增益淹没在数值误差里或者损耗系数没换算成 1/mα 大了上千倍泵浦在初始段就被耗光或者扫频起点正好落在洛伦兹线型的边翼。解决种子功率设到 10⁻⁷ W 量级代表等效自发噪声收敛容差收紧到 10⁻⁸确认为 0.2 dB/km 换算出线性吸收系数扫频中心先对准理论 ν_B。快速诊断方法是先跑一波 0.1 mW 泵浦的小功率对照此时 SBS 增益可忽略泵浦曲线应该只呈现线性损耗——如果这条曲线都不对问题一定在 α 换算先修损耗再看增益。5.4 泵浦耗尽的位置和量级不对现象功率曲线显示泵浦在 z5 km 处就掉到接近零而理论上 20 km 光纤的泵浦衰减主要由线性损耗主导SBS 耗尽应该集中在增益最强的一段。原因把 g_B0 当常数用在整个扫频范围内谱线中心失谐处增益被严重高估另外如果斯托克斯种子功率设成泵浦的十分之一相当于人为制造强反向注入耗尽自然提前。解决每个扫频点必须按洛伦兹线型重新计算 g_B(ν)只有失谐为零时才能用 g_B0。种子功率保持比泵浦小 45 个数量级那是更接近物理实际的注入条件。5.5 换网格数结果大变现象Nz 从 1000 改成 500或从 1000 改成 4000布里渊频移不变但功率曲线明显不同尤其泵浦耗尽深度变化很大。原因耦合方程里增益项和损耗项的特征尺度不同显式欧拉格式的局部截断误差是 O(dz)dz 太大时每个步长的误差累积起来足以改变积分轨迹。损耗项还要求 dz·α 远小于 1α4.6×10⁻⁵ 1/m 时 dz10 m 没问题但 dz 拉到 100 m 就开始劣化。解决养成固定习惯——每次改完参数做一次网格收敛性检查Nz 分别取 1000、2000、4000 跑一遍功率曲线重合才算合格。这类检查花不了几分钟但能挡住最隐蔽的数值误差。6. 进阶把 SBSa.m 扩展成 BOTDA 分布式传感仿真骨架连续波 SBS 仿真跑通后往全分布式传感再走一步就是 BOTDA布里渊光时域分析。改动其实不大把连续泵浦换成脉宽 τ 的脉冲脉冲在光纤里传播时每个位置产生布里渊增益反向探测光在不同时刻接收到的增益对应不同位置。位置映射是 z c·t/(2n)脉宽决定空间分辨率典型 τ10 ns 对应约 1 m 分辨率。仿真上把稳态迭代改成时间-空间二维推进脉冲沿 z 传播每个时间步在该位置处计算 SBS 增益累积到反向探测光的强度里。脚本结构可以在 SBSa.m 外面套一层时间循环核心的耦合方程求解函数不用动。这个阶段我还会额外做两件事一是给增益谱加一个高斯噪声项模拟探测器散粒噪声这样拟合算法在噪声下的精度表现才真实二是把洛伦兹拟合函数从主脚本里独立出来因为 BOTDA 每个位置都要拟合一条谱线成千上万次拟合必须向量化否则 MATLAB 会慢到没法用。从那以后我每次跑 SBS 仿真都强制走一遍完整流程先网格收敛检查再阈值自检最后用已知的 C_T 和 C_ε 做闭环验证三关都过了才敢拿仿真结果去指导实验。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
