2011, Vol. 54
2. 中国地震局地震预测研究所兰州科技创新基地,兰州 730000;
3. 东京大学地震研究所,东京 113-0032
2. Lanzhou Base of Institute of Earthquake Science, China Earthquake Administration, Lanzhou 73000, China;
3. Earthquake Research Institute, University of Tokyo, 1-1-1 Yayoi, Bunkyo-ku, Tokyo 113-0032, Japan
在板块构造理论中，最基本的模型是刚性岩石圈“漂浮"在黏性软流圈之上.随着地震波观测记录的积累，人们得到地球内部信息也越来越多，由此所得速度结构的分辨率也越来越高，对岩石圈和软流圈动力学模式的理解也逐渐由简单变得精细.地质学上的很多证据和地震波层析成像结果都表明，岩石圈结构存在较强的横向不均匀性，且在沟弧盆地作用区域可能发育大量低速异常体.如在我国东海沟弧盆结构体系，层析成像结果[1， 2]就表明软流圈内存在大量低速异常体，其深度约180km.部分强低速异常体在地震学性质上甚至会表现出一种低速界面，其深度的确定，对于了解岩石圈形成､演化和动力学模式具有重要意义.
远震P波接收函数方法是研究地壳和上地幔间断面最为有效的地震学工具之一[3， 4].该方法通过远震资料反褶积分离出间断面的相关震相，进而研究地壳和上地幔间断面，图 1a显示了接收函数中常用震相的射线路径.在进行上地幔410km 和660km深度间断面(后文以410 和660 表示)研究中，接收函数的一次透射转换波发挥了其他地震学手段难以企及的作用，并取得了丰硕的成果[5~9].但是对于100~200km 深度间断面(假设为低速界面，后文简称做L 界面)，其对应的一次透射转换震相(Pls)到时往往和来自地壳多次波的到时比较接近，导致Pls震相辨认较为困难.图 1b 显示了IASP91 地球模型[10]中Moho面和深度为165km 的低速界面对应的一次透射转换波及多次波走时曲线.可以看到L界面的一次透射转换震相Pls到时和Moho界面的多次波非常接近.
|
图 1 接收函数中常用震相的射线路径及走时曲线 (a) -次透射转换震相和多次波的射线路径；实线代表P波，虚线代表S波；(b) IASP91模型中Moho面和深度为165 km的L界面对应的一次透射转换波及多次波走时曲线.同时也标示出了 P410s和P660s的走时曲线. Fig. 1 The ray path of general phases used in the receiver functions and their travel times (a) Sketch of ray paths of first order converted phase and multiples.Solid lines represent the P wave, and dotted lines represent the S wave；(b) The travel time curves of first order converted phase and multiples from Moho and L discontinuity in 165 km depth of IASP91 model. |
为了避免地壳多次波对Pls的干扰，近年来地震学家们发展了S波接收函数方法[11~13].和P波接收函数相比，S 波接收函数上Slp 震相在S 震相到达之前，而多次波在S震相之后，因此在S波接收函数上较容易识别该震相，该方法目前也被应用到了许多地区来探测岩石圈底界面深度(在全球大部分地区，该深度为100~200km)[14~17].但是在计算S波接收函数时，对波形资料信噪比要求较高，能选择出来用于S波接收函数研究的资料不到观测资料的一半.这一要求也限制了该方法的使用.相反，P 波接收函数计算方法简单，结果稳定.如果在P 波接收函数上能明确地确定L 界面多次波的震相，则可有效研究L 界面深度.Tian等[18， 19]利用接收函数的多次波对青藏高原下方的地壳速度结构进行了深入的研究，而利用地壳以下深度间断面的多次波进行地球内部结构探测的工作目前还较少.
全球数字化地震观测起始于20 世纪80 年代末，截止目前，一些历史悠久的固定观测台站已经记录到了大量远震波形资料，使研究者从这些海量资料的接收函数(后文中接收函数如不做特殊说明都指P波接收函数)中提取可靠L 界面的多次波确定台站下方L 界面的深度成为可能.位于上海的SSE台站在我国的地震台站中历史较为久远，是我国首批加入美国地震学研究联合会(IRIS)的台站之一，观测资料较为丰富.本文将提取该台站接收函数上L 界面多次波，研究该台站下方L 界面的深度.
2 观测接收函数上海SSE 台站已有100多年历史，该台站利用佘山台山洞(台基为基岩)放置地震仪，观测环境的背景干扰小，记录资料质量高.自1996年进行数字化观测以来，积累了大量观测资料.文献[8]利用该台站2003~2005年112条观测接收函数，得到台站下方410和660深度分别为410km 和665km.Ma和Zhou[20]也利用该台站同样的资料，结合面波频散，得到台站下方地壳厚度和波速比分别为30.2km和1.77.由于受当时资料数量限制，未对其他间断面进行深入讨论.
Crotwell和Owens[21]基于IRIS网站的数据管理系统开发了接收函数自动计算系统(EARS)，并将全球众多台站远震接收函数经挑选后放在该网站(http://ears.seis.sc.edu/Data/Summary/gauss_2.5/IC/SSE/station.html).本文从EARS 下载了SSE 台站自1996~2009 年以来30°~90°震中距的478条接收函数.接收函数计算中，截取了P 波之前20s和之后120s的资料，使用时间域多次迭代反褶积技术[22]，迭代次数为400 次，选择高斯低通滤波因子为2.5，对应的拐角频率约0.5Hz.图 2显示了台站和地震事件分布.
|
图 2 SSE台站和本研究中所用到的地震事件分布三角表示台站位置，黑点表示地震分布，等值线代表震中距. Fig. 2 Map of the station and events Triangle is the location of the station SSE，the points represent the event locations, and the contour lines represent the epicenter distance. |
震中距接近的接收函数上相同震相到时较为接近，因此将所有接收函数按照1°震中距的间隔进行归并，在每个震中距间隔内将接收函数直接叠加，以压制无关噪声干扰，突出来自间断面的相关震相.图 3a为SSE 台站观测接收函数在不同震中距的归并叠加结果.在观测接收函数上可以较为清晰地看到PmS､P410s 和P660s 震相以及地壳多次波，且P410s和P660s的到时和IASP91 模型预测的较为一致(图 3a中绿线所示).此外，在50s到P660s之间还明显存在振幅为负的信号.
|
图 3 观测接收函数时间､深度序列 (a)不同震中距接收函数归并叠加结果；同时也标示了根据 IASP91模型计算得到的P410s和P660s的理论走时；(b)依据 Pds震相计算得到的不同震中距归并叠加接收函数的深度序列 (左边部分)及叠加(右边部分)，灰色阴影部分表示根据 Bootstrap方法得到的误差，黑线指示了 L界面多次波的变化趋势(由理论模拟得到，详见文字说明)；(c)依据PdpPds震相计算得到的不同震中距归并叠加接收函数的深度序列(左边部分) 及叠加(右边部分)，170km低速界面的震相也被标示. Fig. 3 The time and depth serieses of observed receiver functions (a) Waveforms of the binned receiver functions in the different distances.The travel time curves of P410s and P660s phases rom IASP91 model also marked； (b) Depth series from binned RFs referring the Pds phase of different epicentral distances (lett) and the stacking result (right).The gray shadow indicates the twice of the standard deviation estimated by the bootstrap method.The black line displays the tendency of multiple of L discontinuity copied from synthetic RFs (see in Fig.3b).The black line shows the tendency of the multiples of L discontinuity-PlpPls (see details rn text) ； (c) Same as (b) but referring the PdpPds phase.The phase of relative to the L discontinuity is also marked. |
观测接收函数为时间序列，通过不同震中距间隔内的叠加，虽然可以看到某些相关震相的走时趋势，但人们一般更为关注间断面深度.为此，在射线追踪基础上，参考IASP91模型，依据Pds(d为不同深度间断面，在本研究中d=1，2，3，…，800km)震相到时将不同震中距归并叠加的接收函数时间序列转换为深度序列[8， 9].图 3b为不同震中距归并叠加接收函数深度序列，通过该过程，和间断面相关的一次透射转换震相被对齐，如图中的Pms､P410s和P660s震相.利用Bootstrap 重采样方法[23]对不同震中距深度序列进行叠加，并以标准差的2 倍作为误差来确定间断面深度，如图 3b 右侧.据此得到410和660 的深度分别为410±0.8km 和659±1.0km.由于不同震中距地壳多次波到时差别不大，在图 3b 中100~200km 内对地壳多次波和一次透射转换波不能很好地区分，因此无法据此判断该深度间断面.另外，在该图上还看到500~660km深度之间存在明显负信号(图中为红色)，该信号基本上沿一倾斜直线分布，具有较强相关性，因此是一个震相.该震相不和P660s平行，表明其不是一次透射转换震相，可能是来自浅部某个间断面的多次波.
为突出多次波，依据多次波PdpPds震相将接收函数转换成深度序列，结果如图 3c所示.该结果中类似P410s和P660s的一次透射转换震相沿着倾斜方向，而不同深度多次波PdpPds会对齐，图 3b中倾斜的红色负震相在该图 170km 深度上对齐，且叠加结果中也清晰地看到了该震相.同样利用Bootstrap重采样方法得到该间断面深度为170±1.5km.由于该震相幅度为负，意味着对应间断面上下两侧波速由大变小.依此推断，该台站下方约170km 深度存在一个低速界面.
3 全波场理论接收函数模拟为了进一步解释本文结果，利用反射率法计算了全波场理论地震图[24].所构建地球模型的地壳为单层，各参数来自接收函数和面波联合确定的结果[20]，在170km 深度加入一个速度跃变为4% 的低速界面，其他部分直接采用IASP91模型.如图 4a为所构建的1000km 深度以内地球模型.震源为简单爆炸源，震源深度30km, 这和本研究实际资料的平均震源深度较为接近.如图 4b和4c为全波场理论地震图的垂向和径向分量.在理论地震图中也标示了参考IASP91模型计算得到的常见震相走时曲线.
|
图 4 本研究中构建的SSE台站下方地球模型及全波场理论地震图 (a) SSE台站下方1000km深度以内地球模型及与IASP91的比较；(b)径向理论地震图，所有震相按照P对齐，同时标出了一些常见震相的走时曲线；(c)垂直向理论地震图. Fig. 4 The earth model beneath SSE station and full wavefield synthetic seismograms. (a) Local earth model beneath SSE station and IASP91 model；(b) Radial synthetic seismograms aligned by P phase, some usual phases are also marked； (c) Vertical synthetic seismograms aligned by P phase. |
采用和计算观测接收函数完全相同的参数，根据理论地震图计算理论接收函数.图 5a为不同震中距理论接收函数，其中Pms､PmpPms､PmsSms､P410s和P660s到时都和观测接收函数较为一致.L 界面对应的多次波PlpPls也明显呈现于其中.采取和处理实际资料相同的方法与参数，将理论接收函数依据一次透射转换震相Pds和多次震相PdpPds分别转换到深度序列，如图 5b和5c所示.图 5b 中410和660的深度都被很好地还原到了模型预设值.图 5c确定的170km 深度低速界面也和预设模型一致.170km 深度L 界面多次震相PlpPls在图 5b中沿着一个倾斜的方向随震中距分布(图 5b中黑色直线指示).将指示PlpPls倾斜方向的黑线以完全相同的比例复制到观测接收函数图像上(图 3b).比较图 3b和图 5b可以看到观测接收函数上PlpPls分布趋势和理论模拟结果非常一致.此外，理论接收函数上P660s之后，L 界面的另一个多次震相PlpSls(PlsPls)的分布趋势依然较为清晰，采取同样方式，将指示该震相的黑色直线也复制到观测接收函数上(图 3b).沿该方向，观测接收函数上也明显存在一致震相.所构建L 界面模型，能较合理解释观测接收函数中大部分震相，特别是可以非常自然地解释L 界面对应的两个到时不同的多次波，因此该模型是合理的.
|
图 5 理论接收函数时间､深度序列 (a)不同震中距理论接收函数，P410s和P660s的理论走时曲线也显示于图中；(b)依据Pds震相计算得到的不同震中距理论接收函数的深度序列，同图 3b，黑线指示了 L界面多次波的分布趋势；(c)依据PdpPds震相计算得到的不同震中距理论接收函数的深度序列，同图 3c. Fig. 5 The time and depth serieses of synthetic receiver functions (a) Synthetic receiver functions of different epicentral distances； (b) Same as Fig .2b for synthetic receiver functions.The black line displays the tendency of multiple of L discontinuity；(c) Same as Fig.3b for synthetic receiver functions. |
在利用接收函数确定一个间断面时，最为理想的情况是同时观测到该间断面对应的三个震相，即图 1a中所示的Pds､PdpPds和PdpSds(PdsPds).Zhu和Kanamori[25]提出的接收函数H-k叠加方法就是通过提取观测接收函数这三个震相来确定地壳厚度和波速比.Moho面对应的这三个震相幅度较大，易于识别，因此在确定Moho深度时该方法是高效稳定的.但对较深处低速界面而言，同时确定这三个震相较为困难.如图 1b所示，Pls的到时和Moho的多次波时间较为接近，一般情况下较难识别.在50s到P660s之间，除PlpPls震相外，无其他明显相关震相.PlpSls(PlsPls)一般在P660s之后，理论模拟表明这个时间段接收函数波形较为复杂(图 5a)，会对该震相的识别造成干扰，且只在55°~80°震中距内分布较为清晰.在观测不到一次透射转换波情况下，这两个多次波提供了L 界面非常重要的信息，据此可以确定界面深度.
上述确定间断面深度时，参考了IASP91模型，该模型和SSE 台站下方的速度结构存在一定差异，会导致所确定上地幔间断面和L界面深度存在一定误差.为减小误差，本文参考全球S波层析成像结果[26](该模型结果于2006年进行了升级)来构建该台站下方上地幔模型，地壳部分依旧采用图 4a的单层模型.图 6a为构建的该台站地球模型.在该模型基础上，参考Pds和PdpPds震相将观测接收函数转换为深度序列，如图 6b和6c.根据该结果，确定出410和660 的深度分别为403±0.7km 和655±0.9km, L 界面深度为170±0.6km.由于本文利用了尽可能多的观测资料，因此所确定410和660的深度也应更加可靠.
|
图 6 SSE台站下方改进的地球模型及参照此得到的观测接收函数深度序列 (a)参考S波层析成像结果构建的SSE台站下方800 km之内地球模型；(b)参考S波层析成像模型，依据P办震相得到的观测接收函数深度序列，同图 3b；(c)参考S波层析成像模型，依据PdpP办震相得到的观测接收函数深度序列，同图 3c. Fig. 6 Improved Earth model beneath SSE station and depth serieses of receiver functions (a) Local Earth model beneath SSE station from S tomography (Grand, 2002，updated in 2006) ； (b) Same as Fig.2b but referring the local earth model； (c) Same as Fig.3b but referring the localearth model. |
本文发现的L界面和岩石圈底界面具有相同性质，都是地震波速度降低的间断面，面波层析成像结果表明[1]，该区域岩石圈厚度约70km, 和本研究得到L 界面的深度相差较大，因此本研究得到的L界面可能是非岩石圈底界面.层析成像结果还表明，该区域下方软流圈中，在180km 深度附近低速异常体发育，因此本文所确定低速界面应是由低速异常体所导致.
Debayle和Sambridge[27]利用全球面波频散资料，在Voronoi图基础上根据大圆路径分布自动划分网格，得到了面波全球层析成像结果(http://rses.anu.edu.au/cadi/code/tomoeye).在SSE 台站所在区域，其观测数据有较好的射线覆盖.图 7为沿SSE 台站所处经度的一个剖面(纬度从25°~35°).该结果显示在160~180km 深度之间S 波速度扰动发生了较大变化，也说明该深度区域具备形成低速界面的条件，与本文所得170km 结果相一致.基于此进一步推断，该区域软流圈中低速异常体的发育程度可能比以前人们所想的更强烈，以至于形成明显的低速界面.
|
图 7 面波层析成像模型在SSE台站下方的纵剖面 (Debayle and Sambridge, 2004)该剖面经度和台站经度相等，纬度从25°到 35°.颜色表示了速度的扰动(dR/Vs%). Fig. 7 Surface-wave tomography slice through the SSE station (Debayle and Sambridge, 2004)The longitude of the slice equals to the longitude of SSE，the start and end latitude is 25°N and 35°N respectively.The color represents the perturbation of Vs(dVs/Vs%). |
从SSE台站多年远震观测接收函数中提取了深度170km 低速界面的多次波.全波场理论地震图得到的理论接收函数，合理自然地解释了该观测结果.结合前人层析成像结果，推测该低速界面可能由软流圈中大量发育的低速异常体形成.能有效稳定地提取低速界面的多次波，大量高质量的观测接收函数是必要保证.随着数字化地震波观测资料的不断积累，该方法可以越来越广泛被应用于更多的台站或者台阵.且全球岩石圈的平均深度在100~200km, 因此本方法可望有效的用于岩石圈深度的确定.
致谢中国地震局和国家留学基金委共同资助了作者在东京大学ERI一年的学术访问.东京大学ERI的Hitoshi Kawakatsu 教授和Teh-Ru AlexSong博士在资料的分析中给予了较多意见，德国地学研究中心(GFZ)的马延路博士在理论地震图计算程序调试中给予了大量帮助，和中国地质科学院地质力学所安美建博士的讨论对于本研究结果的解释帮助较多，两位审稿人给出了中肯的意见.图 2绘制中利用了GMT 软件[28]，在此一并表示感谢.
| [1] | 蔡学林, 朱介寿, 曹家敏, 等. 华南地区岩石圈三维结构类型与演化动力学. 大地构造与成矿学 , 2003, 27(4): 299–310. Cai X L, Zhu J S, Cao J M, et al. Three-Dimensional tectonic types and evolutional dynamics of lithosphere of south China region. Geo. Tectonica of Memalloqenia (in Chinese) , 2003, 27(4): 299-310. |
| [2] | 蔡学林, 朱介寿, 曹家敏, 等. 中国东南大陆边缘带岩石圈三维结构—动力学模式. 海洋地质与第四纪地质 , 2005, 25(3): 25–34. Cai X L, Zhu J S, Cao J M, et al. Three-Dimensional and dynamic types of the continental margin lithosphere in southeast China. Marine Geology & Quaternary Geology (in Chinese) , 2005, 25(3): 25-34. |
| [3] | Langston C A. Structure under Mount Rainier, Washington, inferred from teleseismic body waves. J.Geophys. Res. , 1979, 84(B4): 4749-4762. |
| [4] | Ammon C J. The isolation of receiver effects from teleseismic P waveforms. Bull. Seismol. Soc. Am. , 1991, 81: 2504-2510. |
| [5] | Yuan X, Ni J, Kind R, et al. Lithospheric and upper mantle structure of southern Tibet from a seismological passive source experiment. J.Geophys.Res. , 1997, 102(B12): 27491-27500. DOI:10.1029/97JB02379 |
| [6] | Li X, Kind R, Yuan X. Seismic study of upper mantle and transition zone beneath hotspots. Phys. Earth Planet. Inter. , 2003, 136: 79-92. DOI:10.1016/S0031-9201(03)00021-9 |
| [7] | Ai Y, Zheng T, Xu W, et al. A complex 660kmn discontinuity beneath northest China. Earth Planet. Sci. Let. , 2003, 212: 63-71. DOI:10.1016/S0012-821X(03)00266-8 |
| [8] | Shen X, Zhou H, Kawakatsu H. Mapping the upper mantle discontinuities beneath China with teleseismic receiver functions. Earth Planets Space , 2008, 60(7): 713-719. DOI:10.1186/BF03352819 |
| [9] | 沈旭章, 周蕙兰. 接收函数CCP-PWS偏移方法探测东北地区620km深处低速层. 科学通报 , 2009, 54(17): 3067–3075. Shen X, Zhou H. The low-velocity layer at the depth of 620 km beneath Northeast China. Chinese Sci. Bull. (in Chinese) , 2009, 54(17): 3067-3075. DOI:10.1007/s11434-008-0559-z |
| [10] | Kennett B L N, Engdahl R E. Travel times for global earthquake location and phase identification. Geophys. J. Int. , 1991, 105: 429-465. DOI:10.1111/gji.1991.105.issue-2 |
| [11] | Farra V, Vinnik L P. Upper mantle stratification by P and S receiver functions. Geophys. J. Int. , 2000, 141: 699-712. DOI:10.1046/j.1365-246x.2000.00118.x |
| [12] | Oreshin S, Vinnik L, Peregoudov D, et al. Lithosphere and asthenosphere of the Tien Shan imaged by S receiver functions. Geophys. Res. Lett. , 2001, 29: 1029-1031. |
| [13] | Yuan X, Kind R, Li X, et al. The S receiver functions: synthetics and data example. Geophys. J. Int. , 2006, 165(2): 555-564. DOI:10.1111/j.1365-246X.2006.02885.x. |
| [14] | Li X, Kind R, Yuan X, et al. Rejuvenation of the lithosphere by the Hawaiian plume. Nature , 2004, 427: 827-829. DOI:10.1038/nature02349 |
| [15] | Kumar P, Yuan X, Kumar M R, et al. The rapid drift of the Indian tectonic plate. Nature , 2007, 449: 894-897. DOI:10.1038/nature06214 |
| [16] | Chen L, Cheng C, Wei Z G. Seismic evidence for significant lateral variations in lithospheric thickness beneath the central and western North China Craton. Earth Planet Sci. Lett. , 2009, 286: 171-183. DOI:10.1016/j.epsl.2009.06.022 |
| [17] | Kawakatsu H, Kumar P, Takei Y, et al. Seismic Evidence for Sharp Lithosphere-Asthenosphere Boundaries of Oceanic Plates. Science , 2009, 324: 499-502. DOI:10.1126/science.1169499 |
| [18] | Tian X, Wu Q, Zhang Z, et al. Joint imaging by teleseismic converted and multiple waves and its applicaton in the INDEPTH-III passive seismic arry. Geophys. Res. Lett. , 2005, 32: L21315. DOI:10.1029/2005GL023686 |
| [19] | Tian X, Wu Q, Zhang Z, et al. Identification of multiple reflected phases from migration receiver function profile: An example for the INDEPTH-III passive teleseismic P waveform data. Geophys. Res. Lett. , 2005, 32: L08301. DOI:10.1029/2004GL021885 |
| [20] | Ma Y, Zhou H. Crustal thicknesses and Poisson's ratios in China by joint analysis of teleseismic receiver functions and Rayleigh wave dispersion. Geophys. Res. Lett. , 2007, 34: 12304. DOI:10.1029/2007GL029848 |
| [21] | Crotwell H, Owens T. Automated Receiver Function Processing. The Electronic Seismologist column of Seismological Research Letters, 2005 |
| [22] | Ligorria J, Ammon C. Iterative deconvolution and receiver-function estimation. Bull. Seismol. So. Am. , 1999, 89: 1395-1400. |
| [23] | Efron B, Tibshirani R J. An introduction to the Bootstrap. London, Chapman & Hall, 1998 |
| [24] | Wang R. A simple orthonormalization method for the stable and efficient computation of Green's functions. Bull. Seism. Soc. Am. , 1999, 89: 733-741. |
| [25] | Zhu L, Kanamori H. Moho depth variation in southern California from teleseismic receiver functions. J. Geophys. Res. , 2000, 105(B2): 2969-2980. DOI:10.1029/1999JB900322 |
| [26] | Grand S P. Mantle shear-wave tomography and the fate of subducted slabs. Phil. Trans. R. So. Lond A , 2000, 360: 2475-2491. |
| [27] | Debayle E, Sambridge M. Inversion of massive surface wave data sets: model construction and resolution assessments. J. Geophys. Res. , 2004, 109: B02316. DOI:10.1029/2003JB002652 |
| [28] | Wessel P, Smith W. New version of the generic mapping tools released. Eos Trans Am Geophys Union , 1995, 76: 329. |