与叠后地震反演相比,叠前反演能获得多种储层弹性参数反演结果,可以针对不同的地质状况选择不同的参数反演结果来实现储层识别(Sunetal.,2014).随着油气勘探逐步深入,具有复杂地质构造的油气藏逐渐成为开发的重点,常规的线性反演方法受地震数据频带跟子波影响,很难有效地进行储层预测.而以地质统计学反演为代表的随机反演方法充分利用了测井数据在分辨率上的优势,可以有效地识别地下薄储层;王保丽等,2015;孙瑞莹等,2016;沈洪涛等,2017).

叠前地质统计学反演属于多参数反演方法,因此需要考虑反演参数之间的相关性.常规的地质统计学反演方法主要运用了序贯高斯模拟方法来实现数据的模拟(姜晓宇等,2014).序贯高斯模拟属于条件模拟方法,其首先对硬数据进行高斯变换,之后依次获取各个节点的条件累计分布函数并对其进行抽样,最后对样本值进行高斯反变换,得到模拟结果.由于对模拟数据进行高斯变换的影响,当进行多参数模拟时,序贯高斯模拟就难以再现模拟参数之间的统计分布特征(Soares,2001).针对该问题,Soares(2001)提出了直接序贯模拟(DirectSequentialSimulation,DSS)和直接序贯协模拟(DirectSequentialCo-Simulation,Co-DSS)方法.该模拟方法可以直接在原始数据域对数据进行模拟,不需要进行高斯变换,其利用克里金估计得到的局部均值和方差直接对全局条件分布函数进行采样,因此模拟参数之间的相关特性可以得到很好地体现.之后,基于联合分布的直接序贯协模拟方法(HortaandSoares,2010;Nunesetal.,2016)的提出更进一步解决了再现非线性相关变量之间统计特点的问题.

传统地质统计学反演的基础是逐道反演方法(HaasandDubrule,1994),其受地震信噪比影响,反演结果往往存在着误差累计.而基于全局随机反演策略(Azevedoetal.,2013,2015)的地质统计学反演方法将模拟及协模拟方法作为模型扰动技术,将模拟结果正演得到的合成地震记录与实际地震记录之间的相关系数作为模拟结果的优选标准以及模型扰动的约束条件,经过不断迭代,可以得到与实际地震记录最为匹配的模拟结果.然而基于直接序贯模拟及协模拟的全局迭代地质统计学反演方法仅通过比较不同模拟结果所合成的地震记录与观测地震记录之间的局部相关系数大小来实现对模拟结果的优选,存在一定的不稳定性(Pereiraetal.,2016).针对该问题,本文引入Metropolis优化算法(Metropolisetal.,1953;Hasting,1970)以及考虑先验模型约束及横向约束的目标函数来对模拟结果进行优选,其不仅可以有效避免非线性问题中的局部极值问题,还添加了先验模型约束项与横向约束项(HamidandPidlisecky,2015)来提高反演结果的稳定性及横向连续性.同时,我们根据目标函数的表达形式,重新构建了基于目标函数的优化局部相关系数表达形式,提高了反演的稳定性.

本文将基于联合分布的直接序贯协模拟方法与Metropolis优化算法相结合,在贝叶斯框架下构建目标函数,并构建新的优化局部相关系数,提出了基于Metropolis优化算法和优化局部相关系数的全局迭代地质统计学反演方法.首先利用直接序贯模拟及基于联合分布的直接序贯协模拟得到纵横波阻抗及密度的多个模拟结果,并通过正演关系得到合成地震记录;之后计算对应的目标函数值,利用Metropolis准则对模拟结果进行优选;最后根据目标函数计算当前状态下的优化局部相关系数,并将其与当前反演结果作为条件数据,利用直接序贯协模拟及基于联合分布的直接序贯协模拟重新模拟得到新的模拟结果,并再次利用Metropolis准则对模拟结果进行优选;经过一系列迭代,直到得到最终的反演结果.

1、理论与方法

1.1 直接序贯模拟

直接序贯模拟算法可以利用由克里金估计得到的局部均值和方差,直接对模拟参数的全局条件概率分布函数进行抽样.其基础是简单克里金估计(Goovaerts,1997),克里金均值Z(x0)*及克里金方差σ2E(x0)可写作:

(xi)为模拟数据在临近点xi的均值,λi为临近点xi的权重系数,Z(xi)表示模拟数据在临近点xi的变量值,C(0)为模拟数据在待估计点x0的方差,C(xi,x0)为模拟数据在待估计点x0与临近点xi的协方差,N表示参与简单克里金计算的临近点个数.其中,权重系数可通过式(3)计算:

根据硬数据的分布特征,定义一个全局累计高斯分布函数G(y(x0)*,σ2E(x0)),其中y(x0)*为克里金均值Z(x0)*所对应的高斯值.在高斯分布函数G(y(x0)*,σ2E(x0))中产生一个随机的高斯数yG,其所对应的反演参数值即为模拟值Zs(x0)(AzevedoandSoares,2017).

相比传统的序贯高斯模拟,直接序贯模拟不需要原始数据进行任何非线性变换,通过对变量之间的联合概率分布进行抽样,可以有效地再现变量之间的统计关系.

1.2 基于联合分布的直接序贯协模拟

对于叠前地质统计学反演来说,需要模拟得到多个变量(如纵波阻抗、密度、横波阻抗)的模拟结果,这需要考虑变量之间的相关性.传统的协模拟方法只适用于变量之间满足线性相关的情况,当变量之间并不满足完全线性关系时,协模拟方法就难以表征变量之间的统计分布特征(肖亚楠,2018).

针对该问题,Horta和Soares(2010)提出了基于联合分布的直接序贯协模拟方法.在已知其中一个变量的模拟结果的情况下,根据变量之间的联合概率密度分布可以计算其他变量的条件概率分布,通过对条件概率分布进行抽样来再现变量之间的统计关系.

其基础是同位简单协克里金(Xuetal.,1992;AlmeidaandJournel,1994),假设Z2为待模拟变量,Z1为已知变量,则待模拟变量Z2在待估计点x0处的克里金均值Z2(x0)*与克里金方差σ2E(x0)可用式(6)和(7)进行表征:

式中,CZ2(xi,xj)是待模拟变量Z2在临近点xi与xj之间的协方差,CZ1(0)为已知变量Z1在待估计点x0处的方差,CZ1Z2(x0,xi)为待模拟变量Z2与已知变量Z1在临近点xi处与待估计点x0处之间的协方差,CZ1Z2(xj,xi)为待模拟变量Z2与已知变量Z1在临近点xi与xj处之间的协方差,CZ1Z2(0)为已知变量Z1与待模拟变量Z2在待估计点x0处之间的协方差.由于相关系数实际上是未考虑量纲的协方差,待模拟变量Z2与已知变量Z1的协方差CZ1Z2(xi,xj)、CZ1Z2(x0,xi)及CZ1Z2(0)可以利用两者的局部相关系数进行转化,因此局部相关系数可以通过计算克里金权重,应用到随机模拟之中.

基于联合分布的直接序贯模拟方法通过对变量之间的全局联合概率分布函数抽样,不仅可以保证模拟结果满足变量之间的联合概率分布特性,还可以提高模拟的精度.

利用Marmous2模型的一部分进行模拟测试.图1为密度模型(a)、密度伪井数据(b)及已知的纵波阻抗模拟结果(c).图2为纵波阻抗与密度的联合概率分布及某模拟点的条件概率分布.图3为直接序贯模拟、直接序贯协模拟及基于联合概率分布的直接序贯协模拟的一次模拟结果.可以发现,基于联合分布的直接序贯协模拟方法的模拟结果更加精确.

1.3 基于优化局部相关系数的叠前地质统计学反演方法

本文提出的叠前地质统计学反演方法是在贝叶斯理论的框架下,将直接序贯模拟及基于联合分布的直接序贯协模拟作为模型扰动技术,将考虑先验模型约束项及横向约束项的目标函数及Metropolis准则作为模拟结果的优选手段,并将基于目标函数的优化局部相关系数作为条件约束的全局迭代随机反演方法.

(1)Zoeppritz近似方程及正演关系

叠前地震反演的理论基础是Zoeppritz方程,然而精确的Zoeppritz方程极其复杂,不利于反演(李建华等,2016).因此,许多学者(AkiandRichards,1980;Shuey,1985;Gidlowetal.,1992)对其进行了简化,提出了Zoeppritz近似方程.

本文主要根据Fatti近似方程(Fattietal.,1994)来进行叠前反演.该近似公式具体形式为

图1密度模型及已知的纵波阻抗模拟结果

图2纵波阻抗与密度的联合分布及密度条件概率分布直方图

图3不同模拟结果的对比

(3)基于贝叶斯定理的目标函数构建方法及优化局部相关系数的构建

利用Metropolis优化算法对模型进行优选,需要建立似然函数.本文在贝叶斯理论的框架下,构建了反演的似然函数,并引入了由测井数据构建的平滑模型约束信息及横向约束项,来提高反演的稳定性和横向连续性.

1.4 基于全局迭代反演策略的叠前随机反演流程

本文根据全局随机反演策略,提出了基于Metropolis优化算法的全局迭代反演方法.考虑到纵波阻抗反射系数在AVO近似方程中的权重问题,将纵波阻抗作为横波阻抗及密度模拟的条件数据.因此,其具体步骤如下:

(1)将已知的测井数据为硬数据,运用直接序贯模拟方法模拟得到纵波阻抗的NS个实现,并利用纵波阻抗模拟结果计算每个模拟点横波阻抗及密度的条件概率分布函数.

(2)利用同位简单协克里金方法计算每个模拟点横波阻抗与密度的局部均值与方差,并利用局部均值和方差对条件概率分布函数进行抽样,得到横波阻抗及密度的NS个实现.

(3)利用模拟得到的纵波阻抗、横波阻抗和密度,运用Fatti近似方程,正演得到叠前合成地震记录.

(4)利用贝叶斯定理建立目标函数,分区块计算每一区域的目标函数值,利用Metropolis准则判断是否进行状态更新,并利用当前状态下的纵横波阻抗及密度值计算每一区域的优化局部相关系数.

(5)将当前纵横波阻抗、密度模拟结果及优化局部相关系数作为软数据,利用直接序贯协模拟及考虑联合分布的直接序贯协模拟得到新的模拟结果,重复步骤3～5,直到达到最大迭代次数.

2、模型试算

为了验证反演的可行性,利用Marmous2模型的一部分,进行了叠前反演试算.设置最高迭代次数为10000次,每次迭代利用协模拟方法分别模拟得到纵横波阻抗及密度的10次模拟结果,利用Metropolis准则选取其中的最优结果,并将其作为下次迭代的条件数据.为提高算法的计算效率,本文利用了并行策略.图4为纵波阻抗、横波阻抗和密度模型及其纵向横向变差函数.为了分别测试Metropolis优化算法及优化局部相关系数对反演结果的影响,分别利用基于Metropolis优化算法和优化局部相关系数的反演方法、基于Metropolis优化算法和原局部相关系数的反演方法以及只考虑原局部相关系数大小变化的反演方法,前2种反演方法目标函数一致(均为式(14)),仅局部相关系数的计算公式不同(分别为式(16)、式(15)),反演结果分别如图5所示,图6为反演结果与真实模型的相对误差绝对值.图7为模型第18道真实模型与不同方法反演结果的对比,图8为所对应的合成角道集与真实角道集的对比,可以发现在相同迭代次数的情况下,三种反演方法的反演结果均能获得较好的反演结果,但基于优化相关系数与Metropolis优化算法的反演结果相较于其他反演结果稳定性更高,总体误差更小.

图4纵波阻抗、横波阻抗和密度模型数据及纵向横向变差函数

(a)(b)(c）纵波阻抗模型及其对应的纵向横向变差函数；（d)(e)(f）横波阻抗模型及其对应的纵向横向变差函数；（g)(h)(i）密度模型及其对应的纵向横向变差函数．

图9为基于Metropolis优化算法分别利用原相关系数与优化局部相关系数两种反演方法的目标函数值随迭代次数的变化,可以发现利用优化局部相关系数的反演方法目标函数能更加稳定地收敛,且最终结果的目标函数值更小,在并行计算的情况下,共耗时12min17s.而基于原局部相关系数的反演方法在8000次迭代左右,才基本达到平稳状态,且收敛过程具有较大的起伏.图10为某采样点达到平稳状态下纵横波阻抗及密度的后验分布直方图及此时纵波阻抗模拟结果所对应的横波阻抗与密度的条件分布直方图,可以发现横波阻抗和密度的后验概率分布均位于其条件概率分布之中且满足正态分布,因此可以选取达到平稳状态下多个反演结果的均值作为反演结果.

从模型测试中可以看出,基于Metropolis优化算法及优化局部相关系数的叠前地质统计学反演是可以适用于叠前反演的,且稳定性与收敛性相比传统的考虑相关系数变化的迭代反演方法有了一定程度的提高.

图5不同反演方法的反演结果对比

(a)(d)(g）基于Metropolis优化算法及优化局部相关系数的反演结果；（b)(e)(h）基于Metropolis优化算法及原局部相关系数的反演结果；（c)(f)(i）仅考虑原局部相关系数的反演结果．

当构造起伏较大时,等时求取横向变差函数及定义横向约束项具有较大误差,此时可利用层位数据求取沿层横向变差函数,并定义对应的横向约束项,消除地层起伏因素的影响.因此,我们根据模型数据,选取底部1.45～1.5s之间的部分数据,对其进行一定程度的旋转,以此为模型进行大倾角测试.图11为模型数据及反演的结果,可以发现利用沿层横向变差函数及对应横向约束项的反演结果在大倾角状况下具有更好的反演效果,因此在大倾角状况下沿层位求取变差函数及定义横向约束项是较为必要的.

3、实际数据测试

为了验证基于Metropolis抽样及优化局部相关系数的全局迭代叠前地质统计学反演在实际数据中的应用效果,利用实际数据进行反演测试.该数据来自东部某工区,属于致密砂岩储层,储层内部“甜点”表征低孔低渗背景下局部高孔渗的区域,是勘探开发中最为关注的地质目标.经过岩石物理分析得知,纵波阻抗可以很好地表征“甜点”的分布.利用井数据求取纵向变差函数,利用确定性反演结果求取横向变差函数.为了提高反演的稳定性,将叠前角度道集叠加为三个部分角度叠加数据,其中心角度分别为7°、15°、24°.利用部分叠加地震数据进行AVA反演测试.

图6不同反演方法反演结果的相对误差绝对值

(a)(d)(g）基于Metropolis优化算法及优化局部相关系数的反演方法；（b)(e)(h）基于Metropolis优化算法及原局部相关系数的反演方法；（c)(f)(i）仅考虑原局部相关系数的反演方法．

纵横波阻抗及密度的反演结果以及确定性反演纵波阻抗的反演剖面如图12所示,可以发现,井曲线(井位置位于CDP=64)与反演剖面是相吻合的.图13为反演结果所合成的中角度地震记录与实际中角度地震记录的对比,可以发现两者较为接近,说明反演结果所合成的地震记录与实际地震记录具有较好的相关性.图14是目标函数随迭代次数的变化,我们可以发现迭代5000次以后,反演基本收敛,因此我们可以选取收敛后一系列反演结果的均值作为最终反演结果,即为图12所示的剖面.图15展示了全局迭代叠前地质统计学反演及确定性反演在井位置处(CDP=64)的反演结果,可以发现地质统计学反演结果在井位置处反演结果与真实井数据完全吻合,这是由于地质统计学反演利用井数据作为硬数据进行反演,因此相比确定性反演具有更高的分辨率,可以很好反映确定性反演没能表征出的2.13s左右的薄层.因此,基于Metropolis优化算法及优化局部相关系数的全局迭代地质统计学反演可以较好地应用于叠前实际数据反演.

图7模型第18道真实模型与三种反演结果的对比

(a）基于优化局部相关系数与Metropolis优化算法的反演方法；（b）基于原局部相关系数与Metropolis优化算法的反演方法；（c）仅考虑原局部相关系数的反演方法．红色虚线：真实模型，蓝色实线：反演结果．

图8模型第18道的真实角道集与三种反演结果所合成的角道集对比

(a）真实角道集；（b）基于Metropolis优化及优化局部相关系数反演结果所合成角道集；（c）基于Metropolis优化及原局部相关系数反演结果所合成角道集；（d）仅考虑原局部相关系数反演结果所合成角道集．

4、结论

本文将Metropolis优化算法与全局迭代地质统计学反演方法相结合,提出了基于Metropolis优化算法的全局迭代地质统计学反演方法,该方法的优势主要存在以下几点:

(1)在贝叶斯理论的框架下,建立目标函数,并引入横向约束,通过Metropolis优化算法对模拟结果进行优选,可以克服原有全局迭代地质统计学方法仅利用地震记录局部相关系数进行优化所带来的不确定性,提高反演的稳定性.

(2)根据目标函数构建了新的优化局部相关系数,将新的基于目标函数的局部相关系数应用于协模拟之中,提高了叠前地质统计学反演的收敛性与稳定性.

(3)模型测试与实际数据应用表明,基于Metropolis抽样算法与优化局部相关系数的全局迭代地质统计学反演方法可以有效地应用于叠前反演,且能较好地反演较薄储层.

图9基于Metropolis优化算法的地质统计学反演目标函数值随迭代次数变化

图10某采样点纵横波阻抗及密度的后验分布直方图及横波阻抗与密度的条件分布直方图

图11大倾角测试模型及反演结果

图12纵波阻抗、横波阻抗、密度反演剖面

图13反演结果所合成地震记录与实际地震记录的对比

图14目标函数随迭代次数变化

图15井位置处（CDP=64）反演结果与井数据对比

参考文献：

[1]姜晓宇,计智锋,毛凤军等.2014.随机模拟技术在储层物性研究中的应用.地球物理学进展,29(6):2665-2668.

[2]李建华,刘百红,张延庆等.2016.叠前AVO反演在储层含油气性预测中的应用.石油地球物理勘探,51(6):1180-1186.

[3]刘婵娟.2017.基于横向约束的快速随机反演方法研究[硕士论文].青岛:中国石油大学(华东).

[4]沈洪涛,郭乃川,秦童等.2017.地质统计学反演技术在超薄储层预测中的应用.地球物理学进展,32(1):248-253.

[5]孙瑞莹,印兴耀,王保丽等.2016.基于随机地震反演的Russell流体因子直接估算方法.地球物理学报,59(3):1143-1150.

[6]王保丽,印兴耀,丁龙翔等.2015.基于FFT-MA谱模拟的快速随机反演方法研究.地球物理学报,58(2):664-673.

[7]肖亚楠.2018.基于序贯协模拟的地震随机反演方法研究[硕士论文].青岛:中国石油大学(华东).

赵晨,张广智,张佳佳,郗诚,肖亚楠.基于Metropolis优化的叠前全局迭代地质统计学反演方法[J].地球物理学报,2020,63(08):3116-3130.

基金:国家科技重大专项(2016ZX05027004-001);国家自然基金项目(41674130);中央高校基本科研业务费专项资金(18CX06020A)联合资助.