具体实施方式

以下结合附图和实施例具体说明本发明的技术方案。虽然附图中显示了本发明的示例性实施例，然而应当理解，可以以各种形式实现本发明，而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更透彻地理解本发明，并且能够将本发明的范围完整地传达给本领域的技术人员。

本技术领域技术人员可以理解，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)，具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语，应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非被特定定义，否则不会用理想化或过于正式的含义来解释。

随着高频GNSS接收机技术的发展和GNSS定位精度的提高，GNSS历元间差分方法可以在大地震近场区域实时、连续地提供mm/s精度的速度时间序列。为弥补常规地震预警系统中强震仪、地震仪在大地震近场区域因发生仪器倾斜、旋转以及量程超限等导致的数据不可靠，解决现有震源反演策略低普适性、低时效性的问题，本发明公开了一种基于单一矩形模型的地震震源快速确定方法及系统。本发明通过联合使用强震数据和高频GNSS数据，以同时做到P波到时的精确拾取和大震近场区域地震波的准确捕获。并且，本发明充分发挥了GNSS历元间差分方法的时效性及单一矩形模型的普适性，凭借震源参数遍历方法，以实现震源破裂长度、方向的快速确定，以期为地震预警、地震动场预测、震后灾害评估及快速救援等提供快速、可靠的震源参考信息。

图1为本发明实施例提供的一种基于单一矩形模型的地震震源快速确定方法的流程图。参照图1，实施例具体步骤如下：

S11、强震仪或地震仪采集地震台站原始运动信息，GNSS接收机采集GNSS观测站基本导航定位信息，并将采集到的所述地震台站原始运动信息、所述GNSS观测站基本导航定位信息发送到地面计算中心；

S12、所述地面计算中心根据所述地震台站原始运动信息，构建地震台站特征函数时间序列，并根据所述地震台站特征函数时间序列，判断地震台站是否发生地震触发。如果没有发生地震触发，所述地面计算中心继续判断地震台站是否发生地震触发；否则，所述地面计算中心开启地震模式，并根据所述地震台站特征函数时间序列，精确拾取地震台站初动时刻；

S13、所述地面计算中心根据本地存储的地震台站位置信息、所述地震台站初动时刻，确定震中位置；所述地面计算中心根据所述GNSS观测站基本导航定位信息，生成GNSS观测站运动信息；

S14、所述地面计算中心根据本地存储的GNSS观测站位置信息、所述GNSS观测站运动信息以及所述震中位置，估计点源震级，并根据所述点源震级估计单一矩形模型的初始断层尺寸；

S15、所述地面计算中心根据所述初始断层尺寸、所述GNSS观测站位置信息、所述GNSS观测站地面峰值参数观测值以及所述震中位置，反演单一矩形模型参数，并将所述震中位置、所述点源震级、所述单一矩形模型参数等地震预警信息发送到地震预警中心。

进一步地，所述S12中，所述地面计算中心根据所述地震台站原始运动信息，构建地震台站特征函数时间序列，并根据所述地震台站特征函数时间序列，判断地震台站是否发生地震触发，以及精确拾取地震台站初动时刻，通过以下具体流程实现：

第一步：构建同时兼顾频率和振幅特性变化的地震台站特征函数时间序列，具体形式如下：

式中，CF_k表示构建的第k时刻的地震台站特征函数值；x_k、x_k-1分别表示第k时刻和第k-1时刻的地震台站原始运动信息，如果地震台站配置强震仪，则所述地震台站原始运动信息为地震台站加速度；M_k为加权因子，与地震台站环境噪声和仪器采样频率有关，具体实施时可根据经验值给定，通常取3。

第二步：利用长短时平均法(STA/LTA)，判断地震台站是否发生地震触发，并对P波到时进行粗略拾取。标准公式(2)和递归公式(3)如下所示：

式中，STA和LTA分别表示特征函数时间序列在短、长时窗内的平均值，STA_k为第k时刻的特征函数时间序列在短时窗内的平均值，LTA_k为第k时刻的特征函数时间序列在长时窗内的平均值，第k-1时刻的情况类推；n为时窗内数据的个数，n_STA为短时窗内数据的个数，n_LTA为长时窗内数据的个数。标准公式和递归公式，在数值上是等价的。这里，实施例优选采用计算效率更高的递归公式。

当STA与LTA的比值大于预设的经验阈值时，k被判定为P波粗略到时。值得注意的是，长时窗和短时窗的窗口长度影响着经验阈值的选取，一般遵循的取值建议是：短时窗窗口长度相对于长时窗越小，拾取P波的敏感性越高，经验阈值应设置得大一些；短时窗窗口长度相对于长时窗越大，拾取P波的敏感性越低，经验阈值应设置得小一些。然后，以k为中心，向前截取长时窗、向后截取短时窗数据作为精确拾取地震台站初动时刻的时间窗口。通过这种较为合理的时间窗口设置，拾取结果一般具有较好的时效性和较高的精度。

第三步：采用赤池信息量准则(AIC)，精确拾取地震台站初动时刻。数学公式如下：

AIC_k＝k×log10{var(CF[1,k])}+(L-k-1)×log10{var(CF[k+1,L])} (4)

式中，L为精确拾取初动时刻时间窗口的数据长度；k为该时间窗口内从左到右的数据序号；var(CF[1,k])表示该时间段内数据方差。当AIC取得最小值时，第k个数据对应的时刻即为地震台站初动时刻。

进一步地，所述S13中，所述地面计算中心根据本地存储的地震台站位置信息、所述地震台站初动时刻，确定震中位置，具体包括：

至少有3个地震台站发生地震触发，才会开始估计地震震中。并且假设地震波传播过程中，地震波传播速度恒定不变、各向同性，具体步骤如下：

第一步：将触发的地震台站的绝对位置信息(大地坐标)，转化为以首个触发台站为原点的站心坐标。具体地，可以先将大地坐标(BLH)转化为空间直角坐标(XYZ)；再将空间直角坐标转化为站心坐标(ENU)。

第二步：在站心坐标系下，采用最小二乘平差，估计震中位置。根据P波到达各个地震台站的时间差，可以得到如下误差方程：

式中，(E₀,N₀)表示震中坐标，(E_i,N_i)表示地震台站坐标，其中下标i(i＝1,2,3,L,n)代表发生地震触发的地震台站序号，n代表此刻发生地震触发的地震台站总数；D_i表示第i个地震台站的震中距；V和T_i分别表示地震波传播速度和第i个地震台站的台站初动时刻；v_i,j则表示仪器观测误差、初动时刻拾取误差、地震台站坐标误差以及近似误差等的综合影响，其中下标i和j均代表发生地震触发的地震台站序号，且满足i＞j。注意的是，震中的确定还可使用网格搜索等方法，这里只是给出一种可行的实现方案，方便实施本发明时更加深入地、透彻地理解本发明的各个流程细节与特色，而并不是对本发明的一种限制。

第三步：将站心坐标系下的震中位置转化为大地坐标，具体地，可先将站心坐标转化为空间直角坐标，再将空间直角坐标转化为大地坐标。

进一步地，所述S13中，所述地面计算中心根据所述GNSS观测站基本导航定位信息，生成GNSS观测站运动信息，具体包括：

所述地面计算中心，根据载波相位观测值信息、伪距观测值信息、广播星历(所述GNSS观测站基本导航定位信息)，采用高频GNSS历元间差分方法，连续地、实时地获取GNSS观测站速度时间序列(所述GNSS观测站运动信息)。高频GNSS历元间差分方法的数学模型如下所示：

式中，t_k和t_k+1表示两个相邻的GNSS观测历元，其中下标k和k+1为GNSS观测历元标识序号；c为光速；dt_r和dt^s分别表示接收机和卫星钟差；和T_r ^s分别表示电离层和对流层延迟；为多路径误差；为其他误差和观测噪声的综合影响；和分别表示在历元t_k和t_k+1从卫星s到接收机r的单位方向向量；p^s(t_k)和p^s(t_k+1)分别表示在历元t_k和t_k+1时的卫星坐标；p_r(t_k)则表示在历元t_k时的接收机坐标；为载波相位观测值的历元间单差；cΔdt^s(t_k,t_k+1)＝c×[dt^s(t_k+1)-dt^s(t_k)]为卫星钟差变率；历元间接收机位移增量Δp_r(t_k,t_k+1)和接收机钟差变率cΔt_r(t_k,t_k+1)使用最小二乘估计；卫星位置和钟差通过广播星历计算得到，接收机位置使用伪距单点定位方法(Pseudo-rangeSingle Point Positioning)计算，电离层延迟变率对流层延迟变率和多路径的影响非常小，可忽略。Δp_r(t_k,t_k+1)为接收机在相邻历元间(t_k,t_k+1)的位移增量，除以采样间隔即为接收机平均速度。

进一步地，所述S14中，所述地面计算中心根据本地存储的GNSS观测站位置信息、所述GNSS观测站运动信息以及所述震中位置，估计点源震级，通过以下具体流程实现：

第一步：所述地面计算中心根据所述GNSS观测站位置信息、所述震中位置，计算GNSS观测站震中距。震中距计算公式如下：

EPID_i＝arccos(sinB_i sinB₀+cosB_i cosB₀ cos(L_i-L₀))×R_earth (7)

式中，B₀和L₀为震中的纬度和经度；B_i和L_i为GNSS观测站的纬度和经度；R_earth为地球平均半径，一般可取6371km；EPID_i为震中距。

第二步：所述地面计算中心根据GNSS观测站速度时间序列(所述GNSS观测站运动信息)，提取GNSS观测站地面峰值速度(GNSS观测站地面峰值参数观测值)。

式中，PGV表示GNSS观测站三维地面峰值速度；N_v(t)、E_v(t)、U_v(t)分别表示GNSS观测站t时刻在北方向、东方向以及垂直方向上的速度分量。

第三步：所述地面计算中心根据所述GNSS观测站震中距、GNSS观测站地面峰值速度(所述GNSS观测站地面峰值参数观测值)，采用经验公式(9)，估计地震点源震级。

log10(PGV)＝-5.025+0.741×Mw-0.111×Mw×log10(EPID) (9)

其中，Mw为矩震级；log10(*)表示以10为底的对数操作。

进一步地，所述S14中，所述地面计算中心根据所述点源震级，采用经验公式(10)，估计单一矩形模型的初始断层尺寸，具体地，所述初始断层尺寸包括初始断层长度和初始断层宽度。

式中，L0和W0分别为单一矩形模型的初始断层长度和宽度；Mw为矩震级，实际计算时，可用所述点源震级表示。

进一步地，所述S15中，地面计算中心根据所述初始断层尺寸、所述GNSS观测站位置信息、所述GNSS观测站地面峰值参数观测值以及所述震中位置，反演单一矩形模型参数，并将所述震中位置、所述点源震级、所述单一矩形模型参数等地震预警信息发送到地震预警中心，具体包括：

如图4所示，单一矩形模型参数包括震中位置(B₀，L₀)、震源深度(H)、断层长度(L)、主破解方向上震中到矩形断层一端的长度(l)和宽度(W)、震源破裂方向(θ)以及震源相对位置(k)，模型假设整个矩形面与震源保持同一高度。具体参数反演流程如下：

第一步：依据所述单一矩形模型参数的物理含义及先验信息，确定合理的震源参数遍历范围，具体地，所述先验信息包括初始断层长度和宽度。优选建议采用的一种较为合理的震源参数遍历范围如下：

顾及所述单一矩形模型参数反演的精度和效率，确定合理的震源参数遍历步长。优选建议采用的一种较为合理的震源参数遍历步长如下：

式中，ΔH、ΔL、ΔW、Δθ和Δk分别表示震源深度、断层长度和宽度、震源破裂方向以及震源相对位置的遍历步长。

第二步：根据所述震源参数遍历范围和所述震源参数遍历步长，生成震源参数解空间，并根据所述震源参数解空间、所述GNSS观测站位置信息、所述GNSS观测站地面峰值参数观测值以及所述震中位置，逐一为所述震源参数解空间中的每个震源参数可能解，计算GNSS观测站地面峰值参数理论值，具体通过以下流程实现：

首先：根据所述震源参数遍历范围和所述震源参数遍历步长，生成震源参数解空间。

其次：根据所述震源参数解空间、所述GNSS观测站位置信息以及所述震中位置，逐一为每个所述震源参数可能解，计算GNSS观测站断层距。根据GNSS观测站的位置，计算GNSS观测站断层距，分为如下两种情形：

(a)GNSS观测站位于矩形断层在地面上的投影：震源深度即为GNSS观测站断层距；

(b)GNSS观测站不在矩形断层在地面上的投影：根据下式计算GNSS观测站断层距：

式中，D_(i,j)表示GNSS观测站到线段(i,j)的最短距离，i和j为矩形断层两个相邻的端点(见图4)端点标记为1、2、3、4，min(D_(i,j))是从中取最小者；FD_h表示GNSS观测站断层距水平分量；H为震源深度；FD为GNSS观测站断层距。

然后：根据所述GNSS观测站断层距和所述GNSS观测站地面峰值速度(所述GNSS观测站地面峰值参数观测值)，通过使用最小二乘平差，拟合经验衰减公式的系数。简单、合理的经验衰减公式如下：

log10(PGV)＝a+b×log10(FD) (14)

式中，PGV为GNSS观测站地面峰值速度；FD为GNSS观测站断层距；a和b为需要拟合的衰减公式系数。

最后：根据所述GNSS观测站断层距以及所述经验衰减公式，计算GNSS观测站地面峰值速度理论值(所述GNSS观测站地面峰值参数理论值)。

第三步：根据所述GNSS观测站地面峰值速度、所述GNSS观测站地面峰值速度理论值以及目标函数，从所述震源参数解空间中确定震源参数最优解，具体地，实施例把所述震源参数最优解作为当前时刻单一矩形模型参数的反演解。使用的目标函数如下：

式中，表示GNSS观测站地面峰值速度；表示GNSS观测站地面峰值速度理论值；n为当前时刻可使用的GNSS观测站数量，下标i代表可使用的GNSS观测站序号；min表示最小值。

具体地，在实际工作中，至少3个地震台站发生地震触发，才开始计算地震震中；至少有2个GNSS观测站可用，才开始震源参数的反演；并且只要新增地震台站或GNSS观测站，立即重新估计地震震中或反演震源参数。

对于方法实施例，为了简单描述，故将其都表述为一系列的动作组合，但是本领域技术人员应该知悉，本发明实施例并不受所描述的动作顺序的限制，因为依据本发明实施例，某些步骤可以采用其他顺序或者同时进行。其次，本领域技术人员也应该知悉，说明书中所描述的实施例均属于优选实施例，所涉及的动作并不一定是本发明实施例所必须的。

图2示出了本发明一个实施例的基于单一矩形模型的地震震源快速确定系统的结构示意图。参照图2，本发明实施例的基于单一矩形模型的地震震源快速确定系统具体包括地震台站、GNSS观测站201、地面计算中心202、地面预警中心203，所述地震台站配置强震仪或地震仪，所述GNSS观测站配置GNSS接收机；

所述地震台站，用于采集地震台站原始运动信息，并将采集到的所述地震台站原始运动信息发送到地面计算中心202；

所述GNSS观测站，用于采集GNSS观测站基本导航定位信息，并将采集到的所述GNSS观测站基本导航定位信息发送到地面计算中心202；

所述地面计算中心202，用于根据本地存储的地震台站位置信息、本地存储的GNSS观测站位置信息以及所述地震台站原始运动信息、所述GNSS观测站基本导航定位信息，判断地震台站是否发生地震触发、确定地震震中、计算点源震级以及反演单一矩形模型参数，并将地震震中信息、点源震级信息以及单一矩形模型参数信息等地震预警信息发送到地震预警中心203；

所述地震预警中心203，用于根据接收到的所述地震预警信息，做出相应的地震预警举措。

具体实施时，地面计算中心和其他地震台站、GNSS观测站、地面预警中心之间分别建立通信。

本发明实施例中，如图3所示，所述地面计算中心202，具体包括接收单元2021、第一计算单元2022、第二计算单元2023、处理单元2024和发送单元2025，其中：

接收单元2021，用于接收所述地震台站原始运动信息、所述GNSS观测站基本导航定位信息；

第一计算单元2022，用于根据所述地震台站原始运动信息以及本地存储的所述地震台站位置信息，判断地震台站是否发生地震触发以及确定地震震中位置；

第二计算单元2023，用于根据所述GNSS观测站基本导航定位信息、本地存储的所述GNSS观测站位置信息以及所述地震震中位置，生成GNSS观测站地面峰值参数观测值以及估计地震点源震级；

处理单元2024，用于根据所述GNSS观测站地面峰值参数观测值、所述GNSS观测站位置信息、所述初始断层尺寸以及所述地震震中信息，快速反演单一矩形模型参数；

发送单元2025，用于把所述地震震中位置、所述地震点源震级以及所述单一矩形模型参数等地震预警信息发送到地震预警中心。

进一步地，所述第一计算单元2022，具体用于构建地震台站特征函数时间序列；判断地震台站是否发生地震触发；精确拾取地震台站初动时刻；确定地震震中位置。

进一步地，所述第二计算单元2023，具体用于生成GNSS观测站运动信息；提取GNSS观测站地面峰值参数观测值；计算GNSS观测站震中距；估计地震点源震级；确定单一矩形模型的初始断层尺寸。

具体实施时，接收单元2021、第一计算单元2022、第二计算单元2023、处理单元2024和发送单元2025可以依次连接。

对于系统实施例而言，由于其与方法实施例基本相似，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。

具体实施时，本发明技术方案提出的方法可由本领域技术人员采用计算机软件技术实现自动运行流程，实现方法的系统装置例如存储本发明技术方案相应计算机程序的计算机可读存储介质以及包括运行相应计算机程序的计算机设备，也应当在本发明的保护范围内。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。