(1)
摘要 利用新近发展起来的台阵网格自动事件定位方法, 对鄂尔多斯地块及秦岭‒大别造山带区域 136 个宽频带流动台站记录的连续波形进行扫描和定位。通过与地震目录的对比发现, 不仅在常规地震带上分布着大量的事件, 而且在地震目录中极少发生地震的区域出现大量地震事件集中分布的现象。这些事件主要分布于鄂尔多斯地块南部以及秦岭‒大别造山带附近的断层带。经过对定位台站的抽样分析和对定位事件地震波形的仔细人工确认, 验证了这些被扫描并定位的事件的真实性。通过慢度分析, 获得所定位事件的面波传播视速度, 绝大部分事件对应的视速度在 2.5～5.2km/s 之间, 与地震面波群速度相符。所定位的地震事件可用于补充秦岭‒大别造山带的地震目录, 同时可以通过进一步的波形拟合来分析其震源机制。
关键词 识别和定位; 面波; 地震目录; 事件分布; 视速度
地震事件定位是地震学中的经典问题。早期的研究中大多利用体波的到时进行地震定位, 定位方法包括绝对定位和相对定位。绝对定位是经典的定位方法, 其中的典型代表是 Geiger 法[1]。在 Backus等[2]反演论的基础上, Klein[3]提出 HYPOINVERSE方法并得到广泛使用。相对定位法可以提高地震定位的精度, 其中应用最广泛的是 Waldhauser 等[4]提出的双差定位方法。该方法利用不同台站之间的相对到时差, 更加精确地确定地震之间的相对位置。
微震事件震级小, 信号弱, 其地震记录信噪比低, 很难直接从波形中拾取震相到时。对于震级较低的震颤等事件, 地震目录中也鲜有记录, 因此很难了解没有或少有较大地震区域的地下信息[5‒6]。针对这种情况, 模板匹配技术被发展出来, 用于定位能量较低的地震事件[7‒9]。模板匹配技术利用波形互相关来检测具有相似波形的微弱信号, 可以较有效地侦测低信噪比事件。随着地震台站密度的增加, 台阵数据成像技术应用于地震定位工作中。汶川地震后, 一些学者利用反投影方法进行余震定位, 取得较好的效果[10‒12]。
de Groot-Hedlin 等[13‒14]提出一种新的台阵网格自动定位(automated event location using a mesh of arrays, AELUMA)方法, 利用密度足够大的地震台网来探测和定位地球物理事件。该方法利用大阵列数据的相干性, 确定台站得到地震信号的方位角以及到达时间等信息, 并依此进行准确的震源定位。Fan 等[15‒16]利用此方法, 对美国东西海岸的风暴及其引发的风暴地震和墨西哥湾的滑坡地震进行了系统的定位和研究, 取得良好的效果。
秦岭‒大别造山带是华北克拉通与华南克拉通的缝合带。与其南北两边的四川盆地和鄂尔多斯地块相比, 渭河地堑和秦岭‒大别造山带的岩石圈厚度更薄, 是青藏高原上地幔受到印度板块挤压而向东移动的通道[17‒20]。然而, 中国台网正式地震目录(https://data.earthquake.cn/)显示, 2000 年 1 月以来, 秦岭‒大别造山带的地震活动性偏弱。本文基于北京大学和中国地质科学院联合布设的位于鄂尔多斯地块南部以及秦岭‒大别造山带附近的流动地震台站的观测数据, 对该区域进行事件的扫描和定位。
模板匹配技术需要研究区域内具有足够的已知地震来为未知事件的定位提供模板, 但是本文研究区域内已知地震事件很少, 不能提供足够的模板。AELUMA 方法仅仅利用台阵的地震记录, 同样可以进行全局扫描, 不需要更多先验信息和地震模板, 因此本文采用 AELUMA 方法进行地震事件的探测和定位工作。
AELUMA 方法利用长周期(10～40s)连续信号的相干性分析来确定地震波方位角、到时和局部相速度等信息, 从而确定事件的发生时刻和震源位置[15]。该方法不需要了解震源类型, 也不需要精确的地下速度结构模型, 可以直接寻找并定位震源, 并且计算效率较高。同时, 该方法可以利用振幅较小的数据进行扫描定位, 可以提高对释放能量较少的微震事件的探测效果。该方法只在平面内, 根据台站记录的相干性处理数据, 从而减少台站高程对结果的影响, 使问题简单化, 并且因地表面波信号较为丰富, 可以增加所用观测数据的信噪比。但是, 由于我们利用面波进行扫描和定位, 也导致一些震源较深的地震无法扫描到, 造成遗漏。然而, 大陆地震多发生于地壳的中上部, 深度较浅, 因此适合用此方法进行扫描。
AELUMA 方法的基本思想是把台阵划分为较小的三角形网格, 并挑选每个网格检测到的可靠信号, 再根据这些信号的慢度和方位, 定位一个事件的震源位置和时间, 包括阵列离散化、信号检测、分组分析和震源定位 4 个步骤。
第一步, 阵列离散化。利用 Delaunay 三角剖分方法, 将大阵列台站划分成 3 个台站组成的三元组[21]。为了保证最佳定位效果, 仅使用边长在 10～ 600km 之间, 内角在 30°～120°之间的三元组。每个三元组都可以检测到通过该子阵列的相干信号, 并且得到信号的到达时间、局部相速度以及方位角等信息。
第二步, 信号检测。利用 τ-p 方法来检测每个三元组中的相干信号, 并利用相干性来确定重心位置上相干信号的到达时间、传播方向和局部相速度等, 计算公式[22]如下:
(1)
(2)
(3)
其中, θ 和
分别为方位角和相速度; Ti, j 为台站 i 与j 之间的时间延迟;
为水平慢度矢量; 
是台站i与j之间的距离。
处理实际信号时, 将每天的波形分割成 600s的时窗, 以 180s 为步长进行滑动。为保证波形的相干性, 仅保留 3 个互相关系数平均值大于 0.5 的三元组的信号。经过处理后, 得到合格的信号检测结果。
黑色三角形表示台站位置, 蓝色圆点表示地震目录中事件位置, 灰色线条表示断层
图1 本文使用的台阵位置以及2012—2013年地震目录事件位置
Fig. 1 Location of the array used in this paper and the location of events in the 2012‒2013 earthquake catalog
第三步, 分组分析。为保证相干信息的效果, 需要根据台站到时, 将所有三元组检测信号分为不同的簇集。首先得到三元组之间的距离, 然后将不同距离范围的三元组检测信号分为不同的簇集, 再根据每个簇集的探测结果进行震源定位。
第四步, 震源定位, 分两步进行。首先定位每个簇集的探测震源点, 即利用每个簇集的三元组所检测到的传播方位角和到达时间, 估算潜在震源点的大致范围; 然后用更密的网格在此范围内进行搜索, 结合所有簇集, 求得最佳震源探测点。
为了使数据尽可能准确, 需要选取台站数量和台站间距合适的台网数据进行研究。如图 1 所示, 本文选择鄂尔多斯地区 136 个台站组成的台站阵列, 通过对 2012—2013 年期间数据连续性好的 224 天台站地震数据扫描, 进行地震事件的探测和定位。首先对所有台站进行德劳内三角剖分, 剖分结果见图2。对宽频带台站的垂直分量连续数据进行 10～40 s 的带通滤波, 确保我们使用的是面波成分。由于使用的面波频段较低, 波长较长, 因此面波信号受地球介质小尺度横向非均匀性的影响较小, 可以忽略地球速度模型对探测结果的不利影响。
红色三角形为台站位置, 蓝线组成的三角形为台站阵列剖分所得三元组, 本文求得的参数均位于三元组重心处(黑色圆点)
图2 Delaunay三角剖分方法对台站阵列的离散化结果
Fig. 2 Discretization results of station array by delaunay triangulation method
在所用有效数据范围内, 每天定位的事件数目从十几个到六十几个不等。以 2013 年 9 月 22 日为例, 我们将当天 120 个有效台站的连续地震数据分成 1439 个时间窗, 利用 AELUMA 方法得到 1860 个检测点。移除 623 个重复检测点和 217 个孤立检测点后, 最终保留 1020 个有效检测点, 形成 200 个有效的三元组。一天内所有检测点的振幅、相关系数、相速度和方位角等参数与时间之间的对应关系如图 3 所示, 可以看到在一些时间点有很多检测点排列整齐, 因此确定在这些时刻存在地震事件 信号。
根据上述信息, 在台站周围用 AELUMA 方法扫描定位出 41 个事件及其对应的发生时刻, 这些事件的位置如图 4 所示。与图 4 中星号显示的台网中心地震目录对比, 可以看到两者并不重合。这一现象主要是由两者所用台站分布和数据类型不同造成的。值得注意的是, 本文定位的大部分事件都不能在地震目录中查询到相应的记录, 因此本研究聚焦于这些未被地震目录收录的事件。
在共 224 天时间里, 扫描并定位出 9582 个事件(图 5 中篮色圆点所示), 这些事件主要密集地分布在青藏高原东南缘、太行山、鄂尔多斯周缘等地震带。经过与中国台网中心地震目录中地震事件位置(图 1)的对比, 我们发现在地震目录中很少有地震事件记录的一个区域(31.9°—34.4°N, 109°—115°E)存在大量定位事件。这一区域为秦岭‒大别造山带, 其中有多条断层发育。
为了更加直观地观察事件的发生密度, 分别绘制定位事件和地震目录的散点密度图(图 6)。将图 6的区域均分为 60×60 个网格, 每个网格的大小为0.35°×0.22°, 按照每个网格内的事件数量计算事件的密度。对事件密度图进行分析后可以确定, 利用AELUMA 方法在研究区域扫描并定位出的事件密度是足够的, 并非一些事件偶然集中在这一区域, 足以说明有必要对该区域的事件做更多的探索。这些事件集中分布在鄂尔多斯地块南部边缘和秦岭‒大别造山带的断层带, 说明它们与该地区的断层分布联系紧密。另外, 对比地震目录与定位事件的密度分布可以发现, 台站阵列范围内的定位事件非常丰富, 距离台站阵列越远, 则定位事件分布越稀疏, 甚至在一些地震目录中地震事件较密集的区域也鲜有定位事件。因此, 我们认为 AELUMA 方法对台站阵列的分布比较敏感。在阵列内部, 台阵接收的信号覆盖范围广, 定位效果较好; 靠近阵列边缘处, 台阵接收的信号覆盖范围缩小, 定位结果稍为逊色;在远离阵列的地方, 由于信噪比降低, 有用信息易被覆盖, 接收信息的角度也较为单一, 导致定位效果明显下降, 并且离阵列越远, 效果越差。
图3 对2013年9月22日的地震数据处理后得到的所有检测点振幅、相关系数、相速度和反方位角与时间的对应关系
Fig. 3 Corresponding relationship between amplitude, correlation coefficient, phase velocity and inverse azimuth angle of all detection points obtained after processing and analysis of seismic data recorded on September 22, 2013
图4 2013 年 9 月 22 日定位事件和地震目录事件位置分布
Fig. 4 Location distribution of events and seismic catalog events on September 22, 2013
为了观察不同台站对定位事件的影响, 我们先后 4 次随机选取 2013 年 9 月 22 日的 100 个台站数据进行定位, 得到图 7 中 4 个结果。可以发现, 由于完全随机选取台站数据, 台站分布不同导致这 4 次结果不完全重合, 但台站阵列内部的定位结果非常接近, 事件的相对位置小于 5km。同时可以看到, 由于随机选取的台站数目减少, 定位事件的总数也比图 4 中的少。这也从侧面说明台阵分布和台站数量都影响着最终的定位效果。
AELUMA 方法不能直接给出事件的震级, 震级 ML 的确定需要已知区域的地震波衰减参数。因此, 本研究采用相对震级来估算检测到的事件的震级。我们挑选地震台网目录中的一个事件(2012 年8 月 21 日 08:41:20 于 34.092°N, 114.675°E 处发生的2.3 级地震)作为参考事件, 拾取离参考事件最近的5 个台站的振幅, 将其几何扩散校正后的平均值 A0与其他事件相应的振幅 A 进行比较, 利用公式
ML =lg(A/A0 +M0)(4)
图5 研究时段内定位事件的分布
Fig. 5 Distribution of located events in the study time period
(a)定位事件; (b)地震目录事件
图6 本文定位地震事件对数密度分布
Fig. 6 Logarithmic density distribution of located seismic events in this study
图7 4次随机选取2013年9月22日的100个台站数据进行定位的不同结果
Fig. 7 Different results of 100 stations that were selected randomly four times on September 22, 2013
估算其他事件的震级。
利用式(4)估算选取出的 1000 个定位事件的震级, 结果如图 8 所示。可以看出, 事件的震级集中在 1.5～2.5 级之间, 低于 1 级和高于 3 级的事件很少, 说明探测到的主要是低能量的事件, 这也是在本地区正式地震目录中地震数目较少的原因。宽频带流动台阵的观测可以在一定程度上增加本地区地震目录的完备性。
为了进一步验证在秦岭‒大别造山带识别的大量定位事件的真实性, 我们针对研究区域的所有事件的发震时刻, 截取 300s 的波形, 按照依震中距排列的地震波形及相应上午 τ-p 域结果进行分析, 对定位事件进行确认。
图8 估算震级的分布范围
Fig. 8 Estimate the magnitude distribution range
以 2013 年 9 月 22 日 06:05:29 在 33.82°N, 113.33°E 处的一个定位事件(图 4 中蓝色五角星所示)为例。以 06:05:29 为零时刻, 截取所有台站的波形, 经过1～2Hz 滤波, 依震中距排序后, 得到如图 9 所示的波形。可以看出, 在震中距 0～300km 范围内, 存在走时近似一条直线的面波信号。
将图 9 中震中距在 0～300km 范围内的地震波形数据进行 τ-p 变换, 得到 τ-p 域图像(图 10)来进一步确定事件的发震时刻和射线参数。图 10 中在 τ=0, p=0.3s/km 处存在能量峰值, 表明在 0 时刻存在群速度为 3.33km/s 的面波。分析研究区域内 1000 个定位事件, 其中 93%可以在波形图和 τ-p 图像中观察到明显的事件特征。因此, 可以确信这些事件是真实的。
图9 北京时间2013年9月22日的示例事件的波形
Fig. 9 Waveforms of sample events on September 22, 2013 (Beijing time)
图10 图9中波形数据的τ-p域能量分布
Fig. 10 Energy distribution in τ-p domain of waveforms in Fig. 9
图11 事件视速速度统计结果
Fig. 11 Statistics of apparent velocities of located events
根据 τ-p 域分析结果, 进一步统计研究区域定位所得事件地震波传播的视速度(图 11), 可以发现地震波传播的视速度大部分在 2.5～5.2km/s 之间, 其中 70%的事件波速在 3.0～4.0km/s 之间, 出现几率最高的视速度为 3.45 km/s, 符合真实地球介质中面波的群速度范围。极个别事件计算得到很大的视速度值, 符合近似垂直入射远震体波的视速度范围, 不符合区域地震的传播特征, 我们在最终的地震目录中去除了这部分结果。
本文利用近年发展起来的 AELUMA 方法, 对鄂尔多斯地块中南部的连续地震记录进行扫描和定位, 并与已有地震目录进行对比, 得到如下结论。
1)在鄂尔多斯地块南部及秦岭‒大别造山带定位出 1780 个事件, 这些事件多分布在块体的边界处, 与断层分布吻合, 因此多为构造地震事件。
2)在地震目录之外, 鄂尔多斯地块东南的秦岭‒大别地区存在大量事件, 波形分析和 τ-p 变换结果确认这些事件是真实存在的。
3)研究区域内出现的大量事件的速度分布集中在 2.5～5.2km/s 之间, 出现几率最大的视速度为3.45km/s, 处于面波的正常速度范围。
4)本研究在地震活动较少的区域发现的这些真实事件及其发生时刻、发生位置和速度特征可以很好地补充现有的地震目录, 为完善地震目录以及研究地震事件的分布和活动性提供依据。
AELUMA 方法在事件的扫描和定位方面有不错的效果, 但不能给出事件的震源机制和深度信息, 并且定位效果对台阵的覆盖范围和台站数量要求较高, 不适用于所有地区。后续研究中, 可以通过波形拟合等方法进一步研究地震事件更全面的信息。
参考文献
[1]Geiger L. Probability method for the determination of earthquake epicenters from the arrival time only. Bulletin of St Louis University, 1912, 8(1): 56‒71
[2]Backus G, Gilbert F. The resolving power of gross earth data. Geophys J R Astr Soc, 1968, 16: 169‒205
[3]Klein F W. Hypocenter location program HYPOIN-VERSE [EB/OL]. (1978) [202‒10‒16]. https://pubs.er. usgs.gov/publication/ofr78694
[4]Waldhauser F, Ellsworth W L. A double-difference earthquake location algorithm: method and application to the northern Hayward fault, California. Bull Seism Soc Am, 2000, 90(6): 1353‒1368
[5]Dragert H. Episodic tremor and slip in the Cascadia Subduction Zone: a story of discovery // AGU Fall Meeting Abstracts. San Francisco, 2003: G22F-01
[6]Porritt R W, Allen R M, Boyarko D C, et al. Inves-tigation of Cascadia segmentation with ambient noise tomography. Earth & Planetary Science Letters, 2011, 309(1/2): 67‒76
[7]Gibbons S J, Ringdal F. The detection of low mag-nitude seismic events using array-based waveform correlation. Geophysical Journal International, 2006, 165(1): 149‒166
[8]Gibbons S J, Sørensen M B, Harris D B, et al. The detection and location of low magnitude earthquakes in northern Norway using multi-channel waveform correlation at regional distances. Physics of the Earth and Planetary Interiors, 2007, 160(3/4): 285‒309
[9]Peng Z G, Zhao P. Migration of early aftershocks following the 2004 Parkfield earthquake. Nature Geo-science, 2009, 2(12): 877‒881
[10]Kiser E, Ishii M. Hidden aftershocks of the 2011 Mw 9.0 Tohoku, Japan earthquake imaged with the back-projection method. Journal of Geophysical Research Solid Earth, 2013, 118: 5564‒5576
[11]Roten D, Miyake H, Koketsu K. A Rayleigh wave back-projection method applied to the 2011 Tohoku earthquake. Geophysical Research Letters, 2012, 39 (2): L02302
[12]黄媛, 吴建平, 张天中, 等. 汶川8.0级大地震及其余震序列重定位研究. 中国科学: 地球科学, 2008, 38(10): 1242‒1249
[13]de Groot-Hedlin C D, Hedlin M A H, Walker K T. Detection of gravity waves across the US Array: a case study. Earth planet Sci Lett, 2014, 402: 346‒352
[14]de Groot-Hedlin C D, Hedlin M A H. A method for detecting and locating geophysical events using groups of arrays. Geophys Journal International, 2015, 203(2): 960‒971
[15]Fan W, McGuire J, de Groot-Hedlin C D, et al. Stormquakes. Geophysical Research Letters, 2019, 46(12): 909‒918
[16]Fan W, McGuire J, Shearer P M. Abundant spontane-ous and dynamically triggered submarine landslides in the Gulf of Mexico. Geophysical Research Letters, 2020, 47: e2020GL087213
[17]Huang Z, Xu M, Wang L, et al. Shear wave splitting in the southern margin of the Ordos block, North China. Geophys Res Lett, 2008, 35: L19301
[18]Chang L, Wang C, Ding Z. Upper mantle anisotropy in the Ordos Block and its margins. Sci China Earth Sci, 2011, 54: 888‒900
[19]Li J, Wang X J, Niu F L. Seismic anisotropy and implications for mantle de-formation beneath the NE margin of the Tibet plateau and Ordos plateau. Phys Earth Planet Inter, 2011, 189: 157–170
[20]Yu Y, Chen Y J. Seismic anisotropy beneath the southern Ordos block and the Qinling-Dabie orogen, China: eastward Tibetan asthenospheric flow around the southern Ordos. Earth and Planetary Science Letters, 2016, 455: 1‒6
[21]Lee D T, Schachter B J. Two algorithms for cons-tructing a Delaunay triangulation. Int J Comput Sci, 1980, 9(3): 219‒242
[22]Fan W, de Groot-Hedlin C D, Hedlin M A H, et al. Using surface waves recorded by a large mesh of three-element arrays to detect and locate disparate seismic sources. Geophys Journal International, 2018, 215: 942‒958
Application of the Automated Event Location Using a Mesh of Arrays in Southern Ordos and Qinling-Dabie Orogenic Belt
Abstract The recently developed automated event location using a mesh of arrays (AELUMA) method is used to scan and locate the continuous waveforms recorded by 136 broadband mobile seismic stations in Ordos basin and its surrounding region. Compared with the seismic catalog provided by national earthquake data center, besides the events located in conventional seismic zone, a large number of events were located in Qingling-Dabie orogenic belt with few earthquakes in the catalog. The detected events were further confirmed through the bootsctrap analysis and carefully manual examination. The apparent velocities of these events were also obtained from slowness analysis. Most apparent velocities were 2.5‒5.2 km/s which was consistent with the group velocity of surface wave propagation. The events detected and located in this paper can be a supplementary of the earthquake catalog and further used to determine the focal mechanism through waveform fitting.
Key words detect and locate; surface wave; seismic catalog; events distribution; apparent velocity
doi: 10.13209/j.0479-8023.2022.045
国家自然科学基金(41874071)和地震行业科研专项经费(201408013)资助
收稿日期: 2021-05-29;
修回日期: 2021-06-26