CN111914609B - 井震联合叠前地质统计学弹性参数反演方法及装置 - Google Patents
井震联合叠前地质统计学弹性参数反演方法及装置 Download PDFInfo
- Publication number
- CN111914609B CN111914609B CN202010380490.8A CN202010380490A CN111914609B CN 111914609 B CN111914609 B CN 111914609B CN 202010380490 A CN202010380490 A CN 202010380490A CN 111914609 B CN111914609 B CN 111914609B
- Authority
- CN
- China
- Prior art keywords
- data
- seismic
- prestack
- grid node
- grid
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 68
- 239000011159 matrix material Substances 0.000 claims description 73
- 230000006870 function Effects 0.000 claims description 37
- 238000004364 calculation method Methods 0.000 claims description 36
- 238000004422 calculation algorithm Methods 0.000 claims description 15
- 238000004590 computer program Methods 0.000 claims description 13
- 238000006243 chemical reaction Methods 0.000 claims description 5
- 241000239366 Euphausiacea Species 0.000 claims 2
- 238000010586 diagram Methods 0.000 description 28
- 238000004088 simulation Methods 0.000 description 13
- 238000009826 distribution Methods 0.000 description 11
- 238000005070 sampling Methods 0.000 description 10
- 230000000875 corresponding effect Effects 0.000 description 7
- 239000013598 vector Substances 0.000 description 7
- 230000008569 process Effects 0.000 description 5
- 238000012545 processing Methods 0.000 description 4
- 238000011161 development Methods 0.000 description 3
- 230000000149 penetrating effect Effects 0.000 description 3
- 230000035939 shock Effects 0.000 description 3
- 238000003860 storage Methods 0.000 description 3
- 238000000342 Monte Carlo simulation Methods 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 2
- 238000012804 iterative process Methods 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 239000013049 sediment Substances 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000002596 correlated effect Effects 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 230000001351 cycling effect Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 238000005315 distribution function Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- 238000005293 physical law Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000004576 sand Substances 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000002948 stochastic simulation Methods 0.000 description 1
- 230000001502 supplementing effect Effects 0.000 description 1
- 238000003325 tomography Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V11/00—Prospecting or detecting by methods combining techniques covered by two or more of main groups G01V1/00 - G01V9/00
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/08—Feature extraction
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/30—Assessment of water resources
Abstract
本发明提供了一种井震联合叠前地质统计学弹性参数反演方法及装置，涉及地质勘探技术领域，该方法包括：获取测井数据、叠前地震数据、地震子波数据和地质空间结构模型数据；确定访问所有网格节点的路径；根据测井数据、叠前地震数据和地质空间结构模型数据，确定第一网格节点的测井数据点和地震数据点；根据观测数据、地震子波数据、变差函数和初始叠前弹性参数计算第一网格节点的克里金期望和克里金协方差；根据克里金期望和克里金协方差生成第一网格节点的井震联合叠前地质统计学弹性参数，按照路径生成除第一网格节点外剩余网格节点的井震联合叠前地质统计学弹性参数。本发明可以提高反演结果分辨率和效率。
Description
技术领域
本发明涉及地质勘探技术领域，尤其是涉及一种井震联合叠前地质统计学弹性参数反演方法及装置。
背景技术
叠前地震弹性参数反演作为储层描述的工具之一，在获取储层物性信息、评价储层质量方面有着重要应用。根据反问题中解的形式，可以将地震反演分为确定性反演和随机反演两大类。当获得了模型参数的后验概率分布之后，确定性反演是采用最优化理论，寻求一个最大后验概率解，使得观测地震数据与合成地震数据的误差达到最小，这类方法通常分辨率较低；随机反演则是使用采样算法，需经历大量迭代过程，计算效率较低。
发明内容
本发明提供了一种井震联合叠前地质统计学弹性参数反演方法及装置，可以提高叠前弹性参数反演结果分辨率，并能够提升计算效率。
第一方面，本发明实施例提供了一种井震联合叠前地质统计学弹性参数反演方法，该方法包括：获取测井数据、叠前地震数据、地震子波数据和地质空间结构模型数据；所述地质空间结构模型包括多个网格节点；所述地质空间结构模型数据包括变差函数和每个所述网格节点的初始叠前弹性参数；确定访问所有所述网格节点的路径；根据所述测井数据、所述叠前地震数据和所述地质空间结构模型数据，确定第一网格节点的测井数据点和地震数据点；根据所述测井数据点和所述地震数据点确定观测数据，并根据所述观测数据、所述地震子波数据、所述变差函数和所述初始叠前弹性参数计算所述第一网格节点的克里金期望和克里金协方差；根据所述克里金期望和所述克里金协方差生成所述第一网格节点的井震联合叠前地质统计学弹性参数，按照所述路径生成除所述第一网格节点外剩余网格节点的井震联合叠前地质统计学弹性参数。
第二方面，本发明实施例还提供一种井震联合叠前地质统计学弹性参数反演装置，该装置包括：获取模块，用于获取测井数据、叠前地震数据、地震子波数据和地质空间结构模型数据；所述地质空间结构模型包括多个网格节点；所述地质空间结构模型数据包括变差函数和每个所述网格节点的初始叠前弹性参数；路径模块，用于确定访问所有所述网格节点的路径；确定模块，用于根据所述测井数据、所述叠前地震数据和所述地质空间结构模型数据，确定第一网格节点的测井数据点和地震数据点；计算模块，用于根据所述测井数据点和所述地震数据点确定观测数据，并根据所述观测数据、所述地震子波数据、所述变差函数和所述初始叠前弹性参数计算所述第一网格节点的克里金期望和克里金协方差；生成模块，用于根据所述克里金期望和所述克里金协方差生成所述第一网格节点的井震联合叠前地质统计学弹性参数，按照所述路径生成除所述第一网格节点外剩余网格节点的井震联合叠前地质统计学弹性参数。
第三方面，本发明实施例还提供一种计算机设备，包括存储器、处理器，所述存储器中存储有可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现上述叠前弹性参数反演方法。
第四方面，本发明实施例还提供一种具有处理器可执行的非易失的程序代码的计算机可读介质，所述程序代码使所述处理器执行上述叠前弹性参数反演方法。
本发明实施例带来了以下有益效果：本发明实施例提供了一种井震联合叠前地质统计学弹性参数反演方案，该方案首先获取测井数据、叠前地震数据、地震子波数据和地质空间结构模型数据，其中，地质空间结构模型包括多个网格节点，每个网格节点对应一个初始叠前弹性参数，地质空间结构模型数据还包括变差函数，变差函数用于描述地质空间结构模型的空间结构，之后，确定访问所有网格节点的路径，根据测井数据、叠前地震数据和地质空间结构模型数据，确定路径中当前节点，即第一网格节点的测井数据点和地震数据点，根据测井数据点和地震数据点确定观测数据，再根据观测数据、地震子波数据、变差函数和初始叠前弹性参数计算第一网格节点的克里金期望和克里金协方差，从而再基于克里金期望和克里金协方差得到第一网格节点的井震联合叠前地质统计学弹性参数，按照确定好的路径生成除第一网格节点外剩余网格节点的井震联合叠前地质统计学弹性参数。本发明实施例可以提反演结果分辨率和效率。
本发明的其他特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。
为使本发明的上述目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附附图，作详细说明如下。
附图说明
为了更清楚地说明本发明具体实施方式或现有技术中的技术方案，下面将对具体实施方式或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施方式，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的井震联合叠前地质统计学弹性参数反演方法流程图；
图2中的(a)为本发明实施例提供的传统贝叶斯线性化AVO反演结果示意图；
图2中的(b)为本发明实施例提供的井震联合地质统计学反演结果示意图；
图3中的(a)为本发明实施例提供的二维理论模型与测井虚线示意图；
图3中的(b)为本发明实施例提供的小角度(0-10°)叠加地震数据示意图；
图4中的(a)为本发明实施例提供的测井数据第1次序贯高斯模拟实现结果示意图；
图4中的(b)为本发明实施例提供的测井数据第2次序贯高斯模拟实现结果示意图；
图5中的(a)为本发明实施例提供的地震数据最小二乘法反演结果示意图；
图5中的(b)为本发明实施例提供的井震联合反演第1次实现结果示意图；
图5中的(c)为本发明实施例提供的井震联合反演第2次实现结果示意图；
图5中的(d)为本发明实施例提供的井震联合反演第3次实现结果示意图；
图6中的(a)为本发明实施例提供的水平变差函数反演结果地质统计参数示意图；
图6中的(b)为本发明实施例提供的垂直变差函数反演结果地质统计参数示意图；
图7为本发明实施例提供的实际三维工区井震联合叠前地质统计学反演结果示意图；
图8为本发明实施例提供的高分辨率井震联合叠前地质统计学反演实施的装置图；
图9中的(a)为本发明实施例提供的一维理论模型曲线与硬数据点(测井数据)示意图；
图9中的(b)为本发明实施例提供的合成角道集地震数据示意图；
图10为本发明实施例提供的一种井震联合叠前地质统计学弹性参数反演装置结构示意框图；
图11为本发明实施例提供的另一种井震联合叠前地质统计学弹性参数反演装置结构示意框图；
图12为本发明实施例提供的计算机设备结构框图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合附图对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
目前，根据反问题中解的形式，可以将地震反演分为确定性反演和随机反演两大类。事实上，无论是确定性反演还是随机反演，都可以表达为贝叶斯理论框架下的概率推断问题。二者的区别在于，当获得了模型参数的后验概率分布之后，确定性反演是采用最优化理论，寻求一个最大后验概率解，使得观测地震数据与合成地震数据的误差达到最小，这类方法通常受地震资料频带范围的限制，无法识别地震分辨率之外的薄储层信息；随机反演则是使用采样算法，完整的探索后验解空间，每一次实现结果都代表可能的解，不仅能够获得高分辨率反演结果，并且能够对反演结果进行不确定性评价。
在随机地震反演中，根据采样策略的不同，又可将随机反演分为传统蒙特卡洛反演和地质统计学反演。蒙特卡洛方法能够较好的对高维数据进行降维处理，适用于大规模、多参数的反演问题，但是此类方法通常需要经过一系列的“拒绝-接受”重复迭代采样过程，计算效率较低。当地下介质的模型参数能用统计学特征(例如期望、协方差)表述时，地质统计学中的序贯模拟能够提供一种高效的随机采样方法，而且每次采样结果与观测数据和先验信息都相符合。1994年，Haas和Dubrule首次将地质统计学应用到地震反演中，将地震资料作为约束条件，控制测井资料序贯模拟结果的取舍，构成了传统地质统计学反演的基本思想。但是在这种方法中，地震资料仅起到了匹配约束的作用，没有直接参与反演计算，与蒙特卡洛方法类似，该方法仍需经历大量迭代过程，计算效率低。2003年，Buland和Omre提出了空间结构约束下的贝叶斯线性化AVO(Amplitude variation with offset，振幅随偏移距的变化)反演方法，假设纵、横波速度和密度服从对数高斯分布，利用测井资料求取垂向变差函数并且构建先验模型协方差矩阵，获得了后验期望和协方差的解析表达式，从而可以对后验分布进行快速采样实现。此方法中，地震资料是直接参与了反演计算，但是对于测井资料，仅利用其计算了模型的先验协方差矩阵，没有直接参与反演计算。2006年，Hansen等在模型参数服从高斯分布假设下，提出了线性反问题的序贯高斯模拟求解方法，将地质统计学与射线理论相结合，成功应用于井间探地雷达层析成像反演中。在这种方法中，探地雷达旅行时观测数据和测井资料都直接参与了反演计算，提高探地雷达反演的分辨率和降低了反演结果的不确定性。这三篇文章奠定了地质统计学反演的方法基础。近年来的发展主要在于模型参数先验分布假设的拓展和包括训练图像在内的高效复杂先验信息表征，使得地质统计学反演具有更加广泛的适用性。
在地震勘探领域，目前尚未见有将叠前地震数据、测井数据、地质空间结构先验信息等不同类型、不同尺度的数据整合，进行高分辨率井震联合叠前地质统计学反演的有效方法。
基于此，本发明实施例提供的一种井震联合叠前地质统计学弹性参数反演方法及装置，将地质统计学与地震反演理论相结合，整合不同类型、不同尺度的数据，能够提高叠前弹性参数反演结果分辨率，并且降低反演结果的不确定性。
为便于对本实施例进行理解，首先对本发明实施例所公开的一种井震联合叠前地质统计学弹性参数反演方法进行详细介绍。
本发明实施例提供了一种叠前弹性参数反演方法，参见图1所示的一种井震联合叠前地质统计学弹性参数反演方法流程图，该方法包括以下步骤：
步骤S102，获取测井数据、叠前地震数据、地震子波数据和地质空间结构模型数据。
在本发明实施例中，测井数据、叠前地震数据地、震子波数据和地质空间结构模型数据都是预先准备的数据。其中，地质空间结构模型可以用于描述地质空间结构信息，地质空间结构模型中包括多个网格节点，每个网格节点的叠前弹性参数包括纵波速度、横波速度以及密度三个维度的数据。本发明实施例中初始叠前弹性参数可以预先进行设定，例如初始叠前弹性参数可以是反射界面处上下两层介质的平均纵波速度数据、平均横波速度数据和平均密度数据。
地质空间结构模型数据还包括变差函数，变差函数可以用于描述地质空间结构模型的空间结构。
步骤S104，确定访问所有网格节点的路径。
在本发明实施例中，地质空间结构模型中包括多个网格节点，需要对每个网格节点进行访问，对每个网格节点的初始叠前弹性参数进行修正及优化，得到叠前弹性参数反演结果，因此，需要确定对各个网格节点进行访问的先后顺序，并且，确定出的路基需要覆盖所有网格节点。
步骤S106，根据测井数据、叠前地震数据和地质空间结构模型数据，确定第一网格节点的测井数据点和地震数据点。
在本发明实施例中，第一网格节点是多个网格节点中正在被访问的当前网格节点，在每个网格节点的周围都分布着多个网格节点，将测井数据按照地质空间结构模型数据进行网格划分后，第一网格节点附近的网格节点可以作为测井数据点；将叠前地震数据按照地质空间结构模型数据进行网格划分后，第一网格节点附近的网格节点可以作为地震数据点。地震数据点的个数通常大于测井数据点的个数。
步骤S108，根据测井数据点和地震数据点确定观测数据，并根据观测数据、地震子波数据、变差函数和初始叠前弹性参数计算第一网格节点的克里金期望和克里金协方差。
在本发明实施例中，根据测井数据点对应的测井数据和地震数据点对应的地震数据可以确定观测数据。从而，观测数据联合两类观测数据，之后，根据观测数据、地震子波数据、变差函数和初始叠前弹性参数再进行计算，能够有助于提升反演结果的分辨率。通过计算第一网络节点的克里金期望和克里金协方差可以得到第一网格节点对应的叠前弹性参数的数据分布情况。
步骤S110，根据克里金期望和克里金协方差生成第一网格节点的井震联合叠前地质统计学弹性参数，按照路径生成除第一网格节点外剩余网格节点的井震联合叠前地质统计学弹性参数。
在本发明实施例中，在得到了第一网格节点的克里金期望和克里金协方差之后，根据克里金期望和克里金协方差确定数据范围；在数据范围中随机选择目标数据作为第一网格节点的井震联合叠前地质统计学弹性参数。按照确定好的路径，对除第一网格节点外剩余的每个网格节点进行访问，按照步骤S106至步骤S110，进行计算，可以得到每个网格节点的井震联合叠前地质统计学弹性参数，进而，得到井震联合叠前地质统计学弹性参数反演结果。
本发明实施例提供了一种井震联合叠前地质统计学弹性参数反演方案，该方案首先获取测井数据、叠前地震数据、地震子波数据和地质空间结构模型数据，其中，地质空间结构模型包括多个网格节点，每个网格节点对应一个初始叠前弹性参数，地质空间结构模型数据还包括变差函数，变差函数用于描述地质空间结构模型的空间结构，之后，确定访问所有网格节点的路径，根据测井数据、叠前地震数据和地质空间结构模型数据，确定路径中当前节点，即第一网格节点的测井数据点和地震数据点，根据测井数据点和地震数据点确定观测数据，再根据观测数据、地震子波数据、变差函数和初始叠前弹性参数计算第一网格节点的克里金期望和克里金协方差，从而再基于克里金期望和克里金协方差得到第一网格节点的井震联合叠前地质统计学弹性参数，按照确定好的路径生成除第一网格节点外剩余网格节点的井震联合叠前地质统计学弹性参数。本发明实施例可以提高反演结果分辨率和效率。
考虑到为了提升网格节点反演结果的准确性，根据测井数据、叠前地震数据和地质空间结构模型数据，确定第一网格节点的测井数据点和地震数据点，可以按照如下步骤执行：
根据测井数据和地质空间结构模型数据，将第一网格节点的预设距离范围内的网格节点作为测井数据点；根据叠前地震数据和地质空间结构模型数据，将第一网格节点的预设距离范围内的网格节点作为地震数据点。
在本发明实施例中，预设距离范围可以根据实际需求进行设置，例如，可以以第一网格节点为中心，作一个球体，球体内的范围可以作为预设距离范围，球体的半径根据需求进行设置，也可以以第一网格节点为中心，作长方体或正方体等，本发明实施例对此不作具体限定。测井数据和叠前地震数据都是三维体数据，因此，根据地质空间结构模型数据中的多个网格节点，可以分别对测井数据和叠前地震数据进行划分，对于测井数据，第一网格节点的预设距离范围内的网格节点作为测井数据点，对于叠前地震数据，第一网格节点的预设距离范围内的网格节点作为地震数据点。
需要说明的是，可以选择预设距离范围内的所有点作为测井数据点或地震数据点，也可以选择预设距离范围内的部分点作为测井数据点或地震数据点，具体可以根据实际需求进行选择，本发明实施例对此不作具体限定。
考虑到联合井震数据进行反演，为进一步提升反演结果的分辨率，根据测井数据点和地震数据点确定观测数据，可以按照如下步骤执行：
根据测井数据点的测井数据确定测井观测数据；对地震数据点的地震数据进行时深转换，得到地震观测数据；对测井观测数据和地震观测数据进行组合，得到观测数据。
在本发明实施例中，在确定了测井数据点之后，根据测井数据点的位置可以确定测井数据，将测井数据直接作为测井数据点的测井观测数据。由于采集到的地震数据为时间域数据，在确定了地震数据点之后，先对地震数据点的地震数据进行时深转换，得到深度域数据，之后，将深度域的地震数据作为地震数据点的地震观测数据。分别将测井观测数据和地震观测数据记为矩阵的形式，将两个矩阵进行组合，即可得到观测数据，例如，测井观测数据为dobsA，地震观测数据为dobsB，合并可得观测数据为
考虑到为了提升反演的效率，根据观测数据、地震子波数据、变差函数和初始叠前弹性参数计算第一网格节点的克里金期望和克里金协方差，可以按照如下步骤执行：
根据初始叠前弹性参数和地震子波数据计算地震正演算子矩阵；根据初始叠前弹性参数和变差函数计算先验模型协方差矩阵；根据初始叠前弹性参数、观测数据、地震正演算子矩阵和先验模型协方差矩阵计算第一网格节点的克里金期望；根据地震正演算子矩阵和先验模型协方差矩阵计算克里金协方差。
在本发明实施例中，根据初始叠前弹性参数和地震子波数据得到地震正演算子矩阵，根据初始叠前弹性参数和变差函数得到先验模型协方差矩阵，进而根据初始叠前弹性参数、观测数据、地震正演算子矩阵和先验模型协方差矩阵得到第一网格节点的克里金期望，根据地震正演算子矩阵和先验模型协方差矩阵得到克里金协方差。通过逐步计算，实现对地震数据、测井数据和地质空间结构先验信息的整合。上述计算过程具体可以按照如下步骤执行。
根据初始叠前弹性参数和地震子波数据按照如下公式计算地震正演算子矩阵：其中，G为地震正演算子矩阵，W为地震子波数据，A为对角矩阵，D为差分矩阵。
在本发明实施例中，假设平面界面两侧为各向同性、弹性介质，且两侧介质的纵波速度vP、横波速度vS和密度ρ具有较小的差异，则当入射纵波穿过此界面时，反射纵波的时间连续反射系数rPP(t,θ)与入射角θ的关系可以表达为：
其中系数项和aρ(t,θ)由Aki和Richards公式给出：
并且和/>为反射界面处上下两层介质的平均纵、横波速度，/>为反射界面处入射角与透射角的平均，其中透射角根据Snell定律(斯涅尔定律)获得。在给定时间间隔和一系列入射角的情况下，可将式(1)改写为矩阵-向量的离散表达形式：
或进一步简化为
其中r(θi)为nm维列向量，由某一角度下各离散时间点的反射系数构成； 和Aρ(θi)为nm×nm维对角阵，对角线元素分别由式(2)至(4)中/> 和aρ(t,θ)计算得到；D为nm×nm维差分矩阵；m＝[lnvP,lnvS,lnρ]T为3nm维列向量，由模型参数的对数值组成；/>为Kronecker积(克罗内克积)算子；nθ为入射角个数；nm为模型参数vP、vS和ρ的维度。然后根据褶积模型，可将观测地震资料dobs表示为：
其中由不同地震角道集dobs(θi)构成，各角道集与该道的模型参数采样点相等，即各角道集为nm维列向量；W为子波褶积矩阵；e为地震误差项，与dobs维度相同；并且将/>定义为地震正演算子矩阵，维度为nmnθ×3nm。
根据上述公式(1)至(7)，可知，在给定了震子波数据和初始叠前弹性参数的情况下，根据初始叠前弹性参数，可以计算矩阵A，进而，可以根据公式计算地震正演算子矩阵G。
根据初始叠前弹性参数和变差函数按照如下公式计算先验模型协方差矩阵：Cm(p,q)＝C0ν(h)，其中，ν(h)为变差函数，h为模型空间中任意两点p和q间的距离，Cm为先验模型协方差矩阵，C0为时不变的协方差矩阵。
在本发明实施例中，贝叶斯反演方法能够综合观测数据中的多种先验信息，通过概率理论获得完备的后验解空间，并且能够表征反演结果的不确定性。假设模型参数m＝[lnvP,lnvS,lnρ]T服从多维高斯分布，即m服从下标3nm表示多维高斯分布的维度，则先验概率密度函数ρ(m)可以表示为
其中mprior＝[lnvP0,lnvS0,lnρ0]T为3nm维列向量，由纵、横波速度和密度的光滑背景模型组成；Cm为3nm×3nm维的先验模型协方差矩阵：
Cm(p,q)＝C0ν(h) (9)
其中h为模型空间中任意两点p和q间的距离，ν(h)为空间结构协方差矩阵，主要有高斯模型、球状模型、指数模型，以及多种模型的套和结构，C0为时不变的协方差矩阵，由对数域纵、横波速度和密度三个列向量计算得到，其表达式为
其中 且/>和/>为三参数间的相关系数。
根据上述公式(8)至(10)，可知，在给定了初始叠前弹性参数和变差函数的情况下，根据初始叠前弹性参数，可以计算时不变的协方差矩阵C0，进而，可以根据公式Cm(p,q)＝C0ν(h)，计算先验模型协方差矩阵。
根据初始叠前弹性参数、观测数据、地震正演算子矩阵和先验模型协方差矩阵按照如下公式计算第一网格节点的克里金期望：
其中，/>为第一网格节点的克里金期望，mpriorA+B为初始叠前弹性参数，CmA+B为第一网格节点的先验模型协方差矩阵，GA+B为第一网格节点的地震正演算子矩阵，CdA+B为地震协方差矩阵，dobsA+B为观测数据。
根据地震正演算子矩阵和先验模型协方差矩阵按照如下公式计算克里金协方差：其中，/>为第一网格节点的克里金协方差，CmA+B为第一网格节点的地震正演算子矩阵，GA+B为第一网格节点的地震正演算子矩阵，CdA+B为地震协方差矩阵。
在本发明实施例中，假设地震误差项服从零均值高斯分布，即e服从可得地震似然函数L(m)的表达式：
其中Cd为地震协方差矩阵，通过测井合成记录与实际井旁道地震资料估算得到，维度为nθnm×nθnm。根据贝叶斯理论，将式(8)和(11)相结合，可得全概率的后验解空间，其中后验期望为
并且后验协方差为
在地质统计学理论范畴中，式(12)和式(13)分别为在观测资料约束下的简单克里金估计，其中为克里金期望，/>为克里金方差。由此可以将地震反演与地质统计学理论相联系。
在地质统计学理论中，序贯高斯模拟是一种非迭代的高效采样算法，能够充分利用上述反问题中的线性和高斯分布特性，其基本思想为将已模拟点作为条件数据，同原始已知点一起，在给定邻域内参与下一个点的计算，直至随机路径中所有结点模拟完成。借鉴这种思想，在本发明实施例中，将观测数据划分为两类(A类和B类)，其中A类数据为模型参数的直接测量值，包括测井数据和先前模拟的点，也可称为“硬数据”或者“点数据”；B类数据为模型参数的线性平均观测值，即通过线性正演算子获得的测量数据，也可称为“软数据”或者“体数据”，这里指叠前地震角道集数据。联合两类观测数据，式(7)中的正演模型可改写为
或者简化为
dobsA+B＝GA+BmA+B+eA+B (15)
其中dobsA、GA、mA和eA分别代表A类数据的观测数据、正演算子、模型参数和误差项；dobsB、GB、mB和eB分别代表B类数据的观测数据、正演算子、模型参数和误差项；dobsA+B、GA+B、mA+B和eA+B代表相应的两类数据联合物理量。这里需要注意的是，与地震正演算子GB不同，正演算子GA为简单的单位矩阵，不依赖于任何物理规律。
根据式(12)和(13)的后验解形式，对式(15)进行贝叶斯反演，可以得到在测井数据、叠前地震数据和地质空间结构先验信息共同约束下的后验期望为
并且后验协方差为
其中mpriorA+B和dobsA+B均为(3nA+3nB)维列向量，nA和nB分别为邻域内A类和B类数据的采样点数，CmA+B和CdA+B分别为(3nA+nθnB)×(3nA+nθnB)维模型协方差矩阵和数据协方差矩阵，表达式为
其中CmAA和CmBB分别为A类和B类数据的模型协方差矩阵；CmAB为两类数据的互协方差矩阵；CdBB为B类数据的数据协方差矩阵。
式(16)和(17)中的后验期望和后验协方差即为多种数据共同约束下的局部克里金期望和方差，进而可以对后验概率模型进行序贯模拟求解，实现高分辨率地质统计学反演。若观测数据仅包含A类数据，即仅有测井数据，则该算法退化为地纯粹的质统计学序贯高斯模拟；相反，若观测数据仅包含B类数据，即只有地震数据，则该算法退化为传统的贝叶斯线性化AVO反演。因此，本发明中方法既适用于多井的成熟开发区块，同时也适用于无井或少井的初期勘探区块。
考虑到为了进一步提升反演结果的分辨率，根据克里金期望和克里金协方差生成第一网格节点的井震联合叠前地质统计学弹性参数之后，还可以执行如下步骤：
根据第一网格节点的井震联合叠前地质统计学弹性参数确定测井观测数据。
在本发明实施例中，在计算完成第一网格节点的井震联合叠前地质统计学弹性参数后，将该第一网格节点的井震联合叠前地质统计学弹性参数作为测井数据点的测井数据，在确定测井观测数据时使用。因此，在访问后续的网格节点时，可以得到优化后的测井观测数据，有利于提升反演结果的分辨率。
考虑到为了能够评估反演结果的不确定性，该方法可以执行如下步骤：
获取随机实现次数；按照实现次数生成多个叠前弹性参数。
在本发明实施例中，对于一次反演，假设地质空间结构模型包括k个网格节点，那么一次实现，会生成k个网格节点的叠前弹性参数。预先设置随机实现次数，例如，可以设置为n次，那么该方案会生成n次实现结果，每次实现结果中包括k个网格节点的叠前弹性参数。也就是说，对于某一个网格节点，会产生n个叠前弹性参数与之对应，可以根据n个叠前弹性参数的值评估反演结果的不确定性。
为了验证发明内容的有效性，分别设计一维和二维理论模型，以及三维实际资料的反演进行验证。图9为设计的一维理论模型，采用无条件序贯高斯模拟方法获得。在模型参数设置中，设定纵波速度、横波速度和密度的期望分别为3000m/s、1500m/s、2250kg/m3,方差分别0.0074、0.0074、0.0024，三参数间的相关性都设置为0.7，并且设置模型的空间结构为具有15ms变程的高斯变差函数。图9中的(a)中曲线代表真实理论模型，用于与反演结果的参照对比，实心点代表已知的点数据(A类数据)，相当于测井数据。图9中的(b)为使用Zoeppritz方程对理论模型曲线进行AVA正演模拟获得的地震角道集(B类数据)，添加了30％的相关噪声。以图9中的测井数据和地震角道集作为输入数据，分别进行传统贝叶斯线性化AVO反演和本发明井震联合叠前地质统计反演，并对反演结果进行100次随机模拟实现，结果分别如图2中的(a)和图2中的(b)所示。从图中箭头指示可以看出，由于传统贝叶斯线性化AVO反演方法仅利用了地震数据(B类数据)，在硬数据点(测井数据)位置处，反演结果与理论模型存在一定偏差，但是从本发明的反演结果可以看出，反演结果与理论模型在硬数据点处完全重合，以此可以验证测井数据参与了地震反演的计算。并且从随机模拟实现结果可以看出，本发明反演结果的不确定性低于传统贝叶斯线性化AVO反演方法。
图3为二维理论模型，以图3中的(a)中的测井数据和图3中的(b)中的合成地震数据作为输入，分别进行仅使用测井数据的序贯高斯模拟，仅使用地震数据的最小二乘反演，以及本发明的井震联合反演，结果分别如图4中的(a)-(b)、图5中的(a)和图5中的(b)-(d)所示。从图4中可以看出，仅使用测井数据的序贯高斯模拟方法虽然能够获得高分辨率插值结果，但是2次模拟结果都与图3中的(a)中的理论模型相差较大，并且不确定性较高，如图中箭头所示。从图5中的(a)可以看出，受地震资料频带范围的限制，最小二乘反演方法仅利用了地震数据，反演结果分辨率与地震数据相当，分辨率较低。然而，从图5中的(b)-(d)本发明所得到的反演结果看出，反演结果不仅具有高分辨率特征，而且3次实现结果较为一致，都与图3中的(a)中理论模型相似，反映了反演结果的不确定性较低。图6中的(a)和图6中的(b)为不同方法得到的水平和垂直变差函数，可以看出本发明(A+B类反演)的结果更加接近实际模型的地质空间结构信息。
图7为采用本发明方法获得的实际三维工区叠前反演结果，以检验该方法的横向效果。测试工区的面积约100km2，如图中箭头指示，反演结果能够清晰展示两条河道沉积的地质特征，符合研究区三角洲前缘水下分流河道砂体的沉积背景。
综上一维、二维理论模型和实际资料测试结果，可以验证本发明在提高叠前地震反演结果分辨率，降低反演结果多解性，和重建地质先验信息方面的优势。本发明方法既是对地震反演的发展，也是对地质统计学理论的补充，在薄层精细勘探和油气预测中具有较好的应用前景。
参见图8所示的高分辨率井震联合叠前地质统计学反演实施的装置图，在一个实施例中，可以按照如下步骤实施本方案：
步骤1：读取输入数据，包括测井数据和叠前角道集数据，模型参数空间结构信息(变差函数)ν(h)，地震子波W，初始光滑背景模型参数mpriorA+B，以及随机实现次数n；
步骤2：定义一条访问模型空间中所有网格节点的随机路径，节点数为k；
步骤3：在当前访问节点j的邻域内，选择测井数据和已模拟点数据作为A类观测数据dobsA，选择角道集地震数据dobsB，构建混合观测数据dobsB；
步骤4：求解混合观测数据约束下的克里金方程组(式16和17)，获得当前节点的局部克里金期望和协方差/>
步骤5：对局部高斯累积分布函数进行随机采样获得当前节点的随机实现值并将该模拟值添加到已知点数据中；
步骤6：重复步骤3-5，循环直至随机路径中所有节点计算完毕，获得井震联合地质统计学反演结果的1次实现；
步骤7：重复步骤2-6，循环直至达到设定的随机实现次数n，输出最终纵横波速度和密度模型参数的多次实现结果msim。
本发明实施例提供了一种井震联合叠前地质统计学弹性参数反演方法及装置，首先根据Aki和Richards近似公式建立地震AVA正演模型，然后在贝叶斯理论框架下将地震数据、测井数据和地质空间结构先验信息整合为地层参数的后验模型，最后采用序贯模拟方法对后验概率进行随机实现，获得多个既尊重地震数据又符合测井数据和先验地质信息的高分辨率反演结果，同时能够评估反演结果的不确定性。
本发明实施例还提供一种叠前弹性参数反演融合装置，参见图10所示的一种叠前弹性参数反演装置结构框图，该装置包括：
获取模块71，用于获取测井数据、叠前地震数据、地震子波数据和地质空间结构模型数据；地质空间结构模型包括多个网格节点；地质空间结构模型数据包括变差函数和每个网格节点的初始叠前弹性参数；路径模块72，用于确定访问所有网格节点的路径；确定模块73，用于根据测井数据、叠前地震数据和地质空间结构模型数据，确定第一网格节点的测井数据点和地震数据点；计算模块74，用于根据测井数据点和地震数据点确定观测数据，并根据观测数据、地震子波数据、变差函数和初始叠前弹性参数计算第一网格节点的克里金期望和克里金协方差；生成模块75，用于根据克里金期望和克里金协方差生成第一网格节点的井震联合叠前地质统计学弹性参数，按照路径生成除第一网格节点外剩余网格节点的井震联合叠前地质统计学弹性参数。
在一个实施例中，确定模块，具体用于：根据测井数据和地质空间结构模型数据，将第一网格节点的预设距离范围内的网格节点作为测井数据点；根据叠前地震数据和地质空间结构模型数据，将第一网格节点的预设距离范围内的网格节点作为地震数据点。
在一个实施例中，确定模块，具体用于：根据测井数据点的测井数据确定测井观测数据；对地震数据点的地震数据进行时深转换，得到地震观测数据；对测井观测数据和地震观测数据进行组合，得到观测数据。
在一个实施例中，计算模块，具体用于：根据初始叠前弹性参数和地震子波数据计算地震正演算子矩阵；根据初始叠前弹性参数和变差函数计算先验模型协方差矩阵；根据初始叠前弹性参数、观测数据、地震正演算子矩阵和先验模型协方差矩阵计算第一网格节点的克里金期望；根据地震正演算子矩阵和先验模型协方差矩阵计算克里金协方差。
在一个实施例中，计算模块，具体用于：根据初始叠前弹性参数和地震子波数据按照如下公式计算地震正演算子矩阵：其中，G为地震正演算子矩阵，W为地震子波数据，A为对角矩阵，D为差分矩阵。
在一个实施例中，计算模块，具体用于：根据初始叠前弹性参数和变差函数按照如下公式计算先验模型协方差矩阵：Cm(p,q)＝C0ν(h)，其中，ν(h)为变差函数，h为模型空间中任意两点p和q间的距离，Cm为先验模型协方差矩阵，C0为时不变的协方差矩阵。
在一个实施例中，计算模块，具体用于：根据初始叠前弹性参数、观测数据、地震正演算子矩阵和先验模型协方差矩阵按照如下公式计算第一网格节点的克里金期望：
其中，为第一网格节点的克里金期望，mpriorA+B为初始叠前弹性参数，CmA+B为第一网格节点的先验模型协方差矩阵，GA+B为第一网格节点的地震正演算子矩阵，CdA+B为地震协方差矩阵，dobsA+B为观测数据。
在一个实施例中，计算模块，具体用于：根据地震正演算子矩阵和先验模型协方差矩阵按照如下公式计算克里金协方差：
其中，/>为第一网格节点的克里金协方差，CmA+B为第一网格节点的地震正演算子矩阵，GA+B为第一网格节点的地震正演算子矩阵，CdA+B为地震协方差矩阵。
在一个实施例中，确定模块，还用于：根据第一网格节点的井震联合叠前地质统计学弹性参数确定测井观测数据。
在一个实施例中，参见图11所示的一种井震联合叠前地质统计学弹性参数反演装置结构框图，该装置还包括随机模块76，用于：获取随机实现次数；按照实现次数生成多个叠前弹性参数。
本发明实施例还提供一种计算机设备，参见图12所示的计算机设备结构示意框图，该计算机设备包括存储器81、处理器82，存储器中存储有可在处理器上运行的计算机程序，处理器执行计算机程序时实现上述任一种方法的步骤。
所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的计算机设备的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。
本发明实施例还提供一种具有处理器可执行的非易失的程序代码的计算机可读介质，程序代码使处理器执行上述任一种方法的步骤。
本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
最后应说明的是：以上所述实施例，仅为本发明的具体实施方式，用以说明本发明的技术方案，而非对其限制，本发明的保护范围并不局限于此，尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，其依然可以对前述实施例所记载的技术方案进行修改或可轻易想到变化，或者对其中部分技术特征进行等同替换；而这些修改、变化或者替换，并不使相应技术方案的本质脱离本发明实施例技术方案的精神和范围，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应所述以权利要求的保护范围为准。
Claims (22)
1.一种井震联合叠前地质统计学弹性参数反演方法，其特征在于，包括：
获取测井数据、叠前地震数据、地震子波数据和地质空间结构模型数据；所述地质空间结构模型包括多个网格节点；所述地质空间结构模型数据包括变差函数和每个所述网格节点的初始叠前弹性参数；
确定访问所有所述网格节点的路径；
根据所述测井数据、所述叠前地震数据和所述地质空间结构模型数据，确定第一网格节点的测井数据点和地震数据点；
根据所述测井数据点和所述地震数据点确定观测数据，并根据所述观测数据、所述地震子波数据、所述变差函数和所述初始叠前弹性参数计算所述第一网格节点的克里金期望和克里金协方差；
根据所述克里金期望和所述克里金协方差生成所述第一网格节点的井震联合叠前地质统计学弹性参数，按照所述路径生成除所述第一网格节点外剩余网格节点的井震联合叠前地质统计学弹性参数。
2.根据权利要求1所述的方法，其特征在于，根据所述测井数据、所述叠前地震数据和所述地质空间结构模型数据，确定第一网格节点的测井数据点和地震数据点，包括：
根据所述测井数据和所述地质空间结构模型数据，将所述第一网格节点的预设距离范围内的网格节点作为测井数据点；
根据所述叠前地震数据和所述地质空间结构模型数据，将所述第一网格节点的预设距离范围内的网格节点作为地震数据点。
3.根据权利要求1所述的方法，其特征在于，根据所述测井数据点和所述地震数据点确定观测数据，包括：
根据所述测井数据点的测井数据确定测井观测数据；
对所述地震数据点的地震数据进行时深转换，得到地震观测数据；
对所述测井观测数据和所述地震观测数据进行组合，得到观测数据。
4.根据权利要求1所述的方法，其特征在于，根据所述观测数据、所述地震子波数据、所述变差函数和所述初始叠前弹性参数计算所述第一网格节点的克里金期望和克里金协方差，包括：
根据所述初始叠前弹性参数和所述地震子波数据计算地震正演算子矩阵；
根据所述初始叠前弹性参数和所述变差函数计算先验模型协方差矩阵；
根据所述初始叠前弹性参数、所述观测数据、所述地震正演算子矩阵和所述先验模型协方差矩阵计算所述第一网格节点的克里金期望；
根据所述地震正演算子矩阵和所述先验模型协方差矩阵计算所述克里金协方差。
5.根据权利要求4所述的方法，其特征在于，包括：根据所述初始叠前弹性参数和所述地震子波数据按照如下公式计算地震正演算子矩阵：
其中，G为地震正演算子矩阵，W为地震子波数据，A为对角矩阵，D为差分矩阵。
6.根据权利要求4所述的方法，其特征在于，包括：根据所述初始叠前弹性参数和所述变差函数按照如下公式计算先验模型协方差矩阵：
Cm(p,q)＝C0ν(h)
其中，ν(h)为变差函数，h为模型空间中任意两点p和q间的距离，Cm为先验模型协方差矩阵，C0为时不变的协方差矩阵。
7.根据权利要求4所述的方法，其特征在于，包括：根据所述初始叠前弹性参数、所述观测数据、所述地震正演算子矩阵和所述先验模型协方差矩阵按照如下公式计算所述第一网格节点的克里金期望：
其中，为第一网格节点的克里金期望，mpriorA+B为初始叠前弹性参数，CmA+B为第一网格节点的先验模型协方差矩阵，GA+B为第一网格节点的地震正演算子矩阵，CdA+B为地震协方差矩阵，dobsA+B为观测数据。
8.根据权利要求4所述的方法，其特征在于，包括：根据所述地震正演算子矩阵和所述先验模型协方差矩阵按照如下公式计算所述克里金协方差：
其中，为第一网格节点的克里金协方差，CmA+B为第一网格节点的地震正演算子矩阵，GA+B为第一网格节点的地震正演算子矩阵，CdA+B为地震协方差矩阵。
9.根据权利要求1所述的方法，其特征在于，根据所述克里金期望和所述克里金协方差生成所述第一网格节点的井震联合叠前地质统计学弹性参数之后，还包括：
根据所述第一网格节点的井震联合叠前地质统计学弹性参数确定测井观测数据。
10.根据权利要求1-9任一项所述的方法，其特征在于，还包括：
获取随机实现次数；
按照所述实现次数生成多个叠前弹性参数。
11.一种井震联合叠前地质统计学弹性参数反演装置，其特征在于，包括：
获取模块，用于获取测井数据、叠前地震数据、地震子波数据和地质空间结构模型数据；所述地质空间结构模型包括多个网格节点；所述地质空间结构模型数据包括变差函数和每个所述网格节点的初始叠前弹性参数；
路径模块，用于确定访问所有所述网格节点的路径；
确定模块，用于根据所述测井数据、所述叠前地震数据和所述地质空间结构模型数据，确定第一网格节点的测井数据点和地震数据点；
计算模块，用于根据所述测井数据点和所述地震数据点确定观测数据，并根据所述观测数据、所述地震子波数据、所述变差函数和所述初始叠前弹性参数计算所述第一网格节点的克里金期望和克里金协方差；
生成模块，用于根据所述克里金期望和所述克里金协方差生成所述第一网格节点的井震联合叠前地质统计学弹性参数，按照所述路径生成除所述第一网格节点外剩余网格节点的井震联合叠前地质统计学弹性参数。
12.根据权利要求11所述的装置，其特征在于，所述确定模块，具体用于：
根据所述测井数据和所述地质空间结构模型数据，将所述第一网格节点的预设距离范围内的网格节点作为测井数据点；
根据所述叠前地震数据和所述地质空间结构模型数据，将所述第一网格节点的预设距离范围内的网格节点作为地震数据点。
13.根据权利要求11所述的装置，其特征在于，所述确定模块，具体用于：
根据所述测井数据点的测井数据确定测井观测数据；
对所述地震数据点的地震数据进行时深转换，得到地震观测数据；
对所述测井观测数据和所述地震观测数据进行组合，得到观测数据。
14.根据权利要求11所述的装置，其特征在于，所述计算模块，具体用于：
根据所述初始叠前弹性参数和所述地震子波数据计算地震正演算子矩阵；
根据所述初始叠前弹性参数和所述变差函数计算先验模型协方差矩阵；
根据所述初始叠前弹性参数、所述观测数据、所述地震正演算子矩阵和所述先验模型协方差矩阵计算所述第一网格节点的克里金期望；
根据所述地震正演算子矩阵和所述先验模型协方差矩阵计算所述克里金协方差。
15.根据权利要求14所述的装置，其特征在于，所述计算模块，具体用于：根据所述初始叠前弹性参数和所述地震子波数据按照如下公式计算地震正演算子矩阵：
其中，G为地震正演算子矩阵，W为地震子波数据，A为对角矩阵，D为差分矩阵。
16.根据权利要求14所述的装置，其特征在于，所述计算模块，具体用于：根据所述初始叠前弹性参数和所述变差函数按照如下公式计算先验模型协方差矩阵：
Cm(p,q)＝C0ν(h)
其中，ν(h)为变差函数，h为模型空间中任意两点p和q间的距离，Cm为先验模型协方差矩阵，C0为时不变的协方差矩阵。
17.根据权利要求14所述的装置，其特征在于，所述计算模块，具体用于：根据所述初始叠前弹性参数、所述观测数据、所述地震正演算子矩阵和所述先验模型协方差矩阵按照如下公式计算所述第一网格节点的克里金期望：
其中，为第一网格节点的克里金期望，mpriorA+B为初始叠前弹性参数，CmA+B为第一网格节点的先验模型协方差矩阵，GA+B为第一网格节点的地震正演算子矩阵，CdA+B为地震协方差矩阵，dobsA+B为观测数据。
18.根据权利要求14所述的装置，其特征在于，所述计算模块，具体用于：根据所述地震正演算子矩阵和所述先验模型协方差矩阵按照如下公式计算所述克里金协方差：
其中，为第一网格节点的克里金协方差，CmA+B为第一网格节点的地震正演算子矩阵，GA+B为第一网格节点的地震正演算子矩阵，CdA+B为地震协方差矩阵。
19.根据权利要求11所述的装置，其特征在于，所述确定模块，还用于：
根据所述第一网格节点的井震联合叠前地质统计学弹性参数确定测井观测数据。
20.根据权利要求11-19任一项所述的装置，其特征在于，还包括随机模块，用于：
获取随机实现次数；
按照所述实现次数生成多个叠前弹性参数。
21.一种计算机设备，包括存储器、处理器，所述存储器中存储有可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现上述权利要求1至10任一项所述的方法的步骤。
22.一种具有处理器可执行的非易失的程序代码的计算机可读介质，其特征在于，所述程序代码使所述处理器执行上述权利要求1至10任一项所述的方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010380490.8A CN111914609B (zh) | 2020-05-08 | 2020-05-08 | 井震联合叠前地质统计学弹性参数反演方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010380490.8A CN111914609B (zh) | 2020-05-08 | 2020-05-08 | 井震联合叠前地质统计学弹性参数反演方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111914609A CN111914609A (zh) | 2020-11-10 |
| CN111914609B true CN111914609B (zh) | 2023-12-01 |
Family
ID=73238007
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010380490.8A Active CN111914609B (zh) | 2020-05-08 | 2020-05-08 | 井震联合叠前地质统计学弹性参数反演方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111914609B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113608263B (zh) * | 2021-09-02 | 2023-06-20 | 中海石油(中国)有限公司深圳分公司 | 一种叠前反演的可行性评价方法、装置、设备和存储介质 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104614763A (zh) * | 2015-01-19 | 2015-05-13 | 中国石油大学(北京) | 基于反射率法的多波avo储层弹性参数反演方法及系统 |
| WO2017024702A1 (zh) * | 2015-08-11 | 2017-02-16 | 深圳朝伟达科技有限公司 | 一种射线弹性参数的反演系统 |
| WO2018107904A1 (zh) * | 2016-12-12 | 2018-06-21 | 中国石油大学(华东) | 一种精确反演杨氏模量和泊松比的方法 |
| CN110031896A (zh) * | 2019-04-08 | 2019-07-19 | 中国石油天然气集团有限公司 | 基于多点地质统计学先验信息的地震随机反演方法及装置 |
-
2020
- 2020-05-08 CN CN202010380490.8A patent/CN111914609B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104614763A (zh) * | 2015-01-19 | 2015-05-13 | 中国石油大学(北京) | 基于反射率法的多波avo储层弹性参数反演方法及系统 |
| WO2017024702A1 (zh) * | 2015-08-11 | 2017-02-16 | 深圳朝伟达科技有限公司 | 一种射线弹性参数的反演系统 |
| WO2018107904A1 (zh) * | 2016-12-12 | 2018-06-21 | 中国石油大学(华东) | 一种精确反演杨氏模量和泊松比的方法 |
| CN110031896A (zh) * | 2019-04-08 | 2019-07-19 | 中国石油天然气集团有限公司 | 基于多点地质统计学先验信息的地震随机反演方法及装置 |
Non-Patent Citations (2)
| Title |
|---|
| 叠前随机反演方法及其在薄层预测中的应用;孙月成;周家雄;马光克;郭爱华;王聪;;天然气工业(12);全文 * |
| 岩石物理参数高分辨率地质统计学反演;姜文龙;杨锴;;石油物探(06);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111914609A (zh) | 2020-11-10 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| KR101797451B1 (ko) | 상호상관 목적 함수를 통한 해양 스트리머 데이터에 대한 동시 소스 반전 | |
| CN110031896B (zh) | 基于多点地质统计学先验信息的地震随机反演方法及装置 | |
| CA2690373C (en) | Method for velocity analysis using waveform inversion in laplace domain for geophysical imaging | |
| CN104536044B (zh) | 一种地震数据的插值去噪方法以及系统 | |
| Lisitsa et al. | Combination of the discontinuous Galerkin method with finite differences for simulation of seismic wave propagation | |
| CN106526674B (zh) | 一种三维全波形反演能量加权梯度预处理方法 | |
| EP2715405B1 (en) | Method of processing seismic data by providing surface offset common image gathers | |
| US9625593B2 (en) | Seismic data processing | |
| CN108369289A (zh) | 使用全波场反演点扩展函数分析设计地球物理勘测的方法 | |
| CN103119552A (zh) | 同时源编码和源分离作为全波场反演的实际解决方案 | |
| CN102892972A (zh) | 地球物理数据的迭代反演方法中的伪迹减少 | |
| EA032186B1 (ru) | Сейсмическая адаптивная фокусировка | |
| US11733413B2 (en) | Method and system for super resolution least-squares reverse time migration | |
| CN103119472B (zh) | 利用同时和顺序源方法进行全波形反演的混合方法 | |
| CN109738950B (zh) | 基于三维稀疏聚焦域反演的噪声型数据一次波反演方法 | |
| CN109738952A (zh) | 基于全波形反演驱动的被动源直接偏移成像方法 | |
| CN111914609B (zh) | 井震联合叠前地质统计学弹性参数反演方法及装置 | |
| US11397273B2 (en) | Full waveform inversion in the midpoint-offset domain | |
| US20230140656A1 (en) | Method and system for determining seismic processing parameters using machine learning | |
| Copeland* et al. | Fracture surface extraction and stress field estimation from three-dimensional microseismic data | |
| Le Bouteiller | Eulerian approach of Hamilton-Jacobi equation with a discontinuous Galerkin method in heterogeneous anisotropic medium: Application to seismic imaging | |
| CN109061730B (zh) | 一种薄互层模型约束下的地震波阻抗确定性反演方法 | |
| CN113219534B (zh) | 一种叠前深度偏移速度质控方法、装置、介质及电子设备 | |
| US11860324B2 (en) | System and method for 3D seismic velocity estimation | |
| Kim et al. | Resolution enhancement for a seismic velocity model using machine learning |
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 |