三维电场直流电阻率边界单元法地形改正

段长生，王绪本，吕玉增，李长伟，黄修保

（1．成都理工大学，成都 610059；2．桂林理工大学，桂林 541004；3．赣中南地质矿产勘查研究院，南昌 330029）

三维电场直流电阻率边界单元法地形改正

段长生1，3，王绪本1，吕玉增2，李长伟2，黄修保3

（1．成都理工大学，成都 610059；2．桂林理工大学，桂林 541004；3．赣中南地质矿产勘查研究院，南昌 330029）

电法观测易受到地形的影响，给资料解释造成困难。这里对数值计算中地面网格节点的立体角（ω）公式进行了严格地推导，利用边界单元法对三维直流电场的地形影响进行了计算并编写了程序。90°角域山脊地形算例表明，解析解与数值解的电位及视电阻率相对误差小于5%，验证了算法及程序的正确性。实例计算了3D山谷地形2D观测剖面的地形效应，当旁侧剖面与主剖面距离大于24 m时，地形效应最大值小于5%。

边界单元法；地形改正；立体角；球面三角形

0 引言

20世纪50年代以来，电法勘探在矿产资源勘查和工程勘察等方面的应用越来越广泛，但地形起伏较大的地区，不仅采集难度加大，资料畸变也严重，给解释造成困难，因此地形影响引起了有关学者们的重视［1－4］。起伏地形不但可以引起假异常，还会掩盖地下由矿体或目标物引起的真异常，如果不能正确地认识和消除，就会导致完全错误的解释结果。

对起伏地形资料的解释，人们一般采用经验判断、带地形直接反演和比值法消除地形的影响。比值法是用实测视电阻率除以“地形校正因子（地形模型视电阻率／地形模型电阻率）”，从而得到改正后的平地形下的视电阻率。比值法中对地形校正因子的计算的方法很多，早期的有基于物理实验的定性解释法；基于数学公式的解析计算方法（比如保角变换法［5］）。但由于定性解释的局限性以及数学公式的复杂性，这些方法的应用受到一定限制。随着计算机技术的发展和地球物理解释方法的进步，众多学者研究了适合于起伏地形的2D、3D数值模拟和解释技术，从而使地形影响问题的模拟、识别和校正成为可能［6］。主要的数值计算方法有：积分方程法、边界元方法、有限差分法和有限元方法等。其中的边界单元法利用格林函数将三维边值问题转换成二维计算，计算量减少、计算速度提高，有利于大型数值计算以及快速反演的实现。

作者从区域边界点立体角的定义出发，依据计算球面三角形面积的原理，利用地面网格点的三维坐标给出了该地形改正公式中的ω（立体角）的计算方法［7］。由于该法仅利用已有网格节点坐标，并不局限剖分单元的形式，因而可以用于多种剖分方式下节点的立体角计算。据此作者编写了地形改正程序，给出了90°角域山脊地形算例，计算结果表明，解析解与数值解的电位及视电阻率相对误差小于5%，验证了算法及程序的正确性。实例计算了3D山谷地形2D剖面的地形效应，在设计模型下，当测线与主剖面的距离小于10 m时，视电阻率受地形影响仍然较大，中心值大于1．1Ω·m，而当距离大于24 m时，中心值小于1．05Ω·m，地形效应最大值小于5%，几乎可以忽略。

1 边界元模拟法

1．1 边值问题

假定地下介质是均匀的，电阻率为1Ω·m，地表A点电流强度为1 A，则电位满足的基本方程为［8］

其中：ωA是电源点A对地下区域所张的立体角；δ是狄拉克函数。

1．2 边界积分方程

其中：ωp是计算点p的对地下区域所张的立体角；up、u0p分别是计算点p的电位、正常电位；是与计算单元的法向向量的夹角。

1．3 边界元法

1．3．1 单元积分

式（2）中第二项利用边界单元法求解，首先用三角形单元将计算区域划分为m个小区域，单元内采用线性插值则

其中：l是三角元顶点；Nl是形函数；xl、yl、zl是顶点坐标，

1．4 立体角ω的计算方法

根据式（3）求解方程组时还需要确定各网格节点的立体角。对于不光滑边界上的孤立点p，立体角ωp定义为：以p点为球心，作一半径为ε（无穷小）的球体，若区域内球体表面积s，则ωp＝s／ε2［1］。

在水平面（图1中红线）上，x、y将整个区域按照矩形剖分后再剖分为三角形，x方向剖分nx＿n次。其中一内部节点p的剖分情况如图1所示，第n、n－nx＿n＋1、n＋1点的坐标分别为（xn，yn，zn）、（xn－nx＿n＋1，yn－nx＿n＋1，zn－nx＿n＋1）、（xn＋1，yn＋1，zn＋1）。对第n点做一半径为ε（ε趋于无穷小）的球体，此时该点的立体角由地面以下区域内球体的表面积表示。按图1的方法将地面以下表面积的计算分为八个部分（现以图1中第二部分的计算为例），则在水平面上，n－nx＿n＋1、n＋1点到n点及n－nx＿n＋1点到n＋1点的距离分别为a1＝

图1 内部节点网格示意图Fig．1 Schematic diagram for inner node mesh

待求表面积为图2蓝线区域部分，可由θ1、θ2、α、ε来表示。

图2 球面三角形面积计算示意图Fig．2 Schematic diagram for Spherical triangle area calculation

同样，可以求出其他七个部分的表面积，叠加就得到整个的表面积。从而易得p点的立体角为

2 地形校正

利用式（3）可求得网格节点电位

其中：K为装置系数；ΔU为电位差；I为供电电流；ρT五数值计算视电阻率。

当I＝1A，大地电阻率ρ＝1Ω·m，计算得到的ρT在数值上等于校正因子αT，若实测视电阻率为ρs，则经比值法地改后的视电阻率

3 算例

3．1 算例1

计算模型为验证模型，即有解析解的90°角域山脊地形，如图3所示，模型的电阻率为1Ω·m。由于y坐标具有对称性，1 A正负点源供电点的x、z坐标分别为：Axz＝（－40 m、－40 m），Bxz＝（40 m、－40 m）。在供电点中间范围x＝－20．5 m到x＝20．5 m测量。数值解与解析解电位对比如图4所示，数值解与解析解曲线形态一致，数值相近，表1将数值解与解析解电位和视电阻率值进行了统计，21对电位和视电阻率值数值和解析计算结果中，相对误差均较小，其中电位最大相对误差为3．87%，平均相对误差为1．73%；电阻率最大相对误差为2．34%，平均相对误差为0．64%。

图3 90°角域山脊地形Fig．3 90°ridge topography model

图4 数值解与解析解电位对比图Fig．4 Comparison between analytical and numerical solutions

3．2 算例2

模型为地形效应特征分析模型，电阻率为1 Ω·m的山谷地形，中心坐标为（0 m，0 m，－4．631 m）。x、y坐标范围均为－40 m～40 m。沿y方向除主剖面外布置了5条线，距离主剖面依次为5 m、10 m、24 m、30 m、45 m，边界单元法计算的视电阻率异常曲线图示于图5（a），图5（b）是相应的地形示意。由图5可见，总体来讲，距离主剖面越远地形效应越小，当测线与主剖面的距离小于10 m时，视电阻率受地形影响仍然较大，中心值大于1．1Ω· m，而当距离大于24 m时，中心值小于1．05Ω·m，地形效应最大值小于5%，几乎可以忽略。

3．3 算例3

模型为地改效果模型，凹地形下方存在一个低阻异常球体，球体半径4 m，埋深12 m，异常球体电阻率为0．05Ω·m，区域内电阻率为1Ω·m，如图6（a）所示，图中只画出了主剖面的地形，实际地形是绕对称轴旋转面成的凹陷。在主剖面x坐标为－39 m，39 m处分别供正负点电源。首先用有限单元法计算模型主剖面中间梯度法的视电阻率异常，如图6（b）；然后用边界单元法计算由地形引起的视电阻率异常，如图6（c）；最后用比值法进行地形影响改正，并将改正结果与水平大地有限元计算的结果作比较，如图6（d）所示。显然图6（b）中的正异常是由于地形引起的，直接解释将会得到错误的结果。图6（d）表明二者相差较小，经计算的平均误差仅为2．721%。

表1 数值解与解析解电位和视电阻率值统计表Tab．1 Statistics values of the potential and apparent resistivity of numerical solution and the analytical solution

图5 地形及边界单元法视电阻率异常曲线图Fig．5 Topography and abnormal curve of apparent resistivity boundary element method

图6 模型及数值计算对比图Fig．6 Model and the comparison chart of numerical calculation

网格大小为61×27，在Intel（R）Core（TM）2 Duo CPU，2．66 GHz，3．00 GB内存，Windows XP下运行时间不到5 s，数值计算效率较高。

4 结语

作者从点源场的边值问题出发，通过格林公式建立边界积分方程。利用区域剖分网格节点的坐标，采用求解球面三角形面积的方法，推导出了求解边界单元法地形改正公式中立体角的计算方法，将此方法运用到电阻率三维边界单元法数值模拟中，实现了对积分方程的求解。算例说明该方法误差较小，且由于边界元法有类似于降维的效果，可快速处理较大的网格数据，因此可作为数值模拟理论计算或者实际生产中消除地形效应的一种方法。

［1］ 汤洪志，刘庆成，龚育龄，等．边界单元法在高密度电阻率法二维地形改正中的应用效果［J］．物探与化探，2001，25（6）：457－479．

［2］ 陈伯舫，侯作中，范国华，等．有限差分法计算三维地形影响的电磁感应［J］．地震学报，1998，20（5）：541－544．

［3］ 黄兰珍，田宪漠．电阻率法地形改正及其在工程地质勘查总的应用［J］．物探化探计算技术，1997，19（3）：238－241．

［4］ 李正容，何继善，温佩琳．电阻率和CSAMT法地形改正［D］．长沙：中南大学，1991．

［5］ 徐世浙，王柏钧．保角变换在电法勘探中的应用［M］．北京：地质出版社，1977．

［6］ 刘树才，何昭友，刘志新，等．适合地形起伏的二维有限单元数值模拟技术［J］．物化探计算技术，2005，27（2）：131－134．

［7］ 《数学手册》编写组．数学手册［M］．北京：人民教育出版社，1979．

［8］ 徐世浙．地球物理中的边界单元法［M］．北京：科学出版社，1994．

Topographic correction of 3D electric field in DC resistivity by boundary element method

DUAN Chang-sheng1，2，WANG Xu-ben1，2，LV Yu-zeng3，LI Chang-wei3，HUANG Xiu-bao2
（1．Chengdu University of Technology，Chengdu 610059，China；2．Guilin University of Technology，Guilin 541004，China；3．Institute of Geological and Mineral Exploration in Central and South Jiangxi Province，Nanchang 330029，China）

It is difficult to interpret DC electrical prospecting data which were easily Influenced by the terrain．In this paper，the solid angle of the ground grid nodes for the numerical calculation of formula were rigorously deduced，3D topography effect of the DC electric field were calculated by the boundary element method using fortran Language．The results of ridge terrain of 900 angular domain show the relative error of potential and apparent resistivity are less than 5 percent and verify the correctness of the algorithm and program．3D topography effect were analyzed in valley terrain by 2D sections．In the design mode，when the lateral distance from the cross section to the main cross－section greater than 24m，the maximum value of the topographic effect less than 5%．

boundary element method；topographic correction；solid angle；spherical triangle

P 631．3＋22

A

10．3969／j．issn．1001-1749．2014．05．03

1001-1749（2014）05-0529-06

2014-05-12 改回日期：2014-07-25

国家高技术研究发展计划（863计划）项目（2009AA06Z108）；广西自然科学基金项目（2013GXNSFAA019277）

段长生（1974－），男，博士，主要从事电磁法正演模拟和反演成像，E-mail：49119428＠qq．com。