CN112053438A - 基于水平集的成矿构造深部推断建模方法 - Google Patents

Abstract

本发明涉及三维地质建模技术领域，尤其涉及一种基于水平集的成矿构造深部推断建模方法。该方法为：构建成矿构造初始模型并确定地质空间，进行单元划分获得多个立体单元，判断立体单元与成矿构造的空间关系获得初始物性参数；计算多个立体单元至成矿构造初始模型的最短距离，并根据立体单元所在区域赋予符号，获得初始水平集函数；根据地球物理异常数据及其正演模型，构建深部推断目标函数；以初始物性参数和初始水平集函数作为输入，利用水平集方法对目标函数进行迭代求解直至最优值，并获取最优水平集函数；提取最优水平集函数的零水平集，并根据零水平集获得更加准确可靠的成矿构造修正模型。

Description

基于水平集的成矿构造深部推断建模方法

技术领域

本发明涉及三维地质建模技术领域，尤其涉及一种基于水平集的成矿构造深部推断建模方法。

背景技术

三维地质建模理论和技术日益完善，对重现地质体形态特征具有重要作用。利用三维地质建模技术建立的成矿构造模型对于矿产资源预测与勘查工作意义重大。

目前，成矿构造的建模方法主要是基于地质勘查资料或者其他图件资料的显式建模方法，对于图件资料稀疏或者空白的成矿构造深部主要采用人工推断的方式建模，没有充分利用反映深部结构信息的地球物理异常数据，导致建立的成矿构造模型准确性和可靠性不高，对于深部构造特征和规律的启示十分有限，极大地影响了深部找矿工作的质量和成果。

发明内容

基于此，本发明针对上述的问题，提供了一种基于水平集的成矿构造深部推断建模方法，针对现有技术难以准确推断成矿构造深部结构特征的缺陷，旨在解决如何提升成矿构造模型准确性和可靠性的技术问题。

本发明实施例提供的一种基于水平集的成矿构造深部推断建模方法，具体包括：

构建成矿构造初始模型并确定地质空间，将所述地质空间进行单元划分获得多个立体单元，通过判断立体单元与成矿构造的空间关系获得初始物性参数；

计算所述多个立体单元至成矿构造初始模型的最短距离，并根据所述立体单元所在的区域赋予最短距离符号，获得初始水平集函数；

根据地球物理异常数据及其正演模型，构建深部推断的目标函数；

以所述初始物性参数和初始水平集函数作为输入，利用水平集方法对目标函数进行迭代求解至目标函数达到最优值，获得最优水平集函数；以及

提取所述最优水平集函数的零水平集，并根据所述零水平集获得成矿构造的修正模型。

更进一步地，所述通过判断立体单元与成矿构造的空间关系获得初始物性参数步骤，具体包括：

将所述立体单元点p_i(x_i,y_i,z_i)沿铅垂方向向成矿构造初始模型投影，获得投影点坐标

根据z_i与

的大小，赋予对应立体单元点的初始物性参数

所述s_int和s_ext分别表示所述立体单元点p_i位于成矿构造初始模型之上和位于成矿构造初始模型之下的岩体物性参数值。

更进一步地，计算多个立体单元至成矿构造初始模型的最短距离，并根据所述立体单元所在的区域赋予最短距离符号，获得初始水平集函数步骤具体包括：

过所述立体单元点向成矿构造初始模型作垂线，计算获得垂线段长度；

当所述立体单元区域在成矿构造初始模型之上时，最短距离值为垂线段长度，当所述立体单元区域在成矿构造初始模型之下时，最短距离值为垂线段长度的负值；

将上述获得的各立体单元的带符号的最短距离值构成水平集函数。

更进一步地，所述地球物理异常数据为表征不同岩石介质在物理性质(密度、磁性、电性、弹性、放射性等)方面存在的差异所引起的地球物理场上的局部变化值，包括重力异常数据、磁异常数据、电阻率异常数据、介质弹性异常数据、放射性异常数据等。

更进一步地，所述深部推断的目标函数为

所述y为地球物理异常测量值，为N×1向量，y＝{y_i,i＝1,2,…,N}；

所述h(s(p))为该地球物理异常对应的正演模型，将所述多个立体单元坐标和初始物性参数的数据转化为N个地球物理异常计算值；

L⁰(p)为初始水平集函数，L^t(p)为第t次迭代时的水平集函数，λ为与深度成负相关的函数。

更进一步地，以所述初始物性参数和初始水平集函数作为输入，利用水平集方法对目标函数进行迭代求解至目标函数达到最优值，获得最优水平集函数步骤，具体包括：

根据所述目标函数和所述初始物性参数分布计算迭代速度；

以所述迭代速度对水平集函数进行演化，获得水平集函数变化量以及变化后的水平集函数值；

根据所述变化后的水平集函数值更新物性参数，完成一次迭代过程；

将所述更新的物性参数进行所述物理异常数据的正演，获得更新的目标函数；

重复所述迭代过程和目标函数更新至目标函数达到最优值，获得最优水平集函数。

有益效果：

本发明的基于水平集的成矿构造深部推断建模方法通过地质空间单元划分与单元物性参数分配、水平集函数值计算、深部推断目标函数构建、目标函数的水平迭代求解和零水平集的提取及推断后模型构建，在特定的地球物理异常数据驱动下，将水平集表达下的成矿构造初始模型按照由目标函数确定的速度进行演化计算，当基于模型的正演数据计算值与地球物理异常数据值最为接近时，推断过程完成，得到更加准确可靠的成矿构造修正模型。

应当理解的是，以上的一般描述和后文的细节描述仅是示例性和解释性的，并不能限制本发明。

附图说明

图1为本发明实施例提供的一种基于水平集的成矿构造深部推断建模方法流程图；

图2为本发明实施例提供的地质空间单元划分示意图；

图3为本发明实施例提供的单元物性参数分配示意图；

图4为本发明实施例提供的目标函数的水平集迭代求解算法；

图5为本发明实施例提供的招平断裂中段成矿构造初始模型效果图；

图6为本发明实施例提供的招平断裂中段立体单元岩体密度参数分配效果图；

图7为本发明实施例提供的成矿构造初始模型水平集函数可视化效果图；

图8为本发明实施例提供的招平断裂中段重力异常值；

图9为本发明实施例提供的深部推断过程目标函数变化情况；

图10为本发明实施例提供的成矿构造修正模型效果图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。

如图1所示，在一个实施例中，提出了一种基于水平集的成矿构造深部推断建模方法，具体包括以下步骤：

步骤101，构建成矿构造初始模型并确定地质空间，将所述地质空间进行单元划分获得多个立体单元，通过判断立体单元与成矿构造的空间关系获得初始物性参数。

在本发明实施例中，地质空间是由成矿构造的初始模型确定的三维坐标范围，通常表示为Ω＝{(X_min～X_max),(Y_min～Y_max),(Z_min～Z_max)}。其中Ω代表地质空间，X_min,Y_min,Z_min分别为地质空间X,Y,Z方向的坐标最小值，X_max,Y_max,Z_max分别为地质空间X,Y,Z方向的坐标最大值。选取合适的尺寸(dx,dy,dz)对地质空间进行划分，其中dx,dy,dz分别代表X,Y,Z方向划分的单元尺度，则X,Y,Z方向的单元个数分别为nx＝(X_max-X_min)/dx；ny＝(Y_max-Y_min)/dy；nz＝(Z_max-Z_min)/dz，共计nx×ny×nz＝M个立体单元，并以中心点坐标代表各单元，即{p_i＝(x_i,y_i,z_i),i＝1,2,…,M}。其中p_i代表三维地质空间中的任意一点，x_i,y_i,z_i为其三维坐标。如图1所示。

在本发明实施例中，将点p_i沿Z轴向成矿构造模型投影，得到投影点及其坐标

比较点p_i和

的z值大小并赋予对应的岩体物性参数

其中s_int和s_ext分别表示点p_i位于成矿构造之上(上盘区域)和位于成矿构造之下(下盘区域)的岩体物性参数值，s(p_i)则表示点p_i分配的物性参数。如图2所示。

步骤S102，计算所述多个立体单元至成矿构造初始模型的最短距离，并根据所述立体单元所在的区域赋予最短距离符号，获得初始水平集函数。

在本发明实施例中，计算地质空间所有点到成矿构造模型的带符号的最短距离组成水平集函数值：过点p_i向成矿构造模型作垂线，记录垂足点及坐标

则点到成矿构造的最短距离为垂线段

的长度如公式(1)所示：

根据点p_i所在区域(上盘区域或下盘区域)为最短距离赋予符号，如公式(2)所示。

其中d(p_i)代表点p_i到成矿构造模型的带符号的最短距离。空间中所有点的符号最短距离值构成水平集函数L(p)＝{d(p_i),i＝1,2,…,M}。

步骤S103，根据地球物理异常数据及其正演模型，构建深部推断的目标函数。

在本发明实施例中，所述地球物理异常数据为表征不同岩石介质在物理性质(密度、磁性、电性、弹性、放射性等)方面存在的差异所引起的地球物理场上的局部变化值，包括重力异常数据、磁异常数据、电阻率异常数据、介质弹性异常数据、放射性异常数据等，根据特定的地球物理异常数据及其正演模型，构建通用的深部推断目标函数，形式如公式(3)所示。

目标函数第一项中y为地球物理异常测量值，为N×1向量，即y＝{y_i,i＝1,2,…,N}。而h(s(p))为该地球物理异常对应的正演模型，将M个包含地质空间点的坐标和物性参数的数据转化为N个地球物理异常计算值：

目标函数第二项中L⁰(p)代表初始水平集函数值，L^t(p)代表第t次迭代时的水平集函数值，λ则是与深度成负相关的函数。

步骤S104，以所述初始物性参数和初始水平集函数作为输入，利用水平集方法对目标函数进行迭代求解至目标函数达到最优值，获得最优水平集函数。

在本发明实施例中，目标函数的求解过程是不断地对水平集表示下的成矿构造模型进行演化直至目标函数达到最优值。为此，我们确定水平集函数演化所需的沿法线方向的速度如公式(4)所示。

根据演化速度v(p)确定的水平集函数的变化量如公式(5)所示。

其中δL(p)表示水平集函数的变化量，

为水平集函数梯度向量的模。综上所述，水平集求解目标函数的一次运算过程可以表达为以下步骤(假设当前为第t次迭代)：

根据已有的目标函数和物性参数分布计算当前的演化速度v^t(p)；

以当前的速度对水平集函数进行演化，获得其变化量δL^t(p)，并得到当前的水平集函数值L^t(p)＝L^t-1(p)+δL^t(p)。其中L^t-1(p)为第t-1次迭代获得的水平集函数值；

根据水平集函数值更新空间点的物性参数值，即

完成本次迭代过程；

基于当前点的物性参数分配值s^t(p)进行地球物理异常正演h(s^t(p))，并计算当前目标函数值O^t；

重复所述迭代过程和目标函数更新至目标函数达到最优值，获得最优水平集函数。

基于上述步骤，目标函数水平集迭代求解算法构建如图3所示。

步骤S105，提取所述最优水平集函数的零水平集，并根据所述零水平集获得成矿构造的修正模型。

在本发明实施例中，提取迭代推断结束后的水平集函数中值为零的数据集合组成零水平集{p＝(x,y,z)|L(p)＝0}，同时也代表成矿构造模型的散点集合。利用代表成矿构造的散点构建出成矿构造线框模型完成深部推断建模过程。

本方法在特定的地球物理异常数据驱动下，将水平集表达下的成矿构造初始模型按照由目标函数确定的速度进行演化计算，当基于模型的正演数据计算值与地球物理异常数据值最为接近时，推断过程完成，得到更加准确可靠的成矿构造修正模型。

实施例

下面结合附图对本发明具体实施例作进一步详细描述。以胶西北金矿集区招平断裂带中段成矿构造及重力异常数据为例，基于水平集的成矿构造深部推断与修正建模方法实施包含以下六个步骤：

(1)招平断裂中段成矿构造初始模型的构建和地质空间的划分

在地质剖面图中对成矿构造进行圈定得到其二维线串，导入到三维建模GOCAD软件中形成三维线串，依据勘探工程布置和成矿构造展布顺次连接各线串，用TIN(Triangulated Irregular Network，不规则三角网)模拟成矿构造界面，建立成矿构造的线框模型，称之为初始模型(如附图4)。为方便后续数值计算，将成矿构造线框模型导出为由三角形顶点及其三维坐标构成的数据文件。

根据成矿构造初始模型确定地质空间范围为Ω＝{(40528900～40535400),(4110100～4133400),(-3000～50)}，选取立体单元尺寸大小为(dx＝50,dy＝50,dz＝50)，由此得到X,Y,Z方向的单元个数分别为nx＝130,ny＝466,nz＝61，共计3695380个立体单元。在ParmAssign软件中计算立体单元坐标信息，并为其分配岩体密度参数，如附图5所示的密度参数分配效果图。应当理解的是立体单元尺寸的大小可以预先设定也可以有用户根据实际情况确定。

(2)成矿构造深部推断与修正建模软件ParmAssign和LSInve开发

ParmAssign软件基于Visual Studio 2015平台，采用C++语言进行开发，主要用于划分地质空间立体单元、计算立体单元坐标与分配立体单元物性参数，共包括4个模块：①数据读取模块；②立体单元坐标计算模型；③立体单元物性参数分配模块；④结果保存与输出模块。以下步骤中步骤(3)采用该软件完成。

LSInve软件基于MATLAB平台，采用matlab语言进行开发，主要用于生成水平集函数、迭代求解目标函数、提取零水平集等，共包括5个模块：①数据读取模块；②水平集函数生成模型；③地球物理正演计算模块；④水平集演化模块；⑤物性参数更新模块；⑥零水平集提取模块；⑦结果存储与输出模块。

(3)成矿构造初始模型水平集函数计算

利用LSInve软件中的LevelsetReset模块生成成矿构造初始模型的水平集函数值，将其导出后在GOCAD软件中进行可视化，如附图6。

(4)深部推断目标函数构建及求解

选取招平断裂中段重力异常为测量值(N＝33)，如图7所示。在此基础上构建如公式(3)所示的目标函数，并利用LSInve软件中的LevelsetInversAdjust模块对目标函数进行迭代求解。设置最大迭代次数为3000次，当迭代至2408次时目标函数收敛，演化结束，如图8所示，完成深部推断求解过程。

(5)成矿构造修正模型建立与可视化

利用LSInve软件提取推断修正后水平集函数的零水平集并输出，在三维建模软件GOCAD中建立其线框模型，如图9所示。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

本领域技术人员在考虑说明书及实践这里公开的发明后，将容易想到本公开的其它实施方案。本申请旨在涵盖本公开的任何变型、用途或者适应性变化，这些变型、用途或者适应性变化遵循本公开的一般性原理并包括本公开未公开的本技术领域中的公知常识或惯用技术手段。说明书和实施例仅被视为示例性的，本公开的真正范围和精神由权利要求指出。

应该理解的是，虽然本发明各实施例的流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，各实施例中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一非易失性计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。

以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。

Claims (6)

1.一种基于水平集的成矿构造深部推断建模方法，其特征在于，具体包括：

构建成矿构造初始模型并确定地质空间，将所述地质空间进行单元划分获得多个立体单元，通过判断立体单元与成矿构造的空间关系获得初始物性参数；

计算所述多个立体单元至成矿构造初始模型的最短距离，并根据所述立体单元所在的区域赋予最短距离符号，获得初始水平集函数；

根据地球物理异常数据及其正演模型，构建深部推断的目标函数；

以所述初始物性参数和初始水平集函数作为输入，利用水平集方法对目标函数进行迭代求解至目标函数达到最优值，获得最优水平集函数；以及

提取所述最优水平集函数的零水平集，并根据所述零水平集获得成矿构造的修正模型。

2.根据权利要求1所述的基于水平集的成矿构造深部推断建模方法，其特征在于，所述通过判断立体单元与成矿构造的空间关系获得初始物性参数步骤，具体包括：

将所述立体单元点p_i(x_i,y_i,z_i)沿铅垂方向向成矿构造初始模型投影，获得投影点坐标

根据z_i与

的大小，赋予对应立体单元点的初始物性参数

所述s_int和s_ext分别表示所述立体单元点p_i位于成矿构造初始模型之上和位于成矿构造初始模型之下的岩体物性参数值。

3.根据权利要求1所述的基于水平集的成矿构造深部推断建模方法，其特征在于，计算所述多个立体单元至成矿构造初始模型的最短距离，并根据所述立体单元所在的区域赋予最短距离符号，获得初始水平集函数步骤具体包括：

过所述立体单元点向成矿构造初始模型作垂线，计算获得垂线段长度；

当所述立体单元在成矿构造初始模型之上时，最短距离值为垂线段长度，当所述立体单元在成矿构造初始模型之下时，最短距离值为垂线段长度的负值；

将上述获得的各立体单元的带符号的最短距离值构成水平集函数。

4.根据权利要求1所述的基于水平集的成矿构造深部推断建模方法，其特征在于，所述地球物理异常数据为重力异常数据、磁异常数据、电阻率异常数据、介质弹性异常数据或放射性异常数据的任意一种。

5.根据权利要求1所述的基于水平集的成矿构造深部推断建模方法，其特征在于，所述深部推断的目标函数为：

所述y为地球物理异常测量值，为N×1向量，y＝{y_i,i＝1,2,…,N}；

所述h(s(p))为该地球物理异常对应的正演模型，将所述多个立体单元坐标和初始物性参数的数据转化为N个地球物理异常计算值；

L⁰(p)为初始水平集函数，L^t(p)为第t次迭代时的水平集函数，λ为与深度成负相关的函数。

6.根据权利要求1所述的基于水平集的成矿构造深部推断建模方法，其特征在于，以所述初始物性参数和初始水平集函数作为输入，利用水平集方法对目标函数进行迭代求解至目标函数达到最优值，获得最优水平集函数步骤，具体包括：

根据所述目标函数和所述初始物性参数分布计算迭代速度；

以所述迭代速度对水平集函数进行演化，获得水平集函数变化量以及变化后的水平集函数值；

根据所述变化后的水平集函数值更新物性参数，完成一次迭代过程；

将所述更新的物性参数进行所述物理异常数据的正演，获得更新的目标函数；

重复所述迭代过程和目标函数更新至目标函数达到最优值，获得最优水平集函数。

Images