基于MATLAB的受激布里渊散射放大仿真与参数调优 简介这份压缩包提供的是一个基于MATLAB平台的受激布里渊散射放大模拟程序面向光学工程、光纤通信以及非线性光学领域的学生与科研人员主要用于研究布里渊放大和散射过程中的物理现象与规律。资源共包含一个文件即放大程序点m文件整个压缩包占用空间仅约一千字节代码非常精简便于直接运行和调整参数。目前已有五百四十三人学习使用在光纤通信系统、光学放大器设计以及光信号处理等场景中具有实际参考价值。该程序基于受激布里渊散射的物理模型通过数值方法求解光波与声波相互作用的波动方程模拟泵浦光激发声波并产生斯托克斯光的放大过程用户可以观察泵浦功率、光纤长度、信号频率等关键参数对增益效果的影响并借助结果绘图功能直观理解布里渊散射的机理为相关实验方案设计和系统性能优化提供有力支持。1. 受激布里渊散射的MATLAB仿真到底在算什么假设你手头正在搭一套基于布里渊放大器的光纤传感系统输入信号只有微瓦量级探测器灵敏度又卡得很紧又或者你在做非线性光学的课程项目想复现教科书里那条经典的受激布里渊增益曲线。这个压缩包里的放大程序.m要干的事情很直接把泵浦光和反向传播的斯托克斯信号在光纤里的功率演化方程以数值方式从头到尾求解一遍得到沿线功率分布、放大增益、泵浦耗尽点这类设计参数。它不依赖COMSOL或者Photonics集成工具纯MATLAB脚本就能跑适合作为SBS放大链路预研和参数扫描的起点。下面按物理模型、数值求解、参数调优三层拆开讲并给出可以直接改参数运行的核心代码。2. 增益从哪里来SBS三波耦合方程与边界条件2.1 光子-声子散射的频移关系受激布里渊散射可以理解为三个波之间的参量过程泵浦光子、反向传播的斯托克斯光子和声学声子。泵浦光通过电致伸缩效应在光纤内激励起一个以声速传播的密度波这个密度波相当于一个运动的布拉格光栅把后向的泵浦光散射成频率更低的光。能量守恒给出 ω_s ω_p - Ω动量守恒给出 k_s k_p - q其中Ω是声波角频率。换算成频率布里渊频移近似为ν_B 2 n ν_p v_a / c其中n是纤芯折射率v_a是纵向声速ν_p是泵浦光频率。在1550 nm窗口普通石英单模光纤的典型值约为11 GHz这个数值直接决定了放大器和传感系统的失谐量设置。声波在传播过程中会衰减其寿命通常只有几个纳秒对应增益谱是洛伦兹型半高全宽典型为30 MHz左右。这也是为什么仿真里不能把布里渊增益当成常数。在窄带信号条件下信号光和泵浦光的频差只要偏离ν_B几百MHz增益就会掉到峰值的一半以下。实际工程中泵浦激光器与信号激光器的相对频率稳定度直接影响增益波动所以后边的参数扫描那章会专门讨论失谐量怎么扫。2.2 稳态强度耦合方程的适用边界在绝大多数SBS放大预研场景里脉冲宽度远大于声子寿命且信号频谱远窄于布里渊增益谱这时候可以不写包含时间项的瞬态方程直接用稳态强度耦合方程dp_p/dz - (ω_p/ω_s) (g_B / Aeff) P_p P_s - α P_pdP_s/dz - (g_B / Aeff) P_p P_s α P_s这里的z是泵浦正向传播方向P_s沿反向传播所以第二个方程右边是负的耦合项。g_B是布里渊峰值增益系数典型值48e-11 m/W和纤芯成分、掺杂浓度有关。Aeff是有效模场面积常见取值80 μm²换算成平方米是80e-12。α是功率衰减系数注意单位必须用每米而不是每千米。这套方程忽略了相位失配和偏振不匹配适用于普通单模光纤中偏振随机但统计均匀的情况。如果信号功率很强导致泵浦被大量消耗方程左侧的耦合项会自然体现增益压缩效应这正是后面要看重分析的饱和行为。对于更高精度的相位敏感仿真需要在方程里引入复振幅和相位失配项但从放大器增益设计的角度看强度方程已经够用。2.3 为什么这是个边值问题而不是初值问题从输入条件上看泵浦从z0注入信号从zL注入两者在光纤中相向而行。因此P_p的初值在z0已知而P_s的初值在zL已知。常规的ode45只能从单侧积分这里就必须采用打靶法或者MATLAB的bvp系列求解器。打靶法的思路是先猜一个P_s(0)正向积分到zL看P_s(L)是否等于已知的注入功率不相等就修正猜测值。这个做法在小信号情况下收敛很快但一旦泵浦耗尽明显迭代矩阵容易病态。bvp5c则是把整个区间离散成网格用有限差分配合牛顿迭代同时满足边界条件工程上更省心也更容易控制精度。实际写MATLAB时我一般用bvp5c因为它的网格自适应控制比手写打靶法稳。需要注意边界条件的写法泵浦强度在z0固定为Pp0/Aeff信号强度在zL固定为PsL/Aeff。如果把两个边界条件都放在同一端程序会收敛到一个数学上存在但物理上无意义的解仿真出来的增益曲线和解析解对不上。3. 用 MATLAB 数值积分求解放大程序.m 的核心结构3.1 求解前的单位与参数归一化SBS方程里最容易出问题的是单位不统一。比如习惯上光纤损耗写0.2 dB/km而布里渊增益系数g_B是m/W长度是米这时候必须将0.2 dB/km转换成每米的功率衰减系数。0.2 dB/km除以4.343得到0.046 km^-1也就是4.6e-5 m^-1。很多从文献里抄来的参数Aeff是80 μm²写成m²是80e-12。如果直接把μm²代进去耦合系数会相差10^-6量级出来的增益曲线完全不对。另一个需要注意的点是泵浦与信号功率的单位。方程里P_p和P_s应使用同一功率单位且g_B/Aeff之后乘积的单位是1/(W·m)。功率用W面积用m²长度用m这样耦合项的量纲自然消掉。部分参考书使用光强I和光强增益系数g_B关系为IP/Aeff两者本质上是一回事。在放大程序.m里推荐在一个单独的参数初始化段完成所有常数定义再传递给求解器避免在主程序里到处改数字。下面的代码是把方程组直接写成bvp5c能接受的格式并给出完整的求解流程。3.2 用 bvp5c 求解反向SBS方程% sbs_bvp_demo.m % 正向泵浦 z0-L反向信号 zL-0 clear; close all; % 物理常数及光纤参数 c 3e8; n 1.45; va 5945; % 石英声速 m/s nu_p c / 1550e-9; nu_B 2 * n * va * nu_p / c; % 约 11 GHz Gamma 2 * pi * 30e6; % 声子衰减率 rad/s g0 5e-11; % 峰值增益 m/W dnu 0; % 信号-泵浦频差相对nu_B的失谐 Hz dB 1 / (1 (2*pi*dnu/Gamma)^2); % 洛伦兹失谐因子 gB g0 * dB; L 20e3; % 长度 m alpha 0.046e-3; % 损耗 1/m Aeff 80e-12; % 有效面积 m^2 Pp0 0.2; % 泵浦注入功率 W PsL 1e-5; % 信号注入功率 W r 1; % 频率比 omega_p/omega_s % 边界值问题的初值猜测 solinit bvpinit(linspace(0, L, 20), [Pp0/Aeff PsL/Aeff]); sol bvp5c((z,y) sbs_ode(z,y,gB,alpha,r), ... (ya,yb) sbs_bc(ya,yb,Pp0/Aeff,PsL/Aeff), solinit); z linspace(0, L, 200); y deval(sol, z); Pp y(1,:) * Aeff; % 泵浦功率 Ps y(2,:) * Aeff; % 信号功率 fprintf(信号输出功率%.4e W增益%.2f dB\n, ... Ps(end), 10*log10(Ps(end)/PsL)); function dydz sbs_ode(z, y, gB, alpha, r) % y(1)Ip 泵浦强度y(2)Is 信号强度 dydz zeros(2,1); dydz(1) -r*gB*y(1)*y(2) - alpha*y(1); dydz(2) gB*y(1)*y(2) alpha*y(2); end function res sbs_bc(ya, yb, Pp0, PsL) % 泵浦强度在 z0 固定信号强度在 zL 固定 res [ya(1) - Pp0; yb(2) - PsL]; end这段代码的核心是sbs_ode里的两个耦合项。第一行rgBy(1)y(2)表示泵浦由于SBS转移给信号而产生的损耗r是频率比一般近似为1第二行的gBy(1)*y(2)是信号获得增益。两个方程的符号相反正好满足能量守恒。边界函数sbs_bc里ya表示z0处的状态yb表示zL处的状态数组位置和y的定义顺序一致。这里信号在zL处被固定为PsL/Aeff所以程序会自动调整z0处的信号强度使其积分到L端时正好等于注入值。bvp5c很适合这种常系数耦合方程但要注意初始猜测会影响收敛。如果泵浦功率或光纤长度增大一个量级最好先用小信号解析解作为初值或者增加网格点数量。另外求解器的默认容差为1e-3严格仿真时需要显式设置否则增益曲线在高增益区域会出现抖动。3.3 参数表与工程取值下面是二氧化硅单模光纤在1550 nm附近的一组典型参数大部分公开文献都能对得上但不是某个固定版本参数符号典型值单位说明折射率n1.45无纤芯等效折射率声速v_a5945m/s纵向声学速度声子线宽Γ_B/2π30MHz增益谱半高全宽布里渊增益g_B5e-11m/W与光功率和掺杂有关有效面积Aeff80μm²换算为80e-12 m²损耗α0.046km^-10.2 dB/km布里渊频移ν_B11GHz1550 nm处表里最容易被误用的是g_B和Aeff。不少论文给的是g_B/Aeff组合后的系数比如g_B/Aeff约0.6 W^-1m^-1复制进代码时要先拆开还是直接使用取决于你选用的是功率方程还是强度方程。功率方程直接用g_B/Aeff合在一起的耦合系数更不容易错。4. 泵浦耗尽与增益饱和算出来的数据怎么读4.1 从指数增益到泵浦耗尽当信号很弱时泵浦沿光纤的损耗几乎只有线性衰减信号增长的解析表达式为Ps(L) ≈ Ps(0) exp(g_B Pp0 L_eff / Aeff)其中有效长度 L_eff (1 - exp(-αL)) / α当αL很大时约等于1/α。这个公式是校验数值仿真最有用的工具。但随着信号功率增大泵浦能量被明显消耗增益不再是指数关系出现了压缩现象。如果你的仿真里信号增益在提升泵浦功率后上升变缓甚至输出信号功率趋于饱和说明耦合项已经把泵浦拉低这正是物理上泵浦耗尽的表现。在放大程序.m的曲线图上泵浦功率沿z的分布会从指数衰减变成一种在z接近L处快速下跌的形状。这是因为泵浦在光纤末端被反向信号大量吸收能量转移给斯托克斯光。通常把泵浦功率跌落到初始值的一半的那个位置叫作耗尽点设计系统时希望耗尽点出现在光纤末端附近这样整段光纤提供的增益最充分。4.2 从仿真结果提取增益与转换效率运行上面的代码后Ps(end)是信号从光纤zL端输出的功率。放大器的净增益为G_dB 10 log10(Ps(end) / PsL)这里PsL是在zL注入的信号功率不是z0处的功率。如果结果出现负增益先看泵浦功率是否超过了阈值再看信号注入方向是否和代码里的边界条件一致。常见错误是把PsL当成输出直接用Ps(1)算增益那算出来的是z0处的内部信号强度不能反映放大器性能。功率转换效率定义为信号功率增加量与泵浦注入功率的比值eta (Ps(end) - PsL) / Pp0SBS放大器的能量转换效率受限于量子亏损理论上不超过ν_s/ν_p约为1实际上由于损耗和泵浦耗尽能到50%已经很不错。把eta也打印出来和增益一起画成泵浦功率Pp0的函数就能直观看到增益提升和效率回落的不同区间。4.3 失谐频率对增益谱的影响SBS增益谱的洛伦兹形状可以由dnu参数控制dnu是信号光频率与泵浦光频率之差再减去ν_B后的数值。把前面的代码包装成一个循环每个dnu调用一次bvp5c记录输出增益dnu_list -80e6:5e6:80e6; G_list zeros(size(dnu_list)); for k 1:numel(dnu_list) dnu dnu_list(k); delta_w 2*pi*dnu / Gamma; gB g0 / (1 delta_w^2); % 重新调用求解器这里省略重复参数设置 % G_list(k) 10*log10(Ps_end / PsL); end plot(dnu_list/1e6, G_list); xlabel(频率失谐 (MHz)); ylabel(增益 (dB));选择频差时要围绕光纤的实际ν_B展开。比如1550 nm普通光纤ν_B约11 GHz如果你的泵浦光源标称线宽1 MHz信号光源用可调谐激光器先扫描到增益最大位置再在这个基础上做精细步长扫描。这样能同时验证声子线宽Γ_B的取值是否合理。5. 调参小信号解析校验、容差设置和批量扫描5.1 用小信号解析解校验数值结果每次改了代码先用小信号条件下的解析增益做一次快速校验。取PsL1e-8 WPp00.1 WL20 km损耗α0.046e-3 m^-1则L_eff约18 km解析增益约为G_dB_analytic 10 log10( exp(gB Pp0 L_eff / Aeff) )代入gB5e-11、Pp00.1、L_eff18000、Aeff80e-12指数项约1.125增益约4.9 dB。因为这个场景里泵浦还远没有耗尽数值结果应该和它非常接近。如果bvp5c结果和解析值偏差超过0.5 dB检查边界条件是否把信号方向写反以及损耗单位是不是混用了km和m。5.2 设置bvp5c的容差默认RelTol是1e-3对增益波动敏感时设置opts bvpset(RelTol, 1e-6, AbsTol, 1e-8, Stats, on); sol bvp5c((z,y) sbs_ode(z,y,gB,alpha,r), ... (ya,yb) sbs_bc(ya,yb,Pp0/Aeff,PsL/Aeff), ... solinit, opts);网格自适应会承担大部分精度问题但若仍然不收敛先把linspace点数从20改成100再检查初值猜测的数量级是否离真实解太远。Stats参数会在求解结束后报告实际网格点和残差评估次数方便判断是否值得继续加密。5.3 用parfor批量扫描泵浦功率当需要扫描几十个Pp0值时每次调用bvp5c都是独立计算完全可以用parfor。要注意把参数定义放到循环体内或打包成函数避免广播变量把内存撑满。典型做法是写一个sbs_simulate(Pp0, dnu)函数返回增益和效率然后并行循环。设置并行池后一次参数扫描可以从几分钟压到几十秒最终把多条增益-功率曲线叠加在同一张图上观察不同光纤长度下的饱和趋势。本文还有配套的精品资源点击获取