河南某大型裂隙巖溶水源地地下水位動態分析

姜寶良，陳寧寧，李小建，潘 登，張曉偲，姜英博

（1.華北水利水電大學，河南 鄭州 450046；2.濟源市布袋溝供水站，河南 濟源 459000；3.河南省地質環境監測院，河南 鄭州 450016）

河南某大型裂隙巖溶水源地位于濟源克井盆地出口青多村東南，為濟源市第一水廠和第二水廠提供水源。自1994年投入運行以來，對地下水位、開采量、水質等進行了系統的統計和監測工作，積累了豐富的監測資料。隨著氣象、水文以及開采條件的變化，水源地的地下水資源也隨之發生變化，這些變化集中反映在水位動態上。

地下水動態是水文地質學的重要研究內容，早在20世紀40年代國內外相關學者就開始研究地下水動態[1]。目前地下水動態分析方法可分為解析法、數值法和數理統計法3 類。解析法和數值法必須了解邊界條件和相關水文地質參數。裂隙巖溶地下水的水文地質條件復雜，很難取得所需參數，難以建模求解。數理統計的逐步回歸分析方法，可以從眾多影響因素中，智能篩選出較為重要的影響因素，建立水位或泉水流量與相關因素（變量）的相關方程，具有一定實用價值。國內學者對華北平原地下水的動態類型、水位時空變化、下降速率及影響因素等開展研究分析[2-7]，這些研究主要針對華北平原大區域，數據的空間尺度和時間跨度都較大，對局部區域分析深度不夠。姜寶良，趙貴章等[8-9]利用逐步回歸分析方法，研究新鄉百泉泉水流量動態，確定了泉水流量動態的主要影響因素，建立了天然條件下泉水流量動態預測的逐步回歸模型，預測了天然狀態下的泉水流量，評價了百泉泉域的裂隙巖溶水資源量；姜寶良等[10]通過某大型水源地的地下水位動態分析，建立了水位與開采量、降水量等相關因素逐步回歸模型，進行了水位預測和資源評價。國內外針對地下水位動態特征及其影響因子的相關研究多以城市、流域、灌區等較大面積的區域為研究對象[11-12]，較小尺度范圍內的地下水位變化研究較少且結果不夠精確[13-14]。蟒河口水庫建設過程中，部分防滲工程未隨水庫主體工程完工，導致水庫建成后漏水現象嚴重。蟒河口水庫滲漏對水源地的地下水環境影響尚未有人研究。本文根據水源地水位與開采量、降水量、蟒河口水庫的蓄水量、水庫蓄水后水位及水源地與水庫水質等資料，定量確定地下水位動態的影響因素，論述了蟒河口水庫蓄水后對水源地的補給效果。

1 研究區概況

克井盆地北、西和南部邊界由盤古寺斷層、三樊斷層、封門口斷層和五龍口斷層構成。盆地北部為東西向延伸的太行山南麓余脈，西部萬羊山及東南部孔山為低山丘陵，盆地內由蟒河沖洪積扇和山前坡洪積扇裙組成。地形東高西低，北高南低。

克井盆地含水巖組分為3 類：松散巖類孔隙含水巖組、碳酸鹽巖類裂隙巖溶含水巖組和碎屑巖類裂隙含水巖組[15]。因碳酸鹽巖類裂隙巖溶含水巖組為水源地開采目的層，故主要介紹此類含水巖組。根據克井盆地碳酸鹽巖裂隙巖溶含水巖組的埋藏條件可分為：裸露型、覆蓋型和埋藏型。太行山南麓、孔山、萬羊山碳酸鹽巖大面積裸露（裸露型），巖性主要為寒武—奧陶系的灰巖、白云質灰巖、白云巖等，加之斷層交錯，裂隙巖溶發育。五龍口斷層與封門口斷層之間、盆地出口的萬羊山與孔山之間，碳酸鹽巖被第四系松散堆積物覆蓋（覆蓋型），局部奧陶系灰巖與第四系卵礫石直接接觸，形成“天窗”。盆地內石炭—二疊系砂泥巖地層掩埋（埋藏型）（圖1）。克井盆地裂隙巖溶地下水主要接受大氣降水入滲補給，沁河側向徑流補給和蟒河口水庫滲漏補給。地下水排泄以人工開采為主，也有以泉的形式集中排泄。

圖1 克井盆地水文地質條件略圖Fig.1 Simplified hydrogeological map of the Kejing Basin

2 水源地地下水位動態

青多水源地是濟源市唯一的大型裂隙巖溶地下水供水水源地，其多年地下水位動態可分為2 個階段（圖2）。

圖2 水源地水位與開采量、降水量、蟒河口水庫蓄水量動態曲線圖Fig.2 Dynamic curves of groundwater levels,extraction rate,precipitation and water storage of the Manghekou Reservoir

2.1 第一階段（1994年6月—2014年8月）

蟒河口水庫建成前，水源地水位主要受降水和開采影響，水位一般滯后降水3～4 個月。根據水源地開采量的變化又分2 個時段。

2.1.1 第1 時段（1994年6月—2005年12月）

水源地開采量維持在約2×104m3/d，地下水位主要受降水量控制，呈現豐水期（年）水位上升，枯水期（年）水位下降的動態特征，開采量影響較小，基本呈現天然動平衡狀態，水位峰值滯后降水峰值2 個月以上。如：枯水年（1997年）降水量僅348.5 mm，水位持續下降。從年初的162.8 m 下降到年末的153.01 m，下降近10 m，呈負均衡；豐水年（1996年）降水量932 mm，7—9月降水量達742.3 mm，水位從7月的151.99 m 迅速上升到9月的163.09 m，呈正均衡狀態；平水年呈均衡狀態。該時段枯水年少，豐水年和平水年多，開采量小，總體呈正均衡狀態。

2.1.2 第2 時段（2006年1月—2014年8月）

水源地開采量由2.2×104m3/d 逐漸增大到約6.5×104m3/d，地下水位受降水和開采雙重因素控制，呈現持續下降趨勢。僅在豐水年（2011年降水量897.5 mm，7—9月降水量達620.5 mm）水位由7月的146.6 m 急劇上升到11月的157.15 m，之后又持續下降至2014年6月的148.15 m。該時段由于開采量持續增加，水位持續下降，總體呈負均衡狀態。

2.2 第二階段（2014年9月—2019年12月）

2014年9月10 日蟒河口水庫開始蓄水。2015年汛期沒有蓄水。2016—2019年連續4年汛期蓄水。

水源地開采量由2014年8月的6.5×104m3/d 增大到2019年6月的8.8×104m3/d。蟒河口水庫蓄水后雖水源地開采量不斷增加，但水源地水位不降反升，出現新的均衡狀態。2019年7月底，濟源第三水廠（利用河口村水庫水）建成供水，水源地開采量開始減少，至2019年12月開采量減少到約6×104m3/d，水源地水位上升。由于蟒河口水庫距水源地約7 km，其水位滯后蟒河口水庫蓄水高峰2～3 個月。

3 水源地地下水位動態逐步回歸分析

以水源地水位H為因變量，降水量P、開采量Q、蟒河口水庫蓄水量W及前期水位等為自變量，進行逐步回歸計算。

3.1 數據和算法

（1）數據

收集1994年以來水源地地下水位、降雨量、開采量、水庫水位及2014年開始蓄水后的水庫蓄水量數據：

①水源地水位H選用觀測系列長且相距較近的C3、C4、C5 井的逐月平均水位。

②降水量P選用水文部門布置在克井盆地及附近的交地、竹園、黃龍廟3 個雨量站和五龍口、濟源2 個水文站逐月降水量的算術平均值。

③水源地開采量Q采用濟源市自來水公司統計的逐月平均開采量[16-18]。

④蟒河口水庫蓄水量W采用逐月平均蓄水量。

（2）逐步回歸分析方法

逐步回歸分析是在多元回歸分析的基礎上派生出來的一種算法，能從眾多變量（或預先盡可能多地考慮一些變量）中自動挑選重要變量[8，10，16]，并確定其數學表達式的一種統計方法。逐步回歸模型為：

式中：Hi—當月地下水位/m；

x1、x2、···、xn—自變量；

B0—常數項；

bi（i=1，2，···，n）—回歸系數。

該方法常用來分析眾多因素中無法確定究竟是哪些因素對因變量有影響的情況，在數據分析中應用廣泛。其最大特點是在引入或剔除一個變量時都要進行F檢驗。具體流程為根據向前選擇法則引入自變量，然后根據向后剔除法則將回歸方程中F檢驗概率最小且符合剔除條件的自變量剔除，重復進行上述引入和剔除過程直到回歸方程外的自變量均不符合引入條件、回歸方程中的變量均不符合剔除條件。

計算結果中，除式（1）中出現的因素外，還出現各影響因素F檢驗值F，F值越大，該因素對因變量影響越大；復相關系數R，R越大，方程回歸效果越顯著。

3.2 結果與分析

3.2.1 蟒河口水庫蓄水前（1994.6—2014.8）

以當月水位Hi為因變量，以水源地當月開采量Qi、當月及前1～4月降水量（Pi、Pi-1、Pi-2、Pi-3、Pi-4）、前期水位（Hi-1）為自變量，對1994年6月—2014年8月數據進行逐步回歸分析計算，F臨界檢驗值取2.0，計算得到回歸方程為：

該回歸方程的復相關系數R=0.88；剩余標準差SY=1.36；F*=263.70。

由式（2）可知，水源地水位主要受前期水位、前2月降水量和當月開采量影響，其中前期水位的影響最為顯著（F=610.6），其次是前2月降水量（F=46.8）；最后是當月開采量（F=4.9）。

3.2.2 蟒河口水庫蓄水后（2014.9—2019.12）

在求解該階段逐步回歸方程的基礎上，探明前期水位（Hi-1）對水源地地下水位的影響，故分兩種情況進行討論：

①以水源地水位Hi為因變量，以水源地當月開采量Qi、當月及前1～4月降水量（Pi、Pi-1、Pi-2、Pi-3、Pi-4）、當月及前1～4月蟒河口水庫蓄水量（Wi、Wi-1、Wi-2、Wi-3、Wi-4）為自變量，采用不同的F檢驗值分別進行逐步回歸計算。

②以水源地水位Hi為因變量，以水源地當月開采量Qi、當月及前1～4月降水量（Pi、Pi-1、Pi-2、Pi-3、Pi-4）、當月及前1～4月蟒河口水庫蓄水量（Wi、Wi-1、Wi-2、Wi-3、Wi-4）、前期水位Hi-1為自變量，采用不同的F檢驗值分別進行逐步回歸計算。

計算結果見表1。

由表1可知，蟒河口水庫蓄水后，在不考慮前期水位的影響下，相關性較為一般（R=0.595 5～0.746 6）；加入前期水位Hi-1這一因素后，方程相關性明顯提升（R=0.856 6～0.918 3）。由此可見水源地地下水位受前期水位影響最為顯著，其次是前1月降水、前2月蟒河口水庫蓄水量，雖開采量由6.5×104m3/d 增加到最大的8.77×104m3/d，其對水位基本無影響或影響甚微。

表1 水庫蓄水后逐步回歸方程Table 1 Stepwise regression equations after reservoir impoundment

由于水庫蓄水后的逐步回歸方程基本都無開采量的體現，但實際預測過程中必須考慮開采量的影響，故不使用第二階段的回歸方程作為水位預測方程。

4 蟒河口水庫補給效果分析

蟒河口水庫位于克井盆地西北的北蟒河出山口，控制流域面積94 km2，總庫容為1 094×104m3，主要任務是防洪、補給地下水等，蟒河口水庫建成后，其防洪效果明顯，但補給效果備受詬病。

蟒河口水庫壩址區為寒武系中上統的碳酸鹽巖地層，受盤古寺斷層影響，次級斷裂發育，巖體破碎，裂隙巖溶發育，滲漏嚴重。水庫蓄水水位標高為255～310 m；水源地地面標高約173 m，水位埋深10～25 m，水位標高148～163 m，水庫蓄水水位高于水源地水位100 m 以上。蟒河口水庫距水源地直線距離約7 km，其間在石炭—二疊系地層之下隱伏有巨厚層的奧陶—寒武系的碳酸鹽巖裂隙巖溶含水層，蟒河口水庫滲漏補給碳酸鹽巖裂隙巖溶地下水，通過裂隙巖溶通道由北向南徑流補給水源地地下水[19]。

利用蟒河口水庫蓄水前建立的逐步回歸方程式（2），預測在不考慮蟒河口水庫蓄水影響下水源地水位（2014年9月—2019年12月），見圖3。實測水位與預測水位差值見表2。

由圖3和表2可知，蟒河口水庫蓄水期間實測水位均高于不考慮蟒河口水庫蓄水的預測水位，其差值最大1.46～3.10 m，最小0.21～0.89 m，平均0.86～1.57 m。說明蟒河口水庫蓄水對水源地裂隙巖溶地下水補給效果顯著。

圖3 逐步回歸方程（2）預測值與實測水位動態曲線圖Fig.3 Dynamic curves of the predicted values of the stepwise regression equation（2）and the measured groundwater levels

表2 蟒河口水庫蓄水期間水源地實測水位與不考慮蟒河口水庫蓄水預測水位差值Table 2 Differences between the measured groundwater levels near the wellfield during the impoundment period of the Manghekou Reservoir and the predicted groundwater levels without considering the impoundment of the Manghekou Reservoir /m

從2018年蟒河口水庫與青多水源地水質分析結果（表3）可知，蟒河口水庫和水源地的水化學類型均為HCO3·SO4—Ca·Mg 型，水源地主要化學組分除K+外，其它含量均大于蟒河口水庫，這是由于蟒河口水庫水滲漏補給裂隙巖溶地下水后，經過較長距離的運移和溶濾作用，溶解了碳酸鹽巖礦物質，使水源地地下水中的Ca2+、Mg2+等離子的含量增加，總硬度和溶解性總固體也增加。

表3 蟒河口水庫與青多水源地水質分析對比表Table 3 Comparison of hydrochemical analyses between the Manghekou Reservoir and the Qingduo wellfield

5 結論

（1）第一階段水位受前期水位的影響最顯著、其次是前2月降水量、開采量影響最小。其中第1 時段水源地開采量固定，地下水位主要受水源地降水影響，開采量影響很小；第2 時段增大開采量，水位主要受開采影響，其次是降水量。

（2）第二階段水位受前期水位影響最為顯著，前1月降水量、前2月水庫蓄水量次之。雖開采量增加較多，但其對水位的影響很小或基本沒有影響，其水位滯后蟒河口水庫蓄水1～4 個月、降水1～3月。

（3）第二階段（水庫蓄水后）計算得出的回歸方程中，開采量對地下水位無影響或產生微弱影響，無法對實際工作中考慮開采量條件的地下水位進行預測，不具有現實意義，故不采用該階段方程作為預測方程。

通過分析逐步回歸方程的預測水位與水源地實測水位差、水源地與水庫水質發現，蟒河口水庫蓄水后對水源地地下水補給效果顯著，但本研究在此方面只進行了定性研究，今后需要從同位素示蹤方向進一步研究。