2010, Vol. 53
2. 中国地震局地震预测研究所兰州科技创新基地, 兰州 730000;
3. 甘肃省地震局兰州观象台, 兰州 730000;
4. 中国地震局地球物理研究所, 北京 100081
2. Lanzhou Base of Institute of Earthquake Prediction, China Earthquake Administration, Lanzhou 730000, China;
3. Lanzhou Observation Station, China Earthquake Administration, Lanzhou 730000, China;
4. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
远震P波的波场记录中，包含了P波初动信号和P波的尾波.其中P波初动是各台站都可以记录到的而且相干的信号，而其尾波则是一种不相干的信号.邻近台站的相干信号主要和震源、介质的波速以及地球内部间断面的深度有关，因此地震波记录中相干信号部分可用于确定震源机制解和研究地球内部的速度结构，如利用P波初动反演震源机制解的方法[1]，用不同台站P波的到时资料反演地球内部速度结构的层析成像方法[2]，研究地壳速度结构和上地幔间断面的接收函数方法[3]，利用一些特殊震相对核幔边界非均匀体进行探测等[4].而地震波记录中的不相干信号，也称为地震波场扰动（seismic wavefield fluctuations：SWF），认为由震源和台站之间射线路径上的小尺度非均匀体散射引起[5, 6]，通过分析SWF信号，可以提取这些小尺度散射体的一些平均特征.
对SWF进行研究时，选择的最为理想的时间窗口是资料中不包含相干信号干扰的纯净的散射尾波.对于震中距在30°~90°之间的深源地震（深度为几百公里），P波的尾波记录中没有能量较强的深度震相（pP和sP）的干扰，相干信号的主要能量来自于P波初动，而这种相干信号可以通过一些相对简便可行的技术消除[7, 8]，以此得到较为纯净的SWF.因此深远地震的P波波场是研究SWF的理想资料.随着地震散射理论的发展，SWF被广泛地应用于地震台站下方地壳和岩石层小尺度非均匀体的研究[8~19].近来，Shen和Ritter[20]将原来用大孔径台阵资料研究台阵下方地壳和岩石层小尺度非均匀体的方法，成功应用到了《全面禁止核试验条约》组织（CTBTO）用于核监测的在兰州小孔径台阵，提取了远震P波的SWF，研究了该台阵下方地壳中小尺度非均匀体的尺度.
本研究中，我们利用Ritter等[8]提出的方法和海拉尔CTBTO小孔径地震台阵资料，对该台阵下方地壳和岩石层中的不均匀散射体尺度进行研究.
2 方法和资料 2.1 方法本研究中利用的Ritter等的方法[8]（简称RSS98）是Shaoiro & Kneib（简称SK93）[16]和Shapiro等[17]（简称SK96）方法的发展.下面简约介绍该方法的基本原理.
当远波P波穿过非均匀的介质时，记录到的P波波场中包含有相干的初动信号和不相干的SWF，即P的尾波.对于随机介质中的任一点r，在任意时刻t，整个波场可以表示为
|
(1) |
其中ut是整个波场，u是相干波场（也称主波场），uf是不相干波场，即SWF.〈〉代表了空间或者统计的平均.
根据SK93，无量纲参数ε被用来描述SWF的强弱，
|
(2) |
如果假设非弹性的能量耗散可以忽略，那么对相干波场部分的干扰主要来自于散射，则能量守恒可以利用强度表示如下：
|
(3) |
It为总强度，Ic=|〈u〉|2，If=|〈uf〉|2，将其带入上式，可得到
|
(4) |
在弱的SWF区域，〈ε2〉 < < 1（SK93），意味着波场中的能量主要是相干的部分.如果〈ε2〉 > > 1，则SWF（或不相干部分）是波场中的主要部分，意味着波在传播过程中经过了强烈的非均匀体区域.
忽略背景散射和非弹性衰减的情况下，频率f、〈ε2〉和非均匀介质的结构参数满足以下关系
|
(5) |
其中，σ为介质中波速的扰动，a是非均匀体的相关长度，L是散射层的厚度，c是散射层的平均波速.以上公式中用到了Born近似，介质被认为是Gaussian或指数性质的介质（RSS98），并且有如下两个假设. akσ2 < < 1，ka≥1，k为波数.在基于L/c2和σ2基础上，可以计算a.
2.2 台阵和资料海拉尔CTBTO台阵位于内蒙古陈巴尔虎旗的三八牧场，距离海拉尔市38km.该台阵由9个子台组成，大致以内外两个同心圆分布，直径约为3km（图 1）.各子台地震仪均放置于地下40~50m的钻井中，岩层为火山角砾岩，背景干扰较少，是比较理想的观测环境.各子台都安装有英国Guralp公司生产的CMG-3ESPV垂直向地震计，内外同心圆上子台以A和B进行区分（图 1所示），其平坦型幅频特征的频率范围为0.02~16 Hz；每分钟采样点数为40个.在B1子台上增设有一套瑞士三分向甚宽频带STS-2地震计[21].Hao和Zheng [22]利用2002~2006年的远震资料，对该台阵做了慢度-方位角校正，结果表明考虑到台阵校正后，理论慢度和方位角与观测值更为接近.该区域地壳厚度约33km[23]，可分为上地壳、中地壳和下地壳三层.经过多年的连续观测，海拉尔CTBTO地震台阵已经积累了丰富的观测资料，这些资料为台阵下方地壳和地幔内部速度结构的探测奠定了坚实的基础.
|
图 1 海拉尔CTBTO台阵子台和本研究中地震事件分布 （a）海拉尔CTBTO台阵子台分布（图中三角所示）；（b）本研究中所使用到的16个深远地震分布（黑色圆点），地震的详细资料列于表 1.三角表示海拉尔台阵中心位置. Fig. 1 Map of the stations and events (a) Distribution of the recording stations (triangles) of Hailaer CTBTO seismic array; (b) Distribution of the 16 epicenters (dots), detailed information on the events is listed in Table 1, triangle represents the center of the array. |
本研究用了9个CMG-3ESPV垂直向地震计的记录，由于仪器型号完全相同，因此我们没有去除资料中的仪器响应，也未对信号进行滤波.对于这些数据，我们截取了P之前10s和P之后60s的资料.
我们首先以如下标准对数据进行了筛选：（1）震中距范围40°~90°，以保证台阵下方的入射波波前近乎是垂直入射的；（2）震源深度大于200km，以避免pP或sP等震相干扰P的尾波；（3）至少有7道或者以上的波形可用，以便使结果有较好的统计特征.依此标准，共选取了2003~2004年的16个地震事件，其中地震目录参考USGS地震目录（http://neic.usgs.gov/neis/epic/epic_global.html），各地震的具体参数见表 1.
|
表 1 本研究中使用的地震信息 Table 1 Hypocentral parameters of the 16 teleseismic events |
我们以2003年1月17日的地震事件为例说明数据的处理过程.如图 2中的1~9道记录为该地震的原始波形，图中标示出了P波初动，其后紧跟着大于噪声水平（参考P之前的记录）的散射尾波，即SWF.我们截取了P之前10s和之后60s的波形.在该时间窗口内，相干震相只有P，而pP、sP和PcP都不在该时间窗口内.第10道波形为9道单台记录的叠加，叠加之前对数据只根据入射P波的方位角和慢度对各台资料作了P波到时校正[24].通过叠加，相干的P波初动信号被完整地保留，而SWF为随机信号，通过叠加之后会被压制，因此叠加结果中主要包含的是相干信号.每个台站原始波形中减去叠加波形（第10道波形）就可以分离出每个单台的SWF，如图 2中11~19道记录.和原始记录相比，各台站的SWF波形中，相干性和幅度最强的初动P波信号被较好压制.对于各台站的SWF波形，我们也进行了叠加，结果如图 2中的第20道记录.明显可以看到，SWF叠加波形的幅度很小，基本和P波之前的噪声水平接近，这也表明我们分离出的各台的SWF波形是不相干信号.和震源、震源区介质以及地幔深部相干的信号都包含在叠加波形中[25]，因此分离出的SWF中，和震源、震源区介质以及地幔深部相干的信号都被消除，而主要包含了不相干的散射信号.
|
图 2 表 1中1号地震的原始波形和各子台分离出的SWF.1~9曲线道为原始波形，各子台波形按照观测到的初至P波对齐.第10道曲线为1~9道记录的叠加.11~19道曲线为各子台的SWF，由1~9道波形中减去第10道叠加记录得到.第20道曲线为11~19道SWF的叠加，其幅度较弱，和初至P之前的噪声水平相当. Fig. 2 Observed waveforms and residual waveforms (fluctuation field) of the event no.01 (Table 1).The original recordings (traces 1~9) are time-shifted according to the measured P arrival time.The stack (trace 10) represents the mean wavefield.The SWF (traces 11~19) are determined by subtracting the stacked waveform of each event (trace 10) from the individual recordings.The trace 20 is the stack of the 9 SWF seismograms with the same amplitude scaling; it is close to the noise level shown prior to the P onset. |
为了进一步检验SWF是否成功地从原始记录中被分离出来，我们采用了速度谱（Vespagrams）技术[26]对信号作了处理.图 3a为以0.1s/°为步长、在0~10s/°的慢度区间内对图 2中台阵的原始地震波形做的速度谱分析结果.在速度谱分析过程中，为了突出相干性强的信号，我们采用了N次根叠加方法[27]，其中我们选用N=20.通过速度谱分析方法，可以帮助我们较为容易地识别弱震相.在速度谱分析结果中，显示出能量最大的信号（如图 3a中的白叉）对应的到时和P波初动一致，对应的慢度为7.62s/°，根据IASP91模型计算得到的理论慢度为7.54s/°.这一结果进一步说明在该时间窗口内，最主要的相干信号为P波初动信号.对于各单台的SWF（图 2中11~19道记录），我们也采用同样的参数作速度谱分析，结果如图 3b所示，在速度谱分析图中，其强度远小于原始波形速度谱分析结果，没有明显的相干信号.这说明各单台SWF波形中，主要包含不相干的散射信号，而且也没有受到其他震相干扰.
|
图 3 表 1中1号地震原始波形和各子台SWF（见图 2）的速度谱分析 （a）原始观测波形的速度谱分析结果，黑色圆圈表示根据IASP91模型计算得到的初至P波的位置，白色叉号表示由速度谱分析的得到的初至P波的位置，灰度色标表示了波形的叠加幅度；（b）各子台SWF的速度谱分析，结果显示没有明显的相干信号. Fig. 3 Vespagrams of the seismic velocity waveforms and SWF of the event No.01 (Fig. 2) (a) Vespagram of the original waveforms.The black circle represents the P arrival as predicted by IASP91, the white cross represents the maximum intensity which is the observed P arrival.Gray scaling represents the amplitude of the stacked waveforms; (b) Vespagram of the SWF, note there are no coherent arrivals. |
原始波形在频率域中的强度我们用功率谱（振幅谱的平方）来表示.如图 4中的虚线为图 2地震各原始波形功率谱的叠加平均结果Ic，实线为图 2中各道SWF功率谱的叠加平均结果If.从图中可以看到，在频率小于3Hz时，Ic>If，而频率大于3Hz后，随着频率连续增大，SWF的能量逐渐占优.表 1中其余15个地震事件的结果也和该地震事件类似.根据公式（2）我们可以得到ln（ε2）随频率变化的结果，如图 5a为图 2中地震的处理结果.结果显示，在0.5~3Hz的区间内，散射强度随着频率增加而增大，当频率更高时，〈ε2〉达到饱和.其余15个地震事件的处理结果也和该结果相似，图 5b显示了所有16个地震的结果以及平均叠加结果.由于16地震来自不同方位和不同震中距，因此ln（ε2）的平均叠加结果代表了海拉尔台阵下方的平均散射特征.图 5b也表明，在3 Hz以上，〈ε2〉达到饱和，较大能量的散射约从0.5Hz开始，因此对于所有地震事件，我们统一选定0.5~3.0 Hz的频率范围来确定γ的值.
|
图 4 表 1中1号地震的功率谱 点线代表了图 2中1~9道波形谱的叠加，即Ic，实线代表了图 2中11~19道波形谱的叠加，即If.图形显示3Hz以下Ic占优势，高于此频率，Ic和If接近. Fig. 4 Squared amplitude spectra of the velocity recordings of event No.01 The spectrum of the stack (mean wavefield, Ic) is shown as a dotted line, the stack of the squared amplitude spectrum of the residual waveforms (If) as a solid line. Note the higher intensity of the mean field below 3 Hz, which belongs to the primary P wave.At higher frequencies the Ic is close to If. |
|
图 5 ln（ε2）随频率的变化曲线 （a）表 1中1号地震对应的ln（ε2）；（b）表 1中所有地震的ln（ε2）（灰色曲线）及其叠加（黑色曲线）. Fig. 5 Plots of ln (ε2) versus frequency (a) ln (ε2) for a single event (No.01 in Table 1); (b) Stack of ln (ε2) for all analyzed events (black line) given in Table 1.Gray lines represent ln (ε2) of all single events. |
根据公式（5），我们利用最小二乘拟合来确定γ的值，以此估算海拉尔台阵下方非均匀体的性质.图 6a显示了图 2中地震事件ln（ε2 +1）的拟合结果，图 6b显示了所有地震ln（ε2+1）叠加之后的拟合结果.每一个地震事件得到γ值列于表 1中.γ的最大和最小值分别为0.0352Hz-2（表 1中2号地震）和0.0569 Hz-2（4号地震），其均值γ=0.0482±0.0069Hz-2.根据ln（ε2 +1）叠加之后的拟合结果为0.0467Hz-2.γ值的变化在一个较小的范围内，表明我们得到的结果较为稳定，同时也说明该结果和震源以及深部介质无关.
|
图 6 ln（ε2 +1）随频率的变化及其最小二乘拟合γ的结果 （a）表 1中1号地震的拟合结果；（b）表 1中所有地震ln（ε2 +1）平均值的拟合结果. Fig. 6 Plots of ln (ε2 +1) versus frequency with a quadratic LSQ fit for the determination of the parameter γ (seetext) (a) Fit for a single event (No.01 in Table 1); (b) Fit for the average of ln (ε2 +1) for all analyzed events in Table 1. |
在已知L/c2，即介质厚度L和平均速度c的情况下，根据我们计算得到的0.0352 Hz-2 < γ < 0.0482Hz-2，利用公式5，我们可以确定介质的参数σ2a.L和c的值参考前人的研究结果[23, 28]列于表 2中.我们认为，在海拉尔台阵下方散射层厚度的最大值是岩石层的厚度，An和Shi[28]利用地震面波层析成像结果和地热资料，估计出该区域的岩石层厚度约为110km.而地壳速度结构我们参考沈旭章和周蕙兰[23]利用接收函数反演得到的结果，地壳厚度约为34km，而且可以分为上中下三层.在以上结果基础上，我们可以得到不同层厚时的L/c2值，进一步得到σ2a的值.图 7显示了γ不同时，σ2a随L/c2的变化曲线.γ=0.0352 Hz-2和γ=0.0569Hz-2的两条曲线，给出了σ2a的变化范围.譬如，SWF由整个地壳中的散射体产生，则L/c2=0.9701km-1，我们可以计算出σ2a范围为0.00046~0.00074km，假如该台阵下方地壳中速度扰动σ为1%，则散射体的相关长度变化范围为4.6~7.4km.用上述方法计算得到的地壳中各层、整个地壳和整个岩石层中可能存在散射体的参数都列于表 2.本研究中得到的所有σ2a的值都小于0.01km，在这种条件下，数值实验也表明得到的结果是最为合理的[19].
|
|
表 2 海拉尔台阵下方L/c2的可能值和不同散射层速度扰动σ对应的散射体的相关长度a Table 2 Layer parameters L/c2 for possible scattering layers beneath the Hailaer CTBTO array and calculated correlation lengths a of the heterogeneities using different velocity contrasts σ for the fluctuations in a scattering layer |
|
图 7 σ2a随L/c2的变化曲线 图中显示出了对于γ的最大值、最小值和平均值对应的曲线.上、中、下地壳、整个地壳和整个岩石层根据L/c2来进行区分. Fig. 7 Plot of the statistical medium parameters σ2a versus the average quantity L/c2 of the scattering layer. Curves are given for the extremal determined values of γ and for the stacked.Upper, middle, lower, whole crust and whole lithosphere ranges are separated by their values of layer thickness L divided by the squared average seismic velocity c. |
本研究通过分析不同震中距和不同方位角的远震P波的尾波，确定了海拉尔台阵下方存在小尺度的不均匀散射体.根据16个不同的地震事件所确定的〈ε2〉（图 5）和γ的值（图 7和表 1）都在一个较小而稳定的范围内变化，说明了我们结果的可靠性.根据公式（3），我们确定了各种可能的散射层模型（表 2）.
当散射体尺寸a>L/4时，在对应的散射层中散射体数量比较少，可能产生一些相干信号，而不能用来解释具有随机性质的SWF，因此对于a>L/4的结果我们认为是不合理的，表 2中用中括号标示. Huang和Zhao[29]体波层析成像结果表明，在该台阵所在区域，波速扰动σ在地壳和岩石层中不超过3%.在这种情况下，公式（3）中σ2ak < < 1的条件可以完全满足.另外一个条件，ka≥1在a≥2.0km（对应最低频率0.5 Hz）时才能够满足，对于相关尺寸小于2.0 km的散射体，本研究中所用的理论（SK93，SSG96）不能合理地解释，在表 2中以圆括弧标示出.综上所述，本研究中观测到的海拉尔台阵远震P波尾波的SWF可以用下述散射层模型进行合理的解释：散射层的厚度为10~110km，速度的扰动在1%~3%之间，散射体的相关长度在2.0~7.4km之间.
Ritter等[8]利用该方法研究了法国Massif中部的散射体，得到γ的平均值为0.51 Hz-2.他们利用岩石层中由岩浆的浸入引起的的非均匀性来解释观测到的SWF.Rothert和Ritter[18]研究了德国东南部的小尺度非均匀体，得到γ的平均值为0.129 Hz-2，并以地壳（速度扰动为3%~7%）中的散射来对其进行解释.Hock等[19]对本研究中的方法和能量流方法[15, 30]（EFM）确定散射体相干尺度的结果进行了对比，讨论了两种方法的优缺点，并用之研究了欧洲中部和东南部11个小区域中的不均匀散射体，得到的γ值范围为0.16~0.44 Hz-2.以上结果中得到的γ值都大于我们的结果，最主要的原因可能是我们研究中所利用的台阵的孔径（约1.5km）远小于欧洲台阵的孔径（50~200km）.Shen和Ritter[20]将该方法应用到了小孔径的兰州CTBTO台阵，研究了该台阵下方的不均匀散射体，得到的γ范围为0.051~0.084 Hz-2，并以地壳（速度扰动为1%~3%）中可能存在2.4~8.2km的散射模型对其结果进行了解释，该结果和本研究较为接近.
4 结论本研究通过分离海拉尔小孔径地震台阵中16个远震P波的尾波波场，得到了SWF.对SWF分析处理后得到了较为稳定的描述不均匀散射体性质的参数（〈ε2〉和γ）.结合前人在该区域的研究工作，我们推断，在该台阵下方，分层地壳或者岩石层中的散射层模型可以较好地解释我们所观测到的SWF.该散射层模型的厚度为10~110km，速度的扰动在1%~3%之间，散射体的相关长度在2.0~7.4km之间.
致谢所用数据都由中国地震局地球物理研究所提供，两位审稿人给出了中肯和合理的评审意见，资料处理中用到了美国Sandia实验室的Matseis软件，德国Karlsruhe大学Ritter教授在方法和资料分析中给予了很多帮助，图 1绘制中利用了GMT软件[31]，在此一并表示感谢.
| [1] | Kasahars K. Computer program for a fault-plane solution. Bull. Seism. Soc. Am , 1963, 53(1): 113. |
| [2] | Aki K, Christoffersson A, Husebye E S. Determination of the three-dimensional seismic structure of the lithosphere. J. Geophys. Res , 1977, 82: 277-296. DOI:10.1029/JB082i002p00277 |
| [3] | Langston C A. Structure under Mount Rainier, Washington, inferred from teleseismic body waves. J. Geophys. Res. , 1979, 84(B4): 4749-4762. |
| [4] | Shen X Z, Zhou H L. Locating scatters in the bottom of the mantle beneath eastern Tibet with PKIKP precursors. Chinese Science Bulletin , 2009. |
| [5] | Sato H, Fehler M C. Seismic wave propagation and scattering in the heterogeneous earth. Springer-Verlag New York, 1998 |
| [6] | Mack H. Nature of short-period P-wave signal variations at LASA. J. Geophys. Res. , 1969, 74: 3161-3170. DOI:10.1029/JB074i012p03161 |
| [7] | RitterJ R R, Mai P M, Stoll G, Fuchs K. Scattering of teleseismic waves in the lower crust: observations in the Massif Central, France. Phys. Earth. Planet. Int. , 1997, 104: 127-146. DOI:10.1016/S0031-9201(97)00046-0 |
| [8] | Ritter J R R, Shapiro S A, Schechinger B. Scattering parameters of the lithosphere below the Massif Central, France, from teleseismic wavefield records. Geophys. J. Int. , 1998, 134: 187-198. DOI:10.1046/j.1365-246x.1998.00562.x |
| [9] | Aki K. Scattering of P waves under the Montana Lasa. J. Geophys. Res. , 1973, 78: 1334-1346. DOI:10.1029/JB078i008p01334 |
| [10] | Aki K, Chouet B. Origin of coda waves: source, attenuation, and scattering effects. J. Geophys. Res , 1975, 80: 3322-3342. DOI:10.1029/JB080i023p03322 |
| [11] | Flatte S M, Wu R S. Small-scale structure in the lithosphere and asthenosphere deduced from arrival time and amplitude fluctuations at NORSAR. J. Geophys. Res. , 1988, 93: 6601-6614. DOI:10.1029/JB093iB06p06601 |
| [12] | Langston C A. Scattering of teleseismic body waves under Pasadena, California. J. Geophys. Res. , 1989, 94: 1935-1951. DOI:10.1029/JB094iB02p01935 |
| [13] | Korn M. A modified energy-flux model for lithospheric scattering of teleseismic body waves. Geophys. J. Int. , 1990, 102: 165-175. DOI:10.1111/gji.1990.102.issue-1 |
| [14] | Wagner G S, Langston C A. A numerical investigation of scattering effects for teleseismic propagation in a heterogeneous layer over a homogeneous half-space. Geophys. J. Int. , 1992, 110: 486-500. DOI:10.1111/gji.1992.110.issue-3 |
| [15] | Korn M. Determination of site-dependent scattering Q from P-wave coda analysis with an energy-flux model. Geophys. J. Int. , 1993, 113: 54-72. DOI:10.1111/gji.1993.113.issue-1 |
| [16] | Shapiro S A, Kneib G. Seismic attenuation by scattering: theory and numerical results. Geophys. J. Int. , 1993, 114: 373-391. DOI:10.1111/gji.1993.114.issue-2 |
| [17] | Shapiro S A, Schwarz R, Gold N. The effect of random isotropic inhomogeneities on the phase velocity of seismic waves. Geophys. J. Int. , 1996, 127: 783-794. DOI:10.1111/j.1365-246X.1996.tb04057.x |
| [18] | Rothert E, Ritter J R R. Small-scale heterogeneities below the Graefenberg array, Germany, from seismic wavefield fluctuations of Hindu Kush events. Geophys. J. Int. , 2000, 140: 175-184. DOI:10.1046/j.1365-246x.2000.00013.x |
| [19] | Hock S, Korn M, Ritter J R R, Rothert E. Mapping random lithospheric heterogeneities in northern and central Europe. Geophys. J. Int. , 2004, 157: 251-264. DOI:10.1111/gji.2004.157.issue-1 |
| [20] | Shen X Z, Ritter J R R. Small-scale heterogeneities below the Lanzhou CTBTO seismic array, from seismic wavefield fluctuations. J. Seismo. , 2009. |
| [21] | 郝春月, 郑重, 郭燕平, 等. 中国数字地震台网(CDSN)和IMS/PS台阵的监测定位能力评估. 地震地磁观测与研究 , 2006, 27(2): 56–63. Hao C Y, Zheng Z, Guo Y P, et al. The calculation of locating ability of China Digital Seismological Network (CDSN) and IMS/PS seismic array (in Chinese). Seismological and Geomagnetic Observation and Research (in Chinese) , 2006, 27(2): 56-63. |
| [22] | Hao C Y, Zheng Z. Slowness-azimuth corrections of teleseismic events for IMS primary arrays in China. J. Seismo. , 2009. DOI:10.1007/s10950-008-9137-8 |
| [23] | 沈旭章, 周蕙兰. 接收函数近邻反演方法的改进和对海拉尔台下地壳速度结构的研究. 中国科学院研究生院学报 , 2005, 22(3): 322–328. Shen X Z, Zhou H L. Inversion of velocity structure under HLR seismic station with receiver function and NA method (in Chinese). Journal of the Graduate School of the Chinese Academy of Sciencss (in Chinese) , 2005, 22(3): 322-328. |
| [24] | Rost S, Thomas C. Array seismology: Methods and applications. Rev. Geophys. , 2002, 40: 1-27. |
| [25] | Bannister S C, Husebye E S, Ruud B O. Teleseismic P-coda analyzed by three-component and array techniques: deterministic location of topographic P-to-Rg scattering near the NORESS array. Bull. Seism. Soc. Am , 1990, 80: 1969-1986. |
| [26] | Davies D, Kelly E J, Filson J R. Vespa process for the analysis of seismic signals. Nature Physical Science , 1971, 232: 8-13. |
| [27] | Muirhead K J, Datt R. The Nth root process applied to seismic array data. Geophys. J. R. Astron. Soc. , 1976, 47: 197-210. DOI:10.1111/j.1365-246X.1976.tb01269.x |
| [28] | An M J, Shi Y L. Lithospheric Thickness of the Chinese Continent. Phys. Earth Planet. Ints , 2006, 159: 257-266. DOI:10.1016/j.pepi.2006.08.002 |
| [29] | Huang J L, Zhao D P. High-resolution mantle tomography of China and surrounding regions. J. Geophys. Res. , 2006, 111: B09305. DOI:10.1029/2005JB004066 |
| [30] | Korn M. Modelling the teleseismic P coda envelope; depth dependent scattering and deterministic structure, Special issue: Stochastic Seismic Wave Fields and Realistic Media, Neustadt, Germany, 11-15 March 1996, guest eds Korn M, Sato H, Scherbaum F. Phys. Earth. Planet. Int. , 1997, 104: 23-36. DOI:10.1016/S0031-9201(97)00044-7 |
| [31] | Wessel P, Smith W H F. New version of the generic mapping tools released. Eos. Trans. Am. Geophys. Union. , 1995, 76: 329. |