2020, Vol. 40
2. 南方科技大学理学院，深圳市学苑大道1088号，518055
GRACE卫星给出全球2002-05~2017-06共160组高精度的月时变重力场模型，为研究全球大尺度地表质量变化提供了重要的数据支撑。国外许多学者在时变重力场反演地表质量迁移的方法和理论基础上[1]，利用GRACE卫星数据分析了陆地水储量的变化情况[2-4], 国内很多学者[5-8]也开展了利用GRACE卫星监测华北平原、长江流域、黄河流域、黑河流域等地区陆地水储量变化的研究。
珠江流域是中国面积第三、径流量第二的大河流域，也是珠三角地区水资源供给的重要保障，因此，利用GRACE卫星研究该流域的陆地水储量变化显得尤为重要。本文利用CSR05、GFZ05、JPL05、ITSG-Grace2016、CSR06、GFZ06、JPL06及ITSG-Grace2018等8种GRACE时变重力场模型反演珠江流域2005-01~2012-12的陆地水储量变化，并结合GLDAS水文模型、降水数据及实测地下水数据等资料详细分析该地区陆地水储量的季节性特征和长期变化趋势。
1 反演理论与方法利用GRACE卫星时变重力场模型计算全球表面质量变化的公式为[1]：
| $\begin{array}{*{20}{c}} {\Delta \bar \sigma (\theta , \lambda ) = \frac{{a{\rho _{{\rm{ave }}}}}}{3}\sum\limits_{n = 0}^\infty {\sum\limits_{m = 0}^n {{{\bar P}_{{\rm{nm }}}}} } (\cos \theta )\frac{{2n + 1}}{{1 + {k_n}}}{w_n}}\\ {\left( {\Delta {{\bar C}_{nm}}\cos m\lambda + \Delta {{\bar S}_{nm}}\sin m\lambda } \right)} \end{array} $ | (1) |
式中，a为地球的赤道处半径，ρave为地球平均密度，ΔCnm、ΔSnm为规格化球谐系数的变化，由每个月份的球谐系数减去平均值获得，其他参数见文献[1]。
在利用式(1)进行计算时需要注意以下几点：1)由于GRACE卫星获取的C20项精度不高，可用SLR测定的C20项进行替换[9]；2)GRACE月时变模型没有1阶项，可采用Swenson方法等计算的地心改正项[10]；3)为削弱冰后回弹对反演结果的影响，还要对重力场模型进行冰川均衡改正(GIA)；4)GRACE时变重力场模型存在南北方向条带误差及高频噪声[11]，需采用Swenson方法等去相关滤波[11]和高斯滤波[12]等数据处理算法来抑制这2类误差的影响。其中，高斯平滑系数wn可以由递推公式获得[1]：
| $w_{0}=\frac{1}{2 \pi} ; w_{1}=\frac{1}{2 \pi}\left[\frac{1+\mathrm{e}^{-2 b}}{1-\mathrm{e}^{-2 b}}-\frac{1}{b}\right]\\ w_{n+1}=-\frac{2 n+1}{b} w_{n}+w_{n-1} $ | (2) |
式中，
首先对球谐系数变化量ΔCnm和ΔSnm进行Swenson方法去相关处理，并由式(2)计算得到高斯平滑系数wn，再代入式(1)即可获得全球质量变化的等效水柱高[8]:
| $\Delta h=\Delta \bar{\sigma} / \rho_{w} $ | (3) |
式中，ρw为水的平均密度(1 000 kg/m3)。
由于各时变模型扣除了潮汐影响和非潮汐的大气和海洋影响，经上述数据处理即可获得全球陆地水储量变化的等效水柱高格网值，再根据纬度余弦加权平均的方法计算区域质量变化的时间序列[7]。通过对时间序列进行拟合，即可得到区域陆地水储量变化1 a项、0.5 a项的振幅、相位和趋势项等信息，其公式为[13]：
| $\begin{array}{*{20}{c}} {f(t) = A + Bt + C\sin \omega t + D\cos \omega t + }\\ {E\sin 2\omega t + F\cos 2\omega t} \end{array} $ | (4) |
式中，A、B、C、D、E、F为待求常数，ω=2π。
2 研究数据本文研究数据采用2005-01~2012-12的CSR05、GFZ05、JPL05、ITSG-Grace2016、CSR06、GFZ06、JPL06和ITSG-Grace2018等8种GRACE时变重力场模型、GLDAS水文模型、降水数据及地下水测井监测数据，以便综合分析珠江流域的陆地水储量变化。
2.1 GRACE时变重力场时间序列选用ICGEM网站发布的8种GRACE时变重力场模型，其中3大官方机构2005版时变模型和ITSG-Grace2016模型利用GRACE2002版数据反演(称为A组模型)，2006版时变模型和ITSG-Grace2018模型利用GRACE2003版数据反演(称为B组模型)。模型的截断阶次为60，选取2005-01~2012-12共92个月的数据，其中2011-01、2011-06、2012-05、2012-10数据缺失，采用三次样条函数内插获取。
2.2 水文数据分别收集2005-01~2012-12的GLDAS水文模型、全国降雨数据和实测地下水数据等3种水文资料，各数据均为每月1个值。其中，GLDAS全球水文模型由NASA和NCEP共同建立，空间分辨率为1°×1°，包括NOAH、VIC、CLM、MOSAIC四种陆地表面模式[14]；降雨数据来自中国气象数据网(http://data.cma.cn/)的格网数据集，空间分辨率为0.5°×0.5°；地下水测井资料来自《中国地质环境监测地下水位年鉴》。
3 结果分析 3.1 不同GRACE时变重力场模型的比较为分析不同时变重力场模型的精度，对各模型的位系数不作任何滤波处理，直接利用A、B两组共8种时变重力场模型计算2005-01~2012-12全球质量变化趋势(等效水柱高)，结果见图 1。由图可知，各模型的结果都存在不同程度的条带误差，3大官方机构发布的2006版时变模型的条带误差要远小于2005版时变模型，且ITSG-Grace2018模型的条带误差也比ITSG-Grace2016模型小。整体来看，B组模型精度优于A组模型，且在所有模型中ITSG-Grace2018模型的精度最好。
|
图 1 不同时变重力场模型的精度比较 Fig. 1 Accuracy comparison of different temporal gravity field models |
以珠江流域为例，进一步分析8种时变模型之间的差异，图 2为8种重力场模型经Swenson去相关算法与300 km高斯滤波组合方法处理后计算出的珠江流域陆地水储量在空间上的变化趋势。由图可知，各时变模型反演的珠江流域陆地水储量变化的差异不大，但A组模型的反演结果比B组模型的反演结果略大。珠江流域大部分区域的陆地水储量处在增长状态，增速约为3~5 mm/a，其中北部地区可能受三峡水库蓄水的影响，陆地水储量增速较快；而西部地区即云南东部、贵州南部等降雨较少的区域，陆地水储量处于平衡或略呈下降的趋势。
|
图 2 珠江流域陆地水储量变化趋势的空间分布 Fig. 2 Spatial distribution of TWS trend in the PRB |
根据珠江流域的边界信息，利用纬度余弦加权平均的方法算得珠江流域质量变化的时间序列，结果见图 3。由图可知，CSR05、GFZ05、JPL05与ITSG-Grace2016模型时间序列的相关系数分别为0.992 6、0.988 6和0.983 9，而CSR06、GFZ06、JPL06与ITSG-Grace2018模型时间序列的相关系数分别为0.993 4、0.991 8和0.993 2，因此，与A组模型相比，B组模型之间的差异更小，说明不同机构利用GRACE2003版KBR(K-band ranging)等数据反演的时变模型具有更好的一致性。利用式(4)对各时变模型的时间序列进行拟合，获取1 a振幅等分量，结果见表 1，各时变模型的1 a振幅平均值约为65.4 mm，1 a相位平均值约为223.4°。对扣除1 a振幅和0.5 a振幅后的时间序列进行直线拟合，得到8种模型的趋势项分别为4.49 mm/a、4.17 mm/a、3.67 mm/a、4.68 mm/a、3.48 mm/a、3.71 mm/a、3.90 mm/a和4.00 mm/a。与B组模型相比，除了JPL05模型外，A组模型反演的珠江流域陆地水储量变化趋势均稍高一些。考虑到各时变重力场与GLDAS水文模式的相关系数均在0.87左右，且ITSG-Grace2018模型的条带误差小、精度高，因此后续分析主要基于ITSG-Grace2018时变重力场模型展开。
|
图 3 珠江流域陆地水储量变化时间序列 Fig. 3 Time series of TWS variations in the PRB |
|
表 1 珠江流域陆地水储量变化的1 a项、0.5 a项的振幅、相位 Tab. 1 Annual and semiannual amplitude and phase of TWS variations in the PRB |
将GLDAS水文模型的格网数据球谐展开并截断到60阶，采用与时变重力场模型相同的数据处理策略对其进行处理，并与GRACE时变重力场模型反演的结果进行对比。由图 4(a)可知，珠江流域的陆地水储量变化表现出一定的季节性，夏秋盈余，冬春亏损。ITSG-Grace2018模型与GLDAS水文模型确定的珠江流域的陆地水储量变化结果具有较强的相关性，两者的相关系数达到0.87，但GLDAS水文模型的振幅比ITSG-Grace2018模型的振幅略小，原因是GRACE重力场模型的反演结果反映的是珠江流域在综合因素影响下的质量变化情况，而GLDAS水文模型只包含了土壤水和积雪等表层物质的变化情况，反演结果符合2种数据的特性。进一步分析珠江流域的陆地水储量变化可知，在扣除ITSG-Grace2018模型的季节性信号后，对余下的信号进行直线拟合，结果见图 4(b)。结果表明，2005-01~2008-12珠江流域的陆地水以17.4 mm/a的速度增长，2008年末陆地水储量达到最大，2009-01~2009-12以35.4 mm/a的速度回落，2010-01~2012-12又以8.6 mm/a的速度慢慢回升，整体来看，2005-01~2012-12珠江流域陆地水储量以4.0 mm/a的速度增长。
|
图 4 ITSG-Grace2018与GLDAS反演的珠江流域陆地水储量变化 Fig. 4 TWS variations in the PRB derived from ITSG-Grace2018 and GLDAS |
降雨是水文过程的重要组成部分，也是中低纬度沿海地区陆地水主要的补给来源。图 5(a)中柱体为珠江流域的月降雨量情况，降雨的峰值出现在每年的6~8月，11~12月和次年的1~2月降雨量最少。整体来看，珠江流域的降雨量在时间上表现出一定的周期性，蓝色曲线为月降雨量扣除平均降雨量后每月的降雨异常量，红色曲线为ITSG-Grace2018模型反演的珠江流域陆地水储量的等效水柱高，两者具有较好的一致性。降雨较多的时期，珠江流域的陆地水储量也相应增加，降雨异常量与陆地水储量的相关系数为0.54，说明降雨为珠江流域陆地水储量变化的重要影响因素。通过式(4)对降雨异常量进行拟合，获得降雨异常量的1 a相位为175.0°。对比表 1中重力场模型的1 a相位信息可知，珠江流域的陆地水储量达到峰值要比降雨达到峰值的时间滞后约48 d，这是由于GRACE卫星获得的陆地水包括了观测月份之前的降雨部分，且地表水还需经过下渗、径流等水文过程才能转换成陆地水，因此出现陆地水储量变化比降雨滞后的现象[15]。
|
图 5 珠江流域陆地水与降雨的响应关系 Fig. 5 Response relationship between TWS and precipitation in the PRB |
按照春(3~5月)、夏(6~8月)、秋(9~11月)、冬(12月至次年1~2月)四季分别统计每个季度的陆地水储量和降雨异常量，结果见图 5(b)。珠江流域在2006年和2008年夏季的降雨量比其他年份夏季的降雨量要多，分别达59 cm和57 cm，但由图 5(a)可知，2006年夏季3个月的降雨量比较均匀，而2008年夏季的降雨量主要集中在6月。经查询，2008-06中旬华南地区遭遇持续强降雨，造成福建、江西、湖南、广东、广西、贵州等地区发生严重的洪涝灾害。ITSG-Grace2018模型的反演结果中2008年夏季陆地水储量比其他年份夏季的储水量要高，且季节性降雨异常量与ITSG-Grace2018模型反演的陆地水储量季节性变化的相关系数为0.65，说明GRACE卫星能够监测到珠江流域的洪水事件。
3.4 与地下水数据的对比分析地下水是陆地水的重要组成部分，分析区域地下水的变化情况也非常重要，但利用地下水测井观测等传统手段无法全面反映区域地下水的变化情况，可考虑采用流域水储量基本平衡方程[16]来获取珠江流域的地下水储量变化情况。由于珠江流域位于中低纬度地区，降雪量较少，本文将忽略该成分，则可根据公式ΔGW=ΔTWS－ΔSM获得地下水储量的变化情况，结果见图 6，其中绿色曲线为ITSG-Grace2018模型扣除GLDAS水文模型后获得的珠江流域地下水储量的变化情况。由图可知，地下水也具有明显的季节性变化特征，通过对其进行直线拟合后发现，珠江流域的地下水以4.3 mm/a的速度增长。
|
图 6 ITSG-Grace2018与GLDAS估算地下水变化与测井观测的地下水变化 Fig. 6 Groundwater storage variations derived from ITSG-Grace2018, GLDAS and monitoring wells observation |
利用珠江流域4个时段上较为连续的地下水测井数据对GRACE卫星的反演结果进行验证分析。地下水实测数据要与GRACE卫星的反演结果进行比较，需将地下水位变化量转化为水储量变化量，其计算公式为[6]：
| $\Delta \mathrm{GW}=\Delta H \times \mu $ | (5) |
式中，ΔH为地下水埋深变化值，μ为相应地区的给水度，根据地质情况，广东和广西地区的给水度分别为0.030和0.025。图 6中井1和井2均位于广州市白云区，井3和井4分别位于柳州市城中区和鱼峰区，由图可知，珠江流域的地下水测井水位在夏秋上升、冬春下降，具有一定的季节性，与ITSG-Grace2018模型反演的地下水变化情况基本一致。但珠江流域较为广阔，本文仅采用了4个局部地下水井的监测数据对GRACE卫星的反演结果进行验证，具有一定的局限性。
4 结语本文以CSR05、GFZ05、JPL05、ITSG-Grace2016、CSR06、GFZ06、JPL06及ITSG-Grace2018等8种GRACE时变重力场模型作为研究对象，分析不同版本时变模型的精度，最后选择ITSG-Grace2018模型、GLDAS水文模型、降雨数据、实测地下水数据等资料综合分析2005-01~2012-12珠江流域陆地水储量的变化情况，得出以下结论：
1) 3大官方机构发布的2006版GRACE时变模型的精度较2005版有显著改善，ITSG-Grace2018模型精度相对ITSG-Grace2016模型也有一定改善，且新版不同时变模型间的差异更小，说明不同机构利用GRACE2003版KBR等数据反演的时变模型具有更好的一致性。
2) 2005-01~2012-12，ITSG-Grace2018模型与GLDAS水文模型的反演结果具有较强的相关性，两者的相关系数达到0.87，季节性变化特征与降雨量、地下水观测数据都具有较好的一致性。
3) 珠江流域陆地水储量经历了先增长(2005-01~2008-12)-后下降(2009-01~2009-12)-再回升(2010-01~2012-12)3个状态，整体来看，2005-01~2012-12珠江流域陆地水储量以4.0 mm/a的速度增长。
致谢: ICGEM提供时变重力场模型数据，NASA和NCEP提供GLDAS水文数据，中国气象数据网提供降雨数据，冯伟博士提供GRACE_MATLAB_Toolbox软件，在此表示感谢。
| [1] |
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): 30 205-30 229 DOI:10.1029/98JB02844
(  0) |
| [2] |
Rodell M, Velicogna I, Famiglietti J S. Satellite-Based Estimates of Groundwater Depletion in India[J]. Nature, 2009, 460(7 258): 999-1 002
( 0) |
| [3] |
Tiwari V M, Wahr J, Swenson S. Dwindling Groundwater Resources in Northern India, from Satellite Gravity Observations[J]. Geophysical Research Letters, 2009, 36(18): 184-201
( 0) |
| [4] |
Famiglietti J S, Lo M, Ho S L, et al. Satellites Measure Recent Rates of Groundwater Depletion in California's Central Valley[J]. Geophysical Research Letters, 2011, 38(3)
( 0) |
| [5] |
Gong H L, Pan Y, Zheng L Q, et al. Long-Term Groundwater Storage Changes and Land Subsidence Development in the North China Plain (1971-2015)[J]. Hydrogeology Journal, 2018, 26(5): 1 417-1 427 DOI:10.1007/s10040-018-1768-4
( 0) |
| [6] |
冯伟, 王长青, 穆大鹏, 等. 基于GRACE的空间约束方法监测华北平原地下水储量变化[J]. 地球物理学报, 2017, 60(5): 1 630-1 642 (Feng Wei, Wang Changqing, Mu Dapeng, et al. Groundwater Storage Variations in the North China Plain from GRACE with Spatial Constraints[J]. Chinese Journal of Geophysics, 2017, 60(5): 1 630-1 642)
( 0) |
| [7] |
尼胜楠, 陈剑利, 李进, 等. 利用GRACE卫星时变重力场监测长江、黄河流域水储量变化[J]. 大地测量与地球动力学, 2014, 34(4): 49-55 (Ni Shengnan, Chen Jianli, Li Jin, et al. Terrestrial Water Storage Change in Yangtze and Yellow River Basins from Grace Time-Variable Gravity Measurements[J]. Journal of Geodesy and Geodynamics, 2014, 34(4): 49-55)
( 0) |
| [8] |
罗志才, 李琼, 钟波. 利用GRACE时变重力场反演黑河流域水储量变化[J]. 测绘学报, 2012, 41(5): 676-681 (Luo Zhicai, Li Qiong, Zhong Bo. Water Storage Variations in Heihe River Basin Recovered from GRACE Temporal Gravity Field[J]. Acta Geodaetica et Cartographica Sinica, 2012, 41(5): 676-681)
( 0) |
| [9] |
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) DOI:10.1029/2010JB000850
( 0) |
| [10] |
Swenson S, Chambers D, Wahr J. Estimating Geocenter Variations from a Combination of GRACE and Ocean Model Output[J]. Journal of Geophysical Research:Solid Earth, 2008, 113(B8)
( 0) |
| [11] |
Swenson S, Wahr J. Post-Processing Removal of Correlated Errors in GRACE Data[J]. Geophysical Research Letters, 2006, 33(8)
( 0) |
| [12] |
Jekeli C.Alternative Methods to Smooth the Earth's Gravity Field[D].Columbus: Ohio State University, 1981
( 0) |
| [13] |
Ran J, Ditmar P, Klees R, et al. Statistically Optimal Estimation of Greenland Ice Sheet Mass Variations from GRACE Monthly Solutions Using an Improved Mascon Approach[J]. Journal of Geodesy, 2018, 92(3): 299-319 DOI:10.1007/s00190-017-1063-5
( 0) |
| [14] |
Rodell M, Houser P R, Jambor U E A, et al. The Global Land Data Assimilation System[J]. Bulletin of the American Meteorological Society, 2004, 85(3): 381-394 DOI:10.1175/BAMS-85-3-381
( 0) |
| [15] |
冯伟, Lemoine J M, 钟敏, 等. 利用重力卫星GRACE监测亚马逊流域2002~2010年的陆地水变化[J]. 地球物理学报, 2012, 55(3): 814-821 (Feng Wei, Lemoine J M, Zhong Min, et al. Terrestrial Water Storage Changes in the Amazon Basin Measured by GRACE during 2002-2010[J]. Chinese Journal of Geophysics, 2012, 55(3): 814-821 DOI:10.6038/j.issn.0001-5733.2012.03.011)
( 0) |
| [16] |
Rodell M, Chen J L, Kato H, et al. Estimating Groundwater Storage Changes in the Mississippi River Basin (USA) Using GRACE[J]. Hydrogeology Journal, 2007, 15(1): 159-166
( 0) |
2. College of Science, Southern University of Science and Technology, 1088 Xueyuan Road, Shenzhen 518055, China