MT阻抗张量在二维与三维介质中的分析

胡玮哲， 熊高君

(成都理工大学 地球物理学院 成都 610059 )

0 引言

大地电磁测深法(Magnetotelluric ,MT)是20世纪50年代初，由前苏联学者A.T.Tikhonov[1]与法国地球物理学家L.cagnird[2]分别提出，以天然交变电磁场为场源来获取有关地下地质体的电性结构特征的信息，属于频率域测深的一种方法。通过对二维以及三维的阻抗张量的分析来表征地球系统的传输性质以及描述其地质意义。传统的处理方法是将Z(阻抗张量)进行旋转分析，一般可以满足对二维构造问题的研究。但三维的分析只能将其近似的表现为二维介质来研究，所以Z的旋转分析效果是有限的。在阻抗张量的分析方法中，Swift[3]提出了对阻抗张量的传统分解方法,通过对观测坐标的旋转，使得阻抗张量的对角最小或者反对角元素最大，由此求得构造主轴方位以及阻抗大小；Yee[4-5]利用数学上的奇异值分解方法推导出阻抗张量Z的正则分解形式，既结合了传统旋转分析的基础又拓展了对阻抗张量地质意义的认知；国内对大地电磁测深法的真正应用起始于20世纪70年代，于20世纪初王光锷[6]阐述了阻抗张量的正则分解与传统分析方法的等价性条件。由于电磁场于三维介质传播的复杂性，目前国内、外对三维介质阻抗张量的研究正在进行。笔者通过分析二维与三维介质的阻抗张量，来阐述对阻抗张量的意义。

1 二维介质阻抗张量的分析

1.1 Z的旋转分析

(1)

旋转后所得的阻抗张量：

Z(θ)=[u]·Z·[u]T

(2)

即可得

(3)

图1 坐标变换示意图Fig.1 Coordinate transformation

由式(3)可以令：

Zxx(θ)=Z2-Z1cos 2θ+Z4sin 2θ

Zxy(θ)=Z3-Z1sin 2θ+Z4cos 2θ

Zyx(θ)=-Z3+Z1sin 2θ+Z4cos 2θ

Zyy(θ)=Z2+Z1cos 2θ-Z4sin 2θ

(4)

其中

(5)

若经过旋转角θ0观测坐标系与二维介质电性主轴重合，则

Zxx(θ0)=0,Zxy(θ0)=ZTE

Zyx(θ0)=-ZTM,Zyy(θ0)=0

(6)

但在实际观测中Zxx(θ0)、Zyy(θ0)不可能为“0”，所以可构造函数F(θ)使得

(7)

所以在二维介质中，通过对观测坐标系中的阻抗张量进行旋转分析可以获得旋转角度θ0，从而可以确定二维介质的电性主轴所在的正交方位。

1.2 的SVD分解

由Z代数中矩阵的奇异值分解原理可知[8]，对阻抗张量可进行如下分解：

Z=U·∧·V+(+表示转置共轭)

称为Z的奇异值，且λ1≥λ2>0

对于任意的平面单色电磁波，取角频率ω，波数为k，沿Z轴传播，则在地面x-y平面内，一般产生椭圆极化。以电场矢量为例，数学表达式为式(8)。

(8)

其中：A为波的振幅；a为共同相位；φ为相对相位；

(9)

为了更方便地确定电磁波的极化状态，我们引入复极化率P:

(10)

由式(9)可以令

a1=cosθ，a2=sinθ,θ∈[0,π]

所以相应的电场向量与磁场向量对应形式为：

(11)

(12)

即

(13)

(14)

我们可以清晰观察出[U]、[V]分别描述了电场与磁场的极化状态，并且有两个相互正交的主极化状态。

依据上述的分解性质，将共同相位引入∧中并与[U]、[V]结合:

(15)

2 三维介质的阻抗张量的分析

2.1 三维介质阻抗张量的通用分析

(16)

其中：α1、α2、β1、β2、γ1为对应分量的系数；它们可以通过电阻率ρ(Ω·m)与空间坐标(x，y，z)计算得出。

(17)

其中:α3、α4、β3、β4、γ2为对应分量的系数，它们可以通过电阻率ρ(Ω·m)与空间坐标(x，y，z)计算得出。将大地电磁场源分离成两个相互正交的线性极化平面波, 场源Mx、场源My存在如下形式：

(18)

其中电场分量在水平方向的分量示意如图2 所示。

图2 正交场源作用下电场分量在水平方向 的分量示意图Fig.2 Component of the electric field component in the horizontal direction under the action of the orthogonal field source

根据场的叠加原理，对于一般垂直入射平面电磁波，各分量存在如下关系：

(19)

或表示为：

(20)

因为E=Z·H

(21)

(22)

即

(23)

联立式(14)、式(15)解方程组(23)，可得阻抗张量元素解

(24)

(25)

(26)

(27)

式(24)～式(27)为计算阻抗张量各元素的通用公式。 因此三维介质条件下，将场源M分解为电场形式或磁场形式代入即可获得其值。综上对于张量阻抗的计算，是不同场源在水平面上的分量叠加，将不同的场源分量带入通用公式中可以计算出张量阻抗元素的值。

2.2 三维介质阻抗张量扩展

在三维介质的一般情况下，电场矢量存在:

(28)

(29)

(30)

k2=-iωμδ-ω2με

(31)

由式(28)～式(31)可知在三维条件下，电磁场的6个分量Ex、Ey、Ez、Hx、Hy、Hz全部耦合，不能解耦。

现推导出三维介质张量阻抗的形式，由方程可知，在三维介质中6个电磁场分量是相互耦合的，因此，提出如下三阶阻抗张量形式：

(32)

由对张量阻抗的通用分析可知：

(33)

(34)

(35)

在实际的操作过程中可以获得Ex、Ey、Hx、Hy、Hz这5个电磁场分量，结合式(35)计算出的电场垂直分量Ez，则电磁场的6个分量都成为已知。因此,可由电磁场分量求出阻抗张量各元素，由此可计算出相应的视电阻率以及相位特征。

Ex=ZxxHx+ZxyHy+ZxzHz

(36)

讨论阻抗张量的计算方法。

若对某一频率的电磁场作3次独立观测，则应有

(37)

式(37)中，上标表示对电磁场分量的观测次数。

写作矩阵形式：

(38)

式(38)可写为

(39)

同理存在

因地下地质体的相应是非线性的，所以对式(39)有

因为|H～|≠0

其中*表示伴随矩阵。

(40)

同理

(41)

(42)

综上由式(40)～式(42)，将阻抗张量写作：

(43)

2.3 三维介质阻抗张量算法

为求阻抗张量计算方法，可采用算法(代数重构)进行反演计算。通过对给定初值进行计算，再将计算的结果与已知的数据求差，将差值以对应系数所占的权重进行分配，以上一次迭代结果相加，得到下一次的迭代值。

当满足给定误差时，停止迭代，输出阻抗元素。获得在三维介质中的阻抗张量元素值，由此给出了在无噪的情况下三维介质阻抗张量计算的具体方法。

但实际上不能直接这样处理，因为实际观测的数据存在噪音，仅从3组观测数据来确定阻抗张量元素是不严谨的。为尽量减小噪音对数据质量的损坏，可以采用对大量观测数据求平均值的方法，一般根据最小二乘原理求取阻抗张量元素的最佳估计值。

令

Exc=ZxxHx+ZxyHy+ZxzHz

(44)

其中:Hx、Hy、Hz为实测值；Exc为理论计算值，因观测误差与噪音的综合影响，Exc并不等于实际观测值，因此可以定义均方差函数

(45)

式中：*表示共轭复数；N为观测次数(一般取3)。

将式(44)代入式(45)中得：

(46)

为使Ψ→min，则有

(47)

因为式(46)中Zxx、Zxy、Zxz均为复数，所以

iImZxy)Hyi+(ReZxz+(ReZxz+iImZxz)Hzi]}·

(48)

以Zxy为例，对其实部与虚部求偏导数，则有

(49)

(50)

(51)

若用<>符号表示功率谱的平均值，则式(51)可表示为:

(52)

若对Zxx实部与虚部求偏导数，则有

(53)

若对Zxz实部与虚部求偏导数，则有

(54)

将(52)～(54)写作矩阵形式

(55)

也可利用ART算法可以直接反演出阻抗张量。

同理可分析

Ey=ZyxHx+ZyyHy+ZyzHz

Ez=ZzxHx+ZzyHy+ZzzHz

的情况。

以上通过对已观测的数据进行最小二乘处理以及求出观测数据的功率谱平均值，可以得到如式(55)的矩阵形式，采用算法可以求解出阻抗张量元素的值。

3 结论

1)大地电磁测深法利用天然交变的平面电磁波为场源，依据在地表观测到的包含地下地质体信息的电场值与磁场值计算出阻抗张量，并由阻抗张量计算出的视电阻率以及相位曲线特征来分析大地电性结构特征。 在 MT张量阻抗的分析中，对严格的二维介质可以进行旋转分析以及正则分解，若要描述复杂的三维介质必须采用二维近似以及引入二维偏离度，其中正则分解可以描述主阻抗与电场与磁场的关系，可以确定地球的传输性质，较传统的旋转分析有了更明显的地质含义。但自然界地下地质体是以三维形式出现，因此二维近似会产生偏差。而阻抗张量的通用分析则通过对场源的分解计算阻抗张量元素的值，可以直接对三维介质进行表述。

2)对于张量阻抗的分析中，对严格的二维介质可以进行旋转分析以及正则分解，若要描述复杂的三维介质必须采用二维近似以及引入二维偏离度，其中正则分解可以描述主阻抗与电场与磁场的关系，可以确定地球的传输性质，较传统的旋转分析有了更明显的地质含义。但自然界地下地质体是以三维形式出现，因此二维近似会产生偏差。而阻抗张量的通用分析则通过对场源的分解计算阻抗张量元素的值，可以直接对三维介质进行表述。

3)对三维介质阻抗张量的形式进行了更进一步的描述，采用了三阶张量的形式表述阻抗，使得在形式上进行大地电磁测深时获得的电场与磁场的关系统一，可以计算出三维介质的阻抗张量。

4)通过三维介质阻抗张量扩展中的分解方法，可以获得阻抗张量的计算公式，对阻抗张量元素进行分部处理，结合观测数据的功率谱平均值计算出功率谱矩阵，利用算法反演出对应的阻抗张元素值。

参考文献：

[1] А.Н.Тихонв.Об определиии Электрических Характерстих Глубокихслоев Земной Коры，Докл.АН[J].СССР73，1950：295-297.

[2] CAGNIARD L．Basic theory of the magnetotelluric method of geophysical prospecting[J]．Geophysics，1953，18(3):605-635．

[3] SWIFT C W. A magnetotelluric investigation of an electrical conductivity in the south western united states[D].Cambridge:University of Cambridge，1967.

[4] E.YEE,K .V .PAULSON.The canonical decomposition and its relationship to the other forms of magnetotelluric impedance tensor analysis[J].Journal of Geophysics ，1987， 61(3):173-189.

[5] E.YEE,K.V .PAULSON.Canonical decomposition of the telluric transfer tensor[J].Journal of Geophysics，1987，61(3):190-199.

[6] 李伟.一种改进的阻抗张量分解方法及其应用[D].长沙：中南大学，2010.

LI W. An improved impedance tensor decomposition method and its application [D]. Changsha: Central South University, 2010. (In Chinese)

[7] 陈乐寿，王光鄂.大地电磁测深法[M].北京：地质出版社，1990.

CHEN L S, WANG G E. Magnetotelluric sounding method[M]. Beijing: Geological Publishing House, 1990. (In Chinese)

[8] 晋光文，孙杰，王继军.大地电磁(MT)阻抗张量的正则分解及其初步应用[J].地震地质，1998,20(3)：243-249.

JIN G W, SUN J, WANG J J. Magnetotelluric (MT) impedance tensor canonical decomposition and its preliminary application[J]. Seismology and geology, 1998,20 (3): 243-249. (In Chinese)

[9] 尹兵祥，王光鄂.大地电磁的阻抗张量正则分解的地质含义分析[J].地球物理学报,2001(3)：279-288.

YIN B X, WANG G E. Geological implication analysis of impedance tensor decomposition of magnetotelluric[J].Chinese Journal of Geophysics,2001(3):279-288.(In Chinese)

[10] 谭捍东，魏文博，邓明,等. 大地电磁的张量阻抗通用计算公式[J].石油地球物理勘探,2004，39(1)：113-116.

TAN H D, WEI W B, WANG L, et al. Universal calculation formula of tensor impedance of magnetotelluric petroleum geophysical exploration[J].Oil geophysical prospecting,2004, 39 (1): 113-116.(In Chinese)