华北山前平原典型井灌区地下水水位变化影响因素分析

张 婧，马贵宏，高 雅，沈彦俊，刘晓英，邹 妍，戴茂华

(1.河北省气象与生态环境实验室，河北 石家庄 050021；2.河北省气候中心，河北 石家庄 050021； 3.河北省石家庄市高邑县气象局，河北 石家庄 051330；4.河北省水文勘测研究中心，河北 石家庄 050031； 5.中国科学院遗传与发育生物学研究所农业资源研究中心，河北 石家庄 050021；6.天津长城滨银汽车金融有限公司，天津 300457； 7.天津师中师教育信息咨询有限公司，天津 300074；8.河北省农林科学院旱作农业研究所，河北 石家庄 053000)

地下水是世界上最大的淡水资源库，为全球提供了36%的饮用水和42%的灌溉用水，是农业、工业、环境和饮用水必不可少的供水来源[1-2]。华北平原是我国重要的粮食高产区之一，由于半干旱半湿润的气候条件无法满足两熟制(冬小麦和夏玉米)作物的需水量，且地表水资源较匮乏，农业地下水灌溉量占地下水总开采量的80%以上，人类活动已成为影响华北平原地下水环境变化的主导因素[3-4]。由于地下水资源的过度开采，1970—2008年华北平原地下水水位累计下降了8.74 m，山前平原地区地下水水位降幅最大，达到16.34 m，地下水漏斗连区成片引发了一系列地下水资源和环境问题[5-9]。随着华北地区社会经济的快速发展和城市人口的不断增长，城市工业和生活用水量将继续增加，水资源短缺问题将进一步加剧，会严重影响该地区社会经济和生态环境的可持续发展[4-7]。

在前人对华北平原地下水水位变化的归因分析中，降水量减少、气温升高、上游水利工程建设、农业种植结构改变和地下水开采等因素是造成该地区水资源减少的重要原因[10-11]。其中，人为开采尤其是农业灌溉是地下水水位下降的主要原因，灌溉引起的实际蒸散量高于降水量，导致该地区地下水资源持续减少[12-15]。华北山前平原小麦和玉米种植面积与农业地下水开采量关系密切，二者的增减变化趋势较为一致[14]，粮食产量增产的同时，作物消耗大量土壤水，造成补给地下水的入渗量大大减少，也是地下水水位连年下降的一个重要原因[2,15]。

除了人为影响因素外，气候变化尤其是降水量变化对区域水资源影响较为显著[13,16-19]。在影响华北山前平原含水层的气候因素中，除了降水量外，还有来自上游的地下水侧向补给[16-18]。山西作为华北平原地下水产流区和地下水补给区，自20世纪末以来降水量明显减少，导致地表产流和地下水补给量减少，成为太行山山前平原地下水水位下降的另一个主要原因[17]。降水量和气温的波动以及蒸散效应与地下水的补给和排泄类型具有一定的联系，华北平原地下水水位有2～7 a准周期变化，这与ENSO的波动周期较为一致[13,18-19]。一些学者从气候要素和人类活动两个方面开展了综合分析，认为气候变化影响人类活动，进而影响到地下水储量的变化。Taylor等[20]认为全球气候变化对地下水的间接影响大于其对地下水补给的直接影响，间接影响主要体现在气候变化通过影响灌溉需水量和灌溉来源间接影响地下水储量的变化。石家庄平原降水量影响农业地下水开采强度，干旱年份的地下水开采量明显多于湿润年份，高强度的地下水开采常常会掩盖气候变化在地下水动态变化中的作用[14,16,21]。因此，华北平原浅层地下水动态变化是气候与人类活动叠加效应的结果。

由于研究区和分析方法不同，得到的气候变化和人类活动对地下水动态变化影响的结论也存在差异。开展地下水动态影响定量研究的方法较多，主要有水文模型和数学统计方法两大类。利用MODFLOW、VFM等水文模型可以模拟地下水水位动态变化，而主成分分析、灰色理论等统计方法能较客观地分析地下水动态变化中各影响因子的响应。许月卿[22]利用投影寻踪回归技术，从降水、气温等自然因素和地下水超采、水利工程、农田产量、作物结构等人为因素两方面得出河北平原地下水开采是引起地下水水位下降的最主要因素，相对贡献率为54.7%，其次是降水和地表径流，相对贡献率为25.6%和19.7%。张展羽等[23]基于主成分分析与多变量时间序列(controlled auto-regressive，CAR)构建了济南市陡沟灌区地下水水位模型，通过对比主成分回归模型、逐步回归模型和BP神经网络模型的模拟结果，发现主成分-时间序列模型模拟结果与实测值相对误差最小，模拟效果较好。

在已有研究中，Wang等[18]利用同位素示踪法发现华北山前平原的元氏地区上游地下水侧向补给元氏的地下水量占地下水总补给量的50%～88%，虽然上游的侧向补给对华北山前平原地下水影响较大，但目前考虑上游地区对山前平原地下水动态变化的定量影响研究尚未见报道，且华北山前平原浅层地下水埋深较大，地质条件复杂，目前对该地区地下水水位变化原因的定量分析较少。本文选取华北山前平原典型井灌区石家庄栾城为研究区，从气候条件和人类活动两方面对栾城地下水位变化的原因进行分析，以期为合理制定“以水定产”的地下水管理措施和促进高效农业节水灌溉提供参考。

1 研究数据和方法

1.1 研究区概况

华北山前平原位于华北平原的西部和北部的狭长地带，宽度一般在30～60 km，土层深厚，有机质含量较高，地下水自西北流向东南。受季风和地形等因素影响，山前平原年降水波动较大，年内降水时间分布不均，6—9月降水量约占全年的70%～80%[13]。气候变化背景下极端气候事件增多以及人类活动强度日趋增大等因素使该地区地下水水位变化原因复杂，且各影响因素间存在一定的相互作用。

石家庄栾城区为华北山前平原石家庄市东南部的典型井灌区，该地区地势平坦，以冬小麦和夏玉米轮作制农田为主，农田有效灌溉面积比例达到100%，占区域总面积的68%。栾城区冬小麦生育期自10月初至次年6月初，夏玉米生育期自冬小麦收获至9月末，农业用水占全部用水量的88%以上。

研究区的地下水位观测点位于中国科学院栾城农业生态系统试验站内(37.53°N，114.41°E，海拔为50.1 m)，该试验站自1974年开始记录地下水水位数据，其地下水水位波动趋势与华北山前平原区地下水水位变化具有较好的一致性[11,24]，能较好地代表栾城冬小麦和夏玉米轮作系统下的地下水水位变化。

1.2 研究数据

1975—2015年地下水位埋深观测数据来自中国科学院栾城农业生态系统试验站，该数据采用加拿大Solinst公司制造的地下水水位尺(-101型)设备观测，仪器精度为±1 mm。农田实际蒸散月数据(2008—2015年)由涡度相关系统观测并计算得到[25-26]。气象数据(1994—2015年)来自河北省气象局气象信息中心，包括日降水量、日平均气温、日最高气温、日最低气温、日平均风速、日平均相对湿度和日照数据等。社会经济数据(1994—2015年)来自《河北农村统计年鉴》，包括石家庄市区人口数和栾城、正定、平山、井陉4县(区)的小麦和玉米的种植面积、有效灌溉面积和农村用电量。

1.3 研究方法

1.3.1 水分盈亏

本文将年降水量减去蒸散量定义为水分盈亏量：

D=P-ET

(1)

式中：D为水分盈亏量，mm；P为年降水量，mm；ET为蒸散量，mm。当水分盈亏量为正值时，表明该地区水分有盈余；当水分盈亏量为负值时，表明年降水量不足以弥补作物生长的蒸散量，此时出现水量亏缺。当蒸散量为涡度相关系统观测值时，ET为实际蒸散量，本文利用实际蒸散量等数据分析研究区的年内水分盈亏变化。

由于研究区实际蒸散量数据的时间序列较短，在分析1994—2015年地下水水位变化影响因素中的水分盈亏量时，选用1998年联合国粮农组织(FAO)推荐并修订的Penman-Monteith模型计算参考作物的潜在蒸散量：

(2)

式中：ET0为参考作物的潜在蒸散量，mm/d；Δ为饱和蒸汽压与空气温度曲线的斜率，kPa/℃；Rn为参考作物净辐射，MJ/(m2·d)；G为土壤热通量，MJ/(m2·d)；γ为干湿表常数；T为2 m高处日平均气温，℃；u2为2 m高处日平均风速，m/s；es为饱和水汽压，kPa；ea为实际水汽压，kPa。

1.3.2 数据标准化

由于不同类型数据的单位和数量级不同，在进行主成分分析之前，需要先将原始数据标准化。z-score是常用的标准化方法之一，该方法基于原始数据的均值和标准差对数据进行标准化处理，标准化后为无量纲数据，均值为0，方差为1，计算公式为

(3)

式中：x为原始数据；μ为全部数据的均值；σ为标准差。

1.3.3 主成分回归分析法

华北山前平原地下水水位变化原因复杂，在进行变量分析时，为了消除多个影响因子间的相互作用，本文采用主成分回归分析法进行拟合。主成分回归分析是多元统计分析中减少自变量维数的一种较理想的工具[27]，因此被广泛应用于气候变化和水资源等领域。

主成分回归分析的基本思路是将原来的指标体系重新组合成一组(若干个)新的且相互无关的综合指标，从中选取自变量所有线性组合中方差最大者作为主成分，并尽可能多地反映原来指标体系的信息，最后以主成分代替原变量进行回归分析[28-29]。在主成分基础上构建了多元线性回归模型，将回归方程中各主成分系数的权重作为各主要影响因子对地下水水位变化的贡献率。

图1 1975—2015年观测点地下水水位年内和年际变化Fig.1 Annual and interannual change of groundwater level in the observation point during 1975-2015

图2 2008—2015年观测点水平衡项和地下水水位月变化Fig.2 Monthly changes of water balance and groundwater level in observation point during 2008-2015

1.3.4 多变量时间序列分析法

多变量时间序列分析法[30]可从多变量的时间序列中提取有效信息表征复杂系统的动态特征，广泛应用于水文预测，具有较好的预测效果。在主成分回归分析基础上，将主成分回归分析构造的综合影响因子作为输入变量，导入多变量时间序列CAR模型[31]，采用递推最小二乘法对估计模型参数由低阶向高阶进行多变量时间序列模拟，通过拟合优度确定模型的最高阶次。为了检验所构建的多变量时间序列模型，本文从1994—2015年共21个样本中选取前16个构建模型，用剩余样本进行检验。

对m个变量的时间序列组构建n阶多变量时间序列模型，计算公式为

yt=a1yt-1+a2yt-2+…+anyt-n+b10x1,t+b11x1,t-1+…+b1nx1,t-n+b20x2,t+

b21x2,t-1+…+b2nx2,t-n+…+bm0xm,t+bm1xm,t-1+…+bmnxm,t-n+εt

(4)

式中：aj、bij为多变量时间序列模型系数(i=1，2，…，m；j=1，2，…，n)；yt-j、xi,t-j为因变量和自变量的时间序列；t为时间序列，t>1。

1.3.5 模型检验

模型验证过程中采用决定系数(R2)、模拟值与实测值的均方根误差(RMSE)、相对误差(MRE)和纳什系数4个评价指标来表征模型的模拟效果。

2 结果与分析

2.1 地下水水位年变化

图1为1975—2015年观测点地下水水位年内和年际变化，1975—2015年栾城年平均地下水水位呈显著下降趋势(P<0.001)，平均每年下降0.84 m，年平均地下水水位变化幅度在-2.8～2.9 m之间。1975—2015年栾城年平均地下水水位由37.4 m下降至5.9 m，同时，地下水埋深不断增加，由1975年的12.7 m增加至2015年的44.2 m。期间仅有5年(1976年、1977年、1982年、1996年和2003年)年平均地下水水位回升，其中1996年水位上升幅度最大，回升近3 m。1997年和1998年连续两年降水量偏少使地下水水位下降幅度均超过2.5 m，2004—2015年地下水水位出现连续11年持续下降现象，地下水资源持续减少。

2.2 影响地下水水位变化的气象因素

受季风气候影响，2008—2015年栾城年内降水量集中在6—9月，占全年降水量的77.9%，其他时段月降水量大多不足30 mm。由图2可以看出，栾城典型井灌区实际蒸散观测值存在两个峰值，与一年两熟的作物需水变化相对应，分别是冬小麦需水量最大的5月和夏玉米需水量最大的8月。夏玉米生育期与当地气候雨热同季，且降水量高于实际蒸散量，因此在8—9月地下水水位有所回升。冬小麦生育期(3—5月)降水量较少，不能满足作物生长需要，地下水高强度灌溉导致地下水水位在该期间持续下降，多年平均地下水水位累积下降2.9 m左右。由此可见，栾城降水量和水分盈亏是影响其地下水动态变化的主要气象因素，2008—2015年栾城区高强度地下水灌溉是保障该地区农作物稳产、高产的关键[15]。

2.3 影响地下水水位变化的人为因素

2.3.1 当地农业生产

华北平原地下水水位年际变化主要受人为因素影响，农业灌溉是引起当地地下水水位下降的主要影响因素。1967年以来，华北平原冬小麦和夏玉米种植面积明显增加，地下水水位开始逐年下降，种植结构的变化对栾城地下水动态变化影响显著，地下水采补失衡严重[32]。20世纪50年代至21世纪初，石家庄平原农业地下水开采量与小麦和玉米总种植面积存在密切相关性，二者变化趋势较为一致，年际和年代际地下水农业开采量随小麦和玉米总种植面积的变化而变化[11]。因此，本文将栾城小麦和玉米的种植面积作为当地人类活动对地下水动态变化的影响因子。

2.3.2 上游人为因素

栾城农田井灌区在雨季(6—9月)后至次年3月冬小麦春灌前，地下水水位以持续回升为主，除雨季降水直接补给外，上游地区地下水侧向补给量占总补给量的50%以上[18]。20世纪70年代末以来，栾城区境内交河、沙河和泥河上游先后修建中小型水库30余座，使该地区4条季节性河流的年均径流量明显减少，沙河和泥河几乎无水下泄，河流补给地下水量也显著减少[11]。根据《华北平原地下水可持续利用地图集》[33]中华北平原浅层地下水等水位线图，栾城及周边地下水水位自西北向东南方向递减，假设上游距离栾城越近，对栾城地下水补给影响越大。通过对比栾城区西北方向的石家庄市区和邻近3个县(正定、平山和井陉)的小麦和玉米种植面积、有效灌溉面积、农村用电量和人口数量等因子在不同组合下构建的主成分回归模型的显著性差异，发现石家庄市区人口作为上游人为因素时模拟效果最佳，这主要是由于石家庄市区为栾城区地下水上游最邻近地区，近些年石家庄市区城市化进程加快，市区水资源消耗和地下水开采量显著增加，影响了对下游栾城区的地下水侧向补给，因此，本文将石家庄市区人口作为栾城上游的主要人为因素。

2.4 气候和人为因素对地下水动态变化影响的定量分析

表1 1994—2015年潜在自变量的统计特征量

本文在分析栾城及其上游气候因子和人类活动对栾城地下水水位动态变化的定性影响基础上，考虑了栾城及其上游(正定、平山、井陉和石家庄市区)的降水量、水分盈亏量、小麦和玉米种植面积、农村用电量、人口数量等多个影响因子，通过对比不同组合下的拟合模型的模拟结果，最终将栾城地下水水位作为因变量，选取栾城降水量(x1)、栾城水分盈亏量(x2)、栾城小麦和玉米种植面积(x3)、栾城上游降水量(x4)和石家庄市区人口(x5)作为潜在自变量，潜在自变量在1994—2015年的统计特征量见表1。

表2 各主成分的特征值、贡献率和累积贡献率

从表2可以看出，最大的特征值为1.670，其贡献率为55.8%，前3个主成分的累积贡献率高达97.0%，可以代表原指标97.0%的信息，用这3个主成分即可代替原指标，而且他们是相互独立的影响因子。用表3中前3个主成分的特征向量为系数构造各主成分的线性方程：

Z1=-0.591x1-0.584x2-0.551x4

(5)

Z2=-0.781x3-0.194x4-0.585x5

(6)

Z3=-0.557x3-0.180x4+0.810x5

(7)

从主成分的线性方程来看，式(5)中栾城降水量(x1)的

表3 各主成分特征值对应的特征向量

系数最大，对Z1的波动起主导作用，所以将第一主成分Z1视作本地降水因子；式(6)中栾城小麦和玉米的种植面积(x3)的系数最大，对Z2的变化起主导作用，将第二主成分Z2视作本地人为因子；式(7)中上游石家庄市区人口(x5)的系数最大，对Z3的变化起主导作用，可视作上游主要影响因子。这3个因子彼此之间相互独立，不存在共线性，且代表了原指标的大部分信息。以本地降水因子Z1、本地人为因子Z2和上游影响因子Z3为自变量，以地下水水位为因变量，利用最小二乘法进行多元线性拟合后，得到的回归方程为

Y=0.112 6Z1-0.490 5Z2-0.360 6Z3

(8)

回归方程的相关系数R为0.62，决定系数R2为0.39，F检验值为3.79，通过了0.01的显著性检验，说明回归方程拟合效果较理想。从回归方程中3个因子系数的权重得出本地降水因子对栾城地下水水位回升的贡献率为11.7%，本地人为因子贡献率为-50.9%，上游影 响因子贡献率为-37.4%(正负号表示对地下水水位回升的正负作用)。由此可以看出，栾城小麦和玉米种植面积对地下水水位变化的影响最大，人为因子在地下水水位的变化中起主导作用。与人为因子相比，降水因子对地下水水位影响较小，地下水水位随降水量的增加出现不同程度的回升。总体上，在以栾城典型井灌区为代表的太行山山前平原区，其地下水水位变化主要来自人类活动的影响，受自然因素影响的比重仅占11.7%。

2.5 主成分-时间序列模型

地下水水位建模对流域水资源管理具有重要意义和必要性，也是井灌区用水管理的重要依据[34]。束龙仓[35]在分析石家庄地下水水位影响因素和构建回归模型时发现本月地下水水位与本月降水量、本月开采量和前一个月水位关系较密切。本文借鉴随机理论和成因理论的思想，在构建栾城地下水水位模型时，基于上文主成分回归分析所选出的3个主要综合因子，考虑上一年地下水水位对本年度的影响，根据决定系数将模型定阶到上一年(t-1)，剔除不显著因子后，将上一年地下水水位、当年的Z1、上一年的Z2和上一年的Z1作为自变量，以当年地下水水位为因变量，利用1994—2010年数据构建研究区地下水水位主成分-时间序列模型：

Y=-0.185 13+1.178 46Yt-1-0.024Z1，t+0.085 6Z2，t-1-0.081 5Z1，t-1

(9)

该主成分-时间序列模型的模拟值与实测值的相关系数R为0.994，决定系数R2为0.988，平均相对误差为0.08，均方根误差为0.132，纳什系数为0.98。

图3 2011—2015年观测点地下水水位模拟值和观测值Fig.3 Comparison between simulated and measured groundwater level in test station of Luancheng during 2011—2015

为检验上述模型的效果，对比2011—2015年栾城典型井灌区模型的模拟值与地下水水位观测值(图3)发现，验证期的模拟值与观测值变化趋势一致，相对误差在-19.4%～2.1%之间，说明该模型适用于栾城井灌区的地下水水位预测，即利用降水量、水分盈亏量、种植规模等数据可以较好地推测当年的地下水水位。

2.6 讨论

华北平原人口众多，工农业用水量大，地下水是该地区最重要的供水水源和战略资源。在气候变化和人类活动共同影响下，华北山前平原土地利用、种植结构和城镇化的发展和变化显著增加了地下水开采量，地下水水位持续下降。虽然栾城在1994—2015年小麦和玉米种植面积有所减少，但粮食产量和单产总体增加使作物实际耗水量增加，作物生育期的降水量不能满足作物需水量，即使在雨季地下水水位仍不能得到足够的补给，地下水开采量长期大于补给量[32]。由于数据来源、分析方法和对地下水补给排泄认知的差异，加上栾城地下水水位变化原因十分复杂，已有研究中对地下水持续下降的原因分析结论不一。本文针对研究区地下水补给特点，综合考虑了栾城本地及其上游的气候因素和人类活动对栾城地下水水位变化的影响，利用主成分回归分析法排除各影响因子间的相互作用，得出综合影响因子及其贡献率，可为合理开发地下水资源提供科学依据。

由于目前地下水观测数据不足，对地下水补给量和补给来源的认识尚存在一定分歧，限制了我们对地下水和气候之间动态关系的理解[36-37]，加之观测数据样本量较小，时间尺度较短，缺少准确的地下水补径排和开采量数据等原因，给定量分析研究区地下水水位变化原因带来一定的不确定性。因此，地下水水位模型构建还有待进一步考虑其他参数，如地下水水位的不同埋深范围、地下水水位的年变化、前期土壤含水量、土壤性质(尤其是表层土壤性质)、季节以及降水强度等[38]。此外，栾城上游降水量减少和植被覆盖度增加，使自然降水补给减少的同时增加了降水量的截留，进一步减少了对下游地下水的侧向补给量；近些年来，华北山前平原种植类型多样化，土地利用变化对地下水的影响也有待于进一步分析。

3 结 论

a.1975年以来，栾城降水量增加对地下水补给具有正效应，当地人类活动对地下水水位下降影响最大，同时，栾城上游人类活动对栾城地下水水位下降的贡献也不容忽视。1994—2015年，华北山前平原地下水水位变化主要受本地气候因子、本地人类活动和上游人类活动三方面影响因子共同作用，各因子对当地地下水水位回升的贡献率分别为11.7%、-50.9%和-37.4%，这三方面影响因子的主导因素分别为栾城降水量、栾城小麦和玉米种植面积和上游石家庄市区人口。

b.构建地下水水位的主成分-时间序列模型，经检验模型的模拟值与实测值趋势一致，表明该模型能较好地预测华北山前平原典型井灌区多因子影响下的地下水水位变化。

c.在开发利用华北平原地下水资源时，应充分考虑上游对下游的可能影响，尽量减少对下游水资源补给的不利影响。通过构建统计关系模型，可以在一定程度上解决目前对复杂机理问题认知的不确定性。基于主成分-时间序列模型可以根据上一年地下水水位和降水量等数据，合理规划当地农业种植结构和种植规模，统筹考虑预期粮食产量与地下水开采量，从而有效缓解华北平原地下水长期超采的严峻局面。