基于VBA和群井叠加原理的水文地质参数反演和降深预测

刘 军

(水利部新疆维吾尔自治区水利水电勘测设计研究院，新疆 乌鲁木齐 830000)

水文地质参数获取是地下水资源评价中的重点之一，这些参数包括渗透系数、给水度和贮水系数等，而抽水试验是获取水文地质参数的最简单方法之一。将抽水试验按照抽水井数量可分为单井抽水和群井抽水试验；单井定流量抽水试验往往可采用稳定流的Dupuit公式或非稳定流的Theis公式计算含水层参数[1- 3]，在稳定流的单井抽水试验中(观测井小于1个)，需要引入影响半径的经验公式与Dupuit公式联立方程组才可计算出渗透系数，在非稳定流的单井抽水试验可通过Theis公式的Jacob直线图解法求取含水层渗透参数；而群井抽水试验通常采用数值模拟的方法进行参数反演及水位预测[4- 14]，但尚未较好地解决这一问题。

在基坑排水和水源地开采时往往是群井同时抽水，且各孔之间距离较近，这种情况下各井周围含水层的水文地质条件差异较小，可视含水层为均质的，因此各井的水文地质参数也相同，在这种条件下，可以将含水层概化为均质无限含水层。本文基于乌鲁木齐地铁2号线群井抽水试验结合叠加原理提出了一种参数反演的方法，并利用VBA编程预测群井附近含水层在任意时刻由于抽水导致的地下水位降深，再根据抽水前的初始流场即可得到抽水后任意时刻的流场。该方法的提出为水文地质参数获取及基坑排水等问题提供了新思路。

1 材料及研究方法

本文主要利用单井抽水试验和叠加原理提出了水文地质参数反演与水位预测的方法。具体内容如下：

1.1 单井抽水的Theis公式

在均质无限含水层的单井定流量抽水试验中，抽水井周围任意一点i的降深可用Theis公式[15]表示为：

(1)

式中，t—抽水试验时间，min；si—任意点i在t时刻的降深，m；Q—抽水井的涌水量，m3/d；T—含水层导水系数，m2/d；ri—任意点距抽水井的距离，m；μ*—贮水系数。

1.2 叠加原理

由线性偏微分方程和线性定解条件组成的定解问题，可以采用叠加原理[15]，而承压含水层的群井抽水试验过程中，地下水位在承压含水层顶板以上时，地下水运动的偏微分方程满足叠加原理的条件。当承压含水层中同时存在2口井同时定流量抽水时，如图1所示，抽水任意一点i处产生的降深是两口抽水井分别抽水时在i点处产生的降深叠加；将这一理推广到计算N口抽水井同时定流量抽水时在任意一点i处的降深，可表达为：

图1 群井叠加原理示意图[15]

(2)

式中，N—抽水井个数；QN—第N口抽水井的涌水量，m3/d；rN—第i点与第N口抽水井之间的距离，m。

2 参数反演与水位预测

2.1 参数反演方法探讨

单井抽水试验可采用Theis公式的配线法将观测孔数据与标准曲线拟合，从而计算含水层参数。根据这一特性，可将群孔抽水试验所有孔的观测降深与公式(2)的计算值拟合，从而反推导水系数T和贮水系数，这样能够很好的避免由于个别观测孔数据异常而对结果产生较大影响。参数反演的具体方法如下：

(1)确定抽水孔个数N及各抽水孔涌水量QN，绘制各观测孔的降深-时间散点图。

(2)预先给定任意大于零的导水系数T和贮水系数μ*(根据经验赋值能够提高参数反演效率)，根据公式(2)计算各观测孔(含有观测数据的抽水孔)的叠加降深，绘制随参数改变的动态降深-时间曲线。

(3)通过不断调整导水系数T和贮水系数μ*，以达到所有观测孔观测到的降深-时间散点图与计算的降深时间曲线达到最佳的拟合状态，此时导水系数T和μ*就确定了。在Excel中通过公式和函数功能即可实现观测孔曲线拟合，此处不再赘述。

采用这种方法拟合参数，个别观测孔的数据可能无法完全拟合，这是因为在拟合过程中各个观测孔的标准曲线是根据一个综合的平均参数(渗透系数和贮水系数)计算出来的，因此拟合出的参数也是整个含水层的平均值，避免因个别井的含水层结构突变而计算出的参数不具有代表性。

2.2 水位降深预测方法

根据拟合后的参数及公式(2)，可计算群井抽水试验中任意一点i(x，y)在t时刻的降深s(x，y，t)，因此只要给定一个范围，即可得到该范围内在任意时刻的降深或流场，其具体方法如下：

(1)确定一个所需要计算的流场范围，需要给出横坐标范围(xmax,xmin)和纵坐标范围(ymax,ymin)，设定一个水位计算点的间隔距离Δx和Δy，根据范围和间隔采用VBA代码提取出所有计算点的坐标值，计算点的间隔距越小插值出的降深场(流场)精度越高，但计算速度越慢。

(2)选定所需计算流场的时间点t，根据公示(2)采用VBA代码实现计算所有计算点i(x，y)在各时间点t时的降深值。

(3)采用ArcGIS软件可将各时刻的降深计算结果插值得到降深场，若需要绘制流场，将初始流场与降深场栅格化叠加即可得到各时刻的流场。

2.3 VBA程序实现

提取特定范围内的计算点坐标i(x，y)可采用VBA代码实现的具体代码见表1。

表1 提取计算点坐标的代码

提取出各计算点坐标后，给定一个抽水时间t，即可根据公式(2)在Excel中计算出该时刻各计算点的降深s(x，y，t)，通过插值软件即可获取任意时刻的降深场。

3 实例

3.1 含水层概况

本次研究选取了乌鲁木齐市轨道交通2号线某试验场地的承压含水层，该含水层平面面积40000m2，含水层厚度在20～30m，压力水头高出地表2～3 m，含水层岩性单一；上覆隔水顶板厚度在20～25m左右；下层有稳定的隔水底板，含水层周围无特殊构造，因此可将该含水层概化为均质、含水层水平无限延伸的承压含水层，含水层结构示意图如图2所示。

图2 含水层结构示意图

3.2 群井抽水试验概况

该含水层区域内有4口抽水井和4个观测孔，抽水期间4口抽水井同时抽水，抽水孔均采取定流量抽水，各孔的涌水量见表2。本次抽水试验时长为300min，4个观测孔同时观测水位并记录观测数据，各试验井(孔)参数见表2。

表2 抽水试验井(孔)参数表

3.3 参数反演

图3 观测孔降深拟合曲线图

由拟合结果可看出，4个观测孔的观测降深与群井叠加原理的计算降深拟合结果较好，表明该方法能够有效的反演水文地质参数，试算得到的导水系数和注水系数可视为该含水层的平均值，可用该参数进行流场预报。

3.4 水位降深预测

本次预报情景为4口抽水井同时开采，在每口涌水量均为2500m3/d的条件下，预报抽水30d内的水位变化情况，各观测孔30d水位降深变化图如图4所示，抽水后1、10、30d降深场如图5所示。由图4观测孔的降深变化曲线可以看出，随着不断地抽水，水位逐渐下降，水位下降速率逐渐变小。由图5抽水试验的不同时间点的降深场可以看出，抽水井群的降落漏斗连成一片，漏斗降深逐渐变大，等降深线逐渐向外扩散，这种特征满足Theis公式的特点，也符合实际的抽水情况。

图5 预测降深场(单位：m)

图4 观测孔降深变化曲线图

4 结语

结合VBA程序和水文地质学中的群井叠加原理提出一种反演水文地质参数的方法，将该方法应用于乌鲁木齐市地铁2号线某站点的承压含水沉群井抽水试验中，得到以下结论：

(1)基于群井叠加原理反演水文地质参数的方法能够有效利用所有观测孔数据，同时可避免异常孔的数据。该方法计算出的水文地质参数为某一区域含水层的平均值，适用于小范围内的基坑排水和水源地抽水条件下的求参。通过实例验证了群井叠加原理反演水文地质参数的可靠性。

(2)基于VBA程序能够有效提取设定范围和设定密度的计算点坐标，结合群井叠加原理能够预测不同抽水情景(定流量)下、不同时刻的降深场，通过对实例应用，预测结果较好。基于这一特征，还可将该方法用于解决“给定降深求解涌水量的问题”，为地下施工排水问题提供一种新的解决思路。

(3)群井叠加方法反演参数是基于Theis公式反演参数，因此该方法仅适用于水文地质结构较简单的情况，含水层需要满足Theis公式的假设条件。而在含水层结构复杂的条件下，应当采用数值法解决问题。

(4)与数值模拟相比，群井叠加的VBA方法预测流场操作简便，结合其它插值软件可绘制降深场或流场，在预测水位的过程中应当注意，当水位降至承压顶板以后，地下水运动不再遵循承压含水层的Theis公式，因此也不可使用于群井叠加方法预测水位降深。