全文快速搜索:   高级搜索

  中国石油大学学报(自然科学版)  2019, Vol. 43 Issue (6): 48-58  DOI:10.3969/j.issn.1673-5005.2019.06.006
0

引用本文 [复制中英文]

周东红, 黄金强, 李振春, 等. 各向异性介质一阶速度-应力方程平面波最小二乘逆时偏移[J]. 中国石油大学学报(自然科学版), 2019, 43(6): 48-58. DOI: 10.3969/j.issn.1673-5005.2019.06.006.
[复制中文]
ZHOU Donghong, HUANG Jinqiang, LI Zhenchun, et al. Plane-wave least-squares reverse time migration based on the first-order velocity-stress equation in anisotropic media[J]. Journal of China University of Petroleum (Edition of Natural Science), 2019, 43(6): 48-58. DOI: 10.3969/j.issn.1673-5005.2019.06.006.
[复制英文]

基金项目

国家科技重大专项(2016ZX05024-003-011);国家重点基础研究发展计划(2014CB239006);国家自然科学基金项目(41774133);中央高校科研业务费专项(R1401005A);贵州大学引进人才科研项目(贵大人基合字(2018)27号);贵州省人才基地项目(RCJD2018-21)

作者简介

周东红(1968-),男,教授级高级工程师,硕士,研究方向为地震勘探方法。E-mail:zhoudh@cnooc.com.cn

通信作者

黄金强(1991-),男,博士,研究方向为地震波传播、成像与反演等。E-mail:1194998838@qq.com

文章历史

收稿日期:2019-03-21
各向异性介质一阶速度-应力方程平面波最小二乘逆时偏移
周东红1, 黄金强2, 李振春3, 吕丁友1     
1. 中海石油天津分公司渤海石油研究院, 天津 300452;
2. 贵州大学资源与环境工程学院, 贵州贵阳 550025;
3. 中国石油大学华东地球科学与技术学院, 山东青岛 266580
摘要: 为满足各向异性储层精细勘探的需求, 基于各向异性介质拟声波一阶速度-应力方程, 首先推导适用于反演成像的偏移成像算子及线性正演算子, 进而构建出泛函梯度及散射波场; 其次从最小二乘反演框架出发, 发展各向异性介质拟声波最小二乘逆时偏移成像方法; 为了提高反演效率, 进一步引入平面波编码技术, 并结合预条件、正则化及GPU加速等措施, 最终发展一种稳定高效的各向异性介质平面波最小二乘反演成像多级优化策略, 并进行典型模型反演成像测试。结果表明:该方法能够对逆掩断层、高陡倾构造、崎岖浅表层等复杂构造进行相对保幅成像, 成像效率较高, 且迭代过程中自动压制交叉干扰和偏移噪音, 进而提升剖面的整体成像质量。
关键词: 各向异性介质    一阶速度-应力方程    线性正演算子    平面波编码    平面波最小二乘偏移    
Plane-wave least-squares reverse time migration based on the first-order velocity-stress equation in anisotropic media
ZHOU Donghong1 , HUANG Jinqiang2 , LI Zhenchun3 , LÜ Dingyou1     
1. Bohai Petroleum Research Institute, Tianjin Branch, CNOOC, Tianjin 300452, China;
2. School of Resources and Environment Engineering, Guizhou University, Guiyang 550025, China;
3. School of Geosciences in China University of PetroleumEast China, Qingdao 266580, China
Abstract: High-precision imaging of anisotropic media has always been a research hotspot but remains a difficulty in the field of seismic exploration, in which the computational efficiency is one of the key factors restricting the application of inversion imaging method. In order to meet the demand of precise exploration for anisotropic reservoir, this paper first derives the imaging operator and linearized forward modeling operator suitable for inversion imaging based on the first-order velocity-stress pseudo-acoustic-wave equation of anisotropic media. The gradient of functional and scattered wavefields is then constructed. In the least square inversion framework, an anisotropic media pseudo-acoustic-wave least-squares reverse time migration method was developed afterwards.In order to improve the inversion efficiency, plane-wave coding technique, precondition, regularization, and GPU acceleration were introduced, and a stable and efficient plane-wave least-square inversion imaging optimization strategy for anisotropic media was developed. Tests on typical models show that this approach can perform with high imaging efficiency amplitude-preserving imaging on complex structures such as overthrust, high-steep-dip structure, rugged shallow surface layer and so on, and also can automatically suppress crosstalk and migration noise with the iteration to improve the overall imaging quality.
Keywords: anisotropic media    first-order velocity-stress equation    linearized modeling operator    plane-wave coding    plane-wave least-squares reverse time migration    

地震各向异性已经成为地震资料处理中不可忽视的重要因素[1]。中国西部碳酸盐岩地区、页岩油层以及海上探区都发现了大量的各向异性储层, 油气储量极其丰富[2]。然而, 各向异性地震处理仍面临诸多困难[3-5], 开展各向异性储层快速且精确成像方法研制尤为重要。各向异性储层高精度成像方法面临计算量大且多参数建模困难等问题[6], 其中各向异性Kirchhoff偏移被广泛应用于实际资料的早期成像处理中, 但成像精度不高, 且无法解决多波至、焦散区、阴影区以及复杂高陡倾构造等问题[7], 而各向异性高斯束偏移在一定程度上弥补了Kirchhoff偏移的不足[8]。各向异性波动方程成像大致分为3类:第一类是一阶速度-应力方程[9-10]; 第二类是耦合型拟声波方程, 通过设计不同的辅助变量获得多种类型的二阶耦合型拟声波方程及相应的伪横波干扰压制方案[11-12]; 第三类是纯qP波方程, Zhan等[13]简化qP波频散关系, 推导出一种时间-波数域纯qP波方程, Sheng等[14]利用波场梯度代替波矢量, 给出一种类似各向同性的纯qP方程, 此外, 借助低秩分解来构建纯qP波传播算子是当前较为稳定的实现方式[15-16]。综上, 当前各向异性成像研究主要集中在偏移成像范畴, 反演成像研究较少, 当地震数据采集不足或观测系统不规则、地下构造复杂时, 常规偏移仅能实现地下构造的模糊成像[17]。最小二乘逆时偏移(LSRTM)利用合成数据与观测数据之间的匹配程度来矫正常规RTM成像面临中深部成像能量不足、同相轴连续性较差、低频噪音严重、震源效应突出等问题, 同时缓解由于地震采集不规则、炮间能量不均衡及地震数据缺失所引起的RTM成像效果不佳。常规LSRTM大多基于各向同性假设, 其中最为典型的是Dai等[18]建立的时间-空间域LSRTM成像算法, 成像框架基本能够保证计算稳定, 在实施平面波编码、预条件、正则化、及GPU加速等多级优化措施方面趋于成熟[19]。但是, 由于没有考虑地层各向异性引起的地震波走时畸变和振幅误差, 合成数据与观测数据极不匹配, 在成像中经常出现反射波不能正确归位、同相轴发生扭曲以及迭代寻优速度慢或不收敛等问题。李振春等[20]基于二阶耦合型波动方程, 首次实现了VTI介质最小二乘逆时偏移, 并对其进行平面波加速优化, 但二阶拟声波方程存在伪横波干扰、模型参数限制、数值不稳定、及边界反射严重等问题。基于一阶速度-应力方程开展LSRTM能够更好地适应交错网格高阶有限差分, 可在不增加计算量的前提下提升计算精度, 在边界处理方面, 一阶波动方程的PML吸收效果明显优于二阶PML[21]。考虑到LSRTM计算效率低, I/O及存储成本高等问题, 笔者在实现qP波一阶速度-应力方程LSRTM算法的基础上, 考察平面波编码和GPU加速措施来提升LSRTM的运算效率, 降低计算成本。

1 方法原理 1.1 各向异性拟声波一阶速度-应力方程及线性化

一阶速度-应力方程具有不存在模型参数限制、适应变密度和交错网格高价有限差分、利于构建PML完全匹配层吸收边界条件等优势, 因此被广泛应用于各向异性介质成像与反演中[22-23]。基于Alkhalifah声学近似理论[24]的非均匀各向异性拟声波一阶速度-应力方程可描述为

$ \left\{ \begin{array}{l} s_{{\rm{p}}z}^2\frac{{\partial {\sigma _{\rm{H}}}}}{{\partial t}} = \rho \left( {1 + 2\varepsilon } \right)\left( {\frac{{\partial {v_x}}}{{\partial x}} + \frac{{\partial {v_y}}}{{\partial y}}} \right) + \rho \sqrt {1 + 2\delta } \frac{{\partial {v_z}}}{{\partial z}} + {f_x},\\ s_{{\rm{p}}z}^2\frac{{\partial {\sigma _{\rm{V}}}}}{{\partial t}} = \rho \sqrt {1 + 2\delta } \left( {\frac{{\partial {v_x}}}{{\partial x}} + \frac{{\partial {v_y}}}{{\partial y}}} \right) + \rho \frac{{\partial {v_z}}}{{\partial z}} + {f_z},\\ \rho \frac{{\partial {v_x}}}{{\partial t}} = \frac{{\partial {\sigma _{\rm{H}}}}}{{\partial x}},\\ \rho \frac{{\partial {v_y}}}{{\partial t}} = \frac{{\partial {\sigma _{\rm{H}}}}}{{\partial y}},\\ \rho \frac{{\partial {v_z}}}{{\partial t}} = \frac{{\partial {\sigma _{\rm{V}}}}}{{\partial z}}. \end{array} \right. $ (1)

式中, σHσV分别为质点振动的横向及垂向应力分量; vxvyvz分别为质点振动的速度分量; (x, y, z)为空间位置坐标; t为时间坐标; spz为沿着对称轴方向的qP波相慢度, 即纵波速度vp的倒数; εδ为各向异性参数, 用于表征介质的等效各向异性特征; ρ为介质密度; fxfy为体力分量。为了方便数学推导及理论阐述, 将正演方程(1)表达为矩阵形式:

$ {\mathit{\boldsymbol{A}}_0}\left( \mathit{\boldsymbol{M}} \right) \cdot \mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{M}} \right) = {\mathit{\boldsymbol{F}}_0}. $ (2)

其中

$ \begin{array}{l} {\mathit{\boldsymbol{A}}_0}\left( \mathit{\boldsymbol{M}} \right) = \\ \left( {\begin{array}{*{20}{c}} {s_{{\rm{p}}z}^2\frac{\partial }{{\partial t}}}&0&{ - \rho \left( {1 + 2\varepsilon } \right)\frac{\partial }{{\partial x}}}&{ - \rho \left( {1 + 2\varepsilon } \right)\frac{\partial }{{\partial y}}}&{ - \rho \sqrt {1 + 2\delta } \frac{\partial }{{\partial z}}}\\ 0&{s_{{\rm{p}}z}^2\frac{\partial }{{\partial t}}}&{ - \rho \sqrt {1 + 2\delta } \frac{\partial }{{\partial x}}}&{ - \rho \sqrt {1 + 2\delta } \frac{\partial }{{\partial y}}}&{ - \rho \frac{\partial }{{\partial z}}}\\ { - \frac{\partial }{{\partial x}}}&0&{\rho \frac{\partial }{{\partial t}}}&0&0\\ { - \frac{\partial }{{\partial y}}}&0&0&{\rho \frac{\partial }{{\partial t}}}&0\\ 0&{ - \frac{\partial }{{\partial z}}}&0&0&{\rho \frac{\partial }{{\partial t}}} \end{array}} \right), \end{array} $
$ \mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{M}} \right) = {\left( {\begin{array}{*{20}{c}} {{\sigma _{\rm{H}}}}&{{\sigma _{\rm{V}}}}&{{v_x}}&{{v_y}}&{{v_z}} \end{array}} \right)^{\rm{T}}}, $
$ {\mathit{\boldsymbol{F}}_0} = {\left( {\begin{array}{*{20}{c}} {{f_x}}&{{f_z}}&0&0&0 \end{array}} \right)^{\rm{T}}}. $

式中, A0为正演算子, 由各向异性模型参数M决定; u为正演波场向量矩阵, 与模型参数、模型尺寸、以及网格大小等计算参数有关, 并构成一一对应关系; F0为震源函数向量矩阵, 也即应力边界条件; c11c13c33分别为各向异性刚度系数。c11c13与和33与Thomsen弱各向异性参数vpεδ满足如下关系:

$ \left\{ {\begin{array}{*{20}{l}} {{c_{11}} = \rho v_{\rm{p}}^2\left( {1 + 2\varepsilon } \right) = \rho v_{\rm{p}}^2\left( {1 + 2\eta } \right)\left( {1 + 2\delta } \right),}\\ {{c_{13}} = \rho v_{\rm{p}}^2\sqrt {1 + 2\delta } ,}\\ {{c_{33}} = \rho v_{\rm{p}}^2.} \end{array}} \right. $ (3)

其中η=(ε-δ)/(1+2δ)是描述各向异性介质波前面形态的重要参数, 当η=0时, 表征椭圆(或椭球)各向异性特征; 当η>0时, qP波波前面展现出类椭圆(或类椭球)各向异性规律; 当η < 0时, qP波波前面呈现出矩形(或立体)各向异性形态, 对于各向异性二阶波动方程不再适应于该类介质[15]

最小二乘逆时偏移的关键是利用反偏移记录与观测记录的数据残差来求取反射率模型残差, 即反演成像结果的更新量, 这里假定了反射率模型残差与数据残差之间符合线性关系。所谓反偏移, 即通过成像结果来求取刻画反射波信息的反偏移记录, 将其与观测记录进行比较, 用以判定该成像结果是否达到最优, 为此, 反偏移算子或者说线性正演算子的构建至关重要, 定义表征反射率模型的速度扰动$ m = \Delta s_{{\rm{p}}z}^2, \left( {s_{{\rm{p}}z}^2 = s_{{\rm{p}}z0}^2 + \Delta s_{{\rm{p}}z}^2} \right)$, 根据地震波摄动理论, 背景速度与背景波场存在线性关系, 应用Born近似, 即采用背景波场代替总场, 并忽略等号右端的高阶项$ \Delta s_{{\rm{p}}z0}^2\frac{{\partial \delta {\sigma _{\rm{H}}}}}{{\partial t}}, \Delta s_{{\rm{p}}z0}^2\frac{{\partial \delta {\sigma _{\rm{v}}}}}{{\partial t}}$, 可导出线性化波动方程:

$ \left\{ \begin{array}{l} s_{{\rm{p}}z{\rm{0}}}^2\frac{{\partial \delta {\sigma _{\rm{H}}}}}{{\partial t}} = \rho \left( {1 + 2\varepsilon } \right)\left( {\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{\partial \delta {v_y}}}{{\partial y}}} \right) + \\ \rho \sqrt {1 + 2\delta } \frac{{\partial \delta {v_z}}}{{\partial z}} - m\frac{{\partial {\sigma _{\rm{H}}}}}{{\partial t}},\\ s_{{\rm{p}}z{\rm{0}}}^2\frac{{\partial \delta {\sigma _{\rm{V}}}}}{{\partial t}} = \rho \sqrt {1 + 2\delta } \left( {\frac{{\partial \delta {v_x}}}{{\partial x}} + \frac{{\partial \delta {v_y}}}{{\partial y}}} \right) + \\ \rho \frac{{\partial \delta {v_z}}}{{\partial z}} - m\frac{{\partial {\sigma _{\rm{V}}}}}{{\partial t}},\\ \rho \frac{{\partial \delta {v_x}}}{{\partial t}} = \frac{{\partial \delta {\sigma _{\rm{H}}}}}{{\partial x}},\\ \rho \frac{{\partial \delta {v_y}}}{{\partial t}} = \frac{{\partial \delta {\sigma _{\rm{H}}}}}{{\partial y}},\\ \rho \frac{{\partial \delta {v_z}}}{{\partial t}} = \frac{{\partial \delta {\sigma _{\rm{V}}}}}{{\partial z}}. \end{array} \right. $ (4)

用矩阵表示为

$ \mathit{\boldsymbol{A}}\left( \mathit{\boldsymbol{M}} \right) \cdot \delta \mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{M}} \right) = m\mathit{\boldsymbol{F}}. $ (5)

其中

$ \delta \mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{M}} \right) = {\left( {\begin{array}{*{20}{c}} {\delta {\sigma _{\rm{H}}}}&{\delta {\sigma _{\rm{V}}}}&{\delta {v_x}}&{\delta {v_y}}&{\delta {v_z}} \end{array}} \right)^{\rm{T}}}, $
$ \mathit{\boldsymbol{F}} = {\left( {\begin{array}{*{20}{c}} {\frac{{\delta {\sigma _{\rm{H}}}}}{{\partial t}}}&{\frac{{\delta {\sigma _{\rm{V}}}}}{{\partial t}}}&0&0&0 \end{array}} \right)^{\rm{T}}}. $

式中, δu(M)为扰动波场矩阵, 主要用于刻画地震反射波信息; F为由入射波场构建的震源函数向量矩阵; 线性正演算子A的形式与A0相同, 但A是由各向异性背景参数构建。

1.2 伴随方程与梯度更新公式

根据Tarantola的梯度构建方案[25], 需要对数据残差进行反向延拓, 进而计算伴随波场, 为了获取伴随波场, 假定χHχV分别表示横向及垂向的伴随应力分量, Δdx、Δdz分别表示观测记录残差分量, 可由伴随状态法导出如下伴随方程:

$ \left\{ \begin{array}{l} s_{{\rm{p}}z0}^2\frac{{\partial {\chi _{\rm{H}}}}}{{\partial t}} = \frac{{\partial {u_x}}}{{\partial x}} + \frac{{\partial {u_y}}}{{\partial y}} - \Delta {d_x},\\ s_{{\rm{p}}z0}^2\frac{{\partial {\chi _{\rm{V}}}}}{{\partial t}} = \frac{{\partial {u_z}}}{{\partial z}} - \Delta {d_z},\\ \rho \frac{{\partial {u_x}}}{{\partial t}} = \rho \left( {1 + 2\varepsilon } \right)\frac{{\partial {\chi _{\rm{H}}}}}{{\partial x}} + \rho \sqrt {1 + 2\delta } \frac{{\partial {\chi _{\rm{V}}}}}{{\partial x}},\\ \rho \frac{{\partial {u_y}}}{{\partial t}} = \rho \left( {1 + 2\varepsilon } \right)\frac{{\partial {\chi _{\rm{H}}}}}{{\partial y}} + \rho \sqrt {1 + 2\delta } \frac{{\partial {\chi _{\rm{V}}}}}{{\partial y}},\\ \rho \frac{{\partial {u_z}}}{{\partial t}} = \rho \sqrt {1 + 2\delta } \frac{{\partial {\chi _{\rm{H}}}}}{{\partial z}} + \frac{{\partial {\chi _{\rm{V}}}}}{{\partial z}}. \end{array} \right. $ (6)

为表述统一, 将式(6)写为矩阵形式:

$ {\mathit{\boldsymbol{A}}^*}\left( \mathit{\boldsymbol{M}} \right) \cdot {\mathit{\boldsymbol{u}}^*}\left( \mathit{\boldsymbol{M}} \right) = \Delta \mathit{\boldsymbol{d}} = \delta \mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{M}} \right) - {\mathit{\boldsymbol{d}}^{{\rm{obs}}}}. $ (7)

其中

$ \begin{array}{l} {\mathit{\boldsymbol{A}}^*}\left( \mathit{\boldsymbol{M}} \right) = \\ \left( {\begin{array}{*{20}{c}} { - s_{{\rm{p}}z0}^2\frac{\partial }{{\partial t}}}&0&{\frac{\partial }{{\partial x}}}&{\frac{\partial }{{\partial y}}}&0\\ 0&{ - s_{{\rm{p}}z0}^2\frac{\partial }{{\partial t}}}&0&0&{\frac{\partial }{{\partial z}}}\\ {\frac{\partial }{{\partial x}}\rho \left( {1 + 2\varepsilon } \right) \cdot }&{\frac{\partial }{{\partial x}}\rho \sqrt {1 + 2\delta } \cdot }&{ - \rho \frac{\partial }{{\partial t}}}&0&0\\ {\frac{\partial }{{\partial y}}\rho \left( {1 + 2\varepsilon } \right) \cdot }&{\frac{\partial }{{\partial y}}\rho \sqrt {1 + 2\delta } \cdot }&0&{ - \rho \frac{\partial }{{\partial t}}}&0\\ {\frac{\partial }{{\partial z}}\rho \sqrt {1 + 2\delta } \cdot }&{\frac{\partial }{{\partial z}}\rho \cdot }&0&0&{ - \rho \frac{\partial }{{\partial t}}} \end{array}} \right), \end{array} $
$ {\mathit{\boldsymbol{u}}^ * }\left( \mathit{\boldsymbol{M}} \right) = {\left( {\begin{array}{*{20}{l}} {{\chi _{\rm{H}}}}&{{\chi _{\rm{V}}}}&{{u_x}}&{{u_y}}&{{u_z}} \end{array}} \right)^{\rm{T}}}, $
$ \Delta \mathit{\boldsymbol{d}} = {\left( {\begin{array}{*{20}{l}} {\Delta {d_x}}&{\Delta {d_z}}&0&0&0 \end{array}} \right)^{\rm{T}}}. $

式中, A*(M)为伴随正演算子, 与线性正演算子A(M)满足共轭转置特征, 即(A*)-1=(A-1)*; u*(M)为伴随波场向量矩阵; Δd为数据残差。

VTI介质误差泛函可表示为水平和垂直两个分量的残差在L2模条件下同时达到最小, 即目标泛函为

$ \mathit{\boldsymbol{J}}\left( \mathit{\boldsymbol{M}} \right) = \frac{1}{2}\left\| {\delta \mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{M}} \right) - {\mathit{\boldsymbol{d}}^{{\rm{obs}}}}} \right\|_2^2. $ (8)

与全波形反演类似, 但不同之处在于, 最小二乘偏移的目标泛函是反偏移记录与观测记录的残差, 而全波形反演的目标泛函是正演记录与观测记录的残差, 然后根据线性正演公式(5), 可得

$ \frac{{\partial \delta \mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{M}} \right)}}{{\partial m}} = {\mathit{\boldsymbol{A}}^{ - 1}}\left( \mathit{\boldsymbol{M}} \right)\frac{{\partial \left( {m\mathit{\boldsymbol{F}}} \right)}}{{\partial m}} = {\mathit{\boldsymbol{A}}^{ - 1}}\left( \mathit{\boldsymbol{M}} \right)\mathit{\boldsymbol{F}}. $ (9)

因此, 目标函数J(m)对本文中定义的速度模型扰动m的梯度可表示为

$ \begin{array}{l} \frac{{\partial \mathit{\boldsymbol{J}}\left( \mathit{\boldsymbol{M}} \right)}}{{\partial m}} = \left\langle {\frac{{\partial \delta \mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{M}} \right)}}{{\partial m}},\delta \mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{M}} \right) - {\mathit{\boldsymbol{d}}^{{\rm{obs}}}}} \right\rangle = \\ \left\langle {{\mathit{\boldsymbol{A}}^{ - 1}}\left( \mathit{\boldsymbol{M}} \right)\mathit{\boldsymbol{F}},\delta \mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{M}} \right) - {\mathit{\boldsymbol{d}}^{{\rm{obs}}}}} \right\rangle = \\ \left\langle {\mathit{\boldsymbol{F}},{{\left[ {{\mathit{\boldsymbol{A}}^{ - 1}}\left( \mathit{\boldsymbol{M}} \right)} \right]}^*}\left( {\delta \mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{M}} \right) - {\mathit{\boldsymbol{d}}^{{\rm{obs}}}}} \right)} \right\rangle = \\ \left\langle {\mathit{\boldsymbol{F}},{{\left[ {{\mathit{\boldsymbol{A}}^*}\left( \mathit{\boldsymbol{M}} \right)} \right]}^{ - 1}}\left( {\delta \mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{M}} \right) - {\mathit{\boldsymbol{d}}^{{\rm{obs}}}}} \right)} \right\rangle = \left\langle {\mathit{\boldsymbol{F}},{\mathit{\boldsymbol{u}}^*}\left( \mathit{\boldsymbol{M}} \right)} \right\rangle = \\ \int_0^{\rm{T}} {\left( {\frac{{\partial {\sigma _{\rm{H}}}}}{{\partial t}}{\chi _{\rm{H}}} + \frac{{\partial {\sigma _{\rm{V}}}}}{{\partial t}}{\chi _{\rm{V}}}} \right){\rm{d}}t} = {\sigma _{\rm{H}}}{\chi _{\rm{H}}} + {\sigma _{\rm{V}}}{\chi _{\rm{V}}}. \end{array} $ (10)

由此可知, 伴随波场与正传波场的零延迟互相关即为目标泛函梯度, 在求取数据残差偏移结果时, 本文中统一采用震源补偿作为正则化算子来改善收敛速度, 并增强中深部成像能量。最小二乘偏移与FWI的不同之处在于线性化正演算子A(M)与背景模型参数M有关, 但与模型参数扰动m无关。

综上, 各向异性介质平面波LSRTM算法流程可概述为:①通过线性正演算子计算预测波场, 进而求取炮数据残差; ②利用偏移算子及梯度公式计算梯度方向, 成像条件决定于梯度公式; ③通过预条件算子修正梯度方向和步长得到共轭梯度方向和新的步长; ④计算模型更新量, 更新成像结果; ⑤将数据残差或模型残差与预先设定的阈值进行对比, 判断是否需要进行下一次迭代。

1.3 Collino分裂式PML边界条件

考虑到PML边界条件, 在交错网格剖分方式下的频空域波动方程迭代更新公式为

$ \left\{ \begin{array}{l} \left( {{\rm{i}}\omega + {d_x}} \right){{\tilde v}_x} = \frac{1}{\rho }\frac{{\partial \tilde p}}{{\partial x}},\\ \left( {{\rm{i}}\omega + {d_z}} \right){{\tilde v}_z} = \frac{1}{\rho }\frac{{\partial \tilde q}}{{\partial x}},\\ \tilde p = {{\tilde p}^\parallel } + {{\tilde p}^ \bot },\\ \left( {{\rm{i}}\omega + {d_x}} \right){{\tilde p}^\parallel } = {c_{11}}\frac{{\partial {{\tilde v}_x}}}{{\partial x}},\\ \left( {{\rm{i}}\omega + {d_z}} \right){{\tilde p}^ \bot } = {c_{13}}\frac{{\partial {{\tilde v}_z}}}{{\partial z}},\\ \tilde q = {{\tilde q}^\parallel } + {{\tilde q}^ \bot },\\ \left( {{\rm{i}}\omega + {d_x}} \right){{\tilde q}^\parallel } = {c_{13}}\frac{{\partial {{\tilde v}_x}}}{{\partial x}},\\ \left( {{\rm{i}}\omega + {d_z}} \right){{\tilde q}^ \bot } = {c_{33}}\frac{{\partial {{\tilde v}_z}}}{{\partial z}}. \end{array} \right. $ (11)

式中, ω为角频率; $ {{\tilde v}_x}$$ {{\tilde v}_z}$为频率域的速度分量; $ {\tilde p}$$ {\tilde q}$为频率域的应力分量; 上标||为分裂后水平方向的波场分量, 上标⊥为分裂后垂直方向的波场分量; dxdz分别为水平和垂直方向的衰减系数。将式(11)变换到时间域, 质点横向振动速度${{\tilde v}_x} $与横向应力$ {\tilde p}$的离散更新公式可表示为

$ \left\{ \begin{array}{l} {D_t}\left( {{v_x}} \right)_{i + \frac{1}{2},j}^n + {d_{x,i + \frac{1}{2}}}\frac{{\left( {{v_x}} \right)_{i + \frac{1}{2},j}^{n + \frac{1}{2}} + \left( {{v_x}} \right)_{i + \frac{1}{2},j}^{n - \frac{1}{2}}}}{2} = \\ \frac{1}{\rho }{D_x}\left( p \right)_{i + \frac{1}{2},j}^n\left( p \right)_{i,j}^{n + 1} = \left( {{p^\parallel }} \right)_{i,j}^{n + 1} + \left( {{p^ \bot }} \right)_{i,j}^{n + 1},\\ {D_t}\left( {{p^\parallel }} \right)_{i,j}^{n + \frac{1}{2}} + {d_{x,i}}\frac{{\left( {{p^\parallel }} \right)_{i,j}^{n + 1} + \left( {{p^\parallel }} \right)_{i,j}^n}}{2} = {c_{11}}{D_x}\left( {{v_x}} \right)_{i,j}^{n + \frac{1}{2}},\\ {D_t}\left( {{p^ \bot }} \right)_{i,j}^{n + \frac{1}{2}} + {d_{x,i}}\frac{{\left( {{p^ \bot }} \right)_{i,j}^{n + 1} + \left( {{p^\parallel }} \right)_{i,j}^n}}{2} = {c_{11}}{D_x}\left( {{v_x}} \right)_{i,j}^{n + \frac{1}{2}}. \end{array} \right. $ (12)

式中, 下标(i, j)为空间整网格点, $ \left( {i + \frac{1}{2}, j} \right)$$ \left( {i, j + \frac{1}{2}} \right)$$ \left( {i + \frac{1}{2}, j + \frac{1}{2}} \right)$为空间半网格点; 上标n为时间整网格点, $ n + \frac{1}{2}$为时间半网格点。在时空交错网格中, 时间偏倒数Dt采用中心差分、空间偏导数DxDz采用高阶有限差分, 进一步得到振动速度$ {{\tilde v}_x}$与横向应力$ {\tilde p}$的详细迭代格式:

$ \left\{ \begin{array}{l} \frac{{\left( {{v_x}} \right)_{i + \frac{1}{2},j}^{n + \frac{1}{2}} - \left( {{v_x}} \right)_{i + \frac{1}{2},j}^{n - \frac{1}{2}}}}{{dt}} + {d_{x,i + \frac{1}{2}}}\frac{{\left( {{v_x}} \right)_{i + \frac{1}{2},j}^{n + \frac{1}{2}} + \left( {{v_x}} \right)_{i + \frac{1}{2},j}^{n - \frac{1}{2}}}}{2} = \\ \frac{1}{\rho }{D_x}\left( p \right)_{i + \frac{1}{2},j}^n,\\ \left( {\frac{1}{{dt}} + \frac{{{d_{x,i + \frac{1}{2}}}}}{2}} \right)\left( {{v_x}} \right)_{i + \frac{1}{2},j}^{n + \frac{1}{2}} = \left( {\frac{1}{{dt}} - \frac{{{d_{x,i + \frac{1}{2}}}}}{2}} \right)\left( {{v_x}} \right)_{i + \frac{1}{2},j}^{n - \frac{1}{2}} = \\ \frac{1}{\rho }{D_x}\left( p \right)_{i + \frac{1}{2},j}^n,\\ \frac{{\left( {{p^\parallel }} \right)_{i,j}^{n + 1} - \left( {{p^\parallel }} \right)_{i,j}^n}}{{dt}} + {d_{x,i}}\frac{{\left( {{p^\parallel }} \right)_{i,j}^{n + 1} + \left( {{p^\parallel }} \right)_{i,j}^n}}{2} = \\ {c_{11}}{D_x}\left( {{v_x}} \right)_{i,j}^{n + \frac{1}{2}}\left( {\frac{1}{{dt}} + \frac{{{d_{x,i}}}}{2}} \right)\left( {{p^\parallel }} \right)_{i,j}^{n + 1} = \left( {\frac{1}{{dt}} - \frac{{{d_{x,i}}}}{2}} \right)\left( {{p^\parallel }} \right)_{i,j}^n + \\ {c_{11}}{D_x}\left( {{v_x}} \right)_{i,j}^{n + \frac{1}{2}},\\ \frac{{\left( {{p^ \bot }} \right)_{i,j}^{n + 1} - \left( {{p^ \bot }} \right)_{i,j}^n}}{{dt}} + {d_{z,j}}\frac{{\left( {{p^ \bot }} \right)_{i,j}^{n + 1} + \left( {{p^ \bot }} \right)_{i,j}^n}}{2} = {c_{13}}{D_z}\left( {{v_z}} \right)_{i,j}^{n + \frac{1}{2}},\\ \left( {\frac{1}{{dt}} + \frac{{{d_{z,j}}}}{2}} \right)\left( {{p^ \bot }} \right)_{i,j}^{n + 1} = \left( {\frac{1}{{dt}} - \frac{{{d_{z,j}}}}{2}} \right)\left( {{p^ \bot }} \right)_{i,j}^n + {c_{11}}{D_z}\left( {{v_z}} \right)_{i,j}^{n + \frac{1}{2}}. \end{array} \right. $ (13)

其中衰减参数为

$ {d_{x,i}} = d0xd_x^2, $
$ {d_{z,i}} = d0zd_z^2, $
$ d0x = 3.0{v_{{\rm{pmax}}}}\log \left( {\frac{1}{{{R_{{\rm{coef}}}}}}} \right)/\left( {2{N_{{\rm{PML}}}}{\rm{d}}x} \right), $
$ d0z = 3.0{v_{{\rm{pmax}}}}\log \left( {\frac{1}{{{R_{{\rm{coef}}}}}}} \right)/\left( {2{N_{{\rm{PML}}}}{\rm{d}}z} \right), $
$ {d_x} = \frac{{\left| {x - {x_{{\rm{PML}}}}} \right|}}{{{N_{{\rm{PML}}}}{\rm{d}}x}},{d_z} = \frac{{\left| {z - {z_{{\rm{PML}}}}} \right|}}{{{N_{{\rm{PML}}}}{\rm{d}}z}}. $

式中, vpmax为模型的最大纵波速度; 在各向异性介质中, Rcoef为经验系数, 一般取0.000001;NPML为PML扩展层的网格厚度; dxdz分别为纵横向网格间距; dxdz分别由衰减层位置(x, z)到对应衰减边界(xPML, zPML)的距离决定。分析迭代公式可见, 该PML离散格式能够很好地适应GPU异构平台下的算法编写, 下文中统一采用PML边界条件, 扩展层厚度为30个网格点。

1.4 平面波编码技术

平面波编码是多震源混叠数据直接成像技术的一种实现方式, 目前在多震源全波形反演、多震源最小二乘偏移中被广泛使用[20]。相对于常规的逐炮偏移, 平面波偏移是将几十或几百炮地震数据按照线性时移编码合成几个或几十个平面波域数据, 然后进行多炮同时偏移, 进而提高计算效率的一种加速手段, 在x-t域, 可用射线参数p来表示经过不同时移的平面波记录:

$ p = \frac{{{\rm{d}}t}}{{{\rm{d}}x}} = \frac{{\sin \theta }}{{v\left( \theta \right)}}. $ (14)

式中, θ为地表激发倾角; v为在地表位置θ方向的纵波相速度。对于二维观测方式, 平面波记录可表示为

$ \mathit{\boldsymbol{d}}\left( {{x_g},t;p} \right) = \sum\limits_{{x_s}} \mathit{\boldsymbol{d}} \left( {{x_{\rm{g}}},t;{x_{\rm{s}}}} \right) * \mathit{\boldsymbol{\delta }}\left( {t - p \cdot {x_{\rm{s}}}} \right). $ (15)

式中, d(xg, t; xs)为炮域的地震数据; d(xg, t; p)为平面波域地震数据, 即平面波记录; δ(t-p·xs)为时移函数, 时移量为p·xs, 随震源坐标线性变化; *为卷积运算; xs为每一炮震源坐标; xg为接收点坐标。因此平面波编码过程可以表述为:首先将炮域的地震数据进行逐炮时移, 然后进行线性叠加, 最终将多炮地震数据压缩成1个平面波记录。

平面波偏移是使用相同编码方式的平面波震源激发产生震源波场, 将平面波记录作为边界进行反向传播获得检波点波场, 互相关成像即获得平面波成像结果, 可见其计算量仅相当于1炮地震数据的偏移成本。然而单个平面波成像结果通常会遭受严重的串扰噪音, 为了在提高计算效率的同时保证计算精度, 通常将多个不同编码方式的平面波成像结果进行叠加, 以此来增强有效信号, 减弱串扰噪音。平面波域最终成像结果可以表示为

$ {I_{\rm{m}}}\left( x \right) = \sum\limits_{i = 1}^N {{I_{{\rm{m}}i}}\left( {x,{p_i}} \right)} . $ (16)

式中, pi为合成第i个平面波记录所选用的射线参数; Imi为第i个平面波记录作为输入获得的成像结果; Im为最终成像结果; N为总的平面波记录个数, 通常远小于炮数。

2 模型测试

本文中采用简单Salt模型检验方法的有效性, 并在逆掩断层模型上验证方法对复杂构造模型的适应性。

2.1 Salt模型正演

图 1(a)~(d)分别为Salt模型对应的真实纵波速度、偏移纵波速度以及各向异性参数。真实纵波速度用来合成观测记录, 偏移纵波速度是通过对真实纵波速度的平滑获得, 用以表征背景速度。根据定义, 通过背景速度与真实速度能够求出反射系数模型(图 1(e)), 这是反演成像试图通过迭代来逼近的理论值, 也是成像试验中用以比较成像质量优劣的重要参考。需要说明的是, 正演计算中对各向异性参数进行平滑处理可以提高算法的稳定性。模型网格大小为676×201, 空间网格间距Δxz=10m, 因此模型尺寸总计为6750 m×2000 m, 模型速度变化范围为1500~4482 m/s, 权衡稳定性条件与计算效率, 时间步长Δt设计为0.8 ms, 选定主频为20 Hz的雷克子波作为激发震源。

图 1 简单Salt模型 Fig.1 A simple Salt model

观测系统设计如下:激发总炮数338炮, 激发井深20 m, 炮点间距20 m, 首炮位置位于模型左端第一个网格点处, 采用676道接收, 接收点间距10 m, 总的记录长度为5.0 s, 因此每道采样点数为6251。首先按照上述参数选用真实纵波速度模型及各向异性参数, 并采用高阶有限差分线性正演算法(式(4))合成VTI介质散射波地震响应(图 2(a)), 分析可见:①线性化的反偏移记录, 也即散射波记录中不含有直达波, 这是由于直达波属于入射波场; ②散射波记录中蕴含着更加丰富的反射波信息, 这是因为散射波的激发震源为地下介质中所有的反射源, 且无直达波等强能量干扰; ③反射同相轴清晰、不存在严重的数值频散假象、反射波刻画准确, 这些现象初步证实了本文中线性化正演算子的正确性。

图 2 散射波地震响应 Fig.2 Seismic response of linearized scattered-wave

为了测试平面波偏移成像方法的优势。将上述Salt模型的炮域地震记录合成31个含不同射线参数p的平面波域地震记录, p的取值范围是-0.27~0.27 s·km-1, 图 2(b)p=0时的平面波记录, 观察分析可知:①由于采用了平面波线性时移编码策略, 将338炮炮域地震记录压缩成31个平面波记录, 数据量减少为原来的1/10, 释放了大量存储空间, 计算中也明显降低了I/O成本与计算时间; ②合成的平面波记录在地震波形形态上更接近于叠加剖面, 从平面波记录中可以粗略地判断出浅表层的地层分布情况; ③深部地层及复杂地质体所对应的反射波由于相互干涉叠加, 很难从记录中直接辨识, 需要通过进一步的偏移成像处理来使反射波归位、绕射波收敛。

2.2 偏移成像试验

图 3为本文中Salt模型的最小二乘反演成像结果, 其中图 3(a)(b)分别为VTI介质一阶速度-应力方程LSRTM迭代1次和30次的成像剖面。由图 3中可见:在LSRTM迭代1次, 也即RTM的成像结果中(图 3(a)), 出现明显的震源效应和低频噪音, 中深部成像能量较弱, 在覆盖次数不足的区域出现由浅至深呈“斜坡状”分布的弱照明现象, 在盐体边界位置出现成像不清晰, 在盐体下部及顶部构造处出现严重的偏移噪音; LSRTM迭代30次后(图 3(b)), 震源效应得到很好地抑制, 低频噪音基本消除, 中深部同相轴的振幅能量更加均衡, 盐体边界的成像分辨率更高, 成像剖面的信噪比整体上得到大幅度提升。

图 3 反演成像结果 Fig.3 Inversion images

图 3(c)(d)分别为VTI介质一阶速度-应力方程PLSRTM迭代1次和30次的成像剖面。对比可知:在PLSRTM迭代1次, 也即平面波RTM的成像结果中, 由于采用了平面波震源, 没有出现点震源效应, 对于低频噪音严重等现象与图 3(a)几乎一致; 在理论上, 平面波RTM结果可能会存在交叉干扰, 而图 3(c)中没有出现明显的交叉干扰, 这是因为对31个平面波道集的成像剖面进行偏移叠加已经显著地压制掉了交叉干扰, 提升了剖面信噪比; 通过PLSRTM迭代, 图 3(d)的成像质量接近于图 3(b), 然而计算效率较LSRTM明显提升, 该模型的计算时间约为LSRTM的1/10, 这是因为借助平面波编码将338炮域地震记录压缩成了31个平面波记录, 致使成像过程中炮循环次数由338降为31, 而成像算法的其他模块基本保持不变, 提升运行效率是PLSRTM的主要优势, 同时也降低存储成本。

图 3(e)(f)分别为各向同性PLSRTM迭代1次和30次的成像结果, 与图 3(c)(d)对比可见:在1次迭代结果中, 盐丘体上部边界出现成像模糊, 盐丘体下覆介质中噪音严重, 成像杂乱, 总体成像效果变差; 迭代30次之后, 低频噪音得到压制, 深部能量得到增强, 但是断层面等陡倾构造成像不清晰(图(f)虚线箭头所指), 盐丘下部的地层界面同相轴连续性变差(图(f)虚线框所示), 整体成像质量欠佳, 这是由于输入的地震记录存在各向异性效应所致, 与各向同性相比, 地震记录走时及振幅都存在较大差异, 采用各向同性算法进行反演成像处理, 将导致模拟数据与观测数据的走时出现严重不匹配, 即周波跳跃现象, 使反演陷入局部极值; 上述现象也从侧面说明, 最小二乘迭代可以压制低频成像噪音、补偿中深部成像能量等, 但是在地震波走时存在严重畸变的情况下, 仍需发展其他方法来缓解走时不匹配的问题。

为了更直观地展示LSRTM及PLSRTM的成像特点, 图 4给出了上述3组成像试验的数据残差收敛曲线。从收敛曲线变化趋势反映出:①随着迭代次数的增加, 各向异性PLSRTM与各向异性LSRTM的数据残差都能稳定地收敛到5%以下, 表明各向异性PLSRTM与各向异性LSRTM的成像质量相当, 但是各向异性PLSRTM具有更高的计算效率; ②在前1~8次迭代中, 各向异性PLSRTM的收敛效果略差于各向异性LSRTM, 这是由于使用平面波编码将炮域地震数据压缩为平面波记录, 在反演成像中产生了交叉干扰, 降低了成像剖面的信噪比, 映射到地震数据中即表现为较大的数据残差; ③迭代8次之后, 各向异性PLSRTM的收敛效果与各向异性LSRTM几乎完全一致, 这说明各向异性PLSRTM方法能够在迭代过程中自动压制串扰噪音, 提高剖面的信噪比, 改善成像质量; ④各向同性PLSRTM的数据残差只能收敛到60%, 在迭代10次之后, 收敛极其缓慢, 基本不再下降, 这反映出由于数据走时不匹配, 出现了严重的周波跳跃现象, 使反演过程陷入局部极值, 不能获得全局最优解。上述试验现象清楚地证实了本文中反演成像算法的正确性及开展各向异性研究的必要性。

图 4 归一化数据残差收敛曲线 Fig.4 Normalized convergence curves of data residuals
2.3 复杂逆掩断层模型测试

在验证PLSRTM成像方法正确性的基础上, 将该方法应用于复杂模型成像试算, 进一步验证方法对复杂模型的适应性和算法的健壮性。改造的复杂加拿大逆掩断层模型参数如图 5所示, 从真实速度模型(图 5(a))可见, 浅表层横向变速非常剧烈, 地下地层中背斜与向斜构造、高陡构造, 以及逆掩断层发育, 是检验成像算法适应性和健壮性的典型模型。图 5(b)为偏移速度模型, 用于构建反演成像算法, 图 5(c)(d)分别为各向异性Epsilon和Delta模型参数, 是表征速度各向异性的重要参数。

图 5 改造的逆掩断层模型 Fig.5 Modified Canada overthrust model

模型网格尺寸分别为834×500, 纵横向网格间距都是20 m, 因此该模型的横向长度为16660 m, 纵向深度为9980 m, 浅层速度为3599.95 m/s, 模型最大速度为6000 m/s, 选取时间步长为1.0 ms, 子波主频30 Hz。观测系统设计如下:激发总炮数为217炮, 激发井深20 m, 炮点间距40 m, 首炮位置位于模型左端第一个网格点处, 采用834道全接收, 接收点间距20 m, 总的记录长度为5.0 s, 因此每道地震数据的采样点数为5001。

按照上述观测系统和计算参数, 采用VTI介质高阶有限差分法线性化正演模拟算子计算出炮域地震记录, 并将其合成35个平面波记录, 射线参数p=0时的平面波记录如图 6所示。从图中可见, 由于浅表层崎岖变化, 且受构造运动影响剧烈, 造成平面波记录中反射结构杂乱、绕射波发育, 很难直接识别出地层分布情况, 图中没有直达波, 主要的反射波信息清晰可见, 这是线性正演记录的重要标志。

图 6 VTI介质合成平面波记录 Fig.6 Plane-wave records of VTI media

对上述平面波记录进行平面波域最小二乘逆时偏移成像, 迭代1次和30次的成像结果分别如图 7(a)(b)所示。由图 7中可见:①PLSRTM迭代1次的成像结果中交叉干扰严重(如实线箭头所指)、受偏移噪音干扰、剖面信噪比低、分辨率差、复杂构造区域的同相轴连续性欠佳; ②通过PLSRTM迭代, 迭代结果中交叉干扰、偏移噪音基本消失, 信噪比显著提高; 逆掩断层、高陡倾构造、崎岖浅表层等复杂构造准确、清晰地成像; 深部同相轴的连续性更好, 剖面整体分辨率得到提升。上述实验结果验证了本文方法对高陡构造、逆掩断层、崎岖浅表层等复杂构造的成像适应性及本算法的健壮性。

图 7 VTI介质平面波最小二乘逆时偏移成像结果 Fig.7 Plane-wave LSRTM images by using PLSRTM

进一步, 图 8给出了逆掩断层模型VTI介质PLSRTM成像的数据残差收敛曲线, 该曲线更直观地反映出PLSRTM成像方法能够使数据残差迅速稳定地收敛到5%以下; 在前几次迭代中, 收敛速度最快, 随着迭代次数的增加, 数据残差的收敛速度逐渐变得缓慢。

图 8 归一化数据残差收敛曲线 Fig.8 Normalized convergence curves of data residuals

此外, 考虑到地震叠前逆时偏移方法受地震子波、地震数据的品质、背景参数的精度, 以及地球介质的理论假设等因素的综合影响, 今后需要进一步考察该算法对含噪音地震数据的适应能力, 以及对速度模型精度的依赖性等问题。

3 结论

(1) 采用推导的线性正演算子能较精确地重构qP波一阶散射波场, 运用推导的偏移算子作用于数据残差可获取较为准确的误差泛函梯度向量。

(2) 将平面波编码技术等引入到各向异性介质最小二乘偏移中, 发展了各向异性介质平面波最小二乘偏移方法, 能够极大地提高偏移成像效率, 且迭代过程中自动压制交叉干扰和偏移噪音。

(3) PLSRTM成像方法能够对逆掩断层、高陡倾构造、崎岖浅表层等复杂构造进行准确、清晰地成像, 信噪比显著提高, 剖面整体分辨率得到提升。

(4) 考虑地震各向异性, 在速度模型及各向异性参数相对准确的情况下, 可以实现地下构造的准确成像, 通过采用各向异性最小二乘逆时偏移, 能够实现地下介质的相对保幅成像。

参考文献
[1]
HOUGHTON J C. Use of the truncated shifted Pareto distribution in assessing size distribution of oil and gas fields[J]. Mathematical Geology, 1988, 20(8): 907-937. DOI:10.1007/BF00892970
[2]
POON D C, MCCORMACK M, THIMM H F. The application of fractal geostatistics to oil and gas property evaluation and reserve estimates[J]. The Journal of Canadian Petroleum Technology, 1993, 32(10): 24-27.
[3]
BARTON C C, LA POINTE P R. Fractals in the earth sciences[M]. New York: Plenum Press, 1995: 13-33.
[4]
LAHERRERE J. Distribution of field sizes in a petroleum system: parabolic fractal, lognormal or stretched exponential?[J]. Marine & Petroleum Geology, 2000, 17(4): 539-546.
[5]
曾怡. 分形法预测油气储量与资源量[J]. 石油实验地质, 1998(2): 152-154.
ZENG Yi. A new method of petroleum resources estimation: fractal method[J]. Petroleum Geology & Experiment, 1998(2): 152-154.
[6]
郭秋麟, 谢红兵, 米石云, 等. 油气资源分布的分形特征及应用[J]. 石油学报, 2009, 30(3): 379-385.
GUO Qiulin, XIE Hongbing, MI Shiyun, et al. Fractal model for petroleum resource distribution and its application[J]. Acta Petrolei Sinica, 2009, 30(3): 379-385.
[7]
刘晓冬, 徐景祯. 分形方法预测气田数量及其储量[J]. 石油学报, 2000, 21(2): 42-44.
LIU Xiaodong, XU jingzhen. The number and reserves of gas fields are predicted by fractal method[J]. Acta Petrolei Sinica, 2000, 21(2): 42-44. DOI:10.3321/j.issn:0253-2697.2000.02.008
[8]
陈新, 王绪龙, 靳涛. 石油储量分布的分形特征及其预测[J]. 新疆地质, 2001, 19(4): 297-299.
CHEN Xin, WANG Xulong, JIN Tao. The fractal characteristic and its prediction of the distribution of oil reserves[J]. Xinjiang Geology, 2001, 19(4): 297-299. DOI:10.3969/j.issn.1000-8845.2001.04.014
[9]
冯阵东, 戴俊生, 刘景东. 饶阳凹陷油藏储量在空间分布的分形特征[J]. 新疆石油地质, 2010, 31(4): 372-375.
FENG Zhendong, DAI Junsheng, LIU Jingdong. Fractal characters of oil Reserves in spatial distribution in Raoyang Sag[J]. Xinjiang Petroleum Geology, 2010, 31(4): 372-375.
[10]
宋宁, 王铁冠, 刘东鹰, 等. 分形方法在苏北盆地金湖凹陷石油资源评价中的应用[J]. 地质科学, 2006, 41(4): 578-585.
SONG Ning, WANG Tieguan, LIU Dongying, et al. Application of fractal method predicating oil resources in the Jinhu Sag, North Jiangsu Basin[J]. Chinese Journal of Geology, 2006, 41(4): 578-585. DOI:10.3321/j.issn:0563-5020.2006.04.003
[11]
鞠玮, 侯贵廷, 肖芳锋. 墨西哥湾盆地陆棚区油气田数量与储量规模的分形分析[J]. 北京大学学报(自然科学版), 2011, 47(6): 1049-1055.
JU Wei, HOU Guiting, XIAO Fangfeng. Fractals of oil and gas field reserve distribution in outer continental shelf (OCS), Gulf of Mexico Basin[J]. Acta Scientiarum Naturalium Universitatis Pekinensis, 2011, 47(6): 1049-1055.
[12]
肖芳锋, 侯贵廷, 李乐, 等. 世界油气田规模分布的分形研究[J]. 应用基础与工程科学学报, 2011, 19(1): 95-103.
XIAO Fangfeng, HOU Guiting, LI Le, et al. Fractal study on oil and gas field size distribution of the world[J]. Journal of Basic Science and Engineering, 2011, 19(1): 95-103. DOI:10.3969/j.issn.1005-0930.2011.01.011
[13]
柳广弟, 孙明亮. 剩余压力差在超压盆地天然气高效成藏中的意义[J]. 石油与天然气地质, 2007, 28(2): 203-208.
LIU Guangdi, SUN Mingliang. Significance of excess differential pressure in highly efficient gas accumulation in over-pressured basins[J]. Oil & Gas Geology, 2007, 28(2): 203-208. DOI:10.3321/j.issn:0253-9985.2007.02.011
[14]
赵贤正, 蒋有录, 金凤鸣, 等. 富油凹陷洼槽区油气成藏机理与成藏模式:以冀中坳陷饶阳凹陷为例[J]. 石油学报, 2017, 38(1): 67-76.
ZHAO Xianzheng, JIANG Youlu, JIN Fengming, et al. Hydrocarbon accumulation mechanism and model of sub-sags in hydrocarbon-rich sag: a case study of Raoyang sag in Jizhong depression[J]. Acta Petrolei Sinica, 2017, 38(1): 67-76. DOI:10.3969/j.issn.1671-4067.2017.01.022
[15]
蒋有录, 万涛, 林会喜, 等. 成藏期剩余压力与储层排替压力下限耦合恢复油气成藏过程:以济阳坳陷车西洼陷为例[J]. 石油学报, 2011, 32(2): 265-272.
JIANG Youlu, WAN Tao, LIN Huixi, et al. Reconstruction of hydrocarbon accumulation process by the matched relationship between surplus pressure and lower limit of reservoir displacement pressure in the hydrocarbon accumulation period: taking Chexi Sag in Jiyang Depression as an example[J]. Acta Petrolei Sinica, 2011, 32(2): 265-272.
[16]
刘池洋, 张东东. 盆地复杂系统特征与研究思想和方法论[J]. 西北大学学报(自然科学版), 2009, 39(3): 350-391.
LIU Chiyang, ZHANG Dongdong. Characteristics, thought of study and methodology of complex system of basin[J]. Journal of Northwest University(Natural Science Edition), 2009, 39(3): 350-391.
[17]
彭少芳, 张昭. 线性和非线性非平衡态热力学进展和应用[M]. 北京: 化学工业出版社, 2006: 8-16.
[18]
曾丹苓. 工程非平衡热动力学[M]. 北京: 科学出版社, 1991: 102-125.
[19]
李如生. 非平衡态热力学和耗散结构[M]. 北京: 清华大学出版社, 1986: 127-160.
[20]
徐国宾, 练继建. 流体最小熵产生原理与最小能耗率原理(Ⅱ)[J]. 水利学报, 2003(6): 43-47.
XU Guobin, LIAN Jijian. Theories of the minimum rate of energy dissipation and the minimum entropy production of flow (Ⅱ)[J]. Journal of Hydraulic Engineering, 2003(6): 43-47. DOI:10.3321/j.issn:0559-9350.2003.06.007
[21]
赵丽娜, 徐国宾. 基于广义流和广义力的河流能耗率推导[J]. 天津大学学报(自然科学与工程技术版), 2015(12): 1126-1129.
ZHAO Lina, XU Guobin. Derivation of energy dissipation rate of river based on generalized flux and generalized force[J]. Journal of Tianjin University(Science and Technology), 2015(12): 1126-1129.
[22]
王铮. 利用耗散结构理论分析河道演变[J]. 地理科学, 1989, 9(2): 173-180.
WANG Zheng. Analysis of river channel change with theory of dissipative structure[J]. Scientia Geographica Sinica, 1989, 9(2): 173-180.
[23]
张学文. 组成论[M]. 合肥: 中国科学技术大学出版社, 2003: 176-182.
[24]
王志宏, 李建明. 饶阳凹陷异常高压与油气成藏关系[J]. 岩性油气藏, 2014, 26(6): 15-19.
WANG Zhihong, LI Jianming. Abnormal high pressure and its relation to hydrocarbon accumulation in Raoyang Sag[J]. Lithologic Reservoirs, 2014, 26(6): 15-19. DOI:10.3969/j.issn.1673-8926.2014.06.003
[25]
刘华, 张丰荣, 蒋有录, 等. 饶阳凹陷洼槽区地层压力特征及成因机制[J]. 中国石油大学学报(自然科学版), 2016, 40(4): 37-46.
LIU Hua, ZHANG Fengrong, JIANG Youlu, et al. Characteristics and genetic mechanism of formation pressure in sags of Raoyang Depression[J]. Journal of China University of Petroleum(Edition of Natural Science), 2016, 40(4): 37-46. DOI:10.3969/j.issn.1673-5005.2016.04.005