DC勘探有限元法的两点改进

1、修正的点源二维DC勘探有限元模拟韩思旭1欧东新1 李勇2桂林理工大学资源与环境工程系桂林541004；中国地质科学院地球物理地球化学勘查研究所，廊坊065000摘 要 本文是在文献1中所提出的5个波数来近似模拟的点源二维DC勘探总电位法有限元法的基础上提 出的，介绍了一种简单的对计算结果进行修正的方法，以此来消除总电位法中电源点的影响，并且提高了在有 限个波数的情况下大极距情况下值的可靠性，从而进一步的使得在有限个波数近似计算的结果更加准确。关键词 点源二维 DC三角单元修正系数中图分类号：P631.325文献标识码：A在文献1中提出5个波数来近似模拟的点源二维DC勘探的有限单元法后，相继也出

2、现了很多的文章对其 进行研究，由于上面所提的这5个波数的只能满足100米以内的计算精度，而在超出100米以外以及1米以内 的值误差会急剧的变化，使得此方法大大的受到了限制。于是在此提出对上面所提的有限元计算结果进行修正 的方法，使得电源点附近和100米之外的数据都能很好的与理论数据拟合。在不增加波数情况下提高有限单元 法近似计算的准确性。1、波数域中点源二维电场变分问题F (U) = j (V U )2 + -k 2。U 2 - 16 (A)U d QQK (kr ) cos( r , n) 一 K (kr ) cos( r , n)+ 一 J k 1 AA1 BBU 2 d rK (kr )

3、 - K (kr )r”0 A 0 B TOC o 1-5 h z HYPERLINK l bookmark14 o Current Document 6 F (U ) = 0文献1中所取的5个波数及反富士变换权值为表1：12345K0.0047580.04070110.14088550.3932251.088380g0.00994720.3816190.09803270.25115310.7260814 TOC o 1-5 h z 2、有限单元法在三角单元内电导率和电位都进行线性插值，设电位U和电导率a是单元内的线性函数：U = N u + N u + N u(2)11223=N + N +

4、 N (3)其中N1、N2、N3称为形函数：基金项目：国家自然科学基金项目(40774057)、广西地质工程中心重点实验室开放基金项目(桂科能07109011-K003)、广西 高校人才小高地创新团队建设项目资助。作者简介：1N =(a x+ b y + c )-2 A1111N 二=(a x+b y+c )22A2221N=(a x+b y+c )32A333（4）图1三角单元的节点编号图Kleaibiaaaabijmjjbbbabijm(c + c + c )12 Am m(7)a1=yj-ymb=1xm-xjc = x.y-xmyja2=ym-yib=2xi-xmc = x y-xyma

5、3=y.-yjb=3xj-xic3 = xi y j-xjyiA1=一(2a b -a b )（5）积分的第一部分为j -c(VU )2 dQ = -UTK. U(6)e其中：（8）r cCC 111213K =k 2CCC京 60212223CCCL 313233（9）积分的第二部分为:2 k 2C UT j NNTd Q U其中：C=6b+ 2b+ 2bC= 2b+ 2b+b11ijm12ijmC=2b+ b+ 2bC= 2b+ 2b+b13ijm21ijm(10)C= 2b+ 6b+ 2bC=b+ 2b+ 2b22i-jm23ijmC= 2b+b+ 2bC=b+ 2b+ 2b31ijm3

6、2ijmC= 2b+ 2b+ 6b33ijm积分的第三部分为:j 5 (A )d。= - IUQ(11)由5函数的积分性质，只有电源点的Ua对该项积分才有贡献。积分的第四部分(边界积分)为: K (kr ) cos( r , n) K (kr ) cos( r , n)jb k i a a i a a U2 d r2K (kr ) K (kr )80 A0 B-= uTK u(12)设单元的边落在左无穷远边界上，由于无穷远边界离电源较远，可以将8b k K (kr ) cos( r , n) - K (kr ) cos( r ,n)看作是常数，提到积分号之外，经过适当计算，得: =1 aa1b

7、K (kr ) K (kr )3b + b b + b 0K = B b +b b + 3b 00 0 0(13)K (kr ) cos( r , n) K (kr ) cos( r , n)B = Lk 1 a A1 B :jB-12 L K (kr ) K (kr ) J对所有的三角单元求和，即得到以下的泛函F (U )的数值表达式：F (U ) = - UtKU + Ut P(15)其中K为所有三角单元以及边界单元的(K1e+K2e+K3e)相加组成的N*N阶对称带宽矩阵，而U是由所有三角 单元中的N个节点上的转换电位组成的列矢向量。(15)式泛函F (U )对U求变分，并令其为零，就得

8、到了求波数域中各节点上的电位的线性方程组：K u P=0或K U= P(16)解这个方程组就能得到波数域中的各个节点上的电位值U，计算出各个波数所对应的电之后，再进行反富 士变换，相加就得到空间域的电位值。u (r) = U (r, k ) gj=1(17)3、视电阻率的求取由文献3中视电阻的计算方法可以很方便的求的地面点视电阻率。计算方法为:(18)P s = K 其中：K =竺 （19）+AM AN BM BM4、视电阻率的修正由于本文有限单元模拟中只取了表1中的5个波数1,且这五个波数只能在100米以内满足计算精度。从 计算出的视电阻率结果中可以看出曲线的首端和尾端与理论值都不能拟合（本

9、文可见算例），这是由于首端离电 源点比较近，受到电源点的影响比较大，而尾端则受到所取波数以及计算机舍入误差等的影响。鉴于总体考虑， 在此提出对视电阻率进行修正的办法，以使此方法得以改进。下面介绍一下具体方法和理论：我们可以把大地看成是由一个均匀的背景和局部的异常叠加。简单表示如下：大地=背景* AP（20）我们所做的有限单元计算方法受到距离和计算机舍入误差的影响比较大，且这两者的影响完全是由于计算 方法所引起的，与模型无关，所以在计算时可以对此方法进行适当的修正。由于均匀大地是所有模型所共有的 （P背景）特性，最具有一般性，所以可以以均匀水平大地的模型来计算这些修系数。具体计算如下：QX =

10、P，理/（21）S其中：QX 一为修正系数P理论一为均匀大地模型理论计算视电阻率值p；一为均匀大地模型有限单元法计算视电阻率值5、算例5.1、中间层过渡变化模型的计算p s（。m）设计模型如图2所示，顶层电阻率p 1=50Q m，底层电阻率p 2=1000Q m，中间层电阻率按电导率线性变化 规律从p 1变到p 2,厚度为底层的5倍。图3就是利用解析法与有限元法所得测深曲线对比图。b 1=1/50(0 1nl 深=1/501/1000(0 m) h(m)2=1/1000 (0 .m)图2、中间层过渡变化模型示意图(a)、修正前的对比曲线1010001000.11101001000(b)、修正后

11、的对比曲线1一曲线为有限元模拟曲线；2一曲线为有限元模拟修正曲线；3解析计算理论曲线图3、中间层过渡变化模型的测深曲线与理论曲线对比图由图中计算对比的效果看，在做修正之后结果与理论值拟合得更加好了，而且在100之后数据逐渐趋向于 1000。- m (即底层电阻率)，而不是一直向上升，超出1000。 m。5.2、D型层状模型的计算10001000.1110100100010000(a)、修正前的对比曲线10001000.1110100100010000(b)、修正后的对比曲线1曲线为有限元模拟曲线；2曲线为有限元模拟修正曲线；3解析计算理论曲线；4模型曲线图4、D型水平大地的测深曲线与理论曲线对

12、比图5.3、H型的层状模型的计算0.1110100100010000(b)、修正后的对比曲线1一曲线为有限元模拟曲线；2一曲线为有限元模拟修正曲线;3解析计算理论曲线；4模型曲线图5、H型水平大地的测深曲线与理论曲线对比图5.4、K型的层状模型的计算0.11101001000100000.1110100100010000(b)、修正后的对比曲线1曲线为有限元模拟曲线；2曲线为有限元模拟修正曲线；3解析计算理论曲线；4模型曲线图6、K型水平大地的测深曲线与理论曲线对比图5.4、Q型的层状模型的计算0.1110100100010000(a)、修正前的对比曲线1000100100.111010010

13、0010000(b)、修正后的对比曲线1一曲线为有限元模拟曲线；2一曲线为有限元模拟修正曲线;3解析计算理论曲线；4模型曲线图7、Q型水平大地的测深曲线与理论曲线对比图0.1110100100010000(a)、修正前的对比曲线0.11101001000100005.5、HK型的层状模型的计算(b)、修正后的对比曲线1曲线为有限元模拟曲线；2曲线为有限元模拟修正曲线；3解析计算理论曲线；4模型曲线图8、HK型水平大地的测深曲线与理论曲线对比图6、结论通过以上算例可以看出对点源二维有限元总电位法计算的视电阻率进行系数修改时可行的。尤其是电源点附近的点修改之后与 理论值拟合的非常的好，当然在曲线的

14、尾部的几个点由于受到边界的影响也不能达到完全保证拟合，但是总体上也提高了点源二 维电法正演中5个波数的计算时的有效性。参考文献1徐世浙，地球物理中的有限单元法M，科学出版社，1994。2阮百尧，三角单元部分电导率分块连续变化点源二维电场有限元数值模拟J，广西科学，2001，8-13阮百尧、贺玉田、程志平，电法勘探教程M，桂林工学院，2003, 4144。Mallick,K.and Janin,S.C,1979,Resistivity sounding on a layered transitional earth,GeophysicalProspecting,Vol.27,PP.869875.

15、Modified two-dimensional DC Point-source finite element simulationHan Si-xuiOu Dong-xinLi Yong2Guilin University of Technology Resources and Environmental Engineering, Guilin 541004;Chinese Academy of Geological Sciences Institute of Geophysical and Geochemical Exploration, Langfang 065000Abstract：

16、This article introduced one simple method to amend the resault of calculation which based on the wave numbers of 5 to approximate the point source of two-dimensional simulation of DC total exploration potential finite element method which put forward in the literature 1.This method eliminated the impact of power piont in the total potential method and raised the reliability of values in the case of limited wave numbers and large polar distance.Then furtherly make sure that the re