地球物理学报  2017, Vol. 60 Issue (1): 271-282   PDF    
纵向和横向同时约束AVO反演
霍国栋1 , 杜启振1 , 王秀玲2 , 陈刚3 , 郑笑雪1     
1. 中国石油大学(华东)地球科学与技术学院, 青岛 266580;
2. 中国石化石油物探技术研究院, 南京 211103;
3. 中石油新疆油田公司勘探开发研究院, 克拉玛依 834000
摘要: 为了提高AVO(amplitude versus offset)反演结果的精度和横向连续性,本文提出了一种新的AVO反演约束方法,该方法结合贝叶斯原理和卡尔曼滤波算法实现了对反演参数纵向和横向的同时约束.文章首先结合反演参数的纵向贝叶斯先验概率约束和反演参数的横向连续性假设建立了与卡尔曼滤波算法对应的AVO反演系统的数学模型,然后将该数学模型代入卡尔曼滤波算法框架,利用卡尔曼滤波算法实现了双向约束AVO反演.二维模型测试和实际数据测试结果表明,相对于单纯的纵向贝叶斯先验概率约束,双向约束能更准确地刻画参数的横向变化,得到更准确、横向连续性更好的反演结果.
关键词: 贝叶斯先验概率约束      卡尔曼滤波算法      双向约束      横向连续性     
AVO inversion constrained simultaneously in vertical and lateral directions
HUO Guo-Dong1, DU Qi-Zhen1, WANG Xiu-Ling2, CHEN Gang3, ZHENG Xiao-Xue1     
1. School of Geosciences, China University of Petroleum(East China), Qingdao 266580, China;
2. Sinopec Geophysical Research Institute, Nanjing 211103, China;
3. Exploration and Development Research Institution, PetroChina Xinjiang Oilfield Company, Karamay 834000, China
Abstract: The prior probability distribution constraint, which is based on the Bayesian principle, is widely used in AVO (amplitude versus offset) inversion to make the inversion process well-posed. It assumes that the parameters of each trace obey some kind of prior probability distribution in vertical direction. And then the Bayesian principle is used to build the posterior probability density function and to maximize it to obtain the best inversion results in the probability theory sense. However, we cannot obtain satisfying inversion results just depending on the Bayesian constraint. We think that one of the reasons for this less-than-ideal inversion results is that the Bayesian constraint works only in vertical direction but ignores the lateral continuity of parameters. To increase the accuracy and lateral continuity of inversion results, a new AVO inversion method is proposed to constrain the parameters in both vertical and lateral directions by combining the Bayesian principle and Kalman filtering algorithm. We assume that the parameters of two neighboring seismic traces are similar for earth media change gradually in the lateral direction. Based on this hypothesis, we think that the parameters of a seismic trace are determined by two parts:one is the information of parameters contained in the seismic record of that trace, and the other is the similarity with its neighboring trace. Building a mathematical model of the two parts above permits to obtain the measuring equation and the state transition equation. The measuring equation, which is established according to the Bayesian principle, is the maximum posterior probability density function of parameters. The state transition equation and the measuring equation are substituted into the Kalman filtering algorithm frame, and then the Kalman filtering algorithm is used to predict and modify the parameters of each trace. The parameters after predicting and modifying are the final double-direction constrained AVO inversion results. The Marmousi2 model and actual data are used to test the proposed method. The inversion results show that the double-direction constrained method has better lateral continuity as well as accuracy than that with only vertical direction Bayesian constraint. And the double-direction constrained method is based on prior hypothesis, which does not depend on well logging or horizon information. So this method can be widely used even in the area without well logging or horizon information..
Key words: Bayesian prior probability constraint      Kalman filtering algorithm      Double-direction constraint      Lateral continuity     
1 引言

AVO反演主要利用Zoeppritz方程的近似式研究地震反射系数与地下介质参数之间的关系,以褶积模型作为联系反射系数与地震记录的纽带,并在褶积模型的基础上添加测井信息、地质信息以及一些先验信息以约束介质参数,最后形成AVO反演的数学模型.基于上述模型,利用线性、非线性以及概率统计等数学方法从地震记录中提取介质参数(刘洋等,2012; 彭真明等,2008; 邴萍萍等,2012; 陈建江和印兴耀,2007).然后以岩石物理为桥梁将所得参数转化为流体指示因子,用于储层预测与流体识别.作为一种有效的叠前弹性参数反演方法,AVO反演一直是研究热点且被广泛用于实际生产.

AVO反演技术虽然已经相当成熟,但它仍存在一些局限,其中的一个主要原因在于AVO反演是一个病态问题(马劲风等,2000; 闫惠中,2012; 印兴耀等,2014),对反演参数加一些约束是降低其病态程度的重要手段(马劲风等,2000).利用测井信息以及其他地质信息直接对参数做一定程度上的约束能够得到与地质资料较为吻合的反演结果(马劲风等,2000; 张铭等,2012; 田军等,2013).基于概率统计的约束是另一种重要的手段,其中基于贝叶斯理论的先验概率约束在AVO反演中应用广泛(Buland and Omre,2003; 陈建江和印兴耀,2007; Alemie and Sacchi,2011; 侯栋甲等,2014),它假设弹性参数在纵向上服从特定的先验概率分布,通过建立参数的后验概率密度函数并最大化之得到概率论意义下的最优反演结果,以此降低反演的病态程度.但仅通过贝叶斯先验概率约束还不能达到很好的反演效果,首先,反演结果与真实值相差比较大,Buland和Omre(2003)指出仅依靠贝叶斯软约束只能得到较为准确的纵波阻抗,而几乎不能提供密度信息,且反演的不确定性较高.田军等(2013)指出单纯地利用贝叶斯概率化方法做AVO反演,所得结果是带限的,影响了反演结果的相对可靠性.其次,单纯的贝叶斯约束(无低频模型约束)所得反演结果在横向上不具备良好的连续性.通过分析认为造成反演结果不理想的原因之一在于贝叶斯约束只在单个地震道上对反演参数做纵向上的约束,而忽略了参数在横向上的关联性.

为了提高AVO反演结果的精度和横向连续性,本文结合贝叶斯原理和卡尔曼滤波算法提出了一种新的AVO反演约束方法,该方法能在纵向和横向上同时对参数做约束.由于地球介质是连续分布的,在横向上具有渐变性,因此本文假设相邻地震道的弹性参数具有一定的相似性.基于该假设,本文认为一个地震道的弹性参数由两部分共同决定,一是本道的地震记录所包含的弹性参数的信息,二是本道与其相邻道的弹性参数所具备的相似性.分别建立上述两个部分的数学模型得到测量方程和状态转换方程,其中的测量方程为依据贝叶斯原理建立的参数的最大后验概率密度函数,以上两个方程即为与卡尔曼滤波算法对应的AVO反演系统的数学模型.然后将状态转换方程和测量方程代入卡尔曼滤波算法框架,利用卡尔曼滤波算法(Kalman,1960; 彭丁聪,2009)对每一地震道的弹性参数进行预测和修正两个步骤,最后得到每一个地震道最优的弹性参数.上述方法在不加入测井、层位约束条件的情况下,在纵向和横向上同时对参数做约束,以提高反演结果的精度和横向连续性.

2 纵向贝叶斯先验概率约束

纵向约束即假设参数在一个地震道上服从特定的先验概率分布,以此降低反演的病态程度,使一个地震道的参数反演趋于稳定.该部分首先建立AVO褶积正演模型,然后以建立的正演模型为基础,依据贝叶斯原理推导反演参数的最大后验概率密度函数,以此为卡尔曼滤波双向约束AVO反演做准备.

2.1 建立正演模型

本文采用Fatti反射系数近似公式(Fatti et al.,1994)做三参数AVO反演,Fatti公式如(1)式:

(1)

其中A=1+tan2θB=-8k2sin2θ,C=4k2sin2θ-tan2θk为横纵波速度比,θ为入射角.分别为纵、横波阻抗反射率和密度反射率,IP=vPρIS=vSρ.如果θ取为一个定值,则(1)式表达的仅是一个入射角、一个反射界面的情况.为了构建一个CMP道集的AVO正演模型,将(1)式扩展为M个入射角、N个反射界面的情况并表达为矩阵形式,可得:

(2)

其中Ri以及rPrSrdN维列向量,Ri为第i个入射角所有反射界面的反射系数,rPrSrd分别为所有反射界面的纵、横波阻抗反射率和密度反射率.AiBiCiN×N对角系数矩阵(文中粗黑符号代表矩阵或向量).利用褶积模型将反射系数转化为地震记录:

(3)

(3)式中di为第i个入射角的N维地震记录列向量,Wi为第i个入射角的N×N子波矩阵,ni为与di对应的N维噪声向量.将(2)式代入(3)式得:

(4)

将(4)式简记为:

(5)

(5)式即为将子波包含在内的一个CMP道集的AVO褶积正演模型,它刻画了地震记录d与弹性参数r之间的关系.(5)式中d为CMP道集数据向量,G为正演算子,r为参数向量,n为噪声向量.

2.2 贝叶斯原理

贝叶斯原理的核心思想体现在如下的贝叶斯公式:

(6)

其中P(r|d)r的后验概率密度函数,P(d|r)d的似然函数,P(r)r的先验概率分布函数,P(d)为边缘概率,其为一常数,因此(6)式等价于:

(7)

根据(7)式可知贝叶斯原理为:参数的后验概率密度函数由观测数据的似然函数和参数的先验概率分布函数共同决定.具体到本文的三参数AVO反演中即为:纵、横波阻抗反射率和密度反射率的后验概率密度函数由地震记录的似然函数以及三个参数的先验概率分布函数共同决定.

2.3 建立似然函数

通过假设地震记录中的噪声满足特定的概率分布能将地震记录和反演参数联系起来,本文假设各入射角和反射界面的地震噪声满足独立同高斯分布(Buland and Omre,2003),考虑所有入射角和反射界面可得噪声的联合概率密度函数,即似然函数的表达式为:

(8)

其中:ni为第i个采样点的噪声,μn为噪声的均值,σn为噪声的标准差.假设噪声均值为零,并将(5)式代入(8)式可得:

(9)

(9) 式即为地震记录的似然函数,利用极大似然参数估计方法可以由(9)式导出参数的最小二乘解.

2.4 反演参数的先验概率约束

贝叶斯约束的优越之处在于它在(9)式的基础上加入了对参数的先验概率约束,虽然加入的先验概率分布不一定符合实际,但是会在一定程度上降低问题的病态程度,这样将比直接的最小二乘反演更稳定可靠.本文假设反演参数在纵向上服从三元高斯分布(Buland and Omre,2003),则第i个采样点的参数的概率密度函数如下:

(10)

其中ri(3×1)为第i个采样点处的参数向量,它的三个分量分别为纵、横波阻抗反射率和密度反射率.Σ(3×3)为三个分量之间的协方差矩阵,μr(3×1)为所有ri的均值向量.考虑一个地震道上的所有N个采样点就得到参数的联合概率密度函数为:

(11)

为了将(11)式中的ri表示成r的形式,我们引入矩阵si从整个r中选取第i个采样点处的参数向量ri,即ri=sir.si是一个3×3N的矩阵,其元素满足(Alemie and Sacchi,2011):

(12)

ri=sir代入(11)式可得:

(13)

其中.(13)式即为反演参数所满足的先验三元高斯联合概率密度函数.

2.5 最大后验概率密度函数

根据贝叶斯原理将(9)式和(13)式代入(7)式得到后验概率密度函数:

(14)

使(14)式最大即可得到最大后验概率密度函数.由(14)式可知,要使其最大等价于使(15)式最小,即(15)式为目标函数.

(15)

对(15)式关于ri(i=1,2,…,3N)求导并令求导结果为零,得:

(16)

(16)式即为最大后验概率密度函数,求解(16)式可得到最大后验概率解,即最终的贝叶斯约束反演结果.从(16)式可见贝叶斯约束项的权重为观测数据噪声的方差(σn2),如果观测数据中没有噪声(σn2=0),贝叶斯方法就退化为最小二乘法.因此,贝叶斯约束项的权重应该根据地震记录的信噪比调整,当信噪比高时减小权重,反之增大权重.

3 基于卡尔曼滤波算法的双向约束AVO反演

纵向上的贝叶斯先验概率约束能够使单个地震道的反演趋于稳定,但各地震道之间相互独立、缺少必要的联系.由于地球介质在横向上具有渐变性,特别是层状介质具有较强的横向连续性.因此,本文假设相邻两个地震道的弹性参数具有一定的相似性,并建立相应的数学模型.将该数学模型和贝叶斯最大后验概率密度函数(16)式代入卡尔曼滤波算法框架,利用卡尔曼滤波算法实现对参数纵向和横向的同时约束,以此反演得到的参数受纵向先验概率分布以及横向连续性假设的双重约束.

3.1 建立AVO反演系统的数学模型

要利用卡尔曼滤波算法实现地震剖面的双向约束反演,需要建立与卡尔曼滤波算法对应的AVO反演系统的数学模型,该数学模型包括两个方程:状态转换方程和测量方程.

首先,我们将各地震道之间的参数变化看作是状态随空间的转换,相邻道之间状态的转换由如下的状态转换方程描述:

(17)

其中的r对应于(5)式中的rrk+1rk分别代表地震剖面中第k+1道和第k道的反演参数,Ak是状态变换矩阵,即联系第k+1道和第k道反演参数的中间量.本文假设相邻地震道的参数具有相似性,因此对于所有k值,Ak均取为单位矩阵以刻画参数的横向连续性,所以(17)式就变为:

(18)

(18) 式即为AVO反演系统的状态转换方程,由(18)式可见第k+1道的反演参数可由第k道的反演参数预测.由于相邻两道的反演参数不可能完全相等,所以在转换过程中存在转换误差wkwk是与r维数相同的列向量.wk的大小由地质体的横向变化剧烈程度决定,横向变化越剧烈wk越大,反之wk越小.wk服从零均值多元高斯分布,其协方差矩阵为Q,即

(19)

Q为状态转换误差协方差矩阵.由上述分析可知Q反映了地质体横向变化的剧烈程度.由于地质体横向变化的剧烈程度需要根据反演结果确定,因此Q是一个先验量,它是双向约束AVO反演的输入,本文假设Q不随k变化而变化.

其次,一个地震道的参数不仅与其相邻道的参数有关,还与实测地震数据有关.把地震数据与参数之间的关系式称为测量方程,测量方程最简单的形式为前文构建的正演模型(5)式,利用测量方程我们可以对状态转换方程预测的参数做修正,使得最终反演所得参数更加准确.由于(5)式中没有任何约束,所以直接利用(5)式不能得到理想的修正效果,这一点在后文有相应的分析说明.因此,为了提高修正效果前文应用贝叶斯原理在(5)式的基础上对参数施加了先验概率约束,得到了最大后验概率密度函数(16)式,在(16)式右边加一个误差项vk就得到了测量方程(20)式:

(20)

vk体现了利用(20)式反演参数的精度,反演精度越低vk越大,反之vk越小.vk服从零均值多元高斯分布,其协方差矩阵为R,即

(21)

R为测量误差协方差矩阵,本文假设R不随k变化而变化.和状态转换误差协方差矩阵Q一样,R也是一个反演的输入量,由于R体现了利用(20)式反演参数的精度,而影响(20)式反演精度的因素主要包括(20)式所用数学模型与实际情况的吻合程度、(20)式的病态程度以及地震数据的质量.因此,实际情况下R的确定要综合考虑上述几个因素.将(20)式简记为:

(22)

其中:.至此,与卡尔曼滤波算法对应的AVO反演系统的数学模型已经建立,即:状态转换方程(18)式和测量方程(22)式.其中状态转换方程体现了横向约束,而测量方程体现了纵向先验概率约束,下文将把这两个方程代入卡尔曼滤波算法框架,实现双向约束AVO反演.

3.2 利用卡尔曼滤波算法实现双向约束AVO反演

卡尔曼滤波算法是一种自回归算法,它利用状态转换方程和测量方程对参数进行预测和修正两大步骤,通过这两个步骤得到每一个状态的最优结果,这里的状态在本文中即为地震道的参数.

首先利用状态转换方程(18)式由一个地震道的参数预测与之相邻的下一个地震道的参数.假设第k道的参数已经得到,现在要反演第k+1道的参数.根据状态转换方程,第k+1道的参数可由第k道的参数预测为:

(23)

其中rk+1,k是由第k道的参数预测的第k+1道的参数,rk,k是已反演得到的第k道的最优结果.由(23)式可见,第k+1道的预测值就是第k道的最优结果,这样粗糙的预测值必然存在误差,该预测误差由预测误差协方差矩阵Pk+1,k来刻画:

(24)

其中Pk+1,krk+1,k的预测误差协方差矩阵,Pk,krk,k的误差协方差矩阵,Q是状态转换误差协方差矩阵.由(24)式可见,第k+1道的预测误差有两个来源,一是Pk,k:第k道的最优结果rk,k中包含的误差,二是Q:状态转换过程中引入的误差.

完成预测后就有了第k+1道的预测值rk+1,k及相应的预测误差Pk+1,k,现在利用实测地震数据和测量方程(22)式对第k+1道的预测值做修正,修正后的最优结果为:

(25)

rk+1,k+1即为第k+1道最终的双向约束AVO反演结果,其中Kg为卡尔曼增益(Kalman Gain),其作用是使最优结果rk+1,k+1和真实值之间的误差平方和最小,Kg的表达式为:

(26)

其中R为测量误差协方差矩阵.

到现在为止,已经完成了预测和修正两步,并得到了第k+1道的最优结果rk+1,k+1.但是要计算rk+1,k+1的误差协方差矩阵Pk+1,k+1,以此将第k+1道的最优结果中存在的误差过渡到第k+2道,Pk+1,k+1的表达式为:

(27)

其中I为单位矩阵.这样,当系统进入第k+2道时,Pk+1,k+1就对应于(24)式中的Pk,k卡尔曼滤波算法就能自回归地运行下去,直到完成整个地震剖面的反演.

通过上面的分析可见,地震剖面中后面道的反演考虑了前面道反演结果中存在的误差,因此不会造成横向上的误差累积.(23)—(27)式即为实现双向约束AVO反演的五个核心公式,这五个公式以卡尔曼滤波算法为框架,并将参数的纵向贝叶斯先验概率约束和参数的横向连续性约束包含入内,从而实现对参数纵向和横向的双重约束.双向约束AVO反演的卡尔曼滤波算法可表示为如下流程图 1.

图 1 双向约束AVO反演方法的卡尔曼滤波递归算法流程图 Fig. 1 Flow chart of Kalman filtering recursive algorithm for double-direction constrained AVO inversion

对于该算法的实现有三点需要说明.首先,算法的启动需要初始地震道的参数向量r0,以及与之对应的误差协方差矩阵P0,这两个初始值不需要很准确,因为卡尔曼滤波算法对它们不敏感.其次,本文的横向约束和纵向约束都基于先验假设,其约束力度分别由QR控制.因此,QR的选取应该综合考虑地质体的横向变化剧烈程度以及纵向上贝叶斯约束的准确程度,这样才能使纵向约束和横向约束达到平衡,使反演结果既准确又具有良好的横向连续性.最后,反演开始的地震道可以是剖面中的任意一道,如果是剖面两端的地震道即从一端反演到另一端,如果是剖面中间的某一道即从那一道向两边反演直到整个剖面反演完毕.但是,以不同的地震道作为反演起始道能得到非常相似的反演结果.

利用上述方法反演得到参数向量r后,根据道积分公式(28)将每一道的纵、横波阻抗反射率和密度反射率转换为纵、横波阻抗和密度.

(28)

4 模型测试

本文利用Marmousi2模型左上角的一部分对双向约束AVO反演方法做分析测试.图 2是Marmousi2模型的纵、横波阻抗和密度剖面,从图 2可见该模型具有一定的横向连续性但并非标准的水平层状介质,模型中部有一个气藏,可将其视为一个突变,因此通过该模型能够较好地测试双向约束AVO反演算法的反演效果.由于Marmousi2模型给出的是深度域的纵、横波阻抗和密度剖面,因此本文的模型测试部分在深度域进行.由于AVO反演是由一个角道集反演得到一道的参数,因此Marmousi2模型的正演数据体是一个三维数据体.为了便于显示,从正演数据体中抽取300CDP处过气藏的深度域角道集显示如图 3图 3a为无噪声数据,图 3b为加10%高斯噪声的数据.图 3所示角道集是根据褶积正演模型(5)式正演得到,子波主频为35 Hz,入射角范围为0°~30°,角度增量为1.5°.

图 2 Marmousi2模型的纵、横波阻抗和密度剖面 Fig. 2 P-wave impedance,S-wave impedance and density sections of the Marmousi2 model
图 3 300CDP处的合成角道集 (a)无噪;(b)10%高斯噪声. Fig. 3 Synthetic angle gathers of 300CDP (a)Without noise;(b)With 10% Gauss noise.

利用不同反演方法对含噪声正演数据做反演,反演中涉及的均值向量μr、协方差矩阵Σ、状态转换误差协方差矩阵Q以及测量误差协方差矩阵R可由模型参数统计得到,这样能使每种方法均达到理想的反演效果,便于不同方法之间的对比.

图 4是纵向贝叶斯约束方法的反演结果,从左到右依次是纵波阻抗、横波阻抗和密度剖面,图 5图 4所示反演结果与模型之间的差值.图 6是利用卡尔曼滤波算法反演所得结果,但其中的测量方程采用的是正演模型(5)式,而非加入了贝叶斯先验概率约束的(22)式,即只加入了横向约束而没加入纵向约束,图 7图 6所示反演结果与模型之间的差值.图 8是双向约束方法的反演结果,图 9图 8所示反演结果与模型之间的差值.将图 2图 4图 8在300CDP处的参数抽取出来做单道对比,得到图 10.将图 10中的双向约束反演结果和纵向贝叶斯约束反演结果分别与模型值做差得到二者的误差曲线对比图,如图 11所示.

图 4 纵向贝叶斯约束方法的反演结果 Fig. 4 Inversion results of vertical direction Bayesian constrained method
图 5 纵向贝叶斯约束方法的反演结果与模型的差值 Fig. 5 Differences between the inversion results of vertical direction Bayesian constrained method and the Marmousi2 model
图 6 横向卡尔曼滤波约束方法的反演结果 Fig. 6 Inversion results of lateral direction Kalman filtering constrained method
图 7 横向卡尔曼滤波约束方法的反演结果与模型的差值 Fig. 7 Differences between the inversion results of lateral direction Kalman filtering constrained method and the Marmousi2 model
图 8 双向约束方法的反演结果 Fig. 8 Inversion results of double-direction constrained method
图 9 双向约束方法的反演结果与模型的差值 Fig. 9 Differences between the inversion results of double-direction constrained method and the Marmousi2 model
图 10 模型值、双向约束反演结果、纵向贝叶斯约束反演结果在300CDP处的单道对比 Fig. 10 Single trace comparison of Marmousi2 model,inversion results of double-direction constrained method and inversion results of vertical direction Bayesian constrained method at CDP number of 300. Red,green and blue lines respectively represent the Marmousi2 model,inversion results of double directions constrained method and inversion results of vertical direction Bayesian constrained method
图 11 300CDP处双向约束反演结果和纵向贝叶斯约束反演结果与模型的误差曲线对比 Fig. 11 Comparison of the error curves at CDP number of 300. Green curves are the error curves between the double-direction constrained inversion results and the model value. Blue curves are the error curves between the vertical direction Bayesian constrained inversion results and the model value

图 2为标准对比图 4图 8可见,图 4所示反演结果道与道之间相互独立,横向连续性较差,其中的气藏形状较为模糊,而图 8所示的三个剖面较图 4具有更好的横向连续性,其中的气藏形状也更清晰.对比图 5图 9可知,图 8所示反演结果比图 4所示反演结果具有更高的准确度.特别是密度项的准确度有较大的提高,图 4的密度项与模型值相差较大,而图 8的密度项与模型值相差很小.由图 10所示的在300CDP处的单道对比可见,双向约束反演结果比单纯的纵向贝叶斯约束反演结果更接近模型值,从图 11所示的误差曲线能更清楚地看到双向约束反演结果与模型值的差值在零值附近扰动,而单纯的纵向贝叶斯约束反演结果与模型的差值则远离零值,且随着深度的增加误差变大.上述对比结果表明,本文提出的双向约束方法能够得到比纵向贝叶斯约束方法更准确、横向连续性更好的反演结果,且能更准确地刻画地质体横向上的突变.

图 2为标准对比图 6图 8可见,图 6图 8具有比较接近的横向连续性.对比图 7图 9可知,图 8所示反演结果比图 6所示反演结果更准确.上述对比结果表明,双向约束方法能够得到与横向约束方法横向连续性相当但准确度更高的反演结果.双向约束方法能够得到更准确反演结果的原因体现在数学上即纵向贝叶斯约束的加入使得卡尔曼增益Kg的计算变得更稳定.纵向上不加贝叶斯约束时测量方程为(5)式,相应的Kg表达式为:

(29)

其中U为地震数据噪声的协方差矩阵.加入贝叶斯约束后测量方程为(22)式,相应的Kg表达式为(26)式.两种情况下Kg计算的稳定性主要取决于(Gk+1Pk+1,kGk+1T+U)-1和(Bk+1Pk+1,kBk+1T+R)-1计算的稳定性,而它们又取决于Gk+1Pk+1,kGk+1T+UBk+1Pk+1,kBk+1T+R的条件数,二者的条件数如图 12所示.

图 12 条件数对比 Fig. 12 Comparison of condition numbers

图 12可见,单纯的横向约束时其Gk+1Pk+1,kGk+1T+U的条件数远高于双向约束时Bk+1Pk+1,kBk+1T+R的条件数,因此双向约束的Kg比横向约束的Kg计算更稳定且耗时更少,所以双向约束反演结果比单纯的横向约束反演结果更准确.

5 实际数据测试

对某工区一测线的实际叠前CDP道集做超道集处理以提高信噪比,计算超道集时滚动窗口内包含6个相邻的CDP道集.再将得到的超道集转化为角道集,为了便于显示将每个CDP处的角道集做叠加得到如图 13所示的叠后道集,共131道,CDP间距为20 m,时长280 ms.图 12中的每一道都对应于一个角道集,各角道集内均有5道,角度范围为3°~25°,角度增量为5.5°,如图 14所示为CDP号为65的角道集.利用统计方法从角道集中提取的子波如图 15所示,为时间长度100 ms的零相位子波.

图 13 利用角道集叠加得到的剖面 Fig. 13 Section generated by stacking angle gathers
图 14 65CDP处的角道集 Fig. 14 Angle gathers at CDP number of 65
图 15 从角道集提取的地震子波 Fig. 15 Seismic wavelets extracted from angle gathers

以上述角道集和地震子波作为反演的输入,分别测试纵向贝叶斯先验约束方法和双向约束方法,得到图 16图 17所示反演结果,从左到右依次为纵波阻抗、横波阻抗和密度.

图 16 纵向贝叶斯约束方法的反演结果 Fig. 16 Inversion results of vertical direction Bayesian constrained method
图 17 双向约束方法的反演结果 Fig. 17 Inversion results of double-direction constrained method

对比图 16图 17可见,贝叶斯约束反演结果和双向约束反演结果具有比较接近的数值范围.但是双向约束反演结果较贝叶斯约束反演结果更稳定,横向连续性更好,特别是密度的横向连续性得到了较大提高.

6 结论

本文提出的双向约束AVO反演方法以参数的先验概率分布和横向连续性假设为基础,结合贝叶斯原理和卡尔曼滤波算法在纵向和横向上同时对参数做约束.相比于单纯的纵向贝叶斯先验概率约束,双向约束提高了反演结果的精度和横向连续性,能更准确地刻画地质体的横向变化,增强了油气指示性.

由于该双向约束反演方法以先验假设作为基础,因此它不依赖于测井、层位约束信息,具有较广的应用范围.在缺乏上述信息难以构建参数的准确低频模型的地区也能从纵向和横向两个方面去约束反演参数,得到更可靠、横向连续性更好的反演结果.

致谢

感谢Gary S. Martin提供Marmousi2模型,感谢HRS公司提供实际地震资料,感谢评审专家和本文责编反馈的意见和建议,本文曾得到中国石油大学(北京)硕士研究生刘洪星的帮助,谨此致谢.

参考文献
Alemie W, Sacchi M D. 2011. High-resolution three-term AVO inversion by means of a Trivariate Cauchy probability distribution. Geophysics, 76(3): R43-R55. DOI:10.1190/1.3554627
Bing P P, Zao S Y, Lu J T. 2012. Non-linear AVO inversion based on support vector machine. Chinese J. Geophys. (in Chinese), 55(3): 1025-1032. DOI:10.6038/j.issn.0001-5733.2012.03.033
Buland A, Omre H. 2003. Bayesian linearized AVO inversion. Geophysics, 68(1): 185-198. DOI:10.1190/1.1543206
Chen J J, Yin X Y. 2007. Three-parameter AVO waveform inversion based on Bayesian theorem. Chinese J. Geophys. (in Chinese), 50(4): 1251-1260.
Fatti J L, Smith G C, Vail P J, et al. 1994. Detection of gas in sandstone reservoirs using AVO analysis:A 3-D seismic case history using the Geostack technique. Geophysics, 59(9): 1362-1376. DOI:10.1190/1.1443695
Hou D J, Liu Y, Hu G Q, et al. 2014. Prestack multiwave joint inversion for elastic moduli based on Bayesian theory. Chinese J. Geophys. (in Chinese), 57(4): 1251-1264. DOI:10.6038/cjg20140422
Kalman R E. 1960. A new approach to linear filtering and prediction problems. Journal of Basic Engineering, 82(1): 35-45. DOI:10.1115/1.3662552
Liu Y, Zhang J S, Hu G M, et al. 2012. Study of three-term non-Gaussian pre-stack inversion method. Chinese J. Geophys. (in Chinese), 55(1): 269-276. DOI:10.6038/j.issn.0001-5733.2012.01.026
Ma J F, Wang X J, Jia C H, et al. 2000. Study of constraint methodology in constrained impedance inversion. Geophysical Prospecting for Petroleum (in Chinese), 39(2): 52-63.
Peng D C. 2009. Basic principle and application of Kalman filter. Software Guide (in Chinese), 8(11): 32-34.
Peng Z M, Li Y L, Wei W G, et al. 2008. Nonlinear AVO inversion using particle filter. Chinese J. Geophys. (in Chinese), 51(4): 1218-1225.
Tian J, Wu G C, Zong Z Y. 2013. Robust three-term AVO inversion and uncertainty analysis. Oil Geophysical Prospecting (in Chinese), 48(3): 443-449.
Yan H Z. 2012. PP and PS waves joint AVO inversion and fluid prediction[Master's thesis] (in Chinese). Qingdao:China University of Petroleum (East China).
Yin X Y, Cao D P, Wang B L, et al. 2014. Research progress of fluid discrimination with prestack seismic inversion. Oil Geophysical Prospecting (in Chinese), 49(1): 22-34.
Zhang M, Zhang W Q, Xu H X, et al. 2012. Application of AVO constrained sparse spike inversion of reservoir prediction in Pucheng oilfield. Journal of Oil Gas Technology (in Chinese), 34(4): 78-82.
邴萍萍, 曹思远, 路交通. 2012. 基于支持向量机的非线性AVO反演. 地球物理学报, 55(3): 1025–1032. DOI:10.6038/j.issn.0001-5733.2012.03.033
陈建江, 印兴耀. 2007. 基于贝叶斯理论的AVO三参数波形反演. 地球物理学报, 50(4): 1251–1260.
侯栋甲, 刘洋, 胡国庆, 等. 2014. 基于贝叶斯理论的叠前多波联合反演弹性模量方法. 地球物理学报, 57(4): 1251–1264. DOI:10.6038/cjg20140422
刘洋, 张家树, 胡光岷, 等. 2012. 叠前三参数非高斯反演方法研究. 地球物理学报, 55(1): 269–276. DOI:10.6038/j.issn.0001-5733.2012.01.026
马劲风, 王学军, 贾春环, 等. 2000. 波阻抗约束反演中的约束方法研究. 石油物探, 39(2): 52–63.
彭丁聪. 2009. 卡尔曼滤波的基本原理及应用. 软件导刊, 8(11): 32–34.
彭真明, 李亚林, 魏文阁, 等. 2008. 粒子滤波非线性AVO反演方法. 地球物理学报, 51(4): 1218–1225.
田军, 吴国忱, 宗兆云. 2013. 鲁棒性AVO三参数反演方法及不确定性分析. 石油地球物理勘探, 48(3): 443–449.
闫惠中. 2012. 多波联合AVO反演及流体预测[硕士论文]. 青岛:中国石油大学(华东).
印兴耀, 曹丹平, 王保丽, 等. 2014. 基于叠前地震反演的流体识别方法研究进展. 石油地球物理勘探, 49(1): 22–34.
张铭, 张文起, 徐海霞, 等. 2012. AVO约束稀疏脉冲反演技术在濮城油田储层预测中的应用. 石油天然气学报, 34(4): 78–82.