2022, Vol. 42
2. 自然资源部国土卫星遥感应用中心，北京市百胜村1号, 100048
三江源是长江、黄河、澜沧江的发源地，是长江、黄河中下游地区和东南亚国家生态安全和区域可持续发展的生态屏障[1]。三江源地下水储量影响着区域生态环境变化和长江、黄河、澜沧江中下游的地面径流量，对其进行研究有助于了解全球变暖背景下青藏高原地区水储量的变化特征。
传统的地下水监测手段主要是通过地面监测井定期监测地下水的水位变化，难以实现大空间、大尺度下的地下水变化监测[2]。自2002年GRACE卫星发射以来，利用时变重力场模型研究大范围区域地下水的时空变化成为可能。在2017年GRACE卫星停止服务后，美国国家航空航天局(NASA)与德国地球科学研究中心联合研制了GRACE-FO卫星，继续承担全球及区域水资源变化监测的任务。
目前，已有大量的研究人员利用时变重力场数据研究区域水资源的变化[3-6]，但对三江源区域地下水的分析研究较少[7]。本文利用GRACE和GRACE-FO提供的2003~2020年的时变重力场模型、GLDAS模型和湖泊水量变化数据反演三江源区域地下水的时空变化，通过奇异谱分析方法(singular spectrum analysis, SSA)提取时间序列中的趋势信号和周期信号，分别对黄河源、长江源、澜沧江源不同区域的地下水变化进行分析，并结合区域内降水数据分析三江源地下水年度补给规律。
1 研究区概况三江源区域位于青海省南部，地理位置为31°19′~36°12′N、89°45′~102°23′E，总面积为302 500 km2，为长江、黄河和澜沧江的源头汇水区(图 1)。三江源区域是青藏高原的腹地和主体，以山地地貌为主，海拔1 971~6 565 m，河流密布，湖泊、沼泽众多，是世界上海拔最高、面积最大、湿地类型最丰富的地区，也是中国淡水资源的重要补给地。
|
图 1 三江源地区位置 Fig. 1 The location of the Three-River headwater region |
采用美国得克萨斯大学空间研究中心(CSR)提供的GRACE和GRACE-FO卫星RL06数据，模型最大阶数为96，时间范围为2003-01~2020-12。GRACE与GEACE-FO之间11个月(2017-07~2018-05)的数据缺失不作处理。此外，由于GRACE卫星在前期调试运行阶段及后期出现传感器性能下降和供能不足的问题造成20个月的模型缺失，对此，采用三次样条法进行插值[8]。
GLDAS数据选取GLDAS-NOAH水文模型，时间分辨率为1个月，空间分辨率为1°× 1°，时间范围与时变重力场数据时间范围一致。本研究主要提取水文模型中的4层土壤水(0~200 cm)、雪水当量和冠层水变化。
湖泊数据采用国家青藏高原科学数据中心的青藏高原2000~2017年高时间分辨率湖泊水位及水量变化数据集[9]，以及Cooley等[10]提供的2018~2020年全球湖泊水量变化数据集。从数据集中提取三江源地区13个大、中型湖泊的水量数据，湖泊位置如图 1所示。
降雨和温度数据来自国家气象科学数据中心提供的中国地面气温月值0.5°×0.5°格点数据集(V2.0)和中国地面降水月值0.5°×0.5°格点数据集(V2.0)，时间范围为2003-01~2020-12。
2.2 理论方法 2.2.1 地下水储量变化反演由于GRACE和GRACE-FO卫星对重力场低阶项不敏感，且模型中含有高阶噪声、奇偶阶相关误差等[11]，在计算之前需要对时变重力场数据进行预处理，步骤如下：1)用SLR(satellite laser ranging)得到的C20、C30项代替模型中的C20、C30项[12]；2)采用半径为300 km的高斯滤波对模型系数进行处理，抑制高阶噪声[11]；3)用滑动窗进行相关滤波，消除模型系数间的相关误差[13]。得到滤波后的陆地水反演公式为：
| $ \begin{array}{*{20}{c}} {\Delta H\left( {\theta , \lambda } \right) = \frac{{a{\rho _{{\rm{ave}}}}}}{{3{\rho _{\rm{w}}}}}\sum\limits_{l = 0}^N {\sum\limits_{m = 0}^l {{{\bar P}_{lm}}\left( {\cos \theta } \right)W\left( l \right)} } \frac{{2l + 1}}{{1 + {k_l}}} \times }\\ {\left( {\Delta {C_{lm}}\cos \left( {m\lambda } \right) + \Delta {S_{lm}}\sin \left( {m\lambda } \right)} \right)} \end{array} $ | (1) |
式中，ΔH为反演的陆地水储量变化，表示方式为等效水高值；λ为地心经度；θ为地心余纬；ρave和ρw分别为地球平均密度(5 517 kg/m3)与水的密度(1 000 kg/m3)；Plm为完全规格化缔合勒让德函数；W(l)为高斯平滑函数；kl为l阶负荷勒夫数；ΔClm为月时变重力场Clm系数相对于重力场Clm系数平均值的变化量；ΔSlm为月时变重力场Slm系数相对于重力场Slm系数平均值的变化量; N为时变重力场模型的最大阶数。
经式(1)处理后的等效水高值仍受到冰后回弹和误差泄露的影响，可分别利用ICE-6G-C模型和正演模型消除冰后回弹和误差泄露的影响[14]。
依据水量平衡公式，三江源地区的地下水储量变化可表示为：
| $ \begin{array}{*{20}{c}} {{\rm{\Delta GWS}} = {\rm{\Delta TWS}} - {\rm{\Delta SMS}} - }\\ {{\rm{\Delta SWE}} - {\rm{\Delta CWS}} - {\rm{\Delta LAKE}}} \end{array} $ | (2) |
式中，ΔTWS为陆地水储量变化，由时变重力场模型计算得到；ΔSMS为土壤湿度含水量变化，ΔSWE为雪水当量变化，ΔCWS为植被冠层含水量变化，这3个变量都可以从GLDAS的NOAH模型中获得；ΔLAKE为湖泊水量变化，由湖泊水量数据集计算得到。
将13个湖泊分别看作是离散的点质量，收集到范围内13个湖泊的经纬度位置和水量变化速率，得到三江源内13个等效点质量及其位置。根据式(3)得到球谐展开系数，并利用式(4)得到空间域的湖泊水量变化速率[15]：
| $ \left\{ {\begin{array}{*{20}{c}} {\Delta {C_{lm}} = \sum\limits_{i = 1}^n {\frac{{\Delta {M_i}}}{{4{\rm{ \mathsf{ π} }}{a^2}{\rho _{\rm{w}}}}}{{\bar P}_{lm}}\left( {\cos {\theta _i}} \right)\cos m{\lambda _i}} }\\ {\Delta {S_{lm}} = \sum\limits_{i = 1}^n {\frac{{\Delta {M_i}}}{{4{\rm{ \mathsf{ π} }}{a^2}{\rho _{\rm{w}}}}}{{\bar P}_{lm}}\left( {\cos {\theta _i}} \right)\sin m{\lambda _i}} } \end{array}} \right. $ | (3) |
| $ \begin{array}{c} \Delta L\left( {\theta , \lambda } \right) = \sum\limits_{l = 0}^N {\sum\limits_{m = 0}^l {\left( {\Delta {C_{lm}}\cos m\lambda + } \right.} } \\ \Delta {S_{lm}}\sin m\left. \lambda \right){{\bar P}_{lm}}\left( {\cos \theta } \right) \end{array} $ | (4) |
式中，n为点质量总数，n=13；ΔMi为第i个点质量变化速率；θi、λi分别为第i个点的余纬和经度；N为球谐系数截断阶数。
2.2.2 奇异谱分析奇异谱分析(SSA)是一种处理非线性时间序列的方法，通过对所研究时间序列的轨迹矩阵进行分解、重构等操作，提取出时间序列中趋势分量、周期分量、噪声分量等[16]。SSA的简要处理过程为：构建轨迹矩阵、进行奇异值分解(SVD)、分组、对角平均化，具体公式可参考文献[17]。
3 结果与分析 3.1 整体时变分析2003~2020年整个三江源区域地下水储量的平均变化时间序列如图 2所示。考虑到时间序列的连续性，在进行SSA分析时选择2003-01~2017-06的数据，窗口长度选为48个月(图 3)。
|
图 2 2003~2020年三江源区域地下水平均变化时间序列 Fig. 2 The time series of the average change of the groundwater in the Three-River headwater region from 2003 to 2020 |
|
图 3 SSA结果 Fig. 3 The results of SSA |
对前15阶重构成分(reconstruction components, RC)进行w-correlation分析，结果如图 3 (a)所示。从相关性分析结果中可以看出，第5个及之后的RC不能很好地相互分离，说明其中噪声信息占比很大；RC1为趋势项；RC3、RC4的相关性很好，归为同一信号，为周期项；RC2中含有部分趋势项和周期项，二者没有很好地分离。
RC1的趋势分量如图 3(b)所示，可以看出，2013-01之前三江源地下水储量有一个缓慢增加的趋势，2013-01之后有一个缓慢减少的趋势，故以2013-01为界对三江源地下水变化进行空间分析。
RC3和RC4组合的周期分量如图 3(c)所示，可以看出，地下水变化有一个明显的年周期。对周期分量时间序列的每年相同月份相加再取平均，得到地下水储量的年度变化规律如图 4所示。可以看出，在1 a中，三江源区域地下水变化基本能够维持动态平衡，随季节发生周期性波动：1~3月地下水量不断减少，3月地下水储量最少，3~9月地下水量缓慢增加，9月达到峰值，9~12月缓慢减少。
|
图 4 三江源地下水年度变化规律 Fig. 4 Annual variation of groundwater in the Three-River headwater region |
以SSA提取的趋势信息为参考，将2003~2020年划分为2个时间段。2003-01~2012-12三江源区域地下水储量变化率如图 5(a)所示，可以看出，该时间段内三江源区域地下水储量变化大部分为正增长，仅在澜沧江源头东南部分呈现缓慢负增长。地下水储量变化率空间分布整体上呈现东南向西北方向增加的趋势。长江源区域均表现为正增长，西北方向部分变化率为1.2 cm/a左右，东南方向部分变化率为0.7 cm/a左右。黄河源区域地下水增长主要在西南部分，变化率为0.6 cm/a左右，东北部分地下水变化稳定，无明显增加或减少。澜沧江源区域西北部分地下水有轻微增加，变化率为0.5 cm/a左右，东南部分表现为轻微减少，变化率为-0.4 cm/a左右。2013-01~2020-12三江源区域地下水储量变化率如图 5(b)所示。该时间段内三江源区域地下水变化明显的地区有3个，分别为长江源区域西北部、黄河源区域东部、澜沧江源区域。地下水在长江源区域西北部表现为明显的正增长，变化率为1.3 cm/a左右，在黄河源区东部和澜沧江源区域表现为明显的负增长，黄河源区域东部变化率为-1 cm/a左右，澜沧江源区域变化率为-1.1 cm/a左右。
|
图 5 三江源区域地下水储量变化率 Fig. 5 The change rate of the groundwater reserves of the Three-River headwater region |
三江源区域地下水变化的平均时间序列如图 6所示。对时间序列线性拟合后，得到区域内地下水变化的趋势，并对变化趋势进行显著性检验。结果表明，2003-01~2012-12三江源区域整体表现为显著正增长，增长率为0.7 cm/a(P < 0.01)；长江源区域整体表现为显著正增长，增长率为1.0 cm/a(P < 0.01)；黄河源区域整体表现为显著正增长，增长率为0.5 cm/a(P < 0.01)；澜沧江源区域整体稳定，无显著增加或减少。2013-01~2020-12三江源区域整体表现为负增长，增长率为-0.3 cm/a(P < 0.01)；长江源区域整体稳定，无显著增加或减少；黄河源区域整体表现为显著负增长，增长率为-0.6 cm/a(P < 0.01)；澜沧江源区域整体表现为显著负增长，增长率为-1.2 cm/a(P < 0.01)。
|
图 6 三江源区域地下水变化时间序列 Fig. 6 Time series of groundwater changes in the Three-River headwater region |
对中国地面降水月值0.5°×0.5°格点数据集(V2.0)进行提取，分别绘制2003~2012年和2013~2020年的降雨量变化速率图(图 7)。由图可见，2003~2012年降雨量在长江源、黄河源与长江源交接处呈明显正增长，而在黄河源区域和澜沧江源区域呈现较为轻微的正增长。整体来看，这一时间段内降雨量在各个区域均呈现正增长，与这一时间段内地下水储量变化较为一致。2013~2020年降雨量在黄河源区东部和澜沧江区东南部有较为明显的正增长，在长江源区域有轻微的正增长。总体来看，区域内降水在这一时间段内依旧表现为增加，与区域内地下水变化表现为相反的趋势。对区域内降雨量和地下水变化时间序列进行相关性分析，得到相关性仅为0.35，表明三江源地下水变化的长期趋势和降雨量变化的相关性不强。
|
图 7 三江源区域降雨量变化分布 Fig. 7 Distribution of rainfall changes in the Three-River headwater region |
Deng等[18]研究认为，地下水变化的长期趋势受降雨量与温度的共同影响。从中国地面气温月值0.5°×0.5°格点数据集(V2.0)中提取三江源区域的温度月值，发现2003~2020年温度呈缓慢增加的趋势，18 a内平均温度上升0.76 ℃。温度的升高使区域内蒸散发量增加，和降雨量一起影响地下水变化。从图 8可以看出，2013年之前三江源的降雨量多于2013年之后的降雨量，图中阴影区域为降雨量最少的年份，地下水在这些年份中都有明显的下降趋势，表明三江源地下水的短期变化主要受降雨量变化的影响。三江源地区2003~2012年地下水呈现正增长归因于降雨量多、温度低，2013~2020年地下水呈现负增长归因于降雨量少、温度高。
|
图 8 降雨量与地下水变化对比 Fig. 8 Comparison of rainfall and groundwater changes |
对降雨量进行每年相同月份取平均处理，如图 9所示。可以看出，降雨量和地下水一样具有明显的年度变化特征：1月最少，1~7月逐渐增加，7月达到最大值，7~12月逐渐减少。对2个时间序列进行时间滞后相关性分析，得出地下水变化滞后降雨量变化2个月的相关性为0.97。表明地下水年度变化与降雨量年度变化规律一致，地下水变化的年度特征主要受降雨量的影响，三江源区域降雨对地下水的补给滞后期为2个月。
|
图 9 降雨量与地下水年度变化规律 Fig. 9 Annual variation of rainfall and groundwater |
1) 2003~2012年三江源区域地下水储量整体呈现增加趋势，增长率为0.7 cm/a；2013~2020年三江源地下水储量呈现减少趋势，增长率为-0.3 cm/a。
2) 2003~2012年长江源地下水储量呈现增加趋势，增长率为1.0 cm/a；黄河源地下水储量呈现增加趋势，增长率为0.5 cm/a；澜沧江源地下水储量稳定，无明显增加或减少趋势。2013~2020年长江源地下水储量整体稳定；黄河源和澜沧江源地下水储量均呈现明显的减少趋势，增长率分别为-0.6 cm/a和-1.2 cm/a。
3) 三江源区域地下水储量具有明显的年度变化规律：1~3月不断减少，3月达到最低值，3~9月缓慢增加，9月达到峰值，9~12月缓慢减少。地下水年度变化规律与降雨量年度变化规律一致，地下水的年度变化规律主要受降雨补给的影响，三江源区域降雨对地下水的补给滞后期为2个月。
| [1] |
李盈昌, 杨为民, 刘真真. 三江源区土地利用和景观格局变化分析[J]. 西南林业大学学报, 2014, 34(1): 58-66 (Li Yingchang, Yang Weimin, Liu Zhenzhen. Analyses on Change of Land Use and Landscape Pattern in the Three-River Headwater Region[J]. Journal of Southwest Forestry University, 2014, 34(1): 58-66)
(  0) |
| [2] |
刘晨, 许才军, 刘洋, 等. 基于GRACE RL06数据探测三江源地区陆地水储量变化[J]. 大地测量与地球动力学, 2020, 40(10): 1092-1096 (Liu Chen, Xu Caijun, Liu Yang, et al. Terrestrial Water Storage Changes in Three-River Source Region as Detected Based on GRACE RL06 Data[J]. Journal of Geodesy and Geodynamics, 2020, 40(10): 1092-1096)
( 0) |
| [3] |
Hu X G, Chen J L, Zhou Y H, et al. Seasonal Water Storage Change of the Yangtze River Basin Detected by GRACE[J]. Science in China Series D: Earth Sciences, 2006, 49(5): 483-491 DOI:10.1007/s11430-006-0483-5
( 0) |
| [4] |
罗志才, 李琼, 钟波. 利用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) |
| [5] |
贾旭, 文汉江, 刘焕玲, 等. GRACE时变重力场与GPS速度场地表沉降研究[J]. 测绘科学, 2020, 45(3): 15-22 (Jia Xu, Wen Hanjiang, Liu Huanling, et al. Study on the GRACE Time-Varying Gravity Field and GPS Velocity Field Surface Subsidence[J]. Science of Surveying and Mapping, 2020, 45(3): 15-22)
( 0) |
| [6] |
周苗, 常晓涛, 朱广彬, 等. 利用时变重力数据反演地下水储量变化[J]. 测绘科学, 2021, 46(3): 67-72 (Zhou Miao, Chang Xiaotao, Zhu Guangbin, et al. Groundwater Reserves Changes from Time-Varying Gravity Data[J]. Science of Surveying and Mapping, 2021, 46(3): 67-72)
( 0) |
| [7] |
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: 228-239 DOI:10.1016/j.epsl.2016.06.002
( 0) |
| [8] |
陈智伟, 张兴福, 周波阳, 等. GRACE月时变模型缺失值的填补方法与精度分析[J]. 地球物理学进展, 2020, 35(3): 815-822 (Chen Zhiwei, Zhang Xingfu, Zhou Boyang, et al. Filling Method and Accuracy Analysis of Missing Values of Monthly GRACE Temporal Gravity Model[J]. Progress in Geophysics, 2020, 35(3): 815-822)
( 0) |
| [9] |
Li X D, Long D, Huang Q, et al. High-Temporal-Resolution Water Level and Storage Change Data Sets for Lakes on the Tibetan Plateau during 2000-2017 Using Multiple Altimetric Missions and Landsat-Derived Lake Shoreline Positions[J]. Earth System Science Data, 2019, 11(4): 1603-1627 DOI:10.5194/essd-11-1603-2019
( 0) |
| [10] |
Cooley S W, Ryan J C, Smith L C. Human Alteration of Global Surface Water Storage Variability[J]. Nature, 2021, 591(7848): 78-81 DOI:10.1038/s41586-021-03262-3
( 0) |
| [11] |
李琼. 地表物质迁移的时变重力场反演方法及其应用研究[D]. 武汉: 武汉大学, 2014 (Li Qiong. Earth's Surface Mass Transport Recovered from Temporal Gravity Field and Its Applications[D]. Wuhan: Wuhan University, 2014)
( 0) |
| [12] |
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) |
| [13] |
詹金刚, 王勇, 郝晓光. 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) |
| [14] |
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) |
| [15] |
相龙伟, 汪汉胜, 贾路路. GRACE监测青藏高原及邻区陆地水储量变化结果的可变性[J]. 大地测量与地球动力学, 2017, 37(3): 311-318 (Xiang Longwei, Wang Hansheng, Jia Lulu. The Variability of Terrestrial Water Storage Changes in the Tibetan Plateau and Adjacent Areas Retrieved by GRACE Data[J]. Journal of Geodesy and Geodynamics, 2017, 37(3): 311-318)
( 0) |
| [16] |
Elsner J B. Analysis of Time Series Structure: SSA and Related Techniques[J]. Journal of the American Statistical Association, 2002, 97(460): 1207-1208
( 0) |
| [17] |
王鹏程, 郭金运, 王方建, 等. 基于奇异谱分析的汶川地震震前地倾斜数据分析[J]. 地球物理学进展, 2020, 35(2): 467-474 (Wang Pengcheng, Guo Jinyun, Wang Fangjian, et al. Singular Spectrum Analysis on Ground Tilt Data before Wenchuan Earthquake[J]. Progress in Geophysics, 2020, 35(2): 467-474)
( 0) |
| [18] |
Deng H J, Pepin N C, Liu Q, et al. Understanding the Spatial Differences in Terrestrial Water Storage Variations in the Tibetan Plateau from 2002 to 2016[J]. Climatic Change, 2018, 151: 379-393 DOI:10.1007/s10584-018-2325-9
( 0) |
2. Land Satellite Remote Sensing Application Center, MNR, 1 Baisheng Village, Beijing 100048, China