0 引言

次声波在传播过程中不易被吸收，容易发生衍射而绕过障碍物。地震台站记录的连续地震波形数据中也包含次声波，通过频率域分析可以将低频信号提取出来。低频信号异常特征表现在两方面，一方面为频谱特征的频带变窄并向低频偏移，另一方面为地脉动低频持续波动现象，又被称为“前驱波”、“慢滑动”、“间歇性震颤”、“低频率地震”、“震颤异常波”、“相关脉动”和“临震微波动”等(赵根模等，2001；唐林波等，2002；蒋骏等，2009；杨立明等，2018a、2018b)，该现象在2008年汶川8.0级、2010年青海玉树7.1级和2014年新疆于田7.3级地震前均有相关研究(杨立明等，2009、2015)。早在多年前就有学者(冯德益等，1983、1994)认为，大震前短周期面波的频谱可能会在震前显示出一定的前兆异常特征，其物理本质可能是断层的预滑或慢破裂(陈运泰等，1979；郭增建等，1991)，而之后的破裂成核理论(Dieterich，1992)中具有前兆意义的推论也说明在震前会出现低频波。近些年，很多研究工作致力于利用被动源成像研究地壳和岩石圈变形特征，这对了解中强地震深部孕震环境具有重要意义(沈旭章，2013)。目前，探测岩石圈结构成熟的方法为接收函数法(Phinney，1964；Langston，1979)。接收函数法主要有P波接收函数和S波接收函数两种。其中用P波接收函数获得地壳上地幔结构时，由于莫霍面和壳内间断面多次反射波的干扰而无法较好地确定岩石圈边界，而S波接收函数不受间断面多次反射震相的干扰，在确定岩石圈-软流圈边界具有一定优势(Farra et al，2000；Yuan et al，2006)，但在实际应用中，能够获得稳定可靠的P波接收函数的数目却十分有限。近年来发展起一种将勘探地震学方法应用于天然地震学结构研究的新方法——自相关方法，即利用自相关从透射波响应中提取出反射波信号。该方法基本原理早在20世纪70年代就已被提出(Claerbout，1968)，其核心思想是垂直入射的透射波自相关等价于反射波。自相关技术已被广泛用于探索澳大利亚的地震构造，并在不同方面进行了改进(Kennett et al，2018)。

2008年汶川8.0级地震发生后，关于龙门山断裂带的研究工作积极开展，包括从P波、S波远震接收函数偏移成像结果中揭示汶川地震震源区下方地壳岩石圈变形特征(沈旭章，2013)，从地震测深、大地电磁测深、小震双差定位分析、精细三维P波速度模型等方面进行研究，发现震源区下方中下地壳及上地幔顶部存在切割莫霍面向东倾斜的隐伏断裂带或壳幔韧性剪切带(蔡学林等，2004、2008)，并提供出地壳流沿龙门山断裂带上浸等多项地震发生依据(雷建设等，2009；Deng et al，2011)。

本文在利用远震P波及其尾波自相关探测汶川8.0级地震震中所在区域地下结构变化特征的同时，分析连续地震波形中低频异常信号与震中地下结构变形之间的关系。

1 地质构造背景

2008年5月12日14时28分4秒，汶川8.0级地震发生在龙门山断裂带上，作为世界上发生在陆内的最大逆冲型地震之一，其震源深度仅约14km，形成了迄今为止空间上分布最为复杂、长度最大的逆冲型同震地表破裂带。龙门山断层属于逆冲型断层，其中的一个特点就是上盘为主动盘、下盘为被动盘，因此上盘会摇得特别厉害，由此造成近7万人死亡、36万多人受伤、18万多人失踪、500多万人无家可归(李西等，2008)。龙门山构造带位于青藏高原东缘边界山脉，而青藏高原东缘不仅是中国东西部一级地貌边界带，也是中国最重要的地震活动带，被称为南北地震带，龙门山构造带恰好位于我国南北地震带的中段，是最重要的地震活动带之一。从板块构造理论的观点来看，龙门山构造带位于扬子板块的西北部，西部为松潘甘孜褶皱带，北部为秦岭板块俯冲碰撞造山带的米仓山地体和大巴山逆掩构造体，西南部为康滇地轴。从构造带分布特征上看，龙门山构造带北起广元，南至天全，长约500km，宽约30km，呈NE—SW展布，北东与大巴山相交，南西被鲜水河断裂相截，构造带主要由一系列大致平行的NE—SW向的逆冲断裂组成，由西至东分别为：汶川—茂县断裂(也称后龙门山断裂)、映秀—北川断裂(也称龙门山中央断裂)、安县—灌县断裂(也称龙门山前山断裂)(Deng et al，2011)(图 1)。映秀—北川断裂带是汶川地震的发震断裂，形成约270km的地表破裂，是汶川地震中最长的地表破裂带(Xu et al，2008；Fu et al，2011)。

2 原理及方法 2.1 频谱包络

地震计记录的以1h为窗长的连续波形，在未记录到地震、爆破等明显事件时，其频谱基本上应该处于一个稳定的形态，且其相应的频谱包络幅值最大值也基本处于一个水平状态波动。为获取频谱包络幅值最大值的时序变化特征，在对连续波形处理时关键要进行两个步骤，一是通过快速傅里叶变换得到1h连续地震波形频谱，二是寻找频谱包络幅值最大值(图 2)。将所得频谱图中不同频率的振幅最高线连接起来形成的曲线称为频谱包络线，频谱包络线在0~5Hz的最大值简写为MASE(Maximum Amplitude of Spectrum Envelope)，即连续波形数据中需要被提取的第二个特征量，表达式为

$ x_{\mathrm{i}}=\frac{a_0}{2}+\sum\limits_{k=1}^m\left(\sqrt{a_k^2+b_k^2} \cos \left(\frac{2 {\rm{ \mathsf{ π} }} k i}{N}+\arctan \left(-\frac{b_k}{a_k}\right)\right)\right) $

(1)

$ a_0=\frac{2}{N} \sum\limits_{i=0}^{N-1} x_i, a_k=\frac{2}{N} \sum\limits_{i=0}^{N-1} x_i \cos \frac{2 {\rm{ \mathsf{ π} }} k i}{N} $

(2)

$ b_k=\frac{2}{N} \sum\limits_{i=0}^{N-1} x_i \sin \frac{2 {\rm{ \mathsf{ π} }} k i}{N}, k=1, 2, 3, \cdots, m $

(3)

$ c_k=\sqrt{a_k^2+b_k^2}, \phi_k=\arctan \left(-\frac{b_k}{a_k}\right) $

(4)

其中，{x_i}为1h记录到的连续波形数据序列，c_k为振幅谱。

2.2 自相关提取反射波

地震计记录到的远震数据W(ω)包含震源信息S(ω)、震源至间断面的传输路径信息P(ω)、间断面至自由表面的多次波发生的反射和透射信息C(ω)，其中震源信息可忽略不计，则一个无干扰的地震道远震记录在频率域可以表示为

$ W(\omega)=P(\omega) C(\omega) $

(5)

由自相关与卷积的关系可知，对W(ω)做自相关后有

$ \begin{aligned} A(W(\omega)) & =\left[W^{\mathrm{T}}(\omega)\right] *[W(\omega)] \\ & =C^{\mathrm{T}}(\omega) * P^{\mathrm{T}}(\omega) * P(\omega) * C(\omega)=C^{\mathrm{T}}(\omega) * C(\omega) \end{aligned} $

(6)

当远震P波作为透射波在垂直入射至水平层时，震源信息作为点源可以忽略不计，透射响应T(t)与反射响应R(t)的关系为

$ \frac{T(t) * T^{\mathrm{T}}(t)}{\varPi \mathrm{det}}=\left(1+R(t)+R^{\mathrm{T}}(t)\right) $

(7)

构建台站叠加自相关图来提取地壳反射信息的过程，在一定程度上能够提高地壳反射系数的一致性并压制噪声。

采用自相关方法获取以台站为虚拟震源和接收点的地震反射波，处理流程主要有预处理、自相关、动校正和叠加四个步骤(图 3)，由于该方法是利用垂直分量记录的远震地震事件P波震相以及全球震相PKP和PKIKP等及其尾波构建反射波(Ruigrok et al，2012)，其中最关键的步骤是动校正，用来消除由入射慢度引起的地震反射相位的走时差，对射线进行时差校正，使其垂直入射(Kennett，2001)。

根据Taylor公式展开，即

$ \begin{aligned} (1+x)^m & =1+m x+\frac{m(m-1)}{2 !} x^2+\cdots+\frac{m(m-1) \cdots(m-n+1)}{n !} x^n+\cdots \\ & \approx 1+m x+\frac{m(m-1)}{2 !} x^2 \end{aligned} $

(8)

时差动校正量具体确定如下：假设界面深度为h₀，上覆介质速度为v₀，震源位于地表时垂直反射双程时间τ₀=2h₀/v₀，令式(8)中的x=-v₀²p²、m=1/2，得出入射慢度为p的时差动校正量为-v₀²p²/2，τ-P域反射时间可表示为

$ \tau=\tau_0\left(1-v_0^2 p^2\right)^{\frac{1}{2}}=\tau_0\left(1-\frac{1}{2} v_0^2 p^2+\cdots\right) \approx \tau_0\left(1-\frac{1}{2} v_0^2 p^2\right) $

(9)

3 数据处理及结果分析

数据处理过程中使用了两种类型数据，第一种数据是汶川8.0级地震周边的9个固定台站记录的连续地震波形数据，另一种是位于汶川8.0级地震震中附近的中国科学台阵所记录的远震事件波形数据。连续波形数据用来提取频谱包络幅值最大值的时序变化特征，而远震地震事件波形数据用来提取震中下方地下结构特征。

3.1 连续地震波形

选取2008年汶川8.0级地震附近9个地震台站(JYA、JMG、AXI、DFU、JJS、YZP、WCH、MXI、HSH)的连续地震波形数据进行频谱分析(图 4)，以1h为研究时间窗长，发现0~5Hz频谱形态在不同时段发生明显变化(Liang et al，2018)。

将每小时段连续地震波形数据的0~5Hz频谱幅值包络最大值作为异常信号跟踪对象，在这里记作MASE(Maximum Amplitude of Spectrum Envelope)。在2008年4月1日—2008年5月16日研究时段，JYA、JMG、AXI、DFU、JJS、YZP 6个台站的MASE变化形态非常一致(图 4(a)、4(b))。此外，低频异常信号不仅在汶川8.0级地震主震发生前的5~7天(2008年4月17日19时—2008年4月19日12时)出现高幅值脉冲，且在主震发生后持续多天(2008年5月12日14时—2008年5月16日23时，至少4天)的MASE高频现象也非常一致，持续多天呈衰减状态符合地震能量释放过程。但是，由于数据处理截止时间为2008年5月16日23时，无法判断低频异常信号在震后持续的准确时长，但必然是大于4天，也能展示出地震能量逐渐衰减的过程。除此以外，由于WCH、MXI、HSH 3个台站在汶川8.0级地震发生后数据中断，无法看到与JYA、JMG、AXI、DFU、JJS、YZP形态一致的同震响应(图 4(c)、4(d))。然而，WCH、MXI、HSH 3个台站在汶川8.0级地震发生前阶段性(5~7天)出现的低频信号高幅值脉冲与JYA、JMG、AXI、DFU、JJS、YZP 6个台站记录的时间及形态非常吻合。特别指出的是，2008年4月17日19时至2008年4月19日21时约40小时内也出现了MASE先升后降的异常变化形态，该40小时内时间段MASE异常形态可能反映“微破裂、慢滑移”过程特征。

3.2 远震地震事件波形

根据现有慢地震、微破裂等孕震机理，本研究进一步来探究连续地震波形中出现的地震波低频异常信号，是否与监测台站下方的结构变化有关系。为此，在探测汶川8.0级地震震中区域下方地下结构情况时，依然采用与MASE跟踪频率段相一致的数据(0~5Hz的地震波形数据)，数据来源于中国地震科学台阵在汶川8.0级地震震中附近所记录的远震地震波形(图 5(a))。采用距离台阵30°~90°范围，时间尺度为2006年11月—2007年12月、2007年5月—2008年5月、2008年5月—2009年8月、2009年9月—2010年12月4个不同时段的远震地震事件(图 5(b))。具体使用数据为P波前10s至后100s的远震地震波形数据，再进行去仪器响应、提高信噪比等操作后，在频率域采用[0.8，2]Hz、[0.5，4]Hz、[1，5]Hz进行带通滤波后做自相关计算。最后，将不同慢度的远震事件结果进行动校正后叠加来削弱震源对结果的影响程度。

由于2008年汶川8.0级地震发生后震中附近很多地震计遭到破坏，数据停记，选取了一条震前和震后数据记录相对完整的测线进行研究(图 5(a)中震中东侧白色三角形台站组成的条形测线)，对该测线进行地下剖面探测，在主震发生前后不同研究时段查看其地下结构的情况。被选取的台站名称由北至南依次为：KMX03、KMY06、KMY09、KCD01、KCD03、KCD06、KCD10、KLS03、KLS15、KLS14，其中KMY06和KMY09仅在2007年5月—2008年5月之间存在数据。测线剖面结构结果显示在2006年11月—2007年12月、2007年5月—2008年5月、2008年5月—2009年8月、2009年9月—2010年12月4个不同研究时段内，地下结构形态发生变化，特别是2009年9月—2010年12月表现最为显著(图 6)，KLS03与KLS15地下20~40km处在4个不同时段内结构形态差异明显。再对震中周边的几个台站(KCD05、KWC05、KWC03、KYA08、KYA09)进行不同时间段的跟踪，发现距2008年汶川8.0级地震震中最近的KCD05和KWC05(图 7、图 8)在震后2008年5月—2009年8月结构变化非常明显，并且KCD05在震前2007年5月—2008年5月的结构相较于2006年11月—2007年12月也发生了较为明显的变化，特征为莫霍面上移，可以理解为软流圈的热物质上涌，待地震发生后，由于震源处能量释放，使得传输到地表的能量转化为逆断抬高式地表破裂现象，而向下的能量在震后使得岩石圈纵向变形。

汶川8.0级地震地表破裂带由两支NE向同向破裂带和一支NW向同震破裂带组成，其中两支NE向同震破裂带沿中央断裂和前山断裂发育。值得注意的是，在2008年5月—2009年8月时段KCD05比KWC05变形特征明显，结合修正后的汶川8.0级地震地表破裂带分布(图 1)来看，KCD05恰好位于前山断裂靠近震中位置，也在沿灌县—安县断裂发育的汉旺—白鹿破裂带的西端；而KWC05台站位于后山断裂(非发震断裂)，也就是在汶川—茂汶断裂并靠近震中的位置，但汶川8.0级地震发生时并未出现地表同震破裂，虽然KCD05和KWC05均在震中附近，但台站下方地下结构在2008年5月—2009年8月变形程度却完全不同。同样，同在后山断裂(汶川—茂汉断裂)上的KWC03比KWC05在2008年5月—2009年8月时段的变形特征明显，虽然KWC05比KWC03更靠近震中，但KWC03相比KWC05变形反而表现得更强烈，这可能是由于主震发生后引起后山断裂的汶川—茂汉断裂东端结构不稳定，成为余震进一步发生的优势区域，而KWC05即使距汶川8.0级地震震中近，却不在后山断裂余震优势发震区域上。从实际发生的5级以上余震看，也证实即使KWC05距离汶川8.0级地震震中更近一些，但5级以上余震却没有发生在KWC05位置附近，所以主震发生后KWC05没有记录到非常明显的地下结构变形信息，与震后余震分布及地表破裂结果相一致。此外，尽管KWC03与汶川8.0级地震震中相距较远，却在主震8.0级地震发生后，在KWC03所在位置附近先后发生2008年5月12日19时11分汶川M_S5.6、2008年5月13日19时54分汶川M_S5.0、2008年5月14日22时54分汶川M_S5.2地震，余震的发生也解释了8.0级主震发生后，KWC03地下结构处于不稳定状态(图 9)。

最后，位于8.0级地震震中西南侧的KYA08、KYA09在2009年9月—2010年12月结构形态相较于2006年11月—2007年12月、2007年5月—2008年5月、2008年5月—2009年8月三个时段表现得更不稳定。对同一台站，这种不稳定在[0.8，2]Hz、[1，5]Hz、[0.5，4] Hz不同频段均表现一致(图 10、图 11)。

4 结论与讨论

实验室中的“岩石局部失稳—岩石结构破坏”过程，在实际应用中，从密集台阵记录的背景噪声数据中捕捉关于包含震源信息的低频信号及地下结构方面进行验证：含有震源信息的低频异常信号，其产生的原因是由监测区下方结构的变化所引起的，从异常现象到结构变化两方面来佐证地震发生的紧迫性。本文依据岩石破裂声发射实验思路，开展了从2008年汶川8.0级地震波形数据中提取低频异常信号，得到其在震前、震后的形态特征，并通过探查震中附近地下结构在四个时间段的形态特征，寻找低频信号与地下结构的关系，由此得到以下主要结论：

(1) 本文在思路上通过地下结构在不同时段的变化来解释地震前低频异常现象，是现象与地下结构两方面共同开展的震例回溯性研究；在数据使用上，对于同一地震计记录的地震数据，既使用包含震源信息的背景噪声数据，又使用能够反映地下结构的远震数据；在方法选用上，地震波低频异常信号提取技术是将岩石力学实验室探索工作转移到地球内部这样一个大的实验场中，由于地球是一个各向异性的粘弹介质，使用被动源自相关技术时不需要震源信息，对地震数据进行重建就能获取有效的地下介质信息。

(2) 在对2008年汶川8.0级地震回溯性研究中，主震发生后持续出现低频异常信号的同时，从地下结构中也能看到震前、震后的明显变化，但是在2008年4月17日19时至2008年4月19日21时约40小时内，震中附近9个台站提取出的地震波低频信号先升后降的异常变化在本次结构信息中无法看到，这需要使用间隔更加密集的数据进行处理才有可能看到变化。但是由于提取低频信号与探索地下结构所使用的数据是同频段的(0~5Hz)，故震后低频异常信号密集出现的现象可以确定是由地下结构变化引起的，这个变化可能不仅仅是地壳的变化，也有可能是整个更深部岩石圈整体的变化。

天然微震与实验室岩石破裂等实验有相似之处，实验中高应力作用下的岩石出现“微小破裂萌生(微尺度、微能量)—裂缝扩展、贯通(小尺度、小能量的破裂)—岩石局部失稳(中尺度、中能量的冲击)—岩体结构破坏、大范围失稳(大尺度、大能量的释放)”的发展过程在今后地震数据的高效处理下可能得以实现，希望未来在时域、频域双向实时监测条件下开展基于地震前兆现象与地下结构变化的研究。

致谢: 感谢中国地震局地球物理研究所“中国地震科学探测台阵数据中心”为本研究提供地震波形数据(doi：10.12001/ChinArray.Data)(中国地震科学台阵，2006)，感谢四川省地震局提供的连续地震波形数据，感谢各位老师的耐心指导。