接收函数方法(Langston，1979)已成为地壳上地幔地震成像的一个重要方法．但是，迄今为止，接收函数的解释仍主要基于水平分层模型(Owens et al.，1984；Duecker and Sheehan，1997；Zhu and Kanamori，2000)．例如，对于接收函数反演来说，为了适应一维水平分层模型的限制，人们不得不通过不同方位接收函数的叠加来压制数据中横向非均匀地壳介质引发的散射(Liu et al．，2014)．另一方面，近年来，各向异性接收函数的研究日益受到人们的重视(田宝峰等，2008；房立华和吴建平，2009；Schulte-Pelkum and Mahan，2014)，但地壳横向非均匀介质引起的散射给接收函数各向异性参数提取带来了困难，而低通滤波将会降低各向异性参数的分辨能力．这并非是人们所希望的.

因此，对基于地壳水平分层模型的接收函数研究来说，在提取接收函数的过程中，环境噪声已得到很大程度的压制，而影响信息提取的主要噪声来自于介质横向非均匀引起的散射．与环境噪声不同，地壳横向非均匀介质引起的散射属于自激式的次生波场，它通常与试图分析研究的信号具有相近的频带，且在空间取向上具有较大的随机性．但是，理论合成数据的检验表明，各向异性介质接收函数波场在空间上具有很好的连续性，或者称为可追踪性(Nagaya et al．，2008；房立华和吴建平，2009)．本文将基于这一认识，研究压制接收函数中的横向非均匀散射方法，这对推进接收函数研究将具有重要的实际意义．

曲波变换(curvelet transform)是近年来迅速发展的多尺度、多方向高维信号分析工具(Candès and Donoho，2000，2004)．相较于经典的二维数据去噪方法，基于曲波变换的去噪，对信号的损伤最小(张军华等，2005)．在方向性特征识别上，曲波变换优于小波变换(wavelet transform)．它尤其适用于地震同相轴检测与去噪．迄今为止，曲波变换已经被广泛地应用于地震勘探数据的多次波去除(Herrmann and Verschuur，2004；张素芳等，2006)，噪声压制(Neelamani et al．，2008；彭才等，2008)，偏移成像(Douma and De Hoop，2007)以及地震环境噪声的格林函数提取(Stehly et al．，2011，2015)．曲波变换所具有的特性为我们压制接收函数中横向非均匀地壳介质引起的散射提供了新思路．

但是，直接将曲波变换应用于接收函数是困难的，因为曲波变换要求数据在各个维度上必须满足等间距采样．一方面受制于地震全球分布的条带性，另一方面受具体环境条件对地震台站(特别是流动台站)布设的限制，接收函数在空间方向上的采样通常是非等间距的，甚至可能存在大范围(例如，在方位空间)数据缺失．因此，对于曲波变换去噪处理，我们需要对接收函数波场进行正则化(均匀空间采样)．利用距离加权叠加的方法，Chai等(2015)实现了对空间分布相对密集且接近均匀的接收函数波场的正则化．Zhang和Zheng(2015)提出用三次样条插值提高稀疏台站接收函数地震成像的分辨率．作为一种探索，本文将引入由Candès等(2006b)和Donoho(2006)提出的压缩感知(compressive sensing)理论，对接收函数波场进行重建，进而形成可用于曲波变换处理的数据，实现接收函数数据集的二维去噪．为便于理解，本文的讨论将从简要介绍曲波变换和压缩感知的基本概念开始．

曲波可以认为是一类具有方向性的特殊小波，它提供对曲线型边缘的最优稀疏表示(Candès and Donoho，2004)．曲波在空间域和频率域内均满足抛物尺度关系，即宽度约等于长度的平方．在空间域内，曲波的能量在其垂向上呈现振荡衰减，形态接近平面波．在频率域内，曲波的支撑区域是一个围绕原点的楔形区域，它属于紧支撑．为了便于直观理解，我们利用Curvelab工具包绘制了曲波在空间域和频率域内的对应关系，如图 1所示．需要说明的是，1)为了便于观察，我们对空间域内的曲波进行了平移，同时，因为二维信号的频谱关于原点对称，我们仅绘制了一半的支撑域；2)曲波在空间域和频率域的相互正交特征取决于Fourier变换的时频缩放特性．

图 1

Fig. 1

图 1 曲波在空间域和频率域的对应关系 (a) 空间域的4个尺度、方向和位置不同的曲波； (b) 频率域中相应曲波的支撑区域． (a)和(b)中曲波的数字编号相互对应． Fig. 1 Curvelets in the spatial and frequency domain (a) Four curvelets with different scales, directions and positions in the spatial domain； (b) Their corresponding support regions in the frequency domain．Curvelets have been labelled with the number in both domains．

根据Candès等(2006a)，曲波函数定义为

如前所述，曲波变换基为紧框架，其重构公式具有以下形式：

为此，假定由观测数据得到的二维(空间)接收函数d符合线性矩阵方程:

传统的信号采样需要遵循Nyquist-Shannon采样定理，即有限带宽的信号采样频率必须达到其带宽的2倍以上才能够被精确重建．但是，近年来出现的压缩感知理论突破了经典采样定律的限制. 压缩感知理论提出后迅速受到地震勘探领域的关注(Herrmann and Hennenfent，2008；Naghizadeh and Sacchi，2010)．在天然地震领域，压缩感知理论也已被用于大震的震源破裂模式研究(Yao et al．，2011，2013)和天然地震波场的逆时偏移(Shang，2014)．

压缩感知理论认为，当信号是稀疏的或在某个变换域内是稀疏的，利用一个与变换基不相关的观测矩阵(或称为感知矩阵)可将高维信号投影到低维空间，并通过求解一个优化问题，以高概率实现由少量投影重建原始信号(Candès and Wakin，2008).这意味着从远低于Nyquist-Shannon采样频率的样本中也可以重建完整信号．实际上，稀疏性条件意味着信号所携带的信息存在某种冗余度，即信号的自由度远小于信号的长度．如果能够通过选择合适的变换基，使信号变换系数大部分为零或小到足以忽略，则可对该信号进行数据压缩．

原则上，低于Nyquist频率的采样操作将导致信号失真，即等间距的欠采样会引发频谱的重叠效应，而随机欠采样则等效于频率域内的随机噪声．当观测矩阵与变换基不相关时，采样不完整所引发的噪声与原始信号在变换域内可被分离．一般地，随机欠采样与常用的变换基不相关．这样，信号重建问题转变为频率域内的去噪问题．

需要关注的是，观测矩阵相应于对连续信号的采样操作．接收函数在时间上的采样满足采样定理，但在空间上采样是不均匀的. 若我们用d表示实际观测的接收函数矩阵，且矩阵的行代表时间，列代表空间位置，则有

根据Naghizadeh和Sacchi(2010)，我们选择曲波作为变换基，则有

若令A=GC^HW，待求曲波系数c可以通过最小化目标函数

为检验方法的可行性，我们首先考虑单个台站记录的不同方位接收函数．为此，利用作者独立开发的程序，我们计算了合成的各向异性接收函数，并将它作为数值检验的依据．需要说明的是，我们依据的算法与Levin和Park(1997)的方法没有原则上区别．

计算各向异性理论接收函数所依据的地壳模型参数如表 1所示．由表 1可知，为简单起见，我们仅考虑平面波入射单层各向异性地壳的情况．计算中采用的水平慢度为0.0618 s·km^-1．大体上，这相当于震中距为60°的情况．对于合成的接收函数，我们加入高斯白噪声，用来模拟包含噪声的观测数据．加入噪声后的接收函数信噪比为0 dB．

表 1(Table 1)

表 1 数值检验所用的各向异性模型 Table 1 Anisotropic model parameters for the numerical test 深度 (km) α (km·s-1) β (km·s-1) ρ (g·cm-3) B E θ(°) φ(°) 30 6.7 3.8 2.7 10% 10% 45 110 ∞ 7.8 4.5 3.3 - - - - 注：α,β,ρ分别为P波速度、S波速度和密度．B和E分别为P波和S波各向异性强度，θ和φ分别为各向异性的倾角和方位角．

表 1 数值检验所用的各向异性模型 Table 1 Anisotropic model parameters for the numerical test

我们考虑接收函数在空间方位上均匀采样的情况．图 2给出了相应的理论接收函数切向分量及其随反方位角(back azimuth，0°～360°，间隔为1°)的变化．图 2表明，在时间―空间域内，接收函数切向分量的同相轴有很好的可追踪性(图 2a)，而在相应的频率域内，这表现为主要能量仅集中在狭窄范围内(图 2b)，这意味着，若观测数据中存在横向非均匀散射，经曲波变换就可以简单地通过阈值去噪．

图 2

Fig. 2

图 2 理论接收函数切向分量及其频谱 (a) 在时-空域内的接收函数切向分量， 红色和蓝色分别代表正振幅和负振幅； (b) 二维Fourier振幅谱， 纵轴和横轴分别对应时间和反方位角． Fig. 2 The transverse component of synthetic anisotropic receiver functions in the time and frequency domain (a) The transverse component of receiver functions in the time domain．Red color denotes positive amplitude and blue color negative amplitude； (b) The normalized Fourier spectrum of transverse component of receiver functions． The vertical and horizontal axes are relevant to the time and back azimuth，respectively．

图 3给出了对接收函数(图 2a)加入随机白噪声后(图 3a)利用本文阈值去噪方法得到的结果(图 3b)．在进行曲波阈值去噪时，我们仅保留了大于阈值(阈值函数在最高尺度取值为4σ，而在其他尺度取值为3σ，其中σ为加入噪声的标准差)且接近垂直方向的曲波系数(Stehly et al．，2011)．比较图 2a和图 3b可知，经过曲波阈值去噪后，原有信号可以正确恢复．

图 3

Fig. 3

图 3 接收函数(切向分量)的曲波阈值去噪 (a) 加入随机噪声后的接收函数； (b) 去噪后的接收函数． Fig. 3 Denoised receiver functions (transverse component) using curvelet transform (a) Noisy receiver functions； (b) Denoised receiver functions．

如前所述，与通常的人工地震勘探不同，接收函数沿不同方位的采样受控于震源的分布．这意味着我们很难将曲波阈值去噪技术直接用于观测数据．为此，我们通过压缩感知对进行接收函数波场进行重建，以便得到在空间上均匀采样的接收函数波场．为了模拟这种情况，如图 4a所示，通过对图 2a的随机抽样，我们构筑了方位非均匀分布的接收函数波场．图 4b给出了利用本文方法重建的接收函数波场(切向分量)．重建波场的信噪比为36.3 dB，表明接收函数的振幅、极性和相位均可以正确恢复(Herrmann and Hennenfent，2008)．

图 4

Fig. 4

图 4 方位非均匀分布的合成接收函数(切向分量)及其波场重建 (a) 由随机抽样得到的方位非均匀分布的接收函数波场； (b) 重建的接收函数波场．其反方位角的间隔为1°． Fig. 4 Uneven azimuthal wavefield of the synthetic receiver functions (transverse component) and its reconstruction (a) Uneven azimuthal wavefield after being randomly sampled； (b) Wavefield reconstruction．The interval of back azimuth is 1°．

对于单台多事件观测来说，数据中缺失部分的重建精度主要取决于地震事件的分布．通常，固定台站由于记录时间长，较容易获得方位分布密集的地震事件．因此，本文首先选取IRIS全球地震台网在澳大利亚东北部地区的固定台CTAO(-20.0882°N，146.2545°E，台网编号IU). 它所处的地理位置有利于获取方位较完整的地震事件覆盖．所用的地震为该台站记录的151个远震事件(自1999年2月到2015年7月，震中距30°≤Δ≤90°，M_b>5.5)，接收函数的提取采用了基于单台的时间域迭代反褶积方法(Ligorría and Ammon，1999)．

图 5给出了台站CTAO的151个远震接收函数的二维时间-方位Fourier谱．图 5显示了分布在整个空间方位的较弱噪声．比较图 5b和图 2b可知，与横向均匀分层模型不同，由于实际接收区地壳介质的横向非均匀性，接收函数的切向分量存在着不可忽视的横向非均匀散射．图 5a表明，这种横向非均匀散射同样存在于接收函数的径向分量，尽管其散射强度相对较弱，它应是非均匀空间采样导致的结果．

图 5

Fig. 5

图 5 台站CTAO不同方位接收函数的二维Fourier谱 (a) 径向分量(R)； (b) 切向分量(T)． Fig. 5 2D Fourier spectrum of the receiver functions with back azimuths at Station CTAO (a) Radial component (R)； (b) Transverse component (T)．

为了客观评价波场重建在实际数据应用中的效果，我们首先将提取的151组接收函数(径向与切向分量)按照反方位角1°的间隔叠加获得123组接收函数，叠加后的事件分布如图 6c所示．我们以类似jitter欠采样的方式(Hennenfent and Herrmann，2008)从123组接收函数中随机抽取出35组来模拟方位缺失的数据(如图 6d中蓝色和红色圆圈所示). 利用剩余的88组接收函数，以反方位角1°为间隔，进行波场重建和去噪．最后，我们从重建后的接收函数中选择出20组数据(如图 6d中蓝色圆圈所示)，并与实际观测接收函数进行波形对比(如图 6a和6b所示)．

图 6

Fig. 6

图 6 台站CTAO的观测数据与重建数据的波形对比 (a) 径向分量的波形对比； (b) 切向分量的波形对比，其中，黑色实线代表实际观测数据，红色实线代表波场重建数据；(c) 按反方位角1°间隔叠加后的台站CTAO的远震事件分布，红色圆圈代表远震事件，绿色五角星代表台站位置； (d) 从(c)中随机抽取 出的远震事件分布，蓝色圆圈代表(a)和(b)中进行对比的远震事件． Fig. 6 Waveform comparison between the observed and reconstructed data from station CTAO (a) Radial component； (b) Transverse component， black solid traces：the observed data， red solid traces：the reconstructed data；(c) Teleseismic event distribution of station CTAO after azimuthal binning across interval of 1°， red circle： event， green star： station； (d) Randomly extracted events from (c). Blue circles are selected for waveform comparison in (a) and (b)．

观察可知，波场重建很好地恢复了缺失部分的数据．与此同时，波场重建校正了观测数据中的一些振幅和相位的异常，重建后的接收函数随反方位角振幅和相位的变化更为连续和合理，这是因为曲波阈值去噪与波场重建是同时进行的．因此，对于单台多事件观测来说，本文提出的曲波变换去噪方法不仅能够有效压制接收函数中的随机散射噪声，而且在地震事件分布较好的情况下或地震事件相对密集的一些区域内，可以对方位缺失的数据进行重构．

图 7给出了CTAO台站123组接收函数的观测数据和曲波去噪后数据按反方位角10°的间隔分组叠加的结果．对图 7我们可有以下观察：1)利用本文给出的曲波阈值去噪方法，我们不但得到了完整(方位)的接收函数波场，而且由于横向散射波场得到压制，经过曲波去噪的接收函数(径向和切向)相较于原始数据显示了更为清晰、合理的震相追踪效果，表明去噪并未削弱有效信息的能量；2)对比去噪处理前后全方位接收函数叠加的结果(sum)，它们的波形没有明显的不同；3)对方位缺失部分预测的接收函数振幅和相位的追踪与其他部分是协调的，这为后续研究提供更多可用数据．

图 7

Fig. 7

图 7 固定台站CTAO实际接收函数的曲波去噪 (a) 去噪前的径向分量； (b) 去噪前的切向分量； (c) 去噪后的径向分量； (d) 去噪后的切向分量． 各方位的接收函数为 以10°为间隔进行方位分组叠加的结果．‘sum’表示全方位接收函数的叠加． Fig. 7 Noise attenuation of real receiver functions with curvelet transform at permanent station CTAO (a) Radial component before denoise；(b) Transverse component before denoise；(c) Radial component after denoise；(d) Transverse component after denoise．Each trace is the result through azimuthal binning across interval of 10°．‘sum’denotes the summation over all of back azimuths．

相对于固定台站，流动地震台站由于观测时间较短(除地理位置的因素之外)，通常很难获得方位覆盖密集的事件分布．同时，对于中国大陆布设的台站来说，往往存在如图 8a所示的地震事件在西北方向的严重缺失(原因是地中海地震带的地震活动性相对较弱，M_b>5.0地震事件较少)．为此，需要探讨方位严重缺失对波场重建带来的影响．

图 8

Fig. 8

图 8 按实际地震方位抽样的理论接收函数波场重建 (a) 台站S05记录的地震事件分布； (b) 按实际地震方位采样的接收函数波场(切向分量)； (c) 重建的接收函数波场(切向分量)； (d)重建合成接收函数的波场误差．SB：四川盆地； SG：松潘地块； CD：川滇地块； LMS：龙门山断裂带； XSH：鲜水河断裂带； 黄色五 角星：台阵中心位置； 红色三角：观测台站； 红色圆圈：地震事件． Fig. 8 Wavefield reconstruction of synthetic receiver functions with real event distribution (a) Teleseismic event distribution at Station S05； (b) Receiver functions with real event distribution (transverse component)； (c) Reconstructed wavefield of the receiver functions； (d) Errors of reconstructed wavefield．SB：Sichan basin； SG：Songpan terrain； CD： Chuandian fragment； LMS： Longmen Shan fault belt； XSH：Xianshuihe fault belt； Yellow star：array center； Red triangle： station； Red circle：event．

实际上，实际数据方位的缺失导致对实际数据按反方位角的采样难以满足随机性. 为此，作为方法检验，我们按照实际数据的方位分布对理论接收函数波场进行抽样并重建．我们选用了川西台阵中的一个台站(S05)，图 8a给出了它的地理位置及相应的震中分布.有关川西台阵和台站详细的信息可参考相关文章(刘启元等，2008)．

图 8a表明，台站S05的接收函数方位除了空间上的不均匀分布，在西北方向反方位角200°~300°的区间还存在明显的方位缺失．由于其方位缺失超过了曲波长度，它导致缺失的数据难以精确重建(Naghizadeh and Sacchi，2010)．作为对压缩感知理论在方位缺失情况下波场重建能力的检验，图 8d给出了重建的理论接收函数与原有的合成接收函数波场(图 2a)的差异．图 8d表明，利用本文方法可以正确重建方位密集区的波场，但对于方位严重缺失的区域，重建的接收函数存在大约6%～8%的振幅误差，并导致重建波场的信噪比下降(10.4 dB)．

为了进一步检验本文曲波阈值去噪方法能否用于流动观测数据，并对结果进行评价．我们对台站S05记录的280个(M_b>5.0)接收函数进行曲波去噪，接收函数的提取采用了基于台阵观测的三分量接收函数的多道最大或然性反褶积方法(刘启元和Kind，2004)，其反方位角(back azimuth)范围为30°～311°，震中距范围为30° ≤Δ≤ 90°．

图 9给出了台站S05接收函数的观测数据和重建数据按反方位角10°的间隔分组叠加的结果．图 9表明：1)经过曲波去噪的接收函数(径向和切向)相较于原始数据显示了更为清晰、合理的震相追踪效果；2)对于反方位角200°～300°的严重缺失部分(已超过曲波长度)，接收函数的振幅不能精确恢复，但这并不破坏该区域之外的接收函数；3)对比去噪处理前后全方位接收函数叠加的结果(sum)，它们的波形没有明显不同；4)对方位严重缺失部分预测的接收函数相位追踪与其他部分是协调的，表明其相位信息是有效的．

图 9

Fig. 9

图 9 台站S05实际接收函数的曲波去噪 (a) 去噪前的径向分量；(b) 去噪前的切向分量；(c) 去噪后的径向分量；(d) 去噪后的切向分量．各方位的接收函数为 以10°为间隔进行方位分组叠加的结果．‘sum’表示全方位接收函数的叠加． Fig. 9 Noise attenuation of real receiver functions with curvelet transform at Station S05 (a) Radial component before denoise；(b) Transverse component before denoise；(c) Radial component after denoise；(d) Transverse component after denoise．Each trace is the result through azimuthal binning across interval of 10°．‘sum’denotes the summation over all of back azimuths．

下面，我们进一步考虑多台阵列观测的情况．对于二维线性观测剖面来说，接收函数的曲波阈值去噪在一定程度上类似于地震勘探的情况(Herrmann and Hennenfent，2008；彭才等，2008)．两者的不同在于，用于天然地震观测的线性台阵难以保证空间上的等间隔采样．这对于本文考虑的曲波阈值去噪来说并没有原则上的困难．与单个台站的情况不同，这种观测方式吸引人的地方在于可以避免方位缺失带来的问题．

作为一个例子，我们考虑沿川西台阵31°N测线上19个台站的接收函数．台站分布如图 8a所示．其剖面总长度约为420 km．台站间距约为5～40 km(刘启元等，2009)．图 10给出了川西台阵31°N测线接收函数去噪结果及其与原始数据的对比．所用的地震为2006年12月7日的千岛群岛M_S6.4地震(震中位置46.26°N，153.80°E，深度16 km，发震时刻19h 10m 23.2s UTC)．提取接收函数时，仍采用了三分量多道最大或然性反褶积方法(刘启元和Kind，2004)．由于该事件在台站S04和S14的记录质量问题，我们用另外两个地震事件的相应数据加以替换．对台站S04，我们选用了2007年4月9日的千岛群岛M_S5.6地震(震中位置48.38°N，154.84°E，深度44 km，发震时刻10h 18m 2.3s UTC)．对台站S14，我们选用了2007年10月25日的千岛群岛以东M_S5.9地震(震中位置45.96°N，153.63°E，深度10 km，发震时刻为13h 50min 4.0s UTC)．它们的震中距和方位角与被替换数据大体相同．

图 10

Fig. 10

图 10 川西台阵31°N测线19个台站接收函数曲波阈值去噪 (a) 原始的径向分量； (b) 去噪后的径向分量； (c) 原始的切向分量； (d) 去噪后的切向分量． 图的左侧标出了台站代码． Fig. 10 Denoised receiver functions at 19 stations along the 31°N profile in the western Sichuan seismic array (a) Radial component before denoise； (b) Radial component after denoise； (c) Transverse component before denoise； (d) Transverse component after denoise．The station codes are labeled on the left of figures．

对于图 10可以有如下观察：1)接收函数的信噪比得到明显提高，同时震相沿测线的可追踪性增强；2)与松潘—甘孜地块内的台站不同(S05—S13)，四川盆地内各台站(S14—S19)的接收函数去噪后显示了简单的波形特征，这意味着四川盆地具有较为均匀且各向异性较弱的地壳，而松潘—甘孜地块内各台站接收函数较强的能量则与该地壳内的低速层及下地壳的变形方式有关(Liu et al．，2014)；3)观察零时刻前的背景噪声，并与图 9相比，我们发现对于阵列观测方式来说，其去噪效果不如单台观测的方式(例如，台站S05)的好．

通过把曲波变换和压缩感知引入接收函数的数据处理，本文发展了一种二维多道接收函数的噪声压制方法．数值检验和实际观测数据的测试结果表明，本文发展的接收函数曲波变换去噪方法有很好的可行性，它既可以用于单个台站观测的情况，也可以用于多台阵列观测的情况．上述两种情况，其压制横向非均匀散射的效果是明显的．这有利于基于水平分层地壳模型的接收函数研究，这意味着以往所需要的叠加数据量可以减少，从而减小由于叠加造成的多次波能量削弱，并有助于缩短流动观测所需的时间.

对以地壳各向异性研究为目标的接收函数偏振分析来说(齐少华等，2009)，这意味着大幅度减少了散射噪声的干扰，有助于提高接收函数偏振分析的稳定性．由于本文讨论的去噪方法可以同时弥补受地震方位控制的数据缺失，这有助于依据接收函数随方位变化规律的地壳各向异性研究(田宝峰等，2008；Schulte-Pelkum and Mahan，2014)．

但是，本文结果表明，单台观测数据的去噪效果优于多台阵列观测的情况．其原因在于，单个台站不同方位接收函数波场中，横向散射波在空间取向上的随机性强于线性剖面测线的情况．单事件的阵列观测对接收函数波场的空间采样依赖于台站布设方式. 相比之下，单台多方位接收函数波场的空间采样更加接近随机分布．理论上，曲波去噪趋向于对相干性较强信息的保护，而对在空间取向上随机性强的横向散射波趋向于压制.

对于我国大部分地区来说，全球地震分布导致的接收函数方位缺失不利于利用本文方法正确预测缺失方位的数据，但这并不妨碍单纯的去噪目的．如何进一步克服这种困难有赖于今后进一步的深入研究．