壓縮感知高分辨率接收函數(shù)疊加成像及其應(yīng)用

白蘭淑， 吳慶舉， 張瑞青

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

0 引言

巖石圈精細(xì)結(jié)構(gòu)是研究地球深部動(dòng)力學(xué)機(jī)制的重要依據(jù).接收函數(shù)成像對(duì)臺(tái)站下方速度界面縱向分辨能力高，是有效的深部界面結(jié)構(gòu)研究手段.傳統(tǒng)接收函數(shù)共轉(zhuǎn)換點(diǎn)(Common Conversion Point，CCP)疊加成像方法是基于層狀界面假設(shè)的準(zhǔn)二維疊前成像方法，具有可靠、計(jì)算成本低、實(shí)用性強(qiáng)的優(yōu)點(diǎn)，近20年來(lái)得到廣泛的研究和應(yīng)用(吳慶舉和曾融生，1998；Zhu and Kanamori, 2000；吳建平等，2001；陳九輝，2007；劉啟元等, 2008；李永華等，2009；Tian et al.，2011；Zhang et al.，2014；Wang et al.，2016；武巖等，2018).

接收函數(shù)成像方法縱向分辨能力高，而其水平分辨能力很大程度上取決于地表觀測(cè)臺(tái)站分布.臺(tái)站越密集，水平分辨率相對(duì)更高，反過(guò)來(lái)，臺(tái)站越稀疏，由于需采用較大的疊加半徑，水平分辨率有限.當(dāng)臺(tái)站密度不均勻時(shí)，為了兼顧臺(tái)站稀疏區(qū)域的成像效果，常常統(tǒng)一選取較大尺度疊加半徑進(jìn)行疊加成像，以獲得更為光滑、連續(xù)的成像界面.然而，這不僅造成原本臺(tái)站密集區(qū)下方成像界面的水平分辨率的損失，還會(huì)使成像界面產(chǎn)生明顯的階梯狀疊加痕跡.

經(jīng)典CCP成像方法基于層狀速度模型假設(shè)，通過(guò)層狀界面轉(zhuǎn)換波射線追蹤進(jìn)行接收函數(shù)成像，對(duì)高角度、陡傾角界面的成像可能存在一定偏差.基于波動(dòng)理論的波動(dòng)方程疊后偏移方法(Chen et al.，2005)和疊前偏移方法(Shang et al.，2012；Jiang et al.，2019)有效克服了經(jīng)典CCP方法無(wú)法對(duì)劇烈起伏界面、地下復(fù)雜結(jié)構(gòu)進(jìn)行準(zhǔn)確成像的問(wèn)題，但其三維實(shí)際應(yīng)用相比常規(guī)CCP成像方法需要巨大的計(jì)算成本，且對(duì)臺(tái)站密度要求非常高.

如果臺(tái)陣分布稀疏、不均勻，通過(guò)對(duì)所有臺(tái)站接收函數(shù)進(jìn)行規(guī)則化、加密，可以有效提升成像分辨率.Zhang和Zheng(2015)采用三次樣條插值函數(shù)對(duì)接收函數(shù)進(jìn)行了加密、規(guī)則化，并應(yīng)用于鄂爾多斯地區(qū)地震臺(tái)陣觀測(cè)數(shù)據(jù)，結(jié)果表明加密、規(guī)則化的接收函數(shù)對(duì)界面的橫向起伏形態(tài)具有更高的分辨率.Song等(2017)在我國(guó)南部深部結(jié)構(gòu)研究中，首先基于徑向基函數(shù)(Radial Basis Function) 對(duì)原始地震數(shù)據(jù)進(jìn)行了加密、規(guī)則化處理，再提取接收函數(shù)并進(jìn)行成像，成像質(zhì)量獲得有效提升.Hu等(2018)采用Stretching-and-Squeezing插值方法對(duì)接收函數(shù)進(jìn)行了規(guī)則化插值處理，并進(jìn)一步通過(guò)對(duì)插值的接收函數(shù)進(jìn)行CCP成像提升了成像質(zhì)量.然而，關(guān)于地震數(shù)據(jù)或接收函數(shù)的加密、規(guī)則化方法在數(shù)學(xué)、物理含義層面上仍待進(jìn)一步深入討論和研究.

壓縮感知理論認(rèn)為具有稀疏性的信號(hào)可以通過(guò)求解稀疏促進(jìn)算法從少量觀測(cè)數(shù)據(jù)中重構(gòu)出高精度的原始信號(hào)(Donoho，2006).由于地震信號(hào)具有特定的頻率、方向、空間分布特征，因此在一些數(shù)據(jù)變換域表現(xiàn)出稀疏特性.在勘探地震學(xué)領(lǐng)域，壓縮感知理論在數(shù)據(jù)采集、數(shù)據(jù)重建、高分辨率地下結(jié)構(gòu)成像等方面得到廣泛應(yīng)用(Herrmann and Hennenfent，2008；唐剛，2010；白蘭淑等，2014；Mosher et al.，2014；Bai et al.，2017).天然地震學(xué)領(lǐng)域，壓縮感知理論在地震數(shù)據(jù)采集和重建技術(shù)等研究中得到了新的應(yīng)用(Bai et al.，2020).

本文從壓縮感知理論出發(fā)，研究可以有效減少橫向分辨率損失的壓縮感知高分辨率接收函數(shù)成像方法.本方法充分利用地震信號(hào)具備的稀疏特性，首先在稀疏變換域通過(guò)壓縮感知稀疏促進(jìn)算法對(duì)接收函數(shù)進(jìn)行重建，獲得規(guī)則化、加密的“虛擬臺(tái)陣”接收函數(shù)，接著對(duì)其進(jìn)行時(shí)深轉(zhuǎn)換、疊前振幅補(bǔ)償以及小尺度疊加半徑疊加成像，有效減少常規(guī)CCP方法采用大疊加半徑導(dǎo)致的橫向分辨率損失.本文方法不僅可以保證臺(tái)站密集區(qū)成像分辨率，同時(shí)能夠充分利用地震信號(hào)的稀疏性有效提升臺(tái)站稀疏區(qū)成像質(zhì)量.本研究方法結(jié)合壓縮感知理論和經(jīng)典CCP疊加成像方法實(shí)現(xiàn)高分辨率接收函數(shù)成像，具有明確的數(shù)學(xué)理論支撐和相應(yīng)的物理含義，華北克拉通科學(xué)臺(tái)陣數(shù)據(jù)的應(yīng)用結(jié)果驗(yàn)證了其可靠性和有效性.本文方法在地球深部精細(xì)結(jié)構(gòu)研究中具有廣闊的應(yīng)用前景.

1 方法

本研究基于地震信號(hào)具有稀疏性這一假設(shè)，首先根據(jù)實(shí)際臺(tái)陣分布情況建立密集、規(guī)則化分布的“虛擬臺(tái)陣”(以下引號(hào)省略)，利用實(shí)際數(shù)據(jù)計(jì)算得到的接收函數(shù)，通過(guò)壓縮感知稀疏促進(jìn)算法重構(gòu)出密集、規(guī)則化虛擬臺(tái)陣接收函數(shù).然后，按照經(jīng)典CCP疊加成像方法，采用層速度模型，將所有虛擬臺(tái)站接收函數(shù)沿著射線路徑進(jìn)行時(shí)深轉(zhuǎn)換.接著，對(duì)所有成像點(diǎn)的成像值進(jìn)行疊前振幅校正.最后，選取小尺度疊加半徑對(duì)所有地下成像點(diǎn)進(jìn)行共轉(zhuǎn)換點(diǎn)疊加獲得高分辨率成像剖面.

1.1 壓縮感知稀疏促進(jìn)算法重構(gòu)接收函數(shù)

壓縮感知理論認(rèn)為稀疏信號(hào)可以通過(guò)稀疏促進(jìn)算法從少量觀測(cè)數(shù)據(jù)中重構(gòu)出高精度的原始信號(hào)(Donoho，2006).因此，接收函數(shù)是否具有稀疏性是壓縮感知理論能否成功應(yīng)用于接收函數(shù)重構(gòu)和高分辨率成像的重要前提條件.

遠(yuǎn)震事件產(chǎn)生的地震波傳播至臺(tái)站下方，其地震波場(chǎng)在空間上是連續(xù)變化的，地表的地震波場(chǎng)自然也呈空間連續(xù)變化.因此，距離相近的臺(tái)站所記錄的地震信號(hào)通常具有相似的走時(shí)、波形形態(tài)和頻帶特性.數(shù)學(xué)上，接收函數(shù)通過(guò)同一臺(tái)站事件波形的徑向(或切向)分量和垂直分量的反褶積運(yùn)算獲得.因此，接收函數(shù)理論上也與觀測(cè)的地震記錄一樣具有相似的走時(shí)、波形形態(tài)和頻帶特性，同一轉(zhuǎn)換震相對(duì)應(yīng)的接收函數(shù)在剖面上的局部小尺度空間范圍內(nèi)斜率變化穩(wěn)定、近乎呈線性.針對(duì)接收函數(shù)的這一特征，我們認(rèn)為其在一些數(shù)據(jù)變換域可以進(jìn)行稀疏表示，如傅里葉、Radon、Curvelet等變換域.一旦有了接收函數(shù)具有稀疏特性的這一合理假定，就可以通過(guò)壓縮感知稀疏促進(jìn)反演重構(gòu)出加密、規(guī)則化的虛擬臺(tái)陣接收函數(shù).

下面，我們按照壓縮感知理論建立對(duì)原本不規(guī)則、不均勻分布地震臺(tái)陣接收函數(shù)進(jìn)行重構(gòu)的基本框架.

地表水平面為連續(xù)的二維空間，地震波從地下傳播至地表之后在二維地表水平面上的地震波場(chǎng)為三維連續(xù)信號(hào)，其中兩個(gè)維度是地表水平面空間坐標(biāo)，一個(gè)維度是時(shí)間.從該三維連續(xù)地震波場(chǎng)提取的接收函數(shù)也為相同維度上的地震信號(hào)，我們用S表示.如果采用非常密集的、規(guī)則分布的二維“虛擬”地震臺(tái)陣記錄這三維連續(xù)接收函數(shù)信號(hào)S，并對(duì)時(shí)間方向進(jìn)行離散化，可以得到離散化接收函數(shù)信號(hào)f∈RN(N=Nx×Ny×Nt)，我們稱其為原始信號(hào).其中，Nx、Ny分別為水平面兩個(gè)正交方向上的臺(tái)站個(gè)數(shù)，Nt為每個(gè)虛擬臺(tái)站記錄的時(shí)間采樣總數(shù).原始信號(hào)f在稀疏變換域有如下表示：

f=Ψu,

‖u‖0=k,

(1)

其中，Ψ代表稀疏變換矩陣，u代表稀疏變換基對(duì)應(yīng)的系數(shù)向量，‖·‖0為L(zhǎng)0范數(shù)，k表示向量(或矩陣)的非零元素個(gè)數(shù).從而，u中非零元素越少，0范數(shù)越小，u越稀疏.本文采用三維傅里葉變換作為稀疏表達(dá)接收函數(shù)的數(shù)據(jù)變換.由于實(shí)際臺(tái)陣位置基本不落在規(guī)則網(wǎng)格點(diǎn)上，因此我們特別采用了三維非均勻傅里葉變換(Non-uniform Fast Fourier Transform，NUFFT)(nufft3d-v1.3 package 2011. https:∥github.com/zgimbutas/nufft3d).

實(shí)際布設(shè)二維臺(tái)陣時(shí)，地面臺(tái)站往往分布不規(guī)則、不均勻，有些區(qū)域密集，而有些區(qū)域稀疏.我們將從實(shí)測(cè)臺(tái)站記錄的地震數(shù)據(jù)中提取的接收函數(shù)稱為“不完整”數(shù)據(jù)，并用y∈RM(M?N)表示.那么

y=Φf=ΦΨu=Θu,

(2)

其中,Θ∈RM×N(M?N)為測(cè)量矩陣，其必須服從“有限等距性質(zhì)(Restricted Isometry Property，RIP)”(Candès，2008).有限等距性質(zhì)可以描述近似正交矩陣的非正交程度，從理論上保證k-稀疏信號(hào)能由M個(gè)測(cè)量值y重構(gòu)出長(zhǎng)度為N的原始數(shù)據(jù)f.隨機(jī)高斯矩陣與大多數(shù)正交基構(gòu)成的矩陣不相關(guān)，以高概率滿足RIP.

對(duì)于不規(guī)則分布的地震臺(tái)陣而言，可以認(rèn)為其觀測(cè)數(shù)據(jù)在二維水平空間維度上是近隨機(jī)的.實(shí)際布設(shè)中，也存在較多的地震臺(tái)陣采用近等間距的規(guī)則布設(shè).關(guān)于地震觀測(cè)點(diǎn)規(guī)則或不規(guī)則分布時(shí)的地震數(shù)據(jù)重構(gòu)效果，已有大量通過(guò)模擬數(shù)據(jù)或?qū)嶋H數(shù)據(jù)開(kāi)展的研究和分析(尤其是勘探地震學(xué)領(lǐng)域).研究表明只要地震數(shù)據(jù)具有稀疏性，均可以高概率獲得高精度重構(gòu)結(jié)果(Herrmann and Hennenfent，2008；唐剛，2010；白蘭淑等，2014；黃小剛等，2014；Bai et al.，2017，2020).因此，我們假定無(wú)論是相對(duì)等間距、規(guī)則布設(shè)，還是不等間距、隨機(jī)布設(shè)，都可以實(shí)現(xiàn)虛擬臺(tái)陣接收函數(shù)的成功重構(gòu).

我們的目標(biāo)就是從上述不完整接收函數(shù)y中重建出密集、規(guī)則的虛擬臺(tái)陣接收函數(shù)f.在公式(2)中，從低維空間上的觀測(cè)值y中求解獲得原始高維信號(hào)f，未知量個(gè)數(shù)大于已知量個(gè)數(shù)，是一個(gè)欠定方程，必須施加約束條件才可以獲得定解.基于原始信號(hào)的稀疏性假設(shè)，壓縮感知理論提供了對(duì)重構(gòu)信號(hào)施加稀疏約束并通過(guò)稀疏促進(jìn)算法求解上述欠定方程的解決方法.公式(2)的最稀疏解，可以通過(guò)求解如下問(wèn)題獲得：

(3)

求解公式(3)中的稀疏促進(jìn)問(wèn)題是一個(gè)經(jīng)典的NP-hard非凸優(yōu)化問(wèn)題，對(duì)于稀疏信號(hào)而言可以通過(guò)用L1范數(shù)代替L0范數(shù)變成凸優(yōu)化問(wèn)題，并進(jìn)一步轉(zhuǎn)化為如下非約束化的子問(wèn)題：

(4)

公式(4)中，等式右側(cè)第一項(xiàng)為數(shù)據(jù)殘差二范數(shù)約束項(xiàng)，第二項(xiàng)為L(zhǎng)1范數(shù)稀疏約束項(xiàng)，λ為權(quán)重系數(shù).本文用基于Landweber下降法的冷卻閾值迭代技術(shù)求解上述問(wèn)題(Herrmann and Hennenfent，2008；唐剛，2010)，以三維傅里葉域作為稀疏表達(dá)接收函數(shù)的變換域.具體迭代算法如表1所示.

表1中，Tλ(·)為硬閾值函數(shù)，滿足

(5)

1.2 疊前振幅校正

通過(guò)稀疏促進(jìn)算法重構(gòu)虛擬臺(tái)陣接收函數(shù)時(shí)，由于虛擬臺(tái)陣的臺(tái)間距比實(shí)際布設(shè)臺(tái)間距小很多，重構(gòu)的虛擬臺(tái)站接收函數(shù)振幅偏弱，尤其是在臺(tái)站極為稀疏的區(qū)域.如果直接按照經(jīng)典方法進(jìn)行簡(jiǎn)單的疊加平均操作，會(huì)導(dǎo)致成像界面能量不均衡.

針對(duì)上述問(wèn)題，我們提出一種疊前振幅校正公式，具體表達(dá)式如下：

α=(max(e-min(disthori(xima,xista)),c))-1,

(6)

式中，α為校正因子，xima和xista分別代表地下第ima個(gè)成像點(diǎn)和第ista個(gè)臺(tái)站的坐標(biāo)，disthori(xima,xista)代表這兩點(diǎn)之間的水平距離.min(disthori(xima,xista))代表距離第ima個(gè)成像點(diǎn)最近的地震臺(tái)站與該成像點(diǎn)之間的水平距離.該值越大，說(shuō)明該成像點(diǎn)附近臺(tái)站越稀疏，校正因子α取值越大.c為小于1的一個(gè)常數(shù)，用來(lái)設(shè)定振幅校正因子的上界.由于不同地震事件篩選出的地震臺(tái)站分布有所差異，我們對(duì)每個(gè)地震事件接收函數(shù)依次進(jìn)行疊加前的振幅校正.

表1 迭代閾值稀疏促進(jìn)算法Table 1 The iterative thresholding method for the sparsity inversion

1.3 小尺度疊加半徑疊加成像

疊前振幅補(bǔ)償后，我們采用小尺度疊加半徑進(jìn)行傳統(tǒng)CCP疊加成像(ccp1.0. tar package 2006. https:∥www.eas.slu.edu/People/LZhu/downloads/ccp1.0.tar)，獲得高分辨率地下結(jié)構(gòu)成像結(jié)果.疊加半徑大小通常需要綜合實(shí)際臺(tái)陣分布、虛擬臺(tái)陣網(wǎng)格間距等因素進(jìn)行選取，通?？刹捎锰摂M臺(tái)陣臺(tái)間距相當(dāng)?shù)某叨?

(7)

式中，divsta為虛擬臺(tái)站接收函數(shù)時(shí)深轉(zhuǎn)換之后的(單條射線對(duì)應(yīng)的)深度域振幅，xima代表地下第ima個(gè)成像點(diǎn)的坐標(biāo)，ievt、ivsta分別代表地震事件和虛擬臺(tái)站編號(hào)，Nray為以xima為中心的疊加箱內(nèi)總射線條數(shù)，Image為振幅校正后疊加得到的像.

1.4 壓縮感知高分辨率接收函數(shù)疊加成像步驟

本文結(jié)合壓縮感知理論和經(jīng)典CCP接收函數(shù)疊加成像的方法具體實(shí)現(xiàn)流程如下:

(1)提取遠(yuǎn)震事件接收函數(shù)，篩選高信噪比接收函數(shù)臺(tái)站數(shù)較多的遠(yuǎn)震事件接收函數(shù)；

(2)根據(jù)實(shí)際地震臺(tái)陣分布，在矩形區(qū)域內(nèi)建立密集、規(guī)則化虛擬地震臺(tái)陣；

(3)針對(duì)每一個(gè)地震事件，利用篩選出的接收函數(shù)，通過(guò)表1所示的壓縮感知稀疏促進(jìn)反演算法重構(gòu)虛擬臺(tái)陣接收函數(shù)；當(dāng)篩選出的地震事件共Nevt個(gè)，虛擬臺(tái)陣共Nvsta個(gè)，則所有地震事件重構(gòu)的總接收函數(shù)共Nevt×Nvsta條.

(4)綜合區(qū)域構(gòu)造背景、實(shí)際觀測(cè)臺(tái)陣分布、虛擬臺(tái)陣分布等因素，建立三維地下成像點(diǎn)網(wǎng)格；

(5)設(shè)定合適的小尺度長(zhǎng)方體疊加箱，對(duì)所有篩選的遠(yuǎn)震事件依次進(jìn)行(6)—(7)的操作；

(6)在水平層狀界面假設(shè)下，根據(jù)參考速度模型和遠(yuǎn)震Ps轉(zhuǎn)換波對(duì)應(yīng)的射線參數(shù)計(jì)算每個(gè)臺(tái)站對(duì)應(yīng)的轉(zhuǎn)換波射線路徑.對(duì)所有臺(tái)站接收函數(shù)時(shí)間序列從零時(shí)刻開(kāi)始，從地表向下沿著轉(zhuǎn)換波射線路徑放在射線所落入的成像點(diǎn)上，即時(shí)深轉(zhuǎn)換.

(7)以每個(gè)地震事件視為處理單元，根據(jù)公式(6)式對(duì)虛擬臺(tái)陣接收函數(shù)依次進(jìn)行成像點(diǎn)振幅校正；

(8)對(duì)每個(gè)地下成像點(diǎn)圈定疊加箱范圍，將落入該點(diǎn)疊加箱里的成像能量進(jìn)行疊加取平均，獲得該成像點(diǎn)的像.

(9)獲得三維地下結(jié)構(gòu)成像剖面.

經(jīng)典方法通過(guò)共成像點(diǎn)疊加箱進(jìn)行疊加成像主要有兩個(gè)目的.首先，疊加可以有效提高成像結(jié)果信噪比.其次，由于臺(tái)站分布稀疏或者地下結(jié)構(gòu)復(fù)雜，有些成像區(qū)域沒(méi)有射線經(jīng)過(guò)、或者非常稀少.因此，利用較大尺度的疊加箱可以得到更加連續(xù)、光滑的成像剖面.然而，如前所述，采用過(guò)大的疊加半徑，會(huì)嚴(yán)重影響臺(tái)站密集區(qū)的原本較高的成像分辨率.本文與經(jīng)典方法的不同點(diǎn)在于，我們首先在數(shù)據(jù)域?qū)邮蘸瘮?shù)進(jìn)行基于壓縮感知理論下的加密、規(guī)則化重構(gòu)，再在成像域用小疊加半徑進(jìn)行疊加成像，而不像經(jīng)典方法直接在成像域進(jìn)行大疊加窗光滑疊加.從而，本文方法最大的優(yōu)點(diǎn)在于：不僅可以很好地保留臺(tái)站分布密集區(qū)域的原始數(shù)據(jù)含有的細(xì)節(jié)特征，還可以充分利用實(shí)際數(shù)據(jù)中包含的數(shù)據(jù)特征，重構(gòu)出臺(tái)站稀疏區(qū)域臺(tái)站的接收函數(shù)主要特征.

2 方法測(cè)試與實(shí)際應(yīng)用

2.1 數(shù)值模擬測(cè)試

我們?cè)O(shè)計(jì)了如圖1所示的2.5維模型，該模型介質(zhì)參數(shù)沿東西向呈不均勻變化，沿南北方向保持不變. 模型中存在一個(gè)界面，其形態(tài)呈中段傾斜、兩側(cè)水平，界面深度最左端達(dá)70 km，最右端為20 km.上層縱橫波速度和密度分別為5500 m·s-1、2800 m·s-1、2.332 g·cm-3，下層縱橫波速度和密度分別為7289 m·s-1、4217 m·s-1、2.9 g·cm-3.遠(yuǎn)震波場(chǎng)從模型東側(cè)底部以平面波的形式入射，入射角20°、后方位角90°.地面接收臺(tái)站共71×41=2991個(gè)，等間距布設(shè)，東西、南北向臺(tái)間距分別為4.29 km、5 km，臺(tái)站布設(shè)區(qū)域東西、南北向跨度分別達(dá)300 km、200 km.

圖1 數(shù)值測(cè)試中采用的速度模型Fig.1 Velocity model used in the synthetic test

本研究采用SEM-FK方法模擬遠(yuǎn)震事件的地震波場(chǎng)，方法理論介紹詳見(jiàn)Tong等(2014).我們?cè)趫D2a—c中分別展示了正演模擬獲得的CC′測(cè)線上71個(gè)密集、規(guī)則臺(tái)站(圖1中黑色圓點(diǎn))的垂向、東西向、南北向地震波場(chǎng)記錄.記錄中第一個(gè)到達(dá)震相為P波，其后到達(dá)的第二個(gè)震相為Ps轉(zhuǎn)換波.這些規(guī)則、密集臺(tái)站的地震記錄是我們后續(xù)從中隨機(jī)抽取稀疏臺(tái)站記錄并作為“實(shí)測(cè)臺(tái)站”數(shù)據(jù)的基礎(chǔ)數(shù)據(jù)，同時(shí)也是結(jié)果對(duì)比分析中所需的參考數(shù)據(jù).首先對(duì)規(guī)則、密集臺(tái)陣記錄的兩個(gè)水平分量進(jìn)行轉(zhuǎn)換獲得了徑向、切向水平記錄，分別如圖2d—e所示.

接著，利用徑向和垂向分量數(shù)據(jù)，通過(guò)時(shí)間域反褶積計(jì)算，獲得了規(guī)則、密集的接收函數(shù)(此后簡(jiǎn)稱為參考接收函數(shù))，具體結(jié)果展示在圖3中.如圖3a所示的臺(tái)站平面分布圖中，灰色空心圓代表2991個(gè)規(guī)則、密集臺(tái)站，用于記錄規(guī)則、密集的地震波場(chǎng)數(shù)據(jù).我們建立位置與這些規(guī)則、密集臺(tái)站一致的虛擬臺(tái)陣，以便于進(jìn)行波形對(duì)比；紅色圓圈代表從2991條地震臺(tái)站中隨機(jī)抽選出的145個(gè)臺(tái)站，占比約4.8%，在數(shù)值測(cè)試中視其為“實(shí)測(cè)臺(tái)站”，并將這些實(shí)測(cè)臺(tái)站接收的記錄用作實(shí)際觀測(cè)數(shù)據(jù)；虛擬臺(tái)站位置與最初布設(shè)的規(guī)則、密集臺(tái)站位置一致.藍(lán)色實(shí)心圓代表CC′測(cè)線上的71個(gè)虛擬臺(tái)站，紅色實(shí)心圓代表CC′測(cè)線10 km范圍內(nèi)的實(shí)測(cè)臺(tái)站.我們采用1.1節(jié)中介紹的壓縮感知稀疏促進(jìn)重構(gòu)方法對(duì)145個(gè)隨機(jī)臺(tái)站的接收函數(shù)進(jìn)行了密集、規(guī)則化重構(gòu)，最終獲得了2991個(gè)虛擬臺(tái)站接收函數(shù)高精度重構(gòu)結(jié)果.圖3b—d分別展示了CC′測(cè)線10 km范圍內(nèi)的實(shí)測(cè)臺(tái)站(紅色實(shí)心圓)的接收函數(shù)、重構(gòu)的接收函數(shù)以及參考接收函數(shù)，圖中零時(shí)刻對(duì)應(yīng)P波到時(shí).通過(guò)對(duì)比重構(gòu)的接收函數(shù)和參考接收函數(shù)波形，雖然僅利用4.8%臺(tái)站的數(shù)據(jù)量，但由于接收函數(shù)具有良好的稀疏性，因此仍然可以高精度重構(gòu)規(guī)則化、密集虛擬臺(tái)陣的接收函數(shù).

圖2 正演模擬獲得的CC′測(cè)線上的波場(chǎng)記錄(a) Z軸垂向分量； (b) 東西向分量； (c) 南北向分量； (d) 徑向分量； (e) 切向分量.Fig.2 The seismic record of the single line along CC′ obtained through forward modeling(a) Z-direction vertical component; (b) EW-direction component; (c) NS-direction component; (d) Radial component; (e) Tangential component.

圖3 壓縮感知稀疏促進(jìn)算法獲得的接收函數(shù)重構(gòu)結(jié)果(a) 臺(tái)站平面分布圖.灰色空心圓代表2991個(gè)規(guī)則、密集臺(tái)站，用于獲得規(guī)則、密集的地震波場(chǎng)記錄，虛擬臺(tái)站位置與這些臺(tái)站位置一致.紅色圓圈(包括實(shí)心和空心)代表我們從2991個(gè)地震臺(tái)站中隨機(jī)抽選出的145個(gè)臺(tái)站，占比約4.8%，我們視其為“實(shí)測(cè)臺(tái)站”，這些實(shí)測(cè)臺(tái)站接收記錄用作實(shí)際觀測(cè)數(shù)據(jù).藍(lán)色實(shí)心圓代表CC′測(cè)線上的71個(gè)虛擬臺(tái)站，紅色實(shí)心圓代表該條測(cè)線10 km范圍內(nèi)的實(shí)測(cè)臺(tái)站； (b) CC′測(cè)線附近實(shí)測(cè)臺(tái)站(圖3a中紅色實(shí)心圓)的接收函數(shù)； (c) CC′測(cè)線虛擬臺(tái)站處重構(gòu)的接收函數(shù)； (d) CC′ 測(cè)線虛擬臺(tái)站處參考接收函數(shù).Fig.3 Receiver functions reconstructed by CS-based sparsity-promoting algorithm(a) Distribution of the stations. Grey hollow circles represent 2991 dense, regular stations that obtains regular, dense seismic records. The virtual stations are all located in the same positions. The red circles (including solid and hollow circles) represents 145 randomly distributed stations, which are chosen from 2991 seismic stations with the proportion about 4.8%, are regarded as “real” observation stations and the seismic record of these real stations are used as real observed data. The blue solid circles are 71 virtual stations located along the CC′ line, and the red solid circles are the real stations located within 10 km from the CC′ line; (b) Receiver functions of the real stations near the CC′ line; (c) Reconstructed receiver functions of the virtual stations; (d) Reference receiver functions of the virtual station position along the CC′ line.

圖4 壓縮感知接收函數(shù)成像方法的數(shù)值模擬測(cè)試結(jié)果及對(duì)比(a) 經(jīng)典方法利用2991個(gè)規(guī)則、密集臺(tái)站真實(shí)接收函數(shù)獲得的結(jié)果(疊加半徑=2 km)； (b) 壓縮感知接收函數(shù)成像方法利用145個(gè)臺(tái)站數(shù)據(jù)重構(gòu)的2991個(gè)虛擬臺(tái)站接收函數(shù)進(jìn)行成像獲得的結(jié)果(疊加半徑=2 km)； (c—e) 經(jīng)典方法利用145個(gè)臺(tái)站(圖3a中紅色 圓圈所示)接收函數(shù)、并分別采用疊加半徑=2 km、5 km、10 km獲得的結(jié)果.Fig.4 Imaging results of the synthetic test using the Compressive-Sensing-based imaging method and the CCP method(a) Imaging results of the CCP method using the receiver functions of 2991 regular, dense stations (stacking radius=2 km); (b) Imaging results of the Compressive-Sensing-based imaging method using the receiver functions of 2991 virtual stations that reconstructed from 145 real station data (stacking radius=2 km); (c—d) Imaging results of the CCP method using the receiver functions of 145 real stations with stacking radius=2 km, 5 km, 10 km respectively.

圖5 中國(guó)地震局地球物理研究所華北科學(xué)臺(tái)陣與密集、規(guī)則化虛擬臺(tái)陣分布(a) 中國(guó)地震局地球物理研究所科學(xué)臺(tái)陣與虛擬臺(tái)陣分布，分別用紫色、黑色三角形表示； (b) 各臺(tái)站使用的接收函數(shù)條數(shù).顏色越深， 表示所使用接收函數(shù)越多； (c—d) 第40個(gè)和第62個(gè)地震事件所使用的接收函數(shù)對(duì)應(yīng)臺(tái)站分布.Fig.5 Distribution of the ChinArray stations and the dense, regularized virtual stations(a) Distribution of the ChinArray stations and virtual stations that used in the research (marked with purple and black triangles respectively); (b) Number of receiver functions used for each station. The darker the station, the higher the number of receiver functions; (c—d) Distributions of the stations of the receiver functions used for the 40th and 62th event.

圖6 本文選用的遠(yuǎn)震事件震源分布Fig.6 Distribution of teleseismic events chosen in the real data application

進(jìn)一步，我們對(duì)重構(gòu)后的2991條虛擬臺(tái)站接收函數(shù)進(jìn)行成像獲得了高分辨率成像結(jié)果.為了對(duì)比分析，我們也采用經(jīng)典CCP疊加成像方法分別對(duì)2991條參考接收函數(shù)、145條稀疏臺(tái)站接收函數(shù)進(jìn)行了成像.對(duì)于經(jīng)典方法的參考接收函數(shù)疊加成像我們采用了2 km疊加半徑，成像結(jié)果沿CC′測(cè)線的剖面如圖4a所示.成像結(jié)果中轉(zhuǎn)換界面清晰，與圖1中真實(shí)界面的位置、形態(tài)一致；經(jīng)典方法的稀疏臺(tái)站接收函數(shù)的疊加成像分別采用了2 km、5 km、10 km的疊加半徑，成像結(jié)果分別如圖4c—e所示.疊加半徑為2 km時(shí)，由于臺(tái)站分布過(guò)于稀疏，界面很多段沒(méi)有得到成像.當(dāng)疊加半徑增加至5 km時(shí)，斜面的成像結(jié)果相對(duì)更連續(xù)，但成像界面仍然不夠清晰，且疊加半徑的擴(kuò)大導(dǎo)致斜面寬度變寬，水平成像分辨率有所損失.當(dāng)疊加半徑達(dá)10 km時(shí)，斜面的橫向分辨率進(jìn)一步降低.這些結(jié)果表明當(dāng)臺(tái)站過(guò)于稀疏時(shí)，經(jīng)典方法采用大尺度疊加半徑進(jìn)行疊加所能提升的成像質(zhì)量有限，且很容易導(dǎo)致強(qiáng)起伏界面的橫向分辨率損失；我們利用壓縮感知稀疏促進(jìn)算法重構(gòu)的接收函數(shù)，采用與真實(shí)接收函數(shù)成像相同的2 km疊加半徑進(jìn)行了CCP疊加成像，成像結(jié)果剖面如圖4b所示.壓縮感知重構(gòu)接收函數(shù)的成像結(jié)果與真實(shí)接收函數(shù)成像結(jié)果近乎一致，陡傾角斜面得到準(zhǔn)確成像，成像分辨率相比直接進(jìn)行稀疏臺(tái)站接收函數(shù)成像的結(jié)果得到顯著提高.數(shù)值測(cè)試結(jié)果很好地驗(yàn)證了壓縮感知接收函數(shù)成像方法“基于地震波場(chǎng)具有稀疏性這一前提重構(gòu)出密集、規(guī)則化接收函數(shù)，再采用小尺度疊加半徑進(jìn)行共轉(zhuǎn)換點(diǎn)疊加成像”的處理模式可以有效提高成像分辨率，提升成像質(zhì)量.

2.2 實(shí)際應(yīng)用實(shí)例

本研究進(jìn)一步將基于壓縮感知重構(gòu)技術(shù)的接收函數(shù)疊加成像方法應(yīng)用于華北克拉通觀測(cè)地震數(shù)據(jù)，獲得了華北克拉通深部結(jié)構(gòu)成像剖面.需要說(shuō)明的是，我們僅討論本文研究方法對(duì)地下界面成像分辨率和成像質(zhì)量的提升效果，不分析深部結(jié)構(gòu)特征及其所指示的動(dòng)力學(xué)問(wèn)題.

中國(guó)地震局地球物理研究所在2006年11月—2009年9月期間，在我國(guó)華北地區(qū)布設(shè)了200多個(gè)流動(dòng)地震臺(tái)站(圖5a中紫色三角形).該地區(qū)臺(tái)陣布設(shè)采用點(diǎn)面相結(jié)合的方式，即通過(guò)密集的線性測(cè)線獲得較高空間分辨率的二維深部結(jié)構(gòu)剖面，同時(shí)結(jié)合二維平面上相對(duì)稀疏的流動(dòng)臺(tái)陣觀測(cè)研究較大尺度三維深部結(jié)構(gòu).

臺(tái)站密度在兩條線性測(cè)線和東北側(cè)局部區(qū)域非常密集，而其他區(qū)域相對(duì)稀疏.如前所述，如果統(tǒng)一采用大尺度疊加半徑進(jìn)行CCP成像很容易導(dǎo)致臺(tái)站密集區(qū)成像分辨率的損失，因此將壓縮感知高分辨率接收函數(shù)疊加成像方法應(yīng)用于該套觀測(cè)數(shù)據(jù)，可以驗(yàn)證方法可靠性和有效性.

基于深部結(jié)構(gòu)特點(diǎn)，我們建立如圖5a中黑色三角形所表示的矩形虛擬臺(tái)陣，共71×41個(gè)虛擬臺(tái)站，矩形長(zhǎng)邊、短邊分別沿北西西方向、北北東方向，兩個(gè)方向臺(tái)間距分別約8 km、11 km.

本研究使用了2006—2007年期間30°<震中距<90°、MS>5.0的遠(yuǎn)震事件數(shù)據(jù)計(jì)算獲得的Ps震相接收函數(shù).我們對(duì)低信噪比臺(tái)站數(shù)據(jù)進(jìn)行剔除之后，篩選出臺(tái)站數(shù)大于100個(gè)的遠(yuǎn)震事件100個(gè)(圖6)作為試驗(yàn)數(shù)據(jù)，圖5b展示了所有臺(tái)站累計(jì)使用的接收函數(shù)條數(shù)，臺(tái)站顏色越深表示累計(jì)所使用的該臺(tái)接收函數(shù)條數(shù)越多，最大者達(dá)到100條.

每個(gè)地震事件所篩選出的臺(tái)站存在一定差異，圖5c、5d分別展示了第40個(gè)地震事件和第62個(gè)地震事件篩選出的臺(tái)站分布.接著，以每個(gè)地震事件為處理單元，以該事件所有接收函數(shù)視為不完整數(shù)據(jù)，通過(guò)壓縮感知稀疏促進(jìn)反演獲得虛擬臺(tái)陣接收函數(shù).對(duì)每次事件，本文共進(jìn)行了100次稀疏促進(jìn)反演.我們對(duì)重構(gòu)結(jié)果抽取了2個(gè)比較典型的線性虛擬臺(tái)陣排列(CC′測(cè)線和NN′測(cè)線)進(jìn)行展示.圖7a—b分別展示了第40個(gè)地震事件和第62個(gè)地震事件所抽取的臺(tái)站分布，紅色、藍(lán)色實(shí)心圓分別代表實(shí)際臺(tái)站(距所選取的線性虛擬臺(tái)站排列0.3°以內(nèi))和線性虛擬臺(tái)站排列.其中，第40個(gè)地震事件所選的CC′測(cè)線附近實(shí)際觀測(cè)臺(tái)站比較稀疏，距離0.3°以內(nèi)的臺(tái)站非常少，而第62個(gè)地震事件所選取的NN′測(cè)線附近的實(shí)際臺(tái)站分布則較為密集.圖7c—d分別為第40個(gè)地震事件CC′測(cè)線附近實(shí)測(cè)臺(tái)站數(shù)據(jù)提取的接收函數(shù)以及CC′測(cè)線上重構(gòu)的虛擬臺(tái)站接收函數(shù)，而圖7e—f則分別為第62個(gè)地震事件的NN′測(cè)線附近實(shí)際接收函數(shù)和NN′測(cè)線上的虛擬臺(tái)站接收函數(shù).可以看到第40個(gè)事件重構(gòu)的接收函數(shù)由于附近臺(tái)站非常稀少，導(dǎo)致振幅偏弱.第62個(gè)事件重構(gòu)的接收函數(shù)由于附近實(shí)際觀測(cè)臺(tái)站較多，因此振幅較強(qiáng).

獲得所有地震事件的虛擬臺(tái)陣接收函數(shù)之后，我們依次進(jìn)行接收函數(shù)時(shí)深轉(zhuǎn)換、振幅校正和小尺度疊加半徑疊加成像.其中，成像點(diǎn)網(wǎng)格大小沿著北西西、北北東、深度方向分別為1 km、4 km和0.5 km.我們采用了不同的水平方向疊加箱尺寸對(duì)經(jīng)典CCP成像方法以及本文方法成像結(jié)果進(jìn)行對(duì)比分析.水平方向上，經(jīng)典CCP方法分別采用了5 km、10 km、20 km、30 km疊加半徑，本文方法分別采用了5 km、10 km、20 km疊加半徑，垂直方向上，統(tǒng)一采用了1 km疊加半徑.我們分別展示具有代表性的兩條北西西向測(cè)線，其中，南線(CC′)位于中心區(qū)域，臺(tái)站相對(duì)稀疏，北線(NN′)位于北側(cè)，臺(tái)站相對(duì)密集.

圖7 實(shí)際資料應(yīng)用獲得的虛擬臺(tái)站接收函數(shù)重構(gòu)結(jié)果 (a—b) 第40個(gè)地震事件和第62個(gè)地震事件的臺(tái)站分布.紅色圓(包括實(shí)心和空心)代表145個(gè)實(shí)測(cè)臺(tái)站，綠色圓代表2991個(gè)虛擬臺(tái)站，藍(lán)色實(shí)心圓分別代表沿著CC′或NN′測(cè)線上的虛擬臺(tái)站，紅色實(shí)心圓分別代表距離CC′或NN′測(cè)線30 km以內(nèi)的實(shí)測(cè)臺(tái)站； (c—d) 第40個(gè)地震事件觀測(cè)臺(tái)站(圖(a)中的紅色實(shí)心圓)實(shí)際提取的CC′測(cè)線附近的接收函數(shù)和CC′測(cè)線上虛擬臺(tái)站接收函數(shù)展示； (e—f) 第62個(gè)地震事件觀測(cè)臺(tái)站(圖(b)中的紅色實(shí)心圓)實(shí)際提取的NN′測(cè)線附近的接收函數(shù)和NN′測(cè)線上虛擬臺(tái)站接收函數(shù)展示.Fig.7 The reconstructed receiver functions of virtual stations in the real data application (a—b) Distribution of the stations for the 40th and 62th event. The red (including solid and hollow circles), green circles represent 145 real observation stations and 2991 virtual stations respectively. The blue solid circles are virtual stations located along the CC′ or the NN′ line respectively, and the red solid circles are real observation stations located within 30 km from the CC′ or the NN′ line respectively; (c—d) Receiver functions of the real observation stations of the 40th event near the CC′ line and the reconstructed receiver functions of the virtual stations along the CC′ line; (e—f) Receiver functions of the real observation stations of the 62th event near the NN′ line and the reconstructed receiver functions of the virtual stations along the NN′ line.

圖8 經(jīng)典CCP方法和本文方法成像結(jié)果的南線(CC′)剖面對(duì)比 (a—d) 經(jīng)典CCP方法采用水平疊加半徑=5 km、10 km、20 km、30 km獲得的成像結(jié)果的南線(CC′)剖面； (e—g) 本文方法采用水平疊加半徑=5 km、10 km、20 km獲得的成像結(jié)果的南線(CC′)剖面.Fig.8 The CC′ profile of the result obtained by the CCP method and our proposed method (a—d) The CC′ profile of the result obtained by the CCP method with the horizontal stacking radius=5 km，10 km，20 km，30 km respectively; (e—g) The CC′ profile of the result obtained by our proposed method with the horizontal stacking radius=5 km，10 km， 20 km respectively.

圖9 經(jīng)典CCP方法和本文方法成像結(jié)果的北線(NN′)剖面對(duì)比 (a—d) 經(jīng)典CCP方法采用水平疊加半徑=5 km、10 km、20 km、30 km獲得的成像結(jié)果的北線(NN′)剖面； (e—g) 本文方法采用水平疊加半徑=5 km、10 km、20 km獲得的成像結(jié)果的北線(NN′)剖面.Fig.9 The NN′ profile of the result obtained by the CCP method and our proposed method (a—d) The NN′ profile of the result obtained by the CCP method with the horizontal stacking radius=5 km，10 km，20 km，30 km respectively; (e—g) The NN′ profile of the result obtained by our proposed method with the horizontal stacking radius=5 km，10 km， 20 km respectively.

圖8a—d分別展示了經(jīng)典CCP方法對(duì)實(shí)測(cè)臺(tái)站接收函數(shù)采用疊加半徑5 km、10 km、20 km、30 km疊加成像得到的南線(CC′)剖面結(jié)果.可以看出，5 km和10 km疊加成像時(shí)，水平距離=100 km以西區(qū)域出現(xiàn)較為明顯的成像界面，界面深度從50～40 km不等，在剖面的水平距離=100～200 km處出現(xiàn)明顯的中斷.進(jìn)一步，20 km、30 km疊加成像結(jié)果中，西側(cè)中斷處界面變得更為連續(xù)，但產(chǎn)生明顯的水平階梯狀疊加痕跡，且原本西側(cè)較高分辨率成像部位分辨率產(chǎn)生一定損失.

圖8e—g分別展示了本文方法采用疊加半徑5 km、10 km、20 km對(duì)重構(gòu)的虛擬臺(tái)陣接收函數(shù)進(jìn)行成像得到的南線(CC′)剖面結(jié)果.對(duì)比發(fā)現(xiàn)，在5 km疊加成像結(jié)果中，原有數(shù)據(jù)中包含的橫向不均勻性細(xì)節(jié)基本得到很好刻畫(huà)，在水平距離=100 km以西的區(qū)域，界面顯示出起伏而連續(xù)的變化，且整個(gè)剖面沒(méi)有出現(xiàn)明顯的疊加痕跡.10 km疊加后的界面相比5 km要更為平滑，界面能量更為均衡，同樣沒(méi)有明顯的疊加痕跡，但界面局部起伏形態(tài)有所減緩；20 km疊加結(jié)果界面更為光滑、連續(xù)，但原本局部的橫向不均勻性基本全都消失.因此，我們認(rèn)為對(duì)重構(gòu)的虛擬臺(tái)陣接收函數(shù)進(jìn)行疊前的小尺度疊加成像可以獲得更高分辨率的成像結(jié)果.這樣，可以在保證數(shù)據(jù)較密集部位的成像分辨率的同時(shí)，能夠有效提升臺(tái)站稀疏區(qū)下方的成像質(zhì)量.

北線(NN′剖面)附近有一條北西西向線性密集臺(tái)陣，相比南線臺(tái)間距小、連續(xù)性好.圖9a—d分別展示了對(duì)原始接收函數(shù)按5 km、10 km、20 km、30 km疊加半徑進(jìn)行成像后的結(jié)果.由于臺(tái)站相對(duì)密集，經(jīng)典CCP方法進(jìn)行5 km疊加獲得的北線成像界面連續(xù)性相比南線好很多，但仍然在臺(tái)站稀疏部位出現(xiàn)了成像空白區(qū)(圖9a).這一空白部位在采用10 km以上尺度的疊加半徑之后得到了填補(bǔ)(圖9b)，界面連續(xù)、完整，但階梯狀疊加痕跡較明顯.圖9e—g分別展示了對(duì)虛擬臺(tái)陣接收函數(shù)采用5 km、10 km、20 km疊加半徑得到的剖面.其中，5 km疊加的結(jié)果中界面連續(xù)、分辨率高，局部橫向變化清晰，按10 km、20 km疊加后界面變得更為光滑，但損失了一定分辨率.同樣，北線虛擬臺(tái)陣接收函數(shù)的成像結(jié)果沒(méi)有出現(xiàn)明顯的疊加痕跡.總的來(lái)說(shuō)，由于臺(tái)陣更為密集，北線剖面總體成像質(zhì)量相比南線高，本文方法獲得的成像結(jié)果相比經(jīng)典方法成像分辨率更高，有效克服經(jīng)典方法實(shí)際應(yīng)用中存在的水平分辨率損失問(wèn)題.

3 討論

首先，我們討論壓縮感知接收函數(shù)成像方法的成像結(jié)果是否可靠，尤其是臺(tái)站稀疏區(qū)域.壓縮感知稀疏促進(jìn)算法求解接收函數(shù)一范數(shù)(稀疏性)約束下的實(shí)際接收函數(shù)和重構(gòu)接收函數(shù)之差的二范數(shù)最小化問(wèn)題(見(jiàn)(4)式)，即求得與實(shí)際數(shù)據(jù)擬合度高的稀疏解.數(shù)據(jù)殘差二范數(shù)項(xiàng)，其作用在于約束重建的密集、規(guī)則化虛擬臺(tái)陣接收函數(shù)和實(shí)際接收函數(shù)波形之間的擬合程度.接收函數(shù)的一范數(shù)約束項(xiàng)，用來(lái)約束接收函數(shù)重構(gòu)結(jié)果的稀疏性.這一稀疏性的物理意義可以解釋為，深部速度界面在局部空間范圍內(nèi)變化緩慢、特征穩(wěn)定，因此該局部特征即使只有較少的觀測(cè)數(shù)據(jù)也可以得到充分捕捉，并利用該局部特征重構(gòu)出缺失部位的特征信號(hào).界面的局部特征越顯著，接收函數(shù)的稀疏性越好，重構(gòu)精度越高，成像結(jié)果越可靠.本文的數(shù)值模擬測(cè)試中，145條稀疏臺(tái)站接收函數(shù)無(wú)法對(duì)陡傾角斜面進(jìn)行準(zhǔn)確成像，即便嘗試多個(gè)不同尺度的疊加半徑其成像分辨率的提升效果仍然有限.然而，由于接收函數(shù)具有稀疏特性，我們以壓縮感知理論作為理論基礎(chǔ)，僅使用145個(gè)隨機(jī)分布的實(shí)測(cè)臺(tái)站重構(gòu)獲得了2991個(gè)規(guī)則、密集的虛擬臺(tái)站的接收函數(shù)高精度重構(gòu)結(jié)果，最終對(duì)虛擬臺(tái)站接收函數(shù)進(jìn)行小尺度疊加半徑成像并獲得了高分辨率成像結(jié)果.這表明我們對(duì)地震數(shù)據(jù)的稀疏性假設(shè)以及基于壓縮感知理論進(jìn)行接收函數(shù)重構(gòu)和小尺度半徑疊加成像是合理、有效的.需要說(shuō)明的是，如果在局部區(qū)域界面起伏非常劇烈的同時(shí)臺(tái)站分布也極為稀疏，那么觀測(cè)數(shù)據(jù)不足以約束界面的局部特征.這將導(dǎo)致反演問(wèn)題的解的范圍較大，從而影響接收函數(shù)的重構(gòu)精度以及最終成像結(jié)果的可靠程度.但是，稀疏性是一種相對(duì)概念，臺(tái)站稀疏區(qū)域的觀測(cè)數(shù)據(jù)量仍然可以較好地約束較大尺度的地下結(jié)構(gòu)特征.這是因?yàn)榇蟪叨冉Y(jié)構(gòu)對(duì)應(yīng)的特征仍然是顯著的，因此仍可以從稀疏分布臺(tái)站接收函數(shù)中重構(gòu)出能夠反映大尺度結(jié)構(gòu)特征的虛擬臺(tái)陣接收函數(shù)，也就是說(shuō)，這些臺(tái)站稀疏區(qū)域的成像結(jié)果仍然適用于研究相對(duì)較大尺度的橫向不均勻性.

如果臺(tái)站足夠密集、且等間隔規(guī)則分布，那么經(jīng)典CCP方法足以獲得理想的成像結(jié)果.但一旦臺(tái)站分布不均勻、不夠密集，那么統(tǒng)一采用大尺度疊加窗一定會(huì)帶來(lái)分辨率的損失.本文數(shù)值模擬測(cè)試結(jié)果表明當(dāng)臺(tái)站稀疏區(qū)域地下存在起伏較大的界面時(shí)，采用過(guò)大的疊加窗會(huì)造成分辨率損失.壓縮感知接收函數(shù)成像方法能夠在盡可能保留原有數(shù)據(jù)中的細(xì)節(jié)特征的同時(shí)能夠刻畫(huà)臺(tái)站缺失區(qū)的較大尺度的結(jié)構(gòu)變化特征，極大程度地減少經(jīng)典方法存在的上述分辨率損失問(wèn)題，有效提高成像分辨率.另外，稀疏促進(jìn)反演還可以對(duì)稀疏性差的隨機(jī)干擾噪聲有很好的壓制作用，理論上能夠進(jìn)一步提升成像質(zhì)量.因此，壓縮感知接收函數(shù)成像方法在深部結(jié)構(gòu)研究中具有很好的實(shí)用價(jià)值.

其次，我們討論是否有必要在虛擬臺(tái)陣接收函數(shù)成像時(shí)進(jìn)行振幅校正.圖10a展示了對(duì)某個(gè)地震事件的接收函數(shù)進(jìn)行成像時(shí)的振幅校正值在水平面上的分布.圖中，紫色三角形為該地震事件對(duì)應(yīng)的臺(tái)站分布.對(duì)該事件接收函數(shù)進(jìn)行時(shí)深轉(zhuǎn)換后，都按照?qǐng)D10a所示的振幅校正值進(jìn)行振幅校正.成像點(diǎn)與其最近臺(tái)站的水平距離越近，其振幅校正值越接近1，校正量越小，反之，則校正量越大.圖10b—c分別展示了本文方法在振幅校正前、后成像結(jié)果的南線(CC′)剖面.圖10b中剖面右側(cè)對(duì)應(yīng)臺(tái)站稀疏區(qū)，成像能量較弱，在振幅校正后，振幅得到補(bǔ)償(圖10c)，成像界面能量更為均衡.因此，當(dāng)臺(tái)站密度不均勻時(shí)，我們認(rèn)為振幅校正有助于提高虛擬臺(tái)陣接收函數(shù)成像質(zhì)量.

本文的數(shù)值測(cè)試和實(shí)際資料應(yīng)用均采用水平層狀模型進(jìn)行成像，未來(lái)研究中可以進(jìn)一步采用各臺(tái)站下方的層速度模型，以獲得更為可靠的成像結(jié)果.另外，基于波動(dòng)方程的接收函數(shù)疊前偏移成像方法理論上相比本文所采用的基于射線理論的經(jīng)典CCP疊前成像方法對(duì)界面起伏劇烈的復(fù)雜構(gòu)造具有更高的分辨率，但三維波動(dòng)方程疊前成像的計(jì)算成本非常高，其實(shí)用性有待進(jìn)一步攻克.對(duì)于本研究中的應(yīng)用實(shí)例，界面起伏變化相對(duì)緩慢，采用射線類疊前成像方法足以得到可靠的高分辨率成像結(jié)果.

4 結(jié)論

本文將壓縮感知理論應(yīng)用于深部結(jié)構(gòu)成像方法研究，基于深部速度間斷面在局部空間范圍內(nèi)變化緩慢、特征穩(wěn)定從而接收函數(shù)具有稀疏特性這一假定，通過(guò)壓縮感知重構(gòu)技術(shù)獲得密集、規(guī)則化虛擬臺(tái)陣的接收函數(shù)，再通過(guò)成像點(diǎn)振幅校正和小尺度疊加半徑疊加成像，獲得高分辨率的深部轉(zhuǎn)換界面成像剖面.本文方法可以極大程度地減少經(jīng)典CCP成像方法在實(shí)際應(yīng)用中為了兼顧臺(tái)站稀疏區(qū)域而統(tǒng)一采用大尺度疊加半徑所引起的臺(tái)站密集區(qū)成像分辨率損失以及成像界面上階梯狀的疊加痕跡，能夠在保證臺(tái)站密集區(qū)成像結(jié)果的高分辨率的同時(shí)提升臺(tái)站稀疏區(qū)的成像質(zhì)量.數(shù)值模擬測(cè)試驗(yàn)證了壓縮感知接收函數(shù)成像方法可以充分利用接收函數(shù)的稀疏性，從極少的觀測(cè)數(shù)據(jù)中高精度重構(gòu)出規(guī)則、密集的虛擬臺(tái)站接收函數(shù)，并進(jìn)一步通過(guò)小尺度疊加半徑實(shí)現(xiàn)對(duì)地下起伏界面的高分辨率成像.我們將本文方法應(yīng)用于中國(guó)地震科學(xué)探測(cè)臺(tái)陣(ChinArray)華北地區(qū)觀測(cè)數(shù)據(jù)，獲得了華北克拉通深部結(jié)構(gòu)成像剖面.實(shí)際應(yīng)用結(jié)果表明，利用壓縮感知高分辨率接收函數(shù)疊加成像可以有效提升成像結(jié)果的空間分辨率，有效提升成像質(zhì)量.基于壓縮感知重構(gòu)技術(shù)的接收函數(shù)成像方法在地球深部精細(xì)結(jié)構(gòu)研究中具有廣闊的應(yīng)用前景.

圖10 虛擬臺(tái)陣疊加成像振幅校正量分布及校正前后的成像結(jié)果對(duì)比(a) 單個(gè)地震事件成像過(guò)程中虛擬臺(tái)陣振幅校正值平面分布.圖中紫色三角形代表實(shí)測(cè)臺(tái)站； (b—c) 振幅校正前和校正后的 虛擬臺(tái)陣接收函數(shù)成像結(jié)果的南線(CC′)剖面(5 km疊加半徑).Fig.10 Distribution of the amplitude correction value for the virtual stations and comparison of the imaging results before and after the correction(a) Distribution of the amplitude correction value for the virtual stations in single event imaging process. Purple triangles represent the real observed stations; (b—c) The CC′ profile (with the stacking radius = 5 km) of the imaging result using the receiver functions of virtual stations before and after the amplitude correction.

致謝本文實(shí)際資料使用中國(guó)地震科學(xué)探測(cè)臺(tái)陣觀測(cè)數(shù)據(jù)，遠(yuǎn)震波場(chǎng)的正演模擬通過(guò)Specfem3D開(kāi)源軟件中的SEM-FK模塊在國(guó)家超級(jí)計(jì)算天津中心計(jì)算實(shí)現(xiàn).感謝三位審稿人對(duì)本文提出的寶貴建議.本研究開(kāi)展數(shù)值模擬測(cè)試期間，與多倫多大學(xué)何彬博士、應(yīng)急管理部國(guó)家自然災(zāi)害防治研究院劉少林研究員對(duì)SEM-FK方法進(jìn)行了討論，在此表示感謝.