欢迎点评！如果公式乱码，请用浏览器访问，用鼠标右键 -> 加载映像 显示公式。
需要安装的库：Jupyter，Jupyterlab，numpy，scipy，matplotlib，ipympl，mpl_interactions
用 pip 安装：
pip install Jupyter Jupyterlab numpy scipy matplotlib ipympl mpl_interactions
以下内容源自地下水动力学课程教学内容，程序可在 Jupyter notebook 中运行。
非线性拟合 — Theis 公式
Theis 公式是非线性的，所以需要将其转化为初始参数附近的线性形式，如用泰勒级数展开并忽略高阶项，我们要从初始参数开始寻找参数变化的最优步长。据此可使用最小二乘法拟合。
取初始参数 ，
，降深
的泰勒级数展开为
式中，。因为
，得
记
忽略高阶无穷小，(2) 式可写为
式 (4) 为 Theis 公式取参数（）时参数微小变化的线性近似式。
从初始参数出发，确定最优步长 ，可以得到一步优化的参数。
以 为未知系数,
为
的预测值，
根据初始参数（
）及观测时间计算，由此式可构造最小二乘问题。
程序设计需要考虑以下问题
- 最小二乘法得出的
有时为不合理的负值，这是需要收缩步长使参数为正值。可以用二分法缩小步长保证参数为正值；
- 最小二乘法只是在初值的基础上进行了一步优化。为了得到最优参数，需要用得出的参数作为新的参数初值迭代计算；
- 初值对计算效率有影响，简单方法是选取两组观测数据按 Jacob 公式计算参数初值。
直接上代码
%matplotlib widget
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import exp1
# 控制小数的显示精度
np.set_printoptions(precision=4)
def calc_u(p, r, t): # 计算变量 u
S, T = p
return r**2*S/4/T/t
# 定义井函数计算方法，可用 scipy.special.exp1 计算.
# 也可用多项式逼近的方法计算井函数
def theis_drawdown(p, t, Q, r): # 计算降深
S, T = p
u = calc_u(p, r, t)
s_theis = Q/4/np.pi/T*exp1(u)
return s_theis
def theis_Dfun(p, t, s, Q, r): # 计算 ds/dT,ds/S
"""Calculate and return ds/dS and ds/dT for the Theis equation.
The parameters S, T are packed into the tuple p.
"""
S, T = p
u = calc_u(p, r, t)
dsdS = -Q/4/np.pi/T/S*np.exp(-u)
dsdT = Q/4/np.pi/T**2*(-exp1(u) + np.exp(-u))
return np.array((-dsdS, -dsdT)).T # 返回负梯度
def theis_resid(p, t, s, Q, r): # 计算降深预测误差
S, T = p
return s - theis_drawdown(p, t, Q, r)
# 抽水量 Q, 观测孔位置 r
r = 125.0 # m, 15#孔位置
Q = 1440 # m^3/d
t=np.array([10., 20., 30., 40., 60., 80., 100., 120., 150., \
210., 270., 330., 400., 450., 645., 870., 990., 1185.]) / 1440 # 单位 day
s = np.array([0.16, 0.48, 0.54, 0.65, 0.75, 1., 1.12, 1.22, 1.36, 1.55, 1.7, \
1.83, 1.89, 1.98, 2.17, 2.38, 2.46, 2.54]) # m
# 取数组长度
n = len(t)
# 初始参数怎么取？最简单方法是取中间相邻的两点，用 Jacob 两点公式计算。
# 取数组长度
n = len(t)
i1 = int(n/2)
i2 = i1 + 1
t1 = t[i1]
t2 = t[i2]
s1 = s[i1]
s2 = s[i2]
kk = (s1 - s2)/np.log10(t1/t2)
T0 = 0.183*Q/kk
S0 = 2.25*T0*t1/r**2/np.float_power(10, s1/kk)
p = S0, T0
S, T = p
# 循环计数器
j = 0
while True:
c = theis_drawdown(p, t, Q, r)
a = theis_Dfun(p, t, s, Q, r)[:, 0]
b = theis_Dfun(p, t, s, Q, r)[:, 1]
X = np.vstack([a, b]).T # 形成系数矩阵
# 调用 np.linalg.solve
beta = np.linalg.solve(np.dot(X.T, X), np.dot(X.T, c - s))
DS = beta[0]
DT = beta[1]
while True:
if S + DS < 0: # 步子太大扯了蛋，参数为负不合理！
DS = DS/2.0 # 缩回半步
else:
break
while True:
if T + DT < 0: # 步长太大，参数为负数不合理！
DT = DT/2.0 # 缩回半步
else:
break
j = j + 1 # 循环计数器
if j > 1000: # 循环次数多，程序可能有错误
print(' 循环超过 1000 次，请先检查程序再运行！')
break
if (abs(DT) < 1.0E-6) or (abs(DS) < 1.0e-8): # 判断计算误差
break
# 新参数值
p = S + DS, T + DT
S, T = p
# 生成绘图数据
x0 =[i for i in np.arange(-3, 0.1, 0.1)]
x =np.power(10, x0)
y = theis_drawdown(p, x, Q, r)
plt.style.use('default')
fig, ax = plt.subplots(dpi=100)
ax.grid(True)
ax.set_xscale("log")
ax.set_yscale("log")
ax.set_xlim(0.001, 1)
ax.set_ylim(0.1, 10)
ax.set_aspect(1)
ax.plot(t, s, 'o', label='观测值')
ax.plot(x, y, label='拟合曲线')
plt.xlabel(r'$\log t$')
plt.ylabel(r'$\log s$')
plt.legend(
prop={'family': 'Kaiti', 'size': 10},
loc=4,title="图例",
title_fontproperties={'family': 'Simsun'})
# 指定图标题，显示中文
ax.set_title("Theis公式最小二乘法求参", fontproperties={'family': 'Simsun'})
plt.show()
# 显示最终结果
print('循环次数 = {}'.format(j))
print('T = {:.2f} m^2/d,'.format(T), 'S = {:.4e}'.format(S))
print('RSS = {:.4e}'.format(np.sqrt(np.sum((c - s)**2))))
image.png
