直线边界非稳定流抽水试验水文地质参数智能计算方法与流程

1.本发明属于水文地质参数计算领域，具体涉及一种直线边界非稳定流抽水试验水文地质参数智能计算方法。


背景技术：

2.在直线边界附近进行的非稳定流抽水试验，边界可分为直线隔水边界与直线补给边界，水文地质参数的计算可采用直线图解法或标准曲线配线法，应用配线法时，根据井函数绘制的标准曲线不是单一特定的曲线，而需要根据抽水井和观测孔的位置制定特定的标准曲线，标准曲线的计算、绘制及配线计算十分麻烦，有些情况下仅存在理论可行性，人工配线计算的人为误差也较大，因而，其实际应用受到很大影响。
3.针对上述问题，本发明人提出了采用人工智能优化配线计算水文地质参数的想法，计算中，利用人工智能理论模拟人工配线，不仅计算精度高，而且免除了繁杂的井函数计算、标准曲线的绘制、配线的过程，使得该方法的应用可以更加普及、实用。
4.对于智能优化配线来讲，优化理论仅是一种计算方法，所述智能优化算法为差分进化、粒子群算法、模拟退火算法、遗传算法、免疫算法，本示例中选用的智能优化理论为遗传算法和粒子群算法。


技术实现要素：

5.为解决现有技术中存在的问题，本发明提供了直线边界非稳定流抽水试验水文地质参数智能计算方法，免除了以往采用人工进行复杂的井函数计算、绘制标准曲线、人工配线的繁琐过程。
6.本发明采用的技术方案是，直线边界非稳定流抽水试验水文地质参数智能计算方法，包括以下步骤：
7.s1：根据直线边界附近的非稳定流抽水试验，获取观测孔降深
‑
时间(s
‑
t)观测数据，在对数坐标系下绘制s
‑
t试验曲线；
8.s2：根据边界条件、观测井位置的不同，在相同模数的对数坐标下，绘制标准曲线，根据最小二乘法原理，构建试验曲线与试验曲线配线拟合时最小离差平方和目标函数；
9.s3：初始化智能优化算法中的参数，采用智能优化算法求解配线离差平方和目标函数的最小值，完成试验曲线与标准曲线的最优拟合；
10.s4：标准曲线为单一曲线时、计算试验曲线在最优拟合时的纵、横移动距离a、b，当标准曲线为曲线簇时，计算试验曲线在最优拟合时的纵、横移动距离a、b及抽水井至隔水，或补给，边界的距离d；
11.s5：在试验曲线与标准曲线最优拟合时，随机选取一个匹配点，根据匹配点在两个坐标系中的值，根据公式计算水文地质参数。
12.本发明直线边界非稳定流抽水试验水文地质参数智能计算方法的有益效果如下：
13.可以直接给出完整干扰井非稳定流抽水试验计算水文地质参数的具体公式，在生
产实践中，不需要水文地质人员自行推演计算公式，该方法具有统一性和普遍性。
附图说明
14.图1为背景技术中直线隔水边界映射理论示意图。
15.图2为背景技术中直线隔水边界抽水试验示意图。
16.图3为直线隔水边界配线法示意图。
17.图4为本发明直线隔水边界抽水试验优化配线示意图。
18.图5为示例1中标准曲线及s
‑
t观测数据曲线。
19.图6为示例1中第10次优化迭代计算配线结果。
20.图7为示例1中第20次优化迭代计算配线结果。
21.图8为示例1中第50次优化迭代计算配线结果。
22.图9为示例1中第500次优化迭代计算配线结果。
23.图10为示例1中第1000次优化迭代计算配线结果。。
24.图11为示例1中适应度函数进化曲线。
25.图12为示例2补给边界抽水试验示意图。
26.图13为直线补给边界配线法示意图。
27.图14为示例2标准曲线及原始数据曲线。
28.图15为示例2中第3次优化迭代计算配线结果。
29.图16为示例2中第10次优化迭代计算配线结果。
30.图17为示例2中第50次优化迭代计算配线结果。
31.图18为示例2中第100次优化迭代计算配线结果。
32.图19为示例2中第300次优化迭代计算配线结果。
33.图20为示例2中适应度函数变化曲线。
具体实施方式
34.下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。
35.标准曲线配线法求解水文地质参数
36.(1)theis井函数w(u)
37.进行承压完整井非稳定流抽水试验时，theis推导出的井函数w(u)为：
[0038][0039]
其中，w(u)为theis井函数；t为导水系数；t为自抽水开始到计算时刻的时间；r为观测孔至抽水孔的距离；s为含水层的贮水系数；e
i
(
‑
u)为指数积分函数。
[0040]
井函数w(u)没有解析解，为了计算方便，通常将w(u)展开成级数形式：
[0041][0042]
(2)井函数w(u)的近似解
[0043]
井函数w(u)没有解析解，大部分采用级数展开式的近似解，适用范围具有一定的局限性，为此，很多科学工作者对此展开了研究，下面仅列出两种近似解。
[0044]
swamee and ojha(1990)提出了井函数w(u)的全域解近似计算公式：
[0045][0046]
研究表明，与w(u)表列值的比较表明，在整个参数值范围内，公式的计算精度可以控制在1％以内(rajesh srivastava，1995)。
[0047]
allen(1954)and hastings(1955)提出了w(u)的近似计算公式(rajesh srivastava，1995)：
[0048]
如果0＜u≤1
[0049][0050]
|ε(u)|＜2
×
10
‑8[0051]
式中，a0＝
‑
0.57721566；a1＝0.99999193；a2＝
‑
0.24991055；a3＝0.05519968；a4＝
‑
0.00976004；a5＝0.00107857；b0＝0.2677737343；b1＝8.6347608925；b2＝18.0590169730；b3＝8.5733287401；c0＝3.9584969228；c1＝21.0996530827；c2＝25.6329561486；c3＝9.5733223454。
[0052]
直线隔水边界附近的单井非稳定流
[0053]
标准曲线配线法求解水文地质参数
[0054]
(一)标准曲线配线法原理
[0055]
抽水井位于直线隔水边界附近时，观测孔降深应等于抽水量为 q的实际抽水井和边界相对称的另一侧，具有同样流量 q的虚井共同影响的叠加。根据叠加原理，观测孔的降深s是实井引起的降深s1和虚井引起的降深s2之和如图1所示。
[0056][0057]
式中
[0058]
式中，r1为观测孔至实井的距离；r2为观测孔至虚井的距离；d为抽水井至隔水边界的距离。
[0059][0060]
[0061]
令为距离比，记r＝r1，u＝u1，则则u2＝β2×
u
[0062]
公式1可改写为：
[0063]
构建井函数φ(u,β)：φ(u,β)＝[w(u) w(u2)]＝w(u) w(β2u)
[0064]
公式1可改写为：
[0065]
类似于theis标准曲线的情况，给定不同的β值，可以绘制出ф(u,β)
‑
1/u标准曲线。
[0066]
当观测孔位于抽水井和边界相垂直的轴线上图1(b)时，若观测孔至实井的距离r1＝r，则r2＝2d
±
r，
±
号根据观测孔的位置确定，则
[0067][0068]
由此可知公式1的方括号内的数值为u1和r/d的函数，这样，只要已知u1和r/d的值，即可计算出u2，通过查井函数表即可得到w(u2)的值，因而可计算出w(u1) w(u2)的值，如果用引进一个新的井函数则公式1可改写为公式2。
[0069][0070][0071]
(二)计算步骤
[0072]
与常用的配线法类似，在透明双对数坐标纸上作降深—时间(s,t)实际资料曲线，与和标准曲线进行配线拟合，读出两张坐标纸上配合点的坐标值s，t，带入下面公式求出水文地质参数。
[0073][0074][0075]
[0076][0077]
非稳定流求参数的方法，一般采用直线图解法和配线法。直线边界附近的井流试验，无论是隔水边界还是二提到的补给边界，应用配线法时，井函数不是单一特定的曲线，要根据抽水井和观测孔的位置制定特定的标准曲线，因而，其应用受到一定影响。
[0078]
(三)举例说明，实施例1
[0079]
有一砂岩组成的承压含水层，顶板为粘土质砂岩层，底板为粘土岩。根据勘探资料知底板向西南方向迅速抬高，构成隔水边界(图2)。抽水井的直径为0.4m。抽水进行了18d，稳定抽水量为1140m3/d。抽水井和观测孔近似位于与边界垂直的同一轴线上，观测孔距抽水井的距离r＝20m。抽水过程中观测孔的水位降深资料列于表1。
[0080]
表1抽水试验观测数据
[0081][0082]
(1)根据题意，稳定抽水量q＝1140m3/d＝47.5m3/h，r＝20m。
[0083]
一般情况下，试验人员不能确定抽水井距隔水边界的确切距离d，需要绘出大量的标准曲线进行配线计算，配线完成后，根据配线采用的具体标准曲线的r/d值，再根据公式6反算d值，为计算方便，本次计算不再画出过多的标准曲线，d采用直线图解法计算出的距离，d＝159m，r/d＝0.1258。
[0084]
(2)计算参变量u2，根据题意，观测孔位于抽水井的另一侧，远离隔水边界，则：
[0085][0086]
(3)在透明双对数坐标纸上作降深—时间(s
‑
t)实际资料曲线，与和标准曲线进行配线拟合，如图3所示。
[0087]
(4)选取配合点(本例选用的配合点是观测时间t＝50min，沉降s＝1.93m的点)读出两张坐标纸上配合点的坐标值[s]，[t]，如下：
[0088]
[s]＝1.93m，[t]＝50min，
[0089]
(5)代入公式计算导水系数t和贮水系数s分别为：
[0090][0091][0092]
数据可以看出，配线法的计算结果与本例题直线图解法计算结果基本一致。
[0093]
遗传算法(geneticalgorithm，ga)优化配线求解水文地质参数
[0094]
目标函数的建立
[0095]
人工智能配线是采用人工智能算法，模拟人工配线求取水文地质参数，构建目标函数，用人工智能算法求解目标函数最小值是本方法的重要步骤。
[0096]
配线法的实质是在双对数坐标系中，以[w(u1) w(u2)]
‑
1/u标准曲线作为标准值，以抽水试验(t
i
,s
i
)观测值作为原始数据，移动抽水试验观测点(t
i
,s
i
)向[w(u1) w(u2)]
‑
1/u标准曲线靠拢，使得观测点(t
i
,s
i
)向[w(u1) w(u2)]
‑
1/u标准曲线做最优拟合，观测点(t
i
,s
i
)沿着横轴移动距离为10
a
，沿着纵轴移动距离为10
b
，当拟合的离差平方和为最小时，观测点与标准曲线达到最优拟合。为此，可以根据最小二乘法的原理，构建配线离差平方和最小值目标函数f
v
＝minf(a,b)。对于曲线簇的标准曲线，构建配线离差平方和最小值目标函数f
v
＝minf(a,b,d)，则属于高维度目标函数优化问题。
[0097]
设有直线隔水边界附近，进行定流量非稳定流抽水试验，获得n组观测数据为(t
i
,s
i
)，i＝1，2，3，
…
，n。井函数w(u)可由近似计算公式计算或调用matlab内置w(u)井函数直接计算获得，当抽水试验观测值(t
i
,s
i
)与[w(u1) w(u2)]
‑
1/u标准曲线最优配线拟合时，根据最小二乘法原理，可以构建配线离差平方和最小值目标函数。
[0098]
(1)最小二乘法基本原理
[0099]
根据最小二乘法原理构建试验曲线与标准曲线配线拟合时离差平方和目标函数最小值：
[0100][0101]
其中，y
i
为第i次测量时的真实值；f(x
i
)为选取的模型得到的y的预测值；y为模型理论值；n为观测试验数据的个数；
[0102]
(2)构建目标函数
[0103]
如果抽水井到隔水边界距离d和观测孔至实井的距离r为已知值，且观测井位于隔水边界和抽水井之外，根据最小二乘法基本原理，则目标函数可采用公式7：
[0104]
[0105]
如果观测井位于隔水边界和抽水井之间，则目标函数可采用下面公式8：
[0106][0107]
如果抽水井到隔水边界距离d是未知量，观测孔至实井的距离r为已知值，则目标函数稍微复杂一些，未知量变为三个(a，b，d)，同样，目标函数也可采用上面的公式，需要进行多元函数寻优计算，迭代优化完成后，可以计算出s
‑
t曲线配线移动参数a和b以及抽水井到隔水边界距离d。
[0108]
(二)水文地质参数智能优化配线计算
[0109]
本文采用遗传算法(ga)，实现隔水边界抽水试验，降深—时间(s
‑
t)实际资料曲线与和标准曲线优化配线。
[0110]
遗传算法(ga)优化配线法计算步骤如下：
[0111]
(1)初始化：设置进化代数计数器t＝0，设置最大进化代数g，随机生成m个个体作为初始群体p(0)。
[0112]
(2)个体评价：计算群体p(t)中各个个体的适应度。
[0113]
(3)选择运算：将选择算子作用于群体。选择的目的是把优化的个体直接遗传到下一代或通过配对交叉产生新的个体再遗传到下一代。选择操作是建立在群体中个体的适应度评估基础上的。
[0114]
(4)交叉运算：将交叉算子作用于群体。遗传算法中起核心作用的就是交叉算子，从群体中随机选择两个个体，通过两个染色体的交换组合，把父辈的优秀特征遗传给子辈，从而产生新的优秀个体。
[0115]
(5)变异运算：将变异算子作用于群体。即是对群体中的个体串的某些基因座上的基因值作变动，以产生更优秀的个体。
[0116]
(6)终止条件判断：若t＝g，则以进化过程中所得到的具有最大适应度个体作为最优解输出，终止计算，选取匹配点，输出a、b、d，f
v
的计算结果、绘制配线图等图件，输出导水系数t和贮水系数s，否则返回(3)继续搜索。
[0117]
隔水边界遗传算法(ga)智能优化配线法计算实例分析
[0118]
为对比分析，计算实施例1的数据。
[0119]
(1)遗传算法(ga)输入的参数
[0120]
群体规模np＝10、交叉概率p
c
＝0.6、变异概率p
m
＝0.1、遗传运算终止进化代数g＝1000。
[0121]
(2)构建目标函数
[0122]
如果抽水井到隔水边界距离d和观测孔至实井的距离r为已知值，且观测井位于隔水边界和抽水井之外，则目标函数可采用公式7：
[0123][0124]
(3)设计编程
[0125]
①
打开matlab程序，设定隔水边界遗传算法(ga)计算目录路径；
[0126]
②
构建抽水试验观测数据文件data01.mat，形成两行(t
i
,s
i
)、m列(抽水试验观测个数)的数据文件；
[0127]
③
编制目标函数计算程序jscb.m，对a、b、d进行初始赋值。井函数w(u)的计算可采用swamee and ojha(1990)提出的井函数w(u)的全域解近似计算公式。如不采用经验公式，可利用matlab内置w(u)井函数直接调用求解；
[0128]
④
编制遗传算法(ga)隔水边界优化配线计算主程序yccl.m；
[0129]
⑤
编制变异操作getvariationmethod.m、交差操作getcrossmethod.m文件。
[0130]
目标函数jscb.m中，设置了allen(1954)and hastings(1955)提出的井函数w(u)近似计算公式和swamee and ojha(1990)提出的井函数w(u)的全域解近似计算公式，计算时可根据需要修改注释语句选用。
[0131]
本次计算采用的是swamee and ojha(1990)提出的井函数w(u)的全域解近似计算公式。
[0132]
(4)优化计算过程
[0133]
为直观了解遗传算法(ga)优化配线计算过程，运行matlab程序，在当前路径窗口选取遗传算法(ga)配线法计算路径，在当前文件夹窗口，找到遗传算法(ga)优化算法计算程序yccl.m，右键打开yccl.m，可在编辑器中查看、修改程序中的输入设置参数。本例中，群体规模n＝10，遗传运算终止进化代数g＝1000，交叉概率p
c
＝0.6，变异概率p
m
＝0.1。
[0134]
计算中修改了遗传运算终止进化代数g，以了解遗传算法(ga)迭代计算收敛情况，设置进化代数g＝10、20、50、500，1000，其他参数不作调整，g＝10，原始数据曲线已开始向右移动靠近标准曲线(图4～图5)，但计算结果不稳定，每次计算给出的配线离差平方和最小值f
v
都不同，f
v
＝125.9069607122196～1.044255764176940，变化较大。图5～图11为迭代计算过程中的配线状态曲线及适应度函数变化曲线，图中很清楚反映出ga算法快速迭代配线的过程，一般情况下，经过20步迭代计算，配线离差平方和最小值的均值(随机计算十次的均值)有所减小，经过50步迭代计算后，计算结果更加接近最终结果，配线离差平方和最小值的均值进一步缩小，经1000次迭代计算完成后，配线离差平方和最小值基本稳定在f
v
＝0.9839176482603303，反映出该方法随着遗传迭代运算次数的增加，逐渐接近最优计算结果。由于是三元函数的寻优计算，迭代寻优配线计算用时略长。
[0135]
上面各图中未显示出大量的标准曲线，为便于观看，仅示例性的给出了r/d＝0.05、0.5、1.0及配线结束后的最优标准曲线。
[0136]
(4)计算结果
[0137]
采用遗传算法(ga)，计算中选取的匹配点坐标值如下：
[0138]
[s]＝1.930000，[t]＝50min，
[0139]
1000次迭代计算完成后，计算耗时8.7s，计算结果如下：
[0140]
导水系数t＝117.8283150640656m2/d；
[0141]
贮水系数s＝2.010393771504699
×
10
‑3；
[0142]
抽水井到隔水边界距离d＝149.67m；
[0143]
标准曲线匹配参数r/d＝0.1336255；
[0144]
配线离差平方和最小值f
v
＝0.9839176482603303；
[0145]
横向平移值a＝0.1403963787429783、纵向平移值b＝0.9160768545838138。
[0146]
例题1人工配线的计算结果为：
[0147][0148][0149]
抽水井到隔水边界距离d＝159m，配线标准曲线参变量r/d＝0.1258。
[0150]
对比实施例1人工配线的计算结果，遗传算法(ga)优化配线计算结果更符合实际情况，抽水井到隔水边界距离d是未知量，配线法的操作几乎不可完成，因此，采用人工智能配线使得这种计算方法成为可能，各参数的计算更加便利。
[0151]
为了解直接调用matlab内置w(u)井函数与采用经验公式的差别，仅改变一下井函数的调用方法，采用直接调用matlab内置w(u)井函数计算目标函数，其他参数不变，结果计算时间令人心烦意乱，一夜未出结果，后将遗传运算终止代数g从1000降至500，其他参数不变，经过漫长的631分钟才出结果，为了这次得来不易的计算，还是将计算结果显示如下，供大家参考：
[0152]
同样采用遗传算法(ga)，直接调用matlab内置w(u)井函数，参数设置为群体规模np＝10、交叉概率p
c
＝0.6、变异概率p
m
＝0.1、遗传运算终止进化代数g＝500。优化迭代计算完成后，选取匹配点的坐标值如下：
[0153]
[s]＝1.930000，[t]＝50min，
[0154]
500次迭代计算完成后，计算耗时631min，计算结果如下：
[0155]
导水系数t＝116.2090799547340m2/d；
[0156]
贮水系数s＝2.054800742292142
×
10
‑3；
[0157]
抽水井到隔水边界距离d＝148.9m；
[0158]
标准曲线匹配参数r/d＝0.1343142；
[0159]
配线离差平方和最小值f
v
＝0.9919952173046122；
[0160]
横向平移值a＝0.1454970187012802、纵向平移值b＝0.9225340293849773。
[0161]
究其原因，还是在计算隔水边界井函数的过程中，涉及到多项合并调用内置积分，大部分时间均消耗在内置数值积分上，所以，在不影响计算精度的情况下，采用swamee and ojha(1990)提出的井函数w(u)的全域解近似计算公式是最好的解决方法。
[0162]
直线补给边界附近的单井非稳定流抽水实验
[0163]
标准曲线配线法求解水文地质参数
[0164]
标准曲线配线法原理
[0165]
当边界是一个水位固定的补给边界时，边界的影响等于实际抽水井和边界另一侧与抽水井对称的、具有流量
–
q的虚构注水井共同作用的叠加。因为此时虚井的流量是
–
q，则
[0166][0167]
式中
[0168]
式中，r1为观测孔至实井的距离；r2为观测孔至虚井的距离，类似隔水边界的推导：
[0169][0170]
令则u2＝β
′×
u1[0171]
因此，只要给出观测井的具体坐标值(x,y)、抽水井至隔水边界的距离d，已知u1时，即可计算出u2的值，可构建井函数φ(u1,β
′
)：
[0172]
φ(u1,β
′
)＝[w(u1)
‑
w(u2)]＝w(u1)
‑
w(β
′
u1)
[0173]
公式9可改写为：
[0174]
当观测孔位于抽水井和边界相垂直的轴线上(图1b)时，若观测孔至实井的距离r1＝r，则r2＝2d
±
r，
±
号根据观测孔的位置确定，则：
[0175][0176]
引进一个新的井函数
[0177][0178]
式9可以改写为：
[0179]
(二)计算步骤
[0180]
与常用的配线法类似，在透明双对数坐标纸上作降深一时间实际资料曲线，与和标准曲线进行配线拟合，读出两张坐标纸上配合点的坐标值s，t，带入公式求出水文地质参数。
[0181][0182][0183]
(三)实施例2
[0184]
某半无限承压含水层，在距河边190m处打一抽水井，在抽水井到河边的垂直线段上设一观测孔，该孔距河185m，现以29.5m3/h的抽水量进行抽水试验(如图12)，在观测孔中测得的水位降深如表2所示，计算含水层的导水系数和贮水系数(邹立芝，1991)。
[0185]
表2抽水试验观测数据
[0186][0187]
(1)根据题意，稳定抽水量q＝29.5m3/h，r＝5m，为方便看图，未画出过多的标准曲线，d＝190m，r/d＝0.0263；
[0188]
(2)计算参变量u2，根据题意，观测孔位于抽水井与补给边界之间，则：
[0189][0190]
(3)在透明双对数坐标纸上作降深—时间(s
‑
t)实际资料曲线，与和标准曲线进行配线拟合，如图13；
[0191]
(4)选取配合点(本实施例选用的配合点是观测时间t＝90min，沉降s＝1.37m的点)读出两张坐标纸上配合点的坐标值s，t，如下：
[0192]
s＝1.37，t＝90，
[0193]
(5)代入公式计算导水系数t和贮水系数s分别为：
[0194][0195]
[0196]
配线离差平方和最小值f
v
＝0.66134。
[0197]
基于随机权重粒子群(randwpso)优化配线求解水文地质参数
[0198]
(一)目标函数的建立
[0199]
人工智能配线是采用人工智能算法，模拟人工配线求取水文地质参数，构建目标函数，用人工智能算法求解目标函数最小值是本方法的重要步骤。
[0200]
配线法的实质是在双对数坐标系中，以[w(u1)
‑
w(u2)]
‑
1/u标准曲线作为标准值，以抽水试验(t
i
,s
i
)观测值作为原始数据，移动抽水试验观测点(t
i
,s
i
)向[w(u1)
‑
w(u2)]
‑
1/u标准曲线靠拢，使得观测点(t
i
,s
i
)向[w(u1)
‑
w(u2)]
‑
1/u标准曲线做最优拟合，观测点(t
i
,s
i
)沿着横轴移动距离为10
a
，沿着纵轴移动距离为10
b
，当拟合的离差平方和为最小时，观测点与标准曲线达到最优拟合。为此，可以根据最小二乘法的原理，构建配线离差平方和最小值目标函数f
v
＝minf(a,b)。对于曲线簇的标准曲线，构建配线离差平方和最小值目标函数f
v
＝minf(a,b,d)，则属于高维度目标函数优化问题。
[0201]
设有直线补给边界，进行定流量非稳定流抽水试验，获得n组观测数据为(t
i
,s
i
)，i＝1，2，3，
…
，n。井函数w(u)可由近似计算公式计算或调用matlab内置w(u)井函数直接计算获得，当抽水试验观测值(t
i
,s
i
)向[w(u1)
‑
w(u2)]
‑
1/u标准曲线最优配线拟合时，根据最小二乘法原理，可以构建配线离差平方和最小值目标函数。
[0202]
(1)最小二乘法
[0203]
根据最小二乘法原理构建试验曲线与标准曲线配线拟合时离差平方和目标函数最小值：
[0204][0205]
其中，y
i
为第i次测量时的真实值；f(x
i
)为选取的模型得到的y的预测值；y为模型理论值；n为观测试验数据的个数；
[0206]
(2)构建目标函数
[0207]
如果抽水井到隔水边界距离d和观测孔至实井的距离r为已知值，且观测井位于隔水边界和抽水井之外，则目标函数可采用公式：
[0208]
如果观测孔位于抽水井与补给边界之间，则目标函数为公式13：
[0209][0210]
如果观测孔位于补给边界与抽水井之外，则目标函数为公式14：
[0211][0212]
如果补给边界方位已知，抽水井到补给边界距离d是未知量，观测孔至实井的距离r为已知值，则目标函数未知量为三个(a，b，d)，不同于两个未知数的计算，目标函数需要采用上面的公式，进行多元函数寻优计算，迭代优化完成后，可以计算出s
‑
t曲线配线移动参
数a和b以及抽水井到补给边界距离d，确定配线标准曲线的配线参数r/d值。对于补给边界，一般情况下，抽水井到补给边界的距离d是已知的。
[0213]
(二)randwpso算法优化配线法计算步骤
[0214]
将标准pso算法中惯性权重ω设定为服从某种随机分布的随机数，这样一定程度上可克服ω的线性递减所带来的不足。
[0215]
我们采用随机权重粒子群(randwpso)优化配线算法，实现补给边界降深
–
时间(s
‑
t)实际资料曲线与和标准曲线进行优化配线。
[0216]
randwpso算法优化配线法计算步骤如下：
[0217]
(1)设置粒子数目n、学习因子c1、c2、随机权重平均值的最大值mean_max、随机权重平均值的最小值mean_min、随机权重的方差sigma、最大迭代次数m、粒子的维数d；
[0218]
(2)随机初始化种群中各粒子的位置和速度，计算井函数w(u)及适应度值；
[0219]
(3)评价每个粒子的适应度，将当前各粒子的位置和适应度值存储在各粒子的p
best
(个体最优位置)中，将所有p
best
中适应度值最优个体的位置和适应度值存储于g
best
(群体最优)中；
[0220]
(4)更新权重其中，μ
max
、μ
min
为随机权重平均值的最大值和最小值，n(0,1)表示标准正态分布的随机数，rand(0,1)表示0到1之间的随机数，σ为随机权重的方差。
[0221]
(5)用下式更新粒子的速度和位移：
[0222]
v
id
＝ω
×
v
id
c1r1(p
id
‑
x
id
) c2r2(g
id
‑
x
id
)
[0223]
x
id
＝x
id
v
id
[0224]
(6)对每个粒子，将其适应度值与其经历过的最好位置作比较，如果较好，则将其作为当前的最好位置；
[0225]
(7)比较当前所有p
best
和g
best
的值，更新g
best
；
[0226]
(8)若满足停止条件(通常为预设的运算精度或迭代次数)，搜索停止，选取匹配点，输出配线坐标纵横a、b移动值、配线离差平方和最小值f
v
，绘制原始曲线图、配线过程图、最终优化配线图，输出导水系数t和贮水系数s，否则返回(4)继续搜索。
[0227]
基于randwpso算法的补给边界降深—距离优化配线法计算实例
[0228]
为方便对比分析，计算采用实施例2的数据(表2)。
[0229]
(1)randwpso算法输入的参数
[0230]
粒子数目n＝100、学习因子c1＝c2＝2.05，随机权重平均值的最大值mean_max＝0.8、随机权重平均值的最小值mean_min＝0.5、随机权重的方差sigma＝0.2、最大迭代次数m＝300、粒子的维数d＝2。
[0231]
(2)构建目标函数
[0232]
根据题意，观测孔位于抽水井与补给边界之间，则目标函数为：
[0233][0234]
如果补给边界方位已知，抽水井到补给边界距离d是未知量，观测孔至实井的距离r为已知值，则目标函数未知量为三个(a，b，d)，不同于两个未知数的计算，目标函数需要采用上面的公式，进行多元函数寻优计算，迭代优化完成后，可以计算出s
‑
t曲线配线移动参数a和b以及抽水井到补给边界距离d，确定配线标准曲线的配线参数r/d值。对于补给边界，一般情况下，抽水井到补给边界的距离d是已知的。
[0235]
根据题意，本例题目标函数采用公式13，由于抽水井距补给边界的距离d为已知，本例题变为二元函数寻优问题，计算相对简单。
[0236]
(3)设计编程
[0237]
①
打开matlab程序，设定补给边界随机权重粒子群randwpso计算目录路径；
[0238]
②
构建抽水试验观测数据文件data01.mat，形成两行(t
i
,s
i
)、m列(抽水试验观测个数)的数据文件。
[0239]
③
编制目标函数计算程序fitness03.m，对a、b进行初始赋值。井函数w(u)的计算可采用井函数w(u)的计算可采用swamee and ojha(1990)提出的井函数w(u)的全域解近似计算公式。如不采用经验公式，可利用matlab内置w(u)井函数直接调用求解。计算时可根据需要修改注释语句选用。
[0240]
④
编制随机权重粒子群randwpso计算主程序randwpso.m
[0241]
⑤
编制执行程序run0807.m。
[0242]
(4)优化计算过程
[0243]
为了解randwpso算法优化配线计算过程，运行matlab程序，在当前路径窗口选取补给边界randwpso配线法计算路径，在当前文件夹窗口，找到随机权重粒子群算法执行程序run0807.m，右键打开run0807.m，可在编辑器中查看、修改程序中的输入设置参数。粒子数目n＝50，学习因子c1＝c2＝2.05，随机权重平均值的最大值mean_max＝0.8,随机权重平均值的最小值mean_min＝0.5,随机权重的方差sigma＝0.2，最大迭代次数m＝300，问题的维数d＝2。
[0244]
计算中修改迭代次数m值的大小，以了解randwpso算法迭代计算收敛情况，设置m＝3、10、50、300，其他参数不作调整，迭代计算一开始，m＝3，试验曲线即快速搜索靠近标准曲线(图14～图15)，但计算不稳定，每次计算给出的配线离差平方和最小值f
v
都不同，f
v
＝0.768744944309303～26.152251199957075，变化较大。图15～图20为不同迭代次数下的配线曲线及适应度函数变化曲线，图中反映出randwpso算法快速迭代配线的过程，一般情况下，经过10步迭代计算，配线离差平方和最小值f
v
接近最优迭代计算结果，但不稳定，经过50步迭代计算后，计算结果与最终结果基本一致，m＝300，配线离差平方和最小值基本稳定在f
v
＝0.656016610316616，反映迭代计算收敛速度快、计算结果稳定。
[0245]
(5)计算结果
[0246]
经randwpso算法优化配线计算完成后，选取匹配点的坐标值如下：
[0247]
[s]＝1.37，[t]＝90，
[0248]
300次迭代计算完成后，计算耗时3.86s，计算结果如下：
[0249]
导水系数t＝319.7154961911451m2/d；
[0250]
贮水系数s＝9.173710220167313
×
10
‑4；
[0251]
标准曲线匹配参数r/d＝0.0263；
[0252]
配线离差平方和最小值f
v
＝0.656016610316616。
[0253]
实施例2人工配线的计算结果为(邹立芝，1991)：
[0254][0255][0256]
配线离差平方和最小值f
v
＝0.66134。
[0257]
对比实施例2人工配线的计算结果，randwpso算法优化配线计算结果符合更为精确、结果更为符合实际情况，采用人工智能配线使得这种计算方法成为可能，各参数的计算更加便利。

技术特征：
1.直线边界非稳定流抽水试验水文地质参数智能计算方法，其特征在于，包括以下步骤：s1：根据直线边界附近的非稳定流抽水试验，获取观测孔降深
‑
时间s
‑
t观测数据，在对数坐标系下绘制s
‑
t试验曲线；s2：根据边界条件、观测井位置的不同，在相同模数的对数坐标下，绘制标准曲线，根据最小二乘法原理，构建试验曲线与试验曲线配线拟合时最小离差平方和目标函数；s3：初始化智能优化算法中的参数，采用智能优化算法求解配线离差平方和目标函数的最小值，完成试验曲线与标准曲线的最优拟合；s4：标准曲线为单一曲线时、计算试验曲线在最优拟合时的纵、横移动距离a、b，当标准曲线为曲线簇时，计算试验曲线在最优拟合时的纵、横移动距离a、b及抽水井至隔水，或补给，边界的距离d；s5：在试验曲线与标准曲线最优拟合时，随机选取一个匹配点，根据匹配点在两个坐标系中的值，根据公式计算水文地质参数。2.根据权利要求1所述的直线边界非稳定流抽水试验水文地质参数智能计算方法，其特征在于，所述s3的智能优化算法包括以下子步骤：s31：构建目标函数；s32：选取人工智能计算方法，初始化智能优化算法中的参数，求解配线离差平方和目标函数的最小值，完成试验曲线与标准曲线的最优拟合，计算含水层水文地质参数。3.根据权利要求1所述的直线边界非稳定流抽水试验水文地质参数智能计算方法，其特征在于，所述智能优化算法包括蚁群算法、粒子群算法、模拟退火算法、遗传算法、差分进化算法、免疫算法和禁忌搜索算法。4.根据权利要求1所述的直线边界非稳定流抽水试验水文地质参数智能计算方法，其特征在于，所述边界为直线隔水边界时，单井非稳定流抽水试验的井函数为：当观测孔位于抽水井和边界相垂直的轴线上时，若观测孔至实井的距离r1＝r，则r2＝2d
±
r，
±
号则根据观测孔的位置确定，则井函数井函数其中，r1为观测孔至实井的距离；r2为观测孔至虚井的距离；d为抽水井至隔水边界的距离，w(u)为theis井函数；t为导水系数；t为自抽水开始到计算时刻的时间；r为观测孔至抽水孔的距离；s为含水层的贮水系数；e
i
(
‑
u)为指数积分函数；u)为指数积分函数；
标准曲线为试验曲线为s
‑
t，配线离差平方和最小目标函数为：如果抽水井到隔水边界距离d和观测孔至实井的距离r为已知值，且观测井位于隔水边界和抽水井之外，根据最小二乘法基本原理，则目标函数可采用公式：如果观测井位于隔水边界和抽水井之间，则目标函数可采用下面公式：如果抽水井到隔水边界距离d是未知量，观测孔至实井的距离r为已知值，则目标函数稍微复杂一些，未知量变为三个(a，b，d)，同样，目标函数也可采用上面的公式，需要进行多元函数寻优计算，迭代优化完成后，可以计算出s
‑
t曲线配线移动参数a和b以及抽水井到隔水边界距离d，水文地质参数计算公式为：水文地质参数计算公式为：水文地质参数计算公式为：水文地质参数计算公式为：式中，t为导水系数，s为贮水系数，q表示抽水井流量，m表示含水层厚度。5.根据权利要求1所述的直线边界非稳定流抽水试验水文地质参数智能计算方法，其特征在于，所述边界为直线补给边界时，单井非稳定流抽水试验的井函数为：当观测孔位于抽水井和边界相垂直的轴线上时，若观测孔至实井的距离r1＝r，则r2＝2d
±
r，则：式中，d表示抽水井到隔水边界距离，r表示观测孔至实井的距离；井函数：
标准曲线为试验曲线为s
‑
t，配线离差平方和最小目标函数为：如果抽水井到隔水边界距离d和观测孔至实井的距离r为已知值，且观测井位于隔水边界和抽水井之外，则目标函数可采用公式：如果观测孔位于抽水井与补给边界之间，则目标函数为公式：如果观测孔位于补给边界与抽水井之外，则目标函数为公式：如果补给边界方位已知，抽水井到补给边界距离d是未知量，观测孔至实井的距离r为已知值，则目标函数未知量为三个(a，b，d)，不同于两个未知数的计算，目标函数需要采用上面的公式，进行多元函数寻优计算，迭代优化完成后，可以计算出s
‑
t曲线配线移动参数a和b以及抽水井到补给边界距离d，确定配线标准曲线的配线参数r/d值，对于补给边界，一般情况下，抽水井到补给边界的距离d是已知的，水文地质参数计算公式为：水文地质参数计算公式为：式中，t为导水系数，s为贮水系数，q表示抽水井流量，表示直线补给边界时，单井非稳定流抽水试验的井函数。
技术总结
本发明公开了直线边界非稳定流抽水试验水文地质参数智能计算方法，其包括直线隔水边界和直线补给边界，S1：根据直线边界附近的非稳定流抽水试验，获取观测孔降深


技术研发人员：张勇
受保护的技术使用者：张勇
技术研发日：2021.03.03
技术公布日：2021/6/29