2. 地球空间信息技术协同创新中心, 武汉 430079;
3. 武汉大学地球空间环境与大地测量教育部重点实验室, 武汉 430079
2. Collaborative Innovation Center of Geospatial Technology, Wuhan 430079, China;
3. Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, Wuhan 430079, China
使用地震背景噪声来探测地球内部结构是地震学的一大进步,人们逐渐认识到背景噪声中包含有丰富的地下结构信息.许多学者(如Sabra et al., 2005; Shapiro et al., 2005; Yao et al., 2006)通过对两个地震台站的长周期背景噪声进行互相关运算来提取台站间的经验格林函数,从而确定台站之间地壳和上地幔的地球介质信息.目前,背景噪声互相关技术已在全球许多地区得到应用,例如,(1)通过背景噪声互相关函数提取的面波频散曲线来获取区域速度结构.如Yang等(2007)首次利用地震背景噪声对欧洲进行了面波层析成像;Bensen等(2008)通过对美国大陆的地震背景噪声层析成像研究获得了Rayleigh面波和Love面波不同周期的频散图与群速度图.(2)在叠加长周期的互相关函数作为台站间的理论格林函数的基础上,通过提取走时偏移来监测地下介质波速变化.如Brenguier等(2008)分析法属Piton de la Fournaise火山21个台站背景噪声得到的波速变化,发现火山喷发前波速均出现了约0.1%的降低,且在火山喷发期间逐渐恢复正常;Stehly等(2016)利用2008年汶川地震震中区域21个台站的背景噪声,发现主震当天有0.07%的波速下降.
处于北美板块和科科斯板块交界的墨西哥是大地震活跃区域.DeMets等(2010)研究表明科科斯板块以大约69~75 mm·a-1速度向北美板块俯冲,由于地壳压力增速加快,巨大能量骤然增加,使得墨西哥发生了许多超过8级的大地震,如1932年Jalisco MW8.1地震和1995年Colima MW8.0地震,通常这些地震都发生在横跨墨西哥中部和南部海岸的1800 km长的俯冲带区域,并具有相似的发震机制(Melgar et al., 2018).2017年9月19日—1985年墨西哥8.1级特大地震32周年纪念日,在墨西哥中部发生了一次MW7.1地震,美国地质调查局(USGS)远震解显示震中位于(-98.489°,18.550°),深度48 km.由于震中位于内陆,地震能量全部释放在包括首都墨西哥城在内的中南部城市,部分地区建筑遭受严重破坏,造成巨大的经济损失和300多人死亡.值得注意的是2017年9月7日在墨西哥沿岸近海(-93.899°, 15.022°)还发生了一次MW8.2特大地震,两次地震相隔约650 km(Segou and Parsons, 2018).Hayes(2017)研究表明MW7.1地震是一次正断层机制的地震,倾角44°~47°,震源在科科斯俯冲带板内深部,而不是浅层的板块边界,同时其余震极少,与该区域其他地震机制不太相同.
在墨西哥中部有2个宽频带固定台站,它们恰好位于2017年地震的震源区,这为运用互相关函数研究较长时间尺度内震源区介质物理性质的变化提供了较好的条件.本文首先对双台站24个月的波形资料进行不同频带下的互相关处理;然后利用小波变换对互相关函数进行滤波来提高格林函数的稳定性;在此基础上,分别使用Stretching和MWCS方法测量的走时偏移来分析该时段台站间地震波速变化;最后结合已有的发震机理、破裂过程等研究成果,探讨引起这种变化的原因.
1 数据处理本文采用美国地震学研究联合会(Incorporated Research Institutions for Seismology,IRIS)提供的宽频带垂直分量连续波形数据来计算地震波速变化.数据处理分为四步:(1)对单台数据进行预处理得到地震背景噪声;(2)通过计算台站对互相关函数并叠加来提取经验格林函数;(3)计算不同时段的经验格林函数与参考经验格林函数的走时偏移δt/t;(4)计算地震波波速变化δv/v及其误差.
1.1 观测数据从IRIS数据中心查询得知,2017年墨西哥地震震中周围2°范围内仅有2个地震仪,分别为G.UNM和MX.TLIG,其中G.UNM(-99.18°,19.33°)位于墨西哥城,MX.TLIG(-98.57°,17.56°)位于格雷罗州东部,两个台站距震中分别为125 km和85 km,两台站间相距207 km(图 1).为了分析2017年地震前后该区域的地震波速变化,本文选取了G.UNM和MX.TLIG台站从2016年4月1日至2018年4月1日20 Hz的垂直分量连续波形数据,所选时段内该区域没有其他超过MW4.5级的地震,G.UNM站在2016年6月有连续13天的数据缺失,MX.TLIG站没有大的数据缺失,因此数据基本能够反映该时段内的介质波速变化.
单台站数据预处理是为了获取高质量的地震背景噪声(Bensen et al., 2007).在预处理时,首先将连续波形数据按天分割,去除仪器响应、去均值和去线性趋势,并通过0.1~0.2 Hz,0.2~0.5 Hz,0.05~0.1 Hz等不同频段的带通滤波来提取对应深度的波速信息;随后选择One-bit方法(Cupillard et al., 2011)进行时域归一化来降低地震、仪器误差、毛刺等非平稳噪声源的影响,该方法直接将地震波形振幅正值取1,负值取-1,实现简单;最后通过谱白化(频域归一化)压制单频强烈的微震信号,这样就得到了背景噪声.
1.3 背景噪声互相关对背景噪声进行互相关处理可以得到互相关函数,对互相关函数求导可以得到经验格林函数(下文统称为格林函数).在处理过程中,将每日数据划分为长度为1800 s的移动时间窗口,然后对两个台站对应窗口内的背景噪声在频率域内进行互相关运算,保留滞后时间[-240,240] s的格林函数,并将当天所有的移动窗口格林函数叠加,就可以获得该台站对的单日格林函数(图 2a).随后,对前后各15天(合计30天)的互相关函数进行叠加,就得到了该天的互相关函数,称为当前格林函数(图 2b).最终将2016年4月1日至2018年4月1日所有的单日格林函数叠加为参考格林函数(图 2c).
由于噪声源的不均匀分布和噪声成分的不确定性,加上台站间距较远,发现得到的单日格林函数信噪比低,难以分辨波形,而通过叠加得到的当前格林函数与参考格林函数的信噪比逐渐升高,可以见到清晰的直达面波和尾波(图 2).因此可以利用参考格林函数对单日格林函数进行滤波来进一步提高信噪比,从而得到更高质量的单日格林函数.Hadziioannou等(2011)在处理帕克菲尔德地震背景噪声时使用了S变换处理互相关函数;Stehly等(2016)在处理汶川地震背景噪声时使用曲波变换处理互相关函数,结果均表明利用这些变换构建的自适应滤波能够提高单日格林函数信噪比,从而提升时间分辨率.由于小波变换具有多尺度的时频窗口特性,在地震研究中也常用小波变换分离噪声,因此在本研究中,选取了小波变换来构造滤波器,可以发现小波变换后的单日格林函数(图 2d)相比于原单日格林函数(图 2a)有了明显提升.
1.4 走时偏移计算计算走时偏移δt/t需要将当前格林函数与参考格林函数在时间窗[t1,t2]内进行对比,t1和t2一般通过区域速度结构得到直达面波或尾波的理论到时来确定,由于尾波在介质中的传播时间更长,对速度变化更加敏感,因而更受到研究者的青睐.然而本研究中的两个台站相距207 km,远大于一般研究中十几km至几十km的台间距,尾波的走时增加,会使得单日格林函数的信噪比降低,当前格林函数与参考格林函数的相关系数降低(Sabra et al., 2005).为了提高单日格林函数的信噪比,我们选择了包括直达波和尾波的[70, 190] s时间窗来计算走时偏移.
计算时间窗内走时偏移常用的方法主要有Stretching或MWCS(Moving-Window Cross-Spectral)算法.其中,Stretching方法(Hadziioannou et al., 2009)是将当前格林函数正负轴相加平均,乘以因子(1+δt/t),以使当前格林函数与参考格林函数在[t1,t2]段的相关系数最大(图 3),同时相关系数作为相关性高低的指标,也代表着噪声源分布的稳定性,相关系数越高则表示噪声源越稳定.而MWCS方法(Clarke et al., 2011)的原理是在当前格林函数上划分1 s步长的滑动窗,每个滑动窗内计算互相关谱求出两个格林函数的频域相位差得到δt,[t1, t2]内线性拟合所有滑动窗的δt得到斜率,即为当日的δt/t(图 4).
最后,根据介质沿空间均匀变化的理论假设,得出相对波速变化与相对走时偏移的关系(Ratdomopurbo and Poupinet, 1995):
(1) |
对于Stretching方法,设定最大δt/t为0.5%,并在[-0.5%,0.5%]等分1000个值,局部搜索出相关系数最大时的δt/t.同时为了分析地震波速变化的精度水平,还计算了其误差(Weaver, 2011):
(2) |
其中,C为相关系数,T为频带宽的倒数,ω为中心频率.
对于MWCS方法,对走时偏移进行加权线性拟合,并强制过零点来消除仪器误差.其误差计算公式(Clark et al., 2011)为
(3) |
其中,pi表示权重,ti为滞后时间,〈t〉为t的加权平均值.
2 结果与分析 2.1 Stretching方法结果与分析图 5给出了在滞后时间[70,190] s使用Stretching方法得到的墨西哥地震介质速度变化,结果表明在地震发生前,波速基本处于较平稳状态;2017年9月19日地震发生时有明显的波速抖动,意味着地震会引起地震波速度的突然变化,其中0.1~0.2 Hz的结果显示波速上升约0.07%,平均相关系数为0.96,误差为0.77×10-4,说明此波动变化值可信,也意味着G.UNM和MX.TLIG之间确实能监测到比较可靠的地震波速度变化.0.2~0.5 Hz与0.05~0.1 Hz的结果分别显示下降0.05%和上升0.03%,但是波动及误差较大,具有较大的不确定性.
利用小波变换提高互相关函数信噪比可以获取更精确的信息,经过对多种小波基和变换方法的尝试,最终选用db5小波做单级小波变换,变换后结果(图 5(d—f))与原波速变化趋势相似,但是其平均误差更低,整体趋势更加平稳.说明在单日格林函数中包含噪声成分不明确的条件下,利用参考格林函数构建的自适应滤波对于提高信噪比有明显效果.另外0.2~0.5 Hz的震时波速变化仍显示了微小的下降趋势,0.05~0.1 Hz的结果则几乎没有变化.
地震波速变化误差主要分布在三个区间内,分别为第70~90天、第240~320天以及第600~680天,第一个区间主要是因为2016年6月16日开始的13天数据缺失,导致该时段缺少了相应的单日格林函数,信噪比降低而造成的误差;后面两个误差区间则是由于当前格林函数受到了季节性变化的影响(Meier et al., 2010).误差成因还可以在当前格林函数(图 6)中得到印证,由于受到噪声源分布以及台站对方位角的影响,本应正负轴对称的格林函数(图 2和图 6)均表现出负轴波形大于正轴的情形,表明背景噪声主要受到太平洋微震的影响(Sabra et al., 2005),使得MX.TLIG方向微震能量多于G.UNM方向.
选择滑动窗30 s,滑距1 s,最小相关系数0.8,最大误差0.1 s和最大走时偏移0.1 s,基于MWCS方法给出的结果见图 7.对比Stretching方法和MWCS方法的结果,可以发现MWCS方法结果误差稍低,但在震前对于速度波动更加敏感,几处比较明显的尖峰不符合该区域在2017年地震前没有大于4.5级地震发生的情况,但是在0.1~0.2 Hz和0.2~0.5 Hz的结果中仍然分别显示了相同的上升和下降趋势.另外,小波变换对于MWCS方法仍具有降噪平滑的效果,但是不及对Stretching方法的效果显著.因此,与MWCS方法相比,Stretching方法更加稳定,受地震之外的噪声影响较小.由于稀疏的台站分布和数据记录,本研究更倾向于Stretching方法测量得到的结果.
本文所选用的0.1~0.2 Hz背景噪声,其对5~10 km深度的剪切波速更为敏感,震时波速表现为上升趋势,而0.2~0.5 Hz的结果则表明更浅的地层波速有下降趋势.与该区域的历史地震相比,2017年MW7.1地震属于比较特殊的板内地震,其发震断层处于俯冲带由平板转向高倾角的区域.Melgar等(2018)认为该地震是由弯曲应力积累而引发的,并给出最大滑动量约1.3 m的结果(图 8a).Pérez-Campos等(2008)利用接收函数发现地壳浅层有速度层分界,表明该区域地壳结构复杂,存在产生本研究中速度变化差异的构造背景.
由于地震波传播的复杂性,要找到台站对区域的地震波波速变化的物理机制极其困难.如Zaccarelli等(2011)在分析2009年L′Aquila地震时认为浅层地壳结构中裂缝和孔隙密度的增加导致同震波速的变化;而Hadziioannou等(2011)则将2004年Parkfield地震波速变化归因于剪切模量下降导致浅层地壳的破坏;Stehly等(2016)认为导致2008年汶川地震波速变化的原因可能是脆性上地壳的孔隙回弹.刘志坤和黄金莉(2010)总结了造成地震波速变化的几种可能物理机制:地震断层区的破坏与愈合,强地面运动的地表破坏,断层区及地壳介质的应力变化等.由于2017年地震发生在俯冲板块的上边界,属于科科斯板内事件,同时本文的计算深度所反映的是北美板块浅层地壳的波速变化,因此地震断层的影响应该不是直接原因.此次地震产生的地表最大加速度约130 Gal,而这个数值的加速度在其他地震中并不一定引起明显波速变化,同时本文两个深度的波速变化并不一致,因此强地面运动的影响应该也不是其主要因素.
同时,本研究还采用Coulomb软件(Toda et al., 2005)和有限断层模型(Melgar et al., 2018)计算了地震造成的应力应变效应(图 8b),发现体应变在震源中心有膨胀趋势,而震源外围几十公里范围则是压缩趋势,台站对穿过正应变区与负应变区引起的波速变化有正有负,且数值与其他地震相比较小,另外依据0.1~0.2 Hz和0.2~0.5 Hz频段波速变化的差异,我们也比较了对应深度的应力变化,发现变化均不大,所以同震应力应变效应也不是导致波速变化的合理解释.
与其他地震相似,2017年墨西哥地震之后的波速逐渐恢复.其他学者使用GPS监测基线变化或是沿断层走向的位移量,发现GPS位移或震颤活动与震后速波速恢复趋势有吻合性(Brenguier et al., 2008;Rivet et al., 2014).或是断层区不同空间位置的地震波速度变化的时序特性与地震序列的时空分布特征有着较好的对应关系,认为其与震后的余震或应力松弛有关(刘志坤和黄金莉, 2010).根据墨西哥国家地震服务的测定,此次地震的余震较少,在震后5天仅有15个约3级的余震,且分布稀疏,无法看出余震与波速变化之间的明显关系.
上述分析表明本次震源区的地震波波速变化与断层区运动、静态应力应变状态、强地面运动等均不存在显著关系,参考Zaccarelli等(2011)的研究,本文认为造成地震波速变化的最大可能原因是2017年墨西哥地震导致震源区内浅层地壳结构破坏,改变了岩石密度、孔隙度等物理性质,从而导致震源区地震波波速随之发生了变化.
4 结论本文利用2个地震台24个月的背景噪声资料提取了2017年MW7.1墨西哥地震震中区域的地震波速度变化.研究发现虽然较远的台站间距造成互相关函数不稳定,但是通过小波变换能够提升单日格林函数信噪比,使得波速变化结果更精确,未来可以通过对背景噪声性质的深入分析,尝试其他的滤波器,获得更佳结果.同时,通过Stretching方法和MWCS方法分别提取到的地震波速变化比较表明MWCS方法对于噪声的敏感导致结果波动较大,而Stretching方法结合自适应滤波效果更好,更能保持结果稳定性.
两种方法的结果均表明2017年地震震前的波速基本处于较平稳状态,震时5~10 km深度地震波速度升高了约0.07%,而2~5 km深度波速下降约0.05%, 震后波速逐渐恢复至正常水平.其原因可能是地震对不同深度地层产生不同的影响,破坏了脆性岩石的结构,造成了破裂或孔洞,引起了地震波速的变化,之后随着断层和孔隙的愈合,地震波速度又逐渐恢复.
大部分研究使用了较多的台站来研究整个区域的平均波速变化,偶尔会挑选典型的台站进行分析,如Zaccarelli等(2011)仅使用3个台站来研究L′Aquila地震造成的波速变化.由于台站数量和分布的限制,本研究仅提取了双台站间的波速变化,因此该结果仅代表G.UNM和MX.TLIG之间的地下介质信息,不能反应整个区域的波速变化.未来可以通过增加观测台站数量,得到分辨率更高、范围更广的地下介质信息.
此外,本文未发现明显的震前波速变化征兆,表明该时段的区域地震波速度变化主要是由地震引起的构造变化所触发,未来可以结合岩石孔隙水或应力的变化、地温信息和地球动力学理论等来研究地震的孕育和发生.
致谢 感谢两位审稿专家给出的建设性意见.研究中使用的宽频带地震连续波形资料由IRIS提供,数据处理中使用了MSNoise程序包,图件的绘制采用了GMT软件,在此一并致谢.
Bensen G D, Ritzwoller M H, Barmin M P, et al. 2007. Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurements. Geophysical Journal International, 169(3): 1239-1260. DOI:10.1111/j.1365-246X.2007.03374.x |
Bensen G D, Ritzwoller M H, Shapiro N M. 2008. Broadband ambient noise surface wave tomography across the United States. Journal of Geophysical Research:Solid Earth, 113(B5): B05306. DOI:10.1029/2007JB005248 |
Brenguier F, Campillo M, Hadziioannou C, et al. 2008. Postseismic relaxation along the San Andreas fault at Parkfield from continuous seismological observations. Science, 321(5895): 1478-1481. DOI:10.1126/science.1160943 |
Clarke D, Zaccarelli L, Shapiro N M, et al. 2011. Assessment of resolution and accuracy of the Moving Window Cross Spectral technique for monitoring crustal temporal variations using ambient seismic noise. Geophysical Journal International, 186(2): 867-882. DOI:10.1111/j.1365-246X.2011.05074.x |
Cupillard P, Stehly L, Romanowicz B. 2011. The one-bit noise correlation:a theory based on the concepts of coherent and incoherent noise. Geophysical Journal International, 184(3): 1397-1414. DOI:10.1111/j.1365-246X.2010.04923.x |
DeMets C, Gordon R G, Argus D F. 2010. Geologically current plate motions. Geophysical Journal International, 181(1): 1-80. DOI:10.1111/j.1365-246X.2009.04491.x |
Hadziioannou C, Larose E, Coutant O, et al. 2009. Stability of monitoring weak changes in multiply scattering media with ambient noise correlation:laboratory experiments. The Journal of the Acoustical Society of America, 125(6): 3688-3695. DOI:10.1121/1.3125345 |
Hadziioannou C, Larose E, Baig A, et al. 2011. Improving temporal resolution in ambient noise monitoring of seismic wave speed. Journal of Geophysical Research:Solid Earth, 116(B7): B07304. DOI:10.1029/2011JB008200 |
Hayes G P, Wald D J, Johnson R L. 2012. Slab1.0:A three-dimensional model of global subduction zone geometries. Journal of Geophysical Research Solid Earth, 117(B1). DOI:10.1029/2011JB008524 |
Hayes G P. 2017. The finite, kinematic rupture properties of great-sized earthquakes since 1990. Earth and Planetary Science Letters, 468: 94-100. DOI:10.1016/j.epsl.2017.04.003 |
Liu Z K, Huang J L. 2010. Temporal changes of seismic velocity around the Wenchuan earthquake fault zone from ambient seismic noise correlation. Chinese Journal of Geophysics (in Chinese), 53(4): 853-863. DOI:10.3969/j.issn.0001-5733.2010.04.010 |
Meier U, Shapiro N M, Brenguier F. 2010. Detecting seasonal variations in seismic velocities within Los Angeles basin from correlations of ambient seismic noise. Geophysical Journal of the Royal Astronomical Society, 181(2): 985-996. |
Melgar D, Pérez-Campos X, Ramirez-Guzman L, et al. 2018. Bend faulting at the edge of a flat slab:The 2017 MW7.1 Puebla-Morelos, Mexico Earthquake. Geophysical Research Letters, 45(6): 2633-2641. DOI:10.1002/2017GL076895 |
Pérez-Campos X, Kim Y H, Husker A, et al. 2008. Horizontal subduction and truncation of the Cocos Plate beneath central Mexico. Geophysical Research Letters, 35(18): L18303. DOI:10.1029/2008GL035127 |
Ratdomopurbo A, Poupinet G. 1995. Monitoring a temporal change of seismic velocity in a volcano:Application to the 1992 eruption of Mt. Merapi (Indonesia). Geophysical Research Letters, 22(7): 775-778. DOI:10.1029/95GL00302 |
Rivet D, Campillo M, Radiguet M, et al. 2014. Seismic velocity changes, strain rate and non-volcanic tremors during the 2009-2010 slow slip event in Guerrero, Mexico. Geophysical Journal International, 196(1): 447-460. DOI:10.1093/gji/ggt374 |
Sabra K G, Roux P, Kuperman W A. 2005. Emergence rate of the time-domain Green's function from the ambient noise cross-correlation function. The Journal of the Acoustical Society of America, 118(6): 3524-3531. DOI:10.1121/1.2109059 |
Segou M, Parsons T. 2018. Testing earthquake links in Mexico from 1978 to the 2017 M=8.1 Chiapas and M=7.1 Puebla shocks. Geophysical Research Letters, 45(2): 708-714. DOI:10.1002/2017gl076237 |
Shapiro N M, Campillo M, Stehly L, et al. 2005. High-resolution surface-wave tomography from ambient seismic noise. Science, 307(5715): 1615-1618. DOI:10.1126/science.1108339 |
Stehly L, Froment B, Campillo M, et al. 2016. Monitoring seismic wave velocity changes associated with the MW7.9 Wenchuan earthquake:increasing the temporal resolution using curvelet filters. Geophysical Journal International, 201(3): 1939-1949. |
Toda S, Stein R S, Richards-Dinger K, et al. 2005. Forecasting the evolution of seismicity in southern California:Animations built on earthquake stress transfer. Journal of Geophysical Research, 110(B05S16): 361-368. |
Weaver R L. 2011. On the amplitudes of correlations and the inference of attenuations, specific intensities and site factors from ambient noise. Comptes Rendus Géoscience, 343(8-9): 0-622. |
Yang Y, Ritzwoller M H, Levshin A L, et al. 2007. Ambient noise Rayleigh wave tomography across Europe. Geophysical Journal International, 168(1): 259-274. DOI:10.1111/j.1365-246X.2006.03203.x |
Yao H, Hilst R D V D, Hoop M V D. 2006. Surface-wave array tomography in SE Tibet from ambient seismic noise and two-station analysis-I. Phase velocity maps. Geophysical Journal International, 166(2): 732-744. DOI:10.1111/j.1365-246X.2006.03028.x |
Zaccarelli L, Shapiro N M, Faenza L, et al. 2011. Variations of crustal elastic properties during the 2009 L'aquila earthquake inferred from cross-correlations of ambient seismic noise. Geophysical Research Letters, 38(24): L24304. DOI:10.1029/2011GL049750 |
刘志坤, 黄金莉. 2010. 利用背景噪声互相关研究汶川地震震源区地震波速度变化. 地球物理学报, 53(4): 853-863. DOI:10.3969/j.issn.0001-5733.2010.04.010 |