一种页岩储层脆性和各向异性参数预测方法及系统

本发明属于油气、页岩气勘探地震和岩石力学参数预测领域，尤其涉及页岩储层脆性和各向异性参数预测方法及系统。

背景技术：

本部分的陈述仅仅是提供了与本发明相关的背景技术信息，不必然构成在先技术。

地震反演是获得地下介质内部图像、对储层进行精细描述的有效方法，也是高分辨率地震勘探的最终表现形式，地震数据反演很大程度上提高了储层表征的价值。随着地震研究的重点由勘探逐渐向开发转移，以及由常规油气向非常规油气尤其是页岩气的转移，通过地震反演等手段来揭示地下油气藏的精细分布特征，对油气藏储层进行精雕细刻的技术也受到越来越多的关注。

页岩储层的水平薄互层构造可以等效为vti介质。从目前的研究来看，地震反演研究主要涉及基于zoeppritz方程的avo反演方法和全波形反演方法。基于精确zoeppritz方程和全波形反演方法能利用丰富的叠前信息预测模型参数，但其计算量巨大，在反演尺度和计算效率上不能满足实际油藏储层精细表征要求，尤其对于实际三维大偏移距地震数据。

针对vti介质，传统的基于zoeppritz方程近似反射系数是一个包含弹性参数和弱各向异性参数的表达式，虽然通过等价改写能重新表示为5个模型参数的线性表达式，但是5个未知参数造成的病态反演问题限制了方法的应用。并且，常规反演获得的弹性参数难以直接对页岩脆性进行刻画，岩石力学参数的间接计算会引入累积误差。

技术实现要素：

为了解决在vti介质的页岩储层脆性和各向异性参数预测中，传统的基于zoeppritz方程近似反射系数中未知参数过多，造成反演病态，并且岩石力学参数不能直接估计的问题，本发明的第一个方面提供一种页岩储层脆性和各向异性参数预测方法，其能够满足地震叠前反演对vti页岩储层脆性和各向异性特征的预测。

为了实现上述目的，本发明采用如下技术方案：

一种页岩储层脆性和各向异性参数预测方法，包括：

采集工区内的测井数据，并基于采集的测井数据确定地震子波；

构建含三属性参数的vti介质反射系数方程，确定三属性参数；

基于测井数据和三属性参数获取先验分布函数、初始三属性参数模型；

根据初始三属性参数模型、地震子波确定反演残差；

基于先验分布函数和似然函数构建最大后验概率意义下的反演目标函数；

根据反演目标函数和反演残差确定最优三属性参数模型；

依据最优三属性参数模型，预测页岩储层脆性、各向异性参数。

进一步地，所述构建含三属性参数的vti介质反射系数方程的步骤为：

利用弹性波逆散射理论、波恩近似及稳相法，推导出基于纵横波模量的vti介质反射系数方程；

利用纵横波模量与杨氏模量和泊松比间的定量关系，将基于纵横波模量的vti介质反射系数方程等价改写为基于岩石力学参数的vti介质反射系数方程；

通过对基于岩石力学参数的vti介质反射系数方程进行整合化简，获得含三属性参数的vti介质反射系数方程。

进一步地，所述确定地震子波的步骤为：

采集工区内的地震叠前道集、实际井旁角度域地震道集；

基于地震叠前道集和测井数据，采取统计方法提取依赖于入射角度的地震子波；

将测井数据为输入模型，利用vti介质反射系数方程正演模拟角度域的pp波道集；

将pp波道集与实际井旁角度域地震道集对比，得到振幅缩放因子；

将振幅缩放因子应用于所述的地震子波。

进一步地，所述先验分布函数的获取过程具体为：基于测井数据和三属性参数，统计三属性参数的先验信息，形成先验分布函数。

进一步地，所述三属性参数的先验信息包括三属性参数的均值和三属性参数的协方差矩阵。

进一步地，所述初始三属性参数模型是通对测井数据进行平滑处理并沿解释层位外推得到的。

进一步地，所述根据初始三属性参数模型、地震子波确定反演残差的具体步骤为：

获取实际角道域地震道集；

将所述初始三属性参数模型与敏感度矩阵做乘积得到反射系数矩阵，并与地震子波矩阵进行褶积，得到角度域的角道集；

将实际角道域地震道集与所述角度域的角道集作差，得到反演残差。

为了解决上述问题，本发明的第二个方面提供一种页岩储层脆性和各向异性参数预测系统，其能够满足地震叠前反演对vti页岩储层脆性和各向异性特征的预测。

为了实现上述目的，本发明采用如下技术方案：

一种页岩储层脆性和各向异性参数预测系统，包括：

地震子波确定模块，用于采集工区内的测井数据，并基于采集的测井数据确定地震子波；

三属性参数确定模块，用于构建含三属性参数的vti介质反射系数方程，确定三属性参数；

先验分布函数和初始三属性参数模型获取模块，用于基于测井数据和三属性参数获取先验分布函数、初始三属性参数模型；

反演残差确定模块，用于根据初始三属性参数模型、地震子波确定反演残差；

反演目标函数构建模块，用于基于先验分布函数和似然函数构建最大后验概率意义下的反演目标函数；

最优三属性参数模型确定模块，用于根据反演目标函数和反演残差确定最优三属性参数模型；

预测模块，用于依据最优三属性参数模型，预测页岩储层脆性、各向异性参数。

本发明的第三方面提供一种电子设备，包括存储器和处理器以及存储在存储器上并在处理器上运行的计算机指令，所述计算机指令被处理器运行时，完成第一方面所述方法的步骤。

本发明的第四方面提供一种计算机可读存储介质，用于存储计算机指令，所述计算机指令被处理器执行时，完成第一方面所述方法的步骤。

本发明的有益效果是：

第一、本发明所述的地震叠前vti介质avo反演基于zoeppritz方近似方程，相对于基于精确zoeppritz方程和波动方程的非线性反演而言，基于近似方程的线性反演计算效率高，反演结果稳定，反演参数目的性强，具有很好的工业生产应用价值。

第二、本发明建立的含三属性参数的vti介质反射系数方程，相比包含5个模型参数的传统vti介质反射系数方程，该方程能有效提高反演的稳定性和精度。

第三、本发明通过三变量高斯先验分布引入模型参数之间的相关性，在降低反演不确定性的同时，提高了反演的精度。

第四、本发明通过对新三属性参数的交会关系分析，明确新属性参数间物理意义，指导储层脆性、强各向异性发育段的预测及评价。

附图说明

构成本发明的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。

图1是本发明实施例一所述页岩储层脆性和各向异性参数预测方法流程图；

图2是本发明实施例一a地区测井曲线图；

图3是本发明实施例一cdp叠加剖面图；

图4是本发明实施例一新属性参数反演剖面图；

图5是本发明实施例一测井位置反演曲线对比及合成道集与实际道集对比图；

图6是本发明实施例一测井位置反演曲线合成道集与实际道集对比图；

图7是本发明实施例一基于新三属性间交会关系的测井位置脆性参数和各向异性参数曲线计算结果图；

图8是本发明实施例一基于新三属性间交会关系的密度、泊松比和各向异性参数剖面计算结果图。

具体实施方式

下面结合附图与实施例对本发明作进一步说明。

应该指出，以下详细说明都是例示性的，旨在对本发明提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

实施例1

本实施例提供一种页岩储层脆性和各向异性参数预测方法。

针对vti介质的页岩储层脆性和各向异性参数预测，传统的基于zoeppritz方程近似反射系数中未知参数过多，造成反演病态，并且岩石力学参数不能直接估计的问题，提出了一种基于vti介质avo反演的页岩储层脆性和各向异性参数评价方法。利用弹性波逆散射理论、波恩近似及稳相法，推导出基于纵横波模量的vti介质地震波反射系数方程，并利用纵横波模量与杨氏模量和泊松比间的定量关系，将上述方程等价改写为基于岩石力学参数的vti介质反射系数方程。通过对方程进行整合化简，获得含3个新属性参数的线性反射系数方程。在贝叶斯框架下，提出了vti介质新三属性参数线性反演方法，实现了新属性参数的准确预测。通过建立新属性参数间的交会关系，进一步识别页岩储层位置及高脆性、强各向异性储层段，最终实现页岩储层脆性和各向异性特征的准确识别。

如图1所示，本实施例提供的一种页岩储层脆性和各向异性参数预测方法，具体步骤包括：

s1：采集工区内的地震叠前道集、测井数据以及实际井旁角度域地震道集；基于所述的地震叠前道集、测井数据以及实际井旁角度域地震道集确定地震子波。

基于实际的地震叠前道集和测井数据提取角度依赖的子波，得到子波矩阵w；基于测井数据和vti介质反射系数方程正演模拟地震角道集，并结合实际井旁角度域地震道集确定振幅缩放因子。

s101:本发明假设反演之前地震子波已知，因而，需要基于实际的地震叠前道集和测井数据采取统计方法提取子波，子波在传播过程中受地层的影响会发生波形或频率变化，提取依赖于入射角度的地震子波能有效提高振幅匹配程度。

s102:实际的地震振幅往往是相对值，采用vti介质反射系数方程正演模拟的地震数据振幅与实际振幅存在一定数值差异。以测井数据为输入模型利用vti介质反射系数方程(这里的vti介质反射系数方程可以为基于纵横波模量的vti介质地震波反射系数方程、基于岩石力学参数的vti介质反射系数方程、含三属性参数的vti介质反射系数方程中的任意一种，本实施例使用的是含三属性参数的vti介质反射系数方程)正演模拟角度域的pp波道集，与实际井旁角度域地震道集对比，计算振幅缩放因子，并应用于所提取的地震子波，达到模拟记录与实际记录的振幅匹配。当地震数据信噪比较高时，为角道集的每一道使用统一的振幅缩放因子，以保证振幅随偏移距的变化关系；当信噪比低时，可分近、中、远偏移距分别计算振幅缩放因子，保证模拟记录与实际记录的最佳匹配，减少噪声对反演过程的影响。

s2：建立含三属性参数的vti介质反射系数方程；即，推导vti介质三属性参数的反射系数方程，具体详细叙述如下：

s201：利用弹性波逆散射理论、波恩近似及稳相法，推导出基于纵横波模量的vti介质反射系数方程，

其中，ρ是密度，θ为入射角，s是散射函数，具体结果如下，

其中，pμ(θ)＝-2gsin²θ，和和是纵横波模量反射系数，是密度反射系数，δ,ε是弱各向异性参数。

s202：利用纵横波模量与杨氏模量和泊松比间的定量关系，将上述方程等价改写为基于岩石力学参数的vti介质反射系数方程

其中，和是杨氏模量和泊松比反射系数，

s203：通过对方程(3)进行整合化简，获得含三属性参数的vti介质反射系数方程

其中，三属性参数函数包括：第一属性参数为第二属性参数为第三属性参数为d(a)表示即下层与上层属性a的比，e为杨氏模量，ρ为密度，σ为泊松比；可以看出，第一属性参数与杨氏模量、密度和泊松比相关；第二属性参数近似等于杨氏模量；第三属性参数与第一属性参数的相关性依赖于密度和弱各向异性参数。因此，第二属性参数用于预测杨氏模量；联合第二属性参数、第一属性参数和密度计算预测泊松比；联合第一属性参数和第三属性参和密度计算各向异性参数。

s3：基于测井数据，统计新三属性参数的先验信息，形成符合该工区的模型参数先验分布函数；基于工区内所有测井数据统计模型参数的先验信息，所述三属性参数的先验信息包括三属性参数的均值和三属性参数的协方差矩阵；具体详细叙述如下：

本发明采用三变量高斯分布函数作为先验分布函数；基于工区内所有测井数据(纵横波速度(vp,vs)、密度(ρ)、弱各向异性参数ε和δ)统计各模型参数(a，b，c)的均值，求取各参数(a，b，c)的自相关系数和互相关系数，构建三参数相关的协方差矩阵，形成符合该工区的模型参数先验分布函数；在后续的反演目标函数中三变量高斯分布函数相应的规则化表达式，即，先验分布函数为：

其中，m＝[a,b,c]^t为参数模型，μ和cm分别为所统计的模型参数的均值和三参数(a，b，c)协方差矩阵，三参数(a，b，c)的推导过程详见公式(1)-(4)；表示杨氏模量、密度和泊松比σ的t2次方的乘积；表示杨氏模量乘上研究表明在中度和强各向异性情况下，则，b可等于杨氏模量e；可以看出与a的相关性依赖于密度和弱各向异性参数ε；其中，

假设不同时间点的模型参数(a，b，c)互不相关，可以得到cm的表达式为：

其中为第i个时间点属性a的自相关系数，为第i个时间点属性b与c的互相关系数，其它类推；这些系数通过对工区内测井数据进行时间延迟统计得到；由于三变量高斯分布通过协方差矩阵融合属性a、b以及c之间的相关性，降低了三个属性参数之间的不确定性。

s4：对所述测井数据进行平滑处理并沿解释层位外推，建立初始三属性参数模型，所述初始三属性参数模型的建立利用三维空间插值方法。具体详细叙述如下：

建立弹性参数模型主要利用三维空间插值方法，其的技术流程是首先利用散点插值的方法对各个层位(通过商业软件操作人工或自动拾取的具有明显分界面的层位)的数据(拾取的时间数据)进行插值，完成地质层位建模，然后根据地质层位进行测井数据横向插值，即将测井信息进行横向插值，计算得到地下每个点上的新属性参数值，完成初始新属性参数(a，b，c)建模的任务。

s5：根据初始三属性参数模型、地震子波确定反演残差，具体包括：

将所述初始三属性参数模型作为输入，与敏感度矩阵做乘积得到反射系数矩阵，并与地震子波矩阵进行褶积，得到角度域的角道集；

采集实际角道域地震道集；将实际角道域地震道集与所述角度域的角道集作差，得到反演残差。

s6：基于所述的模型参数先验分布函数和似然函数构建最大后验概率意义下的反演目标函数；

s7：根据所述的反演目标函数和反演残差确定最优三属性参数模型。

基于bayesian原理、先验信息和正演算子构建最大后验概率意义下的反演目标函数，通过求目标函数关于模型参数的导数，并求解目标函数的极小值等价于令其导数等于零，获得最优三属性参数模型。

贝叶斯方法是用来计算条件概率的一种概率统计方法，数学上它的表达式为

p(m|d)∝p(d|m)p(m)(7)

其中，p(m|d)是后验概率分布；p(d|m)是从模型参数空间映射到观测数据空间的似然函数；p(m)是在观测之前对模型参数的先验了解，称为先验分布；m为模型参数，d为地震观测数据。

假设噪声相互独立并且服从高斯分布，观测数据的似然函数可以表示为

其中，是噪声的方差，n是数据采样点数,g＝wl为正演算子，w为子波矩阵，l是敏感度矩阵

假设模型参数之间相互独立，则先验分布形式为，

其中，cm是模型参数的方差；

将反演问题转换为直接求最大后验概率解，即求解以下目标函数的极小值，所述反演目标函数为

其中，是数据的方差，通过求目标函数关于模型参数的导数，

等价于令其导数等于零，整理得到，

其中，μh是超参数。

s8：依据最优三属性参数模型，预测vti介质页岩储层脆性、各向异性参数；第二属性参数用于预测杨氏模量；联合第二属性参数、第一属性参数和密度计算预测泊松比；联合第一属性参数和第三属性参和密度计算各向异性参数。即，通过对新属性参数间的交会图分析，评估页岩储层位置及高脆性发育段，具体包括；

s801：利用估计的第二属性参数预测杨氏模量，研究表明在中度和强各向异性情况下，则，反演的b可近似表示vti介质的杨氏模量e。

s802：利用反演得到的第一属性参数b(或杨氏模量)以及估测的第一属性参数获得泊松比；其中，背景弹性特征可从测井数据计算得到，密度的获取方法：一种是根据经验公式，二是由于小角度地震反射对各向异性不敏感，因此密度可以利用小角度道集进行各向同性反演获得；

s803：利用反演得到的第一属性参数a和第三属性参数c，首先计算a/c＝ρ²/e^2ε，接着利用s802获取的密度曲线或剖面可计算各向异性参数。

通过以上具体步骤的处理，实现了vti介质页岩储层脆性、各向异性参数预测及评估的难题。

为了验证本发明所述页岩储层脆性和各向异性参数预测方法的应用效果，下面以四川页岩某研究区实际叠前反演过程为例，进行分析。

如图2所示为a地区测井曲线信息，包括纵横波速度、密度、各向异性参数。测井数据用于计算新属性参数、反演初始模型、正则化模型参数协方差以及新属性参数间交会关系。

如图3所示为cdp叠加剖面图，图3中的竖线为测井位置，从图3可以看出在2210-2255ms的位置具有强反射轴。

如图4所示为新属性参数反演剖面，图4中的曲线为测井数据，从图4可以看出a、b和c三属性在2210-2255ms具有相对低值，与图3的cdp叠加剖面的强反射轴和测井曲线有较好对应。

如图5所示为测井位置反演曲线对比，图5中的实线为反演曲线，点实线为实际测井曲线，虚线为初始平滑模型，从图可以看出反演曲线与测井曲线有很好一致性，证实了方法的有效性。

如图6为井位的反演曲线合成道集与实际道集对比，其中，(a)实际道集，(b)对实际道集进行k-l变换后的拉平道集，(c)反演曲线合成道集，从图可以看出反演曲线合成道集与实际拉平际道集有很好一致性，进一步证实了方法的有效性。

图7为基于三属性间交会关系的测井位置脆性参数和各向异性参数曲线计算结果，其中，(a)反演b属性(实线)与实际b属性(点实线)、实际杨氏模量(虚线)曲线对比；(b)反演泊松比(实线)与实际泊松比(点实线)曲线对比；(c)反演a/c(实线)与实际a/c(点实线)曲线对比；(d)反演各向异性参数(实线)与实际各向异性参数(点实线)曲线对比。从图可以看出反演的杨氏模量、泊松比和各向异性参数曲线与相应的实际测井曲线有很好一致性，从反演曲线图可以看出页岩储层位置，高脆性和强各向异性发育段在2210-2255ms。

图8为基于新三属性间交会关系的泊松比和各向异性参数剖面计算结果，图中曲线为测井数据，从图可以看出泊松比和各向异性参数在2210-2255ms具有相对低值，与图3的cdp叠加剖面的强反射轴和测井曲线有较好对应，可以合理识别页岩储层位置，高脆性和强各向异性发育段。

实施例2

本实施例提供一种页岩储层脆性和各向异性参数预测系统，包括：

地震子波确定模块，用于采集工区内的测井数据，并基于采集的测井数据确定地震子波；

三属性参数确定模块，用于构建含三属性参数的vti介质反射系数方程，确定三属性参数；

先验分布函数和初始三属性参数模型获取模块，用于基于测井数据和三属性参数获取先验分布函数、初始三属性参数模型；

反演残差确定模块，用于根据初始三属性参数模型、地震子波确定反演残差；

反演目标函数构建模块，用于基于先验分布函数和似然函数构建最大后验概率意义下的反演目标函数；

最优三属性参数模型确定模块，用于根据反演目标函数和反演残差确定最优三属性参数模型；

预测模块，用于依据最优三属性参数模型，预测页岩储层脆性、各向异性参数。

实施例3

本实施例还提供了一种电子设备，包括存储器和处理器以及存储在存储器上并在处理器上运行的计算机指令，所述计算机指令被处理器运行时，完成实施例1所述方法的步骤。

实施例4

本实施例还提供了一种计算机可读存储介质，用于存储计算机指令，所述计算机指令被处理器执行时，完成实施例1所述方法的步骤。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。