CN115932960A

CN115932960A - 全矩张量震源机制反演方法、系统、设备及存储介质

Info

Publication number: CN115932960A
Application number: CN202211527642.8A
Authority: CN
Inventors: 文健
Original assignee: University of Science and Technology of China USTC
Current assignee: University of Science and Technology of China USTC
Priority date: 2022-12-01
Filing date: 2022-12-01
Publication date: 2023-04-07

Abstract

本发明公开了一种全矩张量震源机制反演方法、系统、设备及存储介质，方法包括：步骤S1，获取控制参数，控制参数包括：速度模型、震源深度和台站位置；步骤S2，利用步骤S1获取的控制参数通过快速广义反透射系数算法计算不同震源深度的格林函数；步骤S3，将各格林函数代入以矩张量距离为停止准则的邻域算法中进行震源机制反演，得出最终震源深度和震源机制解。通过采用快速广义反透射系数算法计算格林函数，不仅提升了计算效率，也较低了存储量，且能计算任意台站‑震源分布的地震记录，通过采用以矩张量距离为停止准则的邻域算法进行震源机制反演，相比网格搜索提高了震源机制反演的效率。

Description

全矩张量震源机制反演方法、系统、设备及存储介质

技术领域

本发明涉及地震的震源机制反演领域，尤其涉及一种全矩张量震源机制反演方法、系统、装置及存储介质。

背景技术

震源机制反演是揭示地震断层或断裂活动的重要手段。反演得到的地震矩张量可以反映破裂的类型，估算破裂的尺度和能量释放的大小等。对于大部分地震断层，可以用简单位错模型解释(如图1所示)，即地震是发生在断层面上相对滑动。简单位错模型中有三个控制参数，其中断层走向φ和倾角δ控制断层面的方向，滑动角λ控制断层面两侧的相对滑动。但当地震是由火山喷发、水力压裂、塌陷或爆炸等引起时，简单位错模型就不适用了。此时需要引入两个无量纲的量β和γ来决定地震矩张量中非简单位错的比例。因此，准确得到地震矩张量就可以给出地震事件的类型等相关参数。因此，震源机制反演除了用于天然地震观测，在水力压裂油气开采、核爆监测等领域均有广泛应用。例如，如果是水力压裂地震，则可以估算出压裂的尺度，评估压裂效果；如果是爆炸，则可以判断爆炸类型和估算当量。

现有全波形震源机制反演算法有两种，一种是Dreger和Helmberger(1993)年提出的时间域震源机制反演算法(TDMT_INVC，Time-Domain Moment Tensor INVerse Code)，算法流程如图2所示。另一种是Zhu和Helmberger(1996)完善得到的GCAP方法(Cut AndPaste)(算法流程如图3所示)。二者对比发现，格林函数的计算思路相同，都是分别计算三种简单位错源和一种爆炸源在台站引起的三分量位移或速度记录并保存这12个量。在后续反演过程中，TDMT_INVC采用的是将非线性问题线性化后迭代反演的思路，而GCAP方法中是采用网格搜索直接搜索走向、倾角、滑动角等5个参数。

现有的利用波形进行震源机制反演的算法在生成格林函数时，一般采用Minson和Dreger(2008)年提出的分别计算四种特定断层引起的位移作为格林函数的方式计算。这四种特定断层分别为：垂直走滑断层，垂直倾滑断层，45°倾角的倾滑断层以及爆炸源。但这种计算方式要分别计算四种特定的断层引起的位移，这个计算过程有的算法可能需要两次，有的需要四次。例如，FK(Wand&Herrmann,1980；Zhu&Rivera,2002)需要计算两次，并存储四种震源三个分量的记录。

现有的格林函数计算方法至少具有以下缺点：

(1)需要进行两次或四次正演计算才能得到全矩张量反演需要的格林函数，计算效率较低。

(2)需要存储12个分量的记录，存储量大。

(3)计算速度慢，特别是当震源深度和台站埋深接近的时候。

(4)不能计算震源和台站位于同一深度或台站有一定埋深的情况，这类情形主要出现在对水力压裂地震进行观测、塌陷或爆炸引起的地震的观测中。

而在反演中使用的现有网格搜索算法存在反演速度和反演精度受网格大小控制，网格越小，反演速度越慢的缺点。

有鉴于此，特提出本发明。

发明内容

本发明的目的是提供了一种全矩张量震源机制反演方法、系统、设备及存储介质，能提高格林函数计算效率和减少存储，以及提高震源机制反演的效率，进而解决现有技术中存在的上述技术问题。

本发明的目的是通过以下技术方案实现的：

一种全矩张量震源机制反演方法，包括：

步骤S1，从各个地震数据中心或从已获得的原始波形数据中获取控制参数，控制参数包括：速度模型、震源位置和台站位置；

步骤S2，利用所述步骤S1获取的控制参数通过快速广义反透射系数算法计算得出不同震源深度的格林函数；

步骤S3，将计算得出的每个震源深度的格林函数代入以矩张量距离为停止准则的邻域算法中进行震源机制反演，得到不同震源深度的震源机制反演结果及误差函数，选取误差函数最小的震源深度和震源机制解作为该震源的最终震源深度和震源机制解，即完成震源机制的反演。

一种全矩张量震源机制反演系统，包括：

控制参数获取单元、格林函数计算单元和震源机制反演单元；其中，

所述控制参数获取单元，用于获取控制参数，控制参数包括：速度模型、震源位置和台站位置；

所述格林函数计算单元，与所述控制参数获取单元通信连接，用于利用所述控制参数获取单元获取的控制参数通过快速广义反透射系数算法计算得出不同震源深度的格林函数；

所述震源机制反演单元，与所述格林函数计算单元通信连接，用于将所述格林函数计算单元计算得出的每个震源深度的格林函数代入以矩张量距离为停止准则的邻域算法进行震源机制反演，得到不同震源深度的震源机制反演结果及误差函数，选取误差函数最小的震源深度和震源机制解作为该震源的最终震源深度和震源机制解。

一种处理设备，包括：

至少一个存储器，用于存储一个或多个程序；

至少一个处理器，能执行所述存储器所存储的一个或多个程序，在一个或多个程序被处理器执行时，使得所述处理器能实现本发明所述的方法。

一种可读存储介质，存储有计算机程序，当计算机程序被处理器执行时能实现本发明所述的方法。

与现有技术相比，本发明所提供的全矩张量震源机制反演方法、系统、设备及存储介质，其有益效果包括：

通过采用快速广义反透射系数算法计算格林函数，不仅提升了计算效率，也较低了存储量，且能计算任意台站-震源分布的地震记录，通过采用以矩张量距离为停止准则的邻域算法进行震源机制反演，相比网格搜索提高了震源机制反演的效率，且震源和台站可以分布任意位置，适用范围更广。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。

图1为表示大多数地震的简单位错模型示意图。

图2为现有时域震源参数反演方法的流程图。

图3为现有GCAP反演方法的流程图。

图4为本发明实施例提供的全矩张量震源机制反演方法的流程图。

图5为本发明实施例提供的全矩张量震源机制反演方法的具体处理流程图。

图6为本发明实施例提供的全矩张量震源机制反演方法的快速广义反透射系数算法计算不同震源深度格林函数的流程图。

图7为本发明实施例提供的全矩张量震源机制反演方法的以矩张量距离为停止准则的邻域算法的流程图。

图8为本发明实施例提供的快速GRTM算法(即格林函数计算方法)与现有GRTM算法计算结果对比示意图；其中，(a)为当台站和震源都位于地表时，原始GRTM(实线)和本发明快速GRTM(虚线)计算得到的波形对比图；(b)为原始的GRTM和本发明快速GRTM计算对于同一台站(GZ.BJT)同样震中距的震源不同震源深度的震源引起的位移所需要的时间对比，其中，三角为原始GRTM，圆圈为FGRTM。

图9为本发明实施例提供的快速GRTM算法与现有FK算法计算结果对比示意图；其中，(a)为FK和FGRTM计算对于同一台站(GZ.BJT)同样震中距的震源不同震源深度的震源(震源深度范围为0.1km到9km)引起的波形对比图，其中，黑实线为FK，黑虚线为FGRTM；(b)为两种算法计算时间对比，其中，三角为FK，圆圈为FGRTM。

图10为本发明实施例提供的全矩张量震源机制反演系统的结构框图。

具体实施方式

下面结合本发明的具体内容，对本发明实施例中的技术方案进行清楚、完整地描述；显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例，这并不构成对本发明的限制。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

首先对本文中可能使用的术语进行如下说明：

术语“和/或”是表示两者任一或两者同时均可实现，例如，X和/或Y表示既包括“X”或“Y”的情况也包括“X和Y”的三种情况。

术语“包括”、“包含”、“含有”、“具有”或其它类似语义的描述，应被解释为非排它性的包括。例如：包括某技术特征要素(如原料、组分、成分、载体、剂型、材料、尺寸、零件、部件、机构、装置、步骤、工序、方法、反应条件、加工条件、参数、算法、信号、数据、产品或制品等)，应被解释为不仅包括明确列出的某技术特征要素，还可以包括未明确列出的本领域公知的其它技术特征要素。

术语“由……组成”表示排除任何未明确列出的技术特征要素。若将该术语用于权利要求中，则该术语将使权利要求成为封闭式，使其不包含除明确列出的技术特征要素以外的技术特征要素，但与其相关的常规杂质除外。如果该术语只是出现在权利要求的某子句中，那么其仅限定在该子句中明确列出的要素，其他子句中所记载的要素并不被排除在整体权利要求之外。

除另有明确的规定或限定外，术语“安装”、“相连”、“连接”、“固定”等术语应做广义理解，例如：可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本文中的具体含义。

术语“中心”、“纵向”、“横向”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”、“顺时针”、“逆时针”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述和简化描述，而不是明示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本文的限制。

下面对本发明所提供的全矩张量震源机制反演方法、系统、装置及存储介质进行详细描述。本发明实施例中未作详细描述的内容属于本领域专业技术人员公知的现有技术。本发明实施例中未注明具体条件者，按照本领域常规条件或制造商建议的条件进行。本发明实施例中所用试剂或仪器未注明生产厂商者，均为可以通过市售购买获得的常规产品。

如图4所示，本发明实施例提供一种全矩张量震源机制反演方法，包括：

步骤S3，将得出的每个震源深度的格林函数代入以矩张量距离为停止准则的邻域算法中进行震源机制反演，得到不同震源深度的震源机制反演结果及误差函数，选取误差函数最小的震源深度和震源机制解作为该震源的最终震源深度和震源机制解。

参见图6，上述方法的步骤S2中，利用获取的控制参数通过快速广义反透射系数算法按以下方式计算得出不同震源深度的格林函数，包括：

步骤S21，对每一个震源台站对，计算震源到台站之间不同频率和波数的广义反透射系数；

步骤S22，计算不同阶数的Bessel函数在不同频率和波数时的Bessel函数值；

步骤S23，结合步骤S21和S22中计算得到的广义反透射系数和Bessel函数值，计算不同频率和波数时关于格林函数I_i，i＝1…10的被积函数的值；

步骤S24，利用谷峰平均法对步骤S23得到的格林函数I_i，i＝1…10值进行波数积分，然后再进行反傅立叶变换得到其时间域的值，即得出不同震源深度的格林函数。

参见图7，上述方法的步骤S3中，按以下方式将得出的每个震源深度的格林函数代入以矩张量距离为停止准则的邻域算法中进行震源机制反演，包括：

步骤S31，初始化，生成均匀随机分布在模型空间的n_s个初始模型并计算每个初始模型的误差函数，误差函数e为：

其中，

是某一震源台站的观测记录；

是根据地震在该震源台站上引起位移的公式计算得到的合成记录；cc是

和

的相关系数；所述模型空间由断层走向、倾角、滑动角等5个参数构成，每个初始模型均由随机的5个参数构成；

步骤S32，开始迭代求解，令loop＝1；

步骤S33，从前loop步生成的loop×n_s中选取误差函数最小的前n_r个模型；

步骤S34，在由n_r个模型参数构成的Voronoi体中重新生成n_s个新模型并计算各新模型的误差函数，令loop＝loop+1；

步骤S35，判断迭代次数是否到达预设的最大迭代次数或平均矩张量距离的小于预设值，若否，则返回步骤S33继续迭代；若是，则执行步骤S36；

步骤S36，输出误差函数最小的模型参数，即为误差函数最小的震源深度和震源机制解，作为该震源的最终震源深度和震源机制解，即完成震源机制的反演。

参见图10，本发明实施例还提供一种全矩张量震源机制反演系统，包括：

所述控制参数获取单元，用于从各个地震数据中心或从已获得的原始波形数据(即地震原始波形数据)中获取控制参数，控制参数包括：速度模型、震源位置(经纬度)和台站位置；

参见图6，上述系统的格林函数计算单元中，利用获取的控制参数通过快速广义反透射系数算法按以下方式计算得出不同震源深度的格林函数，包括：

步骤S21，对每一个震源台站对，计算震源到震源台站之间不同频率和波数的广义反透射系数；

参见图7，上述系统的震源机制反演单元中，按以下方式将得出的每个震源深度的格林函数代入以矩张量距离为停止准则的邻域算法中进行震源机制反演，包括：

其中，

是某一震源台站的观测记录；

和

的相关系数；

步骤S32，开始迭代求解，令loop＝1；

本发明实施例还提供一种处理设备，包括：

至少一个存储器，用于存储一个或多个程序；

至少一个处理器，能执行所述存储器所存储的一个或多个程序，在一个或多个程序被处理器执行时，使得所述处理器能实现上述的方法。

本发明实施例进一步提供一种可读存储介质，存储有计算机程序，当计算机程序被处理器执行时能实现上述的方法。

综上可见，本发明实施例的方法及系统，通过采用快速广义反透射系数算法计算格林函数，由于节省了修正反透射系数的计算步骤，不仅提升了计算效率，也较低了存储量，并且，引入了峰谷平均法计算波速积分，能够计算任意台站-震源分布的地震记录，即使台站和震源位于同一深度(对于这种情况，其它FK算法是无法计算的)也可以计算，通过采用以矩张量距离为停止准则的邻域算法进行震源机制反演，相比网格搜索提高了震源机制反演的效率。

为了更加清晰地展现出本发明所提供的技术方案及所产生的技术效果，下面以具体实施例对本发明实施例所提供的全矩张量震源机制反演方法及系统进行详细描述。

实施例1

如图4、5所示，本实施例提供一种全矩张量震源机制反演方法，包括：

步骤S1，从各个地震数据中心(或从已获得的原始波形数据(即地震原始波形数据))中获取控制参数，控制参数包括：速度模型、震源位置和震源台站位置；

步骤S3，将得出的每个震源深度的格林函数代入通过以矩张量距离为停止准则的邻域算法中进行震源机制反演，得到不同震源深度的震源机制反演结果及误差函数；

步骤S4，选取误差函数最小的震源深度和震源机制解作为该震源的最终震源深度和震源机制解，即完成震源机制反演。

参见图6，上述步骤S2中，快速广义反透射系数算法计算格林函数I_i的流程如图5所示，其计算原理依据以下公式(1)，根据广义反透射算法的原理(Chen,1999,Eq.77a～cand 78a～c)，对于一个地震，其中任一台站上引起的位移可以改写为：

其中，θ是台站方位角；格林函数I_i,i＝1…10是和Bessel函数和广义反透射系数有关的波数积分函数，公式如下。

其中，J_i和J_i′,i＝0,1,2为Bessel函数及其导数；ω是计算波形的频率；k是波数；r是震源中心到台站的距离；z是震源深度；

和

和第j层的广义反透射系数有关(具体公式参见Chen,1999,Eq.74a～c的内容)；M_pq,p,q＝x,y,z是归一化的地震矩张量，是断层的走向、倾角、滑动角以及β和γ的函数。

因此，基于上述公式(1)就可以进行震源机制反演，而I_i可当作一种等效的格林函数，其只与台站的震中距和震源深度有关，与台站方位角无关。此外，在原始的GRTM算法中，在计算广义反透射系数之前，有一个修正的反透射系数需要计算。而本发明使用的快速广义反透射系数算法将修正反透射系数的计算部分移除，直接计算广义反透射系数，这样本发明的快速广义反透射系数算法相比原始的GRTM计算效率提高了40％(参见图8b)；此外，由于引入了峰谷平均法计算波速积分，该快速广义反透射系数算法能够计算任意台站-震源分布的地震记录，即使台站和震源位于同一深度(对于这种情况，其它FK算法是无法计算的)也可以计算(图8a)。

因为震源机制，即地震矩张量是断层走向0≤φ＜360，倾角0≤δ≤90、滑动角-180＜λ≤180及震源球上的余纬度0≤β≤180和经度-30≤γ≤30的函数，所以在进行震源机制反演时，本发明的方法没有直接反演m_pq而是反演上述5个角度值。反演过程中会设定一个误差函数e：

其中，

是震源台站的观测记录；

是根据公式(1)计算得到的合成记录；cc是

和

的相关系数。

参见图7，上述步骤S3中，将得出的每个震源深度的格林函数代入通过以矩张量距离为停止准则的邻域算法中进行震源机制反演，具体流程为：

步骤S31，初始化，生成均匀随机分布在模型空间的n_s个初始模型并计算每个模型的误差函数；

步骤S32，开始迭代求解，令loop＝1；

步骤S33，从前loop步生成的loop×n_s中选取误差函数最小的n_r个模型；

步骤S34，在由n_r个模型参数构成的Voronoi体中重新生成n_s个新模型并计算相应的误差函数，令loop＝loop+1；

步骤S35，判断迭代次数是否到达预设的最大迭代次数或平均矩张量距离的小于预设值，若否则返回步骤S33继续迭代；若是，则执行步骤36；

本发明的以矩张量距离为停止准则的邻域算法中，前面迭代步骤的结果在最优解附近生成新的模型。相比于有非最优解的模型空间，该算法会更倾向于最优解附近的模型空间。该邻域算法有两个控制参数：每次迭代要生成的模型数目n_s和要重采样的Voronoi体的个数n_r。在该邻域算法中，假设n_s和n_r相等，此时邻域算法更接近一种随机采样算法而不是一种优化算法。此外，在本邻域算法中还引入了两个矩张量距离的概念。两个张量X＝(x_ij)和Y＝(y_ij)的距离定义为

其中，

‖X‖＝X·X。两个张量的距离的取值范围为[0,180]。

综上可见，本发明实施例的方法及系统，通过使用快速广义反透射系数算法计算格林函数，由于省去了修正反透射系数的计算步骤，不仅提升了计算效率，也降低了存储量，引入了峰谷平均法计算波速积分，不仅可以计算任意台站-震源分布的格林函数I_i。而且，在计算精度相同的情况下，还可以大大提高计算效率(图9(a)、图9(b))，此外，对每个台站，传统FK算法需要存储12个分量作为格林函数，而本发明的快速广义反透射系数算法只需要存储10个I_i，减少了大概17％的存储量；本发明步骤3的以矩张量距离为停止准则的邻域算法相比于传统的网格搜索算法，大大减少了误差函数计算的次数。而反演的效率和误差函数计算的次数是成正比的，误差函数计算次数越多反演速度越慢。假如网格搜索的网格大小为10°，那么网格搜索的精度即为10°。因此要完成一次反演，网格搜索需要进行1723680次误差函数计算。而对于该以矩张量距离为停止准则的领域算法，经测试，当n_s为1000时，反演效率和精度都能兼顾。假设最大迭代次数为50，那么邻域算法需要的误差函数计算次数为51000，因此，邻域算法至少比网格搜索要快33倍，大幅提升反演效率。

本领域普通技术人员可以理解：实现上述实施例方法中的全部或部分流程是可以通过程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random Access Memory，RAM)等。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。本文背景技术部分公开的信息仅仅旨在加深对本发明的总体背景技术的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人员所公知的现有技术。