第十一章 逆时偏移

逆时偏移

在本章中，我们推导逆时偏移（RTM）算法。在石油工业生产中，该算法被广泛应用于对盐下反射体的成像。同时，它还被用于迭代地偏移残差LSM和FWI。与其他偏移方法相比，例如第 10 章中讨论的绕射叠加偏移，RTM 考虑了波场中的所有类型的波，如果速度模型足够精确的话其同时包含了初至和多次波。 这可以提升成像分辨率，但它的代价是对偏移速度模型中的误差更敏感。当频率小20Hz时，RTM 是用于偏移盐体等复杂地质构造的首选成像方法。我们将推导出通用成像方程并讨论RTM的实际应用。

11.1 前言

第10章介绍了估计反射率分布（即偏移成像）的绕射叠加偏移（Diffraction-stack migration, DM）方法，其沿着“时差曲线”对地震道的振幅进行加权求和。对于一般速度分布，偏移距-时间坐标中的时差曲线是不规则的、类似双曲线的。标准 DM 存在的问题是它只考虑了初至成像，而无法模拟多次波或者说不能很好地偏移多次波。 事实上，与时差曲线相关的初至能量相对于其他反射可能是相当弱的，所以 DM 成像结果可能质量较差。这一问题在盐下成像时（入射和反射波场defocus）非常明显。
为了从一定程度上解决由于大速度对比导致的散焦（defocusing）问题，Baysal et al. (1983), McMechan (1983)和Whitmore (1983)引入了逆时偏移的概念。相较于使用射线追踪法，逆时偏移中使用波动方程的有限差分解估计格林函数$G(\mathbf x| \mathbf x’)$。如果速度模型是准确已知的，那么所有来自可能的成像点的arrivals都被考虑在内，包括多次波、绕射波和转换波等。另外，如果速度模型足够精确，能够有效提升偏移成像的信噪比，能够对陡倾角成像、也能够对高对比速度下方构造成像。但是，偏移速度模型中小的误差可能会引起偏移成像中较大的defocusing。所以，速度模型精度越高，RTM成像精度越高。下章节将推导通用成像算法下的RTM公式。

11.2 通用成像算法

地震成像和偏移的目标分别是根据记录的地震数据估计地球的慢度和反射率分布。两种算法都源自相同的梯度优化方法，如步骤所示：

1. 损失函数：假设由$\mathbf s$处发出的谐波能量源被地震检波器$\mathbf g$记录。为了简单起见，这里假设声学近似，观测到的压力谱用 $P(\mathbf g| \mathbf s)^{obs}$ 表示，其中未知源小波$w(t)$的频谱用$W(\omega)$表示。压力损失函数$\epsilon$如下所示: $$\epsilon = 1/2 \int_{\omega}d \omega \sum_g \sum_s |\Delta P (\mathbf g | \mathbf s)|^2 \tag{11.1}$$ 其中，$\Delta P(\mathbf g| \mathbf s)=P(\mathbf g| \mathbf s)-P(\mathbf g| \mathbf s)^{obs}$是数据残差，$P(\mathbf g| \mathbf s)$为通过波动方程数值求解得到的预测压力场。假设已知的初始速度模型是一个足够精确到实际数据的模型。对源所有频率进行积分并对源和检波点的所有坐标求和。预测压力场通过有限差分法求解时空域或时频域波动方程得到。
2. 慢度梯度和Frechet导数：使用最速下降法，第$i+1$次迭代的慢度模型$s(\mathbf x)^{(t)}$由以下公式得到： $$s(\mathbf x)^{(t+1)}=s(\mathbf x)^{(t)}-\alpha \gamma(\mathbf x)^{(t)}, \tag{11.2}$$ 其中$\alpha$为迭代步长，$\gamma(\mathbf x)^{(t)}$损失梯度，梯度定义如下： $$\begin{aligned} \gamma(\mathbf{x})^{(i)} & =\frac{\partial \epsilon_i}{\partial s(\mathbf{x})} \\ & =\int_\omega d \omega \sum_g \sum_s \operatorname{Real}\left[\Delta P(\mathbf{g} \mid \mathbf{s})_i \frac{\partial P(\mathbf{g} \mid \mathbf{s})_i^*}{\partial s(\mathbf{x})}\right], \end{aligned}$$ 其中第$i$次数据残差$
   \Delta P(\mathbf{g} \mid \mathbf{s})_i=P(\mathbf{g} \mid \mathbf{s})_i-
   \Delta P(\mathbf{g} \mid \mathbf{s})_{obs}$。算子$Real[ ]$表示求括号中复数的实部。下面表述中实数符号将被删除，因为对频率的求和是对称的，被加数的虚部是在 ω 反对称。
   $P$关于$s$的扰动是Frechet导数（公式10.23），将其带入11.3得到梯度公式： $$\gamma (\mathbf x)^{t} = 2s(\mathbf x) \int_{\omega} d \omega \omega^2 \sum_g \sum_s W(\omega)^* x \Delta P(\mathbf g | \mathbf s)_t G(\mathbf x| \mathbf g)_t^* G(\mathbf x| \mathbf s)_t^* \tag{11.4}$$ 其中，$G(\mathbf x | \mathbf x’)$是第$i$个模型的格林函数。
3. 慢度更新：将公式11.4带入11.2得到慢度更新公式： $$\begin{aligned} s(\mathbf{x})^{(t+1)}= & s(\mathbf{x})^{(t)}-\alpha 2 s(\mathbf{x})^{(t)} \int_\omega d \omega \omega^2 \sum_g \sum_s \\ & \times \Delta P(\mathbf{g} \mid \mathbf{s})_t W(\omega)^* G(\mathbf{x} \mid \mathbf{g})_t^* G(\mathbf{x} \mid \mathbf{s})_t^*, \end{aligned} \tag{11.5}$$ 重复步骤1到3知道收敛到$s(\mathbf x)$。以下五种成像算法都能从以上公式推导得到。

• 反射波偏移。如果平滑的背景慢度模型不产生反射，那么数据残差$\Delta P(\mathbf g | \mathbf s)_t$只包含上行反射，那么第一次迭代的梯度$\delta s(\mathbf{x})=s(\mathbf{x})^{(1)}-s(\mathbf{x})^{(0)}$即为偏移剖面，其为地下反射体分布的近似。如果使用双程波动方程的有限差分解，那么这一近似即为逆时偏移剖面(Baysal et al., 1983; McMechan, 1983; Whitmore,1983)。如果使用单程波格林函数(Gazdag, 1978; Claerbout, 1992)替换公式11.5的双程波格林函数，那么公式11.5变为单程波偏移；如果使用渐近的基于射线的格林函数，其变为到绕射叠加偏移。这种渐近格林函数仅能够模拟初至反射波，但是高斯束法(Hill, 1990; Hill, 2001; Nowack et al., 2003)能够模拟多次波。
• 最小二乘偏移。如果平滑背景慢度在每次迭代中都是相同的，即$G\left(\mathbf{x} \mid \mathbf{x}^{\prime}\right)_1=G\left(\mathbf{x} \mid \mathbf{x}^{\prime}\right)_2=G\left(\mathbf{x} \mid \mathbf{x}^{\prime}\right)_3=\ldots$最后的reflectivity image即为最小二乘偏移剖面（LSM）(Nemeth, 1996; Nemeth et al., 1999; Duquet et al., 2000)。最小二乘偏移又被称为线性化Born反演(Lailly, 1983; Tarantola, 1984b)，15-18章内容。
• 全波形反演。如果格林函数中的高波数反射率和平滑慢度模型在每次迭代后都更新，那么即为全波形反演 (FWI) tomo（arantola 和 Valette，1982；Tarantola，1984a 和 1987；Mora，1987 和 1989），也被称为非线性反演成像（Pica 等，1990）。第 19 章到第 23 章介绍了理其论和示例。
• 旅行时层析成像。如果数据残差除地震道之外的数据，例如旅行时间残差$\tau$，那么 Fréchet 导数的调整导致波动方程层析成像 (Woodward, 1992) 或波动方程
  走时反演 (Luo and Schuster, 1991a, 1991b). 详见第 20 至 21 章。
• 广义FWI。广义FWI算法以以下公式为目标函数： $$\epsilon=\frac{1}{2} \Delta \mathbf{d}^{\dagger} \Delta \mathrm{T}^2 \Delta \mathrm{d} \tag{11.6}$$ 梯度为： $$\frac{\partial \epsilon}{\partial s(\mathbf{x})}=\frac{\partial \Delta \mathbf{d}^{\dagger}}{\partial s(\mathbf{x})} \Delta \mathbf{T}^2 \Delta \mathbf{d}+\Delta \mathbf{d}^{\dagger} \frac{\partial \Delta \mathbf{T}}{\partial s(\mathbf{x})} \Delta \mathrm{T} \Delta \mathbf{d}$$ 为方便期间，假定矢量和矩阵均为实数，$\Delta \mathbf T$是旅行时残差的对角矩阵。如果$\Delta \mathbf{d}=\mathbf{L} \Delta \mathbf{m}$，此时梯度变为： $$\begin{aligned} \frac{\partial \epsilon}{\partial s(\mathbf{x})}= & \overbrace{\frac{\partial \Delta \mathbf{m}^{\dagger}}{\partial s(\mathbf{x})} \mathbf{L}^{\dagger} \Delta \mathbf{T}^2 \Delta \mathbf{d}}^{\text {small-residual FWI }}+\overbrace{\Delta \mathbf{m}^{\dagger} \frac{\partial \mathbf{L}^{\dagger}}{\partial s(\mathbf{x})} \Delta \mathbf{T}^2 \Delta \mathbf{d}}^{\text {large-residual FWI }} \\ & +\overbrace{\Delta \mathbf{d}^{\dagger} \frac{\partial \Delta \mathbf{T}}{\partial s(\mathbf{x})} \Delta \mathbf{T} \Delta \mathbf{d}}^{\text {WT }} \end{aligned} \tag{11.8}$$ 其中，WT表示将在第20章中讨论的波动方程旅行时反演，其用于更新速度模型的低波数分量。小残差 FWI 梯度类似于经典FWI 并且对小数据残差有效。 正如第19到22章所讨论，它可以同时更新速度模型的低波数和高波数分量但是偏向于更深处的高波数分量。当存在大量数据残差时，大残差 FWI 梯度用于高斯牛顿法（参见Gill et al.，1981 和Exercise 12) 。这一流程等同于大残差的高斯-牛顿法应用于像域反演，如第26章所讨论。

当$i=0$时，整理公式15可得到慢度残差：

$$\begin{aligned} m(\mathbf{x})^{\mathrm{mig}}= & s(\mathbf{x})^{(1)}-s(\mathbf{x})^{(0)}=-\alpha 2 s(\mathbf{x}) \int_\omega d \omega \omega^2 \sum_g \sum_s \\ & \times \Delta P(\mathbf{g} \mid \mathbf{s}) W(\omega)^* G(\mathbf{x} \mid \mathbf{g})^* G(\mathbf{x} \mid \mathbf{s})^*, \end{aligned} \tag{11.9}$$

其中，$m(\mathbf x)^{mig}$，即偏移剖面，近似了反射体分布。第 10.1.5 节中提供该公式在绕射叠加偏移中的解释。接下来的两节根据广义绕射偏移和逆时偏移算法解释公式 11.9。

11.2.1 逆时偏移=广义绕射叠加偏移

方程 11.9 描述的偏移图像是纯实数，因此对于一个炮集，复共轭可以得到

$$\begin{aligned} m(\mathbf{x})^{\mathrm{mig}} & =-2 s(\mathbf{x}) \int_\omega d \omega \sum_g \omega^2 W(\omega) \Delta P(\mathbf{g} \mid \mathbf{s})^* G(\mathbf{x} \mid \mathbf{g}) G(\mathbf{x} \mid \mathbf{s}), \\ & =2 \int_t d t \sum_g \overbrace{\ddot{d}(\mathbf{g}, t \mid \mathbf{s}, 0)}^{\text {reflection data }} \overbrace{g(\mathbf{x}, t \mid \mathbf{g}, 0) \star g(\mathbf{x}, t \mid \mathbf{s}, 0)}^{\text {migration kernel }}, \end{aligned}$$

其中，为了符号方便，$W(\omega)=1/2 \pi, \Delta p=(\mathbf g, t | \mathbf s, 0) \rightarrow d(\mathbf g, t | \mathbf s, 0), s(\mathbf x) \rightarrow 1, w(t) =\delta(t)$。
公式11.10与广义绕射叠加偏移（GDM）等价，偏移图像域逆时偏移相同（Schuster, 2002)。如果格林函数只考虑直达波，（即$g(\mathbf{x}, t \mid \mathbf{s}, 0) \approx \delta(t-|\mathbf{x}-\mathbf{s}| / c) /|\mathbf{x}-\mathbf{s}|$），GDM公式变为由公式10.47所描述的DM方法。如图11.1a所示，偏移核$g(\mathbf{x}, t \mid \mathbf{g}, 0) \star g(\mathbf{x}, t \mid \mathbf{s}, 0)=\delta(t-|\mathbf{x}-\mathbf{g}| / c-|\mathbf{x}-\mathbf{s}| / c) /|\mathbf{x}-\mathbf{g}||\mathbf{x}-\mathbf{s}|$对于在$\mathbf x=(0,d)$的特定成像点在$x_g—t$空间为一条双曲线。相反，如果格林函数是通过有限差分求解的波动方程，偏移核在数据空间中为多条双曲线，如图 11.1b所示。
与图 11.1a 中的初至偏移核不同，只有当trail成像点与实际散射位置十分接近时，$x_g—t$空间中预测的多个时差曲线（彩色实线和虚线）与实际时差（黑色双曲线）的点积才可以忽略不计。 因此，多次波偏移图像的分辨率应该比初至波偏移更高。
GDM图像的计算与DM的计算相类似，对于一个特定的trail成像点$\mathbf x$，$m(\mathbf x)^{mig}$是通过炮集记录与偏移核之间的点积$g(\mathbf{x}, t \mid \mathbf{g}, 0) \star g(\mathbf{x}, t \mid \mathbf{s}, 0)$。与DM不同的是，GDM核同时考虑了初至波和所有的多次散射波。详情见第13章。

11.2.2 逆时偏移

公式11.9中的梯度可以整理为一个更易理解的形式：

$$\begin{aligned} & m(\mathbf{x})^{\mathrm{mig}}=-\int_\omega d \omega \omega^2 \sum_g[W(\omega) G(\mathbf{x} \mid \mathbf{s})]^*\left[G(\mathbf{g} \mid \mathbf{x})^* \tilde{d}(\mathbf{g} \mid \mathbf{s})\right] \\ & =\overbrace{[w(t) \star g(\mathbf{x}, t \mid \mathbf{s}, 0)]}^{D(\mathbf{x}, t)=\text { downgoing source }} \\ & \left.\otimes \overbrace{\left[\sum_g g(\mathbf{g},-t \mid \mathbf{x}, 0) \star \ddot{d}(g, t \mid \mathbf{s}, 0)\right]}^{U(\mathbf{x}, t)=\text { backpropagated upgoing reff. }}\right|_{t=0}, \\ & =\sum_t D(\mathbf{x}, t) U(\mathbf{x}, t) \\ & \end{aligned}$$

方便起见，取$\mathbf s(\mathbf x)$为$1/2 \pi$，假设只有一个源$\mathbf s$，$\tilde{d}(\mathbf{g} \mid \mathbf{s})=\Delta P(\mathbf{g} \mid \mathbf{s})$表示只含有上行反射波的地震道数据（已经去掉了直达波和潜波）。这一方程表明， 在trail成像点$\mathbf x$的叠前RTM $m(\mathbf x)^{\text {mig}}$通过下行震源波场$D(\mathbf x, t)$与反传反射波$U(\mathbf x, t)$的零延迟时间互相关计算得到(McMechan, 1983; Whitmore, 1983)。
$D(\mathbf{x}, t)=w(t) \star g(\mathbf{x}, t \mid \mathbf{s}, 0)$表示下行波场（在均质介质中，源波场是严格的下行波，但在非均质介质中它由上行波和下行波组成。 尽管如此，我们还是将其称为“下行”源波。），反传的反射上行波为$U(\mathbf{x}, t)=\sum_g g(\mathbf{g},-t \mid \mathbf{x}, 0) \star \ddot{d}(\mathbf{g}, t \mid \mathbf{s}, 0)$。对于脉冲源，反传的反射只有在$\mathbf x=\mathbf x_{scatt}$时才会与$D(\mathbf x, t)$相遇。理想情况下，这意味着唯一的非零点积$\sum_t U(\mathbf{x}, t) D(\mathbf{x}, t)$出现在散射点的位置。然而在实际应用中，由于震源子波的有限带宽、有限的源-接收器孔径、稀疏的检波器和源采样间隔、偏移速度不准确、多次波和转换波，以及反向传播的反射波场作为散焦波场出现在散射体下方等情况，会导致偏移伪影的出现。
可以在RTM中加预处理项$I(\mathbf x)$来对较弱的偏移振幅进行加权以提升较弱的振幅：

$$m(\mathbf{x})^{\mathrm{mig}}=\frac{\sum_t D(\mathbf{x}, t) U(\mathbf{x}, t)}{I(\mathbf{x})} \tag{11.12}$$

这一项可以视为对Hessian逆的近似（见第11.4章节）且可以补偿减弱的下行源波场和由于几何扩展和衰减引起的上行反射(Beydoun and Mendes, 1989)。它也可以用来对震源子波反褶积，这种情况下通常在时间域用加权的反褶积滤波器或者带通的矩阵来近似Hessian的逆，而非进行频域的划分(Hu et al., 2001; Yu et al., 2006)。如果有多个炮集，还需要对源进行叠加来获得叠加RTM偏移剖面。当然，RTM 成像的质量取决于偏移速度的准确性，因为较大的速度错误将导致由不同炮集叠加得到的$m(\marthbf x)^{mig}$非常模糊。

11.3 逆时偏移数值示例

11.4 逆时偏移实际应用

11.5 总结