朱 琴，田曉華，楊志光

1.河北省水文工程地質(zhì)勘查院(河北省遙感中心)，河北 石家莊 050021；2.河北省地下水資源與生態(tài)環(huán)境安全技術(shù)創(chuàng)新中心，河北 石家莊 050021

0 前言

地下水是自然資源的重要組成部分，其水位變化受到多種復(fù)雜因素的影響。這些影響因素既包括人類活動(如地下水開采)，又包括固有的水文地質(zhì)條件。然而，由于水文地質(zhì)條件與地下水位變化關(guān)系的復(fù)雜性，以及不同地區(qū)地質(zhì)條件的巨大差異，量化地質(zhì)要素對地下水水位變化的影響成為一個難題。盡管現(xiàn)有的地下水?dāng)?shù)值模擬能夠在一定程度上將水文地質(zhì)要素數(shù)據(jù)代入模型中，雖然這些模型可以對地下水水位進(jìn)行預(yù)測，但是仍然難以準(zhǔn)確量化水文地質(zhì)條件對地下水水位變化的實(shí)際影響。

在已有的研究中，關(guān)于漏斗或者地下水水位影響因子的研究分為3種類型，第1種類型為僅做定性分析，認(rèn)為地下水開采是水位變化的主控因子。王淑虹等[1]定性分析了瑪納斯河流域綠洲區(qū)地下水水位變化的影響因素，梁中華[2]定性分析了大連主要灌區(qū)地下水埋深變化影響因子，陶建華[3]定性分析了宿州市地下水漏斗影響因子，馬波等[4]定性分析了平羅縣西大灘區(qū)潛水降落漏斗影響因子，古力米日·買買提[5]定性分析了葉爾羌河主要灌區(qū)地下水埋深變化影響因子，大家普遍認(rèn)為主要影響因子為地下水開采。第2種為采用數(shù)學(xué)分析方法，大部分采用回歸分析，努爾麥麥提·艾麥提[6]采用多元線性回歸分析計算了貢獻(xiàn)率，計算結(jié)果表明氣溫變化對區(qū)域地下水埋深變化影響較大，相對貢獻(xiàn)率高于50%，地下水開采影響相對貢獻(xiàn)率高于30%。劉陽[7]采用多元線性回歸，計算結(jié)果表明地下水開采量是阜新地區(qū)地下水埋深變化的主因，各區(qū)縣影響相對貢獻(xiàn)率均超過30%，其次為降水量，氣溫和蒸發(fā)影響貢獻(xiàn)率相對較低。張婧等[8]采用多元線性回歸，計算結(jié)果表明氣候因子、本地人類活動和上游人類活動對當(dāng)?shù)氐叵滤换厣呢暙I(xiàn)率分別為11.7%、-50.9%和-37.4%。劉靜等[9]采用多元線性回歸分析計算了貢獻(xiàn)率，GDP、工業(yè)用水是遼河平原地下水埋深變化的關(guān)鍵影響因子。陳彬鑫等[10]采用多元線性回歸計算表明地下水開采量已逐漸成為影響莫索灣灌區(qū)地下水埋深的主要因素。文靜等[11]采用多元線性回歸，計算表明人口規(guī)模、上游來水量、地下水開采量是影響河西走廊地下水埋深變化的主要因素。徐月琴等[12]采用多元線性回歸，計算了降水和蒸發(fā)量對濰坊市濰北區(qū)地下水降落漏斗的影響。白宜斐等[13]采用多元線性回歸，計算表明葉爾羌河灌區(qū)地下水埋深加深的主導(dǎo)因子是地下水的超采。其他方法，如史彩霞等[14]采用相關(guān)性分析，發(fā)現(xiàn)地下水過量開采與區(qū)域性干旱氣候是太原盆地內(nèi)地下水降落漏斗形成的主要原因。于曉露等[15]采用Sobol全局敏感性分析方法識別影響內(nèi)蒙古黃旗海盆地地下水資源演變的主要控制因素。翟小艷[16]運(yùn)用灰色關(guān)聯(lián)度法分析表明城區(qū)—宋古漏斗區(qū)主要影響因素為地下水過量開采。楊阿敏等[17]采用灰色關(guān)聯(lián)度法揭示西部地下水位的主導(dǎo)影響因子為降雨量，東部地下水位的主導(dǎo)影響因子為滲漏量和降雨量。石媛媛等[18]基于主成分分析法，計算了民勤荒漠綠洲干旱累計方差貢獻(xiàn)率，氣象干旱的累計方差貢獻(xiàn)率為86.853%；水文干旱受制于民調(diào)水量和出庫水量，累計方差貢獻(xiàn)率為92.190%；農(nóng)業(yè)干旱主要是由耕地面積與有效灌溉面積增大所致，累計方差貢獻(xiàn)率為86.179%。張曼菲等[19]采用主成分分析法和灰色關(guān)聯(lián)分析法確定影響涇惠渠灌區(qū)地下水埋深動態(tài)的主要因素。結(jié)果表明：灌區(qū)1981—2010年地下水位平均埋深基本呈增加趨勢，突變點(diǎn)出現(xiàn)在1994年，1994年之前影響灌區(qū)地下水埋深動態(tài)的主要因素為蒸發(fā)量，1994年之后主要影響因素為蒸發(fā)量、地下水開采量和灌渠引水量。王默涵等[20]采用地理探測器定量分析了邯鄲平原區(qū)的主要影響因子為工業(yè)用水變化和降水量變化。第3種方法是采用當(dāng)前較為熱門的機(jī)器學(xué)習(xí)模型進(jìn)行分析。關(guān)于“冀棗衡”漏斗影響因子的研究以定性研究為主，尹新明[21]認(rèn)為超采是“冀棗衡” 漏斗形成的主要原因。王建恒[22]分析了氣溫、日照、蒸發(fā)量等氣象要素對降水及衡水地區(qū)地下水水位變化的影響。

鑒于目前關(guān)于“冀棗衡”漏斗影響因子的研究以定性分析為主，定量分析不同時期主控因子的貢獻(xiàn)率的研究和研究地質(zhì)因子對漏斗貢獻(xiàn)的相對較少，為此本次采用多元線性回歸這種較為成熟的研究方法，將地下水開采和水文地質(zhì)因子作為影響因子，進(jìn)行地下水水位降落漏斗的影響因子分析，有助于我們更全面地了解“冀棗衡”地區(qū)地下水水位的動態(tài)變化，以期為未來地下水資源的合理開發(fā)和保護(hù)提供科學(xué)依據(jù)。

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

“冀棗衡”漏斗是深層地下水開采形成的漏斗，位于冀州—棗強(qiáng)—衡水一帶，地下水開采深度150～160 m，開采層位為第三含水層組，形成于1969年。1969年衡水市附近深層地下水位埋深約2.9 m，“冀棗衡”地區(qū)地下水位埋深約9.3 m，漏斗中心位于衡水市衡62觀測孔，漏斗面積只有幾平方千米，受工業(yè)和農(nóng)業(yè)集中開采影響，漏斗面積不斷加深和擴(kuò)大，到1973年5月底，水位埋深大于10 m，范圍已達(dá)2 980 km2。近十幾年超量開采深層水，致使水位逐年下降。2004年，漏斗北部已擴(kuò)展至衡水地區(qū)外，東南部已延伸至邢臺、滄州、德州地區(qū)，形成大型的復(fù)合型深層水位降落漏斗，漏斗面積為8 815 km2，覆蓋整個衡水地區(qū)。2016年至今，受景縣及鄰近山東省德州市地下水集中開采的影響，漏斗中心位于景縣留智廟鎮(zhèn)八里莊村，中心水位埋深在120.86～129.27 m之間；2020年漏斗中心移向西南故城縣房莊鄉(xiāng)堤口村，中心水位埋深為130.31 m。

表1是“冀棗衡”漏斗的統(tǒng)計數(shù)據(jù)。圖1是“冀棗衡”漏斗中心水位變化圖，從圖1可以看出1975—2020年漏斗中心水位持續(xù)下降，水位埋深范圍為32.44～130.31 m，下降幅度為97.87 m，水位下降速率為2.11 m/a。漏斗的影響面積自1975年持續(xù)擴(kuò)大，至2004年影響面積擴(kuò)展到整個衡水市，2004年以后影響面積的統(tǒng)計數(shù)據(jù)為8 815 km2(衡水市面積)。圖2為漏斗封閉面積歷史變化圖，從圖中可以看出漏斗封閉面積總體變化趨勢為波動中持續(xù)擴(kuò)大。近幾十年“冀棗衡”漏斗的位置不斷的變化，由西南向東北方向運(yùn)移。漏斗封閉區(qū)域歷史范圍見圖3。從圖3可以看出，1975年漏斗位于衡水市桃城區(qū)；1984年漏斗位于南宮市和棗強(qiáng)縣；1992年漏斗范圍較大，包括武邑縣、景縣、桃城區(qū)、棗強(qiáng)縣和南宮市；2005年漏斗范圍最大，從南到北，包括了邢臺南宮市、衡水市和滄州獻(xiàn)縣、東光縣等地；2016年漏斗范圍較小，西側(cè)桃城區(qū)、武邑縣等地不再位于封閉漏斗范圍內(nèi)，漏斗主要包含棗強(qiáng)縣、固城縣和景縣；2020年漏斗向東北發(fā)展，包含武邑縣、阜城縣、景縣、吳橋縣和固城縣。

圖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 “冀棗衡”漏斗統(tǒng)計表Table 1 Statistical table of “JZH” funnel

“冀棗衡”地區(qū)4期的地下水水位變幅數(shù)據(jù)(1985—1992年、1992—2005年、2005—2016年、2016—2020年)如圖4所示。從圖中可以看出，1984—1992年水位下降較大的地區(qū)位于棗強(qiáng)縣和桃城區(qū)，水位下降幅度為30 m左右，深州市的西北側(cè)水位有一定的回升，上升幅度為20 m左右；1992—2005年水位下降較大的地區(qū)位于景縣的北邊和南邊，水位變幅為44 m左右，整個衡水地區(qū)地下水水位變幅均為負(fù)值，說明區(qū)域水位均下降；2005—2016年水位下降較大的地區(qū)位于景縣東邊，水位變幅為26 m左右，深州市和武強(qiáng)縣水位出現(xiàn)了一定的回升，其中深州市的回升幅度為10 m左右，武強(qiáng)縣為4 m左右；2016—2020年水位下降較大的地區(qū)位于冀州區(qū)和衡水市的西北角，水位變幅為35 m左右，其次為故城縣和阜城縣，水位變幅為25 m左右，與阜城縣相鄰的滄州泊頭市水位出現(xiàn)了一定的回升，回升幅度為20 m左右。

總體來說，1984—2020年“冀棗衡”地區(qū)的地下水水位總體呈下降趨勢。1992—2005年是水位下降幅度最大的時間段，區(qū)域水位整體下降，而其他3個階段均存在局部水位上升。這4個階段水位變幅由大到小依次為：1992—2005年、2016—2020年、1985—1992年、2005—2016年。這4個階段水位變化速率由大到小依次為：2016—2020年、1985—1992年、1992—2005年、2005—2016年。

2 漏斗發(fā)展演化的影響因子分析

已有的研究表明衡水地區(qū)深層地下水的來源主要有側(cè)向徑流補(bǔ)給、越流補(bǔ)給和黏性土壓密釋水[23]。其中淺層水對深層水的越流補(bǔ)給主要發(fā)生在第一含水組和第二含水組之間，第三、第四含水組主要來源于儲存量和側(cè)向徑流補(bǔ)給量[24]。深層地下水的主要排泄項(xiàng)為地下水開采。深層地下水開采是影響漏斗發(fā)展的主要因素，通過對比開采量與水位埋深(圖5)，可知開采與埋深變化趨勢總體一致。兩者之間進(jìn)行相關(guān)性分析的結(jié)果顯示(圖6)，兩者之間具有良好的相關(guān)性，相關(guān)系數(shù)R2為0.804 5。

圖6 開采量與“冀棗衡”漏斗中心埋深相關(guān)性分析圖Fig.6 Correlation analysis between mining volume and burial depth of “JZH” funnel center

一般來說，降雨量比較大的年份開采量相對較小，從而導(dǎo)致地下水水位下降速率較慢，由此說明降雨量對地下水水位是間接影響。對比降雨量與開采量(圖7)和降雨量與漏斗中心水位(圖8)。從圖7可以看出，降雨量與開采量在1990年前有一定的對應(yīng)關(guān)系，1990年前隨著降雨量的增加，開采量有所減少。但是從圖8 可以看出 1990 年后，水位仍舊持續(xù)下降，降雨量相對穩(wěn)定，開采量持續(xù)升高。由此說明，降雨量在1990年前對開采量有所影響，但是影響強(qiáng)度較小，水位持續(xù)下降，1990年后降雨量對開采量影響較弱，水位持續(xù)下降。利用1990年前后降雨量與開采量數(shù)據(jù)進(jìn)行相關(guān)性分析(圖9)，1990年前兩者之間具有良好的負(fù)相關(guān)性，相關(guān)系數(shù)R2為0.76。1990年后兩者之間相關(guān)性較差。由此說明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 開采量與降雨量相關(guān)性分析圖Fig.9 Analysis of correlation between mining and rainfall

3 漏斗空間水位變化的影響因素分析

本次以“冀棗衡”漏斗為研究對象。利用ArcGIS軟件在研究區(qū)范圍內(nèi)隨機(jī)生成620個點(diǎn)，提取這些點(diǎn)的彈性釋水系數(shù)數(shù)據(jù)、滲透系數(shù)數(shù)據(jù)、底板埋深、含水層厚度數(shù)據(jù)、水位變幅數(shù)據(jù)(1985—1992年、1992—2005年、2005—2016年、2016—2020年)。采用SPSS軟件進(jìn)行回歸分析，分析4個時間階段水文地質(zhì)要素對水位變幅的影響強(qiáng)度。結(jié)合開采量數(shù)據(jù)分析回歸分析結(jié)果的合理性。其中彈性釋水系數(shù)數(shù)據(jù)利用河北平原的1 051個鉆孔數(shù)據(jù)計算每個鉆孔彈性釋水系數(shù)，然后利用ArcGIS對點(diǎn)數(shù)據(jù)進(jìn)行克里金插值，利用插值后柵格數(shù)據(jù)提取本次研究點(diǎn)的彈性釋水系數(shù)數(shù)據(jù)；滲透系數(shù)的數(shù)據(jù)利用鉆孔的抽水試驗(yàn)獲得的滲透系數(shù)進(jìn)行克里金插值后獲得；頂板埋深采用第一含水層底板埋深柵格數(shù)據(jù)獲得；含水層厚度數(shù)據(jù)采用第一含水層底板埋深與第四含水層底板埋深進(jìn)行算術(shù)計算獲得；水位變幅數(shù)據(jù)利用1982年、1992年、2005年、2016年和2020年水位柵格數(shù)據(jù)進(jìn)行計算獲得。

3.1 “冀棗衡”地區(qū)各特征因子分布特征

“冀棗衡”地區(qū)的滲透系數(shù)和彈性釋水系數(shù)普遍偏低(圖10、圖11)，滲透系數(shù)數(shù)值范圍為3.99～14.4 m/d，冀州和景縣東北地區(qū)相對較高，其他地區(qū)偏低。彈性釋水系數(shù)的數(shù)值范圍為0.000 1～0.002 5，吳橋和深州等地較高，其他地區(qū)偏低?！凹綏椇狻钡貐^(qū)深層水的底板標(biāo)高-450 m左右，頂板標(biāo)高為20 m左右，深層水厚度400多米。

圖10 研究區(qū)滲透系數(shù)分布圖Fig.10 Spatial distribution of infiltration coefficient in the study area

圖11 研究區(qū)彈性釋水系數(shù)分布圖Fig.11 Spatial distribution of elastic water release coefficients in the study area

圖12 地質(zhì)要素對“冀棗衡”漏斗影響的解釋強(qiáng)度Fig.12 Interpreted intensity of the influence of geological elements on “JZH” funnel

3.2 水文地質(zhì)指標(biāo)對水位變化的貢獻(xiàn)率分析

通過對研究區(qū)620個點(diǎn)的含水層厚度數(shù)據(jù)、彈性釋水系數(shù)數(shù)據(jù)、滲透系數(shù)數(shù)據(jù)、含水層頂板埋深、1984—1992年水位變幅數(shù)據(jù)進(jìn)行分析，采用spss軟件進(jìn)行回歸分析，結(jié)果見表2、表3、表4。4期分析數(shù)據(jù)的Durbin—Watson系數(shù)均接近2，滿足樣本獨(dú)立要求。共線性診斷的顯著性P值為0.000a，小于0.05，說明自變量之間相互獨(dú)立，滿足計算模型要求。共線性統(tǒng)計量VIF均小于5，說明模型計算結(jié)果可靠。

表2 模型匯總b表Table 2 Summaryb of the mldel

表3 Anovab表Table 3 Table of Anovab

表4 系數(shù)a表Table 4 Table of coefficientsa

“冀棗衡”漏斗4個階段的分析可知，1984—1992年、1992—2005年、2005—2016年、2016—2020年水文地質(zhì)指標(biāo)對水位變幅的解釋度R2依次為：0.382、0.212、0.045、0.163(表2)。1984—1992年對水位變幅影響較大的水文地質(zhì)指標(biāo)依次為頂板埋深、彈性釋水系數(shù)；1992—2005年對水位變幅影響較大的水文地質(zhì)指標(biāo)依次為彈性釋水系數(shù)、頂板埋深、含水層厚度；2005—2016年對水位變幅影響較大的水文地質(zhì)指標(biāo)依次為含水層厚度、頂板埋深、滲透系數(shù)；2016—2020年對水位變幅影響較大的水文地質(zhì)指標(biāo)依次為含水層厚度、頂板埋深、滲透系數(shù)、彈性釋水系數(shù)。

一般認(rèn)為漏斗發(fā)展的主控因子是地下水開采，隨著地下水開采量的增加，開采對地下水水位下降的影響強(qiáng)度逐漸增加，從而導(dǎo)致地質(zhì)要素對漏斗的影響強(qiáng)度逐漸降低。結(jié)合圖7可知，1984—1990 年開采量逐漸減少；1990—2005 年是開采量快速增長；2005—2016 年開采量在波動中繼續(xù)增長，2016—2020 年開采量逐漸下降。 與地質(zhì)要素對水位變化解釋率的變化具有良好負(fù)相關(guān)。由此說明，回歸分析結(jié)果符合一般認(rèn)識。

4 結(jié)論

近45年來“冀棗衡”漏斗中心水位持續(xù)下降，多年平均下降速率為2.1 m/a，漏斗面積持續(xù)擴(kuò)大，自2014年漏斗面積已近擴(kuò)大至衡水全境。漏斗演化特征為由西北向東南方向發(fā)展，漏斗中心由早期的衡水市棗強(qiáng)縣，向東南部的景縣、固城縣發(fā)展。4個時期的水位下降幅度顯示，1985—1992年水位下降幅度較大地區(qū)為桃城區(qū)和棗強(qiáng)縣，水位降幅約25 m；1992—2005年為景縣和固城縣，降幅約44 m；2005年后水位下降速率減小，2005—2016年為景縣，降幅約26 m；2016—2020年為固城縣，降幅約25 m。

地下水開采是影響漏斗的主控因子，漏斗中心水位埋深與開采量具有良好的相關(guān)性，1990年前降雨量與開采量具有良好的相關(guān)性，1990年后降雨量與開采量之間的相關(guān)性較小。應(yīng)用回歸分析和GIS技術(shù)分析漏斗區(qū)水位變化影響因子，結(jié)果顯示水文地質(zhì)因子對研究區(qū)地下水水位空間變化的影響解釋率為：1984—1992年為38.2%、1992—2005年為21.2%、2005—2016年為4.5%、2016—2020年為16.3%，該結(jié)果與研究區(qū)開采強(qiáng)度具有良好的對應(yīng)關(guān)系，隨著開采強(qiáng)度的增加地質(zhì)要素對漏斗影響強(qiáng)度在降低，符合一般認(rèn)識。研究成果可以為衡水地區(qū)的地下水修復(fù)提供技術(shù)參考。