2018, Vol. 38
2. 中国测绘科学研究院，北京市莲花池西路28号, 100830;
3. 河北农业大学城乡建设学院，河北省保定市灵雨寺街289号，071001
GRACE重力卫星为研究大尺度的陆地水储量变化提供了可能。Wang等[1]较早利用2002-05~2010-05的GRACE RL-04产品研究了三峡水库蓄水引起的水储量变化；卢飞等[2]利用最新的GRACE RL-05产品对三峡库区水储量的变化开展了相应的研究，但他们对于信号泄漏误差的改正尚存在不足。Chen等[3-4]提出一种正向建模(forward modeling)方法，可以有效地减少GRACE数据后处理过程中信号衰减和泄漏造成的误差。本文主要侧重于对正向建模方法的实现，设计模拟算例验证正向建模对于信号衰减和泄漏的改正作用，并利用该方法恢复改正GRACE反演的三峡库区(图 1)水储量变化，最后与全球水文模型GLDAS、CPC以及全球降水数据GPCC(global precipitation climatology centre)进行对比分析，给出其时空分布特征解释。
|
图 1 三峡地区地形及库区、研究区域示意图 Fig. 1 The DEM of Three Gorges area and diagram of reservoir area, study area |
使用美国得克萨斯大学空间研究中心(Center for Space Research, University of Texas at Austin，CSR)提供的2002-07~2014-12共138个月的GRACE Level-2(Release-05)数据，期间缺失的数据采用样条插值方法补齐，该数据产品已将非潮汐的海洋、大气信号以及各种潮汐影响扣除[5]。对各个月份的球谐系数扣除2004~2009年的平均值，得到球谐系数变化(ΔC、ΔS)。考虑到GRACE的C20项精度不高，采用卫星激光测距(SLR)观测的C20进行替代[6]。根据重力场球谐系数反演以等效水高变化Δh表示的地表质量变化公式为[7]：
| $ \begin{array}{l} \Delta h\left( {\theta, \lambda } \right) = \frac{{a{\rho _e}}}{{3{\rho _w}}}\mathop \sum \limits_{l = 0}^{60} \frac{{2l + 1}}{{1 + {k_l}}}{W_l} \cdot \\ \mathop \sum \limits_{m = 0}^l \left( {\Delta {C_{lm}}\cos m\lambda + \Delta {S_{lm}}\sin m\lambda } \right){P_{lm}}\left( {\cos \theta } \right) \end{array} $ | (1) |
式中，θ和λ分别为计算点的余纬和经度，ρe为地球的平均密度，ρw为水的平均密度，kl为负荷Love数，Wl为滤波因子，Plm(cosθ)为l阶m次的完全规格化勒让德函数。由于受到卫星轨道误差、K波段测距误差等对重力场解算的影响，在使用式(1)计算等效水高时将球谐系数截断至60阶，以减小高阶项误差对计算结果的影响。同时，使用高斯滤波和改进的P3M15去相关[8-9]滤波削弱仍存在的噪声和南北条带误差等的影响。高斯滤波因子Wl可以通过式(2)迭代计算得到[7]：
| $\begin{array}{l} {W_0} = \frac{1}{{2\pi }}\\ {W_1} = \frac{1}{{2\pi }}\left[{\frac{{1 + {e^{-2\alpha }}}}{{1-{e^{-2\alpha }}}} - \frac{1}{\alpha }} \right]\\ {W_{n + 1}} = - \frac{{2n + 1}}{\alpha }{W_n} + {W_{n - 1}}\\ \alpha = \frac{{\ln 2}}{{\left( {1 - \cos \left( {{r_{1/2}}/{\rm{R}}} \right)} \right)}} \end{array}$ | (2) |
式中, r1/2为滤波半径，R为地球平均半径。
1.2 水文模式数据使用的水文模式数据来自美国宇航局哥达航空中心和美国国家环境预报中心共同建立的全球陆地资料同化系统(GLDAS)以及美国国家海洋和大气局气象预报中心提供的陆地同化数据集(CPC)。其中GLDAS采用NOAH、VIC和CLM等多种地表模型，本文采用基于NOAH地面模型的GLDAS，主要反映的是4层土壤水和雪水当量，空间分辨率为1°×1°，时间范围为2002-01~2014-12，每月一值[10]。CPC水文模式是根据全球观测到的降水分布而建立的，主要反映表层土壤水含量和积雪变化，其时间间隔为1个月，空间分辨率为0.5°×0.5°，时间跨度为1948-01至今[11]，本文选用2002~2014年的数据。水文模式数据的处理策略与GRACE相同。
2 正向建模方法在GRACE数据后处理过程中，由于对重力场球谐系数进行了截断及使用了高斯平滑滤波，会不同程度地削弱原始信号，引起信号的“泄漏”误差：一方面，使得研究区域内的信号泄漏到研究区域之外，造成研究区域内的信号衰减(leakage out)；另一方面，研究区域之外的信号会泄漏进研究区域，混入到研究区域内的信号中(leakage in)[12]。上述因素会导致研究区域内信号有一定程度的失真，对进一步的研究带来误差。本文模拟了三峡地区的重力场信号处理过程，以等效水高变化率的形式表示，如图 2所示。
|
图 2 不同条件下模拟信号处理结果 Fig. 2 The results under different condition |
图 2(a)为模拟的初始信号，模拟信号是28°N~32°N、107°E~111°E之间变化率为-80 mm/a、1°×1°的格网点，其他地区的变化率均设为0；图 2(b)为将模拟信号进行球谐展开并将球谐系数截断到60阶的结果；图 2(c)为球谐系数截断到60阶，并使用300 km高斯滤波的结果；图 2(d)为球谐系数截断到60阶，并使用500 km高斯滤波的结果。可以看出，GRACE数据后处理过程明显削弱了原始信号，且信号不同程度地泄漏到周围地区，这种现象随着滤波半径的增大而更加明显，当滤波半径为500 km时，信号几乎完全湮没在背景场中。为了获得研究区域内真实的质量变化时间序列，需要使用一定的方法来减小这种泄漏误差的影响。目前常用的改正方法主要有2种：1)利用水文模式数据估计“leakage in”和“leakage out”的影响，然后从GRACE时间序列中扣除“leakage in”影响，回加信号泄漏的“leakage out”影响，但由于水文模式数据存在较大的不确定性，使得对泄漏误差的估计也存在较大的误差[13]；2)对研究区域内的水文信号采用与GRACE相同的后处理过程，分析信号后处理前后的衰减比例，得到一个尺度因子，然后对于GRACE时间序列先扣除利用水文模式数据估计的“leakage in”影响，再利用尺度因子来恢复研究区域内的真实信号强度[14]。如果进一步假设研究区域内的信号是均一分布的，则可以不依赖于水文模式来估计尺度因子[15]；但如果区域内存在比较大的质量分布差异，如三峡库区蓄水会造成局部地区水储量产生明显变化，那么这种方法也存在较大的不确定性[16]。
针对泄漏误差的改正，Chen等[3-4]设计了一种正向建模方法以恢复原始的真实信号。
1) 对于全球1°×1°的每个格网点，令其初始值等于GRACE计算得到的等效水高值(或等效水高变化)，计算过程采用P3M15去相关滤波和半径为500 km的高斯滤波，结果记为FM，并将其作为FM0的初值。
2) 将FM0进行球谐系数展开，截断到60阶，0阶项和1阶项均设为0，然后采用半径为500 km的高斯滤波，计算全球陆地1°×1°格网点的等效水高值(或等效水高变化)，结果记为FM1。
3) 将每一个格网点上GRACE计算的等效水高值(或等效水高变化)与FM1之间的差值回加到FM0上，形成新的初值，然后重复步骤2)。在回加过程中可以使用1.2的缩放因子来加快收敛速度。
4) 当FM1与GRACE计算的等效水高值(或等效水高变化)之间的差值小于特定的阈值或达到一定的迭代次数时，结束计算。
通过以上建模过程可以使得数据逐渐逼近真值。需要指出的是，由于使用的去相关滤波器是非线性的，无法对其进行量化评估，所以在正向建模过程中只使用高斯滤波。图 3为分别使用正向建模方法和尺度因子方法对三峡地区模拟数据的500 km高斯滤波结果进行改正前后的对比图。
|
图 3 信号经正向建模或尺度因子改正的结果对比 Fig. 3 The correction based on forward modeling and scale factor |
图 3(a)为对图 2(a)中模拟信号进行500 km高斯滤波之后的结果；图 3(b)为正向建模方法恢复后的结果，即最后一次迭代时的FM0，可以看到，经正向建模后，泄漏明显收敛，对信号的恢复效果十分显著；图 3(c)为正向建模方法恢复后的信号经过500 km滤波之后的结果，即最后一次迭代时的FM1，可以看出，图 3(a)和图 3(c)的差值很小，说明恢复的信号在经过相同的滤波处理后与模拟信号的滤波结果基本相同，这时可以认为经正向建模方法恢复的信号是可信的；图 3(d)为使用了基于均一假设估计的尺度因子恢复后的信号结果[15]，可以看到，尺度因子法虽然对信号振幅的衰减有一定恢复作用，但未解决信号的泄漏问题。经过计算，三峡地区采用正向建模方法之前结果与真实信号不符值的RMS为19 mm/a，使用尺度因子方法和正向建模方法改正之后的RMS分别为17 mm/a、7 mm/a，说明正向建模所恢复的真实信号(图 3(b))明显减小了泄漏效应带来的误差。
3 三峡库区水储量反演 3.1 水储量变化时间序列分析利用正向建模方法改正GRACE数据结果，并结合水文模式数据CPC、GLDAS和GPCC全球降水数据，研究三峡库区的水储量变化，研究区域为图 1中黑色矩形部分。其中GRACE反映的是库区总体的水储量变化，主要包括土壤水、积雪水、地表径流和地下水，而水文模型只包含了土壤水和积雪水。反演结果的时间序列如图 4、5所示。
|
图 4 正向建模改正前后GRACE时间序列 Fig. 4 GRACE TWS time series before and after forward modeling |
|
图 5 水文模式与降水数据时间序列 Fig. 5 Time series of hydrological model and precipitation |
由图 4可以发现，GRACE反演的库区总水储量和水文模式数据的地表水储量均表现出明显的季节性变化。利用正向建模改正后，GRACE时间序列振幅的衰减得到一定程度的恢复。经过计算，GRACE在经过正向建模改正前后的年周期项振幅分别为7.0±2.4 cm和7.2±2.7 cm。其中第1次蓄水前等效水高的振幅恢复效果最为明显，详细的分析见后文(后文中，GRACE结果均指经过正向建模改正后的结果)。水文模式数据(CPC与GLDAS的均值)的年周期项振幅为6.8±2.1 cm，水文模式数据主要反映的是土壤水分的季节性变化，这种变化主要受降水影响。结合图 5中GPCC降水的时间序列也可以说明这种相关关系，即夏季降水多，土壤水增加；冬季降水少，土壤水也相应减少。根据水量平衡方程(GRACE前后月份之差=降雨-蒸发-径流(包含水库蓄放水))可知，影响GRACE结果的因素更多一些，其结果除了反映土壤水分的季节性变化外，还反映库区蓄放水、蒸散发等的季节性变化，所以其振幅比水文模式的结果稍大[17]。
从图 4、5中看到，在年际变化趋势上，2006年和2011年的降水量较小，其月平均降水分别较前一年减少20.3 cm和8.4 cm。GRACE和水文模式数据的时间序列也均呈现出比较明显的减小趋势，但水文模式数据的减小趋势更明显，这也说明GRACE结果受水库蓄放水等多种因素的共同影响，降水只是其中重要的因素之一。2007年的月平均降水量较前年增加了32.3 cm，增加趋势十分明显，但GRACE和水文模式数据均没有明显的变化趋势。推测这一时期达到了区域最大有效水储量，因而水储量不会一直增加，盈余部分通过地表径流等流出区域外[18]。这一现象容易诱发洪水灾害，但根据中华人民共和国水利部发布的《2007中国水旱灾害公报》(http://www.mwr.gov.cn/)来看，2007年间长江流域并未有大型洪涝灾害发生，说明三峡工程在汛期削峰上发挥了一定的作用。
通过对研究时段内的数据进行一阶线性拟合发现，GLDAS反演的等效水高变化速率为-0.1±0.4 mm/a，CPC反演的等效水高变化速率为0.2±0.6 mm/a，而GRACE反演的等效水高变化为7.2±0.6 mm/a。拟合结果表明，该地区土壤水分的变化不足以引起陆地水储量在长期趋势上产生较大的变化。三峡水库从2003年开始蓄水至今经历了3个阶段：2003-06开始第1次蓄水，水位从约70 m上升到135 m；2006-10蓄水至156 m；2010-10蓄水至175 m。据此分析，该地区水储量变化主要是由三峡蓄水引起的。为了分析蓄水量的影响，以3次蓄水时间为节点，计算3个时段的等效水高平均值。结果显示，3次蓄水引起的等效水高变化分别为60.2 mm、28.3 mm、20.5 mm。考虑本文监测区域面积为19.9万km3，将上述等效水高变化换算为库容量变化，并与利用水位观测估计的库容量变化进行对比。其中水位观测结果可以通过以下3个步骤得到：1)利用水文站获得某一片区域的水位变化；2)利用遥感影像获得某一片区域的水域面积变化；3)水位变化乘以面积变化得到库容量变化。对比结果如表 1(单位km3)所示。
|
表 1 三峡库区不同蓄水阶段相对蓄水前的库容量变化的估计 Tab. 1 Estimation of the reservoir capacity change before the impoundment of the Three Gorges reservoir in different storage stages |
由表 1可以看出，GRACE监测结果与水文站观测估计的结果比较接近，存在差异的原因是:1)岸边存在斜坡等，水位观测估计所得到的结果也有较大的不确定性；2)GRACE结果不只包含因库区蓄水引起的水储量变化，如2006年降水减小，会引起地表土壤水储量的变化(图 5)，GRACE的监测结果中也包含这些影响；3)正向建模尚不能完全改正泄漏误差的，仍有部分残余。本文结果较卢飞等[2](3次蓄水引起的等效水高总变化分别为52 mm、18 mm和7 mm)和Wang等[1](前2次蓄水引起的库容量变化分别为10.84 km3和4.54 km3)反演的结果偏大，其主要原因是本文使用正向建模改正泄漏误差，削弱了振幅衰减现象，使得结果与水位观测估计得到的库容量变化更为接近，说明利用正向建模改正后的结果更可靠。
3.2 正向建模改正前后水储量空间分布对比§3.1提到利用正向建模改正后的GRACE时间序列振幅衰减在第1次蓄水(2003-06)前的恢复最为明显，本节着重分析这一时间段前后水储量变化率的空间分布、滤波对信号的湮没和正向建模对信号空间分布的恢复作用。GRACE观测数据中噪声较高，为了提高信噪比，得到较为可靠的信号空间分布，就需要使用足够长时段的数据[19]。将第1、2次蓄水引起的等效水高变化一起分析，时间跨度为2003~2009年。图 6、7分别给出了这期间正向建模改正前后三峡地区的水储量变化速度的平均值。
|
图 6 2003~2009年正向建模改正前水储量变化速度 Fig. 6 Terrestrial water variation before forward modeling during 2003~2009 |
|
图 7 2003~2009年正向建模改正后水储量变化速度 Fig. 7 Terrestrial water variation after forward modeling during 2003~2009 |
由图 6可以看到，三峡地区的东南部分以亏损为主，西北部分以盈余为主。使用500 km的高斯滤波之后，呈现出从西北向东南依次递减的阶梯状分布，同时由于滤波半径比较大，水库蓄水引起的水储量变化几乎完全湮没在了库区周边地区向研究区域内的泄漏信号之中，同时水库蓄水的正信号也泄漏出研究区域，如东南方向的亏损信号受这一影响，信号振幅衰减到了0 mm/a附近。GRACE从2002-04才开始提供数据，从水库蓄水前到第1次蓄水这一时段内数据较短，且由于轨道调整等原因，前几个月份的数据精度较差，信噪比低[20]。再加上使用了较大半径的高斯滤波，从时间序列中就很难发现第1次蓄水引起的等效水高变化，基本被湮没在了背景场中。
由图 7可以发现，正向建模改正之后泄漏效应明显收敛，可以清楚地看到库区蓄水引起的水储量变化，库区内的等效水高变化大致为8 mm/a，最大亏损集中在东南区域，达到-4 mm/a左右，西部区域的盈余信号也得到一定程度的恢复，水储量变化大致为4.5 mm/a。但蓄水引起的等效水高变化并未完全收敛于库区，在库区与东南地区亏损信号和西北地区盈余信号交界处仍有部分泄漏现象，其主要原因是：1)GRACE的分辨率大致为300 km[7]，不能使信号精细地收敛于库区之内；2)正向建模时用到了球谐展开，由于球谐函数的性质，在展开时也有信号泄漏，这一点可以从图 2(b)中看出。
4 结语本文使用正向建模方法改正GRACE反演的三峡库区陆地水储量变化的泄漏误差，有效改善了反演结果，并将其与GLDAS和CPC水文模式结果、GPCC降水数据以及水位观测估计的库区容量变化进行比较分析，结果表明：
1) 通过对模拟数据的分析计算发现，泄漏误差改正的量级约为1 cm，因此在利用GRACE反演陆地水储量变化时这一因素不可忽视。采用正向建模方法改正后，被湮没的信号得到有效恢复，一定程度上改善了信号的空间分辨效果。
2) 从三峡库区陆地水储量变化的时间序列上看，正向建模改正前后GRACE年周期项振幅分别为7.0±2.4 cm和7.2±2.7 cm，改正后信号泄漏有所减小。GRACE和水文模式数据均表现出明显的季节性特征，但GRACE的周年项振幅略大于水文模式的周年项振幅(6.8±2.1 cm)。水文模式数据主要由降水数据驱动，而GRACE结果除受降水的影响外，还受库区蓄放水、蒸散发等因素影响，故其振幅稍大。利用GRACE监测得到三峡库区3次蓄水引起的等效水高变化分别为60.2 mm、28.3 mm、20.5 mm，这一结果与利用水位观测估计得到的库区容量变化具有较好的一致性。
3) 对2003~2009年正向建模改正前后GRACE反演得到的三峡地区等效水高变化空间分布进行分析比较发现，改正前信号泄漏现象比较严重，库区蓄水信号几乎被湮没在背景场中；改正后，蓄水信号明显收敛于三峡库区，效果十分明显。但由于受到GRACE分辨率和球谐函数性质等因素的限制，仍有部分信号泄漏现象。
| [1] |
Wang X W, Linage C D, Famiglietti J, et al. Gravity Recovery and Climate Experiment (GRACE) Detection of Water Storage Changes in the Three Gorges Reservoir of China and Comparison with in Situ Measurements[J]. Water Resources Research, 2011, 47(12): 1091-1096
(  0) |
| [2] |
卢飞, 游为, 范东明, 等. 由GRACE RL05数据反演近10年中国大陆水储量及海水质量变化[J]. 测绘学报, 2015, 44(2): 160-167 (Lu Fei, You Wei, Fan Dongming, et al. Chinese Continental Water Storage and Ocean Water Mass Variations Analysis in Recent Ten Years Based on GRACE RL05 Data[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(2): 160-167)
( 0) |
| [3] |
Chen J L, Li J, Zhang Z Z, et al. Long-Term Groundwater Variations in Northwest India from Satellite Gravity Measurements[J]. Global & Planetary Change, 2014, 116(3): 130-138
( 0) |
| [4] |
Chen J L, Wilson C R, Li J, et al. Reducing Leakage Error in GRACE-Observed Long-Term Ice Mass Change: A Case Study in West Antarctica[J]. Journal of Geodesy, 2015, 89(9): 925-940 DOI:10.1007/s00190-015-0824-2
( 0) |
| [5] |
UTCSR. UTCSR Level-2 Processing Standards Document for Level-2 Product Release 0005[Z]. Center for Space Research, the University of Texas at Austin, 2013
( 0) |
| [6] |
Cheng M K, Ries J C, Tapley B D. Variations of the Earth's Figure Axis from Satellite Laser Ranging and GRACE[J]. Journal of Geophysical Research :Solid Earth, 2011, 116(B1)
( 0) |
| [7] |
Wahr J, Molenaar M, Bryan F. Time Variability of the Earth's Gravity Field: Hydrological and Oceanic Effects and Their Possible Detection Using GRACE[J]. Journal of Geophysical Research:Solid Earth, 1998, 103(B12): 30205-30230 DOI:10.1029/98JB02844
( 0) |
| [8] |
Swenson S, Wahr J. Post-Processing Removal of Correlated Errors in GRACE Data[J]. Geophysical Research Letters, 2006, 33(8)
( 0) |
| [9] |
詹金刚, 王勇, 郝晓光. GRACE时变重力位系数误差的改进去相关算法[J]. 测绘学报, 2011, 40(4): 442-446 (Zhan Jingang, Wang Yong, Hao Xiaoguang. Improved Method for Removal of Correlated Errors in GRACE Data[J]. Acta Geodaetica et Cartographica Sinica, 2011, 40(4): 442-446)
( 0) |
| [10] |
Rodell M, Houser P R, Jambor U, et al. The Global Land Data Assimilation System[J]. Bulletin of the American Meteorological Society, 2004, 84(3): 381-394
( 0) |
| [11] |
Fan Y, Dool H V D. Climate Prediction Center Global Monthly Soil Moisture Data Set at 0.5 Resolution for 1948 to Present[J]. Journal of Geophysical Research: Atmospheres, 2004, 109(D10)
( 0) |
| [12] |
Klees R, Zapreeva E A, Winsemius H C, et al. The Bias in GRACE Estimates of Continental Water Storage Variations[J]. Hydrology & Earth System Sciences, 2007, 11(4): 1227-1241
( 0) |
| [13] |
Longuevergne L, Scanlon B R, Wilson C R. GRACE Hydrological Estimates for Small Basins: Evaluating Processing Approaches on the High Plains Aquifer, USA[J]. Water Resources Research, 2010, 46(11): 6291-6297
( 0) |
| [14] |
Feng W, Zhong M, Lemoine J M, et al. Evaluation ofGroundwater Depletion in North China Using the Gravity Recovery and Climate Experiment (GRACE) Data and Ground-Based Measurements[J]. Water Resources Research, 2013, 49(4): 2110-2118 DOI:10.1002/wrcr.20192
( 0) |
| [15] |
Feng W, Lemoine J M, Zhong M, et al. Mass-Induced Sea Level Variations in the Red Sea from GRACE, Steric-Corrected Altimetry, in Situ Bottom Pressure Records, and Hydrographic Observations[J]. Journal of Geodynamics, 2014, 78(3): 1-7
( 0) |
| [16] |
冯伟, 王长青, 穆大鹏, 等. 基于GRACE的空间约束方法监测华北平原地下水储量变化[J]. 地球物理学报, 2017, 60(5): 1630-1642 (Feng Wei, Wang Changqing, Mu Dapeng, et al. Groundwater Storage Variations in the North China Plain from GRACE with Spatial Constraints[J]. Chinese J Geophys, 2017, 60(5): 1630-1642)
( 0) |
| [17] |
Pan Y, Zhang C, Gong H L, et al. Detection of Human-Induced Evapotranspiration Using GRACE Satellite Observations in the Haihe River Basin of China[J]. Geophysical Research Letters, 2017, 44(1): 190-199 DOI:10.1002/2016GL071287
( 0) |
| [18] |
Reager J T, Famiglietti J S. Global Terrestrial Water Storage Capacity and Flood Potential Using GRACE[J]. Radiology, 2009, 36(23): 195-215
( 0) |
| [19] |
Liu Y C, Hwang C, Han J, et al. Sediment-Mass Accumulation Rate and Variability in the East China Sea Detected by GRACE[J]. Remote Sensing, 2016, 8(9): 777 DOI:10.3390/rs8090777
( 0) |
| [20] |
钟敏, 段建宾, 许厚泽, 等. 利用卫星重力观测研究近5年中国陆地水量中长空间尺度的变化趋势[J]. 科学通报, 2009, 54(9): 1290-1294 (Zhong Min, Duan Jianbin, Xu Houze, et al. Trend of China Land Water Storage Redistribution at Medi- and Large-Spatial Scales in Recent Five Years by Satellite Gravity Observation[J]. Chinese Sci Bull, 2009, 54(9): 1290-1294)
( 0) |
2. Chinese Academy of Surveying and Mapping, 28 West-Lianhuachi Road, Beijing 100830, China;
3. Urban and Rural Construction Institute, Hebei Agricultural University, 289 Lingyusi Street, Baoding 071001, China