一次波(一次反射波、一次散射波)地震数据的波动方程偏移成像起源于20世纪70年代初期^[1-3], 是当前地震数据处理中最重要的方法技术之一, 也是获得地下三维构造图像的主要技术手段^[4-7]。然而, 随着油气勘探开发目标的日趋复杂化和精细化(如复杂构造岩性油气藏和非常规油气藏的勘探开发), 需要我们发展高分辨、高保真的反演成像方法。宽方位、宽频带、高密度数据采集方法技术的应用, 为我们利用波形(即振幅、相位)信息进行高分辨、高保真的反演成像方法研究提供了数据基础。地震数据全波形反演方法理论^[8-10]的发展也为我们开展高分辨、高保真的反演成像方法研究提供了借鉴。此外, 地震数据的波动方程偏移成像方法理论也需要从当前的构造成像走向以地层物性参数和储层参数为目标的岩性成像。

地震数据是地下介质变化情况的客观反映, 不仅包含地震波的走时信息, 还含有地震波的波形信息。地震勘探中地震波的频带具有特定的范围, 相应的地震波长也有一定的范围, 因此在利用地震数据研究地下介质变化时, 须考虑地震波长的长短。地震数据的正演表达是构建地震波反演成像方法的基础, 不同的正演表达会导致不同的反演目标和反演方法。从地震数据的全波形反演可知, 要实现地震数据的岩性成像不仅需要利用好地震数据的走时信息, 更要准确利用地震数据的波形信息。为了在地震数据的偏移成像中准确地利用波形信息, 使当前的构造成像走向以地层物性参数和储层参数为目标的岩性成像, 我们必须从波动理论出发为地震数据构建一个严谨的线性正演表达。

基于上述认识, 我们从地震波的控制方程(波动方程)出发, 利用扰动理论, 并根据地下非均匀体的尺寸或其物理性质空间变化特征尺度与地震波长之间的相对大小关系, 将非均匀体划分为产生散射波的散射体和产生反射波的反射体, 推导出控制散射波的散射波方程和控制反射波的反射波方程, 然后在一定的近似条件下分别得到散射地震数据和反射地震数据的线性正演表达式, 再利用线性反演理论构建基于散射地震数据和反射地震数据的线性正演表达式的地震数据波形成像方法理论。上述工作也是对我们近年相关研究工作^[11-19]的补充、完善和系统化。

1 当前地震数据波动方程偏移的主要不足

根据CLAERBOUT^[1]提出的地震数据波动方程偏移成像理论和BERKHOUT^[20-21]提出的有关地震波传播的WRW概念模型, 在给定满足地震波运动学的准确的光滑模型(也称为偏移的宏观模型)下, 当前的地震数据波动方程逆时偏移主要由两步组成: 第一步是波场外推, 即利用宏观模型下的波动方程分别对震源波场和记录波场进行顺时外推和逆时外推; 第二步是波场成像, 即利用CLAERBOUT提出的时间一致性成像原理^[22-26]建立的成像条件(公式)对外推得到的记录波场和震源波场进行成像, 得到由地下反射系数估计的偏移成像结果, 以此实现对反射面或散射(绕射)体空间位置的构造成像。经过近50年的发展, 当前的地震数据偏移成像方法主要包括: ①利用描述波场传播的Kirchhoff积分公式的射线Kirchhoff偏移方法和射线Beam偏移方法^[27-32]; ②在Born近似和WKBJ近似建立的散射数据表达式或Kirchhoff近似建立的基于反射系数的反射数据表达式的基础上, 利用渐进反演理论构建的积分形式的射线Kirchhoff真振幅偏移方法^[33-36]; ③利用单程波波动方程和Claerbout成像条件的单程波偏移方法^[37-40]; ④利用双程波波动方程和Claerbout成像条件的逆时偏移方法^[41-43]; ⑤借助射线Kirchhoff真振幅偏移方法或利用采集照明补偿的波动方程真振幅偏移方法^[44-53]; ⑥利用Born近似的地震数据正演方程建立的波动方程最小二乘偏移方法^[54-55]; ⑦利用反射波方程构建的波动方程最小二乘偏移方法^{[13, 17]}。

由于利用了准确的波动方程进行波场外推, 逆时偏移被认为是当前理论上最为完善的偏移方法。对于标量波方程的地震数据常规逆时偏移方法有以下计算公式。

1) 光滑偏移速度模型下的地下入射波场计算。

$ \left\{\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\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) $

(1)

2) 光滑偏移速度模型下利用逆时外推重建地下反射波场。

$ \left\{\begin{array}{l} \left\{\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\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)=\mathrm{d}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \delta\left(\boldsymbol{x}_{\mathrm{g}}-\boldsymbol{x}\right) \end{array}\right. $

(2)

3) 利用成像条件(反射系数成像条件)进行波场成像。

$ \begin{gathered} I(\boldsymbol{x})=\int\limits_{0}^{T} \mathrm{~d} t u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) / u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \text { 或退化为 } \\ I(\boldsymbol{x})=\int\limits_{0}^{T} \mathrm{~d} t u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \end{gathered} $

(3)

(1) 式和(2)式分别表示利用光滑偏移速度模型下的波动方程进行震源波场的顺时外推和记录波场的逆时外推; (3)式为基于Claerbout成像原理的成像公式, 即小角度入射条件下, 无穷平边界和平面波的镜面反射系数成像公式。在上述公式中, ${\boldsymbol{\nabla} ^2}$为Laplace算子, 有${\boldsymbol{\nabla} ^2} = {\partial ^2}/\partial {x^2} + {\partial ^2}/\partial {y^2} + {\partial ^2}/\partial {z^2}$; x_s表示震源位置坐标矢量; v_b(x)为光滑偏移速度模型; u_i(x, t; x_s)表示光滑偏移速度模型下的地下入射波场, 也称为震源波场; f(t)表示震源时间函数; u_r(x, t; x_s)代表光滑偏移速度模型下的地下反射波场; d(x_g, t; x_s)表示检波点坐标x_g处的地震数据; T表示最大记录时间; I(x)表示偏移成像得到的地下反射系数估计。

当前地震数据波动方程偏移成像方法存在的不足主要有以下3方面。①不区分散射数据和反射数据, 将散射和反射视为同义词, 用同一种公式对散射数据和反射数据进行偏移成像^{[7, 11]}。我们知道, 在波动力学中, 散射波与反射波具有不同的波场特征, 散射波是由局部非均匀体形成的二次源产生的, 而反射波是由空间上具有一定延续度的非均匀体形成的二次源产生的, 可视为散射波的叠加, 所以不能将基于散射概念建立的散射数据表达应用于反射数据的偏移成像, 同样也不能将基于平面波平界面的反射概念的波场关系应用于散射数据的偏移成像。②成像公式没有正确考虑散射波、反射波的不同产生机制, 公式(3)所表示的成像公式仅适合平面波小角度入射到无穷大平反射面的特殊情况。③缺乏与波动方程偏移成像相适应的地震数据严谨正演表达式, 使得偏移成像不能准确地利用地震数据的波形信息, 得到的偏移成像结果会出现相位误差。在当前的偏移成像中, 对于地震数据的正演表达有基于介质模型参数扰动的Born近似表达^{[33, 56-57]}和基于反射面反射率的Kirchhoff近似表达^[58-60]。前者适合相对地震波长为小尺度的局部非均匀体产生的散射数据, 而不适用于相对地震波长为大尺度的非均匀体产生的反射数据, 后者主要用来表达反射面产生的反射数据, 但也存在问题(见下文)。基于地震数据Born近似表达而建立的偏移成像方法得到的成像结果是有关模型参数扰动的近似估计, 如果非均匀体为小尺度的局部散射体, 则该近似估计可有效反映散射体的位置; 如果非均匀体为大尺度的反射体, 则该近似估计难以准确反映反射体的边界, 因此基于模型参数扰动Born近似的地震数据表达是不适合反射数据的构造(反射面位置)成像。根据上述认识, 我们认为在偏移成像中对于地震散射数据和反射数据应该有不同的正演表达公式, 并由此建立其不同的偏移成像计算公式。

在应力与应变间满足线性近似的假定下, 地震震源激发的地震波在地下的运动可以用下述二阶微分算子方程描述:

$ L(m) u=f $

(4)

式中: L为与地下模型物性参数m(也称为弹性参数)有关的二阶偏微分算子, 也称为波动算子; u表示地震波在地下运动状态的状态变量, 也称为地震波场; f为激发地震波的震源函数。方程(4)可描述震源和地下模型m变化产生的各种波, 如直达波、散射波(绕射波)、反射波、透射波和转换波(对于弹性波动方程)等等, 也是标量波方程、声波方程和弹性波方程等的抽象形式, 因此被称为一般(或通用)波动方程, 是地震波场模拟、地震数据全波形反演的基础方程。

给定满足地震波运动学的准确的光滑偏移模型m_b(也称为用于偏移的宏观模型), 逆时偏移中用于波场传播的通用波动方程为:

$ L\left(m_{\mathrm{b}}\right) u_{\mathrm{i}}=f $

(5)

由于m_b的光滑性, 方程(5)和其对应的齐次方程只能描述地震波的传播, 而不能描述地震波的散射、反射和波型转换, 因此入射波场u_i中不包含散射波场、反射波场和转换波场。

逆时偏移中Claerbout成像公式(条件)(3)的反射波场除以入射波场的成像公式是平面波小角度入射到无穷大平反射面的反射系数计算公式。对于无穷大平边界和平面波, 在入射角小于临界角的条件下, 反射系数为与频率无关的实数^[61-62], 有:

$ u_{\mathrm{r}}=R u_{\mathrm{i}} $

(6)

式中: R为与入射角有关的镜面反射系数。如果入射波不是平面波或反射面不是无穷大平边界或入射角大于临界角, 则公式(6)不成立。因为公式(6)中的反射系数R为实数, 则要求式中的反射波与入射波应有相同的波形(即有相同的相位)。我们将(6)式所表示的反射波与入射波之间的关系称为镜面反射方程, 它也被用于反射数据的Kirchhoff近似表达和基于Kirchhoff近似表达的偏移成像中^[58-59]。因此, 我们认为当前广泛应用的Claerbout成像条件和Kirchhoff近似对于非无穷大平边界和非平面波地震数据的偏移成像存在不足, 它无法考虑入射波与反射波、散射波之间的相位差异。

在当前的偏移成像方法中, 将方程(5)和方程(6)作为描述地震波传播和反射(散射)的通用方程, 并将它们联合作为偏移成像方法中地震数据的正演表达。由上述的分析可知, 方程(5)和方程(6)是不同条件下相互独立的波动方程, 它们的联合不能作为非无穷大平边界和非平面波情况下地震数据的正演表达。因此, 我们认为当前的偏移成像缺乏与之适应的地震数据严谨正演表达公式, 致使其不能准确地利用地震数据的波形信息。也由于方程(6)的严格假定条件, 导致当前的偏移成像方法对于散射数据和反射数据的偏移成像通常都会出现相位误差。

为了使偏移成像方法更适合地下实际情况和能准确地利用地震数据波形信息, 弥补当前地震数据偏移成像方法理论中存在的不足和修正当前波动方程偏移成像方法对散射数据和反射数据偏移成像存在的相位误差, 有必要为偏移成像方法提供与之相适应的地震数据严谨正演表达, 以满足地下复杂构造油气藏、岩性油气藏和非常规油气藏勘探开发对高分辨、高保真偏移成像的要求。本文在满足地震波运动学的准确的光滑模型下, 首先将扰动法应用于通用波动方程, 得到描述非均匀体产生的扰动波场的波动方程。依据地震波长与非均匀体大小或非均匀体物理性质空间变化特征尺度之间的相对大小关系, 将非均匀体划分为在空间上具有有限延续度的产生散射波的散射体或空间上具有一定延续度的产生反射波的反射体, 相应的描述扰动波场的波动方程转化为散射波动方程和反射波动方程。利用小扰动假定、Born近似(或一次波近似)和高频近似(地震波长相对于非均匀体大小或其物理性质空间变化特征尺度为小量), 得到基于散射体物性参数相对扰动的散射数据线性正演表达和基于反射体波阻抗相对扰动的反射数据线性正演表达。为了描述反射体边界对反射波的作用, 在高频近似条件下通过定义反射体波阻抗相对扰动沿入射波传播方向的方向导数为反射体边界的局部反射系数, 推出基于局部反射系数的反射数据线性正演表达。最后给出标量波、声波和弹性波方程下散射数据和反射数据的具体线性正演表达式。

2 地震波的散射和反射及其通用方程

地震勘探中的一次地震波在地下传播的物理过程直观上可表述为: 地表震源激发的入射波场向下传播, 入射波场与非均匀体作用形成二次源, 并产生次生波场(散射波场、反射波场等), 次生波场向上传播到地表被检波器接收。在上述的直观表述中, 地震波的传播和散射(反射)是相互独立的, 但在描述地震波场运动状态的通用波动方程(4)中, 地震波的传播与散射(反射)是交织在一起的。这是因为方程(4)中的模型m不仅控制了地震波的传播, 也控制了地震波的散射(反射)。如果方程(4)中的模型m为光滑模型, 则波动方程(4)就不能描述地震波的散射、反射和波型转换, 即波动方程(5)中的波场u_i不包含散射、反射和转换波。为了研究光滑模型下地震波的波场模拟、反演和偏移成像方法, 如反射波形反演、反射波阻抗反演、反射数据的偏移成像(包括最小二乘偏移成像)方法等, 我们首先需要建立光滑模型下描述地震波的散射和反射的波动方程, 并以此得到适合地震数据偏移成像的地震数据严谨正演表达公式。

对通用波动方程(4)应用扰动法^[63], 即分别引入模型m和波场u的扰动表达: m=m_b+δm, u=u_i+δu, 其中, m_b为光滑的背景模型(简称为光滑模型), δm为模型扰动(非均匀体的表示), u_i为光滑模型m_b下的入射波场(也称为背景波场), δu为模型扰动δm产生的扰动波场。将模型与波场的扰动表达代入方程(4), 得到与其对应的扰动方程, 即非均匀体产生的扰动波场所满足的波动方程。

$ L\left(m_{\mathrm{b}}\right) \delta u=-\delta L\left(m_{\mathrm{b}}, \delta m\right) u $

(7)

式中: L(m_b)为光滑模型下的波动算子; δL(m_b, δm)为扰动算子, 有δL(m_b, δm)=L(m)-L(m_b)=L(m_b+δm)-L(m_b)。

地震数据是一种具有特定频带范围的波形数据, 因此利用地震数据研究地下非均匀体结构必须考虑地震波长与非均匀体δm的空间尺寸或者非均匀体δm物性空间变化的特征尺度之间的相对大小关系。如果非均匀体δm的尺寸a或者其物性空间变化特征尺度a相对于地震波场u的波长λ比较小, 通常假定a/λ≤1, 则非均匀体δm相对于波长在空间上只有有限的延续度, 可视为产生散射波的局部散射体, 相应的波动方程(7)化为散射波所满足的波动方程, 即:

$ L\left(m_{\mathrm{b}}\right) u_{\mathrm{s}}=-\delta L_{\mathrm{s}}\left(m_{\mathrm{b}}, \delta m\right) u $

(8)

式中: u_s表示散射体产生的散射波场; δL_s(m_b, δm)表示与扰动δm或相对扰动δm/m_b有关的散射算子, 它与波场u的相互作用是激发散射波的虚源。方程(8)也称为描述散射波的通用方程。由于δm相对波长比较小, 因此方程(8)的右端源可视为点源, 产生的散射波没有特定的传播方向。利用背景模型下波动方程(5)所对应的Green函数, 将方程(8)表示为时空域积分形式, 有:

$ u_{\mathrm{s}}=-\int_{\varOmega} \mathrm{d} v g *{ }_{t} \delta L_{\mathrm{s}}\left(m_{\mathrm{b}}, \delta m\right) u $

(9)

式中: Ω为散射波源的分布空间; g为光滑模型m_b下波动方程(5)的Green函数; “*_t”表示时间褶积。

如果非均匀体δm的尺寸a或者其物性空间变化特征尺度a相对于地震波场u的波长λ比较大, 通常假定a/λ>1, 则非均匀体δm相对于波长在空间上具有一定的延续度, 可视为产生反射波的反射体, 相应的波动方程(7)化为反射波所满足的波动方程, 即:

$ L\left(m_{\mathrm{b}}\right) u_{\mathrm{r}}=-\delta L_{\mathrm{r}}\left(m_{\mathrm{b}}, I_{\mathrm{r}}, \sigma\right) u $

(10)

式中: u_r表示反射体产生的反射波场; δL_r(m_b, I_r, σ)表示与波阻抗相对扰动I_r和角度σ有关的反射算子, 它与波场u的相互作用是激发反射波的虚源; σ表示入射方向与反射方向之间的开角, 它是由波动算子中物理参数的空间导数和波场的空间导数在高频近似条件下的表达式引入的; 波阻抗相对扰动I_r不同于散射体的物性参数相对扰动δm/m_b, 它不仅与m_b, δm有关, 还与反射开角σ有关。方程(10)也称为描述反射波的通用方程。由于δm相对于波长在空间上具有一定延续度, 因此方程(10)的右端源可视为由点源集成的局部平面波源, 产生的反射波可视为具有特定传播方向的局部平面波。利用背景模型下波动方程(5)所对应的Green函数, 同样可将方程(10)表示为积分形式, 有:

$ u_{\mathrm{r}}=-\int_{\varOmega} \mathrm{d} v g *{ }_{t} \delta L_{\mathrm{r}}\left(m_{\mathrm{b}}, I_{\mathrm{r}}, \sigma\right) u $

(11)

式中: Ω为局部平面波源的分布空间。

我们根据非均匀体δm的尺寸a或者其物性空间变化特征尺度a与地震波长λ之间的相对大小关系, 将δm分别划分为产生散射波的散射体和产生反射波的反射体, 并分别得到了描述散射波和反射波的通用方程(公式(8)和公式(10)), 以及它们的积分形式(公式(9)和公式(11))。

利用对δm的小扰动假定和一阶Born近似(忽略多次散射波), 方程(8)退化为描述一次散射波u_s的非齐次波动方程, 即:

$ L\left(m_{\mathrm{b}}\right) u_{\mathrm{s}}=-\delta L_{\mathrm{s}}\left(m_{\mathrm{b}}, \delta m\right) u_{\mathrm{i}} $

(12)

相应的方程(9)退化为:

$ u_{\mathrm{s}}=-\int_{\varOmega} \mathrm{d} v g *{ }_{t} \delta L_{\mathrm{s}}\left(m_{\mathrm{b}}, \delta m\right) u_{\mathrm{i}} $

(13)

利用对I_r的小值假定和一次反射波近似(忽略多次反射波), 方程(10)退化为描述一次反射波u_r的非齐次波动方程, 即:

$ L\left(m_{\mathrm{b}}\right) u_{\mathrm{r}}=-\delta L_{\mathrm{r}}\left(m_{\mathrm{b}}, I_{\mathrm{r}}, \sigma\right) u_{\mathrm{i}} $

(14)

相应的方程(11)退化为:

$ u_{\mathrm{r}}=-\int_{\varOmega} \mathrm{d} v g *{ }_{t} \delta L_{\mathrm{r}}\left(m_{\mathrm{b}}, I_{\mathrm{r}}, \sigma\right) u_{\mathrm{i}} $

(15)

如果反射体δm内部的物性变化很小或不变化, 则反射主要是由反射体边界产生^[64]。为了描述反射体边界对反射的作用和满足对反射体边界成像的需要, 本文将产生反射的波阻抗相对扰动I_r沿入射波传播方向I的方向导数定义为反射体边界局部切平面的局部反射系数, 有:

$ r=\frac{\partial}{\partial I} I_{\mathrm{r}} $

(16)

式中: r表示局部反射系数, 是局部平面波入射到反射体边界位置x处局部切平面上的局部反射系数, 它不仅可以表征反射面的空间位置, 还能反映反射面上的物性变化。I_r在数学形态上是一个阶跃函数, 而r在数学形态上是一个沿反射体边界分布的δ函数。图 1为分界面上点x处局部切平面的反射示意图。图中L_p为点x处的局部切平面; r为点x处的局部反射系数; u_i和u_r分别为具有局部平面波性质的入射波和反射波; σ为点x处的反射开角。

假定地震波长相对于光滑模型m_b空间变化特征尺度满足高频近似条件(即光滑模型m_b的空间变化相对于地震波长尺度可视为缓慢变化甚至不变化), 则公式(16)中的入射波传播方向导数在波数域对应为ik_i, k_i为入射波传播方向的波数, 有k_i=ω/[v_b(x)], 也称为入射波传播方向在x处的局部波数, ω为地震波的圆频率, v_b(x)为点x处光滑模型的速度。因此公式(16)在高频近似条件下的波数域表达式为:

$ r=i k_{\mathrm{i}} I_{\mathrm{r}}=\frac{i \omega}{v_{\mathrm{b}}(\boldsymbol{x})} I_{\mathrm{r}} $

(17)

公式(16)和公式(17)定义的局部反射系数r是地下物性参数相对扰动空间变化的反映。不同于公式(6)中与频率无关的无量纲镜面反射系数R, 本文定义的局部反射系数r与传播方向有关且具有长度倒数的量纲, 它是高频近似条件下, 局部平面入射波作用于局部反射平面的反射系数, 所以与入射波传播方向上的局部波数有关(在形式上与频率有关)。

利用公式(16)和公式(17)定义的局部反射系数r, 可将方程(14)和方程(15)分别改写为:

$ L\left(m_{\mathrm{b}}\right) u_{\mathrm{r}}=-\delta L_{\mathrm{r}}\left(m_{\mathrm{b}}, r, \sigma\right) u_{\mathrm{i}} $

(18)

$ u_{\mathrm{r}}=-\int_{\varOmega} \mathrm{d} v g *{ }_{t} \delta L_{\mathrm{r}}\left(m_{\mathrm{b}}, r, \sigma\right) u_{\mathrm{i}} $

(19)

式中: δL_r(m_b, r, σ)表示与局部反射系数r有关的反射算子。方程(18)的右端源也是一个局部平面波源。

3 地震数据线性正演表达的具体形式

利用前文推导出的基于散射体物性参数扰动的一次散射波通用方程(12)和方程(13)以及基于反射体波阻抗相对扰动的一次反射波通用方程(14)和方程(15)或基于反射体边界局部反射系数的一次反射波通用方程(18)和方程(19), 结合具体的标量波、声波和各向同性弹性波方程, 我们给出散射数据和反射数据线性正演表达的具体形式。

3.1 标量波地震数据的线性正演表达

如果通用波动方程(4)中模型物性参数m仅为速度v(x), 则方程(4)为非齐次标量波动方程, 即:

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

(20)

根据扰动方法, 有v(x)=v_b(x)+δv(x), u(x, t; x_s)=u_b(x, t; x_s)+δu(x, t; x_s), 其中, v_b(x)为光滑速度模型, δv(x)为速度模型扰动, u_b(x, t; x_s)为光滑速度模型下的背景波场, 在下文也称为入射波场u_i(x, t; x_s), δu(x, tx_s)为速度模型扰动产生的扰动波场。

假定速度扰动体δv(x)的大小或其速度空间变化特征尺度相对于地震波长比较小, 即满足上文定义的散射体条件, 速度扰动体δv(x)可视为散射体。根据通用的一次散射波方程(12)和方程(13), 可得到下述标量散射波方程:

$ \left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\nabla^{2}\right] u_{\mathrm{s}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\frac{a_{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}} $

(21)

$ \begin{aligned} u_{\mathrm{s}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=& \int_{\varOmega} \mathrm{d} \boldsymbol{x} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) * \frac{a_{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} $

(22)

式中: u_s(x, tx_s)表示由a_v(x)产生的一次散射波; a_v(x)为速度相对扰动, 有a_v(x)=1－v_b²(x)/v²(x)≈2δv(x)/v_b(x); g(x_g, x, t)为下述背景速度模型下标量波方程(23)对应的Green函数; 对于光滑速度模型下的入射波场u_i(x, tx_s), 有:

$ \left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\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)

方程(22)也称为基于速度相对扰动的标量波地震散射数据的线性正演表达。

假定速度扰动体δv(x)的大小或其速度空间变化特征尺度相对于地震波长比较大, 即满足前文定义的反射体条件, 速度扰动体δv(x)可视为反射体。根据通用的一次反射波方程(14), 可得到下述标量反射波方程:

$ \left[\frac{1}{v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\nabla^{2}\right] u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\frac{a_{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}} $

(24)

式中: u_r(x, tx_s)表示由速度相对扰动a_v(x)产生的一次反射波。

根据方程(15), 由方程(24)可得到下述基于速度相对扰动的积分形式标量波地震反射波数据线性正演表达, 即:

$ \begin{aligned} u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=& \int_{\varOmega} \mathrm{d} \boldsymbol{x} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) *_{t} \frac{a_{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} $

(25)

由于标量波动算子中不含有对速度参数的空间导数, 其对应的反射算子也就不含有角度σ, 因此反射数据的表达式(25)与散射数据的表达式(22)有相同的形式, 但其中a_v(x)的分布范围Ω大小不同。

如果反射体δv(x)内部的速度变化很小或不变化, 则反射可视为由反射体边界产生。根据反射体边界局部切平面的局部反射系数定义式(17), 对于速度相对扰动, 有下述的局部反射系数表示式:

$ r(\boldsymbol{x})=\frac{i \omega}{v_{\mathrm{b}}(\boldsymbol{x})} a_{v}(\boldsymbol{x}) $

(26)

利用公式(26)和基于局部反射系数的通用一次反射波方程(18), 可得到基于局部反射系数的标量反射波方程:

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

(27)

根据方程(19), 由方程(27)可以得到基于局部反射系数r的积分形式标量反射波数据线性正演表达式:

$ \begin{gathered} u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=-\int_{\varOmega} \mathrm{d} \boldsymbol{x} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) *_{t} \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} $

(28)

基于Kirchhoff近似, BLEISTEIN>等^[36]推出了与公式(28)相类似的基于公式(6)中的镜面反射系数R的反射数据表达式, 由于公式(6)中镜面反射系数R的严格条件, 我们认为BLEISTEIN等推导出的反射数据表达式对于非平面入射波和反射面为非无穷大平边界存在局限性, 不适合用于构建实际地震数据的偏移成像。

3.2 声波地震数据的线性正演表达

考虑到地下密度变化对地震波传播的影响, 假定通用波动方程(4)中模型m为包含速度v(x)和密度ρ(x)的声学介质, 则方程(4)为非齐次声波方程, 即:

$ \begin{gathered} \left\{\frac{1}{\rho(\boldsymbol{x}) v^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\nabla \cdot\left[\frac{1}{\rho(\boldsymbol{x})} \nabla\right]\right\} u\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) \\ =f(t) \delta\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{s}}\right) \end{gathered} $

(29)

式中: $\boldsymbol{\nabla} $表示梯度算子, 有$\boldsymbol{\nabla} = \frac{\partial }{{\partial x}}\mathit{\boldsymbol{i}} + \frac{\partial }{{\partial y}}\mathit{\boldsymbol{j}} + \frac{\partial }{{\partial z}}\mathit{\boldsymbol{l}};\boldsymbol{\nabla \cdot} $表示求散度运算。利用扰动法, 有v(x)=v_b(x)+δv(x), ρ(x)=ρ_b(x)+δρ(x), u(x, tx_s)=u_b(x, tx_s)+δu(x, tx_s), 其中, v_b(x)和ρ_b(x)为光滑背景模型, δv(x)和δρ(x)分别为速度和密度扰动, u_b(x, tx_s)为光滑模型下的背景波场, 在下文也同样称为入射波场u_i(x, tx_s), δu(x, tx_s)为速度和密度扰动产生的扰动波场。

假定δv(x)和δρ(x)的大小或其空间变化特征尺度相对于地震波长满足散射体条件, 即δv(x)和δρ(x)对应的非均匀体可视为产生散射波的散射体。根据通用的一次散射波方程(12)和方程(13), 可得到下述声散射波方程:

$ \begin{aligned} &\left\{\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\nabla \cdot\left[\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \nabla\right]\right\} u_{\mathrm{s}}(\boldsymbol{x}, t ; \\ &\left.\boldsymbol{x}_{\mathrm{s}}\right)=\frac{\left(a_{v}+a_{\rho}\right)}{\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}}-\nabla \cdot \\ &{\left[\frac{a_{\rho}}{\rho_{\mathrm{b}}(\boldsymbol{x})} \nabla u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)\right]} \end{aligned} $

(30)

$ \begin{aligned} &u_{\mathrm{s}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\int_{\varOmega} \mathrm{d} \boldsymbol{x}\left\{\frac{a_{v}+a_{\rho}}{\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}}-\nabla \cdot\right. \\ &\left.\left[\frac{a_{\rho}}{\rho_{\mathrm{b}}(\boldsymbol{x})} \nabla u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)\right]\right\} *_{t} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) \end{aligned} $

(31)

式中: u_s(x, tx_s)为速度和密度扰动δv(x)和δρ(x)产生的一次散射波场; a_v为速度相对扰动, 有a_v=1－v_b²(x)/v²(x)≈2δv(x)/v_b(x); a_ρ为密度相对扰动, 有a_ρ=1－ρ_b(x)/ρ(x)≈δρ(x)/ρ_b(x); g(x_g, x, t)为Green函数; 对于光滑模型v_b(x)和ρ_b(x)的入射波场u_i(x, tx_s), 有:

$ \begin{gathered} \left\{\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\nabla \cdot\left[\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \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} $

(32)

对方程(31)应用分步积分法, 有:

$ \begin{aligned} &u_{\mathrm{s}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\int_{\varOmega} \mathrm{d} \boldsymbol{x}\left\{\frac{a_{v}+a_{\rho}}{\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}} *_{t}\right. \\ &\left.g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right)+\frac{a_{\rho}}{\rho_{\mathrm{b}}(\boldsymbol{x})} \nabla u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right) *_{t} \nabla g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right)\right\} \end{aligned} $

(33)

方程(33)也称为基于速度、密度相对扰动的声散射波数据的线性正演表达。

假定δv(x)和δρ(x)的大小或其空间变化特征尺度相对于地震波长满足反射体条件, 即δv(x)和δρ(x)对应的非均匀体可视为产生反射波的反射体。根据方程(15)和方程(33), 可得到下述基于速度和密度相对扰动的声反射波方程, 即:

$ \begin{aligned} &u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\int_{\varOmega} \mathrm{d} \boldsymbol{x}\left[\frac{a_{v}+a_{\rho}}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \cdot\right. \\ &\ \ \ \ \ \ \ \ \frac{\partial^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} *_{t} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right)+\frac{a_{\rho}}{\rho_{\mathrm{b}}(\boldsymbol{x})} \nabla u_{\mathrm{i}}(\boldsymbol{x}, t ; \\ &\ \ \ \ \ \ \ \ \left.\left.\boldsymbol{x}_{\mathrm{s}}\right) *_{t} \nabla g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right)\right] \end{aligned} $

(34)

式中: u_r(x_g, t; x_s)为由a_v和a_ρ产生的一次反射波数据。

将方程(34)变换到频率域, 有:

$ \begin{aligned} &u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, \omega ; \boldsymbol{x}_{\mathrm{s}}\right)=\int_{\varOmega} \mathrm{d} \boldsymbol{x}\left[\frac{-\omega^{2}\left(a_{v}+a_{\rho}\right)}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} u_{\mathrm{i}}\left(\boldsymbol{x}, \omega ; \boldsymbol{x}_{\mathrm{s}}\right)\cdot\right. \\ &\left.G\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, \omega\right)+\frac{a_{\rho}}{\rho_{\mathrm{b}}(\boldsymbol{x})} \nabla u_{\mathrm{i}}\left(\boldsymbol{x}, \omega ; \boldsymbol{x}_{\mathrm{s}}\right) \nabla G\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, \omega\right)\right] \end{aligned} $

(35)

光滑速度模型v_b(x)的空间变化相对于地震波长尺度可视为缓慢变化甚至常数, 因此在光滑模型中, 入射波传播方向的局部波数矢量k_i和反射波传播方向的局部波数矢量k_r分别为k_i=ω/v_b(x)$\mathit{\boldsymbol{\hat k}}$_i和k_r=ω/v_b(x)$\mathit{\boldsymbol{\hat k}}$_r, 其中, $\mathit{\boldsymbol{\hat k}}$_i和$\mathit{\boldsymbol{\hat k}}$_r分别为入射波和反射波的传播方向单位矢量。对于$\mathit{\boldsymbol{\hat k}}$_i和$\mathit{\boldsymbol{\hat k}}$_r, 在局部反射面上有$\mathit{\boldsymbol{\hat k}}$_i·$\mathit{\boldsymbol{\hat k}}$_r=-cosσ。(35)式右端第2项中的$\boldsymbol{\nabla} {u_i}\left( {\mathit{\boldsymbol{x}}, \omega ;{\mathit{\boldsymbol{x}}_s}} \right)$和$\boldsymbol{\nabla} G\left( {{\mathit{\boldsymbol{x}}_g}, \mathit{\boldsymbol{x}}, \omega } \right)$分别表示对u_i(x, ω; x_s)和G(x_g, x, ω)进行梯度运算。u_i(x, ω; x_s)和G(x_g, x, ω)在高频近似条件下的波数域表达式分别为ik_iu_i(x, ω; x_s)=－iω/v_b(x)u_i(x, ω; x_s)$\mathit{\boldsymbol{\hat k}}$_i和ik_rG(x_g, x, ω)=iω/v_b(x)G(x_g, x, ω)$\mathit{\boldsymbol{\hat k}}$_r。利用上述高频近似关系, 反射波方程(35)可改写为:

$ \begin{gathered} u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, \omega ; \boldsymbol{x}_{\mathrm{s}}\right)=-\int_{\varOmega} \mathrm{d} \boldsymbol{x} G\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, \omega\right)\cdot \\ \frac{a_{v}+a_{\rho}(1+\cos \sigma)}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \omega^{2} u_{\mathrm{i}}\left(\boldsymbol{x}, \omega ; \boldsymbol{x}_{\mathrm{s}}\right) \end{gathered} $

(36)

将方程(36)变换到时间域, 有:

$ \begin{gathered} u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\int_{\varOmega} \mathrm{d} \boldsymbol{x} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) *_{t} \\ \frac{a_{v}+a_{\rho}(1+\cos \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} $

(37)

方程(37)对应的微分形式方程为:

$ \begin{aligned} &\left\{\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\nabla \cdot\left[\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \nabla\right]\right\} \cdot \\ &u_{\mathrm{r}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\frac{\left[a_{v}+a_{\rho}(1+\cos \sigma)\right]}{\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{aligned} $

(38)

公式(38)中的a_v+a_ρ(1+cosσ)是一个与开角σ有关的声波阻抗相对扰动。令I_r(x, σ)=a_v+a_ρ(1+cosσ), 如果取σ=0, 即沿反射面法向入射, 则有I_r(x, σ=0)=a_v+2a_ρ=2δv(x)/v_b(x)+2δρ(x)/ρ_b(x)≈2δI(x)/I_b(x), I_b(x)=v_b(x)ρ_b(x), δI(x)=δ[v(x)ρ(x)]=ρ(x)δv(x)+v(x)δρ(x)。由此可知, 小角度反射主要受阻抗变化的控制。如果取σ=π, 则有I_r(x, σ=π)=a_v=2δv(x)/v_b(x)。由此可知, 大角度反射主要受速度变化的控制。由方程(37), 我们可以得到基于反射体声波阻抗相对扰动的声波反射数据线性正演表达公式, 即

$ \begin{gathered} u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\int_{\varOmega} \mathrm{d} \boldsymbol{x} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) *_{t} \\ \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} $

(39)

为了得到基于反射体边界局部反射系数的反射波动方程和反射数据线性正演表达公式, 我们定义声反射的局部反射系数为声波阻抗相对扰动I_r(x, σ)沿入射波传播方向I的方向导数。根据反射体边界局部切平面的局部反射系数定义式(17), 对于声波阻抗相对扰动, 有下述的局部反射系数表达式:

$ \begin{aligned} r(\boldsymbol{x}, \sigma)=\mathrm{i}& \boldsymbol{k}_{\mathrm{i}}\left[I_{\mathrm{r}}(\boldsymbol{x}, \sigma)\right]=\frac{\mathrm{i} \omega}{v_{\mathrm{b}}(\boldsymbol{x})}\left[a_{v}+a_{\rho}(1+\right. \\ &\cos \sigma)] \end{aligned} $

(40)

利用(40)式定义的局部反射系数, 由方程(38)在高频近似条件下可得到基于反射体边界局部反射系数的声反射波方程:

$ \begin{gathered} \left\{\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x}) v_{\mathrm{b}}^{2}(\boldsymbol{x})} \frac{\partial^{2}}{\partial t^{2}}-\nabla \cdot\left[\frac{1}{\rho_{\mathrm{b}}(\boldsymbol{x})} \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} $

(41)

根据方程(19), 由方程(41)可以得到基于反射体边界局部反射系数的积分形式声波反射数据线性正演表达公式, 即

$ \begin{gathered} u_{\mathrm{r}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=-\int_{\varOmega} \mathrm{d} \boldsymbol{x} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) *_{t} \\ \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} $

(42)

如果在声波方程中不考虑密度变化, 则声波方程(29)退化为标量波方程(20), 相应的方程(33)、方程(39)和方程(42)分别退化为方程(22)、方程(25)和方程(28)。

3.3 弹性波地震数据的线性正演表达

为了更真实地描述地下地震波的运动, 假定方程(4)中地下模型为各向同性弹性介质, 其模型m的物性参数为Lamé系数λ(x)、μ(x)和密度ρ(x), 则方程(4)为非齐次弹性波动方程, 即:

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

(43)

式中: u(x, tx_s)为位移矢量波场, 在笛卡尔坐标系下, 有u(x, tx_s)=[u^x(x, tx_s), u^y(x, tx_s), u^z(x, tx_s)], 其中, u^x(x, tx_s), u^y(x, tx_s), u^z(x, tx_s)分别代表u(x, tx_s)的x, y, z方向分量; f(t, x_s)为矢量震源函数; L是一个3×3偏微分算子:

$ \boldsymbol{L}=\left(\begin{array}{lll} L_{x x} & L_{x y} & L_{x z} \\ L_{y x} & L_{y y} & L_{y z} \\ L_{z x} & L_{z y} & L_{z z} \end{array}\right) $

(44)

其中, ${L_{ii}} = \rho \left( \mathit{\boldsymbol{x}} \right)\frac{{{\partial ^2}}}{{\partial {t^2}}} - \frac{\partial }{{\partial i}}\left[ {\lambda \left( \mathit{\boldsymbol{x}} \right) + 2\mu \left( \mathit{\boldsymbol{x}} \right)} \right]\frac{\partial }{{\partial i}} - \sum\limits_{j \ne i} {\frac{\partial }{{\partial j}}\mu \left( \mathit{\boldsymbol{x}} \right)} \frac{\partial }{{\partial j}}\mu \left( \mathit{\boldsymbol{x}} \right)\frac{\partial }{{\partial j}}, {L_{ij}} = \frac{\partial }{{\partial i}}\lambda \left( \mathit{\boldsymbol{x}} \right)\frac{\partial }{{\partial j}} + \frac{\partial }{{\partial j}}\mu \left( \mathit{\boldsymbol{x}} \right)\frac{\partial }{{\partial i}}, j \ne i, i, j = x, y, z$。

将弹性波动方程(43)中的模型参数转化为纵波速度[α(x)]、横波速度[β(x)]和密度[ρ(x)], 有α²(x)=[λ(x)+2μ(x)]/ρ(x)、β²(x)=μ(x)/ρ(x), 相应的纵、横波速度和密度光滑模型[α_b(x)、β_b(x)、ρ_b(x)]满足α_b²(x)=[λ_b(x)+2μ_b(x)]/ρ_b(x)、β_b²(x)=μ_b(x)/ρ_b(x), 对应的纵、横波速度和密度相对扰动为a_α=1－α_b²(x)/α²(x)≈2δα(x)/α_b(x), a_β=1－β_b²(x)/β²(x)≈2δβ(x)/β_b(x), a_ρ=1－ρ_b(x)/ρ(x)≈δρ(x)/ρ_b(x)。地震波在光滑模型中传播时不发生散射、反射和波型转换。

假定扰动量δα(x)、δβ(x)、δρ(x)的大小或其空间变化特征尺度相对于地震波长满足前文定义的散射体条件, 即δα(x)、δβ(x)、δρ(x)对应的非均匀体可视为产生散射波的散射体。根据通用散射波方程(12), 可得到下述通用的弹性散射波方程:

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

(45)

式中: u_s(x, tx_s)为弹性一次散射波的位移矢量波场; L_b为光滑模型下的弹性波传播算子; ΔL_s为弹性波散射算子, 有:

$ \Delta \boldsymbol{L}_{\mathrm{s}}=\boldsymbol{L}-\boldsymbol{L}_{\mathrm{b}}=\left(\begin{array}{ccc} \Delta L_{\mathrm{s}}^{x x} & \Delta L_{\mathrm{s}}^{x y} & \Delta L_{\mathrm{s}}^{x z} \\ \Delta L_{\mathrm{s}}^{y x} & \Delta L_{\mathrm{s}}^{y y} & \Delta L_{\mathrm{s}}^{y z} \\ \Delta L_{\mathrm{s}}^{z x} & \Delta L_{\mathrm{s}}^{z y} & \Delta L_{\mathrm{s}}^{z z} \end{array}\right) $

(46)

其中, $\Delta L_s^{ii} = \rho_{\mathrm{b}} \left[ { - {a_\rho }\frac{{{\partial ^2}}}{{\partial {t^2}}} + \alpha _{\rm{b}}^{\rm{2}}\frac{\partial }{{\partial i}}\left( {{a_\rho } + {a_\alpha }} \right)\frac{\partial }{{\partial i}} + \beta _b^2 \cdot \sum\limits_{j \ne i} {\frac{\partial }{{\partial j}}\left( {{a_\rho } + {a_\beta }} \right)\frac{\partial }{{\partial j}}} } \right], \Delta L_s^{ij} = {\rho _{\rm{b}}}\left[ {\alpha _{\rm{b}}^2\frac{\partial }{{\partial i}}\left( {{a_\rho } + {a_\alpha }} \right) \cdot \frac{\partial }{{\partial j}} - 2\beta _{\rm{b}}^2\frac{\partial }{{\partial i}}\left( {{a_\rho } + {a_\beta }} \right)\frac{\partial }{{\partial j}} + \beta _{\rm{b}}^2\frac{\partial }{{\partial j}}\left( {{a_\rho } + {a_\beta }} \right)\frac{\partial }{{\partial i}}} \right], j \ne i, {\mathit{\boldsymbol{u}}_{\rm{i}}}\left( {\mathit{\boldsymbol{x}}, t;{\mathit{\boldsymbol{x}}_{\rm{s}}}} \right)$为光滑模型下入射波的位移矢量波场, 有:

$ \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) $

(47)

利用光滑模型弹性波方程(47)所对应的矢量Green函数, 可将方程(45)中的弹性散射波数据写成与方程(13)对应的积分形式^[65]。方程(45)也称为基于纵、横波速度、密度相对扰动的弹性散射波数据的线性正演方程。

假定δα(x), δβ(x), δρ(x)的大小或其空间变化特征尺度相对于地震波长满足前文定义的反射体条件, 即δα(x), δβ(x), δρ(x)对应的非均匀体可视为产生反射波的反射体。根据通用反射波方程(14), 利用与前两节类似的推导可得到弹性反射波方程, 即:

$ \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) $

(48)

式中: u_r(x, tx_s)为弹性一次反射波的位移矢量波场; ΔL_r为弹性波反射算子。

由于笛卡尔坐标系下的弹性波场中P波和S波耦合, 在弹性波反射数据的逆时偏移成像中, 为了避免弹性波场分量成像时产生交叉噪声和得到易于物理解释的成像结果, 需要在成像前将弹性波场解耦为P/S波型波场^[66-69]。本文将方程(48)中笛卡尔坐标系下的入射矢量波场u_i(x, tx_s)=(u_i^x(x, tx_s), u_i^y(x, tx_s), u_i^z(x, tx_s))和反射矢量波场u_r(x, tx_s)=[u_r^x(x, tx_s), u_r^y(x, tx_s), u_r^z(x, tx_s)]解耦为以P波、SH波和SV波的极化方向为坐标系的波场u_i^P(x, tx_s)和u_r^PS(x, tx_s), 有u_i^PS(x, tx_s)=[u_i^P(x, tx_s), u_i^SH(x, tx_s), u_i^SV(x, tx_s)], u_r^PS(x, tx_s)=[u_r^P(x, tx_s), u_r^SH(x, tx_s), u_r^SV(x, tx_s)](详细的弹性波场P/S解耦分离见本研究工作的第二部分)。再将光滑模型下的弹性波传播算子L_b和反射算子ΔL_r也变换成P/S波型对应的算子, 则方程(48)可表示为:

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

(49)

式中: ΔL_r^PS为与ΔL_r对应的反射算子^[70], 有:

$ \Delta \boldsymbol{L}_{\mathrm{r}}^{\mathrm{PS}}=\left(\begin{array}{lll} \Delta L_{\mathrm{r}}^{\mathrm{PP}} & \Delta L_{\mathrm{r}}^{\mathrm{PSH}} & \Delta L_{\mathrm{r}}^{\mathrm{PSV}} \\ \Delta L_{\mathrm{r}}^{\mathrm{SHP}} & \Delta L_{\mathrm{r}}^{\mathrm{SHSH}} & \Delta L_{\mathrm{r}}^{\mathrm{SHSV}} \\ \Delta L_{\mathrm{r}}^{\mathrm{SVP}} & \Delta L_{\mathrm{r}}^{\mathrm{sVSH}} & \Delta L_{\mathrm{r}}^{\mathrm{SVSV}} \end{array}\right) $

(50)

ΔL_r^PS中的各个元素代表不同类型入射波与不同类型反射波(或转换波)间的反射(或转换)作用。在地震波长相对于背景模型的空间变化特征变化尺度满足高频近似的条件下, ΔL_r^PS中的元素在频率域满足^{[56, 70]}:

$ \begin{gathered} \Delta L_{\mathrm{r}}^{\mathrm{PP}}=-\rho_{\mathrm{b}} \omega^{2}\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} $

(51)

$ \Delta L_{\mathrm{r}}^{\mathrm{SHSH}}=-\rho_{\mathrm{b}} \omega^{2}\left[a_{\beta} \cos \sigma+a_{\rho}(1+\cos \sigma)\right] $

(52)

$ \Delta L_{\mathrm{r}}^{\mathrm{SVSV}}=-\rho_{\mathrm{b}} \omega^{2}\left[a_{\beta} \cos 2 \sigma+a_{\rho}(\cos \sigma+\cos 2 \sigma)\right] $

(53)

$ \Delta L_{\mathrm{r}}^{\mathrm{SVP}}=-\rho_{\mathrm{b}} \omega^{2} \sin \sigma\left[a_{\rho}\left(\frac{\alpha_{\mathrm{b}}}{\beta_{\mathrm{b}}}+2 \cos \sigma\right)+2 a_{\beta} \cos \sigma\right] $

(54)

$ \Delta L_{\mathrm{r}}^{\mathrm{PSV}}=\rho_{\mathrm{b}} \omega^{2} \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] $

(55)

$ \Delta L_{\mathrm{r}}^{\mathrm{SHP}}=0, \Delta L_{\mathrm{r}}^{\mathrm{PSH}}=0, \Delta L_{\mathrm{r}}^{\mathrm{SHSV}}=0, \Delta L_{\mathrm{r}}^{\mathrm{SVSH}}=0 $

在法向入射时, σ=0, 有:

$ \Delta L_{\mathrm{r}}^{\mathrm{PP}}=-\rho_{\mathrm{b}} \omega^{2}\left(a_{\alpha}+2 a_{\rho}\right) $

(56)

$ \Delta L_{\mathrm{r}}^{\mathrm{SHSH}}=-\rho_{\mathrm{b}} \omega^{2}\left(a_{\beta}+2 a_{\rho}\right) $

(57)

$ \Delta L_{\mathrm{r}}^{\mathrm{SVSV}}=-\rho_{\mathrm{b}} \omega^{2}\left(a_{\beta}+2 a_{\rho}\right) $

(58)

$ \Delta L_{\mathrm{r}}^{\mathrm{SHP}}=0, \Delta L_{\mathrm{r}}^{\mathrm{PSH}}=0, \Delta L_{\mathrm{r}}^{\mathrm{SHSV}}=0, \Delta L_{\mathrm{r}}^{\mathrm{SVSH}}=0, $

$ \Delta L_{\mathrm{r}}^{\mathrm{SVP}}=0, \Delta L_{\mathrm{r}}^{\mathrm{PSV}}=0 $

此时不发生波型转换。

将ΔL_r^PS中各个元素代入(50)式, 然后再将(50)式转换到时间域, 有:

$ \Delta \boldsymbol{L}_{\mathrm{r}}^{\mathrm{PS}}=\rho \boldsymbol{E}_{\mathrm{r}}^{\mathrm{IPS}} \frac{\partial^{2}}{\partial t^{2}} $

(59)

其中, ${\bf{E}}_{\rm{r}}^{{\rm{IPS}}} = \left( {\begin{array}{*{20}{c}} {I_{\rm{r}}^{{\rm{pp}}}}&0&{I_{\rm{r}}^{{\rm{PSV}}}}\\ 0&{I_{\rm{r}}^{{\rm{SHSH}}}}&0\\ {I_{\rm{r}}^{{\rm{SVP}}}}&0&{I_{\rm{r}}^{{\rm{SVSV}}}} \end{array}} \right)$。

弹性阻抗相对扰动矩阵E_r^IPS中的各个元素称为弹性波阻抗相对扰动, 有:

$ 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 $

(60)

$ I_{\mathrm{r}}^{\mathrm{SHSH}}=a_{\beta} \cos \sigma+a_{\rho}(1+\cos \sigma) $

(61)

$ I_{\mathrm{r}}^{\mathrm{SVSV}}=a_{\beta} \cos 2 \sigma+a_{\rho}(\cos \sigma+\cos 2 \sigma) $

(62)

$ 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] $

(63)

$ \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} $

(64)

如果假定横波速度为0, 即β_b=0, 则上面弹性PP波阻抗相对扰动I_r^PP退化为声波方程(39)中的声波阻抗相对扰动I_r=a_α+a_ρ(1+cosσ)。

将(59)式代入方程(49), 有:

$ \boldsymbol{L}_{\mathrm{b}}^{\mathrm{PS}} \boldsymbol{u}_{\mathrm{r}}^{\mathrm{PS}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=\rho \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}} $

(65)

方程(65)称为高频近似条件下基于弹性波阻抗相对扰动矩阵E_r^IPS的P/S波型反射波方程。根据方程(15), 由方程(65)可以得到基于弹性波阻抗相对扰动矩阵E_r^IPS的P/S波型积分形式的弹性反射波数据线性正演表达公式, 即:

$ \begin{aligned} \boldsymbol{u}_{\mathrm{r}}^{\mathrm{PS}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=& \int_{\varOmega} \mathrm{d} \boldsymbol{x} \boldsymbol{g}\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) *_{t} \rho_{\mathrm{b}}(\boldsymbol{x}) \boldsymbol{E}_{\mathrm{r}}^{\mathrm{IPS}}\cdot \\ & \frac{\partial^{2} \boldsymbol{u}_{\mathrm{i}}^{\mathrm{PS}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t^{2}} \end{aligned} $

(66)

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

光滑背景模型的空间变化相对于地震波长尺度可视为缓慢变化甚至常数, 对于P波入射, 其入射波方向的局部波数k_i^P=ω/α_b(x); 对于S波入射, 其入射波方向的局部波数k_i^S=ω/β_b(x)。为了得到基于弹性局部反射系数的反射波方程, 我们定义反射体边界弹性局部反射系数为反射体弹性阻抗相对扰动沿入射波传播方向I的方向导数。根据反射体边界局部切平面的局部反射系数定义式(17), 结合(50)式中各元素的弹性波阻抗相对扰动表达式, 对于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)-\right. \\ \left.2 \frac{\beta_{\mathrm{b}}^{2}}{\alpha_{\mathrm{b}}^{2}}\left(a_{\rho}+a_{\beta}\right) \sin ^{2} \sigma\right] \end{gathered} $

(67)

对于SH波入射SH波反射的局部反射系数, 有:

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

(68)

对于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] $

(69)

对于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] $

(70)

对于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} $

(71)

将上述弹性反射波的局部反射系数定义式代入(50)式的弹性波反射算子ΔL_r^PS, 得到基于局部反射系数的弹性波反射算子ΔL_r^RPS, 即:

$ \begin{aligned} \Delta \boldsymbol{L}_{\mathrm{r}}^{\mathrm{RPS}} &=\\ &\left(\begin{array}{ccc} \rho_{\mathrm{b}} \alpha_{\mathrm{b}} r^{\mathrm{PP}} \mathrm{i} \omega & 0 & \rho_{\mathrm{b}} \beta_{\mathrm{b}} r^{\mathrm{PSV}} \mathrm{i} \omega \\ 0 & \rho_{\mathrm{b}} \beta_{\mathrm{b}} r^{\mathrm{SHSH}} \mathrm{i} \omega & 0 \\ \rho_{\mathrm{b}} \alpha_{\mathrm{b}} r^{\mathrm{SVP}} \mathrm{i} \omega & 0 & \rho_{\mathrm{b}} \beta_{\mathrm{b}} r^{\mathrm{SVSV}} \mathrm{i} \omega \end{array}\right) \\ &=i \omega \rho_{\mathrm{b}} \boldsymbol{E}_{\mathrm{r}}^{\mathrm{RPS}} \end{aligned} $

(72)

$ \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) $

(73)

式中: E_r^RPS称为弹性局部反射系数矩阵。

由方程(65)可得到基于弹性局部反射系数矩阵E_r^RPS的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 u_{\mathrm{i}}^{\mathrm{PS}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} $

(74)

根据通用方程(19), 由方程(74)可以得到基于弹性局部反射系数矩阵E_r^RPS的P/S波型积分形式的弹性反射波数据线性正演表达公式, 即:

$ \begin{gathered} \boldsymbol{u}_{\mathrm{r}}^{\mathrm{PS}}\left(\boldsymbol{x}_{\mathrm{g}}, t ; \boldsymbol{x}_{\mathrm{s}}\right)=-\int_{\varOmega} \mathrm{d} \boldsymbol{x} g\left(\boldsymbol{x}_{\mathrm{g}}, \boldsymbol{x}, t\right) *_{t} \rho_{\mathrm{b}}(\boldsymbol{x}) \boldsymbol{E}_{\mathrm{r}}^{\mathrm{RPS}} \cdot\\ \frac{\partial u_{\mathrm{i}}^{\mathrm{PS}}\left(\boldsymbol{x}, t ; \boldsymbol{x}_{\mathrm{s}}\right)}{\partial t} \end{gathered} $

(75)

由上述推导出的标量波、声波和弹性波的反射数据线性正演表达公式可知, 基于波阻抗相对扰动的反射数据线性正演表达公式的右端包含入射波场的时间二阶导数与波阻抗相对扰动的乘积, 而基于局部反射系数的反射数据线性正演表达公式的右端包含入射波场的时间一阶导数与局部反射系数的乘积。不同于平面波无穷大平边界的无量纲镜面反射系数, 本文基于高频近似和局部平面波定义的局部切平面反射系数是一个具有长度倒数量纲的变量, 是真实模型与光滑模型之间模型参数相对扰动沿入射波传播方向的方向导数, 它不仅依赖于光滑模型, 也与入射波传播方向有关。Zoeppritz方程不仅描述了反射波与入射波之间的关系, 还定义了平面波无穷大平边界镜面反射系数与入射方向和地下模型物理参数之间的定量关系。本文利用反射体近似和高频近似推导的基于反射体弹性波阻抗相对扰动的P/S波型弹性反射波方程与Zoeppritz方程类似, 不仅建立了反射波与入射波之间的数学物理关系, 也定义了反射体弹性波阻抗相对扰动与入射方向和地下模型物理参数之间的定量关系, 可为进一步的基于偏移成像得到的角度域共成像点道集的岩性分析与反演提供理论基础。由于本文定义的局部反射系数与入射波传播方向的局部波数有关, 因此在目前我们还不能给出倾斜入射情况下局部反射系数与地下模型物理参数和入射波传播方向之间的具体函数关系。

4 结论

当前地震数据波动方程偏移成像方法不区分散射波和反射波, 也缺乏与之相应的严谨正演方程, 致使偏移成像不能准确地利用地震数据的波形信息, 得到的偏移成像结果存在相位误差。为了使波动方程偏移成像方法能更好地适应地下实际情况和准确地利用波形信息, 弥补当前偏移成像方法理论的不足, 修正地震数据波动方程偏移成像结果的相位误差, 并为发展地层物性参数和储层参数成像方法打下基础, 本文在给定满足地震波运动学的准确的光滑模型基础上, 利用扰动方法将波动方程转化为描述非均匀体产生的扰动波场的波动方程, 根据地下非均匀体的大小或非均匀体物理性质空间变化特征尺度与地震波长之间的相对大小关系, 提出将非均匀体划分为相对于地震波长在空间上具有有限延续度的产生散射波的散射体或相对于地震波长在空间上具有一定延续度的产生反射波的反射体, 推导出反射波波动方程。为区别于散射波, 本文认为反射波是由点源集成的局部平面波源产生的具有方向性的局部平面波。对于散射波, 利用小扰动假定和一阶Born近似(忽略多次散射波), 可得到基于散射体物性参数相对扰动的地震一次散射数据的线性正演表达。对于反射波, 利用波阻抗相对扰动的小值假定、一次反射波近似(忽略多次反射波)和高频近似, 可得到基于反射体波阻抗相对扰动的地震一次反射数据的线性正演表达。反射体的波阻抗相对扰动不同于散射体的物性参数相对扰动, 它不仅与反射体的物性参数相对扰动有关, 还与反射开角有关。基于反射体弹性波阻抗相对扰动的P/S波型弹性反射波方程, 不仅建立了反射波与入射波之间的数学物理关系, 也定义了反射体弹性波阻抗相对扰动与入射方向和地下模型物理参数相对扰动之间的定量关系。对于高频近似条件下定义的反射体与反射波, 本文将波阻抗相对扰动沿入射波传播方向的方向导数定义为反射体边界局部切平面对于局部平面波的局部反射系数, 得到基于局部反射系数的地震一次反射数据的线性正演表达。不同于平面波入射到无穷大平界面的与频率无关的无量纲镜面反射系数, 本文定义的局部反射系数具有长度倒数的量纲且与入射波传播方向有关(即入射波的波数), 推导出的反射数据表达式具有更广的适应性。本文推导出的地震数据线性正演表达(特别是基于反射体波阻抗相对扰动的地震反射数据线性正演表达), 是开展地震数据波形成像和其它波形线性反演方法研究的基础方程。