三维黏声最小二乘逆时偏移方法模型研究

图1 三维黏滞声介质交错网格示意

表1 黏滞声波波场分量和黏弹性参数的空间位置

网格点	1	2	3	4
黏声波波场分量与对应参数	p, M, τ^P	u, ρ^-1	l, ρ^-1	w, ρ^-1

新窗口打开| 下载CSV

参考图1的参数差分示意图与表1的对应的参数位置与波场分量,对上述三维方程差分处理。具体的差分格式为

(4) $\{\begin{array}{l} u_{i + 1 / 2, j, k}^{m} = u_{i + 1 / 2, j, k}^{m - 1} - \frac{Δ t}{Δ xρ} \{\overset{N}{\sum_{n = 1}} a_{n} [p_{i + 1 / 2 + (2 n - 1) / 2, j, k}^{m - 1 / 2} - p_{i + 1 / 2 - (2 n - 1) / 2, j, k}^{m - 1 / 2}]\}, \\ l_{i, j + 1 / 2, k}^{m} = l_{i, j + 1 / 2, k}^{m - 1} - \frac{Δ t}{Δ yρ} \{\overset{N}{\sum_{n = 1}} a_{n} [p_{i, j + 1 / 2 + (2 n - 1) / 2, k}^{m - 1 / 2} - p_{i, j + 1 / 2 - (2 n - 1) / 2, k}^{m - 1 / 2}]\}, \\ w_{i, j, k + 1 / 2}^{m} = w_{i, j, k + 1 / 2}^{m - 1} - \frac{Δ t}{Δ zρ} \{\overset{N}{\sum_{n = 1}} a_{n} [p_{i, j, k + 1 / 2 + (2 n - 1) / 2}^{m - 1 / 2} - p_{i, j, k + 1 / 2 - (2 n - 1) / 2}^{m - 1 / 2}]\}, \\ u l_{i + 1 / 2, j, k}^{m} = \frac{Δ t}{Δ xρ} \{\overset{N}{\sum_{n = 1}} a_{n} [p_{i + 1 / 2 + (2 n - 1) / 2, j, k}^{m - 1 / 2} - p_{i + 1 / 2 - (2 n - 1) / 2, j, k}^{m - 1 / 2}]\}, \\ l l_{i, j + 1 / 2, k}^{m} = \frac{Δ t}{Δ yρ} \{\overset{N}{\sum_{n = 1}} a_{n} [p_{i, j + 1 / 2 + (2 n - 1) / 2, k}^{m - 1 / 2} - p_{i, j + 1 / 2 - (2 n - 1) / 2, k}^{m - 1 / 2}]\}, \\ w l_{i, j, k + 1 / 2}^{m} = \frac{Δ t}{Δ zρ} \{\overset{N}{\sum_{n = 1}} a_{n} [p_{i, j, k + 1 / 2 + (2 n - 1) / 2}^{m - 1 / 2} - p_{i, j, k + 1 / 2 - (2 n - 1) / 2}^{m - 1 / 2}]\}, \\ p_{i, j, k}^{m + 1} = 2 p_{i, j, k}^{m} - p_{i, j, k}^{m - 1} + (1 + τ_{i, j, k}^{m}) \{\frac{Δ t^{2} v^{2}}{Δ xρ} \{\overset{N}{\sum_{n = 1}} a_{n} [u l_{i + (2 n - 1) / 2, j, k}^{m} - u l_{i - (2 n - 1) / 2, j, k}^{m}]\} + \\ \frac{Δ t^{2} v^{2}}{Δ yρ} \{\overset{N}{\sum_{n = 1}} a_{n} [l l_{i, j + (2 n - 1) / 2, k}^{m} - l l_{i, j - (2 n - 1) / 2, k}^{m}]\} + \frac{Δ t^{2} v^{2}}{Δ zρ} \{\overset{N}{\sum_{n = 1}} a_{n} [w l_{i, j, k + (2 n - 1) / 2}^{m} - w l_{i, j, k - (2 n - 1) / 2}^{m}]\}\}, \\ R_{i, j, k}^{m + 1 / 2} = \frac{1}{(1 + \frac{Δ t}{2 σ_{i, j, k}})} \{(1 - \frac{Δ t}{2 σ_{i, j, k}}) R_{i, j, k}^{m - 1 / 2} - \frac{ρ v^{2} Δ t τ_{i, j, k}^{m}}{2 τ_{σ_{i, j, k}}^{m}} \{\frac{1}{Δ x} \{\overset{N}{\sum_{n = 1}} a_{n} [u l_{i + (2 n - 1) / 2, j, k}^{m} - u l_{i - (2 n - 1) / 2, j, k}^{m}]\} + \\ \frac{1}{Δ y} \{\overset{N}{\sum_{n = 1}} a_{n} [l l_{i, j + (2 n - 1) / 2, k}^{m} - w l_{i, j - (2 n - 1) / 2, k}^{m}]\} + \frac{1}{Δ z} \{\overset{N}{\sum_{n = 1}} a_{n} [w l_{i, j, k + (2 n - 1) / 2}^{m} - w l_{i, j, k - (2 n - 1) / 2}^{m}]\}\} 。 \end{array}$

上式中:u、l、w分别为速度场在x、y、z方向的分量,ul、ll、wl分别是中间变量,P代表应力场,R代表记忆变量。

文中采用的三维情况下的稳定性条件为

(5) $c_{u} Δ t [\overset{m}{\sum_{n = 1}} a_{n} {(- 1)}^{n - 1}] \sqrt[]{\overset{3}{\sum_{i = 1}} \frac{1}{Δ x_{i}^{2}}} \leq 1,$

其中对于Δx_i,当取1、2或3时分别对应Δx、Δy或者Δz;c_u为最高的高频极限速度。

当所有方向的空间网格间距都相等时,条件(5)简化为:

(6) $c_{u} \frac{Δ t}{h} [\overset{m}{\sum_{n = 1}} a_{n} {(- 1)}^{n - 1}] \sqrt[]{3} \leq 1 。$

1.3 三维黏声最小二乘伴随算子的构建

基于Albert Tarantola对数据空间标量乘法的定义,可以得到标量乘法:

(7) $< d_{1} (x, t), d_{2} (x, t) >_{D} = \int_{x_{1}}^{x_{2}} d x \int_{t_{1}}^{t_{2}} [d_{1} (x, t) d_{2} (x, t)] d t 。$

在对地震波正演模拟计算时,首先是将波动方程进行线性化处理,在求解黏声逆时偏移伴随算子时,根据式(7)定义的方程,得到伴随算子满足的标量乘法:

(8) $< p^{*}, Lp >_{D} = < L^{*} p^{*}, p >_{D},$

式中:L是将波动方程线性化的算子,L^*是L的伴随算子,同理,p^*就是波动方程线性化算子L伴随算子L^*下的正演数据的波场值。将最原始的方程代入到式(8)中,得到

(9) $< L^{*} p^{*}, p >_{D} = \int_{x_{1}}^{x_{2}} d x \int_{t_{1}}^{t_{2}} \{p^{*} (x, t) [\frac{1}{v^{2}} \frac{\partial p^{2}}{\partial t^{2}} - \frac{\partial}{\partial t} [(1 - \overset{I}{\sum_{i = 1}} (1 - \frac{τ_{εi}}{τ_{σi}}) e^{- t / t_{σi}}) H (t)] \times ρ \nabla \cdot (\frac{1}{ρ} \nabla p)]\} d t,$

为了求解上述方程,还需要给出相应的初始边界条件,初始边界条件的表达式为

(10) $\begin{array}{l} p (x, t) |_{t = t_{1}} = 0, \end{array}$

(11) $\begin{array}{l} \frac{\partial p}{\partial t} (x, t) |_{t = t_{1}} = 0 。 \end{array}$

则方程(9)被整理为

(12) $\begin{array}{l} \int_{x_{1}}^{x_{2}} d x \int_{t_{1}}^{t_{2}} p^{*} (x, t) \frac{1}{v^{2}} \frac{\partial p^{2}}{\partial t^{2}} d t = \\ \int_{x_{1}}^{x_{2}} d x [p^{*} \frac{1}{v^{2}} \frac{\partial p}{\partial t} |_{t_{1}}^{t_{2}} - \frac{\partial p^{*}}{\partial t} \frac{1}{v^{2}} p |_{t_{1}}^{t_{2}} + \int_{t_{1}}^{t_{2}} p \frac{1}{v^{2}} \frac{\partial^{2} p^{*}}{\partial t^{2}} d t] 。 \end{array}$

伴随波场满足的边界条件为

(13) $\begin{array}{l} p^{*} (x, t) |_{t = t_{2}} = 0, \end{array}$

(14) $\begin{array}{l} \frac{\partial p^{*}}{\partial t} (x, t) |_{t = t_{2}} = 0 。 \end{array}$

通过推导可得到三维黏声波伴随算子满足的波动方程(具体推导见附录A):

(15) $\frac{1}{v^{2}} \frac{\partial^{2} p^{*}}{\partial t^{2}} = \nabla \cdot \{\frac{1}{ρ} \nabla (ρ p^{*} (t) * \frac{\partial}{\partial t} [(1 - \overset{I}{\sum_{i = 1}} (1 - \frac{τ_{εi}}{τ_{σi}}) e^{- t / t_{σi}}) H (t)])\} + f 。$

三维黏声波伴随算子可通过下式所示的离散格式进行求解:

(16) $\{\begin{array}{l} u_{i + 1 / 2, j, k}^{m} = u_{i + 1 / 2, j, k}^{m - 1} - \frac{Δ t}{Δ xρ} \{\overset{N}{\sum_{n = 1}} a_{n} [p_{i + 1 / 2 + (2 n - 1) / 2, j, k}^{* m - 1 / 2} - p_{i + 1 / 2 - (2 n - 1) / 2, j, k}^{* m - 1 / 2}]\}, \\ l_{i, j + 1 / 2, k}^{m} = l_{i, j + 1 / 2, k}^{m - 1} - \frac{Δ t}{Δ yρ} \{\overset{N}{\sum_{n = 1}} a_{n} [p_{i, j + 1 / 2 + (2 n - 1) / 2, k}^{* m - 1 / 2} - p_{i, j + 1 / 2 - (2 n - 1) / 2, k}^{* m - 1 / 2}]\}, \\ w_{i, j, k + 1 / 2}^{m} = w_{i, j, k + 1 / 2}^{m - 1} - \frac{Δ t}{Δ zρ} \{\overset{N}{\sum_{n = 1}} a_{n} [p_{i, j, k + 1 / 2 + (2 n - 1) / 2}^{* m - 1 / 2} - p_{i, j, k + 1 / 2 - (2 n - 1) / 2}^{* m - 1 / 2}]\}, \\ u l_{i + 1 / 2, j, k}^{m} = \frac{Δ t}{Δ xρ} \{\overset{N}{\sum_{n = 1}} a_{n} [p_{i + 1 / 2 + (2 n - 1) / 2, j, k}^{* m - 1 / 2} - p_{i + 1 / 2 - (2 n - 1) / 2, j, k}^{* m - 1 / 2}]\}, \\ l l_{i, j + 1 / 2, k}^{m} = \frac{Δ t}{Δ yρ} \{\overset{N}{\sum_{n = 1}} a_{n} [p_{i, j + 1 / 2 + (2 n - 1) / 2, k}^{* m - 1 / 2} - p_{i, j + 1 / 2 - (2 n - 1) / 2, k}^{* m - 1 / 2}]\}, \\ w l_{i, j, k + 1 / 2}^{m} = \frac{Δ t}{Δ zρ} \{\overset{N}{\sum_{n = 1}} a_{n} [p_{i, j, k + 1 / 2 + (2 n - 1) / 2}^{* m - 1 / 2} - p_{i, j, k + 1 / 2 - (2 n - 1) / 2}^{* m - 1 / 2}]\}, \\ p_{i, j, k}^{* m + 1} = 2 p_{i, j, k}^{* m} - p_{i, j, k}^{* m - 1} + (1 + τ_{i, j, k}^{m}) \{\frac{Δ t^{2} v^{2}}{Δ xρ} \{\overset{N}{\sum_{n = 1}} a_{n} [u l_{i + (2 n - 1) / 2, j, k}^{m} - u l_{i - (2 n - 1) / 2, j, k}^{m}]\} + \\ \frac{Δ t^{2} v^{2}}{Δ yρ} \{\overset{N}{\sum_{n = 1}} a_{n} [l l_{i, j + (2 n - 1) / 2, k}^{m} - l l_{i, j - (2 n - 1) / 2, k}^{m}]\} + \frac{Δ t^{2} v^{2}}{Δ zρ} \{\overset{N}{\sum_{n = 1}} a_{n} [w l_{i, j, k + (2 n - 1) / 2}^{m} - w l_{i, j, k - (2 n - 1) / 2}^{m}]\}\}, \\ R_{i, j, k}^{m + 1 / 2} = \frac{1}{(1 + \frac{Δ t}{2 σ_{i, j, k}})} \{(1 - \frac{Δ t}{2 σ_{i, j, k}}) R_{i, j, k}^{m - 1 / 2} - \frac{ρ v^{2} Δ t τ_{i, j, k}^{m}}{2 τ_{σ_{i, j, k}}^{m}} \{\frac{1}{Δ x} \{\overset{N}{\sum_{n = 1}} a_{n} [u l_{i + (2 n - 1) / 2, j, k}^{m} - u l_{i - (2 n - 1) / 2, j, k}^{m}]\} + \\ \frac{1}{Δ y} \{\overset{N}{\sum_{n = 1}} a_{n} [l l_{i, j + (2 n - 1) / 2, k}^{m} - w l_{i, j - (2 n - 1) / 2, k}^{m}]\} + \frac{1}{Δ z} \{\overset{N}{\sum_{n = 1}} a_{n} [w l_{i, j, k + (2 n - 1) / 2}^{m} - w l_{i, j, k - (2 n - 1) / 2}^{m}]\}\} 。 \end{array}$

1.4 黏声LSRTM实现流程

最小二乘偏移理论的实质就是在地震反问题中求解最优解的过程。在最小二乘意义下最小化模拟数据与观测数据的误差方程,然后使用扰动模型更新初始模型以获得最终的成像结果。笔者主要是在了解梯度法的基础上,利用预条件共轭梯度法求解最小二乘最优解。最小二乘逆时偏移可以使用下面的迭代公式进行计算:

(17) $m^{(k + 1)} = m^{(k)} + β^{(k)} d k^{(k)},$

式中:m是模型参数;dk⁽^k⁾是经过第k次修正之后的梯度方向;β⁽^k⁾是迭代求解的步长。梯度和步长由式(18)~(21)求得:

(18) $\begin{array}{l} g^{(k)} = L^{T} [L m^{(k)} - d_{obs}], \end{array}$

(19) $\begin{array}{l} α^{(k)} = \frac{[g^{(k)}]^{T} [C g^{(k)}]}{[g^{(k - 1)}]^{T} [C g^{(k - 1)}]}, \end{array}$

(20) $\begin{array}{l} d k^{(k)} = C k^{(k)} + βd k^{(k - 1)}, \end{array}$

(21) $\begin{array}{l} β^{(k)} = \frac{[d k^{(k)}]^{T} [Cd k^{(k)}]}{[L k^{(k - 1)}]^{T} [Ld k^{(k - 1)}]} 。 \end{array}$

式中:g⁽^k⁾和α⁽^k⁾是由最速下降法求解获得的梯度方向与步长。假设步长α满足抛物线的性质,利用插值法求解。L是线性正演算子;T表示矩阵的转置;d_obs是观测数据;L^T是偏移算子;C表示预条件算子,通过该预条件算子可以使反演问题更稳定、收敛速度更快。在文中,我们的预条件算子是使用的震源波场照明方法。

图2为黏声波最小二乘逆时偏移成像的技术路线图,图中主要分为三个部分:第一部分是黏声波正演模拟,第二部分是黏声波逆时偏移,第三部分是最小二乘逆时偏移迭代的过程。

图2

图2 黏声最小二乘逆时偏移成像技术路线流程

从流程图上可以看到正演部分的核心技术是构造高精度的正演算子。逆时偏移部分的关键是稳定高效的逆时传播算子的构建,最小二乘逆时偏移部分的关键就是稳定的反偏移算子的构建,解决计算量问题的关键就是引入合理极性编码技术(见附录B)。地震波在黏弹(声)介质中反向传播时,会出现高频不稳定现象,通常采用加正则化算子方法、低频滤波方法和伴随算子方法。本论文采用黏声波伴随算子的方法克服反向传播的高频不稳定,因此,对于本论文提出的三维黏声最小二乘逆时偏移来说,最关键的环节就是黏声波伴随算子的引入和计算。在迭代过程中,为了简化,假设品质因子参数和密度参数都是准确的,也就是说,品质因子和密度的扰动为零,而p参数存在扰动。

2 模型试算

2.1 Marmousi模型

与常规逆时偏移方法一样,黏声波最小二乘逆时偏移成像也需要分为三个部分:第一部分是震源波场的正向延拓,第二部分是检波场的反向延拓,第三部分是利用成像条件对其进行偏移成像。由于黏声最小二乘逆时偏移成像方法需要在避开不稳定性问题的同时对黏滞性的吸收衰减进行补偿,故其关键技术就是如何构建稳定的逆时传播算子。笔者主要运用伴随算子对其进行求解,以确保逆时传播算子的稳定性。为了验证方法的准确性,先运用二维黏声最小二乘逆时偏移方法对Marmousi模型进行方法测试,通过声波逆时偏移与黏声波逆时偏移成像效果对比,来验证黏声波方程对地震波在构造深部高波数能量的补偿作用,通过声波与黏声波最小二乘逆时偏移成像对比,来凸显在反演框架下的地震偏移成像的优势。

图3所用的Marmousi模型大小为4 416 m×2 244 m,横向采样点数为369,纵向采样点数为188,纵横向空间采样间隔12 m。图3a是真实的速度模型,图3b为偏移所用的Q模型,图3c为平滑P波速度模型,图3d为慢度扰动模型。均一的密度为 2.0 g/cm³。地表以下12 m处等间隔放炮148炮,检波器总数为369个,为了确保数值求解的稳定性,采样步长为0.5 ms。图4是对应的成像结果,其中图4a是声波逆时偏移的结果,图4b是黏声波逆时偏移的结果,图4c是声波最小二乘逆时偏移第30次迭代的结果,图4d是黏声波最小二乘逆时偏移第30次迭代的结果。图4a与图4b的偏移成像结果表明,在构造的浅部两种偏移方法都能对构造进行精确的成像,但是在深层,常规声波逆时偏移得到的成像结果同相轴变得模糊不清,分辨率很低,深部能量较弱,而黏声波逆时偏移由于补偿了地震波的高波数成分,所以获得了成像结果清晰、能量均衡的偏移剖面。但是黏声波逆时偏移的偏移剖面上同样存在着低频噪声、对深部储层成像分辨率不高且振幅不均衡的问题。为了解决这些问题,将逆时偏移成像结果作为最小二乘反演的第一步输入,通过求取相应的梯度和步长来不断更新成像结果,最后实现对深部构造成像的目的。图4c展示的是常规声波最小二乘逆时偏移的成像结果。与图4d对比,发现以声波LSRTM成像为理论的成像方法,不仅不能准确校正介质黏滞性引起的速度频散问题,而且还致使反射界面的空间归位不准确,反射界面空间的波场特征、振幅变化和反射系数不准确,引入成像假象。由于该方法没有对地震波的衰减进行补偿,所以其成像振幅与真实的慢度扰动相比,远小于后者。图4d是黏声波LSRTM迭代30次的偏移成像结果,当迭代30次时,成像振幅得到了补偿,能量得到了均衡,衰减了的高波数成分得到了有效恢复,压制了逆时偏移中存在的低频噪声,同时校正了由于弱照明区域照明不均衡所引起的模型横向变化问题。图5所示为x=2 760 m处反射系数与成像结果单道对比。从图中可以看出,第一次迭代的成像结果振幅能量较弱,而经过30次迭代后,成像结果的振幅能量与真实成像结果(反射系数)能量接近,而且因黏声介质造成的相位变化也得到了很好的校正。因此,文中的黏声最小二乘逆时偏移成像算法能够得到高质量的保幅成像结果。

图3

图3 Marmousi模型

a—速度模型;b—Q模型;c—偏移速度模型;d—慢度扰动模型

图4

图4 偏移成像结果

a—声波LSRTM第1次迭代;b—黏声LSRTM第1次迭代;c—声波LSRTM第30次迭代;d—黏声波LSRTM第30次迭代

图5

图5 x=2 760 m处反射系数与成像结果单道对比

为了证明本文黏声伴随算子方法对低频滤波方法和添加正则化项方法的优势,给出了两种算法的对比(低频滤波方法与添加正则化项方法本质上是相同的,都是对高频成分进行滤波,只是实现方式不同,这里只给出低频滤波方法的结果),如图6所示。可以看出,采用低频滤波方法损失了部分有效能量,而本文方法振幅的保幅性更好。

图6

图6 黏声最小二乘逆时偏移成像结果

a—本文方法;b—低通滤波方法

最后,为了探索本文方法对偏移速度模型的精度要求,引入随机误差,并对速度模型进行平滑,图7给出了误差分别为3.84%、5.63%、7.52%和10.13%的偏移速度模型,图8表示采用不同精度偏移速度模型的收敛曲线。从图中可以看出,在一定误差范围内,本文方法可以得到很好的收敛,随着误差增加,收敛效果变差。

图7

图7 不同误差的偏移速度模型

a—误差=3.84%;b—误差=5.63%;c—误差=7.52%;d—误差=10.13%

图8

图8 不同误差的迭代收敛曲线

2.2 三维简单模型测试

随着地震勘探难度的增加,二维地震数据处理已经不能满足勘探需求,如今面临的勘探都是三维地震勘探,即地震数据正从二维向三维过渡,如何发展高精度的三维偏移成像方法,已迫在眉睫^[47]。在二维研究的基础上,增加一维,发展到三维,实现了三维黏声最小二乘逆时偏移。

设计图9所示的3层平层模型,该模型的网格个数为81×51×61,其中网格间距为10 m。时间采样间隔为1 ms,地表等间隔放炮200炮,将200炮设计成一个超道集。模型设计为3个水平层(图9)。品质因子模型为一均匀模型,Q值为30。在黏声最小二乘逆时偏移成像中,我们经常默认反射率模型就代表着真实的偏移成像结果。图10为采用震源极性编码技术,将200炮三维黏声正演模拟得到的炮记录组合成一个超级炮,分别进行三维黏声最小二乘逆时偏移成像。采用OpenMP加速方式进行成像,单次迭代耗时75.5 min,20次迭代总共耗时近1 510 min。

图9

图9 平层模型

a—速度模型;b—慢度扰动模型

图10

图10 平层模型不同迭代次数下黏声最小二乘逆时偏移成像结果

a—第1次迭代;b—第20次迭代

图10a是迭代了1次的成像结果,即相当于常规的三维黏声最小二乘逆时偏移成像结果,图10b是迭代20次的偏移成像结果。通过对比偏移成像剖面发现,当经过1次迭代的时候,其偏移剖面上存在着大量的低频噪声,构造特征模糊不清,高陡构造成像质量较差;经过不断的迭代,当迭代20次的时候,低频噪声得到了有效的压制,振幅能量得到了有效的均衡,构造剖面清晰,实现了真正的三维保幅逆时偏移成像。

2.3 三维盐丘模型测试

为了测试三维黏声最小二乘逆时偏移对复杂构造的成像效果及其该方法在三维成像效果中的优势,下面对三维盐丘模型进行了测试。模型网格数为338×338×201,时间采样间隔为1 ms,空间采样间隔为20 m,炮点和检波点个数都是600,并且均匀的分布于地表,并且利用极性编码技术将600炮组成一个超道集,构成一个超级炮。盐丘的品质因子根据已知的速度模型,利用李氏经验公式进行求取。图11所示的盐丘模型是从三维盐丘模型中截取的一部分,图12所示的是从三维盐丘模型中截取的盐丘等深度切片模型的成像结果。采用OpenMP并行进行成像,单次迭代耗时160 min。

图11

图11 盐丘模型偏移成像结果

a—盐丘模型(从三维盐丘模型中截取);b—反射系数;c—第1次迭代成像结果;d—第20次迭代成像结果

图12

图12 盐丘等深度切片模型偏移成像结果

a—盐丘等深度切片模型(从三维盐丘模型中截取);b—反射系数;c—第1次迭代成像结果;d—第20次迭代成像结果

从图11c和12c偏移结果可以发现,在模型的浅部区域,第1次迭代结果就取得了分辨率较高的成像结果,但是截取部分的成像品质非常不理想,存在明显的低频噪声;当迭代到20次时,已经可以比较好地对构造进行成像,均衡了两侧的振幅,解决了弱照明区域成像质量差的问题。不仅如此,应用三维黏声最小二乘逆时偏移成像技术,还可以补偿地下介质对地震波的吸收衰减,尤其是地震波的高波数成分。

3 结论与认识

1) 常规逆时偏移方法存在固有的低频噪声,在野外实际地震资料采集的过程中,观测系统的不规则性也会对偏移成像结果产生重大影响。因此,为了获得高品质的成像结果,需要将成像引入到反演的框架下,发展黏声最小二乘逆时偏移成像方法。该方法不仅能补偿地震波衰减的能量,压制低频噪声和均衡振幅,还能获得高保真度、高信噪比和高分辨率的偏移成像剖面。

2) 基于伴随状态理论,推导了3D黏声波最小二乘逆时偏移的伴随算子以及伴随算子满足的波动方程;利用黏声最小二乘逆时偏移对Marmousi模型进行测试,验证了方法的正确性;然后,利用简单的三层平层模型和三维盐丘模型截取的一部分来测试方法的有效性,最终得到了高品质的成像效果。

3) 采用编码技术将多炮数据压缩为一个超道集,三维黏声最小二乘逆时偏移的计算量降为一个单炮的数量级,使得三维黏声最小二乘逆时偏移的实际生产应用成为可能。

4) 将三维黏声最小二乘逆时偏移应用于实际生产时,除了要考虑计算量外,还需要获得较为准确的偏移速度和品质因子。采用黏声全波形反演方法对偏移速度和品质因子进行提取是目前最热门的研究方向,子波提取的准确性也是决定该方法成像效果的决定性因素之一。

附录A 三维黏声最小二乘逆时偏移伴随算子

将伴随算子满足的边界条件代入到式(12)中,则式(12)简化为:

(A1) $\int_{x_{1}}^{x_{2}} d x \int_{t_{1}}^{t_{2}} p^{*} (x, t) \frac{1}{v^{2}} \frac{\partial p^{2}}{\partial t^{2}} d t = \int_{x_{1}}^{x_{2}} d x \int_{t_{1}}^{t_{2}} p \frac{1}{v^{2}} \frac{\partial^{2} p^{*}}{\partial t^{2}} d t,$

为了推导公式的方便,将方程中代表黏滞性的这一项用一个变量进行替代,即:

(A2) $A (t) = \frac{\partial}{\partial t} [(1 - \overset{I}{\sum_{i = 1}} (1 - \frac{τ_{εi}}{τ_{σi}}) e^{- t / t_{σi}}) H (t)] 。$

为了研究的方便,首先以x方向为例进行公式的推导,把对空间函数的求导和褶积运算做展开处理,然后变换积分顺序,得到(A3)式:

(A3) $\begin{array}{l} \int_{x_{1}}^{x_{2}} d x \{- p^{*} A (t - τ) ρ \frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial p (τ)}{\partial x})\} = - p^{*} A (t - τ) \frac{\partial p (τ)}{\partial x} |_{x_{1}}^{x_{2}} + \\ \frac{\partial}{\partial x} (p^{*} A (t - τ) ρ) \frac{1}{ρ} p (τ) |_{x_{1}}^{x_{2}} - \int_{x_{1}}^{x_{2}} d x \{p (τ) \frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial}{\partial x} (p^{*} A (t - τ) ρ))\}, \end{array}$

方程(A1)转化成:

(A4) $\begin{array}{l} \int_{t_{1}}^{t_{2}} d t \int_{x_{1}}^{x_{2}} d x \int_{z_{1}}^{z_{2}} d z \int_{t_{1}}^{t_{2}} d τ \{- p^{*} A (t - τ) ρ \frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial p (τ)}{\partial x})\} = \\ \int_{t_{1}}^{t_{2}} d t \int_{x_{1}}^{x_{2}} d x \int_{z_{1}}^{z_{2}} d z \int_{t_{1}}^{t_{2}} d τ \{- p (τ) \frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial}{\partial x} (p^{*} A (t - τ) ρ))\}, \end{array}$

将 $\int_{t_{1}}^{t_{2}} (\int_{t_{1}}^{t} d τ)$ dt= $\int_{t_{1}}^{t_{2}} (\int_{τ}^{t_{2}} d t)$ dτ式代入到(A4)式中,运用分步积分,化简该式为:

(A5) $\begin{array}{l} \int_{t_{1}}^{t_{2}} d t \int_{x_{1}}^{x_{2}} d x \int_{z_{1}}^{z_{2}} d z \int_{t_{1}}^{t} d τ \{- p^{*} A (t - τ) ρ \frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial p (τ)}{\partial x})\} = \\ \int_{t_{1}}^{t_{2}} d t \int_{z_{1}}^{z_{2}} d z \int_{x_{1}}^{x_{2}} d x \{- p (τ) \frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial}{\partial x} (ρ \int_{t_{1}}^{t_{2}} p^{*} A (t - τ) d τ))\} 。 \end{array}$

设与黏滞性有关项的变量满足如下关系,即

(A6) $A^{*} (t) = A (- t),$

故(A5)式最终化简为:

(A7) $\begin{array}{l} \int_{t_{1}}^{t_{2}} d t \int_{x_{1}}^{x_{2}} d x \int_{z_{1}}^{z_{2}} d z \int_{t_{1}}^{t} d τ \{- p^{*} A (t - τ) ρ \frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial p (τ)}{\partial x})\} = \\ \int_{t_{1}}^{t_{2}} d t \int_{x_{1}}^{x_{2}} d x \int_{z_{1}}^{z_{2}} d z \{- p (τ) \frac{\partial}{\partial x} (\frac{1}{ρ} \frac{\partial}{\partial x} (ρ p^{*} (t) * A^{*} (t)))\}, \end{array}$

同理,y方向的表达式为:

(A8) $\begin{array}{l} \int_{t_{2}}^{t_{3}} d t \int_{x_{2}}^{x_{3}} d x \int_{y_{2}}^{y_{3}} d y \int_{z_{2}}^{z_{3}} d z \int_{t_{1}}^{t} d τ \{- p^{*} A (t - τ) ρ \frac{\partial}{\partial y} (\frac{1}{ρ} \frac{\partial p (τ)}{\partial y})\} = \\ \int_{t_{2}}^{t_{3}} d t \int_{x_{2}}^{x_{3}} d x \int_{y_{2}}^{y_{3}} \int_{z_{2}}^{z_{3}} d z \{- p (τ) \frac{\partial}{\partial y} (\frac{1}{ρ} \frac{\partial}{\partial y} (ρ p^{*} (t) * A^{*} (t)))\}, \end{array}$

同理,z方向的表达式为:

(A9) $\begin{array}{l} \int_{t_{2}}^{t_{3}} d t \int_{x_{2}}^{x_{3}} d x \int_{z_{2}}^{z_{3}} d z \int_{t_{1}}^{t} d τ \{- p^{*} A (t - τ) ρ \frac{\partial}{\partial z} (\frac{1}{ρ} \frac{\partial p (τ)}{\partial z})\} = \\ \int_{t_{2}}^{t_{3}} d t \int_{x_{2}}^{x_{3}} d x \int_{z_{2}}^{z_{3}} d z \{- p (τ) \frac{\partial}{\partial z} (\frac{1}{ρ} \frac{\partial}{\partial z} (ρ p^{*} (t) * A^{*} (t)))\} 。 \end{array}$

把式(13)、(14)、(A7)和式(A8)代入到式(11)中,经过一系列的化简,得到了初始波动方程和它的伴随算子满足的波动方程及其边界条件:

(A10) $\begin{array}{l} L^{*} p^{*} = \frac{1}{v^{2}} \frac{\partial^{2} p^{*}}{\partial t^{2}} - \nabla \cdot (\frac{1}{ρ} \nabla (ρ p^{*} (t) * A^{*} (t))) = f \end{array}$

(A11) $\begin{array}{l} p^{*} (x, t) |_{t = t_{1}} = 0, \end{array}$

(A12) $\begin{array}{l} \frac{\partial p^{*} (x, t) |_{t = t_{1}}}{\partial t} = 0 。 \end{array}$

将中间变量A^*(t)代入式(A10),得到伴随方程的显式表达形式:

(A13) $\frac{1}{v^{2}} \frac{\partial^{2} p^{*}}{\partial t^{2}} = \nabla \cdot \{\frac{1}{ρ} \nabla (ρ p^{*} (t) * \frac{\partial}{\partial t} [(1 - \overset{I}{\sum_{i = 1}} (1 - \frac{τ_{εi}}{τ_{σi}}) e^{- t / t_{σi}}) H (t)])\} + f 。$

附录B 极性编码实现步骤

极性编码技术的实现步骤,可以从以下几点出发:

1)构建超道集。对于固定的采集系统,编码的多震源数据可以写成:

(B1) $d = \overset{S}{\sum_{i = 1}} P_{i} d_{i},$

其中:S是组成超道集的炮集数量;P_i是编码时使用的编码函数,当P_i是编码矩阵时,可以确保 ${P^{T}}_{i}$ P_i是理想的单位矩阵。

假设每一道的正演数据d_i与反射率模型m之间满足

(B2) $d_{i} = L_{i} m,$

其中,L_i是线性正演模型的算子。将式(B2)代入到式(B1),有:

(B3) $d = \overset{S}{\sum_{i = 1}} P_{i} L_{i} m = Lm,$

这里将超道集的正演模型算子定义为:

(B4) $L = \overset{S}{\sum_{i = 1}} P_{i} L_{i} 。$

2)多震源数据的偏移成像。由于式(B4)给出了超道集的正演模型算子,因此根据超道集的正演模型算子可以得到超道集的偏移算子,即

(B5) $L^{T} = \overset{S}{\sum_{i = 1}} {L^{T}}_{i} {P^{T}}_{i},$

则组成超道集的偏移结果可以通过下式给出:

(B6) $\begin{array}{l} m_{mig} = L^{T} d = L^{T} \overset{S}{\sum_{i = 1}} P_{i} L_{i} m, \end{array}$

(B7) $\begin{array}{l} \begin{array}{r} m_{mig} & = L^{T} d = L^{T} \overset{S}{\sum_{i = 1}} P_{i} L_{i} m \\ = \overset{S}{\sum_{j = 1}} {L^{T}}_{j} {P^{T}}_{j} \overset{S}{\sum_{i = 1}} P_{i} L_{i} m \\ = \overset{S}{\sum_{i = 1}} \overset{S}{\sum_{i = 1}} {L^{T}}_{j} {P^{T}}_{j} P_{i} L_{i} m \\ = \overset{正确的成像项}{\overset{︷}{\overset{S}{\sum_{i = 1}} {L^{T}}_{i} L_{i} m}} + \overset{串扰噪声项}{\overset{︷}{\sum_{j \neq i} \overset{S}{\sum_{i = 1}} {L^{T}}_{j} {P^{T}}_{j} P_{i} L_{i} m}} 。 \end{array} \end{array}$

观察上面的偏移成像结果发现,最终的成像结果由两部分组成:第一部分是标准的偏移成像结果,第二部分是超道集之间相互混叠的时候产生的串扰噪声;串扰噪声项幅值大小会依据不同的相位编码函数而不同。

文中进行编码的函数为:

(B8) $P_{i} = rand (\pm 1) 。$

3)极性编码之后的偏移成像计算。在研究极性编码最小二乘逆时偏移方法的时候,根据第一部分定义的超道集,将超道集代入到之前定义的最小二乘梯度公式中,通过化简整理发现,随着最小二乘逆时偏移模型的不断迭代更新,串扰噪声不断地得到压制。

The authors have declared that no competing interests exist.

作者已声明无竞争性利益关系。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Carcione J

Wave propagating in anisotropic linear viscoelastic media:theory and simulated wave-fields

[J]. Geophys. J. In., 1990,101(3):739-742.

DOI:10.1111/gji.1990.101.issue-3 URL [本文引用: 1]

[2]

李振春, 王清振 .

地震波衰减机理及能量补偿研究综述

[J]. 地球物理学进展, 2007,22(4):1147-1152.

随着中、浅层油气勘探程度的不断提高，深层油气藏的勘探已经成为当今勘探研究的主要对象，而深层油气勘探要获得突破关键是用于综合研究的地震剖面要有较高的信噪比和分辨率，深层与浅层能量均衡.因此，国内外对地震波衰减机理和补偿方法进行了大量研究.本文综述了基于岩石物理和散射理论的地震波的衰减机理，重点概括了地震波衰减的影响因素，指出衰减与孔隙度、频率、温度、压力、岩性等因素有关，最后分类介绍了各种地震波衰减的能量补偿方法.

[3]

李金丽, 李振春, 管路平 , 等.

地震波衰减及补偿方法

[J]. 物探与化探, 2015,39(3):456-465.

DOI:10.11720/wtyht.2015.3.04 Magsci [本文引用: 1]

地震波在地下介质中传播时,由于地下介质的吸收衰减作用,地震波的部分弹性能量不可逆转的转化为热能而发生耗散,使得地震波的能量发生衰减,相位产生畸变,降低了地震资料的分辨率和信噪比。为了实现对地震波的吸收衰减进行补偿,国内外许多学者在这方面做了大量的研究工作。笔者简单介绍了地震波的衰减机制和影响地震波衰减的主要因素,重点概括了各种地震波衰减补偿的研究方法,如反Q滤波方法、时频分析方法和反Q偏移方法,分析了各种补偿方法的优缺点。最后用黏声波逆时偏移方法对地震波衰减进行了补偿,并预测了地震波衰减补偿研究的发展趋势。

[4]

刘洋, 李炳秀, 刘财 , 等.

局部时频变换域地震波吸收衰减补偿方法

[J]. 吉林大学学报:地球科学版, 2016,46(2):594-602.

[5]

Bickel S

, Natarajan R

Plane-wave Q deconvolution

[J]. Geophysics, 1985,50(9):1426-1439.

DOI:10.1190/1.1442011 URL [本文引用: 1]

[6]

Gelius L

Inverse Q-filtering:a spectral balancing technique

[J]. Geophysical Prospecting, 1987,35(6):656-667.

DOI:10.1111/gpr.1987.35.issue-6 URL [本文引用: 1]

[7]

Varela C

, Rosa

A L R

,Ulrych

T J

Modeling of attenuation and dispersion

[J]. Geophysics, 1993,58(8):1167-1173.

DOI:10.1190/1.1443500 URL [本文引用: 1]

[8]

Zhang X

, Han L

, Zhang F

, et al.

An inverse Q-filter algorithm based on stable wavefield continuation

[J]. Applied Geophysics, 2007,4(4):263-270.

DOI:10.1007/s11770-007-0040-9 URL [本文引用: 1]

[9]

Hargreaves N

Similarity and the inverse Q filter:some simple algorithms for inverse Q filtering

[J]. Geophysics, 1992,57(7):944-947.

DOI:10.1190/1.1443307 URL [本文引用: 1]

[10]

Wang Y

A stable and efficient approach of inverse Q Filtering

[J]. Geophysics, 2002,67(2):657-664.

DOI:10.1190/1.1468627 URL [本文引用: 1]

[11]

Wang Y

Inverse Q-filter migration

[J]. Geophysics, 2008,73(1):S1-S6.

DOI:10.1190/gpysa7.73.1so_1 URL [本文引用: 1]

[12]

裴江云, 陈树民 .

近地表Q值求取及振幅补偿

[J]. 地球物理学进展, 2001,16(4):18-22.

DOI: Magsci [本文引用: 1]

利用面波衰减的Q模型,反演出近地表的Q值.据此,给出消除振幅受近地表影响的振幅静校正,利用时间剩余静校正方法解决了振幅问题.实际计算结果表明本文方法能有效地去除振幅异常.

[13]

高军, 凌云 , 等.

时频域球面发散和吸收补偿

[J]. 石油地球物理勘探, 1996,31(6):865-866.

[14]

Bickel S

, Natarajan R

Plane-wave Q deconvolution

[J]. Geophysics, 1985,50(9):1426-1439.

DOI:10.1190/1.1442011 URL [本文引用: 1]

[15]

姚振兴, 高星, 李维新 .

用于深度域地震剖面衰减与频散补偿的反Q滤波方法

[J]. 地球物理学报, 2003,46(2):229-233.

根据地震波在非弹性介质中的传播规律，本文提出了一种在深度域地震剖面进行反Q滤波的新方法. 在深度域的Q滤波算子符合地震波在衰减介质中的传播规律，不仅考虑了介质吸收对地震波振幅的影响，而且还保证了所造成的波形畸变满足因果规律，即地震体波具有某种频散性质. 利用该方法处理了2条总长度约60km的某浅海地区叠后深度偏移地震剖面. 合成记录剖面的正演和反演计算结果表明，在对地震剖面进行深度域反Q滤波之后，不仅达到了压缩中、深部目的层子波宽度提高分辨率和信噪比的目的，而且保持了原始地震剖面的整体形态和能量分布.

[16]

Yan H

, Liu

Estimation of inverse Q filtering for prestack reflected PP and converted PS waves

[J]. Applied Geophysics, 2009,6(1):59-69.

DOI:10.1007/s11770-009-0009-y URL [本文引用: 1]

[17]

Dai

, West G

Inverse Q Migration

[C]//Expanded Abstracts of the 64 ^th Annual SEG Meeting,Society of Exploration Geophysicists , 1994: 1418-1421.

[18]

Mittet

, Sollie

Prestack depth migration with compensation for absorption and dispersion

[J]. Geophysics, 1995,60(5):1485-1494.

DOI:10.1190/1.1443882 URL [本文引用: 1]

[19]

Deng

, McMechan

G A

True-amplitude prestack depth migration

[J]. Geophysics, 2007,72(3):S155-S166.

DOI:10.1190/1.2714334 URL [本文引用: 1]

[20]

Zhang

, Zhang

, Zhang H

Compensating for visco-acoustic effects in reverse-time migration

[C]//Expanded Abstracts of the 80 ^th Annual SEG Meeting,Society of Exploration Geophysicists , 2010: 3160-3164.

[21]

Fletcher R

, Nichols

, Cavalca

Wavepath-consistent effective Q estimation for Q compensated reserve-time migration

[C]//Expanded Abstracts of 74 ^th Annual Internat EAGE Mtg , 2012.

[22]

, Li

, Huang J

, et al.

The least-squares reverse time migration for viscoacoustic medium based on a stable reverse-time propagator

[C]//Expanded Abstracts of the 85 ^th Annual SEG Meeting,Society of Exploration Geophysicists , 2015a: 3977-3980.

[23]

, Li

, Huang

, et al.

The application of pseudospectral method and a stable reverse-time propagator for viscoacoustic RTM

[C]//Expanded Abstracts of the 85 ^th Annual SEG Meeting,Society of Exploration Geophysicists , 2015b: 3977-3980.

[24]

, Huang

, Li

, et al.

A new reverse time migration method for viscoacoustic VTI medium

[C]//77 ^th EAGE Conference and Exhibition,Expanded Abstracts , 2015c.

[25]

, Li

, Huang

, et al.

Viscoacoustic anisotropic full waveform inversion

[J]. Journal of Applied Geophysics, 2017,136:484-497.

DOI:10.1016/j.jappgeo.2016.12.001 URL [本文引用: 1]

[26]

, Huang

, Li

, et al.

Attenuation compensation in anisotropic least-squares reverse time migration

[J]. Geophysics, 2017,82(6):S411-S423.

DOI:10.1190/geo2016-0677.1 URL [本文引用: 1]

[27]

Zhu T

, Harris J

, Biondi

Q-compensated reverse-time migration

[J]. Geophysics, 2014,73(9):S77-S87.

[28]

Zhu T

, Harris J

, Biondi

Q-compensated reverse-time migration

[J]. Geophysics, 2014,73(9):S77-S87.

[29]

李振春, 孙小东, 刘洪 .

复杂地表条件下共反射面元(CRS)叠加方法研究

[J]. 地球物理学报, 2006,49(6):1794-1801.

在地表地形复杂的情况下，静校正不易做好，这是制约山地资料处理质量的一个很重要的因素.复杂地表共反射面元（CRS）叠加不需对叠前数据做静校正，而且在得到叠加剖面后可以利用叠加得到的波场参数剖面实现基准面重建.地震数据的试算表明，复杂地表CRS叠加得出的剖面与常规处理剖面相比有着较高的信噪比和同相轴连续性.与水平地表CRS叠加不同的是，在复杂地表CRS叠加的时距公式中，波场三参数耦合，难以通过简化CRS道集的方法将它们全部分离并逐个优化.引入模拟退火算法后，有效地解决了这一组合优化的难题.

[30]

李振春

. 地震叠前成像理论与方法[M]. 东营: 中国石油大学出版社, 2011.

[31]

李振春

地震偏移成像技术研究现状与发展趋势

[J]. 石油地球物理勘探, 2014,49(1):1-21.

地震偏移成像是在一定的数学物理模型（声介质、弹性介质等）基础上，应用相应的地球物理理论，将地面观测到的多次覆盖数据反传，消除地震波的传播效应并得到地下介质模型图像的过程。本文首先从不同的分类角度对地震偏移进行了简单评述，较为详细地阐述了如今广泛应用的几类叠前深度偏移方法，分析了计算效率和成像精度、保幅性、起伏地表以及介质复杂性对偏移成像适用性的制约，指出选取偏移成像方法时需要具体问题具体分析，要明确方法的假设条件和适用范围。最后通过总结近年来SEG/EAGE会议及出版物中有关地震偏移成像的讨论专题，展望了地震偏移未来走向和发展趋势，可以概括为：①叠前深度偏移已经成为地震偏移的研究主流；②逆时偏移逐步由理论研究步入工业化应用；③地震偏移成功地从二维走向三维；④起伏地表偏移受到广泛关注；⑤TI、VTI和TTI介质偏移成为研究热点；⑥反演偏移已经作为偏移领域的新宠登上历史舞台。

[32]

曲英铭, 黄建平, 李振春 , 等.

分层坐标变换法起伏自由地表弹性波叠前逆时偏移

[J]. 地球物理学报, 2015,58(8):2896-2911.

[33]

李振春, 李庆洋, 黄建平 , 等.

裂缝型储层时空双变正演模拟研究

[J]. 地球物理学进展, 2014,29(3):1183-1188.

DOI:10.6038/pg20140324 Magsci [本文引用: 1]

本文将深层裂缝型储层视为各向同性背景下的小尺度非均质体，采用时空双变正演模拟算法模拟其地震响应，并将Lanczos滤波算子推广到该算法中，在保证最优的计算效率及内存的前提下，解决了常规变网格算法长时间采样下的不稳定问题.测试结果验证了本文算法对深层微构造的模拟精度及效率.最后将本算法应用到裂缝储层正演模拟中，定量分析了不同裂缝开度、密度和倾角下的地震振幅响应规律，为裂缝型储层的检测和解释提供了可靠的理论指导.

[34]

曲英铭, 黄建平, 李振春 , 等.

一种基于非规则网格的地震波射线追踪方法

[J]. 石油物探, 2014,53(6):627-632.

[35]

曲英铭, 黄建平, 李振春 , 等.

基于单元交错网格的变坐标系正演模拟方法在声—弹介质中的应用

[J]. 石油物探, 2015,54(5):582-591.

[36]

刘洋, 李向阳, 陈双全 .

高精度双吸收边界条件在地震波场模拟中的应用(英文)

[J]. Applied Geophysics, 2015,12(1):111-119.

DOI:10.1007/s11770-014-0463-z URL [本文引用: 1]

[37]

刘洋

波动方程时空域有限差分数值解及吸收边界条件研究进展

[J]. 石油地球物理勘探, 2014,49(1):35-46.

波动方程数值解是波动方程正演、逆时偏移和全波形反演的核心技术之一。本文对波动方程数值求解的有限差分技术和吸收边界条件进行了分析，重点总结了基于时空域频散关系的有限差分、自适应可变空间算子长度有限差分、优化有限差分及混合吸收边界条件等方法，介绍了这些方法在逆时偏移和波形反演中的应用。

[38]

刘洋, 李承楚, 牟永光 .

任意偶数阶精度有限差分法数值模拟

[J]. 石油地球物理勘探, 1998,33(1):1-10.

本文从Taylor级数展开式出发，推导出了任意阶导数的任意偶数阶精度差分格式。文中给出了计算相应差分系数的公式，并计算出了一阶、二阶、三阶和四阶导数的若干阶精度差分系数。最后讨论了声波方程的任意偶数阶差分数值解的稳定性条件，并进行了数值模拟试验。结果表明差分精度越高，频散越小，数值模拟的效果越好。

[39]

赵岩, 刘洋, 任志明 .

基于优化时空域高阶有限差分方法的粘滞声波叠前逆时偏移(英文)

[J]. Applied Geophysics, 2014,11(1):50-62.

DOI:10.1007/s11770-014-0414-8 URL [本文引用: 1]

[40]

陈可洋

地震波逆时偏移方法研究综述

[J]. 勘探地球物理进展, 2010,33(3):153-159.

[41]

陈可洋

基于高阶有限差分的波动方程叠前逆时偏移方法

[J]. 石油物探, 2009,48(5):475-478.

[42]

邓文志, 李振春, 王延光 , 等.

基于稳定逆时传播算子的黏声介质最小二乘逆时偏移

[J]. 物探与化探, 2015,39(4):791-796.

DOI:10.11720/wtyht.2015.4.22 Magsci [本文引用: 1]

基于GSLS模型黏声介质二阶拟微分方程,采用伪谱法进行数值模拟。针对黏声介质逆时传播过程中产生的高频不稳定问题,提出加入规则化算子对其进行消除的方法,构建了稳定的逆时传播算子。在最小二乘反演的基础上,将黏声介质逆时偏移与最小二乘思路相结合,发展了带有振幅补偿的黏声介质最小二乘逆时偏移(LSRTM)。Marmousi模型结果表明:相对于常规最小二乘逆时偏移,黏声介质最小二乘逆时偏移校正了地层的黏滞性,得到了更加精确可靠的保幅成像剖面。

[43]

刘玉金, 李振春 .

扩展成像条件下的最小二乘逆时偏移

[J]. 地球物理学报, 2015,58(10):3771-3782.

[44]

Romero L

, Ghigliaz D

, Ober C

, et al.

Phase encoding of shot records in prestack migration

[J].Geophysics, 65(2):426-436.

[45]

Schuster G

, Wang

, Huang

, et al.

Theory of multisource crosstalk reduction by phase-encoded statics

[J]. Geophys. J. Int., 2011,184:1289-1303.

DOI:10.1111/gji.2011.184.issue-3 URL [本文引用: 1]

[46]

Huang Y

, Schuster G

Multisource least-squares migration of marine streamer and land data with frequency-division encoding

[J]. Geophysical Prospecting, 2012,60:663-680.

DOI:10.1111/gpr.2012.60.issue-4 URL [本文引用: 1]

[47]

陈可洋, 吴清岭, 李来林 , 等.

松辽盆地三维地震资料连片处理关键技术及其应用效果分析

[J]. 岩性油气藏, 2012,24(2):87-91.