石油物探  2022, Vol. 61 Issue (3): 512-520  DOI: 10.3969/j.issn.1000-1441.2022.03.013
0
文章快速检索     高级检索

引用本文 

冉崎, 陈康, 韩嵩, 等. 起伏地表棱柱波逆时偏移成像方法[J]. 石油物探, 2022, 61(3): 512-520. DOI: 10.3969/j.issn.1000-1441.2022.03.013.
RAN Qi, CHEN Kang, HAN Song, et al. Reverse time migration imaging method of prismatic waves for surface topography[J]. Geophysical Prospecting for Petroleum, 2022, 61(3): 512-520. DOI: 10.3969/j.issn.1000-1441.2022.03.013.

基金项目

国家自然科学基金项目(42174138, 42074133)、中国石化地球物理重点实验室开放基金项目(wtyjy-wx2017-01-04)、中石油重大科技专项(ZD2019-183-003)、中央高校基础研究基金项目(19CX02010A)和学校人才引进费项目(20180041)共同资助

第一作者简介

冉崎(1973—), 男, 硕士, 高级工程师, 长期从事地震勘探综合解释、应用及相关研究工作。Email: ranqi@petrochina.com.cn

通信作者

曲英铭(1990—), 男, 博士, 副教授, 博士生导师, 主要从事地震波传播、成像与反演等方面的研究工作。Email: quyingming@upc.edu.cn

文章历史

收稿日期:2021-02-10
起伏地表棱柱波逆时偏移成像方法
冉崎1, 陈康1, 韩嵩1, 马博2, 刘畅2, 曲英铭2    
1. 中国石油天然气股份有限公司西南油气田分公司勘探开发研究院, 四川成都 610021;
2. 中国石油大学(华东)地球科学与技术学院地球物理系, 山东青岛 266580
摘要:山前带地区地表起伏剧烈、地下构造复杂。当地表速度横向变化剧烈或地表高程变化很大时, 基于高程静校正的常规偏移成像方法难以消除起伏地形对偏移的影响。针对山前带高陡构造的成像, 提出了一种曲坐标系下起伏地表棱柱波逆时偏移成像方法。将棱柱波逆时偏移成像方法引入曲坐标系下的声波方程中, 推导出曲坐标系棱柱波逆时偏移方程, 利用棱柱波提高山前带高陡构造的成像精度, 并采用贴体网格剖分技术及坐标变换技术, 在曲坐标系下进行偏移计算。简单起伏地表凹陷模型和起伏地表火山岩模型试算结果表明, 坐标变换后的起伏地表被映射为水平地表, 曲坐标系下起伏地表棱柱波逆时偏移成像方法克服了起伏地表对成像的影响, 得到了准确的山前带高陡构造成像结果。
关键词起伏地表    棱柱波    逆时偏移    贴体网格    高陡构造    曲坐标系    
Reverse time migration imaging method of prismatic waves for surface topography
RAN Qi1, CHEN Kang1, HAN Song1, MA Bo2, LIU Chang2, QU Yingming2    
1. Research Institute of Exploration and Development, PetroChina Southwest Oil & Gas Field Company, Chengdu 610021, China;
2. School of Geosciences, China University of Petroleum (East China), Qingdao 266580, China
Abstract: The surface topography of a piedmont zone fluctuates violently, and its underground structure is complex. When the near-surface lateral velocity or elevation changes drastically, the conventional migration imaging methods based on static correction cannot overcome the influence of undulating surface topography on the migration imaging. To overcome the imaging problem of steeply dipping structures in piedmont areas, a prismatic wave reverse-time migration method based on body-fitted grids is proposed in the curved coordinate system, which makes complete use of prismatic waves to improve the imaging accuracy of steeply dipping structures in the foothill areas. The method was implemented in the curved coordinate system using body-fitted grids and coordinate transformation technology. Numerical examples on simple topographic sag and topographic volcanic rock models demonstrated that the irregular surface was transformed to a horizontal one by coordinate transformation. The proposed prismatic wave reverse-time migration method was able to produce accurate and steeply dipping structural images of the foothill areas and overcome the influence of undulating surface topography.
Keywords: surface topography    prismatic wave    reverse time migration    body-fitted grid    steeply dipping structures    curved coordinate system    

山前带地区地表起伏剧烈, 复杂的高陡构造发育。剧烈变化的起伏地表给常规地震勘探工作带来了极大的挑战。常规起伏地表地震资料处理方法包括静校正、基准面校正以及直接成像法[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 曲坐标系下的声波方程

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

$ \left\{\begin{array}{l} \frac{\partial \sigma}{\partial t}=-\rho v_{\mathrm{P}}^{2}\left(\frac{\partial v_{x}}{\partial x}+\frac{\partial v_{z}}{\partial z}\right) \\ \frac{\partial v_{x}}{\partial t}=\frac{1}{\rho} \frac{\partial \sigma}{\partial x} \\ \frac{\partial v_{z}}{\partial t}=\frac{1}{\rho} \frac{\partial \sigma}{\partial z} \end{array}\right. $ (1)

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

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

$ \left\{\begin{array}{l} \xi_{x x}+\xi_{z z}=P(x, z) \\ \eta_{x x}+\eta_{z z}=Q(x, z) \end{array}\right. $ (2)

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

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

$ \left\{\begin{array}{l} x=x(\xi, \eta) \\ z=z(\xi, \eta) \end{array}\right. $ (3)

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

$ \left\{\begin{array}{l} 1=x_{\xi} \xi_{x}+x_{\eta} \eta_{x} \\ 0=z_{\xi} \xi_{x}+z_{\eta} \eta_{x} \end{array}\right. $ (4)
$ \left\{\begin{array}{l} 0=x_{\xi} \xi_{z}+x_{\eta} \eta_{z} \\ 1=z_{\xi} \xi_{z}+z_{\eta} \eta_{z} \end{array}\right. $ (5)

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

$ \left\{\begin{array}{l} \xi_{x}=\frac{1}{J} z_{\eta} \\ \xi_{z}=\frac{-1}{J} x_{\eta} \\ \eta_{x}=\frac{-1}{J} z_{\xi} \\ \eta_{z}=\frac{1}{J} x_{\xi} \\ J=x_{\xi} z_{\eta}-x_{\eta} z_{\xi} \end{array}\right. $ (6)

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

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

$ \left\{\begin{array}{l} \frac{\partial \sigma}{\partial t}=-\rho v_{\mathrm{P}}^{2}\left(\frac{\partial v_{x}}{\partial \xi} \xi_{x}+\frac{\partial v_{x}}{\partial \eta} \eta_{x}+\frac{\partial v_{z}}{\partial \xi} \xi_{z}+\frac{\partial v_{z}}{\partial \eta} \eta_{z}\right) \\ \frac{\partial v_{x}}{\partial t}=\frac{1}{\rho}\left(\frac{\partial u}{\partial \xi} \xi_{x}+\frac{\partial u}{\partial \eta} \eta_{x}\right) \\ \frac{\partial v_{z}}{\partial t}=\frac{1}{\rho}\left(\frac{\partial u}{\partial \xi_{z}} \xi_{z}+\frac{\partial u}{\partial \eta} \eta_{z}\right) \end{array}\right. $ (7)
1.2 棱柱波逆时偏移

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

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

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

$ \begin{gathered} m_{\text {mig }}\left(x \mid x_{\mathrm{s}}\right)=\sum\limits_{w} \sum\limits_{g} \omega^{2} W^{*}(\omega) G^{*}\left(x \mid x_{\mathrm{s}}\right) \\ G^{*}\left(x \mid x_{\mathrm{g}}\right) d\left(x_{\mathrm{g}} \mid x_{\mathrm{s}}\right) \end{gathered} $ (8)

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

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

$ \begin{aligned} G_{1}\left(x \mid x_{\mathrm{s}}\right)=& \int_{x^{\prime}} \omega^{2} m_{1}\left(x^{\prime}\right) G_{0}\left(x^{\prime} \mid x_{\mathrm{s}}\right) G_{0} \cdot \\ &\left(x^{\prime} \mid x\right) \mathrm{d} x^{\prime} \end{aligned} $ (9)

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

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

$ \begin{gathered} \sum\limits_{\omega} \omega^{2} W^{*}(\omega) G_{0}^{*}\left(x \mid x_{\mathrm{s}}\right) G_{1}^{*}\left(x \mid x_{\mathrm{g}}\right) d_{1}\left(x_{\mathrm{g}} \mid x_{\mathrm{s}}\right) \\ =\sum\limits_{\omega} \omega^{2} W^{*}(\omega) \int\limits_{x^{\prime}} \omega^{2} m_{1}\left(x^{\prime}\right) G_{0}^{*}\left(x^{\prime} \mid x_{\mathrm{s}}\right) \cdot \\ G_{0}^{*}\left(x^{\prime} \mid x\right) \mathrm{d} x^{\prime} G_{0}^{*}\left(x \mid x_{\mathrm{g}}\right) d_{2}\left(x_{\mathrm{g}} \mid x_{\mathrm{s}}\right) \\ =\sum\limits_{\omega} \omega^{2} P_{1}^{*}\left(x \mid x_{\mathrm{s}}\right) Q_{0}\left(x \mid x_{\mathrm{s}}\right) \end{gathered} $ (10)

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

$ \left\{\begin{array}{l} {\left[\nabla^{2}+\omega^{2} v_{0}^{-2}(x)\right] P_{0}(x)=W(\omega) \delta\left(x-x_{\mathrm{s}}\right)} \\ {\left[\nabla^{2}+\omega^{2} v_{0}^{-2} 0(x)\right] P_{1}(x)=\omega^{2} m_{1}(x) P_{0}(x)} \\ {\left[\nabla^{2}+\omega^{2} v_{0}^{-2} 0(x)\right] Q_{0}\left(x \mid x_{\mathrm{s}}\right)=} \\ \quad d_{2}\left(x_{\mathrm{g}} \mid x_{\mathrm{s}}\right) \delta\left(x-x_{\mathrm{g}}\right) \end{array}\right. $ (11)

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

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

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

$ \left\{\begin{array}{l} \rho \frac{\partial v_{x}}{\partial t}=\frac{\partial p}{\partial \xi} \frac{\partial \xi}{\partial x}+\frac{\partial p}{\partial \eta} \frac{\partial \eta}{\partial x} \\ \rho \frac{\partial v_{z}}{\partial t}=\frac{\partial p}{\partial \xi} \frac{\partial \xi}{\partial z}+\frac{\partial p}{\partial \eta} \frac{\partial \eta}{\partial z} \\ \frac{1}{v_{\mathrm{P}}^{2}} \frac{\partial p}{\partial t}=\rho\left(\frac{\partial v_{x}}{\partial \xi} \frac{\partial \xi}{\partial x}+\frac{\partial v_{x}}{\partial \eta} \frac{\partial \eta}{\partial x}+\frac{\partial v_{z}}{\partial \xi} \frac{\partial \xi}{\partial z}+\frac{\partial v_{z}}{\partial \eta} \frac{\partial \eta}{\partial z}\right) \end{array}\right. $ (12)

速度扰动可表示为:

$ v_{\mathrm{P}}=v_{\mathrm{P} 0}+\delta v_{\mathrm{P}} $ (13)

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

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

$ \left\{\begin{array}{l} \rho \frac{\partial v_{x 0}}{\partial t}=\frac{\partial p}{\partial \xi} \frac{\partial \xi}{\partial x}+\frac{\partial p}{\partial \eta} \frac{\partial \eta}{\partial x} \\ \rho \frac{\partial v_{z 0}}{\partial t}=\frac{\partial p}{\partial \xi} \frac{\partial \xi}{\partial z}+\frac{\partial p}{\partial \eta} \frac{\partial \eta}{\partial z} \\ \frac{1}{v_{\mathrm{P} 0}^{2}} \frac{\partial p_{0}}{\partial t}=\rho\left(\frac{\partial v_{x 0}}{\partial \xi} \frac{\partial \xi}{\partial x}+\frac{\partial v_{x 0}}{\partial \eta} \frac{\partial \eta}{\partial x}+\right. \\ \ \ \ \ \ \ \ \ \frac{\partial v_{z 0}}{\partial \xi} \frac{\partial \xi}{\partial z}+\frac{\partial v_{z 0}}{\partial \eta} \frac{\partial \eta}{\partial z} \end{array}\right. $ (14)

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

$ \frac{1}{v_{\mathrm{P}}^{2}} \approx \frac{1}{v_{\mathrm{P} 0}^{2}}-\frac{2 \delta v_{\mathrm{P}}}{v_{\mathrm{P} 0}^{3}} $ (15)

将(15)式代入(12)

式中, 并减去(14)式, 可得:

$ \left\{\begin{array}{l} \rho \frac{\partial \delta v_{x}}{\partial t}=\frac{\partial \delta p}{\partial \xi} \frac{\partial \xi}{\partial x}+\frac{\partial \delta p}{\partial \eta} \frac{\partial \eta}{\partial x} \\ \rho \frac{\partial \delta v_{z}}{\partial t}=\frac{\partial \delta p}{\partial \xi} \frac{\partial \xi}{\partial z}+\frac{\partial \delta p}{\partial \eta} \frac{\partial \eta}{\partial z} \\ \frac{1}{v_{\mathrm{P} 0}^{2}} \frac{\partial \delta p}{\partial t}=\rho\left(\frac{\partial \delta v_{x}}{\partial \xi} \frac{\partial \xi}{\partial x}+\frac{\partial \delta v_{x}}{\partial \eta} \frac{\partial \eta}{\partial x}\right)+\rho\left(\frac{\partial \delta v_{z}}{\partial \xi} \frac{\partial \xi}{\partial z}+\right. \\ \ \ \ \ \ \ \ \ \left.\frac{\partial \delta v_{z}}{\partial \eta} \frac{\partial \eta}{\partial z}\right)+\frac{2 \partial v_{\mathrm{P}}}{v_{\mathrm{P} 0}^{3}} \frac{\partial p_{0}}{\partial t}+O\left(\partial v_{\mathrm{P}}\right) \end{array}\right. $ (16)

反射系数模型可表示为:

$ m\left(v_{\mathrm{P}}\right)=\frac{\delta v_{\mathrm{P}}}{v_{\mathrm{P} 0}} $ (17)

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

$ \left\{\begin{array}{l} \rho \frac{\partial \delta v_{x}}{\partial t}=\frac{\partial \delta p}{\partial \xi} \frac{\partial \xi}{\partial x}+\frac{\partial \delta p}{\partial \eta} \frac{\partial \eta}{\partial x} \\ \rho \frac{\partial \delta v_{z}}{\partial t}=\frac{\partial \delta p}{\partial \xi} \frac{\partial \xi}{\partial z}+\frac{\partial \delta p}{\partial \eta} \frac{\partial \eta}{\partial z} \\ \frac{1}{v_{\mathrm{P} 0}^{2}} \frac{\partial \delta p}{\partial t}=\rho\left(\frac{\partial \delta v_{x}}{\partial \xi} \frac{\partial \xi}{\partial x}+\frac{\partial \delta v_{x}}{\partial \eta} \frac{\partial \eta}{\partial x}\right)+\rho\left(\frac{\partial \delta v_{z}}{\partial \xi} \frac{\partial \xi}{\partial z}+\right. \\ \ \ \ \ \ \ \ \ \left.\frac{\partial \delta v_{z}}{\partial \eta} \frac{\partial \eta}{\partial z}\right)+\frac{2 m\left(v_{\mathrm{P}}\right)}{v_{\mathrm{P} 0}^{2}} \frac{\partial p_{0}}{\partial t} \end{array}\right. $ (18)

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

$ \left\{\begin{array}{l} \rho \frac{\partial v_{x 1}}{\partial t}=\frac{\partial p_{1}}{\partial \xi} \frac{\partial \xi}{\partial x}+\frac{\partial p_{1}}{\partial \eta} \frac{\partial \eta}{\partial x} \\ \rho \frac{\partial v_{z 1}}{\partial t}=\frac{\partial p_{1}}{\partial \xi} \frac{\partial \xi}{\partial z}+\frac{\partial p_{1}}{\partial \eta} \frac{\partial \eta}{\partial z} \\ \frac{1}{v_{\mathrm{P} 0}^{2}} \frac{\partial p_{1}}{\partial t}=\rho\left(\frac{\partial v_{x 1}}{\partial \xi} \frac{\partial \xi}{\partial x}+\frac{\partial v_{x 1}}{\partial \eta} \frac{\partial \eta}{\partial x}\right)+\rho\left(\frac{\partial v_{z 1}}{\partial \xi} \frac{\partial \xi}{\partial z}+\right. \\ \ \ \ \ \ \ \ \ \left.\frac{\partial v_{z 1}}{\partial \eta} \frac{\partial \eta}{\partial z}\right)+\frac{2 I}{v_{\mathrm{P} 0}^{2}} \frac{\partial p_{0}}{\partial t} \end{array}\right. $ (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 用于常规逆时偏移和棱柱波逆时偏移的单炮记录

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

图 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 结论与认识

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

参考文献
[1]
BERRYHILL J R. Wave equation datuming[J]. Geophysics, 1979, 44(8): 1329-1344. DOI:10.1190/1.1441010
[2]
BERRYHILL J R. Wave equation datuming before stack[J]. Geophysics, 1984, 49(11): 2064-2066. DOI:10.1190/1.1441620
[3]
李录明, 罗省贤. 波场延拓表层模型校正[J]. 石油地球物理勘探, 2001, 36(5): 572-583.
LI L M, LUO S X. Surface model correction by wave-field continuation[J]. Oil Geophysical Prospecting, 2001, 36(5): 572-583. DOI:10.3321/j.issn:1000-7210.2001.05.008
[4]
曲英铭. 起伏地表直接成像技术研究进展[J]. 石油物探, 2019, 58(5): 625-644.
QU Y M. Research progress of topographic imaging methods[J]. Geophysical Prospecting for Petroleum, 2019, 58(5): 625-644. DOI:10.3969/j.issn.1000-1441.2019.05.001
[5]
WIGGINS J W. Kirchhoff integral extrapolation and migration of nonplanar data[J]. Geophysics, 1984, 49(8): 1239-1248. DOI:10.1190/1.1441752
[6]
JAGER C, HERTWECK T, SPINER M. Ture-amplitude Kirchhoff migration from topography[J]. Expanded Abstracts of 73rd Annual Internat SEG Mtg, 2003, 909-913.
[7]
秦宁, 王延光, 杨晓东, 等. 非水平地表高斯束叠前深度偏移及山前带应用实例[J]. 石油地球物理勘探, 2017, 52(1): 81-86.
QIN N, WANG Y G, YANG X D, et al. Gaussian beam prestack depth migration for undulating-surface area in piedmont zone[J]. Oil Geophysical Prospecting, 2017, 52(1): 81-86.
[8]
GRAY S H. Gaussian beam migration of common-shot records[J]. Geophysics, 2005, 70(4): S71-S77. DOI:10.1190/1.1988186
[9]
岳玉波, 李振春, 钱忠平, 等. 复杂地表条件下保幅高斯束偏移[J]. 地球物理学报, 2012, 55(4): 1376-1383.
YUE Y B, LI Z C, QIAN Z P, et al. Amplitude preserved Gaussian beam migration under complex topographic conditions[J]. Chinese Journal of Geophysics, 2012, 55(4): 1376-1383. DOI:10.6038/j.issn.0001-5733.2012.04.033
[10]
黄建平, 袁茂林, 段新意, 等. 一种解耦的起伏地表弹性波高斯束偏移方法[J]. 石油地球物理勘探, 2015, 50(3): 460-468.
HUANG J P, YUAN M L, DUAN X Y, et al. Decoupled elastic Gaussian beam migration for rugged topography[J]. Oil Geophysical Prospecting, 2015, 50(3): 460-468.
[11]
杨继东, 黄建平, 王欣, 等. 复杂地表条件下叠前菲涅尔束偏移方法[J]. 地球物理学报, 2015, 58(10): 3758-3770.
YANG J D, HUANG J P, WANG X, et al. Prestack Fresnel beam migration method under complex topographic conditions[J]. Chinese Journal of Geophysics, 2015, 58(10): 3758-3770. DOI:10.6038/cjg20151026
[12]
徐秀刚, 孙道朋, 刘怀山. 基于起伏地表的小波束叠前深度偏移成像技术[J]. 海洋地质前沿, 2015, 31(12): 58-64.
XU X G, SUN D P, LIU H S. Research of prestack depth migration image-making technique with beamlets for topographic relief areas[J]. Marine Geological Frontiers, 2015, 31(12): 58-64.
[13]
黄建平, 杨继东, 李振春, 等. 基于有效邻域波场近似的起伏地表保幅高斯束偏移[J]. 地球物理学报, 2016, 59(6): 2245-2256.
HUANG J P, YANG J D, LI Z C, et al. An amplitude preserved Gaussian beam migration based on wave filed approximation in effective vicinity under irregular topographical conditions[J]. Chinese Journal of Geophysics, 2016, 59(6): 2245-2256.
[14]
陈林谦. 基于起伏地表单程波波动方程的叠前深度偏移成像技术应用与研究[D]. 北京: 中国石油大学, 2017
CHEN L Q. Application and method study of pre-stack one-way equation depth migration based on irregular topography[D]. Beijing: China University of Petroleum, 2017
[15]
段心标, 王华忠, 白英哲, 等. 基于GPU的三维起伏地表单程波叠前深度偏移[J]. 石油物探, 2016, 55(2): 223-230.
DUAN X B, WANG H Z, BAI Y Z, et al. 3D one-way wave equation pre-stack depth migration from topography based on the acceleration of GPU[J]. Geophysical Prospecting for Petroleum, 2016, 55(2): 223-230. DOI:10.3969/j.issn.1000-1441.2016.02.008
[16]
张宇, 徐升, 张关泉, 等. 真振幅全倾角单程波方程偏移方法[J]. 石油物探, 2007, 46(6): 582-587.
ZHANG Y, XU S, ZHANG G Q, et al. True amplitude full tilt one-way equation migration[J]. Geophysical Prospecting for Petroleum, 2007, 46(6): 582-587. DOI:10.3969/j.issn.1000-1441.2007.06.010
[17]
WHITMORE N. Iterative depth migration by backward time propagation[J]. Expanded Abstracts of 53rd Annual Internat SEG Mtg, 1983, 382-385.
[18]
BAYSAL E D, KOSLOFF D, SHERWOOD J W. Reverse time migration[J]. Geophysics, 1983, 48(1): 1514-1524.
[19]
LOEWENTHAL D, MUFTI I. Reverse time migration in spatial frequency domain[J]. Geophysics, 1983, 48(5): 627-635. DOI:10.1190/1.1441493
[20]
刘红伟, 刘洪, 李博, 等. 起伏地表叠前逆时偏移理论及GPU加速技术[J]. 地球物理学报, 2011, 54(7): 1883-1892.
LIU H W, LIU H, LI B, et al. Pre-stack reverse time migration for rugged topography and GPU acceleration technology[J]. Chinese Journal of Geophysics, 2011, 54(7): 1883-1892. DOI:10.3969/j.issn.0001-5733.2011.07.022
[21]
王延光, 匡斌. 起伏地表叠前逆时深度偏移与并行实现[J]. 石油地球物理勘探, 2012, 47(2): 266-271.
WANG Y G, KUANG B. Pre-stack reverse-time depth migration on rugged topography and parallel computation realization[J]. Oil Geophysical Prospecting, 2012, 47(2): 266-271.
[22]
鲁雁翔, 白超英. 起伏层状介质中曲线交错网格有限差分弹性波逆时偏移成像[J]. 地球物理学进展, 2019, 34(2): 605-614.
LU Y X, BAI C Y. Time reverse migration based on curvilinear stagger-grid finite difference method[J]. Progress in Geophysics, 2019, 34(2): 605-614.
[23]
兰海强, 刘佳, 白志明. VTI介质起伏地表地震波场模拟[J]. 地球物理学报, 2011, 54(8): 2072-2084.
LAN H Q, LIU J, BAI Z M. Wave-field simulation in VTI media with irregular free surface[J]. Chinese Journal of Geophysics, 2011, 54(8): 2072-2084. DOI:10.3969/j.issn.0001-5733.2011.08.014
[24]
李庆洋, 黄建平, 李振春. 起伏地表贴体全交错网格仿真型有限差分正演模拟[J]. 石油地球物理勘探, 2015, 50(4): 633-642.
LI Q Y, HUANG J P, LI Z C. Undulating surface body-fitted grid seismic modeling based on fully staggered-grid mimetic finite difference[J]. Oil Geophysical Prospecting, 2015, 50(4): 633-642.
[25]
HAWKINS K. The challenge presented by North Sea Central Graden salt domes to all DMO algorithms[J]. First Break, 1994, 12(6): 327-343.
[26]
FARMER P A, JONES I F, ZHOU H, et al. Application of reverse time migration to complex imaging problems[J]. First Break, 2006, 24(9): 65-73.
[27]
JIN S, XU S, WALRAVEN D. One-return wave equation migration: Imaging of duplex waves[J]. Expanded Abstracts of 76th Annual Internat SEG Mtg, 2006, 2338-2342.
[28]
MALCOLM A E, URSIN B, MAARTEN V. Seismic imaging and illumination with internal multiples[J]. Geophysical Journal International, 2009, 176(3): 847-864. DOI:10.1111/j.1365-246X.2008.03992.x
[29]
MALCOLM A E, DEHOOP M V, URSIN B. Recursive imaging with multiply scattered waves using partial image regularization: A North Sea case study[J]. Geophysics, 2011, 76(2): B33-B42. DOI:10.1190/1.3537822
[30]
QU Y M, LI Z C, HUANG J P, et al. Prismatic and full-waveform joint inversion[J]. Applied Geophysics, 2016, 13(3): 511-518. DOI:10.1007/s11770-016-0568-7
[31]
CAVALCA M, LAILLY P. Towards the tomogtaphic inversion of prismatic reflections[J]. Expanded Abstracts of 71st Annual Internat SEG Mtg, 2001, 726-729.
[32]
MARMALYEVSKYY N, ROGANOV Y, GORNYAK Z, et al. Migration of duplex waves[J]. Expanded Abstracts of 75th Annual Internat SEG Mtg, 2005, 2025-2029.
[33]
MALCOLM A E, URSIN B, MAARTEN V. Seismic imaging and illumination with internal multiples[J]. Geophysical Journal International, 2009, 176(3): 847-864. DOI:10.1111/j.1365-246X.2008.03992.x
[34]
BROTO K, LAILLY P. Towards the tomogtaphic inversion of prismatic reflections[J]. Expanded Abstracts of 71st Annual Internat SEG Mtg, 2001: SEG-2001-0726
[35]
LI Y F, AGNIHOTRI Y, DY T. Prismatic wave imaging with dual flood RTM[J]. Expanded Abstracts of 81st Annual Internat SEG Mtg, 2011, 3290-3294.
[36]
DAI W. Multisource least-squares migration and prism wave reverse time migration[D]. Utah: The University of Utah, 2012
[37]
MAKSYM K, HENNING K. East-squares reverse time migration with prism waves[J]. Expanded Abstracts of 89th Annual Internat SEG Mtg, 2019, 4575-4579.
[38]
YANG J Z, LIU Y Z, ELITA Y Y, et al. Joint least-squares reverse time migration of primary and prismatic waves[J]. Geophysics, 2019, 84(1): S29-S40. DOI:10.1190/geo2017-0850.1
[39]
蒋丽丽, 孙建国. 基于Poisson方程的曲网格生成技术[J]. 世界地质, 2008, 27(3): 298-305.
JIANG L L, SUN J G. Source terms of elliptic system in grid generation[J]. Global Geology, 2008, 27(3): 298-305. DOI:10.3969/j.issn.1004-5589.2008.03.011
[40]
ROBERT H S, ALVIN K B. Seismic migration theory and practice[J]. The Journal of the Acoustical Society of America, 1986, 81(5): 1651-1659.