引 言
河渠水位的变化是影响两岸地下水动态的重要因素．河渠间潜水的非稳定运动规律的研究，对区域地下水资源评价、灌排渠系设计、地下水动态预报都有着重要的意义[1-4]．其中，地表水体附近潜水非稳定流与渗流规律的研究[5-11]已在渠道附近地下水动态过程研究[12-16]中得到广泛的应用．
在之前的研究中，假设河渠完全切割潜水含水层，将河渠视为第一类边界的河渠-半无限潜水非稳定流模型，大体分为两类： 1) 河渠水位变化瞬时完成并保持稳定[1-4]，对于这类简单的边界条件函数，可利用常用的积分变换方法，求出模型的理论解； 2) 河渠水位持续变化[1,3,5-6]，对于这类问题多将水位的变化过程概化为一常规函数，如线性变化，当函数形式较为简单时，也可以通过常用的积分变化方法进行模型求解．对于上述两类问题现已有较为深入的研究，但实际情况中，河渠水位多难以一已知函数形式进行概化，即便某些特定条件下可给出变化过程的具体函数形式，但当函数过于复杂时[10]，其常用的积分变换方法便难以进行模型求解．而且在实际情况中，其河渠水位变化过程多为一未知函数，故对未知边界函数形式条件下的河渠-半无限潜水非稳定流模型问题进行相关研究具有实际意义．
针对上述情况，本文采用不考虑边界条件形式的Laplace变换方法求解，利用Laplace变换中的有关性质，给出模型的解；同时为了进一步讨论该理论解在实际中的运用，本文对河渠水位变化过程采用Lagrange插值进行离散化处理，并借助文献[17-18]中的算例数据进行验算，对其应用展开介绍．
1 经典模型及其解
经典的河渠附近潜水非稳定运动研究做出了如下假设：
① 潜水含水层均质各向同性，隔水底板水平、平面上无限延伸；
② 河渠完整切割含水层，河渠水位变动过程为函数H(t)；
③ 潜水初始水位h(x,0)水平，潜水流为一维流；
④ 忽略上部入渗量，无垂向交换作用．
上述水文地质概念模型，可写成
式中，μ是潜水含水层给水度；h是潜水水位,m；k是含水层渗透系数,m/d；H(t)是河渠水位随时间变化函数,m．H(t)=H0，H0为常数，即“河流水位起始瞬时升高H0后，水位保持不变”．
图1 河渠附近潜水的非稳定运动
Fig. 1 Transient phreatic motion near a river or a canal
当含水层厚度较大时，水位变化h(x,t)-h(x,0)≤hm时，根据Boussinesq方程第一线性化方法，模型
可改写成
式中，u(x,t)=h(x,t)-h(x,0)；a=khm/μ，a为潜水含水层的导压系数,m2/d．
当H(t)=H0时，利用Laplace变换，可获得上述模型的解：
(9)
式中，
是余误差函数．式(9)即为Ferris经典模型的解．
2 理论模型及其解
上述模型通常将河渠水位变化过程概化描述为一具体函数形式，但实际中的河渠水位变化过程H(t)，虽在小范围流域内有一定的规律，但由于外界干扰因素过多，其变化过程更具随机性，再加上人为干扰因素，基本无法对水位变化过程进行统一的具体函数表达．针对这一局限性，本文将河渠水位变化过程直接以通用的H(t)函数形式表达，在不追求H(t)的具体函数形式基础上，利用Laplace变换方法求解这类模型的理论解．
对模型(Ⅱ)求关于t的Laplace变换，获得
式中,
为u关于t的Laplace变换的象函数，s为Laplace算子，L和L-1分别为Laplace变换算符与逆变换算符．
对于模型
，当式(11)过于复杂或未知时，难以采用上述方法直接求解．对此，针对式(11)，不直接求H(t)的象函数，而只将其以算符的形式运用于计算过程中，再利用变换的相关性质和相应的特征函数进行求解．详解过程如下．
模型中式(10)的通解为
(13)
式中，λ1,λ2为待定常数；由边界条件(11)和(12)，模型的特定解为
(14)
值得指出的是，当H(t)未知时，
的Laplace逆变换求解比较困难，根据Laplace变换的“卷积定理”，对式(14)进行逆变换：
(15)
式中，*为卷积算符．
参照常用Laplace变换中的函数：
(16)
对式(16)，根据Laplace变换的“微分定理”，有
(17)
其中
结合式(15)、(17)，有
(18)
对式(18)，将卷积打开，注意h(x,t)和u(x)的关系可得
(19)
对于式(19)，本文对河渠水位变化过程以H(t)进行概化，而之前的研究中多以一已知常规函数进行处理，但在实际中，河渠水位的变化多难以函数形式进行概化，即便存在，也多为一些非常规函数，难以求解，而本文通过对Laplace变换性质的运用，避免了边界函数变换过程的复杂性，进一步推动了经典模型的解在实践中的运用．同时需要指出的是，上述求解过程中，边界条件H(t)虽未进行形式上的变换，但其实际参与了Laplace变换与逆变换的过程，所以H(t)应满足Laplace变换的要求，在实际中河渠水位变化基本是可以满足变化要求的．
3 运 用
3.1 Lagrange线性插值
实际中，显然难以给出复杂的水位变化过程H(t)统一的数学表达式；此时，可采用Lagrange线性插值方法[17]，对H(t)进行离散化处理．
对于时段Ti=ti-ti-1，有
(20)
则有
(21)
δ(t-ti-1)为Heaviside函数，具有如下性质：当t<ti-1时，δ(t-ti-1)=0；当t≥ti-1时，δ(t-ti-1)=1．
3.2 求参
河渠附近潜水非稳定流模型研究的重要目的之一，是以其为工具，利用地下水位动态数据计算含水层参数，也即模型中的参数．直接体现在式(19)中的参数有潜水含水层的导压系数a,m2/d．根据a和μ之间关系，还隐含潜水含水层渗透系数为潜水含水层的导压系数k，m/d，其中k=μa/hm．而一般而言，给水度μ可通过野外简单试验获得；利用非稳定流模型，主要是求导压系数a．
在之前河渠附近潜水非稳定运动的参数求解研究中，吴丹等[17-18]利用配线法和拐点法等传统理论方法对含水层参数进行求解，但由于上述方法对函数解析式构成形式要求较高，处理过程过于复杂，且需要大量的数据计算，在推广使用时具有较大的难度．针对上述情况，本文将利用MATLAB软件进行含水层参数的求解，为相关领域的实践工作提供借鉴．
在对水位变化过程进行插值后，不难发现若将式(21)直接代入模型的解中，其计算式中含有复杂的卷积计算，这会给曲线拟合求参带来较大难度．对此，本文在水位数据处理可行性和卷积定理微分特性的基础之上，考虑对潜水位变动速度和时间进行曲线拟合，已达到含水层参数求解的目的．
设潜水位变动速度φ(x,t)=∂h(x,t)/∂t，注意到u(x,t)=h(x,t)-h(x,0)，结合卷积定理的微分特性，由式(18)，有
(22)
结合式(21)，注意到当t=0时，H(t)=H0及Heaviside函数的性质，得
(23)
令
则上式变为
(24)
对于式(24)，当距离x已知后，即转化为含1个参数的非线性方程．为此，本文将利用MATLAB软件曲线拟合工具箱中的自定义函数拟合功能对上式进行优化拟合，确定含水层参数．
3.3 河渠与潜水之间交换量计算
1) 算式
河渠与潜水之间的补排关系，由河渠水位与潜水水位关系所决定，这也表明，河渠附近潜水非稳定渗流模型，可以用来计算河渠与潜水之间的水量交换．其中，河渠与潜水之间的交换量，可以用单位河渠长度上的交换量表达，即单宽流量q(t)；t时刻的交换强度q(t)，由Darcy定理：
(25)
如前文，在实际问题的运用中，还需将具体的H(t)代入到式中，进一步换算，才能得到河渠与潜水之间的水量交换计算公式．若H(t)采用Lagrange插值进行离散化处理，根据式(25)可得
(26)
式(26)即为河渠与潜水之间在t时刻的交换强度算式；q(t)>0，表示潜水接受河渠补给；q(t)<0，表示潜水排泄河渠．
2) λ段的影响规律
为简便，仅讨论H0，λ1的情形．由式(26)，去掉等式括号内的第2项，转化为Ferris模型的交换强度算式，对应的是河渠水位瞬间变动H0所获得的水量交换强度，记为qH；式(26)去掉等式括号内第1项，对应的是H0=0条件下、仅由河渠水位变化段所获得的交换强度，记为qλ． λ段的影响，可用qλ/qH来反映．由式(26)，有
qλ/qH=2λ1t/H0．
(27)
假设水位在变动段获得的水位升幅也为H0(即λt等于H0)，由式(27)，对交换量的贡献，λt是H0的2倍；也即，水位变动段的水位升幅虽然也为H0，但其对交换量的贡献，是不考虑变化过程影响而计算出效果的2倍；这是水位在t时段内缓慢上升至λt的过程中，被不断抬升的水位持续增强补给作用的过程累积效应所至．该累积效应也是经典模型边界概化方法无法反映的．
4 实 例 研 究
安徽淮北平原中部区域潜水含水层广泛发育，其下分布有不完全连续黏性土隔水层；潜水位埋深多较浅,且农业区内灌渠系统较完善，灌渠之间的间距多为2 km左右，干渠渠首大部分建设有节制闸；区内有一口国家级地下水位自记观测井，距离干渠距离约为60 m．
2014年7月24日12时至2014年7月25日12时，干渠节制闸关闸蓄水．根据水位观测数据，渠道水位变化过程可分为两个过程: 第一阶段,关闸20分钟内，水位迅速上升2.0 m;第二阶段,水位上升速度迅速减缓，至25日12时升幅为0.34 m．对应的地下水观测孔水位按3h摘录，如表1所示．
根据式(24)，第二阶段H0=2.0 m，λ=0.34 m/d．
表1 潜水水位动态数据
Table 1 The dynamic data of the water table
t/h36910111215182124h/m26.80126.81826.85726.87126.88726.9026.96227.02027.07927.135φ(x,t)/(m/h)0.000 20.005 80.012 90.014 70.016 10.017 20.019 00.019 60.019 40.018 9
注 初始水位为26.80 m，附近地面标高30.66 m．
Note The initial water table is 26.80 m and the ground level is 30.66 m．
图2 曲线拟合过程
Fig. 2 The curve fitting process
如图2，在MATLAB软件曲线拟合工具箱的custom equation选项下输入拟合函数，经计算，a值约为35.9 m2/h，即861.6 m2/d．求解结果与相关地区文献[19]中计算的值(854 m2/d)基本一致．
从图2中可以看出，曲线拟合程度较高，其计算结果精度较好，同时，可以看出该方法不涉及复杂的软件操作，对各种水位变化数据都有着较好的适用性，具有较好的推广价值．
5 结 论
通过上述模型求解和实际运用的研究，得到以下结论：
1) 针对半无限域河渠附近潜水非稳定运动经典模型中河渠水位边界条件概化的局限性，将河渠水位变化过程概化为通用函数形式，借助Laplace变换求解，可得到形式较为简单的理论解．此外，Laplace变换在模型求解中可有效避免边界条件函数变化所带来的象函数求逆的困难，对于未知边界函数的模型求解有着重要的意义．
2) 针对河渠水位变化过程，利用Lagrange插值的方法对河渠水位过程进行离散化处理；采用MATLAB软件曲线拟合工具箱中的自定义函数拟合功能进行含水层参数求解，参考相关文献，结果基本一致，且操作简易，对本专业实际工作和实验提供借鉴，具有较好的推广价值．
3) 河渠水位在其变动过程中，被不断抬升的河渠水位(以上升为例)，使河渠对地下水的补给作用呈持续增大并伴生有过程累积效应；在计算河渠与潜水之间水量交换时，采用经典模型的计算方法，将因忽略累积效应而导致计算偏差，偏差幅度为2λ1t/H0～λ1t/(H0+λ1t)．
4) 从数据处理的简便性出发，本文对水位变化过程只进行了简单的插值离散处理．随着计算机技术的发展，可以对长系列的水文数据进行处理，考虑到水位变化在实际情况中的随机性，需进一步研究各情况下理论解的实际应用．
5) 当实际水文地质条件较为简单时，在合理概化的基础之上可以直接采用解析法来求取水文地质参数并预测地下水动态．但需要注意的是解析法的局限性也非常明显，它不能适用于含水层的各种复杂条件，此外模型是建立在一维概化的条件下，因此将它利用于实际计算时不可能得到精确的结果，但作为一种进行理论研究和初步估算的工具来说，具有重要的意义．
[1] 张蔚榛. 地下水非稳定流计算和地下水资源评价[M]. 北京: 科学出版社, 1983.(ZHANG Weizhen. Calculation of Unsteady Flow of Groundwater and Evaluation of Groundwater Resources[M]. Beijing: Science Press, 1983.(in Chinese))
[2] 沙金煊. 农田不稳定排水理论与计算[M]. 北京: 中国水利水电出版社, 2004.(SHA Jinxuan. Theory and Calculation of Unsteady Drainage in Farmland[M]. Beijing: China Water & Power Press, 2004.(in Chinese))
[3] 薛禹群. 地下水动力学[M]. 北京: 地质出版社, 2010.(XUE Yuqun. Dynamics of Groundwater[M]. Beijing: Geological Publishing House, 2010.(in Chinese))
[4] 瞿兴业. 农田排灌渗流计算及其应用[M]. 北京: 中国水利水电出版社, 2011.(QU Xingye. Calculation and Application of Drainage and Irrigation Seepage in Farmland[M]. Beijing: China Water & Power Press, 2011.(in Chinese))
[5] 贾志峰, 李文宾, 贾志锐. 潜水非稳定渗流边界条件处理方法研究[J]. 灌溉排水学报, 2013, 32(3): 44-49.(JIA Zhifeng, LI Wenbin, JIA Zhirui. Boundary conditions treatment method in un-confined unsteady flow seepage process[J]. Journal of Irrigation and Drainage, 2013, 32(3): 44-49.(in Chinese))
[6] 张鸿雁. 河渠水位曲线回水影响半无限含水层河渠附近地下水非稳定流计算[J]. 长春地质学院学报, 1987, 17(3): 319-330.(ZHANG Hongyan. Unsteady groundwater flow calculation in a semi-infinite aquifer near a river or channel affected by backwater level in the river or channel[J]. Journal of Changchun College of Geology, 1987, 17(3): 319-330.(in Chinese))
[7] 阿里木·吐尔逊, 周志芳, 木塔力甫·依明尼亚孜. 河渠附近潜水非稳定运动的一种通解[J]. 河海大学学报(自然科学版), 2003, 31(6): 649-651.(ALIM Tursun, ZHOU Zhifang, MUTALIP Iminniyaz. A universal solution to unstable groundwater movement in vicinity of canals[J]. Journal of Hohai University(Natural Sciences), 2003, 31(6): 649-651.(in Chinese))
[8] 杨红坡, 谢新宇, 张继发, 等. 潜水一维非稳态运动的解析理论及应用[J]. 水科学进展, 2004 , 15(1): 82-86.(YANG Hongpo, XIE Xinyu, ZHANG Jifa, et al. Analytical solution of one-dimensional transient phreatic flow and its application[J]. Advances in Water Science, 2004, 15(1): 82-86.(in Chinese))
[9] 陶月赞, 席道瑛. 垂直与水平渗透作用下潜水非稳定渗流运动规律[J]. 应用数学和力学, 2006, 27(1): 53-59.(TAO Yuezan, XI Daoying. Rule of transient phreatic flow subjected to vertical and horizontal seepage[J]. Applied Mathematics and Mechanics, 2006, 27(1): 53-59.(in Chinese))
[10] LIANG X Y, ZHANG Y K. Analytic solutions to transient groundwater flow under time-dependent sources in a heterogeneous aquifer bounded by fluctuating river stage[J]. Advances in Water Resources, 2013, 8(58): 1-9.
[11] MAHDAVI A. Transient-state analytical solution for groundwater recharge in anisotropic sloping aquifer[J]. Water Resources Management, 2015, 29(10): 1-14.
[12] YOUNGS E G, RUSHTON K R. Dupuit-Forchheimer analyses of steady-state water-table heights due to accretion in drained lands overlying undulating sloping impermeable beds[J]. Journal of Irrigation and Drainage Engineering, 2009, 135(4): 467-473.
[13] SU N H. The fractional Boussinesq equation of groundwater flow and its applications[J]. Journal of Hydrology, 2017, 547(2): 403-412.
[14] BANSAL R K. Approximation of surface-groundwater interaction mediated by vertical stream bank in sloping terrains[J]. Journal of Ocean Engineering and Science, 2017, 2(1): 18-27.
[15] 张建锋, 李国敏, 张元, 等. 塔河下游间歇性输水河道附近地下水位动态响应[J]. 地球物理学报, 2012, 55(2): 622-630 .(ZHANG Jianfeng, LI Guomin, ZHANG Yuan, et al. Responses of groundwater levels to intermittent water transfer in the lower Tarim River[J]. Chinese Journal of Geophysics, 2012, 55(2): 622-630 .(in Chinese))
[16] 张学宏, 李颜, 郝培章, 等. 水文资料插值计算方法探讨[J]. 海洋预报, 2008, 25(1): 5-13.(ZHANG Xuehong, LI Yan, HAO Peizhang, et al. Discussion on the interpolation calculation methods of hydrological data[J]. Marine Forecasts, 2008, 25(1): 5-13.(in Chinese))
[17] 吴丹, 陶月赞, 林飞. 河渠水位线性变化条件下河渠-潜水非稳定流模型及其解[J]. 应用数学和力学, 2018, 39(9): 1043-1050.(WU Dan, TAO Yuezan, LIN Fei. Solution of the transient stream-groundwater model with linearly varying stream water levels[J]. Applied Mathematics and Mechanics, 2018, 39(9): 1043-1050.(in Chinese))
[18] 吴丹, 陶月赞, 林飞. 复杂函数边界控制下的潜水非稳定流模型及解的应用[J]. 水利学报, 2018, 49(6): 725-731.(WU Dan, TAO Yuezan, LIN Fei. Application of unsteady phreatic flow model and its solution under the boundary control of complicated function[J]. Journal of Hydraulic Engineering, 2018, 49(6): 725-731.(in Chinese))
[19] 陶月赞, 曹彭强, 席道瑛. 垂向入渗与河渠边界影响下潜水非稳定流参数的求解[J]. 水利学报, 2006, 37(8): 913-917.(TAO Yuezan, CAO Pengqiang, XI Daoying. Parameter estimation for semi-infinite phreatic aquifer subjected to vertical seepage and bounded by channel[J]. Journal of Hydraulic Engineering, 2006, 37(8): 913-917.(in Chinese))