GeophyAI

Python与地球物理数据处理

0%

反问题

11. 反问题例子:线性化地震数据反演

我们必须对间接测量得到的地震数据使进行处理和解释之后,才能较为合理地推断地下模型。我们需要对地下介质的物理性质进行建模,例如确定岩层之间的界面位置、估计岩层物理性质等。由一个给定数据集到与该数据集相关联的最优地下模型的这一过程被称为求解反问题。在这一节中,我们介绍石油地震学家解决这一问题所采取的一些方法。

分步式反演

反问题的基础要素

反问题存在于诸多领域,例如量子力学、医学成像等。但是,较为广义的反问题理论是由研究地球物理数据的人们发展起来的。地球物理学家试图了解地球内部,但唯一可用的数据是在地球表面、海面、海底或钻孔等位置收集的。
脱离于学科概念,反问题的解决通常要包含以下三个步骤:

  1. 建立模型参数。即通过可用数据构建最佳的模型参数。但是由于我们数据的限制,这一般是不可能的。例如,根据以下$AVA$响应$f_{ava}$确定两个各向同性弹性岩层之间的弹性参数的对比度(constrast)。当$\theta \in [0,30^{\circ}]$时,由于$sin^4\theta$较小,因此以上公式中只有两个参数(即$[\Delta I_P, \Delta I_S-3\Delta\mu]$),而非三个参数(即$[\Delta I_P, \Delta I_S-3\Delta\mu, \Delta \mu]$)。所以,由于数据的局限性,我们无法重构出精确的模型。物理量$\Delta I_P, \Delta I_S, \Delta \mu$分别代表纵波阻抗对比度、横波阻抗对比度和剪切模量对比度。
    11.1
  2. 正演(forward modeling)。对于给定的模型参数,该步骤意味着我们要在一定物理定律的约束下来预测该模型参数下的数据(如图11.1所示)。方程11.1给出了在特定弹性参数对比度下的$AVA$响应。波动方程也是正演的一种基本方式(如第2和6节所示),基于波动方程的正演模拟允许我们使用密度、拉梅系数等弹性参数来预测地震记录。
  3. 反演(inverse problem)。该步骤需要我们使用观测数据来推断模型参数,使得由该模型参数正演得到的预测数据(Predicted data)在特定的准则(Criterion)下能够很好地拟合观测数据(Observed data)。使预测数据和观测数据的误差平方和最小是最为常用的准则(如公式11.2所示)其中,$P_{obs}(x_s, x_r, t)$和$P_{pred}(x_s, x_r, t)$分别代表观测数据和合成数据。$x_s$表示炮点位置($s=1,2,…$);$x_r$表示接受点位置($r=1,2,…$);$t$表示旅行时$t_0\leq t\leq t_1$。

以上步骤中有非常强的反馈机制,在任何一步中的提升一般都可以用到其它两个步骤中。

在本节中我们只关注地震数据的反演。

非唯一性、不稳定性、收敛性、不确定性、成本

石油地震学中反问题的求解通常会遇到五个基础性问题:①非唯一性。即如何确定由给定数据所反演的模型与实际情况最为匹配。②不稳定性。即数据中的微小扰动可能会导致反演结果中的较大扰动。③收敛性。即迭代求解反问题时可能出现。④不确定性。即物理模型不准确、物理模型不完整性以及测量的不确定性等⑤成本。反问题求解过程中的正问题的成本。

11.2
为了使我们刚才提出的反问题解决方案的相关问题更加具体,让我们看几个例子。如图11.2所示,两个完全不同模型,入射角为$30^{\circ}$以下时的亚临界$AVA$响应(振幅随角度变化)几乎完全一样,此即非唯一性问题。对于超过$30^{\circ}$的入射角、$AVA$响应的差异足以区分这两种模型。换言之,石油地震学问题中大量的非均匀性可以通过改进我们的理论或抛弃一些假设来解决,或者通过改进我们采集数据的观测系统,例如,长偏移距,以便采集大于$30^{\circ}$入射角的数据。
11.3
图11.3显示了亚临界PP反射的$AVA$响应来重构横波速度所产生的不稳定性问题。我们可以看到,当对横波速度施加$50%$变化量后,$PP$反射的$AVA$响应几乎没有变化。这种不稳定性的结果就是,我们可能无法通过$PP$反射数据来区分松散沉积物和固结岩层。但是,正如在非唯一性问题中所讲的那样,我们可以改进理论、抛弃一些假设或者改进采集数据的观测系统来解决不稳定性问题。例如,在长偏移距数据能够看到明显的$PS$反射能量。事实上,从图11.3中可以看出,可以通过使用$PP$反射的临界后$AVA$响应或者$PS$的$AVA$响应来对两种模型加以区分。
11.4
接下来让我们看一下物理模型的不准确性问题。假设我们正演中用于预测$PP$反射的物理模型是基于声波方程,而不是弹性方程(见第2章和第6章)。如图11.4所示,尽管声波方程能够很好地预测旅行时,但是振幅确实完全不正确的,在大偏移距处更为明显。同样,这些不准确的结果是,我们可能无法区分某些岩性。

正如我们在第7~10章中所描述的,地震数据含有大量的多余的能量,例如地滚波、涌浪噪声、多次波等等。当反问题迭代求解时,这些不需要的信号(通常比期望的信号大得多)会减慢收敛到最优解的速度,甚至会导致地下模型的不精确。事实上,选择最佳地下模型的大多数标准包括拟合观测数据和为给定地下模型预测的数据。但是,根据这些标准建立的模型往往同时解释了数据中的可能是噪声的大振幅信号,例如地滚波等。解决这个收敛性问题的一种方法是在开始反演之前从数据中尽可能多地去除不必要的能量。

下面考虑石油地震学中正问题的成本。正问题即在一组初始条件、最终条件和边界条件下,求解控制波在地球中传播的微分方程。有限差分模拟是求解波动方程的常用方法,其对偏微分方程中的偏导数作了数值近似(附录C)。对时间和空间进行适当地离散后,就能够较为精确的计算波动方程中的偏微分。在这种情况下,有限差分正演模拟(FDM)技术是迄今为止通过复杂地质模型(如勘探与生产行业目前面临的模型)模拟弹性波传播的最精确的工具。但是,其计算成本非常高。建设我们对一个离散成$1000\times1000\times500$个空间网格的复杂地质模型进行弹性波正演模拟,在时间上递推4000次。我们估计,在SGI Origin 2000(由Silicon Graphics,Inc.提供)上,使用20个CPU进行40000次的小型3D正演计算需要12年以上时间,远远超过了某些石油储层的生产寿命(编者注:本书写于2005年)。

地震反问题的分步式方法

在实际应用中,我们并不是开发单一的反演算法,让其接收野外数据输入然后生成地下模型,而是将其划分为多个步骤,以更好地处理不稳定性、唯一性、不确定性、收敛性,正演成本等问题。这里的基本思想是,我们对数据、问题的物理性质以及其他先验地质和地球物理信息的了解越深入,越有助于引导反问题解远离不稳定性和非均匀性。这些知识甚至可以帮助我们减少正演成本,因为它允许我们采用具有多种人类交互可能性的多步反演方案,而不是单步方案,如图11.5所示。
11.5
当前海洋数据多步反演方案包含以下关键步骤:

  1. 去鬼波(即对震源和接收器鬼波的影响进行校正);
  2. 去多次波(即自由面多次波的衰减);
  3. 宏观模型速度的估计(即估计最适合地震数据中包含的走时的速度模型),这也被称为背景模型;
  4. 线性化反演(即,地震数据反演,假设数据已针对多次反射进行了校正,并且背景模型已知)。

陆地地震中,关键步骤如下:

  1. 地滚波衰减;
  2. 静校正;
  3. 宏观模型速度估计;
  4. 线性化反演。

井中地震关键步骤如下:

  1. 上下分离;
  2. 多次波衰减;
  3. 线性化反演。

注意,在井中地震中,速度模型通常可从测井测量中获得。因此,通常不需要进行速度估计。我们不能用三章内容来描述所有这些处理步骤。我们将重点讨论海洋数据多步反演的四个步骤。对这四个步骤的描述为有助于了解本章内容。

由于成本原因,其中的一些步骤,例如多次波去除,都是从几何光学或者纯信号处理领域发展来的(例如Chapter8中$fk$倾角滤波方法),并不是基于波动方程的。这些工具倾向于只关注旅行时变化而忽略振幅变化。并且,这些方法还会产生假象,并且在数据中留有较大的残差。但是由于这些方法对于计算资源的需求较小,因此应用较为广泛。但随着近年来计算能力的提升,人们不再考虑使用这些几何光学工具来替换基于波动方程的方法。此外,为达到现在地震资料处理中追求的高分辨率和多次波衰减的目的,我们也必须基于运动学和动力学更加精确的波动方程方法。这些方法背后的理论是地震学中用于表示波场的两个积分方程:$Lippmann-Schwinger$积分方程和表示定理。$Lippmann-Schwinger$积分方程导出了所谓的$Born$散射级数,而表示定理导出了$Kirchhoff$散射级数。上述反演的四个步骤(去鬼波、多次波衰减、宏观模型估计和线性化反演)均可以通过$Born$散射级数和$Kirchhoff$散射级数推导得出,正如我们在第10章多次衰减中看到的那样。本章的重点是线性化反演和宏观模型的估计

在前面的章节中,我们提到了NMODMO,偏移等作为地震成像的工具。我们刚才描述的地震反问题中,这些工具从何而来?在本章中,我们将看到偏移只是线性化反演的一个特例。我们还将动校正叠加与线性化反演、以及倾角校正叠加联系起来。

反问题示例的关键假设

在对数据进行去多次波和去鬼波处理(如Cahpters9和10所述)之后,下一步是从包含初至波的数据中构建地下模型。这一章节中我们通过线性化反演解来讲述如何解决反问题,线性化反演非常好地提供了一种从原始数据域构建地震剖面的方法。线性化反演问题的关键假设如下面的章节所述。

Born近似

如前所述,通过正演模拟求解物理规律允许我们以最经济有效的方式获得预测数据。对于给定的地下模型,我们本章节的任务是选择一组物理定律来预测只包含初至波的拖揽数据。
第6节中所述的一阶波动方程以及初始和边界条件,或者是$Lippmann-Schwinger$积分方程是两个解决正问题的工具。

基于有限差分求解正问题

使用波动方程求解只有初至波的正问题与求解包含多次波的完整数据几乎是完全相同的。我们将不施加如Chapter3所述的自由界面边界条件,而是假设水层是无限的。从实用角度看,当用有限差分计算波动方程的解时(见附录C),我们需要用吸收边界条件来代替自由界面边界条件。因此,就计算所消耗的时间而言,使数据仅包含初至波是没有任何收益的。实际上,在某些情况下,仅计算包含初至波的波场要比包含多次波的全波场数据更加费时,因为吸收边界条件的数值模拟要比自由界面边界条件计算成本更高。

基于Born近似求解正问题

正演模拟的另一种解决方案是$Lippmann-Schwinger$积分方程,或与该积分方程相关的$Born$散射级数(见Box 11.1)。这要求我们将实际的介质分解为背景模型(或参考模型、宏观模型)和扰动模型(散射模型)两部分。在这一分解中,震源和接收点均位于背景模型中,如图11.6所示。没有接触到扰动模型,直接从震源到达接收点的波被称为直达波场(非扰动波场),与扰动模型相互作用的波则构成了扰动波场(散射波场)。
11.6
对于给定的地下介质模型,如Box 11.1所示的$Lippmann-Schwinger$积分方程允许我们预测产生只包含初至的数据,这里用$P$来表示:

其中,

这里,$P_0$表示直达波场,$S$表示震源,$G_0$代表背景介质中波动方程的格林函数,$G$代表真实介质中波动方程的格林函数,$W$表示散射势。散射波场$\delta P$可表示为:

公式11.3对于任何背景介质都是有效的,而$G_0$是否准去的描述了震源检波点界面以及散射点并不重要。11.3所示的积分可以通过迭代的方式或者通过如Box 11.1所示的$Born$级数方式求解。
公式11.3还可以改写为:

如果我们假设$G_0$足够近似$G$,即$G_0$能够精确预测从散射点到震源和检波器的波场,有如下公式

式中$I$表示恒等算符——当$W^{‘}$近似等于$W^{‘}$,公式11.6退化为$Born$近似,即

或者

此外,在高于二阶项均可忽略不计的前提下,我们可以通过截断$Born$级数Box 11.1中,只取其前两项,来得到公式11.9和11.10。

因此,我们能够利用公式11.9预测给定散射势的初至波场(即求解正演模拟问题)。然而,与第六章中的$Lippmann-Schwinger$积分方程或一阶波动微分方程组相反,我们在将实际模型分解为背景和扰动模型时,不能简单地选择任意背景模型。只有当背景模型与源-接收器界面和散射点之间的实际介质相同时,方程11.9才有效。另请注意,与11.3中的$Lippmann-Schwinger$积分方程以及有限差分求解波动方程相反,利用$Born$近似得到的方程11.9求解正问题不涉及迭代过程或具有无穷多个项的级数。因此,方程11.9中的$Born$近似不会产生数值不稳定性问题,也不会遇到迭代过程中的收敛问题,其解也不会有无穷多项级数。这些特点使$Born$近似成为解决正演模拟问题的一个非常经济有效的方案,特别是对于背景介质中与波动方程相关的格林函数是明确的或已知解析解时,例如在均匀和一维背景介质中。

此外,其中包含我们最终想要恢复的弹性参数的数据和散射势,在$Born$近似下是线性的。这一特性大大简化了反问题的求解。$Born$近似的缺点是描述波在源-接收器界面和散射点之间传播的背景模型必须是精确的。否则,$Born$近似产生的响应误差较大。

让我们用另一种观点解释$Born$近似,正如第10章所述,初至波要求震源和检波器均位于地下,由于我们在本章中考虑的拖缆数据是在地球表面附近生成和收集的,因此我们必须以某种方式将数据从地球表面附近的实际震源和接收器位置外推到地下的所需位置。这种外推是通过计算描述波在背景介质中传播的格林函数$G_0$来实现的。为了计算这个格林函数,我们需要知道源所在位置和散射点之间介质,以及散射点和接收器位置之间介质的性质。换句话说,我们至少需要知道地下模型的某些部分,以便重建“整个”模型。这一要求被称为基于玻恩近似的地震成像悖论。由于目前大多数地震成像方案都是基于$Born$近似的(即使这种假设在一些出版物中往往不明确),这一要求通常被称为地震成像悖论。

总之,用本章所提到的$Born$近似,我们必须将“从包含初至的地震数据中重建地下模型”这一问题分解为两部分:①找到背景速度模型,以使得我们能够将数据从震源-检波器界面外推到散射点上;②重建散射模型,其中包含引起散射和反射的弹性性质。在接下来的四节中,我们将重点讨论“在背景模型已知前提下重建引起散射的物理属性”这一问题。关于如何估计背景模型的讨论将在标题为“偏移”和“估计背景速度的模型”的章节中进行。

平滑背景介质

如图11.7所示,方程11.9中的$Born$近似使我们能够预测初至以及多次波。由于在大多数成像中都假定数据只包含初至,因此我们需要调整11.9的$Born$近似,使预测数据满足成像假设。如果描述背景介质中波动方程的格林函数$G_0$不包含反射波,此时11.9的$Born$近似将只预测初至波。换句话说,如果我们假设背景模型是平滑的(即没有反射或只有非常弱的反射),那么11.9中的$Born$近似将只预测初至波。实现这种平滑的经典方法是将实际模型分解为两部分:低空间频率分量和高空间频率分量。$a1D$(横向不变)模型的分解示例如图11.8所示,低空间频率(长波数)分量表示背景模型,高空间频率(即短波数)分量引起散射和反射。
11.7~8
当背景模型中没有反射点时,如图11.8所示,描述波在背景介质中传播的直达波场基本上可以以旅行时和几何扩散来描述,而这两者主要取决于速度变化。这也是为什么在石油地震学中背景模型一般被称为背景速度模型

Born近似的局限性

Box 11.1所描述的那样,11.3所示$Lippmann-Schwinger $积分方程的解可以表示成$Born$级数的形式,其中此级数的前两项为$Born$近似,有时我们将这种近似称为单散射点近似。通过截断$Born$级数,只取其前两项,我们假设背景介质与实际介质十分相近,并且高于二阶的$Born$级数项均可以忽略不计。
下面我们以实际的声速介质为例,将$Born$近似的假设公式化。我们从$Box 11.1$中定义的散射势(散射势的表达式源自量子力学,$W$实际上是粒子散射的势,这一表达现在被那些将散射理论应用于地震学的学者们所广泛采用)出发:

其中$K_0({\rm \textbf x})$和$\sigma_0({\rm \textbf x})$分别为背景介质的压缩系数(体积模量的倒数)和比容积(密度的倒数)。$K({\rm \textbf x})$和$\sigma({\rm \textbf x})$分别为实际介质的压缩系数(体积模量的倒数)和比容积(密度的倒数)。在$Born$近似的假设中,背景介质中十分接近于实际介质,即如下公式成立

下面我们通过比较$Born$近似法的预测数据与精确解对应的数据,来更具体地观察公式11.13所提出的假设。
我们必须通过$Born$近似或者利用有限差分求解波动方程的数值解,是因为精确的解析解非常少。为数不多的精确解析解一般仅限于平面波的传播。因此,这里我们将通过一维平面波的传播比较$Born$近似和精确解析解。一维模型如图11.9所示,其由一个嵌入在速度为$V_{P0}=2.0 km/s$的背景介质中的流体介质组成,其中嵌入物速度为$V_{P1}$、厚度为$h$。因此,在这一模型中的散射势可以表示为:

如果我们将公式11.14带入公式11.10,将得到$Born$近似公式:

其中,$D$为$\delta K(\rm \textbf x)$的定义域,即$\delta K(\rm \textbf x)$在$D$外为null。对于一维介质中垂直于纵轴的入射平面波,方程11.15与坐标$x$和$y$无关;即

并有

其中,$z_r$代表检波器深度,$z_s$代表震源深度。假设震源位于表面$z_s=0$。公式11.6变为:

其中,

对于反射问题(即$z_r \leq z_1$),对11.18积分可得

公式11.20中的指数项控制了波的传播。我们可以看出其只受背景介质$V_{P0}$的影响,而与$V_{P1}$无关。换句话说,$Born$解以背景介质的速度在嵌入体中传播。我们稍后将讨论这些问题对Born解的准确性的影响。
11.9所示图的准确压力场可在经典的声学书中看到。所以,我们这里只给出最后的结果:

与公式11.20所示的$Born$解不同的是,在准确解中,控制波传播的指数项不但取决于背景介质速度$V_{P0}$,同时还受嵌入体速度$V_{P1}$影响。
11.9_10

在我们对方程11.20中的$Born$解和方程11.21中的精确解进行数值比较之前,我们先介绍我们即将讨论的三个重要参数:①在背景介质中传播的信号的主波长,$\lambda_0=V_{P0}/f_c$,其中$f_c$是源信号的主频;② 在嵌入体中传播的信号的主波长,$\lambda_1=V_{P1}/f_c$③嵌入体与背景介质的相对扰动,$\Delta K=\delta K/K_0$。对于下面的数值实验,源信号的主频为35Hz(见图11.10)。
11.11~13
图11.11、11.12和11.13比较了不同嵌入体厚度和不同速度下的$Born$解和精确解。这些计算中使用的震源信号频谱如图11.10所示。我们可以看到,对于信号在嵌入体中传播的波长远大于平板厚度的情况(即$h\leq \lambda_1$), 当相对扰动小于0.36时,$Born$近似解与精确解几乎一致。然而,随着板的厚度增加,我们可以看到,在$Born$近似解中,嵌入体底部的反射开始偏离精确解,即使$\Delta K$很小。前面我们提到了$Born$近似解与精确解存在差异的原因,即玻恩近似解是以背景介质的速度而不是以平板速度在平板中传播的。如果背景介质的速度与嵌入体速度之差较大,或者波在嵌入体中传播的持续时间足够长,则$Born$近似无法准确预测嵌入体底部反射的到达时间。

我们注意到,$Born$近似中的反射波到时与嵌入体的速度和厚度无关。然而,这种反射的振幅却受到嵌入体与背景介质之间的相对扰动,即$\Delta K$的影响。当$\Delta K=0.55$时,$Born$近似解误差较大。
另外,对于板厚为$25m$,板速度为$3.0 km/s$的情况,$Born$近似解的形状与精确解的形状似乎有很大不同。形状不同是因为玻恩近似假定波在板中以背景速度传播,它预测的两个同相轴之间有一个小的重叠,而在精确解中,基于波以实际速度传播,由于嵌入体中速度快因此两个同相轴完全重叠。
我们还可以看到,当$\Delta K=0.55$,$h \geq100m$时,由于背景介质是均匀的,所以在精确解中能够看到明显的层间多次波,而$Born$近似中却并不明显。如图11.7所示,如果介质是平滑的,$Born$近似无法预测层间多次波;均匀介质是平滑介质的理想情况。
总结来讲,$Born$近似解需要①在嵌入体中传播的波的波长要远远大于嵌入体厚度;②相对扰动要足够小。我们通过一维介质中的平面波解得出了以上结论。在更一般意义上,其由嵌入体的厚度$h$和散射势平均强度的乘积来决定。对于一维情况,$Born$近似有效的标准是:

以上公式对于二维和三维情况下的球面波同样适用,尽管对于二维和三维介质中的球面波这一结论非常难验证。距离和全散射势的平均强度的乘积才是最重要的。成功应用$Born$近似的关键是,背景介质的选择或估计必须足够准确以捕捉实际模型中包含的所有长波长信息,而扰动模型必须基本上由短波组成。图11.8所示为理想情况下,$Born$近似将介质分解为背景介质和扰动介质。