CN108108579B - 直流电阻率无单元法中耦合有限单元法的边界处理方法 - Google Patents

Abstract

本发明提供了一种直流电阻率无单元法模拟中耦合有限单元法的边界处理方法，包括以下步骤：对二维地电模型建立范围较小的无单元法区域Ω₁，该区域采用无单元法计算，使用规则分布的矩形或者平行四边形背景网格覆盖无单元法区域Ω₁，获得该区域无单元法方程组；在无单元法区域Ω₁外围采用快速扩展的有限单元法网格进行剖分，建立足够大的有限单元法区域Ω₂，使得满足第一类边界条件，并采用有限单元法计算获得该区域有限单元法方程组；将两个区域的方程组合并后进行求解，获得观测点的视电阻率参数。本发明能够基于任意节点分布离散模型，对任意复杂的地电模型适应性强，利用有限单元法处理边界大幅提高了传统直流电阻率无单元法正演的计算效率。

Description

直流电阻率无单元法中耦合有限单元法的边界处理方法

技术领域

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

背景技术

直流电阻率勘探是地球物理勘探中的一种重要方法，被广泛应用于固体矿产资源勘探、水文地质勘察、环境治理与监测、工程地球物理勘查等领域。测量的视电阻率与地下介质的电阻率有着直接的关系，通过人工向地下供电，在地表或者井中观测视电阻率可以对地下电阻率异常体分布进行判断。随着直流电阻率勘探技术的发展，对复杂地形、地下介质复杂形态和分布的地电模型的高精度、高适应性和灵活性的正演方法的需求日益增长，无单元法是新兴的一种数值模拟方法(Belytschko,et al.,1994；Hadinia and Jafari,2015)，其仅需节点信息，不依赖网格链接信息，摆脱了网格的约束而具有高灵活性和适应性的特点，同时由于采用高精度的插值方法其也具有高精度的特点，被广泛研究，目前在直流电阻率正演模拟中已获得了应用(麻昌英等,2017)。然而，在传统的直流电阻率无单元法中，由于无单元法计算效率不高，在边界处理时往往采用第三类边界条件，以尽量缩小计算域范围，减少计算成本，但在多电极情况下，无单元法处理第三类边界较费时。如果采用第一类边界条件，以消除边界计算，需要足够大的计算区域，对任意的节分布，无单元法往往需要节点之间的合理布局，计算区域的扩大将严重增加无单元法的计算成本。

因此，有必要设计一种快速扩展计算域、高效的直流电阻率无单元正演的边界处理方法。

发明内容

本发明所解决的技术问题是，针对现有技术的不足，提供了一种直流电阻率无单元法模拟中耦合有限单元法的边界处理方法，能够基于向外快速扩展的有限单元网格，使用少量的节点和网格极大地扩展计算域范围，使得满足第一类边界条件，大幅提高无单元法的计算效率。

本发明的技术方案为：

步骤1、地电模型建立：

首先，确定二维地电模型中介质电阻率的分布情况，以及异常体的几何形态和地形起伏，并设置好电极布设位置、观测装置、观测点位置、供电电流；

建立范围较小但包含地电模型核心区域的无单元法区域Ω₁，在无单元法区域中将二维地电模型采用一组任意分布的节点进行离散，根据地质异常体位置和几何形态、地形起伏形态、电极位置布置节点分布，节点可不规则任意分布，局部可以任意加密节点；

步骤2、在无单元法区域进行无单元法计算：

首先，使用规则分布的矩形或者平行四边形背景网格覆盖无单元法区域Ω₁，平行四边形背景网格可模拟起伏地形，在每一个背景网格中采用无单元法计算边值问题(1)式对应的2.5维直流电阻率变分问题(2)式(徐世浙,1994)；

其中，Ω为计算域，σ＝(x,z)为介质电导率，U(λ,x,z)为波数域电位，λ为波数，I₀为电流，δ为Kroneckerdelta函数，x＝[x,z]^T为Ω上的任意一点，A为场源点位置，▽为梯度运算符，n为外法线单位向量，δ₀为变分符号；F(U)表示函数U的泛函；

然后，在每一个背景网格中布置n_g个高斯积分点x_g，对每一个高斯积分点x_g构造一个局部支持域，使用其内部包含的n个节点信息构造RPIM形函数，利用该组形函数对高斯积分点x_g处的场值进行插值计算，获得该区域无单元法方程组；

在2.5维直流电阻率无单元法正演模拟中，为无单元法区域Ω₁中任意点x^T＝[x,z]构造局部支持域Ω_q，假设Ω_q中包含有n个节点{x₁,x₂,…,x_n}，它们对应的波数域电位值为{U₁,U₂,…,U_n}，则任意点的电位值波数域U(x)的RPIM近似可表示为

式中，p_j(x)为空间坐标x中的单项式，m为基函数个数，二维情况下选m＝3时p^T(x)＝[1x z]。R^T(x)＝[R(x,x₁)R(x,x₂)…R(x,x_n)]为径向基函数向量，R(x,x_i)为径向基函数，本发明的方法中采用复合二次函数(4)式作为径向基函数。

R(x,x_i)＝(r_i ²+(α_cd_c)²)^q (4)

包含α_c和q两个形状参数，d_c为计算点x支持域中节点的平均间距，r_i为计算点x与节点x_i之间的距离，a_i和b_j为待定系数。最终获得矩阵形式的RPIM方程(麻昌英等，2017)

式中

为(n+m)维向量，G为RPIM系数矩阵

其中径向基函数对应的系数矩阵为

多项式基函数对应的系数矩阵为

RPIM形函数为

其中支持域中节点对应的RPIM形函数为

Φ^T(x)＝[φ₁(x) φ₂(x) … φ_n(x)] (10)

为避免对G求逆，(9)式两边同时乘以G

变换可得到

利用标准线性方程组求解器求解式(12)可得到RPIM形函数。利用(12)式可计算RPIM

形函数的一阶和二阶导数:

更高阶导数同理可得。

如附图1所示，假设计算域Ω被分割为无单元法区域Ω₁和有限单元法区域Ω₂，Ω＝Ω₁∪Ω₂，Ω₁中包含N₁个无单元法节点，Ω₂中包含N₂个有限单元法节点，其中无单元法区域Ω₁边界上的N₃个节点为无单元法和有限单元法共用节点，则计算域Ω中节点总数为N＝N₁+N₂-N₃，Ω₁为模型的核心区域，采用无单元法进行正演计算，Ω₂为外围扩边区域，采用有限单元法进行正演计算，采用快速延伸网格扩大计算域，以满足第一类边界条件，认为边界处场值为0，无需计算边界积分，特别指出的是，无单元法区域被有限单元法区域包围，与截断边界Γ_∞无交集，因此无单元法区域内无需做边界计算。

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

K^EFGU^EFG＝F (15)

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

K_qU_q＝F_q (16)

U_q为局部支持域Ω_q内节点的波数域电位U构成的向量，其中

Φ为局部支持域Ω_q内节点的形函数构成的形函数矩阵，K_q为局部支持域子刚度系数矩阵，它的各个元素的表达式为

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

将所有背景积分单元内的积分点的局部支持域子方程组(16)式组装起来可得到直流电阻率无单元法方程(15)式。

步骤3、有限单元法处理边界

在传统的直流电阻率无单元法中，由于无单元法计算效率不高，在边界处理时往往采用第三类边界条件，以尽量缩小计算域范围，减少计算成本，但在多电极情况下，无单元法处理第三类边界较费时。因此，步骤2中的无单元法区域Ω₁通常适当覆盖地电模型核心区域，范围较小。

为了克服无单元法边界处理的不足，利用有限单元法数值稳定性和计算效率高的特点，单元的快速延伸不会造成计算失败，在无单元法区域Ω₁外围采用快速扩展的有限单元法网格进行剖分，建立足够大的有限单元法区域Ω₂，使得满足第一类边界条件，此时认为边界处场值为0，避免第三类边界积分计算，并采用有限单元法计算获得该区域有限单元法方程组；

如附图1所示，在有限单元法区域Ω₂使用矩形单元剖分，假设矩形单元内电导率σ为常数，矩形单元内任意一点近似场值U^h(x)线性变换，对U^h(x)采用双线性插值，根据文献徐世浙(1994)，在矩形单元子域Ω_e上直流电阻率变分问题(14)式可写为子方程组形式

K_eU_e＝S_e (21)

式中K_e＝K_e1+K_e2为有限单元法单元子刚度系数矩阵，它们的表达式分别为

其中

a和b为矩形单元的长和宽。矩形单元四个节点坐标分别表示为x_l＝(x_l,z_l)，l＝1,2,3,4，下标1、2、3和4为矩形单元顶点对应的节点局部编号；U_e＝(U₁,U₂,U₃,U₄)^T为矩形单元节点上的波数域电位列向量。将有限单元法区域Ω₂所有矩形单元的子方程组(19)式组装得到直流电阻率有限单元法方程组

K^FEU^FE＝S (25)

式中K^FE为由全部矩形单元子域Ω_e的K_e组装得到的N₂×N₂维的有限单元法刚度系数矩阵，U^FE为有限单元法区域内所有矩形单元节点上的波数域电位子向量U_e组合成的N₂×1维列向量，S为N₂×1维的有限单元法方程右端项列向量，K_e为矩形单元子域Ω_e形成的子系数矩阵。由于有限单元法形函数满足Kronecker delta函数性质，其各个元素也可表达为(20)式。

步骤4、耦合无单元法与有限单元法：

在无单元Galerkin-有限单元耦合法中，耦合界面Γ_int上两种方法需满足电位相容性条件(26)式：

I表示耦合界面上两种方法的共同节点，

为耦合界面上共同节点的无单元Galerkin法计算的波数域电位，

为耦合界面上共同节点的有限单元法计算的波数域电位，若耦合界面为介质分界面，即两侧介质电阻率不同，此时耦合界面上还需满足电位自然边界条件(27)式：

传统的无单元法Galerkin法通常采用MLS形函数，由于MLS形函数不具有Kronecker delta函数性质，因此不具有相容性，若将各自得到代数方程直接耦合，则会造成两个子域在耦合界面上场值不连续，产生较大误差，两种方法耦合起来比较了困难，本发明的方法采用具有Kronecker delta函数性质的RPIM形函数，使得(26)式自然得到满足，因此可以克服MLS形函数不具有Kronecker delta函数性质而产生的耦合困难，另外，本发明的方法中选择耦合面两侧介质电阻率相同的区域分割有限单元法区域和无单元法区域，使得计算域内部介质界面上满足自然边界条件(27)式，同时将无单元法中局部支持域限定在无单元法区域，将无单元法区域计算获得的无单元法方程组K^EFGU^EFG＝F和有限单元法区域计算获得的有限单元法方程组K^FEU^FE＝S直接合并组装起来，获得整体计算域Ω(其中Ω＝Ω₁+Ω₂)对应的整体方程组KU＝P，其中K＝K^EFG+K^FE，U＝U^EFG+U^FE，P＝F+S，求解整体方程组获得计算域Ω内节点的电位场值，再根据观测装置信息可计算得到观测点的视电阻率参数。

有益效果：

常规的直流电阻率无单元法计算效率不高，在边界处理时往往采用第三类边界条件，以尽量缩小计算域范围，减少计算成本，但无单元法处理第三类边界较费时。本发明将计算域中模型核心部分分为无单元法区域，使用无单元法计算，以更好地模拟复杂地电模型，在外围区域采用快速扩展的有限单元法网格扩展计算域，以满足第一类边界条件，此时认为边界积分为0，避免边界积分计算，仅增加少量的节点和有限单元网格，大幅提高了直流电阻率无单元法的计算效率，利用了RPIM形函数具有Kronecker delta函数性质，使得两方法区域交界面上独立计算的场值具有相容性，使得两种方法在交界面上可直接耦合，处理简便。

本发明可以对二维任意几何形态、地形起伏、电阻率参数分布复杂的地电模型进行直流电阻率法正演计算，大幅提高了直流电阻率无单元法正演的计算效率。

附图说明：

图1为复杂地形起伏模型无单元法区域示意图，展示了模型地形起伏情况和节点布置情况。

图2为复杂地形起伏模型无单元法区域外围区域有限单元法网格剖分示意图，采用了快速扩展的有限单元法网格扩展计算域，以满足第一类边界条件。

图3为常规的直流电阻率无单元法和本发明的方法对复杂地形起伏模型的视电阻率观测正演模拟结果，(a)为常规的直流电阻率无单元法的正演模拟结果，(b)为本发明的方法的正演模拟结果。

具体实施方式：

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

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

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

步骤2、外围有限单元法剖分文件：在模型外围区域建立有限单元法网格剖分文件，确定有限单元法区域范围，网格分布，节点坐标。

步骤3、进行无单元法计算：在无单元法区域，根据模型设计进行直流电阻率无单元法正演计算。

步骤4、进行有限单元法计算：在外围有限单元法区域，根据网格剖分进行有限单元法正演计算。

步骤5、将两种方法组合求解，获得正演模型直流电阻率观测装置视电阻率。

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

如附图1所示，建立为水平方向宽度140m(X:-70～70m)，垂直方向宽度40m(Z:0～-40m)的复杂地形起伏模型的无单元法区域Ω₁，在X＝-14m处为山脊最高处(最高为7.6m)，X＝10处为山谷最低处(最高为3.4m)。在直流电阻率无单元法正演模拟中，可在局部域任意加密节点提高模拟精度，如在场源附近加密节点可显著提高近场源模拟精度。因此，如附图3所示，针对该复杂地形起伏模型，在电场场值变化剧烈的场源区域(即浅地表区域)和地形变化区域采用密集的节点分布以提高模拟精度，同时利用无单元法的强适应性，采用任意的节点分布以适应任意起伏的地形，若有进一步加密节点需求，可直接加密节点而无需其他处理，相比于有限单元法中加入新的节点需要重新剖分单元网格，这体现了无单元法的便利性和灵活性。随着与场源及地形距离的增加，电场场值变化趋于平缓，节点分布随之设置为逐渐稀疏，以降低计算成本。无单元法区域Ω₁共使用7749个节点。如附图2所示，外围采用与快速扩展的的有限单元网格扩展建立有限单元法区域Ω₂，水平范围超过4千米宽，垂直范围超过2千米高，共740个节点，657个单元。复杂地形起伏模型介质电阻率为100Ω·m。在地表X:-58～58m范围内等间隔2m布置59根供电和观测电极，对模型进行高密度电法温纳装置视电阻率观测正演模拟。首先，在无单元法区域先进行常规的无单元法正演，即在模型区域内采用第三类边界条件进行直流电阻率无单元法正演，正演视电阻率如附图3(a)所示。然后将无单元法区域和有限单元法区域结合，进行本发明的利用有限单元法扩边处理边界的直流电阻率无单元法正演，正演视电阻率如附图3(b)所示。两种方法正演所花费的计算时间如表1所示。

表1复杂起伏地形模型视电阻率观测不同模拟方法正演的计算时间

将附图3(a)和(b)对比分析可知，两种方法的模拟结果几乎一致，表明当采用相同的无单元法模型时，本发明的方法可获得与常规的直流电阻率无单元法一致的模拟精度。表1列出了两种方法对复杂地形起伏模型的正演计算时间，常规的直流电阻率无单元法施加第三类界条件而增加了边界计算时间，由于边界计算量与电极数、波数个数、边界段数、高斯点数等成正比，因此在多电极供电观测条件下往往耗费较多计算时间。本发明的方法中由于有限单元法网格数很少，且有限单元法计算效率高，使得有限单元法计算时间很少，但本发明的方法由于增加了有限单元法区域的节点需要花费更多方程求解时间，然而本发明的方法中仅使用少量的有限单元法节点，因此不会明显增加方程求解时间。根据表1，虽然本发明的方法需要花费更多的方程求解时间，但常规的直流电阻率无单元法需要较多的边界计算时间，使得总体上本发明的方法比常规的直流电阻率无单元法花费计算时间大幅减少。表明本发明的方法相比于常规的直流电阻率无单元法，通过使用快速延伸的有限单元法网格扩展计算域，仅少量增加节点的情况下使得满足第一类边界条件，节省了边界计算时间，从而提高了计算效率。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细说明，本领域的普通技术人员应当理解；其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。

Claims (3)

1.一种直流电阻率无单元法模拟中耦合有限单元法的边界处理方法，其特征在于，包括以下步骤：

步骤1、地电模型建立：

确定二维地电模型中介质电阻率的分布情况，以及异常体的几何形态和地形起伏，并设置好电极布设位置、观测装置、观测点位置、供电电流；

建立包含地电模型核心区域的无单元法区域Ω₁，在无单元法区域中将二维地电模型采用一组任意分布的节点进行离散，节点根据地质异常体位置和几何形态、地形起伏形态、电极布设位置布置；

步骤2、在无单元法区域进行无单元法计算：

首先，使用规则分布的平行四边形背景网格覆盖无单元法区域Ω₁，平行四边形背景网格模拟起伏地形，在每一个背景网格中采用无单元法计算如下式所示的2.5维直流电阻率边值问题对应的变分问题；

其中，Ω为计算域，σ＝(x,z)为介质电导率，U(λ,x,z)为波数域电位，λ为波数，I₀为电流，δ为Kronecker delta函数，x＝[x,z]^T为Ω上的任意一点，A为场源点位置，

为梯度运算符，δ₀为变分符号；F(U)表示函数U的泛函；

然后，在每一个背景网格中布置n_g个高斯积分点x_g，对每一个高斯积分点x_g构造一个局部支持域，使用其内部包含的n个节点信息构造RPIM形函数，利用该组形函数对高斯积分点x_g处的场值进行插值计算，获得该区域无单元法方程组K^EFGU^EFG＝F，其中K^EFG为N₁×N₁维的无单元法刚度系数矩阵，U^EFG为无单元法区域节点波数域电位对应的N₁×1维的列向量，F为N₁×1维的无单元法方程右端项列向量；

步骤3、有限单元法处理边界

在无单元法区域Ω₁外围采用快速扩展的有限单元法网格进行剖分，建立有限单元法区域Ω₂，使得计算域边界处满足认为边界处场值为0的第一类边界条件，采用有限单元法计算获得该区域有限单元法方程组K^FEU^FE＝S，式中K^FE为由全部矩形单元子域Ω_e的K_e组装得到的N₂×N₂维的有限单元法刚度系数矩阵，U^FE为有限单元法区域内所有矩形单元节点上的波数域电位子向量U_e组合成的N₂×1维列向量，S为N₂×1维的有限单元法方程右端项列向量；K_e为矩形单元子域Ω_e形成的子系数矩阵；

步骤4、将无单元法区域计算获得的无单元法方程组K^EFGU^EFG＝F和有限单元法区域计算获得的有限单 元法方程组K^FEU^FE＝S合并组装起来，获得整体计算域Ω对应的整体方程组KU＝P，其中Ω＝Ω₁+Ω₂，其中K＝K^EFG+K^FE，U＝U^EFG+U^FE，P＝F+S，求解整体方程组获得计算域Ω内节点的电位场值，再根据观测装置信息，计算得到观测点的视电阻率参数；

具体为：

在无单元Galerkin-有限单元耦合法中，耦合界面Γ_int上需满足电位相容性条件：

I表示耦合界面上两种方法的共同节点，

为耦合界面上共同节点的无单元Galerkin法计算的波数域电位，

为耦合界面上共同节点的有限单元法计算的波数域电位，若耦合界面为介质分界面，即两侧介质电阻率不同，此时耦合界面上还需满足电位自然边界条件：

其中，n为外法线单位向量。

2.根据权利要求1所述的直流电阻率无单元法模拟中耦合有限单元法的边界处理方法，其特征在于：步骤3中所述的有限单元法扩边的边界处理方法：

假设计算域Ω被分割为无单元法区域Ω₁和有限单元法区域Ω₂，Ω＝Ω₁∪Ω₂，Ω₁中包含N₁个无单元法节点，Ω₂中包含N₂个有限单元法节点，其中无单元法区域Ω₁边界上的N₃个节点为无单元法和有限单元法共用节点，则计算域Ω中节点总数为N＝N₁+N₂-N₃，Ω₁为模型的核心区域，采用无单元法进行正演计算，Ω₂为外围扩边区域，采用有限单元法进行正演计算，采用快速延伸网格扩大计算域，以满足第一类边界条件，认为边界处场值为0，无需计算边界积分。

3.根据权利要求2所述的直流电阻率无单元法模拟中耦合有限单元法的边界处理方法，其特征在于：步骤2中背景积分单元Ω_E内直流电阻率无单元法子方程组中元素的表达：

在背景积分单元Ω_E内，对于每一个积分点在其局部支持域Ω_q内无单元法子方程组可表示为

K_qU_q＝F_q (1)

U_q为局部支持域Ω_q内节点的波数域电位U构成的向量，其中

Φ为局部支持域Ω_q内节点的形函数构成的形函数矩阵，K_q为局部支持域子刚度系数矩阵，它的各个元素的表达式为

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

Images