基于弹性波场矢量分解与低秩分解的地震正演模拟方法

摘要

本发明属于勘探地球物理学领域，具体地，涉及一种弹性地震波场正演模拟方法。首先，该方法基于波场矢量分解原理将弹性地震波场分解为纵波波场和横波波场，并分别计算相应的纵、横波波场传播矩阵，所得到的传播矩阵包含了对介质参数和数值模拟参数的补偿，因而具有较高的计算精度；然后，基于低秩分解原理将传播矩阵进行分解得到矢量纵波和矢量横波的传播算子，从而降低数值模拟的计算量，提高计算效率；最后，将数值模拟得到的更新后的矢量纵波记录和矢量横波记录进行耦合，得到地震波场正演模拟结果。本发明构建的正演方法得到的地震波场记录几乎不存在数值频散，数值模拟稳定性较高，是一种高精度的弹性波场数值模拟方法。

说明书

技术领域

本发明属于勘探地球物理学领域，具体地，涉及一种弹性地震波场正演模 拟方法。

背景技术

正演模拟方法是研究地震波在地球介质中传播规律的重要手段，同时也是 开展地震波成像方法和波动方程反演方法等技术研究的理论基础。就描述地球 介质的方程而言，目前最为常用的方法主要有声波方程正演模拟和弹性波方程 正演模拟。在过去相当长的时间里，基于声波方程的正演方法由于其计算速度 快，计算机内存需求小而受到了工业界和学术界的青睐，发挥了巨大的作用。 但是实际地球介质是一个弹性体，简单的声学近似不足以描述地震波在其中传 播的全部波场信息。特别是随着近年来石油工业的发展，涌现出大量的以弹性 波为基础的地震勘探方法，例如弹性波逆时偏移成像、多分量联合反演等。新 的发展也对弹性波模拟方法提出了更高的要求，因此开展高精度的弹性波正演 模拟方法研究势在必行。

目前普遍使用的弹性波数值模拟方法主要有有限差分方法、有限元方法和 伪谱法。其中，有限差分方法具有计算速度快、内存需求小的优点，但是这种 算法的数值模拟精度较低，越来越难以满足工业界对高精度数值模拟算法的需 求；有限元方法可以模拟任意几何形状的地球介质模型，数值模拟精度很高， 但是这类方法的计算量非常大，并且数值频散较为严重；谱方法通过快速正反 傅里叶变换计算微分，在均匀介质情况下可以认为是一种具有无穷高阶近似的 数值方法，但是不适用于复杂变速介质限制了该类算法的应用范围。近些年来， 在声波数值模拟中，新发展出了一种矩阵低秩分解的谱方法，这类方法对数值 微分进行了一个与介质参数和模拟参数相关的补偿，极大的提高了计算精度， 且计算量较为适中，但是其传播算子并不适用于弹性波方程。

本专利提出了一种基于弹性矢量波场分解和矩阵低秩分解构建弹性波波场 外推算子，该方法是一种弹性波递归时间积分波场外推算法，保持了传统谱方 法计算精度高，同时克服了传统谱方法不适用于复杂变速介质的难题。该方法 对弹性波场正演模拟及其相关领域的研究具有重要意义。

发明内容

为克服现有技术的缺陷，本发明提供一种弹性地震波场正演模拟方法。

为实现上述目的，本发明采用下述方案：

基于波场矢量分解原理将弹性地震波场分解为一个矢量纵波波场和一个矢 量横波波场，并分别计算相应的纵、横波波场传播矩阵；然后，基于矩阵低秩 分解原理将传播矩阵进行分解得到矢量纵波和矢量横波的传播算子；最后，将 数值模拟得到的更新后的矢量纵波记录和矢量横波记录进行耦合，得到地震波 场正演模拟结果，具体包括如下步骤：

步骤1，基于波场矢量分解原理将弹性地震波场分解为一个矢量纵波波场 和一个矢量横波波场，根据递归时间积分原理分别计算矢量纵波和矢量横波的 的地震波场传播矩阵；

步骤2，基于矩阵低秩分解原理，分别将第1步计算得到的纵横波传播矩阵 进行分解得到相应的地震波场传播算子；

步骤3，基于第1步和第2步得到的地震波场传播算子进行弹性波场外推。 在每一个时间步长中，分别应用纵横波的地震波传播算子计算得到更新后的纵、 横波场；

步骤4，基于第1步、第2步和第3步得到的更新后的纵、横波场，通过纵、 横波场耦合得到更新后的弹性波场。

相对于现有技术，本发明的有益效果如下：

1、基于弹性波场矢量分解构建的弹性波场正演模拟方法本质上是一种谱方 法，因此该算法具有计算精度高，几乎不存在数值频散且算法稳定性较好。同 时在算子推导中对介质参数(弹性参数和密度)和模拟参数(空间步长和时间 步长)进行了补偿，在一定程度上解决了谱方法不适用于复杂变速介质正演的 难题；

2、虽然本文提出的波场传播算子依然是作用于单独的纵波或横波，但是通 过求解弹性波矢量分解而得到波场传播矩阵内隐含了弹性波场的全部信息，所 以本文方法是一种弹性波的正演模拟方法。相比较于传统的伪声波或伪横波方 程递归时间积分波场外推算法，本文提出数值模拟方法能够更加真实的模拟地 震波场的动力学特征和运动学特征；

3、由于传播算子矩阵通常具有较低的秩，所以通过矩阵低秩分解得到的两 个子矩阵的规模都远远小于原始传播矩阵。因此，矩阵低秩分解可以极大的降 低迭代过程中的傅里叶反变换次数，从而提高数值模拟的计算效率。

附图说明

图1是基于波场矢量分解与耦合的弹性波正演模拟方法流程图；

图2是三维空间角度定义示意图；

图3是地质模型示意图；

图4(a)是更新后的纵波波场；

图4(b)是更新后的横波波场；

图5(a)是本专利方法得到的正演模拟结果；

图5(b)是传统有限差分方法得到的正演模拟结果。

图6是Sigsbee2a模型示意图；

图7(a)是1050毫秒时刻本专利所提出的正演模拟方法得到的波场快照；

图7(b)是1050毫秒时刻交错网格有限差分方法得到的波场快照；

图8(a)是本专利所提出的正演模拟方法得到的地震剖面；

图8(b)是交错网格有限差分方法得到的地震剖面。

具体实施方式

实施例一，如图1所示，在对一个实际地质模型(如图-3所示)进行模拟 时，基于波场矢量分解原理将弹性地震波场分解为一个矢量纵波波场和一个矢 量横波波场，并分别计算相应的纵、横波波场传播矩阵；然后，基于矩阵低秩 分解原理将传播矩阵进行分解得到矢量纵波和矢量横波的传播算子；最后，将 数值模拟得到的更新后的矢量纵波记录和矢量横波记录进行耦合，得到地震波 场正演模拟结果，具体包括如下步骤：

步骤1，基于波场矢量分解原理将弹性地震波场分解为一个矢量纵波波场 和一个矢量横波波场，根据递归时间积分原理分别计算矢量纵波和矢量横波的 的地震波场传播矩阵。

波场矢量分解的具体方法如下：

第一步，对于地震波场的位移矢量U进行空间傅里叶变换，得到波数域的 地震波场

$\hat{U}=\mathrm{FFT}\left(U\right)$

其中，U表示空间域的位移矢量、表示波数域的位移矢量，FFT表示快速傅 里叶变换；

第二步，对于波数域的位移矢量点乘纵波偏振矢量A^P得到波数域的纵波 位移矢量

${\hat{U}}^P=A^P(A^P·\hat{U})$

其中，表示波数域的位移矢量、表示波数域的纵波位移矢量，A^P表示纵波 偏振矢量；

第三步，对于波数域的位移矢量叉乘纵波偏振矢量A^P得到波数域的横波 位移矢量

${\hat{U}}^s=-A^P\times (A^P\times \hat{U})$

其中，表示波数域的位移矢量、表示波数域的横波位移矢量，A^P表示纵波 偏振矢量；

第四步，对得到的结果进行傅里叶反变换，从而得到分解后的P波矢量和 横波矢量：

$U^P={\mathrm{FFT}}^{-1}\left[{\hat{U}}^P\right]$

$U^S={\mathrm{FFT}}^{-1}\left[{\hat{U}}^S\right]$

其中，表示波数域的纵波位移矢量、表示波数域的横波位移矢量、表示 空间域的纵波位移矢量、表示空间域的横波位移矢量、FFT^-1表示傅里叶反变 换

递归时间积分的具体方法如下，以纵波为例：

第一步，假设该质点传播的地震波场可以近似成一个平面波：

$u^P=u_0^P·e^{\pm j({\tilde{k}}_{x}x+{\tilde{k}}_{y}y+{\tilde{k}}_{z}z)}·e^{-\mathrm{j\omega t}}$

其中，是一个常数，表示P波振幅；表示沿 三维空间不同方向的归一化波数；ω、t和j分别表示圆周频率，时间和复数单 位，并且ω＝|k|v_p；

第二步，平面P波沿空间水平的分量同样是一个平面波，通过坐标轴映射 可以得到该分量的表达式为：

其中，角度θ和的分别表示平面波传播方向与水平方向和垂直方向的夹角，具 体定义如图2所示；

第三步，根据平面波公式，可以得到P波沿水平方向的分量的递归时间积 分波场外推公式：

$\left(\begin{array}{c}{u}^P(n+1)+u^P(n-1)-2u^P\left(n\right)=2(\cos \left(\omega \frac{\mathrm{\Delta t}}{2}\right)-1)·u^P\left(n\right)\\ =2(\cos \left(\right|k|·v_p·\frac{\mathrm{\Delta t}}{2})-1)·u^P\left(n\right)\end{array}\right)$

同理，横波的波场外推公式表示为：

$\left(\begin{array}{c}{u}^S(n+1)+u^S(n-1)-2u^S\left(n\right)=2(\cos \left(\omega \frac{\mathrm{\Delta t}}{2}\right)-1)·u^S\left(n\right)\\ =2(\cos \left(\right|k|·v_S·\frac{\mathrm{\Delta t}}{2})-1)·u^S\left(n\right)\end{array}\right)$

其中，u^P和u^S分别表示纵横波位移场；n+1、n、n-1分别表示下一个时刻、当 前时刻和过去时刻；v_P和v_S分别表示纵横波的相速度；ω、Δt和|k|分别表示圆 周频率、时间步长和波数

步骤2，通过矩阵低秩分解原理，将递归时间积分波场外推矩阵进行分解， 分别得到纵波和横波的外场外推算子。矩阵低秩分解后的子矩阵由原矩阵的部 分行、部分列以及数个常数组成，可以表示为：

$W(x,k)\approx W^{1}A{W}^2\approx \underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}W(x,k_m)·a_{\mathrm{mn}}·W(x_n,k)=M^L·M^R$

其中，矩阵W(x,k)是一个需要进行分解的目标矩阵；W(x,k_m)表示一个由矩阵W某 几列组成的矩阵；W(x_n,k)表示一个由矩阵W某几行组成的矩阵；a_mn表示一个常 数矩阵；M^L＝W(x,k_m)·a_mn，M^R＝W(x_n,k)；

矩阵低秩分解的具体做法为：

第一步，随机选取W(x,k)中的β·r_ε列(通常β＝3或4即可，r_ε表示矩阵W的 秩)构成一个新的矩阵S。对新构建的矩阵进行主元QR分解，其前r_ε个主元列 对应于S矩阵的r_ε行。定义W₁是原始矩阵W(x,k)的一个子矩阵，包含对应与r_ε的 列数据；

第二步，随机选取W(x,k)中的β·r_ε行(通常β＝3或4，r_ε表示矩阵W的秩) 构成一个新的矩阵T。定义W₂是原始矩阵W(x,k)的一个子矩阵，包含对应与r_ε的 行数据；

第三步，选取矩阵W₁和矩阵W₂相互交叉的部分求伪逆，即为中间矩阵 A＝W⁺(x_n,k_m)；

第四步，结合以上步骤就可以得到矩阵的低秩分解W(x,k)≈W¹AW²。

步骤3，在每一个时刻的波场外推时，采用如下的计算流程

第一步，在每一个时间步长中，对地震波场进行傅里叶正变换得到波数域 的地震波场

$\hat{U}=\mathrm{FFT}\left(U\right)$

其中，U表示空间域的位移矢量、表示波数域的位移矢量，FFT表示快速傅 里叶变换；

第二步，对波数域的地震波场分别乘以矢量纵波波场外推算子和矢量 横波波场外推算子：

其中，表示波数域的位移矢量；表示波数域的纵波位移矢量；表示波数 域的横波位移矢量；n+1、n-1和n表示下一个时间时刻、过去时刻和当前时刻； MP^R和MS^R分别表示纵波和横波是地震波传播右矩阵；

第三步，对波数域的纵横波矢量和进行傅里叶反变换，并乘以地震波 传播左矩阵得到更新后的纵横波波场，如图4(a)和图4(b)所示：

其中，和表示波数域的纵横波位移矢量，U^P和U^S表示空间域的纵横波位 移矢量；MP^L和MS^L分别表示纵波和横波是地震波传播左矩阵；IFFT表示傅里 叶反变换。

步骤4，将更新后的纵横波位移矢量耦合得到地震波正演模拟记录，如图5 (a)所示：

U＝U^P+U^S

其中，U^P和U^S表示空间域的纵横波位移矢量，U表示地震波正演模拟记录。

图5(b)为传统有限差分方法得到的波场快照，通过对比，本专利提出的 弹性波数值模拟方法能准确的描述弹性波场的地震波传播特征，通过该算法得 到的地震记录几乎不存在数据频散，具有较高的计算稳定性和数值模拟精度。

在第二个实施例中，为了验证本专利方法对复杂模型的适用性，将本发明提 出的基于弹性波场矢量分解与低秩分解的地震正演模拟方法用于Sigsbee2a模型 (详见文献：Paffenholz J.,McLain B.,Zaske J.,et al.Subsalt Multiple Attenuation  and Imaging:Observations from the Sigsbee2b Synthetic Dataset[A].in 2002 SEG  Annual Meeting[C].Society of Exploration Geophysicists,2002)，该模型纵波速度 如图6所示。通过对比弹性波正演模拟得到的波场快照(图7)和地震剖面(图 8)：本专利提出的正演模拟方法较为稳定，没有出现数值溢出；得到的地震记 录同相轴清晰，几乎没有数值频散。而传统的交错网格正演模拟方法所得到的 波场快照和地震剖面都存在较为明显的数值频散，影响了正演模拟的精度。通 过Sigsbee2a模型测试，证明了本专利所提出的基于弹性波场矢量分解与低秩分 解的地震正演模拟方法的有效性。