CN110187384A - 贝叶斯时移地震差异反演方法及装置 - Google Patents
贝叶斯时移地震差异反演方法及装置 Download PDFInfo
- Publication number
- CN110187384A CN110187384A CN201910529992.XA CN201910529992A CN110187384A CN 110187384 A CN110187384 A CN 110187384A CN 201910529992 A CN201910529992 A CN 201910529992A CN 110187384 A CN110187384 A CN 110187384A
- Authority
- CN
- China
- Prior art keywords
- data
- variable quantity
- seismic
- elastic parameter
- inversion
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 106
- 239000011159 matrix material Substances 0.000 claims abstract description 86
- 238000009795 derivation Methods 0.000 claims abstract description 27
- 238000000605 extraction Methods 0.000 claims abstract description 16
- 230000006870 function Effects 0.000 claims description 64
- 238000012544 monitoring process Methods 0.000 claims description 42
- 238000009826 distribution Methods 0.000 claims description 24
- 238000004422 calculation algorithm Methods 0.000 claims description 16
- 238000005065 mining Methods 0.000 claims description 16
- 238000007619 statistical method Methods 0.000 claims description 14
- 238000004590 computer program Methods 0.000 claims description 10
- 239000000284 extract Substances 0.000 claims description 10
- 238000003860 storage Methods 0.000 claims description 9
- 206010003673 Atrioventricular block complete Diseases 0.000 claims description 8
- 238000011160 research Methods 0.000 claims description 8
- 201000002931 third-degree atrioventricular block Diseases 0.000 claims description 8
- 238000012360 testing method Methods 0.000 claims description 7
- 238000004458 analytical method Methods 0.000 claims description 6
- 241000208340 Araliaceae Species 0.000 claims description 5
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 claims description 5
- 235000003140 Panax quinquefolius Nutrition 0.000 claims description 5
- 230000008859 change Effects 0.000 claims description 5
- 235000008434 ginseng Nutrition 0.000 claims description 5
- 230000015572 biosynthetic process Effects 0.000 claims description 2
- 238000004364 calculation method Methods 0.000 claims description 2
- 230000008569 process Effects 0.000 description 11
- 238000007796 conventional method Methods 0.000 description 4
- 230000000903 blocking effect Effects 0.000 description 3
- 230000008878 coupling Effects 0.000 description 3
- 238000010168 coupling process Methods 0.000 description 3
- 238000005859 coupling reaction Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000004088 simulation Methods 0.000 description 3
- 238000010586 diagram Methods 0.000 description 2
- 238000005315 distribution function Methods 0.000 description 2
- 235000013399 edible fruits Nutrition 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000005484 gravity Effects 0.000 description 2
- 238000012545 processing Methods 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 239000004215 Carbon black (E152) Substances 0.000 description 1
- 241001269238 Data Species 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 230000006854 communication Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000013075 data extraction Methods 0.000 description 1
- 238000013499 data model Methods 0.000 description 1
- 230000001419 dependent effect Effects 0.000 description 1
- 238000013213 extrapolation Methods 0.000 description 1
- 239000004744 fabric Substances 0.000 description 1
- 238000001914 filtration Methods 0.000 description 1
- 238000010304 firing Methods 0.000 description 1
- 229930195733 hydrocarbon Natural products 0.000 description 1
- 150000002430 hydrocarbons Chemical class 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000002360 preparation method Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 230000006641 stabilisation Effects 0.000 description 1
- 238000011105 stabilization Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
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
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明提供了一种贝叶斯时移地震差异反演方法及装置。方法包括：对采集的地震数据进行匹配处理；利用匹配良好的地震数据提取角度依赖的地震子波，确定振幅缩放因子；提取得到多个弹性参数及对应的均值，计算包含多个弹性参数之间统计相关性的协方差矩阵；提取得到注采引起的弹性参数变化量曲线及对应的均值，计算包含多个弹性参数变化量统计相关性的协方差矩阵；根据沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型；估算正演算子矩阵；构建反演目标函数，对其求导并令导数为零，得到弹性参数变化量的求解表达式；对求解表达式进行迭代求解，获得最优弹性参数变化量反演结果。
Description
技术领域
本发明涉及油气田地震勘探开发和储层参数预测技术领域，尤指一种贝叶斯时移地震差异反演方法及装置。
背景技术
借助地球物理反演技术可以对油气开采引起的地震响应变化进行量化，得到储层弹性参数的变化量，进而结合地球物理方法进一步预测得到储层含油气饱和度、压力以及温度等物性参数的变化量，为开发中后期剩余油分布以及储层精细描述提供重要的地震信息。高精度反演方法使反演结果具有重要价值，可以提高地震数据解释的可信度。利用反演结果可以更好地实现储层空间描述，很大程度地提高储层表征的价值。叠前地震数据相对于叠后资料保留了更多的地下地质信息，基于叠前地震数据的反演能够获取更加丰富的弹性参数信息，因此叠前AVO/AVA时移地震反演方法被广泛应用于实际生产中。通常，这些方法都是以精确Zoeppritz方程近似公式为正演基础，然而近似公式推导过程中诸多的假设条件和较低的计算精度会极大地限制这些常规时移地震反演方法的适用性和反演精度，不能满足油藏开发人员对地震反演结果越来越高的精度要求。理论上，利用现有基于精确Zoeppritz方程的高精度反演方法开展时移地震分别反演可以有效克服基于近似公式的反演方法存在的一系列问题，但是分别反演不能有效利用相同位置不同时间采集的地震数据之间的耦合性，因此该思路并非最优时移地震反演策列。基于波动方程的时移地震反演算法虽然能够利用全波场的信息来预测储层弹性参数的变化，但其计算量大，在反演尺度和计算效率上不能满足实际油藏储层精细表征要求。
综上所述，目前基于叠前地震数据的时移地震反演方法研究存在以下问题：1、基于精确Zoeppritz方程近似公式的时移地震反演方法精度低、适用性差。2、基于精确Zoeppritz方程的时移地震分别反演方法不能有效利用基础数据和监测数据之间的耦合性。3、基于波动方程的时移地震反演方法计算量大，大型三维工区实际数据反演应用困难。4、因注采引起的弹性参数变化量具有明显的分块特征，需要有针对性的引入更加合理的先验模型。5、传统时移地震叠前AVA反演方法一般只考虑PP波数据，对其他多波信息的利用不充分。
发明内容
为了解决至少一个上述问题，本发明实施例提供一种贝叶斯时移地震差异反演方法，所述方法包括：
对同一工区不同时间采集的地震数据进行匹配处理，将匹配良好的基础数据和监测数据相减，得到时移地震差异数据；
利用匹配良好的地震数据提取角度依赖的地震子波，结合测井数据、精确Zoeppritz方程以及实际井旁道观测地震数据，确定振幅缩放因子；
根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个弹性参数之间统计相关性的协方差矩阵；
利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵；
利用地震构造解释资料和测井数据，根据沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型；
推导基于精确Zoeppritz方程的时移地震差异数据正演方程，并利用基于Zoeppritz方程的反演方法估算所述时移地震差异数据正演方程的正演算子矩阵；
根据Bayesian理论，同时引入包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式；
采用迭代重加权最小二乘算法对弹性参数变化量的求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果。
可选的，在本发明一实施例中，所述根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个参数之间统计相关性的协方差矩阵包括：通过测井数据统计分析，选取符合研究工区弹性参数分布特征的先验分布函数；结合统计分析得到包含各参数之间统计相关性的协方差矩阵构建用于基础数据对应的弹性参数反演的先验模型。
可选的，在本发明一实施例中，所述利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵包括：利用二次测井或完全非线性反演方法获取目标工区经过一段时间的油气开采之后对应井位的测井曲线；将相同位置不同时间获取的测井曲线相减，得到因油气开采引起的弹性参数变化量的真实曲线；引入服从微分拉普拉斯分布的垂向块约束项，利用曲线统计分析获得包含各参数变化量统计相关性的协方差矩阵，形成符合弹性参数变化量分布特征的包含垂向块约束项的先验模型。
可选的，在本发明一实施例中，所述根据Bayesian理论及包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式包括：根据Bayesian理论，引入包含垂向块约束项的先验模型，结合时移地震差异数据正演方程，构建最大后验概率意义下的反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式。
本发明实施例还提供一种贝叶斯时移地震差异反演装置，所述装置包括：
时移地震差异数据模块，用于对同一工区不同时间采集的地震数据进行匹配处理，将匹配良好的基础数据和监测数据相减，得到时移地震差异数据；
振幅缩放因子模块，用于利用匹配良好的地震数据提取角度依赖的地震子波，结合测井数据、精确Zoeppritz方程以及实际井旁道观测地震数据，确定振幅缩放因子；
第一协方差矩阵模块，用于根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个弹性参数之间统计相关性的协方差矩阵；
第二协方差矩阵模块，用于利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵；
初始模型模块，用于利用地震构造解释资料和测井数据，根据沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型；
正演算子矩阵模块，用于推导基于精确Zoeppritz方程的时移地震差异数据正演方程，并利用基于Zoeppritz方程的反演方法估算所述时移地震差异数据正演方程的正演算子矩阵；
求解表达式模块，用于根据Bayesian理论，同时引入包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式；
反演结果模块，用于采用迭代重加权最小二乘算法对弹性参数变化量的求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果。
可选的，在本发明一实施例中，所述第一协方差矩阵模块包括：先验分布函数单元，用于通过测井数据统计分析，选取符合研究工区弹性参数分布特征的先验分布函数；第一先验模型单元，用于结合统计分析得到包含各参数之间统计相关性的协方差矩阵构建用于基础数据对应的弹性参数反演的先验模型。
可选的，在本发明一实施例中，所述第二协方差矩阵模块包括：测井曲线单元，用于利用二次测井或完全非线性反演方法获取目标工区经过一段时间的油气开采之后对应井位的测井曲线；真实曲线单元，用于将相同位置不同时间获取的测井曲线相减，得到因油气开采引起的弹性参数变化量的真实曲线；第二先验模型单元，用于引入服从微分拉普拉斯分布的垂向块约束项，利用曲线统计分析获得包含各参数变化量统计相关性的协方差矩阵，形成符合弹性参数变化量分布特征的包含垂向块约束项的先验模型。
可选的，在本发明一实施例中，所述求解表达式模块包括：求解表达式单元，用于根据Bayesian理论，引入包含垂向块约束项的先验模型，结合时移地震差异数据正演方程，构建最大后验概率意义下的反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式。
本发明实施例还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现以下步骤：
对同一工区不同时间采集的地震数据进行匹配处理，将匹配良好的基础数据和监测数据相减，得到时移地震差异数据；
利用匹配良好的地震数据提取角度依赖的地震子波，结合测井数据、精确Zoeppritz方程以及实际井旁道观测地震数据，确定振幅缩放因子；
根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个弹性参数之间统计相关性的协方差矩阵；
利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵；
利用地震构造解释资料和测井数据，根据沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型；
推导基于精确Zoeppritz方程的时移地震差异数据正演方程，并利用基于Zoeppritz方程的反演方法估算所述时移地震差异数据正演方程的正演算子矩阵；
根据Bayesian理论，同时引入包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式；
采用迭代重加权最小二乘算法对弹性参数变化量的求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果。
本发明实施例还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现以下步骤：
对同一工区不同时间采集的地震数据进行匹配处理，将匹配良好的基础数据和监测数据相减，得到时移地震差异数据；
利用匹配良好的地震数据提取角度依赖的地震子波，结合测井数据、精确Zoeppritz方程以及实际井旁道观测地震数据，确定振幅缩放因子；
根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个弹性参数之间统计相关性的协方差矩阵；
利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵；
利用地震构造解释资料和测井数据，根据沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型；
推导基于精确Zoeppritz方程的时移地震差异数据正演方程，并利用基于Zoeppritz方程的反演方法估算所述时移地震差异数据正演方程的正演算子矩阵；
根据Bayesian理论，同时引入包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式；
采用迭代重加权最小二乘算法对弹性参数变化量的求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果。
本发明为差异反演提供了更高精度的正演算子，能够有效地克服基于精确Zoeppritz方程近似公式的反演方法在大偏移距信息利用以及方法适用性等方面存在的问题，进一步提升反演结果的稳定性和精度，相比于传统方法而言能够提供更加可靠的弹性参数变化量反演结果。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例一种贝叶斯时移地震差异反演方法的流程图；
图2A-图2C为本发明实施例中输入的叠前AVAPP波角道集地震数据，基础数据(图2A)、监测数据(图2B)及差异数据(图2C)；
图3A-图3C为本发明实施例中输入的叠前AVAPS波角道集地震数据，基础数据(图3A)、监测数据(图3B)及差异数据(图3C)；
图4A-图4C为本发明实施例输入的含噪叠前AVAPP波角道集地震数据，基础数据(图4A)、监测数据(图4B)及差异数据(图4C)；
图5A-图5C为本发明实施例中输入的含噪叠前AVAPS波角道集地震数据，基础数据(图5A)、监测数据(图5B)及差异数据(图5C)；
图6A-图6C为本发明实施例中无噪声情况下反演得到的纵波速度变化量(图6A)、横波速度变化量(图6B)及密度变化量(图6C)；
图7A-图7C为本发明实施例中信噪比为2:1的情况下反演得到的纵波速度变化量(图7A)、横波速度变化量(图7B)及密度变化量(图7C)；
图8为本发明实施例一种贝叶斯时移地震差异反演装置的结构示意图。
具体实施方式
本发明实施例提供一种贝叶斯时移地震差异反演方法及装置。
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
如图1所示为本发明实施例一种贝叶斯时移地震差异反演方法的流程图，图中所示方法包括：
步骤S1，对同一工区不同时间采集的地震数据进行匹配处理，将匹配良好的基础数据和监测数据相减，得到时移地震差异数据；
步骤S2，利用匹配良好的地震数据提取角度依赖的地震子波，结合测井数据、精确Zoeppritz方程以及实际井旁道观测地震数据，确定振幅缩放因子；
步骤S3，根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个弹性参数之间统计相关性的协方差矩阵；
步骤S4，利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵；
步骤S5，利用地震构造解释资料和测井数据，根据沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型；
步骤S6，推导精确Zoeppritz方程的时移地震差异数据正演方程，并利用基于Zoeppritz方程的反演方法估算所述时移地震差异数据正演方程的正演算子矩阵；
步骤S7，根据Bayesian理论，同时引入包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式；
步骤S8，采用迭代重加权最小二乘算法对弹性参数变化量的求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果。
作为本发明的一个实施例，根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个参数之间统计相关性的协方差矩阵包括：通过测井数据统计分析，选取符合研究工区弹性参数分布特征的先验分布函数；结合统计分析得到包含各参数之间统计相关性的协方差矩阵构建用于基础数据对应的弹性参数反演的先验模型。
作为本发明的一个实施例，利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵包括：利用二次测井或完全非线性反演方法获取目标工区经过一段时间的油气开采之后对应井位的测井曲线；将相同位置不同时间获取的测井曲线相减，得到因油气开采引起的弹性参数变化量的真实曲线；引入服从微分拉普拉斯分布的垂向块约束项，利用曲线统计分析获得包含各参数变化量统计相关性的协方差矩阵，形成符合弹性参数变化量分布特征的包含垂向块约束项的先验模型。
在本实施例中，根据Bayesian理论及包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式包括：根据Bayesian理论，引入包含垂向块约束项的先验模型，结合时移地震差异数据正演方程，构建最大后验概率意义下的反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式。
针对常规的时移地震反演方法所存在的精度低，适用性差等问题，本发明是在研究了以下问题的基础之上提出的：(1)基于精确Zoeppritz方程近似公式的时移地震反演方法精度低、适用性差。(2)基于精确Zoeppritz方程的时移地震分别反演方法不能有效利用基础数据和监测数据之间的耦合性。(3)基于波动方程的时移地震反演方法计算量大，大型三维工区实际数据反演应用困难。(4)因注采引起的弹性参数变化量具有明显的分块特征，需要有针对性的引入更加合理的先验模型。(5)传统时移地震叠前AVA反演方法一般只考虑PP波数据，对其他多波信息的利用不充分。
具体的，本发明的基于精确Zoeppritz方程的贝叶斯时移地震差异反演方法，具体包括：
(1)对同一工区不同时间采集的地震数据进行匹配处理并将匹配良好的基础数据和监测数据相减，得到时移地震差异数据：通过反复的时移地震匹配处理试验和效果分析，确定合理的时移地震匹配处理流程，然后基于该流程对同一工区不同时间采集的地震数据进行匹配处理，得到匹配良好的基础数据和监测数据，最后将二者相减得到时移地震差异数据。
(2)利用匹配良好的地震数据提取角度依赖的地震子波，并结合测井数据、精确Zoeppritz方程以及实际井旁道观测地震数据确定振幅缩放因子：基于实际的地震叠前道集和测井数据采取统计方法提取若干个依赖于入射角度的地震子波；以测井数据为输入模型利用精确Zoeppritz方程正演模拟角度域的PP波与PS波道集，与实际井旁角度域地震道集对比，计算振幅缩放因子，并应用于所提取的地震子波。
(3)基于工区内基础数据对应的测井数据提取得到纵波速度、横波速度和密度三个弹性参数及其均值，并计算包含三个参数之间统计相关性的协方差矩阵：通过测井数据统计分析，选取符合研究工区弹性参数分布特征的先验分布函数，结合统计得到包含各参数之间统计相关性的协方差矩阵构建用于基础数据对应弹性参数反演的先验模型。
(4)利用测井手段或是地球物理方法获取监测数据对应井位的实际测井曲线，然后基于工区内所有测井数据提取得到注采引起的弹性参数变化量曲线及其均值，并计算包含三个弹性参数变化量统计相关性的协方差矩阵：利用二次测井或完全非线性反演等方法获取目标工区经过一段时间的油气开采之后对应井位的测井曲线；将相同位置不同时间获取的测井曲线相减，得到因油气开采引起的弹性参数变化量的真实曲线；假设模型参数服从高斯分布的同时引入服从微分拉普拉斯分布的垂向块约束项，利用这些曲线统计分析获得包含各参数变化量统计相关性的协方差矩阵，最终形成符合弹性参数变化量分布特征的先验模型。
(5)利用地震构造解释资料和测井数据，基于沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型：利用地震构造解释资料，基于沉积模式建立地质模型，并将测井资料，按构造模式进行插值和外推，得到每条测线的初始参数模型，包括基础数据对应的弹性参数和差异数据对应的弹性参数变化量模型。
(6)推导基于精确Zoeppritz方程的时移地震差异数据正演方程，然后利用基于精确Zoeppritz方程的贝叶斯高精度反演方法估算该正演方程的正演算子矩阵：从精确Zoeppritz方程出发，借助泰勒级数展开推导得到时移地震差异数据对应的正演方程，然后结合步骤(3)和(5)得到的先验模型和初始模型，利用基于精确Zoeppritz方程的贝叶斯高精度反演方法估算该正演方程的正演算子矩阵。
(7)基于Bayesian理论，同时引入包含垂向块约束项的先验模型，构建最大后验概率意义下的反演目标函数，将该目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式：基于Bayesian反演理论框架，引入步骤(4)构建的包含垂向块约束项的先验模型，结合步骤(6)得到的时移地震差异数据正演方程，构建最大后验概率意义下的反演目标函数，将该目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式。
(8)采用迭代重加权最小二乘算法对弹性参数变化量的求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果：由于步骤(7)得到的弹性参数变化量求解表达式是非线性的，故采用迭代重加权最小二乘算法该求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果。
在本实施例中，图2A-图2C是本发明实施例输入的叠前AVA(Amplitude versusangle，振幅随入射角度的变化)PP波角道集地震数据，基础数据Base data(图2A)、监测数据Monitor data(图2B)、差异数据Difference data(图2C)；图3A-图3C是本发明实施例输入的叠前AVA(Amplitude versus angle，振幅随入射角度的变化)PS波角道集地震数据，基础数据Base data(图3A)、监测数据Monitor data(图3B)、差异数据Difference data(图3C)；图4A-图4C是本发明实施例输入的含噪叠前AVA(Amplitude versus angle，振幅随入射角度的变化)PP波角道集地震数据，基础数据Base data(图4A)、监测数据Monitor data(图4B)、差异数据Difference data(图4C)；图5A-图5C是本发明实施例输入的含噪叠前AVA(Amplitude versus angle，振幅随入射角度的变化)PS波角道集地震数据，基础数据Basedata(图5A)、监测数据Monitor data(图5B)、差异数据Difference data(图5C)；图6A-图6C是本发明实施例无噪声情况下采用基于精确Zoeppritz方程的贝叶斯时移地震差异反演方法反演得到的纵波速度变化量ΔVp(图6A)、横波速度变化量ΔVs(图6B)、密度变化量ΔRho(图6C)；图7A-图7C是本发明实施例信噪比为2:1的情况下采用基于精确Zoeppritz方程的贝叶斯时移地震差异反演方法反演得到的纵波速度变化量ΔVp(图7A)、横波速度变化量ΔVs(图7B)、密度变化量ΔRho(图7C)。
本发明具体技术方案与工作步骤详细叙述如下：
(1)本发明假设反演之前基础测线观测数据和监测测线观测数据具有良好的重复性，因此需要对基础数据和监测数据进行时移地震匹配处理。通常，需要通过反复的时移地震匹配处理试验和效果分析，确定最佳的时移地震匹配处理流程，其中主要包括面元重置、能量均衡处理、时差校正处理、频率均衡、相位校正、基于标志层的匹配滤波等。基于确定的最佳处理流程对同一工区不同时间采集的地震数据进行匹配处理，便可得到匹配良好的基础数据和监测数据，将二者相减可得到反演所需的时移地震差异数据。
(2)本发明假设反演之前地震子波已知，因此，需要基于实际的地震叠前道集和测井数据采取统计方法提取子波，子波在传播过程中受地层的影响会发生波形或频率变化，提取依赖于入射角度的地震子波能有效提高振幅匹配程度。
实际的地震振幅往往是相对值，采用精确Zoeppritz方程正演模拟的地震数据振幅与实际振幅存在一定数值差异。以测井数据为输入模型利用精确Zoeppritz反射系数方程正演模拟角度域的PP波与PS波道集，与实际井旁角度域地震道集对比，计算振幅缩放因子，并应用于所提取的地震子波，达到模拟记录与实际记录的振幅匹配。当地震数据信噪比较高时，为角道集的每一道使用统一的振幅缩放因子，以保证振幅随偏移距的变化关系；当信噪比低时，可分近、中、远偏移距分别计算振幅缩放因子，保证模拟记录与实际记录的最佳匹配，减少噪声对反演过程的影响。
(3)建立模型参数的先验模型。为了降低地震反演的不确定性，提高反演过程的稳定性，需要从其它途径获得地下介质地震弹性参数模型的信息作为先验信息，即建立用于反演的先验模型。本发明需要建立两个先验模型：一个用于基础数据弹性参数反演的先验模型，另一个用于时移地震差异数据反演的先验模型。用于基础数据弹性参数反演的先验模型构建的技术流程是首先对基础数据对应的测井数据进行统计分析，得到包含纵横波速度和密度三个弹性参数之间统计相关性的协方差矩阵，然后假设弹性参数服从特定的分布，并由此构建得到对应的先验模型。用于时移地震差异数据反演的先验模型构建需要已知监测数据对应井位的测井数据，因此需要进行测井数据获取。对于具备二次测井条件的井位采用测井的方法采集得到开发后的纵横波速度和密度曲线，对于不具备二次测井条件的井位采用基于完全非线性算法反演方法进行测井曲线重构。将得到的监测数据(开发后)对应的测井数据与已有基础数据(开发前)对应的测井数据相减，得到因油气开采引起的纵横波速度和密度三个弹性参数的真实变化曲线。基于这些变化曲线统计得到的各模型参数变化量及其均值，求取各参数变化量的自相关系数和互相关系数，构建协方差矩阵。
本发明涉及的时移地震差异反演采用高斯分布函数和微分拉普拉斯分布函数作为先验分布函数，结合已构建的协方差矩阵，形成符合该工区的先验分布函数。在后续的反演目标函数中先验模型相应的规则化表达式为：
F(Δm)＝const+F1(Δm)+F2(Δm) (1)
其中F1(Δm)代表高斯约束项，其表达式如下：
F2(Δm)代表块微分拉普拉斯约束项，其表达式如下：
其中，CΔm为包含三个差异数据相关性的协方差矩阵，N为模型参数的长度，μ为差异数据的均值向量(三个弹性参数需要分别求取)，D是一阶微分算子，kl,l＝1,2,3是尺度因子，三个弹性参数之间可能会不一样。
(4)利用地震构造解释资料和测井数据，基于沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型。
建立参数模型主要利用三维空间插值方法，其的技术流程是首先利用散点插值的方法对各个层位的数据进行插值，完成地质层位建模，然后根据地质层位进行弹性参数横向插值，即将测井信息进行横向插值，计算得到地下每个点上的参数值，完成初始参数建模的任务。
(5)推导基于精确Zoeppritz方程的时移地震差异数据正演方程。时移地震基础数据和监测数据的正演过程分别如下：
d1＝G(m1)+n1 (4)
d2＝G(m2)+n2 (5)
其中，G为Zoeppritz方程正演算子。d1为基础数据，d2为监测数据。m1为基础数据对应的弹性参数，m2为监测数据对应的弹性参数，n1和n2分别为两次数据采集时的噪声数据。
将等式(5)在基础数据对应的弹性参数m1处进行泰勒展开：
将等式(6)和等式(4)相减得：
令e＝n0+n2-n1，则有：
Δd＝LΔm+e (8)
其中，Δd是差异地震数据，Δm是弹性参数的变化量，是基于精确Zoeppritz方程的正演算子矩阵关于纵横波速度以及密度三个弹性参数的一阶偏导数在基础数据对应的弹性参数m1处的取值，e是与差异数据直接相关的误差项，也可以理解为噪声。
由于公式(8)所示正演方程中的正演算子是基于精确Zoeppritz方程的正演算子矩阵关于纵横波速度以及密度三个弹性参数的一阶偏导数在基础数据对应的弹性参数m1处的取值，因此在开展后续时移地震差异反演之前需要借助Zhou等(2017)提出的基于精确Zoeppritz方程的高精度贝叶斯非线性反演方法对基础数据进行反演，得到所需的正演算子矩阵L。
(6)基于Bayesian反演理论框架，引入已构建的包含垂向块约束项的先验模型，结时移地震差异数据正演方程，构建最大后验概率意义下的反演目标函数，将该目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式。
假设公式(8)中的误差项e服从零均值高斯分布，则似然函数可表示为：
上式中，CD是噪声的协方差矩阵，Nd是观测数据的长度。对上式取负对数：
将公式(1)和公式(10)代入贝叶斯反演理论框架，约去常数项，可以得到如下所示目标函数：
在实际处理过程中，通常假设观测数据的噪声是不相关的，则噪声的协方差矩阵会简化为一个对角矩阵，即CD＝σn 2I，其中σn 2表示噪声的方差，I是Nd×Nd的单位矩阵，Nd是观测数据的长度。对于多波地震数据，可以将上述目标函数进行拓展，本发明以PP-PS波联合反演为例进行了阐述。假设PP波差异地震数据的噪声方差为σΔPP，PS波差异地震数据的噪声方差为σΔPS，那么上述目标函数可以表示为：
其中α＝σΔpp/σΔps控制PS波数据的比重，β＝σΔpp控制先验信息的比重。
将上述目标函数关于Δm求导并令导数为零，
整理得，
其中
(7)公式(14)所示的弹性参数变化量求解表达式右端项元素A也包含Δm，因此需要采用迭代重加权最小二乘算法该求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果，如图6A-图6C所示，是本发明实施例基于精确Zoeppritz方程的贝叶斯时移地震差异反演方法反演得到的纵波速度变化量(图6A)、横波速度变化量(图6B)和密度变化量(图6C)，图中纵轴表示时间，单位秒，横轴自左至右表示纵波速度(单位：km/s)、横波速度(单位：km/s)和密度(单位：g/cm3))。基于精确Zoeppritz方程的贝叶斯时移地震差异反演方法能够较高精度的预测到油气开采引起的弹性参数变化量，由于引入了包含密度信息的先验分布，并采用PS横波进行联合反演，其对密度模型也预测准确。如图7A-图7C所示为加入信噪比为2的随机噪声时的反演结果，先验模型的引入对保持反演过程稳定起到了关键作用。
本发明具有以下优点：1、本发明为差异反演提供了更高精度的正演算子。2、本发明所述的时移地震差异反演方法为基于精确Zoeppritz方程的叠前AVA反演方法，能够有效地克服基于精确Zoeppritz方程近似公式的反演方法在大偏移距信息利用以及方法适用性等方面存在的问题。3、本发明联合多波数据进行时移地震差异反演，能够进一步提升反演结果的稳定性和精度。4、本发明引入了包含垂向微分拉普拉斯块约束项的先验模型，能够很好地刻画弹性参数变化量的分块特征，提升反演结果的精度。5、本发明能够提供一种高可靠性的时移地震差异反演方法，相比于传统方法而言能够提供更加可靠的弹性参数变化量反演结果。
如图8所示为本发明实施例一种贝叶斯时移地震差异反演装置的结构示意图，图中所示装置包括：
时移地震差异数据模块10，用于对同一工区不同时间采集的地震数据进行匹配处理，将匹配良好的基础数据和监测数据相减，得到时移地震差异数据；
振幅缩放因子模块20，用于利用匹配良好的地震数据提取角度依赖的地震子波，结合测井数据、精确Zoeppritz方程以及实际井旁道观测地震数据，确定振幅缩放因子；
第一协方差矩阵模块30，用于根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个弹性参数之间统计相关性的协方差矩阵；
第二协方差矩阵模块40，用于利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵；
初始模型模块50，用于利用地震构造解释资料和测井数据，根据沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型；
正演算子矩阵模块60，用于推导基于精确Zoeppritz方程的时移地震差异数据正演方程，并利用基于Zoeppritz方程的反演方法估算所述时移地震差异数据正演方程的正演算子矩阵；
求解表达式模块70，用于根据Bayesian理论，同时引入包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式；
反演结果模块80，用于采用迭代重加权最小二乘算法对弹性参数变化量的求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果。
作为本发明的一个实施例，第一协方差矩阵模块包括：先验分布函数单元，用于通过测井数据统计分析，选取符合研究工区弹性参数分布特征的先验分布函数；第一先验模型单元，用于结合统计分析得到包含各参数之间统计相关性的协方差矩阵构建用于基础数据对应的弹性参数反演的先验模型。
作为本发明的一个实施例，第二协方差矩阵模块包括：测井曲线单元，用于利用二次测井或完全非线性反演方法获取目标工区经过一段时间的油气开采之后对应井位的测井曲线；真实曲线单元，用于将相同位置不同时间获取的测井曲线相减，得到因油气开采引起的弹性参数变化量的真实曲线；第二先验模型单元，用于引入服从微分拉普拉斯分布的垂向块约束项，利用曲线统计分析获得包含各参数变化量统计相关性的协方差矩阵，形成符合弹性参数变化量分布特征的包含垂向块约束项的先验模型。
在本实施例中，求解表达式模块包括：求解表达式单元，用于根据Bayesian理论，引入包含垂向块约束项的先验模型，结合时移地震差异数据正演方程，构建最大后验概率意义下的反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式。
基于与上述一种贝叶斯时移地震差异反演方法相同的申请构思，本发明还提供了上述一种贝叶斯时移地震差异反演装置。由于该一种贝叶斯时移地震差异反演装置解决问题的原理与一种贝叶斯时移地震差异反演方法相似，因此该一种贝叶斯时移地震差异反演装置的实施可以参见一种贝叶斯时移地震差异反演方法的实施，重复之处不再赘述。
本发明为差异反演提供了更高精度的正演算子，能够有效地克服基于精确Zoeppritz方程近似公式的反演方法在大偏移距信息利用以及方法适用性等方面存在的问题，进一步提升反演结果的稳定性和精度，相比于传统方法而言能够提供更加可靠的弹性参数变化量反演结果。
本发明实施例还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现以下步骤：
对同一工区不同时间采集的地震数据进行匹配处理，将匹配良好的基础数据和监测数据相减，得到时移地震差异数据；
利用匹配良好的地震数据提取角度依赖的地震子波，结合测井数据、精确Zoeppritz方程以及实际井旁道观测地震数据，确定振幅缩放因子；
根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个弹性参数之间统计相关性的协方差矩阵；
利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵；
利用地震构造解释资料和测井数据，根据沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型；
推导基于精确Zoeppritz方程的时移地震差异数据正演方程，并利用基于Zoeppritz方程的反演方法估算所述时移地震差异数据正演方程的正演算子矩阵；
根据Bayesian理论，同时引入包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式；
采用迭代重加权最小二乘算法对弹性参数变化量的求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果。
本发明实施例还提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现以下步骤：
对同一工区不同时间采集的地震数据进行匹配处理，将匹配良好的基础数据和监测数据相减，得到时移地震差异数据；
利用匹配良好的地震数据提取角度依赖的地震子波，结合测井数据、精确Zoeppritz方程以及实际井旁道观测地震数据，确定振幅缩放因子；
根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个弹性参数之间统计相关性的协方差矩阵；
利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵；
利用地震构造解释资料和测井数据，根据沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型；
推导基于精确Zoeppritz方程的时移地震差异数据正演方程，并利用基于Zoeppritz方程的反演方法估算所述时移地震差异数据正演方程的正演算子矩阵；
根据Bayesian理论，同时引入包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式；
采用迭代重加权最小二乘算法对弹性参数变化量的求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果。
基于与上述一种贝叶斯时移地震差异反演方法相同的申请构思，本发明还提供了上述一种计算机设备及一种计算机可读存储介质。由于该一种计算机设备及一种计算机可读存储介质解决问题的原理与一种贝叶斯时移地震差异反演方法相似，因此该一种计算机设备及一种计算机可读存储介质的实施可以参见一种贝叶斯时移地震差异反演方法的实施，重复之处不再赘述。
本发明为差异反演提供了更高精度的正演算子，能够有效地克服基于精确Zoeppritz方程近似公式的反演方法在大偏移距信息利用以及方法适用性等方面存在的问题，进一步提升反演结果的稳定性和精度，相比于传统方法而言能够提供更加可靠的弹性参数变化量反演结果。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分步骤可以通过程序来指令相关的硬件来完成，该程序可以存储于一计算机可读取存储介质中，比如ROM/RAM、磁碟、光盘等。
以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (10)
1.一种贝叶斯时移地震差异反演方法，其特征在于，所述方法包括：
对同一工区不同时间采集的地震数据进行匹配处理，将匹配良好的基础数据和监测数据相减，得到时移地震差异数据；
利用匹配良好的地震数据提取角度依赖的地震子波，结合测井数据、精确Zoeppritz方程以及实际井旁道观测地震数据，确定振幅缩放因子；
根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个弹性参数之间统计相关性的协方差矩阵；
利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵；
利用地震构造解释资料和测井数据，根据沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型；
推导基于精确Zoeppritz方程的时移地震差异数据正演方程，并利用基于Zoeppritz方程的反演方法估算所述时移地震差异数据正演方程的正演算子矩阵；
根据Bayesian理论，同时引入包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式；
采用迭代重加权最小二乘算法对弹性参数变化量的求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果。
2.根据权利要求1所述的方法，其特征在于，所述根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个参数之间统计相关性的协方差矩阵包括：
通过测井数据统计分析，选取符合研究工区弹性参数分布特征的先验分布函数；
结合统计分析得到包含各参数之间统计相关性的协方差矩阵构建用于基础数据对应的弹性参数反演的先验模型。
3.根据权利要求1所述的方法，其特征在于，所述利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵包括：
利用二次测井或完全非线性反演方法获取目标工区经过一段时间的油气开采之后对应井位的测井曲线；
将相同位置不同时间获取的测井曲线相减，得到因油气开采引起的弹性参数变化量的真实曲线；
引入服从微分拉普拉斯分布的垂向块约束项，利用曲线统计分析获得包含各参数变化量统计相关性的协方差矩阵，形成符合弹性参数变化量分布特征的包含垂向块约束项的先验模型。
4.根据权利要求3所述的方法，其特征在于，所述根据Bayesian理论及包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式包括：根据Bayesian理论，引入包含垂向块约束项的先验模型，结合时移地震差异数据正演方程，构建最大后验概率意义下的反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式。
5.一种贝叶斯时移地震差异反演装置，其特征在于，所述装置包括：
时移地震差异数据模块，用于对同一工区不同时间采集的地震数据进行匹配处理，将匹配良好的基础数据和监测数据相减，得到时移地震差异数据；
振幅缩放因子模块，用于利用匹配良好的地震数据提取角度依赖的地震子波，结合测井数据、精确Zoeppritz方程以及实际井旁道观测地震数据，确定振幅缩放因子；
第一协方差矩阵模块，用于根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个弹性参数之间统计相关性的协方差矩阵；
第二协方差矩阵模块，用于利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵；
初始模型模块，用于利用地震构造解释资料和测井数据，根据沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型；
正演算子矩阵模块，用于推导基于精确Zoeppritz方程的时移地震差异数据正演方程，并利用基于Zoeppritz方程的反演方法估算所述时移地震差异数据正演方程的正演算子矩阵；
求解表达式模块，用于根据Bayesian理论，同时引入包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式；
反演结果模块，用于采用迭代重加权最小二乘算法对弹性参数变化量的求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果。
6.根据权利要求5所述的装置，其特征在于，所述第一协方差矩阵模块包括：
先验分布函数单元，用于通过测井数据统计分析，选取符合研究工区弹性参数分布特征的先验分布函数；
第一先验模型单元，用于结合统计分析得到包含各参数之间统计相关性的协方差矩阵构建用于基础数据对应的弹性参数反演的先验模型。
7.根据权利要求5所述的装置，其特征在于，所述第二协方差矩阵模块包括：
测井曲线单元，用于利用二次测井或完全非线性反演方法获取目标工区经过一段时间的油气开采之后对应井位的测井曲线；
真实曲线单元，用于将相同位置不同时间获取的测井曲线相减，得到因油气开采引起的弹性参数变化量的真实曲线；
第二先验模型单元，用于引入服从微分拉普拉斯分布的垂向块约束项，利用曲线统计分析获得包含各参数变化量统计相关性的协方差矩阵，形成符合弹性参数变化量分布特征的包含垂向块约束项的先验模型。
8.根据权利要求7所述的装置，其特征在于，所述求解表达式模块包括：求解表达式单元，用于根据Bayesian理论，引入包含垂向块约束项的先验模型，结合时移地震差异数据正演方程，构建最大后验概率意义下的反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式。
9.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现以下步骤：
对同一工区不同时间采集的地震数据进行匹配处理，将匹配良好的基础数据和监测数据相减，得到时移地震差异数据；
利用匹配良好的地震数据提取角度依赖的地震子波，结合测井数据、精确Zoeppritz方程以及实际井旁道观测地震数据，确定振幅缩放因子；
根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个弹性参数之间统计相关性的协方差矩阵；
利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵；
利用地震构造解释资料和测井数据，根据沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型；
推导基于精确Zoeppritz方程的时移地震差异数据正演方程，并利用基于Zoeppritz方程的反演方法估算所述时移地震差异数据正演方程的正演算子矩阵；
根据Bayesian理论，同时引入包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式；
采用迭代重加权最小二乘算法对弹性参数变化量的求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果。
10.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现以下步骤：
对同一工区不同时间采集的地震数据进行匹配处理，将匹配良好的基础数据和监测数据相减，得到时移地震差异数据；
利用匹配良好的地震数据提取角度依赖的地震子波，结合测井数据、精确Zoeppritz方程以及实际井旁道观测地震数据，确定振幅缩放因子；
根据工区内基础数据对应的测井数据，提取得到多个弹性参数及对应的均值，并计算包含多个弹性参数之间统计相关性的协方差矩阵；
利用测井手段或地球物理方法获取监测数据对应井位的实际测井曲线，根据工区内的测井数据提取得到注采引起的弹性参数变化量曲线及对应的均值，并计算包含多个弹性参数变化量统计相关性的协方差矩阵；
利用地震构造解释资料和测井数据，根据沉积模式建立时间域的基础数据对应的初始弹性参数模型和差异地震数据对应的初始弹性参数变化量模型；
推导基于精确Zoeppritz方程的时移地震差异数据正演方程，并利用基于Zoeppritz方程的反演方法估算所述时移地震差异数据正演方程的正演算子矩阵；
根据Bayesian理论，同时引入包含垂向块约束项的先验模型，构建反演目标函数，对所述反演目标函数关于弹性参数变化量求导并令导数为零，得到弹性参数变化量的求解表达式；
采用迭代重加权最小二乘算法对弹性参数变化量的求解表达式进行迭代求解，并通过反演残差控制最大迭代次数，获得最优弹性参数变化量反演结果。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910529992.XA CN110187384B (zh) | 2019-06-19 | 2019-06-19 | 贝叶斯时移地震差异反演方法及装置 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201910529992.XA CN110187384B (zh) | 2019-06-19 | 2019-06-19 | 贝叶斯时移地震差异反演方法及装置 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110187384A true CN110187384A (zh) | 2019-08-30 |
| CN110187384B CN110187384B (zh) | 2021-07-20 |
Family
ID=67722334
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201910529992.XA Active CN110187384B (zh) | 2019-06-19 | 2019-06-19 | 贝叶斯时移地震差异反演方法及装置 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110187384B (zh) |
Cited By (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110471106A (zh) * | 2019-09-20 | 2019-11-19 | 西南石油大学 | 一种基于滤波器设计的时移地震反演方法 |
| CN110661548A (zh) * | 2019-09-13 | 2020-01-07 | 广东石油化工学院 | 一种利用l1模反演的plc信号滤波方法和系统 |
| CN111736208A (zh) * | 2020-06-24 | 2020-10-02 | 重庆大学 | 变权重联合P波和S波初至数据的微震事件Bayes定位方法、系统及介质 |
| CN111948712A (zh) * | 2020-08-10 | 2020-11-17 | 中海石油(中国)有限公司 | 一种基于深度域地震记录的叠前线性反演方法 |
| CN115951405A (zh) * | 2023-03-08 | 2023-04-11 | 中国地震局地球物理研究所 | 一种基于稀疏化卷积和自适应矩估计的速度变化成像方法 |
Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104597490A (zh) * | 2015-01-28 | 2015-05-06 | 中国石油大学(北京) | 基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法 |
| US20160116620A1 (en) * | 2014-10-24 | 2016-04-28 | Ion Geophysical Corporation | Methods and systems for seismic inversion and related seismic data processing |
| CN105954804A (zh) * | 2016-07-15 | 2016-09-21 | 中国石油大学(北京) | 页岩气储层脆性地震预测方法及装置 |
| CN107092029A (zh) * | 2017-04-26 | 2017-08-25 | 中国石油大学(北京) | 一种地震反演方法和装置 |
| WO2018102043A1 (en) * | 2016-12-02 | 2018-06-07 | Exxonmobil Upstream Research Company | Method for estimating petrophysical properties for single or multiple scenarios from several spectrally variable seismic and full wavefield inversion products |
-
2019
- 2019-06-19 CN CN201910529992.XA patent/CN110187384B/zh active Active
Patent Citations (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20160116620A1 (en) * | 2014-10-24 | 2016-04-28 | Ion Geophysical Corporation | Methods and systems for seismic inversion and related seismic data processing |
| CN104597490A (zh) * | 2015-01-28 | 2015-05-06 | 中国石油大学(北京) | 基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法 |
| CN105954804A (zh) * | 2016-07-15 | 2016-09-21 | 中国石油大学(北京) | 页岩气储层脆性地震预测方法及装置 |
| WO2018102043A1 (en) * | 2016-12-02 | 2018-06-07 | Exxonmobil Upstream Research Company | Method for estimating petrophysical properties for single or multiple scenarios from several spectrally variable seismic and full wavefield inversion products |
| CN107092029A (zh) * | 2017-04-26 | 2017-08-25 | 中国石油大学(北京) | 一种地震反演方法和装置 |
Non-Patent Citations (4)
| Title |
|---|
| ARILD BULAND，等: "Bayesian time-lapse inversion", 《GEOPHYSICS》 * |
| LIN ZHOU，等: "Prestack AVA inversion of exact Zoeppritz equations based on modified Trivariate Cauchy distribution", 《JOURNAL OF APPLIED GEOPHYSICS》 * |
| 刘明阳，等: "基于多波时移地震反演的油藏监测研究", 《中国地球物理2013》 * |
| 王守东，等: "时移地震资料贝叶斯AVO波形反演", 《地球物理学报》 * |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN110661548A (zh) * | 2019-09-13 | 2020-01-07 | 广东石油化工学院 | 一种利用l1模反演的plc信号滤波方法和系统 |
| CN110471106A (zh) * | 2019-09-20 | 2019-11-19 | 西南石油大学 | 一种基于滤波器设计的时移地震反演方法 |
| CN110471106B (zh) * | 2019-09-20 | 2022-03-11 | 西南石油大学 | 一种基于滤波器设计的时移地震反演方法 |
| CN111736208A (zh) * | 2020-06-24 | 2020-10-02 | 重庆大学 | 变权重联合P波和S波初至数据的微震事件Bayes定位方法、系统及介质 |
| CN111948712A (zh) * | 2020-08-10 | 2020-11-17 | 中海石油(中国)有限公司 | 一种基于深度域地震记录的叠前线性反演方法 |
| CN111948712B (zh) * | 2020-08-10 | 2021-08-10 | 中海石油(中国)有限公司 | 一种基于深度域地震记录的叠前线性反演方法 |
| CN115951405A (zh) * | 2023-03-08 | 2023-04-11 | 中国地震局地球物理研究所 | 一种基于稀疏化卷积和自适应矩估计的速度变化成像方法 |
| CN115951405B (zh) * | 2023-03-08 | 2023-05-09 | 中国地震局地球物理研究所 | 一种基于稀疏化卷积和自适应矩估计的速度变化成像方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN110187384B (zh) | 2021-07-20 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110187384A (zh) | 贝叶斯时移地震差异反演方法及装置 | |
| CN107589448B (zh) | 一种多道地震记录反射系数序列同时反演方法 | |
| EP0750203B1 (en) | Subsurface modeling from seismic data and secondary measurements | |
| CN104597490B (zh) | 基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法 | |
| Tamburini et al. | Retrieving surface deformation by PSInSAR™ technology: A powerful tool in reservoir monitoring | |
| CN104749624B (zh) | 一种地震岩相识别及其不确定性定量评价同步实现方法 | |
| Wang et al. | Current developments on micro-seismic data processing | |
| CN105388518A (zh) | 一种质心频率与频谱比联合的井中地震品质因子反演方法 | |
| CN108645994A (zh) | 一种基于多点地质统计学的地质随机反演方法及装置 | |
| Li et al. | CO2 storage monitoring based on time-lapse seismic data via deep learning | |
| Pevzner et al. | Feasibility of time-lapse seismic methodology for monitoring the injection of small quantities of CO2 into a saline formation, CO2CRC Otway Project | |
| Witten et al. | Microseismic image-domain velocity inversion: Marcellus Shale case study | |
| CN106950596B (zh) | 一种基于子波迭代估计的有限差分对比源全波形反演方法 | |
| Biswas et al. | Density inversion from seismic using a trans-dimensional approach: A field dataset example | |
| CN107678065B (zh) | 提高地震分辨率的保构造井控空间反褶积方法和装置 | |
| CN105929453B (zh) | 一种具有信道衰落的无穷分布时滞神经网络系统的状态估计方法 | |
| CN108680957B (zh) | 基于加权的局部互相关时频域相位反演方法 | |
| Wang et al. | Geophysical and production data history matching based on ensemble smoother with multiple data assimilation | |
| Huang et al. | Data assimilated time-lapse viscoacoustic full-waveform inversion for monitoring CO2 geological storage | |
| Pereira et al. | Integrating of initial guess models into geostatistical seismic inversion methodologies | |
| Caers | Geostatistics: from pattern recognition to pattern reproduction | |
| Ren et al. | On application issues of automatic dispersion curves picking by machine learning | |
| Huang et al. | New seismic attribute: Fractal scaling exponent based on gray detrended fluctuation analysis | |
| CN110568488A (zh) | 一种基于非线性混沌算法的生物礁储层识别方法 | |
| Maity et al. | Novel hybrid ANN autopicker for hydrofrac data: A comparative study |
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 |