基于Debye方程拟合实验数据：Matlab/Python问题求助

问题背景 #

目标是用Debye模型拟合20℃水的复介电常数实部实验数据，频段覆盖200MHz到6GHz，同时参考公开文献数据。标准Debye模型实部公式应为：
$$E’ = E(∞) + \frac{E(static) - E(∞)}{1 + (ωτ)^2}$$
（注：原问题中公式符号有误，正确形式下低频时$E’$趋近于$E(static)$，高频趋近于$E(∞)$）
其中：

• $E’$：复介电常数实部
• $E(∞)$：光频介电常数（预期值5.2）
• $E(static)$：静态介电常数（预期值≈80.2±0.2）
• $ω=2π×频率$（单位Hz）
• $τ$：弛豫时间（预期值9.36×10⁻¹²s）

当前遇到三个核心问题：

1. Matlab拟合不收敛 #

尝试了两个变体方程，但均未得到符合预期的参数：

• 错误形式：a+(b-a)/1+(2*3.142*x*c)^2（分母缺失括号，公式逻辑完全错误）
• 四参数变体：a+(b-a)/(1+(2*3.142*x*c)^2)^d

2. Python拟合参数偏差大 #

• 使用基础Debye方程拟合时，仅$E(static)$符合预期，其余参数偏差显著
• 使用四参数变体方程拟合得到结果：

  80.29080554191091
  1.0
  1.8969872910331487e-11
  1.0000322206749337
  

  其中$E(∞)$（对应第二个值）仅为1.0，弛豫时间推导值远小于预期的9.36e-12s

3. 数据文件读取报错 #

读取第二组数据时触发ValueError：

File ~\Documents\Akin\Antenna_Design\Dr Young\Dielectric_Extraction\CCT\555\250223\030323\FittingScript3.py:16 in 
experimentdata = np.array(rawdata[1:],dtype=float)
ValueError: could not convert string to float: ''

解决方案与优化建议 #

数据读取修复 #

报错原因是CSV文件中存在空字符串，需先清理数据再转换：

# 读取CSV后添加数据清理步骤
cleaned_data = []
for row in rawdata[1:]:
    # 过滤空值并转换为float
    valid_vals = [float(val.strip()) for val in row if val.strip() != '']
    cleaned_data.append(valid_vals)
experimentdata = np.array(cleaned_data, dtype=float)

同时检查CSV文件：确保没有 trailing commas、空白行或缺失数据列。

Python拟合优化 #

1. 修正模型公式与单位一致性

原模型定义存在两个关键问题：

• 公式符号错误（与标准Debye模型不符）
• 频率单位未转换（输入频率为MHz，需转换为Hz计算ω）

修正后的模型定义：

def debye_model(x, Estatic, Einf, tau, d=2):
    # x: 频率（MHz），需转换为Hz
    omega = 2 * np.pi * x * 1e6
    return Einf + (Estatic - Einf) / (1 + (omega * tau) ** d)

2. 使用物理意义的初始猜测

避免使用[1.0,1.0,1.0,1.0]这类无意义初始值，直接用预期参数作为起点：

# 基于文献预期值的初始猜测
init_guess = [80.2, 5.2, 9.36e-12, 2.0]

3. 添加参数边界约束

通过bounds限制参数在物理合理范围内，防止优化器发散：

bounds = (
    [79.0, 4.0, 1e-12, 1.5],  # 下界：Estatic≥79，Einf≥4，τ≥1e-12，d≥1.5
    [81.0, 6.0, 1e-11, 2.5]   # 上界：Estatic≤81，Einf≤6，τ≤1e-11，d≤2.5
)
# 执行拟合
fit_params, cov = curve_fit(debye_model, frequencydata, rpermdata, 
                           p0=init_guess, bounds=bounds, absolute_sigma=True)

4. 分离介电常数与电导率拟合

电导率的Debye模型与介电常数实部不同，不应复用同一模型。电导率的Debye形式为：
$$σ(ω) = σ_{dc} + \frac{(Estatic - Einf)ε_0 ω² τ}{1 + (ωτ)^2}$$
需单独定义模型进行拟合。

Matlab拟合修复 #

• 修正公式中的括号错误：将a+(b-a)/1+(2*3.142*x*c)^2改为a + (b-a)/(1 + (2*3.142*x*c)^2)
• 使用与Python相同的初始猜测和参数边界约束
• 确保频率单位转换正确（MHz→Hz）

修正后的完整Python代码示例 #

# -*- coding: utf-8 -*-
import numpy as np
import csv
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# 常数定义
e0 = 8.854e-12
pi = np.pi

# 读取并清理CSV数据
with open("water.csv", "r") as i:
    rawdata = list(csv.reader(i, delimiter=","))

cleaned_data = []
for row in rawdata[1:]:
    valid_vals = [float(val.strip()) for val in row if val.strip() != '']
    cleaned_data.append(valid_vals)
experimentdata = np.array(cleaned_data, dtype=float)

frequencydata = experimentdata[:, 0]
rpermdata = experimentdata[:, 1]
conductivitydata = experimentdata[:, 2]

# 计算电导率（单位S/m）
conductivity = 2 * pi * e0 * frequencydata * 1e6 * conductivitydata

# 绘制原始数据
plt.figure(1, dpi=120)
plt.title("20℃水介电常数Debye拟合")
plt.xlabel("频率 (MHz)")
plt.ylabel("相对介电常数/电导率")
plt.scatter(frequencydata, rpermdata, color='yellow', label="实测介电常数")
plt.scatter(frequencydata, conductivitydata, color="orange", label="实测电导率")

# 定义修正后的Debye模型（介电常数实部）
def debye_rperm(x, Estatic, Einf, tau, d):
    omega = 2 * pi * x * 1e6
    return Einf + (Estatic - Einf) / (1 + (omega * tau) ** d)

# 初始猜测与边界
init_guess = [80.2, 5.2, 9.36e-12, 2.0]
bounds = ([79.0, 4.0, 1e-12, 1.5], [81.0, 6.0, 1e-11, 2.5])

# 拟合介电常数
fit_rperm, cov_rperm = curve_fit(debye_rperm, frequencydata, rpermdata, 
                                p0=init_guess, bounds=bounds, absolute_sigma=True)
Estatic_fit, Einf_fit, tau_fit, d_fit = fit_rperm

# 打印拟合结果
print("拟合参数：")
print(f"E(static) = {Estatic_fit:.2f}")
print(f"E(∞) = {Einf_fit:.2f}")
print(f"弛豫时间τ = {tau_fit:.2e} s")
print(f"指数d = {d_fit:.4f}")

# 生成拟合曲线
x_line = np.arange(min(frequencydata), max(frequencydata), 1)
rperm_fit_line = debye_rperm(x_line, *fit_rperm)

# 绘制拟合曲线
plt.plot(x_line, rperm_fit_line, '--', color='black', label="Debye拟合介电常数")
plt.legend()
plt.show()

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