2014, Vol. 29
2. 东华理工大学, 核工程与地球物理学院, 南昌 330013
2. School of Nuclear Engineering and Geophysics, East China Institute of Technology, Nanchang 330013, China
复电阻率法(CR)也称频谱激电(SIP)是以岩、矿石的电化学性质差异为基础，通过观测10-3~103 Hz交变场源激发下大地的复电阻率频谱进行探测的地球物理方法.该方法野外工作与直流电阻率法相似；可用于解释的参数丰富，是较为公认的具有多参数反应异常能力的地球物理方法之一；早期研究学者就指出多频复电阻率法可以用来确定围岩响应、实现含矿异常识别( Zonge and Wynn，1975;Pelton et al., 1978).该方法主要用于金属矿、油气资源、地下水、有(无)机污染等的探测，随着对激发极化机理研究的不断深入，该方法还将不断扩大新的应用领域.
复电祖率正演模拟已经取得不少的研究成果，刘松阐述了复电阻率的理论基础和所面临的问题(刘崧，1998)；蔡军涛等开展了基于总场的二维地质体复电阻率有限元数值模拟(蔡军涛等，2007)；李勇等人开展了基于异常复电位的2.5维有限元数值模拟的方法及复电阻率反演研究(李勇等， 2011a，2011b).张辉开展了同时考虑激电效应、电磁效应的积分方程三维复电阻率正演(张辉等，2006).提高计算效率是复电阻率正演面临的主要问题之一，本文通过边界近似结合基于图论理论的矩阵重 排与填入元分析方法设计的直接解法，有效的提高了正演计算的速度.研究工作在面向对象的正反演统一框架下开展(张志勇和李曼，2012)，通过将Cole-Cole模型引入边值问题，得出复电阻率问题满足的变分，进而应用有限元法进行求解.
1 复电阻率2.5维在不考虑电磁耦合条件下，复电祖率正演计算与直流类似，不同之处需要引入Cole-Cole模型表示随频率变化的复电阻率.波数域的二维复电位边值问题为(Coggon，1971;徐世浙，1994)
 | 图 1 二维结构模拟 Fig. 1 Model of 2D structure |
对(3)求变分并令其为零，得到线性方程组
|
表 1 傅氏变换系数km和gm Table 1 km and gm parameters of Fouriers inverse transform |
地下由岩、矿石的激发极化效应产生的复电阻率的性质可以用Cole-Cole模型(罗延钟和方胜，1986;张桂青等，1987;阮百尧和徐世浙，1998； 阮百尧和罗润林，2003;蔡军涛等，2004)表示
统一框架是根据地球物理正演与反演特征，应用面向对象思想开发的一系列基础类.该框架正演部分主要包含模型剖分、单元插值、钢度矩阵生成及线性方程组求解等功能.其中默认剖为基于二叉树结构的收缩网格(张志勇和刘庆成，2013)(见图 2)，单元采用三角单元，可实现双线性与双二次插值；钢度矩阵合成对常见的单元面积分、边界积分形式进行了开发；线性系统求解分别开发了迭代解法与直接解法，其中迭代解法包括常用预处理(SOR，SSOR，不完全乔利斯基分解等)共轭梯度迭代方法，直接解法主要开发通过对矩阵的重排减少填入元的数量、通过符号分析确定分解矩阵元素分布，充分利用矩阵稀疏性的高效Cholesky分解(刘斌等，2012)及LDLT分解算法.
 | 图 2 有限单元网格剖分示意图 Fig. 2 Sketch map of finite element mesh |
在统一框架下，2.5维复电阻率的开发，只需在二维有限单元基类(CEFModel2D)派生新类(假定为CEF Model2DCR)；在派生类中设定单元物性数为四(即电 阻率、充电率、充电时间、频率相关系数)；加入Cole- Cole模型计算函数，计算单元复电阻率；重载单元场值梯度平方属性值函数(GetElementGradientSquare) 、单元场值平方属性值函数(GetElementFieldSquare)、边界值平方属性值函数(GetElementBoundaryField Square)；加入程序功能需要的其它函数，如电位计算、视电阻率振幅计算、相位计算、频散率计算等.
 | 图 3 复电祖率正演流程图 Fig. 3 Process of complex resistivity modeling |
分析钢度矩阵可知K1e与单元的复电阻率有关、K2e除与单元的复电阻率有关还与波数有关、K3e与单元的复电阻率和波数及供点电到边界的距离等有关，如果能构造近似的K3e使其只与供电点无关，那么对于所有的供电点只需要进行一次矩阵分解(Cholesky或LDLT)就可以通过回代完成所有供电点计算，这无疑将大大提高计算效率.对于有限单元来说，如果剖分形确定，则钢度矩阵的非零元素分布就可以确定，也就是只需通过一次符号分析就可为后续所有计算提供支持.研究工作采用了指数扩边，即将两侧边界及底边界迅速扩展(大于计算区10倍以上)，这样计算区供电点到边界的矢径可以用计算区中心到边界的矢径近似，则近似的K3e与供电点位置无关.
下图 4总结出了两种求解方案，方案1针对完全求解K3e而设计的，首先对频点循环，其次点源循环，然后波数循环，合成钢度矩阵K，对钢度矩阵K进行符号分析，采用Cholesky分解算法求解线性方程组，如图 4右边部分所示；方案2针对近似K3e而设计的，具体流程为：首先频点循环，其次波数循环，然后求解系数矩阵并对系数矩阵进行符号分析，最后点源循环，采用LDLT算法求解线性方程组，其流程如图 4左边所示.对比两方案可见方案2只进行了一次符号分析，对于相同的频率及波束进行了一次矩阵分解，求解效率远远高于方案1.为了比较两种方案，设置模型进行试算，表 2为100个不同供电点单个频点计算时间.本次计算采用的计算机CPU为Inter Core(TM)2主频2.33 GHz，内存3 G.
 | 图 4 有限单元数值解与一维解析解的对比 Fig. 4 Comparision of FEM numerical and one-dimensional analytical solutions |
|
|
表 2 不同算法计算的时间对比 Table 2 Contrast of calculation time |
4.1 均匀半空间模型
为了验证算法的正确性，设计均匀大地电阻率为ρ1=100 Ωm不极化介质模型，计算了理论电位曲线并与2.5D有限元算法计算的结果进行对比，计算结果如图 4，除了尾支由于边界的影响产生的误差较大外(相对误差2.67%)，曲线吻合很好，说明了本文提出的2.5D正演算法可靠.
4.2 单异常体模型建立如图 5的单异常体模型，分别进行高阻高极化与低阻高极化计算.取无极化围岩电阻率为ρ0=100 Ωm；低阻高极化异常体的Cole-Cole参数，ρ1=50 Ωm，m1=0.25，τ1=0.05，c1=0.25；高阻高极化异常体的Cole-Cole参数为：ρ2=200 Ωm，m2=0.25，τ2=0.05，c2=0.25.
 | 图 5 模型断面示意图 Fig. 5 Sketch map of the model |
采用对称四级装置，测量时极距AB=18.0 m(供电极距)、MN=6.0 m(接收极距)，模型体正上方视电阻率振幅和相位随着频率的变化，其计算结果如下图 6所示：
 |
图 6 频谱曲线 (a)振幅谱；(b)相位谱. Fig. 6 Spectrum curve |
从图 6中可以看出，复电阻率的振幅随着的频率的增加而逐渐减小，低阻体曲线变化更快表明低阻体表现出的激电性质比高阻体要高；相位曲线在各个频点上皆为负值，相位存在极小值且位置与充电时间的倒数对应.
4.2.2 Wenner-α装置复电阻率振幅与相位拟断面图 7a为低阻高极化模型四频点振幅拟断面图，图 7b为相位拟断面图，从图中可见，低阻体形成了明显的低阻异常，异常随着频率增加逐渐减弱，异常体的位置与异常极值对应.
 |
图 7 低阻高极化模型 (a)复电阻率振幅拟断面图；(b)复电阻率相位拟断面图. Fig. 7 Polarized low resistivity model |
图 8a为高阻高极化模型四频点振幅拟断面图，图 8b为相位拟断面图，从图中可见，高阻体形成了明显的高阻异常，异常随着频率增加逐渐减弱，异常体的位置与异常极值对应.
 |
图 8 高阻高极化模型 (a)复电阻率振幅拟断面图；(b)复电阻率相位拟断面图. Fig. 8 polarized high resistivity model |
图 9a、b中可以看出，频散率能够很好的反应异常的存在，异常体的位置与极值位置对应，低阻异常体频散率断面图相对高阻异常体频散率断面图异常幅值大.
 |
图 9 频散率断面图 (a)低阻高极化拟断面图；(b)高阻高极化拟断面图. Fig. 9 Pseudo-Section for Ps |
本文通过引入Cole-Cole模型表示地下介质的复电阻率性质，得出了2.5维复电阻率的边值问题，并在面向对象统一框架结构下进行了有限单元正演.通过边界积分的近似计算，利用基于符号分析的直接解法设计了正演的快速计算方案.进行了模型试算，试算表明程序计算正确、算法稳定高效.本文主要取得以下成果：
(1)面向对象的统一框架结构大大缩断了程序开发周期，提高了代码的重用性；
(2)近似边界及基于符号分析的直接解法适合于快速求解复电阻率问题，计算效率高，边界近似对计算精度影响不大；
(3)相同极化条件下，低阻高极化体较高阻阻高极化异常明显，由于介质的激电性质使电阻率振幅谱曲线随频率降低而升高、相位谱极值与时间常数对应；
(4)研究工作为考虑电磁耦合条件下的2.5D复电祖率正演及反演奠定了基础.
| [1] | Cai J T, Ruan B Y, Zhao G Z, et al. 2007. Two-dimensional modeling of complex resistivity using finite element method[J]. Chinese J. Geophys. (in Chinese), 50(6): 1869-1876, doi: 10.6038/cjg20071115. |
| [2] | Cai J T, Ruan B Y, Luo R L. 2004. Inverse calculation on apparent intrinsic induced polarization parametric curves over multi-layered earth[J]. Journal of Central South University: Natural Science (in Chinese), 35(4): 662-666. |
| [3] | Coggon J H. 1971. Electromagnetic and electrical modeling by the finite element method [J]. Geophysics, 36(1): 132-155. |
| [4] | Erdogan E, Demirci I, Candansayar M E. 2008. Incorporating topography into 2D resistivity modeling using finite-element and finite-difference approaches [J]. Geophysics (in Chinese), 73(3): F135-F142. |
| [5] | Johansen H K, Sorensen K. 1979. Fast Hankel transforms[J]. Geophysical Prospecting, 27(4): 876-901. |
| [6] | Li J M. 2005. Geoelectric Field and Electrical Exploration (in Chinese) [M]. Beijing: Geological Publishing House. |
| [7] | Li X B, Piao H R. 1998. Induced polarization and resistivity modeling of three-dimensional bodies in two-layered Earth using integral equation[J]. Acta Geophysica Sinica (in Chinese), 31(3): 342-352, doi:10.6038/cjg19880529. |
| [8] | Li Y, Lin P R, Li T L, et al. 2011a. Finite element method for solving anomalous complex potential of 2.5-D complex resistivity[J]. Journal of Jilin University (Earth Science Edition) (in Chinese), 41(5): 1596-1604, doi:10.13200378/j.cnki.jjuese.20110926. |
| [9] | Li Y, Lin P R, Zheng C J, et al. 2011b. Studies of complex resistivity forward simulation and inversion over a polarization Earth[J]. Journal of Jilin University (Earth Science Edition) (in Chinese), 41(S1): 355-361, doi:10.13278/j.cnki.jjuese.20110926. |
| [10] | Liu B, Li S C, Li S C, et al. 2012. 3D electrical resistivity inversion with least-squares method based on inequality constraint and its computation efficiency optimization[J]. Chinese J. Geophys. (in Chinese), 55(1): 260-268, doi:10.6038/j.issn.0001-5733.2012.01.025. |
| [11] | Liu S. 1998. Spectral Induced Polarization[M]. Wuhan: China University of Geosciences Press. |
| [12] | Luo Y Z, Fang S. 1986. An approximate inversion of the apparent complex resistivity spectrum [J]. Earth Science: Journal of China University of Geosciences (in Chinese), 23(1): 93-102 doi:10.1007/s11430-013-4685-3. |
| [13] | Pelton W H, Ward S H, Hallof P G, et al. 1978. Mineral discrimination and removal of inductive coupling with multifrequency IP [J]. Geophysics, 43(3): 588-609. |
| [14] | Ruan B Y, Xu S Z. 1998. Finite element method for modeling resistivity sounding on 2-D geoelectric model with line variation of conductivity within each block[J]. Earth Science-Journal of China University of Geosciences (in Chinese), 23(3): 303-307. |
| [15] | Ruan B Y, Luo R L. 2003. A new recursive inversion method of the complex-resistivity spectrum [J]. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 25(4): 298-301, doi:10.3969001 /j.issn.1001-1749.2003.11.04. |
| [16] | Wang D Y, Li T L, Li J P, et al. 2010. Forward simulation of a 3D complex resistivity model[J], Progress in Geophys. (in Chinese), 25(1): 266-271. doi:10.3969/j.issn.1004-2903.2010.01.035. |
| [17] | Weller A, Seichter M, Kampke A. 1996. Induced-polarization modelling using complex electrical conductivities[J]. , Geophys. J. Int.127(2): 387-398. |
| [18] | Xu S Z. 1988. Selection of wave number k in Fourier inverse transformation for point source and 2-D electric field problems [J]. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 10(3): 235-239, doi:10.3969/j.issn.1001-1749.1988.09.03. |
| [19] | Xu S Z. 1994. The Finite Element Method in Geophysics (in Chinese) [M]. Beijing: Science Press. |
| [20] | Yang X H, He J S, Tong X Z. 2008. Numerical simulation of frequency-domain IP with FEM[J]. Progress in Geophys.(in Chinese), 23(4): 1186-1189. |
| [21] | Zhang G Q, Cui X W, Luo Y Z. 1987. The determination of intrinsic parameters by inversing IP apparent spectrum [J]. Geology and Prospecting (in Chinese), 23(4): 48-54. |
| [22] | Zhang H, Li T L, Dong R X. 2006. Modeling electromagnetic responses of complex resistivity 3-D body using volume integral equation method[J]. Coal Geology & Exploration (in Chinese), 34(1): 70-74, doi:10.3969/j.issn.1001002-1986.2006.01.021. |
| [23] | Zhang Z Y, Liu Q C. 2013. Study on 2D MT Numerical Simulation Using FEM based on bi-tree Grid[J]. OGP (in Chinese), 48(3): 482-487. |
| [24] | Zhang Z Y, Li M. 2012. The framework designing of geophysical Modeling and inversion base on object-oriented method[J]. Progress in Geophys. (in Chinese), 27(3): 1207-1212, doi:10.6038/j.issn.1004-2903.2012.03.047. |
| [25] | Zonge K L, Wynn J C. 1975. Recent advances and applications in complex resistivity measurements [J]. Geophysics, 40(5): 851-864. |
| [26] | 蔡军涛, 阮百尧, 赵国泽,等. 2007. 复电阻率法二维有限元数值模拟[J]. 地球物理学报, 50(6): 1869-1876. doi: 10.6038/cjg20071115. |
| [27] | 蔡军涛, 阮百尧, 罗润林. 2004. 层状大地视真频参数测深曲线的反演[J]. 中南大学学报: 自然科学版, 35(4): 662-666. |
| [28] | 李金铭. 2005. 地电场与电法勘探[M]. 北京: 地质出版社. |
| [29] | 李晓波, 朴化荣. 1988. 两层大地中三维体的激发极化与电阻率响应的积分方程模拟[J]. 地球物理学报, 31(3): 342-352, doi:10.6038/cjg19880529. |
| [30] | 李勇, 林品荣, 李桐林,等. 2011a. 基于异常复电位2.5维CR有限元数值模拟[J]. 吉林大学学报(地球科学版), 41(5): 1596-1604, doi:10.13278/j.cnki.jjuese.20110926. |
| [31] | 李勇, 林品荣, 郑采君,等. 2011b. 极化大地复电阻率正反演[J]. 吉林大学学报(地球科学版), 41(增刊1): 355-361, doi:10.13278/j.cnki.jjuese.20110926. |
| [32] | 刘斌, 李术才, 李树忱,等. 2012. 基于不等式约束的最小二乘法三维电阻率反演及其算法优化[J]. 地球物理学报, 55(1): 260-268, doi:10.6038/j.issn.0001-5733.2012.01.025. |
| [33] | 刘崧. 1998. 谱激电法[M]. 武汉: 中国地质大学出版社. |
| [34] | 罗延钟, 方胜. 1986. 视复电阻率频谱的一种近似反演方法[J]. 地球科学: 中国地质大学学报, 23(1): 93-102, doi:10.1007/s11430-013-4685-3. |
| [35] | 阮百尧, 徐世浙. 1998. 电导率分块线性变化二维地电断面电阻率测深有限元数值模拟[J]. 地球科学-中国地质大学学报, 23(3): 303-307. |
| [36] | 阮百尧, 罗润林. 2003. 一种新的复电阻率频谱参数的递推反演方法[J]. 物探化探计算技术, 25(4): 298-301, doi:10.3969/j.issn.1001-1749.2003.11.04. |
| [37] | 王大勇, 李桐林, 李建平,等. 2010. 三维复电阻率模型电磁场正演模拟研究[J]. 地球物理学进展, 25(1): 266-271, doi:10.3969/j.issn.1004-2903.2010.01.035. |
| [38] | 徐世浙. 1988. 点电源二维电场问题中付氏反变换的波数k的选择[J]. 物探化探计算技术, 10(3): 235-239, doi:10.3969/j.issn.1001-1749.1988.09.03. |
| [39] | 徐世浙. 1994. 地球物理中的有限单元法[M]. 北京: 科学出版社. |
| [40] | 杨晓弘, 何继善, 童孝忠. 2008. 频率域激电有限元数值模拟[J]. 地球物理学进展, 23(4): 1186-1189. |
| [41] | 张桂青, 崔先文, 罗延钟. 1987. 一种反演频谱激电法视频谱求取真参数的方法[J]. 地质与勘探, 23(4): 48-54. |
| [42] | 张辉, 李桐林, 董瑞霞. 2006, 体积分方程法模拟复电阻率三维体电磁响应[J]. 煤田地质与勘探, 34(1): 70-74, doi:10.3969/j.issn.1001-1986.2006.01.021. |
| [43] | 张志勇, 刘庆成. 2013. 基于收缩二叉树结构网格剖分的大地电磁二维有限单元法正演研究[J]. 石油地球物理勘探, 48(3): 482-487. |
| [44] | 张志勇, 李曼. 2012. 基于面向对象思想的正演反演统一框架结构[J]. 地球物理学进展, 27(3): 1207-1212, doi:10.6038/j.issn.1004-2903.2012.03.047. |