2. 中国地震局地球物理研究所,北京 100081;
3. 台湾大学地质系,台北 10617;
4. 北京大学地球与空间科学学院,北京 100871;
5. 北京科技大学材料科学与工程学院,北京 100083
2. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China;
3. Department of Geosciences, National Taiwan University, Taipei 10617, Taiwan;
4. Department of Geophysics, School of Earth and Space Sciences, Peking University, Beijing 100871, China;
5. School of Materials Science and Engineering, USTB,Beijing 100083, China
台湾是世界上构造复杂和地震活动最强烈的地区之一.欧亚板块和菲律宾海板块在这里相互碰撞.菲律宾海板块在台湾岛东部沿着琉球海沟向北下插到欧亚板块之下.欧亚板块在台湾岛以南的中国南海向东俯冲到菲律宾海板块之下.台湾岛东南部的纵谷(LV)是欧亚板块和菲律宾海板块的缝合带.该缝合带将台湾分为两个构造区.由于地质构造活动,台湾岛大部分地区在西北-东南向压缩下以大约8cm/a的速度消减[1].起始于约4 Ma的台湾造山运动在地质年代上相对比较年轻.剧烈的地壳形变导致较强的地震活动.大量地震的观测资料为研究地球结构和地震震源提供了很好的天然实验室.确定这里的三维应力场方向对于该地区的板块之间动力学研究、地震动力学研究乃至地震危险性分析均具有十分重要的科学意义.对于区分薄壳模型[2, 3]和厚壳模型[4]也具有一定的帮助.早在20 世纪80年代,李锡堤[5]测量了台湾地区断层面及该面上的滑动方向,并反演了台湾地区应力场的大致轮廓,结果表明台湾地区的压应力主要发生在台湾中南部,台中与花莲以北地区的最大主应力轴呈现一顺时针旋转的趋势,北部大屯山区显示区域性的扩张作用,宜兰地区则表现出强烈的南北向扩张作用.Yeh等[5]使用200个台湾及其附近的震源机制解求得的台湾应力状态大致符合菲律宾板块向西北挤压欧亚板块的形式,最大主应力轴为东南-西北走向,在台湾西部主压应力轴呈现扇形分布.Rau 和Wu[6]使用97个小震(2.7≤ML≤5.7)震源机制资料得到台湾全区的主压应力方向为289°,南部地区的主压应力方向为93°呈东西向.Wu等[7]根据中央气象厅台网、台湾强地面运动台网和日本气象厅数据记录的P 波初动资料和基于三维速度结构[8]得到的地震精定位结果[9],得到了1635个高质量的震源机制解,并依此求解了该地区更加详细的应力场方向,进一步验证了前人给出的应力分布结果.之后Wu等进一步对1991~2007 年的更多的地震事件进行研究,得到了4761个震源机制数据[10],虽然采用了大量地震的P 波初动资料,但仅能给出一个区域的平均应力方向,无法得到随深度分布的应力方向.
为了得到台湾地区的三维应力方向,本研究采用大量P波初动资料求解该地区的平均P轴、T轴方向,虽然P轴和T轴与地壳应力场的压缩轴和拉张轴并不完全一致,但在一定程度上预示着应力场的方向.
2 方法目前采用震源机制反演应力场有几种方法[11~13].但我们直接采用P 波极性数据反演应力场而省去了求解震源机制的中间步骤,因为:第一,采用震源机制仅利用了能够根据P 波初动求解震源机制的资料,大量较小地震的P 波初动资料无法利用;第二,由于求解单个地震的震源机制需要覆盖全方位的P波极性数据,大部分地震不能得到严格方位覆盖的P波极性数据,因此导致求解的震源机制具有较大误差[14, 15].因此我们采用大量P 波初动数据求解综合震源机制确定应力场方向的反演方法.采用大量小震P 波初动资料推测应力场方向的研究是由Aki[16]首先提出的.该方法可以利用大量的不能单独确定震源机制的P 波初动数据推断应力方向,因此得到了较为广泛的应用.李钦祖等[17]根据该方法求得了红山台和沙城台两个区域的应力场.许忠淮等[18]进一步提出了采用研究区域的多个地震和多个台站求解区域应力场的方法,并采用合成数据对此方法进行了严格测试,验证该方法求得地壳应力场的能力.许忠淮等总结了该方法得到的结果,并综合该方法得到了中国大陆的应力场方向和动力学模型[19, 20].
我们将台湾地区划分为0.25°×0.25°×2.5km的三维网格,依次求解三维网格点处的应力场方向.这里垂向距离间隔约为水平间隔的1/10,主要是考虑到地表的自由面可能导致地表应力场和深部应力场的方向有较大的差别[21].每个网格点选择周围地震在台站上观测的P 波初动符号进行反演.由于每个地震距网格点的距离不同,对网格点应力方向确定的贡献也不相同,因此我们根据距网格点的距离给予不同地震的P 波初动不同的权重.每个地震的P波初动权重w仿照Shen等[22]的大地测量数据求取应变的方式求取:
(1) |
其中,D为距离衰减常数,本研究取为25km.r为折合距离,按照下式计算:
(2) |
其中,x,y,z为地震的经度、纬度和深度,x0,y0,z0为网格点的经度、纬度和深度.这种加权重的方式意味着距网格点折合距离为20.8km 的地震的P 波初动数据只有50% 的权重,距网格点折合距离为35km的地震的P 波初动数据只有13% 的权重,而折合距离大于50km 地震的P 波初动数据的权重小于2%.为了增加计算速度,对于每个网格点只选择折合距离r≤50km 范围内的地震的P 波初动资料.另外为了保证应力场反演的质量,我们只选择P波数据个数≥100的网格点进行反演.
为了确定每个网格点与P 波初动数据拟合最好的综合震源机制解,进而分析其平均P、T轴方向,我们采用1°×1°×1°的网格尝试方法搜索P、B、T轴的方位.每一种尝试给出与综合震源机制模型不符合的P 波初动符号数(考虑权重)与总P 波初动符号个数(考虑权重)之比,我们称之为矛盾比.选择最小的矛盾比所对应的P、B、T轴方位则得到该网格点的综合震源机制解.由于该方法综合了多个地震的P波初动数据,得到的综合震源机制解的P和T方位大体可以认为是该网格点的最大和最小压应力的方位[20].如果初动数据足够多,而且有很好的方位覆盖,则只能得到一个矛盾比最小的震源机制解.在实际反演中,我们去除了少数几个有两个或两个以上最小矛盾比的综合震源机制解.
3 数据自1991年,台湾的地震观测部门在400km×550km 的范围内每年可以记录到18000 个地震事件.这些地震事件通过三维速度结构给出了精确的震源位置[9].吴逸民[10]研究小组整理了该地区1991~2007年的299210 个地震的503393 个P 波初动资料,根据这些资料和他们改进的遗传算法确定震源机制程序,获得了4761个地震的震源机制,并对该地区的应力场进行分析.研究仅用了其中的约1/5(106891个P波初动资料)数据,其他的资料由于无法确定震源机制而没有用来分析该地区的应力场.我们采用该地区的503393个P 波初动资料求解台湾地区的三维应力方向.由于台湾的区域很小,并且地壳结构复杂,难于观测到Pn波,本研究采用了Pg波的初动,离源角和方位角根据Wu等[7]的三维速度模型计算.
4 结果 4.1 总体应力场对台湾地区21.25~25.75°N,119.75~123°E,0~150km 的范围内进行空间扫描求解综合震源机制,图 2给出反演综合震源机制的P、T轴的方向.可以看出,本反演结果比起前人结果更为详尽,不仅涵盖了前人研究没有涵盖的研究区域,而且在深度方向上有很好的分辨.在浅部(蓝色棒)与Wu等[7, 10]的结果有很大的相似性,如台湾总体的P轴方向为SEE-NWW,在台湾东部P轴呈扇形向南北散开,并且呈现明显的区域应力方向的变化,表明了本文反演方法的正确性.另外,本研究结果还显示了深部P、T轴方向和倾俯角与浅部的不同,并且这些变化基本连续进行.由于菲律宾板块向欧亚板块下俯冲,导致台湾东北部的深部地震比较活跃,可以看到东北部浅部P轴方位与深部并不一致,呈现出一定的旋转,并且这种旋转随深度逐渐进行.台湾GPS 分析的数值模拟显示,自台湾中部到东北部,应力场经历了NW-SE 方向的压缩应力方向到台湾东南部张应力场的转换[23].台湾东北部的地震主要由于欧亚板块与吕宋岛弧在菲律宾海板块上的碰撞所致,其中兰屿地区发生的较深地震(约100km)是中国南海块体向菲律宾海板块之下俯冲所致.这里深部的平均P轴方位与浅部也有很大不同,也呈现随深度P轴逐渐旋转的特性.分析这些特性对于理解台湾地区的地球动力学模型具有重要意义.
很多研究者对台湾地区的3-D 速度模型进行过研究[24~28],然而这些研究采用的大部分资料为中央气象局地震台网的数据,最近Wu等[8]采用台湾强震观测项目(TSMIP)观测的41141 个S-P 波到时数据和中央气象局地震台网数据对台湾的地球速度模型进行反演,极大地改善了P 波速度结果和P/S波速度比的分辨率和覆盖程度.从其P 波速度成像结果中可以清楚地识别出板块边界,在25km 的深度上,自东向西,可以识别出欧亚板块的前陆高速区、中央山脊区的山根低速区和以纵谷为缝合带的东部菲律宾海板块的高速区.本研究采用他们得到的速度成像结果与本研究得到的P、T轴方向进行比较,研究它们之间的可能关系.
台湾的总体压应力方向表现为N110°E,致使该地区的总体构造特征大体垂直于此方向.为了进一步分析P、T轴的三维特征,我们选择剖面上网格点周围靠近该点4.1 节已经计算过的综合震源机制点,将这些点的震源机制转换为矩张量[29, 30],叠加作为剖面网格点的矩张量,求解该网格点矩张量的P、T轴方位(本征向量方向),作为该网格点的综合震源机制P、T轴方向.周围点的综合震源机制的权重按照(1)式计算.为了保证剖面点的局部应力方向特征的分析,本文只选择距剖面点的折合距离小于25km 的范围.
将图 1a 所示的AA′,BB′,CC′,DD′,EE′,FF′的P轴方向与Wu等[8]给出的Vp , Vp 扰动速度剖面叠加,得到的剖面为图 2 和图 3,将图 1b 所示的GG′,HH′,II′,JJ′,KK′,LL′,MM′,NN′的T轴方向与Wu等[8]给出的Vp , Vp 扰动速度剖面叠加,得到的剖面为图 4和图 5.台湾地区Moho面深度一直没有确切结论,为简单起见,我们将Wu等的速度剖面结果选取速度跳跃到7km/s的点,然后将这些点在空间上用FIR 滤波器[31]滤除空间周期为10km 以下的周期进行平滑,作为Moho面显示在图 2~5中.从图 2,3 可以看出,地壳中的应力场相对较为均匀,山脉之下的P轴与周围区域明显不同,特别是AA′,CC′,EE′剖面上,表现为造山运动的影响.海洋地壳中的应力场与陆地较为连续.靠近地幔顶部P轴方向从地壳中的大体水平变为倾斜,表现了地幔与地壳在力学性质上有明显差别.在AA′剖面海洋地壳的浅部,P轴方向接近垂直于地面,表现为正断拉张性质.AA′剖面的70~80km的深处,P轴方向比较凌乱,80~100km 的P轴方向又呈现水平,由此往下P轴变陡.这可能是中国南海板块下插导致的各深度层次应力不同所致.在FF′剖面,海洋地壳中的应力场较为凌乱,深处较为均匀.P轴分布与Vp 扰动(图 3)的相关性不太明显.图 4所表现的T轴方向分布在地壳浅部较为均匀,在壳幔分界处有一定的转换,但没有P轴明显.图 5的大部分剖面15~45km 的深度范围内存在T轴变陡的区域,对于剖面GG′,II′,JJ′表现为Vp 低速扰动,对于KK′,LL′,MM′,NN′表现为Vp 高速扰动区域,且自西向东逐渐加深.对应于Vp 高速扰动区域为本尼奥夫带.此处的T轴变得垂直表明了本尼奥夫带具有更明显的逆冲性质.
根据台湾地区1991~2007 年的P 波初动资料采用综合震源机制解的方法求取了台湾地区地下的三维应力场方向.台湾地区的总体应力场在东部地区呈现扇形分布,更为详细地展现了该地区的总体应力场特征.台湾地区地下深处的应力场与浅部有较为显著的差别,在Moho面附近有较为显著的应力场方向转向变化,说明台湾地区的地壳和地幔的力学性质有明显差别.地壳中的应力场相对比较一致.台湾大部分深部15~45km 范围内的T轴较陡,在东部剖面上与本尼奥夫带有较好的对应,为菲律宾板块下插到欧亚板块之下的证据.在西部的剖面上与Vp 的低速扰动对应.这些结果为台湾地区的动力学模拟及解释提供了新的约束.
本研究采用矛盾比最小的准则搜索综合震源机制的P、T轴方向,数据越多对模型的约束越好,台站的方位覆盖越好对综合震源机制解的约束越好.本文采用最小矛盾比是否唯一来判断解的好坏,如何客观估计解的误差是采用P 波初动符号求解震源机制要解决的问题,这需要在今后的工作中加强该方面的研究.但这里的综合震源机制表现为较好的连续性,说明结果也具有一定的可信度.
严格来讲,本文所用的P 波综合震源机制除了与地球内部应力场的方向有关外,还与应力场的相对大小和摩擦系数有关[31],但通过改变摩擦系数和改变应力相对大小在一定程度下均可得到相同的综合P波辐射花样,表明反演摩擦系数和应力相对大小的困难[32].因此作为一级近似,本研究仅对应力场的方向进行确定.
本文研究只是初步结果,还可能可以从中解译很多与构造有关的信息.下一步的工作将结合更为细致的构造信息和地球物理模拟结果进行综合解译.
致谢该研究得益于与台湾师范大学叶恩肇博士和中央研究院地球所许雅茹博士的讨论,中国地震局地球物理研究所许忠淮研究员提出了宝贵意见,论文评审者提出了修改性建议,谨向他们表示衷心感谢.本文绘图采用GMT 绘图软件[33].
[1] | Yu S B, Chen H Y, Kuo L C. Velocity field of GPS stations in the Taiwan area. Tectonophysics , 1997, 274(1-3): 41-59. DOI:10.1016/S0040-1951(96)00297-1 |
[2] | Suppe J. Mechanics of mountain building and metamorphism in Taiwan. Mem. Geol. Soc. China , 1981, 4: 67-89. |
[3] | Davis D, Suppe J, Dahlen F A. Mechanics of fold-and-thrust belts and accretionary wedges. J. Geophys. Res. , 1983, 88(B2): 1153-1172. DOI:10.1029/JB088iB02p01153 |
[4] | |
[5] | 李锡堤. 大地应力分析与弧陆碰撞对于台湾北部古应力场变迁之影响. 台北: 国立台湾大学地质研究所 , 1986: PP202. Li X D. Crustal Stress Analysis and the Effects of Arc-Continent Collision on North Part of Taiwan Region (in Chinese). Taibei: Institute of Geosciences, National Taiwan University (in Chinese) , 1986: PP202. |
[6] | Yeh Y H, Barrier E, Lin C H, et al. Stress tensor analysis in the Taiwan area from focal mechanisms of earthquakes. Tectonophysics , 1991, 200(1-3): 267-280. DOI:10.1016/0040-1951(91)90019-O |
[7] | Rau R J, Wu F T. Active tectonics of Taiwan orogeny from focal mechanisms of small-to-moderate-sized earthquakes. TAO , 1998, 9(4): 755-778. |
[8] | Wu Y M, Zhao L, Chang C H, et al. Focal mechanism determination in Taiwan by genetic algorithm. Bull. Seism. Soc. Amer. , 2008, 98(2): 651-661. DOI:10.1785/0120070115 |
[9] | Wu Y M, Chang C H, Zhao L, et al. Seismic tomography of Taiwan: improved constraints from a dense network of strong-motion stations. J. Geophys. Res. , 2007, 112: B08312. DOI:10.1029/2007JB004983 |
[10] | Wu Y M, Chang C H, Zhao L, et al. A comprehensive relocation of earthquakes in Taiwan from 1991 to 2005. Bull. Seism. Soc. Amer. , 2008, 98(3): 1471-1481. DOI:10.1785/0120070166 |
[11] | Wu Y M, Hsu Y J, Chang C H, et al. Temporal and spatial variation of stress field in Taiwan from 1991 to 2007: Insights from comprehensive first motion focal mechanism catalog. Earth Planet Sci. Lett. , 2010, 298(3-4): 306-316. DOI:10.1016/j.epsl.2010.07.047 |
[12] | Gephart J W, Forsyth D W. An improved method for determining the regional stress tensor using earthquake focal mechanism data: application to the San Fernando earthquake sequence. J. Geophys. Res. , 1984, 89(B11): 9305-9320. DOI:10.1029/JB089iB11p09305 |
[13] | Michael A J. Determination of stress from slip data: faults and folds. J. Geophys. Res. , 1984, 89(B13): 11517-11526. DOI:10.1029/JB089iB13p11517 |
[14] | Horiuchi S, Russo G, Hasegawa A. Discrimination of fault planes from auxiliary planes based on simultaneous determination of stress tensor and a large number of fault plane solutions. J. Geophys. Res. , 1995, 100(B5): 8327-8338. DOI:10.1029/94JB03284 |
[15] | 许忠淮, 汪素云, 高阿甲, 等. 我国部分早期震源机制解答的重新测定. 地震地磁观测与研究 , 1994, 16(4): 34–42. Xu Z H, Wang S Y, Gao A J, et al. Redetermination of some early focal mechanism solutions of Chinese earthquakes. Seismological and Geomagnetic Observation and Research (in Chinese) (in Chinese) , 1994, 16(4): 34-42. |
[16] | 万永革, 吴忠良, 周公威, 等. 中国地震震源机制测定结果的比较. 地震地磁观测与研究 , 2001, 22(5): 1–15. Wan Y G, Wu Z L, Zhou G W, et al. China seismic mechanism comparison for different determination. Seismological and Geomagnetic Observation and Research (in Chinese) (in Chinese) , 2001, 22(5): 1-15. |
[17] | Aki K. Earthquake generating stress in Japan for the years 1961 to 1963 obtained by smoothing the first motion radiation patterns. Bull. Earthq. Res. Inst. , 1966, 44(2): 447-471. |
[18] | 李钦祖, 王泽皋, 贾云年, 等. 由单台小地震资料所得两个区域的应力场. 地球物理学报 , 1973, 16(1): 49–61. Li Q Z, Wang Z G, Jia Y N, et al. Stress field obtained for two regions from weak earthquake data recorded at a single seismic station. Chinese J. Geophys. (Acta Geophysica Sinica) (in Chinese) (in Chinese) , 1973, 16(1): 49-61. |
[19] | 许忠淮, 阎明, 赵仲和. 由多个小地震推断的华北地区构造应力场的方向. 地震学报 , 1983, 5(3): 268–279. Xu Z H, Yan M, Zhao Z H. Evaluation of the direction of tectonic stress in North China from recorded data of a large number of small earthquakes. Acta Seismologica Sinica (in Chinese) (in Chinese) , 1983, 5(3): 268-279. |
[20] | 许忠淮, 汪素云, 黄雨蕊, 等. 由大量地震的资料推断的我国大陆构造应力场. 地球物理学报 , 1989, 32(6): 636–647. Xu Z H, Wang S Y, Huang Y R, et al. The tectonic stress field of Chinese continent deduced from a great number of earthquakes. Chinese J. Geophys. (Acta Geophysica Sinica) (in Chinese) (in Chinese) , 1989, 32(6): 636-647. |
[21] | Xu Z H, Wang S Y, Huang Y R, et al. Tectonic stress field of China inferred from a large number of small earthquakes. J. Geophys. Res. , 1992, 97(B8): 11867-11877. DOI:10.1029/91JB00355 |
[22] | Bokelmann G H R, Beroza G C. Depth-dependent earthquake focal mechanism orientation: Evidence for a weak zone in the lower crust. J. Geophys. Res. , 2000, 105(B9): 21683-21695. DOI:10.1029/2000JB900205 |
[23] | Shen Z K, Jackson D D, Ge B X. Crustal deformation across and beyond the Los Angeles basin from geodetic measurements. J. Geophys. Res. , 1996, 101(B12): 27957-27980. DOI:10.1029/96JB02544 |
[24] | Rau R J, Ching K E, Hu J C, et al. Crustal deformation and block kinematics in transition from collision to subduction: global positioning system measurements in northern Taiwan, 1995-2005. J. Geophys. Res. , 2008, 113: B09404,. DOI:10.1029/2007JB005414 |
[25] | Rau R J, Wu F T. Tomographic imaging of lithospheric structures under Taiwan. Earth Planet. Sci. Lett. , 1995, 133(3-4): 517-532. DOI:10.1016/0012-821X(95)00076-O |
[26] | Ma K F, Wang J H, Zhao D P. Three-dimensional seismic velocity structure of the crust and uppermost mantle beneath Taiwan. J. Phys. Earth , 1996, 44(2): 85-105. DOI:10.4294/jpe1952.44.85 |
[27] | Shin T C, Chen Y L. Study on the earthquake location of 3-D velocity structure in the Taiwan area. Meteorol. Bull. , 1988, 42: 135-169. |
[28] | Kim K H, Chiu J M, Pujol J, et al. Three-dimensional Vp and Vs structural model associated with the active subduction and collision tectonics in the Taiwan region. Geophys. J. Int. , 2005, 162(1): 204-220. DOI:10.1111/gji.2005.162.issue-1 |
[29] | Aki K, Richards P G. Quantitative Seismology: Theory and Methods. San Francisco: W H Freeman, Vol 1 , 1980: 1-932. |
[30] | 万永革. "地震静态应力触发"问题的研究. 北京: 中国地震局地球物理研究所, 2001. Wan Y G. Research on "Seismic Static Stress Triggering" Problem (in Chinese) . Beijing: Institute of Geophysics, China Seismological Bureau, 2001 |
[31] | 万永革. 数字信号处理的MATLAB实现. 北京: 科学出版社, 2007 : PP326 . Wan Y G. Digital Signal Processing Based on MATLAB (in Chinese). Beijing: Science Press, 2007 : PP326 . |
[32] | 万永革, 盛书中, 许雅儒, 等. 不同应力状态和摩擦系数对综合P波辐射花样影响的模拟研究. 地球物理学报 , 2011, 54(4): 994–1001. Wan Y G, Sheng S Z, Hsu Y J, et al. Effect of stress ratio and friction coefficient on composite P wave radiation patterns. Chinese Journal of Geophysics (in Chinese) (in Chinese) , 2011, 54(4): 994-1001. |
[33] | Wessel P, Smith W M F. New, improved version of generic mapping tools released. EOS Trans. , 1998, 79(4): 579. |