我国的华北平原位于太行山麓以东、黄河以北，包括北京市、天津市、河北省的全部平原以及河南省和山东省的黄河以北平原，是我国重要的粮食生产区、政治和经济中心.自20世纪中叶以来，由于气候变化、工农业和城市的快速发展等因素，华北平原成为中国水资源压力最大的地区之一 (Xia et al., 2012).近年来，华北地区地下水开采量已占当年水资源利用量的70%以上 (Zheng et al., 2010).长期的地下水超采已导致华北平原地下水水位不断下降 (Liu et al., 2001).在农业种植密度较高的太行山山前平原区 (主要包括石家庄、邢台和保定等地)，近四十年来的浅层地下水水位下降达50 m，河北中东部平原区的深层地下水水位下降高达90 m (Foster et al., 2004).在华北中东部平原区，地下水水位的持续下降产生了大面积的地下水降落漏斗，并由此导致地面沉降、海水入侵、地下水污染和土壤盐碱化等一系列地质环境问题，严重制约了社会经济发展 (张兆吉等, 2009; 费宇红等, 2009).因此，开展对华北平原地下水储量变化的监测工作，对合理利用和有效管理区域水资源，具有十分重要的科学意义和社会价值.

近十余年来卫星重力技术的快速发展，特别是美德联合研制的GRACE (Gravity Recovery and Climate Experiment) 重力卫星的成功发射，使得我们可以获取空间中长尺度全球和区域时变重力场变化信息 (Tapley et al., 2004; Wahr et al., 1998; 孙文科, 2002; 宁津生, 2002; 许厚泽, 2001).从GRACE观测的陆地水变化中扣除地表水和土壤水等的贡献，即可得到地下水的时空分布和变化趋势信息.GRACE卫星目前已成为监测全球大中尺度地下水变化的唯一观测手段.国内外学者在利用GRACE卫星重力监测地下水变化方面开展了长期研究工作.早期，Rodell和Famiglietti (2002)通过对美国大平原地区地下水储量变化的模拟研究，验证了利用GRACE研究区域地下水变化的可行性.GRACE卫星发射以后，Strassberg等 (2007)和Rodell等 (2007)从GRACE实际观测的陆地水变化中扣除水文模式输出的土壤水变化，研究了美国大平原和伊利诺伊州地区的地下水储量季节性变化，其估计结果与地下水水井观测的结果符合得很好，从而证明了GRACE重力卫星在监测地下水变化方面的巨大应用价值.Yeh等 (2006)、Swenson等 (2008b)和Strassberg等 (2009)利用GRACE和实测土壤水资料估计了美国伊利诺伊州、俄克拉荷马州和大平原地区的地下水变化，其结果与实测地下水水位资料有很好的一致性.随着GRACE观测时间长度的增加，国内外学者开始进一步关注大中尺度及流域尺度地下水变化的长期趋势.近年来，GRACE在全球多个地区成功探测到由于农业灌溉导致的地下水超采和持续地下水储量减少，如在印度北部地区 (Chen et al., 2014; Long et al., 2016; Rodell et al., 2009; Tiwari et al., 2009)、美国加州中央山谷地区 (Famiglietti et al., 2011; Scanlon et al., 2012) 和中东地区 (Joodaki et al., 2014; Voss et al., 2013).

在我国华北地区，国内外学者也逐步开展了利用GRACE监测陆地水和地下水变化方面的研究工作.钟敏等 (2009)利用2003—2007年的GRACE资料研究了中国大尺度陆地水变化趋势，发现华北地区存在约2.4 cm·a^-1的陆地水减少信号，并推测这一信号是由华北地下水超采导致的.苏晓莉等 (2012)利用2002年至2010年的GRACE数据也发现华北的陆地水储量存在明显的减少趋势.Feng等 (2013)进一步结合GRACE观测数据、水文模式数据和地下水水井实测数据，得到2003年至2010年华北地区的地下水储量正以-8.3±1.1 km³·a^-1的速率减少.Huang等 (2015)则进一步探讨了华北的山前平原和中东部平原地下水亏损速率的差异性.

以上研究工作主要关注于华北平原地下水变化的趋势项，目前对该地区地下水储量年际变化规律的研究以及地下水模型与GRACE结果的对比分析则相对较少.为此，本文在前人的研究基础上，提出基于空间约束的GRACE质量反演方法监测华北平原地下水储量变化，并通过与地下水实测资料和地下水模型结果的对比分析，研究了华北平原地下水年际变化特征与长期变化趋势，并探讨了深层地下水变化的贡献.

本文采用美国德克萨斯大学空间研究中心CSR (Center for Space Research，University of Texas at Austin) 提供的2002年4月至2014年11月共140个月的GRACE level-2(Release 05) 月重力场模型.数据为60阶球谐系数，并已扣除了非潮汐大气、非潮汐高频海洋信号、海潮、固体潮和固体极潮等的影响 (Bettadpur, 2007).受卫星的近极地轨道设计影响，GRACE无法观测到地心运动，同时对重力场的C₂₀项也不敏感.因此本文加回了Swenson等 (2008a)基于大气海洋模型与高阶GRACE数据计算的地心改正项；利用卫星激光测距 (SLR) 观测的C₂₀项替换了GRACE的原始C₂₀项 (Cheng and Tapley, 2004)；并采用400 km高斯平滑来降低高阶球谐系数的噪声 (Jekeli, 1981; Wahr et al., 1998).此外，由于地壳均衡调整 (GIA) 过程导致的重力场长期变化趋势，利用GIA模型进行了扣除 (A et al., 2013; Paulson et al., 2007).

在GRACE数据处理中，由于球谐展开的阶数限制以及滤波技术 (比如高斯平滑) 的使用，会引起信号“泄漏”误差 (leakage)：研究区域内的信号会泄漏到周围地区，并造成研究区域内的信号衰减 (leakage out)；同时，周围的信号也可能泄漏到研究区域内 (leakage in).为了获得研究区域内真实的质量变化时间序列，需充分考虑“泄漏”误差的影响.目前常用的“泄漏”误差改正方法主要有两种.一种方法是利用水文模式资料模拟的陆地水变化信息，通过进行与GRACE相同的后处理过程，来分别估计“leakage in”和“leakage out”的影响，并从GRACE得到的时间序列中扣除“leakage in”效应，并加回信号对外泄漏的“leakage out”效应.但是，该方法严重依赖于所使用的水文模式资料.由于目前的水文模式资料存在较大的不确定性，这导致对“泄漏”误差的估计也可能存在较大误差 (Klees et al., 2007; Longuevergne et al., 2010).另一种方法，首先也是利用水文模式资料估计“leakage in”影响，然后通过比较使用GRACE后处理前后研究区域内的水文信号的衰减比例，估计一个尺度因子 (scaling factor)，该尺度因子代表了研究区域内信号由于GRACE后处理导致的信号衰减的幅度.为获得真实的质量变化时间序列，从GRACE后处理得到的时间序列中扣除“leakage in”影响后，采用尺度因子来恢复研究区域内的真实信号强度.如果进一步假设研究区域内的信号是均一分布的，则可以不依赖于水文模式来估计尺度因子.这种基于均一假设的尺度因子估计方法，被广泛用于研究极地冰盖质量平衡、区域海平面变化和地下水变化等 (Feng et al., 2014; Rodell et al., 2009; Swenson and Wahr, 2007; Velicogna and Wahr, 2006).但是，如果研究区域较大，且区域内存在较为明显的质量分布差异，则基于均一假设估计的尺度因子也可能存在较大不确定性.

鉴于以上两种方法在估计GRACE时间序列方面的不足，近年来发展了基于空间约束的区域质量变化估计方法.Jacob等 (2012)将全球主要山岳冰川与极地冰盖的空间分布信息定义为球谐域的点质量，并与GRACE给出的球谐系数建立线性关系，从而估计了全球冰川与冰盖质量变化对海平面变化的影响.基于类似的思想，除了在球谐域构建点质量模型，在空间域也可以构建相应的点质量模型.基于空间域的约束方法已被用于研究高亚洲冰川和南极冰盖的质量变化 (Tang et al., 2012; Yi and Sun, 2014).该方法基于对信号源所在位置的先验信息，通过空间约束的方法恢复真实的质量变化信号，一定程度上避免了信号外泄问题.考虑到华北的地下水变化主要集中在平原地区，平原以外地区的地下水信号对华北平原内信号的“leakage in”效应相对较小.因此，本文需着重考虑华北平原的地下水信号向周边的信号外泄，即“leakage out”影响.以下将对本文利用GRACE估计华北平原地下水变化时所使用的空间约束方法进行简要阐述.

利用GRACE提供的大地水准面球谐系数估计地表某一点的密度异常的公式为 (Wahr et al., 1998)

(1)

其中，Δσ表示高斯滤波后得到的密度异常，a表示地球半径，ρ_ave表示地球平均密度，λ和θ表示经度和地心纬度，l和m表示球谐系数的阶和次，k_l表示l阶负荷勒夫数，W_l表示高斯平滑函数，P_lm(cosθ) 表示缔合勒让德函数，ΔC_lm和ΔS_lm表示扣除平均后的GRACE给出的大地水准面球谐系数.

而大地水准面的球谐系数可以表示为地表“真实”密度异常的球面积分（Chao and Gross, 1987）:

(2)

对上式做离散化处理后可得

(3)

其中，Δλ_i和Δθ_i表示离散化的格网划分间隔，i代表划分的地表质量单元.

将式 (3) 代入式 (1)，经整理后可得到

(4)

其中，下标i表示待求信号区域的格网点，共计I个点；下标j表示整个研究区域的格网点，共计J个点.记Δσ和Δσ分别为I维向量和J维向量，并记A为J×I的矩阵，且

(5)

则式 (4) 可简写为

(6)

上式建立了GRACE数据后处理得到的地表质量异常与待求的“真实”质量异常之间的函数模型.本文将华北地区按照0.25°的规则格网划分，取华北平原作为待求的信号区域，取经纬度范围 (34°N—43°N，110°E—120°E) 为整个研究区域 (即为华北平原区信号可能的泄漏范围).需特别说明的是，采用0.25°格网划分是为了更好地覆盖华北平原地区，并与实测地下水结果保持一致的格网间隔.GRACE实际的空间分辨率约为200~300 km，本文的空间约束方法则是基于先验信息将地下水信号约束在华北平原地区.由于式 (6) 中设计矩阵呈明显的病态性，导致无法使用经典最小二乘方法进行求解.为此，本文采取了Tikhonov正则化的方法来抑制设计矩阵的病态性，同时结合L曲线确定正则化参数.基于该空间约束方法，本文将“真实”信号 (即需研究的地下水变化信号) 约束在华北平原地区，进而减少高斯平滑导致的信号泄漏.

(7)

其中，κ为正则化参数，I为单位阵.

本文使用的水文模式数据有：美国宇航局哥达航空中心建立的四个版本的GLDAS水文模式 (NOAH、VIC、CLM和MOSAIC)(Rodell et al., 2004)、美国国家海洋和大气管理局气象预报中心提供的CPC水文模式 (Fan and van den Dool, 2004)、德国法兰克福大学提供的WGHM2.2b水文模式 (Döll et al., 2014)、美国国家大气研究中心CLM4.5水文模式 (Oleson et al., 2013) 和中国科学院地理科学与资源研究所构建的中国地区VIC水文模式 (Zhang et al., 2014).本文取八个水文模式输出的土壤湿度与雪量数据的平均，作为最终的估计值.根据模式输出结果发现，华北地区雪量变化对长期陆地水变化的影响可忽略 (＜0.01 km³·a^-1)，而土壤水变化的趋势项为-0.19± 0.08 km³·a^-1.本文对水文模式输出的格网数据做球谐函数展开，并进行了与GRACE数据相同的平滑处理，进而从GRACE观测的陆地水变化中扣除了土壤水与积雪的影响.

为了从陆地水中分离出地下水，还需扣除地表水的贡献.华北地区的大部分河流均建有大量蓄拦大坝，因此表面水体对长期陆地水变化的贡献主要来自于水库蓄水变化 (Han, 2003).基于《海河流域水资源公报》的统计结果，2002—2014年华北地表水库蓄水的变化趋势为0.28 km³·a^-1.本文在利用GRACE估算华北地下水变化时已将其扣除.

本文收集了中国水利部地下水监测中心公布的华北平原浅层地下水水位数据和地下水埋深等值线图.为了将浅层地下水水位变化转化为与GRACE观测等效的地下水储量变化，需将水位数据乘以给水度得到储量变化的等效水柱高.其中，给水度代表了地下水下降一个单位，单位体积含水层所能释放出来的重力水的体积.在华北平原，浅层含水层的给水度从沿海地区 (＜0.05) 向西逐步增加，在太行山山前平原区的给水度最大达到0.2以上 (张兆吉等, 2009).本文参考相关文献，得到华北平原地区的平均给水度约为0.06(张兆吉等, 2009).根据地下水水位变化时间序列和给水度数据，本文计算得到华北平原浅层地下水储量变化结果.

为了与GRACE估计结果和地下水实测资料开展对比和分析，本文进一步收集了考虑地下水变化的全球水文模型结果和华北平原区域地下水模型结果.其中WGHM水文模型给出了全球每月一值的地下水储量变化分布 (Döll et al., 2014)；而PCR-GLOBWB水文模型则只给出了2000年全球地下水储量变化速率分布图 (Wada et al., 2010).此外，本文收集了分别基于MODFLOW和MIKE SHE地下水建模软件开发的华北平原地下水模型结果 (Cao et al., 2013; Qin et al., 2013).与全球水文模型相比，华北平原的区域地下水模型更多地考虑了华北地区的土壤类型、水文地质参数、实测水文气象资料和人类活动影响等因素，并利用实测水井资料进行了模型校正和评估.

如图 1a所示，基于GRACE的空间约束方法反演的华北平原地下水储量变化与地面水井实测的浅层地下水水位变化，均表现出一致的年际变化特征.通过与华北平原每月降雨量资料和年平均降雨量异常的对比可以发现，地下水的年际变化规律与降雨的年际变化特征十分吻合.2003—2004年，GRACE观测的地下水储量增加和水井实测的地下水水位抬升，与该时间段的平均降雨量增加相吻合.特别是在2003年，华北平原降雨量较常年平均高出104 mm，GRACE观测结果与水井实测水位均表现出快速回升.2005—2009年，由于降雨量相对减少，华北平原的地下水储量与地下水水位也持续降低.在随后的2010—2013年，降雨量相对多年平均表现出增加趋势，该时间段GRACE观测的地下水储量则未出现明显减少，而地下水水位则有所回升.在降雨量显著减少的2014年，GRACE与水位资料则表明了地下水的快速亏损.

图 1

Fig. 1

图 1 2002—2014年，(a) GRACE得到的华北平原地下水储量变化时间序列、华北平原平均地下水水位变化时间序列以及每月降雨量数据.其中，虚线表示原始时间序列；实线表示扣除季节项并进行三个月滑动平均后的时间序列；(b) 华北平原年平均降雨量异常.降雨数据来自中国气象科学数据共享服务网 (http://data.cma.cn/data/cdcindex.html[2016-10-31]) Fig. 1 (a) Groundwater storage variations (equivalent water height) in the NCP from GRACE (blue) and groundwater level changes in the NCP from phreatic well observations (red) for the period of 2002—2014. Monthly precipitation data in the NCP are also shown in panel (a). Dashed lines present original time series. Solid lines show the time series in which seasonal cycles were removed and a 3-month boxcar smoothing was applied. (b) Yearly precipitation anomalies in the NCP from 2002 to 2014. Precipitation data from the China Meteorological Data Sharing Service System (CMDSSS) were obtained from http://data.cma.cn/data/cdcindex.html[2016-10-31]

自2002年以来，华北平原地下水储量存在明显的长期亏损趋势.GRACE估计得到的2002—2014年华北平原地下水储量亏损速率为-7.4±0.9 km³·a^-1，这相当于每年等效水高5.6±0.6 cm的质量亏损.基于地面水井资料估计的浅层地下水水位下降速率为-14.4±1.0 cm·a^-1.考虑华北平原的平均给水度与总面积后，基于水井估计的亏损速率为-1.2 km³·a^-1.GRACE结果与水井观测结果之间的巨大差异，说明了华北平原的深层地下水亏损速率可能远大于浅层.深层地下水的贡献将在3.3节详细讨论.

除了地下水变化的时间序列，图 2进一步给出了基于GRACE得到的每年华北平原地下储量变化的空间分布，并已扣除了研究时间段的多年平均.如图所示，从地理分布上看，华北平原的太行山山前平原地区的地下水储量变化特征明显，并存在长期亏损和明显的年际变化特征.例如：与2012年相比，华北平原地下水储量在2013年略有回升，但在2014年则表现出明显的大范围亏损.这与图 1中的GRACE时间序列以及降雨的年际变化特征吻合.如图 3所示，基于地面水井实测结果表明，华北平原浅层地下水的变化主要集中在太行山的山前平原地区.这主要是由于该地区是浅层地下水的主要开采区，且存在高强度的农业灌溉 (Foster et al., 2004; Liu et al., 2001).对比图 2和图 3可以发现，受限于卫星观测资料的空间分辨率，GRACE尚无法提供如水井资料给出的高空间分辨率结果，但其总体时空变化趋势与水井资料有一定的吻合.

图 2

Fig. 2

图 2 2002—2014年基于GRACE得到的每年华北平原地下水储量变化的空间分布 Fig. 2 Spatial patterns of yearly groundwater storage variations in the NCP estimated from GRACE from 2002 to 2014

图 3

Fig. 3

图 3 2002—2014年基于地下水水位数据和给水度数据得到的每年华北平原浅层地下水储量变化的空间分布 Fig. 3 Spatial patterns of yearly groundwater storage variations in shallow aquifers of NCP estimated from phreatic well observations and specific yield data from 2002 to 2014

图 4进一步给出了2002—2014年，华北平原地下水储量长期变化趋势的GRACE监测结果与地面水井观测结果.与图 3中地下水储量年变化的空间分布相一致的是，地面水井结果表现出华北浅层地下水变化趋势的显著空间差异性 (图 4b).在山前平原地区 (如石家庄和邢台等地)，浅层地下水亏损速率较大，达到每年6 cm以上；但在中东部平原区，浅层地下水的亏损并不显著，部分地区还略有增加.基于GRACE监测的地下水储量变化趋势，则主要表现为整个平原区 (特别是山前平原和中部平原) 的地下水亏损，且在数值上高于水井实测结果.这一方面是受目前GRACE的空间分辨率的制约，无法给出较高空间分辨率的地下水变化信息；另一方面，地面水井资料仅能给出浅层地下水变化趋势，而GRACE给出的结果是总的地下水储量变化趋势.因此，两者的差异说明了华北平原存在深层地下水的亏损.

图 4

Fig. 4

图 4 2002—2014年，(a) 基于GRACE得到的华北平原地下水储量变化速率；(b) 基于地下水实测与给水度数据得到的华北平原浅层地下水储量变化速率 Fig. 4 (a) Trend map of GWS variations in the NCP from GRACE from 2002 to 2014; (b) Trend map of GWS variations in shallow aquifers of NCP from phreatic well observations and specific yield data from 2002 to 2014

除地下水实测数据外，本文进一步将GRACE结果与全球和区域地下水模型结果进行了综合比较.如图 5a所示，基于MODFLOW区域地下水模型估计的2002—2008年华北平原地下水储量变化趋势，表现出明显的空间差异性.如图 5b所示，基于WGHM全球水文模型估计的结果也表现出山前平原和中部平原地区的快速亏损，但其估计的数值是区域地下水模型结果的两倍左右 (表 1).虽然基于PCR-BLOBWB水文模型估计的华北平原地下水储量变化速率，与GRACE和区域地下水模型结果更为吻合 (表 1)，但其空间分布上存在较大问题 (图 5c).通过对WGHM和PCR-BLOBWB两个全球水文模型的分析可以发现，两者均存在各自的问题.WGHM虽然较好地模拟了华北平原地下水储量变化趋势的空间分布，但在数值上远远高估了变化趋势；PCR-BLOBWB虽然在数值上与GRACE和区域地下水结果吻合，但其无法很好地模拟地下水变化的空间分布.与上述两种考虑地下水变化的全球水文模型相比，区域地下水模型能够较好地模拟华北地区的地下水储量时空分布特征.而GRACE估计结果无论是在空间分布上还是对总体趋势项的估计方面，均能够与区域地下水模型有很好的吻合.这也说明了GRACE能够有效监测华北平原地下水的时空分布.

图 5

Fig. 5

图 5 (a) MODFLOW地下水模型估计的2002—2008年华北平原地下水储量变化趋势 (Cao et al., 2013)；(b) 基于WGHM全球水文模型估计的2002—2014年华北平原地下水储量变化趋势 (Döll et al., 2014)；(c) 基于PCR-GLOBWB全球水文模型估计的2000年华北平原地下水储量变化趋势 (Wada et al., 2010) Fig. 5 Trend map of GWS variations in the NCP based on (a) MODFLOW model from 2002 to 2008(Cao et al., 2013); (b) WGHM model from 2002 to 2014 (Döll et al., 2014), and (c) PCR-GLOBWB model for the year of 2000(Wada et al., 2010)

表 1(Table 1)

表 1 基于GRACE、地下水实测、全球和区域地下水模型估计的华北平原地下水储量变化趋势的比较 (单位：km3·a-1) Table 1 Comparison of GWS depletion rates in the NCP estimated from GRACE, well observations, global and regional groundwater models (unit: km3·a-1) 数据源 时间段 趋势项 (km3·a-1) GRACE 2002—2014 -7.4±0.9 GRACE 2002—2008 -5.3±2.2 地下水实测 2002—2014 -1.2 地下水实测 2002—2008 -1.8 MODFLOW地下水模型 (Cao et al., 2013) 2002—2008 -4.4 MIKE SHE地下水模型 (Qin et al., 2013) 2002—2008 -6.4 WGHM水文模型 (Döll et al., 2014) 2002—2008 -10.3 WGHM水文模型 (Döll et al., 2014) 2002—2014 -11.2 PCR-GLOBWB水文模型 (Wada et al., 2010) 2000 -4.9

表 1 基于GRACE、地下水实测、全球和区域地下水模型估计的华北平原地下水储量变化趋势的比较 (单位：km3·a-1) Table 1 Comparison of GWS depletion rates in the NCP estimated from GRACE, well observations, global and regional groundwater models (unit: km3·a-1)

表 1综合给出了GRACE、地下水实测、全球和区域地下水模型估计的华北平原地下水储量变化趋势.图 6则更直观地给出了不同研究时间段，实测与模型给出的华北平原地下水储量亏损速率的比较.结合表 1和图 6可以发现，2002—2008年，基于MIKE SHE和MODFLOW的两个地下水模型估计的华北平原地下水储量亏损速率分别为-6.4 km³·a^-1和-4.4 km³·a^-1.这与GRACE同时期的估计值-5.3±2.2 km³·a^-1较为吻合，同时与PCR-GLOBWB水文模型给出的2000年的结果也十分一致.而同期的WGHM全球水文模型则明显高估了地下水亏损速率.

图 6

Fig. 6

图 6 基于GRACE、地下水实测、全球和区域地下水模型估计的华北平原地下水储量变化趋势.直线在横坐标上的长度代表了研究时间段，纵坐标代表该时间段估计的地下水亏损速率.阴影面积代表GRACE估计结果的不确定性 Fig. 6 Compilation of GWS depletion rates in the NCP estimated from GRACE, well observations, global and regional groundwater models. Lengths of lines on the horizontal axis represent time periods, while values on the vertical axis represent the GWS depletion rates during the corresponding time period. Shaded area presents the uncertainties of GRACE estimates

2008年以后，虽然目前尚没有可供参考的区域地下水模型结果，但是考虑到2002—2008年GRACE与区域模型估计结果十分吻合，我们可以推断，2002—2014年的GRACE结果仍可以较好地估计华北平原地下水储量变化趋势.2002—2014年，WGHM水文模型仍明显高估了华北平原地下水的亏损速率 (表 1).

华北平原地下含水层自上而下一般分为四个含水层组：第一含水层组为浅层地下水 (即潜水)，其下的三个含水层组为深层地下水 (即承压水)(Foster et al., 2004; 张兆吉等, 2009).在山前平原地区，由于第一和第二含水层组的混合开采，两者具有较好的水力联系，一般统称为浅层地下水.如图 4b所示，从潜水水井资料得到的华北平原浅层地下水储量变化趋势图可以看出，华北平原的浅层地下水亏损主要发生在山前平原地区.虽然GRACE估计的总的地下水储量变化趋势大于潜水水井结果，但其在空间分布上也表明了山前平原存在较大的地下水亏损 (图 4a).

针对浅层地下水储量变化的估计，可以用潜水水位乘以给水度获得.但是深层地下水储量的变化则又分为两个部分.第一个部分是承压含水层的弹性释水，这可以利用承压水的水井资料乘以弹性贮水系数获得.由于华北地区弹性贮水系数约为0.00125，较给水度小约2个数量级，因此这部分弹性地下水储量变化可忽略不计 (Cao et al., 2013).第二个部分是非弹性的压密释水，这部分地下水亏损会导致含水层的压缩和地表的沉降.与山前平原不同，在华北平原的中东部地区，由于浅层地下水主要为无法使用的咸水，地下水的开采以承压水为主.中东部地区严重的深层地下水超采已经导致大面积的地面沉降，如在天津、沧州和衡水等地 (Liu et al., 2001; 费宇红等, 2009).本文基于GRACE估计的华北平原地下水亏损速率远大于潜水水井观测结果，这说明了华北地区存在较大的深层地下水亏损.

本文使用的两个华北平原地下水模型，同时模拟了浅层和深层地下水的变化，其估计的地下水储量变化趋势结果也远大于潜水水井实测结果 (表 1).2002—2008年，两个地下水模型估计的地下水亏损速率的平均值为-5.4 km³·a^-1，这与同时期本文利用GRACE估计的结果-5.3±2.2 km³·a^-1吻合.从GRACE中扣除同期潜水水井实测的浅层地下水的贡献后，可以得到：2002—2008年和2002—2014年两个时间段，华北平原深层地下水的亏损速率分别为-3.5±2.2 km³·a^-1和-6.2±0.9 km³·a^-1.

如图 7所示，受GRACE空间分辨率的制约以及水文模式的不确定性，GRACE估计的华北平原地下水变化趋势存在较为明显的信号“泄漏”问题，即部分华北平原地下水亏损信号位于平原区以外.在图 7所代表的整个研究区域 (34°N—43°N，110°E—120°E) 内，经300 km和400 km高斯平滑处理后得到的2002—2014年华北地下水亏损速率分别为-7.1±0.9 km³·a^-1和-5.7±0.8 km³·a^-1；如仅以华北平原作为研究区域，则经300 km和400 km高斯平滑处理后的结果分别为-1.7±0.3 km³·a^-1和-1.3±0.2 km³·a^-1.以上结果表明，随着高斯平滑半径的增大，华北平原地下水趋势信号的强度减弱；而仅以华北平原作为研究区域，则会产生更为严重的信号低估.

图 7

Fig. 7

图 7 2002—2014年基于GRACE得到的华北地区地下水储量变化速率 (未使用空间约束).(a)、(b) 分别表示使用300km高斯平滑和400 km高斯平滑的结果 Fig. 7 Trend maps of GWS variations in the North China from GRACE during 2002—2014, with (a) 300 km or (b) 400 km Gaussian filter applied. No spatial constraints are applied

如2.1节所述，使用高斯平滑导致的研究区域内信号强度的减弱，可以通过使用尺度因子或附加空间约束来恢复.对采用300 km和400 km高斯平滑的华北平原GRACE估计结果分别进行尺度因子恢复，得到2002—2014年地下水亏损速率分别为-6.2±1.1 km³·a^-1和-7.0±1.0 km³·a^-1；而对采用300 km和400 km高斯平滑的GRACE数据分别进行空间约束后的结果为-7.0±0.9 km³·a^-1和-7.4±0.9 km³·a^-1.通过与区域地下水模型结果的独立对比后发现，基于均一假设的尺度因子方法仍会不同程度地低估华北平原地下水的真实信号.本文所采用的基于空间约束的GRACE反演方法，则可以更好地揭示华北地下水变化的时空分布和长期趋势.

本文基于GRACE的空间约束方法分析了华北平原地下水储量变化的时空分布特征，并与水井实测资料、全球和区域地下水模型进行了综合比较和分析，研究结果表明：

(1) 华北平原地下水变化存在明显的年际变化特征，并与该地区的降雨异常相关.在降雨较多的2003—2004年，地下水储量显著回升；在降雨较少的2005—2009年，地下水储量显著减少；随着2010年以来降雨量的逐步增加，地下水储量有所回升；2014年，华北平原年降雨量显著减少，地下水的亏损则也相应较大.

(2) 除年际变化外，华北平原地下水存在明显的长期亏损趋势.GRACE估计得到的2002—2014年华北平原地下水储量亏损速率为-7.4±0.9 km³·a^-1.基于水井资料得到的华北平原浅层地下水的亏损速率为-1.2 km³·a^-1，远小于GRACE估计结果.两者之间的巨大差异，说明华北平原存在明显的深层地下水亏损.2002—2014年，结合GRACE与水井资料估算得到的华北平原深层地下水亏损速率为-6.2±0.9 km³·a^-1.

(3) 在空间分布上，地面水井实测资料表明，华北平原的浅层地下水亏损主要发生在太行山的山前平原地区.GRACE估计结果也表明，华北平原的地下水储量减少主要发生在山前平原和中部平原区，与水井资料较为吻合.基于WGHM全球水文模型估计的地下水亏损速率，在空间分布上也与GRACE和水井资料吻合.

(4) 2002—2008年，GRACE估计结果与基于MODFLOW和MIKE SHE地下水建模软件构建的区域地下水模型结果十分吻合.区域地下水模型结果提供了GRACE结果的独立验证，说明GRACE重力卫星可以有效监测华北平原的地下水储量变化趋势.而全球水文模型由于未充分考虑华北地区的实测资料和水文地质信息等因素，无法较好地模拟该地区的地下水变化.其中，WGHM模型高估了华北平原的平均地下水亏损速率，而PCR-BLOBWB模型则无法较好地模拟地下水变化的空间分布.

虽然本文基于空间约束的GRACE质量反演方法有效地估计了华北平原的地下水储量亏损速率，但是与地面水井资料和区域地下水模型结果相比，GRACE重力卫星仍无法提供较高分辨率的地下水储量空间变化信息.这在根本上受制于目前GRACE卫星的观测精度及其轨道设计，且无法通过数据后处理来完全解决.随着未来GRACE Follow-on卫星的发射以及下一代重力卫星计划的实施，有望提供更高时空分辨率的全球重力场模型，进一步拓展重力卫星的水文学应用范围.