地球重力场能够反映出地下物质的分布及其运动状态，是现代地球科学解决人类面临的资源、环境和灾害等紧迫课题的基础地学信息。在地壳运动过程中，由各种原因引起的区域构造单元应力场的强弱变化都会在重力图像中展示^[1]。流动重力观测和连续重力观测的时间序列都能反映出观测点附近区域的重力变化情况，其中，连续重力观测可以精确探测到由地球内部各圈层的物质迁移引起的重力变化效应，包括叠加在一起的长期和短期变化，是研究固体潮汐和地震前兆异常的重要手段^[2]。一般连续重力观测的计算精度会略高于流动重力观测，从理论上来说，流动重力观测得到的地震前重力场变化信息连续重力观测同样也能观测到，而连续重力信号将理论重力潮汐变化剔除后得到的重力非潮汐变化信号，能够更加直观地反映台站周围垂直位移和地下物质的运移情况^[3]。

沈阳台PET型重力仪重力非潮汐变化自2017-06出现速率明显增大的趋势^[4]，变化量值约为5 μGal/月，截止到2017-09累计变化达到15 μGal。从重力时序变化和潮汐变化来看，该仪器观测稳定，而同时段的牡丹江台连续重力观测数据出现与沈阳台相反的变化趋势(图 1, 来自中国地震局地震研究所韦进)。沈阳台重力观测点位于辽宁东北部地区郯庐断裂带的中北段，区域构造断裂发育，历史上曾发生多次中强地震，如1765-03-15沈阳故宫M_S5.5地震、2003-03-30沈阳东陵区M_S4.3地震和2013-01-23辽阳灯塔M_S5.1地震。本文主要探讨沈阳台重力非潮汐的变化趋势与构造活动、地下介质及周边环境变化的关系，以期为相关地震预报提供可靠依据。

1 沈阳台重力观测概况

沈阳重力基准台位于沈阳市东陵公园南麓山洞(41.82°N，123.50°E)，地质上位于辽宁东部平原地区^[5]，年均降雨量为600~1 000 mm，降雨主要集中在每年7~9月份，山洞洞口高于最高洪水位和地下水位。台站周围地层结构简单，为前震旦系鞍山群及第四系上更新统地层，地表混合花岗岩出露较好，重力仪就放置在以混合花岗岩为基岩的摆放墩上。

沈阳重力台的PET型重力仪为美国Micro-g LaCoste公司生产的进口仪器，仪器精度为10 μGal，秒采样，2007-08正式运行。重力仪观测资料的内在精度指标^[6]显示，仪器运行状态稳定，数据质量较好。表 1为仪器的噪声水平、各月份Nakai检验小于2段数的百分比及M₂波(重力固体潮汐的M₂波潮汐因子)的中误差情况。

2 重力非潮汐变化分析 2.1 数据预处理

重力仪器在进行重力场变化观测时，不可避免地会将干扰信息记录下来，同时由于仪器自身或供电问题，长周期连续记录中会混杂数据的突跳或断记，因此在计算地球潮汐常数前应对观测数据进行预处理，在不影响结果精度的前提下尽可能地将各种干扰删除。本文采用Tsoft数据预处理方法，首先利用消去法去除观测残差中的各种干扰信号(即对残差中不同类型的干扰信号使用不同的方法修正)，顾及大气重力负荷效应，利用恢复法将消去干扰信号后的残差恢复(即将修正后的残差加上台站重力潮汐模拟理论值，也叫合成潮)，并对漂移进行拟合，得到消去干扰信号的观测值，即台站重力固体潮非潮汐量。

2.2 重力非潮汐变化

自2007年沈阳台连续重力观测以来，观测数据较为稳定，年变周期波动正常，每年6~8月雨季时处于年变低值，12月到次年2月处于年变高值状态。2015-05沈阳台重力非潮汐出现持续下降的趋势，并于6月开始转折持续上升，至同年7月开始整体呈持续平稳状态，在此次上升下降过程中变化幅度累积高达30 μGal。该现象是由于2015-05-27~05-31电源故障断电，开机后重力仪内温没有达到工作温度，导致仪器读数不稳定，当仪器内温达到工作温度并保持稳定后，仪器的读数恢复正常。而2017-06开始，M₂波非潮汐因子数据明显增加，速率为5 μGl/月，到2017-11连续增加5个月，打破自2015年断电后回归的稳定状态(图 2(a))。而相对应的M₂波潮汐因子的变化一直稳定，在非潮汐变化出现趋势性加速过程中，M₂波仍维持正常量值的小幅波动变化，并未有同步变化(图 2(b))。

2.3 重力非潮汐变化干扰因素分析

沈阳台连续重力自观测以来的资料^[4]显示，由于台站地理位置的特殊性，重力变化基本不受当地降雨、气压和大风、雷电等天气干扰因素的影响。对辽宁地区2016-05~2017-09共4期流动重力观测数据进行0.5 a尺度的分析发现：1)沈阳地区的重力场在2016-05~08这0.5 a间处于正值变化(图 3(a))；2)2016-08~2017-05，重力由将近-30 μGal左右变为20 μGal，正负变化将近50 μGal(图 3(b))；3)2017-05~09，重力又发生正负值反转，正值达10 μGal，在这1.5 a的时间里，沈阳周边近100 km范围内重力由正值转为负值，整体上变化幅值不大，但沈阳地区一直处于正负交替状态(图 3(c))。重力背景场变化与沈阳台连续重力非潮汐的趋势性变化有一定的关联性，用流动观测结果进行约束，重新对2016年初至2017-10沈阳台PET连续重力观测数据进行非线性零漂改正，分析得到的重力残差时间序列是准确的，说明沈阳地区地下介质活动有所起伏，这与沈阳台重力非潮汐变化自2017年下半年开始速率增大相对应，也从流动重力背景场变化说明了沈阳台重力非潮汐趋势变化的可靠性。

3 垂向沉降及抽水对重力非潮汐变化的影响计算

沈阳台距离居民区较远，除台站院内花圃园灌溉用水外，周边无其他工业井、农业井用水等情况。在距离山洞口大约40 m处建有苗圃花棚，该花棚距离山洞内重力仪约70 m。2017-05在原花圃区域种花，并且抽取台站地下水进行浇水灌溉至同年10月末。由于台站基岩为混合花岗岩，岩层表面及裂隙中含有褐铁矿，在花圃园抽水后会使岩石的孔隙水减少，影响介质密度分布，对重力潮汐变化产生一定的影响。

为了准确识别沈阳台PET型重力仪非潮汐变化趋势性加速产生的原因^[7-8]，并对可能产生干扰的各项因素进行详细分析，对同台GPS观测资料的垂向变化及花圃灌溉抽取地下水对局部重力观测产生的影响进行研究，并在分析中采用数值计算将影响信息进行量化，使沈阳台重力非潮汐异常变化的分析结果更加真实。

3.1 GPS基线垂向重力变化量值分析

考虑到冬春两季水负荷对重力的影响不明显^[9]，利用2017-06-01~2017-11-30的时变重力过程(变化总量为-48 μGal)和沈阳台GPS站(沈阳测点)的垂直位移时变过程(变化总量为28 mm)研究两者之间的关系。

沈阳台GPS国家基准站时间序列自2017-05垂直分量出现下降，下降幅度略高于往年(图 4)。该基准站同样受周边用水干扰(图 5)，排除抽水干扰对GPS基准站的影响后发现，GPS垂直位移可以反映该地区垂直方向的变化特征。2017-06开始，重力非潮汐变化速率出现较往年明显增大的趋势，计算变化量值约为5 μGal/月。由于重力变化在地面垂直梯度上发生相对应的变化，根据大地测量中垂直位移和重力变化规律，对应变化量值的理论值为-302 μGal/m，计算2017-06~11的GPS垂向变化数值和同时段重力潮汐变化对应的量值^[10-11]。

根据公式计算2017-06-01~11-30 GPS变化量与重力点值的变化：

$h_{v}=1.27 \times 10^{6} t+0.022 $

(1)

$g_{\mathrm{res}}=4.3 \times 10^{4} t-1.396 $

(2)

计算结果表明变化值为-1.714 μGal/mm，即2017-06~11这0.5 a内，沈阳台重力垂直梯度变化值为-1.714 μGal/mm。

分析该时间段内定点重力潮汐变化量值与垂向位移线性拟合比值及相关系数的对应变化量可知，与大地测量中垂直位移对应的重力变化数值差异较大。然而，该时间段内的短期重力变化和地壳垂直位移之间的明显负相关规律与上述理论规律一致(地壳下沉，重力增加)，进一步表明该地区地壳垂直位移的变化态势为重力量值处于上升趋势，地壳整体处于下沉变化趋势，两者之间的关系与理论值的量化差异说明垂直形变的变化不足以影响到连续重力观测出现异常变化。

3.2 苗圃及花棚用水模拟计算重力变化

沈阳台在2015~2017年未灌溉，仅在2017-04~06进行抽水灌溉，该时段用水量约为260 t，与重力值变化时间段相符，计算该时段抽水造成的重力变化。重力学中正演模拟任何形状的外重力影响公式为：

$\Delta g(p)=g \iiint\limits_{v} \frac{\rho(q)(\zeta-z)}{r_{p q}^{3 / 2}} \mathrm{d} \xi \mathrm{d} \eta \mathrm{d} \zeta $

(3)

建立苗圃和花棚灌溉引起的重力变化模型，总水量为260 t，水的密度为1×10³ kg/m³，重力观测点到抽水灌溉点最大距离约40 m，观测点以北均为山体，多为花岗岩。可将苗圃灌溉看成一个立方体，将其归结为长方体重力正演模型，以x₁、x₂、y₁、y₂、z₁、z₂为边界的长方体重力正演结果为：

$\begin{array}{*{20}{c}} {\Delta g = G\Delta \rho ||| - x\ln (y + r) - y\ln (x + r) + }\\ {\mathit{z}arctan\frac{{xy}}{{zr}}|_{{x_1}}^{{x_2}}|_{{y_1}}^{{y_2}}|_{{z_1}}^{{z_2}}} \end{array} $

(4)

设定(y₁，y₂)为(0 m，40 m)，(x₁，x₂)为(-40 m，40 m)，(z₁，z₂)为(-1 m，0 m)。

根据cuboids_gravity(x₁, x₂, y₁, y₂, z₁, z₂, density)程序计算，且x₁=-40、x₂=40、y₁=0、y₂=40、z₁=-80.081 25、z₂=-80、ρ=10³ kg/m³，最后得到的计算结果为Δg=0.218 1 μGal，即2015~2017年的抽水对重力变化的影响仅能形成0.218 1 μGal的重力变化。经过量值计算，排除抽水影响，重力非潮汐变化值为14.781 9 μGal，仍为趋势性加速变化。

此外，重力场变化受抽水引起的地下水位影响主要表现为局部地区水位变化导致重力变化^[12]，沈阳台出现的地下水位下降应为区域性变化，且抽水引起地下水位下降的幅度较小，故认为抽水引起的地下水位的变化并不能引起15 μGal的重力变化，排除干扰后确定该重力变化为前兆异常现象。

4 结语

1) 2016-09~2017-09三期0.5 a尺度的流动重力数据表明，沈阳地区重力正负反转变化较大，该地区地下介质活动有所起伏，与重力仪的非潮汐变化趋势性上升相一致。

2) 同台站GPS基准时间序列垂直分量出现下降，垂向变化引起的重力变化量远小于重力仪非潮汐变化量，说明该重力非潮汐变化出现异常；经过模拟计算2015~2017年末由花棚灌溉用水引起水位变化产生的重力值变化量为0.218 1 μGal，远小于沈阳台重力非潮汐自2017-06~09共计15 μGal的变化量，排除抽水影响，该非潮汐重力变化量可达14.781 9 μGal。通过GPS垂向对重力影响的计算和抽水干扰模型的计算分析后发现，沈阳台重力非潮汐数据短期加速上升为前兆异常。

3) 对于连续重力潮汐因子的趋势性变化，不仅需要对仪器情况进行分析，还要对GPS垂向沉降引起的重力变化进行分析，并将两者量值进行对比；在对比周边环境产生的干扰时，需选用合适的数学模型进行定量计算及对比分析。经过上述分析和计算，使确定异常的可靠性增强，为地震预报提供了更为准确的判定依据。

致谢: 感谢中国地震局地震研究所韦进老师提供所需的部分资料及在本文撰写中的特别指导，异常落实工作得到了沈阳台有关人员的大力协助，在此一并表示感谢！