一种微地震震源机制反演方法及装置的制造方法
【专利摘要】本发明提供了一种微地震震源机制反演方法及装置，涉及微地震技术领域。方法包括：获取微地震的原始波场记录；根据原始波场记录的特征选取多个基本矩张量；根据弹性波动方程有限差分正演模拟方法获取所述多个基本矩张量的模拟波场记录；根据波场叠加原理，建立所述模拟波场记录和所述原始波场记录的线性方程组，并根据最小二乘法求解所述线性方程组中各个基本矩张量对应的预设系数；根据求解后的预设系数对所述多个基本矩张量进行线性组合，生成反演后微地震震源的矩张量结果。本发明可以解决当前还没有一种适用性更强并且效率较高的方法用于研究微地震的震源机制的问题。
【专利说明】
_种微地震震源机制反演方法及装置
技术领域
[0001] 本发明涉及微地震技术领域，尤其涉及一种微地震震源机制反演方法及装置。
【背景技术】
[0002] 随着当前油田生产与开发的需要，用微地震方法对油田生产进行实时动态监测对 油田开发具有十分重要意义。其中，对微地震震源的研究是微地震监测中的一项主要内容， 目的是阐明微地震破裂的成因和震源机制，以便更好地认识和了解压裂裂缝。震源机制是 地震震源处地球介质的运动方式。震源机制研究的内容包括:确定地震断层面的方位和岩 体的错动方向，研究震源处岩体的破裂和运动特征，以及这些特征和震源所辐射的地震波 之间的关系。
[0003] 当前，常用的震源机制分析方法，可以通过震源机制反演或求解震源机制矩张量 解来实现。一般情况下，现有技术中对震源机制反演的研究主要分为两大类:一类是基于格 林函数反演矩张量的方法，该方法具有计算方便、效率高、对简单模型计算精度高等优点， 但对于抗噪性及非均质反演的误差较大;第二类是基于弹性波场的矩张量反演，该方法能 够适应多种类型的介质，但由于震源类型的研究没有格林函数反演方法充分，如果考虑震 源类型的变化，则会增加反演的变量，使反演更加复杂。因此，当前还没有一种适用性更强 并且效率较高的方法用于研究微地震的震源机制。
【发明内容】
[0004] 本发明的实施例提供一种微地震震源机制反演方法及装置，以解决当前还没有一 种适用性更强并且效率较高的方法用于研究微地震的震源机制的问题。
[0005] 为达到上述目的，本发明采用如下技术方案：
[0006] -种微地震震源机制反演方法，包括：
[0007] 获取微地震的原始波场记录；
[0008] 根据原始波场记录特征选取多个基本矩张量；
[0009] 根据弹性波动方程有限差分正演模拟方法获取所述多个基本矩张量的模拟波场 记录；
[0010] 根据波场叠加原理，建立所述模拟波场记录和所述原始波场记录的线性方程组， 并根据最小二乘法求解所述线性方程组中各个基本矩张量对应的预设系数；
[0011]根据求解后的预设系数对所述多个基本矩张量进行线性组合，生成反演后微地震 震源的矩张量结果。
[0012] 进一步的，该微地震震源机制反演方法，还包括：
[0013] 根据所述反演后微地震震源的矩张量结果进行正演模拟，获取得到第一波场记 录；
[0014] 将所述第一波场记录与所述原始波场记录进行对比，以确定第一波场记录与所述 原始波场记录的误差。
[0015]具体的，所述多个基本矩张量为6个基本矩张量，所述6个基本矩张量为：
[0018] 具体的，根据弹性波动方程有限差分正演模拟方法获取所述多个基本矩张量的模 拟波场记录，包括： -Mu 0 0 1 「.0 MS2 0"
[0019] 获取预先设置的ISO震源0 M:2 0 ，DC震源:M优=姐:1 0 0， _ 0 ：Q M}i\ [ 0 0 0 ~Mn 0 0: _ 以及CLVD震源：0 M:_ 0 ； 0 0 -2.V/ _ 4 J _
[0020] 根据ISO震源的加载公式、DC震源的加载公式和CLVD震源的加载公式确定6个基本 矩张量的模拟波场记录；
[0021] 其中，所述ISO震源的加载公式、DC震源的加载公式和CLVD震源的加载公式分别 为：
[0025] 其中，A]^+V2(? + 1/2,/，幻表示速度场Vx分量在n+1/2步时网格点（i+1/2，j，k)的变 化; AG l/::(/-l/2，./U)表示速度场Vx分量在n+1/2步时网格点（i-1/2，j，k)的变化量； Ar," 1/2 (/，./ +1/2, /〇表示速度场V y分量在n +1 / 2步时网格点（i，j +1 / 2，k )的变化量； AF," ^(/, /-l,/2,A〇表示速度场Vy分量在n+1/2步时网格点（i，j-1/2，k)的变化量； Ar" ^ (/，./，/「+1/2)表示速度场V z分量在n +1 / 2步时网格点（i，j，k +1 / 2 )的变化量； Af+!'2 (/，M -1/2)表示速度场V z分量在n + 1 / 2步时网格点（i，j，k -1 / 2 )的变化量； ^、"|/:(/ + 1/2,./_-1，幻表示速度场￥1分量在11+1/2步时网格点（1 + 1/2，」-1，1〇的变化量； AK"11/2(/ +1，./ -1/2, A')表示速度场Vy分量在n+1/2步时网格点（i+1，j-1/2，k)的变化量;dt为 微地震震源采样间隔，P为微地震受力点处的密度，V为微地震网格单元的体积，V = dxdydz; fii(dt ? n)为Mil对应的子波在dt ? n时亥lj的振幅值；f22(dt ? n)为M22对应的子波在dt ? n时 亥IJ的振幅值;f33(dt ? n)为M33对应的子波在dt ? n时刻的振幅值;f12(dt ? n)为M12对应的子 波在dt ? n时刻的振幅值;f2i(dt ? n)为M21对应的子波在dt ? n时刻的振幅值。
[0026] 具体的，根据波场叠加原理，建立所述模拟波场记录和所述原始波场记录的线性 方程组，并根据最小二乘法求解所述线性方程组中各个基本矩张量对应的预设系数，包括： [0027]根据波场叠加原理，建立所述模拟波场记录和所述原始波场记录的线性方程组：
[0030]根据最小二乘法求解所述线性方程组中各个基本矩张量对应的预设系数ki、k2、 k3、k4、ks、k6 〇
[0031 ] -种微地震震源机制反演装置，包括：
[0032] 微地震震源数据获取单元，用于获取微地震的原始波场记录；
[0033] 基本矩张量选取单元，用于根据原始波场记录的特征选取多个基本矩张量；
[0034] 模拟波场记录获取单元，用于根据弹性波动方程有限差分正演模拟方法获取所述 多个基本矩张量的模拟波场记录；
[0035] 预设系数确定单元，用于根据波场叠加原理，建立所述模拟波场记录和所述原始 波场记录的线性方程组，并根据最小二乘法求解所述线性方程组中各个基本矩张量对应的 预设系数；
[0036] 矩张量结果生成单元，用于根据求解后的预设系数对所述多个基本矩张量进行线 性组合，生成反演后微地震震源的矩张量结果。
[0037] 进一步的，所述的微地震震源机制反演装置，还包括：
[0038]正演模拟单元，用于根据所述反演后微地震震源的矩张量结果进行正演模拟，获 取得到第一波场记录；
[0039] 对比单元，用于将所述第一波场记录与所述原始波场记录进行对比，以确定第一 波场记录与所述原始波场记录的误差。
[0040] 具体的，所述基本矩张量选取单元得到的多个基本矩张量为6个基本矩张量，所述 6个基本矩张量为：
[0043] 此外，所述模拟波场记录获取单元，具体用于： _MU G 0 1 「0 M12 0"
[0044] 获取预先设置的ISO震源:0 M22 0 ，DC震源:M21 0 0， 〇 o mj3\ ooo 0 0 以及CLVD震源：#^加=0 M泣 0 ； 0 0 -2Mm
[0045] 根据ISO震源的加载公式、DC震源的加载公式和CLVD震源的加载公式确定6个基本 矩张量的模拟波场记录；
[0046] 其中，所述ISO震源的加载公式、DC震源的加载公式和CLVD震源的加载公式分别 为：
[0049]其中，Ar7+1/2(/ + l/2,y',幻表示速度场Vx分量在n+1/2步时网格点（i+l/2，j，k)的变 化量；Af+1/2(/ -1/2,/，幻表示速度场Vx分量在n+1 /2步时网格点（i-1 /2，j，k)的变化量； A f v 2 (/，./ +1/2，/c)表示速度场V y分量在n +1 / 2步时网格点（i，j +1 / 2，k )的变化量； 厶7,"|/2(/,./-1/2，幻表示速度场￥7分量在11 + 1/2步时网格点（1，」-1/2，1〇的变化量； M?' *(/,./,t_+_ 1/2)表示速度场Vz分量在n +1 /2步时网格点（i，j，k +1 /2 )的变化量； AfA' -1/2 >表示速度场V z分量在n + 1 / 2步时网格点（i，j，k -1 / 2 )的变化量； AK" 1/:(/ +1/2,- U)表示速度场Vx分量在n+1 /2步时网格点（i +1 /2，j-1，k)的变化量； /^::'|/:(/ + 1../-1/2，々)表示速度场￥7分量在11+1/2步时网格点（1+1，」-1/2，1〇的变化量;(^为 微地震震源采样间隔，P为微地震受力点处的密度，V为微地震网格单元的体积，V = dxdydz; fii(dt ? n)为Mil对应的子波在dt ? n时亥lj的振幅值；f22(dt ? n)为M22对应的子波在dt ? n时 亥IJ的振幅值;f33(dt ? n)为M33对应的子波在dt ? n时刻的振幅值;f12(dt ? n)为M12对应的子 波在dt ? n时刻的振幅值;f2i(dt ? n)为M21对应的子波在dt ? n时刻的振幅值。
[0050]此外，所述预设系数确定单元，具体用于：
[0051]根据波场叠加原理，建立所述模拟波场记录和所述原始波场记录的线性方程组：
~Rf ： k' R\ - 卩.丄 fc 焉 /为预设系数矩阵；；为原始波场记录矩阵； kA \ I VK_ : -災.
[0054]根据最小二乘法求解所述线性方程组中各个基本矩张量对应的预设系数ki、k2、 k3、k4、ks、k6 〇
[0055] 本发明实施例提供的微地震震源机制反演方法及装置，首先获取微地震的原始波 场记录，进而可以根据原始波场记录特征选取多个基本矩张量;之后，根据弹性波动方程有 限差分正演模拟方法获取所述多个基本矩张量的模拟波场记录;从而根据波场叠加原理， 建立所述模拟波场记录和所述原始波场记录的线性方程组，并根据最小二乘法求解所述线 性方程组中各个基本矩张量对应的预设系数;之后，根据求解后的预设系数对所述多个基 本矩张量进行线性组合，生成反演后微地震震源的矩张量结果。本发明能够适应微地震资 料的实时处理，提高水力压裂结果分析的准确性，可以解决当前没有一种适用性更强并且 效率较高的方法用于研究微地震的震源机制的问题。
【附图说明】
[0056] 为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现 有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本 发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可 以根据这些附图获得其他的附图。
[0057]图1为本发明实施例提供一种微地震震源机制反演方法的流程图一；
[0058]图2为本发明实施例提供一种微地震震源机制反演方法的流程图二；
[0059] 图3为本发明实施例中的基本力偶关系示意图；
[0060] 图4为本发明实施例中的两种不同震源强度的记录曲线示意图；
[0061] 图5为本发明实施例中震源强度与记录振幅之间的线性关系曲线示意图；
[0062] 图6为本发明实施例中ISO震源、DC震源和两者混合的震源的记录曲线示意图； [0063]图7为本发明实施例中ISO和DC震源记录相加与混合震源的差异示意图；
[0064] 图8为本发明实施例中经过矩张量反演后的记录与模拟记录差异示意图；
[0065] 图9为本发明实施例中的反演记录和含噪记录的波形示意图；
[0066]图10为本发明实施例中的一种三分量弹性波正演模拟波形示意图；
[0067]图11为本发明实施例中的6个基本矩张量的波形记录示意图；
[0068]图12为本发明实施例中的通过反演得到的矩张量进行弹性波正演模拟的波形记 录示意图和与原始模拟记录的差异示意图；
[0069]图13为本发明实施例中反演前后单道波形对比示意图；
[0070] 图14为本发明实施例提供的一种微地震震源机制反演装置的结构示意图一；
[0071] 图15为本发明实施例提供的一种微地震震源机制反演装置的结构示意图二。
【具体实施方式】
[0072] 下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完 整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于 本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他 实施例，都属于本发明保护的范围。
[0073] 如图1所示，本发明实施例提供一种微地震震源机制反演方法，包括：
[0074] 步骤101、获取微地震的原始波场记录。
[0075] 步骤102、根据原始波场记录的特征选取多个基本矩张量。
[0076]步骤103、根据弹性波动方程有限差分正演模拟方法获取所述多个基本矩张量的 模拟波场记录。
[0077]步骤104、根据波场叠加原理，建立所述模拟波场记录和所述原始波场记录的线性 方程组，并根据最小二乘法求解所述线性方程组中各个基本矩张量对应的预设系数。
[0078] 步骤105、根据求解后的预设系数对所述多个基本矩张量进行线性组合，生成反演 后微地震震源的矩张量结果。
[0079] 本发明实施例提供的微地震震源机制反演方法，首先获取微地震的原始波场记 录，进而根据原始波场记录的特征选取多个基本矩张量;之后，根据弹性波动方程有限差分 正演模拟方法获取所述多个基本矩张量的模拟波场记录;从而根据波场叠加原理，建立所 述模拟波场记录和所述原始波场记录的线性方程组，并根据最小二乘法求解所述线性方程 组中各个基本矩张量对应的预设系数;之后，根据求解后的预设系数对所述多个基本矩张 量进行线性组合，生成反演后微地震震源的矩张量结果。本发明能够适应微地震资料的实 时处理，提高水力压裂结果分析的准确性，可以解决当前没有一种适用性更强并且效率较 高的方法用于研究微地震的震源机制的问题。
[0080] 为了使本领域的技术人员更好的了解本发明，下面列举一个更为详细的实施例， 值得说明的是为了适应微地震资料的实时处理，提高水力压裂结果分析的准确性，进而对 开采方案进行及时调整，本发明是基于波动方程正演理论和弹性波场叠加原理的微地震震 源机制矩张量反演的方法，通过本发明能够对微地震的破裂过程和裂缝空间几何形态进行 定量描述，为压裂工程及储层改造提供基础数据。
[0081] 如图2所示，本发明实施例提供一种微地震震源机制反演方法，包括：
[0082]步骤201、获取微地震的原始波场记录。
[0083]步骤202、根据原始波场记录的特征选取多个基本矩张量。
[0084] 此处，该多个基本矩张量可以为6个基本矩张量，所述6个基本矩张量为： "1 0 〇1 「0 10"! 「0 0 I-
[0085] Mn - 0 0： 0 = 0： Q 0： Ml3 = 0 0 0 0 0 0」 0 0」 L〇 0 0 '()：0： 0] 「0 0 0] 「0 f) () _
[0086] m2% = 0 1 0 M a = 0 0 1 m33 = 0: 0： 0 , 0 0 (:)」 [0 0 0」 [_0 0 1
[0087] 步骤203、根据弹性波动方程有限差分正演模拟方法获取所述多个基本矩张量的 模拟波场记录。
[0088]下面就弹性波动方程有限差分正演模拟方法进行详细阐述：
[0089] -般情况下，在均匀各向异性介质中的弹性波方程可写为：
[0091]式中Ui和Uk为位移场;Cijki是刚度张量，又称为弹性矩阵;P为密度;fi为体力；t是 传播时间。
[0092]在常规地震勘探的正演模拟中，震源函数h是不包含角度、偏振等信息的。但是在 水力压裂裂缝过程中，产生的微地震震源通常是比常规地震复杂的，包含丰富的震源类型 以及裂缝的走滑角，倾向和倾角等角度特征。地震矩张量（Seismic Moment Tensor，简称 SMT)是描述这一复杂震源的有利工具。本发明将地震矩张量理论与弹性波方程结合来描述 三种基本的震源类型：ISO震源类型、DC震源类型和CLVD震源类型。本发明是基于Graves在 1996年提出的震源加载方法，在震源的强度及矩张量的表达方式上有所改变。在震源加载 分类中，主要有两种最为普遍:将震源加载在应力场和速度场。本发明中所用的方法是将震 源加载在速度场中。
[0093 ]如图3所示，一般情况下所有震源都可以由9对基本的力偶来表示，其表达式为： ~Mn M12 Mr：
[0094] M = Mn MZ2 M2i Mn Mm M33_
[0095]其中M为矩张量，Mij( i = 1，2，3; j = 1，2，3)为每对基本力偶。本发明以位于X分量的 力偶fX为例。其包含三个分量:Mn，M12，M13，分别位于X，Y，Z轴(如图3所示）。矩张量中的力矩 为dx，所以所受力要除以dx。因此，力矩Mn的表达式为：
[0097]其中，dt为采样间隔，P为该受力点处的密度，V为网格单元的体积，V = dxdydz，f ⑴为无量纲的单位震源。类似的，M12和M13可分别表示为：
[0100] 分量fy和fZ有类似的表达式，此处不再赘述。
[0101] 通过对这9对基本的力偶进行组合，可以产生全部的基本震源IS0、DC和CLVD类型。 本发明中只描述最具代表性的预先设置的三类震源。该预先设置的ISO震源表达式为： "Mu 0 0 1 「0 Mn O' MIS0 = 0 Mn 0 ，DC震源的表达式为:= _M2, 0 〇，以及CLVD震源的表达 0 0 从泊 0 0 0 ~Mn 0 0 式为0 Mv_ 0 。 0 0 -2M3i
[0102] 根据ISO震源的加载公式、DC震源的加载公式和CLVD震源的加载公式确定6个基本 矩张量的模拟波场记录。
[0103] 其中，所述ISO震源的加载公式、DC震源的加载公式和CLVD震源的加载公式分别 为：
[0107] 其中，AF广1/2(/ +1/2,义幻表示速度场Vx分量在n+1/2步时网格点（i+1/2，j，k)的变 化量；AF^1/2(/-l/2,y,幻表示速度场Vx分量在n+1/2步时网格点（i-1/2，j，k)的变化量； AF,n+1/2 (W +1/2，: &)表示速度场V y分量在n +1 / 2步时网格点（i，j +1 / 2，k)的变化量； A G 1/2 (/,./ -1/2, /〇表示速度场V y分量在n +1 / 2步时网格点（i，j -1 / 2，k )的变化量； 1/2 (/，》+1/2)表示速度场V z分量在n +1 / 2步时网格点（i，j，k +1 / 2 )的变化量； Af心(/，./，/、? -1/2)表示速度场V z分量在n + 1 / 2步时网格点（i，j，k - 1 / 2 )的变化量； Af "(Z + l/X y-U)表示速度场Vx分量在n+1/2步时网格点（i + 1/2, j-l，k)的变化量； AC/:(/ + L./-l/2，/〇表示速度场Vy分量在n+1/2步时网格点（i+1，j-1/2，k)的变化量;dt为 微地震震源采样间隔，P为微地震受力点处的密度，V为微地震网格单元的体积，V = dxdydz; fii(dt ? n)为Mil对应的子波在dt ? n时亥lj的振幅值；f22(dt ? n)为M22对应的子波在dt ? n时 亥IJ的振幅值;f33(dt ? n)为M33对应的子波在dt ? n时刻的振幅值;f12(dt ? n)为M12对应的子 波在dt ? n时刻的振幅值;f2i(dt ? n)为M21对应的子波在dt ? n时刻的振幅值。当计算第(n+ 1 /2)步的速度场时，所施加的外力等价于速度场额外的变化量。随后在计算第n+1步时的应 力场时，施加在速度场的变换量转化为应力场的变换量，从而产生正应力和切应力。根据以 上速度场震源进行弹性波数值模拟，能够得到三种基本震源在各向异性介质中的波场记 录。
[0108] 此处，弹性波方程有限差分正演模拟的主要作用是:矩张量可以用来描述微地震 事件的各种震源机制，能够以具体的数值张量表示裂缝产生时的应力状态;使用矩张量工 具加载到弹性波方程中，以所施加的矩张量作为波场的初始值，进行弹性波场的模拟;最后 记录模拟的波场，观测不同震源机制在不同介质中的响应，为震源机制解释提供理论支持。
[0109] 步骤204、根据波场叠加原理，建立所述模拟波场记录和所述原始波场记录的线性 方程组，并根据最小二乘法求解所述线性方程组中各个基本矩张量对应的预设系数。
[0110] 此处，对波场叠加原理进行详细阐述：
[0111]①震源强度线性叠加原理：
[0112] -般情况下，震源强度越强，得到的地震记录越强。当震源强度在弹性破裂范围 内，且研究模型较小的情况下，可以认为震源强度的大小与记录振幅的大小呈线性关系。即 多个震源的叠加引起的记录与多个震源记录的和是相同的，即波场叠加原理。当波在介质 中传播时，如果存在几列波，则每列波能够保持自己的传播规律而不响应其它波。在波的重 叠范围内，各点的振幅值为各列波在该点引起的振幅的矢量和。此原理成立的条件为介质 为线性，且震动不强烈。而微地震监测中，产生的震动微弱，震源与检波器间的间距一般不 超过1000m，因此微地震事件的记录可以认为符合波场叠加原理。
[0113]震源振幅强度与记录强度的关系可以表示为：
[0114] f(k ? M)=k ? f(M)
[0115] 其中M为矩张量，k为系数。即震源矩张量为kM产生的波场等于k倍的震源矩张量M 产生的波场。如图4所示，模拟了震源强度为100和1000的两个记录，然后将震源强度为100 的记录乘以10并与震源强度为1000的记录做差，发现两者误差很小，波形匹配一致。图5为 震源强度与记录振幅之间的线性关系曲线，也证明了地震波在传播过程中是符合波场的叠 加原理的。
[0116] ②震源类型线性叠加原理
[0117] 根据上述震源强度与记录振幅的关系，也可以推断出震源类型与记录振幅的线性 关系，即不同种类的震源在传播过程中仍然符合波场的叠加原理。其关系可表示为：
[0118] f (M1+M2+M3) =f (Mi)+f (M2)+f (M3)
[0119] 其中施^3为不同的矩张量，代表不同类型的震源。
[0120] 例如，IS0、DC和两者混合的震源的记录如图6所示。ISO震源与DC和ISO的混合震源 记录的差异较小，只在波峰和波谷处可以看出差异。DC震源记录比较平缓。将ISO和DC震源 的记录相加，比较其与混合震源的差异，如图7所示。两种曲线匹配吻合较好。因此，可以认 为公式f (M4M2+M3) = f (M〇+f (M2)+f (M3)在微地震监测中是成立的，即可以将矩张量先分解 再进行波场模拟。
[0121]此外，在步骤204中所建立的模拟波场记录和原始波场记录的线性方程组，可以表 示为： ;可以 筒记为Ax = b。
[0123]其中，n为微地震监测中所用三分量检波器的个数；
为模拟波场记录矩阵， 「之1 : t Ri - pi 可以简记为A; 3为预设系数矩阵，可以简记为X;丨为原始波场记录矩阵，可以简即为 k4 t R,; k-5 ： kJ ，丨 RH y _ b。之后，可以根据最小二乘法求解所述线性方程组中各个基本矩张量对应的预设系数h、 k2、k3、k4、ks、k6 〇
[0124] 步骤205、根据求解后的预设系数对所述多个基本矩张量进行线性组合，生成反演 后微地震震源的矩张量结果。
[0125] 步骤206、根据所述反演后微地震震源的矩张量结果进行正演模拟，获取得到第一 波场记录。
[0126] 步骤207、将所述第一波场记录与所述原始波场记录进行对比，以确定第一波场记 录与所述原始波场记录的误差。
[0127] 对于上述步骤205至步骤207，其具体实现和实验过程及结果可以如下：
[0128] 发明人在研究过程中首先将页岩层简化为均匀介质进行初步的试验论证。在均匀 介质中，受到的外力是均匀传播的，因此需要考虑ISO震源。如果页岩层可近似认为水平层 状介质，则剪切应力产生的震源会发生在水平层，即会产生水平面的DC震源。为方便计算， 在均匀介质中的模拟，只考虑了这两种震源进行模拟。
[0129] 图8中实线为ISO和DC混合震源的单道记录。应用最小二乘反演方法后，得到单道 记录为反演虚线。从图8中可以发现，经过矩张量反演后的记录与模拟记录差异很小，匹配 的较好。原始的矩张量为Mo,反演得到的矩张量为M。结果为： '1 1 0] [0.9977 0.9824 0
[0130] lf0= 1 1 0 , M= 0.9824 0 9977 0 _0 0 IJ [ 0 0 0.9977_
[0131] 如图9所示，反演得到结果(反演记录)与原始记录(含噪记录)的波形得到匹配，求 得的矩张量值误差也很小。说明在该种简化的模型中，本发明所提出的反演方法是成立的。
[0132] 上述反演方法得到的是假设的ISO震源和DC震源的组合矩张量，其算法的准确性 在于假设条件的符合程度。相比上述反演方法，本发明的算法含有更多的未知系数，计算量 有所增加。
[0133] 本发明假设已知震源的矩张量为： ^3 1 2
[0134] ,M0= 1. 3 3 J 3之
[0135] 其三分量弹性波正演模拟波形如图10所示，波形更为复杂。由于三类基本震源的 组合成分是未知的，因此不能应用前文反演算法。通过求取每一个矩张量分量的波场，通过 最小二乘反演方法，可求得各个分量的系数，最后将各个矩张量分量组合求得原始矩张量。
[0136] 运用弹性波方程矩张量模拟方法，首先得到6个基本矩张量的波形记录，如图11所 示。通过得到的模拟记录，可建立上述公式Ax = b等式左边的系数矩阵A。等式右边矩阵b由 原始矩张量正演模拟记录得到，如图1 〇所示。
[0137] 对于公式Ax = b中未知系数x的求取，本发明中所用的反演算法为最小二乘法。通 过对全波形记录的匹配，计算得到的矩张量为： "2.999980 0.999424 2.000005"
[0138] M= 0.999424 2.999970 3.000471 2.000005 3.000471 3.000003
[0139] 其结果与原始矩张量差异很小，无噪情况下，相对误差控制在0.1 %以内，基本满 足了精度要求。
[0140] 根据反演得到的矩张量进行弹性波正演模拟，其波形记录如图12的a所示，与原始 模拟记录的差异如图12的b所示。其波形记录振幅的误差小于0.015%，相对于原始波形振 幅值，可以忽略不计。因此，可以证明反演得到的波形与原始波形差异很小，得到的震源机 制矩张量是十分准确的。同样，通过对反演前后单道波形对比及其残差分析，如图13所示， 也能证明反演结果的精确性。
[0141] 本发明实施例提供的微地震震源机制反演方法，首先获取微地震的原始波场记 录，进而可以根据原始波场记录的特征选取多个基本矩张量;之后，根据弹性波动方程有限 差分正演模拟方法获取所述多个基本矩张量的模拟波场记录;从而根据波场叠加原理，建 立所述模拟波场记录和所述原始波场记录的线性方程组，并根据最小二乘法求解所述线性 方程组中各个基本矩张量对应的预设系数;之后，根据求解后的预设系数对所述多个基本 矩张量进行线性组合，生成反演后微地震震源的矩张量结果。本发明能够适应微地震资料 的实时处理，提高水力压裂结果分析的准确性，可以解决当前没有一种适用性更强并且效 率较高的方法用于研究微地震的震源机制的问题。
[0142] 如图14所示，本发明实施例提供一种微地震震源机制反演装置，包括：
[0143] 微地震震源数据获取单元31，用于获取微地震的原始波场记录。
[0144] 基本矩张量选取单元32,根据原始波场记录的特征选取多个基本矩张量。
[0145] 模拟波场记录获取单元33，用于根据弹性波动方程有限差分正演模拟方法获取所 述多个基本矩张量的模拟波场记录。
[0146] 预设系数确定单元34,用于根据波场叠加原理，建立所述模拟波场记录和所述原 始波场记录的线性方程组，并根据最小二乘法求解所述线性方程组中各个基本矩张量对应 的预设系数。
[0147] 矩张量结果生成单元35,用于根据求解后的预设系数对所述多个基本矩张量进行 线性组合，生成反演后微地震震源的矩张量结果。
[0148] 进一步的，如图15所示，该微地震震源机制反演装置，还包括：
[0149] 正演模拟单元36,用于根据所述反演后微地震震源的矩张量结果进行正演模拟， 获取得到第一波场记录。
[0150] 对比单元37,用于将所述第一波场记录与所述原始波场记录进行对比，以确定第 一波场记录与所述原始波场记录的误差。
[0151] 具体的，所述基本矩张量选取单元32得到的多个基本矩张量为6个基本矩张量，所 述6个基本矩张量为： '1 0 〇1 「0 10] 「0 0 1_
[0152] Mn = 0 0 Q M12= 〇 0 0 Ma = 0 0 0 :0 0 0」 :0 0」 0 0 0: '〇 o 〇1 「〇: .〇: .〇~| 「〇 o
[0153] Af22 - 0 1 0 =00 1 A# 33 = 0 0 0 〇： 0 0 0」 [0 0 0」 [0 0 1
[0154] 此外，所述模拟波场记录获取单元33，具体用于： ~MU 0 0 1 r 0 Mn 0"
[0155] 获取预先设置的130震源：从/15〇= 0 M1Z 0 ，DC震源：财21 0 0, _ 0 0 M,3J L 0 0 0 ^Mu 0 0 - 以及CLVD震源：0 Mri 0 ： 0 0 -2M" _ Si _
[0156] 根据ISO震源的加载公式、DC震源的加载公式和CLVD震源的加载公式确定6个基本 矩张量的模拟波场记录。
[0157] 其中，所述ISO震源的加载公式、DC震源的加载公式和CLVD震源的加载公式分别 为：
[0161]其中，Afn+1/2(f+l/2,y,/:)表示速度场Vx分量在n+1/2步时网格点（i+l/2，j，k)的变 化量；Af 1/::(/ -1/2，./，々）表示速度场Vx分量在n+1 /2步时网格点（i-1 /2，j，k)的变化量； A f v 2 (/，./ +1/2，/c)表示速度场V y分量在n +1 / 2步时网格点（i，j +1 / 2，k )的变化量； 厶7,"|/2(/,./-1/2，幻表示速度场￥7分量在11 + 1/2步时网格点（1，」-1/2，1〇的变化量； M?' *(/,./,t_+_ 1/2)表示速度场Vz分量在n +1 / 2步时网格点（i，j，k +1 /2 )的变化量； Af 1/:(/../，A -1/2)表示速度场Vz分量在n+1/2步时网格点（i，j，k-1/2)的变化量； AK" 1/:(/ + l/2,./-U)表示速度场Vx分量在n+1/2步时网格点（i + 1/2, j-l，k)的变化量； K /_ j/2,/、?)表示速度场Vy分量在n+1/2步时网格点（ i+1 j-1//2，k)的变化量;扣为 微地震震源采样间隔，P为微地震受力点处的密度，V为微地震网格单元的体积，V = dxdydz; fii(dt ? n)为Mil对应的子波在dt ? n时亥lj的振幅值；f22(dt ? n)为M22对应的子波在dt ? n时 亥IJ的振幅值;f33(dt ? n)为M33对应的子波在dt ? n时刻的振幅值;f12(dt ? n)为M12对应的子 波在dt ? n时刻的振幅值;f2i(dt ? n)为M21对应的子波在dt ? n时刻的振幅值。
[0162] 此外，所述预设系数确定单元34,具体用于：
[0163] 根据波场叠加原理，建立所述模拟波场记录和所述原始波场记录的线性方程组：
~Rf ： k' R\ - 卩.丄 fc 焉 /为预设系数矩阵；；为原始波场记录矩阵。 kA \ I VK_ : -災.
[0166] 根据最小二乘法求解所述线性方程组中各个基本矩张量对应的预设系数ki、k2、 k3、k4、ks、k6 〇
[0167] 值得说明的是，本发明实施例提供的微地震震源机制反演装置的具体实现方式可 以参见上述图1和图2对应的方法实施例，此处不再赘述。
[0168] 本发明实施例提供的微地震震源机制反演装置，首先获取微地震的原始波场记 录，进而根据原始波场记录的特征选取多个基本矩张量;之后，根据弹性波动方程有限差分 正演模拟方法获取所述多个基本矩张量的模拟波场记录;从而根据波场叠加原理，建立所 述模拟波场记录和所述原始波场记录的线性方程组，并根据最小二乘法求解所述线性方程 组中各个基本矩张量对应的预设系数;之后，根据求解后的预设系数对所述多个基本矩张 量进行线性组合，生成反演后微地震震源的矩张量结果。本发明能够适应微地震资料的实 时处理，提高水力压裂结果分析的准确性，可以解决当前没有一种适用性更强并且效率较 高的方法用于研究微地震的震源机制的问题。
[0169] 本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序 产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实 施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机 可用存储介质（包括但不限于磁盘存储器、CD-ROM、光学存储器等）上实施的计算机程序产 品的形式。
[0170] 本发明是参照根据本发明实施例的方法、设备（系统）、和计算机程序产品的流程 图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流 程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序 指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产 生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实 现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0171] 这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特 定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指 令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或 多个方框中指定的功能。
[0172] 这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计 算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或 其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一 个方框或多个方框中指定的功能的步骤。
[0173]本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例 的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员， 依据本发明的思想，在【具体实施方式】及应用范围上均会有改变之处，综上所述，本说明书内 容不应理解为对本发明的限制。
【主权项】
1. 一种微地震震源机制反演方法，其特征在于，包括： 获取微地震的原始波场记录； 根据原始波场记录的特征选取多个基本矩张量； 根据弹性波动方程有限差分正演模拟方法获取所述多个基本矩张量的模拟波场记录； 根据波场叠加原理，建立所述模拟波场记录和所述原始波场记录的线性方程组，并根 据最小二乘法求解所述线性方程组中各个基本矩张量对应的预设系数； 根据求解后的预设系数对所述多个基本矩张量进行线性组合，生成反演后微地震震源 的矩张量结果。2. 根据权利要求1所述的微地震震源机制反演方法，其特征在于，还包括： 根据所述反演后微地震震源的矩张量结果进行正演模拟，获取得到第一波场记录； 将所述第一波场记录与所述原始波场记录进行对比，以确定第一波场记录与所述原始 波场记录的误差。3. 根据权利要求2所述的微地震震源机制反演方法，其特征在于，所述多个基本矩张量 为6个基本矩张量，所述6个基本矩张量为：4. 根据权利要求3所述的微地震震源机制反演方法，其特征在于，根据弹性波动方程有 限差分正演模拟方法获取所述多个基本矩张量的模拟波场记录，包括：获取预先远 以及CLVD震源： 根据ISO震源的加载公式、DC震源的加载公式和CLVD震源的加载公式确定6个基本矩张 量的模拟波场记录； 其中，所述ISO震源的加载公式、DC震源的加载公式和CLVD震源的加载公式分别为：其中，4(_|/：!(/ + 1/2,./,幻表示速度场￥1分量在11+1/2步时网格点（1+1/2，」，1〇的变化量； Δ fl/2(/ -1/2，./, λ·)表示速度场V x分量在η +1 / 2步时网格点（i -1 / 2，j，k )的变化量； ΔC +1/2，/〇表示速度场V y分量在η +1 / 2步时网格点（i，j +1 / 2，k )的变化量； 厶7,"|/2(/,./-1/2，幻表示速度场￥7分量在11 + 1/2步时网格点（1，」-1/2，1〇的变化量； ΔΓ?/,./, 1/2)表示速度场V z分量在η +1 / 2步时网格点（i，j，k +1 / 2 )的变化量； (Λ /, A· -1/2)表示速度场Vz分量在η + 1 / 2步时网格点（i，j，k- 1 / 2 )的变化量； ΔΓ |/2(/_+_1/2，./-1,/〇表示速度场￥1分量在11+1/2步时网格点（1 + 1/2，」-1，1〇的变化量； 么广"':(/+1,./-1/2，々)表示速度场￥7分量在11+1/2步时网格点（1+1，」-1/2，1〇的变化量;(^为 微地震震源采样间隔，P为微地震受力点处的密度，V为微地震网格单元的体积，V = dxdydz; fii(dt · η)为Mii对应的子波在dt · η时亥Ij的振幅值；f22(dt · η)为M22对应的子波在dt · η时 亥IJ的振幅值;f33(dt · η)为M33对应的子波在dt · η时刻的振幅值;f12(dt · η)为M12对应的子 波在dt · η时刻的振幅值;f2i(dt · η)为M2i对应的子波在dt · η时刻的振幅值。5.根据权利要求4所述的微地震震源机制反演方法，其特征在于，根据波场叠加原理， 建立所述模拟波场记录和所述原始波场记录的线性方程组，并根据最小二乘法求解所述线 性方程组中各个基本矩张量对应的预设系数，包括： 根据波场叠加原理，建立所述模拟波场记录和所述原始波场记录的线性方程组：为预设系数矩阵I原始波场记录矩阵； 根据最小二乘法求解所述线性方程组中各个基本矩张量对应的预设系数khk^khlu、 ks、1? 〇6. -种微地震震源机制反演装置，其特征在于，包括： 微地震震源数据获取单元，获取微地震的原始波场记录； 基本矩张量选取单元，用于根据原始波场记录的特征选取多个基本矩张量； 模拟波场记录获取单元，用于根据弹性波动方程有限差分正演模拟方法获取所述多个 基本矩张量的模拟波场记录； 预设系数确定单元，用于根据波场叠加原理，建立所述模拟波场记录和所述原始波场 记录的线性方程组，并根据最小二乘法求解所述线性方程组中各个基本矩张量对应的预设 系数； 矩张量结果生成单元，用于根据求解后的预设系数对所述多个基本矩张量进行线性组 合，生成反演后微地震震源的矩张量结果。7. 根据权利要求6所述的微地震震源机制反演装置，其特征在于，还包括： 正演模拟单元，用于根据所述反演后微地震震源的矩张量结果进行正演模拟，获取得 到第一波场记录； 对比单元，用于将所述第一波场记录与所述原始波场记录进行对比，以确定第一波场 记录与所述原始波场记录的误差。8. 根据权利要求7所述的微地震震源机制反演装置，其特征在于，所述基本矩张量选取 单元得到的多个基本矩张量为6个基本矩张量，所述6个基本矩张量为：9. 根据权利要求8所述的微地震震源机制反演装置，其特征在于，所述模拟波场记录获 取单元，具体用于： 获取预先? 以及CLVD震源：根据ISO震源的加载公式、DC震源的加载公式和CLVD震源的加载公式确定6个基本矩张 量的模拟波场记录； 其中，所述ISO震源的加载公式、DC震源的加载公式和CLVD震源的加载公式分别为：其中，Δ(/ + 1/2../U)表示速度场Vx分量在n+1 /2步时网格点（i+1/2，j，k)的变化量； E示速度场Vx分量在n + 1/2步时网格点（i-1/2, j，k)的变化量； :示速度场Vy分量在n + 1/2步时网格点（i，j + l/2，k)的变化量； 示速度场V y分量在η +1 / 2步时网格点（i，j -1 / 2，k )的变化量； 示速度场Vz分量在n + 1/2步时网格点（i，j，k+l/2)的变化量； 亓速度场Vz分量在n + 1/2步时网格点（i，j，k-l/2)的变化量； >表示速度场Vx分量在n+1 /2步时网格点（i +1 /2，j -1，k)的变化量； 良示速度场Vy分量在n+1/2步时网格点（1+1，」-1/2，1〇的变化量;(^为 微地震震源采样间隔，P为微地震受力点处的密度，V为微地震网格单元的体积，V = dxdydz; fii(dt · η)为Mii对应的子波在dt · η时亥Ij的振幅值；f22(dt · η)为M22对应的子波在dt · η时 亥IJ的振幅值;f33(dt · η)为M33对应的子波在dt · η时刻的振幅值;f12(dt · η)为M12对应的子 波在dt · η时刻的振幅值;f2i(dt · η)为M2i对应的子波在dt · η时刻的振幅值。10.根据权利要求9所述的微地震震源机制反演装置，其特征在于，所述预设系数确定 单元，具体用于： 根据波场叠加原理，建立所述模拟波场记录和所述原始波场记录的线性方程组：其中，η为微地震监测中所用三分量检波器的个数；根据最小二乘法求解所述线性方程组中各个基本矩张量对应的预设系数khk^khlu、 ks、k_6 〇
【文档编号】G01V1/28GK105891884SQ201610490599
【公开日】2016年8月24日
【申请日】2016年6月28日
【发明人】王润秋, 李会俭, 陈筱青, 江勇勇, 黄炜霖, 周延鑫, 袁民, 袁一民, 郝学兵
【申请人】中国石油大学(北京)