多尺度陣列嵌套組合反演賓川氣槍源區(qū)橫波速度結(jié)構(gòu)

孫楠, 潘磊, 王偉濤, 葉泵, 王彬, 陳曉非,5*

1 中國科學(xué)技術(shù)大學(xué)地球與空間科學(xué)學(xué)院， 合肥 230026 2 云南省地震局， 昆明 650224 3 南方科技大學(xué)地球與空間科學(xué)系 深圳 518055 4 中國地震局地球物理研究所， 北京 100081 5 深圳市深遠海油氣勘探技術(shù)重點實驗室， 深圳 518055

0 引言

背景噪聲成像是利用背景噪聲分析提取面波頻散信息，進而反演地下介質(zhì)速度結(jié)構(gòu)的方法，是近年來地震學(xué)研究的熱點，無論在工程物理勘探還是地殼上地幔大尺度速度結(jié)構(gòu)研究方面都得到了快速的發(fā)展，其中面波頻散曲線的提取是這些研究的關(guān)鍵問題.基于不同的研究尺度，頻散提取方法可以分為兩種：(1)淺地表面波勘探尺度的空間自相關(guān)法(Aki, 1957)、F-K法(Capon, 1969；周云騰和張致付，2020)、相移法(Park et al., 1998；伍敦仕等，2017)、傾斜疊加法(Xia et al., 2007)以及拉東變換法(McMechan and Yedlin, 1981; Luo et al., 2008)等;(2)大地球物理尺度的雙臺法(Knopoff, 1964; Yao et al., 2006)和基于程函方程的相速度測量方法(Lin et al., 2009)等.隨著面波成像技術(shù)的發(fā)展，高階面波在地下結(jié)構(gòu)分析中的影響越來越被重視，利用多階面波聯(lián)合反演可以降低反演結(jié)果的非唯一性，使反演結(jié)果更接近真實結(jié)構(gòu)(羅銀河, 2008；Xia et al., 2003；何耀峰, 2005)，而傳統(tǒng)方法對高階面波的反應(yīng)不敏感，很難從實際背景噪聲中提取到清晰的高階面波信息.最近陳曉非課題組提出了基于臺陣數(shù)據(jù)的頻率-貝塞爾變換法(Frequency-Bessel, F-J方法)，能夠有效地從背景噪聲中提取分辨率高、頻率范圍廣的基階和高階面波頻散曲線(Wang et al., 2019；Pan et al., 2019；Wu et al., 2020)，無論在淺地表勘探(李雪燕，2019；吳華禮，2019)，還是地殼、上地幔等大尺度的速度結(jié)構(gòu)成像研究中都取得了很好的應(yīng)用成果(王建楠，2019；Zhan et al., 2020; Wu et al., 2020)，對被動源和主動源數(shù)據(jù)研究都有不錯的效果(Li and Chen, 2020)，是一種應(yīng)用前景非常廣闊的新方法.

但是，對于橫向結(jié)構(gòu)變化非常大的區(qū)域，目前背景噪聲成像方法存在一些矛盾.即小尺度陣列對于局部淺層結(jié)構(gòu)的反演成像有效，能夠得到分辨率較高的淺層速度結(jié)構(gòu)，但對深部速度結(jié)構(gòu)不敏感.而大尺度陣列對深部速度結(jié)構(gòu)研究有效，對淺層結(jié)構(gòu)卻約束不夠.因此，為得到分辨率較高地從淺層到深層的速度結(jié)構(gòu)，本文中我們采用了一種多尺度陣列嵌套組合的方式提取Rayleigh面波頻散曲線.即用小尺度陣列資料提取相對高頻部分的頻散曲線，而用大尺度陣列資料提取較低頻率的頻散曲線，融合不同尺度陣列資料的分析結(jié)果，獲得寬頻帶頻散曲線，通過對寬頻帶頻散曲線的反演獲得研究目標區(qū)域從淺至深(0～70 km)的橫波速度結(jié)構(gòu).

云南位于青藏高原東南緣，受板塊推擠作用，區(qū)域內(nèi)橫向結(jié)構(gòu)復(fù)雜，中強地震持續(xù)活躍，地下構(gòu)造特征備受關(guān)注.為探測和監(jiān)測區(qū)域地下結(jié)構(gòu)變化，2012年在云南賓川建成了水庫激發(fā)氣槍震源發(fā)射臺，利用氣槍震源可重復(fù)性、頻率低、傳播距離遠的優(yōu)勢進行了一系列研究(陳颙等，2007；王彬等，2016；王寶善等，2011；陳蒙，2014；孫楠和孫耀充，2020；王偉濤等，2017；向涯等，2021).為得到氣槍源發(fā)射臺周邊區(qū)域的淺層速度結(jié)構(gòu)，2017年3—5月圍繞氣槍發(fā)射臺開展了密集臺陣實驗，記錄到大量的連續(xù)數(shù)據(jù)，本研究利用這一密集臺陣記錄到的連續(xù)背景噪聲數(shù)據(jù)，采用F-J方法提取Rayleigh波頻散曲線.由于該地區(qū)近地表橫向結(jié)構(gòu)變化大、地質(zhì)地貌格局復(fù)雜，直接獲得較寬頻帶頻散曲線的難度較大，因此，我們采用從密集臺陣到云南地區(qū)多尺度臺陣嵌套組合的方式，來拓寬頻散曲線的頻段，進而反演得到研究區(qū)域不同深度的橫波速度結(jié)構(gòu)，為該區(qū)地下結(jié)構(gòu)的探測和監(jiān)測提供基礎(chǔ).

1 臺站數(shù)據(jù)

2017年密集臺陣實驗圍繞賓川氣槍震源發(fā)射臺布設(shè)，布設(shè)時間長度從兩周到兩個月不等，記錄了大量的連續(xù)波形數(shù)據(jù)，臺陣覆蓋面積40 km×40 km，臺間距為2～3 km，共使用了381個短周期地震計，頻帶范圍為5 s～150 Hz，覆蓋區(qū)內(nèi)涵蓋盆地和丘陵，地表起伏最大1.4 km，下文中該區(qū)用BA來表示.本文要反演的淺層區(qū)域為BA區(qū)的一部分，位于密集臺陣內(nèi)部氣槍發(fā)射臺的西南側(cè)，覆蓋區(qū)域長25 km寬10 km，包括57個臺站，下文中該區(qū)域用PA來表示.

為反演更深層的信息，我們將研究區(qū)域不斷擴大，收集了賓川氣槍源臺網(wǎng)記錄的連續(xù)波形數(shù)據(jù)，臺網(wǎng)覆蓋滇西地區(qū)(24.25°N—27.00°N，98.70°E—102.05°E)，共包括60個臺站，臺間距最大能到200 km以上，數(shù)據(jù)記錄時間范圍為2015—2018年，下文中該區(qū)用BC來表示.之后，我們將區(qū)域繼續(xù)擴大到云南地區(qū)(20°N—30°N，96°E—107°E)，收集了云南區(qū)域臺網(wǎng)2016—2019年3月共72臺的連續(xù)波形數(shù)據(jù)，該區(qū)內(nèi)臺間距最大為927 km，下文中該區(qū)用YN來表示.本文所有區(qū)域的臺站數(shù)據(jù)均選取了垂直分量的連續(xù)記錄，采用多尺度陣列嵌套組合的方式進行分析計算，具體臺站分布見圖1，4個陣列具體信息見表1.

圖1 多尺度陣列臺站分布(a) 云南地區(qū)-YN； (b) 滇西地區(qū)-BC； (c) 2017年密集臺陣區(qū)-BA，其中紅色矩形框內(nèi)為本文研究的最小區(qū)域-PA，臺站用三角表示，黃色表示云南區(qū)域地震臺網(wǎng)，綠色表示賓川氣槍源臺網(wǎng)，藍色為2017年密集臺陣，紅色為研究小區(qū)域臺站，黑色五角星為氣槍發(fā)射臺位置.Fig.1 Maps showing station distribution of multi-scale arrays(a) Yunnan area-YN； (b) Western Yunnan area-BC； (c) Dense array area in 2017-BA, where the red rectangle is the smallest study area in this paper-PA. Triangles represent stations, yellow for Yunnan regional seismic network, green for Binchuan airgun source seismic network, blue for dense array station in 2017, red for dense array of the smallest study area in this paper, and black pentagram indicates location of airgun source launcher.

表1 多尺度陣列信息表Tab1e 1 Information of multi-scale arrays

2 數(shù)據(jù)處理

文中數(shù)據(jù)處理主要包括4步：連續(xù)數(shù)據(jù)預(yù)處理、背景噪聲互相關(guān)、提取頻散信息和反演橫波速度，具體流程見圖2，對不同陣列數(shù)據(jù)分別作數(shù)據(jù)預(yù)處理，去除儀器響應(yīng)的影響，經(jīng)過多次篩選比較，YN區(qū)和BC區(qū)連續(xù)波形降采樣到10 Hz，密集臺陣降采樣到100 Hz.背景噪聲互相關(guān)計算采用姚華建提供的程序包(Yao et al., 2006, 2011)進行處理.針對連續(xù)波形記錄中地震事件和氣槍信號的影響，采用時間域one-bit的歸一化操作進行消除(Bensen et al., 2007).

圖2 數(shù)據(jù)處理流程圖Fig.2 Flow chart of data processing

頻散曲線提取方法采用陳曉非課題組(Wang et al., 2019)提出的F-J方法.該方法首先對噪聲互相關(guān)函數(shù)的頻譜C(r,ω)作貝塞爾變換，得到互相關(guān)系數(shù)的F-J頻散譜：

(1)

其中，J0是零階第一類貝塞爾函數(shù).然后利用貝塞爾函數(shù)的同階正交性，得到：

I(ω,k)=A·Im[g(ω,k)]，

(2)

其中，g(ω,k)是核函數(shù)，可以寫為：

(3)

D(ω,k)具有頻散特征，頻散圖上的點接近頻散曲線時，D(ω,k)趨近于零，I(ω,k)會達到極值點，這些極值點就是頻散點，而I(ω,k)可以利用背景噪聲的互相關(guān)函數(shù)近似計算出來，這樣就可以通過拾取頻散譜圖上的極值點來獲得頻散曲線信息.

提取頻散曲線后，我們采用基于BFGS校正的擬牛頓法 (何耀峰，2005；田玥和陳曉非，2006；Pan et al., 2019)進行反演.首先以參考模型為中心，參考前人的研究經(jīng)驗，在±0.4 km·s-1擾動范圍內(nèi)隨機生成50個初始速度模型(吳高雄，2020)，利用擬牛頓迭代進行反演，得到50個速度模型結(jié)果，然后將50個反演結(jié)果利用殘差作為權(quán)重參數(shù)進行加權(quán)平均，計算得到最終的橫波速度模型.

3 結(jié)果及分析

3.1 F-J頻散譜

經(jīng)過上述數(shù)據(jù)處理，得到4個尺度區(qū)域的F-J頻散譜圖(圖3)，顯示臺陣區(qū)域越大，頻散譜分布頻段越低.在低頻段，相同頻率下的YN頻譜分布明顯比BC更窄，極值點的拾取也會更準確，PA中低頻譜段噪聲過大也暫不考慮，在高頻段出現(xiàn)多階信號，我們速度模型正演與頻散譜對比，大致確定PA中基階和1階頻散曲線的位置.基于以上考慮，分別拾取F-J頻散譜上的極大值點，連線得到4個區(qū)域的面波頻散曲線.圖4顯示YN區(qū)頻散曲線頻率范圍為0.008～0.15 Hz，BC區(qū)頻散曲線頻段為0.046～0.33 Hz，BA區(qū)頻段為0.18～0.56 Hz，PA區(qū)基階頻段為0.55～1.13 Hz，1階頻段為0.86～1.1 Hz.

圖3 多尺度陣列的背景噪聲數(shù)據(jù)提取到的F-J頻散譜Fig.3 F-J dispersion spectra of multi-scale arrays from ambient noise data

圖4 多尺度陣列的F-J頻散曲線Fig.4 F-J dispersion curves of multi-scale arrays

將四個區(qū)域的頻散譜按照區(qū)域大小順序進行組合(圖5a)，發(fā)現(xiàn)隨著區(qū)域的逐漸擴大(PA→BA→BC→YN)，多尺度陣列的頻散譜圖能夠很好的連接在一起，頻散譜能量逐漸向低頻端拓寬.圖5b為組合后得到的頻散曲線，更直觀地表現(xiàn)出多尺度陣列的頻散曲線在相同頻段能夠很好地吻合，表明越往深部區(qū)域介質(zhì)橫向均勻尺度越大，可以通過擴大區(qū)域尺度來拓寬低頻信息，從而反演深層速度結(jié)構(gòu).

圖5 多尺度陣列嵌套組合的F-J頻散譜(a)和頻散曲線(b)Fig.5 F-J dispersion spectrum (a) and dispersion curves (b) of combined multi-scale arrays

3.2 橫波速度反演

背景噪聲互相關(guān)結(jié)果經(jīng)過F-J變換后，提取到多尺度陣列的頻散曲線，下面將利用基于BFGS校正的擬牛頓迭代法反演不同深度的橫波速度結(jié)構(gòu).反演中對于參考模型地選取，我們綜合了中國地震局地球物理研究所吳建平研究員提供的速度模型和陳思文(2015)利用走時成像研究得到的滇西地區(qū)速度結(jié)構(gòu)模型(圖6)，參考模型共分42層，前兩層層厚1 km，其余層厚2 km，最深為80 km.

圖6 參考速度模型Fig.6 Reference velocity models

首先利用PA區(qū)頻散曲線信息反演該區(qū)的淺層橫波速度.圖7為以參考模型為中心，在±0.4 km·s-1擾動范圍內(nèi)隨機生成的50個初始速度模型分布.當只用基階頻散曲線反演的時候，50個模型的反演結(jié)果(圖8a)顯示：基階頻散曲線(紅色實線)與實際值(黑點)擬合很好，而一階頻散曲線擬合不好，紅色實線分布離散，圖8c為其反演得到的橫波速度結(jié)構(gòu)，灰色代表50個模型的反演結(jié)果分布，顏色越深表示該速度值的比重越大，紅色實線為50個反演結(jié)果加權(quán)平均得到的最終速度模型，結(jié)果顯示橫波速度的反演深度到4 km，0～2 km的橫波速度變化特別大，從2.33 km·s-1增加到3.15 km·s-1(圖8c)，羅睿潔等(2015)利用衛(wèi)星遙感和地表調(diào)查認為該區(qū)域地層主要為泥盆系灰?guī)r，那么橫波速度0.82 km·s-1的變化梯度可能是從沉積層到基底層的轉(zhuǎn)變.反演加入一階頻散曲線后(圖8b)，50個模型的基階和一階頻散曲線擬均能得到很好地擬合，反演模型結(jié)果(圖8d)顯示：50個速度模型在4 km以上更加收斂，也就是高階面波的加入極大地降低反演結(jié)果的非唯一性，提高反演的穩(wěn)定性，多個不同初始模型的反演結(jié)果更集中，減少了反演過程中淺層結(jié)果對初始模型的依賴，而且相對基階，在同一頻率下高階面波的敏感深度更深，使得PA區(qū)橫波速度的反演深度在8 km左右依然很收斂，表明該頻段多階聯(lián)合反演的約束深度能達到8 km處.接下來的反演中，該區(qū)域淺層(<8 km)速度結(jié)構(gòu)就采用PA區(qū)多階反演結(jié)果，固定淺層反演結(jié)果不變，通過不斷擴大臺站分布區(qū)域來反演得到8 km以下的橫波速度結(jié)構(gòu).

圖7 擾動區(qū)間內(nèi)50個初始速度模型的分布Fig.7 Initial velocity models within the disturbance interval

圖8 分別利用PA區(qū)基階和基階加一階的頻散信息反演得到的擬合頻散曲線和橫波速度模型(a)(c) 利用基階頻散曲線得到的反演結(jié)果; (b)(d) 基階加一階頻散曲線的反演結(jié)果；橫波速度模型圖(c)(d)中，兩側(cè)虛線為以參考模型為中心的初始模型擾動范圍區(qū)間，灰色實線表示不同初始模型反演得到的結(jié)果分布，顏色越深表示該速度值比重越大，紅色實線為最后反演得到的平均速度模型，下同.Fig.8 Fitting dispersion curves and shear wave velocity models inverted with PA region model(a) and (c) Using fundamental dispersion curves； (b) and (d) Fundamental plus first-order dispersion curves. In (c) and (d), dotted lines represent disturbance interval of initial models centered at the reference model, grey solid lines represent inversion results using different initial models, the darker the color, the greater the proportion of velocity value; and red solid lines represent the final average velocity model of inversion, the same below.

利用PA區(qū)頻散信息反演得到淺層橫波速度后，逐漸擴大區(qū)域尺度繼續(xù)反演(圖9).首先擴大到BA區(qū)，獲取PA-BA組合后的頻散曲線，可以發(fā)現(xiàn)基階頻段拓寬，低頻段從0.55 Hz延伸到0.18 Hz，50個模型的頻散曲線擬合較好(圖9a)，橫波速度的反演深度加深，最深達到12 km，深度10 km處出現(xiàn)明顯的橫波速度增加現(xiàn)象，12 km處稍微減小(圖9d).

繼續(xù)擴大區(qū)域到滇西地區(qū)，或取BC-BA-PA組合后的頻散曲線進行反演，基階低頻段拓寬到0.046 Hz(圖9b)，反演模型深度顯著加深，最深能達到45 km左右(圖9e)，前人的研究表明(陳思文，2015；白志明和王椿鏞，2004)賓川地區(qū)莫霍面埋深約為45 km，也就是通過拓寬低頻段信息，能反演得到該區(qū)從地表到地殼的橫波速度結(jié)構(gòu).在深度10 km左右出現(xiàn)明顯的高速-低速轉(zhuǎn)變，這與張云鵬等(2020)利用天然地震研究該區(qū)上地殼速度結(jié)構(gòu)的結(jié)果一致，該區(qū)域?qū)?yīng)張云鵬等(2020)研究中x=-10，y=-5的區(qū)域，根據(jù)其計算結(jié)果該區(qū)域在10 km上下也存在明顯的高低速變化現(xiàn)象.

最后將區(qū)域擴大到云南地區(qū)，獲取YN-BC-BA-PA組合的頻散曲線進行反演，顯示低頻最低能達到0.008 Hz(圖9c)，反演深度能達70 km，70 km以上50個不同初始模型反演結(jié)果相對集中，10～45 km的反演結(jié)果比BC-BA-PA的反演結(jié)果更收斂(圖9f)，表明低頻信息的拓寬不僅能加深反演深度，對反演結(jié)果也具有一定的約束，降低了反演結(jié)果的非唯一性，在45 km左右存在橫波速度迅速增加的現(xiàn)象，與該區(qū)域莫霍面的深度吻合.

圖9 多尺度區(qū)域反演得到的頻散曲線擬合圖和橫波速度模型圖Fig.9 Fitting dispersion curves and shear wave velocity models inverted by multi-scale arrays

3.3 結(jié)果對比分析

將本文最終反演得到橫波速度模型與參考模型進行對比，發(fā)現(xiàn)二者的深部速度結(jié)果基本吻合，而淺層結(jié)果存在明顯差異，尤其是近地表處(圖10a).將兩種速度結(jié)構(gòu)分別進行正演計算，與區(qū)域F-J頻散譜分布進行比較，圖10b中黃色虛線為參考模型的正演結(jié)果，紅色虛線為反演模型的正演結(jié)果，顯示反演模型的正演結(jié)果(紅色虛線)與實際頻散譜在所有頻段都吻合.而參考模型的正演結(jié)果(黃色虛線)與實際頻散曲線在低頻段(<0.08 Hz，由YN區(qū)提取)基本重合，兩者反映的深部的速度結(jié)構(gòu)應(yīng)該一致，圖10a速度模型對比中證實了這一結(jié)果，二者在深部的速度結(jié)構(gòu)基本一致.而0.08 Hz以上參考模型的正演結(jié)果(黃色虛線)與實際頻散譜分布有明顯區(qū)別，對應(yīng)二者的速度結(jié)構(gòu)在淺層存在明顯差異.那么，對于該區(qū)域不同深度的橫波速度結(jié)構(gòu)，這種多尺度陣列嵌套組合的反演結(jié)果相對更準確.

圖10 參考模型和反演模型的橫波速度結(jié)構(gòu)(a)與F-J頻散譜正演模擬分布(b)Fig.10 Shear wave velocity (a) and F-J dispersion spectrum forward modeling (b) of reference model and inversion model

4 結(jié)論

本文基于多尺度陣列嵌套組合的方式，利用頻率-貝塞爾變換方法，提取連續(xù)背景噪聲數(shù)據(jù)中面波頻散信息，將不同尺度的頻散曲線融合進而反演，得到了賓川氣槍發(fā)射臺區(qū)域不同深度的橫波速度結(jié)構(gòu).得到以下主要結(jié)論：

淺層橫波速度反演結(jié)果顯示，多階面波反演結(jié)果更加收斂，反演深度更深.同樣，利用多尺度陣列嵌套組合的方式(小區(qū)域PA—密集臺陣BA—賓川BC—云南YN)，使低頻信息不斷拓寬(0.55 Hz—0.18 Hz—0.046 Hz—0.008 Hz)，對反演結(jié)果也能提供一定的約束，降低反演結(jié)果非唯一性.

從該區(qū)域反演得到0～70 km的橫波速度變化特征可以看出：速度隨深度曲折增長，淺層0～2 km處迅速增加，可能表明該區(qū)域沉積層的厚度在2 km以上的范圍內(nèi).速度在深度10 km附近存在高低速的轉(zhuǎn)變，說明該區(qū)域地殼內(nèi)部存在速度異常體.在深度45 km處橫波速度出現(xiàn)明顯增大，與前人研究所得的該區(qū)域莫霍面的埋深相吻合，從結(jié)果上證實了本文利用多尺度陣列嵌套組合的方法是可行的.將反演得到的速度模型正演結(jié)果與實際頻散譜進行對比，發(fā)現(xiàn)反演得到的區(qū)域橫波速度的結(jié)果相對更準確，更能反映地下的真實介質(zhì)情況，為該區(qū)地下結(jié)構(gòu)的探測和監(jiān)測提供了基礎(chǔ).

本文提出的多尺度陣列嵌套組合方式，為背景噪聲成像提供了一種新的思路和方法.通過將多個陣列提取的頻散曲線組合，得到頻帶更寬的頻散曲線，既可以反演深度范圍更廣的速度結(jié)構(gòu)，又可以對反演結(jié)果提供更好地約束，而且該方法操作簡單方便，實用性極強，在區(qū)域結(jié)構(gòu)研究中具有非常好的應(yīng)用前景.

致謝感謝審稿專家為文章的修改提供的寶貴建議，感謝中國地震局地球物理研究所張云鵬博士、云南地震局彭關(guān)靈工程師和王光明工程師為本文提供資料.