利用GRACE数据联合新型尺度因子校正法提高陆地水储量变化准确性

杨帅， 郑伟,3*， 尹文杰， 刘杰

1 河南理工大学测绘与国土信息工程学院， 焦作河南 454003 2 中国空间技术研究院钱学森空间技术实验室， 北京 100094 3 辽宁工程技术大学测绘与地理科学学院， 阜新辽宁 123000

0 引言

自2002年重力恢复与气候实验(Gravity Recovery and Climate Experiment，GRACE)卫星发射以来，GRACE重力卫星提供了高精度观测地球时变重力场信号的新手段(Tapley et al., 2004).截止到目前，GRACE重力卫星已被广泛应用于陆地水储量(Terrestrial Water Storage，TWS)变化(Rodell et al., 2009; Seyoum and Milewski, 2016)、冰川变化及其质量平衡(Velicogna and Wahr, 2006a,b; Mohajerani et al., 2018)、全球平均海洋质量变化(Chen et al., 2019)、地震(Han et al., 2006; Chen et al., 2007; Zhou et al., 2018)等多领域研究.GRACE重力卫星在2017年10月停止运行，随着GRACE Follow-on重力卫星的成功发射，该卫星反演的地球重力场模型继续在地球科学领域发挥着极其重要的作用(郑伟等, 2014).

基于Wahr等(1998)的理论，地球表面的质量变化可以利用GRACE球谐系数解反演得到.然而，球谐系数中既含有随系数阶数而增加的随机误差，又包含与系统相关的系统误差.因此，许多研究者相继提出了空间平滑滤波(如高斯平滑滤波(Wahr et al., 1998)、维纳滤波(Sasgen et al., 2007)等)和去条带滤波(Swenson和Wahr(2006)的去条带方法、Chen等(2009)的P4M6方法等).这些滤波方法在消除误差影响的同时也会改变真实的物理信号.Chen等(2015)通过实验模拟得出截断及平滑300 km产生的泄露误差在南极洲可达80.95 Gt·a-1(实测值-150.00 Gt·a-1，模拟值-69.05 Gt·a-1)；Jin和Zou(2015)研究表明格陵兰岛冰盖质量损失的信号中，42.4%是由泄漏误差造成，对泄漏误差进行校正后，格陵兰岛的冰川融化速度从-15.58±0.93 Gt·a-2增至-26.19±1.67 Gt·a-2.因此，泄露误差是GRACE反演质量变化中的主要限制因素，校正泄露误差对提高反演结果的准确性至关重要.

目前，GRACE数据泄露误差常用的恢复方法包括：(1)尺度因子法(Scaling Factor，SF).Long等(2015a)利用全球水文模型PCR-GLOBWB得到的尺度因子用于恢复长江流域水储量变化的信号损失，该方法对于提高研究区水储量变化的准确性具有重要价值；Huang等(2019)基于全球不同水文模型对尺度因子泄露误差校正法进行敏感性分析，对中国西南部水储量变化进行了评估.尺度因子法依赖先验水文模型，虽然水文模型已被广泛应用于GRACE数据泄露误差校正的先验信息(Long et al., 2014; Chen et al., 2017; Pan et al., 2017)，但是水文模型存在一些缺点，例如，地表水和地下水表示成分不完整以及没有考虑人类活动影响等.(2)迭代恢复法(Forward Modeling，FM).该方法不需要借助先验模型，主要受GRACE原始观测值的约束.Chen等(2014)利用迭代恢复法恢复了印度西北长期地下水变化，新估算的水储量变化趋势反映迭代恢复法可以有效减少GRACE估计中的泄露误差；吴云龙等(2015)基于迭代恢复法可有效恢复黑河流域水储量变化初始信号，泄露误差减少，提高水储量变化的空间分辨效果.迭代恢复法计算过程相对复杂，在部分区域可能存在信号恢复不全或过量恢复现象.与尺度因子法相比，迭代恢复法是一种完全不同于尺度因子法的泄露误差恢复方法，恢复结果可与尺度因子法校正结果进行验证对比分析.

不同于前人的已用研究，本文首次利用模型与GRACE数据的均方根误差和相关性加权尺度因子，构建新型尺度因子校正法(New Scale Factor Correction Method，NSFCM).该方法充分考虑多种水文模型与GRACE之间的一致性和准确性信息，相比于传统尺度因子法和迭代恢复法，可以提高GRACE陆地水储量变化的准确性.

1 新型尺度因子校正法的构建

传统尺度因子法使用单一水文模型计算的尺度因子依赖水文模型结果.然而，不同的水文模型由于模型的结果以及驱动数据的差异，导致模型输出结果间也存在一定的差异(Long et al., 2013).因此，本文根据多种水文模型计算的陆地水储量与CSR-SHc数据之间的一致性和准确性信息赋予权重，综合多种尺度因子模型提高GRACE陆地水储量变化的准确性.

新型尺度因子校正法主要分为以下几个步骤：(1)水文格网数据转化为球谐数据；(2)高斯平滑滤波和去条带滤波；(3)球谐系数转化为格网数据；(4)基于最小二乘原理计算尺度因子；(5)计算水文模型的陆地水储量与CSR-SHc数据均方根误差和相关性；(6)基于模型权重加权尺度因子.具体方法如下：

(1)水文格网数据转化为球谐数据

根据Wahr等(1998)的理论，可将水文模型陆地水储量转为与之对应的球谐系数：

×sinθdθdλ，

(1)

(2)高斯平滑滤波和去条带滤波

为了保持与GRACE数据处理相同的滤波，本文使用Swenson去条带滤波和300 km的高斯平滑滤波.去条带滤波通常采用以l阶为中心，宽为w的二次多项式滤波窗口，Swenson和Wahr(2006)采用的平滑窗口宽度w和球谐系数次数m的关系为

(2)

其中，max表示取两者的最大值，A=30和K=10表示选择的经验参数.

高斯滤波是指对不同阶次的球谐系数赋予不同权重，从而减小GRACE数据高阶项误差的方法.Jekeli(1981)提出用递归关系计算权重：

(3)

其中，r表示高斯平滑半径，a表示地球平均半径，e表示自然常数，Wl表示权重.

(3)球谐系数转化为格网数据

根据Wahr等(1998)的理论，可将滤波后的球谐系数转为格网数据：

(4)

(4)基于最小二乘原理计算尺度因子

尺度因子是基于水文模型陆地水储量模拟滤波过程造成的信号损失，并以最小二乘减小原始模型时间序列与滤波后模型时间序列之间差异的方法.对于任意格网点，在t时间序列对应的未滤波的陆地水储量为Δhi(t)，滤波后的陆地水储量为Δηi(t)，通过最小二乘原则(Landerer and Swenson, 2012)：

(5)

其中，N表示时间序列个数，SF表示所求尺度因子.

(5)计算水文模型陆地水储量与CSR-SHc数据的均方根误差和相关性

皮尔逊相关系数(PR，Pearson correlation coefficient)通常用来描述两个时间序列之间在一致性，均方根误差(RMSE，Root Mean Square Error)反映了模拟值与“真实”值之间的偏差程度.对于任意格网点相关系数计算如下(Yin et al., 2020)：

(6)

对于任意格网点均方根误差计算如下(Yin et al., 2020)：

(7)

则该格网点分配的权重为

W=PR+RMSE-1,

(8)

其中，W表示权重，N表示时间序列个数，Δh(t)表示t时间序列水文模型的陆地水储量，ΔSH(t)表示t时间序列CSR-SHc的陆地水储量.

(6)基于模型权重加权尺度因子

基于尺度因子模型的权重，可得新型尺度因子校正法的公式如下：

(9)

其中，NSFCM-SF表示新型尺度因子校正法计算的尺度因子，SF表示模型尺度因子，n表示尺度因子模型数目(CLM，MOS，NOAH，VIC，WGHM).

2 研究区域和数据

2.1 研究区域

长江流域地处中国南部，经纬度范围为85°E—122°E和24°N—36°N，横跨中国西部、中部和东部三大区域共计19个省、市和自治区，长江流域位置及其海拔高度如图1所示.长江流域年平均降水量为1067 mm且时空分布不均匀，约70～90%集中在5—10月.长江干流宜昌以上为长江上游，流域面积为100×104km2，宜昌至湖口为长江中游，流域面积为68×104km2，湖口至长江入海口为长江下游，流域面积为12×104km2.鉴于长江下游流域面积较小，将GRACE数据应用于小于GRACE分辨率的区域，反演结果可能存在较大不确定性(Long et al., 2015b).因此，本研究将长江中游和下游统称为长江中下游，如图1虚线以东所示，虚线以西为长江上游.

图1 长江流域位置及其海拔高度Fig.1 Location and elevation of the Yangtze River basin

2.2 GRACE数据

本研究使用的数据包括两种网格数据产品：CSR-SHc和CSR-M(CSR Mascon).CSR-SHc数据是基于美国德克萨斯大学空间研究中心(Center for Space Research，CSR)、美国喷气推进实验室(Jet Propulsion Laboratory，JPL)和德国波兹坦地球科学研究中心(GeoForschungsZentrum，GFZ)的RL05球谐系数得到，数据最大阶次为60阶，时间跨度范围为2002年4月—2017年1月，利用三次样条插值得到月份缺失的数据，每月网格数据表示相对于2004年1月—2009年12月平均值的质量变化.GRACE的C20值比卫星激光测距(Satellite Laser Ranging，SLR)的C20值具有更大的不确定性，因此，用SLR的C20系数解代替GRACE的C20系数(Cheng et al., 2011).使用Swenson等(2008)的方法估计一阶项，应用于地心校正.基于Geruo等(2012)的模型进行冰川均衡调整(Glacial Isostatic Adjustment，GIA)校正.为了减少GRACE月度图中南北条带的相关误差，对数据进行了去条带处理.此外，数据还应用了300 km高斯平滑滤波以抑制球谐系数中的高阶随机误差.

CSR-M数据是基于Level-1B数据解算的，没有应用任何平滑或去条带滤波.因此，CSR-M数据无需进一步处理即可使用.CSR-M数据空间分辨率为0.5°× 0.5°，该数据集并非针对特定应用定制，适用于所有感兴趣的科学领域，例如海洋学、陆地表面水文学、冰冻圈等研究.研究表明，CSR-M数据解适用于恢复全球和大流域尺度(大于约3°× 3°)中不同质量源的质量再分配，而球谐系数解在有效处理时，在较小或局部尺度上优于CSR-M数据解(Zhang et al., 2019).

2.3 水文数据

全球陆地同化系统水文模型(Global Land Data Assimilation System，GLDAS)，是由美国国家航空局和美国国家海洋与大气管理局共同开发(Rodell et al., 2004).该模型通过陆地表面建模和数据同化技术，输出陆地表面各项参数(土壤温度和土壤水分、蒸发量和地表径流等)和水文变量.本文使用的GLDAS数据包括四个陆地表面模型NOAH、CLM、MOS和VIC：NOAH是一种土地表面过程模型，既将土壤表面和植被视为一个整体的模型；CLM包含三个重要的子模型：美国国家大气研究中心的地表模型，生物圈—大气转移计划和中国科学院大气物理研究所的陆地表面模型；MOS是根据单元内植被类型的分布，将每个模型网格单元划分的模型；VIC模型是基于热平衡和物理动力机制的大规模半分布式水文模型.本研究使用的数据空间分辨率为1°×1°，时间分辨率为1个月，选取的数据时间范围与GRACE数据一致.GLDAS模型具有不同深度的土壤水分，VIC模拟3层，最大深度1.90 m；NOAH模拟4层，最大深度2.00 m；CLM模拟10层，最大深度3.43 m；MOS模拟3层，最大深度3.50 m.本研究模型在计算TWS时均考虑了土壤水分、雪水当量和总冠层蓄水量，并且每月数据均扣除了2004年1月—2009年12月数据平均值获得TWS相对变化量.

WGHM(Water GAP Global Hydrology Model)水文模型由德国法拉克福大学自然地理研究所研制(Döll et al., 2014)，不同于GLDAS水文模型，该模型包含地下水成分，时间跨度范围为2000年1月—2016年12月，空间分辨率为0.5°×0.5°，每月一值.本研究选取了2002年4月—2016年12月的数据，为与CSR-SHc数据空间分辨率保持一致，对WGHM水文模型数据进行了插值.

2.4 降水数据

中国地面降水月值数据集是基于国家气象中心基础资料专项最新整编中国地面台站降水资料(赵煜飞等, 2014).该数据利用薄盘样条法进行空间插值，从而生成1961年至今的分辨率0.5°×0.5°的中国降水月值格点数据，数据范围为72°E—136°E和18°N—54°N.本文选取的数据时间范围为2002年4月—2017年1月，与CSR-SHc数据时间范围一致.

3 新型尺度因子校正法的验证

3.1 数据处理

3.1.1 截断及滤波引起的信号衰减

GRACE球谐系数截断以及去条带和高斯平滑滤波处理会导致GRACE反演结果存在泄露误差(冯伟等, 2017)，主要分为两种情况：(1)当研究区内部信号较强，而周围信号较弱时，球谐系数截断以及滤波后，研究区内部的信号会泄露到研究区外，则表明该区域由外泄露误差控制，称为外泄露误差；(2)当研究区内部信号较弱，而周围信号较强时，球谐系数截断以及滤波后，周围的信号会泄露到研究区内部，则表明该区域由内泄露误差控制，称为内泄露误差.

本研究利用GLDAS-NOAH水文模型计算出TWS空间分布趋势变化，扣除研究区域长江流域外的数据，模拟长江流域球谐系数截断及滤波引起的外泄露误差.图2a是由2002年4月—2017年1月的GLDAS-NOAH数据计算的长江流域TWS空间分布趋势变化.将图2a中的数据经过球谐展开并截断至60阶，结果如图2b所示.图2c和图2d分别采用300 km和500 km高斯平滑滤波处理所得.可以看出图2a中A1区域的取值范围为0.50～1.00 cm·a-1，图2b中A2区域的取值范围为0.25～0.50 cm·a-1，图2c中A3区域的取值范围为0.25～0.30 cm·a-1，图2d中A4区域的取值范围为0.10～0.15 cm·a-1，从A1→A2→A3→A4标注的红色区域不断向周围扩散信号，而且该区域的TWS空间分布趋势变化不断变小.可以得出经球谐系数展开及截断后信号向周围扩散，再经过高斯平滑滤波后信号进一步衰减，并且高斯平滑半径越大，信号损失越强.

图2 GLDAS-NOAH模型计算的长江流域TWS空间分布趋势变化Fig.2 The spatial distribution trends of TWS in the Yangtze River basin calculated by the GLDAS-NOAH model

为了评估截断及滤波的外泄露误差效应，分别计算了长江流域内的区域平均趋势变化进行对比分析，定义趋势变化损失率(Rate)为

Rate=(NOAH-Model)/NOAH,

(10)

其中，NOAH为GLDAS-NOAH水文模型的区域平均趋势，Model为水文模型滤波后的区域平均趋势.

据表1可知，初始数据计算长江流域的区域平均趋势为0.08 cm·a-1；经过球谐展开并截断至60阶，不进行高斯平滑的区域平均趋势为0.07 cm·a-1，相对于初始数据信号损失12%；经过球谐展开并截断至60阶，进行300 km的高斯平滑的区域平均趋势为0.06 cm·a-1，相对于初始数据信号损失25%；经过球谐展开并截断至60阶，进行500 km的高斯平滑的区域平均趋势为0.05 cm·a-1，相对于初始数据信号损失37%.这些数据也验证了截断以及滤波造成的信号衰减.因此，恢复GRACE数据由于截断以及滤波造成的信号损失，对于利用GRACE球谐数据反演研究区域质量变化具有重要意义.

表1 不同实验模型计算的区域平均趋势变化Table 1 Regional average trend changes calculated by different experimental models

3.1.2 尺度因子模型

本文使用了GLDAS四种水文模型和WGHM计算的尺度因子模型，基于新型尺度因子校正法得出的模型如图3a所示.全球水文模型CLM、MOS、NOAH、VIC和WGHM数据基于公式(5)应用于1°×1°的所有网格点，扣除边界外数据得到长江流域网格化的尺度因子，如图3b、3d、3e、3f和3g所示.为比较尺度因子结果，本文使用了由官方提供的CLM4.0尺度因子(Landerer and Swenson, 2012)，将CLM4.0尺度因子保留长江流域区域内数据如图3c所示.尺度因子CLM-SF(图3b)、MOS-SF(图3d)、NOAH-SF(图3e))和VIC-SF(图3f)空间分布类似，NSFCM-SF(图3a)、CLM4.0-SF(图3c)和WGHM-SF(图3c)相比于图3b、3d、3e、3f在长江下游的空间分布明显不同.这些差异也反映了不同模型TWS数据空间分布的差异，这是由于不同水文模型数据考虑的组成成分不同造成.

图3a—3f模型的尺度因子网格值的空间分布基本都大于1，而图3g模型的尺度因子网格值的空间分布基本都小于1，这是由于WGHM水文模型与GLDAS水文模型相比包含了地下水成分，不同模型的尺度因子显示出明显的差异和空间异质性.本文利用不同模型的尺度因子计算了长江流域内空间平均尺度因子如表2所示，长江流域的GLDAS模型空间平均尺度因子都大于1，其中CLM4.0-SF空间平均尺度因子最大(1.42)，WGHM-SF空间平均尺度因子最小(0.88)，空间平均尺度因子从一定程度反映了对校正结果周年振幅和长期趋势的影响.

表2 不同模型计算的长江流域空间平均尺度因子Table 2 Spatial mean scale factor of The Yangtze River Basin calculated by different models

图3 不同模型计算的长江流域尺度因子Fig.3 Scale factors of the Yangtze River Basin calculated by different models

3.1.3 迭代恢复法

迭代恢复法由Chen等(2006)首次提出，经过改进后，广泛应用于GRACE数据泄露误差校正.迭代恢复法的具体步骤如下：(1)在1°×1°网格的每个网格点上(图4b)，令初始质量变化率为GRACE的质量变化率(图4a).(2)正演模型质量变化率(图4c)由步骤(1)的1°×1°网格化模型质量变化率(图4b)，进行球谐展开，并将球谐系数截断至60阶，0阶和1阶的系数设为0，然后应用300 km的高斯平滑滤波得到.(3)在每个格网点上，将初始的GRACE的质量变化率(图4a)和正演模型质量变化率(图4c)之间的差值乘以1.20，回加到1°×1°网格化模型质量变化率(图4b)中，作为新的模型质量变化率进行步骤(2)，重复该过程.连续迭代后产生正演模型质量变化率和初始GRACE质量变化率之间的一致性增加(图4a与图4c).(4)当初始的GRACE的质量变化率(图4a)和正演模型质量变化率(图4c)之间的差值小于特定的阈值或者迭代特定次数后，停止迭代.

经过迭代60次后，将得到球谐信号转化为格网信号即为迭代恢复法恢复的信号结果(图4b)，比较图4a和图4c发现具有较强的一致性，他们之间的差值(图4d)较小(色棒值-0.10～0.10 cm·a-1).因此，认为迭代恢复法恢复的结果信号是可信的.

3.2 新型尺度因子校正方法对比分析

3.2.1 长江流域水储量变化空间分布趋势对比分析

根据GLDAS水文模型、WGHM水文模型和新型尺度因子校正法计算的尺度因子和官方提供的CLM4.0尺度因子获得校正后的GRACE数据时变信号如图5所示.可以看出，尺度因子法(WGHM水文模型除外)和迭代恢复法(图4b)校正结果相比于CSR-SHc反演结果(图5e)空间分布趋势信号明显得到增强，尤其是在长江中下游.

为进一步对比不同泄露误差校正方法恢复结果，本文以CSR-M数据校正结果作为CSR-SHc数据的真实信号.不同泄露误差校正方法恢复结果空间分布趋势与CSR-M数据空间分布趋势相比，NSFCM-Trend(图5a)、CLM-Trend(图5b)、MOS-Trend(图5g)、NOAH-Trend(图5h)、VIC-Trend(图5i)与CSR-M-Trend(图5d)一致，而CLM4.0-Trend(图5c)标注区域C与CSR-M-Trend相比(图5d)明显存在过量恢复现象.WGHM-Trend(图5f)标注区域D和FM-Trend(图4b)标注区域B与CSR-M-Trend相比(图5d)明显存在信号恢复不全现象.对于CSR-SHc-Trend(图5e)标注E区域，GLDAS模型增强了该区域的陆地水储量亏损信号，CLM4.0-Trend(图5c)和WGHM-Trend(图5f)在该区域的信号与CSR-M-Trend接近，NSFCM-Trend(图5a)与GLDAS模型相比在该区域的信号更接近CSR-M-Trend.本文也计算了长江流域不同泄露误差校正方法恢复结果的区域平均趋势如表3所示.由区域平均趋势可以得出CSR-M-Trend最大，WGHM-Trend最小，NSFCM-Trend区域平均趋势大于GLDAS模型计算的尺度因子校正结果的区域平均趋势，但是小于CLM4.0-Trend和FM-Trend区域平均趋势.综上所述，NSFCM-Trend避免了FM-Trend、CLM4.0-Trend和WGHM-Trend存在信号恢复不全和过量恢复现象，并且NSFCM-Trend区域平均趋势比GLDAS水文模型尺度因子校正结果的区域平均趋势较好.

图4 迭代恢复法恢复GRACE质量变化信号Fig.4 Forward modeling method to restore GRACE quality change signal

图5 长江流域不同恢复方法的空间分布趋势Fig.5 Spatial distribution trends of different restoration methods in the Yangtze River basin

表3 长江流域不同恢复方法的区域平均趋势Table 3 Regional average trend of different restoration methods in the Yangtze River basin

3.2.2 长江流域水储量变化时间序列对比分析

图6显示了长江流域不同泄露误差校正方法恢复结果得到的水储量变化时间序列.不同泄露误差校正方法恢复结果与CSR-M数据相比周年相位一致性较好，但周年振幅具有明显差异，表明尺度因子法和迭代恢复法主要改变水储量变化时间序列的周年振幅，而对水储量时间序列的周年相位影响较小.为了进一步比较不同泄露误差校正方法恢复结果的差异，本文计算了长江流域水储量变化时间序列的周年振幅和长期趋势如图7所示.

图6 长江流域不同泄露误差校正方法恢复的TWS时间序列Fig.6 TWS time series recovered by different leakage error correction methods in the Yangtze River Basin

图7 长江流域不同泄露误差恢复方法恢复的TWS周年振幅和长期趋势 (a) 周年振幅； (b) 长期趋势.Fig.7 Annual amplitude and long-term trend of TWS recovered by different leakage error recovery methods in the Yangtze River Basin (a) Annual amplitude; (b) Long-term trends.

图7显示了长江流域不同泄露误差校正方法恢复的水储量变化长期趋势和周年振幅，可以得出CSR-SHc数据的周年振幅和长期趋势处于较低水平，存在信号泄露现象.泄露误差校正后，周年振幅和长期趋势均得到增强，信号泄露得到一定程度上的恢复.在周年振幅上，TWS-WGHM、CSR-M数据校正结果小于CSR-SHc数据，GLDAS水文模型和CLM4.0尺度因子校正结果均明显大于CSR-SHc数据.因为GRACE原始观测数据模拟的TWS在精细的空间比例上不敏感(Landerer and Swenson, 2012)，而CSR-M和CSR-SHc主要受GRACE卫星原始观测数据约束，因此CSR-M、CSR-SHc和基于FM法恢复的周年振幅较小；GLDAS水文模型模拟的TWS大多数变化发生在次季节与季节的时间尺度上，因此由模型模拟的TWS计算的尺度因子优化了周年振幅.在长期趋势上，WGHM模型计算的尺度因子降低了CSR-SHc数据的长期趋势；FM法、CLM4.0尺度因子、NSFCM尺度因子和GLDAS水文模型计算的尺度因子优化了CSR-SHc数据的长期趋势.其中，TWS-FM(0.34±0.08 cm·a-1)、TWS-CLM4.0数据(0.34±0.11 cm·a-1)和TWS-NSFCM(0.29±0.09 cm·a-1)较好，更接近于CSR-M数据(0.39±0.07 cm·a-1)的长期趋势.

4 新型尺度因子校正法的应用

4.1 长江子流域水储量变化对比分析

NSFCM-SF尺度因子综合GLDAS水文模型计算的尺度因子校正结果空间分布趋势的优点，避免了CLM4.0-SF、FM和WGHM校正结果空间分布趋势不均匀的现象.在长期趋势上，NSFCM-SF尺度因子校正结果虽然低于CLM4.0-SF和FM，但优于GLDAS水文模型校正结果计算的长期趋势，在周年振幅上，NSFCM-SF尺度因子校正结果明显优于FM和CSR-M的结果.综合考虑，本文以下主要讨论NSFCM-SF尺度因子校正结果.选取CSR-M、CSR-SHc、CLM4.0和TWS-NSFCM模型，进一步分析NSFCM-SF在长江子流域(长江上游和长江中下游)的校正结果.

图8显示了长江子流域CSR-M、CSR-SHc、TWS-CLM4.0和TWS-NSFCM模型的水储量时间序列.据图8可知，长江子流域的TWS-NSFCM和TWS-CLM4.0周年振幅明显优于CSR-SHc的结果，长江上游TWS-NSFCM周年振幅优于TWS-CLM4.0，长江中下游TWS-CLM4.0周年振幅优于TWS-NSFCM.长江中下游的时间序列相比于长江上游更复杂，季节性的TWS变化更大，这可能与长江中下游季节性地表水存储变化有关(Long et al., 2015a).本文计算的长江子流域水储量变化长期趋势和周年振幅如表4所示.长江子流域的长期趋势，TWS-NSFCM与TWS-CLM4.0相比更接近于CSR-M数据.长江上游CSR-SHc长期趋势为0.07 cm·a-1，经NSFCM-SF尺度因子校正后TWS-NSFCM长期趋势为0.14 cm·a-1，长期趋势提高了100%，长江中下游CSR-SHc长期趋势为0.40 cm·a-1，经NSFCM-SF尺度因子校正后TWS-NSFCM长期趋势为0.49 cm·a-1，长期趋势提高了23%，可以得出长江流域水储量上升趋势主要集中在长江中下游.

表4 长江子流域的TWS长期趋势和周年振幅Table 4 Long-term trends and annual amplitudes of TWS in the Yangtze River sub-basin

图8 长江子流域的TWS时间序列 (a) 长江上游； (b) 长江中下游.Fig.8 TWS time series of the Yangtze River sub-basin (a) The upper Yangtze River; (b) The middle and lower Yangtze River.

4.2 长江流域降水与水储量变化的比较

本文基于NSFCM-SF尺度因子校正结果时间序列与国家气象中心的降水月值数据进行了比较，如图9所示.可以看出长江流域以及子流域降水主要集中在5—10月，降水明显增多时，TWS-NSFCM水储量变化明显增加，并且滞后降水数据1—2月.降水数据变化与TWS-NSFCM水储量变化一致，这表明降水是长江流域水储量变化的主要原因.在2002年4月—2007年1月长江流域和子流域的水储量和降水均呈现上升趋势，长江流域水储量上升趋势为0.29 cm·a-1，降水上升趋势为0.06 cm·a-1；上游水储量上升趋势为0.14 cm·a-1，降水上升趋势为0.03 cm·a-1；中下游水储量上升趋势为0.49 cm·a-1，降水上升趋势为0.10 cm·a-1.从而可以得出降水的上升趋势主要集中在中下游，这与长江流域水储量上升趋势主要集中在中下游的特征一致，进一步验证了降水是长江流域水储量变化的主要原因.本文也计算了长江流域年降水异常如图10所示，由于2010年、2011年和2015年降水数据不完整，文中不讨论降水数据不完整的年份.据图10可知，多数年份降水量接近或低于多年平均值，但2012年、2014年和2016年降水量明显高于多年平均值，这可能与近年来全球变暖、气候极端事件频发有关(例如厄尔尼诺等).结合图9，降水较多的年份，TWS-NSFCM水储量年际变化明显增多，因此，长江流域TWS-NSFCM水储量年际变化特征与降水的年际变化规律一致.

图10 长江流域年降水异常Fig.10 Annual precipitation anomalies in the Yangtze River Basin

5 结论

GRACE卫星的出现为研究大尺度区域的TWS变化提供了有效的手段.然而，GRACE数据处理中的泄露误差是反演质量变化的主要限制因素，传统尺度因子法校正泄露误差依赖水文模型结果，并且不同的水文模型由于模型的结果以及驱动数据的差异，导致模型输出结果间也存在一定的差异.因此，本文提出新型尺度因子校正法，旨在综合多种尺度因子模型提高GRACE陆地水储量变化的准确性.

(1)针对GRACE数据在信号处理过程中存在的泄露误差，本文在传统尺度因子法的基础上，根据水文模型的陆地水储量与GRACE-SHc数据的均方根误差和相关性赋予权重，构建了新型尺度因子校正法.相比于传统尺度因子法使用单一水文模型计算的尺度因子，新型尺度因子校正法充分考虑多种模型与GRACE-SHc数据一致性和准确性信息，校正数据处理过程中的泄露误差，恢复真实的水文信号，从而提高GRACE反演区域陆地水储量变化的准确性.这对于更好地理解区域陆地水储量变化，进一步认识洪涝灾害以及对区域水资源管理具有重要意义.

(2)对新型尺度因子校正法在长江流域进行验证，空间分布趋势上，GLDAS-SF和NSFCM-SF校正结果与CSR-M相近，而WGHM-SF、FM和CLM4.0-SF校正结果与CSR-M相比存在信号恢复不全和过量恢复现象.GLDAS-SF、CLM4.0-SF和NSFCM-SF校正结果较好地恢复了区域水储量的周年振幅，但长期趋势恢复较弱，FM方法主要恢复了区域水储量的长期趋势，相比于尺度因子周年振幅恢复结果较弱，而NSFCM-SF与GLDAS-SF相比，长期趋势优于GLDAS-SF校正结果.

(3)基于新型尺度因子校正法校正结果显示，长江流域、上游和中下游在2002年4月—2017年1月水储量呈现上升趋势，分别为0.29 cm·a-1、0.14 cm·a-1和0.49 cm·a-1，相比于长江流域、上游和中下游CSR-SHc上升趋势0.21 cm·a-1、0.07 cm·a-1和0.40 cm·a-1，分别提高了38%、100%和23%.长江流域水储量上升趋势主要集中在中下游.长江流域降水上升趋势为0.06 cm·a-1，上游降水上升趋势为0.03 cm·a-1，中下游降水上升趋势为0.10 cm·a-1，降水的上升趋势主要集中在中下游，这与长江流域水储量上升趋势主要集中在中下游的特征一致，进而表明降水是长江流域水储量变化的主要原因.本研究方法对青藏高原冰川变化亏损信号泄露入长江流域校正结果还存在不足之处，后期会根据一些实测数据进一步改进本研究方法.

致谢感谢喷气推进实验室提供的GRACE每月质量网格陆地产品数据(https:∥grace.jpl.nasa.gov/data/get-data/monthly-mass-grids-land/)；感谢德克萨斯大学空间研究中心提供的GRACE网格数据产品CSR Mascon数据；感谢戈达德地球科学数据和信息服务中心提供的GLDAS水文模型数据(https:∥disc.gsfc.nasa.gov/)；感谢德国法拉克福大学自然地理研究所提供的WGHM水文模型数据；感谢国家气象科学数据中心提供的降水数据(http:∥data.cma.cn/)；感谢GMT开源画图软件和冯伟博士的开源代码.感谢两位匿名审稿人的意见，这些意见有助于提高文章质量.杨帅、郑伟和尹文杰为共同第一作者，郑伟和尹文杰为共同通讯作者.