劉中培，齊明坤，韓宇平，曹潤祥，冷靜

(華北水利水電大學 a.水資源學院；b.河南省黃河流域水資源節(jié)約集約利用重點實驗室，鄭州 450046)

地下水作為重要的水資源，在維持社會經(jīng)濟發(fā)展、生態(tài)環(huán)境和農(nóng)業(yè)生產(chǎn)需求方面發(fā)揮著重要作用[1].地下水開采量的不斷增大導致地下水漏斗的產(chǎn)生及發(fā)展，并在某些地區(qū)出現(xiàn)了難以恢復的永久性降落漏斗[2-4].地下水漏斗的發(fā)展會進一步引發(fā)諸如地面沉降、地裂縫、咸水入侵、地下水污染等環(huán)境地質(zhì)問題[5-6]，影響區(qū)域生態(tài)保護和高質(zhì)量發(fā)展[7].目前關于地下水漏斗的研究主要集中在漏斗演變特征及修復等方面，而對地下水漏斗演變趨勢的預測研究相對較少[8-11].

人民勝利渠灌區(qū)位于黃河下游，是河南省重要的糧食生產(chǎn)區(qū)，作為重要的灌溉水源，地下水在農(nóng)業(yè)發(fā)展中扮演著重要的角色[12].長期大量的開采地下水導致灌區(qū)地下水位普遍下降，出現(xiàn)了地下水降落漏斗并不斷擴展.以往的研究主要集中在水資源管理、利用，以及優(yōu)化配置等方面，有關灌區(qū)地下水漏斗方面的研究較少[13-15].本文采用地理空間分析方法識別了人民勝利渠灌區(qū)地下水漏斗的形成時間，分析了地下水漏斗的中心水位、面積演變特征[16-17]，構建了地下水漏斗的中心水位模型和面積模型對漏斗中心水位和面積進行預測，從而為實現(xiàn)地下水的科學管理與區(qū)域高質(zhì)量發(fā)展提供依據(jù).

1 研究區(qū)概況

人民勝利渠灌區(qū)(113°31′E～114°25′E，35°0′N～35°30′N)位于黃河北岸，總面積為1 486.84 km2.灌區(qū)屬于溫帶大陸性季風氣候，四季分明，春季干旱多沙、夏季炎熱多雨、秋季秋高氣爽、冬季寒冷干燥.多年平均氣溫約14 ℃，無霜期220 d左右.根據(jù)國家氣象中心新鄉(xiāng)站多年的降水資料，灌區(qū)多年平均降水量為620 mm，年內(nèi)降水量分配不均，主要集中在6至9月份.多年平均潛在蒸發(fā)量為1 864 mm.

根據(jù)《河南省新鄉(xiāng)縣區(qū)域水文地質(zhì)調(diào)查》成果，研究區(qū)地層可分為第三系上新統(tǒng)，以棕紅黏土、粉質(zhì)黏土為主；第四系下更新統(tǒng)，為冰積地層；第四系中更新統(tǒng)，為洪積地層；第四系上更新統(tǒng)下部，為沖積地層；全新統(tǒng)一段為沖積、洪積地層.全新統(tǒng)二段為沖積、洪積、風積地層.全新統(tǒng)三段為沼澤堆積地層.依據(jù)含水介質(zhì)及類型，區(qū)域調(diào)查深度范圍內(nèi)地下水可劃分為松散巖類孔隙水，半膠結碎屑巖類孔隙裂隙水兩類.根據(jù)含水層埋深條件可分為淺層含水組、中深層含水組和深層含水組.按5 m的降深計算又可分為強富水區(qū)、中富水區(qū)、弱富水區(qū).其中灌區(qū)西部關堤一帶為強富水區(qū)(>3 000 m3/d)，灌區(qū)中部獲嘉一帶為中富水區(qū)(1 000～3 000 m3/d)，灌區(qū)北部共產(chǎn)主義渠附近為弱富水區(qū)(<1 000 m3/d).研究區(qū)概況見圖1.

2 材料與方法

2.1 數(shù)據(jù)來源及處理

本文分析對象為淺層地下水，地下水埋深資料來源于河南省人民勝利渠管理局.收集了研究區(qū)內(nèi)24個觀測井(由于個別觀測井經(jīng)緯度數(shù)據(jù)缺失，某些年份觀測井不到24個)1953-2017年逐月地下水埋深觀測數(shù)據(jù)，地下水水位為觀測井地面高程與埋深之差；地下水漏斗面積數(shù)據(jù)獲取首先通過ArcGIS軟件繪制研究區(qū)地下水等水位線圖，然后識別地下水漏斗的范圍，最后利用空間分析提取數(shù)據(jù)實現(xiàn).

本文采用基于經(jīng)驗模態(tài)分解(Empirical Mode Decomposition,EMD)處理后的數(shù)據(jù)構建漏斗中心水位自回歸(Autoregression,AR)模型，進行中心水位的預測；結合普通最小二乘(OLS)方法，建立漏斗面積與漏斗中心水位變化的OLS回歸模型，預測漏斗面積.

2.2 研究方法

2.2.1基于EMD的AR模型

EMD是對一種信號進行平穩(wěn)化處理，將信號中的不同尺度的波段或趨勢項提取出來，從而得到一系列的不同特征尺度的經(jīng)驗模態(tài)函數(shù)(Intrinsic Mode Function,IMF)和一個殘余項[18-19].基于EMD的AR模型建模步驟具體如下所示.

步驟1 提取原始時間信號的IMF和殘余項.提取IMF過程中，IMF必須符合兩個條件：I(t)極值的數(shù)量(最大值和最小值數(shù)量之和)與零穿越的數(shù)量必須相等或最多相差1；在I(t)的任意點，局部最大值定義的包絡線的平均值和局部最小值定義的包絡線的平均值應該等于零，提取過程可以通過“篩選”步驟來完成[20-21].原始時間序列D(t)表示為：

(1)

式中：D(t)表示原始時間序列，rn表示最終的殘余項，Ii表示不同頻率的IMF.

步驟2 對所有的經(jīng)驗模態(tài)函數(shù)IMFs分別建立AR模型[22-23].AR模型計算流程如圖2所示，模型數(shù)學表達式為：

Ii,j=μ+φ1(Ii,j-1-μ)+φ2(Ii,j-2-μ)+…+

φp(Ii,j-p-μ)+εt,

(2)

式中：Ii,j為第i個IMF第j時刻序列值，μ為序列均值，φ1,φ2,…,φp為自回歸系數(shù)，εt為白噪聲項.

2.2.2平穩(wěn)性及白噪聲性檢驗

AR模型要求時間序列是平穩(wěn)的非白噪聲序列，因此需對EMD分解的各經(jīng)驗模態(tài)函數(shù)IMFs先進行平穩(wěn)性檢驗再進行白噪聲性檢驗.數(shù)據(jù)的平穩(wěn)性采用單位根來進行檢驗，若不存在單位根，則時間序列平穩(wěn)，若存在單位根，則時間序列不平穩(wěn)[24].對于不平穩(wěn)的時間序列采用差分處理，直到時間序列平穩(wěn)，通過單位根檢驗[24-25].序列的白噪聲性檢驗采用Ljung-Box Q檢驗的方法，用以檢驗時間序列是否存在自相關，若時間序列存在自相關，則序列為非白噪聲性序列，否則為白噪聲性序列[26].上述檢驗過程均由Eviews專業(yè)統(tǒng)計學軟件完成.

2.2.3模型階數(shù)和參數(shù)估計

模型階數(shù)的確定是模型參數(shù)(φi,μ,εt)估計的前提條件，此外模型適用性檢驗的核心就是解決模型的定階問題.本文采用赤池信息準則(Akaike Information Criterion,AIC)確定AR模型階數(shù)[27].

(3)

(4)

2.2.4模型的診斷與修正

對各經(jīng)驗模態(tài)函數(shù)IMFs建立的AR模型進行殘差序列相關檢驗和殘差異方差檢驗，通過一系列檢驗來驗證模型的準確性.對沒有通過檢驗的模型進行修正，使其通過檢驗，以實現(xiàn)對模型的應用.本文采用拉格朗日乘數(shù)檢驗(Lagrange Multiplier Test,LMT)對殘差序列相關性進行檢驗，采用懷特檢驗(White Test,WT)對殘差異方差進行檢驗[28].

2.2.5OLS回歸模型

本文使用OLS分析漏斗面積與中心水位的關系[29].在地下水漏斗演變過程中，漏斗中心水位和漏斗面積關聯(lián)性較大，通過夏莊漏斗處的中心水位和漏斗面積的相關分析，二者呈現(xiàn)高度相關，且二次函數(shù)擬合時R2值最大，擬合程度最好.建立二者的回歸方程，數(shù)學表達式為：

Y/m2=c+b(X/m)+a(X/m)2,

(5)

式中：Y為因變量(漏斗面積)，X為自變量(中心水位)，c為常數(shù)，a,b分別為一次項和二次項系數(shù).

3 結果與分析

3.1 地下水漏斗演變分析

3.1.1漏斗形成及發(fā)展

為分析研究區(qū)內(nèi)地下水漏斗的產(chǎn)生時間及演變特征，根據(jù)1953-2017年地下水位埋深數(shù)據(jù)，利用ArcGIS地統(tǒng)計學模塊克里金(Kriging)插值法對地下水位進行空間插值，并繪制每年的地下水水位等值線圖，地下水位變化明顯的年份地下水水位情況如圖3所示.

由圖3(a)和圖3(b)對比可知，1975年以前灌區(qū)未出現(xiàn)地下水降落漏斗(圖3(a))，1976年在灌區(qū)的中部形成以董莊為中心的地下水降落漏斗，但地下水漏斗范圍較小(圖3(b)).隨后漏斗面積逐漸擴大，1990年董莊漏斗面積達到最大，如圖3(c)所示.圖3(d)和圖3(e)對比發(fā)現(xiàn)，2000年漏斗中心仍然在董莊(圖3(d))，與1990年相比只是漏斗面積有所變化，而至2001年，以夏莊為中心的地下水漏斗開始逐漸形成，但此時董莊漏斗和夏莊漏斗基本相對獨立，未形成統(tǒng)一的漏斗(圖3(e)).至2003年，隨著地下水位的持續(xù)下降，董莊漏斗和夏莊漏斗連通，形成了一個面積更大的地下水降落漏斗(圖3(f)).圖3(g)和圖3(h)表明，隨著時間的推移，降落漏斗中心逐漸轉移至夏莊附近，并且漏斗面積持續(xù)擴大至研究區(qū)外圍.

綜上所述，研究區(qū)地下水漏斗演變大致經(jīng)歷了三個階段.第一個階段為1975年以前，灌區(qū)未出現(xiàn)明顯的地下水降落漏斗，地下水位等值線分布相對均勻，地下水由灌區(qū)西南向東北方向流動.第二個階段為1975-2000年，地下水漏斗開始出現(xiàn)，此階段地下水漏斗中心在董莊附近.當?shù)叵滤┒烦霈F(xiàn)后，地下水流場發(fā)生了明顯變化，形成了以地下水漏斗為中心的局部地下水流動系統(tǒng).灌區(qū)地下水位等值線變得疏密不一，漏斗區(qū)等值線密集，水力梯度大，非漏斗區(qū)等值線稀疏，水力梯度相對較小.第三階段為2001年至現(xiàn)在，隨著地下水位的持續(xù)下降，漏斗面積逐漸增加，董莊漏斗和夏莊漏斗連通，逐漸形成了以夏莊為中心的更大的地下水降落漏斗，至2017年漏斗邊界已擴展至研究區(qū)外圍.

3.1.2漏斗中心水位及面積演變

漏斗中心水位及面積是衡量漏斗發(fā)展的重要表征，地下水降落漏斗中心及其水位變化反映了漏斗的位置及垂向演變特征；漏斗面積反映了地下水降落漏斗水平發(fā)展特征.研究區(qū)地下水漏斗演變特征如圖4所示.

由圖3和圖4可知，1990年之前，董莊地下水水位不斷下降，地下水漏斗面積逐年變大.1990年之后，漏斗中心水位逐年恢復，漏斗面積逐漸變小，整個變化過程較為平緩.1990-2000年董莊漏斗中心水位上漲約2 m，漏斗面積減小了約60 km2.2001年漏斗中心轉移至夏莊，形成以夏莊為中心的地下水降落漏斗.與董莊相比，夏莊漏斗中心水位與漏斗面積的變化情況大致為兩個階段.2001-2005年二者變化幅度較大，漏斗中心水位下降了約5 m，漏斗面積增加了150 km2.2006-2017年二者變化幅度明顯變小，漏斗中心水位下降了約2 m，面積增加了約5 km2.

綜上所述，董莊漏斗變化過程與夏莊漏斗變化過程有所不同，董莊漏斗中心水位先下降后上升，面積先增大后減小，兩者變化幅度相對較小，總體呈現(xiàn)負相關關系.夏莊漏斗中心水位和漏斗面積演變經(jīng)歷兩個階段，第一個階段為2001-2005年，變化明顯，水位迅速下降，面積快速增加；第二階段為2005年之后，兩者變化相對平穩(wěn)，總體呈現(xiàn)明顯負相關關系.

3.2 地下水漏斗模型構建及預測

3.2.1漏斗中心水位及面積模型構建

(1)漏斗中心水位預測模型構建

選取1953-2017年夏莊漏斗中心處的水位數(shù)據(jù)，利用EMD方法，將漏斗中心水位數(shù)據(jù)分解成4個經(jīng)驗模態(tài)函數(shù)IMFs，分解結果如圖5所示.人民勝利渠灌區(qū)的地下水位變化過程是非平穩(wěn)、非線性的，是由多種波動成分共同作用的結果，從中可以提取出4種不同波段的經(jīng)驗模態(tài)分量IMF和一個殘余項：IMF1波動周期不明顯，波動幅度變化比較大；IMF2波動周期與IMF1相比變化明顯，波動周期為準10 a，波動幅度變化??；IMF3波動周期相比于IMF2變大，波動周期為準15 a，波動幅度小于0.5 m；IMF4波動周期最長，波動周期為準27 a，波動幅度在1983年以前變化最小，1983年后變化最大；殘余項表明漏斗中心水位呈現(xiàn)逐年下降趨勢.

對4個經(jīng)驗模態(tài)函數(shù)IMFs進行檢驗，檢驗結果均為平穩(wěn)的非白噪聲序列.由AIC準則確定IMF1,IMF2,IMF3,IMF4的AR模型分別為AR(2),AR(4),AR(6),AR(6).AR模型參數(shù)計算結果見表1.

表1 AR模型參數(shù)計算

采用拉格朗日乘數(shù)檢驗LMT,懷特檢驗WT，分別對殘差序列進行相關性和異方差檢驗[30]，結果表明各經(jīng)驗模態(tài)函數(shù)IMFs的AR模型均通過檢驗.

選取1953-2002年作為模型率定期，2003-2017年作為模型驗證期，模擬預測結果由表2可知，模擬預測結果與實際觀測結果擬合程度較好，相對誤差均小于5%，模型精度高，模擬效果理想，可利用此模型對未來年份漏斗中心水位進行預測.

表2 地下水位模擬預測結果誤差分析

(2)漏斗面積預測模型構建

根據(jù)2001-2017年夏莊漏斗中心水位和漏斗面積數(shù)據(jù)，利用SPSS軟件的OLS回歸分析模塊對中心水位和面積進行分析.通過計算確定了回歸系數(shù)及模型相應參數(shù)，再通過比選，得到了地下水漏斗面積的回歸方程.模型相應參數(shù)見表3.

表3 模型參數(shù)匯總

從表3中看出，二次函數(shù)建立回歸模型時的R2值最大，為0.861，擬合程度最好，因此夏莊漏斗中心水位和漏斗面積回歸模型如(6)式所示：

Y/m2=-3 234.814+114.865(X/m)-0.947(X/m)2,

(6)

式中：Y為漏斗面積，X為中心水位.

3.2.2漏斗中心水位及面積預測

利用EMD-AR模型對夏莊2022-2030年漏斗中心水位進行預測，預測結果如圖6(漏斗中心水位曲線)所示，2022-2030年漏斗中心水位下降約1.5m，水位下降過程分為兩個階段.第一階段為2025年以前，漏斗中心水位約以0.08m/a速度緩慢下降.第二階段為2025年以后，漏斗中心水位約以每年0.25m/a速度下降.

基于水位預測結果，根據(jù)漏斗中心水位和漏斗面積回歸模型對漏斗面積進行預測，預測結果如圖6(漏斗面積柱狀圖)所示，2022-2030年漏斗面積增加約8.5km2，面積增加過程可分為兩個階段.第一階段為2025年以前，漏斗面積約以1.4km2/a速度增加，增長速度較快.第二階段為2025年以后，漏斗約以0.85km2/a速度增加，增長速度較為緩慢.

計算結果表明，在現(xiàn)有條件下，灌區(qū)漏斗中心水位持續(xù)下降，漏斗面積進一步擴大，亟須采取壓采措施控制漏斗的發(fā)展.

需要指出的是，本文使用的EMD-AR方法明顯的優(yōu)點是所需資料相對不多(主要為長系列水位數(shù)據(jù)和漏斗形成后的面積數(shù)據(jù))，目的是利用EMD對時間序列平穩(wěn)化處理后的本征模函數(shù)IMF自身序列實現(xiàn)對未來情況的預測，是根據(jù)數(shù)據(jù)發(fā)展趨勢進行的預測，屬于統(tǒng)計模型的范疇，如果坡度、坡向、地面植被、開采政策等地下水影響因子不發(fā)生突變的話，數(shù)據(jù)系列具有一致性，用此方法預測效果較好.

4 結 論

(1)1976年以前，灌區(qū)地下水流場未發(fā)生較大變化，未出現(xiàn)地下水漏斗.1976年以后，灌區(qū)地下水流場開始發(fā)生改變，地下水漏斗最早出現(xiàn)在董莊，而后由董莊過渡，最終形成以夏莊為中心的地下水降落漏斗.在漏斗形成中，漏斗中心處的水位和漏斗面積在不斷擴大，中心水位與漏斗面積變化呈現(xiàn)高度相關.

(2)2022-2030年夏莊地下水漏斗中心處水位持續(xù)下降.以2025年為分界點，2022-2025年漏斗中心水位下降緩慢，下降速度約為0.08m/a；2025-2030年漏斗中心水位下降迅速，下降速度約為0.25m/a.2022-2030年夏莊地下水漏斗中心水位下降約1.5m.

(3)隨著漏斗中心水位的降低，漏斗面積也不斷擴大.2022-2025年漏斗面積增長較快，增長速度約為1.4km2/a；2025-2030年漏斗面積增長緩慢，增長速度約為0.85km2/a.2022-2030年夏莊漏斗面積增長約8.5km2.