石油物探  2022, Vol. 61 Issue (2): 293-309  DOI: 10.3969/j.issn.1000-1441.2022.02.011
0
文章快速检索     高级检索

引用本文 

陈生昌, 鲁方正, 刘亚楠, 等. 基于地震数据线性正演表达的波形成像(二): 地震数据的波形成像[J]. 石油物探, 2022, 61(2): 293-309. DOI: 10.3969/j.issn.1000-1441.2022.02.011.
CHEN Shengchang, LU Fangzheng, LIU Yanan, et al. Waveform imaging based on linear forward representations of seismic data—Part 2: Waveform imaging of seismic data[J]. Geophysical Prospecting for Petroleum, 2022, 61(2): 293-309. DOI: 10.3969/j.issn.1000-1441.2022.02.011.

基金项目

国家自然科学基金项目(41874133, U19B2008)资助

第一作者简介

陈生昌(1965—), 男, 教授, 博士生导师, 主要从事勘探地球物理和计算地球物理的研究工作。Email: chenshengc@zju.edu.cn

文章历史

收稿日期:2021-01-25
基于地震数据线性正演表达的波形成像(二): 地震数据的波形成像
陈生昌1, 鲁方正1, 刘亚楠1, 周华敏2    
1. 浙江大学地球科学学院, 浙江杭州 310027;
2. 长江水利委员会长江科学院, 湖北武汉 430010
摘要:将第一部分论述中提出的地震数据线性正演表达公式作为地震数据线性反演的正演方程, 利用线性反演理论提出包括散射数据波形成像和反射数据波形成像的地震数据波形成像方法理论, 给出了标量波散射数据、声波反射数据和弹性波反射数据波形成像的具体计算公式。散射数据波形成像是对散射体物性参数相对扰动的线性反演, 反射数据波形成像是对反射体波阻抗相对扰动或反射体边界局部反射系数的线性反演。在波形成像中, 如果将波场传播算子的伴随算子作为波场传播算子的逆, 则波形成像可转化为波形偏移; 如果将波场传播算子的最小二乘逆作为波场传播算子的逆, 则波形成像可以转化为最小二乘波形偏移。散射数据波形偏移可以实现对散射体空间位置的成像, 反射数据波形偏移不仅可以实现对反射体构造的准确成像, 还可以实现对反射体波阻抗相对扰动的成像。受角度域反偏移计算复杂度的限制, 反射数据最小二乘波形偏移仅可实现反射体近法向波阻抗相对扰动的最小二乘反演和反射体边界局部近法向反射系数的最小二乘反演。相较于逆时偏移, 波形偏移结果不仅在理论上具有更加准确的相位和更高的分辨率, 而且还不增加偏移计算量。地震数据波形成像不仅弥补了当前构造成像的不足, 还可以用于地下的岩性成像。合成地震数据的数值试验结果验证了方法的有效性。
关键词散射数据    反射数据    线性正演表达    线性反演    波形偏移    最小二乘波形偏移    
Waveform imaging based on linear forward representations of seismic data—Part 2: Waveform imaging of seismic data
CHEN Shengchang1, LU Fangzheng1, LIU Yanan1, ZHOU Huamin2    
1. School of Earth Sciences, Zhejiang University, Hangzhou 310027, China;
2. Changjiang River Scientific Research Institute, Changjiang Water Resources Commission, Wuhan 430010, China
Abstract: This paper presents the second part of the theoretical work on seismic data waveform imaging.The method relies on the linear forward representation formulas proposed in the first part of the work and makes use of forward equations for the linear inversion of seismic data.Both scattering and reflection data waveform imaging are discussed.Specific calculation formulas are presented for waveform imaging of scalar wave scattering data, acoustic wave reflection data, and elastic wave reflection data.Waveform imaging of scattering data is achieved via a linear inversion of the relative perturbation of the scatterer's physical parameters, whereas imaging of the reflected data stems from a linear inversion of the relative perturbation of the wave impedance of the reflector (or the local reflection coefficient of the reflector boundary).During the process of waveform imaging, if the adjoint operator of the wavefield propagator is used as the inverse of the wavefield propagator, the waveform imaging is transformed into waveform migration.In contrast, if the least-squares inverse of the wavefield propagator is used as the inverse of the wavefield propagator, the waveform imaging is transformed into the least-squares waveform migration.Imaging of the scatterers' spatial position can be obtained via waveform migration of scattering data.Finally, accurate imaging of the structure of the reflector and imaging of the relative perturbation of the reflector wave impedance can both be realized by means of waveform migration of reflection data.Computational complexity of the angle-domain de-migration limits the least-squares waveform migration of reflection data; thus, least-squares inversion can only be performed with respect to the relative perturbation of the near-normal wave impedance of the reflector and the local near-normal reflection coefficient of the reflector boundary.Compared with the current reverse time migration method, waveform migration has a higher resolution and returns a theoretically correct phase without increasing the computational burden.The waveform imaging method proposed in this paper can not only address the shortcomings of current structural imaging but can also be used for subsurface lithology imaging, as demonstrated by the results of synthetic seismic data processing.
Keywords: scattered data    reflection data    linear forward representations    linear inversion    waveform migration    least squares waveform migration    

在波形成像论文的第一篇中, 对于给定的地震波运动学光滑模型, 将扰动法应用于波动方程, 得到了描述非均匀体扰动波场的波动方程。依据地震波长与非均匀体大小或非均匀体物理性质空间变化特征尺度之间的相对大小关系, 将非均匀体划分为在空间上具有有限延续度可产生散射波的局部散射体或空间上具有一定延续度可产生反射波的反射体, 并将相应的扰动波场的波动方程转化为散射波动方程和反射波动方程。对于散射波, 利用小扰动假定和Born近似, 得到基于散射体物性参数相对扰动的散射数据线性正演表达式; 对于反射波, 利用波阻抗相对扰动的小值假定、一次波近似和高频近似(地震波长相对于非均匀体大小或其物理性质空间变化特征尺度为较小量), 得到基于反射体波阻抗相对扰动的反射数据线性正演表达式。为了描述反射体边界对反射波的作用和反射体边界上物理量在反射数据偏移成像中的广泛应用, 在高频近似条件(光滑模型的空间变化相对于地震波长尺度可视为缓慢变化甚至不变化)下定义反射体波阻抗相对扰动沿入射波传播方向的方向导数为反射体边界局部反射系数, 推导出了基于反射体局部反射系数的反射数据线性正演表达式, 并且给出了标量波、声波和弹性波方程下散射数据和反射数据的线性正演表达式。

本文是波形成像论文的第二篇, 将第一篇推导出的地震一次散射数据和一次反射数据线性正演表达式作为地震数据线性反演的正演方程, 利用线性反演方法提出包括散射数据波形成像和反射数据波形成像的地震数据波形成像方法理论。首先根据第一部分论述中提出的地震数据线性正演表达式, 提出用于散射数据波形成像和反射数据波形成像的地震波传播通用方程、散射通用方程和反射通用方程, 再利用线性反演理论建立地震数据波形成像的基本方法理论框架, 最后在波形成像的线性反演中利用地震波传播算子的伴随算子和最小二乘逆算子分别提出波形成像的波形偏移方法和最小二乘波形偏移方法。对于地震反射数据在两种不同线性正演表达下的波形成像, 首先给出面向地层波阻抗相对扰动的岩性成像和反射面局部反射系数的构造成像, 然后针对具体的标量波散射数据、声波反射数据和弹性波反射数据分别给出它们的波形偏移和最小二乘波形偏移具体计算公式, 最后进行合成数据的波形成像数值试验, 试验结果表明波形偏移结果相对于逆时偏移结果具有更加准确的相位和更高的分辨率。

1 地震数据波形成像的方法理论框架

地震数据的正演表达式是地震数据反演成像的基础, 不同的地震数据正演表达式会导致不同的反演目标和反演方法。当前地震数据波动方程偏移方法缺少与之对应的严谨的地震数据正演方程, 使得偏移方法中描述波场传播的方程与描述波场成像的成像公式不匹配, 皆因波场传播方程和成像公式是由不同条件下的波动方程导出的。偏移成像公式应该和产生散射波或反射波的虚源有关, 即应该和散射波或反射波与入射波之间所满足的数学物理方程有关。本文将描述散射波或反射波与入射波之间关系的数学物理方程称之为散射(波)方程或反射(波)方程。在第一部分论述中关于地震数据线性正演表达的基础上, 首先建立与地震数据波形成像有关的地震波传播、散射、反射的通用方程, 然后在此基础上利用反演理论构建地震数据波形成像的方法理论框架。

1.1 用于散射数据波形成像的地震波传播与散射通用方程

当前的地震数据波动方程偏移成像不区分散射和反射, 将无穷大平边界平面波的镜面反射方程(第一部分论述中的(6)式)作为描述散射波与入射波之间关系的散射方程, 用于散射数据偏移成像。第一部分论述中将相对于地震波长尺度可近似为点源的虚源产生的地震波定义为散射波, 因此我们认为在散射波偏移成像中采用第一部分论述中的(6)式描述散射不合适, 它不能准确地描述局部的散射现象和散射波与入射波之间的相位变化。

根据第一部分论述中在光滑模型下建立的散射数据线性正演表达公式, 我们建立用于描述地震波传播与一次散射的通用波动方程, 即:

$ L\left(m_{\mathrm{b}}\right) u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=f(\boldsymbol{x}, t) $ (1)
$ L\left(m_{\mathrm{b}}\right) u_{\mathrm{s}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=-\delta L_{\mathrm{s}}\left(m_{\mathrm{b}}, \delta m\right) u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) $ (2)

式中: mb为已知的光滑模型; L(mb)为光滑模型下的通用波动算子; ui(x, t; xs)为光滑模型下的背景波场或入射波场; f(x, t)为震源函数; xs为震源位置; us(x, t; xs)为一次散射波场; δLs(mb, δm)为散射算子, 它的具体表达式与δm/mb和波动算子L(mb)的具体形式有关; δm为模型扰动, 其空间尺寸大小或其物理性质空间变化特征尺度相对于地震波长为较小量, 满足散射体条件。波动方程(1)描述了光滑模型下的入射波场及其传播, 它所对应的Green函数也是光滑模型下散射波的波场传播算子。波动方程(2)不仅描述了散射波场与入射波场之间应满足的数学物理关系, 也描述了散射波的产生(即产生散射波的虚源)与散射波的传播。我们将(1)式(或其对应的齐次方程)和(2)式统称为用于散射数据波形成像的地震波传播与散射通用方程, 它们是由相同条件下的波动方程推导而来的, 因此二者相互匹配, 由此可以进一步推导出一次散射波数据的线性正演表达, 该方程也是利用一次散射波对模型相对扰动δm/mb或模型扰动δm进行反演成像的基础方程。

1.2 用于反射数据波形成像的地震波传播与反射通用方程

在当前的地震数据波动方程逆时偏移方法中, 将无穷大平边界平面波的镜面反射方程(第一部分论述中的(6)式)作为一般情况下反射波的反射方程, 并用于反射数据偏移成像。这样的镜面反射方程不能准确地描述非平边界、非平面波情况下的反射以及反射波与入射波之间的相位差异。根据第一部分论述中在光滑模型下利用高频近似和局部平面波建立的反射数据线性正演表达公式, 可以得到地震波传播与反射的通用方程。对于基于反射体波阻抗相对扰动产生的一次反射, 可以得到:

$ L\left(m_{\mathrm{b}}\right) u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=f(\boldsymbol{x}, t) $ (3)
$ L\left(m_{\mathrm{b}}\right) u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=-\delta L_{\mathrm{r}}\left(m_{\mathrm{b}}, I_{\mathrm{r}}, \sigma\right) u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) $ (4)

式中: ur(x, t; xs)表示一次反射波场; δLr(mb, Ir, σ)表示随角度σ变化的反射算子, 它的表达式与波阻抗相对扰动Ir和波动算子L(mb)的关系式有关; σ表示入射方向与反射方向之间的开角。(3)式描述了光滑模型下的入射波场及其传播规律, 所对应的Green函数也是光滑模型下反射波的波场传播算子。(4)式不仅反映了反射波场与入射波场之间应满足的数学物理关系, 也描述了基于相对扰动Ir的反射波产生(即产生反射波的虚源)与传播规律。我们将(3)式(或其对应的齐次方程)和(4)式称为基于波阻抗相对扰动Ir的反射数据波形成像地震波传播与反射通用方程, 它们均由相同条件下的波动方程推导而来, 因此二者是相互匹配的, 由此可以推导出基于波阻抗相对扰动Ir的一次反射波数据的线性正演表达公式, 该方程也是利用一次反射波数据对波阻抗相对扰动Ir进行反演成像的基础方程。

考虑基于反射体边界局部反射系数产生的反射, 推导得到如下的地震波传播与反射的通用方程:

$ L\left(m_{\mathrm{b}}\right) u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=f(\boldsymbol{x}, t) $ (5)
$ L\left(m_{\mathrm{b}}\right) u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=-\delta L_{\mathrm{r}}\left(m_{\mathrm{b}}, r, \sigma\right) u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) $ (6)

式中: r表示反射体边界局部反射系数; δLr(mb, r, σ)表示基于r随角度σ变化的反射算子, 它的具体表达式与r和波动算子L(mb)的具体形式有关。(6)式不仅反映反射波场与入射波场之间应满足的数学物理关系, 也描述了基于局部反射系数的反射波的产生与传播规律。我们将(5)式(或其对应的齐次方程)和(6)式称为利用基于局部反射系数的反射数据波形成像的地震波传播与反射通用方程, 它们也是由相同条件下的波动方程推导而来, 因此二者是相互匹配的, 它们的结合可作为基于反射体边界局部反射系数产生的一次反射波数据的线性正演表达, 由此可以推导出利用一次反射波对反射体边界局部反射系数r进行反演成像的基础方程。

由上述分析可知, 对于反射波和散射波而言, 它们具有相同的入射方程和传播方程(或传播算子), 却具有不同的散射方程和反射方程。

1.3 地震数据波形成像的方法理论框架

根据上文中的地震波传播通用方程、散射通用方程和反射通用方程, 再利用(1)式所对应的Green函数, 可得到下述通用的散射数据或反射数据线性积分表达式:

$ u\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\int_{\Omega} \mathrm{d} \boldsymbol{x} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) *{}_{\mathrm{t}} s\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) $ (7)

式中: u(xg, t; xs)代表散射数据或反射数据; xg代表检波点坐标; Ω表示散射波或反射波虚源的空间分布范围; g(xg, x, t)为Green函数, 代表从Ω内xxg的波场传播算子; *t代表时间褶积; s(x, t; xs)为(2)式、(4)式或(6)式右端的散射波或反射波虚源项, 该虚源项包含了光滑模型下的入射波场ui(x, t; xs)、产生散射或反射的模型参数和散射算子或反射算子。

利用(7)式左端项的地震数据u(xg, t; xs), 估计右端积分中的虚源项s(x, t; xs), 即为求解第一类线性积分方程, 也称为反源问题(inverse source problem)[1]。将(7)式代表的积分方程离散为矩阵方程形式, 如下:

$ \mathit{\boldsymbol{d}} = \mathit{\boldsymbol{Lm}} $ (8)

式中: d代表共炮点道集地震数据(即一次散射数据或一次反射数据)矩阵; L代表由波场传播算子组成的波场传播矩阵; m表示待求的地下虚源矩阵。

对于矩阵方程(8), 有下述的最小二乘解:

$ \boldsymbol{m}=\left(\boldsymbol{L}^{*} \boldsymbol{L}\right)^{-1} \boldsymbol{L}^{*} \boldsymbol{d}=\boldsymbol{L}^{-\mathrm{g}} \boldsymbol{d} $ (9)

式中: L-g代表波场传播矩阵L的最小二乘逆矩阵; L*代表矩阵L的伴随矩阵。如果(L*L)为酉矩阵或忽略(L*L)-1的作用, 即用L*作为L的逆矩阵, 则可得到矩阵方程(8)的近似解ma, 即:

$ {\mathit{\boldsymbol{m}}_{\rm{a}}} = {\mathit{\boldsymbol{L}}^*}\mathit{\boldsymbol{d}} $ (10)

对于大规模反演成像而言, 利用最小二乘方法求解(9)式存在逆矩阵计算困难和不稳定性问题, 但得到的解保真性和分辨率高。(10)式的求解方法称为反投影法或反传播法, 具有计算量小、稳定的特点, 所得到的解保真性和分辨率不及最小二乘法所得结果。对于求解(9)式中存在的逆矩阵计算困难和解不稳定问题, 均可由迭代方法解决[2]

求解(9)式或(10)式均可以得到产生散射波或反射波的虚源s(x, t; xs)的解估计。为了得到产生散射波的散射体物性参数扰动或产生反射波的反射体波阻抗相对扰动或反射体边界局部反射系数的成像结果, 我们还需要利用波动方程(1)计算虚源s(x, t; xs)中的入射波场ui(x, t; xs), 然后再从s(x, t; xs)的解估计中消除ui(x, t; xs), 同时还要考虑散射算子或反射算子的具体形式, 并将其从s(x, t; xs)的解估计中消除。这些操作类似于常规逆时偏移中的波场成像运算, 但在常规逆时偏移方法中不区分散射算子和反射算子之间的差异, 并利用无穷大平界面平面波的镜面反射算子R(σ)I近似本文中的散射算子和反射算子, R(σ)为随角度σ变化的反射系数[3], I为单位算子。因此, 偏移成像公式应该与散射或反射虚源中的散射算子或反射算子(或散射或反射虚源的表达式)有关。

根据上述内容可得到如下的地震数据波形成像方法理论框架:

1) 利用观测到的地震数据(散射数据或反射数据)求解反源问题, 得到对地下散射波或反射波虚源的解估计;

2) 利用光滑模型下的波动方程(1)求解地下入射波场;

3) 利用地下入射波场、散射算子或反射算子的具体形式, 根据散射虚源或反射虚源的解估计得到散射体物性参数扰动或反射体波阻抗相对扰动或反射体边界局部反射系数的成像。

在波形成像方法理论框架的反源问题求解中, 利用(10)式得到虚源的解估计, 对应的地震数据波形成像方法称为地震数据波形偏移方法。在波形成像方法理论框架的反源问题求解中, 利用(9)式得到虚源的解估计, 则对应的地震数据波形成像方法称为地震数据最小二乘波形偏移方法。

利用上述最小二乘波形偏移方法, 在理论上可以得到散射体物性参数扰动、反射体波阻抗相对扰动、反射体边界局部反射系数的高分辨率和高保真成像结果。但由于实际地震数据的频带局限性(特别是缺少低频成分)和观测范围的有限性(缺失大偏移距观测), 导致得到的最小二乘波形偏移结果通常缺少低波数成分, 仅仅是一种对散射体物性参数扰动或反射体波阻抗相对扰动的视反演, 因此需要在最小二乘波形偏移中加入有关散射体物性参数扰动或反射体波阻抗相对扰动低波数成分的先验信息, 才能实现对散射体物性参数扰动值或反射体波阻抗相对扰动值的高分辨率和高保真反演。

2 标量波散射数据的波形成像

根据上节给出的地震数据波形成像方法理论框架, 对于不同波动方程所描述的地震数据线性正演表达式, 给出不同的地震波传播方程和散射或反射方程, 进而可以得到不同的地震数据波形成像公式。在给定的光滑速度模型下, 对于标量波动方程所描述的散射线性数据正演表达式, 存在下述具体的地震波传播方程和散射方程, 即:

$ \left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla}^{2}\right] u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=f(t) \delta\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{s}}\right) $ (11)
$ \left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla}^{2}\right] u_{\mathrm{s}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\frac{a_{\mathrm{v}}(\boldsymbol{x})}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} $ (12)

式中: $ \partial^{2}$为Laplace算子, 有$ \partial^{2}$=2/x2+2/y2+2/z2; vb(x)为光滑速度模型; av(x)表示速度相对扰动, 有av(x)=1-vb2(x)/v2(x)≈2δv(x)/vb(x), 其中, δv(x)为满足散射体条件的速度扰动。方程(12)的右端项为产生散射波的虚源, 包含速度相对扰动、散射算子和入射波场, 其中散射算子由光滑速度模型和时间二阶导数组成, 即[1/vb2(x)](2/t2)。相应的散射数据积分表达式为:

$ \begin{aligned} u_{\mathrm{s}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=& \int_{\Omega} \mathrm{d} \boldsymbol{x} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) *{ }_{\mathrm{t}} \frac{a_{\mathrm{v}}(\boldsymbol{x})}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \cdot \\ & \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{aligned} $ (13)

式中: g(xg, x, t)表示光滑速度模型下方程(11)的Green函数。

将上节的地震数据波形成像方法理论框架及其波形偏移的伴随求解方法以及(10)式应用于求解线性积分方程(13), 再经过一些简单的推导可得到下述求取速度相对扰动av(x)近似估计的标量波散射数据波形偏移计算公式。

1) 地下入射波场的计算:

$ \left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla}^{2}\right] u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=f(t) \delta\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{s}}\right) $ (14)

2) 基于边值条件和波场逆时传播的地下散射虚源近似反演:

$ \left\{\begin{array}{l} {\left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla}^{2}\right] u_{\mathrm{s}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=0} \\ u_{\mathrm{s}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=u_{\mathrm{s}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \delta\left(\boldsymbol{x}_{\mathrm{g}}-\boldsymbol{x}\right) \end{array}\right. $ (15)

3) 对av(x)的近似估计:

$ a_{\mathrm{v}}(\boldsymbol{x})=v_{\mathrm{b}}^{2}(\boldsymbol{x}) \int_{0}^{T} \mathrm{~d} t u_{\mathrm{s}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} $ (16)

上述公式中: us(xg, t; xs)表示共炮点道集散射数据; T表示地震数据的最大记录时间。由于我们将产生散射波的散射体相对于地震波长近似为一个点源体, 因此利用(16)式得到的av(x)的近似估计可以实现对散射体空间位置的成像[4]

(16) 式相较于常规逆时偏移成像公式的主要区别是多了对震源波场的时间二阶导数, 因此(16)式得到的成像结果相对于常规逆时偏移结果存在180°的相位差异(即反极性), 而且分辨率也会更高。

为了得到具有高保真性和高分辨率的av(x)反演结果, 根据上节的地震数据波形偏移方法理论框架及最小二乘波形偏移的最小二乘公式, 再结合上述的散射数据波形偏移计算公式, 可通过迭代方式求解av(x)的散射数据最小二乘波形偏移计算公式。

① 初始化。

k=0, avk(x)=av0(x)=0。

入射波场计算: 利用公式(14)计算入射波场。

② 迭代过程

计算数据残差:

$ \begin{gathered} \delta d^{k}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=u_{\mathrm{s}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)-\int_{\Omega} \mathrm{d} \boldsymbol{x} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, \right.\\ t) * {}_{t}\frac{a_{\mathrm{v}}^{k}(\boldsymbol{x})}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{gathered} $ (17)

数据残差的逆时外推:

$ \left\{\begin{array}{l} {\left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla}^{2}\right] \delta u^{k}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=0} \\ \delta u^{k}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\delta d^{k}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \delta\left(\boldsymbol{x}_{\mathrm{g}}-\boldsymbol{x}\right) \end{array}\right. $ (18)

avk(x)的修正, 为:

$ \begin{gathered} a_{\mathrm{v}}^{k+1}(\boldsymbol{x})=a_{\mathrm{v}}^{k}(\boldsymbol{x})+\alpha_{k} v_{\mathrm{b}}^{2}(\boldsymbol{x}) \int_{0}^{T} \mathrm{~d} t \delta u^{k}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)\cdot \\ \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{gathered} $ (19)

k=k+1, 迭代结束。

修正公式(19)中的αk为第k次迭代步长因子(下同), 可采用简单的线性搜索法获得[5], 也可通过震源方向和检波点方向的地震波双向照明强度近似求解[6]得到。

3 声波反射数据的波形成像

根据前文所述的反射数据波形成像的地震波传播与反射通用方程, 针对速度、密度变化的声波动方程所描述的地震反射数据, 在给定光滑速度vb(x)、密度ρb(x)的模型下, 对于基于波阻抗相对扰动产生的反射, 有下述具体的地震波传播方程和反射方程:

$ \begin{gathered} \left\{\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla} \cdot\left[\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \boldsymbol{\nabla}\right]\right\} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \\ =f(t) \delta\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{s}}\right) \end{gathered} $ (20)
$ \begin{gathered} \left\{\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla} \cdot\left[\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \boldsymbol{\nabla}\right]\right\} u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)= \\ \frac{I_{\mathrm{r}}(\boldsymbol{x}, \sigma)}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{gathered} $ (21)

式中: $\boldsymbol{\nabla} $表示梯度算子, 有$\boldsymbol{\nabla} $=(/x)i+(/y)j+(/z)l; $\boldsymbol{\nabla} $·表示散度运算; Ir(x, σ)表示与角度σ有关的声波阻抗相对扰动, 有Ir(x, σ)=av+aρ(1+cosσ), 其中av为满足反射体条件的速度相对扰动, 有av=[2δv(x)]/[vb2(x)]; aρ为满足反射体条件的密度相对扰动, 有aρ=δρ(x)/ρb(x), 其中δρ(x)表示密度扰动, 进一步有: Ir(x, σ)=2(δv/vb)sin2(σ/2)+2(δI/Ib)cos2(σ/2), 其中δI=ρδv+vδρ为声阻抗扰动, Ib=ρbvb为背景声阻抗。方程(20)的右端项就是产生反射波的虚源, 包含声波阻抗相对扰动、反射算子和入射波场, 其中反射算子由光滑速度、密度模型和时间二阶导数组成, 即1/ρb(xvb2(x)(2/t2)。相应的反射数据积分表达为:

$ \begin{gathered} u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\int_{\Omega} \mathrm{d} \boldsymbol{x} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) * {}_\mathrm{t} \frac{I_{\mathrm{r}}(\boldsymbol{x}, \sigma)}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})}\cdot \\ \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{gathered} $ (22)

式中: g(xg, x, t)表示光滑速度、密度模型下声波方程(20)的Green函数。

如果不考虑密度的变化, 即δρ(x)=0, 则与角度有关的声波阻抗相对扰动简化为与角度无关的速度相对扰动, 即Ir(x, σ)=av(x), 方程(20)至方程(22)简化为下述描述基于速度相对扰动产生的反射所对应的方程:

$ \left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla}^{2}\right] u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=f(t) \delta\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{s}}\right) $ (20-1)
$ \begin{gathered} {\left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla}^{2}\right] u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\frac{a_{\mathrm{v}}(\boldsymbol{x})}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \cdot} \\ \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{gathered} $ (21-1)
$ \begin{gathered} u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\int_{\Omega} \mathrm{d} \boldsymbol{x} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) *{}_{\mathrm{t}}\cdot \\ \frac{a_{\mathrm{v}}(\boldsymbol{x})}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{gathered} $ (22-1)

对于基于反射体边界局部反射系数产生的反射, 可以得到具体的地震波传播方程、反射方程和反射数据积分表达式, 如下:

$ \begin{gathered} \left\{\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla} \cdot\left[\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \boldsymbol{\nabla}\right]\right\} \cdot \\ u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=f(t) \delta\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{s}}\right) \end{gathered} $ (23)
$ \begin{gathered} \left\{\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla} \cdot\left[\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \boldsymbol{\nabla}\right]\right\} u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \\ =\frac{-r(\boldsymbol{x}, \sigma)}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}(\boldsymbol{x})} \frac{\partial u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} \end{gathered} $ (24)
$ \begin{gathered} u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=-\int_{\Omega} \mathrm{d} \boldsymbol{x} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) *{}_{\mathrm{t}}\cdot \\ \frac{r(\boldsymbol{x}, \sigma)}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}(\boldsymbol{x})} \frac{\partial u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} \end{gathered} $ (25)

式中: r(x, σ)表示反射体边界的局部反射系数, 是声波阻抗相对扰动Ir(x, σ)沿入射波传播方向的方向导数, r(x, σ)=(iω)/[vb(x)]Ir(x, σ)=(iω)/[vb(x)][av+aρ(1+cosσ)]。方程(24)的右端项是产生反射波的虚源, 包含局部反射系数、反射算子和入射波场, 其中反射算子由光滑速度、密度模型和时间一阶导数组成, 即-1/[ρb(x)vb(x)]/t

对于方程(23)至方程(25), 如果不考虑密度的变化, 则δρ(x)=0, 可简化为如下方程:

$ \left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla}^{2}\right] u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=f(t) \delta\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{s}}\right) $ (23-1)
$ \begin{gathered} {\left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla}^{2}\right] u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\frac{-r(\boldsymbol{x})}{v_{\mathrm{b}}(\boldsymbol{x})} \cdot} \\ \frac{\partial u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} \end{gathered} $ (24-1)
$ \begin{aligned} u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=-& \int_{\Omega} \mathrm{d} \boldsymbol{x} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) * {}_{\mathrm{t}} \frac{r(\boldsymbol{x})}{v_{\mathrm{b}}(\boldsymbol{x})} \cdot \\ & \frac{\partial u_{\mathrm{i}}\left(\boldsymbol{x}, t; \boldsymbol{\boldsymbol { x } _ { \mathrm { s } } )}\right.}{\partial t} \end{aligned} $ (25-1)

式中: r(x)代表反射体边界的局部反射系数, 是速度相对扰动沿入射波传播方向的方向导数, r(x)=iω/vb(x)av(x)。

将上文的地震数据波形成像方法理论框架及波形偏移的伴随求解方法公式(10)分别应用于线性积分方程(22)和方程(25), 再经过一些简单的推导可以得到求取声波阻抗相对扰动Ir(x, σ)和反射体边界局部反射系数r(x, σ)近似估计的声波反射数据波形偏移计算公式, 如下。

1) 地下入射波场的计算:

$ \begin{gathered} \left\{\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla} \cdot\left[\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \boldsymbol{\nabla}\right]\right\} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \\ =f(t) \delta\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{s}}\right) \end{gathered} $ (26)

2) 基于边值条件和波场逆时传播的地下反射波虚源近似反演:

$ \left\{\begin{array}{l} \left\{\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla} \cdot\left[\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \boldsymbol{\nabla}\right]\right\}\cdot \\ \ \ \ \ \ \ \ \ u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=0 \\ u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \delta\left(\boldsymbol{x}_{\mathrm{g}}-\boldsymbol{x}\right) \end{array}\right. $ (27)

3) 对波阻抗相对扰动Ir(x, σ)和局部反射系数r(x, σ)的近似估计:

$ \begin{gathered} I_{\mathrm{r}}(\boldsymbol{x}, \sigma)=2 \frac{\delta v}{v_{\mathrm{b}}} \sin ^{2} \frac{\sigma}{2}+2 \frac{\delta I}{I_{\mathrm{b}}} \cos ^{2} \frac{\sigma}{2}=\rho_{\mathrm{b}}(\boldsymbol{x})\cdot \\ v_{\mathrm{b}}^{2}(\boldsymbol{x}) \int_{0}^{T} \mathrm{~d} t u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{gathered} $ (28)
$ \begin{gathered} r(\boldsymbol{x}, \sigma)=-\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}(\boldsymbol{x}) \int_{0}^{T} \mathrm{~d} t u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \cdot \\ \frac{\partial u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} \end{gathered} $ (29)

式中: ur(xg, t; xs)为共炮点道集反射数据。(28)式和(29)式相较于常规逆时偏移成像公式的主要区别在于多了时间二阶导数和时间一阶导数, 因此它们得到的成像结果相对于常规逆时偏移结果分别存在180°和90°的相位差异, 此外分辨率有所提高。

为了得到σ角度域的Ir(x, σ)和r(x, σ), 在进行(28)式和(29)式右边的运算前需要将波场ur(x, t; xs)和ui(x, t; xs)转换到相应的角度域[7-8]。由(28)式可知, 利用角度域的共成像点道集Ir(x, σ)数据和Ir(x, σ)的表达式可以进行类似于基于Zoeppritz方程的振幅随入射角变化(AVA)速度扰动、密度扰动和波阻抗扰动分析和反演。但是波阻抗相对扰动Ir(x, σ)不同于Zoeppritz方程中随入射角变化的反射系数, 在地震波长相对于反射体的大小或其物理性质空间变化的特征尺度满足高频近似条件下, 其表达式是由波动方程导出的, 而Zoeppritz方程是由波动方程在平边界和平面波条件下导出的; 再者Ir(x, σ)中的角度σ是根据局部切平面和局部平面波的概念定义的, 是x点处的局部角度, 不同于Zoeppritz方程中假定无穷大平边界和平面波条件下定义的角度。由于本文定义的局部反射系数与入射波的局部波数有关(即与频率有关), 因此我们目前还不能推导出r(x, σ)与速度变化、密度变化、波阻抗变化和角度之间的具体函数关系, 因此利用波形偏移得到的角度域共成像点道集r(x, σ)难以进行AVA速度分析、密度分析和反演。

如果在应用(28)式和(29)式进行偏移成像前不进行波场ur(x, t; xs), ui(x, t; xs)的角度分解, 则由(28)式和(29)式得到的成像结果是Ir(x, σ), r(x, σ)在σ角度域的平均值Ir(x), r(x)。由于反射地震勘探中的偏移距一般都比较小, 所以σ角度域的平均值Ir(x), r(x)分别体现了波阻抗的相对扰动δI(x)/Ib(x)以及局部界面近法线方向的局部反射系数的影响。

如果不考虑密度的变化, 即δρ(x)=0, 将上文中的地震数据波形成像方法理论框架及其波形偏移的伴随求解方法公式(10)应用于线性积分方程(22-1)和方程(25-1), 再经过一些简单的推导可分别得到常密度声波反射数据波形偏移计算公式、速度相对扰动av(x)和反射体边界局部反射系数r(x)近似估计公式。

1) 地下入射波场的计算:

$ \left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla}^{2}\right] u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=f(t) \delta\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{s}}\right) $ (26-1)

2) 基于边值条件和波场逆时传播的地下反射波虚源近似反演:

$ \left\{\begin{array}{l} {\left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla}^{2}\right] u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=0} \\ u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \delta\left(\boldsymbol{x}_{\mathrm{g}}-\boldsymbol{x}\right) \end{array}\right. $ (27-1)

3) 对速度相对扰动av(x)和局部反射系数r(x)的近似估计:

$ a_{v}(\boldsymbol{x})=v_{\mathrm{b}}^{2}(\boldsymbol{x}) \int_{0}^{T} \mathrm{~d} t u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} $ (28-1)
$ r(\boldsymbol{x})=-v_{\mathrm{b}}(\boldsymbol{x}) \int_{0}^{T} \mathrm{~d} t u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \frac{\partial u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} $ (29-1)

由于标量波方程(或常密度声波方程)所描述的反射波与入射角无关[9], 所以标量波散射数据的波形成像公式(14)式至(16)式和常密度声波反射数据的波形成像速度扰动公式(26-1)式至(28-1)式在形式上完全一致。

为了得到具有高保真性和高分辨率的声波阻抗相对扰动反演结果, 根据上节的地震数据波形成像方法理论框架及代表最小二乘波形偏移解的(9)式, 再结合上述的声波反射数据波形偏移计算公式, 可得到下述利用迭代方法求解Ir(x, σ)的声波反射数据最小二乘波形偏移计算公式。由于角度域反偏移计算的困难, 本文中的声波阻抗相对扰动最小二乘波形偏移成像只考虑对Ir(x, σ)在σ角度域的平均值Ir(x)进行估计, 具体迭代过程如下。

1) 初始化。

k=0, Trk(x)=Tr0(x)=0。

入射波场计算: 利用(26)式计算入射波场。

2) 迭代过程。

计算数据残差:

$ \begin{array}{l} \delta d^{k}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)-\int_{\Omega} \mathrm{d} \boldsymbol{x} \cdot\\ \ \ \ \ g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) *{}_{\mathrm{t}} \frac{\bar{I}_{\mathrm{r}}^{k}(\boldsymbol{x})}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{array} $ (30)

数据残差的逆时外推:

$ \left\{\begin{array}{l} \left\{\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla} \cdot\left[\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \boldsymbol{\nabla}\right]\right\}\cdot \\ \ \ \ \ \ \ \ \ \delta \boldsymbol{u}_{\mathrm{r}}^{k}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=0 \\ \delta \boldsymbol{u}_{\mathrm{r}}^{k}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\delta d^{k}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \delta\left(\boldsymbol{x}_{\mathrm{g}}-\boldsymbol{x}\right) \end{array}\right. $ (31)

Trk(x)的修正:

$ \begin{gathered} \bar{I}_{\mathrm{r}}^{k+1}(\boldsymbol{x})=\bar{I}_{\mathrm{r}}^{k}(\boldsymbol{x})+\alpha_{k} \rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x}) \int_{0}^{T} \mathrm{~d} t \delta\cdot \\ \boldsymbol{u}_{\mathrm{r}}^{k}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{gathered} $ (32)

k=k+1, 迭代结束。

如果不考虑密度的变化, 即δρ(x)=0, 可以得到如下的速度相对扰动av(x)的反射数据最小二乘波形偏移。

1) 初始化。

k=0, avk(x)=av0(x)=0。

入射波场计算: 利用(26-1)式计算入射波场。

2) 迭代过程。

计算数据残差:

$ \begin{gathered} \delta d^{k}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)-\int_{\Omega} \mathrm{d} \boldsymbol{x} g\cdot \\ \left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) *{ }_{\mathrm{t}} \frac{a_{\mathrm{v}}^{k}(\boldsymbol{x})}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{gathered} $ (33)

数据残差的逆时外推:

$ \left\{\begin{array}{l} {\left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla}^{2}\right] \delta \boldsymbol{u}_{\mathrm{r}}^{k}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=0} \\ \delta \boldsymbol{u}_{\mathrm{r}}^{k}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\delta d^{k}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \delta\left(\boldsymbol{x}_{\mathrm{g}}-\boldsymbol{x}\right) \end{array}\right. $ (34)

avk(x)的修正:

$ \begin{gathered} a_{\mathrm{v}}^{k+1}(\boldsymbol{x})=a_{\mathrm{v}}^{k}(\boldsymbol{x})+\alpha_{k} v_{\mathrm{b}}^{2}(\boldsymbol{x}) \int_{0}^{T} \mathrm{~d} t \delta \boldsymbol{u}_{\mathrm{r}}^{k} \cdot\\ \left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{gathered} $ (35)

k=k+1, 迭代结束。

为了得到具有高保真性和高分辨率的局部反射系数反演结果, 首先根据上节的地震数据波形偏移方法理论框架及代表最小二乘波形偏移的(9)式, 再结合上述利用声波反射数据进行局部反射系数波形偏移计算, 最终可通过迭代方式求解r(x, σ)的声波反射数据最小二乘波形偏移计算公式。同样地, 由于角度域反偏移计算存在困难, 本文的局部反射系数最小二乘波形偏移也只考虑对r(x, σ)在σ角度域的平均值r(x)进行估计, 具体迭代流程如下。

1) 初始化。

k=0, rk(x)=r0(x)=0。

入射波场计算: 利用(26)式计算入射波场。

2) 迭代过程。

计算数据残差:

$ \begin{array}{c} \delta d_{r}{ }^{k}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)+\int_{\Omega} \boldsymbol{x} g \cdot \\ \left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) * \frac{\bar{r}^{k}(\boldsymbol{x})}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}(\boldsymbol{x})} \frac{\partial u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} \end{array} $ (36)

数据残差的逆时外推:

$ \left\{\begin{array}{l} \left\{\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla} \cdot\left[\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \boldsymbol{\nabla}\right]\right\} \cdot \\ \ \ \ \ \ \ \ \ \delta \boldsymbol{u}_{\mathrm{r}}^{k}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=0 \\ \delta \boldsymbol{u}_{\mathrm{r}}^{k}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\delta d_{\mathrm{r}}^{k}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \delta\left(\boldsymbol{x}_{\mathrm{g}}-\boldsymbol{x}\right) \end{array}\right. $ (37)

rk(x)的修正:

$ \begin{gathered} \bar{r}{}^{k+1}(\boldsymbol{x})=\bar{r}{}^{k}(\boldsymbol{x})-\alpha_{k} \rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}(\boldsymbol{x}) \int_{0}^{T} \mathrm{~d} t \delta \boldsymbol{u}_{\mathrm{r}}^{k} \cdot \\ \left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \frac{\partial u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} \end{gathered} $ (38)

k=k+1, 迭代结束。

如果不考虑密度的变化, 即δρ(x)=0, 得到下述局部反射系数下r(x)反射数据最小二乘波形偏移成像结果, r(x)是速度相对扰动沿入射波传播方向的方向导数。

1) 初始化。

k=0, rk(x)=r0(x)=0。

入射波场计算: 利用(26-1)式计算入射波场。

2) 迭代过程。

计算数据残差:

$ \begin{gathered} \delta d_{r}^{k}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)+\int_{\Omega} \mathrm{d} \boldsymbol{x} g\cdot \\ \left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) * \frac{r^{k}(\boldsymbol{x})}{v_{\mathrm{b}}(\boldsymbol{x})} \frac{\partial u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} \end{gathered} $ (39)

数据残差的逆时外推:

$ \left\{\begin{array}{l} \left\{\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\boldsymbol{\nabla}^{2}\right\} \delta \boldsymbol{u}_{\mathrm{r}}^{k}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=0 \\ \delta \boldsymbol{u}_{\mathrm{r}}^{k}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\delta d_{\mathrm{r}}^{k}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \delta\left(\boldsymbol{x}_{\mathrm{g}}-\boldsymbol{x}\right) \end{array}\right. $ (40)

rk(x)的修正:

$ \begin{gathered} r^{k+1}(\boldsymbol{x})=r^{k}(\boldsymbol{x})-\alpha_{k} v_{\mathrm{b}}(\boldsymbol{x}) \int_{0}^{T} \mathrm{~d} t \delta \boldsymbol{u}_{\mathrm{r}}^{k} \cdot \\ \left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \frac{\partial u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} \end{gathered} $ (41)

k=k+1, 迭代结束。

4 弹性波反射数据的波形成像

对于纵、横波速度和密度变化的各向同性弹性波动方程所描述的地震反射数据, 假定地震波在给定的纵、横波速度和密度光滑模型αb(x), βb(x), ρb(x)下不发生散射、反射和波型转换, 利用反射数据波形成像的地震波传播与反射通用方程, 可以得到如下的弹性入射波方程和一次反射波方程:

$ \boldsymbol{L}_{\mathrm{b}} \boldsymbol{u}_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=f\left(t, \boldsymbol{x}_{\mathrm{s}}\right) $ (42)
$ \boldsymbol{L}_{\mathrm{b}} \boldsymbol{u}_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\Delta \boldsymbol{L}_{\mathrm{r}} \boldsymbol{u}_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) $ (43)

式中: Lb为光滑模型下的弹性波传播算子; ui(x, t; xs)为光滑模型下入射波的位移矢量波场, 在笛卡尔坐标系下, 有ui(x, t; xs)=[uix(x, t; xs), uiy(x, t; xs), uiz(x, t; xs)]; f(t, xs)为矢量震源函数; ur(x, t; xs)为弹性一次反射波的位移矢量波场, 在笛卡尔坐标系下, 有ur(x, t; xs)=[urx(x, t; xs), ury(x, t; xs), urz(x, t; xs)]; $\Delta $Lr为弹性波反射算子。由方程(43)可以得到弹性反射波数据线性正演表达公式, 即

$ \boldsymbol{u}_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\int_{\Omega} \mathrm{d} \boldsymbol{x} \boldsymbol{g}\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) *{}_{t} \Delta \boldsymbol{L}_{\mathrm{r}} \boldsymbol{u}_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) $ (44)

式中: g(xg, x, t)为光滑模型下弹性波方程(42)所对应的矢量Green函数。

由于笛卡尔坐标系下的弹性波场中P波和S波是耦合的, 在弹性波反射数据的逆时偏移成像中, 为了避免弹性波场的分量对分量成像产生的交叉噪声以及得到易于物理解释的成像结果, 需要在成像前将弹性波场解耦为P/S波型波场[10-15]。本文将笛卡尔坐标系下的弹性波场转换成以P、SH和SV波极化方向为坐标系下的P/S波型弹性波场, 实现弹性波场的P波和S波解耦。

在各向同性介质中, 弹性波场的传播方向与P波的极化方向平行, 而P、SH和SV波极化方向相互正交。首先利用Ponyting矢量[16]分别获取入射波场ui(x, t; xs)和反射波场ur(x, t; xs)的入射P波极化方向矢量和反射P波极化方向矢量, 然后对这两个极化方向矢量进行单位化处理, 得到入射P波极化方向单位矢量$ \hat{\boldsymbol{e}}{}_{\mathrm{i}}^{\mathrm{P}}$和反射P波极化方向单位矢量$ \hat{\boldsymbol{e}}{}_{\mathrm{r}}^{\mathrm{P}}$。在反射点以$ \hat{\boldsymbol{e}}{}_{\mathrm{r}}^{\mathrm{P}}$$ \hat{\boldsymbol{e}}{}_{\mathrm{i}}^{\mathrm{P}}$定义参考平面, 垂直该参考平面的方向(右手坐标系)为SH波极化方向, 进而得到SH波极化方向单位矢量$ \hat{\boldsymbol{e}}{}^{\mathrm{SH}}$, 即$ \hat{\boldsymbol{e}}{}^{\mathrm{SH}}$=$ \hat{\boldsymbol{e}}{}_{\mathrm{r}}^{\mathrm{P}}$×$ \hat{\boldsymbol{e}}{}_{\mathrm{i}}^{\mathrm{P}}$。由$ \hat{\boldsymbol{e}}{}_{\mathrm{i}}^{\mathrm{P}}$$ \hat{\boldsymbol{e}}{}_{\mathrm{r}}^{\mathrm{P}}$$ \hat{\boldsymbol{e}}{}^{\mathrm{SH}}$可得到入射SV波极化方向单位矢量$\hat{\boldsymbol{e}}{}_{\mathrm{i}}^{\mathrm{SV}}$和反射SV波极化方向单位矢量$\hat{\boldsymbol{e}}{}_{\mathrm{r}}^{\mathrm{SV}}$, 即$\hat{\boldsymbol{e}}{}_{\mathrm{i}}^{\mathrm{SV}}$=$ \hat{\boldsymbol{e}}{}_{\mathrm{i}}^{\mathrm{P}}$×$ \hat{\boldsymbol{e}}{}^{\mathrm{SH}}$, $\hat{\boldsymbol{e}}{}_{\mathrm{r}}^{\mathrm{SV}}$=$ \hat{\boldsymbol{e}}{}_{\mathrm{r}}^{\mathrm{P}}$×$ \hat{\boldsymbol{e}}{}^{\mathrm{SH}}$。利用投影法得到弹性波场在各个极化方向的分量, 对于入射波场ui(x, t; xs), 有${\boldsymbol{u}}{}_{\mathrm{i}}^{\mathrm{P}}$(x, t; xs)=$ \hat{\boldsymbol{e}}{}_{\mathrm{i}}^{\mathrm{P}}$·ui(x, t; xs), uiSH(x, t; xs)=$ \hat{\boldsymbol{e}}{}^{\mathrm{SH}}$·ui(x, t; xs), ${\boldsymbol{u}}{}_{\mathrm{i}}^{\mathrm{SV}}$(x, t; xs)=$\hat{\boldsymbol{e}}{}_{\mathrm{i}}^{\mathrm{SV}}$·ui(x, t; xs); 对于反射波场ur(x, t; xs), 有${\boldsymbol{u}}{}_{\mathrm{r}}^{\mathrm{P}}$(x, t; xs)=$ \hat{\boldsymbol{e}}{}_{\mathrm{r}}^{\mathrm{P}}$·ur(x, t; xs), urSH(x, t; xs)=$ \hat{\boldsymbol{e}}{}^{\mathrm{SH}}$·ur(x, t; xs), ${\boldsymbol{u}}{}_{\mathrm{r}}^{\mathrm{SV}}$(x, t; xs)=$\hat{\boldsymbol{e}}{}_{\mathrm{r}}^{\mathrm{SV}}$·ur(x, t; xs)。经波场解耦后的入射波场和反射波场, 有uips(x, t; xs)=[${\boldsymbol{u}}{}_{\mathrm{i}}^{\mathrm{P}}$(x, t; xs), uiSH(x, t; xs), ${\boldsymbol{u}}{}_{\mathrm{i}}^{\mathrm{SV}}$(x, t; xs)], urps(x, t; xs)=[${\boldsymbol{u}}{}_{\mathrm{r}}^{\mathrm{P}}$(x, t; xs), urSH(x, t; xs), ${\boldsymbol{u}}{}_{\mathrm{r}}^{\mathrm{SV}}$(x, t; xs)]。

将光滑模型下的弹性波传播算子Lb和反射算子$\Delta $Lr转换成与P/S波型对应的算子, 根据波形成像论文的第一篇, 本文(43)式可表示为:

$ \boldsymbol{L}_{\mathrm{b}}^{\mathrm{ps}} \boldsymbol{u}_{\mathrm{r}}^{\mathrm{ps}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\rho_{\mathrm{b}} \boldsymbol{E}_{\mathrm{r}}^{\mathrm{ips}} \frac{\partial^{2} \boldsymbol{u}_{\mathrm{i}}^{\mathrm{ps}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} $ (45)

其中, 弹性阻抗相对扰动矩阵Erips, 可表示为:

$ \boldsymbol{E}_{\mathrm{r}}^{\mathrm{ips}}=\left(\begin{array}{ccc} I_{\mathrm{r}}^{\mathrm{pp}} & 0 & I_{\mathrm{r}}^{\mathrm{psv}} \\ 0 & I_{\mathrm{r}}^{\mathrm{shsh}} & 0 \\ I_{\mathrm{r}}^{\mathrm{svp}} & 0 & I_{\mathrm{r}}^{\mathrm{svsv}} \end{array}\right) $ (46)

Erips中的元素代表不同类型入射波与不同类型反射波或转换波间的弹性波阻抗相对扰动, 具体如下:

$ I_{\mathrm{r}}^{\mathrm{pp}}=a_{\alpha}+a_{\rho}(1+\cos \sigma)-2 \frac{\beta_{\mathrm{b}}^{2}}{\alpha_{\mathrm{b}}^{2}}\left(a_{\rho}+a_{\beta}\right) \sin ^{2} \sigma $ (47)
$ I_{\mathrm{r}}^{\mathrm{shsh}}=a_{\beta} \cos \sigma+a_{\rho}(1+\cos \sigma) $ (48)
$ I_{\mathrm{r}}^{\mathrm{svsv}}=a_{\beta} \cos 2 \sigma+a_{\rho}(\cos \sigma+\cos 2 \sigma) $ (49)
$ I_{\mathrm{r}}^{\mathrm{svp}}=\sin \sigma\left[a_{\rho}\left(\frac{\alpha_{\mathrm{b}}}{\beta_{\mathrm{b}}}+2 \cos \sigma\right)+2 a_{\beta} \cos \sigma\right] $ (50)
$ \begin{aligned} I_{\mathrm{r}}^{\mathrm{psv}} &=-\frac{\beta_{\mathrm{b}}^{2}}{\alpha_{\mathrm{b}}^{2}} \sin \sigma\left[a_{\rho}\left(\frac{\alpha_{\mathrm{b}}}{\beta_{\mathrm{b}}}+2 \cos \sigma\right)+2 a_{\beta} \cos \sigma\right] \\ &=-\frac{\beta_{\mathrm{b}}^{2}}{\alpha_{\mathrm{b}}^{2}} I_{\mathrm{r}}^{\mathrm{svp}} \end{aligned} $ (51)

上述弹性波阻抗相对扰动Irmn(m, n=p, sh, sv)不仅与纵、横波速度、密度扰动有关, 还与反射开角σ有关。由于σ的正负值与炮、检的相对方向有关, 在Irsvp, Irpsv的表达式中包含有关于角度σ的奇函数sinσ, 因此Irsvp, Irpsvσ=0的两侧会出现极性反转。

根据地震数据波形成像论文的第一篇, (45)式可表示为基于弹性局部反射系数矩阵Errps的P/S波型反射波方程, 即:

$ \boldsymbol{L}_{\mathrm{b}}^{\mathrm{ps}} \boldsymbol{u}_{\mathrm{r}}^{\mathrm{ps}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=-\rho_{\mathrm{b}} \boldsymbol{E}_{\mathrm{r}}^{\mathrm{rps}} \frac{\partial \boldsymbol{u}_{\mathrm{i}}^{\mathrm{ps}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} $ (52)

其中, Errps为反射体边界弹性局部反射系数矩阵, 其表达式为:

$ \boldsymbol{E}_{\mathrm{r}}^{\mathrm{rps}}=\left(\begin{array}{ccc} \alpha_{\mathrm{b}} r^{\mathrm{pp}} & 0 & \beta_{\mathrm{b}} r^{\mathrm{psv}} \\ 0 & \beta_{\mathrm{b}} r^{\mathrm{shsh}} & 0 \\ \alpha_{\mathrm{b}} r^{\mathrm{svp}} & 0 & \beta_{\mathrm{b}} r^{\mathrm{svsv}} \end{array}\right) $ (53)

P波入射、P波反射的波场局部反射系数为:

$ \begin{gathered} r^{\mathrm{pp}}=\mathrm{i} k_{\mathrm{i}}^{\mathrm{p}} I_{\mathrm{r}}^{\mathrm{pp}}=\frac{\mathrm{i} \omega}{\alpha_{\mathrm{b}}}\left[a_{\alpha}+a_{\rho}(1+\cos \sigma)-2 \frac{\beta_{\mathrm{b}}^{2}}{\alpha_{\mathrm{b}}^{2}}\left(a_{\rho}+\right.\right. \\ \left.\left.a_{\beta}\right) \sin ^{2} \sigma\right] \end{gathered} $ (54)

SH波入射、SH波反射的波场局部反射系数为:

$ r^{\text {shsh }}=\mathrm{i} k_{\mathrm{i}}^{\mathrm{s}} I_{\mathrm{r}}^{\text {shsh }}=\frac{\mathrm{i} \omega}{\beta_{\mathrm{b}}}\left[a_{\beta} \cos \sigma+a_{\rho}(1+\cos \sigma)\right] $ (55)

SV波入射、SV波反射的波场局部反射系数为:

$ r^{\mathrm{svsv}}=\mathrm{i} k_{\mathrm{i}}^{\mathrm{s}} I_{\mathrm{r}}^{\mathrm{svsv}}=\frac{\mathrm{i} \omega}{\beta_{\mathrm{b}}}\left[a_{\beta} \cos 2 \sigma+a_{\rho}(\cos \sigma+\cos 2 \sigma)\right] $ (56)

P波入射、SV波反射的波场局部反射系数为:

$ r^{\mathrm{svp}}=\mathrm{i} k_{\mathrm{i}}^{\mathrm{p}} I_{\mathrm{r}}^{\mathrm{svp}}=\frac{\mathrm{i} \omega}{\alpha_{\mathrm{b}}} \sin \sigma\left[a_{\rho}\left(\frac{\alpha_{\mathrm{b}}}{\beta_{\mathrm{b}}}+2 \cos \sigma\right)+2 a_{\beta} \cos \sigma\right] $ (57)

SV波入射、P波反射的波场局部反射系数为:

$ \begin{gathered} r^{\mathrm{psv}}=\mathrm{i} k_{\mathrm{i}}^{\mathrm{s}} I_{\mathrm{r}}^{\mathrm{psv}}=-\frac{\mathrm{i} \omega}{\beta_{\mathrm{b}}} \frac{\beta_{\mathrm{b}}^{2}}{\alpha_{\mathrm{b}}^{2}} \sin \sigma\left[a_{\rho}\left(\frac{\alpha_{\mathrm{b}}}{\beta_{\mathrm{b}}}+2 \cos \sigma\right)+\right. \\ \left.2 a_{\beta} \cos \sigma\right]=-\frac{\beta_{\mathrm{b}}}{\alpha_{\mathrm{b}}} r^{\mathrm{psv}} \end{gathered} $ (58)

由上述的rmn(m, n=p, sh, sv)表达式可知, rsvp, rpsv是关于σ的奇函数, 所以rsvp, rpsvσ=0的两侧也会出现极性反转。

利用上述的P波、SH波和SV波极化方向单位矢量和P波、SH波和SV波标量波场(本文的P波、SH波和SV波分量), 还可以得到P、SH和SV波的矢量波场[12-13, 17]。对入射波场ui(x, t; xs)而言, P波的矢量波场为:

$ \begin{aligned} \boldsymbol{u}_{\mathrm{i}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=& \hat{\boldsymbol{e}}_{\mathrm{\textrm {P }}}^{\mathrm{P}} \boldsymbol{u}_{\mathrm{i}}^{\mathrm{P}}\left(\boldsymbol{x}, t \boldsymbol{;} \boldsymbol{x}_{\mathrm{s}}\right)=\hat{\boldsymbol{e}}_{\mathrm{i}}^{\mathrm{P}}\left[\hat{\boldsymbol{e}}_{\mathrm{i}}^{\mathrm{P}} \cdot\right.\\ &\left.\boldsymbol{u}_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)\right] \end{aligned} $ (59)

SH波的矢量波场为:

$ \begin{aligned} \boldsymbol{u}_{\mathrm{i}}^{\mathrm{SH}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=& \hat{\boldsymbol{e}}^{\mathrm{SH}} \boldsymbol{u}_{\mathrm{i}}^{\mathrm{SH}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\hat{\boldsymbol{e}}^{\mathrm{SH}}\left[\hat{\boldsymbol{e}}^{\mathrm{SH}} \cdot\right.\\ &\left.\boldsymbol{u}_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)\right] \end{aligned} $ (60)

SV波的矢量波场为:

$ \begin{aligned} \boldsymbol{u}_{\mathrm{i}}^{\mathrm{SV}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=& \hat{\boldsymbol{e}}_{\mathrm{i}}^{\mathrm{SV}} \boldsymbol{u}_{\mathrm{i}}^{\mathrm{SV}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\hat{\boldsymbol{e}}_{\mathrm{i}}^{\mathrm{SV}}\left[\hat{\boldsymbol{e}}_{\mathrm{i}}^{\mathrm{SV}} \cdot\right.\\ &\left.\boldsymbol{u}_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)\right] \end{aligned} $ (61)

对于反射波场ur(x, t; xs), P波的矢量波场为:

$ \begin{gathered} \boldsymbol{u}_{\mathrm{r}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\hat{\boldsymbol{e}}_{\mathrm{r}}^{\mathrm{P}} \boldsymbol{u}_{\mathrm{r}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\hat{\boldsymbol{e}}_{\mathrm{r}}^{\mathrm{P}}\left[\hat{\boldsymbol{e}}_{\mathrm{r}}^{\mathrm{P}} \cdot\right. \\ \left.\boldsymbol{u}_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)\right] \end{gathered} $ (62)

SH波的矢量波场为:

$ \begin{aligned} \boldsymbol{u}_{\mathrm{r}}^{\mathrm{SH}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=& \hat{\boldsymbol{e}}^{\mathrm{SH}} \boldsymbol{u}_{\mathrm{r}}^{\mathrm{SH}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\hat{\boldsymbol{e}}^{\mathrm{SH}}\left[\hat{\boldsymbol{e}}^{\mathrm{SH}} \cdot\right.\\ &\left.\boldsymbol{u}_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)\right] \end{aligned} $ (63)

SV波的矢量波场为:

$ \begin{gathered} \boldsymbol{u}_{\mathrm{r}}^{\mathrm{SV}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\hat{\boldsymbol{e}}_{\mathrm{r}}^{\mathrm{SV}} \boldsymbol{u}_{\mathrm{r}}^{\mathrm{SV}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\hat{\boldsymbol{e}}_{\mathrm{r}}^{\mathrm{SV}}\left[\hat{\boldsymbol{e}}_{\mathrm{r}}^{\mathrm{SV}} \cdot\right. \\ \left.\boldsymbol{u}_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)\right] \end{gathered} $ (64)

利用地下介质模型为光滑背景模型下入射方程(42)和弹性波线性积分表达式(44), 结合上文的地震数据波形成像方法理论框架及其波形偏移的伴随求解方法公式(10), 将笛卡尔坐标系下的弹性波场转化为P波、SH波和SV波极化方向坐标系下的弹性波场, 再经过一些简单的推导就可以得到对弹性波阻抗相对扰动矩阵Erips中的元素Irmn(m, n=p, sh, sv)和反射体局部弹性反射系数矩阵Errps中的元素rmn(m, n=p, sh, sv)的近似估计。以PP和PSV波的波形偏移为例, 针对P波入射和P波、SV波反射的波场可以得到如下的弹性波阻抗相对扰动Irpp, Irsvp和局部弹性反射系数rpp, rsvp的弹性波反射数据波形偏移计算步骤和计算公式。

1) 地下入射弹性波场的计算:

$ \boldsymbol{L}_{\mathrm{b}} \boldsymbol{u}_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\boldsymbol{f}\left(t, \boldsymbol{x}_{\mathrm{s}}\right) $ (65)

2) 基于边值条件和波场逆时传播的地下反射波虚源近似反演:

$ \left\{\begin{array}{l} \bf{L}_{\mathrm{b}} \boldsymbol{u}_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=0 \\ \boldsymbol{u}_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\boldsymbol{u}_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \delta\left(\boldsymbol{x}_{\mathrm{g}}-\boldsymbol{x}\right) \end{array}\right. $ (66)

式中: ur(xg, t; xs)代表弹性反射波的炮集数据。

3) 对(65)式得到的入射波场ui(x, t; xs)和(66)式得到的反射波场ur(x, t; xs)进行波场转换, 得到它们在P波, SH波和SV波极化方向的入射波场uips(x, t; xs)=[${\boldsymbol{u}}{}_{\mathrm{i}}^{\mathrm{P}}$(x, t; xs), uiSH(x, t; xs), ${\boldsymbol{u}}{}_{\mathrm{i}}^{\mathrm{SV}}$(x, t; xs)]和反射波场urps(x, t; xs)=${\boldsymbol{u}}{}_{\mathrm{r}}^{\mathrm{P}}$(x, t; xs), urSH(x, t; xs), ${\boldsymbol{u}}{}_{\mathrm{r}}^{\mathrm{SV}}$(x, t; xs)。

4) 得到入射P波分量和反射P波、SV波分量间的分量对分量弹性波阻抗相对扰动Irpp, Irsvp的近似估计和局部弹性反射系数rpp, rsvp的近似估计, 即:

$ \begin{gathered} I_{\mathrm{r}}^{\mathrm{pp}}(\boldsymbol{x}, \sigma)=a_{\alpha}+a_{\rho}(1+\cos \sigma)-2 \frac{\beta_{\mathrm{b}}^{2}}{\alpha_{\mathrm{b}}^{2}}\left(a_{\rho}+a_{\beta}\right) \cdot \\ \sin ^{2} \sigma=\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \int_{0}^{T} \mathrm{~d} t u_{\mathrm{r}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \frac{\partial^{2} u_{\mathrm{i}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{gathered} $ (67)
$ \begin{gathered} I_{\mathrm{r}}^{\mathrm{svp}}(\boldsymbol{x}, \sigma)=\sin \sigma\left(a_{\rho}\left(\frac{\alpha_{\mathrm{b}}}{\beta_{\mathrm{b}}}+2 \cos \sigma\right)+2 a_{\beta} \cos \sigma\right)= \\ \frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \int_{0}^{T} \mathrm{~d} t u_{\mathrm{r}}^{\mathrm{SV}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \frac{\partial^{2} u_{\mathrm{i}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{gathered} $ (68)
$ \begin{gathered} r^{\mathrm{pp}}(\boldsymbol{x}, \sigma)= \frac{-1}{\alpha_{\mathrm{b}}(\boldsymbol{x}) \rho_{\mathrm{b}}(\boldsymbol{x})} \int_{0}^{T} \mathrm{~d} t u_{\mathrm{r}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)\cdot \\ \frac{\partial u_{\mathrm{i}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} \end{gathered} $ (69)
$ \begin{gathered} r^{\mathrm{svp}}(\boldsymbol{x}, \sigma)=\frac{-1}{\alpha_{\mathrm{b}}(\boldsymbol{x}) \rho_{\mathrm{b}}(\boldsymbol{x})} \int_{0}^{T} \mathrm{~d} t u_{\mathrm{r}}^{\mathrm{SV}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \cdot \\ \frac{\partial u_{\mathrm{i}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} \end{gathered} $ (70)

为了得到σ角度域的Irpp(x, σ), rpp(x, σ), Irsvp(x, σ), rsvp(x, σ), 在(67)式至(70)式进行右边项运算前, 需要将反射波场和入射波场转换到相应的角度域[18-20]。由上述计算过程可知, Irsvp(x, σ), rsvp(x, σ)的成像结果在σ=0的两侧会出现极性反转。由(67)式和(68)式可知, 利用Irpp(x, σ), Irsvp(x, σ)的角度域共成像点道集数据和它们的表达式可以进行类似于基于Zoeppritz方程的AVA纵、横波速度扰动、密度扰动和弹性波阻抗扰动分析和反演。如果在应用(67)式至(70)式前不将反射波场和入射波场转换到角度域, 则利用(67)式至(70)式得到的成像结果为Irpp(x, σ), rpp(x, σ), Irsvp(x, σ), rsvp(x, σ)在σ角度域的平均值Irpp(x), rpp(x), Irsvp(x), rsvp(x)。对于Irpp(x, σ), rpp(x, σ)而言, 它们的σ角度域的平均值Irpp(x), rpp(x)主要反映了P波阻抗的相对扰动和局部界面近法线方向的局部弹性反射系数。受极性反转的影响, Irsvp(x, σ), rsvp(x, σ)在σ角度域的平均值Irsvp(x), rsvp(x)存在不确定性。

为了得到具有高保真性和高分辨率的弹性波阻抗相对扰动和局部弹性反射系数反演结果, 根据上节的地震数据波形成像方法理论框架及代表最小二乘波形偏移的(9)式, 再结合上述的弹性波反射数据波形偏移计算公式, 可通过迭代方式得到求解Irmn(x, σ)的弹性波反射数据最小二乘波形偏移计算公式。由于角度域反偏移计算的困难, 本文的弹性波阻抗相对扰动最小二乘波形偏移成像只考虑Irpp(x, σ)在σ角度域平均值Irpp(x)的估计。

1) 初始化。

k=0, Irppk(x)=Irpp0(x)=0。

入射波场计算: 利用公式(65)计算入射波场。

反射波场逆时外推:

$ \left\{\begin{array}{l} L_{\mathrm{b}} \boldsymbol{u}_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=0 \\ \boldsymbol{u}_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\boldsymbol{u}_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \delta\left(\boldsymbol{x}_{\mathrm{g}}-\boldsymbol{x}\right) \end{array}\right. $ (71)

对波场ui(x, t; xs)和ur(x, t; xs)进行转换得到波场uips(x, t; xs)和urps(x, t; xs)。

2) 迭代过程。

计算分量位波场残差:

$ \begin{aligned} \delta \boldsymbol{u}_{\mathrm{r}}^{\mathrm{P}_{k}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=& \boldsymbol{u}_{\mathrm{r}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)-\rho_{\mathrm{b}}(\boldsymbol{x}) \bar{I}_{\mathrm{r}}^{\mathrm{pp}_{k}}(\boldsymbol{x}) \cdot \\ & \frac{\partial^{2} \boldsymbol{u}_{\mathrm{i}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{aligned} $ (72)

Irppk(x)的修正:

$ \begin{gathered} \bar{I}_{\mathrm{r}}^{\mathrm{pp}_{k+1}}(\boldsymbol{x})=\bar{I}_{\mathrm{r}}^{\mathrm{pp}_{k}}(\boldsymbol{x})+\frac{\alpha_{k}}{\rho_{\mathrm{b}}(\boldsymbol{x})} \int_{0}^{T} \mathrm{~d} t \delta \boldsymbol{u}_{\mathrm{r}}^{\mathrm{P}_{k}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \cdot \\ \frac{\partial^{2} \boldsymbol{u}_{\mathrm{i}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{gathered} $ (73)

k=k+1, 迭代结束。

同样地, 可通过迭代方式求解rpp(x)的估计的最小二乘波形偏移计算公式, 即:

1) 初始化。

k=0, rppk(x)=rpp0(x)=0。

入射波场计算: 利用公式(65)计算入射波场。

反射波场逆时外推: 利用公式(71)进行反射波场逆时外推。

对波场ui(x, t; xs)和ur(x, t; xs)进行转换得到波场uips(x, t; xs)和urps(x, t; xs)。

2) 迭代过程。

计算分量位波场残差:

$ \begin{gathered} \delta \boldsymbol{u}_{\mathrm{r}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\boldsymbol{u}_{\mathrm{r}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)+\rho_{\mathrm{b}}(\boldsymbol{x}) \alpha_{\mathrm{b}}(\boldsymbol{x})\cdot \\ \bar{r}^{\mathrm{pp}_{k}}(\boldsymbol{x}) \frac{\partial \boldsymbol{u}_{\mathrm{i}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} \end{gathered} $ (74)

rppk(x)的修正:

$ \begin{gathered} \bar{r}^{\mathrm{pp}_{k+1}}(\boldsymbol{x})=\bar{r}^{\mathrm{pp}_{k}}(\boldsymbol{x})-\frac{\alpha_{k}}{\rho_{\mathrm{b}}(\boldsymbol{x}) \alpha_{\mathrm{b}}(\boldsymbol{x})} \int_{0}^{T} \mathrm{~d} t \cdot\\ \delta \boldsymbol{u}_{\mathrm{r}}^{\mathrm{P}_{k}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \frac{\partial \boldsymbol{u}_{\mathrm{i}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} \end{gathered} $ (75)

k=k+1, 迭代结束。

由于Irsvp(x, σ), rsvp(x, σ)在σ角度域存在极性倒转, 因此不考虑它们的最小二乘波形偏移。

如果在上述的弹性波反射数据波形偏移计算步骤3)的波场转换得到的是P波、SH波和SV波的入射波场和反射波场${\boldsymbol{u}}{}_{\mathrm{i}}^{\mathrm{P}}$(x, t; xs), uiSH(x, t; xs), ${\boldsymbol{u}}{}_{\mathrm{i}}^{\mathrm{SV}}$(x, t; xs), ${\boldsymbol{u}}{}_{\mathrm{r}}^{\mathrm{P}}$(x, t; xs), urSH(x, t; xs), ${\boldsymbol{u}}{}_{\mathrm{r}}^{\mathrm{SV}}$(x, t; xs), 则步骤4)进行PP波和PSV波的矢量波场成像, 可以得到Irppv, Irsvpv, rppv, rsvpv。根据弹性波矢量波场的成像公式[21], 和$ \hat{\boldsymbol{e}}{}_{\mathrm{i}}^{\mathrm{P}}$·$ \hat{\boldsymbol{e}}{}_{\mathrm{r}}^{\mathrm{P}}$=-cosσ, $ \hat{\boldsymbol{e}}{}_{\mathrm{i}}^{\mathrm{P}}$·$\hat{\boldsymbol{e}}{}_{\mathrm{r}}^{\mathrm{SV}}$=-sinσ, 得到Irppv, Irsvpv, rppv, rsvpv为:

$ \begin{gathered} I_{\mathrm{r}}^{\mathrm{pp} v}(\boldsymbol{x}, \sigma)=\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \int_{0}^{T} \mathrm{~d} t \boldsymbol{u}_{\mathrm{r}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \cdot \frac{\partial^{2} \boldsymbol{u}_{\mathrm{i}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}{ }_{\mathrm{s}}\right)}{\partial t^{2}} \\ =-I_{\mathrm{r}}^{\mathrm{pp}}(\boldsymbol{x}, \sigma) \cos \sigma \end{gathered} $ (76)
$ \begin{gathered} I_{\mathrm{r}}^{\mathrm{svpv}}(\boldsymbol{x}, \sigma)=\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \int_{0}^{T} \mathrm{~d} t \boldsymbol{u}_{\mathrm{r}}^{\mathrm{sv}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \cdot \frac{\partial^{2} \boldsymbol{u}_{\mathrm{i}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \\ =-I_{\mathrm{r}}^{\mathrm{svp}}(\boldsymbol{x}, \sigma) \sin \sigma \end{gathered} $ (77)
$ \begin{gathered} r^{\mathrm{ppv}}(\boldsymbol{x}, \sigma)=\frac{-1}{\alpha_{\mathrm{b}}(\boldsymbol{x}) \rho_{\mathrm{b}}(\boldsymbol{x})} \int_{0}^{T} \mathrm{~d} t \boldsymbol{u}_{\mathrm{r}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)\cdot \\ \frac{\partial \boldsymbol{u}_{\mathrm{i}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t}=-r^{\mathrm{pp}}(\boldsymbol{x}, \sigma) \cos \sigma \end{gathered} $ (78)
$ \begin{gathered} r^{\mathrm{svpv}}(\boldsymbol{x}, \sigma)=\frac{-1}{\alpha_{\mathrm{b}}(\boldsymbol{x}) \rho_{\mathrm{b}}(\boldsymbol{x})} \int_{0}^{T} \mathrm{~d} t \boldsymbol{u}_{\mathrm{r}}^{\mathrm{SV}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)\cdot \\ \frac{\partial \boldsymbol{u}_{\mathrm{i}}^{\mathrm{P}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t}=-r^{\mathrm{svp}}(\boldsymbol{x}, \sigma) \sin \sigma \end{gathered} $ (79)

由上述的矢量波场成像计算公式可知, 矢量波场的PSV成像结果Irsvpv(x, σ), rsvpv(x, σ)在σ角度域不会出现极性倒转。

由上述的标量波、声波和弹性波数据的波形成像计算公式可知, 波形成像方法中的波形偏移、最小二乘波形偏移的计算过程与逆时偏移、最小二乘逆时偏移的计算过程相同, 计算量相当。对于地震数据的波形成像处理而言, 如果波形偏移成像处理的目标是散射体, 则将地震数据视为散射数据进行本文提出的散射数据波形偏移, 可以得到散射体的波形偏移成像结果; 如果偏移成像处理的目标是反射体(或反射面), 则将地震数据视为反射数据进行本文提出的反射数据波形偏移, 可以得到反射体(或反射面)的波形偏移成像结果。

上述提出的反射数据波形成像不仅能实现反射界面的准确成像, 还能得到地层波阻抗相对变化(岩性)的成像, 综合应用波阻抗相对扰动的角度域共成像点道集数据和波阻抗相对扰动的表达式还可以进行类似于基于Zoeppritz方程的AVA分析和反演。

5 数值试验

为了验证上述波形成像方法理论, 我们首先采用不考虑密度变化的标量波动方程合成地震数据进行标量波地震数据的波形成像试验, 然后再采用考虑密度变化的声波动方程合成地震数据进行声波地震数据的波形成像试验。

5.1 标量波地震数据的波形成像试验 5.1.1 Sigsbee2A模型试验

采用SEG提供的Sigsbee2A模型数据验证本文提出的散射数据波形偏移方法和反射数据局部反射系数波形偏移方法效果。图 1a为Sigsbee2A的地层模型, 该模型中不仅包含高速盐体、精细的地层, 还包含许多局部高速散射体和底部的水平高速反射面, 十分适合用于验证散射数据波形偏移和反射数据波形偏移方法的效果。图 1b为用于波形偏移试验的光滑速度模型, 模型大小为1067×1201个网格点, 水平和垂直方向的网格间距分别为22.86m和7.62m。该模拟数据的观测系统为海上拖缆观测系统, 共有500炮, 炮间距为45.72m, 拖缆长度为7932.42m, 检波点间距为22.86m, 各炮的道数不均, 最少为57道, 最多为348道, 时间采样间隔为8ms, 时间采样点数为1500。

图 1 Sigsbee2A模型 a 真实地层模型; b 用于波形偏移试验的光滑速度模型

分别应用常规逆时偏移方法(波形成像第一部分论述中的(1)式至(3)式)、本文提出的散射数据波形偏移方法计算公式(本文(14)式至(16)式)和反射数据局部反射系数波形偏移方法计算公式(本文(26-1)式至(29-1)式)对Sigbee2A模型数据进行试验。图 2为应用常规逆时偏移方法、反射数据波形偏移方法和散射数据波形偏移方法得到的Sigsbee2A模型偏移结果。

图 2 Sigbees2A模型偏移结果 a 常规逆时偏移方法; b 反射数据波形偏移方法; c 散射数据波形偏移方法

图 2中的偏移结果可以看出, 3种偏移方法都取得了不错的成像结果。对比图 2中的3种偏移结果可知, 对于高速散射体的成像, 散射数据波形偏移方法得到的散射体偏移结果要优于常规偏移方法和反射波形偏移方法的偏移成像结果; 对于反射面的成像, 反射数据波形偏移方法得到的反射面的偏移成像结果要优于常规偏移方法和散射波形偏移方法的偏移成像结果。

为了仔细对比常规逆时偏移、反射数据波形偏移和散射数据波形偏移对散射体和平反射面的成像结果的差异, 我们分别对图 1所示的红色方框中的高速散射体和蓝色方框中的高速平反射面的成像结果进行对比。图 3图 1红色方框中的高速散射体模型和分别利用逆时偏移方法、反射数据波形偏移方法、散射数据波形偏移方法得到的成像结果。由图 3可以清楚地看出, 散射数据波形偏移方法(其偏移结果见图 3d)对高速散射体成像(成像位置和零相位)准确; 逆时偏移方法(其偏移结果见图 3b)对高速散射体成像位置准确, 但出现了180°的相位差(反极性), 将高速散射体成像为低速散射体; 反射数据波形偏移方法(其偏移结果见图 3c)对高速散射体成像出现了90°的相位差, 致使成像位置变浅, 有1/4波长的深度误差, 高速散射体模型位于成像结果的波峰至波谷过渡带的中间位置。由图 3还可以看出, 散射数据波形偏移结果的分辨率最高, 而逆时偏移结果的分辨率最低。图 3的试验结果与本文的理论(见波形成像第一部分论述中逆时偏移方法中的(3)式和本文的(16)式(29-1)式)相符, 波形成像第一部分论述中的(3)式与本文的(16)式相差一个时间二阶导数, 即产生180°的相位差, 本文的(29-1)式与(16)式相差一个负号和时间一阶导数, 即产生-90°的相位差。本文的(16)式和(29-1)式中入射波场的时间导数有助于提高偏移结果的分辨率。

图 3 高速散射体模型(a)、逆时偏移(RTM)方法(b)、反射数据波形偏移(RWFM)方法(c)和散射数据波形偏移(SWFM)方法(d)成像结果

图 4图 1蓝色方框中的高速平反射面模型和分别利用逆时偏移方法、反射数据波形偏移方法、散射数据波形偏移方法得到的成像结果。由图 4可以清楚地看出, 采用反射数据波形偏移方法(其偏移结果见图 4c)得到的高速平反射面成像结果(成像位置和零相位)准确; 采用逆时偏移方法(其偏移结果见)图 4b得到的高速平反射面成像结果存在90°的相位差, 致使成像位置变浅, 并存在1/4波长的深度误差, 高速平反射面模型位于成像结果的波峰至波谷过渡带的中间位置; 采用散射数据波形偏移方法(其偏移结果见图 4d)对高速平反射面成像出现了-90°的相位差, 致使成像位置变深, 并存在1/4波长的深度误差, 高速平反射面模型位于成像结果的波谷至波峰过渡带的中间位置。图 4中的试验结果与本文的理论相符, 波形成像第一部分论述中的(3)式与本文的(28-1)式相差一个一阶时间导数, 即产生90°的相位差, 本文的(16)式与(29-1)式相差一个负号和时间一阶导数, 即产生-90°的相位差。

图 4 高速平反射面模型(a)、逆时偏移(RTM)方法(b)、反射数据波形偏移(RWFM)方法(c)和散射数据波形偏移(SWFM)方法(d)成像结果

对比上述偏移成像结果可知, 将本文提出的散射数据波形偏移方法应用于散射体成像、反射数据波形偏移方法应用于反射面成像, 均可以取得理想的效果, 它们的成像结果相对于逆时偏移结果不仅分辨率高, 而且位置准确、相位(极性)正确。散射数据波形偏移方法有利于散射体的准确成像, 反射数据波形偏移方法有利于反射面的准确成像。

5.1.2 反射数据的局部反射系数最小二乘波形偏移试验

为了验证本文提出的反射数据局部反射系数最小二乘波形偏移方法, 我们截取Sigsbee2A速度模型(图 5a)的一部分进行试验。对于图 5a的速度模型采用标量波动方程的有限差分方法合成地震数据。本次试验合成了300炮地震数据, 炮点范围是0~6km, 炮间距20m, 为中间放炮两边接收, 每炮由601道接收, 道间距10m, 时间采样间隔是0.5ms, 时间采样点数为3000, 震源为主频18Hz的Ricker子波。由于图 5a的速度模型为层状结构, 因此合成地震数据主要为反射数据。图 5b为用于最小二乘波形偏移试验的光滑速度模型, 它是对图 5a的速度模型进行光滑处理得到的。

图 5 部分Sigsbee2A模型及其波形偏移结果、最小二乘波形偏移结果和模型的理论垂向反射率 a 真实速度模型; b 光滑速度模型; c 波形偏移结果; d 35次迭代的最小二乘波形偏移结果; e 模型的真实垂向反射率

在试验中, 我们首先应用反射数据局部反射系数波形偏移方法对合成的炮道集数据进行局部反射系数波形偏移, 图 5c为得到的波形偏移结果。由图 5c可以看出, 将局部反射系数波形偏移方法应用于该模型可获得很好的成像结果。图 5d为应用本文提出的反射数据局部反射系数最小二乘波形偏移方法得到的偏移结果, 共进行了35次迭代。对比图 5c图 5d可以看出, 图 5d的最小二乘波形偏移结果的分辨率相较于图 5c的波形偏移结果明显得到了提高, 并且层位面与断层面的成像更清晰。由图 5e所示的模型真实垂向反射率可以看出, 图 5d的最小二乘波形偏移结果的振幅保真性相较于图 5c的波形偏移结果也存在明显提升。

5.2 声波地震数据的波形成像

为了验证声波地震数据的波形成像方法, 采用考虑密度变化的声波动方程合成地震数据进行声波阻抗相对扰动的声波反射数据波形偏移试验。

图 6a为用于合成地震数据的速度模型, 它是在常速为2000m/s的背景模型上叠加一些理论速度值在(-100m/s, 100m/s)变化的随机层。图 6b为用于合成地震数据的密度模型, 它是在常密度为2000kg/m3的背景模型上叠加一些理论密度值在(-50kg/m3, 50kg/m3)变化的随机层。图 6c为利用图 6a图 6b的速度、密度扰动、背景速度和密度得到的声波阻抗相对扰动模型。利用声波方程的有限差分方法合成地震数据, 本试验合成了50炮地震数据, 炮点范围为0~2km, 炮间距为40m, 为中间放炮两边接收, 每炮由201道接收, 道间距为10m, 时间采样间隔为1ms, 时间采样点数为3000, 震源为主频20Hz的Ricker子波。由于图 6的模型为层状结构, 因此合成的地震数据为反射数据。用于波形成像的速度和密度模型分别为速度为2000m/s的常速模型和密度为2000kg/m3的常密度模型。

图 6 用于波形成像试验的速度模型(a)、密度模型(b)和相对声波阻抗相对扰动模型(c)

利用(26)式至(28)式所表示的波阻抗相对扰动波形偏移计算公式对合成地震反射数据进行波阻抗相对扰动的波形偏移。图 7为波形偏移得到的结果, 其中图 7a为波形偏移的角度域叠加结果, 图 7bx=1000m处的角度域共成像点道集, 图 7cx=1000m处角度域共成像点道集零度角对应的波阻抗相对扰动波形偏移与波阻抗相对扰动模型值对比结果(右上)和角度域共成像点道集中60°角对应的速度相对扰动波形偏移与速度相对扰动模型值对比结果(右下)。对比图 6图 7中的波形偏移结果可以看出, 本文提出的声波阻抗相对扰动波形偏移方法可以很好地实现对声波阻抗相对扰动的成像。由图 6c图 7a验证了在(28)式中不进行波场角度分解而直接得到角度域的平均波阻抗相对扰动是一种快速的波阻抗相对扰动成像方法, 可用来对地层波阻抗变化进行成像。由图 7c中的对比曲线可以看出, 波阻抗相对扰动的角度域共成像点道集中的小角度道和大角度道可分别作为波阻抗相对扰动和速度相对扰动的近似估计。图 7b的波阻抗相对扰动角度域共成像点道集为进一步的角度域叠前岩性分析和反演提供了基础数据。

图 7 波阻抗相对扰动的波形偏移结果 a 角度域叠加结果; b x=1000m处的角度域共成像点道集; c x=1000m处角度域共成像点道集中零度(右上)和60°(右下)角所对应波阻抗相对扰动波形偏移和速度相对扰动波形偏移与模型值对比结果

上述合成数据的数值试验结果验证了本文提出的标量波和声波的散射数据波形偏移方法、反射数据波形偏移方法和最小二乘波形偏移方法的正确性和有效性。根据波形成像方法的计算公式, 我们认为本文提出的波形成像方法具有很强的实用性。弹性波数据的波形成像数值试验结果将在后续的文章中展示。

6 结论

对于波形成像的光滑背景模型, 在地震波不发生散射、反射和波型转换的假定下, 基于地震数据的线性正演表达, 本文首先给出了地震数据波形成像的地震波通用传播方程、散射数据波形成像的通用散射方程和反射数据波形成像的通用反射方程, 然后利用线性反演理论建立了面向散射数据和反射数据的地震数据波形成像方法理论框架。散射数据波形成像是对散射体物性参数扰动的线性反演, 适合相对于地震波长较小尺寸的散射体和绕射体的成像。反射数据波形成像是对反射体波阻抗相对扰动或反射体边界局部反射系数的线性反演, 适合相对于地震波长较大尺寸的反射体物性和地层反射界面的成像。反射数据波形成像不仅能实现对反射界面位置的准确成像, 还能得到对地层波阻抗相对变化的岩性成像。在波形成像中, 如果利用波传播算子的伴随算子进行地震波的反传播, 则波形成像退化为波形偏移, 如果利用波传播算子的最小二乘逆算子进行地震波的反传播, 则波形成像退化为最小二乘波形偏移。波形偏移计算量小稳定性好, 但其偏移成像结果会因数据观测系统和地下复杂构造的影响在保真性和分辨率方面存在不足。最小二乘波形偏移计算量大, 但能有效校正数据观测系统和地下复杂构造的影响, 因此其偏移成像结果保真性好、分辨率高。我们给出了标量波散射数据、声波反射数据和弹性波反射数据波形成像的计算式, 并进行了标量波地震数据和声波地震数据的波形成像数值试验, 取得了与理论相一致的试验结果。我们提出的波形成像方法弥补了当前地震数据偏移成像方法理论上存在的不足, 可以准确地利用地震数据的波形信息, 修正当前逆时偏移方法存在的相位误差, 得到的偏移结果没有相位误差, 分辨率也有所提高。与常规的逆时偏移和最小二乘逆时偏移的计算相比, 波形偏移和最小二乘波形偏移仅增加了对入射波场的时间导数计算。反射体波阻抗相对扰动随反射开角、岩性变化之间的定量表达式和反射数据波形成像得到的波阻抗相对扰动角度域共成像点道集为角度域成像数据的岩性分析和反演提供了方法理论和数据基础, 这是我们下一步研究工作的方向。

参考文献
[1]
DEVANEY J. Mathematical foundations of imaging, tomography and wavefield inversion[M]. Cambridge: Cambridge University Press, 2012: 267-371.
[2]
KIRSCH A. An introduction to the mathematical theory of inverse problems[M]. 2nd ed. Switzerland: Springer, 2011: 1-346.
[3]
AKI K, RICHARD P. Quantitative seismology[M]. 2nd ed. Sausalito: University Science Books, 2002: 145-178.
[4]
BLEISTEIN N, COHEN J, STOCKWELL J. Mathematics of multidimensional seismic imaging, migration and inversion[M]. New York: Springer Verlag, 2001: 367-489.
[5]
NOCEDA J, WRIGHT S. Numerical optimization[M]. 2nd ed. New York: Springer Science & Business Media, 2006: 1023-1156.
[6]
ZHOU H, CHEN S, AND ZHOU L, et al. Reflected wave least squares reverse time migration with angle illumination compensation[J]. Acta Geophysica, 2019, 67(3): 1551-1561.
[7]
SAVA P, FOMEL S. Angle-domain common-image gathers by wavefield continuation methods[J]. Geophysics, 2003, 68(3): 1065-1074.
[8]
XU S, ZHANG Y, TANG B. 3D angle gathers from reverse time migration[J]. Geophysics, 2011, 76(2): S77-S92.
[9]
STOLT R, WEGLEIN B. Seismic imaging and inversion: Application of linear inverse theory[R]. Houston: Cambridge University Press, 2012
[10]
WAPENAAR C, KINNEGING N, BERKHOUT A. Principle of prestack migration based on the full elastic two-way wave equation[J]. Geophysics, 1987, 52(2): 151-173.
[11]
YAN J, SAVA P. Isotropic angle domain elastic reverse time migration[J]. Geophysics, 2008, 73(6): S229-S239.
[12]
DELLINGER J, ETGEN J. Wave-field separation in two-dimensional anisotropic media[J]. Geophysics, 1990, 55(7): 914-919.
[13]
ZHANG Q, MCMECHAN G. 2D and 3D elastic wavefield vector decomposition in the wavenumber domain for VTI media[J]. Geophysics, 2010, 75(3): D13-D26.
[14]
WANG W, MCMECHAN G. Vector-based elastic reverse time migration[J]. Geophysics, 2015, 80(6): S245-S258.
[15]
周熙焱, 常旭, 王一博, 等. 基于纵横波解耦的三维弹性波逆时偏移[J]. 地球物理学报, 2018, 61(3): 1038-1052.
ZHOU X Y, CHANG X, WANG Y B, et al. 3D elastic reverse time migration based on P-and S-wave decoupling[J]. Chinese Journal of Geophysics, 2018, 61(3): 1038-1052.
[16]
YOON K, MARFURT K. Reverse-time migration using the Poynting vector[J]. Exploration Geophysics, 2006, 37(1): 102-107.
[17]
ZHU H. Elastic wavefield separation based on the Helmholtz decomposition[J]. Geophysics, 2017, 82(2): S173-S183.
[18]
ZHANG Q, CHANG M. Direct vector-field method to obtain angle-domain common-image gathers from isotropic acoustic and elastic reverse time migration[J]. Geophysics, 2011, 76(5): WB135-WB149.
[19]
DUAN Y, SAVA P. 3D angle decomposition for elastic reverse time migration[J]. Expanded Abstracts of 86th Annual Internat Mtg, 2016, 2356-2361.
[20]
TANG C, MCME C. Multidirectional-vector-based elastic reverse time migration and angle-domain common-image gathers with approximate wavefield decomposition of P-and S-waves[J]. Geophysics, 2018, 83(1): S57-S79.
[21]
WANG C, CHENG J, ARNTSEN B. Scalar and vector imaging based on wave mode decoupling for elastic reverse time migration in isotropic and transversely isotropic media[J]. Geophysics, 2016, 81(5): S383-S398.