 Advances in Geosciences 地球科学前沿, 2012, 2, 246-251 http://dx.doi.org/10.12677/ag.2012.24037 Published Online December 2012 (http://www.hanspub.org/journal/ag.html) The Finite Element Method for Solving Magnetotelluric Field Boundary-Value Problem on Anisotropy Media* Yu Shi1, Tian ya Luo2, Bin Xiong2#, Ch a ngwei Li2, Hong Yang2, Yezhong Huang2 1Northeast Team, Bureau of Geological Exploration & Development of Jiangxi Province, Shangrao 2College of Earth Sciences , Guilin University of Technology, Guilin Email: shiyusyq@163.com, #xiongbin@msn.com Received: Oct. 10th, 2012; revised: Nov. 4th, 2012; accepted: Nov. 24th, 2012 Abstract: Based on absorption of previous work, this paper firstly derives the variational prob lem of magnetotellu ric in 2-D anisotropic media and provides the numerical solution calculated by finite element method, then computes apparent resistivity. At last, discusses the characteristics of magneto telluric field v arying with frequency. Keywords: Magnetotellurics; FEM; Anisotropy Medium; Forward Modeling 有限元法求解各向异性介质中大地电磁场边值问题研究* 石 宇1, 罗天涯 2，熊 彬2#，李长伟 2，杨 红2，黄业中 2 1江西省地质矿产勘查开发局赣东北大队，上饶 2桂林理工大学地球科学学院，桂林 Email: shiyusyq@163.com, #xiongbin@msn.com 收稿日期：2012 年10 月10 日；修回日期：2012 年11 月4日；录用日期：2012 年11 月24 日 摘 要：在吸收前人工作基础上，文中给出了二维各向异性介质中大地电磁场的变分问题以及有限单元法数值 解，并计算了视电阻率，由其随频率的变化规律讨论电磁场的变化特征。 关键词：大地电磁；有限单元法；各向异性介质；正演 1. 引言 随着地球结构和物质组成的研究发展，人们逐渐 认识到电导率随方位变化引起的地球介质电性各向 异性对于解释大地电磁场数据显得越来越重要[1]，国 内外的学者做了很多的研究，Josef Pek等[2]研究了一 维各向异性层状介质大地电磁阻抗和参数灵敏度；A. M. Osella等[3]研究了各向异性二维结构大地电磁 响 应；电性各向异性的产生已经成为很多文章热议的主 题，例如，深部地壳研究(Everet and Constable, 1999)； 破裂勘查和填图(Le Masne and Vasseur, 1981)；矿 产 勘 查(Al -Garni and Everet t, 2003)；测 井 曲 线 (Lu and Alum- baugh, 2001)；Barber (2004)等证明了怎样通过测井曲 线确定电性各向异性，Ellis(2010)等通过选择晶粒形 状和排列为沉积各向异性的出现提出一个有效的介 质模型等[4]；Reddy 等[5](1971)讨论的浸渍各向异性的 大地电磁测深效应；Reddy 等[6]在1975 年进一步研究 了横向不均匀各向异性介质中的大地电磁测深响应； 赵生凯、徐世浙[7](1983)研究了有限单元法计算良导嵌 入体二维模型的大地电磁测深曲线 ；徐世浙，赵生凯[8] 在1985 年对二维各向异性地电断面大地电磁场的有 限元解法进行了深入地研究；陈乐寿，王光愕[9](1990) 推导出了层状各向异性介质中视电阻率的递推关系 公式；林长佑等[10]探讨了水平层状对称各向异性 *基金项目：国家自然科学基金项目(40974077、41164004)，广西自 然科学基金项目(2011GXNSFA018003)，广西高校优秀人才资助计 划(桂教人[2010]41 号)。 #通讯作者。 Copyright © 2012 Hanspub 246  有限元法求解各向异性介质中大地电磁场边值问题研究 介质的大地电磁资料反演，并用于识别地震的深部电 性各向异性变化前兆；杨长福[11](1997)对二维 MT 各 向异性正演模拟进行了研究；J. Pek等[12]在1997 年讨 论了二维各向异性介质中大地电磁场的有限差分法； 刘云，王绪本在 2010 年对大地电磁二维自适应地形有 限元正演模拟进行研究[13]等。从前人的研究可以看 出，对各向异性介质研究的学术价值较大。基于前人 的研究，文中就各向异性介质中大地电磁场的异常特 征进行讨论。 2. 二维各向异性介质大地电磁场的 边值问题 二维各向异性介质 MT 两种极化模式下满足的边 值问题[14]为 0 uu (1a) 1 ABu (1b) 0 AD,BC u n (1c) 0 CD uku n (1d) 对于 TE 波 z uE, i , 1 i , (2) 研究区域包括空气和地下，上边界 AB取在距离地表 足够远的位置，AD、BC 分别为距离局部不均匀体足 够远的左右边界，CD 为距离地表足够远的下边界(图 1)。对于 TM波 z uH, i , 11 12 21 22 , (3) 22 11 11 sincos ii , (4a) 22 22 11 cos sin ii , (4b) 12 2111 1sin 2 2ii , (4c) 上边界 AB直接取在地表，下边界以及左右边界同 TE 波(图2)式中 为二维哈密顿算符，ki ， 为研究区域， 、 分别为平行层面和垂直层面的 Figure 1. Study region of TE 图1 . TE波研究区域 Figure 2. Study region of TM 图2 . TM波研究区域 电导率， 为介电常数， 为介质的磁导率，i为复数 单位， 为定态电磁波的角频率， 为物性界面与水 平面的夹角。 3. 二维各向异性介质大地电磁场的 变分问题 与边值问题对应的变分问题为 22 CD 111 dd 222 Fuu uugku (5a) AB 1u (5b) 0Fu (5c) 对于 TE 波， 1 gi ，对于 TM 波， 1 gi 。 4. 有限单元法 文中采用的区域剖分方式为，通过引矩形单元的 Copyright © 2012 Hanspub 247  有限元法求解各向异性介质中大地电磁场边值问题研究 Copyright © 2012 Hanspub 248 两条 (6) 这里，下标 i，j，m是三角单元按 列的 的区域积分和外边界积 分分解为 各单 元积 22 CD 2 2 CD CD 111 dd 222 11 dd 22 1d 2 ee Fuu uugku uu u gku (7) 对角线，将其细分为四个小三角形[14-16]如图 3所 示。在三角单元 e上，将u表为线性插值函数[14] iijjmm uNuNu Nu 逆时针方向排 三角节点编号如图4所示，Ni，Nj，Nm为线性插 值的基函数。 将(5a)式中 对(7)式各个单元和边界逐一积分，并消去中间节点， 得到 4 × 4的矩阵，将它们扩展成为全体节点的矩阵， 再对扩展矩阵求和，之后对泛函取极值，得到线性代 数方程组，在解之前，代入 AB 线上的边界值，采用 文献[14]的解法即可得到各节点的 z E或 z H ，计算辅 助场，进而可以得到地表的视电阻率。 分之和 5. 正演结果分析 5.1. 分层均匀介质模型 Figure 3. Triangle element in the rectangle eleme nt 图3. 矩形单元中的三角单元 给出具有解析解的各向同 性层状介质模型 ： 11000 m ，11080 mh ；2100 m ， 21000 mh ；31000 m ，其正演结果如图 5所 示。根据正演结果数据的误差分析得：TE 波、TM 波 与解析解的均方相对误差均为 0.2%，这验证了本算 法对各向同性介质计算的可靠性；由于各向异性大 地电磁异常场的解析解目前还未见到，因此文中仅 对几个典型模型进行数值计算。 i j m Figure 4. Triangle element 图4. 三角单元 (a) (b) Figure 5. Comparison between anaal and numerical solution by FEM 图5. 有限元数值解与解析解对比 lytic  有限元法求解各向异性介质中大地电磁场边值问题研究 5.3. 倾斜高阻体模型 5.2. 向斜良导体模型 ，其正演结果如图7所示。 从图中可知， 1) 在低频段，大地电磁异常场主要受到 深部 倾斜高阻体模型，其视电阻 率拟断面图如图 9所示。从图中可以清楚地看出很明 如图 6所示向斜模型 向斜良导体的影响，随着频率的增大，由异常体 引起的大地电磁异常场慢慢变弱，曲线逐步过渡到反 映向斜良导体上方各向同性均匀介质的视电阻率；2) 在各向同性良导体的电导率等于各向异性良导体的 平均电导率的情况下，二者呈现出不同的异常幅度， 各向同性体的异常幅度约为24 倍，而各向异性体的 异常幅度约为 85 倍。可见，将各向异性介质近似为 各向同性来解释将会造成较大的偏差甚至是错误的 结果。 如图 8所示各向异性 Figure 6. Syncline model 图6. 向斜模型 (a) (b) (c) Figure 7. Comparison of the magnetotelluric profile between isotropic and anisotropic mdia 图7. 各向同性与各向异性介质大地电磁剖面比较 e Copyright © 2012 Hanspub 249  有限元法求解各向异性介质中大地电磁场边值问题研究 Figure 8. Gradient and high resistivity medium 图8. 倾斜高阻体 (a) (b) 图9. 视电阻率拟断面图 的高阻异常体，其与围岩的分界线清晰明了； 示，电阻率变化规律如图11 所示。 从图 用矩形单元再细分为四个三角形的剖分 方式 的模型，采用有限单元法对各向同性与各向异性介质 电磁测深二维正演，同时对比TE 波和 TM 常幅度约为 24倍，而各向异 性体 Figure 9. Pseudo-section of apparent resistivity 显另一 进行大地 方面可以得出，TE 波较准确 地确定了 高阻 体的纵 向 分布，而横向分辨率不足； TM 波具有很好的横向分辨 率，但却造成了高祖体在纵向上的畸变。 5.4. 高低阻模型 模型如图 10所 中可以清晰的看到低阻异常体和高异常体阻体， 它们与围岩的区别明显，同时可以看到 TE 波具有较 好的纵向分辨率，TM 波具有较好的横向分辨率。 6. 结论 文中采 ，在解线性方程组之前，消去矩形单元中心节点， 这样就达到既不增加网格节点，又可以模拟倾斜物性 界面的目的。最后文中通过引用参考文献[14]和[19] 波，得出两点结论。 1) 同一频率下，在各向同性体电阻率与各向异性 体平均电阻率相同的情况下，二者的异常幅度却相差 很大，各向同性体的异 的异常幅度约为 85倍。因此在大地电磁测深工 作中，介质的各向异性是必须考虑的，而且地球内部 结构与物质组成的各向异性现象普遍存在，各向同性 Figure 10. Gradient and high resistivity medium 图10. 倾斜高阻体 Copyright © 2012 Hanspub 250  有限元法求解各向异性介质中大地电磁场边值问题研究 TE -3-2-10 1 2 3 x/k m -3 -2 -1 0 1 2 3 lg(f) / Hz (a) TM -3-2-101 2 3 x/k m -3 -2 -1 0 1 2 3 lg(f) / Hz (b) Figure 11. Pseudo-section of apparent resistivity 图11. 视电阻率拟断面图 结构只是比较理想的地电模型，只有深入地研究大 电磁场在各向异性介质中的特征，才能得到更接近 际情况的解释; 2) TE波具有很好的纵向分辨率，但横向分辨率 不足，TM 波易于造成异常场在纵向上的畸变，但较 好地确定异常场在横向上的分布，因此仅仅依据某 种类型 此， art´ı, P. Queralt, J. Ledo, C. Farquharson. Dimensionality int of electrical anisotrop y in magnetotelluric responses. Phy- sics of the Earth and Planetary Interiors, 2010, 182(3-4): 139- 151. totelluric impedances and para- ropic layered media. Computers & Geosciences, 2002 28(8): 939-950. M. Hoversten, K. Key and J. S. Chen. Resolution of ric effect of dipping ani- us and anisotropic media. Geophysics, 1975, 40(6): 2. 社, 2011. 06, 28(2): 104-107. 地[1 实[17 一 的电磁波得出的结论，可能是错误的解释，因 在实际工作中应该注意两者的结合，避免望“形” 生意。 参考文献 (References) [1] A. M impr [2] J. Pek, F. A. M. Santos. Magne metric sensitivities for1-D anisot [3] A. M. Osella, P. Martinelli. Magnetotelluric response of anisotropic 2-D structures. Geophysical Journal International, 1993, 115(3): 819-828. [4] V. Brown, reservoir scale electrical anisotropy from marine CSEM data. Geophysics, 2012, 77(2): 147-158. [5] I. K. Reddy, Rankin D. Magnetotellu sotropies. Geophys Prospect, 1971, 19(1): 84-97. [6] I. K. Reddy, D. Rankin. Magnetotelluric Response of Laterally in- homogeneo 1035-1045. [7] 赵生凯, 徐世浙. 有限单元法计算良导嵌入体二维模型的大 地电磁测深曲线 物化探电子计[J]. 算技术, 1983, 1: 14-21. [8] 徐世浙, 赵生凯. 二维各向异性地电断面大地电磁场的有限 元解法[J]. 地震学报, 1985, 7(1): 80-90. [9] 陈乐寿, 王光锷. 大地电磁测深法[M]. 北京: 地质出版社, 1990. [10] 林长佑, 武玉霞, 杨长福等. 水平层状对称各向异性介质的大 地电磁资料反演[J]. 地球物理学报, 1996, 39(增刊): 326-33 [11] 杨长福. 大地电磁二维对称各向异性介质的有限元数值模拟 [J]. 西北地震学报, 1997, 19(2): 27-33. [12] J. Pek, T. Verner. Finite-difference modeling of magnetotelluric fields in two-dimensional anisotropic media. Geophysical Jour- nal International, 1997, 128(1): 505-521. [13] 刘云, 王绪本. 大地电磁二维自适应地形有限元正演模拟[J]. 地震地质, 2010, 32(3): 382-391. [14] 徐世浙. 地球物理中的有限单元法[M]. 北京: 科学出版社, 1994. [15] 熊彬, 罗延钟. 电导率分块均匀的瞬变电磁 2.5维有限元模拟 [J]. 地球物理学报, 2006, 49(2): 590-597. 6] 罗延钟, 张桂清. 电子计算机在电法勘探中的应用[M]. 武汉: 武汉地质学院出版社, 1987. ] 吕玉增, 熊彬, 薛霆虓. 地球物理数据处理基础[M]. 北京: 地质出版 [18] 万汉平. 大地电磁测深的 TE和TM 极化模式对比研究[D]. 成 都理工大学, 2010. [19] 徐凯军, 李桐, 张辉, 李建平. 利用积分方程法的大地电磁三 维正演[J]. 西北地震学报, 20 Copyright © 2012 Hanspub 251 |