CN113159321A - 重力约束下断裂面形态的贝叶斯推断方法 - Google Patents

Abstract

本发明提供了一种重力约束下断裂面形态的贝叶斯推断方法，包括：步骤1，将目标断裂面所属的三维地质空间进行离散化表示；步骤2，根据给定地质先验信息的先验函数计算目标断裂面的先验概率；步骤3，根据给定物性分布计算目标断裂面的地球物理正演值。本发明基于贝叶斯模型结合重力和地质先验信息推断断裂面形态，利用马尔科夫蒙特卡洛算法推断断裂面形态后验概率分布，结合信息熵算法可视化断裂面深部形态的不确定性空间分布，有效缓解单纯依赖地球物理反演或地质推断导致的深部形态不确定性，提升断裂面深部三维结构重建的准确性和有效性。

Description

重力约束下断裂面形态的贝叶斯推断方法

技术领域

本发明涉及三维地质建模技术领域，特别涉及一种重力约束下断裂面形态的贝叶斯推断方法。

背景技术

在矿产资源三维定量预测工作中，经常会涉及到研究区深部成矿部位数据稀少的问题。在有地质勘探工程支持的地区，可以获得一些高可靠钻孔数据和地质剖面数据，以进行近地表和浅部的结构推断，然而在深部区域勘探数据较少，且误差较大。这种误差较大的数据可能会对建模精度造成严重影响。传统基于地球物理反演的成矿预测多解性一直存在，也使得结果的真实性和有效性受到限制。

发明内容

本发明提供了一种重力约束下断裂面形态的贝叶斯推断方法，其目的是为了解决传统基于地球物理反演的成矿预测的多解性一直存在，导致预测结果的真实性和有效性受到限制的问题。

为了达到上述目的，本发明的实施例提供了一种重力约束下断裂面形态的贝叶斯推断方法，包括：

步骤1，将目标断裂面所属的三维地质空间进行离散化表示；

步骤2，根据给定地质先验信息的先验函数计算目标断裂面的先验概率；

步骤3，根据给定物性分布计算目标断裂面的地球物理正演值；

步骤4，根据计算出的目标断裂面的地球物理正演值和给定地球物理观测数据集合计算目标断裂面的实测地球物理数据的概率密度分布；

步骤5，根据计算出的目标断裂面的实测地球物理数据的概率密度分布和给定目标断裂面形态下的物性分布差异计算目标断裂面的似然函数；

步骤6，根据计算出的目标断裂面的先验概率和计算出的目标断裂面的似然函数构建贝叶斯模型并计算贝叶斯模型下的目标断裂面的后验概率；

步骤7，对目标断裂面的轮廓线上的关键点进行提取，建立以各关键点为结点，关键点间空间临近度为边的图结构；在图结构和扰动幅度符合高斯分布的约束下对关键点进行随机扰动，对扰动后的关键点重构目标断裂面轮廓线以生成扰动断裂面模型；

步骤8，计算扰动后目标断裂面的后验概率，根据前一个目标断裂面的后验概率和当前扰动后目标断裂面的后验概率构建马尔科夫蒙特卡洛算法接受率，根据马尔科夫蒙特卡洛算法接受率判断是否接受对当前扰动断裂面模型的采样；

步骤9，重复执行步骤7和步骤8，直到得到目标采样数量的扰动断裂面模型时停止，获得多个对目标断裂面形态的采样，实现对目标断裂面形态的贝叶斯推断；

步骤10，获得多个扰动断裂面模型后利用信息熵算法对目标断裂面进行不确定性的可视化分析。

其中，所述步骤1具体包括：

对目标断裂面所属的三维地质空间

进行离散化表示，将目标断裂面所属的三维地质空间

分割为多个立体单元，在目标断裂面所属的三维地质空间

中，利用隐函数F:

对目标断裂面形态进行表达，使得目标断裂面所属的三维地质空间中任意一点c，满足

采用隐函数建立每个立体单元中心点c＝(x,y,z)相对于目标断裂面F的位置关系，如下所示：

其中，所述步骤2具体包括：

根据给定地质先验信息的先验函数计算目标断裂面的先验概率，如下所示：

其中，Z表示正则化项，E(F)表示地质先验信息的先验函数。

其中，所述步骤3具体包括：

根据给定物性分布Θ，计算目标断裂面的地球物理正演值，如下所示：

x_i＝x-ξ_i (4)

y_j＝y-η_j (5)

z_k＝z-ζ_k (6)

μ_ijk＝(-1)ⁱ(-1)^j(-1)^k (8)

其中，G表示重力常量，ρ表示剩余密度，(x,y,z)表示地面观测点坐标，(ξ_i,η_j,ζ_k)表示每个立体单元角点坐标，Θ表示给定物性分布，μ_ijk表示与立体单元角点位置相关的符号参数，r_ijk表示每个立体单元角点与目标点之间的距离，i、j、k分别表示立体单元角点的位置。

其中，所述步骤4具体包括：

根据目标断裂面的地球物理正演值和给定地球物理观测数据集合计算目标断裂面的实测地球物理数据的概率密度分布，如下所示：

其中，P(ε|Θ)表示实测地球物理数据的概率密度分布，ε表示地球物理观测数据集合，ε＝{ε₁,ε₂,...,ε_N}，fwd(Θ)表示目标断裂面的地球物理正演值，ε_n表示第n个观测点的地球物理观测值，σ表示常数。

其中，所述步骤5具体包括：

根据目标断裂面的实测地球物理数据的概率密度分布和目标断裂面形态下的物性分布差异，计算目标断裂面的似然函数，如下所示：

P(ε|F)＝∫P(ε|Θ)P(Θ|F)dΘ (10)

其中，P(ε|F)表示似然函数，P(Θ|F)表示目标断裂面形态下的物性分布差异，P(Θ|F)根据地球物理观测数据集合ε＝{ε₁,ε₂,...,ε_N}采用核密度估计构建得到。

其中，所述步骤6具体包括：

给定地球物理观测数据集合ε＝{ε₁,ε₂,...,ε_N}，利用地球物理观测数据集合ε推断目标断裂面，在贝叶斯模型框架下，根据目标断裂面的先验概率和目标断裂面的似然函数计算目标断裂面的后验概率，如下所示：

其中，P(F)表示目标断裂面的先验概率，P(F)表达目标断裂面形态和产状的先验信息，P(ε)表示地球物理观测数据集合的概率分布，为常数。

其中，所述步骤7具体包括：

步骤71，获得目标断裂面的轮廓线；

步骤72，在目标断裂面的轮廓线上选取关键点p_a，其中，P＝{p₁,p₂,p₃,...,p_A}表示关键点的集合，p_a表示关键点，a＝1,2,3,...,A；

步骤73，对每个关键点p_a，选取每个关键点p_a周边的d个最邻近的关键点p_b，构建关键点p_a分别与邻近关键点p_b的边连接，形成以关键点为结点，邻近关系为边的图结构；

步骤74，根据关键点在x,y方向上的扰动符合高斯分布，对关键点的扰动大小进行随机采样，根据关键点的扰动大小对选取的关键点集合P进行扰动，得到扰动后关键点位置集合

步骤75，通过对扰动后关键点位置集合进行优化，得到保持目标断裂面整体形态的关键点位置

具体为求解目标函数极小化问题，如下所示：

其中，p'_a表示优化后关键点的位置，L表示Laplace矩阵，δ_a为上一次扰动后第a个点的Laplace坐标，λ表示对扰动项的权重系数；

步骤76，基于扰动优化后的关键点集P'，通过线性插值得到目标断裂面各轮廓上的其他顶点，根据扰动优化后的关键点集P'和其他顶点重建目标断裂面轮廓线；

步骤77，根据重建的目标断裂面轮廓线，采用隐式建模方式，构建目标断裂面形态扰动后的扰动断裂面模型。

其中，所述步骤8具体包括：

根据前一个目标断裂面的后验概率和当前扰动后目标断裂面的后验概率构建马尔科夫蒙特卡洛算法接受率，如下所示：

其中，α表示接受率，F表示前一个目标断裂面，F*表示扰动后目标断裂面，P(F|ε)表示前一个目标断裂面的后验概率，P(F^*|ε)表示扰动后目标断裂面的后验概率，q(F*|F)表示扰动后目标断裂面采样与前一次目标断裂面采样的转移概率，满足q(F*|F)～N(F,σ)，q(F*|F)表示前一次目标断裂面采样与扰动后目标断裂面采样的转移概率，q(F*|F)＝q(F|F*)。

其中，所述步骤9具体包括：

通过信息熵算法对目标断裂面进行不确定性的可视化分析，对每个立体单元s_l考虑相对于断裂面上盘F_up和断裂面下盘F_down的位置变化，计算每个立体单元的信息熵值，如下所示：

其中，H(s_l)表示每个立体单元的信息熵值，s_l表示立体单元，

表示第l个立体单元位于目标断裂面上盘的概率，

表示第l个立体单元位于目标断裂面下盘的概率，

和

的值通过马尔科夫蒙特卡洛算法采样获得的多个扰动断裂面模型进行计算，通过计算出的每个立体单元的信息熵值度量和可视化目标断裂面所属的三维地质空间的不确定性。

本发明的上述方案有如下的有益效果：

本发明的上述实施例所述的重力约束下断裂面形态的贝叶斯推断方法，在贝叶斯模型框架下，融合地质先验知识和地球物理数据，利用断裂面上下盘的物性差异性，实现对断裂面深部形态的推断及其不确定性度量，有助于提升断裂面三维结构重建的准确性和有效性。

附图说明

图1为本发明的流程图；

图2为本发明的关键点扰动示意图；

图3为本发明的招平断裂带仅地质约束下的不确定性结果示意图，(a)为全空间信息熵分布示意图；(b)为空间非零信息熵分布示意图；(c)为断裂面上盘非零信息熵分布示意图；(d)为断裂面下盘非零信息熵分布示意图；

图4为本发明的招平断裂带地质约束与重力约束下的不确定性结果示意图，(a)为全空间信息熵分布示意图；(b)为空间非零信息熵分布示意图；(c)为相对于仅地质约束下断裂面信息熵增加的部分示意图；(d)为相对于仅地质约束下断裂面信息熵减少的部分示意图；

图5为本发明的招平断裂带地质约束、重力约束及法向量约束下的不确定性结果示意图，(a)为全空间信息熵分布示意图；(b)为空间非零信息熵分布示意图；(c)为相对于地质约束和重力约束下断裂面信息熵增加的部分示意图；(d)为相对于地质约束和重力约束下断裂面信息熵减少的部分示意图。

具体实施方式

为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。

本发明针对现有的传统基于地球物理反演的成矿预测的多解性一直存在，导致预测结果的真实性和有效性受到限制的问题，提供了一种重力约束下断裂面形态的贝叶斯推断方法。

如图1至图5所示，本发明的实施例提供了一种重力约束下断裂面形态的贝叶斯推断方法，包括：步骤1，将目标断裂面所属的三维地质空间进行离散化表示；步骤2，根据给定地质先验信息的先验函数计算目标断裂面的先验概率；步骤3，根据给定物性分布计算目标断裂面的地球物理正演值；步骤4，根据计算出的目标断裂面的地球物理正演值和给定地球物理观测数据集合计算目标断裂面的实测地球物理数据的概率密度分布；步骤5，根据计算出的目标断裂面的实测地球物理数据的概率密度分布和给定目标断裂面形态下的物性分布差异计算目标断裂面的似然函数；步骤6，根据计算出的目标断裂面的先验概率和计算出的目标断裂面的似然函数构建贝叶斯模型并计算贝叶斯模型下的目标断裂面的后验概率；步骤7，对目标断裂面的轮廓线上的关键点进行提取，建立以各关键点为结点，关键点间空间临近度为边的图结构；在图结构和扰动幅度符合高斯分布的约束下对关键点进行随机扰动，对扰动后的关键点重构目标断裂面轮廓线以生成扰动断裂面模型；步骤8，计算扰动后目标断裂面的后验概率，根据前一个目标断裂面的后验概率和当前扰动后目标断裂面的后验概率构建马尔科夫蒙特卡洛算法接受率，根据马尔科夫蒙特卡洛算法接受率判断是否接受对当前扰动断裂面模型的采样；步骤9，重复执行步骤7和步骤8，直到得到目标采样数量的扰动断裂面模型时停止，获得多个对目标断裂面形态的采样，实现对目标断裂面形态的贝叶斯推断；步骤10，获得多个扰动断裂面模型后利用信息熵算法对目标断裂面进行不确定性的可视化分析。

其中，所述步骤1具体包括：对目标断裂面所属的三维地质空间

进行离散化表示，将目标断裂面所属的三维地质空间

分割为多个立体单元，在目标断裂面所属的三维地质空间中任意一点c，满足

采用隐函数建立每个立体单元中心点c＝(x,y,z)相对于目标断裂面F的位置关系，如下所示：

其中，所述步骤2具体包括：根据给定地质先验信息的先验函数计算目标断裂面的先验概率，如下所示：

其中，Z表示正则化项，E(F)表示地质先验信息的先验函数。

其中，所述步骤3具体包括：根据给定物性分布Θ，计算目标断裂面的地球物理正演值，如下所示：

x_i＝x-ξ_i (4)

y_j＝y-η_j (5)

z_k＝z-ζ_k (6)

μ_ijk＝(-1)ⁱ(-1)^j(-1)^k (8)

其中，G表示重力常量，ρ表示剩余密度，(x,y,z)表示地面观测点坐标，(ξ_i,η_j,ζ_k)表示每个立体单元角点坐标，Θ表示给定物性分布，μ_ijk表示与立体单元角点位置相关的符号参数，r_ijk表示每个立体单元角点与目标点之间的距离，i、j、k分别表示立体单元角点的位置。

其中，所述步骤4具体包括：根据目标断裂面的地球物理正演值和给定地球物理观测数据集合计算目标断裂面的实测地球物理数据的概率密度分布，如下所示：

其中，P(ε|Θ)表示实测地球物理数据的概率密度分布，ε表示地球物理观测数据集合，ε＝{ε₁,ε₂,...,ε_N}，fwd(Θ)表示目标断裂面的地球物理正演值，ε_n表示第n个观测点的地球物理观测值，σ表示常数。

其中，所述步骤5具体包括：根据目标断裂面的实测地球物理数据的概率密度分布和目标断裂面形态下的物性分布差异，计算目标断裂面的似然函数，如下所示：

P(ε|F)＝∫P(ε|Θ)P(Θ|F)dΘ (10)

其中，P(ε|F)表示似然函数，P(Θ|F)表示目标断裂面形态下的物性分布差异，P(Θ|F)根据地球物理观测数据集合ε＝{ε₁,ε₂,...,ε_N}采用核密度估计构建得到。

本发明的上述实施例所述的重力约束下断裂面形态的贝叶斯推断方法，由于断裂面上下盘可能具有不同的岩性，因此剩余密度并非一个恒定的值，并且同一种岩性内部密度也有所变化，因此在计算目标断裂面的地球物理正演值时对断裂面研究区进行分区最优密度反演，根据每个区的最佳密度反演求得目标断裂面的地球物理正演值，将目标断裂面的地球物理正演值与给定的地球物理观测数据集合进行比较，获得目标断裂面形态导致的物性分布差异从而进行似然函数的计算。

其中，所述步骤6具体包括：给定地球物理观测数据集合ε＝{ε₁,ε₂,...,ε_N}，利用地球物理观测数据集合ε推断目标断裂面，在贝叶斯模型框架下，根据目标断裂面的先验概率和目标断裂面的似然函数计算目标断裂面的后验概率，如下所示：

其中，P(F)表示目标断裂面的先验概率，P(F)表达目标断裂面形态和产状的先验信息，P(ε)表示地球物理观测数据集合的概率分布，为常数。

本发明的上述实施例所述的重力约束下断裂面形态的贝叶斯推断方法，目标断裂面的似然函数P(ε|F)反映基于目标断裂面F的地球物理正演值与地球物理观测数据集合的一致性程度。

其中，所述步骤7具体包括：步骤71，获得目标断裂面的轮廓线；步骤72，在目标断裂面的轮廓线上选取关键点p_a，其中，P＝{p₁,p₂,p₃,...,p_A}表示关键点的集合，p_a表示关键点，a＝1,2,3,...,A；步骤73，对每个关键点p_a，选取每个关键点p_a周边的d个最邻近的关键点p_b，构建关键点p_a分别与邻近关键点p_b的边连接，形成以关键点为结点，邻近关系为边的图结构；步骤74，根据关键点在x,y方向上的扰动符合高斯分布，对关键点的扰动大小进行随机采样，根据关键点的扰动大小对选取的关键点集合P进行扰动，得到扰动后关键点位置集合

步骤75，通过对扰动后关键点位置集合进行优化，得到保持目标断裂面整体形态的关键点位置

具体为求解目标函数极小化问题，如下所示：

其中，p'_a表示优化后关键点的位置，L表示Laplace矩阵，δ_a为上一次扰动后第a个点的Laplace坐标，λ表示对扰动项的权重系数；步骤76，基于扰动优化后的关键点集P'，通过线性插值得到目标断裂面各轮廓上的其他顶点，根据扰动优化后的关键点集P'和其他顶点重建目标断裂面轮廓线；步骤77，根据重建的目标断裂面轮廓线，采用隐式建模方式，构建目标断裂面形态扰动后的扰动断裂面模型。

其中，所述步骤8具体包括：根据前一个目标断裂面的后验概率和当前扰动后目标断裂面的后验概率构建马尔科夫蒙特卡洛算法接受率，如下所示：

其中，α表示接受率，F表示前一个目标断裂面，F*表示扰动后目标断裂面，P(F|ε)表示前一个目标断裂面的后验概率，P(F^*|ε)表示扰动后目标断裂面的后验概率，q(F*|F)表示扰动后目标断裂面采样与前一次目标断裂面采样的转移概率，满足q(F*|F)～N(F,σ)，q(F*|F)表示前一次目标断裂面采样与扰动后目标断裂面采样的转移概率，q(F*|F)＝q(F|F*)。

本发明的上述实施例所述的重力约束下断裂面形态的贝叶斯推断方法，由于马尔科夫链需要经过多次状态转移才能达到稳定，因此对之前的采样进行舍弃，使采样结果更加接近真实结果。

其中，所述步骤9具体包括：通过信息熵算法对目标断裂面进行不确定性的可视化分析，对每个立体单元s_l考虑相对于断裂面上盘F_up和断裂面下盘F_down的位置变化，计算每个立体单元的信息熵值，如下所示：

其中，H(s_l)表示每个立体单元的信息熵值，s_l表示立体单元，

表示第l个立体单元位于目标断裂面上盘的概率，

表示第l个立体单元位于目标断裂面下盘的概率，

和

的值通过马尔科夫蒙特卡洛算法采样获得的多个扰动断裂面模型进行计算，通过计算出的每个立体单元的信息熵值度量和可视化目标断裂面所属的三维地质空间的不确定性。

本发明的上述实施例所述的重力约束下断裂面形态的贝叶斯推断方法，通过每个立体单元的信息熵值，以表示其不确定性，最终实现对目标断裂面所属的三维地质空间不确定性的可视化。

本发明上述实施例所述的重力约束下断裂面形态的贝叶斯推断方法，地质先验信息的先验函数包括：目标断裂面形态一致性先验函数，目标断裂面法向量先验函数和目标断裂面平滑先验函数；

目标断裂面形态一致性先验函数E₁(F^t)表示为：

E₁(F^t)＝λ(z)||F^t-F₀||² (15)

其中，F₀表示初始目标断裂面计算得到的隐函数数据集，F^t表示目标断裂面第t次扰动计算得到的隐函数数据集，λ(z)表示与深度相关的权重项，控制高斯分布扰动断裂面模型在深度上的变化程度，z表示对应深度；

在产状先验信息中引入法向量约束，令目标断裂面为F(x,y,z)＝0的三维曲面上某点的法向量为该点的梯度方向，以

表示，目标断裂面法向量先验函数E₂(F^t)，如下所示：

其中，C表示已知梯度数据集，

表示对应已知梯度数据集的每次迭代后的梯度值，n_i表示已知梯度数据集每一点的梯度，s_i表示目标断裂面上的第i个立体单元，λ表示系数，为一常量；

加入拉普拉斯算子维持每次迭代的目标断裂面的平滑特征，目标断裂面平滑先验函数E₃(F^t)，如下所示：

E₃(F^t)＝||ΔF^t||² (17)

其中，Δ表示拉普拉斯算子：

其中，ΔF^t表示加入拉普拉斯算子后的每次迭代的目标断裂面，

表示某一方向上的偏导数；

地质先验信息的先验函数，如下所示：

E(F)＝E₁(F^t)+E₂(F^t)+E₃(F^t) (19)。

本发明上述实施例所述的重力约束下断裂面形态的贝叶斯推断方法，以山东招平主断裂面形态推断为例说明，在此例中，深部断裂面三维结构具有很强的不确定性，为获得更精细的地质三维模型，结合地质认知和地球物理重力数据，利用所述重力约束下断裂面形态的贝叶斯推断方法实现断裂面研究区的不确定性预测建模，具体实施方式按以下步骤描述：1.根据研究区范围确定离散化立体单元中心点坐标，结合工程控制线串对每个立体单元对应上下盘的初始位置进行隐式表达；2.构建贝叶斯框架模型融合地质和地球物理数据信息，利用重力数据计算似然函数，将地质信息作为贝叶斯框架模型的先验约束，以此作为后续马尔科夫蒙特卡洛算法采样的依据；3.在断裂面轮廓线上提取关键点，基于轮廓线关键点，建立关键点之间的图结构，基于图结构信息，构建轮廓线采样目标函数和极小化目标函数得到的线性方程组，4.基于轮廓线采样方法和隐式建模方法，在马尔科夫蒙特卡洛算法框架下开展对断裂面形态的采样，得到断裂面形态采样结果，结合成矿后验概率得到采样结果接受率α，实现对断裂面形态后验概率的近似推断；5.基于断裂面形态采样结果，利用信息熵算法，可视化招平主断裂面形态的不确定性分布，得到不确定性推断结果，不确定性推断结果如图3、图4和图5所示。

本发明上述实施例所述的重力约束下断裂面形态的贝叶斯推断方法，基于贝叶斯框架结合重力和地质先验信息推断断裂面形态，利用马尔科夫蒙特卡洛算法推断断裂面形态后验概率分布，最终结合信息熵可视化断裂面深部形态的不确定性空间分布，从而在有效缓解单纯依赖地球物理反演或地质推断导致的深部形态不确定性，提升断裂面深部三维结构重建的准确性和有效性。

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

Claims (10)

1.一种重力约束下断裂面形态的贝叶斯推断方法，其特征在于，包括：

步骤1，将目标断裂面所属的三维地质空间进行离散化表示；

步骤2，根据给定地质先验信息的先验函数计算目标断裂面的先验概率；

步骤3，根据给定物性分布计算目标断裂面的地球物理正演值；

步骤4，根据计算出的目标断裂面的地球物理正演值和给定地球物理观测数据集合计算目标断裂面的实测地球物理数据的概率密度分布；

步骤5，根据计算出的目标断裂面的实测地球物理数据的概率密度分布和给定目标断裂面形态下的物性分布差异计算目标断裂面的似然函数；

步骤6，根据计算出的目标断裂面的先验概率和计算出的目标断裂面的似然函数构建贝叶斯模型并计算贝叶斯模型下的目标断裂面的后验概率；

步骤7，对目标断裂面的轮廓线上的关键点进行提取，建立以各关键点为结点，关键点间空间临近度为边的图结构；在图结构和扰动幅度符合高斯分布的约束下对关键点进行随机扰动，对扰动后的关键点重构目标断裂面轮廓线以生成扰动断裂面模型；

步骤8，计算扰动后目标断裂面的后验概率，根据前一个目标断裂面的后验概率和当前扰动后目标断裂面的后验概率构建马尔科夫蒙特卡洛算法接受率，根据马尔科夫蒙特卡洛算法接受率判断是否接受对当前扰动断裂面模型的采样；

步骤9，重复执行步骤7和步骤8，直到得到目标采样数量的扰动断裂面模型时停止，获得多个对目标断裂面形态的采样，实现对目标断裂面形态的贝叶斯推断；

步骤10，获得多个扰动断裂面模型后利用信息熵算法对目标断裂面进行不确定性的可视化分析。

2.根据权利要求1所述的重力约束下断裂面形态的贝叶斯推断方法，其特征在于，所述步骤1具体包括：

对目标断裂面所属的三维地质空间

进行离散化表示，将目标断裂面所属的三维地质空间

分割为多个立体单元，在目标断裂面所属的三维地质空间

中，利用隐函数

对目标断裂面形态进行表达，使得目标断裂面所属的三维地质空间中任意一点c，满足

采用隐函数建立每个立体单元中心点c＝(x,y,z)相对于目标断裂面F的位置关系，如下所示：

3.根据权利要求2所述的重力约束下断裂面形态的贝叶斯推断方法，其特征在于，所述步骤2具体包括：

根据给定地质先验信息的先验函数计算目标断裂面的先验概率，如下所示：

其中，Z表示正则化项，E(F)表示地质先验信息的先验函数。

4.根据权利要求3所述的重力约束下断裂面形态的贝叶斯推断方法，其特征在于，所述步骤3具体包括：

根据给定物性分布Θ，计算目标断裂面的地球物理正演值，如下所示：

x_i＝x-ξ_i (4)

y_j＝y-η_j (5)

z_k＝z-ζ_k (6)

μ_ijk＝(-1)ⁱ(-1)^j(-1)^k (8)

其中，G表示重力常量，ρ表示剩余密度，(x,y,z)表示地面观测点坐标，(ξ_i,η_j,ζ_k)表示每个立体单元角点坐标，Θ表示给定物性分布，μ_ijk表示与立体单元角点位置相关的符号参数，r_ijk表示每个立体单元角点与目标点之间的距离，i、j、k分别表示立体单元角点的位置。

5.根据权利要求4所述的重力约束下断裂面形态的贝叶斯推断方法，其特征在于，所述步骤4具体包括：

根据目标断裂面的地球物理正演值和给定地球物理观测数据集合计算目标断裂面的实测地球物理数据的概率密度分布，如下所示：

其中，P(ε|Θ)表示实测地球物理数据的概率密度分布，ε表示地球物理观测数据集合，ε＝{ε₁,ε₂,...,ε_N}，fwd(Θ)表示目标断裂面的地球物理正演值，ε_n表示第n个观测点的地球物理观测值，σ表示常数。

6.根据权利要求5所述的重力约束下断裂面形态的贝叶斯推断方法，其特征在于，所述步骤5具体包括：

根据目标断裂面的实测地球物理数据的概率密度分布和目标断裂面形态下的物性分布差异，计算目标断裂面的似然函数，如下所示：

P(ε|F)＝∫P(ε|Θ)P(Θ|F)dΘ (10)

其中，P(ε|F)表示似然函数，P(Θ|F)表示目标断裂面形态下的物性分布差异，P(Θ|F)根据地球物理观测数据集合ε＝{ε₁,ε₂,...,ε_N}采用核密度估计构建得到。

7.根据权利要求6所述的重力约束下断裂面形态的贝叶斯推断方法，其特征在于，所述步骤6具体包括：

给定地球物理观测数据集合ε＝{ε₁,ε₂,...,ε_N}，利用地球物理观测数据集合ε推断目标断裂面，在贝叶斯模型框架下，根据目标断裂面的先验概率和目标断裂面的似然函数计算目标断裂面的后验概率，如下所示：

其中，P(F)表示目标断裂面的先验概率，P(F)表达目标断裂面形态和产状的先验信息，P(ε)表示地球物理观测数据集合的概率分布，为常数。

8.根据权利要求7所述的重力约束下断裂面形态的贝叶斯推断方法，其特征在于，所述步骤7具体包括：

步骤71，获得目标断裂面的轮廓线；

步骤72，在目标断裂面的轮廓线上选取关键点p_a，其中，P＝{p₁,p₂,p₃,...,p_A}表示关键点的集合，p_a表示关键点，a＝1,2,3,...,A；

步骤73，对每个关键点p_a，选取每个关键点p_a周边的d个最邻近的关键点p_b，构建关键点p_a分别与邻近关键点p_b的边连接，形成以关键点为结点，邻近关系为边的图结构；

步骤74，根据关键点在x,y方向上的扰动符合高斯分布，对关键点的扰动大小进行随机采样，根据关键点的扰动大小对选取的关键点集合P进行扰动，得到扰动后关键点位置集合

步骤75，通过对扰动后关键点位置集合进行优化，得到保持目标断裂面整体形态的关键点位置

具体为求解目标函数极小化问题，如下所示：

其中，p'_a表示优化后关键点的位置，L表示Laplace矩阵，δ_a为上一次扰动后第a个点的Laplace坐标，λ表示对扰动项的权重系数；

步骤76，基于扰动优化后的关键点集P'，通过线性插值得到目标断裂面各轮廓上的其他顶点，根据扰动优化后的关键点集P'和其他顶点重建目标断裂面轮廓线；

步骤77，根据重建的目标断裂面轮廓线，采用隐式建模方式，构建目标断裂面形态扰动后的扰动断裂面模型。

9.根据权利要求8所述的重力约束下断裂面形态的贝叶斯推断方法，其特征在于，所述步骤8具体包括：

根据前一个目标断裂面的后验概率和当前扰动后目标断裂面的后验概率构建马尔科夫蒙特卡洛算法接受率，如下所示：

其中，α表示接受率，F表示前一个目标断裂面，F*表示扰动后目标断裂面，P(F|ε)表示前一个目标断裂面的后验概率，P(F^*|ε)表示扰动后目标断裂面的后验概率，q(F*|F)表示扰动后目标断裂面采样与前一次目标断裂面采样的转移概率，满足q(F*|F)～N(F,σ)，q(F*|F)表示前一次目标断裂面采样与扰动后目标断裂面采样的转移概率，q(F*|F)＝q(F|F*)。

10.根据权利要求9所述的重力约束下断裂面形态的贝叶斯推断方法，其特征在于，所述步骤10具体包括：

通过信息熵算法对目标断裂面进行不确定性的可视化分析，对每个立体单元s_l考虑相对于断裂面上盘F_up和断裂面下盘F_down的位置变化，计算每个立体单元的信息熵值，如下所示：

其中，H(s_l)表示每个立体单元的信息熵值，s_l表示立体单元，

表示第l个立体单元位于目标断裂面上盘的概率，

表示第l个立体单元位于目标断裂面下盘的概率，

和

的值通过马尔科夫蒙特卡洛算法采样获得的多个扰动断裂面模型进行计算，通过计算出的每个立体单元的信息熵值度量和可视化目标断裂面所属的三维地质空间的不确定性。

Images