一种基于马尔可夫决策过程的反射层位追踪方法及装置

技术领域

本发明涉及地球物理勘探技术领域，特别是一种基于马尔可夫决策过程的反射层位追踪方法及装置。

背景技术

在勘探地震中，成像剖面上的反射结构信息对速度反演和地震地质解释都具有重要的作用。在速度反演过程中，合理地应用反射结构信息，能够有效地提高速度反演的收敛性，约束反演解，从而获得具有地质意义的速度模型。在地震地质解释中，依据成像剖面上的明确的反射特征信息可以确定地质层位的相互关系，提高地震地质解释的可靠性和精度，有利于建立准确的储层模型、提高储层预测精度。因此，从复杂多变的成像剖面中自动地、高效地提取期望的反射结构信息是非常重要的科学问题。目前索引到的成像剖面上的自动拾取或追踪同相轴的方法主要有以下几种方式：1.主要是基于地震剖面上的各种属性信息，根据同相轴横向上的相似性，采用各种算法(信号和图像处理类方法)，识别并提取同相轴的位置。比如：采用振幅信息(Howard，1991)，相位信息 (Zeng等，1998；Stark，2003；Wu和Zhong，2012)，波形信息(Marfurt等， 1998；Lacaze和Valding，2009；Zou等，2013)，倾角/斜率/结构信息(Bakker， 2002；Lomask等，2006；Parks，2010；Wu，2017)或者采用多属性信息共同约束(Lou和Zhang，2018；Li等2020；Zhang等，2020)；2.采用智能化的学习算法，主要包括有监督学习的人工神经网络方法(Harrigan等，1992；Huang 等，2005；Peters等，2019)以及卷积神经网络(Waldeland等，2018；Shi等， 2019；Tschannen等，2020；Zhang等，2021)；无监督学习的自组织网络方法 (Huang等，1990；陆文凯和牟永光，1998)和卷积自编码网络方法(Shi等， 2020)。

然而常规的基于属性信息的拾取方法严重依赖于属性，属性的好坏制约了方法的有效。当横向上属性变化(子波形态特征)比较剧烈时，这类方法很难追踪到正确的同相轴位置。采用机器学习类的方法，借助于大数据分析，可以根据数据的统计特征，学习成像剖面上的结构信息具有广大的应用前景。然而有监督网络算法，需要大量的训练数据和标签，训练过程需要消耗大量的时间，无监督网络算法存在没有明确的物理含义，结果评价不可靠等问题。因此有必要探索更合适、高精度的反射结构拾取方法。

若把自动地同相轴拾取问题看成一个最优决策过程，则可以通过模仿人类的行为，利用计算机算法代替人工拾取。在机器学习算法中，除了有监督学习算法和无监督学习算法，还有一类学习算法。该类学习算法通过与环境的交互来学习，不断的试错和收益进行决断。这类算法叫做强化学习，其核心算法依赖于马尔可夫决策过程实现。本方法将马尔可夫决策过程(MDP)引入到同相轴拾取过程中。通过拾取目标与环境产生交互，感知环境的状态，采取动作获得收益，并影响环境的状态，能够从中不断地学习，实现同相轴的最优拾取。

发明内容

本部分的目的在于概述本发明的实施例的一些方面以及简要介绍一些较佳实施例。在本部分以及本申请的说明书摘要和发明名称中可能会做些简化或省略以避免使本部分、说明书摘要和发明名称的目的模糊，而这种简化或省略不能用于限制本发明的范围。

鉴于上述和/或现有的拾取方法中存在的问题，提出了本发明。

因此，本发明所要解决的问题在于如何从复杂多变的成像剖面中自动地、高效地提取期望的反射结构信息。

为解决上述技术问题，本发明提供如下技术方案：

第一方面，本发明实施例提供了一种基于马尔可夫决策过程的反射层位追踪方法，其包括，输入偏移成像剖面，对所述成像剖面进行预处理并设计拾取的参数；

计算所述成像剖面的属性信息，包括包络属性、相位属性、倾角属性、极值点属性和子波质心属性；

确定种子点的个数和位置；

建立马尔可夫决策过程模型，并采用马尔可夫决策过程拾取反射层位信息；

判断拾取结果的有效性，并存储有效的反射层位信息。

作为本发明所述基于马尔可夫决策过程的反射层位追踪方法的一种优选方案，其中：所述拾取参数包括同相轴拾取的最大数目N

作为本发明所述基于马尔可夫决策过程的反射层位追踪方法的一种优选方案，其中：所述确定种子点的个数和位置包括，

根据计算的不同的属性，采用边缘检测的方式，计算不同属性的边缘点；

将不同属性的边缘结果进行聚类，确定对应的中心点位置；

根据同相轴间隔约束和同相轴最多个数约束，依次提取每个同相轴的初始点位置。

作为本发明所述基于马尔可夫决策过程的反射层位追踪方法的一种优选方案，其中：所述马尔可夫决策过程模型包括决策时刻、状态空间、行动空间、状态转移概率和收益函数，其中，在二维的成像剖面中，纵向网格点为nz，横向网格点为nx；状态为成像剖面上纵向的空间位置，状态的集合表示为 S＝[1,...,nz]，动作为从横向位置x

作为本发明所述基于马尔可夫决策过程的反射层位追踪方法的一种优选方案，其中：所述决策时刻为成像剖面上的横向的空间位置，t∈[1,...,nx]，其中，在二维的成像剖面中，纵向深度网格点为nz，横向网格点为nx；

所述状态空间的宽度为

所述行动空间在倾角属性信息的导引下，采用随机性策略实现从状态到动作的映射，公式为：

其中，p为当前t时刻状态s处的倾角，a

所述状态转移概率表示为从t-1时刻状态

所述收益函数是根据相邻道的多属性参数的相似程度确定的，当前决策时刻的收益函数为：

其中，Δf＝||f(s

作为本发明所述基于马尔可夫决策过程的反射层位追踪方法的一种优选方案，其中：所述采用马尔可夫决策过程拾取反射层位信息包括，

当决策时刻到达终点时(剖面的边界)，或者当长期的累计收益小于给定的阈值时(存在断层)，将自动转入到终止状态，停止该条同相轴的追踪，并记录同相轴的长度。

作为本发明所述基于马尔可夫决策过程的反射层位追踪方法的一种优选方案，其中：所述判断拾取结果的有效性，并存储有效的反射层位信息包括，

若追踪的同相轴小于给定的阈值长度，该条同相轴将被舍弃，开始下一个种子点的追踪过程，重新进行追踪；

反之，记录当前追踪的反射结构信息，继续下一个种子点的追踪过程；

当所有的种子点都追踪完毕之后，输出最终的反射结构信息。

第二方面，本发明实施例提供了一种基于马尔可夫决策过程的反射层位追踪系统，其包括：处理模块，用于输入偏移成像剖面，并对所述成像剖面进行预处理；

计算模块，用于设计拾取的参数并计算所述成像剖面的属性信息，根据计算的不同的属性，确定种子点的位置；

建模模块，用于建立马尔可夫决策过程模型，并采用马尔可夫决策过程拾取反射层位信息；

分析模块，用于判断拾取结果的有效性，并存储有效的反射层位信息。

第三方面，本发明实施例提供了一种计算机设备，包括存储器和处理器，所述存储器存储有计算机程序，其中：所述处理器执行所述计算机程序时实现上述方法的任一步骤。

第四方面，本发明实施例提供了一种计算机可读存储介质，其上存储有计算机程序，其中：所述计算机程序被处理器执行时实现上述方法的任一步骤。

本发明有益效果为：

1.成像剖面结果往往无法直接映射到同相轴，通过将地震成像剖面映射到属性空间，可以突显同相轴的差异(纵向的差异和横向的连续性)，从而检测并提取同相轴。然而由于成像剖面包含相干、非相干的噪音，成像剖面变化剧烈。特别是在弱反射情况，出现断层破裂带，不整合面等情况，计算的单一属性结果往往是不正确的，无法追踪到正确的同相轴。多属性可以突出成像剖面性质的不同方面，通过多属性共同约束，可以弥补单一属性的缺点。

2.将成像剖面上的反射同相轴(反射结构)拾取看成是随机环境中的智能决策问题。通过拾取目标与环境交互，感知环境的状态，采取动作获得收益，并影响环境的状态，能够从中不断地学习，实现同相轴的最优拾取。通过考虑最优的马尔可夫最优决策过程，可以避免单道、或局部区域内的波形变化引起的影响，稳健地拾取期望的反射结构信息。

3.自动化的反射层位拾取结果不仅可以应用于速度建模，也可以应用在地震地质层位的解释过程中，具有非常广泛的应用价值。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其它的附图。其中：

图1为基于马尔可夫决策过程的反射结构追踪流程图。

图2为局部偏移成像结果示意图。

图3为局部剖面的包络属性。

图4为局部剖面的相位属性。

图5为局部剖面的倾角属性。

图6为同相轴拾取结果叠加在成像结果之上的拾取结果对比图。

图7为鲸型模型上拾取多个反射结构的结果示意图。

图8为overthrust模型上背景速度下的偏移成像结果。

图9为overthrust模型上拾取的反射层位结果。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合说明书附图对本发明的具体实施方式做详细的说明，显然所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明的保护的范围。

在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是本发明还可以采用其他不同于在此描述的其它方式来实施，本领域技术人员可以在不违背本发明内涵的情况下做类似推广，因此本发明不受下面公开的具体实施例的限制。

其次，此处所称的“一个实施例”或“实施例”是指可包含于本发明至少一个实现方式中的特定特征、结构或特性。在本说明书中不同地方出现的“在一个实施例中”并非均指同一个实施例，也不是单独的或选择性的与其他实施例互相排斥的实施例。

本发明结合示意图进行详细描述，在详述本发明实施例时，为便于说明，表示器件结构的剖面图会不依一般比例作局部放大，而且所述示意图只是示例，其在此不应限制本发明保护的范围。此外，在实际制作中应包含长度、宽度及深度的三维空间尺寸。

同时在本发明的描述中，需要说明的是，术语中的“上、下、内和外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一、第二或第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。

本发明中除非另有明确的规定和限定，术语“安装、相连、连接”应做广义理解，例如：可以是固定连接、可拆卸连接或一体式连接；同样可以是机械连接、电连接或直接连接，也可以通过中间媒介间接相连，也可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。

实施例1

参照图1，为本发明第一个实施例，该实施例提供了一种基于马尔可夫决策过程的反射层位追踪方法，包括：

S100：输入成像剖面，对成像剖面进行一些适当的预处理(主要包括去噪、提升分辨率、振幅均衡等操作)，并设计拾取的参数。

应说明的是，拾取参数主要包括，同相轴拾取的最大数目N

S200：计算成像剖面的属性信息，包括包络属性、相位属性、倾角属性、极值点属性和子波质心属性，共同约束拾取结果。

(1)相位属性

应说明的是，相位属性通过对地震剖面施加希尔伯特变换，生成复信号获得。针对地震剖面，施加希尔伯特变换：

h(z,x)＝h(z,x)+ih'(z,x) (1)

其中，h'(z,x)是h(z,x)的虚部。h(z,x)是复地震道。

其中，env(z,x)为包络，θ(z,x)为瞬时相位。

然而上式求解的瞬时相位的范围为(-180,180)，存在不连续性。因此采用下列计算瞬时相位的余弦值：

进一步地，考虑到env(z,x)为零时，除法会引起不稳定问题，本方法采用最小二乘反演除的方式，优化瞬时相位余弦值的计算：

cos(θ(z,x))＝argmin||h(z,x)-cos(θ(z,x))env(z,x)||

最终，瞬时相位余弦值计算为：

其中，eps为给定的极小阈值，保证计算的稳定性。

(2)倾角属性

应说明的是，倾角属性通过结构张量方法得到。考虑子波的震荡性和噪音的存在，首先对成像剖面进行平滑滤波，然后利用希尔伯特变换，求取包络，利用包络值求取结构张量。

给定成像剖面包络A，其结构张量矩阵可表示为：

其中，

其中λ

(3)极值点属性

应说明的是，极值点的提取主要包括极大值和极小值。可以对波形提取极大值和极小值，也可对其他属性结果(比如包络、相位)提取极值点。求取极值点的方法有很多，比如通过一阶导数，二阶导数等等方法。本方法通过比较在局部窗内的大小来提取极值点。

极大值和极小值分别采用如下策略计算：

其中，f(x)为选择的不同属性值。通过不同的局部窗的多次加权叠加，可以有效避免异常点的存在。另外，根据实际情况，同性轴的拾取要求具有稀疏性，相邻较近的极值点，采用极大值抑制方法进行处理。

(4)子波质心属性

应说明的是，子波(波形)的质心可以表示为该段子波内所有点的一个算术平均：

该方法依赖于局部窗的大小，需要估计出子波的大致范围。因此将该属性应用于拾取过程第二个阶段的状态转移和最优决策过程中，通过第一步的约束，确定局部窗的大小。

S300：建立马尔可夫决策过程的同相轴拾取模型，包括，决策时刻，状态空间，行动空间，状态转移概率和收益函数。

(1)决策时刻

应说明的是，决策时刻定义为成像剖面上的横向的空间位置。t∈[1,...,nx]，(在二维的成像剖面中，纵向深度网格点为nz，横向网格点为nx)。

(2)状态空间

应说明的是，将成像剖面上纵向的空间位置定义为状态，状态的集合表示为S＝[1,...,nz]。考虑到剖面中存在多条同相轴，当拾取某一条同相轴时，纵向上距离较远的位置对当前同相轴的拾取影响较小。本方法将状态空间缩小至指定范围内宽度

(3)行动空间

应说明的是，将动作定义为从横向位置x

其中，p为当前t时刻状态s处的倾角，a

(4)状态转移概率

应说明的是，状态转移概率表示为从t-1时刻状态

(5)收益函数

应说明的是，本方法在高维多属性框架下定义收益函数。根据相邻道的多属性参数的相似程度，确定当前决策时刻的状态收益函数值。因此，当前决策时刻的收益函数可以定义为：

其中，Δf＝||f(s

S400：确定初始种子点的位置，以及状态空间的大小和折扣因子。其中初始种子点对应于状态空间中的初始状态。

种子点的位置确定方式包括：

根据计算的不同的属性，采用边缘检测的方式，计算不同属性的边缘点；

将不同属性的边缘结果进行聚类，确定对应的中心点位置；

根据同相轴间隔约束和同相轴最多个数约束，依次提取每个同相轴的初始点位置。

根据获取的起始时刻的种子点位置

若横向起伏变化不剧烈，亦可采用如下方式定义初始的状态空间大小：

仅计算当前决策时刻状态对后续有限步的影响，将折扣因子定义为：

其中，N

S500：反射同相轴追踪过程并输出反射结构信息。

应说明的是，将同相轴拾取过程定义为高维多属性的马尔可夫状态转移模型之后，就可以根据状态和收益的变化进行追踪求解。

当决策时刻到达终点时(剖面的边界)，或者当长期的累计收益小于给定的阈值时(存在断层)时，自动地转入到终止状态，停止该条同相轴的追踪，并记录同相轴的长度。

应说明的是，若追踪的同相轴小于给定的阈值长度时，该条同相轴将被舍弃，开始下一个种子点的追踪过程，重新进行追踪。反之，记录当前追踪的反射结构信息，重新进行追踪。当所有的种子点都追踪完毕之后，输出最终的反射结构信息。

进一步的，本实施例还提供一种基于马尔可夫决策过程的反射层位追踪系统，包括：

处理模块，用于输入偏移成像剖面，并对所述成像剖面进行预处理；

计算模块，用于设计拾取的参数并计算所述成像剖面的属性信息，根据计算的不同的属性，确定种子点的位置；

建模模块，用于建立马尔可夫决策过程模型，并采用马尔可夫决策过程拾取反射层位信息；

分析模块，用于判断拾取结果的有效性，并存储有效的反射层位信息。

本实施例还提供一种计算机设备，适用于基于马尔可夫决策过程的反射层位追踪方法的情况，包括：

存储器和处理器；存储器用于存储计算机可执行指令，处理器用于执行计算机可执行指令，实现如上述实施例提出的基于马尔可夫决策过程的反射层位追踪方法。

该计算机设备可以是终端，该计算机设备包括通过系统总线连接的处理器、存储器、通信接口、显示屏和输入装置。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的通信接口用于与外部的终端进行有线或无线方式的通信，无线方式可通过WIFI、运营商网络、NFC(近场通信)或其他技术实现。该计算机设备的显示屏可以是液晶显示屏或者电子墨水显示屏，该计算机设备的输入装置可以是显示屏上覆盖的触摸层，也可以是计算机设备外壳上设置的按键、轨迹球或触控板，还可以是外接的键盘、触控板或鼠标等。

本实施例还提供一种存储介质，其上存储有计算机程序，该程序被处理器执行时实现如上述实施例提出的实现基于马尔可夫决策过程的反射层位追踪方法。

本实施例提出的存储介质与上述实施例提出的数据存储方法属于同一发明构思，未在本实施例中详尽描述的技术细节可参见上述实施例，并且本实施例与上述实施例具有相同的有益效果。

实施例2

参照图2～9，为本发明第二个实施例，该实施例提供了一种基于马尔可夫决策过程的反射层位追踪方法，为了验证本发明的有益效果，通过对比实验进行科学论证。

首先测试采用本方法拾取单个同相轴的应用效果。

图2为某一理论模型数据在背景速度不准确情况下偏移成像的局部剖面。

可以看到，由于偏移速度不准，成像结果不聚焦，剖面中出现大量的偏移画弧假象，横向上的子波波形变化较大。

图3～5分别展示了局部剖面的包络、相位和倾角属性。

可以看到，不同的属性能够突出成像剖面中某一部分的特征，如图3中黑框所示，当包络属性横向出现不连续情况时，图4中黑框所示的相位属性却具有良好的连续性，因此采用多属性可以相互补充单一属性的缺陷。

但是也存在所有属性都不好的情况，如图3～5中黑圈所示，此时若基于常规的单道相似性进行追踪，无法追踪到准确的反射结构，无法获得该区域的有效信息。虽然在局部的区域，所有的属性表现都不好，但是在该区域的横向的局部邻域内具有一定的相似性。因此可以利用前后邻域的属性，采用多步的预测方法确定该部分区域的反射结构情况。

图6展示了采用常规的单道互相关拾取和本方法提出的多属性与马尔可夫决策过程结合方法的拾取结果。从中可以看出，横向上的波形剧烈变化和区域上的属性差异导致常规的单道方法无法正确地追求同相轴的位置。而基于高维多属性、多道收益的马尔可夫决策过程方法可以有效地避开局部的噪音干扰，成功地追踪到正确的反射层位。

接下来采用鲸型模型的偏移成像结果测试多个反射结构的拾取。

在本实施例中，种子点的数目最大设置为5，同相轴的间隔最小为0.1km, 同相轴的最小长度为横向点数的1/10。经过多属性计算，进行边缘检测，聚类后获得的种子点如图7中五星所示。在此基础上，经过常规的多属性约束确定状态空间之后，并采用马尔可夫决策方法自动地追踪同相轴，其结果如图7所示。可以看到，采用本方法提取的反射结构与真实的结构相吻合，说明了本方法的有效性。

最后，采用比较复杂的overthrust模型中的2D切面进一步测试本方法的有效性。图8为采用背景速度模型偏移成像结果，图9为采用本方法拾取的反射层位结果。

在本例中，仅提取了8个反射层位。由于背景速度不准，偏移成像剖面中的反射层位出现不聚焦、不连续现象，使得采用常规的方法难以追踪到正确的反射层位信息。而采用本方法提出的马尔可夫决策过程可以追踪到正确的反射层位信息，从而可以有选择性地筛选出期望的反射层位，应用到速度建模和地震地震解释过程中。

应说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。