馬小軍 吳慶舉 潘佳鐵 鐘世軍 徐 薈

1)中國地震局地球物理研究所，北京 100081

2)寧夏回族自治區(qū)地震局，銀川 750001

3)北京市地震局，北京 100080

0 引言

解釋青藏高原地殼縮短增厚成因的動(dòng)力學(xué)模型包括：1)剛性塊體模型，該模型認(rèn)為高原下部的物質(zhì)沿著大型斷裂向E擠出(Molnaretal.，1975；Tapponnieretal.，1982)；2)下地殼管道流模型，即青藏高原周緣中下地殼廣泛分布塑性通道流(Roydenetal.，1997；Clarketal.，2000)。因此，很多地震學(xué)者的研究都是圍繞上述模型試圖解釋高原地殼增厚、殼幔變形的動(dòng)力學(xué)機(jī)制(Shapiroetal.，2004；李永華等，2006；Yuanetal.，2013；Shenetal.，2015；Wangetal.，2016)。地震波走時(shí)成像是研究地球內(nèi)部殼幔結(jié)構(gòu)和動(dòng)力學(xué)特征的重要手段，走時(shí)成像的分辨率與臺(tái)站空間分布和射線交叉程度正相關(guān)。以往由于臺(tái)站空間尺度大、分布不均勻且稀少，所得研究結(jié)果的精度和分辨率受到了限制。現(xiàn)今，開展密集臺(tái)陣觀測(cè)已成為一種研究殼幔精細(xì)結(jié)構(gòu)的重要方法。“中國科學(xué)臺(tái)陣探測(cè)”2期項(xiàng)目于2013—2016年在南北地震帶北段布設(shè)了676個(gè)密集流動(dòng)臺(tái)站，臺(tái)站之間平均距離為35km，所得觀測(cè)結(jié)果為得到區(qū)域更精細(xì)的三維波速結(jié)構(gòu)提供了良好的數(shù)據(jù)支撐。以這些數(shù)據(jù)為基礎(chǔ)，目前已在體波成像(張風(fēng)雪等，2018)、接收函數(shù)成像(王興臣等，2017；Wangetal.，2017)和面波成像(Lietal.，2017；潘佳鐵等，2017；鐘世軍等，2017；鄭晨等，2018；Wangetal.，2020)等方面取得了很多研究成果。

背景噪聲成像方法自創(chuàng)建以來，在殼幔結(jié)構(gòu)研究方面得到了廣泛應(yīng)用，如全球尺度(Sageretal.，2020)、區(qū)域尺度(Sabraetal.，2005；Shapiroetal.，2005)和小尺度(Brenguieretal.，2007；Mordretetal.，2015)剪切波速度結(jié)構(gòu)研究，通過噪聲互相關(guān)提取體波并將其用于地幔轉(zhuǎn)換帶結(jié)構(gòu)的研究(Polietal.，2012a，b)。傳統(tǒng)的噪聲面波成像以射線走時(shí)為基礎(chǔ)，假設(shè)面波沿大圓路徑傳播。實(shí)際在橫向結(jié)構(gòu)十分不均勻的地區(qū)，面波的傳播路徑存在一定程度的彎曲，且傳統(tǒng)的面波射線走時(shí)成像往往采用正則化反演方法，正則化過程中的模型平滑參數(shù)的選擇容易引入虛假異常。

近年來，Lin等(2009)發(fā)展了基于程函方程的噪聲面波成像方法。程函層析成像方法充分考慮射線的彎曲，通過追蹤波前局部相位變化得到波場(chǎng)走時(shí)梯度，進(jìn)而由波場(chǎng)梯度得到慢度矢量。雖然面波程函層析成像具有反演速度快、分辨率高等優(yōu)點(diǎn)，但要求臺(tái)站布局足夠密集，其分辨率主要由臺(tái)站孔徑和格點(diǎn)尺寸決定。Lin等(2009)提出的噪聲程函成像方法首先需提取頻散曲線，目前廣泛采用的頻時(shí)分析(Frequency-time analysis，F(xiàn)TAN)相速度頻散提取方法(Levshinetal.，1972，1992)是根據(jù)高斯窄帶通濾波后解析信號(hào)的瞬時(shí)相位的定義得到的，其中需要相位模糊度參數(shù)，而相位模糊度參數(shù)由理論或經(jīng)驗(yàn)相速度頻散曲線來確定，無疑將帶來系統(tǒng)誤差(Bensenetal.，2007；Linetal.，2008；Bouéetal.，2013)。Jin等(2015)發(fā)展了基于多通道互相關(guān)相位延遲的天然地震面波自動(dòng)成像方法，通過自動(dòng)計(jì)算鄰近所有臺(tái)站對(duì)的相對(duì)相位走時(shí)延遲，采用程函成像方法反演得到面波頻散的二維分布圖像。多通道互相關(guān)相能夠抑制通過相位延遲方法計(jì)算得到的面波頻散的非相干噪聲，并降低因多路徑散射波導(dǎo)致的測(cè)量偏差。

圖1 研究區(qū)域的構(gòu)造背景與臺(tái)陣分布Fig.1 Tectonic setting and stations in the study area.藍(lán)色三角為科學(xué)臺(tái)陣2期流動(dòng)臺(tái)站，灰色實(shí)線為主要的構(gòu)造邊界或斷裂。HYF 海原斷裂；SQS 南祁連縫合帶；KF 昆侖斷裂；HT 河套-吉蘭泰盆地；ALS 阿拉善盆地；OD 鄂爾多斯地塊；YCD 銀川盆地；QL 祁連造山帶； KQ 西秦嶺造山帶；SG 松潘-甘孜地塊

本研究收集了“中國科學(xué)臺(tái)陣探測(cè)” 2 期項(xiàng)目的676個(gè)流動(dòng)臺(tái)陣記錄的連續(xù)波形數(shù)據(jù)(圖1)，通過噪聲互相關(guān)計(jì)算疊加2015年1—12月的互相關(guān)波形，得到ZZ、ZR、RZ、ZT、TZ、RT、TR、TT等9個(gè)分量的經(jīng)驗(yàn)格林函數(shù)。本文將Jin等(2015)的方法首次發(fā)展應(yīng)用到背景噪聲層析成像中，首先采用多通道互相關(guān)計(jì)算相位的相對(duì)走時(shí)延遲；然后基于程函層析成像得到整個(gè)研究區(qū)域Rayleigh面波8～40s周期的相速度高分辨率分布結(jié)果，大部分研究區(qū)域的分辨率在40～100km之間；最后對(duì)地震面波程函成像結(jié)果和傳統(tǒng)噪聲程函成像結(jié)果進(jìn)行了對(duì)比分析，并結(jié)合其他研究成果探討了研究區(qū)域的構(gòu)造變形和動(dòng)力學(xué)過程。

1 數(shù)據(jù)與方法

1.1 波形互相關(guān)計(jì)算相位延遲

Jin等(2015)提出的多通道互相關(guān)地震面波成像方法，是對(duì)同一震源、不同臺(tái)站接收到的面波波形進(jìn)行互相關(guān)，以自動(dòng)測(cè)量相對(duì)相位和波包群的相對(duì)變化。對(duì)于噪聲經(jīng)驗(yàn)格林函數(shù)，我們以一個(gè)臺(tái)站為“震源”，測(cè)量以這個(gè)臺(tái)站為虛源的所有臺(tái)站對(duì)的相位延遲。假設(shè)以臺(tái)站1為震源，臺(tái)站1與臺(tái)站2、臺(tái)站3之間的經(jīng)驗(yàn)格林函數(shù)波形分別為S1、S2，則S1和S2之間的加窗互相關(guān)函數(shù)可表示為含5個(gè)參數(shù)的小波近似：

F(ω)*WsC(t)≈AGa[σ(t-tg)]cos[ω(t-tp)]

(1)

通過非線性最小二乘擬合可以得到群延遲和相延遲tg和tp。其中，F(xiàn)(ω)為對(duì)應(yīng)頻率ω的高斯窄帶通濾波器，Ga為高斯函數(shù)，A為尺度因子，WsC(t)為加窗互相關(guān)函數(shù)，σ為半帶寬。

1.2 程函層析成像

與傳統(tǒng)的基于射線走時(shí)的成像方法不同，程函層析成像是從基于射線理論的程函方程發(fā)展而來的，標(biāo)準(zhǔn)形式的程函方程為(Shearer，2009)

(2)

(3)

式中，δτp為波場(chǎng)相對(duì)相位延遲時(shí)間，ri為2個(gè)臺(tái)站之間的路徑。

與Lin等(2009)通過求解式(2)中的波場(chǎng)梯度求解慢度空間不同，Jin等(2015)采用最小方差方法反演單個(gè)事件的慢度空間分布，反演的目標(biāo)函數(shù)為

(4)

權(quán)重因子λ可控制反演的平滑度，δτ為相對(duì)相位延遲時(shí)間，可對(duì)觀測(cè)波形與參考波形進(jìn)行互相關(guān)計(jì)算得到。

1.3 數(shù)據(jù)處理和質(zhì)量控制

進(jìn)行單臺(tái)數(shù)據(jù)預(yù)處理時(shí)，首先按照Bensen等(2007)提出的過程進(jìn)行處理，即波形截取、去儀器響應(yīng)、去均值和去趨勢(shì)，并在頻率域進(jìn)行歸一化；然后采用Poli等(2012b)提出的地震窗截?cái)喾椒?，使用長(zhǎng)4h的時(shí)間窗在每天的連續(xù)波形上滑動(dòng)，并進(jìn)行50%的重疊，去除每個(gè)時(shí)間窗內(nèi)大于3倍標(biāo)準(zhǔn)差的波形；最后計(jì)算互相關(guān)函數(shù)并經(jīng)線性疊加得到最終的互相關(guān)函數(shù)。我們通過互相關(guān)函數(shù)的Hilbert變換計(jì)算得到經(jīng)驗(yàn)格林函數(shù)，相較于計(jì)算互相關(guān)函數(shù)的負(fù)導(dǎo)數(shù)得到的經(jīng)驗(yàn)格林函數(shù)，Hilbert變換具有不改變波形振幅的優(yōu)點(diǎn)。為了減少噪聲源非均勻分布對(duì)反演結(jié)果的影響，我們對(duì)因果信號(hào)和非因果信號(hào)進(jìn)行反序平均，得到最終的經(jīng)驗(yàn)格林函數(shù)(Linetal.，2008)。

按照J(rèn)in等(2015)的定義，波形互相關(guān)系數(shù)為

(5)

1.4 相速度計(jì)算

本文采用Jin等(2015)提出的程函層析成像方法(式(3))，對(duì)以單臺(tái)為虛源的所有射線的相對(duì)相延遲進(jìn)行反演，得到網(wǎng)格尺寸為0.2°×0.2°的波場(chǎng)相速度分布。

圖2 為以臺(tái)站15581和51555為虛震源的8s周期Rayleigh面波相速度結(jié)果，其中箭頭為波場(chǎng)走時(shí)梯度方向，即波場(chǎng)的傳播方向。8s周期面波的相速度對(duì)淺層地殼的S波速度較為敏感，受地表構(gòu)造約束，鄂爾多斯地塊西緣、銀川盆地、祁連造山帶西部和松潘-甘孜地塊的東北部表現(xiàn)為低速異常，而西秦嶺造山帶東部為高速異常。

圖2 以臺(tái)站15581(a)和51555(b)為虛源(紫色五角星)的8s周期相速度及波場(chǎng)梯度方向Fig.2 Phase velocity and gradient direction at the 8s period of Rayleigh wave with the central station 15581(a) and 51555(b) as the virtual source.

圖3 臺(tái)站64016對(duì)應(yīng)的所有射線的走時(shí)延遲與距離差的關(guān)系Fig.3 Travel time Delay with relative distance difference for the Station 64016.彩色叉為不同周期的走時(shí)延遲與均值之差<8s的結(jié)果，黑色圓圈為超過8s的結(jié)果

圖4 平均頻散曲線Fig.4 The mean dispersion curve.藍(lán)色實(shí)線為利用頻時(shí)分析方法(FTAN)計(jì)算的平均相速度頻散，紅色實(shí)線為本文結(jié)果；誤差棒為測(cè)量的標(biāo)準(zhǔn)差

最終每個(gè)周期的相速度結(jié)果由所有單臺(tái)結(jié)果加權(quán)疊加得到，主要分為3步：首先對(duì)所有單臺(tái)相速度結(jié)果進(jìn)行加權(quán)平均疊加得到平均相速度結(jié)果；之后對(duì)每個(gè)單臺(tái)相速度進(jìn)行選擇，剔除單臺(tái)走時(shí)超過8s的結(jié)果；最后對(duì)剔除后的結(jié)果進(jìn)行加權(quán)平均疊加，得到最終的相速度分布。

1.5 反演分辨率估計(jì)與誤差分析

程函層析成像反演充分考慮了射線的彎曲，反演的分辨率由臺(tái)陣的分布、射線交叉程度和網(wǎng)格大小決定(Linetal.，2009；Qiuetal.，2019)。本文采用檢測(cè)板測(cè)試方法給出了0.4°×0.4°、1°×1°尺度異常格點(diǎn)的檢測(cè)板恢復(fù)結(jié)果。

圖5 8s、25s和40s周期的檢測(cè)板測(cè)試結(jié)果Fig.5 Checkerboard test image at 8s，25s and 40s periods.a 網(wǎng)格尺寸為0.4°×0.4°的結(jié)果；b 網(wǎng)格尺寸為1°×1°的結(jié)果

為檢驗(yàn)結(jié)果的可靠性，我們還進(jìn)行了誤差分析檢驗(yàn)，每個(gè)格點(diǎn)相速度的誤差由所有單臺(tái)源相速度的標(biāo)準(zhǔn)誤差給出(Linetal.，2009)：

(6)

其中，n為單臺(tái)源的數(shù)目。圖6 為8～40s周期的相速度誤差，除邊緣地區(qū)外，研究區(qū)域的相速度誤差都在40m/s以下，表明反演的精度較高，結(jié)果可靠。

圖6 8～40s周期的相速度誤差Fig.6 Phase velocity uncertainty at periods 8～40s.

2 結(jié)果

圖7 不同周期的相速度深度敏感核Fig.7 Sensitivity kernel of the phase velocity at different periods.

短周期(8～12s)相速度主要反映了上地殼(深6～15km)的S波速度結(jié)構(gòu)(圖8)，上地殼的速度結(jié)構(gòu)與區(qū)域地表地貌特征具有很好的相關(guān)性，鄂爾多斯盆地西緣、銀川地塹、河套-吉蘭泰盆地、阿拉善盆地南部表現(xiàn)為顯著的低速異常，另外祁連造山帶西部、松潘-甘孜地塊東北部也顯示出低速異常(約3.1km/s)。高速異常地區(qū)與造山帶、山脈隆起的分布一致，如西秦嶺造山帶、祁連造山帶東部呈現(xiàn)高速特征。以上結(jié)果與前人的研究一致(Shenetal.，2016；Lietal.，2017；潘佳鐵等，2017；鐘世軍等，2017)。

中等周期(18～22s)的相速度主要對(duì)中地殼(深15～30km)的S波速度結(jié)構(gòu)敏感(圖8)。祁連造山帶西部、松潘-甘孜地塊東北部表現(xiàn)出明顯的低速異常(約3.1km/s)，而河套-吉蘭泰盆地呈現(xiàn)低速異常特征，隨著周期增大低速異常逐漸減弱；阿拉善盆地南部呈現(xiàn)微弱的低速異常；蒙古褶皺系、祁連造山帶東部、鄂爾多斯盆地西緣和西秦嶺造山帶為高速特征。

長(zhǎng)周期(25～40s)相速度主要對(duì)下地殼—上地幔頂部的S波速度結(jié)構(gòu)敏感(圖9)，下地殼速度結(jié)構(gòu)的橫向分布與莫霍面的起伏相關(guān)(Wangetal.，2017；王興臣等，2017)。青藏高原內(nèi)部與其周緣地區(qū)的莫霍面深度具有顯著的差異性，青藏地塊內(nèi)部的祁連造山帶西部和松潘-甘孜地塊東北部表現(xiàn)為顯著的低速異常特征，而青藏地塊周緣的阿拉善盆地、鄂爾多斯盆地西緣、西秦嶺造山帶等大范圍區(qū)域顯示出明顯的整體一致的高速特征，但河套-吉蘭泰盆地與周緣相比顯示出顯著的低速異常。先前的面波成像研究也顯示同樣的結(jié)果(Lietal.，2017；潘佳鐵等，2017；鐘世軍等，2017；楊志高等，2019)。

圖8 8s、12s、18s和22s周期的相速度Fig.8 Phase velocity at 8s，12s，18s and 22s periods.

圖9 25s、30s、35s和40s周期的相速度Fig.9 Phase velocity at 25s，30s，35s and 40s periods.

3 討論

3.1 與他人結(jié)果的對(duì)比

程函層析成像的優(yōu)勢(shì)在于能夠反演局部小尺度的結(jié)構(gòu)特征，具有較高的橫向分辨率。Jin等(2015)創(chuàng)立的多通道互相關(guān)面波自動(dòng)成像方法主要用于長(zhǎng)周期(>20s)天然地震面波成像。鐘世軍等(2017)利用科學(xué)臺(tái)陣2期天然地震面波數(shù)據(jù)，同樣采用多通道互相關(guān)面波自動(dòng)成像方法(Jinetal.，2015)得到了長(zhǎng)周期地震面波的相速度頻散。本文將該方法應(yīng)用到短周期的噪聲面波成像中。我們將本文的相速度純路徑頻散結(jié)果與鐘世軍等(2017)的地震面波成像進(jìn)行了對(duì)比。

圖10 本文結(jié)果(a)與鐘世軍等(2017)結(jié)果(b)的差異Fig.10 The phase velocity difference at different periods between this study (a) and the result from ZHONG Shi-jun et al.，2017 (b).

圖10 為本文與鐘世軍等(2017)結(jié)果的比較(14s、20s、40s周期)，左列(圖10a)為本文的結(jié)果，右列(圖10b)為鐘世軍等(2017)的結(jié)果。從圖10 可以看出，對(duì)于14s、20s周期二者在大部分研究區(qū)域的一致性較好，相速度差異在誤差范圍之內(nèi)，均可刻畫局部小尺度速度結(jié)構(gòu)的細(xì)節(jié)，本文結(jié)果的分辨率為20～40km，特別是短周期(8～25s)的分辨率可達(dá)20～30km。但祁連山西部和西秦嶺的14s周期的結(jié)果具有較大的差異，鐘世軍等(2017)給出的祁連山西部的速度值則更低。由于天然地震的短周期面波射線密度較低，因此頻散的可靠性較低。

在40s周期，兩者在大尺度上特征一致，但在小區(qū)域上的差異較大，天然地震面波成像結(jié)果具有更高的分辨率，推測(cè)主要原因在于兩者使用的數(shù)據(jù)存在差異，背景噪聲源的分布特征決定了長(zhǎng)周期的面波頻散數(shù)目較少，而地震面波長(zhǎng)周期頻散數(shù)目更多，因此結(jié)果具有更高的分辨率。

3.2 與傳統(tǒng)噪聲程函層析成像方法的對(duì)比

Lin等(2009)提出的傳統(tǒng)噪聲程函層析成像方法首先需要提取面波頻散曲線，然后再進(jìn)行程函層析成像。目前，可使用自動(dòng)FTAN方法提取Rayleigh波相速度頻散曲線(Levshinetal.，1972，1992)，利用窄帶通高斯濾波器對(duì)噪聲互相關(guān)函數(shù)的包絡(luò)信號(hào)進(jìn)行濾波，測(cè)量中心頻率對(duì)應(yīng)的最大振幅和到時(shí)，進(jìn)而提取面波群速度頻散，通過相慢度、瞬時(shí)相位、群速度和相位模糊度參數(shù)的關(guān)系式導(dǎo)出相速度。其中，需要由參考模型頻散曲線或已有頻散得到相位模糊度參數(shù)(Bensenetal.，2007；Linetal.，2008)。

圖11 本文方法反演的結(jié)果與Lin等(2009)的傳統(tǒng)程函成像方法結(jié)果的對(duì)比(12s和30s周期)Fig.11 Comparison of phase velocity maps for periods of 12s and 30s.a 本文計(jì)算結(jié)果；b 傳統(tǒng)程函成像結(jié)果；c 二者之差(圖a減去圖b的結(jié)果)

3.3 低速層成因的簡(jiǎn)單解釋

3.3.1 河套-吉蘭泰盆地的低速層

地殼淺部的低速層往往與其上的沉積蓋層密切相關(guān)。研究區(qū)內(nèi)的銀川盆地、鄂爾多斯盆地西緣和河套-吉蘭泰盆地的上地殼存在巨厚沉積層，河套-吉蘭泰盆地的上地殼存在顯著的低速異常并持續(xù)到18s周期。已有的地質(zhì)研究表明，河套-吉蘭泰盆地是新生代以來形成的大型斷陷盆地，新生界最大沉積層可厚達(dá)10km以上(鄭孟林等，2006)。而長(zhǎng)周期(35～40s)相速度(相當(dāng)于40～80km深度)結(jié)果表明河套-吉蘭泰盆地的下地殼、上地幔也存在弱低速層。Chen等(2009)的S波接收函數(shù)研究和He等(2017)的面波成像研究表明河套-吉蘭泰盆地的巖石圈厚度只有80km，遠(yuǎn)低于其南部的鄂爾多斯盆地的巖石圈厚度(160km)。同時(shí)，楊嵩等(2013)對(duì)華北地區(qū)的巖石圈厚度及上地幔溫度進(jìn)行了研究，發(fā)現(xiàn)河套-吉蘭泰盆地下方的溫度明顯高于其周緣地區(qū)。本文推測(cè)，“大地幔楔”中的地幔軟流圈熱物質(zhì)上涌可能是下地殼及巖石圈地幔出現(xiàn)弱低速異常的主要原因(鄭晨等，2018；Leietal.，2019)。

3.3.2 松潘-甘孜地塊和祁連山西部的低速層

有研究表明，地殼的徑向各向異性是云母晶體的定向排列所致。正的徑向各向異性與擴(kuò)張型地殼(即水平向物質(zhì)流)對(duì)應(yīng)；反之，負(fù)的徑向各向異性對(duì)應(yīng)垂向的增厚。而松潘-甘孜地塊東北部存在弱的負(fù)徑向各向異性特征(汪洋等，2001；Baoetal.，2015；楊志高等，2019)，表明不太可能存在水平向的下地殼管道流。

4 結(jié)論

本文基于青藏高原東北緣及鄰區(qū)密集臺(tái)陣記錄的連續(xù)波形得到了背景噪聲經(jīng)驗(yàn)格林函數(shù)，進(jìn)而發(fā)展了基于多通道互相關(guān)的面波程函層析成像方法，首次將其應(yīng)用于背景噪聲層析成像中，反演得到了研究區(qū)域不同周期的相速度分布圖像。

經(jīng)對(duì)比，本文所得結(jié)果與已有的天然地震程函層析成像結(jié)果一致，表明本文結(jié)果是可靠的，且具有高分辨率和高精度的特點(diǎn)。與傳統(tǒng)程函層析成像方法的研究結(jié)果對(duì)比，本文結(jié)果降低了長(zhǎng)周期信號(hào)由于信噪比的降低和有效射線數(shù)目減少所帶來的誤差，壓制了非相干多路徑散射面波引起的測(cè)量誤差，提高了反演的精度和穩(wěn)定性。

反演得到的8～40s周期的相速度分布圖像表明，河套-吉蘭泰盆地的下地殼、上地幔顯示出弱低速特征，地幔高溫、軟弱物質(zhì)的侵入可能造成了該區(qū)域下地殼、上地幔出現(xiàn)弱低速。松潘-甘孜地塊東北部中下地殼存在低速體異常，可能與下地殼的部分熔融有關(guān)。祁連造山帶西段中下地殼廣泛存在低速層，分析認(rèn)為，祁連山西段中下地殼可能不具備部分熔融或地殼管道流，低速層可能與其地殼在外力作用下的縮短、增厚和解耦有關(guān)。

致謝中國地震局地球物理研究所“中國地震科學(xué)臺(tái)陣探測(cè)中心”為本文提供了連續(xù)地震波形數(shù)據(jù)；Jin 和Gaherty博士為本研究提供了面波自動(dòng)成像程序(ASWMS)(1)http：∥www.iris.edu/ds/products/aswms/。，馮力理博士為本研究提供了傳統(tǒng)程函層析成像程序(2)https：∥github.com/NoisyLeon。；審稿專家為本文提供了很好的意見和建議。在此一并表示感謝！