技术领域
[0001] 本
发明涉及地震监测和地震预测领域，特别涉及一种快速获取输入地震事件的地震参数估计的方法和系统以及地震实时搜索引擎。
背景技术
[0002] 地震事件是地球内部在某个局部区域产生的介质构造突变而发生的，包括工业活动引发的微震事件，为增加石油与
天然气的产量而采用注
水破裂引起的微震事件，强度较小或者较大的自然地震事件，核爆炸事件或者是炸药爆炸事件。
[0003] 当一个地震事件产生时，该地震事件对应的
地震波形
信号数据会被地震监测台站记录下来(地震波接收仪器)，我们想要尽可能快的推断出该地震事件发生的
位置，
震级以及
震源机制等信息。目前该领域最为先进的技术能够利用纵波和/或横波的传播时间信息在地震发生后
马上推算出地震发生的位置，但是其震源机制信息却不得不花费数小时或者是数天的时间才能够通过分析得到。对于当前地震
定位技术和地震震级判断技术来说，仍然需要一定水平的专家来进行推算才能保证估计结果的准确程度。以发生在2004年12月26日苏
门答腊北部西海岸的9.1级的地震事件为例，因为该地震事件发生时是圣诞节假期的缘故，所以该地震事件没有被及时的报告，海啸的警报没有发出。实时地报告发生的地震事件的震源机制是很重要的。如Bouchon，M.，Karabulut，H.，Aktar，M.，Ozalaybey，S.，Schmittbuhl，J.，and Bouin，M-P.(2011).Extended nucleation of the 1999 Mw 7.6 Izmit Earthquake，Sciences 331，877-88分析了1999年土
耳其伊兹米特发生的7.6级地震的震前征兆，发现一系列前震的震源机制都显示相似的
断层错动，直至
加速到动态断裂。实时地获得前震的震源机制能够帮助我们开发出前期预警系统或者更进一步，甚至是预测地震的发生。然而，获得地震震源机制是需要一些诸如
波形模拟和波形反演等过程的处理，这一处理过程通常在地震发生之后的若干小时或者数天之内进行，并且这一过程通常要求使用高端的模拟与计算技术。
[0004] 需要一种快速获取发生的地震事件的地震参数信息的方法。
发明内容
[0005] 根据本发明
实施例，提出了一种获取地震事件的地震参数估计的方法，包括：
[0006] 输入地震事件的
地震图数据，地震图数据包括地震波形；
[0007] 在预先建立的历史地震图
数据库和
理论地震图数据库中通过波形匹配，采用近似最近邻搜索方法进行并行搜索，按照预定的相似度条件找出与该输入地震事件的地震图数据相似的多个地震图数据，作为搜索结果集合，其中所述历史地震图数据库存储实际发生的地震事件的地震图数据，所述理论地震图数据库存储人工模拟的地震图数据；以及[0008] 从搜索结果集合中确定匹配地震图数据，并从中获取输入地震事件的地震参数估计。
[0009] 根据本发明实施例，提出了一种获取地震事件的地震参数估计的系统，包括：
[0010] 历史地震图数据库，存储实际发生的地震事件的地震图数据；
[0011] 理论地震图数据库，存储人工模拟的地震图数据；
[0012] 输入装置，输入地震事件的地震图数据，地震图数据包括地震波形；
[0013] 搜索引擎，在历史地震图数据库和理论地震图数据库中通过波形匹配，采用近似最近邻方法进行并行搜索，按照预定的相似度条件找出与该输入地震事件的地震图数据相似的多个地震图数据，作为搜索结果集合，以及从搜索结果集合中确定匹配地震图数据，并从中获取输入地震事件的地震参数估计。
[0014] 根据本发明实施例，提出了一种地震图搜索引擎，包括：
[0015] 搜索装置，针对输入的地震事件的包括地震波形的地震图数据，在预先建立的历史地震图数据库和理论地震图数据库中通过波形匹配，采用近似最近邻方法进行并行搜索，按照预定的相似度条件找出与该输入地震事件的地震图数据相似的多个地震图数据，作为搜索结果集合，其中所述历史地震图数据库存储实际发生的地震事件的地震图数据，所述理论地震图数据库存储人工模拟的地震图数据；以及
[0016] 确定装置，在搜索结果集合中，确定匹配地震图数据，并从中获取输入地震事件的地震参数估计。
附图说明
[0017] 通过下面结合附图说明本发明的示例实施例，将使本发明的上述及其它目的、特征和优点更加清楚，其中：
[0018] 图1示出了三分量地震图数据的示例；
[0019] 图2示出了根据本发明实施例的获取地震事件的地震参数估计的系统的示意
框图；
[0020] 图3示出了建立历史地震图数据库的过程的示意
流程图；
[0021] 图4示出了建立理论地震图数据库的过程的示意流程图；
[0022] 图5示出了具有多个时间维度的地震图数据的示例；
[0023] 图6示出了建立多棵随机KD树(MRKD-Tree)结构的过程的示意图；
[0024] 图7示出了利用多棵随机KD树第一轮搜索过程；
[0025] 图8示出了利用多棵随机KD树第二轮搜索过程；
[0026] 图9示出了根据本发明实施例的地震事件搜索引擎的示意框图；
[0027] 图10示出了根据本发明实施例的获取地震事件的地震参数估计的方法的流程图；
[0028] 图11示出了根据本发明实施例的搜索结果的示例；
[0029] 图12示出了三种方法的地震图搜索时间的比较：精确线性搜索(Linear Search)，局部敏感哈希(LSH)方法，以及MRKD-Tree方法；以及
[0030] 图13示出了输入地震事件和使用线性搜索(Linear Search)和MRKD-Tree方法分别得到的最佳匹配地震事件之间的最大互相关系数的比较。
具体实施方式
[0031] 以下参照附图，对本发明的示例实施例进行详细描述，本发明不限于下述示例实施例。为了清楚描述本发明的基本思想，附图中仅示出了与本发明的技术方案密切相关的部件、功能或步骤，并且以下描述中省略了对已知技术、功能、部件或步骤的具体描述。
[0032] 地震事件是地球内部在某个局部区域产生的介质构造突变而发生的。地震事件是由地下的
岩石构造瞬间破裂所造成的。地震破裂机制可以用几个参数描述，例如破裂面的走向
角度，倾角角度，滑动角度，震级和
震源深度。考虑这些参数的结合，以及地震波记录的
震中距离，其相应的地震图有独特的波形记录，如横波，纵波，面波都有不同的特征。地震图是一个时间信号序列。对于给定的地球地质构造，每一个地震图对应着一个特定震源机制参数和震中距的组合。所以从地震图可以反推震源机制参数和震中距。从地震台或地震仪记录到的是地震波传播到记录地点引起的介质位移，或位移的速度(一阶时间导数)，或位移的加速度(二阶时间导数)。一般地震台可采用三分量(垂直，东西，南北方向震动记录)地震图数据，或一分量(垂直分量)地震图数据。图1显示的是三分量地震图数据的示意性示例，其中上面的图表示南北振动分量，中间的图表示东西振动分量，下面的图表示垂直振动分量，横轴为相对时间(秒)，纵轴为地震仪所在处地球介质振动的速度。
[0033] 在地震学领域，地震模拟中存在一个基本假设：如果同一个地震监测台站记录的两个地震图是相同的，那么这两个地震图对应的震源解以及它们传播过程的速度结构也应该是一样的。因此，在地震学领域，采用正演和反演方法来查找波形匹配的地震图是估计地震参数和地球结构的标准方法，例如可以参见Xu，Y.，Herrmann，R.B.，and Keith，D.K.(2010).Source parameters of regional snall-to-moderate earthquakes in the Yunnan-Sichuan region of China，Bull.Seismol.Soc.Am.100，No.5B，pp.2518-2531，以及Bonner，J.，Herrmann. R.，and Harley，B.(2010).Variable-period surface-wave magnitudes：a rapid and robust estimator of seismic moments.Bull.Seismol.Soc.Am.100，No.5A，pp.2301-2309。
[0034] 本发明实施例对于发生的地震事件，通过快速搜索历史记录的地震图数据和人工模拟的地震图数据，通过波形匹配找出与该地震事件的地震图数据相似的地震图数据，并从中获取该地震事件的地震参数估计。根据本发明实施例，历史记录和人工模拟的地震图数据是预先准备好的。此外，该搜索方法不仅仅是找出一个地震波形的解，还可以找出在一定相似度条件下不同的地震事件的波形的解集。根据本发明实施例，对于每个地震事件的波形的解就像是网络搜索引擎的结果中的每个网页链接一样，可以按相似度排序。这样使得有机会检测地震模拟中非唯一的地震事件的波形的解集。根据本发明实施例，采用近似最近邻搜索方法，能够在地震发生后几秒钟内，快速获得这些非唯一的地震事件的波形的解集，从而在极短的时间内对发生的地震事件有全面的认识。考虑到地震信号噪声的影响，模型参数的相互作用，以及地球结构的复杂性，需要获取多个非唯一的解，并且可以从这些解的特征评估搜索和解集的可靠性和解析度。
[0035] 图2示出了根据本发明实施例的获取地震事件的地震参数估计的系统的示意框图。该系统10包括：历史地震图数据库110，存储实际发生的地震事件的地震图数据；理论地震图数据库120，存储人工模拟的地震图数据；输入装置130，输入地震事件的地震图数据，地震图数据包括地震波形；搜索引擎140，在历史地震图数据库110和理论地震图数据库120中通过波形匹配，采用近似最近邻方法进行并行搜索，按照预定的相似度条件找出与该输入地震事件的地震图数据相似的多个地震图数据，作为搜索结果集合，以及从搜索结果集合中确定匹配地震图数据，并从中获取输入地震事件的地震参数估计。
[0036] 历史和理论地震图数据库可以是本地或远程存储的数据库，也可以是集中式或分布式的数据。
[0037] 根据本发明实施例，系统10还可以包括评估装置(未示出)，通过比较所有匹配地震图数据的互相关系数分布，来评估搜索结果的解析度和可靠性。
[0038] 地震图数据可以是如图1所示的三分量波形信号，也可以采用任何其他适合的形式，例如一分量垂直波形信号。地震图数据可以是数字化记录的，采用如标准的SAC，SEG2，SEGD或SEGY格式等存储，在数据文件头存储国际地震组织认定的地震参数，地震参数例如包括震中距、震源深度、震源机制和震级，其中震源机制可以包括走向角度、倾角角度和滑动角度。
[0039] 历史地震图数据库110
[0040] 历史地震图数据库110是预先建立的。为了建立历史地震图数据库110，可以利用单个或多个不同地震监测台站记录的地震事件的数字化记录。例如，可以针对某个特定的地区以及地震监测台网建立数据库，或对于大尺度的区域性以及超远地震距离的地震问题，根据震中距、震源深度和震源机制建立数据库。根据本发明实施例，可以在数据库中分开存储每个地震监测台的数字化记录。可以针对每个地震台获取地震参数估计，并将多个地震台地震参数估计的交集确定为输入地震事件的地震参数估计。根据本发明实施例，历史地震图数据库中的地震图数据是经过索引处理的，可以采用局部敏感哈希(LSH)方法进行搜索或者采用MRKD-Tree方法进行搜索。LSH是计算机科学领域的一种应用于文档检索和图片检索的常用的快速的搜索方法，在此省略对其的详细描述。详细方法描述参见Slaney，M.，and Casey，M.(2008)，Locality-sensitive hashing for finding nearest neighbors，IEEESingal Processing Magzine，March，pp128-131，图3示出了建立历史地震图数据库的一个示例的示意流程图300，在该示例中，应用MRKD-Tree方法，对地震图数据进行索引。如图3所示，在从单个或多个不同地震监测台站获得记录的地震图数据后，在步骤302，进行
信号处理、去噪，拾取P波初至时间。在步骤304，对地震图数据进行分段
频率滤波以分成多个频率分段，例如如选择频率分段[0.01Hz，0.04Hz]，[0.03Hz，0.06Hz]，并拾取滤波后的面波的最大振幅到时。这可以例如去掉不同台站数据中的仪器响应(参见Herrmann，R.(1973).Some aspects of band-pass filtering of surface
waves.Bull.Seismol.Soc.Am.63，No.2，pp.663-671)。在步骤306，针对每个频率分段，分别根据地震图数据的初至到时和最大振幅到时来对齐波形，截取得到两组地震图数据。在步骤308，分别对得到的两组地震图数据应用MRKD-Tree方法，以针对每一组地震图数据建立一棵或多棵KD树。因此，创建了按照MRKD-Tree结构进行索引的历史地震图数据库。根据本发明实施例，可以对具有这种MRKD-Tree结构的数据库应用MRKD-Tree方法与近似最近邻搜索方法的结合进行搜索。
[0041] 根据本发明实施例，地震图数据可以是多维度数据，这里“维度”可以是时间维度，或者也称为时间
采样。对于每一组地震图数据建立KD树可以包括：计算组内所有地震图数据在每个维度上的均值和方差，确定多个维度中具有相对于均值的最大方差的一个或多个维度，并且对于所确定的一个或多个维度中的每一个维度，以该维度上的均值作为中位数，将地震图数据划分成两部分，对每一个部分递归地进行计算、确定和划分步骤，直到每一部分只剩下单个地震图数据，从而建立一棵或多棵KD树。下面参照图5和6进一步描述如何对地震图数据应用MRKD-Tree方法的示例。图5显示的是一个数据库的10个地震图在3个最大方差分布的时间维度以及相应的均方差值(分别为1.67e-004，2.14e-004，1.85e-004)，其中地震图数据按初至到时对齐。本发明不限于此，可以具有多个时间维度。此外，可以在每一个时间维度上计算均方差，选出均方差最大的几个维度，作为地震图数据的时间维度。图6示出了建立MRKD-Tree结构的过程的示意图。如图所示，用三个地震图作数据库，每个地震图有四个维度。表格显示了地震图的幅值(地球介质振动速度)在每个维度的均值，以及其均方差。均方差的计算公式如下：
[0042]
[0043] 其中，Ai(t)是数据库里第i个地震图在维度t的振幅幅值，α(t)是所有地震图在维度t的振幅，N是数据库里地震图的个数。在图6的表格中，最大的均方差是在维度2和4，同为0.04667。那么可以同时在维度2和4开始建树。对于维度2，任何地震图的幅值在维度2小于均值0.3的放在左边，大于这个均值的放在右边，从而进行了划分。然后继续在左边或右边剩余地震图中应用该规则重复这个分树过程，直至最后只剩一个地震图。同时也在维度4开始建另外一棵树，依据同样的划分方法。由此，得到如图6下部所示的MRKD-Tree结构。
[0044] 综上，在针对每一组地震图数据建立MRKD-Tree结构时，根据方差取值最大的维度(时间)来将所有地震图数据分摊成两部分，即地震图数据空间被一个
正交超平面切分成了两个空间。这两部分的数据接着分别递归的进行分摊过程，即空间切分过程，最后创建了平衡二叉树。
[0045] 理论地震图数据库120
[0046] 实际的地震图数据在地球内部的
覆盖范围是有限制的。一个新发生的地震可能或者不可能与以往历史中的任何一个地震事件相似，特别是由于实际地震图的数据库不很完整，常有可能搜索不到非常相似的实际地震图。因此，本发明提出了采用理论模拟地震图数据作为补充，建立更大的地球内部覆盖范围。于是，根据本发明实施例，预先建立理论地震图数据库120。这样，对于新产生的地震事件，搜索过程是搜索历史数据库与理论数据库的并行过程。这两个并行搜索过程的返回结果都能够提供重要的独立信息，从而用于帮助估计当前地震事件的地震参数。根据本发明实施例，由于实际地震图数据是真实的地震记录，没有任何地球简化的假设，而理论地震图数据都是在各种简化的假设条件下计算出来的，因此优选地不将两种数据库混合在一起做搜索。当从历史地震图数据库和理论地震图数据库分别得到的匹配地震图数据与输入地震时间的地震图数据的互相关系数相同或相近时，搜索引擎优先选择历史地震图数据库的匹配地震图数据，来估计地震参数。
[0047] 为了准备一个充足的理论计算的地震图数据库，需要计算数百万的到数千万的人工模拟的地震图。地震图的特征，如横波和纵波的时间差，面波的
频谱特征，地震记录的时间系列长短等，随监测记录距离的不同而不同。如果目标是为了监测一个局部的区域，可以采用三维有限元或者有限差分方法，用三维地球速度模型，网格化的震源参数(如震中距和震源深度等)来计算短
波长(高频)的理论地震图。方法如Zhang，W .，and X.Chen.2006.Traction image method for irregular free surface boundaries in finite difference seismic wave simulation.Geophysical Journal International 167(1)：337-353.对于区域性的和超远距离的地震，可以采用弹性波模拟方法，使用全球平均的一维地球模型来计算长波长(低频)的地震图。一维模型如：Dziewonski，A.M.and Anderson，D.L.，(1981)，Preliminary reference earth model.Phys.Earth Planet.Inter.，25：297-356.无论距离远近，以上所有的地震图的波形数据都会用最初信号到时或者最大振幅到时进行对齐。此外，可以对特定的地震参数范围进行格点化，例如，对于震中距10°至35°(1°≈111公里)，深度0到50公里，按照距离间隔5公里，深度间隔5公里，走向角度间隔10°，倾角角度间隔10°，滑动角度间隔10°，并对每一个格点，人工模拟或合成理论地震图。同样，可以采用SAC格式存储理论地震图。
[0048] 图4示出了建立理论地震图数据库的一个示例的示意流程图400，在该示例中，也应用MRKD-Tree方法，对地震图数据进行索引。在步骤402，对特定的地震参数范围，进行格点化，并对每一个格点，计算理论地震图。在步骤404，将特定的地震台仪器频谱响应与理论地震图做卷积，使得理论数据与实际数据可比，因为实际接收的数据有地震台频谱响应的卷积效应。在步骤406，拾取P波初至到时。在步骤408，做分段频率滤波并拾取滤波后的最大振幅到时。在步骤410，针对每个频率分段，分别根据地震图数据的初至到时和最大振幅到时来对齐波形，截取得到两组地震图数据。在步骤412，分别对得到的两组地震图数据应用MRKD-Tree方法，以针对每一组地震图数据建立一棵或多棵KD树。因此，创建了按照MRKD-Tree结构进行索引的理论地震图数据库。在理论地震图数据库中，针对每一组地震图数据建立MRKD-Tree结构的过程与在历史地震图数据库中执行的过程相同。
[0049] 输入装置130
[0050] 输入装置130可以是用于输入地震事件的地震图数据的
接口，例如可以与单个或多个地震监测台的输出装置或
通信接口相连，以从单个或多个地震监测台接收产生的地震事件的地震图记录。输入装置130也可以包括例如
键盘、
触摸屏等用户接口，用户或操作者可以通过输入装置130输入系统设置、预定条件等，来控制和管理系统的运行。可以采用任何适当结构或形式来实现输入装置130。
[0051] 搜索引擎140
[0052] 搜索引擎140可以是分立实现的设备，也可以实现为与历史地震图数据库110、理论地震图数据库120和输入装置130中的一个或多个形成在同一位置。搜索引擎140可以利用
微处理器、
现场可编程门阵列或者专用计算设备来实现。
[0053] 一个充足的地震事件的数据库包含数以亿计到数以十亿计的地震图。即使是使用最先进的计算机设备，如果采用原始的线性搜索方法，也没有足够的处理能
力来查询如此大量的数据。诸如近似最近邻搜索(例如，参见Muja，M.，Lowe，D.G.(2009).Fast approximate nearest neighbors with automatic algorithm configuration，VISAPP(1)2009：331-340)的计算机搜索技术使得我们可以从大规模数据库中在一定准确率程度下进行快速的查找，从而得到相似的地震图，并且这一技术在互联网的数据量级下，已经成功地应用于图片检索，音乐检索以及视频检索领域中。
[0054] 搜索引擎140采用近似最近邻方法在历史地震图数据库110和理论地震图数据库120中进行并行搜索，按照预定的相似度条件找出与该输入地震事件的地震图数据相似的多个地震图数据，作为搜索结果集合。例如，预定的相似度条件可以是预定的搜索准确率，例如90％，凡是满足这一条件的地震图数据都可以作为搜索结果。根据本发明实施例，搜索引擎140在搜索过程中还进一步结合了局部敏感哈希(LSH)方法或者MRKD-Tree方法，从而进一步加快了搜索过程。图7示出了根据本发明实施例的搜索引擎140的示意框图。搜索引擎140包括搜索装置710和确定装置720。针对通过输入装置130输入的地震事件的包括地震波形的地震图数据，搜索装置710在历史和理论地震图数据库中通过波形匹配，采用近似最近邻方法进行并行搜索，按照预定的相似度条件找出与该输入地震事件的地震图数据相似的多个地震图数据，作为搜索结果集合。确定装置720在搜索结果集合中，确定匹配地震图数据，并从中获取输入地震事件的地震参数估计。根据本发明实施例，系统10还可以包括输出装置(未示出)，连接至搜索引擎140，输出搜索引擎140获取的地震参数估计。输出装置例如可以是显示屏、扬声器等。
[0055] 根据本发明实施例，搜索装置710可以在历史和理论地震图数据库中，针对多组地震图数据，并行地进行近似最近邻搜索，按照预定的相似度条件(例如，预定的搜索准确率，例如90％)，得到每组独立的多个搜索结果，并按照相似度从高到低得出一个整体的优先队列池，取出前面若干个(例如，20个)最近邻，作为最终搜索结果集合。这些搜索结果可以看成是该大规模数据库所有数据的一个子集。这里，相似度可以基于输入地震图的地震图数据与搜索结果集合中每个搜索结果之间的欧式距离，即，两者之间的差值。
[0056] 因为搜索到的地震图在预处理对齐时可能存在时间选择上的偏差，所以搜索结果根据欧式距离大小顺序进行排列就可能不是很准确。根据本发明实施例，相似度可以基于输入地震图的地震图数据与搜索结果集合中每个搜索结果之间的互相关系数。(归一化的)互相关系数可以如下计算：
[0057]
[0058] 其中，x(i)和y(i)是两个时间序列地震图，每个都有N个分立值(i＝1，2，3，...N)，mx是x(i)的平均值，my是y(i)的平均值。相关产生的时间序列r(j)(j＝1，2，...N)值应该在[-1，+1]之间。这个值越大，表明两个地震图越相似。虽然互相关系数的计算开销远大于欧式距离的计算开销，但是从相似度度量来看，其准确程度要欧式距离更为准确。由于只在近似最近邻的搜索结果里计算互相关系数，其数据规模要远小于整个数据库，因此计算互相关系数不会花费很多的时间。
[0059] 根据本发明实施例，确定装置720计算输入地震图的地震图数据与搜索结果集合中每个搜索结果之间的互相关系数，并采用互相关系数最大的搜索结果作为匹配地震图数据，从而确定匹配地震图数据。如果搜索装置710已按照互相关系数的大小形成了整体的优先队列池，则确定装置720可以直接从优先队列池中选择互相关系数最大的搜索结果，作为匹配地震图数据。当从历史地震图数据库和理论地震图数据库分别得到的匹配地震图数据与输入地震时间的地震图数据的互相关系数相同或相近(例如，相似度满足预定
阈值)时，确定装置720选择历史地震图数据库的匹配地震图数据。
[0060] 图9示出了根据本发明实施例的获取地震事件的地震参数估计的方法的流程图。方法800包括在步骤802，输入地震事件的地震图数据。地震图数据可以是在地震监测台记录的，并且是经过预处理的。预处理可以包括：进行信号处理、去噪，拾取P波初至时间；对地震图数据进行分段频率滤波以分成多个频率分段，并拾取滤波后的面波的最大振幅到时；
针对每个频率分段，分别根据地震图数据的初至到时和最大振幅到时来对齐波形，截取得到两组地震图数据。
[0061] 方法800还包括：在步骤804，在历史和理论地震图数据库中通过波形匹配，采用近似最近邻搜索方法进行并行搜索，按照预定的相似度条件找出与该输入地震事件的地震图数据相似的多个地震图数据，作为搜索结果集合。这里结合图6所示示例进行具体描述。如上所述，图6示意性示出了针对数据库中的每一组地震图数据建立MRKD-Tree结构。该图显示的是同时建立两棵树的例子。具体而言，在输入数据库地震图建树时，并行地分别判断每个地震图数据在其维度2和4的值是否大于或小于树顶的均值0.3和0.6。例如对于建立从维度2开始的树，任何地震图的幅值在维度2小于均值0.3的放在左边(ID＝1，2)，大于这个均值的放在右边(ID＝3)。放在左边的所有地震图(ID＝1，2)的最大均值是维度4的0.65。同样，在维度4，均值小于0.65的地震图放在
节点的左边(ID＝1)，均值大于0.65的地震图放在节点的右边(ID＝2)，直到到达了底层的最后一个地震图，这样便建立了一棵MRKD-Tree，至此树1建树结束。用同样的方式我们建立了右侧(树2)从维度4开始的树。
[0062] 图7和图8显示两轮搜索的过程。
[0063] 图7显示的是第一轮搜索过程，即在每棵MRKD-Tree中从根结点开始查询，具体在每棵MRKD-Tree中是按照如下方式进行查询的：首先，从根结点开始，比较根结点的“均值”与待查询地震图数据的“选择维度”的值，如果待查询地震图数据的“选择维度”的值小于等于根结点的“均值”，则进入根结点的左儿子结点，同时计算右儿子与待查询地震图之间的启发式距离信息(hdist)并将右儿子结点插入到优先队列池中，反之，如果大于根结点的“均值”，则进入根结点的右儿子结点，同时计算左儿子相应的启发式距离信息(hdist)并将左儿子结点插入到优先队列池中，如图7中的a)所示，进入树1的“选择维度＝4，均值＝0.65”的结点，并且计算树1的ID＝3结点的启发式距离信息并该结点插入到优先队列池，按照同样的处理方式，递归的比较内部结点的“均值”与待查询地震图在该结点对应的“选择维度”的值，进入左儿子或者右儿子结点，并且计算“未进入”的儿子结点的启发式距离信息并将其插入优先队列池中，直到到达叶子结点为止，此时将该叶子结点对应ID加入到结果队列中去，如图7中的b)所示；如果结果队列中已经包含了该叶子结点对应的ID，那么将跳过此次操作，不会重复插入相同的ID到结果队列中；至此完成了一棵树的第一轮搜索，按照同样的处理方式，搜索每棵MRKD-Tree，如图7中的c)所示，演示了树2的第一轮搜索结果。搜索完每棵MRKD-Tree标志着第一轮搜索的结束。
[0064] 图8显示的是从第二轮开始的搜索流程。首先，从优先队列池中取出队列首部的结点，即启发式距离信息(hdist)最小的结点；然后，根据此结点的信息，查找到对应的树标号，以及结点位置信息，从而，由此结点开始，按照第一轮的自上到下的方式，比较该结点的“均值”与待查询地震图在该结点对应“选择维度”的值，如果待查询地震图的“选择维度”的值小于等于该结点“均值”，则进入左儿子结点，并且计算其右儿子结点的启发式距离信息，将其加入优先队列池，若是大于该结点“均值”，则左右儿子互换，按照同样方式进行；如此反复，按照这种方式递归的进入儿子结点，直到到达叶子结点为止，然后将叶子结点加入结果队列，若队列中已经含有该结点，则不重复加入。这些操作都与第一轮中的相同。需要注意的是，如果我们
指定了结果队列大小的上限，那么结果队列中结点个数一旦达到这个上限，则整个搜索过程结束，否则继续按照如上方式进行查询，直到所有结点都加入结果队列为止。图8中的a)到d)每个图对应的是一次从优先队列池中取出结点并查询的过程。
[0065] 对于这种多棵KD
树搜索，原始的
基础的计算机理论与方法在Silpa-Anan，C.，and Hartley，R.(2008).Optimised KD-trees for fast image descriptor matching，In CVPR论文中有详尽描述。
[0066] 方法800还包括：在步骤806，从步骤804中找出的搜索结果集合中确定匹配地震图数据，并从中获取输入地震事件的地震参数估计。地震参数可以包括如下中的至少一个：震中距、震源深度、震源机制和震级。根据本发明实施例，方法800还可以包括：通过比较所有匹配地震图数据的互相关系数分布，来评估搜索结果的解析度和可靠性。
[0067] 下面描述利用本发明实施例的模拟示例，其中估计一个人工合成地震的地震参数。假定输入人工合成地震图的震中距300公里，震源深度25公里，震源机制：走向角度＝170°，倾角角度＝230°，滑动角度＝170°。对于这个示例，理论数据库可以采用如下参数来进行格点化：震中距100至1000km，以5km为间隔，震源深度5-30公里，以5km为间隔，震源机制：走向角度＝0至360°，倾角角度＝0至90°，滑动角度＝-180°至180°，间隔为10°。将人工合成地震图的波形输入根据本发明实施例的系统10，图11显示了一个搜索结果的示例，包括MRKD-Tree方法搜索的二十个最佳匹配结果。最上方的粗线的地震图是搜索出的与输入地震图完全相同的地震图，以及其它19个人工模拟的相似地震事件的地震图。
[0068] 图12示出了三种方法的地震图搜索时间的比较：精确的线性搜索(Linear Search)，局部敏感哈希(LSH)方法，以及MRKD-Tree方法。在搜索相同数量的地震图的数据库时，线性精确搜索返回100个最佳匹配结果，并且得到了100％的准确率和最长的搜索时间。LSH方法和MRKD-Tree方法得到的结果有90％的准确率，但是MRKD-Tree方法所用时间远远小于LSH方法和线性精确搜索方法。对于实际应用来说，90％的准确率已经足够。可见，MRKD-Tree方法的速度和效率非常高，根据实验结果，在数以百万计的数据量的数据库上，MRKD-Tree方法仅仅花费数秒时间就可以得到搜索结果，实现了快速甚至实时的地震参数估计。
[0069] 图13示出了输入地震事件和使用线性搜索(Linear Search)和MRKD-Tree方法分别得到的最佳匹配地震事件之间的最大互相关系数的比较。从图中可以看出，尽管MRKD-Tree搜索得到的结果中缺失了50个匹配地震图数据中的若干个，但是对于实际应用，MRKD-Tree方法已经能够搜索得到了足够多统计上需要的匹配地震图数据。考虑到MRKD-Tree方法相比于线性精确搜索的非常高的速度和效率，根据本发明实施例的应用MRKD-Tree方法的方法和系统具有几乎实时获取地震参数估计的优势。
[0070] 以上描述了根据本发明实施例的快速获取输入地震事件的地震参数估计的方法和系统以及地震时间搜索引擎，能够在得到输入地震数据后数秒钟内估计出该输入地震事件的地震参数。该方法和系统使用地震事件搜索引擎，通过建立大规模的历史实际地震资料数据库和建立能够完整覆盖所有震源参数的理论地震图的大规模数据库，快速从数据库中搜索出相似于输入地震图的所有地震事件。一旦搜索得到了输入地震事件对应的匹配结果，该系统立即显示出这些匹配地震事件的地震参数，以便获取例如震源机制和其他震源参数的估计。本发明实施例方法和系统能够利用计算机搜索技术，应用于天然地震的实时监测，实时监测石油与天然气开采生产所激发的微震事件，并且也可以用于实时监测核爆或者炸药爆炸等。
[0071] 以上所述是本发明的示例实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。
