华北平原位于中国东部，是中国的第二大平原，地处黄河下游，北依燕山，西傍太行山，东临渤海，包括北京、天津市与河北省的全部平原区，以及河南、山东省黄河以北的平原区，人口稠密，是中国重要的政治经济中心及工农业基地。

随着工农业和城市的快速发展及人口数量的不断增加，水资源的开采和利用规模及强度不断增大，水资源供需矛盾突出，华北地区地下水位持续下降，面临严重的水资源危机，人均水资源量很低，成为中国水资源压力最大的地区之一。水资源短缺成为制约经济社会和生态环境协调发展的关键因素。尤其在华北中东部平原区，地下水水位的长期开发导致了包括地面沉降、海水入侵、地下水污染和土壤盐碱化等一系列环境问题，严重影响社会经济发展和人类生活^[1-2]。所以，研究华北平原陆地水及地下水储量变化的空间分布特征和时间变化规律，对于该地区水资源可持续利用具有重要科学和社会意义。

华北平原的水储量变化问题受到科学家们的长期关注，其地下水变化成为近几年的热点研究问题。国内外诸多学者基于GRACE(gravity recovery and climate experiment)重力场模型在海平面变化、陆地水储量(terrestrial water storage, TWS)及地下水变化方面开展了长期的研究工作^[3-4]。GRACE重力卫星发射后，提供了前所未有的全球地球重力场模型，可以反映大尺度空间上的物质质量迁移信息，已被证实能十分有效地检测水储量变化，从而广泛应用于水文方面的研究^[5-8]。对于华北平原地区，钟敏等^[9]基于GRACE时变重力场模型研究中国的陆地水储量变化，发现在2003—2007年华北地区水储量以2.4cm/a(等效水高)的速率亏损。Feng等^[10]利用GRACE观测数据，结合水文模型和水井实测数据计算得到2003—2010年华北平原的地下水储量亏损速率约为2.2cm/a。Huang等^[11]考虑到地下水浅层与深层分布的不同，采用降尺度的方法将华北平原分为山前平原和中东部平原，进一步讨论两者地下水亏损速率的差异，结果显示两者在2003—2013年地下水变化速率分别为(-46.5±6.8)mm/a和(-16.9±1.9)mm/a，相对水井实测资料，GRACE更能充分地反映深层地下水信息。Gong等^[12]主要研究地下水储量的长期趋势和地面沉降的关系。Zhao等^[13]采用不同的GRACE产品分析华北地区2004—2016年以及2013—2016年地下水的亏损速率，结果表明2013年后随着降水的减少亏损速率加快。上述研究结果都表明华北地区水储量持续减少。但是，不同学者采用不同的观测资料和数据处理方法所得到的结果各有不同，并且他们更多的是关注地下水储量的变化。

为缓解华北地区水资源短缺的问题，中国实施了南水北调工程，将长江的水源向华北地区输送，是世界上最大规模的调水工程，输水线路如图 1所示。该工程包括西、中、东3条线路，目前东线、中线自2014年底已全面通水，使沿线缺水地区的供水保障能力有了大幅提升，北京、天津市及河北、河南的诸多城市，“南水”已经逐渐成为主要供水源，超过1亿人口受益。然而，南水北调工程是否缓解以及多大程度上缓解华北地区的水储量减少趋势，目前尚未看到较为科学的评估报告。

本研究目的是利用重力卫星GRACE与GRACE Follow-on观测数据及水文模型数据，分析华北平原水储量的时空变化特征，并量化南水北调工程对华北平原地区水储量的贡献。把2015年1月作为时间节点，对华北平原地区、北京市和天津市在2003—2015及2015—2019年期间的水储量变化分别进行趋势对比，讨论南水北调工程对这些地区水储量变化的影响。

1 研究区域及南水北调供水量

南水北调中线一期工程于2014年12月12日正式通水，以丹江口水库为起点，经陶岔渠首闸放水，分别流经河南、河北，经过保定市的西黑山分水口后，向北通往北京市，向东流入天津市，从而向沿途的20余座城市提供生产生活用水(图 1)。东线工程从长江下游江苏省扬州市水利枢纽中心抽引长江水，利用京杭大运河及与其平行的河道输水，并串联起调蓄作用的几个湖泊，出东平湖后分两路输水：一路向北，在位山附近经隧洞穿过黄河，自流到天津；另一路向东，通过胶东输水干线经济南输水到烟台、威海^[14]。东线一期工程主要是向江苏和山东两省供水。南水北调工程实施后，据人民日报报道，截至2019年12月12日，中线工程已通水5周年，累计调水258亿m³，且供水量在持续快速增加，在北京城区“南水”占自来水供应的73%。2018年末北京市平原区地下水位与上年同期相比回升1.94m，地下水储量相应增加9.9亿m³，地下水位上升区占37%，维持稳定区占45%^[15]；天津市平原区浅层地下水储量相比上年增加0.62亿m³，水位上升区面积占28.6%，相对稳定区占71%^[16]。南水北调中线总干渠向海河流域补水5.04亿m³，增加了河道水面面积，回补了地下水，改善了地下水超采的情况^[17]。

在本文的研究中，采用各省市水资源公报提供的每年调水量数据(表 1)，估算华北平原内各省份的受水区面积，算出受水区面积与全省(市)面积的比例，得到整个华北平原及北京、天津的南水北调水量。考虑到南水北调的水主要流入平原，本文的研究范围包含北京市(扣除山区)、天津市及河北省的平原地区，在河南和山东省的受水区面积计算中也减去了山区的面积。南水北调的水量本身相当于受水区水储量的变化，虽然有逐年略为增加的现象，但基本上处于较稳定的供水态势。因此，在下面计算比较中，我们把2015—2018年在华北平原供水量的平均值作为年供水量。2014年12月的输水量较少，暂且忽略2014年的全年变化，在处理GRACE观测数据时以2015年1月作为两个不同时期的分段点。南水北调水量的计算结果如表 2所示。

上述南水北调工程向华北平原地区提供了大量的水资源，该工程对华北水储量减少趋势是否有显著的缓解作用，输水后华北平原的水储量变化现状如何等问题是本文的研究内容，研究区域如图 1所示。

2 重力卫星GRACE与GRACE Follow-on观测数据及其处理

本研究采用重力卫星GRACE及其Follow-on观测数据对华北平原水储量进行计算分析。重力卫星GRACE是美国航空航天局和德国航空中心的合作项目，旨在获取高分辨率的时变重力场模型，发射于2002年，直到2017年10月27日才结束其生命周期。为实现与GRACE卫星的数据衔接，2018年GRACE Follow-on卫星发射。主要有3家机构提供GRACE数据，发布60阶和96阶的球谐系数产品。

本文采用2004年1月—2017年6月及2018年7月—2019年4月由美国德克萨斯大学空间研究中心(CSR)发布的GRACE及GRACE Follow-on卫星的RL06数据，其分辨率为300km。具体的数据处理方法是：利用卫星激光测距(SLR)观测的二阶项替换所有的二阶项球谐系数^[18]，利用Swenson等^[19]基于大气海洋模型与高阶GRACE数据提供的地心改正项增加一阶项系数，再通过扣除整段时间的平均值，根据Wahr等^[20]提出的理论得到陆地水储量TWS的变化。为去除噪声的影响，采用300km滤波半径的高斯平滑，结果用等效水高(equivalent water height, EWH)表示。对TWS时间序列做去周期处理，分别获取2004—2015年及2015—2019年的线性趋势，如图 2(a)所示。结果表明，华北平原水储量存在长期亏损趋势。

图 2还给出这两段时间内水储量变化趋势(斜率)的空间分布，可以看出太行山的山前平原地区水储量亏损特征明显，这与该地区以往发现的地下水储量变化特征相吻合^[21]。该地区水储量剧烈下降的主要原因是大规模开采地下水用于农业灌溉^[22-23]。2015年后亏损最大区域向山东西部方向转移。由于2015年后GRACE观测时间较短，对这段时间观测值进行趋势拟合存在较大误差。总体而言，南水北调后华北平原内陆地水储量亏损呈现缓解趋势。

3 基于尺度因子恢复的华北平原水储量变化

上述水储量变化数据(图 2)是在300km平滑滤波的条件下得到的结果，存在信号泄露问题，不够准确，需要恢复。在GRACE球谐系数展开过程中，受60阶截断的限制及高斯平滑的影响，会使研究区域内的信号和周边地区信号相互泄露、彼此干扰^[24-25]，造成研究区内的信号衰减，产生信号泄露误差。为此，用尺度因子的方法进行数据恢复。Landerer和Swenson^[26]分别检测在某个流域的单一尺度因子^[27]，以及流域内每个网格点的时间域尺度因子。前一种方法因单个因子而无法反映小区域或区域边缘处的实际质量变化，后一种方法是对水文模型在每个网格滤波与未滤波的质量变化时间序列进行拟合来确定尺度因子，相对而言恢复信号的精度更高，但严重依赖于所采用的水文模型。目前，在华北平原使用水文模型仍存在很大的不确定性，故本文采用一种新的空间域网格尺度因子方法。图 3展示空间域网格尺度因子方法的具体实现过程。首先假设研究区域内覆盖一层1cm的等效水高，为理想的原始信号，如图 3(a)所示。经过60阶截断和300km的高斯平滑后，信号衰减，数值明显变小。信号形态整体上呈近似圆盘扩散形状，在边界处的衰减、畸变尤为严重(图 3(b))。

考虑到对整个华北平原的水储量进行估计的同时，还需要对南水北调主要的受水区域(北京和天津市)进行估计。而受边界效应的影响，这两个地区的信号在GRACE处理后衰减、畸变较为严重。因此，采用理想核函数的思想，把研究区域内的信号权重处理为1，而区外权重处理为0，通过计算每个网格GRACE处理前后信号的衰减比率，得到所有网格的尺度因子，如图 3(c)所示，并应用到设定的研究区域上。

将上述网格尺度因子与上一节得到的GEACE每月网格数据对应相乘，并对其时间序列做与图 2相同的处理后，得到恢复后华北平原水储量变化的时间序列和两段时期内的趋势分布(图 4)。2004—2015年陆地等效水高变化趋势为(-4.49±0.40)cm/a，相当于水储量变化速率(-5.92±0.53)Gt/a。2015年后趋势向零值移动，变为(-2.25±2.12)cm/a，相当于水储量变化速率(-2.96±2.79)Gt/a。该结果显示，基于重力卫星观测得到的水储量变化趋势在南水北调前后的差值为2.96Gt/a，与南水北调的年平均输水量(3.08Gt/a，表 2)基本吻合。

结果表明，南水北调缓解了大约50%的华北水储量变化趋势，由原来的-5.92Gt/a减少为-2.96Gt/a。从空间分布形态看，信号恢复后趋势的相对分布特征(图 4(b)、4(c))与之前(图 2(b)、2(c))基本一致。从变化幅度的量级看，2004—2015年间在亏损最为严重的地区速率达到-10.94cm/a，2015—2019年亏损最大的地方只有-5.80cm/a，水储量亏损速率明显减弱。

为检验本研究结果的正确性，我们考察了前人所开展的研究和结果。Su等^[28]反演2002—2010年华北地区陆地水储量，其变化速率为-1.1cm/a。Moiwo等^[29]估计华北平原2002—2009年间的TWS变化，其亏损速率为-1.68cm/a。Feng等^[30]利用GRACE数据得到的2002—2014年间华北平原陆地水亏损结果表明，变化速率为-1~-2cm/a。相比这些结果而言，本文结果(-0.93cm/a)略小；但是考虑到研究所采用数据的长度不同及结果的离散度，这些结果总体上还是体现了较为一致的变化趋势。结果之间的差异主要来自以下原因：我们采用观测数据的时间段长度以及所划定的研究区域与前人有所不同。我们舍去了2002—2003年间的GRACE数据，因为这期间数据明显不稳定，对长年速率的估计影响较大；我们选取的观测数据时间段相对较长，便于研究南水北调前后的变化。此外，本文所采用的利用GRACE恢复陆地水储量的方法在天山及长江中下游区域的储水变化研究中已经得到验证^[31-32]。

4 讨论 4.1 气象因素对华北平原水质量平衡的影响

为进一步评估华北水储量变化的成因，本节讨论气象因素对华北平原水储量变化的影响。采用美国宇航局哥达航空中心建立的GLDAS水文模型^[33]，提取气象方面的数据。它包含4个版本的数据(NOAH、MOS、VIC、CLM)。采用GLDAS的子模型CLM给出的降水、蒸发、径流的逐月解计算华北平原水储量，以此确定气象因素对总水储量的影响。所用的水均衡公式为

$ \Delta {\rm{TWS}} = P - E - R, $

其中，P、R、E分别表示逐月的降水量、径流量与蒸发量。对于这几个气象因素分别计算其累积量，分段拟合趋势，计算这两段时间内它们的变化速率，结果总结在表 2中。两段时间内的TWS变化速率分别为-2.53Gt/a与-2.04Gt/a。

这说明气象因素对陆地水储量的贡献也为负值，且量级上具有重要的影响。2015—2019年的负变化趋势相对于2004—2015年有所减缓，在一定程度上也对华北平原水储量亏损起到了缓解作用。但是，由表 2可以看出，相对于南水北调的输水速率(3.08Gt/a)，由气象因素造成的水储量亏损减缓趋势(0.49Gt/a)只能解释2015年后总水储量速率变化的一小部分，不是主要因素。根据气象数据得到的ΔTWS与GRACE的结果存在较大差异，这是因为GLDAS水文模型的气象数据不太精确，尤其是径流量存在很大的不确定性，这或许可以解释表 2前两列的数据。此外，它们在前后两段时间的差值，能够反映气象因素对于水储量的影响，具有一定的参考价值。

另外，在南水北调输水前，GRACE观测的水储量亏损速率大约是气象因素造成的水储量亏损速率的2倍，它们之差可能体现了人为抽水这个因素。南水北调输水后，基于气象数据得到的TWS亏损速率基本保持不变，GRACE观测的总水储量亏损速率却降低了约一半，说明南水北调减缓了华北平原水亏损的情况，而气象因素导致的水储量亏损有所减缓但不是主要因素。

另一方面，假设南水北调输水前后，人们抽取地下水用水的速率不变，GRACE观测到的前后两个时段的水储量变化速率的差异(表 2最后一列)应该大于等于气象因素和南水北调输水速率变化的总和，因为增加的水量都应该被卫星探测到。但是表 2中GRACE总水储量速率差异却略小于南水北调年输水速率。其原因可能是“南水”送来后人们用水量略有增加，以及水储量的变化速率相较于输水有所滞后。

最后一个因素是，GLDAS数据产品是根据有限的气象数据集和土壤模型而得到的，其本身包含了模型误差；更重要的是它没有考虑人类活动如华北大面积地抽取地下水用于农业灌溉等因素。

4.2 华北平原地下水储量变化

为进一步分析南水北调的水对于华北平原地区水储量的补充是在地表还是地下，进而计算了地下水储量变化。在华北地区，总水储量变化包含土壤湿度(SM)变化、等效雪水(SWE)变化以及地下水(GWS)变化3个部分，可用公式表示为

$ \Delta {\rm{TWS}} = \Delta {\rm{SM}} + \Delta {\rm{SWE}} + \Delta {\rm{GWS}}. $

华北平原地表水体如人工湖和蓄水库等的变化对于长期的陆地水储量变化较小，将其忽略^[34]。对于土壤湿度及雪量，同样用GLDAS模型进行估计，取NOAH、VIC、MOS、CLM 4个模型的平均值。这些模型取的土壤层的最大厚度分别为200、190、300和343cm，远小于地下水的埋深。基于重力卫星的陆地水储量减去土壤湿度及雪量后可计算得到地下水储量变化。

如图 5(a)所示，华北平原内土壤及雪水储量在2004—2015年的变化趋势为(-0.28±0.08)Gt/a，呈现略微减少的态势，2015—2019年变化速率增加，变为(-0.10±0.37)Gt/a，两者相差0.18Gt/a。图 5(b)为地下水储量变化的时间序列，两段时间内变化趋势分别为(-5.64±0.48)Gt/a和(-2.87±2.67)Gt/a，说明南水北调通水后对地下水的恢复贡献达到2.77Gt/a，大约缓解了50%的地下水储量亏损趋势。因此，南水北调对总水储量的贡献在土壤水和地下水两方面均有显著的体现，它更多地补充了地下水的亏损，对地下水的恢复有重要作用^[35]。

4.3 北京市与天津市的水储量变化

南水北调中线的终点为北京的团城湖和天津的外环河。考虑到北京市和天津市受水区面积小而供水量大，南水北调对其影响效果显著，因而对这两个地区进行单独分析。结果如图 6所示。北京市2004—2015年陆地水储量变化速率为(-0.19±0.03)Gt/a，2015—2019年变化速率为(0.00±0.13)Gt/a，两者相差0.19Gt/a。天津市在这两段时间内陆地水储量变化趋势分别为(-0.29±0.05)Gt/a和(-0.02±0.21)Gt/a，两者相差0.27Gt/a。而南水北调的水量在这两市的变化速率分别为1.05和0.84Gt/a，表明南水北调水量与重力卫星所反映的水储量变化具有趋势上的一致性，但是GRACE所得到的结果只分别体现了北京18%与天津30%的南水北调贡献量。

至于为什么由GRACE得到的水储量变化小于南水北调的水量，可能有两方面原因：一是受GRACE空间分辨率限制，在小区域重力变化信号的恢复中信号泄露的现象更加明显，对北京、天津市水储量变化可能有一定的低估；二是GRACE数据2017年后部分月份缺失，特别是GRACE以及Follow-on观测数据间断1年多，二者的数据衔接可能不够准确。这个问题还值得进一步研究。

值得说明的是，北京和天津地区的空间尺度比较小，而GRACE的空间分辨率大约为300km，本质上反演如此小区域的水质量变化具有一定的困难。然而，由于本文采用变尺度因子方法，可以在一定程度上提高对小区域的真实信号恢复效果。另外，这两个区域位于华北平原的边缘处，与研究区域外的信号泄露因素可能会造成加大影响，北京和天津地区的结果应该含有一定误差。但是，本研究结果基本上反映了这两个地区的变化趋势，具有一定的参考价值。

5 结论与展望

本文利用卫星重力GRACE及GRACE Follow-on观测数据计算华北平原的水储量变化，并与南水北调工程供水量数据进行对比，评估南水北调工程对华北平原总体水储量、土壤水和地下水储量的影响。结果表明，南水北调工程的中、东线自2014年底通水后，其输送的总水量缓解了大约50%的华北平原陆地水储量亏损趋势，表明南水北调工程在缓解华北平原水储量亏损中起到了重要且积极的作用，同时也为华北地下水长期超采后的逐步恢复提供了有利条件。

该研究表明，相比于传统水井监测手段受限于空间观测点位覆盖率的问题，连续的重力卫星观测数据更能反映大尺度水资源迁移的时空变化特征，尤其在研究华北平原等大流域水资源变化问题上具有明显的优势，因为重力卫星具有连续观测并且精度高的特点。在GRACE数据处理过程中，恢复信号所采用的是空间域尺度因子方法，该方法是在假定研究区域内存在质量而区域外无质量分布情况下进行的，也就是假设研究区域外质量不发生变化。实际上，研究区域外也存在一定的质量变化，这些质量变化在数据处理过程中会泄漏至研究区域内，使得研究区域的计算结果有一定误差。考虑到华北平原远远大于GRACE空间分辨率，区域外的信号泄露误差应该比较小，其具体的量化结果有待进一步探讨和分析。还应该指出的是，由于GRACE Follow-on的观测时间短且与GRACE观测时段之间存在1年多的间断，2015—2019年间华北平原水储量变化趋势的计算存在较大不确定性。随着GRACE Follow-on观测数据的积累及其观测数据的进一步评估和优化，相信能更准确地估计华北平原水储量变化趋势，进而更准确地评估南水北调工程对华北平原水危机的缓解作用。