CN105629303A - 基于岩石物理的叠前裂缝定量预测方法及系统 - Google Patents

Abstract

本发明提供了一种基于岩石物理的叠前裂缝定量预测方法及系统，所述叠前裂缝定量预测方法包括：根据待测储层的叠前CMP道集中的方位角及入射角的大小，将所述叠前CMP道集划分为多个方位角道集：方位角道集1，方位角道集2，…，方位角道集n；对各所述方位角道集进行叠加、偏移处理，计算各方位角道集的纵波反射系数；根据Ruger方程精确式以及各所述方位角道集及其纵波反射系数，计算待测储层中裂缝的切向各异性系数Δ_T；根据所述切向各异性系数Δ_T计算待测储层的裂缝密度e。本发明能够准确地实现裂缝密度的定量估算，进而准确地进行油气藏预测和指导开发井位的部署，达到了降低油气藏预测风险的目的。

Description

基于岩石物理的叠前裂缝定量预测方法及系统

技术领域

本发明涉及石油、天然气地震勘探反演和定量解释领域，尤其是涉及了一种基于岩石物理的叠前裂缝定量预测方法及系统。

背景技术

AVOAz(Amplitudeversusoffsetandazimuth)反演就是利用反射纵波振幅随入射角和方位角的变化规律反演地下裂缝密度以及裂缝发育方向。大量地球物理学家曾对各向异性介质中反射、透射近似关系式进行了研究，当前最为出名且被大家广泛应用的是Ruger基于弱各向异性的概念，结合Thomsen各向异性参数，提出的HTI介质中纵波反射系数随方位角和入射角变化的公式，即Ruger方程精确式，具体形式见式(1)。Ruger方程精确式是入射角i、方位角φ、Thomsen各向异性三参数(ε、δ、γ)以及介质弹性参数(纵波速度V_p、横波速度V_s、密度ρ、纵波阻抗Z和横波切向模量G等)的函数。

式中，i和φ分别表示纵波的入射角和方位角；φ_sym为HTI介质对称轴方向的方位角；(i,φ)为与入射角i和方位角φ有关的界面纵波反射系数；Z＝ρα为纵波波阻抗；ρ为介质密度，g/cm³；α为纵波速度，m/s；为界面上下的波阻抗之差与平均波阻抗的比值；G＝ρβ²为横波切向模量，β为横波速度；ε、δ和γ为Thomsen的各向异性参数；Δ[·]表示上、下界面物理量之差；[·]表示上、下界面物理量的均值。上标V代表HTI介质。由公式(1)可以得到不同参数条件下HTI介质的反射系数随入射角和方位角的变化情况，进而可进一步探讨和研究HTI介质的各向异性。

纵波各向异性裂缝预测方法作为目前应用比较广泛、效果较好的一种裂缝预测方法。大量国内外学者都进行了深入的探索，其中最早出现的较为经典的是传统的椭圆拟合法。该方法原理简单、易于实现，被广泛应用于各种商业化软件中，并成功用于油田实际资料的解释中。椭圆拟合法是在小入射角前提下，忽略Ruger方程精确式中的高阶项部分，从而对Ruger方程进一步简化，得到随方位角变化的AVO梯度项，它们的具体表达式分别为：

B(φ_k)＝B^iso+B^anicos²(φ_k-φ_sym)(2)

式(2)中的φ_k为第k次的地震观测方位角。事实上，式(2)可以近似为一个椭圆，并且地下介质中的裂缝强度越大，由各向异性拟合出的方位椭圆的扁率越大，其长轴或短轴方向代表裂缝走向。因此，理论上只需知道三个或三个以上方位的反射地震数据，就可以实现对目的地层中任意一点的裂缝发育密度和方位的预测。但是，传统的椭圆拟合方法的准确性会受到一系列因素的影响，当入射角较大、信噪比较低、方位角分布不均匀、盖层的物性和各向异性横向变化、上覆介质的透射各向异性较强时，该方法的裂缝预测结果误差较大。

随后本领域又发展了多种预测方法，包括算法上的改进和最终反演参数的选取，都尝试提高裂缝预测的精度，但大部分方法都受到入射角的限制且精度不高，因为这些方法都是在忽略高阶项的基础上发展而来的。再后来，少量学者提出基于Ruger方程三项的反演方法，大大提高了预测精度且克服了入射角的影响，但这样得到的各向异性梯度项、甚至Thomsen的各向异性参数和裂缝密度之间并没有直接的关系。事实上，地下的各向异性并非都只是由裂缝引起的。另一方面，现有的这些方法都不能研究裂缝中充填的流体性质，而这对于后面的开发布井工作的指导至关重要。

在实际的裂缝预测中，由于各种因素的影响，裂缝的密度及发育方向很难准确估计，尤其是裂缝密度的定量预测。传统的椭圆拟合方法由于受到各种因素的制约，不能准确得到裂缝密度和方向，而后发展的基于三项方程的反演方法虽然弥补了椭圆拟合法的多种不足，但其仍不能实现裂缝密度的定量预测。

发明内容

为解决上述技术问题，本发明提供了一种基于岩石物理的叠前裂缝定量预测方法及系统。

本发明一方面提供了一种基于岩石物理的叠前裂缝定量预测方法，所述叠前裂缝定量预测方法包括：

根据待测储层的叠前CMP道集中的方位角及入射角的大小，将所述叠前CMP道集划分为多个方位角道集：方位角道集1，方位角道集2，…，方位角道集n；

对各所述方位角道集进行叠加、偏移处理，计算各方位角道集的纵波反射系数；

根据Ruger方程精确式以及各所述方位角道集及其纵波反射系数，计算待测储层中裂缝的切向各异性系数Δ_T；

根据所述切向各异性系数Δ_T计算待测储层的裂缝密度e。

在一实施例中，根据Ruger方程精确式以及各所述方位角道集及其纵波反射系数，计算待测储层中裂缝的切向各异性系数Δ_T，包括：

根据HTI介质中各向异性参数与裂缝岩石物理参数之间的关系，对Ruger方程精确式进行变换，得到纵波反射系数简化方程：

将所述方位角道集及其对应的纵波反射系数代入所述纵波反射系数简化方程得到矩阵方程：

Ax＝b；

根据所述矩阵方程建立一用于反演所述x的最优化目标函数：

对所述最优化目标函数进行求解，得到所述切向各异性系数Δ_T；

其中，x₄＝Δ_T；x₅＝Δ_N；i和φ分别表示纵波的入射角和方位角；r为与入射角i和方位角φ有关的纵波反射系数；z＝ρα为纵波波阻抗；ρ为介质密度，g/cm³；α为纵波速度，m/s；G＝ρβ²为横波切向模量，β为横波速度；ε、δ和γ为Thomsen的各向异性参数；Δ_T为裂缝岩石物理参数中的切向各异性系数；Δ_N为裂缝岩石物理参数中的法向各异性系数；g＝μ/(λ+2μ)为横纵波速度比的平方，λ、μ为不含裂隙各向同性岩石的拉梅参数；

R_k、P_k及W_k为中间变量，R_k＝sin²φ_ktan²i_k，P_k＝cos²φ_ksin²i_k，W_k＝sin²i_ktan²i_k，r_k＝r_k(i_k,φ_k)为方位角道集k对应的纵波反射系数，b＝[r₁(i₁,φ₁)；r₂(i₂,φ₂)；…；r_n(i_n,φ_n)]^T，x＝[x₁,x₂,x₃,x₄,x₅]^T，k∈[1，n]，n为方位角道集的数量。

在一实施例中，对所述最优化目标函数进行求解，得到所述切向各异性系数Δ_T，包括：

对所述最优化目标函数进行Tikhonov正则化处理，经过Tikhonov正则化处理后的所述最优化目标函数为：

计算经过Tikhonov正则化处理后的所述最优化目标函数的最小二乘解x＝(A^TA+τI)^-1A^Tb，即得所述切向各向异性系数Δ_T；

其中，τ为正则化参数。

在一实施例中，根据所述切向各异性系数Δ_T计算待测储层的裂缝密度e，包括：

根据所述切向各异性系数Δ_T与裂缝密度e之间的关系，按照下式计算所述裂缝密度e：

在一实施例中，所述切向各异性系数Δ_T与裂缝密度e之间的关系是根据Hudson币状裂隙模型及Schenberg线性滑动模型确定的。

本发明另一方面还提供了一种基于岩石物理的叠前裂缝定量预测系统，所述基于岩石物理的叠前裂缝定量预测系统包括：

方位角道集划分单元，用于根据待测储层的叠前CMP道集中的方位角及入射角的大小，将所述叠前CMP道集划分为多个方位角道集：方位角道集1，方位角道集2，…，方位角道集n；

纵波反射系数获取单元，用于对各所述方位角道集进行叠加、偏移处理，计算各方位角道集的纵波反射系数；

岩石物理参数计算单元，用于根据Ruger方程精确式以及各所述方位角道集及其纵波反射系数，计算待测储层中裂缝的切向各异性系数Δ_T；

裂缝密度计算单元，用于根据所述切向各异性系数Δ_T计算待测储层的裂缝密度e。

在一实施例中，所述岩石物理参数计算单元包括：

Ruger方程预处理模块，用于根据HTI介质中各向异性参数与裂缝岩石物理参数之间的关系，对Ruger方程精确式进行变换，得到纵波反射系数简化方程：

矩阵方程生成模块，用于将所述方位角道集及其对应的纵波反射系数代入所述纵波反射系数简化方程得到矩阵方程：

Ax＝b；

最优化目标函数建立模块，用于根据所述矩阵方程建立一用于反演所述x的最优化目标函数：

求解模块，用于对所述最优化目标函数进行求解，得到所述切向各异性系数Δ_T；

其中，x₄＝Δ_T；x₅＝Δ_N；i和φ分别表示纵波的入射角和方位角；r为与入射角i和方位角φ有关的纵波反射系数；z＝ρα为纵波波阻抗；ρ为介质密度，g/cm³；α为纵波速度，m/s；G＝ρβ²为横波切向模量，β为横波速度；ε、δ和γ为Thomsen的各向异性参数；Δ_T为裂缝岩石物理参数中的切向各异性系数；Δ_N为裂缝岩石物理参数中的法向各异性系数；g＝μ/(λ+2μ)为横纵波速度比的平方，λ、μ为不含裂隙各向同性岩石的拉梅参数；

R_k、P_k及W_k为中间变量，R_k＝sin²φ_ktan²i_k，P_k＝cos²φ_ksin²i_k，W_k＝sin²i_ktan²i_k，r_k＝r_k(i_k,φ_k)为方位角道集k对应的纵波反射系数，b＝[r₁(i₁,φ₁)；r₂(i₂,φ₂)；…；r_n(i_n,φ_n)],x＝[x₁,x₂,x₃,x₄,x₅]^T，k∈[1，n]，n为方位角道集的数量。

在一实施例中，所述求解模块具体用于：

对所述最优化目标函数进行Tikhonov正则化处理，经过Tikhonov正则化处理后的所述最优化目标函数为：

计算经过Tikhonov正则化处理后的所述最优化目标函数的最小二乘解x＝(A^TA+τI)^-1A^Tb，即得所述切向各向异性系数Δ_T；

其中，τ为正则化参数。

在一实施例中，所述裂缝密度计算单元具体用于：

根据所述切向各异性系数Δ_T与裂缝密度e之间的关系，按照下式计算所述裂缝密度e：

在一实施例中，所述切向各异性系数Δ_T与裂缝密度e之间的关系是根据Hudson币状裂隙模型及Schenberg线性滑动模型确定的。

本发明方法在基于Ruger方程精确式的基础上，结合岩石物理模型，能够准确地实现裂缝密度的定量估算，进而准确地进行油气藏预测和指导开发井位的部署，达到了降低油气藏预测风险的目的。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1A为本发明实施例在不同入射角时，忽略Ruger方程精确式中高阶项和不忽略高阶项两种情况下反射系数随方位角变化的曲线对比图，图中箭头指示入射角增大方向；

图1B为本发明实施例sin²(i)和sin²(i)tan²(i)分别随入射角变化的曲线；

图2A为本发明实施例纵波各向异性指数对拟合椭圆扁率的影响；

图2B为本发明实施例纵波变异系数对拟合椭圆扁率的影响；

图2C为本发明实施例横波各向异性指数对拟合椭圆扁率的影响；

图3A为本发明实施例拟合椭圆扁率随盖层纵波速度及横波速度变化的曲面；

图3B为本发明实施例拟合椭圆扁率随反射界面上下介质的阻抗比及速度比变化的曲面；

图4为本发明实施例各向异性梯度项随盖层纵波速度和横波速度变化的曲面；

图5A为本发明实施例裂缝饱含水的情况下，法向各向异性系数Δ_N随横纵波速度比和裂缝密度变化的曲线，图中箭头指示裂缝密度增大方向；

图5B为本发明实施例裂缝饱含水的情况下，切向各向异性系数Δ_T随横纵波速度比和裂缝密度的变化曲线，图中箭头指示裂缝密度增大方向；

图5C为本发明实施例裂缝含油水混合物的情况下，法向各向异性系数Δ_N随横纵波速度比和裂缝密度变化的曲线，图中箭头指示裂缝密度增大方向；

图5D为本发明实施例裂缝含油水混合物的情况下，切向各向异性系数Δ_T随横纵波速度比和裂缝密度变化的曲线，图中箭头指示裂缝密度增大方向；

图5E为本发明实施例裂缝含气水混合物的情况下，法向各向异性系数Δ_N随横纵波速度比和裂缝密度变化的曲线，图中箭头指示裂缝密度增大方向；

图5F为本发明实施例裂缝含气水混合物的情况下，切向各向异性系数Δ_T随横纵波速度比和裂缝密度变化的曲线，图中箭头指示裂缝密度增大方向；

图6A为本发明实施例裂缝中含油水混合物时，Δ_N随含水饱和度和裂缝纵横比变化的曲线，图中箭头指示裂缝纵横比增大方向；

图6B为本发明实施例裂缝中含气水混合物时，Δ_N随含水饱和度和裂缝纵横比变化的曲线，图中箭头指示裂缝纵横比增大方向；

图7为本发明实施例基于岩石物理的叠前裂缝定量预测方法的流程示意图；

图8为本发明实施例计算裂缝的岩石物理参数的流程示意图；

图9A为本发明实施例一模型的真实裂缝密度的分布情况；

图9B为对图9A所示模型进行反演得到的反射界面上裂缝密度的分布情况；

图9C为本发明实施例反演出的裂缝密度与真实裂缝密度的残差；

图10A为本发明实施例另一模型的反射界面上真实裂缝密度的分布情况；

图10B为对图10A所示模型进行反演得到的反射界面上裂缝密度的分布情况；

图10C为本发明实施例反演出的裂缝密度与真实裂缝密度的残差；

图11为本发明实施例基于岩石物理的叠前裂缝定量预测系统的结构示意图；

图12为本发明实施例岩石物理参数计算单元3的结构示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

背景技术中提到的裂缝预测技术都忽略了Ruger方程精确式中高阶项，而这样会受到入射角的限制，并且影响最终预测的精度。为了直观显示高阶项在Ruger方程精确式中所起的作用，本发明实施例设置了三组大小不同的入射角：20°、30°、40°，根据忽略高阶项的Ruger方程精确式计算得到反射系数大小如图1A所示。从图1A中可以看出，在入射角较小的情况下，反射系数的真实值与忽略高阶项后的反射系数基本重合，因此在入射角较小的情况下可以忽略高阶项的影响。但当入射角增大时，反射系数的大小严重偏离真实值，如此就必须考虑高阶项的影响，追其根源是由高阶项的系数sin²(i)tan²(i)决定的。

图1B显示了sin²(i)和sin²(i)tan²(i)分别随入射角变化的曲线。由图1B可以看出，随着入射角的增大，高阶项占反射系数真实值的比重越来越大，当入射角超过45°时，高阶项在反射系数中占据主导地位。而在实际的勘探生产中，随着深层致密油气的勘探需求和采集技术的不断发展，大偏移距的数据越来越多，并且大偏移距数据体现的方位各向异性更强，传统的两项反演不能达到要求。此外，高阶项还能对裂缝方位的预测起到指导作用，这就有必要开展基于Ruger方程精确式的AVOAz反演方法。

从最早的传统椭圆拟合方法发展到现在的Ruger方程三项反演的方法，人们都尝试得到各种参数因子用以表征地下裂缝的密度，主要包括椭圆扁率、各向异性梯度项和Thomsen各向异性三参数等，但这些参数因子甚至得到的各向异性参数都不能直接代表地下的裂缝密度。换句话说，上述参数因子都只是定性的代表地下裂缝密度的相对大小，而不是对裂缝密度进行准确的定量描述。

综上所述，传统的椭圆拟合方法由于受到各种因素的制约，不能准确得到裂缝密度和方向，而后发展的基于三项方程的反演方法虽然弥补了椭圆拟合法的多种不足，但其仍不能实现裂缝密度的定量预测。

本发明结合岩石物理，对用于表征裂缝密度因子的影响因素进行了分析。

首先，本发明分析了用椭圆扁率来表征裂缝密度的准确性。假定Thomsen各向异性三参数可以代表地下的裂缝密度。计算得到单一Thomsen各向异性参数中的纵波各向异性指数ε、纵波变异系数δ及横波各向异性指数γ对拟合椭圆扁率的影响，分别如图2A～图2C所示。

图2A为在纵波变异系数δ及横波各向异性指数γ固定不变的情况下，拟合椭圆扁率随纵波各向异性指数ε变化的曲线。图2B为在纵波各向异性指数ε及横波各向异性指数γ固定不变的情况下，拟合椭圆扁率随纵波变异系数δ变化的曲线。图2C为在纵波各向异性指数ε及纵波变异系数δ固定不变的情况下，拟合椭圆扁率随横波各向异性指数γ变化的曲线。从图2A～图2C可以得出以下结论：1)介质的各向异性越强(即表示地下介质的裂缝密度越大)时，拟合出的各向异性椭圆的扁率越大，这与通常的认识相吻合。但是，本实施例只研究了单一变量对拟合椭圆扁率的影响，如果当多个变量同时变化时，情况将变得十分复杂，此时拟合椭圆的扁率不一定随着介质各向异性的增强而变大。2)当反射界面上下介质的各向异性差为正各向异性和弱各向异性时，曲线更加接近线性化，导致预测结果要好于负各向异性和强各向异性两种情况，同时这也决定了使用椭圆扁率只能定性而非定量代表地下的裂缝密度。

当忽略反射系数公式—Ruger方程精确式中的高阶项时，可以得出拟合椭圆扁率的近似式：

式中，上标V代表HTI介质。从式(5)可以看出，拟合出椭圆的扁率不仅与界面上下介质的各向异性参数有关，还与界面上下介质的纵横波速度、纵波速度差、纵波阻抗差、横波切向模量差等因素有关，如图3A及图3B所示。图3A为本发明实施例拟合椭圆扁率随盖层纵波速度α及横波速度β变化的曲面，图3B为本发明实施例拟合椭圆扁率随反射界面上下介质的阻抗比及速度比变化的曲面。从图3A及图3B中可以看出，拟合椭圆扁率大并不一定是由各向异性引起的。因此，盖层的物性参数也会影响裂缝预测结果。

其次，本发明分析了用各向异性梯度项来表征裂缝密度的准确性。假定Thomsen各向异性三参数可以代表地下的裂缝密度，各向异性梯度项的表达式如下所示：

从式(6)可以看出，各向异性梯度项是各向异性参数的线性函数，这在一定程度上要优于用拟合椭圆扁率来表征裂缝密度的情况。同时也可以由式(6)看出，各向异性梯度项也会受到盖层速度的影响，各向异性梯度项随盖层纵波速度和横波速度变化的情况如图4所示。相较于利用拟合椭圆扁率来表征裂缝密度而言，利用Thomsen各向异性三参数表征地下裂缝密度消除了盖层密度和纵波阻抗等因素的影响，表明利用各向异性梯度项来表征裂缝密度在很大程度上要优于拟合椭圆扁率。

最后，本发明分析了Thomsen各向异性三参数与裂缝密度的关系。关于裂缝岩石物理模型的研究很多，较为常用的是Hudson的薄币状形状的裂隙模型及Schenberg的线性滑动模型。通过研究两个模型的关系，得到岩石物理参数中法向各向异性系数Δ_N和切向各向异性系数Δ_T与裂缝参数之间的关系：

式(8)中，g＝μ/(λ+2μ)为横纵波速度比的平方，μ为不含裂隙各向同性岩石的拉梅参数，K'和μ'分别为裂隙中充填物的体积模量和剪切模量，e为裂缝的体密度，χ为裂缝的纵横比。

Bakulin等又推导了HTI介质中各向异性参数与上面提到的裂缝岩石物理参数之间的关系：

ε^(V)＝-2g(1-g)Δ_N(9)

δ^(V)＝-2g[(1-2g)Δ_N+Δ_T](10)

这样就建立了裂缝参数(裂缝密度、裂缝纵横比及裂缝充填流体类型)和Thomsen各向异性参数之间的关系。下面本发明对裂缝岩石物理参数中的法向各向异性系数Δ_N和切向各向异性系数Δ_T与裂缝参数的变化特征进行了分析。

图5A～图5F给出了裂缝饱含水、含油水混合物和含气水混合物三种情况下，裂缝岩石物理参数随横纵波速度比和裂缝密度变化的特征。从图5A～图5F中可以看出，无论裂缝饱含水、含油水混合物和含气水混合物，Δ_N和Δ_T都随裂缝密度的增大呈上升趋势。分析图5A、图5C及图5E中的Δ_N曲线变化可发现，裂缝含气、含油和含水时，Δ_N和Δ_T变化幅度均有差异，但是对于Δ_T来说，不论裂缝含气、含油还是含水，其数值变化幅度均相同，这说明Δ_T对裂缝充填物的类型并不敏感，可以认为Δ_T的数值变化主要是由裂缝密度引起的。如果能从地震数据中估测得到Δ_T，则可以直接预测岩石的裂缝密度。而估测出Δ_N以及结合预测的裂缝密度，则可以估测裂缝充填物的类型，为裂缝型储层的勘探和开发提供支持。

从式(8)中可以看出，如果裂缝中充填了流体，其剪切模量为0，那么Δ_T就不受裂缝纵横比的影响。从前面的分析可知，Δ_T的数值也不受充填物类型的影响，因此只需分析Δ_N随裂缝纵横比的变化特征，如图6A及图6B所示。图6A为本发明实施例裂缝中含油水混合物时，Δ_N随含水饱和度和裂缝纵横比变化的曲线，图6B为本发明实施例裂缝中含气水混合物时，Δ_N随含水饱和度和裂缝纵横比变化的曲线，图中箭头指示裂缝纵横比增大方向。如图6A及图6B所示，Δ_N的数值随含水饱和度的增大而减小，这也说明了Δ_N确实与充填物类型有关。同时，对于含油水混合物裂缝来说，Δ_N在裂缝纵横比较大的情况下随含水饱和度的变化明显；对于含气水混合物裂缝来说，Δ_N在裂缝纵横比较小的情况下随含水饱和度的变化明显。

通过前面的分析可知，法向各向异性系数Δ_N不仅与裂缝密度有关，还会受到充填物类型和裂缝纵横比的影响，而切向各向异性系数Δ_T只与裂缝密度有关。考虑到实际中裂缝的纵横比很小，因此可将包含纵横比的项忽略不计，对式(7)和(8)进一步化简，得到裂缝中含油水混合物的情况下Δ_N和Δ_T的表达式：

Δ_N＝0(12)

对式(7)和(8)进一步化简，得到裂缝中含气情况下Δ_N和Δ_T的表达式：

上述两组等式进一步证实了前面的分析结果，即只有切向各向异性系数Δ_T与裂缝密度有关，同时由式(9)、(10)和(11)可知，只有横波各向异性指数γ仅与裂缝密度有关系，但这也只是定性的，进一步说明了用Thomsen各向异性参数来代表裂缝密度具有一定的误差。所以要想准确定量求取裂缝密度，就需要反演得到切向各向异性系数Δ_T的数值。

由以上分析可知，只有计算得到切向各向异性系数Δ_T的值，才能实现裂缝的定量预测，这就启发用裂缝岩石物理参数替代Thomsen各向异性参数来表示方位反射系数以求取切向各向异性系数Δ_T的值。在充分考虑了现有方法的局限性后，本发明提出了一种基于岩石物理的叠前裂缝定量预测方法。

图7为本发明实施例基于岩石物理的叠前裂缝定量预测方法的流程示意图。如图7所示，基于岩石物理的叠前裂缝定量预测方法主要包括以下步骤：

步骤S1、根据待测储层的叠前CMP道集中的方位角及入射角的大小，将该叠前CMP道集划分为多个方位角道集：方位角道集1，方位角道集2，…，方位角道集n。

首先根据方位角的大小将待测储层的叠前CMP道集初步划分为多个数据体，然后在各数据体中按照入射角的大小对本数据体再次划分，得到最终的方位角道集1～方位角道集n。

步骤S2、分别对上述各方位角道集进行叠加、偏移处理，计算各方位角道集的纵波反射系数。

步骤S3、根据Ruger方程精确式以及上述各方位角道集及其对应的纵波反射系数，计算待测储层中裂缝的切向各异性系数Δ_T。

步骤S4、根据上述切向各异性系数Δ_T计算待测储层的裂缝密度e。

本发明分析了高阶项在Ruger方程精确式中不可忽略的重要作用，并详细分析了各种用于表征裂缝密度因子的影响因素，验证了前人方法的局限性，同时给出了能够准确定量求得裂缝密度的一种新的方案和思路，提升了裂缝预测的精度，达到了降低油气藏预测风险的目的。

利用Ruger方程精确式求解切向各异性系数Δ_T时，可按照图8所示步骤进行：

步骤31、根据HTI介质中各向异性参数与裂缝岩石物理参数之间的关系，对Ruger方程精确式进行变换，得到纵波反射系数简化方程。

HTI介质中各向异性参数与上面提到的裂缝岩石物理参数之间的关系如式(9)～式(11)所示，当利用裂缝岩石物理参数替代Thomsen各向异性参数来表示方位反射系数时，将式(9)、(10)和(11)代入式(1)中即可，整理Ruger方程精确式可以得到新的反射系数方程：

对式(16)进行简化，可得如下所示纵波反射系数简化方程：

其中，x₄＝Δ_T；x₅＝Δ_N。

i和φ分别表示纵波的入射角和方位角；r为与入射角i和方位角φ有关的纵波反射系数；z＝ρα为纵波波阻抗；ρ为介质密度，g/cm³；α为纵波速度，m/s；G＝ρβ²为横波切向模量，β为横波速度；ε、δ和γ为Thomsen的各向异性参数；Δ_T为裂缝岩石物理参数中的切向各异性系数；Δ_N为裂缝岩石物理参数中的法向各异性系数；g＝μ/(λ+2μ)为横纵波速度比的平方，λ、μ为不含裂隙各向同性岩石的拉梅参数。

步骤S32、将步骤S1中得到的各方位角道集及步骤S2中得到的纵波反射系数代入式(17)中，得到如下矩阵方程：

Ax＝b(18)

在上述矩阵方程(18)中，各参数矩阵分别如下所示：

R_k、P_k及W_k为中间变量，R_k＝sin²φ_ktan²i_k，P_k＝cos²φ_ksin²i_k，W_k＝sin²i_ktan²i_k，r_k＝r_k(i_k,φ_k)为方位角道集k对应的纵波反射系数，b＝[r₁(i₁,φ₁)；r₂(i₂,φ₂)；…；r_n(i_n,φ_n)]^T,x＝[x₁,x₂,x₃,x₄,x₅]^T，k∈[1，n]，n为方位角道集的数量，向量x为要求解的变量，其包含了两个裂缝岩石物理参数Δ_T及Δ_N。

步骤S33、根据上述矩阵方程(18)建立一用于反演上述变量x的最优化目标函数，以求解矩阵方程(18)。其中最优化目标函数如式(19)所示：

步骤S34、对式(19)所示最优化目标函数进行求解，得到切向各异性系数Δ_T。

为了使反演方程(17)更加具有适定性和抗噪性，本发明实施例利用Tikhonov正则化方法对式(20)所示的最优化目标函数进行正则化处理，以对目标函数的解进行约束。经过正则化处理后的最优化目标函数为：

式中，τ为正则化参数，用于控制残差项与正则化约束项之间的权重。

计算式(20)所示的经过Tikhonov正则化处理后的最优化目标函数的最小二乘解x＝(A^TA+τI)^-1A^Tb，即得切向各向异性系数Δ_T。

在通过计算上述最小二乘解得到切向各向异性系数Δ_T的值后，即可根据式(13)或者式(15)所示的切向各向异性系数Δ_T与裂缝密度e之间的关系，计算裂缝密度e。其中，根据式(13)或者式(15)可得：

在本发明中，切向各向异性系数Δ_T与裂缝密度e之间的关系是根据Hudson币状裂隙模型及Schenberg线性滑动模型确定的，这样得到的裂缝密度既不受地层各向异性、地应力各向异性等因素的影响，也不会因裂缝充填流体类型的变化而变化，它就是地下裂缝真实的定量描述，是常规裂缝预测方法所不能做到的。

由于在推导反演方程时保留了Ruger方程精确式中的高阶项，因此本发明比传统的裂缝反演方法和椭圆拟合法更加精确。

为了更好地理解本发明的基于岩石物理的叠前裂缝定量预测方法及其有益效果，下面结合具体的例子进行说明。

本发明实施例使用了两个叠前方位数据对本发明的基于岩石物理的叠前裂缝定量预测方法进行了测试。在本实施例中，采用HTI介质来模拟实际地层中的裂缝分布情况。在合成方位地震数据时，使用了20°、30°和40°三个入射角，假设每个入射角生成5个方位角的数据，这5个方位角分别为18°、54°、90°、126°和162°。此时，可以得到15个方位角道集，代表每个CDP点具有15次覆盖。

图9A为地下真实裂缝密度的分布情况，图9B显示了反演出的反射界面上裂缝密度的分布情况，图9C为反演出的裂缝密度与真实裂缝密度的残差。从图9A～图9C中可以看出，利用本发明提供的预测方法反演得到的结果与真实模型十分匹配，表明本发明的可靠性和准确性较高。并且，本发明实施例中用于反演的数据中含有40°入射角的数据，由于本发明在推导反演方程中没有忽略高阶项，因此本发明在入射角较大时仍然适用，得到的反演结果也依然很精确。

图10A显示了另一模型反射界面上真实裂缝密度的分布情况，图10B为反演出的反射界面上裂缝密度的分布情况，图10C为反演出的裂缝密度与真实裂缝密度的残差。从图10A～图10C中可以看出，反演结果与真实模型十分匹配。

考虑到实际中大偏移距数据的存在，并且为了利用大偏移距各向异性强以及能够进行深层致密油气勘探的优势，本发明可以稳定地定量求取裂缝发育信息。

现有的裂缝预测方法在求取裂缝密度时，都是计算得到某种表征因子来定向表征地下裂缝密度的，而选用本发明的基于岩石物理的叠前裂缝定量预测方法能够定量得到裂缝密度，从而消除了其他各向异性对裂缝密度的干扰。在测井和岩石物理分析前提下，能够进一步计算得到裂缝流体指示因子，这是常规裂缝反演所不能做到的，特别是对含气裂缝有较强的敏感性，从而为后期的开发布井给予一定的指导，可降低油气藏预测的风险。

基于与图7所示的基于岩石物理的叠前裂缝定量预测方法相同的发明构思，本发明实施例还提供了一种基于岩石物理的叠前裂缝定量预测系统，如下面实施例所述。由于该基于岩石物理的叠前裂缝定量预测系统解决问题的原理与图7中基于岩石物理的叠前裂缝定量预测方法相似，因此该基于岩石物理的叠前裂缝定量预测系统的实施可以参见图7中基于岩石物理的叠前裂缝定量预测方法的实施，重复之处不再赘述。

图11为本发明实施例基于岩石物理的叠前裂缝定量预测系统的结构示意图。如图11所示，基于岩石物理的叠前裂缝定量预测系统包括：方位角道集划分单元1、纵波反射系数获取单元2、岩石物理参数计算单元3及裂缝密度计算单元4。

其中，方位角道集划分单元1，用于根据待测储层的叠前CMP道集中的方位角及入射角的大小，将所述叠前CMP道集划分为多个方位角道集：方位角道集1，方位角道集2，…，方位角道集n。

纵波反射系数获取单元2，用于对各所述方位角道集进行叠加、偏移处理，计算各方位角道集的纵波反射系数。

岩石物理参数计算单元3，用于根据Ruger方程精确式以及各所述方位角道集及其纵波反射系数，计算待测储层中裂缝的切向各向异性系数Δ_T。

裂缝密度计算单元4，用于根据所述切向各向异性系数Δ_T计算待测储层的裂缝密度e。

如图12所示，岩石物理参数计算单元3通常包括Ruger方程预处理模块31、矩阵方程生成模块32、最优化目标函数建立模块33及求解模块34。Ruger方程预处理模块31用于根据HTI介质中各向异性参数与裂缝岩石物理参数之间的关系，对Ruger方程精确式进行变换，得到纵波反射系数简化方程(17)。矩阵方程生成模块32用于将所述方位角道集及其对应的纵波反射系数代入纵波反射系数简化方程(17)得到矩阵方程(18)。最优化目标函数建立模块33用于根据矩阵方程(18)建立一用于反演上述变量x的最优化目标函数，如式(19)所示。求解模块34用于对式(19)所示最优化目标函数进行求解，得到切向各向异性系数Δ_T。

一般地，上述的求解模块34通常对式(19)所示的最优化目标函数进行Tikhonov正则化处理，经过Tikhonov正则化处理后的所述最优化目标函数如式(20)所示，然后计算经过Tikhonov正则化处理后的所述最优化目标函数的最小二乘解x＝(A^TA+τI)^-1A^Tb，即得所述切向各向异性系数Δ_T。

在一实施例中，裂缝密度计算单元4具体用于根据裂缝岩石物理参数Δ_T与裂缝密度e之间的关系，按照式(21)计算待测储层的裂缝密度e。

在一实施例中，上述切向各向异性系数Δ_T与裂缝密度e之间的关系是根据Hudson币状裂隙模型及Schenberg线性滑动模型确定的。

本发明方法在基于Ruger方程精确式的基础上，结合岩石物理模型，能够准确地实现裂缝密度的定量估算，进而准确地进行油气藏预测和指导开发井位的部署，达到了降低油气藏预测风险的目的。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种基于岩石物理的叠前裂缝定量预测方法，其特征在于，所述叠前裂缝定量预测方法包括：

根据待测储层的叠前CMP道集中的方位角及入射角的大小，将所述叠前CMP道集划分为多个方位角道集：方位角道集1，方位角道集2，…，方位角道集n；

对各所述方位角道集进行叠加、偏移处理，计算各方位角道集的纵波反射系数；

根据Ruger方程精确式以及各所述方位角道集及其纵波反射系数，计算待测储层中裂缝的切向各异性系数Δ_T；

根据所述切向各异性系数Δ_T计算待测储层的裂缝密度e。

2.根据权利要求1所述的基于岩石物理的叠前裂缝定量预测方法，其特征在于，根据Ruger方程精确式以及各所述方位角道集及其纵波反射系数，计算待测储层中裂缝的切向各异性系数Δ_T，包括：

根据HTI介质中各向异性参数与裂缝岩石物理参数之间的关系，对Ruger方程精确式进行变换，得到纵波反射系数简化方程：

x₁+x₂sin²i+x₃sin²itan²i+x₄[-g(1+sin²φtan²i)+2g]cos²φsin²i

+x₅gcos²φsin²i[(gsin²φ+g-1)tan²i+2g-1]＝r；

将所述方位角道集及其对应的纵波反射系数代入所述纵波反射系数简化方程得到矩阵方程：

Ax＝b；

根据所述矩阵方程建立一用于反演所述x的最优化目标函数：

$\min E{\left(x\right)}_{obj}=\left|\right|Ax-b|{|}_2^2;$

对所述最优化目标函数进行求解，得到所述切向各异性系数Δ_T；

其中， $x_1=\frac{1}{2}\frac{\mathrm{\Δ}z}{\stackrel{\‾}{z}};$ $x_2=\frac{1}{2}\left\{\frac{\mathrm{\Δ}\mathrm{\α}}{\stackrel{\‾}{\mathrm{\α}}}-{\left(\frac{2\stackrel{\‾}{\mathrm{\β}}}{\stackrel{\‾}{\mathrm{\α}}}\right)}^2\frac{\mathrm{\Δ}G}{\stackrel{\‾}{G}}\right\};$ $x_3=\frac{1}{2}\frac{\mathrm{\Δ}\mathrm{\α}}{\stackrel{\‾}{\mathrm{\α}}};$ x₄＝Δ_T；x₅＝Δ_N；i和φ分别表示纵波的入射角和方位角；r为与入射角i和方位角φ有关的纵波反射系数；z＝ρα为纵波波阻抗；ρ为介质密度，g/cm³；α为纵波速度，m/s；G＝ρβ²为横波切向模量，β为横波速度；ε、δ和γ为Thomsen的各向异性参数；Δ_T为裂缝岩石物理参数中的切向各异性系数；Δ_N为裂缝岩石物理参数中的法向各异性系数；g＝μ/(λ+2μ)为横纵波速度比的平方，λ、μ为不含裂隙各向同性岩石的拉梅参数；

$A=\left[\begin{array}{ccccc}1& {\sin }^2{i}_1& W_1& \[-g(1+R_1)+2g\]P_1& P_1\[({\mathrm{gsin}}^2{\mathrm{\φ}}_1+g-1){\tan }^2{i}_1+2g-1\]g\\ 1& {\sin }^2{i}_2& W_2& \[-g(1+R_2)+2g\]P_2& P_2\[({\mathrm{gsin}}^2{\mathrm{\φ}}_2+g-1){\tan }^2{i}_2+2g-1\]g\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ 1& {\sin }^2{i}_n& W_n& \[-g(1+R_n)+2g\]P_n& P_n\[({\mathrm{gsin}}^2{\mathrm{\φ}}_n+g-1){\tan }^2{i}_n+2g-1\]g\end{array}\right],$

R_k、P_k及W_k为中间变量，R_k＝sin²φ_ktan²i_k，P_k＝cos²φ_ksin²i_k，W_k＝sin²i_ktan²i_k，r_k＝r_k(i_k,φ_k)为方位角道集k对应的纵波反射系数，b＝[r₁(i₁,φ₁)；r₂(i₂,φ₂)；…；r_n(i_n,φ_n)]^T，x＝[x₁,x₂,x₃,x₄,x₅]^T，k∈[1，n]，n为方位角道集的数量。

3.根据权利要求2所述的基于岩石物理的叠前裂缝定量预测方法，其特征在于，对所述最优化目标函数进行求解，得到所述切向各异性系数Δ_T，包括：

对所述最优化目标函数进行Tikhonov正则化处理，经过Tikhonov正则化处理后的所述最优化目标函数为：

$\min E{\left(x\right)}_{obj-reg}=\left|\right|Ax-b|{|}_2^2+\mathrm{\τ}|\left|X\right|{|}_2^2;$

计算经过Tikhonov正则化处理后的所述最优化目标函数的最小二乘解x＝(A^TA+τI)^-1A^Tb，即得所述切向各向异性系数Δ_T；

其中，τ为正则化参数。

4.根据权利要求1所述的基于岩石物理的叠前裂缝定量预测方法，其特征在于，根据所述切向各异性系数Δ_T计算待测储层的裂缝密度e，包括：

根据所述切向各异性系数Δ_T与裂缝密度e之间的关系，按照下式计算所述裂缝密度e：

$e=\frac{3(3-2g){\mathrm{\Δ}}_T}{16}.$

5.根据权利要求4所述的基于岩石物理的叠前裂缝定量预测方法，其特征在于，所述切向各异性系数Δ_T与裂缝密度e之间的关系是根据Hudson币状裂隙模型及Schenberg线性滑动模型确定的。

6.一种基于岩石物理的叠前裂缝定量预测系统，其特征在于，所述基于岩石物理的叠前裂缝定量预测系统包括：

方位角道集划分单元，用于根据待测储层的叠前CMP道集中的方位角及入射角的大小，将所述叠前CMP道集划分为多个方位角道集：方位角道集1，方位角道集2，…，方位角道集n；

纵波反射系数获取单元，用于对各所述方位角道集进行叠加、偏移处理，计算各方位角道集的纵波反射系数；

岩石物理参数计算单元，用于根据Ruger方程精确式以及各所述方位角道集及其纵波反射系数，计算待测储层中裂缝的切向各异性系数Δ_T；

裂缝密度计算单元，用于根据所述切向各异性系数Δ_T计算待测储层的裂缝密度e。

7.根据权利要求6所述的基于岩石物理的叠前裂缝定量预测系统，其特征在于，所述岩石物理参数计算单元包括：

Ruger方程预处理模块，用于根据HTI介质中各向异性参数与裂缝岩石物理参数之间的关系，对Ruger方程精确式进行变换，得到纵波反射系数简化方程：

x₁+x₂sin²i+x₃sin²itan²i+x₄[-g(1+sin²φtan²i)+2g]cos²φsin²i

+x₅gcos²φsin²i[(gsin²φ+g-1)tan²i+2g-1]＝r；

矩阵方程生成模块，用于将所述方位角道集及其对应的纵波反射系数代入所述纵波反射系数简化方程得到矩阵方程：

Ax＝b；

最优化目标函数建立模块，用于根据所述矩阵方程建立一用于反演所述x的最优化目标函数：

$\begin{array}{cc}\min & E{\left(x\right)}_{obj}=\left|\right|Ax-b|{|}_2^2\end{array};$

求解模块，用于对所述最优化目标函数进行求解，得到所述切向各异性系数Δ_T；

其中， $x_1=\frac{1}{2}\frac{\mathrm{\Δ}z}{\stackrel{\‾}{z}};$ $x_2=\frac{1}{2}\left\{\frac{\mathrm{\Δ}\mathrm{\α}}{\stackrel{\‾}{\mathrm{\α}}}-{\left(\frac{2\stackrel{\‾}{\mathrm{\β}}}{\stackrel{\‾}{\mathrm{\α}}}\right)}^2\frac{\mathrm{\Δ}G}{\stackrel{\‾}{G}}\right\};$ $x_3=\frac{1}{2}\frac{\mathrm{\Δ}\mathrm{\α}}{\stackrel{\‾}{\mathrm{\α}}};$ x₄＝Δ_T；x₅＝Δ_N；i和φ分别表示纵波的入射角和方位角；r为与入射角i和方位角φ有关的纵波反射系数；z＝ρα为纵波波阻抗；ρ为介质密度，g/cm³；α为纵波速度，m/s；G＝ρβ²为横波切向模量，β为横波速度；ε、δ和γ为Thomsen的各向异性参数；Δ_T为裂缝岩石物理参数中的切向各异性系数；Δ_N为裂缝岩石物理参数中的法向各异性系数；g＝μ/(λ+2μ)为横纵波速度比的平方，λ、μ为不含裂隙各向同性岩石的拉梅参数；

$A=\left[\begin{array}{ccccc}1& {\sin }^2{i}_1& W_1& \[-g(1+R_1)+2g\]P_1& P_1\[({\mathrm{gsin}}^2{\mathrm{\φ}}_1+g-1){\tan }^2{i}_1+2g-1\]g\\ 1& {\sin }^2{i}_2& W_2& \[-g(1+R_2)+2g\]P_2& P_2\[({\mathrm{gsin}}^2{\mathrm{\φ}}_2+g-1){\tan }^2{i}_2+2g-1\]g\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ 1& {\sin }^2{i}_n& W_n& \[-g(1+R_n)+2g\]P_n& P_n\[({\mathrm{gsin}}^2{\mathrm{\φ}}_n+g-1){\tan }^2{i}_n+2g-1\]g\end{array}\right],$

R_k、P_k及W_k为中间变量，R_k＝sin²φ_ktan²i_k，P_k＝cos²φ_ksin²i_k，W_k＝sin²i_ktan²i_k，r_k＝r_k(i_k,φ_k)为方位角道集k对应的纵波反射系数，b＝[r₁(i₁,φ₁)；r₂(i₂,φ₂)；…；r_n(i_n,φ_n)]^T，x＝[x₁,x₂,x₃,x₄,x₅]^T，k∈[1，n]，n为方位角道集的数量。

8.根据权利要求7所述的基于岩石物理的叠前裂缝定量预测系统，其特征在于，所述求解模块具体用于：

对所述最优化目标函数进行Tikhonov正则化处理，经过Tikhonov正则化处理后的所述最优化目标函数为：

$\min E{\left(x\right)}_{obj-reg}=\left|\right|Ax-b|{|}_2^2+\mathrm{\τ}|\left|x\right|{|}_2^2;$

计算经过Tikhonov正则化处理后的所述最优化目标函数的最小二乘解x＝(A^TA+τI)^-1A^Tb，即得所述切向各向异性系数Δ_T；

其中，τ为正则化参数。

9.根据权利要求6所述的基于岩石物理的叠前裂缝定量预测系统，其特征在于，所述裂缝密度计算单元具体用于：

根据所述切向各异性系数Δ_T与裂缝密度e之间的关系，按照下式计算所述裂缝密度e：

$e=\frac{3(3-2g){\mathrm{\Δ}}_T}{16}.$

10.根据权利要求9所述的基于岩石物理的叠前裂缝定量预测系统，其特征在于，所述切向各异性系数Δ_T与裂缝密度e之间的关系是根据Hudson币状裂隙模型及Schenberg线性滑动模型确定的。