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

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

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

0 引言

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

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

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

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

1 研究區概況

華北平原位于我國東部，分布范圍為E113.0°～ 119.5°，N34.5°～ 40.5°。北靠燕山，西傍太行山，南抵黃河下游，東臨渤海，包括京津2市與河北省的全部平原區，以及山東、河南2省黃河以北的平原區，面積約為14萬km2，總人口約1.33億。華北平原地勢較低，最高海拔為100 m左右，自北、西、南3個方向向渤海灣傾斜，海拔高度逐漸降低[9]。該區氣候屬中緯度大陸性半干旱季風氣候，具有四季分明的特點，多年平均降水量約為550 mm，年內降水多集中在7—9月，占全年降水量的75%左右，冬季降水較少。華北平原有黃河、海河及灤河等大河流，此外還有徒駭河、馬頰河及河北沿海諸河等直接入海的小河流，共有大小河流近60條。研究區位置如圖1所示。

圖1 研究區位置示意圖

2 數據與方法

2.1 數據源

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

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

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

ΔGWS=ΔTWS-ΔSMS。

(1)

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

2.2 EOF分析法

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

(2)

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

然后計算矩陣X的協方差陣C，即

(3)

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

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

(4)

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

(5)

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

(6)

式中PC中每行數據就是對應每個特征向量的時間系數。

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

(7)

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

(8)

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

3 結果與討論

3.1 反演結果驗證

計算華北平原內各個觀測井從2005—2013年間的實測地下水數據均值，與GRACE計算得到的相應時段的GWS進行比較，檢驗數據精度。對比結果如圖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組數據的皮爾遜系數(pearson correlation coefficient)達到了0.754，說明實測地下水數據和計算得到的地下水數據具有較強相關性，表現在數據的周期和年度振幅變化一致性較好。表1中，GRACE反演得到地下水年度振幅均值為135.1 mm，方差為33.7，實測地下水周年振幅均值為 149.2 mm，方差為20.5。

3.2 GWS變化EOF分解

表2為華北平原GWS變化進行EOF分析后得出的前5個模態的方差貢獻率。

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

Tab.2Accumulatedvariancecontributionrateofthefirst5EOFmodelsofGWSinNorthChinaPlain(%)

參數EOF模態12345方差貢獻率79.9311.764.661.680.81累計79.9391.6996.3598.0398.84

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

在95%置信度水平下，按照式(5)計算特征根誤差(圖3)。從特征根誤差范圍看，第1—3模態特征根誤差范圍沒有重疊，存在顯著差別。而第4個特征根及以后的特征根誤差范圍基本上重疊，沒有顯著差別。綜上考慮，本文僅選取前3個模態所對應的典型場，進一步分析華北平原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特征向量場在華北平原內符號相同，表明華北平原地下水儲量變化在空間變化上具有很好的一致性，表現為區域內地下水儲量的增加或者減少，但南北地下水儲量的變化幅度不同，南強北弱。這種全區一致的特性占總體方差的79.93%，說明影響華北平原地下水儲量變化的因素比較單一，第1模態基本上表達了華北平原水儲量變化場的主要結構。正值中心位于冀魯豫3省交界區域，表明該區域的地下水變化幅度最大。從圖4(b)可以看出，第1模態空間向量場對應的權重在逐漸減少，華北平原地下水呈現減少趨勢。但這種多年的趨勢性減少又表現出明顯的年內季節性，即1—6月份表現為儲量減少，而6—12月份表現為儲量增加。這種趨勢性減少、季節性變化與華北平原的地下水開采利用、降水量年內分布特征密切相關。由于長期過量超采[5]，華北平原地下水儲量在多年變化上表現為持續下降。這種下降主要是由于冬小麥的灌溉引起的[5]并導致了上半年地下水儲量減少。在下半年，經過6—9月份的雨季補給，地下水儲量又得到一定程度的恢復。從第1模態的方差貢獻率看，這種自然-人為疊加的變化模式是華北平原GWS變化的主要特征，約占其總體變化的80%。第2特征向量場的方差貢獻率為11.76%，表現為東北—西南格局，反映了沿海與內陸地下水變化呈現相反的趨勢。正值中心在河南北部，0值線出現在研究區域中央，負值沿渤海灣分布。時間序列沒有明顯的趨勢性變化，但年內波動大。這種時空變化可能與沿海-內陸的水文條件變化有關。第3特征向量場方差貢獻率為4.66%，呈現出西北—東南相反的變化趨勢。負值中心集中在太行山脈和燕山山脈的山前平原區，正值集中在渤海灣北部和南部地區。時間序列同樣沒有明顯的趨勢性特征。該模態可能受山前-平原的用水結構、水文地質條件影響，與Huang 等[13]研究結果相似，在山前地區地下水開采強度大，但補給量也大，地下水動態強烈。

4 結論

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

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

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

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

參考文獻(References)：

[1] 張利平,夏 軍,胡志芳.中國水資源狀況與水資源安全問題分析[J].長江流域資源與環境,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] 劉昌明.建設節水型社會 緩解地下水危機[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] 張光輝,費宇紅,劉春華,等.華北平原灌溉用水強度與地下水承載力適應性狀況[J].農業工程學報,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].長江流域資源與環境,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] 張兆吉,雒國中,王 昭,等.華北平原地下水資源可持續利用研究[J].資源科學,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] 魏鳳英.全國夏季降水區域動態權重集成預報試驗[J].應用氣象學報,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] 施 能,章愛國,余錦華.氣象學中使用統計檢驗的幾個重要注記[J].氣象科學,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.