如何无需初始猜测拟合质心与峰高随机分布的多高斯峰信号

问题描述 #

我是相关领域新手，发帖前已尽可能完成前置调研，若有无意疏漏还请包涵。
我正在从示波器采集电压-时间序列数据，时间bin宽度为0.8纳秒。我会运行多次「数据捕获」周期，单次捕获包含固定数量的样本，存在5_{15个高斯峰，峰的具体数量未知。这些高斯峰的半高全宽（FWHM）范围相对固定，在2}3纳秒之间，峰高不定，到达时间随机（即质心位置非周期性分布）。
我目前使用Python对数据进行高斯拟合，通过scipy.optimise库和astropy库已经取得了一定成果，下方为使用scipy.optimise的代码。目前我可以完成多高斯拟合，但代码中的关键步骤是需要手动提供峰数量的「初始猜测」，同时还要给出每个峰的质心位置、峰高、峰宽的估计值。有没有方法可以将代码通用化，无需手动提供初始猜测？如果放宽初始猜测的约束条件，拟合质量会明显下降。我已知所有峰都是宽度约束明确的高斯峰，希望能将代码通用化，可对任意单次捕获的数据拟合得到峰质心和峰高参数。

import ctypes
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

#Get data from file
with open('test3.txt') as f:
    w, h = [float(x) for x in next(f).split()]
    print(w, h)
    array = [[float(x) for x in line.split()] for line in f]

#Separate    
x, z = zip(*array)
#Change sign since fitting routine seems to
#prefer positive numbers
y=[ -p for p in z]

def func(x, *params):
    y = np.zeros_like(x)
    for i in range(0, len(params), 3):
        ctr = params[i]
        amp = params[i+1]
        wid = params[i+2]
        y = y + amp * np.exp( -((x - ctr)/wid)**2)
    return y

#Guess the peak positions, heights, and widths
guess = [16, 5, 2, 75, 5, 2, 105, 5, 2, 139, 5, 2, 225, 5, 2, 315, 5, 2, 330, 5, 2]

#Fit and print parameters to screen
popt, pcov = curve_fit(func, x, y, p0=guess)
print(popt)
fit = func(x, *popt)

#Plot stuff
plt.plot(x, y)
plt.plot(x, fit , 'r-')
plt.show()

解决方案 #

核心思路 #

利用已知的峰宽约束，结合自动峰检测生成初始参数，同时给拟合过程加边界限制，避免结果跑飞：

1. 用scipy.signal.find_peaks自动识别原始数据中的峰位置、峰高，过滤噪声产生的假峰
2. 根据高斯FWHM和sigma的换算关系（FWHM = 2.3548 * sigma），2_{3ns的FWHM范围对应sigma为0.85}1.27ns，直接将峰宽参数限制在该区间
3. 给质心、峰高参数也加合理的边界约束，保证拟合结果符合物理意义

完整修改代码 #

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

# 读取数据
with open('test3.txt') as f:
    w, h = [float(x) for x in next(f).split()]
    print(w, h)
    array = [[float(x) for x in line.split()] for line in f]

x, z = zip(*array)
x = np.array(x)
y = np.array([-p for p in z])  # 取反适配正峰检测

def multi_gauss(x, *params):
    y = np.zeros_like(x)
    for i in range(0, len(params), 3):
        ctr = params[i]
        amp = params[i+1]
        wid = params[i+2]
        y += amp * np.exp( -((x - ctr)/wid)**2)
    return y

# ------------------- 自动生成初始猜测 -------------------
# 峰检测参数可根据实际噪声水平调整
# width参数单位为采样点，0.8ns/bin对应2ns为2.5个bin，3ns为3.75个bin，所以设为(2,5)
# prominence参数用于过滤基线噪声，可根据实际数据调整
peaks, properties = find_peaks(y, width=(2, 5), prominence=0.5)
peak_centers = x[peaks]
peak_amps = properties['peak_heights']

# 控制峰数量在5~15之间，超过15个时保留高度最高的15个
if len(peak_centers) > 15:
    sort_idx = np.argsort(peak_amps)[::-1]
    peak_centers = peak_centers[sort_idx[:15]]
    peak_amps = peak_amps[sort_idx[:15]]
if len(peak_centers) <5:
    # 峰数不足时可降低prominence阈值重新检测
    pass

# 生成初始猜测：每个峰对应[质心, 峰高, 峰宽（默认取中间值1.06ns）]
guess = []
for ctr, amp in zip(peak_centers, peak_amps):
    guess.extend([ctr, amp, 1.06])

# 生成参数边界约束，格式为(所有参数的下限, 所有参数的上限)
lower_bounds = []
upper_bounds = []
for ctr, amp in zip(peak_centers, peak_amps):
    lower_bounds.extend([ctr - 1, amp * 0.5, 0.85])  # 质心±1ns，峰高±50%，峰宽0.85~1.27ns
    upper_bounds.extend([ctr + 1, amp * 1.5, 1.27])
bounds = (lower_bounds, upper_bounds)
# -----------------------------------------------------------

# 拟合时传入bounds参数
popt, pcov = curve_fit(multi_gauss, x, y, p0=guess, bounds=bounds)
print("拟合参数（按[质心, 峰高, 峰宽]循环排列）：\n", popt)
fit = multi_gauss(x, *popt)

# 绘图
plt.plot(x, y, label='原始数据')
plt.plot(x, fit , 'r-', label='拟合结果')
plt.scatter(peak_centers, peak_amps, c='g', marker='*', label='检测到的初始峰位置')
plt.legend()
plt.show()

调整说明 #

如果出现假峰漏检/多检的情况，可调整find_peaks的参数：

• 降低prominence值可以检测到更矮的峰
• 调整width范围适配实际峰宽
• 增加distance参数（单位为采样点）可以避免相邻过近的重叠峰被多次检测

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