CN113238280A - 一种基于格林函数的地震监测方法 - Google Patents

一种基于格林函数的地震监测方法 Download PDF

Info

Publication number
CN113238280A
CN113238280A CN202110704678.8A CN202110704678A CN113238280A CN 113238280 A CN113238280 A CN 113238280A CN 202110704678 A CN202110704678 A CN 202110704678A CN 113238280 A CN113238280 A CN 113238280A
Authority
CN
China
Prior art keywords
correlation
cross
earthquake
green
waveform
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.)
Granted
Application number
CN202110704678.8A
Other languages
English (en)
Other versions
CN113238280B (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.)
Chengdu Univeristy of Technology
Original Assignee
Chengdu Univeristy of Technology
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 Chengdu Univeristy of Technology filed Critical Chengdu Univeristy of Technology
Priority to CN202110704678.8A priority Critical patent/CN113238280B/zh
Publication of CN113238280A publication Critical patent/CN113238280A/zh
Application granted granted Critical
Publication of CN113238280B publication Critical patent/CN113238280B/zh
Active 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/28—Processing seismic data, e.g. analysis, for interpretation, for correction
    • G01V1/288—Event detection in seismic signals, e.g. microseismics

Landscapes

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

Abstract

本发明公开了一种基于格林函数的地震监测方法，包括：步骤一，对研究区域划分网格。准备格林函数，并以二进制文件的形式保存在计算机中。步骤二，格林函数与观测数据进行互相关，得到互相关格林函数，并以二进制文件的形式保存在计算机中。第三步，假定一个震源机制，依据震源机制对互相关格林函数进行加权叠加得到平均互相关波形。第四步，重复第二、三步得到所有的虚拟震源的平均互相关波形。第五步，依据全部的平均互相关波形的振幅值检测地震。本发明的优点是：使需要进行的互相关次数与震源机制的种数无关，从而减少互相关次数，具备更高的效率。

Description

一种基于格林函数的地震监测方法
技术领域
本发明涉及地震监测技术领域，特别涉及一种基于格林函数的地震识别、定位和震源机制反演方法。
背景技术
地震仪组成地震观测台网。通过地震观测台网，可实现对地振动的物理量，如加速度、速度和位移的测量和记录，这便是地震监测。地震监测是地震学研究和防震减灾的基础。相较于中强地震，小微地震的地震波具有能量弱，频率高，衰减快的特点。这导致小微地震是地震监测的难点。
波形匹配法是一种在既有观测硬件条件的基础上，加强监测能力的方法。波形匹配法使用已知地震的波形去和连续波形做互相关计算来寻找地震信号。这个已知地震称为模版地震。波形匹配法首先将模版地震的波形和连续观测记录做互相关。然后，将得到的互相关波形叠加、平均得到平均互相关波形。当平均互相关波形的振幅达到一定阈值的时候，即可识别到一个疑似地震。多个疑似地震可能对应的是同一个地震，所以需要假设一个时间窗。在这个时间窗内，具备最高互相关值的疑似地震记录会被程序确认，而其他记录则被移除。
上述的波形匹配法需要使用者提供模版地震。这就限制了波形匹配法的应用范围。比如，地震频发的地震带往往存在地震极少的地震空区。地震空区在未来存在发生强震的可能性，是涉及公众安全的重大课题，是地震学研究的热点。然而，地震空区往往就没有足够的地震可供用作模版。另外，前期缺乏研究的地区也有可能无法提供足够的模版地震。
因为上述的原因，近年产生了使用虚拟地震为模版的做法。虚拟地震的波形是通过人为假定震源参数，使用计算机计算出来的。利用这些计算出来的合成波形代替真实地震的波形去和连续观测记录进行互相关。这一做法可以解决没有模版地震的问题。但是，因为可能的震源参数很多，尤其是可能的震源机制解很多，这导致计算量巨大，使这一做法缺乏实践能力。
发明内容
本发明针对现有技术的缺陷，提供了一种基于格林函数的地震监测方法。
为了实现以上发明目的，本发明采取的技术方案如下：
一种基于格林函数的地震监测方法，包括以下步骤：
步骤一：地震监测区域划分为经度、纬度和深度的三维网格，网格所在位置为可能的震源位置。
步骤二：准备格林函数文件。依据一维速度模型、观测台站的位置坐标和步骤一的网格点的经纬度、深度坐标，准备好全部格林函数文件gi(i∈[1,3])。计算出格林函数的归一化系数。归一化系数是
Figure BDA0003131724460000021
2g1g2、2g2g3和2g1g3。
步骤三：格林函数与观测波形进行互相关。假定震源在某一个三维网格位置上，将实际的连续观测记录从ENZ坐标旋转到RTZ坐标。将格林函数与实际的连续观测做互相关，并得到半归一化的互相关格林函数：
Figure BDA0003131724460000022
u是某一分量的实际的观测波形；每一个分量都需要各自进行步骤三。
步骤四：对步骤三得到的半归一化的互相关格林函数进行加权叠加，并完成归一化。假定一个震源机制
Figure BDA0003131724460000023
(
Figure BDA0003131724460000024
是断层面的走向，δ是断层面的倾角，λ是断层面的滑动角)，对格林函数进行加权叠加，即可得平均互相关波形。首先，单分量的互相关波形可以由下面的公式(1)给出：
Figure BDA0003131724460000025
公式(1)中的
Figure BDA0003131724460000031
已经由步骤二完成。公式(1)中的ai是权重系数，可在表1中查得:
表1加权系数查询表
Figure BDA0003131724460000032
公式(1)中的v是理论波形，|v|由公式(2)得到：
Figure BDA0003131724460000033
在单分量的互相关波形得到后，叠加后再平均即可得平均互相关波形。将平均互相关波形中振幅达到了既定阈值的采样点的相关信息保存，平均互相关波形不保存。
步骤五：步骤四完成后在完成了一个震源位置上一个震源机制的疑似地震的探测。不断重复步骤四，完成同一个震源位置上不同震源机制的疑似地震的探测。
步骤六：重复步骤三至五，完成所有震源位置上所有震源机制的疑似地震的探测并保留相关信息。
步骤七：依据所有的疑似地震的记录，挑选出确定的探测到的地震的信息。
进一步地，步骤中三维网格覆盖需要进行地震监测的区域。三维网格的网格点的分布应该均匀覆盖。如果有需要特别监测的区域需要加密网格点。网格点的间距选择经纬度方向上2km，深度方向上5km。
进一步地，步骤七中的挑选方法是：首先选择出互相关值最大的结果作为确定的探测到的地震，然后找到互相关值第二大的结果，如果该结果的发震时刻和前面已经确定的地震的发震时刻不小于固定时窗，则也确定为探测到的地震，否则舍弃。不断重复这个过程，完成对疑似地震的挑选。
与现有技术相比，本发明的优点在于：
使需要进行的互相关次数与震源机制的种数无关，从而减少互相关次数，具备更高的效率。
附图说明
图1是本发明实施例地震监测方法流程图；
图2是本发明实施例研究区域的网格划分俯瞰图；
图3是本发明实施例一个观测数据由ENZ坐标旋转到RTZ坐标的对比图；
图4是本发明实施例一个波形记录与格林函数做互相关得到的互相关格林函数；
图5是本发明实施例平均互相关波形图。
具体实施方式
为使本发明的目的、技术方案及优点更加清楚明白，以下根据附图并列举实施例，对本发明做进一步详细说明。
如图1所示，一种基于格林函数的地震监测方法，包括以下步骤：
步骤一：网格划分。图2展示了该地区的大致情况。图中的三角形为既有的地震台站。蓝色小点为已观测到的地震。沙滩球表示一主震的震源机制解。红色圆点为划分的网格点(经纬度方向上是2km间隔，深度方向上5km间隔)。可以看到，这些网格应该均匀地覆盖了需要进行地震监测的区域。
步骤二：准备格林函数文件。依据上述的地震台站的经纬度坐标和网格点的经纬度坐标和深度，即可计算所有可能的震中距。再依据当地地区的一维速度模型，准备好可能的全部格林函数文件gi。最后，计算出所有格林函数的归一化系数。这些归一化系数是
Figure BDA0003131724460000051
2g1g2、2g2g3和2g1g3。这些归一化系数应当保存到相应的格林函数的文件的自定义头段变量上。
步骤三：格林函数与观测波形进行互相关。在上面的网格点中，选择任一一个作为假定的震源位置。依据假定的震中位置，将实际的连续观测波形记录从ENZ坐标旋转到RTZ坐标(图3)。Z分量不需要变换。对这些连续观测波形，找到与之对应的格林函数，并进行滑动互相关。在这里，互相关系数的具体表达式是：
Figure BDA0003131724460000052
(u是某一分量的实际的观测波形；gi是对应的格林函数)。所有连续观测波形和与之对应的所有格林函数都需要进行这样的操作，得到的互相关结果称为半归一化的互相关格林函数(图4)。
步骤四：对步骤三得到的半归一化的互相关格林函数进行加权叠加，并完成归一化。假定一个震源机制
Figure BDA0003131724460000053
对互相关格林函数进行加权叠加，即可得互相关波形，再得到平均互相关波形。下面说明，加权叠加和归一化的的具体方式。加权叠加和归一化应该同时进行。一个分量的互相关波形是
Figure BDA0003131724460000054
相关表达式的含义和计算方法请见发明内容部分。在单分量的互相关波形得到后，叠加后再平均即可得平均互相关波形(图5)。将平均互相关波形中振幅达到了既定阈值的采样点的相关信息保存在一个文本文件之中。一个地震的相关信息的震源位置和震源机制和之前假定的一致。互相关波形中时间即为检测到的地震的发震时刻。这里相关信息中的一条地震记录称为疑似地震。信息保留后，互相关格林函数、互相关波形和平均互相关波形全部删除。
步骤五：步骤四完成了在完成了一个震源位置(网格点)上一个震源机制的疑似地震的探测。不断重复步骤四，完成同一个震源位置上不同震源机制的疑似地震的探测。
步骤六：经历步骤三、四、五之后，完成了一个震源位置(网格点)上所有震源机制的疑似地震的探测。不断重复步骤三、四、五，完成所有震源位置上所有震源机制的疑似地震的探测并保留相关信息。
步骤七：依据前面的所有的疑似地震的记录，挑选出确定的探测到的地震的信息。挑选的方法是：首先选择出互相关值最大的结果作为确定的探测到的地震，然后找到互相关值第二大的结果，如果该结果的发震时刻和前面已经确定的地震的发震时刻不小于固定时窗，则也确定为探测到的地震，否则舍弃。不断重复这个过程，完成对疑似地震的挑选。
本发明的主要意义是相对现有技术，在保证效果不变的前提下，减少了互相关次数，增加了效率。下面分析本发明对互相关次数减少的情况：在上述的步骤三中知道，在本发明中互相关的具体公式为
Figure BDA0003131724460000061
和震源机制解无关，而只与格林函数gi有关。这样，本发明在有Nst个三分量地震仪时，互相关次数为8*Nst(因为参见表1，非零的加权系数是8个)。而传统方法是理论波形与实际波形做互相关。理论波形和震源机制相关，所以传统方法需要进行的互相关次数为Nfc*3*Nst。通常Nfc*3是远大于8的。因此，本发明相对于传统方法减少了互相关次数，从而提高了效率。
本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的实施方法，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (3)

1.一种基于格林函数的地震监测方法，其特征在于，包括以下步骤：
步骤一：地震监测区域划分为经度、纬度和深度的三维网格，网格所在位置为可能的震源位置；
步骤二：准备格林函数文件；依据一维速度模型、观测台站的位置坐标和步骤一的网格点的经纬度、深度坐标，准备好全部格林函数文件gi(i∈[1,3])；计算出格林函数的归一化系数；归一化系数是
Figure FDA0003131724450000011
2g1g2、2g2g3和2g1g3；
步骤三：格林函数与观测波形进行互相关；假定震源在某一个三维网格位置上，将实际的连续观测记录从ENZ坐标旋转到RTZ坐标；将格林函数与实际的连续观测做互相关，并得到半归一化的互相关格林函数：
Figure FDA0003131724450000012
u是某一分量的实际的观测波形；每一个分量都需要各自进行步骤三；
步骤四：对步骤三得到的半归一化的互相关格林函数进行加权叠加，并完成归一化；假定一个震源机制
Figure FDA0003131724450000013
(
Figure FDA0003131724450000014
是断层面的走向，δ是断层面的倾角，λ是断层面的滑动角)，对格林函数进行加权叠加，即可得平均互相关波形；首先，单分量的互相关波形可以由下面的公式(1)给出：
Figure FDA0003131724450000015
公式(1)中的
Figure FDA0003131724450000016
已经由步骤二完成；公式(1)中的ai是权重系数，可在表1中查得:
表1 加权系数查询表
Figure FDA0003131724450000017
Figure FDA0003131724450000021
公式(1)中的v是理论波形，|v|由公式(2)得到：
Figure FDA0003131724450000022
在单分量的互相关波形得到后，叠加后再平均即可得平均互相关波形；将平均互相关波形中振幅达到了既定阈值的采样点的相关信息保存，平均互相关波形不保存；
步骤五：步骤四完成后在完成了一个震源位置上一个震源机制的疑似地震的探测；不断重复步骤四，完成同一个震源位置上不同震源机制的疑似地震的探测；
步骤六：重复步骤三至五，完成所有震源位置上所有震源机制的疑似地震的探测并保留相关信息；
步骤七：依据所有的疑似地震的记录，挑选出确定的探测到的地震的信息。
2.根据权利要求1所述的一种基于格林函数的地震监测方法，其特征在于：步骤中三维网格覆盖需要进行地震监测的区域；三维网格的网格点的分布应该均匀覆盖；如果有需要特别监测的区域需要加密网格点；网格点的间距选择经纬度方向上2km，深度方向上5km。
3.根据权利要求1所述的一种基于格林函数的地震监测方法，其特征在于：步骤七中的挑选方法是：首先选择出互相关值最大的结果作为确定的探测到的地震，然后找到互相关值第二大的结果，如果该结果的发震时刻和前面已经确定的地震的发震时刻不小于固定时窗，则也确定为探测到的地震，否则舍弃；不断重复这个过程，完成对疑似地震的挑选。
CN202110704678.8A 2021-06-24 2021-06-24 一种基于格林函数的地震监测方法 Active CN113238280B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110704678.8A CN113238280B (zh) 2021-06-24 2021-06-24 一种基于格林函数的地震监测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110704678.8A CN113238280B (zh) 2021-06-24 2021-06-24 一种基于格林函数的地震监测方法

Publications (2)

Publication Number Publication Date
CN113238280A true CN113238280A (zh) 2021-08-10
CN113238280B CN113238280B (zh) 2023-02-24

Family

ID=77140685

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110704678.8A Active CN113238280B (zh) 2021-06-24 2021-06-24 一种基于格林函数的地震监测方法

Country Status (1)

Country Link
CN (1) CN113238280B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114563820A (zh) * 2022-03-07 2022-05-31 中国矿业大学(北京) 地球物理监测方法、装置及系统
CN115184995A (zh) * 2022-06-23 2022-10-14 中国地震局地震预测研究所 一种基于测震数据确定走滑型地震发震断层方向的方法
CN117687094A (zh) * 2023-11-14 2024-03-12 中国地震局地球物理研究所 一种估计情景地震高概率宽频带地震动的方法

Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007040949A (ja) * 2005-07-07 2007-02-15 Taisei Corp リアルタイム地震情報を利用した地震動の予測システム
US20100067328A1 (en) * 2008-09-17 2010-03-18 Andrew Curtis Interferometric directional balancing
US20110320180A1 (en) * 2010-06-29 2011-12-29 Al-Saleh Saleh M Migration Velocity Analysis of Seismic Data Using Common Image Cube and Green's Functions
KR20120019518A (ko) * 2010-08-26 2012-03-07 서울대학교산학협력단 송신원 추정을 통한 시간 영역 역시간 구조보정 방법 및 장치
CN104749627A (zh) * 2015-03-23 2015-07-01 中国石油集团川庆钻探工程有限公司地球物理勘探公司 基于相似性的微地震信号凸显方法
CN105929444A (zh) * 2016-04-08 2016-09-07 中国科学院地质与地球物理研究所 一种基于互相关偏移与最小二乘思想的微地震定位方法
CN106291684A (zh) * 2015-06-06 2017-01-04 中国石油化工股份有限公司 一种盲源地震波场的地震响应恢复与虚源道集构建方法
CN106646609A (zh) * 2017-01-13 2017-05-10 成都理工大学 多次扫描的微地震多参数联合快速反演方法
CN106908835A (zh) * 2017-03-01 2017-06-30 吉林大学 带限格林函数滤波多尺度全波形反演方法
CN107807385A (zh) * 2017-09-18 2018-03-16 成都理工大学 一种区域应力场时空变化反演方法
CN110058299A (zh) * 2018-09-14 2019-07-26 南方科技大学 地震定位方法、装置及终端设备
CN112068198A (zh) * 2020-08-24 2020-12-11 西南科技大学 基于地震波全波形特征的裂缝破裂尺度的描述方法

Patent Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2007040949A (ja) * 2005-07-07 2007-02-15 Taisei Corp リアルタイム地震情報を利用した地震動の予測システム
US20100067328A1 (en) * 2008-09-17 2010-03-18 Andrew Curtis Interferometric directional balancing
US20110320180A1 (en) * 2010-06-29 2011-12-29 Al-Saleh Saleh M Migration Velocity Analysis of Seismic Data Using Common Image Cube and Green's Functions
KR20120019518A (ko) * 2010-08-26 2012-03-07 서울대학교산학협력단 송신원 추정을 통한 시간 영역 역시간 구조보정 방법 및 장치
CN104749627A (zh) * 2015-03-23 2015-07-01 中国石油集团川庆钻探工程有限公司地球物理勘探公司 基于相似性的微地震信号凸显方法
CN106291684A (zh) * 2015-06-06 2017-01-04 中国石油化工股份有限公司 一种盲源地震波场的地震响应恢复与虚源道集构建方法
CN105929444A (zh) * 2016-04-08 2016-09-07 中国科学院地质与地球物理研究所 一种基于互相关偏移与最小二乘思想的微地震定位方法
CN106646609A (zh) * 2017-01-13 2017-05-10 成都理工大学 多次扫描的微地震多参数联合快速反演方法
CN106908835A (zh) * 2017-03-01 2017-06-30 吉林大学 带限格林函数滤波多尺度全波形反演方法
CN107807385A (zh) * 2017-09-18 2018-03-16 成都理工大学 一种区域应力场时空变化反演方法
CN110058299A (zh) * 2018-09-14 2019-07-26 南方科技大学 地震定位方法、装置及终端设备
CN112068198A (zh) * 2020-08-24 2020-12-11 西南科技大学 基于地震波全波形特征的裂缝破裂尺度的描述方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
S. DONNER 等: "Seismic moment tensors from synthetic rotational and translational ground motion: Green’s functions in 1-D versus 3-D", 《GEOPHYSICAL JOURNAL INTERNATIONAL》, 30 June 2020 (2020-06-30), pages 161 - 179 *
李世林等: "利用背景噪声互相关观测近对跖点台站间体波格林函数", 《地球物理学报》, vol. 59, no. 06, 15 June 2016 (2016-06-15), pages 2028 - 2038 *
王冲鹏 等: "地震虚源法边界效应研究", 《地球物理学进展》, vol. 36, no. 3, 31 March 2021 (2021-03-31), pages 1111 - 1121 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114563820A (zh) * 2022-03-07 2022-05-31 中国矿业大学(北京) 地球物理监测方法、装置及系统
CN115184995A (zh) * 2022-06-23 2022-10-14 中国地震局地震预测研究所 一种基于测震数据确定走滑型地震发震断层方向的方法
CN115184995B (zh) * 2022-06-23 2024-04-19 中国地震局地震预测研究所 一种基于测震数据确定走滑型地震发震断层方向的方法
CN117687094A (zh) * 2023-11-14 2024-03-12 中国地震局地球物理研究所 一种估计情景地震高概率宽频带地震动的方法

Also Published As

Publication number Publication date
CN113238280B (zh) 2023-02-24

Similar Documents

Publication Publication Date Title
CN113238280B (zh) 一种基于格林函数的地震监测方法
CN111164462B (zh) 一种人工源面波勘探方法、面波勘探装置及终端设备
Toomey et al. Tomographic imaging of the shallow crustal structure of the East Pacific Rise at 9° 30′ N
CN106093864A (zh) 一种麦克风阵列声源空间实时定位方法
JP5650353B2 (ja) 地震イベントパラメータの推定を取得する方法及びそのシステム、地震イベント探索エンジン
CN107831542B (zh) Ddw高精度深度域井震匹配方法
CN101893698B (zh) 噪声源测试分析方法及其装置
CN109001801B (zh) 基于多次迭代蚁群算法的断层变尺度识别方法
CN108957403B (zh) 一种基于广义互相关的高斯拟合包络时延估计方法及系统
CN105116448B (zh) 一种转换波方位各向异性校正方法及装置
CN111983676A (zh) 一种基于深度学习的地震监测方法及装置
CN110389377B (zh) 基于波形互相关系数相乘的微震偏移成像定位方法
CN103926623A (zh) 一种压制逆时偏移低频噪音的方法
CN105093319A (zh) 基于三维地震数据的地面微地震静校正方法
CN105629299A (zh) 角度域叠前深度偏移的走时、角度表获取方法及成像方法
CN104570116A (zh) 基于地质标志层的时差分析校正方法
CN105425299A (zh) 确定地层裂缝分布的方法和装置
CN112946743B (zh) 区分储层类型的方法
CN101825722B (zh) 一种鲁棒的地震信号瞬时频率的估计方法
CN106257309A (zh) 叠后地震数据体处理方法及装置
CN112114358A (zh) 一种基于三维地震资料表征的地下火山通道识别方法
CN111474580B (zh) 一种基于炮检距矢量片的方位角道集提取方法和系统
CN110873895A (zh) 一种变网格微地震逆时干涉定位方法
CN112630840B (zh) 基于统计特征参数的随机反演方法及处理器
CN112305591B (zh) 隧道超前地质预报方法、计算机可读存储介质

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