压缩感知高分辨率接收函数叠加成像及其应用

白兰淑， 吴庆举， 张瑞青

中国地震局地球物理研究所， 北京 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方法进行了讨论，在此表示感谢.