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

朱 琴，田曉華，楊志光

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%，該結果與研究區開采強度具有良好的對應關系，隨著開采強度的增加地質要素對漏斗影響強度在降低，符合一般認識。研究成果可以為衡水地區的地下水修復提供技術參考。