CN106932820A - 基于时域伪谱方法的声波方程逆时偏移成像方法 - Google Patents

Abstract

基于时域伪谱方法的声波方程逆时偏移成像方法，涉及地球物理勘探。根据初始模型的具体情况结合现有计算资源状况，合理划分炮孔径及炮间距，利用计算资源、计算时间和精度，设定每炮的间距相同、炮孔径相同、震源位置及深度相同；对每一炮执行处理；所有炮单独处理完毕后，把所有炮的成像结果按具体划分情况，依次进行叠加，形成整体成像结果；在成像结果上应用高通滤波器对低频噪声进行过滤，提高成像分辨率。数值仿真是逆时偏移成像的最重要组成部分。逆时偏移成像需要对双程波动方程进行波场延拓，求解过程中需要同时实现震源波场和检波器波场的数值迭代和互相关，这都需要使用某种数值方法来具体实现。

Description

基于时域伪谱方法的声波方程逆时偏移成像方法

技术领域

本发明涉及地球物理勘探，尤其是涉及基于时域伪谱方法的声波方程逆时偏移成像方法。

背景技术

逆时偏移(Reverse time migration,RTM)方法是近年来油气勘探领域的研究热点之一，也是目前成像精度最高的地震波叠前深度偏移方法之一。逆时偏移的基本思想最早在1978年由Hemon提出[1]。随后，Baysal[2]，Whitmore[3]及McMECHAN[4]分别提出了对逆时偏移思想的具体应用。与其他成像方法相比，逆时偏移成像方法具有成像精度高、适应性广等特点，尤其是在高倾角、水平速度变化大的复杂地质结构成像中具有明显的优势，能够对常规地震成像方法中无法成像的复杂地质结构进行高精度成像。

逆时偏移成像中所需的庞大资源仍然是制约着逆时偏移成像技术的广泛应用的主要瓶颈之一，迫切需要提高逆时偏移的效率，降低成本。要缓解计算瓶颈，可以考虑使用高效快速算法来更快更好的求解波动方程，从而减少所需的计算资源。时域伪谱方法(Thepseudospectral time-domainalgrithm，PSTD)是用于求解偏微分方程常用数值方法之一[5],[6]。时域伪谱方法的基本思想是：时间域导数利用差分近似，空间域导数则使用快速傅里叶变换求取，并通过时间的递推来模拟波的传播过程。时域伪谱方法的原理简单，网格剖分不需要借助辅助工具，介质参数亦不需要做平均化处理，傅里叶变换借助现有的开源傅里叶变换库可以快捷高效实现。本发明使用时域伪谱方法对地下复杂介质进行逆时偏移成像，与传统的二阶差分和高阶差分方法相比，不仅能保证计算所需的精度，而且能有效的降低空间采样率，从而节省存储空间、提高计算效率[7]。

IO带宽瓶颈也是逆时偏移成像中需要考虑的重要问题。传统的逆时偏移成像需要把其中一个波场在所有时刻的快照保存到磁盘，在成像过程中再读取出来，这很大程度上成为整个成像过程的瓶颈。为了解决这个问题，采用Clapp[8]和Dussaud[9]提出波场保存方法，该方法先将源波场进行外推，每个时间步只保存源波场中靠近计算区域边界的若干层，有效地缓解了IO瓶颈。

发明内容

本发明的目的在于提供高效快速求解逆时偏移成像问题的基于时域伪谱方法的声波方程逆时偏移成像方法。

本发明包括以下步骤：

1)根据初始模型的具体情况结合现有计算资源状况，合理划分炮孔径及炮间距，利用计算资源、计算时间和精度，设定每炮的间距相同、炮孔径相同、震源位置及深度相同；

2)对每一炮执行处理的具体方法可为：

(1)在现场实验中获得检波器数据或使用真实模型仿真产生检波器数据，对于一阶声波方程，需要获得声压信号；

(2)使用初始模型，在指定炮点放置震源信号，并应用一阶声波方程进行数值仿真，对震源波场进行正向波场延拓直到最大时间步，即从t₀到t_max，同时保存每个时间步的边界波场，根据模型大小，每个方向各保存3～5层不等。

在步骤2)第(2)部分中，所述一阶声波方程，是采用关于质点速度矢量v和压力p的运动方程和质量守恒方程；采用时域伪谱方法作为数值仿真方法，作用在函数u(j_η)(η＝x,y,z)上的伪谱求导算子定义如下：

其中F_η是正快速傅里叶变换，是逆快速傅里叶变换，k_η是空间波数。

(3)使用初始模型，并利用保存的边界波场作为边界条件，使用数值方法仿真，对震源波场进行逆时外推，即从t_max到t₀；同时使用对检波器波场作为边界条件注入，使用数值方法仿真，对检波器波场进行逆时外推，即从t_max到t₀；仿真使用的公式和仿真方法同步骤(2)；

(4)对步骤(2)、(3)中的震源波场和检波器波场应用震源补偿归一化互相关成像条件；所述震源补偿归一化互相关成像条件如下：

其中，S(x,y,z,t)是震源波场，R(x,y,z,t)是检波器波场，time从t₀到t_max是时间步数；

3)所有炮单独处理完毕后，把所有炮的成像结果按步骤1)的具体划分情况，依次进行叠加，形成整体成像结果；

4)在成像结果上应用高通滤波器对低频噪声进行过滤，提高成像分辨率。

数值仿真是逆时偏移成像的最重要组成部分。逆时偏移成像需要对双程波动方程进行波场延拓，求解过程中需要同时实现震源波场和检波器波场的数值迭代和互相关，这都需要使用某种数值方法来具体实现。因此数值方法的实现很大程度上决定了逆时偏移成像的效率和质量。

与现有的统逆时偏移成像方法相比，本发明具有以下技术效果：

1)针对求解双程波动方程所需的庞大计算时间和存储资源，用时域伪谱方法(Thepseudospectral time-domain algrithm，PSTD)取代传统的时域有限差分方法，大幅度降低了所需的计算时间和所需内存。时域伪谱方法的基本思想是：时间域导数利用差分近似，空间域导数则使用快速傅里叶变换求取，并通过时间的递推来模拟波的传播过程。时域伪谱方法的原理简单，网格剖分不需要借助辅助工具，介质参数亦不需要做平均化处理，同时傅里叶变换借助现有的开源傅里叶变换库可以快捷高效实现。本发明使用时域伪谱方法对地下复杂介质进行逆时偏移成像，与现有的二阶差分和高阶差分方法相比，不仅能保证计算所需的精度，而且能有效的降低空间和时间采样率，从而节省存储空间、提高计算效率。

2)针对成像过程中所需的存储资源和输入输出带宽，先将源波场进行外推，每个时间步只保存源波场中靠近计算区域边界的若干层，然后再把源波场逆推，并把之前保存的边界波场注入以恢复源波场。这种方法有效的减少了所需的存储空间并缓解了IO瓶颈。

附图说明

图1为Marmousi II模型(仅P波)。

图2为Marmousi II初始模型。

图3为Marmousi II模型在空间采样率等于2时的成像结果。

图4为SEG/EAGE 3D Salt模型切片x＝7000m。

图5为SEG/EAGE 3D Salt模型切片y＝7000m。

图6为SEG/EAGE 3D Salt模型成像结果切片x＝7000m(空间采样率等于2)。

图7为SEG/EAGE 3D Salt模型成像结果切片y＝7000m(空间采样率等于2)。

具体实施方式

以下将结合附图对本发明的技术方案、成像效果和效率作进一步说明。

1)根据初始模型的具体情况结合现有计算资源状况，合理划分炮孔径及炮间距。既要考虑充分的利用计算资源，又要考虑计算时间和精度。所以炮间距不能过大，并满足必要的多次覆盖率。为了简便起见，我们设定每炮的间距相同，炮孔径相同，震源位置及深度也相同。

2)对每一炮执行以下处理步骤:

(1)在现场实验中获得检波器数据或使用真实模型仿真产生检波器数据。对于一阶声波方程，需要获得声压信号。

(2)使用初始模型，在指定炮点放置震源信号，并应用一阶声波方程进行数值仿真，对震源波场进行正向波场延拓直到最大时间步，即从t₀到t_max，同时保存每个时间步的边界波场(根据模型大小，每个方向各保存3-5层不等)。

(3)使用初始模型，并利用保存的边界波场作为边界条件，使用数值方法仿真，对震源波场进行逆时外推，即从t_max到t₀；同时使用对检波器波场作为边界条件注入，使用数值方法仿真，对检波器波场进行逆时外推，即从t_max到t₀。

(4)对步骤(2)、(3)中的震源波场和检波器波场应用震源补偿归一化互相关成像条件。

3)所有炮单独处理完毕后。把所有炮的成像结果按步骤1)的具体划分情况，依次进行叠加，形成整体成像结果。

4)在成像结果上应用高通滤波器对低频噪声进行过滤，提高成像分辨率。

本发明可应用于均匀各向同性介质中的声波逆时偏移成像。本发明能够在极低采样下(每最小波长2个采样点)得到高精度图像，因此能够节省大量存储空间和计算时间。

图1和图2为Marmousi II的P波速度分布。图1是真实模型；图2为仿真时使用的初始模型，由真实模型的每个介质点在附近100米矩形区域内取平均得到。

图3是Marmousi II模型采用本发明中的逆时偏移成像方法，采用时域伪谱方法在空间采样等于2时得到的成像结果，从结果来看，各主要的反射界面都得到了较清晰的图像。仿真使用时图2所示的速度模型，并使用中心频率为20Hz的Ricker子波作为震源信号。成像结果显示，主要反射界面均得到了较好的反映。基于时域伪谱方法的逆时偏移成像能够在极低采样下达到高精度的成像结果。

对于这个二维问题，采用时域伪谱方法，与传统的12阶差分方法相比，每个维度采样点数降低到二分之一，两个维度叠加就下降到了四分之一，即节省了75％的存储空间；同时，时域伪谱方法所需的计算时间步数降低到12阶差分方法的59％。

图4和图5为SEG/EAGE 3D Salt模型的速度分布。图4、图5是真实模型在不同角度的切片。

图6和图7是SEG/EAGE 3D Salt模型采用本发明中的逆时偏移成像方法，采用时域伪谱方法并在空间采样等于2时得到的成像结果在不同位置的切片。仿真使用中心频率为13.63赫兹的雷克子波作为震源信号。从成像结果的切片图6和图7中我们可以看到，即使在最低的采样率下，盐丘的各个位置都得到了清晰的成像，包括高倾角的区域。

对于这个三维问题，采用时域伪谱方法，与传统的12阶差分方法相比，每个维度采样点数降低到二分之一，三个维度叠加就下降到了八分之一，即节省了87.5％的存储空间；同时，时域伪谱方法所需的计算时间步数降低到12阶差分方法的59％。

经实践证明，采用时域伪谱方法实现的逆时偏移成像能大幅度降低所需的存储资源和计算时间，并能得到高精度的成像结果。

本发明涉及一种新的基于时域伪谱方法的声波方程逆时偏移成像方法，以时域伪谱方法替代传统的二阶中心差分和高阶差分方法。时域伪谱方法使用快速傅里叶变换来求取空间导数，再结合完美匹配层技术解决傅里叶变换隐含的周期性限制，从而能够快速高精度地求解逆时偏移成像问题。与现在的二阶中心差分和高阶差分方法相比，在同样精度下，采用时域伪谱方法可以大幅度减少存储空间并节省大量的计算时间。在存储空间上，采用时域伪谱方法与十二阶差分方法相比，在二维空间问题中可最多节省75％的存储空间，在三维空间问题中则可最多节省87.5％的存储空间；同时，时域伪谱方法所需的计算时间步数降低到十二阶差分方法中心差分方法的59％。本发明适用于复杂地下结构下高精度成像，并能够高效快速求解大尺度问题。在多个地震模型中应用了本发明中提出的逆时偏移方法，并验证了其精确、快速、计算成本低等特点。

Claims (2)

1.基于时域伪谱方法的声波方程逆时偏移成像方法，其特征在于包括以下步骤：

1)根据初始模型的具体情况结合现有计算资源状况，合理划分炮孔径及炮间距，利用计算资源、计算时间和精度，设定每炮的间距相同、炮孔径相同、震源位置及深度相同；

2)对每一炮执行处理；

3)所有炮单独处理完毕后，把所有炮的成像结果按步骤1)的具体划分情况，依次进行叠加，形成整体成像结果；

4)在成像结果上应用高通滤波器对低频噪声进行过滤，提高成像分辨率。

2.如权利要求1所述基于时域伪谱方法的声波方程逆时偏移成像方法，其特征在于在步骤2)中，所述对每一炮执行处理的具体方法为：

(1)在现场实验中获得检波器数据或使用真实模型仿真产生检波器数据，对于一阶声波方程，需要获得声压信号；

(2)使用初始模型，在指定炮点放置震源信号，并应用一阶声波方程进行数值仿真，对震源波场进行正向波场延拓直到最大时间步，即从t₀到t_max，同时保存每个时间步的边界波场，根据模型大小，每个方向各保存3～5层不；

所述一阶声波方程，是采用关于质点速度矢量v和压力p的运动方程和质量守恒方程；采用时域伪谱方法作为数值仿真方法，作用在函数u(j_η)(η＝x,y,z)上的伪谱求导算子定义如下：

$D_{\mathrm{\η},PS}\left(f\right(\mathrm{\η}\left)\right)=F_{\mathrm{\η}}^{-1}\{{\mathrm{ik}}_{\mathrm{\η}}{F}_{\mathrm{\η}}\[f\left(\mathrm{\η}\right)\]\},---\left(1\right)$

其中F_η是正快速傅里叶变换，是逆快速傅里叶变换，k_η是空间波数；

(3)使用初始模型，并利用保存的边界波场作为边界条件，使用数值方法仿真，对震源波场进行逆时外推，即从t_max到t₀；同时使用对检波器波场作为边界条件注入，使用数值方法仿真，对检波器波场进行逆时外推，即从t_max到t₀；仿真使用的公式和仿真方法同步骤(2)；

(4)对步骤(2)、(3)中的震源波场和检波器波场应用震源补偿归一化互相关成像条件；所述震源补偿归一化互相关成像条件如下：

$I_{RTM}(x,y,z)=\frac{\underset{time}{\mathrm{\Σ}}S(x,y,z,t)R(x,y,z,t)}{\underset{time}{\mathrm{\Σ}}{S}^2(x,y,z,t)}---\left(2\right)$

其中，S(x,y,z,t)是震源波场，R(x,y,z,t)是检波器波场，time从t₀到t_max是时间步数。