基于自动相关判别先验的叠前同时反演方法研究

纪永祯,朱立华,林正良,胡华锋

(中国石油化工股份有限公司石油物探技术研究院,江苏南京211103)

叠前振幅随角度变化(AVA)反演技术是利用叠前地震数据的振幅随角度变化特征,采用反演的手段获取岩石弹性参数的技术[1-4]。相较于叠后反演,由于利用了不同入射角数据的信息,可以获得更多样的弹性参数估计,因而在储层描述中起着至关重要的作用[5]。但由于利用的叠前地震数据与弹性参数的映射关系更加复杂,信噪比更低,因而叠前反演方法需要应对更多难题[6-8]。其中一个难题就是降低解的非唯一性,通常的解决办法是在目标函数中加入正则化项或在贝叶斯框架下加入待估计参数的先验信息约束,在许多可能的解中选择一个最优或最符合假设的解,且在贝叶斯框架下加入待求解参数的先验信息约束,在一定条件下与正则化约束具有相似的内涵,可以推导成相同形式[9]。

基于贝叶斯理论框架的叠前反演研究将似然函数与先验地质信息结合,利用求解最大后验概率来建立反演的目标函数,进而反演弹性参数。其中先验信息的正则约束作用可以较好地解决叠前反演中所存在的非唯一性和病态问题[9]。在构建反演目标函数时,似然函数的构建一般体现反演参数的正演结果与地震数据的匹配关系,但先验信息的构建并没有统一的假设及论述。

BULAND等[10]通过概率统计分析认为弹性参数(特别是声波速度)大体符合高斯分布,基于此构建出贝叶斯反演方法并被广泛应用。DOWNTON等[11]将长尾分布作为先验分布,获得了稀疏的反射系数估计,并认为常规高斯分布先验(这类先验约束与最小二乘本质相同)获得的结果比较光滑,不利于高分辨率分析和解释,发现柯西分布先验比高斯分布先验分辨率更高。此后,利用“长尾分布”作为先验受到关注,其中较为常用的是柯西分布及其变体(如修正柯西分布等)。THEUNE等[12]分析了弹性参数反射系数稀疏假设下的两个先验模型,发现可微拉普拉斯分布先验可获得一个凸的目标函数,而柯西分布先验获得的是非凸函数,因此可微拉普拉斯先验更具优势。ALEMIE等[13]将多元柯西分布先验引入贝叶斯反演,这种多元先验分布可以合并模型参数之间的相关性,并获得高分辨率的结果。随后,Huber分布、t分布等先验先后被引入到贝叶斯反演中,均获得了较好的效果[9,14-16]。以上先验均假设反射系数是随机变量,其分布满足某一特定分布,然而实际弹性参数反射系数的分布受独特的地质条件影响,常表现出独特的特征,无法用某一特定分布准确描述;不同弹性参数的反射系数往往具有不同的分布,对拟反演的不同弹性参数的反射系数进行同时反演时,采用同一分布进行约束也不够合理;并且,特定分布一般由测井资料统计获取,一经确定则不随地震数据变化,而地质特征的横向变化引起的反射系数分布变化难以得到反映。因此,弹性参数反射系数的反演结果可能出现误差。自动相关性判别(automatic relevance determination,ARD)先验可以解决上述问题,这种先验信息不再假设反射系数是符合特定分布的随机变量,而是给每一个反射系数分配一个参数,通过估计该参数获得对反射系数的最终估计。该参数作为一个待估计的特定值,由迭代学习自适应确定,并随着不同地震道的特征和不同弹性参数反射系数的特征变化。这种先验会保持反射系数的稀疏性[17-18],从而使反演结果具有较高的分辨率,其在文字处理、图像分类和动态光散射分析的应用中显现了较好的潜力[19-21]。

本文通过引入ARD先验,将弹性参数反射系数的先验信息融合到反演中,降低反演的不适定性和非唯一性。为了降低噪声对弹性参数反射系数的振幅保持的影响,在反演过程中,引入低频趋势约束。并且,考虑到时间域叠前反演(利用地震数据的所有频率成分)常受数据中5～100Hz以外低信噪比信息的影响,采用了频域叠前反演策略。另外,考虑到常规基于弹性阻抗的叠前反演方法尽管具备较为稳定的优势[22-23],但反演过程中需要先获取分角度弹性阻抗估计,再转换为弹性参数的估计,会存在一定的转换误差,因此,采用同时反演策略,在反演过程中利用各角度信息直接同时获取不同弹性参数的估计,可以减少弹性阻抗转换为弹性参数时的转换误差。

1 方法原理

具有不同入射角度的地震叠前数据可以看成是地层界面处的反射系数与不同角度子波的褶积,地层界面处的反射系数随入射角度的变化规律可以用Zoeppritz方程及其线性近似式描述[24],其中Fatti近似式如下[25]:

(1)

式中:RPP(θ)是随角度变化的地层界面反射系数,由入射角和弹性参数反射系数决定;RIP,RIS,Rρ分别是纵波反射系数、横波反射系数和密度反射系数;γ是横波速度与纵波速度的比值;θ是入射角。考虑到大入射角地震数据较难获取,导致密度反演不准确,进而会影响到纵横波阻抗的反演稳定性,且Fatti公式可以在省略密度项的同时保持对地层反射系数的正演精度[25],因此本文将采用两项Fatti公式作为正反演中采用的近似式。

将两项Fatti公式写成矩阵形式:

RPP(θ)=[A(θ)B(θ)][RIPRIS]T

(2)

式中:A(θ)=(1+tan2θ)/2,B(θ)=-4γ2sin2θ。在等式两边同时进行傅里叶变换,可以获得单一入射角叠前数据在频率域的正演矩阵形式:

(3)

式中:

如果采用N个入射角的叠前地震数据,即θ∈[θ1,θN],那么正演矩阵形式如下:

(4)

式中:

其中,diag[·]代表对角矩阵。

为了提高反演过程的稳定性,在以上正演矩阵中加入趋势约束,形式如下:

(5)

式中:Ψ=F-1ΛFC∈RK×K,Λ=diag[H],Λ为以汉宁窗函数向量(H)作为对角元素的对角矩阵,C是积分矩阵,F和F-1分别为正、反傅里叶变换矩阵。趋势约束矩阵Ψ的物理意义是,首先利用积分矩阵C将弹性参数反射系数转换成对应的阻抗对数,然后利用正反傅里叶变换矩阵F和F-1以及汉宁窗矩阵Λ将低频趋势提取出来,与模型趋势进行匹配约束[3]。LowIP=F-1ΛF[log(IP)]∈RK×1,LowIS=F-1ΛF[log(IS)]∈RK×1,分别是从模型中获取的纵波阻抗和横波阻抗的趋势,IP和IS分别为纵波阻抗和横波阻抗;λ1,λ2分别是趋势约束的权重系数。

为便于公式推导和理解,将公式(5)写成如下简单形式:

dMN+2K=G(MN+2K)×2Km2K+n

(6)

式中:d代表公式(5)中等号左边的项;G包含傅里叶变换矩阵、角度相关矩阵和低频趋势提取矩阵;m代表待估计的纵横波阻抗参数反射系数向量;n代表噪声向量。

在贝叶斯框架下,假设噪声向量符合均值为0,方差为σ2的高斯分布,那么地震数据的似然函数有如下形式[26]:

p(d|m,σ2,θ)=(2πσ2)-MNexp·

(7)

引入ARD先验信息,形式如下[18,27]:

(8)

式中:h=[h1,h2,…,h2K]T,其包含2K个独立参数,这些参数决定了每一个时间采样点处的反射系数的幅值,物理意义是每一个反射系数mk均服从零均值高斯分布,但每一个反射系数的方差均不相同,由独立参数hk表示。TIPPING[27]和WIPF等[18]证实了该先验可以有效保持待估计弹性参数反射系数的稀疏特征。根据贝叶斯理论,待估计参数,即弹性参数反射系数的后验概率有如下形式:

(9)

Σ=(H+σ-2GTG)-1

(10)

μ=σ-2ΣGTd

(11)

其中,C为常数,H=diag(h1,h2,…,h2K)是一个对角矩阵,μ是待估计弹性参数反射系数的后验概率均值,即由弹性参数反射系数组成的反演结果向量。由公式(10)和公式(11)可见,为了获取μ的估计,需要估计H和σ2。本文利用输入的叠前地震数据估计方差σ2;采用第二型最大似然估计方法,将(12)式的边缘似然函数最大化获得H的估计。

p(d|h,σ2,θ)=(2π)MN|Q|exp(dTQ-1d)

(12)

式中:Q=σ2I+GH-1GT。最大化边缘似然函数的求解方法可以分为3类,具体可以参考TIPPING[27],MACKAY[28]和YUAN等[29]文献。最后,将获得的H和σ2的估计代入公式(11)获得μ,即弹性参数反射系数[RIP,RIS]T的反演结果。

获得了纵波阻抗和横波阻抗反射系数的估计后,采用公式(13),可以获得纵波阻抗和横波阻抗的估计[3]。

(13)

式中:IP0和IS0分别是纵波阻抗和横波阻抗的初值。

从方法的原理介绍中可以得出:引入的ARD先验,将反射系数先验信息融合到反演中,不再假设反射系数是符合特定分布的随机变量,而是给每一个反射系数分配一个参数,通过估计该参数获得反射系数的最终估计。该方法对反射系数的估计不依赖初始模型,完全由地震数据驱动获得。在反演过程中加入的趋势约束虽然在一定程度上依赖于初始模型,但也可以从偏移速度等数据中获取,使得反演结果不会过多地受到井数据影响。采用在频率域进行反演的策略,可以通过选取优势频带避开叠前数据中的低信噪比成分,提高反演的准确度。由叠前地震数据直接获得弹性参数的估计,避免了传统基于弹性阻抗的反演方法中,从叠前数据到弹性阻抗再到弹性参数的转换误差。

2 模拟数据测试

设计的弹性参数模型和采用的趋势约束(参数模型的0～6Hz成分)如图1a所示,基于该模型,利用公式(1)获得了不同入射角的反射系数,并与一个主频为30Hz,采样间隔为1ms的雷克子波进行褶积得到模拟的不同入射角的叠前数据,如图1b所示。加入20%的随机噪声得到含噪声的模拟数据,如图1c所示。图1d和图1e分别给出了未加趋势约束的无噪声和含噪声数据反演结果;图1f和图1g分别给出了加趋势约束的无噪声和含噪声数据反演结果。反演时,分别选取5～120Hz和10～70Hz作为无噪声和含噪声数据的频带输入。由图1d和图1f可见,在无噪声情况下,是否加入趋势约束并未对反演结果产生过多影响,反演结果均非常好地重现了模型特征。图1g的反演结果优于图1e,体现了趋势约束具有一定的抗噪性。加趋势约束的反演结果(图1f、图1g)较好地重现了模型特征,展示了趋势约束的有效性。由模拟数据的反演实验可见,本文提出的方法具有可行性和抗噪性;ARD先验使得反演结果(无论加入趋势约束与否)分界面清晰,“块化”效果明显,具有保持边界的优势,可提供分辨率更高,更有利于后续储层和地质解释的反演结果;趋势约束的加入可以有效减小由于随机噪声引起的弹性参数反射系数的幅值畸变,有利于反演结果保持稳定并且使反演结果的准确性更高。

图1 无噪和含噪模拟数据对本文方法的测试结果a 弹性参数模型和趋势约束; b 无噪声模拟数据; c 含噪声模拟数据; d 未加趋势约束的无噪声数据反演结果; e 未加趋势约束的含噪声数据反演结果; f 加趋势约束的无噪声数据反演结果; g 加趋势约束的含噪声数据反演结果

采用经典的推覆体2D模型进一步测试本文方法,利用30Hz雷克子波并加入10%的随机噪声获得了含噪声的、采样间隔为1ms的不同角度模拟数据(入射角度分别为10°、20°和30°),模拟数据的纵横波阻抗参数模型如图2a和图2b所示;利用ARD先验,但未加入趋势约束(参数模型的0～6Hz成分)的反演结果如图2c和图2d所示;加入趋势约束的反演结果如图2e和图2f所示。总体上看,反演获得的纵波阻抗结果优于横波阻抗结果;趋势约束的加入提高了反演结果的横向连续性,并且减轻了随机噪声引起的“挂面条”现象,这一提高对于稳定性相对较差的横波阻抗反演结果更为明显。2D模型的测试结果进一步证实了本文方法在获得高分辨率反演结果和保持反演结果横向稳定性方面的有效性和优势。

图2 2D推覆体模型模拟数据对本文方法的测试结果a 纵波阻抗参数模型; b 横波阻抗参数模型; c 未加趋势约束的纵波阻抗反演结果; d 未加趋势约束的横波阻抗反演结果; e 加趋势约束的纵波阻抗反演结果; f 加趋势约束的横波阻抗反演结果

3 实际数据应用

实际数据来自川西地区,该地区目的层发育有多期河道储层。首先,选取了一个井旁道地震数据作为方法的实验数据。不同入射角的叠前井旁道地震数据如图3a所示,采用的入射角度分别为5°、15°和25°。反演时选取10～70Hz作为输入数据的频带范围。对应的趋势约束曲线(弹性参数测井曲线的低频成分)、弹性参数测井曲线与反演结果如图3b和图3c 所示。由图3b和图3c可见,反演结果展示了与测井曲线较好的匹配关系,并且反演结果具有较好的边界特征,较好地反映了地层弹性参数的变化特征,其中蓝色虚线框中较薄的储层也被较好地展现出来,有利于储层识别和后续地质解释。尽管在180～220ms薄层的连续变化特征没有被完美地重现,反演结果仍然对后续的储层和地质解释具有重要意义。

图3 实际井旁地震道数据反演结果a 不同入射角井旁地震道数据; b 叠合了测井曲线和趋势约束的纵波阻抗反演结果; c 叠合了测井曲线和趋势约束的横波阻抗反演结果

选取实际数据中某二维剖面作为实验数据。不同入射角的叠前部分叠加地震数据如图4所示,采用的入射角分别为7°、14°和21°。采用的低频趋势约束由提取测井弹性参数曲线插值模型的低频成分(0～6Hz)获得。通过井震标定获得不同角度子波进行反演。反演时选取10～70Hz作为输入数据的频带范围,反演后将结果中高于150Hz的成分滤除。图5a和图5b分别给出了采用本文方法反演得到的纵波阻抗和横波阻抗。由图5a和图5b可见,本文方法得到的反演结果分辨率较高,地层界面清晰,与地质认识较为吻合,为后续储层和地质解释提供了较好的基础资料。图5c和图5d分别给出了采用某商业软件得到的纵波阻抗和横波阻抗反演结果。由图5可见,本文方法具有更高的分辨率,尤其是横波阻抗的反演结果优势明显,具有更好的可解释性,有利于后续储层与地质解释。不足的是,尽管横波阻抗反演结果的改善很大,但剖面局部仍存在一些由于噪声和振幅畸变造成的不稳定现象,后续研究可以重点尝试利用多道策略加强反演结果的横向稳定性。本文方法的计算效率在可接受的范围内,但略逊于商业软件,后续研究可以进行针对性的算法优化以提高效率。

图4 具有不同入射角的叠前部分叠加地震数据a 入射角为7°; b 入射角为14°; c 入射角为21°

图5 实际地震数据采用本文方法和某商业软件得到的横波阻抗和纵波阻抗反演结果a 本文方法反演得到的纵波阻抗; b 本文方法反演得到的横波阻抗; c 某商业软件反演得到的纵波阻抗; d 某商业软件反演得到的横波阻抗

4 结论

通过引入自动相关判别先验,在频率域提出了一种新的叠前同时反演方法。自动相关判别先验假设每个反射系数都由一个待估计的参数决定,避免了对反射系数的总体分布进行假设,并保证了反演结果的较高分辨率和对岩性分界面的识别能力。同时,该参数可以随着不同地震道和弹性参数反射系数的特征自适应变化,提高了反演的准确性。反演中加入的趋势约束有效提高了反演的稳定性和对噪声的适应能力。本文方法利用优势频带获取反演结果,避免了利用信噪比低的频带信息,有利于提高反演结果的准确性。模拟数据和实际数据实验结果证实了本文方法的可行性、优势和应用潜力。需要注意的是,基于离散傅里叶变换的低频约束会在一定程度上降低反演的效率,在实际应用中,可以考虑采用频谱分析缩减反演利用的频率区间、精细划分层位控制参与反演的采样点数等方面进一步优化算法、提高效率。