
简介一份面向光纤传感与非线性光学领域的 MATLAB 仿真脚本围绕受激布里渊散射三波耦合方程编写适合研究生、高年级本科生及工程师用于科研入门、课程设计和算法调试。程序涵盖泵浦光、信号光与布里渊散射光三类光波的耦合建模支持设置泵浦功率、光纤长度、频率失谐等参数求解稳态或瞬态三波耦合方程组输出散射强度、增益谱等结果从而直观理解 SBS 频移与阈值特性。压缩包共一个文件为独立脚本整体大小约2KB轻量简洁无需额外数据文件在 MATLAB 中即可直接运行和修改。已有 609 人浏览学习。借助该脚本可观察三波耦合过程中光功率的转移规律对比不同参数组合下的布里渊散射变化并可在其基础上扩展温度、应变等分布式传感参数仿真为非线性光纤光学建模和光纤传感系统设计提供基础工具。1. 受激布里渊散射为什么绕不开三波耦合做光纤传感的人早晚会撞上受激布里渊散射SBS。它在应变和温度分布式测量里是核心物理效应但在仿真和数据处理时你看到的不是一句频移与应变线性相关那么简单而是一组三波耦合方程。泵浦光、斯托克斯光、声波这三者之间不是独立的传播关系而是互相喂能量、彼此制约的耦合系统。只看稳态增益系数、不碰方程最多能算个阈值想做分布式传感的定量反演或者评估泵浦耗尽的影响就必须把方程摆出来解。这也是为什么像 SBSa 这类数据归档出现在项目里时第一步不是直接画谱线而是先确认你手上有没有一套能对上实验条件的耦合方程参数。这篇文章就沿着SBSa.rar → scattering fiber → 三波耦合方程 → 光纤传感这条线把从方程到拟合的路径走一遍适合已经做过 BOTDA/BOTDR 但想在理论上往深走一点的人。2. 从散射光纤到三波耦合物理图像的建立2.1 为什么是三波而不是两波普通散射只涉及光与介质但 SBS 里有一个关键差异介质本身变成了一个动态参与者。泵浦光通过电致伸缩在光纤内产生一个以声速移动的密度波这个声波又反过来通过弹光效应调制折射率让泵浦光布拉格衍射成斯托克斯光。斯托克斯光与泵浦光的干涉场继续驱动声波形成一个反馈闭环。所以方程里至少要有三个场泵浦波、斯托克斯波、声波。如果只写光强耦合像dIp/dz -gB Ip Is那已经是对声波做绝热消除之后的结果工程上够用但解释不了声子寿命、增益谱的洛伦兹轮廓、以及脉冲光的瞬态行为。式子里声波不是静态光栅而是一个带有阻尼的受迫振动系统。声场的弛豫时间决定了 SBS 增益谱的带宽这个带宽反过来限制了 BOTDA 里脉冲宽度与频移分辨率的关系。所以做光纤传感时三波耦合不是故弄玄虚而是让你在频移之外还能推算出声子寿命和增益系数这些真正决定传感距离与精度的参数。2.2 慢变振幅近似下的三波耦合方程常见的三波耦合方程写法很多但都基于慢变振幅近似忽略二阶导数项同时假设泵浦与斯托克斯都是准单色波。下面这组是在光纤传感里用得最多的形式z 沿光纤长度t 是时间声场用复数振幅表示∂Ep/∂z (1/vg) ∂Ep/∂t i κ Es Q - (α/2) Ep ∂Es/∂z - (1/vg) ∂Es/∂t -i κ Ep Q* (α/2) Es ∂Q/∂t Γ Q i Λ Ep Es*其中 Ep 和 Es 分别是泵浦和斯托克斯缓变包络Q 是声波振幅Γ 是声子阻尼率数值上等于 1/声子寿命 τBκ 和 Λ 是耦合系数。不同文献的 κ 与 Λ 可能差一个共轭或符号但只要最后推出的增益系数 gB 一致就行。这个方程组要解的是三场沿光纤的演化泵浦从 z0 注入斯托克斯通常从 zL 反向注入或由自发噪声产生声场在时间上随时跟随但存在惯性。工程计算里最常用的简化是把时间项去掉直接看稳态。这时声场可以消掉方程退化成两个光强方程。但如果做 BOTDA 掺铒光纤放大器的瞬态分析时间项就得保留。2.3 参数怎么定gB、声子寿命和有效面积解方程之前必须先定参数。下面的表给了 1550 nm 单模光纤的典型量级实际值随光纤掺杂、纤芯材料和温度变化查 SBSa 数据时可以用这些做初始猜测参数符号典型值说明布里渊频移νB约 10.8 GHz 1550nm随应变与温度近似线性变化增益系数gB1~5 × 10^-11 m/W与光纤掺杂和声波限制有关声子寿命τB5~20 ns决定增益谱半宽增益谱半高全宽ΔνB15~70 MHz与 τB 成反比损耗系数α0.2~0.3 dB/km用线性单位代入方程有效面积Aeff80 μm^2 量级影响泵浦光强这几个参数的关联很紧密ΔνB 1/(π τB)gB 又和声波模式重叠度挂钩。多数三波耦合方程实现里并不会直接面对 Aeff而是把它折进耦合系数里。你要知道的是当你从某个光纤传感实验里反推出 gB 比表里大很多时先检查是不是把光强换算成了功率、把连续光当成了聚焦后的局域光强。3. 三波耦合方程的求解稳态解、阈值与数值实现3.1 稳态方程组的无量纲化在光纤传感的场景里脉冲宽度远大于声子寿命时可以忽略时间微分项使用稳态近似。把声场消去后得到泵浦与斯托克斯的光强耦合方程组dIp/dz -gB Ip Is - α Ip dIs/dz -gB Ip Is α Is注意方向约定泵浦沿 z斯托克斯沿 -z。所以在数值积分时Is 的初始条件在 zL 处给出。直接对这个方程组积分会碰上一个数值问题如果 Ip 和 Is 数量级差十个量级乘法容易溢出。常见做法是先无量纲化把光强归一化到各自的参考值比如用泵浦输入功率对应的光强 I0p Ip / I0 s Is / I0 Z z / L代入后得到dp/dZ -γ p s - αL p ds/dZ -γ p s αL s其中 γ gB I0 LαL α L。这样就把乘积项的尺度压下来数值积分会稳定很多。用 Any 语言做 RK4 时边界条件是 p(0)1s(1)给定的输入斯托克斯光强。3.2 Python 解稳态三波耦合的最小实现下面的 Python 例子用scipy.integrate.solve_bvp处理两点边界值问题因为泵浦和斯托克斯不是一个方向不允许直接单次积分。代码里明确了边界条件和参数换算。import numpy as np from scipy.integrate import solve_bvp # 参数设置L为光纤长度(m)gB为布里渊增益(m/W) # I0为泵浦输入光强(W/m^2)Is0为斯托克斯输入光强(W/m^2) L 1000.0 gB 3e-11 alpha 4.6e-5 # 0.2 dB/km 换算成线性单位/m I0 1e9 # 1W 有效面积 1e-9 m^2 得到的光强 Is0 1e3 # 微瓦量级作为种子光输入 gamma gB * I0 * L alphaL alpha * L r0 Is0 / I0 # 无量纲化斯托克斯初值 def fun(z, y): p, s y dpdz -gamma * p * s - alphaL * p dsdz -gamma * p * s alphaL * s return np.vstack((dpdz, dsdz)) def bc(ya, yb): # 泵浦在z0处为1斯托克斯在z1处为r0 return np.array([ya[0] - 1.0, yb[1] - r0]) z np.linspace(0, 1, 500) y_init np.zeros((2, z.size)) y_init[0] 1.0 y_init[1] r0 sol solve_bvp(fun, bc, z, y_init) z_phys sol.x * L p sol.y[0] s sol.y[1] Ip p * I0 Is s * I0 # 输出泵浦耗尽信息 print(输出端泵浦光强占比: {:.2f}.format(p[-1])) print(输入端斯托克斯光强: {:.2e} W/m^2.format(Is[0]))fun里直接写无量纲方程组bc把边界条件收紧ya[0]-1.0表示泵浦在 z0 归一化后等于 1yb[1]-r0表示斯托克斯在 z1 处等于输入种子。solve_bvp的初值猜测需要大致合理否则可能不收敛所以代码里把泵浦初猜成常数 1、斯托克斯初猜成常数r0。如果光纤长度特别长或增益特别大建议把z的网格加密到 1000 点以上否则在 L 端附近出现边界层时解会振荡。3.3 阈值条件与声子寿命的关系当斯托克斯没有种子光、只靠自发噪声起振时三波耦合稳态方程没办法给出噪声驱动的解。工程上用的是布里渊阈值判定式gB P_th L_eff / A_eff ≈ 21其中 L_eff (1 - exp(-αL)) / α。这个条件是泵浦在光纤长度上对斯托克斯的净增益与腔反馈达到平衡的经验值对单模光纤通常取 20~25但具体到散射光纤会因为模式分布不同略有偏移。你可以用上面的代码做一次验证把 Is0 设置成极小值例如 1e-12 量级然后看输入端斯托克斯光强是否随 P0 出现指数倍增到接近泵浦量级。如果出现陡峭上升那个拐点位置就是阈值。阈值和声子寿命在这里体现为时间尺度的关系只有在声子寿命远小于斯托克斯与泵浦的拍周期时稳态方程才成立。拿 BOTDA 为例脉冲宽度一般要求大于声子寿命的数倍否则增益谱会展宽、频移峰位失真。所以读 SBSa 数据时先看它的时域分辨率与脉冲宽度如果脉冲宽度比 5 ns 短三波耦合方程就必须保留 ∂Q/∂t 项。4. 用 SBSa.rar 里的数据做三波耦合参数拟合4.1 拿到 SBSa.rar 后第一步做什么这类压缩包通常包含的是布里渊散射实验的频域谱线可能按不同泵浦功率、不同光纤段组织。不要假设里面已经是干净的频移表常见的处理流程是解包、目录梳理、数据抽样、拟合、反演参数。先做目录树检查避免脚本把二进制文件和文本文件混着读mkdir sbsa unzip SBSa.rar -d sbsa find sbsa -type f | head -50unzip对 .rar 不适用时改用unar或7z x。看到文件命名规律后用 Python 批量读取。如果里面是 CSV 或 MATLAB 的 .mat 文件直接np.loadtxt或scipy.io.loadmat如果是二进制 raw 数据那就需要先找读取说明或按固定字节长度试解析。4.2 从频谱数据拟合布里渊频移与增益谱宽SBS 的增益谱在稳态下是洛伦兹形状这是三波耦合方程对连续波做洛伦兹解的直接结果。所以拟合函数用洛伦兹而不是高斯如果数据明显不对称再考虑卷积声子线型与脉冲光谱。下面的代码用curve_fit从一段频率扫描数据中提取中心频率nu_B、峰值增益g0和半高全宽df_Bimport numpy as np from scipy.optimize import curve_fit def lorentz(f, nu0, g0, df): return g0 * (df / 2)**2 / ((f - nu0)**2 (df / 2)**2) # 假设读取了两个数组freq(MHz), gain(归一化) # 实际使用前先对 gain 做基线扣除 freq np.linspace(10700, 10900, 2001) gain lorentz(freq, 10798.5, 1.0, 35.0) 0.002 * np.random.randn(freq.size) p0 [np.mean(freq), 0.8, 30.0] bounds ([10700, 0, 5], [10900, 5, 80]) popt, pcov curve_fit(lorentz, freq, gain, p0p0, boundsbounds) nu_B, g_peak, df_B popt tau_B 1.0 / (np.pi * df_B * 1e6) * 1e9 # 声子寿命(ns) print(布里渊频移: {:.2f} MHz.format(nu_B)) print(峰值增益: {:.4f}.format(g_peak)) print(半高全宽: {:.2f} MHz.format(df_B)) print(声子寿命: {:.2f} ns.format(tau_B))注意df_B与物理线宽的换算如果数据频率单位是 MHzcurve_fit返回的df_B直接以 MHz 为单位。在把半高全宽换算成声子寿命时要用1/(π Δν)如果改成 GBps 单位后算出来的寿命落在 20 ns 以上说明数据可能来自脉冲展宽后的表观线宽不是真实的声子寿命。这是三波耦合方程与实验数据对照时的经典陷阱。4.3 拟合结果的坑基线与泵浦耗尽用洛伦兹拟合 SBS 谱最大的隐患是数据没有做基线扣除。受激布里渊散射实验里光电探测器上总会叠加 Rayleigh 散射、放大自发辐射或干涉仪泄漏的基底信号。如果基底不是常数而是缓慢随频率起伏直接拟合的nu0会偏移几个 MHz。做法是用频谱两端各 10% 的窗口算线性基线先减掉再拟合。泵浦耗尽会表现为主峰值下降或频谱不对称。在三波耦合方程里当斯托克斯光强增长到与泵浦可比时泵浦沿 z 方向被显著消耗此时洛伦兹拟合的峰值增益偏低频移反而不影响太大。如果你发现某个功率点拟合出的g0随泵浦功率不是线性增长而是趋于饱和说明实验数据已经进入了强耦合区。这时应该改用第三章的稳态方程做正向仿真来拟合整个谱线而不是单个洛伦兹峰。5. 验证三波耦合结果的两个快速技巧5.1 用阈值和增益谱宽互相校验拟合完 SBSa 数据得到的τB可以直接用来校验实验条件。如果实验中脉冲宽度为 T_pulse增益谱的等效宽度会从稳态的ΔνB展宽为约sqrt((ΔνB)^2 (1/T_pulse)^2)。你可以把拟合出的表观线宽代入这个式子反推出真实的ΔνB再和泵浦功率阈值做一次闭环验证。数值上很多散射光纤在这两个独立测量下能得到一致的声子寿命不一致通常说明三波耦合方程里忽略了声波的空间走离walk-off效应。5.2 用频移-应变标定的正反算校验数据切割最后一个小技巧是用频移对应变/温度系数做自洽检查。1550 nm 单模光纤的布里渊频移温度系数约 1 MHz/°C应变系数约 0.05 MHz/με。从 SBSa 谱线拟合出的频移变化量除以对应环境变化量如果偏离这两个量级太多先怀疑数据切割位置是否错位再怀疑洛伦兹拟合是否被旁瓣干扰。拿一组已知温度差的剖面数据做一次线性回归截距就是参考频移斜率就是温度系数拿它与标称值对照比单独看一个单点谱可靠得多。这个步骤做完你的三波耦合参数才算真正和实验对上了。本文还有配套的精品资源点击获取