大地電磁二維各向異性反演及其在青藏高原北部的應用

喻國， 肖騎彬， 李滿

1 中國地震局地質(zhì)研究所， 北京 100029 2 中國地震局地質(zhì)研究所地震動力學國家重點實驗室, 北京 100029

0 引言

地球介質(zhì)的電性特征受裂縫以及礦物的定向排列、含水量、部分熔融等的影響(Wannamaker，2005；Wang et al.，2006；Yoshino et al.，2006；Dai and Karato，2014)，在不同層次上(地殼以及上地幔)都可能表現(xiàn)出各向異性.而分辨以及解釋包含各向異性效應的數(shù)據(jù)也一直以來是大地電磁測深(MT)方法領域的難點，在針對大地電磁實際數(shù)據(jù)的各向異性研究中，大多數(shù)是應用維性工具、正演模擬等定性分析手段去解釋各向異性結(jié)構(gòu)(Wannamaker，2005；Martí，2014).在較早的研究中，人們觀察到大地電磁阻抗張量非對角元素的相位曲線有時在長周期存在明顯的差異，認為這種相位分離現(xiàn)象是地球深部電各向異性的表現(xiàn)，并有研究將MT數(shù)據(jù)中的相位分離現(xiàn)象與地震學中的橫波分離現(xiàn)象進行對比分析(Kurtz et al.，1993；Simpson，2001；Bahr and Simpson，2002；Leibecker et al.，2002；Eaton et al.，2004; Eaton and Jones，2006；Padilha et al.，2006).然而，隨著模型研究的深入，Heise等(2006)指出相位分離現(xiàn)象是由電導率的空間差異或梯度分布引起的，這種現(xiàn)象并不能指征電導率是各向同性的或各向異性的本質(zhì)特征，因而，MT數(shù)據(jù)的相位分離與地震橫波分離是完全不同類別的現(xiàn)象.相位分離的一種特殊現(xiàn)象是低頻段的相位超象限現(xiàn)象(Chouteau and Tournerie，2000)，它可以由各向異性感應產(chǎn)生，并被作為各向異性結(jié)構(gòu)可能存在的印跡(Heise and Pous，2003；Kumar and Manglik，2012；Martí，2014；Liddell et al.，2016).但是，相位超象限現(xiàn)象同樣也可以通過三維各向同性結(jié)構(gòu)來模擬(Lezaeta and Haak，2003；Weckmann et al.，2003； Ichihara and Mogi，2009； Ichihara et al.，2013；Xiao et al.，2016；邵貴航和肖騎彬，2016).另外，理論模型研究表明，含有本征各向異性體的模型響應可以用較為復雜的各向同性結(jié)構(gòu)(如宏觀各向異性結(jié)構(gòu))進行模擬(Eisel and Haak，1999；Heise and Pous，2001；Martí，2014)，這種等效性主要是因為在相應的探測深度上MT方法分辨率不足而導致的(Weidelt，1999)，這也提供了利用特殊的等效各向同性結(jié)構(gòu)解釋各向異性的可能性.但是，這種等效性的存在也會使得我們通過各向同性反演獲得的電阻率結(jié)構(gòu)發(fā)生明顯的扭曲(Miensopust and Jones，2011；L?wer and Junge，2017)，從而對研究區(qū)的結(jié)構(gòu)產(chǎn)生誤導.此外，各向異性結(jié)構(gòu)并不是都存在等效的各向同性結(jié)構(gòu)，例如，簡單的均勻各向異性半空間的響應無法通過構(gòu)建等效的各向同性結(jié)構(gòu)來獲得(Martí，2014).這些研究結(jié)果表明，探測各向異性是極其復雜的過程，僅僅依靠資料的定性分析以及正演模擬等手段來判斷地下深部的電各向異性結(jié)構(gòu)還存在著較大的爭議，而為了尋找更可靠的各向異性模型來合理解釋地下可能存在的結(jié)構(gòu)，我們有必要完成相應的反演算法，同時積累和總結(jié)實際數(shù)據(jù)應用中的經(jīng)驗和認識.

大地電磁各向異性正、反演方法的研究最早可追溯到20世紀60年代，接下來直到70年代末期關于各向異性的大地電磁一維正、反演和二維正演也一直是研究熱點.Reddy和Rankin(1971)討論了層狀模型中，傾斜各向異性對大地電磁響應的影響；Reddy和Rankin(1975)利用有限元算法最早實現(xiàn)了二維各向異性模型響應的計算；Abramovici(1974)完成了一維任意各向異性的正演算法；Abramovici和Shoham(1977)應用矩陣廣義逆的方法完成針對一維垂直各向異性模型的反演.而在80年代至90年代中期關于大地電磁各向異性問題的研究成果發(fā)表量明顯減少(徐世浙和趙生凱，1985).90年代后期至今，Pek和Vener(1997)利用有限差分法實現(xiàn)了計算二維任意各向異性模型的正演算法；Yin(2003)討論了大地電磁各向異性反演固有的多解性；Pek和Santos(2006)將標準的一維奧卡姆反演拓展到各向異性情況；Pek等(2011)完成二維任意各向異性反演算法，并研究了不同參數(shù)對反演的影響.隨著計算機技術的快速革新，關于各向異性的大地電磁正、反演方法再次成為研究熱點(楊長福，1997；Martinelli and Osella，1997；Weidelt，1999；Pek and Santos，2002； Li，2002；Li et al.，2003；Li and Pek，2008；楊長福等，2005；Baba et al.，2006；Mandolesi and Jones，2012；Chen and Weckmann，2012；Plotkin，2012；秦林江和楊長福，2012；胡祥云等，2013；霍光譜等，2012，2015；熊彬等，2013；Qin et al.，2013；Yang et al.，2015；Montahaei and Oskooi，2014；L?wer and Junge，2017；Cao et al.，2018；Kong et al.，2018；Liu et al.，2018；Yu et al.，2018；Xiao et al.，2020；喻國等，2019).

然而，雖然基于各向異性介質(zhì)的大地電磁正演技術已經(jīng)可以進行復雜的模型計算，但是各向異性反演在實際數(shù)據(jù)中的應用依然進展緩慢，這其中一個關鍵因素在于對各向異性參數(shù)的恢復上存在很大的局限(Yin，2003；Pek et al.，2011；Chen and Weckmann，2012).Baba等(2006)首先將二維各向異性反演應用在穿過東太平洋南部洋中脊的大地電磁剖面的解釋中，算法在Rodi和Mackie(2001)的二維各向同性反演代碼的基礎上修改完成了針對三個主軸電阻率的垂直各向異性反演.Baba等(2006)在對反演結(jié)果的分析中指出，相比于各向同性的反演結(jié)果，各向異性反演并沒有明顯的改善對數(shù)據(jù)的擬合程度，因此無法單獨依靠大地電磁數(shù)據(jù)分析來確定各向異性模型為最優(yōu)模型；然而，通過其他地球物理觀測資料與大地電磁的各向異性模型之間的相互印證，才表明得到的各向異性模型是可靠的結(jié)果.Pek和Santos(2006)將一維各向異性反演得到的模型作為二維結(jié)構(gòu)試錯過程中的初始模型.Pek等(2011)也將他們完成的二維各向異性反演代碼應用于葡萄牙南部采集的大地電磁數(shù)據(jù)，他們認為二維各向異性反演結(jié)果有可能受到觀測數(shù)據(jù)中存在的非二維效應的影響，但沒有對數(shù)據(jù)作深入探究；其中，反演待解模型參數(shù)為三個主軸電阻率、各向異性走向角以及傾斜角，但反演最終模型整體上呈現(xiàn)較強的各向同性.林長佑等(2004)將一維大地電磁垂直各向異性反演用在天?！赖堑貐^(qū)實際資料的定性解釋中.楊長福等(2005)開發(fā)了二維大地電磁垂直各向異性反演算法(針對三個主軸電阻率)，認為當電導率沿垂向和傾向表現(xiàn)相同的各向異性時，可以分別用TE和TM極化反演結(jié)果來代表兩個主軸方向電阻率模型.霍光譜等(2012)基于馬夸特方法，完成了大地電磁層狀任意各向異性的視電阻率和相位反演，并指出需要將數(shù)據(jù)揭示的電各向同性與各向異性現(xiàn)象區(qū)分，否則無法得到準確合理的地層信息，甚至會出現(xiàn)自相矛盾的情況.秦林江和楊長福(2012)應用廣義逆矩陣方法完成一維大地電磁各向異性反演算法，指出對實測資料反演時，如果數(shù)據(jù)是受各向異性影響的，那么必須用各向異性的處理方法才能得到更客觀真實的結(jié)果.Cao等(2018)應用有限內(nèi)存擬牛頓方法完成了針對三個主軸電阻率的大地電磁反演算法，但還沒在實際中應用，并且目前還沒有其他的大地電磁三維各向異性反演算法完成.通過以上國內(nèi)外大地電磁各向異性反演應用的情況可以知道，大地電磁各向異性反演算法還不完善，一維各向異性反演在實際應用中存在明顯的局限性；二維各向異性反演的實例較少；三維反演目前還沒有對數(shù)據(jù)應用的實例.因此，合理全面地認識和理解地下存在的電各向異性結(jié)構(gòu)，需要進一步發(fā)展各向異性反演算法和更多地開展在大地電磁實際數(shù)據(jù)中的應用.

現(xiàn)今基于三維各向同性介質(zhì)的反演技術近年來發(fā)展十分迅速并已經(jīng)在實際資料解釋中廣泛應用(Siripunvaraporn et al.，2005a，b；Egbert and Kelbert，2012；Kelbert et al.， 2014).Kelbert等(2014)完成的ModEM程序?qū)⒎囱葸^程封裝成不同功能的單獨模塊，為MT方法的推廣應用創(chuàng)造了條件.本文在已經(jīng)實現(xiàn)的二維任意各向異性大地電磁正演算法的基礎上(Yu et al.，2018)，構(gòu)建基于各向異性介質(zhì)的反演目標函數(shù)，實現(xiàn)對目標函數(shù)值及其梯度的計算.參考ModEM的非線性共軛梯度法(NLCG)反演模塊，實現(xiàn)對各向異性介質(zhì)的二維大地電磁反演，通過理論模型驗證反演結(jié)果的穩(wěn)定性.在此基礎上，我們對東昆侖—柴達木盆地過渡帶一條MT剖面數(shù)據(jù)進行各向異性反演解釋，該剖面所在的測區(qū)的地下地震各向異性背景豐富，但還沒有相關電各向異性結(jié)構(gòu)的研究.我們結(jié)合該地區(qū)以往的地震各向異性研究成果以及本文各向異性反演得到的結(jié)果來判斷深部各向異性結(jié)構(gòu)的存在性，并進一步分析各向異性的成因及其構(gòu)造意義.

1 正演模擬

在大地電磁方法中，電場E和磁場H滿足麥克斯韋方程組，通過基本的數(shù)學代數(shù)運算可以得到電場矢量方程或磁場矢量方程.在二維各向異性介質(zhì)中，因為似穩(wěn)假設的Maxwell方程組無法再解耦成單獨的TE或TM模式，因此這里我們采用一般形式的電場矢量方程

方程組具體表達如下：

(1)

(2)

其中，Rz表示沿Z軸旋轉(zhuǎn)的旋轉(zhuǎn)矩陣，Rx表示沿X軸旋轉(zhuǎn)的旋轉(zhuǎn)矩陣，T表示矩陣轉(zhuǎn)置.

在右旋笛卡爾坐標系下，對模型進行矩形網(wǎng)格剖分，采用經(jīng)典的交錯網(wǎng)格格式來定義場的位置(Yee，1966).對方程(1)進行數(shù)值離散，得到線性方程組

AE=B，

(3)

其中，A為系數(shù)矩陣，B為右端項，E為所有的未知電場.方程(3)采用直接法求解器PARDISO(Schenk and G?rtner，2004；Kuzmin et al.，2013)進行求解.在求得電場值E后，構(gòu)建插值向量egi和hgi(下標i=x，y)，統(tǒng)一用解向量E計算地表指定點的電場egiE和磁場hgiE(Newman and Alumbaugh，2000).

計算相互垂直的兩組初始水平極化條件下的電場E1和E2后，求得大地電磁響應函數(shù)，以在某一觀測點的阻抗張量元素Zxy為例，

(4)

2 反演理論

大地電磁各向異性反演問題，可以歸結(jié)為求解以下目標函數(shù)的極小值問題(Tikhonov and Arsenin，1977)：

φ=φd+λsφs+λaφa.

(5)

相比較于各向同性反演，構(gòu)建的目標函數(shù)中除了包含數(shù)據(jù)擬合項和模型粗糙度項外，還引入了各向異性懲罰項(Pain et al.，2003；Pek and Santos，2006；Pek et al.，2011；Chen and Weckmann，2012).

φa=mTKm.

φa=mTKm=∑[(lnρx-lnρy)2+(lnρx-lnρz)2+(lnρy-lnρz)2]

其一般形式為

φa=mTKm=∑[a12(lnρx-lnρy)2+a13(lnρx-lnρz)2+a23(lnρy-lnρz)2]

其中a12，a13和a23為正數(shù).

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

使得

(6)

(7)

(8)

(9)

該式表示依次繞Z軸旋轉(zhuǎn)α角和γ角，α和γ表示相同的旋轉(zhuǎn)角，所以可以將走向角α和偏向角γ合并為同一個參數(shù)α，如下：

(10)

因此，綜合考慮大地電磁各向異性反演的特點，我們在反演中設置待求解的模型參數(shù)只包括ρx、ρy、ρz和α.

(6)—(8)式中，計算量主要集中在(8)式中數(shù)據(jù)對模型的偏導數(shù)計算上.這里N是Np(頻率點數(shù))、Ns(測點數(shù))和Nr(每個測點的響應函數(shù)數(shù)量)的乘積.取d包含全部阻抗張量元素(Zxx，Zxy，Zyx，Zyy)，那么(8)式可以展開為

(11)

其中，k代表迭代次數(shù)，α為當前搜索步長，p為模型搜索方向.反演計算中實際需要調(diào)整的主要參數(shù)包括：初始正則化參數(shù)λs和λa，正則化參數(shù)更新因子k，初始線搜索步長α0，數(shù)據(jù)誤差門檻，以及與模型協(xié)方差相關的光滑參數(shù).

在反演迭代的計算過程中，當數(shù)據(jù)擬合均方根誤差(RMS)減小量低于某個閥值時，正則化參數(shù)更新為λs/k和λa/k.而初始線搜索步長α0在這里采取人為設定的方式，根據(jù)大量反演算例的經(jīng)驗，初始線搜索步長的大小會顯著影響反演的收斂過程，當步長過大，模型在初始階段易發(fā)生較大調(diào)整而產(chǎn)生冗余構(gòu)造，在反演后期難以校正，并且這個影響在各向異性反演中更為突出.因此，在實際操作中，我們需要盡量避免在首次迭代中出現(xiàn)很大的模型參數(shù)調(diào)整.

3 理論模型的反演檢驗

模型1(Md1)為一包含三個各向異性體的二維模型(圖1a).模型背景電阻率為100 Ωm，在8.65 km以上包含兩個不同規(guī)模的各向異性體A1和A2，它們的頂部埋深都是0.135 km，寬度同為3.2 km.兩個異常體之間相距0.8 km，底部埋深分別為8.65 km和1.36 km.兩個異常體具有相同的各向異性參數(shù)：ρx=1 Ωm、ρy=200 Ωm、ρz=1 Ωm、α=30°、β=0°、γ=0°.異常體A1下伏一個各向異性層，層厚14.42 km，模型參數(shù)為ρx=5 Ωm、ρy=50 Ωm、ρz=5 Ωm、α=120°、β=0、γ=0.將該模型在水平方向剖分為46個網(wǎng)格，地表向下剖分為39個網(wǎng)格，地表定義26個觀測點，計算周期范圍0.015627～15627 s內(nèi)29個周期點的阻抗響應數(shù)據(jù).

如圖1b所示，兩次反演中擬合RMS偏差相對較大的點是位于異常體A1和A2之間的測點13和14，這表明二維電性結(jié)構(gòu)橫向的劇烈變化與各向異性效應一起增加了數(shù)據(jù)擬合的難度.從模型恢復程度來看(圖1c—j)，ρz出現(xiàn)了一些細微的構(gòu)造，但依然無法實現(xiàn)對異常結(jié)構(gòu)的有效分辨，這主要是因為ρz在二維反演中幾乎無法還原，其結(jié)構(gòu)的變化主要是受目標函數(shù)中各向異性約束項的影響(Pek et al.，2011).

圖1c和d表明，在λs=λa=0.001時，反演恢復A1和A2區(qū)域的各向異性參數(shù)特點表現(xiàn)為，主軸電阻率ρx為200 Ωm左右的相對高阻ρmax，ρy為1 Ωm左右的相對低阻ρmin，這與真實模型中水平主軸電阻率的分布特點(ρx=1 Ωm表現(xiàn)為低阻ρmin，ρy=200 Ωm表現(xiàn)為高阻ρmax)相反，而在圖1f中，A1和A2區(qū)域恢復的走向角α約為-60°，與真實的α=30°相差90°.從單獨的各向異性參數(shù)來看，還原的參數(shù)與真實參數(shù)不同，但是反演還原的參數(shù)組(ρx=ρmax，ρy=ρmin，α=-60°)與真實的參數(shù)組(ρx=ρmin，ρy=ρmax，α=30°)代表的是相同的各向異性體，如(12)式所示，兩組參數(shù)表述的是相同的電導率張量，因此，初始λs=λa=0.001計算的反演結(jié)果圖1c、d和f基本還原了A1和A2區(qū)域真實的各向異性異常體，還原的異常體參數(shù)與真實異常參數(shù)等效.同樣，對圖1g、h和j也可以得出相同的結(jié)論.

(12)

對于上部異常體A1和A2的尺度恢復，異常中心越集中在淺部越有利于準確地把握異常體的尺度，這也符合大地電磁方法淺部分辨率高的特點(Oldenburg et al.，2013)，走向方位角也能較好恢復異常體的尺度.對于最下方層狀異常體，兩次反演結(jié)果都不能實現(xiàn)對下部異常體的恢復，但是在反演的電性結(jié)構(gòu)中，模型淺部兩側(cè)出現(xiàn)了相對異常結(jié)構(gòu)(如圖1c中的C1和C2、圖1d中的B1和B2).結(jié)合測點RMS擬合偏差值觀察，兩側(cè)的測點擬合程度高，因此，兩側(cè)的冗余結(jié)構(gòu)是反演的多解性所導致的.

圖1 二維各向異性理論模型及其反演結(jié)果(a)二維理論模型Md1示意圖及其視電阻率和相位響應擬斷面圖；圖中小倒三角表示地表觀測點位置，A1、A2、A3表示各向異性體.(b)取不同正則化參數(shù)反演得到的各測點RMS值.紅色實點表示λs=0.001和λa=0.001時的擬合RMS值；藍色空心點表示λs=0.1和λa=0.1時的擬合RMS值.(c)—(f)當取正則化參數(shù)λs=0.001和λa=0.001時反演得到的模型與理論模型對比，黑框為理論模型中異常體的范圍；(c)—(f)依次為模型參數(shù)ρx、ρy、ρz和α的剖面圖.(g)—(j)當取正則化參數(shù)λ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.

圖2 不同正則化參數(shù)組合下對應數(shù)據(jù)擬合RMS值的變化Fig.2 The RMS values calculated by different regularization parameter combinations

4 東昆侖—柴達木過渡帶巖石圈電各向異性研究

4.1 研究區(qū)域各向異性背景

東昆侖—柴達木過渡帶東西跨度約1000 km，橫亙于青藏高原的北部(圖3a).它們之間以柴達木南緣斷裂(SQDF)相隔(圖3b)，東昆侖造山帶大陸地殼主要形成于古元古代晚期，自顯生宙以來，東昆侖經(jīng)歷了多期洋陸演化，同時，東昆侖造山帶也是一條巨型構(gòu)造巖漿巖帶(莫宣學等，2007).左行走滑的東昆侖斷裂帶與東昆侖山相伴而生，它繼承了早期縫合帶特征，并在新生代重新活動.該斷裂帶現(xiàn)今的走滑速率在中段約11 mm·a-1，并在印度—亞洲板塊匯聚過程中吸收北東向擠壓而產(chǎn)生縮短變形，對青藏高原的形成及演化發(fā)揮了重要作用(Tapponnier et al.，2001).

在青藏高原中部，觀察到的地震各向異性可以通過中下地殼的弱流動層解釋(Shapiro et al.，2004)，同時該軟弱流動層在電性上也表現(xiàn)出較強的各向異性(Le Pape et al.，2015)，但關于它是否越過東昆侖進入到柴達木盆地之下，這仍存在爭議(Zhu and Helmberger，1998；Unsworth et al.，2004；Shi et al.，2009；Karplus et al.，2011；Le Pape et al.，2012；Wei et al.，2014).如圖3b，我們先后跨過東昆侖造山帶與柴達木盆地之間的過渡帶的不同位置采集了5條大地電磁測深剖面的數(shù)據(jù)，并在前期的研究工作中發(fā)現(xiàn)，在剖面L15中祁漫塔格下方的上地幔頂部和下地殼底部存在低阻異常體(圖3c中C1和C2異常體)(Xiao et al.，2013，2016，2018).同時，該上地幔頂部低阻異常體對應的地面向南約150 km處的地表存在新生代的火山巖(Jolivet et al.，2003；Wang et al.，2005)，Wang等(2005)認為這是增厚的松潘—甘孜地殼下部發(fā)生部分熔融的證據(jù).

在這一地區(qū)的地震快波偏振方向顯示，從阿爾金斷裂帶到柴達木盆地，快波偏振方向從與阿爾金斷裂帶走向平行向近東西向轉(zhuǎn)變，發(fā)生了順時針旋轉(zhuǎn)，而少部分的快波偏振方向與地表山脈結(jié)構(gòu)走向不一致，其中，在祁漫塔格內(nèi)部的一個測點顯示最大左行剪切方向為北東向，與祁漫塔格山脈走向不同(Herquel et al.，1999； Soto et al.，2012)，這也為青藏高原北部的動力學過程提供了新的認識.Le Pape等(2012)在重新獲取大地電磁剖面的各向異性模型時，發(fā)現(xiàn)各向異性的存在使東昆侖與柴達木接觸帶的地殼結(jié)構(gòu)發(fā)生了明顯的變化，為松潘—甘孜中下地殼低阻流動層繼續(xù)向北穿透提供了新的證據(jù).這一節(jié)中，我們利用新的大地電磁各向異性反演程序?qū)15剖面數(shù)據(jù)進行重新反演，進一步分析在祁漫塔格與柴達木過渡部位低阻通道的各向異性特征，并在下一節(jié)的討論中利用電各向異性參數(shù)的變化樣式來分析青藏高原北部軟物質(zhì)流的存在及區(qū)域應力環(huán)境.

4.2 二維各向異性結(jié)構(gòu)

圖3 東昆侖—柴達木盆地過渡帶和鄰區(qū)構(gòu)造背景以及前期大地電磁測點分布圖(a)研究區(qū)大地構(gòu)造背景，圖中黃色粗箭頭表示區(qū)域應力方向；ATF：阿爾金斷裂帶；KLF：東昆侖斷裂帶；HYF：海原斷裂帶；RRF：紅河斷裂帶；XXF：鮮水河—小江斷裂帶；JLF：嘉黎斷裂帶.(b)東昆侖—柴達木盆地地質(zhì)圖，圖中三角形表示大地電磁測深點；其中紅色為本次研究的剖面；SQDF：柴達木南緣斷裂.(c)沿L15剖面早期二維各向同性反演結(jié)果，剖面方向從SW至NE(Xiao et al.，2018).斷裂構(gòu)造參考鄧起東等(2003). 地質(zhì)圖據(jù)1:250萬中國地質(zhì)圖，在線瀏覽地址: 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/.

反演的模型網(wǎng)格為84×240，84為橫向網(wǎng)格剖分，240是地下分層數(shù)，模型加入地形的變化，初始模型電阻率為100 Ωm.初始線搜索步長α0=20.

首先，我們進行三組數(shù)據(jù)分類的各向同性反演：對Zxy和Zyx數(shù)據(jù)的反演、對單獨Zxy分量的反演以及對單獨Zyx分量的反演.通過各向異性反演程序?qū)崿F(xiàn)各向同性反演，只需設定模型參數(shù)ρx=ρy=ρz，α=β=γ=0°，這樣處理后，目標函數(shù)中各向異性懲罰項φa始終為0，不再對反演過程產(chǎn)生作用.然后，選取全部阻抗張量數(shù)據(jù)(Zxx、Zxy、Zyx、Zyy)，進行反演，待還原的參數(shù)仍為ρx、ρy、ρz和α.

從圖4a對Zxy和Zyx數(shù)據(jù)的各向同性反演結(jié)果來看，模型的結(jié)構(gòu)特征與Xiao等(2018)中各向同性反演的二維電性結(jié)構(gòu)相似，在剖面最南端的上地幔頂部出現(xiàn)了低阻異常C1，并向北延伸，與柴達木盆地下地殼層次的低阻體C2相連，形成低阻帶iso1(圖4a)，與圖3c中的異常體C1和C2對應(Xiao et al.，2018).然而，這種結(jié)構(gòu)在單獨的Zxy分量數(shù)據(jù)以及單獨Zyx分量數(shù)據(jù)的反演結(jié)果中并不明顯(圖4b和c).

其次，我們在進行阻抗全張量的各向異性反演計算并測試了大量不同正則因子組合后，選擇粗糙度正則因子為λs=0.1，其得到的各向異性模型與各向同性模型的結(jié)構(gòu)特征更接近.另外，我們在圖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的結(jié)果模型對比可以看出，λa在0.1到10.0范圍內(nèi)得到的各向異性結(jié)構(gòu)比較相似，隨著λa增大，模型中電各向異性程度降低，當λa達到100.0時，各向異性結(jié)構(gòu)不再明顯，在λa=1000.0時，ρx、ρy以及ρz的電阻率結(jié)構(gòu)幾乎一致，并且沒有明顯的各向異性走向角α的結(jié)構(gòu)變化，模型接近于各向同性結(jié)構(gòu)；從圖5可以看出，當λa從小到大變化時，最終RMS值在λa=1.0時達到最小，然后呈上升趨勢，其中在λa=10.0時，數(shù)據(jù)擬合最終RMS開始顯著增大，對應的目標函數(shù)值也明顯變大，各向異性懲罰項的目標函數(shù)值明顯變小.

根據(jù)對圖4和圖5的分析，我們選擇RMS值最小即數(shù)據(jù)擬合最好的結(jié)果作為各向異性反演的最終模型(圖4h—k)，正則因子組合為λs=0.1和λa=1.0.圖6a—d也同為該各向異性模型，ρx和ρy分量的電性結(jié)構(gòu)具有較大的差異，表現(xiàn)出明顯的各向異性特征，其中，ρx剖面中的低阻體C1所在位置對應ρy剖面中的相對高阻體R1，C1核心區(qū)域各向異性走向角α約為45°；ρx剖面中相對高阻體R2所在位置對應ρy剖面中的低阻體C2，C2核心區(qū)域各向異性走向方位角α大約為-45°.這里x軸垂直于測線并指向西北，根據(jù)兩個異常體的各向異性走向角的旋轉(zhuǎn)效應，兩者的低阻電阻率分量C1和C2所在軸的空間方位是在與x軸指向方向夾角45°(順時針旋轉(zhuǎn)45°)的軸線上，即兩個各向異性異常體的最大和最小主軸電阻率的空間展布一致.因此，可以將圖6a中C1和R2的位置歸為一個各向異性帶ani1，它的模型參數(shù)可以表示為ρx～3 Ωm，ρy～300 Ωm，α～45°.因為前文提到，已經(jīng)將觀測系統(tǒng)逆時針旋轉(zhuǎn)55°至測線方向，那么再將觀測系統(tǒng)順時針旋轉(zhuǎn)55°至正南北-正東西坐標系統(tǒng)，就可以知道ani1區(qū)域里低阻主軸電阻率所在軸的方位呈北偏西10°，如圖7所示.

ani1與圖4a中iso1的形態(tài)接近，前者的整體埋深略淺于后者.另外，在柴達木盆地中心部位的中下地殼，也存在明顯的各向異性體ani2(圖6a中C3所在區(qū)域)，它的模型參數(shù)可以表示為ρx～3 Ωm，ρy～500 Ωm，α～ 0°.圖6a中ρz分量的電性結(jié)構(gòu)在淺層次有所變化，但變化程度小，這主要也與前文提到的原因相關，即二維反演幾乎無法分辨ρz.

在數(shù)據(jù)擬合RMS方面，二維各向同性反演的總體RMS值要小于各向異性反演的結(jié)果，這與各向異性反演中考慮了全部阻抗張量元素有關.如圖8a所示，對全阻抗張量各向異性反演的結(jié)果(圖6a—d)另外計算每個測點對阻抗張量非對角元Zxy和Zyx的RMS，其擬合程度與針對非對角元Zxy和Zyx的各向同性反演(圖4a)的RMS擬合程度非常接近，這一定程度說明，各向異性反演過程在對Zxy和Zyx充分擬合的同時，也考慮到了主對角元素Zxx和Zyy的信息.從圖8b來看，各向異性反演結(jié)果中，對Zxx和Zyy的整體擬合效果不如Zxy和Zyx，但是在測線兩端和中間部位，即測點1—3、8—14以及18—27的區(qū)間，二者擬合效果相當，而且這三個位置大致與圖6a—d中的各向異性結(jié)構(gòu)ani1和ani2的位置對應，因此，Zxx和Zyy的加入有利于反演過程中還原各向異性結(jié)構(gòu).

為進一步驗證各向異性結(jié)構(gòu)的穩(wěn)健性，我們以圖4a中的各向同性結(jié)構(gòu)模型為基礎，將該模型大于10.9 km深度的電阻率重新置換為100 Ωm，形成新的初始模型，采用與圖6a—d相同的各向異性反演參數(shù)設置.在反演開始時，數(shù)據(jù)擬合RMS為3.5369，明顯高于圖4a的1.3860.在經(jīng)過37次迭代后，RMS降為1.7133.反演結(jié)果如圖6e—h所示，ρx和ρy剖面中還原的ani1和ani2區(qū)域的各向異性體與圖6a—d同樣位置的異常特征基本一致.因此，這個結(jié)果說明各向異性帶ani1以及ani2能夠被觀測數(shù)據(jù)所分辨，并且不能被各向同性體所替代.

圖4 剖面L15在不同正則化參數(shù)組合下的各向異性反演模型(a)—(c)分別為對阻抗Zxy和Zyx分量、單獨的Zxy分量以及單獨的Zyx分量的二維各向同性反演結(jié)果；(d)—(g)正則化參數(shù)組合為λs= 0.1和λa=0.1的各向異性反演結(jié)果；(h)—(k)正則化參數(shù)組合為λs=0.1和λa=1.0的各向異性反演結(jié)果；(l)—(o)正則化參數(shù)組合為λs=0.1和λa=10.0的各向異性反演結(jié)果；(p)—(s)正則化參數(shù)組合為λs=0.1和λa=100.0的各向異性反演結(jié)果；(t)—(w)正則化參數(shù)組合為λs=0.1和λa=1000.0的各向異性反演結(jié)果.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.

圖5 測線L15在不同各向異性正則化參數(shù)下的反演結(jié)果統(tǒng)計，其中粗糙度正則因子固定為λs=0.1(a) RMS值； (b) 總目標函數(shù)值； (c) 粗糙度項目標函數(shù)值； (d) 各向異性約束項目標函數(shù)值.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.

5 討論

(1)Pek等(2011)在二維各向異性反演研究中指出，因為ρz對大地電磁正演響應的影響非常微弱，ρz分量也幾乎不能被分辨，在數(shù)據(jù)中加入傾子等與垂直磁場相關的傳輸函數(shù)也無法改善對這個分量的還原情況，各向異性反演中成像的ρz結(jié)構(gòu)主要受目標函數(shù)中各向異性約束項的控制.本文中二維理論模型的反演算例對ρz還原也很不理想，因此，即使我們在實際的反演計算中得到較滿意的數(shù)據(jù)擬合程度，但要進一步還原和解釋ρz分量結(jié)構(gòu)還需尋找其他對ρz敏感的地球物理資料或者地質(zhì)結(jié)構(gòu)信息來支持.

(2)電各向異性是通過電導率張量來表征，而電導率張量可以通過不同的各向異性參數(shù)組合來形成，例如，方位各向異性參數(shù)組(ρx=10 Ωm，ρy=1000 Ωm，α=-45°)和(ρx=1000 Ωm，ρy=10 Ωm，α=45°)實際上代表的是相同的各向異性結(jié)構(gòu)，這可能會導致在反演中，我們還原出了正確的各向異性結(jié)構(gòu)，但是卻由等效的各向異性參數(shù)組表述.在對二維理論模型的還原結(jié)果中，也出現(xiàn)了這種各向異性參數(shù)等效效應.因此，在對各向異性結(jié)構(gòu)進行解釋時，要將各向異性參數(shù)當作一個整體看待，避免結(jié)構(gòu)認識上的偏差.

(3)在各向異性反演中，淺部冗余結(jié)構(gòu)會壓制對深部各向異性體的分辨.在對理論二維各向異性模型的恢復中，可以觀察到，隨深度增加，對各向異性體的分辨能力逐漸下降，很難得到完整的恢復，同時，淺部會呈現(xiàn)出冗余結(jié)構(gòu)，這種淺層冗余構(gòu)造同樣也是各向同性反演中的難題，它的出現(xiàn)主要是因為反演問題的多解性，而考慮各向異性又顯著提升了多解性的程度，因此，通過壓制淺部冗余結(jié)構(gòu)來提升分辨深部異常體的能力，是將來各向異性反演中需要深入研究的方向.

(4)對剖面L15的觀測數(shù)據(jù)重新進行大地電磁各向異性反演，發(fā)現(xiàn)在祁漫塔格山下方50～100 km的區(qū)域存在方位各向異性體ani1.該異常體的水平低阻主軸電阻率空間展布在北偏西10°的軸線上，水平高阻主軸電阻率空間展布在北偏東80°的軸線上.在剖面L15鄰近區(qū)域已有的遠震寬頻橫波分裂研究結(jié)果(Herquel et al.，1995，1999； Soto et al.，2012；Wu et al.，2015)顯示，在阿爾金斷裂、昆侖山斷裂以及柴達木盆地周邊，測量的大部分快波偏振方向與主要的地質(zhì)構(gòu)造走向相吻合(McNamara et al.，1994；Guilbert et al.，1996).已有實例表明，大多數(shù)的快波偏振方向是沿著山脈和其分支走向的(Silver，1996)，而快波偏振方向可以用來指示剪切帶方向，這說明大多數(shù)深部剪切帶與地表構(gòu)造帶可能受到相同變形機制的控制，少數(shù)快波偏振方向與地表構(gòu)造不一致則暗示可能兩者受到更復雜的應力環(huán)境影響.

圖6 驗證剖面L15深部各向異性結(jié)構(gòu)存在性，結(jié)果均為利用全部阻抗張量數(shù)據(jù)對模型參數(shù)ρx、ρy、ρz和α進行各向異性反演獲得(a)—(d)初始模型為100 Ωm半空間；(e)—(h)初始模型是在圖4a模型的基礎上，圖5a深度10.9 km以上保持不變，將圖5a深度10.9 km以下的區(qū)域用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區(qū)域的各向異性異常水平主軸電阻率的空間展布實線表示正南北-正東西的坐標系統(tǒng)；虛線表示沿測線方向的坐標系統(tǒng)，其中y軸是沿測線展布方向；點虛線表示ani1區(qū)域各向異性體的水平主軸電阻率ρx和ρy的空間展布，最終低阻主軸電阻ρx在北偏西10°的軸線上.Fig.7 Spatial distribution of the horizontal principal resistivities of the anomaly in the ani1 area of Fig.6The 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.

圖8 不同二維反演結(jié)果的測點RMS對比(a) 各向同性與各向異性反演結(jié)果中次對角元素擬合RMS的對比；三角形：圖4a各向同性模型中對Zxy和Zyx計算的RMS值；圓形：圖6a—d各向異性模型中對Zxy和Zyx計算的RMS值. (b)圖6a—d各向異性反演結(jié)果中，對不同阻抗張量元素組合的擬合RMS對比.方形：對Zxy和Zyx計算的RMS值；菱形：對Zxx和Zyy計算的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.

關于地震各向異性與電各向異性間的聯(lián)系，現(xiàn)在仍然存在爭議，當測量的快波偏振方向與電各向異性主軸方向吻合時，有研究認為兩者產(chǎn)生于相同的原因(Collins et al.，2006；Frederiksen et al.，2006；Padilha et al.，2006；Carcione et al.，2007；Jones et al.，2009)；而也有一些研究實例認為這兩者之間沒有明顯的關系(Hamilton et al.，2006).地震各向異性與電各向異性之間是否存在關系，這還需在將來通過更多的研究工作去尋找和驗證.

剖面L15成像的祁漫塔格山下各向異性體ani1所在位置屬于上地幔頂部，Caricchi等(2011)在解釋上地幔電各向異性成因時指出，巖石圈下部的部分熔融通道在簡單剪切環(huán)境下會產(chǎn)生電各向異性.其中，在剪切作用的垂直方向上由于存在應力梯度會降低熔融連通性，產(chǎn)生高阻；而沿剪切方向產(chǎn)生相對高連通性，表現(xiàn)出低阻.

而我們在前期對此次區(qū)域的三維各向同性反演研究中(Xiao et al.，2018；Yu et al.，2020)，將剖面L15下的上地幔低阻體解釋為熔融的弱物質(zhì)流，是源自松潘—甘孜下地殼流穿過祁漫塔格山脈向柴達木盆地下方運動的新證據(jù)，并且該處鄰近區(qū)域位于多個剪切帶交匯處(阿爾金斷裂帶，祁漫塔格造山帶，昆侖山斷裂帶等).根據(jù)Caricchi等(2011)對地幔電各向異性成因的解釋，圖6中的各向異性體ani1則可能是該弱物質(zhì)流運動穿過祁漫塔格山脈下方時與剪切運動相互作用而形成.

這里ani1指示的剪切作用位于北偏西10°的軸線上，該軸線方位與Soto等(2012)在幾乎相同位置的地震各向異性結(jié)果指示的北偏東45°左右的剪切帶方位不一致，這很可能與兩種各向異性結(jié)構(gòu)在深度上的變化和跨度不一致相關；另外，ani1指示的剪切作用方位與祁漫塔格山脈走向也不一致，根據(jù)這兩者間的關系進行推論：深部剪切作用方位與地表造山帶地勢的不一致暗示，祁漫塔格山脈地下應力環(huán)境并不單一，ani1指示的深部剪切作用可能是鄰近區(qū)域剪切運動的綜合體現(xiàn)，祁漫塔格在內(nèi)的鄰近區(qū)域深部剪切應力環(huán)境可能還受到阿爾金斷裂帶中部較大的左行走滑運動的影響.

綜上，剖面L15反演還原的電各向異性結(jié)構(gòu)與地震各向異性結(jié)果指示了不同的剪切作用方位，并且這兩種各向異性的空間結(jié)構(gòu)差異較大，所以無法確認兩者成因之間的關系；此外，上地幔電各向異性體ani1可能是祁漫塔格山脈下方的弱物質(zhì)流受到剪切作用的影響而產(chǎn)生的，其表征的深部剪切作用方位與祁漫塔格地表造山帶走向不一致則暗示該處應力環(huán)境復雜，可能還受鄰近阿爾金斷裂帶左行走滑的影響.

6 結(jié)論

基于大地電磁任意各向異性有限差分正演算法以及成熟的NLCG反演模塊，本文實現(xiàn)了二維大地電磁各向異性反演，理論模型響應的反演算例也驗證了算法的穩(wěn)定性.

利用該反演程序，我們對穿過東昆侖祁漫塔格山—柴達木盆地的一條大地電磁剖面進行了重新反演.該剖面早期的各向同性結(jié)果顯示，在祁漫塔格山脈下方的上地幔頂層存在連續(xù)的低阻帶.在本文中，我們對該剖面進行各向異性反演，并在祁漫塔格下方相同層次的位置發(fā)現(xiàn)了明顯的各向異性體，其水平低阻主軸電阻率所在軸的方位呈北偏西10°，水平高阻主軸電阻率所在軸的方位呈北偏東80°.我們認為這里的各向異性異常是上地幔弱物質(zhì)流在祁漫塔格下方受剪切運動的作用而產(chǎn)生沿垂直于剪切作用的方向和沿剪切作用的方向的熔融物質(zhì)導電性差異而形成的；也暗示這里的應力環(huán)境除了印度—歐亞板塊碰撞帶來的遠程擠壓效應外，還可能受到鄰近阿爾金走滑斷裂帶左行剪切運動的影響.

致謝感謝Gary Egbert教授提供了ModEM的反演代碼.感謝三位匿名審稿人提出的寶貴意見，使本文得以進一步改善.