带地形的大地电磁各向异性二维模拟及实例对比分析

霍光谱, 胡祥云, 黄一凡, 韩波

1 河南省地质调查院， 郑州 4500002 中国地质大学地球物理与空间信息学院，地球内部多尺度成像湖北省重点实验室， 武汉 430074



带地形的大地电磁各向异性二维模拟及实例对比分析

霍光谱1, 胡祥云2*, 黄一凡2, 韩波2

1 河南省地质调查院， 郑州 4500002 中国地质大学地球物理与空间信息学院，地球内部多尺度成像湖北省重点实验室， 武汉 430074

带地形的大地电磁二维正演数值模拟多数基于电性各向同性理论，由于地球内部电性各向异性现象的普遍存在，基于电性各向异性理论研究地形起伏情况下大地电磁二维正演数值模拟就显得非常迫切.本文首先由麦克斯韦方程出发，引入张量电导率，求得一组关于平行走向的电场分量Ex和磁场分量Hx的二阶偏微分方程，使用有限差分法求解出Ex和Hx的近似解，并以此求得其他场分量；其次，引入地形因素，改变变量在网格节点中的排列方式，选择交错排列方式从而给有限差分系数矩阵的最大带宽分配合理的存储空间；最后，使用Weaver的方法解决TM模式下，在地-空分界面垂直于构造走向的一些区域存在不同电导率的问题.通过对带地形的二维电性各向异性结构做正演模拟，研究地形因素对大地电磁响应的影响；以电性各向异性理论为基础，将地形因素引入对实测大地电磁资料的处理中，通过做二维正演拟合和未引入地形因素的结果做对比，说明电性各向异性现象的普遍存在，认识地形因素对观测大地电磁场的影响，为今后分析解释实测大地电磁资料包含地形因素和电性各向异性情况提供理论基础和技术指导.

地形； 电性各向异性； 大地电磁； 有限差分； 张量电导率

1 引言

20世纪50年代初Tikhonov（1950）和Cagnird（1953）分别提出利用天然电磁场进行勘探的方法，从70年代开始大地电磁法进入飞速发展时期，到目前已经成为国内外学者（Weiss and Newman，2002；赵国泽等，2004；Wannamaker et al，2005；魏文博等，2006,2010；胡祥云等，2013）研究地球内部电性结构最常用的方法之一.Ku等（1973）发现地形对大地电磁场的影响，20世纪80年代开始，不断有学者将地形因素引入大地电磁数值模拟中， Wannamaker等（1986）用有限元法对起伏地形下二维大地电磁响应进行数值模拟.Chouteau和Bouehard（1988）研究了二维大地电磁响应的地形校正问题.徐世浙等（1985,1992,1997）用有限元和边界元法计算起伏地形下水平层状介质的大地电磁响应，并对畸变后的大地电磁响应做地形改正，这些研究成果均是基于电性各向同性理论.

20世纪60年代末，Parkhomenko（1967），Keller（1968）和Hill（1972）分别讨论了岩石中的微各向异性问题，与此同时，电性各向异性现象的理论研究也逐渐发展起来.O′Brienh和 Morrison（1967）研究层状各向异性介质中电磁场的递推公式；Reddy和Rankin （1971）研究倾斜各向异性对大地电磁场的影响.Reddy和Rankin（1975）用有限元法对偏微分方程做数值计算，较早对二维电性各向异性问题进行了理论研究.虽然在其后的20多年间（直到2000之后才有学者将其理论改进后用于大地电磁实测资料处理）未能广泛运用在大地电磁实测资料处理中，但其所做研究在理论意义上的影响却尤为深远，不仅让世人重视电性各向异性现象，更为今后的二维电性各向异性研究奠定了理论基础.Pek和Toh（1997）在Reddy和Rankin（1975）研究的基础上，并将前人（Haak，1972；Brewitt-Taylor and Wearver，1976；Cerv et al.，1978）在二维电性各向同性介质中使用的有限差分法引入求解二维电性各向异性问题中，用有限差分法对偏微分方程做数值模拟计算，是正演计算二维电性各向异性介质大地电磁响应方法中影响最为深远和广泛的一种求解方法.Li（2002，2008）细述用有限元法计算二维电性各向异性结构感应电磁场，并于2008年用自适应网格剖分法代替以往的网格剖分法，是一次成功的改进，取得显著的效果.霍光谱等（2012）通过引入权因子成功改进马奎特反演理论，将基于电性各向同性理论的反演方法扩展应用到电性各向异性理论.

基于电性各向异性理论带地形的大地电磁数值模拟在最近20年间才逐渐发展起来，用于实际资料处理的更为少见.Pek和Toh（2000）将公式扩展应用在大地含地形和海洋测深中.Key和Ovall（2011）使用有限元法采用自适应非结构性网格剖分方式，系统研究二维大地电磁及海洋可控源电磁法，是最近较为系统而完善的方法.

2 带地形的大地电磁电性各向异性二维正演理论

本文采用有限差分法参照Pek等（1997,2000,2002）的方法研究带地形的大地电磁二维电性各向异性正演问题.

假设在笛卡尔坐标系中，有一个二维电性结构，它的构造走向平行于x轴，y垂直于x轴，z轴垂直于xy平面，且正向向下.模型由一个有限区域组成，在其中的电性各向异性区域是二维有限结构.假设地表为水平，对应z=0.在地表上空（z<0）为绝缘空气层.来自z→-∞的平面波垂直向下传播，ω=2π/T，T为周期，单位为秒（s）.

谐变电磁场的麦克斯韦方程组为

（1）

（2）

上式为描述电磁场分布的微分方程.e-iω t为时谐因子，σ为张量电导率.

张量电导率σ为三阶对称矩阵，可分解为由3个电导率σ1,σ2,σ3和一个旋转矩阵，旋转矩阵又可通过3个Euler′s空间旋转元素依次旋转角度为αL,αD,αS求得，角度αS,αD,αL分别定义为各向异性介质的走向，倾角，倾向（图1所示）.

（3）

其中R为旋转矩阵元素：

（4）

通过公式（3）可以表示空间中任何位置的电性各向异性结构.

由2-D性质∂/∂x=0，将公式简化为只含Ex和Hx的一组二阶偏微分方程（胡祥云等，2013）.

电场沿走向时（即TE模式）：

（5）

磁场沿走向时（即TM模式）：

图1 各向异性基本参数示意图（依次运用三个Euler′s旋转元素αS,αD,αL可将导电面变化到空间任何位置）

（6）

公式（5）—（6）中，

A=（σxyσyz-σxzσyy）/D,B=（σxyσyz-σxzσzz）/D,

用有限差分法求解公式（5）—（6），求得网格节点（j,k）上TE模式的近似差分方程，这个近似差分方程代表公式（5）的有限差分形式：

（7）

公式（7）中，

（p,q）∈{（j-1,k-1）,（j-1,k+1）,（j+1,

k-1）,（j+1,k+1）}

同样可以求得在网格节点（j,k）上TM模式下的线性差分方程：

（8）

在公式（8）中，

（p,q）∈{（j-1,k-1）,（j-1,k+1）,（j+1,

k-1）,（j+1,k+1）}.

用有限差分法对公式（5）—（6）做数值近似后，须将线性代数公式（7）—（8）安排在一个系统，以便进一步处理.因为包含Ex和Hx两个变量，且这两个变量多数时候并不同时出现在同一个网格节点，TE模式的方程需要在全部网格节点做近似计算，而TM模式的方程只需在导电底层内部做近似计算.

引入地形因素时，对应的有限差分网格可通过起伏地形对应的阶梯函数来近似.并不影响用有限体积法在中心网格节点的任何位置对公式（5）—（6）做近似计算.将变量排入一个系统有两种方式：顺序排列和交错排列.本文选取交错排列方法.

图2a是水平地形，对应的存储矩阵是2c；2b是起伏地形，对应的存储矩阵是2d.由图2b可以看到，在网格中，不同的列，Hx的个数是可变的，所以图2d中有限差分系数矩阵的带宽是可变的：带宽的最宽部分由地形的海拔高度决定；带宽的最窄部分由地形的盆地深度决定.这也是起伏地形相对于水平地形不同的地方，带地形情况下，就必须考虑对公式做相应的有效修改，即：在用高斯消元法对有限差分方程组求解时，须为有限差分系数矩阵的最大带宽分配相应的存储空间.

与水平地表相比，含起伏地形时，TM模式在地-空分界面（y方向上）存在一些不同电导率的区域，如：σj,kσj+1,k+1≠σj+1,kσj,k+1.在这些区域的节点上，没有连续的电流和连续的切向电流.为解决这个问题，须将Weaver等（1985,1986）的光滑方法应用在这些网格节点的局部导电区域.

Weaver等（1985,1986）认识到用矩形网格剖分模型时，有时网格节点四周的4个单元格会存在不同的电导率，并给出解决此问题的方法：

1） 对选定的网格节点周围的局部电导率重新塑造为光滑的电导率分布，这并不影响用有限差分法近似求解偏微分方程.

2） 使用同样的近似方法（有限差分法）求原始偏微分方程的近似值.

在用有限差分方法对二阶偏微分方程（5）—（6）做近似计算后，可以得到有限差分线性方程组，并求出整个模型中各个节点上的Ex和Hx的近似值.进而可以求出磁场分量Hy和Hz.

求解过程为：

1） 将场分量Ex在中心节点的两个负方向上，分别在各自的网格步长中用泰勒公式展开到二阶.

2） 用二阶导数替换这些展开项，它们等于原来的二阶方程.

3） 在中心节点的两边求出合适的电导率平均值.

4） 通过从先前步骤中减去适当比例方程，消除Ex剩下的二阶导数.

在电性各向异性情况中，磁场分量Hx的影响也同样必须重新考虑.它会从本质上影响计算电导率平均值的过程.

在积分单元上对电导率或电阻率求平均值是可以重复多次进行的，但之间的差别很微小，在电性各向异性结构中，电流密度的各个分量不是对应的电场分量乘以标量电导率，而是所有电场分量乘以一个电导率张量.在分界面的边界条件，须对各个区域的不同电导率张量进行合适的修改.

下面将表明通过使用有限体积近似法，Weaver等（1986）的求导公式可用于电性各向异性情况.例如：Ex在z方向上的导数∂Ex（yj,zk）/∂z.对公式

图2 变量在网格节点中的交错排列方式

（5）用体积积分法分别计算积分单元Gj,k的两个子单元，即：积分区域的上半部分（z>zk）和下半部分（z

（9）

式中，

（10）

公式（10）前两项只包含电场分量Ex，这与Weaver等（1986）用于计算电性各向同性情况下改进的求导公式一致.公式（10）中的第一项代表一个导数值，它由网格节点（j,k-1）,（j,k）,（j,k+1）对Ex做抛物差值计算求出.第二项是一个修改量，它是原始偏微分方程的结果，代表了二阶偏导数∂2Ex/∂z2的值.这个修改量与中心节点上下平均电导率的差成正比.对于存在较大电导率差值时，忽视这个修改量会造成导出的场值精度快速下降，即使场分量Ex非常准确并有着很高的精度.

公式（10）的后两项包含了磁场分量Hx，表示了在电性各向异性情况下，H极化模式的影响.这些项由电导率张量的非对角元素由集合参数A和B及它们的平均值B-和B+所决定.

同样的，重复以上过程，可以容易求得积分区域Gj,k的左右积分子区域IE（Gj,k,yj-）和IE（Gj,k,yj+），并得到∂Ex（j,k）/∂y的近似值.

随后可以计算由麦克斯韦方程组导出的电场分量Ey和Ez，将有限体积积分用于公式（6）中，求出∂Hx/∂y和∂Hx/∂z.

（11）

3 带地形的二维电性各向异性模型计算与分析

在数据采集过程中，数据质量是基础是关键，大地电磁方法应用成功的与否在很大程度上取决于数据采集质量.我国地域辽阔，地形复杂多变，工业生产、人文因素引起的强干扰区较为普遍.

大地电磁法多数是研究大尺度构造，虽然在多数情况下，地形的起伏相对于测线长度变化较小，但也存在一些情况，例如在构造带附近，地形会有较为剧烈的变化，此时就要求依据实际情况建立含有起伏地形的地电模型用于处理解释.对含地形因素的二维电性各向异性地电模型做正演研究，可以为分析研究存在二维电性各向异性介质情况下，地形因素对大地电磁场产生的影响提供理论支撑.

为说明地形的影响，首先对均匀半空间包含地垒和地堑两种地形的情况分别做正演计算.均匀半空间ρ=100Ωm（图3所示）.

图3 含地形的均匀半空间

正演计算时将模型剖分为69×47的剖分单元（垂向47个剖分单元），测点数为69个（即在地表的网格节点处），采用21个周期点（0.1, 0.3, 0.5, 0.7, 1, 3, 5, 7, 10, 30,50,70,100,300,500,700,1000,1300,1500,1700,2000）.计算结果如图4、图5所示.

图4、图5分别是均匀半空间中地垒地形和地堑地形的大地电磁响应结果.可以直观看出：无论是地垒还是地堑，地形的影响主要表现在地形起伏区域，突出的畸变出现在地形剧烈变化的交汇处（平地与山坡交界及山坡与山顶交界处）.在地垒地形的起伏区域，TE模式的视电阻率值和阻抗相位值呈现出小-大-小的变化，TM模式的视电阻率值和阻抗相位值呈现出大-小-大的变化.而在地堑地形的起伏区域，TE模式和TM模式中视电阻率值和阻抗相位值的变化正好与地垒情况相反.此外，无论是地垒还是地堑，地形的影响对于视电阻率值来说主要集中在TE模式的中高频段和TM模式的低频段，而对于阻抗相位来说影响几乎都是在中高频段.最后，地形的影响程度和地形倾角有关.倾角越大，影响就越大；反之，则越小.

图4 含地垒地形的均匀半空间的大地电磁响应拟断面图

图5 含地堑地形的均匀半空间的大地电磁响应拟断面图

网格剖分情况（图7）及周期点与均匀半空间地垒相同，计算结果如图8所示.

图7 网格剖分图

将图8的结果同之前（胡祥云等，2013）未含地形的正演结果对比可知：含地垒情况下TE模式的视电阻率拟断面图中低阻异常区域在水平方向上的延展有所扩大，图8a中低阻区域展示了二维电性各向异性介质的存在范围，受地形影响在33～35km，lg（f）=0的低阻区域有明显凹陷；阻抗相位拟断面图的形态基本一致，细微区别在于含地垒情况下阻抗相位值最大值和最小值区域水平方向上的延伸稍有增加，在33～35km，lg（f）=-1.5的阻抗相位值有所增大；受地形影响TE模式下的视电阻率最小值有所增大，阻抗相位值几乎无变化.TM模式的视电阻率值和阻抗相位值受地形的影响较为明显，在山脚处（24km和44km处）视电阻率值和阻抗相位值均发生了畸变，表现为视电阻率值跳跃增大，阻抗相位值跳跃降低；在山顶处（33～35km）视电阻率值有明显降低，阻抗相位值有明显增大；受地形影响TM模式下的视电阻率最小值有所降低，阻抗相位最大值有所增大.

在含地堑地形的层状介质中放置一个二维电性各向异性介质，对地电模型做正演计算，并与之前（胡祥云等，2013）未含地形的正演结果做对比分析，介质1、2、3参数同地垒地电模型.地堑的顶部宽度为20km（从24km到44km），底部为2km，海拔高度-3km.如图9所示.

网格剖分情况及周期点与均匀半空间地堑情况相同，计算结果如图10所示.

图8 含地垒地形的层状电性各向异性地电模型的大地电磁响应拟断面图

图9 含地堑地形的层状电性各向异性地电模型

将图10的结果同之前未含地形的正演结果（胡祥云等，2013）对比可知：含地堑情况下TE模式的视电阻率拟断面图中低阻异常区域在水平方向上的延展有所扩大，图1a中低阻区域展示了二维电性各向异性介质的存在范围，受地形影响在33～35km，lg（f）=0的低阻区域有明显凸起；阻抗相位拟断面图的形态基本一致，细微区别在于含地堑情况下阻抗相位值最大值和最小值区域水平方向上的延伸稍有增加，在33～35km，lg（f）=-1位置阻抗相位值有所减小.TM模式的视电阻率值和阻抗相位值受地形的影响较为明显，在山脚处（24km和44km处）视电阻率值和阻抗相位值均发生了畸变，表现为视电阻率值跳跃减小，阻抗相位值跳跃增大；在山顶处（33～35km）视电阻率值表现为增大-减小-增大的跳跃式畸变，阻抗相位值表现为减小-增大-减小的跳跃式畸变.

总体来说，地形因素对于大地电磁响应的影响主要存在于地形的突变位置（平地与山脚交汇处及山坡与山顶交汇处），而在平稳过渡位置（山坡）对于大地电磁响应的影响较为平缓.其次，不论是地垒地形还是地堑地形，对于TE模式下的视电阻率值和阻抗相位值影响都较小；对于TM模式的影响则较大.此外，地堑地形的影响要强于地垒地形的影响，最后，地形的影响程度和地形倾角有直接关系，地形倾角越大，影响越强烈；反之则越小.

4 大地电磁实测资料解释及电性结构对比分析

选取的大地电磁测线位于贵州某地，测线位于贵州山区，地形起伏较大，加入地形的正演拟合结果展示了地形对大地电磁响应的影响.

图11是测线位置图，测线由西南向东北首先穿过山谷，随后翻过山脊进入山坡平缓区域.测线长测线的测点距为40m，对于长度为1400m的测线来说，数据量是十分充足的.图12中的4幅拟断面图都不同程度的出现了许多闭合圆圈，即使通过数据插值成图，也不能使这些区域的高/低值平滑.说明实测资料中在同一周期点，水平方向上的数据值跳动十分剧烈.本文通过构建一个含地形的二维电性各向异性地电模型，并对其做正演计算，根据正演拟合的结果分析该区域的电性结构并研究地形对大地电磁响应的影响.

图10 含地堑地形的层状电性各向异性地电模型的大地电磁响应拟断面图

图11 测线位置图

1400m，地形起伏高差240m，测点距40m.测线的视电阻率拟断面图和阻抗相位拟断面图如图12所示.

图14中TE和TM模式的正演结果受地形的影响都较大.正演计算的结果和实测数据有很好的拟合.图14a中，在lg（f）=1Hz，水平方向的200～400m，800～1400m区域出现了对应的高阻区，在X=1000m和1200m出现的高阻闭合圆圈也有很好的对应，并且高阻最大值也保持一致.图14c和图12c的形态基本保持一致，在阻抗相位的低值区域（lg（f）=2-1Hz，X=100～300m及800～1400m）都表现出高阻介质（5、6、7）在深度方向上的延展范围.图14b中，在lg（f）=0.5Hz，水平方向的100～300m，1000～1200m区域，出现高阻闭合圆圈，与图12b中的情况一致，高阻最大值也保持一致.图14d和图12d形态稍有差异，主要表现在水平方向600～800m范围内.图14d的低值区域在水平方向上很好的呈现出5、6、7三个电性各向异性介质的延展范围.根据正演拟合的结果，推断地形的突变点（X=200、700、1200m）及5、6、7三个介质的变化点（X=240、320、700、900、1140）是造成数据在同一周期点呈锯齿状起伏的主要原因.总体来说，正演拟合的结果基本上与实测数据一致，不但在位置上一一对应出高阻闭合圈的范围，而且在高阻区域的最大值也能够保持一致.

图12 带地形的二维大地电磁响应拟断面图

图13 带地形的二维电性各向异性地电模型

实测资料的正演拟合结果表明：在大尺度下地形起伏较缓时可以忽略地形因素对大地电磁响应造成的影响，而在地形起伏较大区域（如山区和构造带区域），地形变化会对大地电磁响应造成显著的影响，使得视电阻率值及阻抗相位值在地形突变位置产生畸变，此时就必须重视地形变化带来的影响，必要时需对地形造成的畸变进行校正.

5 结论

本文基于Pek（1997,2000）的研究理论，系统而详实的阐述了含地形的二维电性各向异性正演理论；并通过理论研究说明地形因素对大地电磁响应的影响，实测资料的对比研究说明地形因素对大地电磁响应的影响程度.

基于水平地形二维大地电磁各向异性研究成果，引入地形因素，改变变量在网格节点中的排列方式，选择交错排列方式从而给有限差分系数矩阵的最大带宽分配合理的存储空间，随后，使用Weaver等（1985,1986）的方法解决TM模式下，在地-空分界面Y方向上的一些区域存在不同电导率的问题，对带地形的二维大地电磁电性各向异性正演问题进行系统研究.正演计算地垒及地堑情况下的地电模型，研究地形对大地电磁响应的影响.

图14 带地形的二维电性各向异性地电模型的大地电磁响应拟断面图

以本文的研究成果为基础，将带地形的二维大地电磁电性各向异性理论引入对实测大地电磁资料的处理解释中，说明电性各向异性现象的普遍存在，同时验证理论的正确性及程序的实用性，并且对山区地形起伏剧烈情况下地形变化对大地电磁响应的影响程度深入分析，为今后处理解释大地电磁资料中包含地形因素的电性各向异性现象提供理论依据和技术指导，并开拓了对大地电磁实测资料处理的思路和方法.

致谢 感谢捷克科学院地球物理研究所JosefPek教授对本文的帮助.

Brewitt-Taylor C R, Weaver J T. 1976. On the finite difference solution of two-dimensional induction problems.Geophys.J.Int., 47（2）: 375-396.

Cagniard L. 1953. Basic theory of the magneto-telluric method of geophysical prospecting.Geophysics, 18（3）: 605-635.

Chouteau M, Bouehard K. 1988. Two-dimensional terrain correction in magnetotelluric surveys.Geophysics, 53（6）: 854-862.

Haak V. 1972. Magnetotellurics: Determination of the transfer functions in regions with lateral changes of the electrical conductivity.Z.Geophys. （in German）, 38: 85-102.

Hill D G. 1972. A laboratory investigation of electrical anisotropy in Precambrian rocks.Geophysics, 37（6）: 1022-1038.

Hu X Y, Huo G P, Gao R, et al. 2013. The magnetotelluric anisotropic two-dimensional simulation and case analysis.ChineseJ.Geophys. （in Chinese）, 56（12）: 4268-4277, doi: 10.6038/cjg20131229.

Huo G P, Hu X Y, Fang H, et al. 2012. Magnetotelluric joint inversion for anisotropic conductivities in layered media.ActaPhys.Sin. （in Chinese）, 62（12）: 129101.

Huo G P, Hu X Y, Liu M. 2011. Review of the forward modeling of magnetotelluric in the anisotropy medium research.ProgressinGeophys. （in Chinese）, 26（6）: 1976-1982,doi: 10.3969/j.issn.1004-2903.2011.06.011.

Jin G W, Zhao G Z, Xu C F, et al. 1998. The affection and correction on magnetotelluric response data for inclination two-dimension terrain.SeismologyandGeology（in Chinese）, 20（4）: 454-458.

Keller G V. 1968. Electrical prospecting for oil.∥ Quarterly of the Colorado School of Mines, v. 63, no. 2. Golden: Colorado School of Mines.Key K, Ovall J. 2011. A parallel goal-oriented adaptive finite element method for 2.5-D electromagnetic modelling.Geophys.J.Int., 186（1）: 137-154.

Ku C C, Hsieh M S, Lim S H. 1973. The topographic effect in electromagnetic fields.CanadianJ.EarthSci., 10（5）: 645-656.

Li Y G. 2002. A finite-element algorithm for electromagnetic induction in two-dimensional anisotropic conductivity structures.Geophys.J.Int., 148（3）: 389-401.

Li Y G, Pek J. 2008. Adaptive finite element modelling of two-dimensional magnetotelluric fields in general anisotropic media.Geophys.J.Int., 175（3）: 942-954.

O’Brien D P, Morrison H F. 1967. Electromagnetic fields in an N-layer anisotropic half-space.Geophysics, 32（4）: 668-677. Parkhomenko E I. 1967. Electrical Properties of Rocks. US: Springer.

Pek J, Verner T. 1997. Finite-difference modelling of magnetotelluric fields in two-dimensional anisotropic media.Geophys.J.Int., 128（3）: 505-521.

Pek J, Toh H. 2000. Numerical modeling of MT fields in 2-D anisotropic structures with topography and bathymetry considered.∥ Protokolluber das 18, Kolloquium “Electromagnetische Tiefenforschung”. Altenberg, Germany, Hordt, A., and J. Stoll, eds., DGG, 190-199.

Pek J, Santos F A M. 2002. Magnetotelluric impedances and parametric sensitivities for 1-D anisotropic layered media.ComputersandGeosciences, 28（8）: 939-950.

Reddy I K, Rankin D. 1971. Magnetotelluric effect of dipping anisotropies.GeophysicalProspecting, 19（1）: 84-97.

Reddy I K, Rankin D. 1975. Magnetotelluric response of laterally inhomogeneous and anisotropic media.Geophysics, 40（6）: 1035-1045.

Tikhonov A N. 1950. On determining electrical characteristics of the deep layers of the earth′s crust.Deki.Akud.Nuck., 73: 295-297. Wang X B, Li Y N, Gao Y C. 1999. Two dimensional topographic responses in magneto telluric sounding and it′s correction methods.ComputingTechniquesforGeophysicalandGeochemicalExploration（in Chinese）, 21（4）: 327-332.

Wannamaker P E, Stodt J A, Rijo L. 1986. Two-dimensional topographic responses in magnetotellurics modeled using finite elements.Geophysics, 51（11）: 2131-2144.

Wannamaker P E. 2005. Anisotropy versus heterogeneity in continental solid Earth electromagnetic studies: fundamental response characteristics and implications for physicochemical state.Surv.Geophys., 26（6）: 733-765.

Weaver J T, Le Quang B V, Fischer G. 1985. A comparison of analytic and numerical results for a two-dimensional control model in electromagnetic induction-I. B-polarization calculations.Gephys.J.Int., 82（2）: 263-277.

Weaver J T, Le Quang B V, Fischer G. 1986. A comparison of analytical and numerical results for a 2-D control model in .

J.Int., 87（3）: 917-948.

Wei W B, Jin S, Ye G F, et al. 2006. Conductivity structure of crust and upper mantle beneath the northern Tibetan Plateau: Results of super-wide band magnetotelluric sounding.ChineseJ.Geophys. （in Chinese）, 49（4）: 1215-1225.

Wei W B, Jin S, Ye G F, et al. 2010. On the conductive structure of Chinese continental lithosphere-experiment on “standard monitoring network” of continental EM parameters （SinoProbe-01）.ActaGeologicaSinica（in Chinese）, 84（6）: 788-800.

Weiss C J, Newman G A. 2002. Electromagnetic induction in a fully 3-D anisotropic earth.Geophysics, 67（4）: 1104-1114.

Xu S Z, Zhao S K. 1985. The topographic effects on magnetotelluric response.NorthwesternSeismologicalJournal（in Chinese）, 7（4）: 69-78.

Xu S Z, Wang Q Y, Wang J. 1992. Modeling 2-D terrain effect on MT by the boundary element method.ActaGeophysicaSinica（in Chinese）, 35（3）: 380-388.

Xu S Z, Li Y G, Liu B. 1997. The method and efficiency of 2-D terrain correction for Hx-polarization of MT.ChineseJ.Geophys. （in Chinese）, 40（6）: 842-846.

Zhao G Z, Tang J, Zhan Y, et al. 2004. Study on the crustal electrical structure and block deformation in Northeast margin of Tibetan Plateau.ScienceinChina（Ser.D） （in Chinese）, 34（10）: 908-918.

附中文参考文献

胡祥云, 霍光谱, 高锐等. 2013. 大地电磁各向异性二维模拟及实例分析. 地球物理学报, 56（12）: 4268-4277, doi: 10.6038/cjg20131229.

霍光谱, 胡祥云, 方慧等. 2012. 层状各向异性介质大地电磁联合反演研究. 物理学报, 61（12）: 129101.

霍光谱, 胡祥云, 刘敏. 2011. 各向异性介质中大地电磁正演研究综述. 地球物理学进展, 26（6）: 1976-1982, doi: 10.3969/j.issn.1004-2903.2011.06.011.

魏文博, 金胜, 叶高峰等. 2006. 藏北高原地壳及上地幔导电性结构-超宽频带大地电磁测深研究结果. 地球物理学报, 49（4）: 1215-1225.

魏文博, 金胜, 叶高峰等. 2010. 中国大陆岩石圈导电性结构研究-大陆电磁参数“标准网”试验（SinoProbe-01）. 地质学报, 84（6）: 788-800.

徐世浙, 赵生凯. 1985. 地形对大地电磁勘探的影响. 西北地震学报, 7（4）: 69-78.

徐世浙, 王庆乙, 王军. 1992. 用边界单元法模拟二维地形对大地电磁场的影响. 地球物理学报, 35（3）: 380-388.

徐世浙, 李予国, 刘斌. 1997. 大地电磁Hx型波二维地形改正的方法与效果. 地球物理学报, 40（6）: 842-846.

赵国泽, 汤吉, 詹艳等. 2004. 青藏高原东北缘地壳电性结构和地块变形关系的研究. 中国科学（D辑）, 34（10）: 908-918.

（本文编辑 汪海英）

MT modeling for two-dimensional anisotropic conductivity structure with topography and examples of comparative analyses

HUO Guang-Pu1, HU Xiang-Yun2*, HUANG Yi-Fan2, HAN Bo2

1HenanInstituteofGeologicalSurvey,Zhengzhou450000,China2HubeiSubsurfaceMulti-scaleImagingKeyLaboratory,InstituteofGeophysicsandGeomatics,ChinaUniversityofGeosciences,Wuhan430074,China

The magnetotelluric （MT） method is a technique for probing electrical conductivity structure of the Earth, from the near-surface down to upper mantle. MT observations can be significantly influenced by topography. Most of existing forward modeling algorithms for 2-D MT with topography are based on the electrical isotropic theory. However, it has been well established that the electrical anisotropy is widely present in the earth interior. Since the electrical anisotropy in crust and upper mantle is the connection between electric structure and the geological structure, it is vital to account for anisotropic effects while modeling 2-D MT fields with topography.We present a 2-D MT modeling approach for anisotropic media with topography. The solution is based on a finite-difference （FD） discretization of the second-order Maxwell′s equation in terms of electric fields parallel to strike for TE mode and magnetic fields parallel to strike for TM mode. The topography effect is accounted by changing the way in which the sampling electromagnetic fields are arranged, i.e. the sampling fields are in a staggered-order to make sure that the maximum bandwidth of the FD coefficient matrix can be assigned moderate memory spaces. The Weaver's approach is used to deal with the problem that the conductivity of some area on the air-earth interface may vary in the direction perpendicular to strike. Once the primary field is solved, the dual field can be obtained very easily by applying a discrete differential operator to the primary field.Two synthetic models for modeling the horst structure and graben structure, respectively, are used to evaluate the effects of topography. The responses of models with topography are compared with that of models without topography. It is found that most significant differences occur in the regions with sharp topography boundaries, such as the boundary between the foot of a mountain and the mountainslope and the boundary between the mountain slope and the mountaintop, while the minor differences appear in flat regions. The topography has much greater impact on the TM responses than TE responses, either for the horst model or the graben model. The graben structure can have more effects than the horst structure. The sharper the slope, the greater the influence is.We demonstrate the validity of the algorithm for 2-D MT anisotropic forward modeling with topography by numerical tests on both synthetic models and real datasets. The effects of topography are assessed by analyzing the modeling results of the horst model and the graben model. Finally, the impacts of the topography on MT responses for both large scale models and sharp slope models are evaluated by fitting the observed MT data through forward modeling.

Topography； Electrical anisotropy; Magnetotelluric; Finite difference; Tensor conductivity

国家深部探测专项第3项目（SinoProbe-03）,“十二五”国家科技支撑计划课题（2011BAB04B01），国家重点基础研究发展计划（2013CB733203）和国家自然科学基金（41274077、41474055）联合资助.

霍光谱,男,1983年生,博士,主要从事电磁法正反演研究.E-mail：huoguangpu1218@163.com

*通讯作者 胡祥云,男,1966年生,教授,博士生导师,主要从事电磁法理论应用教学与研究工作.E-mail:xyhu@cug.edu.cn

10.6038/cjg20151230.

10.6038/cjg20151230

P631

2015-05-15，2015-10-30收修定稿

霍光谱, 胡祥云, 黄一凡等. 2015. 带地形的大地电磁各向异性二维模拟及实例对比分析.地球物理学报，58（12）:4696-4708,

Huo G P, Hu X Y, Huang Y F, et al. 2015. MT modeling for two-dimensional anisotropic conductivity structure with topography and examples of comparative analyses.ChineseJ.Geophys. （in Chinese）,58（12）:4696-4708,doi:10.6038/cjg20151230.