对称非各向同性体大地电磁法二维异常特征

图1 层状岩石的电阻率分布模型

a—实际情况下的岩石;b—等效模型

Fig.1 Resistivity model of the layered rock

a—rock in the actual situation;b—equivalent model

表1 部分岩石和地层的非各向同性系数范围

Table 1 Anisotropy coefficient range of some rocks and strata

类型	λ范围
冲积层	1.02~1.10
干燥页岩,固结页岩	1.10~1.60
煤—无烟煤(沥青质的)	2.00~2.55(1.70~2.60)
辉长岩	1.10~2.00
花岗岩	1.05~1.50
石墨板岩	2.00~2.80
硬石膏和页岩互层	4.00~7.50
页岩、砂岩互层	1.05~1.15
石灰岩	1.00~1.14
含磁铁矿角页岩	1.20
磁铁矿	1.65
页岩层	1.02~1.05
砂岩层	1.10~1.60
绿泥石片岩—黑云母	1.12
板岩	1.10~2.25
火山岩	1.10~1.20

新窗口打开| 下载CSV

2 对称非各向同性体二维有限单元法

2.1 基本方程

在大地电磁测深法当中,通常假设磁导率和介电常数为真空当中的数值,即μ≈μ₀和ε≈ε₀,所使用的时间因子为e^-i^ωt,则相应的电磁场方程可以表示为^[15]

(2)

在式(2)当中,当地下介质为各向同性,其电导率为标量,而当地下介质为非各向同性,其电导率为张量。

图2a所示为在TM模式下的研究区域,定义x轴为走向方向,在二维各向同性介质中放置一个非各向同性体,该非各向同性体的电性主轴与测量轴之间存在倾角α,此时与异常体层面方向平行的电导率为σ_t=σ_y',与层面方向垂直的电导率为σ_n=σ_z'。

图2

图2 研究区域和单元分析

a—研究区域;b—母单元和子单元

Fig.2 Survey region and element analysis

a—survey region;b—parent element and subunit

在图2a所示的电性主轴坐标系x'y'z'中,电导率可以表示为:

(3)

σ' = [\begin{array}{l} σ_{x'} & 0 & 0 \\ 0 & σ_{y'} & 0 \\ 0 & 0 & σ_{z'} \end{array}],

在地面坐标系xyz中,电导率可以表示为式(4):

(4)

\begin{array}{l} σ = Rσ' R^{T} = [\begin{array}{l} 1 & 0 & 0 \\ 0 & \cos α & - \sin α \\ 0 & \sin α & \cos α \end{array}] [\begin{array}{l} σ_{x'} & 0 & 0 \\ 0 & σ_{y'} & 0 \\ 0 & 0 & σ_{z'} \end{array}] \cdot \\ [\begin{array}{l} 1 & 0 & 0 \\ 0 & \cos α & \sin α \\ 0 & - \sin α & \cos α \end{array}], \end{array}

其中R为旋转矩阵,用于连接坐标系xyz和x'y'z'。将式(4)展开则有:

(5)

σ - i ω ε_{0} = [\begin{array}{l} σ_{x'} - i ω ε_{0} & 0 & 0 \\ 0 & σ_{y'} co s^{2} α + σ_{z'} si n^{2} α - i ω ε_{0} & \frac{1}{2} (σ_{y'} - σ_{z'}) \sin 2 α \\ 0 & \frac{1}{2} (σ_{y'} - σ_{z'}) \sin 2 α & σ_{y'} si n^{2} α + σ_{z'} co s^{2} α - i ω ε_{0} \end{array}] 。

因x轴为走向方向,因此对x方向上的偏导数为0,将电磁场方程(2)进一步展开,可以得到TE和TM两种极化模式,其表达式如下:

TE极化模式(E_x型):

(6)

\{\begin{array}{l} \frac{\partial H_{z}}{\partial y} - \frac{\partial H_{y}}{\partial z} = (σ_{x'} - i ω ε_{0}) E_{x} \\ \frac{\partial E_{x}}{\partial z} = i ω μ_{0} H_{y} \\ - \frac{\partial E_{x}}{\partial y} = i ω μ_{0} H_{z} \end{array}

TM极化模式(H_x型):

(7)

\{\begin{array}{l} \frac{\partial E_{z}}{\partial y} - \frac{\partial E_{y}}{\partial z} = i ω μ_{0} H_{x} \\ \frac{\partial H_{x}}{\partial z} = (σ_{y'} co s^{2} α + σ_{z'} si n^{2} α - i ω ε_{0}) E_{y} \\ + \frac{1}{2} (σ_{y'} - σ_{z'}) \sin 2 α E_{z} \\ - \frac{\partial H_{x}}{\partial y} = \frac{1}{2} (σ_{y'} - σ_{z'}) \sin 2 α E_{y} \\ + (σ_{y'} si n^{2} α + σ_{z'} co s^{2} α - i ω ε_{0}) E_{z} \end{array}

通过式(6)可以看出,在TE极化模式下,非各向同性等同于电导率为σ_x'的各向同性,可以用各向同性来研究。从式(7)可以看出,在TM极化模式下,非各向同性跟电性主轴y'和z'方向上的电导率以及倾角有关系,不同于各向同性,需单独考虑。

2.2 变分问题

通过电磁场方程的展开式可以得到相应的变分问题为:

(8)

内边界条件为自然边界条件,不会对边值问题产生影响,因此不在考虑范围内。对于式(8):Ω为研究区域,Γ为区域的外边界,▽为二维哈密顿算符。在TM极化模式情况下,上边界AB为地表,不包含空气层,此时u=H_x,并且在地表的场值为常量1,τ可以表示为:

(9)

\begin{array}{l} τ = [\begin{array}{l} τ_{11} & τ_{12} \\ τ_{21} & τ_{22} \end{array}], \\ \{\begin{array}{l} τ_{11} = \frac{1}{σ_{y'} - i ω ε_{0}} si n^{2} α + \frac{1}{σ_{z'} - i ω ε_{0}} co s^{2} α \\ τ_{22} = \frac{1}{σ_{y'} - i ω ε_{0}} co s^{2} α + \frac{1}{σ_{z'} - i ω ε_{0}} si n^{2} α \\ τ_{12} = τ_{21} = - \frac{1}{2} (\frac{1}{σ_{y'} - i ω ε_{0}} - \frac{1}{σ_{z'} - i ω ε_{0}}) \sin 2 α \end{array} \end{array}

其中:

$λ = i ω μ_{0}, k = \sqrt[]{- i ω μ_{0} σ_{y'}}, ω = 2 π f 。$

下边界CD为离地表面足够远、左边界DA和右边界BC均为离异常体足够地远。

2.3 有限单元法

采用有限单元法作为数值模拟方法,求解步骤如下^[15]:

2.3.1 网格剖分

在TM极化模式下,研究区域为地表以下,将研究区域剖分为规则的矩形网格,在异常体周围和接近地表的地方适当加密网格,按照有限单元法的剖分原则来进行剖分。

2.3.2 单元分析

采用图2a中的自然坐标系来研究,采用图2b当中的母单元和子单元来剖分。节点编号如子单元中所示,将子单元映射到母单元(为等参单元)上,用等参变量(ξ,η)来描述,每一个等参变量的定义域为[-1,1],设置子单元的中点坐标为(y₀,z₀),单元宽为a,高为b。则有母单元和子单元之间的坐标转换公式为:

(10)

\{\begin{array}{l} y = y_{0} + \frac{a}{2} ξ \\ z = z_{0} + \frac{b}{2} η \end{array}

其微分关系为:

(11)

\{\begin{array}{l} d y = \frac{a}{2} ξ d ξ \\ d z = \frac{b}{2} η d η \\ d y d z = \frac{ab}{4} d ξ d η \end{array}

通过母单元可以构造如下形函数:

(12)

\{\begin{array}{l} N_{1} = \frac{1}{4} (1 - ξ) (1 + η) \\ N_{2} = \frac{1}{4} (1 - ξ) (1 - η) \\ N_{3} = \frac{1}{4} (1 + ξ) (1 - η) \\ N_{4} = \frac{1}{4} (1 + ξ) (1 + η) \end{array}

则单元场值可以表示为:

(13)

\begin{array}{r} u = N_{1} u_{1} + N_{2} u_{2} + N_{3} u_{3} + N_{4} u_{4} = \overset{4}{\sum_{i = 1}} N_{i} u_{i} = N^{T} u_{e} = {u^{T}}_{e} N, \end{array}

其中u_i(i=1,…,4)是矩形单元四个顶点的待定场值,u_e=(u₁,…,u₄)^T为矩形单元的场值列向量,N_e=(N₁,…,N₄)^T是形函数向量。将式(8)离散化为研究区域内所有单元的组合,则有:

(14)

其中单元积分有:

其中,K₁_e、K₂_e和K₃_e为4×4的矩阵,均为对称矩阵。

2.3.3 系数总体合成和求变分

将上面得到的单元积分扩展成全体节点,将各个单元的扩展矩阵相加,则可以得到如下公式:

(15)

F (u) = \frac{1}{2} u^{T} (K_{1} - K_{2} + K_{3}) u = \frac{1}{2} u^{T} Ku,

其中:u为全体节点的场值组成的列向量,K₁、K₂和K₃为各个单元矩阵相应的单元系数扩展叠加后得到的总体系数矩阵,K=K₁-K₂+K₃为正定、对称和稀疏矩阵。对式(15)求变分可以得到:δF(u)=δu^TKu=0,由于δu的任意性,则有:Ku=0,将AB上的边界条件代入上式,用乔列斯基分解法解线性方程组后,得各节点的场值,即u=H_x(TM模式)值。

2.4 辅助场、视电阻率和阻抗相位的计算

在TM极化模式下,通过解方程得到的H_x为主场值,依据电磁场方程展开式可以得到相应地面的辅助场值为

$E_{y} = \frac{1}{σ_{y'} - i ω ε_{0}} \frac{\partial H_{x}}{\partial z} 。$

通过主场值和辅助场值可以得到阻抗张量的表达式,从而可以得到相应的视电阻率和阻抗相位,其表达式为^[15]:

$Z = \frac{E_{y}}{H_{x}}, ρ_{s} = \frac{1}{ω μ_{0}} | Z |^{2}, φ = ta n^{- 1} (\frac{Im (Z)}{Re (Z)}) 。$

3 模型响应异常特征分析

3.1 模型一

如图3所示,设置的模型为厚板模型,异常体的顶部埋深为1 000 m,异常体的尺寸宽为2 000 m,高为8 000 m;围岩为均匀半空间各向同性介质,电阻率为1 000 Ω·m;将研究区域划分为44×86,其中水平方向上的间隔为44,垂直方向上的间隔为86。

图3

图3 厚板模型

Fig.3 Thick plate model

3.1.1 计算精度检验

为验证程序的正确性,设置图3非各向同性体在三个电性主轴方向上的电阻率分别为:10、10、250 Ω·m,此时其非各向同性系数为λ=5,设置非各向同性层理倾角(以下简称倾角)为45°。分别用Pek和文章程序计算该模型,对比频率范围为从0.01 Hz到1 000 Hz,对比位置为异常体中心点,将两种方法计算出来的视电阻率和阻抗相位分别在同一图形中表示出来。如图4所示,其中实线表示为Pek有限差分计算结果,圆圈表示为文章所使用的有限单元计算结果。可以看出从低频到高频在TM模式下视电阻率和阻抗相位曲线都拟合得很好,验证了程序的正确性。

图4

图4 中心点不同频率对计算结果的影响

Fig.4 The influence of different frequencies of the center point on the calculation results

3.1.2 倾角对非各向同性体正演响应的影响

设置图3异常体在三个电性主轴方向上的电阻率分别为:10、10、490 Ω·m,此时的非各向同性系数λ=7,设置倾角分别为0°、30°、60°、90°。对比频率为0.1 Hz,对比范围为所有测点,将4个角度视电阻率和阻抗相位计算结果分别在同一图形中表示出来(图5)。

图5

图5 不同倾角对计算结果的影响

Fig.5 The influence of different dip angles on the calculation results

从图5可以看出,在异常体中心点上方,视电阻率数值以0°、30°、60°、90°的顺序逐渐增大,阻抗相位曲线以0°、30°、60°、90°的顺序进行变化,说明当非各向同性系数不变时,倾角的变化将对异常体中心点的视电阻率和阻抗相位呈现出一定规律的影响。

3.1.3 非各向同性系数对非各向同性体正演响应的影响

设置图3非各向同性层理倾角为45°,根据一般岩石和地层非各向同性系数的范围1~7.5选择对比的非各向同性系数分别为1、3、5、7。对比频率为0.1 Hz,对比范围为所有测点,将四个非各向同性系数的视电阻率和阻抗相位计算结果分别在同一图形中表示出来(图6)。可以看出,在异常体中心点上方,视电阻率曲线以7、5、3、1的顺序变化,阻抗相位曲线以1、3、5、7的顺序进行变化,说明当倾角不变时,随着非各向同性系数的增大,非各向同性效应更加明显。

图6

图6 不同非各向同性系数对计算结果的影响

Fig.6 The influence of different anisotropy coefficients on the calculation results

3.1.4 断面图特征

为了更加直观和形象地展示非各向同性系数和倾角的变化对非各向同性体正演响应的影响,设置图3非各向同性层理倾角分别为0°、30°、60°、90°,当异常体为各向同性时在三个电性主轴方向上的电阻率分别为:10、10、10 Ω·m,设置非各向同性系数分别为1、3、5、7,即其他电性主轴方向上的电阻率不变,z'方向上的电阻率分别为10、90、250、490 Ω·m。对比频率范围为0.1~1 000 Hz,对比范围为所有测点,将异常体在4个倾角和4个非各向同性系数的视电阻率和阻抗相位正演拟断面图分别展示在同一色棒下(图7)。

图7

图7 厚板不同倾角和不同非各向同性系数拟断面对比

Fig.7 Section comparison figures of different dip angles and anisotropy coefficients of thick plate model

可以看出,图7a和图7b的(1)号图为倾角为0°、非各向同性系数为1,此时为各向同性。对比随着倾角和非各向同性系数变化正演响应与(1)号图的变化规律。从(1)~(4)可以看出当非各向同性系数为1时,随着倾角的变化,正演响应不会发生变化;从(1)、(5)、(9)、(13)可以看出当倾角为0°时,正演响应呈现比较微弱的变化;从(5)~(8)、(9)~(12)、(13)~(16)三组可以看出,在每一个非各向同性系数不变的情况下,对于本算例中异常体所取角度范围,非各向同性层理倾角越大,非各向同性效应越强;从(2、6、10、14)、(3、7、11、15)、(4、8、12、16)三组可以看出,在每一个倾角不变的情况下,随着非各向同性系数的增大,非各向同性效应增强;当倾角为30°和60°,且非各向同性系数不为1时正演响应表现出不对称性。

3.2 模型二

如图8所示M形模型,异常体总跨度长为5 200 m,高为950 m,厚为200 m,顶部埋深为100 m,将异常体放置在围岩电阻率为100 Ω·m的各向同性介质当中。将研究区域剖分为120×66,水平方向上的间隔为120(包括左右扩展的各18个网格),垂直方向上的间隔为66。对比频率范围0.1~1 000 Hz,选择41个频点,当异常体为各向同性时在三个电性主轴方向上的电阻率分别为:1、1、1 Ω·m,设置非各向同性系数分别为1、3、5、7,即只改变z'方向上的电阻率,分别为1、9、25、49 Ω·m,设置非各向同性层理倾角,即左翼倾角分别为0°、30°、60°、90°,相对应的右翼倾角分别为180°、150°、120°、90°。通过同一色棒下的正演拟断面图说明该异常体响应特征随不同倾角和非各向同性系数的变化规律。

图8

图8 M形异常体模型

Fig.8 M type abnormal body model

图9中倾角的变化是指左翼倾角的变化,同时伴随着右翼倾角的变化,标识上为左翼倾角变化。可以看出,图9a和图9b的(1)号图为倾角为0°、非各向同性系数为1,此时为各向同性。从(1)~(4)可以看出当非各向同性系数为1时,随着倾角的变化,正演响应不会发生变化,从视电阻率和阻抗相位图均可以模拟出M形异常体;从(5)~(8)、(9)~(12)、(13)~(16)三组可以看出,在每一个非各向同性系数不变的情况下,对于本算例中M形异常体所取角度范围,随着左翼非各向同性层理倾角增大,非各向同性效应增强;从(1、5、9、13)、(2、6、10、14)、(3、7、11、15)、(4、8、12、16)四组可以看出,在每一个倾角不变的情况下,随着非各向同性系数的增大,非各向同性效应增强;当倾角为30°和60°时正演响应没有表现出不对称性。

图9

图9 M形异常体不同倾角和不同非各向同性系数拟断面对比

Fig.9 Section comparison figures of different dip angles and anisotropy coefficients of M type abnormal body

4 结论与建议

使用有限单元的数值模拟方法,实现了对对称非各向同性体的大地电磁二维研究,正演模型依据一般岩石和矿石的电学性质选取非各向同性系数,对比分析了不同非各向同性系数和倾角对正演异常特征的影响。

1) 通过对比文章所采用的有限单元法与前人有限差分法的计算结果验证了正演程序的正确性。

2) 通过固定非各向同性系数,改变倾角的模型正演曲线对比,说明对于文章中垂直板状异常体所取角度范围,在异常体中心点处视电阻率值会随着非各向同性层理倾角的变大而增加。

3) 通过固定倾角,改变非各向同性系数的模型正演曲线对比,说明在异常体中心点处非各向同性效应随着非各向同性系数的变大而增强。

4) 通过简单的厚板模型和复杂的M形模型,对异常体选取不同的非各向同性系数和倾角,分别对比视电阻率和阻抗相位拟断面图。对比结果说明,当非各向同性系数为1时,正演响应不会随倾角而变;非各向同性效应在非各向同性系数不变时,对于文章中的垂直板状和M形异常体所取角度范围,会随非各向同性层理倾角增加而增强,在倾角不变时,会随非各向同性系数变大而更加凸显。这说明无论对于简单还是复杂的地质体,都会受到非各向同性的影响。

通过以上对称非各向同性体随着倾角和非各向同性系数变化的异常特征,可以看出不可以忽略非各向同性效应的存在,为更好地模拟地下地质情况,在进行正演模拟时需建立非各向同性模型。

(本文编辑:沈效群)

致谢:

感谢捷克学者Josef Pek对二维各向异性有限差分正演代码的开源提供。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

李金铭

. 地电场与电法勘探[M]. 北京: 地质出版社, 2005.

Li J

Electric and electrical prospecting[M]. Beijing: Geology Publishing House, 2005.

[2]

内吉 J G, 萨拉夫 P D . 大地介质电磁各向异性问题[M]. 邹永辉, 陈德志, 译. 北京: 地质出版社, 1992.

[本文引用: 2]

Negi J

, Saraf

P D

(translated by Zou Y H, Chen D Z). Anisotropy in geoelectromagnetism[M]. Beijing: Geology Publishing House, 1992.

[本文引用: 2]

[3]

O’Brueb D

, Morrison H

Electromagnetic fields in an N-layered anisotropic half-space

[J]. Geophysics, 1967,32(4):668-677.

[4]

Pek

, Santos

F A M

Magnetotelluric impedances and parametric sensitivities for 1-D anisotropic layered media

[J]. Computers & Geosciences, 2002,28:939-950.

[5]

霍光谱, 胡祥云, 方慧 , 等.

层状各向异性介质大地电磁联合反演研究

[J]. 物理学报, 2012, 61(12): 129101(1-10).

Huo G

, Hu X

, Fang

, et al.

Magnetotelluric joint inversion for anisotropic conductivities in layered media

[J]. Acta Physica Sinica, 2012, 61(12):129101(1-10).

[6]

Reddy I

, Rankin

Magnetotelluric response of latterlly inhomogeneous and anisotropic media

[J]. Geophysics, 1975,40(6):1035-1045.

[7]

徐世浙, 赵生凯 .

二维各向异性地电剖面的大地电磁场的有限元解法

[J]. 地震学报, 1985,7(1):80-90.

Xu S

, Zhao S

Solution of magnetotelluric field equations for a two-dimensional, anisotropic geoelectric section by the finite element method

[J]. Acta Seismologica Sinica, 1985,7(1):80-90.

[8]

Pek

, Verner

Finite-difference modeling of magnetotelluric fields in two-dimensional anisotropic media

[J]. Geophysics Journal International, 1997,128(1):505-521.

[9]

A finite-element algorithm for electromagnetic induction in two-dimensional anisotropic conductivity structures

[J]. Geophysics Journal International, 2002,148:389-401.

[10]

, Key

2D marine controlled-source electromagnetic modeling, part 1: an adaptive finite-element algorithm

[J]. Geophysics, 2007, 72(2): WA51-WA62.

[11]

, Pek

Adaptive finite element modelling of two-dimensional magnetotelluric fields in general anisotropic media

[J]. Geophysics Journal International, 2008,175:942-954.

DOI:10.3969/j.issn.1004-2903.2011.06.011 Magsci [本文引用: 1]

[12]

林长佑, 武玉霞, 杨长福 , 等.

水平层状对称各向异性介质的大地电磁资料反演

[J]. 地球物理学报, 1996,39(S1):342-348.

DOI: Magsci [本文引用: 1]

论述了某些条件下地球介质呈现出明显的电导率各向异性的情形，研讨了均匀水平层状对称各向异性介质大地电磁正、反演有关问题，试验表明所给出的反演方法是适用和有效的．给出了利用大地电磁资料了解地壳应力状态及研究膨胀各向异性地震前兆的一种途径．

Lin C

, Wu Y

, Yang C

, et al.

Magnetotelluric inversion for symmetrically anisotropic layered medium

[J]. Chinese Journal of Geophysics, 1996,39(S1):342-348.

Magsci [本文引用: 1]

[13]

霍光谱, 胡祥云, 刘敏 .

各向异性介质中大地电磁正演研究综述

[J]. 地球物理学进展, 2011,26(6):1976-1982.

在地球物理资料处理解释中,各向异性现象是一个无法忽视的重要因素.各向异性介质中大地电磁(MT)的正演数值模拟计算一直是国内外电磁感应领域研究的前沿课题.文中简要介绍了正演数值模拟计算的发展历程和几种主要的研究方法,探讨了未来各向异性介质中大地电磁正演研究方向、需要解决的难点和面临的挑战.

Huo G

, Hu X

, Liu

Review of the forward modeling of magnetotelluric in the anisotropy medium research

[J]. Progress in Geophysics, 2011,26(6):1976-1982.

Magsci [本文引用: 1]

[14]

胡祥云, 霍光谱, 高锐 , 等.

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

[J]. 地球物理学报, 2013,56(12):4268-4277.

DOI:10.6038/cjg20131229 Magsci [本文引用: 1]

经过半个多世纪的发展，国内外利用大地电磁法研究地球内部电性结构取得了令人瞩目的成就，这些研究成果多数是基于电性各向同性理论.然而地球内部普遍存在电性各向异性现象，地壳和上地幔中存在的电性各向异性是地电模型、地下结构和构造模型间一个重要的联系因素.本文首先由麦克斯韦方程出发，引入张量电导率，根据二维电性各向异性结构的特点，得到一组关于平行走向的电场分量Ex和磁场分量Hx的偏微分方程.使用有限差分法求解偏微分方程，求出Ex和Hx的近似解，并以此求得其它场分量；随后，通过对普通及特殊的二维电性各向异性结构做正演模拟，研究其对观测大地电磁场的影响，从而认识在普通及某种特定地质条件下的电磁传播特性，为其后对大地电磁实测资料的处理解释奠定理论基础；最后，以本文的研究成果为基础，将电性各向异性理论引入对实测大地电磁资料的处理解释中，通过对新疆某地的大地电磁资料做二维正演拟合解释，说明了电性各向异性现象的普遍存在，也验证了理论的正确性及算法的实用性，为今后分析解释大地电磁资料中的电性各向异性现象提供理论依据和技术指导，并开拓了对大地电磁实测资料处理的思路和方法.

Hu X

, Huo G

, Gao

, et al.

The magnetotelluric anisotropic two-dimensional simulation and case analysis

[J]. Chinese Journal of Geophysics, 2013,56(12):4268-4277.

Magsci [本文引用: 1]

[15]

徐世浙

. 地球物理中的有限单元法[M]. 北京: 科学出版社, 1994.

Xu S

The finite element method in geophysics[M]. Beijing: Science Press, 1994.