远震P波波形数据包含大量地震台站下方地壳和上地幔速度间断面所产生的Ps转换波及其多次反射波的信息，接收函数是利用其垂直分量对径向分量和切向分量作反褶积处理后得到的时间序列，代表台站下方地壳上地幔速度结构远震平面波的响应，基本与震源及传播路径无关，是研究地球内部物理、地壳上地幔结构的有效方法，并已得到了广泛应用(吴庆举和曾融生，1998；吴庆举等, 2004; Kosarev et al., 1999; 刘启元等, 2000；杨毅和周蕙兰, 2001; 段永红等, 2005; Tian et al. 2005a, 2005b; 胡家富等, 2005; Chen and Ai, 2009; 房立华和吴建平, 2009; 查小惠和雷建设, 2013).接收函数的提取原理简单，难点在于如何使用一种稳定、可靠的算法来获得.在接收函数方法应用发展过程中，多种提取的算法被提出，但每一种方法都有自身的局限.

频率域反褶积( Helmberger and Wiggins, 1971)是其中常用的一种方法，即用远震P波的水平分量与垂直分量的谱振幅相除并反变回时间域得到接收函数.由于地震仪弹性系统的固有周期是有限的，因此实际地震资料是有限带宽的，这就导致垂直分量中一般都存在近零的频率成分，直接在频率域作除法运算往往不稳定，为确保频率域反褶积的稳定性，故引入了‘水准量’对垂直分量作预白化处理，以压制近零的频谱成分(Langston et al., 1979)，但实际应用中，‘水准量’的选择往往存在着困难，过小难以保证其稳定性，过大就不可避免地压制了能量较弱的转换波震相，特别是上地幔间断面的微弱震相；另一方面，快速傅里叶变换假定了时间信号是周期无限的，但我们在计算接收函数时往往只截取其中的一部分，这就隐含着时窗之外数据为零或呈周期性变化的不合理假设，并且单一窗函数的选择不得不在频率泄露和频谱估计方差之间进行权衡(Park et al., 1987)，这影响了接收函数的计算精度和分辨率.为此，一系列提取接收函数的方法被提出：包括吴庆举等提出的时间域Wiener滤波(2003a)、最大熵谱反褶积(2003b)频率域多道反褶积(2007)，刘启元和Kind(2004)的最大或然性反褶积方法以及Park (1987)和Levin (2000) 提出的多椭球窗频率相关估计方法.时间域反褶积主要受信号中最大振幅频率成分的影响，而且在应用中通常将数据频率限制在0.5Hz以下，虽然能够获得较好的410km 和660km间断面信息，但对于地壳内结构存在着困难(Dueker and Sheehan, 1998).Park和Levin的多椭球窗频率相关估计用一簇数据窗代替单一数据窗，很好地抑制了频率泄露、保证了数据信息的饱和度，而且把同等长度的震前数据作为计算接收函数的一部分，避免了在使用‘水准量’过程中产生的困难.与频率域反褶积和时间域反褶积方法相比，获得的接收函数精确度更高、稳定性更强，但用此方法获得的接收函数，可用的信息仅仅包含在前10s内，10s 以后接收函数的信息被压制(Helffrich et al., 2006)，这阻碍了利用接收函数方法对深部结构(例：地幔过渡带)的研究.为克服这个缺点，一种时间延长的多椭球窗频率相关估计方法被提出(Helffrich et al., 2006; Takuo et al., 2008)，对上述方法进行了改进，在保证精确度和稳定性的前提下，可以计算任意时间长度的接收函数.

多椭球窗频率相关估计方法的关键在于椭球序列(Slepian, 1978)的确定.椭球序列依赖于两个参数：带宽W和时间带宽积P(详见2.2节)，在计算接收函数时要选择合适的带宽W，来获得高的频率分辨率以及小的方差，以获得高质量的接收函数，但其没有确定的表达式，改变参数W值，椭球序列必须重新计算，给实际应用带来了困难.而Reidel提出的多正弦窗序列 (详见2.3节) ，在抑制频谱泄露方面，与多椭球窗取得同样的效果，并且具有更小的局部偏差(Riedel and Sidorenko, 1995).另外，多正弦窗有确切的表达式，其带宽W的改变只需增加或减少多正弦窗的个数就可以实现，在实际应用中更为简便.

在本文中我们用多正弦窗代替多椭球窗序列，并且采用最优等效震源时间函数长度(张洪双等, 2009)来计算接收函数，在保证了接收函数的稳定性和精确度的前提下，提取更加方便.

接收函数的频率域表达式为

将Ē_R(ω)和Ē_T(ω)反变回时间域就得到了径向和切向接收函数.

离散椭球窗(Slepian, 1978)(又称Slepian窗)是1978年由Slepian提出的一簇相互正交的离散椭球序列ω^(k)，已经成为多窗谱分析中最常用的一簇数据窗，得到广泛的应用(陈赟等, 2005; Anandan et al., 2004; Zhu et al., 1989).

假定有一实数序列v(n)，n值为[0,N-1]，其能量为1，即

离散椭球窗有两个重要参数W和P，W为所选择的中心波瓣的宽度，即半带宽，要求数据窗在带宽(-W，W)内具有最大的能量，即半带宽W外具有最小频率泄露；P为时间带宽积(P=NW，N为数据窗的长度)，P决定了具有最小频率泄露的数据窗的个数K，K=2P-1(Abhey et al., 2006).图 1a 显示了时间带宽积NW为2.5时它的前四个Slepian 窗函数在时域形态，当k=0时，其形状与传统的余弦窗非常相似，而其他高阶的情况则与传统的各种时窗有很大不同，偶数阶时窗是关于中心对称的，而奇数阶不关于中心对称；最低阶k=0的权值为正，而其他各阶权值均有正有负，利用多重窗函数可以实现对有限时间序列的正交采样(陈赟等, 2005)，每个数据窗对时间序列完成不同的采样，第一个数据窗丢失的统计信息可以被第二个或其他的数据窗恢复.相对于传统的时窗分析方法而言，这种方法具有频谱光滑、畸变小，一致性好的特点.图 1(c、d、e和f中实线)显示了这四个Slepian窗函数在频域的形态，可以看出，随着k的增加，Slepian窗主瓣的宽度基本不 变，但旁瓣的幅度明显增加，即频谱的泄露增加了.

图 1 离散椭球窗在时域与频域的形态(a) 时间带宽积NW为2.5时它的前四个Slepian窗函数在时域形态；(b) 前四个多正弦窗在时域的形态； (c)、(d)、(e)和(f) 中实线和虚线分别表示椭球窗和正弦窗函数在频域的形态. Fig.1 Shapes of discrete prolate tapers in time-domain and frequency-domain (a) The four lowest-order Slepian tapers (time-bandwidth NW=2.5) in time-domain; (b) The four lowest-order sinusoidal taper in time-domain; (c)、(d)、(e) and (f) show the Slepian tapers (solid line) and sinusoidal tapers (dash line) in frequency-domain respectively.

多正弦窗是Riedel和Sidorenko (1995)提出的两种最小偏差正交窗函数之一，相对于离散椭球窗，正弦窗的个数不受频率带宽的限制，可以通过增加或减少窗的个数，以实现平滑的的多窗谱估计.

多正弦窗函数的定义为

如同Slepian窗，多正弦窗也是对数据进行正交采样.k反映了窗函数与横轴交点的个数，只有第 一个窗函数是通常意义下的数据窗，后面的窗函数用于强调窗口边缘数据的作用，更好地体现谱峰的包络(吴红卫等, 2007).

在本方法中，用一系列短的多正弦窗 h_k(n)， (n=0,1，…,N₁-1) 对一个时间域数据序列d(n)， 时间间隔为Δt，(n=0,1，…,N₂-1，N₂>N₁)移动重复加窗，直到覆盖其全部数据.对每一步加窗的数据通过傅里叶变换转换到频率域，然后将一系列子窗数据在频率域内叠加，其保留了每一个子窗数据的相位特征.其具体的公式为

这样，接收函数的表达式为

D^(k)_V(ω)、D^(k)_R(ω)和D^(k)_T(ω)分别表示径向、垂向和切向分量通过第k阶正弦窗得到的频率域形式； ^(k)_V(ω) 为D^(k)_V(ω)的共轭； N^(k)_V(ω) 为垂向分量信号前的噪声估计；K表示多正弦窗的最大阶数；α为控制高斯频谱带宽的常数.

我们分别用合成地震图和实际观测到的远震资料对本文提出的多正弦窗提取接收函数方法进行数值检验，以考察实际应用效果.

首先，我们设计一个理论的接收函数时间序列(图 2b)，并选取某个实际远震记录的垂向分量作为等效震源时间函数(图 2a)，随机噪声自然也包括其中，使之与之前的理论接收函数褶积形成径向分量(图 2c).然后分别应用本文研究的方法和频率域反褶积方法对实际记录的垂向分量和合成的径向分量作反褶积运算来提取接收函数(如图 2d所示).可以发现，多正弦窗提取接收函数方法能得到与理论值基本相同的接收函数，表明了此方法的有效性.

图 2 多正弦窗接收函数提取方法的模型检验 (a) 实际远震记录的垂向分量；(b) 理论的接收函数；(c) 合成的远震记录径向分量；(d) 本文方法(灰线)与 传统频率域( Helmberger and Wiggins (1971))方法(黑线)提取的接收函数. Fig.2 The model test of the multiple sinusoidal tapers to calculate receiver function (a) Vertical component of the field teleseismic event; (b) The theoretical reciver function; (c) Radial component of teleseismic event synthesized from (a) and (b); (d) Receiver functions estimated by the method of this paper(grey line) and Helmberger and Wiggins (1971) (black line).

我们也用全球地震台网的台站TLY(张建利等, 2012; Si et al., 2013)所记录的1991—2012年间的地震数据(图 3显示了台站TLY的位置以及记录的地震事件)对传统频率域反褶积、多椭球窗以及多正弦窗频率相关估计三种方法进行对比.三种方法在计算接收函数时，对原始数据采用同样的处理，滤波频带为0.03 Hz到 2 Hz的带通滤波，反褶积时采用的高斯系数为1.0.从对比结果看，多椭球窗和多正弦窗频率相关估计获得的接收函数明显好于传统的频率域反褶积，大部分多椭球窗和多正弦窗频率相关估计获得的接收函数相似，少部分多正弦窗频率相关估计得到的接收函数优于多椭球窗频率相关估计方法.图 4和图 5表示了两个实例，图中的a是远震事件的三分量P波波形数据，对其选取最优等效震源时间函数长度(张洪双等, 2009)，b、c、d分别表示用频率域反褶积、多椭球窗频率相关估计和多正弦窗法得到的接收函数.从图 4、5都能看出，多椭球窗和多正弦窗频率相关估计方法明显优于传统的频率域反褶积方法，在图 4中，多椭球窗和多正弦窗频率相关估计方法得到的接收函数基本相同，而图 5显示多正弦窗频率相关估计方法的结果则比多椭球窗频率相关估计更好.

图 3 台站TLY所处位置以及记录的地震事件分布(a)显示台站所在的位置，红色三角形代表台站TLY；(b)显示了地震事件相对于台站的方位，圆圈代表所记录的地震事件，圆圈大小代表震级，颜色越深代表震源越深. Fig.3 The position of the station TLY and the seismic event recorded (a) shows the position of the station TLY and the red triangle represents the station. (b) shows the positions of earthquakes, circles represent the earthquakes and their sizes show the earthquake magnitude. The darker the circle, the deeper the focus of the earthquake.

图 4 (实例一)三种方法对台站TLY实际数据提取的接收函数(a)实际数据垂向分量，南北分量和东西分量；(b)、(c)、(d)分别表示用传统频率域反褶积方法、多椭球窗以及多正弦窗方法提取的接收函数. Fig.4 (Example 1) Receiver functions from the station TLY using three methods(a) show vertical、north-south and east-west component respectively; (b)、(c)、(d) show receiver functions estimated by traditional deconvolution in frequency-domain、multiple prolate tapers and multiple sinusoidal tapers methods respectively.

图 5 (实例二)三种方法对台站TLY实际数据提取的接收函数(a)实际数据垂向分量、南北分量和东西分量；(b)、(c)、(d)分别表示用传统频率域反褶积方法、多椭球窗以及多正弦窗方法提取的接收函数. Fig.5 (Example 2) Receiver functions from the station TLY using three methods(a) show vertical, north-south and east-west component respectively; (b)、(c)、(d) show receiver functions estimated by traditional deconvolution in frequency-domain, multiple prolate tapers and multiple sinusoidal tapers methods respectively.

多正弦窗的带宽可以通过增加或减少正弦窗的个数来实现，正弦窗频率相关估计获取接收函数的质量与所采用正弦窗个数有关，适当的正弦窗个数能够取得高质量的接收函数.正弦窗个数的选择一般为3～5个，具体根据不同的数据再做调整.图 6 显示了不同的正弦窗个数对获得的接收函数的影响.

图 7展示了对TLY台站数据用三种方法获得的接收函数的叠加效果.a、b、c分别展示了传统频 率域反褶积、多椭球窗以及多正弦窗频率相关估计 三种方法得到的震中距在40°~50°，反方位角在 100°~110°内的19条接收函数，由于得到的接收函数在相同震中距及反方位角内，故具有较好的一致性，有利于比较.从这三幅图可以看出，相比于传统的频率域反褶积，多正弦窗和多椭球窗频率相关估计方法得到的接收函数更加干净，消除了低频干扰，能够得到明显的莫霍转换波震相，而且P波之前的振幅明显更小，具有更高的可靠性.我们拾取了三种方法获得的19条接收函数在莫霍面的转换波P_ms以及直达P波的振幅(接收函数中0时刻P波振幅)，得到其振幅比(图 7d)，蓝色圆圈、绿色三角以及红色五角星分别代表传统频率域、多椭球窗以及多正弦窗频率相关估计三种方法计算的接收函数的结果，无法清晰分辨出莫霍转换波震相的，没有标注在图上.从图中可以看出，多正弦窗频率相关估计获得的接收函数的转换波信号最强，传统频率域反褶积次之，多椭球窗频率相关估计最弱.对台站TLY所有的接收函数，依据IASP91模型进行时深转换后叠加(图 7e)，叠加方法采用相位加权叠加，从此图可以看出，莫霍转换波信号的强度由高到低是多正弦窗频率域相关估计、多椭球窗频率相关估计和传统频率域反褶积.综合图 7d和7e，从单个接收函数看，多椭球窗频率相关估计获得的莫霍转换波信号受到了压制，其能量低于传统频率域反褶积的结果；而从叠加结果看，多椭球频率相关估计方法的接收函数信号得到增强，其能量高于传统频率域反褶积，说明其一致性更好；但多正弦窗频率相关估计得到的接收函数莫霍转换波震相不仅信号强而且具有很好的一致性.图 7f是将图 7e中的上地幔过渡带界面信号进行放大.可以看出，其结果与莫霍面转换波表现一致.

图 6 不同的正弦窗个数对接收函数的影响 (a)实际数据的垂向分量、南北分量以及东西分量；(b) 用不同的正弦窗个数所得的接收函数. 绿色、红色和蓝色分别表示所采用的正弦窗个数为3、6以及10个. Fig.6 Effect on the receiver function using different number of the sinusoidal tapers (a) shows vertical, north-south and east-west component; (b) shows receiver function estimated using different number of the sinusoidal tapers. Green、red and blue show the number of the sinusoidal tapers is 3、6 and 10 respectively.

图 7 三种方法对台站TLY实际数据提取的接收函数叠加以及拾取莫霍转换波与直达P波的振幅比 (a)、(b)和(c)分别表示用传统频率域反褶积、多椭球窗和多正弦窗方法对震中距在40°~50°，反方位角在100°~110°内的19个地震记录计算的接收函数；(d)对三种方法计算的19条接收函数拾取的莫霍转换波与直达P波的振幅比，其中红色五角星、蓝色圆圈和绿色三角分别对应多正 弦窗、多椭球窗方法和传统频率域方法的结果； (e)所有台站TLY提取的接收函数的叠加； (f)表示对图(e)中黑色线框标记的区域进行放大. Fig.7 Stack of the receiver function from station TLY with three methods and the ratio of the converted wave at Moho to the direct P wave (a)、(b)、(c) represent 19 receiver functions within epicentral distance 40°~50°, reverse azimuth 100°~110°, estimated by traditional deconvolution in frequency-domain, multiple prolate tapers and multiple sinusoidal tapers methods respectively. (d) show the ratio of the converted wave at Moho to the direct P wave from 19 receive function estimated by three methods. Red star, blue circle and green triangle represent the results obtained by multiple sinusoidal tapers method, multiple prolate tapers method and traditional deconvolution in frequency-domain; (e) shows the stack all the receiver function from the station TLY; (f) shows the enlargement of the area delimited by black line in (e).

对三种方法计算的接收函数，我们用bootstrap 算法(Press et al.,1992; Efron and Tibshirani,1986) 抽样叠加100次，拾取每一次叠加后的410 km和660 km界面转换波P410s、P660s与直达P波的振幅比，统计结果分别展示于图 8a和图 8b.图中黑色、灰色和灰白色分别代表多正弦窗、多椭球窗以及传统频率域的结果.可以看出，多正弦窗频率相关估计得到的地幔过渡带转换波震相振幅最大，优于另外两种方法，传统频率域方法的信号最弱.

图 8 用Bootstrap方法抽样叠加100次拾取的上地幔间断面转换波与直达P波振幅比的统计结果 (a) 410 km间断面；(b) 660 km间断面.其中灰白色、灰色和黑色分别代表传统频率域反褶积、多椭球窗和多正弦窗方法的结果. Fig.8 Distribution of the ratio of converted wave at the upper discontinuities to the direct P wave, which obtained 100 stack using bootstrap method (a) 410 km discontinuity；(b) 660 km discontinuity. Off-white, gray and white represent result from traditional disconvolution in frequency-domain, multiple prolate tapers and multiple sinusoidal tapers method.

接收函数方法是研究地壳、上地幔速度界面的重要方法之一，其研究内容主要包括接收函数的提取和接收函数的反演两大部分，前者是用接收函数方法进行研究的基础.提取的接收函数越准确、精度越高，得到的研究结果越可靠.在接收函数发展的过程中，有多种方法被提出，但每一种方法都有自身的局限.时间域反褶积主要受信号中最大振幅频率成分的影响；传统的频率域反褶积方法存在“水准量”选择的困难；多椭球窗频率相关估计方法中的多椭球序列的确定选择复杂且不方便.而多正弦窗是一种最小偏差正交窗函数，其有明确的表达式，应用方便.本文提出用多正弦窗法代替多椭球窗来提取接收函数，分别用合成资料和实际资料与传统的频率域反褶积和多椭球窗频率相关估计进行对比，多正弦窗法是一种有效地提取接收函数的方法，而且得到了质量更高的接收函数，有效地提了高接收函数的分辨率，对于背景噪音较大的信号尤为明显.