CN112130199A - 一种优化的最小二乘逆时偏移成像方法 - Google Patents

Abstract

本发明公开了一种优化的最小二乘逆时偏移成像方法，具体按照以下步骤实施：步骤1、采用波动方程对地震数据进行偏移与反偏移；步骤2、采用共轭梯度法构建最小二乘框架，并引入优化因子，对共轭梯度法的收敛速度进行提升；步骤3、采用GPU的CUDA计算架构进一步提升优化的LSRTM成像的计算效率。采用优化的最小二乘逆时偏移成像方法能够在保持相同成像精度下提高算法的综合计算效率，综合计算效率平均提升20％，有效验证了优化的最小二乘逆时偏移成像能够较好的提高传统最小二乘逆时偏移成像方法的计算效率。

Description

一种优化的最小二乘逆时偏移成像方法

技术领域

本发明属于信号处理技术领域，涉及一种优化的最小二乘逆时偏移成像方法。

背景技术

地震勘探的核心目标是尽可能定量地、精确地描述油气藏，地震波成像由定位反射(散射)点位置发展到当前的估计(角度)反射系数是地震勘探的核心需求。随着实际需求的不断提高和大量资源开采殆尽的现状，对复杂构造成像需求急剧增加，如何更精准更快速的成像成为研究热点。然而，复杂地质构造具有无规律起伏地表、复杂地下构造、处理数据量大等特点，这给偏移成像带来了巨大的挑战。

基于反演理论的最小二乘逆时偏移(Least-Square Reverse Time Migration，LSRTM)将成像问题当作一个反问题来处理，通过比较由偏移剖面所产生的合成数据与实际采集数据之间的相关性来判定成像结果是否准确，通过多次自动修正成像结果来提升相关性，以寻求更接近于真实的地下反射系数，从而更好地进行岩性储层成像和储层参数反演。该方法是实现地震成像理论由常规地下岩石的几何结构描述向保幅成像的推进和发展，具有更高的成像精度，也是实现高精度储层参数反演的关键。最小二乘逆时偏移运算成本非常高，除了因为双程波方程波场传播算子计算量大以外，另一个更重要的原因是迭代收敛速度慢。LSRTM通过迭代函数使误差函数达到最小，得到数据空间或者模型空间求解，模型空间方法不需要显式求解Hessian 矩阵，而数据空间方法需要显式求解Hessian矩阵，然后对偏移结果应用 Hessian矩阵的逆可得到反演结果。然而，LSRTM成像在进行1次迭代时包含1次逆时偏移成像和1次残差的反偏移过程，即1次LSRTM成像的时间为常规RTM成像时间的2.5倍，而且随着迭代次数的增加计算量将成倍增加，这一问题严重限制LSRTM成像方法在实际地震数据处理中的应用。

在对LSRTM成像方法的效率提升方面，Cerdan等提出用LSRTM方法计算迭代求解低秩最小二乘问题的稀疏预条件，其更新正则化问题的不完全分解，用以得到原始解，有效的加速了迭代收敛的速度；Yang提出限制背景波场传播时间至停止时间而不是最大记录时间，并给出可以选取最大记录时间的一半作为停止时间，从而减少计算量提高计算效率；Guo等提出一种基于一阶速度-应力波动方程的最小二乘逆时偏移成像方法，其通过线性化波动方程为速度-应力方程，导出相应的速度扰动方程和伴随方程，从而更新梯度公式改善成像质量和收敛速度。

发明内容

本发明的目的是提供一种优化的最小二乘逆时偏移成像方法，具有通过优化因子加快共轭梯度法的收敛速度，从而提高最小二乘逆时偏移成像算法效率的特点。

本发明所采用的技术方案是，一种优化的最小二乘逆时偏移成像方法，具体按照以下步骤实施：

步骤1、采用波动方程对地震数据进行偏移与反偏移；

步骤2、采用共轭梯度法构建最小二乘框架，并引入优化因子，对共轭梯度法的收敛速度进行提升；

步骤3、采用GPU的CUDA计算架构进一步提升优化的LSRTM成像的计算效率。

步骤1具体按照以下步骤实施：

步骤1.1、利用常密度声波方程，获得波场递推结果；

步骤1.2、用扰动法对速度模型进行表达，得到反偏移方程；

步骤1.3、对炮集d(x_s,x_g,t)的偏移通过两次有限差分的模拟，实现偏移和反偏移。

步骤1.1具体为：

通过式

其中，v为速度场，f为震源项，u为全波场；

对该式进行时间2阶空间12阶差分变换，获得波场递推结果，如公式 (2)所示：

步骤1.2具体为：

通过式

式中，v₀ ²为假定的背景(或参考)速度模型；α(x)表示相对于背景速度模型的速度不均匀分布，也称为速度扰动函数；

同样的，全波场也可分为背景波场u₀和扰动波场δu：

u(x_s,x,t)＝u₀(x_s,x,t)+δu(x_s,x,t) (4)

把式(3)和式(4)带入式(1)，应用Born近似可得扰动波场δu的方程：

方程(5)为波场传播是线性传播，是由入射波场与速度不均匀相互作用产生的散射波场，也就是反偏移方程。

步骤1.3具体为：

一次是与炮点有关的正传波场，一次与是检波点有关的反传波场，如下所示：

其中，q(x_s,x,t)是与检波点有关的反传波场；

用矢量矩阵符号写成算子形式，Born正演的过程可描述为：

d＝Lm (8)

其中，d代表了反偏移数据，

是反射系数模型，L表示Born 正演算子，相应的，偏移过程可描述为：

m_mig＝L^*d (9)

其中，m_mig是偏移成像结果，L*是偏移算子，从而实现偏移和反偏移。

步骤2具体为：

步骤2.1、利用共轭梯度法，求解修正因子β；

步骤2.2、引入优化因子τ_k进行对参数β^PRP进行改进，在保证其成像精度的条件下加快收敛速度。

步骤2.1具体为：

搜索方向是负梯度方向和上一次搜索方向的一个组合，其公式如下：

其中g_k为梯度，α为搜索步长，z^k为搜索方向，β为搜索方向z^k的修正因子；不同的共轭梯度法有不同的β；PRP(Polak-Ribière-Polyak)共轭梯度法具有较好的数值表现，其参数β的求取如公式(11)所示：

该参数β影响最小二乘的收敛速度。

步骤2.2具体为：

优化的公式如下：

式中：α为计算步长；β为共轭梯度修正因子；z为共轭梯度；k为迭代次数；τ_k为调节因子，可以自动调节参数β，使其满足|β|≤β^PRP，优化后的方法不仅保持了PRP共轭梯度法的性能，而且具有良好收敛性能和数值表现。

步骤3具体为：

最小二乘逆时偏移算法需进行迭代计算，计算量随迭代次数线性增加，采取GPU的并行策略为迭代计算时，先对计算任务进行划分，对每个任务独立地进行波场模拟并行计算，计算时每个线程对应波场网格的一个点，每个线程计算一个波场值即每个线程进行一次差分方程计算。

本发明的有益效果是：本发明提出了一种优化的最小二乘逆时偏移成像方法，采用优化的最小二乘逆时偏移成像方法能够在保持相同成像精度下提高算法的综合计算效率，综合计算效率平均提升20％，有效验证了优化的最小二乘逆时偏移成像能够较好的提高传统最小二乘逆时偏移成像方法的计算效率。

附图说明

图1是本发明优化的最小二乘逆时偏移成像方法流程图；

图2是波场的GPU多线程计算流程图；

图3是盐体速度模型；

图4是优化的最小二乘逆时逆时偏移成像结果；

图5是改进的PRP共轭梯度法与传统PRP共轭梯度法收敛速度的对比图。

具体实施方式

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

本发明一种优化的最小二乘逆时偏移成像方法，如图1所示，具体按照以下步骤实施：

步骤1、采用波动方程对地震数据进行偏移与反偏移；

具体按照以下步骤实施：

步骤1.1、利用常密度声波方程，获得波场递推结果；

步骤1.2、用扰动法对速度模型进行表达，得到反偏移方程；

步骤1.3、对炮集d(x_s,x_g,t)的偏移通过两次有限差分的模拟，实现偏移和反偏移。

步骤1.1具体为：

通过式

其中，v为速度场，f为震源项，u为全波场；

对该式进行时间2阶空间2N阶差分变换，获得波场递推结果，如公式 (2)所示：

步骤1.2具体为：

通过式

式中，v₀ ²为假定的背景(或参考)速度模型；α(x)表示相对于背景速度模型的速度不均匀分布，也称为速度扰动函数；

同样的，全波场也可分为背景波场u₀和扰动波场δu：

u(x_s,x,t)＝u₀(x_s,x,t)+δu(x_s,x,t) (4)

把式(3)和式(4)带入式(1)，应用Born近似可得扰动波场δu的方程：

方程(5)为波场传播是线性传播，是由入射波场与速度不均匀相互作用产生的散射波场，也就是反偏移方程。

步骤1.3具体为：

一次是与炮点有关的正传波场，一次与是检波点有关的反传波场，如下所示：

其中，q(x_s,x,t)是与检波点有关的反传波场；

用矢量矩阵符号写成算子形式，Born正演的过程可描述为：

d＝Lm (8)

其中，d代表了反偏移数据，

是反射系数模型，L表示Born 正演算子，相应的，偏移过程可描述为：

m_mig＝L^*d (9)

其中，m_mig是偏移成像结果，L*是偏移算子，从而实现偏移和反偏移。

步骤2、采用共轭梯度法构建最小二乘框架，并引入优化因子，对共轭梯度法的收敛速度进行提升。

步骤2.1、利用共轭梯度法，求解修正因子β；

步骤2.2、引入优化因子τ_k进行对参数β^PRP进行改进，在保证其成像精度的条件下加快收敛速度。

步骤2.1具体为：

搜索方向是负梯度方向和上一次搜索方向的一个组合，其公式如下：

其中g_k为梯度，α为搜索步长，z^k为搜索方向，β为搜索方向z^k的修正因子；不同的共轭梯度法有不同的β；PRP(Polak-Ribière-Polyak)共轭梯度法具有较好的数值表现，其参数β的求取如公式(11)所示：

该参数β影响最小二乘的收敛速度。

步骤2.2具体为：

引入优化因子τ_k进行对参数β^PRP进行改进，在保证其成像精度的条件下加快收敛速度，从而节约计算成本。优化的公式如下：

式中：α为计算步长；β为共轭梯度修正因子；z为共轭梯度；k为迭代次数；τ_k为调节因子，可以自动调节参数β，使其满足|β|≤β^PRP，优化后的方法不仅保持了PRP共轭梯度法的性能，而且具有良好收敛性能和数值表现。

步骤3、采用GPU的CUDA计算架构进一步提升优化的LSRTM成像的计算效率；

步骤3具体为：地震波逆时偏移波场值的计算是沿着时间线进行的，同一时刻网格每个点计算时，各点之间没有相互联系，因此GPU加速技术适合应用在逆时偏移成像。最小二乘逆时偏移算法需进行迭代计算，计算量随迭代次数线性增加，采取GPU的并行策略为迭代计算时，先对计算任务进行划分，对每个任务独立地进行波场模拟并行计算；这样节约了存储空间，避免一次传输过多数据超出显存导致计算无法进行，降低了设备的硬件要求。

计算时每个线程对应波场网格的一个点，每个线程计算一个波场值即每个线程进行一次差分方程计算。波场的多线程计算流程如图2所示。

优化的LSRTM成像方法对综合计算效率的影响为：

本申请选用盐体数据对算法进行测试。图3为盐体速度模型，图4为优化的LSRTM成像结果，图5为传统共轭梯度法和优化的共轭梯度法收敛速度的对比图。如图4所示，优化的LSRTM能够对盐体精确成像；由图5可见采用优化的LSRTM采用的改进的PRP公认梯度法的收敛速度要由于传统的PRP共轭梯度法的收敛速度，其在迭代48次时已经收敛，而传统的PRP方法在55次时才收敛，故改进的PRP方法速度要由于传统的PRP，从而使得优化的LSRTM成像方法速度要由于传统的LSRTM成像方法。

Claims (9)

1.一种优化的最小二乘逆时偏移成像方法，其特征在于，具体按照以下步骤实施：

步骤1、采用波动方程对地震数据进行偏移与反偏移；

步骤2、采用共轭梯度法构建最小二乘框架，并引入优化因子，对共轭梯度法的收敛速度进行提升；

步骤3、采用GPU的CUDA计算架构进一步提升优化的LSRTM成像的计算效率。

2.根据权利要求1所述的一种优化的最小二乘逆时偏移成像方法，其特征在于，所述步骤1具体按照以下步骤实施：

步骤1.1、利用常密度声波方程，获得波场递推结果；

步骤1.2、用扰动法对速度模型进行表达，得到反偏移方程；

步骤1.3、对炮集d(x_s,x_g,t)的偏移通过两次有限差分的模拟，实现偏移和反偏移。

3.根据权利要求1所述的一种优化的最小二乘逆时偏移成像方法，其特征在于，所述步骤1.1具体为：

通过式

其中，v为速度场，f为震源项，u为全波场；

对该式进行时间2阶空间2N阶差分变换，获得波场递推结果，如公式(2)所示：

4.根据权利要求1所述的一种优化的最小二乘逆时偏移成像方法，其特征在于，所述步骤1.2具体为：

通过式

式中，v₀ ²为假定的背景(或参考)速度模型；α(x)表示相对于背景速度模型的速度不均匀分布，也称为速度扰动函数；

同样的，全波场也可分为背景波场u₀和扰动波场δu：

u(x_s,x,t)＝u₀(x_s,x,t)+δu(x_s,x,t) (4)

把式(3)和式(4)带入式(1)，应用Born近似可得扰动波场δu的方程：

方程(5)为波场传播是线性传播，是由入射波场与速度不均匀相互作用产生的散射波场，也就是反偏移方程。

5.根据权利要求1所述的一种优化的最小二乘逆时偏移成像方法，其特征在于，所述步骤1.3具体为：

一次是与炮点有关的正传波场，一次与是检波点有关的反传波场，如下所示：

其中，q(x_s,x,t)是与检波点有关的反传波场；

用矢量矩阵符号写成算子形式，Born正演的过程可描述为：

d＝Lm (8)

其中，d代表了反偏移数据，

是反射系数模型，L表示Born正演算子，相应的，偏移过程可描述为：

m_mig＝L^*d (9)

其中，m_mig是偏移成像结果，L*是偏移算子，从而实现偏移和反偏移。

6.根据权利要求1所述的一种优化的最小二乘逆时偏移成像方法，其特征在于，所述步骤2具体为：

步骤2.1、利用共轭梯度法，求解修正因子β；

步骤2.2、引入优化因子τ_k进行对参数β^PRP进行改进，在保证其成像精度的条件下加快收敛速度。

7.根据权利要求1所述的一种优化的最小二乘逆时偏移成像方法，其特征在于，所述步骤2.1具体为：

搜索方向是负梯度方向和上一次搜索方向的一个组合，其公式如下：

其中g_k为梯度，α为搜索步长，z^k为搜索方向，β为搜索方向z^k的修正因子；不同的共轭梯度法有不同的β；PRP(Polak-Ribière-Polyak)共轭梯度法具有较好的数值表现，其参数β的求取如公式(11)所示：

该参数β影响最小二乘的收敛速度。

8.根据权利要求1所述的一种优化的最小二乘逆时偏移成像方法，其特征在于，所述步骤2.2具体为：

优化的公式如下：

式中：α为计算步长；β为共轭梯度修正因子；z为共轭梯度；k为迭代次数；τ_k为调节因子，可以自动调节参数β，使其满足|β|≤β^PRP，优化后的方法不仅保持了PRP共轭梯度法的性能，而且具有良好收敛性能和数值表现。

9.根据权利要求1所述的一种优化的最小二乘逆时偏移成像方法，其特征在于，所述步骤3具体为：

最小二乘逆时偏移算法需进行迭代计算，计算量随迭代次数线性增加，采取GPU的并行策略为迭代计算时，先对计算任务进行划分，对每个任务独立地进行波场模拟并行计算，计算时每个线程对应波场网格的一个点，每个线程计算一个波场值即每个线程进行一次差分方程计算。

Images