基于GRACE 的黄河流域陆地水储量时空变化研究

杨 飞，张成业，李 军，宋子恒，郭添玉

(1.中国矿业大学(北京) 煤炭资源与安全开采国家重点实验室，北京 100083；2.中国矿业大学(北京) 地球科学与测绘工程学院，北京 100083；3.煤炭开采水资源保护与利用国家重点实验室，北京 102211)

黄河流域是我国重要的经济地带，构成了我国重要的生态屏障。推进黄河流域水资源节约集约利用是黄河流域生态保护和高质量发展的主要目标任务之一。彭苏萍等[1]指出，水资源是支撑黄河流域煤炭开发、区域经济发展和生态环境协调的基础资源，而水资源短缺也是黄河流域中上游生态环境治理的瓶颈。因此，了解并掌握黄河流域水资源及其变化显得尤为重要，是当前黄河流域及其矿区生态环境研究的热点方向。传统监测手段大多利用地面水资源动态监测网点及相应的监测系统来实现，会受到监测站点分布和数目的影响。

GRACE(Gravity Recovery and Climate Experiment)重力卫星计划由美国国家宇航局和德国空间飞行中心联合实施，卫星于2002 年发射并于2004 年公开其重力数据[2]。GRACE 卫星采用极地近圆轨道低卫星跟踪技术，通过获取星间距离的变化及其残差来恢复时变重力场信息[3-4]。在消除天体引潮力对海洋、大气和相关地球动力过程引起的时变重力场影响后，GRACE 时变重力场反映的主要是两极冰盖、山岳冰川、陆地水储量以及海平面变化信息[5-6]。因此，利用GRACE 时变重力场可以反演陆地水储量变化信息，成为当前进行中长尺度陆地水资源时空变化分析的重要手段[7-9]。J.Wahr 等[10-11]推导了GRACE 数据计算陆地水储量变化的基础理论公式，首次利用GRACE卫星11 个月的数据反演密西西比河流域、亚马孙河流域和孟加拉湾地区的储水量变化，证明GRACE 具有在月尺度上探测到这些流域储水量变化的潜力。随后，众多学者提出多种基于GRACE 的陆地水储量反演的改进方法[12-17]，并成功应用于世界范围内不同区域陆地水储量的变化研究中[18-21]。

当前，利用GRACE 卫星数据计算陆地水储量变化的方法主要包括球谐系数方法和Mascons 方法[22-23]。邹贤才等[24]详细对比了上述2 种方法并发现Mascons 方法存在巨大优势；其他学者的研究结果同样证明Mascons 方法在一些流域的表现要优于球谐系数方法[25-26]。2015 年以来，美国得克萨斯州州立大学空间研究中心(CSR)和美国喷气推进实验室(JPL)发布用Mascons 方法计算陆地水储量的产品，分别称为CSRM 和JPL-M。2 种产品在不同区域陆地水储量变化的比较和研究中表明，CSR-M 产品可以更好地反映整个区域地表陆地水迁移的实际情况[27-28]。

基于此，利用CSR 提供的RL06 Mascons 水储量变化产品，开展黄河流域上、中、下游近20 a 陆地水储量变化的研究。通过纬圈长度加权平均计算每个子流域的水储量变化均值，采用时间序列分解方法得到黄河流域水储量变化的趋势、年周期等整体规律，利用空间分析比较不同区域水储量变化的特征，以期为更好地了解和掌握黄河流域水资源及其变化规律提供依据，为流域矿区的生态保护与可持续发展提供参考。

1 黄河流域概况

黄河发源于青海省巴颜喀拉山脉，横跨青藏高原、内蒙古高原、黄土高原和黄淮海平原，先后流经青海、四川、甘肃、宁夏、内蒙古、陕西、山西、河南和山东9 个省/自治区。其东西、南北跨度分别长约1 900 km和1 100 km，经纬度介于96°E-119°E 和32°N-42°N，流域总面积79.5 万km2。受大气环流和季风环流的影响作用比较复杂，黄河流域内不同地区的气候差异显著，主要属于南温带、中温带和高原气候区。根据黄河水利委员会的划分方案[29]，黄河上游指的是从源头到内蒙古河口镇，河长3 472 km，流域面积38.5 万km2；黄河中游指的是从内蒙古河口镇到河南郑州桃花峪，河长1 206 km，流域面积34.3 万km2；黄河下游指的是从郑州桃花峪到渤海，河长786 km，流域面积2.3 万km2。黄河流域中上游分布有晋北、晋东、晋中、黄陇、陕北、神东和宁东7 个国家大型煤炭基地，是我国目前主要煤炭生产区和调出区；而黄河下流分布有河南、鲁西等煤炭基地。黄河流域及其分界线如图1 所示。

图1 黄流流域分布Fig.1 Map of the Yellow River Basin

2 数据及分析方法

本文采用美国德克萨斯大学空间研究中心CSR提供的2002 年4 月到2017 年6 月的最新版本RL06 Mascons 产品。较RL05 版本，RL06 使用了改进的背景模型、新的地球系统参考框架和新的数据处理策略，显著提升了数据质量。CSR RL06 Mascons 产品是对GRACE 卫星的高精度星间观测数据Level1B 进行直接解算，C20 项采用卫星激光测月(SLR)数据的C20项替换，同时考虑地心一阶项和冰后回弹改正，采用Tikhonov 正则化进行方程组的求解。最终，CSR RL06 Mascons 数据采用NetCDF 格式，以0.25°×0.25°分辨率格网的形式发布。

为分析2002-2017 年黄河流域上中下游各个区域的水储量变化，需要计算每个区域的平均水储量变化值。采用对每个格网点的水储量变化数据直接取平均值的传统方法，相当于认为每个格网点等权。然而，CSR RL06 Mascons 产品为0.25°×0.25°格网点，数据在地球不同纬度之间代表的实际面积差异较大，需要考虑地球曲面的影响。本文采用纬圈长度加权方法，根据地理纬度为每个格网点赋予权重来获取不同区域的水储量变化均值。计算公式如下：

式中：Wr为某个区域的水储量变化均值；Δh为某一个格网点的水储量变化值；λi和φi分别为格网点经度和纬度；n为该区域总格网点数量。

本文采用时间序列分解方法来分析黄河流域上中下游不同区域水储量变化的时间序列特征。将各个区域水储量变化均值时间序列W(t)分解为趋势项、年周期项和半年周期项。利用最小二乘法进行参数拟合，拟合形式如下：

式中：t为时间；a为常数项；b为趋势项；c1、c2和d1、d2分别为年周期变化和半年周期变化。

3 结果分析

通过纬圈长度加权方法得到每个月黄河流域上中下游的水储量变化均值，采用等效水高表示[3-4,6-7]，其时间序列及线性拟合趋势如图2 左侧所示，可以看到，2002-2017 年，黄河上游水储量变化均值的波动范围较小，最大和最小值分别为32.8、-56.9 mm，线性拟合趋势项为-2.3 mm/a，水储量变化递减的趋势不明显。而黄河流域中游和下游的水储量变化均值表现出明显的递减趋势，其趋势项分别达到-10.3、-17.5 mm/a。在黄河流域中游，水储量变化均值的最大和最小值分别为108.9、-168.6 mm；在黄河流域下游，它们分别达到了143.8、-301.9 mm。上述结果表明近20 a 黄河流域3 个流域段水储量变化的差异性，上游水储量变化稳定，月均值差异较小，趋势性弱；中游的水储量变化呈递减趋势；下游的水储量变化范围最大且递减趋势尤为明显。

去趋势后，黄河流域上中下游水储量变化均值的时间序列如图2 右侧所示，3 个流段都表现出周年性及季节性变化。值得注意的是，黄河中下游水储量变化在2003 年秋季有显著的提升。这是因为在2003 年上半年，黄河流域来水量较常年减少40%～70%，为实测记录以来最小值；进入汛期以后，黄河流域大部分地区降水量增加，较常年增多10%～70%，洪水时水库得以蓄水；所以在2003 年秋季，黄河流域水储量变化迅速从负值变化到正值。

图2 黄河流域上中下游水储量变化均值的时间序列及趋势和去趋势后的时间序列Fig.2 Time series and trend of the mean water storage anomalies(left),and the time series after detrending(right) in the upper,middle and lower Yellow River Basin

为进一步分析黄河流域不同流段水储量变化的月度规律以及数据的变化幅度，本文统计了2002 年4 月到2017 年6 月3 段流域各个月份的水储量变化均值并绘制箱形图，如图3 所示。可以发现，黄河上中下游的水储量变化月度范围有显著的区域差异，上游整体处于-40～40 mm，而中游和下游分别处于-100～100 mm和-200～200 mm。黄河上中下游每个月的数据都在一个相对稳定的范围内波动，春冬季节的波动明显大于夏秋季节，这与该流域“冬干春旱，夏秋多雨”的气候特点密切相关。另外，有少数波动表现为异常值，这些异常值集中于夏季，受到夏季降水和高温蒸发等因素影响。3 个流域段多年平均水储量变化在9 月份达到最大值，且在整个夏季和秋季维持较大值。

图3 黄河流域上中下游水储量变化月数据箱形图Fig.3 The boxplots of monthly mean data of the terrestrial water storage anomalies in the upper,middle and lower Yellow River Basin

通过对黄河流域上中下游水储量变化均值的时间序列进行分解，得到其趋势、年周期和半年周期等信息，见表1。可以看到，时间序列分解得到的趋势项与线性拟合的趋势差别很小。黄河上游的趋势项只有-2.17 mm/a，年振幅和半年振幅都较小，说明该区域水储量变化在近20 a 一直处于比较稳定的状态。相较于黄河中游，黄河下游的趋势项、年振幅及半年振幅的增幅分别达到73%、38%和98%。上述结果表明，黄河流域上中下游水储量变化存在不同的趋势项和周期性。确定系数的范围为0～1，表征最小二乘法参数拟合的可靠程度的统计指标。上中下游的确定系数均大于0.5，表明时间序列分解得到较好的模型拟合参数。其中，中游确定系数为0.70，参数拟合结果优于上游和下游。

表1 黄河流域上中下游水储量变化均值的时间序列分解统计Table 1 Statistics of the water storage anomalies in the upper,middle and lower Yellow River Basin through time series decomposition

为表现黄河流域水储量变化的空间分布情况，本文对该区域0.25°×0.25°格网点数据进行时间序列分解，获取每个格网点2002 年6 月到2017 年6 月的时间序列拟合参数。基于此，开展黄河流域水储量变化空间分析。图4 表示水储量变化趋势(图4a)和确定系数(图4b)的空间分布。可以看到，黄河流域水储量变化为增加趋势的地区主要集中在黄河上游源头附近区域，其年增加量在5 mm 左右。黄河下游及中游靠近黄河的区域表现出明显的递减趋势，年减少量基本达到15 mm 以上。其他一些区域趋势项大多处于-10 mm/a 以下。黄河流域水储量变化存在随经度变化的规律，由西向东表现出越来越明显的递减趋势。从确定系数的空间分布可以看到黄河流域拟合结果的南北差异，流域北侧拟合可靠度整体优于流域南侧。整个流域确定系数大于0.5 的格网点占64%，大于0.3的格网点达到85%。

图4 黄河流域水储量参数的空间分布Fig.4 Spatial distribution of water storage factors in the Yellow River Basin

图5 表示黄河流域水储量变化的年周期振幅和半年周期振幅的空间分布。可以看到，整个流域年周期振幅都明显大于其半年周期振幅。黄河上游源头附近区域、中游桃花峪及黄河下游区域年周期振幅较大，达到30 mm 及以上。在上述区域，高山融雪、降水量等季节性差异是导致其年周期振幅相对较大的原因。黄河流域中上游其他区域的年周期振幅较小，处于10～20 mm。黄河流域水储量变化的半年周期振幅表现出由西北向东南增大的趋势，在黄河上游小于10 mm，黄河中游处于10～20 mm，黄河下游达到了30 mm。这与不同区域所处的气候环境等因素有直接关系。

图5 黄河流域水储量不同变化周期的空间分布Fig.5 Spatial distribution of water storage anomalies in different period in the Yellow River Basin

4 结 论

a.黄河上游水储量变化的趋势性不明显，只有-2.3 mm/a；中下游表现出显著的递减趋势，分别达到-10.3 mm/a 和-17.5 mm/a。GRACE 水储量变化数据很好地反映出2003 年黄河流域来水量、降水量的特殊性。

b.GRACE 水储量变化数据表明黄河流域水储量变化在9 月达到最大值，而春冬季节波动较大，与该流域“冬干春旱，夏秋多雨”的气候特点密切相关。水储量变化的异常值集中于夏季，与夏季降水和高温蒸发等因素相关。

c.黄河上游源头附近区域的水储量变化呈增长趋势，增长率为5 mm/a 左右。黄河流域水储量变化存在随经度变化的趋势，由西向东表现出越来越明显的递减趋势，最大递减率超过20 mm/a。

d.黄河流域年周期振幅明显大于其半年周期振幅。受高山融雪、降水量等季节性差异的影响，年周期振幅在上游源头附近区域、中游桃花峪及下游区域达到较大值。半年周期振幅表现出由西北向东南增大的趋势，可能与区域的气候环境等因素相关。