Python中使用curve_fit进行双高斯拟合的问题

问题分析 #

原拟合结果不佳的核心原因包括：

• find_peaks参数设置不合理，峰值定位存在偏差
• 初始参数猜测精度不足，峰位置和宽度的估算受数据离散性影响较大
• curve_fit未设置参数约束，拟合过程中参数易偏离合理范围
• 模型未考虑基线噪声，导致整体拟合偏移

优化步骤 #

• 精准定位峰值：调整find_peaks的distance和height参数，确保两个峰被正确识别分离。
• 优化初始参数：用二次插值细化峰位置，结合稳健的半高宽计算宽度参数，同时加入基线初始值（取数据基线区域均值）。
• 添加参数边界：通过curve_fit的bounds参数限制参数范围（如高度为正、宽度为正、峰位置在数据区间内），防止拟合发散。
• 引入基线项：在双高斯模型中加入常数基线，适配数据的背景噪声。

修改后的代码 #

import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import find_peaks
from scipy.optimize import curve_fit
from scipy.interpolate import interp1d

# 原始数据
y_data = np.array([
    1.500e-04, 1.500e-04, 1.500e-04, 1.500e-04, 1.700e-04, 1.600e-04, 1.800e-04,
    1.600e-04, 1.700e-04, 2.300e-04, 2.500e-04, 3.200e-04, 3.200e-04, 3.800e-04,
    4.000e-04, 5.000e-04, 5.600e-04, 5.600e-04, 6.500e-04, 7.500e-04, 9.100e-04,
    1.180e-03, 1.550e-03, 2.110e-03, 2.880e-03, 3.850e-03, 5.200e-03, 6.780e-03,
    8.950e-03, 1.123e-02, 1.403e-02, 1.723e-02, 2.031e-02, 2.330e-02, 2.495e-02,
    2.433e-02, 2.171e-02, 1.725e-02, 1.231e-02, 8.080e-03, 4.980e-03, 2.960e-03,
    1.970e-03, 1.830e-03, 2.220e-03, 2.880e-03, 3.700e-03, 4.650e-03, 5.820e-03,
    7.150e-03, 8.450e-03, 9.510e-03, 1.006e-02, 9.660e-03, 8.560e-03, 6.910e-03,
    5.100e-03, 3.380e-03, 2.170e-03, 1.230e-03, 7.000e-04, 3.600e-04, 2.200e-04,
    1.700e-04, 1.500e-04, 1.600e-04, 1.100e-04, 1.100e-04, 1.200e-04, 1.000e-04,
    1.200e-04, 1.100e-04, 1.300e-04, 1.500e-04, 1.200e-04, 1.600e-04, 1.200e-04,
    1.200e-04, 1.200e-04, 1.200e-04, 1.100e-04, 1.500e-04, 1.500e-04, 1.300e-04,
    1.100e-04, 8.000e-05, 1.200e-04, 1.200e-04, 1.100e-04, 1.100e-04, 1.500e-04])

x_data = np.array([
    6555.101, 6555.201, 6555.301, 6555.401, 6555.501, 6555.601, 6555.701, 6555.801,
    6555.901, 6556.001, 6556.101, 6556.201, 6556.301, 6556.401, 6556.501, 6556.601,
    6556.701, 6556.801, 6556.901, 6557.001, 6557.101, 6557.201, 6557.301, 6557.401,
    6557.501, 6557.601, 6557.701, 6557.801, 6557.901, 6558.001, 6558.101, 6558.201,
    6558.301, 6558.401, 6558.501, 6558.601, 6558.701, 6558.801, 6558.901, 6559.001,
    6559.101, 6559.201, 6559.301, 6559.401, 6559.501, 6559.601, 6559.701, 6559.801,
    6559.901, 6560.001, 6560.101, 6560.201, 6560.301, 6560.401, 6560.501, 6560.601,
    6560.701, 6560.801, 6560.901, 6561.001, 6561.101, 6561.201, 6561.301, 6561.401,
    6561.501, 6561.601, 6561.701, 6561.801, 6561.901, 6562.001, 6562.101, 6562.201,
    6562.301, 6562.401, 6562.501, 6562.601, 6562.701, 6562.801, 6562.901, 6563.001,
    6563.101, 6563.201, 6563.301, 6563.401, 6563.501, 6563.601, 6563.701, 6563.801,
    6563.901, 6564.001, 6564.001])

# 优化峰值检测：设置distance确保两个峰被正确分离
peaks, _ = find_peaks(y_data, height=0.001, distance=20)

# 定义带基线的双高斯函数
def double_gaussian(x, a1, b1, c1, a2, b2, c2, baseline):
    g1 = a1 * np.exp(-((x - b1) ** 2) / (2 * c1 ** 2))
    g2 = a2 * np.exp(-((x - b2) ** 2) / (2 * c2 ** 2))
    return g1 + g2 + baseline

# 更稳健的FWHM计算函数
def find_fwhm(x, y, peak_idx):
    half_max = y[peak_idx] / 2
    # 左半部分插值找半高位置
    left_points = y[:peak_idx+1][::-1]
    left_xs = x[:peak_idx+1][::-1]
    left_idx = np.argmax(left_points <= half_max)
    if left_idx == 0:
        left_x = left_xs[0]
    else:
        x0, x1 = left_xs[left_idx-1], left_xs[left_idx]
        y0, y1 = left_points[left_idx-1], left_points[left_idx]
        left_x = x0 + (half_max - y0) * (x1 - x0)/(y1 - y0)
    
    # 右半部分插值找半高位置
    right_points = y[peak_idx:]
    right_xs = x[peak_idx:]
    right_idx = np.argmax(right_points <= half_max)
    if right_idx == 0:
        right_x = right_xs[0]
    else:
        x0, x1 = right_xs[right_idx-1], right_xs[right_idx]
        y0, y1 = right_points[right_idx-1], right_points[right_idx]
        right_x = x0 + (half_max - y0) * (x1 - x0)/(y1 - y0)
    return right_x - left_x

# 用二次插值细化峰位置
interp_func = interp1d(x_data, y_data, kind='quadratic')
x_fine = np.linspace(x_data[peaks[0]-2], x_data[peaks[0]+2], 1000)
y_fine = interp_func(x_fine)
peak1_x = x_fine[np.argmax(y_fine)]

x_fine2 = np.linspace(x_data[peaks[1]-2], x_data[peaks[1]+2], 1000)
y_fine2 = interp_func(x_fine2)
peak2_x = x_fine2[np.argmax(y_fine2)]

# 峰高取插值后的最大值
peak1_a = np.max(y_fine)
peak2_a = np.max(y_fine2)

# 计算FWHM并转换为高斯宽度
fwhm1 = find_fwhm(x_data, y_data, peaks[0]) / 2.355
fwhm2 = find_fwhm(x_data, y_data, peaks[1]) / 2.355

# 基线初始值：取数据前10个点的均值
baseline_guess = np.mean(y_data[:10])

# 初始参数猜测
initial_guess = [peak1_a, peak1_x, fwhm1, peak2_a, peak2_x, fwhm2, baseline_guess]

# 设置参数边界，避免拟合发散
bounds = (
    [0, min(x_data), 0.01, 0, min(x_data), 0.01, 0],  # 参数下界
    [max(y_data)*1.2, max(x_data), 2, max(y_data)*1.2, max(x_data), 2, max(y_data)*0.1]  # 参数上界
)

# 执行带边界约束的拟合
params, covariance = curve_fit(double_gaussian, x_data, y_data, p0=initial_guess, bounds=bounds)

# 生成拟合曲线数据
x_fit = np.linspace(min(x_data), max(x_data), 10000)
fitted_data = double_gaussian(x_fit, *params)

# 分离单个高斯分量和基线
g1 = params[0] * np.exp(-((x_fit - params[1]) ** 2) / (2 * params[2] ** 2))
g2 = params[3] * np.exp(-((x_fit - params[4]) ** 2) / (2 * params[5] ** 2))
baseline = params[6] * np.ones_like(x_fit)

# 绘图展示
plt.scatter(x_data, y_data, label='原始数据', s=10)
plt.plot(x_fit, g1 + baseline, 'g--', label='高斯峰1')
plt.plot(x_fit, g2 + baseline, 'm--', label='高斯峰2')
plt.plot(x_fit, fitted_data, 'r-', label='双高斯拟合曲线', linewidth=2)
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

优化效果说明 #

修改后的代码能够：

• 精准定位两个峰的位置，避免拟合错位
• 拟合曲线与原始数据的匹配度显著提升
• 高斯分量的高度、宽度和位置与实际数据特征对齐

内容的提问来源于stack exchange，提问作者bougab