張 婧，馬貴宏，高 雅，沈彥俊，劉曉英，鄒 妍，戴茂華

(1.河北省氣象與生態(tài)環(huán)境實(shí)驗(yàn)室，河北 石家莊 050021；2.河北省氣候中心，河北 石家莊 050021； 3.河北省石家莊市高邑縣氣象局，河北 石家莊 051330；4.河北省水文勘測(cè)研究中心，河北 石家莊 050031； 5.中國(guó)科學(xué)院遺傳與發(fā)育生物學(xué)研究所農(nóng)業(yè)資源研究中心，河北 石家莊 050021；6.天津長(zhǎng)城濱銀汽車金融有限公司，天津 300457； 7.天津師中師教育信息咨詢有限公司，天津 300074；8.河北省農(nóng)林科學(xué)院旱作農(nóng)業(yè)研究所，河北 石家莊 053000)

地下水是世界上最大的淡水資源庫，為全球提供了36%的飲用水和42%的灌溉用水，是農(nóng)業(yè)、工業(yè)、環(huán)境和飲用水必不可少的供水來源[1-2]。華北平原是我國(guó)重要的糧食高產(chǎn)區(qū)之一，由于半干旱半濕潤(rùn)的氣候條件無法滿足兩熟制(冬小麥和夏玉米)作物的需水量，且地表水資源較匱乏，農(nóng)業(yè)地下水灌溉量占地下水總開采量的80%以上，人類活動(dòng)已成為影響華北平原地下水環(huán)境變化的主導(dǎo)因素[3-4]。由于地下水資源的過度開采，1970—2008年華北平原地下水水位累計(jì)下降了8.74 m，山前平原地區(qū)地下水水位降幅最大，達(dá)到16.34 m，地下水漏斗連區(qū)成片引發(fā)了一系列地下水資源和環(huán)境問題[5-9]。隨著華北地區(qū)社會(huì)經(jīng)濟(jì)的快速發(fā)展和城市人口的不斷增長(zhǎng)，城市工業(yè)和生活用水量將繼續(xù)增加，水資源短缺問題將進(jìn)一步加劇，會(huì)嚴(yán)重影響該地區(qū)社會(huì)經(jīng)濟(jì)和生態(tài)環(huán)境的可持續(xù)發(fā)展[4-7]。

在前人對(duì)華北平原地下水水位變化的歸因分析中，降水量減少、氣溫升高、上游水利工程建設(shè)、農(nóng)業(yè)種植結(jié)構(gòu)改變和地下水開采等因素是造成該地區(qū)水資源減少的重要原因[10-11]。其中，人為開采尤其是農(nóng)業(yè)灌溉是地下水水位下降的主要原因，灌溉引起的實(shí)際蒸散量高于降水量，導(dǎo)致該地區(qū)地下水資源持續(xù)減少[12-15]。華北山前平原小麥和玉米種植面積與農(nóng)業(yè)地下水開采量關(guān)系密切，二者的增減變化趨勢(shì)較為一致[14]，糧食產(chǎn)量增產(chǎn)的同時(shí)，作物消耗大量土壤水，造成補(bǔ)給地下水的入滲量大大減少，也是地下水水位連年下降的一個(gè)重要原因[2,15]。

除了人為影響因素外，氣候變化尤其是降水量變化對(duì)區(qū)域水資源影響較為顯著[13,16-19]。在影響華北山前平原含水層的氣候因素中，除了降水量外，還有來自上游的地下水側(cè)向補(bǔ)給[16-18]。山西作為華北平原地下水產(chǎn)流區(qū)和地下水補(bǔ)給區(qū)，自20世紀(jì)末以來降水量明顯減少，導(dǎo)致地表產(chǎn)流和地下水補(bǔ)給量減少，成為太行山山前平原地下水水位下降的另一個(gè)主要原因[17]。降水量和氣溫的波動(dòng)以及蒸散效應(yīng)與地下水的補(bǔ)給和排泄類型具有一定的聯(lián)系，華北平原地下水水位有2～7 a準(zhǔn)周期變化，這與ENSO的波動(dòng)周期較為一致[13,18-19]。一些學(xué)者從氣候要素和人類活動(dòng)兩個(gè)方面開展了綜合分析，認(rèn)為氣候變化影響人類活動(dòng)，進(jìn)而影響到地下水儲(chǔ)量的變化。Taylor等[20]認(rèn)為全球氣候變化對(duì)地下水的間接影響大于其對(duì)地下水補(bǔ)給的直接影響，間接影響主要體現(xiàn)在氣候變化通過影響灌溉需水量和灌溉來源間接影響地下水儲(chǔ)量的變化。石家莊平原降水量影響農(nóng)業(yè)地下水開采強(qiáng)度，干旱年份的地下水開采量明顯多于濕潤(rùn)年份，高強(qiáng)度的地下水開采常常會(huì)掩蓋氣候變化在地下水動(dòng)態(tài)變化中的作用[14,16,21]。因此，華北平原淺層地下水動(dòng)態(tài)變化是氣候與人類活動(dòng)疊加效應(yīng)的結(jié)果。

由于研究區(qū)和分析方法不同，得到的氣候變化和人類活動(dòng)對(duì)地下水動(dòng)態(tài)變化影響的結(jié)論也存在差異。開展地下水動(dòng)態(tài)影響定量研究的方法較多，主要有水文模型和數(shù)學(xué)統(tǒng)計(jì)方法兩大類。利用MODFLOW、VFM等水文模型可以模擬地下水水位動(dòng)態(tài)變化，而主成分分析、灰色理論等統(tǒng)計(jì)方法能較客觀地分析地下水動(dòng)態(tài)變化中各影響因子的響應(yīng)。許月卿[22]利用投影尋蹤回歸技術(shù)，從降水、氣溫等自然因素和地下水超采、水利工程、農(nóng)田產(chǎn)量、作物結(jié)構(gòu)等人為因素兩方面得出河北平原地下水開采是引起地下水水位下降的最主要因素，相對(duì)貢獻(xiàn)率為54.7%，其次是降水和地表徑流，相對(duì)貢獻(xiàn)率為25.6%和19.7%。張展羽等[23]基于主成分分析與多變量時(shí)間序列(controlled auto-regressive，CAR)構(gòu)建了濟(jì)南市陡溝灌區(qū)地下水水位模型，通過對(duì)比主成分回歸模型、逐步回歸模型和BP神經(jīng)網(wǎng)絡(luò)模型的模擬結(jié)果，發(fā)現(xiàn)主成分-時(shí)間序列模型模擬結(jié)果與實(shí)測(cè)值相對(duì)誤差最小，模擬效果較好。

在已有研究中，Wang等[18]利用同位素示蹤法發(fā)現(xiàn)華北山前平原的元氏地區(qū)上游地下水側(cè)向補(bǔ)給元氏的地下水量占地下水總補(bǔ)給量的50%～88%，雖然上游的側(cè)向補(bǔ)給對(duì)華北山前平原地下水影響較大，但目前考慮上游地區(qū)對(duì)山前平原地下水動(dòng)態(tài)變化的定量影響研究尚未見報(bào)道，且華北山前平原淺層地下水埋深較大，地質(zhì)條件復(fù)雜，目前對(duì)該地區(qū)地下水水位變化原因的定量分析較少。本文選取華北山前平原典型井灌區(qū)石家莊欒城為研究區(qū)，從氣候條件和人類活動(dòng)兩方面對(duì)欒城地下水位變化的原因進(jìn)行分析，以期為合理制定“以水定產(chǎn)”的地下水管理措施和促進(jìn)高效農(nóng)業(yè)節(jié)水灌溉提供參考。

1 研究數(shù)據(jù)和方法

1.1 研究區(qū)概況

華北山前平原位于華北平原的西部和北部的狹長(zhǎng)地帶，寬度一般在30～60 km，土層深厚，有機(jī)質(zhì)含量較高，地下水自西北流向東南。受季風(fēng)和地形等因素影響，山前平原年降水波動(dòng)較大，年內(nèi)降水時(shí)間分布不均，6—9月降水量約占全年的70%～80%[13]。氣候變化背景下極端氣候事件增多以及人類活動(dòng)強(qiáng)度日趨增大等因素使該地區(qū)地下水水位變化原因復(fù)雜，且各影響因素間存在一定的相互作用。

石家莊欒城區(qū)為華北山前平原石家莊市東南部的典型井灌區(qū)，該地區(qū)地勢(shì)平坦，以冬小麥和夏玉米輪作制農(nóng)田為主，農(nóng)田有效灌溉面積比例達(dá)到100%，占區(qū)域總面積的68%。欒城區(qū)冬小麥生育期自10月初至次年6月初，夏玉米生育期自冬小麥?zhǔn)斋@至9月末，農(nóng)業(yè)用水占全部用水量的88%以上。

研究區(qū)的地下水位觀測(cè)點(diǎn)位于中國(guó)科學(xué)院欒城農(nóng)業(yè)生態(tài)系統(tǒng)試驗(yàn)站內(nèi)(37.53°N，114.41°E，海拔為50.1 m)，該試驗(yàn)站自1974年開始記錄地下水水位數(shù)據(jù)，其地下水水位波動(dòng)趨勢(shì)與華北山前平原區(qū)地下水水位變化具有較好的一致性[11,24]，能較好地代表欒城冬小麥和夏玉米輪作系統(tǒng)下的地下水水位變化。

1.2 研究數(shù)據(jù)

1975—2015年地下水位埋深觀測(cè)數(shù)據(jù)來自中國(guó)科學(xué)院欒城農(nóng)業(yè)生態(tài)系統(tǒng)試驗(yàn)站，該數(shù)據(jù)采用加拿大Solinst公司制造的地下水水位尺(-101型)設(shè)備觀測(cè)，儀器精度為±1 mm。農(nóng)田實(shí)際蒸散月數(shù)據(jù)(2008—2015年)由渦度相關(guān)系統(tǒng)觀測(cè)并計(jì)算得到[25-26]。氣象數(shù)據(jù)(1994—2015年)來自河北省氣象局氣象信息中心，包括日降水量、日平均氣溫、日最高氣溫、日最低氣溫、日平均風(fēng)速、日平均相對(duì)濕度和日照數(shù)據(jù)等。社會(huì)經(jīng)濟(jì)數(shù)據(jù)(1994—2015年)來自《河北農(nóng)村統(tǒng)計(jì)年鑒》，包括石家莊市區(qū)人口數(shù)和欒城、正定、平山、井陘4縣(區(qū))的小麥和玉米的種植面積、有效灌溉面積和農(nóng)村用電量。

1.3 研究方法

1.3.1 水分盈虧

本文將年降水量減去蒸散量定義為水分盈虧量：

D=P-ET

(1)

式中：D為水分盈虧量，mm；P為年降水量，mm；ET為蒸散量，mm。當(dāng)水分盈虧量為正值時(shí)，表明該地區(qū)水分有盈余；當(dāng)水分盈虧量為負(fù)值時(shí)，表明年降水量不足以彌補(bǔ)作物生長(zhǎng)的蒸散量，此時(shí)出現(xiàn)水量虧缺。當(dāng)蒸散量為渦度相關(guān)系統(tǒng)觀測(cè)值時(shí)，ET為實(shí)際蒸散量，本文利用實(shí)際蒸散量等數(shù)據(jù)分析研究區(qū)的年內(nèi)水分盈虧變化。

由于研究區(qū)實(shí)際蒸散量數(shù)據(jù)的時(shí)間序列較短，在分析1994—2015年地下水水位變化影響因素中的水分盈虧量時(shí)，選用1998年聯(lián)合國(guó)糧農(nóng)組織(FAO)推薦并修訂的Penman-Monteith模型計(jì)算參考作物的潛在蒸散量：

(2)

式中：ET0為參考作物的潛在蒸散量，mm/d；Δ為飽和蒸汽壓與空氣溫度曲線的斜率，kPa/℃；Rn為參考作物凈輻射，MJ/(m2·d)；G為土壤熱通量，MJ/(m2·d)；γ為干濕表常數(shù)；T為2 m高處日平均氣溫，℃；u2為2 m高處日平均風(fēng)速，m/s；es為飽和水汽壓，kPa；ea為實(shí)際水汽壓，kPa。

1.3.2 數(shù)據(jù)標(biāo)準(zhǔn)化

由于不同類型數(shù)據(jù)的單位和數(shù)量級(jí)不同，在進(jìn)行主成分分析之前，需要先將原始數(shù)據(jù)標(biāo)準(zhǔn)化。z-score是常用的標(biāo)準(zhǔn)化方法之一，該方法基于原始數(shù)據(jù)的均值和標(biāo)準(zhǔn)差對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，標(biāo)準(zhǔn)化后為無量綱數(shù)據(jù)，均值為0，方差為1，計(jì)算公式為

(3)

式中：x為原始數(shù)據(jù)；μ為全部數(shù)據(jù)的均值；σ為標(biāo)準(zhǔn)差。

1.3.3 主成分回歸分析法

華北山前平原地下水水位變化原因復(fù)雜，在進(jìn)行變量分析時(shí)，為了消除多個(gè)影響因子間的相互作用，本文采用主成分回歸分析法進(jìn)行擬合。主成分回歸分析是多元統(tǒng)計(jì)分析中減少自變量維數(shù)的一種較理想的工具[27]，因此被廣泛應(yīng)用于氣候變化和水資源等領(lǐng)域。

主成分回歸分析的基本思路是將原來的指標(biāo)體系重新組合成一組(若干個(gè))新的且相互無關(guān)的綜合指標(biāo)，從中選取自變量所有線性組合中方差最大者作為主成分，并盡可能多地反映原來指標(biāo)體系的信息，最后以主成分代替原變量進(jìn)行回歸分析[28-29]。在主成分基礎(chǔ)上構(gòu)建了多元線性回歸模型，將回歸方程中各主成分系數(shù)的權(quán)重作為各主要影響因子對(duì)地下水水位變化的貢獻(xiàn)率。

圖1 1975—2015年觀測(cè)點(diǎn)地下水水位年內(nèi)和年際變化Fig.1 Annual and interannual change of groundwater level in the observation point during 1975-2015

圖2 2008—2015年觀測(cè)點(diǎn)水平衡項(xiàng)和地下水水位月變化Fig.2 Monthly changes of water balance and groundwater level in observation point during 2008-2015

1.3.4 多變量時(shí)間序列分析法

多變量時(shí)間序列分析法[30]可從多變量的時(shí)間序列中提取有效信息表征復(fù)雜系統(tǒng)的動(dòng)態(tài)特征，廣泛應(yīng)用于水文預(yù)測(cè)，具有較好的預(yù)測(cè)效果。在主成分回歸分析基礎(chǔ)上，將主成分回歸分析構(gòu)造的綜合影響因子作為輸入變量，導(dǎo)入多變量時(shí)間序列CAR模型[31]，采用遞推最小二乘法對(duì)估計(jì)模型參數(shù)由低階向高階進(jìn)行多變量時(shí)間序列模擬，通過擬合優(yōu)度確定模型的最高階次。為了檢驗(yàn)所構(gòu)建的多變量時(shí)間序列模型，本文從1994—2015年共21個(gè)樣本中選取前16個(gè)構(gòu)建模型，用剩余樣本進(jìn)行檢驗(yàn)。

對(duì)m個(gè)變量的時(shí)間序列組構(gòu)建n階多變量時(shí)間序列模型，計(jì)算公式為

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

(4)

式中：aj、bij為多變量時(shí)間序列模型系數(shù)(i=1，2，…，m；j=1，2，…，n)；yt-j、xi,t-j為因變量和自變量的時(shí)間序列；t為時(shí)間序列，t>1。

1.3.5 模型檢驗(yàn)

模型驗(yàn)證過程中采用決定系數(shù)(R2)、模擬值與實(shí)測(cè)值的均方根誤差(RMSE)、相對(duì)誤差(MRE)和納什系數(shù)4個(gè)評(píng)價(jià)指標(biāo)來表征模型的模擬效果。

2 結(jié)果與分析

2.1 地下水水位年變化

圖1為1975—2015年觀測(cè)點(diǎn)地下水水位年內(nèi)和年際變化，1975—2015年欒城年平均地下水水位呈顯著下降趨勢(shì)(P<0.001)，平均每年下降0.84 m，年平均地下水水位變化幅度在-2.8～2.9 m之間。1975—2015年欒城年平均地下水水位由37.4 m下降至5.9 m，同時(shí)，地下水埋深不斷增加，由1975年的12.7 m增加至2015年的44.2 m。期間僅有5年(1976年、1977年、1982年、1996年和2003年)年平均地下水水位回升，其中1996年水位上升幅度最大，回升近3 m。1997年和1998年連續(xù)兩年降水量偏少使地下水水位下降幅度均超過2.5 m，2004—2015年地下水水位出現(xiàn)連續(xù)11年持續(xù)下降現(xiàn)象，地下水資源持續(xù)減少。

2.2 影響地下水水位變化的氣象因素

受季風(fēng)氣候影響，2008—2015年欒城年內(nèi)降水量集中在6—9月，占全年降水量的77.9%，其他時(shí)段月降水量大多不足30 mm。由圖2可以看出，欒城典型井灌區(qū)實(shí)際蒸散觀測(cè)值存在兩個(gè)峰值，與一年兩熟的作物需水變化相對(duì)應(yīng)，分別是冬小麥需水量最大的5月和夏玉米需水量最大的8月。夏玉米生育期與當(dāng)?shù)貧夂蛴隉嵬荆医邓扛哂趯?shí)際蒸散量，因此在8—9月地下水水位有所回升。冬小麥生育期(3—5月)降水量較少，不能滿足作物生長(zhǎng)需要，地下水高強(qiáng)度灌溉導(dǎo)致地下水水位在該期間持續(xù)下降，多年平均地下水水位累積下降2.9 m左右。由此可見，欒城降水量和水分盈虧是影響其地下水動(dòng)態(tài)變化的主要?dú)庀笠蛩兀?008—2015年欒城區(qū)高強(qiáng)度地下水灌溉是保障該地區(qū)農(nóng)作物穩(wěn)產(chǎn)、高產(chǎn)的關(guān)鍵[15]。

2.3 影響地下水水位變化的人為因素

2.3.1 當(dāng)?shù)剞r(nóng)業(yè)生產(chǎn)

華北平原地下水水位年際變化主要受人為因素影響，農(nóng)業(yè)灌溉是引起當(dāng)?shù)氐叵滤幌陆档闹饕绊懸蛩亍?967年以來，華北平原冬小麥和夏玉米種植面積明顯增加，地下水水位開始逐年下降，種植結(jié)構(gòu)的變化對(duì)欒城地下水動(dòng)態(tài)變化影響顯著，地下水采補(bǔ)失衡嚴(yán)重[32]。20世紀(jì)50年代至21世紀(jì)初，石家莊平原農(nóng)業(yè)地下水開采量與小麥和玉米總種植面積存在密切相關(guān)性，二者變化趨勢(shì)較為一致，年際和年代際地下水農(nóng)業(yè)開采量隨小麥和玉米總種植面積的變化而變化[11]。因此，本文將欒城小麥和玉米的種植面積作為當(dāng)?shù)厝祟惢顒?dòng)對(duì)地下水動(dòng)態(tài)變化的影響因子。

2.3.2 上游人為因素

欒城農(nóng)田井灌區(qū)在雨季(6—9月)后至次年3月冬小麥春灌前，地下水水位以持續(xù)回升為主，除雨季降水直接補(bǔ)給外，上游地區(qū)地下水側(cè)向補(bǔ)給量占總補(bǔ)給量的50%以上[18]。20世紀(jì)70年代末以來，欒城區(qū)境內(nèi)交河、沙河和泥河上游先后修建中小型水庫30余座，使該地區(qū)4條季節(jié)性河流的年均徑流量明顯減少，沙河和泥河幾乎無水下泄，河流補(bǔ)給地下水量也顯著減少[11]。根據(jù)《華北平原地下水可持續(xù)利用地圖集》[33]中華北平原淺層地下水等水位線圖，欒城及周邊地下水水位自西北向東南方向遞減，假設(shè)上游距離欒城越近，對(duì)欒城地下水補(bǔ)給影響越大。通過對(duì)比欒城區(qū)西北方向的石家莊市區(qū)和鄰近3個(gè)縣(正定、平山和井陘)的小麥和玉米種植面積、有效灌溉面積、農(nóng)村用電量和人口數(shù)量等因子在不同組合下構(gòu)建的主成分回歸模型的顯著性差異，發(fā)現(xiàn)石家莊市區(qū)人口作為上游人為因素時(shí)模擬效果最佳，這主要是由于石家莊市區(qū)為欒城區(qū)地下水上游最鄰近地區(qū)，近些年石家莊市區(qū)城市化進(jìn)程加快，市區(qū)水資源消耗和地下水開采量顯著增加，影響了對(duì)下游欒城區(qū)的地下水側(cè)向補(bǔ)給，因此，本文將石家莊市區(qū)人口作為欒城上游的主要人為因素。

2.4 氣候和人為因素對(duì)地下水動(dòng)態(tài)變化影響的定量分析

表1 1994—2015年潛在自變量的統(tǒng)計(jì)特征量

本文在分析欒城及其上游氣候因子和人類活動(dòng)對(duì)欒城地下水水位動(dòng)態(tài)變化的定性影響基礎(chǔ)上，考慮了欒城及其上游(正定、平山、井陘和石家莊市區(qū))的降水量、水分盈虧量、小麥和玉米種植面積、農(nóng)村用電量、人口數(shù)量等多個(gè)影響因子，通過對(duì)比不同組合下的擬合模型的模擬結(jié)果，最終將欒城地下水水位作為因變量，選取欒城降水量(x1)、欒城水分盈虧量(x2)、欒城小麥和玉米種植面積(x3)、欒城上游降水量(x4)和石家莊市區(qū)人口(x5)作為潛在自變量，潛在自變量在1994—2015年的統(tǒng)計(jì)特征量見表1。

表2 各主成分的特征值、貢獻(xiàn)率和累積貢獻(xiàn)率

從表2可以看出，最大的特征值為1.670，其貢獻(xiàn)率為55.8%，前3個(gè)主成分的累積貢獻(xiàn)率高達(dá)97.0%，可以代表原指標(biāo)97.0%的信息，用這3個(gè)主成分即可代替原指標(biāo)，而且他們是相互獨(dú)立的影響因子。用表3中前3個(gè)主成分的特征向量為系數(shù)構(gòu)造各主成分的線性方程：

Z1=-0.591x1-0.584x2-0.551x4

(5)

Z2=-0.781x3-0.194x4-0.585x5

(6)

Z3=-0.557x3-0.180x4+0.810x5

(7)

從主成分的線性方程來看，式(5)中欒城降水量(x1)的

表3 各主成分特征值對(duì)應(yīng)的特征向量

系數(shù)最大，對(duì)Z1的波動(dòng)起主導(dǎo)作用，所以將第一主成分Z1視作本地降水因子；式(6)中欒城小麥和玉米的種植面積(x3)的系數(shù)最大，對(duì)Z2的變化起主導(dǎo)作用，將第二主成分Z2視作本地人為因子；式(7)中上游石家莊市區(qū)人口(x5)的系數(shù)最大，對(duì)Z3的變化起主導(dǎo)作用，可視作上游主要影響因子。這3個(gè)因子彼此之間相互獨(dú)立，不存在共線性，且代表了原指標(biāo)的大部分信息。以本地降水因子Z1、本地人為因子Z2和上游影響因子Z3為自變量，以地下水水位為因變量，利用最小二乘法進(jìn)行多元線性擬合后，得到的回歸方程為

Y=0.112 6Z1-0.490 5Z2-0.360 6Z3

(8)

回歸方程的相關(guān)系數(shù)R為0.62，決定系數(shù)R2為0.39，F(xiàn)檢驗(yàn)值為3.79，通過了0.01的顯著性檢驗(yàn)，說明回歸方程擬合效果較理想。從回歸方程中3個(gè)因子系數(shù)的權(quán)重得出本地降水因子對(duì)欒城地下水水位回升的貢獻(xiàn)率為11.7%，本地人為因子貢獻(xiàn)率為-50.9%，上游影 響因子貢獻(xiàn)率為-37.4%(正負(fù)號(hào)表示對(duì)地下水水位回升的正負(fù)作用)。由此可以看出，欒城小麥和玉米種植面積對(duì)地下水水位變化的影響最大，人為因子在地下水水位的變化中起主導(dǎo)作用。與人為因子相比，降水因子對(duì)地下水水位影響較小，地下水水位隨降水量的增加出現(xiàn)不同程度的回升?？傮w上，在以欒城典型井灌區(qū)為代表的太行山山前平原區(qū)，其地下水水位變化主要來自人類活動(dòng)的影響，受自然因素影響的比重僅占11.7%。

2.5 主成分-時(shí)間序列模型

地下水水位建模對(duì)流域水資源管理具有重要意義和必要性，也是井灌區(qū)用水管理的重要依據(jù)[34]。束龍倉[35]在分析石家莊地下水水位影響因素和構(gòu)建回歸模型時(shí)發(fā)現(xiàn)本月地下水水位與本月降水量、本月開采量和前一個(gè)月水位關(guān)系較密切。本文借鑒隨機(jī)理論和成因理論的思想，在構(gòu)建欒城地下水水位模型時(shí)，基于上文主成分回歸分析所選出的3個(gè)主要綜合因子，考慮上一年地下水水位對(duì)本年度的影響，根據(jù)決定系數(shù)將模型定階到上一年(t-1)，剔除不顯著因子后，將上一年地下水水位、當(dāng)年的Z1、上一年的Z2和上一年的Z1作為自變量，以當(dāng)年地下水水位為因變量，利用1994—2010年數(shù)據(jù)構(gòu)建研究區(qū)地下水水位主成分-時(shí)間序列模型：

Y=-0.185 13+1.178 46Yt-1-0.024Z1，t+0.085 6Z2，t-1-0.081 5Z1，t-1

(9)

該主成分-時(shí)間序列模型的模擬值與實(shí)測(cè)值的相關(guān)系數(shù)R為0.994，決定系數(shù)R2為0.988，平均相對(duì)誤差為0.08，均方根誤差為0.132，納什系數(shù)為0.98。

圖3 2011—2015年觀測(cè)點(diǎn)地下水水位模擬值和觀測(cè)值Fig.3 Comparison between simulated and measured groundwater level in test station of Luancheng during 2011—2015

為檢驗(yàn)上述模型的效果，對(duì)比2011—2015年欒城典型井灌區(qū)模型的模擬值與地下水水位觀測(cè)值(圖3)發(fā)現(xiàn)，驗(yàn)證期的模擬值與觀測(cè)值變化趨勢(shì)一致，相對(duì)誤差在-19.4%～2.1%之間，說明該模型適用于欒城井灌區(qū)的地下水水位預(yù)測(cè)，即利用降水量、水分盈虧量、種植規(guī)模等數(shù)據(jù)可以較好地推測(cè)當(dāng)年的地下水水位。

2.6 討論

華北平原人口眾多，工農(nóng)業(yè)用水量大，地下水是該地區(qū)最重要的供水水源和戰(zhàn)略資源。在氣候變化和人類活動(dòng)共同影響下，華北山前平原土地利用、種植結(jié)構(gòu)和城鎮(zhèn)化的發(fā)展和變化顯著增加了地下水開采量，地下水水位持續(xù)下降。雖然欒城在1994—2015年小麥和玉米種植面積有所減少，但糧食產(chǎn)量和單產(chǎn)總體增加使作物實(shí)際耗水量增加，作物生育期的降水量不能滿足作物需水量，即使在雨季地下水水位仍不能得到足夠的補(bǔ)給，地下水開采量長(zhǎng)期大于補(bǔ)給量[32]。由于數(shù)據(jù)來源、分析方法和對(duì)地下水補(bǔ)給排泄認(rèn)知的差異，加上欒城地下水水位變化原因十分復(fù)雜，已有研究中對(duì)地下水持續(xù)下降的原因分析結(jié)論不一。本文針對(duì)研究區(qū)地下水補(bǔ)給特點(diǎn)，綜合考慮了欒城本地及其上游的氣候因素和人類活動(dòng)對(duì)欒城地下水水位變化的影響，利用主成分回歸分析法排除各影響因子間的相互作用，得出綜合影響因子及其貢獻(xiàn)率，可為合理開發(fā)地下水資源提供科學(xué)依據(jù)。

由于目前地下水觀測(cè)數(shù)據(jù)不足，對(duì)地下水補(bǔ)給量和補(bǔ)給來源的認(rèn)識(shí)尚存在一定分歧，限制了我們對(duì)地下水和氣候之間動(dòng)態(tài)關(guān)系的理解[36-37]，加之觀測(cè)數(shù)據(jù)樣本量較小，時(shí)間尺度較短，缺少準(zhǔn)確的地下水補(bǔ)徑排和開采量數(shù)據(jù)等原因，給定量分析研究區(qū)地下水水位變化原因帶來一定的不確定性。因此，地下水水位模型構(gòu)建還有待進(jìn)一步考慮其他參數(shù)，如地下水水位的不同埋深范圍、地下水水位的年變化、前期土壤含水量、土壤性質(zhì)(尤其是表層土壤性質(zhì))、季節(jié)以及降水強(qiáng)度等[38]。此外，欒城上游降水量減少和植被覆蓋度增加，使自然降水補(bǔ)給減少的同時(shí)增加了降水量的截留，進(jìn)一步減少了對(duì)下游地下水的側(cè)向補(bǔ)給量；近些年來，華北山前平原種植類型多樣化，土地利用變化對(duì)地下水的影響也有待于進(jìn)一步分析。

3 結(jié) 論

a.1975年以來，欒城降水量增加對(duì)地下水補(bǔ)給具有正效應(yīng)，當(dāng)?shù)厝祟惢顒?dòng)對(duì)地下水水位下降影響最大，同時(shí)，欒城上游人類活動(dòng)對(duì)欒城地下水水位下降的貢獻(xiàn)也不容忽視。1994—2015年，華北山前平原地下水水位變化主要受本地氣候因子、本地人類活動(dòng)和上游人類活動(dòng)三方面影響因子共同作用，各因子對(duì)當(dāng)?shù)氐叵滤换厣呢暙I(xiàn)率分別為11.7%、-50.9%和-37.4%，這三方面影響因子的主導(dǎo)因素分別為欒城降水量、欒城小麥和玉米種植面積和上游石家莊市區(qū)人口。

b.構(gòu)建地下水水位的主成分-時(shí)間序列模型，經(jīng)檢驗(yàn)?zāi)Ｐ偷哪M值與實(shí)測(cè)值趨勢(shì)一致，表明該模型能較好地預(yù)測(cè)華北山前平原典型井灌區(qū)多因子影響下的地下水水位變化。

c.在開發(fā)利用華北平原地下水資源時(shí)，應(yīng)充分考慮上游對(duì)下游的可能影響，盡量減少對(duì)下游水資源補(bǔ)給的不利影響。通過構(gòu)建統(tǒng)計(jì)關(guān)系模型，可以在一定程度上解決目前對(duì)復(fù)雜機(jī)理問題認(rèn)知的不確定性?；谥鞒煞?時(shí)間序列模型可以根據(jù)上一年地下水水位和降水量等數(shù)據(jù)，合理規(guī)劃當(dāng)?shù)剞r(nóng)業(yè)種植結(jié)構(gòu)和種植規(guī)模，統(tǒng)籌考慮預(yù)期糧食產(chǎn)量與地下水開采量，從而有效緩解華北平原地下水長(zhǎng)期超采的嚴(yán)峻局面。