油气地震勘探的根本目标是由地表观测和人工激发的地震波场推演或估计地下岩石介质的几何结构及弹性参数信息, 并基于这些信息进行油气藏描述与评价。这个根本目标从油气地震勘探起始至今从未改变, 相信今后也不会改变。

当前, 在“两宽一高”地震数据采集技术、全波形反演(FWI)和最小二乘逆时深度偏移(LS_RTM)成像技术、高性能计算机技术快速进步的促进下, 油气地震勘探逐渐进入了精确的反演成像阶段。大数据、机器学习与人工智能领域的技术进步, 更进一步提高了业界对高精度地震波成像和油藏描述的预期。

但是, 精确地震成像的代表性技术, 即全波形反演技术和最小二乘逆时深度偏移成像技术经过近十几年的发展依然很难实用化, 根源在于待反演弹性参数与实测波场之间的关系不是线性的。很多学者试图解决如何求解强非线性问题^[1-8]。另外, 勘探地震数据采集非常特别, 仅在地表进行激发和观测, 而近地表介质环境异常复杂, 激发与接收环境变化剧烈, 使得接收的地震子波空变严重, 与地表相关的噪声极强, 数据信噪比很低, 而且空间采样不规则; 空间采样孔径和地震子波的频带有限; 地震波场包含各种复杂的矢量波场等。地下岩石介质的弹性参数变化更复杂, 不同探区、甚至同一探区不同区域之间弹性参数都不一致。这就使得待反演参数与实测波场之间的关系异常复杂, 即便用双相介质波动方程也难以反映这种复杂关系。事实上, 对于模型驱动的参数估计问题, 选择复杂的正问题试图解释或描述数据中更多的波现象, 必然会引入更多的待反演参数, 从而显著增加参数空间的规模, 反问题的多解性和难以收敛性会更加严重。

从地震勘探学的技术发展历史可以看出, 业界正在使用的地震波成像技术有其独到的特征。这种独特性首先来自于地震勘探面对的特殊地下介质情况, 即以沉积形成的层状介质为主, 中间夹杂着构造运动和火山活动等地质原因引起的介质局部异常变化。地表激发的地震波, 在这样的介质中传播且在地表接收, 具有以反射波为主的特征。地震勘探学独创的多次覆盖观测方式进一步催生了经典的地震数据成像处理流程及相关关键技术, 其中关键的技术环节包括: 初至波或折射波射线层析; 静校正、动校正及均方根速度估计; 叠前时间偏移; 基于共成像点道集(common image gather, CIG)剩余时差的反射波射线层析和叠前深度偏移。常规的地震波成像技术体系本质上是奠基在水平层状介质和横向缓变介质中波传播理论基础上的, 各种波现象的时距关系、高频近似下的射线理论和Zoeppritz方程共同建起了走向地震波成像的桥梁。正是这样的处理流程和关键技术支撑了石油工业发展到今天, 至今还在发挥主体作用。

复杂介质情形下的地震波成像问题, 本质上就是基于实测波场的参数估计问题。探区地下介质被视为一个系统, 在地表激发人工震源激励地下介质系统, 产生系统输出, 仅在地表用检波器检测系统输出。这代表了对地震勘探系统的物理抽象。波动方程建立起了从对系统激励到系统响应的数学物理关系的抽象描述。基于波动方程建立的数学物理关系, 勘探介质系统的弹性参数估计问题(或地震波反演成像问题)可以被视为贝叶斯(Bayes)参数估计理论在勘探地震中的一个应用实例。贝叶斯参数估计理论是一个很完美的理论框架, 可以处理任意复杂介质情况下的参数估计问题。但是, 只有在对应的反问题是凸问题时, 参数估计反问题的求解才是收敛和稳定的, 解的意义才是明确的。遗憾的是, 实际地震观测得到的复杂波场和待反演参数之间的关系是强非线性的, 对应的参数估计反问题也不再是凸的, 可能有非常多的极值点, 而对于强非线性反问题目前并没有有效的数学工具来求解。这是当前全波形反演技术(也包括迭代最小二乘逆时深度偏移成像技术)难以在地震勘探中实用化的根本原因。

从地震勘探的物理实质出发, 提出若干凸问题的组合来逼近强非线性反问题的解是正确的逻辑路线。实际上, 这也是经典的地震数据成像处理流程隐含的思想逻辑。在当前油气勘探对地震波成像精度要求越来越高的情况下, 如何将贝叶斯估计理论下的全波形反演技术实用化, 既有很高的理论意义, 也能有效提升地震波成像精度, 进而改善油藏描述精度。

基于上述认识, 我们提出了特征波反演成像(characteristic wave inversion, CWI)的基本思想^[9], 即不追求对实测波场中全部波现象的模拟, 而是模拟其中的部分特征波场。特征波场的定义为: 时空局部的、单震相的和带方向的带限波场。还可以进一步地将特征波场退化为波形包络, 甚至到达时。特征波场与待反演参数之间有近似线性关系, 基于这种关系可以定义一个较凸的反问题。特征波反演成像的本质是基于特征波场且满足同因之果差异度量要求定义一组更凸的泛函, 将强非线性的全波形反演问题转化为一组更凸的反问题。更进一步地, 基于空间信息融合理论, 发展新的宽带波阻抗建模方法。

总之, 利用“两宽一高”地震观测数据, 做好精确的子波处理, 基于现代信号分析最新方法理论开展保频带和保振幅的噪声压制, 基于特征波反演成像技术流程和关键技术实现高精度速度建模和宽带反射系数成像, 基于空间信息融合方法发展出宽带波阻抗成像方法和技术, 实现油藏精细描述, 是在更高层次上走向精确地震勘探的道路。

1 上一代高精度地震勘探道路

李庆忠^[10]在《走向精确勘探的道路-高分辨率地震勘探系统工程剖析》一书中明确指出: “根据国内外发展动向, 从技术上来说, 我认为今后地震勘探将立足于三维地震及高分辨率勘探两大支柱上”。他强调, 高分辨率问题的研究主要集中于地震子波的优势频带宽度和相对波阻抗的频带宽度上, 它们是高分辨率地震勘探的核心问题, 高分辨率勘探是一个系统工程。

20世纪90年代以前, 地震数据采集以二维观测系统为主。原则上, 这样的观测系统只能对水平层状介质或横向缓变速介质进行有效成像, 并基于成像结果进行油藏描述和刻画。而且近地表条件不能过于复杂, 满足地表一致性假设条件。

事实上, 地表一致性假设成立、地下介质横向缓变速, 与此对应的折射波(和回转波)时距关系、反射波及散射波时距关系、高频近似下的波传播射线理论、Zoeppritz方程及其各种简化形式、一维波动方程, 以及据此发展出的各种地震波成像处理技术均代表了上一代地震勘探的技术道路。

在上一代地震勘探的技术路线中, 基于折射波时距关系的折射波层析、基于回转波时距关系的初至波层析、静校正、基于反射波时距关系的动校正及倾角时差校正(NMO/DMO)及均方根速度分析、基于散射波时距关系的Kirchhoff积分叠前时间偏移(PSTM)构成了地震波成像的基础。2000年前后逐渐普及的三维地震数据采集以及基于高频近似下射线理论的Kirchhoff积分叠前深度偏移(PSDM)和共成像点道集射线理论反射波层析速度反演, 将地震波成像精度提高了一大步。基于Zoeppritz方程及其各种简化形式的AVA分析与反演、基于1D波动方程的宽带波阻抗反演以及基于1D褶积模型的道积分, 共同促进了油藏描述技术的进步。各种属性提取技术提升了油藏含油气性评价结果的可靠性。总之, 2000年后十几年的时间, 石油工业的繁荣正是奠基在这样的地震波成像处理技术的基础上。

事实上, 到目前为止, 上述各种描述波场规律的数学物理方程, 以及基于它们发展的地震波成像方法, 依然是实际地震数据处理中常用的核心技术。

上一代高精度地震勘探技术系列中, 子波处理是极为受重视的领域。横向缓变速介质情形下的油气地震勘探中, 高分辨率地震勘探主要体现在垂向分辨率的提升上, 垂向分辨率的提升主要取决于激发与接收地震子波的分辨率。

实际中, 检波器接收到的同相轴上的地震子波的影响因素很多, 这可以视为一个综合褶积过程, 一般用如下模型表示。

$ w(t)=S(t) * c(t) * a(t) * m(t) $

(1)

式中: S(t)是震源激发后稳定的波前面上的波形; c(t)既包含各种近地表因素, 也包括震源本身和检波器本身等的影响; a(t)是吸收衰减因素的影响; m(t)是鬼波、鸣震的影响。子波处理的核心问题是消除各种影响, 仅保留S(t)。最后, 还要把S(t)变成零相位子波。

最困难的是对c(t)的消除。到目前为止, 地表一致性反褶积依然是消除c(t)影响的主要方法。事实上, 地表一致性振幅校正和相位(时差)校正是基于反射子波的统计结果进行的, 不是基于波传播的物理进行的。它只是按照统计规律消除了道间反射子波的振幅跃变和相位(时差)跃变。一般地, 自然界不偏好跃变, 这些高频变化往往被认为是非物理因素导致的。这种基于地表一致性假设和统计规律的处理方法是一种无奈的选择, 吸收衰减因素和鬼波、鸣震影响则按照物理规律消除。

本质上, 在横向缓变速介质假设下, 估计地下介质的(角度)反射系数才是高分辨油气勘探最本质的目标。当然, 直接估计横向缓变速介质情形下的宽带(绝对)波阻抗更有利于高精度油藏描述。其中: 最本质的困难包括近地表相关的影响因素、震源子波的带限和观测孔径的带限。很多情况下, 与近地表相关的影响因素的消除成为高精度地震勘探的真正瓶颈。

与近地表相关的各种因素, 除了对地震子波有显著的影响外, 对高分辨率地震勘探的最大障碍就是如何保真地压制随机噪声和不能用于成像的与近地表相关的波现象。面波以及面波散射波是最主要的与近地表相关的波现象, 多次折射波和鸣震也是很显著的噪声。海上油气勘探中, 深水和浅水情形下, 噪声特点存在显著不同。主要表现为: ①深水情形下, 鬼波和其它自由地表相关多次波是主要噪声波场; ②浅水情形下, 噪声波场要复杂得多, 鬼波、水体相关多次波及长程多次波是主要的噪声波场; ③极浅水情形下, 会有自由表面、极浅水体及海底附近岩石介质系统产生的复杂波现象。实际上, 只有一次反射波可以用来成像。各种去噪方法与技术, 不仅要去除各种近地表相关的噪声, 还要保持低频段和高频段等非优势频段的有效信号成分, 尽可能提升全频带的信噪比。地震子波的分辨率取决于全频带的频谱特征, 尤其是全频带的振幅谱特征^[11]。为了提高分辨率, 仅提升优势频段的信噪比不是理论上的最佳追求, 往往是工程实践上的无奈选择。

从上述分析可以清楚地看出, 李庆忠^[10]提出的上一代高精度地震勘探的技术要点可概括为, 消除共成像点道集中各角度(及各偏移距)上反射子波的所有影响因素, 保持反射子波的宽频带, 保持AVA(AVO)关系。

据此, 可以总结出上一代高精度地震勘探的技术流程及关键技术点为:

1) 宽带地震子波的激发与接收, 尤其是重视低频的激发。期望的反射地震子波是具有2~80Hz频宽、5个以上倍频程的有效频宽;

2) 消除近地表影响的地震子波处理, 包括地表一致性振幅校正、相位校正和子波整形校正等;

3) 消除近地表相关噪声且保频带的去噪处理, 包括鸣震(鬼波)的压制;

4) 吸收衰减因素的消除;

5) 迭代的叠前偏移与成像道集速度分析;

6) 成像道集后处理, 主要是剩余时差校正、子波零相位校正和最优叠加。期望成像子波具有2~80Hz、5个以上倍频程的有效频宽;

7) 1D(宽带)波阻抗反演或道积分。

2 新一代高精度地震波成像的理论方向

2000年后, 勘探地震技术发展显著加速。勘探地震学的核心问题逐渐得以聚焦, 即将三维油气探区视为一个弹性介质系统, 人工震源作为激励该系统的输入, 地表观测的波场作为该系统的输出(理论上讲, 应该在三维油气探区的所有边界上布设地震波震源和接收器, 才能采集到系统的完整激励和输出), 该实际物理系统被认为可以由人为建立的数学物理方程(弹性波方程)来描述, 然后在上述正过程描述的基础上, 基于贝叶斯估计理论对实际物理系统的弹性参数(主要是纵波速度、横波速度和密度, 实际上更多的是角度反射系数和纵波阻抗)进行估计。在弹性参数估计(通常称地震波反演成像)结果的基础上, 结合岩石物理关系, 进行岩石物性参数的估计, 进而对储层含油气性进行评价, 最后进入油气开发阶段。这是当今油气地震勘探的基本理论和技术逻辑。

据此, 油气地震勘探的核心问题可以描述为: 将地表(包括井中)观测的波场(叠前数据)视为随机信号, 一定的数学物理方程(波动方程)可以模拟(预测或解释)实测波场中的波现象(部分), 预测误差满足一定的概率分布(高斯分布), 假设有来自井中或其它来源的关于要反演参数的先验信息, 将该先验信息也视为随机的(实质上这仅仅反映了先验信息的不确定性), 基于贝叶斯公式获得要估计参数满足的后验概率密度函数, 后验概率密度函数最大(MAP)时, 对应的参数估计结果(MAP原则下的反演解)就被认为是地震波反演问题的解。这就是全波形反演(以及线性层析成像和最小二乘逆时深度偏移成像技术)的理论基础。也可以说是整个地震数据分析的理论基础。

地震波成像的本质可以用如下贝叶斯公式进一步阐述。

$ \rho_{\boldsymbol{M}}(\boldsymbol{m} \mid \boldsymbol{d})=\frac{\rho(\boldsymbol{d}, \boldsymbol{m})}{\rho(\boldsymbol{d})}=\frac{\rho(\boldsymbol{d} \mid \boldsymbol{m}) \rho(\boldsymbol{m})}{\rho(\boldsymbol{d})} $

(2)

后验概率密度函数表达了基本的贝叶斯公式。式中: ρ(d|m)是反映反演参数与数据(或波场)之间预测关系不确定性的先验概率函数, 既包括了波动方程不合理导致的对波场预测的不确定性, 也包括了仪器对波场观测的不确定性; ρ(m)代表了对反演参数先验信息的不确定性; ρ(d)代表对实测波场的先验信息。一般地, 认为没有关于实测波场的先验信息, 因此, 假设ρ(d)是均匀分布的, ρ_M(m|d)显然是已知观测数据情形下, 反演参数所要满足的条件后验概率密度函数。

假若可以得到条件后验概率密度函数ρ_M(m|d), 它规定的条件期望及条件方差就给出了反演参数的估计结果$\mathit{\boldsymbol{\hat m}} $和精度评价C_M。

$ \left\{\begin{array}{l} \hat{\boldsymbol{m}}=\int \boldsymbol{m} \rho_{\boldsymbol{M}}(\boldsymbol{m} \mid \boldsymbol{d}) \mathrm{d} m \\ \boldsymbol{C}_{\boldsymbol{M}}=\int(\boldsymbol{m}-\hat{\boldsymbol{m}})(\boldsymbol{m}-\hat{\boldsymbol{m}})^{\mathrm{T}} \rho_{\boldsymbol{M}}(\boldsymbol{m} \mid \boldsymbol{d}) \mathrm{d} \boldsymbol{m} \end{array}\right. $

(3)

但是, 在极高维情况下, 公式(3)中的积分是很难计算的。因此, 用条件后验概率密度函数ρ_M(m|d)最大对应的参数估计结果作为要反演的结果, 即:

$ \hat{\boldsymbol{m}}_{\mathrm{MAP}}=\underset{\boldsymbol{m}}{\operatorname{argmax}} \rho_{\boldsymbol{M}}(\boldsymbol{m} \mid \boldsymbol{d}) $

(4)

在条件后验概率密度函数ρ_M(m|d)是单峰的情况下, 取这样的结果作为反演解是很合理的选择。

但是, 条件后验概率密度函数ρ_M(m|d)很难确定。为此, 假设先验概率密度函数ρ(d|m)和ρ(m)都是高斯分布, 条件后验概率密度函数ρ_M(m|d)也就是高斯的。因此, 就有条件后验概率密度函数ρ_M(m|d)的最大化等价于如下误差泛函的最小化, 即:

$ \begin{array}{c} \hat{\boldsymbol{m}}=\underset{\boldsymbol{m}}{\operatorname{argmin}}\{S(\boldsymbol{m})\}=\underset{\boldsymbol{m}}{\operatorname{argmin}}\left\{( \boldsymbol{m} - \boldsymbol{m} _ { \text { prior } } ) ^ { \mathrm { T } } \boldsymbol { C } _ { \boldsymbol { M } } ^ { - 1 } \left(\boldsymbol{m}-\right. \right.\\ \left.\left.\boldsymbol{m}_{\text {prior }}\right)+\left[\boldsymbol{g}(\boldsymbol{m})-\boldsymbol{d}_{\mathrm{obs}}\right]^{\mathrm{T}} \boldsymbol{C}_{\boldsymbol{D}}^{-1}\left[\boldsymbol{g}(\boldsymbol{m})-\boldsymbol{d}_{\mathrm{obs}}\right]\right\} \end{array} $

(5)

式中: d_obs是观测波场; g(m)代表波场模拟算子, 体现待反演参数与波场之间的关系; m_prior是初始参数; C_M^－1代表系统参数的协方差矩阵的逆, 体现系统参数分量之间的关联性; $\boldsymbol{C}_{\boldsymbol{D}}^{-1} $代表测量数据的协方差矩阵的逆, 体现数据分量之间的关联性。据此, 参数估计问题就退化为公式(5)定义的误差泛函求极小的最优化问题。很显然, 公式(5)就是最经典的全波形反演。此为统计学家定义的参数估计理论框架。

事实上, 数学分析学家从变分原理的角度直接定义某种误差泛函, 同样可以给出公式(5)定义的最优化问题, 进而解决系统参数估计问题。但是, 统计学家给出的贝叶斯框架下的系统参数估计理论更深入地阐述了参数估计的本质。

公式(5)定义的参数估计反问题, 奠基在待反演参数与系统输出之间存在非线性关系的基础上。公式(5)可能有非常多的局部极小值, 对应非线性反问题解的严重不唯一性。

基本的解决方案是在某一个初始解(期望它靠近真解)附近, 对待反演参数与系统输出之间的非线性关系进行泰勒(Taylor)展开并在初始值m_B点做局部线性化处理。

$ \begin{array}{l} \boldsymbol{d}_{\mathrm{obs}}=\boldsymbol{g}\left(\boldsymbol{m}_{B}+\Delta \boldsymbol{m}\right)=\boldsymbol{g}\left(\boldsymbol{m}_{B}\right)+\frac{\partial \boldsymbol{g}}{\partial \boldsymbol{m}} \Delta \boldsymbol{m}+ \\ \quad \quad\quad\frac{1}{2} \frac{\partial^{2} \boldsymbol{g}}{\partial \boldsymbol{m}^{2}} \Delta \boldsymbol{m}^{2}+\cdots \\ \Delta \boldsymbol{d}_{\mathrm{obs}} \approx \boldsymbol{g}\left(\boldsymbol{m}_{B}+\Delta \boldsymbol{m}\right)-g\left(\boldsymbol{m}_{B}\right)=\frac{\partial \boldsymbol{g}}{\partial \boldsymbol{m}} \Delta \boldsymbol{m} \\ \Delta \boldsymbol{d}_{\mathrm{obs}}=G \Delta \boldsymbol{m} \end{array} $

(6)

待反演参数与系统输出之间的非线性关系进行局部线性化后, 得到参数扰动Δm与波场扰动Δd_obs之间的线性关系, 则公式(5)变成:

$ \begin{aligned} \hat{\boldsymbol{m}}=& \underset{\boldsymbol{m}}{\operatorname{argmin}}\left\{\left(\boldsymbol{m}-\boldsymbol{m}_{\text {prior }}\right){ }^{\mathrm{T}} \boldsymbol{C}_{\boldsymbol{M}}^{-1}\left(\boldsymbol{m}-\boldsymbol{m}_{\text {prior }}\right)+\right.\\ &\left.\left(G \Delta \boldsymbol{m}-\Delta \boldsymbol{d}_{\mathrm{obs}}\right)^{\mathrm{T}} \boldsymbol{C}_{\boldsymbol{D}}^{-1}\left(\boldsymbol{G} \Delta \boldsymbol{m}-\Delta \boldsymbol{d}_{\mathrm{obs}}\right)\right\} \end{aligned} $

(7)

这是经典的线性反演问题。公式(7)是地震数据分析的基础之一。地震数据分析中涉及到的很多种信号处理方法、图像处理方法以及成像方法可以说都在解决线性反演问题。地震成像结果或反演结果的分辨率由参数估计结果的后验方差决定。

由于上述参数反演的表述过于数学化, 尽管它很本质, 但地震勘探中关于高精度成像的讨论并没有按照上述理论框架进行, 而是转而讨论成像子波的频带、倍频程以及整个振幅谱特征, 以此界定是否实现了高精度地震成像。但是, 无论偏数学化的讨论还是偏工程化的讨论, 都只是视角不同而已。高精度地震勘探的核心依然是在高精度地震成像的基础上实现对油气藏的精确描述。

贝叶斯框架下的系统参数估计理论及基于此发展出的方法技术, 具有高度的包容性和广泛的适用性。所有类似的基于物理系统外的观测数据推测物理系统内的物理参数专业学科都遵循同样的贝叶斯参数估计理论。待反演参数与系统输出(波场)之间的关系(即波动方程, 尤其是积分形式表达的波动方程)是上述参数反演理论的基石。

全波形反演和最小二乘逆时深度偏移成像技术是贝叶斯参数估计理论在地震勘探中的应用范例。这样的理论方法并不假设地下介质分布是层状的、波场也不必以反射波为主。理论上, 这样的反演成像理论与方法适用性更强, 应该能得到更高分辨率的成像结果。但是, 十几年来, 全波形反演和迭代法最小二乘逆时深度偏移成像技术在实用化过程中举步维艰, 虽然投入了大量人力、物力, 但并没有在油气勘探工业界得到大规模实际应用。

从地震勘探学历史发展看, 地震波成像有其独特的观点和做法。地震波反演成像方法与技术的独特性来自于地震勘探面对的特殊地下介质情况, 即以沉积形成的层状介质为主, 中间夹杂着构造运动和火山活动等地质原因引起的介质局部畸变。波在这样的介质中传播具有自己的特点(即以一次反射波为主要的波现象)。地震勘探中的地表观测方式也深深地影响着地震波反演成像的方式与方法。因此, 地震波反演成像不能完全照搬简洁而完美的贝叶斯估计理论。地震波成像从一开始就分成背景速度估计与叠前偏移成像, 就佐证了上述观点。这不是在否认贝叶斯估计理论的普适性, 而是强调具体问题(地震波反演成像)的特殊性。

最近几年来, 我们一直在总结提炼地震勘探中地震波反演成像问题的提法和解法, 据此给出了地震波反演成像的总体定位: 是一个信息不充分情形下的系统参数估计问题。并形成了解决地震波反演成像问题的总体认知, 即不是硬解一个强非线性反演问题, 而是提出并求解一组邻近的且凸的线性反演问题。据此, 提出了特征波反演成像理论^[9], 其中包括如下核心方法技术: 特征波提取方法技术、波动理论特征透射波层析成像、波动理论特征反射波层析成像、最小二乘(方位角度)反射系数反演成像和宽带波阻抗成像方法。图 1是特征波反演成像流程和关键技术环节的示意。图 1a为(类)CMP道集初始背景速度建模, (类)CMP道集是叠前时间偏移成像道集, 尤其是经深时转换后的叠前深度偏移成像道集, 经反动校后形成的, 经过(类)CMP道集的速度分析可以得到较好的背景速度。图 1b为波动理论透射波到达时层析成像结果。图 1c为特征反射波波动理论到达时层析成像结果。图 1d为特征反射波波动理论波形层析成像结果。图 1e为常规的全波形反演结果。图 1f为最小二乘逆时深度偏移成像结果。图 1g为宽带波阻抗建模结果。

需要说明的是, 宽带波阻抗建模是我们提出的新一代高精度地震波成像的期望结果。但与老一代的精确地震勘探^[10]不同, 我们提出的宽带波阻抗建模并没有把它归结为利用波动方程的参数估计问题, 而是转化成空间信息融合问题。利用基于透射波和反射波的高精度速度反演方法估计0~8Hz的背景速度成分; 利用井中信息、偏移成像的结构信息、岩石物理信息、乃至地质构造和沉积信息, 建立背景密度模型, 二者融合成背景波阻抗。最小二乘叠前深度偏移或逆时偏移, 结合稀疏约束, 得到尽可能宽带的零角度反射系数。在空间信息融合的理论框架下, 将两部分信息融合在一起得到宽带波阻抗建模结果。宽带波阻抗反演这样一个典型的贝叶斯参数估计问题, 具有很强的非线性性(比单纯的速度反演非线性性更强), 就被转化成了稳定的空间信息融合建模问题。

图 2是宽带波阻抗建模示意, 忽略了背景密度的建模过程。图 2a是图 1中各种方法速度建模的抽道显示结果, 显示模型中间一道的速度变化情况。其中曲线的颜色与图 1中各图边框颜色一致。从图 1a至图 1e可以看出, 各种速度反演方法的精度逐步提升。图 2b是图 1f所示的最小二乘逆时深度偏移成像结果的最中间一道的宽带反射系数。图 2c展示了逐渐精细的背景阻抗与宽带反射系数融合后产生的宽带波阻抗的波数谱。最后的宽带波阻抗的波数谱用红线表示^[11]。很显然, CLAERBOUT ^[12]提出的中波数缺失被补齐。

图 2c中的中波数部分能否补齐, 很大程度上取决于野外地震数据采集情况。“两宽一高”的地震数据采集是新一代高精度地震勘探能否实现的基础。

综上所述, 可以总结出新一代高精度地震勘探的技术要点为: 基于“两宽一高”地震数据, 按照特征波成像技术流程和技术要点, 实现宽带波阻抗建模。

新一代高精度地震勘探的技术流程及关键技术点主要有: ① “两宽一高”地震数据观测。反射地震子波具有2~80Hz频宽、5个以上倍频程的有效频宽。②消除近地表影响的地震子波处理。③消除近地表相关噪声且保频带的去噪处理。④基于透射波和特征反射波的波动理论各向异性速度(Q值)反演与建模。⑤带吸收衰减补偿的宽带反射系数反演成像。⑥成像道集后处理, 主要是剩余时差校正、子波零相位校正和最优叠加。成像子波具有2~80Hz、5个以上倍频程的有效频宽。⑦基于空间信息融合的宽带波阻抗建模。

3 精确地震波成像对“两宽一高”数据采集的需求

当前, 地震勘探的目标地质体越来越小、复杂、隐蔽且埋深越来越大, 对地震波采集与处理都提出了更高的要求。不仅要求地质结构成像结果精确, 更希望有效提高储层描述与刻画的精度。“两宽一高”地震数据采集是达到上述目标的基本要求。

下面的理论分析展示了“两宽一高”地震数据采集对于高精度地震波成像必要性的本质。

考虑平面波入射到一个异常体或散射体上, 同时在远离该异常体或散射体上的接收点接收散射波场, 如图 3所示。当背景介质为常速时, 可以假设虚拟平面波沿直线穿过目标体, 产生观测平面波(投影)。目标体的像(或其波数谱)和散射场的关系可以被推导出来, 其类似于傅里叶(Fourier)投影定理。

此处仅考虑标量波方程, 速度由v(x)来描述。对于无源区域中的压力场$\tilde{\boldsymbol{u}}(\boldsymbol{x} ; \omega) $频域形式的标量波动方程(Helmholtz方程)为:

$ \nabla^{2} \tilde{\boldsymbol{u}}(\boldsymbol{x} ; \omega)+\left[\frac{\omega}{v(\boldsymbol{x})}\right]^{2} \tilde{\boldsymbol{u}}(\boldsymbol{x} ; \omega)=0 $

(8)

式中: x表示空间位置; ω表示频率。将目标函数m(x)定义成如下形式的速度扰动, 则:

$ \boldsymbol{m}(x)=1-\left[\frac{v_{0}}{v(\boldsymbol{x})}\right]^{2} $

(9)

式中: v₀为背景速度, 将公式(9)代入公式(8)得到:

$ \nabla^{2} \tilde{\boldsymbol{u}}(\boldsymbol{x} ; \omega)+k^{2} \tilde{\boldsymbol{u}}(\boldsymbol{x} ; \omega)=k_{0}^{2} \boldsymbol{m}(\boldsymbol{x}) \tilde{\boldsymbol{u}}(\boldsymbol{x} ; \omega) $

(10)

式中: $ \boldsymbol{k}=\omega / v(\boldsymbol{x}) ; \boldsymbol{k}_{0}=\omega / v_{0}$。

将总波场分为入射波场$\tilde{\boldsymbol{u}}_{0}(\boldsymbol{x} ; \omega) $和散射波场$ \tilde{\boldsymbol{u}}_{S}$ (x; ω)两部分:

$ \tilde{\boldsymbol{u}}(\boldsymbol{x} ; \omega)=\tilde{\boldsymbol{u}}_{0}(\boldsymbol{x} ; \omega)+\tilde{\boldsymbol{u}}_{S}(\boldsymbol{x} ; \omega) $

(11)

并将公式(11)代入公式(10)得到:

$ \nabla^{2} \tilde{\boldsymbol{u}}_{S}(\boldsymbol{x} ; \omega)+\boldsymbol{k}^{2} \tilde{\boldsymbol{u}}_{S}(\boldsymbol{x} ; \omega)=\boldsymbol{k}_{0}^{2} \boldsymbol{m}(\boldsymbol{x}) \tilde{\boldsymbol{u}}(\boldsymbol{x} ; \omega) $

(12)

注意: $\nabla^{2} \tilde{\boldsymbol{u}}_{0}(\boldsymbol{x} ; \omega)+\boldsymbol{k}^{2} \tilde{\boldsymbol{u}}_{0}(\boldsymbol{x} ; \omega)=0 $。公式(12)将散射场和目标函数及总波场联系在一起。利用格林(Green)定理, 可以得到公式(12)的解, 即:

$ \tilde{\boldsymbol{u}}_{S}(\boldsymbol{x} ; \omega)=-\boldsymbol{k}^{2} \int_{\varOmega} G\left(\boldsymbol{x}_{S}, \boldsymbol{x} ; \omega\right) \boldsymbol{m}(\boldsymbol{x}) \tilde{\boldsymbol{u}}\left(\boldsymbol{x}, \boldsymbol{x}_{R} ; \omega\right) \mathrm{d} x $

(13)

公式(13)中的积分是对目标体分布范围Ω进行的。公式(13)将目标体外表记录数据和未知的速度扰动m(x)联系在一起。散射场是由实际介质和背景介质间的波阻抗差异引起的。总场$\tilde{\boldsymbol{u}}(\boldsymbol{x} ; \omega) $作用到速度扰动m(x)上产生相互作用$ \boldsymbol{m}(\boldsymbol{x}) \tilde{\boldsymbol{u}}(\boldsymbol{x} ; \omega)$, 形成二次源, 散射场由此产生。理论上, 背景速度场可以任意选择。但背景场的选择直接影响散射势的大小。事实上, 速度扰动不能变化过大, 这是玻恩(Born)近似的理论要求。

对公式(13)引入玻恩(Born)近似得到

$ \begin{array}{c} \tilde{\boldsymbol{u}}_{S}(\boldsymbol{x} ; \omega)=-\boldsymbol{k}^{2} \int_{\varOmega} G\left(\boldsymbol{x}, \boldsymbol{x}_{R} ; \omega\right) \boldsymbol{m}(\boldsymbol{x}) \tilde{\boldsymbol{u}}_{0}\left(\boldsymbol{x}_{S},\right. \\ \boldsymbol{x} \boldsymbol{;} \omega) \mathrm{d} \boldsymbol{x} \end{array} $

(14)

公式(13)中的背景速度可以任意选, 但近似式(14)只有在速度扰动m(x)较小时才比较准确。格林函数G(x_S, x; ω)只有在背景速度分布很简单时(譬如常速或层状介质时), 可以解析地写出。一般地, 我们用数值近似或渐近近似(例如几何近似)来描述格林函数。

常背景速度情况下, 入射平面波可以表示成:

$ \tilde{\boldsymbol{u}}_{0}(\boldsymbol{x} ; \omega)=\mathrm{e}^{\mathrm{i} \boldsymbol{k} \cdot\left(\boldsymbol{x}_{S}-\boldsymbol{x}\right)} $

(15)

式中: k是入射波传播方向的波矢量。同样地, 常背景速度情况下, 三维自由空间格林函数为:

$ G\left(\boldsymbol{x}, \boldsymbol{x}_{R} ; \omega\right)=\frac{1}{\left|\boldsymbol{x}_{R}-\boldsymbol{x}\right|} \mathrm{e}^{\mathrm{i} \boldsymbol{k}_{R} \cdot\left(\boldsymbol{x}-\boldsymbol{x}_{R}\right)} $

(16)

将公式(15)和(16)代入公式(14)得到:

$ \begin{aligned} &\tilde{\boldsymbol{u}}_{S}(\boldsymbol{x} ; \omega)=-\boldsymbol{k}^{2} \int_{\varOmega} \frac{\mathrm{e}^{\mathrm{i} \boldsymbol{k}_{R} \cdot\left(\boldsymbol{x}-\boldsymbol{x}_{R}\right)}}{\left|\boldsymbol{x}-\boldsymbol{x}_{R}\right|} \boldsymbol{m}(\boldsymbol{x}) \mathrm{e}^{\mathrm{i}\boldsymbol{k}_{S} \cdot\left(\boldsymbol{x}_{S}-\boldsymbol{x}\right)} \mathrm{d} \boldsymbol{x} \\ &=-\boldsymbol{k}^{2} \int_{\varOmega} \frac{1}{\left|\boldsymbol{x}-\boldsymbol{x}_{R}\right|} \boldsymbol{m}(\boldsymbol{x}) \mathrm{e}^{\mathrm{i}\left[\boldsymbol{k}_{S} \cdot\left(\boldsymbol{x}_{S}-\boldsymbol{x}\right)+\boldsymbol{k}_{R} \cdot\left(\boldsymbol{x}-\boldsymbol{x}_{R}\right)\right]} \mathrm{d} \boldsymbol{x} \\ &=-\boldsymbol{k}^{2} \int_{\varOmega} \frac{1}{\left|\boldsymbol{x}-\boldsymbol{x}_{R}\right|} \boldsymbol{m}(\boldsymbol{x}) \mathrm{e}^{\mathrm{i} \omega\left[T\left(\boldsymbol{x}_{S}, \boldsymbol{x}\right)+T\left(\boldsymbol{x}, \boldsymbol{x}_{R}\right)\right]} \mathrm{d} \boldsymbol{x} \end{aligned} $

(17)

式中: T(x_S, x)代表源点到像点的旅行时; T(x_R, x)代表像点到检波点的旅行时。从源点到像点再到检波点的旅行时为:

$ T\left(\boldsymbol{x}_{S}, \boldsymbol{x}, \boldsymbol{x}_{R}\right)=T\left(\boldsymbol{x}_{S}, \boldsymbol{x}\right)+T\left(\boldsymbol{x}, \boldsymbol{x}_{R}\right) $

(18)

总的旅行时有如下的Taylor展开形式:

$ \begin{aligned} &\mathrm{e}^{\mathrm{i} \omega\left[T\left(\boldsymbol{x}_{S}, \boldsymbol{x}\right)+T\left(\boldsymbol{x}, \boldsymbol{x}_{R}\right)\right]}\\ &\approx \mathrm{e}^{\mathrm{i} \omega\left[T\left(\boldsymbol{x}_{S}, \boldsymbol{x}_{0}\right)+T\left(\boldsymbol{x}_{0}, x_{R}\right)\right]} \mathrm{e}^{\mathrm{i} \omega\left(\boldsymbol{x}-\boldsymbol{x}_{0}\right)\left[p_{S}\left(\boldsymbol{x}_{S}, \boldsymbol{x}_{0}\right)+\boldsymbol{p}_{R}\left(\boldsymbol{x}_{0}, \boldsymbol{x}_{R}\right)\right]}\\ &=\mathrm{e}^{\mathrm{i} \omega\left[T\left(\boldsymbol{x}_{S}, \boldsymbol{x}_{0}\right)+T\left(\boldsymbol{x}_{0}, \boldsymbol{x}_{R}\right)\right]} \mathrm{e}^{\mathrm{i} \bar{k} \cdot\left(\boldsymbol{x}-\boldsymbol{x}_{0}\right)} \end{aligned} $

(19)

将公式(19)代入公式(17)得到:

$ \begin{aligned} &\tilde{\boldsymbol{u}_{S}}(\boldsymbol{x} ; \omega)=-\boldsymbol{k}^{2} \int_{\varOmega} \frac{1}{\left|\boldsymbol{x}-\boldsymbol{x}_{R}\right|} \boldsymbol{m}(\boldsymbol{x}) \mathrm{e}^{\mathrm{i} \omega\left[T\left(\boldsymbol{x}_{S}, \boldsymbol{x}\right)+T\left(\boldsymbol{x}, \boldsymbol{x}_{R}\right)\right]} \mathrm{d} \boldsymbol{x} \\ &\approx-\boldsymbol{k}^{2} \int_{\varOmega} \frac{1}{\left|\boldsymbol{x}-\boldsymbol{x}_{R}\right|} \boldsymbol{m}(\boldsymbol{x}) \mathrm{e}^{\mathrm{i}\omega\left[T\left(\boldsymbol{x}_{S}, \boldsymbol{x}_{0}\right)+T\left(\boldsymbol{x}_{0}, \boldsymbol{x}_{R}\right)\right]} \mathrm{e}^{\mathrm{i}\bar{k}\left(\boldsymbol{x}-\boldsymbol{x}_{0}\right)} \mathrm{d} \boldsymbol{x} \\ &\approx-\boldsymbol{k}^{2} \frac{\mathrm{e}^{\mathrm{i} \omega}\left[T\left(\boldsymbol{x}_{S}, \boldsymbol{x}_{0}\right)+T\left(\boldsymbol{x}_{0}, \boldsymbol{x}_{R}\right)\right]}{\left|\boldsymbol{x}_{0}-\boldsymbol{x}_{R}\right|} \int_{\varOmega} \boldsymbol{m}(\boldsymbol{x}) \mathrm{e}^{\mathrm{i} \bar{k} \cdot\left(\boldsymbol{x}-\boldsymbol{x}_{0}\right)} \mathrm{d} \boldsymbol{x} \end{aligned} $

(20)

公式(19)中: 波矢量$ \overline{\boldsymbol{k}}=\omega\left[\boldsymbol{p}_{S}\left(\boldsymbol{x}_{S}, \boldsymbol{x}_{0}\right)+\boldsymbol{p}_{R}\left(\boldsymbol{x}_{0}\right.\right.$, $\left.\left.\boldsymbol{x}_{R}\right)\right]=(2 \omega \cos \theta) /\left[v\left(\boldsymbol{x}_{0}\right)\right] \boldsymbol{n}维波动理论下的傅里叶切片定 $。

针对公式(20), 记:

$ \widetilde{\boldsymbol{m}}(\bar{\boldsymbol{k}})=\int_{\varOmega} \boldsymbol{m}(\boldsymbol{x}) \mathrm{e}^{\mathrm{i} \bar{\boldsymbol{k}} \cdot\left(\boldsymbol{x}-\boldsymbol{x}_{0}\right)} \mathrm{d} \boldsymbol{x} $

(21)

很显然, 公式(21)是关于目标函数(即速度扰动)的三维傅里叶变换, 可以认为公式(21)是三维波动理论下的傅里叶切片定理。只有波矢量$ \overline{\boldsymbol{k}}=(2 \omega \cos \theta) /$ [v(x₀)]n定义的所有波数谱成分都存在时, 或者说只有所有角度的炮检组合及所有频率成分都存在时, 才能由上式的傅里叶反变换得到速度扰动量的精确解。很显然, 这取决于入射波与散射波的张角展布和频率成分的展布。波动理论层析中的傅里叶切片是由波矢量$ \overline{\boldsymbol{k}}=(2 \omega \cos \theta) /\left[v\left(\boldsymbol{x}_{0}\right)\right] \boldsymbol{n}$决定的, 也就是由观测系统决定的。

至此, 理论上证明了精确的地震波反演成像必须要有宽带和宽方位的地震观测数据。至于高密度采集, 在上述理论分析中体现不出来。高密度采集对高精度成像的贡献主要体现在无假频地压制各种噪声和无假频的波动方程反演成像两个方面。噪声对地震波反演成像的影响与对地震波偏移成像的影响是非常不同的, 不满足高斯分布的强噪声的存在完全可以使得反演算法不收敛。噪声压制在全波形反演和最小二乘逆时深度偏移成像为代表的成像技术阶段具有更特别的意义和重要性。

真正实现“两宽一高”地震数据采集是十分困难的。陆上, 尤其是山前带(或山地)地震勘探更是如此。最核心的是当前的采集装备水平很难适应复杂的地表条件。节点(单点)地震数据采集技术可以确定地认为是陆上地震勘探中“两宽一高”数据采集的必要技术。智能无线节点检波器技术的逐渐成熟会大幅提升“两宽一高”地震数据采集的可实现性。但是, 经济高效的宽频震源技术的发展始终比较缓慢。

陆上地震勘探中, 高精度成像的最大障碍在于复杂近地表条件变化引起的强噪声和地震子波的空间时间变化剧烈。近地表相关的噪声主要由震源以及震源附近岩性变化和地表高程变化而引起。检波器端对噪声的贡献相对较小。复杂近地表条件下的噪声在炮集上往往表现为强振幅、低频、短连续和低视速度的特点, 假频现象非常明显。如此特征的噪声目前缺乏有效的压制技术, 对后续的地震波成像, 无论偏移成像或是反演成像都会造成严重干扰。针对这样的复杂波现象(或称为噪声), 首先, 要尽可能无假频地记录下来, 然后, 构造合适的去噪方案尽可能地压制。复杂地表条件下, 有效地采集噪声在“两宽一高”地震数据采集中至关重要。如何弱化不同炮检对间地震子波的不一致也是陆上地震数据采集时要仔细考虑的问题。

事实上, “两宽一高”地震勘探技术是一项综合性的技术, 是今后油气地震勘探领域居于中心地位的技术系列。它不仅是指“两宽一高”地震数据采集技术, 挖掘出海量观测数据中包含的与储层有关的信息, 对油气储层更精确地描述, 并拓展到储层开发阶段, 提高整个勘探过程的经济效益, 才是“两宽一高”地震勘探技术的真正目标。

4 新一代精确地震勘探的分辨率

李庆忠^[10]判断“三维地震和高分辨率勘探”是当时的技术支柱是非常有远见的。不过, “两宽一高”与“全波形反演及最小二乘逆时偏移成像技术”时代, 高分辨率的概念需要进一步拓展。

本质上, 地震波反演成像结果的分辨率是其描述地下介质能力的体现。地震数据分析的核心理论与方法基础是线性反演。无偏与方差最小的参数估计结果是最好的理论结果。公式(7)的反演结果可以表达为:

$ \hat{\boldsymbol{m}}=\left(\boldsymbol{G}^{\mathrm{T}} \boldsymbol{C}_{D}^{-1} \boldsymbol{G}+\boldsymbol{C}_{M}^{-1}\right)^{-1}\left(\boldsymbol{G}^{\mathrm{T}} \boldsymbol{C}_{D}^{-1} \boldsymbol{d}^{\mathrm{obs}}+\boldsymbol{C}_{M}^{-1} \boldsymbol{m}_{\text {prior }}\right) $

(22)

而$ \mathit{\boldsymbol{\hat m}}$的精度(也可以称为反演解的分辨率)由后验方差$ \hat{\boldsymbol{C}}_{M}=\left(\boldsymbol{G}^{\mathrm{T}} \boldsymbol{C}_{D}^{-1} \boldsymbol{G}+\boldsymbol{C}_{M}^{-1}\right)^{-1}$决定。可以看出, 它们分别是高斯型条件后验概率密度函数$ \rho_{M}\left(\boldsymbol{m} \mid \boldsymbol{d}^{\mathrm{obs}}\right) \sim N$$ \left(\hat{\boldsymbol{m}}, \hat{\boldsymbol{C}}_{M}\right)$的条件期望与方差。这是当前地震勘探中地震波成像结果及其分辨率的最根本解释。

一般地, 方差被解释为海森(Hessian)矩阵的逆。也就是说, 海森矩阵决定了反演成像的分辨率。

对于玻恩近似意义下的最小二乘逆时深度偏移成像而言, 海森矩阵的每一行对应地下任一成像点的点扩散函数(point spreading function, PSF)。该函数表现为带通滤波性质。正是点扩散函数算子褶积在真像上产生了模糊的像(常规偏移成像结果)。因此, 点扩散函数反褶积形成了提高成像分辨率的基本理论基础。

上述地震波反演成像分辨率的理论概念讨论很难在勘探地震实践中得到广泛应用。为此, 在地震勘探中提出激发地震子波、观测地震子波和成像子波的概念, 便于指导高分辨率勘探的实践。水平或缓横向变化介质情形下, 用激发地震子波和观测地震子波的概念就可以讨论它们对薄层的分辨能力。复杂介质情形下, 依据偏移成像后的成像子波分辨率定义垂向分辨率, 各种子波的完整谱特征决定了它们的分辨率和对地层的分辨能力^[11], 而不是优势频带或主频。尤其是子波中低频成分的存在, 对油气勘探中岩性地层的分辨有着举足轻重的地位。因此, 总体而言, 2~80Hz(5个以上倍频程)的成像子波可以被作为实际情况下高分辨率地震勘探的追求。

以“两宽一高”和波动理论反演成像为代表的地震勘探时代, 地震波成像的目标应从高分辨率的角度成像道集延伸到宽带波阻抗成像。宽带波阻抗才是进行高精度油藏描述更合适的成像结果。

图 4对比了缺低频情形下的归一化带限波阻抗与带限反射系数。可以看出, 带限反射系数在缺低频的情形下, 分辨率降低不明显, 可以很好地描述反射层的反射系数变化。但缺低频情形下的带限波阻抗分辨率降低十分明显, 与真实的波阻抗之间差异十分显著。说明在上一代的精确地震勘探中, 带限地震子波的研究备受关注。从激发阶段就提出用小井深、小药量、高频检波器保证观测到包含高频的反射子波; 成像处理阶段, 各种反褶积方法进一步提升处理后子波的高频, 各种高精度的静校正和动校正方法保证不伤害高频成分。这些做法符合时代的工程实际。

以“两宽一高”和波动理论反演成像为代表的地震勘探时代从数据采集到子波处理、反演成像有新的逻辑。贝叶斯框架下的参数反演已经对高精度地震波成像的本质理论阐述得十分清楚。由于它难以直接指导高精度勘探的工程实践, 转而将2~80Hz(5个以上倍频程)的成像子波作为实际情况下高分辨率地震勘探的追求, 更进一步地, 将0~45Hz的(绝对)宽带波阻抗成像作为新一代高分辨率地震勘探的目标。

图 5至图 9对比了不缺低频情形下的归一化带限波阻抗与带限反射系数。可以明显地看出, 在不缺失低频的情况下, (绝对)宽带波阻抗对地层岩性参数变化有可靠的直接指示作用。而带限反射系数依然和缺低频时的带限反射系数差异不大。因此, 将(绝对)宽带波阻抗成像作为目标可以大幅提升油藏描述的精度。

图 10至图 13展示了(绝对)宽带波阻抗在描述岩性变化中的优势。可以看出, 0~45Hz的(绝对)宽带波阻抗已经可以很好地反映地下介质的空间变化规律。将0~45Hz的(绝对)宽带波阻抗成像作为当前及今后精确地震勘探的道路是合理的选择。

5 新一代精确地震勘探中地震波成像处理关键技术

与其它广义遥感方法不同, 地震勘探的特殊性体现在以下3个重要方面: ①人工震源在复杂的近地表介质环境中激发; ②地震波在地质历史中形成的复杂介质中传播; ③检波器在复杂地表环境下以空间非规则的布设方式进行地震波场采集。脉冲源激发到形成稳定的震源子波, 构成波前面之间经历了复杂的物理过程, 地震波在近地表介质中传播形成了不能用于地下介质成像的强干扰波场, 对弱能量的低频端和高频端有效反射信号形成压制性干扰, 严重制约成像分辨率的提升。用于成像的地震波方程不能描述脉冲源激发到形成稳定的震源子波, 构成波前面之前的复杂的物理过程, 也不能描述检波器本身对地震波场的滤波作用以及检波器与介质相互作用的物理过程。以几何地震学为核心的上一代地震勘探主要利用地震波到达时信息, 尽管也很重视源端和检端存在的问题, 但还没有上升到理论认识的高度。

新一代精确地震勘探中地震波成像处理奠基在波动方程和贝叶斯参数估计理论基础上, 从地震波激发、传播到接收全过程的物理实质要有新的认识, 才能发展出有效的成像处理方法技术。

基本的成像处理方法分为4个部分考虑是合理的: ①激发和接收端的各种因素, 即震源本身及震源和介质的相互作用、检波器本身(包括检波器组合效应)及检波器和介质的相互作用, 不是波动方程可以描述的。这些因素引起炮与炮及道与道之间的反射子波在振幅、相位和频带上的变化。这部分因素的压制或消除的基本思想是地表一致性假设下的、基于统计的均衡方法, 是非物理的、数据驱动的解决方案。一般地, 被称为子波处理。目的是消除炮与炮及道与道之间的反射子波在振幅、相位和频带上的不一致。②与近地表相关的波现象(噪声)的消除以及对环境噪声的压制。近地表相关的噪声波场主要是面波及面波散射波, 折射波及多次折射波, 导波和回转波也都是不能用于中深层目标成像的噪声波场。目前, 面波、折射波和回转波对近地表速度建模的贡献得到充分重视。噪声压制的本质思想是对信号进行最佳的建模预测以达到去噪目的, 即希尔伯特(Hilbert)空间中的基函数(字典)的线性组合形成信号模型, 在贝叶斯估计理论下实现最佳信号模型的建立, 然后进行数据中所含信号的最佳预测, 进而压制噪声。高维、多域、稀疏表达和基函数字典是当今解决去噪问题的关键点。③高维地震数据规则化。宽方位地震数据的偏移距和方位角分布更不均匀, 尤其是近偏移距时更是如此。高维地震数据规则化在保真成像中有重要意义, 而仅仅重视构造成像时凸显不出真正的价值。④(各向异性)速度与Q值建模与各向异性和吸收衰减介质中的叠前深度偏移成像。折射波与回转波在大偏移距和宽方位情形下存在明显的波现象, 主要体现为初至波与早至波。利用这两种波现象进行(波动理论)层析成像获取近地表速度是非常重要的环节, 能有效地提高(各向异性)速度建模和成像的精度。全方位角度道集反射波层析成像会是当前及今后很长一段时间内宽方位地震数据(各向异性)速度建模的核心工具。能产生方位角度道集的、至少针对TTI介质的、基于OVT数据或共角度(或射线参数)数据Kirchhoff积分叠前深度偏移是宽方位高密度数据成像的首选成像方法。主要原因还是其很高的计算效率。⑤成像道集的后处理。叠前成像的本质是共成像点处不同炮检对反射子波的同相叠加, 这是高精度成像本质基础。(各向异性)偏移速度的正确性是实现同相叠加的最重要保障。另外, 共成像点道集上的去噪声、剩余时差处理和子波相位处理是保障同相位叠加的重要环节。基于上述成像后处理结果进行的相对及绝对波阻抗反演是当前开展储层描述的必要步骤。

关于与地表相关的噪声处理及地表一致性处理, 包括高维地震数据规则化技术, 在(各向异性)速度建模方法技术逐步完善的情况下, 将成为进一步提高成像保真性的瓶颈。陆上地震数据采集中, 主要的噪声是面波、直达波散射波和面波散射波。特殊情况下也会有导波散射波。这些都是很强的噪声, 与后续的反射波叠加在一起, 对保真成像的结果造成严重影响。对全波形反演及最小二乘逆时深度偏移成像的影响更是举足轻重, 极容易导致根本不收敛。这也是反演成像时代, 为什么业界对噪声特别关注的根本原因。以前, 仅仅进行构造成像的时代, 噪声的影响不是特别受关注, 因为偏移成像考虑的是叠加成像, 只要能进行同相位叠加, 就可以凸显出反射(绕射及散射)点的像。但是, 反演成像强调的是预测(或逼近), 如果不能彻底地消除预测或逼近的波现象, 就会导致反演成像的不收敛。因此, 直达波、浅层折射波、面波和导波遇到近地表异常体都会导致显著的散射波, 它们的压制需要特别关注。本质上讲, 对于上述波现象, 最好的预测器是波动方程, 信号理论中的线性及非线性预测器无论如何也不如波动方程预测器。但问题是用波动方程预测地表相关的噪声需要较精确的弹性参数模型及精确的地表形态。尽管如此, 今后高精度地去除与地表相关的噪声都需要波动方程模型预测器的参与。更高精度的地表一致性校正技术在反演成像逐渐普及的时代需要再度深入研究, 毕竟复杂而空变的震源函数是反演成像的重大障碍。高维数据规则化技术在数据规则化的同时, 将数据由直角坐标系转换到极坐标系是非常有价值的做法。与极坐标系中的成像方法结合在一起, 会较大幅度提高宽方位地震成像的保真度。

总之, 在“两宽一高”采集技术逐渐普及的时代, 成像处理技术的核心没有发生本质变化, 依然是得到高精度的地下弹性参数场(反射系数或宽带波阻抗), 满足储层解释的要求。但是反演成像对噪声的要求有本质的不同, 对子波的空变也有更严格的要求。简而言之, 近地表因素对高精度成像的影响更为凸显。

6 机器学习和人工智能在新一代精确地震勘探中的可能应用

从图 14中可以看出, 油气地震勘探问题可以抽象概括为: 利用各种不同来源和不同可信度的信息, 对油气藏进行综合评价的问题。全局来看, 应该是一个典型的人工智能和大数据分析问题。但也是一个长链条、综合性的数据分析和决策问题, 其中很多关键环节并不是大数据分析问题, 例如非常核心的地震波偏移成像与反演成像问题, 也包括1D波阻抗反演和叠前AVA反演问题等都不属于大数据分析问题。

无论如何, 统一在贝叶斯最佳决策框架下, 利用统计分析、信息融合和信息综合应用的观点, 开展当今的油气勘探, 尤其是油气藏描述及含油气性判断是毋庸置疑的。基于大数据和高性能计算技术, 机器学习与人工智能正在深刻地改变很多行业的工作及经营模式, 油气勘探行业也不会例外。

但是, 我们对机器学习与人工智能在地震数据分析领域、油气藏智能识别和含油气判断领域到底能发挥什么作用以及发挥多大作用还是要有清晰认识。

现代信号与图像分析是地震数据分析的基础, 地震数据分析的每一个环节中几乎都会用到现代信号与图像分析的知识, 譬如高维数据规则化、高维数据去噪声、包括全波形反演和最小二乘逆时深度偏移成像都建立在现代信号分析中信号预测和估计理论基础上。现代信号与图像分析与地震数据分析及反演成像目前已经紧密地联系在一起, 可以说统一在贝叶斯估计理论的基础之上。

现代信号分析理论已经发展到针对高维时空信号、基于统计估计理论和建立字典基、进行稀疏特征表达。当前备受关注的各种机器学习算法是现代信号分析理论的高级发展阶段。当今广受重视的深度神经网络学习算法, 也是在深入吸收信号稀疏特征表达思想的基础上, 由基本的神经网络算法发展而来的。

因此, 机器学习与人工智能在地震勘探中的应用, 更确切地讲, 是在地震数据分析中的应用。可以说, 在很多环节上都可以发挥作用(图 14)。当前的研究重点也应集中在某些代表性的技术点上, 譬如: 自动与智能去噪、自动与智能初至波提取、自动与智能的初始速度建模和自动与智能构造解释等方面。至于智能的油藏解释、智能的储层含油气性评价技术可能还要很长时间, 等到人工智能业态在油气勘探领域发展比较成熟, 算法也有了更高阶段的发展的时候, 实用性会有很大提高。我们认为, 创建一个能有效应用的油气勘探机器(人), 代替人类油藏评价和钻井决策, 估计还需要很长时间, 但这应该成为不懈追求的长远目标。

机器学习与人工智能算法只是提供了一种数据分析的新工具。无论是否进入大数据和人工智能阶段, 油气地震勘探的核心目标依然是精确地刻画和评价油气藏。“两宽一高”数据采集技术、广义的地震波反演成像技术(推进到宽带波阻抗成像)和多种信息融合的油藏描述与评价技术, 始终是3个核心问题。

7 结论与讨论

李庆忠提出“精确地震勘探的道路”至今已近30年。尽管绿色能源、碳中和和碳达峰等先进的能源战略理念不断地冲击油气资源的发展方向, 但是油气资源的战略地位并没有发生明显的改变。高精度地震成像理论、方法和技术发展的必要性和紧迫性并没有降低。大数据和人工智能技术的异军突起也没有改变必须对油气藏开展精确描述的根本需求。

到目前为止, 地震勘探中地震波成像处理的核心问题依然是消除近地表因素对地震子波的影响、消除近地表噪声影响和消除观测系统不规则的影响, 建立精确的背景(各向异性)速度模型和Q值模型, 进行保真和高分辨率的方位角度反射系数成像以及成像道集后处理。然后, 进入油藏描述阶段。但是, “两宽一高”地震数据采集技术、更高性能的计算机技术、全波形反演和最小二乘逆时深度偏移成像技术, 也包括机器学习和人工智能, 预示着地震勘探已经进入全新的阶段。

走向精确勘探的道路需要在新的阶段给出新描述。无论是地震波反演成像问题或是油藏描述问题, 本质上都是信息不足情况下的最佳判决问题(或最佳决策问题)。贝叶斯判决理论下基于各种信息综合得到最佳的(或可靠性和精度最高的)决策(或判断、或反演结果)是本质要求。基于以上认识, 提出具备可实现性的、新一代的“精确地震勘探的道路”是非常有必要的。我们认为, 如下技术点组合可以代表新一代的“精确地震勘探的道路”, 即“两宽一高”地震数据采集技术; 更高精度子波处理技术; 基于高维、字典基和稀疏特征表达的信号处理技术(解决去噪、数据规则化、数据压缩和去混叠等问题); 建立更精确速度和Q值模型以及估计宽带反射系数的特征波反演成像技术; 宽带波阻抗成像技术以及人工智能油藏描述技术。