CN104459773A - 基于交错网格Lowrank分解的无条件稳定地震波场延拓方法 - Google Patents

Abstract

本发明涉及一种基于交错网格Lowrank算子的无条件稳定有限差分地震波场延拓方法，包括以下步骤：利用一阶速度-应力方程的Fourier积分解构建交错网格上的Lowrank算子；利用衰减函数对交错网格Lowrank算子进行衰减约束；利用加权最小二乘方法计算有限差分系数；利用得到的有限差分系数实现无条件稳定的有限差分波场延拓。本发明基于Lowrank分解近似波数-空间域算子，通过衰减约束和加权最小二乘方法得到计算稳定且精度高的优化有限差分系数，构建的有限差分计算格式能实现无条件稳定的地震波场延拓。

Description

基于交错网格Lowrank分解的无条件稳定地震波场延拓方法

技术领域

本发明属于勘探地球物理学领域，具体地，涉及一种基于交错网格Lowrank算子的无条件稳定有限差分波地震场延拓方法。

背景技术

石油和天然气是关系到国计民生的重要能源。地震勘探是寻找油气资源，解决油气勘探、开发问题的有效方法。地震正演、地震成像和地震反演是地震勘探中的重要技术，而这三种技术的计算效率和精度都依赖于所采用的时间域地震波场延拓方法。目前最常用的两类时间域波场延拓方法包括有限差分方法和谱方法。有限差分利用差分代替微分，具有计算量小且容易实现的优点，被广泛应用与地震勘探技术中。有限差分的系数可以通过Taylor级数展开或是优化算法来求取。前者可以看作是对伪谱算子级数展开的截断，而后者可以看作是对伪谱算子在某个或某些特定频率上的最小二乘拟合。有限差分本质上是对伪谱算子的一种近似，这种近似使得有限差分存在计算精度低，频散严重以及计算不稳定等问题。在实际计算过程中，为了提高有限差分的计算精度，往往采用高阶的空间差分算子，但由于计算机存储等限制，在时间上通常还是采用二阶精度的差分，因此波场延拓在时间方向上精度较低，这也导致了频散误差严重。另外，为了保证计算的稳定，需要采用比地震采集数据小得多的时间步长进行波场延拓，增加了计算量。

近年来随着计算机技术的发展，在波数-空间域进行波场延拓成为可能。在波数-空间域构建的波场延拓算子能够补偿时间离散引起的误差，即便采用大时间步长延拓，仍能保持极高的精度和稳定性，对于依赖于波场延拓方法的地震正演、成像和反演技术的发展具有重要意义。但是由于实际地震数据量巨 大，直接利用波数-空间域算子处理实际数据依然受到计算速度和计算存储的限制。将空间波数域算子和有限差分方法相结合，发挥空间波数域的波场延拓算子精度高、稳定性好的优点以及有限差分计算速度快的优点，是解决波场延拓问题的一种新思路。交错网格Lowrank有限差分方法就是基于这种思想提出的一种新的波场延拓方法。该方法利用Lowrank分解处理波数-空间域的算子，得到有限差分形式的计算格式，在保证精度的前提下，大大节省了计算量，并且发挥交错网格的特点，在不增加计算量的前提下，提高计算精度。虽然该方法精度很高，但时间步长的选择依然受到稳定性条件的限制。一种在大时间步长延拓时，能保证计算稳定的高精度波场延拓方法对于实际地震资料的处理具有重要价值。

发明内容

为了得到适用于实际地震资料处理且计算稳定的高精度地震波场延拓方法，本发明提出了一种基于交错网格Lowrank算子的无条件稳定有限差分地震波场延拓方法，通过对交错网格Lowrank算子加入衰减约束提高稳定性，通过考虑了频率权重的加权最小二乘拟合提高计算精度，最后得到对任意实际步长都计算稳定的有限差分系数，并将其用于地震波场的延拓。

为实现上述目的，本发明的技术方案如下：

基于交错网格Lowrank算子的无条件稳定有限差分波场延拓方法，其特征在于，包括以下步骤：

步骤1：利用一阶速度-应力方程的Fourier积分解构建交错网格上的Lowrank算子

步骤2：利用衰减函数对交错网格Lowrank算子进行衰减约束

步骤3：利用加权最小二乘方法计算有限差分系数

步骤4：利用得到的有限差分系数实现无条件稳定的有限差分波场延拓。

相对于现有技术，本发明的有益效果如下：能够实现无条件稳定的高精度的地震波场延拓，有效压制数值频散，在达到相同计算精度的条件下，预期能提高2至3倍的计算效率。

附图说明

图1是基于交错网格Lowrank算子的无条件稳定有限差分波场延拓方法的流程图。

图2是衰减函数曲线。

图3是频率权重函数曲线。

图4是常规有限差分方法得到的地震波场。

图5是加入衰减约束后的交错网格Lowrank有限差分方法得到的地震波场。

图6是同时加入衰减约束和频率权重后的交错网格Lowrank有限差分方法得到的地震波场。

具体实施方式

如图1所示，基于交错网格Lowrank算子的无条件稳定有限差分波场延拓方法，其特征在于，包括以下步骤：

步骤1：利用一阶速度-应力方程的Fourier积分解构建交错网格上的Lowrank算子

根据地下介质的速度和密度模型以及波场延拓参数，利用一阶速度-应力方程的Fourier积分解在交错网格上构建波数-空间域的波场延拓算子；对波数-空间域波场延拓算子应用Lowrank分解，得到交错网格Lowrank算子。具 体方法如下：

地震波场在地下介质中的传播可以近似的用声波方程描述。声波方程有多种形式，其中一阶速度-应力声波方程考虑了速度和密度变化，适用于地震波场的延拓。在速度和密度变化的介质中，一阶速度-应力方程为

$[\frac{\∂}{\∂t}+A]u(x,t)=0---\left(1\right)$

其中 $u(x,t)=\left[\begin{array}{c}u(x,t)\\ p(x,t)\end{array}\right]$ 为矢量波场，u表示质点震动的速度，p表示压力，A为

$A=\left[\begin{array}{cc}0& \frac{1}{\mathrm{\ρ}\left(x\right)}\▿\\ \mathrm{\ρ}\left(x\right)v^2\left(x\right)\▿\·& 0\end{array}\right]---\left(2\right)$

其中利用一阶速度-应力方程(1)的Fourier积分解，可以构建如下的波数-空间域波场延拓算子

式(3)给出了在波数-空间域实现波场延拓的一般形式。如果当前时刻为t，那么下一时刻的地震波场u(x,t+Δt)可以通过对当前时刻的波场作用一个波数-空间域的积分算子再加上一时刻的波场u(x,t-Δt)得到。因此实现式(3)给出的波场延拓的核心是构建波数-空间域的积分算子。本发明的步骤1利用一阶速度-应力方程的Fourier积分解给出交错网格上波数-空间域的波场延拓算子，这是后续步骤的基础。

式(1)中时间和空间上的偏导数对应的离散形式为

$\frac{\∂}{\∂t}p(x,t)\≈\frac{p(x,t+\mathrm{\Δt}/2)-p(x,t-\mathrm{\Δt}/2)}{\mathrm{\Δt}}---\left(4\right)$

$\frac{\∂}{\∂x}p(x+\mathrm{\Δx}/2,t)\≈D_x^{+}P(x,t)---\left(5\right)$

$\frac{\∂}{\∂x}p(x-\mathrm{\Δx}/2,t)\≈D_x^{-}P(x,t)---\left(6\right)$

$\frac{\∂}{\∂z}p(x+\mathrm{\Δz}/2,t)\≈D_z^{+}P(x,t)---\left(7\right)$ $\frac{\∂}{\∂z}p(x-\mathrm{\Δz}/2,t)\≈D_z^{-}P(x,t)---\left(8\right)$

将式(4)-(8)代入(1)式，可得交错网格上一阶速度-应力方程的离散形式

$\frac{u_x(x_1,t^{+})-u_x(x_1,t^{-})}{\mathrm{\Δt}}=-\frac{1}{\mathrm{\ρ}\left(x_1\right)}{D}_x^{+}p(x,t)$

$\frac{u_z(x_2,t^{+})-u_x(x_2,t^{-})}{\mathrm{\Δt}}=-\frac{1}{\mathrm{\ρ}\left(x_2\right)}{D}_z^{+}p(x,t)---\left(9\right)$

$\frac{p(x,t+\mathrm{\Δt})-p(x,t)}{\mathrm{\Δt}}=-\mathrm{\ρ}\left(x\right)v^2\left(x\right)[D_z^{-}{u}_x(x_1,t^{+})+D_z^{-}{u}_z(x_2,t^{+})]$

其中

x₁＝(x+Δx/2,z)，x₂＝(x,z+Δz/2)，

t⁺＝t+Δt/2,t^-＝t-Δt/2

下面推导计算x和z方向的空间偏导数的波数-空间域算子，即式(3)中W(x,k)涉及到的一阶偏导数算子。对应于方程(1)的二阶波动方程为

$\▿\·\frac{1}{\mathrm{\ρ}\left(x\right)}\▿p(x,t)-\frac{1}{\mathrm{\ρ}\left(x\right)v^2\left(x\right)}\frac{{\∂}^{2}p(x,t)}{{\∂t}^2}=0---\left(10\right)$

对于常速度和常密度的情况，方程(10)可以写成如下波数域的形式

$\frac{{\∂}^2\hat{p}(k,t)}{{\∂t}^2}=-v_0^2{k}^2\hat{p}(k,t)---\left(11\right)$

其中是波场p(k,t)的空间傅里叶变换。方程(11)的解析解为 

$\hat{p}(k,t+\mathrm{\Δt})=e^{\±{\mathrm{iv}}_0\left|k\right|\mathrm{\Δt}}\hat{p}(k,t)---\left(12\right)$

将此解析解带入二阶时间导数微分格式中，可得

$\begin{array}{c}\frac{\hat{p}(k,t+\mathrm{\Δt})-2\hat{p}(k,t)+\hat{p}(k,t-\mathrm{\Δt})}{{\mathrm{\Δt}}^2}\\ =-{\left(v_0\right|k\left|\right)}^2{\sin c}^2\left(v_0\right|k|\mathrm{\Δt}/2)\hat{p}(k,t)\end{array}---\left(13\right)$

其中sinc(x)＝sin(x)/x。

通常情况下，速度和密度都是随着空间位置变化的，在速度梯度和时间步 长较小的情况下，可以将(13)式中的常速度v₀替换为变速度v(x)，从而得到近似式

$\begin{array}{c}\frac{\hat{p}(k,t+\mathrm{\Δt})-2\hat{p}(k,t)+\hat{p}(k,t-\mathrm{\Δt})}{\mathrm{\Δ}{t}^2}\\ \≈-{\left(v\left(x\right)\right|k\left|\right)}^2\sin c^2\left(v\left(x\right)\right|k|\mathrm{\Δt}/2)\hat{p}(k,t)\end{array}---\left(14\right)$

对上式两端应用逆傅里叶变换，可以得到波场在时间-空间域的递推格式

$\begin{array}{c}\frac{p(x,t+\mathrm{\Δt})-2p(x,t)+p(x,t-\mathrm{\Δt})}{\mathrm{\Δ}{t}^2}\\ \≈-v^2\left(x\right)F^{-1}\left[\right|k{|}^{2}s{\mathrm{inc}}^2\left(v\left(x\right)\right|k|\mathrm{\Δt}/2)F\left[p\right](x,t)\end{array}---\left(15\right)$

其中F表示空间傅里叶变换。上式右端可看作是一个依赖于空间位置x和波数k的算子作用于波场，记为

$[{\▿}^{v\left(x\right)\mathrm{\Δt}}{]}^{2}p(x,t)\≡-F^{-1}[\left|k{|}^2{\sin c}^2\left(v\left(x\right)\right|k|\mathrm{\Δt}/2)F\right[p(x,t)]---\left(16\right)$

利用(16)式，(15)式可以表示为 

$\frac{p(x,t+\mathrm{\Δt})-2p(x,t)+p(x,t-\mathrm{\Δt})}{\mathrm{\Δ}{t}^2}\≈{-v}^2\left(x\right)[{\▿}^{v\left(x\right)\mathrm{\Δt}}{]}^{2}p(x,t)---\left(17\right)$

在(16)和(17)式中，算子相当于标准的梯度算子，但是该算子还依赖于参数v(x)和Δt。当速度为常数时，(17)式严格成立，这意味着算子补偿了由于时间离散而导致的误差。事实上，对于变化的速度场，算子也能在一定程度上补偿时间离散引起的误差。

对于地震波场这样的有限带宽信号，通过对波场应用傅里叶变换可以准确求得波场在空间上的偏导数。类似于常规的梯度算子，可以定义

$[{\▿}^{v\left(x\right)\mathrm{\Δt}}{]}^2\≡D_x^{+}{D}_x^{-}+D_z^{+}{D}_z^{-}---\left(18\right)$

其中x，z方向的偏导数算子为

$D_x^{+}p(x,t)\≡F^{-1}\left[i{k}_x{e}^{i{k}_x\mathrm{\Δx}/2}\sin c\left(v\left(x\right)\right|k|\mathrm{\Δt}/2)F\right[p(x,t)\left]\right]$

$D_x^{-}p(x,t)\≡F^{-1}\left[i{k}_x{e}^{-i{k}_x\mathrm{\Δx}/2}\sin c\left(v\left(x\right)\right|k|\mathrm{\Δt}/2)F\right[p(x,t)\left]\right]---\left(19\right)$

$D_z^{-}p(x,t)\≡F^{-1}\left[i{k}_z{e}^{i{k}_z\mathrm{\Δz}/2}\sin c\left(v\left(x\right)\right|k|\mathrm{\Δt}/2)F\right[p(x,t)\left]\right]$

$D_z^{-}p(x,t)\≡F^{-1}\left[i{k}_z{e}^{i{k}_z\mathrm{\Δz}/2}\sin c\left(v\left(x\right)\right|k|\mathrm{\Δt}/2)F\right[p(x,t)\left]\right]$

其中F表示空间傅里叶变换，_kx和k_z为波数分量，满足式中指数项使得波场在计算时沿着x轴的正向或负向交错Δx/2个网格步长，或者沿着z轴的正向或负向交错Δz/2个网格步长。 

直接利用公式(19)计算空间偏导数，需要的计算量为，这对于实际地震数据的处理是难以承受的。下面采用Lowrank分解对公式(19)进行近似处理，已降低计算量。以式(19)中x方向的空间偏导数为例，定义波数-空间域算子矩阵

W_x(x,k)≡k_xsinc(v(x)|k|Δt/2)   (20) 

对式(20)应用Lowrank分解，将算子矩阵分解为三个矩阵的乘积

其中W_x(x,k_m)和W_x(x_n,k)分别是从矩阵W_x(x,k)选取M列和N行组成的子矩阵，a_mn是中间矩阵的系数，这里M和N是依赖于矩阵秩的整数，通常M和N远小于计算点数N_x。因此，可以得到计算空间偏导数的交错网格Lowrank算子

$D_x^{+}p(x,t)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{W}_x(x,k_m)a_{\mathrm{mn}}{F}^{-1}\left[{\mathrm{ie}}^{{\mathrm{ik}}_x\mathrm{\Δx}/2}{W}_x(x_n,k)\mathrm{taper}\left(k\right)F\right[p(x,t)\left]\right]---\left(22\right)$

类似的可以得到其余计算x和z方向的空间偏导数的交错网格Lowrank算子

$D_x^{-}p(x,t)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{W}_x(x,k_m)a_{\mathrm{mn}}{F}^{-1}\left[{\mathrm{ie}}^{{\mathrm{ik}}_x\mathrm{\Δx}/2}{W}_x(x_n,k)\mathrm{taper}\left(k\right)F\right[p(x,t)\left]\right]$

$D_z^{+}p(x,t)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{W}_z(x,k_m)a_{\mathrm{mn}}{F}^{-1}\left[{\mathrm{ie}}^{{\mathrm{ik}}_x\mathrm{\Δx}/2}{W}_z(x_n,k)\mathrm{taper}\left(k\right)F\right[p(x,t)\left]\right]---\left(23\right)$

$D_z^{-}p(x,t)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{W}_z(x,k_m)a_{\mathrm{mn}}{F}^{-1}\left[{\mathrm{ie}}^{{\mathrm{ik}}_x\mathrm{\Δx}/2}{W}_z(x_n,k)\mathrm{taper}\left(k\right)F\right[p(x,t)\left]\right]$

步骤2：对交错网格Lowrank算子进行衰减约束

构造波数相关的光滑衰减函数；利用衰减函数对步骤1构建的交错网格Lowrank算子进行衰减约束。具体方法如下：

为了提高大时间步长波场延拓的稳定性，我们给出对传播算子进行衰减约束的方法。交错网格Lowrank算子式(22)-(23)是一个光滑震荡的函数，根据Von-Neumann稳定性分析易知，计算稳定的充分条件是交错网格Lowrank算子的取值对于任意情况都位于区间[-1,+1]内。但在实际计算过程中，数值误差使得交错网格Lowrank算子在时间步长选择较大时，往往不能满足稳定性条件。在为此需对交错网格Lowrank算子施加衰减约束，施加衰减约束之后的的交错网格Lowrank算子为

$D_x^{+}p(x,t)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{W}_x(x,k_m)a_{\mathrm{mn}}{F}^{-1}\left[{\mathrm{ie}}^{{\mathrm{ik}}_x\mathrm{\Δx}/2}{W}_x(x_n,k)\mathrm{taper}\left(k\right)F\right[p(x,t)\left]\right]$

$D_x^{-}p(x,t)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{W}_x(x,k_m)a_{\mathrm{mn}}{F}^{-1}\left[{\mathrm{ie}}^{{\mathrm{ik}}_x\mathrm{\Δx}/2}{W}_x(x_n,k)\mathrm{taper}\left(k\right)F\right[p(x,t)\left]\right]---\left(24\right)$

$D_z^{+}p(x,t)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{W}_z(x,k_m)a_{\mathrm{mn}}{F}^{-1}\left[{\mathrm{ie}}^{{\mathrm{ik}}_x\mathrm{\Δx}/2}{W}_z(x_n,k)\mathrm{taper}\left(k\right)F\right[p(x,t)\left]\right]$

$D_z^{-}p(x,t)=\underset{m=1}{\overset{M}{\mathrm{\Σ}}}\underset{n=1}{\overset{N}{\mathrm{\Σ}}}{W}_z(x,k_m)a_{\mathrm{mn}}{F}^{-1}\left[{\mathrm{ie}}^{{\mathrm{ik}}_x\mathrm{\Δx}/2}{W}_z(x_n,k)\mathrm{taper}\left(k\right)F\right[p(x,t)\left]\right]$

其中taper(k)为和波数相关的衰减函数。这里衰减函数的选取应该满足如下三个条件，

①衰减函数是关于波数k的光滑函数； 

②能够自动识别出交错网格Lowrank算子的稳定临界条件；

③能在临界值附近加入可控的衰减。

我们选用取如下形式的衰减函数，

$\mathrm{taper}\left(k\right)=\left\{\begin{array}{cc}1-\mathrm{\&beta;exp}[1-\frac{1}{1-{\left(\frac{k-k_0}{\mathrm{\&alpha;}/2}\right)}^2}],& \mathrm{if}|k-k_0|<\mathrm{\&alpha;}/2\\ 1-\mathrm{\&beta;exp}[1-\frac{1}{1-{\left(\frac{k+k_0}{\mathrm{\&alpha;}/2}\right)}^2}],& \mathrm{if}|k+k_0|<\mathrm{\&alpha;}/2\\ 1,& \mathrm{else}\end{array}\right.$

易知taper(k)关于波数k的任意阶导函数存在(即该函数是光滑的)，通常k的取值范围是从负Nyquist波数到正Nyquist波数。式中k₀，α和β是三个控制参数，k₀为特征波数，是使得W(x,k)达到稳定条件所限制的临界值的波数，α控制衰减的宽度，β控制衰减的大小。图2给出了取不同控制参数的衰减函数图像，k₀取为0.175，实线所示的衰减函数α＝300，β＝0.02，虚线所示的衰减函数α＝150，β＝0.04。

步骤3：利用加权最小二乘方法计算有限差分系数

构造频率加权函数；利用加权最小二乘拟合步骤2提出的加衰减约束交错网格Lowrank算子(24)，计算得到有限差分系数。具体方法如下：

采用式(24)进行波场延拓，需要在每一个时刻上做多次空间Fourier变换，这使得该方法的并行实现较为困难。下面根据式(24)给出易用于分布式内存并行计算平台的有限差分方法计算格式。

以式(24)中x方向的偏导数为例，利用Fouirer变化的相移性质，可将计算偏导数的Fouirer积分算子转化到空间域求解，在空间域其求解格式具有如下的有限差分形式，

$D_x^{+}p(x,t)=\frac{1}{2}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}G(x,l)[p(x_R,t)-p(x_L,t)]---\left(26\right)$

其中p(x,t)为定义在空间位置x和时间t上的波场，x_R＝(x+l¹Δx,z+l²Δz)，x_L＝(x-(l¹-1)Δx,z+l²Δz)，G(x,l)为有限差分的差分系数，可由下面的最小二乘问题求得

$\min \left|\right|[\mathrm{taper}\left(k\right)W(x,k)-\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}G(x,l)B(l,k)]|{|}_2---\left(27\right)$

其中W(x,k)＝sin(v(x)|k|Δt/2)，有限差分系数G(x,l)是依赖于空间位置的，随着不同位置模型参数的不同自适应地变化。

通过在波数k的整个取值范围内拟合交错网格Lowrank算子求取有限差分系数，在实际的计算中k的取值范围是从负Nyquist频率到正Nyquist频率，Nyquist频率的大小由空间步长和网格点数通过采样定理确定。但是地震波场是带限的，地震波场中有效信号的绝大部分能量都集中在主频附近的一个小的频带范围内，这个有效的频带范围仅是正、负Nyquist频率之间很小的一个区间，在这个区间内波场延拓算子的精度决定了波场延拓过程中的计算精度。因此为了提高计算精度，在总体拟合精度不变的前提下，可以让主频附近的拟合精度尽可能高一些，而容许在远离主频的波数误差可大一些，因为这部分的有效能量很小，即便是其误差较大对整体的计算精度的影响也不会很大。为此要求在做最小二乘拟合时，对主频附近的波数设定较大的权重，以保证这部分算子的精度。

基于加入频率权重和衰减约束的交错网络Lowrank算子的优化有限差分系数可以通过求解加权最小二乘问题得到，

$\min \left|\right|[\mathrm{taper}\left(k\right)W(x,k)-\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}G(x,l)B(l,k)]|{|}_2---\left(28\right)$

其中w(k)为Gauss型的双峰权函数

$w\left(k\right)=\frac{1}{2}[e^{-\frac{{(k-k_0)}^2}{a}}+e^{-\frac{{(k+k_0)}^2}{a}}]---\left(29\right)$

其中f₀为地震子波主频，v为介质速度。图3给出了两个Gauss型权函数的图像，介质速度为4.5 km/s，其中实线表示的是主频为20Hz，a＝20000时的Gauss权函数，虚线表示主频为30Hz，a＝2000时的Gauss权函数。

式(28)给出的最小二乘问题的解为

其中为一个N_k×N_k的对角矩阵，其元素为权重函数。

通过衰减约束保证式(30)给出的有限差分系数G对于任意大的时间步长计算都是稳定的，求取G时所用的加权最小二乘方法保证了有限差分系数的精度。

步骤4：利用得到的有限差分系数实现无条件稳定的有限差分波场延拓

将步骤3得到的优化差分系数(30)用于波场延拓可实现无条件稳定的地震波场延拓。采用步骤3中得到优化有限差分系数，用下式实现波场延拓过程中沿着x轴正方向空间偏导数的计算，

$D_x^{+}p(x,t)=\frac{1}{2}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}G(x,l)[p(x_R,t)-p(x_L,t)]---\left(31\right)$

其中x_R＝(x+l¹Δx,z+l²Δz)，x_L＝(x-(l¹-1)Δx,z+l²Δz)，G(x,l)为步骤3中得到优化有限差分系数，l¹,l²,l³＝1,L,L表示有限差分的结构，类似的对于沿着x轴负方向空间偏导数，其有限差分计算格式为

$D_x^{-}p(x,t)=\frac{1}{2}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}G(x,l)[p(x_R,t)-p(x_L,t)]---\left(32\right)$

其中x_R＝(x+(l¹-1)Δx,z+l²Δz)，x_L＝(x-l¹Δx,z+l²Δz)。对于沿着z轴正方向交错的空间偏导数，其有限差分计算格式为

$D_z^{+}p(x,t)=\frac{1}{2}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}G(x,l)[p(x_U,t)-p(x_D,t)]---\left(33\right)$

其中x_U＝(x+l¹Δx,z+l²Δz)，x_D＝(x+l¹Δx,z-(l²-1)Δz)。对于沿着z轴负方向交错的空间偏导数，其有限差分计算格式为

$D_z^{-}p(x,t)=\frac{1}{2}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}G(x,l)[p(x_U,t)-p(x_D,t)]---\left(34\right)$

其中x_U＝(x+l¹Δx,z+(l²-1)Δz)，x_D＝(x+l¹Δx,z-l²Δz)。

将式(31)-(34)式带入式(9)，可得到无条件稳定的有限差分波地震场延拓计算式，

$u_x(x_1,t^{+})=u_x(x_1,t^{-})-\frac{\mathrm{\&Delta;t}}{2\mathrm{\&rho;}\left(x_1\right)}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}G(x,l)[p(x_R,t)-p(x_L,t)]$

$u_z(x_2,t^{+})=u_x(x_2,t^{-})-\frac{\mathrm{\&Delta;t}}{2\mathrm{\&rho;}\left(x_2\right)}\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}G(x,l)[p(x_U,t)-p(x_D,t)]$

$\begin{array}{c}p(x,t+\mathrm{\&Delta;t})=p(x,t)-\frac{\mathrm{\&rho;}\left(x\right)v^2\left(x\right)\mathrm{\&Delta;t}}{2}\left[\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}G(x,l)\right[u_x(x_R,t^{+})-u_x(x_L,t^{+})]\\ +\underset{l=1}{\overset{L}{\mathrm{\&Sigma;}}}G(x,l)[u_z(x_U,t^{+})-u_z(x_D,t^{+})]]\end{array}---\left(35\right)$

其中

x₁＝(x+Δx/2,z)，x₂＝(x,z+Δz/2)，

t⁺＝t+Δt/2,t^-＝t-Δt/2

地震勘探通常采用4ms的时间采用间隔记录地震信号，常规方法计算不稳定采用4ms时间步长进行数值波场延拓时计算不稳定，如图4所示。图5所示为仅采用本发明步骤2中提出的衰减约束得到的波场，衰减约束能够使得计算稳定，但当传播时间较大时，会出现数值误差。图6所示为在计算过程中同时加入本发明步骤2中提出的衰减约束和步骤3中提出的频率权重得到的波场 记录，可以看出当使用4ms时间步长时，计算稳定，而且得到的波场精确，没有数值噪音。

Claims (5)

1.一种基于交错网格Lowrank算子的无条件稳定有限差分地震波场延拓方法，包括以下步骤：

步骤1：利用一阶速度-应力纵波方程的Fourier积分解构建交错网格上的Lowrank算子

步骤2：利用衰减函数对交错网格Lowrank算子进行衰减约束

步骤3：利用加权最小二乘方法计算有限差分系数

步骤4：利用得到的有限差分系数实现无条件稳定的有限差分波场延拓。

2.根据权利要求1所述的基于交错网格Lowrank算子的无条件稳定有限差分波场延拓方法，其特征在于，步骤1为根据地下介质的速度和密度模型以及波场延拓参数，利用一阶速度-应力纵波方程的Fourier积分解在交错网格上构建波数-空间域的波场延拓算子；对波数-空间域波场延拓算子应用Lowrank分解，得到交错网格Lowrank算子；具体方法如下：

一阶速度-应力纵波方程可以描述地震波场在声波介质中的传播。随速度和密度变化而变化的一阶速度-应力纵波方程为

其中矢量波场u表示质点振动的速度，p表示压力，A为

上式中时间和空间上的偏导数对应的离散形式为

在交错网格上，一阶速度-应力纵波方程的离散形式为

其中

x₁＝(x+Δx/2,z)，x₂＝(x,z+Δz/2)，

t⁺＝t+Δt/2,t^-＝t-Δt/2

利用Lowrank分解，用于计算偏导数的交错网格Lowrank算子为

。

3.根据权利要求1所述的基于交错网格Lowrank算子的无条件稳定有限差分波场延拓方法，其特征在于，步骤2为构造波数相关的光滑衰减函数；利用衰减函数对步骤1构建的交错网格Lowrank算子进行衰减约束；具体方法如下：

利用衰减函数，对交错网格Lowrank算子中进行衰减约束，加入衰减约束后的偏导数算子为

其中taper(k)为构造衰减函数，形式如下

其中中k₀，α和β是三个控制参数，k₀为特征波数，是使得W(x,k)达到稳定条件所限制的临界值的波数，α控制衰减的宽度，β控制衰减的大小。

4.根据权利要求1-3所述的基于交错网格Lowrank算子的无条件稳定有限差分波场延拓方法，其特征在于，步骤3为利用加权最小二乘拟合步骤2提出的加衰减约束交错网格Lowrank算子，得到优化有限差分系数；具体方法如下：

以式中x方向的偏导数为例，利用Fouirer变化的相移性质，可将计算偏导数的Fouirer积分算子转化到空间域求解，在空间域其求解格式具有如下的有限差分形式，

其中p(x,t)为定义在空间位置x和时间t上的波场，x_R＝(x+l¹Δx,z+l²Δz)，x_L＝(x-(l¹-1)Δx,z+l²Δz)，G(x,l)为有限差分的差分系数，可由下面的最小二乘问题求得

其中W(x,k)＝sin(v(x)|k|Δt/2)，G(x,l)为有限差分系数，是依赖于空间位置的，随着不同位置模型参数的不同自适应地变化。w(k)为Gauss型的双峰权函数，取为

其中f₀为地震子波主频，v为介质速度。该最小二乘问题的解为

其中为一个N_k×N_k的对角矩阵，其元素为权重函数。

通过衰减约束保证G给出的有限差分系数对于任意大的时间步长计算都是稳定的，求取G时所用的加权最小二乘方法保证了有限差分系数的精度。

5.根据权利要求1-4所述的基于交错网格Lowrank算子的无条件稳定有限差分波场延拓方法，其特征在于，步骤4为:将步骤3得到的优化差分系数用于波场延拓可实现无条件稳定的地震波场延拓，具体方法如下，

采用步骤3中得到优化有限差分系数，用下式实现波场延拓过程中沿着x轴正方向空间偏导数的计算，

其中x_R＝(x+l¹Δx,z+l²Δz)，x_L＝(x-(l¹-1)Δx,z+l²Δz)，G(x,l)为步骤3中得到优化有限差分系数，l¹,l²＝1,…,L表示有限差分的结构，类似的对于沿着x轴负方向空间偏导数，其有限差分计算格式为

其中x_R＝(x+(l¹-1)Δx,z+l²Δz)，x_L＝(x-l¹Δx,z+l²Δz)。对于沿着z轴正方向交错的空间偏导数，其有限差分计算格式为

其中x_U＝(x+l¹Δx,z+l²Δz)，x_D＝(x+l¹Δx,z-(l²-1)Δz)。对于沿着z轴负方向交错的空间偏导数，其有限差分计算格式为

其中x_U＝(x+l¹Δx,z+(l²-1)Δz)，x_D＝(x+l¹Δx,z-l²Δz)。

将上述计算偏导数的格式带入步骤1中的一阶速度-应力纵波方程的离散形式，得到无条件稳定的有限差分波地震场延拓计算式，

其中

x₁＝(x+Δx/2,z)，x₂＝(x,z+Δz/2)，

t⁺＝t+Δt/2,t^-＝t-Δt/2。