简介这套MATLAB仿真资源以SBS受激布里渊散射三维谱为切入点面向光通信、光子学方向的学生与研究人员帮助零基础快速理解泵浦光、声子与斯托克斯光之间的非线性耦合机制。压缩包内仅含1个.m脚本大小仅1KB结构非常精简适合直接运行或在此基础上二次修改用于不同场景下的SBS谱形模拟。代码涵盖模型建立、光场传播模拟、声子动力学、SBS增益计算与频谱分析等关键步骤用户可灵活调整泵浦功率、光纤长度、折射率等参数动态观察三维受激布里渊散射谱的变化趋势从而深入掌握SBS效应的物理图像并为光纤通信系统中抑制SBS、优化传输性能提供直观的仿真支持。已有287人学习下载既可作为相关课程设计的参考代码也是科研入门阶段理解受激散射现象的实用工具。1. 从 sbs.rar 到 SBS 谱数据到物理量的解析链路把 sbs.rar 解压出来里面往往不是一张现成的图而是一组频率扫描曲线或二维矩阵。文件名里的 SBS 指受激布里渊散射Stimulated Brillouin Scattering谱线记录的是泵浦光在介质中驱动声学声子后产生的频移光强分布。做光纤传感的人拿到这套数据第一反应是找布里渊频移随温度和应变的变化做激光系统的人则会关心增益谱的线宽和阈值。下面这条链路从洛伦兹线型讲起再到 Python 提取频移和线宽最后落到分布式传感标定和数据验证。2. 受激布里渊散射谱的物理图像洛伦兹线型与关键参数2.1 三波耦合与声学声子SBS 谱从哪里来2.1.1 布里渊频移的守恒关系受激布里渊散射在背向散射几何里最典型一束泵浦光角频率 \omega_p与一束反向传播的斯托克斯光角频率 \omega_s \omega_p - \Omega通过电致伸缩效应在介质中共同驱动一个声学声子场角频率 \Omega。能量守恒和动量守恒同时成立时声子场被共振增强散射效率达到峰值。这个共振频率就是布里渊频移。对于体波声速为 V_a、折射率为 n 的介质背向布里渊频移可以写成 \nu_B 2 n V_a / \lambda_p其中 \lambda_p 是泵浦波长。代入 1550 nm 波段石英光纤的典型值n 约 1.45V_a 约 5945 m/s\nu_B 约 11.1 GHz。这就是为什么你在 SBS 谱上看到的峰总在 10~12 GHz 附近而不是在 THz 量级。2.1.2 稳态近似下的洛伦兹增益谱把三波耦合方程在泵浦未耗尽、声子场达到稳态的条件下做简化声子场振幅对驱动项正比于泵浦场与斯托克斯场共轭的乘积的响应是一个阻尼谐振子的频率响应。其模平方的谱形状为g(\nu) g_B \frac{(\Delta\nu_B / 2)^2}{(\nu - \nu_B)^2 (\Delta\nu_B / 2)^2}其中 \Delta\nu_B 是增益谱的半高全宽FWHMg_B 是峰值增益系数。也就是说理想 SBS 增益谱是单峰洛伦兹线型不是高斯线型。这个区别直接决定了后面的拟合函数形式任何用高斯去拟合 SBS 谱的做法至少在稳态受激散射场景下是缺少物理依据的。2.2 SBS 谱必须看准的四个参数拿到谱线之后真正要提取和关心的量就四个它们各自对应不同的物理过程和实验设定。参数符号1550 nm 石英光纤典型量级决定因素在测量中的角色布里渊频移\nu_B约 11 GHz折射率、声速、泵浦波长温度/应变测量的主读数增益线宽\Delta\nu_B20~50 MHz声子寿命、泵浦功率影响拟合精度与空间分辨率峰值增益系数g_B约 5 \times 10^{-11} m/W材料电致伸缩系数决定阈值功率和信号强度阈值功率P_{th}米级光纤数百 mW公里级可低于 10 mWg_B、有效面积、有效作用长度判断系统是否进入受激散射线宽 \Delta\nu_B 与声子寿命 T_B 有粗略关系 \Delta\nu_B \approx 1 / (\pi T_B)。常规单模光纤的声子寿命在 10 ns 量级所以线宽落在几十 MHz 是正常的。增益系数 g_B 则主要用于估算阈值P_{th} \approx 21 A_{eff} / (g_B L_{eff})A_{eff} 是有效模场面积L_{eff} 是有效作用长度。实验上如果发现谱峰异常宽先不要怀疑声子寿命变了多数情况是泵浦功率过高导致的增益展宽占主导。2.3 从 sbs.rar 解出来的数据怎么对应物理量这类数据文件最常见的组织方式有两种一种是每个文件单条谱线两列分别为频率和光强另一种是整个测量存成一个二维矩阵行是频率扫描点列是空间位置或时间序列。打开文件第一件事不是画图而是确认频率轴的单位是 MHz 还是 GHz以及扫描方向是从低频到高频还是反过来。如果数据是从光谱仪采的波长轴需要先做大范围换算。1550 nm 附近 1 nm 对应的频率间隔约 125 GHz换算关系 \Delta\nu \approx (c / \lambda^2) \Delta\lambdaimport numpy as np wl_nm 1549.5 # 中心波长单位 nm wl_step_nm 0.01 # 波长步进 freq_step_ghz 299792.458 * wl_step_nm / wl_nm**2 print(f波长步进对应的频率步进约 {freq_step_ghz:.3f} GHz)换算完成后还要检查扫描窗口和步进。我的经验是扫描窗口至少覆盖 3~5 个线宽也就是常规光纤里留 200~300 MHz频率步进要小于线宽的 1/5否则拟合出的 \nu_B 会出现量化误差。一个 3 GHz 扫描范围、1 MHz 步进的设定足够覆盖温度和应变同时变化时的频移漂移量也够把洛伦兹峰的肩部形状采清楚。3. 用 Python 从 SBS 谱提取布里渊频移和线宽的拟合流程3.1 加载 sbs.rar 解出的光谱矩阵解压后最常见的存储格式是 CSV 或纯文本表格。读取时把第一列当作频率轴其余列当作不同测量点的散射谱。可以用 pandas 一次性读入避免手写逐行解析import numpy as np import pandas as pd data pd.read_csv(sbs_spectrum.csv, headerNone).values freq data[:, 0] # 频率轴单位假定为 MHz spec_mat data[:, 1:] # 形状 (N_freq, N_point) spec spec_mat[:, 0] # 先取第一个点看形态 print(f频率范围 {freq[0]:.1f} ~ {freq[-1]:.1f} MHz, 点数 {freq.size})逻辑说明列数含义取决于仪器导出的约定读取前先看文件头部几行避免把第一列之外的标识列误当成数据。若频率轴单位是 GHz先统一乘 1000 转成 MHz否则后面拟合初值和边界都要跟着改。这一步如果发现矩阵里有 NaN需要按列处理不要整行丢弃否则会破坏频率轴对齐。3.2 单条谱线的洛伦兹拟合与参数初值有了频率轴和一条谱线直接用 scipy 的 curve_fit 做最小二乘拟合。拟合函数就是前文那个洛伦兹线型多一个常数项用来吸收探测器暗电流或 ASE 基底from scipy.optimize import curve_fit def lorentz(f, amp, fB, df, offset): SBS 增益谱洛伦兹线型df 为半高全宽 FWHM return offset amp * (df / 2) ** 2 / ((f - fB) ** 2 (df / 2) ** 2) # 初值用质心法粗略定位峰位避免 argmax 被脉冲噪声带偏 mask (freq freq[spec.argmax()] - 100) (freq freq[spec.argmax()] 100) fB0 np.sum(freq[mask] * spec[mask]) / np.sum(spec[mask]) p0 [spec.max() - spec.min(), fB0, 50.0, spec.min()] bounds ([0, freq[0], 1.0, -np.inf], [np.inf, freq[-1], 500.0, np.inf]) popt, pcov curve_fit(lorentz, freq, spec, p0p0, boundsbounds) perr np.sqrt(np.diag(pcov)) print(f布里渊频移 {popt[1]:.2f} ± {perr[1]:.2f} MHz, 线宽 {popt[2]:.2f} MHz)参数说明amp 是峰值增益幅度fB 是布里渊频移df 是谱线 FWHMoffset 是光谱基底。边界设置的逻辑是df 下限取 1 MHz 避免除零和退化上限取 500 MHz 覆盖高泵浦展宽情形fB 限制在扫描窗口内防止拟合跳出物理范围。质心初值比 argmax 初值稳健因为 SBS 谱叠加噪声后最大值点会随机抖动而质心对整个峰形的加权更稳定。3.3 批量拟合多位置、多时刻的 SBS 谱分布式传感系统里一个采集周期会产生几百到几千条谱线。逐条调用 curve_fit 时要处理两类失败拟合不收敛以及收敛到窗口边界上的假峰。批量循环里我会同时记录残差均方根作为后续筛选依据fits [] for i in range(spec_mat.shape[1]): y spec_mat[:, i] if np.isnan(y).any(): fits.append((np.nan, np.nan, np.nan)) continue # 用前一条谱的结果限制当前初值保证沿光纤连续变化 fB0 freq[y.argmax()] if i 0 else fits[-1][0] p0 [y.max() - y.min(), fB0, 50.0, y.min()] try: popt, _ curve_fit(lorentz, freq, y, p0p0, boundsbounds, maxfev20000) residual y - lorentz(freq, *popt) rmse np.sqrt(np.mean(residual ** 2)) fits.append((popt[1], popt[2], rmse)) except RuntimeError: fits.append((np.nan, np.nan, np.nan))逻辑说明用上一个点的频移作为当前初值利用了光纤中温度应变连续变化这一先验能显著减少跳变。rmse 保留下来后以第一条干净谱的 rmse 为基准超过其 3 倍的点视为坏点。坏点常见于熔接点、断纤处或泵浦功率骤降的位置这些位置的谱可能不是完整洛伦兹形强行保留会污染后续的标定。3.4 拟合结果的物理解读与错误边界识别拟合完成不等于解调完成。频移和线宽要回到物理场景里做合理性检验相邻采样点的频移差在短段普通光纤中通常只有零点几 MHz 到几 MHz如果出现几十 MHz 的突变优先怀疑拟合窗口偏移或噪声尖峰而不是真的发生了剧烈应变。现象可能原因处理方式频移跳到扫描窗口边缘初值给错拟合收敛到边界改用质心初值或追踪前一位置结果线宽拟合结果接近边界上限泵浦功率过高或谱含双峰降低泵浦功率检查是否为双峰叠加残差明显大于噪声本底窗口截断、基底非线性扩大扫描窗口或加二次基底项频移曲线出现单点毛刺该位置 SNR 过低剔除该点后用相邻点插值4. 分布式传感中的 SBS 谱解调与温度应变标定4.1 BOTDA 与 BOTDR 的谱形态差异同样是布里渊散射受激和自发的谱数据差异很大处理流程的侧重点也不同。BOTDR 是单端自发散射信号弱需要大量平均BOTDA 是双端受激过程增益谱信噪比高但系统里同时存在泵浦和探测两束光。项目BOTDABOTDR物理过程受激放大/损耗自发散射 相干探测谱线特征洛伦兹增益峰或损耗谷洛伦兹峰叠加噪声幅度极小实验配置双端接入光纤单端即可测量平均次数需求较低通常需数千次以上典型空间分辨率米级米级BOTDA 的谱是增益峰拟合时峰值朝上如果系统工作在损耗谱模式探测器看到的反而是一个凹陷这时拟合函数可以取负幅度或对原始数据做翻转。BOTDR 的数据则要注意噪声基底是否被低通滤波压平滤波过重会把洛伦兹肩部削掉导致线宽偏窄。4.2 逐点解调频移并拼接成距离曲线沿光纤每个位置都有一条独立的谱解调目标是把整条光纤的 \nu_B(z) 曲线拼出来。逐个点独立拟合最容易出的问题是相邻点频移抖动因为每个点的 SNR 不同。用追踪式拟合能显著改善fit_list [] fb_prev None for i in range(spec_mat.shape[1]): y spec_mat[:, i] if fb_prev is None: fB0 freq[y.argmax()] else: fB0 fb_prev lo, hi fB0 - 200, fB0 200 # 只保留峰附近 400 MHz 窗口 m (freq lo) (freq hi) p0 [y[m].max() - y[m].min(), fB0, 50.0, y[m].min()] local_bounds ([0, lo, 1.0, -np.inf], [np.inf, hi, 500.0, np.inf]) popt, _ curve_fit(lorentz, freq[m], y[m], p0p0, boundslocal_bounds) fb_prev popt[1] fit_list.append(fb_prev)逻辑说明把每个点的拟合窗口限制在上一位置频移附近 400 MHz 内既排除了远处假峰的干扰也能容忍相邻点几十 MHz 的真实频移差对应上千度的温度变化或上万个微应变。窗口范围可以按光纤最大应变梯度折算高压输电线或海底缆场景适当放宽普通光缆 200~400 MHz 足够。拼接完成后频移曲线如果出现整体跳变先看是不是解包问题而不是传感事件。沿着光纤的 \nu_B 变化应该连续突变点要回到原始谱上去人工确认峰形。4.3 从 MHz 到物理量温度与应变标定系数布里渊频移对温度和应变的响应在石英光纤中近似线性这也是 SBS 传感能定量测量的基础。常用经验系数如下物理量典型系数单位说明温度系数 C_T0.9~1.1MHz/°C随光纤掺杂浓度略有差异应变系数 C_\varepsilon0.048~0.058MHz/\mu\varepsilon随光纤涂覆层和纤芯成分变化解调公式为 \Delta\nu_B C_T \Delta T C_\varepsilon \Delta \varepsilon。做标定时最稳妥的做法是在被测光纤旁边放一段同型号的参考光纤段已知温度变化作为对照不要直接照抄厂商给出的常数不同批次光纤的实际系数可能有百分之几的偏差。4.4 温度应变交叉敏感的边界单根普通单模光纤上温度变化 1°C 和应变变化 20 \mu\varepsilon 造成的频移量相当因此仅凭一个 \nu_B 无法同时区分温度和应变。这是 SBS 传感的固有瓶颈不是拟合算法能解决的。常见做法是参考段补偿温度或者同时测量频移和线宽做双参量解耦但线宽同时受泵浦功率影响实用中要做功率校正。工程上多数场景会假定一段光纤内温度均匀只解应变如果你拿到的数据是单端 BOTDR 且没有温度参考通道就不要宣称同时测出了高精度的温度和应变。5. SBS 谱拟合质量的进阶验证残差、峰型畸变与双峰处理拟合做完下一步不是直接发布结果而是验证每个拟合峰是不是真的洛伦兹形。我的常规检查有三步第一把拟合残差做出来与谱线噪声本底的标准差比较残差峰峰值小于 3 倍噪声标准差说明峰形描述合理第二正反方向扫描的两组叠加数据分别拟合fB 差值应小于 1 MHz如果差出几十 MHz说明扫描步进或仪器滞后引入了不对称第三把质心法和拟合法的 fB 对比两者差超过线宽的 1/3 时叠加了双峰或非对称畸变的可能性很大。双峰叠加是保偏光纤中很常见的情况。两个偏振轴对应不同的有效折射率和声速布里渊频移出现两个峰间距从几 MHz 到几十 MHz。这时单峰洛伦兹拟合会把结果落在两个峰的中间某个加权位置物理意义不明确。处理方式是双洛伦兹叠加from scipy.signal import find_peaks def double_lorentz(f, a1, f1, d1, a2, f2, d2, c0): return (c0 a1 * (d1 / 2) ** 2 / ((f - f1) ** 2 (d1 / 2) ** 2) a2 * (d2 / 2) ** 2 / ((f - f2) ** 2 (d2 / 2) ** 2)) peaks, _ find_peaks(spec, prominence0.1 * (spec.max() - spec.min())) if len(peaks) 2: p0 [spec[peaks[0]] - spec.min(), freq[peaks[0]], 40, spec[peaks[1]] - spec.min(), freq[peaks[1]], 40, spec.min()] popt, _ curve_fit(double_lorentz, freq, spec, p0p0, maxfev50000)双峰拟合的初值不要靠猜先用 find_peaks 找到局部极大值分别作为两个 f 的初值幅度初值取各自峰值相对基底的高度。如果两个峰间距小到 5 MHz 以内拟合会高度相关这时需要固定一个峰的线宽先拟合出峰位再放开所有参数。另一个常见的峰型畸变来自高泵浦功率下的增益展宽和泵浦耗尽峰顶会被压平甚至凹陷。判断方法是固定本底噪声比较低功率和高功率两组数据中拟合线宽的变化若线宽随功率明显增大说明谱已进入功率展宽区——这种情况下测到的线宽不再等于声子寿命对应的本征线宽。本文还有配套的精品资源点击获取