Theis井流广义群井抽水试验参数反演最优化方法

Optimization Method for Parameter Inversion of Generalized Group Well Pumping Test in Theis Well Flow

摘要: 群井抽水试验、阶梯流量抽水试验、间断性阶梯流量试验、抽水 + 恢复试验，本文统称为广义群井抽水试验。传统的手工配线法无法解决广义群井抽水试验的参数反演；简化的s-logt分析法，不能完整利用全时程试验数据。现有的研究，侧重于利用抽水阶段或恢复阶段的试验数据进行参数反演，不考虑前期抽水试验对恢复水位的影响，割裂了恢复过程与抽水过程的有机联系。本文尝试利用现代计算工具Mathcad及1stOpt，对Theis井流广义群井抽水试验进行全时程最优化反演，并考虑抽水井的井损效应。利用4个经典实例对广义群井抽水试验最优化参数反演进行了实例验证，获得了十分满意的效果。最优化求解，水文地质参数T、S可以不进行条件约束，但井损系数c需要进行条件约束：c ≥ 0。Mathcad的最小误差解函数Minerr()，是搜寻目标方程E = 0的解，即使方程无解，也会给出与方程解误差最小的一组近似解。最小误差解并不一定严格满足约束条件，井损系数c有可能解出负值，不符合井损与流量呈正相关的预期。为此，可将井损系数c改用绝对值|c|进行最优化解算。

Abstract: The group well pumping test, stepped flow pumping test, intermittent stepped flow test, and pumping + recovery test are collectively referred to as the generalized group well pumping test in this article. The traditional manual wiring method cannot solve the parameter inversion of generalized group well pumping tests; The simplified s-logt analysis method cannot fully utilize the full time experimental data. The existing research focuses on parameter inversion using experimental data from the pumping or recovery stages, without considering the impact of previous pumping tests on the recovery water level, thus cutting off the organic connection between the recovery process and the pumping process. This article attempts to use modern computing tools Mathcad and 1stOpt to perform full time optimization inversion on the generalized group well pumping test of Theis well flow, considering the well loss effect of the pumping well. Four classic examples were used to verify the optimization parameter inversion of generalized group well pumping test, and satisfactory results were obtained. Optimization solution, hydrogeological parameters T and S can be unconstrained, but the well loss coefficient c needs to be constrained: c ≥ 0. The minimum error solution function Minerr() in Mathcad is used to search for the solution of the objective equation E = 0. Even if the equation has no solution, it will provide a set of approximate solutions with the minimum error from the equation solution. The minimum error solution may not strictly satisfy the constraint conditions, and the well loss coefficient c may result in negative values, which does not meet the expected positive correlation between well loss and flow rate. For this purpose, the well loss coefficient c can be optimized by using the absolute value |c|.

1. 引言

Theis井流是地下水非稳定井流模型的经典模型。

水文地质参数通常是通过抽水试验对井流模型进行反演解算获得。抽水试验，按抽水井数量，通常分为单井及群井抽水试验(多个抽水井)；按抽水降次，可分为单降次及多降次抽水试验。多降次抽水试验，各降次之间不进行恢复试验时，通常也称之为阶梯流量抽水试验。

阶梯流量抽水试验、间断性阶梯流量抽水试验、抽水 + 恢复试验，为广义阶梯流量抽水试验，均可转换为等价的同孔群井抽水试验。

本文将群井抽水试验、广义阶梯流量抽水试验，统称为广义群井抽水试验。

地下水水文地质参数的反演解算，传统的经典方法是利用单孔抽水阶段或恢复阶段的试验资料，采用logs-logt手工配线法以及简化的s-logt分析法(图解法)进行求解。对恢复试验，是将“准稳定”后的恢复试验仿照独立的抽水试验进行反演，不考虑前期抽水试验对恢复水位的影响，割裂了恢复过程与抽水过程的有机联系。

上世纪80年代，开始利用计算机及袖珍计算机进行水文地质参数的反演[1]-[5]，目前的研究热点是最优化反演[6]-[10]。有作者利用2点降深比值法及其变形算法解算导水系数T及释水系数S [4] [5]，张勇等(2019)深入研究了水文地质参数反演的智能算法[9]，洪声亮等(2025)使用Mathcad及1stOpt软件分别对Theis井流模型的抽水阶段数据及恢复阶段数据进行了参数反演解算[10]。

目前，对广义群井抽水试验全时程数值反演研究较少，需进一步探索研究。

本文尝试使用现代计算工具Mathcad及1stOpt，对Theis井流广义群井抽水试验进行全时程最优化反演，并考虑抽水井的井损效应。利用4个经典实例对广义群井抽水试验的群井抽水试验、阶梯流量抽水试验、间断性阶梯流量抽水试验、抽水 + 恢复试验等4种试验模型进行了参数反演实例验证，获得了十分满意的效果。

2. Theis井流单井抽水试验模型

2.1. Theis井流单井模型

经典的Theis非稳定井流模型，含水层为无界承压含水层，假定含水层为水平无界、均质各向同性、等厚、无补给，抽水量为含水层存储水量的弹性释放。Theis井流单井模型详见图1。

Figure 1. Theis well flow single well pumping test model

图1. Theis井流单井抽水试验模型

Theis给出了该模型的水位降深解[11]-[13]：

$s=\frac{Q}{4\pi T}\cdot w\left(u\right)$ (1)

$w\left(u\right)={\displaystyle {\int }_u^{\infty }\frac{e^{-x}}{x}}dx=Ei\left(u\right)$ (2)

$u=\frac{r^2}{4at}=\frac{S\cdot r^2}{4T\cdot t}$ (3)

$T=k\cdot M$ (4)

$a=\frac{T}{S}$ (5)

式中：

$w\left(u\right)$ ：井函数；

$Ei\left(u\right)$ ：指数积分函数；

$s$ ：降深/回升(m)；

$Q$ ：井出水量(m³/d)；

$M$ ：含水层厚度(m)；

$r$ ：井距(m)；

$t$ ：抽水延续时间/停抽恢复时间(min)；

$k$ ：渗透系数(m/d)；

$a$ ：导压系数(m²/d)；

$T$ ：导水系数(m²/d)；

$S$ ：弹性释水系数。

2.2. Theis井函数$w\left(u\right)$ 的数值计算

井函数$w\left(u\right)$ 的计算，本文引用文献[10]的研究成果，采用变量代换后的积分式：

$w\left(u\right)={\displaystyle {\int }_{\ln u}^4\exp \left(e^y\right)}dy$ (2-1)

3. Theis井流广义群井抽水试验模型

群井抽水试验(详见图2)，是由2个(含2个)以上抽水井进行的抽水试验，各井的抽水流量及开抽时间可以不同。

Figure 2. Pumping test diagram of group wells, (a) Well layout plan; (b) Q-t curve

图2. 群井抽水试验图，(a) 井孔布置图；(b) Q-t曲线

设各抽水井到观察井的井距$r$ 、流量$Q$ 及开抽时间$tp$ 为：

$\{\begin{array}{l}r=\left[r_1,r_2,\cdot \cdot \cdot ,r_j,\cdot \cdot \cdot ,r_n\right]\\ Q=\left[Q_1,Q_2,\cdot \cdot \cdot ,Q_j,\cdot \cdot \cdot ,Q_n\right]\\ tp=\left[t{p}_1,t{p}_2,\cdot \cdot \cdot ,t{p}_j,\cdot \cdot \cdot ,t{p}_n\right]\end{array}$ (6)

推论：广义阶梯流量抽水，等价于同孔群井抽水。

设抽水试验的阶梯流量$Qt$ 及其开抽时间$tp$ 为：

$\{\begin{array}{l}Qt=\left[Q{t}_1,Q{t}_2,\cdot \cdot \cdot ,Q{t}_j,\cdot \cdot \cdot ,Q{t}_n\right]\\ tp=\left[t{p}_1,t{p}_2,\cdot \cdot \cdot ,t{p}_j,\cdot \cdot \cdot ,t{p}_n\right]\end{array}$ (7)

则，由叠加原理，可将阶梯流量抽水试验等价转换为同孔群井抽水试验(详见图3)：

$\{\begin{array}{l}r=\left[r_0,r_0,\cdot \cdot \cdot ,r_0,\cdot \cdot \cdot ,r_0\right]\\ Qt=\left[Q{t}_1,Q{t}_2,\cdot \cdot \cdot ,Q{t}_j,\cdot \cdot \cdot ,Q{t}_n\right]\\ Q=\left[Q_1,Q_2,\cdot \cdot \cdot ,Q_j,\cdot \cdot \cdot ,Q_n\right]\\ tp=\left[t{p}_1,t{p}_2,\cdot \cdot \cdot ,t{p}_j,\cdot \cdot \cdot ,t{p}_n\right]\\ Q_1=Q{t}_1\\ Q_j=Q{t}_j-Q{t}_{j-1}\end{array}$ (8)

Figure 3. Step flow pumping test diagram, (a) Qt-t curve diagram of stepped flow pumping; (b) Equivalent group well pumping Q-t curve graph

图3. 阶梯流量抽水试验图，(a) 阶梯流量抽水Qt-t曲线图，(b) 等价群井抽水Q-t曲线图

3.1. 广义群井降深通解

广义群井抽水试验，利用叠加原理，得到水位降深通解s [11]：

$s={\displaystyle \sum _{j=1}^n\frac{Q_j}{4\pi T}\cdot w\left(u_j\right)}$ (9)

当考虑井损时，叠加井损项$\Delta s$ [11]：

$\Delta s={\displaystyle \sum _{j=1}^{n}c\cdot Q{t}_j^2}c\ge 0$ (10)

$u_j=\frac{r_j^2}{4\cdot a\cdot t_j}$ (11)

$t_j=t-t{p}_j$ (12)

其中：

$Q_j$ ：为群井第j井流量(或等价群井第j井流量)；

$Q{t}_j$ ：为第j阶阶梯流量；

$t_j$ ：为$Q_j$ 的抽水延续时间；

$t{p}_j$ ：为第j阶阶梯流量开抽时间；

$t$ ：为抽水试验总延续时间；

$c$ ：为井损系数；

3.2. 广义群井降深简化解

当$u_j<0.05$ 时，由式(9)得到水位降深简化解[11]：

$s={\displaystyle \sum _{j=1}^n\left[\frac{Q_j}{4\pi T}\cdot \ln \left(2.25\cdot a\right)+\frac{Q_j}{4\pi T}\cdot \ln \left(\frac{t_j}{r_j^2}\right)\right]}$ (9-1)

4. Theis井流广义群井抽水试验参数最优化反演方法

本文仅讨论广义群井抽水试验的通解反演，不涉及简化解反演。

广义群井抽水试验，水位降深$s$ 与水文地质参数$\left(T,S\right)$ 为非线性关系，难以使用人工配线法求取水文地质参数。

由式(9)，水位降深为水文地质参数$(T,S)$ 的函数:

$s\left(T,S,r,t\right)={\displaystyle \sum _{j=1}^n\frac{Q_j}{4\pi T}\cdot w\left(T,S,r_j,t_j\right)}$ (13)

当考虑抽水井井损时，叠加井损项：

$s\left(T,S,r,t\right)={\displaystyle \sum _{j=1}^n\frac{Q_j}{4\pi T}\cdot w\left(T,S,r_j,t_j\right)}+{\displaystyle \sum _{j=1}^{n}c\cdot Q{t}_j^2}c\ge 0$ (14)

基于最小二乘法的参数求解，就是使实测降深$sg$ 与理论降深$s$ 的二乘目标函数E最小：

$\min E\left(T,S\right)={\displaystyle \sum _{i=1}^n{\left(s{g}_i-s{}_i\right)}^2}$ (15)

目标函数E值，也即是残差平方和SSE。

本文探索使用2款非常实用的数学软件Mathcad及最优化专业软件1stOpt进行反演解算。

Mathcad优化解可使用最小误差函数Minerr()及最小优化函数Minimize()，为局部最优解。

Minerr()函数，专用于求解最小二乘法目标函数的最优化问题，其优点在于即使初值偏离最优解较远或无严格解时，仍可寻找一组局部近似解–最小误差解。Minimize()函数，是更为一般的通用最优化方法，寻找约束条件下的最小优化解。两者算法不同，其解也可能有一定差异。实际求解时，可对参数赋不同初值，将解算结果进行对比，找到符合水文地质条件的一组局部最优解。

1stOpt，是国内七维高科有限公司开发的最优化软件，解算速度快，代码非常简洁，其最大特色是不依赖于参数初值即可大概率地找到全局最优解。

4.1. Mathcad反演代码

使目标函数最小的最优化求解，Mathcad是通过Given定义求解块实现的。

1) 输入试验数据：

$r：={\left(r_1,r_2,\cdot \cdot \cdot ,r_i,r_n\right)}^{\text{T}}$

$Qt：={\left(Q{t}_1,Q{t}_2,\cdot \cdot \cdot ,Q{t}_i,Q{t}_n\right)}^{\text{T}}$

$Q：={\left(Q_1,Q_2,\cdot \cdot \cdot ,Q_i,Q_n\right)}^{\text{T}}$

$tp：={\left(t{p}_1,t{p}_2,\cdot \cdot \cdot ,t{p}_i,t{p}_n\right)}^{\text{T}}$

$tg：={\left(t{g}_1,t{g}_2,\cdot \cdot \cdot ,t{g}_i,t{g}_m\right)}^{\text{T}}$

$sg：={\left(s{g}_1,s{g}_2,\cdot \cdot \cdot ,s{g}_i,s{g}_m\right)}^{\text{T}}$

2) 定义函数

$u\left(T,S,r,t\right)：=\frac{1440\cdot S\cdot r^2}{4\cdot T\cdot t}$

$w\left(T,S,r,t\right)：={\displaystyle {\int }_{\ln \left(u\left(T,S,r,t\right)\right)}^4\exp \left(-e^y\right)}dy$

$sc\left(T,S,c,t\right)：={\displaystyle \sum _{j=1}^n\frac{Q_j}{4\cdot \pi \cdot T}\cdot w\left(T,S,r_j,t\right)}+{\displaystyle \sum _{j=1}^{n}c\cdot Q{t}_j^2}c\ge 0$

不考虑井损时，删除井损参数c及井损叠加项，下同。

3) 参数赋初值

$\left(\begin{array}{l}T\\ S\\ c\end{array}\right)：=\left(\begin{array}{l}{T}_0\\ S_0\\ c_0\end{array}\right)$

可根据情况给出不同的初值$\left(\begin{array}{l}{T}_0\\ S_0\\ c{}_0\end{array}\right)$

4) 定义二乘目标函数

$Ec\left(T,S,c\right)：={\displaystyle \sum _{i=0}^{m-1}{\left(sc\left(T,S,c,t{g}_i\right)-s{g}_i\right)}^2}$

5) 局部优化解

5a. Minerr解

$\begin{array}{l}\text{Given }//定义求解块\\ c\ge 0//约束条件\\ Ec\left(T,S,c\right)：=0//目标方程\end{array}$

$\left(\begin{array}{l}T\\ S\\ c\end{array}\right)：=\mathrm{Minerr}\left(T,S,c\right)//最小误差解$

5b. Minimize解

$\begin{array}{l}\text{Given }//定义求解块\\ c\ge 0//约束条件\\ \left(\begin{array}{l}T\\ S\\ c\end{array}\right)：=\mathrm{Minimize}\left(E,T,S,c\right)//最小优化解\end{array}$

6) 计算残差平方和

$SSE：=E\left(T,S,c\right)//残差平方和越小越好$

4.2. 1stOpt反演代码

1stOpt最优化反演代码也很简洁，由7部分构成：

1) 定常量赋值

Constant n= ; //广义群井数n

2) 变常量赋值(快捷方式数组赋值)

Constant r (n) = [r₁, r₂ , r_i, …, r_n];

Constant Qt (n) = [Qt₁, Qt₂ , Qt_j, …, Qt_n];

Constant Q (n) = [Q₁, Q₂ , Q_j, …, Q_n];

Constant tp (n) = [tp₁, tp₂, tp_j, …, tp_n];

3) 定义变量

Variable t, s;

4) 参数赋初值与约束(参数赋初值非必需)

Parameters TT = [0,];

Parameters SS = [0,];

Parameters c = [0,]; //不考虑井损时，删除本行

5) 字符串常量赋字(井函数w)

ConstStr w = Expint(u); //1stOpt直接用指数积分函数计算井函数w(u)

6) 定义拟合函数(Theis水位降深函数s，pascal编程方式)

Start Program [Pascal];

Procedure Main Model;

var i, j: integer;

var tj, u: double;

Begin

for i: = 0 to DataLength - 1 do begin //i数据长度循环

s [i]: = 0.0;

for j: = 1 to n do begin //j群井数循环

tj: = t[i]-tp[j] ; //tj群井流量Qj作用时间

if tj > 0 then begin

u: = 1440*r[j]^2*SS/(4*TT*tj);

s[i]: = s[i]+Q[j]/(4*pi*TT)*w ; //层流降深

Ds: = c*(Qt[j])^2 ; //井损项

s[i]: = s[i]+ Ds ; //总降深

end; //end if

end; //end j

end;// end i

End;// end begin

EndProgram;

7) 抽水试验t-s数据(m组数据)

Data.

t1, t2 , ti … tm ;

s1, s2 , ti … sm ;

注：因1stOpt不区分大小写字母，代码中TT为导水系数系数T，SS为弹性释水系数S。

5. Theis井流广义群井抽水试验参数最优化反演实例

5.1. 群井抽水试验

【实例1】选自：经典习题集，文献[13]，p.119，习题16。

某承压含水层布置了井1、井2、井3及观测孔4进行群井抽水试验。各井抽水流量分别为Q₁ = 6000 m³/d，Q₂ = 4000 m³/d，Q₃ = 2000 m³/d。观测孔至各抽水井的井距分别为r₁ = 150 m，r₂ = 600 m及r₃ = 200 m。抽水试验开始时，2号井先开泵，5 min后3号井先开泵，10 min后井1号井才开泵。井孔平面布置及Q-t曲线见图4，抽水试验s-t数据见表1。试利用上述资料求解含水层的导水系数T及弹性释水系数S。(注：原文献为习题，无题解答案)

Figure 4. Layout of group well pumping test boreholes, (a) Well hole plan; (b) Q-t curve

图4. Theis井流群井抽水试验，(a) 井孔平面图；(b) Q-t曲线

Table 1. Theis well group pumping test data

表1. Theis井流群井抽水试验数据表

	群井			观测孔
井号	井1	井2	井3	G1
井距r/m	150	600	200
延续时间	Q1	Q2	Q3	s	备注
t /min	m^3/d	m^3/d	m^3/d	m	tp/min
0		4000			第2井开抽
2		4000		0.05
5		4000	2000	0.20	第3井开抽
10	6000	4000	2000	0.32	第1井开抽
15	6000	4000	2000	0.38
20	6000	4000	2000	0.44
30	6000	4000	2000	0.50
40	6000	4000	2000	0.56
60	6000	4000	2000	0.65
90	6000	4000	2000	0.70
120	6000	4000	2000	0.73
180	6000	4000	2000	0.79
240	6000	4000	2000	0.84
300	6000	4000	2000	0.90
360	6000	4000	2000	0.93
420	6000	4000	2000	0.96
480	6000	4000	2000	0.98
600	6000	4000	2000	1.01
720	6000	4000	2000	1.04
840	6000	4000	2000	1.05
1020	6000	4000	2000	1.09	停抽

Table 2. Theis well group pumping test parameter calculation table

表2. Theis井流群井抽水试验参数解算表

解算 工具	解算 方法	赋参数初值		最优解	参数解算值
		T m^2/d	S	Minerr /Minimize	T m^2/d	S	c d^2/m^5	SSE m^2
Mathcad	局部最优	100	0.01		2086.058	2.693E−03		1.093750
					4753.586	3.604E−04		0.174538
		1000	0.001		6363.952	1.159E−04		0.077033
					6972.551	7.557E−05		0.071702
		10000	0.0001		6973.595	7.527E−05		0.071699
					6973.591	7.527E−05		0.071699
				局部最优解	6973.595	7.527E−05		0.071699
1stOpt	全局最优			全局最优解	6973.593	7.527E−05		0.071699

Figure 5. The s-t fitting curve of the pumping test for the Theis well flow group

图5. Theis井流群井抽水试验s-t拟合曲线图

本例，抽水量较大，但观测井距离抽水井较远，不考虑观测井水位的井损效应。根据群井抽水试验资料，最优化解算结果详见表2，s-t拟合曲线图见图5。Mathcad软件局部最优解与1stOpt的全局最优解一致。

s-t拟合曲线图，在前10 min拟合不好，水位下降值大于理论值，可能受到了其它因素的影响。

5.2. 阶梯流量抽水试验

【实例2】选自：经典教材，文献[11]，p.147，例题，表4-7-1。

在某地某承压井进行了4个阶梯流量的抽水试验。抽水井井径0.109 m，阶梯流量Qt分别为：3000,6000,8000,9000 m³/d；各阶梯流量的开抽时间tp为：0,120,240,360 min。阶梯流量抽水试验资料详见表3及图6。原文献考虑了井损问题，用s-logt分析法进行了解答。

Table 2. Theis well group pumping test parameter calculation table

表2. Theis井流群井抽水试验参数解算表

延续时间	恢复时间	阶梯流量	延续时间	主井	备注
t/min	t′/min	Qt/m^3/d	t/min	s/m	tp/min
0		3000	0	0.000	第1阶开抽
1		3000	1	3.630
2		3000	2	3.790
5		3000	5	4.010
10		3000	10	4.180
20		3000	20	4.340
50		3000	50	4.560
120		6000	120	4.770	第2阶开抽
121		6000	121	8.940
122		6000	122	9.110
125		6000	125	9.330
130		6000	130	9.510
140		6000	140	9.690
170		6000	170	9.950
240		8000	240	10.250	第3阶开抽
241		8000	241	13.330
242		8000	242	13.450
245		8000	245	13.600
250		8000	250	13.730
260		8000	260	13.870
290		8000	290	14.070
360		8000	360	14.350	第4阶开抽
361		9000	361	15.980
362		9000	362	16.040
365		9000	365	16.130
370		9000	370	16.190
380		9000	380	16.290
410		9000	410	16.420
480		9000	480	16.630	停抽

Figure 6. Ladder flow pumping test diagram and equivalent group well flow pumping test diagram, (a) Step flow pumping; (b) Equivalent group well pumping

图6. 阶梯流量抽水及等价群井流量抽水试验图，(a) 阶梯流量抽水；(b) 等价群井抽水

根据阶梯抽水试验资料，主井开抽流量较大，为3000 m/d，宜优先考虑井损效应。解算结果详见表4，s-t拟合曲线详图7。

Table 4. Parameter calculation table for Theis well flow step flow pumping test (considering well loss)

表4. Theis井流阶梯流量抽水试验参数解算表(考虑井损)

解算 工具	解算 方法	赋参数初值			最优解	参数解算值
		T m^2/d	S	c min^2/m^5	Minerr /Minimize	T m^2/d	S	c d^2/m^5	SSE m^2
Mathcad	局部 最优解	10	0.1	0		443.384	3.879E−01	0.000E+00	6.565778
						469.568	2.472E−01	0.000E+00	5.140931
		100	0.01	0		638.761	8.291E−03	0.000E+00	1.315364
						634.797	8.783E−03	0.000E+00	1.344217
		1000	0.001	0		971.268	3.662E−05	7.206E−09	0.346288
						721.797	1.532E−03	0.000E+00	1.233325
					局部最优解	971.268	3.662E−05	7.206E−09	0.346288
1stOpt	全局最优解				全局最优解	919.558	8.730E−05	6.513E−09	0.333317
原教材	图解法				s-logt图解法	1016.000	1.016E−04	3.00E−08	75.452973

Figure 7. s-t fitting curve of Theis well flow step flow pumping test (considering well loss), (a) 1stOpt Global Optimal Solution; (b) Original manual s-logt graphical solution

图7. Theis井流阶梯流量抽水试验s-t拟合曲线图(考虑井损)，(a) 1stOpt全局最优化解；(b) 原文手工s-logt图解法解

从残差平方和SSE来看，Mathcad的局部解与1stOpt的全局解大致相当，1stOpt解好于Mathcad解。手工s-logt分析法的SSE远远大于1stOpt及Mathcad，参数解算精度低。

从s-t曲线拟合来看，1stOpt及athcad解整体拟合较好(图7(a))；手工s-logt分析解整体拟合差，100 min钟后，s-t曲线越往后偏离越来越大(图7(b))。

5.3. 间断性阶梯流量抽水试验

【实例3】选自：论文，文献[5]。

抽水试验井位于某城郊北部平原区，岩性较单一，上部12 m~14 m亚粘土。静止水位埋深7 m左右，含水层厚120 m~130 m，具承压性。抽水试验井井径0.15 m，在具抽水井中心16 m一观测孔。抽水进行了3个阶梯流量的抽水试验Qt₁ = 1500 m³/d，Qt₂ = 0 m³/d，Qt₃ = 2000 m³/d，试验过程中采用间断性阶梯流量抽水。600 min时Qt₂停抽，1800 min时刻Qt₃开抽。抽水试验资料见表5及图8。

Table 5. Data table of intermittent step flow pumping test for Theis well flow

表5. Theis井流间断性阶梯流量抽水试验数据表

日期	时间	延续时间	阶梯抽水 延续时间	阶梯流量	延续时间		主井 降深	阶梯抽水 起始时间	备注
y-m-d	h: min	t/min	t′/min	Qt/m^3/d	t/min		s/m	tp/min
		0		1500	0		0.000	0	第1阶起始
		20	20	1500	20		4.172
		25	25	1500	25		4.658
		30	30	1500	30		5.061
		40	40	1500		40	5.710
		50	50	1500		50	6.220
		60	60	1500		60	6.641
		70	70	1500		70	6.998
		80	80	1500		80	7.309
		90	90	1500		90	7.585
		100	100	1500		100	7.832
		120	120	1500		120	8.262
		150	150	1500		150	8.790
		180	180	1500		180	9.222
		210	210	1500		210	9.589
		240	240	1500		240	9.904
		300	300	1500		300	10.423
		360	360	1500		360	10.860
		480	480	1500		480	11.551
		600		0		600	12.807	600	第2阶起始
		720		0		720	4.263
		900		0		900	2.639
		1200		0		1200	1.665
		1500		0		1500	1.227
		1800	0	2000		1800	0.974	1800	第3阶起始
		2100	300	2000		2100	14.705
		2400	600	2000		2400	16.808
		3000	1200	2000		3000	18.872		试验结束

Figure 8. Theis intermittent stepped flow pumping test, (a) Intermittent staircase flow; (b) Equivalent group well

图8. Theis间断性阶梯流量抽水试验，(a) 间断性阶梯流量；(b) 等价群井

间断性阶梯流量抽水，开抽流量大，为1500 m³/d，宜优先考虑井损。

将间断性阶梯抽水作为整体进行全时程解算，最优解算结果见表6，s-t拟合曲线见图9。

解算结果表明，井损可忽略不计，Mathcad的局部最优解与1stOpt的全局最优解完全相同。1stOpt及Mathcad的SSE不足原文解的一半，最优化解算值好于原文解。

Table 6. Parameter calculation table for Theis well flow step flow pumping test (considering well loss)

表6. Theis井流间断性阶梯流量抽水试验参数解算表(考虑井损)

解算 工具	解算 方法	赋参数初值		最优解	参数解算值
		T m^2/d	S	Minerr /Minimize	T m^2/d	S	c d^2/m^5	SSE m^2
Mathcad	局部最优	10	0.1		1.174	1.705E−01	2.864E−06	266.752400
					49.060	9.668E−04	0.000E+00	0.563540
		100	0.01		52.575	1.019E−03	2.032E−07	2.158688
					49.060	9.668E−04	0.000E+00	0.563540
		1000	0.001		862.518	8.102E−08	2.459E−06	181.450700
					49.063	9.666E−04	0.000E+00	0.563540
				局部最优解	49.060	9.668E−04	0.000E+00	0.563540
1stOpt	全局最优			全局最优解	49.060	9.668E−04	0.000E+00	0.563540
原文解	2点降深比值法				51.770	8.050E−04		1.189440

Figure 9. The s-t fitting curve of intermittent stepped flow pumping test for Theis

图9. Theis间断性阶梯流量抽水试验s-t拟合曲线图

5.4. 抽水 + 恢复试验

【实例4】选自：经典教材，文献[12]，p.136，例题，抽水，表4-2；p.146，例题，恢复，表4-5。

1976年在江苏丰县进行了承压含水层的非稳定流抽水试验。抽水试验场地位于江苏丰县王沟公社曹楼大队冲洪积平原，承压含水层为粉细砂及亚砂土(粉土)，抽水井井深120 m，直径0.20 m，自1976年11月6日9时40分开始抽水，抽水后进行了恢复观测。抽水持续97 h，流量为22.60 m³/h = 542.4 m³/d。观测孔1距抽水井117.85 m。抽水试验数据详见表7。

Table 7. Theis well flow pumping + recovery test data table

表7. Theis井流抽水 + 恢复试验数据表

日期	时间	延续时间	恢复时间	阶梯 流量	观1	观1	延续时间	恢复时间	观1	观1
y-m-d	h: min	t/min	t′/min	Qt/m^3/d	降深水位 s/m	回升水位 s′/m	t/min	t′/min	降深水位 s/m	回升水位 s′/m
1976/11/6	9:40	0		542.4	0.000	开抽	5825	5	1.729	0.001
		8		542.4	0.002		5826	6	1.728	0.002
		10		542.4	0.005		5828	8	1.727	0.003
		12		542.4	0.006		5830	10	1.725	0.005
		15		542.4	0.007		5832	12	1.721	0.009
		20		542.4	0.008		5835	15	1.719	0.011
		25		542.4	0.011		5840	20	1.715	0.015
		30		542.4	0.020		5845	25	1.710	0.020
		40		542.4	0.029		5850	30	1.706	0.024
		50		542.4	0.038		5855	35	1.702	0.028
		60		542.4	0.050		5860	40	1.698	0.032
		80		542.4	0.093		5865	45	1.692	0.038
		100		542.4	0.130		5870	50	1.685	0.045
		120		542.4	0.165		5880	60	1.670	0.060
		150		542.4	0.220		5900	80	1.640	0.090
		180		542.4	0.270		5920	100	1.605	0.125
		210		542.4	0.330		5940	120	1.570	0.160
		240		542.4	0.370		5970	150	1.515	0.215
		300		542.4	0.465		6000	180	1.465	0.265
		360		542.4	0.530		6030	210	1.400	0.330
		480		542.4	0.655		6060	240	1.360	0.370
		600		542.4	0.755		6120	300	1.273	0.457
		720		542.4	0.880		6180	360	1.200	0.530
		900		542.4	1.000		6300	480	1.080	0.650
		1200		542.4	1.150		6420	600	0.980	0.750
		1500		542.4	1.220		6540	720	0.895	0.835
		1800		542.4	1.320		6720	900	0.794	0.936
		2100		542.4	1.390		7020	1200	0.650	1.080
		2400		542.4	1.450		7320	1500	0.560	1.170
		3000		542.4	1.510		7620	1800	0.452	1.278
		3600		542.4	1.670		7920	2100	0.390	1.340
		4200		542.4	1.710		8520	2700	0.260	1.470
		4800		542.4	1.720		10040	4220	0.090	1.640
		5820		542.4	1.730	停抽， 恢复

观测孔1距抽水井117.85 m，不需考虑井损。将抽水试验 + 恢复试验作为一个整体，转化为等价的群井抽水，对其进行全时程最优化解算。解算结果见表8及图10。

解算结果表明，Mathcad的局部最优解与1stOpt的全局最优解结果完全相同。Mathcad与1stOpt的SSE，不足手工配线法的1/4，解算结果好于原教材手工配线法。从s-t曲线拟合来看，最优化解的s-t曲线整体拟合好，仅在恢复末期(停抽1800 min，30 h后)拟合较差(图10)，水位回升快于理论值，此时可能受到外界或抽排水渗漏引起水位回升的影响。

Table 8. Calculation table of Thetis well flow pumping + recovery test parameters

表8. Theis井流抽水 + 恢复试验参数解算表

解算 工具	解算 方法	赋参数初值			最优解	参数解算结果
		T m^2/d	S	c d^2/m^5	Minerr /Minimize	T m^2/d	S	c d^2/m^5	SSE m^2
Mathcad	局部最优	10	0.1			98.163	1.211E−03		0.373405
						98.163	1.211E−03		0.373405
		100	0.01			98.163	1.211E−03		0.373405
						98.163	1.211E−03		0.373405
		1000	0.001			140.975	2.064E−04		6.041790
						98.163	1.211E−03		0.373405
					局部最优解	98.163	1.211E−03		0.373405
1stOpt	全局最优				全局最优解	98.163	1.211E−03		0.373405
原教材	手工配线法				logs-logt配线法	80.160	1.540E−03		1.682171

注：手工配线法解算值为纯抽水试验段lgs-lgt手工配线解算值。

Figure 10. The s-t fitting curve of Theis well flow pumping and recovery test, (a) 1stOpt global solution; (b) logs-logt manual wiring solution

图10. Theis井流抽水 + 恢复试验s-t拟合曲线图，(a) 1stOpt全局解；(b) logs-logt手工配线解

6. 最优化解算的几点说明

1. 关于最优化算法

1stOpt非线性最优化算法，有通用全局最优化算法(UGO)、稳健全局优化法、简面体爬山法、差分进化法、最大继承法、经典局部优化算法、遗传算法、拟退火法、粒子群法等。其独特的通用全局优化算法UGO，极大地提高了非线性寻优能力，无需对参数赋初值即可大概率地找到全局最优解。一般只需采用默认算法UGO。

Mathcad的局部优化算法，一般也只需采用默认设置。最小误差解函数Minerr()可供选择的优化算法有列文伯格法、共轭梯度法及拟牛顿法。最小优化解函数Minimize()可供选择的优化算法只有共轭梯度法及拟牛顿法。

2. 关于求解参数的约束

层流降深式中，参数的定义域为$u>0$ ，$T>0$ ，$S>0$ 。因此，在最优化求解时，可以不对其进行约束。

考虑井损时，井损与流量呈正相关，因此，最优化求解时需要对井损系数进行约束：$c\ge 0$ 。

3. Mathad的最优解

Mathcad的最小误差解函数 Minerr()，是搜寻目标方程$E=0$ 的解，即使方程无解，也会给出与方程解误差最小的一组近似解。这个解值并不一定严格满足约束条件，有可能解出井损系数c为负值，不符合井损与流量呈正相关预期。为次，可将井损项的井损系数$c$ 改用绝对值$\left|c\right|$ 进行优化求解。

最小优化解Minimize()函数，严格按约束条件给出满足约束条件的数值解。

Minerr解与Minimize解并不一定相等，有时甚至相差较大。

4. 关于模型的适用性

Theis井流模型是一假定的理想模型，实际情况与Theis井流模型不一定完全吻合。在抽水试验中，也还可能存在各种自然的和非自然的干扰因素，导致流量、时间及水位等的观测数据受到一定程度的影响。因此，实际的s-t曲线与理论s-t曲线不一定能完全拟合。在实际解算时，应注意这些因素的影响。

7. 结语

广义群井抽水试验的参数反演，常规的人工配线法难以解算；简化的s-logt分析法，只能利用部分试验资料，求解精度不高，s-t曲线整体拟合较差。

利用观测孔资料反演水文地质参数，一般不需要考虑井损问题。利用抽水井资料解算水文地质参数，当抽水量较大时，由于三维流、紊流、及过滤器水头损失等影响，一般需优先考虑井损问题。

本文使用现代数学工具Mathcad及最优化工具1stOpt，利用4个经典实例对广义群井抽水试验的群井抽水试验、阶梯流量抽水、间断性阶梯流量抽水试验、以及抽水 + 恢复试验等4种试验模型进行了参数反演验证，获得了十分满意的效果。

Mathcad局部最优解与1stOpt全局最优解，两者结果基本相同，1stOpt的全局最优解更为精准。

最优化求解时，水文地质参数T，S可以无约束求解，但井损系数$c$ 应进行约束求解。

Mathcad的最小误差解函数Minerr()，是搜寻目标方程$E=0$ 的解，即使方程无解，也会给出与方程解误差最小的一组近似解。这个解值并不一定严格满足约束条件，有可能解出负值，不符合井损与流量呈正相关预期。为此，可将井损系数$c$ 改用绝对值$\left|c\right|$ 进行优化求解。

Theis井流模型是一假定的理想模型，实际情况与Theis井流模型不一定完全吻合。在抽水试验中，有可能存在各种自然的和非自然的干扰因素，导致观测数据受到一定程度的影响，从而使实际的s-t曲线与理论s-t曲线不一定能完全拟合。在实际解算时，应注意这些因素的影响。

致 谢

七维高科有限公司张伟博士对1stOpt软件的使用给予了大力支持，对本文给予了热情指导，在此深以致谢！

NOTES

^*第一作者。

^#通讯作者。