一种地震共反射点道集数据的叠前反演方法

张繁昌,印海燕,翁 斌,王保丽

(1.中国石油大学 地球资源与信息学院,山东东营 257061;2.中国海洋石油研究总院,北京 100027)

一种地震共反射点道集数据的叠前反演方法

张繁昌1,印海燕2,翁 斌2,王保丽1

(1.中国石油大学 地球资源与信息学院,山东东营 257061;2.中国海洋石油研究总院,北京 100027)

叠前地震资料含有地层的纵波、横波速度和密度等信息。利用叠前反演获得隐藏在地震数据中的这些基本参数后,即可揭示大量岩性及孔隙流体性质的信息。这里推导了平面波在层状弹性介质中传播的正演算子,提出了一种基于B rent正交搜索方向组的叠前三参数反演方法,该方法不需要求解庞大而复杂的导数矩阵。通过自适应退火因子和罚函数来处理约束条件,提高了算法的稳定性。将K-L变换引入到方向置换过程,有效防止了搜索方向组的线性相关。经理论模型和油田实际数据的反演结果表明,该反演方法是一种利用叠前地震数据进行储层预测的有效手段。

叠前反演;地震道集;B rent算法;正演算子;K-L变换

0 前言

纵波的传播速度随着岩石中孔隙流体类型的变化而变化,而横波速度则主要取决于岩石骨架的矿物成份。因此,同时利用纵波、横波信息,可以有效地进行储层识别和烃类检测[1、2]。我们知道,地震波在沿着非法向路径传播时,每遇到一个地层分界面,就会产生反射纵波和反射横波,最后被地面检波器接收[3]。换句话说,叠前地震数据不仅包含有纵波和密度信息,而且携带了地层横波信息。叠前反演的目的,就是利用叠前地震数据来得到地下的纵波、横波等信息[4]。

叠前反演主要有波动方程反演和AVO反演二个方向。叠前波动方程反演在理论上虽然比较成熟,但由于其运算量极大,对噪音敏感,在实际中并没有得到应用。而基于Zoepp ritz方程的叠前道集AVO反演由于运算简单,抗噪性好而被广泛应用[5]。由于三参数AVO反演问题是病态的[6],业界普遍采用二参数方程进行反演,以达到稳定反演问题的目的,但同时也引入了较大误差[7]。另外,一般反演方法都要利用梯度信息,即在反演过程中须求解目标函数的一阶导数矩阵。由于叠前反演数据量大,待求解的地层参数较多,梯度矩阵将占用相当大的内存空间,从而限制了叠前反演的数据规模。

作者在本文中,将采用一种基于B rent正交搜索方向组[8]的方法,以实现叠前道集的三参数同时反演。该方法的优势是在反演过程中不需要导数信息。叠前反演目标函数采用一种退火罚函数方法构建,克服了三参数反演问题的病态性质,提高了该算法的稳定性。在反演迭代过程中,将KL变换[9]构造的正交方向作为新的搜索方向组,则有效地避免了搜索方向的线性相关,提高了算法的效能。

1 叠前反演问题描述

在非线性叠前地震反演中,实际地震道集dobs与由地质模型参数m合成的道集dmod=A(m),见

式(1)。

其中 矩阵A为正演算子;e为残差噪音。

可以建立式(2)的叠前反演目标函数:

其中 dobsij表示第i次覆盖、第j个样点的实际观测地震数据;dmodij表示对应位置的合成道集数据。

实际上,式(2)定义的是预测误差向量e的欧几里的长度平方,可以用任意范数来求解地层模型参数。由于合成数据与实际数据差异很大的奇异点对低阶范数影响不大,而对高阶范数影响甚大,为了减小少数奇异点对目标函数的影响,将式(2)用L1模代替,如式(3)所示:

2 模型参数与正演算子的表示

作者在进行常规地震叠后反演认为,地震波是法向入射到界面上的,但事实上常规叠加道是不同入射角(或不同偏移距)地震记录的水平叠加,不能真正代表法线入射的地震记录。所以,叠加已经破坏了真实的振幅关系。而利用叠前道集反演,要基于Zoepp ritz方程的三项线性简化方程[10、11],通常用地层弹性参数的变化率表示(见式4):

其中 Rpp(θj)为纵波反射系数,其大小随入射角θj而变化,产生AVA效应;Vp、Vs及ρ分别表示地层分界面二侧的局部平均纵波、横波速度和密度;γ=4(Vs/Vp)2;纵波、横波速度和密度的相对变化率分别用ΔVp/Vp、ΔVs/Vs、Δρ/ρ表示,这些参数的相对变化量称该参数在地层界面产生的反射系数。

由式(4)可以看出,非零偏移距地震道的反射振幅,除了与纵波速度及密度有关外,还含有横波速度信息,从而得到泊松比等其它地层弹性参数。记:

将方程(4)重新整理,并写成矩阵形式(5)。

其中

方程(5)描述的是单个地层界面,在某一入射角或偏移距下的情形。将其扩展到N个地层界面、K个偏移距的一般情形,可以用分块矩阵表示为式(6)。

式中 向量Rpk是偏移距为k的所有界面的反射系数;rp、rs、rρ分别是纵波、横波速度和密度反射系数向量;分块矩阵Ek、Fk及Gk均为对角矩阵,如

由于带限叠前地震数据是随偏移距变化的反射系数和子波的褶积,设Wk是偏移距为k时的子波矩阵,dk为相应的地震数据,则有式(7)。

作为解释人员,更加关心的是地层参数信息,而不是地层分界面的反射系数。由于反射系数和地层参数的对数值之间存在微分关系,如rp=∂lnVp/2∂t,以差分形式表示即为:

以矩阵形式表示为rp=D lnVp,其中D为差分矩阵:

类似地:

最后得到叠前道集的正演方程,如式(8)。

其中 系数矩阵即为式(1)中的正演算子A;地层参数m=[VpVsρ]T;地震道集数据为dmod=[d1,d2,…,dK]T。

在地震数据有N个时间采样点,K次覆盖的情形下,向量m有N×3个元素,向量d mod有N×K个元素,正演算子A为N×K行、N×3列的大型矩阵。

3 罚函数的构建

由于叠前地震资料普遍信噪比较低,且观测孔径不足,导致了问题的不适定性。仅依赖合成道集与实际道集的误差大小得到的地层参数是非唯一的,甚至不具有实际物理意义和地质意义,这也是造成三参数反演病态性的主要原因。为了解决这个问题,得到具有地质意义的解,所言反演的地层模型参数Vp、Vs和ρ,除了要使式(3)最小外,还要满足约束条件:

这样,引入先验信息对解空间进行约束,可以使模型参数在先验值附近进行搜索,减小寻优范围,加快收敛速度。

在引入先验信息后,无约束反演问题就变成了有约束反演问题。为了能够利用B rent方法进行求解,作者设计了一种退火罚函数法,将有约束问题转化为无约束问题。设模型参数的罚函数为:

由式(9)可以看出,当反演的模型参数值偏离先验值太多,大于门槛值η时,该约束项才起作用;而小于η时,不予惩罚。惩罚因子σ吸取模拟退火[12]的思想:令σ=T(k),T(k)是模拟退火算法中的冷却进度表,其管理方程为T(k)=αT(k-1),其中0<α<1,k为迭代次数。这样,在反演初期阶段,先验约束的作用比较大,能够很好地防止反演参数的过度调整而产生的剧烈抖动。随着反演过程的不断进行,T逐渐下降,罚函数的作用逐渐减小,最终得到具有地质意义的解。

4 反演问题的求解

在引入先验信息约束后,叠前反演所要求解的问题为式(10):

对式(10)的反演有多种方法,业界普遍使用的是基于梯度信息的反演方法[13],在反演过程中需要计算式(10)的一阶导数。由于叠前反演数据量大,待求解的地层参数较多,由一阶导数构成的梯度矩阵将占用相当大的内存。

在不使用导数信息的反演方法中,B ren t方法是一种相当有效的方法,其基本思想是用任意n个线性无关的方向ξ1、…、ξn作为初始搜索方向组,经过n次迭代,陆续构造出n个互相共扼的搜索方向组u1、…、un。在每次迭代中更新搜索方向组时,都保留已经得到的共扼方向,而用新构造出的方向替换使目标函数值下降最大的方向。

在目标函数有多个极值点的情况下,B ren t方法会出现这样的问题:在第k次迭代中,由于搜索方向的改变,使u1、…、uk变得线性相关或接近线性相关[14],从而严重限制了该方法在多极值反演问题中的应用。为避免搜索方向出现线性相关,作者在方向置换策略过程中引入K-L变换,以保证搜索方向线性无关。

在叠前反演中,设地层参数解空间的初始搜索方向组为ui(i=1,...,n),其中u1、…、uk是n维空间中k个关于A共轭的方向,令:

则有:

(1)置初始迭代次数k=0,地层参数的初始估计mk∈Rn,初始搜索方向组ui=ξi∈Rn(i=1,...,n)。

(2)记mk,0=mk为Rn空间中的初始点,按次序i=1、…、n沿各搜索方向寻优,即求 βi,使Q(mk,i-1+βiui)满足式(10),并置:

(3)找到使目标函数下降最快的p个方向,记录在集合Φ={l1,l2,…,lp}⊂{1,2,…,n}中,利用K-L变换对下标集中的方向组进行置换。

(4)置mk+1=mk,n,若满足式(10)规定的门槛值,则mk+1为最优解;否则,置k=k+1,返回步骤(2)继续迭代。

5 方法应用

5.1 方法测试

图1 无噪声时的道集对比Fig.1 Gather com parisonw ithou t no ise

模型试验所采用的叠前道集(见图1(a)),由实际井曲线根据Zoepp ritz方程(如下页图2中的虚线所示)产生,地震子波为主频35 Hz的R icker子波,共188m s,1m s采样,最小偏移距350m,最大偏移距3 200m,共二十次覆盖。将采用K-L变换进行方向置换的新算法与原算法进行了对比,搜索方向组的线性相关指标(越大表示越线性相关)如下页图3所示。由此可见,原算法中搜索方向组随着迭代次数的增加,变得越来越线性相关(如图3中实线所示),致使算法不收敛。而新的方向置换算法,有效地防止了线性相关(如图3中虚线所示)现象。

利用本文中的反演算法,对图1(a)的叠前道集进行纵横波速度及密度三参数,同时反演的结果如下页图2中的实线所示。通过与实际井曲线的对比可见,二者基本重合。需说明的是,本反演方法是按照地震数据的采样间隔进行反演的,事先并不知道地层界面的位置。反演结果合成的叠前地震道集如图1(b)所示,图1(c)为前二者的残差,可见反演结果的合成道集与真实模型参数的道集极为相似,残差基本为零。

为测试地震数据中随机噪声对反演结果的影响,在图1(a)的地震道集中,加入不同程度的随机噪声。

后面图4(a)是加入30%随机噪声时的地震道集,用本文中的方法对其进行叠前反演,反演结果如后面图5所示。

由此可见,纵波、横波曲线与实际井曲线能够很好地吻合,虽然密度曲线的质量有所下降,但仍具有较好的可比性。反演结果的合成道集如图4(b)所示(见下页),图4(c)所示(见下页)的残差数据基本上是随机噪声。

当地震道集中的随机噪声达到50%时(见后面图6(a)),有效信号已淹没在噪声中,后面图6(b)所示反演结果的合成道集与真实地层模型的地震道集(如图1(a)所示)已差别很大,后面图6(c)的残差记录中还隐约残留有效信号同相轴的影子。图6(a)地震道集的反演结果如后面图7所示。

与前面的反演结果相比,各参数与实际井曲线的偏差明显增大,但从总体上看,纵波速度与原测井曲线仍具有较好的相似性,横波和密度参数的反演质量下降较大,可见纵波速度较其它二个参数的反演具有更强的抗噪能力,原因可能是由于纵波对地震AVO响应的贡献最大,这从方程(4)中各参数前面的系数不难看出。

5.2 实际数据的应用

后面的图8(a)是丽水凹陷过L2井的叠前地震道集,其中在1.87 s附近(实心箭头所指)为气层,1.95 s(空心箭头所指)处为假亮点。L2井已有实测的纵波速度、横波速度和密度曲线,对该位置进行叠前反演,是为了检验方法的有效性,然后再对工区内没有横波速度的井位,利用本方法得到相应的横波速度。

图2 无噪声时反演结果与实际井曲线的对比Fig.2 Comparison between inverted and realwell logw ithoutnoise(Dashed line is realwell log,so lid line is inverted resu lt)

图3 线性相关指标值随迭代次数的变化Fig.3 Changesof linear-dependen t ind icato rw ith iteration

图9(见后面)从左到右是L2井位置处反演的纵波速度、横波速度及密度曲线,与该井实际测井曲线的对比,其中实线为反演曲线,虚线为L2的实测井曲线。

图8(b)(见后面)是由反演结果生成的道集,与图8(a)的同相轴位置和振幅变化关系都很一致。

图4 含较小噪声时的道集对比Fig.4 Gather comparisonw ith less noise

图5 图4(a)所示地震道集的反演结果与实际井曲线的对比Fig.5 Comparison between inverted resultof Fig 4(a)and realwell log

图6 含较大噪声时的道集对比Fig.6 Gather comparisonw ith large noise

图7 图6(a)所示地震道集的反演结果与实际井曲线的对比Fig.7 Comparison between inverted resultof Fig 6(a)and realwell log

见下页,综合图8和图9可以看出,在1.87 s处,纵波速度降低,密度明显减小,而横波速度变化微弱,因此在图8中实心箭头所指的同相轴是由储层含气所致。而在1.95 s处的纵波、横波速度均减小,密度值却较大,故在图8中空心箭头所指的同相轴是由岩性变化引起的,并不是由于储层含气引起的。

6 结论

图8 L2井反演合成道集与实际地震道集的对比Fig.8 Comparison between synthetic and real gatherofwellL2

图9 L2井三参数反演结果与实际井曲线的对比Fig.9 Comparison between 3-term inverted resultand wellL2 real log

叠前地震共反射点道集中含有地层的纵波速度、横波速度和密度等信息,在利用叠前反演获得这些参数后,可揭示岩性及储层所含流体性质的信息。因此,叠前反演已成为岩性油气藏勘探和烃类检测的有力手段。

作者在本文中详细推导了平面波在层状弹性介质中传播的正演算子,进而提出了基于B rent正交搜索方向组的叠前三参数反演方法,该方法不需要求解和存储复杂的导数矩阵。通过自适应退火因子和罚函数来构建约束条件,克服了三参数反演的病态性质,提高了算法的稳定性。方向置换则利用K-L变换实现,有效地防止了搜索方向组的线性相关现象。经理论模型检验和实际数据的应用表明,本方法反演过程稳定,效果良好,是一种利用叠前道集数据反演地层纵波、横波及密度信息的有效手段。

[1]ANDERSON P F,GRAY FD.U sing LMR for dual attribute litho logy identification[J].71 th Ann.Internat M tg.,SEG Expanded Abstracts,2001:201.

[2]周中彪.基于岩石物理模型的测井约束横波速度计算方法研究[J].物探化探计算技术,2010,32(5):536.

[3]印兴耀,韩文功,李振春,等.地震技术新进展(下)[M].山东东营:中国石油大学出版社,2006.

[4]肖思和,李曙光,许多,等.叠前弹性波阻抗反演在储层预测中的应用[J].物探化探计算技术,2010,32(5):476.

[5]王永刚.地震资料综合解释方法[M].山东东营:中国石油大学出版社,2007.

[6]DEBSK IW,TARANTOLA A.Inform ation on elastic param eters obtained from the amp litudes of reflected waves[J].Geophysics,1995,60(5):1426.

[7]陈建江,印兴耀.基于贝叶斯理论的AVO三参数波形反演[J].地球物理学报,2007,50(4):1251.

[8]BRENTR P.A lgorithm form inim izationw ithoutderivatives[M].Englewood C liffs:Prentice-Hall,1973.

[9]张贤达.矩阵分析与应用[M].北京:清华大学出版社,2004.

[10]AK IK,R ICHARDSPG.Quantitative seismo logy:theory andm ethods[M].Cam bridge:W.H.Freem an and Co.,1980.

[11]CONNOLLY P.Elastic impedance[J].The Leading Edge,1999,18(4):438.

[12]PEID,LOU IE JN,SATISH K P.App lication of sim ulated annealing inversion on high-frequency fundam ental-mode Rayleigh wave dispersion curves[J].Geophysics,2007,72(5):R77.

[13]施光燕,董加礼.最优化方法[M].北京:高等教育出版社,2005.

[14]高旅端.正交程度及其在B rent方法中的应用[J].北京工业大学学报,1997,23(2):42.

P 631.4

A

1001—1749(2011)01—0011—09

国家“863”项目资助(2006AA 09A 102)

2010-06-17 改回日期:2010-11-02

张繁昌(1972-),男,副教授,博士,主要从事地震反演、储层解释方面的研究。