对于复杂地质构造区域的叠前深度偏移成像来说, 除了选取适当的偏移算法之外, 高精度的速度模型是获得良好成像质量的另一关键因素, 而如何准确估算地层的各向异性参数则一直是速度建模中的难点及研究重点。传统的基于层析成像反演的地震速度建模方法的基本原理是通过数学优化过程调整初始速度模型使得实际观测数据与模拟数据之间的差异(即优化过程的目标函数)最小化, 进而求出某种数学意义下与实际数据吻合的最佳速度模型。对于不同的层析成像算法及不同的已知数据, 目标函数的选择往往也不尽相同, 例如, 如果已知数据是拾取的地震波旅行时, 则目标函数可定义为拾取的旅行时与正演模拟的旅行时之差; 如果数据是叠前深度偏移后输出的共成像点道集, 则目标函数可定义为相应成像道集中同相轴的拉平程度。众所周知, 与大多数反问题类似, 地震层析成像反演通常是欠定和病态的^[1-2], 而且该问题在反演各向异性介质模型参数的情况下尤为严重。另外, 由于地震资料采集的技术、环境以及成本等诸多因素的限制, 实际地震数据往往只是空间有限采集范围内的地表数据, 这也造成层析成像反演不同模型参数之间, 以及模型参数和地层深度之间存在一定的不确定性, 即存在多种模型参数组合可以使得观测数据与模拟数据在一定误差范围内相吻合^[3-5]。由于这些理论和实践的限制, 在实际地震资料处理中经常会发现, 尽管层析成像反演确实使选定的目标函数减小, 但反演求出的速度模型与其它地质和地球物理资料, 如声波测井曲线、测井标记等并不一致, 有时甚至差距很大。

为减少速度建模过程中的这种不确定性并改进层析成像反演的结果, 实际生产中往往需要综合利用已知的其它地质、地球物理、岩石物理和声波测井等资料, 如深度偏移成像的地层深度与测井得到的地层深度之间的井震深度差(以下简称井震差), 对速度模型进行进一步校正^[5-6]。工业界中一种常规做法是, 在深度偏移后对偏移速度(即层析成像反演求出的速度)进行拉伸以使地层的成像深度与测井得到的地层深度相匹配, 尽量使井震差为0, 然后将由此估算得到的拉伸系数(亦即校正系数)由井位处通过某种内插或外推算法扩展到整个模型^[7]。在实际地震资料处理中, 如果研究区中地层大致水平, 则偏移成像后的速度拉伸校正方法有效, 然而当研究区域存在倾斜的复杂地质构造时, 则需要经过层析成像反演-深度偏移-拉伸校正之间的大量反复迭代才能建立较为满意的各向异性速度模型。另外需要指出的是, 此类方法的成功与否在很大程度上取决于算法中插值或外推方法的具体选择。

地层的井震差也可作为约束条件引入到层析成像反演^[5-8]。前人的相关研究中, BAKULIN等^[8]将VTI介质层析反演的范围限定在井位附近, 并在反演过程中引入与井标记(Well marker)有关的线性约束, 然后在此基础上对地表资料和井资料进行联合反演, 并通过预处理方法将测井信息拓展到井附近的局部反演区域中, 由此使反演模型的更新遵从地质层位并受空间平滑性的约束。SPADAVECCHIA等^[5]也提出了类似方法, 并利用BP2007人工模拟数据的测试, 证明了所提方法不仅提高了叠前深度偏移的地震地层深度与井标记之间的一致性, 而且改善了偏移后共成像点道集的拉平程度。

在前人研究工作的基础上, 本文提出一种井约束条件下的各向异性层析成像反演方法, 与前人工作不同之处在于, 本文通过引入了一个所谓的“桥”函数, 将地层的井震差信息由井位(通常是稀疏的)延拓到整个模型空间, 提高建模的精度。首先对理论算法的要点进行了讨论, 之后利用模拟数据和实际数据验证了方法的有效性。

1 方法原理

本文讨论的层析成像是基于地震波高频近似, 即射线追踪的旅行时反演方法, 并在模型空间引入了正则化。如果已知井震差, 则层析成像反演目标函数C(m)可表示为:

$ C\left(\mathit{\boldsymbol{m}} \right)=\|\mathit{\boldsymbol{Am}}-{{S}_{\text{pick}}}{{\|}^{{{\gamma }_{1}}}}+{{\lambda }_{1}}\|{{\nabla }^{2}}\mathit{\boldsymbol{m}}{{\|}^{{{\gamma }_{2}}}}+{{\lambda }_{2}}\|{{W}_{\text{cal}}}-{{W}_{\text{obs}}}{{\|}^{{{\gamma }_{3}}}} $

(1)

式中:右端第一项中A即为文献[9]中公式(2)里面的层析矩阵; m为待求的速度模型参数; S_pick为叠前深度偏移共成像点道集中拾取的同相轴深度差。右端第二项由拉普拉斯算子(∇²)构成, 是对模型的正则化以避免0空间问题。由于层析反演得到的速度模型存在一定误差, 即地层的测井深度W_obs与相应的偏移成像深度W_cal并不一致, 存在深度误差, 为此在(1)式右端引入第三项, 以在优化过程中减少此误差。λ₁, λ₂和λ₃为加权系数, ‖·‖^γ₁, ‖·‖^γ₂和‖·‖^γ₃表示向量的某种范数, 通常选取为L₂^[10]。

如何在层析反演中根据地层的井震差实现对速度模型的校正？为此, 首先讨论一种简单情况, 即井位附近区域内的校正。其基本原理, 如图 1所示, 将观测的地层深度误差(即某一测井标记与偏移地震剖面上相应地层深度之差)Δz转化为零炮检距射线的旅行时误差。图中|ΔP|=Δz·cosθ为Δz在零炮检距射线上的投影, 射线由成像点P₀出发并指向地表, θ为P₀点处射线方向与测井方向之间的夹角。为便于讨论, 假设P₀为深度偏移像点位置, m₀为相应的偏移速度, 亦即层析成像反演的初始速度模型, P为地层的真实位置, 亦即测井深度, m为待求的真实速度模型。对P₀处的像进行逆偏移, 则可恢复其真实旅行时, 且满足:

$ T\left( \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{m}} \right)=T({{\mathit{\boldsymbol{P}}}_{0}}, {{\mathit{\boldsymbol{m}}}_{0}}) $

(2)

对(2)式做泰勒展开并略去高阶项, 可得:

$ \begin{align} T\left( \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{m}} \right)&=T({{\mathit{\boldsymbol{P}}}_{0}}+\Delta \mathit{\boldsymbol{P}}, {{\mathit{\boldsymbol{m}}}_{0}}+\Delta \mathit{\boldsymbol{m}}) \\ &\approx T({{\mathit{\boldsymbol{P}}}_{0}}, {{\mathit{\boldsymbol{m}}}_{0}})+\frac{\partial T({{\mathit{\boldsymbol{P}}}_{0}}, {{\mathit{\boldsymbol{m}}}_{0}})}{\partial \mathit{\boldsymbol{P}}}\Delta \mathit{\boldsymbol{P}}+ \\ &\frac{\partial T({{\mathit{\boldsymbol{P}}}_{0}}, {{\mathit{\boldsymbol{m}}}_{0}})}{\partial \mathit{\boldsymbol{P}}}\Delta \mathit{\boldsymbol{m}}+O({{\Delta }^{2}}) \\ \end{align} $

(3)

将(3)式代入(2)式, 可得井震差约束条件下的速度校正为:

$ \frac{\partial T({{\mathit{\boldsymbol{P}}}_{0}}, {{\mathit{\boldsymbol{m}}}_{0}})}{\partial \mathit{\boldsymbol{m}}}\Delta \mathit{\boldsymbol{m}}=-\frac{\partial T({{\mathit{\boldsymbol{P}}}_{0}}, {{\mathit{\boldsymbol{m}}}_{0}})}{\partial \mathit{\boldsymbol{P}}}\Delta \mathit{\boldsymbol{P}} $

(4)

式中:Δm为待求的模型参数扰动; ΔP为与已知的井震差有关的量(图 1)。如果在研究区域内井位足够稠密并且空间均匀分布, 则由(1)式及(4)式可直接求出区域内速度模型参数。

(1) 式和(4)式描述的是理想条件下的井约束反演。然而实际应用中层析成像反演解是否有效与很多因素有关, 例如前面已经提到, 如果研究区域内测井数量有限且位置分布稀疏, 直接应用(1)式和(4)式只能使井周围的局部区域的速度模型得到更新。此外, 基于大规模线性系统求解的计算效率的要求, 实际计算中通常采用共轭梯度法等数值方法进行求解, 再加上迭代次数的限制, 有时得到的速度更新量可能呈尖刺状或振荡状, 这些都是由反演的病态性导致的。

为降低反演的病态性并获得稳定的层析反演结果, 本文引入一个所谓的“桥”函数f, 可将(1)式改写为:

$ \begin{align} &C\left( \mathit{\boldsymbol{m}} \right)=\|A\mathit{\boldsymbol{m}}-{{S}_{\text{pick}}}{{\|}^{{{\gamma }_{1}}}}+{{\lambda }_{1}}\|{{\nabla }^{2}}\mathit{\boldsymbol{m}}{{\|}^{{{\gamma }_{2}}}}+ \\ &\ \ \ \ \ \ {{\lambda }_{2}}\|Df({{W}_{\text{cal}}}({{x}_{0}})-{{W}_{\text{obs}}}({{x}_{0}}), x){{\|}^{{{\gamma }_{3}}}} \\ \end{align} $

(5)

函数f(x₀, x)的作用是将井位x₀处的井震差加权外推到井周围的任意一点x, 算子D的作用是强制反演过程中地层的成像深度与测井深度保持一致。通常情况下函数f难以直接求出解析解, 但通过求解一个线性偏微分方程组可以得到其数值解。“桥”函数与井的位置以及局部地质构造有关, 如地震解释的层位、地层倾角和方位角场等。尽管理论上可将测井信息拓展到整个模型空间, 但实际应用中延拓范围则是根据井的位置分布、地震数据质量、同相轴拾取及地质解释的可信度等人工设定。

2 理论模型实验

选取SEAM模拟数据进行试算, 结果见图 2。研究区域以纵测线3225为中心, 纵测线方向长度17km, 横测线方向宽度5km, 深度8km。准确速度模型和各向异性参数ε如图 2a和图 2d所示。我们还在准确模型中加入人为扰动, 如图 2b和图 2e所示, 以模拟有误差的反演初始模型。为模拟地层的井震差, 分别使用准确速度模型和初始模型进行深度偏移, 并在偏移剖面上拾取地层的深度(图 3a和图 3b), 然后将得到的深度差作为井震差, 如表 1中速度校正前各井标记的深度误差所示。图 2c和图 2f为应用本文方法反演得到的初始模型的校正量, 可以看出, 反演结果大体上恢复了准确的速度模型。使用反演校正后的模型进行深度偏移, 结果如图 3c所示, 可以发现, 井震差在3个测量位置处基本为0, 如表 1和图 3c所示, 相应道集中的同相轴的拉平程度也得到改进, 如图 4所示。

3 实际数据应用

实际数据是在墨西哥湾采集的某一广角数据。建立的初始速度模型如图 5所示, 图中黄线表示测井位置, 位于纵测线5 800与横测线40 977交会处。在初始速度模型的深度偏移成像剖面上, 盐丘顶部的深度为2 503 m, 比相应的井标记(2 474 m)深约29 m, 井震差如图 6a所示。

对v₀(沿各向异性对称轴的速度)和各向异性参数ε进行井约束条件下的联合反演, 得到速度模型校正量。采用校正后的速度模型进行深度偏移, 结果显示盐丘的顶部深度为2 477 m(图 6b)。这意味着, 经过井约束条件下层析成像, 地层的成像深度误差从29 m减小到3 m(小于深度采样间隔)。图 7表示速度修正量以及各向异性参数修正量。图 8和图 9给出了利用初始速度模型以及更新后的速度模型偏移的共成像点道集(图中红线表示盐丘顶部(TOS)的井标记), 由图 8可见, 道集中有一些同相轴在大炮检距处不平, 这表明初始速度模型中存在误差, 而图 9上共成像点道集的平坦度得到改善。

4 结论

本文提出了一种井震差约束条件下TTI介质速度建模的方法, 通过构造一个“桥”函数, 将井震差信息从空间稀疏分布的井位延伸到其周边区域。这样的处理方式有助于减少反演算子的不适应性并提高解的稳健性。SEAM模拟数据和实际数据的应用表明, 该方法不仅提高了深度偏移的地层深度与测井得到的地层的一致性, 同时也提高了偏移后共成像点道集中同相轴的拉平程度。理论上可将测井信息拓展到整个模型空间, 但实际应用中扩展的范围是有限的, 取决于井的位置分布、数据质量、同相轴的拾取以及地质解释的可信度等因素。

致谢: 感谢Statoil公司的Gentiana Ionescu, Teresa Szydlik, Aaron DeNosaquo和Mike Cogan等同事在速度建模和深度偏移等资料处理过程中有益的讨论和重要贡献。同时感谢Statoil公司对本文工作的支持和发表许可。