張明剛 魏長(zhǎng)壽

1 山東科技大學(xué)測(cè)繪與空間信息學(xué)院，青島市前灣港路579號(hào)，266590 2 內(nèi)蒙古科技大學(xué)礦業(yè)與煤炭學(xué)院，內(nèi)蒙古自治區(qū)包頭市阿爾丁大街7號(hào)，014010

近20 a來，華北平原農(nóng)業(yè)迅猛發(fā)展，地下水超量開采，造成地下水資源稀缺。因此，獲取華北平原地區(qū)的地下水儲(chǔ)量變化，探究影響其變化的因素對(duì)當(dāng)?shù)剞r(nóng)業(yè)的良好運(yùn)營(yíng)非常重要。

對(duì)華北平原地下水儲(chǔ)量進(jìn)行監(jiān)測(cè)的研究較多[1-3]，但一般側(cè)重于獲取地下水水儲(chǔ)量的變化趨勢(shì)，減少后處理誤差，很少關(guān)注農(nóng)業(yè)用水量變化、降水量變化的影響，缺乏對(duì)地下水減少影響因素的具體分析。本文使用GRACE衛(wèi)星數(shù)據(jù)質(zhì)量較好的2003～2015年的數(shù)據(jù)反演華北平原地下水儲(chǔ)量的時(shí)空變化，并利用經(jīng)驗(yàn)正交函數(shù)(empirical orthogonal function,EOF)具體分析影響華北平原地下水儲(chǔ)量時(shí)空變化的主要因素。

1 研究區(qū)概況

華北平原位于31°36′～40°29′N、112°13′～120°53′E，屬半干旱、半濕潤(rùn)季風(fēng)性氣候，春季干旱少雨、蒸發(fā)強(qiáng)，夏季高溫多雨，冬季寒冷干燥，地區(qū)年均降水量在空間上分布不均。冬小麥和夏玉米的一年兩熟制是該地區(qū)的典型農(nóng)業(yè)種植模式[4]。其中，玉米的灌漿期在7月份左右，此時(shí)華北平原氣溫較高，雨水供應(yīng)充足，地下水開采對(duì)地下水儲(chǔ)量的變化影響較??；華北平原冬小麥的產(chǎn)量占全國(guó)的60%～80%，高密度的冬小麥種植需要消耗大量的水資源，且其灌漿期為3～5月，恰逢干旱少雨的春季，導(dǎo)致水資源嚴(yán)重短缺，此時(shí)超過60%的用水來自地下水[5]。

圖1 研究區(qū)位置Fig.1 Location of the study area

2 數(shù)據(jù)與方法

2.1 監(jiān)測(cè)井?dāng)?shù)據(jù)

從中國(guó)地質(zhì)環(huán)境監(jiān)測(cè)研究院獲取2003～2013年華北平原地區(qū)61口監(jiān)測(cè)井水位數(shù)據(jù)。為了將水位數(shù)據(jù)轉(zhuǎn)化為與GRACE等效的水儲(chǔ)量變化，參考文獻(xiàn)[6]，取華北平原的平均給水度為0.06，將水位數(shù)據(jù)乘以區(qū)域給水度求得水儲(chǔ)量變化的等效水柱高。

2.2 GRACE重力衛(wèi)星數(shù)據(jù)獲取陸地水儲(chǔ)量變化

使用CSR提供的GRACE level-2 RL06數(shù)據(jù)產(chǎn)品，時(shí)間跨度為2003-01～2015-12，分辨率為1個(gè)月，使用三次樣條差值方法補(bǔ)齊缺失數(shù)據(jù)[7]。對(duì)GRACE數(shù)據(jù)進(jìn)行如下預(yù)處理：1)用SLR對(duì)應(yīng)項(xiàng)替代原始的1階項(xiàng)和C20項(xiàng)；2)對(duì)生成的全球格網(wǎng)數(shù)據(jù)進(jìn)行去平均處理；3)采用P4M15去相關(guān)濾波去除“南-北條帶誤差”；4)采用300 km的高斯濾波對(duì)去條帶之后的球諧系數(shù)進(jìn)行平滑，降低高階項(xiàng)球諧系數(shù)的噪聲；5)采用正向建模法改正信號(hào)的泄露誤差[8]。最后反演獲得陸地水儲(chǔ)量(terrestrial water storage，TWS)變化。

2.3 GLDAS水文模型數(shù)據(jù)

選取2003～2015年共156個(gè)月GLDAS Noah陸地模型輸出的0～200 cm土壤含水量數(shù)據(jù)和雪水當(dāng)量月解數(shù)據(jù)來估計(jì)地表水儲(chǔ)量的變化，其空間分辨率為1°×1°，并采用與GRACE數(shù)據(jù)相同的處理方式獲取土壤水儲(chǔ)量(soil moisture storage，SMS)變化與雪水當(dāng)量(snow water equivation,SWE)變化。

2.4 地下水儲(chǔ)量變化ΔGWS的獲取

簡(jiǎn)化后的區(qū)域水量平衡方程為[9]：

ΔGWS=ΔTWS-ΔSMS-ΔSWE

利用GRACE獲取的ΔTWS，減去GLDAS水文模型數(shù)據(jù)獲取的ΔSWS和ΔSWE，即可得地下水儲(chǔ)量(ground water storage，GWS)的變化ΔGWS。

2.5 TRMM衛(wèi)星降雨數(shù)據(jù)

降雨對(duì)地下水儲(chǔ)量研究具有重要意義。本文使用熱帶降雨測(cè)量任務(wù)(tropical rainfall measuring mission,TRMM)衛(wèi)星數(shù)據(jù)(https:∥disc.gsfc.nasa.gov/datasets/TRMM_3B43_7/summary)估算華北平原的年均降水量。

2.6 小麥藍(lán)水足跡

小麥藍(lán)水足跡WFblue指在小麥生長(zhǎng)過程中用于灌溉消耗的地表與地下水的總量[10-11]。由于華北平原降水量在季節(jié)上分布不均，使得冬小麥在生長(zhǎng)期間主要依靠地下水補(bǔ)充灌溉。因此，本文利用冬小麥WFblue數(shù)據(jù)來說明華北平原年均灌溉用水量的空間變化及其對(duì)地下水儲(chǔ)量的影響。冬小麥WFblue數(shù)據(jù)可以從水足跡網(wǎng)站(https://waterfootprint.org/en/resources/waterstat/wf-crop-production-and-consumption-china/)獲取，因網(wǎng)站數(shù)據(jù)僅截止到2009年，故本文使用的冬小麥WFblue數(shù)據(jù)時(shí)間跨度為2003～2009年。

2.7 EOF的原理

EOF的基本原理是對(duì)包含p個(gè)空間點(diǎn)的變量隨時(shí)間進(jìn)行分解。假設(shè)有n個(gè)點(diǎn)的資料，則任一空間點(diǎn)i和任一時(shí)間點(diǎn)j的變量值(xij)可以看成由p個(gè)空間函數(shù)EOFik和時(shí)間函數(shù)tkj(k=1,2,…,p)的線性組合。具體原理參見文獻(xiàn)[12]，本文不再贅述。分解后使用North檢驗(yàn)對(duì)各模態(tài)進(jìn)行檢核，以保證各模態(tài)相互獨(dú)立[13]。

3 結(jié)果與分析

3.1 反演結(jié)果驗(yàn)證

圖2、3分別為2003～2013年月尺度和季尺度GRACE獲取的ΔGWS與監(jiān)測(cè)井實(shí)測(cè)ΔGWS的時(shí)變序列，二者振幅和趨勢(shì)大致相同。從圖2可知，月尺度下，2個(gè)時(shí)間序列之間的相關(guān)系數(shù)為0.725 3；從圖3可知，季尺度下，2個(gè)時(shí)變序列之間的相關(guān)系數(shù)為0.793 3。由此可得，基于GRACE計(jì)算得到的地下水儲(chǔ)量變化與實(shí)測(cè)數(shù)據(jù)之間在時(shí)間上有較強(qiáng)的相關(guān)性，用GRACE研究華北平原的地下水儲(chǔ)量變化是可行的，且時(shí)間尺度越大，結(jié)果越準(zhǔn)確。

圖2 月尺度GRACE反演和實(shí)測(cè)ΔGWS時(shí)間序列Fig.2 Monthly time series of ΔGWS inversed by GRACE and measured ΔGWS

圖3 季尺度GRACE反演和實(shí)測(cè)ΔGWS時(shí)間序列Fig.3 Seasonal time series of ΔGWS inversed by GRACE and measured ΔGWS

在2003～2005年初，GRACE獲取的地下水儲(chǔ)量變化相比監(jiān)測(cè)井實(shí)測(cè)數(shù)據(jù)整體偏低，這可能是由于現(xiàn)有的監(jiān)測(cè)井?dāng)?shù)量有限、有些月份缺失、空間分辨率不足所致。

3.2 地下水儲(chǔ)量變化影響因素分析

圖4為2003～2013年監(jiān)測(cè)井實(shí)測(cè)ΔGWS、2003～2015年GRACE獲取的ΔGWS、降水以及2003～2009年冬小麥WFblue隨時(shí)間的變化。其中，地下水對(duì)降水的響應(yīng)延遲2～3個(gè)月，故將降水時(shí)序后移3個(gè)月。

圖4 監(jiān)測(cè)井實(shí)測(cè)ΔGWS、GRACE反演的ΔGWS、冬小麥WFblue和降水隨時(shí)間的變化Fig.4 ΔGWS obtained from monitoring wells,ΔGWS inversed by GRACE,winter wheat WFblue and precipitation vary with time

可以看出，2003～2015年華北平原地下水儲(chǔ)量以2.62±0.40 cm/a的速率銳減，可分為2個(gè)階段：1)2003～2008年地下水儲(chǔ)量處于上升狀態(tài)，速率為3.50±1.47 cm/a；2)2009年開始，地下水儲(chǔ)量由增轉(zhuǎn)減，2009～2015年地下水儲(chǔ)量的下降速率為4.44±0.54 cm/a，2014～2015年下降最快，達(dá)到8.96±4.52 cm/a。

圖4還表明，2003～2005年降水較多，冬小麥WFblue處于較低值，冬小麥灌溉對(duì)地下水的依賴程度低，因此地下水儲(chǔ)量在此期間處于上升階段；2006年降水較少，不能滿足冬小麥灌溉的需要，地下水被大量用于灌溉生產(chǎn)，冬小麥WFblue處于2003～2009年的最大值；到了干旱少雨的2009年，降水峰值未超過200 mm，灌溉超采情況進(jìn)一步加劇，地下水儲(chǔ)量開始急劇下降；2012年和2013年降雨量有所上升，但因農(nóng)業(yè)灌溉常年超采地下水，導(dǎo)致地下水資源補(bǔ)給、排泄的平衡難以恢復(fù)，地下水儲(chǔ)量并未回升；2014年和2015年降水量最少，使得地下水儲(chǔ)量進(jìn)一步減少，地下水儲(chǔ)量下降速率達(dá)到峰值。由以上分析可知，冬小麥的灌溉活動(dòng)和降水與華北平原地下水儲(chǔ)量下降關(guān)系密切。

3.3 年際ΔGWS的EOF分解

為清晰地展現(xiàn)華北平原地下水儲(chǔ)量的年際變化，進(jìn)一步分析影響地下水儲(chǔ)量變化的主要因素，將GRACE獲取的年際ΔGWS進(jìn)行EOF分解，表1為EOF分解后前5個(gè)模態(tài)的方差貢獻(xiàn)率?？梢钥闯觯?個(gè)模態(tài)的累計(jì)方差貢獻(xiàn)率達(dá)到93.09%，基本可以解釋總體的地下水儲(chǔ)量變化特征。

表1 華北平原年際ΔGWS EOF分解后前5個(gè)模態(tài)方差貢獻(xiàn)率Tab.1 Accumulated variance contribution rate of the frist five EOF modes of ΔGWS in North China plain

表2為進(jìn)行North檢驗(yàn)的前5個(gè)模態(tài)的特征值λ、特征值的誤差范圍ej和檢驗(yàn)值λj-λj+1。從表2可以看出，EOF分析結(jié)果均符合North檢驗(yàn)標(biāo)準(zhǔn)，所對(duì)應(yīng)的EOF是有價(jià)值的。綜上考慮，本文選取前2個(gè)整體模態(tài)，結(jié)合2003～2009年冬小麥年均WFblue數(shù)據(jù)與2003～2015年年均降水?dāng)?shù)據(jù)進(jìn)一步分析華北平原ΔGWS的時(shí)空變化特征。

表2 North檢驗(yàn)結(jié)果Tab.2 Results of North test

圖5為第1特征向量與冬小麥年均WFblue對(duì)比。從圖5(a)中可以看出，2003～2009年第1特征向量場(chǎng)與研究區(qū)冬小麥年均WFblue空間分布相關(guān)系數(shù)為-0.69，呈現(xiàn)良好的負(fù)相關(guān)性。第1特征向量場(chǎng)在黃河流域與海河流域的大部分地區(qū)為負(fù)值，在邯鄲、邢臺(tái)、鶴壁、濮陽(yáng)、新鄉(xiāng)一帶ΔGWS長(zhǎng)期處于低值區(qū)，地下水虧損嚴(yán)重，與此區(qū)域較高的冬小麥WFblue相對(duì)應(yīng)；只有在淮河流域的南部第1特征向量場(chǎng)為正值，地下水儲(chǔ)量增加，與農(nóng)業(yè)灌溉開采地下水較少有關(guān)。不同區(qū)域地下水儲(chǔ)量的變化幅度不同，黃河流域地下水儲(chǔ)量下降最快，以黃河為界，向南北兩側(cè)下降速率減小。從圖5(b)中可以看出，2003～2009年第1時(shí)間序列與冬小麥年均WFblue時(shí)變序列呈現(xiàn)良好的一致性。這種與冬小麥年均WFblue一致的特性占總體方差的80.04%，基本表達(dá)了地下水儲(chǔ)量變化場(chǎng)的主要結(jié)構(gòu)，表明冬小麥的灌溉用水是影響華北平原地下水儲(chǔ)量變化的主要因素。

圖5 第1特征向量與冬小麥年均WFblue對(duì)比Fig.5 Comparison of the first feature vector and annual average WFblue of the winter wheat

圖6(a)為2003～2015年第2特征向量場(chǎng)與年均降水的空間分布?？梢钥闯?，華北平原的年均降水量由北向南隨緯度遞減，在南部的淮河流域降雨量較大，而中部的黃河下游平原和北部的京津冀地區(qū)則較少?？傮w而言，降水量空間分布不均，以黃河為界，呈現(xiàn)出南多北少的格局，在黃河以北地區(qū)，太行山山前平原整體上降水量相對(duì)較少。這與第2特征向量場(chǎng)有著強(qiáng)一致性，二者之間的相關(guān)系數(shù)為0.93。結(jié)合圖6(b)可知，二者的時(shí)變序列趨勢(shì)保持一致，表明降水是影響華北平原地下水儲(chǔ)量年際變化的次要因素。

圖6 第2特征向量與年均降水對(duì)比Fig.6 Comparison of the second fecture vector and annual average precipitation

4 結(jié) 語(yǔ)

本文利用GRACE重力衛(wèi)星數(shù)據(jù)反演了華北平原2003～2015年地下水儲(chǔ)量的年際變化，使用EOF分解方法進(jìn)行時(shí)空分解，并將第1特征向量和第2特征向量結(jié)合冬小麥年均WFblue和降水進(jìn)行對(duì)比分析，得到如下結(jié)論：

1)GRACE反演結(jié)果與監(jiān)測(cè)井實(shí)測(cè)數(shù)據(jù)相關(guān)性好，證明了使用GRACE與GLDAS水文模型結(jié)合反演地下水儲(chǔ)量時(shí)空變化的可行性。

2)2003～2015年華北平原地下水儲(chǔ)量以2.62±0.40 cm/a的速率下降，使用EOF將華北平原地下水儲(chǔ)量的年際變化分解為2個(gè)主要模態(tài)，第1模態(tài)占總體變化的80.04%，第2模態(tài)占總體變化的13.05%，二者對(duì)總體變化的解釋率可以達(dá)到93.09%。

3)第1特征向量場(chǎng)呈現(xiàn)以黃河流域?yàn)橹行摹⑾蚰媳眱蓚?cè)逐漸變緩的趨勢(shì)，與2003～2009年冬小麥年均WFblue空間分布之間有良好的負(fù)相關(guān)性，相關(guān)系數(shù)為-0.69，可見冬小麥的灌溉用水是華北平原地下水儲(chǔ)量減少的主要因素。

4)第2特征向量場(chǎng)方差貢獻(xiàn)率為13.05%，在時(shí)間上沒有明顯的趨勢(shì)性變化，由北向南隨緯度遞減，與年均降水的空間分布之間的相關(guān)系數(shù)為0.93，可知降水是影響華北平原地下水儲(chǔ)量變化的第二大因素。