基于EMD-AR模型的黄河下游灌区地下水漏斗演变预测

刘中培，齐明坤，韩宇平，曹润祥，冷静

(华北水利水电大学 a.水资源学院；b.河南省黄河流域水资源节约集约利用重点实验室，郑州 450046)

地下水作为重要的水资源，在维持社会经济发展、生态环境和农业生产需求方面发挥着重要作用[1].地下水开采量的不断增大导致地下水漏斗的产生及发展，并在某些地区出现了难以恢复的永久性降落漏斗[2-4].地下水漏斗的发展会进一步引发诸如地面沉降、地裂缝、咸水入侵、地下水污染等环境地质问题[5-6]，影响区域生态保护和高质量发展[7].目前关于地下水漏斗的研究主要集中在漏斗演变特征及修复等方面，而对地下水漏斗演变趋势的预测研究相对较少[8-11].

人民胜利渠灌区位于黄河下游，是河南省重要的粮食生产区，作为重要的灌溉水源，地下水在农业发展中扮演着重要的角色[12].长期大量的开采地下水导致灌区地下水位普遍下降，出现了地下水降落漏斗并不断扩展.以往的研究主要集中在水资源管理、利用，以及优化配置等方面，有关灌区地下水漏斗方面的研究较少[13-15].本文采用地理空间分析方法识别了人民胜利渠灌区地下水漏斗的形成时间，分析了地下水漏斗的中心水位、面积演变特征[16-17]，构建了地下水漏斗的中心水位模型和面积模型对漏斗中心水位和面积进行预测，从而为实现地下水的科学管理与区域高质量发展提供依据.

1 研究区概况

人民胜利渠灌区(113°31′E～114°25′E，35°0′N～35°30′N)位于黄河北岸，总面积为1 486.84 km2.灌区属于温带大陆性季风气候，四季分明，春季干旱多沙、夏季炎热多雨、秋季秋高气爽、冬季寒冷干燥.多年平均气温约14 ℃，无霜期220 d左右.根据国家气象中心新乡站多年的降水资料，灌区多年平均降水量为620 mm，年内降水量分配不均，主要集中在6至9月份.多年平均潜在蒸发量为1 864 mm.

根据《河南省新乡县区域水文地质调查》成果，研究区地层可分为第三系上新统，以棕红黏土、粉质黏土为主；第四系下更新统，为冰积地层；第四系中更新统，为洪积地层；第四系上更新统下部，为冲积地层；全新统一段为冲积、洪积地层.全新统二段为冲积、洪积、风积地层.全新统三段为沼泽堆积地层.依据含水介质及类型，区域调查深度范围内地下水可划分为松散岩类孔隙水，半胶结碎屑岩类孔隙裂隙水两类.根据含水层埋深条件可分为浅层含水组、中深层含水组和深层含水组.按5 m的降深计算又可分为强富水区、中富水区、弱富水区.其中灌区西部关堤一带为强富水区(>3 000 m3/d)，灌区中部获嘉一带为中富水区(1 000～3 000 m3/d)，灌区北部共产主义渠附近为弱富水区(<1 000 m3/d).研究区概况见图1.

2 材料与方法

2.1 数据来源及处理

本文分析对象为浅层地下水，地下水埋深资料来源于河南省人民胜利渠管理局.收集了研究区内24个观测井(由于个别观测井经纬度数据缺失，某些年份观测井不到24个)1953-2017年逐月地下水埋深观测数据，地下水水位为观测井地面高程与埋深之差；地下水漏斗面积数据获取首先通过ArcGIS软件绘制研究区地下水等水位线图，然后识别地下水漏斗的范围，最后利用空间分析提取数据实现.

本文采用基于经验模态分解(Empirical Mode Decomposition,EMD)处理后的数据构建漏斗中心水位自回归(Autoregression,AR)模型，进行中心水位的预测；结合普通最小二乘(OLS)方法，建立漏斗面积与漏斗中心水位变化的OLS回归模型，预测漏斗面积.

2.2 研究方法

2.2.1基于EMD的AR模型

EMD是对一种信号进行平稳化处理，将信号中的不同尺度的波段或趋势项提取出来，从而得到一系列的不同特征尺度的经验模态函数(Intrinsic Mode Function,IMF)和一个残余项[18-19].基于EMD的AR模型建模步骤具体如下所示.

步骤1 提取原始时间信号的IMF和残余项.提取IMF过程中，IMF必须符合两个条件：I(t)极值的数量(最大值和最小值数量之和)与零穿越的数量必须相等或最多相差1；在I(t)的任意点，局部最大值定义的包络线的平均值和局部最小值定义的包络线的平均值应该等于零，提取过程可以通过“筛选”步骤来完成[20-21].原始时间序列D(t)表示为：

(1)

式中：D(t)表示原始时间序列，rn表示最终的残余项，Ii表示不同频率的IMF.

步骤2 对所有的经验模态函数IMFs分别建立AR模型[22-23].AR模型计算流程如图2所示，模型数学表达式为：

Ii,j=μ+φ1(Ii,j-1-μ)+φ2(Ii,j-2-μ)+…+

φp(Ii,j-p-μ)+εt,

(2)

式中：Ii,j为第i个IMF第j时刻序列值，μ为序列均值，φ1,φ2,…,φp为自回归系数，εt为白噪声项.

2.2.2平稳性及白噪声性检验

AR模型要求时间序列是平稳的非白噪声序列，因此需对EMD分解的各经验模态函数IMFs先进行平稳性检验再进行白噪声性检验.数据的平稳性采用单位根来进行检验，若不存在单位根，则时间序列平稳，若存在单位根，则时间序列不平稳[24].对于不平稳的时间序列采用差分处理，直到时间序列平稳，通过单位根检验[24-25].序列的白噪声性检验采用Ljung-Box Q检验的方法，用以检验时间序列是否存在自相关，若时间序列存在自相关，则序列为非白噪声性序列，否则为白噪声性序列[26].上述检验过程均由Eviews专业统计学软件完成.

2.2.3模型阶数和参数估计

模型阶数的确定是模型参数(φi,μ,εt)估计的前提条件，此外模型适用性检验的核心就是解决模型的定阶问题.本文采用赤池信息准则(Akaike Information Criterion,AIC)确定AR模型阶数[27].

(3)

(4)

2.2.4模型的诊断与修正

对各经验模态函数IMFs建立的AR模型进行残差序列相关检验和残差异方差检验，通过一系列检验来验证模型的准确性.对没有通过检验的模型进行修正，使其通过检验，以实现对模型的应用.本文采用拉格朗日乘数检验(Lagrange Multiplier Test,LMT)对残差序列相关性进行检验，采用怀特检验(White Test,WT)对残差异方差进行检验[28].

2.2.5OLS回归模型

本文使用OLS分析漏斗面积与中心水位的关系[29].在地下水漏斗演变过程中，漏斗中心水位和漏斗面积关联性较大，通过夏庄漏斗处的中心水位和漏斗面积的相关分析，二者呈现高度相关，且二次函数拟合时R2值最大，拟合程度最好.建立二者的回归方程，数学表达式为：

Y/m2=c+b(X/m)+a(X/m)2,

(5)

式中：Y为因变量(漏斗面积)，X为自变量(中心水位)，c为常数，a,b分别为一次项和二次项系数.

3 结果与分析

3.1 地下水漏斗演变分析

3.1.1漏斗形成及发展

为分析研究区内地下水漏斗的产生时间及演变特征，根据1953-2017年地下水位埋深数据，利用ArcGIS地统计学模块克里金(Kriging)插值法对地下水位进行空间插值，并绘制每年的地下水水位等值线图，地下水位变化明显的年份地下水水位情况如图3所示.

由图3(a)和图3(b)对比可知，1975年以前灌区未出现地下水降落漏斗(图3(a))，1976年在灌区的中部形成以董庄为中心的地下水降落漏斗，但地下水漏斗范围较小(图3(b)).随后漏斗面积逐渐扩大，1990年董庄漏斗面积达到最大，如图3(c)所示.图3(d)和图3(e)对比发现，2000年漏斗中心仍然在董庄(图3(d))，与1990年相比只是漏斗面积有所变化，而至2001年，以夏庄为中心的地下水漏斗开始逐渐形成，但此时董庄漏斗和夏庄漏斗基本相对独立，未形成统一的漏斗(图3(e)).至2003年，随着地下水位的持续下降，董庄漏斗和夏庄漏斗连通，形成了一个面积更大的地下水降落漏斗(图3(f)).图3(g)和图3(h)表明，随着时间的推移，降落漏斗中心逐渐转移至夏庄附近，并且漏斗面积持续扩大至研究区外围.

综上所述，研究区地下水漏斗演变大致经历了三个阶段.第一个阶段为1975年以前，灌区未出现明显的地下水降落漏斗，地下水位等值线分布相对均匀，地下水由灌区西南向东北方向流动.第二个阶段为1975-2000年，地下水漏斗开始出现，此阶段地下水漏斗中心在董庄附近.当地下水漏斗出现后，地下水流场发生了明显变化，形成了以地下水漏斗为中心的局部地下水流动系统.灌区地下水位等值线变得疏密不一，漏斗区等值线密集，水力梯度大，非漏斗区等值线稀疏，水力梯度相对较小.第三阶段为2001年至现在，随着地下水位的持续下降，漏斗面积逐渐增加，董庄漏斗和夏庄漏斗连通，逐渐形成了以夏庄为中心的更大的地下水降落漏斗，至2017年漏斗边界已扩展至研究区外围.

3.1.2漏斗中心水位及面积演变

漏斗中心水位及面积是衡量漏斗发展的重要表征，地下水降落漏斗中心及其水位变化反映了漏斗的位置及垂向演变特征；漏斗面积反映了地下水降落漏斗水平发展特征.研究区地下水漏斗演变特征如图4所示.

由图3和图4可知，1990年之前，董庄地下水水位不断下降，地下水漏斗面积逐年变大.1990年之后，漏斗中心水位逐年恢复，漏斗面积逐渐变小，整个变化过程较为平缓.1990-2000年董庄漏斗中心水位上涨约2 m，漏斗面积减小了约60 km2.2001年漏斗中心转移至夏庄，形成以夏庄为中心的地下水降落漏斗.与董庄相比，夏庄漏斗中心水位与漏斗面积的变化情况大致为两个阶段.2001-2005年二者变化幅度较大，漏斗中心水位下降了约5 m，漏斗面积增加了150 km2.2006-2017年二者变化幅度明显变小，漏斗中心水位下降了约2 m，面积增加了约5 km2.

综上所述，董庄漏斗变化过程与夏庄漏斗变化过程有所不同，董庄漏斗中心水位先下降后上升，面积先增大后减小，两者变化幅度相对较小，总体呈现负相关关系.夏庄漏斗中心水位和漏斗面积演变经历两个阶段，第一个阶段为2001-2005年，变化明显，水位迅速下降，面积快速增加；第二阶段为2005年之后，两者变化相对平稳，总体呈现明显负相关关系.

3.2 地下水漏斗模型构建及预测

3.2.1漏斗中心水位及面积模型构建

(1)漏斗中心水位预测模型构建

选取1953-2017年夏庄漏斗中心处的水位数据，利用EMD方法，将漏斗中心水位数据分解成4个经验模态函数IMFs，分解结果如图5所示.人民胜利渠灌区的地下水位变化过程是非平稳、非线性的，是由多种波动成分共同作用的结果，从中可以提取出4种不同波段的经验模态分量IMF和一个残余项：IMF1波动周期不明显，波动幅度变化比较大；IMF2波动周期与IMF1相比变化明显，波动周期为准10 a，波动幅度变化小；IMF3波动周期相比于IMF2变大，波动周期为准15 a，波动幅度小于0.5 m；IMF4波动周期最长，波动周期为准27 a，波动幅度在1983年以前变化最小，1983年后变化最大；残余项表明漏斗中心水位呈现逐年下降趋势.

对4个经验模态函数IMFs进行检验，检验结果均为平稳的非白噪声序列.由AIC准则确定IMF1,IMF2,IMF3,IMF4的AR模型分别为AR(2),AR(4),AR(6),AR(6).AR模型参数计算结果见表1.

表1 AR模型参数计算

采用拉格朗日乘数检验LMT,怀特检验WT，分别对残差序列进行相关性和异方差检验[30]，结果表明各经验模态函数IMFs的AR模型均通过检验.

选取1953-2002年作为模型率定期，2003-2017年作为模型验证期，模拟预测结果由表2可知，模拟预测结果与实际观测结果拟合程度较好，相对误差均小于5%，模型精度高，模拟效果理想，可利用此模型对未来年份漏斗中心水位进行预测.

表2 地下水位模拟预测结果误差分析

(2)漏斗面积预测模型构建

根据2001-2017年夏庄漏斗中心水位和漏斗面积数据，利用SPSS软件的OLS回归分析模块对中心水位和面积进行分析.通过计算确定了回归系数及模型相应参数，再通过比选，得到了地下水漏斗面积的回归方程.模型相应参数见表3.

表3 模型参数汇总

从表3中看出，二次函数建立回归模型时的R2值最大，为0.861，拟合程度最好，因此夏庄漏斗中心水位和漏斗面积回归模型如(6)式所示：

Y/m2=-3 234.814+114.865(X/m)-0.947(X/m)2,

(6)

式中：Y为漏斗面积，X为中心水位.

3.2.2漏斗中心水位及面积预测

利用EMD-AR模型对夏庄2022-2030年漏斗中心水位进行预测，预测结果如图6(漏斗中心水位曲线)所示，2022-2030年漏斗中心水位下降约1.5m，水位下降过程分为两个阶段.第一阶段为2025年以前，漏斗中心水位约以0.08m/a速度缓慢下降.第二阶段为2025年以后，漏斗中心水位约以每年0.25m/a速度下降.

基于水位预测结果，根据漏斗中心水位和漏斗面积回归模型对漏斗面积进行预测，预测结果如图6(漏斗面积柱状图)所示，2022-2030年漏斗面积增加约8.5km2，面积增加过程可分为两个阶段.第一阶段为2025年以前，漏斗面积约以1.4km2/a速度增加，增长速度较快.第二阶段为2025年以后，漏斗约以0.85km2/a速度增加，增长速度较为缓慢.

计算结果表明，在现有条件下，灌区漏斗中心水位持续下降，漏斗面积进一步扩大，亟须采取压采措施控制漏斗的发展.

需要指出的是，本文使用的EMD-AR方法明显的优点是所需资料相对不多(主要为长系列水位数据和漏斗形成后的面积数据)，目的是利用EMD对时间序列平稳化处理后的本征模函数IMF自身序列实现对未来情况的预测，是根据数据发展趋势进行的预测，属于统计模型的范畴，如果坡度、坡向、地面植被、开采政策等地下水影响因子不发生突变的话，数据系列具有一致性，用此方法预测效果较好.

4 结 论

(1)1976年以前，灌区地下水流场未发生较大变化，未出现地下水漏斗.1976年以后，灌区地下水流场开始发生改变，地下水漏斗最早出现在董庄，而后由董庄过渡，最终形成以夏庄为中心的地下水降落漏斗.在漏斗形成中，漏斗中心处的水位和漏斗面积在不断扩大，中心水位与漏斗面积变化呈现高度相关.

(2)2022-2030年夏庄地下水漏斗中心处水位持续下降.以2025年为分界点，2022-2025年漏斗中心水位下降缓慢，下降速度约为0.08m/a；2025-2030年漏斗中心水位下降迅速，下降速度约为0.25m/a.2022-2030年夏庄地下水漏斗中心水位下降约1.5m.

(3)随着漏斗中心水位的降低，漏斗面积也不断扩大.2022-2025年漏斗面积增长较快，增长速度约为1.4km2/a；2025-2030年漏斗面积增长缓慢，增长速度约为0.85km2/a.2022-2030年夏庄漏斗面积增长约8.5km2.