地震震源监测方法、装置、计算机设备、存储介质与流程

1.本发明涉及地震监测技术领域，特别是涉及地震震源监测方法、装置、计算机设备、存储介质。


背景技术：

2.当震级较大的地震发生时，它的破裂路径和速度是决定其破坏能力的重要因素。基于大量自由参数紧密拟合的观测值，利用有限断层反演方法来确定破裂的运动学特征。然而，要建立一个运动学模型包含有断层结构的先验信息，其特点是固有的非唯一性，并且不能确保地震动力学方面的机械一致性。不断增加的计算资源，允许开发观测约束动态破裂模型，以补充数据驱动的分析。这种破裂场景提供了对复杂断层滑动的物理上自洽的描述，而它们的复杂性限制了可行数值实验的总数。大规模、密集的地震阵列仪器的兴起，使得在空间和时间上跟踪地震的互补技术成为可能。这种方法需要非常有限的先验知识，以简单和快速的方式对相干高频能量辐射成像。
3.旋转可以从平移运动的空间梯度导出，阵列观测法是通过平动分量计算旋转分量的一种测量方法。此方法利用大量平动传感器组成传感器阵列，在局部均匀结构上，利用有限差分方法作用于校准良好的传感器测量的平移分量，计算得到对应的地震旋转信号。该方法要求地震仪本身对旋转运动敏感，虽然所需的最少台站数量是三个，但通常需要更多的几个方位覆盖相当好的台站才能进行稳定的观测。另外，基于阵列观测法对旋转分量测量的质量控制非常困难。参考台站处的相关空间梯度，可以通过一定数量的密集台阵来估计(oliveiraandbolt,1989；spudichetal.,1995)。frohlichandpulliam(1999)指出，与基于旅行时的方法相比，经典的单站方法存在一些模糊性，例如存在180
°
反方位角误差问题。平移和旋转运动数据的联合分析可以帮助克服这些缺点。
4.使用阵列数据对地震特性进行成像的最常用技术可以分为两类，它们都基于分析p波的相位信息。第一类方法基于传统的阵列测量。称为反投影方法，地震能量辐射通过应用阵列波束形成技术成像。反投影首次在2004年苏门答腊-安达曼地震源(震中)得到成功证明。方向性效应被用来表征断层机制(kr
ü
gerandohrnberger,2005；ishiietal.,2005)。第二类地震破裂追踪方法是通过单站估计入射波的反方位角(backazimuth，baz)。在偏振分析中，标准地震仪的三个平移分量可用于估计入射波的后方位角和入射角。cochardetal.(2006)提出：在自由表面上，x、y、z方向的三个旋转分量可以被表示为：
[0005][0006]
其中表示位移波场的时间倒数。因此旋转可以由横向运动的空间梯度推导出来。
[0007]
拜耳等(2012)开发了一种通过局部和区域p波到达的偏振分析来跟踪移动源的单
站方法。


技术实现要素：

[0008]
本技术提供了一种地震震源监测方法、装置、计算机设备和存储介质。
[0009]
第一方面提供了一种地震震源监测方法，所述方法包括：
[0010]
获取监测点在预设历史时间段内采集的地震数据，所述地震数据包括：立体空间中相互垂直的第一坐标轴上的第一平动分量和第一旋转分量，第二坐标轴上的第二平动分量和第二旋转分量，第三坐标轴上的第三平动分量和第三旋转分量；
[0011]
确定震源至所述监测点的候选反方位角；
[0012]
根据所述第一平动分量、所述第二平动分量和各所述候选反方位角，得到各所述候选反方位角对应的横向加速度；
[0013]
确定各所述横向加速度和所述第三旋转分量之间的零滞后互相关系数；
[0014]
根据所述候选反方位角和所述零滞后互相关系数，确定所述震源至所述监测点的最终反方位角。
[0015]
在一些实施例中，所述确定震源至所述监测点的候选反方位角，包括：
[0016]
将0
°
到360
°
反方位角以预设度数间隔划分为多份，得到多个候选反方位角。
[0017]
在一些实施例中，所述根据所述候选反方位角和所述零滞后互相关系数，确定所述震源至所述监测点的最终反方位角，包括：
[0018]
选择与最大的零滞后互相关系数对应的候选反方位角为所述最终反方位角。
[0019]
在一些实施例中，所述确定震源至所述监测点的候选反方位角，包括：
[0020]
根据幅值最大的第一旋转分量和幅值最大的第二旋转分量，得到估算反方位角；
[0021]
判断所述估算反方位角是否小于0，若所述估算反方位角不小于0，则以所述估算反方位角为所述候选反方位角；若所述估算反方位角小于0，则使所述估算反方位角加上180
°
，得到所述候选反方位角。
[0022]
在一些实施例中，所述确定震源至所述监测点的候选反方位角，包括：
[0023]
判断所述零滞后互相关系数是否大于0，若大于0，则使所述候选反方位角加上180
°
，得到所述最终反方位角；若小于或等于0，则所述候选反方位角为所述最终反方位角。
[0024]
在一些实施例中，所述根据幅值最大的第一旋转分量和幅值最大的第二旋转分量，得到估算反方位角，包括：
[0025]
将所述幅值最大的第一旋转分量和所述幅值最大的第二旋转分量代入估算反方位角计算公式，得到所述估算反方位角，其中，所述估算反方位角计算公式为：
[0026][0027]
式中，θ
baz
为估算反方位角；re为幅值最大的第一旋转分量；rn为幅值最大的第二旋转分量。
[0028]
在一些实施例中，所述地震数据还包括：地震波中的p波的波速和s波的波速以及所述p波和所述s波的传播时差；
[0029]
所述方法还包括：根据所述p波的波速、所述s波的波速以及所述传播时差，确定所述地震波的震源至监测点的距离。
[0030]
第二方面提供了一种地震震源监测装置，所述装置包括：
[0031]
数据输入单元，用于获取监测点在预设历史时间段内采集的地震数据，所述地震数据包括：立体空间中相互垂直的第一坐标轴上的第一平动分量和第一旋转分量，第二坐标轴上的第二平动分量和第二旋转分量，第三坐标轴上的第三平动分量和第三旋转分量；
[0032]
候选确定单元，用于确定震源至所述监测点的候选反方位角；
[0033]
横向加速度计算单元，用于根据所述第一平动分量、所述第二平动分量和各所述候选反方位角，得到各所述候选反方位角对应的横向加速度；
[0034]
互相关系数计算单元，用于确定各所述横向加速度和所述第三旋转分量之间的零滞后互相关系数；
[0035]
结果输出单元，用于根据所述候选反方位角和所述零滞后互相关系数，确定所述震源至所述监测点的最终反方位角。
[0036]
第三方面提供了一种计算机设备，包括存储器和处理器，所述存储器中存储有计算机可读指令，所述计算机可读指令被所述处理器执行时，使得所述处理器执行上述所述地震震源监测方法的步骤。
[0037]
第四方面提供了一种存储有计算机可读指令的存储介质，所述计算机可读指令被一个或多个处理器执行时，使得一个或多个处理器执行上述所述地震震源监测方法的步骤。
[0038]
上述地震震源监测方法、装置、计算机设备和存储介质，获取监测点采集的立体空间的三个方向的平动分量和旋转运动分量；根据立体空间的三个方向的平动分量和旋转分量，确定震源至监测点的反方位角；根据振波中的p波和s波的波速以及传播时差，确定震源至监测点的距离。因此，单台法没有场地的限制以及人工布置的成本也得到大大减少，反方位角信息可以通过单台旋转分量和平动分量的结合来获取。互相关方法侧重于扭转运动(垂向轴的旋转)，偏振分析方法注重摇摆运动(水平轴的旋转)。基于水平旋转分量偏振分析的反方位角估计，不受微小但有影响的异步化的影响，因为只涉及一个数据记录器。偏振分析方法在使用起来更加灵活，且可以用于p、sv、sh、rayleigh波以及love波中。
附图说明
[0039]
图1为本技术一个实施例中地震震源监测方法的流程图；
[0040]
图2为本技术一个实施例中地震震源监测方法的利用互相关方法实现震源追踪流程图；
[0041]
图3为本技术一个实施例中地震震源监测方法的在反方位角为280
°
情况下，横向加速度与第三旋转分量波形图；横坐标对应时间(单位，秒)，纵坐标对应归一化振幅；实线表示第三旋转分量，虚线表示转化的横向加速度；
[0042]
图4为本技术一个实施例中地震震源监测方法的在270s到300s这个时间窗口内，互相关方法对每个反方位角求得的互相关系数，横坐标对应每个反方位角，纵坐标为对应的互相关系数；
[0043]
图5为本技术一个实施例中地震震源监测方法的利用偏振分析方法实现震源追踪流程图；
[0044]
图6为本技术一个实施例中地震震源监测方法的在整个时间窗口里，偏振分析方
法预测得到的反方位角和理论反方位角；其中黑色星号为每个时间窗口预测的反方位角，水平实线为理论的反方位角数值；横坐标为采样时间(单位秒)，纵坐标为反方位角的角度；
[0045]
图7为本技术一个实施例中地震震源监测装置的结构框图；
[0046]
图8为一个实施例中计算机设备的内部结构框图。
具体实施方式
[0047]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0048]
可以理解，本技术所使用的术语“第一”、“第二”、第三”等可在本文中用于描述各种元件，但这些元件不受这些术语限制。这些术语仅用于将第一个元件与另一个元件区分。
[0049]
一个实施例中提供的地震震源监测方法的实施环境图，在该实施环境中，可以包括计算机设备以及地震监测仪。
[0050]
计算机设备具有接口，例如可以为接口是api(applicationprogramminginterface，即应用程序接口)。地震监测仪为设置在检测点的采集的地震数据的装置，计算机设备根据地震监测仪采集的数据进行地震震源监测(如确定震源的反方位角)。
[0051]
需要说明的是，终端以及计算机设备可为智能手机、平板电脑、笔记本电脑、台式计算机等，但并不局限于此。计算机设备以及终端可以通过蓝牙、usb(universalserialbus，通用串行总线)或者其他通讯连接方式进行连接，本发明在此不做限制。
[0052]
如图1至6所示，在一个实施例中，提出了一种地震震源监测方法，该地震震源监测方法可以应用于上述的计算机设备中，该方法可以基于python程序实现。如图1所示，具体可以包括以下步骤：
[0053]
步骤101、获取监测点在预设历史时间段内采集的地震数据，地震数据包括：立体空间中相互垂直的第一坐标轴上的第一平动分量和第一旋转分量，第二坐标轴上的第二平动分量和第二旋转分量，第三坐标轴上的第三平动分量和第三旋转分量；
[0054]
可以理解的是，该步骤中，获取监测点采集的立体空间的三个坐标轴的平动分量和旋转分量，包括：
[0055]
采用六分量地震仪采集立体空间的三个方向的平动分量和旋转分量；该六分量地震仪，可以包括：三个光纤陀螺、三个加速度计和六分量信号处理模块，其中，三个光纤陀螺的敏感轴相互正交，三个加速度计的敏感轴相互正交；每个光纤陀螺的敏感轴分别与一个加速度计的三个敏感轴一一平行或重合；六分量信号处理模块的输入端连接每个光纤陀螺的检测信号输出端以及以及连接每个加速度计的检测信号输出端，六分量信号处理模块用于生成光纤陀螺所需的调制信号，以及根据光纤陀螺的检测信号输出端输出的检测信号，得到检测角速度，并对检测角速度进行误差补偿六分量信号处理模块还用于根据三轴加速度计的检测信号输出端输出的检测信号，得到检测平移加速度，并对检测平移加速度进行误差补偿。其中，建立空间直角坐标系，该空间直角坐标系包括相互正交的第一坐标轴(x轴)、第二坐标轴(y轴)和第三坐标轴(z轴)，光纤陀螺和加速度计的敏感轴分别在第一坐标
轴(x轴)、第二坐标轴(y轴)和第三坐标轴(z轴)上。
[0056]
进一步地，六分量地震仪包含三个正交的光纤陀螺和三个同样正交的加速度计，光纤陀螺用于测量物体的角速度信息，三轴加速度计用于测量搭载体的线运动信息，六分量信号处理模块通过处理光纤陀螺和加速度计的输出，能够提供完整的地震波场信息。
[0057]
步骤102、确定震源至监测点的候选反方位角；
[0058]
可以理解的是，候选反方位角可以是根据经验确定的，也可以是根据地震数据确定的。
[0059]
在一些实施例中，采用互相关方法实现震源最终反方位角的确定，则在互相关方法中上述步骤102中，根据经验设置候选反方位角可以包括：
[0060]
步骤1021a、将0
°
到360
°
反方位角以预设度数间隔划分为多份，得到多个候选反方位角。
[0061]
在一些实施例中采用偏振分析方法实现震源最终反方位角的确定，则在偏振分析方法中上述步骤102中根据地震数据确定候选反方位角可以包括：
[0062]
步骤1022a、根据幅值最大的第一旋转分量和幅值最大的第二旋转分量，得到估算反方位角；
[0063]
在该步骤中，将幅值最大的第一旋转分量和幅值最大的第二旋转分量代入估算反方位角计算公式，得到估算反方位角，其中，估算反方位角计算公式为：
[0064][0065]
式中，θ
baz
为估算反方位角；re为幅值最大的第一旋转分量；rn为幅值最大的第二旋转分量。
[0066]
步骤1022b、判断估算反方位角是否小于0，若估算反方位角不小于0，则以估算反方位角为候选反方位角；若估算反方位角小于0，则使估算反方位角加上180
°
，得到候选反方位角。
[0067]
步骤103、根据第一平动分量、第二平动分量和各候选反方位角，得到各候选反方位角对应的横向加速度；
[0068]
该步骤中，将第一平动分量、第二平动分量和各候选反方位角代入横向加速度公式中，得到各候选反方位角的横向加速度。
[0069]
其中，横向加速度计算公式为：
[0070][0071]
式中，ae为第一平动分量；an为第二平动分量；θ
baz
为候选反方位角；a
t
为横向加速度。
[0072]
步骤104、确定各横向加速度和第三旋转分量之间的零滞后互相关系数；
[0073]
该步骤中，将各横向加速度和第三旋转分量代入相关计算公式中，得到横向加速度和第三旋转分量之间的零滞后互相关系数。其中，相关计算公式为
[0074][0075]
式中，a
t
各个采集时刻对应的横向加速度；rz为各个采集时刻采集的垂向旋转分量；为横向加速度的平均值；为垂向旋转分量的平均值；xcorr(a
t
,rz)为横向加速度和垂向旋转分量之间的零滞后互相关系数。
[0076]
可以理解的是，六分量地震仪以预设的时间间隔(如1秒采集一次)采集地震数据，各个采集时刻对应的横向加速度(a
t
)是将各个采集时刻采集的为第一平动分量、第二平动分量、反方位角(候选反方位角或最终反方位角)代入横向加速度计算公式得到的。横向加速度的平均值是一段时间(一个时间窗口，如3秒)的横向加速度的平均值。
[0077]
步骤105、根据候选反方位角和零滞后互相关系数，确定震源至监测点的最终反方位角。
[0078]
该步骤中，横向加速度和第三旋转分量之间的零滞后互相关系数与震源的反方角之间有一定的联系，因此，根据候选反方位角和零滞后互相关系数，可以确定震源至监测点的最终反方位角。
[0079]
在一些实施例中，采用互相关方法实现震源最终反方位角的确定，则在互相关方法中上述步骤105中，根据候选反方位角和零滞后互相关系数，确定震源至监测点的最终反方位角可以包括：
[0080]
步骤1051a、选择与最大的零滞后互相关系数对应的候选反方位角为最终反方位角。
[0081]
在一些实施例中采用偏振分析方法实现震源最终反方位角的确定，则在偏振分析方法中上述步骤105中，根据候选反方位角和零滞后互相关系数，确定震源至监测点的最终反方位角可以包括：
[0082]
步骤1052a、判断零滞后互相关系数是否大于0，若大于0，则使候选反方位角加上180
°
，得到最终反方位角；若小于或等于0，则候选反方位角为最终反方位角。
[0083]
进一步地，可以理解的是，获取监测点在预设历史时间段内采集的地震数据，可以是获取当前时刻以往12s内的地震数据，监测点每1秒采集一次的频率采集地震数据，根据该地震数据确定震源至监测点的最终反方位角，可以将上述地震数据划分为3个时间窗口(每个时间窗口时长为4s)，以一个时间窗口为最小单位，采用上述实施例中根据每个时间窗口内的地震数据确定该时间窗口内的候选反方位角，进而确定各时间窗口内的最终反方位角，最后根据多个时间窗口的最终反方位角确定震源至监测点的反方位角。也就是，预测的最终反方位角大都分布在理论反方位角附近，个别时间窗口的结果误差较大，可能受到其它外界因素的影响，整体上预测结果与理论数值比较吻合，因此，将多个时间窗口的最终反方位角线型拟合，得到震源至监测点的反方位角。当然，划分时间窗口后，上述实施例中“幅值最大”、“最大的零滞后互相关系数”等是在单位时间窗口内最大。
[0084]
如图2至4，在一应用场景中，利用2021年5月21日武汉仪器采集到的云南大理发生6.4级地震数据，来测试互相关方法的可行性。其理论反方位角为280
°
，震源(震中)距测试台站距离约为1583km。
[0085]
如图2所示，上述方法可以包括：
[0086]
步骤a1、从旋转地震仪上获得三个平动分量以及旋转分量地震数据，其中截取记录数据长度为390s，从地震发生时开始。再对数据做移除仪器响应等预处理；
[0087]
步骤a2、将0
°
到360
°
反方位角以5
°
间隔划分为72份，对于给定的每一个反方位角，将第一水平分量和第二水平分量根据横向加速度公式，得到横向加速度。图3为在反方位角为280
°
情况下，横向加速度与第三旋转分量波形图，从图上可以看出波形能较好的拟合；
[0088]
步骤a3、将记录到的390s地震数据按30s为一个时间窗口均匀划分为13份。在每个时间窗口内，对所有候选反方位角和每个时间窗口进行循环。将得到的横向加速度a
t
与第三旋转分量rz代入相关计算公式中，横向加速度和第三旋转分量之间的零滞后互相关系数。图4是270s到300s时间窗口内，每个反方位角对应的互相关系数；
[0089]
步骤a4、选取互相关系数最大值位置对应的反方位角作为该时间窗口的最后结果。
[0090]
当所有的时间窗口循环结束后，可以得到整个时间道的反方位角预测结果。
[0091]
如图5、6，在一应用场景中，利用2021年1月12日在广州进行的甲烷爆炸实验来测试偏振分析方法的可行性。其理论反方位角为190
°
，甲烷源距测试台站距离约为305m。如图5所示，上述方法可以包括：
[0092]
步骤b1、从地震仪下载获得三个平动分量以及旋转分量地震数据，其中截取记录数据长度为15s，从甲烷源爆炸时开始。再对下载数据做移除仪器响应等预处理；
[0093]
步骤b2、将记录到的15s地震数据按1s为一个时间窗口均匀划分为15份。在每个时间窗口内，选择最大幅值的第一旋转分量和第二旋转分量来计算候选反方位角。
[0094]
步骤b3、在每一个1s时间窗口内，根据计算得到的候选反方位角将水平分量转化成横向加速度。
[0095]
步骤b4、将第一平动分量、第二平动分量和候选反方位角代入横向加速度公式中，得到各候选反方位角的横向加速度。
[0096]
步骤b5、如果横向加速度与第三旋转分量为正相关，则得到的反方位角则加上180
°
，反之不变。
[0097]
步骤b6、对每个时间窗口进行循环，得到整个时间道的反方位角预测结果，如图6所示。可以看到预测的反方位角比较均匀分布在理论方位角附近，但是也有部分时间窗口由于受到其它外界因素的影响，结果误差较大。
[0098]
在一些实施例中，上述地震波中的p波的波速和s波的波速以及p波和s波的传播时差；
[0099]
上述方法还包括：根据p波的波速、s波的波速以及传播时差，确定地震波的震源至监测点的距离。具体包括根据如下公式计算：
[0100][0101]
式中，其中s为震源至监测点的距离；v
p
为p波的波速；和vs分别为s波的波速；t1和t2分别为根据波形得到的p波和s波到时。
[0102]
在一个实施例中，提供了一种确定地震波相速度的装置，该确定地震波相速度的装置可以集成于上述的计算机设备中，如图7所示，该装置具体可以包括：
[0103]
数据输入单元711，用于获取监测点在预设历史时间段内采集的地震数据，地震数据包括：立体空间中相互垂直的第一坐标轴上的第一平动分量和第一旋转分量，第二坐标轴上的第二平动分量和第二旋转分量，第三坐标轴上的第三平动分量和第三旋转分量；
[0104]
候选确定单元712，用于确定震源至监测点的候选反方位角；
[0105]
横向加速度计算单元，用于根据第一平动分量、第二平动分量和各候选反方位角，得到各候选反方位角对应的横向加速度；
[0106]
互相关系数计算单元713，用于确定各横向加速度和第三旋转分量之间的零滞后互相关系数；
[0107]
结果输出单元714，用于根据候选反方位角和零滞后互相关系数，确定震源至监测点的最终反方位角。
[0108]
如图8所示，在一个实施例中，提出了一种计算机设备，该计算机设备可以包括通过系统总线连接的处理器、存储介质、存储器和网络api接口。其中，该计算机设备的存储介质存储有操作系统、数据库和计算机可读指令，数据库中可存储有控件信息序列，该计算机可读指令被处理器执行时，可使得处理器实现一种地震震源监测方法。该计算机设备的处理器用于提供计算和控制能力，支撑整个计算机设备的运行。该计算机设备的存储器中可存储有计算机可读指令，该计算机可读指令被处理器执行时，可使得处理器执行一种地震震源监测方法。该计算机设备的网络api接口用于与终端连接通信。本领域技术人员可以理解，图8中示出的结构，仅仅是与本技术方案相关的部分结构的框图，并不构成对本技术方案所应用于其的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
[0109]
处理器执行计算机程序时实现以下步骤：获取监测点在预设历史时间段内采集的地震数据，地震数据包括：立体空间中相互垂直的第一坐标轴上的第一平动分量和第一旋转分量，第二坐标轴上的第二平动分量和第二旋转分量，第三坐标轴上的第三平动分量和第三旋转分量；确定震源至监测点的候选反方位角；根据第一平动分量、第二平动分量和各候选反方位角，得到各候选反方位角对应的横向加速度；确定各横向加速度和第三旋转分量之间的零滞后互相关系数；根据候选反方位角和零滞后互相关系数，确定震源至监测点的最终反方位角。
[0110]
在一个实施例中，提出了一种存储有计算机可读指令的存储介质，该计算机可读指令被一个或多个处理器执行时，使得一个或多个处理器执行以下步骤：获取监测点在预设历史时间段内采集的地震数据，地震数据包括：立体空间中相互垂直的第一坐标轴上的第一平动分量和第一旋转分量，第二坐标轴上的第二平动分量和第二旋转分量，第三坐标轴上的第三平动分量和第三旋转分量；确定震源至监测点的候选反方位角；根据第一平动分量、第二平动分量和各候选反方位角，得到各候选反方位角对应的横向加速度；确定各横向加速度和第三旋转分量之间的零滞后互相关系数；根据候选反方位角和零滞后互相关系数，确定震源至监测点的最终反方位角。
[0111]
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，该计算机程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，前述的存储介质可为磁碟、光盘、只读存储记忆体(read-onlymemory，rom)等非易失性存储介质，或随机存储记
忆体(randomaccessmemory，ram)等。
[0112]
以上实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
[0113]
以上实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。