大地电磁各向异性二维模拟及实例分析_胡祥云

第５６卷第１２期２０１３年１２月地　球　物　理　学　报ＣＨＩＮＥＳＥ　ＪＯＵＲＮＡＬ　ＯＦ　ＧＥＯＰＨＹＳＩＣＳＶｏｌ．５６，Ｎｏ．１２Ｄｅｃ．，２０１３胡祥云，霍光谱，高锐等．大地电磁各向异性二维模拟及实例分析．地球物理学报，２０１３，５６（１２）：４２６８－４２７７，ｄｏｉ：１０．６０３８／ｃｊｇ２０１３１２２９．Ｈｕ　Ｘ　Ｙ，Ｈｕｏ　Ｇ　Ｐ，Ｇａｏ　Ｒ，ｅｔ　ａｌ．Ｔｈｅ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ａｎｉｓｏｔｒｏｐｉｃ　ｔｗｏ－ｄｉｍｅｎｓｉｏｎａｌ　ｓｉｍｕｌａｔｉｏｎ　ａｎｄ　ｃａｓｅ　ａｎａｌｙｓｉｓ．Ｃｈｉｎｅｓｅ　Ｊ．Ｇｅｏｐｈｙｓ．（ｉｎ　Ｃｈｉｎｅｓｅ），２０１３，５６（１２）：４２６８－４２７７，ｄｏｉ：１０．６０３８／ｃｊｇ２０１３１２２９．大地电磁各向异性二维模拟及实例分析胡祥云１，霍光谱１，３＊，高　锐２，王海燕２，黄一凡１，张云霞１，左博新１，蔡建超１１中国地质大学地球物理与空间信息学院，武汉　４３００７４２中国地质科学院地质研究所岩石圈研究中心，北京　１０００３７３河南省地质调查院，郑州　４５００００摘　要　经过半个多世纪的发展，国内外利用大地电磁法研究地球内部电性结构取得了令人瞩目的成就，这些研究成果多数是基于电性各向同性理论．然而地球内部普遍存在电性各向异性现象，地壳和上地幔中存在的电性各向异性是地电模型、地下结构和构造模型间一个重要的联系因素．本文首先由麦克斯韦方程出发，引入张量电导率，根据二维电性各向异性结构的特点，得到一组关于平行走向的电场分量Ｅｘ和磁场分量Ｈｘ的偏微分方程．使用有限差分法求解偏微分方程，求出Ｅｘ和Ｈｘ的近似解，并以此求得其它场分量；随后，通过对普通及特殊的二维电性各向异性结构做正演模拟，研究其对观测大地电磁场的影响，从而认识在普通及某种特定地质条件下的电磁传播特性，为其后对大地电磁实测资料的处理解释奠定理论基础；最后，以本文的研究成果为基础，将电性各向异性理论引入对实测大地电磁资料的处理解释中，通过对新疆某地的大地电磁资料做二维正演拟合解释，说明了电性各向异性现象的普遍存在，也验证了理论的正确性及算法的实用性，为今后分析解释大地电磁资料中的电性各向异性现象提供理论依据和技术指导，并开拓了对大地电磁实测资料处理的思路和方法．关键词　大地电磁，有限差分，电性各向异性，张量电导率ｄｏｉ：１０．６０３８／ｃｊｇ２０１３１２２９中图分类号　Ｐ６３１收稿日期２０１２－１２－１９，２０１３－０７－０８收修定稿基金项目　国家深部专项（ＳｉｎｏＰｒｏｂｅ－０２－０１），国家自然科学基金（批准号：４１２７４０７７）和湖北省自然科学基金（ＮＯ．２０１１ＣＤＡ１２３）联合资助．作者简介　胡祥云，男，１９６６年生，教授，博士生导师，主要从事电磁法理论应用教学与研究工作．Ｅ－ｍａｉｌ：ｘｙｈｕ＠ｃｕｇ．ｅｄｕ．ｃｎ＊通讯作者　霍光谱，男，１９８３年生，博士，主要从事电磁法正反演研究．Ｅ－ｍａｉｌ：ｈｕｏｇｕａｎｇｐｕ１２１８＠１６３．ｃｏｍＴｈｅ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ａｎｉｓｏｔｒｏｐｉｃ　ｔｗｏ－ｄｉｍｅｎｓｉｏｎａｌ　ｓｉｍｕｌａｔｉｏｎ　ａｎｄ　ｃａｓｅ　ａｎａｌｙｓｉｓＨＵ　Ｘｉａｎｇ－Ｙｕｎ１，ＨＵＯ　Ｇｕａｎｇ－Ｐｕ１，３＊，ＧＡＯ　Ｒｕｉ　２，ＷＡＮＧ　Ｈａｉ－Ｙａｎ２，ＨＵＡＮＧ　Ｙｉ－Ｆａｎ１，ＺＨＡＮＧ　Ｙｕｎ－Ｘｉａ１，ＺＵＯ　Ｂｏ－Ｘｉｎ１，ＣＡＩ　Ｊｉａｎ－Ｃｈａｏ１１　Ｉｎｓｔｉｔｕｔｅ　ｏｆ　Ｇｅｏｐｈｙｓｉｃｓ　ａｎｄ　Ｇｅｏｍａｔｉｃｓ，Ｃｈｉｎａ　Ｕｎｉｖｅｒｓｉｔｙ　ｏｆ　Ｇｅｏｓｃｉｅｎｃｅｓ，Ｗｕｈａｎ　４３００７４，Ｃｈｉｎａ２　Ｌｉｔｈｏｓｐｈｅｒｅ　Ｒｅｓｅａｒｃｈ　Ｃｅｎｔｅｒ，Ｉｎｓｔｉｔｕｔｅ　ｏｆ　Ｇｅｏｌｏｇｙ，ＣＡＧＳ，Ｂｅｉｊｉｎｇ１０００３７，Ｃｈｉｎａ３　Ｈｅｎａｎ　Ｉｎｓｔｉｔｕｔｅ　ｏｆ　Ｇｅｏｌｏｇｉｃａｌ　Ｓｕｒｖｅｙ，Ｚｈｅｎｇｚｈｏｕ４５００００，ＣｈｉｎａＡｂｓｔｒａｃｔ　Ａｆｔｅｒ　ｍｏｒｅ　ｔｈａｎ　ｈａｌｆ　ｃｅｎｔｕｒｙ　ｄｅｖｅｌｏｐｍｅｎｔ，ｒｅｍａｒｋａｂｌｅ　ａｃｈｉｅｖｅｍｅｎｔｓ　ｈａｖｅ　ｂｅｅｎ　ｍａｄｅｏｎ　ｔｈｅ　ｓｔｕｄｙ　ｏｆ　ｅａｒｔｈ′ｓ　ｉｎｔｅｒｉｏｒ　ｅｌｅｃｔｒｉｃａｌ　ｓｔｒｕｃｔｕｒｅ　ｂｙ　ｕｓｉｎｇ　ｔｈｅ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｍｅｔｈｏｄ　ａｔ　ｈｏｍｅａｎｄ　ａｂｒｏａｄ．Ｔｈｅ　ｍａｊｏｒｉｔｙ　ｏｆ　ｔｈｅｓｅ　ｒｅｓｅａｒｃｈ　ｒｅｓｕｌｔｓ　ａｒｅ　ｂａｓｅｄ　ｏｎ　ｔｈｅ　ｅｌｅｃｔｒｉｃａｌ　ｉｓｏｔｒｏｐｉｃ　ｔｈｅｏｒｙ．Ｈｏｗｅｖｅｒ，ｔｈｅ　ｅｌｅｃｔｒｉｃａｌ　ａｎｉｓｏｔｒｏｐｙ　ｐｈｅｎｏｍｅｎｏｎ　ｉｓ　ｐｒｅｖａｌｅｎｔ　ｉｎ　ｔｈｅ　ｅａｒｔｈ　ｉｎｔｅｒｉｏｒ，ａｎｄ　ｔｈｅｐｒｅｓｅｎｃｅ　ｏｆ　ｅｌｅｃｔｒｉｃａｌ　ａｎｉｓｏｔｒｏｐｙ　ｉｎ　ｔｈｅ　ｃｒｕｓｔ　ａｎｄ　ｕｐｐｅｒ　ｍａｎｔｌｅ　ｉｓ　ａｎ　ｉｍｐｏｒｔａｎｔ　ｌｉｎｋ　ｆａｃｔｏｒ　ａｍｏｎｇｇｅｏｅｌｅｃｔｒｉｃ　ｍｏｄｅｌ，ｕｎｄｅｒｇｒｏｕｎｄ　ｓｔｒｕｃｔｕｒｅ　ａｎｄ　ｔｅｃｔｏｎｉｃ　ｍｏｄｅｌ．Ｉｎ　ｔｈｉｓ　ｐａｐｅｒ，ｔｈｅ　ｔｅｎｓｏｒｃｏｎｄｕｃｔｉｖｉｔｙ　ｈａｓ　ｂｅｅｎ　ｆｉｒｓｔ　ｉｎｔｒｏｄｕｃｅｄ　ｂｙ　ｓｔａｒｔｉｎｇ　ｆｒｏｍ　ｔｈｅ　Ｍａｘｗｅｌｌ　ｅｑｕａｔｉｏｎｓ，ｔｈｅｎ　ｇｅｔｔｉｎｇ　ａ　ｓｅｔ　１２期胡祥云等：大地电磁各向异性二维模拟及实例分析ｏｆ　ｐａｒｔｉａｌ　ｄｉｆｆｅｒｅｎｔｉａｌ　ｅｑｕａｔｉｏｎｓ　ｗｉｔｈ　ｒｅｓｐｅｃｔ　ｔｏ　Ｅｘａｎｄ　Ｈｘａｃｃｏｒｄｉｎｇ　ｔｏ　ｔｈｅ　ｃｈａｒａｃｔｅｒｉｓｔｉｃ　ｏｆ　ｔｗｏ－ｄｉｍｅｎｓｉｏｎ　ｅｌｅｃｔｒｉｃａｌ　ａｎｉｓｏｔｒｏｐｉｃ　ｓｔｒｕｃｔｕｒｅ．Ｂｙ　ｕｓｉｎｇ　ｔｈｅ　ｐａｒｔｉａｌ　ｄｉｆｆｅｒｅｎｔｉａｌ　ｅｑｕａｔｉｏｎｓ，ｔｈｅａｐｐｒｏｐｒｉａｔｅ　ｓｏｌｕｔｉｏｎ　ｏｆ　Ｅｘａｎｄ　Ｈｘｈａｓ　ｂｅｅｎ　ｓｏｌｖｅｄ，ａｎｄ　ｗｏｒｋｉｎｇ　ｏｕｔ　ｏｔｈｅｒ　ｆｉｅｌｄ　ｃｏｍｐｏｎｅｎｔｓ　ｏｎ　ｔｈｅｂａｓｉｓ　ｏｆ　Ｅｘａｎｄ　Ｈｘ．Ｔｈｅｎ，ｉｎ　ｏｒｄｅｒ　ｔｏ　ｒｅｃｏｇｎｉｚｅ　ｔｈｅ　ｃｈａｒａｃｔｅｒｉｓｔｉｃｓ　ｏｆ　ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃ　ｐｒｏｐａｇａｔｉｏｎｕｎｄｅｒ　ｃｏｍｍｏｎ　ａｎｄ　ｃｅｒｔａｉｎ　ｇｅｏｌｏｇｉｃａｌ　ｃｏｎｄｉｔｉｏｎｓ，ｔｈｅ　ｉｎｆｌｕｅｎｃｅ　ｏｎ　ｔｈｅ　ｍｅａｓｕｒｅｄ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃｈａｓ　ｂｅｅｎ　ｒｅｓｅａｒｃｈｅｄ　ｂｙ　ｔｈｅ　ｆｏｒｗａｒｄ　ｓｉｍｕｌａｔｉｏｎ　ｏｎ　ｏｒｄｉｎａｒｙ　ａｎｄ　ｓｐｅｃｉａｌ　ｔｗｏ－ｄｉｍｅｎｓｉｏｎ　ａｎｉｓｏｔｒｏｐｉｃｅｌｅｃｔｒｉｃａｌ　ｓｔｒｕｃｔｕｒｅ，ｗｈｉｃｈ　ｈａｓ　ｌａｉｄ　ｔｈｅｏｒｅｔｉｃａｌ　ｆｏｕｎｄａｔｉｏｎ　ｆｏｒ　ｔｈｅ　ｓｕｂｓｅｑｕｅｎｔ　ｐｒｏｃｅｓｓｉｎｇ　ａｎｄｉｎｔｅｒｐｒｅｔａｔｉｏｎ　ｏｆ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｍｅａｓｕｒｅｄ　ｄａｔａ．Ｆｉｎａｌｌｙ，ｔｈｅ　ｅｌｅｃｔｒｉｃａｌ　ａｎｉｓｏｔｒｏｐｙ　ｔｈｅｏｒｙ　ｉｓｉｎｔｒｏｄｕｃｅｄ　ｉｎ　ｔｈｅ　ｐｒｏｃｅｓｓｉｎｇ　ａｎｄ　ｉｎｔｅｒｐｒｅｔａｔｉｏｎ　ｏｆ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｍｅａｓｕｒｅｄ　ｄａｔａ　ｂａｓｅｄ　ｏｎ　ｔｈｅｒｅｓｅａｒｃｈ　ｒｅｓｕｌｔｓ　ｉｎ　ｔｈｅ　ｐａｐｅｒ，ａｎｄ　ｔｈｅ　ｗｉｄｅｓｐｒｅａｄ　ｅｘｉｓｔｅｎｃｅ　ｏｆ　ｅｌｅｃｔｒｉｃａｌ　ａｎｉｓｏｔｒｏｐｙ　ｐｈｅｎｏｍｅｎｏｎｈａｓ　ｂｅｅｎ　ｓｔａｔｅｄ　ａｎｄ　ｅｘａｃｔｎｅｓｓ　ｏｆ　ｔｈｉｓ　ｔｈｅｏｒｙ　ａｎｄ　ｐｒａｃｔｉｃａｂｉｌｉｔｙ　ｏｆ　ａｌｇｏｒｉｔｈｍ　ｈａｖｅ　ｂｅｅｎ　ｖｅｒｉｆｉｅｄ　ｂｙｍｅａｎｓ　ｏｆ　ｔｈｅ　ｔｗｏ－ｄｉｍｅｎｓｉｏｎ　ｆｏｒｗａｒｄ　ｆｉｔｔｉｎｇ　ｉｎｔｅｒｐｒｅｔａｔｉｏｎ　ｏｎ　ｔｈｅ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｄａｔａ　ｉｎｓｏｍｅｐｌａｃｅ　ｉｎ　Ｘｉｎｊｉａｎｇ，ｗｈｉｃｈ　ｃａｎ　ｐｒｏｖｉｄｅ　ｔｈｅｏｒｅｔｉｃａｌ　ｂａｓｉｓ　ａｎｄ　ｔｅｃｈｎｉｃａｌ　ｇｕｉｄａｎｃｅ　ｆｏｒ　ａｎａｌｙｚｉｎｇａｎｄ　ｅｘｐｌａｉｎｉｎｇ　ｔｈｅ　ｅｌｅｃｔｒｉｃａｌ　ａｎｉｓｏｔｒｏｐｙ　ｐｈｅｎｏｍｅｎｏｎ　ｉｎ　ｔｈｅ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｄａｔａ，ａｎｄ　ｍｏｓｔｉｍｐｏｒｔａｎｔｌｙ，ｉｔ　ｏｐｅｎｓ　ｕｐ　ｔｈｅ　ｉｄｅａ　ａｎｄ　ｓｏｌｕｔｉｏｎ　ｉｎ　ｔｈｅ　ｐｒｏｃｅｓｓｉｎｇ　ｏｆ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｍｅａｓｕｒｅｄ　ｄａｔａ．Ｋｅｙｗｏｒｄｓ　Ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ，Ｆｉｎｉｔｅ　ｄｉｆｆｅｒｅｎｃｅ，Ｅｌｅｃｔｒｉｃａｌ　ａｎｉｓｏｔｒｏｐｙ，Ｔｅｎｓｏｒ　ｃｏｎｄｕｃｔｉｖｉｔｙ１　引　言大地电磁法是研究地壳及上地幔电性结构的主要手段，近年来，我国学者取得了显著的进展［１－４］．Ｗａｎａｍａｋｅｒ［５］在２００５年详细阐述了电性各向异性在地电模型和地球动力学解释中的重要作用，同时说明了电性各向异性现象在地球内部是普遍存在的．二维大地电磁电性各向异性正演模拟是通过解偏微分方程或积分方程的数值计算方法实现的，被广泛使用的两种方法是：有限元法和有限差分法［６］．早在１９７５年，Ｒｅｄｄｙ［７］等采用有限元法对偏微分方程做了数值计算，较早的开展了二维电性各向异性问题的理论研究工作．虽然在其后的二十多年间（直到２０００年之后才有学者将其理论改进后用于大地电磁实测资料处理）未能广泛运用在大地电磁实测资料处理中，但其所做研究在理论意义上的影响却尤为深远．徐世浙等［８］用有限元法计算了二维电性各向异性地电断面的大地电磁场，是我国较早见于公开发表的理论研究文章．李予国［９］公布了用于计算二维电性各向异性结构感应电磁场的有限元公式．该公式采用有限元法求解表示平行于走向的场分量Ｅｘ和Ｈｘ的一组耦合方程组得到了有限元的线性方程，使用预条件共轭梯度法求得这个线性方程的结果；通过样条插值对Ｅｘ和Ｈｘ进行数值微分求得垂直于走向的场分量Ｅｙ和Ｈｙ．随后，与Ｐｅｋ等［１０］数值计算的结果做了对比，证明了理论的正确性．李予国等［１１］采用自适应网格剖分法对之前的有限元法正演计算的二维电性各向异性结构的大地电磁响应做了一次成功改进，并取得了很好的效果．Ｐｅｋ等［１０］在Ｒｅｅｄｙ等［７］研究基础上，将前人（Ｈａａｋ［１２］，Ｂｒｅｗｉｔｔ－Ｔａｙｌｏｒ　＆Ｗｅａｖｅｒ［１３］，ˇＣｅｒｖ　＆Ｐｒａｕｓ［１４］等）在二维电性各向同性介质中使用的有限差分法引入求解二维电性各向异性问题中，用有限差分法对偏微分方程做了数值模拟计算，是正演模拟计算二维电性各向异性介质中大地电磁响应问题影响较为深远的一种求解方法．２　大地电磁电性各向异性二维正演理论本文采用有限差分法对大地电磁二维电性各向异性正演问题进行研究，参照了Ｐｅｋ等［１０］的研究方法．假设在笛卡尔坐标系中，有一个二维电性结构，它的构造走向平行于ｘ轴，ｙ轴垂直于ｘ轴，ｚ轴垂直于ｘｙ平面，且正向向下．模型由一个有限区域组成，在其中的电性各向异性区域是二维有限结构．假设地表为水平，对应ｚ＝０．在地表上空（ｚ＜０）为绝缘空气层．来自ｚ→－∞的平面波垂直向下传播，ω＝２π／Ｔ，Ｔ为周期，单位为秒（ｓ）．谐变电磁场的麦克斯韦方程组为：Δ×Ｅ＝ｉωμ０Ｈ，（１）Δ×Ｈ＝σＥ．（２） 上式为描述电磁场分布的微分方程．ｅ－ｉωｔ为时９６２４地球物理学报（Ｃｈｉｎｅｓｅ　Ｊ．Ｇｅｏｐｈｙｓ．）５６卷　谐因子，σ为张量电导率．ｘ轴平行于构造走向，由２－Ｄ性质 ／ｘ＝０，公式（１）—（２）分解为： Ｅｚ ｙ－ Ｅｙ ｚ＝ｉωμ０Ｈｘ，（３） Ｅｘ ｚ＝ｉωμ０Ｈｙ，（４）－ Ｅｘ ｙ＝ｉωμ０Ｈｚ，（５） Ｈｚ ｙ－ Ｈｙ ｚ＝σｘｘＥｘ＋σｘｙＥｙ＋σｘｚＥｚ，（６） Ｈｘ ｚ＝σｙｘＥｘ＋σｙｙＥｙ＋σｙｚＥｚ，（７）－ Ｈｘ ｙ＝σｚｘＥｘ＋σｚｙＥｙ＋σｚｚＥｚ．（８） 在公式（３）—（８）中，如果Ｅｘ和Ｈｘ已知，那么Ｅｙ，Ｅｚ，Ｈｙ和Ｈｚ均可求得；将公式简化为只含Ｅｘ和Ｈｘ的一组二阶偏微分方程．电场沿走向时（即ＴＥ模式）： ２　Ｅｘ ｙ２＋ ２　Ｅｘ ｚ２＋ｉωμ０（σｘｘ＋Ａσｚｘ＋Ｂσｘｙ）Ｅｘ ＋ｉωμ０Ａ Ｈｘ ｙ＋ｉωμ０Ｂ Ｈｘｚ＝０，（９）磁场沿走向时（即ＴＭ模式）： ｙσｙｙＤ Ｈｘ （）ｙ＋ ｚσｚｚＤ Ｈｘ （）ｚ＋ｙσｙｚＤ Ｈｘ （）ｚ＋ｉωμ０Ｈｘ＋ ｚσｙｚＤ Ｈｘ （）ｙ－ （ＢＥｘ） ｙ＋ （ＡＥｘ） ｚ＝０，（１０）公式（９）—（１０）中：Ａ＝（σｘｙσｙｚ－σｘｚσｙｙ）／Ｄ，Ｂ＝（σｘｙσｙｚ－σｘｚσｚｚ）／Ｄ，Ｄ＝σｙｙσｚｚ－σ２ｙｚ．用有限差分法求解公式（９）—（１０）：首先将二维结构投影到一个数值网格（是一个有限网格区域），在电性各向异性区域，再分为一个矩形单元的电性均匀体系．网格通常是不规则的，首先应该符合模型的几何形式，同时满足数值运算的一般规则，即设计的数值网格应该满足数值模拟的需要．使用ˇＣｅｒｖ和Ｐｒａｕｓ［１４］的方法，将积分插值法用于导出各个网格节点的有限差分方程．由图１所示，黑色圆实点代表节点（ｊ，ｋ），阴影区域代表该节点四周的矩形积分单元Ｇｊｋ，通过该积分单元对方程（９）和（１０）求积分．矩形积分单元Ｇｊｋ表示为：Ｇｊｋ＝［ｙ－ｊ，ｙ＋ｊ］×［ｚ－ｋ，ｚ＋ｋ］＝ｙｊ－１２Δｙｊ，ｙｊ＋１２Δｙｊ＋［］１×ｚｋ－１２Δｚｋ，ｚｋ＋１２Δｚｋ＋［］１．图１　节点（ｊ，ｋ）周围的有限差分网格Ｆｉｇ．１　Ｆｉｎｉｔｅ　ｄｉｆｆｅｒｅｎｃｅ　ｍｅｓｈｅｓ　ａｒｏｕｎｄ　ｔｈｅ　ｎｏｄｅ（ｊ，ｋ）求解这些积分，需用各自剖分网格四周的四个节点做近似计算．公式（６）对应的积分形式为：∮Γ（Ｇｊｋ）Ｈ·ｄｇ＝∫Ｇｊｋ（σＥ）·ｄＧ，（１１）公式（１１）中Γ（Ｇｊｋ）是积分单元Ｇｊｋ的定向边界，ｄｇ是沿着此边界的积分路径的元素，ｄＧ＝ｄＧｅｘ，ｅｘ是在ｘ方向的单位矢量，ｄＧ是Ｇｊｋ的区域元素，表示为： Ｅｘ（）ｚＫＬＫＬ－Ｅｘ （）ｙＬＮＬＮ＋Ｅｘ（）ｚＮＭＮＭ－ Ｅｘ（）ｙＭＫＭＫ－ｉωμ０Ｓ（Ｇｊｋ）（ｊｘ）Ｇｊｋ＝０，（１２）公式（１２）中 Ｅｘ（）ｚＫＬ是指函数 Ｅｘｚ穿过积分单元边界ＫＬ的值，相应的 Ｅｘ （）ｙＬＮ，Ｅｘ（）ｚＮＭ， Ｅｘ（）ｙＭＫ意义可类推．ＫＬ是积分单元点Ｋ、Ｌ之间的距离（从ｘ正向看，沿逆时针方向为正），Ｓ（Ｇｊｋ）是积分单元的面积，ｊｘ是电流密度的ｘ分量，ｊｘ＝（σＥ）ｘ＝σｘｘＥｘ＋σｘｙＥｙ＋σｘｚＥｚ，（ｊｘ）Ｇｊｋ是ｊｘ穿过积分单元Ｇｊｋ的值．基于公式（１２）的近似计算，得到了一个在网格节点（ｊ，ｋ）对于ＴＥ模式的近似差分方程，这个近似差分方程代表了公式（９）的有限差分形式：　∑ｊ＋１ｐ＝ｊ－１∑ｋ＋１ｑ＝ｋ－１ＣＥＥｊｋ（ｐ，ｑ）Ｅｘ（ｐ，ｑ）＋∑ｊ＋１ｐ＝ｊ－１∑ｋ＋１ｑ＝ｋ－１ＣＥＨｊｋ（ｐ，ｑ） ×Ｈｘ（ｐ，ｑ）＝０，（１３）公式（１３）中ｐ可以是ｊ－１，ｊ和ｊ＋１，ｑ可以是ｋ－１，ｋ和ｋ＋１．Ｃ表示各个剖分节点上是否有Ｅｘ和Ｈｘ分量，其上标的第一个Ｅ表示为ＴＥ模式，第二个Ｅ和Ｈ分别为电场分量和磁场分量；且当（ｐ，ｑ）０７２４　１２期胡祥云等：大地电磁各向异性二维模拟及实例分析∈｛（ｊ－１，ｋ－１），（ｊ－１，ｋ＋１），（ｊ＋１，ｋ－１），（ｊ＋１，ｋ＋１）｝时，有：ＣＥＥｊｋ（ｐ，ｑ）＝ＣＥＨｊｋ（ｐ，ｑ）＝０，对应的ＴＭ模式积分形式，由公式（３）在积分单元Ｇｊｋ上的积分表示为：∮Γ（Ｇｊｋ）Ｅ·ｄｇ＝∫Ｇｊｋ（σＥ）·ｄＧ，（１４）同样可以获得在网格节点（ｊ，ｋ）对于Ｈ极化下的线性差分方程：　∑ｊ＋１ｐ＝ｊ－１∑ｋ＋１ｑ＝ｋ－１ＣＨＥｊｋ（ｐ，ｑ）Ｅｘ（ｐ，ｑ）＋∑ｊ＋１ｐ＝ｊ－１∑ｋ＋１ｑ＝ｋ－１ＣＨＨｊｋ（ｐ，ｑ） ×Ｈｘ（ｐ，ｑ）＝０，（１５）与ＴＥ模式的线性差分公式（１３）相比，在公式（１５）中，当（ｐ，ｑ）∈｛（ｊ－１，ｋ－１），（ｊ－１，ｋ＋１），（ｊ＋１，ｋ－１），（ｊ＋１，ｋ＋１）｝时，只有：ＣＨＥｊｋ（ｐ，ｑ）＝０ 在这些节点，相应的磁场系数ＣＨＨｊｋ（ｐ，ｑ）一般不为零，它代表了电性各向异性的影响．系数矩阵Ｃαβｊｋ（ｐ，ｑ）的元素并不是每一个都有数值，也就是说在任意情况下，都不用计算全部９个节点的系数元素（１８个值）．对于ＴＥ模式，有５个节点（１０个值）；对于ＴＭ模式，有全部９个节点，但只有１４个非零系数元素．在全电性各向异性情况下，系数参数的个数会达到最大值．但随着网格单元包含的电性各向异性的减少，非零的系数参数也会逐渐减少．用有限差分法对公式（９）和公式（１０）做数值近似后，须将线性代数公式（１３）和公式（１５）安排在一个系统，以便做进一步的处理．因为包含了Ｅｘ和Ｈｘ两个变量，且这两个变量不在相同的网格节点，ＴＥ模式的方程需要在全部网格节点做近似计算，而ＴＭ模式的方程只需在导电底层内部做近似计算．由图２可知，系数矩阵是一个对称、带宽限定矩阵，可以使用高斯消元对其进行特定的修改，从而达到求解有限差分线性代数方程组的目的．为了存储矩阵的上半频带，需要对高斯消元公式做一次改进来实现，复数（ＮＳＴＯＲ）需被放置在内存中，ＮＳＴＯＲ为：ＮＳＴＯＲ＝（Ｎ－１）（２　ＭＥ＋ＭＡ－２烐烏烑）ｎｕｍｂｅｒ　ｏｆ　ｅｑｕａｔｉｏｎｓ（２　ＭＥ＋ＭＡ＋１烐烏烑）ｂａｎｄ－ｗｉｄｔｈ，（１６）公式（１６）中Ｎ是水平网格步长的数量，ＭＥ、ＭＡ分别是导电地层和空气中垂向网格步长的数量．３　二维电性各向异性模型计算与分析大地电磁法被普遍应用于地壳与上地幔大尺度构造研究中，在一般情况下，地形的起伏相对于测量剖面的长度较小，所以用水平层状二维电性各向异性结构的正演模拟，可满足多数地质条件下电磁传图２　变量在网格节点中的交错排列方式符号表示了有限差分线性代数方程组系数矩阵在各个节点的存在情况．（ａ）中，空心圆节点表示只有电场分量，实心圆节点表示同时有电场和磁场分量；（ｂ）中，实心圆代表电场分量，实心方块代表磁场分量，空心符号表示在矩阵模式中由电性各向异性存在造成的系数．Ｆｉｇ．２　Ｔｈｅ　ｓｔａｇｇｅｒｅｄ　ｍａｎｎｅｒ　ｏｆ　ｖａｒｉａｂｌｅｓ　ｉｎ　ｇｒｉｄ　ｎｏｄｅｓＴｈｅ　ｓｙｍｂｏｌｓ　ｉｎｄｉｃａｔｅ　ｔｈｅ　ｐｒｅｓｅｎｃｅ　ｏｆ　ｅａｃｈ　ｎｏｄｅ　ｏｆ　ｃｏｅｆｆｉｃｉｅｎｔ　ｍａｔｒｉｘ　ｏｆ　ｆｉｎｉｔｅ　ｄｉｆｆｅｒｅｎｃｅ　ｌｉｎｅａｒ　ａｌｇｅｂｒａｉｃ　ｅｑｕａｔｉｏｎｓ．（ａ）ｈｏｌｌｏｗ　ｃｉｒｃｌｅ　ｎｏｄｅｓｉｎｄｉｃａｔｅ　ｏｎｌｙ　ｔｈｅ　ｅｌｅｃｔｒｉｃ　ｆｉｅｌｄ　ｃｏｍｐｏｎｅｎｔ，ｓｏｌｉｄ　ｃｉｒｃｌｅ　ｎｏｄｅｓ　ｉｎｄｉｃａｔｅ　ｂｏｔｈ　ｔｈｅ　ｅｌｅｃｔｒｉｃ　ａｎｄ　ｍａｇｎｅｔｉｃ　ｆｉｅｌｄ　ｃｏｍｐｏｎｅｎｔｓ；（ｂ）ｓｏｌｉｄ　ｃｉｒｃｌｅｓ　ａｎｄｓｏｌｉｄ　ｓｑｕａｒｅｓ　ｉｎｄｉｃａｔｅ　ｅｌｅｃｔｒｉｃ　ａｎｄ　ｍａｇｎｅｔｉｃ　ｆｉｅｌｄ　ｃｏｍｐｏｎｅｎｔｓ，ｒｅｓｐｅｃｔｉｖｅｌｙ，ｈｏｌｌｏｗ　ｓｙｍｂｏｌｓ　ｉｎｄｉｃａｔｅ　ｔｈｅ　ｃｏｅｆｆｉｃｉｅｎｔｓ　ｄｕｅ　ｔｏ　ｔｈｅ　ｐｒｅｓｅｎｃｅ　ｏｆｅｌｅｃｔｒｉｃｉｔｙ　ａｎｉｓｏｔｒｏｐｙ　ｉｎ　ｔｈｅ　ｍａｔｒｉｘ　ｍｏｄｅ．１７２４地球物理学报（Ｃｈｉｎｅｓｅ　Ｊ．Ｇｅｏｐｈｙｓ．）５６卷　播特性的研究目的．在两层电性各向同性介质中放置一个二维电性各向异性介质，对地电模型做正演计算，分析研究介质中大地电磁场的形态特征和分布情况．第一层称为介质１，为各向同性，在第二层各向同性介质（称为介质３）中放置各向异性介质２．模型参数为：介质１：ρ１＝３００Ωｍ；介质２：ρ１２＝１０Ωｍ，ρ２２＝１００Ωｍ，ρ３２＝１００Ωｍ，α＝３０°，其中α表示测量轴与电性主轴的夹角；介质３：ρ３＝１０００Ωｍ．介质１深度范围是０ｋｍ到－４ｋｍ；介质２的大小为２０×４ｋｍ２，在水平方向上的位置为２４ｋｍ至４４ｋｍ之间、深度方向上的位置为－６ｋｍ至－１０ｋｍ之间．如图３所示：图３　层状电性各向同性介质中电性各向异性介质地电模型Ｆｉｇ．３　Ｅｌｅｃｔｒｉｃａｌ　ａｎｉｓｏｔｒｏｐｉｃ　ｍｅｄｉｕｍ　ｇｅｏｅｌｅｃｔｒｉｃｍｏｄｅｌ　ｉｎ　ｔｈｅ　ｌａｙｅｒｅｄ　ｉｓｏｔｒｏｐｉｃ　ｍｅｄｉｕｍ 正演计算时将模型剖分为４５×３５的剖分单元（垂向３５个剖分单元中包含１０个空气层的剖分单元），测点数为４５个（即在地表的网格节点处），采用２１个周期点（０．１，０．３，０．５，０．７，１，３，５，７，１０，３０，５０，７０，１００，３００，５００，７００，１０００，１３００，１５００，１７００，２０００，单位为ｓ）．计算结果如图４所示（图中视电阻率拟断面图单位为Ωｍ，相位拟断面图单位为°）：图４展示了正演计算得到的大地电磁响应结果．正演计算的结果能够很好的反映出电性各向异性介质的存在范围和区域：ＴＥ模式的拟断面图清晰地反映了电性各向异性介质在纵向上的顶部埋深和延展范围；而ＴＭ模式的拟断面图反映电性各向异性介质在横向上的延续范围和长度比较准确．图４也很好的反映了地电模型的层状性质，在频率ｌｏｇ（ｆ）＝－１以上部分，视电阻率值表现出了覆盖层的相对低阻特性；在频率ｌｏｇ（ｆ）＝－１以下部分，视电阻率值表现出了背景场值的高阻特性．对图３的地电模型做一个微小的调整，将介质２的参数改为：ρ１２＝１０Ωｍ，ρ２２＝１００Ωｍ，ρ３２＝１００Ωｍ，α＝０°．即：令测量轴和电性主轴相互重合，但介质电阻率值保持不变．通过正演计算，研究两种地电模型中大地电磁场分布的差异，计算结果如图５所示（图中视电阻率拟断面图单位为Ωｍ，相位拟断面图单位为°）：图５为修改后模型正演计算得到的大地电磁响应结果．图５与图４在形态上大体相似，但细节处又有明显不同，主要表现在：ＴＥ模式下图５（ａ）的低阻区域较图４（ａ）有明显扩大，图５（ｃ）的高值区域有明显扩大，而低值区域的延展范围有了明显缩小；ＴＭ模式下图５（ｂ）的低阻区域较图４（ｂ）在水平方向上略有收窄，图５（ｄ）中高值区域的延展范围有了明显缩小；ＴＥ模式下图５（ａ）的视电阻率最小值较图４（ａ）有了明显降低，而ＴＭ模式下图５（ｂ）的视电阻率最小值较图４（ｂ）有了明显升高．这些细微的不同表明了观测角对于大地电磁响应结果的影响和作用：观测角的存在会令两个电性主轴的电阻率值有平均的效果，即具有低值电性主轴表现出来的视电阻率值不会太低，具有高值电性主轴表现出来的视电阻率值不会太高；但当电性主轴为０时，具有低值电性主轴表现出来的视电阻率值会降低，具有高值电性主轴表现出来的视电阻率值会升高．对图３的地电模型再做一个微小的调整，将介质２的参数改为：ρ１２＝１０Ωｍ，ρ２２＝１０００Ωｍ，ρ３２＝１０００Ωｍ，α＝０°．即：不仅令测量轴和电性主轴相互重合，而且使ρ２２、ρ３２的电阻率值等于背景场值．计算结果如图６所示（图中视电阻率拟断面图单位为Ωｍ，相位拟断面图单位为°）．由图６的正演计算结果可知：在ＴＭ模式下，无法识别该模型的二维结构，只表现出了层状介质的电性各向同性性质．原因是当观测角为０时，测量轴与电性主轴方向是重合的，在电阻率值和背景值是相同的情况下，该方向具有层状同性电性特征．在ＴＥ模式下，与图５一样，可以明显的看到一个低阻异常体的存在，并且，低阻异常的范围较初始模型计算得到的结果（图４（ａ））有了非常明显的扩大，特别是在深度方向上异常的范围比水平方向上的异常范围要延展很多．通过对二维电性各向异性介质的正演计算，并对正演计算结果进行研究分析，能够做出以下结论：在电性各向异性主轴同测量轴存在夹角的情况下，ＴＥ、ＴＭ模式都可以有效的分辨出异常区域，并对二维结构在顶部埋深、深度和宽度的延展上都有非常好的反映；观测角的存在会令两个电性主轴的电阻率值有平均的效果，在电性各向异性主轴同测量２７２４　１２期胡祥云等：大地电磁各向异性二维模拟及实例分析图４　层状电性各向同性介质中电性各向异性介质地电模型的大地电磁响应拟断面图（ａ）ＴＥ模式视电阻率拟断面图；（ｂ）ＴＭ模式视电阻率拟断面图；（ｃ）ＴＥ模式阻抗相位拟断面图（ｄ）ＴＭ模式阻抗相位拟断面图．Ｆｉｇ．４　Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｒｅｓｐｏｎｓｅ　ｏｆ　ｅｌｅｃｔｒｉｃａｌ　ａｎｉｓｏｔｒｏｐｉｃｍｅｄｉｕｍ　ｇｅｏｅｌｅｃｔｒｉｃ　ｍｏｄｅｌ　ｉｎ　ｔｈｅ　ｌａｙｅｒｅｄ　ｉｓｏｔｒｏｐｉｃ　ｍｅｄｉｕｍ（ａ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ａｐｐａｒｅｎｔ　ｒｅｓｉｓｔｉｖｉｔｙ　ｉｎ　ＴＥ　ｍｏｄｅ；（ｂ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ａｐｐａｒｅｎｔ　ｒｅｓｉｓｔｉｖｉｔｙ　ｉｎ　ＴＭ　ｍｏｄｅ；（ｃ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ｉｍｐｅｄａｎｃｅ　ｐｈａｓｅ　ｉｎ　ＴＥ　ｍｏｄｅ；（ｄ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ｉｍｐｅｄａｎｃｅ　ｐｈａｓｅ　ｉｎ　ＴＭ　ｍｏｄｅ．图５　层状电性各向同性介质中电性各向异性介质地电模型的大地电磁响应拟断面图（ａ）ＴＥ模式视电阻率拟断面图；（ｂ）ＴＭ模式视电阻率拟断面图；（ｃ）ＴＥ模式阻抗相位拟断面图；（ｄ）ＴＭ模式阻抗相位拟断面图．Ｆｉｇ．５　Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｒｅｓｐｏｎｓｅ　ｏｆ　ｅｌｅｃｔｒｉｃａｌ　ａｎｉｓｏｔｒｏｐｉｃｍｅｄｉｕｍ　ｇｅｏｅｌｅｃｔｒｉｃ　ｍｏｄｅｌ　ｉｎ　ｔｈｅ　ｌａｙｅｒｅｄ　ｉｓｏｔｒｏｐｉｃ　ｍｅｄｉｕｍ（ａ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ａｐｐａｒｅｎｔ　ｒｅｓｉｓｔｉｖｉｔｙ　ｉｎ　ＴＥ　ｍｏｄｅ；（ｂ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ａｐｐａｒｅｎｔ　ｒｅｓｉｓｔｉｖｉｔｙ　ｉｎ　ＴＭ　ｍｏｄｅ；（ｃ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ｉｍｐｅｄａｎｃｅ　ｐｈａｓｅ　ｉｎ　ＴＥ　ｍｏｄｅ；（ｄ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ｉｍｐｅｄａｎｃｅ　ｐｈａｓｅ　ｉｎ　ＴＭ　ｍｏｄｅ．３７２４地球物理学报（Ｃｈｉｎｅｓｅ　Ｊ．Ｇｅｏｐｈｙｓ．）５６卷　轴不存在夹角的情况下，若测量轴其中的一个方向上的电阻率值同背景场值相同，必然会造成对应的极化模式无法分辨出二维电性各向异性结构．４　大地电磁实测资料解释及电性结构分析选取的大地电磁测区位于新疆北部，该测线（图７所示）位于准噶尔盆地的古尔班通古特沙漠边缘，地势平坦，整条测线高程起伏较小，海拔高度在３００～４００ｍ．测点均在沙漠边缘，干扰较小．中心测点（Ｃ）位置为Ｎ４５°００′，Ｅ８６°００′，西一测点（Ｗ１）距中心点约１５ｋｍ，西二测点（Ｗ２）距中心点约４０ｋｍ，东一测点（Ｅ１）距中心点约２０ｋｍ，东二测点（Ｅ２）距中心点约４０ｋｍ．测线各个测点的大地电磁视电阻率和阻抗相位的曲线形态都非常平滑，仅在个别频点上存在飞点，只需做简单的圆滑处理．对此测线的大地电磁数据选取了频率为１１．２～０．０００５Ｈｚ之间的数据进行分析研究．该测线的视电阻率拟断面图和阻抗相位拟图６　层状电性各向同性介质中电性各向异性介质地电模型的大地电磁响应拟断面图（ａ）ＴＥ模式视电阻率拟断面图；（ｂ）ＴＭ模式视电阻率拟断面图；（ｃ）ＴＥ模式阻抗相位拟断面图；（ｄ）ＴＭ模式阻抗相位拟断面图．Ｆｉｇ．６　Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｒｅｓｐｏｎｓｅ　ｏｆ　ｅｌｅｃｔｒｉｃａｌ　ａｎｉｓｏｔｒｏｐｉｃｍｅｄｉｕｍ　ｇｅｏｅｌｅｃｔｒｉｃ　ｍｏｄｅｌ　ｉｎ　ｔｈｅ　ｌａｙｅｒｅｄ　ｉｓｏｔｒｏｐｉｃ　ｍｅｄｉｕｍ（ａ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ａｐｐａｒｅｎｔ　ｒｅｓｉｓｔｉｖｉｔｙ　ｉｎ　ＴＥ　ｍｏｄｅ；（ｂ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ａｐｐａｒｅｎｔ　ｒｅｓｉｓｔｉｖｉｔｙ　ｉｎ　ＴＭ　ｍｏｄｅ；（ｃ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ｉｍｐｅｄａｎｃｅ　ｐｈａｓｅ　ｉｎ　ＴＥ　ｍｏｄｅ；（ｄ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ｉｍｐｅｄａｎｃｅ　ｐｈａｓｅ　ｉｎ　ＴＭ　ｍｏｄｅ．图７　测线位置图Ｆｉｇ．７　Ｔｈｅ　ｌｏｃａｔｉｏｎ　ｇｒａｐｈ　ｏｆ　ｓｕｒｖｅｙ　ｌｉｎｅ４７２４　１２期胡祥云等：大地电磁各向异性二维模拟及实例分析图９　二维电性各向异性地电模型Ｆｉｇ．９　Ｔｈｅ　ｔｗｏ－ｄｉｍｅｎｓｉｏｎａｌ　ｅｌｅｃｔｒｉｃａｌ　ａｎｉｓｏｔｒｏｐｉｃ　ｇｅｏｅｌｅｃｔｒｉｃ　ｍｏｄｅｌ断面图如图８所示（图中视电阻率拟断面图单位为Ωｍ，相位拟断面图单位为°）．因为测点较少，必须通过插值成图，造成了阻抗相位拟断面图的成图效果较差，但视电阻率拟断面图的成图效果还是能够满足研究要求的．在以往的研究中，基本上是通过对整条测线的大地电磁数据做ＴＥ、ＴＭ或ＴＥ＋ＴＭ反演，从而研究测区的地下电性结构，反演方法也多种多样（Ｏｃｃａｍ［１５－１６］、ＲＲＩ［１７］、ＮＬＣＧ［１８－１９］等），这些研究的理论基础都是建立在电性各向同性理论上．本文尝试通过建立一个二维电性各向异性地电模型，基于前文所阐述的电性各向异性理论，对其做数值模拟计算，将计算结果同图８做对比，从而分析该测区的电性结构．建立的二维电性各向异性地电模型如图９所示．由前文所述，整条测线的高程起伏较小，故图９中的地电模型不包含地形起伏．图９中该地电模型包含６个区域（分别由圆圈中的数字所指示），其中２、３、６、８区域为电性各向异性区域，４、５为电性各向同性区域．模型参数如下所示：介质２：ρ１２＝３６Ωｍ，ρ２２＝３０Ωｍ，ρ３２＝３６Ωｍ，α＝２０°；介质３：ρ１３＝３０Ωｍ，ρ２３＝２０Ωｍ，ρ３３＝３０Ωｍ，α＝３０°；介质６：ρ１６＝２Ωｍ，ρ２６＝１０００００Ωｍ，ρ３６＝２Ωｍ，α＝０°；介质８：ρ１８＝１０Ωｍ，ρ２８＝７Ωｍ，ρ３８＝１０Ωｍ，α＝０°；介质４：ρ４＝２Ωｍ；介质５：ρ５＝１５Ωｍ．各个介质的区域范围如图９所示，这里不再一一叙述．正演计算时将模型剖分为８１×５５的剖分单元（垂向５５个剖分单元中包含１０个空气层的剖分单元），测点数为８１个（即在地表的网格节点处），采用２８个周期点（０．１，０．３，０．５，０．７，１，３，５，７，１０，２０，３０，４０，５０，６０，７０，８０，９０，１００，２００，４００，５００，６００，７００，９００，１０００，１３００，１５００，２０００，单位为ｓ）．计算结果如图１０所示（图中视电阻率拟断面图单位为Ωｍ，相位拟断面图单位为°）．为方便对比，图１０使用了同图８一致的色标，可以直观看到，不论是ＴＥ模式还是ＴＭ模式，视电阻率拟断面图都有很好的拟合；阻抗相位拟断面图在整体上差强人意，主要的差别出现在ｌｏｇ（ｆ）≥０和５０～７０ｋｍ之间的区域．对于整个低阻背景场（２Ωｍ）来说，存在高阻异常区域的３个主要部分分别是：ｌｏｇ（ｆ）≥０和０～５０ｋｍ之间的区域；ｌｏｇ（ｆ）≥０和７５～８０ｋｍ之间的区域；－１≤ｌｏｇ（ｆ）＜０和２０～５０ｋｍ之间的区域．无论是高阻异常区域的视电阻率值，还是异常区域的空间位置以及在水平和深度方向的延伸，图１０中的视电阻率拟断面图中都有很好的反映．由此可以判断，所建立的二维电性各向异性地电模型能够反映出测区的电性结构．由图９可知，该区域的电性结构相当复杂，存在多个电性各向异性区域，各个电性区域间的相互影响也非常大，特别是在－１≤ｌｏｇ（ｆ）＜０和２０～５０ｋｍ之间出现的高阻异常区域，此区域在ＴＥ模式中没有任何异常出现，而在ＴＭ模式中变为一个非常明显的高阻异常区（相对于背景场值），这一现象本文在前文中有过详细阐述，属于电性各向异性结构的特例．在电磁法中，低阻背景场值会掩盖很多有用的电性信息，对高阻区域的影响尤为明显．在建立模型时，需要适当的扩大各个高阻区域的范围．在数值模拟计算时，对于此类背景场值和各个区域的电阻率值等级相差很大的情况，需要细分剖分网格，并加密频率，才能够将计算结果在成图后，很好的展现出各个不同电性区域的边界和范围．对于此测区所做的二维电性各向异性结构的研究，不仅验证了二维电性各向异性基础理论及算法的正确性，也说明了在大地电磁资料中普遍存在电性各向异性现象，同时为今后的实际大地电磁资料处理提供了一种新的思路和方法．５７２４地球物理学报（Ｃｈｉｎｅｓｅ　Ｊ．Ｇｅｏｐｈｙｓ．）５６卷图８　测线的视电阻率拟断面图和阻抗相位拟断面图（ａ）ＴＥ模式视电阻率拟断面图；（ｂ）ＴＭ模式视电阻率拟断面图；（ｃ）ＴＥ模式阻抗相位拟断面图；（ｄ）ＴＭ模式阻抗相位拟断面图．Ｆｉｇ．８　Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ａｐｐａｒｅｎｔ　ｒｅｓｉｓｔｉｖｉｔｙ　ａｎｄ　ｉｍｐｅｄａｎｃｅ　ｐｈａｓｅ　ｏｆ　ｓｕｒｖｅｙ　ｌｉｎｅ（ａ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ａｐｐａｒｅｎｔ　ｒｅｓｉｓｔｉｖｉｔｙ　ｉｎ　ＴＥ　ｍｏｄｅ；（ｂ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ａｐｐａｒｅｎｔ　ｒｅｓｉｓｔｉｖｉｔｙ　ｉｎ　ＴＭ　ｍｏｄｅ；（ｃ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ｉｍｐｅｄａｎｃｅ　ｐｈａｓｅ　ｉｎ　ＴＥ　ｍｏｄｅ；（ｄ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ｉｍｐｅｄａｎｃｅ　ｐｈａｓｅ　ｉｎ　ＴＭ　ｍｏｄｅ．图１０　二维电性各向异性地电模型的大地电磁响应拟断面图（ａ）ＴＥ模式视电阻率拟断面图；（ｂ）ＴＭ模式视电阻率拟断面图；（ｃ）ＴＥ模式阻抗相位拟断面图；（ｄ）ＴＭ模式阻抗相位拟断面图．Ｆｉｇ．１０　Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｒｅｓｐｏｎｓｅ　ｏｆ　ｔｗｏ－ｄｉｍｅｎｓｉｏｎ　ｅｌｅｃｔｒｉｃａｌ　ａｎｉｓｏｔｒｏｐｉｃ　ｇｅｏｅｌｅｃｔｒｉｃ　ｍｏｄｅｌ（ａ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ａｐｐａｒｅｎｔ　ｒｅｓｉｓｔｉｖｉｔｙ　ｉｎ　ＴＥ　ｍｏｄｅ；（ｂ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ａｐｐａｒｅｎｔ　ｒｅｓｉｓｔｉｖｉｔｙ　ｉｎ　ＴＭ　ｍｏｄｅ；（ｃ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ｉｍｐｅｄａｎｃｅ　ｐｈａｓｅ　ｉｎ　ＴＥ　ｍｏｄｅ；（ｄ）Ｔｈｅ　ｐｓｅｕｄｏ－ｓｅｃｔｉｏｎ　ｏｆ　ｉｍｐｅｄａｎｃｅ　ｐｈａｓｅ　ｉｎ　ＴＭ　ｍｏｄｅ．５　结　论本文基于Ｐｅｋ等（１９９７）［１０］的研究成果，系统而详实的阐述了二维电性各向异性正演理论；并将电性各向异性理论成功的应用于大地电磁实测资料的处理解释中，为今后大地电磁资料的处理解释提供了新的思路和方法，主要取得了以下成果：１）从麦克斯韦方程组出发，采用矩形剖分方法，使用有限差分法为正演研究工具，对二维大地电磁电性各向异性正演问题进行了系统研究．通过对二维电性各向异性地电模型的正演计算，研究了二维电性各向异性结构对观测大地电磁场的影响，为今后大地电磁实测资料的处理解释提供理论依据．６７２４　１２期胡祥云等：大地电磁各向异性二维模拟及实例分析２）以本文的研究成果为基础，将电性各向异性理论引入对实测大地电磁资料的处理解释中，通过对新疆古尔班通古特沙漠边缘大地电磁资料做二维正演拟合解释，说明了电性各向异性现象是普遍存在的，也验证了理论的正确性及算法的实用性，能够为今后分析和解释大地电磁资料中的电性各向异性现象提供理论依据和技术指导，开拓了对大地电磁实测资料处理解释的思路和方法．致　谢　感谢捷克科学院地球物理研究所ＪｏｓｅｆＰｅｋ教授对本文的帮助．参考文献（Ｒｅｆｅｒｅｎｃｅｓ）［１］　王家映．我国大地电磁测深研究新进展．地球物理学报，１９９７，４０（增刊）：２０６－２１６． Ｗａｎｇ　Ｊ　Ｙ．Ｎｅｗ　ｄｅｖｅｌｏｐｍｅｎｔ　ｏｆ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｓｏｕｎｄｉｎｇ　ｉｎＣｈｉｎａ．Ｃｈｉｎｅｓｅ　Ｊ．Ｇｅｏｐｈｙｓ（ｉｎ　Ｃｈｉｎｅｓｅ），１９９７，４０（Ｓｕｐｐｌ）：２０６－２１６．［２］　魏文博．我国大地电磁测深新进展及瞻望．地球物理学进展，２００２，１７（２）：２４５－２５４． Ｗｅｉ　Ｗ　Ｂ．Ｎｅｗ　Ａｄｖａｎｃｅ　ａｎｄ　Ｐｒｏｓｐｅｃｔ　ｏｆ　ＭａｇｎｅｔｏｔｅｌｌｕｒｉｃＳｏｕｎｄｉｎｇ（ＭＴ）ｉｎ　Ｃｈｉｎａ．Ｐｒｏｇｒｅｓｓ　ｉｎ　Ｇｅｏｐｈｙｓｉｃｓ（ｉｎ　Ｃｈｉｎｅｓｅ），２００２，１７（２）：２４５－２５４．［３］　赵国泽，汤吉，詹艳等．青藏高原东北缘地壳电性结构和地块变形关系的研究．中国科学（Ｄ辑），２００４，３４（１０）：９０８－９１８． Ｚｈａｏ　Ｇ　Ｚ，Ｔａｎｇ　Ｊ，Ｚｈａｎ　Ｙ，ｅｔ　ａｌ．Ｓｔｕｄｙ　ｏｎ　ｔｈｅ　ｃｒｕｓｔａｌｅｌｅｃｔｒｉｃａｌ　ｓｔｒｕｃｔｕｒｅ　ａｎｄ　ｂｌｏｃｋ　ｄｅｆｏｒｍａｔｉｏｎ　ｉｎ　Ｎｏｒｔｈｅａｓｔｍａｒｇｉｎ　ｏｆ　Ｔｉｂｅｔａｎ　Ｐｌａｔｅａｕ．Ｓｃｉｅｎｃｅ　ｉｎ　Ｃｈｉｎａ（Ｓｅｒ．Ｄ）（ｉｎＣｈｉｎｅｓｅ），２００４，３４（１０）：９０８－９１８．［４］　Ｇｕｏｚｅ　Ｚｈａｏ　ｅｔ　ａｌ，Ｃｒｕｓｔａｌ　ｓｔｒｕｃｔｕｒｅ　ａｎｄ　ｒｈｅｏｌｏｇｙ　ｏｆ　ｔｈｅＬｏｎｇｍｅｎｓｈａｎ　ａｎｄ　Ｗｅｎｃｈｕａｎ　Ｍｗ７．９ｅａｒｔｈｑｕａｋｅ　ｅｐｉｃｅｎｔｒａｌａｒｅａ　ｆｒｏｍ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｄａｔａ，ＧＥＯＬＯＧＹ，Ｄｅｃｅｍｂｅｒ　２０１２；４０（１２）：１１３９－１１４２．［５］　Ｗａｎｎａｍａｋｅｒ　Ｐ．Ａｎｉｓｏｔｒｏｐｙ　Ｖｅｒｓｕｓ　Ｈｅｔｅｒｏｇｅｎｅｉｔｙ　ｉｎ　ＣｏｎｔｉｎｅｎｔａｌＳｏｌｉｄ　Ｅａｒｔｈ　Ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃ　Ｓｔｕｄｉｅｓ：Ｆｕｎｄａｍｅｎｔａｌ　ＲｅｓｐｏｎｓｅＣｈａｒａｃｔｅｒｉｓｔｉｃｓ　ａｎｄ　Ｉｍｐｌｉｃａｔｉｏｎｓ　ｆｏｒ　Ｐｈｙｓｉｃｏｃｈｅｍｉｃａｌ　Ｓｔａｔｅ．ＳｕｒｖＧｅｏｐｈｙｓ，２００５，２６（６）：７３３－７６５．［６］　霍光谱，胡祥云，刘敏．各向异性介质中大地电磁正演研究综述．地球物理学进展，２０１１，２６（６）：１９７６－１９８２． Ｈｕｏ　Ｇ　Ｐ，Ｈｕ　Ｘ　Ｙ，Ｌｉｕ　Ｍ．Ｒｅｖｉｅｗ　ｏｆ　ｔｈｅ　ｆｏｒｗａｒｄ　ｍｏｄｅｌｉｎｇｏｆ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｉｎ　ｔｈｅ　ａｎｉｓｏｔｒｏｐｙ　ｍｅｄｉｕｍ　ｒｅｓｅａｒｃｈ．Ｐｒｏｇｒｅｓｓ　ｉｎ　ｇｅｏｐｈｙｓｉｃｓ（ｉｎ　Ｃｈｉｎｅｓｅ），２０１１，２６（６）：１９７６－ １９８２．［７］　Ｒｅｄｄｙ　ＩＫ，Ｒａｎｋｉｎ　Ｄ．Ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｒｅｓｐｏｎｓｅ　ｏｆ　ｌａｔｅｒａｌｌｙｉｎｈｏｍｏｇｅｎｅｏｕｓ　ａｎｄ　ａｎｉｓｏｔｒｏｐｉｃ　ｍｅｄｉａ．Ｇｅｏｐｈｙｓｉｃｓ，１９７５，４０（６）：１０３５－１０４５．［８］　徐世浙，赵生凯．二维各向异性地电断面大地电磁场的有限元法解法．地震学报，１９８５，０７（０１）：８０－９０． Ｘｕ　Ｓ　Ｚ，Ｚｈａｏ　Ｓ　Ｋ．Ｓｏｌｕｔｉｏｎ　ｏｆ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｆｉｅｌｄ　ｅｑｕａｔｉｏｎｓｆｏｒ　ａ　ｔｗｏ－ｄｉｍｅｎｓｉｏｎａｌ，ａｎｉｓｏｔｒｏｐｉｃ　ｇｅｏｅｌｅｃｔｒｉｃ　ｓｅｃｔｉｏｎ　ｂｙ　ｔｈｅｆｉｎｉｔｅ　ｅｌｅｍｅｎｔ　ｍｅｔｈｏｄ．Ａｃｔａ　Ｓｅｉｓｍｏｌｏｇｉｃａ　Ｓｉｎａｃａ（ｉｎＣｈｉｎｅｓｅ），１９８５，０７（０１）：８０－９０．［９］　Ｌｉ　Ｙ．Ａ　ｆｉｎｉｔｅ－ｅｌｅｍｅｎｔ　ａｌｇｏｒｉｔｈｍ　ｆｏｒ　ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃ　ｉｎｄｕｃｔｉｏｎｉｎ　ｔｗｏ－ｄｉｍｅｎｓｉｏｎａｌ　ａｎｉｓｏｔｒｏｐｉｃ　ｃｏｎｄｕｃｔｉｖｉｔｙ　ｓｔｒｕｃｔｕｒｅｓ．Ｇｅｏｐｈｙｓ　ＪＩｎｔ，２００２，１４８（３）：３８９－４０１．［１０］　Ｐｅｋ　Ｊ，Ｖｅｒｎｅｒ　Ｔ．Ｆｉｎｉｔｅ－ｄｉｆｆｅｒｅｎｃｅ　ｍｏｄｅｌｌｉｎｇ　ｏｆ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃｆｉｅｌｄｓ　ｉｎ　ｔｗｏ－ｄｉｍｅｎｓｉｏｎａｌ　ａｎｉｓｏｔｒｏｐｉｃ　ｍｅｄｉａ．Ｇｅｏｐｈｙｓ　Ｊ　Ｉｎｔ，１９９７，１２８（３）：５０５－５２１．［１１］　Ｌｉ　Ｙ，Ｐｅｋ　Ｊ．Ａｄａｐｔｉｖｅ　ｆｉｎｉｔｅ　ｅｌｅｍｅｎｔ　ｍｏｄｅｌｌｉｎｇ　ｏｆ　ｔｗｏ－ｄｉｍｅｎｓｉｏｎａｌ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｆｉｅｌｄｓ　ｉｎ　ｇｅｎｅｒａｌ　ａｎｉｓｏｔｒｏｐｉｃｍｅｄｉａ．Ｇｅｏｐｈｙｓ　Ｊ　Ｉｎｔ，２００８，１７５（３）：９４２－９５４．［１２］　Ｖ　Ｈａａｋ．Ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃｓ：Ｄｅｔｅｒｍｉｎａｔｉｏｎ　ｏｆ　ｔｈｅ　ｔｒａｎｓｆｅｒｆｕｎｃｔｉｏｎｓ　ｉｎ　ｒｅｇｉｏｎｓ　ｗｉｔｈ　ｌａｔｅｒａｌ　ｃｈａｎｇｅｓ　ｏｆ　ｔｈｅ　ｅｌｅｃｔｒｉｃａｌｃｏｎｄｕｃｔｉｖｉｔｙ．Ｚ　Ｇｅｏｐｈｙｓ（ｉｎ　Ｇｅｒｍａｎ），１９７２，３８：８５－１０２．［１３］　Ｂｒｅｗｉｔｔ－Ｔａｙｌｏｒ　Ｃ　Ｒ，Ｗｅａｖｅｒ　Ｊ　Ｔ．Ｏｎ　ｔｈｅ　ｆｉｎｉｔｅ　ｄｉｆｆｅｒｅｎｃｅｓｏｌｕｔｉｏｎ　ｏｆ　ｔｗｏ－ｄｉｍｅｎｓｉｏｎａｌ　ｉｎｄｕｃｔｉｏｎ　ｐｒｏｂｌｅｍｓ．Ｇｅｏｐｈｙｓ　Ｊ．Ｒ．ａｓｔｒ．Ｓｏｃ．，１９７６，４７（２）：３７５－３９６．［１４］　ˇＣｅｒｖ　Ｖ，Ｐｒａｕｓ　Ｏ，Ｈｖｏˇｚｄａｒａ　Ｍ．Ｎｕｍｅｒｉｃａｌ　ｍｏｄｅｌｌｉｎｇ　ｉｎｌａｔｅｒａｌｌｙ　ｉｎｈｏｍｏｇｅｎｅｏｕｓ　ｇｅｏｅｌｅｃｔｒｉｃａｌ　ｓｔｒｕｃｔｕｒｅｓ．Ｓｔｕｄｉａ　Ｇｅｏｐｈｙｓ．ｅｔ　Ｇｅｏｄａｅｔ．，１９７８，２２（１）：７４－８１．［１５］　Ｃｏｎｓｔａｂｌｅ　Ｓ　Ｃ，Ｐａｒｋｅｒ　Ｒ　Ｌ，Ｃｏｎｓｔａｂｌｅ　ＣＧ．Ｏｃｃａｍ′ｓｉｎｖｅｒｓｉｏｎ；ａｐｒａｃｔｉｃａｌ　ａｌｇｏｒｉｔｈｍ　ｆｏｒ　ｇｅｎｅｒａｔｉｎｇ　ｓｍｏｏｔｈ　ｍｏｄｅｌｓｆｒｏｍ　ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃ　ｓｏｕｎｄｉｎｇ　ｄａｔａ．Ｇｅｏｐｈｙｓｉｃｓ，１９８７，５２（３）：２８９－３００．［１６］　ｄｅ　Ｇｒｏｏｔ－Ｈｅｄｌｉｎ　ＣＤ，Ｃｏｎｓｔａｂｌｅ　ＳＣ．Ｏｃｃａｍ＇ｓ　ｉｎｖｅｒｓｉｏｎ　ｔｏｇｅｎｅｒａｔｅ　ｓｍｏｏｔｈ，ｔｗｏ－ｄｉｍｅｎｓｉｏｎａｌ　ｍｏｄｅｌｓ　ｆｒｏｍ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃｄａｔａ．Ｇｅｏｐｈｙｓｉｃｓ，１９９０，５５（１２）：１６１３－１６２４．［１７］　Ｓｍｉｔｈ　ＪＴ，Ｂｏｏｋｅｒ　ＪＲ．Ｒａｐｉｄ　Ｉｎｖｅｒｓｉｏｎ　ｏｆ　Ｔｗｏ－ａｎｄ　Ｔｈｒｅｅ－Ｄｉｍｅｎｓｉｏｎａｌ　Ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　Ｄａｔａ．Ｊ　Ｇｅｏｐｈｙｓ　Ｒｅｓ，１９９１，９６（Ｂ３）：３９０５－３９２２．［１８］　Ｆｌｅｔｃｈｅｒ　Ｒ，Ｒｅｅｖｅｓ　ＣＭ．Ｆｕｎｃｔｉｏｎ　ｍｉｎｉｍｉｚａｔｉｏｎ　ｂｙ　ｃｏｎｊｕｇａｔｅｇｒａｄｉｅｎｔｓ．Ｃｏｍｐｕｔｅｒ　Ｊ．，１９６４，７（２）：１４９－１５４．［１９］　Ｒｏｄｉ　Ｗ，Ｍａｃｋｉｅ　ＲＬ．Ｎｏｎｌｉｎｅａｒ　ｃｏｎｊｕｇａｔｅ　ｇｒａｄｉｅｎｔｓ　ａｌｇｏｒｉｔｈｍｆｏｒ　２－Ｄ　ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ　ｉｎｖｅｒｓｉｏｎ．Ｇｅｏｐｈｙｓｉｃｓ，２００１，６６（１）：１７４－１８７．（本文编辑　刘少华）７７２４。