2018, Vol. 61
2. 中国科学院测量与地球物理研究所大地测量与地球动力学国家重点实验室, 武汉 430077;
3. 中国地震局地球物理勘探中心, 郑州 450002
2. State Key Laboratory of Geodesy and Earth's Dynamics, Institute of Geodesy and Geophysics, Chinese Academy of Sciences, Wuhan 430077, China;
3. Geophysical Exploration Center, China Earthquake Administration, Zhengzhou 450002, China
从远震三分量波形记录提取的接收函数，有效压制了震源复杂性，为地壳上地幔结构研究提供了关键信息.接收函数中各震相的到时以及幅值信息对台站下方的速度结构非常敏感(Langston, 1977, 1979；吴庆举等，2007；刘启元和Kind，2004；许卫卫和郑天愉，2002；贺传松等，2003)，具有良好的垂向分辨率，是分析地球内部结构和揭示其动力学过程的重要手段，被广泛地应用于Moho、LAB(lithosphere- asthenosphere boundary)、上下地幔过渡带间断面(410 km、660 km界面)的研究中.接收函数的计算大多基于地表台站的波形资料，近年来数量日益增多的井下地震计也为此类研究提供了新的数据来源.
基于地表台站的接收函数在Moho面、岩石圈结构的研究中发挥了重要作用.Zhu和Kanamori(2000)从南加州地震台的地震记录中提取远震P波接收函数，并在H-K域叠加获得了该区域Moho面的分布.Chen等(2006)利用基于远震P波接收函数的波动方程叠后偏移成像方法，研究了中国东部郯庐断裂带的岩石圈结构.张广成等(2013)利用中国东北布设的116个流动地震台站以及121个固定宽频台站(国家台网和区域台网)所记录的824个远震事件数据，采用远震P波接收函数CCP叠加和H-K叠加两种方法，获得了研究区详尽的地壳厚度图像.在LAB的研究中，接收函数相较于地震面波方法具有更高的横向分辨率，特别是对于S波接收函数，其具有不受间断面多次反射震相干扰的特点.Chen等(2008)提取华北地区东北部S波接收函数，利用波动方程偏移成像方法研究了该区域岩石圈厚度的横向变化.在上下地幔过渡带界面(410 km、660 km界面)的研究中，地表台站的接收函数能够精确描绘界面的横向变化特征，为研究地幔温度分布和变化，以及地幔对流提供了重要依据.Shen等(1996)通过接收函数研究了410 km和660 km的不连续面及其横向变化，揭示了冰岛下方以及毗邻的大西洋中脊处400~700 km深度上涌的地幔流.
近年来井下地震观测取得了长足进展，其数据也可用于接收函数计算.井下地震记录受地表噪声干扰小、数据信噪比高、可提高对地震资料的利用率.例如：Dreger等利用Hayward fault台网井下地震计记录的高质量地震信号研究M3+的微小地震破裂过程(Taira et al., 2015).对于有沉积层的地区，地表的地震数据往往会受到沉积盆地放大效应影响，而井下地震计通常布置在几百米深度，甚至在基岩面以下，基本不受盆地的影响，可以近似作为基岩记录.例如Fukushima等利用该特性，将日本KiK-net的井下地震数据近似作为基岩记录，对沉积盆地地表的地震记录作反卷积来研究沉积盆地的放大效应(Fukushima et al., 2016).以前虽然井下地震计较少，但仍在地震定位、地震事件监测中发挥了重要作用.现在，井下地震观测范围有了极大的扩展，特别是在地震重点区(例如我国的首都圈地区，日本的Hi-Net地震台网等)，为了更好地监测地震活动性，布设了大量的地震仪，获得了大量的井下地震记录(Okada et al., 2004).首都圈地区(包括北京市、天津市及河北省)是我国防震减灾重点示范区，为了加强对首都圈地区的地震监测，提高地震应急快速反应能力，从1999年到2001年，建设了实时传输的首都圈测震台网(http://www.csndmc.ac.cn/network/base/area/bases/main).该台网由107个台站组成，台站之间的距离约50 km.其中53个台站配备有井下地震计，埋深最浅的为赵县台(ZHX，150 m)，埋深最大的为唐海台(TAH，480 m).王林瑛等(2008)综合利用首都圈地区地表和井下台站波形数据，研究了文安地震前后该地区的波速比(VP/VS)分布特征的变化情况.Wu等(2008)、高原和吴晶(2008)通过地表和井下台站波形数据分析了首都圈地区地壳中剪切波速度的横向不均匀性.
然而，除上述优点外，井下地震计不仅能观测到与地表地震记录相同的震相，还能观测到较强的地表反射波(图 1a).地表是自由边界，在井下地震记录垂向分量上，反射P波和入射P波振幅相近、相位相反，会在特定频率上干涉相消，使得垂向分量的振幅谱出现能量极小，甚至接近于零的频谱极小(spectrum null)现象.这种频谱极小现象在Rayleigh波的振幅谱中已有研究(Tsai and Aki, 1970)，但本文则着重于研究井下地震计波形中P波的频谱极小现象.在已有的正演接收函数方法中，多数针对地表接收台站，基本不适用于井下地震计的情形.本文首先展示首都圈地区井下地震计波形垂向分量的频谱极小现象，随后，结合首都圈地区已有的地壳模型，发展了一套井下地震计接收函数的正演计算方法.并在此基础上解释频谱极小现象，研究频谱极小现象对井下地震计接收函数稳定性的影响.
|
图 1 射线路径示意图及白家疃台和唐海台波形垂向分量振幅谱观测比较图 (a)地表及井下接收台站地震射线示意图；(b)三次地震事件，在白家疃台(BJT)和唐海台(TAH)地震波形垂向分量；(c, d, e) 2013年5月24日5时44分Okhotsk海(54.892°N, 153.221°E)MW8.3地震、2013年7月7日18时35分New Ireland附近(51.849°N, 178.735°E)MW7.9地震、2014年6月23日20时53分Aleutian Island附近(51.849°N, 178.735°E)MW7.3地震在白家疃地表台(BJT，上)及唐海井下台(TAH，下)地震波形垂向分量平滑前(红色)及平滑后(黑色)的振幅谱；(f)三次地震事件，白家疃台(BJT，左)及唐海台(TAH，右)地震波形垂向分量振幅谱平滑结果; (g)三次地震事件，白家疃台对唐海台地震波形垂向分量振幅谱平滑结果的比值(BJT/TAH).其中红色曲线、绿色曲线、蓝色曲线依次为Okhotsk Sea地震、New Ireland附近地震、Aleutian Island附近地震的结果；箭头指示了spectrum null出现的频率. Fig. 1 Ray paths, amplitude spectrum of seismograms on vetical component at BJT and TAH station and its comparison (a) Ray paths of the receiver at surface and in borehole; (b, c, d) The amplitude spectrum before smoothing (red line) and after smoothing (black line) from seismograms on vetical component at BJT station on surface (up) and TAH station in borehole (down), which are from the earthquakes occuring at 05 : 44 hours on 24 May 2013 in Okhotsk Sea (MW8.3), at 18 : 35 hours on 07 July 2013 near the New Ireland (MW7.9) and at 05 : 44 hours on 23 June 2013 near the Aleutian Island (MW7.3) separately; (e) Vertical component of seismogram from three events at BJT station (left) and TAH station (right); (f) The result of amplitude spectrum being smoothed from three events at BJT station (left) and TAH station (right); (g) The ratio of amplitude spectrum being smoothed at BJT station and that at TAH station (BJT/TAH).Among them, the red, green and blue curve are the results of the Okhotsk Sea earthquake, the earthquake near New Ireland and the earthquake near Aleutian Island separately.The arrows indicate the frequency at which spectrum null appears. |
基岩地表台的波形基本反映了震源的频谱信息，将其振幅谱与井下地震计数据的振幅谱进行对比，可以研究波形被改造的情况.这里选取了相同地震事件，同属首都圈地区的地表白家疃台(BJT)和井下地震计埋深最大的唐海台(TAH，480 m)的地震数据进行对比(郑秀芬等，2009).深源地震相比于浅源地震干扰较小，在研究中，选取了震源较深的以下三个地震事件(表 1)，得到白家疃台(BJT)和唐海台(TAH)波形数据垂向分量的振幅谱，并将结果进行平滑(图 1(b、c、d)).由图可见，在白家疃台(BJT)各地震事件波形数据垂向分量的振幅谱中，随着频率的增高，能量逐渐降低.而在唐海台(TAH)的结果中，除上述特点外，还在1.1 Hz、3.3 Hz及5.5 Hz的附近频段内出现能量极小的频谱极小现象.该结果说明这种频谱极小现象不是由震源造成的，而是由于接收台站处的结构造成的.而接收函数主要是通过径向分量和垂向分量反卷积得到，频谱极小会造成反卷积不稳定的情况.提取接收函数可以通过时间域反卷积和频率域反卷积两类方法计算获得，其中频率域water-level反卷积方法是通过引入water-level (水准量)的方式来避免频谱极小造成的反卷积不稳定，以提取稳定的接收函数(张洪双等，2009).
|
表 1 观测地震事件表 Table 1 Catlog of events |
井下地震计接收函数的正演包括模型的选取(本文选取首都圈地区模型)、理论波形的正演计算、从理论地震波形中利用反卷积提取接收函数三个部分.
2.1 首都圈地区模型根据前人的研究成果(Zheng et al., 2005；Chong and Ni, 2009；罗艳等，2008)构建了首都圈地区地壳速度模型如图 2所示.该地区上地壳由沉积盖层和其下的结晶基底两部分组成.沉积盖层的厚度为4.9 km，顶部为新生代低速沉积层，顶部至结晶基底为梯度速度层.在结晶基底之下，P波、S波速度均为线性增加.Moho面深度为35 km，Moho面之下是均匀半空间.
|
图 2 首都圈地区地壳模型图 实线代表S波速度，标有VP的虚线代表P波速度，标有density的点画线代表密度分布曲线. Fig. 2 The crustal model of captial area The solid line represents the S wave velocity, the dashed line labeled VP represents the P wave velocity, the dashed-dotted line labeled density represents the density distribution curve. |
井下地震计往往架设在沉积盆地地区，浅层的地震波速横向变化不大，具有较好的一维结构特征，对于井下地震计波形的正演，可以使用水平分层模型.在水平分层介质中，常用的波形正演方法有“广义射线法”、“反射率法”、“全波形理论”、“WKBJ”以及“广义反透射系数法”等(Helmberger，1968；Fuchs and Muller，1971；Richards, 1973, 1976；Chapman, 1976, 1978；Kennett and Kerry, 1979；Kennett，1980；Kennett and Clarke, 1983；Chen，1993).其中Kennett方法能够计算井下理论地震图，但是无法计算井下地震波上下行波的各分量，而上行波分量的计算可为后续研究压制频谱极小的波形分解方法提供基础.Tao等(2014)给出了由自由地表波形记录计算地下波形的方法.本文参考Tao等(2014)的推导方法，从波动方程开始推导，得出同时计算井下理论地震波及上下行波分量的公式.
本文选取入射波只有P波的P-SV系统进行计算分析远震P波接收函数.模型由n个水平层和半空间组成，第一层的上表面为自由表面(图 3).在推导中，坐标轴设定与Aki和Richards(2002)中的规定一致.对于水平分层介质，每一层均满足波动方程(Haskell，1960、1962)：
|
图 3 分层结构及半空间模型图 Fig. 3 Layered structure and half space model |
|
其中
每一层的波动方程，可表示成位移、应力向量的微分形式(Aki and Richards, 2002).
|
(1) |
其中U(j)，V(j)表示水平方向和垂直方向的谐波展开位移系数；P(j)，S(j)表示水平方向和垂直方向的谐波展开应力系数.(1)式可归纳为
|
(2) |
其中，f(j)(z)表示位移应力向量；M(j)是一个4×4的系数矩阵.在地表采用了自由边界条件，即地表的应力分量均为0；最后一层采用均匀半空间.根据Aki和Richards的分析(Aki and Richards, 2002)，微分方程的解为
|
(3) |
其中
|
|
为了方便表示，(3)式可以写为如下形式：
|
(4) |
这里，E11(j)，E12(j)，E21(j)和E22(j)是E(j)的2×2的子矩阵，
|
对于前n层介质，任意一层介质的上下界面由方程(2)的解(4)知
|
(5) |
|
(6) |
比较(5)和(6)，可得到如下关系
|
(7) |
其中，th(j)=z(j)-z(j-1)表示任意一层介质的厚度.
为了简化，可令
|
则有
|
(8) |
对于整个模型内部任意两层介质的界面，由连续边界条件知
|
其中，
|
(9) |
对于地表z(0)处，根据波动方程有
|
(10) |
其中
|
(11) |
根据地表z(0)处自由边界条件，应力分量为0，有
|
(12) |
式(11)(12)代入(10)，有
|
(13) |
由于只考虑P波入射的情况，故有
|
(14) |
联立式(3)、(8)、(9)、(13)、(14)可得
|
(15) |
|
(16) |
(15)式中，E(1)A(1)(E(1))-1… E(n)A(n)(E(n))-1E(n+1)是4×4的矩阵，
令
|
(17) |
可得
|
(18) |
由此，地表的位移和井下各界面处的位移应力向量均可以解得.
以此为基础，通过传播矩阵对于井下地震计，给定深度后的位移可由下式获得.
|
其中c, d由
在计算地表及井下各深度理论地震波形后，可用反卷积方法提取接收函数.Oldenburg对反卷积方法进行了详细的讨论(Oldenburg，1981)，可以分为频率域和时间域两大类方法.其中频率域比较常见的反卷积方法有Clayton和Wiggins等采用的水准量(water-level)反卷积方法，刘启元等采用的最大似然原理的复谱比法(Helmberger and Wiggins, 1971；Clayton and Wiggins, 1976；刘启元等，1996)；时间域的方法有Gurrola等利用多道滤波的维纳(winner)滤波反卷积方法，Ammon等的迭代反卷积方法和吴庆举等的最大熵法反褶积方法等(Gurrola et al., 1995；Ligorría and Ammon, 1999；吴庆举等，2003).本文选取了频率域的水准量反卷积方法和时间域的迭代反卷积方法同时进行接收函数的计算，并将两种方法的结果进行对比，以此来验证结果的稳定性.
2.3 首都圈地区模型的井下地震计理论接收函数利用首都圈地区模型，计算P波平面波入射时井下地震计的理论波形.首都圈地区井下地震台站最大埋深约480 m，因此本文研究的深度主要集中在0.5 km内.在这里分别计算0~0.5 km，每相隔0.1 km处的理论地震图.对于震源时间函数，由于地下接收点埋深较浅，入射波震相与地表反射波震相到时相差较小，为便于区分，使用持续时间为0.1 s的三角波，得到了理论地震图如图 4所示.
|
图 4 首都圈模型下地表及井下接收点的理论地震图 (a)径向分量；(b)垂向分量. Fig. 4 The synthetic seismogram of receiver at surface and in borehole with capital area crustal model (a) Radial component; (b) Vertical component. |
由图可知，对于径向分量，在地表的接收点记录到比较明显的震相依次有：直达p波、4.9 km结晶基底界面产生的转换波pbs、4.9 km界面处反射的一重多次波ppPbs、Moho面产生的转换Pms波.地表理论地震波形上出现的震相，在井下地震计理论波形上则以一组震相的形式出现，同一组内各震相之间到时差随井下地震计埋深的增大而变大.例如，直达p波对应上行P1波和下行Pfp波和Pfs波，转换波Pbs对应上行S1波和下行Sfs波.多次波ppPbs和转换Pms波与转换波Pbs情况类似.对于垂向分量，在地表接收点比较明显的震相依次是直达p波、4.9 km界面产生的一重多次波pPp、二重多次波pPpPp，而没有S波震相.对于井下地震计，直达p波对应上行P1波和Pfp波，特别是P1波的能量相对理论波形上其他震相较强，且由于自由地表是强反射界面，使得Pfp波和P1波振幅相近，偏振方向一致.
对于垂向分量，其主要震相是入射P波和地表反射的P波，两者振幅相近，但存在一定的相位差，因此会在某些频率上出现干涉相消，即表现为振幅谱上出现频谱极小现象，该现象可能对接收函数造成影响.
以上分析表明计算得到的理论地震图是可靠的.在得到地表及地下不同深度处的理论地震波形后，可以计算其接收函数.在反卷积提取接收函数过程中，对于频率域的水准量方法，选择的时间窗口为从理论地震图的0 s开始，至P波之后80 s的窗口，该窗口可完整覆盖模型介质中两个速度界面产生的ps震相，以及相应的多次波(Langston, 1977, 1979).设置水准因子为0.01，高斯因子选择2.5.对于时域迭代反卷积方法，选择的时间窗口和高斯因子与频率域的水准量方法相同，迭代次数5000次，该迭代次数可以保证后续每迭代一次的结果与前一次结果的变化率不超0.01%(Ligorría and Ammon, 1999).最终得到地表及不同深度接收点的理论接收函数(图 5).由图可知，当高斯因子为2.5(约1.25 Hz低通滤波)时，不同深度接收函数的特征和前面各深度理论地震波形特征类似.地表接收函数中的每个震相，在井下地震计接收函数中会分裂成一组震相，使其更加复杂化.对同一接收点，地表、0.1 km和0.2 km深度处的两种反卷积方法提取的接收函数的互相关系数都在99%以上，可以认为接收函数的结果是稳定的.然而在0.3 km、0.4 km和0.5 km深度处，两种方法提取的接收函数的互相关系数开始下降，且接收点深度越深差异越大.这种互相关系数降低的现象不是由数值误差引起的，而是由于接收函数结果不稳定造成的，且不稳定情况随接收点深度增大而加剧.
|
图 5 不同深度的时间域及频率域方法得到的接收函数(高斯因子2.5) 黑色表示频率域水准量water-level方法结果，灰色表示时间域迭代反卷积方法结果. Fig. 5 The receiver functions from time domain and frequency domain method at different depths (Gaussian factor 2.5) Black lines indicate the result of the frequency domain water- level method, and gray lines indicate the result of the time domain iterative deconvolution method. |
在探究井下地震计接收函数随着接收点深度增大，不稳定性加剧的原因之前，先来解释井下地震计数据垂向分量的频谱极小现象.分别对地表及不同深度接收点的理论地震波形的垂向分量做傅里叶变换，获得其振幅谱，并与震源时间函数的振幅谱进行比较(图 6).由图可见，各振幅谱总体特征相似，低频能量最大，随着频率升高，能量减弱.这说明各深度的理论地震波形总体保持了震源时间函数的特征.但除地表之外，地下各接收点处的接收函数振幅谱在一定频率处出现振幅能量的极小，即频谱极小现象；且频谱极小现象呈周期性变化特征，随着接收点深度的加深，周期变小，频谱极小开始出现的频率值变低.对此有如下解释：对于自由地表，当入射P波的入射角度较小时，P1与Pfp的振幅基本相等，且P1与Pfp的偏振方向相同，相位不同.当相位差为π的奇数倍时，P1与Pfp干涉相消，出现频谱极小现象.其具体的频率值推导如下：
|
图 6 震源时间函数及不同深度的理论地震波形垂向分量的振幅谱 Fig. 6 amplitude spectrum on vertical component from source time function and synthetic seismogram at different depth |
当P波穿过一个厚度为l的层后，相位变化为
|
其中，ω=2πf，i是入射角，且有
|
需要说明的是，此公式计算的结果主要受模型介质的速度变化的影响，而对射线参数不敏感，具体将在后面的讨论中说明.对于首都圈模型，射线参数为0.06 s·km-1时在不同深度，频谱极小出现的频率公式见表 2.
|
|
表 2 频率公式 Table 2 The formulas for frequency |
特别地，h=0.48 km时的结果与在唐海台(TAH)观测到的垂向分量波形振幅谱在1.1 Hz、3.3 Hz、5.5 Hz附近能量急剧减小的特征一致(图 1).
3.2 井下地震计接收函数的稳定性分析在了解井下地震计的波形垂向特征后，可解释井下地震计接收函数结果的不稳定情况随接收点深度增大而加剧现象.在做出了所有接收函数的振幅谱(图 7)后可知，与接收函数结果在时域的互相关结果(图 5)类似，地表、0.1 km和0.2 km深度处，两种方法提取的接收函数的振幅谱一致性较好，而0.3 km、0.4 km和0.5 km深度处接收函数的振幅谱大体上比较一致，只是在部分频段上不一致，具体为：在0.3 km处，振幅谱在1.64 Hz附近不同，在0.4 km处，振幅谱在1.25 Hz附近不同，在0.5 km处，振幅谱在1.01 Hz附近不同.由此知，两种方法接收函数结果的互相关系数随接收点深度降低是由于在某些频段内能量不同造成的，且这些频段随着接收点深度的增大向低频端移动，能量的差异随接收点深度而增大.此外，不同深度处能量差异频段的中心点与其理论地震波形垂向分量的频谱极小初次出现的频率一一对应，说明这种现象是地震波垂向分量的频谱极小现象造成的.对此可以这样解释：频率域和时间域的反卷积，接收函数可以理解为径向分量与垂向分量在频率域的商，只是反卷积的实现方式不同.在频谱极小出现的频段内，作为分母的垂向分量的能量快速减小，远小于作为分子的径向分量能量，使比值的结果不稳定，从而使得两种接收函数在频谱极小出现的频段内的能量不同，造成结果不稳定.
|
图 7 不同深度的时间域及频率域方法的接收函数振幅谱(高斯因子2.5) 灰色表示频率域水准量(water-level)方法结果，黑色表示时间域迭代反卷积方法结果，黑框指示两者振幅谱有差异频段. Fig. 7 Amplitude spectrum of receiver function at different depth from time domain and frequency domain method (Gaussian factor 2.5) Gray lines indicate the result of the frequency domain water-level method, black lines indicate the result of the time domain iterative deconvolution method, the black box indicates the difference in the amplitude spectrum between results from two methods. |
高斯因子也是一个重要参数，下面来研究不同的高斯因子对接收函数稳定性的影响.这里分别计算高斯因子为1、1.5、2时不同深度的接收函数，并和前面高斯因子为2.5的结果进行对比，为了便于分析，也给出了各接收函数结果的振幅谱(图 8).由图可知，当高斯因子为1时(约0.5 Hz低通滤波)，对于地表和地下接收点的接收函数，两种方法的结果一致性很好，各深度接收函数结果的互相关系数在99.9%以上，且振幅谱的一致性也较好.这说明选择高斯因子为1时，各深度可以提取稳定的接收函数.当高斯因子为1.5时(约0.75 Hz低通滤波)，在0.5 km深度，两种方法接收函数结果的互相关系数略有降低，从振幅谱上可以看出，其差异是由1 Hz附近两种结果的能量不一致造成的，这恰好与0.5 km深度的地震波形垂向分量在1.01 Hz开始出现频谱极小相对应；其他深度处的两种方法接收函数结果一致性与高斯因子为1时的情况类似.当高斯因子为2时(约1 Hz低通滤波)，在0.5 km深度振幅谱中的1 Hz附近频段，两种方法接收函数能量差异变大；0.4 km深度处两振幅谱能量在1.25 Hz附近频段出现差异.当高斯因子为2.5时(约1.25 Hz低通滤波)，0.5 km和0.4 km深度处，在两种接收函数能量不同的频段内能量差异继续变大；0.3 km深度处1.65 Hz附近频段出现能量差异.由此可见，当接收点的深度确定时，出现频谱极小的频段也就确定了，选择适当的高斯因子可以压制垂向分量频谱极小对接收函数稳定性的影响，获得稳定的接收函数.
|
图 8 高斯因子及接收点深度对接收函数及振幅谱影响图 黑色表示频率域水准量反卷积方法结果，红色表示时间域迭代反卷积方法结果. (a)高斯因子为1时，接收函数(左)及其振幅谱(右); (b)高斯因子为1.5时，接收函数(左)及其振幅谱(右); (c)高斯因子为2时，接收函数(左)及其振幅谱(右); (d)高斯因子为2.5时，接收函数(左)及其振幅谱(右). Fig. 8 The influence of Gaussian factor and receiver depth on receiver function and amplitude spectrum Black lines indicate the result of the frequency domain water-level deconvolution method and red lines indicateresult of the time domain iterative deconvolution method. (a) When the Gaussian factor is 1, the receiver function(left) and its amplitude spectrum(right); (b) When the Gaussian factor is 1.5, the receiver function(left) and its amplitude spectrum(right); (c) When the Gaussian factor is 2, the receiver function(left) and its amplitude spectrum(right); (d) When the Gaussian factor is 2.5, the receiver function(left) and its amplitude spectrum (right). |
水准因子也可能影响接收函数的稳定性，以下进行探讨.通过上述研究可知，0.5 km深度处台站的接收函数稳定性受频谱极小影响更大，因此我们选定这个深度处的接收函数研究水准因子的影响.在高斯因子为2.5时，分别使用时间域迭代反卷积和频率域水准因子反卷积方法提取接收函数，并计算了振幅谱(图 9).其中频率域水准因子和反卷积中水准因子分别选择0.01和0.25.图 9a中1.1 Hz附近频段内，两种方法提取的接收函数振幅谱一致性较差，表明提取的接收函数受垂向分量的频谱极小现象的影响，在该频段内不稳定.在图 9b中可看到，0~0.8 Hz频段范围内，不同方法提取的接收函数振幅谱有较明显差别.这表明使用频率域水准因子反卷积提取接收函数时，水准因子会对接收函数造成误差，且水准因子越大误差也越大.
|
图 9 水准因子对井下地震计的接收函数的影响 灰色表示频率域水准量反卷积提取的接收函数的振幅谱，黑色表示时间域迭代反卷积提取的接收函数的振幅谱. (a)高斯因子为2.5时，0.5 km处接收函数振幅谱，其中水准量反卷积中水准因子为0.01; (b)高斯因子为2.5时，0.5 km处接收函数振幅谱，其中水准量反卷积中水准因子为0.25. Fig. 9 Effect of water-level on the receiving function from borehole seismometer Gray line and black line represent the amplitude spectrum of the receiver function from frequency domain level deconvolution method and time domain iterative deconvolution method respectively. (a) With the Gaussian factor 2.5, the amplitude spectrum of receiver function at 0.5 km, where the water-level is 0.01 for water-level deconvolution; (b) With the Gaussian factor 2.5, the amplitude spectrum of receiver function at 0.5 km, where the water-level is 0.25 for water-level deconvolution. |
反卷积方法对于频谱极小造成的接收函数不稳定的情况也有影响，这里作进一步的分析.对于一个稳定、正确的接收函数，可以将接收函数与垂向分量波形进行卷积，并将结果与径向分量进行比较，如果两者一致则说明接收函数的结果正确.对于高斯因子2.5的情况，0.5 km深度处的接收函数受频谱极小的影响相对于浅部更大，此处选取该接收函数的结果进行分析(图 10).图 10a中蓝色波形表示理论地震图，黑色波形是频率域方法提取的接收函数与理论地震图垂向分量卷积的结果，红色波形是时域方法提取的接收函数与理论地震图垂向分量卷积的结果.图 10b是图 10a中三个波形的振幅谱，不同的颜色与前面相对应.将频率域的接收函数与理论地震图垂向分量的卷积结果和理论地震图的径向分量相比较可见，在初至P波之前存在信号，这是不正确的.两波形的振幅谱在1.01 Hz附近频段明显不同，表明这是由频谱极小使得接收函数不稳定造成的.将时间域的接收函数与理论地震图垂向分量的卷积结果和理论地震图的径向分量相比较可见，虽然在波形上相似度很高，但是两波形的振幅谱在1.01 Hz附近频段仍然存在一定的不同，说明频谱极小仍然会使得接收函数在频谱极小出现的频率及其附近频段不稳定.由此可以得出结论：时间域的迭代反卷积或频率域的水准量反卷积方法难以克服频谱极小现象造成的接收函数不稳定.
|
图 10 频谱极小对时间域迭代反卷积与频率域水准量反卷积影响图 (a)接收函数与理论地震波形垂向分量卷积结果，其中蓝色是理论地震波形径向分量，黑色为频率域方法提取的接收函数与理论地震图垂向分量卷积的结果，红色是时域方法提取的接收函数与理论地震图垂向分量卷积的结果; (b)为(a)中三个波形分别对应的振幅谱，红框指示了三者幅值差异较大频段，右上为差异处放大图. Fig. 10 The influence of spectrum null on time domain iterative deconvolution and frequency domain water-level deconvolution (a) The convolution results between the receiver function and the vertical component of synthetic waveform, blue line is the radial component of synthetic waveform, black line is the result of convolution between the receiver function from the frequency domain method and the vertical component of synthetic seismogram, red line is the result of convolution between the receiver function from the time domain method and the vertical component of synthetic seismogram; (b) The amplitude spectrum corresponding to the three waveforms in (a), the red frame indicates the frequency band in which the amplitude of three waveforms before are different, the figure upper right is the enlarged image of red frame area. |
频谱极小影响接收函数的稳定性，而频谱极小的具体频率值分布与射线参数、井下地震计的埋深、震源持续时间等因素有关.
从前面频谱极小的计算公式
|
图 11 不同深度接收点频谱极小与慢度(入射角)关系图 从上到下依次为井下0.2 km、0.3 km、0.4 km、0.5 km接收点结果. Fig. 11 The relationship between spectrum null at different depth and slowness (angle of incidence) The lines from top to bottom are results at 0.2 km、0.3 km、0.4 km、0.5 km in borehole |
除了接收点深度和入射角度外，震源持续时间可能对频谱极小有影响.对于5.5级地震，震源持续时间约为3 s，6.0级地震的震源持续时间约为5 s，6.5级地震的震源持续时间约为7 s.使用与前面研究射线参数时一样的半空间模型，分别计算震源持续时间为3 s、5 s、7 s时0.5 km深度处接收函数，并将结果作傅里叶变换得到振幅谱(图 12).图中黑色的振幅谱是频率域水准量方法的结果，红色的振幅谱是时间域迭代反卷积的结果.从图中可以发现，由不同震源持续时间的理论地震图计算得到的接收函数，其振幅谱的整体特征是类似的，说明接收函数反映的是接收点下方的结构信息，而与震源特征关系不大.对于相同震源持续时间的理论地震图得到的接收函数的振幅谱，在1.5 Hz附近频段内两种方法的结果相差较大，由计算可知，由于在1.51 Hz处存在频谱极小现象，会在1.51 Hz附近的频段内导致接收函数不稳定.由此可见，频谱极小现象的出现不受震源持续时间的影响，进一步来说，一般用于提取远震接收函数的5.5级以上的地震，在井下观测时会存在频谱极小现象.
|
图 12 不同震源持续时间的接收函数的振幅谱 图中黑色的振幅谱是频率域water-level方法的结果，灰色的振幅谱是时间域迭代反卷积的结果.黑色方框中指示出两结果中振幅谱能量相差较大的频段. Fig. 12 Amplitude spectrum of receiver functions from the source with different duration time Black lines are results of amplitude spectrum from frequency domain water-level deconvolution method, gray lines are results of amplitude spectrum from time domain iterative deconvolution method. The black square indicates the frequency band in which the amplitude energy is different. |
本文利用基于平面波入射的传播矩阵方法正演计算了井下接收点的理论地震图，发现其垂向分量振幅谱存在频谱极小现象.对于相同的模型，频谱极小开始出现的频率值随深度的增加而减小，随着入射角的增大(射线参数增大)变化不大，且对深度变化的敏感性明显大于对于入射角度的敏感性.在提取接收函数时，垂向分量的频谱极小现象会使得接收函数的结果在频谱极小出现的频率及其附近的频段内不稳定，而且频率域的水准量方法或时间域的迭代反卷积方法均难以克服这种不稳定的情况.当接收点上方的P波速度结构未确定时，对于利用井下数据提取接收函数，推荐的做法是：首先分析垂向分量的振幅谱，分析频谱极小首次出现的频率，在反卷积过程中选择合适的高斯因子压制频谱极小对反卷积的影响，从而获取稳定的接收函数.
本文井下地震计的理论地震波形的算法中，除了能够计算各深度的位移外，还可计算出各层中上行P波、S波和下行P波、S波，即对波场进行一种分解.利用该波场分解的方法预期可以压制频谱极小对井下地震计接收函数的影响，需要开展研究.
致谢 中国地震局地球物理研究所“国家数字测震台网数据备份中心”为本研究提供固定台站地震波形数据，中国地震台网中心提供了地震目录数据，在此深表感谢.感谢审稿专家提供宝贵的修改意见和建议.
Aki K, Richards P G. 2002. Quantitative Seismology. 2nd ed. California: University Science Books: 74-75.
|
Chapman C H. 1976. A first-motion alternative to Geometrical Ray Theory. Geophysical Research Letters, 3(3): 153-156. DOI:10.1029/gl003i003p00153 |
Chapman C H. 1978. A new method for computing synthetic seismograms. Geophysical Journal International, 54(3): 481-518. DOI:10.1111/j.1365-246x.1978.tb05491.x |
Chen L, Zheng T Y, Xu W W. 2006. A thinned lithospheric image of the Tanlu Fault Zone, eastern China:constructed from wave equation based receiver function migration. Journal of Geophysical Research:Solid Earth, 111(B9): B09312. DOI:10.1029/2005jb003974 |
Chen L, Tao W, Zhao L, et al. 2008. Distinct lateral variation of lithospheric thickness in the Northeastern North China Craton. Earth and Planetary Science Letters, 267(1-2): 56-68. DOI:10.1016/j.epsl.2007.11.024 |
Chen X F. 1993. A systematic and efficient method of computing normal modes for multilayered half-space. Geophysical Journal International, 115(2): 391-409. DOI:10.1111/j.1365-246x.1993.tb01194.x |
Chong J J, Ni S D. 2009. Near surface velocity and QS structure of the Quaternary sediment in Bohai basin, China. Earthquake Science, 22(5): 451-458. DOI:10.1007/s11589-009-0451-1 |
Clayton R W, Wiggins R A. 1976. Source shape estimation and deconvolution of teleseismic bodywaves. Geophysical Journal International, 47(1): 151-177. DOI:10.1111/j.1365-246x.1976.tb01267.x |
Fuchs K, Müller G. 1971. Computation of synthetic seismograms with the reflectivity method and comparison with observations. Geophysical Journal International, 23(4): 417-433. DOI:10.1111/j.1365-246X.1971.tb01834.x |
Fukushima R, Nakahara H, Nishimura T. 2016. Estimating S-wave attenuation in sediments by deconvolution analysis of KiK-net Borehole seismograms. Bulletin of the Seismological Society of America, 106(2): 552-559. DOI:10.1785/0120150059 |
Gao Y, Wu J. 2008. Compressive stress field in the crust deduced from shear-wave anisotropy:an example in capital area of China. Chinese Science Bulletin, 53(18): 2840-2848. DOI:10.1007/s11434-008-0310-9 |
Gurrola H, Baker G E, Minster J B. 1995. Simultaneous time-domain deconvolution with application to the computation of receiver functions. Geophysical Journal International, 120(3): 537-543. DOI:10.1111/j.1365-246x.1995.tb01837.x |
Haskell N A. 1960. Crustal reflection of plane SH waves. Journal of Geophysical Research, 65(12): 4147-4150. DOI:10.1029/JZ065i012p04147 |
Haskell N A. 1962. Crustal reflection of plane P and SV waves. Journal of Geophysical Research, 67(12): 4751-4768. DOI:10.1029/JZ067i012p04751 |
He C S, Wang C Y, Wu Q J. 2003. The receiver function method and its new progress. Progress in Geophysics (in Chinese), 18(2): 224-228. DOI:10.3969/j.issn.1004-2903.2003.02.008 |
Helmberger D V. 1968. The crust-mantle transition in the Bering Sea. Bulletin of the Seismological Society of America, 58(1): 179-214. |
Helmberger D V, Wiggins R A. 1971. Upper mantle structure of midwestern United States. Journal of Geophysical Research, 76(14): 3229-3245. DOI:10.1029/jb076i014p03229 |
Kennett B L N, Kerry N J. 1979. Seismic waves in a stratified half space. Geophysical Journal International, 57(3): 557-583. DOI:10.1111/j.1365-246x.1979.tb06779.x |
Kennett B L N. 1980. Seismic waves in a stratified half space-Ⅱ. Theoretical seismograms. Geophysical Journal International, 61(1): 1-10. DOI:10.1111/j.1365-246x.1980.tb04299.x |
Kennett B L N, Clarke T J. 1983. Seismic waves in a stratified half-space-IV:P-SV wave decoupling and surface wave dispersion. Geophysical Journal International, 72(3): 633-645. DOI:10.1111/j.1365-246x.1983.tb02824.x |
Langston C A. 1977. The effect of planar dipping structure on source and receiver responses for constant ray parameter. Bulletin of the Seismological Society of America, 67(4): 1029-1050. |
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 |
Ligorría 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, Rainer K, Li S C. 1996. Maximal likelihood estimation and nonlinear inversion of the complex receiver function spectrum ratio. Acta Geophysica Sinica (in Chinese), 39(4): 500-511. |
Liu Q Y, Kind R. 2004. Multi-channel maximal likelihood deconvolution method for isolating 3-component teleseismic receiver function. Seismology and Geology (in Chinese), 26(3): 416-425. DOI:10.3969/j.issn.0253-4967.2004.03.006 |
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. DOI:10.3321/j.issn:0001-5733.2008.04.022 |
Okada Y, Kasahara K, Hori S, et al. 2004. Recent progress of seismic observation networks in Japan-Hi-net, F-net, K-NET and KiK-net-. Earth, Planets and Space, 56(8): xv-xxviii. DOI:10.1186/BF03353076 |
Oldenburg D W. 1981. A comprehensive solution to the linear deconvolution problem. Geophysical Journal International, 65(2): 331-357. DOI:10.1111/j.1365-246X.1981.tb02716.x |
Richards P G. 1973. Calculation of body waves, for caustics and tunnelling in core phases. Geophysical Journal International, 35(1-3): 243-264. DOI:10.1111/j.1365-246x.1973.tb02426.x |
Richards P G. 1976. On the adequacy of plane-wave reflection/transmission coefficients in the analysis of seismic body waves. Bulletin of the Seismological society of America, 66(3): 701-717. |
Shen Y, Solomon S C, Bjarnason I T, et al. 1996. Hot mantle transition zone beneath Iceland and the adjacent Mid-Atlantic Ridge inferred from P-to-S conversions at the 410-and 660-km discontinuities. Geophysical Research Letters, 23(24): 3527-3530. DOI:10.1029/96gl03371 |
Taira T, Dreger D S, Nadeau R M. 2015. Rupture process for micro-earthquakes inferred from borehole seismic recordings. International Journal of Earth Sciences, 104(6): 1499-1510. DOI:10.1007/s00531-015-1217-8 |
Tao K, Liu T Z, Ning J Y, et al. 2014. Estimating sedimentary and crustal structure using wavefield continuation:theory, techniques and applications. Geophysical Journal International, 197(1): 443-457. DOI:10.1093/gji/ggt515 |
Tsai Y B, Aki K. 1970. Precise focal depth determination from amplitude spectra of surface waves. Journal of Geophysical Research, 75(29): 5729-5744. DOI:10.1029/jb075i029p05729 |
Wang L Y, Guo Y X, Liu F, et al. 2008. Temporal vP/vS variation characteristics in different zones of China's Capital Circle area before and after Wen'an earthquake. Acta Seismologica Sinica (in Chinese), 30(3): 240-253. DOI:10.3321/j.issn:0253-3782.2008.03.003 |
Wu J, Gao Y, Chen Y T. 2008. Crustal seismic anisotropy in southeastern Capital area, China. Acta Seismologica Sinica, 21(1): 1-10. DOI:10.1007/s11589-008-0001-2 |
Wu Q J, Tian X B, Zhang N L, et al. 2003. Receiver function estimated by maximum entropy deconvolution. Acta Seismologica Sinica (in Chinese), 25(4): 382-389. DOI:10.3321/j.issn:0253-3782.2003.04.005 |
Wu Q J, Li Y H, Zhang R Q, et al. 2007. Receiver function estimated by multi-channel deconvolution. Chinese Journal of Geophysics (in Chinese), 50(3): 791-796. DOI:10.3321/j.issn:0001-5733.2007.03.018 |
Xu W W, Zheng T Y. 2002. The receiver function method and its progress. Progress in Geophysics (in Chinese), 17(4): 605-613. DOI:10.3969/j.issn.1004-2903.2002.04.007 |
Zhang G C, Wu Q J, Pan J T, et al. 2013. Study of crustal structure and Poisson ratio of NE China by H-K 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. Research of stable method to estimate receiver function in frequency-domain. Chinese Journal of Geophysics (in Chinese), 52(10): 2483-2490. DOI:10.3969/j.issn.0001-5733.2009.10.007 |
Zheng T Y, Zhao L, Chen L. 2005. A detailed receiver function image of the sedimentary structure in the Bohai Bay Basin. Physics of the Earth and Planetary Interiors, 152(3): 129-143. DOI:10.1016/j.pepi.2005.06.011 |
Zheng X F, Ouyang B, Zhang D N, et al. 2009. Technical system construction of Data Backup Centre for China Seismograph Network and the data support to researches on the Wenchuan earthquake. Chinese Journal of Geophysics, 52(5): 1412-1417. DOI:10.3969/j.issn.0001-5733.2009.05.031 |
Zhu L P, Kanamori H. 2000. Moho depth variation in southern California from teleseismic receiver functions. Journal of Geophysical Research:Solid Earth, 105(B2): 2969-2980. DOI:10.1029/1999jb900322 |
高原, 吴晶. 2008. 利用剪切波各向异性推断地壳主压应力场:以首都圈地区为例. 科学通报, 53(23): 2933-2939. DOI:10.3321/j.issn:0023-074X.2008.23.015 |
贺传松, 王椿镛, 吴庆举. 2003. 接收函数方法及其新的进展. 地球物理学进展, 18(2): 224-228. DOI:10.3969/j.issn.1004-2903.2003.02.008 |
刘启元, Rainer K, 李顺成. 1996. 接收函数复谱比的最大或然性估计及非线性反演. 地球物理学报, 39(4): 500-511. DOI:10.3321/j.issn:0001-5733.1996.04.008 |
刘启元, Kind R. 2004. 分离三分量远震接收函数的多道最大或然性反褶积方法. 地震地质, 26(3): 416-425. DOI:10.3969/j.issn.0253-4967.2004.03.006 |
罗艳, 崇加军, 倪四道, 等. 2008. 首都圈地区莫霍面起伏及沉积层厚度. 地球物理学报, 51(4): 1135-1145. DOI:10.3321/j.issn:0001-5733.2008.04.022 |
王林瑛, 郭永霞, 刘芳, 等. 2008. 文安地震前后首都圈分区波速比时变特征. 地震学报, 30(3): 240-253. DOI:10.3321/j.issn:0253-3782.2008.03.003 |
吴庆举, 田小波, 张乃铃, 等. 2003. 计算台站接收函数的最大熵谱反褶积方法. 地震学报, 25(4): 382-389. DOI:10.3321/j.issn:0253-3782.2003.04.005 |
吴庆举, 李永华, 张瑞青, 等. 2007. 用多道反褶积方法测定台站接收函数. 地球物理学报, 50(3): 791-796. DOI:10.3321/j.issn:0001-5733.2007.03.018 |
许卫卫, 郑天愉. 2002. 接收函数方法及研究进展. 地球物理学进展, 17(4): 605-613. DOI:10.3969/j.issn.1004-2903.2002.04.007 |
张广成, 吴庆举, 潘佳铁, 等. 2013. 利用H-K叠加方法和CCP叠加方法研究中国东北地区地壳结构与泊松比. 地球物理学报, 56(12): 4084-4094. DOI:10.6038/cjg20131213 |
张洪双, 田小波, 滕吉文. 2009. 稳定的频率域提取接收函数方法研究. 地球物理学报, 52(10): 2483-2490. DOI:10.3969/j.issn.0001-5733.2009.10.007 |
郑秀芬, 欧阳飚, 张东宁, 等. 2009. "国家数字测震台网数据备份中心"技术系统建设及其对汶川大地震研究的数据支撑. 地球物理学报, 52(5): 1412-1417. DOI:10.3969/j.issn.0001-5733.2009.05.031 |