摘要
本文就目前复电阻率法中的一些问题,开展了同时考虑激电和电磁感应效应的带地形的复电阻率2.5维电磁场正反演理论研究。首先将整个求解区域剖分为许多小单元,引入Cole-Cole模型表示各剖分单元的复电阻率,针对分块均匀的介质,从电磁场满足的麦克斯韦方程出发,推导出波数域中电场和磁场分量耦合的微分方程,利用基于有限元的双二次插值基函数和Galerkin法实现了起伏地表条件下电偶源的复电阻率2.5维正演模拟,并分别分析了地形参数和异常体复电阻率参数对正演结果的影响特征;然后推导出观测点处电磁场对于地下剖分各单元复电阻率参数的偏导数公式,利用偏导数公式并结合复电阻率2.5维电磁场正演算法及阻尼最小二乘最优化拟合技术,实现了起伏地表情况下电偶源的复电阻率2.5维电磁场反演。理论算例的结果表明,本文中所提出的带地形的复电阻率2.5维电磁场正反演方法在理论上是完全可行的。
The resistivity and induced polarization of geologic body generally show the characteristics of complex resistivity. Some scholars have expressed the complex resistivity by some mathematical models, for example the Cole-Cole model. It has been proved that the electrical anomaly of geologic body can be quantitative estimated by Cole-Cole model parameters (direct resistivityρ0, chargeability m , time constantτand frequency dependencec ). It is very important to study the methods to obtain the complex parameters of geologic body.
There are a lot of mountains in china, mountainous area take part in two over three. In order to enhance the veracity and reasonably of the complex resistivity data, we must consider the effect of topography. At present among data interpretation of electromagnetic method, the computation quantity of 3D forward problem and inversion is enormous. In the mountainous area, the efficiency of 3D forward problem and inversion cannot obtain the enhancement but the cost increase. The 2D forward problem and inversion have small computation quantity and achieve many achievements of research. But because of the source is 3D or finite, the 2D forward problem and inversion are only approximate simulation. Therefor, 2.5D mode(l3D or finite source over an earth whose conductivity varies in 2D)simulation have particular predominance. Compared with 3D model, 2.5D model reduces the computation quantity, which can be fulfilled easily on common computers. Compared with 2D model, it is easier for 2.5D model with topography to approach field geologic condition. so the research on electromagnetic forward and inversion of 2.5D model with topography has important theoretical and practical significance. In practical, electromagnetic effect and induced polarization exist at the same time and influence each other. The electromagnetic coupling will become the interference factor on induced polarization method, while the induced polarization will become the interference factor on frequency domain and time domain electric methods. so it is hard to separate them. Considering all these problems, we develop the research on electromagnetic forward and inversion of 2.5D complex resistivity with topography.
In this paper, we use the finite element method to compute the 2.5D electromagnetic response with topography, In finite element method solution, the solve region must be divided into small units, in this paper, we use the anomalous quadrangle grid to divide the region, and the realistic topography is realized by different altitude of node, furthermore the grid are small and denser in the middle area where the source and anomalous body distribute, while the grid are big and sparser at the edge of area. As a result, not only there is enough space for dividing area, but also complex parts of the geoelectric cross section gather in the center of dividing area. We deduce the frequency domain differential equation of electric field and magnetic field coupling according to Maxwell equation for this block homogeneous medium, then isoparametric elements and Galerkin method are used to disperse the differential equation, we obtain the linear equations which various nodes function value satisfies. Introducing Cole-Cole model into the 2.5D electromagnetic model,and introducing the pseudo-delta function to solve the source place irregularity. Moreover storage of system matrix employ the univariate nonzero element contraction which can greatly reduce the memory occupancy, the Cholesky preconditioned biconjugate gradient (ICBCG) method is used to solve the global system matrix equations. The integration is carried by continued fraction method in Fourier inverse transform.
By model calculating, this article analyses the effect of source frequency, topography, and Cole-Cole model parameters (ρ0, m ,τ,c ) on electric and magnetic field at survey site. It concludes that: 1. the value of horizontal electric field is changing from positive to negative gradually with source frequency; the real component of electric field is a positive anomaly above the anomalous body when the direct resistivity of anomalous body is high resistivity, the performance is a negative anomaly whenρ0is low resistivity, the imaginary part of electric field is opposite to real part; The real component of horizontal electric field decreases with the chargeability above anomalous body, the imaginary part of electric field increases with chargeability at left side of anomalous body, but decreases at the right; The horizontal electric field decrease with time-constant above anomalous body; The real component of horizontal electric field increases with the frequency dependence above anomalous body, but the imaginary part decreases with it; The abnormity of terrain increases with the frequency; When identical frequency the real part of electric field is bigger than the imaginary part. The amplitude of horizontal electric field at the middle of topography almost don’t change with direct resistivity, chargeability, time constant and frequency dependence. 2. the vertical magnetic field is a negative anomaly above the anomalous body when the direct resistivity of anomalous body is high resistivity, the performance is a positive anomaly whenρ0is low resistivity; The real component of vertical magnetic field decreases with the chargeability at left side of anomalous body, but increases at the right; the imaginary part of vertical magnetic field increases with chargeability at left side of anomalous body, but decreases at the right; The real part of vertical magnetic field increases with time-constant above anomalous body, while imaginary part decreases with time-constant; The real component of vertical magnetic field increases with the frequency dependence above anomalous body, but the imaginary part is opposite to real part. The amplitude of vertical magnetic field at the middle of topography almost don’t change with direct resistivity, chargeability, time constant and frequency dependence of the anomalous body.
As for the inversion of induced polarization parameter, this paper develops the partial derivative matrix that reflects the rate of real measurement field to the Cole-Cole model parameters by reciprocity theory. It is achieved that the complex resistivity 2.5D electromagnetic inversion of dipole source with topography by using classical damped least square method. In the inversion, the convergence rate of dc resistivity is fastest , the convergence rate of dependence is the second, and the convergence rate of time constant is slowest. On the other side, due to strong correlativity between time constant and chargeability, the time constant can’t be accurate inversed when all the complex resistivity parameters take part in the inversion. For this problem, considering the chargeability can be measured or estimated before the inversion calculation. So we fixed the values of the chargeability in the inversion, and then inversed other resistivity parameters. Through the method we can obtain accurate result of inversion. The results of theoretical calculation example indicate that the method about inversing resistivity parameters of 2.5D geologic body from electromagnetic information is feasible.
The production of research further developed the theory of forward and inversion for complex resistivity method. It has important theoretical and practical significance.
引文
[1] Snyder D D. A method for modeling the resistivity and IP response of two-dimensional bodies [J].Geophysics,1976,41(5):997-1015.
[2] Fox R C,Hohmann G W,Killpack T J,et al.Topographic effects in resistivity and induced-polarization surveys[J].Geophysics,1980,45(1):75-93.
[3]周熙襄,钟本善.电法勘探数值模拟技术[M].成都:科学技术出版社,1986.
[4]徐世浙.点电源二维电场问题中付氏变换的波数K的选择[J].物探化探计算技术,1988,10(3):235- 239.
[5]罗延钟.关于用有限元法对二维构造作电阻率法模拟的几个问题[J].地球物理学报,1986,29(6): 613-621.
[6] Coggon J H. Electromagnetic and electrical modeling by the finite-element method[J]. Geophysics, 1971,36(2):132-155.
[7] Stoyer C H, Greenfield R J. Numerical solutions of the response of a two-dimensional earth to an oscillating magnetic dipole source[J].Geophysics,1976,41(3):519-530.
[8] Lee K H. Electromagnetic scattering by a two-dimensional in homogeneity due to an oscillating magnetic dipole[D].Berkeley:Univ. of California,1978.
[9] Lee K H, Morrison H F. A numerical solution for the electromagnetic scattering by a two-dimensional in homogeneity[J].Geophysics,1985,50(3):466-472.
[10] Doherty J. EM modeling using surface integral equations[J].Geophys.Prosp.,1988,36:644-668.
[11] Unsworth M J, Traves B J, Chave A D. Electromagnetic induction by a finite electric dipole source over a 2-D earth[J].Geophysics,1993,58(2):198-214.
[12] Torres V C, Habashy T M. Rapid 2.5D forward modeling and inversion via a new nonlinear scattering approximation[J].Radio science,1994,29(3):1051-1079.
[13] Mitsuhata Y. 2-D electromagnetic modeling by finite-element method with a dipole source and topography[J].Geophysics,2000,65(2):465-475.
[14] Yu guo Li, Kerry Key. Two-dimensional marine controlled-source electromagnetic modeling:Part 1 An adaptive finite-element algorithm[J].Geophysics,2007,72(2):WA51-WA62.
[15]孟永良,罗延钟.电偶源CSAMT法二维正演的有限单元算法[A].见:罗延钟编.勘查地球物理勘查地球化学文集,第20集[C].北京:地质出版社,1996,103-114.
[16]毛先进,鲍光淑. 2.5维问题电阻率正演的新方法[J].中南工业大学学报,1997,28(4):307-310.
[17]底青云,MartynUnsworth,王妙月.复杂介质有限元法2.5维可控源音频大地电磁法数值模拟[J].地球物理学报,2004,47(4):723-730.
[18]底青云,MartynUnsworth,王妙月.有限元法2.5维CSAMT数值模拟[J].地球物理学进展,2004,19(2): 317-324.
[19] Newman G A, Hohmann G W. Transient electromagnetic responses of high-contrast prisms in alayered earth[J].Geophysics,1988,53:691-706.
[20] Leppin M. Electromagnetic modeling of 3D sources over 2D in homogeneties in the time domain[J]. Geophysics, 1992, 57: 994-1003.
[21] Everett M E, Edwards R N. Transient marine electromagnetics: the 2.5D forward problem[J]. Geophys J Int,1993,113:545-561.
[22] Mitsuhata Y.Computation of the electrical responses of mid-ocean ridged structures[D].Toronto: University of Toronto,1994.
[23]孟永良,罗延钟.二维地电构造上谐变电偶极子电磁场有限元算法[A].1994年中国地球物理学会第十届学术年会论文集[C],1994.
[24]王华军,罗延钟.中心回线瞬变电磁法2.5维有限单元算法[J].地球物理学报, 2003,46(6):855 -862.
[25]熊彬,罗延钟.电导率分块均匀的瞬变电磁2.5维有限元数值模拟[J].地球物理学报, 2006,49 (2): 590-597.
[26] Kong F N, Johnstad S E, Rosten T,et al. A 2.5D finite-element-modeling difference method for marine CSEM modeling in stratified anisotropic media[J]. Geophysics, 2008,73(1): F9-F19.
[27] Newman G A. Croswell electromagnetic inversion using integral and differential equations[J]. Geophysics, 1995,60:899-911.
[28] Torres V C, Habashy T M. Rapid 2.5D forward modeling and inversion via a new nonlinear scattering approximation[J]. Radio science,1994, 29(3): 1051-1079.
[29] Unsworth M J, Oldenburg D W. Subspace inversion of electromagnetic data: applicatinon to mid-oceanridge exploration[J]. Geophys J.Int,1995, 123: 261-168.
[30] Lu X. Rapid relaxation inversion of CSAMT data [J]. Geophys J Int, 1999, 138(1): 381-392.
[31]毛先进,鲍光淑. 2.5维电阻率成像的新方法[J].物探与化探, 1999,23:150-152.
[32] Mitsuhata Y,Uchida T,Amano H. 2.5-D inversion of frequency-domain electromagnetic data generated by a grounded-wire source[J].Geophysics, 2002,67(6): 1753-1768.
[33]底青云,MartynUnsworth,王妙月.2.5维有限元法CSAMT数值反演[J].石油地球物理勘探, 2006,41(1):100-106.
[34] Shen J S, Sun W B, Zhao W J, et al. A 2.5D cross-hole electromagnetic modeling and inversion method and its application to survey data from the Gudao oil field,east China[J].Journal of Geophysics and engineer,2008,5(4):401-411.
[35] Pelton W H. Interpretation of complex resistivity and dielectric data[D].Salt Lake City,Univ. of Utah, 1977.
[36] Pelton W H, Ward S H, Hallof P G., et al. Mineral discrimination and removal of inductive coupling with multrfrequency IP [J]. Geophysics, 1978,43(3): 588-609.
[37] Zonge K L, Hughes L J, Carlson N R. Hydrocarbon exploration using induced polarization, apparent resistivity and electromagnetic scattering[C]. The fifty-first annual international meeting and exposition SEG, 1981.
[38] Wong J. An electrochemical model of the IP phenomenon in disseminated sulfide ores [J]. Geophysics,1979,44:1245-1265.
[39] Major J, Silic J. Restrictions on the use of Cole-Cole dispersion models in complex resistivity interpretation[J]. Geophysics, 1981,46(6): 916-931.
[40] Sunde E D. Earth conduction effects in transmission systems[M].Van Nostrand, New York, N.Y.,1949.
[41] Dey A, Morrison H F. Electromagnetic coupling in frequency and time-domain induced polarization surveys over a multilayered earth [J]. Geophysics, 1973,38(2): 380-405.
[42] Guptasarma D. Effect of surface polarization on resistivity modeling[J]. Geophysics, 1983,48:98-106.
[43] Guptasarma D. True and apparent spectra of buried polarizable targets[J].Geophysics,1984,49:1534- 1540.
[44] Soininen H. The behavior of the apparent resistivity phase spectrum in the case of a polarizable prism in an unpolarizable half-space[J].Geophysics,1984,49:1534-1540.
[45] Soininen H. The behavior of the apparent resistivity phase spectrum in the case of two polarizable media [J]. Geophysics, 1985,50:810-819.
[46] Brown R J. EM coupling in multi-frequency IP and a generalization of the Cole-Cole impedance model[J]. Geophysical prospecting,1985,33:282-302.
[47] Dias C A. Developments in a model to describe low-frequency electrical polarization of rocks [J].Geophysics, 2000, 65:437-451.
[48]张赛珍,李英贤,张树椿,等.我国几个金属矿区岩、矿石的低频相位频率特性及影响因素[J].地球物理学报, 1984,27:176-189.
[49]罗延钟,张桂青.频率域激电法原理[M].北京:地质出版社,1988.
[50]罗延钟,吴之训.研究剩余电磁效应的理想参数[J].物探与化探,1992,16:370-375.
[51]罗延钟,吴之训.谱激电法中频率相关系数的应用[J].地球物理学报, 1992,35(4):490-499.
[52]柯马罗夫B A著,阎立光等译.激发极化电法勘探[M].北京:北京地质出版社,1983.
[53]罗延钟,许妙立,吴之训.面极化球体上复电阻率的频谱特征和模拟相似性准则[J].桂林冶金地质学院学报,1986,6:145-155.
[54] Xiong Z, Luo Y, Wang S,et al. Induced-polarization and electromagnetic modeling of a three-dimensional body in a two-layer anisotropic earth[J]. Geophysics,1986,51: 2235-2246.
[55] Meng Y L, Luo Y Z. Calculation of harmonic electromagnetic field and scale modeling laws forsurface polarizable bodies[J], Chinese journal of Geophysics, 1989,32:145-158.
[56] Liu Song, Vozoff K. The complex resistivity spectra of models consisting of two polarizable media of different intrinsic properties [J]. Geophysical prospecting, 1985, 33:1029-1062.
[57]罗延钟,方胜.极化椭球体上频谱激电法的异常形态[J].桂林冶金地质学院学报,1985, 5:49-63.
[58]傅良魁.复电阻率法异常的频谱及空间分布规律[J].地质与勘探,1981,2:46-53.
[59]傅良魁,张虎豹.频谱激电法中的若干理论研究[J].物探与化探,1985,9:410-419.
[60]傅良魁,李金铭,陈兆洪.埋藏极化体激电时间谱视参数的试验研究结果[J].地质与勘探, 1985,12:38-43.
[61]李金铭,陈兆洪.埋藏极化体时间域激电谱的理论研究[J].桂林冶金地质学院学报,1987,7: 295-305.
[62]张辉,李桐林,董瑞霞.积分方程法模拟复电阻率三维电磁响应[J].煤田地质与勘探,2006,34(1):70- 74.
[63]张辉,李桐林,董瑞霞,等.利用高斯求积和连分式展开计算电磁张量格林函数积分[J].地球物理学进展,2005, 19(3):667-670.
[64]张辉,李桐林,董瑞霞.体积分方程法模拟电偶源三维电磁响应[J].地球物理学进展,2006,21(2):386- 390.
[65]张辉.复电阻率三维电磁场正反演研究[D].长春:吉林大学地球探测科学与技术学院,2006.
[66]罗延钟,方胜.视复电阻率频谱的一种近似反演方法[J].地球科学,1986,11:93-102.
[67]张桂青,崔先文,罗延钟.一种反演频谱激电法视谱直接求取真参数的方法[J].地质与勘探, 1985,23:48-54.
[68]张桂青,罗延钟,崔先文.由视谱直接反演真谱参数的实用方法[J].地球科学,1990,16: 229-239.
[69]张桂青.Cole-Cole模型复频谱的最优化反演[A].见:傅良魁主编.电法勘探文集[C].北京:地质出版社,1986,233-245.
[70]王自力,张赛珍.一种真复电阻率谱参数的求解方法[J].地球物理学报,1990,33(6): 712-727.
[71]刘崧,官善友,高鹏飞.求极化椭球体真Cole-Cole参数的联合谱激电反演[J].地球物理学报,1994,37(S2):542-551.
[72]刘崧,薛继安,徐建华.求地下真柯尔-柯尔参数的研究[J].物探与化探,2000,24:51-61.
[73] Xiang J, Jones N B, Cheng D,et al. Direct inversion of the apparent complex-resistivity spectrum [J]. Geophysics, 2001,66(5): 1399-1404.
[74]阮百尧,罗涧林.一种新的复电阻率频谱参数递推反演方法[J].物化探计算技术,2003,25(4):298- 301.
[75]徐凯军.2.5维复电阻率电磁场正反演研究[D].长春:吉林大学地球探测科学与技术学院,2007.
[76]李建平.带地形的复电阻率三维电磁场正反演研究[D].长春:吉林大学地球探测科学与技术学院,2008.
[77]齐基威茨O C,邱Y K.结构和连续力学中的有限单元法[M].北京:国防工业出版社,1967.
[78] Rodi W L. A technique for improving the accuracy of finite element solutions for magnetotelluric data[J]. Geophys.J.R.astr.SOC,1976,44:483-506.
[79] Rijo L. Modeling of electric and electromagnetic data[D]. Salt Lake City, Univ. of Utah, 1977.
[80] Kaikkonen P. Numerical VLF modeling[J].Geophysical Prospecting, 1979,27: 106-136.
[81] Dey A, Morrison H F. Resistivity modeling for arbitrarily shaped three-dimensional structures[J]. Geophysics,1979,44:753-780.
[82] Pridmore D F, Hohmann G W, Ward S H,et al. An investigation of finite-element modeling for electrical and electromagnetic data in three dimensions[J].Geophysics,1981,46(7):1009-1024.
[83] Wannamaker P E, Stodt J A, Rijo L. Two-dimensional topographic response in magnetotellurics modeled using finite element[J].Geophysics,1986,51(11):2131-2144.
[84]朱伯芳.有限单元法原理与应用[M].北京:水利出版社,1979.
[85]罗延钟,张桂青.电子计算机在电法勘探中的应用[M].武汉:武汉地质学院出版社,1987.
[86]徐世浙,赵生凯.二维各向异性地电剖面的大地电磁场的有限单元解法[J].地震学报,1985a, 7(7): 80-90.
[87]徐世浙.点源二维各向异性地电断面直流电场的有限元解法[J].中东海洋学院学报, 1988a,18(1): 81-90.
[88]徐世浙.电导率分层线性变化的水平层的点电源电场的数值解[J].地球物理学报, 1986,29 (1):84- 90.
[89]徐世浙.地球物理中的有限单元法[M].北京:科学出版社,1994.
[90]徐世浙.二维分块均匀物体重力异常的计算[J].中国科技大学学报,1984,14: 126-132.
[91]徐世浙.有限元法及其在物探中的应用简介[J].物探化探计算技术, 1982,4:86-103.
[92]金建铭著,王建国译.电磁场有限元方法[M].西安:西安电子科技大学出版社,1998.
[93]李大潜.有限元素法在电法测井中的应用[M].北京:石油工业出版社,1980.
[94]周刚.人工源电磁法二维半问题等参有限元法数值模拟研究[D].北京:中国地质大学,2006.
[95]王胜阁.地形起伏条件下可控源音频大地电磁法2.5维正演研究[D].北京:中国地质大学,2007.
[96]熊彬,罗延钟.电导率分块均匀的瞬变电磁2.5维有限元数值模拟[J].地球物理学报,2006,49(2):590- 597.
[97]孟永良,罗延钟,昌彦君.时间谱电阻率法的二维正演算法[J].地球科学-中国地质大学学报, 2000, 25(6):656-662.
[98]阮百尧,熊彬,徐世浙.三维地电断面电阻率测深有限元数值模拟[J].地球科学,2001,26 (1):73-77.
[99]阮百尧,熊彬.电导率连续变化的三维电阻率测深有限元数值模拟[J].地球物理学报,2002,45: 131-138.
[100]黄俊革.三维电阻率/极化率有限元正演模拟与反演成像[D].长沙:中南大学,2003.
[101]林莘,刘志刚.ICCG算法在SF6罐式高压断路器三维电场有限元计算中的应用[J].中国电机工程学报,2001,21(2):21-24.
[102]宁石英.ICCG算法在解复线性方程组中的应用[J].哈尔滨电工学院学报,1990,13(2):213- 218.
[103]文代刚.求解有限元复代数方程组的实型ICCG法[J].电工技术学报,1996,11(4): 61-64.
[104]吴小平,徐果明.不完全Cholesky共轭梯度法及其在地电场计算中的应用[J].石油地球物理勘探, 1998,33(1):90-94.
[105]蒋豪贤,何志伟,陈强,等.基于加权余量法的地网接地电阻的数值计算[J].华南理工大学学报(自然科学版),1996,24(1),173-180.
[106]徐士良.FORTRAN常用算法程序集[M].北京:清华大学出版社,1995.
[107]徐世浙.点电源二维电场问题中付氏变换的波数K的选择[J].物探化探计算技术, 1988,10(3): 235-239.
[108]柳建新,刘海飞.计算最优化离散波数的优化算法[J].物探化探计算技术,2005,27(1):34-38.
[109]张献民,陈文华.点源二维势场问题的边界元法中点源处理方法[J].地质与勘探,1989,25(10): 37-42.
[110]刘崧.谱激电法[M].武汉:中国地质大学出版社,1998.
[111]Wait J R. Overvoltage research and geophysical applications[M].London,Pergamen Press,1959.
[112]Ward S H, Fraser D C. Conduction of electricity in rocks[J].In soc. Expel. Geophys. Ed. Comm. Mining Geophys, 1967, 197~223.
[113]Madden T R, Cantwell T. Induced polarization:A review[J]. In soc. Expel. Geophys. Ed. Comm. Mining Geophys, 1967, 373~400.
[114]Dias C A. Analytical model for a polarizable medium at radio and lower frequencies [J]. Geophys. Res.,1972,77:4945-4956.
[115]Zonge K L,Hughes L J. Comparison of time,frequency,and phase measurements in IP[J].Geophysical Prospecting,1972,20:626-648.
[116]刘崧.激电法在寻找深埋金属矿中的应用前景及方法理论中值得研究的一些问题[J].地质科技情报,1985.
[117]张辉,李桐林.利用cole-cole模型组合得到SIP真参数的联合频谱最优化反演[J].西北地震学报,2004,26(2):108-112.
[118]Marquardt D W. An algorithm for least-squares estimation of non-linear parameters [J]. J. Soc. Indust. Appl. Math,1963,11:431-441.
[119]Marquardt D W. Generalized inverses ridge regression, biased linear estimation and nonlinear estimation,Technometrics,1970,12:591-612.
[120]姚姚.地球物理反演基本理论与应用方法[M].北京:中国地质大学出版社,2002.
[121]Patricia P L, Wannamaker P E. Calculating the two-dimesional magnetotelluric jacobian in finteelements using reciprocity[J]. Geophys.J.Int,1996,127:806-810.
[122]McGillivray P R., Oldenburg D W, Ellis R G,et al. Calculation of sensitivities for the frequency- domain electromagnetic problem[J].Geophys.J.Int.,1994,116:1-4.
[123]Farquharson C G, Oldenburg D W. Approximate sensitivities for the electromagnetic inverse problem [J].Geophys.J.Int.,1996,126:235-252.
[124]Parasnis D S. Reciprocity theorem in geoelectric and geoelectromagnetic work[J]. Geoexpl.,1988, 25(3):177-198.
[125]Zienkiewicz O C,Taylor R L. The finite element, 6th ed.[M].New York:McGraw-Hill,2005.