接收函数是地震记录水平分量与垂直分量反褶积得到的由间断面转换波和多次反射波组成的时间域波形序列.接收函数采用“等效震源假定”，与震源和地震波传播路径无关，是台站下方介质对近垂直入射地震体波的脉冲响应(Langston, 1979).远震接收函数被广泛用来求取间断面和速度结构特征，目前成为研究区域壳幔结构常规的数据分析工具(Chen et al., 2006; Wang et al., 2010; He et al., 2014; Si et al., 2017; Tian et al., 2014).远震P波接收函数主要利用30°~90°震中距和5.5~7.0震级的地震记录，反演壳/幔间断面深度(如Moho, 410 km, 660 km)和壳-上地幔顶部S波速度，其垂向分辨率可达到~公里级(徐强和赵俊猛，2008；Duan et al., 2017; Hu et al., 2015)；远震S波接收函数主要利用55°~85°(S)/85°~120°(SKS)震中距和5.5~7.0级的地震记录来研究Moho、LAB或410~660 km间断面深度或岩石圈地幔速度，其垂向分辨率可达到~10 km级(王峻，2016；Farra and Vinnik, 2000).由于地球内部结构的非弹性衰减作用，远震P波通常主要包含1s以上的信号，高频信号不足.与远震相比，近震(＜5°)数据通常能到达较高频率.然而，由于近震P波之后往往伴随多种震相，一般难以满足接收函数的震相单一假设条件.因此，借用接收函数概念的近震浅层结构研究不多(例如Li et al., 2014; Ni et al., 2014).然而，地表浅层(~m-km级)信息能为获得准确的深部结构以及防震减灾等提供先验信息，因此采用地震学方法(比如近震接收函数)约束浅层结构十分有必要(张宝龙等，2016；晏鹏等，2017; Su et al., 2015).而我国西部强震发生后，往往伴随滑坡次生灾害，可靠的浅层速度结构信息可为滑坡减灾研究提供重要的参考.

例如，2008年的汶川大地震在龙门山断层附近引发了数量众多的滑坡，对离震中较远的西山村马崩组渔湾组滑坡也造成了显著的影响.该坡体位于四川省阿坝州境东南部理县通化镇杂谷脑河左岸(图 1)，纵长约4200 m，横宽约1700 m，前后缘高程分别约为1500 m和2910 m.坡体为特大型土质滑坡，由第四系全新统崩坡积层和第四系残坡积层组成，厚度小于100 m，下伏基岩主要为泥盆系强风化~中风化千枚岩.滑坡体主要由处于强烈变形~蠕动变形状态且呈圈椅状的3个大滑坡组成，其中h1滑坡主滑方向为178°，总体坡度20.5°~32°；h2滑坡主滑方向为208°，坡度平均约22°；h3滑坡主滑方向为208°，坡度平均约25°(湖南省地质建设工程(集团)总公司，2010).在汶川地震中，该滑坡体上(与震中相距~50 km)部分居民房屋产生裂缝甚至倒塌.震后，滑坡变形范围进一步扩大，变形尺度进一步加剧，直接威胁滑坡上158户586人以及坡下的岷江一级支流杂古脑河及317国道安全，潜在经济损失约3000~5000万.马崩组渔湾组滑坡治理是一项非常紧迫的工程，也是国家重点基础研究发展计划(973计划)“西部山区大型滑坡致灾因子识别、前兆信息获取与预警方法研究”一个关键的研究区域.采用地震学方法，比如接收函数，研究马崩组渔湾组滑坡体的结构成为该973计划的一个重要部分.

图 1

Fig. 1

图 1 研究区域、地震台站和钻孔分布 (a)中黑色圆点标注钻孔位置(参考b)，蓝色三角形为地震台站，绿线为滑坡体边界；(b)为滑坡体h1中一条近南北向的工程地质剖面图(来自成都理工大学地质灾害防治与地质环境保护国家重点实验室). Fig. 1 Study region, seismic-station and borehole distribution Black dots in (a) show the location of the boreholes (b), blue triangles show the seismic-station location, green lines show the landslide boundaries. b shows a nearly N-S trending engineering geological cross-section in landslide h1 (from the State Key Laboratory of Geohazard Prevention and Geoenvironment Protection in the Chengdu University of Technology).

马崩组渔湾组滑坡本身的坡度特征和其与下覆岩石物性显著差异，使近震接收函数方法的近垂直入射假定和滑坡体底界面转换震相的形成得以满足，进而使利用该方法研究其结构特征成为可能.本文以滑坡体h1上4个钻孔资料得到的厚度作为约束，采用新发展的接收函数V_P-k叠加方法和传统的H-k叠加方法(Zhu and Kanamori, 2000)得到了台站下方滑坡体的平均P波速度和V_P/V_S比；基于上述结果，采用接收函数波形反演方法(Li et al., 2010)获得了台站下方滑坡体的S波速度结构，并探讨了滑坡体h1的结构特征及其相关的意义.

为了测试近震高频接收函数获取理县西山村滑坡体内部结构特征的可行性，项目团队于2014年7月在滑坡体h1上布设了4个相邻台间距~400 m的宽频带地震仪，并进行为期3个月的连续观测(图 1).这些地震台站紧邻该区已有的钻孔位置，便于对比地震学的研究结果和钻孔得到的滑坡体厚度和岩性结构.台站布设期间记录到大量高质量的地震事件，震源深度浅于15 km，震中主要在汶川县和理县，为提取近震接收函数提供了数据保证.基于研究区滑坡体和地震特征，参照中国地震局台网中心提供的地震目录，我们选取震中距小于0.5°，反方位角在0~100°和260~360°之间的1~3.5级地震，使得入射P波以近垂直入射.对这些地震，我们截取了P波到时前0.2 s和后0.6 s的三分量波形，以避免后续震相的干扰，然后人工挑选去除P波初至不清晰以及信噪比低的数据，便于下一步的接收函数计算.

传统远震P波接收函数的提取假设垂直分量与地表垂直，一般主要包括两个步骤(图 2a)：第一步把地震记录从原始的ENZ(东西、南北、垂向)三分量旋转到RTZ(径向、切向、垂向)三分量；第二步用垂直分量对径向和切向分量进行反褶积来提取接收函数.对于本研究区域的近震接收函数来说，由于滑坡体具有一定的坡度，地震数据波形的垂直分量不是垂直于滑坡体表面，其水平分量也并非平行于滑坡体表面，因此需要对原始地震三分量进行旋转.除此外，对于近震接收函数，由于震中距只有几十公里，地震定位误差也可能会影响接收函数计算.因此，我们在提取接收函数时，首先参考已有的研究，选取滑坡体倾角为26°和走向为-2°，把原始地震三分量(ENZ)旋转到平行于滑坡体表面(LTQ)，然后采用时间域迭代反褶积计算接收函数(Zhu and Kanamori, 2000，图 2b).同时，我们在PREM模型表层增加一层50 m厚和700 m·s^-1的P波速度和250 m·s^-1的S波速度的薄层，理论测试了发生在5~15 km深和20~35 km震中距的地震P波和Ps波的入射角，发现都小于10°，表明本文所用地震的定位误差对接收函数影响很小.图 2c展示了台站S01记录到的一个1.4级，震中距和反方位角分别为25.4 km和307.1°的地震波形原始三分量到提取接收函数的示意图.从提取的径向接收函数(RI)可以看出，在直达P波后的0.07 s左右有一个明显的Ps转换震相，暗示本文提取接收函数的可靠性.

图 2

Fig. 2

图 2 传统远震接收函数(a)和滑坡体近震接收函数三分量旋转示意图(b)，台站S01记录到的一个地震事件从原始三分量(E, N, Z)旋转到平行于滑坡体的径向分量(L)和垂直滑坡体分量(Q)以及径向接收函数波形图(RI)示意图 Fig. 2 Three components rotation schematic diagrams for traditional teleseismic receiver function (a) and the landslide local earthquake receiver function (b). (c) shows the diagram from original three-component (E, N, Z) for a seismic event recorded by station S01 to radial component (L) and vertical component (Q) and to radial receiver function (RI).

对于提取的径向接收函数，我们进行人工挑选，去除P和Ps震相不清晰的接收函数.S03由于仪器故障，没能获得有效接收函数.最终，我们保留了由3个台站记录到的共计46个地震提取到的63个径向接收函数(图 3a).这些事件的震中距主要分布在0.2°~0.3°，反方位角主要分布在75°~100°和295°~345° (图 3b).值得注意的是，与台站S04不同，台站S01和S02不同反方位角的Ps转换波到时和振幅都存在差异.前人研究表明(Aki and Richards, 1980)，波阻抗，倾斜特性，慢度，各向异性，局部流体等都会影响接收函数Ps震相振幅.研究区台站处于人类活动剧烈的滑坡体h1上，上述结构特征在观测台站的75°~100°和295°~345°反方位角范围内都有可能存在差异，进而导致台站不同方位角的Ps震相振幅存在差异，表明研究区结构存在小尺度的强烈横向变化.接收函数纵向分辨率影响反演的滑坡体结构的可靠性.我们截取包含Ps和PpPs震相的窗口(0.05~0.3 s)做频谱分析(图 3b).结果表明，S波除了在5 Hz处具有较高的能量外，一直到12 Hz甚至15 Hz都还有比较明显的能量.假定滑坡平均S波速度约300 m·s^-1(参考下文研究)以及主频为10 Hz，可得到滑坡体主能量波长约30 m，其垂直分辨率极限小于10 m(参考1/4波长分辨率标准).鉴于研究区滑坡体厚度都大于30 m，甚至超过60 m，我们认为本研究的数据能用来分辨滑坡体的界面和速度结构信息.

图 3

Fig. 3

图 3 本文用到的地震事件(a)和高斯系数为40的接收函数(b)信息. (b)中矩形框内展示了S波震相(0.05~0.3 s)的振幅谱 Fig. 3 Information for events (a) and receiver functions with a Gaussian value of 40 (b) for this study. The amplitude spectrums of the S-wave phases (0.05~0.3 s) are shown in the rectangular frames in (b)

在接收函数波形记录中，对于某一特定间断面，紧随初至P波之后主要有间断面Ps转换波和PpPs及PpSs+PsPs地表多次波三个震相.对于一维水平单层模型，上述三个震相与初至P波的到时差可表示为

(1)

(2)

(3)

其中，h为间断面到地表的深度，V_P和V_S分别为平均P波和S波速度，p为射线参数.Zhu和Kanamori等(2000)提出了利用Moho面Ps转换波和PpPs及PpSs+PsPs地表多次波共同约束地壳结构的接收函数H(地壳厚度)-k(V_P/V_S)叠加方法，被广泛用来获取Moho面深度与V_P/V_S信息(李永华等，2006；龚辰等，2016；Wang et al., 2014; Wei et al., 2011, 2015).该方法在给定的地壳平均P波速度前提下，搜索H-k平面内对应最佳地壳厚度与V_P/V_S的最大能量值(危自根等, 2011, 危自根和陈凌，2012).基于相似原理，当给定台站下方的间断面深度h时，同样可以获得地表到间断面的平均P波和S波速度.图 4表明，在给定滑坡体厚度下，方程(1)，(2)，(3) 在V_P-V_P/V_S平面上为三条斜率不同，相交于同一点的曲线，而这个相交点对应于真实的平均V_P与V_P/V_S比.在实际数据处理过程中，我们构造一个单层水平模型下V_P-k(V_P/V_S)平面内的叠加函数S(V_P, k)，

(4)

图 4

Fig. 4

图 4 间断面不同震相VP-VP/VS关系示意图 其中，间断面深度h=0.06 km，射线参数为p=0.15，Ps, PpPs和PpSs+PsPs震相与直达P波的到时差由平均P波速度0.08 km·s-1和VP/VS比3.0计算得出. Fig. 4 Schematic diagram of VP-VP/VS relation for different phases from discontinuity Discontinuity depth h is 0.06 km. Ray parameter p is 0.15.The time delays between the direct P wave and the Ps, PpPs and PpSs+PsPs phases are calculated by a average P velocity of 0.08 km·s-1 and VP/VS ratio of 3.0.

其中，r(t)为径向接收函数，t_Ps，t_PpPs，t_PpSs+PsPs为不同震相的理论到时，w_Ps，w_PpPs，w_PpSs+PsPs为满足三者之和等于1的加权系数.在V_P-k域内，当给定层厚度h时，V_P-k平面上最大叠加能量对应最佳的V_P和V_P/V_S.

钻孔资料提供的滑坡体厚度信息为我们采用上述接收函数V_P-k叠加方法获取台站下方滑坡体的平均P波和S波速度提供了保证.基于滑坡体的坡度范围和接收函数波形特征，我们采用钻孔厚度的90%作为初始输入h和0.7, 0.3, 0.0分别作为Ps, PpPs和PpSs+PsPs震相的加权因子，采取上述发展的V_P-k方法反演了3个台站东西两侧的平均P波速度和V_P/V_S比(图 5).图 5a为三个台站不同反方位角叠加后的接收函数波形，黑色和红色圆圈分别标明与V_P-k叠加结果(图 5b，白色十字叉)相对应的Ps和PpPs震相.图 5b为V_P-k叠加方法搜索能量分布图，白色十字叉标明最终的V_P和V_P/V_S比.值得注意的是，在反方位角295°~345°情况下，台站S01和S04最大的叠加能量对应于图中白色五角星处，其对应的Ps和PpPs震相到时分别为0.16 s和0.26 s，以及0.26 s和0.47 s.鉴于五角星对应的Ps震相都不是紧邻直达P波后的震相，我们认为其对应的V_P-k值为一假象.上述的V_P-k叠加结果表明，三个台站下方滑坡体平均P波速度在500~800 m·s^-1之间，V_P/V_S比分布在2.4~3.1之间.高的V_P/V_S比与浅表结构高富集的含水量和孔隙比从而导致其比地壳平均1.78的波速比偏高一致(高武平等，2014).进一步分析发现，平均P波和V_P/V_S比在研究区东西两侧(不同反方位角)都存在较大差异.对台站S04，尽管东西两侧平均P波和V_P/V_S比差异较小(＜5%)，但是东侧Ps/P比值几乎是西侧的7倍，表明滑坡体东西两侧底界面性质可能存在显著差异.对台站S02，尽管东西两侧平均P波速度较为接近，但西侧V_P/V_S比比东侧高~21%；对台站S01，东西平均P波速度和V_P/V_S比差异分别达到39%和5%.值得注意的是，除了速度之外，台站S01和S02东西Ps和PpPs震相大的到时差还有可能是由滑坡体底界面非水平特征引起.图 5c展示了当平均P波速度为800 m·s^-1时，采用H-k叠加方法得到的台站东西两侧的厚度与V_P/V_S比分布.可以看出，对台站S01，当平均P波速度相同时，反方位角75°~100°下方的滑坡体厚度比295°~345°下方厚21 m时，同样可以得到如图 5a中S01台站所示的接收函数.基于上述因素和本文获得的台站S01速度信息，我们采用TauP软件计算了台站S01在34 m深度的P和Ps穿透点位置，发现其横向采样范围要小于4 m.也就是说，滑坡体底界面倾角要超过65°才能导致本文观测到的台站S01的接收函数波形.综上所述，虽然不能排除台站S01下方滑坡体底界面的倾斜性质，但东西两侧速度的不同是造成不同方位角接收函数波形差异的一个重要原因.

图 5

Fig. 5

图 5 不同反方位角叠加后的接收函数(a)和VP-k叠加结果(b)以及H-k叠加结果(c).(a)中圆圈标明和(b)图VP-k叠加结果一致(白色十字叉)的Ps和PpPs震相到时. (c)图为假设P波为800 m·s-1，采用H-k叠加得到的滑坡体厚度与VP/VS比 Fig. 5 Receiver functions stacked by different back azimuth (a) and VP-k stacking (b) and H-k stacking results (c). Circles in (a) show the arrival times for Ps and PpPs phases corresponding to the VP-k stacking results (b, white crosses). (c) shows the thickness and VP/VS ratio for landslide using the H-k stacking method with (a) P wave velocity of 800 m·s-1

上述采用接收函数V_P-k方法获得了滑坡体的一些深部结构信息，然而主要局限于滑坡体平均的P波速度和V_P/V_S比，要获得详细的不同深度S波速度结构还需要采用接收函数波形反演方法(Chu et al., 2010).接收函数差异演化非线性全局优化反演方法采用两组随机产生的参数向量之差作为产生下一代个体，直接搜索不同参数组合使目标函数收敛于极小值(Li et al., 2010)，能较好的约束5层以下的S波速度结构，被Wei等(2016)和危自根等(2016)用来反演四川盆地和土耳其地壳结构，也被本文用来反演滑坡体下的S波速度结构.滑坡体被分为4层进行反演，厚度由钻孔资料约束，每层的P波速度和V_P/V_S比参考接收函数V_P-k叠加结果，并分别在±0.2 km·s^-1和±1.0范围内变化.图 6展示了在0~0.3 s的时窗内实际(实线)和拟合(虚线)的接收函数以及反演得到的S波速度结构.从接收函数波形拟合程度来看，除了台站S04(295°~345°)实际和拟合波形的Ps震相到时存在较大差异外，其余台站波形都拟合较好，尤其是Ps和PsPs震相.对于台站S04(295°~345°)，其有效的接收函数数量只有2个(图 3)，V_P-k叠加方法存在双极值现象.上述特征，结合较大差异的实际和拟合的接收函数特点(图 6b)，我们认为其反演得到的S波速度结构不一定反映S04(295°~345°)台站下方滑坡体的真实结构信息.

图 6

Fig. 6

图 6 接收函数波形反演示意图 黑实线和虚线分别展示了原始和拟合的接收函数波形，圆圈标注了原始数据中Ps和PpPs震相，矩形框内展示了反演的S波速度. Fig. 6 Schematic diagrams for waveform inversion of receiver functions The black solid and dotted curves show the original and recovered receiver functions, respectively, circles show the Ps and PpPs phases in original receiver functions, inverted S-wave velocity models are shown in rectangular boxes.

已有研究表明，S波速度结构与电阻率存在一定的关系，相对低的速度往往对应着低的电阻率特征.Xu等(2016)采用时移电阻率方法反演了h1滑坡体一条近N-S向剖面下方的电阻率结构，获得了台站S01和S04以及钻孔ZK07(与台站S02距离不超过30 m，图 1)的电阻率信息，为本文进行S波速度和电阻率对比提供了数据基础(图 7).对台站S04(75~100°)，在45~63 m深度S波速度为125 m·s^-1, 是上层S波速度35%，与滑坡体底层低的电阻率特征一致.对台站S02，在滑坡体最下层速度只有78 m·s^-1 (75°~100°)和143 m·s^-1(295°~345°), 分别都不到上层速度的40%，与滑坡体底层的低电阻率对应.对台站S01，其75°~100°反方位角内滑坡体底层不是低速特征，而在295°~345°反方位角内，滑坡体底部存在95 m·s^-1的低速层，只有上层的38%，这与其下方滑坡体底层高低电阻率剧烈的横向变化一致.本文观测到的上述三个台站滑坡体底部的低速特征与前人观测到的底部低电阻率特征良好的一致性间接支持本文与前人结果的可靠性.工程地质研究表明(四川省阿坝州理县通化乡马崩组渔湾组滑坡应急治理工程可行性研究报告)，滑坡体h1目前处于强变形~蠕动变形阶段，人类活动强烈，生活和饮用水乱排乱放，地面水流常年不息，不停入渗滑体，在基覆面附近富水.上述特征与本文观测到的2.4~3.1的高V_P/V_S比，以及78~143 m·s^-1的滑坡体底部低的S波速度结构特征是一致的.本文观测到的台站S01，S02和S04下方滑坡体的S波速度和V_P/V_S结构，结合其富水特征，表明滑坡体h1底界面的抗剪强度相对较弱，稳定性较差，是潜在的滑坡危险区域.

图 7

Fig. 7

图 7 台站下方电阻率(Xu et al., 2016)和S波速度结构示意图 Fig. 7 Schematic diagrams for resistivity(Xu et al., 2016) and S-wave velocity structure for the stations

本文把在理县马崩组渔湾组滑坡h1上4个宽频带地震仪记录到的三个月的近震数据旋转到平行于滑坡体的径向、切向和垂直滑坡体垂向的三分量波形，采用时间域迭代反褶积方法计算接收函数.基于挑选后的接收函数，采用钻孔资料得到的滑坡体厚度做为约束，我们发展了接收函数V_P-k叠加方法，并结合接收函数H-k叠加方法得到了滑坡体3个台站下方平均P波速度和V_P/V_S比，并进一步采用接收函数波形反演方法获得了台站下方滑坡体S波速度结构特征.结果表明，与台站S02不同，台站S01和S04的接收函数Ps和PpPs震相无论振幅还是到时存在着明显的东西差异(75°~100°和295°~345°).三个台站平均P波速度在500~800 m·s^-1之间，V_P/V_S比分布在2.4~3.1之间，且在滑坡体底层都存在78~143 m·s^-1左右的S波低速层.本文观测到的台站下方滑坡体的低速特征与前人观测到的低的电阻率结果一致.上述高的V_P/V_S比和底层低的S波速度结构，与工程地质调查得到的强变形~蠕动变形特征和富水的滑坡体底部结果一致，表明滑坡体h1底界面的抗剪强度相对较弱，稳定性较差，发生滑坡坍塌的可能性相对要大.

由于滑坡体地震流动台站施工困难等原因，本文只获得了h1滑坡体布设的4个台站中3个台站下方的速度结构信息，不足以全面揭示滑坡体结构特征.在滑坡体上布设更加密集和更高规格的台站进行长时间的记录，反演高分辨率的滑坡体速度结构信息，甚至探索结构随时间的变化特征，可为滑坡体的治理提供基础信息.然而，本文的方法要求滑坡体附近有较多数量的近震，对于地震活动性低的地区(例如福建、广东)的适用程度有限.而在我国西部地区的强震区，余震或较高的背景地震活动性可提供丰富的地震事件波形，为近震接收函数提供基础.同时，本文研究的案例为土质滑坡，覆盖层与滑床介质差异明显，且覆盖层厚度变化较为平缓.对于岩质滑坡或者覆盖成三维结构复杂的情况，还需开展深入的研究.

中国科学院测量与地球物理研究所提供流动地震数据资料，成都理工大学地质灾害防治与地质环境保护国家重点实验室提供钻孔资料，中国科学院测量与地球物理研究所倪四道研究员提供了诸多有益的建议，审稿人提出建设性的建议，在此一并致谢！