CN112307640A - 基于自然单元-无限元的三维多源自然电位数值模拟方法 - Google Patents

Abstract

本发明涉及自然电位数值模拟方法技术领域，尤其涉及基于自然单元‑无限元的三维多源自然电位数值模拟方法。包括：S1、构建三维多源自然电位地电模型；S2、构建模型的微分方程和积分方程；S3、基于Laplace插值，构建三维自然单元插值函数及其导数的基本方程；S4、构建三维单向映射无限元插值函数及其导数的基本方程；S5、基于步骤S2‑S4，对模型的全域积分，得到各积分单元的刚度矩阵；S6、将各积分单元的刚度矩阵按节点编号累加，得到总刚度矩阵；S7、构建源向量；S8、基于总刚度矩阵和源向量，求解大型稀疏方程组得到自然电位分布。解决了多源复杂模型中复杂异常体剖分困难、截断边界问题处理繁琐等技术问题。

Description

基于自然单元-无限元的三维多源自然电位数值模拟方法

技术领域

本发明涉及自然电位数值模拟方法技术领域，尤其涉及一种基于自然单元-无限元的三维多源自然电位数值模拟方法。

背景技术

在地球物理勘探领域中，自然电位法因其使用的仪器设备简便、工作效率高、无二次污染、成本低廉等优点得到了广泛应用。尤其是自然电位法对地下溶质运移、氧化还原反应、自然极化、热扩散等过程具有特定的电信号响应，已被广泛应用于矿产油气勘探、污染羽流的圈定、地热资源探测以及土壤参数估计等领域。

数据反演和解释是解决地球物理勘探实际问题的关键，而数值模拟算法是开展反演解释的有效数值工具。有针对性地研究适用于自然电位数值模拟的数值算法有助于提升自然电位法的应用效果。

自然电位场源具有多源性、区域性、动态性分布特点，常规数值模拟算法如有限单元法、有限差分法等受限于网格剖分，难以有效开展复杂模型的数值模拟工作。自然单元法依托于求解域内离散节点的Voronoi结构，采用无网格方式构造插值函数，整个模型空间不受网格单元约束，能有效剖分极端复杂的几何目标体，能有效贴合起伏地形。而对于多源性，有限单元法等常规算法通常依托场源位置对边界施加特定的边界条件，但在多源模型中场源位置只能近似选择多个场源的几何中心或其他可近似替代的位置，影响了计算精度。

发明内容

(一)要解决的技术问题

鉴于现有技术的上述缺点、不足，本发明提供一种基于自然单元-无限元的三维多源自然电位数值模拟方法，其解决了多源复杂模型中复杂异常体剖分困难、截断边界问题处理繁琐等技术问题。

(二)技术方案

为了达到上述目的，本发明采用的主要技术方案包括：

本发明实施例提供一种基于自然单元-无限元的三维多源自然电位数值模拟方法，包括以下步骤：

S1、构建三维多源自然电位地电模型，三维多源自然电位地电模型由多个积分单元组成；

三维多源自然电位地电模型包括自然单元区域和无限单元区域；

无限单元区域包括多个单向映射的无限单元；

S2、构建三维多源自然电位地电模型的微分方程，基于微分方程得到三维多源自然电位地电模型的积分方程；

S3、基于Laplace插值，构建自然单元区域的插值函数及其导数的基本方程；

S4、基于单向映射的无限单元，构建无限单元区域的插值函数及其导数的基本方程；

S5、基于步骤S2-S4，对三维多源自然电位地电模型的全域进行积分，得到各积分单元的刚度矩阵；

S6、将各积分单元的刚度矩阵按节点编号累加，得到总刚度矩阵；

S7、将场源按分布区域、幅值、正负性赋值于相应的空间节点处，构建源向量；

S8、基于总刚度矩阵和源向量，求解大型稀疏方程，得到三维多源自然电位地电模型的自然电位分布。

本发明实施例提出的基于自然单元-无限元的三维多源自然电位数值模拟方法，在自然单元区域边界处耦合无限单元，能够有效解决截断边界问题，提高数值模拟精度。

可选地，步骤S1中，依据模型尺度、异常体几何参数和起伏地形，在自然单元区域布设离散自然节点。

可选地，自然单元区域内分布有自然节点，无限单元区域分布有无限单元，自然单元区域与无限单元区域的耦合处设置有连接点，连接点既是自然节点也是无限单元的节点。

可选地，步骤S2中，三维多源自然电位地电模型的微分方程满足以下公式：

式中，σ为电导率，V为待求自然电位，j_s为场源；

三维多源自然电位地电模型基于自然单元-无限元的积分方程满足以下公式：

式中，Ω为积分单元区域，e为积分单元，V_e为积分单元中各自然邻点的电位值或各无限单元节点的电位值，

为V_e的转置，k_i,j为积分单元刚度矩阵的元素值，k_i,j满足以下公式：

式中，(x,y,z)为空间坐标，φ_i、φ_j均为插值函数，当积分网格在自然单元区域时，其为自然单元不同节点对应的插值函数；当积分网格在无限单元区域时，其为无限单元不同节点对应的插值函数。

可选地，步骤S3中，三维自然单元插值形函数满足以下公式：

式中，I为自然邻点编号，φ_I(x)为自然邻点I对应于高斯积分点的插值函数，α_I(x)、α_J(x)为表征插值函数的中间参数，J为自然邻点循环，n为自然邻点数量。

可选地，α_I(x)满足以下公式：

式中，s_I(x)为二阶Voronoi单胞中对应于自然邻点I的平面面积，h_I(x)为高斯积分点到自然邻点I的距离的一半。

可选地，步骤S4中，三维单向映射无限元插值函数满足以下公式：

式中，

为三维单向映射无限元插值函数，i_N＝1，2，…，8，(ξ,η,ζ)为三维单向映射无限元的母单元局部坐标。

可选地，步骤S5包括：

对自然单元区域的背景积分网格进行循环，对自然单元区域的背景积分网格中的高斯积分点进行循环，搜索自然单元区域的背景积分网格中的高斯积分点的自然邻点，得到自然单元区域的背景积分网格中的自然邻点相关于高斯积分点的形函数及其导数，进而得到自然单元区域的背景积分网格中的高斯积分点的子刚度矩阵；

对无限单元区域的无限单元进行循环、计算，得到无限单元区域的子刚度矩阵。

可选地，步骤S6中，将自然单元区域的所有子刚度矩阵和无限单元区域的所有子刚度矩阵按节点编号累加到耦合法的总刚度矩阵中；

自然单元区域和无限单元区域的耦合满足以下公式：

式中，

为耦合法的总刚度矩阵的元素值，

为自然单元区域刚度矩阵的元素值，

为无限单元区域刚度矩阵的元素值，i、j为节点编号。

可选地，步骤S8中，基于总刚度矩阵和源向量，求解大型稀疏方程，得到三维多源自然电位地电模型的自然电位分布，满足以下公式：

KU＝P

式中，K为总刚度矩阵，P为源向量，U为三维多源自然电位地电模型的自然电位分布。

(三)有益效果

本发明的有益效果是：本发明的基于自然单元-无限元的三维多源自然电位数值模拟方法，将基于Laplace自然邻点插值的自然单元区域的边界处耦合无限单元，相对于现有技术而言，其能够有效剖分复杂地电模型，提高数值模拟精度。

附图说明

图1为本发明提供的基于自然单元-无限元的三维多源自然电位数值模拟方法的流程图；

图2为本发明中自然单元-单向映射无限元的三维空间耦合关系示意图；

图3为本发明中三维多源自然电位地电模型的复杂起伏地形示意图；

图4为本发明中六面体中高斯积分点及其自然邻点的示意图；

图5为本发明中高斯积分点与其自然邻点所构成的Delaunay四面体的示意图；

图6为本发明中高斯积分点与其自然邻点所构成的Delaunay四面体的外接球及球心的示意图；

图7为本发明中高斯积分点的二阶Voronoi单胞的示意图；

图8为本发明中三维单向映射无限元的子单元的示意图；

图9为本发明中三维单向映射无限元的母单元的示意图；

图10为本发明中三维多源自然电位地电模型的地表自然电位分布示意图。

具体实施方式

为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。

本发明实施例提出的基于自然单元-无限元的三维多源自然电位数值模拟方法，在自然单元区域边界处耦合无限单元，能够有效解决截断边界问题，有效剖分复杂几何异常体，提高数值模拟精度。

为了更好的理解上述技术方案，下面将参照附图更详细地描述本发明的示例性实施例。虽然附图中显示了本发明的示例性实施例，然而应当理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了能够更清楚、透彻地理解本发明，并且能够将本发明的范围完整的传达给本领域的技术人员。

实施例1

如图1所示，为本发明提供的基于自然单元-无限元的三维多源自然电位数值模拟方法的流程图，该方法包括以下步骤：

S1、构建三维多源自然电位地电模型，三维多源自然电位地电模型包括自然单元区域和无限单元区域。具体地，依据模型尺度、异常体几何参数、起伏地形等因素，在自然单元区域布设离散自然节点，并在自然单元区域除地表外的5个边界处扩展单向映射的无限单元区域充当边界单元以处理截断边界问题。

如图2所示，为自然单元-单向映射无限元的三维空间耦合关系示意图，其中A区域表示自然单元区域，B区域表示无限单元区域，C为单向放射状网格，表示单向映射无限元，作为三维多源自然电位地电模型的边界单元，五角星标记点表示自然单元区域和无限单元区域的耦合连接点。自然单元区域作为三维多源自然电位地电模型的核心研究区域，1个面为地表，其余5个面与无限单元区域进行耦合。

在自然单元区域内分布有自然节点，在无限单元区域分布有无限单元，在自然单元区域与无限单元区域的耦合处设置有连接点，连接点既是自然节点也是无限单元的节点，是自然单元区域和无限单元区域的耦合点。

如图3所示，为三维多源自然电位地电模型的复杂起伏地形示意图，自然单元区域的尺度为50m×50m×50m，在三维多源自然电位地电模型中心距地表10m深度处，设置有五个邻近离散分布且幅值为1A的场源点。为简化三维多源自然电位地电模型，电阻率设为200Ω·m，即均匀模型。自然节点布设132651个，无限单元数量为12500个。对于无限单元在映射方向上的尺度，取水平地表中心为映射起点，设定映射起点到无限单元边界的距离为映射起点到自然单元区域边界距离的3.5倍。对于靠近自然单元区域边界的自然节点，布设时应保证当该区域的高斯积分点在计算刚度矩阵时，所参与计算的自然邻点不包括除连接点之外的无限单元节点。有效的处理办法是，靠近自然单元区域边界的自然节点按例如立方体顶点等简易方式进行布设，以实现和无限单元区域的有效耦合。

上述中，自然节点为所有自然单元区域中呈离散分布的点，高斯积分点为自然单元区域或无限单元区域的背景积分网格中的离散积分点，自然邻点为自然单元区域中的背景积分网格中的当前高斯积分点的自然邻点。

S2、构建三维多源自然电位地电模型的微分方程，基于微分方程得到三维多源自然电位地电模型的积分方程。

三维多源自然电位地电模型的微分方程满足以下公式：

式中，σ为电导率，V为待求自然电位，j_s为场源。

将公式(1)转换为三维多源自然电位地电模型的积分方程，满足以下公式：

式中，Ω为积分单元区域，e为积分单元，V_e为积分单元中各自然邻点的电位值或各无限单元节点的电位值，

为V_e的转置，k_i,j为积分单元刚度矩阵的元素值，k_i,j满足以下公式：

式中，(x,y,z)为空间坐标，φ_i、φ_j均为插值函数，当积分网格在自然单元区域时，其为自然单元不同节点对应的插值函数；当积分网格在无限单元区域时，其为无限单元不同节点对应的插值函数。

S3、基于Laplace插值，构建自然单元区域的插值函数及其导数的基本方程。

为便于理解，以六面体为例，如图4所示，为高斯积分点x(x,y,z)及其8个自然邻点的示意图，高斯积分点的自然邻点为该六面体的8个顶点。如图5所示，为高斯积分点与其自然邻点所构成的Delaunay四面体的示意图。如图6所示，为高斯积分点与其自然邻点所构成的Delaunay四面体的外接球及球心的示意图。如图7所示，为高斯积分点的二阶Voronoi单胞的示意图。

以自然邻点I为例，其对应于高斯积分点的插值函数满足以下公式：

式中，φ_I(x)为自然邻点I对应于高斯积分点的插值函数，α_I(x)、α_J(x)为表征插值函数的中间参数，J为自然邻点循环，n为自然邻点数量。

α_I(x)满足以下公式：

式中，s_I(x)为二阶Voronoi单胞中对应于自然邻点I的平面面积，h_I(x)为高斯积分点到自然邻点I的距离的一半。

插值函数φ_I(x)的导数满足以下公式：

式中，

为φ_I(x)的导数，

为α_I(x)的导数，j_xyz分别取x、y和z。

上述公式(4)、(6)都可转化为坐标运算，进一步地可得到自然单元区域的各积分网格所对应的刚度矩阵(如公式(3)所示)。

S4、基于单向映射的无限单元，构建无限单元区域的插值函数及其导数的基本方程。

如图8-9所示，分别为三维单向映射无限元的子单元和母单元的示意图。在无限单元中，任意点的空间坐标满足以下公式：

式中，(ξ,η,ζ)为母单元局部坐标，

为子单元的空间坐标。

为映射函数，其中i_M＝1，2，…，8，映射函数满足以下公式：

无限单元中任意点的电位值满足以下公式：

式中，

为插值函数，i_N＝1，2，…，8，

满足以下公式：

综合上述公式(7)-(10)都可转化为坐标运算，可得到无限单元区域的各积分网格所对应的刚度矩阵(如公式(3)所示)。

S5、基于步骤S2-S4，对三维多源自然电位地电模型的全域进行积分，得到各积分单元的刚度矩阵。具体地，对自然单元区域的背景积分网格进行循环，对自然单元区域的背景积分网格中的高斯积分点进行循环，搜索自然单元区域的背景积分网格中的高斯积分点的自然邻点，得到自然单元区域的背景积分网格中的自然邻点相关于高斯积分点的形函数及其导数，进而得到自然单元区域的背景积分网格中的高斯积分点的子刚度矩阵；对无限单元区域的无限单元进行循环、计算，得到无限单元区域的子刚度矩阵。

S6、将自然单元区域的所有子刚度矩阵和无限单元区域的所有子刚度矩阵按节点编号累加到耦合法的总刚度矩阵K中，从而实现自然单元区域和无限单元区域的有效耦合。

自然单元区域和无限单元区域的耦合满足以下公式：

式中，

为耦合法的总刚度矩阵K的元素值，

为自然单元区域刚度矩阵的元素值，

为无限单元区域刚度矩阵的元素值，i、j为节点编号。

S7、处理电流源信息，将场源按分布区域、幅值、正负性赋值于相应的空间节点处，构建源向量P。

假设场源所在空间位置对应的节点编号为i_k、i_l、i_m、i_n，场源幅值为1A，则源向量有

S8、基于总刚度矩阵K和源向量P，求解大型稀疏方程组，得到三维多源自然电位地电模型的自然电位分布U，满足以下公式：

KU＝P (12)

如图10所示，为三维多源自然电位地电模型的地表自然电位分布示意图。

综合上述，本发明提供的基于自然单元-无限元的三维多源自然电位数值模拟方法，实现基于Laplace自然邻点插值的自然单元区域与单向映射无限单元区域的有效耦合。其中，自然单元区域中自然节点布置灵活，不受剖分网格单元的约束，能够适用于复杂几何目标体及起伏地形模型。同时，对于自然单元区域的计算最终都可转化为简单的代数运算，编程实现简便且便于理解。在自然单元区域边界处耦合无限单元，能够有效解决截断边界问题，提高数值模拟精度。

本发明提供的基于自然单元-无限元的三维多源自然电位数值模拟方法，能推动自然电位等多源、动态源复杂地电模型数值模拟工作的开展，为相应地球物理方法应用到污染监测、检测以及其他工程与环境地球物理问题中提供数值计算基础，同时也可为反演计算提供正演算法基础。

本发明提供的基于自然单元-无限元的三维多源自然电位数值模拟方法，采用离散点源近似模拟呈区域性分布的自然电位场源，在相应空间节点位置赋点源幅值；同时由于无限单元的引入，不需要考虑源分布对边界条件的影响，适用于多源、动态源模型。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例，或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。

应当注意的是，在权利要求中，不应将位于括号之间的任何附图标记理解成对权利要求的限制。词语“包含”不排除存在未列在权利要求中的部件或步骤。位于部件之前的词语“一”或“一个”不排除存在多个这样的部件。本发明可以借助于包括有若干不同部件的硬件以及借助于适当编程的计算机来实现。在列举了若干装置的权利要求中，这些装置中的若干个可以是通过同一个硬件来具体体现。词语第一、第二、第三等的使用，仅是为了表述方便，而不表示任何顺序。可将这些词语理解为部件名称的一部分。

此外，需要说明的是，在本说明书的描述中，术语“一个实施例”、“一些实施例”、“实施例”、“示例”、“具体示例”或“一些示例”等的描述，是指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或多个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

尽管已描述了本发明的优选实施例，但本领域的技术人员在得知了基本创造性概念后，则可对这些实施例作出另外的变更和修改。所以，权利要求应该解释为包括优选实施例以及落入本发明范围的所有变更和修改。

显然，本领域的技术人员可以对本发明进行各种修改和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也应该包含这些修改和变型在内。

Claims (10)

1.一种基于自然单元-无限元的三维多源自然电位数值模拟方法，其特征在于，包括以下步骤：

S1、构建三维多源自然电位地电模型，所述三维多源自然电位地电模型由多个积分单元组成；

所述三维多源自然电位地电模型包括自然单元区域和无限单元区域；

所述无限单元区域包括多个单向映射的无限单元；

S2、构建三维多源自然电位地电模型的微分方程，基于所述微分方程得到三维多源自然电位地电模型的积分方程；

S3、基于Laplace插值，构建自然单元区域的插值函数及其导数的基本方程；

S4、基于单向映射的无限单元，构建无限单元区域的插值函数及其导数的基本方程；

S5、基于步骤S2-S4，对三维多源自然电位地电模型的全域进行积分，得到自然单元区域和无限单元区域中各积分单元的刚度矩阵；

S6、将自然单元区域和无限单元区域中各积分单元的刚度矩阵按节点编号累加，得到总刚度矩阵；

S7、将场源按分布区域、幅值、正负性赋值于相应的空间节点处，构建源向量；

S8、基于总刚度矩阵和源向量，求解大型稀疏方程组，得到三维多源自然电位地电模型的自然电位分布。

2.如权利要求1所述的三维多源自然电位数值模拟方法，其特征在于，步骤S1中，依据模型尺度、异常体几何参数和起伏地形，在自然单元区域布设离散自然节点。

3.如权利要求2所述的三维多源自然电位数值模拟方法，其特征在于，自然单元区域内分布有自然节点，无限单元区域分布有无限单元，自然单元区域与无限单元区域的耦合处设置有连接点，连接点既是自然节点也是无限单元的节点。

4.如权利要求3所述的三维多源自然电位数值模拟方法，其特征在于，步骤S2中，三维多源自然电位地电模型的微分方程满足以下公式：

式中，σ为电导率，V为待求自然电位，j_s为场源；

三维多源自然电位地电模型的积分方程满足以下公式：

式中，Ω为积分单元区域，e为积分单元，V_e为积分单元中各自然邻点的电位值或各无限单元节点的电位值，

为V_e的转置，k_i,j为积分单元刚度矩阵的元素值，k_i,j满足以下公式：

式中，(x,y,z)为空间坐标，φ_i、φ_j均为插值函数，当积分网格在自然单元区域时，其为自然单元不同节点对应的插值函数；当积分网格在无限单元区域时，其为无限单元不同节点对应的插值函数。

5.如权利要求4所述的三维多源自然电位数值模拟方法，其特征在于，步骤S3中，三维自然单元插值形函数满足以下公式：

式中，I为自然邻点编号，φ_I(x)为自然邻点I对应于高斯积分点的插值函数，α_I(x)、α_J(x)为表征插值函数的中间参数，J为自然邻点循环，n为自然邻点数量。

6.如权利要求5所述的三维多源自然电位数值模拟方法，其特征在于，α_I(x)满足以下公式：

式中，s_I(x)为二阶Voronoi单胞中对应于自然邻点I的平面面积，h_I(x)为高斯积分点到自然邻点I的距离的一半。

7.如权利要求6所述的三维多源自然电位数值模拟方法，其特征在于，步骤S4中，三维单向映射无限元插值函数满足以下公式：

式中，

为三维单向映射无限元插值函数，i_N＝1，2，…，8，(ξ,η,ζ)为三维单向映射无限元的母单元局部坐标。

8.如权利要求7所述的三维多源自然电位数值模拟方法，其特征在于，步骤S5包括：

对自然单元区域的背景积分网格进行循环，对自然单元区域的背景积分网格中的高斯积分点进行循环，搜索自然单元区域的背景积分网格中的高斯积分点的自然邻点，得到自然单元区域的背景积分网格中的自然邻点相关于高斯积分点的形函数及其导数，进而得到自然单元区域的背景积分网格中的高斯积分点的子刚度矩阵；

对无限单元区域的无限单元进行循环、计算，得到无限单元区域的子刚度矩阵。

9.如权利要求8所述的三维多源自然电位数值模拟方法，其特征在于，步骤S6中，将自然单元区域的所有子刚度矩阵和无限单元区域的所有子刚度矩阵按节点编号累加到耦合法的总刚度矩阵中；

自然单元区域和无限单元区域的耦合满足以下公式：

式中，

为耦合法的总刚度矩阵的元素值，

为自然单元区域刚度矩阵的元素值，

为无限单元区域刚度矩阵的元素值，i、j为节点编号。

10.如权利要求9所述的三维多源自然电位数值模拟方法，其特征在于，步骤S8中，基于总刚度矩阵和源向量，求解大型稀疏方程组，得到三维多源自然电位地电模型的自然电位分布，满足以下公式：

KU＝P

式中，K为总刚度矩阵，P为源向量，U为三维多源自然电位地电模型的自然电位分布。

Images