作者简介: 袁茂林(1990-),男,硕士,2016年获中国石油大学(华东)工学硕士学位,主要从事高斯束正演与反演成像工作。 Email:yml3290432@163.com
Born近似描述了介质扰动引起的一次散射现象,是线性化地震反演的重要基础。笔者基于Born近似,利用高斯束传播算子表征格林函数,推导得到基于地下慢度扰动的一次反射波场理论公式,称这种反射波场的模拟为高斯束线性正演。在编程实现方法的基础上,通过对绕射体模型、复杂砂砾岩体模型以及起伏地表模型进行数值模拟,并对模拟产生的多炮记录进行偏移试验,结果表明:①高斯束线性正演能够准确地模拟出不含直达波和多次波的一次反射波信息;②对高斯束线性正演产生的多炮记录进行偏移能够恢复模型构造特征;③高斯束线性正演易于实现面向目标的数值模拟。
Born approximation describes the primary scattering response of perturbed medium,which is the basis of linearized seismic inversion.Based on Born approximation,we derive the primary reflected wavefiled propagation formula by expressing the Green's function as a Gaussian-beam summation.The simulation of primary reflected waves is called Gaussian-beam linear forward modeling.On this basis,we perform numerical test by using the diffractors model,the glutenite model and the irregular surface model, and compare our results with those generated by wave equation modeling.Besides,we implement migration with corresponding synthetic seismograms generated by Gaussian-beam linear forward modeling.Numerical results indicate:First,Gaussian-beam linear forward modeling can accurately produce recorded data with multiple and the direct wave muted;Second,migration with synthetic seismograms generated by Gaussian-beam linear forward modeling can image corresponding subsurface structures;Third,target-oriented modeling can be carried out naturally.
基于Born近似的正演模拟方法(Born正演)将慢度扰动和背景速度场作为输入, 得到一次反射的地震数据, 这部分能量在勘探地震中占有主要位置, 因此, Born正演是线性化地震反演的重要基础[1,2]。目前, 国内外学者广泛而深入地研究了基于Born近似的线性正演模拟方法, 根据传播算子类型可以大致分为三类:射线类, 单程波动方程类以及双程波动方程类。de Hoop和Bleistein[3]以及Burridge等[4]先后通过等时线叠加实现了各向异性介质的Born正演算法。Jaramillo和Bleistein[5]则基于等时线叠加理论, 推导了声波介质下的Born正演模拟算子。Nemeth等[6]以及黄建平等[7,8]先后将射线类Born正演算子运用到线性波形反演中实现了保幅成像。Kü hl和Sacchi[9]则将裂步双平方根Born正演算子应用到线性波形反演框架下, 实现了AVA/AVP反演。Kaplan等[10]通过单程波动方程的裂步算子来计算适应横向变速的Green函数, 进而实现了Born正演算法。Wu和de Hoop[11]在波动方程Green函数解法的基础上, 借助局部Born近似等数学手段, 发展了广义屏算法。周华敏等[12]发展了稳定Born近似广义屏正演算子并将其运用到线性波形反演中, 解决了深层目标体保幅成像的难题。Wang等[13]则采用傅里叶有限差分算子求解Born散射波场。此外, 基于双程波动方程的Born正演算子[14,15,16]因其较高的精度被广泛地应用于线性波形反演框架中, 以实现“ 三高” 成像。需要指出的是, 射线类Born正演算子模拟地震波场的精度有限, 而波动方程类Born正演算子虽然精度较高, 但往往不够灵活。
高斯束传播算子通过计算复值动力学射线追踪参量, 使其振幅为处处正则的, 不存在波场的奇异性区域, 被认为是Kirchhoff传播算子和单程波动方程传播算子之间的折中方法。Popov[17]以及
笔者基于Born近似实现了高斯束线性正演模拟, 主要工作包括:①介绍了高斯束传播算子的基本理论及其数值求解方法; ②基于Born近似, 采用高斯束叠加积分表征格林函数, 推导了高斯束线性正演算子, 并提出了其算法实现流程; ③采用绕射体模型、复杂砂砾岩体模型以及起伏地表模型进行高斯束线性正演模拟, 并且将模拟产生的多炮记录进行偏移测试, 除此之外, 将高斯束线性正演模拟结果与波动方程正演模拟结果进行了对比分析, 验证了本文方法的正确性、有效性和适应性; ④讨论了高斯束线性正演模拟实现面向目标正演的策略, 并分析了相应的数值试验结果。
模拟地震波场的高斯束传播算子, 可以在射线中心坐标系中讨论。射线中心坐标系如图1所示, 假设M为中心射线Ω 上一点, e1、e2和e3为M点处射线中心坐标系(q1, q2, s)的基矢量, 其中, e3为单位正切矢量, e1和e2垂直于中心射线, 并且使射线中心坐标系为正则的。q1和q2分别为N点在射线中心坐标系中沿基矢量e1和e2的坐标, s为N点沿中心射线到参考点s0的距离。
高斯束可以描述局部波场传播, 是波动方程的高频渐进解。在N点, 三维高斯束的公式表示为[28]
式中:qT=(q1, q2); 矩阵M(s)=
根据高斯束的初始位置和初始方向, 利用运动学射线追踪方程组来求取中心射线的路径和走时:
式中:xi(s)为物理空间坐标(x, y, z)中的射线坐标, pi(s)为射线慢度矢量, τ 为沿射线的走时, dτ 为求积步长。
P(s)、Q(s)满足三维动力学射线追踪方程组:
V(s)为速度场沿二阶导数矩阵:
P(s)、Q(s)初始值可以设定为[29]:
式中:w0为初始束宽; ω r为参考频率, 与地震数据的频带宽度有关。
根据运动学和动力学射线追踪求得的走时、振幅等信息, 可以求取中心射线附近有效宽度内的波场。有效宽度可以按照下述方式来定义:由于高斯束的振幅随着离中心射线的距离而衰减, 可以定义振幅大于中心射线振幅1%的范围为高斯束的有效宽度。为了使高斯束成为波动方程较为精确的高频渐进解, 需要保证速度场在束宽范围内缓慢变化, 因此, 在高斯束的传播过程中, 束宽越窄, 地震波场越精确。但是, 如果初始束宽过小, 高斯束在传播过程中束宽会迅速增大, 降低地震波场的精度。通常情况下, 当初始束宽大于一个波长时, 束宽不会急剧地增大[29], 但是, 如果初始束宽过大, 地震波场则不准确。
地震散射总波场可以由Lippmann-Swinger方程表示, 包含直达波场, 一次散射波场以及多次散射波场, 用Neumann级数展开可以分别表示为零阶项、一阶项以及高阶项。如图2所示, 在背景慢度场中存在一个扰动区域, 扰动慢度场与背景慢度场之间存在轻微差异。那么, 地面记录波场可以由地下所有扰动场产生的一次散射波场叠加形成, 即Born散射波场[1]:
上式中, G(x, xs, ω )为震源格林函数, 描述从震源xs=(xs, ys, 0)到地下散射点x=(x, y, z)的地震波场; G(xr, x, ω )为接收点格林函数, 描述从散射点x到接收点xr=(xr, yr, 0)的地震波场; m(x)=
格林函数G(x, x0, ω )由高斯束uGB(x, x0, p, ω )的叠加积分来表征[29]:
式(7)所示的格林函数描述由点源x0到地下任意一点x的点源波场, 可以由点源x0处按不同初始射线参数p=(px, py, pz)出射的高斯束在点x处叠加形成。初始射线参数间隔必须足够小以保证高斯束叠加积分得到的格林函数是一个带限的, 且无假频的点源波场。上述条件意味着两条邻近高斯束在点源附近点的走时变化必须要小于半个周期。初始射线参数间隔为[29]
式中, ω h为最高频率, 与地震数据的频带宽度有关。
因此, 式(6)中的震源格林函数和接收点格林函数可分别表示为
G(x, xs, ω )≈
G(xr, x, ω )≈
将式(9)和式(10)代入式(6), 得到高斯束线性正演公式:
式(11)描述了由高斯束模拟的一次散射波场在接收点处叠加形成的记录波场。可以看出, 基于Born近似的高斯束正演是一种关于慢度扰动的线性化算子。
根据式(11)所表达的高斯束线性正演算子, 设计了如图3所示的工作流程来实施数值试验。从此工作流程图中可以看出:①高斯束线性正演将背景速度场和慢度扰动作为输入, 得到合成地震记录; ②高斯束线性正演模拟流程描述了震源波场向下延拓, 当遇到慢度扰动场时, 产生一次散射波场被检波器接收并叠加的过程。由此可见, 高斯束线性正演能够模拟一次散射产生的地震波场。
为了验证文中提出的高斯束线性正演算子的正确性、有效性和适应性, 在编程实现高斯束线性正演算法的基础上, 采用绕射体模型、复杂砂砾岩体模型以及起伏地表模型进行数值模拟, 并对模拟产生的合成地震记录进行偏移测试。
为了验证文中高斯束线性正演算子的正确性和有效性, 首先采用绕射体模型进行试算。绕射体速度模型如图4a所示, 模型网格点数为101× 81, 纵横向网格间距均为10 m, 背景速度场为2 000 m/s的常速介质, 其中包含3个速度均为4 000 m/s的绕射体, 绕射体直径均为10 m。图4b为慢度扰动模型。
此模型可以测试高斯束传播算子模拟一次散射波的效果, 按照图3所示的高斯束线性正演流程, 可以得到图5所示的合成地震记录, 震源在地表x=500 m处, 检波器间距为10 m, 总共101道接收, 每道1 000个采样点, 采样间隔为1 ms。从图5可以看出:①高斯束线性正演模拟的地震记录同相轴走时准确, 绕射体响应特征清晰; ②高斯束线性正演模拟的地震记录中不存在直达波和多次散射波, 这是因为高斯束线性正演基于Born近似, 只能模拟一次散射波场。以上数值试验和理论分析验证了文中高斯束线性正演算子的正确性和有效性。
为了从另一个角度说明高斯束线性正演方法的正确性和有效性, 接下来对高斯束线性正演模拟产生的多炮地震记录进行高斯束偏移[30]试算, 输入数据为51炮高斯束线性正演模拟的地震记录, 第一炮的位置为x=0处, 炮间隔为20 m, 101道接收, 每道1 000个采样点, 采样间隔为0.1 ms。从图6可以看出:3个绕射体的成像位置正确, 且聚焦性较好, 但由于偏移算子与高斯束线性正演算子不是互逆的关系, 因此偏移结果只是真实慢度扰动的模糊成像。
为了进一步验证高斯束线性正演算法的适应性, 笔者对复杂的砂砾岩体模型进行了试处理。砂砾岩体速度模型如图7a所示, 此模型速度变化范围较大, 从2 900 m/s到6 000 m/s, 包含高陡倾等复杂构造, 模型网格点数为550× 200, 纵横向网格间距均为10 m。本文数值模拟将平滑后的速度场作为背景速度场。图7b为慢度扰动模型。
基于高斯束传播算子, 可以模拟出如图8所示的点源波场, 点源位置在x=2 740 m处, 可以看到第一层均匀介质两侧的波场出现了扰动现象(黑色箭头所指), 这是因为高斯束经过浅层速度扰动界面(黑色方框所示)时, 射线路径发生回转后高斯束波场叠加所致。总体上, 高斯束传播算子模拟的点源波场不存在奇异性区域。
图9a为高斯束线性正演模拟的单炮记录, 震源在地表x=2 740 m处, 总共550道接收, 检波器间距为10 m, 每道2 000个采样点, 采样间隔为1 ms。图9b为波动方程有限差分正演模拟记录, 对比图9a和图9b, 可以看到:①高斯束线性正演模拟的地震记录同相轴形态及位置与波动方程正演模拟结果相近, 有效地描述了复杂模型中反射体和绕射体的响应特征; ②因为高斯束线性正演模拟的地震波场是由地下所有扰动源产生的一次散射波场叠加形成, 所以高斯束线性正演结果中存在更多的绕射“ 画弧” 现象(黑色箭头所示), 此现象可以解释为更丰富的绕射波信息; ③由于高斯束模拟地震波场的精度有限, 导致合成地震记录中波场信息并不完全(黑色圆圈所示)。
图10a、b分别为高斯束偏移和逆时偏移[31]测试结果, 输入数据为90炮高斯束线性正演模拟的地震记录, 第一炮位置为x=70 m处, 炮间隔为60 m, 每炮550道接收, 检波器间距为10 m, 每道2 000个采样点, 采样间隔为1 ms。从图10可以看出:无论是高斯束偏移结果还是逆时偏移结果, 各个反射层成像位置准确, 较好地恢复了高陡倾等复杂构造特征。但是, 两个偏移成像中存在不同程度的噪声。一方面是偏移噪声所致; 另一方面, 高斯束线性正演产生数据噪声, 逆时偏移对数据信噪比的要求更高, 产生更多的偏移假象。由于偏移算子与高斯束线性正演算子不是互逆的关系, 因此偏移结果只是真实慢度扰动的模糊成像。从另一个角度反映了本文高斯束线性正演模拟方法对复杂模型具有较强的适应性。
为了测试高斯束线性正演方法面向目标正演的效果, 本文对砂砾岩体模型进行了数值模拟, 目标区域设定在x:1 600~2 800 m, z:1 000~1 600 m(图7a中虚线方框所示区域)。在面向目标正演时, 只需将工作流程(见图3)中“ 判断是否存在慢度扰动m(x)?” 改为“ 判断在目标区域是否存在慢度扰动m(x)?” 。
图11a为面向此目标区域进行高斯束线性正演模拟得到的单炮记录, 震源在地表处x=2 740 m, 总共550道接收, 检波器间距为10 m, 每道2 000个采样点, 采样间隔为1 ms。相较于图9a, 图11a显示了对目标区域进行正演模拟的波场信息。图11b为偏移测试结果, 输入数据为面向此目标区域进行高斯束线性正演模拟得到的90炮合成地震记录, 第一炮位置为x=70 m处, 炮间隔为60 m, 每炮550道接收, 检波器间距为10 m, 每道2 000个采样点, 采样间隔为1 ms。从图11b可以看出, 成像结果准确地恢复了目标区域的地质构造特征, 说明了高斯束线性正演方法面向目标正演模拟的正确性、有效性和适应性。在时移地震的实际应用中, 由于背景速度场不变化, 那么基于高斯束传播算子的格林函数只需要计算一次, 针对随时间变化的偏移成像剖面, 可以面向目标储层进行正演模拟, 从而检测目标储层随时间的变化情况。这也是高斯束线性正演模拟相较于常规波动方程有限差分正演模拟的优势之处。
起伏地表模型为了测试本文高斯束线性正演模拟算法对起伏地表的适应性, 将图7a所示的复杂砂砾岩体模型从分别采用黄建平等[32]提出的双复杂条件下高斯束偏移方法和逆时偏移方法进行成像测试, 结果如图14所示, 输入数据为61炮高斯束线性正演模拟的地震记录, 第一炮位置为x=340 m处, 炮间隔为80 m, 每炮550道接收, 检波器间距为10 m, 每道2 000个采样点, 采样间隔为1 ms。从图14可以看出:①由于高斯束线性正演模拟的地震记录中不包含直达波信息, 所以偏移结果中没有起伏地表成像。②无论是高斯束偏移结果还是逆时偏移结果, 各个反射层成像准确, 较好地恢复了高陡倾等复杂构造特征。但由于偏移算子与高斯束线性正演算子不是互逆的关系, 因此偏移结果只是真实慢度扰动的模糊成像。以上分析说明了高斯束线性正演模拟方法对起伏地表模型具有较好的适应性。
笔者提出了一种新的正演模拟方法, 称为高斯束线性正演模拟。基于Born近似, 通过高斯束传播算子模拟点源波场, 构建了高斯束线性正演理论体系及算法流程, 数值测试验证了高斯束线性正演模拟方法的正确性、有效性和适应性, 并且得到以下几点认识:
1) 高斯束线性正演方法以背景速度场和慢度扰动作为输入数据, 模拟得到合成地震记录。由于高斯束线性正演基于Born近似, 因此不能模拟直达波以及多次波等波场信息。
2) 相较于波动方程有限差分正演模拟, 高斯束线性正演模拟方法具有灵活且易于实现面向目标正演模拟的特点, 尤其适用于时移地震中。在时移地震中, 由于背景速度场是不变的, 变化的是偏移成像, 因此, 可以应用高斯束线性正演方法对随时间变化的目标储层进行正演模拟。
3) 常规高斯束偏移仅仅对地下慢度扰动进行模糊成像, 其振幅往往是欠估计的规高斯束偏移仅仅对地下慢度扰动进行模糊成像[32-34]。在实际地震资料处理中, 可以将偏移速度场视为背景速度, 将偏移成像视为慢度扰动。因此, 可以发展一对互为共轭关系的高斯束线性正演/偏移算子, 将其运用到线性波形反演框架下, 实现真振幅成像; 当然, 也可以将高斯束线性正演/偏移算子运用到数据规则化以及去噪处理中。
致谢:感谢中国石油大学(华东)SWPI(Seismic Wave Propagation and Imaging)课题组对本项研究的支持。
The authors have declared that no competing interests exist.
[1] |
|
[2] |
|
[3] |
|
[4] |
|
[5] |
|
[6] |
|
[7] |
|
[8] |
|
[9] |
|
[10] |
|
[11] |
|
[12] |
|
[13] |
|
[14] |
|
[15] |
|
[16] |
|
[17] |
|
[18] |
|
[19] |
|
[20] |
|
[21] |
|
[22] |
|
[23] |
|
[24] |
|
[25] |
|
[26] |
|
[27] |
|
[28] |
|
[29] |
|
[30] |
|
[31] |
|
[32] |
|
[33] |
|
[34] |
|