歡迎點評！如果公式亂碼袁辈，請用瀏覽器訪問马篮，用鼠標右鍵 -> 加載映像 顯示公式秃励。
需要安裝的庫：Jupyter傍念，Jupyterlab，numpy筏餐，scipy开泽，matplotlib，ipympl魁瞪，mpl_interactions
用 pip 安裝：
pip install Jupyter Jupyterlab numpy scipy matplotlib ipympl mpl_interactions
以下內(nèi)容源自地下水動力學課程教學內(nèi)容穆律，程序可在 Jupyter notebook 中運行。
非線性擬合 — Theis 公式
Theis 公式是非線性的导俘，所以需要將其轉(zhuǎn)化為初始參數(shù)附近的線性形式峦耘，如用泰勒級數(shù)展開并忽略高階項，我們要從初始參數(shù)開始尋找參數(shù)變化的最優(yōu)步長旅薄。據(jù)此可使用最小二乘法擬合辅髓。
取初始參數(shù) ，
少梁，降深
的泰勒級數(shù)展開為
式中洛口，。因為
凯沪，得
記
忽略高階無窮小第焰，(2) 式可寫為
式 (4) 為 Theis 公式取參數(shù)（）時參數(shù)微小變化的線性近似式。
從初始參數(shù)出發(fā)妨马，確定最優(yōu)步長 挺举，可以得到一步優(yōu)化的參數(shù)。
以 為未知系數(shù),
為
的預測值烘跺，
根據(jù)初始參數(shù)（
）及觀測時間計算湘纵，由此式可構造最小二乘問題。
程序設計需要考慮以下問題
- 最小二乘法得出的
有時為不合理的負值滤淳，這是需要收縮步長使參數(shù)為正值梧喷。可以用二分法縮小步長保證參數(shù)為正值；
- 最小二乘法只是在初值的基礎上進行了一步優(yōu)化铺敌。為了得到最優(yōu)參數(shù)绊困，需要用得出的參數(shù)作為新的參數(shù)初值迭代計算；
- 初值對計算效率有影響适刀，簡單方法是選取兩組觀測數(shù)據(jù)按 Jacob 公式計算參數(shù)初值。
直接上代碼
%matplotlib widget
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import exp1
# 控制小數(shù)的顯示精度
np.set_printoptions(precision=4)
def calc_u(p, r, t): # 計算變量 u
S, T = p
return r**2*S/4/T/t
# 定義井函數(shù)計算方法煤蹭，可用 scipy.special.exp1 計算.
# 也可用多項式逼近的方法計算井函數(shù)
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
# 取數(shù)組長度
n = len(t)
# 初始參數(shù)怎么缺屎怼？最簡單方法是取中間相鄰的兩點硝皂，用 Jacob 兩點公式計算常挚。
# 取數(shù)組長度
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
# 循環(huán)計數(shù)器
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 # 形成系數(shù)矩陣
# 調(diào)用 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: # 步子太大扯了蛋，參數(shù)為負不合理稽物！
DS = DS/2.0 # 縮回半步
else:
break
while True:
if T + DT < 0: # 步長太大奄毡，參數(shù)為負數(shù)不合理！
DT = DT/2.0 # 縮回半步
else:
break
j = j + 1 # 循環(huán)計數(shù)器
if j > 1000: # 循環(huán)次數(shù)多贝或，程序可能有錯誤
print(' 循環(huán)超過 1000 次吼过，請先檢查程序再運行！')
break
if (abs(DT) < 1.0E-6) or (abs(DS) < 1.0e-8): # 判斷計算誤差
break
# 新參數(shù)值
p = S + DS, T + DT
S, T = p
# 生成繪圖數(shù)據(jù)
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()
# 顯示最終結(jié)果
print('循環(huán)次數(shù) = {}'.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
