运用亮度函数SSA扫描连续波形资料方法确定玉树7.1级地震前辐射源*
2014-12-25沙成宁张晓清袁伏全
沙成宁,张晓清,袁伏全
(青海省地震局,青海西宁810001)
0 引言
开展强震孕震环境研究,揭示地下介质受力作用的过程,分析强震前震源体的空间分布情况,了解震源体周围受力状态对准确判断强震震情具有重要的意义。过去10多年来,青海省及其周边先后发生了1997年11月8日玛尼7.5级、2001年11月14日昆仑山口西8.1级、2008年3月21日于田7.3级、2008年5月12日汶川8.0级、2010年4月14日玉树7.1级、2014年2月12日于田7.3级地震。针对这些强震孕震环境及前兆特征,很多学者运用不同方法和不同的数据进行研究分析,开展了大量的卓有成效的研究工作,从不同的侧重点给出强震孕育发生特征。祝意青等(2008)运用重力观测资料对汶川8.0级地震孕震环境进行了分析研究,结果表明其孕震区域位于重力正负分区的梯度带内。程佳等 (2011)运用应力触发方法研究了玛尼、昆仑山口西、汶川地震所形成的同震和震后形变场的变化过程与特征,研究结果表明先后发生的强震是有关联的。陈祖安等 (2011)运用 (DDA+FEM)理论和形变资料,研究了2001年昆仑山8.1级地震对2008年汶川8.0级地震孕育产生的影响。郑现等 (2011)运用宽频带数字地震台网资料,开展中国大陆中东部地区背景噪声的瑞利波群速度层析成像,给出了该区域地壳厚度的差异、空间分布以及中强地震的空间分布特征。
在数字地震波形资料分析处理方面,一些学者开展了震源矩张量反演、震源破裂过程等方面的研究。近年来,中长周期 (20~40 s)资料也逐步得到开发利用,李文军等 (2006)使用SSA(Source Scanning Algorithm)方法进行地震定位,研究结果表明该方法在不用精确读取到时和计算理论地震图的情况下达到比较理想的定位效果。丁宁霞和张晓清 (2013)运用SSA方法结合青海区域地震台网观测资料确定2008年11月10日大柴旦MS6.3地震震源破裂面。杨海波 (2010)开展了基于SSA算法的震源破裂过程反演方法系统整理研究。陈长云和任金卫 (2013)研究巴颜喀拉块体东部活动块体的划分、形变特征及构造,从构造本和形变的方面给出了巴颜喀拉块体的东向运动存在自西向东的差异性。邓津 (2012)通过对低频事件的全球分布及汶川震前实例分析,其研究得出低频滑移事件在各大板块周围普遍存在。张晓清 (2000)在青海东部及其邻近地区的地震层析成像工作中给出了研究区域的P波速度结构。杜海林 (2007)将聚束能量分析方法和迁移叠加方法应用于宽频带台阵,分析了2004年12月26日苏门答腊—安达曼地震的能量辐射源,得到了震源破裂长度、方向、持续时间和破裂速度。本文在上述工作的基础上,使用青海及周边省份地震台网记录的连续数字地震波形资料中的中长周期成份进行亮度函数扫描辐射源研究,探索SSA方法确定中长周期辐射源的空间位置,为确定强震孕震区提供技术支持。
1 SSA方法和原理
亮度函数SSA扫描算法年由 Kao和 Shan(2004)在北Cascadia俯冲带的慢地震研究中提出,他们通过SSA方法,根据相应的“亮度”函数扫描,找出慢地震的震源位置和发生时间。SSA扫描算法是在整个时空寻找可能辐射源的方法,其基本原理如图1所示。图中,空间上的点η在τ时刻的“亮度”可以由相应台站的理论到时 (即τ加上相应的走时 taη,tbη,tcη)计算而来,如果某处“亮度”大 (图中空心五角星)表明这里有事件发生,如果“亮度”小 (实心五角星)表示这里没有事件发生。
图1 SSA扫描法原理示意图Fig.1 Sketch map of principle of SSA scanning algorithm
假设一个事件由N个台站记录到 (图1中的A、B、C点),首先分别归一化每个数字地震台站的记录un,然后计算某个点 (η)某个时刻 (τ)的“亮度”函数为
式中,un为归一化地震记录,tηn为从点η处到台站n计算的某个最大震相的走时,整个时空中“亮度”函数取最大值时是震源发生的位置和事件发生的时刻。由于使用区域台网的资料,对速度模型的精度要求较高,目前对速度模型认识程度还达不到要求,造成走时计算有一些系统偏差,可以通过平滑的方法来预先处理记录,得到每个台站光滑的“亮度”函数序列为
smooth表示一种平滑方法。对于某一空间点 (η)某个时刻 (τ)的“亮度”函数可以表示为
2 选取资料处理
使用青海及其周边的甘肃、新疆、西藏数字地震台网共65个台站的数据,青海台网地震台站空间分布如图2所示。选取2010年4月14日玉树7.1级地震前,2010年1月1日至4月11日青海省数字地震台网记录的信噪比较高的连续波形,并挑选出每天0~2时连续无地震记录波形、无高频干扰、无突跳点,仪器工作正常的垂直分量波形进行研究。并对地震记录仪幅频特性进行了核实,所使用资料的地震记录仪大部分为BBVS-60和JCZ-1,带宽分别为40 Hz~60 s和50 Hz~360 s,幅频特性曲线如图3所示。
图2 青海台网地震台站空间分布Fig.2 Spatial distribution of stations in Qinghai Seismic Network
图3 大柴旦BBVS-60(a)和格尔木JZC-1(b)仪器幅频特性曲线Fig.3 Amplitude-frequency characteristic of BBVS-60 instrument at Dachaidan Station(a)and JZC-1 instrument at Golmud Station(b)
图4 垂直分量连续波形资料FFT频谱Fig.4 FFT frequency spectrum of continuous seismic waveform in vertical component
3 结果及分析
运用亮度函数SSA扫描连续波形资料方法,确定玉树7.1级地震前30 s辐射源的空间位置变化和随时间推移的演化过程。
首先,通过对所选取资料做FFT频谱分析发现 (图4),在无地震发生时,大部分地震台站连续波形资料均存在0.03 Hz附近的优势频率。再对连续波形进行滤波,将其中的对应信号提取出来,并通过亮度函数SSA进行扫描计算,确定相应周期的辐射源的空间分布。由于SSA方法是在所给定时空范围内,以亮度函数达到最大作为判据确定辐射源的空间位置,因此,只需要给定所要扫描的的空间范围 (经度、纬度、深度)、相应频段波的波速值、地震台站的空间坐标以及空间扫描的步长,即可由亮度函数的最大值确定辐射源的空间位置。本文给出的初始空间范围为 (88.5°~108°E,30°~50°N),地壳深度为0 ~70 km。参考余大新等 (2014)利用Rayleigh相速度和群速度联合反演青藏高原东北缘速度结构,和郑现等(2012)给出的中国大陆中东部地区基于背景噪声的瑞利波层析成像等研究结果,根据现有计算设备的计算能力,为了确保研究区内计算结果的空间分辨率达到最好,确定本文所采用的波速值范围2.8~3.2 km/s;计算网格步长水平向为0.1°,垂直向为5 km。
运用Butterworth滤波器提取青海台网连续波形中频率为0.2~0.4 Hz的波形资料并进行滤波,如图5所示。对提取的资料进行亮度函数扫描,亮度函数值的空间分布如图6所示。从图中可以看出,2010年1月10日亮度函数最高的区域分布在柴达木盆地的西部及西南部,位于东昆仑断裂带的西段 (图6a),最高亮度值达到0.74,深度15 km;2月10日和3月10日亮度函数高值区域的主体仍然在东昆仑断裂带的西段,明显向东昆仑断裂带以东扩展,亮度函数最高值分别为0.69和0.68(图6b、c),深度分别为15 km和17 km;4月10日 (玉树地震前4天)亮度函数高值区移至东昆仑断裂带的东段,并且明显向南扩散,亮度函数最高值达到0.81,深度为19 km;4月14日玉树7.1级地震发生在亮度函数高值区附近 (图6d)。
图5 垂直分量连续波形和滤波后的波形Fig.5 Continuous waveform and waveform after filtering in vertical component
本文计算的亮度函数结果的时空分布特征呈现出最高亮度函数的中心位置基本处于东昆仑断裂带内,在玉树7.1级地震前沿东昆仑断裂逐步由西向东移的现象。陈祖安等 (2011)对青藏高原及其东侧四川盆地鄂尔多斯块体地区构造块体的研究表明,2001年11月14日昆仑山口西8.1级地震,最大水平错距4.5 m,地表最大位移为6 m左右;杨国华等 (2008)根据GPS观测资料,震源区附近的青藏高原向北运动的量级在南端超过了20 mm/a,面应变以压性为主的方式积累应变能,最大绝对值位于玉树地震和汶川地震震中的周围地区,上述研究结果表明在昆仑山口西8.1级地震之后,玉树7.1级地震之前东昆仑断裂带东部及其周围区域存在较高水平应变,从构造和形变方面给出孕震证据。此外,从岩石力学实验角度,根据邓志辉等 (1995)孤立型大事件常发生于高能量向低能量的突变带或高能量背景区内的相对低能量区,高能量不是发震的充分条件,能量空间分布的差异也是发震的重要条件,失稳可能从弱的地方开始;马瑾等 (2008)对断层变形过程中热场与应变场的关系岩石力学实验研究结果表明,岩石断层破裂前存在热红外辐射的亮度温度场和温度场的变化。给出的亮度函数结果的时空分布特征为,空间上表现为最高亮度函数值与玉树7.1级地震震中不重合,玉树地震震中位于辐射源中心的南部,亮度函数最高值的边缘;时间上表现为随时间推移辐射源逐步向震中移动。本文给出的结果与近几年对2008年5月12日汶川8.0级地震、2010年4月14日玉树7.1级地震及昆仑山口西8.1级地震的研究结果相一致,玉树地震之前确实存在30 s左右的辐射源。
图6 亮度函数值的空间分布Fig.6 Spatial distribution of brightness function value
4 结论及讨论
本文使用SSA扫描方法计算了2010年4月14日玉树7.1级地震前青海区域台网连续波形资料,发现在青海东昆仑断裂带及其附近存在30 s周期的辐射源,辐射源随空间和时间的变化对未来的中强地震发震地点有一定指示意义,期望为今后的孕震区的判定提供一种定量判定方法。
由于本文仅对辐射源扫描开展了二维扫描计算,辐射源的大小与地震震级的关系、辐射源时空特征随辐射周期的变化、未来中强地震的发震地点位于辐射源的方位等问题有待于做更进一步的研究。此外,进行频谱分析时,注意对资料的挑选,尽量不使用有地震波形记录的资料,减少提取相应频率的难度。
陈长云,任金卫.2013.巴颜喀拉块体东部活动块体的划分、形变特征及构造意义[J].地球物理学报,56(12):4126-4140.
陈祖安,林邦慧,白武明,等.2011.2001年8.1级昆仑山大震破裂过程及对2008年汶川8.0级大震孕育发生影响的研究[J].地球物理学报,54(1):108-119.
程佳,刘杰,甘卫军,等.2011.1997年以来巴颜喀拉块体周缘强震之间的黏弹性触发研究[J].地球物理学报,54(8):1997-2010.
邓津.2012.低频事件的全球分布及汶川震前实例分析[D].北京:中国地震局地球物理研究所.
邓志辉,马胜利,马瑾,等.1995.粘滑失稳及其物理场时空分布的实验研究[J].地震地质,17(4):305-315.
丁宁霞,张晓清.2013.利用震源扫描方法计算2008年大柴旦6.3级地震震源破裂过程[J].高原地震,25(4):26-30.
杜海林.2007.2004年苏门答腊_安达曼大地震能量辐射源的时间域台阵技术分析[D].北京:中国地震局地球物理研究所.
李文军,李丽,陈棋福.2008.用震源扫描算法(SSA)研究列车源的运动[J].地球物理学报,51(4):1147-1151.
马瑾,马少鹏,刘培洵,等.2008.识别断层活动和失稳的热场标志[J].地震地质,30(2):364-382.
杨国华,张晓东,沈舞春,等.2008.昆仑山口西8.1级地震震后中国西部地壳水平位移场的变化特征[J].地震研究,31(1):77-82.
杨海波.2010.基于震源扫描算法的震源破裂过程反演方法研究[D].成都:成都理工大学.
余大新,李永华,吴庆举,等.2014.利用Rayleigh波相速度和群速度联合反演青藏高原东北缘S波速度结构[J].地球物理学报,57(3):800-811.
张晓清.2000.青海东部及其邻近地区的地震层析成像[D].合肥:中国科学技术大学.
郑现,赵翠萍,周连庆,等.2012.中国大陆中东部地区基于背景噪声的瑞利波层析成像[J].地球物理学报,55(6):1919-1928.
祝意青,梁伟锋,徐云马.2008.重力资料对2008年汶川MS8.0地震的中期预测[J].国际地震动态,(7):36-39.
Kao H.,Shan S.J,2004.The Source-Scanning Algorithm:mapping the distribution of seismic sources in time and space[J].Geophys.J.Int.,157(2):589-594.