CN112904419A - 一种微地震成像方法及终端设备 - Google Patents

一种微地震成像方法及终端设备 Download PDF

Info

Publication number
CN112904419A
CN112904419A CN202110104645.XA CN202110104645A CN112904419A CN 112904419 A CN112904419 A CN 112904419A CN 202110104645 A CN202110104645 A CN 202110104645A CN 112904419 A CN112904419 A CN 112904419A
Authority
CN
China
Prior art keywords
model
probability
velocity
updated
new
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
CN202110104645.XA
Other languages
English (en)
Other versions
CN112904419B (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.)
Southwest University of Science and Technology
Original Assignee
Southwest University of Science and 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 Southwest University of Science and Technology filed Critical Southwest University of Science and Technology
Priority to CN202110104645.XA priority Critical patent/CN112904419B/zh
Publication of CN112904419A publication Critical patent/CN112904419A/zh
Application granted granted Critical
Publication of CN112904419B publication Critical patent/CN112904419B/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. 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/288—Event detection in seismic signals, e.g. microseismics
    • 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/303—Analysis for determining velocity profiles or travel times

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • Acoustics & Sound (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Business, Economics & Management (AREA)
  • Emergency Management (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Measurement Of Velocity Or Position Using Acoustic Or Ultrasonic Waves (AREA)

Abstract

本发明公开了一种微地震成像方法及终端设备，其中，方法包括步骤：基于声波测井数据建立一维速度模型；采用基于贝叶斯理论和可逆跳马尔科夫链蒙特卡罗算法的变维方法对速度模型结构和微地震事件位置同时进行约束反演，对所述一维速度模型进行更新校正，得到更新速度模型；基于所述更新速度模型进行微地震成像。本发明充分利用了微地震事件空间分布广泛，对速度模型约束能力好的特点，在更新速度模型的同时获得了微地震事件的精确位置，一举两得的解决了速度模型校正和微地震成像的问题。

Description

一种微地震成像方法及终端设备
技术领域
本发明涉及微地震成像技术领域，尤其涉及一种微地震成像方法及终端设备。
背景技术
利用微地震监测技术指导水力压裂开采油气资源已经成为一种重要的实用方法和商业手段。微地震监测过程中速度模型对成像精度具有重要影响。由于获得的地质信息有限，传统井中微地震监测速度模型往往设置为一维层状结构，进而利用射孔对层内速度值进行校正以达到精确定位射孔的目的。然而该方法存在诸多缺陷：1.水力压裂阶段每个压裂段位可以提供的射孔个数有限，同时井中监测可以利用的检波器个数及空间方位也受到制约，这导致只有介于射孔到检波器之间的速度模型受到约束，远离该约束段的速度模型无法得到校正，发生在偏离校正段的微地震事件或出现较大的定位误差；2.在一些较复杂的压裂区域，由于前期的沉积构造作用导致真实的地层并非水平，同时有可能存在断层等结构，如果速度模型校正过程中层位个数和深度保持不变，获得的更新模型可能仍然无法精确定位水力压裂产生的微地震事件位置。因此，必须要增加更多的有效信息获得更合理的速度模型，从而可以达到微地震监测水力压裂精确成像的目的。
因此，现有技术还有待于改进和发展。
发明内容
鉴于上述现有技术的不足，本发明的目的在于提供一种微地震成像方法及终端设备，旨在解决基于现有技术更新的速度模型无法达到微地震监测水力压裂精确成像的目的问题。
本发明的技术方案如下：
一种微地震成像方法，其中，包括步骤：
基于声波测井数据建立一维速度模型；
采用基于贝叶斯理论和可逆跳马尔科夫链蒙特卡罗算法的变维方法对速度模型结构和微地震事件位置同时进行约束反演，对所述一维速度模型进行更新校正，得到更新速度模型；
基于所述更新速度模型进行微地震成像。
所述的微地震成像方法，其中，所述基于声波测井数据建立一维速度模型的步骤包括：
获得研究区域的声波测井曲线，根据所述声波测井曲线表征的速度值大小按照深度划分不同的层位，并取同一层位的声波测井数值的平均值作为建立的一维速度模型的速度值，获得该研究区域的一维速度模型。
所述的微地震成像方法，其中，采用基于贝叶斯理论和可逆跳马尔科夫链蒙特卡罗算法的变维方法对速度模型结构和微地震事件位置同时进行约束反演的步骤包括：
确定反演模型的模型参数表达式为：m＝[n，D，VP，VS，H，Z]，其中，n表示反演模型层位个数，D表示每层的深度，VP＝[VP1，VP2，...，VPn]T和VS＝[VS1，VS2，...，VSn]T分别表示每层的P波速度值和S波速度值，H＝[H1，H2，...，Hk]T和Z＝[Z1，Z2，...，Zk]T分别表示微地震在二维定位情况下的水平位置和垂直位置；
确定贝叶斯理论表达式为
Figure BDA0002916883930000021
其中，d表示观测数据，m是模型参数矢量，p(m)是先验模型信息，p(d|m)是似然函数，p(m|d)是后验模型概率，p(d)是观测数据在模型空间中的整体概率，为一常数；
基于所述贝叶斯理论表达式，利用观测数据来预测各模型参数的概率，选取概率较大的模型参数构建新的速度模型和作为微地震事件定位的位置。
所述的微地震成像方法，其中，所述基于所述贝叶斯理论表达式，利用观测数据来预测各模型参数的概率，选取概率较大的模型参数构建新的速度模型和作为微地震事件定位的位置的步骤包括：
在贝叶斯理论中，采用先验信息代表所有模型参数已知信息的总和，其表达式为：p(m)＝p(Z)p(H)p(Vs|n，D)p(VP|n，D)p(D|n)p(n)，其中，p(n)表示所有的有可能层位个数的概率，p(D|n)表示在层位个数n下层位深度的概率，p(VP|n,D)和p(VS|n,D)表示在层位个数为n，层位深度为D下P波和S波速度模型的概率，p(H)和p(Z)表示在给定范围内微地震事件位置的概率；
将所有的先验信息按照均匀分布或者高斯分布设计，基于所述先验信息使反演快速收敛到最优解。
所述的微地震成像方法，其中，所述模型层位个数n是一个变量，所有可能的结果服从均匀分布：p(n)＝1/Δn，其中，Δn＝(nmax-nmin).nmax和nmin代表着最大和最小的可能层位个数；
在n层模型中，深度D用概率表示为：
Figure BDA0002916883930000031
N代表着所有的可能层位深度；
第i层的P波速度值用概率表示为：
Figure BDA0002916883930000032
其中，ΔvP＝(vmax-vmin)P；
第i层的S波速度值用概率表示为：
Figure BDA0002916883930000033
其中，ΔvS＝(vmax-vmin)S；
微地震事件的位置[H,Z]用概率表示为：
Figure BDA0002916883930000034
其中，Δh＝(hmax-hmin)，Δz＝(zmax-zmin)。
所述的微地震成像方法，其中，基于贝叶斯理论和可逆跳马尔科夫链蒙特卡罗算法对所述一维速度模型进行更新校正的步骤包括：
采用可逆跳马尔科夫链蒙特卡罗算法迭代产生后验模型，在每次迭代过程中，一些参数会得到更新从而产生一个新的模型，所述新的模型被用来计算后验似然函数值，然后产生接收概率：
Figure BDA0002916883930000041
其中，mold表示更新前的模型，mnew表示更新后的模型,p(mnew)和p(d|mnew)表示更新模型的先验信息及其似然函数，p(mold)和p(d|mold)分别是更新前的模型的先验信息及似然函数,q(mnew|mold)是更新前模型转换为更新后模型的概率，q(mold|mnew)是更新后模型转换为更新前模型的概率,J是从更新前模型转换为更新后模型的雅克比转换矩阵；
将计算的接收概率α(mnew|mold)与一个服从均匀分布[0,1]的随机数r相对比，如果α≥r,更新后的模型mnew将会被接受，如果α＜r,更新的模型被拒绝，现在的模型mold将会进入下一个循环。
所述的微地震成像方法，其中，基于贝叶斯理论和可逆跳马尔科夫链蒙特卡罗算法对所述一维速度模型进行更新校正的步骤还包括：
输入模型参数mj；
产生随机数a，如果a是奇数，则选择更新速度参数，如果a是偶数，则选择更新微地震事件的位置；
计算接收概率α(mnew|mj)，mnew是新产生的模型参数，如果新的模型被接受，则第j+1个模型mj+1＝mnew,否则,mj+1＝mj。
所述的微地震成像方法，其中，对于更新速度参数，其包括出生、死亡、移动以及改变这四种选择，其中，
出生选择是指随机的产生一个层位，其界面深度服从概率分布：
Figure BDA0002916883930000051
死亡选择是指随机的在已有的层位中选择一个并且删除它，选择概率为
Figure BDA0002916883930000052
移动选择是指随机的在已有的层位中选择一个层位，并且扰动它的深度，扰动概率为
Figure BDA0002916883930000053
改变选择是指随机的选择一个层位中的P波速度值或S波速度值，改变该速度值的大小：
Figure BDA0002916883930000054
其中，N是所有的层位个数,n是模型
Figure BDA0002916883930000055
的层位个数，D是更新的层位深度值，Dj是目前的层位深度值，v是更新的速度值，vj是目前的速度值，u是服从均匀分布[0,1]的一个随机数，σ1和σ2是深度扰动或者速度扰动的标准差。
所述的微地震成像方法，其中，对于微地震事件位置更新，每次选中一个微地震事件，它的位置[H,Z]各被以1/2的概率扰动，扰动的函数表示为：
Figure BDA0002916883930000056
H＝Hi+g×σ3
Figure BDA0002916883930000057
Z＝Zi+g×σ4，其中，H是更新后的水平位置，Hi是更新前的水平位置，Z是更新后的垂直位置，Zi是更新前的垂直位置，g是一个服从均匀分布[0,1]的随机变量，σ3和σ4是水平位置扰动和垂直位置扰动函数的标准差。
一种终端设备，其中，包括处理器，适于实现各指令；以及存储介质，适于存储多条指令，所述指令适于由处理器加载并执行本发明所述的微地震成像方法中的步骤。
有益效果：本发明提出变维思想更新速度模型，本质上是不再依赖测井等一孔之见限制速度模型的结构，本发明在已知地质层位先验信息较少的情况下，利用射孔和微地震事件更新速度模型，同时反演校正微地震事件的位置，最终获得的速度模型不能反映真实的地质信息，但是可以作为等效速度模型对该区域水力压裂产生的微地震事件进行精确定位。因此，该方法不被较片面的先验信息所误导，保证了微地震检测的基本目的以及水力压裂精确成像任务的完成。
附图说明
图1为本发明一种微地震成像方法较佳实施例的流程图。
图2为本发明提供的一种微地震成像方法中同时校正速度模型结构和微地震事件位置的流程图。
图3为本发明终端设备的原理框图。
具体实施方式
本发明提供一种微地震成像方法及终端设备，为使本发明的目的、技术方案及效果更加清楚、明确，以下对本发明进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
本技术领域技术人员可以理解，除非特意声明，这里使用的单数形式“一”、“一个”、“所述”和“该”也可包括复数形式。应该进一步理解的是，本发明的说明书中使用的措辞“包括”是指存在所述特征、整数、步骤、操作、元件和/或组件，但是并不排除存在或添加一个或多个其他特征、整数、步骤、操作、元件、组件和/或它们的组。应该理解，当我们称元件被“连接”或“耦接”到另一元件时，它可以直接连接或耦接到其他元件，或者也可以存在中间元件。此外，这里使用的“连接”或“耦接”可以包括无线连接或无线耦接。这里使用的措辞“和/或”包括一个或更多个相关联的列出项的全部或任一单元和全部组合。
本技术领域技术人员可以理解，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)，具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语，应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样被特定定义，否则不会用理想化或过于正式的含义来解释。
微地震是一种小型的地震(mine tremor or microseismic)，在地下矿井深部开采过程中发生岩石破裂和地震活动，常常是不可避免的现象。由开采诱发的地震活动，微地震事件通常定义为在开采坑道附近的岩体内因应力场变化导致岩石破坏而引起的那些地震事件。在一些较复杂的压裂区域，由于前期的沉积构造作用导致真实的地层并非水平，同时有可能存在断层等结构，如果速度模型校正过程中层位个数和深度保持不变，获得的更新模型可能仍然无法精确定位水力压裂产生的微地震事件位置，从而导致无法精确成像。
基于此，本发明提供了一种微地震成像方法，如图1所示，其包括步骤：
S10、基于声波测井数据建立一维速度模型；
S20、采用基于贝叶斯理论和可逆跳马尔科夫链蒙特卡罗算法的变维方法对速度模型结构和微地震事件位置同时进行约束反演，对所述一维速度模型进行更新校正，得到更新速度模型；
S30、基于所述更新速度模型进行微地震成像。。
与常规方法中利用声波测井构建初始速度模型，继而利用射孔更新速度值以达到准确定位射孔的速度校正方法不同，本实施例对声波测井的依赖程度较低，利用射孔和微地震事件同时对速度模型结构和微地震事件位置进行约束反演，达到获得的等效速度模型能够精确定位射孔以及微地震事件位置的目的。本实施例利用变维(trans-dimensional)反演思想对速度模型结构和微地震事件位置同时进行更新，利用Bayesian inference(贝叶斯理论)概率思想和可逆跳马尔科夫链蒙特卡罗(reversible jump Monte MarkovMonte Carlo，rjMCMC)算法对构建的一维速度模型进行更新校正，允许产生不同层位的模型以达到准确匹配观测数据的目的。在反演过程中，射孔和微地震事件位置同时更新，有效降低了速度模型和微地震位置“耦合”现象的产生。同时本实施例充分利用了微地震事件空间分布广泛，对速度模型约束能力好的特点，在更新速度模型的同时获得了微地震事件的精确位置，一举两得的解决了速度模型校正和微地震成像的问题。
在一些实施方式中，利用射孔信息和遍布于水力压裂各个位置的微地震事件同时校正速度模型结构和确定微地震事件位置，确定反演模型的模型参数表达式为：m＝[n，D，VP，VS，H，Z]，其中，n表示反演模型层位个数，D表示每层的深度，VP＝[VP1，VP2，...，VPn]T和VS＝[VS1，VS2，...，VSn]T分别表示每层的P波速度值和S波速度值，H＝[H1，H2，...，Hk]T和Z＝[Z1，Z2，...，Zk]T分别表示微地震在二维定位情况下的水平位置和垂直位置。本实施例利用贝叶斯理论作为反演方法实现对所述反演模型的变维推导和所述微地震事件位置的精确定位，获得一维速度模型。
具体来讲，确定贝叶斯理论表达式为
Figure BDA0002916883930000081
其中，d表示观测数据，m是模型参数矢量，p(m)是先验模型信息，p(d|m)是似然函数，p(m|d)是后验模型概率，p(d)是观测数据在模型空间中的整体概率，为一常数；基于所述贝叶斯理论表达式，利用观测数据来预测各模型参数的概率，选取概率较大的模型参数构建新的速度模型和作为微地震事件定位的位置。
在本实施例中，所述贝叶斯理论表达式可以写为：p(m|d)∝p(d|m)p(m)，由此可见，后验模型概率仅仅与先验模型信息和似然函数有关。
似然函数p(d|m)通过对比计算数据和观测数据，提供了一种定量评价反演模型与真实模型相似度的度量标准。它与观测目标函数Φ(m)有关，其关系表达式为：
Figure BDA0002916883930000091
其中，gt(m)为依据模型参数m正演计算的数据，Ct为噪声的协方差矩阵，T为转秩，计算的数据与观测数据拟合程度越好，则认为反演模型越接近需求的模型。
在贝叶斯理论中，采用先验信息代表所有模型参数已知信息的总和，其表达式为：p(m)＝p(Z)p(H)p(VS|n，D)p(VP|n，D)p(D|n)p(n)，其中，p(n)表示所有的有可能层位个数的概率，p(D|n)表示在层位个数n下层位深度的概率，p(VP|n,D)和p(VS|n,D)表示在层位个数为n，层位深度为D下P波和S波速度模型的概率，p(H)和p(Z)表示在给定范围内微地震事件位置的概率；为了防止错误的先验信息干扰反演结果，将所有的先验信息按照均匀分布或者高斯分布设计，基于所述先验信息使反演快速收敛到最优解。
在本实施例中，所述模型层位个数n是一个变量，所有可能的结果服从均匀分布：p(n)＝1/Δn，其中，Δn＝(nmax-nmin).nmax和nmin代表着最大和最小的可能层位个数；
在n层模型中，深度D用概率表示为：
Figure BDA0002916883930000101
N代表着所有的可能层位深度；
不同层位的速度值可以相同，因此在第i层的P波和S波速度值可以用概率表示为p(Vi)＝1/Δv，其中，Δv＝(vmax-vmin)，vmax和vmin是最大和最小的可能的速度值；
对于所有的n个层位：
Figure BDA0002916883930000102
因此，第i层的P波速度值用概率表示为：
Figure BDA0002916883930000103
其中，ΔvP＝(vmax-vmin)P；
第i层的S波速度值用概率表示为：
Figure BDA0002916883930000104
其中，ΔvS＝(vmax-vmin)S；
微地震事件的位置[H,Z]用概率表示为：
Figure BDA0002916883930000105
其中，Δh＝(hmax-hmin)，Δz＝(zmax-zmin)，hmax和hmin,zmax和zmin是水平方向和垂直方向中最大和最小的位置。
最终，模型m所有参数的先验信息可以表示为：
Figure BDA0002916883930000106
为了获得准确的后验概率密度，均匀采样方法由于计算模型数量太多，导致计算效率较低。因此，我们采用了马尔科夫链蒙特卡罗方法(Markov chain Monte Carlo，MCMC)来计算后验模型概率。由于需要在不同维度之间变换，我们选择采用reversible jump MCMC(rjMCMC)算法迭代产生后验模型。在每次迭代过程中，一些参数会得到更新从而产生一个新的模型，该模型会被用来计算后验似然函数值，然后产生接收概率：
Figure BDA0002916883930000111
其中，mold表示更新前的模型，mnew表示更新后的模型,p(mnew)和p(d|mnew)表示更新模型的先验信息及其似然函数，p(mold)和p(d|mold)分别是更新前的模型的先验信息及似然函数,q(mnew|mold)是更新前模型转换为更新后模型的概率，q(mold|mnew)是更新后模型转换为更新前模型的概率,J是从更新前模型转换为更新后模型的雅克比转换矩阵；将计算的接收概率α(mnew|mold)与一个服从均匀分布[0,1]的随机数r相对比，如果α≥r,更新后的模型mnew将会被接受，如果α＜r,更新的模型被拒绝，现在的模型mold将会进入下一个循环。
在rjMCMC反演过程中，前面产生的数十个反演模型将会被忽略，原因是避免受初始模型的影响。剩余的模型，将会作为选取样本构成后验概率密度。只要给与足够的迭代次数，反演获得的模型参数计算的后验概率模型就会非常接近真实的后验概率密度。利用该后验概率密度的一些统计学参数，例如平均值和方差，就可以评价最终的速度模型。微地震事件联合反演的位置的标准差，可以作为衡量定位精确性的标准。
在一些实施方式中，在计算反演模型的参数过程中，我们利用rjMCMC算法在不同的层位之间进行跳跃，不同层位速度结构的模型均会被评价是否可以接受，最终高概率的模型会被接受。如图2所示，详细的反演策略可以表示为：
输入模型参数mj；
产生随机数a，如果a是奇数，则选择更新速度参数，如果a是偶数，则选择更新微地震事件的位置；
计算接收概率α(mnew|mj)，mnew是新产生的模型参数，如果新的模型被接受，则第j+1个模型mj+1＝mnew,否则,mj+1＝mj。
在本实施例中，对于更新速度参数，其包括出生、死亡、移动以及改变这四种选择，其中，
出生选择是指随机的产生一个层位，其界面深度服从概率分布：
Figure BDA0002916883930000121
死亡选择是指随机的在已有的层位中选择一个并且删除它，选择概率为
Figure BDA0002916883930000122
移动选择是指随机的在已有的层位中选择一个层位，并且扰动它的深度，扰动概率为
Figure BDA0002916883930000123
改变选择是指随机的选择一个层位中的P波速度值或S波速度值，改变该速度值的大小：
Figure BDA0002916883930000124
其中，N是所有的层位个数,n是模型
Figure BDA0002916883930000125
的层位个数，D是更新的层位深度值，Dj是目前的层位深度值，v是更新的速度值，vj是目前的速度值，u是服从均匀分布[0,1]的一个随机数，σ1和σ2是深度扰动或者速度扰动的标准差。
在本实施例中，对于微地震事件位置更新，每次选中一个微地震事件，它的位置[H,Z]各被以1/2的概率扰动，扰动的函数表示为：
Figure BDA0002916883930000131
H＝Hi+g×σ3
Figure BDA0002916883930000132
Z＝Zi+g×σ4，其中，H是更新后的水平位置，Hi是更新前的水平位置，Z是更新后的垂直位置，Zi是更新前的垂直位置，g是一个服从均匀分布[0,1]的随机变量，σ3和σ4是水平位置扰动和垂直位置扰动函数的标准差。
本实施例提供了同时更新速度结构和微地震事件位置的流程，同时通过利用rjMCMC算法，反演的模型层位可以在不同的纬度之间跳跃，最终可以获得能够拟合观测数据的后验模型。
本发明将变维思想引入井中微地震检测速度模型校正过程中，不仅仅依靠声波测井数据来构建速度模型，同时利用射孔和微地震时间约束速度模型层位信息。在模型更新过程中，不限制模型的层位个数，允许产生不同层位的速度模型匹配观测数据。当更新的速度模型可以拟合观测数据时，该模型既可作为等效模型对微地震事件进行定位，没有必要反映真实的速度模型结构信息。本发明将速度模型结构和微地震事件位置同时进行更新校正，即利用了微地震事件空间分布广泛，射线覆盖程度高，对速度约束好的特点，同时也更新了微地震事件的位置，保证了在优化速度模型结构的同时对微地震事件进行精确定位。
在一些实施方式中，还提供一种存储介质，其中，所述存储介质存储有一个或者多个程序，所述一个或者多个程序可被一个或者多个处理器执行，以实现如本发明所述的微地震成像方法中的步骤。
在一些实施方式中，还提供一种终端设备，如图3所示，其包括至少一个处理器(processor)20；显示屏21；以及存储器(memory)22，还可以包括通信接口(CommunicationsInterface)23和总线24。其中，处理器20、显示屏21、存储器22和通信接口23可以通过总线24完成相互间的通信。显示屏21设置为显示初始设置模式中预设的用户引导界面。通信接口23可以传输信息。处理器20可以调用存储器22中的逻辑指令，以执行上述实施例中的方法。
此外，上述的存储器22中的逻辑指令可以通过软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。
存储器22作为一种计算机可读存储介质，可设置为存储软件程序、计算机可执行程序，如本公开实施例中的方法对应的程序指令或模块。处理器20通过运行存储在存储器22中的软件程序、指令或模块，从而执行功能应用以及数据处理，即实现上述实施例中的方法。
存储器22可包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需的应用程序；存储数据区可存储根据终端设备的使用所创建的数据等。此外，存储器22可以包括高速随机存取存储器，还可以包括非易失性存储器。例如，U盘、移动硬盘、只读存储器(Read-Only Memory，ROM)、随机存取存储器(Random Access Memory，RAM)、磁碟或者光盘等多种可以存储程序代码的介质，也可以是暂态存储介质。
此外，上述存储介质以及终端设备中的多条指令处理器加载并执行的具体过程在上述方法中已经详细说明，在这里就不再一一陈述。
综上所述，本发明提供的变维反演模型考虑到真实的地质岩层存在倾斜、断层等现象，不再仅仅依赖声波测井构建速度模型，摒弃了传统方法中仅仅校正P波速度值和S波速度值，而对层位信息鲜有更新的现象。本发明充分利用了射孔和大量微地震事件的约束，在校正速度模型的同时允许速度模型结构发生改变，速度层位个数可以获得增加或减少，速度层位个数也可以获得改变。虽然校正获得的一维速度模型不再能精确反映真实的地质层位信息，但是可以作为等效速度模型获得精确的微地震事件位置，保证了水力压裂成像的准确性。该速度模型同时可以为后续二维以及三维复杂模型的反演提供依据，在反演过程中作为初始模型，保证了微地震事件定位的精确性。
应当理解的是，本发明的应用不限于上述的举例，对本领域普通技术人员来说，可以根据上述说明加以改进或变换，所有这些改进和变换都应属于本发明所附权利要求的保护范围。

Claims (10)

1.一种微地震成像方法，其特征在于，包括步骤：
基于声波测井数据建立一维速度模型；
采用基于贝叶斯理论和可逆跳马尔科夫链蒙特卡罗算法的变维方法对速度模型结构和微地震事件位置同时进行约束反演，对所述一维速度模型进行更新校正，得到更新速度模型；
基于所述更新速度模型进行微地震成像。
2.根据权利要求1所述的微地震成像方法，其特征在于，所述基于声波测井数据建立一维速度模型的步骤包括：
获得研究区域的声波测井曲线，根据所述声波测井曲线表征的速度值大小按照深度划分不同的层位，并取同一层位的声波测井数值的平均值作为建立的一维速度模型的速度值，获得该研究区域的一维速度模型。
3.根据权利要求2所述的微地震成像方法，其特征在于，采用基于贝叶斯理论和可逆跳马尔科夫链蒙特卡罗算法的变维方法对速度模型结构和微地震事件位置同时进行约束反演的步骤包括：
确定反演模型的模型参数表达式为：m＝[n，D，VP，VS，H，Z]，其中，n表示反演模型层位个数，D表示每层的深度，VP＝[VP1，VP2，...，VPn]T和VS＝[VS1，VS2，...，VSn]T分别表示每层的P波速度值和S波速度值，H＝[H1，H2，...，Hk]T和Z＝[Z1，Z2，...，Zk]T分别表示微地震在二维定位情况下的水平位置和垂直位置；
确定贝叶斯理论表达式为
Figure FDA0002916883920000011
其中，d表示观测数据，m是模型参数矢量，p(m)是先验模型信息，p(d|m)是似然函数，p(m|d)是后验模型概率，p(d)是观测数据在模型空间中的整体概率，为一常数；
基于所述贝叶斯理论表达式，利用观测数据来预测各模型参数的概率，选取概率较大的模型参数构建新的速度模型和作为微地震事件定位的位置。
4.根据权利要求3所述的微地震成像方法，其特征在于，所述基于所述贝叶斯理论表达式，利用观测数据来预测各模型参数的概率，选取概率较大的模型参数构建新的速度模型和作为微地震事件定位的位置的步骤包括：
在贝叶斯理论中，采用先验信息代表所有模型参数已知信息的总和，其表达式为：p(m)＝p(Z)p(H)p(VS|n，D)p(VP|n，D)p(D|n)p(n)，其中，p(n)表示所有的有可能层位个数的概率，p(D|n)表示在层位个数n下层位深度的概率，p(VP|n，D)和p(VS|n，D)表示在层位个数为n，层位深度为D下P波和S波速度模型的概率，p(H)和p(Z)表示在给定范围内微地震事件位置的概率；
将所有的先验信息按照均匀分布或者高斯分布设计，基于所述先验信息使反演快速收敛到最优解。
5.根据权利要求4所述的微地震成像方法，其特征在于，所述模型层位个数n是一个变量，所有可能的结果服从均匀分布：p(n)＝1/Δn，其中，Δn＝(nmax-nmin)·nmax和nmin代表着最大和最小的可能层位个数；
在n层模型中，深度D用概率表示为：
Figure FDA0002916883920000021
N代表着所有的可能层位深度；
第i层的P波速度值用概率表示为：
Figure FDA0002916883920000022
其中，ΔvP＝(vmax-vmin)P；
第i层的S波速度值用概率表示为：
Figure FDA0002916883920000023
其中，ΔvS＝(vmax-vmin)S；
Figure FDA0002916883920000031
微地震事件的位置[H，Z]用概率表示为：
Figure FDA0002916883920000032
其中，
Figure FDA0002916883920000034
Δz＝(zmax-zmin)。
6.根据权利要求5所述的微地震成像方法，其特征在于，基于贝叶斯理论和可逆跳马尔科夫链蒙特卡罗算法对所述一维速度模型进行更新校正的步骤包括：
采用可逆跳马尔科夫链蒙特卡罗算法迭代产生后验模型，在每次迭代过程中，一些参数会得到更新从而产生一个新的模型，所述新的模型被用来计算后验似然函数值，然后产生接收概率：
Figure FDA0002916883920000033
其中，mold表示更新前的模型，mnew表示更新后的模型，p(mnew)和p(d|mnew)表示更新模型的先验信息及其似然函数，p(mold)和p(d|mold)分别是更新前的模型的先验信息及似然函数，q(mnew|mold)是更新前模型转换为更新后模型的概率，q(mold|mnew)是更新后模型转换为更新前模型的概率，J是从更新前模型转换为更新后模型的雅克比转换矩阵；
将计算的接收概率α(mnew|mold)与一个服从均匀分布[0，1]的随机数r相对比，如果α≥r，更新后的模型mnew将会被接受，如果α＜r，更新的模型被拒绝，现在的模型mold将会进入下一个循环。
7.根据权利要求6所述的微地震成像方法，其特征在于，基于贝叶斯理论和可逆跳马尔科夫链蒙特卡罗算法对所述一维速度模型进行更新校正的步骤还包括：
输入模型参数mj；
产生随机数a，如果a是奇数，则选择更新速度参数，如果a是偶数，则选择更新微地震事件的位置；
计算接收概率α(mnew|mj)，mnew是新产生的模型参数，如果新的模型被接受，则第j+1个模型mj+1＝mnew，否则，mj+1＝mj。
8.根据权利要求7所述的微地震成像方法，其特征在于，对于更新速度参数，其包括出生、死亡、移动以及改变这四种选择，其中，
出生选择是指随机的产生一个层位，其界面深度服从概率分布：
Figure FDA0002916883920000041
死亡选择是指随机的在已有的层位中选择一个并且删除它，选择概率为
Figure FDA0002916883920000042
移动选择是指随机的在已有的层位中选择一个层位，并且扰动它的深
Figure FDA0002916883920000043
度，扰动概率为D＝Dj+u×σ1，；
改变选择是指随机的选择一个层位中的P波速度值或S波速度值，改
Figure FDA0002916883920000044
变该速度值的大小：v＝vj+u×σ2，，其中，N是所有的层位个数，n是模型
Figure FDA0002916883920000045
的层位个数，D是更新的层位深度值，Dj是目前的层位深度值，v是更新的速度值，vj是目前的速度值，u是服从均匀分布[0，1]的一个随机数，σ1和σ2是深度扰动或者速度扰动的标准差。
9.根据权利要求7所述的微地震成像方法，其特征在于，对于微地震事件位置更新，每次选中一个微地震事件，它的位置[H，Z]各被以1/2的概率扰动，扰动的函数表示为：
Figure FDA0002916883920000051
H＝Hi+g×σ3
Figure FDA0002916883920000052
Z＝Zi+g×σ4，其中，H是更新后的水平位置，Hi是更新前的水平位置，Z是更新后的垂直位置，Zi是更新前的垂直位置，g是一个服从均匀分布[0，1]的随机变量，σ3和σ4是水平位置扰动和垂直位置扰动函数的标准差。
10.一种终端设备，其特征在于，包括处理器，适于实现各指令；以及存储介质，适于存储多条指令，所述指令适于由处理器加载并执行权利要求1-9任意一项所述的微地震成像方法中的步骤。
CN202110104645.XA 2021-01-26 2021-01-26 一种微地震成像方法及终端设备 Active CN112904419B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110104645.XA CN112904419B (zh) 2021-01-26 2021-01-26 一种微地震成像方法及终端设备

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110104645.XA CN112904419B (zh) 2021-01-26 2021-01-26 一种微地震成像方法及终端设备

Publications (2)

Publication Number Publication Date
CN112904419A true CN112904419A (zh) 2021-06-04
CN112904419B CN112904419B (zh) 2023-01-13

Family

ID=76120281

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110104645.XA Active CN112904419B (zh) 2021-01-26 2021-01-26 一种微地震成像方法及终端设备

Country Status (1)

Country Link
CN (1) CN112904419B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115730424A (zh) * 2022-10-17 2023-03-03 南方科技大学 基于多源大地测量数据的有限断层反演方法、装置及终端
CN116466390A (zh) * 2023-02-17 2023-07-21 南方海洋科学与工程广东省实验室(广州) 一种大型水库诱发地震实时监测定位方法
CN117075212A (zh) * 2023-10-16 2023-11-17 吉林大学 一种隧道磁共振裂隙结构成像方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105589100A (zh) * 2014-10-21 2016-05-18 中国石油化工股份有限公司 一种微地震震源位置和速度模型同时反演方法
CN107703540A (zh) * 2017-06-26 2018-02-16 河海大学 一种微地震定位及层析成像方法
CN111736208A (zh) * 2020-06-24 2020-10-02 重庆大学 变权重联合P波和S波初至数据的微震事件Bayes定位方法、系统及介质
CN112213768A (zh) * 2020-09-25 2021-01-12 南方科技大学 一种联合震源机制反演的地面微地震定位方法及系统

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105589100A (zh) * 2014-10-21 2016-05-18 中国石油化工股份有限公司 一种微地震震源位置和速度模型同时反演方法
CN107703540A (zh) * 2017-06-26 2018-02-16 河海大学 一种微地震定位及层析成像方法
CN111736208A (zh) * 2020-06-24 2020-10-02 重庆大学 变权重联合P波和S波初至数据的微震事件Bayes定位方法、系统及介质
CN112213768A (zh) * 2020-09-25 2021-01-12 南方科技大学 一种联合震源机制反演的地面微地震定位方法及系统

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
T. BODIN,ET AL.: "Transdimensional inversion of receiver functions and surface wave dispersion", 《JOURNAL OF GEOPHYSICAL RESEARCH》 *
THOMAS BODIN AND MALCOLM SAMBRIDGE: "Seismic tomography with the reversible jump algorithm", 《GEOPHYS. J. INT.》 *
XIAO TIAN,ET AL.: "Cross double-difference inversion for microseismic event location using data from a single monitoring well", 《GEOPHYSICS》 *
XIAO TIAN,ET AL.: "Cross double-difference inversion for simultaneous velocity model update and microseismic event location", 《GEOPHYSICAL PROSPECTING》 *
李承瑾等: "跨维贝叶斯反演在地球物理中的研究进展", 《工程地球物理学报》 *
田宵: "井下微地震监测方法研究", 《中国博士学位论文全文数据库 工程科技Ⅰ辑》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115730424A (zh) * 2022-10-17 2023-03-03 南方科技大学 基于多源大地测量数据的有限断层反演方法、装置及终端
CN115730424B (zh) * 2022-10-17 2023-08-04 南方科技大学 基于多源大地测量数据的有限断层反演方法、装置及终端
CN116466390A (zh) * 2023-02-17 2023-07-21 南方海洋科学与工程广东省实验室(广州) 一种大型水库诱发地震实时监测定位方法
CN116466390B (zh) * 2023-02-17 2023-11-03 南方海洋科学与工程广东省实验室(广州) 一种大型水库诱发地震实时监测定位方法
CN117075212A (zh) * 2023-10-16 2023-11-17 吉林大学 一种隧道磁共振裂隙结构成像方法
CN117075212B (zh) * 2023-10-16 2024-01-26 吉林大学 一种隧道磁共振裂隙结构成像方法

Also Published As

Publication number Publication date
CN112904419B (zh) 2023-01-13

Similar Documents

Publication Publication Date Title
CN112904419B (zh) 一种微地震成像方法及终端设备
Jin et al. Multi-objective optimization-based updating of predictions during excavation
AU2013325186B2 (en) Propagating fracture plane updates
EP2491516B1 (en) Method for optimization with gradient information
US10519766B2 (en) Reservoir modelling with multiple point statistics from a non-stationary training image
EP2880592B1 (en) Multi-level reservoir history matching
AU2012369158B2 (en) Systems and methods for selecting facies model realizations
US10620339B2 (en) Static earth model calibration methods and systems using tortuosity evaluations
FR3043227A1 (zh)
US20230384470A1 (en) Method for three-dimensional velocity geological modeling with structures and velocities randomly arranged
Christou et al. Effective sampling of spatially correlated intensity maps using hazard quantization: Application to seismic events
CN105706105A (zh) 静态地球模型网格单元缩放和性质重新取样方法及系统
CN113945966B (zh) 人工压裂裂缝网络构建方法及装置
WO2019231681A1 (en) Inverse stratigraphic modeling using a hybrid linear and nonlinear algorithm
CN112925011B (zh) 一种单井微地震监测方法、存储介质及终端设备
US20150193707A1 (en) Systems and Methods for Estimating Opportunity in a Reservoir System
CN113987095A (zh) 一种面向地下非法采矿识别的gis时空数据模型
US11125905B2 (en) Methods for automated history matching utilizing muon tomography
CN113219534B (zh) 一种叠前深度偏移速度质控方法、装置、介质及电子设备
CN117572498A (zh) 地震初至波拾取方法、装置、存储介质及计算机设备
CN113219533A (zh) 一种叠前时间偏移速度建模方法、装置、介质及电子设备
Riffault et al. Quantifying EGS Permeability Development Using Induced Seismicity

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