2012, Vol. 55
常规大地电磁(Magnetotelluric, MT)测深资料的反演解释主要是基于各向同性介质进行的.近年来，随着地球物理观测技术和理论方法的进展，对地球的认知正逐渐深入，地球各向异性问题越来越引起人们的关注.越来越多的证据[1-6]也表明:地壳和上地幔通常是各向异性的，地球介质的物理性质(如弹性､导电性､磁性､导热性和密度等)存在明显的各向异性，地球内部介质的各向异性是地球科学的难点与前沿课题.MT 法与其他方法(如地震)相比，其成本低､探测深度大(从几百米到几百公里深)，是地壳上地幔探测的主要方法之一.MT 响应函数对不同方向的导电性很敏感，也即对地球介质的电导率各向异性很敏感.通过对地表观测的MT 响应函数进行各向异性反演，就能很容易地求出不同方向的电导率，从而揭示出电性各向异性.由此就可以根据不同方向的电导率或者导电性，推导出与各向异性密切相关的岩石内部介质裂隙优势取向､矿物晶体排列和应力场及形变带特征等信息.此外，利用MT法研究各向异性无需特殊的观测手段.因此利用MT 法进行地球各向异性研究具有其独特的优势.
已有的研究表明，利用传统的各向同性大地电磁反演方法去解释本质上是各向异性的MT 资料，结果势必会带来很大的误差，甚至使解释结果陷入歧途[7].理论上，对某一给定频率，一维各向异性介质会产生全部的4 个MT 阻抗张量响应函数(即Zxx､Zxy､Zyx､Zzz或它们相应的视电阻率及相位)，各向同性介质只能产生一种有意义的视电阻率及相位(即ρxx=ρyy=0，ρxy=ρyx，φxy=φyx)，因此原则上无法用一维各向同性反演方法反演各向异性地电结构.在实际情况下，当各向异性不明显时，人们经常不考虑各向异性，而用各向同性反演得到近似的结果.但是，如果在地电结构是各向异性(即各向异性很明显，不能忽略)的情况下，那么利用各向同性反演只能得到各向同性地电结构的反演结果，这样势必带来很大的误差或错误的解释.而且，现实工作中也经常发现一些用各向同性反演无法拟合观测数据的异常现象，而利用各向异性反演则可以对这些数据做出合理的解释.因此，利用MT 各向异性反演方法去解释大地电磁资料是很有必要且非常有意义的探索.
但是，长期以来MT 各向异性反演一直进展缓慢.目前国内外有关MT 各向异性反演研究的文献报道并不多见.Abramovici等最早尝试以各向异性模型进行了MT 反演，获得了关于各向异性地球的电性参数的阻抗张量导数值[8].Abramovici 和Shoham 于1977年基于广义反演技术开发了一维各向异性电导率MT 反演程序[9].随后，Regis 和Rijo利用相同的反演程序论证了结合来自于测井､地球物理及地质解释的先验信息反演电磁数据对稳定各向异性模型至关重要[10-11].Li等发展了一种用于解决一维CSAMT 各向异性方位角问题的高效非线性反演算法[12].Yin研究了一维层状各向异性模型MT 反演的内在非唯一性[13].Pek 等把标准Occam 一维反演推广到层状各向异性导体的大地电磁反演中[14].国内林长佑等[15]开发了一维层状对称各向异性介质MT 反演方法，但其反演中舍弃了辅助阻抗而仅仅拟合两主阻抗张量元素观测资料.杨长福等[16]分析了某些条件下电导率各向异性行为，并应用林长佑等[15]的一维各向异性MT 反演方法，对甘肃及其邻近地区实测MT 资料进行解释，但其也未考虑辅助阻抗.
事实上，与主阻抗张量一样，辅助阻抗同样蕴含着重要的信息，甚至在某些情况下(如辅助阻抗较大时)，仅考虑主阻抗元素无法较好地还原出目标体，这就需要同时利用主阻抗与辅助阻抗元素进行反演以更好地还原出目标体.鉴于此，本文尝试MT 全张量响应(即所有的阻抗张量的视电阻率和相位)的一维各向异性反演研究，同时考虑主阻抗与辅助阻抗，在一维各向异性反演中同时拟合4 对视电阻率及阻抗相位，开发了相应的算法并应用于一维层状各向异性介质模型的MT 反演，最后把这种方法应用于宁夏固原附近的MT 实测资料的反演，并进行了部分有关一维各向异性反演的讨论.
2 反演方法MT 测深可以获得4对张量阻抗视电阻率及阻抗相位资料，即8个张量，而在各向同性介质一维反演中仅仅考虑1对视电阻率及阻抗相位，即使用2n个数据(设n为频点个数)，对于L(一般取3~4)层各向同性介质，需反演求解的模型参数个数为2L-1，而本文中同时考虑8个张量进行各向异性介质一维反演，需反演求解的模型参数个数为4L-1，是各向同性介质情形下的近2 倍，因而反演问题也变得更加复杂.但利用的数据量为各向同性情形下的4倍(即8n个数据)，因此也就相对降低了解的非唯一性.基于此，非常有必要进行MT 资料的各向异性反演.由于各向异性反演问题非常复杂，为便于实现，本文主要研究相对简单的层状各向异性介质一维反演问题，并采用较为成熟的广义逆矩阵反演法[17-19].
对于L层水平层状各向异性介质(图 1)，需反演求解的模型参数个数为4L-1，是L层各向同性介质的近2倍，要得到更可靠的反演结果，就必然要求更加充分的资料信息量，这也正是本文进行全张量响应一维各向异性介质反演研究的必要.关于广义逆矩阵反演法的详细过程，请参考有关文献[17-19]，本文不再赘述，下面仅给出反演的基本过程.
|
图 1 层状各向异性介质模型几何参数及坐标系 x, y, z以为测量坐标系，xi'(i=1，2，3)为第i层的各向异性主轴，似θi(i=1，2，3)为第i层的各向异性角，ρip(i=1，2，3 p=1，2) 为第i层各向异性主轴 > 上的电阻率. Fig. 1 The geometry and parameters of the layered anisotropic model and coordinate system used x, y, z is themeasured coordinate system, xi'(i= 1，2，3) is the anisotropic principal axis of the zth layer, while θi(i= 1，2，3) is the anisotropic angle and ρip(i= 1,2,3； p=1,2) is the resistivity on the pth anisotropic principal axis. |
对于本文所用的反演方法，观测资料可以写成如式(1)所示的向量形式:
|
(1) |
其中，n为观测频点数；ρaikl与φaikl(i=1，…，n；k，l=x，y)分别代表视电阻率和阻抗相位的观测值，上角标T 代表矩阵的转置.
在MT 法中，由于视电阻率(Ωm)与阻抗相位(rad)具有不同的量纲，因此在进行视电阻率-阻抗相位的联合反演时，必须首先通过某种方法(如单位换算､对数化､加权等)使它们处于相同的量级上，从而使它们在目标函数中的地位相同，以便在反演过程中对模型参数改正量的贡献相当.晋光文[20]在进行各向同性介质的视电阻率-阻抗相位资料联合解释时，采用一种表征视电阻率及阻抗相位在反演过程中相对重要性的权因子，把两种观测资料写入到同一目标函数之中.本文直接对观测视电阻率取常用对数，对相位根据实际情况(一般情况下加权因子β 可取1；当相位资料比较差时，可使β 取值稍微小些)加权，从而使得视电阻率值与阻抗相位在目标函数中的地位相当.相比各向同性介质，全张量阻抗的视电阻率-相位一维各向异性反演目标函数右端项增加近4倍(式(2)).
|
(2) |
式中，ρaikl与φaikl的含义与(1)式中相同；ρcikl与φcikl(i=1，…，n；k，l=x，y)分别代表理论模型的响应函数值；m= (m1，…，mm)T 表示模型参数向量，m为模型参数个数，在本文中所用的模型参数向量可具体表示为m= (ρ11，ρ12，θ1，h1；ρ21，ρ22，θ2，h2；ρ31，ρ32，θ3)T，Δm=m-m0 是模型参数的修改量；β为加权因子，根据实际情况对相位加权；α 为阻尼因子，表示对模型参数修改量Δm进行约束，Δm/mj为模型参数的相对变化量.
反演问题的求解即使目标函数取极小值，根据极值条件，即满足:
|
(3) |
由此可得正规方程:
|
(4) |
式中，α 为阻尼因子，I为单位矩阵，A= (A1kl ，A2kl)T(k，l=x，y)为偏导数矩阵，其中
|
|
Q= (Q1kl ，Q2kl)T(k，l=x，y)为观测值与计算值的差，其中
|
利用奇异值分解(SVD)法，求解方程(5)，可得Δm，进而可以获得求解地球电性模型参数的迭代公式:
|
(5) |
式中，mk为第k次迭代后的模型参数.反复求解(4)式，逐次修改模型参数，直到满足拟合标准(一般将各观测响应拟合到平均相对观测误差水平).
3 正演和偏导数矩阵计算通常情况下，水平层状对称各向异性介质任一层内的电磁场水平分量可由下一层内的电磁场水平分量递推求取.对于下半空间标号为0 的N+1 层介质(地表层标号为N)，地球表面阻抗张量可用(6)式求取(顶层主坐标方向下)，此处仅列出计算结果，有关详细步骤请参考相关文献[15， 21]:
|
(6) |
式中，Fn和Gn为地球表面处的二阶传递矩阵，可由半空间向上逐层递推求取；
|
m4L为模型参数向量m中的第4L(L为模型层数)个元素，u(m)为与频率有关的变量.
由此可以求得视电阻率及阻抗相位:
|
式中Im(Z(m))与Re(Z(m))分别为阻抗张量Z(m)的虚部和实部.
进而可由(7)式求出偏导数矩阵中的元素:
|
(7) |
实际计算偏导数值时，可采用近似数值计算方法按(7)式计算.
4 理论模型反演试验为了检验上文所述反演方法的正确性，首先把一个一维各向异性介质模型退化为一维各向同性介质模型，并用本方法对其进行反演计算，然后用现有的一维各向同性反演方法对同一模型进行反演计算，经比较，两者的结果基本一致，初步验证了本方法的正确性.
为进一步验证其效果，我们利用多种水平层状各向异性介质模型对上述MT 各向异性反演方法进行了检验，限于篇幅，本文仅列出三层理论模型的反演结果.该模型每一层的模型参数用2 个主方向电阻率､1个各向异性角(各向异性主轴xi′与测量坐标x的夹角)与层厚度表示(图 1).图 2为该模型的MT 函数理论响应值与反演计算值的拟合结果，从图中可以看到，当初始模型计算值与理论模型响应值的差异较大时，最终反演模型计算值可以很好地拟合理论模型响应值.考虑到野外MT 观测资料通常受到各种噪声的干扰，为更加符合实际，在理论模型响应值中加入10%的随机噪声后，再次对算法进行检验，经过多次迭代反演计算，基本上仍然能够较好地拟合理论模型响应值(图 3).比较理论模型与反演结果模型(图 4)，可以发现，当观测值(即理论模型响应值)不含噪声时，即使初始模型与理论模型差异较大，经过适当次数的迭代，理论模型也可以很好地恢复出来，而且当在观测值中加入噪声时，基本上仍然可以有效地恢复出理论模型.由此可见，上文所述方法是有效可行的.
|
图 2 三层各向异性介质模型理论观测值(散点线)与初始模型响应值(虚线)及反演拟合结果(实线) 图中，GRkl与GPhskl(k,l= x,y)分别代表理论模型观测视电阻率及相位响应值，ORkl与OPhskl(k,l= x,y)分别代表初始模型视电阻率及相位响应值，NRkl与NPhskl(k,l= x,y)分别代表反演结果拟合视电阻率及相位响应值(下同). Fig. 2 The response for three-layer anisotropic structure (dot line) and the initial model (dashed line) and inversion fitting result (solid line) In the graph, GRkl； and GPhskl； (k,l= x,y) represent the observed apparent resistivity and phase responses to the theoretical model respectively, and ORkl and OPhskl(k,l= x,y) denote the responses to the initialmodel, while NRkl and NPhskl(k,l= x,y) mark theresponses to the rnversion results (hereinafter the same). |
|
图 3 三层各向异性介质模型理论观测值加人10%的随机噪声(散点线)与初始模型响应值(虚线)及相应的反演拟合结果(实线) Fig. 3 The response for three-layered anisotropic structure after adding 10% random noise to the data (dot line) and the initial model (dashed line) and inversion fitting result (solid line) |
|
图 4 一维各向异性模型(实线)与初始模型(虚点线)及反演结果虚线､点线分别为不含噪声与含噪声响应值的反演结果. Pp2分别为各向异性主轴方向1，的电阻率(下同). Fig. 4 Original one-dimensional structure (solid line) and the initial model (dash-dotted line) and results of the inversionThe dashed line is the result of the response exclusive of noise, while the dotted line is the result of the inversion after adding noise to the data. p1 and p2 represent the resistivity on the 1st and 2nd anisotropic principal axis respectively (hereinafter the same). |
为了检验所述方法的实用性，我们对宁夏固原附近的一些MT 测点资料进行反演计算，限于篇幅，此处仅列出一个测点的反演计算结果(图 5).资料处理过程中除剔除个别“飞点"外，未对原始资料作过多的修饰性处理.由图 5 可见，整体而言，反演计算的8个MT 响应函数基本上拟合了观测值，在某些频段反演计算的MT响应函数不能与观测值完全拟合，这可能由于观测资料中包含较大的噪声.总体上认为该方法是可以应用于实测资料处理的，具有一定的实用性.这为今后的进一步研究提供了基础，也为MT 资料解释增加了一种新的手段.图 6为采用本文所述方法对MT 实测资料的反演结果.本文的实例仅是为了说明上述方法是能够用于实测资料反演解释，有关该方法对实测资料的处理效果和地质解释以及与各向同性反演方法的异同已超出本文的讨论范围，作者将另行撰文阐述.
|
图 5 宁夏固原附近的一个MT测点资料的一维各向异性反演结果 Fig. 5 The one-dimensional anisotropic inversion results of the MT field data in the vicinity of Gu'yuan, Ningxia |
|
图 6 实测数据的全张量一维各向异性反演结果 Fig. 6 Models given by the whole impedance one-dimensional anisotropic inversion for the MT field data in the vicinity of Gu'yuan, Ningxia |
本文提出一种全张量视电阻率-相位的一维各向异性反演方法，并应用于理论模型响应值和实测数据的反演计算中.理论模型研究验证了所述方法的正确性和有效性，实例分析计算表明，这种方法是可以用于实测资料的处理解释，具有一定的实用性.
目前MT 一维反演在MT 资料的解释中仍然占有重要的地位，各向同性一维反演仅仅利用1 对视电阻率和阻抗相位，这无疑是对MT 资料众多信息的浪费.本文所述方法采用4对视电阻率和相位，利用的信息量更多，可以相对减小解的非唯一性.反演结果可以为地质构造解释提供更多的信息(如不同深度不同方向上的电阻率信息､各向异性角等)，为地质解释提供更客观的选择.但在进行实测资料反演时，必须首先对资料进行各向异性的分析和识别，如果资料是各向异性的，必须利用各向异性的方法才能得到更客观真实的结果.但是关于MT 资料的各向异性识别问题仍然是一个比较复杂的研究课题，作者将在今后进行相关方面的研究.
| [1] | Keen C, Tramontini C. A seismic refraction survey on the mid-Atlantic ridge. Geophys. J. R. Astr. Soc , 1970, 20(5): 473-491. DOI:10.1111/j.1365-246X.1970.tb06087.x |
| [2] | Minster J B, Jordan T H, Molnar P, et al. Numerical modelling of instantaneous plate tectonics. Geophys. J. R. Astr , 1974, 36(3): 541-576. DOI:10.1111/j.1365-246X.1974.tb00613.x |
| [3] | Christensen N I. The magnitude, symmetry and origin of upper mantle anisotropy based on fabric analyses of ultramafic tectonites. Geophys. J. R. Astr. Soc. , 1984, 76(1): 89-111. DOI:10.1111/j.1365-246X.1984.tb05025.x |
| [4] | Crampln S, Chesnokov E M, Hipkin R G. Seismic anisotropy-the state of the art II. Geophys. J. R. Astr. Soc. , 1984, 76(1): 1-16. DOI:10.1111/j.1365-246X.1984.tb05017.x |
| [5] | Hilley G E, Bürgmann R, Zhang P Z, et al. Bayesian inference of plastosphere viscosities near the Kunlun Fault, northern Tibet. Geophys. Res. Lett. , 2005, 32: L01302. DOI:10.1029/2004GL021658 |
| [6] | 金淑燕. 上地幔的岩石组构和各向异性. 地质科技情报 , 1993, 12(3): 32–38. Jin S Y. Rock fabric and anisotropy of the upper mantle. Geological Science and Technology Information (in Chinese) , 1993, 12(3): 32-38. |
| [7] | 杨长福. 大地电磁二维对称各向异性介质的有限元数值模拟. 西北地震学报 , 1997, 19(2): 27–33. Yang C F. MT numerical simulation of symmetrically 2D anisotropic media based on the finite element method. Northwestern Seismological Journal (in Chinese) , 1997, 19(2): 27-33. |
| [8] | Abramovici F, Landisman M, Shoham Y. Partial derivatives for the one-dimensional magnetotelluric problem. Geophys. J. Roy. Astr. Soc. , 1976, 44(2): 359-378. DOI:10.1111/j.1365-246X.1976.tb03661.x |
| [9] | Abramovici F, Shoham Y. Inversion of anisotropic magnetotelluric data. Geophys. J. R. Astr. Soc., , 1977, 50(1): 55-74. DOI:10.1111/j.1365-246X.1977.tb01324.x |
| [10] | Regis C R T, Rijo L. 1-D inversion of anisotropic magnetotelluric data // Extended Abstracts Book, 50th IC Soc. Bras. Geof. Ed., Brazil, vol. II. 1997: 673-674. |
| [11] | Regis C R T, Rijo L. Approximate equality constraints in the inversion of anisotropic MT data: Abstracts Book // 15th Workshop on Electromagnetic Induction in the Earth. Cabo Frio, Brazil, 2000: 47. |
| [12] | Li X B, Oskooi B, Pedersen L B. Inversion of controlled source tensor magnetotelluric data for a layered Earth with azimuthal anisotropy. Geophysics , 2000, 65(2): 452-464. DOI:10.1190/1.1444739 |
| [13] | Yin C C. Inherent nonuniqueness in magnetotelluric inversion for 1D anisotropic models. Geophysics , 2003, 68(1): 138-146. DOI:10.1190/1.1543201 |
| [14] | Pek J, Santos F A M. Magnetotelluric inversion for anisotropic conductivities in layered media. Physics of the Earth and Planetary Interiors , 2006, 158(2-4): 139-158. DOI:10.1016/j.pepi.2006.03.023 |
| [15] | 林长佑, 武玉霞, 杨长福, 等. 水平层状对称各向异性介质的大地电磁资料反演. 地球物理学报 , 1996, 39(Suppl): 326–332. Lin C Y, Wu Y X, Yang C F, et al. Magnetotelluric inversion for symmetrically anisotropic layered medium. Chinese J. Geophys. (in Chinese) , 1996, 39(Suppl): 326-332. |
| [16] | 杨长福, 林长佑. 甘肃及邻近地区的各向异性介质MT反演及该区地壳深部应力场和形变带. 内陆地震 , 1997, 11(2): 143–147. Yang C F, Lin C Y. MT inversion made for anisotropic medium in Gansu and its neighbouring areas and deep-crustal stress field and deformation belt in the same areas. Inland Earthquake (in Chinese) , 1997, 11(2): 143-147. |
| [17] | Backus G E, Gilbert J F. Numerical applications of a formalism for geophysical inverse problems. Geophys. J. R. Astr. Soc. , 1967, 13(1-3): 247-276. DOI:10.1111/gji.1967.13.issue-1-3 |
| [18] | Backus G E, Gilbert J F. The resolving power of gross Earth data. Geophys. J. R. Astr. Soc. , 1968, 16(2): 169-205. DOI:10.1111/gji.1968.16.issue-2 |
| [19] | Backus G E, Gilbert J F. Uniqueness in the inversion of inaccurate gross Earth data. Phil. Trans. Roy. Soc. London. Ser. A , 1970, 266(1173): 123-192. DOI:10.1098/rsta.1970.0005 |
| [20] | 晋光文. 大地电磁视电阻率和阻抗相位资料的联合解释. 地震地质 , 1987, 9(4): 81–86. Jin G W. Joint interpretation of the apparent resistivities and impedance phases from magnetotelluric data. Seismology and Geology (in Chinese) , 1987, 9(4): 81-86. |
| [21] | 陈乐寿, 王光鄂. 大地电磁测深法. 北京: 地质出版社, 1990 : 85 -109. Chen L S, Wang G E. Magnetotelluric Sounding Method (in Chinese). Beijing: Geological Publishing House, 1990 : 85 -109. |