不同去相关方法对GRACE反演长江流域陆地水储量变化的影响

季新宇，李廷豪

(南京海河测绘科技有限公司，江苏 南京 211101)

0 引 言

长江总长约6 300 km，是仅次于亚马孙河和尼罗河的世界第三大长河，长江流域总面积约为190万平方公里。粮食产量约占中国的一半，人口约占中国的三分之一，在中国社会发展进程中发挥着重要的作用。随着长江流域社会经济的快速发展，水资源问题越来越成为一个十分突出的因素[1]。

GRACE(Gravity Recovery and Climate Experiment)重力卫星计划是由美国宇航局和德国航空航天中心联合发射，提供了一种基于空间测量来监测全球时变重力场的手段。GRACE重力卫星对陆地水储量的变化同样敏感，可探测平均小于1 cm的陆地水储量变化(TerrestrialWater Storage Change,TWSC)。在流域尺度研究方面，胡小工等利用2002年4月至2003年12月之间的数据分析了长江流域水储量的变化，结果表明长江流域水储量周年变化幅度可达3.4 cm的等效水高，最大值出现在春季和秋季[2]。Chen等通过研究成功地监测了亚马孙流域2005年的干旱和2009年的洪水[3-4]。冯伟等研究了亚马孙流域在2010年经历的一次干旱的过程[5]。尼胜楠等利用GRACE时变引力场数据并结合水文模型数据研究了长江流域和黄河流域的陆地水储量变化，发现长江、黄河流域的水储量变化与降雨量在不同尺度和频段存在较好的时空变化一致性[6]。

由于GRACE卫星轨道的设计和传感器误差以及解算时背景模型误差等一系列因素的影响，利用GRACE数据解算的原始时变引力场模型中包含明显的南北方向条带误差，在一定程度上掩盖了真实的陆地水储量变化信号。为了研究不同去相关方法对反演长江流域陆地水储量变化的影响，本文分别利用SWENSON、P3M6、P4M15三种方法对2003年1月至2014年5月之间长江流域的陆地水储量变化进行反演，与不进行去条带处理结果及GLDAS水文模型计算的结果进行比较，从而分析不同去条带误差方法对GRACE反演长江流域陆地水储量变化的空间和时间特征的影响。

1 反演陆地水储量基本原理

利用引力场模型的球谐系数变化计算地球表面质量变化(以等效水高的形式表示)的模型为：

(1)

(2)

式中，a指地球平均半径，r表示高斯滤波半径，本文r取值500 km。

由于卫星轨道设计和观测精度的制约，直接利用位模型反演的陆地水储量变化网格值存在明显的南北向条带误差[8]，对于某固定次，偶(奇)数阶之间存在明显的相关性，并提出利用滑动多项式消除相关性。以Cl,m为例，对Cl-2α,m,…Cl-2,m,Cl,m,Cl+2,m,…Cl+2α,m进行4阶多项式拟合，然后从原始值中扣除拟合值，从而得到去相关之后的值。关于平滑窗口ω(用于多项式拟合的球谐系数的个数)和α之间的关系为ω=2α+1。(Swenson & Wahr,2006)对5次以上的球谐系数进行了去条带处理，其平滑窗口宽度与球谐系数的次数m之间的关系为：

(3)

式中，max()表示两者之的最大值，文献[8]里选择的经验参数值为A=30和K=10。上述方法(下文简称为SWENSON法)是利用阶数来确定滑动窗口的宽度，并进行滑动多项式拟合，而Chen等[3]在计算Sumatra-Andaman地震引起的重力变化时，采用了P3M6方法，即未选择滑动窗口，而是直接对同一次的所有偶(奇)数阶进行多项式拟合，即，保持球谐系数的前m阶保持不变，然后对m阶及m阶以上的系数用n阶多项式进行拟合，奇数项和偶数项分开拟合，将拟合后的值从原始值扣除之后便得到曲线管制后的值，这种方法简称为PnMm方法。同样的，Chambers等[9]在利用GRACE RL04数据和RL05数据比较研究海底压力时使用了P4M15方法。

2 实验验证与分析

2.1 数据来源

实验数据采用美国德克萨斯大学空间研究中心(CSR)发布的GRACE GSM RL05引力场模型数据，研究时段为2003年1月至2014年5月(缺失9个月的数据，后续计算出的陆地水储量变化时间序列对缺失月份数据进行线性插值)。该产品已经扣除了大气、海洋、海潮、固体潮等等信号，主要反映了由于地表质量重新分布引起的时变重力信号。SLR获得的C20项和一阶项代替GRACE数据的C20项和一阶项[10-11]。

为了验证GRACE反演结果的有效性，采用全球陆地数据同化系统(Global Land Data Assimilation System，GLDAS)中的Noah陆地表面模型进行验证。该系统是美国宇航局和美国国家环境预报中心(National Centers for Environmental Prediction，NCEP)国家海洋大气局联合项目，时间分辨率为一个月，1°×1°的空间分辨率，选取四层土壤湿度数据、雪深水当量数据模拟陆地水储量。

2.2 时间序列分析

根据式(1)，分别采用SWENSON、P3M6、P4M15、NONE(不进行去相关)和500 km的高斯滤波组合的方法计算全球陆地水储量变化的等效水高格网值，再依据长江流域轮廓，对流域内的格网值数据按照纬度余弦加权的方法获得长江流域平均陆地水储量变化的时间序列：

(4)

为了更合理地比较GRACE反演结果与水文模型计算结果，首先对GLDAS-Noah的格网数据做球谐展开，并截断至60阶，然后进行半径为500 km的高斯滤波处理，最后同样采用纬度余弦加权的方法计算长江流域的TWSC时间序列。结果如图1所示。

由图1可知，GRACE反演的长江流域陆地水储量变化与GLDAS水文模型计算的结果均表现出明显的季节性变化，水储量变化极小值出现在春季，极大值出现在秋季，在季节尺度上，GRACE反演结果与GLDAS模拟结果具有较好的一致性，SWENSON、P3M6、P4M15、NONE 4种情况下反演的TWSA与GLDAS模拟结果之间的相关系数分别达到0.886、0.878、0.883、0.869。另外，从上图明显看出，利用GRACE位模型反演的陆地水储量变化值均高于GLDAS模拟结果，这是由于GLDAS-Noah模拟的陆地水储量变化只包含了雪深水当量变化值和土壤深度为0～200 cm之间水储量变化，而GRACE监测的陆地水储量变化不仅包含了土壤水和冰雪水变化，还包含了地下水、地表水等一系列变化[12]。

为了进一步比较GRACE不同去条带误差方法反演的陆地水储量变化与GLDAS模拟结果在季节上的变化规律[13]，对时间序列进行了拟合分析，结果如表1所示。

表1 GRACE不同去相关方法反演结果与GLDAS模拟长江流域TWSC的周年项和半周年项的振幅、相位

由表1可知，基于GRACE数据，利用SWENSON、P3M6、P4M15去相关方法计算的长江流域陆地水储量变化的周年振幅关系是：P4M15> P3M6>NONE>SWENSON，半周年振幅：P3M6>P4M15>NONE>SWENSON，GRACE反演结果的周年、半周年振幅均大于GLDAS模拟结果。图2为采用以上几种方法得到的长江流域陆地水储量年平均变化空间分布趋势图。

由图2可知，在长江流域中下游，利用SWENSON方法反演的陆地水储量变化趋势的空间分布与其余3种情况相差较大，这可能是由于SWENSON方法在使用时对于A与K两个经验值在选择上，采用的是文献[8]里的原始参数，并未进行改动，容易产生较大误差。另外，就P3M6与P4M15比较而言，在长江流域上游，P4M15反演出的陆地水储量变化趋势幅值明显大于P3M6的反演结果，而在中下游地区，P3M6反演出的陆地水储量变化趋势幅值明显大于P4M15的反演结果。

3 结 语

(1)GRACE反演的陆地水储量变化时间序列的与GLDAS-Noah模拟的结果较为一致，相关系数均达到0.8以上，另外，GRACE反演结果的幅值均大于GLDAS模拟的幅值，这是由于GRACE反演的陆地水储量变化不仅包括了土壤水和雪水当量变化，还包含了地下水、地表水等一些列水储量变化。

(2)在时间序列上，不同去相关方法反演的陆地水储量变化周年项的半周年项的振幅接近，在空间分布上，不同去相关方法反演的陆地水储量变趋势的差异较大，但P3M6和P4M15反演结果的分布十分接近，而与SWENSON相差较大。在长江流域上游，P4M15方法反演的陆地水储量变化趋势信号强度最强，而在中下游地区，SWENSON方法反演的陆地水储量变化趋势信号强度最强。