CN108710156B - 一种直流电阻率无单元法模拟的支持域快速构造方法 - Google Patents

Abstract

本发明提供了一种直流电阻率无单元法模拟的支持域快速构造方法，包括以下步骤：确定二维地电模型的介质、地形、电极等分布情况，建立直流电阻率无单元法计算域，利用不规则节点离散地电模型；先根据节点的x坐标再根据z坐标由小到大对所有节点进行排列编号，给定初始支持域参数；以背景单元中心点为中心，根据参数构造一级局部支持域，在该支持域中，将背景单元分为四个子域，以子域中心点为中心，构造二级支持域，在该支持域中，为子域中每一个高斯点选取距离其最近的N_g个节点作为支持域；进行无单元法计算，获得观测视电阻率参数。本发明能够基于任意节点分布离散模型，快速构造支持域，提高了常规直流电阻率无单元法正演的计算效率。

Description

一种直流电阻率无单元法模拟的支持域快速构造方法

技术领域

本发明涉及一种勘探地球物理领域的直流电阻率正演方法，特别涉及复杂地电模型的高精度、高灵活性和高适应性无单元正演方法。

背景技术

直流电阻率勘探是地球物理勘探中的一种重要方法，被广泛应用于固体矿产资源勘探、水文地质勘察、环境治理与监测、工程地球物理勘查等领域。测量的视电阻率与地下介质的电阻率有着直接的关系，通过人工向地下供电，在地表或者井中观测视电阻率可以对地下电阻率异常体分布进行判断。随着直流电阻率勘探技术的发展，对复杂地形、地下介质复杂形态和分布的地电模型的高精度、高适应性和灵活性的正演方法的需求日益增长，无单元法是新兴的一种数值模拟方法(Belytschko,et al.,1994；Hadinia and Jafari,2015)，其仅需节点信息，不依赖网格链接信息，摆脱了网格的约束而具有高灵活性和适应性的特点，同时由于采用高精度的插值方法其具有高精度的特点，被广泛研究，目前在直流电阻率正演模拟中已获得了应用(麻昌英等,2017)。然而，直流电阻率无单元法常规的支持域构造方法中，对每一个高斯积分点均在全局节点内独立进行支持域构造，需对所有节点按照节点编号由小到大逐一搜索，当支持域内节点数不满足设定要求时，需扩大或减少支持域尺寸重新构造支持域，直至满足支持域设定要求，耗费大量的时间，限制了无单元法在直流电阻率正演模拟中的应用。

因此，有必要设计一种快速的支持域构造方法，提高直流电阻率无单元正演计算效率。

发明内容

本发明所解决的技术问题是，针对现有技术的不足，提供了一种直流电阻率无单元法模拟的支持域快速构造方法，给定支持域相关参数，首先分别构造一级和二级局部支持域，在根据高斯积分点位置，在对应的二级局部支持域中，根据节点至高斯积分点距离由近到远选取给定数量节点作为高斯积分点支持域，加快支持域构造速度，大幅提高无单元法的计算效率。

本发明的技术方案为：

步骤1、地电模型建立：

首先，根据二维地电模型中介质电阻率的分布、异常体的几何形态和地形起伏情况建立无单元法区域Ω，并设置好电极布设位置、观测装置、观测点位置，在无单元法区域中将二维地电模型采用一组任意分布的节点进行离散，根据地质异常体位置和几何形态、地形起伏形态、电极位置布置节点分布情况人工局部可以任意加密节点，如地形变化区域，供电电极附近，异常体及其周围区域，并在根据正演模拟需求，在场值变化不大或者远离场源电性不变的区域使用稀疏的节点分布，以使节点合理分布，减少计算成本；

步骤2、对节点排序编号及给定初始参数：

如附图1所示，先根据节点的x坐标再根据z坐标由小到大对所有节点进行排列编号；

给定初始支持域尺寸d_x和d_z，一级局部支持域节点数N_max1，一级局部支持域节点数N_max2，高斯积分点支持域节点数N_g，计算域x方向中点的x坐标X_cent；其中，初始支持域尺寸d_x和d_z一般相对较大，可设计为背景单元尺寸的3至10倍，N_max1一般设置为N_max2的2至5倍，N_max2一般设置为N_g的2至5倍；

步骤3、支持域构造设计：

首先，在计算域内采用矩形背景单元覆盖计算域，在背景单元内布置高斯积分点；

其次，对每一个背景单元进行循环：

计算背景单元中心点坐标，以背景单元中心点为中心，使用初始支持域尺寸，若中心点x坐标小于等于X_cent，则由节点编号由小到大搜索包含在内的节点，当所搜索的节点x坐标与中心点坐标差值大于d_x时退出搜索，若中心点x坐标打大于X_cent，则由节点编号由大到小搜索包含在内的节点，当所搜索的节点x坐标与中心点坐标差值小于-d_x时退出搜索；当所包含的节点数小于N_max1时，则将初始支持域尺寸扩大两倍，重复上述步骤，直至所包含的节点数大于N_max1；

计算该组节点至中心点的距离，根据它们的距离使用倒序冒泡法对节点进行排列，当排列好最近的N_max1个节点时退出节点排列，选取该N_max1个节点作为一级局部支持域；

接着，将背景单元根据其四条边中点划分为四个子域，以子域中心点为中心，在一级局部支持域中，选择距离中心最近的N_max2个节点作为二级局部支持域；

最后，对该背景单元内每一个高斯积分点构造支持域：

在高斯积分点对应的二级局部支持域中，计算N_max2个节点至高斯积分点的距离，根据它们的距离使用倒序冒泡法对节点进行排列，当排列好最近的N_g个节点时退出节点排列，选取该N_g个节点作为该高斯积分点的支持域；

对于边界积分计算，同理计算边界中心点，采用如上步骤构造边界高斯积分点支持域。

步骤3、对地电模型进行无单元法计算：

对每一个高斯积分点，利用其支持域内N_g个节点信息构造RPIM形函数(Liu andGu,2007；麻昌英等,2017)，利用该组形函数对高斯积分点处的场值进行插值计算，使用高斯积分计算2.5维直流电阻率边值问题(1)式对应的变分问题(2)式(徐世浙,1994)；

其中，Γ为边界符号，Γ_S为地表边界，Γ_T为截断边界，σ＝(x,z)为介质电导率，U(λ,x)为波数域电位，λ为波数，I₀为电流，δ为Kronecker delta函数，x＝[x,z]^T为Ω上的任意一点，A为场源点位置，为梯度运算符，r_A为点源与边界上任意一点的直线距离，n为外法线单位向量，cos(r,n)为r_A与外法向n的夹角余弦，K₀、K₁分别为第二类一阶、零阶修正贝塞尔函数，δ为变分符号；

采用Galerkin全局弱式法可导出变分问题(2)式的直流电阻率无单元法方程：

KU＝F (3)

其中K为N×N维的无单元法刚度系数矩阵，N为计算域中节点总数，U为无单元法区域节点波数域电位对应的N×1维的列向量，F为N×1维的无单元法方程右端项列向量，在背景网格Ω_e内，对于每一个积分点在其局部支持域Ω_q内(3)式可写为子方程组形式

K_qU_q＝F_q (4)

其中

为局部支持域子刚度系数矩阵，为体积分项子刚度矩阵，为采用第三类边界条件时的边界积分项子刚度矩阵，它的各个元素的表达式为

其中i,j＝1,2,…,n，n为局部支持域内包含的节点总数，φ_i和φ_j分别为支持域中第i和j个节点的形函数；由于RPIM形函数具有Kronecker delta函数性质，右端项F_q的各个元素表达式为(10)式；

将所有背景网格内积分点的局部支持域子方程组(4)式组装起来可得到直流电阻率无单元法方程(3)式，对无单元法方程(3)式进行求解，获得节点电场场值，根据观测装置计算获得观测点的视电阻率参数。

有益效果：

直流电阻率无单元法常规的支持域构造方法中，对每一个高斯积分点均在全局节点内独立进行支持域构造，需对所有节点按照节点编号由小到大逐一搜索，当支持域内节点数不满足设定要求时，需扩大或减小支持域尺寸重新构造支持域，直至满足支持域设定要求，耗费大量的支持域构造时间，大幅增加了无单元法的计算量，导致计算效率低。

本发明的支持域构造方法，根据背景单元中心点分布情况，选择节点搜索方向，根据节点的排序原则适时退出搜索，使用较大的初始支持域尺寸，尽可能避免二次或多次搜索，以背景单元中心点为中心构造一级局部支持域，将背景单元划分为四个子域，以子域中心点为中心，，在一级局部支持域内构造二级局部支持域，然后以子域内高斯积分点为中心，在二级局部支持域内构造高斯积分点的支持域，避免了为每一个高斯积分点在全局域节点内构造支持域，大幅缩减了支持域构造时间，提高直流电阻率无单元法正演计算效率。

本发明可以大幅缩减支持域构造时间，提高了直流电阻率无单元法正演的计算效率。

附图说明：

图1为本发明中的支持域构造示意图。其中，1、地表边界，2、节点，3、截断边界，4、背景单元，5、背景单元中心点，6、一级局部支持域，7、高斯积分点，8、子域中心点，9，二级局部支持域，10、高斯积分点支持域。

图2为基岩起伏地电模型及节点分布示意图。

图3为基岩起伏模型正演观测视电阻率拟断面图,(a)～(c)为本法发明方法计算结果，(d)～(f)为常规支持域构造方法计算结果；(a)～(c)分别对应参数N_max1＝32/48/64，N_max2＝16/24/32，N_g＝4/8/16；(d)～(f)分别对应参数N_g＝4/8/16。

具体实施方式：

以下结合附图和具体实施方式对本发明作进一步的说明。

本发明涉及的直流电阻率观测计算方法包括以下步骤：

步骤1、正演地电模型参数文件的设计：根据二维地电模型中介质电阻率的分布、异常体的几何形态和地形起伏情况，设置模型离散节点信心文件，并设置电极布设、观测装置和无单元法相关参数。

步骤2、节点排序编号及给定初始参数设计：先根据节点的x坐标再根据z坐标由小到大对所有节点进行排列编号：给定初始支持域尺寸d_x和d_z，一级局部支持域节点数N_max1，一级局部支持域节点数N_max2，高斯积分点支持域节点数N_g。

步骤3、支持域构造设计：以背景单元中心点为中心，根据其位置确定搜索方向，使用较大的初始支持域尺寸和给定参数构造初始支持域，在初始支持域内构造背景单元中所有高斯积分点支持域。

步骤4、进行无单元法计算：在计算域上，根据模型设计进行直流电阻率无单元法正演计算，获得直流电阻率无单元法方程组，对获得的直流电阻率无单元法方程组进行求解，获得正演模型节点电场场值，根据观测装置计算获得观测视电阻率。

以下为本发明计算一个基岩起伏地电模型的高密度温纳装置视电阻率观测的实例。

建立一个宽200m(X:-100～100m)，高100m(Z:0～100m)基岩地电模型，其电阻率和节点分布如附图2所示。在地表X:-58～58m范围内等间隔2m布置59根供电和观测电极，对模型进行高密度电法温纳装置视电阻率观测正演模拟。采用无单元法模拟时，直流电阻率控制方程的计算最终转化为方程(8)和(9)的计算，而计算方程(8)和(9)首先要求对积分点出场值进行插值，其关键为通过为积分点构造局部支持域，利用局部支持域节点信息对积分点场值进行插值，局部支持域以积分点为中心，包含距离积分点最近的N_q个节点，N_q的选取较为灵活，一般根据具体模型和模拟要求选取，但较大的N_q会增加计算量，因此一般N_q取4～16。本发明方法中，给定初始支持域尺寸d_x和d_z分别为背景单元长和高的3倍，一级局部支持域节点数N_max1＝32/48/64，二级局部支持域节点数N_max2＝16/24/32，高斯积分点支持域节点数N_g＝4/8/16，常规支持域构造方法中，积分点支持域内包含节点数也设计为N_g＝4/8/16。

采用不同方法和参数的视电阻率模拟结果如附图3所示，附图3(a)～(c)为采用本发明方法获得的视电阻率拟断面图，附图3(d)～(f)为采用常规支持域构造方法获得的视电阻率拟断面图，对比附图3(a)～(c)和附图3(d)～(f)可知，两种方法获得的计算结果几乎一样，表明本发明方法可以获得与常规支持域构造方法相同的模拟精度。表1给出同一台计算机上两种方法使用的节点数、背景单元数、高斯积分点数和不同参数时的支持域构造时间，对比分析可知，相同的节点数、背景单元数、高斯积分点数和支持域构造参数相同情况下，本发明的方法的支持域构造时间相比于常规支持域构造方法大幅减少，从而提高了直流电阻率无单元法计算效率。

表1两种方法的节点数、背景单元数、高斯积分点数、支持域参数和支持域构造时间

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本

发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种直流电阻率无单元法模拟的支持域快速构造方法，其特征在于，包括以下步骤：

步骤1、地电模型建立：

首先，根据二维地电模型中介质电阻率的分布、异常体的几何形态和地形起伏情况建立无单元法区域Ω，并设置好电极位置、观测装置、观测点位置，根据地质异常体位置和几何形态、地形起伏形态、电极位置布置节点分布情况，在无单元法区域中将二维地电模型采用一组任意分布的节点进行离散；

步骤2、对节点排序编号及给定初始参数：

先根据节点的x坐标再根据z坐标由小到大对所有节点进行排列编号，给定初始支持域尺寸d_x和d_z，一级局部支持域节点数N_max1，二级局部支持域节点数N_max2，高斯积分点支持域节点数N_g，计算域x方向中点的x坐标X_cent；

步骤3、支持域构造设计：

首先，在计算域内采用矩形背景单元覆盖计算域，在背景单元内布置高斯积分点；其次，以背景单元中心点为中心，使用初始支持域尺寸，搜索包含在内的节点，选取距离中心点最近的N_max1个节点作为一级局部支持域；将背景单元划分为四个子域，以子域中心点为中心，在一级局部支持域中选取距离中心最近的N_max2个节点作为二级局部支持域；最后，在二级局部支持域中，为在对应的子域内的每一个高斯积分点选取距离其最近的N_g个节点作为支持域；

步骤4、对地电模型进行无单元法计算：

对每一个高斯积分点，利用其支持域内N_g个节点信息构造RPIM形函数，利用该组形函数对高斯积分点处的场值进行插值计算，使用高斯积分计算2.5维直流电阻率边值问题对应的变分问题；

对所有背景单元均进行无单元法计算后，获得该无单元法方程组KU＝F，其中K为N×N维的无单元法刚度系数矩阵，N为计算域中节点总数，U为无单元法区域节点波数域电位对应的N×1维的列向量，F为N×1维的无单元法方程右端项列向量；对无单元法方程进行求解，获得节点电场场值，根据观测装置计算获得观测点的视电阻率参数。

2.权利要求1所述的一种直流电阻率无单元法模拟的支持域快速构造方法，其特征在于：给定初始支持域尺寸d_x和d_z，一级局部支持域节点数N_max1，二级局部支持域节点数N_max2，高斯积分点支持域节点数N_g，计算域x方向中点的x坐标X_cent；其中，初始支持域尺寸d_x和d_z相对较大，设计为背景单元尺寸的3至10倍，N_max1设置为N_max2的2至5倍，N_max2设置为N_g的2至5倍；

步骤3支持域构造设计：

对每一个背景单元进行循环：

计算背景单元中心点坐标，以背景单元中心点为中心，使用初始支持域尺寸，若中心点x坐标小于等于X_cent，则由节点编号由小到大搜索包含在内的节点，当所搜索的节点x坐标与中心点坐标差值大于d_x时退出搜索，若中心点x坐标大于X_cent，则由节点编号由大到小搜索包含在内的节点，当所搜索的节点x坐标与中心点坐标差值小于-d_x时退出搜索；当所包含的节点数小于N_max1时，则将初始支持域尺寸扩大两倍，重复上述步骤，直至所包含的节点数大于N_max1；

计算该组节点至中心点的距离，根据它们的距离使用倒序冒泡法对节点进行排列，当排列好最近的N_max1个节点时退出节点排列，选取该N_max1个节点作为一级局部支持域；

接着，将背景单元根据其四条边中点划分为四个子域，以子域中心点为中心，在一级局部支持域中，选择距离中心最近的N_max2个节点作为二级局部支持域；

最后，对该背景单元内每一个高斯积分点构造支持域：

在高斯积分点对应的二级局部支持域中，计算N_max2个节点至高斯积分点的距离，根据它们的距离使用倒序冒泡法对节点进行排列，当排列好最近的N_g个节点时退出节点排列，选取该N_g个节点作为该高斯积分点的支持域；

对于边界积分计算，同理计算边界中心点，构造边界高斯积分点支持域。