二维电导率各向异性介质电磁感应效应有限元算法

1、二维电导率各向异性介质电磁感应效应的有限元算法李予国 摘要 二维电性各向异性大地电磁场通过有限元方法计算获得。模型由包含各向异性块的层状介质构造组成。每个块或者层可能是由3×3的电导率张量限定。正演问题可以归结为两个耦合的关于垂直-平行于场分量的和。它们由有限元法数值地解出。线性有限元的解系是用预条件共轭梯度法求取的。之后，地表垂直-平行的场分量和，利用样条插值通过对和进行数值微分得到。 二维有限元算法通过对比二维有限差分的解证明是正确的。三个类型的模型的大地电磁响应用来证明各向异性的影响：水平，垂直，倾斜各向异性。第四个模型用来模拟在剪切和俯冲带各向异性的影响。这些模型的响应模拟了

2、电阻率曲线在长周期和具有明显的主对角线元素的张量主抗存在室的分离，正如以前观测到的。1、简介 近些年来，越来越多的关注投入到了对电性各向异性的电磁感应方面上的研究，尤其是为了试图完全理解对更长周期的大地电磁的观测数据。加拿大地区的大地电磁测量显示了电性的各向异性在下地壳和上地幔（Kellett等，1992；Mareschal等1995）。大的来自围绕德国深钻点（KTB）各向异性大地电磁曲线解释为上地壳高度的电各向异性（Eisel 和 Haak 1999）。Rasmussen (1988)使用一个深地壳层内的各向异性模型来解释瑞典南部的大地电磁横断面数据。层状构造的各向异性效应首先被OBrien

3、 和Morrison (1967)研究。Reddy 和 Rankin (1975)首先开始研究二维各向异性模型，他们只考虑水平各向异性的工作。最近，Osella和Martinelli (1993)计算了一些圆滑不规则边界有一定的主轴旋转角的大地电磁响应。Schmucker (1994)提出一个计算不均匀薄的覆盖体在层状半空间的电磁感应效应，其中可以含有一个或者两个各向异性的电导层。用有限元方法，Pek 和 Verner (1997) 还有 Weidelt (1996) 模拟了广义二维、三维各向异性构造，分别具有任意的主轴旋转角。 本文中，二维模型被重新研究，但是使用有限元（FE）的方法。首先，

4、我们详细描述有限元程序的数值近似。之后，我们展示有限元程序模拟不同的简单测试模型的响应。我们的结果和Pek& Verner (1997)的有限差分程序进行了对比。最后，我们计算三种类型的各向异性介质的大地电磁响应：水平，垂直，以及倾斜各向异性。我们最后以模拟一个地质剪切带或者板块俯冲带的地质背景下的模型结束。本文这样的构造是为了展示一个有限元策略去实现二维各向异性介质的电磁感应。尽管有限差分对于这类问题是有效的，有限元解仍然是需要的，因为有限元和有限差分都有自己的特殊的优势还可以彼此互相核对。另外，有限元方法可以模拟非矩形线性构造的真实地球构造。2、边界值问题 考虑到图1中的二维模型。

5、阴影异常区嵌入到一个简单的层状构造中，这个构造包括n层而且底部的层延伸到无穷远。为了简单，异常区域没有显示内边界作为一个定向独立的电导率，它至少和周边的第j层的电导率张量不同。然后我们的程序允许更一般的模型，异常区域可以细分成不同的不均匀块而且可以接触不同的层。 异常区域延伸到y正向到无穷远，例如模型可以出现到一个新的正常构造。在这种状态下算法排除源的感应。模型在走向方向x方向是不变的。感应电磁场在x方向也是不变的，尽管如此场的矢量也有三个分量对于极化的磁源平行或者垂直走向，由于各向异性。因此，通常区分TE和TM模式的各向同行构造就变得不合理了。由此就有结果，当主磁场矢量垂直轴向，电场将有主场

6、方向的分量，导致电荷产生于界面（除了倾斜各向异性的例子）。 假设一个时谐变量，在准静态近似的电磁场的情况下 (1) 其中是自由空间的磁导率，且 为电性电导率张量。它是对称的，当旋转到主轴方向（x',y',z'）上去，可以表示为 图1.本文中考虑的二维各向异性模型在特别的二维例子中，方程（1）消退为常数电导率的均匀块 很明显如果垂直-平行分量和得到了，剩下的分量和，和可由和的空间导数得到。方程（2）-（7）可以得到两个关于分量和的二阶偏微分方程： 从方程（8）和（9）很明显看出各向异性和独立的垂直-平行分量和通过一阶偏微分方程耦合在一起。因此，在这种情况下不存在独立的TE

7、、TM模式对于广义各向异性情况。因此，这些方程必须同时利用和求解。 之后叙述了不同的空间形式的各向异性的模拟结果。如图1一个均匀的异常区域具有各向同性的正常结构的围岩，一个或者两个平行于参考主轴（x，y，z）坐标系的主轴被考虑了进去。2.1 水平各向异性 当，主轴Z'是垂直的，另外两个主轴x'和y'在水平面（x,y）具有相对于x轴的垂直角。方程（8）和（9）就消退为 感应方程变得比方程（8）和（9）更加简单，但是垂直-平行分量和仍然耦合在一起通过关于z的一阶偏导。所以方程（10）和（11）必须联立在一起去为了得到和。2.2 倾斜各向异性 当，电导率张量主轴x'是

8、水平的在走向方向，另外两个主轴y'和z'在垂直面（y,z）具有相对于y轴的倾斜角。现在方程（8）和（9）就解耦为两个独立的模式方程（12）对应于E极化模式可以通过二维各向同行程序解得-只需要电导率替换为作为标量电导率。然而，TM模式的表达式子仍然比较复杂。2.3 垂直各向异性 如果，电导率张量的三个主轴真好和坐标系（x，y，z）是一致的。如果再有，这种结果关于z轴轴对称，水平方向的电导率完全不同于垂直方向的电导率。这样（8）和（9）仍然是耦合在一起的，我们有 尽管B极化模式感应方程比方程（9）更加简单，它不能通过二维各向同性程序解得，因为电导率在水平方向和垂直方向是不同的。如果

9、两个方向的电导率相同，方程（15）消退为各向同性的例子。之后的边界条件应用：模型的外边界条件，设定为狄里克雷边界条件，由对应的一维层状模型的左右边界条件。在模型的顶部和底部由左右边界的值进行线性插值得到（Pek 和 Verner 1997）。在内边界上，切向电场和磁场的分量，和必须连续。根据图1，和为 重构方程式子之间的关系至关重要。设n为不均匀区域的外侧的单位矢量，且有 和来自方程（8）和（9），和为y轴和z轴的单位向量相应的。这样，实质性的偏导数后，表达式为和为发现的切向分量。方程（17）是自动满足的，对自后的偏导是很重要的。结合这些表达式基本的积分-微分方程（23）和（24）可以得到浓缩

10、的形式，这些形式可以直接正演数值模拟对待。更多细节看李（2000）。3 有限元方法 由方程（8）和（9）组成的数值近视问题，基于有限元方法。垂直-水平分量和跟电导构造相切，所以处处连续。于是，在有限元中假设所有分量穿过单元都是满足连续。这种近似建立在模型域，整天包围二维不均匀体，区域延伸的足够远使得异常场在边界处很小。为了避免绝缘的空气层在地面上产生奇异性，我们假设空气层具有很小的电导率，但是非零，实际上小于在我们的计算中。数值计算结果显示电导率到。 加权余量方法用来获得从偏微分方程（8）和（9）得到的积分方程。方程(8)乘以一个关于电场的任意小的的变量，在模型域里的积分： 上述方程中，第一项

11、包含二阶偏导数，可以用格林函数简化。 其中指示模型域的边界。之后方程（19）可以写成等价形式 类似地，方程（9）乘以一个磁场分量任意小的变量,之后用高斯公式改变一下 加上方程（16）这样得到积分方程 这里我们仍然用这个公式模型域可被分割为矩形或者是三角形单元。矩形单元已被李（2000）年描述过。以下的部分，展示三角形单元。公式（20）和（21）分解到每个单元带有指示标e=1,2,ne其中表示为三角形单元e的边界。在形成方程（23）的时候，我们利用方程（18）时替换被积函数为.用一种简单的方法方程（24）由方程（17）获得。 现在边界条件必须被利用。在内部边界切向电场和磁场（）是被要求的。因为每

12、个边界在集成的过程中被遍历两次且方向相反，则在单元边界的线积分的总和为零。在外边界，因为设置的狄利克雷边界条件，电场和磁场的变量和等于0.因此鲜鸡粉也等于0.这样以来，方程（23）和（24）最终消退为根据有限元的线性近似方法，我们假设在每个三角形单元内，和都为关于y,z的线性函数：其中和为第i个顶点的电场和磁场对应的坐标为（），i=1,2,3为三角形（图2），为线性基函数。它们由一下定义： 其中 方程（25）和（26）在一个单元内的面积分用方程（27）-（32）估算。所有单元的积分都可以组装成两个线性方程系。详细内容见附录A。联立这些方程系，我们可以最终写下方程以矩阵的形式： 其中 和是阶方阵

13、（是结点的数目），和不对称的方阵，满足因此矩阵K是一个2的对称方阵。它是稀疏的复元素。U是2列向量，它是包含电场和磁场的未知数。方程（33）代入外边界条件，在内部的点上和可以用共轭梯度方法解出来。图2三角形单元切面 共轭梯度法最早由Hestens 和Stiefel (1952)提出，只适合求实数对称的正定的矩阵。Jacobs (1981)提出一种B共轭梯度法使得共轭梯度法应用至复矩阵非正定系。广义的共轭梯度法通过对系数矩阵k的预处理得到改善。在我们的例子中简单的雅克比扩展是很好的预处理。本文中我们估算的模型，网格包含2500多个点，方程的数目大约5000.我们总是用0向量U=0作为初始值，我们

14、总是采用趋于一致的解。 在共轭梯度方法中，仅稀疏对称的系数矩阵对角线上和下的元素是需要的，它是储存到一维数组中根据Schwarz (1991)的方法。利用这种共轭梯度法求解，计算内存大大降低。3 有限元方法 求解线性方程组（23），我们获得和在每个节点上的值。另外的场分量和很容易从方程（3）和（4）中获得，而分量和是根据（6）、（7）中一些额外的计算求导获得的。总结为 在我们的算法中，需要的导数是用样条插值数值计算得到的。 阻抗张量元素可以用两个正交的线性源极化模式的电场和磁场计算得来。例如主磁场方向在（极化模式一）。或者在y方向上（极化模式2）.主场为z<0处源的感应场和正常构造在处的

15、诱导场的和。 我们得到阻抗张量元素 电阻率和阻抗相位为 5 数值测试 为了测试电性各向异性有限元的适应性，我们的计算结果对于两种模式和已经获得的有限差分方法进行了比较。图3所示第一个测试模型首先由Reddy &和Rankin (1975)提出。我们计算这个模型在10s时候的响应，且和Pek 和Verner (1997)的计算结果对比。结果展示在图4中，有限元和有限差分吻合很好。图3水平各向异性岩脉模型Reddy &和Rankin (1975) 图4图三中的视电阻率（上部）和相位（下部）曲线。方块为Pek 和 Verner (1997)的有限差分结果，实线，为本文描素的有限元结果

16、。图5 有限元和有限差分测试模型的对比，模型包括一个出露地表的岩石块，覆盖在一个各向异性层上面。图6 图5中的模型的视电阻率（上部）和相位（下部）图。方块为Pek 和 Verner (1997)的有限差分结果，实线，为本文描素的有限元结果。 图5展示第二个模型有Pek 和 Verner (1997)年提出。一个水平各向异性层位于一个二维厨楼地表的各向异性场下。所涉及到的两个构造的走向彼此平行，但是不平行于二维模型的走向。模型可以用来展示MT数据的扰动，可以用一个复杂的各向异性例子计算得到。在图6中，有限元计算结果对于30s和Pek & Verner (1997)的对比。电阻率的相对误差

17、低于0.5%，阻抗的差异不超过1°。6 各向异性效应 已经注意横向非均匀各向同性地球模型不足以构成大地电磁的阻抗张量，特别是各向 异性阻抗保持或多或少的不变在延伸区域。这在长周期已经很好的研究，最好在一天之内的变化的周期，曲线最大值和最小值能差两个数量级或更多。的最小值在10（Schmucker 1998）。在结晶岩的表面，类似的观测已经到短周期内，我们的模型已应用到。响应的在x'主轴电阻率我们选择100-500，x'主轴电阻,10,z'主轴电阻率100-500. 在这部分里，一个简单的岩脉模型，在一个各向同性的半空间内，用来展示水平，垂直和倾斜各向异性的大地

18、电磁响应。所有的计算都在t=10s，除非有另外的叙述。对于均匀半空间10s的屈服深度为50km，5km对于10在y的主轴方向。值得一提的是维度的板状体，低边界在9km比半空间的屈服深度小很多。6.1 倾斜各向异性图7 一个具有倾斜异常的岩脉在各向同性的均匀半空间。岩脉的主轴电阻为，且具有变化的倾向角。图8 图7中10s时不同的变化角视电阻率图图9 显示图7中的电场在yz平面内的分布。产生的主场为，计算的周期为10s。 图7展示主轴面向切面的岩脉。各向异性岩脉的电阻率分别为。图8显示视电阻率对于变化的角在10s。这张图表明： （1）视电阻率独立于倾角，大地电磁场仅单向取决于。在这种情况下，且不受

19、各向异性的影响。 （2）电阻率受的影响。曲线不再对称关于模型的中心。模型的最小值偏离中心，偏向的一侧取决于倾角的一侧。偏离会变大随着水平或者垂直方向的角度偏差变大。我们解释这种现象通过看电场的分布。图9展示电场的分布在yz平面对于三个倾角（0°,30°,60°,90°）在极化模式1，产生主场的方向是。箭头的起点是有限元点。箭头长度和方向代表节点处场的大小和方向。图很明显地显示场的大小和方向变化随着各向异性岩脉倾角变化。（3）=0°时候，对应于一个y方向10，x和z方向500的模型响应。类似的=90°，对于500在x，y方向，z方向10

20、.视电阻率和相位对应于变化的角（0°,30°,60°,90°）在中心c（y=0）在图（a）和（b）分别展示。视电阻率在非常低的周期接近均匀本空间的真实电阻率，相位接近45°。然而，随着周期增长，和分离很明显，和同样，偏离大小取决于倾角。然而和因不同而不同，和曲线在长周期也适合均匀半空间的电阻率，相位和近似45°，但仍然受各向异性岩脉的影响。关于c点对称的a点和b点的和在10（c）和10（d）上分别显示。从这些图中可以看出来，视电阻率和相位曲线在a点，和不同于b点。然而当=0或90°时候，却是相同的。图10 图7中A、B、C地

21、表点的视电阻率和相位曲线。6.2 水平各向异性 图11 一个2-D水平各向异性岩脉在均匀各向同性围岩中。岩脉的主轴电阻为，且具有变化的各向异性走向角。图12 大地电磁主轴具有变化的走向角。叉表示值大地电磁张量非对角元素和，棒代表主轴元素的值和，参考Swift-rotated坐标系（x''，y''）。图13电场实部Re（Ex）和Re（Ey），在模型11的表面。产生的主场为，计算的周期为10s。 图11显示了主轴方向，它们现在水平方向倾斜角。不均匀的主轴电阻率分别。图12显示大地电场阻抗用西门子的表述方法沿着一个剖面表面，对于变化的各向异性角在10s。从图12我们可

22、以得到： （1）各向异性块上。最大和最小轴显示最大和最小的电导率。 （2）张量点到率主对角线元素和存在（除了30°和90°），且随着角增加而增加。 （3）远离不均匀体，异常场消退，非对角线元素变的相等和，对角元素消失：。当时。当从各向异性场通过各向同性场的时候，阻抗会在过渡区有一定的减少，这里的方向可能会严重扭曲，特别是涉及到各向异性体和浅部异常体。图13显示电场实部Re（Ex）和Re（Ey），在剖面的表面具有不同的角（=0°，30°，45°，60°，90°）在第二种极化模式，产生的磁场。可以从图中看出随着各向异性角的改变。

23、场的大小和方向都在改变，y方向电场，Re（Ey）存在=0°或=90°；当在两个角之间时Re（Ey）变大。6.3 倾斜岩脉模型 图14 图7中的一个二维模型=90°，岩脉旋转了45°模拟方向独立的剪切和俯冲带，具有不同的平行和垂直电阻率。 图14显示了一个图7中的二维模型=90°，但是岩脉旋转了45°，模拟方向独立的剪切和俯冲带。由于剪切作用，岩脉在剪切面和垂向有不同的电导率。如果没有不同存在，假设整个岩脉作为一个整体比围岩的导电性好，一个关于中心不对称的最小值在走向方向被观测到。如果平行电导率10倍阻性围岩。对称性得到了保证单在岩脉上

24、最小值不可见。如果再剪切方向更加导电，最小值解更加显著，正如期待的（图15a），关于分歧。图15（b），1s一下开始，在周期更长时编程一个静的偏移。图15视电阻率在岩脉的剖面上对于10s图14中的模型（a）,测量点位于y=0（b）。电阻率垂直于岩脉边缘（各向同性例子），或者十倍小于。视电阻率曲线对于在10-1000之间变化，曲线紧挨着。7 结论 我们展示了一个二维大地电磁数值迷你程序。其主要特点如下：（1）每个各向异性块的电导率张量都由一个3乘3的矩阵构成，这样保证任意的各向异性构造都能被考虑，包括空间中的水平、垂直、倾斜的各向异性。（2）数值近似问题由有限元方法得到，这种方法很适合具有坡度的

25、边界模型和地形，如图14中所举的例子。（3）假设一个非常小，但正的空气介质，空气层可以和导电模型成为一个整体。两种模式的方程在整个模型区域均匀近似，这样就简单一些相对于区分空气和大地或者仅大地变量如Pek 和 Verner(1997)，它们把空气的电导率设为0.线性有限元方程利用了快速迭代的技术，增加的空气层中的变量Hx不是过多的问题。数值试验显示场分量是稳定的，实际上不熟宽的空气电导率的影响限制，特殊的从10-12开始，这对偃师的电导率已经足够低，已经低于大地上面的空气电导率。 （4）方程（33）是利用改变的共轭梯度法求的，是一个复的对称矩阵。我们发现简单的雅克比扩展是一种有效的预处理方法。

26、 （5）垂直走向场Ey和Hy用通过样条插样数值得到。致谢 作者对U. Schmucker博士和J. Pek博士对本文的讨论和改进表示感谢。J. Pek博士也提供了MT一维各向异性构造的正演程序和两个二维的有限差分的结果。同时也感谢Martyn Unsworth, Chester Weiss和另外一位匿名的编辑的建设性意见。作者感谢德国交流服务提供的博士奖学金。本文由格延根大学地球物理许愿管理。参考文献 Eisel, M. & Haak, V., 1999. Macro-anisotropy of the electrical conductivity of the crust: a m

27、agneto-telluric study from the German Continental Deep Drilling site (KTB), Geophys. J. Int., 136, 109122. Hestens, M.R. & Stiefel, E., 1952. Method of conjugate gradients for solving linear systems, J. Res. Nat. Bur. Stand., 49, 409436. Jacobs, D.A.H., 1981. The exploitation of sparsity by iter

28、ative methods,in Sparse Matrics and their Uses, pp. 191222, ed. Duff, I.S.,Springer, Berlin. Kellett, R.L., Mareschal, M. & Kurtz, R.D., 1992. A model of lower crustal electrical anisotropiy for the Pontiac Subprovince of the Canadian Shield, Geophys. J. Int., 111, 141150. Li, Y., 2000. Finite e

29、lement modeling of electromagnetic elds in two- and three-dimensional anisotropic conductivity structures, PhD thesis, University of Go ¨ ttingen. Mareschal, M., Kellett, R.L., Kurtz, R.D., Ludden, J.N., Ji, S.,& Bailey, R.C., 1995. Archaean cratonic roots, mantle shear zones and deep elect

30、rical anisotropy, Nature, 375, 6527, 134137. OBrien, D.P. & Morrison, H.F., 1967. Electromagnetic elds in an Nlayer anisotropic half space, Geophysics, 32, 668677. Osella, A.M. & Martinelli, P., 1993. Magnetotelluric response of anisotropic 2-D structures, Geophys. J. Int., 115, 819828. Pek,

31、 J. & Verner, T., 1997. Finite-difference modelling of magneto-telluric elds in two-dimensional anisotropic media, Geophys. J. Int.,128, 505521. Rasmussen, T.M., 1988. Magnetotellurics in southwestern Sweden;evidence for electrical anisotropy in the lower crust?, J. geophys. Res.,93, 78977907. Reddy, I.K. & Rankin, D., 1975. Magnetotelluric response of laterally inhomogeneous and anisotropic media, Geophysics, 40, 10351045. Schwarz, H.R., 199