CN113447981B - 一种基于共成像点道集的反射全波形反演方法 - Google Patents

Abstract

本发明涉及一种基于共成像点道集和成像域偏移速度分析的反射全波形反演方法，包括以下步骤：1)对原始地震数据进行预处理；2)对预处理过的地震数据进行傅里叶变换，得到频率域反射炮记录d(ω，s，g)，其中，ω为圆频率，s为炮点位置，g为检波点位置；3)采用层析成像方法建立初始速度模型，并设定反演参数，包括成像点坐标与数量、频率范围和频率间隔；4)进行频率域声波方程反射全波形反演，获得最终速度剖面。与现有技术相比，本发明具有全自动无需交互、反演精度高、收敛效果好、实施简单等优点。

Description

一种基于共成像点道集的反射全波形反演方法

技术领域

本发明涉及勘探地震学速度建模领域，尤其是涉及一种基于共成像点道集的反射全波形反演方法。

背景技术

地震偏移成像是通过将地表观测到的地震数据重新归位到产生它的位置上，进而获取地球内部的结构特征。利用主动源的勘探地震观测系统通过滚动采集实现对地下成像点的多次覆盖，具有相当的冗余性，进而在成像过程中可以提取不同炮检对数据对同一成像位置的共成像点道集(CIG)。

早在上世纪八十年代Al-Yahya(1989)就提出“无论地下构造如何，如果宏观速度场正确，那么反射地震波偏移后得到的共成像点道集，同相轴就应该被拉平”。基于此原理，大量学者开展了利用层析成像将共成像点道集的非水平程度转化为背景速度修正量的偏移速度分析工作，取得了显著的效果，奠定了利用反射波进行中深层速度建模的基础。但这类方法需要交互拾取操作，影响效率与反演准确性。之后，建立在聚焦能量最大化或最佳相关性框架下的DSO(Symes&Carazzone,1991)与波动方程偏移速度分析(Sava&Biondi,2004)被提了出来。虽然避免了交互拾取，取得了明显的效果，但这类方法需要在成像扩展域实现，计算量非常大。随着全波形反演(FWI)研究突飞猛进，利用反射波形构建中深层速度模型的反射全波形反演方法近十年来成为研究热点之一(Xu et al.,2012；Chi et al.,2015)。这种方法无需交互拾取，计算量与FWI相同，但需要模型高波数先验信息，同时利用反偏移合成反射数据也很难保证动力学特征的准确性。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种基于共成像点道集的反射全波形反演方法。

本发明的目的可以通过以下技术方案来实现：

一种基于共成像点道集的反射全波形反演方法，包括以下步骤：

1)对原始地震数据进行预处理；

2)对预处理过的地震数据进行傅里叶变换，得到频率域反射炮记录d(ω，s，g)，其中，ω为圆频率，s为炮点位置，g为检波点位置；

3)采用层析成像方法建立初始速度模型，并设定反演参数；

4)进行频率域声波方程反射全波形反演，获得最终速度剖面。

所述的步骤1)中，预处理包括去噪、滤波、道均衡和去透射波。

所述的步骤3)中，反演参数包括成像点坐标与数量、频率范围和频率间隔

所述的步骤4)具体包括以下步骤：

41)根据观测系统统计不重复的炮点和检波点的数量及位置；

42)执行频率循环，在每一个频率循环内部，进行炮点端的单频正传波场u₀(s)、检波点端的单频格林函数G₀(g)以及成像点单频格林函数G₀(r)的计算，其中，r为成像点位置；

43)计算目标函数值并判断是否满足退出条件，如满足则保存模型退出，否则执行步骤44)；

44)计算目标函数梯度、照明补偿算子与步长；

45)采用预条件最速下降法计算模型更新量并更新速度模型；

46)重复执行步骤42-46)。

所述的步骤42)具体包括以下步骤：

421)根据9点有限差分格式获取二维频率域声波方程的数值；

422)对于炮点端模拟，震源项则为子波频谱；对于检波点端格林函数，震源项则为1；

423)采用LU分解算法对有限差分方程组求解得到正传波场和格林函数并保存。

所述的步骤43)中，目标函数的表达式为：

其中，ω为圆频率、G₀(r，g)为在背景介质v₀(r)中检波点g激发成像点r处的格林函数，u₀(r，s)为炮点s激发成像点r处的正传波场，d(g，s)为炮检对(g，s)所对应的反射数据，Re表示取复数的实部，v_r为成像点r处的速度。

所述的步骤43)中，目标函数梯度的计算式为：

其中，i为速度模型空间内的梯度更新点，v_i为更新点i的速度，u₀(i，s)为炮点s激发更新点i处的正传波场，G₀(i，r)为成像点r激发更新点i处的格林函数，G₀(i，g)为检波点g激发更新点i处的格林函数，sgn(x)为符号函数。

所述的步骤43)中，退出条件为达到设定的迭代次数上限。

所述的步骤45)中，利用预条件最速下降法计算模型更新量并更新模型的方式为：

其中，H₀为照明补偿算子，λ为防止求逆操作不稳定的阻尼因子，t为标量步长，v^k+1为第k+1次迭代的速度模型，v^k为第k次迭代的速度模型，I为累加像。

所述的照明补偿算子H₀的表达式为：

H₀＝diag{K^TK}

其中，K^T为核函数矩阵。

与现有技术相比，本发明具有以下优点：

一、无需交互拾取，无需模型高波数先验信息，无需反偏移合成反射数据：该发明在理论上建立波场表达的共成像点道集，然后通过一范数定义CIG水平程度目标函数，进而在最优化框架下实现中深层速度模型的自动更新。

二、计算效率高，易于并行：该方法无需扩展域成像，计算量与FWI相同，而且在频率域实现，与现有的波动方程MVA相比，计算时间大大减少。

三、计算稳定，反演精度高：该发明在最优化框架下实现中深层速度模型的自动更新，因此计算稳定，另外，利用波形提取走时信息，同时可以方便地施加预条件，克服现有方法无法实现精确双向照明补偿的难题，因此反演精度高。

附图说明

图1为本发明的方法流程图。

图2为本发明的硬件结构示意图。

图3为本发明目标函数性态分析结果。

图4为实施例1的真实理论模型图。

图5为实施例1的本发明方法反演结果图。

图6为实施例2的真实理论模型图。

图7为实施例2的初始模型图。

图8为实施例2的本发明方法反演结果图。

图中标记说明：

1、数据采集器，2、处理器，3、显示器，4、输入设备。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

如图1所示，本发明提供一种基于共成像点道集的反射全波形反演方法，本发明方法的原理为：

根据Woodward(1992)，某一个炮检对(g，s)对地下r点的像可以表达为

其中，ω为圆频率，v表示速度，Re表示取复数的实部，u₀(r，s)表示炮点s在r处的正传波场，G₀(r，g)表示检波点g引起的r处的格林函数，d(g，s)表示观测到的炮检对(g，s)对应的反射数据。

根据A1-Yahya准则，如果模型是准确的，那么共成像点道集应该是平的，所以，模型的最优估计应该使得如下目标函数取最大值：

进而可以导出其梯度为：

其中，sgn(x)为符号函数，满足：

然后，采用预条件最速下降法按照公式(5)即可实现模型的更新：

其中，λ为防止求逆操作不稳定的阻尼因子，t为标量步长。

为了定量说明发明方法对初始模型的依赖性，进行了目标函数性态分析，结果如图3所示。说明发明方法对初始模型有较弱的依赖性，目标函数凸性较好。

实施例1：

本实施例以二维Inclusion模型作为真实模型(如图4所示)，该模型共有401×201个网格，网格间距为10米×10米，速度最小最大值分别为2000米/秒和2500米/秒，在该模型上进行声波正演模拟，共正演201炮，炮点均匀分布于地下10米深度，炮间距为20米，第一炮在水平位置0米处，检波点以10米间隔均匀分布于地表，每炮的最大偏移距为4000米，最小偏移距为0米，初始模型取2000米/秒的均匀背景模型，反演过程中将2-22Hz的频率域数据以0.2Hz为间隔分为100个频率，反演结果如图5所示。参照真实模型可以发现，本发明提出的基于共成像点道集的反射全波形反演方法可以获得与真实模型相比拟的结果，球形异常体和底部的水平反射层无论在形态上还是数值上都与真实模型比较接近。

具体实施步骤如下：

(1)数据采集器1采集地震波信号，并对原始地震数据进行滤波、道均衡等预处理，并切除直达波；

(2)数据采集器1将步骤(1)预处理后的地震波数据逐道输入到处理器中进行傅里叶变换，以获得每个炮检对的频率域记录；

(3)输入设备4根据地下先验信息，利用层析成像等技术建立初始速度模型，并设定反演频率范围与间隔等参数；

(4)处理器2获得当前观测系统下非重复的炮检点数量与位置；

(5)处理器2计算并存储所有频率的炮点端波场、检波点端格林函数与成像点格林函数；

(6)处理器2进入频率循环，在每个单一频率下执行步骤(7)、(8)；

(7)处理器2根据步骤(5)的波场和格林函数计算单频目标函数、目标函数梯度与照明补偿算子；

(8)处理器2结束频率循环，对所有的频率结果进行叠加获得总的目标函数、目标函数梯度、照明补偿算子与步长；

(9)处理器2根据目标函数值判断是否满足退出条件，如满足则继续执行以下步骤，否则利用梯度与补偿对模型进行更新，并继续执行步骤(5)～(9)；

(10)输出设备3通过显示器显示模型，结束退出；

实施例2：

本实施例以二维Overthrust理论模型作为真实模型(如图6所示)，该模型共有801×187个网格，网格间距为25米×25米，速度最小最大值分别为2300米/秒和6500米/秒，在该模型上进行声波正演模拟，共正演401炮，炮点均匀分布于地下25米深度，炮间距为50米，第一炮在水平位置0米处，检波点以25米间隔均匀分布于地表，双边非规则放炮，每炮的最大偏移距限制为4000米，最小偏移距为0米，初始模型为真实模型高斯平滑结果(如图7所示)，反演过程中将2-22Hz的频率域数据以0.2Hz为间隔分为100个频率，反演结果如图8所示，可见本发明基本准确重建了整个模型的低波数信息。

本发明利用Al-Yahya准则，参照逆时偏移成像(RTM)原理，在理论上建立波场表达的共成像点道集，然后通过一范数定义CIG水平程度目标函数，进而在最优化框架下实现中深层速度模型的自动更新，该发明无需交互拾取，无需模型高波数先验信息，无需反偏移合成反射数据，计算量与FWI相同，计算稳定，反演效果好。

Claims (3)

1.一种基于共成像点道集的反射全波形反演方法，其特征在于，包括以下步骤：

1)对原始地震数据进行预处理；

2)对预处理过的地震数据进行傅里叶变换，得到频率域反射炮记录d(ω，s，g)，其中，ω为圆频率，s为炮点位置，g为检波点位置；

3)采用层析成像方法建立初始速度模型，并设定反演参数；

4)进行频率域声波方程反射全波形反演，获得最终速度剖面，具体包括以下步骤：

41)根据观测系统统计不重复的炮点和检波点的数量及位置；

42)执行频率循环，在每一个频率循环内部，进行炮点端的单频正传波场u₀(s)、检波点端的单频格林函数G₀(g)以及成像点单频格林函数G₀(r)的计算，其中，r为成像点位置，具体包括以下步骤：

421)根据9点有限差分格式获取二维频率域声波方程的数值；

422)对于炮点端模拟，震源项则为子波频谱；对于检波点端格林函数，震源项则为1；

423)采用LU分解算法对有限差分方程组求解得到正传波场和格林函数并保存；

43)计算目标函数值并判断是否满足退出条件，退出条件为达到设定的迭代次数上限，如满足则保存模型退出，否则执行步骤44)，目标函数的表达式为：

其中，ω为圆频率、G₀(r，g)为在背景介质v₀(r)中检波点g激发成像点r处的格林函数，u₀(r，s)为炮点s激发成像点r处的正传波场，d(g，s)为炮检对(g，s)所对应的反射数据，Re表示取复数的实部，v_r为成像点r处的速度；

44)计算目标函数梯度、照明补偿算子与步长，目标函数梯度的计算式为：

其中，i为速度模型空间内的梯度更新点，v_i为更新点i的速度，u₀(i，s)为炮点s激发更新点i处的正传波场，G₀(i，r)为成像点r激发更新点i处的格林函数，G₀(i，g)为检波点g激发更新点i处的格林函数，sgn(x)为符号函数；

45)采用预条件最速下降法计算模型更新量并更新速度模型，利用预条件最速下降法计算模型更新量并更新模型的方式为：

其中，H₀为照明补偿算子，λ为防止求逆操作不稳定的阻尼因子，t为标量步长，v^k+1为第k+1次迭代的速度模型，v^k为第k次迭代的速度模型，I为累加像，所述的照明补偿算子H₀的表达式为：

H₀＝diag{K^TK}

其中，K^T为核函数矩阵；

46)重复执行步骤42-46)。

2.根据权利要求1所述的一种基于共成像点道集的反射全波形反演方法，其特征在于，所述的步骤1)中，预处理包括去噪、滤波、道均衡和去透射波。

3.根据权利要求1所述的一种基于共成像点道集的反射全波形反演方法，其特征在于，所述的步骤3)中，反演参数包括成像点坐标与数量、频率范围和频率间隔。

Images