一种加速微震数值模拟的方法及装置
阅读：1057发布：2020-07-03
专利汇可以提供一种加速微震数值模拟的方法及装置专利检索，专利查询，专利分析的服务。并且本 发明 涉及一种微震数值模拟方法及装置，属于微震监测领域，具体是涉及一种 加速 微震数值模拟的方法及装置。本发明结合微震监测中 震源 与接收 检波器 数量特点，依据源检互换原理，在 现有技术 基础 上，将原震源点作为接收点、原接收点作为震源点，计算 地震 波 场，能够快速得到微震数值模拟中各个接收点的波场信息，与现有技术相比，大大提高了计算效率。，下面是一种加速微震数值模拟的方法及装置专利的具体信息内容。
1.一种加速微震数值模拟的方法，其特征在于，包括：
正演参数给定步骤，根据待模拟对象的地质条件，建立网格速度模型，并指定检波器坐标；
时间场计算步骤，根据震源点个数Ns和检波点个数Nr，确定地震波场计算方式；若Ns＞Nr，进行源检互换及波场计算，即将每个检波点当做震源点，分别计算每个检波点所在点作为震源的地震波场；
波场抽取步骤，抽取每个检波点的地震波场，正演结束。
2.根据权利要求1所述的一种加速微震数值模拟的方法，其特征在于，所述时间场计算步骤具体包括：
地震波场确定子步骤，计算当前震源点在已知速度模型下地震波场的程函方程；
程函方程计算子步骤，采用上风差分法对程函方程进行离散；
地震走时计算子步骤，利用目标点周围八个方向节点处的走时值计算其地震走时的二阶差分计算公式得到目标点走时的多值解；
目标走时确定子步骤，按照窄带扩展基本原理，在求出的目标点走时的多值解中选取最小走时点，即为目标点走时。
3.根据权利要求2所述的一种加速微震数值模拟的方法，其特征在于，所述地震波场确定子步骤中，基于下式确定地震波场的程函方程：
其中t是地震波走时，s是慢度，即速度的倒数，x、z代表二维直角坐标系的2个坐标轴。
4.根据权利要求2所述的一种加速微震数值模拟的方法，其特征在于，所述地震走时计算子步骤中，采用上风差分法对上式进行离散，具体为：
表示走时t在点(i,j)处x方向上的向前、向后的差分算子，
表示走时t在点(i,j)处z方向上的向前、向后的差分算子；
其中：
式中Δx、Δz分别为X、Z方向的网格尺寸。
5.根据权利要求2所述的一种加速微震数值模拟的方法，其特征在于，所述地震走时计算子步骤中，基于下式利用目标点周围八个节点处的走时值计算目标点走时的多值解：
式中，t是地震波走时，ti,j是在点(i,j)处的地震波走时，s是地震波的慢度，h＝Δx＝Δz为正方形网格尺寸，
利用上式计算出多个走时，按照窄带扩展基本原理，在求出的多值解中选取最小走时点，即为目标点走时，计算公式如下：
ti,j＝min(t1i,j,t2i,j,t3i,j,t4i,j,t5i,j,t6i,j)。
一种加速微震数值模拟的方法及装置
技术领域
[0001] 本发明涉及一种微震数值模拟方法及装置，属于微震监测领域，具体是涉及一种加速微震数值模拟的方法及装置。
背景技术
[0002] 微震监测技术是近几年在石油、煤田领域应用较广的一项物探技术，被广泛应用于压裂监测、冲击地压监测等领域。微震监测技术主要涉及正演和反演两大技术问题。微震监测正演问题，即已知速度模型、观测系统、震源坐标等参数，求解微震波场(动力学)及地震波走时(运动学)等信息，为微震监测观测系统设计、定位方法正确性验证等提供依据，是微震监测中的基础问题，也是反演的基础。微震监测反演问题，是正演的逆过程，即已知观测系统、地震波场和走时等信息，求取速度模型、震源参数等信息的过程。
[0003] 微震正演方法主要包括：基于波动方程理论的有限差分法、有限元方法、伪谱法、波前快速推进法等，基于积分方程的边界积分法、边界元法等，基于射线追踪理论的试射法、弯曲法等等。
[0004] 常规微震监测正演数值模拟算法的基本流程是：1.给定模拟区域的三维网格速度模型、震源及检波点坐标等参数；2.计算每个震源点在三维模型空间的地震波场(运动学或者动力学)；3.在地震波场中提取每个检波点处的信息。由于微震监测正演是一个三维问题，模型空间较大，计算三维模型空间的地震波场是一个非常耗时的过程，其主要影响因素包括模型空间大小、网格划分尺寸等。若模型空间较大、网格划分较细，计算一个震源点的三维地震波场需要几十分钟甚至几百分钟时间。假设有M个微震震源，计算1个震源需要N分钟时间，完成所有震源点正演模拟的时间大致是M×N分钟。因此，上述第2步是制约微震正演模拟效率的主要因素。
发明内容
[0005] 本发明主要是解决现有技术所存在的微震数值模拟(正演)速度较慢的技术问题，提供了一种加速微震数值模拟的方法及装置。该方法及装置结合微震监测中震源与检波点数量特点，依据源检互换原理，能够快速得到微震数值模拟中各个检波点的波场信息，与现有技术相比，大大提高了计算效率。
[0006] 本发明的上述技术问题主要是通过下述技术方案得以解决的：
[0007] 一种加速微震数值模拟的方法，具体包括以下步骤：
[0009] 波场计算步骤，根据震源点个数Ns和检波点个数Nr，确定地震波场计算方式。若Ns＞Nr，进行源检互换及波场计算，即将每个检波点当做震源点，分别计算每个检波点所在点作为震源的地震波场；若Ns＜Nr，则按照常规方法，计算每个震源点地震波场。微震监测中，一般属于前者。
[0010] 波场抽取步骤，基于上述计算得到的三维地震波场，抽取每个震源点的地震波场，正演结束。
[0011] 优选的，所述正演模拟参数给定步骤执行以下子步骤：
[0012] 根据待模拟区域建立XYZ三维直角坐标系；基于所述三维直角坐标系对待模拟区域进行网格划分；基于所述网格划分结果，对网格点进行速度赋值，建立速度模型；指定每个检波点的空间位置。其中，网格速度模型、检波点空间坐标及微震震源点分别记为：
[0013] 网格速度模型：vi,j,k(1≤i≤Nx,1≤j≤Ny,1≤k≤Nz)
[0014] 检波点坐标:(xri,yri,zrri)(1≤i≤Nr)
[0015] 震源点坐标：(xsi,ysi,zsri)(1≤i≤Ns)
[0016] 变量说明：vi,j,k为网格节点(i,j,k)对应的速度，Nx为模型X方向网格数，Ny为模型Y方向网格数，Nz为模型Z方向网格数，xri为第i个检波点的X坐标，yri为第i个检波点的Y坐标，zri为第i个检波点的Z坐标，Nr为检波点个数，xsi为第i个震源点的X坐标，ysi为第i个震源点的Y坐标，zsi为第i个震源点的Z坐标，Ns为震源点个数。
[0017] 优选的，所述波场计算模块执行以下子步骤：
[0018] 首先，根据震源点个数Ns和检波点个数Nr，确定地震波场计算方式。若Ns＜Nr，按常规遍历震源点，若Ns＞Nr，依次遍历各个检波点，将检波点当做震源点，原震源点标记为检波点。
[0019] 其次，计算当前震源点，在已知速度模型下的地震波场。以快速推进算法为例，求解下式二维直角坐标系(x-z)下的程函方程：
[0020]
[0021] 其中t是地震波走时，t(x,z)表示在坐标点(x,z)的地震波走时，s是慢度，即速度的倒数，s2(x,z)表示在坐标点(x,z)的慢度值的平方，x、z代表二维直角坐标系的2个坐标轴。
[0023]
[0024] 其中D为差分算子，i、j分别为x方向和z方向的网格节点索引号，max为最大值函数，min为最小值函数。上式可简化为：
[0025]
[0026] 表示走时t在点(i,j)处x方向上的向前、向后的差分算子，表示走时t在点(i,j)处z方向上的向前、向后的差分算子。
[0027] 为保证计算精度，采用二阶上风差分算子。即：
[0028]
[0029]
[0030] 式中Δx、Δz分别为X、Z方向的网格尺寸,ti,j表示在点(i,j)处的走时值。
[0031] 最终，利用目标点(i,j)周围八个节点(i,j±1)，(i,j±2)，(i±1,j)，(i±2,j)处的走时值，共有6种组合方式计算(i,j)点的走时值，6种组合方式的二阶差分计算公式如下：
[0032]
[0033] 式中，t是地震波走时，s是地震波的慢度，h＝Δx＝Δz为正方形网格尺寸。
[0034] 利用上式计算出多个走时，按照窄带扩展基本原理，在求出的多值解中选取最小走时点，即为目标点走时，计算公式如下：
[0035] ti,j＝min(t1i,j,t2i,j,t3i,j,t4i,j,t5i,j,t6i,j)
[0036] 遍历完所有节点，即可获得地震波时间场。
[0037] 优选的，所述波场抽取步骤执行以下子步骤：
[0038] 抽取每个检波点的地震波场，即在上步中求得的时间场t(x,z)中找到检波点位置对应的时间，就是原震源点到检波点的地震走时，正演结束。
[0039] 一种加速微震数值模拟的装置，具体包括以下步骤：
[0040] 正演模拟参数给定装置，根据待模拟对象的地质条件，建立网格速度模型，并指定检波点坐标。
[0041] 波场计算装置，根据震源点个数Ns和检波点个数Nr，确定地震波场计算方式。若Ns＞Nr，进行源检互换及波场计算，即将每个检波点当做震源点，分别计算每个检波点所在点作为震源的地震波场；若Ns＜Nr，则按照常规方法，计算每个震源点地震波场。微震监测中，一般属于前者。
[0042] 波场抽取装置，基于上述计算得到的三维地震波场，抽取每个检波点的地震波场，正演结束。
[0043] 优选的，所述正演模拟参数给定模块执行以下子步骤：
[0044] 根据待模拟区域建立XYZ三维直角坐标系；基于所述三维直角坐标系对待模拟区域进行网格划分，基于所述网格划分结果，对网格点进行速度赋值，建立速度模型；并指定各个检波点的空间位置。其中，网格速度模型、检波点空间坐标及潜在震源点分别记为：
[0045] 网格速度模型：vi,j,k(1≤i≤Nx,1≤j≤Ny,1≤k≤Nz)
[0046] 检波点坐标:(xri,yri,zrri)(1≤i≤Nr)
[0047] 震源点坐标：(xsi,ysi,zsri)(1≤i≤Ns)
[0048] 变量说明：vi,j,k为网格节点(i,j,k)对应的速度，Nx为模型X方向网格数，Ny为模型Y方向网格数，Nz为模型Z方向网格数，xri为第i个检波点的X坐标，yri为第i个检波点的Y坐标，zri为第i个检波点的Z坐标，Nr为检波点个数，xsi为第i个震源点的X坐标，ysi为第i个震源点的Y坐标，zsi为第i个震源点的Z坐标，Ns为震源点个数。
[0049] 优选的，所述波场计算模块执行以下子步骤：
[0050] 首先，根据震源点个数Ns和检波点个数Nr，确定地震波场计算方式。若Ns＞Nr，按常规遍历震源点，若Ns＜Nr，依次遍历各个检波点，将检波点当做震源点，原震源点标记为检波点。
[0051] 其次，计算当前震源点，在已知速度模型下的地震波场。以快速推进算法为例，求解下式二维直角坐标系(x-z)下的程函方程：
[0052]
[0053] 其中t是地震波走时，s是慢度，即速度的倒数，x、z代表二维直角坐标系的2个坐标轴。
[0054] 采用上风差分格式代替式程函方程中的偏微分格式，即采用上风差分法对上式进行离散，其离散式可表达如下：
[0055]
[0056] 其中D为差分算子，i、j分别为x方向和z方向的网格节点索引号，max为最大值函数，min为最小值函数。上式可简化为：
[0057]
[0058] 表示走时t在点(i,j)处x方向上的向前、向后的差分算子，表示走时t在点(i,j)处z方向上的向前、向后的差分算子。
[0059] 为保证计算精度，采用二阶上风差分算子。即：
[0060]
[0061]
[0062] 式中Δx、Δz分别为X、Z方向的网格尺寸。
[0063] 最终，利用目标点(i,j)周围八个节点(i,j±1)，(i,j±2)，(i±1,j)，(i±2,j)处的走时值，共有6种组合方式计算(i,j)点的走时值，6种组合方式的二阶差分计算公式如下：
[0064]
[0065] 式中，t是地震波走时，s是地震波的慢度，h＝Δx＝Δz为正方形网格尺寸。
[0066] 利用上式计算出多个走时，按照窄带扩展基本原理，在求出的多值解中选取最小走时点，即为目标点走时，计算公式如下：
[0067] ti,j＝min(t1i,j,t2i,j,t3i,j,t4i,j,t5i,j,t6i,j)
[0068] 遍历完所有节点，即可获得地震波时间场。
[0069] 优选的，所述波场抽取装置执行以下子步骤：
[0070] 抽取每个检波点的地震波场，即在上步中求得的时间场t(x,z)中找到检波点位置对应的时间，就是原震源点到检波点的地震走时，正演结束。
[0071] 因此，本发明具有如下优点：
[0072] 结合微震监测中震源点数量远远大于检波点数量的特点，依据地震波场模拟中源检互换原理，在现有技术基础上，将原震源点作为检波点、原检波点作为震源点，计算地震波场，最终获得每个源检对的地震波场信息。通过源检互换后，大大提高了微震数值模拟效率。附图说明
[0073] 图1为本发明的流程图。
[0074] 图2是实施例速度模型(二维)及观测系统。
具体实施方式
[0075] 下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。本发明的保护范围不受以下实例的限制。
[0076] 根据Vermeer提出的对称采样原理，地震勘探中源点(或炮点)和检波点互换，可以获得相同的地震信号，即在源点处激发-检波点处接受、检波点处激发-源点处接受，二者接收到的信号是完全对等的，被称为源检互换原理(即炮检互换原理)。
[0077] 本发明基于这一原理，首先根据待模拟对象的地质条件，建立网格速度模型，并指定检波点坐标；然后根据震源点和检波点个数确定波场计算方式，当震源点数小于检波点数时，基于源检互换原理，依次将各个检波点作为震源点，计算模型的地震波场；最后，基于上述计算得到的三维地震波场，抽取每个检波点(原震源点)的地震波场，正演结束。
[0078] 下面结合图1，进一步介绍本实施例的加速微震数值模拟方法。
[0079] 如图1所示，本发明采用如下技术方案：
[0080] 首先，给定正演模拟参数，即根据待模拟对象的地质条件，建立网格速度模型，并指定检波点坐标。
[0081] 具体的，在本实施例中根据待模拟区域建立XYZ三维直角坐标系；基于所述三维直角坐标系对待模拟区域进行网格划分；基于所述网格划分结果，对网格点进行速度赋值，建立速度模型；指定每个检波点的空间位置。其中，网格速度模型、检波点空间坐标及微震震源点分别记为：
[0082] 网格速度模型：vi,j,k(1≤i≤Nx,1≤j≤Ny,1≤k≤Nz)
[0083] 检波点坐标:(xri,yri,zrri)(1≤i≤Nr)
[0084] 震源点坐标：(xsi,ysi,zsri)(1≤i≤Ns)
[0085] 变量说明：vi,j,k为网格节点(i,j,k)对应的速度，Nx为模型X方向网格数，Ny为模型Y方向网格数，Nz为模型Z方向网格数，xri为第i个检波点的X坐标，yri为第i个检波点的Y坐标，zri为第i个检波点的Z坐标，Nr为检波点个数，xsi为第i个震源点的X坐标，ysi为第i个震源点的Y坐标，zsi为第i个震源点的Z坐标，Ns为震源点个数。
[0086] 优选的，所述波场计算模块执行以下子步骤：
[0087] 首先，根据震源点个数Ns和检波点个数Nr，确定地震波场计算方式。若Ns＞Nr，按常规遍历震源点，若Ns＜Nr，依次遍历各个检波点，将检波点当做震源点，原震源点标记为检波点。微震监测中，一般属于前者。
[0088] 其次，计算当前震源点，在已知速度模型下的地震波场。以快速推进算法为例，求解下式二维直角坐标系(x-z)下的程函方程：
[0089]
[0090] 其中t是地震波走时，s是慢度，即速度的倒数，x、z代表二维直角坐标系的2个坐标轴。
[0091] 采用上风差分格式代替式程函方程中的偏微分格式，即采用上风差分法对上式进行离散，其离散式可表达如下：
[0092]
[0093] 其中D为差分算子，i、j分别为x方向和z方向的网格节点索引号，max为最大值函数，min为最小值函数。上式可简化为：
[0094]
[0095] 表示走时t在点(i,j)处x方向上的向前、向后的差分算子，表示走时t在点(i,j)处z方向上的向前、向后的差分算子。
[0096] 为保证计算精度，采用二阶上风差分算子。即：
[0097]
[0098]
[0099] 式中Δx、Δz分别为X、Z方向的网格尺寸。
[0100] 最终，利用目标点(i,j)周围八个节点(i,j±1)，(i,j±2)，(i±1,j)，(i±2,j)处的走时值，共有6种组合方式计算(i,j)点的走时值，6种组合方式的二阶差分计算公式如下：
[0101]
[0102] 式中，t是地震波走时，s是地震波的慢度，h＝Δx＝Δz为正方形网格尺寸。
[0103] 利用上式计算出多个走时，按照窄带扩展基本原理，在求出的多值解中选取最小走时点，即为目标点走时，计算公式如下：
[0104] ti,j＝min(t1i,j,t2i,j,t3i,j,t4i,j,t5i,j,t6i,j)
[0105] 遍历完所有节点，即可获得地震波时间场。
[0106] 优选的，所述波场抽取步骤执行以下子步骤：
[0107] 接下来，抽取每个检波点的地震波场，即在上步中求得的时间场t(x,z)中找到检波点位置对应的时间，就是原震源点到检波点的地震走时，正演结束。
[0108] 本实施例极大提高了微震监测正演模拟的效率。在微震监测中，震源数量远远大于检波点数(Ns＞＞Nr),如在煤田水平井压裂微震监测中，检波点数(Nr)一般十几个，而微震事件(每个事件对应一个源点)(Ns)往往几百个甚至上千个，二者相差1～2个数量级。对这个过程进行正演模拟时，若按现有方法，求每个源点的地震波场，假设一个源点的模拟时间是T，则模拟整个过程所需时间为Ns×T。而按本专利所述方法，源检互换后，只需求每个检波点的地震波场，由于速度模型不变，因此每个检波点的模拟时间也是T，则模拟整个过程所需时间为Nr×T；两种方案的时差为(Ns-Nr)×T，可以看出，Ns与Nr差值越大，加速比越大。
[0109] 在具体实施时，需对比较Ns与Nr差值，若Ns＞Nr，则计算每个检波点地震波场，若Ns＜Nr，则计算每个震源点地震波场，流程如图1所示，微震监测正演模拟一般适于前者。
[0110] 本发明中所述正演方法与常规正演方法的对比效果。本例仅以二维情况说明，具体的物理问题描述为：边长为100m的正方形模型的4个顶点各放置1个检波点，共4个检波点，模型速度为常速度1000m/s，网格尺寸为0.1m*0.1m,网格大小为1001*1001，正方形内随机布置10个震源，计算每个源检对的走时。利用本发明方法运算时间为3058ms，常规方法运算时间为7898ms，加速比为2.58倍。
[0111] 本实施例仅以二维情况简单说明本发明对计算效率的提高，实际应用中均是三维模型，且震源数量远远大于检波点数量，加速比会更大。
| 标题 | 发布/更新时间 | 阅读量 |
|---|---|---|
| 红冲机自动调节定位抽芯组件 | 2021-03-07 | 5 |
| 一种超柔无人机线束 | 2022-02-25 | 5 |
| 2055A提升器内提升臂装置 | 2020-09-17 | 4 |
| 大树移栽高成活率技术 | 2020-06-01 | 3 |
| 充电装置及充电站 | 2021-06-24 | 3 |
| 粒状混凝土掺合料及其造粒方法 | 2021-10-14 | 2 |
| 一种轧钢生产工艺 | 2021-07-06 | 5 |
| 基于航空遥感的森林火灾探测方法 | 2020-07-19 | 5 |
| 预防破坏性大地震和火山大喷发的方法 | 2020-06-22 | 3 |
| 一种钢丝及其制造方法 | 2021-03-04 | 5 |
高效检索全球专利专利汇是专利免费检索,专利查询,专利分析-国家发明专利查询检索分析平台，是提供专利分析，专利查询，专利检索等数据服务功能的知识产权数据服务商。
我们的产品包含105个国家的1.26亿组数据，免费查、免费专利分析。
分析报告
专利汇分析报告产品可以对行业情报数据进行梳理分析，涉及维度包括行业专利基本状况分析、地域分析、技术分析、发明人分析、申请人分析、专利权人分析、失效分析、核心专利分析、法律分析、研发重点分析、企业专利处境分析、技术处境分析、专利寿命分析、企业定位分析、引证分析等超过60个分析角度，系统通过AI智能系统对图表进行解读，只需1分钟，一键生成行业专利分析报告。
现有技术热门专利
QQ群二维码
意见反馈
意见反馈