地球物理学报  2010, Vol. 53 Issue (4): 853-863   PDF    
利用背景噪声互相关研究汶川地震震源区地震波速度变化
刘志坤 , 黄金莉     
中国地震局地震预测研究所, 北京 100036
摘要: 利用2007年3月至2009年3月四川数字地震台网的宽频带连续波形资料,通过计算地震背景噪声互相关提取台站对间的经验格林函数,在0.1~0.5 Hz频带下测量每天经验格林函数与参考经验格林函数的走时偏移, 进而得到各台站对在该时段内的相对地震波速度变化.结果表明,2008年5月12日汶川Ms8.0级地震造成了震源区地震波速度的急剧降低,最大降幅达0.4%;大致以安县为界,余震带西南部地区在汶川主震后波速降即达到最大值,而东北部地区的最大波速降一般出现在主震后的1~4个月,相对地震波速度变化的这种分段特性与地震序列的时空分布特征有较好的对应关系;在震源区外围的四川盆地也观测到了震后波速降低,而川西高原内部则没有出现显著的波速变化.进一步的分析和计算结果表明主震的静态应力变化和强地面运动引起的地表破坏都不能很好地解释震后波速的急剧降低,地震导致的断层区内部结构破坏和周边介质应力状态改变可能是波速变化的主要原因.
关键词: 地震背景噪声      经验格林函数      地震波速度变化      汶川地震     
Temporal changes of seismic velocity around the Wenchuan earthquake fault zone from ambient seismic noise correlation
LIU Zhi-Kun, HUANG Jin-Chi     
Institute of Earthquake Science, China Earthquake Administration, Beijing 100036, China
Abstract: We detected the variations of crustal seismic velocity around the Ms8.0 Wenchuan earthquake fault zone during March 2007-March 2009 from ambient noise cross-correlation by using continuous seismic data recorded by the broadband stations of Sichuan digital seismic network. The data processing procedure was divided into three phases: (1) single station data preparation including time domain normalization, spectral whitening, (2) cross-correlation and temporal stacking of 31-day daily cross-correlation functions (CCF), and (3) estimation of relative velocity change by measuring travel time shifts in frequency 0.1~0.5 Hz between the 31-day stacks and the reference empirical Green function, which was the stack of 14 months of CCF before the Wenchuan earthquake. The preliminary results revealed a sudden post-seismic velocity drop for the station pairs across the Longmenshan fault, the largest of which was more than 0.4%. The seismic velocity changes with time exhibited a spatial difference: in the southwest segment of the Longmenshan fault, the maximum value of velocity reduction occurred at the time immediately after the Wenchuan earthquake, whereas in the northeast segment of the fault velocity reduction appeared 1~4 months after the main earthquake. Furthermore, the spatial extent and magnitude of the post-seismic velocity change in the Sichuan Basin exceeded that in the mountainous regions west of the fault zone. The static stress change and near-surface physical damage caused by strong ground motion could not entirely explain the measurements in this study. We think the temporal changes of seismic velocity maybe related to the damages from fault ruptures and stress changes around the fault zones..
Key words: Ambient noise correlation      Empirical Green function      Seismic velocity change      Wenchuan earthquake     
1 引言

2008年5月12日发生在四川盆地与川西高原交界处的Ms8.0级汶川大地震是带有少量右旋走滑分量的逆冲构造型地震,发震断层沿龙门山断裂带向北东方向延伸,破裂长度超过300km,地表产生了规模宏大的破裂带,地震发生后,沿发震断层的狭长条带又发生了大量的余震.此次地震造成了重大的人员伤亡和经济损失,因而引起了国内外地学工作者的极大关注,在地震的发震机理[1, 2]、破裂过程[3, 4]、地壳结构[5, 6]及形变[7, 8]等方面都进行了研究,取得了很多有意义的成果.

地震孕育及发生前后断层区的应力状态和介质属性会发生明显的变化[9~11],研究这些变化对于了解断层演化规律及地震危险性分析有着极其重要的意义.地震学家很早就致力于地震前后地壳介质的时变性质的研究,提出了包括测量波速比的变化[12]、剪切波分裂的变化[13]、尾波Q值的变化[14]及利用重复的人工[15]或天然震源[16]研究波速变化等地震学方法,但对汶川地震引起的地壳介质性质随时间变化的研究还很少.丁志峰等[17]研究了汶川地震序列剪切波分裂的时空变化现象,观测到靠近余震带南端的地震台站的快慢波到时差在主震后持续减小,而北端台站的到时差比主震前明显增大.但由于主震前龙门山断裂带的地震活动水平不高,他们的研究只用到了主震前的一个地震记录.

利用背景噪声互相关方法研究地壳介质地震波速度变化是近两年发展很快的一种新方法.对两个地震台站的背景噪声进行互相关运算可提取其间的经验格林函数,当噪声记录足够长时,不同时段的经验格林函数波形往往可以保持较高的相似性,如果介质中地震波速度发生变化,则这种变化可以通过测量不同时段经验格林函数波形的微小走时偏移探测得到.该方法比依赖地震发生的上述多种地震学方法有更高的时间采样率,较人工震源的方法经济得多.通过此方法,地震学家观测到了火山区波速的季节性变化[18]和火山喷发前数周内波速的下降[19],也研究了2004年日本新潟县中越地震[20, 21]、2003年圣西蒙和2004年帕克菲尔德地震[22]以及苏门答腊地区大地震[23]的同震和震后地震波速度变化.在国内,到目前为止还未见到利用背景噪声研究与强震相关的波速变化的任何结果发表.本文利用背景噪声互相关方法研究了汶川地震前后2007年3月至2009年3月震源区及周边的相对地震波速度变化,详细分析了波速随时间和空间的变化规律,并探讨了引起这种变化的原因.

2 数据和方法

本次研究应用了四川数字地震台网宽频带地震台的垂直分向连续波形资料.数据处理分为三步:(1)对单台数据进行预处理得到地震背景噪声;(2)通过计算台站对间背景噪声的互相关函数来提取经验格林函数;(3)计算不同时段的经验格林函数与参考经验格林函数的走时偏移,由此得到相对地震波速度的时间变化.

2.1 观测数据

本文主要研究汶川地震震源区及周边的地震波速度变化,故选用四川数字地震台网位于龙门山断裂带附近的台站(图 1)记录的连续波形资料.所用资料的时间范围为2007年3月到2009年3月,但离主震较近的汶川台(WCH)由于在地震中被损坏造成三个多月没有观测资料,错过了地震波速度变化的主要时段,故本次研究没有用该台站的资料;茂县台(MXI)和青川台(QCH)2008年5月12日到2008年6月2日缺测,也给研究结果带来一定影响.

图 1 所用数字地震台站分布 图中背景为研究区地形,三角形代表地震台站,五角星代表汶川主震,白色小圆圈表示重新定位的余震[24],黑色虚线为实测地表破裂位置[25],黑色实线为断层,其中,F1.茂县一汶川断裂,F2.北川一映秀断裂,F3.灌县一江油断裂,F4.青川断裂,F5.新津一成都一德阳断裂,F6.龙泉山断裂,F7.马尔康断裂,F8.岷江断裂,F9.虎牙断裂. Fig. 1 Distribution of seismic stations used in this study The grey base map shows the surface topography. The triangles denote seismic stations. The big star and small circles denote the main earthquake and aftershocks relocated[24], respectively. The black dash lines denote the co-seismic surface rupture zones[25] and the black solid lines denote active faults. F1. Maoxian-Wenchuan fault; F2. Beichuan-Yingxiu fault; F3. Guanxian-Jiangyou fault; F4. Qingchuan fault; F5. Xinjin-Chengdu-Deyang fault; F6. Longquanshan fault; F7. Markam fault; F8. Minjiang fault; F9. Huya fault.
2.2 单台数据预处理

单台数据预处理的目的是消除天然地震和仪器本身异常信号获得背景噪声[26].对于每个台站,只选用每天记录时间超过10h的连续波形资料,用0.05s时间间隔进行重采样并去除均值和线性趋势,再采用滑动窗口绝对值平均的方法对波形进行时域归一化,最后在0.08~1 Hz频带下进行谱白化处理,这样就得到了背景噪声.

2.3 每天经验格林函数与参考经验格林函数的计算

单台数据预处理后,对研究区内各台站对每天的背景噪声在相同时间内作互相关运算,其结果即为台站对间的经验格林函数[27].为提高信噪比,将各台站对当天及其前、后15天(共31天)的互相关函数叠加来代表该天经验格林函数,如果可用来叠加的互相关函数少于20天,则该天的经验格林函数不予计算.为了避免由于叠加而平滑掉5.12地震前、后的地震波速度变化,地震前每天经验格林函数只由地震前的互相关函数用上述方式叠加得到,地震后的经验格林函数计算方式亦然,由此得到了各台站对2007年3月到2009年3月每天的经验格林函数.我们发现通过这种方式叠加的绝大多数台站对间经验格林函数都比较稳定,但其信噪比会随台间距的增大有所降低,故本文一般只选择台间距小于150km的台站对进行研究.图 2a是由上述方式得到的AXI-PWU台站对2007年3月到2009年3月间的经验格林函数的图像,从图中可以看到清晰的面波信号和稳定的尾波信号,各相位在研究时间段内有很好的连续性,这说明从地震背景噪声信号中提取的经验格林函数是稳定可靠的.2008年5月12日的汶川大地震使经验格林函数发生了明显变化,特别是走时变化.为了更清楚地看到这种变化,将图 2a中流逝时间τ在25~45s部分进行了放大,得到图 2b,从该图中可以看到汶川地震(图中白色粗线所示)后,流逝时间40s处的相位走时有较大滞后,反映了地震造成了地壳介质属性的变化.

图 2 2007年3月至2009年3月AXI-PWU台站对每天经验格林函数及参考经验格林函数 (a)中横轴为日期,纵轴为经验格林函数的流逝时间,红色和蓝色分别表示每天经验格林函数归一化振幅的正值和负值,图的右边给出了归一化振幅的色标,彩色图中的波形是横轴所示日期(007年9月15日)的经验格林函数,彩图右侧的波形是由2007年3月1日至2008年5月11日的互相关函数叠加得到的参考经验格林函数.图(b)是对图(a)中流逝时间25〜45 s部分的放大. Fig. 2 The daily and reference EGF from March 2007 to March 2009 for station pair AXI-PWU (a) The horizontal and vertical axes represent the date and the lapse times of EGF (empirical Green function), respectively. The colors indicate the normalized amplitudes of daily EGF as shown in color bar. The waveform in the figure represents a current EGF (Sep 15, 2007). The waveform in the right box represents the reference EGF stacked using cross-correlation functions from March 1, 2007 to May 11, 2008. (b) Same as (a) for an enlarged lapse time window between 25 and 45 s.

为了获得台站对间经验格林函数随时间变化的定量信息,须确定各台站对的参考经验格林函数.研究表明背景噪声源的季节性变化会引起由互相关方法重构的面波信号的走时变化[28],所以应该尽量用较长时间范围(至少1年)的资料来确定参考经验格林函数.另外,汶川地震后研究区内地震活动性的增强在一定程度上降低了从背景噪声中提取的经验格林函数的信噪比和相似度,因此根据可以收集到资料的时段,本文将各台站对2007年3月至2008年5月地震前一年多较稳定的互相关函数进行叠加作为参考经验格林函数,图 2右侧的波形给出了AXI-PWU台站对的参考经验格林函数.

2.4 相对地震波速度变化的计算

假设台站对间地壳介质的相对地震波速度(Δv/v)是空间均匀变化的,则可以通过测量每天经验格林函数与参考经验格林函数的相对走时偏移(Δτ/τ)来计算相对波速的变化[19],即

(1)

因此测量走时偏移(Δτ)是本文数据处理中重要的步骤.

经验格林函数是双边时间函数(如图 2所示),在数据处理过程中我们保存了流逝时间-1000s~1000s的数据,从该时间长度中确定测量走时偏移的时间范围时,我们兼顾经验格林函数的高信噪比以及每天经验格林函数与参考经验格林函数的高相似度这两个条件,用人工判别的方式确定了测量各台站走时偏移的经验格林函数的时间范围约为1倍瑞雷波到时至3~5倍瑞雷波到时.

在计算相对走时偏移前,需要进行仪器时钟误差校正(简称钟差校正).地下介质物性的变化引起经验格林函数走时的变化应该是对称的,即经验格林函数波形的正、负半支同时压缩或拉伸,而仪器的时钟误差会使波形整体向左或向右移动[29].所以我们先计算了各台站对每天的经验格林函数与参考经验格林函数在前述计算时间范围内的整体走时偏移,再用该偏移值对每天的经验格林函数进行钟差校正.计算结果表明多数台站对每天经验格林函数的最大钟差为数十毫秒,只有少数台站对存在大于100ms的钟差,而这些较大钟差也并非发生在各台站对的地震波速度变化较大的时段.

本次研究应用移动窗口互相关谱(moving-windowcross-spectrum)的方法测量相对地震波速度变化,在文献[19]中对该方法进行了详细的介绍,因此在这里只作简要说明.在前面确定的计算时间范围内,用1s的步长取以流逝时间为中心的一系列时间窗,通过计算互相关谱的方法在0.1~0.5 Hz频带内求出时间窗内当前波形与参考波形的走时偏移(Δτ[16, 19],如果当前波形与参考波形相关系数太小,则该时刻的走时偏移计算值舍去.然后由考虑Δτ测量误差的加权线性拟合得到相对走时偏移(Δτ/τ),相对走时偏移的测量误差为

(2)

其中,

最后根据公式(1),由相对走时偏移即可求得相对波速变化(Δv/v).需要说明的是,此处的线性拟合应为强制过零点的拟合,这样可以减小仪器误差造成的走时偏移对计算结果的影响[19].图 3以AXI-PWU台站对为例示出了计算地震前、后20天的经验格林函数走时偏移和相对走时偏移的过程,可以看到用互相关谱方法计算的走时偏移值误差很小,而通过线形拟合得到的相对走时偏移的误差比计算值本身小一个数量级,所以由此得到的相对地震波速度变化是非常可靠的.

图 3 以AX-PWU台站对为例说明相对走时偏移的计算过程 (a)AXI-PWU台站对2008年4月21日(地震前20天)的经验格林函数(红色波形)及参考经验格林函数(蓝色波形);(b)2008年4月21日的走时偏移计算,圆点为走时偏移的计算值,误差棒示出其误差范围,其中,灰色圆点表示计算走时偏移窗口内当前波形与参考波形的相关系数小于0.3的计算值,这些计算值未参与线性拟合,图中灰色虚线的斜率即为通过线性拟合得到的相对走时偏移值;(c)、(d)分别与(a)、(b)类似,只是针对2008年6月1日(地震后20天)的波形和计算过程. Fig. 3 Themeasurements of relative travel-time shifts for station pairAXI-PWU (a) The red and blue waveforms represent current (April 21, 2008) and reference EGF of station pair AXI-PWU, respectively. (b) The dots denote the measurements of travel-time shifts at the lapse times corresponding to horizontal axis. The gray dots denote the calculation valuewhich correlation coefficients of currentand referencewaveforms are smaller than 0.3 and these dots are notused in further process; the slope of the gray dash line is the relative travel -time shitt which is estimated by a linear regression to all black dots. (c) and (d) for June 1, 2008 same as (a) and (b), respectively.
2.5 经验格林函数信噪比对地震波速度变化的影响

按照上述方法我们可以计算每个台站对每天的相对地震波速度变化及其误差,但这些结果是否会受到经验格林函数信噪比的影响是需要考虑的问题,为此,我们也计算了各台站对每天经验格林函数的信噪比,图 4a给出了AXI-PWU台站对经验格林函数信噪比随时间变化的曲线.可以看到主震前经验格林函数的信噪比存在一定起伏,这可能是受到噪声源的影响,但信噪比均大于7;震后信噪比水平在短时间略有降低,但从总体上看与震前水平大体相当.图 4b为AXI-PWU台站对的相对地震波速度变化,对比图 4a可以看出,震前信噪比的变化并没有引起较大的波速变化;而震后观测到波速有明显下降的时段,经验格林函数的信噪比一般大于10.由此说明,在我们的结果中每天经验格林函数的信噪比对地震波速度变化影响不大.另外,信噪比与相对地震波速度的测量误差(图 4c)也没有明显的相关性,在其他台站对上也是类似的情况,这说明我们的观测结果主要是来自地壳介质的变化,并未引入明显的其他干扰(如噪声源变化,余震活动性等).

图 4 AXI-PWU台站对经验格林函数的信噪比与相对地震波速度变化及误差 (a)每天经验格林函数的信噪比,信噪比定义为信号窗口内振幅最大值与噪声窗口内振幅的均方根的比值[26]; (b)和(c)分别为相对地震波速度变化及其误差,图中灰色竖线为汶川主震的发生时间. Fig. 4 The SNR of daily EGF and the measurements and errors of relative seismic velocity changes for station pair AXI-PWU (a) The SNR (signal to noise ratio) of daily EGF for station pair AXI-PWU. SNR is defined as the ratio of maximum amplitude in signal window to root-mean-square amplitude in noise window. (b) and (c) reprensent the measurements and errors of relative seismic velocity changes for station pair AXI-PWU, respectively.
3 结果及解释

本次研究计算了四川区域台网130个台站对的相对地震波速度变化,为了研究汶川地震震源区的波速变化,我们主要分析震源区及周边台间距在150km范围内的41个台站对的结果,图 5是一些典型台站对2007年3月至2009年3月的相对地震波速度变化曲线,左边的台站对穿过或靠近汶川余震带(图 5(a~e)),而右边的台站对则位于余震带的外围(图 5(f~j)).结果显示,汶川Ms8.0级地震前各台站对间波速变化起伏较小,其均值约为0,标准差为0.03%~0.05%,而主震后大部分台站对的波速都发生了不同程度的变化,虽然不同台站对变化持续的时间有长有短,但我们发现在汶川地震发生后的两个月内相对地震波速度变化较为显著,所以计算了各台站对主震后两个月内的相对波速变化的平均值(图 5中各曲线图上黑粗线所示),41个台站对平均相对地震波速度变化的空间分布示如图 6中,这些结果(图 56)显示了汶川地震前后震源区及周边相对波速的时空变化特征.

图 5 典型台站对相对地震波速度变化曲线 图中左边的曲线(a~e)和右边的曲线(f~j)分别表示穿过或靠近余震带和余震带外围台站对2007年3月至2009年3月相对地震波速度变化曲线.曲线上点的颜色代表相对波速变化的误差,蓝色表示误差最小,红色表示误差最大,图的下边给出了误差色标.各曲线图的左下侧为台站对代码,台站的位置参见图 6, 灰色竖线为汶川主震发生的时间.短横线为主震后两个月内的平均波速变化. Fig. 5 The relative seismic velocity changes for several typical station pairs (a)~(e) and (f)~(j) represent the relative seismic velocity changes from March 2007 to March 2009 for station pairs whose great-circle paths connecting the stations cross or are outside the aftershock zone, respectively. The colors of points indicate the errors of measurements and the error scale is shown at the bottom. The gray vertical line denote the time of the Wenchuan earthquake. Short horizontal lines indicate the average seismic velocity change for a period of two months after the main earthquake. The names of each station pair are shown at the lower lett corner of each panel.
图 6 汶川Ms8.0级地震后两个月内平均相对地震波速度变化的空间分布图 黑色的三角形代表台站,连接台站对的线条颜色代表台站对间平均相对地震波速度变化的大小,深红色表示波速减小最大,深蓝色表示波速增加最大,图的底部给出了平均相对波速变化的色标,红色五角星代表汶川主震,绿色小点代表重新定位后的余震. Fig. 6 The spatial distribution of mean values of relative seismic velocity changes for a period of two months after the Wenchuan Ms8.0 earthquake The black triangles denote seismic stations. The colors of lines connecting two stations indicate the magnitude of average seismic velocity changes as the color bar (bottom). The red star and green dots denote the main earthquake and aftershocks, respectively.
3.1 相对地震波速度随时间变化的主要特征

图 5可以看到,跨过余震带的台站对都观测到了汶川Ms8.0级地震后相对波速的急剧下降,如:XJI-YZP、MXI-YZP、AXI-PWU及JMG-PWU台站对(图 5(a,b,d,e))主震后波速降均在0.2%以上,其中AXI-PWU台站对(图 5d)的波速降最大,达到0.4%,该台站对跨过北川的擂鼓镇至平武的平通镇,此地区是地震破裂过程反演中滑动量较大的区域[3],也是地表破裂的同震位移的高值区[25].有的台站对尽管未跨过余震带,但是构成台站对的两个台站距余震带都很近,如AXI-ZJG(图 5c),也产生了较大的波速降.地震后在余震带外围的四川盆地和川西高原的台站对也都很快出现了不同程度的相对地震波速度随时间的变化(图 5(f~j)),但两个区内的台站对比跨过余震带或离余震带较近的台站对(图 5(a~e))的波速变化都明显减小,且两个区的时变特征有所不同.四川盆地内部的一些台站对的地震波速度在主震后迅速下降,如图 5h中的AXI-JJS台站对和图 5i中的JJS-YGD台站对,但盆地边缘的台站对都没有观测到汶川地震后相对波速的急剧变化,如盆地西南缘的EMS-MDS台站对(图 5g)和盆地东北缘的BZH-ZJG台站对(图 5j).位于余震带西北侧的川西高原上,台站对间相对地震波速度随时间的变化较小,如HSH-SPA台站对(图 5f).

3.2 主震后平均相对地震波速度变化的空间分布特征

图 6展示了主震后平均相对地震波速度变化的空间分布,从该图同样看到较大的平均波速降出现在那些跨过余震带或离余震带较近的台站对上,除图 5(a~e)所示台站对外,还包括YZP-MXI、MXI-ZJG、ZJG-PWU及JMG-QCH等台站对,而位于四川盆地和川西高原的台站对的平均相对波速变化则比跨过余震带的台站对的波速变化明显减小,但相对而言,四川盆地比川西高原台站对的波速变化要明显一些,如图 6中四川盆地的YZP-JJS、JJS-XCO、AXI-XCO及XCO-ZJG等台站对都观测到了较大的波速变化.而川西高原的PWU-SPA、HSH-MXI、HSH-SPA、HSH-XJI和GZA-XJI等台站对的平均相对波速变化都小于0.1%,与震前波速变化水平相当,因此,不能确认这些台站对存在可靠的波速变化量.

3.3 震源区内相对地震波速度随时间变化的空间分段性

从相对地震波速度随时间变化曲线可以看到,跨过震源区或离震源区较近的台站对波速都显示出明显的变化,但在从西南至东北约300km的余震带内却呈现出明显的空间分段特征.大致以离安县较近的ZJG台为界,西南部的台站对在主震发生后立即达到最大波速降,如图 5a中的XJI-YZP台站对和图 5b中的MXI-YZP台站对;而东北部的台站对最大波速降出现的时间几乎都滞后于主震1~4个月,如图 5e中的JMG-PWU台站对,最大波速降在2008年的9月初出现,达到0.35%,明显高于震后两个月的平均波速降0.18%,东北部的AXI-ZJG台站对(图 5c)、JMG-QCH和PWU-QCH等台站对也分别在震后的1~4个月出现了较大的波速降.

3.4 主震后相对地震波速度的恢复

震后大多数台站对的相对地震波速度经历了约5个月的恢复期,如XJI-YZP(图 5a)和JMG-PWU(图 5e)台站对等,但也有个别台站对地震波速度恢复时间较长,如MXI-YZP(图 5b)和JJS-YGD(图 5i)台站对的地震波速度直到2009年3月仍然没有恢复到震前的水平,其变化趋势还需要对近期的观测资料收集处理后进行分析.

4 讨论

近年来应用地震波进行地壳介质速度随时间变化方面的研究迅速发展,不少学者通过重复的人工、天然源及背景噪声资料都观测到了地震前后波速的明显变化,但对于这种变化的影响因素目前都还不太确定,归纳起来提出了以下几种可能的物理机制:地震中断层区的破坏与震后愈合[30, 31],强地面运动引起的地表破坏[21, 22, 32~34],断层区及周边地壳介质的应力变化[21~23, 35],地下水位的影响[18]等.

本文计算相对波速变化过程中同时应用了面波和尾波信号,而尾波一般被认为是多次散射的结果[36],从图 3可以看到,由尾波得到的相对走时变化的线性趋势与由面波得到的结果基本一致,这意味着二者可能反映大体相同空间位置的介质属性变化.由于所用地震波传播路径的复杂性,它们的影响范围是台站对周围一定的三维空间,要准确地确定我们观测到的地震波速度变化的物理机制是非常困难的,但可以结合本文的结果和已有的地质与地球物理观测排除一些影响因素,并进一步分析引起这种变化的可能机制.

(1)强地面运动引起的地表破坏

本文应用的经验格林函数频带为0.1~0.5Hz,此频带范围的面波对深度为约1~10km内的地壳剪切波速度变化较为敏感,故地表介质的变化可能对本次研究得到的波速变化结果贡献不大.另外从观测结果上看,在主震的发震断层两侧且与断层距离大体相当的川西高原和四川盆地台站对,其附近的强地面运动的峰值加速度相近[37],但是震后地震波速度变化却存在明显差异,例如位于川西高原的HSH-XJI台站对(图 5f)与位于四川盆地JJS-YGD台站对(图 5i),前者在震后波速没有明显变化,而后者却产生了较大的波速降.这也说明了浅部地壳介质在大地震的强地面运动中的破坏不是影响波速变化的主要因素.

(2)汶川主震的静态应力变化

为了探究地震波速度变化是否由地震的静态应力变化引起的,我们采用Coulomb3.0软件[38]计算了汶川地震的同震体应变.计算中采用了王卫民等[3]的有限地震断层模型,设泊松比为0.25,剪切模量为32GPa,用Okada[39]给出的均匀半空间的解析表达式得到了5km深度的体应变(图 7).从图中可以看出,体应变在余震带外的四川盆地和川西高原均呈膨胀趋势,而在余震带两端以压缩为主.震后平均相对波速变化的空间分布(图 6)表明在以张性应变为主的四川盆地一侧地震波速度下降显著,而川西高原未观测到明显的波速下降,特别是在余震带两端的地震波速度变化与压性应变不一致.由此可见主震的静态应力变化不能很好地解释本文得到的地震波速度变化.

图 7 汶川Ms8.0级地震的同震体应变 体应变的计算深度为5 km,红色区域为体应变膨胀区,蓝色区域为体应变压缩区,白色网格为断层模型的地表投影[3]. Fig. 7 The volumetric strain changes caused by the Ms8.0 Wenchuan earthquake The volumetric strain changes are calculated at the depth of 5 km. The red and blue colors indicate the region of dilatation and compression, respectively. The white grids show the surface projection of fault mode[3].

(3)断层区的破坏与震后地震波速度的恢复

反演得到的震源破裂过程显示汶川主震破裂长度达300km[3],而实地调查的地表破裂带没有这么长,主要是龙门山断裂带的东北段(见图 1)并未发现地表破裂[25].从图 6可以看出,不仅跨过地表破裂带的台站对有明显的波速变化,如MXI-YZP、MXI-ZJG、AXI-PWU和PWU-ZJG等台站对,在未产生地表破裂的龙门山断裂带东北段的台站对(如图 6中的JMG-PWU、JMG-QCH台站对)也观测到了明显的波速降低,这表明本文得到的波速变化对断层深部的性质变化更为敏感,而受地表破裂影响不大.通过对台站对所在位置及其波速变化的进一步分析,我们认为主震后尽管茂县-汶川断裂(F1)与青川断裂(F4)没有产生新的地表破裂,但是穿过这两个断裂的MXI-PWU台站对在汶川主震后有较大的波速降,这可能意味着主震使这两个断层内岩石遭受破坏、结构出现弱化或者断层及周边的应力状态发生了改变;同样,在四川盆地内的许多台站对都存在波速降,如AXI-JJS,JJS-YZP、JJS-YGD和JJS-XCO等台站对(图 6),这些应为龙门山山前隐伏断裂、成都盆地内部的新津-成都-德阳隐伏断裂(F5)和龙泉山断裂(F6)等在地震中物理属性和应力状态发生变化的反映;而川西高原内仅有的少数台站对如:HSH-SPA和HSH-XJI等台站对都没有发现明显的波速变化,可能说明在川西高原一侧的断层,如马尔康断裂(F7)和岷江断裂(F8),在台站对间的面波灵敏度核范围内[40]断层性质没有明显改变.

我们注意到主震后大多数台站对的相对地震波速度恢复到震前水平一般需要5个月左右时间,且余震带东北段与西南段的波速变化形态存在差异,我们还发现这种断层区不同空间位置的地震波速度变化的时序特性与地震序列的时空分布特征有着较好的对应关系(图 8).余震带西南部的XJI-YZP台站对(图 8a)在主震后波速急剧下降,而此时段也正是该台站对间面波灵敏度核范围内频繁的强余震活动时段,随着余震活动的减少,地震波速度在主震后迅速回升.而余震带东北部的JMG-PWU台站对(图 8b)因距主震破裂的两个最大滑动量区域都较远[3],主震造成的波速降幅度不大,但主震后该台站对间面波灵敏度核范围内频繁的强余震活动一直持续到2008年9月,地震波速度下降也相应地持续到这个时段,之后波速才开始回升,两者的变化特征非常同步.这一对应关系表明地震波速度变化与强余震的时空分布有关,可能是强余震的动态和静态应力造成断层区及周边介质物理属性的改变影响了地震波速度.

图 8 相对地震波速度变化与强余震(M≥4)时序分布 (a)上部为XJI-YZP台站对相对地震波速度变化曲线,同图 5a; 下部为该台站对间面波灵敏度核范围内M≥4余震时序分布.(b)上部为JMG-PWU台站对相对地震波速度变化曲线,同图 5e; 下部为该台站对间面波灵敏度核范围内M≥4余震时序分布.(c)中圆圈表示汶川地震M≥4余震(从2008年5月12日到2008年11月20日),两个粗虚线椭圆分别表示XJI-YZP和JMG-PWU台站间0.1~0.5 Hz频带下的面波灵敏度核范围[40],其中的黑点为图(a)、(b)中用到的余震分布. Fig. 8 The relative seismic velocity changes and M-t maps of the Wenchuan earthquake sequence with M≥4 (a) The relative seismic velocity changes for station pair XJI-YZP (top) and M-t maps of earthquakes in the sensitive kernels of surface wave (bottom). (b) Same as (b), but for station pair JMG-PWU. (c) The circles denote the aftershocks with M≥4 (from May 12, 2008 to November 20, 2008). Twodash ellipses represent the sensitive kernels of surfacewave[40] in frequency 0.1~0.5 Hz for stationpairs XJI-YZP and JMG-PWU, respectively. The black dots in ellipses denote the earthquakes plotted in (a) and (b).
5 结论

本文利用四川数字地震台网的连续波形资料,研究了汶川地震震源区及周边41个台站对2007年3月至2009年3月相对地震波速度变化.结果表明,2008年5月12日汶川Ms8.0级地震造成震源区及周边地壳介质速度的急剧下降,跨过余震带或离余震带较近的台站对波速降低较为明显,余震带外围的四川盆地内各台站对间波速也出现了不同程度的下降,但川西高原上的台站对没有观测到显著的波速变化.在从西南至东北约300km的余震带内波速变化呈现空间分段特征,大致以ZJG台为界,西南部地区主震后很快达到最大波速降,而东北部地区的最大波速降出现在主震后的1~4个月之间,相对地震波速度变化的这种分段特性与余震的时序特征有很好的同步关系.

通过背景噪声互相关方法得到的地震波速度变化反映了台站对间较大范围内地壳介质状态及属性的变化,这种变化的准确空间位置及物理机制目前还不太清楚,我们结合其他地质与地球物理资料对本研究的观测结果进行分析,认为强地面运动引起的地表破坏和地震的静态应力变化都不能很好地解释本文得到的地震波速度变化,这种变化可能反映了一定深度范围的断层区内部结构破坏和断层周围应力变化.

致谢

本次研究所用连续波形资料由中国地震台网中心和四川省地震局监测预报研究所提供.中国地震局地震预测研究所郑斯华研究员和四川地震局闻学泽研究员分别在数据处理和地质解释方面提供了帮助,在与美国罗德岛大学YangShen教授、中国地震局地震预测研究所杜建国教授的讨论中得到了有益的启示,三位匿名审稿专家提出了宝贵的意见和建议,在此一并致谢.

参考文献
[1] 张培震, 徐锡伟, 闻学泽, 等. 2008年汶川8.0级地震发震断裂的滑动速率、复发周期和构造成因. 地球物理学报 , 2008, 51(4): 1066–1073. Zhang P Z, Xu X W, Wen X Z, et al. Slip rates and recurrence intervals of the Longmen Shan active fault zone, and tectonic implications for the mechanism of the May 12 Wenchuan earthquake, 2008, Sichuan, China. Chinese J.Geophys. (in Chinese) , 2008, 51(4): 1066-1073.
[2] Royden L H, Burchiel B C, Van der Hilst R D, et al. The geological evolution of the Tibetan Plateau. Science , 2008, 321: 1054-1058. DOI:10.1126/science.1155371
[3] 王卫民, 赵连锋, 李娟, 等. 四川汶川8.0级地震震源过程. 地球物理学报 , 2008, 51(5): 1403–1410. Wang W M, Zhao L F, Li J, et al. Rupture process of the Ms8.0 Wenchuan earthquake of Sichuan. China.Chinese J.GeoPhys. (in Chinese) , 2008, 51(5): 1403-1410.
[4] 张勇, 许力生, 陈运泰. 2008年汶川大地震震源机制的时空变化. 地球物理学报 , 2009, 52(2): 379–389. Zhang Y, Xu L S, Chen Y T. Spatio-temperal variation of the source mechanism of the 2008 great Wenchuan earthquake. Chinese J.Geophys. (in Chinese) , 2009, 52(2): 379-389. DOI:10.1002/cjg2.v52.2
[5] 吴建平, 黄媛, 张天中, 等. 汶川Ms8.0级地震余震分布及周边区域P波三维速度结构研究. 地球物理学报 , 2009, 52(2): 320–328. Wu J P, Huang Y, Zhang T Z, et al. Aftershock distribution of the Ms8.0 Wenchuan earthquake and three dimension P-wave velocity structure in and around source region. Chinese J.Geophys (in Chinese) , 2009, 52(2): 320-328.
[6] Wang Z, Fukao Y, Pei S P. Structural control of rupturing of the Mw 7.92008 Wenchuan Earthquake. China.Earth Planet.Sci.Lett. , 2009, 279: 131-138. DOI:10.1016/j.epsl.2008.12.038
[7] "中国地壳运动观测网络"项目组. GPS测定的2008年汶川Ms8.0级地震的同震位移场. 中国科学(D辑) , 2008, 38(10): 1195–1206. The Project of Crustal Movement Observation Network of China. The coseismic displacement fields of Wenchuan Ms8.0earthquake occurrence in 2008 using GPS data. Science in China (Series D) (in Chinese) , 2008, 38(10): 1195-1206.
[8] 江在森, 方颖, 武艳强, 等. 汶川8.0级地震前区域地壳运动与变形动态过程. 地球物理学报 , 2009, 52(2): 505–518. Jiang Z X, Fang Y, Wu Y Q, et al. The dynamic process of regional crustal movement and deformation before Wenchuan Ms8.0 earthquake. Chinese J.Geophys. (in Chinese) , 2009, 52(2): 505-518.
[9] Mooney W D, Ginzburg A. Seismic measurements of the intemal properties of fault zones. Pure Appl.Geophys. , 1986, 124: 141-157. DOI:10.1007/BF00875723
[10] Scholz C H. The Mechanics of Earthquakes and Faulting. New York: Cambridge University Press, 2002 .
[11] Kanamori H. Mechanics of earthquakes. Ann.Rev.Earth Planet Sci. , 1994, 22: 207-237. DOI:10.1146/annurev.ea.22.050194.001231
[12] Semenov A N. Variation in the travel time of transverse and longitudinal waves before violent earthquakes. Bull.Acad.Sci.USSR, Phys.Solid Earth , 1969, 3: 245-248.
[13] Crampin S, Booth D C, Evans R, et al. Changes in shear wave splitting at Anza near the time of the North Palm Springs earthquake. J.Geophys.Res. , 1990, 95(B7): 11197-11212. DOI:10.1029/JB095iB07p11197
[14] Jin A, Aki K. Temporal change in coda Q before the Tangshan earthquake of 1976 and the Haicheng earthquake of 1975. J.Geophys.Res. , 1986, 91(B1): 665-673. DOI:10.1029/JB091iB01p00665
[15] Reasenberg P, Aki K. A precise, continuous measurement of seismic velocity for monitoring in situ stress. J.Geophys.Res. , 1974, 79(2): 399-406. DOI:10.1029/JB079i002p00399
[16] Poupinet G, Ellsworth W L, Frechet J. Monitoring velocity variations in the crust using earthquake doublets:An application to the Calaveras Fault, California. J.Geophys.Res. , 1984, 89(B7): 5719-5731. DOI:10.1029/JB089iB07p05719
[17] 丁志峰, 武岩, 王辉, 等. 2008年汶川地震震源区横波分裂的变化特征. 中国科学(D辑) , 2008, 38(12): 1600–1604. Ding Z F, Wu Y, Wang H, et al. Variations of shear wave splitting in the 2008 Wenchuan earthquake region. Science in China (Series D) (in Chinese) , 2008, 38(12): 1600-1604.
[18] Sens-Schonfelder C, Wegler U. Passive image interferometry and seasonal variations of seismic velocities at Merapi volcano, Indonesia. Geophys.Res.Lett. , 2006, 33(L21302).
[19] Brenguier F, Shapiro N M, Campillo M, et al. Towards forecasting volcanic eruptions using seismic noise. Nature Geoscience , 2008, 1: 126-130. DOI:10.1038/ngeo104
[20] Wegler U, Sens-Schdnfelder C. Fault zone monitoring with Passive Image Interferometry. Geophys.J.Int. , 2007, 168: 1029-1033. DOI:10.1111/gji.2007.168.issue-3
[21] Wegler U, Nakahara H, Sens-Schonfelder C, et al. Sudden drop of seismic velocity after the 2004 M.6.6 mid-Niigata earthquake, Japan, observed with Passive Image Interferometry. J.Geophys.Res. , 2009, 114(B06305).
[22] Brenguier F, Campillo M, Hadziioannou C, et al. Postseismic relaxation along the San Andreas fault at Parkfield from continuous seismological observations. Science , 2009, 321: 1478-1481.
[23] Xu Z J, Song X D. Temporal changes of surface wave velocity associated with major Sumatra earthquakes from ambient noise correlation. PNAS , 2009, 106(34): 14207-14212. DOI:10.1073/pnas.0901164106
[24] 吕坚, 苏金蓉, 靳玉科, 等. 汶川8.0级地震序列重新定位及其发震构造初探. 地震地质 , 2008, 30(4): 917–925. Lu J, Su J R, Jin Y K, et al. Discussion on relocation and seismo-tectonics of the Ms8.0 Wenchuan earthquake sequences. Seismology and Geology (in Chinese) , 2008, 30(4): 917-925.
[25] 徐锡伟, 闻学泽, 叶建青, 等. 汶川Ms8.0地震地表破裂带及其发震构造. 地震地质 , 2008, 30(3): 597–629. Xu X W, Wen X Z, Ye J Q, et al. The Ms8.0 Wenchuan earthquake surface ruptures and its seismogenic structure. Seismology and Geology (in Chinese) , 2008, 30(3): 597-629.
[26] Bensen G D, Ritzwoller M H, barmin M P, et al. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements. Geophys.J.Int. , 2007, 169: 1239-1260. DOI:10.1111/gji.2007.169.issue-3
[27] Shapiro N M, Campillo M, Stehly L, et al. High resolution surface wave tomography from ambient seismic noise. Science , 2005, 307: 1615-1618. DOI:10.1126/science.1108339
[28] Stehly L, Campillo M, Shapiro N M. A study of the seismic noise from its long-range correlation properties. J.Geophys.Res. , 2006, 111(B10306).
[29] Stehly L, Campillo M, Shapiro N M. Traveltime measurements from noise correlations:stability and detection of instrumental time-shifts. Geophys.J.Int. , 2007, 171: 223-230. DOI:10.1111/gji.2007.171.issue-1
[30] Li Y G, Vidale J E, Aki K, et al. Evidence of shallow fault zone strengthening after the 1992 M7.5 Landers, California, earthquake. Science , 1998, 279: 217-219. DOI:10.1126/science.279.5348.217
[31] Vidale J E, Li Y G. Damage to the shallow Landers fault from the nearby Hector Mine earthquake. Nature , 2003, 421: 524-526. DOI:10.1038/nature01354
[32] Rubinstein J L, Beroza G C. Evidence for widespread nonlinear strong ground motion in the Mw6.9 Loma Prieta earthquake. Bull.Seismol.Soc.Am. , 2004, 94(5): 1595-1608. DOI:10.1785/012004009
[33] Rubinstein J L, Beroza G C. Depth constraints on nonlinear strong ground motion from the 2004 Parkfield earthquake. GeoPhys.Res.Lett. , 2005, 32(L14313).
[34] Peng Z G, Ben-Zion Y. Temporal changes of shallow seismic velocity around the karadere-duzce branch of the north anatolian fault and strong ground motion. Pure Appl.Geophys. , 2006, 163: 567-600. DOI:10.1007/s00024-005-0034-6
[35] Niu F L, Silver P G, Daley T M, et al. Preseismic velocity changes observed from active source monitoring at the Parkfield SAFOD drill site. Nature , 2008, 454: 204-208. DOI:10.1038/nature07111
[36] Snieder R, Gret A, Douma H, et al. Coda Wave Interferometry for estimating nonlinear behaviour in seismic velocity. Science , 2002, 295: 2253-2255. DOI:10.1126/science.1070015
[37] Li X J, Zhou Z H, Huang M, et al. Preliminary analysis of strong-motion recordings from the magnitude 8.0 Wenchuan, China, Earthquake of 12 May 2008. Seismol.Res.Lett. , 2008, 79(6): 844-854. DOI:10.1785/gssrl.79.6.844
[38] Lin J, Stein R S. Stress triggering in thrust and subduction earthquakes, and stress interaction between the southem San Andreas and nearby thrust and strik-slip faults. J.Geophys.Res. , 2004, 109(B02303).
[39] Okada Y. Internal deformation due to shear and tensile faults in a half-space. Bull.Seismol.Soc.Am. , 1992, 82: 1018-1040.
[40] Zhou Y, Nolet G F, Dahlen A, et al. Global upper-mantle structure from finite-frequency surface-wave tomography. J.Geophys.Res. , 2006, 111(B04304).