地球物理学进展  2014, Vol. 29 Issue (4): 1926-1931   PDF    
引用本文
彭艳华, 柳建新, 刘海飞, 孙丽影. 直流电阻率2.5维有限元数值模拟[J]. 地球物理学进展, 2014, 29(4): 1926-1931, doi: 10.6038/pg20140461  
PENG Yan-hua, LIU Jian-xin, LIU Hai-fei, SUN Li-ying. 2.5-D direct current resistivity modelling by finite-element method. Progress in Geophysics, 2014, 29(4): 1926-1931, doi: 10.6038/pg20140461   
直流电阻率2.5维有限元数值模拟
彭艳华1,2, 柳建新1 , 刘海飞1, 孙丽影3    
1. 中南大学地球科学与信息物理学院, 长沙 410083;
2. 深圳市水务规划设计院, 深圳 518000;
3. 浙江省工程物探勘察院, 杭州 310005
收稿日期: 2013-9-1 修回日期: 2014-2-22 .
基金项目：国家自然科学基金资助项目（41174103）和国家“十二五”科技支撑计划项目（2011BAB04B08）联合资助.
作者简介：彭艳华，男，1990年生，江西樟树人，中南大学地球探测与信息技术硕士研究生.主要从事应用地球物理研究.（E-mail：pengyanhuajx@126.com）
通讯作者：柳建新，男，1962年生，湖南岳阳人，博士，教授，主要从事应用地球物理及应用研究.（E-mail：ljx6666@126.com）
摘要：从稳定电流场满足的基本微分方程出发，推导了地面点源电场电位边值问题对应的等价变分问题.用网格离散积分区域时，以四边形为最小研究单元.针对电场在电源附近衰减快，变化梯度大，数值模拟结果在电源附近位置有较大误差问题，采用边界/奇异校正算法改善模拟结果在电源附近的奇异性问题，同时减少边界条件的影响，提高模拟结果精度.从波数域电位求取空间域电位，利用最优化离散波数进行反付氏变换.模型算例结果验证了该算法的可行性与有效性.
关键词： 数值模拟     最优化离散波数     边界/奇异校正    
2.5-D direct current resistivity modelling by finite-element method
PENG Yan-hua1,2, LIU Jian-xin1 , LIU Hai-fei1, SUN Li-ying3    
1. School of Geoscience and Info-Physics, Central South University, Changsha 410083, China;
2. water planning and design institute of Shenzhen, Shenzhen 518000, China;
3. engineering geophysical investigation institute of zhejiang province, Hangzhou, 310005, China
Abstract: Based on the governing equation of stable electrical current flow field, we obtain the variational equations corresponding to boundary value problems. The integral area is discretized into rectangle elements. With respect to problem of electrical potential falls off rapidly near a source, we apply correction for boundary conditions and the source singularities to improve the accuracy of the simulated results. Optimization wavenumber is used to solve inverse Fourier-cosine transformation. The example model case demonstrate the feasibility and efficiency of the forward routine.
Key words: numerical modeling     optimization wavenumber     correction for boundary condition and sources singularities    

0 引 言

在直流电阻率数值模拟方面国内外许多学者进行了大量的相关研究(罗延钟和孟永良，1986；熊斌等，2003；强建科和罗延钟，2007；徐凯军，2007； Zhou et al., 2009).关于二维地电构造条件下点源场的数值模拟问题，许多研究者利用有限单元法(Mao and Bao， 1997)、有限差分法(罗延钟和张桂青，1987)、边界单元法(徐世浙，1994)、自适应(王飞燕，2009；Tang et al., 2010； Tang et al., 2011)有限元法实现2.5维正演数值模拟.积分方程法的优点易于处理三维局部异常体问题，有限差分法是一种比较老的数值模拟方法，计算原理和编制程序都比较简单，适用于计算平面规则形体(如板状体，层状或似层状体)的异常，有限单元法是50年代中期发展起来的一种数值模拟方法，适用于模拟形状不规则的地质体和地形起伏，但计算程序编制较复杂，尤其对三维地电条件下的模拟更为复杂.在2.5维数值模拟的中涉及到反付氏变换问题，徐世浙计算了一组通用的波数及反付氏变换系数(徐世浙，1994)，(Xu等2000；柳建新和刘海飞，2005)提出了最优化离散波数的优化算法，Cardarelli(2006)利用艾尔米特插值进行反傅里叶变换，提高正演结果精度.在电阻率正演模拟中存在电源奇异性问题，利用异常场电位法或在电源位置附近减小网格剖分间距，增加剖分单元数目改善电源附近位置数值模拟结果奇异性问题.本文利用边界/奇异校正(Pidlisecky et al., 2007；Pidlisecky and Knight， 2008)法改善电源及附近位置数值模拟结果奇异性问题，该校正算法同时减小边界条件对数值模拟结果的影响，提高数值模拟结果精度.

1 有限单元法基本原理 1.1 稳定电流场电位边值问题

点源位于三维无限半空间Ω的表面Γs，电位满足微分方程

其中σ为介质电导率，u为电位，I为电流强度，rs为供电点位置，δ为跟电源供电点出有关的狄拉克拉函数.

电场的边值问题

对方程式(1)进行付氏和反付氏变换可进行方程求解.付氏变换、反付氏变换形式为

其中U(x，k，z)代表波数域电位值，u(x，y，z)代表空间域电位值，k为波数.

当电性分布沿某一方向不变(假设沿y方向不变)，对式(1)进行付氏变换后将上述三维问题变成波数域二维问题：

1.2 网格剖分

采用四边形单元对求解区域进行网格剖分，将(6)式中对区域Ω和边界Γ∞的积分分解为对各单元e和Γe的积分之和，即：

总电位边值问题(5)对应的变分问题为：

1.3 线性插值

矩形单元内电位和电导率均采用双线性插值，即：

式中Ni(i=1，2，3，4)为形函数，可表示为

其中ξi，ηi是点i(i=1，2，3，4)的坐标.根据图 1可知，ξ，η与x，z的对应关系为

图 1 区域Ω网格剖分、子单元和母单元的对应关系 (a)有限单元法剖分示意图；(b)单元. Fig. 1 Regional Ω grid discretion，subunits and parent unit
其中x0，z0是子单元中点的坐标；a，b为单元的两个边长；微分关系为

1.4 单元分析

对式(8)中的各积分项进行单元分析：

其中，K 1e=(k1ij)为4×4阶对称矩阵，k1ij=k1ji，i，j=1，2，3，4.

1.5 总体合成

在单元内，将(13)，(14)，(15)和(16)相加，再扩展成由全体节点的矩阵或列阵，由全部单元相加，得：

对(17)式求变分并令其为零，得线性方程组

解线性方程组，得各节点的 U.

2 最优化离散波数算法

地表有两个点电源，地下为均匀半空间各项同性介质，则偶极子源在地表处产生的电场值可表示为(Pidlisecky and Knight， 2008)：

式(19)中r是观测点位置，r+，r-分别是供电电极位置.对式(19)进行傅里叶变换，得到

式(20)中，U为傅里叶变换后的电位值，K0代表零阶贝塞尔函数.

对式(20)进行反傅里叶变换并进行简化近似计算，得到

选取一系列的r，采用矩阵的方式得到

式(22)中，

为使计算波数在设置的测量电极序列上的精度达到最高，采用最小二乘算法使得

式(23)为波数k与加权系数g的函数，其中 I 为单位矩阵，上式对加权系数g进行求导，得到以k为变量的g的表达式

式(24)代入式(23)得到

上式(25)为波数k的非线性函数，通过迭代法解(25)式可求解最优化波数，将所求波数代入式(24)便可得到加权系数值.由于M是变量r的函数，上式(25)是超定函数，且贝塞尔函数衰减极快，导致式(25)在求解过程中会存在不收敛的情况.基于以上原因，将上式修改为

J 为灵敏度矩阵，W为模型加权矩阵，β为阻尼因子，通常较小，一般为10~9.

将上式对k求导，得到第k次更新修正公式

新的迭代波数序列为

重复上述过程，直到误差满足要求停止.

3 方程组右端项边界/奇异校正算法

数值模拟是用有限空间去求解无限大空间，如何有效处理优化边界问题对模拟结果精度有很大的影响，扩大求解区范围增大网格剖分空间能一定程度上降低边界条件对模拟结果的影响，提高数值模拟结果的精度.网格剖分过程中，为改善模拟结果在电源位置附近奇异性问题，在电源位置附近需减小网格步长(采用较密的网格)，增大单元剖分数目，减小单元剖分步长.本文利用校正算法(Lowry et al.，1989； Pidlisecky et al.，2007；Pidlisecky and Knight， 2008)，对(18)右端含电源项P进行校正，无需增加电源附近单元剖分数目和扩大求解剖分区域，减少边界条件对模拟结果的影响，提高电源及附近位置数值模拟结果精度，由可得

将(29)代入式(4)，令y=0且对积分离散化，得到

联合式(30)与式(29)得

如果已知空间域电位 u值，利用式便可求出右端含源项P .边界校正/源校正算法如下：

(1)求取电导率值为 σ H均匀半空间下电位解析解 u H，其中 σ H等于地下真实模型电导率的平均值.

(2)利用(32)式求出校正项 P corr= L－1u H .

将 P corr代入式(31)，得

利用上式对电场 u 进行边界/奇异校正.

假设一均匀半空间，电阻率值100 Ωm，模型尺寸40 m×40 m，剖分区域不外延，水平位置±20 m为网格剖分的边界，矩形剖分单元尺寸大小1 m×1 m.模拟结果如图 2所示.

从图 2看出，观测点离剖分区域边界距离越小，模拟结果精度降低，利用边界/奇异校正处理后，数值模拟结果精度得到改善.图 2表明通过利用边界/奇异校正可以有效改善电源奇异性问题，同时还可以减小边界条件的影响，提高模拟结果的精度.

图 2 利用边界/奇异校数值模拟结果曲线图 (analytical代表解析解，NoBC/singularity correction代表未利用边界/奇异校正，BC/singularity correction代表利用边界/奇异校正.(a)电位值曲线；(b)相对误差曲线.) Fig. 2 Curves of simulated data obtained by using BC/singularity correction

4 模型算例

4.1 水平层状模型

图 4是三层水平层状模型与对称四级测深装置示意图.A1、B1、A2、B2 …是供电电极，M、N是接收电极.模型参数第一层电阻率100 Ωm，厚度12 m，第二层电阻率500 Ωm，厚度15 m，第三层电阻率10 Ωm.该层状模型的视电阻率有解析解.网格剖分参数：剖分数目个数46×52，网格大小不等，模型尺寸大小为1090 m×3200 m.图 3是双频激电对称四级测深装置模拟的视电阻率数值解与解析解曲线对比图，从表 1、图 3可以看出，与没有经过边界/奇异校正处理后的数值解精度相比，经过边界/奇异校正处理后的数值解具有更高的精度.

图 4 三层水平层状模型与对称四级测深装置示意图 Fig. 4 Sketch map of three-layered model and symmetric quadrupole array

图 3 解析解与数值解对比图 (analytical代表解析解，BC代表利用边界/奇异校正后的数值解，No BC代表未利用边界/奇异校正) Fig. 3 Comparison of analytical and simulated results

表 1 对称四级装置电阻率测深数值模拟结果与解析解

(analytical代表解析解，BC利用边界/奇异校正的数值解，NoBC没有利用边界/奇异校正的数值解，rms_BC利用边界/奇异校正的相对误差，rms_NoBC没有利用边界/奇异校正的相对误差) Table 1 analytical solution and simulated data of symmetric quadrupole array

4.2 均匀半空间低阻体模型

图 5为均匀半空间中一低阻体模型.异常体埋深2 m，规模1.5 m×3 m，数据模拟用温纳装置，电极总数36根，模拟数据共195个.图 6正演结果电阻率断面等值线图.

图 5 均匀半空间中一低阻体模型 Fig. 5 A low resistivity block embedded in the half-space

图 6 正演结果电阻率断面等值线图 Fig. 6 Pseudosection of the apparent resistivity of modelling result

4.3 均匀半空间两个异常体模型

图 7为地电断面图.均匀半空间中同时存在低阻体跟高 阻体.围岩电阻率值100 Ωm，低阻体值10 Ωm，高阻体值500 Ωm，规模1.5 m×3 m，间距2 m，采用温纳装置，电极总数36根，模拟数据共195个.图 8为温纳装置正演结果电阻率断面等值线图.

图 7 地电断面图 Fig. 7 Two blocks embedded in the half-space

图 8 正演结果电阻率断面等值线图 Fig. 8 Pseudosection of the apparent resistivity of modelling result
5 小 结

从稳定电流场满足的基本微分方程出发，推导了点源电流场边值问题对应的变分问题，实现了基于边界/奇异校正的直流电阻率2.5D有限元正演数值模拟.算例结果表明边界/奇异校正算法能改善源附近位置奇异性问题，同时还能减小边界条件对模拟结果的影响，提高模拟结果精度.将波数域电位转换成空间域电位时，采用最优化离散波数进行反付氏变换，提高数值模拟结果的精度.


致 谢 在研究工作期间，得到柳建新教授的亲切指导，在此表示感谢！

参考文献
[1] Cardarelli E, Fischanger F. 2006. 2D data modeling by electrical resistivity tomography for complex subsurface geology[J]. Geophysical Prospecting, 54(2): 121-133, doi: 10. 1111/j. 1365-2478. 2006. 00522. x.
[2] Liu J X, Liu H F. 2005. The optimum algorithm for the calculation of the optimized discrete wave-number[J]. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 27(1): 34-38.
[3] Lowry T, Allen M B, Shive P N. 1989. Singularity removal, a refinement of resistivity modeling techniques[J]. Geophysics, 54(6): 766-774, doi: 10. 1190/1. 1442704.
[4] Luo Y Z, Meng Y L. 1986. Some Problems on resistivity modeling for two dimensional structures by finite element method[J]. Chinese Journal of Geophysics (in Chinese), 29(6): 613-621.
[5] Luo Y Z, Zhang G Q. 1987. Application of electrical Computers in Electrical Exploration (in Chinese)[M]. Wuhan: Wuhan College of Geology Press.
[6] Mao X J, Bao G S. 1997. 2. 5-D resistivity tomography using boundary integral equations[J]. J. Cent. South Univ. Technol., 4(2): 104-107.
[7] Pidlisecky A, Knight R. 2008. FW2_5D: A MATLAB 2. 5-D electrical resistivity modeling code[J]. Computers & Geosciences, 34(12): 1645-1654, doi: 10. 1016/j. cageo. 2008. 04. 001.
[8] Pidlisecky A, Haber E, Knight R. 2007. RESINVM3D: A 3D resistivity inversion package[J]. Geophysics, 72(2): h1-h10, doi: 10. 1190/1. 2402499.
[9] Qiang J K, Luo Y Z. 2007. The resistivity FEM numerical modeling on 3-D undulating topography[J]. Chinese Journal of Geophysics (in Chinese), 50(5): 1606-1613.
[10] Ruan B Y. 2001. 2-D electrical modeling due to a current point by fem with variation of conductivity within each triangular element[J]. Guangxi Science (in Chinese), 8(1): 1-3.
[11] Tang J T, Wang F Y, Xiao X. 2011. 2. 5-D DC resistivity modeling considering flexibility and accuracy[J]. Journal of Earth Science, 22(1): 124-130. doi: 10. 1007/s12583-011-0163-z.
[12] Tang J T, Wang F Y, Ren Z Y. 2010. 2. 5-D DC resistivity modeling by adaptive finite-element method with unstructured triangulation[J]. Chinese J. Geophys., 53(3): 708-716.
[13] Wang F Y. 2009. 2. 5-D DC resistivity modeling by adaptive finite-element method with unstructured triangulation (in Chinese)[D]. Changsha: Central South University.
[14] Xiong B, Ruan B Y, Luo Y Z. 2003. 3-D numerical simulation study of DC resistivity anomaly under complicated terrain[J]. Geology and Exploration (in Chinese), 39(4): 60-64.
[15] Xu K J. 2007. Study on 2. 5D Complex Resistivity Electramagnetic Forward and Inversion (in Chinese)[Ph. D. thesis]. Jilin: Jilin University.
[16] Xu S Z, Duan B C, Zhang D H. 2000. Selection of the wavenumbers k using an optimization method for the inverse Fourier transform in 2. 5D electrical modeling[J]. Geophysical Prospecting, 48(5): 789-796, doi: 10. 1046/j. 1365-2478. 2000. 00210. x.
[17] Xu S Z. 1984. the boundary element method solving electric field of point source on two-dimensional geoelectrical section[J]. Journal of the Guilin College of Geology (in Chinese), 4(4): 119-133.
[18] Xu S Z. 1994. The finite element method in geophysics (in Chinese)[M]. Beijing: Science Press.
[19] Zhou B, Greenhalgh M, Greenhalgh M A. 2009. 2. 5-D/3-D resistivity modeling in anisotropic media using Gaussian quadrature grids[J]. Geophys. J. Int., 176(1): 63-80, doi 10. 1111/j. 1365-246X. 2008. 03950. x.
[20] 柳建新, 刘海飞. 2005. 计算最优化离散波数的优化算法[J]. 物化探计算与技术, 27(1): 34-38.
[21] 罗延钟, 孟永良. 1986. 关于用有限单元法对二维构造做电阻率法模拟的几个问题[J]. 地球物理学报, 29(6): 613-621.
[22] 罗延钟, 张桂青. 1987. 电子计算机在电法勘探中的应用[M]. 武汉: 武汉地质学院出版社.
[23] 强建科, 罗延钟. 2007. 三维地形直流电阻率有限元法模拟[J]. 地球物理学, 50(5): 1606-1613.
[24] 阮百尧. 2001. 三角单元部分电导率分块连续变化电源二维电场有限元数值模拟[J]. 广西科学, 8(1): 1-3.
[25] 王飞燕. 2009. 基于非结构化网格的2. 5-D直流电阻率法自适应有限元数值模拟[D]. 长沙: 中南大学.
[26] 熊斌, 阮百尧, 罗延钟. 2003. 复杂地形条件下直流电阻率异常三维数值模拟研究[J]. 地质与勘探, 39(4): 60-64.
[27] 徐凯军. 2007. 2. 5维复电阻率场正反演研究[博士论文]. 吉林: 吉林大学.
[28] 徐世浙. 1984. 点源二维地电剖面的边界单元解法[J]. 桂林冶金地质学院学报, 4(4): 119-113.
[29] 徐世浙. 1994. 地球物理中的有限单元法[M]. 北京: 科学出版社.