Theis井函数计算方法及井模型参数优化计算研究

1
背景
Thies 井 模 型是水文 地 质 和 地 下 水 研究 领域 中 的
重要模型方法之一， 主要用于承压含水层抽水试验资 料建模， 并进而确定水文地质参数如渗透系数 K 和导 是用何种方法计算出模型参数㊂ 水系数 T 等， 其在国内外都有广泛的应用㊂ 该模型的 应用有两个要点， 一是如何 计算模 型 中 的 井函 数 ， 二 解， 该方法简单㊁易于理解和实现， 但误差大的缺点也 关于第一点， 使用最广的方式是级数展开近似求
式 中 ：S 为 实 测 降 深 (m)；Q 为 抽 水 流 量 (m3/d)；T 为 承 压含水层导水系数(m2/d)；W (u)为 Theis 井函数；us 为承
ʏ
exp(-x) dx,u= r2 • us x 4 • T • t
Q • W (u) 4 • π • T
(1)
压含水层贮水系数；t 为抽水时间(d)；r 为观测 孔到抽 水井的距离(m)㊂ T 和 us， 难点之一是公式中的井函数 W (u)无解析解， 只 学上也叫指数积分函数 (ExpInt)， 在水文地质中则常 根据抽 水 试 验 资料 和 模 型公式 (1)确定 模 型参数
N n
由于 Theis 井函数为指数积分式, 数学上 无法 直
W (u)=
i -ln(u)+ða• i u i=0 4
5
式 中 ：γ =0.5772156649015329， 是 欧拉 常 数 ；N 理论 上 趋近于无穷㊂ 展开 项 N㊂ 国内 绝 大 多数 研究 者 都 是采用的这种方 似 数 值 计算 方法会 因 为 级 数 项 的取 舍 而导 致实际应 5~11 时的 u-W (u)图， 可以很明显地看出当 u 大于 2 时 大,其后果是模型参数最终求解结果出现较大误差㊂
-i
2.2 Barry 方法 法,概要如下：
[5]
3,i
-i
Barry 于 2000 年提出了一种估计井函数 W (u)的方 G 1-G exp(-u)ln 1+ u ( h + b • u)2 W (u)= -u G +(1-G ) • exp 1-G
式中：p 1㊁p 2㊁p 3㊁q 1㊁q 2 和 q 3 为常数项系数, 具体值可参 考文献[7]㊂ 2.5 拟合方程替代 巩彦文[8]等提出 了 利用 拟 合 方程 来近 似计算 井 函 W (u)=18.109 • u0.256-31.310 • u0.16+13.482 (6)
2.5 1.5 0.5 -0.5 -1.5 2 1 0 N=5 N=7
(2)
1 • i=0 4 u • exp(u)
ðb•u ðc•u
i=0 i
i
u>1
(4)
实际应用中, 可以根据不同的精度需要选取级数
序号 a b c
法,级数项 N 一般都取 3~5,如文献[1-4]等㊂ 然而该近 用 中 产 生较大 误差 问题 , 如 图 1 为 级 数 项 N 分 别 取 W (u)值 即处于 震荡 发 散 不 稳 定 状态 ,N 值 越 小 误差 越
8.63476
Hale Waihona Puke Baidu9.57332
8.57333
0.00108
5
2.4 契比雪夫（Chebyshev）有理式逼近 的数值近似解㊂ 对于井函数 W (u)， Cody[7]给出的方法如下：
⎧ | | | | | | | | | | | | | | | | | | | | | ⎨ | | | | | | | | | | | | | | | | | | | | | ⎩
误差 小 于 1E-10 的 高 精 度 数 值 解 , 以 此 结果 作 为参 表 3㊂ 照 目 标 值 检 验 上 述 不同 井 函 数 计算 方法的 精 度 , 见
u 高精度 级数 展开 W (u) 切比雪夫 Srivastava 拟合 方程 式(6) 式(7) Barry N=5 N=30
表3 不同方法井函数计算精度对比 Table3 Accurate comparsion of various methods for well function 10.9357198 10.9357195 10.9357195 10.9357198 10.9357198 10.9357198 9.46056389 10.9357198 1E-5 4.03792958 4.03792924 4.03792924 4.03783522 4.03792958 4.03787773 4.03792515 4.03793067 0.01 0.21938393 0.21958956 0.21938360 0.21927152 0.21938393 0.21938351 0.21938407 0.221 1
⎞ | | | | | | | | | | | | | ⎠
u>4
﹙
其中：G =exp(-γ),γ=0.5772156649015328606 ；
ɿ 26 1 • q ,q = 20 • h= + hɕ u 1+ q 47 1+u •ɿ u 2.3 Srivastava 方法[6]
31
b=
ɿ
﹙
﹚
﹚
(3)
n2
2,i
i ðq 2,• i u i=0
• u-i
1<uɤ4
n3
(5)
图 1 级数展开方式的井函数计算示意图 Fig.1 Thewell function with various series expansion
exp(-u)• 1+ 1 • i=0 n u u
3
ðp ðq
i=0
3,i
• u • u
0.564514611673438
5 6 7 8
Fig.2
图 2 拟合替代公式井函数计算示意图 The regression function for well function
u
-0.0336932239452858
2.6 不同方法精度比较 Matlab ㊁1stOpt 等 已 内 置 了 指 数 积 分 函 数 W (u), 可得到 当 今 市面 上 的一 些 著名 数学数 值 计算 软 件 包 如
第3期
程先 云等：Theis 井 函数计算方法及井 模型 参数优 化计算研究
⎧ | | | | | | | | | | | ⎨ | | | | | | | | | | | ⎩
9 uɤ1
i i
2.1 级数展开
接得出其积分表达式,为便于数值计算,通常将其以无 穷级数形式展开,如下式： • u W (u)=-γ-ln(u)-ð (-1)n n • n! n=1
p 5 p• • p 8) u>1 1 exp(1+p• 2 u 3+p• 4 u +p• 6 ln(u))+p• 7 exp(u
uɤ1
(7)
10
4.3 3.8 3.3 2.8 2.3 1.8 1.3 0.8 0.3 -0.2
水 文
表2 拟合替代公式 (7) 系数表 Table2 The coefficient of regression function 序号 1 2 3 4
-0.57722 3.95850 0.2677
0
表1 Srivastava方法系数表 Table1 Coefficient of Srivastava method 0.99999 1 -0.24991 25.63296 18.05902 2 0.05519 3 -0.00976 1 1 4
21.09965
2
Theis 井 模 型
表示为：
定流量非稳定流承压水完整井 Thies 井模 型可以 S=
ɕ i
其中：W (u)=
显 而 易 见 ， 虽然国内外 还 有不 少 其 它 替 代 计算 方式 ， 但都缺乏较为全面和系统的对比介绍㊂ 第二点是模型 参数估算， 虽然自行编程亦 可实现 ， 但 由 于涉及 到 优 化和数值算法以及编程技巧， 对一般水文工作者而言 有相当难度， 而随着当今计算 软硬 件 的 快速 发 展 ， 利 用成熟优化建模计算平台不失为一可用之道， 这类平 台林林总总㊁ 性能各异， 有的优化能力强但无内置井 函数计算函 数 ， 如 Lingo㊁GAMS 等； 有的 井函 数为内 应用效果已被国内外各行业广大研究者认可㊂ 本研究 1stOpt 却可很好满足两者的要求， 使用上也非常简单， 置 函 数 ， 但 优 化 能力 却 相对 较 弱 ， 如 Matlab 等 ； 而
能通过数值计算方法近似求解㊂ 关于井函数 W (u)， 在数 称 为 Theis 井 函 数或井 函 数， 其 数 值 求解 方式国内外 学者已进行了众多的探讨和研究， 除了传统的误差较 主要可归纳总结如下㊂ 大且 实施繁琐 的配 线 法和 Jacob 直 线图 解 法外， 其 它
收稿日期:2014-04-26 基金项目:北京市西郊地区地下水战略储备关键技术研究与示范(Z141100003614060)；水利部公益性行业科研专项经费项目 山洪灾害监测预警 系统标准化研究 (201201058) 作者简介:程先云 (1965-),男， 四川宜宾人， 工学博士， 高级工程师， 主要从事水信息研究㊂ E-mail:xycheng@iwhr.com
数 W (u)的方法,其替代公式如下： 的误差㊂
2 • (1-G ) ,h = (1-G ) • (G 2-6 • G +12) ɕ G• (2-G ) 3 • G• (2-G )2 • b
如图 2 示,拟合式(6)当 u 大于 2 时 W (u)出现明显 笔者 基 于 最小 二乘原理 , 得出 如 下 井 函 数 W (u)替
程先云 1， 郑凡东 2， 杨
浩 2， 曹大岭 1,杨
勇2
优缺点及对模型最终参数计算的影响， 提出了新的简易计算方法； 以经典实际案例为对象， 对比研究 数优化计算问题， 并可广泛用于其它水文模型参数求解㊂ 关键词:Theis 井模型；Theis 井函数；参数优化；1stOpt 中图分类号:P345 文献标识码: A 应慎用广为采纳的级数展开井函数计算方式；(2)1stOpt 通用优化计算平台可很好处理 Theis 井模型参 文章编号:1000-0852(2015)03-0008-06 即采用该平台进行模型参数优化计算㊂ 了优化软件平台在 Thies 井模型参数计算中的实际应用效果；计算结果表明：(1)Thies 井模型参数估算
⎧ | | | | ⎨ | | | | ⎩
解几乎完全重合： W (u)=
代公式,具体系数见表 2,效果图见图 2,和高精度数值 b 1+b 2 • ub 3+b 4 • ln(u)+b 5 • exp(b 6 • u)
p
系数如表 1㊂
Srivastava 于 1998 年提出的方法如公式(4)示,其
第35卷第3期 2015年6月
水 文 JOURNAL OF CHINA HYDROLOGY
vol.35 No.3 Jun.， 2015
Theis 井函数计算方法及井模型参数优化计算研究
(1. 中国水利水电科学研究院， 北京 100038 ；2. 北京市水科学技术研究院， 北京 100044) 摘 要:归纳总结了 Theis 井模型中井函数的不同数值计算方法， 探讨﹑分析和比较了 各 方法的 利弊 ﹑
契比雪夫有理式逼近可用于高精度计算复杂函数
-ln(u)+
ðp •u ðq •u
i=0 1,i i=0 n1 1,i
n1
i
i
0<uɤ1
-1 0 1 2 u 3 4 5
N=9 N=11 N=10 N=8 N=6
• W (u)= exp(-u)
ðp
i=0 n2
⎛ | | | | | | | | | | | | | ⎝
0 2 4 6 8 10
第35卷
高精度 拟合式一 拟合式二
b
p
0.645517319421736 0.220523369822589 -1.00000014554241 -0.637030075359308 -1.22273508901879 0.99922210236199
-0.189854016260644 -1.73604166880858 -0.234855517244515 -2.06712062991927 0.693543466102803 1.29172191628888
4.156969E-6 -2.624277E-5 4.156478E-6 4.156969E-6 4.156969E-6 4.112166E-6 0.7894027 100.17575
10
9.835525E-11 -52247.83580 4027.53816