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

修正的点源二维DC勘探有限元模拟 韩思旭 欧东新1 李 勇2 1.桂林理工大学 资源与环境工程系 桂林 541004； 2.中国地质科学院地球物理地球化学勘查研究所，廊坊 065000 摘 要 本文是在文献[1]中所提出的5个波数来近似模拟的点源二维DC勘探总电位法有限元法的基础上提出的，介绍了一种简单的对计算结果进行修正的方法，以此来消除总电位法中电源点的影响，并且提高了在有限个波数的情况下大极距情况下值的可靠性，从而进一步的使得在有限个波数近似计算的结果更加准确。 关键词 点源二维 DC 三角单元 修正系数 中图分类号：P631.325 文献标识码：A 在文献[1]中提出5个波数来近似模拟的点源二维DC勘探的有限单元法后，相继也出现了很多的文章对其进行研究，由于上面所提的这5个波数的只能满足100米以内的计算精度，而在超出100米以外以及1米以内的值误差会急剧的变化，使得此方法大大的受到了限制。于是在此提出对上面所提的有限元计算结果进行修正的方法，使得电源点附近和100米之外的数据都能很好的与理论数据拟合。在不增加波数情况下提高有限单元法近似计算的准确性。 1、波数域中点源二维电场变分问题 （1） 文献[1]中所取的5个波数及反富士变换权值为表1： 1 2 3 4 5 K 0.004758 0.0407011 0.1408855 0.393225 1.088380 g 0.0099472 0.381619 0.0980327 0.2511531 0.7260814 2、有限单元法 在三角单元内电导率和电位都进行线性插值，设电位U和电导率σ是单元内的线性函数： （2） （3） 其中N1、N2、N3称为形函数： （4） 图1 三角单元的节点编号图 （5） 积分的第一部分为 （6） 其中： （7） 积分的第二部分为： （8） 其中： （9） （10） 积分的第三部分为： （11） 由函数的积分性质，只有电源点的对该项积分才有贡献。 积分的第四部分（边界积分）为： （12） 设单元的边落在左无穷远边界上，由于无穷远边界离电源较远，可以将 看作是常数，提到积分号之外，经过适当计算，得： （13） （14） 对所有的三角单元求和，即得到以下的泛函的数值表达式： （15） 其中K为所有三角单元以及边界单元的（K1e+K2e+K3e）相加组成的N*N阶对称带宽矩阵，而U是由所有三角单元中的N个节点上的转换电位组成的列矢向量。（15）式泛函对U求变分，并令其为零，就得到了求波数域中各节点上的电位的线性方程组： （16） 解这个方程组就能得到波数域中的各个节点上的电位值U，计算出各个波数所对应的电之后，再进行反富士变换，相加就得到空间域的电位值。 (17) 3、视电阻率的求取 由文献[3]中视电阻的计算方法可以很方便的求的地面点视电阻率。计算方法为： （18） 其中： （19） 4、视电阻率的修正 由于本文有限单元模拟中只取了表1中的5个波数[1]，且这五个波数只能在100米以内满足计算精度。从计算出的视电阻率结果中可以看出曲线的首端和尾端与理论值都不能拟合（本文可见算例），这是由于首端离电源点比较近，受到电源点的影响比较大，而尾端则受到所取波数以及计算机舍入误差等的影响。鉴于总体考虑，在此提出对视电阻率进行修正的办法，以使此方法得以改进。下面介绍一下具体方法和理论： 我们可以把大地看成是由一个均匀的背景和局部的异常叠加。简单表示如下： （20） 我们所做的有限单元计算方法受到距离和计算机舍入误差的影响比较大，且这两者的影响完全是由于计算方法所引起的，与模型无关，所以在计算时可以对此方法进行适当的修正。由于均匀大地是所有模型所共有的（）特性，最具有一般性，所以可以以均匀水平大地的模型来计算这些修系数。具体计算如下： （21） 其中： —为修正系数 —为均匀大地模型理论计算视电阻率值 —为均匀大地模型有限单元法计算视电阻率值 5、算例 5.1、中间层过渡变化模型的计算ρs(Ωm) 设计模型如图2所示,顶层电阻率ρ1=50Ωm，底层电阻率ρ2=1000Ωm，中间层电