基于精确Zoeppritz方程的贝叶斯叠前地震随机反演

牛丽萍,胡华锋,周单,郑晓东,耿建华

(1.中石化石油物探技术研究院有限公司,江苏 南京 211103;2. 中国石油勘探开发研究院,北京 100083;3. 同济大学海洋地质国家重点实验室,上海 200092;4.同济大学 海洋与地球科学学院,上海 200092;5.同济大学 海洋资源研究中心,上海 200092)

0 引言

油气地震勘探的目的是探明地下介质的几何形态、岩性和物性甚至流体分布,进而为钻井井位部署提供依据。目前,石油工业界广泛采用叠前AVO反演获得油气储层弹性参数,在岩石物理分析的基础上,对油气储层进行定量刻画。叠前AVO反演的理论基础是描述平面波在地下介质分界面处发生反射和透射的Zoeppritz方程[1]。20世纪,许多学者提出了大量的AVO线性近似公式[2-8],以此来解决Zoeppritz方程求解的非线性问题,并提高反演的计算效率。然而,这些线性AVO近似基于弱反差介质与小入射角的假设,在强反差介质及大角度入射情形下会产生较大的计算误差[9-11]。随着油气勘探的深入,勘探目标转向非均质性强的复杂油气藏,而长排列地震采集以及角度域成像技术的应用也使得有效地震数据的入射角越来越大,因此迫切需要发展适应强反差介质以及大入射角情形下的AVO反演技术。

近些年,很多学者开始研究基于精确Zoeppritz方程的拟线性及非线性反演方法。刘福平等[12]首次给出了纵波反射系数对纵横波速度一阶偏导数的精确表达式;Zhi等[13]利用精确Zoeppritz方程和PP-PS联合反演得到了界面两侧的密度比和3个速度比;Zong等[14]推导了Zoeppritz方程所描述的纵波反射系数的精确表达式,并结合非线性反演方法估计出了反射界面上下层介质的6个弹性参数;Zhi等[15]提出了基于Zoeppritz方程的迭代正则化Levenberg-Marquardt反演方法;Yin等[16]利用纵波速度反射率、横波速度反射率和密度反射率重新表示Zoeppritz方程,然后基于逆算子估计实现三参数的非线性反演;张凌远等[17]利用快速模拟退火算法实现了基于Zoeppritz方程的横波速度、纵横波速度比反演。尽管基于精确Zoeppritz方程的叠前反演方法获得了较大的发展,但是这类反演方法的计算复杂度高,在实际应用中,特别是针对低信噪比数据,反演结果的准确性以及稳定性等依然存在较大挑战[18-20]。

MCMC模拟是一种启发式的全局寻优算法,可以直接模拟任意复杂形式的概率分布。利用MCMC算法进行地震反演,不仅可以避免假设储层弹性参数的分布形式服从高斯分布,而且也可以避免对弹性参数与地震数据之间的物理关系做线性假设。此外,基于MCMC模拟所获得的后验概率分布能够完整地刻画参数估计的不确定性,从而为更科学地利用反演结果对储层进行刻画提供了依据。因此,MCMC算法是一种在复杂储层与大入射角地震数据条件下实现AVO弹性参数反演的理想选择[21-22]。

尽管MCMC算法具有诸多优点,但是在将其应用于叠前弹性参数反演时仍然存在以下两方面的问题:①对于实际复杂地下介质,弹性参数的真实概率分布往往十分复杂,很难寻找一个合适的概率密度函数来描述,常规利用高斯分布刻画弹性参数的统计特征存在较大的局限性;②由于地下模型参数空间巨大以及地震数据中噪声等因素的影响,MCMC对弹性参数后验概率分布的搜索过程极易受到局部极值的影响,这使得MCMC采样过程难以稳定收敛,并且增加了反演结果的不确定性,特别是在低信噪比资料情况下,这种问题更加明显。为了避免这些问题,大多数学者利用Zoeppritz方程的线性近似公式来描述地震正演过程[23-27],并未将MCMC算法应用于求解基于精确Zoeppritz方程的叠前地震反演问题。Pan 等[28]提出了基于精确Zoeppritz方程和MCMC-DRAM(delayed rejection adaptive Metropolis)算法的叠前反演方法,但是DRAM算法无法有效地遍历具有很多局部极值的复杂高维参数空间,而且也难以判断马尔科夫链何时达到收敛[29]。为了降低后验概率分布的复杂度以提高MCMC采样的稳定性与收敛性,Pan等假设弹性参数相互独立且服从高斯分布,然而,该假设条件在实际复杂地下储层中是难以满足的。

本文针对实际复杂储层及低信噪比资料条件下基于精确Zoeppritz方程的叠前地震反演问题,从以下3个方面来提高MCMC叠前AVO反演的稳定性与收敛性:①利用低频模型约束,将待反演参数转换为模型参数的扰动量,从而使得先验概率分布由实际复杂概率分布转换为高斯分布,以降低后验概率分布的复杂度;②通过对似然函数取对数,并利用低频模型来约束地震正演过程;③利用基于随机子空间采样的多链模拟算法对叠前非线性反演问题进行全局寻优,以避免采样过程过早地收敛到局部极值。低信噪比的模拟数据和实际数据测试表明,本文所提方法能够有效提高反演结果的准确性以及稳定性。

1 方法原理

1.1 基于精确Zoeppritz方程的贝叶斯叠前地震反演

贝叶斯理论通过将模型参数视为随机变量,把模型参数的先验分布与似然函数结合起来,从而获得同时反映模型合理性与数据拟合程度的后验概率分布。当给定地震观测数据d时,模型参数x的后验概率分布可以写为[30]:

P(x|d)∝P(x)P(d|x),

(1)

式中:先验分布P(x)表示独立于实际地震观测数据d的模型参数的已知信息;似然函数P(d|x)描述了实际地震观测数据与由模型参数正演所得模拟地震数据之间的拟合程度。

对于实际地下介质而言,弹性参数的先验概率分布形式通常十分复杂,特别是对于复杂岩性储层[31]。为了充分发挥先验低频模型在叠前随机反演过程中的约束作用,我们将待反演参数由弹性参数本身转换为弹性参数的扰动量。由于弹性参数的扰动量围绕弹性参数的低频趋势m0而变化,因此,不管真实的弹性参数m服从简单单峰分布还是复杂多峰分布,弹性参数扰动量Δm一般表现为单峰分布,其先验概率分布可以用高斯分布来表示:

P(Δm)=Ν(Δm;μΔm,ΣΔm),

(2)

式中:μΔm和ΣΔm分别为扰动量的先验期望和先验协方差。根据贝叶斯理论,弹性参数扰动量的后验概率分布可以表示为:

P(Δm|d)∝P(Δm)P(d|Δm),

(3)

假设地震数据中的噪声服从高斯分布,同时考虑到数值计算的稳定性,我们采用对数似然函数:

(4)

dsyn=Sc(m)+e,

(5)

式中:S为地震子波矩阵;e为服从零均值高斯分布的噪声;c(m)为利用精确Zoeppritz方程(见附录A)计算得到的纵波反射系数,且m=m0+Δm。

1.2 基于随机子空间采样的多链模拟

假设待反演地震道的时间采样点数为T,则待反演的弹性参数扰动量个数d=3×T,d即为模型参数空间的维数。为了避免采样过程过早地收敛到局部极值,本文通过同时运行多条初始状态不同的马尔科夫链来搜索模型空间。此外,在对高维模型参数空间的随机搜索过程中,如果对所有模型参数同时进行更新,通常情况下生成的大量候选点会被拒绝,导致马尔科夫链收敛非常缓慢。为了提高候选点的接受率以加速马尔科夫链的收敛过程,在每次迭代时,随机选择某些模型参数组成子空间,只对子空间内的模型参数进行更新,而对子空间之外的模型参数保留原采样点。

1)随机生成弹性参数扰动量的初始样本xi,结合该道的低频模型得到相应的3个弹性参数(纵波速度、横波速度和密度),然后基于精确Zoeppritz方程和褶积模型正演得到该道的模拟地震记录,根据式4计算该样本所对应的似然函数,并将其与弹性参数扰动量的先验分布(式2)相结合得到该样本xi所对应的后验概率密度p(xi);

首先，农村的“空心化”引发的人才短缺，在城镇化、工业化不断发展的当下，大量农村青壮年劳力不断涌入城市，农村常住人口多是一些“老弱妇幼”人员，并且随着农村人口逐渐空心化的同时，逐渐演变发展为农村土地、产业以及各种基础设施的整体“空心化”，当下农村严重缺乏能够满足现代服务业要求的高素质劳动力人口、以及相应基础设施。

2)根据原采样点xi计算第i条链的候选点zi:

zi=xi+dxi。

(6)

式中,dxi为模型参数的更新量,可以表示为:

(7)

式中:B为由随机选择的d′个模型参数组成的子空间;δ为马尔科夫链的组数;μ和ν均为从{1,…,i-1,i+1,…,N}中随机取出δ个整数所组成的向量,且数值无重复;C为尺度因子;ε为从高斯分布Nd′(0,c)中生成的随机数,通常可令c=10-6;

3)计算候选点所对应的后验概率密度p(zi),然后根据Metropolis准则计算接受概率α(xi,zi)[32]:

(8)

4)判断是否接受转移。如果候选点被接受,则xi=zi;否则,xi=xi;

5)对每一个时刻k重复执行步骤2)～4)直至完成N条马尔科夫链的采样,然后进入时刻k+1进行采样直至达到最大迭代次数K。

图1 对弹性参数扰动量的随机采样流程Fig.1 Flow chart of stochastic sampling for the elastic parameter perturbations

2 模拟数据测试

图2 弹性参数的测井曲线及其低频模型Fig.2 Real well logs of elastic parameters and the corresponding LFMs

图3 弹性参数的交会及利用高斯混合分布拟合得到的概率密度等值线Fig.3 Crossplots of the elastic parameters and the corresponding probability density contours obtained by fitting with Gaussian mixture distributions

图4 弹性参数扰动量的交会及利用高斯分布拟合得到的概率密度等值线Fig.4 Crossplots of elastic parameter perturbations and the corresponding probability density contours obtained by fitting with Gaussian distributions

利用精确Zoeppritz方程计算该测井数据的纵波反射系数,并将其与主频为45 Hz的Ricker子波相褶积得到无噪声的地震角道集(图5a),然后加入高斯随机噪声得到最终的模拟地震角道集(图5b),模拟地震数据的信噪比为3,入射角范围为3°～48°,角度采样间隔为3°。图6对比了本文所提改进后的反演方法与改进前的基于Zoeppritz方程反演的效果,可以看出,本文所提方法反演得到的弹性参数与实际测井曲线的拟合度更好。图7为马尔科夫链的收敛性指标随迭代次数的变化,在迭代约32 000次后,马尔科夫链已经达到收敛。图8和图9分别为1 630 ms和1 672 ms处3个弹性参数扰动量的后验统计直方图,不同位置处3个弹性参数扰动量均稳定收敛到后验高斯分布。

a—无噪声;b—信噪比为3a—no noise;b—S/N is 3

图6 含噪声地震数据的弹性参数反演结果Fig.6 The inversion results of the elastic parameters with noisy seismic data

图7 收敛性指标随迭代次数的变化Fig.7 The variation of convergence diagnostic with the number of iterations

图8 时间为1 630 ms处弹性参数扰动量的后验统计直方图Fig.8 Statistical histogram of the sampling results for the elastic parameter perturbations at time 1630 ms

图9 时间为1 672 ms处弹性参数扰动量的后验统计直方图Fig.9 Statistical histogram of the sampling results for the elastic parameter perturbations at time 1672 ms

3 实际数据应用

实际二维地震数据的叠加剖面及测井位置如图10所示,该数据的时间范围为2 100～2 350 ms,采样率是2 ms,地震子波主频为20 Hz。图11展示了不同CDP位置处的地震角道集,角度范围为1°～31°,角度间隔为2°。从图10和图11可以看出,该剖面左侧地震数据的信噪比较低,而在剖面的中间及右侧地震数据的信噪比较高。

图10 叠加地震剖面Fig.10 Poststack seismic section

图11 不同CDP位置处的实际地震角道集Fig.11 Actual seismic angle gathers at different CDP

在地震层位的约束下,我们首先根据测井数据构建弹性参数的低频模型(图12),然后通过对测井位置处弹性参数扰动量的统计分析得到其先验期望和先验协方差,图13展示了弹性参数扰动量的交会及其高斯分布拟合结果,可以看出,高斯分布能够很好地拟合弹性参数扰动量的统计分布特征。

a—纵波速度;b—横波速度;c—密度a—P-velocity;b—S-velocity;c—density

图14为利用本文所提方法反演得到的弹性参数剖面,图15和图16分别展示了W1井和W2井处的反演结果。从两口井的测井曲线可以看出,该工区的弹性参数在2 200 ms附近存在一个非常明显的突变界面,在该界面以下,两口井的测井曲线存在明显的差异,这反映了该工区地层横向变化较大。二维剖面的反演结果与该先验认识一致。从井旁道的反演结果来看,本文所提方法的反演结果能够较好地拟合实际测井曲线,而且基于反演结果的模拟地震角道集也与实际地震角道集相吻合。

a—纵波速度;b—横波速度;c—密度a—P-velocity;b—S-velocity;c—density

a—纵波速度;b—横波速度;c—密度;d—井旁地震角道集a—P-velocity;b—S-velocity;c—density;d—the seismic angle gather at well

a—纵波速度;b—横波速度;c—密度;d—井旁地震角道集a—P-velocity;b—S-velocity;c—density;d—the seismic angle gather at well

a—纵波速度; b—横波速度; c—密度

图18 不同CDP位置处收敛性指标随迭代次数的变化Fig.18 The variation of convergence diagnostic with the number of iterations at different CDP

4 结论

本文针对实际复杂储层及低信噪比资料条件下基于精确Zoeppritz方程的叠前地震反演问题,通过低频模型约束先验概率分布以及地震正演模拟过程,并结合对数似然函数和基于随机子空间采样的多链模拟算法,提出了一种改进的MCMC叠前地震反演方法。低信噪比的模拟数据测试表明,相较于改进前的反演方法,本文所提方法能够获得更加符合实际测井曲线的弹性参数反演结果。二维实际资料的应用效果表明,尽管本文所提方法为逐道反演,但是剖面上相邻两道的反演结果之间具有良好的一致性,这很好地说明了本方法的稳定性;此外,当地震资料的信噪比越低时,弹性参数反演结果的不确定性越大,本文所提方法能够给出可信的、定量的不确定性估计。

附录A Zoeppritz方程

在各向同性的弹性介质中,当平面纵波入射到两层半无限介质的分界面处时,反射系数、透射系数与介质弹性参数和入射角之间的关系可以表示为[1]:

式中:VP,VS和ρ分别表示纵波速度、横波速度和密度,其下标1和2分别表示上层和下层介质;θ1和φ1分别为P波和SV波的反射角;θ2和φ2分别为P波和SV波的透射角;RPP和RPS分别为P波和SV波的反射系数;TPP和TPS分别为P波和SV波的透射系数。

附录B 收敛性检验

对于N条采样长度为T的马尔科夫链,每条链内的方差可以表示为:

其中,⎣·」表示取整运算;而所有链之间的方差可以表示为:

那么,马尔科夫链的收敛性检验指标可以表示为: