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

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

背景技术：

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

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

技术实现要素：

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

本发明的技术方案为：

步骤1、地电模型建立：

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

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

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

给定初始支持域尺寸dx和dz，一级局部支持域节点数nmax1，一级局部支持域节点数nmax2，高斯点支持域节点数ng，计算域x方向中点的x坐标xcent；其中，初始支持域尺寸dx和dz一般相对较大，可设计为背景单元尺寸的3至10倍，nmax1一般设置为nmax2的2至5倍，nmax2一般设置为ng的2至5倍；

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

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

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

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

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

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

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

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

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

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

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

其中，γ为边界符号，γs为地表边界，γt为截断边界，σ＝(x,z)为介质电导率，u(λ,x)为波数域电位，λ为波数，i0为电流，δ为kroneckerdelta函数，x＝[x,z]^t为ω上的任意一点，a为场源点位置，为梯度运算符，ra为点源与边界上任意一点的直线距离，n为外法线单位向量，cos(r,n)为ra与外法向n的夹角余弦，k0、k1分别为第二类一阶、零阶修正贝塞尔函数，δ为变分符号；

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

ku＝f(3)

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

kquq＝fq(4)

其中

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

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

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

有益效果：

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

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

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

附图说明：

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

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

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

具体实施方式：

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

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

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

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

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

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

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

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

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

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

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。