基于广义极值分布的叠前非高斯AVO反演方法

技术领域

本发明属于地震数据反演技术领域，具体涉及一种基于广义极值分布的叠 前非高斯AVO反演方法。

背景技术

在油气储层预测、油藏描述工作中，我们常常使用纵横波速度、体积密度、 泊松比、杨氏模量、体积模量、剪切模量等弹性参数来描述地下介质的岩性特 征。叠前地震反演可以在钻井资料较少的情况下，较为准确的反演出目标层位 的弹性参数。通过添加约束条件等手段我们可以进一步的提高叠前地震反演的 精度，有利于对页岩气储层的精细刻画，对促进页岩气低成本的开发利用具有 重要的意义。AVO(Amplitude variation withoffset,振幅随偏移距的变化)是叠前 反演最常用的方法之一,其基本原理是通过反射振幅与炮间距的函数关系来推 算地层的弹性参数(最为常见的是纵横波速度和体积密度)。传统的反演方法 为了计算简便，通常在两个假设的前提下简化建模过程和进行反演实施：(1) 假设地震数据中的噪声服从高斯分布、柯西分布、T分布等；(2)假设地层反射 系数为白噪。然后现有诸多研究成果表明全球不同地区的大量测井资料显示, 反射系数也不满足高斯分布，其叠前地震数据中包含的噪声也不服从高斯分布。 在反演过程中,如果不考虑这些非高斯性,反演的结果必将和地下实际情况存在 较大的差异。

相比于常规油气藏描述，需要实施水力压裂技术才能满足开发要求的页岩 油气藏描述的要求更高。理论分析和实际资料分析结果表明，常规的叠前地震 反演结果难以满足页岩气“甜点区”精细刻画的需求，其主要原因在于常规的 叠前地震反演均是基于地震数据中包含噪声满足高斯分布的假设，以及未考虑 待反演弹性参数的分布特征。

地下的岩性及岩石孔隙中所含的流体性质的改变都会对地震反射波振幅 产生一定的影响。AVO正是利用该理论进行岩性识别和油气预测。AVO理论基 于射线理论的Zoeppritz方程建立了反射系数、入射角以及弹性参数之间的关 系。后来很多地球物理学家对Zoeppritz方程进行各种不同方式的简化，其中 Shuey对泊松比与反射系数的关系进行了研究，指出泊松比的变化会对反射系 数随入射角变化产生巨大的影响，当入射角小于30°时，所归纳出的近似方程 为 R(θ)≈A+Bsin²θ

其中θ是入射 角度，R(θ)表示PP波反射系数，数值随角度θ变化，A、B分别为截距和梯 度属性，Δα、Δρ、Δσ分别为界面两侧纵横波速度、密度及泊松比的相对变化 量(视泊松比差反射率)，分别代表反射界面两侧介质的纵波速度、 密度及泊松比的平均值，在方程中参数A代表纵波垂直入射时的反射系数，D 为截距项，参数B为梯度项，假定地层σ的背景值为(纵波速度是横波速 度的2倍)，代入梯度项B进行改写得到：

进一步化简可得：

上式表明，当界面两侧介质的纵波阻抗为定值时，视泊松比差反射率Δσ会 对反射振幅随入射角的变化产生较大的影响，在砂岩含气的情况下，其泊松比 会明显减小，引起两侧介质的视泊松比差反射率Δσ增大，从而导致AVO出现 异常。上式表明Δσ变化可反映在AVO截距和梯度属性上，可通过提取截距和 梯度属性进行求和得到视泊松比差反射率。

传统的AVO反演方法为了简化建模过程便于求解，常常假设噪声是高斯白 噪的。叠后地震数据需要由同一地层反射点的多个地震道数据水平叠加。与同 一地层反射点的多个地震道构成的叠前地震道集数据相比，叠后数据的信噪比 更高，其噪声分布基本满足高斯白噪分布。对于叠前地震资料而言，特别在那 些构造很复杂和目的层为非常规储层地区采集的地震数据，噪声是高斯白噪的 假设是不成立的。近年来，越来越多的地球物理学者注意到了地震资料的统计 分布具有非高斯性。如果在地震反演过程中假设其是服从高斯分布的，那么获 得的结果将与地下实际情况存在较大的差异。尽管学者们发现使用高斯分布假 设存在的问题，但是使用非高斯分布去拟合实际数据在实际应用中仍然存在着 许多的问题：首先，需要选择恰当的分布去表征地震数据就较为困难，如何求 取分布的参数也可能会给地震反演带来很大的计算量；其次，在非高斯分布的 假设下，代价函数中需要使用除了l₂范数以外的其他范数，例如l₀、l₄范数等，>

Fisher和Tippett(1928)在研究独立分布的极大值渐进分布时，进一步的提 出了3种新的极值分布：Frechet、Weibull和Gumbel分布。Jekkinson又将这 三种分布统一成一个极值分布——广义极值分布(GEV)，该极值分布由3个 参数来控制分布的形态、位置和尖锐程度，但是该理论存在一个前提条件，当 极值的渐进分布存在且为非退化时该理论才成立。

Kennedy和Eberhart(1995年)在Reynolds(1987年)提出的鸟群集聚 模型的启发下，提出了粒子群算法(Particle Swarm Optimization,PSO)。 PSO算法修正了鸟群集聚模型，使得粒子群能够向最优解移动。粒子群算法(PSO) 是一种模拟鸟群捕食行为的智能算法，是目前求解最优值最为常用的方法之一。 粒子群算法通过每个粒子的信息与群体的信息协同寻找最优值。由于在计算过 程中并不需要计算目标函数梯度等信息，加之该方法需要调整的参数较少，因 此本发明选择此方法求解方程组。

我们使用粒子i去表示一只鸟，每个粒子具有两个信息：位置矢量Xi和速 度矢量v_i，通过以x_i为自变量的适应度函数来评价粒子i在x_i处的好坏。每个>

v_i＝v_i+c₁×rand()×(Pbest_i-x_i)+c₂×rand()×(Gbest_i-x_i)

x_i＝x_i+v_i

其中，i＝1,2,…,N，N是此群中粒子的总数；v_i是粒子的速度；rand()是介于(0,1)之间的随机数；x_i是粒子的当前位置；c₁和c₂是学习因子，通常c₁＝c₂＝2；v_i的>max(大于0)，如果v_i大于v_max，则v_i＝v_max。

上式为PSO的标准形式。公式的第一部分称为“记忆项”，表示上次速度大 小和方向的影响；公式的第二部分称为“自身认知项”，是从当前点指向粒子 自身最好点的一个矢量，表示粒子的动作来源于自己经验的部分；公式的第三 部分称为“群体认知项”，是一个从当前点指向种群最好点的矢量，反映了粒 子间的协同合作和知识共享。粒子就是通过自己的经验和同伴中最好的经验来 决定下一步的运动。

v_i＝ω×v_i+c₁×rand()×(pbest_i-x_i)+c₂×rand()×(gbest_i-x_i)

其中，ω为惯性因子，其值为非负值。当ω较大时，全局寻优能力强，局 部寻优能力弱；当ω较小时，全局寻优能力弱，局部寻优能力强。引入了惯性 因子ω，大大地提高了PSO算法的性能。针对不同的最优化问题，我们可以对 惯性因子ω进行调整，进而对全局寻优和局部寻优能力做出调整，能较大地提 高收敛速度。

以上这些基础粒子群算法仅仅能对单峰且非病态的目标函数求出最优解， 对于多峰或者病态函数，无论增加迭代次数还是种群数目，都很难取得满足精 度要求的解。为了进一步优化PSO结果，刘霞等(2007年)提出一种抑制局部 最优的粒子群算法(RPSO)。该算法引入了变异算子，在PSO算法陷入局部最 优时，给粒子的某些维变量添加一个扰动，使粒子能够跳出局部最优值，到达 全局最优值。但是，当目标函数很复杂时，使用该算法，仍然很难跳出局部最 优点。

发明内容

为了解决现有的基于贝叶斯框架的AVO反演方法中仅仅考虑了地震数据 中噪声的先验分布，未考虑反射系数的先验分布等问题，本发明采用广义极值 分布来表征反射系数的非高斯特性，并基于贝叶斯反演框架，提出了一种基于 广义极值分布的叠前非高斯AVO反演方法，较基于高斯假设的方法的精度有 了明显的提高，获得的结果能够大幅度地提高油气勘探和开发过程中钻井的成 功率。

为实现上述目的，本发明提供一种基于广义极值分布的叠前非高斯AVO 反演方法，包括以下步骤：

S1、根据页岩气钻井、测井数据和叠前地震数据分析页岩储层弹性参数、 叠前地震数据中噪声的分布特征；

S2、在基于贝叶斯理论的AVO反演框架下，构造包含反射系数统计分布 约束的AVO反演目标函数；

S3、采用具有参数自适应调整特征的广义极值分布来描述反射系数的统计 分布；

S4、采用改进的粒子群算法对AVO反演目标函数进行求解，得到AVO反 演结果。

进一步地，所述步骤S2中，将含噪声的地震数据表示为：

d＝Gr+n

其中，d为野外观测地震数据，G为地震子波，r为地层反射系数，n为地 震数据观测噪声。

进一步地，所述步骤S2中，包含反射系数统计分布约束的AVO反演目标 函数表示为：

其中，ε、μ、δ为广义极值分布的3个参数，r_i为测井数>n为噪声的方差。

进一步地，所述步骤S3中，广义极值分布的概率密度函数表示为：

进一步地，所述步骤S3中，假设随机变量x₁,x₂,…x_n为广义极值分布的一>

进一步地，所述步骤S3中，采用极大似然法计算广义极值分布的参数， 表示为：

进一步地，所述步骤S4具体包括以下分步骤：

S41、初始化粒子群大小，初始位置x₀，初始速度v₀，初始温度T，截止>min，衰减因子α；

S42、计算粒子群中各粒子的适应值

S43、判断粒子的适应值是否小于设定阈值ep；若是，则操作结束； 若否，则进行步骤S44；

S44、判断粒子的适应值是否大于粒子下一次迭代后的适应值若是，则进行步骤S47；若否，则进行步骤S45；

S45、判断初始温度T是否大于截止温度T_min；若是，则进行步骤S46；若>

S46、判断若是，则更新当前温度，并进行 步骤S47；若否，则进行步骤S48；

S47、更新粒子在迭代过程中的最优解Pbest_i；

S48、更新所有粒子的全局最优值Gbest_i；

S49、更新下一迭代过程的粒子位置和粒子速度并返回步骤S42。

本发明的有益效果是：本发明通过分析页岩储层弹性参数、叠前地震数据 中噪声的分布特征，在噪声分布特征分析的基础上，基于贝叶斯理论框架，分 别构建考虑叠前地震数据包含噪声分布特征和弹性参数分布特征的叠前地震 反演目标函数，以及综合考虑两种分布的叠前地震反演目标函数，最后结合最 优化反演理论和方法，实现用于实际页岩储层精细刻画的弹性参数反演的目的， 进而指导利用叠前地震资料识别页岩储层甜点区，以及为压裂方案提供科学依 据。

附图说明

图1是本发明的基于广义极值分布的叠前非高斯AVO反演方法流程图；

图2是本发明中反射系数概率密度函数示意图；

图3是本发明中纵波反射率非高斯性图；

图4是本发明中横波反射率非高斯性图；

图5是本发明中密度反射率非高斯性图；

图6是本发明中改进的粒子群算法流程图；

图7是本发明中密度反射率、纵波反射率和横波反射率的反演结果对比图；

图8是本发明中密度、纵波速度、横波速度的反演结果对比图；

图9是本发明中合成地震记录对比图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实 施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅 用以解释本发明，并不用于限定本发明。

如图1所示，一种基于广义极值分布的叠前非高斯AVO反演方法，包括 以下步骤：

S1、根据页岩气钻井、测井数据和叠前地震数据分析页岩储层弹性参数、 叠前地震数据中噪声的分布特征；

S2、在基于贝叶斯理论的AVO反演框架下，构造包含反射系数统计分布 约束的AVO反演目标函数；

S3、采用具有参数自适应调整特征的广义极值分布来描述反射系数的统计 分布；

S4、采用改进的粒子群算法对AVO反演目标函数进行求解，得到AVO反 演结果。

在上述步骤S1中，本发明以实际页岩气钻井、测井资料(包含纵横波速 度、密度)、叠前地震资料为依托，以叠前AVO方法理论，以及AVO近似 方程等为理论基础，借鉴统计分析、信号分析、地震波场正演等若干手段，分 析页岩储层弹性参数、叠前地震数据中噪声的分布特征。

在上述步骤S2中，本发明将含噪声的地震数据表示为 d＝Gr+n

其中，d为野外观测地震数据；G为地震子波；r为地层反射系数；n为地 震数据观测噪声。

贝叶斯理论表示为：

P(r|d)∝P(r)P(d|r)

其中，P(r|d)为对于给定的野外地震数据d，地层反射系数r(模型)正 确的概率，P(r)为地层反射系数的似然函数；P(d|r)为似然函数。

假设噪声是高斯白噪的，方差为则噪声n可表示为：

噪声的似然函数可表示为

根据上面两式可得：

由于似然函数可以用噪声表示，可得

理论上，广义极值分布可以表征任意分布的序列，故使用广义极值分布来 表征反射系数。反射系数的概率密度函数可表示为

其似然函数可表示为：

对上式两边取对数，可得：

对上式两边同时取对数

lnP(r|d)∝lnP(r)+lnP(d|r)

令

K₂＝-Nlnδ+lnK₁

得

将包含反射系数统计分布约束的AVO反演目标函数表示为

其中，ε、μ、δ为广义极值分布的3个参数，r_i为测井数据所得的反射>n为噪声的方差

对于给定的后验概率P(r|d)，通过寻找使后验概率最大的r计算反射系数r； 求P(r|d)得最大值可转化为求-ln(P(r|d))的最小值，忽略常数项，等价于求J的 最小值。

在上述步骤S3中，本发明分析工区内多口井的地层反射系数的统计分布 特征，分析发现地层反射系数呈现明显的厚拖尾特性，使用高斯分布是无法拟 合这种特性的，而广义极值分布可以很好地表征这种非高斯特性。图2为其中 4口井的统计特征，其中红色曲线为使用广义极值分布拟合的概率密度函数， 蓝色曲线为使用高斯分布拟合的概率密度函数。图2中的4个子图显示广义极 值分布可以更好的拟合非高斯分布的厚拖尾特征。

进一步对工区的一口井的纵波反射率、横波反射率和密度反射率的数据进 行统计分析其非高斯性，每个类型的数据都选取中间3000个数据点。

对每口井的非高斯性分析主要分为四个方面：

1)偏度：统计数据分布偏斜方向和程度的度量，是统计数据分布非对称 程度的数字特征。

2)峰度：是研究数据分布陡峭或平滑的统计量，通过对峰度系数的测量， 我们能够判定数据相对于正态分布而言是更陡峭还是平缓。

3)直方图：是表示统计数据变化情况的一种主要工具。用直方图比较直观 地看出数据的分布状态，对于数据分布状况一目了然，便于判断其分布情况。

4)正态概率图：正态概率图用于检查一组数据是否服从正态分布。是实数 与正态分布数据之间函数关系的散点图。如果这组实数服从正态分布，正态概 率图将是一条直线。

图3为纵波反射率γ(n)＝0.293197442084245，χ(n)＝7.574218345262363； 图4为横波反射率γ(n)＝0.495611893417731，χ(n)＝9.427969959545656；图5 为密度反射率γ(n)＝0.081983551852173，χ(n)＝9.829847271142070；其中左图 为数据的直方图，右图为数据的正态概率图，γ(n)为偏度，χ(n)为峰度。从图 3，4和5中，我们可以观察、认识到测井数据的纵波反射率、横波反射率和密 度的偏度均不为0，这表明三个测井数据均不对称；三个测井数据的峰度显然 均与0存在着较大的差距，这表明这些数据均不满足高斯分布；从右图中，这 些数据的正态概率图均偏离了对角线，这进一步证明了测井数据不满足高斯分 布。

通过测井数据可以计算反射系数，假定该反射系数可以代表整个工区的反 射系数的特征。使用广义极值分布去表征反射系数的统计特征是该方法优于传 统的反演方法的重要原因，因此如何尽可能准确的计算广义极值分布的参数是 提高反演精度的关键。

本发明通过调整广义极值分布的三个参数，拟合反射系数和噪声的非高斯 性，从而更加准确地描述地震数据、弹性参数等的非高斯特性。选择不同的广 义极值分布的参数，可以拟合任意分布特征的数据，能够避免使用单一分布的 不足。广义极值分布的概率密度函数表示为：

本发明采用极大似然法计算广义极值分布的参数，从而可以很好地表现样 本的统计特性。

假设随机变量x₁,x₂,…x_n为广义极值分布的一个随机样本，当ε≠0时，其>

两边同时取对数，可以得到对数似然函数：兙俥

式1-6中要求

等式两边同时对ε，μ和δ取偏导，并令其等于零，可得：兙

设

在上述步骤S4中，粒子群算法(Particle Swarm Optimization,PSO)通过 每个粒子的信息与群体的信息协同寻找最优值。由于在计算过程中并不需要计 算目标函数梯度等信息，加之该方法需要调整的参数较少，因此本发明选择此 方法求解方程组。

假定粒子i去表示一只鸟，每个粒子具有两个信息：位置矢量Xi和速度 矢量v_i，通过以x_i为自变量的适应度函数来评价粒子i在x_i处的好坏。每个粒>

由于基于广义极值分布的叠前非高斯AVO反演的目标函数极为复杂，粒 子群算法很容易陷入局部最优值，单纯的调整参数也很难提高算法性能，因此 考虑将模拟退火算法引入其中。对改进的PSO-SA算法进行多次实验，结果显 示在每次更新粒子的最优值时模拟退火的思想能最大程度上的提高算法性能。

如图6所示，上述步骤S4具体包括以下分步骤：

S41、初始化粒子群大小N，初始位置x₀，初始速度v₀，初始温度T，截>min，衰减因子α；

S42、计算粒子群中各粒子的适应值

S43、判断粒子的适应值是否小于设定阈值ep；若是，则操作结束； 若否，则进行步骤S44；

S44、判断粒子的适应值是否大于粒子下一次迭代后的适应值若是，则进行步骤S47；若否，则进行步骤S45；

S45、判断初始温度T是否大于截止温度T_min；若是，则进行步骤S46；若>

S46、判断若是，则更新当前温度，并进行 步骤S47；若否，则进行步骤S48；

S47、更新粒子在迭代过程中的最优解Pbest_i；

S48、更新所有粒子的全局最优值Gbest_i；

S49、更新下一迭代过程的粒子位置和粒子速度并返回步骤S42。

模拟退火算法具有能跳出局部最优解能力，在计算每个粒子的最优值时， 当前粒子具有一定的概率接受一个劣解，增强了粒子的多样性，增强了粒子收 敛到全局最优值的能力。

本发明从井数据中截取采样率为2ms的200个点，分别使用本发明方法和 常规方法执行参数反演，反演结果。

图7分别为密度反射率、纵波反射率和横波反射率的反演结果对比，图8 分别为密度、纵波速度、横波速度的反演结果对比。每个子图中，蓝色曲线为 测井数据，黑色曲线为使用广义极值分布拟合反射系数的反演结果，红色数据 为传统的使用高斯分布拟合反射系数的反演结果。从图7和图8中可得，两者 均能够表现地层的弹性参数变化趋势，但是基于高斯假设的反演结果(红色曲 线)更加的平缓，而基于广义极值分布假设的反演结果(黑色曲线)在弹性参 数突变的地方表现的更好。

图9为分别为叠前地震道集数据(a)、使用广义极值分布拟合反射系数 的反演结果的合成地震记录(b)、使用高斯分布拟合发射系数的反演结果的 合成地震记录(c)。可以观察到基于广义极值分布的方法获得参数反演结果 具有更高的分辨率和精度。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理 解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和 实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种 不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明 的保护范围内。