2018, Vol. 61
2. 国土资源部应用地球物理重点实验室, 长春 130026;
3. 油页岩地下原位转化与钻采技术国家地方联合工程实验室, 长春 130026
2. Key Laboratory of Applied Geophysics, Changchun 130026, China;
3. National-Local Joint Engineering Laboratory of In-situ Conversion, Drilling and Exploitation Technology for Oil Shale, Changchun 130026, China
自Langston(1979)提出在等效震源假设下从远震P波可以获得台站下方地壳上地幔的脉冲响应(接收函数)后，接收函数方法发展迅速(刘启元等，2010)，国内外许多学者提出了不同的利用接收函数探测地下构造特征的方法，其中比较常用的方法有：接收函数反演地壳上地幔S波速度结构(Owens，1984；Ammon，1990；Julià，2000；Tian et al, 2014)、P波或S波接收函数CCP叠加确定地壳、岩石圈和地幔转换带结构(Yuan et al, 1997, 2006；Dueker and Sheehan, 1997；Zhu，2000)、H-κ叠加扫描地壳厚度和平均波速比(Zhu and Kanamori, 2000)以及利用接收函数Ps转换波震相研究地壳各向异性结构(Mcnamara and Owens, 1993；Savage，1998).近些年来，随着宽频带地震仪的普及，这些接收函数方法被广泛地应用于地壳上地幔结构的研究中.
H-κ叠加(Zhu and Kanamori, 2000)是接收函数中最常用的方法之一，其原理是利用莫霍面转换波(Ps)和多次波(PpPs，PpSs+PsPs)震相扫描地壳厚度H和平均波速比VP/VS，该方法可以避免人工读取转换波或多次波时差而引起的误差，并且可以通过Bootstrap重采样(Efron and Tibshirani, 1986)等方法，对扫描结果的误差进行估计，H-κ叠加方法可以为推断研究区的地壳成分和构造演化提供有效约束.由于传统的H-κ叠加方法是基于水平单层各项同性介质的假设，所以在具有复杂地壳结构地区，如地表巨厚沉积层地区、地壳各向异性较强地区、或具有倾斜莫霍面和渐变莫霍面结构地区，传统的H-κ叠加方法并不能获得准确的地壳结构信息.针对这一问题，许多学者探讨了地壳结构对H-κ叠加结果的影响，并提出了一系列改进的H-κ叠加方法.Lombardi等(2008)和Wang等(2010)通过对倾斜界面的接收函数H-κ叠加结果进行分析，得出利用界面上倾方向的方位角范围内的接收函数的叠加，可以得到更接近真实地壳结构的估计.房立华和吴建平(2009)分析了倾斜莫霍面和地壳各向异性对接收函数的影响，认为利用不同方位角的径向和切向接收函数的直达P波和Ps转换波振幅随到时随方位角的变化规律，可以有效地分辨倾斜莫霍面或地壳各向异性对叠加结果的影响.针对具有倾斜界面结构地区的接收函数H-κ叠加，张洪双等(2009)提出了同时利用径向和切向接收函数震相信息扫描地壳厚度H、平均波速比κ和倾斜角θ的H-κ-θ方法，并利用该方法对青藏高原东侧碰撞造山带地区地壳结构进行了研究.
除了复杂的莫霍面结构和地壳各向异性对H-κ扫描结果有较强影响外，沉积盆地地区地表巨厚低速沉积层对H-κ叠加也有很强的影响(罗艳等，2008；Yu et al, 2015).如图 1所示，由于沉积层内地震波速度较低，所以在沉积层底部和地表之间极易产生多次波，形成多次波混响效应，这些多次波主要为S波，并且具有较强能量，会干扰甚至覆盖来自莫霍面的转换波和多次波震相，从而影响H-κ叠加结果的准确性(Zelt and Ellis, 1999；危自根等，2016).针对沉积盆地地区地壳结构准确估计这一问题，一些学者提出对H-κ叠加方法进行改进，如Yeck等(2013)提出的双层模型H-κ叠加方法，以适应沉积盆地地区台站下方地壳结构特点，但是这种处理方法具有较大局限性，只适用于某些特定的地区的地壳结构的估计.一些学者认为，若要更准确地估计沉积盆地地区地壳结构，则需要对沉积盆地地区接收函数进行特殊处理，去除沉积层多次波的影响，因此，他们根据多次波的不同特性提出了一些消除低速沉积层多次波的方法.Langston(2011)提出利用给定的低速层模型信息进行波场延拓，从远震P波响应中分离出低速沉积层产生的上行S波场，消除低速层上行S波对接收函数的影响.Leahy等(2012)利用低速层多次波与莫霍面转换波和多次波在频谱中分布的差异，通过带通或带阻滤波器对带有混响效应的接收函数进行滤波，消除沉积层对接收函数的影响，但是这种方法的难点是滤波器的频带范围的确定.Tao等(2014a)改进了Langston(2011)的方法，通过最小化上行S波能量，从而避免了参考台站信息的使用.最近，Yu等(2015)等提出了利用接收函数归一化自相关构建共振滤波器来消除多次波对接收函数的影响，并且利用加入时差校正的H-κ叠加技术估计低速层下覆地壳厚度和平均波速比，此外，还给出了利用处理后的接收函数中的沉积层底部转换震相(PbS)和莫霍面多次波PpPs、PpSs震相扫描地壳上覆沉积层厚度和波速比的方法.
|
图 1 接收函数沉积层内多次波混响示意图 (a，c，e, g)分别为沉积层底部转换波PbS震相、莫霍面转换波Ps震相、莫霍面多次波PpPs震相和PpSs(PsPs)震相; (b，d，f, h)分别为对应的沉积层内多次波混响. Fig. 1 Schematic diagrams of near-surface reverberation in receiver function The (a, c, e, g) are the converted PbS wave from the bottom of the sediment, converted Ps wave from the Moho discontinuity, and the multiple-wave of PpPs, PpSs (PsPs) in the crust, respectively, (b, d, f, h) are the corresponding near-surface reverberations in the sedimentary layer. |
Yu等(2015)给出的低速沉积层多次波混响模型(图 1)显示，沉积层内的多次波主要是S波，并且呈很好的周期性，这与海上地震勘探中常出现的“鸣震”效应原理近似.预测反褶积方法是勘探地震中常用的去除多次波的方法之一(张军华等，2009)，由于多次波具有明显的周期性，属于可预测量，而一次波是不可预测量，预测反褶积方法可以消除可预测部分，只保留不可预测部分，即一次波，从而达到消除多次波的目的.此外，该方法具有不需要速度模型等辅助信息，参数设置简单，运行速度快等特点.本文将预测反褶积方法引入到接收函数中，消除沉积盆地地区低速沉积层对接收函数的影响，并利用改进的H-κ叠加方法(Yu et al, 2015)估计沉积盆地地区地壳结构.
1 方法 1.1 低速沉积层对接收函数的影响P波接收函数由宽频带地震仪径向分量(R)对垂向分量(Z)反褶积获得，由于其对S波速度和间断面结构较为敏感，所以可以用来对台站下方构造进行研究.利用分布在基岩地区的台站提取的接收函数，可以较为准确地估计地壳结构，然而，如果地表覆盖有巨厚的低速沉积层，其对接收函数产生明显的影响，表现为直达P波震相宽度增加，振幅减小，其后续震相为低速沉积层引起的高振幅、低频率、周期性的多次波(Zelt and Ellis, 1999；Yu et al, 2015)，这些干扰震相覆盖了直达P波以及莫霍面产生的转换波和多次波震相，从而影响对台站下方地壳、地幔结构的准确估计.
在沉积盆地地区，计算得到的接收函数除直达P波震相外，还包括沉积层底部一次转换波震相PbS、莫霍面转换波震相Ps，莫霍面多次波震相PpPs、PsPs+PpSs和沉积层内多次波震相等，而沉积层内多次波主要为S波，由于沉积层S波速度较小，射线路径近垂直于地表，所以沉积层底部的转换波震相和莫霍面转换波震相以及莫霍面多次波震相产生的多次波具有相近的周期性，其接收函数表达式为
|
(1) |
式中，R(t)为沉积盆地地区接收函数表达式，F(t)为不包含沉积层多次波的归一化接收函数，n为沉积层内多次波反射次数，Δt为沉积层内S波双程走时，An为第n次反射的多次波振幅，它与沉积层反射面的反射强度有关.当n=0时，R(t)=A0F(t)，表示接收函数不包含多次波.Yu等(2015)将公式(1)转换到频率域，并构建共振滤波器在频率域中消除接收函数多次波混响.在本研究中，我们利用预测反褶积方法在时间域中消除沉积层内多次波对接收函数的影响.
1.2 预测反褶积基本原理预测反褶积方法是由Peacock和Treitel在1969(Peacock and Treitel, 1969)年提出，该方法在地震资料处理中的得到广泛的应用，其本质为一个反滤波问题.对于任意一时间序列，其成分在可预测性上可以分为可预测量和不可预测量，凡是可以通过过去值和现在值对未来值进行准确预测的部分叫做可预测量，不可通过过去值和现在值对未来值进行预测的部分叫做不可预测量.一个时间序列之所以具有可预测性，是因为该序列各个时刻的值是相关的，反之，不可预测量各个时刻的值是互不相关的随机量.通常认为，实际观测到的时间序列既包括可预测成分，也包括不可预测成分，不可预测成分可以通过实际观测序列和该序列可预测成分之差获得，表达式为
|
(2) |
ε(t)为预测误差序列，即不可预测量，g(t)为实际观测序列，
|
(3) |
(3) 式化简为
|
(4) |
写成线性方程组形式为
|
(5) |
式中rgg(t)为实际观测序列自相关函数，m为预测滤波因子长度.其系数矩阵是一个托布里兹矩阵，该方程组可以使用莱文森递推算法快速求解(吴庆举等，2003).值得注意的是，在求解该方程组过程中需要进行预白噪化处理，即在托布里兹矩阵中主对角线上用(1+b)rgg(0)代替rgg(0), 其中b称为白噪系数，通常是一个很小的正数，在本次研究中全部取为0.001.
将式(2)写为褶积形式为
|
(6) |
式中的d(t)为反滤波因子，它由式(5)中预测滤波因子c(t)构建，表达式为
|
(7) |
将实际观测序列g(t)褶积该反滤波因子，即可求取预测误差ε(t).
本研究将该方法引入到消除沉积盆地地区的接收函数混响效应中，由于沉积层中多次波具有明显的周期性，所以认为它是可以预测部分，而用来估计地壳结构的一次波成分(图 1)，则是不可预测部分，它通过实测接收函数和其可以预测的多次波成分的差获得(公式(2)).在对实际接收函数资料进行处理时，以带有多次波混响的实际观测接收函数作为式(4)中的g(t)，ε(t)为去除沉积层多次波后的接收函数，通过公式(5)、(7)构建实际接收函数的反滤波因子d(t)，将实际接收函数g(t)和反滤波因子d(t)褶积，获得不包含多次波混响的接收函数.利用公式(5)计算接收函数预测多次波滤波因子时，预测步长a的选取非常关键，直接影响接收函数中多次波的压制效果，如果a值选取过大，则多次波成分压制的不够完全，如果a值选取过小，则会造成波形产生畸变，通常认为当预测步长等于多次波周期时，对多次波的压制效果较好(张军华等，2009).对于沉积盆地地区的接收函数，如前所述，其各个主要震相的多次波在沉积层内都近似为近垂直传播的S波，所以可以认为它们都具有相近的周期，则预测步长a应选取为沉积层内S波的双程走时.由于我们并不知道沉积层的厚度，所以难以直接计算得到沉积层内S波的双程走时，为解决这一问题，Yu等(2015)利用接收函数的归一化自相关函数来确定沉积层内S波的双程走时，并定义自相关函数的第一个波峰(通常为t=0时刻)和第一个波谷之间的时间间隔为沉积层内S波的双程走时，本文在用预测反褶积方法处理接收函数时，同样使用该方法确定预测步长a.
利用预测反褶积消除沉积盆地地区接收函数混响的处理步骤总结如下.
① 利用实测接收函数归一化自相关确定预测步长a;
② 利用公式(5)求取多次波的预测滤波因子c(t)，并根据公式(7)构建求取预测误差的反滤波因子d(t);
③ 将实测接收函数g(t)和反滤波因子d(t)进行褶积，得到消除沉积层多次波影响后的接收函数，用于估计沉积盆地地区沉积层下覆地壳厚度和平均波速比.
1.3 H-κ叠加估计含有沉积层地区地壳结构H-κ叠加(Zhu and Kanamori, 2000)方法是接收函数研究中最常用的估计地壳厚度和平均波速比的方法之一，但是利用传统的H-κ叠加方法对沉积盆地地区的地壳厚度和平均波速比进行估计时，在存在较厚沉积层的地区，会对地壳平均波速比和地壳厚度的估计产生较大影响(Yeck et al, 2013).为准确估计沉积层下覆地壳厚度和平均波速比，Yu等(2015)采用双层模型思想，消除上覆沉积层影响，准确计算沉积层下覆地壳厚度和波速比.通过分析利用共振滤波器去除多次波后的接收函数各主要震相的分布特点，Yu等(2015)给出了加入时差校正的H-κ叠加方法，估计沉积层下覆地壳厚度H和平均波速比κ，叠加公式为
|
(8) |
其中莫霍面转换波和多次波走时tPs，tPpPs，tPpSs+PsPs表达式为
|
(9) |
式中δti为第i个接收函数中PbS波相比于P波的时间延迟，可以直接从接收函数中估计，Δti为第i个接收函数沉积层内S波双程走时，或者称为预测步长，其确定方法如1.2部分所述，ri(t)为该台站第i条接收函数莫霍面转换波或多次波Ps、PpPs、PpSs+PsPs与P波的时差所对应的幅值，λ1、λ2、λ3分别为三种震相叠加的权系数，它们绝对值的和为1，N为台站对应的接收函数个数.此外，利用去除多次波后的接收函数，Yu等(2015)同时提出了估计沉积层厚度和平均波速比的方法，详见Yu等(2015)，在此不再赘述.本文将预测反褶积去除多次波混响后的接收函数，利用公式(8)进行叠加，来估计沉积盆地地区地壳厚度和平均波速比.
2 合成数据测试为验证预测反褶积方法应用于接收函数去沉积层多次波混响的有效性，我们利用合成接收函数进行了试算.我们采用和Yu等(2015)相同的带有沉积层的地壳模型，模型共分为三层，如图 2a所示，上层为巨厚沉积层，中间为地壳介质，下部为地幔介质，对各层介质分别设定P波速度、S波速度、密度和层厚度，具体参数如表 1所示.利用反射率法(Kennett, 1983, 1986)合成理论地震图，截取直达P波前10 s，后90 s数据，并利用频率域反褶积方法(Ammon，1991)提取接收函数，高斯因子选取为2.0，射线参数范围为0.04~0.08 s·km-1，间隔为0.001 s·km-1.从合成接收函数中我们可以看出(图 2b)，莫霍面转换波PS震相被沉积层多次波覆盖，同时莫霍面多次波震相也难以辨认.未经处理的接收函数H-κ扫描得到的地壳厚度和平均波速比分别为37 km和1.67(图 2c)，与实际模型存在较大偏差.
|
图 2 合成接收函数测试 (a)计算合成接收函数的模型示意图，各层参数如表 1所示; (b)依射线参数排列的合成接收函数，射线参数范围为0.04~0.08 s·km-1，间隔为0.001 s·km-1; (c)为(b)中接收函数H-κ叠加结果; (d)为(b)中接收函数归一化自相关; (e)预测反褶积处理后的接收函数; (f)为图(e)中接收函数利用改进的H-κ叠加方法(Yu et al, 2015)处理结果. Fig. 2 Test using the synthetic receiver functions (a) The calculating model for the synthetic receiving functions, each layer parameters as shown in Table 1; (b) Synthetic receiver functions plotted against epicentral distance, the ray parameters range from 0.04 s·km-1 to 0.08 s·km-1, with an interval of 0.001 s/km; (c) is the H-κ stacking result using the receiver function in (b); (d) Autocorrelations of each receiver functions in (a); (e) Resulting receiver functions after the application of the predictive deconvolution method; (f) Result of improved H-κ stacking (Yu et al, 2015) using the receiver functions in (e). |
|
表 1 上覆沉积层的地壳模型参数 Table 1 Parameters of crustal model with a sedimentary layer |
图 2d为合成接收函数归一化自相关函数，对于射线参数小于0.074 s·km-1的接收函数，归一化自相关函数中第一个波峰和第一个波谷之间时差为2.025 s，即预测步长选取为2.025 s，其他射线参数的接收函数预测步长为2.05 s，预测滤波因子长度选取为30，利用预测反褶积方法去除接收函数沉积层多次波混响后，我们可以看到比较清晰的Ps、PbS、PpPs、PpSs+PsPs震相(图 2e)，对处理后的接收函数利用公式(8)，(9)进行H-κ叠加处理，波速比κ扫描范围为1.4~1.9，间隔为0.01，地壳厚度H扫描范围为20~50 km，三种震相的权重因子λ1、λ2、λ3分别取0.5、0.3、-0.2，并利用Bootstrap重采样方法(Efron and Tibshirani, 1986)估计扫描结果的标准差，重采样次数为100次，得到的扫描结果地壳厚度为34.9±0.1 km，波速比为1.75±0.0(图 2f)，与真实值35 km，1.75基本相同，说明预测反褶积方法可以有效地去除接收函数沉积层多次波混响，并可以准确地估计沉积层下覆地壳结构.
图 3为利用本文提出的方法和Yu等(2015)的方法对图 2b中射线参数为0.06 s·km-1的合成接收函数处理结果对比图.利用原始接收函数归一化自相关函数设定共振滤波器反射系数为0.75，双程走时为2.025 s，预测反褶积方法中预测步长同样为2.025 s，预测滤波因子长度选取为30.图 3中黑色实线为利用共振滤波器滤波后的接收函数，红色实线为本文提出的预测反褶积方法处理后的接收函数，从图中可以看出，对于预测反褶积方法去除沉积层多次波混响的接收函数，直达P波、沉积层底部转换波PbS波、莫霍面转换波Ps、莫霍面多次波PpPs、PpSs+PsPs震相振幅均强于共振滤波器处理后的接收函数，说明本文提出的预测反褶积方法处理后的接收函数，更有利于利用H-κ叠加方法估计地壳结构.
|
图 3 预测反褶积方法(红)和共振滤波器方法(黑)处理结果对比图 Fig. 3 The contrast figure of predictive deconvolution method (red line) and resonance filter (black line) processing receiver function in sedimentary basin |
合成数据测试结果显示，利用预测反褶积方法可以较好地消除沉积盆地地区接收函数多次波混响，获得清晰的莫霍面转换波和多次波震相，同时利用改进的H-κ叠加方法，可以较为准确地估计沉积层下覆地壳厚度和平均波速比.位于中国东北地区的松辽盆地是典型的构造拉伸盆地，受中生代以来的西太平洋俯冲板块构造作用，该区域地壳厚度较小，平均值约为31 km左右(张广成等，2013；Zhang et al, 2014；Tao et al, 2014b；康健等，2016；朱洪翔等，2017)，并且地表具有巨厚沉积层(杨宝俊等，1996；Tao et al, 2014b；Li et al, 2016).为了测试预测反褶积方法在实际接收函数处理中的应用效果，我们利用本文提出的方法对位于松辽盆地的两个观测台站：NECESSArray台网中的NE68台站和放置在吉林大学深部实验室的CBab台站，对这两个台站的实测接收函数进行处理，并用改进的H-κ叠加方法估计两个台站下方地壳厚度和平均波速比.
NE68台站布设时间为2009年9月到2011年8月.我们选取震中距范围为30°~90°，震级大于5.8级的远震事件，截取P波前15 s，P波后85 s的三分量数据，在径向-切向-垂向(R-T-Z)坐标系中利用频率域反褶积方法(Ammon，1991)提取接收函数，高斯因子选取为2.0，并利用巴特沃斯带通滤波器进行滤波，频带范围为0.03~1 Hz，经过严格挑选，共获得88条接收函数，NE68台站位置及88个远震事件分布如图 4(a，b)所示.图 5a为最终挑选的接收函数，从图 5a中可以看出，沉积层多次波能量较强，基本覆盖了莫霍面转换波Ps和多次波PpPs、PpSs+PsPs的震相，难以直接判断莫霍面转换波和多次波震相位置，图 5c为直接利用图 5a中接收函数进行H-κ扫描的结果，其中地壳厚度为35.2 km，波速比为1.52，地壳厚度与前人的接收函数研究结果35 km(Tao et al, 2014b)和35.2 km (Yu et al, 2015)较为接近，但是波速比结果偏小，Tao等(2014b)和Yu等(2015)的波速比结果分别为1.74和1.73.
|
图 4 (a) NE68和CBab台站位置；(b) NE68台站远震事件分布；(c) CBab台站远震事件分布图 Fig. 4 (a) The location of NE68 and CBab stations; (b) Distribution of teleseismic events recorded by NE68 station; (c) Distribution of teleseismic events recorded by CBab station |
|
图 5 松辽盆地地区NE68台站实际接收函数测试结果 (a)实际提取的NE68台站接收函数; (b)预测反褶积去除沉积层多次波混响后的接收函数; (c)为(a)中接收函数H-κ叠加结果; (d)为(a)中接收函数自相关归一化; (e)各个接收函数在沉积层内的S波双程走时，由(d)中归一化自相关函数的估计得到; (f)利用改进的H-κ叠加方法处理去沉积层多次波混响后的接收函数的扫描结果. Fig. 5 Test using actual receiver functions obtaining from NE68 station in Songliao basin (a) The receiver functions from NE68 station plotted against the epicentral distance; (b) Resulting receiver functions after the application of the predictive deconvolution method; (c) The H-κ stacking result using the receiver functions in (a); (d) Autocorrelations of each receiver functions in (a); (e) The S wave two-way travel-time in sedimentary layer corresponding to each receiver functions, which is estimated by the autocorrelations in (d); (f) Result of improved H-κ stacking using the receiver functions after removing the near-surface reverberation. |
从图 5d为图 5a中接收函数的归一化自相关函数可以估计各个接收函数对应的沉积层内S波的双程走时，经计算，该台下方沉积层S波双程走时变化范围在1.125~1.35 s之间，平均值约为1.2 s(图 5e)，假设S波在沉积层内传播方向垂直于水平地表，且速度为0.7 km·s-1，则NE68台站下方沉积层厚度大约为0.42 km，这与前人的研究结果0.31 km(Tao et al, 2014b)和0.35 km(Yu et al, 2015)较为接近，说明利用归一化自相关可以较为准确地估计台站下方沉积层的S波双程走时.图 5b为去除沉积层多次波混响后的接收函数，从图中可以清晰地观察到莫霍面转换波和多次波震相，将图 5b中接收函数利用改进的H-κ叠加方法估计地壳结构，根据已有资料(傅维洲等，1998)，设定东北地区地壳平均P波速度为6.3 km·s-1，深度扫描范围为20~50 km，间隔为0.1 km，波速比扫描范围为1.4~1.9, 间隔为0.01，三种震相的权重因子λ1、λ2、λ3分别取0.5、0.3、-0.2.由公式(8)，(9)计算得到NE68台站下方地壳厚度为34±0.4 km，波速比为1.74±0.01，该结果和Tao等(2014b)和Yu等(2015)的结果较为接近，说明预测反褶积方法对沉积盆地地区实际数据也有较好的处理效果.
CBab台站布设时间为2012年到2014年.我们选取震中距范围为30°~90°，震级大于5.5级的远震事件，截取P波前10 s，P波后90 s的三分量数据，在径向-切向-垂向(R-T-Z)坐标系中利用时间域迭代反褶积方法(Ligorria and Ammon, 1999)提取接收函数，高斯因子选取为2.0，经过严格挑选，共获得77条接收函数，CBab台站位置及远震事件分布如图 4(a，c)所示.图 6a为CBab台站挑取的77条接收函数，图 6c为直接利用图 6a中接收函数进行H-κ扫描的结果，其中地壳厚度为32.2 km，波速比为1.59.图 6b为预测反褶积方法处理后的接收函数，将图 6b中接收函数利用改进的H-κ叠加方法估计台站下方地壳厚度和平均波速比，计算得到CBab台站下方地壳厚度为28.9±0.5 km，波速比为1.71±0.08.
4 讨论与结论
本文利用合成接收函数和沉积盆地地区实际接收函数对我们提出的预测反褶积方法进行了测试，结果表明对于接收函数中包含沉积层多次波混响的情况，预测反褶积方法可以有效地去除低速沉积层产生的多次波混响，获得清晰的莫霍面转换波和多次波震相，经预测反褶积方法处理后的接收函数，同时也可以观察到来自沉积层底部的转换波PbS震相，这也为利用接收函数方法估计沉积层结构提供了可能.相比于接收函数中其他去除沉积层多次波的方法，本文提出的预测反褶积方法具有以下特点：
(1) 只需要根据实测接收函数确定沉积层S波双程走时，即预测步长这一个参数.利用共振滤波器方法(Yu et al, 2015)对沉积盆地地区接收函数进行处理时，需要从实际接收函数归一化自相关函数中同时确定滤波器所需要的沉积层反射系数和S波双程走时，由于实际数据往往具有较强的噪声干扰，同时确定两个滤波器参数时，增大了误差产生的可能性，从而影响接收函数沉积层多次波的滤除效果.而本文提出的预测反褶积方法避免了对沉积层反射系数的估计，降低了参数的估计产生误差的可能性.
(2) 计算效率高，参数确定简单，适用于处理大批量数据.本文提出的方法只需要在时间域对接收函数进行处理，与Yu等(2015)的方法相比，避免了对接收函数在时间域和频率域的相互转换，计算稳定性更好，此外，相比于Langston(2011)和Tao等(2014a)的波场延拓的方法，原理更为简单.
(3) 预测反褶积方法去除沉积层多次波混响后接收函数各个震相比较清晰，在H-κ叠加方法中增强了对沉积层下覆地壳厚度和平均波速比估计的约束.
(4) 预测反褶积方法实质为共振滤波器方法的时间域实现，所以两种方法去除沉积层多次波混响原理一致.本文提出的预测反褶积对短周期的多次波消除能力较强，但是对于沉积层较厚地区并不适用，这可能与不同沉积层的厚度的多次波产生的机制不同有关，针对这一问题，还需要后续更细致的研究.
致谢 特别感谢两位评审专家对本文提出的修改完善意见，本文使用的NECESSArray台网部分地震波形数据由IRIS数据中心提供.感谢中国科学院地质与地球物理研究所的田小波研究员对本论文提出了建设性修改意见.
Ammon C J. 1991. Isolation of receiver effects from teleseismic P waveforms. Bulletin of the Seismological Society of America, 81(6): 2504-2510. |
Ammon C J, Randall G E, Zandt G. 1990. On the nonuniqueness of receiver function inversions. Journal of Geophysical Research:Solid Earth, 95(B10): 15303-15318. DOI:10.1029/JB095iB10p15303 |
Dueker K G, Sheehan A F. 1997. Mantle discontinuity structure from midpoint stacks of converted P to S waves across the Yellowstone hotspot track. Journal of Geophysical Research:Solid Earth, 102(B4): 8313-8327. DOI:10.1029/96JB03857 |
Fang L H, Wu J P. 2009. Effects of dipping boundaries and anisotropic media on receiver functions. Progress in Geophysics (in Chinese), 24(1): 42-50. |
Fu W Z, Yang B J, Liu C, et al. 1998. Study on the seismology in Manzhouli-Suifenhe geoscience transect of China. Journal of Changchun University of Science and Technology (in Chinese), 28(2): 206-212. |
Efron B, Tibshirani R. 1986. Bootstrap methods for standard errors, confidence intervals, and other measures of statistical accuracy. Statistical Science, 1(1): 75-77. DOI:10.1214/ss/1177013816 |
Julià J, Ammon C J, Herrmann R B, et al. 2000. Joint inversion of receiver function and surface wave dispersion observations. Geophysical Journal International, 143(1): 99-112. DOI:10.1046/j.1365-246x.2000.00217.x |
Kang J, Wei Q H, Zhou L, et al. 2016. Deep structure based on seismic array observations in Daqing Area. Journal of Jilin University (Earth Science Edition) (in Chinese), 46(3): 900-910. |
Kennett B L N. 1983. Seismic Wave Propagation in Stratified Media. Cambridge: Cambridge University Press.
|
Kennett B L N. 1986. Wavenumber and wavetype coupling in laterally heterogeneous media. Geophysical Journal International, 87(2): 313-331. DOI:10.1111/j.1365-246X.1986.tb06626.x |
Langston C A. 1979. Structure under mount rainier, Washington, inferred from teleseismic body waves. Journal of Geophysical Research:Solid Earth, 84(B9): 4749-4762. DOI:10.1029/JB084iB09p04749 |
Langston C A. 2011. Wave-field continuation and decomposition for passive seismic imaging under deep unconsolidated sediments. Bulletin of the Seismological Society of America, 101(5): 2176-2190. DOI:10.1785/0120100299 |
Leahy G M, Saltzer R L, Schmedes J. 2012. Imaging the shallow crust with teleseismic receiver functions. Geophysical Journal International, 191(2): 627-636. DOI:10.1111/gji.2012.191.issue-2 |
Li G L, Chen H C, Niu F L, et al. 2016. Measurement of Rayleigh wave ellipticity and its application to the joint inversion of high-resolution S wave velocity structure beneath northeast China. Geophysical Research:Solid Earth, 121(2): 864-880. DOI:10.1002/2015JB012459 |
Ligorria J P, Ammon C J. 1999. Iterative deconvolution and receiver-function estimation. Bulletin of the Seismological Society of America, 89(5): 1395-1400. |
Liu Q Y, Li Y, Chen J H, et al. 2010. Joint inversion of receiver function and ambient noise based on Bayesian theory. Chinese Journal of Geophysics (in Chinese), 53(11): 2603-2612. DOI:10.3969/j.issn.0001-5733.2010.11.008 |
Lombardi D, Braunmiller J, Kissling E, et al. 2008. Moho depth and Poisson's ratio in the Western-Central Alps from receiver functions. Geophysical Journal International, 173(1): 249-264. DOI:10.1111/gji.2008.173.issue-1 |
Luo Y, Chong J J, Ni S D, et al. 2008. Moho depth and sedimentary thickness in Capital region. Chinese Journal of Geophysics (in Chinese), 51(4): 1135-1145. |
Mcnamara D E, Owens T J. 1993. Azimuthal shear wave velocity anisotropy in the basin and range province using Moho Ps converted phases. Journal of Geophysical Research:Atmospheres, 98(B7): 12003-12017. DOI:10.1029/93JB00711 |
Owens T J, Zandt G, Taylor S R. 1984. Seismic evidence for an ancient rift beneath the Cumberland Plateau, Tennessee:A detailed analysis of broadband teleseismic P waveforms. Journal of Geophysical Research:Solid Earth, 89(B9): 7783-7795. DOI:10.1029/JB089iB09p07783 |
Peacock K L, Treitel S. 1969. Predictive deconvolution:Theory and practice. Geophysics, 34(2): 155-169. DOI:10.1190/1.1440003 |
Savage M K. 1998. Lower crustal anisotropy or dipping boundaries? Effects on receiver functions and a case study in New Zealand. Journal of Geophysical Research:Solid Earth, 103(B7): 15069-15087. DOI:10.1029/98JB00795 |
Tao K, Liu T Z, Ning J Y, et al. 2014a. Estimating sedimentary and crustal structure using wavefield continuation:Theory, techniques and applications. Geophysical Journal International, 197(1): 443-457. DOI:10.1093/gji/ggt515 |
Tao K, Niu F L, Ning J Y, et al. 2014b. Crustal structure beneath NE China imaged by NECESSArray receiver function data. Earth and Planetary Science Letters, 398: 48-57. DOI:10.1016/j.epsl.2014.04.043 |
Tian X B, Liu Z, Si S K, et al. 2014. The crustal thickness of NE Tibet and its implication for crustal shortening. Tectonophysics, 634: 198-207. DOI:10.1016/j.tecto.2014.07.001 |
Wang P, Wang L S, Mi N, et al. 2010. Crustal thickness and average VP/VS ratio variations in southwest Yunnan, China, from teleseismic receiver functions. Journal of Geophysical Research:Solid Earth, 115(B11): B11308. DOI:10.1029/2009JB006651 |
Wei Z G, Chu R S, Chen L, et al. 2016. Analysis of H-κ stacking of receiver functions beneath crust with complex structure:taking the Anatolia Plate as an example. Chinese Journal of Geophysics (in Chinese), 59(11): 4048-4062. DOI:10.6038/cjg20161110 |
Wu Q J, Tian X B, Zhang N L, et al. 2003. Receiver function estimated by Wiener filtering. Earthquake Research in China (in Chinese), 19(1): 41-47. |
Yang B J, Mu S M, Jin X, et al. 1996. Synthesized study on the geophysics of Manzhouli-Suifenhe geoscience transect, China. Chinese Journal of Geophysics (in Chinese), 39(6): 772-782. |
Yeck W L, Sheehan A F, Schulte-Pelkum V. 2013. Sequential H-stacking to obtain accurate crustal thicknesses beneath sedimentary basins. Bulletin of the Seismological Society of America, 103(3): 2142-2150. DOI:10.1785/0120120290 |
Yu Y Q, Song J G, Liu K H, et al. 2015. Determining crustal structure beneath seismic stations overlying a low-velocity sedimentary layer using receiver functions. Journal of Geophysical Research:Solid Earth, 120(5): 3208-3218. DOI:10.1002/2014JB011610 |
Yuan X H, Kind R, Li X Q, et al. 2006. The S receiver functions:Synthetics and data example. Geophysical Journal International, 165(2): 555-564. DOI:10.1111/gji.2006.165.issue-2 |
Yuan X H, Ni J, Kind R, et al. 1997. Lithospheric and upper mantle structure of southern Tibet from a seismological passive source experiment. Journal of Geophysical Research:Solid Earth, 102(12): 27491-27500. |
Zelt B C, Ellis R M. 1999. Receiver-function studies in the Trans-Hudson Orogen, Saskatchewan. Canadian Journal of Earth Sciences, 36(4): 585-603. DOI:10.1139/e98-109 |
Zhang G C, Wu Q J, Pan J T, et al. 2013. Study of crustal structure and Poisson ratio of NE China by H-κ stack and CCP stack methods. Chinese Journal of Geophysics (in Chinese), 56(12): 4084-4094. DOI:10.6038/cjg20131213 |
Zhang H S, Tian X B, Teng J W. 2009. Estimation of crustal VP/VS with dipping Moho from receiver functions. Chinese Journal of Geophysics (in Chinese), 52(5): 1243-1252. DOI:10.3969/j.issn.0001-5733.2009.05.013 |
Zhang J H, Miao Y S, Zheng X G, et al. 2009. Discussion of several theoretical questions to remove seismic multiples using predictive deconvolution. Computing Techniques for Geophysical and Geochemical Exploration (in Chinese), 31(1): 6-10. |
Zhang R Q, Wu Q J, Sun L, et al. 2014. Crustal and lithospheric structure of Northeast China from S-wave receiver functions. Earth and Planetary Science Letters, 401: 196-205. DOI:10.1016/j.epsl.2014.06.017 |
Zhu H X, Tian Y, Liu C, et al. 2017. High-resolution crustal structure of Northeast China revealed by teleseismic receiver function. Chinese Journal of Geophysics (in Chinese), 60(5): 1676-1689. DOI:10.6038/cjg20170506 |
Zhu L P. 2000. Crustal structure across the San Andreas Fault, southern California from teleseismic converted waves. Earth and Planetary Science Letters, 179(1): 183-190. DOI:10.1016/S0012-821X(00)00101-1 |
Zhu L P, Kanamori H. 2000. Moho depth variation in southernCalifornia from teleseismic receiver functions. Journal of Geophysical Research:Solid Earth, 105(B2): 2969-2980. DOI:10.1029/1999JB900322 |
房立华, 吴建平. 2009. 倾斜界面和各向异性介质对接收函数的影响. 地球物理学进展, 24(1): 42-50. |
傅维洲, 杨宝俊, 刘财, 等. 1998. 中国满洲里-绥芬河地学断面地震学研究. 长春科技大学学报, 28(2): 206-212. |
康健, 韦庆海, 周琳, 等. 2016. 利用地震台阵观测资料研究大庆地区深部构造. 吉林大学学报(地球科学版), 46(3): 900-910. |
刘启元, 李昱, 陈九辉, 等. 2010. 基于贝叶斯理论的接收函数与环境噪声联合反演. 地球物理学报, 53(11): 2603-2612. DOI:10.3969/j.issn.0001-5733.2010.11.008 |
罗艳, 崇加军, 倪四道, 等. 2008. 首都圈地区莫霍面起伏及沉积层厚度. 地球物理学报, 51(4): 1135-1145. |
危自根, 储日升, 陈凌, 等. 2016. 复杂地壳接收函数H-κ叠加——以安纳托利亚板块为例. 地球物理学报, 59(11): 4048-4062. |
吴庆举, 田小波, 张乃铃, 等. 2003. 用Wiener滤波方法提取台站接收函数. 中国地震, 19(1): 41-47. DOI:10.3969/j.issn.1001-4683.2003.01.005 |
杨宝俊, 穆石敏, 金旭, 等. 1996. 中国满洲里——绥芬河地学断面地球物理综合研究. 地球物理学报, 39(6): 772-782. |
张广成, 吴庆举, 潘佳铁, 等. 2013. 利用H-K叠加方法和CCP叠加方法研究中国东北地区地壳结构与泊松比. 地球物理学报, 56(12): 4084-4094. DOI:10.6038/cjg20131213 |
张洪双, 田小波, 滕吉文. 2009. 接收函数方法估计Moho倾斜地区的地壳速度比. 地球物理学报, 52(5): 1243-1252. DOI:10.3969/j.issn.0001-5733.2009.05.013 |
张军华, 缪彦舒, 郑旭刚, 等. 2009. 预测反褶积去多次波几个理论问题探讨. 物探化探计算技术, 31(1): 6-10. DOI:10.3969/j.issn.1001-1749.2009.01.003 |
朱洪翔, 田有, 刘财, 等. 2017. 中国东北地区高分辨率地壳结构:远震接收函数. 地球物理学报, 60(5): 1676-1689. DOI:10.6038/cjg20170506 |