基于GRACE和MODIS數(shù)據(jù)的長江流域陸地水儲量變化分析及預測

，，，

(河海大學 水利水電學院,南京 210098)

1 研究背景

陸地水儲量(Terrestrial Water Storage，TWS)是水文循環(huán)的一個重要組成部分，是地表水、地下水、土壤含水量、冰雪和生物體含水量變化的綜合體現(xiàn)。陸地水儲量變化(Terrestrial Water Storage Anomaly，TWSA)的準確估計和趨勢分析，對理解全球和局域水循環(huán)過程、研究和預報氣候的變化、農(nóng)業(yè)生產(chǎn)的安排、洪澇以及其他自然災害的防治都具有重要意義[1]。然而，由于土壤水、地下水等要素的地面觀測站點稀疏，可見光-紅外及微波遙感數(shù)據(jù)只能得到表層土壤含水量，利用傳統(tǒng)方法難以對各分量進行有效觀測，水儲量變化研究進展緩慢。在此背景下，2002年3月份發(fā)射的GRACE(Gravity Recovery and Climate Experiment)重力衛(wèi)星，通過監(jiān)測全球重力場的變化，綜合反映了氣候波動、全球變化和人類活動對水資源儲量的影響，為陸地水儲量變化研究提供了新方法[2]。目前，我國對GRACE重力衛(wèi)星數(shù)據(jù)反演水儲量變化的流域尺度研究，已涉及到長江、黃河、海河、黑河以及雅魯藏布江等流域[3-9]，其中部分研究通過與全球陸面數(shù)據(jù)同化系統(tǒng)(Global Land Data Assimilation System，GLDAS)數(shù)據(jù)進行對比分析，證明GRACE數(shù)據(jù)能較精確地反映各流域水儲量的時空分布變化。

近年來，隨著人類對水資源的高強度開發(fā)利用，流域內(nèi)水分布失衡現(xiàn)象頻現(xiàn)，誘發(fā)了一系列生態(tài)環(huán)境問題。植被作為陸地生態(tài)系統(tǒng)的重要組成部分，對全球能量平衡、生物化學循環(huán)以及水循環(huán)都起著至關(guān)重要的作用。遙感數(shù)據(jù)生成的歸一化植被指數(shù)(Normalized Difference Vegetation Index, NDVI)能很好地反映地表植被的繁茂程度，與植被密度呈正相關(guān)，即NDVI值越大，反映植被覆蓋情況越好[10]。在歸一化植被指數(shù)的研究中，早期相關(guān)文獻所用遙感數(shù)據(jù)多為空間分辨率不高(8 km)的NOAA/AVHRR數(shù)據(jù)，隨著1999年Terra衛(wèi)星的升空，其搭載的中分辨率成像光譜儀MODIS數(shù)據(jù)產(chǎn)品與NOAA/AVHRR數(shù)據(jù)產(chǎn)品相比，具有更高的靈敏度和精度，數(shù)據(jù)質(zhì)量得到明顯改善并且描述植被信息時所受干擾較小，極大提升了植被指數(shù)的監(jiān)測能力[11]。

水是植物生長發(fā)育的關(guān)鍵生態(tài)因子，區(qū)域水資源的分布決定著植被的分布，而植被通過對水分的吸收、蒸騰、滯流以及林地的滲透、貯蓄作用，調(diào)節(jié)降水、蒸發(fā)、徑流和土壤水分的增減，對區(qū)域水儲量有重大影響。近年來，許多學者應(yīng)用遙感數(shù)據(jù)研究了植被與多項氣候因子在不同時空尺度上的響應(yīng)關(guān)系[12-14]，但植被變化與水儲量變化關(guān)系的研究并不多見。

本文利用MODIS數(shù)據(jù)和GRACE重力衛(wèi)星數(shù)據(jù)，分別獲得長江流域歸一化植被指數(shù)和陸地水儲量變化，從趨勢性、相關(guān)性和時空變化特征3方面分析長江流域NDVI以及TWSA多年時空動態(tài)變化特征及二者的響應(yīng)關(guān)系。通過灰色預測方法對NDVI的變化趨勢進行定量分析，并利用長江流域NDVI和TWSA的最小二乘擬合關(guān)系預測流域未來水儲量的變化。

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

2.1 研究區(qū)域概況

長江流域幅員遼闊，有著十分優(yōu)越的自然條件，開發(fā)潛力巨大，但長期以來，由于各種自然因素和人為因素的干擾，流域生態(tài)問題日益突出。長江流域全年季風活躍，氣候條件復雜，年降水量及暴雨時空分布不均，夏季常有強降水出現(xiàn)，使得水儲量變化具有很強的季節(jié)周期，因此，長江流域水儲量變化的研究對管理中國水資源具有重要作用。

2.2 數(shù)據(jù)來源與處理

2.2.1 GRACE重力衛(wèi)星數(shù)據(jù)

GRACE數(shù)據(jù)主要由美國德克薩斯大學空間研究中心(CSR)、德國地學研究中心(GFZ)及美國噴氣推進實驗室(JPL)進行原始數(shù)據(jù)的解算并向公眾發(fā)布，其值為原值減去2004—2009年陸地水儲量平均值后的距平值。本文將3家機構(gòu)發(fā)布數(shù)據(jù)的平均值作為陸地水儲量變化(TWSA)，3份數(shù)據(jù)的空間分辨率均為1°×1°。根據(jù)流域地形圖，將其進行面積加權(quán)并計算平均數(shù)，得到2002—2015年間長江流域各月的平均水儲量變化。

2.2.2 歸一化植被指數(shù)(NDVI)

NDVI的取值范圍在-1～1之間：當-1≤NDVI<0時，表示地面被云、水、雪等對可見光高反射的物體所覆蓋；當NDVI=0時，地表為巖石或裸土；當0

本文的NDVI數(shù)據(jù)來自NASA MODIS陸地產(chǎn)品組開發(fā)的2000—2015年MODIS植被指數(shù)產(chǎn)品MOD13A2(全球空間分辨率為1 km，時間分辨率為16 d)。將長江流域的MODIS數(shù)據(jù)，采用國際通用的MVC最大化合成法，使用遙感圖像處理程序ENVI和MRT批處理工具將其處理為逐月NDVI數(shù)據(jù)。

3 數(shù)據(jù)分析

3.1 趨勢性分析

長江流域2002—2015年間的NDVI和TWSA變化如圖1所示，由圖可知二者的年際變化具有較好的一致性。選取Mann-Kendall秩次相關(guān)檢驗對長江流域2002—2015年NDVI和TWSA數(shù)據(jù)進行趨勢性檢驗，結(jié)果見表1。由表1結(jié)果可看出，NDVI數(shù)據(jù)有顯著增加趨勢(在0.05顯著性水平下)，TWSA無顯著趨勢，但二者總體趨勢都為上升狀態(tài)。

圖1 長江流域TWSA和NDVI年際變化過程Fig.1 Interannual variabilities of TWSA and NDVI inYangtze River basin

變量變化率τ顯著性水平UNDVI0.53852.6825**TWSA0.33331.5086

注：τ值為正代表上升趨勢；U值的**表示在0.05顯著性水平下存在明顯趨勢

圖2 長江流域TWSA和NDVI頻譜圖Fig.2 Power spectrum of TWSA and NDVI in Yangtze River basin

3.2 相關(guān)性分析

為直觀分析NDVI與TWSA之間的關(guān)系，利用Pearson相關(guān)系數(shù)計算二者的相關(guān)性，結(jié)果表明，二者的相關(guān)系數(shù)達到0.815(P<0.01)，由此可知，TWSA與NDVI之間有強相關(guān)性。利用快速傅里葉變換將長江流域的TWSA和NDVI數(shù)據(jù)從時域轉(zhuǎn)到頻域，如圖2所示，由圖2可知，二者之間的相關(guān)性更加顯著了，同時可以看出長江流域TWSA和NDVI在1 a左右的低頻存在明顯的同周期。

圖3表現(xiàn)了二者間的幅值平方相關(guān)系數(shù)，對于低頻部分可明顯看出各要素之間的相關(guān)性接近1，而對于高頻部分，由于高頻固有噪聲等因素的影響，其相關(guān)性較差。

圖3 長江流域TWSA和NDVI幅值平方相關(guān)圖Fig.3 Magnitude-squared coherence of TWSA and NDVI in Yangtze River basin

通過以上分析，可得出結(jié)論：長江流域GRACE陸地水儲量變化與NDVI之間具有高度的相關(guān)性，并存在一定的線性關(guān)系。

3.3 長序列時空變化分析

3.3.1 多年月均變化特征

為進一步分析陸地水儲量變化與歸一化植被指數(shù)之間的相關(guān)性，對二者2002—2015年間多年月均變化進行對比分析，如圖4所示。

圖4 長江流域TWSA和NDVI多年月均變化曲線Fig.4 Curves of monthly average changes of TWSA and NDVI in Yangtze River basin for many years

由圖4可知，植被和水儲量的年內(nèi)變化都較為明顯，二者變化趨勢基本一致。1，2月份植被處于冬季的休眠期，由多年月平均NDVI的變化曲線可知，該時期NDVI處于全年月份最低值，與之相對應(yīng)的TWSA多年月平均變化曲線也有全年1，2月份為最低值的規(guī)律。3月份植被慢慢復蘇進入返青期，但此時植被還處于休眠期末期，因此生長仍較緩慢。4，5月份流域植被進入生長期，植被NDVI有明顯增幅，且增長速度較快，由圖4能明顯看出4，5月份NDVI的增速大于TWSA的增加速度。7，8月份植被進入茂盛時期，此時大部分植被生長達到頂峰狀態(tài)，植物內(nèi)葉綠素大量累積，8月份NDVI值達到全年峰值。8月末至9月初，植被漸漸枯黃，落葉樹種開始落葉，NDVI值出現(xiàn)明顯下降趨勢。10月份，植被指數(shù)持續(xù)下降。11月末，隨著霜降及雨雪的來臨，NDVI值不斷降低直至植被停止生長進入冬季休眠期。

長江流域水資源豐富，降水主要集中于6—8月份。從GRACE曲線上可以看出，12月份至次年5月份，水儲量均處于虧損狀態(tài)；之后持續(xù)上升到8月份，水儲量達到全年峰值；9月份之后出現(xiàn)明顯下降，全年最低值出現(xiàn)在2月份，與NDVI變化過程曲線基本一致。

3.3.2 多年季變化特征

歸一化植被指數(shù)和陸地水儲量不僅存在年際變化，同時也有顯著的季節(jié)變化，對長江流域2002—2015年間NDVI和TWSA的多年季平均變化進行分析，變化過程如圖5所示。

圖5 長江流域TWSA和NDVI多年季平均變化過程Fig.5 Seasonal variation process of TWSA and NDVI in Yangtze River basin for many years

由圖5可知，NDVI和TWSA有顯著季節(jié)變化，二者變化規(guī)律一致，均為夏季>秋季>春季>冬季。春季土壤溫度和土壤水含量較低，部分耐寒植被在嚴冬惡劣的天氣條件下逐漸衰亡，只有適應(yīng)性強的地表植物存活生長；而夏季冰川融水充沛，地表徑流較為豐富，人為種植活動頻繁，使得NDVI達到最大值；秋季降水減少，農(nóng)作物進入收獲階段，耐旱耐鹽植被能夠繼續(xù)存活；11,12月份耐寒耐旱植被依然處于生長狀態(tài)，但大部分植被進入冬眠期，因此冬季植被指數(shù)明顯降低。同時陸地水儲量在不同季節(jié)也有顯著差異，受季風影響，長江流域年降水量和暴雨的時空分布很不均勻，并且由于長江是雨洪河流，大部分暴雨發(fā)生在5—10月份，因此水儲量變化表現(xiàn)出很強的季節(jié)性。夏天常出現(xiàn)強降水，水儲量達到最大值，冬季由于進入枯水期，降水較少，整個流域水儲量為虧損狀態(tài)。

3.3.3 多年時空變化特征

利用NDVI和GRACE數(shù)據(jù)，計算了長江流域2002—2015年的NDVI和TWSA多年平均空間變化，如圖6所示。

圖6 長江流域TWSA和NDVI多年平均空間變化Fig.6 Annual average spatial variations of TWSA and NDVI in Yangtze River basin

從多年空間變化分布上看，長江中游地區(qū)的NDVI均值最高，下游次之，上游最低，這與長江流域的地理環(huán)境相符。長江流域中游地區(qū)以平原為主，地勢較為平坦，土壤資源肥沃，光照條件適宜植被的生長，并且在中游地區(qū)有大量常綠闊葉林原生植被，受人類影響較少，故長江流域中游地區(qū)的NDVI值較高，GRACE所反演的中游地區(qū)平均水儲量值在全流域中最高。源區(qū)植被指數(shù)與水儲量的盈余狀態(tài)吻合度不高，這是因為受源頭地區(qū)的氣候條件影響，該地區(qū)氣溫較低，降低了植被的生長活力，并且長江流域上游的天然林受到破壞，林地面積大幅減少，導致上游地區(qū)尤其是源頭區(qū)域的植被覆蓋情況較差，因此植被指數(shù)最低。而高寒生態(tài)系統(tǒng)對長江源區(qū)水循環(huán)和水源涵養(yǎng)能力具有促進作用，近10 a來，長江源區(qū)水資源量明顯增多[16]，源區(qū)的多年水儲量處于盈余狀態(tài)。對于下游地區(qū)，近年來，下游地區(qū)的一些省市，尤其上海、浙江沿海地帶的快速城市化現(xiàn)象，對植被的生長造成了一定的影響，因此，下游植被指數(shù)較中游稍低，與之相對應(yīng)，下游的水儲量處于盈余狀態(tài)，但較之中游仍較少。

4 長江流域陸地水儲量預測

GRACE數(shù)據(jù)的發(fā)布需要一個較長的解算過程，目前發(fā)布的GRACE數(shù)據(jù)時間跨度為2002年4月份至2016年3月份，不僅時間序列較短，發(fā)布周期也較長，且最新發(fā)布的2015年數(shù)據(jù)中10，11月份的數(shù)據(jù)缺失，12月份的數(shù)據(jù)異常，對分析其未來變化具有一定的局限性。而MODIS數(shù)據(jù)產(chǎn)品每月實時發(fā)布，最新的數(shù)據(jù)已更新到2016年10月份，因此可以由2000—2015年的NDVI年植被指數(shù)，根據(jù)灰色預測模型定量探究其未來變化趨勢，并利用NDVI與GRACE之間的線性關(guān)系，通過最小二乘擬合，預測出GRACE陸地水儲量的未來變化。

4.1 植被指數(shù)預測

灰色預測具有原理簡單、樣本數(shù)量要求低、算法簡便、短期預測精度高、可進行檢驗等優(yōu)點,本文選擇灰色預測方法來預測植被指數(shù)。GM(1，1)模型作為灰色預測理論的核心，水文上已被廣泛用于預測降雨量、徑流量等，主要適用于預測時間短，波動不大的系統(tǒng)對象，可通過灰色預測建立NDVI預測模型。

4.1.1 GM(1，1)模型的建立

建立灰色模型(Grey Model，GM)使用的是微分方程動態(tài)建模方法[17]，常用的GM(1，1)模型是一個擬微分方程的動態(tài)系統(tǒng)，其建模實質(zhì)是對原始數(shù)據(jù)先進行一次累加，使生成的數(shù)據(jù)序列呈現(xiàn)一定規(guī)律，然后通過建立一階微分方程模型，求得擬合曲線，用以對系統(tǒng)進行預測。設(shè)原始數(shù)據(jù)序列(x(0)(1)，x(0)(2)，…，x(0)(n))滿足x(0)(k)≥0，k=1，2，…，n，GM(1，1)模型預測公式為

k=2，3,…，n。

(1)

式中：ɑ為發(fā)展灰數(shù)；u為內(nèi)性控制灰數(shù)。精度檢驗可利用相對誤差Δk，其表達式為

(2)

參照預測精度檢驗等級表進行結(jié)果評價，見表2。

表2 預測精度檢驗等級Table 2 Ratings of prediction accuracy

4.1.2 NDVI灰色預測

選取2000—2014年的數(shù)據(jù)資料作為預測的原始數(shù)據(jù)序列，建立灰色預測GM(1，1)模型，預測2015年的植被指數(shù)，并將模擬值、預測值與實際值作對比檢驗精度。

建立GM(1，1)模型，其表達式為

x′(0)(k)=0.506 4e0.005 9(k-1)，k=2，3,…，n。(3)

模型計算結(jié)果見表3。

表3 GM(1，1)模型模擬值和預測值與實際值的比較Table 3 Comparison between values simulated byGM (1, 1) model and measured values

注：帶*號數(shù)據(jù)為2015年預測值

進一步計算模型的模擬誤差與預測誤差，得出2001—2014年模擬值的平均相對誤差為0.015，2015年的預測值相對誤差為0.019，均通過精度檢驗。因此，可將2000—2015年的NDVI年植被指數(shù)作為原始數(shù)據(jù)序列，建立GM(1，1)模型，即

x′(0)(k)=0.505 1e0.006 4(k-1)，k=2，3，…，n。(4)

將k= 17代入式(4)，可得2016年NDVI的預測值為0.559 8。

4.2 最小二乘擬合

Pearson相關(guān)系數(shù)表明長江流域TWSA與NDVI之間高度相關(guān)，并存在明顯的線性關(guān)系，因此可利用NDVI與TWSA進行最小二乘擬合，從而通過預測NDVI的變化來推導長江流域未來陸地水儲量變化。由于TWSA的值為距平值，故將NDVI的值減去2004—2009年的平均值，得到NDVI距平值后再建立其與TWSA的線性關(guān)系,如圖7所示。利用通過灰色預測得到的長江流域2016年NDVI數(shù)值，最小二乘擬合關(guān)系式推算出2016年長江流域陸地水儲量變化值為1.839 7 cm。

圖7 TWSA與NDVI最小二乘擬合結(jié)果Fig.7 Result of least square fitting between TWSA and NDVI

4.3 陸地水儲量變化豐枯劃分

根據(jù)國家標準《水文基本術(shù)語和符號標準》(GB/T 50095—98)中的徑流豐枯劃分準則，將陸地水儲量變化的豐枯水年細分為特豐水年、偏豐水年、平水年、偏枯水年和特枯水年。通?？烧J為特豐、偏豐水年為對應(yīng)豐水年；特枯、偏枯水年對應(yīng)枯水年。

一般可認為水文數(shù)據(jù)服從P-Ⅲ型分布，本文采用傳統(tǒng)的頻率分析法來劃分陸地水儲量變化的豐枯，統(tǒng)計確定劃分長江流域年均水儲量變化的豐、枯標準(見表4)，其中P為頻率，S為陸地水儲量變化值。根據(jù)劃分表，由預測結(jié)果可知，長江流域2016年應(yīng)為偏豐水年，且數(shù)值接近特豐區(qū)間。

表4 陸地水儲量變化豐平枯劃分標準Table 4 Classification standard of TWSA

5 結(jié) 論

(1)GRACE重力衛(wèi)星數(shù)據(jù)反演的長江流域2002—2015年陸地水儲量變化和同時期歸一化植被指數(shù)(NDVI)在年內(nèi)變化及季節(jié)變化上有較高的吻合度，均呈現(xiàn)2月份達最低值，8月份達峰值，季節(jié)變化為夏季>秋季>春季>冬季的規(guī)律，二者總體均呈微弱上升趨勢，TWSA與NDVI之間有強相關(guān)性。

(2)利用2000—2015年長江流域NDVI值和灰色模型預測其2016年NDVI值，精度滿足預測模型精度要求。再通過NDVI與TWSA之間的強相關(guān)性，結(jié)合最小二乘線性擬合得出長江流域2016年的陸地水儲量變化值。

(3)根據(jù)陸地水儲量豐枯標準判定長江流域2016年為偏豐水年，且數(shù)值接近特豐區(qū)間，表明長江流域未來水儲量為豐水狀態(tài)，基于水儲量變化對極端氣候的高度敏感性，可知2016年長江流域防洪形勢嚴峻。

參考文獻：

[1] 王 文，王 鵬，崔 巍．長江流域陸地水儲量與多源水文數(shù)據(jù)對比分析[J]．水科學進展，2015，26(6)：759-768．

[2] RIEGGER J，TOURIAN M J． Characterization of Runoff-storage Relationships by Satellite Gravimetry and Remote Sensing[J]．Water Resources Research，2014，50(4)：3444-3466．

[3] 胡小工，陳劍利，周永宏，等．利用GRACE空間重力測量監(jiān)測長江流域水儲量的季節(jié)性變化[J]．中國科學D輯 地球科學，2006，36(3)：225-232．

[4] 曹艷萍，南卓銅．利用GRACE重力衛(wèi)星監(jiān)測黑河流域水儲量變化[J]．遙感技術(shù)與應(yīng)用，2011，26(6)：719-727．

[5] 翟 寧，王澤民，伍 岳，等．利用GRACE反演長江流域水儲量變化[J]．武漢大學學報(信息科學版)，2009，34(4)：436-439．

[6] 尼勝楠，陳劍利，李 進，等．利用GRACE衛(wèi)星時變重力場監(jiān)測長江、黃河流域水儲量變化[J]．大地測量與地球動力學，2014，34(4)：49-55．

[7] 羅志才，李 瓊，鐘 波．利用GRACE時變重力場反演黑河流域水儲量變化[J]．測繪學報，2012，41(5)：676-681．

[8] 任永強，潘 云，宮輝力．海河流域地下水儲量空間變化趨勢分析[J]．首都師范大學學報(自然科學版)，2014，35(2)：89-98．

[9] 許朋琨，張萬昌．GRACE 反演近年青藏高原及雅魯藏布江流域陸地水儲量變化[J]．水資源與水工程學報，2013,(1)：23-29．

[10] 李生勇，王曉卿，李 彪．基于MODIS 數(shù)據(jù)的科爾沁區(qū)植被覆蓋時空變化分析[J]．長江科學院院報，2016，33(2)：121-127．

[11] 譚德寶，劉良明，鄢俊潔，等．MODIS 數(shù)據(jù)的干旱監(jiān)測模型研究[J]．長江科學院院報，2004，21(3)：11-15．

[12] SCHMIDTH，GITELSON A．Temporal and Spatial Vegetation Cover Change in Israeli Transition Zone：AVHRR-based Assessment of Rainfall Impact[J]．International Journal of Remote Sensing，2002，21(5)：997-1010．

[13] 張井勇，董文杰，葉篤正，等．中國植被覆蓋對夏季氣候影響的新證據(jù)[J]．科學通報，2003，48(1)：91-95．

[14] 李秀花，師慶東，郭 娟，等．中國西北干旱區(qū)1981—2001年NDVI對氣候變化的響應(yīng)分析[J]．干旱區(qū)資源與環(huán)境，2009，23(2)：12-16．

[15] 張鈦仁，張佳華，申彥波，等．1981—2001年西北地區(qū)植被變化特征分析[J]．中國農(nóng)業(yè)氣象，2010，31(3)：586-590．

[16] 吳珊珊，姚治君，姜麗光，等．基于MODIS的長江源植被NPP時空變化特征及其水文效應(yīng)[J]．自然資源學報，2016，31(1)：39-51．

[17] 王志強，汪潔華，王式功，等．灰色模型在環(huán)渤海地區(qū)降水過程預報中的應(yīng)用[J]．干旱氣象，2012，30(2)：272-275.