河流水沙调节影响下河水地下水相互转化滞后响应时间的确定方法

1.本发明涉及水文地质研究技术领域，尤其涉及一种河流水沙调节影响下河水地下水相互转化滞后响应时间的确定方法。


背景技术：

2.地下水是世界上最大的淡水资源，在维持人类健康、能源和粮食安全、社会经济发展、生态系统服务和生物多样性等方面发挥着重要作用。然而，受自然和人为因素的影响，如，降雨、土地利用、蒸发蒸腾、大坝运行、洪水频率和地下水抽取等，河流与周围地下水之间的水、物质(固体和溶质)和能量的横向和纵向交换发生了显著变化。特别是，河流水沙调节影响下河道径流和水位波动不仅会改变河道形态和河流的边界条件(河床和河岸)，亦能促进河流与地下水之间的水和溶质的运移。例如：河流水位的升高可以改变河流和地下水之间的流动方向，填充非饱和带，驱动地下水运动。这些变化进而影响河流及其沿岸生态系统的健康和安全。
3.为此，国内外针对河流水沙调节影响下地下水对河水的滞后响应展开了诸多有益的研究。然而，河水地下水的相互转化过程，受河流与监测井的距离、河流边界条件(半渗透和部分渗透河床)、河水和地下水之间的水头等多种机制影响，地下水对河流水沙变化的响应具有时间上的延迟或滞后，反之亦然。研究表明，这种滞后效应可能给水生和河岸生态系统带来严重的或潜在的生态后果，例如：春旱、植被退化、土壤和地下水盐碱化、地面沉降和荒漠化等。因此，在量化河流水沙调节影响下地下水与河水转化关系时，必须考虑其滞后效应。但是，现有描述滞后响应的技术方法主要为滑动平均法、广义线性模型、广义加性模型和脉冲响应函数等，此类方法多为定性描述河水地下水相互转化关系，难以定量化确定二者相互转化的滞后时长。


技术实现要素：

4.为了能够对地下水与河水相互转化的滞后响应时间进行量化，本发明提供一种河流水沙调节影响下河水地下水相互转化滞后响应时间的确定方法。本发明首次将分布式滞后非线性模型(distributed lag non-linear model，dlnm)引入到水文地质领域，并对dlnm的暴露-响应因子进行驯化，同时结合机器学习算法，以量化河水地下水相互转化的滞后响应时间，可为水文地质研究及其建模提供新的技术思路。
5.本发明提供一种河流水沙调节影响下河水地下水相互转化滞后响应时间的确定方法，该方法包括：
6.步骤1：收集与处理河流与地下水数据；
7.步骤2：基于极端随机树模型确定影响河道外地下水位变化的主要水沙因素；
8.步骤3:：对原始dnlm模型进行驯化，利用驯化后模型分析得到地下水位对河流主要水沙因素的滞后响应时间。
9.进一步地，步骤1具体包括：获取第一水文站和第二水文站的水沙参数数据，以及第一水文站至第二水文站之间的各个监测井的地下水位数据；所述水沙参数数据包括水位、径流量、含沙量和来沙系数数据。
10.进一步地，步骤2具体包括：
11.基于极端随机树模型计算水位、径流量、含沙量和来沙系数分别对地下水位变化的贡献率以确定影响河道外地下水位变化的主要水沙因素。
12.进一步地，步骤2具体包括：
13.步骤2.1：将水位、径流量、含沙量和来沙系数设置为特征值，将gwl设置为目标值；
14.步骤2.2：将河流与地下水数据分为训练集和测试集；
15.步骤2.3：在极端随机树分类器中，设置树为500棵，最大特征值为4，节点数为2，生成极端随机树模型；
16.步骤2.4：在生成的所述极端随机树模型中，通过对训练集和测试集的实测值的模拟分别预测训练集和测试集的预测值；
17.步骤2.5：将训练集和测试集的实测值和预测值相结合，计算auc值，如果auc《0.5，则返回至步骤2.3重新调整树的棵数、特征值和节点数；如果auc》0.5，则输出特征值对目标值的贡献率。
18.进一步地，步骤3具体包括：
19.步骤3.1：确定影响地下水位变化的主要水沙因素，联合多项式ploy、自然立方样条ns、最大滞后时间lag、横向特征维度argvar和纵向特征维度arglag建立交叉基函数crossbasis()；其中，设置最大滞后时间lag为6个月，自由度df为3和样条最大阶数degree为2；
20.步骤3.2：利用生成的交叉基函数crossbasis()将暴露-响应函数mkbasis()和滞后函数mklagbasis()联系在一起，通过对原始的河流与地下水数据进行特征转换生成新的数据；
21.步骤3.3：采用glm模型对生成的新的数据进行建模，通过设置目标值gwl，主要水沙因素和次要水沙因素生成dlnm模型；
22.步骤3.4：对生成的dlnm模型计算aic值：若aic值不是最小，则返回至步骤3.1，调整参数自由度df和样条最大阶数degree；若aic值最小，则执行步骤3.5；
23.步骤3.5：通过dlnm模型中的crosspred预测函数模拟预测生成目标值gwl对主要水沙因素的响应程度rr值。
24.进一步地，步骤3.2中，所述特征转换是指通过张量积将暴露-响应函数生成的基矩阵和滞后函数生成的基矩阵在交叉基函数的基础上生成新的交叉基矩阵。
25.进一步地，步骤3.3中，生成的dlnm模型如公式(2)所示：
26.e(y
t
)＝α β1w
tl
(wl) ns(time,df
×
year) ns(fr,df) ns(ssc,df) ns(isc,df) (2)
27.其中，t是观测日，y
t
为第t天地下水位，e(y
t
)为第t天预期地下水位；α是基线风险的截距，β1是w
tl
的系数向量，w
tl
是将平均水位应用于dlnm得到的交叉基函数，df是自由度，l表示滞后天数，time表示时间序列，ns表示自然立方样条，wl、fr、ssc和isc分别表示水位、径流量、含沙量和来沙系数。
28.进一步地，步骤3.5中，rr值是通过地下水受水沙调节影响的天数与水沙调节总天数的比率与地下水受非水沙调节影响的天数与非水沙调节总天数的比率的比值得出。
29.本发明的有益效果：
30.基于河流水沙数据和河道外监测井地下水位数据，通过dlnm模拟建坝前后地下水位对水沙变化的滞后响应。dlnm可以为解决现实操作困难提供一种重要的方法，而且操作便捷、耗时短、误差小、可控性高。此外，dlnm对数据限制较小，也能应用在水文领域内其他的滞后效应的研究，应用范围较广。
31.基于机器学习法确定影响地下水位变化主要因素，并联合dlnm分析了地下水位与河流水沙变化的滞后响应关系。但是，就该研究事件及其数据本身的特殊性，与以前的dlnm相比，本发明在驯化中通过直接量化最大滞后时间简化dlnm的运算过程，避免了数据在模型运算过程中的二次或再次处理，减少了运算时间和错误出现的概率，提高了模型的运算效率和准确性。
附图说明
32.图1为本发明实施例提供的河流水沙调节影响下河水地下水相互转化滞后响应时间的确定方法的流程示意图；
33.图2为本发明实施例提供的黄河下游地下水监测井和水文站的位置；
34.图3为本发明实施例提供的ert模型驯化流程图；
35.图4为本发明实施例提供的小浪底调水调沙中期和后期水位(wl)、径流量(fr)、含沙量(ssc)和来沙系数(isc)在建坝前后分别对监测井(yy1、zm18、fq12和kf10)地下水位的贡献率；a1、b1、c1和d1分别表示建坝前水沙变化对yy1、zm18、fq12和kf10地下水位的贡献率；a2、b2、c2和d2分别表示建坝后水沙变化对yy1、zm18、fq12和kf10地下水位的贡献率；
36.图5为本发明实施例提供的dlnm模型驯化流程图；
37.图6为本发明实施例提供的建坝前监测井yy1(a)、zm18(b)、fq12(c)和kf10(d)地下水位季节变化对河道水沙的滞后响应；ssc和wl分别表示含沙量和水位；rr(relative risk)即相对危险度，表示暴露因子(ssc和wl)对地下水位的影响强度；
38.图7为本发明实施例提供的建坝后监测井yy1(a)、zm18(b)、fq12(c)和kf10(d)地下水位季节变化对河道水沙的滞后响应；fr、wl和isc分别表示径流量、水位和来沙系数；rr(relative risk)即相对危险度，表示暴露因子(fr、wl和isc)对地下水位的影响强度。
具体实施方式
39.为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
40.实施例1
41.如图1所示，本发明提供一种河流水沙调节影响下河水地下水相互转化滞后响应时间的确定方法，包括以下步骤：
42.s101：收集与处理河流与地下水数据；
43.具体地，获取第一水文站和第二水文站的水沙参数数据，以及第一水文站至第二水文站之间的各个监测井的地下水位数据；所述水沙参数数据包括水位、径流量、含沙量和来沙系数数据。
44.s102：基于极端随机树模型确定影响河道外地下水位变化的主要水沙因素；
45.具体地，基于极端随机树模型计算水位、径流量、含沙量和来沙系数分别对地下水位变化的贡献率以确定影响河道外地下水位变化的主要水沙因素。
46.作为一种可实施方式，具体包括以下子步骤：
47.s1021：将水位、径流量、含沙量和来沙系数设置为特征值，将gwl设置为目标值；
48.s1022：将河流与地下水数据分为训练集和测试集；
49.s1023：在极端随机树分类器中，设置树为500棵，最大特征值为4，节点数为2，生成极端随机树模型；
50.s1024：在生成的所述极端随机树模型中，通过对训练集和测试集的实测值的模拟分别预测训练集和测试集的预测值；
51.s1025：将训练集和测试集的实测值和预测值相结合，计算auc值，如果auc《0.5，则回到步骤s1023重新调整树的棵数、特征值和节点数；如果auc》0.5，则输出特征值对目标值的贡献率。
52.s103：对原始dnlm模型进行驯化，利用驯化后模型分析得到地下水位对河流主要水沙因素的滞后响应时间。
53.具体地，本步骤具体包括以下子步骤：
54.s1031：确定影响地下水位变化的主要水沙因素，联合多项式ploy、自然立方样条ns、最大滞后时间lag、横向特征维度argvar和纵向特征维度arglag建立交叉基函数crossbasis()；其中，设置最大滞后时间lag为6个月，自由度df为3和样条最大阶数degree为2；
55.s1032：利用生成的交叉基函数crossbasis()将暴露-响应函数mkbasis()和滞后函数mklagbasis()联系在一起，通过对原始的河流与地下水数据进行特征转换生成新的数据；
56.s1033：采用glm模型对生成的新的数据进行建模，通过设置目标值gwl，主要水沙因素和次要水沙因素生成dlnm模型；
57.具体地，原始dnlm模型要求输入数据为时间序列数据，其基本模型如下：
[0058][0059]
其中，μ
t
＝e(y
t
)，g是链接函数，y
t
是t＝1,...,n的一系列结果。函数sj指定变量xj和线性预测变量之间的关系，由参数向量βj定义。变量uk包括其他具有由相关系数γk指定的线性效应的预测变量。βj和γk是相应的参数。
[0060]
驯化后生成的dlnm模型如公式(2)所示：
[0061]
e(y
t
)＝α β1w
tl
(wl) ns(time,df
×
year) ns(fr,df) ns(ssc,df) ns(isc,df)
ꢀꢀ
(2)
[0062]
其中，t是观测日，y
t
为第t天地下水位，e(y
t
)为第t天预期地下水位；α是基线风险的截距，β1是w
tl
的系数向量，w
tl
是将平均水位应用于dlnm得到的交叉基函数，df是自由度，l
表示滞后天数；time表示时间序列，用于控制长期趋势；ns表示自然立方样条，wl、fr、ssc和isc分别表示水位、径流量、含沙量和来沙系数。
[0063]
s1034：对生成的dlnm模型计算aic值：若aic值不是最小，则返回至步骤s1031，调整参数自由度df和样条最大阶数degree；若aic值最小，则执行步骤s1035；
[0064]
s1035：通过dlnm模型中的crosspred预测函数模拟预测生成目标值gwl对主要水沙因素的响应程度rr值。
[0065]
具体地，rr值是通过地下水受水沙调节影响的天数与水沙调节总天数的比率与地下水受非水沙调节影响的天数与非水沙调节总天数的比率的比值得出。
[0066]
实施例2
[0067]
1)数据收集与处理
[0068]
本案例选取黄河下游1983-2015年花园口水文站(hyk)和夹河滩(jht)水文站的水沙参数数据，包括水位(wl)、径流量(fr)、含沙量(ssc)和来沙系数(isc)数据；以及花园口水文站至夹河滩水文站之间的河道外的4个监测井yy1、zm18、fq12和kf10的地下水位数据。以上的水沙参数据和地下水位数据均从黄河水利委员会的水情日报网和黄河统计年鉴获得。通过将花园口水文站的水沙参数数据与yy1和zm18的地下水位数据匹配结合，以及将夹河滩水文站的水沙参数数据与fq12和kf10的地下水位数据匹配结合，分别分析了小浪底大坝运行前后黄河下游地下水位在不同时间尺度上对河道水沙变化的滞后响应。案例研究区如图2所示。
[0069]
2)确定影响河道外地下水位变化的主要水沙因素
[0070]
基于极端随机树模型(ert)计算所选择的水沙参数(wl、fr、ssc、isc)对地下水位变化的贡献大小，确定影响地下水位变化的主要水沙因素。如图3所示，具体过程如下：
[0071]
步骤a1：输入河流水沙和地下水数据；本实施例中，是指1983-2015年水沙(wl、fr、ssc和isc)和地下水位数据；
[0072]
步骤a2：将输入数据中的自变量wl、fr、ssc和isc设置为特征值x，因变量gwl设置为目标值y；
[0073]
步骤a3：将河流水沙和地下水数据按照7:3比例分别划分为训练集和测试集；
[0074]
步骤a4：在极端随机树分类器中，设置树为500棵，最大特征值为4，节点数为2，生成极端随机树模型ert；
[0075]
步骤a5：在生成的ert模型中，通过对训练集和测试集的实测值的模拟分别预测训练集和测试集的预测值；
[0076]
步骤a6：将训练集和测试集的实测值和预测值相结合，计算auc值，判断模型精准度：如果auc《0.5，模型需重新回到步骤a4调整树的棵树、特征值和节点数；如果auc》0.5，说明模型模拟结果具有可信度，此时将输出特征值x对目标值y的贡献率结果，并通过柱状图呈现。
[0077]
具体地，当模型的识别性能最优时，auc近似1，而auc值小于0.5则表示分类不准确。在本案例中，训练集的模型auc大于0.85，测试集的模型auc大于0.6，所构建的模型能够准确捕获水-沉积物变量与地下水位之间的关系。该方法的运行是在python2021中利用sklearn.ensemble中的extratreesregressor包实现的。
[0078]
运算结果如图4所示，图4为小浪底调水调沙中期和后期水位(wl)、径流量(fr)、含
沙量(ssc)和来沙系数(isc)在建坝前后分别对监测井(yy1、zm18、fq12和kf10)地下水位的贡献率。a1、b1、c1和d1分别表示建坝前水沙变化对yy1、zm18、fq12和kf10地下水位的贡献率。a2、b2、c2和d2分别表示建坝后水沙变化对yy1、zm18、fq12和kf10地下水位的贡献率。本案例是在汇总了花园口和夹河滩水沙分别对yy1、zm18、fq12和kf10的贡献率数据的基础上，利用origion2018软件组合绘制了图4所示的饼状图。
[0079]
3)滞后时长计算
[0080]
dlnm本质是研究暴露因子对其因变量的滞后响应关系，可推广到任何时间序列相关的暴露因子对因变量的滞后响应分析中。该模型以广义加性模型为基础，利用交叉基做基础，并通过组合2组函数解释暴露因子与因变量的非线性关系，阐述暴露因子导致的滞后效应。其中暴露-响应函数表达沿暴露因子维度的潜在的非线性暴露-响应曲线。另一个滞后响应函数代替时间维度。在dlnm基础上，为暴露和滞后维度选择多项式和样条函数，以提供足够的灵活性来捕捉暴露因子与其对因变量的延迟影响之间潜在的复杂联系。设置最大滞后时间lag，利用赤池信息准则(aic)调整模型参数自由度df，值越小则模型参数调整和模拟结果越好。同时，引入相对风险(rr，为暴露组和非暴露组的风险比)表征暴露因子对其因变量的响应强度。当rr》1时，暴露因子对其因变量产生响应，rr越大，响应越强。dlnm不仅可以确定暴露因子是否对自变量具有滞后响应，还能弥补确定滞后响应开始和结束时间方面的不足。而发明人认为水文地质学研究事实上正以长时间序列数据为基础，恰恰符合dlnm对数据的要求，因此本发明基于dlnm并对其进行驯化，进而厘定地下水对河流水沙变化的滞后时长。本发明可为水文地质学滞后效应研究提供新的方法。
[0081]
本发明在确定水沙调节前后影响地下水位变化的主要水沙特征的基础上，对dlnm的基本模型(公式1)进行驯化。其基本模型如下：
[0082][0083]
其中μ
t
＝e(y
t
)，g是链接函数，y
t
是t＝1,...,n的一系列结果。函数sj指定变量xj和线性预测变量之间的关系，由参数向量βj定义。变量uk包括其他具有由相关系数γk指定的线性效应的预测变量。βj和γk是相应的参数。
[0084]
如图5所示，驯化过程如下：
[0085]
步骤b1：输入河流水沙和地下水数据；
[0086]
步骤b2：确定影响地下水位变化的主要因素ssc，联合多项式ploy、自然立方样条ns、最大滞后时间lag、横向特征维度argvar和纵向特征维度(滞后维度)arglag建立交叉基函数crossbasis()。其中设置最大滞后时间lag为6个月，自由度df为3和样条最大阶数degree为2。
[0087]
步骤b3：已经生成的交叉基函数将暴露-响应函数mkbasis()和滞后函数mklagbasis()联系在一起，通过对原始河流水沙和地下水数据进行特征转换形成新的数据；
[0088]
具体地，特征转换是指通过张量积将暴露-响应函数生成的基矩阵和滞后函数生成的基矩阵在交叉基函数的基础上生成新的交叉基矩阵。
[0089]
步骤b4：以glm模型为基础对转换生成的新的数据建模，通过设置目标值gwl，主要因素ssc和次要因素fr、wl和isc生成dlnm。
[0090]
步骤b5：通过已经生成的dlnm计算aic值。如果aic值不是最小，则重新回到步骤b2调整参数自由度df和样条最大阶数degree。
[0091]
步骤b6：当输出的aic值最小，通过dlnm中的crosspred预测函数模拟预测生成目标值gwl对主要因素ssc的响应程度rr值，并通过2d图输出结果，如图6和7所示。其中rr值是通过地下水受水沙调节影响的天数与水沙调节总天数的比率与地下水受非水沙调节影响的天数与非水沙调节总天数的比率的比值得出。
[0092]
图6为建坝前监测井yy1(a)、zm18(b)、fq12(c)和kf10(d)地下水位季节变化对河道水沙的滞后响应。ssc和wl分别表示含沙量和水位；rr(relative risk)即相对危险度，表示暴露因子(ssc和wl)对地下水位的影响强度。图7为建坝后监测井yy1(a)、zm18(b)、fq12(c)和kf10(d)地下水位季节变化对河道水沙的滞后响应。fr、wl和isc分别表示径流量、水位和来沙系数；rr(relative risk)即相对危险度，表示暴露因子(fr、wl和isc)对地下水位的影响强度。
[0093]
本发明弥补了河水地下水相互转化的滞后效应用中的不足，明确了地下水位对河流水沙是否产生滞后响应，以及明确产生滞后响应的起止时间。
[0094]
本发明并没有像以往的研究直接选定经验认定的主要影响因素作为dlnm的主要影响因素。而是在ert方法的基础上，根据各个黄河下游各个水文站水沙和监测井的地下水位的具体数据分析了其各自的主要影响因素。运用了具体问题具体分析的原则。这个改变弥补了dlnm以往运用过程中经验和主观选定主要影响因素的不足，为数据滞后效应的合理化分析奠定了基础。
[0095]
最初dlnm模型通常需要定义极端值或者通过设置不同滞后时间长度来观察滞后效应。而本发明本身就以特殊事件(小浪底大坝调水调沙)发生后的地下水位变化对河道水沙变化的滞后响应为研究对象。所以，为了适应该研究内的水文数据特征，方便对地下水位与水沙变化响应关系的观察与分析，在模型驯化过程中本案例将dlnm模型进行了简化。该模型在本发明不再定义极端值和不同滞后时间，而是直接通过观察最大滞后时间段内不同滞后时间点的地下水位与水沙特征变化的响应程度，判断地下水位对河流水沙变化的整体滞后响应情况。这个改变简化了模型，避免了对数据的二次或再次处理。同时减少运算过程、时间和错误出现概率，能提高模型的运算效率和准确性。
[0096]
最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。