空间变量结构约束下AVO贝叶斯反演

1、空间变量结构约束下AVO贝叶斯反演摘要引入来源不同的地下先验信息约束地球物理反演，降低反演解的不确定性程度是地球物理反演研究中一个重要的课题。地质变量是在空间(时间)上展布的变量。地质变量除了具有统计意义上的随机性外，还依赖于空间(时间)位置,具有一定的相关性和结构性。本文引入地质统计学里变差函数的思想，建立关于地下空间变量结构的描述作为先验信息，结合贝叶斯推断理论，通过似然函数建立空间变量和地球物理观测的联系，得到AVO反演解估计后验概率密度的解析形式。选取东海某油田一口实际测井数据进行反演试算应用区间估计和随机模拟的方法对解估计的不确定性程度进行评价，结果表明在空间变量结构约束下，解估计的

2、不确定性程度低于未受到约束的反演解，验证了此文方法的有效性。关键词地质统计学，AVO反演，贝叶斯理论，空间变量结构AVOBayesianinversionconstrainedbyspatialvariablestructureAbstractAnimportantsubjectingeophysicalinversionresearchistodecreasetheuncertaintydegreeofinversesolutionbyintroducingsubsurfaceprioriinformationfromdifferentsourcesasaconstraintofgeophy

3、sicalinversion.Geologicalvariableexhibitsspatial(time)variability.Apartfromtherandomnessintermsofstatistics,geologicalvariabledependsonspatial(time)locationaswell,possessingcertaincorrelationandstructuralproperty.Thepaperintroducesthetheoryofvariogramfromgeostaticstocreateadescriptionaboutthestructu

4、reofsubsurfacespatialvariableasprioriinformation.CombinedwithBayesianinferencetheory,weestablisharelationshipbetweenspatialvariableandgeophysicalobservationvialikelihoodfunction.Finally,wederivetheanalyticalformofAVOinversesolutionestimationsposteriorprobabilitydensity.Practicalloggingdatafromawelll

5、ocatedinanoilfieldinEastChinaSeaisselectedtoconductpilotinversecalculation.Afterwards,uncertaintydegreeofitssolutionestimationisevaluatedbyperformingintervalestimationandstochasticsimulationmethods.Theconsequencesuggeststhatwithspatialvariablestructureasconstraint,theuncertaintydegreeofsolutionestim

6、ationislowerthanitofinversesolutionwhichiswithoutconstraint.Thus,methodproposedinthispaperisprovedtobeeffective.KeywordsGeostatistics,AVOinversion,Bayesiantheory,Spatialvariablestructure0引言众所周知，地球物理反问题是一个多解且不适定的问题。对于地震反演来说，其根本目标是利用观测到的地震数据集，求解地下介质的结构与物理参数1。除了地震数据外，如何让引入更多不同来源的信息约束地震反演，降低反问题的多解性，获得更真

7、实描述地下结构与性质的解一直是地震反演的关键问题之一。由于反问题具有多解性，所以在反问题研究中，重要的不是构造出一种广义的解估计，而是要对各种可能的解估计进行评价，这正是由于反问题解的非唯一性所决定的。贝叶斯方法是一个解决上述问题很好的工具，Tarantola(1987,2005)使用关于模型空间的先验信息和似然函数构造了反问题解的后验概率密度，完整的描述了解估计的不确定性问题；2001年，Scales和Tenorio关于如何理解在反问题研究中的先验信息与不确定性关系，在其论文中给出了精彩的阐述肌Buland(2003,2006)等人应用贝叶斯原理在子波估计、AVO线性反演、时移地震反演方面给

8、出了研究方法和应用实例。2009年，陈建江等人专门论述了AVO反演的不确定性分析a】。在统计学概念下，将反演的地质变量认为是随机变量，一次观测(如测井数据)均为随机变量的一次实现。经典统计学中，通常假定某随机变量采集的样本是完全随机的或在空间上是完全独立的。但通常地，需要通过反演求得的地质变量在时空域上并不完全独立，而是呈现出一定的空间(时间)相关性。19世纪60年代，法国学者G.Matheron发展起来的地质统计学，考虑到空间变量的结构性，以变差函数、克里金估值、随机模拟为基本手段，是研究空间中某些既显示出随机性又显示出结构性的变量的数学地质方法11。1994年，Haas和Dubrule结合

9、地质统计学中的随机模拟思想以及地震数据的反演，提出了现在所谓的地质统计学反演12。随机反演(地质统计学反演)要求对地下模型空间产生多个实现，比较常用的方法是MonteCarlo等计算代价高昂的方法13，使得在实际应用中受到限制。本文在贝叶斯理论框架下，假设地下介质参数服从对数高斯分布，引入地质统计学变差函数的概念，以地质变量的空间分布结构作为先验对AVO反演进行约束。通过先验的地质变量空间分布和描述地球物理观测和地质变量关系的似然函数，建立了反演解的后验概率密度分布的解析形式。应用随机模拟和区间估计的方法，来定性和定量地评价反演解的不确定性程度。数值算例表明，相比于常规的贝叶斯反演方法，考虑地

10、质变量在空间上的结构性来约束反演，更符合实际情况，不确定性程度降低。1方法原理在一般贝叶斯理论框架下，地球物理反演问题的解以及解的不确定性，以后验概率密度的形式来描述：b(mId)L(dIm)p(m)(1)JL(dIm)p(m)dmp(m)为关于模型m分布的先验概率密度；L(dIm)为似然函数，是模型空间m向数据空间d的映射的概率表达；分母JL(dIm)p(m)dm在给定数据情况下为一个常数，也即是概率归一化因子；(mId)是所要求的后验概率密度分布。方便下文讨论，简写L(dIm)、1)式简写为：似然函数L(m)的建立，以(mId)为L(m)、b(m)，忽略掉常数JL(dIm)p(m)dm,b

11、(m)二L(m)p(m)(2)接下来，本节将分三个部分分别讨论模型先验分布p(m)、及后验概率密度分布的求取。1.1模型先验概率分布在各向同性介质中，地下介质的弹性性质可由纵波速度、横波速度、密度来完整地描述,同样地也可以由纵波阻抗、横波阻抗、密度来描述。在此本文选择关于地下介质的弹性参数:纵波速度Vp、横波速度Vs、密度P来作为反演目标参数。关于模型的先验分布在贝叶斯理论中所起的角色，在理论上一直是概率统计学家争论的焦点14。尽管理论上关于模型先验信息的理解还存在争论，但在地球物理实际应用中，在已知地质背景规律或钻井、测井资料指导下反演，已经是大势所趋。关于贝叶斯AVO反演，本文以测井数据作

12、为地下模型的先验信息。但不同于常规的测井约束反演，本文不仅从测井数据里提取关于地下介质的低频信息来约束反演，还进一步挖掘了测井数据里包含的地下介质弹性参数的空间相关性(也即地下空间变量的结构)来进一步约束指导反演。地质统计学里的变差函数是描述空间变量结构的天然数学工具，可以证明在二阶平稳条件下，变差函数与协方差函数等价，在描述空间变量结构方面下文将不再区分协方差函数和变差函数。关于地下空间变量的分布，大量井资料分析表明，弹性介质三参数(纵波速度、横波速度、密度)服从对数高斯分布形式。由以上讨论，本文定义关于空间变量结构的先验分布为：TOC o 1-5 h z HYPERLINK l bookm

13、ark22 o Current Document m(t)二lnVp(t),lnVs(t),lnp(t)tNg,E)(3)mm其中t为地震双程旅行时，g为模型期望，E为模型协方差矩阵：mm HYPERLINK l bookmark24 o Current Document g(t)=Em(t)=卩(t),卩(t),卩(t)t(4)mVpVsE二Covm(t),m(t)二E(t,t)(5)mijmij在二阶平稳假设下(5)式可简化为:E二Co讥m(t),m(t)二E(t)(6)mijm其中t=|ttI，12(C(t)Vp,VpC(t)Vs,VpIC(t)VpC(t)Vp,VsC(t)Vs,VsC

14、(t),VsC(t)Vp,PC(t)Vs,C(t)丿,7)关于先验模型均值gm、协方差函数Em矩阵元的计算，是通过测井数据数据的统计分析获得。gm反映了地质变量的低频信息，Em则分反映了地质变量的空间结构信息。由此建立了空间变量的先验信息，可以用来约束反演。注意的是(7)式先验协方差函数的计算稍显复杂，但是它是描述空间变量信息的载体,是进行空间变量结构约束的关键。为了使本文的方法理论简洁明晰，在这里不讨论(7)式的具体计算，其详细的计算方式与讨论见反演算例部分。1.2似然函数在概率统计上的似然函数描述的是估计的模型响应与真实数据的相似程度，由模型向数据的映射，也就是我们称之为的动力系统或正演方

15、程。在一定的假设下，可以建立地下弹性参数映射到地震数据的线性关系。1919年Zoeppritz建立了平面波入射时反射系数、透射系数的理论表达式15，但其数学表达式复杂、物理意义不直观，本文采用1980年由Aki和Richard提出来的弱纵波反射系数近似式16，近似式表达形式与Stolt与Welein(1985)所采用的一样：dInVp(t)dInVs(t)dInp(t)TOC o 1-5 h zR(t,0)=a(t,0)+a(t,0)+a(t,0)(8)vpdtVsdtPdt其中：1a(t,0)=(1+tan20)(9)Vp2a(t,0)=-4Vs(t)sin20(10)VsVp(t)2a(t

16、,0)-1(1-4Vs(t尸sin20)(11)p2Vp(t)2这里Vp(t)、Vs(t)可以认为是弹性参数的低频背景值，可以通过模型先验均值A获得，则m系数项a、a、a可以先验求得。(8)式对应的矩阵表达式为：VpVsPR-Am(12)系数矩阵A-a(t,0),a(t,0),a(t,0)，m是m关于时间的一阶导数。VpVsP则地震记录d表示为子波S与反射系数R的褶积加上噪音项e:d-SR+e-SAm+e(13)假设噪音服从零均值高斯分布：eN(0,E)(14)e将褶积运算SA记作算子G，由褶积元算的线性性质可知G为线性算子。由此建立了线性正演方程，将(13)式写为：d-SAm+e-Gm+e(

17、15)1.3后验概率分布由模型先验概率分布与似然函数，在线性高斯条件下可得后验概率分布为高斯分布4：TOC o 1-5 h zmIdN(p,E)(16)m|dm|d其中：g=p+EGt(GEGt+E)-i(d-Gp)(17)m|dmmmemE=E-EGt(GEGt+E)-iGE(18)m|dmmmem6力2其中，二学E(t)，=-E(t)。m6Tmm6T2m式(17)完整地描述了反演解的概率性质，包括了反演解的期望与不确定性。由于后验概率分布是高斯分布，反问题的极大后验概率解也就是后验概率的条件期望MAPm=g。这里取得的反演解是真实弹性参数的对数ln(Vp)、ln(Vs)、ln(Den)，m

18、id必须经过指数变换才可得到Vp、Vs、Den。2数值反演算例为了验证本文方法的可行性，选择东海某油田一口包含纵波速度、横波速度、密度测井数据的勘探井进行试验。测井数据经过标定，深度域变换到时间域(双程旅行时)，井数据选取时窗为3500ms3700ms，采样率1ms。井曲线及其选用的低频背景见图1。图1纵波速度、横波速度、密度(红色实线)；模型低频背景(黑色虚线)Fig.1P-wavevelocity,S-wavevelocity,density(redsolidline);Modellow-frequencybackground(blackdotline)由井曲线通过Aki-Richard纵

19、波近似反射系数公式计算得到的反射系数与主频50Hz雷克子波褶积获得角道集地震记录，图2所示。记录时长与井数据相同3500ms3700ms,采样率1ms。计算得到的合成记录作为反演时地震数据的输入。3500352535503575600SE姑36503675370037253750合成地震记录5101520253035Angle图2部分角度合成地震记录Fig.2Partialanglesyntheticseismogram对测井数据Vp、Vs、Den取对数，采用正态概率纸来检验对数弹性参数是否符合高斯分布。由图3可得，认为认为测井数据的对数弹性参数的高斯分布假设是可以接受的。关于模型先验信息，模

20、型均值取图1所示低频背景的对数，二阶平稳假设下认为模型参数在每个采样点上的对数方差为常数b2（t）=0.0162，b2（t）二0.0117，b2（t）=0.000403。关于模VpVsp型参数的空间结构，也即先验协方差函数计算见下一部分。99:99979-H-5999o.n.IJ._u_uoonun-9075502597985o999999o.go.gooDuo552o.o0.100.05.020.010.0030.001ln(Vp)(m/s)0.100.05.20.010.003.1ln(Vs)(m/s)0.003.0010.100.05.20.017.87.847.88ln(Den)(Kg

21、/ms)97985nu999999.93_UO._uon-5_u760_u52o图3对数弹性参数在正态概率纸上的分布Fig.3Distributionoflogarithmelasticparametersonnormalprobabilitypaper2.1先验空间变量结构计算地质统计学里，关于空间变量的结构是利用变差函数工具来进行描述的，变差函数反映了变量自身在空间上的结构性以及变量之间在空间上的依赖关系I。在本文中的空间变量就是需要反演的弹性参数，纵波速度Vp、横波速度Vs、密度Den。一般地，在应用地震数据反演前，我们已经通过钻井获得了关于井眼处地下介质的一些信息（如岩性测井、物性测井

22、、电性测井）。通过井数据的统计分析，可以计算地下介质弹性参数的变差函数。将井上计算到的变差函数引入到反演中来进行空间变量结构的约束，更符合地下地质客观实际。由于变差函数需要满足正定性条件，通常的方式是运用一些已知的、满足正定条件的变差函数的线性组合来建立变差函数模型。一些被证明有效的、常用的、基本的变差函数包括有球形模型、指数模型、高斯模型、DeWijsian模型等，具体的表达式及应用见参考文献11。图4是三种常用的变差函数模型的形态曲线。图4三种常用变差函数模型的形态曲线Fig.4Curvesofthreecommonvariogramfunctionmodel型型型状数斯球指高可以证明，在

23、二阶平稳下，变差函数与协方差函数之和为一常数，所以关于(7)式协方差矩阵元的计算可以通过变差函数的概念来求取。在二阶平稳假设下，变差函数与协方差函数在反应空间结构性方面具有等价的意义，在本文接下来，为了与上述反演理论一致，将使用协方差函数这个术语。在井上计算ln(Vp)、ln(Vs)、ln(Den)的自协方差函数以及相互之间的互协方差函数，通过与理论协方差函数的套合，本文选择高斯模型作为理论协方差函数，协方差函数的高斯模型的理论表达式见式(19),套合参数与效果见表1与图5。C(h)=Ce-(h/a)2(19)其中C为空间变量方差，即分离距离为0时的协方差，a为变程。值得注意的是，在选择理论模

24、型以及求取理论模型参数时，最小二乘法的曲线拟合一般来说不是一种最佳的，一般来说甚至不是合适的选择变差函数模型的方法，必须通过人为的观察和套合来产生理论模型。在套合过程中，应将更多的注意力集中在较小分离距离的协方差函数值。在这里我们选择前15个分离距离来套和协方差函数。U.U20.0150.010.0060C12468101214ln(Vp)自协方差函数图5协方差函数计算套和Fig.5Nestedanalysisofcovariancefunctionln(Vp)自协方差ln(Vs)自协方差ln(Den)自协方差ln(Vp).ln(Vs)互协方差ln(Vp).ln(Den)互协方差ln(Vs).

25、ln(Den)互协方差方差C0.01620.01174.01e-40.01120.00170.0016变程a8.205811.00716.295711.12875.94267.0830(ms)表1套和协方差函数参数Table1.Covariancefunctionparametersofnestedanalysis2.2反演结果不确定性分析这里主要讨论在空间变量结构约束下与没有空间结构约束情况下的反演结果比较。在这里选用极大化后验概率解作为反演解，也即后验验概率的数学期望。图7是空间结构约束下的反演，图6是不受空间结构约束的常规反演。从图上可知，对于两种反演结果的精度来说,两种反演方法都取得了

26、很好的解估计。受空间结构约束与否，似乎并不影响反演的准确度。正如前文的讨论可知，衡量一个解估计的好坏，除了反演精度之外，另一重要的方面是反演结果的不确定性程度。对于不确定性程度评估，此文应用了两种方法，一种方法是是对后验解空间进行随机模拟，定性地探索后验解估计空间的不确定性程度；另一种方法是，应用统计推断上置信区间的概念来进行分析，这是一种定量化的不确定性分析方法。juoO图6不受空间结构约束Vp、Vs、Den反演结果（黑色点线为先验低频背景值，红色阶梯线为实际井数据，蓝色线为反演结果）Fig.6InverseresultsofVp,Vs,Denwithoutspatialstructural

27、constraint(Blackdotlinereferstopriorilow-frequencygroundvalue,redstep-linereferstopracticalwelldataandbluelinereferstoinverseresult)图7空间结构约束Vp、Vs、Den反演结果黑色点线为先验低频背景值，红色阶梯线为实际井数据，蓝色线为反演结果）Fig.7InverseresultsofVp,Vs,Denwithspatialstructuralconstraint(Blackdotlinereferstopriorilow-frequencygroundvalue,

28、redstep-linereferstopracticalwelldataandbluelinereferstoinverseresult)2.2.1随机模拟不确定性分析事实上，由（16）式得到的是关于解空间的概率表达，它给出了每一个可能的解估计的概率。由（17）式和（18）式计算的期望和方差，唯一地决定了关于解空间后验高斯分布的所有性质。由于具有后验分布的解析表达，所以可以应用随机模拟这种方法来部分地探索解空间，给出定性的关于解空间的不确定性评价。在这里后验解空间进行50次随机模拟实现，实现结果如图8和图9。相比于在空间结构约束下的随机模拟实现，不受约束的解波动性更强，变化范围较大，说明解的

29、不确定性程度更强。而受约束的解，则很好地限制在离真实井数据较近的范围内波动。由这50次的随机模拟结果，可以定性地说明此文的方法确实一定程度上减小了解的不确定性程度。35003525355035753600362536503675370037253750Den24002500260027002800Eg/mJ图8不受空间结构约束解的随机模拟实现灰色线为随机模拟解，红色阶梯线为实际井数据）Fig.8Stochasticsimulationrealizationwithoutspatialstructuralconstraint(Greylinereferstostochasticsimulatio

30、nsolution;redstep-linerepresentspracticalwelldata)juFl-3VPjo535253525350535U53飞5.375530500256663336/3ju375237050025666333jo37523720002500300035004000nVsjo51324002500260027002800Kg/mJ5o5553v5Eh3jlj6352635063jo37635237图9空间结构约束下解的随机模拟实现灰色线为随机模拟解，红色阶梯线为实际井数据）Fig.8Stochasticsimulationrealizationconstrain

31、edbyspatialstructure(Greylinereferstostochasticsimulationsolution;redstep-linerepresentspracticalwelldata)2.2.2统计推断不确定性分析统计推断上研究解估计的可靠性程度，主要内容包括点估计和区间估计。本文采用区间估计里置信区间的概念来说明解的不确定性程度。具体的表述方式是，在给定同样解估计、同等置信水平下，应用解估计所处的置信区间的大小来表达表不确定性程度,置信区间越小，表示解的波动性越小，解的不确定性程度越小。对于高斯分布来说，在给定置信水平为1-Q的情况下，估计量的置信区间为：X土zb

32、(20)0/2其中X为估计量，b为估计量标准差，z为0分为点。这里取0.05，置信水平01-0二0.95，则z=1.96。（21）式变为X士1.96b。将式（18）后验期望作为解估计02量，式（19）可以计算解估计在每个采样点上的后验方差，所以可以得到Vp、Vs、Den的后验置信区间。图10与图11给出了不受空间结构约束和在空间结构约束下的先验置信区间和后验置信区间。图中明显看到,在空间结构约束下的置信区间相比于没受约束的要减小很多,从这个意义上讲认为受空间结构约束的获得的反演结果是一个更好的解估计，因为解估计被很好地限制了在两条窄的绿线之间，波动性更小，可靠程度更高。一个量化的解估计评价见表

33、2,列出了两种解估计对Vp、Vs、Den后验置信区间相对于先验置信区间的平均减小程度。350035035753600367636503675370040006000m/s35003526355035763600362636503G7637003253750Vs200030004000m/sDen35003625.35503600362536503675;70037-25-37502400-.2500260Q-2700图10不受空间结构约束的解估计先验、后验置信区间（黑色点线为先验置信区间，红色阶梯线为实际井数据，绿色线为后验置信区间）Fig.10Prioriandposteriorconfid

34、enceintervalofsolutionestimationwithoutspatialstructuralconstraint(Blackdotlinereferstoprioriconfidenceinterval,redstep-linereferstopracticalwelldataandgreenlinereferstoposteriorconfidenceinterval)ju53525350355.3050025666333(SUI)auIIHju37_L523735ju35%.50L!zl3ju63o563ju372537200030004U00nVsjo53525350

35、53下-5.3JLI_til352635063ju376237图11空间结构约束下解估计先验、后验置信区间(黑色点线为先验置信区间，红色阶梯线为实际井数据，绿色线为后验置信区间)Fig.11(BlackPrioriandposteriorconfidenceintervalofsolutionestimationwithspatialstructuralconstraintdotlinereferstoprioriconfidenceinterval,redstep-linereferstopracticalwelldataandgreenlinereferstoposteriorconfid

36、enceinterval)VpVsDen未约束28.5126.6517.19空间结构约束66.8163.0135.47表2解估计后验置信区间平均减小百分比（%）Table2.Averagedecreasepercentofsolutionestimationposteriorconfidenceinterval3结论文章假设地下介质弹性参数的对数符合高斯分布情况，在贝叶斯理论框架下，建立先验信息和似然函数的概率表达，得到了反演解估计的理论表达式，该理论表达式由一个后验概率密度函数的形式给出，完整地描述了反演解的位置信息(期望)和不确定性程度(方差)。在建立先验信息时，引入地质统计学里变差函数的

37、概念，以协方差函数的形式建立反演参数的空间结构信息，以此约束地震反演。通过某实际工区的单井反演算例，以置信区间的概念和随机模拟的方法评估解的不确定程度，验证了在空间结构约束下，相比于未约束的反演解估计不确定性程度要低。进一步研究的建议：本文只是通过单井的数据建立了垂向上的空间变量结构信息，如何引入横向上甚至全三维方向上空间变量结构信息约束地震反演，是下一步要研究且富有意义的课题。参考文献(References)1Cooke,D.A.,andW.A.Schneider,Generalizedlinearinversionofreflectionseismicdata:Geophysics,198

38、3,48(6),665676,杨文采.地球物理反演的理论与方法.北京:地质出版社,1997,1-13.YangWencai,Theoryandmethodofgeophysicalinversion,Beijing:ChinaLandPress,1997,1-13Tarantola,A.,Inverseproblemtheory:Methodsfordatafittingandmodelparameterestimation,1987:Elsevier,ScientificPubl.Co.Inc.Tarantola,A.,Inverseproblemtheoryandmethodsformod

39、elparame-terestimation,2005:SocietyforIndustrialandAppliedMathematics.JohnA.ScalesandLuis,Priorinformationanduncertaintyininverseproblems.Geophysics,2001,66(2):389-397Buland,A.,andH.Omre,BayesianlinearizedAVOinversion:Geophysics,2003a,68(1),185-198.Buland,A.,andH.Omre,Bayesianwaveletestimationfromseismicandwe