引言

地震定位是地震学中最基本、最经典的问题之一,其核心内容是最大限度地提高地震定位精度。影响地震定位精度的因素很多,比如记录台站的分布、震相拾取的精度、速度模型的准确性等,其中影响最大的是速度模型误差和震相到时读取误差。基于波形互相关技术的双差定位方法能够很好地克服这2种误差(黄媛,2008),提高了数字地震波形应用范围。使用互相关技术对P波和S波走时进行校正,减少到时读数误差的影响,很大程度上降低由于直接读取震相到时所造成的误差,并解决因速度模型不准确造成的定位结果分散问题,提高了地震定位的精度。

波形互相关分析分为时域和频域2种计算方法(Poupinet等,1984)。时域互相关分析操作简单,获得了较广泛的应用(Schaff等,20042005Waldhauser等,2000),针对垂直分量互相关计算的局限性,发展了时域多通道相关检测函数并用于计算波形互相关走时差(王清东等,2015)。在时域中选取互相关系数的阈值需要有一定的经验,如果阈值设定过高,可用数据就会很少,达不到高精度定位的要求,但如果阈值取太低,则会得到一些非真实震相信息。一个比较稳妥的方法便是使用双谱法(Du等,2004)对互相关系数结果进行验证,可以更有效地增强互相关的可用信息并提高互相关数据的可信度。双谱法能有效地压制高斯互相关噪声,在三阶谱域内同时计算原始和滤波后波形的时间延时,并用这2个时间延时来验证互相关系数的可靠性。因其具有精度高的特点,自提出以来便得到了一定的应用(Bannister等,2011Bourguignon等,2015Hansen等,2013张广伟等,2015赵翠萍,2006)。

本文首先利用双谱法软件包(BCSEIS)对三峡库区地震事件的波形数据进行互相关分析,并用双差定位方法对地震事件进行精定位,对不同类型数据的定位结果进行比较,获得了三峡库区高精度微震定位结果。

1 地震波形数据处理

三峡水库地震监测台网由26个地震台站组成,包括21个英国L-22型三分向短周期速度摆和5个美国Guralp CMG-3ESPC宽频地震计,数采均为Reftek130且采样率为200Hz(马文涛等,2010),在2009年3月至2010年12月的观测时间范围内一共获得了5275次地震事件观测报告及其波形数据。

在本研究中,我们首先对地震波形数据进行1.0—10Hz带通滤波,得到原始和滤波后波形数据,然后分别在P波和S波到时前后截取一段波形进行互相关分析。我们取P波时窗长128个采样点,P波到时前30个采样点,之后97个采样点;取S波时窗长192个采样点,S波到时前50个采样点,之后141个采样点。图 1显示了1个事件对的波形记录,2个地震事件的编号分别为242和244,图中只绘制了发震时刻到S波时窗范围内的波形记录,此外的波形不参与计算波形互相关系数,事件242和244的详细信息见表 1。波形互相关结算结果说明事件对之间具有高度相似,2次地震事件的波形互相关系数大于0.9,震相延时小于0.04s。将地震目录震相的到时差和波形互相关延时差进行比较(图 2),波形互相关延时可以降低由于震相拾取误差带来的偏差,能获得更高精度的到时差数据,说明用双谱验证方法提取的波形互相关系数具有较高的精度。


(a)地震242的波形记录,(b)地震244的波形记录
图 1 三峡水库加密监测台网记录到的事件对波形(表 1事件对) Fig. 1 Waveforms of event pairs from Three Gorges Reservoir (TGR) intensive seismic observation network (details are in Table 1)
表 1 图 1中地震事件对详细信息 Table 1 Detailed information of event pairs as shown in Fig. 1

黑色波形为事件242波形,红色为244波形,前两行中的P表示该台站地震目录P波走时,第三行中的CC为波形互相关系数,tdt为波形互相关P波到时差,tdiff为地震目录P波到时差
图 2 地震事件对在台站DTP(a)、LPT(b)的波形互相关示意图 Fig. 2 Waveforms showing correlation of event pair #1 from station DTP (a) and LPT (b)

在双谱计算过程中我们使用3个互相关系数阈值来检验计算结果:中央范围系数CClim、下限CClim (l)和上限CClim (u)。一般情况下$ C{C^{\lim (l)}} \le C{C^{\lim }} \le C{C^{\lim (u)}}$CClim相当于传统互相关延时计算中的阈值,在实际检测时我们使用以下原则:

(1)如果事件对 (ij) 最大互相关系数$ CC_{ij}^{\max } \ge C{C^{\lim (u)}} $,我们将接受互相关系数大于CClim (l)的所有台站的互相关计算值;反之,如果小于CClim (l)则否定该台站。即对于一个具有高度相似波形的地震事件对,我们将选择所有互相关系数高于CClim (l)的互相关延时评估结果通过双谱检验,即使有些互相关系数小于CClim,而这样的低于阈值标准的结果一些研究者是不予采用的。

(2)如果$ C{C^{\lim }} \le CC_{ij}^{\max } \le C{C^{\lim (u)}} $,我们将只对互相关系数高于CClim的台站进行验证。

(3)如果$CC_{ij}^{\max } $低于CClim,则地震事件对的互相关延时评估结果将会被丢弃。

我们使用下面的参数设置来验证互相关延时:$ C{C^{\lim \left(l \right)}} = 0.5 $$ C{C^{\lim }} = 0.7 $$ C{C^{\lim \left(u \right)}} = 0.9$${\Delta ^{\lim }} = 10 $,其中${\Delta ^{\lim }}$为双谱验证时带通滤波和天然波形数据之间最大延迟点数。为了与只用传统的互相关系数阈值估计方法相对比,我们选定互相关系数阈值为0.7,获得了传统的互相关延时数据,结果如图 3所示。用传统的互相关系数阈值方法估计得到1492711个P波震相对和749109个S波震相对,使用双谱验证法对互相关延时进行验证后,一共获得了1375372个P波震相对和624299个S波震相对,其中互相关系数大于0.7的P震相对1273958个,S波震相对596212个。通过双谱法验证后获得震相对数目要略少于传统的互相关震相对数据,一个主要的原因可能是双谱法验证剔除了延时较大的波形互相关数据。最后获得的S波的震相对个数少于P波,原因之一是辨识S波比P波难,此外S波受到P波尾波的干扰,降低了S波时窗内计算的互相关系数的质量。


(a)P波传统互相分析结果统计图,纵坐标括号内数值为互相关记录条数,互相关系数统计间隔为0.01;(b)P波双谱验证互相关分析结果统计图,纵坐标括号内的第一个数值为互相关系数大于0.7的记录个数,第二个为总记录个数,3条竖线表示了CClim (l)=0.5,CClim=0.7,CClim (u)=0.9;(c)S波传统互相关分析结果统计,其余同图(a);(d)S波双谱验证互相关分析结果统计图,其余同图(b)
图 3 三峡库区2种互相关系数评估柱状统计图 Fig. 3 Histogram of two kinds of cross-correlation coefficients in the Three Gorges Reservoir Area

在双差定位或双差层析成像方法中,地震对之间的最大距离对反演结果有很大的影响。如果该值取得越小,则建立联系的地震对之间距离越小,定位精度越高,但是太小则能够建立联系的地震就越少(黄媛等,2006)。在进行双差定位或双差层析成像之前,需要评价地震对之间的距离。一般可以计算地震对之间的距离并统计深度分布直方图,从而可以快速地确定合适的地震距,我们计算了研究范围内地震事件之间的距离(图 4),表明地震对之间的距离取20km较为合理,这个值与程序ph2dt给出的“强连接”事件对的最大距离也是一致的。但是该地震距是否与互相关地震对之间的距离匹配,需要进行详细的分析。根据双谱验证法获得的波形互相关数据,统计所有台站记录到的事件对震源距和相应的波形互相关系数,得到了互相关系数与事件对震源距分布图(图 5)。结果表明,单纯地设置地震对之间距离阈值为20km会降低互相关数据的使用质量,定位结果受横向不均匀性的影响较大,因此需要适当地减小互相关地震对之间距离的阈值,这里我们选取6km,以保证大多数互相关系数和距离的一致性,并确保互相关数据的使用效率和定位的精度。


图 4 三峡库区地震目录的震源距离统计图 Fig. 4 Plot of total number of event pairs and separation distance of hypocenters

(a)P波数据分布图;(b)P波互相关系数与震源距密度分布图,横坐标为各个台站记录到的同一地震对的震源距离,纵坐标为相应地震对的互相关系数,其中密度统计图统计间隔震源距为2km,互相关系数为0.1;(c)和(d)为S波分布图,其余同(a)和(b)。密度为对应矩形色块内的数据量与总数据量的比值,衡量的是色块内数据占总数据的比重,无量纲。
图 5 台站地震对和震相互相关系数分布图 Fig. 5 Distribution of event pairs and cross-correlation coefficient for all stations
2 三峡库区双差定位实现

根据2009年3月至2010年12月中国地震局地质研究所三峡水库加密台网观测数据结果,一共得到5275个地震事件,根据地震目录获得P波走时数据43538条和S波走时数据43385条,平均每个地震事件有7个清晰台站记录。为排除震相判读错误和其他信息的干扰,使用最小二乘拟合走时曲线,剔除误差较大的数据,同时使用和达曲线(即TS-TP曲线)对结果进行检验(蔺永等,2014)。图 6表明数据拟合前存在一系列的干扰数据,对数据进行走时曲线最小二乘拟合,剔除误差较大的震相后,得到误差较小的地震事件信息。经过挑选后,得到5244个地震事件,共41423条P波和42354条S波震相数据。图 7为5244个地震事件的信息分布图,分别给出了沿着经度、纬度和深度的地震事件分布图。结果显示,地震事件主要分布在8km深度以上的范围,且地震分布相对集中,长江水系外围分布有离散地震事件。


(a)为最小二乘剔除误差较大震相后走时曲线;(b)为初始震相走时曲线;(c)为剔除误差较大震相后的TS-TP曲线,红线为直线拟合结果vP/vS为1.73;(d)为初始Ts-Tp曲线
图 6 地震走时曲线和和达曲线图 Fig. 6 The travel time curves and Wadati plot

(a)地震事件初始震中分布图;(b)震源深度沿着纬度方向的剖面图;(c)震源深度沿经度方向的剖面图;(d)震源深度统计图
图 7 5244个地震事件的初始震中分布图、深度剖面图和统计图 Fig. 7 Plots of distribution of space, depth and statistical results of 5244 seismic events

根据三峡库区已有的速度模型研究结果(李强等,2009赵旭等,2007),选定速度模型的水平深度分别取0、2、5、8、11、14和20km,对应的P波速度为4.8、5.4、5.65、5.8、6、6.15和6.5km/s,vP/vS值根据图 6的拟合结果取为1.73。

3 三峡水库地震精定位结果及分析

根据BCSEIS波形互相关分析的结果和地震目录数据,我们使用hypoDD双差定位软件对三峡库区水库地震事件进行3组测试:第一组,仅利用地震目录数据进行定位(CAT);第二组,利用标准的波形互相关方法获得的波形延时和地震目录到时差数据(CC+CAT)进行定位;第三组,利用双谱验证获得的波形延时数据和地震目录到时差数据(WCC+CAT)进行定位。为了对比定位结果的差异性,对3组测试数据设置相同的参数和速度模型,最终获得不同数据类型的地震平面位置(图 8)和深度误差分布图(图 9)。图 8表明,第二组和第三组的定位结果明显优于第一组数据的结果,具有较小的定位误差和优势深度(图 9)。


(a)初始地震事件震中分布;(b)仅地震目录数据的hypoDD定位结果;(c)结合传统互相关分析与地震目录数据的hypoDD定位结果;(d)结合双谱验证互相关分析与地震目录数据的hypoDD定位结果。左上角数字为使用不同数据时获得的重定位地震事件数
图 8 不同数据类型的hypoDD定位结果 Fig. 8 Map view of relocation results from different type data by HypoDD

(a)仅使用地震目录数据;(b)结合传统互相关分析与地震目录数据;(c)双谱法互相关分析与地震目录数据hypoDD定位结果,上图为经度方向的深度分布图,下图为经度方向的误差分布图
图 9 不同数据类型的hypoDD定位结果的深度和RMS分布图 Fig. 9 The distribution of depth and RMS from different data by hypoDD relocation

重定位后的地震震中分布更加集中,呈线性化分布,显示出地震事件向内收敛的趋势。在巴东神龙溪两岸表现出3条明显的近东西向条状分布,震源深度较浅,平均深度在5km左右,震中分布与地层走向基本一致,进一步验证了水库蓄水后库水从神龙溪等地下暗河渗入而诱发地震的解释(马文涛等,2010)。中间一条地震分布线东口于2013年12月16日发生了MS5.1级地震,从2009年3月—2010年12月的微震定位结果可以说明,该地震可能是在库水的作用下,贯穿整个东西向的小规模断裂引发的较大地震活动。在泄滩西区域存在的一个陡立东西向的微震活动带,其可能和长江水系的渗透有关。而在秭归香溪河口附近,地震分布在仙女山断裂带上,震源深度向下延深到10km,揭示了库水渗透的最大深度。该处于2014年3月发生2次4级以上地震,可能也是库水渗透引起断裂局部活动的结果。

3组不同类型数据的定位误差统计结果(图 10)表明使用双谱验证获得波形互相关延时与地震目录到时差数据的hypoDD定位精度明显要高于其他2种数据的精度,地震震源位置的平均误差在东西向为3.2m、南北向为3.9m、垂直向为6.2m。


(a)—(c)为地震目录定位结果,(d)—(f)为结合传统互相关分析和地震目录定位结果,(g)—(i)为结合双谱验证波形互相关和地震目录定位结果
图 10 不同数据类型的双差定位结果误差统计 Fig. 10 The histograms of statistical errors of different type data by hypoDD relocation
4 结论

本文应用双谱分析方法,对三峡库区2009年3月至2010年12月的地震观测数据进行了波形互相关分析,获得了高质量的双谱验证波形互相关数据,并使用双差定位方法对地震进行了重定位,获得了高精度的定位结果。

基于三阶谱域的双谱验证波形互相关算法能获得比传统阈值判定更高质量的波形延时数据,得到的互相关地震事件对的距离更加符合实际情况,能为双差定位或者双差层析成像提供更高质量的数据资料。

根据双谱验证方法获得的传统阈值判定和双谱验证的波形互相关数据,分别使用双差定位程序进行了定位,结果表明结合双谱验证的波形互相关数据的定位结果精度高于其他数据的结果,其东西向震源位置平均误差为3.2m、南北向为3.9m、垂直向为6.2m。

对三峡水库地震重定位结果表明,巴东库区神龙溪两岸微震分布明显呈现出3条近东西向的线性分布特征,与新构造时期产生的小规模断裂和碳酸盐岩地层走向一致,揭示了库水主要沿着溶洞或地下暗河渗透而诱发地震活动,其中中间1条地震分布线端点正是发生2013年12月16日湖北巴东5.1级地震的震中位置,说明它是在库水作用下沿着东西向小规模断裂的一次较大岩体错动事件;泄滩西区东西向垂直条带状分布预示地震活动与长江库水渗透之间的关系;而在秭归香溪河口附近,地震分布在仙女山断裂带上,震源深度向下延深到10km,揭示了库水渗透的最大深度。

致谢: 评审专家对本文的修改提出了宝贵的意见,文中大部分图件使用GMT绘图软件包绘制(Wessel等,1995),在此一并表示感谢。
参考文献
黄媛, 2008. 结合波形互相关技术的双差算法在地震定位中的应用探讨. 国际地震动态(4): 29–34.
黄媛, 杨建思, 张天中, 2006. 2003年新疆巴楚-伽师地震序列的双差法重新定位研究. 地球物理学报, 49(1): 162–169.
李强, 赵旭, 蔡晋安, 等, 2009. 三峡水库坝址及邻区中上地壳P波速度结构. 中国科学D辑:地球科学, 39(4): 427–436.
蔺永, 马文涛, 2014. 基于Matlab的三峡水库地震数据处理与分析. 震灾防御技术, 9(3): 447–453. DOI:10.11899/zzfy20140311
马文涛, 徐长朋, 李海鸥, 等, 2010. 长江三峡水库诱发地震加密观测及地震成因初步分析. 地震地质, 32(4): 552–563.
王清东, 朱良保, 苏有锦, 等, 2015. 2012年9月7日彝良地震及余震序列双差定位研究. 地球物理学报, 58(9): 3205–3221. DOI:10.6038/cjg20150916
张广伟, 雷建设, 2015. 2011年云南腾冲5.2级双震发震机理. 地球物理学报, 58(4): 1194–1204. DOI:10.6038/cjg20150409
赵翠萍, 2006. 1997-2003年新疆伽师震源区特征的地震学方法研究. 北京: 中国地震局地球物理研究所.
赵旭, 李强, 蔡晋安, 2007. 三峡库首区最小一维速度模型研究. 大地测量与地球动力学, 27(专刊): 1–7.
Bannister S., Fry B., Reyners M., et al, 2011. Fine-scale Relocation of Aftershocks of the 22 February MW6.2 Christchurch Earthquake using Double-difference Tomography. Seismological Research Letters, 82(6): 839–845.
Bourguignon S., Bannister S., Henderson C. M., et al, 2015. Structural heterogeneity of the midcrust adjacent to the central Alpine Fault, New Zealand:Inferences from seismic tomography and seismicity between Harihari and Ross. Geochemistry, Geophysics, Geosystems, 16(4): 1017–1043. DOI:10.1002/2014GC005702
Du W. X., Thurber C. H., Eberhart-Phillips D., 2004. Earthquake relocation using cross-correlation time delay estimates verified with the bispectrum method. Bulletin of the Seismological Society of America, 94(3): 856–866. DOI:10.1785/0120030084
Hansen S. E., DeShon H. R., Moore-Driskell M. M., et al, 2013. Investigating the P wave velocity structure beneath Harrat Lunayyir, northwestern Saudi Arabia, using double-difference tomography and earthquakes from the 2009 seismic swarm. Journal of Geophysical Research:Solid Earth, 118(9): 4814–4826. DOI:10.1002/jgrb.50286
Poupinet G., Ellsworth W., Frechet J., 1984. Monitoring velocity variations in the crust using earthquake doublets:An application to the Calaveras Fault, California. Journal of Geophysical Research, 89(B7): 5719–5731. DOI:10.1029/JB089iB07p05719
Schaff D. P., Bokelmann G. H. R., Ellsworth W. L., et al, 2004. Optimizing correlation techniques for improved earthquake location. Bulletin of the Seismological Society of America, 94(2): 705–721. DOI:10.1785/0120020238
Schaff D. P., Waldhauser F., 2005. Waveform cross-correlation-based differential travel-time measurements at the Northern California Seismic Network. Bulletin of the Seismological Society of America, 95(6): 2446–2461. DOI:10.1785/0120040221
Waldhauser F., Ellsworth W. L., 2000. A double-difference earthquake location algorithm:Method and application to the northern Hayward fault, California. Bulletin of the Seismological Society of America, 90(6): 1353–1368. DOI:10.1785/0120000006
Wessel P., Smith W. H. F., 1995. New version of the generic mapping tools. Eos, Transactions American Geophysical Union, 76(33): 329.