基于GRACE和MODIS数据的长江流域陆地水储量变化分析及预测

，，，

(河海大学 水利水电学院,南京 210098)

1 研究背景

陆地水储量(Terrestrial Water Storage，TWS)是水文循环的一个重要组成部分，是地表水、地下水、土壤含水量、冰雪和生物体含水量变化的综合体现。陆地水储量变化(Terrestrial Water Storage Anomaly，TWSA)的准确估计和趋势分析，对理解全球和局域水循环过程、研究和预报气候的变化、农业生产的安排、洪涝以及其他自然灾害的防治都具有重要意义[1]。然而，由于土壤水、地下水等要素的地面观测站点稀疏，可见光-红外及微波遥感数据只能得到表层土壤含水量，利用传统方法难以对各分量进行有效观测，水储量变化研究进展缓慢。在此背景下，2002年3月份发射的GRACE(Gravity Recovery and Climate Experiment)重力卫星，通过监测全球重力场的变化，综合反映了气候波动、全球变化和人类活动对水资源储量的影响，为陆地水储量变化研究提供了新方法[2]。目前，我国对GRACE重力卫星数据反演水储量变化的流域尺度研究，已涉及到长江、黄河、海河、黑河以及雅鲁藏布江等流域[3-9]，其中部分研究通过与全球陆面数据同化系统(Global Land Data Assimilation System，GLDAS)数据进行对比分析，证明GRACE数据能较精确地反映各流域水储量的时空分布变化。

近年来，随着人类对水资源的高强度开发利用，流域内水分布失衡现象频现，诱发了一系列生态环境问题。植被作为陆地生态系统的重要组成部分，对全球能量平衡、生物化学循环以及水循环都起着至关重要的作用。遥感数据生成的归一化植被指数(Normalized Difference Vegetation Index, NDVI)能很好地反映地表植被的繁茂程度，与植被密度呈正相关，即NDVI值越大，反映植被覆盖情况越好[10]。在归一化植被指数的研究中，早期相关文献所用遥感数据多为空间分辨率不高(8 km)的NOAA/AVHRR数据，随着1999年Terra卫星的升空，其搭载的中分辨率成像光谱仪MODIS数据产品与NOAA/AVHRR数据产品相比，具有更高的灵敏度和精度，数据质量得到明显改善并且描述植被信息时所受干扰较小，极大提升了植被指数的监测能力[11]。

水是植物生长发育的关键生态因子，区域水资源的分布决定着植被的分布，而植被通过对水分的吸收、蒸腾、滞流以及林地的渗透、贮蓄作用，调节降水、蒸发、径流和土壤水分的增减，对区域水储量有重大影响。近年来，许多学者应用遥感数据研究了植被与多项气候因子在不同时空尺度上的响应关系[12-14]，但植被变化与水储量变化关系的研究并不多见。

本文利用MODIS数据和GRACE重力卫星数据，分别获得长江流域归一化植被指数和陆地水储量变化，从趋势性、相关性和时空变化特征3方面分析长江流域NDVI以及TWSA多年时空动态变化特征及二者的响应关系。通过灰色预测方法对NDVI的变化趋势进行定量分析，并利用长江流域NDVI和TWSA的最小二乘拟合关系预测流域未来水储量的变化。

2 数据与方法

2.1 研究区域概况

长江流域幅员辽阔，有着十分优越的自然条件，开发潜力巨大，但长期以来，由于各种自然因素和人为因素的干扰，流域生态问题日益突出。长江流域全年季风活跃，气候条件复杂，年降水量及暴雨时空分布不均，夏季常有强降水出现，使得水储量变化具有很强的季节周期，因此，长江流域水储量变化的研究对管理中国水资源具有重要作用。

2.2 数据来源与处理

2.2.1 GRACE重力卫星数据

GRACE数据主要由美国德克萨斯大学空间研究中心(CSR)、德国地学研究中心(GFZ)及美国喷气推进实验室(JPL)进行原始数据的解算并向公众发布，其值为原值减去2004—2009年陆地水储量平均值后的距平值。本文将3家机构发布数据的平均值作为陆地水储量变化(TWSA)，3份数据的空间分辨率均为1°×1°。根据流域地形图，将其进行面积加权并计算平均数，得到2002—2015年间长江流域各月的平均水储量变化。

2.2.2 归一化植被指数(NDVI)

NDVI的取值范围在-1～1之间：当-1≤NDVI<0时，表示地面被云、水、雪等对可见光高反射的物体所覆盖；当NDVI=0时，地表为岩石或裸土；当0

本文的NDVI数据来自NASA MODIS陆地产品组开发的2000—2015年MODIS植被指数产品MOD13A2(全球空间分辨率为1 km，时间分辨率为16 d)。将长江流域的MODIS数据，采用国际通用的MVC最大化合成法，使用遥感图像处理程序ENVI和MRT批处理工具将其处理为逐月NDVI数据。

3 数据分析

3.1 趋势性分析

长江流域2002—2015年间的NDVI和TWSA变化如图1所示，由图可知二者的年际变化具有较好的一致性。选取Mann-Kendall秩次相关检验对长江流域2002—2015年NDVI和TWSA数据进行趋势性检验，结果见表1。由表1结果可看出，NDVI数据有显著增加趋势(在0.05显著性水平下)，TWSA无显著趋势，但二者总体趋势都为上升状态。

图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 相关性分析

为直观分析NDVI与TWSA之间的关系，利用Pearson相关系数计算二者的相关性，结果表明，二者的相关系数达到0.815(P<0.01)，由此可知，TWSA与NDVI之间有强相关性。利用快速傅里叶变换将长江流域的TWSA和NDVI数据从时域转到频域，如图2所示，由图2可知，二者之间的相关性更加显著了，同时可以看出长江流域TWSA和NDVI在1 a左右的低频存在明显的同周期。

图3表现了二者间的幅值平方相关系数，对于低频部分可明显看出各要素之间的相关性接近1，而对于高频部分，由于高频固有噪声等因素的影响，其相关性较差。

图3 长江流域TWSA和NDVI幅值平方相关图Fig.3 Magnitude-squared coherence of TWSA and NDVI in Yangtze River basin

通过以上分析，可得出结论：长江流域GRACE陆地水储量变化与NDVI之间具有高度的相关性，并存在一定的线性关系。

3.3 长序列时空变化分析

3.3.1 多年月均变化特征

为进一步分析陆地水储量变化与归一化植被指数之间的相关性，对二者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可知，植被和水储量的年内变化都较为明显，二者变化趋势基本一致。1，2月份植被处于冬季的休眠期，由多年月平均NDVI的变化曲线可知，该时期NDVI处于全年月份最低值，与之相对应的TWSA多年月平均变化曲线也有全年1，2月份为最低值的规律。3月份植被慢慢复苏进入返青期，但此时植被还处于休眠期末期，因此生长仍较缓慢。4，5月份流域植被进入生长期，植被NDVI有明显增幅，且增长速度较快，由图4能明显看出4，5月份NDVI的增速大于TWSA的增加速度。7，8月份植被进入茂盛时期，此时大部分植被生长达到顶峰状态，植物内叶绿素大量累积，8月份NDVI值达到全年峰值。8月末至9月初，植被渐渐枯黄，落叶树种开始落叶，NDVI值出现明显下降趋势。10月份，植被指数持续下降。11月末，随着霜降及雨雪的来临，NDVI值不断降低直至植被停止生长进入冬季休眠期。

长江流域水资源丰富，降水主要集中于6—8月份。从GRACE曲线上可以看出，12月份至次年5月份，水储量均处于亏损状态；之后持续上升到8月份，水储量达到全年峰值；9月份之后出现明显下降，全年最低值出现在2月份，与NDVI变化过程曲线基本一致。

3.3.2 多年季变化特征

归一化植被指数和陆地水储量不仅存在年际变化，同时也有显著的季节变化，对长江流域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有显著季节变化，二者变化规律一致，均为夏季>秋季>春季>冬季。春季土壤温度和土壤水含量较低，部分耐寒植被在严冬恶劣的天气条件下逐渐衰亡，只有适应性强的地表植物存活生长；而夏季冰川融水充沛，地表径流较为丰富，人为种植活动频繁，使得NDVI达到最大值；秋季降水减少，农作物进入收获阶段，耐旱耐盐植被能够继续存活；11,12月份耐寒耐旱植被依然处于生长状态，但大部分植被进入冬眠期，因此冬季植被指数明显降低。同时陆地水储量在不同季节也有显著差异，受季风影响，长江流域年降水量和暴雨的时空分布很不均匀，并且由于长江是雨洪河流，大部分暴雨发生在5—10月份，因此水储量变化表现出很强的季节性。夏天常出现强降水，水储量达到最大值，冬季由于进入枯水期，降水较少，整个流域水储量为亏损状态。

3.3.3 多年时空变化特征

利用NDVI和GRACE数据，计算了长江流域2002—2015年的NDVI和TWSA多年平均空间变化，如图6所示。

图6 长江流域TWSA和NDVI多年平均空间变化Fig.6 Annual average spatial variations of TWSA and NDVI in Yangtze River basin

从多年空间变化分布上看，长江中游地区的NDVI均值最高，下游次之，上游最低，这与长江流域的地理环境相符。长江流域中游地区以平原为主，地势较为平坦，土壤资源肥沃，光照条件适宜植被的生长，并且在中游地区有大量常绿阔叶林原生植被，受人类影响较少，故长江流域中游地区的NDVI值较高，GRACE所反演的中游地区平均水储量值在全流域中最高。源区植被指数与水储量的盈余状态吻合度不高，这是因为受源头地区的气候条件影响，该地区气温较低，降低了植被的生长活力，并且长江流域上游的天然林受到破坏，林地面积大幅减少，导致上游地区尤其是源头区域的植被覆盖情况较差，因此植被指数最低。而高寒生态系统对长江源区水循环和水源涵养能力具有促进作用，近10 a来，长江源区水资源量明显增多[16]，源区的多年水储量处于盈余状态。对于下游地区，近年来，下游地区的一些省市，尤其上海、浙江沿海地带的快速城市化现象，对植被的生长造成了一定的影响，因此，下游植被指数较中游稍低，与之相对应，下游的水储量处于盈余状态，但较之中游仍较少。

4 长江流域陆地水储量预测

GRACE数据的发布需要一个较长的解算过程，目前发布的GRACE数据时间跨度为2002年4月份至2016年3月份，不仅时间序列较短，发布周期也较长，且最新发布的2015年数据中10，11月份的数据缺失，12月份的数据异常，对分析其未来变化具有一定的局限性。而MODIS数据产品每月实时发布，最新的数据已更新到2016年10月份，因此可以由2000—2015年的NDVI年植被指数，根据灰色预测模型定量探究其未来变化趋势，并利用NDVI与GRACE之间的线性关系，通过最小二乘拟合，预测出GRACE陆地水储量的未来变化。

4.1 植被指数预测

灰色预测具有原理简单、样本数量要求低、算法简便、短期预测精度高、可进行检验等优点,本文选择灰色预测方法来预测植被指数。GM(1，1)模型作为灰色预测理论的核心，水文上已被广泛用于预测降雨量、径流量等，主要适用于预测时间短，波动不大的系统对象，可通过灰色预测建立NDVI预测模型。

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

建立灰色模型(Grey Model，GM)使用的是微分方程动态建模方法[17]，常用的GM(1，1)模型是一个拟微分方程的动态系统，其建模实质是对原始数据先进行一次累加，使生成的数据序列呈现一定规律，然后通过建立一阶微分方程模型，求得拟合曲线，用以对系统进行预测。设原始数据序列(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)

式中：ɑ为发展灰数；u为内性控制灰数。精度检验可利用相对误差Δk，其表达式为

(2)

参照预测精度检验等级表进行结果评价，见表2。

表2 预测精度检验等级Table 2 Ratings of prediction accuracy

4.1.2 NDVI灰色预测

选取2000—2014年的数据资料作为预测的原始数据序列，建立灰色预测GM(1，1)模型，预测2015年的植被指数，并将模拟值、预测值与实际值作对比检验精度。

建立GM(1，1)模型，其表达式为

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

模型计算结果见表3。

表3 GM(1，1)模型模拟值和预测值与实际值的比较Table 3 Comparison between values simulated byGM (1, 1) model and measured values

注：带*号数据为2015年预测值

进一步计算模型的模拟误差与预测误差，得出2001—2014年模拟值的平均相对误差为0.015，2015年的预测值相对误差为0.019，均通过精度检验。因此，可将2000—2015年的NDVI年植被指数作为原始数据序列，建立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相关系数表明长江流域TWSA与NDVI之间高度相关，并存在明显的线性关系，因此可利用NDVI与TWSA进行最小二乘拟合，从而通过预测NDVI的变化来推导长江流域未来陆地水储量变化。由于TWSA的值为距平值，故将NDVI的值减去2004—2009年的平均值，得到NDVI距平值后再建立其与TWSA的线性关系,如图7所示。利用通过灰色预测得到的长江流域2016年NDVI数值，最小二乘拟合关系式推算出2016年长江流域陆地水储量变化值为1.839 7 cm。

图7 TWSA与NDVI最小二乘拟合结果Fig.7 Result of least square fitting between TWSA and NDVI

4.3 陆地水储量变化丰枯划分

根据国家标准《水文基本术语和符号标准》(GB/T 50095—98)中的径流丰枯划分准则，将陆地水储量变化的丰枯水年细分为特丰水年、偏丰水年、平水年、偏枯水年和特枯水年。通常可认为特丰、偏丰水年为对应丰水年；特枯、偏枯水年对应枯水年。

一般可认为水文数据服从P-Ⅲ型分布，本文采用传统的频率分析法来划分陆地水储量变化的丰枯，统计确定划分长江流域年均水储量变化的丰、枯标准(见表4)，其中P为频率，S为陆地水储量变化值。根据划分表，由预测结果可知，长江流域2016年应为偏丰水年，且数值接近特丰区间。

表4 陆地水储量变化丰平枯划分标准Table 4 Classification standard of TWSA

5 结 论

(1)GRACE重力卫星数据反演的长江流域2002—2015年陆地水储量变化和同时期归一化植被指数(NDVI)在年内变化及季节变化上有较高的吻合度，均呈现2月份达最低值，8月份达峰值，季节变化为夏季>秋季>春季>冬季的规律，二者总体均呈微弱上升趋势，TWSA与NDVI之间有强相关性。

(2)利用2000—2015年长江流域NDVI值和灰色模型预测其2016年NDVI值，精度满足预测模型精度要求。再通过NDVI与TWSA之间的强相关性，结合最小二乘线性拟合得出长江流域2016年的陆地水储量变化值。

(3)根据陆地水储量丰枯标准判定长江流域2016年为偏丰水年，且数值接近特丰区间，表明长江流域未来水储量为丰水状态，基于水储量变化对极端气候的高度敏感性，可知2016年长江流域防洪形势严峻。

参考文献：

[1] 王 文，王 鹏，崔 巍．长江流域陆地水储量与多源水文数据对比分析[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空间重力测量监测长江流域水储量的季节性变化[J]．中国科学D辑 地球科学，2006，36(3)：225-232．

[4] 曹艳萍，南卓铜．利用GRACE重力卫星监测黑河流域水储量变化[J]．遥感技术与应用，2011，26(6)：719-727．

[5] 翟 宁，王泽民，伍 岳，等．利用GRACE反演长江流域水储量变化[J]．武汉大学学报(信息科学版)，2009，34(4)：436-439．

[6] 尼胜楠，陈剑利，李 进，等．利用GRACE卫星时变重力场监测长江、黄河流域水储量变化[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 数据的科尔沁区植被覆盖时空变化分析[J]．长江科学院院报，2016，33(2)：121-127．

[11] 谭德宝，刘良明，鄢俊洁，等．MODIS 数据的干旱监测模型研究[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]．科学通报，2003，48(1)：91-95．

[14] 李秀花，师庆东，郭 娟，等．中国西北干旱区1981—2001年NDVI对气候变化的响应分析[J]．干旱区资源与环境，2009，23(2)：12-16．

[15] 张钛仁，张佳华，申彦波，等．1981—2001年西北地区植被变化特征分析[J]．中国农业气象，2010，31(3)：586-590．

[16] 吴珊珊，姚治君，姜丽光，等．基于MODIS的长江源植被NPP时空变化特征及其水文效应[J]．自然资源学报，2016，31(1)：39-51．

[17] 王志强，汪洁华，王式功，等．灰色模型在环渤海地区降水过程预报中的应用[J]．干旱气象，2012，30(2)：272-275.