三维地形直流电阻率有限元法模拟

1、第50卷第5期2007年9月地球物理学报CHINESEJOURNALOFGEOPHYSICSVol.50,No.5Sep.,2007强建科,罗延钟.三维地形直流电阻率有限元法模拟.地球物理学报,2007,50(5):16061613QiangJK,LuoYZ.TheresistivityFEMnumericalmodelingon32Dundulatingtopography.ChineseJ.Geophys.(inChinese),2007,50(5):16061613三维地形直流电阻率有限元法模拟强建科1,2,罗延钟21中南大学信息物理工程学院,长沙4100832中国地质大学地球物理与空间

2、信息学院,武汉430074摘要基于稳定电流场的基本方程、下电阻率的有限单元数值模拟算法.离散积分区域时,性插值型函数以及单元刚度矩阵.,能够节约内存.利用Cholesky,当求解有多个供电点的测深问题,在水平层状介质模型上,三维计算结果与解析解或二维数值解十分吻合,计算精度满足误差要求在二维山脊上的二极剖面或三维山谷上的中间梯度剖面上,其三维计算结果与相应模型的土槽实验结果或边界元法计算结果也非常接近.关键词电阻率法,有限单元法,三维数值模拟,起伏地形文章编号0001-5733(2007)05-1606-08中图分类号P631收稿日期2006-12-20,2007-05-17收修定稿There

3、sistivityFEMnumericalmodelingon32DundulatingtopographyQIANGJian2Ke1,2,LUOYan2Zhong21SchoolofInfo2physicsandGeomaticsEngineering,CentralSouthUniversity,Changsha410083,China2InstituteofGeophysics&Geomatics,ChinaUniversityofGeosciences,Wuhan430074,ChinaAbstractBasedonstableelectricityflowfieldequat

4、ion,three2dimensionboundaryquestionanditscorrespondingvariationquestion,thealgorithmofresistivityFEMnumericalmodelingona32Dtopographyisstudied.Theintegralareaisdividedintomanytri2prismunitswithanomalousshapesandlandformcharacteristics.Thefunctionoftri2linearinterpolationandthematrixoftri2prismunitar

5、ethendeduced.Usingvariedbandmatrixandone2dimensionarraytocollectthelargesparsematrix,theintegralcoefficientofallthenodes,wecouldsavememoryspace.Beforesolvingthesystemoflinearequations,theCholeskyalgorithmwasadoptedtodecomposethelargesparsematrixonlyonce,andthenweobtainedthepotentialsofallnodesofhund

6、redsofdifferentcurrentelectrodesbybackwardsubstitutions.Thisisanimportantsteptoimprovethespeedofforwardmodelingwhentherearemoresounding2points.Theresultsofthismodelingshowthatthe32Dnumericalsolutionisinagreementwiththeanalyticsolutiononlevelstratifiedmedium,on22Dmountainridgewithpoletopolearray,oron

7、32Dmountainvalleywithcentralgradientarray.The32Dnumericalsolutionisalsoinagreementwiththeresultsfromphysicalexperimentorboundaryelementmethod.KeywordsResistivitymethod,Finiteelementmethod(FEM),32Dnumericalmodeling,Undulatingtopography基金项目国家自然科学基金项目(60672042),中国地质大学(北京)地下信息探测技术与仪器教育部重点实验室开放基金(GDL0502

8、),中南大学博士后基金联合资助.作者简介强建科,男,1967年生,中南大学在读博士后,主要从事地球物理电磁法正反演算法研究.E2mail:qiangjianke1631com5期强建科等:三维地形上直流电阻率有限元法模拟1607式做数学变换,再依据半空间函数的积分性质得1引言在进行电法勘探时,起伏地形严重影响直流电阻率法或电磁法的应用,有时甚至使异常面目全非.因此,用数值方法计算地形对视电阻率的影响并设法消除它,是提高资料地质解释效果的关键之一.国内外专家已经作过一些有关起伏地形三维电1阻率正演模拟问题的研究工作,Fox等.阐述了地形对电阻率和极化率的影响问题;TrumanHT和2George

9、RJ研究了三维地形对电阻率的影响并进行校正;YuguoL和KlausSpitzer比较了三维直流电阻率有限单元和有限差分算研;兰47812等、徐世浙等,1314时,.Sasaki、熊彬等采用四面体单元剖分目标区的有限单元法实现了三维起伏地形电阻率的正演模拟计算,缺点是单元数增大15176倍,不利于进一步做反演研究.吴小平等利用有限差分法、有限单元法开展了起伏地形下电阻率激发极化32D正、反演技术研究,但该算法的正演计算建立在以长方体为单元的基础上,对实际地形的模拟存在明显不足,因为地表四个点不能惟一地确定一个面.本文研究以三棱柱为单元的有限单元正演算法,可以实现任意起伏地形的电阻率模拟计算.三

10、棱柱具有以下优点:只要给定地面上任何三点的高程,就可以惟一确定一块小面积,而且对高程的变化没有特别限制,这就意味着三棱柱单元可以模拟任意坡度的地形.求解变分方程时,采用有限单元法把目标区域离散成一系列三棱柱单元.含有地形特征的不规则三棱柱,有五个面六个节点,单元刚度矩阵含有6×6个元素,其中有21个非零元素,每个元素由30多项含有5个变量的三重积分组成,若手工计算非常复杂且容易出错,本文通过编制积分算法程序,只要计算一次就可得出21项非零元素的数学表达式,可节约推导时间.3到三维微分方程:+9x99y99z9=-2I(xA)(yA)(zA).三维电场的边值问题是=0,9n+u=0,9

11、nr(2)S(地表)(地下边界).(3)+F(u)=013有限单元法311网格剖分为较好地模拟不平地形,同时又尽可能使网格简单化和规则化,并考虑以后反演方便,正演模拟中采用了三棱柱体网格单元.三维有限单元模拟的三棱柱网格由中心区网格和扩展区网格两个部分组成.考虑到电法的野外工作方式、计算精度、计算量以及使程序简化,采用图1形式的单元剖分方法:在每一个长方体单元内,进一步剖分出交叉对称的三棱柱单元.其中,每个三棱柱单元中的电导率是均匀的.三棱柱单元个数为M×N×L×2,三棱柱网2稳定电流场电位满足的基本方程为了后面讨论方便,(1)式列出有关稳定电流场点电源电位满足的

12、基本方程,详细讨论见文献18.(1)(u)=-j,其中,为电导率,u为电位,j为电流密度.对(1)图1网格剖分示意图Fig.1Schematicdiagramforfiniteelementmesh(u-2I(A)udrud,()(4)1608地球物理学报(ChineseJ.Geophys.)50卷格节点的个数为(M+1)×(N+1)×(L+1).在三维电阻率法勘查时(如高密度),通常采用正方形测网,即测线距和测点距相同.以水平方向呈正方形分布的地面测点,作为三维有限单元模拟中心区网格的地面节点;由其向下引铅垂线,并取铅垂线上按测量电极距等距分布的点,作为中心区网格的地下节

13、点.当地面水平时,这样形成的中心区网格,为节点均匀分布的正方体网格;当地面起伏不平时,水平方向呈正方形分布的四个相邻地面节点可能不在一个平面上,因而不能组成正方体网格.为避免出现上述情况,连接上述四个相邻地面测点的两个对角点(约定连接左前方到右后方的对角点),角形,体单元.,上底和下底为相互平行的、随地面起伏变化的倾斜直角三角形.为减小网格边界的影响,从中心区网格分别向左、右、前、后和下方,按公比值113的等比数列取步长,布置1015个节点,形成扩展区网格.312任意三棱柱单元分析按前一节所述,三维有限单元模拟的网格被剖分为一系列直立三棱柱体.图2示出了其中的一个有代表性的三棱柱体单元.单元的

14、三个棱142536OZ,即都沿铅垂方向,平行于Z轴;它们的长313局部坐标按文献19“上、下底三角形可用面积坐标”,而沿高度方向则用距离坐标来构成局部坐标,式(5)的详细推导见文献18.=,x=,yh(5)=x+y+z-h12,zh,x,y,z为实际坐标,xx21,y=3,z4-z1,hij=hi-hj(i,j=,ij),hi是节点i上方的地面高程.314三线性插值型函数)表示的三文献19给出了用局部坐标(,棱柱体单元型函数的表达式,将(5)式代入其中,便可写出以实际坐标(x,y,z)表示的三棱柱体单元型函数的表达式:1(P)=1-xy1-+zxzyzz,h12h13h122(P)=1-+xz

15、xzyzz,3(P)=1-+yzxzyzz,度相同(即为单元沿Z方向的步长),设为z.三个棱面都是平行四边形,其中,平面1254在坐标平面XOZ上,棱12和54在X轴上的投影长度(即单元沿X方向的步长)设为x;而平面1364在坐标平面YOZ上,棱13和64在Y轴上的投影长度(即单元沿Y方向的步长)设为y.上底和下底是与地形起伏4(P)=1-x-y+-zxzyzzh12h13h12,一致的、相互平行的倾斜三角形.5(P)=+-xzxzyzz,h12h13h126(P)=+-yzxzyzz.基于上面公式给出的型函数,便可在每个三棱柱体单元中,按公式(6),由三棱柱体单元六个节点的电位U1,U2,U

16、3,U4,U5和U6,用线性插值计算给定的坐标点(x,y,z)的电位U(x,y,z),进而实现对电位泛函极值问题的离散化.6U(x,y,z)=图2不平地形三棱柱单元Fig.2Tri2prismcellonundulatingtopographyi=1(x,y,z)U.ii(6)315单元积分为了计算变分方程(4)式中的积分,首先计算各5期强建科等:三维地形上直流电阻率有限元法模拟yh-a12yh1-abh1216091单元e上的积分,把该积分分为三部分:Fe(U)、Fe(U)和Fe(U),则第一部分积分为e23=e=TU(k)U2eijeT1(7)UKU.2ee1其中Ke=(kij),kij=

17、kji,是三棱柱体单元e的电导率,Ue=(U1,U2,U6),T为了便于计算积分,如图2点编号,假设a=,=1,=14,其中i,j=1,2,kij:kij=edxdydz+9x9x9y9y9z9zk11k21k22k31k32k33k41k42k43k44k51k52k53k54k55k61k62k63k64k65k依据上面系数我们可以组成单元刚度矩阵:Fe(U)=1=2(U)2d2h12abc+h1-92+92+9dxdydz,dzdydx,+9x9x9y9y9z9z(8)把对x,y,z的一阶偏导数代入(8)式即可求出系数矩阵(9)式.计算(9)式积分系数是一件非常繁琐的运算,整理时一定要细

18、心.由(8),其中右向12ah13bh12ah13,ac12b12ac1222T,前三项为水平地形时的系数变量,当第四和第六项或第五和第七项为零时,7个系数变量就变成5个,这是二维起伏地形时的系数变量.4-402-2040-22=000004-4040040-420-2000004-20240-40042-1-1-211431-4-341-3-42-1-14340330-3-3-6-330303-300-3-36300330-3-30-330-3-63300-3-30362-1-1-211431-4-341-3-42-1-14342-1-1-211431-412a12b12cbh-312a4a

19、h131-3-42-1-143412bbh12ac1222,(9)1610地球物理学报(ChineseJ.Geophys.)50卷k11k21Ke=1k12k22k32k42k52k62k13k23k33k43k53k63k14k24k34k44k54k64k15k25k35k45k55k65k16k26k36k46k56k矩阵相加,得F(u)=F(u)ek31k41k51k61,=13T2UeKe+KeUe-KeuAI22UTU2TeKeUe-KeuAI2第二项积分为Fe(u)=2TTU KeU-UP2KeU-UPT2I(A)Ud=KeUAI,2(10)=1(10)式积分只与电源点有关,I为

20、电流.第三项积分与边界有关,cos(r,n)23Fe(u)=Udr12542TTUKU-P,23(12)=e+Ke,K=Ke,()=2r6i=1iijjA.对(12)式求变分,并令其为零,)(xd,(11)得线性方程组KU=P,T3=UeKeUe.2(13)式中,K为系数矩阵,U为电位向量,P为与点源有关的向量,解方程组得各节点的电位值U.在集成系数矩阵K时,本文采用变带宽、一维数组方式只存储系数矩阵中非零元素,同时记录非若单元e的一个面1254落在无穷远边界上时,其系数矩阵为42Ke1254=340120000042040零元素在矩阵中的位置,这种做法能够节约内存.;360210317算例依

21、据上述理论基础,选择总电位为研究对象,以三棱柱为计算单元,编制了起伏地形三维正演模拟程序.试算结果表明,在CPU主频为219GHz的计算机上,当依次在125个地面节点上供电,计算18630个节点上的电位时所花费的时间只有6min.首先用水平地形一维介质模型进行检验.同理,当单元e的面1463落在无穷远边界时,系数矩阵为40Ke1463=3000004102402002201;36模型1:三层一维介质,K型曲线,第一层电阻率为50m,厚度为5m;第二层电阻率为200m,厚度为5m;第三层电阻率为10m.解析解采用数字滤波的方法计算,测量装置采用对称四极.在18个极距中,最大误差为2195%,该精

22、度已经可以满足要求了,见图3a.模型2:为了检验电阻率差异较大(10量级)情3当单元e的面123落在向下无穷远边界时,系数矩阵为2-1Ke123=343000400000000123况下算法的稳定性和有效性.三层一维介质,K型曲线,第一层电阻率为50m,厚度为5m;第二层电.24-1000阻率为50000m,厚度为5m;第三层电阻率为m.解析解采用与模型1相同的算法,见图3b.10可以看出,在视电阻率曲线出现最大值以前,数值解与解析解很接近,随着极距加大,尽管曲线形态变化一致,但对应极距的视电阻率误差却逐渐变大了,说明本算法对于具有很大电阻率差异的模型能够计算且基本稳定.316单元集成在单元内

23、,将(7),(10)和(11)式中Ke,Ke,Ke扩展成由全体节点组成的矩阵,再把全体单元系数5期强建科等:三维地形上直流电阻率有限元法模拟1611图3解析解与数值解对比图(a)21=4;(b)21=10001Fig.3Comparisonbetweenanalyticalandnumericalsolutions模型3:为了检验本算法对起伏地形的计算效果,我们选择二维山脊地形,采用二极视电阻率剖面m,AO=法,山坡倾角45°,均匀介质,电阻率为11m.图4是二维有限元、三维有限元数值模拟和土成镜像,这与理论分析结果一致.图中实线是徐世浙(1985)用边界单元法计算的中间梯度装置的视

24、电阻率,虚线是本算法计算结果,可以看出两种结果非常相似,产生误差的主要原因是由于有限单元模型不能完全与边界单元模型一致造成的.槽物理模拟视电阻率曲线对比图.可以看出,三种结果基本吻合,证明本算法是有效可行的.图5三维纯地形中间梯度视电阻率曲线对比Fig.5ComparisonbetweennumericalsolutionsbyFEMand图4二、三维有限元数值模拟与土槽模拟视电阻率曲线对比Fig.4ComparisonbetweennumericalsolutionsbyFEMandexperimentalsolutionsformodelridgebyBEMformodelvalley模型

25、4:图5是一个三维山谷地形(图中只绘出剖面,实际地形是绕对称轴旋转而成的凹陷),在图上面是视电阻率曲线,其曲线变化与地形凹凸基本模型5:图6是一个山脊模型上用二极测深法计算的视电阻率异常断面图(只给出中间一条断面).围岩电阻率为200m,山脊正下方有一个6mm.采用二×4m×4m大小的低阻体,其电阻率为1极测深法,记录点为AM中点,最小间隔为1,最大1612地球物理学报(ChineseJ.Geophys.)50卷图6(Fig.6The232Drange间隔为14.从断面图看出,特征,本来在水平地形上规则变化的异常,由于地形的作用,在低阻异常背景上叠加了一个高阻信息,因此,如

26、果类似的测量结果不做地形校正或反演计算,必然得出错误的结论.HuangLZ,TianXM.Topographiccorrectionforresistivitymethodanditsapplicationtoengineeringgeologicalprospecting.ComputingTechniquesforGeophysicalandGeochemicalExploration(inChinese),1997,19(3):2382415黄兰珍,田宪漠.点源场电阻率法二维地形改正边界元法.物探化探计算技术,1986,8(3):201207HuangLZ,TianXM.Thebound

27、aryelementmethodof22Dtopographiccorrectionfortheresistivityinthepointsourcefield.ComputingTechniquesforGeophysicalandGeochemicalExploration4结论本文从稳定电流场的基本方程出发,推导出含有地形特征的三棱柱单元型函数、系数积分式及单元的系数矩阵,通过集成并求解线性方程组,实现了三维起伏地形条件下电阻率法的有限单元数值模拟算法.算例表明,该算法模拟结果与一维解析解十分吻合,与二维山脊土槽实验、二维有限单元模拟结果和三维边界元结果也非常接近,因此可以证明该算法是正

28、确的,而且速度快,精度高,可以满足对起伏地形的模拟要求,也为进一步研究反演算法奠定了基础.参考文献(References)1FoxRC,HohmannGW,KillpackTJ,etal.Topographiceffectsinresistivityandinduced2polarizationsurveys.57932TrumanHT,GeorgeRJ.Three2dimensionalterraincorrectionsinresistivitysurveys.Geophysics,1984,49(4):4394523YuguoL,KlausSpitzer.Three2dimensiona

29、lDCresistivityforwardmodelingusingfiniteelementsincomparisonwithfinite2differencesolutions.GeophysicalJournalInternational,2002,151(3):9249344黄兰珍,田宪漠.电阻率法地形改正及其在工程地质勘查中Geophysics,1980,45:(inChinese),1986,8(3):2012076田宪漠,黄兰珍.点源场电阻率法三维地形改正的边界元法.成都地质学院学报,1986,13(3):170175TianXM,HuangLZ.Theboundaryeleme

30、ntmethodofthe32Dtopographiccorrectionfortheresistivitymethodinthepointsourcefield.ComputingTechniquesforGeophysicalandGeochemicalExploration(inChinese),1986,13(3):1701757黄兰珍,田宪漠,徐荣.模式勘查及其在工程地质调查中的应用.物探化探计算技术,1994,l6(4):294302HuangLZ,TianXM,XuR.Modellingexplorationanditsapplicationtotheengineeringgeo

31、logicinvestigation.1994,16(4):2943028徐世浙,赵生凯.三维地形上点电源电场的边界单元解法.桂ComputingTechniquesforGeophysicalandGeochemicalExploration(inChinese),林冶金地质学院学报,1985,5(2):163168XuSZ,ZhaoSK.Theboundaryelementmethodcalculatingelectricfieldofapointsourceonthree2dimensiontopography.JournalofGuilinCollegeofGeology(inChin

32、ese),1985,5(2):1631689徐世浙,倪逸.复杂地电条件下点源三维电阻率模拟的新方法.物探化探计算技术,1991,11(1):1320XuSZ,NiY.Anewmethodformodelingtheresistivitysurveyswithpointsourceon32Darbitrarynonhomogeneousgeoelectricalsection.ComputingTechniquesforGeophysicalandGeochemicalExploration(inChinese),1991,11(1):132010XuShizhe,GaoZucheng,Zhao

33、Shengkai.Anintegralformulationforthree2dimensionalterrainmodelingforresistivitysurveys.Geophysics,1988,53:546552的应用.物探化探计算技术,1997,19(3):2382415期强建科等:三维地形上直流电阻率有限元法模拟161311徐世浙.用边界单元解法求解起伏地形上的215维电场值.桂geoelectricfieldofpointsourcebyincompleteCholeskyconjugategradientmethod.ChineseJ.Geophys(inChinese),1998,41(6):84885516吴小平.非平坦地形条件下电阻率三维反演.地球物理学报,2005,48(4):932936林冶金地质学院学报,1984,4:119113XuSZ.Theboundaryelementmethodforsolving2152Delectricalfield