李仙岳 崔佳琪 史海濱 孫亞楠 邢進平

(1.內(nèi)蒙古農(nóng)業(yè)大學(xué)水利與土木建筑工程學(xué)院， 呼和浩特 010018；2.河套灌區(qū)永濟灌域管理局， 巴彥淖爾 015000)

0 引言

地下水作為人類賴以生存的水資源，特別是在我國西北干旱地區(qū)，在維持區(qū)域社會經(jīng)濟和生態(tài)可持續(xù)發(fā)展方面起到重要作用[1-2]。由于氣候變化、不合理的農(nóng)業(yè)灌溉、排水不暢等因素，導(dǎo)致部分區(qū)域地下水位超過臨界水位，土壤鹽堿化問題突顯[3]，影響了灌區(qū)生態(tài)環(huán)境的健康發(fā)展和農(nóng)業(yè)的可持續(xù)發(fā)展。另外，一些區(qū)域由于大規(guī)模節(jié)水改造工程的實施，導(dǎo)致地下水位下降，從而使作物利用地下水量下降、灌溉水量增加[4]。準確預(yù)測灌區(qū)地下水埋深對于農(nóng)田灌溉研究以及灌區(qū)合理用水都具有重要意義[5]。由于地下水埋深受自然條件和人類活動雙重作用的影響，其變化過程非常復(fù)雜[6]。

目前，對地下水位預(yù)測的研究方法較多，主要包括基于物理模型的數(shù)值模擬和基于數(shù)據(jù)驅(qū)動的機器學(xué)習(xí)算法[7]。其中，數(shù)值模擬主要是基于MODFLOW、FEFLOW等模型進行的各種地下水動態(tài)模擬。如楊洋等[8]應(yīng)用MODFLOW模型對井渠結(jié)合后內(nèi)蒙古河套灌區(qū)地下水動態(tài)變化進行模擬，并進行了驗證；楊廣等[9]基于MODFLOW模型對瑪納斯河流域下游地下水位進行了預(yù)測，發(fā)現(xiàn)通過調(diào)整種植結(jié)構(gòu)可減少地下水開采量，以減少地下水負均衡量；LI等[10]應(yīng)用FEFLOW模型對黑河流域綠洲水文過程進行模擬，結(jié)果表明，在現(xiàn)狀和節(jié)水措施下，未來10年的地下水位將持續(xù)下降。數(shù)值模擬方法對邊界條件、水文地質(zhì)參數(shù)等要求嚴格，參數(shù)較難獲取，而且由于邊界及含水層的設(shè)置不當容易導(dǎo)致預(yù)測精度降低、模擬失真等[11-12]。

地下水預(yù)測的機器學(xué)習(xí)算法主要基于BP神經(jīng)網(wǎng)絡(luò)[13-14]、灰色系統(tǒng)模型[15]、支持向量機[16]、馬爾科夫鏈[17-18]等模型，這些模型各具特點，且不需要各類物理參數(shù)，因此在不同學(xué)科中均得到了廣泛應(yīng)用。但這些預(yù)測模型的建模方法要求時間序列具有平穩(wěn)性、獨立性和正態(tài)性，如果時間序列過于復(fù)雜，則得到的預(yù)測結(jié)果不精確。多變量時間序列CAR(Controlled auto-regressive)模型結(jié)合了回歸分析和一維時間序列分析兩種方法的優(yōu)點，能夠較好地模擬和預(yù)報，避免了因處理復(fù)雜時間序列而引起的誤差。目前已廣泛用于水文領(lǐng)域的徑流、水量及人口增長、氣候變化趨勢等預(yù)測上[19-24]。管孝艷等[25]基于時間序列CAR模型對河套灌區(qū)沙壕渠灌域地下水埋深進行了預(yù)測，預(yù)測效果良好，證明該模型具有較好的適用性。ZENG等[26]基于多變量時間序列建立了洞庭湖地區(qū)地下水資源預(yù)測模型，并對不同方案下的地下水資源量進行了預(yù)測。YAO等[27]基于CAR模型預(yù)測了氣候變化及人類活動對艾比湖流域3條河流年徑流量的影響，其中博爾塔拉河和徑河的年徑流量呈增加趨勢，而奎屯河呈輕微減少趨勢。另外，基于多變量時間序列模型和其他模型的耦合使用也得到了廣泛應(yīng)用。如張展羽等[28]基于主成分分析與多變量時間序列進行耦合，對濟南市陡溝灌區(qū)的地下水位進行預(yù)測，結(jié)果表明，適當引入地表水灌溉和減少地下水開采，灌區(qū)地下水位將逐步回升；MANZIONE等[29]基于時間序列結(jié)合地質(zhì)統(tǒng)計學(xué)構(gòu)建了隨機模型，用于巴西東南部Guarani含水層出露地區(qū)的地下水位預(yù)測，為地下水管理提供了重要決策依據(jù)。

多變量時間序列模型在地下水預(yù)測中得到了廣泛應(yīng)用，但北方灌區(qū)冬季氣溫低，存在凍融期，如河套灌區(qū)有超過5個月的凍融期，凍融期內(nèi)地下水埋深受到氣溫的約束[30]，所以采用年尺度數(shù)據(jù)源構(gòu)建CAR模型勢必會忽略凍融造成的誤差，而采用月尺度數(shù)據(jù)源構(gòu)建CAR模型又難以刻畫引水等因素的滯后效應(yīng)，同時在凍融地區(qū)CAR模型中是否增加氣溫這個輸入變量也會對預(yù)測結(jié)果產(chǎn)生影響。目前針對凍融區(qū)基于不同時間尺度數(shù)據(jù)源CAR模型的地下水預(yù)測研究較少。本文以河套灌區(qū)永濟灌域為研究對象，研究不同時間尺度(月、季、年)數(shù)據(jù)源，并考慮不同影響因子輸入變量，應(yīng)用多變量時間序列CAR模型進行地下水動態(tài)變化預(yù)測，并對其差異性進行分析，優(yōu)選適合凍融灌區(qū)的CAR地下水預(yù)測模型，以期為該凍融灌區(qū)地下水的準確預(yù)測提供理論依據(jù)。

1 材料與方法

1.1 研究區(qū)概況

永濟灌域(107°13′～107°42′E，40°36′～41°13′N)地處河套灌區(qū)中部，蒸發(fā)強烈、干燥少雨、日溫差大、日照時間長，年平均降水量133 mm，年平均蒸發(fā)量2 327 mm，蒸發(fā)量是降水量的10～16倍，年平均氣溫7.5～9.7℃，屬于典型的溫帶大陸性干旱、半干旱氣候。農(nóng)業(yè)用水幾乎全部引于黃河水，1999—2017年引黃水量在7.4～10.5億m3之間。灌域地形西南高，東北低，地下水總的流向自西向東，大部分的淺層地下水靠灌域內(nèi)數(shù)條自南而北的排水干渠匯入北部由西向東的總排水干渠，然后再經(jīng)總排水干渠流入東部的烏梁素海，地下水的運動以垂直交替為主。地下水的主要補給來自引黃的入滲補給，其次為降水入滲補給。蒸發(fā)消耗為灌域主要排泄途徑，地下水循環(huán)屬于典型的入滲-蒸發(fā)類型。11月下旬至翌年4月為灌域的凍融期，凍結(jié)至融通歷時180 d。按水文地質(zhì)特點灌域內(nèi)水文地質(zhì)結(jié)構(gòu)分為第一含水層組和第二含水層組。第二含水層組埋深較大，與外界水量交換較少，且以咸水為主，不宜開采利用，故本文以第一含水層組為研究對象。

1.2 數(shù)據(jù)采集和處理

在永濟灌域內(nèi)共布設(shè)了44眼長期地下水觀測井(圖1)，采用人工法在每月1、5、11、21、26日測地下水位埋深，并在干渠設(shè)流量日監(jiān)測點，數(shù)據(jù)來源于永濟灌域管理局，監(jiān)測時間為1999—2017年。氣象數(shù)據(jù)來源中國氣象科學(xué)數(shù)據(jù)共享網(wǎng)(http:∥data.cma.cn/)，包括1999—2017年逐月平均氣溫、逐月累計降雨量、逐月累計水面蒸發(fā)量。

圖1 永濟灌域地下水采樣點分布Fig.1 Distribution of groundwater sampling points in Yongji irrigation field

2 基于多變量時間序列CAR的地下水預(yù)測模型

2.1 CAR模型建模原理

多變量時間序列CAR模型主要采用遞推最小二乘法進行模型參數(shù)評估的方法建模[31-32]。由m個輸入變量和1個輸出變量組建成n階CAR模型，其形式為

yt=a1yt-1+a2yt-2+…+anyt-n+b10x1,t+b11x1,t-1+…+b1nx1,t-n+b20x2,t+b21x2,t-1+…+b2nx2,t-n+…+bm0xm,t+bm1xm,t-1+…+bmnxm,t-n+εt

(1)

式中 {an}、{bmn}——系數(shù)，其中m為正整數(shù)，n為非負整數(shù)

t——時間序列編號，t>n

yt——輸出變量

xm,t-n——輸入變量

εt——殘差

2.2 CAR模型建模步驟

(1)CAR模型的參數(shù)估計

采用遞推最小二乘法參數(shù)估計

θ=[a1,a2,…,an;b10,b11,…,b1n;b20,b21,…,b2n]Tzt=[yt-1,yt-2,…,yt-n;x1,t-0,x1,t-1,…,x1,t-n;x2,t-0,x2,t-1,…,x2,t-n]T

則CAR模型的一般形式可以寫為

(2)

根據(jù)式(2)，在時刻t，θ的遞推最小二乘估計值為

(3)

式中β——遺忘因子，一般取0.9～1.0

I——單位矩陣

參數(shù)Pt初值為λI，其中λ取104，則可利用N組觀察值得到的CAR(n)模型計算殘差平方和

(4)

(2)CAR模型最高階n的判定

根據(jù)已知的N樣本(xj,t,yt)(t=1,2,…,N;j=1,2,…,m)，由低階到高階遞增地對系統(tǒng)進行擬合，并且依次對相鄰的兩個CAR模型采用F檢驗的方法判斷模型階次增加是否合適，得出合適的CAR模型。

(5)

取置信度α=0.05，大均方自由度為m，小均方自由度為N-mn-(m+1)，求出相應(yīng)的臨界值Fα。若F

(3)CAR模型真實階及時滯的檢驗

為避免以上步驟后某些參數(shù)可能為0，對CAR(n)模型中某些系數(shù)是否為0進行F檢驗，將接近于0的參數(shù)剔除后，重新應(yīng)用遞推最小二乘法建立較少參數(shù)的CAR(n)模型，并繼續(xù)采用F檢驗的方法進行檢驗。若不顯著，則較少參數(shù)的CAR(n)模型是真實模型；若顯著，則原來的CAR(n)模型是真實模型。用以決定模型的真實階及其時滯，得到真實模型的參數(shù)估計。

經(jīng)過上述3個步驟，即可建立只保留對系統(tǒng)影響較大因素的較少參數(shù)CAR(n)模型。

(4)模型的率定和驗證

由于永濟灌域有超過5個月的凍融期，考慮到不同時期氣溫因素對地下水埋深的影響不同，本文以永濟灌域1999—2017年月均、季均(生育期-秋澆期-凍融期)、年均引水量X1(億m3)和蒸發(fā)量X2(mm)、降雨量X3(mm)、氣溫X4(℃)作為輸入變量，地下水埋深Y(mm)作為輸出變量，建立不同時間尺度數(shù)據(jù)源(月、季、年)不考慮氣溫、考慮氣溫和只考慮凍融期(12月到翌年4月)氣溫的地下水埋深預(yù)測模型。建模及因子檢驗的顯著性水平為0.05，建模所用遞推最小二乘法的遺忘因子為1.0。為驗證模型的可行性，以不同時間尺度考慮氣溫因素的CAR(T)模型為例，月尺度數(shù)據(jù)源、季尺度數(shù)據(jù)源和年尺度數(shù)據(jù)源CAR模型分別記為CAR(T)M、CAR(T)Q和CAR(T)Y，其中1999—2013年數(shù)據(jù)用于模型的率定，2014—2017年數(shù)據(jù)用于模型的驗證。

2.3 模型評價指標

采用決定系數(shù)(R2)、均方根誤差(RMSE)和Nash-Sutcliffe系數(shù)(Ens)3種指標[33-34]對各模型模擬效果進行評價。其中R2和Ens越接近1，RMSE越小，說明預(yù)測值與實測值相差越小，擬合精度越高。反之，模型擬合精度較差。同時，通過對實測值與預(yù)測值的擬合誤差進行分析，進一步評價模型的模擬精度。

3 結(jié)果與分析

3.1 地下水埋深影響因素分析

隨著灌區(qū)續(xù)建配套及節(jié)水改造的大規(guī)模實施，總體上灌域地下水埋深呈下降趨勢，1999—2017年期間，永濟灌域年平均地下水埋深為2.16 m。按照灌域的灌水特征，通常在每年4月中旬開閘引水，11月中旬停灌，其中大規(guī)模秋澆時間為10月和11月，而在12月至翌年4月河套灌區(qū)土壤和地下水處于封凍-消融階段。為了便于季尺度地下水預(yù)測，將全年分為3個階段，其中12月—翌年4月為凍融期，5—9月為生育期，10—11月為秋澆期。由圖2可知，年內(nèi)地下水埋深的變動主要受引水及生育期耗水的影響呈季節(jié)性周期變化，土壤于11月中旬開始凍結(jié)，凍層厚1～1.5 m，隨著凍層的逐漸加厚，地下水埋深隨之下降，土壤蒸發(fā)降至最小，至翌年3月出現(xiàn)全年最大埋深。此期埋深的下降是由于凍層上下溫差較大，地下水在溫差作用下，向凍層運移，使凍層不斷增厚，使地下水埋深加大。翌年3月無降水與灌溉入滲補給，在凍結(jié)影響下為全年最大埋深2.68 m，土壤消融始于3月初，于4月末完全解凍，凍層以下融凍水回補地下水，埋深開始上升。凍融期年均地下水埋深為2.42 m，其主要影響因素是土壤溫度，而土壤溫度受外界氣溫影響[35-36]。5月中旬夏灌開始，地下水位大幅度上升，同時蒸發(fā)作用劇烈，地下水消耗于蒸發(fā)，埋深出現(xiàn)峰谷交替變化，生育期間年均地下水埋深波谷出現(xiàn)在6月，為1.78 m，波峰出現(xiàn)在9月，為2.49 m，月均埋深為2.06 m。10月開始大規(guī)模秋澆，灌水量大，時間短，這一階段地下水埋深變淺，由于存在一定的滯后性，至11月出現(xiàn)全年最小埋深1.57 m。非凍融期(生育期和秋澆期)，入滲與蒸發(fā)之間的平衡關(guān)系是地下水埋深變化的主導(dǎo)因素[37]。

圖2 永濟灌域地下水埋深及其影響因素Fig.2 Groundwater depth and its influencing factors in Yongji irrigation field

通過對近20年(1999—2017年)灌域地下水埋深與引黃水量、蒸發(fā)量、降雨量以及氣溫進行相關(guān)分析(表1)，結(jié)果顯示這幾項環(huán)境因子與地下水埋深均有較好的相關(guān)性，其中引黃水量和蒸發(fā)量均達到極顯著相關(guān)，可見該區(qū)域水分入滲和耗水是其地下水變動的主導(dǎo)因子，而且降雨量和氣溫也與地下水埋深呈顯著相關(guān)。另外總體上不同尺度(月、季、年)地下水埋深與環(huán)境因子的相關(guān)性比較相近，不同時間尺度數(shù)據(jù)源相關(guān)性由大到小總體上依次為：年度、季度、月度，這可能是由于年尺度通過數(shù)據(jù)平均弱化了變異點的影響。

表1 地下水埋深相關(guān)因子分析Tab.1 Analysis of correlation factors of groundwater depth

3.2 不同時間尺度數(shù)據(jù)源地下水CAR模型的構(gòu)建與驗證

圖3 不同時間尺度地下水CAR(T)模型率定及驗證Fig.3 Calibration and verification of groundwater CAR(T) model at different time scales

由圖3可知，率定期作為率定模型的數(shù)據(jù)來源時期，其擬合效果最好，驗證期為驗證模型的準確性，其擬合效果略差，以季尺度數(shù)據(jù)源考慮氣溫的CAR(T)Q模型為例，率定期R2、Ens和RMSE分別為0.902、0.893和0.045 m，在驗證期模型各指標精度分別下降了9.87%、11.49%和增加了25.00%，但預(yù)測曲線均很好地跟蹤了觀測曲線的變化趨勢，3種時間尺度數(shù)據(jù)源模型預(yù)測精度均較好，模擬效果較差的月尺度數(shù)據(jù)源CAR(T)M模型的率定期R2和Ens均不小于0.771、RMSE為0.069 m，驗證期預(yù)測結(jié)果相對差些，Ens為0.695。但是月度地下水埋深數(shù)據(jù)量多，能刻畫年內(nèi)的地下水波動，季度也能大體反映出地下水的波動特征，而年度則不能反映地下水年內(nèi)的變化?？傮w來說3種尺度的地下水埋深模型預(yù)測性能都較好，模型較為穩(wěn)定，可進一步應(yīng)用于地下水埋深的預(yù)測研究。

3.3 不同時間尺度數(shù)據(jù)源地下水CAR模型的差異性分析

圖4 不同時間尺度CAR(T)模型差異性分析Fig.4 Difference analysis of CAR(T) models at different time scales

為了探索不同時間尺度數(shù)據(jù)源地下水CAR模型的差異，并篩選出最優(yōu)模型，將1999—2017年期間各因子數(shù)據(jù)按照月、季、年3種時間尺度建立CAR模型，通過對比不同時間尺度考慮氣溫的CAR(T)模型之間的精度，得出最優(yōu)時間尺度模型。由圖4可知，月尺度CAR(T)M模型的相對誤差在-13%～11%之間，季尺度CAR(T)Q模型的相對誤差在-4.8%～8.8%之間，年尺度CAR(T)Y模型的相對誤差在-3.1%～4.0%之間?；?種模型的擬合精度比較可知，擬合效果最好的CAR(T)Q模型的R2、Ens和RMSE較擬合效果較差的CAR(T)M模型分別提高了11.30%、11.86%和降低了32.35%。由圖4d可知，灌域在月尺度的相對誤差明顯高于年尺度和季尺度，部分月份相對誤差超過10%。即CAR(T)Q模型的擬合精度明顯高于CAR(T)M和CAR(T)Y模型，CAR(T)Y模型次之，CAR(T)M相對最差。可能是因為地下水埋深變化在時間序列上呈現(xiàn)滯后性，比如氣溫因素，有研究表明氣溫的變化需要滯后46.5 d才能經(jīng)土體傳至潛水面[38]，此外，部分月份降雨量和引黃水量為零，控制因素減少，擬合結(jié)果出現(xiàn)誤差，這也是基于月尺度的CAR模型擬合結(jié)果相對較差的原因。從而驗證了應(yīng)用季尺度進行地下水埋深動態(tài)變化預(yù)測擬合精度最高。

3.4 氣溫對凍融區(qū)地下水CAR模型的影響

為了探索氣溫對凍融灌區(qū)地下水埋深的影響，將不考慮氣溫CARQ、考慮氣溫CAR(T)Q和只考慮凍融期(12月—翌年4月)氣溫CAR(TF)Q的季尺度模型進行精度比較。由表2可知，CAR(TF)Q模型的R2、Ens、RMSE分別為0.941、0.940和0.044 m，CAR(T)Q模型和CARQ模型的R2和Ens分別降低了0.53%、0.64%和2.98%、3.09%，RMSE分別提高了4.55%和11.36%，可以看出，CAR(TF)Q模型略優(yōu)于CAR(T)Q模型，明顯優(yōu)于CARQ模型。從而驗證了僅考慮凍融期氣溫的CAR(TF)Q模型為最優(yōu)模型。其主要由于凍融期氣溫對地下水埋深有很強的約束作用，加之引黃水量和降雨量對其約束作用減弱，所以加入氣溫對凍融期模型的擬合精度有所提高，而非凍融期氣溫對其約束力不強，主要影響因素為引黃水量、蒸發(fā)量和降雨量。可知此模型在引入數(shù)據(jù)時應(yīng)注意其物理背景，而不是相關(guān)影響因子考慮的越多越好。

表2 季尺度數(shù)據(jù)源條件下不同氣溫因子CAR模型比較Tab.2 Comparison of different temperature conditions CAR model using quarter data

4 結(jié)論

(1)灌域地下水埋深與引黃水量、蒸發(fā)量、降雨量和氣溫呈顯著相關(guān)性，不同時間尺度數(shù)據(jù)源的相關(guān)性由大到小依次為：年度、季度、月度。

(2)季尺度數(shù)據(jù)源CAR模型擬合效果明顯優(yōu)于月尺度數(shù)據(jù)源CAR模型和年尺度數(shù)據(jù)源CAR模型，年尺度數(shù)據(jù)源CAR模型擬合精度次之，月尺度數(shù)據(jù)源CAR模型相對較差。季尺度數(shù)據(jù)源CAR模型的R2、Ens和RMSE較月尺度數(shù)據(jù)源CAR模型分別提高了11.30%、11.86%和降低了32.35%。

(3)僅考慮凍融期氣溫的CAR模型優(yōu)于考慮氣溫的CAR模型，不考慮氣溫的CAR模型擬合結(jié)果相對較差。凍融灌區(qū)最優(yōu)地下水預(yù)測模型為僅考慮凍融期氣溫的季尺度CAR模型，模擬值與實測值的R2為0.941，Ens為0.940，RMSE為0.044 m。