Research on Regional Stacking of Broadband Seismic Records at Hongshan Station in Hebei Province
-
摘要: 可靠的震相走时是地震预警技术中精确测定震源位置和发震时刻的基础，本文运用STA/LTA震相识别技术，针对单台（河北红山台）2009—2021年共计12年积累的地震记录进行叠加计算，得到了红山台记录到的区域地震各震相走时曲线。结果显示，震中距0°～50°范围内红山台共成像7种震相的走时曲线，分别为P、S、PP、SS、PcS、ScS以及R面波震相，且随着组合参数变化，叠加成像的震相种类、震中距范围、清晰度均有所不同。此外，通过绘制各震相走时曲线发现，震中距0°～15°范围内，P波、S波及R波走时曲线基本呈线性变化，震中距0°～15°范围内计算得到红山台区域地震P波传播速度为7.5 km/s左右，S波传播速度为4.2 km/s左右，R波传播速度为3.5 km/s左右，介于P波和S波之间存在一个震相的走时痕迹，波速为5.4 km/s左右。本工作对于提升红山台震中距≤1000 km的地震预警定位精度有指导意义。Abstract: Reliable seismic phase travel time plays a fundamental role in the accurate determination of source position and occurrence time in earthquake early warning technology. In this paper, the STA/LTA phase identification technology is used for the first time to stack the seismic records accumulated by a single station (Hebei Hongshan station) for 12 years from 2009 to 2021, and the travel time curves of each seismic phase recorded by Hongshan station are imaged. The results show that the travel time curves of seven known seismic phases in Hongshan station are P, S, PP, SS, PcS, ScS, and R surface wave phases in the range of 0° to 50° epicenter distance. With the change of combination parameters, the types of seismic phases, epicentral distance range, and definition of superimposed imaging are different. In addition, by drawing the travel time curves of each seismic phase, it is found that the travel time curves of P wave, S wave, and R wave basically change linearly within the epicenter distance of 0° to 15°. Within the epicenter distance of 0° to15°, it is calculated that the seismic P-wave propagation velocity in Hongshantai area is about 7.5 km/s, S-wave propagation velocity is about 4.2 km/s and R-wave propagation velocity is about 3.5 km/s. There is a trace of seismic phase travel time between P wave and S wave, and the wave propagation velocity is about 5.4 km/s.
-
Key words:
- Seismic phase /
- Superposition /
- STA / LTA /
- Seismic record
-
引言
地震波中各震相走时曲线对于精确测定震源位置和发震时刻具有基础性作用。随着地震台网的快速发展，地震监测设备的种类和地震事件记录越来越丰富，如何更好的运用地震数据识别震相到时，充分减少噪声对地震信号识别的干扰成为地震工作中面临的一个问题。地震资料叠加可以有效提高地震记录中的信噪比，从而更好识别震相。
1996年Astiz等人（1996）以STA/LTA震相识别技术为基础，通过对全球地震台网记录到的MS5.5以上地震记录的STA/LTA值进行叠加，成像了全球地震记录的各震相走时图。Shearer（1991）通过运用STA/LTA震相识别技术对三分量长周期地震记录进行叠加，成像了全球震相走时图，通过对比长周期叠加成像与国际地震中心（ISC）短周期数据列出的走时，发现PcP、PKKP和PKPPKP等震相可以在ISC数据中观察到，但在长周期叠加中并不明显。Earle（1999）利用叠加技术对32000幅三分量宽频带地震记录进行了地球远震波场的偏振成像，发现低频和高频记录剖面显示了远震波场的不同特征。国外学者（Walck等，1984；Okal等，1998；Shearer等，2011）运用叠加方法针对天然地震进行了研究，并产出了很多有意义的成果。
秦满忠等（2014）使用兰州小孔径地震台阵记录的近10年地震观测垂直分量波形数据，利用STA/LTA值叠加出适用于青藏高原东北缘地区的观测走时曲线。其他学者运用叠加技术研究天然地震的情况较少，且没有运用该技术研究单台震相记录情况的先例，叠加技术通常应用于地球物理勘探中，如岳玉波等（2020）针对共中心点（CMP）叠加成像方法不适用震源和接收点之间存在高程差的问题，提出了一种基于时变数据映射的共反射点（CRP）叠加成像方法。蒋振武等（1995）为了对有倾角地层的井间反射波进行叠加成像，同时减少对初始模型的依赖，提出了一种新方法（DLCDP）。部分学者通过运用不同的叠加方法解决了科研中的实际问题（冯太林等，2001；芦俊等，2018；李希元等，2020）。
目前，用于震相分析定位的地震走时曲线是以全球平均速度模型为基础，适合大区域范围的地震观测资料，无法顾及小区域地质结构的特殊性。本文首次运用STA/LTA震相识别技术，针对单台（河北红山台）多年积累的宽频带观测数据进行叠加计算，成像了该测震台0°～50°范围内记录到的各震相的走时曲线，这些走时曲线不再按照全球或区域的平均走时绘制，而是按照单台的实际到时绘制，客观上更能代表红山台记录各震相的走时特征，为后期研究单台震相到时提供参考。
1. 波形数据选取及预处理
1.1 波形数据选取
本文地震波形数据截取的条件为：（1）震中距红山台为0°～50°（约0～5550 km）；（2）考虑到红山台测震数据的观测质量，选取2009～2021年的地震波形；（3）波形截取时，依据的地震目录为中国地震台网公布的正式观测报告（MS1.5～MS5.4）和IRIS（美国地震学联合研究会）公布的地震目录（MS5.5以上），从中选取红山台记录清晰且震源深度≤40 km的地震波形。（4）波形截取长度为震前1 min至震后30 min。最终截取地震记录共计4331条，地震分布如图1所示。
图 1 区域主要地震空间分布图（2009—2021）Figure 1. Spatial distribution of regional main earthquakes （2009—2021）1.2 数据分析方法原理及预处理
考虑到地震的原始记录叠加可能会有因波动方向不同而相互抵消的情况，本文参考Astiz等（1996）的研究方法，在数据处理过程中不对原始数据进行叠加，而是对原始数据进行一定的预处理，利用STA/LTA震相识别技术（图2），将计算得到的预处理数据STA/LTA值进行叠加（程序采用IRIS提供的修改后的globalstack计算脚本）。宽频带地震计共计观测三个分量（BHE、BHN、BHZ），数据处理过程中，对不同分量数据单独进行叠加计算，最后进行对比分析。计算之前，先将水平分量旋转成径向分量和横向分量。
图 2 STA/LTA震相识别技术原理图Figure 2. Schematic diagram of STA / LTA phase identification technology1.2.1 STA/LTA震相识别技术的基本原理
STA/LTA 震相识别法是由Stevenson于1976年提出并应用于地震波初至到时判断的一种方法（Stevenson，1976），很多学者在后续科研工作中对该方法进行了改进与和发展（刘晓明等，2017；乔红丽等，2018；赵哲，2018；孙印等，2018；李启成等，2019）。STA/LTA 法的基本原理为用Short Time Average（简称“STA”）和Long Time Average（简称“LTA”）的比值来反应地震序列幅值的变化。其中，STA反映了地震信号振幅的变化趋势，LTA反映了采集到的背景噪声信号振幅的变化趋势，如图2所示。
假设在记录到的地震波形上取1个测试点i，分别计算i点对应的长度为l1的短时窗振幅平均值（STA）和长度为l2的长时窗振幅平均值（LTA），通过计算测试点i所对应的STA与LTA的比值即可反映出地震信号强弱的变化，从而判断有无地震信号出现。
$$ LTA = \frac{1}{{{l_2}}}\sum\limits_{i - {i_2}}^i {\left| {A(i)} \right|} $$ (1) $$ STA = \frac{1}{{{l_1}}}\sum\limits_{i - {i_1}}^i {\left| {A(i)} \right|} $$ (2) 式中，
$A(i) $ 是测震数据采集器记录到的地震信号。预设阈值k用于判断有无有效地震发生，当发生1个有效的地震事件并被测震数据采集器采集到时，短时窗振幅平均值与长时窗振幅平均值的比值 STA/LTA 会出现非常明显的增加现象， STA/LTA 值大于判断阈值 k的时刻，定义为本次有效地震事件的地震波初至时间。1.2.2 数据预处理
预处理的过程主要包括：
（1）降采样。原始数据采样率为100 Hz，且截取的波形时长和数据量过大，为减少计算难度，在不影响计算结果的前提下，将采样率降为10 Hz。利用SAC软件中的decimate命令进行降采样，在降采样过程中，decimate命令会对数据自动做低通滤波，避免产生混叠现象。
（2）数据集存储。将波形记录按照震中距间隔0.2°作为1个数据集进行存储，共计250个数据集，个别震中距范围内缺少震例，如6.2°～6.4°、7.6°～7.8°、13.4°～13.6°、20°～20.2°，具体如图3所示。
图 3 距红山台0°～50°范围内各数据集地震波形数统计Figure 3. Statistics of seismic waveform number of each data set in the range of 0° to 50° from Hongshan station（3）滤波。因不同震相频率不同，为显示红山台记录到的所有震相，叠加前使用不同滤波器对各数据集内的波形数据进行高通或低通滤波，本文参考Astiz等（1996）在全球地震叠加并绘制全球各震相走时曲线时所用的参数组合，多次尝试选择了最优的滤波参数，如表1所示。
表 1 数据处理参数组合Table 1. Data processing parameter combination组合序号 滤波 STA/s LTA/s 1 高通0.5 Hz 1.0 9 2 高通0.167 Hz 2.0 20 3 低通0.1 Hz 3.0 30 4 低通0.033 Hz 4.5 45 （4）叠加。将数据集内所有波形数据计算得到的STA/LTA值进行叠加（图4）。
如图4所示，曲线1表示发生于北京时间2015年11月11日08:26:58，震中位于密克罗尼西亚加罗林群岛，距红山台43.3°的MS5.6地震的部分原始波形数据；曲线2表示对该原始曲线进行低通0.1 Hz滤波后的曲线；曲线3表示采用参数STA=3 s、LTA=30 s对滤波后的曲线进行计算得到的STA/LTA曲线；曲线4表示对43.2°～43.4°震中距范围内所有波形数据进行叠加后的结果。
2. 叠加计算
运用绘图脚本，将各数据集内STA/LTA叠加值绘制成震相走时图，绘制时间范围0～30 min，如图5～图8，分别为按照表1所示4种参数组合计算得到的叠加结果及可识别震相的IASP91模型理论走时曲线。。因各数据集内可用的地震波形数量差别较大，导致各震相走时曲线在不同震中距范围内色度和清晰度有所差别。由图3可知，在0°～4°、14.8°～15.6°、21.6°～23.4°震中距范围内波形数量较多，相应地，图5～图8在该震中距范围内震相走时曲线相较其它震中距颜色深，且清晰度更高。为了更清楚地展示红山台测震数据叠加后记录到的震相，本文在单独绘制垂向、径向、横向3个分量叠加图的基础上，同时绘制了可识别震相的IASP91速度模型计算理论走时曲线，震中距范围为0°～50°，震源深度为20 km，时间长度为30 min。
图 5 参数组合1叠加计算结果及IASP91模型理论走时曲线Figure 5. Superposition calculation results of parameter combination 1 and theoretical travel time curves of IASP91 model
图 6 参数组合2叠加计算结果及IASP91模型理论走时曲线Figure 6. Superposition calculation results of parameter combination 2 and theoretical travel time curves of IASP91 model
图 7 参数组合3叠加计算结果及IASP91模型理论走时曲线Figure 7. Superposition calculation results of parameter combination 3 and theoretical travel time curves of IASP91 model
图 8 参数组合4叠加计算结果及IASP91模型理论走时曲线Figure 8. Superposition calculation results of parameter combination 4 and theoretical travel time curves of IASP91 model由图5可知，在0°～50°震中距范围内，垂向、径向、横向分量均能清楚地显示出P波走时曲线，其中垂向记录相较于其他2个分量显示的更为明显；在0°～20°震中距范围内，3个分量均显示出R面波震相走时曲线，其中显示最清楚的是横向分量；在14°～50°震中距范围内，3个观测分量都能看出S波震相走时曲线痕迹，其中横向分量较垂向和径向分量更明显；在0°～14°震中距范围内，S波走时曲线显示不明显，分析原因，除此震中距范围波形数量较少之外，可能与STA/LTA震相识别技术针对初至波识别更好，对距离较近的后至波识别相对较弱有关。
由图6可知，相较于参数组合1计算结果，垂向、径向、横向均更清楚地显示出P波、S波、R波震相的走时曲线，且3个震相痕迹显示的震中距范围分别为0°～50°、14°～50°和0°～20°。此外，还记录到核面反射波震相ScS和PcS 2个震相的痕迹，震中距范围分别为34°～48°和41°～49°。
由图7可知，在0°～9° 震中距范围内，垂向、径向、横向分量P波痕迹显示不明显；S波走时曲线痕迹出现在14°～50°震中距范围内；横向分量R波走时曲线相对垂向和径向分量更为明显，震中距范围为0°～20°；在35°～48°震中距范围内显示出PP波和SS波走时曲线痕迹，且在SS震相走时上显示有面波走时痕迹；在35°～48°震中距范围内显示出R波走时痕迹。
由图8可知，在0°～10°震中距范围内，各震相显示不明显，在10°～50°震中距范围内，垂向分量P波走时曲线痕迹相对径向和横向分量更为清楚，3个测向分量均较明显的显示S波走时曲线，在36°～48°震中距范围内显示出SS波和PP波震相痕迹。在10°～48°震中距范围内，SS震相上方显示有面波走时痕迹，且呈明显频散现象。
3. 区域各震相波速计算
为了更好地显示区域地震主要震相走时曲线并测算其传播速度，选取红山台0°～20°震中距范围内数据，按照表1参数组合进行叠加计算，并绘制该震中距范围内各震相走时曲线图。限于篇幅，每个参数组合仅选取垂直分量叠加图和3个分量显示的震相综合图，如图9～图12。由于滤波参数不同，导致叠加后同一震相显示出的震中距范围不同，故在图中将未显示出的走时曲线用虚线表示。
图 9 参数组合1叠加计算结果（0°～20°）Figure 9. Superposition calculation results of parameter combination 1（0°～20°）
图 10 参数组合2叠加计算结果（0°～20°）Figure 10. Superposition calculation results of parameter combination 2（0°～20°）
图 11 参数组合3叠加计算结果（0°～20°）Figure 11. Superposition calculation results of parameter combination 3（0°～20°）由图9～图12可以看出，在震中距0°～15°范围内，区域地震各震相走时曲线基本成线性变化。综合各参数组合内选取的易于识别的典型地震震中距和走时，本文在震中距0°～15°范围内计算得到红山台区域地震P波传播速度约为7.5 km/s，S波传播速度约为4.2 km/s，R波传播速度约为3.5 km/s。此外，由图10（垂向）可知，在P波和S波走时曲线之间，震中距0°～15°范围内存在1个波速介于P波和S波之间的震相走时痕迹，该震相波速约为5.4 km/s。
图 12 参数组合4叠加计算结果（0°～20°）Figure 12. Superposition calculation results of parameter combination 4（0°～20°）4. 结论与讨论
本文首次运用STA/LTA震相识别技术，针对河北红山台多年（2009—2021）积累的地震记录进行叠加计算，成像了红山台记录到的区域地震各震相走时曲线。结果显示，震中距范围内波形数量越多，在该震中距范围内震相走时曲线颜色相对越深，清晰度越高。震中距0°～50°范围内，红山台共成像7种已知震相的走时曲线，分别为P、S、PP、SS、PcS、ScS以及R面波震相。选取不同的组合参数，叠加成像的震相种类有所不同，相同震相显示出的震中距范围和清晰度也不一致。例如仅在参数组合2中可以观察到PcS、ScS 2个震相的痕迹，SS、PP震相在参数组合3中比参数组合4中显示的更为清楚。
在震中距0°～15°范围内，P波、S波及R波走时曲线基本呈线性变化，综合各参数组合内选取的易于识别的典型地震震中距和走时，本文在震中距0°～15°范围内计算得到红山台区域地震P波传播速度约为7.5 km/s，S波传播速度约为4.2 km/s，R波传播速度约为3.5 km/s。此外介于P波和S波走时曲线之间，在震中距0°～15°范围内存在1个波速约为5.4 km/s的震相走时痕迹。
本工作对于提升红山台震中距≤1000 km的地震预警定位精度有积极的指导意义。但由于地震记录样本数量有限，震相走时曲线连续性和清晰度有所欠缺。下一步工作中，可选取河北地区其他宽频带地震记录进行联合计算，研究探讨各震相震中距记录范围及其波速。
致谢 感谢审稿专家提出的宝贵修改意见！
-
图 1 区域主要地震空间分布图（2009—2021）
Figure 1. Spatial distribution of regional main earthquakes （2009—2021）
图 2 STA/LTA震相识别技术原理图
Figure 2. Schematic diagram of STA / LTA phase identification technology
图 3 距红山台0°～50°范围内各数据集地震波形数统计
Figure 3. Statistics of seismic waveform number of each data set in the range of 0° to 50° from Hongshan station
图 5 参数组合1叠加计算结果及IASP91模型理论走时曲线
Figure 5. Superposition calculation results of parameter combination 1 and theoretical travel time curves of IASP91 model
图 6 参数组合2叠加计算结果及IASP91模型理论走时曲线
Figure 6. Superposition calculation results of parameter combination 2 and theoretical travel time curves of IASP91 model
图 7 参数组合3叠加计算结果及IASP91模型理论走时曲线
Figure 7. Superposition calculation results of parameter combination 3 and theoretical travel time curves of IASP91 model
图 8 参数组合4叠加计算结果及IASP91模型理论走时曲线
Figure 8. Superposition calculation results of parameter combination 4 and theoretical travel time curves of IASP91 model
图 9 参数组合1叠加计算结果（0°～20°）
Figure 9. Superposition calculation results of parameter combination 1（0°～20°）
图 10 参数组合2叠加计算结果（0°～20°）
Figure 10. Superposition calculation results of parameter combination 2（0°～20°）
图 11 参数组合3叠加计算结果（0°～20°）
Figure 11. Superposition calculation results of parameter combination 3（0°～20°）
图 12 参数组合4叠加计算结果（0°～20°）
Figure 12. Superposition calculation results of parameter combination 4（0°～20°）
表 1 数据处理参数组合
Table 1. Data processing parameter combination
组合序号 滤波 STA/s LTA/s 1 高通0.5 Hz 1.0 9 2 高通0.167 Hz 2.0 20 3 低通0.1 Hz 3.0 30 4 低通0.033 Hz 4.5 45 
下载: 导出CSV
-
冯太林, 张学工, 李衍达等, 2001. 折射波地震记录叠加成像方法研究. 地球物理学报, 44（1）: 129—134Feng T. L. , Zhang X. G. , Li Y. D. , et al. , 2001. Research on methodology of stack imaging of refractive seismic recording. Chinese Journal of Geophysics, 44(1): 129—134. (in Chinese) 蒋振武, 吴律, 1995. 井间地震反射波叠加成像的DLCDP法. 石油地球物理勘探, 30（4）: 495—504Jiang Z. W. , Wu L. , 1995. DLCDP method for the stacking and imaging of crossborehole seismic reflection waves. Oil Geophysical Prospecting, 30(4): 495—504. (in Chinese) 李启成, 何书耕, 2019. 用振幅变化长短时均值比实现地震预警中P波自动拾取. 地震工程学报, 41（1）: 138—146Li Q. C. , He S. G. , 2019. Automatic picking up of P-wave first arrival in earthquake early warning using STA/LTA method. China Earthquake Engineering Journal, 41(1): 138—146. (in Chinese) 李希元, 胡望水, 张楠等, 2020. 连续子波反射叠加合成地震记录方法. 大庆石油地质与开发, 39（2）: 133—138Li X. Y. , Hu W. S. , Zhang N. , et al. , 2020. Synthetic seismic recording method by the continuous wavelet-reflection superposition. Petroleum Geology & Oilfield Development in Daqing, 39(2): 133—138. (in Chinese) 刘晓明, 赵君杰, 王运敏等, 2017. 基于改进的STA/LTA方法的微地震P波自动拾取技术. 东北大学学报（自然科学版）, 38（5）: 740—745Liu X. M. , Zhao J. J. , Wang Y. M. , et al. , 2017. Automatic picking of microseismic events P-wave arrivals based on improved method of STA/LTA. Journal of Northeastern University (Natural Science), 38(5): 740—745. (in Chinese) 芦俊, 王赟, 季玉新等, 2018. 多分量地震数据的成像技术. 地球物理学报, 61（8）: 3499—3514Lu J. , Wang Y. , Ji Y. X. , et al. , 2018. Imaging techniques of multi-component seismic data. Chinese Journal of Geophysics, 61(8): 3499—3514. (in Chinese) 乔红丽, 张常在, 2018. 云计算环境下震前震源异常次声波自动识别方法. 地震工程学报, 40（6）: 1331—1336Qiao H. L. , Zhang C. Z. , 2018. Automatic identification of anomalous infrasonic waves prior to earthquake in cloud computing environment. China Earthquake Engineering Journal, 40(6): 1331—1336. (in Chinese) 秦满忠, 沈旭章, 张元生等, 2014. 利用兰州小孔径地震台阵资料叠加观测走时曲线. 地震学报, 36（1）: 59—69Qin M. Z. , Shen X. Z. , Zhang Y. S. , et al. , 2014. Observed travel-time curves by stacking records from Lanzhou small aperture seismic array. Acta Seismologica Sinica, 36(1): 59—69. (in Chinese) 孙印, 潘素珍, 刘明军, 2018. 天然地震识别与震相自动拾取技术进展. 中国地震, 34（4）: 606—620Sun Y. , Pan S. Z. , Liu M. J. , 2018. Seismic phase identification associated with the automatic pickup technology and its application. Earthquake Research in China, 34(4): 606—620. (in Chinese) 岳玉波, 张建磊, 张超阳等, 2020. 基于时变数据映射的地震叠加成像方法. 石油地球物理勘探, 55（2）: 331—340Yue Y. B. , Zhang J. L. , Zhang C. Y. , et al. , 2020. Seismic data stacking based on time-varying mapping. Oil Geophysical Prospecting, 55(2): 331—340. (in Chinese) 赵哲, 2018. 基于STA/LTA法的地震波初至时间提取方法. 中国科技信息, （9）: 83—84. Astiz L. , Earle P. , Shearer P. , 1996. Global stacking of broadband seismograms. Seismological Research Letters, 67(4): 8—18. doi: 10.1785/gssrl.67.4.8 Earle P. S. , 1999. Polarization of the earth’s teleseismic wavefield. Geophysical Journal International, 139(1): 1—8. doi: 10.1046/j.1365-246X.1999.00908.x Okal E. A., Cansi Y., 1998. Detection of PKJKP at intermediate periods by progressive multi-channel correlation. Earth and Planetary Science Letters, 164(1—2): 23—30. Shearer P. M. , 1991. Imaging global body wave phases by stacking long-period seismograms. Journal of Geophysical Research: Solid Earth, 96(B12): 20353—20364. doi: 10.1029/91JB00421 Shearer P. M. , Rychert C. A. , Liu Q. Y. , 2011. On the visibility of the inner-core shear wave phase PKJKP at long periods. Geophysical Journal International, 185(3): 1379—1383. doi: 10.1111/j.1365-246X.2011.05011.x Stevenson P. R. , 1976. Microearthquakes at Flathead Lake, Montana: a study using automatic earthquake processing. Bulletin of the Seismological Society of America, 66(1): 61—80. doi: 10.1785/BSSA0660010061 Walck M. C. , Clayton R. W. , 1984. Analysis of upper mantle structure using wave field continuation of P waves. Bulletin of the Seismological Society of America, 74(5): 1703—1719. doi: 10.1785/BSSA0740051703 -
下载:
点击查看大图
计量
- 文章访问数: 108
- HTML全文浏览量: 49
- PDF下载量: 6
- 被引次数: 0
