壓縮感知高分辨率接收函數疊加成像及其應用

白蘭淑， 吳慶舉， 張瑞青

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

0 引言

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

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

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

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

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

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

1 方法

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

1.1 壓縮感知稀疏促進算法重構接收函數

壓縮感知理論認為稀疏信號可以通過稀疏促進算法從少量觀測數據中重構出高精度的原始信號(Donoho，2006).因此，接收函數是否具有稀疏性是壓縮感知理論能否成功應用于接收函數重構和高分辨率成像的重要前提條件.

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

下面，我們按照壓縮感知理論建立對原本不規則、不均勻分布地震臺陣接收函數進行重構的基本框架.

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

f=Ψu,

‖u‖0=k,

(1)

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

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

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

(2)

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

對于不規則分布的地震臺陣而言，可以認為其觀測數據在二維水平空間維度上是近隨機的.實際布設中，也存在較多的地震臺陣采用近等間距的規則布設.關于地震觀測點規則或不規則分布時的地震數據重構效果，已有大量通過模擬數據或實際數據開展的研究和分析(尤其是勘探地震學領域).研究表明只要地震數據具有稀疏性，均可以高概率獲得高精度重構結果(Herrmann and Hennenfent，2008；唐剛，2010；白蘭淑等，2014；黃小剛等，2014；Bai et al.，2017，2020).因此，我們假定無論是相對等間距、規則布設，還是不等間距、隨機布設，都可以實現虛擬臺陣接收函數的成功重構.

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

(3)

求解公式(3)中的稀疏促進問題是一個經典的NP-hard非凸優化問題，對于稀疏信號而言可以通過用L1范數代替L0范數變成凸優化問題，并進一步轉化為如下非約束化的子問題：

(4)

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

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

(5)

1.2 疊前振幅校正

通過稀疏促進算法重構虛擬臺陣接收函數時，由于虛擬臺陣的臺間距比實際布設臺間距小很多，重構的虛擬臺站接收函數振幅偏弱，尤其是在臺站極為稀疏的區域.如果直接按照經典方法進行簡單的疊加平均操作，會導致成像界面能量不均衡.

針對上述問題，我們提出一種疊前振幅校正公式，具體表達式如下：

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

(6)

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

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

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

疊前振幅補償后，我們采用小尺度疊加半徑進行傳統CCP疊加成像(ccp1.0. tar package 2006. https:∥www.eas.slu.edu/People/LZhu/downloads/ccp1.0.tar)，獲得高分辨率地下結構成像結果.疊加半徑大小通常需要綜合實際臺陣分布、虛擬臺陣網格間距等因素進行選取，通常可采用虛擬臺陣臺間距相當的尺度.

(7)

式中，divsta為虛擬臺站接收函數時深轉換之后的(單條射線對應的)深度域振幅，xima代表地下第ima個成像點的坐標，ievt、ivsta分別代表地震事件和虛擬臺站編號，Nray為以xima為中心的疊加箱內總射線條數，Image為振幅校正后疊加得到的像.

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

本文結合壓縮感知理論和經典CCP接收函數疊加成像的方法具體實現流程如下:

(1)提取遠震事件接收函數，篩選高信噪比接收函數臺站數較多的遠震事件接收函數；

(2)根據實際地震臺陣分布，在矩形區域內建立密集、規則化虛擬地震臺陣；

(3)針對每一個地震事件，利用篩選出的接收函數，通過表1所示的壓縮感知稀疏促進反演算法重構虛擬臺陣接收函數；當篩選出的地震事件共Nevt個，虛擬臺陣共Nvsta個，則所有地震事件重構的總接收函數共Nevt×Nvsta條.

(4)綜合區域構造背景、實際觀測臺陣分布、虛擬臺陣分布等因素，建立三維地下成像點網格；

(5)設定合適的小尺度長方體疊加箱，對所有篩選的遠震事件依次進行(6)—(7)的操作；

(6)在水平層狀界面假設下，根據參考速度模型和遠震Ps轉換波對應的射線參數計算每個臺站對應的轉換波射線路徑.對所有臺站接收函數時間序列從零時刻開始，從地表向下沿著轉換波射線路徑放在射線所落入的成像點上，即時深轉換.

(7)以每個地震事件視為處理單元，根據公式(6)式對虛擬臺陣接收函數依次進行成像點振幅校正；

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

(9)獲得三維地下結構成像剖面.

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

2 方法測試與實際應用

2.1 數值模擬測試

我們設計了如圖1所示的2.5維模型，該模型介質參數沿東西向呈不均勻變化，沿南北方向保持不變. 模型中存在一個界面，其形態呈中段傾斜、兩側水平，界面深度最左端達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.遠震波場從模型東側底部以平面波的形式入射，入射角20°、后方位角90°.地面接收臺站共71×41=2991個，等間距布設，東西、南北向臺間距分別為4.29 km、5 km，臺站布設區域東西、南北向跨度分別達300 km、200 km.

圖1 數值測試中采用的速度模型Fig.1 Velocity model used in the synthetic test

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

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

圖2 正演模擬獲得的CC′測線上的波場記錄(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 壓縮感知稀疏促進算法獲得的接收函數重構結果(a) 臺站平面分布圖.灰色空心圓代表2991個規則、密集臺站，用于獲得規則、密集的地震波場記錄，虛擬臺站位置與這些臺站位置一致.紅色圓圈(包括實心和空心)代表我們從2991個地震臺站中隨機抽選出的145個臺站，占比約4.8%，我們視其為“實測臺站”，這些實測臺站接收記錄用作實際觀測數據.藍色實心圓代表CC′測線上的71個虛擬臺站，紅色實心圓代表該條測線10 km范圍內的實測臺站； (b) CC′測線附近實測臺站(圖3a中紅色實心圓)的接收函數； (c) CC′測線虛擬臺站處重構的接收函數； (d) CC′ 測線虛擬臺站處參考接收函數.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 壓縮感知接收函數成像方法的數值模擬測試結果及對比(a) 經典方法利用2991個規則、密集臺站真實接收函數獲得的結果(疊加半徑=2 km)； (b) 壓縮感知接收函數成像方法利用145個臺站數據重構的2991個虛擬臺站接收函數進行成像獲得的結果(疊加半徑=2 km)； (c—e) 經典方法利用145個臺站(圖3a中紅色 圓圈所示)接收函數、并分別采用疊加半徑=2 km、5 km、10 km獲得的結果.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 中國地震局地球物理研究所華北科學臺陣與密集、規則化虛擬臺陣分布(a) 中國地震局地球物理研究所科學臺陣與虛擬臺陣分布，分別用紫色、黑色三角形表示； (b) 各臺站使用的接收函數條數.顏色越深， 表示所使用接收函數越多； (c—d) 第40個和第62個地震事件所使用的接收函數對應臺站分布.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 本文選用的遠震事件震源分布Fig.6 Distribution of teleseismic events chosen in the real data application

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

2.2 實際應用實例

本研究進一步將基于壓縮感知重構技術的接收函數疊加成像方法應用于華北克拉通觀測地震數據，獲得了華北克拉通深部結構成像剖面.需要說明的是，我們僅討論本文研究方法對地下界面成像分辨率和成像質量的提升效果，不分析深部結構特征及其所指示的動力學問題.

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

臺站密度在兩條線性測線和東北側局部區域非常密集，而其他區域相對稀疏.如前所述，如果統一采用大尺度疊加半徑進行CCP成像很容易導致臺站密集區成像分辨率的損失，因此將壓縮感知高分辨率接收函數疊加成像方法應用于該套觀測數據，可以驗證方法可靠性和有效性.

基于深部結構特點，我們建立如圖5a中黑色三角形所表示的矩形虛擬臺陣，共71×41個虛擬臺站，矩形長邊、短邊分別沿北西西方向、北北東方向，兩個方向臺間距分別約8 km、11 km.

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

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

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

圖7 實際資料應用獲得的虛擬臺站接收函數重構結果 (a—b) 第40個地震事件和第62個地震事件的臺站分布.紅色圓(包括實心和空心)代表145個實測臺站，綠色圓代表2991個虛擬臺站，藍色實心圓分別代表沿著CC′或NN′測線上的虛擬臺站，紅色實心圓分別代表距離CC′或NN′測線30 km以內的實測臺站； (c—d) 第40個地震事件觀測臺站(圖(a)中的紅色實心圓)實際提取的CC′測線附近的接收函數和CC′測線上虛擬臺站接收函數展示； (e—f) 第62個地震事件觀測臺站(圖(b)中的紅色實心圓)實際提取的NN′測線附近的接收函數和NN′測線上虛擬臺站接收函數展示.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 經典CCP方法和本文方法成像結果的南線(CC′)剖面對比 (a—d) 經典CCP方法采用水平疊加半徑=5 km、10 km、20 km、30 km獲得的成像結果的南線(CC′)剖面； (e—g) 本文方法采用水平疊加半徑=5 km、10 km、20 km獲得的成像結果的南線(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 經典CCP方法和本文方法成像結果的北線(NN′)剖面對比 (a—d) 經典CCP方法采用水平疊加半徑=5 km、10 km、20 km、30 km獲得的成像結果的北線(NN′)剖面； (e—g) 本文方法采用水平疊加半徑=5 km、10 km、20 km獲得的成像結果的北線(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分別展示了經典CCP方法對實測臺站接收函數采用疊加半徑5 km、10 km、20 km、30 km疊加成像得到的南線(CC′)剖面結果.可以看出，5 km和10 km疊加成像時，水平距離=100 km以西區域出現較為明顯的成像界面，界面深度從50～40 km不等，在剖面的水平距離=100～200 km處出現明顯的中斷.進一步，20 km、30 km疊加成像結果中，西側中斷處界面變得更為連續，但產生明顯的水平階梯狀疊加痕跡，且原本西側較高分辨率成像部位分辨率產生一定損失.

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

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

3 討論

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

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

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

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

4 結論

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

圖10 虛擬臺陣疊加成像振幅校正量分布及校正前后的成像結果對比(a) 單個地震事件成像過程中虛擬臺陣振幅校正值平面分布.圖中紫色三角形代表實測臺站； (b—c) 振幅校正前和校正后的 虛擬臺陣接收函數成像結果的南線(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.

致謝本文實際資料使用中國地震科學探測臺陣觀測數據，遠震波場的正演模擬通過Specfem3D開源軟件中的SEM-FK模塊在國家超級計算天津中心計算實現.感謝三位審稿人對本文提出的寶貴建議.本研究開展數值模擬測試期間，與多倫多大學何彬博士、應急管理部國家自然災害防治研究院劉少林研究員對SEM-FK方法進行了討論，在此表示感謝.