0 引言

天然地震的远震波形由震源破裂过程时震源区混响效应、传播过程中介质响应、仪器响应以及接收区下方介质响应共同决定.Langston(1979)基于等效震源假定从中长周期远震体波中分离出接收函数.Owens等(1984)将其扩展到宽频带记录，得到宽频带接收函数.所谓接收函数就是台站下方接收区速度间断面产生转换波到时的时间序列.Ammon等(1990)发展了保有接收函数径向分量绝对振幅的接收函数分离方法.刘启元等(1996)提出了接收函数复谱比的最大或然性估计方法，随之发展至基于多台的多道最大或然性反褶积方法，在复杂构造环境条件下成功的改善接收函数的估计(刘启元和Kind，2004).在时间域内，通过使用Wiener滤波、最大熵谱反褶积等方法也可获取高质量的接收函数(Ligorría and Ammon, 1999；吴庆举等, 2003a, 2003b).由于接收函数的易提取、稳定和高分辨率特点，该方法在研究地壳上地幔速度结构中已经得到广泛应用(Zhu and Kanamori, 2000；Julià et al., 2000；Vinnik et al., 2004；Liu et al., 2014).

偏移作为地震数据处理的主要技术，在地震勘探应用中已成功的获取复杂的地下结构.该方法通过确定散射点位置并且恢复波形与振幅特征来实现(郑海山，2000；刘喜武和刘洪，2002；李振春，2014；撒利明等，2015).随着勘探中偏移技术的理论发展愈发成熟，众多学者开始将其应用于天然地震的观测数据中，利用偏移技术研究地壳上地幔的间断面结构.接收函数深度域偏移成像和时间域偏移叠加技术被率先发展(Dueker and Sheehan, 1997；Yuan et al., 1997).因地壳内部横向不均匀性较强，没有地下水平分层假设的波动方程偏移、逆时偏移和基于散射核(scattering kernel)成像等方法也相继被提出(Bostock and Rondenay, 1999；Ryberg and Weber, 2000；Wilson and Aster, 2005；Chen et al., 2005；吴庆举等，2007；Shang et al., 2012；Mancinelli and Fischer, 2017).

在上地幔过渡带间断面即410 km、660 km速度界面的研究中，由于该深度范围内间断面横向变化较弱，基于水平分层介质假设的共转换点叠加(CCP)法更容易获取其结构，故在上地幔过渡带中CCP作为偏移成像的普遍方法(Zhu, 2000；Zhang et al., 2010；Yue et al., 2012).在该方法中，参与叠加的接收函数频率决定成像的垂向分辨率，但由于P₄₁₀s和P₆₆₀s的振幅在噪声干扰下相对较弱，信噪比较低，通常需要低通滤波至5 s甚至更低来对其进行成像.使用接收函数的频率降低，意味着成像分辨率随之下降，从而对上地幔过渡带的厚度估计产生直接影响.一般认为，上地幔410 km和660 km间断面是由上地幔岩石相变产生的速度间断面(Anderson, 1989).其厚度直接决定上地幔过渡带的温度异常与否，所以消除接收函数中随机噪声对上地幔转换震相识别的干扰，给出分辨率更高的偏移成像结果，对研究地球内部构造及深部动力学有着重要的意义.

曲波变换是发展迅速的多尺度、多方向信号分析工具(Candès and Donoho, 2000，2004).相比传统二维数据去噪方法，对信号的损伤最小.基于更好的方向性特征识别，曲波变换尤其适用于地震同相轴检测与去噪(Candès et al., 2006；Herrmann and Verschuur, 2004；Herrmann et al., 2007).由地壳介质横向不均匀引起的散射噪声经曲波变换去噪后得到压制，接收函数信噪比提高，震相的可追踪性改善明显(齐少华等，2016).因此，利用曲波变换对接收函数去噪是压制接收函数随机噪声行之有效的方法.本文通过对川西地区部分流动台站的接收函数使用曲波去噪后进行偏移成像，以检验该方法在提高偏移成像分辨率的可行性.

1 方法 1.1 接收函数的动校正与叠加

对大量可用的接收函数叠加可加强信号，但动校正处理是接收函数上地幔转换震相叠加的基础，也是接收函数叠后偏移的基础.在地壳深度范围，且震中距变化范围小于20°时，可直接通过接收函数叠加得到平均接收函数(Dueker and Sheehan, 1997；Yuan et al., 1997).但在上地幔深度，由于转换波Ps射线路径较长，在叠加前必须要通过动校正(moveout correction)来消除慢度差异.我们在地壳上地幔深度范围内使用三维速度模型(Liu et al., 2014)，100 km以下的深度范围内使用IASP91速度模型进行动校正计算.

假设地球介质均匀分层，V_P、V_S分别为P波和S波波速，p为P波的水平慢度.在深度z处间断面上形成的PS转换波在接收函数上相对直达P的走时为：

(1)

如果以p₀参考水平慢度，那么在深度z产生的Ps转换震相相对P波走时应做以下校正(Yuan et al., 1997)：

(2)

本文研究使用的参考慢度是p₀=6.4 s/(°).

1.2 曲波去噪

曲波变换(curvelet transform)是近来发展迅速的一种数学变换，曲波被认为是具有方向性的小波，可以对曲线形边缘提供最佳的稀疏表示(Candès and Donoho, 2004).所以，曲波作为理想的手段来检测地震数据的波前和压制噪声(Herrmann and Verschuur, 2004).它在空间域和频率域内均满足抛物尺度关系，即宽约为长的平方.在空间域内，形态接近平面波.在频率域内，曲波表现为一个围绕原点的离散楔形区域，组成紧支撑框架.图 1中我们利用curvelet工具包在空间域和频率域内给出四个不同尺度、方向和位置的曲波.曲波函数可定义为(Candès et al., 2006)：

(3)

式中曲波φ_{j, l, k}对应三个参量(j, l, k)：尺度、角度和位置.尺度变化满足2^－j，角度变化满足θ_{j, l}=2πl·2^[j/2]，且0≤θ_{j, l}≤2π.其中0≤l≤2^[j/2]－1，[j/2]代表j/2的整数部分；位置参数满足x_k^{j, l}=R_{θj, l}^－1(k₁·2^－j, k₂·2^－j/2)；R_θ为旋转矩阵，k₁和k₂为平移参数. Candès等(2006)发展了两类实现离散曲波变换的途径：USFFT算法(Unequispaced FFT)和本文使用的Wrap算法(Wrapping-based transform).

基于Wrapping的快速离散curvelet变换的实现主要有四个步骤：(1)通过二维傅里叶变换获取接收函数在频域系数, －n/2≤n₁, n₂≤n/2；(2)每个尺度、方向参数(j, k)有，在频率域获取角度楔形；(3)围绕原点对角度楔形进行装配；(4)对进行二维FFT逆变换，获得离散的曲波系数.

基于曲波变换去噪的步骤如图 2所示，在单个台站数据处理中，假设C为接收函数矩阵，且该矩阵的行、列分别为时间序列和数量信息.R作为与C相同大小的矩阵，行、列也同为时间与数量信息.上文提到过，叠加是压制噪声、增强信噪比的步骤.所以，R矩阵中每道都为接收函数动校正后叠加所得的参照接收函数, N为相同规模的噪声矩阵，假设：

对矩阵C、R进行曲波变换，得其曲波变换系数，分别为c_{j, l, k}和r_{j, l, k}.对于曲波的每个尺度j和角度l，对系数r_{j, l, k}进行归一化处理后得到n_{j, l, k}.最后对c_{j, l, k}·(n_{j, l, k})^m进行曲波逆变换后得到去噪后的接收函数矩阵，m值的大小直接决定接收函数去噪效果，所以我们通过测试决定m值的大小.在环境噪声研究中，也有学者通过类似方法对互相关函数进行去噪处理(Stehly et al., 2015).

1.3 共转换点叠加(CCP)

接收函数偏移的本质是将Ps转换波震相与直达P波的延时映射到发生转换的深度.而共转换点叠加(CCP)是假设地下为水平层状介质或弱横向非均匀介质，通过参考速度模型进行射线路径追踪，然后把接收函数上每个Ps转换波振幅相对应的被校正后的走时归位于发生转换的位置，也就是速度界面所在的地方(Zhu, 2000).对于给定地壳速度模型中, 若假定接收函数径向分量只包含Ps转换波产生的震相, 在均匀分层的地球介质中，x_S为Ps转换波横向传播距离，t_Ps(z)是在深度为z时接收函数中Ps转换波相对直达P波走时.时间域地震事件j的接收函数r_j(t)可以通过

(4)

映射到深度域.这里, I(z, x)为当地下深度z和水平偏移距离为x_S(z)时的成像强度, 求和被限定在深度z处PS转换点处于区域S内的所有接收函数r_j在时间点t=t_Ps(z)的值.本文将通过动校正后使用曲波变换去噪得到的接收函数，进行共转换点叠加，在0~100 km深度范围内同样使用接收函数与背景噪声联合反演的三维速度模型(Liu et al., 2014)，获取31°N剖面的地壳上地幔的速度间断面结构.

2 数据及去噪测试 2.1 数据

本文使用的数据是沿川西台阵31°N测线附近31个台站的2195条P波接收函数，该接收函数的获取基于刘启元和Kind(2004)提出的三分量接收函数的多道最大或然性反褶积方法.台站分布如图 3所示.剖面长度约为670 km，台站间距20~40 km，观测时间为2006年12月—2009年4月.

2.2 去噪测试

为检验接收函数曲波变换去噪的效果，本文选取了KCD08(30.7°N，104.5°E)台站的数据进行测试.通过对接收函数滤波、动校正、曲波去噪，达到去噪目的.

上文提到曲波变换去噪中m值的选取至关重要.我们通过对实际数据进行去噪测试，试图确定合适的m值.图 4a中给出KCD08台站低通滤波至2 s的接收函数, 以及分别选取m为0.05、0.15、0.25、0.35和0.45时去噪后的数据.由图 4b、c、d、e和f可知数据得到不同程度的降噪，选取不同m值去噪后数据的410 km和660 km的转换振幅都相对增强.且并未改变转换波震相的走时信息，同时P₄₁₀s和P₆₆₀s有很好的可追踪性.我们通过比对后选取m=0.25，在实现去噪的同时保留较多原始数据的走时信息.

在前人的研究中，一般选取4 s或5 s的接收函数进行偏移成像.虽然周期为1 s、2 s的接收函数含有更多信息，因有噪声较多或无关震相干扰，使得上地幔速度间断面转换震相并不清晰.将接收函数滤波至0.33 Hz、0.25 Hz，可获得清晰且可追踪的410 km与660 km间断面的转换震相，可是使用低频接收函数伴随着偏移成像垂向分辨率的降低.我们通过使用曲波变换压制高频接收函数的噪声，试图提高参与CCP的接收函数频率.图 5a为KCD08台站滤波至1 s的接收函数图像，仅Moho面转换震相清晰可追踪，其转换波到时在5 s前.图 5b为曲波去噪后的接收函数，经去噪之后，主要速度界面转换震相并未被削弱，而无关震相被压制明显，可直接从接收函数中识别主要转换震相，且其有较好的可追踪性，在叠加后的接收函数中可清楚识别在44 s、69 s存在410 km与660 km速度间断面的转换震相.通过对1 s接收函数去噪，其信噪比显著提高，获取更为直观的震相信息，相信我们的去噪方法在处理相对高频的接收函数中也可有较好前景(Liu et al., 2017).在本文研究工作中将1 s接收函数用于最终结果的0~100 km偏移成像.图 6a为滤波至0.5 Hz的接收函数，相比图 5a的1 s周期的接收函数，可知通过滤波也可实现简单去噪.图 6b同样对2 s接收函数进行曲波去噪处理，上地幔过渡带速度界面的转换震相得到明显加强.本文使用去噪后2 s的接收函数用于100~800 km的偏移成像.

2.3 接收函数偏移成像

通过曲波变换去噪获得明显改善的接收函数，并进行共转换点叠加，以检验去噪对偏移成像的影响.图 7a、b分别为最小周期4 s的接收函数在去噪前后的成像结果，结果比对可知去噪对偏移成像也有明显的改善，410 km、660 km的速度界面主要特征保持一致，对界面识别造成干扰的噪声得到有效去除.

图 8a为周期2 s的接收函数叠加偏移的结果，通过使用更高频率的接收函数，使得图 7a成像分辨率提高，但伴随整体噪声的增强.图 8b是使用1 s、2 s接收函数去噪后分别得到0~100 km和100~800 km的速度间断面图像，相比图 8a界面形态仍保持一致，而且去噪效果显著，成像分辨率明显提高，对速度间断面的识别更加准确.

通过对图 7、图 8比对可知，曲波变换去噪会显著增强成像效果，保证有效信息高度一致的同时提高偏移使用的接收函数频率，安全地将偏移使用的接收函数主周期扩展到2 s，获得更高的垂向分辨率的偏移图像.

3 31°N剖面Moho、410 km和660 km速度界面

川西地区位于青藏高原东缘与华南地块交汇区域，被鲜水河断裂与龙门山断裂分割为川滇地块、松潘—甘孜地块和四川盆地，而且是中国内陆地震频发的主要地区之一.2008年5月12日汶川地震，便位于该区域的龙门山断裂带之上.汶川地震之后，众多国内外地球物理学家对川西区域进行深部探测研究，并获得丰富的成果.通过接收函数反演、地震环境噪声以及走时层析成像等方法探测表明，四川盆地中上地壳有明显高速异常，而在松潘—甘孜地块中下地壳有低速物质且有向东流动趋势，但受到坚硬的四川盆地阻挡，应力缓慢累积，最终导致龙门山断裂带破裂(郭飙等，2009；吴建平等，2009；雷建设等，2009；李昱等，2010；Liu et al., 2014). Zhang等(2009)和Robert等(2009)分别对横跨龙门山断裂的地壳间断面结构进行了接收函数偏移成像研究，结果给出了龙门山下方Moho面的横向变化和中下地壳分界的迹象.因川西地区在地壳深度范围内有很强的横向不均匀性，使用一维速度模型偏移必然会在上地幔间断面成像结果中有所偏差.所以本文通过使用接收函数与背景噪声联合反演的三维速度模型(Liu et al., 2014)和IASP91速度模型(＞100 km)，利用经曲波变换去噪后的接收函数进行叠加成像，试图获取川西地区更为准确的上地幔间断面结构.

图 9给出在川西地区沿31°N剖面0~800 km深度范围内的速度间断面结构图像.地表至100 km深度使用1 Hz接收函数偏移成像，在100~800 km深度范围采用主周期为2 s的接收函数.由图 9可见，Moho速度界面自东向西随地势增高逐渐变厚，在四川盆地深度约为35~45 km，而在川滇地块将近60~70 km，且在龙门山断裂带下方Moho面深度存在明显的突变.在青藏高原东缘下方70~80 km深度附近，连续的Moho界面下方可发现一条连续的正震相界面，标注其为M1.Liu等(2014)给出相同区域的S波速度结构中，在该剖面的M1位置附近S波速度明显增加，这与我们文中偏移结果正震相界面相互对应，这从另一个角度提供了青藏高原东缘Moho面具有10~20 km厚的过渡性Moho的证据.

410 km速度间断面形态清晰，沿整个剖面位于400~405 km深度，在青藏高原东缘下方深度稍浅，而在四川盆地下方略深.410 km间断面成像强度和间断面厚度比较均匀，只是在松潘甘孜下方经度约为103°位置有所减弱并变宽，并且其深度有稍稍变浅的迹象，但与剖面其他部位相比深度变化不超过5 km.660 km间断面同样沿整个剖面没有明显的深度变化，除在2°~3°的区间外均清晰地出现在660 km深度.在2°~3°范围及松潘—甘孜地块下方，660 km间断面则出现较明显的变深的迹象，其深度变化约为10 km.由410 km间断面和660 km间断面决定的上地幔过渡带的厚度在本文研究的东西向大约650 km范围内表现出高原下方稍厚、四川盆地下方稍薄的横向变化特征，上地幔过渡带厚度在高原下方比标准地球模型厚约10 km，而在四川盆地或扬子地块下方则稍厚约5 km.我们在进行接收函数偏移的过程中采用了更加符合川西地区实际速度结构的参考速度模型，在最终的上地幔间断面成像结果中并未发现在青藏高原东缘和扬子地块之间的明显的间断面深度差异(Zhang et al., 2010).

一般认为，上地幔410 km和660 km间断面是由上地幔岩石相变产生的速度间断面(Anderson, 1989).其中，410 km间断面由橄榄石(olivine)的α-到β-相变产生，而660 km间断面则是由γ-尖晶石到钙钛矿(perovskite)的相变产生.根据上地幔岩石相变的Clapeyron方程，相变的发生受介质环境温度和压力的影响.在410 km深度，相变发生的温度-压力曲线的斜率为正，而在约660 km深度处发生的γ-尖晶石到钙钛矿的相变温度-压力曲线斜率为负值.410 km深度和660 km深度不同的相变特性造成的综合结果是，当上地幔过渡带深度温度存在正异常时，上地幔过渡带厚度将减小；当上地幔过渡带深度温度低于正常温度时，上地幔过渡带厚度增加.基于这一认识，本文结果给出青藏高原和扬子地块碰撞交接区域的上地幔过渡带较标准地球模型厚约10 km，所以其上地幔具有低于全球平均值的温度.如果只考虑温度对过渡带厚度影响，每100 K的温度异常导致约10 km的厚度变化(Helffrich，2000；Kind et al., 2002), 同时0.43%的P波速度异常对应于-100 K的温度变化(Cammarano et al., 2003).东亚地区上地幔层析成像结果(Li et al., 2010)给出的川滇地区上地幔具有高速异常，这与本文得到的结果是一致的.

在龙门山断裂带下方410 km(103°)和660 km(101°~103°)间断面表现出了比较复杂的构造，在这一位置410 km和660 km速度间断面宽度较大，并分别呈现出变浅和变深的迹象.在该区域上地幔过渡带由东向西逐渐变厚.造成这一现象的一个可能的原因是由于龙门山断裂带浅层结构过于复杂，造成410 km和660 km间断面成像不好.但在410 km和660 km深度，台站之间接收函数数据的横向覆盖范围(约为150~300 km)已经远远超过了龙门山断裂带的宽度，因此，浅层复杂速度结构不大可能是上地幔过渡带的这种间断面深度变化的决定因素.造成这个小区域过渡带厚度变化的另一个可能的原因则是过渡带附近的温度变化，根据Li等(2010)给出川滇地区上地幔过渡带的高速异常的结果，我们推测该区域岩石圈深度范围的拆沉可能已进入上地幔过渡带，引起该区域局部温度降低，导致过渡带增厚，但需要更多不同学科证据进一步论证.

4 结论

本文通过对接收函数曲波去噪后进行叠加偏移，得到川西区域31°N剖面的速度间断面结构.通过研究结果有以下认识：

(1) 沿31°N纬度线的接收函数偏移成像结果证明了曲波变换滤波的实用性，有效压制了随机噪声和横向散射, 改善接收函数转换震相的可追踪性，数据信噪比显著提高.

(2) 4 s周期的接收函数经曲波变换去噪后，偏移叠加成像效果明显增强.同时上地幔深度的偏移成像可以扩展到使用0.5 Hz的接收函数.采用短周期接收函数使得成像结果的垂向分辨率具有明显提升.

(3) 偏移叠加结果与接收函数和环境噪声联合反演得到的S波速度结构取得的Moho面形态认识基本一致，从另一个角度证实了青藏高原东缘地壳底部到上地幔顶部具有10~20 km厚的过渡性界面，这可能意味着在川滇块体下方上地幔顶部和中下地壳具有物质或热交换.

(4) 间断面成像结果发现在青藏高原东缘和扬子地块交界地区上地幔过渡带厚度较全球平均过渡带厚度厚5~10 km，表现为青藏高原东缘下方厚、四川盆地较薄的特征，显示在该地区上地幔过渡带附近具有较低的温度或较高的地震波速度.而在龙门山断裂带下方约150 km宽度范围内具有更厚的过渡带，这一结果与上地幔层析成像结果一致.在龙门山断裂带以西的小范围内有可能存在下地壳或上地幔物质的拆沉.

值得说明的是，由于上地幔间断面一般并非简单的一级间断面，使用短周期接收函数进行上地幔偏移成像会造成深部界面成像变弱，我们认为综合考察各种不同频带的成像结果对正确理解结构特征是重要的.

致谢  作者感谢两位评阅人认真的评审和提出宝贵建议，使得本文得以进一步完善.