十三陵地震台gPhone重力仪的仪器性能与水文响应分析

佘雅文 付广裕 韦 进

1 中国地震局地震预测重点实验室（中国地震局地震预测研究所），北京市复兴路63号，100036

2 中国地震局地震研究所（地震大地测量实验室），武汉市洪山侧路40号，430071

地下水变化与降雨是影响重力观测的重要因素之一。近年来水文影响对连续重力观测的研究已经从经验性转变为物理模型［1］。Tanaka使用两台gPhone获取的重力残差时间序列较好地反映了重力变化和降雨的响应关系，并用Kazama建立的水文模型较好地模拟了重力变化［2－3］。利用降雨和静水位的变化来确定长期定点重力观测的水文响应模型，对去除干扰、获取当地地壳内部物质迁移以及深部构造活动具有重要意义［4－5］。

本文利用十三陵地震台的两台gPhone重力仪（序号109、118）2013－04～08 连续观测数据进行潮汐分析和提取重力残差处理，从而对仪器的性能进行判定，并在此基础上分析gPhone连续重力观测的水文响应特征。

1 连续重力观测

十三陵地震台位于山区向平原过渡的斜坡地带。两台gPhone重力仪并排放置在近似水平的平整瓷砖上，仪器旁放置数据采集和网络通信设备，仪器和地面的接触较为平整。利用两台重力仪在同一观测条件下得到的数据来验证数据的可靠性，排除仪器自身原因造成的数据失真。由于断电等原因，2013－06－19～07－04数据缺失。

2 数据处理

仪器原始数据为秒采样，利用加汉宁窗的数字低通滤波器将秒采样数据滤波为分钟采样，使用Tsoft软件以480为窗口数、12cpd为采样率将分钟采样降为小时采样。采用Tsoft预处理软件［6］对观测曲线中的尖峰、地震和掉格等干扰信号进行手工修正，对较短的记录中断用三次多项式进行插值填补，保留长时间中断，并把秒采样数据转换成1min至1h采样，对小时值数据进行极移校正和水平校正。

利用BAYTAP－G 软件进行数据处理分析。BAYTAP－G 软件可以将重力数据分离为随机噪声、固体潮汐、响应数据项和漂移项4个部分，同时给出对数据潮汐因子的估计值［7］。使用辅助气压、传感器温度、水平校正值和相对重力数据一起作为输入，其中，辅助气压、传感器温度和水平校正值为辅助数据。气压和传感器温度是重力观测的重要影响因素，气压的影响大于传感器温度的影响。传感器温度虽然是自相关的，但其突变会引起观测数据μGal级的变化，即每0.001 ℃可能会引起2.6μGal的重力变化［3］。选取水平校正作为辅助项是因为，虽然仪器本身的预处理中包括水平校正，但是仪器的倾斜角度是随时间不断变化的，由仪器记录的水平校正数据来看，影响在几十μGal，且变化无规律。BAYTAP－G 在分离数据的过程中并不会因为辅助数据在原始数据中不存在而产生错误结果，因而忽略辅助数据的影响［7］。因为gPhone的漂移不是线性的［3］，所以选择二阶最小二乘法拟合观测数据。对分离出的漂移信号作二阶最小二乘法拟合，用漂移项数据减去拟合数据即可得到重力残差［1］。

3 数据分析

3.1 两台gPhone重力仪观测数据的对比分析

图1给出了两台gPhone重力仪在2013－04－20～06－19观测到的预处理数据、辅助观测值和残差时间序列，其中1（a）为gPhone重力仪观测到的原始数据下采样后的小时值数据，1（b）为仪器传感器温度，1（c）为gPhone重力仪水平校正值，1（d）为仪器围压值，1（e）为重力残差值。水平校正值和传感器温度作为BAYTAP－G的辅助项得到的重力漂移再减去其二阶最小二乘法拟合值，得到重力残差时间序列（05－09 23：00～05－11 12：00之间gPhone109无测量数据）。从图1（b）和图1（c）可以看出，gPhone109和gPhone118虽然放在同一个观测室内，但两台仪器的传感器温度和水平校正项有较大的差异。前者在2013－04－20～06－20时段内传感器温度仅在53.463 ℃～53.465 ℃变化，且大部分时间维持恒温，而后者传感器温度在53.094 ℃～53.103 ℃之间呈阶梯状变化，可见118号仪器的传感器温度变化没有109号仪器稳定。109和118号仪器的水平校正项变化幅度分别为10μGal和35μGal左右，水平稳定性上109 号仪器也强于118 号。使用BAYTAP－G 分离出的109号和118号仪器漂移项通过一阶最小二乘法拟合后计算的漂移速率分别约为5.2μGal／d和5.9μGal／d，109号仪器的漂移速率比118号低，造成这种情况的原因可能是109号仪器的水平性能和传感器温度变化都比118号稳定。总体而言，109号gPhone重力仪的性能优于118号重力仪。

图2为两台重力仪在2013－04－20～06－20的观测数据处理得到的重力残差变化曲线，其中图2（a）为BAYTAP－G 处理两台重力仪数据得到漂移项减去其二阶最小二乘法拟合得到的重力残差，图2（b）为两台仪器二阶重力残差的三次样条插值。从图2（b）可以看出，重力变化的趋势总体是一致的，局部略有差异。由于两台仪器外部环境相同，所以差异来源于仪器本身。从图2（a）中可以看到，109号重力仪比118 号重力仪的重力残差变化要平缓，且基本看不到周期性变化，而118重力残差中可看到明显的周期变化成分。

图1 gPhone109和118号重力仪2013－04－20～06－20预处理数据、辅助观测值和残差时间序列Fig.1 Time series of preprocessing data，auxiliary data and gravity residual of gPhone meters（serial number109and 118）from 2013－04－20～06－20

图2 2013－04－20～06－20两台重力仪残差时间序列比较Fig.2 Comparison of time series of the gravity residuals of two gPhone gravimeters from 2013－04－20～06－20

3.2 两台重力仪记录数据的调和分析

选取109和118号gPhone重力仪2013－04－20～06－20的连续观测记录进行处理和分析。表1是109和118号gPhone重力仪记录的重力数据通过下采样为小时值后，使用BAYTAP－G［8］以及ETERNA33［9］进行潮汐分析的结果。

表1 BAYTAP－G、ETERNA33和DEHANT理论潮汐因子对比Tab.1 Comparison of tidal factors of BAYTAP－G，ETERNA33and DEHANT

从表1可见，两台仪器的数据使用不同处理软件得到的潮汐振幅因子及其中误差较为接近，证明了本文数据处理的正确性。为了提供对比参考，将BAYTAP－G 和ETERNA33的振幅因子与DEHANT 理论模型振幅因子［10］进行比较。可以看出，BAYTAP－G 处理得到的潮汐振幅因子整体较为接近理论值，其原因可能是BAYTAP－G可使用辅助数据抑制数据校正不充分在潮汐分析中的干扰。此外，109号仪器与118号仪器的潮汐振幅因子相比整体更接近于理论值。

综上，109号gPhone仪器整体性能优于118号仪器。由于这个原因，之后有关水文活动造成重力变化的讨论都将围绕109号gPhone重力仪的观测数据进行。

3.3 降雨造成的重力变化分析

图3（a）和图3（b）分别给出了2013－05－25～06－19和2013－07－05～08－04的109号gPhone重力仪重力残差变化和降雨的时间序列图。从图中可以看到，大部分重力变化与降雨存在较好的对应关系。大部分情况下，它们之间的对应关系满足降雨时重力变小、降雨后重力逐渐变大的规律。但是也有少数情况不符合这种模式，这可能是由于其他未知因素的重力变化引起的，其原因尚待进一步研究。

为解释这种重力随降雨变化的特征，建立简单的物理模型（图4（a））。图4（b）中箭头所指位置为109和118号gPhone重力仪的观测室，观测室西侧紧邻海拔为161m 的龙山，观测室的海拔为101m。图4（a）中，降雨开始时水层覆盖在地表，由于龙山海拔高于观测室，在龙山上的雨水对重力仪产生引力倾斜向上，与重力仪处于同一高度的雨水对重力仪的引力相互抵消，这时重力值变小。降雨结束后，由于山体的坡度，雨水更多地流向地势低的地方，渗透到低于重力仪海拔的地下，这时重力值增大。Naujoks使用超导重力仪也观测到了类似现象［4］。

图3 gPhone109二阶最小二乘法拟合残差与降雨对应的时间序列Fig.3 Correlation of time series of 2nd－order curve fitting residual gravity of gPhone meter（serial number 109）and rainfall

图4 重力随降雨变化的物理模型示意图和仪器位置与周边环境Fig.4 Response model of gravity to rainfall and the position and surrounding environment of gPhone meter

3.4 重力变化与地下水静水位变化的对应关系

图5是119号gPhone重力仪观测结果的二阶最小二乘法拟合残差和沙河地震台静水位关系图。沙河地震台与十三陵地震台的距离约为15 km。由图5可知，连续重力观测残差曲线与降雨量之间存在较好的相关性，水位下降时重力减小，水位上升时重力增大。可见在对连续重力观测进行分析时，地下水的影响不可忽视，因此建议在全国的连续重力观测网中配备地下水位观测。

4月左右的沙河静水位下降在重力观测信号中无显示，可能是由于gPhone重力仪在03－24才启动观测，尚未完全稳定所致。由图5 可知，gPhone重力观测大约需要1～2个月时间来稳定信号，之后才能进行较好的分析。109号gPhone重力仪观测到的重力变化和沙河地震台的静水位存在4d左右的延迟对应关系。延迟的原因可能是地下水流动造成的，由于地下结构比较复杂，地下水的活动受很多地质环境的影响，如地下水流动、断层的阻断等，都会造成重力变化对水文活动的延迟［4，11］。

图5 gPhone二阶最小二乘法拟合残差与沙河地震台静水位Fig.5 Time series of 2nd－order curve fitting residual gravity and static water level of Shahe seismic station

4 结 语

本文对十三陵地震台的两台gPhone重力仪（109和118号）在2013－04～08的连续观测数据进行潮汐分析和提取重力残差处理，对仪器的性能进行判定，并在此基础上分析gPhone重力观测的水文响应。结果显示：1）影响gPhone重力仪观测精度的仪器自身因素主要来自传感器温度和仪器的水平稳定性。2）利用BAYTAP－G 软件和ETERNA33软件获得的潮汐分析结果大体一致，证实了两套处理软件的一致性。相对而言，BAYTAP－G 得到的主要振幅因子更接近DEHANT 给出的理论值。3）依据109号gPhone重力仪观测结果给出的二阶重力残差变化和十三陵地震台的降雨有μGal量级的对应关系，并且降雨开始时重力值减小，降雨结束后重力值增大。该现象可根据观测站周围特定的地形条件定性解释。4）109号gPhone重力仪获取的二阶重力残差变化和距离十三陵地震台15km 的沙河地震台静水位变化存在4d左右延迟的正相关对应关系，该对应关系显示了地下水位观测对重力观测解释的必要性，而4d左右的延迟则是由于重力观测与地下水位观测点存在15km 的差异所致。重力变化与地下水位变化之间的正相关对应关系表明，在全国性连续重力观测网中，有必要进一步配备地下水位观测，以提高连续重力观测的数据处理精度。

［1］Tanaka T，Miyajima R，Asai H，et al.Hydrological Gravity Response Detection Using agPhone Below－and Aboveground［J］.Earth，Planets and Space，2013，65（2）：59－66

［2］Kazama T，Okubo S.Hydrological Modeling of Groundwater Disturbances to Observed Gravity：Theory and Application to Asama Volcano，Central Japan［J］.Journal of Geophysical Research：Solid Earth，2009，114（B8）

［3］Tanaka，T.Introduction／Adjustment of gPhone and Its Present Condition of the Data［J］.Rep Tono Res Inst Earthq Sci，2010，25：75－81

［4］Naujoks M，Kroner C，Weise A，et al.Evaluating Local Hydrological Modelling by Temporal Gravity Observations and a Gravimetric Three－Dimensional Model［J］.Geophysical Journal International，2010，182（1）：233－249

［5］韦进，李辉，刘子维，等.武汉九峰地震台超导重力仪观测分析研究［J］.地球物理学报，2012（6）：1 894－1 902（Wei Jin，Li Hui，Liu Ziwei，et al.Observation of Superconducting Gravimeter at Jiufeng Seismic Station［J］.Chinese J Geophys，2012（6）：1 894－1 902）

［6］Camp M V，Vauterin P.Tsoft：Graphical and Interactive Software for the Analysis of Time Series and Earth Tides［J］.Computers ＆Geosciences，2005，31（5）：631－640

［7］Tamura Y，Sato T，Ooe M，et al.A Procedure for Tidal Analysis with a Bayesian Information Criterion［J］.Geophysical Journal International，1991，104（3）：507－516

［8］Akaike H.Likelihood and the Bayes Procedure［J］.Trabajos de Estadística Y de Investigación Operativa，1980，31（1）：143－166

［9］Wenzel H G.The Nanogal Software：Earth Tide Data Processing Package ETERNA 3.30［J］.Bull Inf Marées Terrestres，1996，124：9 425－9 439

［10］Dehant V.Tidal Parameters for an Inelastic Earth［J］.Physics of the Earth and Planetary Interiors，1987，49（1），97－116

［11］Pool D R，Schmidt W.Measurement of Ground－Water Storage Change and Specific Yield Using the Temporal－Gravity Method near Rillito Creek，Tucson，Arizona［R］.Tucson：US Geological Survey，1997