冉 崎,陈 康,韩 嵩,马 博,刘 畅,曲英铭

(1.中国石油天然气股份有限公司西南油气田分公司勘探开发研究院,四川成都610021;2.中国石油大学(华东)地球科学与技术学院地球物理系,山东青岛266580)

山前带地区地表起伏剧烈,复杂的高陡构造发育。剧烈变化的起伏地表给常规地震勘探工作带来了极大的挑战。常规起伏地表地震资料处理方法包括静校正、基准面校正以及直接成像法[1-4],但当地表纵横向速度变化剧烈时,利用常规静校正技术和基准面校正方法难以得到准确的成像结果,因此,起伏地表叠前深度偏移方法得到了快速发展。

起伏地表深度偏移成像方法包括射线类、单程波动方程类和双程波动方程类偏移成像三大类方法。射线类偏移方法主要包括Kirchhoff偏移[5-6]和束偏移[7-13]方法。该类偏移方法运算效率高且对于起伏地表具有较好的适应性。WIGGINS[5]提出了面向起伏地表地震数据的Kirchhoff偏移方法,在此基础上,JAGER等[6]提出了保幅Kirchhoff起伏地表偏移方法。秦宁等[7]针对复杂山前带成像难题提出了一种高斯束叠前深度偏移技术。GRAY[8]将束偏移方法引入Kirchhoff偏移方法,既保留了高计算效率,又提升了偏移剖面的信噪比。岳玉波等[9]基于平面波分解,提出面向起伏地表地震数据的保幅高斯束成像方法。黄建平等[10]将岳玉波等[9]提出的面向起伏地表地震数据的保幅高斯束成像方法推广到了弹性介质。杨继东等[11]提出了一种适用于陆地复杂地表条件的叠前菲涅尔束偏移方法。徐秀刚等[12]基于“逐步累加”的“直接下延”法,将小束源与Fourier传播算子相结合,提出了基于起伏地表的小波束叠前深度偏移成像技术。黄建平等[13]基于有效邻域波场近似理论提出了一种成像精度更高且适用于复杂起伏地表条件的叠前保幅高斯束偏移方法。虽然射线理论凭借较高的运算效率,在工业界已经得到了广泛应用,但受方法本身的限制,射线类偏移方法只能反映运动学特征,不能反映动力学特征,且无法刻画高陡构造等复杂地质体。因此,基于单程波动方程的叠前深度偏移方法得到了快速发展。陈林谦[14]推导并分析了带误差补偿的频率空间域有限差分法单程波延拓方法,在此基础上,提出了带误差补偿的起伏地表叠前深度偏移技术流程。段心标等[15]提出了基于GPU的三维起伏地表裂步傅里叶叠前深度偏移技术。张宇等[16]提出了真振幅全倾角单程波方程偏移方法。近些年基于单程波方程的起伏地表直接成像方法发展较为缓慢,因为就成像精度而言,单程波成像方法不如逆时偏移方法;就计算效率而言,单程波成像方法不如射线类方法。随着计算机技术的快速发展,基于双程波方程的起伏地表叠前逆时偏移技术得到了发展。逆时偏移成像方法因无倾角限制,在高陡构造成像中发挥了重要作用。WHITMORE等[17]、BAYSAL等[18]和LOEWENTHAL等[19]均提出了成像效果良好的逆时偏移成像技术。刘红伟等[20]给出了针对起伏地表直接进行叠前逆时偏移的实现过程。王延光等[21]提出利用逆时偏移实现基于起伏地表的波动方程叠前深度偏移。鲁雁翔等[22]推导出了弹性波正向传播和逆时传播的曲线网格差分格式,实现了起伏层状介质中的波场逆时偏移。兰海强等[23]推导出了一种曲线坐标系下稳定的、显式二阶精度有限差分方法的弹性VTI方程。李庆洋等[24]提出了曲线坐标系下基于仿真型有限差分的贴体全交错网格波场延拓算法。

山前带地区高陡构造的成像难度极大,逆时偏移成像方法虽可以改善高陡构造的成像精度,但因观测系统的局限性,高陡构造的成像能量依旧较弱。传统的地震勘探中,通常将棱柱波视为噪声而将其压制[25]。事实上,棱柱波携带了许多一次波所不包含的地下高陡构造的信息。因此,棱柱波可以被用来提升照明度和高陡构造的成像效果[26-30]。CAVALCA等[31]阐明了棱柱波包含两个反射点和3段反射路径。MARMALYEVSKYY等[32]提出了一种基于克希霍夫的棱柱波成像方法来描述高陡盐丘侧翼。MALCOLM等[33]采用一种迭代单程波算法逐步对一次波和多次波成像。BROTO等[34]通过有限差分正演分析了棱柱波的运动学特征并提出了一种鉴别棱柱波的方法。LI等[35]提出了针对复杂盐丘构造的棱柱波成像流程。DAI[36]提出了一种无需提前拾取反射界面且对初始速度模型的精度要求较低的棱柱波逆时偏移成像方法。MAKSYM等[37]提出了最小二乘棱柱波逆时偏移方法;YANG等[38]提出了一次波与棱柱波联合最小二乘逆时偏移方法。

本文分别介绍了曲坐标系下的声波方程、棱柱波逆时偏移方程以及曲坐标系棱柱波逆时偏移方程,提出了适用于山前带高陡构造成像的起伏地表棱柱波逆时偏移方法,并利用两个典型模型分别对该方法进行了测试,证明了方法的正确性与有效性。

1 方法原理

1.1 曲坐标系下的声波方程

二维声波介质一阶应力-速度方程表示为:

(1)

式中:σ为应力;vx和vz分别为质点在x方向和z方向的速度分量;ρ为密度;vP为纵波传播速度。

本文采用蒋丽丽等[39]提出的正交贴体网格技术对速度模型进行剖分,通过调节网格线与边界的间距、网格线与边界的夹角,缩小了实际值与期望值的差距。基于泊松方程的贴体网格表示为:

(2)

式中:P(x,z),Q(x,z)为不同的“源项”,其作用是调节网格间距与网格形状;x,z为原始笛卡尔坐标系下不同网格节点的坐标值;ξ,η为转换后的曲坐标系下的不同网格节点的坐标值;下角标代表求二阶导数。

物理域(笛卡尔坐标系)与计算域(曲坐标系)关系表示为:

(3)

(3)式两边分别对x和z求导,表示为:

(4)

(5)

整理(4)式和(5)式得到的坐标变换关系表示为:

(6)

式中:J为雅可比行列式。

联立(1)式和(6)式并结合链式法则得到曲坐标系下的波动方程为:

(7)

1.2 棱柱波逆时偏移

棱柱波作为一种多次波,携带了更多的地质体信息,特别是来自地下高陡构造的信息。在描述高陡构造方面,棱柱波逆时偏移可以发挥很大的作用。根据棱柱波的产生机理可将其分为两种类型,第一种棱柱波是地震波先在低倾角反射界面上发生反射,随后又在高陡体侧翼发生反射,传播至地面被检波器接收;第二种棱柱波是地震波先在高陡体侧翼发生反射,随后又在低倾角沉积界面上发生反射,传播至地面被检波器接收。

棱柱波逆时偏移是针对高陡构造而发展出的一种多次波成像技术。该方法包括两个部分:①首先输入偏移速度场以及炮记录,应用常规逆时偏移成像得到低倾角反射层位置;②对第①步得到的成像结果继续应用棱柱波逆时偏移及其特定成像条件,得到棱柱波逆时偏移结果。

频率域共炮域逆时偏移公式可表示为[24]:

(8)

式中:mmig(x|xs)为逆时偏移波场值;W*(ω)为震源共轭频谱;xg为检波点位置;G*为共轭格林函数;d(xg|xs)为炮记录;xs为震源点位置。

对于给定的偏移速度场,可以求得震源波场的格林函数G0,进而可求得反射波场的格林函数G1,即[40]:

(9)

式中:m1(x′)为水平反射层在常规逆时偏移处理后得到的反射系数模型。

由此可得笛卡尔坐标系下的频率域棱柱波逆时偏移公式,可表示为:

(10)

式中:d1(xg|xs)为一阶散射反射波;d2(xg|xs)为二阶散射棱柱波;P1,Q0可由如下公式求得:

(11)

式中:ω为频率;v0(x)为偏移速度;P0(x)为震源波场;P1(x)为反射波场;W(ω)为震源频谱;Q0(x)为检波点波场;m1(x)为水平反射层经常规逆时偏移处理后得到的反射系数模型。

1.3 曲坐标系棱柱波逆时偏移

(7)式还可以写成如下形式,即:

(12)

速度扰动可表示为:

vP=vP0+δvP

(13)

式中:vP0为背景速度;δvP为扰动速度。

应力场p0可通过下式求解:

(14)

由泰勒公式可知,背景速度与扰动速度的关系可表示为:

(15)

将(15)式代入(12)式中,并减去(14)式,可得:

(16)

反射系数模型可表示为:

(17)

将(17)式代入(16)式,可得:

(18)

本文中,扰动场(δp,δvx,δvz)为反射波场(p1,vx1,vz1)。采用常规逆时偏移成像结果I代替m(vP),(18)式可表示为:

(19)

曲坐标系棱柱波逆时偏移的实现主要包括如下9个步骤:①输入笛卡尔坐标系下的地震炮记录和偏移速度模型;②计算曲网格并将偏移速度转换到曲坐标系下;③计算正向延拓的震源波场和逆时延拓的检波点波场;④得到曲坐标系下的常规逆时偏移;⑤在步骤③的正向延拓震源波场的基础上,再一次进行波场正传计算以正向延拓反射波场;⑥在步骤③的逆时延拓检波点波场的基础上,再一次进行波场逆传计算逆时延拓的反射波场;⑦得到曲坐标系下的棱柱波逆时偏移成像结果;⑧将棱柱波逆时偏移的成像结果反变换到笛卡尔坐标系下;⑨输出棱柱波成像结果。

2 模型试算

2.1 简单起伏地表凹陷模型

笛卡尔坐标系下的简单起伏地表凹陷速度模型如图1a所示,3层的速度分别为3000m/s,4000m/s,4500m/s,曲坐标系下的速度模型如图1b所示。该速度模型x方向与z方向均为2400m,网格数目均为301,网格间距均为8m,笛卡尔坐标系下的贴体网格模型如图2a所示,图2b为其局部放大图。

图1 简单起伏地表凹陷速度模型a 笛卡尔坐标系;b 曲坐标系

图2 贴体网格模型(a)及其局部放大结果(b)

采用主频为25Hz的雷克子波,在地表布置30个均匀分布的震源,炮间距为80m,每炮由301道接收,道间距为8m,总记录时间为2.4s,每道采样间隔为0.6ms。第15炮的单炮记录如图3所示。

对图3所示的地震数据分别进行常规逆时偏移和棱柱波逆时偏移处理。图4为540ms时正向延拓的波场快照,图4a和图4c分别为540ms的曲坐标系下正向延拓的p0和p1波场快照,笛卡尔坐标系下的波场快照分别如图4b和图4d所示,可以看出,p0波场快照中来自水平层的地震波能量(图4a和图4b 空心箭头所示)比来自高陡构造的地震波能量(图4a和图4b实线箭头处)强得多;而p1波场快照中两者的能量更均衡(图4c和图4d)。由此可知,在棱柱波逆时偏移过程中应用p1信息进行成像,相比于传统逆时偏移方法应用p1信息进行成像,高陡构造的能量分布更均衡。图5a和图5b分别为常规逆时偏移和棱柱波逆时偏移的成像结果,从图5中可以看出,采用常规逆时偏移方法得到的高陡构造成像效果较差,采用棱柱波逆时偏移方法得到的水平构造成像效果虽然不及常规逆时偏移方法效果,但高陡构造的成像效果得到了明显改善(图5白色箭头处)。

图3 用于常规逆时偏移和棱柱波逆时偏移的单炮记录

图4 540ms时正向延拓的波场快照a 曲坐标系下的p0波场;b 笛卡尔坐标系下的p0波场;c 曲坐标系下的p1波场;d 笛卡尔坐标系下的p1波场

图5 简单起伏地表凹陷模型偏移结果a 常规逆时偏移;b 棱柱波逆时偏移

2.2 起伏地表火山岩模型

本文根据川西北部龙门山前带地区实际模型建立了起伏地表火山岩速度模型。该模型中广泛发育具有高陡倾角的火山通道构造。笛卡尔坐标系下的起伏地表火山岩速度模型如图6a所示,速度范围为2800～6400m/s,曲坐标系下的速度模型如图6b所示。该速度模型x方向为10km,z方向为5km,x方向网格点数为1001,z方向网格点数为501,网格间距均为10m。采用主频为25Hz的雷克子波,在地表布置125个均匀分布的震源,炮间距为80m,每炮由1251道接收,总记录时间为4.4s,每道采样间隔为0.5ms。图7a和图7b分别为采用常规逆时偏移和采用棱柱波逆时偏移得到的成像结果。从图7a中可以看出,起伏地表条件下常规逆时偏移可对低倾角的火山岩体及沉积层清晰成像,但对于地震、地质解释非常重要的火山通道构造并没有得到较好的刻画。对比图7a和图7b发现,起伏地表棱柱波逆时偏移方法在成像时增强了棱柱波的能量,因此可以对含高陡倾角的火山岩通道清晰成像,为地震、地质解释提供了直观的可靠依据。

图6 起伏地表火山岩速度模型a 笛卡尔坐标系;b 曲坐标系

图7 起伏地表火山岩模型偏移结果a 常规逆时偏移;b 棱柱波逆时偏移

3 讨论

利用本文方法进行实际数据处理时,不需要进行高程静校正,只需根据地表高程函数生成贴体网格,本文方法与常规逆时偏移一样都需要较为准确的偏移速度场。在应用本文方法进行偏移前,首先采用常规偏移方法得到低倾角沉积层成像结果,输入的数据包括炮域地震数据、偏移速度场、贴体网格文件、观测系统文件和常规成像结果。虽然本文方法能够改善高陡构造的成像结果,但损失了低倾角沉积层的成像效果,因此,在地震、地质解释时,需要与传统成像方法的成像结果进行联合分析。

4 结论与认识

地表起伏剧烈、高陡构造广泛发育导致山前带地区成像困难,针对这一难题,将基于贴体网格的曲坐标系声波方程与棱柱波逆时偏移技术结合,提出并实现了基于贴体网格剖分的曲坐标系棱柱波逆时偏移成像方法,该方法克服了山前带地表剧烈起伏对地震成像的影响,同时充分利用包含大量高陡构造信息的棱柱波进行成像,改善了山前带高陡构造的成像效果。