地球物理学报  2021, Vol. 64 Issue (6): 2108-2126   PDF    
大地电磁二维各向异性反演及其在青藏高原北部的应用
喻国1,2, 肖骑彬1,2, 李满1,2     
1. 中国地震局地质研究所, 北京 100029;
2. 中国地震局地质研究所地震动力学国家重点实验室, 北京 100029
摘要:地球介质的电性特征在不同深度上都可能表现出各向异性,识别电各向异性有助于深入理解地球内部物质状态及变形环境,促进对动力学模型的理解.本研究在已实现的任意电各向异性大地电磁有限差分正演算法的基础之上,通过构建常规的各向异性反演目标函数,完成对目标函数及其梯度的计算.借助成熟的非线性共轭梯度反演技术,实现了对大地电磁各向异性介质的反演过程.在此基础上,构建了二维方位各向异性理论模型,来验证反演过程的稳定性.理论模型的反演结果直观地显示了大地电磁各向异性反演对模型参数中垂直电导率恢复上的局限性;在各向异性结构分析上需要将电导率张量作为整体考虑.对青藏高原北部东昆仑—柴达木盆地西段的实测大地电磁数据进行反演,并对比早期二维各向同性反演结果发现,各向同性结构中位于祁漫塔格山脉下方上地幔顶部存在低阻异常带,在相同的位置,各向异性结构中表现出明显的方位各向异性,其各向异性低阻主轴电阻率指示的剪切带走向与地表造山带走向不一致,这暗示该处的应力环境复杂,除了受印度—欧亚板块碰撞控制外,还可能受到邻近阿尔金走滑断裂带左行剪切运动的影响.
关键词: 大地电磁      各向异性反演      东昆仑断裂      柴达木盆地      阿尔金断裂     
Two-dimensional magnetotelluric inversion with anisotropy and its application in the northern Tibetan Plateau
YU Guo1,2, XIAO QiBin1,2, LI Man1,2     
1. Institute of Geology, China Earthquake Administration, Beijing 100029, China;
2. State Key Laboratory of Earthquake Dynamics, Institute of Geology, China Earthquake Administration, Beijing 100029, China
Abstract: Electrical anisotropy exists at all scales in earth. Identifying electrical anisotropy contributes to understanding of the physical state in depth, deformable environment and dynamic models within the Earth. In this study, based on the realized finite-difference magnetotelluric forward algorithm considering arbitrary electrical anisotropy, we construct the objective function for anisotropic inversion and calculate its value and gradient. Then, we complete the anisotropic magnetotelluric inversion algorithm with mature non-linear conjugate gradient inversion technology. A two-dimensional synthetic model is used to examine the stability of the inversion algorithm. The synthetic examples directly show the difficulty in recovering vertical resistivity. And, the elements of conductivity tensor should be considered simultaneously as a whole in the anisotropic interpretation. Along a magnetotelluric profile across the east Kunlun to the west Qaidam Basin in northern Tibetan Plateau, the previous two-dimensional isotropic inversion results have imaged a low resistive anomaly in the uppermost mantle beneath Qiman Tagh range. An apparent azimuthal anisotropy structure appears at the same place in the anisotropic inversion results. The shear direction indicated by the minimum principal resistivity is not parallel to the surface orogenic belt, which implies that, in addition to the control of the collision of India-Eurasia plate, the stress environment is possibly affected by left-lateral shear in the strike-slip Altyn Tagh Fault.
Keywords: Magnetotelluric    Anisotropic inversion    East Kunlun Fault    Qaidam Basin    Altyn Fault    
0 引言

地球介质的电性特征受裂缝以及矿物的定向排列、含水量、部分熔融等的影响(Wannamaker,2005Wang et al., 2006Yoshino et al., 2006Dai and Karato, 2014),在不同层次上(地壳以及上地幔)都可能表现出各向异性.而分辨以及解释包含各向异性效应的数据也一直以来是大地电磁测深(MT)方法领域的难点,在针对大地电磁实际数据的各向异性研究中,大多数是应用维性工具、正演模拟等定性分析手段去解释各向异性结构(Wannamaker,2005Martí,2014).在较早的研究中,人们观察到大地电磁阻抗张量非对角元素的相位曲线有时在长周期存在明显的差异,认为这种相位分离现象是地球深部电各向异性的表现,并有研究将MT数据中的相位分离现象与地震学中的横波分离现象进行对比分析(Kurtz et al., 1993Simpson,2001Bahr and Simpson, 2002Leibecker et al., 2002Eaton et al., 2004; Eaton and Jones, 2006Padilha et al., 2006).然而,随着模型研究的深入,Heise等(2006)指出相位分离现象是由电导率的空间差异或梯度分布引起的,这种现象并不能指征电导率是各向同性的或各向异性的本质特征,因而,MT数据的相位分离与地震横波分离是完全不同类别的现象.相位分离的一种特殊现象是低频段的相位超象限现象(Chouteau and Tournerie, 2000),它可以由各向异性感应产生,并被作为各向异性结构可能存在的印迹(Heise and Pous, 2003Kumar and Manglik, 2012Martí,2014Liddell et al., 2016).但是,相位超象限现象同样也可以通过三维各向同性结构来模拟(Lezaeta and Haak, 2003Weckmann et al., 2003Ichihara and Mogi, 2009Ichihara et al., 2013Xiao et al., 2016邵贵航和肖骑彬,2016).另外,理论模型研究表明,含有本征各向异性体的模型响应可以用较为复杂的各向同性结构(如宏观各向异性结构)进行模拟(Eisel and Haak, 1999Heise and Pous, 2001Martí,2014),这种等效性主要是因为在相应的探测深度上MT方法分辨率不足而导致的(Weidelt,1999),这也提供了利用特殊的等效各向同性结构解释各向异性的可能性.但是,这种等效性的存在也会使得我们通过各向同性反演获得的电阻率结构发生明显的扭曲(Miensopust and Jones, 2011Löwer and Junge, 2017),从而对研究区的结构产生误导.此外,各向异性结构并不是都存在等效的各向同性结构,例如,简单的均匀各向异性半空间的响应无法通过构建等效的各向同性结构来获得(Martí,2014).这些研究结果表明,探测各向异性是极其复杂的过程,仅仅依靠资料的定性分析以及正演模拟等手段来判断地下深部的电各向异性结构还存在着较大的争议,而为了寻找更可靠的各向异性模型来合理解释地下可能存在的结构,我们有必要完成相应的反演算法,同时积累和总结实际数据应用中的经验和认识.

大地电磁各向异性正、反演方法的研究最早可追溯到20世纪60年代,接下来直到70年代末期关于各向异性的大地电磁一维正、反演和二维正演也一直是研究热点.Reddy和Rankin(1971)讨论了层状模型中,倾斜各向异性对大地电磁响应的影响;Reddy和Rankin(1975)利用有限元算法最早实现了二维各向异性模型响应的计算;Abramovici(1974)完成了一维任意各向异性的正演算法;Abramovici和Shoham(1977)应用矩阵广义逆的方法完成针对一维垂直各向异性模型的反演.而在80年代至90年代中期关于大地电磁各向异性问题的研究成果发表量明显减少(徐世浙和赵生凯,1985).90年代后期至今,Pek和Vener(1997)利用有限差分法实现了计算二维任意各向异性模型的正演算法;Yin(2003)讨论了大地电磁各向异性反演固有的多解性;Pek和Santos(2006)将标准的一维奥卡姆反演拓展到各向异性情况;Pek等(2011)完成二维任意各向异性反演算法,并研究了不同参数对反演的影响.随着计算机技术的快速革新,关于各向异性的大地电磁正、反演方法再次成为研究热点(杨长福,1997Martinelli and Osella, 1997Weidelt,1999Pek and Santos, 2002Li,2002Li et al., 2003Li and Pek, 2008杨长福等,2005Baba et al., 2006Mandolesi and Jones, 2012Chen and Weckmann, 2012Plotkin,2012秦林江和杨长福,2012胡祥云等,2013霍光谱等, 2012, 2015熊彬等,2013Qin et al., 2013Yang et al., 2015Montahaei and Oskooi, 2014Löwer and Junge, 2017Cao et al., 2018Kong et al., 2018Liu et al., 2018Yu et al., 2018Xiao et al., 2020喻国等,2019).

然而,虽然基于各向异性介质的大地电磁正演技术已经可以进行复杂的模型计算,但是各向异性反演在实际数据中的应用依然进展缓慢,这其中一个关键因素在于对各向异性参数的恢复上存在很大的局限(Yin,2003Pek et al., 2011Chen and Weckmann, 2012).Baba等(2006)首先将二维各向异性反演应用在穿过东太平洋南部洋中脊的大地电磁剖面的解释中,算法在Rodi和Mackie(2001)的二维各向同性反演代码的基础上修改完成了针对三个主轴电阻率的垂直各向异性反演.Baba等(2006)在对反演结果的分析中指出,相比于各向同性的反演结果,各向异性反演并没有明显的改善对数据的拟合程度,因此无法单独依靠大地电磁数据分析来确定各向异性模型为最优模型;然而,通过其他地球物理观测资料与大地电磁的各向异性模型之间的相互印证,才表明得到的各向异性模型是可靠的结果.Pek和Santos(2006)将一维各向异性反演得到的模型作为二维结构试错过程中的初始模型.Pek等(2011)也将他们完成的二维各向异性反演代码应用于葡萄牙南部采集的大地电磁数据,他们认为二维各向异性反演结果有可能受到观测数据中存在的非二维效应的影响,但没有对数据作深入探究;其中,反演待解模型参数为三个主轴电阻率、各向异性走向角以及倾斜角,但反演最终模型整体上呈现较强的各向同性.林长佑等(2004)将一维大地电磁垂直各向异性反演用在天祝—永登地区实际资料的定性解释中.杨长福等(2005)开发了二维大地电磁垂直各向异性反演算法(针对三个主轴电阻率),认为当电导率沿垂向和倾向表现相同的各向异性时,可以分别用TE和TM极化反演结果来代表两个主轴方向电阻率模型.霍光谱等(2012)基于马夸特方法,完成了大地电磁层状任意各向异性的视电阻率和相位反演,并指出需要将数据揭示的电各向同性与各向异性现象区分,否则无法得到准确合理的地层信息,甚至会出现自相矛盾的情况.秦林江和杨长福(2012)应用广义逆矩阵方法完成一维大地电磁各向异性反演算法,指出对实测资料反演时,如果数据是受各向异性影响的,那么必须用各向异性的处理方法才能得到更客观真实的结果.Cao等(2018)应用有限内存拟牛顿方法完成了针对三个主轴电阻率的大地电磁反演算法,但还没在实际中应用,并且目前还没有其他的大地电磁三维各向异性反演算法完成.通过以上国内外大地电磁各向异性反演应用的情况可以知道,大地电磁各向异性反演算法还不完善,一维各向异性反演在实际应用中存在明显的局限性;二维各向异性反演的实例较少;三维反演目前还没有对数据应用的实例.因此,合理全面地认识和理解地下存在的电各向异性结构,需要进一步发展各向异性反演算法和更多地开展在大地电磁实际数据中的应用.

现今基于三维各向同性介质的反演技术近年来发展十分迅速并已经在实际资料解释中广泛应用(Siripunvaraporn et al., 2005a, bEgbert and Kelbert, 2012Kelbert et al., 2014).Kelbert等(2014)完成的ModEM程序将反演过程封装成不同功能的单独模块,为MT方法的推广应用创造了条件.本文在已经实现的二维任意各向异性大地电磁正演算法的基础上(Yu et al., 2018),构建基于各向异性介质的反演目标函数,实现对目标函数值及其梯度的计算.参考ModEM的非线性共轭梯度法(NLCG)反演模块,实现对各向异性介质的二维大地电磁反演,通过理论模型验证反演结果的稳定性.在此基础上,我们对东昆仑—柴达木盆地过渡带一条MT剖面数据进行各向异性反演解释,该剖面所在的测区的地下地震各向异性背景丰富,但还没有相关电各向异性结构的研究.我们结合该地区以往的地震各向异性研究成果以及本文各向异性反演得到的结果来判断深部各向异性结构的存在性,并进一步分析各向异性的成因及其构造意义.

1 正演模拟

在大地电磁方法中,电场E和磁场H满足麦克斯韦方程组,通过基本的数学代数运算可以得到电场矢量方程或磁场矢量方程.在二维各向异性介质中,因为似稳假设的Maxwell方程组无法再解耦成单独的TE或TM模式,因此这里我们采用一般形式的电场矢量方程

方程组具体表达如下:

(1)

假设时间因子为,磁导率采用真空磁导率μ=4π×10-7H·m-1ω表示角频率;σ表示电导率张量,是一个对称正定的3×3矩阵.因此,σ有6个独立参数,并且可以通过3个主轴电导率σxσyσz以及3个各向异性角——走向角α,倾斜角β,偏向角γ来表征(Martí,2014),

(2)

其中,Rz表示沿Z轴旋转的旋转矩阵,Rx表示沿X轴旋转的旋转矩阵,T表示矩阵转置.

在右旋笛卡尔坐标系下,对模型进行矩形网格剖分,采用经典的交错网格格式来定义场的位置(Yee,1966).对方程(1)进行数值离散,得到线性方程组

(3)

其中,A为系数矩阵,B为右端项,E为所有的未知电场.方程(3)采用直接法求解器PARDISO(Schenk and Gärtner,2004Kuzmin et al., 2013)进行求解.在求得电场值E后,构建插值向量egihgi(下标i=xy),统一用解向量E计算地表指定点的电场egiE和磁场hgiE(Newman and Alumbaugh, 2000).

计算相互垂直的两组初始水平极化条件下的电场E1E2后,求得大地电磁响应函数,以在某一观测点的阻抗张量元素Zxy为例,

(4)

2 反演理论

大地电磁各向异性反演问题,可以归结为求解以下目标函数的极小值问题(Tikhonov and Arsenin, 1977):

(5)

相比较于各向同性反演,构建的目标函数中除了包含数据拟合项和模型粗糙度项外,还引入了各向异性惩罚项(Pain et al., 2003Pek and Santos, 2006Pek et al., 2011Chen and Weckmann, 2012).

ϕd为目标函数中的数据拟合项,dnobs表示观测数据;dnrsp表示模型的响应;εn为观测数据误差;N为观测数据总数;ϕs为模型粗糙度项,λs为与模型粗糙度相关的正则化参数,m0表示初始模型向量,m为反演迭代过程中某一次的模型向量,Cm为模型协方差矩阵,模型协方差矩阵Cm的引入可以在反演中灵活地加入先验信息(Siripunvaraporn and Egbert, 2000Siripunvaraporn et al., 2005aEgbert and Kelbert, 2012Kelbert et al., 2008, 2014);φa为模型各向异性惩罚项,λa为对应的正则化参数,K是约束各向异性模型中主轴电导率的矩阵,其作用是使主轴电导率之间的差异最小(Pain et al., 2003),具体表达式如下:

其一般形式为

其中a12a13a23为正数.

在ModEM中,引入新的模型向量

使得

在反演中求取目标函数的梯度时,针对的是模型向量 ,避免了对Cm-1/2的求逆过程.在实际反演过程中,不用对Cm-1/2显式计算,只需要求取Cm-1/2与向量的乘积.

(6)

(7)

(8)

这里,Δdn=(dnobs-dnrsp)/εn2,其中上标*表示复数的共轭,ml为模型m中的第l个参数,nm为总的模型参数.假设模型网格数为ng,每个网格拟反演的各向异性参数个数为na,那么nm=ng×na.反演中,各向异性参数分别定义为主轴电导率的自然对数和各向异性角的弧度值.Yin(2003)指出一维大地电磁各向异性反演无法分辨垂直主电阻率ρz和各向异性倾角β.Pek等(2011)通过模型研究提出,二维大地电磁各向异性反演也几乎不能还原垂直主电阻率ρz,而且类似于一维反演,倾角β也不能被分辨.当我们在反演中不考虑倾角β时,走向角α和偏向角γ具有相同的旋转对称轴Z轴,偏向角γ也等同于走向角α,公式(2)可写为

(9)

该式表示依次绕Z轴旋转α角和γ角,αγ表示相同的旋转角,所以可以将走向角α和偏向角γ合并为同一个参数α,如下:

(10)

因此,综合考虑大地电磁各向异性反演的特点,我们在反演中设置待求解的模型参数只包括ρxρyρzα.

(6)—(8)式中,计算量主要集中在(8)式中数据对模型的偏导数计算上.这里NNp(频率点数)、Ns(测点数)和Nr(每个测点的响应函数数量)的乘积.取d包含全部阻抗张量元素(ZxxZxyZyxZyy),那么(8)式可以展开为

(11)

根据(11)式,对于某一个频率点,在求解数据对模型的偏导数时,它包含一次正演求解响应数据dnrsp,以及一次拟正演过程求解jLiTA-1,其中i代表频率点,j代表不同极化模式,这里的一次正演和拟正演均包含两种极化模式的计算.梯度的求解过程与Newman和Alumbaugh(2000)文中的处理方式一致,其他相关的细节这里不再重复赘述.

当求取了目标函数值及其梯度后,参考ModEM中的NLCG反演模块,我们实现了对各向异性模型的二维大地电磁反演.在反演中对模型向量进行更新:

其中,k代表迭代次数,α为当前搜索步长,p为模型搜索方向.反演计算中实际需要调整的主要参数包括:初始正则化参数λsλa,正则化参数更新因子k,初始线搜索步长α0,数据误差门槛,以及与模型协方差相关的光滑参数.

在反演迭代的计算过程中,当数据拟合均方根误差(RMS)减小量低于某个阀值时,正则化参数更新为λs/kλa/k.而初始线搜索步长α0在这里采取人为设定的方式,根据大量反演算例的经验,初始线搜索步长的大小会显著影响反演的收敛过程,当步长过大,模型在初始阶段易发生较大调整而产生冗余构造,在反演后期难以校正,并且这个影响在各向异性反演中更为突出.因此,在实际操作中,我们需要尽量避免在首次迭代中出现很大的模型参数调整.

3 理论模型的反演检验

模型1(Md1)为一包含三个各向异性体的二维模型(图 1a).模型背景电阻率为100 Ωm,在8.65 km以上包含两个不同规模的各向异性体A1和A2,它们的顶部埋深都是0.135 km,宽度同为3.2 km.两个异常体之间相距0.8 km,底部埋深分别为8.65 km和1.36 km.两个异常体具有相同的各向异性参数:ρx=1 Ωm、ρy=200 Ωm、ρz=1 Ωm、α=30°、β=0°、γ=0°.异常体A1下伏一个各向异性层,层厚14.42 km,模型参数为ρx=5 Ωm、ρy=50 Ωm、ρz=5 Ωm、α=120°、β=0、γ=0.将该模型在水平方向剖分为46个网格,地表向下剖分为39个网格,地表定义26个观测点,计算周期范围0.015627~15627 s内29个周期点的阻抗响应数据.

图 1 二维各向异性理论模型及其反演结果 (a)二维理论模型Md1示意图及其视电阻率和相位响应拟断面图;图中小倒三角表示地表观测点位置,A1、A2、A3表示各向异性体.(b)取不同正则化参数反演得到的各测点RMS值.红色实点表示λs=0.001和λa=0.001时的拟合RMS值;蓝色空心点表示λs=0.1和λa=0.1时的拟合RMS值.(c)—(f)当取正则化参数λs=0.001和λa=0.001时反演得到的模型与理论模型对比,黑框为理论模型中异常体的范围;(c)—(f)依次为模型参数ρxρyρzα的剖面图.(g)—(j)当取正则化参数λs=0.1和λa=0.1时反演得到的模型与理论模型对比. Fig. 1 Results obtained by two-dimensional anisotropic inversion of simulated data generated from Md1 (a) Complete parameters of Md1 and apparent resistivity and phase response pseudosection generated by Md1. The inverted triangles denote the observation locations; A1, A2 and A3 represent anisotropic anomalies. (b) Comparison of RMS value at each observation point under different regularization parameter sets. The red solid points denote the misfit RMS values with λs=0.001 and λa=0.001. The blue empty circles denote the misfit RMS values with λs=0.1 and λa=0.1. (c)—(f) Comparison of final inverted and true models of anisotropic parameter, ρx, ρy, ρz and α, respectively. The regularization parameters are set as λs=0.001 and λa=0.001. The black boxes show anomalous zones. (g)—(j) Comparison of final inverted and true models of anisotropic parameter, ρx, ρy, ρz and α, respectively. The regularization parameters are set as λs=0.1 and λa=0.1.

在对Md1的正演响应数据进行反演时,采用了全部阻抗张量数据,数据ZxxZyy的误差给定为;数据ZxyZyx的误差则分别为各自的幅值的5%.反演模型网格和正演网格一致,初始模型电阻率为100 Ωm.模型的光滑参数(τδyδz)分别为10,0.1和0.1;反演中待求解的模型参数包括ρxρyρzα.初始线搜索步长α0=20,初始正则化参数λs分别取为0.001,0.01,0.1;λa分别选为0.001、0.01、0.1、1.0、10.0、100.0和1000.0进行计算并对比反演结果.图 2表示不同正则因子选择下的最终RMS值,其中,λs=0.1,λa=0.001的参数组合在经过122次反演迭代后得到了最小的RMS值1.1137,此外,在粗糙度正则因子λs初始取值范围为0.001到0.1,以及各向异性惩罚正则因子取值范围为0.001到0.1时,这两个正则因子区间的参数组合得到的最终RMS值差别并不明显,都非常接近1.0,我们在图 1中展示了正则因子选择这两个取值范围边值情况的最终模型(图 1c—fλs=0.001和λa=0.001的情况,图 1g—jλs=0.1和λa=0.1的情况).

图 2 不同正则化参数组合下对应数据拟合RMS值的变化 Fig. 2 The RMS values calculated by different regularization parameter combinations

图 1b所示,两次反演中拟合RMS偏差相对较大的点是位于异常体A1和A2之间的测点13和14,这表明二维电性结构横向的剧烈变化与各向异性效应一起增加了数据拟合的难度.从模型恢复程度来看(图 1c—j),ρz出现了一些细微的构造,但依然无法实现对异常结构的有效分辨,这主要是因为ρz在二维反演中几乎无法还原,其结构的变化主要是受目标函数中各向异性约束项的影响(Pek et al., 2011).

图 1cd表明,在λs=λa=0.001时,反演恢复A1和A2区域的各向异性参数特点表现为,主轴电阻率ρx为200 Ωm左右的相对高阻ρmaxρy为1 Ωm左右的相对低阻ρmin,这与真实模型中水平主轴电阻率的分布特点(ρx=1 Ωm表现为低阻ρminρy=200 Ωm表现为高阻ρmax)相反,而在图 1f中,A1和A2区域恢复的走向角α约为-60°,与真实的α=30°相差90°.从单独的各向异性参数来看,还原的参数与真实参数不同,但是反演还原的参数组(ρx=ρmaxρy=ρminα=-60°)与真实的参数组(ρx=ρminρy=ρmaxα=30°)代表的是相同的各向异性体,如(12)式所示,两组参数表述的是相同的电导率张量,因此,初始λs=λa=0.001计算的反演结果图 1cdf基本还原了A1和A2区域真实的各向异性异常体,还原的异常体参数与真实异常参数等效.同样,对图 1ghj也可以得出相同的结论.

(12)

对于上部异常体A1和A2的尺度恢复,异常中心越集中在浅部越有利于准确地把握异常体的尺度,这也符合大地电磁方法浅部分辨率高的特点(Oldenburg et al., 2013),走向方位角也能较好恢复异常体的尺度.对于最下方层状异常体,两次反演结果都不能实现对下部异常体的恢复,但是在反演的电性结构中,模型浅部两侧出现了相对异常结构(如图 1c中的C1和C2、图 1d中的B1和B2).结合测点RMS拟合偏差值观察,两侧的测点拟合程度高,因此,两侧的冗余结构是反演的多解性所导致的.

4 东昆仑—柴达木过渡带岩石圈电各向异性研究 4.1 研究区域各向异性背景

东昆仑—柴达木过渡带东西跨度约1000 km,横亘于青藏高原的北部(图 3a).它们之间以柴达木南缘断裂(SQDF)相隔(图 3b),东昆仑造山带大陆地壳主要形成于古元古代晚期,自显生宙以来,东昆仑经历了多期洋陆演化,同时,东昆仑造山带也是一条巨型构造岩浆岩带(莫宣学等,2007).左行走滑的东昆仑断裂带与东昆仑山相伴而生,它继承了早期缝合带特征,并在新生代重新活动.该断裂带现今的走滑速率在中段约11 mm·a-1,并在印度—亚洲板块汇聚过程中吸收北东向挤压而产生缩短变形,对青藏高原的形成及演化发挥了重要作用(Tapponnier et al., 2001).

图 3 东昆仑—柴达木盆地过渡带和邻区构造背景以及前期大地电磁测点分布图 (a)研究区大地构造背景,图中黄色粗箭头表示区域应力方向;ATF:阿尔金断裂带;KLF:东昆仑断裂带;HYF:海原断裂带;RRF:红河断裂带;XXF:鲜水河—小江断裂带;JLF:嘉黎断裂带.(b)东昆仑—柴达木盆地地质图,图中三角形表示大地电磁测深点;其中红色为本次研究的剖面;SQDF:柴达木南缘断裂.(c)沿L15剖面早期二维各向同性反演结果,剖面方向从SW至NE(Xiao et al., 2018). 断裂构造参考邓起东等(2003). 地质图据1:250万中国地质图,在线浏览地址: http://www.ngac.org.cn/. Fig. 3 Tectonic setting of the adjacent area of east Kunlun to Qaidam Basin tranistion zone and distribution of previous MT observation points (a) Tectonic setting of surveying area. The yellow thick arrows denote the direction of the regional stress. ATF: Altyn Tagh Fault; KLF: East Kunlun Fault; HYF: Haiyuan Fault; RRF: Red River Fault; XXF: Xianshuihe-Xiaojiang Fault; JLF: Jiali Fault. (b) East Kunlun-Qaidam Basin geological map. The triangles display MT observation sites, the red ones are the profile we study in this paper. SQDF: South Qaidam Fault. (c) Previous isotropic inversion result along profile L15. The direction of profile is from southwest to northeast (Xiao et al., 2018). The fault structure from Deng et al. (2003). The geology data of 1:2.5 million scale can be browsed at the website: http://www.ngac.org.cn/.

在青藏高原中部,观察到的地震各向异性可以通过中下地壳的弱流动层解释(Shapiro et al., 2004),同时该软弱流动层在电性上也表现出较强的各向异性(Le Pape et al., 2015),但关于它是否越过东昆仑进入到柴达木盆地之下,这仍存在争议(Zhu and Helmberger, 1998Unsworth et al., 2004Shi et al., 2009Karplus et al., 2011Le Pape et al., 2012Wei et al., 2014).如图 3b,我们先后跨过东昆仑造山带与柴达木盆地之间的过渡带的不同位置采集了5条大地电磁测深剖面的数据,并在前期的研究工作中发现,在剖面L15中祁漫塔格下方的上地幔顶部和下地壳底部存在低阻异常体(图 3c中C1和C2异常体)(Xiao et al., 2013, 2016, 2018).同时,该上地幔顶部低阻异常体对应的地面向南约150 km处的地表存在新生代的火山岩(Jolivet et al., 2003Wang et al., 2005),Wang等(2005)认为这是增厚的松潘—甘孜地壳下部发生部分熔融的证据.

在这一地区的地震快波偏振方向显示,从阿尔金断裂带到柴达木盆地,快波偏振方向从与阿尔金断裂带走向平行向近东西向转变,发生了顺时针旋转,而少部分的快波偏振方向与地表山脉结构走向不一致,其中,在祁漫塔格内部的一个测点显示最大左行剪切方向为北东向,与祁漫塔格山脉走向不同(Herquel et al., 1999Soto et al., 2012),这也为青藏高原北部的动力学过程提供了新的认识.Le Pape等(2012)在重新获取大地电磁剖面的各向异性模型时,发现各向异性的存在使东昆仑与柴达木接触带的地壳结构发生了明显的变化,为松潘—甘孜中下地壳低阻流动层继续向北穿透提供了新的证据.这一节中,我们利用新的大地电磁各向异性反演程序对L15剖面数据进行重新反演,进一步分析在祁漫塔格与柴达木过渡部位低阻通道的各向异性特征,并在下一节的讨论中利用电各向异性参数的变化样式来分析青藏高原北部软物质流的存在及区域应力环境.

4.2 二维各向异性结构

对测线L15的数据进行二维各向异性反演前,需先将观测数据坐标系统与正演中的坐标系统匹配.在野外观测的测点都是按正南北-正东西方向接收水平电、磁场,分别与直角坐标系中的xy方向对应,而正演过程中,x方向指示结构走向,垂直于测线,y方向是沿测线展布方向.考虑到测线展布并非正南北或者正东西方向,而是垂直于地质构造走向的,因此,这里要将观测数据的坐标系统逆时针旋转55°,使数据的坐标系统的x轴与测线方向垂直(Xiao et al., 2016, 2017, 2018).反演中,我们选取了全阻抗张量的数据,数据误差门槛比例设置为5%,其中,数据ZxxZyy的误差门槛值设定为 ×0.05,数据ZxyZyx的误差门槛值则分别为各自幅值的5%.

反演的模型网格为84×240,84为横向网格剖分,240是地下分层数,模型加入地形的变化,初始模型电阻率为100 Ωm.初始线搜索步长α0=20.

首先,我们进行三组数据分类的各向同性反演:对ZxyZyx数据的反演、对单独Zxy分量的反演以及对单独Zyx分量的反演.通过各向异性反演程序实现各向同性反演,只需设定模型参数ρx=ρy=ρzα=β=γ=0°,这样处理后,目标函数中各向异性惩罚项ϕa始终为0,不再对反演过程产生作用.然后,选取全部阻抗张量数据(ZxxZxyZyxZyy),进行反演,待还原的参数仍为ρxρyρzα.

图 4aZxyZyx数据的各向同性反演结果来看,模型的结构特征与Xiao等(2018)中各向同性反演的二维电性结构相似,在剖面最南端的上地幔顶部出现了低阻异常C1,并向北延伸,与柴达木盆地下地壳层次的低阻体C2相连,形成低阻带iso1(图 4a),与图 3c中的异常体C1和C2对应(Xiao et al., 2018).然而,这种结构在单独的Zxy分量数据以及单独Zyx分量数据的反演结果中并不明显(图 4bc).

图 4 剖面L15在不同正则化参数组合下的各向异性反演模型 (a)—(c)分别为对阻抗ZxyZyx分量、单独的Zxy分量以及单独的Zyx分量的二维各向同性反演结果;(d)—(g)正则化参数组合为λs= 0.1和λa=0.1的各向异性反演结果;(h)—(k)正则化参数组合为λs=0.1和λa=1.0的各向异性反演结果;(l)—(o)正则化参数组合为λs=0.1和λa=10.0的各向异性反演结果;(p)—(s)正则化参数组合为λs=0.1和λa=100.0的各向异性反演结果;(t)—(w)正则化参数组合为λs=0.1和λa=1000.0的各向异性反演结果. Fig. 4 The inversion models of profile L15 under different regularization parameter combinations (a)—(c) are the isotropic inversion results from impedance components Zxy and Zyx, component Zxy, and component Zyx; (d)—(g) are the anisotropic inversion results under regularization parameters λs=0.1 and λa=0.1; (h)—(k) are the anisotropic inversion results under regularization parameters λs=0.1 and λa=1.0; (l)—(o) are the anisotropic inversion results under regularization parameters λs=0.1 and λa=10.0; (p)—(s) are the anisotropic inversion results under regularization parameters λs=0.1 and λa=100.0; (t)—(w) are the anisotropic inversion results under regularization parameters λs=0.1 and λa=1000.0.

其次,我们在进行阻抗全张量的各向异性反演计算并测试了大量不同正则因子组合后,选择粗糙度正则因子为λs=0.1,其得到的各向异性模型与各向同性模型的结构特征更接近.另外,我们在图 4图 5中展示了λa分别选取为0.1、1.0、10.0、100.0、1000.0时,反演得到最终模型和RMS值的对比.图 4中,λa=0.1、1.0、10.0、100.0以及1000.0的结果模型对比可以看出,λa在0.1到10.0范围内得到的各向异性结构比较相似,随着λa增大,模型中电各向异性程度降低,当λa达到100.0时,各向异性结构不再明显,在λa=1000.0时,ρxρy以及ρz的电阻率结构几乎一致,并且没有明显的各向异性走向角α的结构变化,模型接近于各向同性结构;从图 5可以看出,当λa从小到大变化时,最终RMS值在λa=1.0时达到最小,然后呈上升趋势,其中在λa=10.0时,数据拟合最终RMS开始显著增大,对应的目标函数值也明显变大,各向异性惩罚项的目标函数值明显变小.

图 5 测线L15在不同各向异性正则化参数下的反演结果统计,其中粗糙度正则因子固定为λs=0.1 (a) RMS值;(b) 总目标函数值;(c) 粗糙度项目标函数值;(d) 各向异性约束项目标函数值. Fig. 5 The inversion results of profile L15 versus different anisotropic regularization parameters. The roughness regularization parameter is λs=0.1 (a) Final RMS values; (b) Full objective function values; (c) Roughness objective function values; (d) Anisotropic objective function values.

根据对图 4图 5的分析,我们选择RMS值最小即数据拟合最好的结果作为各向异性反演的最终模型(图 4h—k),正则因子组合为λs=0.1和λa=1.0.图 6a—d也同为该各向异性模型,ρxρy分量的电性结构具有较大的差异,表现出明显的各向异性特征,其中,ρx剖面中的低阻体C1所在位置对应ρy剖面中的相对高阻体R1,C1核心区域各向异性走向角α约为45°;ρx剖面中相对高阻体R2所在位置对应ρy剖面中的低阻体C2,C2核心区域各向异性走向方位角α大约为-45°.这里x轴垂直于测线并指向西北,根据两个异常体的各向异性走向角的旋转效应,两者的低阻电阻率分量C1和C2所在轴的空间方位是在与x轴指向方向夹角45°(顺时针旋转45°)的轴线上,即两个各向异性异常体的最大和最小主轴电阻率的空间展布一致.因此,可以将图 6a中C1和R2的位置归为一个各向异性带ani1,它的模型参数可以表示为ρx~3 Ωm,ρy~300 Ωm,α~45°.因为前文提到,已经将观测系统逆时针旋转55°至测线方向,那么再将观测系统顺时针旋转55°至正南北-正东西坐标系统,就可以知道ani1区域里低阻主轴电阻率所在轴的方位呈北偏西10°,如图 7所示.

图 6 验证剖面L15深部各向异性结构存在性,结果均为利用全部阻抗张量数据对模型参数ρxρyρzα进行各向异性反演获得 (a)—(d)初始模型为100 Ωm半空间;(e)—(h)初始模型是在图 4a模型的基础上,图 5a深度10.9 km以上保持不变,将图 5a深度10.9 km以下的区域用100 Ωm均匀半空间替代. Fig. 6 Examining the deep anisotropic structure of the profile L15. The results are obtained from inverting full impedance elements data for recovering parameter ρx, ρy, ρz and α (a)—(d) The initial model is a 100 Ωm half-space. (e)—(f) The initial model is modified from the model in Fig. 4a. The area above 10.9 km depth remains unchanged in Fig. 5a. The area below 10.9 km depth of model in Fig. 5a is replaced with 100 Ωm half-space.
图 7 图 6中ani1区域的各向异性异常水平主轴电阻率的空间展布 实线表示正南北-正东西的坐标系统;虚线表示沿测线方向的坐标系统,其中y轴是沿测线展布方向;点虚线表示ani1区域各向异性体的水平主轴电阻率ρxρy的空间展布,最终低阻主轴电阻ρx在北偏西10°的轴线上. Fig. 7 Spatial distribution of the horizontal principal >resistivities of the anomaly in the ani1 area of Fig. 6 The solid lines denote the coordinate axes in west-east and south-north directions; the dotted lines denote the coordinate system of which the y axis is in the direction along the profile; the dot-dashed lines denote the spatial distribution of the horizontal principal resistivity ρx and ρy of the anomaly in ani1 area, the minimum principal resistivity is in the north to west 10° axis.

ani1与图 4a中iso1的形态接近,前者的整体埋深略浅于后者.另外,在柴达木盆地中心部位的中下地壳,也存在明显的各向异性体ani2(图 6a中C3所在区域),它的模型参数可以表示为ρx~3 Ωm,ρy~500 Ωm,α~ 0°.图 6aρz分量的电性结构在浅层次有所变化,但变化程度小,这主要也与前文提到的原因相关,即二维反演几乎无法分辨ρz.

在数据拟合RMS方面,二维各向同性反演的总体RMS值要小于各向异性反演的结果,这与各向异性反演中考虑了全部阻抗张量元素有关.如图 8a所示,对全阻抗张量各向异性反演的结果(图 6a—d)另外计算每个测点对阻抗张量非对角元ZxyZyx的RMS,其拟合程度与针对非对角元ZxyZyx的各向同性反演(图 4a)的RMS拟合程度非常接近,这一定程度说明,各向异性反演过程在对ZxyZyx充分拟合的同时,也考虑到了主对角元素ZxxZyy的信息.从图 8b来看,各向异性反演结果中,对ZxxZyy的整体拟合效果不如ZxyZyx,但是在测线两端和中间部位,即测点1—3、8—14以及18—27的区间,二者拟合效果相当,而且这三个位置大致与图 6a—d中的各向异性结构ani1和ani2的位置对应,因此,ZxxZyy的加入有利于反演过程中还原各向异性结构.

图 8 不同二维反演结果的测点RMS对比 (a) 各向同性与各向异性反演结果中次对角元素拟合RMS的对比;三角形:图 4a各向同性模型中对ZxyZyx计算的RMS值;圆形:图 6a—d各向异性模型中对ZxyZyx计算的RMS值. (b)图 6a—d各向异性反演结果中,对不同阻抗张量元素组合的拟合RMS对比.方形:对ZxyZyx计算的RMS值;菱形:对ZxxZyy计算的RMS值. Fig. 8 Comparison between the data misfit RMS values of different two-dimensional inversion at each site (a) Comparison between the off-diagonal elements data misfit RMS values of isotropic and anisotropic inversion. Triangle shows the RMS value of impedance elements Zxy and Zyx of the model in Fig. 4a. Circle shows the RMS value of impedance elements Zxy and Zyx of the model in Figs. 6a—d. (b) Comparison between the data misfit RMS values with different combinations of impedance elements of anisotropic model in Figs. 6a—d. Square shows the RMS value of impedance elements Zxy and Zyx. Diamond shows the RMS value of impedance elements Zxx and Zyy.

为进一步验证各向异性结构的稳健性,我们以图 4a中的各向同性结构模型为基础,将该模型大于10.9 km深度的电阻率重新置换为100 Ωm,形成新的初始模型,采用与图 6a—d相同的各向异性反演参数设置.在反演开始时,数据拟合RMS为3.5369,明显高于图 4a的1.3860.在经过37次迭代后,RMS降为1.7133.反演结果如图 6e—h所示,ρxρy剖面中还原的ani1和ani2区域的各向异性体与图 6a—d同样位置的异常特征基本一致.因此,这个结果说明各向异性带ani1以及ani2能够被观测数据所分辨,并且不能被各向同性体所替代.

5 讨论

(1) Pek等(2011)在二维各向异性反演研究中指出,因为ρz对大地电磁正演响应的影响非常微弱,ρz分量也几乎不能被分辨,在数据中加入倾子等与垂直磁场相关的传输函数也无法改善对这个分量的还原情况,各向异性反演中成像的ρz结构主要受目标函数中各向异性约束项的控制.本文中二维理论模型的反演算例对ρz还原也很不理想,因此,即使我们在实际的反演计算中得到较满意的数据拟合程度,但要进一步还原和解释ρz分量结构还需寻找其他对ρz敏感的地球物理资料或者地质结构信息来支持.

(2) 电各向异性是通过电导率张量来表征,而电导率张量可以通过不同的各向异性参数组合来形成,例如,方位各向异性参数组(ρx=10 Ωm,ρy=1000 Ωm,α=-45°)和(ρx=1000 Ωm,ρy=10 Ωm,α=45°)实际上代表的是相同的各向异性结构,这可能会导致在反演中,我们还原出了正确的各向异性结构,但是却由等效的各向异性参数组表述.在对二维理论模型的还原结果中,也出现了这种各向异性参数等效效应.因此,在对各向异性结构进行解释时,要将各向异性参数当作一个整体看待,避免结构认识上的偏差.

(3) 在各向异性反演中,浅部冗余结构会压制对深部各向异性体的分辨.在对理论二维各向异性模型的恢复中,可以观察到,随深度增加,对各向异性体的分辨能力逐渐下降,很难得到完整的恢复,同时,浅部会呈现出冗余结构,这种浅层冗余构造同样也是各向同性反演中的难题,它的出现主要是因为反演问题的多解性,而考虑各向异性又显著提升了多解性的程度,因此,通过压制浅部冗余结构来提升分辨深部异常体的能力,是将来各向异性反演中需要深入研究的方向.

(4) 对剖面L15的观测数据重新进行大地电磁各向异性反演,发现在祁漫塔格山下方50~100 km的区域存在方位各向异性体ani1.该异常体的水平低阻主轴电阻率空间展布在北偏西10°的轴线上,水平高阻主轴电阻率空间展布在北偏东80°的轴线上.在剖面L15邻近区域已有的远震宽频横波分裂研究结果(Herquel et al., 1995, 1999Soto et al., 2012Wu et al., 2015)显示,在阿尔金断裂、昆仑山断裂以及柴达木盆地周边,测量的大部分快波偏振方向与主要的地质构造走向相吻合(McNamara et al., 1994Guilbert et al., 1996).已有实例表明,大多数的快波偏振方向是沿着山脉和其分支走向的(Silver,1996),而快波偏振方向可以用来指示剪切带方向,这说明大多数深部剪切带与地表构造带可能受到相同变形机制的控制,少数快波偏振方向与地表构造不一致则暗示可能两者受到更复杂的应力环境影响.

关于地震各向异性与电各向异性间的联系,现在仍然存在争议,当测量的快波偏振方向与电各向异性主轴方向吻合时,有研究认为两者产生于相同的原因(Collins et al., 2006Frederiksen et al., 2006Padilha et al., 2006Carcione et al., 2007Jones et al., 2009);而也有一些研究实例认为这两者之间没有明显的关系(Hamilton et al., 2006).地震各向异性与电各向异性之间是否存在关系,这还需在将来通过更多的研究工作去寻找和验证.

剖面L15成像的祁漫塔格山下各向异性体ani1所在位置属于上地幔顶部,Caricchi等(2011)在解释上地幔电各向异性成因时指出,岩石圈下部的部分熔融通道在简单剪切环境下会产生电各向异性.其中,在剪切作用的垂直方向上由于存在应力梯度会降低熔融连通性,产生高阻;而沿剪切方向产生相对高连通性,表现出低阻.

而我们在前期对此次区域的三维各向同性反演研究中(Xiao et al., 2018Yu et al., 2020),将剖面L15下的上地幔低阻体解释为熔融的弱物质流,是源自松潘—甘孜下地壳流穿过祁漫塔格山脉向柴达木盆地下方运动的新证据,并且该处邻近区域位于多个剪切带交汇处(阿尔金断裂带,祁漫塔格造山带,昆仑山断裂带等).根据Caricchi等(2011)对地幔电各向异性成因的解释,图 6中的各向异性体ani1则可能是该弱物质流运动穿过祁漫塔格山脉下方时与剪切运动相互作用而形成.

这里ani1指示的剪切作用位于北偏西10°的轴线上,该轴线方位与Soto等(2012)在几乎相同位置的地震各向异性结果指示的北偏东45°左右的剪切带方位不一致,这很可能与两种各向异性结构在深度上的变化和跨度不一致相关;另外,ani1指示的剪切作用方位与祁漫塔格山脉走向也不一致,根据这两者间的关系进行推论:深部剪切作用方位与地表造山带地势的不一致暗示,祁漫塔格山脉地下应力环境并不单一,ani1指示的深部剪切作用可能是邻近区域剪切运动的综合体现,祁漫塔格在内的邻近区域深部剪切应力环境可能还受到阿尔金断裂带中部较大的左行走滑运动的影响.

综上,剖面L15反演还原的电各向异性结构与地震各向异性结果指示了不同的剪切作用方位,并且这两种各向异性的空间结构差异较大,所以无法确认两者成因之间的关系;此外,上地幔电各向异性体ani1可能是祁漫塔格山脉下方的弱物质流受到剪切作用的影响而产生的,其表征的深部剪切作用方位与祁漫塔格地表造山带走向不一致则暗示该处应力环境复杂,可能还受邻近阿尔金断裂带左行走滑的影响.

6 结论

基于大地电磁任意各向异性有限差分正演算法以及成熟的NLCG反演模块,本文实现了二维大地电磁各向异性反演,理论模型响应的反演算例也验证了算法的稳定性.

利用该反演程序,我们对穿过东昆仑祁漫塔格山—柴达木盆地的一条大地电磁剖面进行了重新反演.该剖面早期的各向同性结果显示,在祁漫塔格山脉下方的上地幔顶层存在连续的低阻带.在本文中,我们对该剖面进行各向异性反演,并在祁漫塔格下方相同层次的位置发现了明显的各向异性体,其水平低阻主轴电阻率所在轴的方位呈北偏西10°,水平高阻主轴电阻率所在轴的方位呈北偏东80°.我们认为这里的各向异性异常是上地幔弱物质流在祁漫塔格下方受剪切运动的作用而产生沿垂直于剪切作用的方向和沿剪切作用的方向的熔融物质导电性差异而形成的;也暗示这里的应力环境除了印度—欧亚板块碰撞带来的远程挤压效应外,还可能受到邻近阿尔金走滑断裂带左行剪切运动的影响.

致谢  感谢Gary Egbert教授提供了ModEM的反演代码.感谢三位匿名审稿人提出的宝贵意见,使本文得以进一步改善.
References
Abramovici F. 1974. The forward magnetotelluric problem for an inhomogeneous and anisotropic structure. Geophysics, 39(1): 56-68. DOI:10.1190/1.1440412
Abramovici F, Shoham Y. 1977. Inversion of anisotropic magnetotelluric data. Geophysical Journal International of the Royal Astronomical Society, 50(1): 55-74. DOI:10.1111/j.1365-246X.1977.tb01324.x
Baba K, Chave A D, Evans R L, et al. 2006. Mantle dynamics beneath the East Pacific Rise at 17°S: Insights from the Mantle Electromagnetic and Tomography (MELT) experiment. Journal of Geophysical Research: Solid Earth, 111(B2): B02101. DOI:10.1029/2004jb003598
Bahr K, Simpson F. 2002. Electrical anisotropy below slow and fast-moving plates: paleoflow in the upper mantle?. Science, 295(5558): 1270-1272. DOI:10.1126/science.1066161
Cao H, Wang K P, Wang T, et al. 2018. Three-dimensional magnetotelluric axial anisotropic forward modeling and inversion. Journal of Applied Geophysics, 153: 75-89. DOI:10.1016/j.jappgeo.2018.04.015
Carcione J M, Ursin B, Nordskag J I. 2007. Cross-property relations between electrical conductivity and the seismic velocity of rocks. Geophysics, 72(5): E193-E204. DOI:10.1190/1.2762224
Caricchi L, Gaillard F, Mecklenburgh J, et al. 2011. Experimental determination of electrical conductivity during deformation of melt-bearing olivine aggregates: implications for electrical anisotropy in the oceanic low velocity zone. Earth and Planetary Science Letters, 302(1-2): 81-94. DOI:10.1016/j.epsl.2010.11.041
Chen X, Weckmann U. 2012. Constraint inversion of 2D magnetotelluric data with anisotropic conductivities. //Extended Abstract 21st Workshop. Darwin, Australia.
Chouteau M, Tournerie B. 2000. Analysis of magnetotelluric data showing phase rolling out of quadrant (PROQ). //70th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts.
Collins J L, Everett M E, Johnson B. 2006. Detection of near-surface horizontal anisotropy in a weathered metamorphic schist at Llano Uplift (Texas) by transient electromagnetic induction. Physics of the Earth and Planetary Interiors, 158(2-4): 159-173. DOI:10.1016/j.pepi.2006.05.008
Dai L D, Karato S I. 2014. High and highly anisotropic electrical conductivity of the asthenosphere due to hydrogen diffusion in olivine. Earth and Planetary Science Letters, 408: 79-86. DOI:10.1016/j.epsl.2014.10.003
Deng Q D, Zhang P Z, Ran Y K, et al. 2003. Active tectonics and earthquake activities in China. Earth Science Frontiers (in Chinese), 10(S1): 66-73. DOI:10.3321/j.issn:1005-2321.2003.z1.012
Eaton D W, Jones A G, Ferguson I J. 2004. Lithospheric anisotropy structure inferred from collocated teleseismic and magnetotelluric observations: Great Slave Lake shear zone, Northern Canada. Geophysical Research Letters, 31(19): L19614. DOI:10.1029/2004GL020939
Eaton D W, Jones A. 2006. Tectonic fabric of the subcontinental lithosphere: evidence from seismic, magnetotelluric and mechanical anisotropy. Physics of the Earth and Planetary Interiors, 158(2-4): 85-91. DOI:10.1016/j.pepi.2006.05.005
Egbert G D, Kelbert A. 2012. Computational recipes for electromagnetic inverse problems. Geophysical Journal International, 189(1): 251-267. DOI:10.1111/j.1365-246X.2011.05347.x
Eisel M, Haak V. 1999. Macro-anisotropy of the electrical conductivity of the crust: a magnetotelluric study of the German Continental Deep Drilling site (KTB). Geophysical Journal International, 136(1): 109-122. DOI:10.1046/j.1365-246X.1999.00707.x
Frederiksen A W, Ferguson I J, Eaton D, et al. 2006. Mantle fabric at multiple scales across an Archean-Proterozoic boundary, Grenville Front, Canada. Physics of the Earth and Planetary Interiors, 158(2-4): 240-263. DOI:10.1016/j.pepi.2006.03.025
Guilbert J, Poupinet G, Jiang M. 1996. A study of azimuthal P residuals and shear-wave splitting across the Kunlun range (Northern Tibetan Plateau). Physics of the Earth and Planetary Interiors, 95(3-4): 167-174. DOI:10.1016/0031-9201(95)03120-0
Hamilton M P, Jones A G, Evans R L, et al. 2006. Electrical anisotropy of South African lithosphere compared with seismic anisotropy from shear-wave splitting analyses. Physics of the Earth and Planetary Interiors, 158(2-4): 226-239. DOI:10.1016/j.pepi.2006.03.027
Heise W, Pous J. 2001. Effects of anisotropy on the two-dimensional inversion procedure. Geophysical Journal International, 147(3): 610-621. DOI:10.1046/j.0956-540x.2001.01560.x
Heise W, Pous J. 2003. Anomalous phases exceeding 90° in magnetotellurics: anisotropic model studies and a field example. Geophysical Journal International, 155(1): 308-318. DOI:10.1046/j.1365-246X.2003.02050.x
Heise W, Caldwell T G, Bibby H M, et al. 2006. Anisotropy and phase splits in magnetotellurics. Physics of the Earth and Planetary Interiors, 158(2-4): 107-121. DOI:10.1016/j.pepi.2006.03.021
Herquel G, Wittlinger G, Guilbert J. 1995. Anisotropy and crustal thickness of Northern-Tibet. New Constraints for tectonic modelling. Geophysical Research Letters, 22(14): 1925-1928. DOI:10.1029/95GL01789
Herquel G, Tapponnier P, Wittlinger G, et al. 1999. Teleseismic shear wave splitting and lithospheric anisotropy beneath and across the Altyn Tagh fault. Geophysical Research Letters, 26(21): 3225-3228. DOI:10.1029/1999GL005387
Hu X Y, Huo G P, Gao R, et al. 2013. The magnetotelluric anisotropic two-dimensional simulation and case analysis. Chinese Journal of Geophysics (in Chinese), 56(12): 4268-4277. DOI:10.6038/cjg20131229
Huo G P, Hu X Y, Fang H, et al. 2012. Magnetotelluric joint inversion for anisotropic conductivities in layered media. Acta Physica Sinica (in Chinese), 61(12): 129101. DOI:10.7498/aps.61.129101
Huo G P, Hu X Y, Huang Y F, et al. 2015. MT modeling for two-dimensional anisotropic conductivity structure with topography and examples of comparative analyses. Chinese Journal of Geophysics (in Chinese), 58(12): 4696-4708. DOI:10.6038/cjg20151230
Ichihara H, Mogi T. 2009. A realistic 3-D resistivity model explaining anomalous large magnetotelluric phases: the L-shaped conductor model. Geophysical Journal International, 179(1): 14-17. DOI:10.1111/j.1365-246X.2009.04310.x
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. Tectonophysics, 603: 114-122. DOI:10.1016/j.tecto.2013.05.020
Jolivet M, Brunel M, Seward D, et al. 2003. Neogene extension and volcanism in the Kunlun Fault Zone, northern Tibet: New constraints on the age of the Kunlun Fault. Tectonics, 22(5): 1052. DOI:10.1029/2002TC001428
Jones A G, Evans R L, Muller M R, et al. 2009. Area selection for diamonds using magnetotellurics: Examples from southern Africa. Lithos, 112: 83-92. DOI:10.1016/j.lithos.2009.06.011
Karplus M S, Zhao W, Klemperer S L, et al. 2011. Injection of Tibetan crust beneath the south Qaidam Basin: Evidence from INDEPTH IV wide-angle seismic data. Journal of Geophysical Research: Solid Earth, 116(B7): B07301. DOI:10.1029/2010JB00791
Kelbert A, Egbert G D, Schultz A. 2008. Non-linear conjugate gradient inversion for global EM induction: resolution studies. Geophysical Journal International, 173(2): 365-381. DOI:10.1111/j.1365-246X.2008.03717.x
Kelbert A, Meqbel N, Egbert G D, et al. 2014. ModEM: A modular system for inversion of electromagnetic geophysical data. Computers & Geosciences, 66: 40-53. DOI:10.1016/j.cageo.2014.01.010
Kong W X, Lin C H, Tan H D, et al. 2018. The effects of 3D electrical anisotropy on magnetotelluric responses: synthetic case studies. Journal of Environmental and Engineering Geophysics, 23(1): 61-75. DOI:10.2113/JEEG23.1.61
Kumar P G, Manglik A. 2012. Electrical anisotropy in the Main Central Thrust Zone of the Sikkim Himalaya: Inference from anomalous MT phase. Journal of Asian Earth Sciences, 57: 120-127. DOI:10.1016/j.jseaes.2012.06.017
Kurtz R D, Craven J A, Niblett E R, et al. 1993. The conductivity of the crust and mantle beneath the Kapuskasing Uplift: electrical anisotropy in the upper mantle. Geophysical Journal International, 113(2): 483-498. DOI:10.1111/j.1365-246X.1993.tb00901.x
Kuzmin A, Luisier M, Schenk O. 2013. Fast methods for computing selected elements of the Greens function in massively parallel nanoelectronic device simulations. //Wolf F, Mohr B, Ney D eds. Euro-Par 2013 Parallel Processing. Berlin Heidelberg: Springer, 533-544.
Le Pape F, Jones A G, Vozar J, et al. 2012. Penetration of crustal melt beyond the Kunlun fault into northern Tibet. Nature Geoscience, 5(5): 330-335. DOI:10.1038/ngeo1449
Le Pape F, Jones A G, Unsworth M J, et al. 2015. Constraints on the evolution of crustal flow beneath Northern Tibet. Geochemistry Geophysics Geosystems, 16(12): 4237-4260. DOI:10.1002/2015GC005828
Leibecker J, Gatzemeier A, Hönig M, et al. 2002. Evidence of electrical anisotropic structures in the lower crust and the upper mantle beneath the Rhenish Shield. Earth and Planetary Science Letters, 202(2): 289-302. DOI:10.1016/S0012-821X(02)00783-5
Lezaeta P, Haak V. 2003. Beyond magnetotelluric decomposition: Induction, current channeling, and magnetotelluric phases over 90°. Journal of Geophysical Research: Solid Earth, 108(B6): 2305. DOI:10.1029/2001jb000990
Li Y, Pek J, Brasse H. 2003. Magnetotelluric inversion for 2D anisotropic conductivity structures. //Jordt A, Stoll J eds. Protokoll uber das 20 Kolloquium Elektromagnetische Tiefenforschung. Königstein, Germany, 250-259.
Li Y G. 2002. A finite-element algorithm for electromagnetic induction in two-dimensional anisotropic conductivity structures. Geophysical Journal International, 148(3): 389-401. DOI:10.1046/j.1365-246x.2002.01570.x
Li Y G, Pek J. 2008. Adaptive finite element modelling of two-dimensional magnetotelluric fields in general anisotropic media. Geophysical Journal International, 175(3): 942-954. DOI:10.1111/j.1365-246X.2008.03955.x
Liddell M, Unsworth M, Pek J. 2016. Magnetotelluric imaging of anisotropic crust near Fort McMurray, Alberta: implications for engineered geothermal system development. Geophysical Journal International, 205(3): 1365-1381. DOI:10.1093/gji/ggw089
Lin C Y, Wang S M, Sun C C, et al. 2004. One-dimensional anisotropic inversion for magnetotelluric data in Tianzhu-Yongdeng region. Northwestern Seismological Journal (in Chinese), 26(1): 72-77.
Liu Y, Xu Z H, Li Y G. 2018. Adaptive finite element modelling of three-dimensional magnetotelluric fields in general anisotropic media. Journal of Applied Geophysics, 151: 113-124. DOI:10.1016/j.jappgeo.2018.01.012
Löwer A, Junge A. 2017. Magnetotelluric Transfer functions: phase tensor and tipper vector above a simple anisotropic three-dimensional conductivity anomaly and implications for 3D isotropic inversion. Pure and Applied Geophysics, 174(5): 2089-2101. DOI:10.1007/s00024-016-1444-3
Mandolesi E, Jones A G. 2012. Magnetotelluric inversion in a 2D anisotropic environment. //EGU General Assembly Conference Abstracts. Vienna, Austria: EGU.
Martí A. 2014. The role of electrical anisotropy in magnetotelluric responses: from modelling and dimensionality analysis to inversion and interpretation. Surveys in Geophysics, 35(1): 179-218. DOI:10.1007/s10712-013-9233-3
Martinelli P, Osella A. 1997. MT forward modeling of 3-D anisotropic electrical conductivity structures using the Rayleigh-Fourier method. Journal of Geomagnetism and Geoelectricity, 49(11-12): 1499-1518.
McNamara D E, Owens T J, Silver P G, et al. 1994. Shear wave anisotropy beneath the Tibetan Plateau. Journal of Geophysical Research, 99(B7): 13655-13665. DOI:10.1029/93JB03406
Miensopusτ M P, Jones A G. 2011. Artefacts of isotropic inversion applied to magnetotelluric data from an anisotropic Earth. Geophysical Journal International, 187(2): 677-689. DOI:10.1111/j.1365-246x.2011.05157.x
Mo X X, Luo Z H, Deng J F, et al. 2007. Granitoids and crustal growth in the East-Kunlun orogenic belt. Geological Journal of China Universities (in Chinese), 13(3): 403-414.
Montahaei M, Oskooi B. 2014. Magnetotelluric inversion for azimuthally anisotropic resistivities employing artificial neural networks. Acta Geophysica, 62(1): 12-43. DOI:10.2478/s11600-013-0164-7
Newman G A, Hoversten G M. 2000. Solution strategies for two-and three-dimensional electromagnetic inverse problems. Inverse Problems, 16(5): 1357-1375. DOI:10.1088/0266-5611/16/5/314
Oldenburg D W, Haber E, Shekhtman R. 2013. Three-dimensional inversion of multisource time domain electromagnetic data. Geophysics, 78(1): E47-E57. DOI:10.1190/GEO2012-0131.1
Padilha A L, Vitorello Í, Pádua M B, et al. 2006. Lithospheric and sublithospheric anisotropy beneath central-southeastern Brazil constrained by long period magnetotelluric data. Physics of the Earth and Planetary Interiors, 158(2-4): 190-209. DOI:10.1016/j.pepi.2006.05.006
Pain C C, Herwanger J V, Saunders J H, et al. 2003. Anisotropic resistivity inversion. Inverse Problems, 19(5): 1081-1111. DOI:10.1088/0266-5611/19/5/306
Pek J, Verner T. 1997. Finite-difference modelling of magnetotelluric fields in two-dimensional anisotropic media. Geophysical Journal International, 128(3): 505-521. DOI:10.1111/j.1365-246X.1997.tb05314.x
Pek J, Santos F A M. 2002. Magnetotelluric impedances and parametric sensitivities for 1-D anisotropic layered media. Computers & Geosciences, 28(8): 939-950. DOI:10.1016/S0098-3004(02)00014-6
Pek J, Santos F A M. 2006. Magnetotelluric inversion for anisotropic conductivities in layered media. Physics of the Earth and Planetary Interiors, 158(2-4): 139-158. DOI:10.1016/j.pepi.2006.03.023
Pek J, Santos F A M, Li Y G. 2011. Non-linear conjugate gradient magnetotelluric inversion for 2-D anisotropic conductivities. //Smith Y ed. Schmicker-Weidelt-Kolloquium, Nustadt and der Weinstraße, 187-206.
Plotkin V V. 2012. Inversion of heterogeneous anisotropic magnetotelluric responses. Russian Geology and Geophysics, 53(8): 829-836. DOI:10.1016/j.rgg.2012.06.010
Qin L J, Yang C F. 2012. A magnetotelluric inversion method of the whole tensor impedance responses to one-dimensional anisotropic structure. Chinese Journal of Geophysics (in Chinese), 55(2): 693-701. DOI:10.6038/j.issn.0001-5733.2012.02.033
Qin L J, Yang C F, Chen K. 2013. Quasi-analytic solution of 2-D magnetotelluric fields on an axially anisotropic infinite fault. Geophysical Journal International, 192(1): 67-74. DOI:10.1093/gji/ggs018
Reddy I K, Rankin D. 1971. Magnetotelluric effect of dipping anisotropies. Geophysical Prospecting, 19(1): 84-97. DOI:10.1111/j.1365-2478.1971.tb00586.x
Reddy I K, Rankin D. 1975. Magnetotelluric response of laterally inhomogeneous and anisotropic media. Geophysics, 40(6): 1035-1045. DOI:10.1190/1.1440579
Rodi W, Mackie R L. 2001. Nonlinear conjugate gradients algorithm for 2-D magnetotelluric inversion. Geophysics, 66(1): 174-187. DOI:10.1190/1.1444893
Schenk O, Gärtner K. 2004. Solving unsymmetric sparse systems of linear equations with PARDISO. Journal of Future Generation Computer Systems, 20(3): 475-487. DOI:10.1016/j.future.2003.07.011
Shao G H, Xiao Q B. 2016. Model for phase rolling out of quadrant magnetotelluric data-an example from north to the Hexi corridor. Progress in Geophysics (in Chinese), 31(4): 1480-1491. DOI:10.6038/pg20160410
Shapiro N M, Ritzwoller M H, Molnar P, et al. 2004. Thinning and flow of Tibetan crust constrained by seismic anisotropy. Science, 305(5681): 233-236. DOI:10.1126/science.1098276
Shi D N, Shen Y, Zhao W J, et al. 2009. Seismic evidence for a Moho offset and south directed thrust at the easternmost Qaidam-Kunlun boundary in the Northeast Tibetan plateau. Earth and Planetary Science Letters, 288(1-2): 329-334. DOI:10.1016/j.epsl.2009.09.036
Silver P G. 1996. Seismic anisotropy beneath the continents: Probing the depths of geology. Annual Review of Earth and Planetary Sciences, 24: 385-432. DOI:10.1146/annurev.earth.24.1.385
Simpson F. 2001. Resistance to mantle flow inferred from the electromagnetic strike of the Australian upper mantle. Nature, 412(6847): 632-635. DOI:10.1038/35088051
Siripunvaraporn W, Egbert G. 2000. An efficient data-subspace inversion method for 2-D magnetotelluric data. Geophysics, 65(3): 791-803. DOI:10.1190/1.1444778
Siripunvaraporn W, Egbert G, Lenbury Y, et al. 2005a. Three-dimensional magnetotelluric inversion: data-space method. Physics of the Earth and Planetary Interiors, 150(1-3): 3-14. DOI:10.1016/j.pepi.2004.08.023
Siripunvaraporn W, Egbert G, Uyeshima M. 2005b. Interpretation of two-dimensional magnetotelluric profile data with three-dimensional inversion: synthetic examples. Geophysical Journal International, 160(3): 804-814. DOI:10.1111/j.1365-246X.2005.02527.x
Soto G L, Sandvol E, Ni J F, et al. 2012. Significant and vertically coherent seismic anisotropy beneath eastern Tibet. Journal of Geophysical Research: Solid Earth, 117(B5): B05308. DOI:10.1029/2011JB008919
Tapponnier P, Xu Z Q, Roger F, et al. 2001. Oblique Stepwise Rise and Growth of the Tibet Plateau. Science, 294(5547): 1671-1677. DOI:10.1126/science.105978
Tikhonov A N, Arsenin V Y. 1977. Solutions of Ill-Posed Problems. Winston, Washington: John Wiley & Sons, Inc.
Unsworth M, Wei W, Jones A G, et al. 2004. Crustal and upper mantle structure of northern Tibet imaged with magnetotelluric data. Journal of Geophysical Research: Solid Earth, 109(B2): B02403. DOI:10.1029/2002JB002305
Wang D J, Mookherjee M, Xu Y S, et al. 2006. The effect of water on the electrical conductivity of olivine. Nature, 443(7114): 977-980. DOI:10.1038/nature05256
Wang Q, McDermott F, Xu J F, et al. 2005. Cenozoic K-rich adakitic volcanic rocks in the Hohxil area, northern Tibet: Lower-crustal melting in an intracontinental setting. Geology, 33(6): 465-468. DOI:10.1130/G21522.1
Wannamaker P E. 2005. Anisotropy versus heterogeneity in continental solid earth electromagnetic studies: fundamental response characteristics and implications for physicochemical state. Surveys in Geophysics, 26(6): 733-756. DOI:10.1007/s10712-005-1832-1
Weckmann U, Ritter O, Haak V. 2003. A magnetotelluric study of the Damara Belt in Namibia. Physics of the Earth and Planetary Interiors, 138(2): 91-112. DOI:10.1016/S0031-9201(03)00079-7
Wei W B, Le Pape F, Jones A G, et al. 2014. Northward channel flow in northern Tibet revealed from 3D magnetotelluric modelling. Physics of the Earth and Planetary Interiors, 235: 13-24. DOI:10.1016/j.pepi.2014.07.004
Weidelt P. 1999. 3D conductivity models: implications of electrical anisotropy. //Oristaglio M, Spies B eds. Three-Dimensional Electromagnetics. Salt Lake City: Society of Exploration Geophysicists Publishers, 119-137.
Wu C L, Xu T, Badal J, et al. 2015. Seismic anisotropy across the Kunlun fault and their implications for northward transforming lithospheric deformation in northeastern Tibet. Tectonophysics, 659: 91-101. DOI:10.1016/j.tecto.2015.07.030
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. Tectonophysics, 601: 125-138. DOI:10.1016/j.tecto.2013.05.003
Xiao Q B, Shao G H, Yu G, et al. 2016. Electrical resistivity structures of the Kunlun-Qaidam-Qilian system at the northern Tibet and their tectonic implications. Physics of the Earth and Planetary Interiors, 255: 1-17. DOI:10.1016/j.pepi.2016.03.011
Xiao Q B, Yu G, Liu Z J, et al. 2017. Structure and geometry of the Aksay restraining double bend along the Altyn Tagh Fault, northern Tibet, imaged using magnetotelluric method. Geophysical Research Letters, 44(9): 4090-4097. DOI:10.1002/2017GL072581
Xiao Q B, Yu G, Shao G H, et al. 2018. Lateral rheology differences in the lithosphere and dynamics as revealed by magnetotelluric imaging at the northern Tibetan Plateau. Journal of Geophysical Research: Solid Earth, 123(9): 7266-7284. DOI:10.1029/2017JB015285
Xiao T J, Wang Y, Huang X Y, et al. 2020. MT responses of three-dimensional conductive and magnetic anisotropic anomalies. Geophysical Prospecting, 68(3): 1016-1040. DOI:10.1111/1365-2478.12886
Xiong B, Luo T Y, Li C W, et al. 2013. Modelling of magnetotelluric responses of 2-D electrical anomalous body in anisotropic media using finite element method. Acta Scientiarum Naturalium Universitatis Sunyatseni (in Chinese), 52(4): 143-148.
Xu S Z, Zhao S K. 1985. Solution of magnetotelluric field equations for a two-dimensional, anisotropic geoelectric section by the finite element method. Acta Seismologica Sinica (in Chinese), 7(1): 80-90.
Yang C F. 1997. MT numerical simulation of symmetrically 2-D anisotropic media based on the finite element method. Northwestern Seismological Journal (in Chinese), 19(2): 27-33.
Yang C F, Lin C Y, Sun C C, et al. 2005. Inversions for MT data in 2-D symmetrical anisotropic media. Acta Seismologica Sinica (in Chinese), 27(3): 339-345. DOI:10.3321/j.issn:0253-3782.2005.03.013
Yang M X, Tan H D, Wang M, et al. 2015. Two-dimensional magnetotelluric forward research for the vertical anisotropy. International Journal of Geosciences, 6(10): 1166-1172. DOI:10.4236/ijg.2015.610091
Yee K. 1966. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Transactions on Antennas and Propagation, 14(3): 302-307. DOI:10.1109/TAP.1966.1138693
Yin C C. 2003. Inherent nonuniqueness in magnetotelluric inversion for 1D anisotropic models. Geophysics, 68(1): 138-146. DOI:10.1190/1.1543201
Yoshino T, Matsuzaki T, Yamashita S, et al. 2006. Hydrous olivine unable to account for conductivity anomaly at the top of the asthenosphere. Nature, 443(7114): 973-976. DOI:10.1038/nature05223
Yu G, Xiao Q B, Zhao G Z, et al. 2018. Three-dimensional magnetotelluric responses for arbitrary electrically anisotropic media and a practical application. Geophysical Prospecting, 66(9): 1764-1783. DOI:10.1111/1365-2478.12690
Yu G, Xiao Q B, Li M. 2019. Anisotropic model study for the phase roll out of quadrant data in magnetotellurics: with examples of upper-lower structure. Chinese Journal of Geophysics (in Chinese), 62(2): 763-778. DOI:10.6038/cjg2019L0661
Yu G, Li M, Shao G H, et al. 2020. Electrical resistivity variations of lithospheric mantle beneath the northern Tibetan Plateau with tectonic implications. Journal of Asian Earth Sciences, 198: 104273. DOI:10.1016/j.jseaes.2020.104273
Zhu L P, Helmberger D V. 1998. Moho offset across the northern margin of the Tibetan Plateau. Science, 281(5380): 1170-1172. DOI:10.1126/science.281.5380.1170
邓起东, 张培震, 冉勇康, 等. 2003. 中国活动构造与地震活动. 地学前缘, 10(S1): 66-73. DOI:10.3321/j.issn:1005-2321.2003.z1.012
胡祥云, 霍光谱, 高锐, 等. 2013. 大地电磁各向异性二维模拟及实例分析. 地球物理学报, 56(12): 4268-4277. DOI:10.6038/cjg20131229
霍光谱, 胡祥云, 方慧, 等. 2012. 层状各向异性介质大地电磁联合反演研究. 物理学报, 61(12): 129101. DOI:10.7498/aps.61.129101
霍光谱, 胡祥云, 黄一凡, 等. 2015. 带地形的大地电磁各向异性二维模拟及实例对比分析. 地球物理学报, 58(12): 4696-4708. DOI:10.6038/cjg20151230
林长佑, 王书明, 孙崇赤, 等. 2004. 天祝-永登地区大地电磁资料的一维各向异性反演. 西北地震学报, 26(1): 72-77.
莫宣学, 罗照华, 邓晋福, 等. 2007. 东昆仑造山带花岗岩及地壳生长. 高校地质学报, 13(3): 403-414. DOI:10.3969/j.issn.1006-7493.2007.03.010
秦林江, 杨长福. 2012. 大地电磁全张量响应的一维各向异性反演. 地球物理学报, 55(2): 693-701. DOI:10.6038/j.issn.0001-5733.2012.02.033
邵贵航, 肖骑彬. 2016. 相位超象限大地电磁观测数据的模型研究——以河西走廊北侧为例. 地球物理学进展, 31(4): 1480-1491. DOI:10.6038/pg20160410
熊彬, 罗天涯, 李长伟, 等. 2013. 用有限单元法模拟各向异性介质中二维电性异常体的大地电磁响应. 中山大学学报(自然科学版), 52(4): 143-148.
徐世浙, 赵生凯. 1985. 二维各向异性地电断面大地电磁场的有限元法解法. 地震学报, 7(1): 80-90.
杨长福. 1997. 大地电磁二维对称各向异性介质的有限元数值模拟. 西北地震学报, 19(2): 27-33.
杨长福, 林长佑, 孙崇赤, 等. 2005. 二维对称各向异性介质大地电磁反演. 地震学报, 27(3): 339-345. DOI:10.3321/j.issn:0253-3782.2005.03.013
喻国, 肖骑彬, 李满. 2019. 大地电磁相位超象限现象的各向异性模型研究——以上下结构为例. 地球物理学报, 62(2): 763-778. DOI:10.6038/cjg2019L0661