地震偏移方法的研究始于20世纪70年代, 主要为叠后时间域偏移, 主流方法是单程波偏移。90年代, 国内学者普遍认为三维叠前深度偏移还只是理论上的探索, 距离应用为时尚早。而国外已经在实际资料处理中进行三维叠前深度偏移试验性处理。此时生产中应用的算法还只有Kirchhoff偏移, 越是有经验的专家越是不相信单程波偏移在这个新领域有任何实用的可能性。随后的10年是偏移理论和方法快速发展时期。仅仅两年的时间, 在深海区勘探地震资料处理市场上Kirchhoff偏移逐步被单程波偏移方法所替代, 而且人们也充分认识到了单程波偏移方法的不足。2006年前后, 逆时偏移应运而生, 极大地提高了地震资料处理质量和效果。

本文总结了叠前深度偏移发展过程中所遇到的一些重要的理论和生产问题, 以及相应发展起来的方法和技术。第1节到第3节的主体是作者和同事们在叠前深度偏移方面的工作总结。特别地, 用真振幅偏移的理论和方法作为主线来贯串。真振幅偏移理论是在生产软件研发过程中发展起来的理论工具。最初要回答的问题是:①偏移是什么？②如何考察偏移软件的正确性？③不同种类的偏移方法和结果如何对应与翻译？④如何在生产中合理地使用偏移结果？理论不仅仅是软件研发的指导, 其效力还需要通过产品质量和应用效果来体现。由于篇幅和版权的限制, 本文主要阐述理论和算法, 但是引用了大量的文献。希望读者能够将理论和文献中的应用结果相比照, 以求得到印证。第3节总结了生产实践对逆时偏移所提出的要求和从中发展出来的一些算法, 包含了“叠前道集的输出”, “各向异性”, “非弹性衰减”, “组合震源、鬼波和多次波的处理”等比较繁杂的内容。在整理过程中, 作者尽量在统一的符号框架下给予介绍说明, 采取文献中的思路, 对重要公式重新做了推导, 所导出的公式和文献相比或有细节上的差异。第4节到第6节讨论了当下受到更多关注的最小二乘偏移和全波形反演方法, 指出真振幅偏移的理论和方法是联系成像与反演的纽带。最小二乘偏移在技术上还不很成熟, 全波形反演是当前发展最快的方向, 有关的文章报告数不胜数。第6节对全波形反演的一些重要问题, 如“低频缺失”, “振幅和走时”, “反射、折射与透射”, “多参数反演”等方面做了粗略的总结。最后在第7节“总结与展望”中, 作者对全波形反演寄予了厚望。

1 从真振幅叠前Kirchhoff偏移谈起

偏移成像和速度反演都可以看作波动方程正演模拟的反问题。

对于最简单的三维常密度声波方程, 在零时刻(t=0)和地表炮点位置x_s=(x_s, y_s, 0)激发一个单脉冲震源δ(t)δ(x－x_s)。该震源产生的压力波场满足如下方程:

$ \left[ {\frac{1}{{{v^2}\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - \Delta } \right]p\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \delta \left( t \right)\delta \left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) $

(1)

式中:x=(x, y, z)为任意空间位置; v(x)是地下介质速度; $\Delta = {\partial ^2}/\partial {x^2} + {\partial ^2}/\partial {y^2} + {\partial ^2}/\partial {z^2}$是Laplace算子。在地面x_r=(x_r, y_r, 0)处压力检波器所记录下的声波信号:

$ D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = p\left( {\mathit{\boldsymbol{x}} = {\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) $

(2)

就是我们所说的人工地震信号。一般情况下, 它是一个具有五维变量(x_r, y_r; t; x_s, y_s)的数据体。

传统的偏移成像可以定义为这样的反问题:给定地震记录D(x_r; t; x_s)和一个背景速度v₀(x), 构建地下反射系数模型R₀(x; θ; φ)。注意, 此处的反射系数R₀除了随空间变量x变化之外, 还随反射角θ和反射方位角φ变化。这样, 叠前偏移成像的输入和输出都是五维数据体。

为了保证成像结果能够直接反映地下真实构造, 要求背景速度v₀(x)尽量接近真实速度v(x)。特别地, 要求由这两个速度函数所定义的波传播走时在较大尺度下尽可能一致, 这等价于要求v₀(x)和v(x)在各个传播方向上的平均慢度要很接近。

给定背景速度v₀, 相应的波场p₀满足:

$ \left[ {\frac{1}{{v_0^2\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - \Delta } \right]{p_0}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \delta \left( t \right)\delta \left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) $

(3)

在小扰动假设下, 将地震波场p分解为两部分:背景波场p₀和扰动波场δp。BLEISTEIN等在文献[1]的附录中给出了δp和反射系数之间的关系:

$ \begin{array}{*{20}{c}} {{\rm{ \mathsf{ δ} }}\hat p\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \int {{\mathit{\boldsymbol{R}}_0}\left( {\mathit{\boldsymbol{x}};\theta ;\varphi } \right) \cdot \nabla \left( {{{\hat p}_0}\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;\mathit{\boldsymbol{x}}} \right) \cdot } \right.} }\\ {{{\hat p}_0}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){\rm{d}}\mathit{\boldsymbol{x}}} \end{array} $

(4)

在(4)式中引入了Fourier变换

$ \hat p\left( \omega \right) = \int_{ - \infty }^{ + \infty } {p\left( t \right){{\rm{e}}^{{\rm{i}}\omega t}}{\rm{d}}t} $

(5)

式中:∇是对变量x求梯度; R₀(x; θ; φ)=R₀(x; θ; φ)n, n是与反射面相垂直的单位向量(正方向指向地面)。由公式(3)可知p₀就是波动方程的Green函数, 它在高频渐近意义下的表达式为:

$ {{\hat p}_0}\left( {\mathit{\boldsymbol{y}};\omega ;\mathit{\boldsymbol{x}}} \right) = A\left( {\mathit{\boldsymbol{y}};\mathit{\boldsymbol{x}}} \right){{\rm{e}}^{{\rm{i}}\omega \tau \left( {\mathit{\boldsymbol{y}};\mathit{\boldsymbol{x}}} \right)}} $

(6)

式中:τ(y; x)和A(y; x)是从x到y的走时函数和射线振幅, 分别满足程函方程

$ {\left| {\nabla \tau } \right|^2} = \frac{1}{{{v^2}}} $

(7)

和输运方程

$ 2\Delta \tau \cdot \nabla A + A\Delta \tau = 0 $

(8)

为简便起见, 我们记A_s=A(x; x_s), A_r=A(x_r; x), τ_s=τ(x; x_s)和τ_r=τ(x_r; x)。将射线表示(6)式代入(4)式就得到了扰动波场的高频渐近表达式:

$ \begin{array}{*{20}{c}} {{\rm{ \mathsf{ δ} }}\hat p\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \int {{R_0}\left( {\mathit{\boldsymbol{x}};\theta ;\varphi } \right)K\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};{\mathit{\boldsymbol{x}}_{\rm{s}}};\mathit{\boldsymbol{x}};\omega } \right) \cdot } }\\ {{{\rm{e}}^{{\rm{i}}\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)}}{\rm{d}}\mathit{\boldsymbol{x}}} \end{array} $

(9)

式中:

$ K\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};{\mathit{\boldsymbol{x}}_{\rm{s}}};\mathit{\boldsymbol{x}};\omega } \right) = {\rm{i}}\omega {A_{\rm{s}}}{A_{\rm{r}}}\nabla \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right) \cdot \mathit{\boldsymbol{n}} $

(10)

如果定义公式(9)中的算子为:

$ {\rm{ \mathsf{ δ} }}\hat p = K\left( {{R_0}} \right) $

(11)

并且定义目标函数:

$ \begin{gathered} C\left( R \right) = \frac{1}{2}\iiint {\left| {\hat D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) - } \right.} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;{\left. {K\left( {{R_0}} \right)\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)} \right|^2}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}{\rm{d}}{\mathit{\boldsymbol{r}}_{\rm{r}}}{\rm{d}}\omega \hfill \\ \end{gathered} $

(12)

求解上面的最小二乘优化问题可以得到比较一般的真振幅偏移公式^[2-3]

$ \begin{array}{*{20}{c}} {{R_0}\left( {\mathit{\boldsymbol{x}};\theta ;\varphi } \right) \sim \frac{1}{{8{{\rm{ \mathsf{ π} }}^3}}}\iint {{\rm{i}}\omega \frac{{H\left( {\xi ;\mathit{\boldsymbol{x}}} \right)}}{{{A_{\rm{s}}}{A_{\rm{r}}}\left| {\nabla \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)} \right|}} \cdot }} \\ {{{\rm{e}}^{ - {\rm{i}}\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)}}\hat D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){\rm{d}}\xi {\rm{d}}\omega } \end{array} $

(13)

其中,

$ H\left( {\xi ;\mathit{\boldsymbol{x}}} \right) = \det \left\{ \begin{array}{l} \nabla \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)\\ \frac{\partial }{{\partial {\xi _1}}}\nabla \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)\\ \frac{\partial }{{\partial {\xi _2}}}\nabla \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right) \end{array} \right\} $

(14)

称作Beylkin行列式, 它描述了地面积分坐标到地下成像坐标之间的转换关系; ξ是由x_s和x_r所定义的某种观测系统。实际应用中, 习惯用反射界面的峰值来考察反射系数的强度。为此, 引入归一化的反射系数函数:

$ \begin{array}{l} R\left( {\mathit{\boldsymbol{x}};\theta ,\varphi } \right) = \frac{{2\cos \theta }}{v}{R_0}\left( {\mathit{\boldsymbol{x}};\theta ,\varphi } \right)\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \left| {\nabla \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)} \right|{R_0}\left( {\mathit{\boldsymbol{x}};\theta ,\varphi } \right) \end{array} $

(15)

相应的真振幅偏移公式修改为:

$ \begin{array}{*{20}{c}} {R\left( {\mathit{\boldsymbol{x}};\theta ;\varphi } \right) \sim \frac{1}{{8{{\rm{ \mathsf{ π} }}^3}}}\iint {{\rm{i}}\omega \frac{{H\left( {\xi ;\mathit{\boldsymbol{x}}} \right)}}{{{A_{\rm{s}}}{A_{\rm{r}}}{{\left| {\nabla \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)} \right|}^2}}} \cdot }} \\ {{{\rm{e}}^{ - {\rm{i}}\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)}}\hat D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){\rm{d}}\xi {\rm{d}}\omega } \end{array} $

(16)

对于常用的共炮检距偏移, 炮点与检波器点之间的矢量h=x_r－x_s是固定的, ξ=(x_r+x_s)/2定义为炮点和检波点的中点。对于任意三维速度函数v(x), 共炮检距情况下的Beylkin行列式有比较复杂的射线表达形式, 参见文献[4]。文献[5]和文献[6]给出了适用于时间偏移的层状介质(v(z))假设的解析公式。

如下两种情况下的真振幅公式有比较简单的表达形式。

1) 共炮偏移。

此时偏移输出真振幅共炮道集^[7]:

$ \begin{array}{*{20}{c}} {R\left( {\mathit{\boldsymbol{x}};{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) \sim \frac{1}{{\rm{ \mathsf{ π} }}}\iint {{\rm{i}}\omega \frac{{\cos {\alpha _{\rm{r}}}}}{{{v_{\rm{r}}}}}\frac{{{A_{\rm{r}}}}}{{{A_{\rm{s}}}}}}{{\rm{e}}^{ - {\rm{i}}\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)}} \cdot } \\ {\hat D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}{\rm{d}}\omega } \end{array} $

(17)

式中:v_r=v(x_r); α_r是反射射线在检波点处与垂直方向的夹角。

2) 共反射角偏移。

此时偏移输出真振幅共反射角道集^{[1, 8]}:

$ \begin{gathered} R\left( {\mathit{\boldsymbol{x}};\theta ,\varphi } \right) \sim \frac{{16{\rm{ \mathsf{ π} }}v\left( x \right)}}{{\sin \theta }}\iint {\iiint {{\rm{i}}\omega \frac{{\cos {\alpha _{\rm{s}}}\cos {\alpha _{\rm{r}}}}}{{{v_{\rm{s}}}{v_{\rm{r}}}}}}} \cdot \hfill \\ {A_{\rm{s}}}{A_{\rm{r}}}{{\rm{e}}^{ - {\rm{i}}\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)}}\hat D\delta \left( {\varphi - \varphi '} \right)\delta \left( {\theta - \theta '} \right){\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}{\rm{d}}\omega {\rm{d}}\theta '{\rm{d}}\varphi ' \end{gathered} $

(18)

式中:脉冲函数δ(φ－φ′)和δ(θ－θ′)表示从观测系统(x_s, x_r)到地下角度域(φ, θ)的变换。

文献[5]和文献[9]给出了一些真振幅偏移的计算实例。概括认为, 真振幅偏移方法具有如下功能:

1) 通过精确的走时计算保证构造成像位置准确;

2) 通过正确的相位处理(如三维的90°相位移)保证反射极性正确;

3) 通过正确的频谱操作(如三维的ω滤波器)保证子波在成像过程中不发生畸变, 从而最大限度地保持构造的成像精度;

4) 通过正确的振幅加权补偿背景速度下的几何球面扩散效应, 使得偏移输出剖面具有定量的物理意义。

一般地, 认为由真振幅偏移所输出的道集在局部范围内与反射系数的强弱成正比。这是满足生产中AVO^[10]分析的最为可靠的成像方法。

2 走向逆时偏移

自20世纪90年代开始, 地震资料处理迅速向叠前深度偏移发展。Kirchhoff偏移虽然广泛应用于各个地区的地震资料处理, 但是对于以墨西哥湾为代表的复杂地质构造成像的效果很不理想。究其原因在于走时函数的计算。主要有两个问题:首先, 在复杂构造条件下, 连接两点x和y的射线不止一条, 可以有很多条, 甚至无穷条。这时的走时函数τ(x; y)不是一个单值函数。也就是所谓的多次到达(multi-arrival)问题。为了节省计算时间, 商业化的Kirchhoff偏移软件通常事先计算并存储走时函数, 将τ(x; y)简化成一个单值函数, 从而造成走时函数的误选(没有选取到可以成像的走时分支)和间断(从一个走时分支跳跃到另一个)。前者妨碍了有效成像, 后者导致了叠加剖面的高频噪声。第二个问题更为严重, 我们注意到高频渐近表达式(6)中的走时和振幅函数均与频率无关, 这实际上假设了入射波场的频率相对于背景速度的波数变化更为高频。换句话说, Kirchhoff偏移只适用于背景速度变化相对平滑的构造成像。因此, 工业界开始研发更为先进的用于复杂地质构造成像的方法, 并相继开发出了商业化的叠前Gauss射线束偏移和叠前单程波偏移技术。

Gauss射线束偏移^[11]可以看成是Kirchhoff偏移的升级版, 它以射线束代替射线。在Gauss射线束的框架下, 波动方程的格林函数在频率域可以表示为:

$ \hat G\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{y}};\omega } \right) = \frac{{{\rm{i}}\omega }}{{2{\rm{ \mathsf{ π} }}}}\int {{{\hat u}_{{\rm{GB}}}}\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{y}};\mathit{\boldsymbol{p}};\omega } \right)\frac{{{\rm{d}}{p_x}{\rm{d}}{p_y}}}{{{p_z}}}} $

(19)

式中:û_GB(x; y; p; ω)是在起始点归一化之后的Gauss束函数; p=(p_x, p_y, p_z)是慢度向量, 定义了地表射线的出射方向, 确定了这个方向之后, 从地面一点到达地下点的走时曲线就是唯一的, 这样射线束偏移很自然地解决了多次到达问题。

Ross Hill博士是射线束偏移方法生产化应用的主要推动者之一。经过不断的研发, 至21世纪初, 形成了从去除多次波、速度建模到三维叠前深度偏移的全套生产应用流程, 对墨西哥湾复杂构造成像的效果优于当时流行的Kirchhoff偏移。GRAY等^[12]发展了真振幅射线束偏移理论。由于射线束偏移可以只针对局部或主能量进行成像, 因而具有快速、灵活、信噪比高等特性, 研究人员相继开发出了一些改进的射线束偏移技术, 如射线束叠前深度偏移(BPSDM)^[13]以及控制射线束偏移(CBM)^[14], 等等。

严格的射线束偏移是解决多次走时问题的一种有效途径。但在实际应用中为了节省计算量, 经常通过稳相原理来减少对各个出射方向分支的搜索, 因此降低了计算精度和适用性。最重要的问题是这种方法仍然基于射线追踪, 离不开高频渐近假设, 因而对于复杂构造的成像效果仍然不够理想。

单程波方程偏移肇始于20世纪70年代CLAERBOUT的15°方程分解^[15]和GAZDAG的相位移偏移^[16]。后来的近40年间, 这种方法一直是偏移研究的主要课题之一。21世纪初, 单程波方程曾经被广泛应用于叠前深度偏移。与Kirchhoff积分偏移方法相比, 单程波方程偏移无需求解射线追踪, 避免了高频渐近逼近假设, 更适用于复杂介质和横向速度变化剧烈条件下的构造成像。与逆时偏移方法相比, 单程波方程偏移具有内存小和计算效率高的特点, 适合于大规模偏移资料处理。但是, 传统的单程波方程偏移本身存在以下问题:

1) 其偏移方法建立在直观的成像概念上, 缺乏严格的反演理论支持。因而不像Kirchhoff偏移技术那样, 具有完整的真振幅偏移理论基础。

2) 经过单程波分解后, 波场只能向下或向上传播, 因而波的传播方向受到90°的倾角限制, 无法对回转波进行成像。

为解决上述两个问题, 需要将传统的单程波分解理论建立在严格的数学理论基础之上。张关泉^[17]给出了如下单程波分解:

$ \left\{ \begin{array}{l} \left( {\frac{\partial }{{\partial z}} + {\rm{i}}\mathit{\Lambda } - \mathit{\Gamma }} \right){{\hat p}_D}\left( {\mathit{\boldsymbol{x}};\omega } \right) = \mathit{\Gamma }{{\hat p}_U}\left( {\mathit{\boldsymbol{x}};\omega } \right)\\ \left( {\frac{\partial }{{\partial z}} - {\rm{i}}\mathit{\Lambda } - \mathit{\Gamma }} \right){{\hat p}_U}\left( {\mathit{\boldsymbol{x}};\omega } \right) = \mathit{\Gamma }{{\hat p}_D}\left( {\mathit{\boldsymbol{x}};\omega } \right) \end{array} \right. $

(20)

并且证明波场p=p_D+p_U的走时和振幅函数分别满足程函方程((7)式)和输运方程((8)式)。式中算子Λ和Γ均为拟微分算子:

$ \begin{array}{l} \mathit{\Lambda } = \frac{\omega }{v}\sqrt {I + \frac{{\Delta T}}{{{\omega ^2}}}} = \frac{\omega }{v}\left[ {I - \frac{1}{{\rm{ \mathsf{ π} }}}\int_{ - 1}^1 {\sqrt {1 - {s^2}} \left( {{\omega ^2} + } \right.} } \right.\\ \;\;\;\;\;\left. {{{\left. {{s^2}{\Delta _T}} \right)}^{ - 1}}{\Delta _T}{\rm{d}}s} \right] \end{array} $

(21)

$ \mathit{\Gamma } = \frac{{{v_z}}}{{2v}}{\left( {{\omega ^2} + {\Delta _T}} \right)^{ - 1}}{\omega ^2} $

(22)

其中,

$ {\Delta _T} = {\left( {v{\partial _x}} \right)^2} + {\left( {v{\partial _y}} \right)^2} $

(23)

在此基础上, ZHANG等证明了适当选择单程波方程形式、边界条件和成像条件, 可以得到共炮域和共反射角域的真振幅偏移算法^{[9, 18-19]}。为了使单程波方程突破90°倾角的限制, ZHANG等引入并且实现了具有回转波机制的单程波偏移算法, 使单程波偏移方法具有了大倾角成像能力^[20-21]。

综上所述, 有关单程波方程偏移的问题似乎都解决了, 但在实际生产中, 更为关键的是如何数值计算单平方根算子Λ。自20世纪70年代以来, 地球物理学家在这方面做了长期不懈的努力, 提出了各种各样的单程波求解算法(见文献[22]综述), 但并未能找到一种普适的快速求解单程波方程的方法。事实上, 单平方根算子Λ是一个拟微分算子, 需要用一系列微分、积分算子的求和来逼近。其在频率-波数域的特征是:

$ \begin{array}{l} \lambda = \frac{\omega }{v}\sqrt {1 - \frac{{{v^2}\left( {k_x^2 + k_y^2} \right)}}{{{\omega ^2}}}} \\ \frac{{{\omega ^2}}}{{{v^2}}} > k_x^2 + k_y^2 \end{array} $

(24)

当ω²/v²＜k_x²+k_y²时, 特征λ为虚数, 进入倏逝波区域; 分界面ω²/v²=k_x²+k_y²是算子的奇点。由于带有奇性的算子很难找到快速一致的逼近算式, 因而给单程波方程的数值求解带来了本质困难。基于这个事实, 从单程波方程偏移向逆时偏移发展具有必然性。

1983年, WHITMORE^[23], BAYSAL等^[24]和MCMECHAN^[25]先后发表论文, 提出了一种新的偏移方法, 之后约定俗成, 定名为逆时偏移(Reverse Time Migration)。这种方法还可以向前追溯到1978年HEMON在Geophysical Prospecting上发表的文章^[26]。

逆时偏移直接求解波动方程(1)。与单程波方程相比, 波动算子

$ L = \frac{1}{{{v^2}\left( \mathit{\boldsymbol{x}} \right)}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - \Delta $

(25)

是一个无奇性的偏微分算子, 直接差分求解即可以产生出复杂丰富的波动现象, 包括了反射波、透射波、散射波、折射波、回转波和棱柱波等, 而且没有传播倾角限制。可以认为逆时偏移是目前最为先进的叠前偏移方法。

我们可以将真振幅单程波偏移算法推广到逆时偏移, 从而将这种从物理直观得到的成像方法改造为定量的反演算法。

1) 共炮域逆时偏移算法。

通常的逆时偏移是在炮域实现的。首先从已知炮点坐标x_s正向传播波场p_s:

$ \left( {\frac{1}{{{v^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \delta \left( t \right)\delta \left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) $

(26)

然后从地面检波点出发, 反向传播所采集到的共炮记录D(x_r; t; x_s)。即从最后接收到的信号开始逐渐减少时间来反向求解下面的方程并得到反向波场p_r:

$ \left\{ \begin{array}{l} \left( {\frac{1}{{{v^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = 0\\ {p_{\rm{r}}}\left( {{x_{\rm{r}}},{y_{\rm{r}}},{z_{\rm{r}}} = 0;t} \right) = D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) \end{array} \right. $

(27)

真振幅炮域偏移道集可以通过如下的反褶积成像条件得到^[27]:

$ R\left( {\mathit{\boldsymbol{x}};{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \frac{1}{{\rm{ \mathsf{ π} }}}\iint {\frac{{{{\hat p}_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{{{\hat p}_{\rm{s}}}\left( {\mathit{\boldsymbol{x}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}{\rm{d}}\omega } $

(28)

2) 共反射角域逆时偏移算法。

此时反向波场p_r的方程(27)保持不变, 正向波场p_s的方程(26)的右端震源需要作如下修改^[27]:

$ \left( {\frac{1}{{{v^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \chi \left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) $

(29)

其中χ定义为:

$ \begin{array}{*{20}{c}} {\chi \left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \frac{1}{{{{\left( {2{\rm{ \mathsf{ π} }}} \right)}^3}}}\iiint {\frac{1}{{{v_{\rm{s}}}}}\sqrt {1 - \frac{{v_{\rm{s}}^2}}{{{\omega ^2}}}\left( {k_x^2 + k_y^2} \right)} } \cdot } \\ {{{\rm{e}}^{ - {\rm{i}}\left[ {{k_x}\left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) + {k_y}\left( {\mathit{\boldsymbol{y}} - {\mathit{\boldsymbol{y}}_{\rm{s}}}} \right)} \right]}}{{\rm{e}}^{ - {\rm{i}}\omega t}}{\rm{d}}\omega {\rm{d}}{k_x}{\rm{d}}{k_y}} \end{array} $

(30)

为了得到地下反射角道集, 需要如下的关系式:

$ \left\{ \begin{array}{l} \cos \theta = \frac{{\left( {{\mathit{\boldsymbol{k}}_{\rm{s}}} + {\mathit{\boldsymbol{k}}_{\rm{r}}}} \right) \cdot {\mathit{\boldsymbol{k}}_{\rm{r}}}}}{{\left| {{\mathit{\boldsymbol{k}}_{\rm{s}}} + {\mathit{\boldsymbol{k}}_{\rm{r}}}} \right|\left| {{\mathit{\boldsymbol{k}}_{\rm{r}}}} \right|}}\\ \cos \varphi = \frac{{\left( {{\mathit{\boldsymbol{k}}_{\rm{s}}} \times {\mathit{\boldsymbol{k}}_{\rm{r}}}} \right) \cdot \left[ {{\mathit{\boldsymbol{n}}_x} \times \left( {{\mathit{\boldsymbol{k}}_{\rm{s}}} + {\mathit{\boldsymbol{k}}_{\rm{r}}}} \right)} \right]}}{{\left| {{\mathit{\boldsymbol{k}}_{\rm{s}}} \times {\mathit{\boldsymbol{k}}_{\rm{r}}}} \right|\left| {{\mathit{\boldsymbol{n}}_x} \times \left( {{\mathit{\boldsymbol{k}}_{\rm{s}}} + {\mathit{\boldsymbol{k}}_{\rm{r}}}} \right)} \right|}} \end{array} \right. $

(31)

式中:n_x=(1, 0, 0);k_s和k_r分别是在每个成像点处正向波场和反向波场的波数向量。可以采用如下的互相关型成像条件输出反射角道集^[28]:

$ \begin{gathered} R\left( {\mathit{\boldsymbol{x}};\theta ,\varphi } \right) = 32{{\rm{ \mathsf{ π} }}^2}\frac{{v\left( \mathit{\boldsymbol{x}} \right)}}{{\sin \theta }}\iint {\iiint {{p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){p_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) \cdot }} \hfill \\ \;\;\;\;\;\;\delta \left( {\theta - \theta '} \right)\delta \left( {\varphi - \varphi '} \right){\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}{\rm{d}}t{\rm{d}}\theta '{\rm{d}}\varphi ' \end{gathered} $

(32)

CLAERBOUT^[29]提出了反褶积型和互相关型两种成像条件, 并且指出两种成像条件可以得到一致的构造成像, 不同的偏移振幅。35年后, 我们从真振幅偏移的角度出发, 证明了两种成像条件经过改造后可以分别产生炮域和反射角域的保振幅偏移结果。

3 叠前深度偏移的实践与发展

偏移成像是一项实用生产技术。高精度、高效率的算法和软件编制技术是提高生产质量、节约工作成本的重要保证。叠前逆时偏移于2005年前后应用于实际地震资料处理并迅速得到推广, 在精度、效率和实用性等方面都取得了重要突破。现阶段全波场逆时偏移和Kirchhoff积分偏移各占半壁江山, 而射线束偏移和单程波偏移一般只起到辅助作用。逆时偏移是深海复杂构造勘探中深度成像的主要方法之一。而Kirchhoff偏移是一种灵活的成像方法, 在背景速度比较简单的情况下可以快速地提供高精度道集和成像结果, 在陆地和浅海勘探中仍然发挥着重要作用。理论上讲, 逆时偏移是一种更理想的成像技术, 或许在不远的将来它会得到广泛的应用, 这有待于勘探技术的演化和计算能力的提高。本节侧重介绍过去十多年来发展起来的一些行之有效的逆时偏移理论和算法, 并且探讨一些至今尚未完全解决的问题。

3.1 叠前道集的输出

叠前深度偏移需要输出道集来做速度分析、偏移质量监控、成像后处理以及AVO分析。Kirchhoff偏移算法能够分解为单地震道的输入和输出, 可以很自然地按照地面坐标标定输出道集, 或者附加一点计算量在成像过程中按照射线参数输出反射角和成像倾角道集^[30]。

对于逆时偏移来说, 每次偏移一组共炮记录, 最自然的输出方式是共炮道集。但是一个共炮剖面成像范围很窄, 同相轴在炮道集上表现得很不规则, 而且炮域偏移有严重多解性。XU等^[31]论证了反射角域的共成像点道集是复杂构造成像最为有效的工具:能够大幅度降低偏移解的非唯一性, 从而在每个共反射角剖面上减少偏移假象; 反射角域的真振幅偏移对应的互相关成像条件, 避免了由反褶积成像条件所带来的除法不稳定性; 由反射角道集信息与偏移剖面上的反射面倾角信息相结合, 可以直接确定层析成像所需的入射射线和反射射线, 提高了速度建模效率; 保振幅的反射角道集能够直接用于AVA(AVO)反演和分析; 这种输出方法不依赖于采集观测系统, 适用于陆地、海上拖缆、OBN、VSP以及多次波成像等。

但是, 从逆时偏移中生成反射角道集的效率问题一直没有得到完满解决。严格地由公式(31)来做道集转化非常昂贵, 计算成本达到N⁹量级, 也就是说这个过程的软件实现至少要包含9个嵌套在一起的时间或空间循环。而一个完整的叠前逆时偏移的计算量也只有N⁶量级。常用的Poynting向量法^[32]为了降低计算成本需要假设波场在某个时间和某个地点只有几个主方向。当地下构造复杂时, 各种波场经常重叠在一起, 很难区分开来。这时输出道集的成像质量和稳定性都受到影响。当然, 此类方法在使用中也不断得到改进, 例如文献[33]-[35]中所介绍的方法。另一方面, 虽然真振幅共反射角道集采用了稳定的互相关成像条件, 但是公式(32)中还需要除以sinθ, 所以在0°反射角附近仍然具有奇性。避免这种奇性的一种可行方案是采用SAVA等提出的地下炮检距道集概念^[36]。例如:

$ \begin{gathered} I\left( {\mathit{\boldsymbol{x}};\mathit{\boldsymbol{h}}} \right) = 8\iint {v\left( \mathit{\boldsymbol{x}} \right){p_F}\left( {x - {h_x},y - {h_y},z;t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) \cdot } \hfill \\ \;\;\;\;\;\;{p_B}\left( {x + {h_x},y + {h_y},z;t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){\rm{d}}t{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}} \hfill \\ \end{gathered} $

(33)

利用关系^[36-37]

$ \left\{ \begin{array}{l} \tan \theta = \sqrt {\frac{{k_z^2\left( {k_{hx}^2 + k_{hy}^2} \right) + {{\left( {{k_x}{k_{hx}} + {h_y}{k_{hy}}} \right)}^2}}}{{k_z^2\left( {k_x^2 + k_y^2 + k_z^2} \right)}}} \\ \cos \varphi = \\ \;\;\;\;\;\frac{{{k_{hx}}\sqrt {k_z^2\left( {k_x^2 + k_y^2 + k_z^2} \right)} }}{{\sqrt {k_z^2 + k_y^2} \sqrt {k_z^2\left( {k_{hx}^2 + k_{hy}^2} \right){{\left( {{k_x}{k_{hx}} + {k_y}{k_{hy}}} \right)}^2}} }} \end{array} \right. $

(34)

将地下炮检距道集变换为反射角道集:

$ \hat R\left( {\mathit{\boldsymbol{k}};\theta ,\varphi } \right) = \dot {\hat I}\left[ {\mathit{\boldsymbol{k}}\left( {\theta ,\varphi } \right);{\mathit{\boldsymbol{k}}_h}\left( {\theta ,\varphi } \right)} \right]\frac{1}{{\sin \theta }}\left| {\frac{{\partial \left( {{k_{hx}},{k_{hy}}} \right)}}{{\partial \left( {\theta ,\varphi } \right)}}} \right| $

(35)

可以证明这样所得到的共成像点道集和真振幅偏移公式(32)等价。因为

$ \begin{array}{l} \frac{1}{{\sin \theta }}\left| {\frac{{\partial \left( {{k_{hx}},{k_{hy}}} \right)}}{{\partial \left( {\theta ,\varphi } \right)}}} \right| = \frac{1}{{\sin \theta }}\frac{{\sin \theta \left| {{k_z}} \right|\left| \mathit{\boldsymbol{k}} \right|}}{{{{\cos }^3}\theta }}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = \frac{{\left| {{k_z}} \right|\left| \mathit{\boldsymbol{k}} \right|}}{{{{\cos }^3}\theta }} \end{array} $

(36)

公式(32)中的奇性项1/sinθ可以自然地吸收到反射角道集转换的过程中。该方法实现成本也比较高, 主要是生成地下炮检距道集的计算量是N⁸量级。

关于反射角道集, MONTEL等^[38]发现了一个有趣的现象:当偏移速度不正确的时候, 由下面3种方法所得到的反射角道集上, 共反射点同相轴的曲率不一致。

1) 由Kirchhoff偏移单道输入, 在成像时按照射线参数分解输出的角道集^[30-31];

2) 由共炮域的逆时偏移或者单程波偏移按照Poynting向量分解所产生的角道集^{[28, 33-35]};

3) 由地下炮检距道集转化生成的反射角道集^[36-37]。

因此, 在用叠前偏移道集进行层析成像速度分析之前还需要了解道集是如何生成的。不同方法所得到的共成像点道集也很有可能表现出不同的曲率特征。必须知道输入数据的来龙去脉, 采取相对应的射线追踪方法求取稳相射线, 这样才能保证快速、正确地建立速度模型。

早在单程波偏移时期, ZHOU等^[39]提出了一种通过分解共炮记录输出地面炮检距道集的波动方程偏移方法。此方法的优点是可以直接使用针对Kirchhoff偏移所开发的层析成像软件来进行速度分析, 因而在工业界受到欢迎。需要指出的是, 无论通过怎样的波场组合, 都无法从逆时偏移中得到真振幅三维炮检距道集。

3.2 各向异性

实际勘探中, 完全各向同性的地质构造几乎不存在。使用的速度网格一般是十几米至几百米, 分辨率比较低, 所反演的只能是较大尺度的等效速度。理论上可以证明, 即使是水平均匀介质, 如果采用较大尺度的速度平均(平滑), 同时尽量保持各个方向的走时, 都不可避免地需要引进各向异性的速度概念^[40]。特别是随着采集技术的发展和勘探工作的深化, 多方位角、宽方位角以及全方位角数据越来越普遍, 使得我们能够直接从数据上观测到地震波沿不同方向传播的速度差异。因此, 各向异性建模与井控数据联合反演是当前叠前深度偏移的流行趋势。

最简单的各向异性是VTI(vertical transverse isotropy)。VTI介质以垂直方向为对称轴, 沿水平平面的波传播仍然保持各向同性, 而沿垂直方向的波传播速度(v_n)则通常较慢一些。有很多文献描述了VTI的拟声波方程。其推导过程一般是在弹性波方程的基础上, 将一个S波速度置零而得到P波的频散关系式:

$ \begin{array}{l} \frac{{{v^2}\left( \psi \right)}}{{v_n^2}} = \frac{1}{2} + \varepsilon {\sin ^2}\psi + \\ \;\;\;\;\;\;\;\;\;\;\frac{{\sqrt {{{\left( {1 + 2\varepsilon {{\sin }^2}\psi } \right)}^2} - 8\left( {\varepsilon - \delta } \right){{\sin }^2}\psi {{\cos }^2}\psi } }}{2} \end{array} $

(37)

继而写出满足以上频散关系的声波方程组^[41-42]。(37)式中v(ψ)表示与对称轴方向夹角为ψ的方向上的相速度, v_n≡v(0), ε和δ是引进的Thomsen各向异性参数^[43]。这种方法得到的所谓各向异性拟声波方程虽然简化了偏移计算, 但是只能保证走时的正确性, 方程的振幅大多是错误的, 甚至于方程本身不能保守任何能量形式, 从而造成相应的逆时偏移计算结果的不稳定性。

严格的各向异性方程是在弹性波框架下给出的。下面的VTI方程组^[44]是稳定的, 而且可以证明它与对应的弹性波方程组等价:

$ \begin{array}{l} \frac{1}{{v_n^2}}\frac{{{\partial ^2}}}{{\partial {t^2}}}\left( {\begin{array}{*{20}{c}} p\\ r \end{array}} \right) = \\ \;\;\;\;\left( {\begin{array}{*{20}{c}} {1 + 2\varepsilon }&{\sqrt {1 + 2\delta } }\\ {\sqrt {1 + 2\delta } }&1 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {\frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{{{\partial ^2}}}{{\partial {y^2}}}}&0\\ 0&{\frac{{{\partial ^2}}}{{\partial {z^2}}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} p\\ r \end{array}} \right) \end{array} $

(38)

VTI模型实际上将地质构造简化为理想的水平层状介质。更一般的各向异性是TTI(tilted transverse isotropy), 速度对称轴可以随地质构造处处变化。DUVENECK等^[45]和MACESANU^[46]分别推导出与TTI弹性波方程等价的声波方程形式:

$ \begin{array}{l} \frac{1}{{v_n^2}}\frac{{{\partial ^2}}}{{\partial {t^2}}}\left( {\begin{array}{*{20}{c}} p\\ r \end{array}} \right) = \\ \;\;\;\;\left( {\begin{array}{*{20}{c}} {1 + 2\varepsilon }&{\sqrt {1 + 2\delta } }\\ {\sqrt {1 + 2\delta } }&1 \end{array}} \right)\sum\limits_{i = 1}^3 {\left( {\begin{array}{*{20}{c}} {\mathit{\boldsymbol{A}}_i^{\rm{T}}{\mathit{\boldsymbol{A}}_i}}&{\mathit{\boldsymbol{A}}_i^{\rm{T}}{\mathit{\boldsymbol{B}}_i}}\\ {\mathit{\boldsymbol{B}}_i^{\rm{T}}{\mathit{\boldsymbol{A}}_i}}&{\mathit{\boldsymbol{B}}_i^{\rm{T}}{\mathit{\boldsymbol{B}}_i}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} p\\ r \end{array}} \right)} \end{array} $

(39)

其中, 矩阵A_i和B_i(i=1, 2, 3)由下式定义

$ \begin{array}{l} {\mathit{\boldsymbol{A}}_i} = \left( \begin{array}{l} \frac{\partial }{{\partial {x_1}}}\left( {{R_{1i}}{R_{11}} + {R_{2i}}{R_{21}}} \right)\\ \frac{\partial }{{\partial {x_2}}}\left( {{R_{1i}}{R_{12}} + {R_{2i}}{R_{22}}} \right)\\ \frac{\partial }{{\partial {x_3}}}\left( {{R_{1i}}{R_{13}} + {R_{2i}}{R_{23}}} \right) \end{array} \right)\\ {\mathit{\boldsymbol{B}}_i} = \left( {\begin{array}{*{20}{c}} {\frac{\partial }{{\partial {x_1}}}{R_{3i}}{R_{31}}}\\ {\frac{\partial }{{\partial {x_2}}}{R_{3i}}{R_{32}}}\\ {\frac{\partial }{{\partial {x_3}}}{R_{3i}}{R_{33}}} \end{array}} \right) \end{array} $

(40)

R_ij是如下旋转矩阵的分量:

$ \left( \begin{array}{l} {x'}\\ {y'}\\ {z'} \end{array} \right) = \left( {\begin{array}{*{20}{c}} {\cos \varphi \cos \theta }&{\sin \varphi \cos \theta }&{ - \sin \theta }\\ { - \sin \varphi }&{\cos \varphi }&0\\ {\cos \varphi \sin \theta }&{\sin \theta \sin \varphi }&{\cos \theta } \end{array}} \right)\left( {\begin{array}{*{20}{c}} x\\ y\\ z \end{array}} \right) $

(41)

公式(39)计算起来比较复杂, ZHANG等提出对方程(38)做坐标旋转, 并且保证算子的自共轭性质。所得到的TTI方程组:

$ \begin{array}{l} \frac{1}{{v_n^2}}\frac{{{\partial ^2}}}{{\partial {t^2}}}\left( {\begin{array}{*{20}{c}} p\\ r \end{array}} \right) = \\ \;\;\;\;\left( {\begin{array}{*{20}{c}} {1 + 2\varepsilon }&{\sqrt {1 + 2\delta } }\\ {\sqrt {1 + 2\delta } }&1 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {\partial _{x'}^T{\partial _{x'}} + \partial _{y'}^T{\partial _{y'}}}&0\\ 0&{\partial _{z'}^T{\partial _{z'}}} \end{array}} \right)\left( {\begin{array}{*{20}{c}} p\\ r \end{array}} \right) \end{array} $

(42)

是保能量的, 因而是稳定的^[47]。FLETCHER等对文献[48]中的方程做了推广, 通过引入非零的S波速度来降低方程的不稳定性^[49], 实现了TTI逆时偏移的生产应用。

2008年底以来, TTI各向异性逆时偏移算法在实际资料处理中得到了应用。虽然TTI的偏移计算成本和模型反演的不确定性都大幅度增加, 但应用起来更加灵活, 特别是在墨西哥湾地区的地震资料处理中取得了良好的效果^[50], 并且迅速推广。至此, TTI建模和成像几乎成为叠前深度偏移处理的基本要求。

如果地质构造沿某个方向沉积的同时又受到横向拉伸, 在与沉积方向垂直的平面上, 地震波传播的各向同性也会受到影响, 这就可能形成旋转坐标下的倾斜正交各向异性(Tilted Orthorhombic Anisotropy)。这种情形下, 上述各类拟声波方程组都可以从TTI推广到倾斜正交各向异性的正演或偏移方程^[51-52]。例如公式(42)可以推广为^[51]:

$ \frac{1}{{v_n^2}}\frac{{{\partial ^2}}}{{\partial {t^2}}}\left( {\begin{array}{*{20}{c}} \begin{array}{l} p\\ q \end{array}\\ r \end{array}} \right) = \mathit{\boldsymbol{N}}{\mathit{\boldsymbol{D}}^{\rm{T}}}\mathit{\boldsymbol{D}}\left( {\begin{array}{*{20}{c}} \begin{array}{l} p\\ q \end{array}\\ r \end{array}} \right) $

(43)

其中,

$ \begin{array}{l} \mathit{\boldsymbol{N}} = \\ \left( {\begin{array}{*{20}{c}} {1 + 2{\varepsilon _2}}&{\left( {1 + 2{\varepsilon _2}} \right)\sqrt {1 + 2{\delta _3}} }&{\sqrt {1 + 2{\delta _2}} }\\ {\left( {1 + 2{\varepsilon _2}} \right)\sqrt {1 + 2{\delta _3}} }&{1 + 2{\varepsilon _1}}&{\sqrt {1 + 2{\delta _1}} }\\ {\sqrt {1 + 2{\delta _2}} }&{\sqrt {1 + 2{\delta _1}} }&1 \end{array}} \right) \end{array} $

(44)

$ \mathit{\boldsymbol{D}} = \left( {\begin{array}{*{20}{c}} {{\partial _{x'}}}&0&0\\ 0&{{\partial _{y'}}}&0\\ 0&0&{{\partial _{z'}}} \end{array}} \right) $

(45)

各向异性地震资料处理的奠基人Leon Thomsen博士认为, 现实世界中的各向异性至少是倾斜正交类型, 所以此类偏移有着更加广阔的应用前景。特别是随着高密度全方位角地震采集的广泛实施, 倾斜正交偏移比TTI偏移更容易在地下角度域将各个方位角的共成像点道集拉平, 因而可以更充分地利用全方位角采集的数据信息以提高成像质量和精度^[53]。

前面所述的各种方程都是从弹性波方程出发, 推导出相应的各向异性拟声波方程表达式。以TTI方程为例, 它们在应用中有如下缺陷:

1) 尽管在拟声波方程的推导过程中将S波速度置零, 但是所得到的方程组((39)式或(42)式)的解除去P波之外仍然有一个S波的增根。因而在正演模拟和逆时偏移过程中能够看到拟S波的噪声干扰^[54];

2) 由于前述的拟S波的传播速度与$\sqrt {\varepsilon - \delta }$成正比, 因此, 为了保证方程的稳定性就需要ε≥δ的附加条件;

3) 将TTI方程组(42)与各向同性的方程(1)相比较, 我们不但要多求解一个方程, 而且要附加很多计算量来求解各种交叉导数。直接实现TTI偏移的计算量大约是各向同性偏移的3倍, 增加了生产成本。

针对以上问题, XU等^[55]另辟蹊径, 提出了只有纯P波的TTI拟声波方程:

$ \left[ {\frac{{{\partial ^2}}}{{\partial {t^2}}} - \nabla \cdot \left( {v_n^2S\nabla } \right)} \right]p\left( {\mathit{\boldsymbol{x}};t} \right) = 0 $

(46)

其中的算子

$ \begin{array}{l} S = \left\{ {1 + 2\varepsilon \left( {n_x^2 + n_y^2} \right) + } \right.\\ \left. {\sqrt {{{\left[ {1 + 2\varepsilon \left( {n_x^2 + n_y^2} \right)} \right]}^2} - 8\left( {\varepsilon - \delta } \right)\left( {n_x^2 + n_y^2} \right)n_z^2} } \right\}/2 \end{array} $

(47)

由波场的传播方向

$ \mathit{\boldsymbol{n}} = \left( {{n_x},{n_y},{n_z}} \right) = \frac{\mathit{\boldsymbol{k}}}{{\left| \mathit{\boldsymbol{k}} \right|}} $

(48)

来定义, 可以按照下式

$ \mathit{\boldsymbol{n}} = \frac{{\nabla p}}{{\left| {\nabla p} \right|}} $

(49)

从波场中推算出来。为了增加计算精度, 他们又进一步提出方程

$ \begin{array}{*{20}{c}} {\left\{ {\frac{1}{{v_n^2}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - \left( {1 + \Delta {S_e}} \right)\left[ {\frac{{{\partial ^2}}}{{\partial {z^2}}} - \left( {1 + 2\varepsilon } \right) \cdot } \right.} \right.}\\ {\left. {\left. {\left( {\frac{{{\partial ^2}}}{{\partial {x^2}}} + \frac{{{\partial ^2}}}{{\partial {y^2}}}} \right)} \right]} \right\}p\left( {\mathit{\boldsymbol{x}};t} \right) = 0} \end{array} $

(50)

并且引入了椭圆型修正量^[56]

$ \Delta {S_e} = \frac{1}{2}\left\{ {\sqrt {1 - \frac{{8\left( {\varepsilon - \delta } \right)\left( {n_x^2 + n_y^2} \right)n_z^2}}{{{{\left[ {1 + 2\varepsilon \left( {n_x^2 + n_y^2} \right)} \right]}^2}}}} - 1} \right\} $

(51)

对于Kirchhoff偏移, 从各向同性向各向异性的过渡较为直接。只需将射线方程作相应的修改^[57], 几乎不附加任何偏移计算量。参考文献[58]给出了一个较早的陆上TTI Kirchhoff深度偏移结果的例子。

3.3 非弹性衰减

地下介质并非是一个能量守恒的纯弹性介质, 从实际接收的反射地震波信号中总能够观测到吸收衰减效应, 在气云区附近甚至能够看到由吸收衰减而导致的明显的弱照明效应。这种波在传播中的非弹性吸收衰减效应不仅逐渐削弱高频能量, 而且造成走时随频率改变的频散效应, 在影响偏移能量聚焦的同时也扭曲了成像界面上的子波波形和相位。需要在叠前偏移中加以校正。

根据KJARTANSSON等^[59]的研究成果, 线性粘性效应可以用下面的复频散关系来描述

$ \frac{{{\omega ^2}}}{{{v^2}}}{\left( {\frac{{{\rm{i}}\omega }}{{{\omega _0}}}} \right)^{ - 2\gamma }} = {\mathit{\boldsymbol{k}}^2} $

(52)

其中,

$ \gamma \left( \mathit{\boldsymbol{x}} \right) = \frac{1}{{\rm{ \mathsf{ π} }}}{\tan ^{ - 1}}\left[ {\frac{1}{{Q\left( \mathit{\boldsymbol{x}} \right)}}} \right] \approx \frac{1}{{{\rm{ \mathsf{ π} }}Q\left( \mathit{\boldsymbol{x}} \right)}} $

(53)

$ v = {v_0}\cos \left( {\frac{{{\rm{ \mathsf{ π} }}\gamma }}{2}} \right) $

(54)

式中:Q是空变品质因子; v₀是相对于参考频率ω₀的波速。

公式(52)是一个随频率变化的频散关系, 在频率域写出相应的微分方程相对比较容易, 例如, 其频域单程波形式较早就有人研究并且应用于生产^[60-61]。XIE等^[62]介绍了如何在Kirchhoff偏移中实现吸收衰减补偿的方法:先利用层析成像来反演Q模型, 再采用具有Q补偿功能的Kirchhoff偏移进行深度域成像。这两种技术的结合在亚太地区的浅水勘探中获得了很多成功^[63]。

与单程波偏移通常在频率域求解波动方程不同, 通行的逆时偏移都是在时间域进行波场延拓。这就需要在时间域直接写出或间接逼近公式(52)所对应的拟微分算子。传统的做法是在时间域逼近算子(iω)^－2γ[64]。例如, 可以采用分式逼近

$ {\left( {{\rm{i}}\omega } \right)^{ - 2\gamma \left( x \right)}} \approx 1 + \sum\limits_{l = 1}^L {\frac{{{c_l}\left( \mathit{\boldsymbol{x}} \right)}}{{{\rm{i}}\omega + {\omega _l}\left( \mathit{\boldsymbol{x}} \right)}}} $

(55)

于是由频散关系式(52)就得到近似的带有粘性的声波方程组

$ \left\{ \begin{array}{l} \left[ {\frac{1}{{v_0^2}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - \frac{{{{\cos }^2}\left( {\frac{{{\rm{ \mathsf{ π} }}\gamma }}{2}} \right)}}{{\omega _0^{2\gamma }}}\Delta } \right]p\left( {\mathit{\boldsymbol{x}};t} \right) + \frac{1}{{v_0^2}}\sum\limits_{l = 1}^L {{e_l}\left( {\mathit{\boldsymbol{x}};t} \right)} = 0\\ \left( {{\omega _l} + \frac{\partial }{{\partial t}}} \right){e_l}\left( {\mathit{\boldsymbol{x}};t} \right) = {c_l}p\left( {\mathit{\boldsymbol{x}};t} \right) \end{array} \right. $

(56)

方程(56)引入了多个记忆函数e_l(x; t), 增加了一些计算量。为了与已有的逆时偏移应用程序兼容, ZHANG等^[65]对频散关系式((52)式)做了调整, 近似写为:

$ {\omega ^2} - \frac{{{\rm{i}}\omega }}{Q}{\left( {\frac{{{v_0}\left| \mathit{\boldsymbol{k}} \right|}}{{\omega _0^\gamma }}} \right)^{\frac{1}{{1 - \gamma }}}} - {\left( {\frac{{{v_0}\left| \mathit{\boldsymbol{k}} \right|}}{{\omega _0^\gamma }}} \right)^{\frac{2}{{1 - \gamma }}}} = 0 $

(57)

并且提出了如下的粘性声波方程:

$ \left( {\frac{{{\partial ^2}}}{{\partial {t^2}}} + \frac{\mathit{\Phi }}{Q}\frac{\partial }{{\partial t}} + {\mathit{\Phi }^2}} \right)p = 0 $

(58)

其中, Φ为拟微分算子, 可以定义为:

$ \mathit{\Phi } = {\left( {\frac{{\sqrt { - v_0^2\Delta } }}{{\omega _0^\gamma }}} \right)^{\frac{1}{{1 - \gamma }}}} $

(59)

方程(58)有比较清晰的物理意义。括号中的第2个算子$[(\mathit{\Phi /Q})(\partial /\partial t)]$是衰减项, 第3个算子(Φ²)是频散项。ZHU等^[66]采取了类似的方式, 将频散关系式(52)近似改写为:

$ \begin{array}{*{20}{c}} {\frac{{{\omega ^2}}}{{{c^2}}} - \left( {{\rm{i}}\omega } \right)v_0^{ - 1 + \gamma }\omega _0^{ - \gamma }\sin \left( {\pi \gamma } \right){{\left| \mathit{\boldsymbol{k}} \right|}^{1 + \gamma }} - }\\ {v_0^{2\gamma }\omega _0^{ - 2\gamma }\cos \left( {\pi \gamma } \right){{\left| \mathit{\boldsymbol{k}} \right|}^{2\left( {1 + \gamma } \right)}} = 0} \end{array} $

(60)

并且写出粘性声波方程:

$ \left[ {\frac{1}{{{v^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - \tau \frac{\partial }{{\partial t}}{{\left( { - \Delta } \right)}^{\frac{{1 + \gamma }}{2}}} - \eta {{\left( { - \Delta } \right)}^{\left( {1 + \gamma } \right)}}} \right]p = 0 $

(61)

其中,

$ \left\{ \begin{array}{l} \eta \left( \mathit{\boldsymbol{x}} \right) = - v_0^{2\gamma }\omega _0^{ - 2\gamma }\cos \left( {{\rm{ \mathsf{ π} }}\gamma } \right)\\ x\left( \mathit{\boldsymbol{x}} \right) = - v_0^{ - 1 + \gamma }\omega _0^{ - \gamma }\sin \left( {{\rm{ \mathsf{ π} }}\gamma } \right) \end{array} \right. $

(62)

方程(61)和方程(58)非常相近。方程(61)括号中的第2个算子$\tau (\partial /\partial t){( - \Delta )^{(1 + \gamma )/2}}$是衰减项, 第3个算子η(－Δ)^1+γ是频散项。同时, ZHU等^[66]还给出了等价的方程组形式, 例如:

$ \left\{ \begin{array}{l} \frac{{\partial \mathit{\boldsymbol{u}}\left( {\mathit{\boldsymbol{x}};t} \right)}}{{\partial t}} = \nabla p\left( {\mathit{\boldsymbol{x}};t} \right)\\ \frac{{\partial \varepsilon \left( {\mathit{\boldsymbol{x}};t} \right)}}{{\partial t}} = \nabla \cdot \mathit{\boldsymbol{u}}\left( {\mathit{\boldsymbol{x}};t} \right)\\ p\left( {\mathit{\boldsymbol{x}};t} \right) = {v^2}\left[ {\eta \left( \mathit{\boldsymbol{x}} \right){{\left( { - \Delta } \right)}^\gamma } + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\left. {\tau \left( \mathit{\boldsymbol{x}} \right)\frac{\partial }{{\partial t}}{{\left( { - \Delta } \right)}^{\frac{{\gamma - 1}}{2}}}} \right]\varepsilon \left( {\mathit{\boldsymbol{x}};t} \right) \end{array} \right. $

(63)

式中:u是质点的速度向量; p和ε分别是压力场和应变场。

由于深度偏移生产中至少需要TTI模型, 因此, 如何将粘性声波方程推广到TTI介质是一个重要研究课题。XIE等^[67]将粘性方程组(56)和TTI拟声波方程组(42)相结合, 提出如下的TTI粘性声波方程:

$ \left\{ \begin{array}{l} \frac{{{\partial ^2}}}{{\partial {t^2}}}\left( \begin{array}{l} p\\ r \end{array} \right) + \sum\limits_{l = 1}^L {\left( \begin{array}{l} {e_l}\\ {d_l} \end{array} \right)} = \frac{{v_0^2{{\cos }^2}\left( {\frac{{{\rm{ \mathsf{ π} }}\gamma }}{2}} \right)}}{{\omega _0^{2\gamma }}} \cdot \\ \left( {\begin{array}{*{20}{c}} {1 + 2\varepsilon }&{\sqrt {1 + 2\delta } }\\ {\sqrt {1 + 2\delta } }&1 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {\partial _{x'}^T{\partial _{x'}} + \partial _{y'}^T{\partial _{y'}}}&0\\ 0&{\partial _{z'}^T{\partial _{z'}}} \end{array}} \right)\left( \begin{array}{l} p\\ r \end{array} \right)\\ \left( {{\omega _l} + \frac{\partial }{{\partial t}}} \right)\left( \begin{array}{l} {e_l}\\ {d_l} \end{array} \right) = {c_l}\left( \begin{array}{l} p\\ r \end{array} \right) \end{array} \right. $

(64)

而ZHU等^[68]的方法((63)式)则可以推广到弹性波的表达式。按照他们的思路, 可以将公式(42)改写为带粘性的TTI拟声波方程组:

$ \begin{array}{l} \frac{1}{{v_0^2}}\frac{{{\partial ^2}}}{{\partial {t^2}}}\left( \begin{array}{l} p\\ r \end{array} \right) = \left[ {\eta {{\left( { - \mathit{\Omega }} \right)}^\gamma } + \tau \frac{\partial }{{\partial t}}{{\left( { - \mathit{\Omega }} \right)}^{\frac{{\gamma - 1}}{2}}}} \right] \cdot \\ \left( {\begin{array}{*{20}{c}} {1 + 2\varepsilon }&{\sqrt {1 + 2\delta } }\\ {\sqrt {1 + 2\delta } }&1 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {\partial _{x'}^T{\partial _{x'}} + \partial _{y'}^T{\partial _{y'}}}&0\\ 0&{\partial _{z'}^T{\partial _{z'}}} \end{array}} \right)\left( \begin{array}{l} p\\ r \end{array} \right) \end{array} $

(65)

其中, Ω是一个拟微分算子, 具有特征表达式:

$ \begin{array}{l} \kappa = \left\{ {\left( {1 + 2\varepsilon } \right)\left( {k_x^2 + k_y^2} \right) + k_z^2 + } \right.\\ \left. {\sqrt {{{\left[ {\left( {1 + 2\varepsilon } \right)\left( {k_x^2 + k_y^2} \right) + k_z^2} \right]}^2} - 8\left( {\varepsilon - \delta } \right)\left( {k_x^2 + k_y^2} \right)k_z^2} } \right\}\\ /2 \end{array} $

(66)

使用带粘性的拟声波方程做逆时偏移可以在成像过程中自动校正由吸收作用所导致的相位频散和能量衰减。但是在偏移的反向传播过程中波场的能量会不断增加, 从而导致计算不稳定。XIE等^[67]提出了一种替代的能量补偿方法。将其思路应用于方程(58), 下面以从检波点出发的反向波场为例。

首先, 求解一个沿着时间反向传播、具有等效能量衰减作用的方程

$ \left\{ \begin{array}{l} \left( {\frac{{{\partial ^2}}}{{\partial {t^2}}} - \frac{\mathit{\Phi }}{Q}\frac{\partial }{{\partial t}} + {\mathit{\Phi }^2}} \right){q_r}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = 0\\ {q_r}\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}},{y_r},{z_r} = 0;t} \right) = D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) \end{array} \right. $

(67)

其次, 求解一个反向传播能量不衰减, 同时具有等效频散补偿效应的方程:

$ \left\{ \begin{array}{l} \left( {\frac{{{\partial ^2}}}{{\partial {t^2}}} + {\mathit{\Phi }^2}} \right){r_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = 0\\ {r_r}\left( {{x_{\rm{r}}},{y_r},{z_r} = 0;t} \right) = D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) \end{array} \right. $

(68)

所需要的能量补偿后的反向波场p_r可以由上面两个波场组合得到:

$ {{\hat p}_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = {{\hat q}_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){\left| {\frac{{{{\hat r}_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{{q_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}} \right|^2} $

(69)

这样就将微分方程的不稳定性转化到反褶积当中, 比较易于处理和控制。

由于粘性声波方程引入了空变的非整数阶导数算子(iω)^－2γ(x), 因此, 不可避免地需要数值计算一些比较复杂的拟微分算子, 例如:公式(58)中的Φ, 公式(61)中的(－Δ)^(1+γ)/2和(－Δ)^1+γ。FOMEL等在这方面做了大量研究工作, 如采用低秩(low-rank)方法逼近求解高维复杂的波动算子, 包括各向同性^[69-71]、各向异性^[72]、粘性方程^[73-74]等, 值得关注。

3.4 逆时偏移中的组合震源、鬼波和多次波的处理

前述偏移公式(26)~公式(28)基于以下假设:①激发震源是单点脉冲; ②每次接收到的信号是一个单炮记录; ③震源和检波器都置于地表或海面; ④所偏移的地震波信号是地下一次反射波。

实际生产中这4个假设往往都不成立。以海上拖缆资料为例:

1) 海上震源通常采用气枪组合, 所合成的子波并非点脉冲, 而是有一定频宽的波形, 并且随传播角和方位角变化^[75-76];

2) 为了提高效率、降低成本, 海上采集开始采用多炮同时激发的方式, 即相距不太远的多个气枪在同一反射时窗相继激发, 得到的是几个炮记录的混合剖面^[77-78];

3) 实际资料采集时气枪和检波器都置于水面以下。由于海面的镜像反射作用, 所接收到的地震信号包含了炮点、检波点、以及两个方向混合的鬼波(随传播方向变化), 因此限制了地震资料的频带;

4) 由于水面的自由面反射, 所接收到的地震记录都是一次反射波和多次波的混响。

为了得到高分辨率、高质量的叠前深度偏移结果, 在成像之前需要对地震资料进行预处理。

1) 震源子波校正(Designature):利用近场和远场所记录到的直达波信号估计阵列震源沿着不同方向的综合子波效应^[79], 并且用反褶积去除震源子波在地震记录中的影响^[80];

2) 去混叠效应(Deblending):将多炮同时激发的资料通过去混响的流程恢复为多个单炮记录^[81-83];

3) 去鬼波处理(Deghosting):通过去鬼波处理消除其陷波影响, 展宽资料频带^[84-89];

4) 去多次波(Demultiple):通过预测多次波和自适应减法以及其它辅助手段压制资料中的多次波^[90-93]。

上述各项预处理所需要解决的问题都具有多维性和不适定性, 原则上都不是在单道上能够简单完成的。地震资料质量受采集成本和设备的局限, 只在数据域做处理未必能够得到最优的成像结果。比如, 去掉的一些噪声可能本来就不影响偏移结果, 浪费了生产时间和计算资源; 由于数据中的假频干扰, 某些传播方向上的鬼波算子计算不正确, 去鬼波之后反而影响了偏移效果, 等等。

逆时偏移不仅是偏移成像的一种有效工具, 而且能够比较自然地进行一部分预处理, 从而简化整个叠前深度偏移的生产流程。例如, 如果我们知道气枪阵列的位置(x^l_s, l=1, 2, …, L)和每个气枪的子波(s^l(t), l=1, 2, …, L), 就可以对正向波场的计算公式(26)作相应修改

$ \left( {\frac{1}{{{v^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_s}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \sum\limits_{l = 1}^L {{s^l}\left( t \right)\delta \left( {\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{x}}_s^l} \right)} $

(70)

并且通过反褶积型成像条件((28)式)来自动达到在偏移中去除三维子波效应的目的。

逆时偏移并不限于对单炮记录作偏移, 也可以用于组合炮记录。原则上可以直接偏移多炮同时激发的混响记录。虽然这样得到的成像结果受到交叉成像噪声的干扰^[94], 而且需要联合应用最小二乘偏移等其它方法以提高成像质量, 但这至少是简化多炮同时激发处理流程的一个思路。

ZHANG等^[95]提出适当修改逆时偏移的震源或者地面边界条件, 去鬼波的过程可以在偏移中自动完成。例如, 如果假设炮点和检波点的水下深度分别为z_s和z_r, 可以通过在镜像炮点x_s=(x_s, y_s, -z_s)加入反极性震源来模拟震源鬼波:

$ \begin{array}{l} \left( {\frac{1}{{{v^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_s}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \\ \;\;\;\;\;\;\;\;\;\;\;s\left( t \right)\left[ {\delta \left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) - \delta \left( {\mathit{\boldsymbol{x}} - {{\mathit{\boldsymbol{\bar x}}}_{\rm{s}}}} \right)} \right] \end{array} $

(71)

对于检波点, 可以在反向波场传播之前沿着波场的初始传播方向做去鬼波处理:

$ \left\{ \begin{array}{l} \left( {\frac{1}{{{v^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = 0\\ {{\hat p}_{\rm{r}}}\left( {{x_{\rm{r}}},{y_{\rm{r}}},{z_{\rm{r}}} = 0;\omega } \right) = \frac{{\hat D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\hat G\left( {\omega ;{\alpha _{\rm{r}}}} \right)}} \end{array} \right. $

(72)

式中

$ \hat G\left( {\omega ;{\alpha _{\rm{r}}}} \right) = - 2{\rm{i}}\sin \frac{{\omega \cos \left( {{\alpha _{\rm{r}}}{z_{\rm{r}}}} \right)}}{{{v_0}}} $

(73)

是沿着与垂直方向成α_r角度的传播方向上的鬼波表达式。α_r可以通过下面频率波数域的公式(74)计算得到:

$ \cos {\alpha _{\rm{r}}} = \sqrt {1 - v_0^2\frac{{k_x^2 + k_y^2}}{{{\omega ^2}}}} $

(74)

最后, 利用反褶积型成像条件((29)式)即可得到真振幅偏移结果。上面描述的偏移公式自动去掉了震源效应、震源鬼波和接收点鬼波。而且, 算法保证了这些步骤的逆算子都是三维的, 不受地震资料中假频的干扰。

在应用波动方程直接偏移多次波方面, WHITMORE等^[96]和LU等^[97-98]开展了一系列研究工作, 在浅海和深海勘探中取得了很多实际应用成果。具体的做法是将地震数据分别赋予正向波场和反向波场作为地表边界条件。在成像过程中一次反射和一阶多次波、多次波和高一阶的多次波都能够进行互相关成像。这样不仅将多次波变成有用信号, 而且在特殊情况下(例如:浅海拖缆采集的crossline方向)还弥补了一次波采样稀疏和照明不足的问题。这可以看成是组合炮偏移的一个特例。对于多次波偏移, 可以输出共反射角道集提供速度分析、照明分析和AVO分析^[99]。LU等^[100]提供了一些压制交叉成像噪声干扰的实用方法。

叠前深度偏移在生产实践中还存在很多重要理论问题和算法问题, 例如:假频、数值频散、计算稳定性、分辨率分析、边界条件、起伏地表等问题。围绕这些问题发展出了很多行之有效的方法和手段, 限于篇幅不一一讨论。

4 最小二乘偏移

既然真振幅偏移已经是波动方程正演问题的逆算子了, 为什么还要做最小二乘偏移？这是因为前面讨论的真振幅偏移对数据采集做了理想的假设, 即:规则采样, 无数据假频, 无数据缺失, 全方位角覆盖, 炮检距足够大, 等等。这些条件在实际生产中几乎都不能够满足。最小二乘偏移^[101-102]是直接将偏移当作一个反演过程来求解, 以解决因数据采集、波动传播和成像条件限制所造成的偏移失真、不聚焦、振幅不平衡和噪声干扰等问题。因而, 理论上最小二乘偏移是一种更符合生产实际要求的成像方法。

最小二乘偏移在生产中的应用还处于试验探索阶段。文献中对一些关键步骤的商业软件实现方法描述得并不详细。这里从理论上对该方法进行分析。

常规的最小二乘偏移是定义一个正演的波动算子M, 由给定的偏移结果R(x)(反射系数)来预测观测数据D_syn(x_r; t; x_s), 即

$ \mathit{\boldsymbol{M}}\left( R \right) = {D_{{\rm{syn}}}} $

(75)

然后定义一个目标函数

$ \begin{gathered} E\left( R \right) = \frac{1}{2}{\left\| {{D_{{\rm{syn}}}} - D} \right\|^2} = \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}\iiint {{{\left( {{D_{{\rm{syn}}}} - D} \right)}^2}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}{\rm{d}}t = \min } \hfill \\ \end{gathered} $

(76)

最小二乘偏移的目的就是要找到一个最优的偏移剖面R(x), 使得在目标函数(76)式中, 由计算得到的数据D_syn与实际观测数据D最为接近。这样, 偏移问题被转化为一个最优化问题。

由于逆时偏移还不能以较低的计算成本输出三维角道集, 因此, 通常假定(75)式中的反射系数模型是与反射角无关的全叠加剖面。这时可以定义正演算子如下。

首先, 从已知的炮点坐标x_s传播波场p_s:

$ \left( {\frac{1}{{{v^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_s}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = s\left( t \right){\rm{ \mathsf{ δ} }}\left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) $

(77)

反射波场p_r是由入射波场p_s与反射系数R作用后所产生的次生波场:

$ \left( {\frac{1}{{{v^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_{\rm{r}}}\left( {x;t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = R\frac{\partial }{{\partial t}}{p_s}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) $

(78)

在检波点处得到预测数据

$ {D_{{\rm{syn}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \mathit{\boldsymbol{M}}\left( R \right) = {p_{\rm{r}}}\left( {\mathit{\boldsymbol{x}} = {\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) $

(79)

公式(78)和公式(79)也经常被称作反偏移。求解优化问题((76)式)需要计算变分

$ \frac{{{\rm{ \mathsf{ δ} }}E}}{{{\rm{ \mathsf{ δ} }}R}} = {\mathit{\boldsymbol{M}}^{\rm{T}}}\left( {{D_{{\rm{syn}}}} - D} \right) = {\mathit{\boldsymbol{M}}^{\rm{T}}}{\rm{ \mathsf{ δ} }}D $

(80)

上面的正演共轭算子可以由下面的方程组来计算^[103]

$ \left\{ \begin{gathered} \left( {\frac{1}{{{v^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_s}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = s\left( t \right){\rm{ \mathsf{ δ} }}\left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) \hfill \\ \left( {\frac{1}{{{v^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\; - \int {\frac{{\partial {\rm{ \mathsf{ δ} }}D}}{{\partial t}}\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){\rm{ \mathsf{ δ} }}\left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{\rm{r}}}} \right){\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}} \hfill \\ \frac{{{\rm{ \mathsf{ δ} }}E}}{{{\rm{ \mathsf{ δ} }}R}} = \iint {{p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){p_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){\rm{d}}t{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}} \hfill \\ \end{gathered} \right. $

(81)

(81)式中的第1个方程是沿时间方向正向求解, 第2个方程是沿时间方向反向求解, 第3个方程是互相关成像条件。实际上等价于求解一个逆时偏移。

由于在反演过程中需要进行波场匹配, 3.4中所讨论的组合震源和多炮同时激发得到的数据都可以在最小二乘逆时偏移的框架下直接进行成像, 而且比较自然地解决了不同的炮记录之间交叉成像的噪声干扰^[103]。关于海上鬼波的问题, 也可以通过在正演过程中加入镜像的虚拟震源和虚拟检波器来预测相应的鬼波, 并且在逐次迭代的过程中逐渐消除鬼波^[103-104]。这些都是逆时偏移特有的优点。需要指出的是, 对反偏移公式不能直接强制自由边界条件, 也就是说, 下面的边界条件提法会导致反偏移算法的不稳定性^[105]

$ \left\{ \begin{gathered} \left( {\frac{1}{{{v^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_{\rm{s}}}\left( {x;t;{x_{\rm{s}}}} \right) = f\left( t \right)\delta \left( {x - {x_{\rm{s}}}} \right) \hfill \\ \left( {\frac{1}{{{v^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_{\rm{r}}}\left( {x;t;{x_{\rm{s}}}} \right) = R\frac{\partial }{{\partial t}}{p_{\rm{s}}}\left( {x;t;{x_{\rm{s}}}} \right) \hfill \\ {p_{\rm{s}}}\left( {x,y,z = 0;t;{x_{\rm{s}}}} \right) = 0 \hfill \\ {p_{\rm{r}}}\left( {x,y,z = 0;t;{x_{\rm{s}}}} \right) = 0 \hfill \\ \end{gathered} \right. $

(82)

究其理论, 是因为反偏移不是一个自然的物理过程, 一般而言能量不守恒。所以我们还不能简单地利用最小二乘逆时偏移来同时偏移一次波和多次波。文献[105]提出了一个思路, 通过修改反偏移中的反射系数可以得到保能量的反偏移公式。这需要附加求解一个反演问题, 由叠加剖面(反射系数)得到伪密度函数。

前面介绍的最小二乘偏移在生产应用中还受到很多限制:

首先, 我们并不一定知道震源子波的真实面目, 而且每次震源的响应都会有些变化。地层是一个复杂的粘弹性体, 不可能简单地通过求解声波方程来准确预测反射数据的振幅。这些实际问题使得人工模拟数据和实际数据之间很难直接匹配。据此, 文献[103]把目标函数改为互相关型:

$ \begin{gathered} E\left( R \right) = \hfill \\ - \iiint {\frac{{{D_{{\rm{syn}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\sqrt {\int {{{\left| {{D_{{\rm{syn}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)} \right|}^2}{\rm{d}}t} } \sqrt {\int {{{\left| {D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)} \right|}^2}{\rm{d}}t} } }}} \cdot \hfill \\ {\rm{d}}t{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}} = \min \hfill \\ \end{gathered} $

(83)

文献[106]则将目标函数(76)式修改为

$ E\left( R \right) = \frac{1}{2}{\left\| {{\mathit{\boldsymbol{P}}_f}\left( {{D_{{\rm{syn}}}}} \right) - {\mathit{\boldsymbol{P}}_{\rm{m}}}\left( D \right)} \right\|^2} = \min $

(84)

并且解释P_f是用来匹配数据的滤波算子, 而P_m是对数据所做的预处理。当偏移速度有一定误差时, 甚至于需要借助时间域动态整形(dynamic warping)的手段来匹配数据^[107-108]。

另外, 我们注意到这里定义的正演算子((77)式~(79)式)与波动方程所定义的反射原理((9)式)不等价。为了节省计算量, 忽略了反射系数随反射角变化的问题, 所以这里的最小二乘偏移只是为了提高叠加剖面上构造的聚焦效果, 不具备真振幅成像的理论根据。

最小二乘偏移每次迭代需要做一次正演和一次成像, 至少相当于两次偏移的计算量, 因此, 计算成本比较昂贵。于是, 研究人员寻求了一些替代方案。比较著名的是偏移反褶积方法(migration decon)^[109-113]。对于第一次迭代, 如果我们直接写出泛函(76)式的解

$ R \approx - {\left( {{\mathit{\boldsymbol{M}}^{\rm{T}}}\mathit{\boldsymbol{M}}} \right)^{ - 1}}{\mathit{\boldsymbol{M}}^{\rm{T}}}\left( D \right) $

(85)

式中:－M^T(D)就是常规偏移结果。M^TM是在确定了偏移速度和观测系统后, 随着每个成像点空间变化的聚焦算子, 与数据无关。具体的计算方法是在成像空间的粗网格(x_k, k=1, 2, …, N)上排布空间点列, 作为反射系数, 即$ R(\mathit{\boldsymbol{x}}) = \sum\limits_k {\delta (\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_k})} $。通过正演模拟M计算其在实际观测系统下的接收响应, 再进行偏移(M^T)来计算每个预设点脉冲(δ(x－x_k))的成像效果。这里所得到的像反映了在不完善的观测系统下偏移方法对于每个空间点脉冲的聚焦扭曲程度。近似地表达了算子M^TM在每个空间点上的响应。我们利用这个响应函数对偏移结果逐点做空间反褶积计算((M^TM)^－1), 就能够在一定程度上校正偏移过程中聚焦不充分、振幅不平衡的缺陷。这类方法可以看作是对常规偏移方法的改进。由于缺乏数据匹配流程, 偏移反褶积在加强聚焦效果的过程中无法分辨哪些是真正的构造, 哪些是数据采集的噪声和成像过程中的假象, 所以这种方法并没有真正达到最小二乘偏移的效果。

2014年, SALOMONS等^[114]在EAGE年会上报告了最小二乘Kirchhoff偏移生产软件算法, 该方法比较接近反演理论的要求。具体的目标函数定义为:

$ \begin{gathered} E\left( R \right) = \frac{1}{{{{\left\| D \right\|}^2}}}{\left\| {M\left[ {R\left( {\mathit{\boldsymbol{x}};a} \right)} \right] - D} \right\|^2} + \hfill \\ \;\;\;\;\;\;\;\;\;\;\;{\varepsilon _{{\rm{off}}}}{\left\| {\frac{{\partial R\left( {\mathit{\boldsymbol{x}};a} \right)}}{{\partial a}}} \right\|^2} + {\varepsilon _{{\rm{dip}}}}\left[ {{{\left\| {\frac{{\partial R\left( {\mathit{\boldsymbol{x}};a} \right)}}{{\partial x}}} \right\|}^2} + } \right. \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\left. {{{\left\| {\frac{{\partial R\left( {\mathit{\boldsymbol{x}};a} \right)}}{{\partial y}}} \right\|}^2}} \right] = \min \hfill \\ \end{gathered} $

(86)

从上面的泛函定义中我们看到, 最小二乘偏移的优化过程是在道集上定义的, 其目的是提高共成像点道集的成像质量, 并非只是为了聚焦全叠加剖面。公式(86)中的第1项是对于每个共成像指标(a)剖面都按照正演公式(9)来预测观测数据, 并且和实际数据做匹配; 第2项强制在共成像点上道集的相似性; 第3项强制偏移剖面的横向连续性。SALOMONS等在文中介绍了3个陆上资料处理的实例。由于陆上采集受地形、环境和成本的限制, 数据稀疏而不规则。最小二乘Kirchhoff偏移的结果大幅改善了小炮检距剖面和靠近采集区域边界附近的画弧效应, 消除了很多浅层的采集痕迹, 道集上的成像一致性得到了明显的改进, 而且拓宽了有效频带, 提高了分辨率、连续性和信噪比。

5 全波形反演与叠前深度偏移

全波形反演(full waveform inversion)^[115]是近年来发展较快的地球物理前沿技术之一。它与最小二乘偏移的提法很相似。其中的正演算子

$ \mathit{\boldsymbol{F}}\left( v \right) = {D_{{\rm{syn}}}} $

(87)

由波动方程(1)来定义。与最小二乘偏移不同的是正演算子((87)式)相对于速度变化并不是一个线性算子。这里我们要寻找速度函数v, 使得下面的目标函数最小

$ \begin{gathered} E\left( v \right) = \frac{1}{2}{\left\| {\mathit{\boldsymbol{F}}\left( v \right) - D} \right\|^2} = \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\frac{1}{2}\iiint {{{\left( {{D_{{\rm{syn}}}} - D} \right)}^2}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}{\rm{d}}t} = \min \hfill \\ \end{gathered} $

(88)

即计算得到的数据D_syn与实际观测数据D最为接近。

给定背景速度v₀, 对波动方程(1)做速度扰动δv, 可以得到扰动波场δp满足如下方程:

$ \left( {\frac{1}{{v_0^2}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - \Delta } \right){\rm{ \mathsf{ δ} }}p\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = 2\frac{{{\rm{ \mathsf{ δ} }}v}}{{v_0^3}}\frac{{{\partial ^2}{p_0}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}} $

(89)

式中背景波场p₀(x_r; t; x)由公式(3)定义。

用常规方法求解优化问题(88)式需要计算梯度算子:

$ \frac{{{\rm{ \mathsf{ δ} }}E}}{{{\rm{ \mathsf{ δ} }}v}} = {\left( {\frac{{{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{F}}}}{{{\rm{ \mathsf{ δ} }}v}}} \right)^{\rm{T}}}\left( {{D_{{\rm{syn}}}} - D} \right) = {\left( {\frac{{{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{F}}}}{{{\rm{ \mathsf{ δ} }}v}}} \right)^{\rm{T}}}\left( {{\rm{ \mathsf{ δ} }}D} \right) $

(90)

由方程(89), 我们可以得到梯度算子的具体计算方法:首先求解两个方程

$ \left\{ \begin{gathered} \left( {\frac{1}{{v_0^2}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = f\left( t \right)\delta \left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) \hfill \\ \left( {\frac{1}{{v_0^2}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \hfill \\ \;\;\;\;\;\;\;\;\int {{\rm{ \mathsf{ δ} }}D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)\delta \left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{\rm{r}}}} \right){\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}} \hfill \\ \end{gathered} \right. $

(91)

式中的第一个方程是沿时间方向正向求解, 第二个方程是沿时间方向反向求解。梯度算子由上面两个波场的互相关计算得出:

$ \frac{{{\rm{ \mathsf{ δ} }}E}}{{{\rm{ \mathsf{ δ} }}v}}\left( \mathit{\boldsymbol{x}} \right) = \iint {\frac{2}{{v_0^3}}\frac{{{\partial ^2}{p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}}{p_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){\rm{d}}t{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}} $

(92)

实际上, 全波形反演梯度算子的计算类似于求解一个逆时偏移。由最速下降方法, 速度更新可以近似由下面的公式迭代求解:

$ \begin{gathered} {\rm{ \mathsf{ δ} }}v\left( \mathit{\boldsymbol{x}} \right) = - \alpha \left( \mathit{\boldsymbol{x}} \right)\frac{{{\rm{ \mathsf{ δ} }}E}}{{{\rm{ \mathsf{ δ} }}v}}\left( \mathit{\boldsymbol{x}} \right) = \hfill \\ - 2\alpha \left( \mathit{\boldsymbol{x}} \right) = \iint {\frac{1}{{v_0^3}}\frac{{{\partial ^2}{p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{\partial {t^2}}}{p_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){\rm{d}}t{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}} \hfill \\ \end{gathered} $

(93)

α(x)是通过搜索得到的最优系数。该速度更新在频率域的Kirchhoff型表达式为:

$ \begin{array}{*{20}{c}} {{\rm{ \mathsf{ δ} }}v\left( \mathit{\boldsymbol{x}} \right) = 4{\rm{ \mathsf{ π} }}\alpha \left( \mathit{\boldsymbol{x}} \right)\iiint {\frac{1}{{v_0^3}}{\omega ^2}{A_{\rm{s}}}{A_{\rm{r}}}{{\rm{e}}^{{\rm{i}}\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)}}} \cdot } \\ {{\rm{ \mathsf{ δ} }}D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}{\rm{d}}\omega } \end{array} $

(94)

另一方面由公式(89), 我们可以得到

$ \begin{array}{*{20}{c}} {{\rm{ \mathsf{ δ} }}\hat p\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = - \int {2{\omega ^2}\frac{{{\rm{ \mathsf{ δ} }}v}}{{v_0^3}} \cdot {{\hat p}_0}\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;\mathit{\boldsymbol{x}}} \right) \cdot } } \\ {{{\hat p}_0}\left( {\mathit{\boldsymbol{x}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right){\rm{d}}\mathit{\boldsymbol{x}}} \end{array} $

(95)

由公式(3)所定义的背景波场p₀(x_r; t; x)也可以看作是声波方程的Green函数。将其高频渐近表达式((6)式)代入(95)式得到:

$ {\rm{ \mathsf{ δ} }}\hat p\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = - \int {2{\omega ^2}\frac{{{\rm{ \mathsf{ δ} }}v}}{{v_0^3}}{A_{\rm{s}}}{A_{\rm{r}}}{{\rm{e}}^{{\rm{i}}\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)}}{\rm{d}}\mathit{\boldsymbol{x}}} $

(96)

实际上显式地给出了导算子δF/δv定义:

$ \left( {\frac{{{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{F}}}}{{{\rm{ \mathsf{ δ} }}v}},{\rm{ \mathsf{ δ} }}v} \right) = - \int {2{\omega ^2}\frac{{{\rm{ \mathsf{ δ} }}v}}{{v_0^3}}{A_{\rm{s}}}{A_{\rm{r}}}{{\rm{e}}^{{\rm{i}}\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)}}{\rm{d}}\mathit{\boldsymbol{x}}} $

(97)

所以, 优化问题(88)式的解近似地表示为:

$ {\rm{ \mathsf{ δ} }}v \approx {\left[ {{{\left( {\frac{{{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{F}}}}{{{\rm{ \mathsf{ δ} }}v}}} \right)}^{\rm{T}}}\frac{{{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{F}}}}{{{\rm{ \mathsf{ δ} }}v}}} \right]^{ - 1}}{\left( {\frac{{{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{F}}}}{{{\rm{ \mathsf{ δ} }}v}}} \right)^{\rm{T}}}\left( {{\rm{ \mathsf{ δ} }}\hat D} \right) $

(98)

仿照第一节反演反射系数的方法, 可以得到反演公式(98)在反射角度域的渐进表达式^[116-118]:

$ \begin{gathered} \frac{{{\rm{ \mathsf{ δ} }}v}}{{{v_0}}} \sim \frac{{32{\rm{ \mathsf{ π} }}{v_0}}}{{\sin \theta }}{\cos ^2}\theta \iint {\iiint {\frac{{\cos {\alpha _{\rm{s}}}\cos {\alpha _{\rm{r}}}}}{{{v_{\rm{s}}}{v_{\rm{r}}}}}{A_{\rm{s}}}{A_{\rm{r}}}{{\rm{e}}^{{\rm{i}}\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)}}{\rm{ \mathsf{ δ} }}\hat D}} \cdot \hfill \\ \;\;\;\;\;{\rm{ \mathsf{ δ} }}\left( {\varphi - \varphi '} \right){\rm{ \mathsf{ δ} }}\left( {\theta - \theta '} \right){\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}{\rm{d}}\omega {\rm{d}}\theta '{\rm{d}}\varphi ' \end{gathered} $

(99)

上式的左边与角度θ和φ无关。如果我们知道地下每个点在角度域的成像范围[θ_min, θ_max]×[φ_min, φ_max], 并计算它的面积L_{θ, φ}(x), 对(99)式两边积分, 就得到:

$ \begin{gathered} \frac{{{\rm{ \mathsf{ δ} }}v}}{{{v_0}}} \sim \frac{{32{\rm{ \mathsf{ π} }}}}{{{L_{\theta ,\varphi }}\left( \mathit{\boldsymbol{x}} \right)}}\iiint {\frac{{\cos {\alpha _{\rm{s}}}\cos {\alpha _{\rm{r}}}}}{{{v_{\rm{s}}}{v_{\rm{r}}}}}\frac{{{v_0}{{\cos }^2}\theta }}{{\sin \theta }}{A_{\rm{s}}}{A_{\rm{r}}}{{\rm{e}}^{ - {\rm{i}}\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)}}} \cdot \hfill \\ \;\;\;\;\;\;\;\;\;{\rm{ \mathsf{ δ} }}\hat D{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}{\rm{d}}\omega \hfill \\ \end{gathered} $

(100)

也可以类似推导共炮域的反演公式^[117]:

$ \frac{{{\rm{ \mathsf{ δ} }}v}}{{{v_0}}} \sim \frac{2}{{{L_{{\mathit{\boldsymbol{x}}_x}}}\left( x \right){\rm{ \mathsf{ π} }}}}\iint {\frac{{\cos {\alpha _{\rm{r}}}}}{{{v_{\rm{r}}}}}\frac{{{A_{\rm{r}}}}}{{{A_{\rm{s}}}}}{{\cos }^2}{{\rm{e}}^{ - {\rm{i}}\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)}}{\rm{ \mathsf{ δ} }}\hat D{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}{\rm{d}}\omega } $

(101)

其中, L_xx(x)是对成像点x有贡献的所有地面炮点的面积。

我们可以将基于梯度算子的反演公式(94)和基于逆算子的反演公式(100)或公式(101)作一个对比。这些公式都说明通过适当修改成像过程, 可以建立全叠加成像剖面与速度扰动之间的联系。公式(94)是基于梯度下降法, 只给出了速度扰动的变化方向。而公式(100)或公式(101)则是基于真正的反演, 近似求出了逆算子, 至少可以用来指导全波形反演的叠代优化策略, 或者用来作为反演求解过程中的预条件子以加速迭代收敛。

类似于第2节的方法, 我们可以在逆时偏移的框架下近似实现公式(100)或公式(101)。以公式(101)为例^[119], 我们需要分别计算正向波场

$ \left( {\frac{1}{{v_0^2}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_{\rm{s}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = f\left( t \right)\delta \left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) $

(102)

和反向波场

$ \left\{ \begin{gathered} \left( {\frac{1}{{v_0^2}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - {\nabla ^2}} \right){p_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = 0 \hfill \\ {p_{\rm{r}}}\left( {{x_{\rm{r}}},{y_{\rm{r}}},{z_{\rm{r}}} = 0;t} \right) = {\rm{ \mathsf{ δ} }}D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) \hfill \\ \end{gathered} \right. $

(103)

此时的成像条件需要修改为:

$ {\rm{ \mathsf{ δ} }}v\left( \mathit{\boldsymbol{x}} \right) = \frac{{v_0^3\left( \mathit{\boldsymbol{x}} \right)}}{{2{\rm{ \mathsf{ π} }}{L_{{\mathit{\boldsymbol{x}}_x}}}\left( \mathit{\boldsymbol{x}} \right)}}\Delta \left( {\iint {\frac{{{{\hat p}_{\rm{r}}}\left( {\mathit{\boldsymbol{x}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{{{{\left( {{\rm{i}}\omega } \right)}^3}{{\hat p}_\mathit{\boldsymbol{s}}}\left( {\mathit{\boldsymbol{x}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)}}{\rm{d}}\omega {\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}}} \right) $

(104)

式中引入了Laplace滤波器来表示反射角的余弦函数^[27]:

$ {\cos ^2}\theta = \frac{{{v^2}{{\left| \mathit{\boldsymbol{k}} \right|}^2}}}{{4{\omega ^2}}} $

(105)

本节讨论了全波形反演和偏移之间的关系。我们既可以认为全波形反演是迭代型的偏移过程, 也可认为有两类真振幅偏移^{[28, 118]}:一类是输出偏移道集, 用来反演反射系数的成像过程, 如公式(17)和公式(18), 其结果可以进行地下岩性的振幅分析(AVO, AVA); 另一类是输出全叠加剖面, 用来反演速度等物性参数的反演过程, 如公式(100)和公式(101)。一些计算实例参见文献[28]和文献[118]。

6 全波形反演的问题与探讨

全波形反演是在2005年前后, 几乎与逆时偏移同时进入工业界的视野, 至今已有十余年。继逆时偏移全面投入叠前深度偏移生产之后, 全波形反演一直是地球物理界的攻关重点之一。在许多探区, 特别是海上三维地震数据处理中取得了很多实际应用效果。该方法特别擅长精细刻画浅层速度模型^[120]。理想的全波形反演是一套自动化的地震资料处理过程, 能够直接从原始数据出发, 给出高精度的速度模型。但是在实际生产中, 现阶段全波形反演扮演一个速度建模的辅助角色, 还有待于理论和实用方法上的突破。本节探讨全波形反演在实际应用中面临的主要困难, 并且介绍文献中所提出的一些可能的解决方法。

6.1 低频缺失

经典的全波形反演算法建立在小扰动理论基础上, 即要求模拟数据与观测数据充分接近。简单地说, 就是要求ωΔτ要足够小。这里的Δτ是模拟数据与实际观测数据之间的走时差。当初始模型严重失真时, Δτ比较大, 需要依靠数据中的低频成分对模型做出校正。由于海上数据有效的低频信号通常在3~4Hz, 这就使得全波形反演必须依靠其它速度建模手段得到一个与真实速度比较接近的初始模型, 在此基础上才能反演出更细致的结构。这个小扰动假设极大地限制了全波形反演在生产中的应用。文献中相应出现了一些应对手段。

6.1.1 将数据改造为低频信号

这类方法很多。例如, 可以比较模拟数据与观测数据的规范能量模之差^[121-122]:

$ \begin{gathered} \mathop {\min }\limits_v \frac{1}{2}\iiint {\left[ {\frac{{\int_0^t {{{\left| {{D_{{\rm{syn}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)} \right|}^2}{\rm{d}}s} }}{{\int_0^{{T_{\max }}} {{{\left| {{D_{{\rm{syn}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)} \right|}^2}{\rm{d}}s} }} - } \right.} \hfill \\ \;\;\;\;\;\left. {\frac{{\int_0^t {{{\left| {D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)} \right|}^2}{\rm{d}}s} }}{{\int_0^{{T_{\max }}} {{{\left| {D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};\omega ;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)} \right|}^2}{\rm{d}}s} }}} \right){\rm{d}}t{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}} \hfill \\ \end{gathered} $

(106)

也可以比较两者的包络函数之差^[123-124]:

$ \begin{array}{*{20}{c}} {\mathop {\min }\limits_v \frac{1}{2}\iiint {\left\{ {\sqrt {{{\left[ {{D_{{\rm{syn}}}}\left( t \right)} \right]}^2} + {{\left\{ {H\left[ {{D_{{\rm{syn}}}}\left( t \right)} \right]} \right\}}^2}} - } \right.}} \\ {{{\left. {\sqrt {{{\left[ {D\left( t \right)} \right]}^2} + {{\left\{ {H\left[ {D\left( t \right)} \right]} \right\}}^2}} } \right\}}^2}{\rm{d}}t{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}} \end{array} $

(107)

式中:H[D(t)]表示函数D(t)的Hilbert变换。文献[125]对几种目标函数作了对比, 说明经过上面的修改后, 目标函数的局部极值陷阱减少, 可以部分地补偿信号中缺少低频的限制。

6.1.2 拓宽数据频带

信号处理人员都知道, 在引进了稀疏反射系数的假设下, 可以将有限带宽信号延拓到未能记录到的低频段和高频段中去。此类方法在全波形反演中也得到了应用。例如, 可以认为接收到的数据D(t)是有限频带的震源子波s(t)和单脉冲响应的反射信号G(t)的褶积:

$ \hat D\left( \omega \right) = \hat s\left( \omega \right)\hat G\left( \omega \right) $

(108)

为了延拓D(t)的频谱, 可以假设G(t)为稀疏的反射序列, 于是可以寻找函数G(t), 在公式(108)约束条件下, 使得下面的L₁模泛函最小:

$ \int {\left| {G\left( t \right)} \right|{\rm{d}}t} = \min $

(109)

上面的算法实际上是传统的稀疏反褶积。由于全波形反演是一个高维的计算过程, 因此, 对输入数据做单道处理是不够的, 频谱延拓的结果还需要满足检波域和炮域的连续性^[126]。LI等^[127]引进了更为复杂的相位追踪方法, 假设地震数据具有稀疏性和在空间、频率两个方向上的连续性

$ \hat D\left( {\omega ,x} \right) = \hat f\left( \omega \right)\sum\limits_{j = 1}^N {{a_j}\left( {\omega ,x} \right){{\rm{e}}^{{\rm{i}}bj\left( {\omega ,x} \right)}}} $

(110)

在强制振幅a_j(ω, x)和相位b_j(ω, x)光滑性的条件下可以实现地震信号的频谱延拓。

6.1.3 限制解空间

当波场缺少低频的时候, 目标函数(88)式的解很可能不唯一。由于反演过程中的速度模型修正量是从波场的互相关计算得来的, 如果每次迭代所得到的反演算子梯度场都是高频震荡, 就无法用来反演出低频的速度变化。ESSER等^[128]认为, 可以通过限制反演问题的解空间来缩小解的不唯一性, 并且补充数据中的低频缺失。例如, 可以在要求数据匹配:

$ {\left\| {\mathit{\boldsymbol{F}}\left( v \right) - D} \right\|^2} \leqslant \varepsilon $

(111)

的前提下, 强制所反演的速度函数的全变差在L₁模下极小:

$ {\left\| {v\left( x \right)} \right\|_{TV}} = \int {\left( {\left| {\frac{{\partial v}}{{\partial x}}} \right| + \left| {\frac{{\partial v}}{{\partial y}}} \right| + \left| {\frac{{\partial v}}{{\partial z}}} \right|} \right){\rm{d}}\mathit{\boldsymbol{x}}} = \min $

(112)

这样在反演的过程中, 由新的约束泛函公式(111)和公式(112)所确定的梯度场中会自动补充一些低频的速度变化成分, 以减少不必要的高频变化。采用这种方法甚至能够反演出复杂而巨大的盐丘速度异常。此类方法见效的关键是根据先验知识提出适当泛函空间来定义速度特征, 并且降低反演问题解的不确定性。

6.1.4 直接反演大尺度速度模型

一些研究人员尝试着直接反演较大尺度的速度变化, 使得全波形反演在缺少低频数据的情况下也能够摆脱对初始速度的依赖。最简单方法是在给定的网格和速度变化范围内, 由Monte Carlo法随机生成很多背景速度, 从中选出最优解。但是这需要大量的正演计算, 很难付诸于现实。于是, 人们提出在可变的粗网格上由浅到深逐步做这样的测试。甚至于更进一步, 随着反演过程的演化, 在每次偏移过程中根据地下构造来划分网格或构造区块^[129], 以减小大尺度速度模型的自由度。为了加速收敛, 可以引进模拟淬火法^[130]或遗传基因算法^[131-132]。此外, 还有作者提出了采用水平集合(level-set)算法来反演巨型盐丘体的建议^[133-134]。

6.2 振幅和走时

常规的速度建模方法, 如层析成像, 是根据走时误差来调整速度模型。而经典全波形反演中的速度建模是通过比较模拟数据与观测数据的振幅差异来考查速度的准确性。只有当模拟数据和观测数据的相位差在半个波长范围内, 才可以近似认为振幅的差别与走时的差别是成正比的。否则, 全波形反演就出现跳周现象(cycle skipping), 陷入局部极小的陷阱, 从而误导速度反演。地球物理学家一直都在寻求一种能够从数据中直接提取走时信息来反演速度的全波形反演方法。

LUO等^[135]于1991年提出了通过模拟数据和观测数据互相关来计算走时误差的方法, 并且建立了速度扰动和走时误差之间的反演关系。ZHANG等^[136]将波场从检波器反向传播到0时刻, 通过考察反传后的波场在炮点附近是否聚焦来衡量走时误差, 并更新速度模型。这些方法都试图避免直接用振幅信息来反演速度。

由于观测到的地震数据是反射、透射、折射等各种波的复合叠加, 并且受到假频和噪声干扰, 因此, 如何从中准确地提取走时信息是一个难题。随着图像处理技术的进步, 近期出现了很多新的技术。例如, 可以对模拟数据做局部时间和空间移动来比较其与观测数据的相似度, 以动态整形的方法^[137]计算每个数据点上的走时误差, 并应用于反演^[138-140]。也可以仿照运输问题, 用最优运输距离来衡量将模拟数据中叠加的地震波形分别移动至观测数据上的对应位置所需要付出的时空搬移成本^[141-142]。这比直接在最小二乘意义下求数据残差的方法更能够反映速度误差对于数据中走时信息的影响。

WARNER等^[143]于2016年提出了一种自适应波形反演方法:在迭代的每一步构造一个匹配滤波器w(t), 使得

$ \left( {{D_{{\rm{syn}}}} * w} \right)\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) \approx D\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) $

(113)

同时引进一个权函数T(t)。例如, 可以假定T(t)随时间递增。在此基础上, 把反演问题提为寻找速度模型使得下面的目标函数最小:

$ \mathop {\min }\limits_v \frac{1}{2}\frac{{{{\left\| {T\left( t \right)w\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)} \right\|}^2}}}{{{{\left\| {w\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)} \right\|}^2}}} $

(114)

这实际上是要求在理想的速度模型下匹配滤波器尽量接近一个单脉冲, 从而使得模拟数据和观测数据最大限度地相似。文章中对各类速度反演方法作了探讨与评述, 值得一读。

6.3 反射、折射与透射

地面所观测到的地震数据主要包含了由地下构造变化而产生的反射波、在临界角附近出现的折射波(或称首波, head wave)和由于速度梯度增加从地下返回地面的透射波(或称潜波, diving wave)。人们普遍认为反射波的振幅至少由速度和密度变化来确定, 比较复杂。而折射波和透射波的振幅则主要取决于速度的变化。特别地, 折射波和透射波对速度建模的长波长分量有更大的贡献^[144]。所以在全波形反演进入生产应用初期, 人们倾向于仅仅利用折射波和透射波来反演速度。这需要在反演的预处理中采用空变的时窗切除方法将输入数据中的反射波切除掉。但是, 在地震数据里各种波信号经常重叠在一起, 生硬地切除一部分会导致信号扭曲, 影响反演效果。另一方面, 这3种波现象并非泾渭分明, 当速度变化梯度较小时它们之间的过渡变得很模糊。总之, 只把全波形反演应用于某一部分观测数据, 理论上是合理的, 但在实践上难以操作, 而且也浪费了观测数据中大量和反射有关的地质信息。

在全波形反演过程中, 速度的波长(|k|)和反射角(θ)、波场频率(ω)、声波速度(v)之间有如下关系^[145]:

$ \left| \mathit{\boldsymbol{k}} \right| = \frac{{2\omega \cos \theta }}{v} $

(115)

也就是说当观测数据缺失低频的时候, 长波长的速度修正量可以从大角度的反射信号得到。在KHALIL等^[146]研究成果的基础上, ALKHALIFAH^[147]采用延时道集的方法, 将反演过程的梯度场分解为不同的反射角道集, 用来控制由反射信号所反演得到的速度尺度。

XU等^[148]于2012年提出了反射型全波形反演方法。这种方法的实质是在全波形反演的过程中加入反偏移, 用模拟生成的反射数据来匹配观测的反射数据, 从而指导速度反演。在此基础上又发展出了一系列新的方法, 可以将反射波和透射波都纳入到全波形反演当中^[149-150]。

SUN等^[151]于2016年提出了应用于实际数据处理的基于反射波的全波形反演方法。该方法在每次反演迭代过程中, 给定背景速度v(x), 可以由逆时偏移成像算法得到反射系数R(x)。再通过反偏移算法得到所预测的地面反射数据, 这个过程记为:

$ \mathit{\boldsymbol{F}}\left[ {v,R\left( v \right)} \right] = {D_{{\rm{syn}}}}\left( {{\mathit{\boldsymbol{x}}_{\rm{r}}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) $

(116)

下一步对于每一条地震道找到一个类似于(113)式所定义的匹配滤波器w(x_r; t; x_s)。理想情况下, 我们要求滤波器w(x_r; t; x_s)对于每一条地震道都是相似的, 并且接近于一个点脉冲函数。即求解速度函数v(x), 使得下面的泛函极小

$ \mathop {\min }\limits_v \frac{1}{2}{\left\| {A\left( w \right)} \right\|^2} $

(117)

式中的函数A既可以用来衡量匹配滤波器随接收点位置的变化

$ A\left( w \right) = \frac{1}{{\left\| w \right\|}}\left( {\left| {\frac{{\partial w}}{{\partial {\mathit{\boldsymbol{x}}_{\rm{r}}}}}} \right| + \left| {\frac{{\partial w}}{{\partial {y_{\rm{r}}}}}} \right|} \right) $

(118)

也用来作为罚函数促使匹配滤波器接近单点脉冲

$ A\left( w \right) = \frac{{\left| t \right|w}}{{\left\| w \right\|}} $

(119)

很显然, 该方法和前面所介绍的自适应波形反演方法有相通之处, 也受到了William Symes教授所提出的微分同相优化方法^[152-154](differential semblance optimization, DSO)的影响。

前述反演方法大多基于数据匹配, 即利用在地面观测到的波动现象来反推地下的速度。受限于观测系统, 很多的波场现象, 特别是绝大部分的透射波未能被记录到。因而完全靠地面的接收数据还不足以可靠地重建地下波场。VAN LEEUWEN等^[155]因此提出了波场重建反演方法(wavefield reconstruction inversion)。首先定义正演算子:

$ \mathit{\boldsymbol{F}}\left( v \right)p = s $

(120)

式中:F(v)是波动算子; p是压力波场; s是震源子波。同时定义了在检波点的观测算子。

$ \mathit{\boldsymbol{ \boldsymbol{\varGamma} }}p = D $

(121)

波场重建反演的目的是在反演过程中同时寻找速度和压力波场, 使它们尽量满足上面的两个方程, 即目标函数为^[156]:

$ \mathop {\min }\limits_{v,p} \frac{1}{2}{\left\| {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}p - D} \right\|^2} + \frac{\lambda }{2}{\left\| {\mathit{\boldsymbol{F}}\left( v \right)p - f} \right\|^2} $

(122)

该方法的实质是从震源和观测记录两个方向来控制反演, 尽量减少全波形反演结果中的非物理性质, 从而得到一些有意义的结果。原方法是在频率域提出的, WANG等^[157]将其改进并推广到了时间域。

6.4 多参数反演

前面讲到反射波的振幅不仅受速度变化的影响, 而且也与密度的变化有关。这就有必要在正演模拟中考虑加入密度函数(ρ)

$ \left( {\frac{1}{{{v^2}}}\frac{{{\partial ^2}}}{{\partial {t^2}}} - \rho \nabla \frac{1}{\rho } \cdot \nabla } \right)p\left( {\mathit{\boldsymbol{x}};t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) = \delta \left( t \right)\delta \left( {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}_{\rm{s}}}} \right) $

(123)

类似于第5节的讨论, 可以得到反演公式^[118]

$ \begin{gathered} \frac{{{\rm{ \mathsf{ δ} }}v}}{{{v_0}}} + {\cos ^2}\theta \frac{{{\rm{ \mathsf{ δ} }}\rho }}{{{\rho _0}}} \sim \frac{{32{\rm{ \mathsf{ π} }}}}{{{L_{\theta ,\varphi }}\left( x \right)}}\iiint {\frac{{\cos {\alpha _{\rm{s}}}\cos {\alpha _{\rm{r}}}}}{{{v_{\rm{s}}}{v_{\rm{r}}}}}\frac{{{v_0}{{\cos }^2}\theta }}{{\sin \theta }}} \cdot \hfill \\ \;\;\;\;\;\;\;{A_{\rm{s}}}{A_{\rm{r}}}{{\rm{e}}^{ - {\rm{i}}\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)}}{\rm{ \mathsf{ δ} }}\hat D{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}{\rm{d}}\omega \hfill \\ \end{gathered} $

(124)

这实际上说明如果在全波形反演过程中产生反射角道集, 就可以从振幅随角度的变化关系来反演速度和密度函数。QIN等^[119]推导了共炮域的速度和密度反演公式, 并且根据反射波重建速度和波阻抗函数。

由于全波形反演依靠透射波和较大角度的反射波来修正长波长的速度分量, 而且这两类波场主要沿水平方向传播, 因此, 与层析成像主要依靠近角度的同相反射相比, 全波形反演受到更多各向异性因素的影响。注意反演公式(124)与SHUEY^[158]所提出的声波形式的AVO公式是相容的。根据这个线索, 我们可以按照各向异性下的AVO公式将公式(124)推广到VTI情形^[159]

$ \begin{gathered} \frac{{{\rm{ \mathsf{ δ} }}v}}{{{v_0}}} + {\cos ^2}\theta \frac{{{\rm{ \mathsf{ δ} }}\rho }}{{{\rho _0}}} + \Delta {\rm{ \mathsf{ δ} }}{\sin ^2}\theta {\cos ^2}\theta + \Delta \varepsilon {\sin ^4}\theta \sim \hfill \\ \;\;\;\;\;\;\frac{{32{\rm{ \mathsf{ π} }}}}{{{L_{\theta ,\varphi }}\left( \mathit{\boldsymbol{x}} \right)}}\iiint {\frac{{\cos {\alpha _{\rm{s}}}\cos {\alpha _{\rm{r}}}}}{{{v_{\rm{s}}}{v_{\rm{r}}}}}\frac{{{v_0}{{\cos }^2}\theta }}{{\sin \theta }}}{A_{\rm{s}}}{A_{\rm{r}}}{{\rm{e}}^{ - {\rm{i}}\omega \left( {{\tau _{\rm{s}}} + {\tau _{\rm{r}}}} \right)}} \cdot \hfill \\ \;\;\;\;\;\;{\rm{ \mathsf{ δ} }}\hat D{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{r}}}{\rm{d}}{\mathit{\boldsymbol{x}}_{\rm{s}}}{\rm{d}}\omega \hfill \\ \end{gathered} $

(125)

有大量文献讨论了如何利用全波形反演来求取密度函数^[160-161]、波阻抗^{[119, 162]}、各向异性参数^[163-166]以及吸收衰减因子^[167], 甚至推导了如何反演TTI介质中各向异性的对称轴方向^[168]。

本文只是在扩展的声波方程框架下讨论成像与反演方法。因为受限于技术条件和计算能力, 这些方法还是当前工业化生产和应用的主流。但是, 粘弹性介质模型是对地下构造更好的近似。多分量地震勘探和陆地的数据处理促使我们在粘弹性介质模型的框架下去思考叠前深度偏移的技术发展。全波形反演至少从理论上为我们提供了反演各种粘弹性介质参数的方法, 可以充分利用观测数据得到高精度的储层构造与岩性描述, 其应用前景非常广阔。

7 总结与展望

逆时偏移在生产中的应用是叠前深度偏移发展中的一个里程碑。它避免了Kirchhoff、射线束、单程波等成像技术的很多理论假设和限制, 使得偏移算法的研发和应用更加简明、丰富。原则上讲, 逆时偏移能够直接对多震源进行成像, 同时也可以在偏移过程中比较自然地处理由地表边界条件所产生的鬼波和多次波信号。经过多年的努力, 各向异性和粘性吸收等复杂的波传播现象都可以在逆时偏移的生产应用软件中得到有效的处理。角度域共成像点道集的概念将起伏地表、海上拖缆、OBN、VSP, 甚至多次波成像等各类偏移纳入到统一的应用框架当中。

当前, 逆时偏移在深海区复杂构造勘探中起着决定性的作用, 但是在浅海和陆上资料处理中的应用还远不如Kirchhoff偏移普及。这主要有3个方面的原因:①较浅目的层需要高精度的构造成像结果。当频率增加1倍时, 理论上Kirchhoff偏移的计算量增加到8倍, 而逆时偏移的成本需要增加到16倍, 高精度的逆时偏移快速算法仍然是有待深入研究的课题。②在共成像点道集输出方面, 逆时偏移远不及Kirchhoff偏移灵活。反射角道集的输出在方法上和生产成本上还存在困难。③在具有复杂构造的浅海和陆地区域, 资料处理和速度建模能力还有待于提高。浅海的强烈多次波和陆地上的复合噪声干扰严重降低了资料品质, 也为层析成像建模带来困难。起伏地表情况下的高精度近地表模型构建一直是个难题, 浅海的复杂气云构造也只是在引进全波形反演后才开始得到解决。如果偏移速度既简单又平滑, 逆时偏移的优势也就无从发挥。

现在看来, 逆时偏移只是一个重要的过渡技术, 将来会融入到最小二乘偏移和全波形反演当中去。但是最小二乘偏移至今尚未规模化应用于生产。可以认为, 最小二乘偏移是常规偏移的加强版, 同时又具有和全波形反演类似的算法困难, 例如, 走时和振幅的匹配问题。在速度模型精度不够时, 成像质量和速度误差经常搅在一起。如果在全波形反演中能够实现同时反演速度和波阻抗, 那么也就相当于将最小二乘偏移并入到了全波形反演当中。

在过去的半个世纪, 受数据采集、理论方法和计算能力的限制, 地球物理工作者不得不将一个复杂的叠前深度偏移过程拆分为很多小的专门课题和小的流程步骤并分而治之, 寻求渐进式的发展。随着计算机技术的发展和采集能力的不断提高, 或许我们应该思考能否将以前分解的步骤逐步组合起来, 最大限度地发挥计算效力, 实现人工和时间的节约。全波形反演当然是沿着这个思路的重要代表技术。之所以提出这样的想法是基于以下的观察与思考:

1) 传统处理流程中的震源反褶积、多炮资料分离、去鬼波、多次波压制等信号处理流程在原理上都是高维算子, Q补偿则本质上是一个深度域算子。受资料质量和处理成本的限制, 在时间域正确地实现这些流程是有困难的。特别是地表多次波和层间多次波的预测, 在时间域的计算量就超过了一次甚至多次全工区的逆时偏移。而这些问题都是可以由波动方程加上适当的边界和震源条件来直接预测, 在全波形反演框架下可以比较自然地解决。

2) 现阶段资料采集中的假频、不规则性、稀疏性和数据缺失问题主要由高维数据规则化来解决。这种时间域的解决方案往往忽略了波传播的特性, 而且处理质量随着资料品质、算法能力和实现手段不同而参差不齐。这些问题实际是最小二乘偏移所应该解决的问题, 也可以纳入到全波形反演中来。

3) 尽管全波形反演的建模质量受到低频缺失、跳周等问题的困扰, 生产中已经出现了一些有效的克服方法。实际上, 传统的速度建模手段, 例如层析成像, 通过同相轴拾取来克服低频缺失、跳周等问题。这就要求我们不要拘泥于传统的共轭算子的反演求解算法, 而是应该借鉴其它技术中行之有效的手段, 使全波形反演方法兼收并蓄, 最终臻于大成。

4) 现阶段全波形反演在浅海和深海勘探的广泛应用已经毋庸置疑, 能够解决一些浅层、高精度、复杂构造的速度建模问题。但是总的来说, 还是在传统建模方法的基础上作一些局部化和精细化的补充。在2016年的SEG年会上, Exxon公司报告了全波形反演在深海窄方位角数据的应用结果。信号中缺失0~5Hz频率成分信息, 最高反演频率达到40Hz, 得到了高精度的速度模型和波阻抗模型。Schlumberger公司提出了新的海上生产处理流程:第一步直接用全波形反演在原始数据上构建浅层速度模型; 同时对数据做预处理, 然后采用层析成像与反射波全波形反演相结合构造出高精度的速度模型。这些工作和结果扩展了全波形反演在生产中的作用, 向前迈出了重要的一步。

5) 陆上资料的全波形反演已经有一些应用实例, 虽然还不具有广泛代表性, 但是想要完全解决近地表建模问题, 除去全波形反演, 目前似乎还找不到其它理论上合理而又适用范围广泛的方案。

6) 过去的10年中, 与信息化计算和处理有关的技术快速发展, 大数据挖掘和深度学习取得了巨大的进展。地震资料处理本质上也是从海量数据中寻求地质信息。压缩感知、模式识别、图像处理、特征提取、高维时频分析甚至人工智能学习等技术都和反演问题有着密切的联系。即使是最为复杂的数据噪声, 也可以在反演中通过限制解空间、选择适当变换, 并且在迭代过程中根据波场的可预测性来逐渐剔除。这需要对资料中的信息有更深入的理解, 更统一有效的处理工具。反演技术的发展将使得勘探地球物理在更广阔的领域与其它学科相结合, 成为广义信息技术革命的一部分。

还能列举出很多的理由。一旦目标确定, 我们就应该全力以赴。

致谢: 文章中的很多工作是在前Veritas和CGG两个公司完成的, 在此向所有合作者和同事表示感谢, 特别要感谢徐升、Gilles Lambare、Robert Soubaras、Samuel Gray、James Sun、Carl Notfors、周洪波、汪道柳、Leon Chernis、David Yingst、谢毅、段炼、秦伯涛、唐兵、Andrew Ratcliffe、Vetle Vinje、Gram Robert、Gordon Poole、Adel Khalil、陈丰、张厚柱、张钋、Barry Hung、Matthew Karazincir等人。今年是我的博士生导师张关泉教授逝世5周年, 作为学生, 我的工作深受张老师的影响, 总希望不辜负导师的教诲与期望。同时也祝愿我的另一位学术前辈和合作者, Norman Bleistein教授身体健康。