Spatial domain division based on superposition pattern of seismic waves and spatial variation of motion：Case of SV wave incidence
-
摘要: 针对地震动空间差异问题，以半无限空间内平面SV波入射为例，采用波动理论和叠加原理相结合的方法研究地震动的空间特征。首先，分析入射波和反射波在空间内形成质点运动的叠加模式差异，并依此对空间域进行划分；其次，针对不同泊松比和SV波入射角情形研究叠加区与分离区的分界线控制情况；最后，对比分析不同空间域内的质点运动在峰值、持时等方面的特征。研究结果显示：当SV波入射时，叠加区与分离区的分界线通常由z3 （反射P波与反射SV波的分离线）控制；同时存在两种特殊情况，当SV波垂直入射时，分界线由z1 （入射SV波与反射SV波的分离线）控制，当反射SV波幅值为零时，分界线由z2 （入射SV波与反射P波的分离线）控制。在入射波和反射波的影响下，质点运动时程的形状具有水平不变性。三波贡献时段只出现在深度小于z1的质点的运动时程中，且持时随着深度的增加线性减少；双波贡献时段出现在位于叠加区内的质点的运动时程中，持时沿深度先增加后减少；单波贡献时段随着深度的增加而逐渐加长，在分离区达到最大值。质点运动总持时随深度逐渐增加，在分离线z2和z3深度处存在两个拐点。在质点运动峰值方面，靠近自由面的叠加区质点运动峰值变化较大，深度较大的叠加区和分离区的质点运动峰值一般不变。Abstract: For the investigation of spatial variation of seismic motion, the wave theory and superposition principle are used to illustrate the spatial variation of motions caused by seismic waves. And plane SV wave’s incidence in semi-infinite space is used as a case in this paper. Firstly the difference in superposition pattern of waves in different position is investigated and the space region is divided into different domains according to superposition pattern. Then the border line between superposition domain and separation domain with different Poisson’s ratio and SV wave’s incident angle is discussed. Finally the difference in peak value and duration time of particle motion in different domains is investigated. The result shows that, the border line is usually controlled by z3 (separation line generated by reflected P wave and reflected SV wave). Particularly, there exist two special cases, one is the vertical incidence of SV wave in which the border line is controlled by z1 (separation line generated by incident SV wave and reflected SV wave), the other is the zero amplitude of reflected SV wave in which the border line is controlled by z2 (separation line generated by incident SV wave and reflected P wave). The time history of different particle motion in one depth has the same shape with the influence of incident and reflected waves. In time history of particle motion, the period contributed by three waves only exists in depth lower than z1 and the duration time decreases with depth. The period contributed by two waves exists in all superposition domain and the duration time increases firstly and then decreases with depth. The period contributed by one single wave increases with depth and reaches its maximum in separation domain. The increase of particle motion duration has two slope change in depths of z2 and z3. The peak value of particle motion changes dramatically in superposition domain near free surface and remains constant in deep superposition domain and separation domain.
-
引言
地震震相信息是地震学研究的重要基础数据，广泛应用于地震定位、震源机制分析、走时层析成像等研究，震相自动识别也是地震预警技术的重要基础. 随着地震观测技术的不断发展，利用密集地震台网和大规模流动地震台阵开展地震监测研究受到广泛重视，与此同时观测产生的地震波形数据急剧增加. 为了获取观测区内的微小地震，需要对大量的连续地震观测数据进行仔细分析，而传统方法是依靠人工一帧一帧地浏览地震波形识别地震，并逐一分析读取震相，工作量巨大，整个过程单调枯燥且效率极低，因此地震震相自动拾取和地震自动检测技术的发展，一直备受关注.
震相拾取方法众多，每种方法均有其优势和局限性，而有效性和计算效率是衡量震相拾取方法优劣的两个重要指标，而且这两个指标会受到诸如数据质量、算法简洁性等因素的影响（Chen，Holland，2016）. 应用于实时监测的震相自动识别方法还要求在保证一定效果的情况下，注重算法的计算效率. 与之区别的常规测震分析则需注重有效性和精度，常用于微震识别和地震定位等分析.
震相到时拾取算法可分为在时域、频域中进行两类，包括能量瞬时方法（Allen，1978，1982；Withers et al，1998 ；高淑芳等，2008）、偏振分析（Flinn，1965；Vidale，1986；Bai，2000；Jurkevics，2002）、自回归赤池信息准则（Akaike information criterion，简写为AIC）方法（Akaike，1974；Kitagawa，Akaike，1978；Maeda，1985； Leonard，Kennett，1999；Leonarda，2000；Zhang et al，2003 ；赵大鹏等，2012）、人工神经元（Wang，Teng，1997；Gentili，Michelini，2006；Dai，MacBeth，2007）、高阶统计（Saragiotis et al，2002 ；Panagiotakis et al，2008 ； Baillard et al，2014 ；Hibert et al，2014 ）、小波变换（Anant，Dowla，1997；刘希强等，1998；刘代志等，2005；Akansu et al，2010 ；Vassallo et al，2012 ；Bogiatzis，Ishii，2015）、奇异值分解（Rosenberger，2010；Kurzon et al，2014 ）等方法. 刘希强等（2009）的研究表明，根据信号与噪声的单一差异特征而提出的方法均具有其应用的局限性，特别是低噪声水平条件下震相自动识别的准确率仍不理想，而综合多个信号与噪声的差异特征有利于提高震相自动识别的准确性.
在震相拾取时，利用多台信息能够降低噪声的干扰，提高震相识别的可靠性. 马强等（2013）将德洛内（Delaunay）三角用于地震预警的P波震相拾取，其基本原理是利用短长时窗均值之比（short term average to long term average，简写为STA/LTA）方法得到初步到时，再用AIC方法得到P波的精确到时，通过德洛内三角判别该到时是否为干扰信号. 该方法现已应用于福建省地震预警试验系统中，在剔除干扰信号方面具有良好的效果，但并未将各台站触发到时归类到相同的地震事件中.
针对以上问题，本文拟对拾取的到时信息，利用台间距与台站相邻关系作为约束来减少噪声干扰，并在快速地震定位的基础上，根据台站理论到时选取合理的时间窗，采用STA/LTA方法、偏振分析与基于方差的AIC （var-AIC）的组合方法进行更精细的分析，以期通过该方法获取更可靠的P波震相到时信息.
1. 方法
目前震相拾取主要分为3步：① 数据预处理. 主要对波形进行滤波或时频转换，包括小波变换、奇异值分解、傅里叶变换等；② 地震信号检测. 以STA/LTA方法为主，但由于STA/LTA方法具有一定的延迟性，因此只能得到地震信号所在的大致到时；③ 拾取震相到时. 主要通过AIC方法、偏振分析方法以及人工神经元等方法拾取P波或S波到时.
以上步骤首先基于STA/LTA的检测到时结果，但是地震观测数据经常会受到台站周边人文环境、气候等因素的干扰，例如车辆、工厂、动物活动以及大风、雷电等产生的干扰，由此检测到虚假震相；但是由于这些干扰的传播距离较小，往往仅能被邻近台站观测到，利用相同地震事件中震相到时具有的同源性，可滤除大部分这类噪声干扰. 利用震源位置信息获得预测的P波到时选取窗口，再使用偏振分析和AIC函数等方法进行更精细的震相到时拾取. 同时，对于一些STA/LTA方法获得的信噪比较低的台站，在P波理论到时窗口内进一步分析，从而判定是否能获得有效的到时信息.
1.1 地震事件检测
首先，本文使用STA/LTA方法（Allen，1978，1982）快速获得低精度的地震波初至到时. 计算短时窗N2与长时窗N1内的特征函数CF的平均值比值STA（t）/LTA（t），其中：
${{CF}}(t) = {Z^2}(t),$
(1) ${{LTA}}(i) = {{LTA}}(i - 1)+ \frac{{{{CF}}(i) - {{CF}}(i - {N_1})}}{{{N_1}}},$
(2) ${{STA}}(i) = {{STA}}(i - 1) + \frac{{{{CF}}(i) - {{CF}}(i - {N_2})}}{{{N_2}}},$
(3) 式中，Z表示地震记录的垂直分量. 当STA（t）/LTA（t）超过阈值C1时，该时刻tP记为地震波初至到时；当STA（t）/LTA（tP）低于阈值C2时，该时刻记为信号恢复时刻te，此时获得的tP到时相较P波实际到时有一定滞后. Maeda（1985）提出基于方差的AIC方法，利用该方法来计算AIC函数的最小值，能得到更精确的P波到时，其中AIC函数为
${ AIC}(t) = t\lg [{\mathop{ var}} (D[1,t])]+ (N - t - 1)\lg [{\mathop{ var}} (D[t + 1,N])],$
(4) 式中，N为信号长度，D为截取信号. var-AIC方法需要预先设置好窗长才能计算，在STA/LTA初步拾取后，本文采用一个长度为N3的小窗对低精度的初至到时进行适当校正.
将各台站记录到的触发到时按时间先后顺序排列，并将震相到时和台站信息按照地震事件进行分组归类. 与台站附近的局部噪声干扰不同，一个地震事件若能被多个地震台站观测到，则相邻台站记录的初至震相到时应该具有较好的相关性. 根据地震波传播的特点可知，相邻台站间震相的最大到时差应小于台间距除以波速并考虑一定的误差量，而台站附近的局部噪声干扰具有随机性，相关性较弱. 根据这一特征，通过引入台网德洛内三角划分，并考虑台间距的影响，即可将震相到时差大于约束条件的台间距排除在外. 对于相继触发的两个台站1和2，根据以下3个条件判断二者触发的到时信号是否可能属于同一事件：① 台站2与台站1记录的到时差是否小于台间距除以波速vP；② 台站2与台站1的台间距是否小于台间距阈值C3；③ 台站2与台站1是否在同一德洛内三角内. 满足条件①且满足条件②或③的两台站所触发的到时认为在同一事件内，其它情况则排除在同一事件外.
通过以上步骤，将时空上具有相关性的相邻到时记录暂时归类为同一事件. 但考虑到后续触发的台站，其震中距不断增大，同一事件触发的台站虽然在时间上相邻，但是其既可能在空间上不相邻，也可能不在台间距阈值内，此时可能会造成同一事件记录的重复拾取. 为避免这种情况，将各个相邻并且互相有所关联的事件合并：两个相邻的事件A和B，取事件B触发的前5个台站与事件A触发的最后5个台站，分别计算前一组中台站与后一组每个台站的距离和到时差（事件A和B不足5台时均计算全部），若有两个台站满足上述3个条件，则将事件A与事件B合并为同一个事件. 最后提取其中触发台站数大于一定数值的事件记为地震事件，其它事件则作为干扰噪声予以剔除.
滤除无相关性的噪声并获得地震事件对应的各台站到时后，利用偏振特性计算各台站的线性偏振度P，从而获取S波震相到时. 首先选取一个长度为N的滑动视窗，计算三分向协方差矩阵的3个特征值λ1，λ2和λ3 （λ1>λ2>λ3），则可得到线性偏振度P为
$P(t) = 1 - \frac{{{\lambda _2} + {\lambda _3}}}{{{\lambda _1}}},$
(5) This page contains the following errors:
error on line 1 at column 1: Start tag expected, '<' not foundBelow is a rendering of the page up to the first error.
$D{_{{ b}t}} = \frac{{{\lambda _{t + }}}}{{{\lambda _{t - }}}};$
(6) 同样取一短窗长Ls，计算t短窗的极化度Dst. 拾取震相的特征函数Pol为
${ Pol}(t) = D{_{{ b}t}}D{_{{ s}t}}.$
(7) 当Pol取最大值时，t点即为震相初至. 由于P波和S波均具有偏振特性，实际波形中P波的尾波会影响到S波的偏振性，因此若Pol函数的值低于阈值C4，则认为S波到时的计算结果不可靠，不记录该S波的到时.
1.2 地震定位判断
首先计算某事件所有触发台站的平均坐标，取平均坐标为中心的正方形区域，按照经纬度和深度进行三维网格化，并确定网格步长.
This page contains the following errors:
error on line 1 at column 1: Start tag expected, '<' not foundBelow is a rendering of the page up to the first error.
${t_0} = \sum\limits_{i = 1}^n {\frac{{t_2^i - t_1^i}}{n}} ,$
(8) ${t_3} = \sqrt {\sum\limits_{i = 1}^n {\frac{{{{(t_2^i - t_1^i - {t_0})}^2}}}{n}} } ,$
(9) 式中n为事件所包含的震相到时数. 取所有网格的标准差t3的最小值，其所处位置为震源位置的可能性最大，此时各台站的残差为
$t_4^i = \left| {t_2^i - t_1^i - {t_0}} \right|.$
(10) This page contains the following errors:
error on line 1 at column 1: Start tag expected, '<' not foundBelow is a rendering of the page up to the first error.
1.3 震相到时拾取的组合函数方法
在时窗内拾取震相到时之前，须判断该时窗内是否存在比较清晰的震相. 为此，本文设计了一个简单的震相存在判别方法：首先查看该台站是否为地震检测过程中触发的台站，若是，则认为该台站记录到震相并进行下一步计算，否则计算时窗内的最大偏振度以及窗内长短时触发位置处的信噪比，若二者之一大于各自的阈值C6或C7，则判定该时窗内有震相.
考虑到AIC函数和STA/LTA拾取到时的结果随信噪比的变化较为强烈，当信噪比较低时，AIC函数或STA/LTA可能表现不稳定，从而使得二者拾取到时的差值较大；而偏振分析方法虽然受到信噪比的影响程度相对较小，但其在高信噪比时拾取精度不如前两者高，因此我们设计了一种新的组合方法，提高震相拾取的可靠性，即利用AIC函数与STA/LTA两方法拾取到时的差值作为信噪比高低的判别依据，当差值较小即信噪比较高时，选择较准确的AIC函数的极小值作为可靠结果，信噪比较低时设计新的拾取函数
${ Pick}\left( {{i}} \right) = \left[ {1 - {{AIC}}\left( {{i}} \right)} \right] { Pol}\left( {{i}} \right),$
(11) 该函数的最大值对应P波到时. Pick函数利用归一化后的AIC函数压制了P波到达后其它信号对Pol函数造成的影响. 通过上述方法在时窗内进一步分析P波、S波到时，主要步骤为：① 先在预测P波、S波到时前后取窗长N4，根据震相存在判别方法判断该时窗内是否具有震相，若有则在该时窗内同时使用STA/LTA方法和var-AIC方法检测到时；② 若AIC函数与STA/LTA检测到的到时差小于C8，则将该到时记为P波到时；若大于C8或STA/LTA未检测到任何到时信号，则对时窗内的波形数据使用偏振分析方法并对AIC函数归一化，取拾取函数Pick的最大值为P波到时；③ 以P波到时与AK135计算得到的S波到时的中间值为时窗的起始点，S波到时向后N4范围时窗结束，在该时窗内计算Pol函数的最大值，若最大值大于C4，则记为S波到时.
2. 计算实例
2013年1月，房立华等（2014）在安宁河断裂带附近布设了由32个流动地震台站组成的密集流动地震台阵，即西昌地震观测台阵，开展针对断层带的地震观测研究. 台阵的台间距介于17.7—60.5 km，平均台间距为27.8 km. 地震观测仪的频带为50 Hz—1 200 s，采用GPS授时，观测使用的采样率为100 Hz.
本文以西昌地震台阵记录到的2013年1月28日21时4分ML1.2地震为例，此时共有28个台站正常运行. 首先使用德洛内三角剖分算法将台站划分为若干三角区（图1），接着使用本文提出的方法对该地震事件进行震相到时自动拾取试验，采用的相关参数列于表1. 首先对波形数据进行1—30 Hz的带通滤波处理，然后计算STA/LTA，并拾取震相到时. 在STA/LTA的计算过程中，长时窗N1设为30 s，短时窗N2设为0.05 s，触发阈值C1设为10，考虑到地震波到来后短时间内无法恢复至触发前的噪声水平，因此STA/LTA的恢复阈值C2设为1.5. 在利用STA/LTA方法拾取震相到时的基础上，进一步通过var-AIC和偏振分析方法拾取P波和S波到时. 图2给出了针对XC03台站的地震波形记录获得的分析结果. 由STA/LTA，var-AIC，偏振分析方法以及组合函数方法计算得到的P波初至结果分别为7.52，7.49，7.48和7.48 s，人机交互拾取结果为7.44 s. 与人工方法拾取的到时相比，var-AIC，偏振分析方法和组合函数方法拾取的到时分别相差0.05 s和0.04 s，STA/LTA方法的拾取到时则相差0.08 s.
![]()
图 1 西昌地震台阵的德洛内三角剖分图正方形为快速定位时网格搜索区域，星形为地震实例发生位置，三角形为台站位置Figure 1. Delaunay triangulation of Xi’chang seismic arrayThe square represents the search area，the star stands for the hypocentral location，and triangles for stations表 1 本文方法采用的相关参数Table 1. Related parameters used in the method proposed in this paper参数名 数值 单位 意义 C1 10 长短时方法获得到时的触发阈值 C2 1.5 长短时方法信号恢复时的触发阈值 C3 30 km 多台检测事件台间距的上限 C4 20 偏振分析Pol函数的触发阈值 C5 5 s 定位时的残差阈值 C6 10 震相存在判别的信噪比阈值 C7 1 000 震相存在判别的Pol函数阈值 C8 0.5 s AIC与STA/LTA结果接近的到时差范围 N1 30 s 长时平均值计算的滑动窗长 N2 0.5 s 短时平均值计算的滑动窗长 N3 1 s 为纠正STA/LTA拾取到时的AIC函数计算窗长 N4 4 s 代入AK135后拾取P波和S波到时的窗长 Lb 0.5 s 通过Pol函数计算S波到时设置的长窗长 Ls 0.2 s 通过Pol函数计算S波到时设置的短窗长 vP 5.8 km/s 检测地震事件时所取的P波速度 ![]()
图 2 西昌地震台阵记录的2013年1月28日21时4分ML1.2地震（a） 垂直分量波形；（b） 短长时窗均值之比；（c） 归一化的var-AIC函数；（d） Pol函数拾取P波；（e） Pol函数拾取S波；（f） 改良后的组合函数Figure 2. Event of ML1.2 recorded by Xi’chang seimic array at 21:04 on January 28，2013（a） Vertical component of waveform；（b） STA/LTA；（c） Normalized var-AIC；（d） P-wave picking by Pol-function；（e） S-wave picking by Pol-function；（f） Improved compound function图3给出了2013年1月28日21时4分ML1.2地震在15个观测台站的垂直向地震记录. 由于噪声干扰等因素的影响，使用STA/LTA方法与var-AIC方法初步拾取的震相到时出现了许多与地震信号无关的虚假震相到时（图3a）. 利用地震信号在多台记录中满足同源性的特点，采用德洛内三角和台间距与传播时间关系的约束，在很大程度上消除了台站附近局部噪声干扰的影响，很好地保留了基本满足同源传播特征的地震震相到时（图3b），成功地检测出地震事件. 值得注意的是，本文设置判断到时相关性的台间距阈值C3略大于平均台间距，主要是考虑到台间距阈值C3的主要作用是约束远处台站干扰信号的影响，而平均台间距能够作为台站远近的一项评价指标.
在初步完成地震事件检测后，以同一地震事件所有触发台站的平均坐标为中心，南北向、东西向的长度均取为1°，深度取0—60 km，以0.05°×0.05°×5 km的三维网格进行快速地震定位. 图1中正方形区域为该地震事件的搜索区域，图4给出了该事件的定位结果. 定位结果显示，与人工分析的正式定位结果相比，震中仅相差5.39 km. 在地震定位过程中，对于一些走时残差过大的震相数据进行适当剔除，可以进一步消除各种干扰的影响.
![]()
图 3 地震观测记录及拾取的震相到时图（a） STA/LTA方法拾取的各种可能的P波震相到时；（b） 利用本文地震检测方法分析处理后获得的P波震相到时；（c） 用AK135模型计算的理论到时，其中紫线为P波预测到时，黄线为S波预测到时，蓝线为地震检测方法得到的P波震相到时，红线和绿线分别为取窗拾取的P波到时和S波到时，其中剔除了Pol函数小于C4的S波到时Figure 3. Waveform records and arrival picking in different steps（a） Probable P-waves picking by STA/LTA；（b） P-waves processing by the earthquake detection method in this paper；（c） Arrival time calculated by AK135 model，where purple lines denote arrival time of P-wave，yellow lines denote arrival time of S-wave，blue lines denote arrival time of P-wave processing by earthquake detection method，red and green lines denote P-wave and S-wave picking by setting windows in waveforms，S-wave of which Pol-function below C4 is deleted根据地震定位结果，本文采用AK135模型计算各个台站的P波和S波预测到时（如图3c中的紫色和黄色竖线所示），并在预测到时附近设定适当宽度的时间窗，进一步拾取更可靠的震相到时. 由于采用的速度模型与真实模型存在差异、震相到时拾取存在误差以及定位结果本身也存在误差等多种因素的影响，理论到时与实际震相的到时之间存在一定的差异. 在实际资料的处理中，以P波预测到时为中心前后取4 s时窗，拾取P波到时. 取P波到时与S波预测到时的中间值与S波预测到时后4 s作为时窗计算Pol函数的最大值，即为S波到时. 信噪比较高时（图2），AIC函数与STA/LTA之差很小，此时较准确的AIC函数的极小值可作为可靠结果；若信噪较低或受到干扰信号影响，AIC函数或短长时之比可能表现得并不稳定，二者差值也会较大，此时可通过本文拾取函数拾取P波到时，如图5所示. 图3c为AK135模型计算的预测到时与取窗分析到时对比图，可见相较于预测到时，取窗分析的结果已有很大提高.
![]()
图 4 实际震源位置与算法计算结果的对比图中数字表示事件触发顺序Figure 4. Comparison of hypocentral location with calculation resultThe numbers represent the triggered order of events![]()
图 5 各震相拾取方法精度对比（a） 波形计算窗口；（b） 垂直分量波形；（c） 短长时比值；（d） var-AIC函数；（e） Pol函数；（f） 改良后的组合函数Figure 5. Accuracy comparison of different phase picking methods（a） Calculation window；（b） Vertical component of waveform；（c） STA/LTA；（d） var-AIC；（e） Pol-function；（f） Improved compound function3. 实际资料处理试验与讨论
为验证本文方法的有效性，本文选用2013年西昌地震观测台阵1—8月记录到的304个地震事件，震级分布于ML−0.5—2.2，震源深度最大为37 km. 这些地震是通过人工一帧一帧地浏览观测波形仔细分析后得到的，对许多微弱的地震信号均进行了分析. 根据人工分析获得的地震目录，将每一个地震截取为180 s的单事件波形数据，利用本文方法进行震相自动拾取，并将结果与人工拾取结果进行对比，以确认其有效性.
首先采用台间距阈值约束方法、德洛内三角判别的方法和本文方法对304个地震事件进行检测，并以人工拾取结果作为参照. 表2列出了3种检测方法的结果对比，可以看出，使用本文提出的改进方法检测，遗漏的地震事件明显少于其它两种方法，误拾取的地震事件数与德洛内三角判别检测方法相近.
由于待检测的地震大多震级较小，对于平均台间距为28 km的流动地震台阵而言，一些微小地震事件的地震波形往往只能被少量台站清晰地记录到. 图6给出了检测到的地震数与地震震级及设定的震相触发最少台站数的关系，其中触发震相要满足改进方法的相关约束条件. 可以看出：设定的最少触发台站数越小，能检测到的地震数越大；当ML>0.9时，采用最少3台触发能检测到所有的地震；当ML<0.6时，设定的最小触发台站数对遗漏的地震事件数具有明显的影响；对于ML>0.6地震，设定的最少触发台站数为3，4和5时，地震的正确识别率分别为95%，88%和81%；对于ML>0.3地震，相应的正确识别率则分别为90%，77%和69%. 为了检验震源深度可能对地震检测效果的影响，本文将不同深度的地震检测遗漏情况进行统计，结果如图7所示，可以看出不同深度的遗漏事件占比并无显著的相关性.
检测完地震事件后，为尽可能地保证定位效果，选取触发4台或4台以上的地震事件用于定位，共包含186个事件. 其中，有2个事件在定位过程中因剔除较大残差的到时而导致触发台站少于4个而放弃，最终共有184个事件用于计算. 地震定位采用AK135模型作为理论走时计算的速度模型. 由于本文在通过定位结果获取的预测到时前后4 s的时窗内分析实际到时，因此定位误差小于20 km足以满足进一步分析的要求，将定位结果与人工定位结果作比较，共有80%的地震事件定位误差小于20 km，89%的ML>1.0地震事件的定位误差小于20 km.
根据地震定位结果，采用AK135速度模型计算各台站的预测到时，并以预测到时为中心设定时窗，然后根据震相存在判别方法判断各台站是否存在震相，若存在则拾取该台站时窗内的震相. 为验证震相存在判别方法的准确性，用该方法判别184个地震事件中存在震相的台站数并与地震检测过程以及人工拾取结果进行对比分析，结果如图8所示，可以看出：该方法判别有震相记录的台站与人工拾取到震相记录的台站高度重合，人工拾取震相的85%能够被检测到，而之前不考虑理论走时的震相检测方法仅有77%的震相能被检测出来；对于ML>0.6地震，两种方法的检测能力分别是人工拾取震相的87%和79%；对于ML>1.5地震，分别能检测到89%和73%的震相. 基于理论到时的震相识别，更容易检测到一些人工识别遗漏的震相.图9给出了新方法震相检测情况与人工震相识别结果对比的两个实例. 可以看出，由于受多种因素的影响，人工识别既可能遗漏掉一些微弱的震相，也可能遗漏掉个别清晰的震相，而这些震相有可能被新的自动识别方法很好地识别出来（图9a）；自动识别方法也可能受短时强噪声干扰的影响错误地将噪声干扰识别为震相（图9b），需要通过其它方法尽可能排除这些干扰.
表 2 3种地震检测方法对304个地震事件的检测效果对比Table 2. Comparison of detection results for 304 events by three event detection methods台站数 拾取正确事件数 拾取错误事件数 拾取重复事件数 遗漏事件数 拾取事件数 本文
方法台间距 德洛内三
角判别本文
方法台间距 德洛内三
角判别本文
方法台间距 德洛内三
角判别本文
方法台间距 德洛内三
角判别本文
方法台间距 德洛内三
角判别3 239 186 214 19 3 21 32 19 29 65 118 90 290 208 264 4 186 135 167 8 0 7 12 5 14 118 169 137 206 140 188 5 149 107 135 1 0 2 10 3 10 155 197 169 160 110 147 ![]()
图 6 触发台站数不同时所拾取事件数与震级的关系Figure 6. Relationship between the amount of events with magnitude with different amount of triggered station![]()
图 7 遗漏事件占比与深度关系图Figure 7. Relationship between the proportion of omitted events with focal depth![]()
图 9 本文方法与人工拾取结果对比图横坐标为相对预测到时的时间，纵坐标已对齐预测到时。红线为本文方法拾取到的P波震相，蓝线为人工拾取P波震相Figure 9. Comparison of picking results by manual and the method proposed in this paperHorizontal ordinate is relative time. Vertical ordinate has been lined up with forecasted arrival time. Red lines represent P-waves picking by the method in this paper，and blue lines represent P-waves picking by manual基于地震定位后的震相预测到时设定的时窗，通过震相检测判别后，分别采用var-AIC方法、偏振分析方法和本文组合函数方法进行震相自动拾取，并将结果与人工的拾取结果进行对比分析，拾取效果如图10所示，可见：与人工拾取的震相到时相比，误差在0.2 s内的震相数量，本文方法要优于var-AIC方法和偏振分析方法；在误差大于0.2 s且信噪比低于4 dB时，本文方法优于var-AIC方法但稍弱于偏振分析方法.
![]()
图 10 不同方法震相拾取率与信噪比SNR的关系图（a） 拾取误差小于0.1 s；（b） 拾取误差小于0.2 s；（c） 拾取误差小于0.3 s；（d） 拾取误差小于0.4 sFigure 10. Relationship between phase picking rate with SNR by different methods（a） Deviation below 0.1 s；（b） Deviation below 0.2 s；（c） Deviation below 0.3 s；（d） Deviation below 0.4 s4. 讨论与结论
本文基于var-AIC和偏振分析方法的特点，提出了一种改进的震相拾取方法. 实际资料处理试验揭示，该方法总体上明显优于原有的两种方法，只是在信噪比低于4 dB左右时稍弱于偏振分析方法.
采用STA/LTA方法、var-AIC以及偏振分析方法拾取可能的震相到时，利用地震具有同源性的特征，通过台站的德洛内三角形划分、与台间距有关的最大走时约束以及地震定位约束等，可以较好地排除噪声干扰的影响，提高地震事件检测的可靠性. 西昌台阵2013年1—8月的实际资料检验表明，当作为地震事件判别的最小触发台站数选为4时，本文提出的方法对于ML>0.5地震的识别率为82%，ML>0.8地震的识别率为89%，ML>1.0地震的识别率为96%，比单独使用台间距阈值约束方法或使用德洛内三角判别方法得到的地震识别率分别提高21%和3%—4%.
以地震定位结果为基础，在计算的震相理论到时前后设定合理的震相搜索时间窗，可以从观测资料中拾取更多的震相，提高震相到时拾取的可靠性. 本研究仅采用简单的AK135模型进行地震定位，对于观测震相数量较少、台站分布较差的地震，定位精度可能不高，会在一定程度上影响最终震相拾取的可靠性. 在今后的研究中，可以结合观测区域的实际情况，建立合适的区域速度模型或三维速度模型开展地震定位和地震震相到时计算，有可能获得更加可靠的地震检测和震相拾取结果.
-
![]()
图 1 地震波叠加模式与质点运动空间差异示意图
（a） 问题说明；（b） A点和B点的水平向运动
Figure 1. Superposition pattern of seismic waves and spatial variation of particle motion at different point
（a） Description of problem；（b） Horizontal motion at points A and B
![]()
图 3 不同泊松比时叠加区范围系数a1 （a），a2 （b）和a3 （c）随平面SV波入射角的变化
Figure 3. Variation of superposition range coefficient a1 （a），a2 （b） and a3 （c） with incident angle of SV wave
![]()
图 4 泊松比为0.15 （a），0.25 （b）和0.35 （c）时叠加区范围系数对比
Figure 4. Comparison of superposition range coefficient with Poisson’s ratio of 0.15 （a），0.25 （b） and 0.35 （c）
![]()
图 5 分界线由Si与Pr分离线控制的对应情况
Figure 5. Border controlled by separation line generated by Si and Pr
![]()
图 7 不同深度h下特征位置质点的水平向运动时程
Figure 7. Time history of horizontal particle motions in different depth h
![]()
图 8 不同地震波叠加类型下持时沿深度分布
Figure 8. Duration of different superposition pattern of particle motion in different depth
![]()
图 10 质点运动水平向（a）、垂直向（b）峰值沿深度变化
Figure 10. Peak displacement of particle motion in horizontal （a） and vertical （b） displacement along depth
表 1 不同深度质点地震波贡献时段
Table 1 Contribution time interval of waves on particle motion at different depth
质点 深度/km 单波贡献时段/s 双波贡献时段/s 三波贡献时段/s Si Pr Sr （Si，Pr） （Pr，Sr） （Si，Pr，Sr） D1 0 0—2.00 D2 1.0 0—0.89 2.89—3.17 0.89—1.17 2.00—2.89 1.17—2.00 D3 1.8 0—1.61 2.00—2.11 3.61—4.11 1.61—2.00 2.11—3.61 D4 3.0 0—2.00 2.68—3.52 4.68—5.52 3.52—4.68 D5 9.0 0—2.00 8.05—10.05 10.57—12.57 
下载: 导出CSV
-
陈厚群,徐泽平,李敏. 2008. 汶川大地震和大坝抗震安全[J]. 水利学报,39(10):1158–1167. doi: 10.3321/j.issn:0559-9350.2008.10.002 Chen H Q,Xu Z P,Li M. 2008. Wenchuan earthquake and seismic safety of large dams[J]. Journal of Hydraulic Engineering,39(10):1158–1167 (in Chinese).
贺春晖,王进廷,张楚汉. 2017. 基于震源—河谷波场数值模拟的坝址地震动参数确定方法[J]. 地球物理学报,60(2):585–592. He C H,Wang J T,Zhang C H. 2017. Determination of seismic parameters for dam sites by numerical simulation of the rupture-canyon wave field[J]. Chinese Journal of Geophysics,60(2):585–592 (in Chinese).
何卫平,周宜红,赵春菊. 2019. 平面P波入射半无限介质空间域划分及运动偏振特征[J]. 岩石力学与工程学报,38(2):321–331. He W P,Zhou Y H,Zhao C J. 2019. Spatial division of semi-infinite space and polarization characteristics of particle motion for incident plane primary wave[J]. Chinese Journal of Rock Mechanics and Engineering,38(2):321–331 (in Chinese).
金星. 1993. 地震动场研究状况简介[J]. 世界地震工程,(3):19–22,52. Jin X. 1993. Introduction of research on earthquake motion field[J]. World Earthquake Engineering,(3):19–22,52 (in Chinese).
金星,廖振鹏. 1994. 地震动随机场的物理模拟[J]. 地震工程与工程振动,14(3):11–19. Jin X,Liao Z P. 1994. Physical simulation on the random field of seismic motion[J]. Earthquake Engineering and Engineering Vibration,14(3):11–19 (in Chinese).
林皋. 2017. 地下结构地震响应的计算模型[J]. 力学学报,49(3):528–542. Lin G. 2017. A computational model for seismic response analysis of underground structures[J]. Chinese Journal of Theoretical and Applied Mechanics,49(3):528–542 (in Chinese).
缪惠全,李杰. 2018. 基于物理机制的工程场地地震动相干函数模型[J]. 中国科学:技术科学,48(2):209–216. Miao H Q,Li J. 2018. The coherence function model of ground motion in engineering sites based on physical mechanism[J]. Scientia Sinica:Technologica,48(2):209–216 (in Chinese).
张楚汉,金峰,王进廷,徐艳杰,潘坚文. 2016. 高混凝土坝抗震安全评价的关键问题与研究进展[J]. 水利学报,47(3):253–264. Zhang C H,Jin F,Wang J T,Xu Y J,Pan J W. 2016. Key issues and developments on seismic safety evaluation of high concrete dams[J]. Journal of Hydraulic Engineering,47(3):253–264 (in Chinese).
张翠然,陈厚群,李敏. 2011. 采用随机有限断层法生成最大可信地震[J]. 水利学报,42(6):721–728. Zhang C R,Chen H Q,Li M. 2011. Generation of the maximum credible earthquake by using the stochastic finite fault method[J]. Journal of Hydraulic Engineering,42(6):721–728 (in Chinese).
钟万勰,林家浩,吴志刚,孙东科,张亚辉. 2000. 大跨度桥梁分析方法的一些进展[J]. 大连理工大学学报,40(2):127–135. doi: 10.3321/j.issn:1000-8608.2000.02.001 Zhong W X,Lin J H,Wu Z G,Sun D K,Zhang Y H. 2000. Development of analytical method for long-span bridges[J]. Journal of Dalian University of Technology,40(2):127–135 (in Chinese).
Hao H,Oliveira C S,Penzien J. 1989. Multiple-station ground motion processing and simulation based on SMART-1 array data[J]. Nucl Eng Des,111(3):293–310. doi: 10.1016/0029-5493(89)90241-0
Loh C H,Yeh Y T. 1988. Spatial variation and stochastic modelling of seismic differential ground movement[J]. Earthq Eng Struct Dyn,16(4):583–596. doi: 10.1002/eqe.4290160409
Luco J E,Wong H L. 1986. Response of a rigid foundation to a spatially random ground motion[J]. Earthq Eng Struct Dyn,14(6):891–908. doi: 10.1002/eqe.4290140606
Shearer P M. 2009. Introduction to Seismology[M]. 2nd ed. Cambridge: Cambridge University Press: 215−240.
Zerva A. 2009. Spatial Variation of Seismic Ground Motions: Modeling and Engineering Applications[M]. Boca Raton: CRC Press: 1–64.
Zerva A,Harada T. 1997. Effect of surface layer stochasticity on seismic ground motion coherence and strain estimates[J]. Soil Dyn Earthq Eng,16(7/8):445–457.
-
期刊类型引用(12)
1. 杨敬松，廖府，李智涛，卢泽葳. 一种基于Delaunay三角剖分的天然地震波初至实时拾取算法. 地球物理学进展. 2024(05): 1788-1799 . 
百度学术
2. 邱磊，李彩华. STA/LTA方法拾取天然地震波初至及其改进分析. 地球物理学进展. 2023(04): 1497-1506 . 百度学术
3. 裴玮来，周仕勇. 基于最大顺序统计量的P波初动极性自动读取方法的小震震源机制反演及在云南小江地区的应用研究. 地球物理学报. 2022(03): 992-1005 . 百度学术
4. 颜利君，庞瑶，唐淋. 四川地震编目自动处理系统产出结果分析. 四川地震. 2022(01): 20-23 . 百度学术
5. 张珂，张帆，张晖，王鑫，郝美仙，翟浩. 内蒙古敖汉旗震群的微震匹配定位及其发震断层研究. 地震学报. 2022(04): 594-607 . 本站查看
6. 刘磊，王怀富，王志峰. 基于长短时窗比的地震检测方法应用. 电声技术. 2022(09): 93-96 . 百度学术
7. 金星，林彬华，蔡辉腾，韦永祥，李军，郭阳，于伟恒. 基于气枪信号偏振特征的震相自动判别研究. 地球物理学进展. 2021(01): 8-18 . 百度学术
8. 王旭林，李玉海，王鹏，王晓阳，吕明慧. 基于改进Delaunay三角剖分算法的微地震数据SRV估算与应用. 地球物理学进展. 2021(02): 689-695 . 百度学术
9. 钟媛，胡星星，滕云田，陈波，何朝博，沈晓宇. STA/LTA算法在MEMS烈度计中的适用性分析. 地震研究. 2021(02): 208-215 . 百度学术
10. 周本伟，范莉苹，张龙，李珀任，房立华. 利用卷积神经网络检测地震的方法与优化. 地震学报. 2020(06): 669-683 . 本站查看
11. Haichao Ma，Risheng Chu，Minhan Sheng，Ziye Yu. Template matching for simple waveforms with low signalto-noise ratio and its application to icequake detection. Earthquake Science. 2020(Z1): 256-263 . 必应学术
12. XU Zhen，WANG Tao，XU Shanhui，WANG Baoshan，FENG Xuping，SHI Jing，YANG Minghan. Active Source Seismic Identification and Automatic Picking of the P-wave First Arrival Using a Convolutional Neural Network. Earthquake Research in China. 2019(02): 288-304 . 必应学术
其他类型引用(10)
计量
- 文章访问数: 1547
- HTML全文浏览量: 951
- PDF下载量: 60
- 被引次数: 22
