地球物理学报  2011, Vol. 54 Issue (1): 245-256   PDF    
利用大地电磁三维反演方法获得二维剖面附近三维电阻率结构的可行性
林昌洪1,2,3, 谭捍东1,2,3, 佟拓1,2,3     
1. 中国地质大学地质过程与矿产资源国家重点实验室,北京 100083;
2. 中国地质大学地下信息探测技术与仪器教育部重点实验室,北京 100083;
3. 中国地质大学(北京)地球物理与信息技术学院,北京 100083
摘要: 大地电磁野外实测数据目前大多为二维剖面数据.如何反演这些二维剖面数据获得较为接近实际地电情况的结果,是多数大地电磁工作者关心的问题.我们通过对理论模型的三维响应进行分析和对合成数据及实测资料的反演结果进行对比研究,讨论了利用三维反演的方法来获得大地电磁二维剖面附近三维电阻率结构的可行性.结果表明:可用三维反演的方法来解释二维剖面数据.对大地电磁二维剖面的张量数据进行三维反演,不仅可以沿剖面获得较好的二维断面结果,还能够得到二维反演所不能获得的剖面附近的三维电阻率结构信息.合成数据的反演算例表明:对二维剖面数据进行三维反演时,对角元素对于圈定剖面附近三维异常体的空间分布具有独特作用,应尽量反演所有的张量元素.
关键词: 大地电磁反演      二维      三维      张量数据     
The possibility of obtaining nearby 3D resistivity structure from magnetotelluric 2D profile data using 3D inversion
LIN Chang-Hong1,2,3, TAN Han-Dong1,2,3, TONG Tuo1,2,3     
1. State Key Laboratory of Geological Processes and Mineral Resources, China University of Geosciences, Beijing 100083, China;
2. Key Laboratory of Geo-detection (China University of Geosciences), Ministry of Education, Beijing 100083, China;
3. School of Geophysics and Information Technology, China University of Geosciences, Beijing 100083, China
Abstract: Currently, most of MT(magnetotelluric) data are still collected on 2D profiles. The issue that most MT researchers concern is how to invert these 2D profile data to get better results, which are closer to the real structure. Based on the analysis of 3D tensor impedance response generated from the test models and on the study of the inversion results of synthetic data and field data, we discuss the possibility of obtaining nearby 3D resistivity structure from magnetotelluric 2D profile data using 3D inversion. The results show that it is possible to interpret 2D profile data using 3D inversion method. Not only a reasonable image beneath the profile but also reasonable pictures of nearby 3D structure which can not be got by 2D inversion can be obtained using all tensor elements of the 2D profile data in the 3D inversion. The synthetic examples show that the on-diagonal elements have special effect on recovering the distribution of 3D abnormity near the profile. Therefore, all the tensor elements are suggested to be used in 3D inversion to interpret the profile data.
Key words: Magnetotelluric inversion      2D      3D      Tensor elements     
1 引言

大地电磁测深法是电法勘探中发展比较快的一种方法,已被广泛应用于能源与资源勘探、工程与环境勘察、地震与火山灾害监测、地壳上地幔结构探测、海洋地球物理等领域[1, 2].使用大地电磁方法进行野外地质勘探,为了获得可靠的地下介质电性参数的分布,可在工区布置面积性数据采集工作,然后用三维反演技术对这些数据进行处理解释.大地电磁数据的三维采集和对三维数据进行三维反演也是电磁领域的前沿课题之一.然而,由于勘探工作成本和工区环境等因素,当前大多数实测数据仍然是二维剖面采集.

如何反演这些二维剖面数据获得较为接近实际情况的地电断面是多数大地电磁工作者关心的问题.近年来,随着计算机内存和速度的倍增,二维反演算法的研究已经趋于成熟,大地电磁实测数据的解释主要采用二维反演方法.使用二维反演方法的前提是假设地下电阻率结构是二维分布的,并且测线与二维电阻率结构走向垂直.但是,实际的地下地质情况比较复杂,地下电阻率结构可能是三维分布的,测线也可能与主体电阻率结构走向不垂直,若对这样的实测资料进行二维反演解释,很可能得不到可靠的地电模型,有时甚至会得到严重畸变的结果[3~7].因此,对大地电磁二维剖面数据的有效反演方法技术进行研究是很有必要的.那是否可以用三维反演方法[8~15]对二维剖面数据进行反演? 同时如何反演这些二维剖面数据获得一些三维电阻率结构分布情况也是大地电磁工作者关心的课题.

本文的目的是从理论模型合成数据和实测资料的反演算例出发,探讨利用三维反演的方法来获得大地电磁二维剖面附近三维电阻率结构的可行性和效果.文中,首先用三维数值模拟程序计算出理论模型产生的一些二维剖面张量数据,并分析这些二维剖面数据的响应特征和数据中是否含有三维电阻率结构的信息.接着,用二维共轭梯度反演程序对两种理论模型合成的二维剖面数据进行二维反演,再用大地电磁三维共轭梯度反演程序[8]对同样的二维剖面数据进行了三维反演.然后,我们给出实测二维剖面资料的二维反演和三维反演结果.通过对比这些二维剖面数据的二维反演结果和三维反演结果之间的差异,以及分析二维剖面数据的三维反演结果,获得了有实用价值的新认识,发展了大地电磁二维剖面数据的三维反演解释手段.

2 理论模型及其张量阻抗特征分析 2.1 模型一:单棱柱体模型

为了获得二维剖面合成数据,我们设计了一个三维单棱柱体模型(见图 1).电阻率为10Ωm, 大小为6km×6km×3km, 顶面埋深为3km 的棱柱体埋藏于电阻率为100Ωm 的地下均匀半空间.图 1中黑色虚线分别表示剖面A、B、C 的位置.三维网格剖分为38×38×35(含10 个空气层).用三维交错采样有限差分算法程序[16, 17]计算出单棱柱体模型在地表处产生的频率为3.3Hz、1Hz、0.33Hz和0.1 Hz的张量阻抗数据.

图 1 单棱柱体模型的平面视图 黑色虚线分别表示剖面A、B、C的位置.坐标原点位于棱柱体中心.剖面A 在Y方向上的位置为Y= -2.875km;剖面B 在Y方向上的位置为Y=-0.25km;剖面C在Y方向上的位置为Y=3.75km. Fig. 1 Plan view of the prism model The black dashed lines show the position of profile A,B,C.The coordinate origin is at the center of the prism.Profile A is at Y= -2.875 km along X direction ; Profile Bsat Y = -0.25 km along X direction; Profile C s at Y=3.75 km along X direction.

图 2给出了频率为0.1 Hz时单棱柱体模型在地表处产生的张量阻抗(ZxxZxyZyxZyy)的平面等值线图,该图显示了ZxxZxyZyxZyy响应的平面分布特征.从图中可以看出非对角元素ZxyZyx的异常最大值都位于棱柱体中心,并且异常值往两侧边缘方向逐渐减弱;而对角元素ZxxZyy的异常最大值则分布于棱柱体的四个角点处,棱柱体中心处的对角元素值最小(与背景值一致).对于位于棱柱体的中心剖面B,非对角元素含棱柱体的异常信息最强,而对角元素则没有棱柱体的异常.对于位于棱柱体边缘处的剖面A 和剖面C,对角元素含有很强棱柱体的异常,而非对角元素所含的棱柱体异常相对较弱甚至没有异常(例如,图中剖面C 的Zxy数据关于棱柱体的异常几乎与背景值一致).

图 2 频率为0.1 Hz时单棱柱体模型产生的ZxxZxyZyxZyy响应的平面等值线图黑色实线表示棱柱体的边界,黑色虚线分别表示剖面A、B、C的位置.图中,第一行显示实部;第二行显示虚部;从左到右分别对应ZxyZyxZxxZyy. Fig. 2 Impedance terms ZxxZxyZyx and Zyy generated fromthe prism model

The black solid lines show the prism margins, the black dashed lines show the position of profile A,B,C.The first row denotes the real part, and the second row denotes the imaginary part.The left to the right are the results for ZxxZxyZyx and Zyy.

2.2 模型二:三棱柱体模型

为了模拟复杂地质情况,我们设计了一个较为复杂的三棱柱体模型(见图 3),并且考虑了三个棱柱体对二维剖面数据的不同影响程度,让三个棱柱体分别位于剖面下方的不同位置.在深度0~10km以内,半空间在X方向上被分成三部分,电阻率分别为50Ωm、100Ωm、200Ωm, 其中100Ωm 部分的宽度为10km(-5km~5km),50Ωm 和200Ωm的两部分往两侧无限延伸,这三个部分在Y方向上都无限延伸.10km 以下是电阻率为100Ωm 的均匀半空间.电阻率为1Ωm、大小不同的三个棱柱体分别埋藏于上述电阻率不同的三部分围岩中.三个棱柱体的大小和位置等数据见表 1.图 3 中黑色虚线表示剖面A 的位置,剖面A 在X方向上的位置为Y=-0.1km.三维网格剖分为30×71×40(含10个空气层).用三维交错采样有限差分算法程序计算出剖面A 的频率为10Hz、3.3Hz、1Hz、0.33Hz和0.1Hz的张量阻抗数据.

表 1 三个棱柱体的大小和位置 Table 1 The dimension and position of the three prisms
图 3 三棱柱体模型示意图 黑色虚线表示剖面A 的位置.(a)沿Z方向俯视图;(b)沿Y负方向侧视图. Fig. 3 Plan view of the three prism model The black dashed line show the position of profile A.(a) See from negative to positive along Z direction;(b) See from positive to negative along Y direction.

图 4 给出了剖面A 上三棱柱体模型产生的ZxxZxyZyxZyy响应的拟断面图.图中,周期取对数,ZxxZxyZyxZyy取绝对值的对数.从Zxy的拟断面图分析可知,只有棱柱体ρ1 在剖面A 上有异常,而棱柱体ρ2ρ3 则在剖面上没有反映出异常信息.从Zyx的拟断面图分析可知,棱柱体ρ1ρ2ρ3在剖面上都有异常信息,并且ρ1 的异常最强,ρ3 的异常最弱.从ZxxZyy的拟断面图分析,棱柱体ρ1ρ2ρ3 在剖面上都反映出了异常信息,并且ρ3 的异常最强,ρ1 的异常最弱.

图 4 剖面A 上三棱柱体模型产生的ZxxZxyZyxZyy响应的拟断面图 第一行表示lg|Zxy|;第二行表示lg|Zyx|;第三行表示lg|Zxx|;第四行表示lg|Zyy|. Fig. 4 Impedance terms ,Zy,2y工 and 2yy generated from three prism model on profile A The first row shows lg |Zxy| ,the second row shows lg |Zyx|,the third row shows lg |Zxx|,the fourth rowshows lg |Zyy|.
3 理论模型合成数据的二维反演 3.1 单棱柱体模型合成剖面数据的二维反演

用二维共轭梯度反演程序对合成的剖面A、剖面B和剖面C数据进行二维反演,所有合成数据中都加入1%的高斯随机误差(下同).初始模型为100Ωm, 反演的结果如图 5所示.在剖面A 和剖面B 的反演结果中,TM 模式(对应Zxy数据)反演结果最为合理,棱柱体的大小和形状都还原得比较好.剖面C的TM 模式反演结果对棱柱体几乎没有反映,这验证了剖面C 的Zxy数据中不含棱柱体的异常信息(见图 2).剖面B和剖面C 的TE 模式(对应Zyx数据)和TM+TE 模式反演结果由于在低频段受到三维异常体的影响,棱柱体下方出现低阻异常.从这些二维反演结果可以看出:当地下异常体为三维结构时,通过二维反演方法仅能够获得所反演剖面下方异常体的分布情况(受三维异常体的影响,所获得的二维断面甚至可能产生多余的异常结构),不能获取异常体的平面分布范围.因此,我们考虑对这些二维剖面数据进行三维反演.

图 5 剖面A、B、C 数据的二维反演结果 第一列为TM 模式反演结果,第二列为TE 模式反演结果,第三列为TM+TE 模式联合反演结果. Fig. 5 Results from 2D inversion of profile A,B and C data The first column shows the results of TM mode, the second column shows the results of TE mode, the third coiumn shows the results of TM + TE joint mode.
3.2 三棱柱体模型合成剖面数据的二维反演

用二维共轭梯度反演程序对合成的剖面A 数据进行二维反演.初始模型为100Ωm, 反演的结果如图 6所示.相对而言,在三种二维反演模式的反演结果中,TM+TE 模式反演结果对三部分围岩以及10km 以下均匀半空间的反映最接近真实模型.对于深度0~10km 内50Ωm 和100Ωm 围岩,TM+ TE 模式和TM 模式的反演结果都较为合理,而TE模式的反演结果中两部分围岩的底边界都变浅并且在100Ωm 围岩中出现高阻假异常.对于深度0~10km内200Ωm 围岩,受低阻棱柱体ρ3 的影响,三种模式的反演结果都与真实模型有较大差别.对于10km 以下半空间,TM 模式和TE 模式反演结果分别出现低阻和高阻假异常.由于剖面A 位于低阻棱柱体ρ1 中心上方(见图 3),三种模式的反演结果都较合理的还原出ρ1.由于剖面A 位于低阻棱柱体ρ2 的边缘上方,TM+TE 模式和TM 模式的反演结果对ρ2 反映较弱,TE 模式的反演结果对ρ2 的反映较强.由于剖面A 位于低阻棱柱体ρ3 的侧上方并且与ρ3Y方向上具有一定距离,TM+TE 模式和 TM 模式的反演结果对ρ3 几乎没有反映,只有TE模式的反演结果对ρ3 有较弱的反映.

图 6 剖面A 数据的二维反演结果 左上图为真实模型,左下图为TM 模式反演结果,右上图为TM+TE 模式联合反演结果,右下图为TE 模式反演结果. Fig. 6 Results from 2D inversion of profile A data The top left figure shows the test model, the bottom left figure shows the resutt of TM mode, the top right figure shows the result ofTM+TE joint mode, the bottom right figure shows the result ofTE mode.
4 理论模型合成数据的三维反演 4.1 单棱柱体模型合成剖面数据的三维反演

对单一剖面而言,最典型的剖面是位于棱柱体中心上方的剖面B 和位于棱柱体边缘上方的剖面 A.初始模型为100Ωm, 在PC 机上用张量阻抗三维共轭梯度反演程序对剖面B 和剖面A 数据进行三维反演.PC 机的配置(下同)为:Intel(R)core(TM)2 Duo处理器,主频2.2GHz, 内存2.0G.为了分析对角元素的重要性,反演分四种情况进行(下同):(1)只反演单元素(Zxy),(2)只反演单元素(Zyx),(3)反演两个非对角元素(ZxyZyx),(4)反演四个张量元素(ZxxZxyZyxZyy).

4.1.1 剖面B的反演结果

图 7的第二列为Y=-0.25km 处沿X方向的垂直断面,即剖面B 下方的垂直断面.从该断面四种不同的反演结果看,无论是双元素反演还是四个元素反演得到的地电模型都比较合理地反映了真实模型的分布,而单元素反演的两种结果都对低阻棱柱体没有反映.与该剖面的二维TM 模式反演结果相比,双元素和四个元素的三维反演结果所获得的棱柱体电阻率值更接近模型真实值.

图 7 对剖面B数据进行三维反演的结果 图中第一行表示真实模型,第二行为单元素(Zxy)的反演结果,第三行为单元素(Zyx)的反演结果,第四行为双元素(ZxyZyx)的反演结果,第五行为四个元素(ZxxZxyZyxZyy)的反演结果.黑色虚线表示棱柱体的边界.第一列为深度4.25km 处的水平截面图,第二列为Y=-0.25km 处沿X方向的垂直断面图,第三列为Y=-2.875km 处沿X方向的垂直断面图,第四列为X=-0.25km 处沿Y方向的垂直断面图. Fig. 7 Results from the 3D inversion of 2D data on profile B The top row shows the test model, the second row shows the results from the inversion of only ,the third row shows the results from the inversion of only ,the fourth row shows the results from the inversion of both and ,and the bottom row shows the results from the inversion of all impedance tensor elements.The black dashed lines show the prism margins.The first column shows the horizontal slices at 4.25 km depth, the second column shows the vertical slices at Y= - 0.25 km along the X axis, the third column shows the vertical slices at Y=-2.875 km along the X axis, and the fourth column shows the vertical slices at X=-0.25 km along the Y axis.

图 7的第二行,单元素(Zxy)的反演结果对低阻棱柱体没有反映.图 7 的第三行,单元素(Zyx)的反演结果虽然对低阻棱柱体有反映,但出现了两个目标体:在X方向上,棱柱体的位置相对正确;在深度方面,棱柱体从地表开始就有反映.四个元素(ZxxZxyZyxZyy)同时反演的结果与双元素(ZxyZyx)的反演结果相似:在XYZ方向上,较合理地还原出棱柱体的形状(见图 7 的第四和第五行).分析四个元素反演的结果与双元素的反演结果相似的原因在于:剖面B 的对角元素数据不含棱柱体的异常信息(见图 2),因此对角元素数据对于还原棱柱体没有帮助.

4.1.2 剖面A 的反演结果

图 8的第三列为Y=-2.875km 处沿X方向的垂直断面,即剖面A 下方的垂直断面.从该断面四种不同的反演结果看,虽然低阻棱柱体的形状和电阻率值与真实模型相比有些变化,无论是双元素反演还是四个元素反演得到的地电模型都比较合理地反映了真实模型的分布,并且反演结果与二维 TM+TE 模式联合反演的结果类似.而单元素反演的两种结果都对低阻棱柱体没有反映.

图 8 对剖面A 数据进行三维反演的结果 图中第一行为单元素(Zxy)的反演结果,第二行为单元素(Zyx)的反演结果,第三行为双元素(ZxyZyx)的反演结果,第四行为四个元素(ZxxZxyZyxZyy)的反演结果.黑色虚线表示棱柱体的边界.第一列为深度4.25km 处的水平截面图,第二列为Y=-0.25km 处沿X方向的垂直断面图,第三列为Y=-2.875km 处沿X方向的垂直断面图,第四列为X=-0.25km 处沿Y方向的垂直断面图. Fig. 8 Results from the 3D inversion of 2D data on profile A The top row shows the results from the inversion of only Zxy ,the second row shows the results from the inversion of only Zyx ,the third row shows the results from the inversion of bothZxy and Zyx ,and the bottom row shows the results from the inversion of all impedance tensor elements.The black dashed lines show the prism margins.The first column shows the horizontal slices at 4.25 km depth, the second column shows the vertical slices at Y= -0.25 km along the X axis, the third column shows the vertical slices at Y= -2.875 km along the X axis, and the fourth column ^hows the vertical slices at X= -0.25 km along the Y axis.

图 8的第一行,单元素(Zxy)的反演结果对低阻棱柱体没有反映.图 8 的第二行,单元素(Zyx)的反演结果虽然对低阻棱柱体有反映,但位置偏差比较大:在X方向上,棱柱体的位置相对正确;在Y方向上,棱柱体往Y负方向偏离较大;在深度方面,棱柱体从地表开始就有反映.图 8的第三行,双元素(ZxyZyx)的反演结果比单元素(ZxyZyx)的反演结果合理得多:在XZ方向上,棱柱体的位置都相对正确;在Y方向上,棱柱体往Y负方向偏离.在这四种反演结果中,四个元素(ZxxZxyZyxZyy)同时反演的结果最好:在XZ方向上,棱柱体的位置都相对正确;在Y方向上,虽然棱柱体的位置仍然有些偏差,但偏差在四种结果中最小(见图 8的最后一行).分析四个元素反演的结果优于双元素的反演结果的原因在于:剖面A的对角元素数据含有很强的棱柱体异常信息(见图 2),因此对角元素数据对于还原棱柱体具有重要的作用.

4.1.3 三维反演中选择反演数据的分析

综合分析以上单棱柱体模型合成二维剖面数据的反演结果,无论二维剖面是位于棱柱体中心部位上方还是位于棱柱体边缘部位上方,单元素(ZxyZyx)的反演结果都不能准确地还原出棱柱体,反演结果不是对棱柱体几乎没有反映(见剖面A 和剖面 B的Zxy反演结果),就是棱柱体的位置偏差很大(见剖面A 的Zyx反演结果),甚至得到错误的棱柱体信息,出现两个棱柱体(见剖面B 的Zyx反演结果).因此,对大地电磁二维剖面数据进行三维反演,应该尽可能地避免仅使用单元素数据(ZxyZyx).

对于双元素(ZxyZyx)的反演结果而言,当二维剖面位于棱柱体中心部位上方时,双元素的反演结果能够较合理地还原出棱柱体(见剖面B的Zxy+Zyx反演结果);当二维剖面位于棱柱体边缘部位上方时,双元素的反演结果还原出的棱柱体位置会出现一些偏差(见剖面A 的Zxy+Zyx反演结果).因此,对大地电磁二维剖面数据进行三维反演时,只使用非对角元素(ZxyZyx)的反演结果不是非常可靠.

虽然当二维剖面位于棱柱体中心部位上方时,对角元素数据对于还原棱柱体没有帮助,但是当二维剖面位于棱柱体边缘部位上方时,只有包含对角元素的四个元素(ZxxZxyZyxZyy)的反演结果才能较合理地还原出棱柱体,这也证明了对角元素中包含了许多的三维电阻率结构信息.在实际工作中,二维剖面位于异常体上方的哪个位置是无法事先判定的.因此,在对二维剖面数据进行三维反演时,和非对角元素一样,对角元素在反演中也具有非常重要的作用,应尽量反演所有的张量元素.

4.2 三棱柱体模型合成剖面数据的三维反演

用张量阻抗三维共轭梯度反演程序对剖面A的5个频率(10~0.1 Hz)张量阻抗数据(使用ZxxZxyZyxZyy四个元素)进行三维反演,反演的结果如图 9所示.

图 9 对剖面A 数据进行三维反演的结果 黑色虚线表示棱柱体的边界.第一行为深度4.5km 处的水平截面图,第二行为Y= -1.5km 处沿X方向的垂直断面图,第三行为Y=-0.1km处沿X方向的垂直断面图,第四行为Y=1.5km 处沿X方向的垂直断面图,第五行为X=-12km 处沿Y方向的垂直断面图,第六行为X=0km 处沿Y方向的垂直断面图,第七行为X=12km 处沿Y方向的垂直断面图. Fig. 9 Results from the 3D inversion of 2D data on profile A The black dashed lines show the prism margins.The top row shows the horizontal slices at 4.5 km depth, the second row shows the vertical slices at Y= -1.5 km along the X axis, the third row shows the vertical slices at Y= -0.1 km along the X axis, the fourth row ^hows the vertical slices at Y=1.5 km along the X axis, the fifth r ow ^hows the vertical slices at X= 一12 km along the Y axis, the sixth row show s the vertical slice s at X=0 km along the Y axis, and the bottom row show s the vertical slice s at X=12 km along the Y axi s .

图 9的第一列为真实模型.图 9的第二列是取该剖面二维TM+TE 模式联合反演结果(见图 6)为初始模型的三维反演结果:0~10km 深度范围内200Ωm 的围岩向下和向100Ωm 围岩方向有些延伸,剩余围岩部分无论是边界分布还是电阻率值都基本和真实模型一致;低阻棱柱体ρ1XZ方向上位置相对正确,而在Y方向上边界拉长,电阻率值也接近真实值;低阻棱柱体ρ2XZ方向上位置相对正确,而在Y方向上边界有些拉长,电阻率值还原至20Ωm 左右;低阻棱柱体ρ3Z方向上位置相对正确,而在XY方向上边界缩小,电阻率值还原至50Ωm 左右.图 9的第三列为初始模型取100 Ωm 的三维反演结果:在剖面A 两侧6km范围内的围岩基本与真实模型一致,而6km范围外的围岩电阻率仍为初始模型值100Ωm;三个棱柱体的位置都基本与真实模型一致;低阻棱柱体ρ1 的电阻率值接近于真实值,而低阻棱柱体ρ2ρ3的电阻率值分别还原至20Ωm 和50Ωm 左右.

在三维反演初始模型的选择方面,上文的单棱柱体模型合成剖面数据反演时都选择与围岩一致的100Ωm 为初始模型.在以100Ωm 为初始模型对三棱柱体模型合成剖面数据反演时,三棱柱体模型有三部分围岩,初始模型只与中间的围岩一致.对于另两部分围岩而言,相当于选择了不同的初始模型.这种情况下,虽然三个目标棱柱体还原得比较好,但在200Ωm 和50Ωm 的两部分围岩方面,剖面A 附近的围岩能够较好地还原,而离剖面A 较远处的围岩则不能够恢复.在以二维反演结果为初始模型的三维反演结果中,三部分围岩因为有二维反演结果为基础都还原得比较好,而三个目标棱柱体却形状拉长.分析得到以上结果的原因在于:只使用单一的二维剖面A 数据进行三维反演所获得三维电阻率结构信息仍然是有限的,只能得到剖面附近的三维构造信息,要想更好地还原出围岩和棱柱体(包括形状和电阻率值),需要更多的控制数据.然而,对剖面A的数据进行三维反演至少获得了较合理的剖面A附近的三维电阻率结构分布情况,而且这些信息是二维反演所不能获得的.

5 实测资料的反演

我们用二维共轭梯度反演程序对某地区的实测大地电磁剖面资料进行了二维反演.选取该地区测线1条,测点28个,测线长度6km, 点距200m, 采用8个频点,范围从97~0.35Hz, 初始模型为200Ωm.图 10为该测线数据的二维TM+TE模式联合反演结果.

图 10 实测资料的二维TM+TE 模式联合反演结果 Fig. 10 Results from 2D TM+TE mode inversion of the field data

同样,我们用张量阻抗三维共轭梯度反演程序对该测线的张量阻抗数据进行了三维反演.三维网格剖分为45×75×33(含10个空气层),选取该测线8个频率的ZxxZxyZyxZyy四个元素数据,初始模型为200Ωm, 经过34次迭代反演,耗时42小时9分钟,数据的拟合差从86.9 降到23.5 迭代结束.图 11是该测线数据的三维反演结果.其中,第一行为X=0km 处沿测线方向的断面图,第二行为Y=3km 处垂直测线方向的断面图,第三行为深度1.2 km 处的水平截面图.该地区布置大地电磁工作的目的主要是为了探明地下基岩的起伏情况.从测线下方的断面图看,三维反演结果的基岩界面在沿测线方向2.5~4km 处与二维反演结果存在一定的差异,其余位置与二维反演结果大体一致.从Y=3km处垂直测线方向的断面图和深度为1.2km 处的水平截面图看,测线两侧2km 外地电模型的电阻率基本为初始模型值200Ωm, 也就是说三维反演的结果不能够获得离测线较远处的地质体信息.但是,测线附近的地下基岩的三维起伏情况还是可以通过三维反演结果获得的.

图 11 实测资料的三维反演结果 Fig. 11 Results from 3D inversion of the field data

图 12是实测资料中Zxy数据与三维反演模型阻抗响应的拟断面对比图.可以看出,实测数据Zxy与反演模型的阻抗相应拟断面图吻合得比较好.另外ZxxZyxZyy数据和Zxy数据类似,实测资料和反演模型的阻抗响应都吻合得比较好.因此,可以说得到的地电模型能够代表实测资料所反映的真实地电构造所具有的特征.

图 12 实测资料中Zxy数据与三维反演模型阻抗响应的拟断面对比图 ReZxy表示阻抗Zxy的实部,ImZxy表示阻抗Zxy的虚部. Fig. 12 The comparison between input Zxy impedance data and the impedance response of result model Re Zxy denotes the real part of Zxy and Im Zxy denotes imaginarypart of Zxy
6 结论和建议

在本文中,我们对理论模型的三维阻抗张量响应进行了分析,并且用我们开发的大地电磁张量阻抗三维共轭梯度反演程序对合成数据和实测资料进行了反演,通过分析这些数据的三维反演结果,并将三维反演结果与二维反演的结果相对比,讨论了利用三维反演的方法来获得大地电磁二维剖面附近三维电阻率结构的可行性.合成二维剖面数据和实测资料的三维反演结果表明:用三维反演的方法对大地电磁二维剖面数据进行反演解释是可行的,利用二维剖面的张量数据,不仅可以沿剖面获得较好的二维断面结果,还能够得到剖面附近较可靠的三维电阻率结构信息,这些信息至少是通过二维反演所不能获得的.

在合成数据和实测资料的算例中,我们也发现剖面数据的三维反演结果所能提供的三维构造信息是有限的,在离剖面较远处的三维构造信息是不可靠的.因此,使用二维剖面数据的三维反演结果进行地质解释时需要谨慎使用,在实际工作中,应在重要的关键地区使用.

同时,通过对比理论模型合成二维剖面数据中的不同阻抗数据的三维反演结果,我们发现阻抗张量数据中的对角元素对于圈定剖面附近三维异常体的空间分布具有独特作用.因此,建议在对二维剖面数据进行三维反演时,应尽量反演所有的张量元素.

参考文献
[1] 陈乐寿, 王光锷. 大地电磁测深法. 北京: 地质出版社, 1990 . Chen L S, Wang G E. Magnetotelluric Sounding Method (in Chinese). Beijing: Geological Publishing House, 1990 .
[2] 魏文博. 我国大地电磁测深新进展及瞻望. 地球物理学进展 , 2002, 17(2): 245–254. Wei W B. The progress and looking forward of domestic magnetotelluric sounding. Progress in Geophysics (in Chinese) , 2002, 17(2): 245-254.
[3] 胡祖志, 胡祥云, 何展翔. 三维大地电磁数据的二维反演解释. 石油地球物理勘探 , 2005, 40(3): 353–359. Hu Z Z, Hu X Y, He Z X. Using 2-D inversion for interpretation of 3-D MT data. Oil Geophysical Prospecting (in Chinese) , 2005, 40(3): 353-359.
[4] 陈小斌, 赵国泽, 马霄. 大地电磁三维模型的一维二维反演近似问题研究. 工程地球物理学报 , 2006, 3(1): 9–15. Chen X B, Zhao G Z, Ma X. Research on inversion of MT 3D model approximately by 1D, 2D inversion method. Chinese Journal of Engineering Geophysics (in Chinese) , 2006, 3(1): 9-15.
[5] Ledo J. 2-D versus 3-D magnetotelluric data interpretation. Surveys in Geophysics , 2005, 26: 511-543. DOI:10.1007/s10712-005-1757-8
[6] Ledo J, Queralt P, Marti A, et al. Two-dimensional interpretation of three-dimensional magnetotelluric data: an example of limitation and resolution. Geophys. J. Int. , 2002, 150: 127-139. DOI:10.1046/j.1365-246X.2002.01705.x
[7] Wanamaker P E, Hohmann G W, Ward S H. Magnetotelluric responses of three-dimensional bodies in layered earth. Geophysics , 1984, 49: 1517-1533. DOI:10.1190/1.1441777
[8] Lin C H, Tan H D, Tong T. Three-dimensional conjugate gradient inversion of magnetotelluric sounding data. Applied Geophysics , 2008, 5(4): 314-321. DOI:10.1007/s11770-008-0043-1
[9] Smith J T, Booker J R. Rapid inversion of two- and three- dimensional magnetotelluric data. J. Geophys. Res , 1991, 96: 3905-3922. DOI:10.1029/90JB02416
[10] Mackie R L, Madden T R. Three-dimensional magnetotelluric inversion using conjugate gradients. Geophys. J. Int. , 1993, 115: 215-229. DOI:10.1111/gji.1993.115.issue-1
[11] Newman G A, Alumbaugh D L. Three-dimensional magnetotelluric inversion using non-linear conjugate gradients. Geophys. J. Int. , 2000, 140: 410-424. DOI:10.1046/j.1365-246x.2000.00007.x
[12] Spichak V, Popova. Artificial neural network inversion of magnetotelluric data in terms of three-deimensional earth macroparameters. Geophys. J. Int. , 2000, 142: 15-26. DOI:10.1046/j.1365-246x.2000.00065.x
[13] 谭捍东, 余钦范, JohnBooker, 等. 大地电磁法三维快速松弛反演. 地球物理学报 , 2003, 46(6): 850–855. Tan H D, Yu Q F, Booker J, et al. Three-dimensional rapid relaxation inversion for the magnetotelluric method. Chinese J. Geophys. (in Chinese) , 2003, 46(6): 850-855.
[14] Siripunvaraporn W, Egbert G, Lenbury Y, et al. Three-dimensional magnetotelluric inversion: data-space method. Physics of The Earth and Planetary Interiors , 2005, 150(1-3): 3-14. DOI:10.1016/j.pepi.2004.08.023
[15] 胡祖志, 胡祥云, 何展翔. 大地电磁非线性共轭梯度拟三维反演. 地球物理学报 , 2006, 49(4): 1226–1234. Hu Z Z, Hu X Y, He Z X. Pseudo-three-dimensional magnetotelluric inversion using nonlinear conjugate gradients. Chinese J. Geophys. (in Chinese) , 2006, 49(4): 1226-1234.
[16] 谭捍东, 余钦范, JohnBooker, 等. 大地电磁法三维交错采样有限差分数值模拟. 地球物理学报 , 2003, 46(5): 705–711. Tan H D, Yu Q F, Booker J, et al. Magnetotelluric three-dimensional modeling using the staggered-grid finite difference method. Chinese J. Geophys. (in Chinese) , 2003, 46(5): 705-711.
[17] Lin C H, Tan H D, Tong T. Parallel rapid relaxation inversion of 3D magnetotelluric data. Applied Geophysics , 2009, 6(1): 77-83. DOI:10.1007/s11770-009-0010-5