地球物理学进展  2016, Vol. 31 Issue (4): 1480-1491   PDF    
相位超象限大地电磁观测数据的模型研究——以河西走廊北侧为例
邵贵航, 肖骑彬     
地震动力学国家重点实验室, 中国地震局地质研究所, 北京 100029
摘要: 在大地电磁实测数据中阻抗张量次对角元素的相位值超出正常象限(一、三或二、四象限)的现象称为相位超象限.这种现象无法利用各向同性介质的一、二维反演来解释,而与特定形态的三维各向同性介质有关的电流槽效应有可能是引起相位超象限现象的原因之一.本文首先重建了四个能引起相位超象限的三维各向同性介质理论模型,对模型结构和响应特征进行了总结.随后结合麦克斯韦方程组介绍了电流槽效应产生的原理,并通过近似解析公式分析了电流槽效应对电、磁场的影响.最后,我们针对最近几年来在青藏高原北缘观测到的相位超象限大地电磁测深数据,重点对河西走廊北侧的三个连续相位异常数据进行了三维反演研究,反演模型能够很好地拟合实测的相位超象限数据.通过对模型中存在于上地壳的不同规模高导异常结构进行三维正演检验,最终确定位于测点北东方向的巴丹吉林沙漠构成了区域高导体,它与紧邻测点东侧的局部高导体相连,在平面上形成了“L”型的高导构造,是引起河西走廊北侧相位超象限大地电磁观测数据的主要原因.除此之外,紧邻测点西侧的孤立高导体也会对相位异常的分布范围产生影响.
关键词大地电磁测深     相位超象限     电流槽效应     三维反演     河西走廊    
Model for phase rolling out of quadrant magnetotelluric data—an example from north to the Hexi corridor
SHAO Gui-hang , XIAO Qi-bin     
State Key Laboratory of Earthquake Dynamics, Institute of Geology, China Earthquake Administration, Beijing 100029, China
Abstract: In observed magnetotelluric data, phases of off-diagonal impedance tensors may exceed the normal quadrants (as the first, third or second, forth quadrants), this phenomenon is called Phase Rolling Out of Quadrant (PROQ). The PROQ data cannot be explained with one-dimensional or two-dimensional isotropic methods. The current channeling, which is caused by complex three-dimensional (3-D) isotropic media, is one explanation for the PROQ data. In this study, we reconstructed four synthetic 3-D isotropic models according to former studies and analyzed their structure features and response features. Then we introduced the theory of the current channeling combining with the Maxwell equations, and analyzed the influence of current channeling on electric and magnetic fields with approximate analytical formulas. Finally, 3-D isotropic inversions were applied to three sequential magnetotelluric sites with PROQ at north to the Hexi Corridor, the responses of the final models fit the observed PROQ data very well. With 3-D forward modeling, we tested the contribution of the high-conductivity anomalies in the upper crust to the PROQ data. The regional (large) conductor northeast to the PROQ MT sites is corresponding to the Badain Jaran Desert, it connects a local (small) conductor on the east side of the PROQ MT sites; these two conductors, forming a L-shaped high-conductivity structure in plane view, are the main contributors to the PROQ data north to the Hexi Corridor. Beside that, there is another isolated conductor west to the PROQ sites that can influence the distribution of abnormal phases.
Key words: magnetotelluric sounding     phase rolling out of quadrant     current channeling     3-D inversion     Hexi corridor    
0 引 言

在大地电磁测深的实际观测数据中,阻抗张量Z是一个2×2的复数矩阵.通常情况下,阻抗张量次对角元素(ZxyZyx)的相位值位于一、三或二、四象限内;但是也有极少数测点的ZxyZyx元素的相位会超出正常的象限范围,这种相位超出正常象限的现象称为相位超象限(PROQ,Phase Rolling Out of Quadrant).

Jones等(1988)在文章中首次正式提到大地电磁相位超象限现象,他们在测量科迪勒拉山脉电性结构时,发现位于Nelson基岩(约50 km宽)上的测点阻抗Zyx分量的相位值超出第三象限.此后,在加拿大、智利、纳米比亚、西班牙、印度、日本、芬兰等地的大地电磁观测资料中也见有相位超象限的研究报导(Livelybrooks et al.,1996; Chouteau and Tournerie,2000; Heise and Pous,2003; Lezaeta and Haak,2003; Weckmann et al.,2003; Ichihara and Mogi,2009; Kumar and Manglik,2012; Vaittinen et al.,2012; Kühn et al.,2014; Cherevatova et al.,2015).我们近年在青藏高原北缘获得的大地电磁观测数据中,也发现有部分测点的相位在低频段超出了正常的范围(Xiao et al.,2013; Xiao et al.,2015).由于基于各项同性介质的二维反演技术无法对这些出现在低频段的异常相位进行拟合,我们在作剖面二维反演时不得不剔除这部分异常相位数据.这种处理方式在测线中只有少数几个测点表现为相位超象限时,能够在一定程度上给出区域的深部电性结构,但是所获得的电性结构模型也必然无法合理解释相位超象限现象.

Egbert(1990)Lezaeta和Haak(2003)Weckmann等(2003)Ichihara和Mogi(2009)分别提出了不同的各向同性介质三维理论模型来解释实际观测数据中的相位超象限现象.随着各向同性介质的三维反演技术日趋成熟(Smith and Booker,1991; Mackie and Madden,1993; Newman and Alumbaugh,2000; 谭捍东等,2003; Siripunvaraporn et al.,2005a2005b; 林昌洪等,2011; Egbert and Kelbert,2012; 胡祥云等,2012; Kelbert et al.,2014),关于三维反演技术的应用实例也越来越多(Ingham et al.,2009; Xiao et al.,2011; Tietze and Ritter,2013; Meqbel et al.,2014),因此直接采用三维反演技术对引起相位超象限现象的地下结构进行恢复是可能的,Kühn等(2014)就采用三维反演技术对智利北部的相位超象限实测数据进行了模拟.此外,Pek和Verner(1997)Heise和Pous(2003)、Chen et al.(2009)也通过构建二维各向异性介质模型来解释相位超象限现象.目前国内虽然没有相关研究报导,但是,霍光谱等(2011)秦林江和杨长福(2012)胡祥云等(2013)等人分别开展了基于各向异性介质的大地电磁研究工作.

本文基于三维各向同性介质反演程序Wsinv3dMT(Siripunvaraporn et al.,2005a2005b)和ModEM(Egbert and Kelbert,2012; Kelbert et al.,2014),首先根据前人研究重建各向同性介质三维理论模型,对模型结构及响应特征进行总结,并从理论上分析引起相位超象限现象的原因;然后重点对河西走廊北侧相位超象限大地电磁实测数据进行三维反演解释,分析总结引起相位超象限的地下结构特征.

1 资料采集与数据特征分析

我们所获得的相位超象限大地电磁观测数据位于青藏高原北缘、河西走廊的北侧(图 1中蓝色三角形).野外大地电磁资料采集于2012年6-8月间,观测仪器为德国Metronix公司生产的ADU-07系统.在每个测点采集五个分量的电、磁场信号(即时间序列):三个磁场分量Hx(正南北向)、Hy(正东西向)、Hz(垂直地面)以及两个电场分量Ex(正南北向)、Ey(正东西向).在数据采集过程中设置了4096 Hz、128 Hz和4 Hz三个不同的采样频率,其中低频4 Hz连续采样时间24小时左右.通过数据处理最终得到频率范围在1000~0.0001 Hz的视电阻率和相位曲线,其中测点723、827、838、841、919的yx分量相位以及测点827~829的xy分量相位在低频段超出了正常的象限范围.

图 1 青藏高原北缘河西走廊及邻区地质简图与MT测点分布 图中红色实线表示活动断裂(邓启东等,2003),新生代地层分布据12500000地质图(http://gsd.cgs.cn/download.asp). Q:第四纪;Qh:全新世;Qp:更新世;N:晚第三纪; N2:上新世;N1:中新世;E:早第三纪;E3:渐新世;E2:始新世;E1:古新世;QL:祁连;HC:河西走廊;AL:阿拉善地块;ATF:阿尔金断裂带;NQLF:祁连山北缘断裂;NDXF:大雪山北缘断裂;TLSF:托莱山断裂;KHF:宽滩山-黑山断裂;SLSF:龙首山南缘断裂;YMSF:榆木山断裂;SBSF:北山南缘断裂;EGF:恩格乌苏断裂;NALF:阿拉善北缘断裂.三角形表示MT测点,其中蓝色三角形表示相位超象限测点. Figure 1 Geology and MT site distribution around the Hexi Corridor at north margin of the Tibetan Platean The red lines mean active faults(Deng et al.,2003),the geology data of 1:2.5 M scale from http://gsd.cgs.cn/download.asp. Q: Quaternary; Qh: Holocene; Qp: Pleistocene; N: Neogene; N2: Pliocene; N1: Miocene; E: Paleogene; E3: Oligocene; E2: Eocene; E1: Paleocene; QL: Qilian; HC: the Hexi Corridor; AL: the Alashan Block; ATF: the Altyn Tagh Fault; NQLF: the Northern Qilian Fault; NDXF: the Northern Daxueh Fault; TLSF: the Tuolai Shan Fault; KHF: the Kuantanshan-Heishan Fault; SLSF: the Southern Longshou Shan Fault; YMSF: the Yumu Shan Fault; SBSF: the Southern Beishan Fault; EGF: the Engeerwusu Fault; NALF: the Northern Alashan Fault. The triangles show magnetotelluric sites and the blue triangles show PROQ MT sites.

本文重点针对L8测线中三个连续的相位超象限测点827~829进行解释,图 2所示为这三个测点的实测视电阻率和相位曲线.总体上看,三个测点的视电阻率曲线比较光滑;827号测点的xy分量的阻抗相位在10~1000 s超出了90°,yx分量的阻抗相位在1000 s附近超出了-90°;xy分量的视电阻率曲线在10 s之后的低频段下降显著,与828和829号点存在较大差异.828和829号测点的视电阻率和相位曲线的形态比较接近,两个测点的阻抗相位的xy分量在10 s附近开始超出90°且随周期增大偏离角度越来越大.

图 2 测点L8中827~829三点实测视电阻率和相位曲线 (a)张量阻抗次对角分量ZxyZyx对应的曲线;(b)张量阻抗主对角分量ZxxZyy对应的曲线. Figure 2 Apparent resistivity and phase curves of MT sites 827~829 in profile L8 (a)curves of XY and YX component; (b)curves of XX and YY component.
2 相位超象限的理论模型与理论分析

Egbert(1990)根据地下介质中电流表现出向低阻体汇聚而偏离高阻体的特性,提出了一个概念模型,模型垂向上分为三层:上层为由“S”型高阻体分隔的高导体,中层为高阻体,下层为均匀半空间.这种特殊的电性结构在上层形成迂回的高导通道并引起汇聚电流方向反转,进而导致在低频段电场的相位倒转但磁场相位不变,产生阻抗相位超象限现象.Lezaeta和Haak(2003)为了解释智利安第斯俯冲带大地电磁观测中的相位超象限数据,构建了区域高导异常体和局部脉状高导异常体斜交的三维模型,其中脉状异常体镶嵌于高阻背景之中.该模型中区域异常体模拟智利西侧的太平洋海水,局部异常体模拟安第斯山脉的走滑断层.Weckmann等(2003)在对纳米比亚WF/OL剪切带的剖面大地电磁测深中,发现在位于椭圆形构造体附近的测线上出现了相位超象限现象,该构造体的岩性为高导的含石墨大理石,而围绕该构造体的背景体为高阻的变质岩,因此他们构建了在高阻半空间中嵌入环状高导三维体的模型,并模拟出了相位超象限现象.Ichihara and Mogi(2009)提出了一个在区域高阻背景中存在由大小不等的高导体组成的L型高导异常体的三维模型,Ichihara等(2013)又利用这一模型特征对北海道北部地区的相位超象限大地电磁观测数据分布特征进行了详细分析.

我们依据Egbert(1990)Lezaeta和Haak(2003)Weckmann等(2003)Ichihara和Mogi(2009)提出的模型,重新构建了4个理论模型,如图 3左侧所示.四组模型使用了相同的网格框架,图中显示的模型核心部分平面图中网格数为20×20,每个网格大小为1 km×1 km;在地面以下垂直方向上分为45层,总厚度大于1000 km.我们使用Wsinv3dMT程序中的正演子过程对每个模型的电、磁场响应以及阻抗数据进行了计算.

图 3 三维各向同性介质理论模型及视电阻率和相位响应曲线 模型(a)参考Egbert(1990)的模型描述,模型垂向上分为三层:上层厚度为0.932 km,由1 Ω·m的高导体和1000 Ω·m的S型高阻体组成;中层厚度为0.965 km,由1000 Ω·m的高阻层组成;下层为100 Ω·m的均匀半空间.模型(b)参考Lezaeta和Haak(2003)的模型,该模型分为两层:上层厚度为12.174 km,由0.3 Ω·m的高导体和1000 Ω·m的高阻体组成;下层为1000 Ω·m的均匀半空间.模型(c)参考Weckmann等(2003)的模型,该模型包含两层:上层厚度1.183 km,为电阻率为1 Ω·m的环状高导体嵌于电阻率为1000 Ω·m的背景高阻体中;下层为1000 Ω·m的均匀半空间.模型(d)参考Ichihara和Mogi(2009)提出的模型,分为两层:上层厚度为4.823 km,包含两个大小不等的高导异常体,区域(大的)高导异常体电阻率为0.3 Ω·m,局部(小的)异常体电阻率为3 Ω·m,背景高阻体电阻率为1000 Ω·m;下层为1000 Ω·m的均匀半空间. Figure 3 3-D synthetic models and calculated apparent resistivity and phase curves Model(a)is referred to Egbert(1990)’s model,including three layers: the upper layer,with a thickness of 0.932 km,is composed of is a 1 Ω·m conductive layer divided by 1000 Ω·m S-shaped high-resistivity bodies; the middle layer,with a thickness of 0.965 km,is of 1000 Ω·m; the lower layer is 100 Ω·m homogeneous half-space. Model(b)is referred to Lezaeta and Haak(2003)’s model,including two layers: the upper layer is composed of 0.3 Ω·m conductors and 1000 Ω·m high-resistivity body,with a thickness of 12.174 km; the lower layer is 1000 Ω·m homogeneous half-space. Model(c)is referred to Weckmann et al.(2003)’s model,including two layers: the upper layer includes a 1 Ω·m ring-shaped conductor embedded in 1000 Ω·m high-resistivity body,with a thickness of 1.184 km; the lower layer is 1000 Ω·m homogeneous half-space. Model(d)is referred to Ichihara & Mogi(2009)’s model,including two layers: the upper layer is composed of L-shaped conductors embedded in 1000 Ω·m high-resistivity body,with a thickness of 4.832 km; the L-shaped conductors include a regional conductor(0.3 Ω·m)and a local conductor(3 Ω·m); the lower layer is 1000 Ω·m homogeneous half-space.

图 3a模型在给定点的视电阻率和相位响应曲线来看,e2点的阻抗相位xy分量在频率低于1 Hz时超出了正常象限,出现相位异常的位置位于“反S”型高导异常体的正上方;从阻抗相位在1000 s的切片来看,异常相位只出现在xy分量上(图 4a),在阻抗Zxy相位超出正常象限的位置,电场Ex的相位出现异常,而磁场Hy的相位未出现异常,这说明电场异常是引起相位超象限的主要原因.

图 4 三维各向同性介质理论模型在周期为1000 s时的阻抗相位切片及对应的电、磁场复数矢量,其中箭头长短对应电、磁场的幅值,箭头方向对应电、磁场相位.(a)~(d)分别对应图 3中理论模型的响应 Figure 4 Phase slices with the electic and magnetic field complex vectors at period 1000s of the synthetic 3-D models, the arrow length and direction are corresponding to amplitude and phases of the field. (a)~(d)are the responses of each models in Fig. 3

图 3b模型在给定点的视电阻率和相位曲线来看,l3点的阻抗相位yx分量在频率低于0.017 Hz时出现了超象限,阻抗相位出现异常的位置处于区域异常体和局部异常体交汇的内侧位置,且位于紧贴局部高导异常体的高阻体之上(图 4b).同样的,在阻抗Zyx相位异常的位置,电场Ey相位同样出现异常,而磁场Hx相位未受影响,而且在高导异常体上方的电场幅值小于区域背景值,磁场幅值大于背景值.

图 3c模型在给定点的视电阻率和相位曲线来看,w2点的阻抗相位yx分量在频率低于10 Hz时出现了超象限,w4点的阻抗相位xy分量同样在频率低于10 Hz时出现了超象限.阻抗相位xy分量和yx分量出现异常的位置都处于环状体内侧,且位于紧贴环状高导异常体的高阻体之上(图 4c).在阻抗Zxy相位异常的位置,电场Ex相位出现异常;在阻抗Zyx相位异常的位置,电场Ey相位同样出现异常,而对应的磁场HxHy相位未受影响,同样在高导异常体上方的电场幅值小于区域背景值,磁场幅值大于背景值.

图 3d模型在给定点的视电阻率和相位曲线来看,i2点的阻抗相位yx分量在频率低于1 Hz时出现了超象限.阻抗相位出现异常的位置处于局部高导异常体与区域高导异常体交汇的内侧位置(图 4d).在阻抗Zyx相位异常的位置,电场Ey相位同样出现异常,而磁场Hx相位未受影响.

通过对前人四个理论模型的重构和三维正演计算,我们发现除了模型图 3a,其余模型及其响应具有以下共同特征:

(1) 引起相位异常的模型具有相交的高导异常结构,并且高导异常体镶嵌于高阻体之中;

(2) 阻抗相位出现异常的位置都在高导异常体旁侧的高阻体上方;如果存在大小不等的高导体,异常相位出现在局部(较小的)高导体附近;

(3) 阻抗相位异常出现在低频段;

(4) 阻抗相位异常主要是电场相位变化引起的.

模型图 3a所产生的异常相位出现在高导体的上方,异常阻抗相位对应的电场相位发生了180°反转.Egbert(1990)认为浅层反“S”型高导异常引起的电流槽效应(Current channelling)会使浅部的电流发生集中流动,受反“S”型的高导体形态影响,流动方向反转的电流对应的低频时的电场相位也发生约180°的变化但磁场相位未受影响(如图 4a),从而导致在高导异常体上方出现相位超象限现象,这是对该模型很好的定性解释.然而Egbert提出的电流汇聚引起的电流槽效应无法解释其他三个模型的响应特征,因为其他三个模型发生相位超象限的点并不位于高导体上方,这需要从其他方面来分析.

Price(1973)指出当区域电流流经电导率不同的介质分界面时,由于界面两侧介质的电导率不同,有一部分传导电流(如图 5中的jc)被阻隔在间断面处形成表面电荷,如图 5.这种间断面上表面电荷的累积可以看成一种电容器模型,由表面电荷引起的位移电流对磁场的影响可以忽略,但表面电荷产生的静电场则不能忽略.Jones(1983)将这种由表面电荷引起的不可忽略的静电场也称之为电流槽效应.

图 5 地下介质电性间断面处表面电荷累积示意图 (根据Price(1973)Jones(1983)重画),其中电导率σ12 Figure 5 Surface charge distributions near a vertical conductivity discontinuity in conductivity with σ12 (redraw from Price(1973)and Jones(1983))

为了研究电流槽效应的影响范围,Le Mouel和Menvielle(1982)分析了在均匀半空间中存在长条状三维高导异常体的模型,其中三维异常体的特征长度为l.引入代表磁场的矢量势A和代表电场的标量势Φ两个物理量,并将实测电场E分为正常场En和异常场Ea,正常电场由均匀半空间所产生,不存在表面电荷引起的标量势Φ及静电场.而在异常体附近,由异常体引起的电、磁场(Ea、Ha)表达式为

(1)
(2)
(3)
(4)

其中Ea表示异常电场,Aa表示异常体中的磁场矢量势,Ha表示异常磁场,μ为磁导率,Ja表示异常体中电流密度,σ为异常体电导率.

将公式(1)(2)带入式(3)(4)再合并,可得:

(5)

为了简化计算,对磁场矢量势的两次旋度计算替换为对特征长度的两次相除,将对时间的偏导替换为对特征时间T(近似于周期)的相除.由此,式5可简化为近似等式:

(6)

其中l表示地下异常体的特征长度,T表示特征时间,‘~’号表示符号两侧的式子不是严格相等而是在同一数量级上.

式6进一步整理可得:

(7)

其中δ表示趋肤深度,当异常体的宽度相对于长度不可忽略时,可用特征面积S来代替l2.

由公式7可以总结出,当电、磁场的趋肤深度δ远大于异常体的特征长度l(或δ2远大于S)时,磁场感应部分()要远小于静电场部分(),此时表面电荷引起的电流槽效应是引起异常电场的主要因素.

这种表面电荷引起的电流槽效应相较于电流在高导体中汇聚流动引起的电流槽效应更能解释图 3bd中的理论模型响应特征,因为在电性间断面两侧,高阻体和高导体中的电场都受到电流槽效应影响,当频率相同时,电场在高阻体中的趋肤深度大于高导体中的,因此更容易满足“趋肤深度远大于异常体特征长度”这个条件,即电流槽效应对高阻体的电场影响更大,所以相位异常出现在高阻体上方.根据以上分析,我们认为引起相位超象限的电流槽效应分为两种情况:一种是区域电流受高导体形态影响而发生汇聚流动和方向变化,这时引起的相位超象限发生在高导体正上方;另一种是区域电流在高低阻介质的界面处产生电荷累积并形成静电场,这时引起的相位超象限发生在高阻体上方.

3 河西走廊北侧相位超象限数据反演与模型特征 3.1 三维反演与模型特征

以相位超象限测点827~829为中心,我们根据L8测线上825~831七个连续测点的坐标位置构建了三维初始模型,模型的网格数为X×Y×Z=25×26×82,其中X代表南北方向网格数,Y代表东西方向网格数,Z代表垂向地心方向网格数(包含7层空气层),模型的初始电阻率为100 Ω·m.图 6a所示为模型的核心部分水平网格和垂直分层.反演中采用了全部阻抗张量元素(Zxx、Zxy、Zyx、Zyy),数据的频率范围从809.45~0.000046 Hz.

图 6 (a)三维初始模型核心部分的平面网格剖分图及模型的垂向分层图;(b1)Wsinv3dMTV程序多阶段反演的RMS值变化图;(b2)ModEM程序反演的RMS变化图;(c1)Wsinv3dMT程序反演结果模型切片图;(c2)ModEM程序反演结果模型切片图 Figure 6 (a) Plane view of core section and vertical layers from surface down of the 3-D initial model;(b1)Plots of RMS misfit for each iteration in Wsinv3dMT inversion;(b2)Plots of RMS misfit for each iteration in ModEM inversion;(c1)the 3-D resisticity slices in W-E direction of the Wsinv3dMT inversion model; (c2)the 3-D resisticity slices of the ModEM inversion model

反演中分别采用了Wsinv3dMT(Siripunvaraporn et al.,2005a,2005b)和ModEM(Egbert and Kelbert,2012; Kelbert et al.,2014)程序.其中在进行Wsinv3dMT反演时,空间光滑参数(τ,δx,δy,δz)分别为5,0.1,0.1,0.1,正则化参数选为2,改变步长为0.5.为了提高数据拟合效果,在使用Wsinv3dMT程序时采用了多阶段的反演过程,即选取当前反演过程中RMS最小的模型作为初始模型,在其他参数保持不变的情况下重新进行反演,反演过程中RMS值变化如图 6b1所示.这一反演过程已经在实际应用中被证实可以有效降低模型粗糙度并提高数据拟合效果(Xiao et al.,2012).我们采取了两个阶段,并选取在第二个阶段反演中具有最小RMS值的模型作为最终模型(图 6c1).在利用ModEM程序进行三维反演时,我们同样选取正则化参数2,反演过程中RMS值会随迭代次数增加而降低,我们选取了最后一次迭代的结果作为最终模型(图 6c2).由图 6c1、6c2可以看出,Wsinv3dMT程序和ModEM程序反演结果模型中都在相同位置出现了高导异常体,且高、低阻体之间的位置关系基本一致,只是后者的模型中高导异常体比前者的尺度要小一些,总体上看两个程序的反演结果可以相互对比.

图 7为采用这两种反演程序得到的最终模型的视电阻率和阻抗相位拟合曲线图.可以看出,两种反演方法都能拟合相位超象限数据,ModEM程序对观测数据的拟合程度总体要高.

图 7 三维反演模型响应曲线与原始观测曲线对比 (a)Wsinv3dMT程序对应结果;(b)ModEM程序对应结果. Figure 7 Comparisons of 3-D model responses with the observed curves (a)Results of Wsinv3dMT program;(b)Results of Mod3dEM program.

图 8为两种反演程序的结果模型核心部分在不同深度处的电阻率切片图.可以看出,Wsinv3dMT程序和ModEM程序的反演结果模型的核心部分形态上十分相近,都是在浅层存在大小不一的高导异常体,随深度增加高阻体开始大面积的出现.两个模型都反映出了高导异常体镶嵌于高阻背景之中的结构特征;在高导异常体规模上,后者比前者要小.

图 8 (a)Wsinv3dMT程序反演结果模型不同深度处电阻率切片,(b)ModEM程序反演结果模型不同深度处电阻率切片 Figure 8 (a)The resistivity slices of 3-D inversion model calculated by Wsinv3dMT program at different depth; (b)The resistivity slices of 3-D inversion model calculated by Mod3dEM program at different depth
3.2 引起相位超象限的异常结构分析

图 1来看,测点825~831所在位置全部被第四纪沉积物所覆盖,这些松散沉积物有可能在浅部形成规模不等的高导异常.在测点东北方为巴丹吉林沙漠,有可能形成区域上的低阻构造;在图 8所示在5 km深度以内的电阻率切片以及图 6c1图 6c2中,测点的东北方位存在较大范围的低阻体(即图 8中出现的B1).为了研究引起相位超象限的模型特征,我们以Wsinv3dMT的反演结果模型为基础,对图 8a中的(相对)高导块体A1、A2、B1、B2、C、D1、D2进行验证.将高导异常块体依次用100 Ω·m的均匀介质代替,再通过Wsinv3dMT中的三维正演子过程对修改后的模型进行三维正演计算,观察模型对应的阻抗相位变化.

图 9上半部分为根据上述方法计算得到的阻抗相位切片图.由图中可以看到,几个相对高导(电阻率在几欧姆米到十几欧姆米)异常块体A1、C、D1、D2被代替后计算得到的阻抗相位切片几乎未受到影响,所以可以认为异常块体A1、C、D1、D2对相位异常的产生是没有贡献的.然而,将异常块体A2替换掉后,Zxy分量相位异常的分布范围发生了一定变化,但测点827~829三点仍为相位异常点;同时测点826~828的Zyx分量阻抗也受到了影响,这说明A2块体对相位异常的产生是有一定贡献的.当替换掉块体B1后,测点827的Zxy分量相位不再出现异常;当替换掉块体B2后,827、828的Zxy分量相位不再出现异常,而且Zxy分量的相位异常分布范围明显减小.这说明异常块体A2、B1、B2在引起Zxy分量相位异常上都起到了一定作用,而B1和B2的影响更为明显.

图 9 不同三维模型在周期为854 s的ZxyZyx相位切片.其中orig代表未经改变的反演模型对应的相位切片,A1、A2、B1、B2、C、D1、D2、A1A2、B1B2、A2B2、A1A2B1B2、CD1D2、A1CD1D2、A1A2CD1D2、A1B1CD1D2分别代表将对应异常块体或组合用100 Ω·m的均匀介质代替后计算得到的相位切片 Figure 9 Phase slices of Zxy and Zyxat period 854s calculated from different models. Orig means the phase slices of the unchanged model; A1,A2,B1,B2,C,D1,D2,A1A2,B1B2,A2B2,A1A2B1B2,CD1D2,A1CD1D2,A1A2CD1D2,A1B1CD1D2 mean the phase slices of the modified models with the corresponding abnormal bodies replaced by 100 Ω·m homogeneous media

下面再对A2、B1、B2这三个块体进行更加细致的研究,对多个组合块体进行替换,这些组合包括:A1和A2(简称A1A2,下同)、B1B2、A2B2、A1A2B1B2、CD1D2、A1CD1D2、A1A2CD1D2、A1B1CD1D2.计算结果如图 9下半部分所示,当替换掉A1CD1D2时,计算得到的阻抗相位基本不受影响,说明这四个块体不管是单独存在还是共同存在对相位异常的产生都是没有贡献的.当替换的块体中包含有A2、B1、B2时,相位或多或少都会受到影响,说明这三个块体都对相位异常的产生有贡献,缺少任何一个都会改变相位异常的范围.尤其是当B2与A2、B1之间两两组合时,Zxy分量的相位超象限现象消失,进一步说明了B2的大小和形态影响了相位异常的范围.

3.3 抽象模型

以上正演验证过程显示高导异常块体A2、B1、B2是引起河西走廊北侧地区相位超象限的主要贡献者,其中B2对测点827~829中阻抗相位的影响最为明显.依据图 8中的电阻率切片,可以看到这三个块体的形态位置特征:B1是三个块体中最大的,可以对应于Ichihara and Mogi(2009)的理论模型中的区域异常体,同时它可以与测点东北侧的巴丹吉林沙漠对应;B2对应于局部异常体,但从深度分布上B2比B1位置更浅,B2底部与B1相连(见图 8中2066~2363 m切片),这是与Ichihara模型不同的地方;A2的存在会对相位异常的范围有一定作用.我们据此抽象出了两个相对简单的理论模型:第一个模型如图 10a所示,主要高导异常体只包括B1和B2;第二个模型在此基础上增加了A2块体(图 10b).对这两个抽象模型我们分别计算了它们的三维响应.

图 10 依据河西走廊北侧的三维反演结果分析抽象出的理论模型及其响应 模型(a)包含两个高导异常体,分别是代表B1的区域块体和代表B2的局部块体;模型分为三层:上层厚度为1.183 km,为100 Ω·m的均匀介质;中层厚度为23.146 km,在深度0.932~2.398 km处为电阻率1 Ω·m、范围 10 km×10 km的局部异常体,在深度1.897~6.082 km处为电阻率1 Ω·m、范围30 km×95 km的区域异常体,两者之间存在连接通道,背景电阻率为1000 Ω·m;下层为100 Ω·m的均匀半空间.模型(b)与模型(a)不同之处是包含三个高导异常体,在深度1.897~6.082 km增加了一个电阻率3 Ω·m、范围10 km×20 km的局部异常体.其中指定点a1、a2、a3及b1、b2、b3分别对应于反演模型中827~829三个相位异常测点.(c)和(d)为两个抽象模型在周期为1000 s时的阻抗相位切片及对应的电、磁场复数矢量. Figure 10 Plane view and profile view of the core parts of two synthetic models together with their 3-D responses Model(a)contains two conductors which represent the regional conductor B1 and the local conductor B2,it includes three layers: the upper layer is of 100 Ω·m with a thickness of 1.183 km; the middle layer is 23.146 km in thick,including a 10 km×10 km local conductor(of 1 Ω·m)at depth from 0.0932 km to 2.398 km and a 30 km×95 km regional conductor(of 1 Ω·m)at depth from 1.897 km to 6.082 km,the two conductors are connecting with each other; the lower layer is 100 Ω·m homogeneous half-space. Different from model(a),Model(b)has three conductors,a 10 km×20 km local conductor(of 3 Ω·m)at depth from 1.897 km to 6.082 km is added in this model. The sites a1,a2,a3 and b1,b2,b3 are corresponding to sites 827~829 in the inversion models.(c)and(d)show off-diagonal impedance phase slices together with electric and magnetic field complex vectors of the synthetic models.

图 10左侧所示为抽象模型的核心部分网格图,网格数为20×20,每个网格大小5 km×5 km;垂向方向分45层,总深度大于1000 km.使用Wsinv3dMT程序中正演子过程对每个模型的电、磁场响应以及阻抗数据进行了计算.在理论模型中只考虑B1和B2时(图 10a),所对应的相位响应曲线在a1点并没有超出正常象限范围;而在理论模型中同时考虑B1、B2和A2时(图 10b),与a1点位置对应的b1点的阻抗相位曲线则超出了正常象限范围,这说明异常块体A2确实对阻抗相位的异常范围起到了一定的控制作用.图 10c10d为抽象出的模型在周期为1000 s时的阻抗相位切片及相应的电、磁场复数矢量.通过对比xy分量相位切片可以清楚的看到块体A2对相位异常范围的影响.

4 总结与讨论

4.1 为了能够对青藏高原北缘、河西走廊北侧实测的相位超象限数据进行解释,我们首先对前人提出的能引起相位超象限的理论模型进行了重新构建,并计算了它们的阻抗相位和电、磁场响应.通过分析异常相位产生的位置及其对高导异常体的关系,我们总结了在三维各向同性介质下引起大地电磁相位超象限的结构特征:即模型中具有相交的高导异常结构或高导体呈迂回状,并且高导异常体镶嵌于高阻体之中.这种特殊结构能够产生电流槽效应,电流槽效应能够在两种情况下引起相位超象限:

一种是区域电流受高导体形态影响而发生汇聚流动和方向变化,这时引起的相位超象限发生在高导体的正上方;

另一种是区域电流在高低阻介质的界面处产生累积电荷并形成静电场,这时引起的相位超象限发生在紧邻高导异常体的高阻体上方.

4.2 本文重点就第二种情况从理论公式角度对累积电荷影响电、磁场的过程进行了探究.当然,文中并没有就宏观电流流动方向变化对电、磁场的影响进行更加深入的研究分析.在对累积电荷引起的电流槽效应理论分析中,我们没有得出可以将电流槽效应与相位超象限直接联系起来的近似解析公式,这些需要在今后做进一步分析研究.

4.3 针对河西走廊北侧实际的相位超象限数据,我们同时采用了Wsinv3dMT和ModEM这两种三维反演程序进行计算,所得到的模型不仅能够拟合相位超象限数据,而且具有很强的对比性.通过对模型中高导异常体进行逐一验证,我们发现引起河西走廊北侧相位超象限的高导异常体的结构与Ichihara and Mogi(2009)提出的理论模型有相似之处,即在平面上都存在“L”型的高导构造.但细节上又存在一定差异:我们的反演模型中区域异常体(图 8中B1,可以与巴丹吉林沙漠对应)和局部异常体(图 8中B2,与测点附近的新生代沉积层对应)不在同一深度范围,并且在“L”型异常体之外、测点的西侧还存在一个相对高导的异常体(图 8中A2),它同样对相位异常的范围有一定影响.

致谢 本文受国家自然科学基金(41274080)资助.感谢赵凌强、秦作权、裴海雄等在野外工作上的帮助.感谢审稿专家和编辑对本文提出的修改意见.感谢Weerachai Siripunvaraporn等提供的Wsinv3dMT以及Gary Egbert等提供的ModEM大地电磁三维反演程序.
参考文献
[1] Cherevatova M, Smirnov M Y, Jones A G, et al.2015. Magnetotelluric array data analysis from north-west Fennoscandia[J]. Tectonophysics, 653 : 1–19. DOI:10.1016/j.tecto.2014.12.023
[2] Chouteau M, Tournerie B. 2000. Analysis of magnetotelluric data showing phase rolling out of quadrant (PROQ)[C].//2000 SEG Annual Meeting. Expanded Abstracts, 334-336.
[3] Deng Q D, Zhang P Z, Ran Y K, et al.2003. Active tectonics and earthquake activities in China[J]. Earth Science Frontiers (in Chinese), 10 (S1) : 66–73.
[4] Egbert G D.1990. Comments on'concerning dispersion relations for the magnetotelluric impedance tensor' by E[J]. Yee and K. V. Paulson[J]. Geophysical Journal International, 102 (1) : 1–8.
[5] Egbert G D, Kelbert A.2012. Computational recipes for electromagnetic inverse problems[J]. Geophysical Journal International, 189 (1) : 251–267. DOI:10.1111/j.1365-246X.2011.05347.x
[6] Heise W, Pous J.2003. Anomalous phases exceeding 90°in magnetotellurics: Anisotropic model studies and a field example[J]. Geophysical Journal International, 155 (1) : 308–318. DOI:10.1046/j.1365-246X.2003.02050.x
[7] Hu X Y, Li Y, Yang W C, et al.2012. Three-dimensional magnetotelluic parallel inversion algorithm using data space method[J]. Chinese Journal of Geophysics (in Chinese), 55 (12) : 3969–3978. DOI:10.6038/j.issn.0001-5733.2012.12.009
[8] Hu X Y, Huo G P, Gao R, et al.2013. The magnetotelluric anisotropic two-dimensional simulation and case analysis[J]. Chinese Journal of Geophysics (in Chinese), 56 (12) : 4268–4277. DOI:10.6038/cjg20131229
[9] Huo G P, Hu X Y, Liu M.2011. Review of the forward modeling of magnetotelluric in the anisotropy medium research[J]. Progress in Geophysics (in Chinese), 26 (6) : 1976–1982. DOI:10.3969/j.issn.1004-2903.2011.06.011
[10] Ichihara H, Mogi T.2009. A realistic 3-D resistivity model explaining anomalous large magnetotelluric phases: The L-shaped conductor model[J]. Geophysical Journal International, 179 (1) : 14–17. DOI:10.1111/j.1365-246X.2009.04310.x
[11] Ichihara H, Mogi T, Yamaya Y.2013. Three-dimensional resistivity modelling of a seismogenic area in an oblique subduction zone in the western Kurile arc: Constraints from anomalous magnetotelluric phases[J]. Tectonophysics, 603 : 114–122. DOI:10.1016/j.tecto.2013.05.020
[12] Ingham M R, Bibby H M, Heise W, et al.2009. A magnetotelluric study of Mount Ruapehu volcano, New Zealand[J]. Geophysical Journal International, 179 (2) : 887–904. DOI:10.1111/gji.2009.179.issue-2
[13] Jones A G.1983. The problem of current channelling: A critical review[J]. Geophysical surveys, 6 (1-2) : 79–122. DOI:10.1007/BF01453996
[14] Jones A G, Kurtz R D, Oldenburg D W, et al.1988. Magnetotelluric observations along the lithoprobe southeastern canadian cordilleran transect[J]. Geophysical Research Letters, 15 (7) : 677–680. DOI:10.1029/GL015i007p00677
[15] Kelbert A, Meqbel N, Egbert G D, et al.2014. ModEM: A modular system for inversion of electromagnetic geophysical data[J]. Computers & Geosciences, 66 : 40–53. DOI:10.1016/j.cageo.2014.01.010
[16] Kumar G P, Manglik A.2012. Electrical anisotropy in the main central Thrust Zone of the Sikkim Himalaya: Inference from anomalous MT phase[J]. Journal of Asian Earth Sciences, 57 : 120–127. DOI:10.1016/j.jseaes.2012.06.017
[17] Kühn C, Küster J, Brasse H.2014. Three-dimensional inversion of magnetotelluric data from the Central Andean continental margin[J]. Earth, Planets and Space, 66 (112) : 1–13.
[18] Le Mouel J L, MenVielle M.1982. Geomagnetic variation anomalies and deflection of telluric currents[J]. Geophysical Journal of the Royal Astronomical Society, 68 (3) : 575–587. DOI:10.1111/j.1365-246X.1982.tb04916.x
[19] Lezaeta P, Haak V.2003. Beyond magnetotelluric decomposition: Induction, current channeling, and magnetotelluric phases over 90°[J]. Journal of Geophysical Research: Solid Earth, 108 (B6) : 2305. DOI:10.1029/2001jb000990
[20] Lin C H, Tan H D, Tong T.2011. The possibility of obtaining nearby 3D resistivity structure from magnetotelluric 2D profile data using 3D inversion[J]. Chinese Journal of Geophysics (in Chinese), 54 (1) : 245–256. DOI:10.3969/j.issn.0001-5733.2011.01.026
[21] Livelybrooks D, Mareschal M, Blais E, et al.1996. Magnetotelluric delineation of the Trillabelle massive sulfide body in Sudbury, Ontario[J]. Geophysics, 61 (4) : 971–986. DOI:10.1190/1.1444046
[22] Mackie R L, Madden T R.1993. Three-dimensional magnetotelluric inversion using conjugate gradients[J]. Geophysical Journal International, 115 (1) : 215–229. DOI:10.1111/gji.1993.115.issue-1
[23] Meqbel N M, Egbert G D, Wannamaker P E, et al.2014. Deep electrical resistivity structure of the northwestern U[J]. S. derived from 3-D inversion of USArray magnetotelluric data[J]. Earth and Planetary Science Letters, 402 : 290–304. DOI:10.1016/j.epsl.2013.12.026
[24] Newman G A, Alumbaugh D L.2000. Three dimensional magnetotelluric inversion using non-linear conjugate gradients[J]. Geophysical Journal International, 140 (2) : 410–424. DOI:10.1046/j.1365-246x.2000.00007.x
[25] Pek J, Verner T.1997. Finite-difference modelling of magnetotelluric fields in two-dimensional anisotropic media[J]. Geophysical Journal International, 128 (3) : 505–521. DOI:10.1111/gji.1997.128.issue-3
[26] Price A T.1973. The theory of geomagnetic induction[J]. Physics of the Earth and Planetary Interiors, 7 (3) : 227–233. DOI:10.1016/0031-9201(73)90049-6
[27] Qin L J, Yang C F.2012. A magnetotelluric inversion method of the whole tensor impedance response to one-dimensional anisotropic structure[J]. Chinese Journal of Geophysics (in Chinese), 55 (2) : 693–701. DOI:10.6038/j.issn.0001-5733.2012.02.033
[28] Siripunvaraporn W, Egbert G, Lenbury Y, et al.2005a. Three-dimensional magnetotelluric inversion: Data-space method[J]. Physics of the Earth and Planetary Interiors, 150 (1-3) : 3–14. DOI:10.1016/j.pepi.2004.08.023
[29] Siripunvaraporn W, Egbert G, Uyeshima M.2005b. Interpretation of two-dimensional magnetotelluric profile data with three-dimensional inversion: Synthetic examples[J]. Geophysical Journal International, 160 (3) : 804–814. DOI:10.1111/j.1365-246X.2005.02527.x
[30] Smith J T, Booker J R.1991. Rapid inversion of two- and three-dimensional magnetotelluric data[J]. Journal of Geophysical Research, 96 (B3) : 3905–3922. DOI:10.1029/90JB02416
[31] Tan H D, Yu Q F, Booker J, et al.2003. Three-dimensional rapid relaxation inversion for the magnetotelluric method[J]. Chinese Journal of Geophysics (in Chinese), 46 (6) : 850–854. DOI:10.3321/j.issn:0001-5733.2003.06.019
[32] Tietze K, Ritter O.2013. Three-dimensional magnetotelluric inversion in practice-the electrical conductivity structure of the San Andreas Fault in Central California[J]. Geophysical Journal International, 195 (1) : 130–147. DOI:10.1093/gji/ggt234
[33] Vaittinen K, Korja T, Kaikkonen P, et al.2012. High-resolution magnetotelluric studies of the Archaean-Proterozoic border zone in the Fennoscandian Shield, Finland[J]. Geophysical Journal International, 188 (3) : 908–924. DOI:10.1111/j.1365-246X.2011.05300.x
[34] Weckmann U, Ritter O, Haak V.2003. A magnetotelluric study of the Damara Belt in Namibia: 2[J]. MT phases over 90° reveal the internal structure of the Waterberg Fault/Omaruru Lineament[J]. Physics of the Earth and Planetary Interiors, 138 (2) : 91–112. DOI:10.1016/s0031-9201(03)00079-7
[35] Xiao Q B, Shao G H, Liu-Zeng J, et al.2015. Eastern termination of the Altyn Tagh Fault, western China: Constraints from a magnetotelluric survey[J]. Journal of Geophysical Research: Solid Earth, 120 (5) : 2838–2858. DOI:10.1002/2014jb011363
[36] Xiao Q B, Zhao G Z, Dong Z Y.2011. Electrical resistivity structure at the northern margin of the Tibetan Plateau and tectonic implications[J]. Journal of Geophysical Research, 116 (B12) : B12401. DOI:10.1029/2010jb008163
[37] Xiao Q B, Zhang J, Wang J J, et al.2012. Electrical resistivity structures between the Northern Qilian Mountains and Beishan Block, NW China, and tectonic implications[J]. Physics of the Earth and Planetary Interiors : 200-201–92-104. DOI:10.1016/j.pepi.2012.04.008
[38] Xiao Q B, Zhang J, Zhao G Z, et al.2013. Electrical resistivity structures northeast of the Eastern Kunlun Fault in the Northeastern Tibet: Tectonic implications[J]. Tectonophysics, 601 : 125–138. DOI:10.1016/j.tecto.2013.05.003
[39] 邓起东, 张培震, 冉勇康, 等.2003. 中国活动构造与地震活动[J]. 地学前缘, 10 (S1) : 66–73.
[40] 胡祥云, 李焱, 杨文采, 等.2012. 大地电磁三维数据空间反演并行算法研究[J]. 地球物理学报, 55 (12) : 3969–3978. DOI:10.6038/j.issn.0001-5733.2012.12.009
[41] 胡祥云, 霍光谱, 高锐, 等.2013. 大地电磁各向异性二维模拟及实例分析[J]. 地球物理学报, 56 (12) : 4268–4277. DOI:10.6038/cjg20131229
[42] 霍光谱, 胡祥云, 刘敏.2011. 各向异性介质中大地电磁正演研究综述[J]. 地球物理学进展, 26 (6) : 1976–1982. DOI:10.3969/j.issn.1004-2903.2011.06.011
[43] 林昌洪, 谭捍东, 佟拓.2011. 利用大地电磁三维反演方法获得二维剖面附近三维电阻率结构的可行性[J]. 地球物理学报, 54 (1) : 245–256. DOI:10.3969/j.issn.0001-5733.2011.01.026
[44] 秦林江, 杨长福.2012. 大地电磁全张量响应的一维各向异性反演[J]. 地球物理学报, 55 (2) : 693–701. DOI:10.6038/j.issn.0001-5733.2012.02.033
[45] 谭捍东, 余钦范, BookerJ, 等.2003. 大地电磁法三维快速松弛反演[J]. 地球物理学报, 46 (6) : 850–854. DOI:10.3321/j.issn:0001-5733.2003.06.019