基于GRACE陆地水储量降尺度的塔里木河流域干旱特征及驱动因子分析

魏光辉，杨 鹏，周海鹰，夏 军，陈 杰，桂东伟，巴音达拉

(1.新疆塔里木河流域管理局，新疆 库尔勒 841000；2.中国地质大学(武汉)，地理与信息工程学院，武汉 430074；3.武汉大学 水资源与水电工程科学国家重点实验室，武汉 430072；4.中国科学院地理科学与资源研究所，北京 100101;5.中国科学院新疆生态与地理研究所，乌鲁木齐 830011；6.新疆农业大学水利与土木工程学院，乌鲁木齐 830052)

0 引 言

塔里木河流域地处亚欧大陆腹地，远离海洋，气候干旱，是中国最大的内陆河流域，绿洲灌溉农业为其主要的农业发展模式，由于水资源的短缺，其间存在土壤盐渍化、水资源流失、生态系统退化和破坏等严重的环境问题[1-3]。在气候变化的背景下，水文气象极端事件更加频繁发生，对于靠冰雪融水补给为主的流域而言，非常不利于水资源的开发利用[3]，因此，气候变化也加剧了塔里木河流域水安全问题[4]。虽然塔里木河流域降水和出山口径流在近期有增加趋势，但是从上游向中下游输送的径流却明显减少[3]。因而，分析塔里木河流域水资源量变化十分重要。2013年9月“一带一路”经济发展战略的提出，以实现跨越亚非欧大陆国家和地区的区域经济一体化为目标，推进国际合作，打造陆海内外联动、东西双向开放的全面开放新格局，也为我国西北干旱区经济发展提供契机[5]。因此，分析和研究塔里木河流域水资源的变化对于该区生态环境保护和社会经济的发展具有至关重要的作用。

陆地水储量是全球水、能量以及生物化学循环的重要组成部分[6]，它主要包括地表水、土壤水和地下水，其中地表水包括河川径流湖泊水、雪水当量以及地表植被冠层水，也是水资源量表示的重要方法[7]。研究陆地水储量的时空变化对于有效管理水资源和可持续发展决策具有十分重要的意义。然而，由于地形地貌、海拔高度、地质结构的复杂，加上技术的局限，人工观测陆地水储量的变化存在较多困难，尤其是大尺度区域的陆地水储量的观测更是很难实施。然而， 2002年重力卫星GRACE的发射为大尺度陆地水储量的观测提供了重要工具和手段，且已经被广泛地应用到水文与水资源领域的研究中。例如：Velicogna和Wahr[8,9]基于重力卫星GRACE数据研究发现格陵兰岛南部地区的冰盖在加速减少；Longuevergne等[10]和Voss等[11]利用重力卫星GRACE数据分别研究了中东地表水和地下水的变化；周旭华等[12]基于重力卫星GRACE数据反演了全球陆地水储量的变化，并结合全球陆地同化系统GLDAS模型数据进行了对比分析和研究；钟敏等[13]利用GRACE数据发现中国华北平原、三峡库区、青藏高原、新疆阿尔金山和青川甘交界处陆地水储量变化幅度较大，且京津冀地区、青藏高原地区陆地水储量呈现明显的减少趋势。但是，重力卫星GRACE固有的空间分辨率较低及时间序列较短的特点严重限制了其应用的推广。尽管Long等[14]和Zhang等[15]应用人工智能方法结合水文气象数据对陆地水储量进行重建，延长了陆地水储量的时间序列，其较低的空间分辨率依然未得到有效解决。重力卫星GRACE陆地水储量空间分辨率的提高，需要通过科学合理的数学物理方法对其进行降尺度处理，即通过科学合理的数学物理方法在不影响重力卫星GRACE陆地水储量信息的前提下，提高其空间分辨率，以便其在中小尺度流域的应用与推广。除此之外，在水文领域重力卫星GRACE大多数情况都被用作分析陆地水储量的时空变化，其水文研究应用的潜力还亟待提高和加强。

为了研究流域尺度陆地水储量的时空变化及其影响下的流域干旱特征，本文尝试利用重力卫星GRACE数据和全球陆地同化系统模型数据GLDAS，结合统计学方法，分别以数据格点和等效水高为权重，对塔里木河流域的陆地水储量进行空间降尺度，分析塔里木河流域陆地水储量的时空变化，随后进行降尺度后的陆地水储量分析塔里木河流域的干旱特征分析，简要分析其驱动因素。

1 研究区概况

塔里木河流域是中国最大的内流河流域，其面积达102 万km2，坐落于73°E～93°E，35°N～45°N之间(如图1)，属于典型的温带大陆性气候，因此流域内降水极少，蒸散发却极其强烈。另外流域内温度变化幅度较大，常年冬季温度跨度在-20和10 ℃之间，夏季温度跨度在20和30 ℃之间[16]。就流域内的降水的空间分布而言也存在较大差异，例如：流域内山区降水较多，而平原区降水却在较少[17-19]。另外，流域内的年温度也表现出较大的年际差异[3]。

随着20世纪70年代以来气温的上升[20]，流域内固态降水减少且冰雪融水增加。然而，对于干旱区而言山区降水及来自于山区的冰雪融水在水资源系统和水文过程起着十分重要的作用[21]。塔里木河流域的年均径流量约为3.989 × 1010m3, 其中来自山区的冰雪融水大约占48.2%[22]。作为流域内的主要源流，阿克苏河、和田河及叶尔羌河的径流量分别占总径流量的73.2%，23.2%和3.6%[23]。同时，以灌溉农业为主要农业发展方式的塔里木河流域的社会经济的可持续发展深受流域内水资源波动和供给能力的影响[23]。

2 数据与方法

2.1 数 据

(1)GRACE数据。重力卫星 GRACE数据源自http:∥grace.jpl.nasa.gov。本文所涉及的陆地水储量的等效水高是从第二版重力卫星球谐系数RL05转换得到。同时，在等效水高转换过程中还有一些其他处理细节，如用卫星激光测量数据(SLR)替代C20项[24]，重新估算1阶系数[25]，冰川后回弹系数矫正[26]等。另外， 考虑到重力卫星GRACE球谐系数数据在2004年1月至2009年12月之间不缺失[27]，因此，为了得到球谐系数的距平值，本文从原始的球谐系数中去掉了该时间段的平均值。尽管遗落误差、测量误差等误差存在于重力卫星GRACE的数据中，然而Chen等[25]却证明了GRACE数据具有一定的适用性和可靠性。

(2)GLDAS数据。全球陆地同化系统模型(GLDAS)数据旨在精确反演陆地表面状况[28]。同时，在该系统由3个陆面模型(Noah, Mosaic, Community Land Model-CLM)和一个水文模型(Variable Infiltration Capacity-VIC)组成。本文所用的数据来自于网站(http://mirador.gsfc.nasa.gov/)，文中所涉及的土壤水来自于Noah模型中的时间分辨率为月，空间分辨率为1°×1°和0.25°×0.25°两种数据。该数据被广泛的用到全球及各个流域的水文科学的研究领域中，具有较强的可靠性和普适性[29]。为了能够与GRACE数据进行对比，也将2004年1月至2009年12月的平均值从原始时间序列中扣除。

(3)塔里木河流域水资源数据。本文采用的塔里木河流域水资源数据主要来自于由塔里木河流域管理局编写的《塔里木河流域水资源公报》，其中主要涉及塔里木河流域地下水资源量数据、塔里木河流域供水量数据以及塔里木河流域用水量数据，该数据由塔里木河流域管理局出版和提供，已被生态和水文领域中的很多专家和课题组应用，具有可靠性。

(4)夏季0 ℃层高度数据及MODIS积雪产品数据。夏季(6-8月)0 ℃层高度主要基于库车、喀什和和田3个探空站数据，来源于怀俄明大学气象数据中心(http：// www. weather. uwyo. edu/ upperair/ sounding.html)，可免费获取，时间分辨率为天，另外，本文是基于MODIS积雪产品提取的塔里木河流域夏季的积雪覆盖率，该数据时间分辨率为月，空间分辨率为500 m，并且该数据可以从NASA相关网站获取，具体为http：// modis-land. gsfc. nasa. gov/ snow. html。国内外许多专家和学者基于该数据对塔里木河流域进行过气候变化的分析和研究[30]。

2.2 方 法

2.2.1 陆地水储量估算方法

本文根据Warh等[9]提出的方法，将重力卫星GRACE提供的重力场球谐系数数据转化为等效水高数据，等效水高的具体计算方法如公式(1)：

ΔSnmsin(mφ)]Pnmsinθ

(1)

式中：H为陆地水储量的等效水高，mm；θ为纬度；φ为经度；α为赤道半径，km；ρave为地球密度，kg/m3；kn为勒夫数(常量)；Cnm和Snm分别为重力场球谐系数的余弦和正弦项(阶)；Pnmsinθ为第n阶、第m次标准化后的勒让德函数，其最高展开阶数是60阶。

但是，通常在计算陆地水储量的同时还得考虑重力卫星GRACE中的误差(遗漏误差和测量误差)和尺度因子，其中遗漏误差可以通过经滤波处理前后的月均GLDAS-Noah陆地水储量的均方根误差(RMSD)获得[31]，测量误差可由研究区域与同纬度海洋区域的重力差值对应的等效水高获得[31]，因此，经误差和尺度因子校正后的等效水高为：

Hc=[H-(L+M)]S

(2)

式中：Hc为校正后的等效水高，mm；H为经球谐系数直接求得的等效水高，mm；L和M分别为遗漏误差和测量误差，mm；S为尺度因子(标量)，可以从JPL数据分发网站获得(https:∥grace.jpl.nasa.gov/data/get-data/jpl_global_mascons/)。

值得注意的是重力卫星GARCE三级产品(即GRACE-Level 3)已经发布，它充分的考虑了GRACE数据处理成陆地水储量过程中会出现的各种误差和尺度因子，能很好地反映全球各区域陆地水储量的变化[32],可以从网上免费获得(http:∥www2.csr.utexas.edu/grace/RL05_mascons.html)。

2.2.2 陆地水储量降尺度方法

全球陆面模式的等效水柱包含地表水和土壤水，然而这些陆面模式并未考虑地下水通量的变化，考虑到GRACE陆地水储量的实际空间分辨率约为4°，首先对GRACE陆地水储量数据进行升尺度处理，将空间分辨率1°×1°的GRACE陆地水储量升尺度至4°×4°陆地水储量，并且将0.25°的全球四个陆面模式的等效水柱升尺度至4°×4°。然后根据Wan等[20]提出的方法，以数据格点间的面积为权重因子，对GRACE陆地水储量进行降尺度处理，具体方法如下：

sM=∑(Si×ai)/∑ai=∑(Si×ai)/∑A

(3)

式中：SM是4°×4°的陆面模式的陆地水储量距平值的等效水高，mm；ai是0.25°×0.25°数据格点i的面积，m2；A是包含0.25°×0.25°数据格点i的4°×4°数据格点的总面积，m2。

假定GRACE陆地水储量(mm)是可靠的，因此空间分辨率为4°×4°陆面模式水储量距平值(SM)与4°×4°GRACE陆地水储量距平值(SG)之间的差异可以认为是二者之间的偏差Bias(B)：

B=SM-SG

(4)

因此，基于格点等效水柱权重，进一步将4°×4°陆面模式水储量距平值(SM)与4°×4°GRACE陆地水储量距平值(SG)之间的总水柱差异(B×A)分解到0.25°×0.25°数据格点上：

bi={B×A[(Soi×ai)/∑(Soi×ai)]}/ai=

(B×A×Soi)/∑(Soi×ai)

(5)

式中：bi为空间分辨率0.25°×0.25°等效水高的偏差，mm；Soi为0.25°×0.25°求距平值之前的等效水高，mm。

(6)

2.2.3 基于陆地水储量的水文干旱计算方法

首先，本文将陆地水储量的距平值转化为陆地水储量变化，然后，根据陆地水储量变化与多年同期陆地水储量变化平均值的差值确定干旱期，具体公式表示如下：

TWSCj=TWSAj-TWSAj-1

(7)

式中：TWSCj为陆地水储量的变化，mm；TWSAj表示整个时间序列中第j个月的陆地水储量的距平值，mm；TWSAj-1表示整个时间序列中第j-1个月的陆地水储量距平值，mm。

(8)

随后，根据游程定理确定陆地水储量干旱特征(干旱月数、干旱频率、干旱次数及干旱强度)，TWSCA为负值确定为干旱发生，干旱月数即TWSCA为负值时的月份数量累加[公式(9)]，随后，干旱频率为干旱月份数量占总时间段的百分比[公式(10)]，干旱次数为所有干旱发生事件次数的累加[公式(11)]，然后，干旱强度为所有干旱事件中TWSCA累加的相反数[公式(12)]，具体表示为如下公式：

Dmonth=∑monthTWSCA<0

(9)

Dfrequency=Dmonth/Tmonth×100%

(10)

DTimes=∑eventTWSCA<0

(11)

Dintensity=-∑TWSCA, (TWSCA<0)

(12)

式中：Dmonth为干旱月数，月；Dfrequency为干旱频率，%；DTimes为干旱次数，次；Dintensity为干旱强度，标量；monthTWSCA<0为干旱发生时的月份；Tmonth为总月数；eventTWSCA<0为干旱事件。

3 结果与讨论

3.1 流域陆地水储量变化

根据GRACE数据反演出塔里木河流域的陆地水储量，研究发现流域北部为亏损状态，南部为盈余状态，亏损状态最大达到-15 mm，盈余状态最大达到10 mm(图2)，这与Yang等[4]的研究结果较为一致。2002-2015年间整个时段内塔里木河流域陆地水储量呈下降趋势，下降速度达-1.6±1.1 mm/a，陆地水储量最小值出现在2008年年底和2009年年初，最大值出现在2005年和2010年，在2005、2006、2010和2011年间陆地水储量呈现盈余状态，然而在2007、2008、2009和2014年间陆地水储量呈现亏损状态。另外，塔里木河流域内的GLDAS和GRACE反演的土壤水和陆地水储量变化趋势、相位和振幅较为一致，线性相关关系R2达到0.35。同时，研究发现由CSR、GFZ和JPL发布的球谐系数反演得到的陆地水储量变化较为一致。

图2 塔里木河流域陆地水储量的时空分布

对塔里木河流域内的陆地水储量季节变化研究发现(如图3)，流域季节陆地水储量分别为春季>夏季>冬季>秋季。对其空间分布研究发现，春季流域大多数地区陆地水储量呈盈余状态，只有在流域的东北部区域陆地水储量为亏损状态[图3(a)]；夏季流域北部天山地区陆地水储量为亏损状态，流域南部的昆仑山地区陆地水储量呈现盈余状态[图3(b)]；秋季流域内陆地水储量几乎全大部分区域都为出亏损状态，而在流域东南部的局部区域陆地水储量为盈余状态[图3(c)]；冬季流域东北部和西南部陆地水储量呈现亏损状态，流域东南部陆地水储量呈现盈余状态[图3(d)]。

图3 塔里木河流域的GRCAE陆地水储量的季节空间分布

3.2 流域陆地水储量变化的降尺度分析

由于重力卫星GRACE反演得到的水储量的空间分辨存在明显的局限性，且考虑到土壤水是陆地水储量的重要组成部分，因此根据公式(1)～(6)及结合全球陆地水储量GLDAS数据对重力卫星GRACE陆地水储量进行降尺度分析(如图4)，降尺度到0.25°×0.25°的空间分辨率。研究发现，降尺度前后的陆地水储量在空间分布上存在一定的相似性，即在塔里木河流域北部陆地水储量呈现亏损状态，流域南部陆地水储量呈现盈余状态，这与Chen 等[33]及Deng和Chen[27]的研究结果一致，与此同时，降尺度前后的陆地水储量在时间分布上较为一致，其振幅和相位也存在一致性，两者的线性相关性R2达0.98，相关关系达到了显著水平(P<0.001)。

图4 塔里木河流域的陆地水储量降尺度分析

为了进一步分析降尺度之后的陆地水储量的空间分布情况，对其进行季节变化分析(如图5)。研究发现，春季流域大多数地区陆地水储量呈盈余状态，只有在流域的北部和西部的部分区域陆地水储量为亏损状态[图5(a)]；夏季北部天山地区陆地水储量呈现亏损状态，流域南部的昆仑山地区陆地水储量呈现盈余状态[图5(b)]；秋季流域内陆地水储量几乎全大部分区域都呈现出亏损状态，而在流域东南部的局部区域陆地水储量呈现盈余状态[图5(c)]；冬季流域东北部和西南部陆地水储量呈现亏损状态，流域东南部陆地水储量呈现盈余状态[图5(d)]。其空间分布与降尺度之前的陆地水储量的季节空间分布较为一致。

图5 塔里木河流域的降尺度陆地水储量的季节空间分布

3.3 基于陆地水储量的干旱分析

为了更科学的探究塔里木河流域的干旱情况，本文基于降尺度的陆地水储量来分析塔里木河流域的干旱特征(如图6)。基于公式(7)～(12)研究发现2002年4月至2015年12月间干旱总时长最长及干旱频率最高的区域主要集中在天山南坡的东段、昆仑山北坡的东段，最长干旱时长可达112个月及相对应的干旱频率为68.2%，相比较而言，期间干旱总时长较短及干旱频率较低的区域主要集中在和田河流域的中下游，最短干旱时长为56个月及对应的干旱频率为34.1%；从干旱次数方面来看，天山南坡及昆仑山北部的部分区域干旱次数发生的较少，期间干旱次数最小为2次，开都-孔雀河流域发生干旱次数较多，最多达28次；从干旱强度方面来看，天山南坡及和田河流域上游区域干旱强度较大，最大可达5 383 mm，塔里木盆地中部干旱强度最小，最小几乎为零，这可能与高山冰雪资源消融和人类水资源开发利用有关。

图6 基于陆地水储量的塔里木河流域干旱情况分析

3.4 流域陆地水储量变化的驱动因素分析

为了研究塔里木河流域陆地水储量驱动因素的变化，本文选取了与陆地水储量相关的水文变量进行研究(即MODIS夏季积雪覆盖率、夏季0 ℃层高度及水资源量)(如图7和表1)，研究发现2003-2012年间塔里木河流域年平均夏季积雪覆盖率和0 ℃层高度分别为3.12%和4 723.7 m，二者大致呈现出一个负相关的关系，如：2003-2006年间塔里木河流域夏季积雪覆盖率较高，但是夏季0 ℃层高度却较低，2006-2008年间塔里木河流域夏季0 ℃层高度升高，相应的积雪覆盖率反而降低，直到2009年夏季0 ℃层高度达到最低值，并且出现了较低的积雪覆盖率，势必会影响整个流域的水储量变化，直到2010年时夏季0 ℃层高度达到次高值，积雪覆盖率也达到了较高值，这在一定程度上可以解释Wang等[4]发现的塔里木河流域在2009年出现极端干旱和在2010年出现极端湿润的情况，这与Chen等[33]杨鹏等[30]的研究结果一致。

图7 2003-2012年塔里木河流域夏季积雪覆盖率及夏季0 ℃层高度的变化

表1 塔里木河流域水资源相关变量距平值 亿m3

另外，结合GLDAS的陆地水储量分量研究发现，塔里木河流域冠层水和雪水当量呈现不显著的降低趋势，与陆地水储量的变化趋势较为一致，而土壤水却呈上升趋势，因此，说明该区域的地下水呈现明显的下降趋势(如图8)。同时，为了验证该结论本研究结合《塔里木河流域水资源公报数据》进行了分析和研究，研究发现2004、2007、2008、2009、2011和2012年间流域水资源量盈余为负距平状态，2002、2003、2005、2006、2010和2013年间流域水资源量盈余呈现正距平状态。另外，2004、2005、2009、2010、2011和2013年间流域地下水资源量为负距平状态，2002、2003、2006、2007、2008和2012年间流域地下水资源量呈现正距平状态。除此之外，流域水资源盈余量的距平值和地下水资源量呈下降趋势，上述水资源相关变量基本与流域内的陆地水储量变化较为一致，因此，塔里木河流域陆地水储量的变化主要是因为流域内地下水变化引起的，这一研究结果与Yang等[4]的发现一致。关于陆地水储量时空变化研究对于流域水资源管理和水旱灾害预报具有十分重要的意义。

图8 2002-2015年塔里木河流域GLDAS陆地水储量分量的时空分布

4 结 语

本文利用重力卫星GRACE数据、全球陆地同化系统模型GLDAS数据和塔里木河水资源数据，对塔里木河流域陆地水储量变化、陆地水储量变化的降尺度、塔里木河流域干旱情况及驱动因素进行了研究分析，主要结论如下。

(1)根据GRACE数据反演出塔里木河流域的陆地水储量，研究发现流域北部为亏损状态，南部为盈余状态，亏损状态最大达到-15 mm，盈余状态最大达到10 mm。2002-2015年间整个时段内塔里木河流域陆地水储量呈下降趋势，下降速度达-1.6±1.1 mm/a。

(2)结合全球陆地同化模型数据和统计方法可以很好地对重力卫星GRACE反演的陆地水储量进行降尺度分析，这使得重力卫星GRACE在反演中小尺度流域陆地水储量变化上的应用具有十分重要的意义。

(3)尽管在气候变化背景下塔里木河流域山区来水量有增加趋势，但是塔里木河流域水储量赤字情况依然很严峻，干旱问题突出，水资源科学合理地开发和利用十分有必要。

(4)塔里木河流域陆地水储量受夏季0 ℃层高度及夏季积雪覆盖率的影响较为明显，另外，流域内水资源盈余量的距平值和地下水资源距平值的变化与塔里木河流域内的陆地水储量变化较为一致，因此，流域地下水资源的变化是流域内陆地水储量的变化的主要因素。