基于最大顺序统计量的P波初动极性自动读取方法的小震震源机制反演及在云南小江地区的应用研究

裴玮来，周仕勇

北京大学地球与空间科学学院，北京 100871

0 引言

对地震波形记录中的地震事件进行准确检测并拾取地震震相的到时与初动极性，是地震资料分析的基本内容之一.目前在地震信号采集之中，最常用且被广泛认同的单台地震事件检测与震相到时读取方法主要有以下几类：

(1)基于长短窗信号比的震相到时检测方法.通过计算长短两个时窗的信号的能量比值来判定信号的初至到时(Stevenson,1976).短时能量(STA)显著高于长时能量(LTA)时，设定阈值判定地震信号到达的时间.长短窗算法是最经典的地震波信号拾取算法之一，常被拿来与其他方法进行对比.为了增加算法的灵敏度和精确度，许多学者先后采用不同的特征函数替代地震序列值作为输入(Allen,1978,1982;Baer and Kradolfer,1987)，来提高震相检测的精确性.

(2)基于熵的概念进行震相检测的方法.较有代表性的方法包括以赤池信息准则构建目标函数的震相拾取算法以及互信息方法两类.前者通过熵的概念权衡模型复杂度和数据拟合优良性，构建目标函数(AIC函数)，并寻找给定时窗中AIC函数的最小值的位置来决定时窗内地震信号的到时(Maeda,1985).后者通过计算随机变量之间的互信息值，或以其他信息论概念为基础构建的目标函数来判定地震震相到时.互信息的概念如下(Shannon,1948)：

(1)

其中p(x,y)表征随机变量X与Y的取值的联合概率密度分布；p(x)与p(y)表征随机变量X和Y的边缘概率密度分布；Pmi(x;y)表征点互信息，衡量两个随机事件之间的相关性.一般而言，基于熵的概念进行震相检测的算法在选取的时窗较为合适时能够给出非常精确的结果.

(3)基于高阶统计量如偏度、峰度等参数开发的震相拾取方法，通过计算地震波形的峰度和偏斜度的极值点，寻找地震波震相到时(Nakamula et al.,2007;Ross et al.,2016;Saragiotis et al.,2002).其中偏斜度描述随机变量相对于均值的不对称程度，而峰度描述随机变量分布曲线的尖峭程度.

(4)以机器学习为代表的利用人工智能搭建神经网络进行地震波震相拾取的方法.通过模拟人类的学习行为来训练神经网络，从而达到完成特定任务的目的(Wang and Teng,1997;Zhao and Takano,1999;Hara et al.,2019;Ross et al.,2018;Zhou et al.,2019;Zhu and Beroza,2019;Wong et al.,2021).基于人工智能的方法虽然能给出较高的准确率，但需要提供大量标记样本数据来训练模型，因此在不同地区和数据间的迁移和应用存在一定程度的困难且较为耗时.

(5)各类杂交方法与其他方法.Zhang等(2003)对单分量波形取一系列滑动时窗进行离散小波变换,并在每个时窗内应用AIC准则自动拾取判断是否有P波震相；Bai和Kennett(2001)结合了包括长短窗、平均-短期-长期比值、互相关、PCA在内的多种不同方法进行震相拾取；蒋策等(2018)结合AIC方法与长短窗方法的结果自动选择合适的值作为震相到时；赵大鹏等(2012)利用直达P波信号的峰度和Kurtosis-AIC两个指标，采用峰度触发时间窗、AIC挑震相的思想来达到区域地震事件实时检测和直达P波到时精细拾取的目的.

以上的几种方法主要是针对单台、单分量记录通过振幅或者能量的变化来检测P波震相.虽然它们基本实现了快速、自动化的P波震相到时检测，但是以上的方法大多难以估计拾取到时的不确定度，给定量估计初至到时的误差带来了困难.

作为小震震源机制解反演的重要资料，P波初动的极性是地震资料分析中的又一项基础性工作.P波初动极性的判定在一定程度上基于P波初至到时的准确检测.因此现有的P波初动极性的自动判定方法大多数使用了上文提到的一种或多种P波震相到时读取算子.比如Pugh等(2016)依赖上文中提到的长短窗函数法.Ross等(2018)依赖人工智能拾取P波震相到时的方法.Chen和Holland(2016)依赖基于峰度、偏度与赤池信息准则的联合到时拾取算法.当P波到时被准确挑取时，P波初动的判定本身没有过多的争议和可开发的空间，多数情况下我们都能得到准确的判定结果，但是当其所依赖的P波初至到时拾取算子不够准确或者出现错误时，微小的到时拾取误差就有可能导致不一样的结果.虽然大多数前人提出的方法都经过了检验并取得了不错的结果，但在将其应用到具体问题中时，还是有一定的犯错的几率.在目前主流的P波初动极性判定方法中，仅有Pugh等(2016)具备对初动极性误差进行定量分析的能力.然而由于该方法采用的初至拾取算子的误差定量估计不够严谨，导致其对P波初动极性的错误率分析不够严谨可靠.

现有方法在P波初至与P波初动极性拾取方面存在的问题，归根结底是因为初至到时的不确定度难以估计.因此，一种数学上更加严密、能够考虑结果的不确定度，且到时拾取更加精确的自动化到时拾取算法是有必要且意义重大的.本文采用的新的自动化到时拾取算法，POI(Probability of arrival time and polarity based on Order statistics and Information theory)算法，是基于最大顺序统计量和互信息理论构建的.该方法能够定量计算给定时窗中震相到时的概率密度分布，使P波到时与初动极性的不确定度得到约束.该方法并在自动化反演震源参数等方面可以得到很好的应用.

1 方法原理

1.1 基于最大顺序统计量的P波到时与极性读取方法

POI算法主要由两步构成.第一步是对给定时窗中的地震波数据，使用互信息理论，构造合适的特征函数来拾取地震波的初至时刻.地震波到达往往导致波形记录的振幅突然变大.因此通过构造一个体现振幅变化的随机变量，观察其与表征地震波到时的随机变量的相关性，就可以准确地拾取出地震波初至到时.假设地震波形记录Ai，定义一个振幅值ε作为阈值(图1).其中阈值ε是运用互信息算法对波形的振幅值Ai进行区分从而作为确认初动到时的参数.其值表征噪音和信号的振幅的绝对值大小的分界值.超过该阈值的数据点或数据段被认为其来自地震信号，而小于该阈值的数据点或数据段来自背景噪音，即有

图1 互信息到时判定利用互信息构建的目标函数Z(t,ε)，取目标函数值最大的时刻tinit作为到时.定义式中不同的事件用不同颜色区分出来.此处阈值约为300.Fig.1 Arrival time picking using pointwise mutual informationWhen amplitude threshold is given,arrival time is determined by maximizing objective function Z(t,ε).Here threshold is about 300.

(2)

由互信息理论，与表征振幅值变化的随机变量Y相关性最高的到时随机变量X，标识了震相到时t.在阈值ε确定时，到时t的取值完全由下列目标函数决定：

(3)

其中m和n取0或1,Pm n代表联合概率密度分布，其涵义为

P(X=x0,Y=y0)=P00(t,ε)=Pr{i≤t∧|Ai|≤ε},

P(X=x0,Y=y1)=P01(t,ε)=Pr{i≤t∧|Ai|>ε},

P(X=x1,Y=y0)=P10(t,ε)=Pr{i>t∧|Ai|≤ε},

P(X=x1,Y=y1)=P11(t,ε)=Pr{i>t∧|Ai|>ε}.

(4)

用tinit指代使互信息目标函数Z(t,ε)最大的t的取值，tinit即为阈值已知时认定的真实地震波到时：

(5)

若已知ε的概率密度分布，则可直接推出到时tinit的概率密度分布：

(6)

上式中I为逻辑函数，其括号内条件成立则取值为1，否则取值为0.

方法的第二步是通过引入最大顺序统计量的概念来计算合理的参数阈值ε的概率密度分布，从而使得到时tinit的概率密度分布可以直接通过等式(6)计算出来.

理论上阈值的取值范围可以是从0到正无穷.给出不同的阈值，或者阈值的概率密度分布，所推导出的到时，或到时的概率密度分布往往是不同的.我们希望给出的阈值的概率密度分布能够使得最终推导出的到时的概率密度分布接近真实的到时.为此需要推导能够计算出真实到时tinit的阈值ε所满足的条件.

我们提出最大顺序统计量准则：阈值ε应当与其对应到时tinit前置的背景噪音Asamp，在正态分布假设下的最大顺序统计量相近，即

max(|Asamp|)≈ε,(7)

上式左端表示在正态分布假设下，前置噪音样本Asamp的最大顺序统计量.其中前置噪音样本定义为

Asamp=ξ(ε)={Ai:|Ai|≤ε,i≤η(ε)}.

(8)

这里引入正态分布假设是因为真实的噪音往往是服从正态分布的.因此真正的阈值所对应的前置噪音应当也是大致服从正态分布.在正态分布假设下，前置的背景噪声可以认为是服从如下总体的随机采样：

Asamp～N(0,σ2),(9)

其中σ代表总体的噪声水平，可以用前置噪音样本估计：

(10)

此时前置噪音的最大顺序统计量的概率密度分布可以由下式给出：

(11)

上式中erf代表误差函数，card函数为返回元素个数的函数.考虑所有的阈值取值的情况，若阈值取值的概率密度分布已知，则前置噪音样本的最大顺序统计量的概率密度即为边缘分布：

=x|ε)×p(ε)dε,(12)

式(12)代表参数ε边缘化的过程，等式左边为将阈值边缘化后最大顺序统计量的概率密度分布.根据最大顺序统计量准则，合理的阈值的概率密度分布应该使其推出的最大顺序统计量的总体的概率密度分布接近自身，即

p(ε=x)≈p(max(|Asamp|)=x),(13)

因此有

式(14)中，待求的概率密度分布p(ε)为方程中的唯一未知量.因此可以通过解方程的形式将p(ε)解出，再由式(6)即可求得初至到时的概率密度分布.

初至波峰A定义为初至到时后第一个极值点的振幅值，其概率密度分布可以由初至到时的概率密度直接求得.当求解的阈值ε高于所有振幅绝对值|Ai|时，在所选时窗内无法定义初至到时与初至波峰A，初动的极性被判定为未知.除此之外，初至到时与初至波峰A的概率密度分布存在且可求.Pugh等(2016)中提出了完备的的贝叶斯极性判定公式，公式支持正负两种极性的计算.我们沿用其中的公式，计算初动的极性向上(Y=1)或向下(Y=-1)的概率：

(15)

将Y=1或Y=-1代入公式(15)中，所求得的概率即为初动极性为正或为负的概率.

1.2 基于P波初动的震源机制解与区域应力场反演流程

图2显示了本文工作设计的震源机制解与区域应力场反演流程.这套反演流程采用本文提出的POI P波初动极性判定方法，并综合了已有的HASH算法和SATSI算法.流程包括初动极性拾取、SH/P波振幅比计算、震源机制反演与应力场反演4个主要步骤.由于所有采用的方法本身都是能够实现自动化的，因此这套反演流程不需要人工的介入.

图2 本研究采用的自动化反演流程图图中的每一项计算过程都实现了自动化.其中POI极性拾取是本文采用的新的初动极性拾取方法.Fig.2 Automatic inversion process used in this studyEach step of inversion is independent of human involvement.POI method is the new polarity determination method adopted in the inversion process.

反演流程的第一个步骤是通过原始地震目录截取波形数据.对于给定的地震目录和连续波形记录，选取合适震中距的台站记录(120 km以内)，并通过区域一维速度结构估算P波震相与S波震相的理论到时，而后截取P波到时前10 s至后30 s的时间窗.在经过1～15 Hz的滤波后，选取P波初至到时前后各3 s的速度记录数据作为P波震相拾取时窗，为后续POI极性拾取做准备.分别选用P波到时前5 s到后2 s的时间窗，和S波理论到时前后共7 s的时间窗，分别计算P波与S波的位移记录的最大振幅值，为计算振幅比做准备.

反演流程的第二个步骤是通过波形数据分别同步计算P波初动极性和振幅比.P波初动极性的计算采用1.1节中提出的POI算法，设置质量控制参数(概率阈值，如：0.99)，保留所有概率阈值之上读取的极性数据.反演过程中使用的振幅比资料则需要在初步振幅比值的基础上进行了自由界面校正和台站校正.

由于放置在地表的地震仪器记录到的位移记录为入射波与反射波的叠加，而反演震源机制解的过程中使用的振幅比资料事实上只有入射波的部分，因此需要对初步振幅比资料进行自由界面校正.在众多种类的振幅比资料中，SH/P振幅比资料稳定性较好，在自由界面校正的过程中简单易算，不需要考虑P与SV的互相转化的问题.而SV/P以及S/P在计算过程中由于SV在自由界面表面反射的过程中存在大于临界角的情况，需要对入射角度进行严格分类计算.因此本研究的振幅比资料主要采用SH/P振幅比资料.尽管吴大铭等(1989)通过研究不同入射角范围的反透射系数，认为总体而言针对P波的自由界面校正是不必要的，考虑到本研究中入射角的展布范围，实际反演中还是对振幅比值进行了自由界面校正.自由界面校正的公式参考Aki和Richards(2002).由于定位结果中包含一定的误差，而这些误差可能会对振幅比的自由界面校正量带来影响，我们因此根据hypoinverse程序给出的深度误差计算了校正量的上下界.对于深度定位误差过大导致的振幅比的自由界面校正系数难以确定的情况，舍弃该振幅比值，仅保留在校正量误差在以2倍以内的数据，进行了自由界面校正.

台站校正是指通过计算理论辐射花样，与实际台站记录到的振幅比资料做对比，在一定程度上修正浅层地壳，特别是台站下方的场地效应以及衰减效应.台站校正整体上是一个经验性的校正.本研究中计算了所有台站的台站校正，采用的处理流程遵循Hardebeck和Shearer(2003)的处理流程.将经过自由界面校正与台站校正后得到的振幅比资料用于震源机制解反演.

反演流程的第三个步骤是通过HASH算法，使用第二个步骤计算的初动极性和振幅比数据对所有地震的震源机制进行反演和筛选.HASH程序的源代码是计算S波与P波的振幅比.算法考虑了地震深度、一维速度结构、P波初动极性等数据资料的不确定度，对同一地震事件采用不同的地震深度与不同的速度结构组合，反演地震震源机制，并在所有结果中进行聚类，然后返回可能的震源机制解.我们将HASH源码中计算S/P的目标函数改为计算SH/P的目标函数，其他流程不变.

反演流程的最后一步是通过SATSI算法，使用第三个步骤计算的震源机制解作为数据资料，反演区域应力场的应力主轴方向以及三个主应力的大小关系.SATSI算法是基于法向应力反演，默认断层滑动沿断层面上的剪切应力方向发生的应力张量求解技术.SATSI算法可以在保持连续性的前提下同时求解多个子区域的应力场，从而观察研究区域的应力主方向的空间分布与时间变化.

2 研究区域与研究资料

小江断裂带位于青藏高原东南缘的川滇菱形块体东侧，毗邻华南地块，走向近南北，全长约400 km，分为北、中、南三段，是川滇地区的主要活动断裂之一.多个研究结果显示，由于川滇块体向东南方向远离青藏高原，并围绕东喜马拉雅构造顺时针旋转，小江断裂带的主要运动形式以左旋走滑运动为主.小江断裂带主断层周边分布多条次级断裂，断层的几何形态与运动状态均比较复杂.为了研究小江断裂带中北段以及南段的地震活动性及断层结构，2015年至2019年间，北京大学周仕勇课题组联合中国地震局地球物理所许力生研究员沿小江断裂带走向布设了一个以宽频带仪器为主的流动地震台阵(图3).台阵覆盖小江断裂带的主断层结构，台站间距平均小于20 km.流动台阵的全部连续三分量记录以及昆明、东川两个固定台的连续波形记录构成了本研究的主要资料来源.

本研究中待反演震源机制的地震目录由自动化监测算法结合人工处理得到.利用自动化震相拾取聚类算法对近震的P、S初至震相进行初步的拾取和聚类，得到了包含1953个地震事件的初始地震目录.为了进一步提高震相目录的准确性和震相数量，增加求解资料，通过ObsPy(Beyreuther M et al，2010)和crazyseismic软件(Yu et al.，2017)对初步地震目录进行了人工筛选以及震相增补，去掉误识别并且补充初步拾取过程中遗漏的P波初至震相.经过人工筛选和扩充后的震相文件包含约15000个P波震相到时.它们来自1866个地震事件的所有震中距小于120 km的台站的可识别记录.平均每一个事件的P波震相记录数量大于8.

对于1866个通过筛选的地震事件，采用hypoinverse程序对所有事件进行重定位(Klein，2002).重定位所采用的速度模型参考吴建平等(2006).经过重定位后，选取P波记录数大于5，到时残差小于0.5 s，水平定位误差小于3 km，垂直定位误差小于5 km的地震事件，共674个符合条件的事件被挑选出来(图3).空间上看这些地震事件大多数分布在小江断裂带主断层的附近，少数事件有丛集特征的现象.

图3 本研究选取的研究区域，使用的地震台与选取的674个地震空间分布图倒三角代表固定台站(昆明台、东川台)，正三角代表宽频带流动地震台站，圆点代表地震位置.Fig.3 Demonstration of stations and earthquake catalog used in this studyInverted triangles represent permanent stations,i.e.,Dongchuan station and Kunming station.Triangles represent mobile stations.The catalog contains 674 events and they distribute evenly along the main fault.

3 结果分析

3.1 初动极性判定结果

为了检测方法的有效性，以上节确定的研究区674个定位精度较高的地震事件为基础资料，按照1.2节中的反演流程首先进行P波初动的极性判定.研究区台阵共获得了来自这674个地震的6688条地震记录，将这6688条地震记录截取的P波到时前3 s至后3 s时窗作为输入数据进入POI，开展P波初至到时与极性的概率计算.对于6688个输入数据，计算了每一条的P波初动极性概率以及时窗中的震相到时.图5至图8显示了其中4组计算案例.我们选取不确定度小于0.01的数据(初动极性向上的概率大于0.99或初动极性向下的概率大于0.99)作为有效数据(impulsive)，用作下一步震源机制解反演的输入资料，而将其他数据作为可靠性较低(emergent)或者极性未知(unknown)的初动数据舍弃.在6688组数据中，共有4084条记录被选中.

另一方面，为了验证新方法拾取初动极性的准确性，我们对6688组Z分量数据的初动极性进行了人工拾取，以便与自动拾取的初动极性做比较.人工拾取时采用传统的分类方法，将极性分为向上、向下、未知三类，并且对于向上和向下两类判定出极性的数据，按照质量又将每一类极性的波形分成了impulsive和emergent两类.6688组数据中，人工辨别极性清晰可读，质量较好的数据(impulsive)有2490条，极性可读但质量一般的数据(emergent)有849条，而极性未知的数据(unknown)有3349条.图4给出人工读取与POI自动读取结果的比较图.如果能假设人工能清晰读取的2490条极性结果是正确的，我们可以看到，POI自动读取的正确率在95%左右.此外POI还能将大量人工难以判断的极性高概率地识别出来.

图4 POI算法初动极性判定结果与人工初动极性判定结果的比较图中的每一行代表一类人工拾取初动结果，而每一行中三种极性的占比该类别数据被POI自动算法拾取为向上、未知与向下3种情况的概率.从图中可以看出Impulsive质量的数据判定结果一致性高，对人工判定极性为上、下的数据分别有高达96.5%和92.5%的拾取正确率.对于极性质量一般的数据判定的正确率平均超过75%.Fig.4 Comparison between automatic determined polarities and manual picking polaritiesEach line represents a category in manual picking.For each category,the probability of up,unknown and down polarities decided by POI algorithm is demonstrated.It′s obvious that our algorithm works well for impulsive polarities,with precision higher than 96.5% and 92.5% for up and down polarities,respectively.Our algorithm also has an accuracy above 75% for emergent polarities.

图5 初动极性判定案例人工标记该条记录极性向上，质量等级为Impulsive.POI算法计算其初动极性100%的概率向上.图中左上部分展示的是阈值的概率密度分布，右下图展示的是与之对应的初至到时的概率密度分布.主图中的红线和绿线展示两者之间是如何对应的.本例中到时几乎百分之百集中在3.048 s左右.Fig.5 An example of polarity determinationManual polarity picking classification of the seismic trace is up-impulsive.Polarity probability calculated by POI is 100%upward.The top left picture shows the pdf of threshold value.The down right picture shows the pdf of arrival time.Their correspondence is demonstrated in main graph.In this case pdf of arrival time is almost concentrated at 3.048 s.

图6 初动极性判定案例人工标记该条记录极性向上，质量等级为Emergent.POI算法计算其初动极性100%的概率向下.图中左上部分展示的是阈值的概率密度分布，右下图展示的是与之对应的初至到时的概率密度分布.主图中的红线和绿线展示两者之间是如何对应的.本例中到时几乎百分之百集中在3.271 s左右.Fig.6 An example of polarity determinationManual polarity picking classification of the seismic trace is down-emergent.Polarity probability calculated by POI is 100% downward.The top left picture shows the pdf of threshold value.The down right picture shows the pdf of arrival time.Their correspondence is demonstrated in main graph.In this case pdf of arrival time is almost concentrated at 3.271 s.

图7 初动极性判定案例人工标记该条记录极性未知.POI算法计算其初动极性100%的概率未知.图中左上部分展示的是阈值的概率密度分布，右下图展示的是与之对应的初至到时的概率密度分布.主图中的红线和绿线展示两者之间是如何对应的.本例中到时几乎百分之百大于6 s，意味着POI算法判定在本段波形数据中不存在初动.因此极性为未知.Fig.7 An example of polarity determinationManual polarity picking classification of the seismic trace is unknown.Polarity probability calculated by POI is 100% unknown.The top left picture shows the pdf of threshold value.The down right picture shows the pdf of arrival time.Their correspondence is demonstrated in main graph.In this case pdf of arrival time is almost 100% above 6 s,which means POI method is almost 100% sure that there is no seismic arrival in this waveform section.Thus the polarity is unknown.

图8 初动极性判定案例人工标记该条记录极性未知.POI算法计算其初动极性100%的概率向上.图中左上部分展示的是阈值的概率密度分布，右下图展示的是与之对应的初至到时的概率密度分布.主图中的红线和绿线展示两者之间是如何对应的.本例中到时几乎百分之百集中在2.956 s左右.Fig.8 An example of polarity determinationManual polarity picking classification of the seismic trace is unknown.Polarity probability calculated by POI is 100% upward.The top left picture shows the pdf of threshold value.The down right picture shows the pdf of arrival time.Their correspondence is demonstrated in main graph.In this case pdf of arrival time is almost concentrated at 2.956 s.

3.2 震源机制解反演结果

采用HASH程序(Hardebeck and Shearer,2002,2003)，利用POI读取的P波初动极性与SH/P波的振幅比资料，反演了上述674个地震事件的震源机制解.

考虑到HASH算法的多解性，需要对所得结果进行质量评估.本研究利用多种参数综合指标判断结果的质量，并对震源机制解进行分级与筛选.筛选标准包含断层面角度的不确定度数(Varest)、该组解占所有可能解的比例(Prob)、符合条件的解的总数量(Q-Number)、错误初动极性百分比(Bad Pol)、台站程度覆盖百分比(STA cover)、振幅比误差(SH/P error)等多个方面(表1).结合HASH程序中的基础的质量评定标准，只有在两次质量标准评级均在C或C以上，或者第一次评为A，第二次评级为D或D以上的解才会被保留.

表1 震源机制解质量分级表Table 1 Quality grading of focal mechanisms

图9展示了发震时刻于2016年9月7日的地震事件经过HASH反演算法所得到的震源机制解.图中标注了该事件在反演过程中使用到的初动波形数据以及POI算法计算出的极性向上的概率.从图中可以看出POI算法对初动极性的拾取是基本正确的，采用自动反演流程得到的震源机制解也与初动极性本身大致吻合.

图9 发生于2016年9月7日07时17分51.47秒的地震事件的P波初动极性拾取情况与震源机制解该事件的走向、滑动角及倾角分别为330°，75°，-17°.Fig.9 The polarities and focal mechanism of event occurred at 2016-09-07T07:17:51.47The strike,dip and rake angle are 330°，75°，-17°,respectively.

图10展示了筛选后得到227个质量较高解的地震事件的震源机制结果.其中走滑类型、正断层类型、逆冲断层类型和混合型地震事件的数量分别为87，84，31，25，它们占总事件数量的百分比分别为38.5%，37%，13.5%.11%.走滑类型是占比最多的地震事件类型，其次是正断层型，这与前人计算出的震源机制解类型大致一致(严川,2015)，也与区域的构造背景和GPS观测数据相一致.

图10 4种震源机制解的空间分布情况(a)正断层事件分布;(b)逆冲事件分布;(c)走滑事件分布;(d)混合类型事件分布.Fig.10 Spatial distribution of four type of focal mechanisms(a)Normal type events;(b)Thrust type events;(c)Strike-slip type events;(d)Hybrid type events.

分类型看，走滑型震源机制解的断层面展布仍然以南北和东西方向为主，展示了对小江主断裂带的南北走向特征的刻画.虽然偶尔有滑动方向不严格遵从左旋走滑的地震事件发生，大多数沿小江断裂带主断层分布的事件的震源机制解仍然是左旋走滑，也与整个区域的构造背景符合.

3.3 应力场反演结果

我们进一步统计分析了这227个地震事件的P轴和T轴的空间方位分布.结果见图11.P轴的方位分布整体集中在西北-南东方向，T轴的方位角分布整体更为集中，其优势走向仍为北东-南西.P轴与T轴的优势展布方向与前人的结果一致(吴建平等，2004)，与区域构造背景和GPS观测数据相符合(张培震等，2002).

图11 PT轴展布(a)P轴方位角展布;(b)T轴方位角展布.Fig.11 Strike distribution of P and T axes(a)Strike distribution of P axes;(b)Strike distribution of T axes.

更进一步，采用Hardebeck和Michael(2006)提出的基于法向应力的SATSI反演算法，使用Martínez-Garzón等(2014)开发的matlab程序包，利用227个小震的震源机制解对小江断裂带的区域应力场进行了反演.

由于研究区域较小，地震活动分布不均匀.考虑到程序反演出稳定结果对资料的数量有限制，仅将研究区域划分成南北两块，以25.7°为分界线，分别反演南北两个区域的平均应力张量.反演过程中，为了客观不加入任何先验条件，我们将每个地震事件的两个节面都作为可能的断层面输入程序，并将正确断层面数量的比例保持在默认值0.5，阻尼值采用程序默认值2000.反演结果由图12显示.

图12 应力主轴方向分布图红、绿、蓝依次代表最大至最小应力主轴.(a)南区应力主轴取向;(b)北区应力主轴取向.Fig.12 Principal axes of stress fieldOrientation of maximum,medium and minimum principal axes are represented by red,green and blue dots,respectively.(a)Orientation of axes of south region;(b)Orientation of axes of north region.

从图12中可以看出，在南北两个区域中，应力主轴的方向几乎没有改变.最大压应力轴与竖直方向有一定的偏离，次大压应力轴与最大压应力轴的走向相同，均为北西-南东方向，与上节中P轴的反演结果一致，而最小压应力轴则是与T轴的展布方向基本一致.最大压应力轴与次大压应力轴的倾角不确定度较大，但水平最大主应力方向的不确定非常小.与之对应的，最小主应力轴的方位角和倾角都被约束的非常好，且与前面得到的T轴方向相近，为北东-南西的近似水平方向.

R值表征了三个方向主应力的相对大小，其定义为：

横向对比两个区域的R值，北区的R值为0.1，南区的R值为0.13，二者非常接近.从R值中可以推出，最大压应力轴所对应的应力大小与中间压应力轴相近.这说明P轴方向与竖直方向的应力较为接近，而他们都远远大于最小压应力轴T轴的主应力值.

4 讨论与结论

传统的人工震相到时读取与P波初动极性判断不仅费时，且结果严重依赖于人工的经验性专业知识，效率低、主观性强，已不适合处理基于密集台阵的丰富的小震记录.本文采用的POI算法是从互信息理论出发的.互信息方法本身是成熟的信息学理论，在数学、物理等众多学科中有广泛的应用.已有科研工作者在地震信号的拾取中使用互信息理论，以往的互信息方法往往使用互信息值或者其他基于互信息的函数值作为目标函数，然后设立判据寻找初至到时(李辉等，2007).而本文的方法提供了一个完全不同的解决问题的思路.该方法首次在使用互信息方法，将到时拾取问题转化为模型参数抉择问题后，引入最大顺序统计量准则的概念，来解算参数的概率密度分布，并最终推导出初至到时的概率密度.POI算法没有局限于互信息本身，原理简单，使用的参数少，数学推理严密，同时在计算的过程中自然地引入了正态分布假设，减少了结果的随机性和主观设定先验信息可能造成的误识别.同时POI算法还可以给出不确定度的估计.对于信噪比较低的数据，其他方法往往只能通过设定阈值，对数据进行取舍，而本研究采用的POI算法能够判断该段数据包括P波初动的概率，从而给出概率化的解，这在目前的初动极性检测算法中是没有的.

目前流行的地震目录完备震级检测和基于海量样本训练的人工智能方法(Zhou et al.,2018;Zhou et al.,2019;Zhu and Beroza,2019)虽然有效地解决了地震事件的自动检测与震相到时的自动读取问题，但由于地震记录的区域特性，通常需要针对研究区寻找包含研究区区域特征的新地震记录样本集进行重新训练，给基于海量样本训练的人工智能方法的移植性造成困难.相对人工智能类的方法，本文提出的方法无需搜集海量训练样本，因此有更高的工作效率与可移植性.

本文的第二部分主要内容是依据POI算法，联合已有的成熟的反演程序HASH和MSATSI，开发了一套自动化反演震源机制解及区域应力场的工作流程.该工作流程所需要的输入数据包括经过定位的地震目录，区域速度结构与连续波形记录.本文采用北京大学与中国地震局地球物理研究所在小江区域布置的2015—2019年台阵数据，以及674个定位误差较小的可靠事件目录作为资料，运用POI算法，对P波初动到时与极性进行了自动拾取，计算了振幅比资料，并反演了小震的震源机制解以及区域平均应力场.反演过程中得到了4084条P波初动极性记录以及227个震源机制解.经过与不同结果的比对，P波初动极性的拾取结果、震源机制的反演结果以及应力场推断结果均与前人研究较为吻合，符合基本的区域构造特征和其他数据观测资料.因此POI自动拾取算法，以及其联合HASH与SATSI构建的自动化震源机制与应力场反演流程具有较高的可靠性和应用性和可推广性.我们希望更多地震学家能够试用本文提出的方法与反演流程，将其应用到不同区域和不同观测条件下，不断完善该方法，寻找最优反演流程，并拓展其应用范围.

致谢中国地震局地球物理研究所许力生课题组慷慨地给我们提供他们的地震台阵记录资料.3位匿名评审专家给本文提出了许多很好的建设性修改意见，在此衷心感谢他们.