CN116958470B - 一种融合马尔可夫链和多点统计学的地质建模方法和装置 - Google Patents

Abstract

本发明公开了一种融合马尔可夫链和多点统计学的地质建模方法和装置，所述方法包括：利用指示函数提取出三维网格图像中的转移面中的若干个非平稳地质特征和平稳地质特征，并基于马尔可夫链模型和反距离加权法得到在所述待模拟区域中各模拟点的概率分布，构建出所有非平稳特征对应的非平稳地质结构三维模型；将非平稳地质结构三维模型作为硬约束数据，并基于多尺度策略并根据最大期望算法对初始三维地质模型进行迭代优化，直至累计迭代次数达到预设迭代次数时，得到最终的三维地质模型。本发明通过融合马尔可夫链和反距离加权思想来明确地质剖面结构信息之间的耦合机制，并且能够更有效地重建出非平稳地质特征的空间结构。

Description

一种融合马尔可夫链和多点统计学的地质建模方法和装置

技术领域

本发明涉及地质建模技术领域，尤其涉及一种融合马尔可夫链和多点统计学的地质建模方法和装置。

背景技术

在多点统计学方法中，地质结构的统计特征是随机函数根据训练图像中具有相同的样本值的概率分布来推断的。受控于构造、沉积机制和沉积环境等多重因素影响，地质体的实际空间结构具有典型的、绝对的非平稳性特征。因此，处理非平稳性的能力对于多点统计学重建复杂地质结构至关重要。

为此，国内外许多学者研究开发出不同的方法来模拟非平稳地质现象，现有的方法主要分为三种：一是几何变化法，该几何变化法利用平稳的训练图像，通过划分区域进行相应的几何变化，对数据事件进行旋转和缩放变换，建立对应的模式库，以重建非平稳地质特征，然而，此类方法难以模拟出局部区域地质结构存在显著差异的非平稳特征时对应的地质结构重建，尤其是具有方向延展性特征的非平稳地质结构，如断裂或透镜体等，因此现有技术难以明确地质剖面结构信息之间的耦合机制，使得基于地质剖面结构的三维地质模型不能有效重建出非平稳地质特征的空间结构。

发明内容

本发明实施例提供了一种融合马尔可夫链和多点统计学的地质建模方法和装置，能有效解决现有技术中难以明确地质剖面结构信息之间的耦合机制，使得基于地质剖面结构的三维地质模型不能有效重建出非平稳地质特征的空间结构的问题。

本发明一实施例提供一种融合马尔可夫链和多点统计学的地质建模方法，包括：

获取目标二维地质剖面；

将目标二维地质剖面导入到三维模拟网格中，得到所述目标二维地质剖面在三维空间中表示的二维网格图像，并通过三维模拟网络将所述二维网格图像转换为三维网格图像；

利用指示函数对三维网格图像中的非平稳地质结构和平稳地质结构进行分别标识；

对于所述三维网格图像中的每一与三维网格图像平行的转移面，提取出转移面中的若干个非平稳地质特征和若干个平稳地质特征，并设置所有非平稳特征对应的待模拟区域的尺寸，基于马尔可夫链模型和反距离加权法得到在所述待模拟区域中各模拟点的概率分布，采用随机抽样方法对各模拟点进行抽样得到若干选定模拟点，基于若干选定模拟点在所述待模拟区域中构建出所有非平稳特征对应的非平稳地质结构三维模型；

基于多点统计学算法，从三维网格图像中提取各地质结构的空间分布模式，并基于各空间分布模式建立地质对象空间模式数据库，根据所述空间模式数据库对所述初始地质模型进行序贯模拟，得到初始三维地质模型；

将非平稳地质结构三维模型作为硬约束数据，基于多尺度策略并根据最大期望算法对初始三维地质模型进行迭代优化，直至累计迭代次数达到预设迭代次数时，得到最终的三维地质模型。

优选地，所述通过三维模拟网格将所述二维网格图像转换为三维网格图像，具体包括：

对于二维网格图像中每一待扩展网格单元，获取待扩展网格单元的各待扩展点；

对于每一待扩展点，在与待扩展点处于相同水平位置上的若干相邻点的属性值中随机选出待扩展点的属性值；

根据所有待扩展点的属性值对各待扩展点进行扩展，得到每一待扩展网格单元对应的已扩展网格单元；

根据每一待扩展网格单元对应的已扩展网格单元，得到三维网格图像。

优选地，所述利用指示函数对三维网格图像中的非平稳地质结构和平稳地质结构进行分别标识，具体包括：

利用指示函数将非平稳地质结构的属性值标识为1，将平稳地质结构的属性值都标识为0。

优选地，所述提取出转移面中的若干个非平稳地质特征和若干个平稳地质特征，并设置所有非平稳特征对应的待模拟区域的尺寸，具体包括：

提取出转移面中的若干个非平稳地质特征和若干个平稳地质特征；

在判断转移面中所有提取的特征均为非平稳地质特征时，则计算所有非平稳地质特征的空间位置坐标，并选出最大坐标值和最小坐标值，根据最大坐标值和最小坐标值来计算出外接长方体的尺寸，并将外接长方体的尺寸作为所有非平稳特征对应的待模拟区域的尺寸；

在判断转移面中所有提取的特征至少有一个平稳地质特征时，则根据预设数值来设定所有非平稳特征对应的待模拟区域的尺寸。

优选地，所述基于马尔可夫链模型和反距离加权法得到在所述待模拟区域中各模拟点的概率分布，具体包括：

基于马尔可夫链模型，根据转移面中非平稳地质特征和平稳地质特征的个数之和得到转移计数矩阵，并根据转移计数矩阵得到各非平稳地质特征的转移概率矩阵；

根据各非平稳地质特征的转移概率矩阵，得到在转移面中沿着转移方向的起始非平稳特征和终止非平稳特征的转移概率；

根据如下公式计算出转移面的转移概率：

其中，M_TP(i)表示在预设步长i下转移面的转移概率，dis表示起始非平稳特征和终止非平稳特征之间的距离，TP₁和TP₂分别表示起始非平稳特征和终止非平稳特征的转移概率；

基于反距离加权法，根据各非平稳地质特征的形态拟合出转移面的中线，并遍历所述转移面内的所有点，计算每个点到非平稳特征中线的距离，选出最远的距离；

根据选出的最远的距离，通过以下公式计算待模拟区域中各模拟点到转移面的中线的距离权重：

其中，Z为模拟点到中线的距离权重，d(x,y)为模拟点到中线的距离，d_max为每个点到非平稳特征中线的距离中最远的距离；

根据转移面的转移概率和转移面的中线的距离权重计算得到待模拟区域中各模拟点的概率分布。

优选地，所述根据转移面的转移概率和转移面的中线的距离权重计算得到待模拟区域中各模拟点的概率分布，具体包括：

根据转移面的转移概率和转移面的中线的距离权重，并通过以下公式计算得到待模拟区域中各模拟点的概率分布：

F_Tran(x,y)＝M_TP(i)×I_(k)×Z

其中，F_Tran(x,y)表示模拟点(x,y)的概率分布，M_TP(i)表示在预设步长i下转移面的转移概率，Z表示模拟点到转移面的中线的距离权重。

优选地，所述根据所述空间模式数据库对所述初始地质模型进行序贯模拟，得到初始三维地质模型，具体包括：

按照预设的顺序初始化模拟网格，从未模拟边缘区域空网格开始，从外向内逐层建模，并设定重叠区域的大小，以使得移动到每个位置的模型都有与已建模的区域重合的部分，对于重叠区域中的每个待模拟点，从模拟网格中提取出当前待模拟点的位置的模式，然后搜索空间模式数据库中相应的模式来进行序贯模拟，直到得到初始三维地质模型。

在上述的方法实施例的基础上，本发明对应提供了装置项实施例。

本发明一实施例提供了一种融合马尔可夫链和多点统计学的地质建模装置，包括：数据获取模块、三维网格图像生成模块、标识模块、非平稳地质结构三维模型构建模块、三维地质模型第一生成模块和三维地质模型第二生成模块；

所述数据获取模块，用于获取目标二维地质剖面；

所述三维网格图像生成模块，用于将目标二维地质剖面导入到三维模拟网格中，得到所述目标二维地质剖面在三维空间中表示的二维网格图像，并通过三维模拟网络将所述二维网格图像转换为三维网格图像；

所述标识模块，用于利用指示函数对三维网格图像中的非平稳地质结构和平稳地质结构进行分别标识；

所述非平稳地质结构三维模型构建模块，用于对于所述三维网格图像中的每一与三维网格图像平行的转移面，提取出转移面中的若干个非平稳地质特征和若干个平稳地质特征，并设置所有非平稳特征对应的待模拟区域的尺寸，基于马尔可夫链模型和反距离加权法得到在所述待模拟区域中各模拟点的概率分布，采用随机抽样方法对各模拟点进行抽样得到若干选定模拟点，基于若干选定模拟点在所述待模拟区域中构建出所有非平稳特征对应的非平稳地质结构三维模型；

所述三维地质模型第一生成模块，用于基于多点统计学算法，从三维网格图像中提取各地质结构的空间分布模式，并基于各空间分布模式建立地质对象空间模式数据库，根据所述空间模式数据库对所述初始地质模型进行序贯模拟，得到初始三维地质模型；

所述三维地质模型第二生成模块，用于将非平稳地质结构三维模型作为硬约束数据，基于多尺度策略并根据最大期望算法对初始三维地质模型进行迭代优化，直至累计迭代次数达到预设迭代次数时，得到最终的三维地质模型。

通过实施本发明具有如下有益效果：

本发明实施例提供了一种融合马尔可夫链和多点统计学的地质建模方法和装置，本发明的方法首先通过指示函数对三维网格图像中的非平稳地质结构和平稳地质结构进行分别标识；然后对于所述三维网格图像中的每一与三维网格图像平行的转移面，提取出转移面中的若干个非平稳地质特征和平稳地质特征，并基于马尔可夫链模型和反距离加权法得到在所述待模拟区域中各模拟点的概率分布，采用随机抽样方法对各模拟点进行抽样得到若干选定模拟点，基于若干选定模拟点在所述待模拟区域中构建出所有非平稳特征对应的非平稳地质结构三维模型；因此本发明可以借用马尔可夫链的基本思想保证非平稳地质体的方向延展性，而又引入反距离加权的概念来考虑非平稳地质体的空间位置信息，从而可以建立出有效表达空间关联关系的非平稳地质结构三维模型，将得到的非平稳地质结构三维模型作为硬约束数据，并基于多尺度策略并根据最大期望算法对初始三维地质模型进行迭代优化，直至累计迭代次数达到预设迭代次数时，得到最终的三维地质模型，而最终的三维地质模型因为基于非平稳地质结构三维模型作为硬约束数据得到的，因此最终的三维地质模型中的非平稳地质特征的方向延展性会更好，且非平稳特征的模拟结果更平滑和准确。与现有技术相比，本发明可以通过融合马尔可夫链和反距离加权思想来明确地质剖面结构信息之间的耦合机制，实现了地质结构全局特征的重建，使得基于地质剖面结构的最终的三维地质建模模型能够更有效地重建出非平稳地质特征的空间结构。

附图说明

图1是本发明一实施例提供的一种融合马尔可夫链和多点统计学的地质建模方法的流程示意图。

图2是本发明一实施例提供的将二维网格图像转换为三维网格图像的示意图。

图3是本发明一实施例提供的对模拟点进行反距离加权的示意图。

图4是本发明一实施例提供的模拟点到转移面的中线示意图。

图5是本发明一实施例提供的序贯模拟过程中重叠区的类型示意图。

图6是本发明一实施例提供的一种融合马尔可夫链和多点统计学的地质建模装置的结构示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，是本发明一实施例提供的一种融合马尔可夫链和多点统计学的地质建模方法的流程示意图；

本发明一实施例提供了一种融合马尔可夫链和多点统计学的地质建模方法，包括：

步骤S1：获取目标二维地质剖面；

步骤S2：将目标二维地质剖面导入到三维模拟网格中，得到所述目标二维地质剖面在三维空间中表示的二维网格图像，并通过三维模拟网络将所述二维网格图像转换为三维网格图像；

步骤S3：利用指示函数对三维网格图像中的非平稳地质结构和平稳地质结构进行分别标识；

步骤S4：对于所述三维网格图像中的每一与三维网格图像平行的转移面，提取出转移面中的若干个非平稳地质特征和若干个平稳地质特征，并设置所有非平稳特征对应的待模拟区域的尺寸，基于马尔可夫链模型和反距离加权法得到在所述待模拟区域中各模拟点的概率分布，采用随机抽样方法对各模拟点进行抽样得到若干选定模拟点，基于若干选定模拟点在所述待模拟区域中构建出所有非平稳特征对应的非平稳地质结构三维模型；

步骤S5：基于多点统计学算法，从三维网格图像中提取各地质结构的空间分布模式，并基于各空间分布模式建立地质对象空间模式数据库，根据所述空间模式数据库对所述初始地质模型进行序贯模拟，得到初始三维地质模型；

步骤S6：将非平稳地质结构三维模型作为硬约束数据，基于多尺度策略并根据最大期望算法对初始三维地质模型进行迭代优化，直至累计迭代次数达到预设迭代次数时，得到最终的三维地质模型。

对于步骤S1，在一个优选的实施例中，获取的目标二维地质剖面中包括了复杂地质体几何形态，比如具有断裂或透镜体这些非平稳地质特征。

对于步骤S2，在一个优选的实施例中，将目标二维地质剖面导入到三维模拟网格中，得到所述目标二维地质剖面在三维空间中表示的二维网格图像，并通过三维模拟网络将所述二维网格图像转换为三维网格图像，具体包括：

将二维地质剖面导入到三维模拟网格中，得到三维空间中表示的二维网格图像，并在三维模拟网格中将二维地质剖面转换为三维训练图像，从而得到三维网格图像。

其中，通过三维模拟网格将所述二维网格图像转换为三维网格图像，具体包括：

对于二维网格图像中每一待扩展网格单元，获取待扩展网格单元的各待扩展点；

对于每一待扩展点，在与待扩展点处于相同水平位置上的若干相邻点的属性值中随机选出待扩展点的属性值；

根据所有待扩展点的属性值对各待扩展点进行扩展，得到每一待扩展网格单元对应的已扩展网格单元；

根据每一待扩展网格单元对应的已扩展网格单元，得到三维网格图像。

具体的，在将所述二维网格图像转换为三维网格图像时，如图2所示，(a)左图为二维网格图像的俯视图，右图为扩展后的三维网格图像俯视图；(b)为不同位置的待扩展网格单元的扩展过程示意图，空白区域为扩展区域，P和Q为待扩展网格单元，且P和Q的地质属性值不同，即P和Q的地质结构不同；

本发明在将待扩展网格单元从二维逐层扩展到三维的过程中，由于二维待扩展网格单元内大多数地层是水平分布的，为了避免地层层序错乱，所以待扩展网格单元的属性只从水平方向上相邻的属性值中选择。

具体的，若当前扩展点是顶点，如图2中(b)的当前扩展点P点，则从待扩展网格单元中同一水平位置上与P相邻的两个点的属性值中随机选取点P的属性；

若当前扩展点不是顶点，如图2中(b)的Q点，Q点的属性值是从待扩展网格单元中相同水平位置上三个最邻近点的属性值中随机选择的。

如果待扩展网格单元处于中间(如图2(b)左图)需要向两侧扩展；如果待扩展网格单元处于最外层(如图2(b)右图)需要向内进行扩展，填充完靠近待扩展网格单元的那一层后填充下一层，直到达到需要填充的层数。从而使得最后获得的三维网格图像，其厚度大于或等于一个预设模板的大小，如图2(a)右图所示。

将待扩展网格单元从二维逐层扩展到三维，可以获得三维的空间模式，以进行后续的序贯模拟和迭代过程。

对于步骤S3，在一个优选的实施例中，利用指示函数对三维网格图像中的非平稳地质结构和平稳地质结构进行分别标识，具体包括：

利用指示函数将非平稳地质结构的属性值标识为1，将平稳地质结构的属性值都标识为0。

因为非平稳地质结构具有局部各向异性，为了获得更好的模拟结果，就不能与其他平稳地层一起模拟，所以需要在扩展完成后的三维网格图像中提取出非平稳特征，单独对其进行三维重构。而本发明通过引入指示函数I(k)来标识每个与非平稳特征有关的单元体(即非平稳地质结构)，本发明通过指示函数可以把非平稳地质结构的属性值标识为1，其余地层或地质结构(即平稳地质结构)的属性值都标识为0，从而可以使得后续的步骤S4中可以实现识别并提取出非平稳地质结构特征。

对于步骤S4，在一个优选的实施例中，对于所述三维网格图像中的每一与三维网格图像平行的转移面，提取出转移面中的若干个非平稳地质特征和若干个平稳地质特征，并设置所有非平稳特征对应的待模拟区域的尺寸，具体包括：

提取出转移面中的若干个非平稳地质特征和若干个平稳地质特征；需要说明的是，转移面一般情况下是平行于三维网格图像的一个纵剖面；

在判断转移面中所有提取的特征均为非平稳地质特征时，则计算所有非平稳地质特征的空间位置坐标，并选出最大坐标值和最小坐标值，根据最大坐标值和最小坐标值来计算出外接长方体的尺寸，并将外接长方体的尺寸作为所有非平稳特征对应的待模拟区域的尺寸；

在判断转移面中所有提取的特征至少有一个平稳地质特征时，则根据预设数值来设定所有非平稳特征对应的待模拟区域的尺寸。

因此，为了提高非平稳结构的模拟效率和准确性，如果模拟使用的转移面上都存在非平稳特征，则计算所有转移面上非平稳特征的空间位置坐标，选出最大和最小坐标值，据此计算出非平稳特征分布范围的外接长方体，在后续对非平稳特征进行模拟时即在这个长方体中进行；

如果模拟使用的转移面中至少有一个不存在非平稳特征，导致转移面间无法建立非平稳特征的空间对应关系，可以根据先验经验和所得的数据资料设定非平稳特征的延展长度及其所影响的模拟网格，从而得到外接长方体，即得到待模拟区域的尺寸大小。

进一步的，在获得非平稳特征对应的待模拟区域的尺寸后，基于马尔可夫链模型和反距离加权法得到在所述待模拟区域中各模拟点的概率分布，然后采用随机抽样方法对各模拟点进行抽样得到若干选定模拟点，基于若干选定模拟点在所述待模拟区域中构建出所有非平稳特征对应的非平稳地质结构三维模型；其中，基于马尔可夫链模型和反距离加权法得到在所述待模拟区域中各模拟点的概率分布，具体包括：

马尔可夫链是指具有马尔可夫性质的随机变量的一个数列，描述了一种状态序列，其每个状态值取决于前面有限个状态；马尔可夫链因其具有普遍意义额随机过程，被广泛地应用于近代物理、生物医学、地质学、农学等。

马尔可夫过程采用数学方式可描述为：在时间序列中，假设某一随机过程{x_n}，n＝1,2,…i构成该过程的状态集，在时刻t处于状态x_i，那么在下个时刻t+1处于x_j的概率是(schwarzacherW,2000)：

P_ij(t+1)＝P(X_t+1＝x_j|X₀＝x₀,X₁＝x₁,…,X_t-1＝x_i-1,X_t＝x_i)

＝P(X_t+1＝x_j|X_t＝x_i)

在马尔可夫过程中，与不同状态之间的转移或相同状态之间的转移相关的概率叫做转移概率。使用马尔可夫链模型进行预测模拟，最关键就是要确定或近似计算出马尔可夫状态的转移概率矩阵，数学表达式为：

在随机模拟过程中，时间参数也可以换作其他的量例如空间位置。例如，四个状态{A,B,C,D}为从左向右转移的一维马尔可夫链，转移计数矩阵是从四个状态分布的短序列中分别统计出当前状态转移到不同状态的次数，计算出转移计数矩阵后，根据转移概率矩阵各行元素之和为1的性质，可以得到该序列的四个状态的转移概率矩阵。

而本发明技术借鉴马尔可夫链在状态转移过程具有方向的特性，来预测具有方向延展性的非平稳结构的空间展布，在使用马尔可夫链模型进行非平稳特征模拟过程中，根据转移面中非平稳地质特征和平稳地质特征的个数之和得到转移计数矩阵，并根据转移计数矩阵得到各非平稳地质特征的转移概率矩阵；

根据各非平稳地质特征的转移概率矩阵，得到在转移面中沿着转移方向的起始非平稳特征和终止非平稳特征的转移概率；

根据如下公式计算出转移面的转移概率：

其中，M_TP(i)表示在预设步长i下转移面的转移概率，dis表示起始非平稳特征和终止非平稳特征之间的距离，TP₁和TP₂分别表示起始非平稳特征和终止非平稳特征的转移概率；

示意性的，在使用马尔可夫链模型进行非平稳特征模拟过程中，首先根据指示函数标识出来的非平稳地质特征的属性值和平稳地质特征(如围岩等)属性值，统计转移面内非平稳地质特征的单元体和平稳地质特征的单元体的个数，然后建立转移计数矩阵，根据转移概率各行和为1的性质，得到转移面的初始转移概率矩阵。

接着计算转移面之间的距离，也就是转移距离，可以采用等间隔的转移步长，所以需要对转移面中的非平稳地质特征的距离进行等分，每份距离间隔即为一步的转移步长，示意性的，预设的转移步长为1；

最后可以沿着转移方向根据起始非平稳TI到终止非平稳TI上非平稳属性和围岩属性的变化趋势，计算出每步转移面内，即计算得出转移面的转移概率。

在反距离加权中，每个模拟点都被认为会受到其邻近点的影响，离模拟点越近，那么对模拟点的影响越大；距离越远，影响越小。因此，可以通过对模拟点附近的样本点的值加权平均，估计出模拟点的值。

基于反距离加权法，根据各非平稳地质特征的形态拟合出转移面的中线，并遍历所述转移面内的所有点，计算每个点到非平稳特征中线的距离，选出最远的距离；

根据选出的最远的距离，通过以下公式计算待模拟区域中各模拟点到转移面的中线的距离权重：

其中，Z为模拟点到中线的距离权重，d(x,y)为模拟点到中线的距离，d_max为每个点到非平稳特征中线的距离中最远的距离；

根据各非平稳地质特征的形态拟合出转移面的中线，具体包括：

如图3所示，P₀为模拟点，其周围分布有n个已知点，测量模拟点P₀到附近所有点的距离d_m，利用以下公式计算出附近点所占的权重：

其中，d_m距离通常使用以下的欧式距离计算公式来计算：

拟合过程如图4所示，图4中的1为非平稳特征，图4中的2为非平稳特征的最小外接矩形，3为拟合的中线；

在计算出一非平稳特征的最小外接矩形后，可以通过规则的矩形计算得到拟合中线。按同样的方法沿着转移方向根据起始非平稳特征和终止非平稳特征的变化趋势，拟合出转移面的中线。具体的，在转移面中，沿着转移方向根据起始非平稳特征的中线L1到终止位置上的终止非平稳特征的中线L2的变化趋势，计算得出每步转移面内非平稳特征中线L；在一个优选的实施例中，终止位置可以是另一个转移剖面的终止位置；

然后遍历该转移面内所有点，计算每个点到转移面的中线的距离，选出最远的距离d_max；

再依次计算每个模拟点到中线的距离d(x,y)占最远距离d_max比重的倒数，得到模拟点到转移面的中线的距离权重，计算公式如下所示：

其中，Z为模拟点到中线的距离权重，d(x,y)为模拟点到中线的距离，d_max为每个点到非平稳特征中线的距离中最远的距离；其中，点到线的距离计算公式为：

马尔可夫链有利于实现非平稳特征的延展性，反距离加权法包括非平稳特征的空间位置信息，在待模拟区域，每个空间点上的属性概率分布不仅与所在周边点有关系，而且还与两侧约束数据的位置有关。为此，本发明技术提出了融合马尔可夫链和反距离加权思想的非平稳空间结构的转换函数，即根据转移面的转移概率和转移面的中线的距离权重，并通过以下转换函数的公式计算得到待模拟区域中各模拟点的概率分布：

F_Tran(x,y)＝M_TP(i)×I_(k)×Z

其中，F_Tran(x,y)表示模拟点(x,y)的概率分布，M_TP(i)表示在预设步长i下转移面的转移概率，Z表示模拟点到转移面的中线的距离权重。

本发明在进行非平稳特征模拟过程中，采用马尔可夫链来约束非平稳特征延展的方向，但是仅使用马尔可夫链模拟，模拟结果中非平稳特征分布碎片、离散。反距离加权思想可以把点的空间位置考虑进去，且算法简单易于实现，从而以使非平稳特征的模拟结果更平滑、准确。即本发明借用马尔可夫链的基本思想保证非平稳地质体的方向延展性，引入反距离加权的概念来考虑非平稳地质体的空间位置信息，实现了非平稳地质特征的三维重构，从而可以建立出有效表达空间关联关系的非平稳地质结构三维模型。

对于步骤S5，在一个优选的实施例中，基于多点统计学算法，从三维网格图像中提取各地质结构的空间分布模式，并基于各空间分布模式建立地质对象空间模式数据库，根据所述空间模式数据库对所述初始地质模型进行序贯模拟，得到初始三维地质模型；

在多点统计学方法中，地质结构的统计特征是随机函数根据训练图像中具有相同的样本值的概率分布来推断的，受控于构造、沉积机制和沉积环境等多重因素影响，地质体的实际空间结构具有典型的以及绝对的非平稳性特征。因此，处理非平稳性的能力对于多点统计学重建复杂地质结构至关重要；

以二维地质剖面为建模数据源，本发明建立融合转移概率(即通过基于马尔可夫链模型和反距离加权法得到的)和多点统计学的非平稳地质结构三维重建技术方法，首先是从三维网格图像中提取各地质结构的空间分布模式，并基于各空间分布模式建立地质对象空间模式数据库，根据所述空间模式数据库对所述初始地质模型进行序贯模拟，得到初始三维地质模型，具体包括：

按照预设的顺序初始化模拟网格，从未模拟边缘区域空网格开始，从外向内逐层建模，并设定重叠区域的大小，以使得移动到每个位置的模型都有与已建模的区域重合的部分，对于重叠区域中的每个待模拟点，从模拟网格中提取出当前待模拟点的位置的模式，然后搜索空间模式数据库中相应的模式来进行序贯模拟，直到得到初始三维地质模型。

示意性的，序贯模拟是按照指定路径从已模拟区域向未知区域推进的过程，本发明采用序贯模拟进行初始化，沿着从外到内逐层模拟的路径填充模拟网格中属性空白区域，这种路径有利于地质结构特征受周边转移面的约束，从而确保地层层序、空间分布和形态特征的一致性及地质意义。

为了提高序贯模拟过程中的模型检索性能，在建模之前先把二维网格图像拓展为三维网格图像，然后从生成的三维网格图像中提取所有的三维模式库，并对所有的三维模板库进行分类，建立一个三维模式库。分类时，需要根据建模路径考虑所有可能的重叠区域的形状类型(如图5所示)，然后用汉明距离计算候选模式与已经模拟网格的重叠区域之间的相似度，根据相似度值来对所有模式进行分类。

在本发明中，通过序贯模拟按照一定顺序来初始化模拟网格，从未模拟边缘区域空网格开始，从外向内逐层建模，自定义重叠区域大小，确保移动到每个位置的模型都有与已经建模区域重合的部分。对于每个待模拟点，从模拟网格中提取出当前位置的模式，然后搜索模式库中相应的模式。

当搜索与当前模拟点处数据事件时，根据重叠区域，将该数据事件与模式库中每个类别的代表性模式进行比较，如果相似度达到阈值，则从模式库中随机选择一个该类别的模式，粘贴到当前数据事件中；如果在遍历完所有模式类别之后，没有找到与阈值相似的模式，则降低阈值，并与每个类别的代表性模式中再次进行比较。

由于依据重叠区相似度在模式库中搜索模式时，多点统计学算法有时难以找到与重叠区相似度满足阈值、且较为吻合的模式。这说明三维模式库中所包含的模式信息不够，我们需要降低阈值搜索模式。如此一来，拼贴时选择的模式和已模拟区域之间的连续性受到破坏。所以，序贯模拟得到的结果中存在大量明显拼贴痕迹，不能作为有效的最终模型。但通过多尺度EM迭代可以有效修正模型中的不连续特征来优化模型，使得模拟网格中的三维地质特征更贴近实际。

对于步骤S6，在一个优选的实施例中，将非平稳地质结构三维模型作为硬约束数据，基于多尺度策略并根据最大期望算法(EM算法)对初始三维地质模型进行迭代优化，直至累计迭代次数达到预设迭代次数时，得到最终的三维地质模型，具体包括：

本发明利用GOSIM算法所使用的多尺度策略进行迭代，首先在得到使用序贯模拟获得的初始模型后，对于每个转移面，从最粗尺度开始模拟，并使得尺度逐步在每个维度增加一倍，当模拟进行到最精细尺度时，把模拟好的非平稳地质结构三维模型载入最精细模拟网格中，作为硬数据约束，根据EM算法对初始三维地质模型进行优化，且在每一个尺度模拟中，都会执行若干次EM算法进行迭代优化，且当前尺度模拟完成后，其结果会作为下一尺度进行模拟的初始模型，直至累计迭代次数达到预设迭代次数时，得到最终的三维地质模型。

本发明通过将模拟的非平稳结构作为迭代过程中的硬约束数据，可以使模拟结果不会出现地质语义错乱，如地层混乱的问题，并且使得最终的三维地质模型因为基于非平稳地质结构三维模型作为硬约束数据得到的，而非平稳地质结构三维模型又是通过借用马尔可夫链的基本思想保证非平稳地质体的方向延展性，引入反距离加权的概念来考虑非平稳地质体的空间位置信息进行建立的，因此最终的三维地质模型中的非平稳地质特征的方向延展性会更好，且非平稳特征的模拟结果更平滑和准确。

与现有技术相比，本发明可以通过融合马尔可夫链和反距离加权思想来明确地质剖面结构信息之间的耦合机制，实现了地质结构全局特征的重建，使得基于地质剖面结构得到的最终的三维地质建模模型能够更有效地重建出非平稳地质特征的空间结构。

如图6所示，在上述各种融合马尔可夫链和多点统计学的地质建模方法的实施例的基础上，本发明对应提供了装置项实施例；

本发明一实施例提供了一种融合马尔可夫链和多点统计学的地质建模装置，包括：数据获取模块、三维网格图像生成模块、标识模块、非平稳地质结构三维模型构建模块、三维地质模型第一生成模块和三维地质模型第二生成模块；

所述数据获取模块，用于获取目标二维地质剖面；

所述三维网格图像生成模块，用于将目标二维地质剖面导入到三维模拟网格中，得到所述目标二维地质剖面在三维空间中表示的二维网格图像，并通过三维模拟网络将所述二维网格图像转换为三维网格图像；

所述标识模块，用于利用指示函数对三维网格图像中的非平稳地质结构和平稳地质结构进行分别标识；

所述非平稳地质结构三维模型构建模块，用于对于所述三维网格图像中的每一与三维网格图像平行的转移面，提取出转移面中的若干个非平稳地质特征和若干个平稳地质特征，并设置所有非平稳特征对应的待模拟区域的尺寸，基于马尔可夫链模型和反距离加权法得到在所述待模拟区域中各模拟点的概率分布，采用随机抽样方法对各模拟点进行抽样得到若干选定模拟点，基于若干选定模拟点在所述待模拟区域中构建出所有非平稳特征对应的非平稳地质结构三维模型；

所述三维地质模型第一生成模块，用于基于多点统计学算法，从三维网格图像中提取各地质结构的空间分布模式，并基于各空间分布模式建立地质对象空间模式数据库，根据所述空间模式数据库对所述初始地质模型进行序贯模拟，得到初始三维地质模型；

所述三维地质模型第二生成模块，用于将非平稳地质结构三维模型作为硬约束数据，基于多尺度策略并根据最大期望算法对初始三维地质模型进行迭代优化，直至累计迭代次数达到预设迭代次数时，得到最终的三维地质模型。

需说明的是，以上所描述的装置实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，既可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。另外，本发明提供的装置实施例附图中，模块之间的连接关系表示它们之间具有通信连接，具体可以实现为一条或多条通信总线或信号线。本领域普通技术人员在不付出创造性劳动的情况下，即可以理解并实施。

所属领域的技术人员可以清楚地了解到，为的方便和简洁，上述描述的装置的具体工作过程，可参考前述方法实施例中对应的过程，在此不再赘述。

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。

Claims (7)

1.一种融合马尔可夫链和多点统计学的地质建模方法，其特征在于，包括：

获取目标二维地质剖面；

将目标二维地质剖面导入到三维模拟网格中，得到所述目标二维地质剖面在三维空间中表示的二维网格图像，并通过三维模拟网络将所述二维网格图像转换为三维网格图像；

利用指示函数对三维网格图像中的非平稳地质结构和平稳地质结构进行分别标识；

对于所述三维网格图像中的每一与三维网格图像平行的转移面，提取出转移面中的若干个非平稳地质特征和若干个平稳地质特征，并设置所有非平稳特征对应的待模拟区域的尺寸，基于马尔可夫链模型和反距离加权法得到在所述待模拟区域中各模拟点的概率分布，采用随机抽样方法对各模拟点进行抽样得到若干选定模拟点，基于若干选定模拟点在所述待模拟区域中构建出所有非平稳特征对应的非平稳地质结构三维模型；

基于多点统计学算法，从三维网格图像中提取各地质结构的空间分布模式，并基于各空间分布模式建立地质对象空间模式数据库，根据所述空间模式数据库对初始地质模型进行序贯模拟，得到初始三维地质模型；

将非平稳地质结构三维模型作为硬约束数据，基于多尺度策略并根据最大期望算法对初始三维地质模型进行迭代优化，直至累计迭代次数达到预设迭代次数时，得到最终的三维地质模型；

其中，所述基于马尔可夫链模型和反距离加权法得到在所述待模拟区域中各模拟点的概率分布，具体包括：

基于马尔可夫链模型，根据转移面中非平稳地质特征和平稳地质特征的个数之和得到转移计数矩阵，并根据转移计数矩阵得到各非平稳地质特征的转移概率矩阵；

根据各非平稳地质特征的转移概率矩阵，得到在转移面中沿着转移方向的起始非平稳特征和终止非平稳特征的转移概率；

根据如下公式计算出转移面的转移概率：

其中，M_TP(i)表示在预设步长i下转移面的转移概率，dis表示起始非平稳特征和终止非平稳特征之间的距离，TP₁和TP₂分别表示起始非平稳特征和终止非平稳特征的转移概率；

基于反距离加权法，根据各非平稳地质特征的形态拟合出转移面的中线，并遍历所述转移面内的所有点，计算每个点到非平稳特征中线的距离，选出最远的距离；

根据选出的最远的距离，通过以下公式计算待模拟区域中各模拟点到转移面的中线的距离权重：

其中，Z为模拟点到中线的距离权重，d(x,y)为模拟点到中线的距离，d_max为每个点到非平稳特征中线的距离中最远的距离；

根据转移面的转移概率和转移面的中线的距离权重计算得到待模拟区域中各模拟点的概率分布。

2.如权利要求1所述的一种融合马尔可夫链和多点统计学的地质建模方法，其特征在于，所述通过三维模拟网格将所述二维网格图像转换为三维网格图像，具体包括：

对于二维网格图像中每一待扩展网格单元，获取待扩展网格单元的各待扩展点；

对于每一待扩展点，在与待扩展点处于相同水平位置上的若干相邻点的属性值中随机选出待扩展点的属性值；

根据所有待扩展点的属性值对各待扩展点进行扩展，得到每一待扩展网格单元对应的已扩展网格单元；

根据每一待扩展网格单元对应的已扩展网格单元，得到三维网格图像。

3.如权利要求1所述的一种融合马尔可夫链和多点统计学的地质建模方法，其特征在于，所述利用指示函数对三维网格图像中的非平稳地质结构和平稳地质结构进行分别标识，具体包括：

利用指示函数将非平稳地质结构的属性值标识为1，将平稳地质结构的属性值都标识为0。

4.如权利要求1所述的一种融合马尔可夫链和多点统计学的地质建模方法，其特征在于，所述提取出转移面中的若干个非平稳地质特征和若干个平稳地质特征，并设置所有非平稳特征对应的待模拟区域的尺寸，具体包括：

提取出转移面中的若干个非平稳地质特征和若干个平稳地质特征；

在判断转移面中所有提取的特征均为非平稳地质特征时，则计算所有非平稳地质特征的空间位置坐标，并选出最大坐标值和最小坐标值，根据最大坐标值和最小坐标值来计算出外接长方体的尺寸，并将外接长方体的尺寸作为所有非平稳特征对应的待模拟区域的尺寸；

在判断转移面中所有提取的特征至少有一个平稳地质特征时，则根据预设数值来设定所有非平稳特征对应的待模拟区域的尺寸。

5.如权利要求4所述的一种融合马尔可夫链和多点统计学的地质建模方法，其特征在于，所述根据转移面的转移概率和转移面的中线的距离权重计算得到待模拟区域中各模拟点的概率分布，具体包括：

根据转移面的转移概率和转移面的中线的距离权重，并通过以下公式计算得到待模拟区域中各模拟点的概率分布：

F_Tran(x,y)＝M_TP(i)×I(k)×Z

其中，F_Tran(x,y)表示模拟点(x,y)的概率分布，M_TP(i)表示在预设步长i下转移面的转移概率，Z表示模拟点到转移面的中线的距离权重，I(k)为用于标识不同地质结构的属性值，当地质结构为非平稳地质时则I(k)的属性值为1，当地质结构为平稳地质时则I(k)的属性值为0。

6.如权利要求1所述的一种融合马尔可夫链和多点统计学的地质建模方法，其特征在于，所述根据所述空间模式数据库对初始地质模型进行序贯模拟，得到初始三维地质模型，具体包括：

按照预设的顺序初始化模拟网格，从未模拟边缘区域空网格开始，从外向内逐层建模，并设定重叠区域的大小，以使得移动到每个位置的模型都有与已建模的区域重合的部分，对于重叠区域中的每个待模拟点，从模拟网格中提取出当前待模拟点的位置的模式，然后搜索空间模式数据库中相应的模式来进行序贯模拟，直到得到初始三维地质模型。

7.一种融合马尔可夫链和多点统计学的地质建模装置，其特征在于，包括：数据获取模块、三维网格图像生成模块、标识模块、非平稳地质结构三维模型构建模块、三维地质模型第一生成模块和三维地质模型第二生成模块；

所述数据获取模块，用于获取目标二维地质剖面；

所述三维网格图像生成模块，用于将目标二维地质剖面导入到三维模拟网格中，得到所述目标二维地质剖面在三维空间中表示的二维网格图像，并通过三维模拟网络将所述二维网格图像转换为三维网格图像；

所述标识模块，用于利用指示函数对三维网格图像中的非平稳地质结构和平稳地质结构进行分别标识；

所述非平稳地质结构三维模型构建模块，用于对于所述三维网格图像中的每一与三维网格图像平行的转移面，提取出转移面中的若干个非平稳地质特征和若干个平稳地质特征，并设置所有非平稳特征对应的待模拟区域的尺寸，基于马尔可夫链模型和反距离加权法得到在所述待模拟区域中各模拟点的概率分布，采用随机抽样方法对各模拟点进行抽样得到若干选定模拟点，基于若干选定模拟点在所述待模拟区域中构建出所有非平稳特征对应的非平稳地质结构三维模型；

所述三维地质模型第一生成模块，用于基于多点统计学算法，从三维网格图像中提取各地质结构的空间分布模式，并基于各空间分布模式建立地质对象空间模式数据库，根据所述空间模式数据库对初始地质模型进行序贯模拟，得到初始三维地质模型；

所述三维地质模型第二生成模块，用于将非平稳地质结构三维模型作为硬约束数据，基于多尺度策略并根据最大期望算法对初始三维地质模型进行迭代优化，直至累计迭代次数达到预设迭代次数时，得到最终的三维地质模型；

其中，所述基于马尔可夫链模型和反距离加权法得到在所述待模拟区域中各模拟点的概率分布，具体包括：

基于马尔可夫链模型，根据转移面中非平稳地质特征和平稳地质特征的个数之和得到转移计数矩阵，并根据转移计数矩阵得到各非平稳地质特征的转移概率矩阵；

根据各非平稳地质特征的转移概率矩阵，得到在转移面中沿着转移方向的起始非平稳特征和终止非平稳特征的转移概率；

根据如下公式计算出转移面的转移概率：

其中，M_TP(i)表示在预设步长i下转移面的转移概率，dis表示起始非平稳特征和终止非平稳特征之间的距离，TP₁和TP₂分别表示起始非平稳特征和终止非平稳特征的转移概率；

基于反距离加权法，根据各非平稳地质特征的形态拟合出转移面的中线，并遍历所述转移面内的所有点，计算每个点到非平稳特征中线的距离，选出最远的距离；

根据选出的最远的距离，通过以下公式计算待模拟区域中各模拟点到转移面的中线的距离权重：

其中，Z为模拟点到中线的距离权重，d(x,y)为模拟点到中线的距离，d_max为每个点到非平稳特征中线的距离中最远的距离；

根据转移面的转移概率和转移面的中线的距离权重计算得到待模拟区域中各模拟点的概率分布。