地下水储量作为最为宝贵的自然资源之一，对人类社会生活、工农业发展以及环境保护具有重要作用。全球约有50%的饮用水和40%的工业用水来自地下水^[1]，开采地下水资源便成为维持经济建设和生活用水需求的关键途径。但是，地下水资源的过度开采将会引发地面沉降、水土流失和生态恶化等问题^[2]。在现有的地下水开采方案下，华北平原鲁北地区地面沉降在未来十年仍会继续发展^[3]。因此，找到一种合适的方式来预测地下水储量变化对地下水的合理开发与利用具有显著的现实意义。

地下水储量监测主要依靠实际水位观测数据，但是稀疏的地下水观测数据导致地下水储量的分析受限。重力恢复与气候实验（gravity recovery and climate experiment，GRACE）卫星观测数据与水文模型相结合^[4-5]，革新了地下水储量变化的监测、分析和预测方式。文献[6]从GRACE卫星反演得到的陆地水储量中扣除全球陆地数据同化系统（global land data assimilation system，GLDAS）模拟输出的土壤水和雪水，发现2003—2016年土耳其地下水月消耗速率约6 cm等效水厚（equivalent water thickness，EWT）。文献[7]则将GRACE卫星和WaterGAP全球水文模型（WaterGAP global hydrology model，WGHM）结合，探测到中国鄂尔多斯和山西地区在2003—2014年间，地下水每年分别减少约0.65 cm和1.16 cm。文献[8]将GRACE数据与水文模型进行融合，提高了对澳大利亚东南部和古尔本河流域地下水储量的估计精度。

目前，结合GRACE卫星和水文模型或实际观测值分析地下水储量变化的研究较多^[9-11]，但利用GRACE卫星相关数据预测流域地下水储量变化的研究较少。一方面，传统的地下水储量预测依赖地下水埋深数据、水文循环以及含水层地质构造等复杂信息，同时地下水储量的变化受到人类活动、水岩作用过程、土壤以及降雨调节等自然因素的共同影响^[12]。另一方面，传统的预测模型，如自回归（autoregressive，AR）模型，由于其线性性质导致地下水储量预测精度有限^[13]。人工神经网络的出现提高了地下水储量预测的可靠性。文献[14]使用一种机器学习的方法即支持向量机，验证了GRACE卫星预测地下水储量的可行性；文献[15]将降雨和温度作为驱动因子，采用带外部输入的非线性自回归（non-linear autoregressive with exogenous input，NARX）神经网络成功预测了3种不同类型含水层的地下水储量；文献[16]利用温度、降雨和GRACE反演的陆地水储量作为输入变量，证实了多层感知机（multilayer perceptron，MLP）神经网络可以提高对地下水水位的预测精度。

随着全球气候变化加剧，研究和预测流域地下水储量变化对调节区域水文循环和改善生态环境具有重要意义^[17]。尽管人工神经网络可以提高预测地下水储量的精度，但是不同的网络类型和网络结构具有不同的预测表现。MLP神经网络属于静态神经网络，预测能力有限；而NARX神经网络虽然具有反馈连接，但是外部输入或者驱动因子通常难以确定^[18]；非线性自回归（non-linear autoregressive，NAR）神经网络具有动态机制，同时无需外部输入，仅仅依赖目标序列进行预测^[19]，当外部输入未知或者驱动数据难以获得时，NAR神经网络体现出其优越性。

本文从流域降雨分布出发，提出一种季节调整和NAR神经网络结合的流域地下水储量预测方法，在提高地下水储量预测精度的同时，减少网络训练参数，从而加快网络的收敛速度。首先对流域的降雨和地下水储量进行季节调整，然后利用NAR神经网络对季节调整后的地下水储量进行预测，并与未经季节调整的NAR神经网络、AR模型以及季节性自回归差分移动平均（seasonal autoregressive integrated moving average，SARIMA）模型进行对比。以长江流域、勒拿河流域、鄂毕河流域和叶尼塞河流域为例，分析和验证地下水储量的预测结果。

1 方法和数据

1.1   降雨分布与地下水分布

如果某一流域降雨是独立分布的，满足关系式^[20]：

式中,$P_t$是流域降雨时间序列；${\mu }_P$是其均值；$t$为观测时段；${\epsilon }_t$是独立分布的时间序列,其均值为0,方差为${\sigma }_{\epsilon }^{\mathrm{2}}$。

流域降雨量与径流量和地下水储量关系如图 1所示，该流域当月的地下水储量$G_t$一部分来自地表渗透$a{P}_t$,$a$为渗透系数,另一部分来自上个月滞留的地下水储量$(\mathrm{1}-c)G_{t-\mathrm{1}}$,$c$为地下径流系数。考虑流域降雨均值${\mu }_P$与地下水储量均值${\mu }_G$的影响,则有：

图  1  流域水文循环模型（改自文献[20]）

Figure  1.  Environmental System Model of a Watershed (Modified from Reference [20])

下载: 全尺寸图片 幻灯片

将式（1）代入式（2）可得：

由式（3）可知，如果降雨是独立分布的，那么地下水储量服从一阶自回归模型，即AR（1）模型。

1.2   季节调整

由于流域降雨通常具有季节性变化特征，地下水储量必定也是季节性变化的，因此采用稳定季节性滤波的季节调整方法^[21]，将原始的降雨序列和地下水序列分别分解成非季节性分量和季节分量，结果如图 2、图 3所示。其中，地下水和降雨序列的时间跨度均为2003-01—2015-06，时间分辨率为月。该方法假定季节分量是稳定的，即每年的季节性变化是一致的，因此对经过季节调整后的地下水进行预测时，加回季节分量便得到原始地下水储量的预测。本文还利用单位根检验对经过季节调整后的地下水进行检验^[22]，结果表明，在显著性水平为5%的情况下，经过季节调整后的地下水序列为平稳序列。

图  2  原始降雨序列、非季节性分量和季节分量

Figure  2.  Time Series of Precipitation, Deseasonalized Component and Seasonal Component

下载: 全尺寸图片 幻灯片

图  3  原始地下水序列、非季节性分量和季节分量

Figure  3.  Time Series of Groundwater, Deseasonalized Component and Seasonal Component

下载: 全尺寸图片 幻灯片

1.3   模型评价指标

本文采用4种指标评价模型的预测能力^[18]。

1）均方根误差（root mean square error，RMSE），用来衡量预测模型的全局拟合度：

式中,$y_t\mathrm{、}{o}_t\mathrm{、}n$分别为模型输出值、实际观测值及其样本数量。

2）相关系数R,用来衡量模型的预测能力：

式中,$\overline{y}$和$\overline{o}$分别为模型输出值和实际观测值的均值。

3）标准化均方根误差$R^{\mathrm{*}}$：

式中,${\sigma }_o$为观测值标准方差。

4）纳什系数（Nash-Sutcliff efficiency, NSE）：

根据文献[18]，笔者认为：当NSE > 0.75，$R^{\mathrm{*}}$<0.5时,模型表现为优；当NSE>0.65,$R^{\mathrm{*}}$<0.6时,模型表现为良；当NSE>0.5,$R^{\mathrm{*}}$<0.7时,模型表现为中；否则模型表现属于一般。特别地,当NSE取值为负时,模型预测能力不如均值预测。

1.4   数据

本文使用的降雨数据取自全球降雨测量（global precipitation measurement，GPM）数据产品，作为热带降雨测量任务（tropical rainfall measuring mission，TRMM）卫星的继代品^[23]，时间分辨率为月，空间分辨率为0.1°×0.1°。地下水储量估计是从GRACE卫星反演的陆地水储量中扣除GLDAS水文模型模拟输出的水分量^[24]，即：

式中,$\mathrm{G}\mathrm{W}\mathrm{S}\mathrm{A}$（groundwater storage anomalies）是地下水储量变化；$\mathrm{S}\mathrm{M}\mathrm{S}\mathrm{A}$（soil moisture storage anomalies）、$\mathrm{S}\mathrm{W}\mathrm{E}\mathrm{A}$（snow water equivalent anomalies）和$\mathrm{C}\mathrm{W}\mathrm{S}\mathrm{A}$（cover water storage anomalies）分别是GLDAS水文模型模拟输出的土壤水、雪水和植物冠层水变化,空间分辨为1°×1°；陆地水储量（terrestrial water storage anomaly,$\mathrm{T}\mathrm{W}\mathrm{S}\mathrm{A}$）采用由喷气推进实验室（jet propulsion laboratory, JPL）最新发布的Mascon Level-3解,空间分辨为0.5°×0.5°。具体处理如低阶项替换、泄露误差改正等参见文献[25]。

2 结果和讨论

2.1   自相关函数和偏自相关函数

图 4和图 5分别为经过季节调整后的降雨和地下水的样本自相关系数和偏自相关系数。由图 4可以看出，虽然个别滞后次数的自相关系数超出置信区间（可能归咎于数据的不确定性和自相关函数计算的偶然性^[20]），但是大部分自相关系数取值均落于95%置信区间内，因此可以推定4个流域的非季节性降雨是独立分布的。由图 5可以看出，4个流域的自相关系数都呈现出衰减形式，而偏自相关系数都在滞后次数为1处显著自相关，因此可以推定4个流域的非季节性地下水储量变化均满足一阶自回归模型，即AR（1）模型。

图  4  非季节性降雨样本的自相关函数和偏自相关函数

Figure  4.  Sample Autocorrelation Function and Partial Autocorrelation Function of Deseasonalized Precipitation

下载: 全尺寸图片 幻灯片

图  5  非季节性地下水样本的自相关函数和偏自相关函数

Figure  5.  Sample Autocorrelation Function and Partial Autocorrelation Function of Deseasonalized Groundwater

下载: 全尺寸图片 幻灯片

为了核实AR（1）模型对4个流域的非季节性地下水储量进行拟合的正确性，本文对残差依次进行正态性检验、异方差性检验和独立性检验。正态性检验采用核密度估计，异方差性检验是通过计算平方残差自相关函数，独立性检验则是直接计算残差的自相关函数^[21]。各流域的检验结果类似，因篇幅所限，图 6仅展示勒拿河流域的检验结果。从图 6可以看出，残差满足正态、独立和方差一致性条件。因此，结合图 5结果可知，经过季节调整后的地下水储量满足AR（1）模型。

图  6  勒拿河流域非季节性地下水残差检验

Figure  6.  Residual Test of Deseasonalized Groundwater in the Lena Basin

下载: 全尺寸图片 幻灯片

2.2   地下水储量的预测

考虑到AR模型的线性性质导致其预测能力有限，因此结合季节调整技术和NAR神经网络对4个流域地下水储量进行预测，并与未经季节调整的NAR神经网络、AR（1）以及SARIMA模型的预测结果进行对比。两种NAR神经网络均采用单一隐藏层，隐藏层的激活函数为双曲正切函数，输出层的激活函数采用线性函数。经过季节调整后的NAR神经网络（记为NAR（后））的时延数根据§2.1结果确定为1，隐藏神经元的数目则通过1~24逐个递增的方式来确定不同流域的最佳神经元数目。而季节调整前的NAR神经网络（记为NAR（前））的神经元数目和时延数均采用1~24逐个递增的方式确定。输入数据按照70%、15%、15%划分为训练集、验证集和测试集，同时采用随机初始权值方法对每个NAR神经网络训练100次，取100次训练结果的平均值和2倍标准差分别作为最终的预测结果和对应的95%置信区间。NAR（后）与AR（1）模型都是首先对经过季节调整后的地下水储量进行预测，然后加上季节分量得到地下水储量的预测结果；SARIMA模型则是对原始序列进行季节性差分，季节性差分阶数和非季节性自回归阶数都为1。本文输入数据集是2003-01—2014-12的地下水序列，预测时间为2015-01—2015-06。

图 7是不同方法在不同流域预测结果对比（为了更直观地进行比较，横轴限制在2014-01—2015-06）。由图 7可以看出，在长江流域，NAR（后）神经网络表现明显优于NAR（前）、AR（1）和SARIMA模型；在勒拿河流域，NAR（后）和NAR（前）模型预测结果与实际观测值的符合程度优于AR（1）和SARIMA模型；在鄂毕河流域，NAR（前）神经网络和AR（1）模型表现明显不如NAR（后）和SARIMA模型；在叶尼塞河流域，SARIMA模型和NAR（前）在峰值处的符合程度优于NAR（后）与AR（1）模型。图 7还展示了NAR（后）预测结果的95%置信区间，可以看出，在显著性水平为5%的情况下，NAR（后）预测结果落在地下水实际观测值范围内。

图  7  地下水储量的预测结果对比

Figure  7.  Comparison of Groundwater Prediction Results

下载: 全尺寸图片 幻灯片

表 1统计了不同预测方法的4种评价指标，由表 1可以看出，AR（1）模型在长江流域、鄂毕河流域以及叶尼塞河流域表现一般，且在长江流域和叶尼塞河流域NSE取值为负，说明AR（1）模型预测能力在这两条流域内不如简单的均值预测。SARIMA模型仅在鄂毕河流域表现为优，而在长江流域和勒拿河流域表现一般。NAR（后）与NAR（前）神经网络预测结果的RMSE均在1 cm以内，小于目前GRACE卫星反演地下水储量的不确定性。NAR（后）神经网络在4个流域的表现结果均为优，相关系数超过0.96；而NAR（前）在长江和鄂毕河流域表现一般。在勒拿河和鄂毕河流域，尽管NAR（前）神经网络的RMSE值小于NAR（后），但是由于NAR（后）神经网络的神经元数目和时延数减少，显著提高了网络的收敛速度。

表  1  预测结果的评价指标

Table  1.  Performance Indexes of Prediction Results

预测方法	流域	RMSE/cm	R*	NSE	R	模型评价
AR（1）	长江	1.407	2.958	-9.502	0.732	一般
	勒拿河	0.717	0.469	0.736	0.964	良
	鄂毕河	1.831	0.703	0.406	0.973	一般
	叶尼塞河	1.790	0.963	-0.113	0.969	一般
SARIMA	长江	2.896	1.241	-0.848	0.936	一般
	勒拿河	1.320	0.864	0.105	0.862	一般
	鄂毕河	0.678	0.261	0.919	0.973	优
	叶尼塞河	0.910	0.490	0.712	0.941	良
NAR（前）	长江(19, 23)	1.267	0.859	0.114	0.940	一般
	勒拿河(19，14)	0.580	0.380	0.827	0.914	优
	鄂毕河(20，17)	2.971	1.142	-0.564	0.942	一般
	叶尼塞河(2，12)	0.562	0.302	0.890	0.947	优
NAR（后）	长江(8, 1)	0.480	0.392	0.815	0.966	优
	勒拿河(11, 1)	0.604	0.396	0.812	0.968	优
	鄂毕河(17, 1)	0.787	0.303	0.890	0.990	优
	叶尼塞河(2, 1)	0.810	0.436	0.772	0.973	优
注：流域名称括号内的第一个数字代表神经元数目，第二个数字代表时延数

下载: 导出CSV 

| 显示表格

勒拿河流域NAR模型训练过程如图 8所示。由图 8可以看出，在勒拿河流域，NAR（前）神经网络在第5个历元收敛，而NAR（后）仅需要一个历元就完成模型训练。因此，结合季节调整和NAR神经网络不仅有效提高了流域地下水储的预测精度，同时减小了网络训练参数和收敛时间。

图  8  勒拿河流域NAR模型训练过程

Figure  8.  Training Procedure of NAR Models in the Lena Basin

下载: 全尺寸图片 幻灯片

2.3   讨论

首先，本文提出的预测方法的前提是假定季节分量是稳定的，因为非季节性分量通常受到极端气候的影响。因此预测误差除了数据本身的不确定性外，还可能来自于季节分量的变化。其次，本文预测的时间长度为6个月，这是因为地下水储量变化受到人类活动、降雨以及各种复杂的气候变化影响，超过半年的预测可靠性降低^[26]。NAR神经网络预测的峰值略低于实际观测值（图 7），这是因为本文采用的损失函数是RMSE，而更合理的损失函数应该顾及误差的季节性特征^[20]，这将在未来的研究中持续改进。SARIMA预测结果在鄂毕河流域略优于NAR神经网络（表 1），这是因为SARIMA模型对原始序列进行了一阶季节性差分处理。这也间接证明了对原始序列进行预处理，如季节性差分或者季节调整等，可以提高流域地下水储量的预测精度。同时，季节调整处理也为NAR神经网络时延数的确定提供了有效途径，从而降低神经网络收敛时间。各模型在不同流域的表现不同，除了与模型结构和实验方式有关外，可能与流域地理气候和位置等特征有关，这将在未来的研究中继续完善。最后，从流域降雨分布信息出发，推导出诸如地表水、径流等分布信息，不仅具有明确的物理背景，也会提高人们对区域水文循环的理解，为水资源的合理利用提供决策支持。

3 结语

本文将季节调整技术和NAR神经网络结合，分别对长江流域、勒拿河流域、鄂毕河流域和叶尼塞河流域的地下水储量进行预测，主要结论有：（1）经过季节调整后的降雨和地下水分别满足独立分布和一阶自回归模型；（2）季节调整技术在有效确定NAR神经网络的时延数的同时，减少网络训练参数，可以提高对地下水储量的预测精度和收敛速度；（3）经过季节调整后的NAR神经网络的预测结果优于传统的AR模型和SARIMA模型，预测结果的RMSE小于目前GRACE卫星反演地下水储量的不确定性，相关系数超过0.96。

尽管地下水储量变化模式复杂，受到诸如气候变化和人类活动等影响，但是从流域水文循环的物理背景出发，经过季节调整后使得对地下水储量的预测变得简单、可行，同时也体现了神经网络在对水文时间序列预测的表现优于传统模型。