2. 中国地质大学地下信息探测技术与仪器教育部重点实验室, 北京 100083
3. 武警黄金第七支队, 烟台 264000
2. Key Laboratory of Geo-detection (China University of Geosciences), Ministry of Education, Beijing 100083, China
3. The 7th Gold Detachment of Chinese Armed Police Force, Yantai 264004, China
可控音频大地电磁测深法 (Controlled Source Audio-frequency Magnetotellurics) 简称CSAMT,在地质普查、勘探地热、金属矿产、水文、环境等各个领域应用广泛,发挥着重要的作用 (何继善,1990;吴璐苹等,1996;于昌明,1998;龚飞和底青云,2004;Fu et al., 2013).该方法具有工作效率高、勘探深度范围大、信噪比高等优点 (何继善,1990).根据源和采集方式的不同又可以分为标量、矢量和张量CSAMT等,理论研究和实际生产多以CSAMT标量采集处理为主.随着研究和应用的深入,CSAMT标量采集方式的缺陷也渐渐突出:该方法使用的是单一方向的有限长线源,只能在地下产生单一方向的电流矢量,适合探测水平均匀层状一维地质情况,但是实际条件大多很难满足这种假设,因此在数据处理和反演解释的时候会产生较大的偏差 (Boerner et al., 1993;Wannamakzr, 1997;Garcia et al., 2003;王显祥等,2014).
虽然大地电磁法 (MT) 采用天然场源,能够得到全阻抗张量信息,但是天然场源信号具有随机性,信噪比较低,在干扰较强的地区无法运用MT,而且天然场源在1~10 Hz和1000~3000 Hz频率范围活动水平很低,无法探测浅部地电特征 (陈乐寿和王光锷, 1990),为此有必要人工源主动激发.Caldwell等 (2002)指出,为了用电磁响应完全描述一个三维的地下电性信息,必须有至少两个极化方向的场源激发.
CSAMT张量数据采集、处理和解释很早就被提出,并应用于地质复杂构造区域的地质勘探,取得了良好效果.Li和Pedersen (1991)基于CSAMT方法首次明确提出了这种方法,证明了对于任意的源和接收机位置,任意的地球结构,阻抗张量和倾子矢量可以被唯一的确定,并研究了一维的张量和倾子向量的响应特征;Boerner (1993)等对加拿大纽芬兰布切恩斯矿山进行了CSAMT张量数据采集应用研究,识别出了帕拉因断层为布切恩斯地区逆冲构造史中重要的异常构造;Wannamaker (1997)在新墨西哥州Valles Caldera地热区域展开CSAMT张量数据实际测量和研究,结合钻井和地质资料进行了二维反演.
对单场源CSAMT的研究目前已经取得了不少进展 (Routh and Oldenburg, 1999;王若和王妙月,2003;Di et al., 2004;Weiss and Constable, 2006;张继峰等,2009;林昌洪等,2012),而多场源CSAMT张量采集由于发射和采集系统较复杂,涉及到多个场源的发射和多个场分量的接收,考虑源的CSAMT张量数据正反演的研究进展缓慢,研究较少.Li和Pedersen (1992)完成了各向异性张量可控源一维正演和反演 (Li et al., 2014),孟庆奎等 (2013)对张量CSAMT进行了一维数值模拟分析,分析了张量阻抗和倾子向量随频率的变化规律,胡英才等 (2015)实现了基于矢量有限元张量CSAMT三维正演.此外在硬件方面,德国Metronix公司2011年研制成功了首台满足张量阻抗观测的可控源大地电磁发射、采集系统①.
① 地球物理方法、仪器、应用 (第一辑), 北京欧华联科技有限责任公司.
为了满足电法勘探需求,在复杂地质构造地区实现高精度电法勘探,开发与CSAMT张量数据野外采集硬件系统匹配的数值模拟算法,本文采用三维有限差分法实现了CSAMT三维数值模拟,通过算例说明了CSAMT张量数据采集相比标量采集的优势所在.
1 CSAMT阻抗张量的计算阻抗张量是电磁测深中最重要的概念之一,它描述了地球系统的传输特性,包含了关于地下电性结构的信息 (尹兵祥和王光锷,2001).倾子向量主要反映了大地电性结构的水平横向变化 (陈小斌等,2004).张量CSAMT通过采集地表场的所有水平分量和垂直磁场分量,确定研究区域的阻抗张量和倾子向量.
CSAMT张量数据的获取需要两个或者两个以上不同极化方向的场源激发,场源可以为正交“十”字型源、“L“型源或者多个场源成一定角度同时激发 (何继善,1990).图 1为正交“十”字型场源张量CSAMT的激发和接收系统示意图,首先由第一个场源激发,在接受区域采集5个电磁场分量记为Ex1、Ey1、Hx1、Hy1和Hz1,然后由第二个场源激发,采集同样的5分电磁场分量记为Ex2、Ey2、Hx2、Hy2和Hz2,本文的研究基于此装置.
在CSAMT张量系统中,电场和磁场满足关系式为
(1) |
对于以上关系式,由采集的十个电磁场分量建立两个方程组,可以解得:
(2) |
按照下面的关系式可进一步计算出对应的卡尼亚电阻率ρ和相位ϕ,公式为
(3) |
其中i=x, y; j=x, y.
2 CSAMT三维正演地球物理中的数值模拟方法主要有有限单元法、有限差分法、积分方程法等,他们在解决三维问题中各有优缺点,本文考虑到计算效率、计算精度、地质构造复杂度等问题选择有限差分法进行正演,有限差分法也是目前解决三维问题应用最为广泛的一种方法 (谭捍东等,2003).
2.1 控制方程在CSAMT研究的频率范围内,可忽略位移电流的影响,地下介质的导磁率μ可近似为空气的导磁率μ0,在电性源激发情况下,电场强度矢量E和磁场强度矢量H相互作用的关系满足麦克斯韦方程组,取时间因子e-iωt,麦克斯韦方程组的微分形式为
(4a) |
(4b) |
(4c) |
(4d) |
上式中,i为虚数单位,ω为角频率,单位为rad·s-1,J为电流密度,单位为A·m-2,对 (4a) 式取散度,将 (4b) 式代入,考虑场源的影响,可以得到:
(5) |
上式中,将总场分解为一次场 (背景场) 和二次场 (散射场):
(6) |
其中一次场Ep为均匀或者层状背景介质下有限电偶源产生的电场,本文采用Kerry Key公开的CSAMT一维代码计算一次场 (Key,2009),二次场Es为异常体在一次场下感应产生的电磁场,则有
(7) |
其中σ0为背景电导率,整理 (5) 式有
(8) |
(8) 式即为含源的电场控制方程.
2.2 三维交错网格有限差分离散将控制方程 (8) 式在如图 2所示的Yee网格上进行离散化 (Yee,1966).Yee网格又被称为交错网格,它将电场分量定义在每个网格单元棱边的中心,将磁场分量定义在每个网格单元面的中心.Yee网格的电磁场分量的交错采样方式符合法拉第电磁感应定律和安培环路定律的自然结构,同时,这种电磁场的各分量的空间相对分布位置也非常适合麦克斯韦方程的差分计算,恰当地描述了电磁波的传播特性 (葛德彪和闫玉波,2011).
将矢量方程 (8) 式沿x、y、z三个方向分解,可以得到以下三个标量方程为
(9a) |
(9b) |
(9c) |
以y方向的标量方程 (9b) 式为例,用差分代替微分进行有限差分离散,如图 3,它是与y方向标量方程相关的电场分量分布关系 (Weiss and Constable, 2006),则 (9b) 式离散为
(10) |
上式中,
同样我们在x和z方向上离散就可以得到关于电场的三个离散方程组,除开边界上的电场分量,整个区域剖分的单元棱边上的所有电场组成未知数向量,可以得到一个大型线性方程组为
(11) |
K为大型稀疏系数矩阵,Es为待求二次电场值,S为与场源项和边界条件相关的右端项,无限远的二次电场值取零,因此右端项只与一次场有关.采用拟最小残差法 (QMR) 解该方程组,在迭代过程中引入散度校正.
由于CSAMT求得阻抗张量至少需要两个场源,本文均采取两个源计算,因此一次正演一个频率要解两个上述方程,公式为
(12) |
由三维正演解得的电场值通过插值得到研究区域的磁场值.则有:
(13) |
代入 (2) 式实现CSAMT三维阻抗张量的求解.
2.3 数值模拟精度分析在实现了张量CSAMT的三维正演之后,需要对该正演程序的精度进行评价.本文设计了一个二维模型,利用Kerry Key公开的2.5DCSEM有限单元正演程序计算结果与本文开发的三维正演计算结果进行对比验证.
模型如图 4所示,二维地质体沿x轴无限延伸,将两个正交有限长线源置于原点,源的长度为1000 m,分别平行于x轴和y轴.发射频率为:1、2.15、4.64、10、21.54、46.42、100、215.44、464.16、1000、2154.44、4641.59、10000 Hz,接收点r位于二维异常体上方,坐标为 (2000, 7880, 0),异常体的几何参数和电阻率参数如图.将两种方法正演得到的结果利用阻抗张量和倾子的计算公式 (2) 整理,得到的结果如图 5所示.
从图 5可以看出,四个分量Zxy、Zyx、Tzx和Tzy的实部和虚部吻合地非常好,本文的张量CSAMT三维数值模拟程序精度良好.
3 阻抗张量响应特征分析为了检验本三维程序的正演结果以及分析三维CSAMT阻抗张量数据的响应特征,本文设计模型进行三维正演,低阻异常体和有限长激发源的相对位置及尺寸如图 6所示,场源长度为500 m,场源1(TX1) 和场源2(TX2) 分别与x轴成45°和135°角.
利用式 (2) 整理得到各个频率各个方向的电阻率和相位的响应,以16 Hz和256 Hz为例,分析在地表z=0的平面等值线图特征如图 7和图 8所示.
从图 7和图 8中可以看出,对于低阻模型,不同的阻抗张量分量的响应特征具有较大的差异.阻抗张量反对角分量Zxy和Zyx,显示了低阻异常体的位置,反映了低阻模型在水平两个方向上低阻体与围岩的边界.Zxy的实部和虚部在整体上呈现高值,在异常体区域呈现低值,而且低值区域沿y方向拉伸,从图中标注的低阻异常体的范围可以看出其在x方向上具有较锐利的边界,而在y方向的边界则非常模糊;而Zxy的实部和虚部在整体上呈现低值,在异常体区域呈现高值,高值区域沿x方向拉伸,对比异常体可以看出其在x方向上的边界很模糊,而在y方向上的边界则非常锐利,这两个分量恰好地限制了异常体两个方向上得边界.
阻抗张量的主对角分量Zxx和Zyy的值和反对角的值相比很小,在研究区域内,Zxx和Zyy的实部、虚部在以异常体中心的四个象限内呈现高值或者低值的异常,如图 7所示.这两个阻抗张量分量显示了异常体在平面上的四个角点的位置.
倾子Tzx和Tzy反映了异常体的边界,如图 7所示,倾子Tzx的实部和虚部相比倾子Tzy的实部和虚部,其绝对值数量级较小,为10-2数量级,后者为10-1数量级,随着频率的增加,如图 8所示,倾子Tzy的总体数量级到10-2,而倾子Tzx的则不变,倾子Tzx的实虚部值在x方向上呈现了低阻异常体的边界,以异常体沿x轴的正方向和负方向的边界为界呈现高值和低值区域,其中实部的高值异常和低值异常呈闭合状,而虚部则不闭合;倾子Tzy的实虚部在y方向上呈现了低阻异常体的边界,以异常体沿y轴的正方向和负方向上的边界为界呈现高值和低值异常,其中实部的异常不闭合,虚部的异常呈现闭合形态.随着频率的增大,如图 8所示,两者的虚部分量异常形态也呈现闭合形态.
因此本文认为,增加一个源采集CSAMT张量信息,除了能够在x和y方向上限定异常体的边界外,张量数据还有其他的四个分量提供了其他的地电信息. Zxx和Zyy显示出异常体在水平方向上四个边界交点的位置,倾子Tzx和Tzy突出了异常体的在水平方向上的边界,通过CSAMT采集张量数据能够得到丰富的地电信息.
4 张量数据与标量数据对比分析对于单一场源CSAMT数据采集系统,只能采集关于一个场源的不同频率的电磁场值,通过整理可以得到标量数据,即卡尼亚电阻率和相位,其计算方法由式 (14) 和 (15) 表示;对于多个场源的CSAMT数据采集系统,能够采集到关于不同场源的电磁场值,通过整理可以得到张量数据,其计算方法由式 (2) 和 (3) 表示.公式 (14) 和 (15) 为
(14) |
(15) |
为了对比CSAMT张量采集和标量采集两种方式的优劣,本文根据图 6模型正演的结果,分别以源一、源二正演结果并依据式 (14)、(15) 整理出CSAMT标量数据卡尼亚电阻率和相位数据,同时根据式 (3) 整理出CSAMT阻抗张量对应的电阻率和相位,以16 Hz为例画出平面等值线图如图 9、图 10和图 11所示.
从图 9和图 10中我们可以看出,在异常体上方,两个源的四个电阻率均呈低阻异常,电阻率可低至50 Ω·m,四个相位呈现高值异常,可高至47°,由Ex和Hy场值整理的卡尼亚电阻率ρxy和相位ϕxy的异常形态与延伸趋势相似,而对于由Ey和Hx场值整理的卡尼亚电阻率ρyx和相位ϕyx也是如此.对比图 9和图 10可以发现,两者的异常走向是不同的,两者的异常走向恰好与激发该场的有限长源走向近似平行,因此对于CSAMT标量数据中整理出来的不论是卡尼亚电阻率或者相位,其异常的走向受源的激发方向的影响,而不仅仅取决于接收区域异常体的形态和走向.如果对比图 9或者图 10中两个不同方向上的视电阻率和相位,可以发现卡尼亚电阻率ρyx和相位ϕyx受源的激发方向影响更大.
图 11是由CSAMT张量数据中整理得到视电阻率和相位,从中可以发现从张量数据中整理得到的视电阻率和相位,其平面等值线异常则对应着地下的低阻异常,其异常形态受源的方向的影响较小.同时卡尼亚电阻率ρxy、ρyx及其对应的相位分别在水平x和方向上显示了低阻异常体的边界,因此对于16 Hz这个频率上来说,张量数据很好得反映了低阻异常体在水平面上的位置.
5 结论本文基于含源电场的控制方程,通过交错网格有限差分方法实现CSAMT张量数据的三维正演.设计二维模型进行正演与Kerry Key公开的二维CSAMT正演结果对比,验证了三维正演的正确性.设计具有代表性的三维模型进行正演,分析了单一低阻模型的CSAMT张量数据响应特征,在CSAMT的张量阻抗数据中,反对角阻抗主要在x和y两个方向上限定异常体的水平面上的边界,而主对角阻抗则辅助反对角阻抗提供平面上四个边界交点的位置,倾子向量含有电阻率沿水平变化信息.CSAMT标量数据由于只从一个源中得到,在水平面上只能限定异常体一个方向上的边界,通过采用张量源,能够取得丰富的地电信息;同时对于未知区域,张量采集能够尽可能地减少源的布设方向对于数据处理和解释的影响,提高解释精度.
致谢 感谢审稿专家提出的修改意见和编辑部的大力支持![] | Boerner D E, Wright J A, Thurlow J G, et al. 1993. Tensor CSAMT studies at the buchans mine in central newfoundland[J]. Geophysics, 58(1): 12–19. DOI:10.1190/1.1443342 |
[] | Caldwell G T, Bibby M H, Brown C. 2002. Controlled source apparent resistivity tensors and their relationship to the magnetotelluric impedance tensor[J]. Geophysical Journal International, 151(3): 755–770. DOI:10.1046/j.1365-246X.2002.01798.x |
[] | Chen L S, Wang G E. 1990. Magnetotellurics[M]. Beijing: Geology Publishing House. |
[] | Chen X B, Zhao G Z, Zhan Y, et al. 2004. Analysis of tipper visual vectors and its application[J]. Earth Science Frontiers, 11(4): 626–636. |
[] | Di Q Y, Unsworth M, Wang M Y. 2004. 2^5-D CSAMT modeling with finite element method over 2-D complex earth media[J]. Chinese Journal of Geophysics, 47(4): 825–829. DOI:10.1002/cjg2.v47.4 |
[] | Fu C M, Di Q Y, An Z G. 2013. Application of the CSAMT method to groundwater exploration in a metropolitan environment[J]. Geophysics, 78(5): B201–B209. DOI:10.1190/geo2012-0533.1 |
[] | Garcia X, Boerner D, Pedersen L B. 2003. Electric and magnetic galvanic distortion decomposition of tensor CSAMT data. Application to data from the Buchans Mine (Newfoundland, Canada)[J]. Geophysical Journal International, 154(3): 975–969. |
[] | Ge D B, Yan Y B. 2011. Finite-Difference Time-Domain Method for Electromagnetic Waves[M]. 3rd ed. Xi'an, China: Xi'an University of Electronic Science and Technology Press. |
[] | Gong F, Di Q Y. 2004. 1-D Simulation on model apparent resistivity curve of CSAMT method in certain mine[J]. Progress in Geophysics, 19(3): 631–634. DOI:10.3969/j.issn.1004-2903.2004.03.023 |
[] | He J S. 1990. Controlled Source Audio-frequency Magnetotellurics[M]. Changsha: Central South University of Technology Press. |
[] | Hu Y C, Li T L, Fan C S, et al. 2015. Three-dimensional tensor controlled-source electromagnetic modeling based on the vector finite-element method[J]. Applied Geophysics, 12(1): 35–46. DOI:10.1007/s11770-014-0469-1 |
[] | Key K. 2009. 1D inversion of multicomponent, multifrequency marine CSEM data:Methodology and synthetic studies for resolving thin resistive layers[J]. Geophysics, 74(2): F9–F20. DOI:10.1190/1.3058434 |
[] | Li X B, Pedersen B L. 1991. Controlled source tensor magnetotellurics[J]. Geophysics, 56(9): 1456–1461. DOI:10.1190/1.1443165 |
[] | Li X B, Pedersen L B. 1992. Controlled-source tensor magnetotelluric responses of a layered earth with azimuthal anisotropy[J]. Geophysical Journal International, 111(1): 91–103. DOI:10.1111/gji.1992.111.issue-1 |
[] | Li X B, Oskooi B, Pedersen L B. 2014. Inversion of controlled-source tensor magnetotelluric data for a layered earth with azimuthal anisotropy[J]. Geophysics, 65(2): 452–464. |
[] | Lin C H, Tan H D, Shu Q, et al. 2012. Three-dimensional conjugate gradient inversion of CSAMT data[J]. Chinese Journal of Geophysics, 55(11): 3829–3838. DOI:10.6038/j.issn.0001-5733.2012.11.030 |
[] | Meng Q K, Lin P R, Xu B L, et al. 2013. Study of one-dimensional numerical simulation of tensor CSAMT[J]. Computing Techniques for Geophysical and Geochemical Exploration, 35(4): 435–441. |
[] | Routh S P, Oldenburg W D. 1999. Inversion of controlled source audio-frequency magnetotellurics data for a horizontally layered earth[J]. Geophysics, 64(6): 1689–1697. DOI:10.1190/1.1444673 |
[] | Tan H D, Yu Q F, Booker J, et al. 2003. Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method[J]. Chinese Journal of Geophysics, 46(5): 705–711. DOI:10.3321/j.issn:0001-5733.2003.05.019 |
[] | Wang R, Wang M Y. 2003. Inversion method of controlled source audio-frequency magnetotelluric data[J]. Progress in Geophysics, 18(2): 197–202. DOI:10.3969/j.issn.1004-2903.2003.02.004 |
[] | Wang X X, Di Q Y, Xu C. 2014. Characteristics of multiple sources and tensor measurement in CSAMT[J]. Chinese Journal of Geophysics, 57(2): 651–661. DOI:10.6038/cjg20140228 |
[] | Wannamakzr P E. 1997. Tensor CSAMT survey over the Sulphur Springs thermal area, Valles Caldera, New Mexico, United States of America, Part Ⅰ:Implications for structure of the western caldera[J]. Geophysics, 62(2): 451–465. DOI:10.1190/1.1444156 |
[] | Weiss C J, Constable S. 2006. Mapping thin resistors and hydrocarbons with marine EM methods, Part Ⅱ modeling and analysis in 3D[J]. Geophysics, 71(6): G321–G332. DOI:10.1190/1.2356908 |
[] | Wu L P, Shi K F, Li Y H, et al. 1996. Application of CSAMT to the search for groundwater[J]. Acta Geophysica Sinica, 39(5): 712–717. |
[] | Yee K. 1966. Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media[J]. IEEE Transactions on Antennas and Propagation, 14(3): 302–307. DOI:10.1109/TAP.1966.1138693 |
[] | Yin B X, Wang G E. 2001. Canonical decomposition of magnetotelluric impedance tensor and geological meaning of the parameters[J]. Chinese Journal of Geophysics, 44(2): 279–288. DOI:10.3321/j.issn:0001-5733.2001.02.016 |
[] | Yu C M. 1998. The application of CSAMT method in looking for hidden gold mine[J]. Acta Geophysica Sinica, 41(1): 133–138. |
[] | Zhang J F, Tang J T, Yu Y, et al. 2009. Three dimensional controlled source electromagnetic numerical simulation based on electric field vector wave equation using finite element method[J]. Chinese Journal of Geophysics, 52(12): 3231–3141. DOI:10.3969/j.issn.0001-5733.2009.12.023 |
[] | 陈乐寿, 王光锷. 1990. 大地电磁测深法[M]. 北京: 地质出版社. |
[] | 陈小斌, 赵国泽, 詹艳, 等. 2004. 磁倾子矢量的图示分析及其应用研究[J]. 地学前缘, 11(4): 626–636. |
[] | 葛德彪, 闫玉波. 2011. 电磁波时域有限差分方法[M]. 3版.西安: 西安电子科技大学出版社. |
[] | 龚飞, 底青云. 2004. 某煤矿典型CSAMT法视电阻率曲线的一维模拟[J]. 地球物理学进展, 19(3): 631–634. DOI:10.3969/j.issn.1004-2903.2004.03.023 |
[] | 何继善. 1990. 可控源音频大地电磁法[M]. 长沙: 中南工业大学出版社. |
[] | 胡英才, 李桐林, 范翠松, 等. 2015. 基于矢量有限元法的三维张量CSAMT正演模拟[J]. 应用地球物理, 12(1): 35–46. |
[] | 林昌洪, 谭捍东, 舒晴, 等. 2012. 可控源音频大地电磁三维共轭梯度反演研究[J]. 地球物理学报, 55(11): 3829–3838. DOI:10.6038/j.issn.0001-5733.2012.11.030 |
[] | 孟庆奎, 林品荣, 徐宝利, 等. 2013. 张量CSAMT一维数值模拟分析[J]. 物探化探计算技术, 35(4): 435–441. |
[] | 谭捍东, 余钦范, BookerJ, 等. 2003. 大地电磁法三维交错采样有限差分数值模拟[J]. 地球物理学报, 46(5): 705–711. DOI:10.3321/j.issn:0001-5733.2003.05.019 |
[] | 王若, 王妙月. 2003. 可控源音频大地电磁数据的反演方法[J]. 地球物理学进展, 18(2): 197–202. DOI:10.3969/j.issn.1004-2903.2003.02.004 |
[] | 王显祥, 底青云, 许诚. 2014. CSAMT的多偶极子源特征与张量测量[J]. 地球物理学报, 57(2): 651–661. DOI:10.6038/cjg20140228 |
[] | 吴璐苹, 石昆法, 李荫槐, 等. 1996. 可控源音频大地电磁法在地下水勘查中的应用研究[J]. 地球物理学报, 39(5): 712–717. |
[] | 尹兵祥, 王光锷. 2001. 大地电磁阻抗张量正则分解的地质含义分析[J]. 地球物理学报, 44(2): 279–288. DOI:10.3321/j.issn:0001-5733.2001.02.016 |
[] | 于昌明. 1998. CSAMT方法在寻找隐伏金矿中的应用[J]. 地球物理学报, 41(1): 133–138. |
[] | 张继锋, 汤井田, 喻言, 等. 2009. 基于电场矢量波动方程的三维可控源电磁法有限单元法数值模拟[J]. 地球物理学报, 52(12): 3231–3141. DOI:10.3969/j.issn.0001-5733.2009.12.023 |