衡水地区“冀枣衡”漏斗演化特征及主控因子研究

朱 琴，田晓华，杨志光

1.河北省水文工程地质勘查院(河北省遥感中心)，河北 石家庄 050021；2.河北省地下水资源与生态环境安全技术创新中心，河北 石家庄 050021

0 前言

地下水是自然资源的重要组成部分，其水位变化受到多种复杂因素的影响。这些影响因素既包括人类活动(如地下水开采)，又包括固有的水文地质条件。然而，由于水文地质条件与地下水位变化关系的复杂性，以及不同地区地质条件的巨大差异，量化地质要素对地下水水位变化的影响成为一个难题。尽管现有的地下水数值模拟能够在一定程度上将水文地质要素数据代入模型中，虽然这些模型可以对地下水水位进行预测，但是仍然难以准确量化水文地质条件对地下水水位变化的实际影响。

在已有的研究中，关于漏斗或者地下水水位影响因子的研究分为3种类型，第1种类型为仅做定性分析，认为地下水开采是水位变化的主控因子。王淑虹等[1]定性分析了玛纳斯河流域绿洲区地下水水位变化的影响因素，梁中华[2]定性分析了大连主要灌区地下水埋深变化影响因子，陶建华[3]定性分析了宿州市地下水漏斗影响因子，马波等[4]定性分析了平罗县西大滩区潜水降落漏斗影响因子，古力米日·买买提[5]定性分析了叶尔羌河主要灌区地下水埋深变化影响因子，大家普遍认为主要影响因子为地下水开采。第2种为采用数学分析方法，大部分采用回归分析，努尔麦麦提·艾麦提[6]采用多元线性回归分析计算了贡献率，计算结果表明气温变化对区域地下水埋深变化影响较大，相对贡献率高于50%，地下水开采影响相对贡献率高于30%。刘阳[7]采用多元线性回归，计算结果表明地下水开采量是阜新地区地下水埋深变化的主因，各区县影响相对贡献率均超过30%，其次为降水量，气温和蒸发影响贡献率相对较低。张婧等[8]采用多元线性回归，计算结果表明气候因子、本地人类活动和上游人类活动对当地地下水水位回升的贡献率分别为11.7%、-50.9%和-37.4%。刘静等[9]采用多元线性回归分析计算了贡献率，GDP、工业用水是辽河平原地下水埋深变化的关键影响因子。陈彬鑫等[10]采用多元线性回归计算表明地下水开采量已逐渐成为影响莫索湾灌区地下水埋深的主要因素。文静等[11]采用多元线性回归，计算表明人口规模、上游来水量、地下水开采量是影响河西走廊地下水埋深变化的主要因素。徐月琴等[12]采用多元线性回归，计算了降水和蒸发量对潍坊市潍北区地下水降落漏斗的影响。白宜斐等[13]采用多元线性回归，计算表明叶尔羌河灌区地下水埋深加深的主导因子是地下水的超采。其他方法，如史彩霞等[14]采用相关性分析，发现地下水过量开采与区域性干旱气候是太原盆地内地下水降落漏斗形成的主要原因。于晓露等[15]采用Sobol全局敏感性分析方法识别影响内蒙古黄旗海盆地地下水资源演变的主要控制因素。翟小艳[16]运用灰色关联度法分析表明城区—宋古漏斗区主要影响因素为地下水过量开采。杨阿敏等[17]采用灰色关联度法揭示西部地下水位的主导影响因子为降雨量，东部地下水位的主导影响因子为渗漏量和降雨量。石媛媛等[18]基于主成分分析法，计算了民勤荒漠绿洲干旱累计方差贡献率，气象干旱的累计方差贡献率为86.853%；水文干旱受制于民调水量和出库水量，累计方差贡献率为92.190%；农业干旱主要是由耕地面积与有效灌溉面积增大所致，累计方差贡献率为86.179%。张曼菲等[19]采用主成分分析法和灰色关联分析法确定影响泾惠渠灌区地下水埋深动态的主要因素。结果表明：灌区1981—2010年地下水位平均埋深基本呈增加趋势，突变点出现在1994年，1994年之前影响灌区地下水埋深动态的主要因素为蒸发量，1994年之后主要影响因素为蒸发量、地下水开采量和灌渠引水量。王默涵等[20]采用地理探测器定量分析了邯郸平原区的主要影响因子为工业用水变化和降水量变化。第3种方法是采用当前较为热门的机器学习模型进行分析。关于“冀枣衡”漏斗影响因子的研究以定性研究为主，尹新明[21]认为超采是“冀枣衡” 漏斗形成的主要原因。王建恒[22]分析了气温、日照、蒸发量等气象要素对降水及衡水地区地下水水位变化的影响。

鉴于目前关于“冀枣衡”漏斗影响因子的研究以定性分析为主，定量分析不同时期主控因子的贡献率的研究和研究地质因子对漏斗贡献的相对较少，为此本次采用多元线性回归这种较为成熟的研究方法，将地下水开采和水文地质因子作为影响因子，进行地下水水位降落漏斗的影响因子分析，有助于我们更全面地了解“冀枣衡”地区地下水水位的动态变化，以期为未来地下水资源的合理开发和保护提供科学依据。

1 “冀枣衡”漏斗演化特征

“冀枣衡”漏斗是深层地下水开采形成的漏斗，位于冀州—枣强—衡水一带，地下水开采深度150～160 m，开采层位为第三含水层组，形成于1969年。1969年衡水市附近深层地下水位埋深约2.9 m，“冀枣衡”地区地下水位埋深约9.3 m，漏斗中心位于衡水市衡62观测孔，漏斗面积只有几平方千米，受工业和农业集中开采影响，漏斗面积不断加深和扩大，到1973年5月底，水位埋深大于10 m，范围已达2 980 km2。近十几年超量开采深层水，致使水位逐年下降。2004年，漏斗北部已扩展至衡水地区外，东南部已延伸至邢台、沧州、德州地区，形成大型的复合型深层水位降落漏斗，漏斗面积为8 815 km2，覆盖整个衡水地区。2016年至今，受景县及邻近山东省德州市地下水集中开采的影响，漏斗中心位于景县留智庙镇八里庄村，中心水位埋深在120.86～129.27 m之间；2020年漏斗中心移向西南故城县房庄乡堤口村，中心水位埋深为130.31 m。

表1是“冀枣衡”漏斗的统计数据。图1是“冀枣衡”漏斗中心水位变化图，从图1可以看出1975—2020年漏斗中心水位持续下降，水位埋深范围为32.44～130.31 m，下降幅度为97.87 m，水位下降速率为2.11 m/a。漏斗的影响面积自1975年持续扩大，至2004年影响面积扩展到整个衡水市，2004年以后影响面积的统计数据为8 815 km2(衡水市面积)。图2为漏斗封闭面积历史变化图，从图中可以看出漏斗封闭面积总体变化趋势为波动中持续扩大。近几十年“冀枣衡”漏斗的位置不断的变化，由西南向东北方向运移。漏斗封闭区域历史范围见图3。从图3可以看出，1975年漏斗位于衡水市桃城区；1984年漏斗位于南宫市和枣强县；1992年漏斗范围较大，包括武邑县、景县、桃城区、枣强县和南宫市；2005年漏斗范围最大，从南到北，包括了邢台南宫市、衡水市和沧州献县、东光县等地；2016年漏斗范围较小，西侧桃城区、武邑县等地不再位于封闭漏斗范围内，漏斗主要包含枣强县、固城县和景县；2020年漏斗向东北发展，包含武邑县、阜城县、景县、吴桥县和固城县。

图1 “冀枣衡”漏斗中心埋深变化图Fig.1 Variation of burial depth in the center of the “JZH” funnel

图2 “冀枣衡”漏斗封闭面积变化图Fig.2 Variation of closed area of “JZH” funnel

图3 “冀枣衡”漏斗封闭范围历史变化图Fig.3 Historical variation of the closure range of the “JZH” funnel

表1 “冀枣衡”漏斗统计表Table 1 Statistical table of “JZH” funnel

“冀枣衡”地区4期的地下水水位变幅数据(1985—1992年、1992—2005年、2005—2016年、2016—2020年)如图4所示。从图中可以看出，1984—1992年水位下降较大的地区位于枣强县和桃城区，水位下降幅度为30 m左右，深州市的西北侧水位有一定的回升，上升幅度为20 m左右；1992—2005年水位下降较大的地区位于景县的北边和南边，水位变幅为44 m左右，整个衡水地区地下水水位变幅均为负值，说明区域水位均下降；2005—2016年水位下降较大的地区位于景县东边，水位变幅为26 m左右，深州市和武强县水位出现了一定的回升，其中深州市的回升幅度为10 m左右，武强县为4 m左右；2016—2020年水位下降较大的地区位于冀州区和衡水市的西北角，水位变幅为35 m左右，其次为故城县和阜城县，水位变幅为25 m左右，与阜城县相邻的沧州泊头市水位出现了一定的回升，回升幅度为20 m左右。

总体来说，1984—2020年“冀枣衡”地区的地下水水位总体呈下降趋势。1992—2005年是水位下降幅度最大的时间段，区域水位整体下降，而其他3个阶段均存在局部水位上升。这4个阶段水位变幅由大到小依次为：1992—2005年、2016—2020年、1985—1992年、2005—2016年。这4个阶段水位变化速率由大到小依次为：2016—2020年、1985—1992年、1992—2005年、2005—2016年。

2 漏斗发展演化的影响因子分析

已有的研究表明衡水地区深层地下水的来源主要有侧向径流补给、越流补给和黏性土压密释水[23]。其中浅层水对深层水的越流补给主要发生在第一含水组和第二含水组之间，第三、第四含水组主要来源于储存量和侧向径流补给量[24]。深层地下水的主要排泄项为地下水开采。深层地下水开采是影响漏斗发展的主要因素，通过对比开采量与水位埋深(图5)，可知开采与埋深变化趋势总体一致。两者之间进行相关性分析的结果显示(图6)，两者之间具有良好的相关性，相关系数R2为0.804 5。

图6 开采量与“冀枣衡”漏斗中心埋深相关性分析图Fig.6 Correlation analysis between mining volume and burial depth of “JZH” funnel center

一般来说，降雨量比较大的年份开采量相对较小，从而导致地下水水位下降速率较慢，由此说明降雨量对地下水水位是间接影响。对比降雨量与开采量(图7)和降雨量与漏斗中心水位(图8)。从图7可以看出，降雨量与开采量在1990年前有一定的对应关系，1990年前随着降雨量的增加，开采量有所减少。但是从图8 可以看出 1990 年后，水位仍旧持续下降，降雨量相对稳定，开采量持续升高。由此说明，降雨量在1990年前对开采量有所影响，但是影响强度较小，水位持续下降，1990年后降雨量对开采量影响较弱，水位持续下降。利用1990年前后降雨量与开采量数据进行相关性分析(图9)，1990年前两者之间具有良好的负相关性，相关系数R2为0.76。1990年后两者之间相关性较差。由此说明1990年后随着开采量的变化增加，受降雨量影响较小。

图7 “冀枣衡”漏斗降雨量与开采量变化图Fig.7 Variation of rainfall and mining volume in “JZH” funnel

图8 降雨量与漏斗中心水位变化图Fig.8 Rainfall versus water level in the center of the funnel

图9 开采量与降雨量相关性分析图Fig.9 Analysis of correlation between mining and rainfall

3 漏斗空间水位变化的影响因素分析

本次以“冀枣衡”漏斗为研究对象。利用ArcGIS软件在研究区范围内随机生成620个点，提取这些点的弹性释水系数数据、渗透系数数据、底板埋深、含水层厚度数据、水位变幅数据(1985—1992年、1992—2005年、2005—2016年、2016—2020年)。采用SPSS软件进行回归分析，分析4个时间阶段水文地质要素对水位变幅的影响强度。结合开采量数据分析回归分析结果的合理性。其中弹性释水系数数据利用河北平原的1 051个钻孔数据计算每个钻孔弹性释水系数，然后利用ArcGIS对点数据进行克里金插值，利用插值后栅格数据提取本次研究点的弹性释水系数数据；渗透系数的数据利用钻孔的抽水试验获得的渗透系数进行克里金插值后获得；顶板埋深采用第一含水层底板埋深栅格数据获得；含水层厚度数据采用第一含水层底板埋深与第四含水层底板埋深进行算术计算获得；水位变幅数据利用1982年、1992年、2005年、2016年和2020年水位栅格数据进行计算获得。

3.1 “冀枣衡”地区各特征因子分布特征

“冀枣衡”地区的渗透系数和弹性释水系数普遍偏低(图10、图11)，渗透系数数值范围为3.99～14.4 m/d，冀州和景县东北地区相对较高，其他地区偏低。弹性释水系数的数值范围为0.000 1～0.002 5，吴桥和深州等地较高，其他地区偏低。“冀枣衡”地区深层水的底板标高-450 m左右，顶板标高为20 m左右，深层水厚度400多米。

图10 研究区渗透系数分布图Fig.10 Spatial distribution of infiltration coefficient in the study area

图11 研究区弹性释水系数分布图Fig.11 Spatial distribution of elastic water release coefficients in the study area

图12 地质要素对“冀枣衡”漏斗影响的解释强度Fig.12 Interpreted intensity of the influence of geological elements on “JZH” funnel

3.2 水文地质指标对水位变化的贡献率分析

通过对研究区620个点的含水层厚度数据、弹性释水系数数据、渗透系数数据、含水层顶板埋深、1984—1992年水位变幅数据进行分析，采用spss软件进行回归分析，结果见表2、表3、表4。4期分析数据的Durbin—Watson系数均接近2，满足样本独立要求。共线性诊断的显著性P值为0.000a，小于0.05，说明自变量之间相互独立，满足计算模型要求。共线性统计量VIF均小于5，说明模型计算结果可靠。

表2 模型汇总b表Table 2 Summaryb of the mldel

表3 Anovab表Table 3 Table of Anovab

表4 系数a表Table 4 Table of coefficientsa

“冀枣衡”漏斗4个阶段的分析可知，1984—1992年、1992—2005年、2005—2016年、2016—2020年水文地质指标对水位变幅的解释度R2依次为：0.382、0.212、0.045、0.163(表2)。1984—1992年对水位变幅影响较大的水文地质指标依次为顶板埋深、弹性释水系数；1992—2005年对水位变幅影响较大的水文地质指标依次为弹性释水系数、顶板埋深、含水层厚度；2005—2016年对水位变幅影响较大的水文地质指标依次为含水层厚度、顶板埋深、渗透系数；2016—2020年对水位变幅影响较大的水文地质指标依次为含水层厚度、顶板埋深、渗透系数、弹性释水系数。

一般认为漏斗发展的主控因子是地下水开采，随着地下水开采量的增加，开采对地下水水位下降的影响强度逐渐增加，从而导致地质要素对漏斗的影响强度逐渐降低。结合图7可知，1984—1990 年开采量逐渐减少；1990—2005 年是开采量快速增长；2005—2016 年开采量在波动中继续增长，2016—2020 年开采量逐渐下降。 与地质要素对水位变化解释率的变化具有良好负相关。由此说明，回归分析结果符合一般认识。

4 结论

近45年来“冀枣衡”漏斗中心水位持续下降，多年平均下降速率为2.1 m/a，漏斗面积持续扩大，自2014年漏斗面积已近扩大至衡水全境。漏斗演化特征为由西北向东南方向发展，漏斗中心由早期的衡水市枣强县，向东南部的景县、固城县发展。4个时期的水位下降幅度显示，1985—1992年水位下降幅度较大地区为桃城区和枣强县，水位降幅约25 m；1992—2005年为景县和固城县，降幅约44 m；2005年后水位下降速率减小，2005—2016年为景县，降幅约26 m；2016—2020年为固城县，降幅约25 m。

地下水开采是影响漏斗的主控因子，漏斗中心水位埋深与开采量具有良好的相关性，1990年前降雨量与开采量具有良好的相关性，1990年后降雨量与开采量之间的相关性较小。应用回归分析和GIS技术分析漏斗区水位变化影响因子，结果显示水文地质因子对研究区地下水水位空间变化的影响解释率为：1984—1992年为38.2%、1992—2005年为21.2%、2005—2016年为4.5%、2016—2020年为16.3%，该结果与研究区开采强度具有良好的对应关系，随着开采强度的增加地质要素对漏斗影响强度在降低，符合一般认识。研究成果可以为衡水地区的地下水修复提供技术参考。