基于EMD-AR模型的黃河下游灌區(qū)地下水漏斗演變預(yù)測

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

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

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

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

1 研究區(qū)概況

人民勝利渠灌區(qū)(113°31′E～114°25′E，35°0′N～35°30′N)位于黃河北岸，總面積為1 486.84 km2.灌區(qū)屬于溫帶大陸性季風(fēng)氣候，四季分明，春季干旱多沙、夏季炎熱多雨、秋季秋高氣爽、冬季寒冷干燥.多年平均氣溫約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)二段為沖積、洪積、風(fēng)積地層.全新統(tǒng)三段為沼澤堆積地層.依據(jù)含水介質(zhì)及類型，區(qū)域調(diào)查深度范圍內(nèi)地下水可劃分為松散巖類孔隙水，半膠結(jié)碎屑巖類孔隙裂隙水兩類.根據(jù)含水層埋深條件可分為淺層含水組、中深層含水組和深層含水組.按5 m的降深計(jì)算又可分為強(qiáng)富水區(qū)、中富水區(qū)、弱富水區(qū).其中灌區(qū)西部關(guān)堤一帶為強(qiáng)富水區(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個(gè)觀測井(由于個(gè)別觀測井經(jīng)緯度數(shù)據(jù)缺失，某些年份觀測井不到24個(gè))1953-2017年逐月地下水埋深觀測數(shù)據(jù)，地下水水位為觀測井地面高程與埋深之差；地下水漏斗面積數(shù)據(jù)獲取首先通過ArcGIS軟件繪制研究區(qū)地下水等水位線圖，然后識(shí)別地下水漏斗的范圍，最后利用空間分析提取數(shù)據(jù)實(shí)現(xiàn).

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

2.2 研究方法

2.2.1基于EMD的AR模型

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

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

(1)

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

步驟2 對所有的經(jīng)驗(yàn)?zāi)B(tài)函數(shù)IMFs分別建立AR模型[22-23].AR模型計(jì)算流程如圖2所示，模型數(shù)學(xué)表達(dá)式為：

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

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

(2)

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

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

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

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

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

(3)

(4)

2.2.4模型的診斷與修正

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

2.2.5OLS回歸模型

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

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

(5)

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

3 結(jié)果與分析

3.1 地下水漏斗演變分析

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

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

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

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

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

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

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

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

3.2 地下水漏斗模型構(gòu)建及預(yù)測

3.2.1漏斗中心水位及面積模型構(gòu)建

(1)漏斗中心水位預(yù)測模型構(gòu)建

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

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

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

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

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

表2 地下水位模擬預(yù)測結(jié)果誤差分析

(2)漏斗面積預(yù)測模型構(gòu)建

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

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

從表3中看出，二次函數(shù)建立回歸模型時(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漏斗中心水位及面積預(yù)測

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

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

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

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

4 結(jié) 論

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

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

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