許樂紅, 譚捍東*, 吳萍萍, 彭淼，王帥軍，姜枚

1 中國地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院，北京 100083 2 中國地震局地球物理研究所，北京 100081 3 中國地震局地球物理勘探中心，鄭州 450002 4 中國地質(zhì)科學(xué)院礦產(chǎn)資源所，北京 100037

0 引言

龍門山斷裂帶形成于印支期松潘—甘孜地塊與揚子地體的碰撞以及逆沖推覆作用(林茂炳，1992)，三條主斷裂呈近南北向分布，其南段位于華南地塊、松潘—甘孜地塊和川滇地塊交界處，緊鄰東西向的鮮水河斷裂(圖1b).區(qū)域內(nèi)地形復(fù)雜、演化歷史及板塊動力學(xué)機(jī)制等問題仍存在很大爭議，因此該區(qū)域成為研究青藏高原物質(zhì)向東運動及青藏高原東緣隆升機(jī)制的重要場所(許志琴等, 1991；Worley and Wilson,1996; Hubbard et al.,2009; Xu et al.,2009).目前青藏高原劇烈隆升的主要機(jī)制有：中下地殼管道流模型(Royden et al.,1997,2008; Clark and Royden,2000)、側(cè)向逃逸模式(Tapponnier et al.,2001)與上地殼縮短模型(Hubbard et al.,2009).近年來隨著地球物理方法的不斷改進(jìn)和多資料綜合分析研究的深入，對傳統(tǒng)經(jīng)典模型的質(zhì)疑也逐漸增加.Clark和Royden(2000)認(rèn)為青藏高原東緣的地殼增厚無法用斷層與褶皺作用縮短上地殼的理論來解釋，大部分抬升區(qū)上地殼沒有沿著大型逆沖斷層方向的顯著縮短(Burchfiel et al.，1995)，高原東緣的地勢升高可能是深部地殼增厚與中央高原下方下地殼疏散的直接結(jié)果.朱介壽等(2017)認(rèn)為由青藏高原腹地羌塘高原流出的黏滯性地殼流，帶動上地殼塊體水平移動，遇四川盆地阻擋后的分層形成了上地殼的隆升與下地殼的增厚.He等(2014)認(rèn)為印度板塊持續(xù)向北和印度—亞洲碰撞引起的東西向的相對壓縮是龍門山逆沖斷裂帶和松潘—甘孜地塊地殼增厚的主導(dǎo)因素.

大地電磁法(Magnetotelluric sounding，MT)，探測深度大、不受高阻層屏蔽，廣泛應(yīng)用于殼幔結(jié)構(gòu)探測領(lǐng)域.橫跨龍門山斷裂帶已有多條大地電磁剖面的成果.Zhao等(2012)利用MT數(shù)據(jù)與其他地質(zhì)和地球物理觀測相結(jié)合，提出松潘—甘孜地塊隆起是由四川盆地和龍門山下方穩(wěn)定巖石圈阻擋的東南向地殼流引起的.張樂天等(2012)分析了青藏高原、龍門山斷裂帶及四川盆地下方電性結(jié)構(gòu)的差異及分布特征.詹艷等(2013)在龍門山中段、南段布設(shè)兩條大地電磁剖面，分析了中段與南段的電性差異.王緒本等(2018)結(jié)合長周期大地電磁數(shù)據(jù)分析了龍門山南段電性結(jié)構(gòu)與中、北段的區(qū)別，推斷蘆山地震形成于深部隱伏殼幔韌性剪切帶向上擴(kuò)展的影響.不同電阻率剖面均顯示松潘—甘孜地塊下方存在高導(dǎo)層，龍門山斷裂帶下方揚子地體呈現(xiàn)整塊高阻特征，但大地電磁常規(guī)帶地形反演，垂向分辨率不高，受體積效應(yīng)影響邊界不夠清晰，對深部構(gòu)造分析可能造成誤差.

近年來大地電磁數(shù)據(jù)與地震數(shù)據(jù)之間的約束、聯(lián)合反演算法已日趨成熟并在實際資料中驗證了其有效性與實用性.吳萍萍等(2021)利用VP/VS波速比地震學(xué)模型約束二維大地電磁反演，并應(yīng)用于張渤地震帶的大地電磁資料處理中，分析了區(qū)域內(nèi)地質(zhì)構(gòu)造及孕震環(huán)境.陳高等(2010)利用基于不同物性參數(shù)隨機(jī)分布網(wǎng)格模型的大地電磁與地震數(shù)據(jù)同步聯(lián)合反演，對黔中隆起地區(qū)古生界分布及構(gòu)造特征進(jìn)行了分析.同時，多種地球物理資料綜合分析在龍門山斷裂帶構(gòu)造解釋中也有進(jìn)展.李大虎等(2015)利用臺網(wǎng)及流動臺站天然地震數(shù)據(jù)，P波到時方法結(jié)合重力、航磁數(shù)據(jù)反演，得到了鮮水河斷裂區(qū)三維P波速度結(jié)構(gòu)，認(rèn)為松潘—甘孜地塊地殼物質(zhì)較軟弱，盆地西緣密度較高、物質(zhì)堅硬.徐嘯等(2016)利用寬角反射、折射地震數(shù)據(jù)與重力數(shù)據(jù)聯(lián)合反演，結(jié)合地殼結(jié)構(gòu)與密度分布，分析了龍門山及鄰近區(qū)域上中下地殼厚度、速度及密度的分布規(guī)律，推斷四川盆地剛性基底西緣因擠壓產(chǎn)生的彎曲應(yīng)力是該地區(qū)抬升的重要條件.目前，大地電磁方法與速度結(jié)構(gòu)綜合分析解釋還未在龍門山南段開展.對于地球深部結(jié)構(gòu)及其地球動力學(xué)研究而言，多資料聯(lián)合解釋、互相驗證很可能是未來的發(fā)展趨勢.

本文實現(xiàn)了一種速度模型約束下的帶地形大地電磁二維反演算法，并應(yīng)用于龍門山南段大地電磁實測數(shù)據(jù)中.基于交叉梯度結(jié)構(gòu)約束方法，利用人工地震高分辨率寬角反射/折射P波的速度結(jié)構(gòu)約束大地電磁的電阻率模型，獲得了相較于常規(guī)反演方法更可靠的電性結(jié)果，為分析龍門山斷裂帶南段的深部結(jié)構(gòu)提供了更準(zhǔn)確的約束.

1 區(qū)域地質(zhì)構(gòu)造和大地電磁數(shù)據(jù)

1.1 區(qū)域地質(zhì)構(gòu)造

龍門山斷裂帶沿走向大致以臥龍—懷遠(yuǎn)、北川—安縣一線為界將其劃分為南、中、北三段(李智武等, 2008; 劉樹根等, 2009).此次研究區(qū)域為龍門山南段，區(qū)域包括整個造山帶及其前陸盆地和松潘—甘孜地塊，橫跨區(qū)域內(nèi)主要構(gòu)造單元，且相對于中北段，其形變歷史相對較晚，現(xiàn)今活動性強(qiáng)(李智武等，2008).區(qū)域包含了三個主構(gòu)造單元，松潘—甘孜褶皺帶、龍門山斷裂帶以及四川盆地，是龍門山南段構(gòu)造發(fā)育的代表性地區(qū).

龍門山三大主斷裂包括：北川—映秀斷裂(YBF)、安縣—灌縣斷裂(GAF)和茂縣—汶川斷裂(WMF).北川—映秀斷裂在龍門山南段為耿達(dá)—隴東斷裂(LDF)；安縣—灌縣斷裂在南段為鹽井—五龍斷裂(WLF)；茂縣—汶川斷裂南段是雙石—大川斷裂(SSF).其中北川—映秀斷裂形成于三疊世早期松潘—甘孜造山帶與揚子地體碰撞期，三疊世晚期松潘—甘孜地塊與上揚子地體發(fā)生陸內(nèi)俯沖，三疊世末期松潘—甘孜地塊左旋走滑逆沖于揚子地體之上形成了松潘—甘孜山與前陸盆地(喬秀夫等, 2012).安縣—灌縣斷裂與茂縣—汶川斷裂的逆沖活動發(fā)生于侏羅紀(jì)-早白堊世，古龍門山形成一個由三條逆沖斷裂帶形成的山體.現(xiàn)今的龍門山經(jīng)歷了中生代地震造山與平靜期山脈剝蝕降低的過程，以及新生代多次古地震瞬時造山活動(喬秀夫等, 2016).龍門山南段的前緣及其縱深地區(qū)——寶興地區(qū)的寶興雜巖，存在著多次推覆和多層滑覆現(xiàn)象，以及推覆滑覆疊加樣式(林茂炳，1992).印支期和喜馬拉雅期的造山運動，使得這一地區(qū)地震頻發(fā)、構(gòu)造極不穩(wěn)定.先后發(fā)生了2013年4月20日蘆山地震以及2014年11月22日康定地震.蘆山地震震源精定位的震中為30.28°N，103.98°E，震源深度16.38 km(王小娜等，2015).初步調(diào)查認(rèn)為，主震源位于雙石—大川斷裂附近，多數(shù)余震分布在主斷裂西南側(cè)區(qū)域(劉杰等，2013).龍門山斷裂南段的速度結(jié)構(gòu)、地殼電性結(jié)構(gòu)分布、斷裂帶產(chǎn)狀和深部延伸特征均與大地震活動密切相關(guān).

1.2 大地電磁數(shù)據(jù)來源

雅安地震后由汶川科鉆項目的“橫穿龍門山地震反射剖面探測及深部結(jié)構(gòu)的地震解釋”課題支持，于2013年至2015年間分三次穿越龍門山斷裂帶進(jìn)行大地電磁測深野外數(shù)據(jù)采集工作.測線南起雅安市嚴(yán)橋鎮(zhèn)，沿北西向經(jīng)過雅安市蘆山縣、寶興縣，瀘定縣，穿過小金縣，向北延伸至阿壩縣大寨村，橫穿四川盆地以及龍門山三大主斷裂帶并深入松潘—甘孜地塊，直線距離近190 km.從四川盆地到青藏高原西南緣，高程變化約4500 m，測點設(shè)計間距約5 km,由于地形復(fù)雜斷裂帶多破碎嚴(yán)重，山地居多伴隨河流且植被豐茂，沿途有水電站、高壓塔等強(qiáng)干擾因素的影響，部分區(qū)域測點較稀疏(測點位置如圖1中黑色三角形所示).采集儀器為加拿大鳳凰公司的MTU-5寬頻大地電磁儀，采集頻率范圍為320～0.0005 Hz，電極采用標(biāo)準(zhǔn)“十”字形布設(shè)，平均采集時間為24 h.

圖1 龍門山斷裂帶南段大地電磁測深剖面位置圖及區(qū)域構(gòu)造情況簡圖(a) 青藏高原及周緣塊體分布簡圖.圖中HLB—喜馬拉雅地塊；LSB—拉薩地塊；QTB—羌塘地塊；S-GB—松潘—甘孜地塊；QDB—柴達(dá)木盆地；QL—祁連地塊; ODB—鄂爾多斯盆地; SCB—四川盆地；YTT—揚子地體；SCFS—華南褶皺系.黑色方框為(b)圖范圍，紅色方框為(c)圖范圍；(b)龍門山斷裂帶南段區(qū)域板塊構(gòu)造簡圖.圖中塊體邊界數(shù)據(jù)來自張培震等(2003)、張國民等(2005)和王輝等(2003).青藏高原及周緣板塊主要運動特征據(jù)Tapponnier等(2001)和李英強(qiáng)(2018)修改；(c)大地電磁和地震測點位置圖.圖中紅色藍(lán)色圓點為該區(qū)域歷史地震位置(2000-11-01—2020-01-01)，藍(lán)色代表MS3.0～4.9之間的地震，共計744次；紅色代表MS5.0～6.9之間的地震，共計52次.地震目錄來自USGS(United States Geological Survey，美國地質(zhì)勘探局).黑色三角形為大地電磁測點位置；黑色圓點為寬角折/反射地震接收臺站位置；紅色五角星為寬角折/反射地震炮點位置.紫色直線為圖12中所參考人工反射地震資料(李英強(qiáng)，2018)測線位置.汶川MS8.0地震、蘆山MS7.0地震、康定MS6.3地震的震源機(jī)制解，數(shù)據(jù)來自GCMT(Global Centroid-Moment-Tensor，全球質(zhì)心矩張量).斷裂帶位置(紅色實線)來自于鄧起東等(2003)，張量分解玫瑰圖顯示電性主軸方位為近北東向.Fig.1 Topography map showing MT stations layout and tectonic structures in the south part of the Longmenshan fault zone(a)Distribution diagram of Qinghai-Tibet Plateau and surrounding plates. HLB—Himalaya Block; LSB—Lhasa Block；QTB—Qiangtang Block；S-GB—Songpan-Garzê Block； QDB—Qaidam Basin；QL—Qilian Block; ODB—Ordos Basin; SCB—Sichuan Basin；YTT—Yangtze Terrane；SCFS—South China Fold System. The black box shows the scope of figure (b), the red box shows the area of figure (c); (b) Regional plate tectonics in the south part of the Longmenshan fault zone. The block boundary data according to Zhang P Z et al.,2003; Zhang G M et al.,2005; Wang et al.,2003. The main movement characteristics of the Qinghai-Tibet Plateau and its surrounding plates were modified from Tapponnier et al., 2001 and Li et al.,2018;(c) Location map of MT and seismic profile. Red and blue circles are the position of the historical earthquake in this region (2000-11-01—2020-01-01). Blue circles: 744 earthquakes between MS3.0～4.9; red circles: 52 earthquakes between MS5.0～6.9. The earthquake catalog comes from USGS. Black triangles represent the MT sites; Black circles represent receivers of the wide-angle seismic reflection/refraction profile. The red stars represent the shot points of the wide-angle seismic reflection/refraction profile; The straight purple line is the artificial seismic reflection profile shown in Fig.12, according to Li et al., 2018. The mechanism solutions of Wenchuan MS8.0, Lushan MS7.0, and Kangding MS6.3 earthquakes are obtained from GCMT. The location of faults (red solid lines) refers to Deng et al.(2003). The rose diagrams of MT tensor decomposition show that the electrical strike direction orients nearly northeast.

1.3 大地電磁數(shù)據(jù)預(yù)處理

野外采集到的大地電磁數(shù)據(jù)，通過預(yù)處理過程剔除信噪比低、誤差較大的頻點及數(shù)據(jù)并將時間序列轉(zhuǎn)換至阻抗張量數(shù)據(jù).首先通過快速傅里葉變換將時間域數(shù)據(jù)轉(zhuǎn)換到頻率域.利用ROBUST估計等技術(shù)對數(shù)據(jù)進(jìn)行篩選，計算得到初始阻抗張量，手工剔除個別畸變的頻點數(shù)據(jù).對于個別頻點或連續(xù)頻點值分布不正常，或者儀器不穩(wěn)的情況，運用RHOPLUS理論對其進(jìn)行視電阻率和相位曲線的一致性校正(譚捍東等, 2004).校正后的數(shù)據(jù)采用MT-pioneer軟件分析電性主軸方向(蔡軍濤等, 2010)，利用相位張量方法計算最佳主軸方向角為北東向16°，部分測點的單臺全頻段玫瑰圖(見圖1)顯示電性主軸方向與龍門山斷裂構(gòu)造走向相符，四川盆地區(qū)域測點及小金地區(qū)個別測點三維性較強(qiáng)，其余測點二維性較好，測線整體二維性較強(qiáng).利用MT-editor軟件將所有測點旋轉(zhuǎn)角度至電性主軸方向上.典型測點視電阻率與相位曲線如圖2所示，其中龍門山至松潘—甘孜地塊區(qū)域，地形起伏變化劇烈，數(shù)據(jù)靜位移現(xiàn)象明顯.

2 速度結(jié)構(gòu)約束下的帶地形大地電磁二維反演方法

自從Gallardo和Meju(2003, 2004)提出基于兩種物性參數(shù)梯度的叉乘獲取兩種物性的結(jié)構(gòu)相似性以來，該方法在各種地球物理反演約束或聯(lián)合反演中廣泛應(yīng)用(彭淼等，2013；Peng et al., 2019; Wu et al., 2020；閆政文等, 2020).但已有的研究成果主要針對不帶地形的交叉梯度結(jié)構(gòu)約束，而四川盆地至青藏高原測區(qū)抬升劇烈，地形起伏高差達(dá)5 km，地形是影響反演效果和精度的主要因素之一(王緒本等，1999).因此本文根據(jù)龍門山區(qū)域?qū)嶋H地形反演的需求開發(fā)了帶地形的交叉梯度約束算法，并將已有的速度模型作為固定約束加到算法中，進(jìn)行約束反演.以下介紹算法的實現(xiàn).

圖2 部分典型大地電磁臺站校正后的實測視電阻率和阻抗相位曲線紅色曲線為TE模式；藍(lán)色曲線為TM模式.Fig.2 Apparent resistivity and impedance phase curves after correction at several typical MT stationsThe red curve represents the TE mode and the blue curve represents the TM mode.

2.1 帶地形的交叉梯度項離散方法

根據(jù)交叉梯度函數(shù)的定義，可以將兩種物性的交叉梯度項寫成：

(1)

mr表示電阻率模型，ms表示速度模型.在x-z坐標(biāo)(二維坐標(biāo)系)中，交叉梯度項分量中只有ty有值，tx和tz為零(吳萍萍，2019).根據(jù)矢量叉積求解公式，可將ty寫成：

(2)

采用中心差分法將公式(2)離散化成：

(3)

其中i表示z軸第i個網(wǎng)格點，j表示x軸方向第j個網(wǎng)格點.Δxj表示x軸方向第j個網(wǎng)格的網(wǎng)格間距，Δzi表示z軸方向第i個網(wǎng)格的網(wǎng)格間距.

在考慮地形的交叉梯度函數(shù)的離散化過程中，主要考慮非邊界網(wǎng)格、地形邊界、兩個側(cè)邊界和底界面，跟不帶地形的交叉梯度函數(shù)離散化主要區(qū)別在于地形界面的處理方式.非邊界網(wǎng)格、兩個側(cè)邊界和底界面的計算公式可以參考吳萍萍(2019).以下主要介紹地形邊界，公式(3)的計算方式.

在二維地形界面中，以網(wǎng)格(C)為中心，該網(wǎng)格四周有4個網(wǎng)格，分別為上(U)、下(D)、左(L)、右(R).本文主要將地形邊界網(wǎng)格分成3種情況：

(1)當(dāng)C網(wǎng)格周邊只有一個空氣網(wǎng)格時，則對應(yīng)的情況有圖3(a)、(b)和(c)三種情況.以(a)為例，簡要說明該網(wǎng)格的ty的計算方式，圖中(i-1,j)為空氣網(wǎng)格，則將公式(3)中有關(guān)(i-1,j)的物性和網(wǎng)格參數(shù)用中心網(wǎng)格(i,j)代替，則公式可以改寫成：

(4)

同理，(b)和(c)的公式可以寫成如下：

(5)

(6)

(2)當(dāng)C網(wǎng)格周邊有兩個空氣網(wǎng)格時，則對應(yīng)的情況有圖3(d)、(e)和(f)三種情況.以(d)為例，當(dāng)(i,j-1)和(i-1,j)為空氣時，公式(3)可以改寫成：

(7)

同理，(e)和(f)的公式可以寫成如下：

(8)

圖3 地形邊界網(wǎng)格示意圖Fig.3 Sketch map of terrain boundary grid

ty(i,j)=0.

(9)

(3)當(dāng)C網(wǎng)格周邊有三個空氣網(wǎng)格時，只有一種情況(圖3(g))，即(i,j-1)和(i-1,j)和(i+1,j)為空氣網(wǎng)格，則該網(wǎng)格的ty可以寫成：

ty(i,j)=0.

(10)

在約束反演過程中，ms是已知的速度模型，為常數(shù).將交叉梯度函數(shù)的所有網(wǎng)格組合成一個線性矩陣，

ty=Wcg,smr,

(11)

其中mr中只包含非空氣的介質(zhì)電阻率，Wcg,s為常系數(shù)矩陣.

取ty的平方和作為結(jié)構(gòu)耦合約束項，則約束項可表示為：

(12)

2.2 基于速度結(jié)構(gòu)約束的帶地形大地電磁反演算法

帶地形大地電磁反演目標(biāo)函數(shù)由數(shù)據(jù)項和模型項組成，公式如下：

f(·)是二維大地電磁正演響應(yīng)函數(shù)，λ為拉格朗日因子，Cd為數(shù)據(jù)協(xié)方差矩陣，Cm為模型光滑矩陣，dobs為觀測數(shù)據(jù)矩陣.

將交叉梯度項通過權(quán)重因子η加到無約束反演目標(biāo)函數(shù)中，形成速度模型約束下的帶地形大地電磁二維反演目標(biāo)函數(shù)，具體表示如下：

+η‖Wcg,smr‖2.

(14)

對目標(biāo)函數(shù)(14)求導(dǎo)數(shù)，可得目標(biāo)函數(shù)的梯度：

(15)

其中J為雅可比矩陣,它是正演響應(yīng)函數(shù)f(mr)對模型mr的一階偏導(dǎo)數(shù)矩陣.采用非線性共軛梯度法對加有交叉梯度約束的目標(biāo)函數(shù)(14)進(jìn)行最小化.在求取公式(15)所示的梯度時，該方法避免了直接求解和存儲雅可比矩陣J，先通過解一次擬正演問題求雅可比矩陣的轉(zhuǎn)置和一個向量的乘積，即JT(dobs-f(mr))，進(jìn)而再按照公式(15)去計算梯度，避免了每次迭代都計算完整的大型矩陣，減少了計算時間和內(nèi)存需求.圖4為本文采用的速度結(jié)構(gòu)約束下大地電磁二維帶地形反演算法流程圖.

本文在傳統(tǒng)的大地電磁二維非線性共軛梯度(NLCG)反演(Rodi and Mackie, 2001)基礎(chǔ)上，設(shè)計了帶地形的網(wǎng)格，將地震資料插值到網(wǎng)格得到約束矩陣，引入了帶地形的交叉梯度項，形成了速度結(jié)構(gòu)約束的帶地形大地電磁二維反演算法.

圖4 速度結(jié)構(gòu)約束下大地電磁二維帶地形反演算法流程圖Fig.4 Flow chart of the 2D MT inversion algorithm in consideration of both the terrain and the velocity structure constraints

3 理論模型合成數(shù)據(jù)的反演算例

為了驗證本文算法的可靠性和有效性，設(shè)計了一個長度為280 km的地塹-地壘推覆構(gòu)造模型(圖5)，地形起伏8 km，沿地表均勻布設(shè)46個測點，測點間距6 km，圖中以黑色三角形表示.電阻率模型共分為五部分，A部分為高阻高速異常體，A 塊體ρA=1000 Ωm，vA=6.0 km·s-1；B部分為低阻低速異常體，B塊體ρB=1 Ωm，vB=4.0 km·s-1；C部分為高阻低速異常體，C塊體ρC=400 Ωm，vC=4.0 km·s-1；D部分為低阻高速異常體，D塊體ρD=10 Ωm，vD=5.5 km·s-1.速度模型中單獨設(shè)計了一塊高速異常體E，vE=6.5 km·s-1，此區(qū)域無電阻率異常.A塊體上涌推覆至B塊體之上，背景場值ρ=100 Ωm，v=5.0 km·s-1.設(shè)計帶地形的網(wǎng)格，正演計算出頻段0.00057～320 Hz之間40個頻點的TE、TM數(shù)據(jù)，加入5%的高斯誤差作為觀測數(shù)據(jù).

反演網(wǎng)格橫向采用等間隔剖分，網(wǎng)格間距2 km，縱向為不等間隔剖分，地形部分采用200 m間隔剖分，地形以下按1.5倍間距逐漸增大.初始模型選用電阻率為100 Ωm的均勻半空間.反演結(jié)果如圖6所示，其中圖6a為大地電磁單方法結(jié)果，可見帶地形的反演方法效果較好，無明顯由地形起伏引入的假異?，F(xiàn)象，推覆體構(gòu)造部分形態(tài)、位置與模型基本一致.但深部體積效應(yīng)明顯，高低阻分界面及輪廓圓滑不清晰，高阻體相鄰兩塊體反演結(jié)果有連通的假象.圖6b為速度結(jié)構(gòu)約束下的大地電磁帶地形交叉梯度反演結(jié)果，相較于單方法，高阻體輪廓恢復(fù)較好，與模型區(qū)域基本一致，階梯狀推覆構(gòu)造清晰，高低阻體分界面分辨率有所提高.其中對于塊體C與A，常規(guī)帶地形反演結(jié)果有部分連通趨勢，通過約束反演后，分離為單個獨立塊體.對于低阻體部分，塊體D輪廓有所改善，塊體B頂部輪廓恢復(fù)較好.對于速度異常體E(圖5b)，此處大地電磁結(jié)果無異常，約束后的結(jié)果在該區(qū)域內(nèi)仍然未出現(xiàn)異常，說明約束反演只作用于大地電磁有異常區(qū)域，不會引入虛假異常.電阻率異常體對應(yīng)的速度異常體無論是高速還是低速，均有較好的約束效果.

圖5 電阻率及速度結(jié)構(gòu)的理論模型圖(a) 電阻率模型圖； (b) 速度模型圖.Fig.5 The theoretical model plots of resistivity structure and velocity structure(a) The resistivity model; (b) The velocity model.

圖6 理論模型合成數(shù)據(jù)的大地電磁二維反演結(jié)果對比(a) 常規(guī)帶地形反演； (b) 本文算法進(jìn)行的約束反演結(jié)果.Fig.6 The comparison of 2D MT inversion results of synthetic data from theoretical models(a) The regular inversion result with terrain; (b) The constrained inversion result of the algorithm in this study.

圖7a為大地電磁常規(guī)帶地形反演的L-curve曲線圖，采用了多個正則化因子τ進(jìn)行多次試算，計算所得RMS與粗糙度分布規(guī)律如圖所示，最終選定處于曲線拐點位置τ=10時的反演結(jié)果.圖7b為約束反演過程中，權(quán)重系數(shù)η=1×1016時，約束反演與無約束反演的RMS下降曲線對比，約束反演迭代次數(shù)增加，最終RMS值略高于常規(guī)帶地形反演.圖7c為無約束反演電阻率模型與速度模型的交叉梯度值空間分布圖，可見模型邊界與兩種方法沒有共同異常體的區(qū)域E，這兩部分交叉梯度值較大，即常規(guī)帶地形反演結(jié)果與速度模型在邊界區(qū)域有較大的非零值.圖7d為約束反演與速度模型的交叉梯度值空間分布圖，圖中交叉梯度值在異常體邊界部分明顯減小.約束反演有效的提升了大地電磁反演在輪廓邊界位置的精確度.同時，參考Moorkamp 等(2011)論文，對于反演模型分別做了模型還原度求解，其公式如下：

(16)

圖7 反演RMS曲線及交叉梯度值空間分布圖(a) 大地電磁反演L-curve曲線； (b) 無約束和有約束反演RMS迭代曲線對比； (c) 無速度模型約束的反演結(jié)果交叉梯度分布圖; (d) 速度結(jié)構(gòu)約束的反演結(jié)果交叉梯度分布圖.Fig.7 RMS curves in the inversion and the spatial distribution of the cross-gradient values(a) The L-curve subplot for MT inversion; (b) The RMS iteration curves for no-constrained and constrained inversion; (c) The spatial distribution map of the cross-gradient values from the inversion results without the velocity structure constraints; (d) The spatial distribution map of the cross-gradient values from the inversion results with the velocity structure constraints.

4 小金—雅安剖面實測數(shù)據(jù)的反演

4.1 P波速度結(jié)構(gòu)

利用P波速度模型約束實測數(shù)據(jù)的大地電磁反演中充分考慮了兩種方法在分辨率以及勘探深度、測線區(qū)域的匹配度.對比分析后認(rèn)為寬角折/反射地震資料與大地電磁方法更匹配.我們選擇了中國地震局地球物理勘探中心在蘆山震后布設(shè)的一條長約410 km的人工地震寬角折/反射探測剖面的走時反演結(jié)果(王帥軍等, 2015)，作為約束大地電磁數(shù)據(jù)的速度結(jié)構(gòu).地震觀測系統(tǒng)，由沿剖面的8個炮點激發(fā)和268臺三分量PDS-2型地震儀器同步接收組成，接收臺間距在0.8～2 km之間.對觀測資料進(jìn)行了二維非均勻介質(zhì)動力學(xué)射線追蹤、走時擬合計算，最終獲得二維殼幔速度結(jié)構(gòu)模型.其分辨率約為1 km×1 km，勘探深度達(dá)到70 km，與大地電磁深度大致相當(dāng).由于地震測線較長，故截取其中長度約240 km的P波速度數(shù)據(jù)(圖8)，保證全部覆蓋了大地電磁測點區(qū)域且向測線兩側(cè)各有幾十公里延伸.依據(jù)大地電磁方法的反演網(wǎng)格，采用不等間隔插值，將速度剖面數(shù)據(jù)插值到大地電磁反演網(wǎng)格中，約束大地電磁反演.

如圖8所示，速度模型呈層狀結(jié)構(gòu)，Moho界面清晰.松潘—甘孜地塊與揚子地體速度結(jié)構(gòu)迥異，并以龍門山三大主斷裂為界.松潘—甘孜地塊上覆蓋層已部分出露地表，20 km以上至地表附近有一個低速異常區(qū)，且中下地殼速度以階梯狀垂向增大，其Moho界面深度在60 km左右.揚子地體地表下沉積蓋層較厚，約7～8 km，呈低速特征；Moho界面深度淺，約50 km左右.波速異常區(qū)域與大地電磁電阻率異常區(qū)域?qū)?yīng)較好.

圖8 沿大地電磁測線的P波速度模型截面Moho面深度及各區(qū)塊分布位置參考王帥軍等(2015).Fig.8 The section of VP model along the MT profileThe Moho depth and the distribution of each block refers to Wang et al.(2015).

4.2 速度結(jié)構(gòu)約束下的大地電磁二維反演

反演網(wǎng)格部分，縱向地形區(qū)域以50 m等間隔剖分，不參與計算的地形網(wǎng)格，采取將實測高程數(shù)據(jù)插值到網(wǎng)格節(jié)點的方法來剔除(圖9中紅色區(qū)域)，地形以下區(qū)域網(wǎng)格間距以1.5倍逐漸增大.橫向網(wǎng)格剖分采用不等間隔，測點附近加密，測點兩側(cè)以1.5倍逐漸加大.在保證了模型計算邊界縱向區(qū)域足夠深、橫向跨度足夠大的前提下，大大減少了網(wǎng)格塊的數(shù)量，從而提高了正反演計算效率，減少了內(nèi)存需求.反演的初始模型參考同區(qū)域內(nèi)大地電磁常規(guī)帶地形反演結(jié)果(張樂天等，2012; 王緒本等，2018)，并通過反復(fù)試算后確定選用電阻率為100 Ωm的均勻半空間模型.

加入地形網(wǎng)格后的單獨大地電磁反演結(jié)果如圖10b所示，可見區(qū)域內(nèi)三大主體構(gòu)造走向、位置清晰，電阻率值過渡光滑.松潘—甘孜地塊下方40 km以上普遍存在低阻異常，低阻體以逆沖形態(tài)延伸至龍門山斷裂帶下方大塊高阻體之上，龍門山深部高阻異常分布于70 km以上，四川盆地淺層5 km左右為低阻異常帶，深部20～60 km左右存在低阻異常區(qū).三大塊體具有明顯“低-高-低”的宏觀電性特征，常規(guī)帶地形反演結(jié)果與王緒本等(2018)在該區(qū)域內(nèi)長周期大地電磁反演結(jié)果比對，60 km以上部分主體構(gòu)造分布基本一致，說明單方法反演結(jié)果可靠.加入P波速度模型，進(jìn)行交叉梯度約束反演后的結(jié)果如圖10c所示，約束反演后的結(jié)果有效克服了大地電磁單方法邊界收斂不好的情況，高低阻異常體分布區(qū)域及輪廓、走向均有變化.松潘—甘孜塊體中下地殼低阻體從50 km至地表均有顯示，低阻區(qū)域深度有所增加，兩塊異常體有連通通道.低阻體向龍門山下方高阻體頂部逆沖并深入地表蓋層之下，部分已出露至地表.龍門山斷裂帶高阻體范圍較單方法有了明顯的收斂，底部邊界收斂至50 km左右深度，Moho界面上部，等值線傾斜方向顯示與人工地震資料中龍門山三條主斷裂的深部延伸產(chǎn)狀一致(圖12)，更加符合真實地質(zhì)情況.四川盆地淺地表低阻層約束后薄層特征明顯，淺部低阻體與深部低阻體中間有一條形高阻帶，高低阻過渡清晰，斷裂帶構(gòu)造明顯.中下地殼大范圍低阻體分布，其中西側(cè)緊鄰龍門山高阻體的低阻異常區(qū)域向深部延伸至Moho界面以下，東側(cè)低阻體下部為高阻區(qū)，低阻分布于40 km以上.

圖9 大地電磁反演網(wǎng)格剖分藍(lán)色三角形為測點位置，紅色填充部分為地形網(wǎng)格.Fig.9 Mesh discretization for the MT inversionBlue triangles represent the MT sites，and the red fills are terrain grids.

圖10 實測數(shù)據(jù)大地電磁二維反演結(jié)果(a) 實際地形剖面； (b) 常規(guī)帶地形反演結(jié)果； (c) 本文算法進(jìn)行的約束反演結(jié)果.Fig.10 2D MT inversion results of the measured data(a) The actual topographic section; (b) The regular inversion results with terrain; (c) The constrained inversion results of the algorithm in this study.

圖11 反演RMS曲線及交叉梯度值空間分布圖(a) 大地電磁反演L-curve曲線，紅色空心圓表示不同τ的RMS值，藍(lán)色五角星表示最終選擇值； (b) 無約束和有約束反演RMS迭代曲線； (c) 無速度模型約束的反演結(jié)果交叉梯度分布圖; (d) 速度結(jié)構(gòu)約束的反演結(jié)果交叉梯度分布圖.Fig.11 RMS curves in the inversion and the spatial distribution of the cross-gradient values(a) The L-curve subplot for MT inversion; The red hollow circle curve represents the RMS value for different τ, the blue pentagram represents the optimal value; (b) The RMS iteration curves for no-constrained and constrained inversion; (c) The spatial distribution map of the cross-gradient values from the inversion results without the velocity structure constraints; (d) The spatial distribution map of the cross-gradient values from the inversion results with the velocity structure constraints.

圖12 綜合解釋圖地形數(shù)據(jù)來自GMT(Generic Mapping Tools)軟件，全球地形數(shù)據(jù)包earth_relief；斷裂帶深部產(chǎn)狀及延伸位置依據(jù)人工反射地震結(jié)果(李英強(qiáng)，2018)修改；地表斷裂帶位置(紅線)來自于鄧起東等(2003)；沙灘球為蘆山MS7.0地震震源機(jī)制解，數(shù)據(jù)來自于GCMT.Fig.12 Integrated interpretation mapTerrain data comes from the global terrain data package (earth_relief) of GMT software. The fault geometry and the extended position underground are modified based on the results of artificial seismic reflection profiles (Li,2018); The location of faults (red lines) refers to Deng et al.(2003); The beachball is the focal mechanism solutions of Lushan MS7.0 earthquake obtained from GCMT.

反演選定的參數(shù)、迭代次數(shù)與對應(yīng)的RMS值如圖11a、b所示，常規(guī)帶地形反演試算多個正則化因子，最終選取τ=10的結(jié)果，因為曲線沒有明顯拐點，τ=10與τ=5相比，RMS差別不大但粗糙度相差較大，故選擇τ=10的結(jié)果更光滑.約束前后RMS曲線對比如圖11b所示，RMS值略有增加但差別不大，說明擬合情況較好.圖11c與圖11d分別為約束反演前、后電阻率模型與速度模型的交叉梯度值的空間分布圖，可見約束后整體梯度值都有下降，其中上地殼淺層結(jié)構(gòu)改善顯著.

5 構(gòu)造解釋分析討論

龍門山南段區(qū)域斷裂帶展布較北部更寬，分支斷裂多.本次研究測線自西向東穿過的主要斷裂帶有：磽磧斷裂(QQF)、隴東斷裂(LDF)、五龍斷裂(WLF)、小關(guān)子組斷裂(XGZF)、雙石—大川斷裂(SSF)、新開店斷裂(XKDF)、大邑隱伏斷裂(DYF).以雙石—大川斷裂為界，西側(cè)為韌性較強(qiáng)的基底卷入構(gòu)造，東側(cè)為基底卷入變形隱伏于沉積蓋層之下(李英強(qiáng)，2018)，磽磧鄉(xiāng)至小金段分布多條弧形斷裂.

5.1 松潘—甘孜褶皺帶逆沖推覆構(gòu)造形成機(jī)制

松潘—甘孜褶皺帶位于青藏高原東部，由北部的勞亞板塊、西部的昌都—羌塘微板塊向上仰沖，東部的揚子地體向下俯沖而成(許志琴等，1991).眾多學(xué)者研究表明，松潘—甘孜地塊巖石圈普遍存在低阻異常帶，但低阻層分布不均.對比龍門山中、北段大地電磁結(jié)果(詹艷等，2013；彭淼等，2015)，低阻層位置及規(guī)模與南段相比差異較大，北段低阻層更接近于地表且深度淺，分布在15 km以上，地表出露較多.中段低阻層規(guī)模大，深度更深，位置大約在20～50 km左右，地表出露幾乎沒有.南段低阻層位置較淺，規(guī)模更大，淺地表至50 km范圍內(nèi)均有分布，西側(cè)小金附近有部分出露地表.中部、南部低阻體逆沖趨勢明顯，覆蓋于揚子地體表層并向東側(cè)有一定延伸.

深部地球物理探測中，脫水和部分熔融可用于解釋殼內(nèi)低速-高導(dǎo)層的存在(Meissner and Wever,1986).巖石實驗結(jié)果證實，少量的熔融體便可以大量吸收彈性波能量，造成地震波速的明顯下降.Sato和Sacks(1989)發(fā)現(xiàn)地幔橄欖巖中只要有5%的熔融體存在就可以導(dǎo)致地震波速下降約5%.結(jié)合該區(qū)域較高的縱橫波速度比以及存在的殼內(nèi)負(fù)極性震相，說明青藏高原東緣地殼內(nèi)部存在低黏滯性通道(Zhang et al., 2009)，支持“下地殼流”或“構(gòu)造逃逸”模型.郭曉玉等(2014)推斷北段揚子板塊邊界位于龍日壩斷裂帶附近，王緒本等(2018)推測南段揚子板塊西緣位于大地電磁結(jié)果中的高低阻過渡區(qū).從大地電磁反演結(jié)果來看，南段邊界符合位于高低阻分界面的判斷，但并未有中上地殼物質(zhì)深入揚子地體下方的趨勢，另外一支熱流應(yīng)轉(zhuǎn)向南流入三江塊體(熊熊等，2001).綜合分析后，推斷龍門山南段逆沖推覆構(gòu)造，形成于下地殼流向東運移，造成了青藏高原隆升，在側(cè)向擠壓揚子地體過程中，向東側(cè)運移的熱流分支逆沖至揚子地體之上；另一分支向南移動形成弧形斷裂帶.龍門山斷裂帶的產(chǎn)生源于印度板塊與歐亞大陸發(fā)生陸陸碰撞造成青藏高原深部物質(zhì)運移，運移方向如圖1b所示.

5.2 龍門山斷裂帶南段構(gòu)造特征

龍門山三條主斷裂延伸至南段，寬度逐漸加大，次級斷裂增多，破碎嚴(yán)重.大地電磁反演結(jié)果顯示，薄皮蓋層下方高阻塊體厚約50 km，三大斷裂均分布于高阻體內(nèi)部，沒有高、低阻過渡區(qū).寬角折/反射資料顯示，轉(zhuǎn)換帶內(nèi)部震相紊亂、不清晰，速度結(jié)構(gòu)異常呈現(xiàn)低速異常特征(王帥軍，2015).人工地震資料顯示，斷裂帶呈疊瓦狀排列，交匯于地下20 km處逆沖構(gòu)造滑脫面，符合逆沖推覆特征.

地震精定位結(jié)果顯示，蘆山地震震中位于蘆山向斜下方，雙石斷裂與新開店斷裂附近，但震源深度更深，電阻率結(jié)果中，該區(qū)域為靠近高阻塊體邊緣的位置，附近沒有汶川地震類似的低阻區(qū)，人工地震資料顯示大邑斷裂帶東側(cè)的隱伏斷裂從震源附近穿過.穩(wěn)定的揚子地體，阻擋了松潘—甘孜地塊內(nèi)的物質(zhì)運移.下地殼流沿大地電磁低阻異常帶運移，推測西側(cè)淺層滑脫面沿高低阻過渡帶縱向深入至滑脫層.在高原物質(zhì)運移過程中，擠壓揚子地體，驅(qū)動了蘆山地震的發(fā)生.

5.3 四川盆地西南緣構(gòu)造特征

測線位置穿過四川盆地南端邊緣，此處的沉積規(guī)律與北部基本相符卻更復(fù)雜.該區(qū)域分布有多處地?zé)豳Y源，推測深部有熱源上涌通道.反演結(jié)果顯示淺地表低阻體分布范圍約5～10 km.低阻異常區(qū)域內(nèi)，P波、S波速度結(jié)構(gòu)均為明顯的低速異常，Liu等(2014)推斷該區(qū)域為四川盆地上三疊紀(jì)-第四紀(jì)沉積物的反映.深部低阻帶分為兩個區(qū)域，靠近龍門山的低阻帶有向深部延伸的趨勢，推斷為川西殼幔低阻帶；測線東側(cè)接近邊緣部分的低阻體下方為高阻體，推測為川中高阻體邊緣.地表低阻體與中上地殼低阻帶之間有一層高阻過渡區(qū)域，依據(jù)人工地震結(jié)果劃分的隱伏斷裂位置，可見多條隱伏斷裂存在于此，斷裂帶沿電阻率變化界面分布，推斷蘆山地震可能由隱伏斷裂驅(qū)動產(chǎn)生.

6 結(jié)論

本文利用人工地震資料反演獲得的P波速度約束大地電磁反演，有效提高了大地電磁方法的縱向分辨率.新算法應(yīng)用于龍門山南段地區(qū)大地電磁實測資料中，得到了該區(qū)域更為可靠的電性結(jié)構(gòu)分布.本次研究表明：

(1)大地電磁方法體積效應(yīng)明顯，縱向分辨率低，利用速度模型層狀結(jié)構(gòu)及邊界清晰的特點，可以有效改善大地電磁方法的體積效應(yīng)，更清晰刻畫異常體輪廓；

(2)帶地形的約束反演方法，尤其適用于地形起伏變化劇烈的研究區(qū)域開展多資料反演和解釋；

(3)推斷龍門山南段逆沖推覆構(gòu)造，形成于下地殼流向東運移，造成了青藏高原隆升并側(cè)向擠壓揚子地體，向東側(cè)逃逸的熱流分支逆沖至揚子地體之上，另一分支向南移動形成附近弧形斷裂帶；

(4)推斷松潘—甘孜地塊低阻異常區(qū)下方為深部滑脫層，高原物質(zhì)運移過程中，擠壓揚子地體，驅(qū)動了蘆山地震的發(fā)生.

本文提出的基于地震速度約束的大地電磁反演方法尤其適用于類似龍門山等地形起伏劇烈研究區(qū)域.反演結(jié)果結(jié)合已有物探資料分析了龍門山南段與北段電性結(jié)構(gòu)的差異，推斷出揚子板塊前緣位于高低阻過渡帶附近，分析了蘆山地震由隱伏斷裂驅(qū)動，為龍門山南段的地質(zhì)構(gòu)造解釋提供了新資料.