多尺度陣列嵌套組合反演賓川氣槍源區橫波速度結構

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

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

0 引言

背景噪聲成像是利用背景噪聲分析提取面波頻散信息，進而反演地下介質速度結構的方法，是近年來地震學研究的熱點，無論在工程物理勘探還是地殼上地幔大尺度速度結構研究方面都得到了快速的發展，其中面波頻散曲線的提取是這些研究的關鍵問題.基于不同的研究尺度，頻散提取方法可以分為兩種：(1)淺地表面波勘探尺度的空間自相關法(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)等.隨著面波成像技術的發展，高階面波在地下結構分析中的影響越來越被重視，利用多階面波聯合反演可以降低反演結果的非唯一性，使反演結果更接近真實結構(羅銀河, 2008；Xia et al., 2003；何耀峰, 2005)，而傳統方法對高階面波的反應不敏感，很難從實際背景噪聲中提取到清晰的高階面波信息.最近陳曉非課題組提出了基于臺陣數據的頻率-貝塞爾變換法(Frequency-Bessel, F-J方法)，能夠有效地從背景噪聲中提取分辨率高、頻率范圍廣的基階和高階面波頻散曲線(Wang et al., 2019；Pan et al., 2019；Wu et al., 2020)，無論在淺地表勘探(李雪燕，2019；吳華禮，2019)，還是地殼、上地幔等大尺度的速度結構成像研究中都取得了很好的應用成果(王建楠，2019；Zhan et al., 2020; Wu et al., 2020)，對被動源和主動源數據研究都有不錯的效果(Li and Chen, 2020)，是一種應用前景非常廣闊的新方法.

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

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

1 臺站數據

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

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

圖1 多尺度陣列臺站分布(a) 云南地區-YN； (b) 滇西地區-BC； (c) 2017年密集臺陣區-BA，其中紅色矩形框內為本文研究的最小區域-PA，臺站用三角表示，黃色表示云南區域地震臺網，綠色表示賓川氣槍源臺網，藍色為2017年密集臺陣，紅色為研究小區域臺站，黑色五角星為氣槍發射臺位置.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 數據處理

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

圖2 數據處理流程圖Fig.2 Flow chart of data processing

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

(1)

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

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

(2)

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

(3)

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

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

3 結果及分析

3.1 F-J頻散譜

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

圖3 多尺度陣列的背景噪聲數據提取到的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

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

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

3.2 橫波速度反演

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

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

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

圖7 擾動區間內50個初始速度模型的分布Fig.7 Initial velocity models within the disturbance interval

圖8 分別利用PA區基階和基階加一階的頻散信息反演得到的擬合頻散曲線和橫波速度模型(a)(c) 利用基階頻散曲線得到的反演結果; (b)(d) 基階加一階頻散曲線的反演結果；橫波速度模型圖(c)(d)中，兩側虛線為以參考模型為中心的初始模型擾動范圍區間，灰色實線表示不同初始模型反演得到的結果分布，顏色越深表示該速度值比重越大，紅色實線為最后反演得到的平均速度模型，下同.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區頻散信息反演得到淺層橫波速度后，逐漸擴大區域尺度繼續反演(圖9).首先擴大到BA區，獲取PA-BA組合后的頻散曲線，可以發現基階頻段拓寬，低頻段從0.55 Hz延伸到0.18 Hz，50個模型的頻散曲線擬合較好(圖9a)，橫波速度的反演深度加深，最深達到12 km，深度10 km處出現明顯的橫波速度增加現象，12 km處稍微減小(圖9d).

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

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

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

3.3 結果對比分析

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

圖10 參考模型和反演模型的橫波速度結構(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 結論

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

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

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

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

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