基于GRACE的華北平原地下水儲量時空變化分析

束秋妍， 潘 云， 宮輝力， 黃志勇， 鄭龍群

(1.首都師范大學(xué)資源環(huán)境與旅游學(xué)院，北京 100048； 2.首都師范大學(xué)城市環(huán)境過程與數(shù)字模擬國家重點實驗室培育基地，北京 100048)

0 引言

華北平原為我國第二大平原，也是我國的政治中心所在，人口稠密，是以旱作為主的農(nóng)業(yè)區(qū)，擁有全國最大的小麥、玉米種植區(qū)，是全國重要的糧食基地。華北平原也是我國北方缺水較嚴(yán)重的地區(qū)，人均水資源量僅有450 m3，低于國際公認(rèn)的極度缺水的標(biāo)準(zhǔn)(500 m3)[1]。地下水是區(qū)內(nèi)主要供水源，隨著工農(nóng)業(yè)的生產(chǎn)發(fā)展及人口數(shù)量的不斷增加，地下水的開采強度和開采規(guī)模不斷增大。據(jù)統(tǒng)計，目前華北平原有地下水開采井約200萬眼，從1985年到20世紀(jì)末，多年平均開采量已超過1010m3[2]，且大部分用于農(nóng)業(yè)灌溉。研究華北平原地區(qū)地下水儲量變化的空間分布特征及時間變化規(guī)律，對于該區(qū)地下水資源可持續(xù)利用具有重要意義。

傳統(tǒng)監(jiān)測地下水儲量變化一般利用監(jiān)測井進(jìn)行檢測。這種實地測量地下水的方法僅依靠單點的地下水位，難以充分反映地下水儲量在空間上和時間上的變化特征。GRACE(gravity recovery and climate experiment)重力衛(wèi)星作為一種監(jiān)測地表水儲量變化的新型手段，打破了傳統(tǒng)地面觀測在時間和空間上的局限性，并已在國內(nèi)外地下水水儲量變化研究中得到了應(yīng)用[3-5]。

經(jīng)驗正交函數(shù)(empirical orthogonal function，EOF)分析方法，通過分析矩陣數(shù)據(jù)中的結(jié)構(gòu)特征，提取數(shù)據(jù)主要特征量，在水儲量研究方面應(yīng)用廣泛。Awange等[6]利用EOF方法研究了埃塞俄比亞地區(qū)的陸地水儲量、土壤濕度和降水的時空變化特征，并分析了三者的相關(guān)性及陸地水儲量變化對降水的滯后性； Kang等[7]利用重力衛(wèi)星數(shù)據(jù)結(jié)合EOF分析了中國2003—2010年間陸地水儲量變化特征，發(fā)現(xiàn)水儲量季節(jié)性變化與降水關(guān)系密切，其中梅雨季節(jié)和重大干旱時期在EOF分解的模態(tài)中都有體現(xiàn)； 閻福禮等[8]運用EOF方法分解了基于GRACE得到的長江流域水儲量變化，較好地解釋了該流域水儲量變化的空間規(guī)律。

本文基于GRACE重力衛(wèi)星反演華北平原地下水儲量變化，并用實測數(shù)據(jù)加以驗證。運用EOF對華北平原2003—2015年間地下水儲量變化進(jìn)行分解，分析地下水儲量變化的時間和空間特征，并探討控制其時空變化的主要因素。

1 研究區(qū)概況

華北平原位于我國東部，分布范圍為E113.0°～ 119.5°，N34.5°～ 40.5°。北靠燕山，西傍太行山，南抵黃河下游，東臨渤海，包括京津2市與河北省的全部平原區(qū)，以及山東、河南2省黃河以北的平原區(qū)，面積約為14萬km2，總?cè)丝诩s1.33億。華北平原地勢較低，最高海拔為100 m左右，自北、西、南3個方向向渤海灣傾斜，海拔高度逐漸降低[9]。該區(qū)氣候?qū)僦芯暥却箨懶园敫珊导撅L(fēng)氣候，具有四季分明的特點，多年平均降水量約為550 mm，年內(nèi)降水多集中在7—9月，占全年降水量的75%左右，冬季降水較少。華北平原有黃河、海河及灤河等大河流，此外還有徒駭河、馬頰河及河北沿海諸河等直接入海的小河流，共有大小河流近60條。研究區(qū)位置如圖1所示。

圖1 研究區(qū)位置示意圖

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

2.1 數(shù)據(jù)源

將德國地學(xué)研究中心發(fā)布的GRACE重力場模型產(chǎn)品減去一個平均重力場(或背景重力場)，便可獲取重力異?；驎r變部分。本文減去的是2005年1月至2010年12月平均重力場。所采用核函數(shù)格網(wǎng)分辨率為1°，經(jīng)球諧展開并截斷至60階，去條帶、高斯濾波處理(半徑為200 km)后得到2003年1月—2015年8月期間152個月的陸地水儲量變化(terrestrial water storage, TWS)數(shù)據(jù)，缺省月份的數(shù)據(jù)在實驗處理的過程中一般取相鄰數(shù)據(jù)的均值來替代。

由美國航空航天局戈達(dá)德空間飛行中心和美國海洋和大氣局國家環(huán)境預(yù)報中心聯(lián)合開發(fā)的全球陸面數(shù)據(jù)同化系統(tǒng) (global land data assimilation system version1，GLDAS-1)，可獲取高精度、高空間分辨率的陸表水儲量和能量信息。該系統(tǒng)驅(qū)動4個陸面過程模型： 通用陸面模型(common land model，CLM)，Noah，Mosaic以及可變?nèi)霛B能力模式(variable infiltration capacity，VIC)。本文采用2003年1月—2015年8月期間CLM模型的每月1°×1°土壤水儲量(soil moisture storage，SMS)數(shù)據(jù)集。

從GRACE得到的TWS，包括地表水(湖泊、水庫)、土壤水、地下水、雪水當(dāng)量和植被冠層含水量等。反演地下水儲量需考慮研究區(qū)域水文地質(zhì)條件，扣除非地下水量。對于干旱半干旱區(qū)，土壤水儲量是主要的非地下水部分。本文結(jié)合GRACE監(jiān)測的總陸地水儲量變化ΔTWS和水文模型得到的土壤水儲量變化ΔSMS分離出地下水儲量(ground water storage，GWS)變化，即

ΔGWS=ΔTWS-ΔSMS。

(1)

本文采用的實測地下水位數(shù)據(jù)來源于中國地質(zhì)監(jiān)測院編著的地下水位年鑒，共搜集了2005年1月—2013年12月期間年鑒中連續(xù)監(jiān)測地下水位觀測井?dāng)?shù)據(jù)，共64個； 采用的降水?dāng)?shù)據(jù)來自中國氣象數(shù)據(jù)網(wǎng)(http: //data.cma.cn/)，選取了2003年1月—2014年12月期間華北平原境內(nèi)18個連續(xù)監(jiān)測降水站點的觀測數(shù)據(jù)。

2.2 EOF分析法

EOF分析也稱經(jīng)驗正交函數(shù)分解，基本原理是對包含若干個空間點的變量隨時間進(jìn)行分解。設(shè)樣本容量為m個點的資料，則空間中任意一點i和任一時間點j的變量值xij可看成由n個空間函數(shù)eofik和時間函數(shù)tkj(k=1 , 2 , … ,n)的線性組合[10]。將得到的GWS距平場以矩陣的形式表示如下

(2)

式中：m為空間點的個數(shù)，即空間網(wǎng)格點和測站等；n為時間樣本數(shù)；X中第j列xj=(x1j，x2j，…，xmj)T即為第j個空間場。

然后計算矩陣X的協(xié)方差陣C，即

(3)

計算C的特征根(λ1，…，λm)和特征向量Vm×m，二者滿足

Cm×mVm×m=Vm×mΛm×m,

(4)

式中Λ是m×m維對角陣，即

(5)

將特征值λ從大到小順序排列，λ1>λ2>…>λm，每個非0的特征根對應(yīng)一列特征向量值，也稱eof。將特征向量投影到原始資料矩陣X上，就得到所有空間特征向量對應(yīng)的時間系數(shù)(即主成分)，表達(dá)式為

(6)

式中PC中每行數(shù)據(jù)就是對應(yīng)每個特征向量的時間系數(shù)。

通常情況下，用前幾個模態(tài)就可以近似地反映變量場的主要時空變化。模態(tài)的重要性可以用特征根的大小來判斷，即特征根越大則其所對應(yīng)的模態(tài)也越重要，同樣表明對總方差的貢獻(xiàn)率越高。第i個模態(tài)對整個變量場的貢獻(xiàn)率P為

(7)

實際資料數(shù)據(jù)通過EOF分解后得到的空間模態(tài)是否是隨機(jī)的，是否有物理意義，需要進(jìn)行顯著性檢驗[11]。在95%置信度水平下特征根的誤差可通過式(8)得到檢驗，即

(8)

式中：Δλ為特征根誤差；λ為各模態(tài)特征根；N*為自由度。當(dāng)相鄰的2個特征值滿足λj=λj+1≥Δλ時，就認(rèn)為這2個特征根可以分離，它們之間存在顯著差別，其所對應(yīng)的EOF才具有物理意義，否則將被排除。

3 結(jié)果與討論

3.1 反演結(jié)果驗證

計算華北平原內(nèi)各個觀測井從2005—2013年間的實測地下水?dāng)?shù)據(jù)均值，與GRACE計算得到的相應(yīng)時段的GWS進(jìn)行比較，檢驗數(shù)據(jù)精度。對比結(jié)果如圖2和表1所示。

圖2 GRACE反演和實測GWS變化時間序列

周年GWS均值2005年2006年2007年2008年2009年2010年2011年2012年2013年均值方差GRACE反演值126.5189.2136.288.1108.2106.2124.6172.5164.5135.133.7觀測井實測值182.1148.6173.8111.1152.4142.0145.7150.7136.7149.220.5

從圖2可以看出，2組數(shù)據(jù)的皮爾遜系數(shù)(pearson correlation coefficient)達(dá)到了0.754，說明實測地下水?dāng)?shù)據(jù)和計算得到的地下水?dāng)?shù)據(jù)具有較強相關(guān)性，表現(xiàn)在數(shù)據(jù)的周期和年度振幅變化一致性較好。表1中，GRACE反演得到地下水年度振幅均值為135.1 mm，方差為33.7，實測地下水周年振幅均值為 149.2 mm，方差為20.5。

3.2 GWS變化EOF分解

表2為華北平原GWS變化進(jìn)行EOF分析后得出的前5個模態(tài)的方差貢獻(xiàn)率。

表2華北平原地下水儲量變化EOF分解后前5個模態(tài)累積方差貢獻(xiàn)率

Tab.2Accumulatedvariancecontributionrateofthefirst5EOFmodelsofGWSinNorthChinaPlain(%)

參數(shù)EOF模態(tài)12345方差貢獻(xiàn)率79.9311.764.661.680.81累計79.9391.6996.3598.0398.84

可以從表2中看出，模態(tài)1的方差貢獻(xiàn)最大達(dá)到79.93%，遠(yuǎn)大于其他特征向量的方差貢獻(xiàn)率，第2模態(tài)和第3模態(tài)方差貢獻(xiàn)急劇收斂，分別為11.76%和4.66%，前3個模態(tài)的累積貢獻(xiàn)率達(dá)到96.35%，基本可解釋大部分變化特征，其余模態(tài)所占比重較小。而且當(dāng)相鄰的2個特征根過于接近時，EOF采樣變異性會明顯增大[12]。

在95%置信度水平下，按照式(5)計算特征根誤差(圖3)。從特征根誤差范圍看，第1—3模態(tài)特征根誤差范圍沒有重疊，存在顯著差別。而第4個特征根及以后的特征根誤差范圍基本上重疊，沒有顯著差別。綜上考慮，本文僅選取前3個模態(tài)所對應(yīng)的典型場，進(jìn)一步分析華北平原GWS變化時空變化特征。

圖3 華北平原地下水儲量變化EOF分解特征根

3.3 GWS變化主要時空特征

前3個特征向量場及其分析如圖4所示。

(a) GWS變化第1特征向量場 (b) 第一時間序列13點平滑值

(c) GWS變化第2特征向量場 (d) 第二時間序列13點平滑值

(e) GWS變化第3特征向量場 (f) 第三時間序列13點平滑值

圖42003—2015年間華北平原GWS變化3個特征向量分析

Fig.4ThreefeaturevectorfieldsofvarietyofGWSfrom2003to2015inNorthChinaPlain

從圖4中可以看出，第1特征向量場在華北平原內(nèi)符號相同，表明華北平原地下水儲量變化在空間變化上具有很好的一致性，表現(xiàn)為區(qū)域內(nèi)地下水儲量的增加或者減少，但南北地下水儲量的變化幅度不同，南強北弱。這種全區(qū)一致的特性占總體方差的79.93%，說明影響華北平原地下水儲量變化的因素比較單一，第1模態(tài)基本上表達(dá)了華北平原水儲量變化場的主要結(jié)構(gòu)。正值中心位于冀魯豫3省交界區(qū)域，表明該區(qū)域的地下水變化幅度最大。從圖4(b)可以看出，第1模態(tài)空間向量場對應(yīng)的權(quán)重在逐漸減少，華北平原地下水呈現(xiàn)減少趨勢。但這種多年的趨勢性減少又表現(xiàn)出明顯的年內(nèi)季節(jié)性，即1—6月份表現(xiàn)為儲量減少，而6—12月份表現(xiàn)為儲量增加。這種趨勢性減少、季節(jié)性變化與華北平原的地下水開采利用、降水量年內(nèi)分布特征密切相關(guān)。由于長期過量超采[5]，華北平原地下水儲量在多年變化上表現(xiàn)為持續(xù)下降。這種下降主要是由于冬小麥的灌溉引起的[5]并導(dǎo)致了上半年地下水儲量減少。在下半年，經(jīng)過6—9月份的雨季補給，地下水儲量又得到一定程度的恢復(fù)。從第1模態(tài)的方差貢獻(xiàn)率看，這種自然-人為疊加的變化模式是華北平原GWS變化的主要特征，約占其總體變化的80%。第2特征向量場的方差貢獻(xiàn)率為11.76%，表現(xiàn)為東北—西南格局，反映了沿海與內(nèi)陸地下水變化呈現(xiàn)相反的趨勢。正值中心在河南北部，0值線出現(xiàn)在研究區(qū)域中央，負(fù)值沿渤海灣分布。時間序列沒有明顯的趨勢性變化，但年內(nèi)波動大。這種時空變化可能與沿海-內(nèi)陸的水文條件變化有關(guān)。第3特征向量場方差貢獻(xiàn)率為4.66%，呈現(xiàn)出西北—東南相反的變化趨勢。負(fù)值中心集中在太行山脈和燕山山脈的山前平原區(qū)，正值集中在渤海灣北部和南部地區(qū)。時間序列同樣沒有明顯的趨勢性特征。該模態(tài)可能受山前-平原的用水結(jié)構(gòu)、水文地質(zhì)條件影響，與Huang 等[13]研究結(jié)果相似，在山前地區(qū)地下水開采強度大，但補給量也大，地下水動態(tài)強烈。

4 結(jié)論

本文利用GRACE數(shù)據(jù)反演了華北平原2003年1月—2015年8月地下水儲量變化，采用EOF方法對此進(jìn)行時空分解，得到反映數(shù)據(jù)特征的前3個特征向量場，并對特征向量場和時間系數(shù)進(jìn)行分析，得到如下結(jié)論：

1)華北平原地下水儲量變化可以分解為3個主要模態(tài)，其對總體變化的解釋率達(dá)到96.35%。

2)第1模態(tài)空間變化一致，表現(xiàn)出多年趨勢性減少與年內(nèi)季節(jié)性變化相結(jié)合的特征，可能由研究區(qū)內(nèi)地下水開采、年內(nèi)降水分布共同作用導(dǎo)致。這也是目前華北平原地下水儲量減少的主要模式，對整體變化的解釋率約為80%。

3)第2和第3模態(tài)分別表現(xiàn)出了東北—西南和西北—東南2種變化相反的空間格局，對總體變化的解釋率分別約為12%和5%，在時間上沒有明顯的趨勢性變化，推測可能主要受沿海-內(nèi)陸、山前-平原的水循環(huán)和水文地質(zhì)條件控制。

參考文獻(xiàn)(References)：

[1] 張利平,夏 軍,胡志芳.中國水資源狀況與水資源安全問題分析[J].長江流域資源與環(huán)境,2009,18(2):116-120.

Zhang L P,Xia J,Hu Z F.Situation and problem analysis of water resource security in China[J].Resources and Environment in the Yangtze Basin,2009,18(2):116-120.

[2] 劉昌明.建設(shè)節(jié)水型社會 緩解地下水危機(jī)[J].中國水利,2007(15):10-13.

Liu C M.Building water-saving society and alleviating groundwater crisis[J].China Water Resources,2007(15):10-13.

[3] Yeh P J F,Swenson S C,Famiglietti J S,et al.Remote sensing of groundwater storage changes in Illinois using the gravity recovery and climate experiment (GRACE)[J].Water Resources Research,2006,42(12)W 12203:1-7.doi:10.1029/2006WR005374.

[4] Feng W,Zhong M,Lemoine J M,et al.Evaluation of groundwater depletion in North China using the gravity recovery and climate experiment (GRACE) data and ground-based measurements[J].Water Resources Research,2013,49(4):2110-2118.

[5] 張光輝,費宇紅,劉春華,等.華北平原灌溉用水強度與地下水承載力適應(yīng)性狀況[J].農(nóng)業(yè)工程學(xué)報,2013,29(1):1-10.

Zhang G H,Fei Y H,Liu C H,et al.Adaptation between irrigation intensity and groundwater carrying capacity in North China Plain[J].Transactions of the Chinese Society of Agricultural Engineering,2013,29(1):1-10.

[6] Awange J L,Gebremichael M,Forootan E,et al.Characterization of Ethiopian mega hydrogeological regimes using GRACE,TRMM and GLDAS datasets[J].Advances in Water Resources,2014,74:64-78.

[7] Kang K X,Li H,Peng P,et al.Low-frequency variability of terrestrial water budget in China using GRACE satellite measurements from 2003 to 2010[J].Geodesy and Geodynamics,2015,6(6):444-452.

[8] 閻福禮,李書明,王世新,等.基于EOF方法長江流域2002—2013年GRACE水儲量時空變化研究[J].長江流域資源與環(huán)境,2015,24(S1):131-137.

Yan F L,Li S M,Wang S X,et al.Temporal and spatial variations research of GRACE water storage changes over the Yangtze River Bisin[J].Resources and Environment in the Yangtze Basin,2015,24(S1):131-137.

[9] 張兆吉,雒國中,王 昭,等.華北平原地下水資源可持續(xù)利用研究[J].資源科學(xué),2009,31(3):355-360.

Zhang Z J,Luo G Z,Wang Z,et al.Study on sustainable utilization of groundwater in North China Plain[J].Resources Science,2009,31(3):355-360.

[10] 魏鳳英.全國夏季降水區(qū)域動態(tài)權(quán)重集成預(yù)報試驗[J].應(yīng)用氣象學(xué)報,1999,10(4):402-409.

Wei F Y.Regional consensus forecast method with dynamic weighting for summer precipitation over China[J].Quarterly Journal of Applied Meteorlolgy,1999,10(4):402-409.

[11] 施 能,章愛國,余錦華.氣象學(xué)中使用統(tǒng)計檢驗的幾個重要注記[J].氣象科學(xué),2009,29(5):670-673.

Shi N,Zhang A G,Yu J H,Some important problems of the statistics test in meteorology[J].Scientia Meteorologica Sinica,2009,29(5):670-673.

[12] North G R,Bell T L,Cahalan R F,et al.Sampling errors in the estimation of empirical orthogonal functions[J].Monthly Weather Review,1982,110(7):699-699.

[13] Huang Z Y,Pan Y,Gong H L,et al.Subregional-scale groundwater depletion detected by GRACE for both shallow and deep aquifers in North China Plain[J].Geophysical Research Letters,2015,42(6):1791-1799.