背景噪聲互相關(guān)函數(shù)多尺度成像研究

黎昌成，陳曉非2,,3*

1 南方科技大學(xué)地球與空間科學(xué)系，廣東深圳 518055 2 深圳市深遠(yuǎn)海油氣勘探技術(shù)重點實驗室(南方科技大學(xué)), 廣東深圳 518055 3 南方海洋科學(xué)與工程廣東省實驗室(廣州)，廣州 511458

0 引言

背景噪聲成像方法已被廣泛用于不同尺度地下介質(zhì)結(jié)構(gòu)成像中(Shapiro et al.,2005; Yao et al.,2006；Xia et al.,1999,2003;Harmon et al.,2007;Shen and Ritzwoller,2016;Pan et al.,2019).通常我們會對背景噪聲進行處理，然后獲得有效信息:如走時(Poli et al., 2012a; Mroczek and Tilmann, 2021)，振幅(Knapmeyer-Endrun et al., 2016,2017)，波形(Tromp et al., 2010; Sager et al., 2018a, 2018b)，頻散曲線(Shapiro et al., 2005; Wang et al.,2019; Xi et al., 2021)等.傳統(tǒng)基于頻散曲線成像方法，通常從兩個臺站間的互相關(guān)函數(shù)中提取基階頻散曲線來進行直接計算(Shapiro et al., 2005;Yao et al.,2006,2008).近年來，為了減少反演多解性，從互相關(guān)函數(shù)中提取高階頻散曲線來反演地下介質(zhì)結(jié)構(gòu)的方法受到越來越多的關(guān)注(Xia et al.,2003; Ikeda et al.,2013; Wang et al.,2019; Pan et al.,2019; Wu et al.,2020).但頻散曲線方法主要基于層狀假設(shè)，推廣到二維三維復(fù)雜介質(zhì)成像有降低空間分辨率風(fēng)險.因此，近年來直接利用互相關(guān)函數(shù)來進行成像的方法也得到了發(fā)展，基于背景噪聲互相關(guān)函數(shù)的反演方法已在全球尺度結(jié)構(gòu)成像中證實可以提高反演模型分辨率(Tromp et al., 2010; Sager et al., 2019)，但是該方法依賴于噪聲源的全球分布，導(dǎo)致無法推廣到局部區(qū)域尺度地下介質(zhì)結(jié)構(gòu)成像上(Sager et al.,2019).因此，本文提出了不依賴于噪聲源的互相關(guān)函數(shù)成像方法.另一方面，不同多尺度類反演方法和策略被廣泛應(yīng)用于不同地震成像領(lǐng)域(Gorbatov and Kennett,2003; Lebedev and van der Hilst,2008;Tian et al.,2009；Fishwick et al., 2005; Zhou et al., 2006; Nettles and Dziewoński,2008；Fichtner et al.,2013; 徐義賢和王家映,1998；潘紀(jì)順等,2009; 張文生等,2015;孫楠等,2021)；為了減小計算量增加反演穩(wěn)定性，Bunks等(1995)基于模型重采樣方法，對大尺度模型進行了反演.在地震勘探領(lǐng)域，多尺度反演主要集中于用低頻反演得到大尺度初始模型，然后直接加入高頻信號，恢復(fù)更加精確的結(jié)構(gòu)(Virieux and Operto,2009;Fu et al.,2020;朱寶衡和謝春麗,2021)，但這種策略受網(wǎng)格大小限制并不能從模型尺度上對分辨率進行提高.雖然可以通過劃分更加精細(xì)的網(wǎng)格來提高分辨率，但精細(xì)網(wǎng)格對計算量、反演穩(wěn)定性等方面均會有影響.因此，需要提出新的多尺度反演方法權(quán)衡反演深度、穩(wěn)定性、以及計算量.在背景噪聲成像中，實際數(shù)據(jù)采集，信號處理(從背景噪聲數(shù)據(jù)中提取有效信號)，正演，反演是其重要組成部分.因此，本文將系統(tǒng)的從互相關(guān)函數(shù)處理以及正反演三個方面對互相關(guān)函數(shù)多尺度成像方法進行研究，并提出尺度隨頻率變化的多尺度反演方法和策略，為獲得更深大尺度和淺部精細(xì)結(jié)構(gòu)提供新的方法和策略.

1 理論方法

在背景噪聲成像理論方法中，數(shù)據(jù)處理、正反演方法均是成像的重要組成部分.

對于背景噪聲信號處理，如何獲得穩(wěn)定的噪聲互相關(guān)數(shù)據(jù)，并且衡量實際數(shù)據(jù)誤差的不確定性對于反演十分重要.前人基于各種有效方法從背景噪聲中獲得了穩(wěn)定的互相關(guān)函數(shù)(Benson et al.,2007；Weaver and Yoritomo 2018；Li et al.,2020；Xie et al.,2020；Yang et al.,2020).為了確定互相關(guān)函數(shù)中不同時刻有效信號權(quán)重，本文采用如下方法對噪聲互相關(guān)進行處理(Aster et al.,2012).

對于不同時間段獲得的互相關(guān)函數(shù)將其寫成矩陣形式：

d=[d1,d2,d3,…,dn]T,

(1)

X=[d1-dm,d2-dm,d3-dm,…,dn-dm]T,

(2)

協(xié)方差矩陣為

C=diag(XTX)，

(3)

任意選取若干段實際噪聲數(shù)據(jù)，對不同時間段背景噪聲數(shù)據(jù)處理得到的互相關(guān)函數(shù)，對其求平均可得反演所需數(shù)據(jù)和數(shù)據(jù)協(xié)方差矩陣對角線元素.結(jié)果分別如圖1和圖2所示，為避免不必要噪聲誤差，數(shù)據(jù)處理時，對確定噪聲誤差段添加了零權(quán)重因子.

圖1 不同時間背景噪聲合成的互相關(guān)函數(shù)實部信息

圖2 由不同時間互相關(guān)函數(shù)所求得的協(xié)方差矩陣對角線元素

除背景噪聲信號處理外，本文方法成立的前提條件之一便是正演方法的發(fā)展與成熟.因此，首先建立地下模型參數(shù)與噪聲互相關(guān)函數(shù)之間的關(guān)系式，然后求解該方程.前人給出了背景噪聲互相關(guān)函數(shù)與格林函數(shù)之間的關(guān)系(Sánchez-Sesma and Campillo，2006; Sato et al., 2012; Haney et al., 2012; Haney and Nakahara,2014)：

(4)

互相關(guān)函數(shù)實部與格林函數(shù)虛部成正比，其中A在單一頻率下為常數(shù)，可以通過歸一化的方法消除影響.C為互相關(guān)系數(shù)，i，j從1到3表示不同Z，R，T分量，G為格林函數(shù)，其為單位點力源在一個臺站激發(fā)，然后在另一個臺站接收到的波形.r為兩臺站之間的間距，z=0表示震源和接收點都位于地表.w表示角頻率.Im表示虛部.上述公式正確性，已被廣泛驗證(Wang et al., 2019；Li and Chen,2020；Hu et al., 2020)，當(dāng)然該公式從頻率域變換到時間域具有一定近似性，通常選擇較窄頻段范圍，使頻率近似為常數(shù).

對于一個層狀彈性半空間模型，基于點力源格林函數(shù)垂直分量可以用公式(5)計算(Chen, 1993; Wang et al., 2019; Li et al., 2020)：

(5)

其中：J0表示零階第一類貝塞爾函數(shù)，k表示波數(shù)，r為兩臺站間距，dk表示對波數(shù)積分，g表示核函數(shù)，其是基于各向同性層狀假設(shè)，由廣義反射透射系數(shù)計算得到，其虛部主要表示頻散曲線，其包含地下介質(zhì)物性參數(shù)信息，Chen(1993)給出了具體表達公式.將公式(5)代入到公式(4)可得

(6)

求解式(6)便可由地下模型參數(shù)得到互相關(guān)函數(shù)實部，然后利用互相關(guān)函數(shù)實部信息反演地下介質(zhì)結(jié)構(gòu).為了確定合適的反演策略和方法，首先建立相應(yīng)的目標(biāo)函數(shù)，然后進行敏感性分析.從互相關(guān)函數(shù)中的核函數(shù)g的計算可知，影響互相關(guān)函數(shù)的因素較多，主要有層厚和層位深度，地下介質(zhì)物性參數(shù)(主要指縱橫波速度以及密度)，震源與接收點距離，頻率，以及濾波器種類等.對于資源探測，地下介質(zhì)物性參數(shù)及其所在位置較為重要，因此本研究主要考慮物性參數(shù)和層厚敏感性.

反演目標(biāo)函數(shù)為

(7)

接著考慮目標(biāo)函數(shù)對模型參數(shù)導(dǎo)數(shù)為

(8)

其中，

(9)

表示目標(biāo)函數(shù)梯度.

對于模型m0，敏感性計算公式為

(10)

其中，

d0=G(m0)，

(11)

表示基于模型m0正演計算得到的互相關(guān)函數(shù)，協(xié)方差矩陣應(yīng)由實際數(shù)據(jù)計算得到，在理論分析中，選擇單位矩陣.

接著基于一個遞增模型分析互相關(guān)函數(shù)對模型參數(shù)的敏感性，為后面反演地下介質(zhì)結(jié)構(gòu)提供理論基礎(chǔ).單一頻率敏感性計算頻帶寬度為0.02 Hz.計算歸一化結(jié)果如圖3b所示.

圖3a為橫波模型參數(shù)，縱波與密度由經(jīng)驗公式給出(Shen and Ritzwoller,2016;Wu et al.,2020).圖3b顯示的是0.05～0.3 Hz頻段范圍的互相關(guān)函數(shù)對模型參數(shù)m(依次為縱波速度，密度，橫波速度)敏感性歸一化結(jié)果，從圖3b中可知，互相關(guān)函數(shù)主要對橫波速度敏感，對縱波及密度敏感性較低.當(dāng)同時反演三個參數(shù)時，密度與縱波速度容易陷入局部極值，因此，在反演時僅反演橫波速度，而用經(jīng)驗關(guān)系給出大致的縱波速度與密度值(Pan et al.,2019;Wu et al.,2020).具體經(jīng)驗公式為(Brocher,2005;Shen and Ritzwoller,2016)：

(12)

(13)

本文方法是根據(jù)各向同性層狀假設(shè)，而實際互相關(guān)函數(shù)的提取并無這一假設(shè)，當(dāng)頻率變高時，半空間層狀介質(zhì)假設(shè)已不能滿足實際波場傳播特征，此時需要用更加精細(xì)的方法如有限差分等方法對地下介質(zhì)進行剖分.圖3c和3d敏感性分析表明，高頻對淺部小尺度結(jié)構(gòu)較為敏感，而低頻主要對深度大尺度模型較為敏感，當(dāng)然敏感性具體與模型有關(guān).影響敏感性因素較多，當(dāng)層厚較小，反演層數(shù)和參數(shù)較多時，會使得反演結(jié)果易收斂到局部極值，因此，對于同一深度的模型可采取增大模型厚度的方法來進行反演.本文所要介紹的多尺度反演方法是對同一個深度的模型，在不同反演階段進行不同層厚的劃分，然后用不同頻率的數(shù)據(jù)進行反演，為反演不同尺度結(jié)構(gòu)和深度的模型以及權(quán)衡計算效率提供方法.

圖3 (a)物性參數(shù)模型；(b)不同物性參數(shù)敏感性結(jié)果；(c)不同頻段數(shù)據(jù)敏感性結(jié)果；(d)不同頻段數(shù)據(jù)在深部敏感性結(jié)果

2 數(shù)值試驗

2.1 合成數(shù)據(jù)模型測試

為了驗證上述理論，首先進行理論合成數(shù)據(jù)試驗.設(shè)計圖4a所示背景噪聲模擬和采集系統(tǒng)，采集背景噪聲，并合成互相關(guān)函數(shù)進行反演，初步驗證互相關(guān)函數(shù)反演可行性.圖4a噪聲合成模型中兩臺站分別位于(50 km，0 km)和(-50 km，0 km)處，點力噪聲源均勻分布于60～260 km圓環(huán)內(nèi).噪聲源在一段時間內(nèi)均勻隨機振動，然后合成相應(yīng)的背景噪聲，共60段，其中一段噪聲如圖4b所示.接著對噪聲進行處理，合成互相關(guān)函數(shù)，并得到互相關(guān)函數(shù)實數(shù)部分，然后利用共軛梯度算法反演地下模型，在反演中初始模型選擇速度遞增背景速度場，反演結(jié)果如圖5所示.

圖4 (a)噪聲源與臺站之間的相對位置關(guān)系；(b)模擬的背景噪聲

圖5a為L曲線法確定正則化參數(shù)，由于數(shù)據(jù)誤差中包含權(quán)重因子，因此其值較大.從圖5b、5c、5d可以看出，反演結(jié)果雖然在具體數(shù)值上與真實結(jié)果有一定偏差，但形態(tài)特征基本與真實模型一致.其結(jié)果大體上可以反映出真實的地下結(jié)構(gòu).

2.2 實際數(shù)據(jù)模型測試

接著對一段實際資料進行處理研究，實際資料來源于中國永豐—撫州斷裂地區(qū)(趙延娜等,2015).該區(qū)及附近區(qū)域地下巖石圈軟流圈邊界主要分布范圍為85～150 km范圍之間(Hu and Wang,2000;Priestley and McKenzie,2013;Ouyang et al.,2014;Shan et al.,2017)，不同觀測和定義方法可能會有一定誤差，但地下軟流圈的存在基本已無爭議.此處選擇臺站SC09和SC19噪聲數(shù)據(jù)進行反演，臺站間距約為100 km，使用頻段范圍為0.06～0.11 Hz，共劃分0.06～0.08 Hz，0.07～0.09 Hz，0.09～0.11 Hz三個頻段，初始模型層厚3 km是一遞增模型.將反演結(jié)果與頻散曲線結(jié)果進行對比，驗證方法正確性.反演結(jié)果如圖6所示，其中圖6a為正則化參數(shù)，圖6b為反演結(jié)果，圖7—9為三個頻段下，互相關(guān)函數(shù)比較結(jié)果，其中(a)表示實際互相關(guān)與由初始模型計算得到互相關(guān)比較結(jié)果，圖(b)表示實際互相關(guān)與由反演結(jié)果計算得到互相關(guān)比較結(jié)果，圖(c)表示真實互相關(guān)與初始模型誤差，圖(d)表示真實互相關(guān)與反演結(jié)果誤差.圖10為初始模型計算得到的頻散曲線(藍(lán)線)和反演結(jié)果計算得到的頻散曲線(紅線)與FJ方法(圖10a)和雙臺時頻分析(圖10b)提取的頻散曲線之間的比較結(jié)果，與初始模型相比本文反演結(jié)果計算得到的頻散曲線在對應(yīng)頻段能更好的與實際頻散曲線擬合.從互相關(guān)函數(shù)和頻散曲線比較結(jié)果中可以看出，反演結(jié)果誤差相比于初始誤差是減小的，反演結(jié)果是朝著真實結(jié)果收斂的.從而在實際資料中驗證了本文方法的可行性.

圖6 (a)L曲線示意圖；(b)初始模型與反演結(jié)果示意圖

圖7 頻段為0.06～0.08 Hz時(a)真實數(shù)據(jù)與初始模型合成數(shù)據(jù)比較；(b)真實數(shù)據(jù)與反演結(jié)果合成數(shù)據(jù)比較；(c)真實數(shù)據(jù)與初始模型合成數(shù)據(jù)差；(d)真實數(shù)據(jù)與反演結(jié)果差

圖8 頻段為0.07～0.09 Hz時(a)真實數(shù)據(jù)與初始模型合成數(shù)據(jù)比較;(b)真實數(shù)據(jù)與反演結(jié)果合成數(shù)據(jù)比較;(c)真實數(shù)據(jù)與初始模型合成數(shù)據(jù)差;(d)真實數(shù)據(jù)與反演結(jié)果差

圖9 頻段為0.09～0.11 Hz時(a)真實數(shù)據(jù)與初始模型合成數(shù)據(jù)比較；(b)真實數(shù)據(jù)與反演結(jié)果合成數(shù)據(jù)比較；(c)真實數(shù)據(jù)與初始模型合成數(shù)據(jù)差；(d)真實數(shù)據(jù)與反演結(jié)果差

圖10 本文反演結(jié)果(紅線)和初始模型結(jié)果(藍(lán)線)與基于頻率貝塞爾函數(shù)提取頻散曲線(a)和時頻分析提取頻散曲線(b)對比說明圖

本文敏感性分析表明，由于不同頻段互相關(guān)函數(shù)對不同尺度和深度模型敏感性不同，因此我們可以選擇新的反演方案.在同一臺站間距尺度下，給定一個層厚較大，深度較深的初始模型，接著用較低頻段進行反演，此時可以得到大尺度模型反演結(jié)果.接著對反演結(jié)果重新劃分層厚，并將其作為初始模型，用高頻段信息進行反演，從而得到更加精細(xì)的結(jié)果.與前人多尺度反演方法的不同之處在于本文對層厚進行了重新劃分，可以得到更加精細(xì)的模型結(jié)果.

首先，對原模型劃分層厚，層厚為6 km，用較低頻段進行反演，反演頻段為0.05～0.07 Hz，0.06～0.08 Hz，0.07～0.09 Hz，大尺度模型反演結(jié)果如圖11所示.

從圖11a可以看出本文反演結(jié)果在110 km處左右可見明顯低速異常，基本在該地區(qū)軟流圈范圍內(nèi).頻散曲線對比結(jié)果中紅色點表示由初始模型計算得到的結(jié)果，藍(lán)色點表示由反演結(jié)果計算得到的結(jié)果.從圖11b中可以看出，藍(lán)點與實際提取得到的頻散曲線能較好的擬合.從而說明反演結(jié)果可靠性.

圖11 (a)反演結(jié)果與初始模型示意圖；(b)反演結(jié)果(藍(lán)線)和初始模型(紅線)與頻散曲線比較結(jié)果

為了驗證反演深度可靠性，我們對上面反演結(jié)果(不包含最底層無窮遠(yuǎn)半空間)在對應(yīng)反演頻段下進行敏感性分析.此處協(xié)方差矩陣由實際資料處理得到.從圖12敏感性結(jié)果可以看出，本文所選擇頻段對深部異常結(jié)構(gòu)具有敏感性，說明了反演深度可靠性.在實際反演中反演結(jié)果通常取決于實際資料中的有效信號，而協(xié)方差矩陣可以有效衡量有效信號.

圖12 反演結(jié)果模型不同層厚在三個頻段下深部結(jié)構(gòu)敏感性結(jié)果

2.3 多尺度反演

接著根據(jù)不同頻率對不同層厚和深度敏感性不同，驗證本文模型尺度隨頻率變化而變化的多尺度反演方法的有效性.將前面反演結(jié)果作為初始模型，對原6 km層厚進行重新劃分，劃分為3 km層厚，劃分過程類似于圖14a到圖14b，保留對應(yīng)深度具體數(shù)值不變，然后用較高頻率進行反演.反演結(jié)果如下：

從圖13中可以看出，當(dāng)加入較高頻率后，高頻主要對淺層結(jié)構(gòu)敏感，可以反演得到更加精細(xì)的淺層結(jié)果.由于在反演時我們無法確定基于層狀假設(shè)合適頻段，因此，此處給出了兩個頻段進行反演.由于頻段不一樣，且淺部結(jié)構(gòu)較為復(fù)雜，層狀近似誤差較大，因此在數(shù)值上略有區(qū)別，但兩個高頻段數(shù)據(jù)對淺層反演結(jié)果都可見淺部低速異常結(jié)構(gòu).此低速構(gòu)造可能是該區(qū)斷裂帶的存在導(dǎo)致油氣水或地下巖漿等低速體的聚集導(dǎo)致.通常在一維結(jié)構(gòu)反演中我們可以給定一系列初始模型從而求得相對穩(wěn)定的反演結(jié)果.在二維、三維反演中由于計算量的限制，通常很少采用該方法.本方法可以為復(fù)雜介質(zhì)成像構(gòu)建相對可靠初始模型.

圖13 (a)0.06～0.1 Hz中三個頻段數(shù)據(jù)反演結(jié)果；(b)0.07～0.11 Hz中三個頻段數(shù)據(jù)反演結(jié)果

接著對上述結(jié)果進行解釋.

如圖14所示,對于地下10～20 km處一個異常體，保持臺站間距不變，如圖14b對層厚進行重新劃分，層厚較大時，異常體相對影響較小，層厚較小時，異常體相對影響較大.因此，層厚和臺站間距等均會對結(jié)果產(chǎn)生相應(yīng)影響，而本文多尺度方法可以得到相對準(zhǔn)確的異常體層位信息.另一方面，如圖14b所示，對于低頻長波長信息，較小的層厚占波長比重較小，因此對模型敏感性較小，只有通過增加層厚才能反演得到相對精確的結(jié)果，同時也可以將低頻反演結(jié)果作為高頻反演初始模型，從而反演得到更加精細(xì)的結(jié)果.圖14主要說明合理劃分層厚，選擇合適的臺站間距、頻率對于研究不同尺度地下介質(zhì)真實結(jié)構(gòu)具有較為重要作用.這一結(jié)論和假設(shè)，即成像頻率逐漸變大，尺度逐漸變小，廣泛存在于彈性波成像中，如大尺度地震成像，地震勘探，超聲醫(yī)學(xué)成像中.

圖14理論分析表明，對于各向同性層狀介質(zhì)假設(shè)方法，僅在地下介質(zhì)完全均勻情況下才有可能反演出相對真實的參數(shù)結(jié)果.當(dāng)?shù)叵潞挟惓ｓw時，反演參數(shù)值是地下介質(zhì)參數(shù)的一個平均結(jié)果，層厚的劃分對模型分辨率結(jié)果具有一定影響.

圖14 (a)層厚為20 km示意圖；(b)層厚為10 km不同頻率波長數(shù)據(jù)對地下介質(zhì)響應(yīng)示意圖

3 結(jié)論

本文在各向同性層狀介質(zhì)假設(shè)下基于背景噪聲互相關(guān)函數(shù)成像提出了層厚和反演深度隨頻率變化的多尺度成像方法.本文反演主要是在單一臺站間距下，利用低頻反演大尺度模型，然后對低頻反演結(jié)果重新劃分層厚添加高頻成分反演得到更精細(xì)的結(jié)果，這樣可以權(quán)衡計算量、模型分辨率、反演穩(wěn)定性和反演深度之間的關(guān)系.本文方法為今后二維、三維反演獲得相對可靠初始速度模型以及為多尺度反演提供了相應(yīng)策略.同時對僅能布設(shè)少量地震臺站的地區(qū)如月球火星等地外天體內(nèi)部結(jié)構(gòu)探測具有較為重要的意義.

致謝感謝2020北大暑期課程以及課題組提供的噪聲處理開源程序，同時感謝課題組陳舉慶、閆英偉、康佳語、徐建寬和李正波等人在文章修改過程中提供的文字潤色和校正服務(wù).