2. 洛伦兹拟合的布里渊频移估计

任务说明

本案例针对分布式光纤传感中布里渊频移（BFS）的估计问题。输入为一段实测布里渊增益谱（BGS）数据（频率-增益对），目标是通过洛伦兹曲线拟合估计谱峰对应的频移值，并定量比较不同拟合初值对收敛结果的影响。你需要完成：加载数据、实现洛伦兹拟合、输出BFS估计值、对比不同初值下的估计差异。

洛伦兹拟合模型

布里渊增益谱近似满足洛伦兹线型：

g(f) = g0 * (ΔfB/2)² / [(f - fB)² + (ΔfB/2)²] + noise

其中待估参数为：g0（峰值增益）、fB（布里渊频移）、ΔfB（线宽，即半高全宽）。采用非线性最小二乘（Levenberg-Marquardt）进行拟合。

代码实现

import numpy as np
from scipy.optimize import curve_fit
import json

def lorentzian(f, g0, fB, dfB):
    """洛伦兹线型函数"""
    return g0 * (dfB/2)**2 / ((f - fB)**2 + (dfB/2)**2)

def load_bgs(filepath):
    """加载实测BGS数据，返回频率数组(MHz)和增益数组(线性)"""
    data = np.loadtxt(filepath, delimiter=',')
    return data[:,0], data[:,1]

def estimate_bfs(freq, gain, p0):
    """
    洛伦兹拟合估计布里渊频移
    参数:
        freq: 频率数组 (MHz)
        gain: 增益数组 (线性)
        p0: 初值 [g0, fB, dfB]
    返回:
        popt: 拟合参数 [g0, fB, dfB]
        pcov: 协方差矩阵
    """
    # 剔除异常点（增益非正或NaN）
    mask = np.isfinite(gain) & (gain > 0)
    freq, gain = freq[mask], gain[mask]
    
    # 约束参数范围：g0>0, fB在数据范围内, dfB>0
    bounds = ([0, freq.min(), 0], [np.inf, freq.max(), np.inf])
    
    popt, pcov = curve_fit(lorentzian, freq, gain, p0=p0, bounds=bounds, maxfev=5000)
    return popt, pcov

def main():
    # 加载实测数据（示例数据文件格式：每行"频率,增益"）
    freq, gain = load_bgs('bgs_data.csv')
    
    # 从数据中提取粗略初值
    idx_max = np.argmax(gain)
    f_peak = freq[idx_max]
    g_peak = gain[idx_max]
    
    # 计算半高宽作为线宽初值
    half_max = g_peak / 2
    above_half = gain > half_max
    if np.sum(above_half) >= 2:
        f_above = freq[above_half]
        df_rough = f_above.max() - f_above.min()
    else:
        df_rough = 50.0  # 默认线宽50MHz
    
    # 定义三组不同初值
    p0_list = [
        [g_peak, f_peak, df_rough],           # 初值1：基于数据峰值
        [g_peak*0.8, f_peak - 30, df_rough*1.5],  # 初值2：频移偏小30MHz
        [g_peak*1.2, f_peak + 30, df_rough*0.5],  # 初值3：频移偏大30MHz
    ]
    
    results = []
    for i, p0 in enumerate(p0_list):
        try:
            popt, pcov = estimate_bfs(freq, gain, p0)
            fB_est = popt[1]
            fB_err = np.sqrt(pcov[1,1]) if np.all(np.isfinite(pcov)) else np.nan
            results.append({
                'init': i+1,
                'p0_fB': p0[1],
                'fB_estimate': fB_est,
                'fB_std': fB_err,
                'converged': True
            })
        except Exception as e:
            results.append({
                'init': i+1,
                'p0_fB': p0[1],
                'fB_estimate': None,
                'fB_std': None,
                'converged': False,
                'error': str(e)
            })
    
    # 输出结果
    print("=== 洛伦兹拟合BFS估计结果 ===")
    for r in results:
        if r['converged']:
            print(f"初值{r['init']}: 初值fB={r['p0_fB']:.2f} MHz -> 估计fB={r['fB_estimate']:.4f} MHz (标准差={r['fB_std']:.4f} MHz)")
        else:
            print(f"初值{r['init']}: 初值fB={r['p0_fB']:.2f} MHz -> 拟合失败: {r.get('error','未知错误')}")
    
    # 保存结果
    with open('bfs_estimation_results.json', 'w') as f:
        json.dump(results, f, indent=2)
    print("\n结果已保存至 bfs_estimation_results.json")

if __name__ == '__main__':
    main()

数据准备与运行方式

准备实测BGS数据文件 bgs_data.csv，每行包含两个数值：频率（MHz）和增益（线性幅度），以逗号分隔。若暂无实测数据，可用以下模拟数据生成脚本创建测试文件：

import numpy as np
# 生成模拟BGS数据（真实fB=10600 MHz, 线宽=80 MHz）
f = np.linspace(10400, 10800, 200)
true_fB = 10600.0
true_df = 80.0
true_g0 = 1.0
gain = true_g0 * (true_df/2)**2 / ((f - true_fB)**2 + (true_df/2)**2)
gain += np.random.normal(0, 0.01, size=f.shape)  # 加入1%噪声
np.savetxt('bgs_data.csv', np.column_stack([f, gain]), delimiter=',')
print("模拟数据已生成: bgs_data.csv")

运行环境需安装 numpy 和 scipy。执行主程序：

python bfs_lorentzian_fit.py

验证方法与判定标准

运行程序后，观察终端输出。判定标准如下：

• 收敛性：三组初值均应成功收敛（converged=True），若某组初值导致拟合失败，说明初值选取对算法稳定性有影响，需记录失败原因。
• 估计一致性：三组初值估计的fB应接近（差异小于1 MHz），表明拟合结果对初值不敏感；若差异显著（大于5 MHz），说明存在局部最优解问题。
• 与真实值对比：若使用模拟数据（真实fB=10600 MHz），估计值应在10600±2 MHz范围内。
• 标准差合理性：fB_std应远小于估计值（通常小于1 MHz），反映拟合精度。

通过比较不同初值下的估计结果，可以评估算法对初值的依赖程度：若所有初值收敛到同一值，说明拟合鲁棒；若部分初值不收敛或收敛到不同值，则需要改进初值选取策略（如使用峰值搜索或谱矩方法提供更优初值）。

结果解读

程序输出三组初值对应的BFS估计值及标准差。若三组结果差异小于1 MHz，说明拟合算法对初值不敏感，BFS估计可靠；若某组初值导致不收敛或收敛到明显不同的峰值，则表明初值选取需谨慎。建议以数据峰值附近作为初值（初值1），并对比偏移初值的结果来评估拟合稳定性。最终BFS估计值取三组收敛结果的平均值，标准差反映拟合不确定度。