一种各向异性流体识别因子反演方法及系统与流程

本发明涉及地球物理反演领域，尤其涉及一种各向异性流体识别因子反演方法及系统。
背景技术：
：随着勘探程度的加深，在裂缝型油气藏评价中要求对裂缝含流体情况进行评价。目前由于裂缝各向异性参数比较难以获得，因此大多数研究是基于Ruger进行简化，得到流体识别因子的近似解。然而采用Ruger进行简化的方法，在得到的流体识别因子的近似解中容易引入较多误差，并且基于叠前反演的原理及叠前资料的直接影响，使的反演结果的准确性和可信度大幅降低。技术实现要素：本发明的实施例提供一种各向异性流体识别因子反演方法及系统，能够以叠后方位数据为基础，提高了反演结果的准确性和可信度。为达到上述目的，本发明的实施例采用如下技术方案：第一方面，本发明提供了一种各向异性流体识别因子反演方法，包括：选取三个不同的方位角，分别对选取的方位角的叠前数据按照入射角方向进行叠加，并分别将每个方位角对应的叠后地震数据体进行反演，得到与每个方位角对应的反射系数；对叠前方位角道集进行分析得到入射角信息，并根据入射角信息计算出中间参数；根据选取的三个不同的方位角、反演得到的反射系数和计算出的中间参数，计算得到各向异性参数和变量；根据各向异性参数和变量进行计算得到各向异性流体识别因子。较佳的，所述反演为稀疏脉冲反演。较佳的，所述入射角信息包括最大入射角、最小入射角和入射角角度的间隔。较佳的，每个选取的所述方位角的角度均大于30度。第二方面，本发明提供了一种各向异性流体识别因子反演系统，包括选取叠加模块、分析模块、第一计算模块和第二计算模块；所述选取叠加模块，用于选取三个不同的方位角，分别对选取的方位角的叠前数据分别按照入射角方向进行叠加，并分别将每个方位角对应的叠后地震数据体进行反演，得到与每个方位角对应的反射系数；所述分析模块，用于对叠前方位角道集进行分析得到入射角信息，并根据入射角信息计算出中间参数；所述第一计算模块，用于根据选取的三个不同的方位角、反演得到的反射系数和计算出的中间参数，计算得到各向异性参数和变量；所述第二计算模块，用于根据各向异性参数和变量进行计算得到各向异性流体识别因子。较佳的，所述反演为稀疏脉冲反演。较佳的，所述入射角信息包括最大入射角、最小入射角和入射角角度的间隔。较佳的，每个选取的所述方位角的角度均大于30度。附图说明为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。图1为本发明实施例提供的一种各向异性流体识别因子反演方法的流程图；图2为本发明实施例提供的一种各向异性流体识别因子反演系统的结构图；图3为本发明实施例提供的正演得到的60°、75°、90°方位的叠前数据体示意图；图4为本发明实施例提供的反演得到的60°、75°、90°方位对应的反射系数示意图。具体实施方式下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。本发明实施例提供了一种各向异性流体识别因子反演方法，如图1所示，所述方法包括：101、选取三个不同的方位角，分别对选取的方位角的叠前数据按照入射角方向进行叠加，并分别将每个方位角对应的叠后地震数据体进行反演，得到与每个方位角对应的反射系数。所述步骤101在本发明实施例中可具体实现为，选取任意三个不同方位角并对的叠前数据，分别按照入射角方向进行叠加，然后对叠后方位地震数据体分别进行反演，所述反演为稀疏脉冲反演。在本发明实施例中，对叠后方位地震数据体进行稀疏脉冲反演分为以下三步：(1)反射系数反演；采用最大似然反褶积进行反射系数的反演，最大似然反褶积对地层的假设认为：地层的反射系数是由较大的反射界面的反射和具有高斯背景的小反射叠加组合而成，导出一个最小目标函数J：J=1R2ΣK=1Lr2(K)+1N2ΣK=1Ln2(K)-2Mln(λ)-2(L-M)ln(1-λ)]]>式中，R2和N2分别为反射系数和噪音的均方值，r(K)和n(K)表示第K个采样点的反射系数和噪音，M表示反射层数，L表示采样总数，λ表示给定反射系数的似然值。通过多次迭代，求取反射系数。(2)根据反射系数的反演结果结合阻抗趋势计算一个初始的波阻抗；根据最大似然反褶积计算得到的反射系数，结合初始阻抗模型，采用递推算法，反演得到初始的波阻抗模型：Z(i)=Z(i-1)1+R(i)1-R(i)]]>式中，Z(i)为第i层的波阻抗值，R(i)为第i层的反射系数。(3)结合井的约束条件进行波阻抗反演；约束稀疏脉冲反演对每一道依据目标函数对计算出的初始波阻抗进行调整，包括对反射系数的调整。目标优化函数为：F＝Lp(r)+λLq(s-d)+α-1L1ΔZ式中，r为反射系数序列，Δz为与阻抗趋势的差序列，d为地震道序列，s为合成地震道序列，λ为残差权重因子，α为趋势权重因子，p、q为L模因子。具体的，右式第一项反映了反射系数的绝对值和，第二项反映了合成声波记录与原始地震数据的差值，第三项为趋势约束项。基于以上三步，得到了调整后的与每个方位角的反射系数需要说明的是，每个选取的所述方位角的角度均大于30度，选择较大方位角的叠前数据，进行计算可以减小简化误差，如果选择的方位角太小，比如选取方位角为0,15,30的效果就会不好。102、对叠前方位角道集进行分析得到入射角信息，并根据入射角信息计算出中间参数。所述步骤102在本发明实施例中可具体实现为，所述入射角信息包括最大入射角、最小入射角和入射角角度的间隔。在本发明实施例中，首先采用Ruger(1998)借助弱各向异性的概念，推导出的HTI介质的纵波近似反射系数关系式：式中:和分别为HTI介质上下两层的纵波、横波速度平均值以及密 度平均值；Δδ(V)，Δε(V)和Δγ分别为上下两层各向异性参数差值，与裂缝密度和缝隙充填流体有关；θ和分别代表入射角和方位角，与裂缝的发育方向有关。然后，将HTI介质的纵波近似反射系数关系式按照入射角方向叠加消掉θ，得到：其中A=Σθ1θ212sin2θ,B=Σθ1θ212sin2θtan2θ,D=Σθ1θ2sin2θ,]]>C=Σθ1θ212sin2θtan2θ,E=12(Δρρ‾+Δαα‾)+12{Δαα‾-(2β‾α‾)2(Δρρ‾+2Δββ‾)}sin2θ+12Δαα‾sin2θtan2θ,]]>0°≤θ1<θ2≤30°。最后，根据A=Σθ1θ212sin2θ,B=Σθ1θ212sin2θtan2θ,D=Σθ1θ2sin2θ,]]>0°≤θ1<θ2≤30°计算出中间参数A、B和D。103、根据选取的三个不同的方位角、反演得到的反射系数和计算出的中间参数，计算得到各向异性参数和变量；所述步骤103在本发明实施例中可具体实现为，此时，反射系数分别用中间参数A、B和D表示为：解反射系数的方程组，得到：另外，根据Bakulin(2000)推导出Thomsen各向异性参数与ΔN，ΔT之间的关系：ε(V)＝-2g(1-g)ΔNδ(V)＝-2g[(1-g)ΔN+ΔT]γ＝-2ΔT/2以及，根据Schoenberg和Sayers(1995)根据线性滑动理论，推导出裂缝流体识别因子：Fluidindicator≈gΔNΔτ]]>其中，g=(βα)2.]]>根据Thomsen各向异性参数与ΔN，ΔT之间的关系和裂缝流体识别因子，可以得到：Fluidindicator≈gΔNΔτ=ϵ(V)-δ(V)4gγ+1]]>在Ruger推导出的HTI介质的纵波近似反射系数关系式中，当方位角较大时近似为0，因此解反射系数的方程组得到的和Fluidindicator变为以下形式：F1uidindicator≈gΔNΔT=1-δ(V)4gγ≈1δ(V)(2β‾α‾)2Δγ]]>将方位角反射系数和中间参数A、B、D带入最后得到的Δδ(V)、两个方程中，计算得到各向异性参数Δδ(V)，变量104、根据各向异性参数和变量进行计算得到各向异性流体识别因子。所述步骤103在本发明实施例中可具体实现为，将各向异性参数Δδ(V)，变量带入到Fluidindicator中，得到各向异性流体识别因子，所述各向异性流体识别因子，可以用于分析裂缝内充填流体性质，在裂缝预测中起着至关重要的作用。在本发明实施例的另一种实现方式中，如图2所示，本发明提供了一种各向异性流体识别因子反演系统，包括选取叠加模块201、分析模块202、第一计算模块203和第二计算模块204。所述选取叠加模块201，用于选取三个不同的方位角，分别对选取的方位角的叠前数据分别按照入射角方向进行叠加，并分别将每个方位角对应的叠后地震数据体进行反演，得到与每个方位角对应的反射系数。所述选取模块201在本发明实施例中可具体实现为，选取任意三个不同方位角并对的叠前数据，分别按照入射角方向进行叠加，然后对叠后方位地震数据体分别进行反演，所述反演为稀疏脉冲反演。在本发明实施例中，对叠后方位地震数据体进行稀疏脉冲反演分为以下三步：(1)反射系数反演；采用最大似然反褶积进行反射系数的反演，最大似然反褶积对地层的假设认为：地层的反射系数是由较大的反射界面的反射和具有高斯背景的小反射叠加组合而成，导出一个最小目标函数J：J=1R2ΣK=1Lr2(K)+1N2Σk=1Ln2(K)-2Mln(λ)-2(L-M)ln(1-λ)]]>式中，R2和N2分别为反射系数和噪音的均方值，r(K)和n(K)表示第K个采样点的反射系数和噪音，M表示反射层数，L表示采样总数，λ表示给定反射系数的似然值。通过多次迭代，求取反射系数。(2)根据反射系数的反演结果结合阻抗趋势计算一个初始的波阻抗；根据最大似然反褶积计算得到的反射系数，结合初始阻抗模型，采用递推算法，反演得到初始的波阻抗模型：Z(i)=Z(i-1)1+R(i)1-R(i)]]>式中，Z(i)为第i层的波阻抗值，R(i)为第i层的反射系数。(3)结合井的约束条件进行波阻抗反演；约束稀疏脉冲反演对每一道依据目标函数对计算出的初始波阻抗进行调整，包括对反射系数的调整。目标优化函数为：F＝Lp(r)+λLq(s-d)+α-1L1ΔZ式中，r为反射系数序列，Δz为与阻抗趋势的差序列，d为地震道序列，s为合成地震道序列，λ为残差权重因子，α为趋势权重因子，p、q为L模因子。具体的，右式第一项反映了反射系数的绝对值和，第二项反映了合成声波记录与原始地震数据的差值，第三项为趋势约束项。基于以上三步，得到了调整后的与每个方位角的反射系数需要说明的是，每个选取的所述方位角的角度均大于30度，选择较大方位角的叠前数据，进行计算可以减小简化误差，如果选择的方位角太小，比如选取方位角为0,15,30的效果就会不好。所述分析模块202，用于对叠前方位角道集进行分析得到入射角信息，并根 据入射角信息计算出中间参数。所述分析模块202在本发明实施例中可具体实现为，所述入射角信息包括最大入射角、最小入射角和入射角角度的间隔。在本发明实施例中，首先采用Ruger(1998)借助弱各向异性的概念，推导出的HTI介质的纵波近似反射系数关系式：式中:和分别为HTI介质上下两层的纵波、横波速度平均值以及密度平均值；Δδ(V)，Δε(V)和Δγ分别为上下两层各向异性参数差值，与裂缝密度和缝隙充填流体有关；θ和分别代表入射角和方位角，与裂缝的发育方向有关。然后，将HTI介质的纵波近似反射系数关系式按照入射角方向叠加消掉θ，得到：其中A=Σθ1θ212sin2θ,B=Σθ1θ212sin2θtan2θ,D=Σθ1θ2sin2θ,]]>C=Σθ1θ212sin2θtan2θ,E=12(Δρρ‾+Δαα‾)+12{Δαα‾-(2β‾α‾)2(Δρρ‾+2Δββ‾)}sin2θ+12Δαα‾sin2θtan2θ,]]>0°≤θ1<θ2≤30°。最后，根据A=Σθ1θ212sin2θ,B=Σθ1θ212sin2θtan2θ,D=Σθ1θ2sin2θ,]]>0°≤θ1<θ2≤30°计算出中间参数A、B和D。所述第一计算模块203，用于根据选取的三个不同的方位角、反演得到的反射系数和计算出的中间参数，计算得到各向异性参数和变量；所述第一计算模块203在本发明实施例中可具体实现为，此时，反射系数分别用中间参数A、B和D表示为：解反射系数的方程组，得到：另外，根据Bakulin(2000)推导出Thomsen各向异性参数与ΔN，ΔT之间的关系：ε(V)＝-2g(1-g)ΔNδ(V)＝-2g[(1-g)ΔN+ΔT]γ＝-2ΔT/2以及，根据Schoenberg和Sayers(1995)根据线性滑动理论，推导出裂缝流体识别因子：Fluidindicator≈gΔNΔτ]]>其中，g=(βα)2.]]>根据Thomsen各向异性参数与ΔN，ΔT之间的关系和裂缝流体识别因子，可以得到：F1uidindicator≈gΔNΔτ=ϵ(V)-δ(V)4gγ+1]]>在Ruger推导出的HTI介质的纵波近似反射系数关系式中，当方位角较大时近似为0，因此解反射系数的方程组得到的Δε(V)-Δδ(V)、和Fluidindicator变为以下形式：F1uidindicator≈gΔNΔT=1-δ(V)4gγ≈1δ(V)(2β‾α‾)2Δγ]]>将方位角反射系数和中间参数A、B、D带入最后得到的Δδ(V)、两个方程中，计算得到各向异性参数Δδ(V)，变量所述第二计算模块204，用于根据各向异性参数和变量进行计算得到各向异性流体识别因子。所述第二计算模块204在本发明实施例中可具体实现为，将各向异性参数Δδ(V)，变量带入到Fluidindicator中，得到各向异性流体识别因子，所述各向异性流体识别因子，可以用于分析裂缝内充填流体性质，在裂缝预测中起着至关重要的作用。在本发明实施例的另一种实现方式中，基于表1各向异性介质模型数据进行试算。分层αβρε(v)δ(v)γ124009122.170002280010642.250.040.180.12表1首先对60°，75°，90°方位的叠前数据(如图3所示)分别按照入射角方向进行叠加，然后对叠后方位地震数据体进行稀疏脉冲反演，选择较大方位角的叠前数据进行计算可以减小简化Ruger公式所带来的误差，计算得到对应的反射系数体RPP(60°)，RPP(75°)，RPP(90°)，如图4所示。然后分析叠前方位道集的入射角信息，本次叠前地震数据沿入射角叠加的范围选取为θ1＝0°，θ2＝30°，入射角间隔为1°，因此根据分别计算出A＝1.3607，B＝0.2636，D＝2.7214。在实际资料应用中入射角叠加的范围θ1θ2的选取可以根据叠前方位道集的资料品质进行调整，比如目前大多数叠前地震道集中存在着小入射角和大入射角能量分布不均衡的现象，可以选取θ1＝15°，θ2＝25°来避开资料品质低的道集，选取资料品质高并且各向异性强的道集进行计算。接下来将已知的方位角信息反射系数体RPP(60°)，RPP(75°)，RPP(90°)和A、B、D的值，带入公式：和中，计算得到各向异性参数Δδ(V)，变量通过方程：Fluidindicator≈gΔNΔT=1-δ(V)4gγ≈1δ(V)(2β‾α‾)2Δγ]]>计算得到各向异性流体识别因子，可以用于分析裂缝内充填流体性质，在裂缝预测中起着至关重要的作用。通过以上的实施方式的描述，所属领域的技术人员可以清楚地了解到本发明可借助软件加必需的通用硬件的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在可读取的存储介质中，如计算机的软盘，硬盘或光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本
技术领域：
的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。当前第1页1 2 3