2018, Vol. 38
2. 北海道大学理学研究院，日本札幌市北八西5街，0600808
印度板块和欧亚大陆碰撞形成世界第三极——青藏高原[1]，活跃的地质构造、独特的地理气候环境给高原带来复杂的变化特征，包括构造运动、湖泊、冰川、土壤水、冻土和地下水变化等[2-4]。发生在青藏高原的这些地球物理现象都会带来质量迁移而引起重力变化，因此重力观测提供青藏高原质量变化的整体约束。美、德联合研制的GRACE重力卫星于2002-03发射升空，可用来监测全球大尺度时变重力[5-6]。
关于青藏高原质量变化研究，其数据主要包括遥感、卫星测高、GRACE卫星重力、GPS和水文模型等。Jacob等[7]利用GRACE观测给出2003~2010年青藏高原内部区域质量变化为7 Gt/a，错误地归结为这部分质量增加来自冰川增长，与地面实际观测不符。青藏高原区域冰川分为3类，南部冰川受夏季印度洋带来的强降水影响，帕米尔高原冰川受冬季西风的影响，青藏高原内部冰川受大陆气候控制。地面长期观测资料表明，青藏高原南端和喀拉昆仑山的冰川融化最快，青藏高原内部的冰川融化很慢，不会出现如此迅速的质量变化[8]。Zhang等[9]利用2003~2009年ICESat卫星观测数据得到青藏高原内部湖泊处于快速增长阶段，其质量增长率为8.06 Gt/a，认为能够解释青藏高原区域的质量增加。Yi等[3]在对亚洲高山冰川变化进行评估时发现，青藏高原内部区域质量增长速率为30 Gt/a，意味着湖泊水位的增长并不足以解释青藏高原区域明显的质量增加信号。Liang等[10]给出的现今青藏高原地壳运动三维速度场显示，印度板块向北挤压速率为4 cm/a，青藏高原整体隆升速率为2 mm/a，南部上升速率约为3~4 mm/a，意味着印度板块挤压欧亚板块带来的构造运动可能是导致青藏高原重力增加的一个因素。Yi等[11]联合GRACE和GPS研究表明，青藏高原东缘深部地壳流可能带来0.32 μGal/a的重力增加，但青藏高原区域冻土和地下水状况仍然缺乏大范围的实地观测资料。一般来说，一个更高的湖面需要一个更高的地下水位来维持，青藏高原内部大部分河流是内流湖，水量平衡测量表明渗流起很大作用[12]，Xiang等[4]认为大量残余信号全部来自地下水增长，其估计的地下水增长速率为18.6 Gt/a。由此可见，地下水估计精度取决于已知变量的精度，尤其是水循环系统中的重要成分——土壤水储量的精度。为获取可靠的青藏高原区域质量变化，探究各种影响要素的占比及来源，深刻认识青藏高原区域的动力学变化，有必要利用新的、更准确的数据进行分析。
现今广泛使用的水文模型GLDAS(global land data assimilation system)是一组地表水和能量场的数据集，由GSFC和NCEP联合研制[13]，已有相关实验表明GLDAS模型在估算土壤含水量上的有效性[14-15]。GLDAS水文模型已从GLDAS_1更新至GLDAS_2，GLDAS_2有GLDAS_2.0和GLDAS_2.1两个版本。
Bi等[16]将观测站数据与GLDAS_2.0、GLDAS_1模型数据比较表明，GLDAS_2.0相对于GLDAS_1水文模型数据在青藏高原区域没有显著的精度提升。由于GLDAS_2.0未提供2010年之后的数据，且新的GLDAS_2.1数据相较于GLDAS_2.0存在一定的改进，为探究GLDAS_2.1相对于GLDAS_1的精度提升，以青藏高原区域为例，分别利用GRACE数据、GLDAS_2.1版水文数据及GLDAS_1版水文数据计算青藏高原区域水储量变化，并对比分析3种数据得到的趋势、年振幅、年相位及相关性，说明新的GLDAS水文模型在青藏高原区域的精度提升。
1 数据处理 1.1 GRACE月均等效水高形式的地表质量变化数据由JPL提供[17-19]，该数据采用GRACE RL-05数据集，空间分辨率截断到60阶，融合CSR、GFZ和JPL 3种球谐系数集，用SLR结果替换二阶球谐系数[20]，用Swenson等[21]地心改正模型进行一阶项改正，并进行去相关滤波和300 km高斯平滑。此外，该数据采用三维冰后回弹(PGR)模型Geruo13[22]改正GRACE观测值中的相关信号。
1.2 GLDAS采用GLDAS NOAH水文模型GLDAS_NOAH_M.001(以下称GLDAS_1)及最新GLDAS_NOAH_M.2.1(以下称GLDAS_2.1)月均1°×1°格网的产品，数据范围为180° W ~180° E、60° S~90° N。通过水文模型中的4层土壤(0~0.1 m、0.1~0.4 m、0.4~1 m、1~2 m)湿度和积雪等数据来计算青藏高原区域陆地水储量变化。GLDAS_1采用的驱动数据不连续，2000年采用NOAA的全球同化系统(global data assimilation system，GDAS)大气场分析数据，2001年至今综合采用NOAA/GDAS大气数据、NOAA气候预测中心降水融合数据(climate prediction center merged analysis of precipitation，CMAP)以及美国空军气象局农业气象模拟系统(air force weather agency’s agricultural meteorological modeling system，ARGMET)生成的长短波辐射数据等。GLDAS_2.1融合GDAS、全球降水气候计划(global precipitation climatology project, GPCP)和ARGMET模型，并以实测数据作为约束，尽可能保证与长期气候学上的一致性，GLDAS_2.1以替换掉GLDAS_1.0为目标[23]。
GRACE数据进行去相关滤波及300 km高斯平滑，为将水文模型数据与GRACE数据进行比较，本文将GLDAS_1和GLDAS_2.1数据展开并截断至60阶球谐系数，再进行300 km高斯平滑。为削弱截断和平滑对结果的影响，将得到的各区域质量变化趋势结果乘以一个比例因子(s)，该比例因子通过以下方式获得：构造与GLDAS数据相似的格网点数据f(x, y)，需要计算的区域赋值1，其他为0；f(x, y)展开成球谐系数，并截断至60阶；顾及300 km高斯平滑计算g(x, y)，比例因子s=f(x, y)/g(x, y)。
2 结果与分析 2.1 质量变化趋势、年振幅和年相位利用2003-01~2016-06共162个月的GRACE、GLDAS_1、GLDAS-2.1数据，根据式(1)，通过最小二乘拟合计算青藏高原区域整体质量变化及水储量变化速度的空间分布、区域平均变化速度，并分析比较3种结果的振幅、相位及相关性。考虑到GRACE与GLDAS数据时间不一致，将GRACE数据插值到与GLDAS数据相同的时间序列:
| $ \begin{array}{l} y = a + bx + c\cos \left( {\frac{{2{\rm{ \mathsf{ π} }}}}{{{T_1}}}x + {\varphi _1}} \right) + \\ \;\;\;\;\;d\cos \left( {\frac{{2{\rm{ \mathsf{ π} }}}}{{{T_2}}}x + {\varphi _2}} \right) + \varepsilon \end{array} $ | (1) |
式中，a为常数项，b为趋势项，c、d为年振幅及半年振幅，T1、T2代表年周期及半年周期，φ1、φ2为初始相位，ε代表残差。
图 1(a)、(b)、(c)分别为由GRACE、GLDAS_1及GLDAS_2.1数据得到的研究区域以等效水高形式表示的质量变化速度。GRACE结果表明，青藏高原南部区域存在明显的负信号，主要由印度北部地下水亏损[24-25]和喜马拉雅冰川融化导致[2-3]。青藏高原内部的正信号，由于信号源复杂且缺乏足够的观测资料，不能确定所有信号的构成[3, 7, 9]，但可确定并非由冰川增长导致。研究表明，青藏高原内部冰川在退缩，尽管比周边区域慢得多[8, 26-27]。
|
图 1 2003-01~2016-06三种数据结果 Fig. 1 Results of three kind of data from January 2003 to June 2016 |
图 1(d)、(e)、(f)分别为由GRACE、GLDAS_1及GLDAS_2.1数据得到的研究区域质量变化时间序列的年相位空间分布，图 1(g)、(h)、(i)分别为由GRACE、GLDAS_1及GLDAS_2.1数据得到的研究区域质量变化时间序列的年振幅空间分布。整体而言，GLDAS_2.1与GRACE的年振幅、年相位符合较好或相当，即陆地水储量变化主导青藏高原区域质量季节性变化；GLDAS_1年周期结果空间差异较大，尤其在北部区域存在明显的圆形区域极值，而GLDAS_2.1则不存在。该极值很可能与模型相关，GLDAS_1只是将不同时段的数据集组合在一起，GLDAS_2.1尽可能保证在长期气候学上的一致性。此外，ARGMET短波数据某些年份在北半球会出现明显异常，GLDAS_2.1数据对此异常已进行相关改进。为更好地分析新的GLDAS_2.1模型相较于GLDAS_1模型的精度提升，将研究区域分为A、B、C、D、E、F等6个子区域，如图 2所示。
|
图 2 研究区域划分(B、C、D、E为青藏高原主体，A为帕米尔高原，F为印度北部) Fig. 2 Division of research area |
为比较GLDAS_2.1、GLDAS_1波动，对6个子区域的时间序列分别进行去趋势处理，顾及年周期、半年周期项进行最小二乘拟合，3种数据在6个子区域的去趋势时间序列、年振幅、年相位分别如图 3、图 4、图 5所示。可见，在整个青藏高原区域，GLDAS_2.1与GRACE的年振幅较为一致，但GLDAS_1与GRACE结果在B、C、D、E和F区域存在较大差异。在B区域，GLDAS_1年振幅远大于GRACE结果，D、E和F区域年振幅小于GRACE结果；在整个青藏高原区域，GLDAS_2.1与GRACE年相位较一致，GLDAS_1仅在B区域存在异常，与年相位空间分布类似。年振幅、年相位结果表明，在整个青藏高原区域GLDAS_2.1与GRACE结果较一致；在青藏高原中部、西南部及南部区域GLDAS_1与GRACE年振幅的一致性较差，在青藏高原中部区域GLDAS_1年相位存在异常，即GLDAS_2.1在青藏高原中部、西南部及南部区域相较于GLDAS_1水文模型数据有一定的精度提升。
|
图 3 6个区域质量变化时间序列(从上往下分别为区域A、B、C、D、E、F) Fig. 3 Mass change series of six regions |
|
图 4 3种数据的年振幅 Fig. 4 Year phase of three kinds of data |
|
图 5 3种数据的年相位 Fig. 5 Year amplitude of three kinds of data |
GRACE、GLDAS_1和GLDAS_2.1数据在各子区域的相关性如图 6和表 1所示。可见，GRACE与GLDAS_2.1的水储量时间序列具有强相关性，尤其在D、E、F区域；与GLDAS_2.1相似，GLDAS_1与GRACE在D、E、F区域也具有较高的相关性，但不及GLDAS_2.1，但在B区域GLDAS_1与GRACE相关性不大；GLDAS_1和GLDAS_2.1在B、C区域(青藏高原北部)相关性较差。结果表明，GLDAS_2.1能很好地反映青藏高原区域季节性质量变化，尤其在青藏高原中部区域、东北及东南区域；GLDAS_1数据基本能体现，但在中北部存在异常；GLDAS_2.1较GLDAS_1能更好地反映研究区域质量的季节性变化。
|
图 6 相关性空间分布 Fig. 6 Spatial distribution of correlation |
|
表 1 3种数据的相关性结果 Tab. 1 Correlation of three kinds of data |
对6个子区域A、B、C、D、E、F数据进行加权平均，得到各区域GRACE、GLDAS_1、GLDAS_2.1的质量变化速率，如图 7和表 2(2倍中误差)所示。GLDAS只包括土壤水储量变化，而GRACE反映的是整体的质量变化，还包括冰川、湖泊等的贡献，二者不一定具有相同的趋势性变化，GRACE扣除GLDAS趋势得残差变化速率。GRACE在青藏高原区域的整体质量变化与前人研究结果一致，即青藏高原南部区域为明显的负信号[24-25]，内部为正信号[3, 7]。GLDAS_2.1比GLDAS_1具有更高的精度，其提升主要体现在周期性和相关性两方面。在6个子区域上两个模型的长趋势存在较大差异，主要体现在GLDAS_1的差异被引入相关改正模型中，如估计冰川融化速率[3]、地下水变化[4]。在青藏高原主体区域，GLDAS_2.1比GLDAS_1结果小10.4 Gt/a，若Xiang等[4]采用GLDAS_2.1估算地下水变化，其结果18.6 Gt/a变为29.0 Gt/a，即GLDAS_1对结果产生57%的影响。
|
图 7 各区域质量变化趋势 Fig. 7 Mass change trend of different regions |
|
|
表 2 各区域质量估计差异 Tab. 2 Regional mass estimation and calculation of mutual deviation |
利用GRACE、GLDAS_2.1及GLDAS_1数据计算青藏高原区域水储量变化。结果表明，GLDAS_2.1水储量变化的年振幅、年相位与GRACE基本一致，而GLDAS_1在北部区域存在明显异常；GLDAS_2.1与GRACE年振幅、年相位具有较好的相关性，在中部、东北及东南部区域尤为明显，而GLDAS_1与GRACE在北部区域相关性不大；GLDAS_2.1比GLDAS_1能更好地反映青藏高原区域质量的季节性变化。
| [1] |
Tapponnier P, Xu Z Q, Roger F, et al. Oblique Stepwise Rise and Growth of the Tibet Plateau[J]. Science, 2001, 294(5547): 1671-1677 DOI:10.1126/science.105978
(  0) |
| [2] |
Matsuo K, Heki K. Time-Variable Ice Loss in Asian High Mountains from Satellite Gravimetry[J]. Earth and Planetary Science Letter, 2010, 290(1-2): 30-36 DOI:10.1016/j.epsl.2009.11.053
( 0) |
| [3] |
Yi S, Sun W K. Evaluation of Glacier Changes in High-Mountain Asia Based on 10 year GRACE RL05 Models[J]. Journal of Geophysical Research: Solid Earth, 2014, 119(3): 2504-2517 DOI:10.1002/2013JB010860
( 0) |
| [4] |
Xiang L W, Wang H S, Steffen H, et al. Groundwater Storage Changes in the Tibetan Plateau and Adjacent Areas Revealed from GRACE Satellite Gravity Data[J]. Earth and Planetary Science Letters, 2016, 449(1): 228-239
( 0) |
| [5] |
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-30229 DOI:10.1029/98JB02844
( 0) |
| [6] |
Tapley B D, Bettadpur S, Watkins M, et al. The Gravity Recovery and Climate Experiment: Mission Overview and Early Results[J]. Geophysical Research Letter, 2004, 31(9): 4
( 0) |
| [7] |
Jacob T, Wahr J, Pfeffer W T, et al. Recent Contributions of Glaciers and Ice Caps to Sea Level Rise[J]. Nature, 2012, 482(7386): 514-518 DOI:10.1038/nature10847
( 0) |
| [8] |
Yao T, Thompson L, Yang W, et al. Different Glacier Status with Atmospheric Circulations in Tibetan Plateau and Surroundings[J]. Nature Climate Change, 2012, 2(9): 663-667 DOI:10.1038/nclimate1580
( 0) |
| [9] |
Zhang G Q, Yao T D, Xie H J, et al. Increased Mass over the Tibetan Plateau: From Lakes or Glaciers?[J]. Geophysical Research Letters, 2013, 40(10): 2125-2130 DOI:10.1002/grl.50462
( 0) |
| [10] |
Liang S M, Gan W J, Shen C Z, et al. Three-Dimensional Velocity Field of Present-Day Crustal Motion of the Tibetan Plateau Derived from GPS Measurements[J]. Journal of Geophysical Research:Solid Earth, 2013, 118(10): 5722-5732 DOI:10.1002/2013JB010503
( 0) |
| [11] |
Yi S, Freymueller J T, Sun W K. How Fast is the Middle-Lower Crust Flowing in Eastern Tibet? A Constraint from Geodetic Observations[J]. Journal of Geophysical Research: Solid Earth, 2016, 121(9): 6903-6915 DOI:10.1002/2016JB013151
( 0) |
| [12] |
Zhou S Q, Kang S C, Chen F, et al. Water Balance Observations Reveal Significant Subsurface Water Seepage from Lake Nam Co, South-Central Tibetan Plateau[J]. Journal of Hydrology, 2013, 491(1): 89-99
( 0) |
| [13] |
Rodell M, Houser P R, Jambor U, 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) |
| [14] |
Zaitchik B F, Rodell M, Olivera F. Evaluation of the Global Land Data Assimilation System Using Global River Discharge Data and a Source-to-Sink Routing Scheme[J]. Water Resources Research, 2010, 46(6): 2840-2849
( 0) |
| [15] |
Spennemann P C, Rivera J A, Saulo A C, et al. A Comparison of GLDAS Soil Moisture Anomalies Against Standardized Precipitation Index and Multisatellite Estimations over South America[J]. Journal of Hydrometeorology, 2015, 16(1): 158-171 DOI:10.1175/JHM-D-13-0190.1
( 0) |
| [16] |
Bi H Y, Ma J W, Zheng W J, et al. Comparison of Soil Moisture in GLDAS Model Simulations and Insitu Observations over the Tibetan Plateau[J]. Journal of Geophysical Research:Atmospheres, 2016, 121(6): 2658-2678 DOI:10.1002/2015JD024131
( 0) |
| [17] |
Swenson S C. GRACE Monthly Land Water Mass Grids Netcdf Release 5. 0[EB/OL]. http://dx.doi.org/10.5067/TEMSC-OCL05, 2017-02-21
( 0) |
| [18] |
Landerer F W, Swenson S C. Accuracy of Scaled GRACE Terrestrial Water Storage Estimates[J]. Water Resources Research, 2012, 48(4): 4531
( 0) |
| [19] |
Swenson S C, Wahr J. Post-Processing Removal of Correlated Errors in GRACE Data[J]. Geophysical Research Letters, 2006, 33(8)
( 0) |
| [20] |
Cheng M K, Ries J C, Tapley B D. Variations of the Earth's Figure Axis from Satellite Saser Ranging and GRACE[J]. Journal of Geophysical Research:Solid Earth, 2011, 116(B1)
( 0) |
| [21] |
Swenson S C, 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) |
| [22] |
Geruo A. Computations of the Viscoelastic Response of a 3-D Compressible Earth to Surface Loading: An Application to Glacial Isostatic Adjustment in Antarctica and Canada[J]. Geophysical Journal International, 2013, 192(2): 557-572 DOI:10.1093/gji/ggs030
( 0) |
| [23] |
Rui H L, Beaudoing H K. Readme Document for Global Land Data Assimilation System Version 2 (GLDAS-2) Products[EB/OL]. http://dis.sic.gsfc.nasa.gov/services/grds-gds/gldas, 2014-01-10
( 0) |
| [24] |
Rodell M, Velicogna I, Famiglietti J S. Satellite-Based Estimates of Groundwater Depletion in India[J]. Nature, 2009, 460(7258): 999-1002 DOI:10.1038/nature08238
( 0) |
| [25] |
Tiwari V M, Wahr J, Swenson S C. Dwindling Groundwater Resources in Northern India, from Satellite Gravity Observations[J]. Geophysical Research Letters, 2009, 36(18): 184-201
( 0) |
| [26] |
Bolch T, Yao T, Kang S, et al. A Glacier Inventory for the Western Nyainqentanglha Range and the Nam Co Basin, Tibet, and Glacier Changes 1976-2009[J]. Cryosphere, 2010, 4(3): 419-433 DOI:10.5194/tc-4-419-2010
( 0) |
| [27] |
Yao T D. Glacial Fluctuations and Its Impacts on Lakes in the Southern Tibetan Plateau[J]. Chinese Science Bulletin, 2010, 55(20): 2071-2071 DOI:10.1007/s11434-010-4327-5
( 0) |
2. Faculty of Science, Hokkaido University, 5 Nishi, 8 Kita-Ku, Sapporo 0600808, Japan