0 引言

地下水是重要淡水资源.地下水储量变化常规监测方法包括压力计测量与地面监测网测量.压力计测量覆盖点到局部空间尺度，对实现大尺度地下水位变化监测存在困难.地面监测网测量根据研究区给水度，依靠实测水位数据可估计地下水储量在中长时空尺度上的动态变化特征(Rodell and Famiglietti, 2001)，但这种方式需在局部测量之间插值或非采样区域外推计算，并严格依赖样点测量的精度与可靠性，尤其在监测网动态维护与水位统测方面消耗大量人力物力(Tiwari et al., 2009; 罗志才等, 2012).除了上述两种观测手段之外，地下水文模型也可用于估计地下水储量变化，但在模型刻画及数据获取方面存在很大的局限性(Tiwari et al., 2009; Famiglietti et al., 2011).

重力反演与气候实验卫星(Gravity Recovery and Climate Experiment, GRACE)计划的出现，为监测地球表层质量变化提供了新手段(Tapley et al., 2004; Wahr et al., 2004).在较短时间尺度上，地球重力场主要敏感于表层大气、海洋及陆地水之间的质量交换，而陆地时变重力场主要由陆地水质量变化引起(Tapley et al., 2004).目前GRACE卫星重力技术在陆地水储量变化研究等方面得到了广泛应用(Rodell et al., 2001; Swenson and Wahr, 2002; Swenson et al., 2003; Wahr et al., 2004; 胡小工等, 2006; 罗志才等, 2012).

陆地水储量变化信号中混叠GRACE重力场模型的高阶项噪声(Wahr et al., 2006)，以及特定阶数之间的相关误差(Swenson and Wahr, 2006).球谐系数高阶次截断与滤波方法使用泄露了研究区水质量变化的真实信号.目前已有研究针对GRACE滤波后的结果进行尺度因子标定，以减弱信号的损失程度(Landerer and Swenson, 2012).Klees等(2006)提出了依靠水文模型原始格网时间序列的多时刻尺度因子法，该方法对水文模型的依赖程度较高；冯伟等(2012)则进一步研究了GRACE提取亚马逊流域陆地水变化的单一尺度因子法，该方法可对GRACE结果误差进行有效改正.

陆地水变化包括地下水，相应的时空变化可利用GRACE结合GLDAS、CPC、CLM4.0或CLM4.5等地表水文模型进行最新的定量估算(Tiwari et al., 2009; 冯伟等, 2017).对于地下水超采地区，如印度北部、加利福尼亚中谷、巴基斯坦及华北平原等地，GRACE技术可实现地下水储量的变化监测，监测结果可与常规监测手段及地下水文模型进行比对，结果表明两者一致性较好(Tiwari et al., 2009; Famiglietti et al., 2011; Iqbal et al., 2016; 冯伟等, 2017).

关中地区对地下水资源开发利用频繁，已引发局部多处地面沉降(陶虹等, 2013).

监测地下水动态变化对实现该地区水资源可持续利用具有重要意义.目前针对关中地区地下水的相关研究主要集中在地下水位变化规律及化学性质演变等方面(Yuan et al., 2011; Ma et al., 2012; 陶虹等, 2013; 贺旭波等, 2015; Duan et al., 2016).关于该地区地下水储量变化特征的研究则相对较少.

本文利用GRACE时变重力场及单一尺度因子法反演了关中地区地下水储量变化，与地下水实测资料及全球地下水文模型进行对比分析，并探讨降雨对地下水储量变化的影响.

1 数据与方法 1.1 GRACE数据后处理

本文采用美国德克萨斯大学空间研究中心(Center for Space Research, University of Texas at Austin, CSR)提供的GRACE level-2(RL05)月重力场模型数据，该数据为正则化球谐系数，扣除了固体潮、海潮、固体极潮、非潮汐大气和海洋影响，以及由于日月等其他行星引起的重力扰动(Bettadpur et al., 2007).GRACE时变重力场反演的陆地水储量变化以等效水高变化Δh_w(φ, λ)表示为(Jekeli et al., 1981; Wahr et al., 1998)：

(1)

式(1)中，θ为地心余纬，(φ, λ)为地面点的地心纬度和经度；(ΔC_nm, ΔS_nm)为n阶m次规格化位系数变化；P_nm(·)为经规格化的n阶m次缔合勒让德函数；k_n为n阶负荷勒夫数；ρ_w≈10³ kg·m^-3为水的密度；ρ_e≈5.5×10³ kg·m^-3为固体地球平均密度；R为地球平均半径.

直接利用公式(1)反演的水储量变化存在较明显的高频噪声，需截断高阶项球谐系数以降低高频误差影响，通常采用60阶重力场模型进行反演计算(Swenson and Wahr, 2006; Longuevergne et al., 2010; Feng et al., 2013; 吴云龙等, 2015; Belda et al., 2015).对接近60阶的高阶高次项位系数采用平滑半径为300 km的扇形滤波进行处理(Han et al., 2005; Zhang et al., 2009).奇偶阶位系数间各自序列的相关误差表现为“南-北”条带，利用P3M15反向边界延拓去相关滤波进行相关误差改正，同时调整固定滑动窗口宽度为5点(詹金刚等，2011).因单一滤波有其局限性，现将两种滤波方法进行组合使用，组合滤波计算模型如下：

(2)

式(2)中，逐阶、次的滤波核函数W_n、W_m通过递归计算得到，即

当l为阶数n或次数m时即为W_n和W_m，若W_m=1时，公式(2)演化成各向同性高斯滤波，若W_n=1, W_m=1，公式(2)退化成公式(1)；r为平滑半径.

本文采用的去相关滤波方法，基于最小二乘估计对l≥15阶的球谐系数拟合第l阶的相关误差，以第l阶为滑动窗中心向边界两端各反向延拓1/2窗宽的球谐系数，使得该方法在重力场模型系数的更高次数m上得到应用.建立法方程组求解三次多项式的未知参数a₀、a₁、a₂、a₃；b₀、b₁、b₂、b₃，构建相关误差模型为：

(3)

经过去相关误差改正后的球谐系数为：

1.2 单一尺度因子方法

在陆地水变化的时空尺度上，可认为研究区陆地水变化信号具有空间均一分布特点，且忽略区域外部信号变化影响(Chen et al., 2006).因此本文考虑基于区域整体及每个格网点的不同时段陆地水变化时间序列进行相同尺度的信号恢复.通过比较地表水文模型资料处理前后研究区的水文信号强度，估计一个信号振幅衰减比例因子即尺度因子，该尺度因子同样也代表了GRACE后处理所导致的信号衰减的幅度.这里采用单一尺度因子法对结果进行标定，使得经振幅恢复后的格点时间序列与真实时间序列尽可能接近.具体处理过程为：

(1) 提取全球水文模型GLDAS(NOAH)研究区格网数据，通过空间平均计算研究区陆地水变化时间序列ΔM_T.

(2) 同时将原始格网数据球谐展开至与GRACE相同阶次，按照GRACE后处理过程得到时间变化序列ΔS_T.

(3) 依据与原始时序的残差平方和最小，即

(4)

计算尺度因子scale，陆地水与地表水储量变化的振幅通过该尺度因子scale进行恢复.

1.3 水文数据

(1) GLDAS水文模型

GLDAS全球水文模型由美国宇航局哥达航空中心(NASA)和美国国家环境预报中心(NCEP)共同建立(Rodell et al., 2004)，反映陆地表面土壤水和冰雪变化，空间分辨率为1°×1°，每月一值，包括NOAH、VIC、CLM、MOSAIC四种陆地表面模式.

(2) WGHM水文模型

WGHM(WaterGAP Global Hydrology Model)全球地下水文模型(2.2c)由德国法拉克福大学自然地理研究所(IPG, University of Frankfurt)研制，该模型联合月时序参数以及气候和地理数据建立而成，时段为2003—2013年，空间分辨率为0.5°×0.5°，每月一值(Döll et al., 2014).

(3) 地下水测井资料

地下水测井资料来自陕西省地下水位年鉴，为关中地区监测井观测的地下水位埋深数据，观测时隔为5天，水位单位为m.为方便计算，本文直接引用年鉴的年平均水位数据.

2 结果与分析 2.1 GRACE与水文模型结果分析

利用2003年1月至2014年12月的月重力场模型对关中地区水储量变化进行反演.GRACE因其自身较高的轨道倾角，对重力场低阶项C₂₀ 敏感程度较低，可使用卫星激光测距SLR数据进行替换(Cheng and Tapley, 2005).扣除多年平均获得月重力场相对变化量，计算中缺失月份经插值得到.结合去相关滤波P3M15与300 km Fan滤波估计该地区0.5°×0.5°的陆地水储量变化.最后采用单一尺度因子法对数据结果进行标定，乘以按公式(4)计算得到的尺度因子1.588，作为研究区的最终估计.

陆地水储量变化ΔTWS实际上包括土壤水、地下水、地表水、冰雪及生物量的质量变化总和(Strassberg et al., 2009).根据陕西省水资源公报结果估算近年来关中地区主要河流年平均径流量长期变化趋势约为0.006 km³·a^-1(＜0.01 km³·a^-1).对长期陆地水变化影响来讲，地表径流变化可忽略不计.加之研究区的冰雪与生物量的变化相对较小，因此，本文可认为陆地水变化主要受土壤水与地下水变化影响.

(5)

式(5)中，Δ表示为水储量的相对变化，SM为土壤水含量，GWS为地下水储量.当ΔTWS和ΔSM分别依靠GRACE和陆地表面水文模型获得时，地下水储量变化ΔGWS可通过ΔGWS=ΔTWS-ΔSM进行估计.

同时将GLDAS格网数据球谐展开至与GRACE相同阶次，通过扣平均、组合滤波及尺度因子恢复法估算研究区地表水变化，为减小每种水文模型的不确定性，这里取GLDAS四种模型均值.对水储量变化格网经空间平均得到研究区整体变化结果.对比图 1a两种时间序列可知：季节性变化在GRACE与GLDAS估算的水储量时间变化中表现明显，两者整体变化趋势较为一致，但由于GLDAS模型未考虑深层土壤水及地下水影响，因此模型振幅小于GRACE.图 1b为基于最小二乘拟合计算研究区水储量的季节变化，振幅变化在±5.7 cm，GRACE估计的水储量最大亏损发生在每年的5—6月，最大盈余在每年的10月出现，与水文模型结果一致.

为了定量分析GRACE与GLDAS模型之间的相关程度，对两者相关系数进行计算，如表 1所示，相关系数多数高于0.7，GLDAS模型取平均后与GRACE相关系数明显增大，达到0.8，说明GLDAS平均结果与GRACE反演的陆地水储量变化更为吻合.

在利用尺度因子法获取区域时变重力场信号时，需对尺度因子改正前后GRACE及GLDAS结果的泄露误差进行估算.研究区水储量变化时间序列通过格网点空间平均得到，采用NOAH水文模型改正后结果均方根与滤波后的GRACE均方根RMS_GRACE及NOAH均方根RMS_model比值的乘积，如公式(6)，计算改正后的GRACE泄露误差E_{leak_scaled}^GRACE(Landerer and Swenson, 2012)：

(6)

式(6)中，ΔM_T、scale·ΔS_T在公式(4)中已定义.GRACE基于研究区的泄露误差为5.0 cm，经单一尺度因子改正后GRACE泄露误差减小为4.6 cm.同理GLDAS模型平均结果经改正前后的泄露误差分别为3.7 cm、3.4 cm.

2.2 地下水储量变化分析 2.2.1 与地下水模型对比分析

GRACE与GLDAS结果的相关性分析主要针对水储量随时间的变化，但评估地下水与模型结果之间的一致性程度应综合时空特征的比较.

直接从陆地水变化中扣除GLDAS均值获得研究区地下水储量变化等效水高，经线性回归得到地下水长期变化趋势的空间分布，如图 2a所示，同时将WGHM依照GRACE处理过程估计地下水储量变化.

2003—2013年地下水储量变化率呈现从区域西南方向到东北方向依次递减的阶梯分布，表现出地下水变化趋势的显著空间差异性，这与该地区降雨由南向北、由西至东逐渐递减的地理分布一致(冯晶和钱会, 2015).区域范围内地下水表现出不同程度的亏损，变化幅度在-1.2~-0.2 cm·a^-1，在东北区域达到最大，南部区域略有减少.但在WGHM估计结果中，西部及北部区域大幅亏损信号减弱，整个东南部区域与GRACE结果相比保持稳定，可以看出整个分布格局仍然从西南至东北依次减少，振幅变化在-0.7~-0.2 cm·a^-1，如图 2b所示.GRACE与WGHM估计结果在整体空间变化趋势上均表现为从西南至东北方向逐渐递减，但在局部长期变化速率上仍存在较大的差异性，这种差异可能由于WGHM在进行模型构建时采用的全球气候强度数据在局部空间上的不确定性导致模型本身存在一定误差，目前WGHM模型仍无法准确估计局部地下水变化的空间分布和亏损速率(Döll et al., 2014; 冯伟等, 2017).

除了空间分布总体上一致之外，研究区地下水储量变化整体均呈现长期减少趋势，如图 4所示，2003—2013年GRACE估计结果-0.37 cm·a^-1与WGHM结果-0.35 cm·a^-1十分吻合.2003—2008年期间，GRACE监测到研究区地下水整体逐年上升，这与陶虹等(2013)给出的同期水位埋深年值数据变化结果较为吻合，该年值数据经线性拟合后也表现为地下水的小幅回升.但在WGHM估计的地下水变化中难以体现，这可能由于WGHM地下水模型覆盖全球空间尺度，对区域整体部分时段的结果估算精度有待进一步改善.

2.2.2 地下水实测资料的对比分析

研究区地下水位监测井位置分布如图 3所示.观测数据经时空平均得到地下水位年变化，地下水位升降反映在地球重力场随时间的演变，采用2003—2008年地下水位数据对GRACE反演结果进行验证.已有相关研究表明，近10年来该地区地下水位持续下降局面因受控采限制及引水工程实施得到遏制.图 4a显示了GRACE观测2003—2008年的地下水储量相对多年平均有所增加，增加速率为0.25 cm·a^-1.地下水位在相同时间内也处于逐年上升状态，上升速率约为4.39 cm·a^-1，并在2004年表现出相比前一年的明显抬升，与GRACE观测结果比较一致，如图 4b所示.

若将GRACE观测的地下水储量变化与地下水位变化进行等效对比，需将地下水位结果乘以含水层给水度(冯伟等, 2017).因此根据两者数值结果估计关中地区给水度约为0.06.

2.2.3 降雨影响分析

降雨对关中地区地下水储量变化具有重要影响(贺旭波等, 2015)，本文进一步分析了两者之间的关系.采用中国气象数据网(http://data.cma.cn/data)提供的2003年1月—2014年12月的气象站月值降雨数据，得到关中地区整体降雨变化(分别以降水量实测值及距平值表示，单位mm).

图 5a为2003—2014年地下水储量变化及同期降雨实测值，对原始时序扣除季节项并进行三个月滑动平均得到地下水年际变化，研究区地下水等效水高长期变化趋势为-0.44 cm·a^-1，这与降雨长期时间分布较为一致，如图 5b所示.

图 5a显示了该地区每年7—10月存在大幅持续降雨，GRACE监测结果较降雨数据滞后，这主要因为降雨对地下水补给存在延迟所致.结合图 5b可以发现，地下水储量年变化特征与降雨的年变化规律比较吻合.特别是在2003年，关中地区降雨量较多年平均高出286 mm，GRACE观测的地下水表现出快速回升；2007—2008年地下水先小幅上升而后发生下降突变，与同期降雨先小幅增多，而后大幅减少相一致；2012—2013年地下水处于亏损状态，同期降雨呈逐年减少趋势；2014年降雨偏多，地下水储量衰减趋于平稳.

3 结论

本文利用GRACE卫星重力数据可监测到关中地区地下水储量变化，基于单一尺度因子法研究关中地区地下水储量的时空变化特征，并与全球地下水文模型及实测资料进行了综合比较与分析，研究结论如下：

(1) GRACE反演的关中地区陆地水储量变化与GLDAS模型平均结果在时间序列上具有很强的相关性.

(2) GRACE观测的地下水长期亏损速率存在明显的空间差异性特征，与WGHM结果在时空分布上吻合较好；GRACE结果与地下水实测水位在局部时段内均出现正增长.

(3) 地下水储量变化与年降雨异常密切相关，在年际变化及长期趋势上均表现出较好的一致性.

(4) 本文给出了关中地区给水度约为0.06，该结果可为水文地质学等相关研究提供一定参考.

致谢

感谢德国法拉克福大学自然地理研究所Döll P教授提供的WGHM 2.2c全球地下水文模型.