Python中stijghoogteverlaging函数曲线拟合异常问题求助

看起来你在使用curve_fit进行抽水试验数据拟合时，遇到了函数逻辑和效率的双重问题，我来帮你拆解并修正：

核心问题诊断 #

你之前的两个写法都踩了curve_fit的使用误区：

1. 分支内直接return：循环第一次迭代就触发return，函数只返回第一个元素的计算值，而curve_fit需要函数返回和输入数组长度一致的结果，这直接导致拟合曲线完全错误。
2. 循环后单个return：循环里反复覆盖s变量，最后只返回最后一个元素的值，同样不符合curve_fit的要求；而且Python原生循环处理数组效率极低，加上curve_fit会多次调用函数，直接导致计算耗时爆炸。

修正后的完整代码 #

我用numpy向量化操作重构了你的核心函数，既解决逻辑问题又大幅提升效率：

from scipy.special import exp1
import numpy as np
from scipy.optimize import curve_fit
import pandas as pd
import matplotlib.pyplot as plt

# 数据读取与预处理（保留你的原始逻辑）
data = pd.read_csv("pompproef_data.csv", sep=';')
pb1 = data.iloc[1:,1].values - 1.87
pb2 = data.iloc[1:,2].values - 1.86
t1 = data.iloc[1:,0].values / (60*24)  # 转换为小时单位
volume = 10/1000  # m3
duur = [128,136, 150, 137, 143, 141]  # 秒
totaal = np.sum(duur)
debiet = (((len(duur) * volume)/totaal)) * (60*60*24)  # 转换为m3/d
print(debiet)
print(t1)
print(pb1)

# 全局参数定义
tuit = 15/(24*60)  # 转换为小时
D = 2.0
Q = debiet
r = 4.0  # 观测井到抽水井的距离（米）

# 重构后的水位降深计算函数
def stijghoogteverlaging(t, k, S):
    # 提取公共计算因子，避免重复运算
    common_factor = Q / (4 * np.pi * k * D)
    # 计算基础项（适用于所有t < tuit的情况）
    term1 = exp1(S * r**2 / (4 * k * D * t))
    s = common_factor * term1
    
    # 对t >= tuit的元素，计算差值项并更新结果
    mask = t >= tuit
    if np.any(mask):
        term2 = exp1(S * r**2 / (4 * k * D * (t[mask] - tuit)))
        s[mask] = common_factor * (term1[mask] - term2)
    
    return s

# 执行参数拟合
poptpb1, pcovpb1 = curve_fit(stijghoogteverlaging, t1, pb1, p0=[100, 1e-5], maxfev=10000000)
print('optimale waarde van k voor peilbuis1:', poptpb1[0])
print('optimale waarde van S voor peilbuis1:', poptpb1[1])

# 生成拟合曲线并绘图
tijd = data.iloc[1:,0].values
# 基于原始时间范围生成拟合点（保持小时单位）
t = np.linspace(t1.min(), t1.max(), 1000)
s1 = stijghoogteverlaging(t, poptpb1[0], poptpb1[1])

plt.plot(tijd, pb1, 'r.', label='Gemeten bij 4 meter')
# 将拟合时间转换为分钟，与原始数据单位对齐
plt.plot(t * 24 * 60, s1, 'b', label='fitted bij 4 m')
plt.xlabel('Tijd (minuten)')
plt.ylabel('Stijghoogteverlaging (m)')
plt.legend()
plt.show()

关键修正点说明 #

1. 向量化替代循环：用numpy数组操作替代Python原生循环，效率提升几个数量级，同时自然满足curve_fit对输出长度的要求。
2. 正确使用函数参数：原函数错误地引用全局变量t1，现在改为使用curve_fit传入的t参数，这是拟合逻辑正确的核心。
3. 布尔掩码处理分支：用mask = t >= tuit筛选需要计算差值项的元素，避免提前return的问题，保证所有元素都被正确计算。
4. 单位统一优化：绘图时将拟合时间从小时转换为分钟，与原始数据的时间单位对齐，避免曲线错位。
5. 初始参数调整：把S的初始值从1e-25调整为1e-5，更符合水文地质参数的常规范围，帮助拟合更快收敛。

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