禤键豪 张兴福 陈鑑华 马资颖

1 广东工业大学土木交通工程学院，广州市外环西路100号，510006

GRACE卫星可为陆地水储量变化监测、海平面变化监测等提供重要数据支撑。CSR、JPL、GFZ、TU Graz等科研机构陆续发布了多个版本的GRACE月时变重力场模型，但由于不同机构发布的模型所采用的反演方法和数据处理策略不同，导致各模型精度略有差异，特别是在GRACE任务末期，不同模型间差异较大。Meyer等[1]采用方差分量估计方法融合AIUB-RL02、CSR-RL06、GFZ-RL06、GRGS-RL04和ITSG-Grace2018共5个月时变模型得到COST-G GRACE模型，理论上该模型噪声水平较小，可为利用GRACE时变模型反演区域陆地水储量变化提供更好的数据支持。

黄土高原位于我国中北部，是我国四大高原之一。随着黄土高原植被种植面积的增大，植被蒸腾作用会使该区域的陆地水储量减少[2]。此外，降水、气温、农业用水、工业用水等因素也会影响区域陆地水储量的变化[3]。本文综合利用COST-G GRACE月时变重力场模型、实测降水、实测气温、GLDAS(global land data assimilation system)模型、实测浅层地下水、原煤开采量、归一化植被指数(NDVI)等多源数据对黄土高原陆地水储量变化进行分析，最后利用偏最小二乘回归法分析该区域陆地水储量变化的驱动因素。

1 原理与方法

1.1 GRACE模型反演陆地水储量变化

利用时变重力场球谐模型计算区域质量变化对应的等效水柱高，公式为[4]：

(1)

式中，a为地球赤道处半径，ρave为地球平均密度，ρ为水密度，lm为归一化缔合勒让德多项式，kl是l阶勒夫数，ωlm为平滑核函数，本文采用DDK3滤波[5]，Δlm与Δlm为球谐系数变化量。

通常陆地水储量变化可近似表达为地表水储量变化与地下水储量变化之和[6]：

ΔTWS=ΔSW+ΔGW

(2)

式中，ΔTWS为陆地水储量变化，ΔSW为地表水储量变化，ΔGW为地下水储量变化。利用GRACE时变重力场模型时间序列可以估计ΔTWS，利用GLDAS模型可以计算ΔSW，再利用式(2)即可得到ΔGW。

1.2 实测浅层地下水变化

将实测浅层地下水井监测数据转换为水储量变化，计算公式为[7]：

(3)

式中，ΔGWwell为实测浅层地下水储量变化，Si为子区平均给水度，本文取黄土高原平均给水度为0.03[7]，Wi为子区面积，Δi为相应子区监测井的平均地下水位变化。

1.3 采煤耗水量

煤炭开采会影响区域陆地水储量变化，利用式(4)可将原煤开采量转换为采煤耗水量[3]：

TWScoal=μMc

(4)

式中，TWScoal为采煤用水量，μ为耗水系数，取为0.87 m3/t，Mc为原煤开采量。黄土高原内含陕西、山西、青海、宁夏、内蒙古、河南、甘肃等7个省或自治区，由于部分省或自治区并不完全位于黄土高原内，本文使用面积加权法获取该省或自治区位于黄土高原区域内的原煤开采量。

1.4 偏最小二乘回归

在偏最小二乘回归中，自变量对因变量影响的重要程度可由变量投影重要性(VIP)测定。一般来说，当VIP大于0.8时，自变量对因变量的变化具有重要解释意义，VIP计算公式为[3]：

(5)

式中，p为自变量个数，m为提取主成分个数，th为第h个主成分，Y为因变量，R(Y,th)为Y与th的相关系数，Whi为自变量在主成分th中的权重。本文取降水、气温、NDVI、采煤用水量、农业用水、工业用水、生活用水作为自变量，陆地水储量作为因变量计算VIP值。

2 陆地水储量变化时空分析

2.1 研究数据及处理流程

表1为用于研究黄土高原陆地水储量时空变化特征及其驱动因素的有关数据，共收集到74口浅层地下水监测井的水位数据，其中，45口井收录于《陕西省地下水监测年鉴》，其余收录于《中国地质环境监测地下水位年鉴》，其点位见图1。

表1 数据来源Tab.1 Data sources

黄土高原边界矢量数据来源于文献[8]图1 黄土高原浅层地下水监测井位置分布Fig.1 Distribution of shallow groundwater monitoring wells in the Loess plateau

将COST-G GRACE月时变重力场模型(下文简称为GRACE模型)截断至60阶次，采用图2中数据处理流程反演区域陆地水储量变化。

图2 GRACE模型数据处理流程Fig.2 Data processing flowchart of GRACE model

2.2 陆地水储量时间变化分析

图3为GRACE模型反演的陆地水储量变化与气温异常量(当月值减去平均值)、降水量对比。结果显示：1)气温、降水存在季节性变化特征，但陆地水储量季节性变化特征不明显；2)GRACE模型反演的陆地水储量变化在2003年末、2013年夏季、2016年秋季出现明显峰值，而在2003年全年山西、陕西降水较常年多20%以上[13]，2013年夏季降水较多，导致汾川河发生超历史洪水、渭河发生超警洪水[14]，2016-07中下旬黄土高原内多地降水较常年偏多，导致山西等地发生暴雨洪涝灾害[15]；3)GRACE模型反演的陆地水储量变化在2007年、2009年、2015年和2016年夏季出现较为显著的谷值，2007年夏季黄土高原部分地区由于降水较少使得旱情偏重且持续时间长[15]，2009年冬小麦主产区发生冬春连旱、西北部分地区发生伏旱和秋旱等[15]，2015年全年降水分布南多北少，特别是西北东部和黄淮等区域偏少10%～30%[14]，2016-01～05北方冬麦区发生冬旱，其中山西降水较常年同期偏少60%[15]。由此可见，黄土高原区域陆地水储量变化与气温、降水变化密切相关，GRACE模型反演的陆地水储量变化时间序列中较为明显的峰值与谷值往往与当时或之前的洪旱灾害相关。

图3 GRACE模型反演的黄土高原ΔTWS与降水、气温变化对比Fig.3 Comparison of ΔTWS derived by GRACE and precipitation and temperature in the Loess plateau

图4为黄土高原GRACE模型反演的陆地水储量变化、GLDAS模型计算的地表水储量变化、实测浅层地下水储量变化以及GRACE模型估计的地下水储量变化。由图4可知，4条曲线变化较为一致。2002-04～2016-12黄土高原陆地水储量变化与地表水变化、浅层地下水变化以及GRACE估计的地下水变化的相关系数分别为0.53、0.68、0.89，表明该区域陆地水储量变化与地下水变化相关性最高。2002-04～2016-12黄土高原陆地水储量变化大致可分为3个阶段：2002-04～2003-12以3.64±3.51 cm/a的速率上升(由于该时间段较短，下文不作定量分析)；2004-01～2009-12以1.64±0.25 cm/a的速率快速下降；2010-01～2016-12以0.57±0.29 cm/a的速率缓慢下降。在陆地水储量上升阶段，地表水、实测浅层地下水和GRACE估计的地下水均具有明显的上升趋势；在陆地水储量亏损的2个阶段，地表水分别以-0.57±0.18 cm/a和0.08±0.17 cm/a的速率变化，实测浅层地下水下降速率分别为0.35±0.08 cm/a和0.46±0.06 cm/a，GRACE估计的地下水下降速率分别为1.07±0.20 cm/a和0.65±0.24 cm/a。2002-04～2016-12 GRACE估计的地下水与实测浅层地下水的相关系数达到0.71，而RMSE为2.86 cm，两者具有较好的符合度，2种方法得到的地下水储量均在2004年年初达到峰值。结合图3和图4可知，地表水、降水也在2003年中后期出现极大值，而气温较低，由此推测地下水储量上升的主要原因为降水增加且蒸散量减少，经过一段时间的渗透作用使地下水出现滞后的上升信号，因此陆地水储量也出现相似的上升趋势。结合上述实测浅层地下水下降速率与GRACE估计的地下水亏损速率推断，2004-01～2009-12黄土高原深层地下水也可能存在明显的亏损。

图4 GRACE模型反演的黄土高原ΔTWS、GLDAS模型计算的ΔSW、实测浅层ΔGWwell和GRACE模型估计的ΔGW变化对比Fig.4 Comparison of ΔTWS derived by GRACE model, ΔSW calculated by GLDAS model, measured shallow ΔGWwell and ΔGW estimated by GRACE model in the Loess plateau

2.3 陆地水储量空间变化分析

图5为GRACE反演的陆地水储量、GRACE估计的地下水、GLDAS模型计算的地表水、NDVI、降水、气温等数据变化的空间趋势图。可以看出，GRACE模型反演的2002-04～2016-12陆地水储量在山西存在较大亏损，亏损中心位于黄土高原中东部，速率超过1 cm/a。GLDAS模型计算的地表水在研究时间内并无明显的空间变化趋势。GRACE探测到的2002-04～2016-12陆地水储量亏损信号与该区域地下水下降关联较大，且地下水亏损的空间位置与陆地水储量下降的位置相似。NDVI在黄土高原区域内具有较大的上升趋势，在上升最大的黄土高原中东部区域，降水也具有较大的上升趋势，但由于植被增多之后，地下水与地表水的补给会由于植被蒸腾耗水与冠层截留而降低，且该区域气温具有较大的上升趋势，结合地下水开采等因素，该区域陆地水储量显示为下降趋势。气温在研究区大部分区域均呈现上升趋势，这会增大蒸散作用，对陆地水储量亏损存在一定影响。

图5 黄土高原2002-04～2016-12陆地水储量、地下水、地表水、NDVI、降水、气温变化趋势Fig.5 Variation trend of TWS, GW, SW, NDVI, precipitation and temperature in the Loess plateau from April 2002 to December 2016

3 陆地水储量变化驱动因素分析

GRACE反演黄土高原陆地水储量的亏损大致分为2个阶段：2004-01～2009-12(时段1)呈现较大的亏损趋势，2010-01～2016-12(时段2)表现为较为缓慢的亏损趋势。因此，基于这2个阶段分析陆地水储量变化的驱动因素。其中，时段2原煤开采量数据大部分缺失，故暂不考虑该时段该因素的影响。

降水、气温、NDVI、采煤用水、农业用水、工业用水、生活用水等可直接或间接影响黄土高原陆地水储量变化，以上述7个因素作为自变量，以陆地水储量变化作为因变量，进行偏最小二乘回归分析，表2为各驱动因子的VIP值。在时段1中，采煤用水、生活用水、NDVI的VIP值大于0.8，均为陆地水储量下降的重要驱动因素，其中，采煤用水与生活用水的VIP值为1.60，这2个因素均与人类活动有直接关系。由于退耕还林工程的推进，人工林面积急剧增加[2]，NDVI在时段1中以每年0.003±0.002的速率上升，导致该区域被植被吸收的地表水和地下水增加。因此，在时段1，人类活动与植被作用对陆地水储量的下降起到推进作用。在时段2中，VIP值大于0.8的因素为工业用水、生活用水、气温，该时段中人类活动仍对陆地水储量的下降起到重要作用。但值得注意的是，气温上升也为黄土高原陆地水储量亏损的驱动因素：时段2中黄土高原地区的气温以0.16±0.14 ℃/a的速率上升，而NDVI以每年0.004±0.002的速率上升，黄土高原气温显著上升使得蒸发作用增大，加上植被增多会促进蒸腾作用，最终导致陆地水储量减少。

表2 偏最小二乘回归计算的驱动因子的VIP值Tab.2 VIP values of driving factors obtained by partial least squares regression

4 结 语

本文利用COST-G GRACE月时变重力场模型反演黄土高原陆地水储量变化，并与实测降水、实测气温、GLDAS模型计算的地表水变化、实测浅层地下水变化等进行对比，分析陆地水储量的时空变化特征，最后使用偏最小二乘回归分析陆地水储量变化的驱动因素，得到以下结论：

1)2002-04～2016-12黄土高原陆地水储量具有上升-下降-平缓下降的变化特征。其中，山西陆地水储量亏损最为严重，最大亏损速率超过1 cm/a。对比实测浅层地下水与GRACE估计的地下水可知，两者相关系数达到0.71，符合度较高，通过分析认为黄土高原陆地水储量变化与地下水变化密切相关。

2)经偏最小二乘回归分析可知，2004～2009年采煤用水和生活用水等人类活动以及植被作用是导致黄土高原陆地水储量亏损的重要驱动因素；2010～2016年工业用水、生活用水以及气温变化是黄土高原陆地水储量亏损的重要驱动因素。

致谢：感谢陈剑利教授、冯伟教授和ICGEM提供模型与软件。