长期而稳定的高精度重力潮汐观测作为探测各种地球物理信号的有力工具，已经有效应用于地球内部结构(Hinderer and Crossley，2000)、板块构造(Francis et al.，2004; Sato et al.，2012)、冰后期回弹(Lambert et al.，2001; Sato et al.，2006; Steffen et al.，2009)、地球内部动力学(徐建桥等，2001; 孙和平等，2004; Shiomi，2008)等研究中.而超导重力仪(Superconducting Gravimeter，SG)是目前精度最高、连续性和稳定性最好的重力观测设备，观测精度可达0.01 μGal(1 μGal=1×10^-8 m·s^-2). 1997年全球地球动力学计划(Global Geodynamics Project，GGP)的组织实施(Crossley et al.，1999)，其目的之一就是研究与重力变化相关的各种地球物理现象.

然而，超导重力仪观测的重力信号除了包含地球物理和地球动力学信息之外，也会受到海洋潮汐、气压变化与地下水变化等环境因素的影响.因此，在重力观测的处理中进行环境因素的改正显得越来越重要.对于最主要的海洋与大气负荷效应问题，国内外许多学者都进行了深入的研究(Farrell，1972; Agnew，1997; 孙和平，1997; Boy et al.，2002; 罗少聪等，2005; Zhou et al.，2013)，改正的方法和流程都趋于成熟.但是，地下水变化对重力场观测的影响研究还不够成熟.一些学者曾研究了地下水变化对重力观测的影响(Pool，2008; Leiriãoet al.，2009)，简单的理论研究表明，在布格平板模型中，1 m厚的水层质量变化约能导致42 μGal的重力扰动(Crossley et al.，1998; Knudsen and Andersen，2002).而实际上影响地下水文运动的因素有很多，包括降雨、降雪、蒸发、土壤渗透、地表径流等自然因素，也包括地下水抽取、农作物灌溉、水库蓄水等人为因素(张人权等，2011).如果只针对其中的单个因素着手研究，往往不能准确估计该因素对重力变化的影响.如Harnisch和Harnisch(2006)在研究中用到了降雨、地下水位和土壤湿度数据，但只简单分析了各个影响因素与重力观测变化的相关性，没能 将这些数据结合起来描述地下水的流动，因而缺乏物理基础.而Kazama等(2012)从基本的地下水渗流方程(Richards方程)出发，解决了这样的问题，利用降雨和地下水位数据从物理上模拟出地下水的垂直分布情况，进而估计地下水变化导致的重力效应.

1997年11月，经厂家更新换代后的SG C032成功安装在位于武汉市郊九峰的大地测量与地球动力学国家野外科学观测研究站(简称九峰站)，其高质量的连续重力观测资料为研究该区域重力潮汐和非潮汐变化特征以及地球动力学问题奠定了良好的资料基础(许厚泽等，2000；孙和平等，2005；徐建桥等，2014).同时，九峰站是中国大陆唯一的国际重力潮汐基准站，是正在实施的国家重大科技基础设施建设项目“中国大陆构造环境监测网络”的基准站之一，是不同类型重力仪对比观测以及检验中国大陆的固体潮和海潮模型的重要实验和研究基地.张为民和王勇(2007)曾利用绝对重力仪研究了武汉九峰站的测点环境噪声干扰水平和重力变化速率，指出地下水变化对该站的重力观测影响不容忽视.徐建桥等(2008)也曾利用全球陆地水同化模型LaD和GLDAS研究了全球陆地水循环在武汉九峰站导致的重力变化，同时也指出由于没有考虑台站局部地区的地下水实际分布，超导重力仪观测的重力季节变化与全球陆地水循环导致的重力变化之间存在着较大的差异和大约55天的时间延迟.

本文目的是研究九峰站局部地区地下水变化对超导重力观测的影响，根据武汉国际气象交换站降雨量与九峰站地下水位记录，利用水动力学法从物理上模拟计算了地下水的渗透，得出台站附近区域土壤含水率随深度的分布，进而获得地下水变化对重力观测的影响.利用高精度的超导重力仪观测残差与模拟的地下水重力效应对比，可以检验模拟计算的实际效果；也可以对超导重力残差进行地下水重力效应的改正，去掉与降雨和季节性水循环有关的短周期信号，提高我们检测与地壳运动有关的长期重力变化的能力.

九峰站位于武汉市洪山区石门峰北面的山脊上，山体相对高度差35 m左右.上覆地层为中更新世残积、冲洪积物，下伏基岩由南向北分别为二叠系栖霞组(P_1q)灰岩、二叠系孤峰组(P_2g)硅质岩、三叠系(T)灰岩，地质构造较为稳定.该站始建于20世纪80年代末，在重力观测室建造时，为了减小环境温度变化对仪器的影响，将山体靠北一面挖开，采用整体框架结构，然后进行回填的方式建造，保证观测资料的可靠性和高精度要求(张为民和王勇，2007; 徐建桥等，2008).

九峰站所在武汉市属亚热带季风性湿润气候，总体气候环境良好.雨量充沛，近30年来，年均降雨量1269 mm，且多集中在6—8月.日照充足，年均气温16.7 ℃(据中国气象局公共气象服务中心1971—2000年资料统计).

为了研究地下水变化对九峰站重力变化观测的影响，台站钻取了一口水井，井内安装一台地下水位监测仪，用以记录地下水位的变化.水井距超导重力仪大约130 m，由于研究区域的水文环境相对简单，因此水井的水位可以近似看作台站局部区域的地下水位.图 1的黑色曲线表示水井的水位观测记录，起始面选在地面以下10 m的位置.观测从2008年5月开始至2012年11月结束，原始采样间隔为10 min，重采样间隔为1 h.由图可以清楚看到，4年多的数据足以反映出地下水位的季节性和年际变化特征.

日降雨量和蒸发量数据来自中国气象科学数据共享服务网的《中国地面国际交换站气候资料日值数据集》，该数据集包含了中国194个地面气象站1951年以来的气象资料，包括日平均气压、平均气温、蒸发和降水量等8个要素的观测值.图 1的灰色柱条表示武汉气象国际交换站的日降雨量时间分布，如图所示其降雨量季节性变化明显，多集中在6—8月.

图 1

Fig. 1

图 1 台站观测的地下水位和日降雨量 Fig. 1 Groundwater level measured at station and daily precipitation

九峰站SG C032从1997年11月开始，连续观测至2012年7月，考虑到与地下水位数据之间的时间一致性，本研究选取了2008年5月至2012年7月的超导重力仪观测记录，包括重力和气压数据.重力观测以1 s的采样率进行数据采集，观测精度达到0.01 μGal，仪器的漂移量为2.28 μGal/a(陈晓东，2003).超导重力仪的重力观测数据经过格值标定、修正尖峰间隔跳跃、重采样、扣除潮汐信号、气压改正、极移改正和去漂移项的处理过程之后，剩余的重力残差中主要包含的就是地下水重力效应.

SG数据的预处理采用国际地球潮汐中心推荐的重力固体潮观测资料预处理软件Tsoft(Van Camp and Vauterin，2005)，通过移去\|恢复的人机对话方式，直观剔除并改正观测资料中的干扰信号(如尖峰、间隔、跳跃、大地震干扰等)，并采用低通数字滤波器将原始每秒采样资料变换为每小时采样.然后，利用Eterna3.30标准分析软件(Wenzel，1996)对其进行调和分析，精密确定九峰站的重力潮汐参数和大气重力导纳值(-0.36 μGal/hPa)，进而在SG观测中精密剔除重力潮汐信号和局部大气负荷效应.极移改正用到的极坐标来自国际地球自转服务(International Earth Rotation Service，IERS)网站(http://www.iers.org).仪器漂移的去除采取了线性拟合的方法，因为安装一段时间(4至6个月)之后仪器的漂移基本上会变成时间的线性函数(Van Camp and Francis，2007).

格值转换后的原始超导重力仪的重力观测数据如图 2a所示，其中2012年5月1日出现一处46.79 μGal的跳跃，可能是由仪器掉格所致，采用Tsoft的跳跃修正器加以修正.图 2b的实线所示是经过尖峰间隔跳跃修正、扣除潮汐信号、气压改正、极移改正之后的重力残差曲线，虚线所示是漂移线性拟合项.去掉漂移项之后，得到的最终残差曲线如图 2c所示，为方便起见，未作特殊说明的情况下，后文出现的超导重力残差均指图 2c所示的残差.

图 2

Fig. 2

图 2 超导重力残差 (a)超导重力仪重力观测值；(b)未经漂移改正的重力残差值(实线)和线性漂移项(虚线)；(c)最终的重力残差值. Fig. 2 Superconductivity gravity(SG)residuals (a)SG observations；(b)SG residuals without drift correction(solid)and linear drift(dash)；(c)Final SG residuals.

研究表明，重力仪的观测对近区的地下水运动尤为敏感(Abe et al.，2006; Hasan et al.，2006; Kazama and Okubo，2009; Hector et al.，2013).Kazama和Okubo(2009)的研究指出，Asama火山观测站重力仪观测到的地下水重力效应能被70 m左右范围内的地下水质量变化较好地解释，而如此小的主要影响范围与台站所处的地形和地下水的埋深较浅有关.因此，本文在不考虑地形的水平起伏和地下水的水平流动的情况下，以监测水井的水位代表重力仪局部地区的地下水位，并模拟了非饱和层(地面至地下水位)的渗透过程，以便更精确地确定局部地区地下水流动对重力观测的影响.

假设土壤性质是均匀的，在非饱和区域内取一边长分别为Δx、Δy和Δz的均衡单元体，根据质量守恒原理，地下水渗透满足如下的连续性方程(薛禹群，1986；薛禹群和谢春红，2007；Bear，1979)：

(1)

式中：ρ_w为土壤水密度；v_x、v_y和v_z为三个方向上的渗流速度分量；θ为单元体的含水率

不考虑土壤的弹性形变和水的密度变化，Δx、Δy、Δz和ρ_w为常数，且不考虑水平方向上的流动，则有垂向一维运动的非线性扩散方程，即Richards方程：

(2)

式中：D和K分别为垂向渗透系数和扩散系数，均是含水率θ的函数.

将降雨量数据和地下水位数据分别作为上、下边界条件，利用有限差分法求解上述Richards方程，便可以得到各个时刻的地下水分布.而地下水变化对重力观测的影响则由如下的公式给出(Van Camp et al.，2006; Kazama et al.，2012)：

(3)

式中：g_w(t)为总的地下水重力效应，Δz是水平含水层的厚度.在实际计算中我们取万有引力常数G=6.67× 10^-11 m³·kg^-1·s^-2，水的密度ρ_w=1000 kg·m^-3.

注意到上述地下水渗透的模拟过程和地下水重力效应的计算完全没有用到重力观测值来作为任何程度上的约束条件，因此这里所采用的计算可以说是不依赖于重力观测的.

我们将地表至地下10 m的土壤划分为100个厚度Δz=0.1 m的等厚水平层，模拟了各层含水率随时间的变化，时间间隔Δt=3600 s.其中，根据常见的土壤类型和常见松散岩土孔隙度参考值(张人权等，2011)，并参考张为民等(2001)在武汉地区的土壤样品试验结果，取土壤的孔隙度为n=0.24.图 3显示了几个典型深度(地面以下1~9 m处)的含水率时间变化，9 m以下的土壤一直处于饱和含水状态，因此不作绘制.各深度处的含水率变化趋势一致，都是在冬季达到最低峰值(2011年除外)，在夏季达到最高峰值，这与图 1中的降雨量时间分布和监测水井水位的变化相吻合.而2011年的土壤含水率时间变化的最低峰值和水井最低水位出现在夏季，这主要是因为2011年武汉遭遇冬春连旱，1—5月降雨总量与历年同期比较偏少60%，一度达到严重干旱标准.6月发生强降雨，单月平均降雨量与历年同期比较偏多80%，形成旱涝急转态势(2011年武汉市水资源公报).

对比图 1和图 3可以发现，当大降雨发生时，各层的含水率随即显著增加，之后一段时间内，各层的含水率先后出现指数下降.但是各处含水率的变化并不是同步进行的，相对浅层而言，深层的含水率变化存在着一定的时间延后，这很可能是跟地下水的渗透有关.

图 3

Fig. 3

图 3 模拟计算的地下1~9 m处含水率变化 Fig. 3 Simulated soil water content variations at 1~9 m depth

图 4显示了三个典型时刻土壤含水率随深度的 变化.其中图 4a和图 4c处于水位相对较低的枯水期，水位埋深分别为8.15 m和7.23 m，图 4b处于图 4a和图 4c之间的丰水期，水位埋深为1.36 m.

图 4

Fig. 4

图 4 含水率随深度的变化 Fig. 4 Variation of soil water content with depth

图 5比较了模拟计算的地下水重力效应和超导重力残差，黑线显示的是2008年5月至2012年7月的超导重力残差，蓝线显示的是模拟计算的2008年5月至2012年11月的地下水重力效应.对比两者的变化趋势即可发现，模拟计算可以很好地重现出观测重力残差的季节和年际变化，在降雨较少的冬季达到最低峰值，在降雨较多的夏季达到最高峰值(2011年除外).经过快速傅里叶变换计算两者的年际变化振幅谱和相位谱(图 6)，图 6的振幅谱结果表明两者有着非常相似的频率\|振幅分布.两者在时域和频域上趋势的一致性说明经过文中各项改正之后的重力残差的主要信号来源是局部区域地下水的变化，同时也验证了本文模拟结果的正确性.

图 5

Fig. 5

图 5 模型估计的地下水重力效应(蓝)和SG重力残差时间序列(黑) Fig. 5 Simulated groundwater hydrological effects on gravity(blue)and the time series of SG residuals(black)

图 6

Fig. 6

图 6 模型估计的地下水重力效应(蓝)和SG重力残差(黑)振幅谱及相位谱 Fig. 6 Amplitude spectrum and phase spectrum of the simulated groundwater hydrological effects on gravity(blue) and the SG residuals(black)

模拟计算的结果表明武汉九峰站局部区域地下水变化最大可以引起15.94 μGal的重力变化，在冬夏降雨量差别较大的两季之间，地下水重力效应的变化通常会达到10 μGal以上.因此，在进行精度较高的重力观测时，特别是在时间跨度超过几个月或者几个季度的情况下，必须将局部区域地下水的影响考虑进去.

虽然，模拟计算的地下水重力效应在时域和频域上都表现出了与超导重力残差一致的变化趋势，但是，两者的振幅和相位也存在一定差异.如整个时段内，超导重力残差的变化幅度为12.73 μGal，而模拟的地下水重力效应的变化幅度为15.94 μGal，两者相差3.21 μGal.在两者振幅最大的周年频段上，超导重力残差的振幅和相位分别为：2.96 μGal和-69.0°；模拟的地下水重力效应振幅和相位分别为：4.18 μGal和-46.5°.模拟的年际变化振幅比观测到的大1.22 μGal，相位上提前22.5°.

能够引起两者差异出现的可能原因有许多，如超导重力残差中可能含有非环境因素的地球物理信号；超导重力资料的处理过程中也可能带入不确定因素.这都是与模拟计算无关的影响因素.与模拟计算相关的影响因素包括：(1)地形起伏的影响.由于缺少台站区域的地形资料，我们对地下水渗透模型作了地形的简化处理，即不考虑台站附近区域的地形起伏，将其作为平坦地区进行研究.(2)模型简化的影响.文中采用模型只考虑了垂直方向上的一维地下水渗透，忽略了地下水在横向上的流动，且不考虑土壤的不均匀性，即假设不同深度的土壤性质相同，这显然与实际的情况有所差异.(3)水文资料的不足.地下水的分布和流动对时间和空间都有很强的依赖性，比如地形的起伏、植被的繁茂程度和地表建筑物的遮挡等，都会对地下水的分布产生影响，仅凭一口水井的观测资料难以反映出台站附近区域的地下水详细分布.(4)与邻近区域的地下水交换.地下水系统中土壤含水层除了受到来自降雨的补给和源于蒸发作用的消耗之外，还会与邻近区域的含水层发生地下水交换.

本文基于一维地下水渗透模型，以地面降雨量和地下水位等观测资料为约束条件，模拟获得了武汉九峰台站局部区域土壤含水率随深度变化的分布特征，并进一步估计了其导致的重力变化.从结果来看，局部区域地下水变化最大可以引起15.94 μGal的重力变化，因此，在进行时间跨度较长的精密重力观测时，必须对其作相应的地下水改正.

在时域和频域分别对模拟计算的地下水重力效应和超导重力残差进行对比，发现两者在趋势上均有很好的一致性，这说明超导重力残差的主要信号来源是局部区域地下水的变化，同时这种一致性也验证了本文模拟结果的正确性.此外，超导重力残差能够很好地反映出地下水重力效应随时间变化的季节性和年际特征，因此，采用超导重力残差可以获得台站局部地区的平均地下水分布，可以为了解地下水的渗透过程提供重要约束.

然而，模拟值和观测值之间依然存在着振幅和相位上的偏差，这需要加入地形和更加详细的水文观测资料来提高模拟结果的精度.

感谢中国科学院测量与地球物理研究所大地测量国家野外科学观测研究站和中国气象局的数据支持，以及两位审稿人给文章初稿提出宝贵的评论意见和修改建议.