如何在Python中用有限能域Kramers-Kronig变换迭代拟合DAFS数据

我正在衍射异常精细结构（DAFS）实验中，尝试用Kramers-Kronig（KK）算法转换反常散射因子的实部（f'）与虚部（f''）贡献。

目前已用lmfit包实现实验数据的平滑曲线拟合，拟合函数如下：

def intensity(en, phi, beta, I0=1, slope=0, Ioff=0, fprime=-1, fsec=1):
    costerm = np.cos(phi) + beta*fprime
    sinterm = np.sin(phi) + beta*fsec
    return I0 * (costerm**2 + sinterm**2) + slope*en + Ioff

通过最小化参数，我用以下函数得到了f'的初始值：

def f1_guess(f2, I, I0, phi, beta, Ioff):
    f1_guess = (1/beta) * (np.sqrt(((I-Ioff)/I0) - (math.sin(phi)+(beta*f2))) - math.cos(phi)) 
    return f1_guess

现在我需要通过KK关系将f'的初始值转换为f''，再代入lmfit流程进行强度建模，进而计算精细结构。尝试过scipy.fftpack、scipy.integrate及sympy.integrate工具，但遇到两个核心问题：

1. 必须在有限能量范围内计算KK变换——需要移除已知原子线形f''a(E)，仅在实验能域内积分，而非标准的无限域积分。
2. 无法在脚本中准确区分积分式中的E与E'（或ω与ω'），不知道如何设置变量完成积分计算。

请问如何实现f'到f''的转换，并迭代lmfit模型，直到I0、Ioff、beta、phi这4个参数及f''/f'精细结构函数达到稳定？

一、有限域Kramers-Kronig变换的实现（f'→f''） #

首先明确有限域下的KK变换公式（针对f'转f''）：
对于实验能量范围[E_min, E_max]，f''(E)的表达式为：
$$f''(E) = f''a(E) + \frac{2E}{\pi} \text{PV} \int{E_{\text{min}}}^{E_{\text{max}}} \frac{f'(E') - f'_a(E')}{E'^2 - E^2} dE'$$
其中：

• $f'_a(E)$、$f''_a(E)$是原子的本征反常散射因子（可从XAFS数据库获取）
• PV表示柯西主值积分，用于处理E'=E时的奇点

代码实现步骤： #

1. 准备数据：

   • 把实验得到的f'初始值减去本征f'_a，得到精细结构部分$\Delta f'(E) = f'(E) - f'_a(E)$
   • 确保能量点E是均匀采样的（如果不是，先做插值处理）

2. 处理积分变量E与E'：

   • 用numpy的meshgrid生成E（目标能量点）和E'（积分变量）的二维网格，这样可以对每个E点计算对应的积分
   • 处理奇点：当E'≈E时，用极限值$\frac{d\Delta f'}{dE'}/(2E)$代替被积函数，避免除以0

3. 柯西主值积分计算：
   用scipy.integrate.quad结合主值选项实现：

import numpy as np
from scipy.integrate import quad

def kk_fprime_to_fsec(E, E_prime, delta_fprime, fsec_a):
    fsec = np.zeros_like(E)
    # 先对delta_fprime做插值，确保任意E'都能取值
    delta_fprime_interp = np.interp(E_prime, E, delta_fprime)
    
    for i, e in enumerate(E):
        # 定义被积函数
        def integrand(ep):
            return delta_fprime_interp[np.argmin(np.abs(E_prime - ep))]/(ep**2 - e**2)
        
        # 计算柯西主值积分，wvar指定奇点位置e
        pv_integral, _ = quad(integrand, E_prime.min(), E_prime.max(), weight='cauchy', wvar=e)
        fsec[i] = fsec_a[i] + (2*e/np.pi)*pv_integral
    return fsec

二、迭代lmfit拟合流程 #

要实现参数与f''/f'的稳定迭代，按以下步骤循环执行：

from lmfit import Model, Parameters

# 初始化模型
model = Model(intensity)
params = Parameters()
params.add('phi', value=0.5, min=0, max=np.pi)
params.add('beta', value=0.1, min=1e-5)
params.add('I0', value=1000, min=100)
params.add('Ioff', value=10, min=0)
params.add('slope', value=0, vary=False)  # 根据实验需求设置是否可变

# 加载实验数据（示例，替换为你的实际数据路径）
en_exp = np.loadtxt('energy.dat')
I_exp = np.loadtxt('intensity.dat')
fprime_a = np.loadtxt('fprime_a.dat')  # 本征f'
fsec_a = np.loadtxt('fsec_a.dat')      # 本征f''

# 初始f'（来自你的f1_guess函数）
fprime_initial = f1_guess(fsec_a, I_exp, params['I0'].value, params['phi'].value, params['beta'].value, params['Ioff'].value)
delta_fprime_initial = fprime_initial - fprime_a

# 迭代控制参数
max_iter = 15
tolerance = 1e-4
prev_chisqr = np.inf

for iter_num in range(max_iter):
    # 步骤1：用当前f'计算f''
    fsec_current = kk_fprime_to_fsec(en_exp, en_exp, delta_fprime_initial, fsec_a)
    # 固定f'和f''，拟合全局参数
    model.set_param_hint('fprime', value=fprime_initial, vary=False)
    model.set_param_hint('fsec', value=fsec_current, vary=False)
    
    # 步骤2：拟合I0、Ioff、beta、phi
    result = model.fit(I_exp, params, en=en_exp)
    
    # 步骤3：用拟合得到的全局参数重新计算f'
    new_fprime = f1_guess(fsec_current, I_exp, result.params['I0'].value, result.params['phi'].value, result.params['beta'].value, result.params['Ioff'].value)
    delta_fprime_new = new_fprime - fprime_a
    
    # 步骤4：检查收敛
    current_chisqr = result.chisqr
    if np.abs(current_chisqr - prev_chisqr) < tolerance:
        print(f"迭代收敛，次数：{iter_num+1}")
        break
    
    # 更新参数进入下一轮迭代
    prev_chisqr = current_chisqr
    fprime_initial = new_fprime
    delta_fprime_initial = delta_fprime_new
    params = result.params

# 输出最终拟合结果
print(result.fit_report())

收敛判断说明 #

• 用拟合的卡方值（chisqr）变化作为核心收敛依据，当两次迭代的卡方差小于设定阈值时停止
• 也可以额外监控f'或f''的全局均方误差，当误差小于阈值时辅助判断收敛

注意事项 #

• 确保能量点采样密度足够，避免积分误差过大；若采样不均，先做线性插值处理
• 本征散射因子$f'_a$和$f''_a$必须与实验能量范围完全匹配，可通过插值调整长度
• 如果积分结果噪声较大，可先对$\Delta f'$做平滑处理（比如用scipy.signal.savgol_filter）

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