CN113009559B - 基于多类型地球物理数据的地震评价方法 - Google Patents

基于多类型地球物理数据的地震评价方法 Download PDF

Info

Publication number
CN113009559B
CN113009559B CN202110228947.8A CN202110228947A CN113009559B CN 113009559 B CN113009559 B CN 113009559B CN 202110228947 A CN202110228947 A CN 202110228947A CN 113009559 B CN113009559 B CN 113009559B
Authority
CN
China
Prior art keywords
grid
data
seismic
target area
gridded
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN202110228947.8A
Other languages
English (en)
Other versions
CN113009559A (zh
Inventor
底青云
田飞
杨永友
王中兴
郑文浩
裴仁忠
张文秀
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Geology and Geophysics of CAS
Original Assignee
Institute of Geology and Geophysics of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Geology and Geophysics of CAS filed Critical Institute of Geology and Geophysics of CAS
Priority to CN202110228947.8A priority Critical patent/CN113009559B/zh
Publication of CN113009559A publication Critical patent/CN113009559A/zh
Application granted granted Critical
Publication of CN113009559B publication Critical patent/CN113009559B/zh
Priority to US17/679,222 priority patent/US11360225B1/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/01—Measuring or predicting earthquakes
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282—Application of seismic models, synthetic seismograms
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/301—Analysis for determining seismic cross-sections or geostructures
    • G01V1/302—Analysis for determining seismic cross-sections or geostructures in 3D data cubes
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/34—Displaying seismic recordings or visualisation of seismic data or attributes
    • G01V1/345—Visualisation of seismic data or attributes, e.g. in 3D cubes
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V11/00—Prospecting or detecting by methods combining techniques covered by two or more of main groups G01V1/00 - G01V9/00
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/64—Geostructures, e.g. in 3D data cubes

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供了一种基于多类型地球物理数据的地震评价方法。该方法包括：获取目标区域已发生板内天然地震数据、地表高程数据、布格重力异常数据、航磁异常数据、地震波层析成像、地壳GPS运动速率和震源机制解；构建目标区域平面网格，基于此将板内天然地震数据、地表高程数据、布格重力异常数据和航磁异常数据分别进行平面网格化处理；基于平面网格化的上述数据，确定目标区域强震区的平面分布图；将地表高程图、目标区域强震区的平面分布图与地震波层析成像进行叠加，结合地壳GPS运动速率和震源机制解，分析强震区两侧地壳的运动情况、分析强震区与地壳地质结构的匹配关系，从而确定强震带的地质成因。

Description

基于多类型地球物理数据的地震评价方法
技术领域
本发明属于板内天然地震评价技术领域，特别涉及一种基于多类型地球物理数据的地震评价方法。
背景技术
板内地震是由大陆板块运动引起的发生于板块内部的地震，其地点比较分散，频度较低，危害性大且震源机制复杂。近百年来的大陆地质、地震地质研究表明，看似稳定的大陆，内部不是完全的刚体、也不是完全的弹性体。板块间相对运动产生的、作用在板块边界的构造力被传递到板块内部的远程效应，以及多期构造运动形成的、现为板块内部不同地块之间分界线的先存构造的构造复活，都可导致现今看似稳定的大陆内部不同地块间存在着相对运动，诱发灾难性的地震。因此，板内地震引起巨大的破坏性，受到越来越多的重视。
尽管板内地震发生的并不频繁，但是在板块内部地震能量随距离的衰减通常很低，这会导致即使中等规模的地震也可以是毁灭性的。例如，美国新马德里(1811年)、印度的Latur(1992年)和Gujarat(2001年)等地震都发生在板块内部；特别是，2008年中国汶川7.9级板内地震导致大约70,000人死亡和总损失超过1000亿美元。然而，板内地震具有发生时间、发育深度和地震震级随机性强的特点，简单的投点显示只能定性描述地震的发生特征，基于现有板内地震数据的分布，很难确定哪些区域容易发生板内地震；并且同板间地震相比，板内地震发生在现今应变率非常低的地区，且深部先存断裂难以识别。因此板内地震的精确分析是一项很大的挑战。
发明内容
本发明的目的在于提供一种能够实现板内天然地震量化评价分析的方法。
为了实现上述目的，本发明提供了一种基于多类型地球物理数据的地震评价方法，其中，该方法包括：
获取目标区域已发生板内天然地震数据、地表高程数据、布格重力异常数据、航磁异常数据、地震波层析成像、地壳GPS运动速率和震源机制解；其中，地震波层析成像包括地震面波层析成像和/或地震体波层析成像；
将目标区域平面网格化，构建得到目标区域的平面网格；根据目标区域的平面网格，将板内天然地震数据、地表高程数据、布格重力异常数据和航磁异常数据分别进行平面网格化处理，得到平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据；
基于平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据，确定目标区域强震区的平面分布图(能够体现目标区域中强震区在平面上的分布特征的图，可以选用仅能体现强震区的平面分布特征的图也可以选用既能体现强震区也能体现其他类型震区的平面分布特征的图)；
将地表高程图、目标区域强震区的平面分布图与地震波层析成像进行叠加，结合地壳GPS运动速率和震源机制解，分析强震区两侧地壳的运动情况、分析强震区与地壳地质结构的匹配关系，从而确定强震带的地质成因。
上述基于多类型地球物理数据的地震评价方法基于多属性地质-地球物理数据对板内天然地震进行评价。
在上述基于多类型地球物理数据的地震评价方法中，优选地，所述将目标区域平面网格化，构建得到目标区域的平面网格包括：东西向按照0.01°-0.1°，南北向按照0.01°-0.1°建立网格；
例如，所述平面网格化处理包括：东西向按照0.05°(即宽度为5.55×cosθkm，θ为该点的纬度)，南北向按照0.05°(即长度5.55km)建立网格。
在上述基于多类型地球物理数据的地震评价方法中，优选地，根据目标区域的平面网格，将板内天然地震数据进行平面网格化处理包括：
基于目标区域已发生板内天然地震数据，确定目标区域中各网格的板内地震震级得到目标区域内板内地震震级分布，从而获得目标区域平面网格化的天然地震数据，实现定量表征目标区域内地震能量分布；其中，在确定目标区域中各网格的板内地震震级过程中，对于目标区域中每一个网格，将同一网格内的板内天然地震按照能量进行叠加并基于叠加所得的能量确定该网格的板内地震震级；
更优选地，当网格内叠加所得的能量为零时，将该网格叠加所得的能量定义为1。
在上述基于多类型地球物理数据的地震评价方法中，优选地，根据目标区域的平面网格，将地表高程数据进行平面网格化处理包括：
基于目标区域内地表高程数据，确定目标区域中各网格的地表高程值，从而获得目标区域平面网格化的地表高程数据；其中，在确定目标区域中各网格的地表高程值过程中，对于目标区域中每一个网格，将同一网格内的地表高程数据进行平均得到该网格的平均地表高程值作为该网格的地表高程值；
更优选地，当网格内没有地表高程值数据时，则取该网格东、西、南、北的4个邻近的网格的地表高程值计算几何平均值作为该网格的地表高程值。
在上述基于多类型地球物理数据的地震评价方法中，优选地，根据目标区域的平面网格，将布格重力异常数据进行平面网格化处理包括：
基于目标区域内地表高程数据，确定目标区域中各网格的布格重力异常值，从而获得目标区域平面网格化的布格重力异常数据；其中，在确定目标区域中各网格的布格重力异常值过程中，对于目标区域中每一个网格，将同一网格内的布格重力异常数据进行平均得到该网格的平均布格重力异常值作为该网格的布格重力异常值；
更优选地，当网格内没有布格重力异常数据时，则取该网格东、西、南、北的4个邻近的网格的布格重力异常值计算几何平均值作为该网格的布格重力异常值。
在上述基于多类型地球物理数据的地震评价方法中，优选地，根据目标区域的平面网格，将航磁异常数据进行平面网格化处理包括：
基于目标区域内地表高程数据，确定目标区域中各网格的航磁异常值，从而获得目标区域平面网格化的航磁异常数据；其中，在确定目标区域中各网格的航磁异常值过程中，对于目标区域中每一个网格，将同一网格内的航磁异常数据进行平均得到该网格的平均航磁异常值作为该网格的航磁异常值；
更优选地，当网格内没有航磁异常数据时，则取该网格东、西、南、北的4个邻近的网格的航磁异常值计算几何平均值作为该网格的航磁异常值。
在上述基于多类型地球物理数据的地震评价方法中，优选地，基于平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据，确定目标区域强震区的平面分布图包括：
对平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据分别进行数据归一化处理，得到归一化的平面网格化的天然地震数据、归一化的平面网格化的地表高程数据、归一化的平面网格化的布格重力异常数据和归一化的平面网格化的航磁异常数据；
通过主成分分析方法对归一化的平面网格化的天然地震数据、归一化的平面网格化的地表高程数据、归一化的平面网格化的布格重力异常数据和归一化的平面网格化的航磁异常数据进行降维，得到第一主成分数据和第二主成分数据；
对各组主成分数据进行聚类分析，确定强震区的第一主成分门槛值和第二主成分门槛值，将目标区域分为强震区和非强震区，得到目标区域强震区的平面分布图；
更优选地，在对各组主成分数据进行聚类分析过程中，进一步包括：确定中震区和弱震区的第一主成分门槛值和第二主成分门槛值，将目标区域非强震区进一步分为中震区、弱震区和无震区，进而得到目标区域不同分类震区的平面分布图。
在上述基于多类型地球物理数据的地震评价方法中，优选地，第一主成分数据＝0.91×归一化的平面网格化的天然地震数据-0.878×归一化地表高程+0.843×归一化的平面网格化的布格重力异常数据+0.835×归一化的平面网格化的航磁异常数据；
第二主成分数据＝0.081×归一化的平面网格化的天然地震数据+0.0395×归一化地表高程-0.011×归一化的平面网格化的布格重力异常数据+0.196×归一化的平面网格化的航磁异常数据；
在一具体实施方式中，强震区的第一主成分门槛值为1.52、第二主成分门槛值为0.91；即当第一主成分大于等于1.52且第二主成分大于等于0.91时为强震区。
在一具体实施方式中，对各组主成分数据进行聚类分析通过下述方式进行：利用K-means算法，对两组主成分数据的交会图进行无监督聚类。
在上述基于多类型地球物理数据的地震评价方法中，优选地，将地表高程图、目标区域强震区的平面分布图与地震波层析成像进行叠加包括：
基于目标区域强震区的平面分布图确定目标区域中的强震区分布；
基于目标区中的强震区分布，结合地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图；
将目标区域强震区的平面分布图投影到地表高程图上，获取地表高程与强震区平面分布的叠合图；
将所述地表高程与强震区平面分布的叠合图与所述地震波层析成像剖面图进行叠加。
在上述基于多类型地球物理数据的地震评价方法中，优选地，该方法还包括：将目标区域不同分类震区的平面分布图投影到目标区域的地表高程图上，得到目标区域地表高程与不同分类震区的叠合图。
在上述基于多类型地球物理数据的地震评价方法中，优选地，该方法还包括：根据目标区域强震区的平面分布图，确定经过目标区域强震区的剖面，并构建剖面网格；根据剖面网格，将板内天然地震数据进行剖面网格化处理，得到剖面网格化的天然地震数据(可以体现强震区的板内天然地震纵向分布情况)；
更优选地，所述构建剖面网格将板内天然地震数据进行剖面网格化处理包括：取剖面两侧各1-50km范围内的区域作为该剖面的控制范围，在该剖面上按照长度为1-10km、纵向深度为0.5-5km建立网格；例如，取剖面两侧各20km范围内的区域作为该剖面的控制范围，在该剖面上按照长度为5km、纵向深度为1km建立网格；
进一步优选地，根据剖面网格，将板内天然地震数据进行剖面网格化处理包括：
基于目标区域已发生天然地震记录，确定剖面控制范围内的天然地震记录，将剖面控制范围内的天然地震记录投点到剖面上，进而确定目标区域中各网格的板内地震震级；从而获得目标区域强震区剖面网格化的天然地震数据，实现定量表征目标区域内地震能量分布；其中，在确定目标区域中各网格的板内地震震级过程中，对于目标区域中每一个网格，将同一网格内的板内天然地震按照能量进行叠加并基于叠加所得的能量确定该网格的板内地震震级；
再优选地，当网格内叠加所得的能量为零时，将该网格叠加所得的能量定义为1。
在上述基于多类型地球物理数据的地震评价方法中，优选地，该方法还包括：
将目标区域三维网格化，构建得到目标区域的三维网格；根据目标区域的三维网格，将板内天然地震数据进行三维网格化处理，得到三维网格化的天然地震数据；
确定三维网格化的天然地震数据的强震区的截止值；
基于三维网格化的天然地震数据，利用三维网格化的天然地震数据的强震区的截止值，获取目标区域强震区的立体分布图(能够体现出强震区的立体分布特征)；
更优选地，该方法还包括：将目标区域强震区的立体分布图与地震波层析成像进行叠加，分析强震区与地壳地质结构的匹配关系。
在一具体实施方式中，所述目标区域三维网格化，构建得到目标区域的三维网格包括：东西向按照1-50km、南北向按照1-50km、纵向深度按照0.1-10km建立网格；例如，东西向按照20km、南北向按照20km、纵向深度按照1km建立网格。
在上述基于多类型地球物理数据的地震评价方法中，优选地，根据目标区域的三维网格，将板内天然地震数据进行三维网格化处理包括：基于目标区域已发生板内天然地震数据，确定目标区域中各网格的板内地震震级得到目标区域内板内地震震级分布，从而获得目标区域三维网格化的天然地震数据，实现定量表征目标区域内地震能量分布；其中，在确定目标区域中各网格的板内地震震级过程中，对于目标区域中每一个网格，将同一网格内的板内天然地震按照能量进行叠加并基于叠加所得的能量确定该网格的板内地震震级；
更优选地，当网格内叠加所得的能量为零时，将该网格叠加所得的能量定义为1；
在一具体实施方式中，三维网格化的天然地震数据的强震区的截止值为4。
在上述基于多类型地球物理数据的地震评价方法中，优选地，将目标区域强震区的立体分布图与地震波层析成像进行叠加包括：
基于目标区域强震区的立体分布图确定目标区域中的强震区分布；
基于目标区中的强震区分布，结合地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图；
将目标区域强震区的立体分布图与所述地震波层析成像剖面图进行叠加。
在上述基于多类型地球物理数据的地震评价方法中，优选地，基于目标区中的强震区分布，结合地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图包括：
将目标区域三维网格化，构建得到目标区域的三维网格；根据目标区域的三维网格，将地震波层析成像进行三维网格化处理，得到三维网格化的地震波层析成像；
基于目标区中的强震区分布，结合三维网格化的地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图；
更优选地，根据目标区域的三维网格，将地震波层析成像进行三维网格化处理包括：基于目标区域内地震波层析成像数据，确定目标区域中各网格的地震波层析成像值，从而获得目标区域三维网格网格化的地震波层析成像数据；其中，在确定目标区域中各网格的地震波层析成像值过程中，对于目标区域中每一个网格，将同一网格内的地震波层析成像数据进行平均得到该网格的平均地震波层析成像值作为该网格的地震波层析成像值；
进一步优选地，当网格内没有地震波层析成像数据时，则取该网格上、下、左、右、前、后6个邻近的网格的地震波层析成像值计算几何平均值作为该网格的地震波层析成像值。
在一具体实施方式中，所述目标区域三维网格化，构建得到目标区域的三维网格包括：东西向按照1-50km、南北向按照1-50km、纵向深度按照0.1-10km建立网格。
在上述基于多类型地球物理数据的地震评价方法中，优选地，所述并基于叠加所得的能量确定该网格的板内地震震级通过下述公式实现：
lgE＝4.8+1.5M
式中，E为叠加所得的能量，J；M为地震里氏震级。
本发明提供的技术方案基于多属性地质-地球物理数据进行板内天然地震定量评价、级别划分、成因分析。具体而言，该方法将地震能量分布、地表高程、布格重力异常和航磁异常有效结合，进行震区级别划分，确定地震发育带；进一步与地震波层析成像、地壳GPS运动速率和震源机制解相结合，进行成因分析；有效实现了板内天然地震量化评价分析，为板内地震的精确分析提供了一种可能，推动了板内地震预测工作的进行。
附图说明
为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要使用的附图作简单的介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图：
图1为本发明一实施例提供的基于多类型地球物理数据的地震评价方法的流程示意图。
图2为本发明实施例1中目标区域内天然地震数据投影图。
图3A为本发明实施例1中目标区域内平面网格化的地震震级平面分布图。
图3B为本发明实施例1中目标区域内平面网格化的地表高程数据图。
图3C为本发明实施例1中目标区域内平面网格化的布格重力异常数据图。
图3D为本发明实施例1中目标区域内平面网格化的航磁异常数据图。
图4为本发明实施例1中目标区域不同分类震区的平面分布图。
图5为本发明实施例1中目标区域强震区的立体分布图。
图6A为本发明实施例1中A-A’剖面位置图。
图6B为本发明实施例1中A-A’剖面地震记录图。
图6C为本发明实施例1中A-A’剖面地震震级剖面分布图。
图7A为本发明实施例1中强震区平面分布与成因分析图。
图7B为本发明实施例1中强震区立体分布与成因分析图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚完整的描述。显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明的保护范围。
下面参考本发明的若干代表性实施方式，详细阐述本发明的原理和精神。
参见图1，本发明的一实施例提供了一种基于多类型地球物理数据的地震评价方法，该方法包括：
步骤S1：获取目标区域已发生板内天然地震数据、地表高程数据、布格重力异常数据、航磁异常数据、地震波层析成像、地壳GPS运动速率和震源机制解；其中，地震波层析成像包括地震面波层析成像和/或地震体波层析成像；
步骤S2：将目标区域平面网格化，构建得到目标区域的平面网格；根据目标区域的平面网格，将板内天然地震数据、地表高程数据、布格重力异常数据和航磁异常数据分别进行平面网格化处理，得到平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据；
步骤S3：基于平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据，确定目标区域强震区的平面分布图；
步骤S4：将地表高程图、目标区域强震区的平面分布图与地震波层析成像进行叠加，结合地壳GPS运动速率和震源机制解，分析强震区两侧地壳的运动情况、分析强震区与地壳地质结构的匹配关系，从而确定强震带的地质成因。
其中，目标区域强震区的平面分布图指能够体现目标区域强震区在平面上的分布特征的图；可以选用仅能体现强震区的平面分布特征的图，也可以选用既能体现强震区也能体现其他类型震区的平面分布特征的图(例如目标区域不同分类震区的平面分布图)。
其中，板内天然地震数据：可以源自世界地震中心(ISC)报告、美国哈佛大学和美国地质调查局的CMT解、和/或、中国地震局及其各省局解析处理结果。
其中，地表高程数据：地表高程可以反映地表上固定性物体的变化起伏；根据不同的等高线曲线形状，可以实现表面形貌的可视化。地表高程数据可以直接来自于美国地球物理中心(U.S.National Geophysical Data Center，NGDC)发布的全球地形起伏模型ETOPO；该模型整合了众多全球和区域数据集，根据陆地地形和海洋测深法构建而成，采用UTM/WGS84坐标系统，精度为1′×1′。
其中，布格重力异常数据：布格重力异常是经过纬度改正、地形改正及布格改正后获得的异常，主要是反映地球内部异常质量对重力测量结果的影响，其大范围内的变化主要反映的是莫霍界面的起伏。布格重力异常数据可以收集自中国地质调查局的全国区域重力数据库，数据精度为±2×10-5m/s2。
其中，航磁异常数据：岩石圈磁场主要反映岩石圈内岩石磁化强度(磁化率或磁矩)的分布；由于岩石的磁化环境(指岩石被磁化时的外部磁场及其随时间的变化)、磁性载体及其经历的地质作用与构造演化存在差异，可以得到岩石圈磁场携带着岩石矿物成分、岩石结构、温度与压力状态、深部岩石变质作用及其经历的构造运动与演化过程等地质构造与地球动力学信息。航磁异常数据可以收集自中国国土资源部航空物探遥感中心，数据的网格为2km×2km。
其中，地震面波层析成像：利用背景噪声研究地下结构，即对两个台站长时间的地震噪声记录进行互相关计算，提取台站间的格林函数，获取面波频散特征，并进一步层析成像或测量面波走时变化，可获得壳幔速度结构和对地下介质监测。Shen et al.(Shen etal.,2016)利用分布在中国各地(CEArray，China Array，NECESS，PASSCAL，GSN)及周边地区(Korean earthquake Network，F-Net，KNET)的2000多个地震台站的数据，绘制了覆盖整个中国，向东延伸穿过朝鲜半岛，进入边缘海域，直到日本的8-50s周期的各向同性瑞利波群和相位速度图；这套模型在地壳和地幔岩石圈中都显示出极大的多样性和相当丰富性的速度结构，例如清晰的显示了华南陆块四川盆地的低速异常及其下方的高速异常。本发明使用的地震面波层析成像可以直接采用Shen et al.(Shen et al.,2016)公开的中国面波层析成像的三维数据，数据采样为0.5°×0.5°，垂向采样间隔为1km，纵深从地表到地下150km深度。
其中，地震体波层析成像：地震层析成像使用地震数据逐层剖析法绘制地球内部图像，反演地球内部地质体的速度结构；地质体速度结构的详细信息可以为地球深部地质结构和动力学演化提供一个重要的视角；高分辨率体波层析成像能够研究地球地幔结构，显示更详细的上地幔结构信息，可以更好地反映地球内部的变形和整体地球动力学演化；地震体波层析成像纵深可以达到从地表到地下1000km深度。
其中，地壳GPS运动速率：GPS连续观测资料使用精度高达亚厘米甚至毫米量级的观测数据，可提取地壳形变信息，分析地壳位移、应力、应变场等动力学特征，为定量研究地壳现今水平运动和构造变形提供了革命性的技术手段；例如，Wang and Shen(Wang andShen，2020)对中国境内近25年的GPS观测资料进行了严格的处理，得到了站点的长期速度，得到高空间分辨率和一致性的速度解，能够揭示青藏高原向西挤压、青藏高原东南边界的旋转和华南陆块内部变形等模式，可以获得震前、同震及震后形变，为板内地震研究提供重要依据。本发明使用的地壳GPS运动速率可以直接采用Wang and Shen(Wang and Shen，2020)提供的这套中国境内最新的GPS数据，包含2403个数据点，精度为0.1mm/yr。
其中，震源机制解：震源机制描述了破裂过程的性质和来源，是了解震源区和发震断层应力状态的重要基础；当断层相关事件发生时，震源机制是指滑移的断层平面的方向和滑移矢量，也称为断层平面解。本发明使用的震源机制解可以采用世界地震中心(ISC)报告、美国哈佛大学和美国地质调查局、中国地震局及其各省局解析处理的震源机制解。
其中，分析强震区两侧地壳的运动情况、分析强震区与地壳地质结构的匹配关系，从而确定强震带的地质成因可以通过下述方式实现：
通过强震区分布带及其两侧的地表高程变化以及地震波层析成像变化情况，确定强震区的地壳地质结构从而得到强震区与地壳地质结构的匹配关系；结合地壳GPS运动速率确定强震区两侧地壳的运动情况；基于强震区与地壳地质结构的匹配关系、强震区两侧地壳的运动情况结合震源机制解，确定强震带的地质成因。
进一步，所述将目标区域平面网格化，构建得到目标区域的平面网格包括：东西向按照0.01°-0.1°，南北向按照0.01°-0.1°建立网格；例如，所述平面网格化处理包括：东西向按照0.05°(即宽度为5.55×cosθkm，θ为该点的纬度)，南北向按照0.05°(即长度5.55km)建立网格。
进一步，根据目标区域的平面网格，将板内天然地震数据进行平面网格化处理包括：
基于目标区域已发生板内天然地震数据，确定目标区域中各网格的板内地震震级得到目标区域内板内地震震级分布，从而获得目标区域平面网格化的天然地震数据，实现定量表征目标区域内地震能量分布；其中，在确定目标区域中各网格的板内地震震级过程中，对于目标区域中每一个网格，将同一网格内的板内天然地震按照能量进行叠加并基于叠加所得的能量确定该网格的板内地震震级；
更进一步，当网格内叠加所得的能量为零时，将该网格叠加所得的能量定义为1。
进一步，根据目标区域的平面网格，将地表高程数据进行平面网格化处理包括：
基于目标区域内地表高程数据，确定目标区域中各网格的地表高程值，从而获得目标区域平面网格化的地表高程数据；其中，在确定目标区域中各网格的地表高程值过程中，对于目标区域中每一个网格，将同一网格内的地表高程数据进行平均得到该网格的平均地表高程值作为该网格的地表高程值；
更进一步，当网格内没有地表高程值数据时，则取该网格东、西、南、北的4个邻近的网格的地表高程值计算几何平均值作为该网格的地表高程值。
进一步，根据目标区域的平面网格，将布格重力异常数据进行平面网格化处理包括：
基于目标区域内地表高程数据，确定目标区域中各网格的布格重力异常值，从而获得目标区域平面网格化的布格重力异常数据；其中，在确定目标区域中各网格的布格重力异常值过程中，对于目标区域中每一个网格，将同一网格内的布格重力异常数据进行平均得到该网格的平均布格重力异常值作为该网格的布格重力异常值；
更进一步，当网格内没有布格重力异常数据时，则取该网格东、西、南、北的4个邻近的网格的布格重力异常值计算几何平均值作为该网格的布格重力异常值。
进一步，根据目标区域的平面网格，将航磁异常数据进行平面网格化处理包括：
基于目标区域内地表高程数据，确定目标区域中各网格的航磁异常值，从而获得目标区域平面网格化的航磁异常数据；其中，在确定目标区域中各网格的航磁异常值过程中，对于目标区域中每一个网格，将同一网格内的航磁异常数据进行平均得到该网格的平均航磁异常值作为该网格的航磁异常值；
更进一步，当网格内没有航磁异常数据时，则取该网格东、西、南、北的4个邻近的网格的航磁异常值计算几何平均值作为该网格的航磁异常值。
进一步，步骤S3包括：
对平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据分别进行数据归一化处理，得到归一化的平面网格化的天然地震数据、归一化的平面网格化的地表高程数据、归一化的平面网格化的布格重力异常数据和归一化的平面网格化的航磁异常数据；
通过主成分分析方法对归一化的平面网格化的天然地震数据、归一化的平面网格化的地表高程数据、归一化的平面网格化的布格重力异常数据和归一化的平面网格化的航磁异常数据进行降维，得到第一主成分数据和第二主成分数据；
对各组主成分数据进行聚类分析，确定强震区的第一主成分门槛值和第二主成分门槛值，将目标区域分为强震区和非强震区，得到目标区域强震区的平面分布图；
更进一步，在对各组主成分数据进行聚类分析过程中，进一步包括：确定中震区和弱震区的第一主成分门槛值和第二主成分门槛值，将目标区域非强震区进一步分为中震区、弱震区和无震区，进而得到目标区域不同分类震区的平面分布图；
此时，步骤S3包括：
对平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据分别进行数据归一化处理，得到归一化的平面网格化的天然地震数据、归一化的平面网格化的地表高程数据、归一化的平面网格化的布格重力异常数据和归一化的平面网格化的航磁异常数据；
通过主成分分析方法对归一化的平面网格化的天然地震数据、归一化的平面网格化的地表高程数据、归一化的平面网格化的布格重力异常数据和归一化的平面网格化的航磁异常数据进行降维，得到第一主成分数据和第二主成分数据；
对各组主成分数据进行聚类分析，确定强震区、中震区和弱震区的第一主成分门槛值和第二主成分门槛值，将目标区域分为强震区、中震区、弱震区和无震区四类得到目标区域不同分类震区的平面分布图(能够体现出不同分类震区在平面的分布特征，包括强震区在平面的分布特征)；
其中，对各组主成分数据进行聚类分析优选通过下述方式进行：利用K-means算法，对两组主成分数据的交会图进行无监督聚类；
其中，归一化处理为将数据归一化至0-1范围内；
其中，确定的到的第一主成分、第二主成分的表达式优选为：
第一主成分数据＝0.91×归一化的平面网格化的天然地震数据-0.878×归一化地表高程+0.843×归一化的平面网格化的布格重力异常数据+0.835×归一化的平面网格化的航磁异常数据；
第二主成分数据＝0.081×归一化的平面网格化的天然地震数据+0.0395×归一化地表高程-0.011×归一化的平面网格化的布格重力异常数据+0.196×归一化的平面网格化的航磁异常数据；
在一具体实施方式中，强震区的第一主成分门槛值为1.52、第二主成分门槛值为0.91；即当第一主成分大于等于1.52或者第二主成分大于等于0.91时为强震区；
在一具体实施方式中，强震区的第一主成分门槛值为1.52、第二主成分门槛值为0.91；中震区的第一主成分门槛值为0.47、第二主成分门槛值为0.43；弱震区的第一主成分门槛值为-1.32、第二主成分门槛值为-0.21；即当第一主成分大于等于1.52或者第二主成分大于等于0.91时为强震区；当第一主成分大于等于0.47小于1.52或者第二主成分大于等于0.43小于0.91时为中震区；当第一主成分大于等于-1.32小于0.47或者第二主成分大于等于-0.21小于0.43时为中震区。
进一步，将地表高程图、目标区域强震区的平面分布图与地震波层析成像进行叠加包括：
基于目标区域强震区的平面分布图确定目标区域中的强震区分布；
基于目标区中的强震区分布，结合地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图；
将目标区域强震区的平面分布图投影到地表高程图上，获取地表高程与强震区平面分布的叠合图
将所述地表高程与强震区平面分布的叠合图与所述地震波层析成像剖面图进行叠加。
进一步，该方法还包括：将目标区域不同分类震区的平面分布图投影到目标区域的地表高程图上，得到目标区域地表高程与不同分类震区的叠合图。
进一步，该方法还包括：根据目标区域强震区的平面分布图，确定经过目标区域强震区的剖面，并构建剖面网格；根据剖面网格，将板内天然地震数据进行剖面网格化处理，得到剖面网格化的天然地震数据(可以体现强震区的板内天然地震纵向分布情况)；
更进一步，所述构建剖面网格将板内天然地震数据进行剖面网格化处理包括：取剖面两侧各1-50km范围内的区域作为该剖面的控制范围，在该剖面上按照长度为1-10km、纵向深度为0.5-5km建立网格；例如，取剖面两侧各20km范围内的区域作为该剖面的控制范围，在该剖面上按照长度为5km、纵向深度为1km建立网格；
更进一步，根据剖面网格，将板内天然地震数据进行剖面网格化处理包括：
基于目标区域已发生天然地震记录，确定剖面控制范围内的天然地震记录，将剖面控制范围内的天然地震记录投点到剖面上，进而确定目标区域中各网格的板内地震震级；从而获得目标区域强震区剖面网格化的天然地震数据，实现定量表征目标区域内地震能量分布；其中，在确定目标区域中各网格的板内地震震级过程中，对于目标区域中每一个网格，将同一网格内的板内天然地震按照能量进行叠加并基于叠加所得的能量确定该网格的板内地震震级；
更进一步，当网格内叠加所得的能量为零时，将该网格叠加所得的能量定义为1。
进一步，该方法还包括：
将目标区域三维网格化，构建得到目标区域的三维网格；根据目标区域的三维网格，将板内天然地震数据进行三维网格化处理，得到三维网格化的天然地震数据；
确定三维网格化的天然地震数据的强震区的截止值；
基于三维网格化的天然地震数据，利用三维网格化的天然地震数据的强震区的截止值，获取目标区域强震区的立体分布图(能够体现出强震区的立体分布特征)；
更进一步，该方法还包括：将目标区域强震区的立体分布图与地震波层析成像进行叠加，分析强震区与地壳地质结构的匹配关系；
更进一步，所述目标区域三维网格化，构建得到目标区域的三维网格包括：东西向按照1-50km、南北向按照1-50km、纵向深度按照0.1-10km建立网格；例如，东西向按照20km、南北向按照20km、纵向深度按照1km建立网格；
更进一步，根据目标区域的三维网格，将板内天然地震数据进行三维网格化处理包括：基于目标区域已发生板内天然地震数据，确定目标区域中各网格的板内地震震级得到目标区域内板内地震震级分布，从而获得目标区域三维网格化的天然地震数据，实现定量表征目标区域内地震能量分布；其中，在确定目标区域中各网格的板内地震震级过程中，对于目标区域中每一个网格，将同一网格内的板内天然地震按照能量进行叠加并基于叠加所得的能量确定该网格的板内地震震级；
更进一步，当网格内叠加所得的能量为零时，将该网格叠加所得的能量定义为1；
更进一步，三维网格化的天然地震数据的强震区的截止值为4。
进一步，将目标区域强震区的立体分布图与地震波层析成像进行叠加包括：
基于目标区域强震区的立体分布图确定目标区域中的强震区分布；
基于目标区中的强震区分布，结合地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图；
将目标区域强震区的立体分布图与所述地震波层析成像剖面图进行叠加。
进一步，基于目标区中的强震区分布，结合地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图包括：
将目标区域三维网格化，构建得到目标区域的三维网格；根据目标区域的三维网格，将地震波层析成像进行三维网格化处理，得到三维网格化的地震波层析成像；
基于目标区中的强震区分布，结合三维网格化的地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图；
更进一步，所述目标区域三维网格化，构建得到目标区域的三维网格包括：东西向按照1-50km、南北向按照1-50km、纵向深度按照0.1-10km建立网格；
更进一步，根据目标区域的三维网格，将地震波层析成像进行三维网格化处理包括：基于目标区域内地震波层析成像数据，确定目标区域中各网格的地震波层析成像值，从而获得目标区域三维网格网格化的地震波层析成像数据；其中，在确定目标区域中各网格的地震波层析成像值过程中，对于目标区域中每一个网格，将同一网格内的地震波层析成像数据进行平均得到该网格的平均地震波层析成像值作为该网格的地震波层析成像值；再进一步，当网格内没有地震波层析成像数据时，则取该网格上、下、左、右、前、后6个邻近的网格的地震波层析成像值计算几何平均值作为该网格的地震波层析成像值。
进一步，基于叠加所得的能量确定该网格的板内地震震级通过下述公式实现：
lgE＝4.8+1.5M
式中，E为叠加所得的能量，J；M为地震里氏震级。
实施例1
本实施例提供了一种基于多类型地球物理数据的地震评价方法，该方法用以对中国境内板内天然地震进行分析评价。
具体而言，本实施例提供的基于多类型地球物理数据的地震评价方法包括：
1、获取中国境内已发生板内天然地震数据、地表高程数据、布格重力异常数据、航磁异常数据、地震面波层析成像、地壳GPS运动速率和震源机制解；
其中，中国已发生板内天然地震数据采集自世界地震中心(ISC)报告、美国哈佛大学和美国地质调查局的CMT解、和中国地震局及其各省局解析处理结果。总共收集了自1965以来在中国境内发生的65万个地震记录(如图2所示)；
地表高程数据采集自美国地球物理中心(U.S.National Geophysical DataCenter，NGDC)发布的全球地形起伏模型ETOPO1；该模型整合了众多全球和区域数据集，根据陆地地形和海洋测深法构建而成，采用UTM/WGS84坐标系统，精度为1′×1′；
布格重力异常数据采集自中国境内地区，重力数据源来自中国地质调查局的全国区域重力数据库，数据精度为±2×10-5m/s2；
航磁异常数据采集自中国国土资源部航空物探遥感中心，数据的网格为2km×2km；
地震面波层析成像采集自Shen et al.(Weisen Shen,Michael H.Ritzwoller,Dou Kang,YoungHee Kim,Fan-Chi Lin,Jieyuan Ning,Weitao Wang,Yong Zheng,Longquan Zhou.A seismic reference model for the crust and uppermost mantlebeneath China from surface wave dispersion.Geophysical Journal International,2016,206(2):954-979)公开的中国面波层析成像的三维数据，数据采样为0.5°×0.5°，垂向采样间隔为1km，纵深从地表到地下150km深度；
地壳GPS运动速率采集自Wang and Shen(Min Wang,Shen,Zheng-KangShen.Present-Day Crustal Deformation of Continental China Derived From GPSand Its Tectonic Implications.Journal of Geophysical Research:Solid Earth,2020,125(2):1-22)提供的中国境内最新的GPS数据，包含2403个数据点，精度为0.1mm/yr；
震源机制解采集自世界地震中心(ISC)报告、美国哈佛大学和美国地质调查局、中国地震局及其各省局解析处理的震源机制解，获得了中国境内3115个震源机制解结果。
2、平面网格化处理：
构建中国境内的平面网格：东西向按照0.05°(即宽度为5.55×cosθkm，θ为该点的纬度)，南北向按照0.05°(即长度5.55km)建立网格；
根据中国境内的平面网格，将板内天然地震数据、地表高程数据、布格重力异常数据和航磁异常数据分别进行平面网格化处理，得到平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据；
其中，将板内天然地震数据进行平面网格化处理包括：基于中国境内已发生板内天然地震记录，确定中国境内的各网格的板内地震震级得到中国境内板内地震震级平面分布；其中，在确定各网格的板内地震震级过程中，对于每一个网格，将同一网格内的地震记录按照能量进行叠加并基于叠加所得的能量确定该网格的板内地震震级；其中，当网格内叠加所得的能量为零时，将该网格叠加所得的能量定义为1；地震震级通过下述公式确定：lgE＝4.8+1.5M；式中，E为叠加所得的能量，J；M为地震里氏震级；结果如图3A所示；
地表高程数据进行平面网格化处理通过下述方式实现：对于中国境内的每一个网格，将同一网格内的地表高程数据进行平均得到该网格的平均地表高程值作为该网格的地表高程值，从而获得中国境内平面网格化的地表高程数据；其中，当网格内没有地表高程值数据时，则取该网格东、西、南、北的4个邻近的网格的地表高程值计算几何平均值作为该网格的地表高程值；结果如图3B所示；
布格重力异常数据进行平面网格化处理通过下述方式实现：对于中国境内的每一个网格，将同一网格内的布格重力异常数据进行平均得到该网格的平均布格重力异常值作为该网格的布格重力异常值，从而获得中国境内平面网格化的布格重力异常数据；其中，当网格内没有布格重力异常值数据时，则取该网格东、西、南、北的4个邻近的网格的布格重力异常值计算几何平均值作为该网格的布格重力异常值；结果如图3C所示；
航磁异常数据进行平面网格化处理通过下述方式实现：对于中国境内的每一个网格，将同一网格内的航磁异常数据进行平均得到该网格的平均航磁异常值作为该网格的航磁异常值，从而获得中国境内平面网格化的航磁异常数据；其中，当网格内没有航磁异常值数据时，则取该网格东、西、南、北的4个邻近的网格的航磁异常值计算几何平均值作为该网格的航磁异常值；结果如图3D所示。
3、对平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据分别进行数据归一化处理，得到归一化的平面网格化的天然地震数据、归一化的平面网格化的地表高程数据、归一化的平面网格化的布格重力异常数据和归一化的平面网格化的航磁异常数据；
其中，归一化处理采用如下公式进行，将平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据分别归一化至0-1范围内；
Zi归一化＝(Zi-Zmin)÷(Zmax-Zmin)；
式中，Zi归一化为Zi归一化后的数据；Zi为第i个Z数据；Zmax为Z数据最大值；Zmin为Z数据最小值；其中，Z数据为平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据或平面网格化的航磁异常数据。
通过主成分分析方法对归一化的平面网格化的天然地震数据、归一化的平面网格化的地表高程数据、归一化的平面网格化的布格重力异常数据和归一化的平面网格化的航磁异常数据进行降维，得到第一主成分数据Y1和第二主成分数据Y2；
得到的第一主成分数据＝0.91×归一化的平面网格化的天然地震数据-0.878×归一化地表高程+0.843×归一化的平面网格化的布格重力异常数据+0.835×归一化的平面网格化的航磁异常数据；第二主成分数据＝0.081×归一化的平面网格化的天然地震数据+0.0395×归一化地表高程-0.011×归一化的平面网格化的布格重力异常数据+0.196×归一化的平面网格化的航磁异常数据。
利用K-means算法，对Y1和Y2的交会图进行无监督聚类，确定强震区、中震区和弱震区的Y1门槛值和Y2门槛值，将目标区域分为强震区、中震区、弱震区和无震区四类得到目标区域不同分类震区的平面分布图；
确定的强震区的第一主成分门槛值为1.52、第二主成分门槛值为0.91；中震区的第一主成分门槛值为0.47、第二主成分门槛值为0.43；弱震区的第一主成分门槛值为-1.32、第二主成分门槛值为-0.21；即当第一主成分大于等于1.52或者第二主成分大于等于0.91时为强震区；当第一主成分大于等于0.47小于1.52或者第二主成分大于等于0.43小于0.91时为中震区；当第一主成分大于等于-1.32小于0.47或者第二主成分大于等于-0.21小于0.43时为中震区；
目标区域不同分类震区的平面分布图如图4所示。
4、三维网格化处理：
构建中国境内的三维网格：东西向按照20km、南北向按照20km、纵向深度按照1km建立网格；基于中国境内现有天然地震记录，确定中国境内各网格的板内地震震级得到中国境内板内地震震级三维立体分布；其中，在确定中国境内各网格的板内地震震级过程中，对于中国境内每一个网格，将同一网格内的地震记录按照能量进行叠加并基于叠加所得的能量确定该网格的板内地震震级；其中，当网格内叠加所得的能量为零时，将该网格叠加所得的能量定义为1；地震震级通过下述公式确定：lgE＝4.8+1.5M；式中，E为叠加所得的能量，J；M为地震里氏震级；
确定三维网格化的天然地震数据的强震区的截止值为4；
基于三维网格化的天然地震数据，利用三维网格化的天然地震数据的强震区的截止值，获取中国境内强震区的立体分布图(能够体现出强震区的立体分布特征)；结果如图5所示；
由图4、图5可以看出，在四川盆地西部频发板内地震，出现连续强震区。
构建中国境内的三维网格：东西向按照20km、南北向按照20km、纵向深度按照1km建立网格；基于目标区域内地震波层析成像数据，确定目标区域中各网格的地震波层析成像值，从而获得目标区域三维网格网格化的地震波层析成像数据；其中，在确定目标区域中各网格的地震波层析成像值过程中，对于目标区域中每一个网格，将同一网格内的地震波层析成像数据进行平均得到该网格的平均地震波层析成像值作为该网格的地震波层析成像值；优选地，当网格内没有地震波层析成像数据时，则取该网格上、下、左、右、前、后6个邻近的网格的地震波层析成像值计算几何平均值作为该网格的地震波层析成像值。
5、剖面网格化：
在中国境内内确定强震区的剖面A-A’的位置(如图6A所示)，取剖面两侧各20km范围内的区域作为该剖面的控制范围，在该剖面上按照长度为5km、纵向深度为1km建立网格；基于中国境内已发生板内天然地震记录，确定A-A’剖面控制范围内的板内天然地震记录，将A-A’剖面控制范围内的天然地震记录投点到A-A’剖面上(如图6B所示)，进而确定中国境内A-A’剖面中各网格的板内地震震级得到中国境内A-A’剖面板内地震震级剖面分布(结果如图6C所述)；其中，在确定中国境内A-A’剖面中各网格的板内地震震级过程中，对于A-A’剖面上的每一个网格，将同一网格内投点的地震记录按照能量进行叠加并基于叠加所得的能量确定该网格的板内地震震级；其中，当网格内叠加所得的能量为零时，将该网格叠加所得的能量定义为1；地震震级通过下述公式确定：lgE＝4.8+1.5M；式中，E为叠加所得的能量，J；M为地震里氏震级。
由图6C可知，龙门山地震带、华蓥山地震带和台湾地震带具有较强的地震能量异常地震多发生在30公里以下的深度。
6、板内地震成因的地质分析
(1)地震平面分布与地壳结构的关系
将地表高程图、不同分类震区的平面分布图与三维网格化的地震面波层析成像进行叠加，结合地壳GPS运动速率和震源机制解，分析强震区两侧地壳的运动情况、分析强震区与地壳地质结构的匹配关系，从而确定强震带的地质成因；其中，地表高程图、目标区域强震区的平面分布图与三维网格化的地震波层析成像进行叠加通过下述方式实现：
基于目标区域强震区的平面分布图确定目标区域中的强震区分布；基于目标区中的强震区分布，结合三维网格化的地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图；将目标区域强震区的平面分布图投影到地表高程图上，获取地表高程与强震区平面分布的叠合图；将所述地表高程与强震区平面分布的叠合图与所述地震波层析成像剖面图进行叠加；
具体而言，选择四川盆地西部的地表高程图、强震区分布图与层析成像剖面图相叠合(如图7A所示)，可以发现该区是扬子克拉通与青藏高原相互碰撞形成的；
(2)地震立体分布与地壳结构的关系
将目标区域强震区的立体分布图与三维网格化的地震面波层析成像进行叠加，分析强震区与地壳地质结构的匹配关系；其中，将目标区域强震区的立体分布图与三维网格化的地震波层析成像进行叠加通过下述方式实现：
基于目标区域强震区的立体分布图确定目标区域中的强震区分布；基于目标区中的强震区分布，结合三维网格化的地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图；将目标区域强震区的立体分布图与所述地震波层析成像剖面图进行叠加；
具体而言，选择四川盆地西部的强震区与层析成像剖面图相叠合(如图7B所示)，可以发现该区是扬子克拉通与青藏高原相互碰撞形成的。
本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (17)

1.一种基于多类型地球物理数据的地震评价方法，其中，该方法包括：
获取目标区域已发生板内天然地震数据、地表高程数据、布格重力异常数据、航磁异常数据、地震波层析成像、地壳GPS运动速率和震源机制解；其中，地震波层析成像包括地震面波层析成像和/或地震体波层析成像；
将目标区域平面网格化，构建得到目标区域的平面网格；根据目标区域的平面网格，将板内天然地震数据、地表高程数据、布格重力异常数据和航磁异常数据分别进行平面网格化处理，得到平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据；
基于平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据，确定目标区域强震区的平面分布图；
将地表高程图、目标区域强震区的平面分布图与地震波层析成像进行叠加，结合地壳GPS运动速率和震源机制解，分析强震区两侧地壳的运动情况、分析强震区与地壳地质结构的匹配关系，从而确定强震带的地质成因；
其中，根据目标区域的平面网格，将板内天然地震数据进行平面网格化处理包括：基于目标区域已发生板内天然地震数据，确定目标区域中各网格的板内地震震级得到目标区域内板内地震震级分布，从而获得目标区域平面网格化的天然地震数据；其中，在确定目标区域中各网格的板内地震震级过程中，对于目标区域中每一个网格，将同一网格内的板内天然地震按照能量进行叠加并基于叠加所得的能量确定该网格的板内地震震级；
其中，根据目标区域的平面网格，将地表高程数据进行平面网格化处理包括：基于目标区域内地表高程数据，确定目标区域中各网格的地表高程值，从而获得目标区域平面网格化的地表高程数据；其中，在确定目标区域中各网格的地表高程值过程中，对于目标区域中每一个网格，将同一网格内的地表高程数据进行平均得到该网格的平均地表高程值作为该网格的地表高程值；
其中，根据目标区域的平面网格，将布格重力异常数据进行平面网格化处理包括：基于目标区域内地表高程数据，确定目标区域中各网格的布格重力异常值，从而获得目标区域平面网格化的布格重力异常数据；其中，在确定目标区域中各网格的布格重力异常值过程中，对于目标区域中每一个网格，将同一网格内的布格重力异常数据进行平均得到该网格的平均布格重力异常值作为该网格的布格重力异常值；
其中，根据目标区域的平面网格，将航磁异常数据进行平面网格化处理包括：基于目标区域内地表高程数据，确定目标区域中各网格的航磁异常值，从而获得目标区域平面网格化的航磁异常数据；其中，在确定目标区域中各网格的航磁异常值过程中，对于目标区域中每一个网格，将同一网格内的航磁异常数据进行平均得到该网格的平均航磁异常值作为该网格的航磁异常值。
2.根据权利要求1所述的方法，其中，所述将目标区域平面网格化，构建得到目标区域的平面网格包括：东西向按照0.01°-0.1°，南北向按照0.01°-0.1°建立网格。
3.根据权利要求1所述的方法，其中，将板内天然地震数据进行平面网格化处理过程中，当网格内叠加所得的能量为零时，将该网格叠加所得的能量定义为1。
4.根据权利要求1所述的方法，其中，将地表高程数据进行平面网格化处理过程中，当网格内没有地表高程值数据时，则取该网格东、西、南、北的4个邻近的网格的地表高程值计算几何平均值作为该网格的地表高程值。
5.根据权利要求1所述的方法，其中，将布格重力异常数据进行平面网格化处理过程中，当网格内没有布格重力异常数据时，则取该网格东、西、南、北的4个邻近的网格的布格重力异常值计算几何平均值作为该网格的布格重力异常值。
6.根据权利要求1所述的方法，其中，将航磁异常数据进行平面网格化处理过程中，当网格内没有航磁异常数据时，则取该网格东、西、南、北的4个邻近的网格的航磁异常值计算几何平均值作为该网格的航磁异常值。
7.根据权利要求1所述的方法，其中，基于平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据，确定目标区域强震区的平面分布图包括：
对平面网格化的天然地震数据、平面网格化的地表高程数据、平面网格化的布格重力异常数据和平面网格化的航磁异常数据分别进行数据归一化处理，得到归一化的平面网格化的天然地震数据、归一化的平面网格化的地表高程数据、归一化的平面网格化的布格重力异常数据和归一化的平面网格化的航磁异常数据；
通过主成分分析方法对归一化的平面网格化的天然地震数据、归一化的平面网格化的地表高程数据、归一化的平面网格化的布格重力异常数据和归一化的平面网格化的航磁异常数据进行降维，得到第一主成分数据和第二主成分数据；
对各组主成分数据进行聚类分析，确定强震区的第一主成分门槛值和第二主成分门槛值，将目标区域分为强震区和非强震区，得到目标区域强震区的平面分布图；
其中，在对各组主成分数据进行聚类分析过程中，进一步包括：确定中震区和弱震区的第一主成分门槛值和第二主成分门槛值，将目标区域非强震区进一步分为中震区、弱震区和无震区，进而得到目标区域不同分类震区的平面分布图；
其中，第一主成分数据=0.91×归一化的平面网格化的天然地震数据-0.878×归一化地表高程+0.843×归一化的平面网格化的布格重力异常数据+0.835×归一化的平面网格化的航磁异常数据；
第二主成分数据=0.081×归一化的平面网格化的天然地震数据+0.0395×归一化地表高程-0.011×归一化的平面网格化的布格重力异常数据+0.196×归一化的平面网格化的航磁异常数据。
8.根据权利要求1所述的方法，其中，将地表高程图、目标区域强震区的平面分布图与地震波层析成像进行叠加包括：
基于目标区域强震区的平面分布图确定目标区域中的强震区分布；
基于目标区中的强震区分布，结合地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图；
将目标区域强震区的平面分布图投影到地表高程图上，获取地表高程与强震区平面分布的叠合图；
将所述地表高程与强震区平面分布的叠合图与所述地震波层析成像剖面图进行叠加。
9.根据权利要求1-8任一项所述的方法，其中，该方法还包括：根据目标区域强震区的平面分布图，确定经过目标区域强震区的剖面，并构建剖面网格；根据剖面网格，将板内天然地震数据进行剖面网格化处理，得到剖面网格化的天然地震数据；
其中，所述构建剖面网格将板内天然地震数据进行剖面网格化处理包括：取剖面两侧各1-50km范围内的区域作为该剖面的控制范围，在该剖面上按照长度为1-10km、纵向深度为0.5-5km建立网格；
其中，根据剖面网格，将板内天然地震数据进行剖面网格化处理包括：基于目标区域已发生天然地震记录，确定剖面控制范围内的天然地震记录，将剖面控制范围内的天然地震记录投点到剖面上，进而确定目标区域中各网格的板内地震震级；从而获得目标区域强震区剖面网格化的天然地震数据；其中，在确定目标区域中各网格的板内地震震级过程中，对于目标区域中每一个网格，将同一网格内的板内天然地震按照能量进行叠加并基于叠加所得的能量确定该网格的板内地震震级；其中，当网格内叠加所得的能量为零时，将该网格叠加所得的能量定义为1。
10.根据权利要求1-8任一项所述的方法，其中，该方法还包括：
将目标区域三维网格化，构建得到目标区域的三维网格；根据目标区域的三维网格，将板内天然地震数据进行三维网格化处理，得到三维网格化的天然地震数据；
确定三维网格化的天然地震数据的强震区的截止值；
基于三维网格化的天然地震数据，利用三维网格化的天然地震数据的强震区的截止值，获取目标区域强震区的立体分布图；
其中，所述目标区域三维网格化，构建得到目标区域的三维网格包括：东西向按照1-50km、南北向按照1-50 km、纵向深度按照0.1-10km建立网格；其中，根据目标区域的三维网格，将板内天然地震数据进行三维网格化处理包括：基于目标区域已发生板内天然地震数据，确定目标区域中各网格的板内地震震级得到目标区域内板内地震震级分布，从而获得目标区域三维网格化的天然地震数据；其中，在确定目标区域中各网格的板内地震震级过程中，对于目标区域中每一个网格，将同一网格内的板内天然地震按照能量进行叠加并基于叠加所得的能量确定该网格的板内地震震级；其中，当网格内叠加所得的能量为零时，将该网格叠加所得的能量定义为1；
其中，三维网格化的天然地震数据的强震区的截止值为4。
11.根据权利要求9所述的方法，其中，该方法还包括：
将目标区域三维网格化，构建得到目标区域的三维网格；根据目标区域的三维网格，将板内天然地震数据进行三维网格化处理，得到三维网格化的天然地震数据；
确定三维网格化的天然地震数据的强震区的截止值；
基于三维网格化的天然地震数据，利用三维网格化的天然地震数据的强震区的截止值，获取目标区域强震区的立体分布图；
其中，所述目标区域三维网格化，构建得到目标区域的三维网格包括：东西向按照1-50km、南北向按照1-50 km、纵向深度按照0.1-10km建立网格；其中，根据目标区域的三维网格，将板内天然地震数据进行三维网格化处理包括：基于目标区域已发生板内天然地震数据，确定目标区域中各网格的板内地震震级得到目标区域内板内地震震级分布，从而获得目标区域三维网格化的天然地震数据；其中，在确定目标区域中各网格的板内地震震级过程中，对于目标区域中每一个网格，将同一网格内的板内天然地震按照能量进行叠加并基于叠加所得的能量确定该网格的板内地震震级；其中，当网格内叠加所得的能量为零时，将该网格叠加所得的能量定义为1；
其中，三维网格化的天然地震数据的强震区的截止值为4。
12.根据权利要求10所述的方法，其中，该方法还包括：
将目标区域强震区的立体分布图与地震波层析成像进行叠加，分析强震区与地壳地质结构的匹配关系；
其中，将目标区域强震区的立体分布图与地震波层析成像进行叠加包括：
基于目标区域强震区的立体分布图确定目标区域中的强震区分布；
基于目标区中的强震区分布，结合地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图；
将目标区域强震区的立体分布图与所述地震波层析成像剖面图进行叠加。
13.根据权利要求11所述的方法，其中，该方法还包括：
将目标区域强震区的立体分布图与地震波层析成像进行叠加，分析强震区与地壳地质结构的匹配关系；
其中，将目标区域强震区的立体分布图与地震波层析成像进行叠加包括：
基于目标区域强震区的立体分布图确定目标区域中的强震区分布；
基于目标区中的强震区分布，结合地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图；
将目标区域强震区的立体分布图与所述地震波层析成像剖面图进行叠加。
14.根据权利要求8、12、13中任一项所述的方法，其中，基于目标区中的强震区分布，结合地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图包括：
将目标区域三维网格化，构建得到目标区域的三维网格；根据目标区域的三维网格，将地震波层析成像进行三维网格化处理，得到三维网格化的地震波层析成像；
基于目标区中的强震区分布，结合三维网格化的地震波层析成像，获取目标区域强震区及其两侧的地震波层析成像剖面图；
其中，所述目标区域三维网格化，构建得到目标区域的三维网格包括：东西向按照1-50km、南北向按照1-50 km、纵向深度按照0.1-10km建立网格；
其中，根据目标区域的三维网格，将地震波层析成像进行三维网格化处理包括：基于目标区域内地震波层析成像数据，确定目标区域中各网格的地震波层析成像值，从而获得目标区域三维网格网格化的地震波层析成像数据；其中，在确定目标区域中各网格的地震波层析成像值过程中，对于目标区域中每一个网格，将同一网格内的地震波层析成像数据进行平均得到该网格的平均地震波层析成像值作为该网格的地震波层析成像值；其中，当网格内没有地震波层析成像数据时，则取该网格上、下、左、右、前、后6个邻近的网格的地震波层析成像值计算几何平均值作为该网格的地震波层析成像值。
15.根据权利要求1或11所述的方法，其中，所述并基于叠加所得的能量确定该网格的板内地震震级通过下述公式实现：
Figure 288479DEST_PATH_IMAGE002
式中，E为叠加所得的能量，J；M为地震里氏震级。
16.根据权利要求9所述的方法，其中，所述并基于叠加所得的能量确定该网格的板内地震震级通过下述公式实现：
Figure DEST_PATH_IMAGE004A
式中，E为叠加所得的能量，J；M为地震里氏震级。
17.根据权利要求10所述的方法，其中，所述并基于叠加所得的能量确定该网格的板内地震震级通过下述公式实现：
Figure DEST_PATH_IMAGE004AA
式中，E为叠加所得的能量，J；M为地震里氏震级。
CN202110228947.8A 2021-03-02 2021-03-02 基于多类型地球物理数据的地震评价方法 Expired - Fee Related CN113009559B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202110228947.8A CN113009559B (zh) 2021-03-02 2021-03-02 基于多类型地球物理数据的地震评价方法
US17/679,222 US11360225B1 (en) 2021-03-02 2022-02-24 Earthquake evaluation method based on multi-type geophysical data

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110228947.8A CN113009559B (zh) 2021-03-02 2021-03-02 基于多类型地球物理数据的地震评价方法

Publications (2)

Publication Number Publication Date
CN113009559A CN113009559A (zh) 2021-06-22
CN113009559B true CN113009559B (zh) 2022-02-18

Family

ID=76402157

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110228947.8A Expired - Fee Related CN113009559B (zh) 2021-03-02 2021-03-02 基于多类型地球物理数据的地震评价方法

Country Status (2)

Country Link
US (1) US11360225B1 (zh)
CN (1) CN113009559B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113609729B (zh) * 2021-07-30 2023-09-01 中国科学院大学 一种东亚地区地幔动力学的数值模拟方法
CN113759425B (zh) * 2021-09-13 2022-04-01 中国科学院地质与地球物理研究所 井震联合评价深层古岩溶储层充填特征的方法与系统
CN114563820B (zh) * 2022-03-07 2023-03-24 中国矿业大学(北京) 地球物理监测方法、装置及系统
CN117572530B (zh) * 2024-01-17 2024-04-05 自然资源部第二海洋研究所 一种重力反演莫霍面和海底地震联合厘定洋陆边界的方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104360399A (zh) * 2014-12-10 2015-02-18 中国科学院地质与地球物理研究所 长导线接地源地下瞬变电磁探测方法及装置
CN105277976A (zh) * 2015-09-14 2016-01-27 山东科技大学 基于岩石露头雷达探测的地震正演模拟方法
CN110941021A (zh) * 2019-11-30 2020-03-31 西南交通大学 基于网格点格架函数的重力异常及梯度异常的正演方法

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
GB1109601A (en) * 1963-09-16 1968-04-10 Seismograph Service England Correlation of variables not previously known
US20090064602A1 (en) * 2007-09-11 2009-03-12 Il-Won Tech Co., Ltd Earthquake-resistant flat expansion joint using hinge
US8279707B2 (en) 2010-04-23 2012-10-02 Chevron U.S.A. Inc. Fusing geophysical data representing a geophysical space
CN102368094B (zh) * 2010-10-22 2013-10-16 中国石油天然气股份有限公司 一种地球物理勘探中的矩形窗带通保幅滤波数据处理方法
US20120271552A1 (en) * 2011-04-22 2012-10-25 Dashevsky Yuliy A Increasing the resolution of vsp ava analysis through using borehole gravity information
TWI636276B (zh) * 2014-05-16 2018-09-21 財團法人國家實驗研究院 人工智能地震判斷方法及地震偵測系統
CN107066481B (zh) 2016-12-21 2020-09-29 浙江大学 基于WebGIS的考古地球物理数据可视化平台的构建方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104360399A (zh) * 2014-12-10 2015-02-18 中国科学院地质与地球物理研究所 长导线接地源地下瞬变电磁探测方法及装置
CN105277976A (zh) * 2015-09-14 2016-01-27 山东科技大学 基于岩石露头雷达探测的地震正演模拟方法
CN110941021A (zh) * 2019-11-30 2020-03-31 西南交通大学 基于网格点格架函数的重力异常及梯度异常的正演方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
论上海市地热资源物探勘查;陆衍;《上海国土资源》;20151231;第36卷(第1期);67-72 *

Also Published As

Publication number Publication date
US11360225B1 (en) 2022-06-14
CN113009559A (zh) 2021-06-22

Similar Documents

Publication Publication Date Title
CN113009559B (zh) 基于多类型地球物理数据的地震评价方法
Gu et al. Shallow crustal structure of the Tanlu Fault Zone near Chao Lake in eastern China by direct surface wave tomography from local dense array ambient noise analysis
Lin et al. High-resolution 3D shallow crustal structure in Long Beach, California: Application of ambient noise tomography on a dense seismic array
CN104020492B (zh) 一种三维地震资料的保边滤波方法
CN110045432A (zh) 基于3d-glq的球坐标系下重力场正演方法及三维反演方法
Zhen et al. Segmentation of crustal velocity structure beneath the Shandong‐Jiangsu‐Anhui segment of the Tanlu Fault Zone and adjacent areas and its geological interpretations
Li et al. Crustal velocity structure of the Shaowu-Nanping-Pingtan transect through Fujian from deep seismic sounding-tectonic implications
Shomali et al. Crustal structure of Damavand volcano, Iran, from ambient noise and earthquake tomography
CN103926619A (zh) 一种三维vsp数据的逆时偏移方法
Irwansyah et al. Earthquake hazard zonation using peak ground acceleration (PGA) approach
Gang et al. Crust and upper mantle electrical resistivity structure in the Panxi region of the eastern Tibetan Plateau and its significance
Koulakov et al. Slab interactions in the Taiwan region based on the P-and S-velocity distributions in the upper mantle
Yang et al. Mantle flow in the vicinity of the eastern edge of the Pacific‐Yakutat slab: Constraints from shear wave splitting analyses
Wang Joint inversion of P-wave velocity and Vp-Vs ratio: imaging the deep structure in NE Japan
Seredkina Surface wave tomography of the Arctic from Rayleigh and Love wave group velocity dispersion data
Qin et al. Distribution of intra-crustal low velocity zones beneath Yunnan from seismic ambient noise tomography
Droznina et al. S-wave velocity model for several regions of the Kamchatka Peninsula from the cross correlations of ambient seismic noise
Aboud et al. Using a 3D gravity inversion technique to image the subsurface density structure in the Lunayyir volcanic field, Saudi Arabia
Xu et al. Insight into urban faults by wavelet multi-scale analysis and modeling of gravity data in Shenzhen, China
CN106646607A (zh) 一种提高ct反演分辨率和效率的自适应不等间距网格划分方法
Yong et al. Three-dimensional high-resolution velocity structure imaging and seismicity study of Yangbi Ms6. 4 earthquake
Levshin et al. Surface wave studies of the Greenland upper lithosphere using ambient seismic noise
Sun et al. Crustal velocity, density structure, and seismogenic environment in the southern segment of the North-South Seismic Belt, China
Mala et al. Earthquake Microzonation Study on Batubesi Dam of Nuha, East Luwu, South Sulawesi, Indonesia
Sobolev et al. Towards real‐time tsunami amplitude prediction

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20220218