高斯束层析偏移速度建模方法及应用

邵荣峰,方伍宝,蔡杰雄,倪　瑶,李　辉,郭立鹏

(1.中国石油化工股份有限公司石油物探技术研究院,江苏南京211103;2.中国石油大学(华东)地球科学与技术学院,山东青岛266580;3.同济大学海洋与地球科学学院波现象与反演成像研究组,上海200092)



高斯束层析偏移速度建模方法及应用

邵荣峰1,2,方伍宝1,2,蔡杰雄1,倪瑶1,李辉3,郭立鹏2

(1.中国石油化工股份有限公司石油物探技术研究院,江苏南京211103;2.中国石油大学(华东)地球科学与技术学院,山东青岛266580;3.同济大学海洋与地球科学学院波现象与反演成像研究组,上海200092)

摘要:常规地震层析反演分为射线层析和波动方程层析两类，基于射线类的常规层析方法由于理论假设的限制难以对复杂地区进行准确的速度建模;波动方程层析需要建立高精度的初始模型,其应用受到实际资料品质的限制。为此，开展了高斯束层析偏移速度建模方法研究,利用高斯束核函数构建病态性更小的灵敏度矩阵以提高层析反演的稳定性;采用高斯束叠前深度偏移进行高陡构造成像,并直接提取高分辨率的角度域共成像点道集(ADCIGs)用于高斯束层析反演,进而获得精确的偏移速度模型。为了避免在高陡构造地区人工拾取的繁杂,采用自动拾取技术获取反射点坐标与地层局部构造倾角。模型测试和实际资料应用结果显示,该方法具有较高的反演精度和稳定性。

关键词：高陡构造;高斯束层析;角度域共成像点道集;核函数;自动拾取

20世纪70年代,BOIS等[1]将层析反演从医学领域引入地球物理学,开启了地震层析反演研究的序幕。根据理论基础不同,通常将地震层析反演分为射线层析和波动方程层析两类。射线层析由于理论简单并且计算效率高,已经在生产实践中得到广泛的应用。但是传统的射线层析基于费马原理,认为波沿着介质中的一条射线传播,利用传播的射线路径和地震波走时来反演地下介质速度；同时传统的射线仅能覆盖有限的区域,构建的层析矩阵非常稀疏[2],存在较大的零空间,使得反演方程的求解存在极大的不适定性[3],影响成像分辨率。这些问题严重制约着传统射线层析方法的发展。

为了弥补射线层析方法存在的不足,DEVANVEY[4],TARANTOLA[5]等于20世纪80年代提出了波动方程层析成像方法。波动方程层析[6]基于Born近似和Rytov近似,从理论波场与观测波场的最佳匹配入手,直接反演速度扰动,理论上具有更高的成像分辨率,但它需要计算正向传播波场与反向传播的剩余波场以及两者的相关,非常耗时,同时反演的目标泛函与速度扰动之间存在强烈的非线性关系,需要建立一个高精度的初始模型,而这恰恰是实际生产中很难实现的目标。这些问题如今阻碍着波动方程层析的生产应用。

针对以上问题,胖射线层析[7]、菲涅尔体层析[8]、高斯束层析[9-10]等一些折中的方法被相继提出。其中POPOV等[9]、SEMTCHENOK等[10]提出的高斯束层析方法兼容了传统射线层析计算速度快和相对稳定的优点,通过对一个炮检对建立多个层析方程来减少层析方程的病态性,在一定程度上提高了反演精度和分辨率,但是该方法建立灵敏度矩阵的基础仍然是射线,反演过程中还需要加入正则化。

本文基于高斯束核函数建立的层析反演方程进行偏移速度分析[11]。首先通过高斯束叠前深度偏移[12-13]得到偏移剖面和角度域共成像点道集(ADCIGs);然后根据全自动拾取技术拾取反射点和局部地层倾角,计算出剩余曲率并转化为旅行时残差,利用高斯束核函数计算得到的灵敏度矩阵与旅行时残差建立层析反演方程组,求解并经过多次迭代最终得到满足要求的速度模型。

1方法原理

1.1高斯束层析理论

近年来,POPOV等[9]、SEMTCHENOK等[10]利用高斯束进行旅行时层析成像取得了很好的效果。其层析方程如下:

(1)

式中:si是网格内第i个高斯束中心射线长度;c是网格内慢度差;x0是地下反射点;xg是接收点;U是由震源点激发、接收点xg观测的波场;Uk是第k个高斯束射线计算的观测波场;τg是接收点的旅行时差;φk代表第k个高斯束射线。

根据一个炮检对存在一个旅行时差,传统射线层析只能建立一个层析方程。而从(1)式可以看出,位于接收点附近的射线都可以建立层析方程,这里的“附近”是指接收点位于射线的高斯束范围之内。利用高斯束对波场的贡献大小建立多个层析方程(方程的个数k取决于射线的密度,即k等于2倍高斯束宽度内射线的个数),减少了层析反演的稀疏性,从而降低了层析反演的病态性。但是,这里的高斯束层析灵敏度矩阵仍然建立在射线的基础上,求解过程中往往不收敛。为了克服这一缺陷,本文利用高斯束体建立核函数，构建灵敏度矩阵,增加了层析反演的稳定性。

中心射线附近没有焦散现象时,高斯束层析方程可以写成空间积分[14]的形式:

(2)

其中，

(3)

式中:K(x)是高斯束核函数;p(x)是慢度更新量;τ是走时残差;A[xg,xs,φ]是高斯束振幅,是空间点x的函数;∫ΨA[xg,xs,φ]dφ是高斯束振幅的积分,是与x无关的常数,Ψ是加权系数;xs是震源点;φ是高斯束中心射线标记。从(3)式可以看出,高斯束核函数不仅与连接观测点和接收点的射线(中心射线)路径有关,还与中心射线附近一定宽度范围的高斯束有关。所以高斯束层析更加符合地震波传播的物理规律,改善了传统射线路径对高、低速区域灵敏度的差异所产生的问题。

(4)

式中:φ0是连接炮检点的射线标记;v(φ)是观测点xg在射线φ上垂足x0处的投影速度；Q(φ)是x0处的动力学射线追踪参数；q0是观测点xg到x0的距离；M(φ0)是x0处旅行时二阶导数,ω是圆频率。

在中心射线附近没有焦散现象的前提下,对射线的积分可以转换为:

(5)

将(4)式和(5)式代入(3)式,可得到高斯束层析核函数:

(6)

其中,

(7)

式中:q是观测点xg到射线的垂直距离;qmax是观测点到射线的最大垂直距离;r是中心射线上的点到假想射线的距离;rmax是中心射线上的点到最远射线的距离。根据(6)式可以计算出(3)式所示的高斯束核函数。图1是传统射线路径与高斯束核函数覆盖范围对比图,其中背景速度场是常梯度速度模型,黑色虚线为传统射线路径，湖兰色区域为高斯束核函数覆盖范围。

图1　传统射线路径与高斯束核函数覆盖范围对比

由图1可见,传统的射线路径只是一条射线,通过层析更新速度模型时,只更新射线到达的区域,覆盖范围小,增加了灵敏度矩阵的稀疏性,因而得不到准确的反演结果。高斯束核函数不是单一射线,而是可以覆盖更多区域的波束体,降低了高斯束层析反演矩阵的稀疏性,因而能得到更加精确的速度模型。

1.2反射点、地层倾角自动拾取以及剩余曲率自动求取

在高斯束层析反演过程中,每迭代一次都需要进行偏移并拾取偏移剖面上的反射点坐标,再利用角道集(ADCIGs)中剩余曲率与反射点处局部地层角度的关系式拟合求取剩余曲率。为了快速拾取成像剖面上的反射点和求取成像道集上的剩余曲率,本文引入一种自动拾取技术。

根据反射点落在波峰(或波谷)上的基本特性,我们将断层检测方法[20-22]应用到反射点的自动拾取中。图2给出了波峰(波谷)的识别准则，即：

(8)

图2　波峰、波谷识别准则

根据反射点位于同相轴上,即反射点处的图像在横向邻域范围内具有局部线性性,可以进一步筛选潜在反射点。本文采用图像处理中的结构张量算法[23-24]来计算反射点处的局部线性性与局部地层倾角。设I为二维地震图像,I中表示空间方向信息的结构张量由图像梯度值定义。结构张量表示区域的变化方向和沿变化方向的变化量大小,地震地层纹理和断层纹理则由局部各点的方位信息变化关系确定。引入高斯函数模糊局部细节,可以使得结构张量突出显示区域内信号的复杂性。对于二维图像,结构张量是一个2×2的矩阵:

(9)

其中,gx与gy代表地震图像沿水平方向和垂直方向的梯度,<·>代表二维高斯光滑滤波。对于半正定矩阵G,特征值与特征向量可由|G-λI|=0得到,即:

(10)

式中:λ1是最大特征值,表示结构张量在第一个特征向量方向V1的能量;λ2是最小特征值,表示结构张量在第二个特征向量方向V2的能量。(λ1-λ2)/λ1是线性度,反映局部方向的一致性。特征向量描述了图像局部线性结构的方向性,针对图像的每个点,特征向量V1正交于图像的主结构方向,特征向量V2平行于图像的主结构方向。根据结构张量算法的物理意义,可以计算得到图像中任一点的局部线性性指标(λ1-λ2)/λ1与局部图像切向(即地层局部倾角方向)方向单位向量。值得注意的是,结构张量算法能适应低信噪比地震数据,因此可用于稳健地拾取地下反射点与局部地层倾角。

自动拾取反射点和地层倾角后,利用角道集(ADCIGs)中剩余曲率与反射点坐标位置的关系式拟合求取Δz,计算公式如下[25]:

(11)

其中,

(12)

式中:z0是自动拾取的反射点深度;γ是拾取的反射点深度与真实深度的比值;β是角道集的入射角度;M是与角道集对应的控制点数;zai是M个控制点的偏移深度;Δz是剩余曲率。

图3是公式(12)在z0=1km时不同γ值的曲线分布,可以看出,当偏移速度正确时,成像道集被拉平,此时γ=1,偏移深度与真实深度吻合;当偏移速度不正确时,γ≠1,不同角度的偏移深度与角度之间存在一定的曲率关系。当偏移速度偏小时,γ<1,成像深度曲线成上翘趋势;当偏移速度偏大时,γ>1,成像深度曲线成下弯趋势。通过这种曲率关系在ADCIGs中进行一系列曲率值扫描就可以得到剩余曲率Δz,最后根据地下成像点局部位置关系求取旅行时残差,计算公式[25]如下:

Δt=2sΔzcosαcosβ

(13)

式中:Δt是旅行时残差;s是地下成像点局部慢度值;α是反射层倾角。

自动求取剩余曲率的步骤是:①通过自动拾取的反射点深度z0,扫描不同的(z0,γ)参数对,对每一组参数对应用公式(12)所示的共成像点道集的曲线关系式进行叠加,将同相叠加能量最大时对应的参数对(z0,γ)视为最佳拟合的实际共成像点道集参数;②将最佳拟合得到的γ值代入(11)式求得该共成像点道集每个反射角度对应的剩余曲率Δz,利用(13)式将剩余曲率Δz转化为旅行时残差,代入(2)式可求取慢度更新量。

图3　z0=1km时不同γ值的曲线分布

1.3高斯束层析偏移速度分析流程

基于以上分析可知,高斯束层析偏移速度建模方法的核心在于高斯束层析核函数推导以及全自动拾取技术。具体实现步骤如下:

1) 根据初始的速度场进行高斯束叠前深度偏移,获取偏移剖面和相应的角度域共成像点道集ADCIGs。由于ADCIGs是在高斯束叠前深度偏移过程中提取的,因此能够准确地反映速度与深度的耦合关系,使得反演结果具有较高的精度。

2) 利用全自动拾取技术在偏移剖面上拾取层位界面以及局部倾角,根据(11)式求取深度残差,并转化为旅行时残差。

3) 利用高斯束层析核函数((6)式)构建灵敏度矩阵,将灵敏度矩阵和旅行时残差代入(2)式,利用联合迭代重建(SIRT)方法[26]求解反演方程组(需加入正则化进行约束)得到慢度更新量,进而更新速度模型。

4) 用更新得到的速度模型再次进行高斯束叠前深度偏移,并提取角道集ADCIGs,根据ADCIGs拉平程度以及速度精度要求判断是否进行下一次迭代,是则重复上述步骤,直到速度模型满足精度要求为止。

2模型试算

采用复杂理论模型(图4a)对高斯束层析反演方法的有效性进行了测试。该模型断层发育,地层界面起伏变化较大,断层处速度横向变化剧烈。横纵向采样点数为731×550,纵向采样间隔5m,横向采样间隔10m。首先正演叠前地震记录(图4b),然后将常梯度速度模型作为层析初始速度模型(图5a),进行高斯束叠前深度偏移(图5b)。可以看出,复杂模型中的各个层位均已经归位,但是绕射波没有收敛,且层位界面也不在正确的深度位置。图6a是对图5b进行自动拾取的结果(绿色点线),图6b,图6c,图6d是CDP200,CDP400,CDP600处角道集自动拾取结果。由图6a 可以看出自动拾取的反射点与反射层位吻合,断层、陡倾角等复杂区域也得到了很好的拾取。图7是在CDP200处提取的初始角道集和层析偏移后角道集对比结果,初始角道集同相轴存在上翘现象(图7a),说明初始速度偏小;而经过层析反演得到的速度模型偏移后,角道集被拉平(图7b,图7c),且高斯束层析方法比常规射线层析方法提取的角道集更加接近真实速度模型提取的角道集。图8 是经过5次高斯束层析迭代得到的速度模型以及相应的高斯束叠前深度偏移剖面,可以看到绕射波收敛,反射界面被归位到正确的深度位置,地堑处的同相轴聚焦更好,说明高斯束层析反演在横向速度变化剧烈的区域优于常规射线层析方法。图9是经过8次常规射线层析得到的速度模型以及相应的叠前偏移剖面，将图8与图9对比可以发现,高斯束层析能够提供更加丰富的地质信息,迭代次数也远小于常规射线层析。图10为真实深度、高斯束层析偏移深度以及常规射线层析偏移深度对比,可以看出,高斯束层析偏移结果更加逼近真实速度模型。

图4　复杂理论模型(a)及其正演地震记录(b)

图5　常梯度初始速度模型(a)及高斯束叠前深度偏移剖面(b)

图6　图5b反射点自动拾取结果(a)以及CDP200(b),400(c),600(d)的角道集自动拾取结果

图7　CDP200处角度域共成像点道集a 初始角道集; b 正确速度模型偏移后角道集; c 高斯束层析偏移后角道集; d 常规射线层析偏移后角道集

图8　高斯束层析更新的速度场(a)及高斯束叠前深度偏移剖面(b)

图9　常规射线层析更新的速度场(a)及高斯束叠前深度偏移剖面(b)

图10　真实深度、高斯束层析偏移深度以及常规射线层析偏移深度对比

3实例应用

中国南方和西部含油气区地下发育复杂高陡构造，速度横向变化剧烈,构造成像及速度建模均面临着巨大的挑战。利用某探区一条二维测线对高斯束层析偏移速度建模方法进行了测试。该测线共264炮,每炮240道,道间距为40m,偏移距范围20～4780m。炮记录采样点数1500,采样间隔4ms。CDP范围1～1278,CDP间隔20m,最大深度8.0km。ADCIGs范围0～45°,角度间隔1°。图11为该探区一炮集记录,图12是初始速度模型以及对应的高斯束叠前深度偏移剖面,可见构造形态的成像模糊不清。图13是对图12b中的反射界面进行自动拾取的结果(绿色点和线),可以看到拾取的反射点与偏移剖面上的层位相对应。图14是角道集自动拾取结果,图15为经过10次高斯束层析更新后的速度场以及对应的高斯束叠前深度偏移剖面,高陡构造形态清晰可见。图16是高斯束层析更新速度前后在CDP900处提取的角道集。更新前角道集同相轴存在下弯现象,说明初始速度偏大(图16a);更新后速度逼近真实速度,所以角道集同相轴被拉平(图16b)。

图11　某探区实际炮集记录

图12　某探区初始速度模型(a)及高斯束叠前深度偏移剖面(b)

图13　图12b中反射界面自动拾取的反射点

图14　角道集自动拾取结果a CDP900处; b CDP901处

图15　层析更新后的速度模型(a)及对应的高斯束叠前深度偏移剖面(b)

图16　速度更新前(a)后(b)CDP900处提取的角道集

4结论

1) 基于高斯束核函数建立的层析反演方程能大大减弱层析矩阵的病态性,提高反演的稳定性。

2) 高斯束叠前深度偏移不仅有利于高陡构造成像,而且能够直接提取高分辨率的角度域共成像点道集用于高斯束层析反演,获得精确的偏移速度模型。

3) 自动拾取反射点和局部地层倾角的方法不仅准确度高,而且避免了人工拾取的繁琐,提高了层析速度分析的效率。

致谢:本研究得到了中国石油化工股份有限公司石油物探技术研究院地震成像技术研究所的软硬件支持以及中国石油大学(华东)地震波传播与成像课题组李振春、杨国权和曹文俊老师的理论指导,在此表示衷心的感谢。

参考文献

[1]BOIS P,LA PORTE M,LAVERGNE M.Well-to-well seismic measurements[J].Geophysics,1972,37(3):471-480

[2]刘玉柱,董良国.初至波层析影响因素分析[J].石油地球物理勘探,2007,42(5):544-553

LIU Y Z,DONG L G.The analysis of first wave tomography[J].Oil Geophysical Prospecting,2007,42(5):544-553

[3]王华忠,冯波,王雄文,等.地震波反演成像方法与技术核心问题分析[J].石油物探,2015,54(2):115-125

WANG H Z,FENG B,WANG X W,et al.Analysis of seismic inversion imaging and its technical core issues[J].Geophysical Prospecting for Petroleum,2015,54(2):115-125

[4]DEVANVEY A J.Geophysical diffraction tomography[J].IEEE Transactions on Geoscience and Remote Sensing,1984,22(1):3-13

[5]TARANTOLA A.Inverse problem theory,methods for data fitting and model parameter estimation[M].New York:Elsevier,1987:350-351

[6]WOODWARD M J.Wave-equation tomography[J].Geophysics,1992,57(1):15-26

[7]MICHELENA R J,HARRIS J M.Tomographic traveltime inversion using natural pixels[J].Geophysics,1991,56(5):635-644

[8]VASCO D W,PETERSON J E,MAJER E L.Beyond ray tomography:wavepaths and Fresnel volumes[J].Geophysics,1995,60(6):1790-1804

[9]POPOV M M,SEMTCHENOK N M,VERDEL A R,et al.Reverse time migration with Gaussian beams and velocity analysis applications[J].Expanded Abstracts of 70thEAGE Conference & Exhibition,2008:F048

[10]SEMTCHENOK N M,POPOV M M,VERDEL A R.Gaussian beam tomography [J].Expanded Abstracts of 71stEAGE Conference & Exhibition,2009:U032

[11]WOODWARD M J,NICHOLS D,ZDRAVEVA O,et al.A decade of tomography[J].Geophysics,2008,73(5):VE5-VE11

[12]HILL R N.Gaussian beam migration[J].Geophysics,1990,55(11):1416-1428

[13]蔡杰雄,方伍宝,杨勤勇.高斯束深度偏移的实现与应用研究[J].石油物探,2012,51(5):465-475

CAI J X,FANG W B,YANG Q Y.Realization and application of Gaussian beam depth migration[J].Geophysical Prospecting for Petroleum,2012,51(5):465-475

[14]LI H,WANG H Z.An effective tomography algo-rithm with Gaussian beam[J].CPS/SEG International Ge-ophysical Conference,2014:701-704

[20]TANER M T,KOEHLER F.Velocity spectra-digital computer derivation and applications of velocity functions[J].Geophysics,1969,34(6):859-881

[21]LUO S,HALE D.Velocity analysis using weighted semblance[J].Geophysics,2012,77(2):U15-U22

[22]MARFURT K J,FARMER S L,BAHORICH M S.et al.3-D seismic attributes using a semblance-based coherency algorithm[J].Geophysics,1998,63(4):1150-1165

[23]VLIET L J V,VERBEEK P W.Estimators for orientation and anisotropy in digitized images[J].Proc First Conference of the Advanced School for Computing & Imaging,1995:442-450

[24]LUO Y,WANG Y E,ALBINHASSAN N M,et al.Computation of dips and azimuths with weighted structural tensor approach[J].Geophysics,2006,71(5):V119-V121

[25]秦宁,李振春,杨晓东.高斯波束角道集的旅行时层析速度分析[J].石油地球物理勘探,2013,48(3):366-371,442

QING N,LI Z C,YANG X D.The travel time tomography based on Gaussian beam angle-domain common imaging gathers[J].Oil Geophysical Prospecting,2013,48(3):366-371,442

[26]和锐,杨建思,张翼.地震层析成像方法综述[J].理论与应用研究,2007,16(1):35-48

HE R,YANG J S,ZHANG Y.Review of seismic tomography methods[J].Theory and Applications,2007,16(1):35-48

(编辑：戴春秋)

A method of migration velocity analysis based on Gaussian Beam tomography and its application

SHAO Rongfeng1,2,FANG Wubao1,2,CAI Jiexiong1,NI Yao1,LI Hui3,GUO Lipeng2

(1.SinopecGeophysicalResearchInstitute,Nanjing211103,China;2.SchoolofGeosciences,ChinaUniversityofPetroleum,Qingdao266580,China;3.WavePhenomenaandInversionImage(WPI),SchoolofOceanandEarthScience,TongjiUniversity,Shanghai200092,China)

Abstract:Conventional seismic tomography includes ray tomography and wave-equation tomography.As a result of the limitation of theoretical assumptions,the conventional tomographic methods based on raypath are unable to build accurate velocity model in area with complicated geological conditions.Moreover,wave-equation tomography fails to get a good application result,because it needs to establish the initial model with high precision.The method of migration velocity analysis based on Gaussian beam tomography (GBT) depends on the alternately iterative of migration and tomography to achieve velocity update,by which we can obtain the relatively correct velocity model.This method can improve the stability of inversion by constructing less ill-posed Frechet matrix using Gaussian beams’ kernels.The pre-stack Gaussian beam depth migration not only has great advantages in high and steep structure zones,but also can extract angle-domain common imaging gathers (ADCIGs) with high resolution directly.By inputting the ADCIGs to Gaussian beam tomography,we will be able to obtain accurate velocity model.In order to avoid complex artificial pickup in the high-steep structure area,we proposed automatic picking technique to pick up reflection points and the local stratigraphic structural dip.The model test and field data prove that the Gaussian beam tomographic method has high inversion accuracy and stability.

Keywords:high and steep structure,Gaussian beam tomography,angle-domain common imaging gathers (ADCIGs),kernels,automatic picking

文章编号：1000-1441(2016)01-0091-09

DOI:10.3969/j.issn.1000-1441.2016.01.012

中图分类号：P631

文献标识码：A

基金项目：国家科技重大专项(2011ZX05014-001-002)项目资助。

作者简介：邵荣峰(1987—),男,硕士在读,主要从事地震成像与速度反演研究。

收稿日期：2015-05-22;改回日期：2015-10-10。

邵荣峰,方伍宝,蔡杰雄,等.高斯束层析偏移速度建模方法及应用[J].石油物探,2016,55(1):-99

SHAO Rongfeng,FANG Wubao,CAI Jiexiong,et al.A method of migration velocity analysis based on Gaussian Beam tomography and its application[J].Geophysical Prospecting for Petroleum,2016,55(1):-99

This research is financially supported by the National Science and Technology Major Project of China (Grant No.2011ZX05014-001-002).