超导重力仪是重力测量最敏感、稳定的仪器，由于其具有秒采样速率，能够观测到1纳伽（1.0× 10^-11 m·s^-2）的重力变化，已成为探测地球内部重力场的重要工具，广泛用来确定引起重力场变化的环境因素，如降水、地震等^[1].通过对超导重力数据的分析，已经探测到一系列地震引起的同震变化^[1-3].国内外的研究表明，秒采样的重力仪还经常记录到大地震前数天之内的扰动异常^[4-8]，其变化过程有快有慢，持续或断续直至地震发生.例如，2001年昆仑山口西M8.1级地震和2004年12月26日印尼M9.3级地震前数天内都记录到重力扰动^[4].Lan等^[8]分析了日本的三台超导重力仪和新竹的超导重力仪数据，检测到与近10年间27个地震有关的35组震前重力扰动异常现象.汶川M_s8.0级地震前48h，武汉、成都等重力仪以及全国范围内的地震仪也检测到了这种震前的扰动异常^{[5-7, 9]}.与重力仪相比，宽频带地震仪采样率更高，二者完全独立观测，所获得的结果进一步说明了震前扰动异常的真实性、可靠性^{[5-7, 10]}.

目前，大多数研究所用数据来自单一台站或区域网，台站分布在某种程度上具有一定的方位性，若要获得对震前扰动特征更全面的认识，需从全球尺度通过对更多可利用数据的分析.2011年3月11日，日本M_w9.0级地震发生前数天内，中国重力台网的gPhone弹簧重力仪观测到类似的重力异常扰动^[11-12].与此同时，全球地球动力学项目（Global Geodynamic Project，GGP，http://www.eas.slu.edu/GGP）发布了3月份17个台站共20组超导重力仪的秒采样数据（图 1（a，b），其中BH有3组，WE有2组），为从全球尺度探测震前重力扰动并研究其特征其提供了重要数据.本研究对日本M_w9.0级地震前后7天（3月6-12日）的数据进行分析，检验震前是否存在重力的震前扰动现象，与此期间的地震活动进行对比（图 1（a，b）；图 2），进一步探讨该现象的形成机理.

图 1

Fig. 1

图 1 超导重力仪台站位置及3月6日-3月12日地震活动性 （a）全球台站及地震（M≥4）分布；（b）欧洲台站（图 1a中蓝色框部分）及地震（1.2≤M≤4.4）分布.红色五角星是UTC时间2011年3月11日05:46:23在日本东北发生的Mw9.0级地震. Fig. 1 Superconducting gravimeter stations and seismicity from 6 March to 12 March (a) Global stations and seismicity (M≥4);(b) European stations (inset region of Fig. 1a with frame in blue) and regional seismicity (1.2≤M≤4.4).Red star is the epicenter of Mw9.0 Tohoku-Oki earthquake that occurred at UTC05:46:23 on 11 March, 2011.

图 2

Fig. 2

图 2 地震活动与重力时间序列（以PE台站为例）的对比 （a）全球地震目录M-t图（黑色为日本及其周围的地震活动，即图 1a红色框部分）；（b）超导重力数据、重力固体潮及大气效应；（c）重力残差；（d）欧洲地震活动M-t图（图 1b）.Eq是指日本东北Mw9.0级地震的发震时刻，即UTC时间2011年3月11日05:46:23（下同）. Fig. 2 Comparison between seismicity and gravity time series (take PE station for example) (a)M-t diagram for global seismicity (earthquakes occurring in Japan and its surroundings, i.e.red rectangular region in Fig. 1a, are shown in black); (b) Superconducting gravity data, synthetic solid tide and pressure effect; (c) Residual gravity; (d)M-t diagram for European seismicity (Fig. 1b).'Eq'of abscissa denotes the original time of Mw9.0 earthquake, that is UTC05:46:23 on 11 March, 2011.

为检测震前扰动信息，需扣除固体潮效应及同步观测的大气压重力效应^[1]，却不需要进行潮汐分析，故而在计算中本研究没有采用如T-SOFT、ETERNA等程序^[13-14]进行处理、调和展开，而是采用GGP提供的扣除固体潮、气压的g1sres.f代码（http://www.eas.slu.edu/GGP/snm.html）^[15].该代码更加简洁明了，使用Xi^[15]完全调和展开的高达3070阶的潮汐分析结果，计算得到理论重力固体潮汐值.我们对其进行了改进并用来计算理论重力固体潮，由于高阶分量的影响微小，可忽略不计，实际计算中仅选取了383阶的分量，精确度达0.04×10^-8 m·s^-2.从超导重力数据中扣除理论重力固体潮汐值，获得固体潮汐改正重力.从固体潮汐改正重力中减去大气压效应，即实现大气重力效应改正.为简化计算，直接利用气压导纳（-0.3×10^-8 m·s^-2/ mbar）与气压的乘积确定大气重力效应；虽然可以通过计算获得更为可靠的气压导纳值，其结果并不影响最终的输出结果^[15].利用9次多项式拟合改正去除仪器的漂移及残余潮汐效应，得到重力残差时间序列（图 3）.

图 3

Fig. 3

图 3 地震活动与20组重力残差时间序列的对比 带圆圈、倒三角的短线分别为日本及周边、其他区域5级以上地震，图中震级放大5倍；重力残差的单位：1×10-8 m·s-2；左侧字符是台站名称，其中H3、H4、H5在德国的BH（Bad Homberg），W3、W5在德国的WE（Wettzell）（图 1）. Fig. 3 Comparison between seismicity and residual gravity time series of 20 SGs Lines attached with circle and inverse triangle denote M≥5 earthquakes occurring in Japan and its surroundings, other region, respectively, and the magnitude is enlarged by a factor of 5;residual gravity is in 1×10-8 m·s-2.

图 2a为3月6-12日期间地震（M≥4）的M-t图，图 2b展示了捷克共和国Pecny台站（PE）在地震前5天的观测记录、理论重力固体潮及大气压的重力效应，图 2c为经过固体潮改正、大气压改正的重力残差.图 2a的地震时间序列与图 2b重力观测曲线及图 2c重力残差曲线相比，该台站能观测到全球范围内大部分M≥6级地震所引起的重力高频波动信号.在PE台站周围发生了几个小震（图 1b），但结果显示该台站对近场（欧洲地区）3级小震的响应并不明显（图 2d）.

图 3为20组SG的重力残差变化与地震活动性的对比.除了NY台站在3月10日缺少数据没有记录到对日本3月11日大地震的响应外，其他19组数据都记录到3月11日地震引发的显著的高频波动信号.这20组超导重力仪数据对此前的6级以上地震有响应，如3月7日0:09:37在（-10.31°N，160.75°E）发生6.4级地震，3月6-7日在日本及其周边地区都发生了多次5级以上地震，20组数据也都有不同程度的响应.3月9日2:45:17在（38.52°N，142.84°E）发生7.3级地震，其响应特征在20组记录上都非常显著（图 3）.

图 4a是对PE台站震前5天数据的放大显示，重力残差曲线在震前3-4天（3月7日）开始明显增粗，即出现震前重力扰动现象，振幅最高达~10×10^-8 m·s^-2（图 4a）.对其进行功率谱密度分析，结果（图 4b）显示重力残差的能量峰值主要集中在0.01~0.1Hz频段，同时还有小部分能量集中在0.01Hz以下，表明该信号还包含潮汐成分.通过高通（f≥0.01Hz）滤波滤掉周期大于100s的成分，重力残差不再包含显著的潮汐效应（图 4c），其能量集中的0.02~0.10Hz频段是地震波的优势频段^[1]（图 4d），不过，该优势频段的截止频率不很明显.通过0.10Hz的高通滤波，能消除地震波的主要成分，更清晰地显示非地震波信号分量^[1]（图 5a）.不过，除了在0.10~0.15 Hz范围内的主要能量频段外（图 5b），还存在0.15~0.25Hz频段的次一级的能量峰值（图 5d），表明该频段信号的成分较为丰富（图 5c）.

图 4

Fig. 4

图 4 PE台站震前5天重力变化及其功率谱密度分析 （a）重力残差（图 2c时间截取Eq之前并放大显示）；（b）重力残差功率谱密度；（c）重力残差的高通（f≥0.01Hz）滤波结果；（d）重力残差高通滤波后的功率谱密度. Fig. 4 Gravity changes during 5 days before Mw9.0 earthquakes and the power spectral density (PSDs) (a) Residual gravity (enlarged for the time span before Eq as Fig. 2c); (b) PSD of residual gravity; (c) High pass filter (f≥0.01Hz) for residual gravity; (d) PSD of high pass filter.

图 5

Fig. 5

图 5 震前5天重力变化带通（0.10 Hz≤f≤0.25 Hz）滤波结果及其功率谱密度分析 （a）高通滤波（f≥0.1Hz）；（b）高通滤波功率谱密度；（c）带通（0.15~0.25Hz）滤波；（d）带通滤波功率谱密度. Fig. 5 Gravity changes during 5 days before main shock using band pass filtering and PSDs (a) High passed (f≥0.1Hz) gravity; (b) PSD of high passed gravity; (c) Band pass filter (0.15Hz≤f≤0.25Hz) for residual gravity; (d) PSD of band pass filter.

近年来，郝晓光课题组的研究结果^{[5-6, 16-17]}表明震前扰动的优势频段在0.10~0.18 Hz之间，而与台风有关的扰动的频段在0.20~0.25 Hz之间^[1].根据这一结论（或推断），对原始数据再进行两个频段的带通滤波分解.经反复尝试，带通（0.118~0.18Hz）滤波的结果能够较大程度地压制同震信号、并同时保留异常扰动信息；为保证信号的完整性，第二个频段取0.18~0.25 Hz.两组带通滤波结果分别如图 6a和图 6c所示，图 6b和图 6d分别为其相应的功率谱密度结果.这些结果表明两组带通滤波信号分量的能量至少相差一个量级.前者代表震前扰动信息，可能与震前缓慢的应力积累过程有关.相比而言，0.18~0.25Hz频段的信号，除了仍保留对个别地震的响应之外，变化比较平缓，其功率谱密度特征也很接近白噪声，即便不能排除与台风或海潮有关，至少与震前应力缓慢变化的关系不大.

图 6

Fig. 6

图 6 震前5天PE台站重力变化带通（0.118 Hz≤f≤0.25 Hz）滤波结果及其功率谱密度分析 （a）带通（0.118~0.18Hz）滤波Ⅰ；（b）带通滤波Ⅰ的功率谱密度；（c）带通（0.18~0.25Hz）滤波Ⅱ；（d）带通滤波Ⅱ的功率谱密度. Fig. 6 Gravity changes during 5 days before main shock using 2 band pass filtering and their PSDs (a) Band passed (0.118Hz≤f≤0.18Hz) gravity; (b) PSD of high passed gravity; (c) Band pass filter (0.18Hz≤f≤0.25Hz) for residual gravity; (d) PSD of band pass filter.

相对于M_w9.0级地震，有8个台站（11组）数据经过带通（0.118~0.18 Hz）滤波后能够观测到显著的震前扰动现象（图 7），其振幅从3月7日10时左右开始逐渐增大，到3月9日2:45:17在（38.52°N，142.84°E）发生7.3级地震时振幅增大更明显，波动起伏持续约12h逐渐下降，至3月9日23:37:00在（38.48°N，143.14°E）发生5.4级地震时，逐渐变平缓并一直持续到主震发生.对于其他台站，结果显示除了对个别地震的高频扰动信号较大外，其他时段的变化都非常平缓.

图 7

Fig. 7

图 7 震前5天20组SG重力变化带通（0.118 Hz≤f≤0.18 Hz）滤波结果与地震活动性对比（其他说明同图 3） Fig. 7 Band passed (0.118 Hz≤f≤0.18 Hz) gravity changes during 5 days before main shock and comparison to the seismicity (Others is the same as Fig. 3)

由于其较高的敏感度，超导重力仪的观测数据不可避免地受多种因素的影响，包括地震激发的地球自由振荡、地震背景噪声、水文现象、大气质量变化等等^{[1-3, 8, 15]}.在利用超导重力数据研究长期构造形变时，需要尽可能地将这些影响因素扣除.超导重力仪通常能够检测到频段在0.3~6 mHz、由地震激发的自由振荡^{[1, 18]}，在滤波处理中已去掉.

为扣除与大气压变化有关的重力效应，通常采用经典负荷理论，将大气压变化的负荷效应表示为相应的负荷Green函数与负荷压强的全球褶积积分^{[1, 8, 19-20]}.由于缺少秒采样的全球大气模型，本研究中计算大气压效应时仅以大气压导纳作乘积，简化了全球气压积分过程，可能造成了一定程度的偏差.韦进等^[21]利用九峰地震台SGC053超导重力仪540天的记录与气压观测数据获得气压导纳值，发现固定气压导纳值与频率依赖导纳值进行改正的结果相差不超过2纳伽.另外，从气压改正的幅度来看，大气压变化造成的重力效应通常在小于1微伽的量级，计算结果表明其在5天内的变幅不超过5微伽，与超导重力仪的观测值振幅相比显得很微小^[1].利用负荷格林函数区域积分的形式进行改正也不会有一个量级的变化，故而不足以影响观测值的趋势.

重力异常的主要影响因素是地下水的变化^[1]，根据空间尺度的不同，可分为局部的、区域的乃至全球的变化.降雨资料^[22]显示在地震前5天内欧洲出现过一次大范围的降水过程，不过主要集中在大西洋.由于该资料是日均降雨量^[22]，该降雨过程对重力仪数据造成的影响程度，还有待于对更高频采样数据的分析.一般而言，局部精确的地下水效应对超导重力仪观测值影响很大，强烈的降雨过程将造成重力的阶跃变化，而不是如图 7所示的震前振荡.另外，对台站周围（以及对数天内的数据）水文重力效应的改正是非常复杂的，需要对降水、土壤湿度、地下水位变化以及人类用水等造成的影响进行逐一排除.此外，考虑到所用台站是全球分布的，而海潮效应受台站位置的影响较大，并且已有研究表明海洋潮汐效应并不影响地震阶段的信噪比^[15]，在本研究中未予考虑.

M-t图显示，在主震发生前日本及其周边地区的5级地震显著增多，这些超导重力仪探测到大多数地震引起的高频波动信号.值得关注的是出现震前扰动的这8个重力仪都位于欧洲（图 1a），而在此期间欧洲范围内除了3月6日20:50:27发生的4.4级地震以外，其余地震的震级都小于4.0（图 1b，图 2d）.要探讨大震之前的扰动，需要确定观测到的高频波动信号是否与大地震之前的小震有关.

由于远场引起的密度变化非常微弱，地震引起的高频波动信号的幅值可以通过地震波引起的地表位移近似估计，故而信号的幅值近似与地震波能量成正比、与震中距呈反比^[8].以能量与震中距的比值作为标准，比较近场小震与远场大震所产生的扰动的幅值大小.震中距最大相差4个量级（1~20000 km），而震级对应的能量是指数变化的，震级相差两级则能量相差3个量级，于是震中距10km的3级地震与10000km的5级地震造成的重力扰动振幅相当，故而震中距小于10000km的5级地震要比震中距大于10km的3级地震造成的扰动强烈得多.由于欧洲小震几乎都在3.0级以下，其他地区5级地震到欧洲台站的震中距都小于10000km，故而这些5级以上地震对欧洲各台站的高频波动信号更显著.

在汶川地震发生前48h内，西太平洋上空恰好形成了一次台风（国际命名为Rammasun）.震前扰动似乎可以通过Longuet-Higgins^[23]提出的洋底压强变化激发地球微震模型进行解释^[24-25].后续的多数研究认为汶川地震之前观测到的扰动可能与这次台风引起的异常地脉动密切相关^[9].此后，Tanimoto^[25]利用简正模公式计算进行对比，发现该模型只对海洋深度浅于1km的有效.最近，通过对扰动信号时频特征的细致分析，郝晓光课题组^{[5, 6, 16-17]}最近的研究表明，中国大陆宽频地震仪在汶川地震前记录到异常扰动信号，既包含台风扰动也包含非台风扰动两种优势频率；其中频段在0.2~0.25Hz的扰动主要与台风Rammasum有关，在沿海地区较强而在内陆较弱；而另一种优势频率0.1~0.18Hz的扰动则与台风无关，其在震前10h出现急剧增强，靠近震中的地区较强而在沿海较弱^{[6, 15-16]}.利用多种仪器、长时间的连续观测结果表明，这种扰动信号是真实的并且是普遍的^[10].

最近的研究表明，日本M_w9级地震前的地震活动经历了加速-减速-积累-加速的迁移过程^[26].利用波形相关技术获得了精确的地震目录，Kato等^[26]确定了主震前的两个不同的前震序列（阶段）.它们分别以每天2km和10km的速度沿着海沟轴线向主震震中迁移，最终都在主震震中附近停止，期间还伴随着重复的小震.准静态滑动的时间过程表明，两个序列发生了慢滑移，并向主震的初始破裂点瞬态传播；第二个序列滑动速率高，可能造成了一定的应力加载，促进了主震不稳定破裂.前者开始于2月中旬止于2月底且持续约8天的平静期，后者开始于M_w7.3级前震之后且在12h后出现峰值.这一过程中可能经历了减速矩释放到加速矩释放的转变^[27]，是应力积累过程中两次成核过程的体现^[26].

Suito等^[28]利用连续GPS数据的反演结果，发现M_w9.0级主震前的M7.3级地震产生的震后形变及其总的地震矩都比同震响应大得多，这与通常的认识（即震后滑动量小于等于同震滑动量）正好相反，此前的数值计算结果也获得了类似的结果^[29].这个加速矩释放的时间段与Kato等^[26]的第二个序列的时间段相一致，不过GPS不能确定是否向震中传播.结合GPS、地震活动性的结果，推测该过程可能包含着与主震有关的预滑动^{[26, 30]}.

这些关于主震前的GPS反演结果以及地震活动序列的研究与我们获得的关于震前扰动的特征有许多相似之处，尤其是在慢滑移时间点上.比如，震前扰动大约从3月7日UTC10:13:08在（29.13°N，129.43°E）发生的M5.0级地震开始（图 7），这与Kato等^[26]第一序列后的平静期接近；震前扰动至3月9日发生的M7.3级地震时显著增强，此后波动变化持续约12h，这也与Kato等^[26]第二个序列的峰值点发生在该前震后的半天后一致.重力震前扰动（图 7）到3月9日UTC23:37:00发生5级地震后变为相对均匀的扰动，持续至主震发生，这与Kato等^[26]的第二阶段快速传播相一致.即便如此，这种震前扰动是否的确是由震前的慢滑移序列造成的，还有待进一步检验.

实际上，实验室及理论研究很早就已经发现地震（破裂）前存在着类似的成核过程，断层失稳前发生预滑并逐渐加速，当预滑速度达到动态破裂速度时裂纹扩展直至破裂，即在断裂带的围限区内，稳定的、慢速破裂生长会逐渐发展成不稳定的、高速破裂.大量的岩石破裂实验记录到破裂前的声发射异常扰动^[31]、温度场及热红外亮温场异常等现象^[32-33]，与本研究获得的震前扰动特征有许多相似之处，都反映了岩石（地壳）从成核到破裂经历了由慢到快的成核过程^[27].另一方面，孕震过程是一个极其复杂的过程，多数情况下，震前扰动还受构造带几何形态、介质性质、应力分布等多种条件的控制.近年来，根据对不同几何结构的断层系变形破坏过程中应变、断层位移和声发射事件的时空分布，对典型失稳事件的特征进行分析后，发现具不同几何结构的断层系有不同的变形物理场演化图象，其震前扰动特征有明显差异，并非所有失稳都能观测到扰动，对于同一失稳也并非所有的部位都能观测到震前扰动^{[31, 34-36]}.最近，Lan等^[8]的观测似乎更进一步证实了这种观点.这可能部分地解释本研究只在欧洲台站检测到震前扰动现象.

本研究首先对超导重力数据与地震活动性进行对比，结果表明超导重力仪能观测到全球范围内大部分M_w≥6级地震所引起的重力高频波动，欧洲地区发生的小震对附近超导重力仪的影响不及远场6级以上地震的显著.

通过滤波处理及功率谱密度分析，发现重力残差的地震波信息主要集中在0.01~0.1Hz频段，滤除地震波主要成分后，在0.10~0.25 Hz范围内仍存在两个能量相差一个量级的频段.根据前人的结果及反复尝试，确定了两组优势频段分别为0.18~0.25Hz和0.118~0.18Hz.前者的频谱特征与地震、台风的关系都不大；而后者主要反映了震前扰动信息，不过该扰动信号是否与孕震断层的亚失稳阶段的预滑动或裂纹扩展过程有关，仍需进一步研究.

目前对震前扰动的认识还受观测条件的限制，尤其是对异常发生时间及其持续时间的确定还未达成一致意见^[26-27]，这有赖于更多观测手段、更高分辨率数据的获取及分析.台风等大区域事件容易识别，其他局部或潜在的自然事件却很难排除，这无疑增加了对震前扰动解释的不确定性^[10].相比而言，在加载和观测设备等可以控制的实验过程中，也仍能观测到这种现象，至少说明震前扰动与破裂成核过程关系密切^[32-37].因此，若能结合实验记录、地震仪、连续GPS观测等多手段进行综合研究，有可能获得更多的新的认识.

本研究采用的超导重力数据以及合成重力固体潮代码均来自国际地球动力学合作计划（GGP，http://www.eas.slu.edu/GGP/tohoku2011.html），地震目录来自欧洲-地中海地震中心（www.emsc-csem.org），两位匿名审稿人提出了建设性的修改意见，在此一并表示感谢.