CN105093278B - 基于激发主能量优化算法的全波形反演梯度算子提取方法 - Google Patents

Abstract

本发明公开了一种基于激发主能量优化算法的全波形反演梯度算子的提取方法，通过建立全波形反演参数模型；采用常规的叠加前单炮数据作为观测数据，建立正反演观测系统；将震源置于地表炮点处，利用优化有限差分算法计算该震源产生的正传波场；确定震源波场的正向传播过程中地下各点震源波场最大值及其对应的时刻，并保存地下各点在激发主能量时窗内的波场；再利用优化有限差分模拟算法计算波场残差逆时传播波场，判断地下各点不同时刻的波场残差逆时传播波场是否位于该点激发主能量时窗之内，通过基于激发主能量优化算法得到全波形反演梯度算子。本发明在保证精度的前提下，大大减少了I/O操作，提升了计算效率。

Description

基于激发主能量优化算法的全波形反演梯度算子提取方法

技术领域

本发明涉及石油地震勘探速度建模技术领域，尤其涉及一种基于激发主能量优化算法的全波形反演梯度算子的提取方法。

背景技术

地震波的传播速度是地震勘探中的一个基本参数，在地震资料处理和地震解释的诸多环节中扮演着重要的角色。影响地震波传播速度的因素有很多，诸如地下介质的岩性、密度、孔隙度、埋藏深度等等，这些因素最终体现在地震波的走时、波形、振幅和相位等属性上变化。只有充分了解速度及影响因素的相互关系，才有可能更为准确地估计出地下的结构和构造。我国大多数勘探区域越来越复杂，呈现出地表构造和地下构造双复杂特征。传统的速度分析方法已经很难达到所需的精度要求，也会直接影响成像质量。常用的速度分析方法有：叠加速度分析、偏移速度分析、层析速度反演和全波形反演。这四种速度分析方法的准确性依次升高，复杂度也依次升高。

全波形反演技术是目前精度很高的速度分析和反演工具，其利用叠前地震记录与正演模拟炮记录之间的寻求最佳匹配来获取最优反演速度场。在反演过程中不仅考虑了地震波传播的旅行时等运动学信息，而且加入了振幅、相位、波形等动力学信息，能够适应强横向变速介质和各向异性介质的速度反演。

全波形反演的发展大致可以分为三个阶段，第一个阶段是从20世纪80年代到90年代，这一时期是全波形反演理论形成和初步验证阶段，Lailly、Tarantola等人建立了基于波动方程的波形反演的理论构架，之后一批学者分别利用简单模型对反射和透射波形反演效果进行了理论验证。Gauthier等人率先对理论模型合成的多偏移距数据进行了声波和弹性波全波形反演测试，反演结果既展示了全波形反演方法的应用潜力，也暴露出其根本缺陷——对初始模型的严重依赖性。

第二个阶段是从上世纪90年代初到21世纪初，为全波形稳定发展阶段，全波形反演方法形成了多个分支和变种，其中以Pratt等提出的频率域波形反演方法最具代表性，Pratt提出频率域声波和弹性波波形反演理论，并利用实际井间地震数据对方法进行验证。该方法利用稀疏矩阵直接解法求解频率域的离散化声波或弹性波波动方程，只要频率与物性参数不改变，通过一次稀疏矩阵LU分解和多次高斯消元法中的前向、后向回代运算，即可完成所有炮点的单频波场响应计算；应用研究方面以井间地震波形反演和海洋天然气水合物的叠前波形反演最具影响力。此外为解决全波形反演因初始模型不准而陷入局部极小问题，Zhou等人先后给出了波动方程走时和波形联合反演方法，以充分利用走时和波形两类数据对不同尺度模型参数变化的敏感度的优势，使井间数据全波形反演更为稳定。

自21世纪初至今，全波形反演达到了蓬勃发展的阶段，全波形反演逐步从二维反演走向三维反演，从声波走向弹性波反演，从单参数反演到多参数反演，从各向同性到各向异性反演。此外还有Laplace/Laplace-Fourier域全波形反演、各向异性全波形反演、先验信息约束的全波形反演和基于超级炮和相位编码的全波形反演。此外随着计算机水平的进步，全波形反演的计算效率进一步得到了提升，其实现平台逐步从单机单核走向了多核多节点并行平台、GPU并行平台。在一定程度上，全波形反演的计算效率完全取决于大数据的I/O交换。

发明内容

为了提高时间域全波形反演的计算效率，本发明提供了一种基于激发主能量优化算法的全波形反演梯度算子的提取方法，其利用有限差分算法计算震源波场和波场误差反传波场，再根据激发主能量优化算法计算模型更新量的梯度。该方法可以在保障计算精度的同时，减少大量I/O操作，大大提高全波形反演的计算效率，相比于常规梯度算法，在保持精度的前提下，提高了近3倍的效率。

本发明解决其技术问题所采用的技术方案是：

一种基于激发主能量优化算法的全波形反演梯度算子的提取方法，包括：

建立全波形反演参数模型；

利用常规的叠前单炮数据作为观测数据，建立正反演观测系统；

将震源置于地表炮点处，利用优化有限差分算法计算该震源产生的正传波场；

确定震源波场的正向传播过程中地下各点震源波场最大值及其对应的时刻，并保存地下各点在激发主能量时窗内的波场；

利用优化有限差分模拟算法计算波场残差逆时传播波场，判断地下各点不同时刻的波场残差逆时传播波场是否位于该点激发主能量时窗之内，通过基于激发主能量优化算法得到全波形反演梯度算子。

前述的基于激发主能量优化算法的全波形反演梯度算子的提取方法，具体包括如下步骤：

(1)建立全波形反演参数模型；

(2)采用常规的叠加前单炮数据作为观测数据，建立正反演观测系统；

(3)将震源置于地表炮点处，利用下式表示的优化有限差分算法计算该震源产生的正传波场：

其中U为震源产生的正传波场，c_m为优化有限差分系数，v为介质速度，Δx和Δz为空间采样间隔，Δt为时间采样间隔，i,k为x和z空间坐标，n为时间坐标；

(4)确定震源波场的正向传播过程中地下各点最大值及其对应的时刻Ts；

(5)定义激发主能量时窗长度为正传波场中一个子波的时间长度w，保存地下各点在激发主能量时窗内的波场，即以地下各点的最大值对应时刻Ts为时窗中心的整个主能量时窗之内的波场值；

(6)利用步骤(3)中所述优化有限差分正演模拟算法得到叠前单炮合成数据RU，并其与观测地震数据d_obs求差得到波场残差RU-d_obs；

(7)利用下式表示的优化有限差分模拟算法计算波场残差逆时传播波场：

其中RU-d_obs为波场残差，L^*为反传算子，c_m为优化有限差分系数，U_b为波场残差逆时传播波场，v为介质速度，Δx和Δz为空间采样间隔，Δt为时间采样间隔，i,k为x和z空间坐标，n为时间坐标；

(8)判断地下任意一点的不同时刻波场残差的逆时传播波场是否位于该点激发主能量时窗之内，如果位于，则按下式提取基于激发主能量优化算法的用于全波形反演的梯度算子：

其中RU-d_obs为波场残差，L^*为反传算子，为正演波场对时间的二阶偏导，v为空间速度分布，Ts为激发最大振幅对应时间步，w为激发主能量时窗长度，g(v)为全波形反演的梯度。

上述方案更进一步包括：

步骤(1)中全波形反演的模型参数求取采用迭代方法求解，迭代公式为

m^k+1＝m^k+Δm (1)

其中m^k+1，m^k分别表示第k次迭代模型参数和第k+1次迭代模型参数，Δm是参数模型的更新量，其表达式为：

Δm＝-αg (2)

其中g是参数模型更新量的梯度算子；α是计算步长。所以利用先验信息建立全波形反演参数的初始模型；

步骤(3)中将震源置于地表炮点处，采用优化有限差分算法计算震源激发的正传波场；

正传波场U满足如下的波动方程：

其中U(x,z,t)为震源激发的正传波场；v(x,z)为介质速度；

利用Taylor展开进行时间和空间差分，将时间离散形式：

和空间差分的形式

带入波动方程，得到各向同性介质正传波场数值计算公式

其中c_m为优化有限差分系数，v为介质速度，Δx和Δz为空间采样间隔，Δt为时间采样间隔，i,k为x和z空间坐标，n为时间坐标；

步骤(4)是在震源波场的正向传播过程中，根据正传波场U计算地下各点的最大值及其对应的时刻Ts；

步骤(7)中利用优化有限差分模拟算法计算波场残差逆时传播波场，其中逆时传播波场满足如下方程：

其中U_b(x,z,t)为波场误差的逆时传播波场；v(x,z)为介质速度，

利用Taylor展开进行时间和空间差分，即将时间离散形式：

和空间差分的形式

带入逆时传播波场满足的方程，则得出波场误差的逆时传播波场数值计算公式：

其中RU-d_obs为波场残差，L^*为反传算子，c_m为优化有限差分系数，U_b为波场残差逆时传播波场，v为介质速度，Δx和Δz为空间采样间隔，Δt为时间采样间隔，i,k为x和z空间坐标，n为时间坐标。

本发明的有益效果是：传统的全波形算法在构造梯度时需要保存地震震源波场全部时刻波场信息，基于激发主能量优化梯度算法定义了构成梯度的主能量，只需保存构成激发主能量时窗的正传波场，并与相应反传波场计算即可得到模型更新量的梯度，大大减少了I/O操作，在保持精度的前提下，可以提高3倍的效率。

附图说明

图1是基于激发主能量优化算法的全波形反演梯度算子提取方法的工作流程图。

图2是真实参数模型图。

图3是初始参数模型图。

图4是模型各点激发最大振幅图。

图5是模型各点激发最大振幅所对应时刻。

图6是模型数据常规梯度算法计算的梯度。

图7是模型数据激发主能量优化算法计算的梯度。

图8是实际资料常规梯度算法计算的梯度。

图9是实际资料激发主能量优化算法计算的梯度。

图10是实际资料常规梯度算法反演结果。

图11是实际资料激发主能量优化算法反演结果。

图12是实际资料单道反演结果对比图。

图13是实际资料反演效率对比图。

具体实施方式

为使本发明的目的、特征和优点能更明显易懂，下文特举出较佳实施例，并配合所附图式，作详细说明如下：

实施例1

结合附图1，一种基于激发主能量优化算法的全波形反演梯度算子的提取方法，包括如下步骤：

建立全波形反演参数模型；

采用常规的叠前单炮数据作为观测数据，建立正反演观测系统；

将震源置于地表炮点处，利用优化有限差分算法计算该震源产生的正传波场；

确定震源波场的正向传播过程中地下各点震源波场最大值及其对应的时刻，并保存地下各点在激发主能量时窗内的波场；

利用优化有限差分模拟算法计算波场残差逆时传播波场，判断地下各点不同时刻的波场残差逆时传播波场是否位于该点激发主能量时窗之内，通过基于激发主能量优化算法得到全波形反演梯度算子。

具体步骤为：

步骤1：全波形反演的模型参数求取采用迭代方法求解，迭代公式为

m^k+1＝m^k+Δm (1)

其中m^k+1，m^k分别表示第k次迭代模型参数和第k+1次迭代模型参数，Δm是参数模型的更新量，其表达式为：

Δm＝-αg (2)

其中g是参数模型更新量的梯度算子；α是计算步长。所以利用先验信息建立全波形反演参数的初始模型。

步骤2：因为本发明是针对叠前单炮数据，所以采用常规的叠前单炮数据作为观测数据，建立正反演观测系统。

步骤3：将震源置于地表炮点处，计算正传波场U。正传波场U满足如下的波动方程：

其中U(x,z,t)为震源激发的正传波场；v(x,z)为介质速度。

本发明采用优化有限差分算法计算震源激发的正传波场。利用Taylor展开进行时间和空间差分，将时间离散形式：

和空间差分的形式

带入波动方程，得到各向同性介质正传波场数值计算公式

其中c_m为优化有限差分系数，v为介质速度，Δx和Δz为空间采样间隔，Δt为时间采样间隔，i,k为x和z空间坐标，n为时间坐标；

步骤4：在震源波场的正向传播过程中，根据正传波场U计算地下各点的最大值及其对应的时刻Ts。

步骤5：定义激发主能量时窗长度为正传波场中一个子波的时间长度w，保存以地下各点的最大值对应时刻Ts为时窗中心的整个主能量时窗之内的波场值，即保存地下各点在激发主能量时窗内的波场，为计算梯度奠定基础。

步骤6：利用优化有限差分正演模拟算法得到叠前单炮合成数据RU，并将其与观测地震数据d_obs求差得到波场残差RU-d_obs，

步骤7：利用优化有限差分模拟算法计算波场残差逆时传播波场。逆时传播波场满足如下方程：

其中U_b(x,z,t)为波场误差的逆时传播波场；v(x,z)为介质速度。

本发明采用优化有限差分算法计算波场误差的逆时传播波场。利用Taylor展开进行时间和空间差分，即将时间离散形式：

和空间差分的形式

带入逆时传播波场满足的方程，则可以得出波场误差的逆时传播波场数值计算公式：

其中RU-d_obs为波场残差，L^*为反传算子，c_m为优化有限差分系数，U_b为波场残差逆时传播波场，v为介质速度，Δx和Δz为空间采样间隔，Δt为时间采样间隔，i,k为x和z空间坐标，n为时间坐标。

步骤8：因为常规的时间域内基于时间相关条件的梯度算子表达式为：

其中RU-d_obs为波场残差，L^*为反传算子，为正演波场对时间的二阶偏导，v为空间速度分布。

而本发明的激发主能量优化算法是在激发主能量时窗内利用正传波场和波场误差逆时波场来计算梯度，其表达式为：

其中RU-d_obs为波场残差，L^*为反传算子，为正演波场对时间的二阶偏导，v为空间速度分布，Ts为激发最大振幅对应时间步，w为激发主能量时窗长度。

所以在步骤7之后，判断地下各点不同时刻的波场残差逆时传播波场是否位于该点激发主能量时窗之内，如果位于，则按公式(12)提取基于激发主能量优化算法的用于全波形反演的梯度算子。

激发主能量优化梯度算法在波场的正向传播过程中只需保存构成激发主能量时窗的正传波场，并与相应反传波场计算即可得到模型更新量的梯度，大大减少了I/O操作，在保持精度的前提下，可以提高3倍的效率。

应用实验例1

本发明基于激发主能量优化算法计算的梯度算子的提取方法，应用于模型数据，取得了理想的计算效果。图2为真实参数模型；图3为初始参数模型；首先将震源置于地表炮点处，利用优化有限差分算法计算该震源产生的正传波场；确定震源波场的正向传播过程中地下各点震源波场最大值(图4)及其对应的时刻(图5)，并保存地下各点在激发主能量时窗内的波场；再利用优化有限差分模拟算法计算波场残差逆时传播波场，判断地下各点不同时刻的波场残差逆时传播波场是否位于该点激发主能量时窗之内，通过基于激发主能量优化算法得到全波形反演梯度算子(图7)。与常规算法计算梯度结果(图6)相比，基于激发主能量优化算法计算的梯度结果精度与其相当，但是计算速度更快，计算效率大大提高。

应用实验例2

应用本发明基于激发主能量优化算法计算的梯度算子的提取方法，应用于油田某块区的地震数据，取得了精确的反演效果和明显的效率提高。图8为常规梯度算法计算的梯度，图9为基于激发主能量优化算法计算的梯度，二者计算精度相当。图10为常规梯度算法的反演结果，图11为基于激发主能量优化算法的反演结果，图12为单道反演结果对比图，反演结果表面主能量梯度算法和常规梯度算法误差在5％以下，即在全波形反演过程中，基于激发主能量优化算法的梯度提取方法满足精度的需求。图13为反演效率对比图。常规梯度算法进行一次迭代的时间为320min，激发主能量优化梯度算法迭代一次的时间为90min。总结，在多炮的全波形反演算法中，在相同的参数情况下，激发主能量优化梯度算法相比于常规梯度算法，在保持精度的前提下，提高了近3倍的效率。

Claims (3)

1.一种基于激发主能量优化算法的全波形反演梯度算子提取方法，包括如下步骤：

建立全波形反演参数模型；

利用常规的叠前单炮数据作为观测数据，建立正反演观测系统；

将震源置于地表炮点处，利用优化有限差分正演模拟算法计算该震源产生的正传波场；

确定震源波场的正向传播过程中地下各点震源波场最大值及其对应的时刻，并保存地下各点在激发主能量时窗内的波场；

利用优化有限差分模拟算法计算波场残差逆时传播波场，判断地下各点不同时刻的波场残差逆时传播波场是否位于该点激发主能量时窗之内，通过基于激发主能量优化算法得到全波形反演梯度算子。

2.依据权利要求1所述的基于激发主能量优化算法的全波形反演梯度算子提取方法，具体包括如下步骤：

(1)建立全波形反演参数模型；

(2)采用常规的叠加前单炮数据作为观测数据，建立正反演观测系统；

(3)将震源置于地表炮点处，利用下式表示的优化有限差分正演模拟算法计算该震源产生的正传波场：

其中U为震源产生的正传波场，c_m为优化有限差分系数，v为介质速度，Δx和Δz为空间采样间隔，Δt为时间采样间隔，i,k为x和z空间坐标，n为时间坐标；

(4)确定震源波场的正向传播过程中地下各点最大值及其对应的时刻Ts；

(5)定义激发主能量时窗长度为正传波场中一个子波的时间长度w，保存地下各点在激发主能量时窗内的波场，即以地下各点的最大值对应时刻Ts为时窗中心的整个主能量时窗之内的波场值；

(6)利用步骤(3)中所述优化有限差分正演模拟算法得到叠前单炮合成数据RU，并其与观测地震数据d_obs求差得到波场残差RU-d_obs；

(7)利用下式表示的优化有限差分模拟算法计算波场残差逆时传播波场：

其中RU-d_obs为波场残差，L^*为反传算子，c_m为优化有限差分系数，U_b为波场残差逆时传播波场，v为介质速度，Δx和Δz为空间采样间隔，Δt为时间采样间隔，i,k为x和z空间坐标，n为时间坐标；

(8)判断地下任意一点的不同时刻波场残差的逆时传播波场是否位于该点激发主能量时窗之内，如果位于，则按下式提取基于激发主能量优化算法的用于全波形反演的梯度算子：

其中RU-d_obs为波场残差，L^*为反传算子，为正演波场对时间的二阶偏导，v为空间速度分布，Ts为激发最大振幅对应时间步，w为激发主能量时窗长度，g(v)为全波形反演的梯度。

3.依据权利要求2所述的基于激发主能量优化算法的全波形反演梯度算子提取方法，更进一步包括：

步骤(1)中全波形反演的模型参数求取采用迭代方法求解，迭代公式为

m^k+1＝m^k+Δm (1)

其中m^k+1，m^k分别表示第k次迭代模型参数和第k+1次迭代模型参数，Δm是参数模型的更新量，其表达式为：

Δm＝-αg (2)

其中g是参数模型更新量的梯度算子；α是计算步长；所以利用先验信息建立全波形反演参数的初始模型；

步骤(3)中将震源置于地表炮点处，采用优化有限差分正演模拟算法计算震源激发的正传波场；

正传波场U满足如下的波动方程：

其中U(x,z,t)为震源激发的正传波场；v(x,z)为介质速度；

利用Taylor展开进行时间和空间差分，将时间离散形式：

和空间差分的形式

带入波动方程，得到各向同性介质正传波场数值计算公式

其中c_m为优化有限差分系数，v为介质速度，Δx和Δz为空间采样间隔，Δt为时间采样间隔，i,k为x和z空间坐标，n为时间坐标；

步骤(4)是在震源波场的正向传播过程中，根据正传波场U计算地下各点的最大值及其对应的时刻Ts；

步骤(7)中利用优化有限差分模拟算法计算波场残差逆时传播波场，其中逆时传播波场满足如下方程：

其中U_b(x,z,t)为波场误差的逆时传播波场；v(x,z)为介质速度，

利用Taylor展开进行时间和空间差分，即将时间离散形式：

和空间差分的形式

带入逆时传播波场满足的方程，则得出波场误差的逆时传播波场数值计算公式：

其中RU-d_obs为波场残差，L^*为反传算子，c_m为优化有限差分系数，U_b为波场残差逆时传播波场，v为介质速度，Δx和Δz为空间采样间隔，Δt为时间采样间隔，i,k为x和z空间坐标，n为时间坐标。