如何用Scipy优化带参数约束的光谱图多高斯峰拟合？

问题背景 #

现有代码用于拟合XRF能谱数据，需满足以下约束：

• 固定4个高斯峰的mean值（[24.210, 24.002, 27.276, 27.238] keV）
• 固定3个高斯峰的振幅与第一个峰的比例（[1, 0.5538, 0.1673, 0.1673*0.5185]）
• 仅保留两个拟合自由度：第一个高斯的振幅、所有高斯的公共sigma
  原代码拟合效果偏差，需优化实现上述约束并提升拟合精度。

优化后代码 #

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from google.colab import drive

# 挂载Google Drive
drive.mount('/content/drive')

# 读取数据文件
def read_numeric_data(file_path):
    with open(file_path, 'r', encoding='latin1') as f:
        data = []
        for line in f:
            stripped_line = line.strip()
            if all(char.isdigit() or char in "-+.e" for char in stripped_line):
                data.extend([float(x) for x in stripped_line.split()])
    return np.array(data)

# 加载数据
file_path = '/content/drive/MyDrive/ProjetoXRF_Manuel/April2024/20240402_In.mca'
data = read_numeric_data(file_path)

# 能量刻度参数
a = -0.0188026396003431
b = 0.039549044037714

# 固定参数配置
MEANS = np.array([24.210, 24.002, 27.276, 27.238])
AMP_RATIOS = np.array([1, 0.5538, 0.1673, 0.1673 * 0.5185])

# 多高斯拟合函数（仅保留两个自由度：主振幅A、公共sigma）
def multi_gaussian(x, A, sigma):
    eps = 1e-10
    # 计算每个高斯分量并求和
    gaussians = A * AMP_RATIOS * np.exp(-(x[:, np.newaxis] - MEANS)**2 / (2 * sigma**2 + eps))
    return np.sum(gaussians, axis=1)

# 数据预处理
working_data = data[:885]
max_counts = working_data.max()
norm_data_y = working_data / max_counts
data_x = a + b * np.arange(len(norm_data_y))

# 初始参数设置：主振幅设为1（归一化后），sigma初始值0.1
p0 = [1.0, 0.1]

# 设置参数边界：振幅>0，sigma>0
bounds = ([0, 1e-3], [2, 1.0])

# 加权拟合：XRF数据误差为sqrt(counts)，归一化后误差为sqrt(counts)/max_counts
data_err = np.sqrt(working_data) / max_counts

# 执行拟合
params, params_cov = curve_fit(
    multi_gaussian,
    data_x,
    norm_data_y,
    p0=p0,
    bounds=bounds,
    sigma=data_err,
    absolute_sigma=True
)

# 提取拟合结果
fit_A, fit_sigma = params
print(f"拟合主振幅: {fit_A:.4f}")
print(f"公共sigma: {fit_sigma:.4f}")

# 绘制结果
fig, ax = plt.subplots(figsize=(10, 6))
# 原始数据
ax.scatter(data_x, norm_data_y, color='black', marker='o', s=10, label='原始数据')
# 误差棒
ax.errorbar(data_x, norm_data_y, yerr=data_err, fmt='none', ecolor='black', capsize=2)
# 拟合曲线
x_fine = np.linspace(22, 28, 2000)
y_fit = multi_gaussian(x_fine, fit_A, fit_sigma)
ax.plot(x_fine, y_fit, color='red', linewidth=1.5, label='拟合曲线')

# 图表配置
ax.set_xlabel('能量 (keV)')
ax.set_ylabel('归一化计数')
ax.set_title('In元素XRF能谱多高斯拟合')
ax.set_xlim(22, 28)
ax.legend()
plt.show()

关键优化点 #

• 重构拟合函数：移除原函数冗余参数，仅保留A（主振幅）和sigma（公共标准差）两个自由参数，严格遵循约束条件
• 加权拟合：利用XRF数据的统计误差（sqrt(counts)）作为拟合权重，提升拟合精度
• 参数边界约束：设置振幅和sigma的合理范围，避免拟合过程中出现非物理的负参数
• 简化初始参数：仅传入两个初始值，减少拟合算法的搜索空间，提升收敛稳定性

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