基于GRACE RL06 数据探测华北平原水储量变化

何国庆

（1. 上海市测绘院，上海 200063）

GRACE 卫星由美国国家航空航天局（NASA）和德国航空中心（DLR）共同研制、部署，其科学意义是探测地球时变重力场信号分析地球质量密度变化与地球重力场的时变特征[1-9]。本文以华北平原为研究对象，利用最新的GRACE RL06重力卫星数据反演华北平原近年来水储量变化，分析其时空变化特征，旨在实现华北平原水储量时空连续监测，为华北平原合理开发利用水资源提供依据。

1 数据与方法

1.1 GRACE重力数据

GRACE之所以能反演水储量变化，是基于John[10]的理论依据，利用GRACE 时变重力场模型数据反演水储量变化时，通常以等效水高的形式表示：

式中，θ为地心余纬；λ为地心经度；a为地球半径；ρave为地球平均密度（ρave=5 517 kg/m3）；kl为l阶负荷勒夫数；ρw为水的密度（ρw=1 000 kg/m3）；l、m为球谐系数的阶数和次数；Pˉlm(cosθ)为完全规格化的缔合勒让德函数；ΔClm、ΔSlm为完全规格化的球谐系数变化量；Wlm为平滑函数。

本文采用GSR公布的自2002-08～2016-08共154个月的RL06 版本GSM 数据，其中存在部分月份数据丢失的情况。RL06 数据相比于RL05 采用了新的背景模型，改善了数据处理方法，使得数据精度明显提高。在利用式（1）进行计算时需要注意：受卫星轨道及星载仪器等影响，GRACE球谐系数中低阶项（一阶项和二阶项）观测精度较差，需要对其进行替换；GRACE时变重力场数据中高阶项存在着较大的误差，并且误差随着阶数的升高而增大，对于该项误差常用高斯滤波或各向异性滤波来处理；GRACE球谐系数奇（偶）数阶项数之间存在一定的相关性，这种相关性导致了反演结果在图上表现为“南北条带”误差，针对该误差处理方法，通常采用去相关滤波进行处理。

1.2 Mascon数据

Mascon 数据由CSR、GFZ、JPL 三家机构利用GRACE Level-1 数据各自解算的地球重力场模型数据。Watkins[11]等国外学者通过实验发现Mascon 数据较球谐系数反演的重力场模型结果有更高的时空分辨率和信噪比。Mascon数据是重力场模型的另一种表达方式，数据已经替换过了低阶项、改正了GIA（冰川均衡调整）等步骤，不用进行后处理和空间滤波平滑。与球谐系数不同的是，球谐系数没有特定的位置信息，而Mascon数据中每个质量集中块都具有特定的地球物理位置。

1.3 水文数据

收集了研究期间内的水文模型数据和降水数据。其中水文模型采用GLDAS 数据，GLDAS 由美国国家航空航天局NASA 和美国国家环境预报中心NCEP 共同建立，作为高精度的水文模型，能通过卫星与地面观测数据约束模型，提供最优化的近实时的地球表面质量变化信息[12]。降水数据来自GPCC，它是由德国气象局DWD 成立的，根据地面观测数据整合了全球陆地表面降水数据集，并且还接受了世界气象组织和全球电信传输系统的每天地面天气和月气候数据[13]。GLDAS 提供的水文模型数据有Noah、VIC、Mosaic和CLM；其时间分辨率有3 h和1个月，空间分辨率提供0.25°×0.25°和1°×1°；GPCC数据提供了4种空间分辨率：0.25°×0.25°、0.5°×0.5°、1°×1°、2.5°×2.5°。本文均采用1°×1°、时间分辨率为1个月的水文数据，时间跨度选取与GRACE数据相同。

1.4 基于尺度因子的恢复

由于GRACE 重力模型数据提供的球谐系数是有限的，无法完全展开，这就造成了部分信号无法反映的截断误差；另一方面，使用滤波方法对球谐系数进行处理时，也会产生误差，使得GRACE 估计的结果是有偏差的。本文将采用尺度因子方法来恢复上述内容中“泄露”的信号，即利用水文模型数据来模拟泄露误差的影响，再利用尺度因子恢复GRACE 反演信号。具体操作步骤是：①对水文模型数据进行球谐展开，并展开到和GRACE 数据一样的阶数；②对展开后的水文模型球谐系数进行和GRACE数据同样的后处理操作，得到研究区域的质量变化时间序列ΔS′；③计算原始水文数据得到的研究区域质量变化的时间序列ΔS；④将处理后的时间序列ΔS′和原始数据的时间序列ΔS进行最小二乘拟合，计算出尺度因子k；⑤对GRACE反演的研究区域质量变化时间序列信号乘上尺度因子k。其中尺度因子k满足处理后的时间序列ΔS′要和原始数据得到的时间序列ΔS的残差平方和最小。即：

2 结果分析

2.1 数据处理

本文采用CSR提供的GRACE RL06 GSM重力模型数据产品，截取时间跨度为2002-08～2016-08 共计154 个月份。对球谐系数中一阶项利用Chambers 解算出的一阶项数据进行替代，对球谐系数中二阶项数据利用SLR 提供的数据进行替代。对于南北条带误差，采用Duan 滑动窗口去相关方法进行去条带处理。采用平滑半径为300 km 的扇形滤波方法抑制高阶项噪声。最后，对得到的全球水储量变化格网数据进行按纬度的余弦值定权，并根据研究区域的选择，对研究区域的格网数据进行加权平均计算出该研究区域的平均水储量变化的时间序列。

利用式（2）计算出在Duan 滑动窗口去相关滤波和300 km扇形滤波方法处理下的尺度因子k，这里计算出的尺度因子k=1.573，本文利用此值乘以GRACE计算出的结果，作为GRACE反演水储量的最终估计值。

为了分析华北平原水储量变化在时间尺度上的特征，通过最小二乘拟合对华北平原水储量时间序列进行以下形式拟合：

式中，b为趋势项；Ai为振幅；φi为相位；Ti为周期（T1为周年，T2为半周年）。通过以上拟合公式可以求出华北平原水储量变化的年际变化趋势、周年和半周年项等变化特征。

2.2 GRACE反演结果

华北平原是我国第二大平原，坐落于我国东部，西到太行山，北达燕山，东邻渤海，南至黄河，是我国主要的经济、文化、政治中心。利用GRACE时变重力模型数据计算出的研究期间2002-08～2016-08 全球及华北平原水储量变化速率，从图1可以看出华北平原在研究期间内水储量整体上呈递减趋势。在空间上，华北平原西南部递减的速率最大，东北部递减速率相对较小，递减速率呈北缓南快的空间分布特征。

图1 2002年8月至2016年8月华北平原水储量变化速率分布

图2 给出了由GRACE 反演的华北平原2002-08～2016-08 水储量变化的时间序列图，以便对华北平原年际变化特征进行具体分析。如图所示，在2002-08～2016-08 期间，华北平原水储量变化总体趋势上呈现下降现象，下降速率为1.38 cm/a。2002年，我国北方遭遇旱灾，黄河遭遇特枯、京杭大运河山东济宁段断航50 d、河北辽河干流断流158 d[14]，受2002年极端干旱气候的影响,华北平原水储量在2002—2003 年处于较低值。2004年我国雨水丰富，部分中小河流发生较大洪水，所以在2002—2004年期间，华北水储量变为上升趋势，上升速率达到11.13 cm/a，2004 年华北平原水储量到达研究期间内的峰值；2004—2008 年间，华北平原水储量以2.87 cm/a的速率下降；2009—2013年，华北平原水储量表现为缓慢增长的趋势，增长速率为0.28 cm/a；2013—2016年，华北平原水储量减少速率最快为4.29 cm/a；2016 年7 月我国发生大面积的洪涝灾害，华北平原水储量在该期间内增幅较大。

图2 华北平原2002-08～2016-08水储量变化

由图2 可以看出，华北平原水储量变化存在明显的季节性变化，分别在冬季达到最大值，在夏季达到最小值。了解周年和半周年变化振幅和相位有利于了解华北平原水储量季节性变化特征，图3 计算出了由GRACE 解算出的2002 年8 月至2016 年8 月期间华北平原水储量变化的周年、半周年项。

图3 华北平原2002-08～2016-08水储量变化周年与半周年项

图3 表明华北平原南部属于黄河下游区域，雨水量较充分，因此华北平原南部周年和半周年振幅大于北方区域；周年振幅由华北平原西南向东北逐渐减小，最大值出现在黄河下游河南境内达到31.5 mm；华北平原周年相位表现为东部略大于西部，整个华北平原平均周年相位大约为325 d，说明华北平原水储量峰值一般出现在10～11 月附近，这一结论与图2 华北平原水储量变化时间序列所给出的信息相一致；与华北平原水储量变化的周年振幅相比，半周年振幅相对较小，其空间分布特征与周年振幅分布特征相似，也是由西南向东北方向递减，半周年振幅最大值达16.1 mm；华北平原半周年相位在空间上呈现西部大于东部的现象，说明华北平原半周年变化振幅东部快于西部，整个华北平原半周年相位平均约为130 d。

为了更好地揭示华北平原年内水储量变化的时空特征，本文利用GRACE 数据计算了华北平原2002-08～2016-08期间内各个月份的平均水储量，如图4所示，给出了华北平原研究期间内多年各月平均水储量的空间分布。

从月尺度上分析，有利于剔除时间上的复杂性和空间上的异质性，由图4可知从1月份至6月份华北平原水储量不断减少，6 月份至9 月份水储量不断增加，即华北平原水储量冬春之际处于不断亏损状态，夏季属于盈余补充状态。其中，1、2月份水储量空间分布差异不大，都是东部水储量高于西部，这是因为华北平原东邻渤海，东部区域气候相对湿润的原因。整个华北平原1 月份水储量略高于2 月水储量，水储量变化范围在9～17 mm 之间，整个华北平原3 月份水储量总值大概为0 左右，表明了华北平原3 月份水储量达到了1 a内的均值。4～7月华北平原水储量处于负值，是1 a 中最缺水的几个月份，且河南境内缺水表现最为严重，北京天津地区缺水情况相对较轻，其中4～6 月华北平原水储量空间分布呈现为西南地区低于东北地区，六月份整个华北平原水储量达到1 a当中最小值，其均值约为-50 mm；7～9月份华北平原水储量不断增加，由亏损状态转为盈余状态，9月份水储量达到了1 a之中最大值，且空间分布上表现为华北平原西南地域盈余最大达到34 mm，东北部相对较小，最小地区只有16 mm；10～11月份华北平原水储量空间分布情况与8月份恰好相反，西南地区水储量高于东北部。总体来说，华北平原春夏水储量较少，秋冬水储量较多；且冬春水储量整体不断减少，在6 月份达到谷值，夏季水储量不断补充在9月份达到峰值。

图4 华北平原水储量各月空间分布

2.3 GRACE与GLDAS和Mascon结果对比

本文采用GLDAS 水文模型中的Noah 数据和CSR提供的Mascon 数据与GRACE 观测结果进行比较分析。前文已经介绍了GLDAS可提供降雨、降雪、蒸发量、大气压、风速等多种数据，本文将选取其提供的地表积雪数据、土壤水含量及植被冠层水含量数据来模拟研究区域的水储量变化，其中土壤水含量分为4 层，包含了地下0～10 cm、10～40 cm、40～100 cm 及100～200 cm的土壤水。

图5给出了由GRACE、GLDAS Noah和Mascon 3种数据得到的华北平原水储量变化时间序列。由图可知，通过3 种数据计算得到的华北平原水储量长期变化趋势基本一致。经计算，GRACE 与GLDAS 反演结果相关系数为0.76；Mascon与GLDAS结果相关系数为0.71；GRACE 与Mascon 相关性最高为0.94，3 种数据反演结果呈强相关。从反演结果振幅上来看，GRACE反演的结果振幅大于Mascon 估算的结果，而Mascon估算的结果振幅又大于GLDAS水文模型的结果，这是因为GRACE 反演的结果相较于Mascon 进行了信号恢复，GRACE和Mascon反演结果又相较于GLDAS水文模型不仅包含了土壤表层水、地表积雪和植被冠层水，同时还包括了深层地下水和各种形式存在的水元素。与GRACE 和Mascon 结果相比，GLDAS 水文模型估算的结果在2011 年以后高于GRACE 和Mascon 结果，原因是近年来华北平原地下水遭到了严重开采，农业灌溉和生活用水主要来源于地下水，资料显示[15]2012年华北平原地下水开采量已达当年水资源利用量的6成以上，地下水又是GRACE反演结果中主要的组成部分，所以地下水的超采骤减会相应导致GRACE反演结果相较于平均水平骤缩，对GLDAS估算结果没有直接影响。GRACE、Mascon 和GLDAS 反演结果都显示出2004 年和2016 年水储量增幅最大，2002 年和2009年水储量变化偏小，这说明华北平原在2002年、2004 年、2009 年及2016 年发生了水文事件。查阅中国水旱灾害公报发现华北平原在2002 年和2009 年出现严重的旱灾；2004年和2016年7月全国出现大面积的洪涝灾害，华北平原受灾严重。

图5 GRACE、GLDAS和Mascon估算的华北平原水储量变化

2.4 与GPCC降水数据对比

区域水储量的变化是由区域内降水、蒸发、径流共同作用的，其水平衡方程表达式可写作为：

式中，ΔS为区域水储量变化量；P为降水量；E为蒸散发；R为地表、地下径流。我们可以看出降水量是影响区域水储量变化中的唯一输入量，因此降水量对区域水储量的影响至关重要。

因为GRACE 观测的某月区域水储量也包括了研究区域内该月份之前的降水部分，因此GRACE 观测结果会较降水结果出现一定的滞后现象，为了解决这一问题，使GRACE 观测结果与降水结果更具有直接相关性，本文对GRACE 得到的华北平原水储量变化时间序列做进一步处理：首先针对GRACE 中缺失的月份数据利用三次样条插值进行补充，然后将2002-08～2016-08 水储量变化数据用该月份数据减去前一个月份的数据，这样就得到了该月的水储量变化量。通过上述步骤计算得到的华北平原各个月份水储量变化量去和降水数据进行比较发现，由GRACE 计算出的月水储量变化量与GPCC月降水量峰值和谷值对应较好，峰值出现在1 a 的夏秋季；谷值出现在1 a的冬春季，说明华北平原水储量流量的季节性变化与降水的季节性规律有较好的一致性；水储量流量的变化走势与降水量的走势较为吻合，说明降水量是影响华北平原水储量变化的重要因素，但水储量的变化还受蒸散发和径流的影响，所以两者走势并不会完全相同。

3 结 语

本文利用2002-08～2016-08 的GRACE RL06 版本GSM 数据反演了华北平原水储量变化，同时结合了GLDAS水文模型数据、Mascon数据及GPCC降水数据对华北平原区域水储量时空变化进行了详细的分析，得到如下结论。

1）华北平原水储量变化呈现出明显的季节性和周年性变化，春夏水少、秋冬水多，水储量在1 a中6月份达到最小值，9月份达到最大值。研究期间内华北平原水储量以1.38 cm/a 的速度减少，其中2004—2008 年间，以2.87 cm/a 的速率减少；2009—2013 年，以0.28 cm/a的速率缓慢增长；2013—2016年，下降速率达4.29 cm/a。空间分布整体表现为由东北部向西南部减少速率不断加快。

2）GRACE 反演结果与GLDAS 和Mascon 计算结果具有较强的一致性，GRACE 观测的季节性变化振幅高于水文模型模拟的结果，与Mascon计算结果相关度最高。三者都能探测到华北平原在研究期间内发生了干旱和洪涝灾害。

3）GRACE 反演水储量流量变化与GPCC 降水数据具有较好的相关性，2016年华北平原水储量的增加和2002年华北平原的旱情，都在同时期降水数据上有着良好的反映，说明降雨是华北平原水储量变化的重要因素之一。