2016, Vol. 59
2. 防灾科技学院, 河北三河 065201;
3. 中国地震局地球物理研究所, 北京 100081
2. Institute of Disaster Prevention, Sanhe Hebei 065201, China;
3. Institute of Geophysics, China Earthquake Administration, Beijing 100081, China
There are two steps in this research. Firstly, by using the global seismicity EHB catalog (Mb≥4.7) from 1960 to 2008, we investigate the Tonga-Kermadec subduction zone (18.5°S-28.5°S) geometry by plane fitting. The plane fitting result indicates that the strike of slab is nearly 196°and the dip of slab is nearly 48° in the range of 18.5°S-28.5°S. In the following stress field study, the Tong-kermadec subduction zone is divided into 9 sub-regions to solve stress field based on the strike and dip we obtained. Using focal mechanisms from the Global Centroid Moment Tensor (1993-2012), the present-day stress tensor for different sub-region and depth interval is solved based on grid search method.
In shallow part (70~300 km), stress inversion results indicate that the stress regime has strong heterogeneity with the variety of plunges, and the azimuth of maximum principal compression stress is rotated near 24°S, which is associated with the subducting Louisville Seamount Chain. The ridge push within the slab is subjected to resistance due to the Louisville Seamount Chain, which leads to the stress rotation. As the motion of the Chain is southward along trench, the "stress shadow zone" is generated in the southwest of Louisville Seamount Chain.
In intermediate part (300~500 km), the azimuths of maximum and minimum principal stress are rotated from north to south gradually. This phenomenon is probably associated with compression due to north-south mantle flow for the slab. The compression for the slab is greater near the northern edge of the slab beneath the Lau Basin. But the compression induced from north-south mantle flow is weakening gradually in southward direction from the northern edge.
In deep part (500~700 km), the distribution of principal stress is down dip compression. It may be caused by the resistance of its downward motion due to increases of depth and mantle viscosity.
Our results also show different stress plunge between main subduction zone and the western offset slab. The horizontal or vertical maximum compression stress axis in the offset slab shows no down dip compression induced from the upper subduction zone, which confirms the offset slab separation from main subduction zone.
全球俯冲带动力学作为板块构造学、地球动力学、地震学等学科重要组成部分，一直是国内外学者研究的热点领域(臧绍先等，1992；Hayes et al.，2012).Wadati(1935)和Benioff(1954)先后发现了俯冲板片的深震地震带，后人称之为Wadati-Benioff带，该地震带的发现为俯冲带动力学研究提供了天然试验场.近年来，随着全球地震台网的不断增加，观测技术不断发展，全球地震数据的观测数量和质量均有很大改善，为全球俯冲板片的构造形态及动力学研究提供了很好的数据支持(Chen et al.，2004; Syracuse and Abers，2006；Fukao and Obayashi，2013；张浪平等，2013，2014).
随着应力反演技术的不断发展，有关利用地震观测资料反演俯冲带应力场的研究，同样成为了俯冲带动力学研究的热点问题之一.Christova研究小组曾运用Gephart和Forsyth(1984)的网格搜索法反演得到了琉球—九州、伊豆—小笠原、千岛、瓦鲁阿图等俯冲带不同深度的应力场(Christova and Scholz，2003; Christova et al.，2004; Christova，2004，2005，2015)；Wada等(2010)采用Arnold和Townend(2007)的贝叶斯参数估计方法和Michael(1987)的应力张量线性反演方法得到了卡斯卡迪亚俯冲带的应力状态；Meighan等(2013)采用Hardebeck和Michael(2006)提出的时空应力线性反演方法(SATSI)证实了波多黎各海沟板片拆离现象.
汤加—克马德克俯冲带作为澳大利亚—太平洋板块边界的一部分，是全球板块大地构造学研究的热点区域.该俯冲带是全球地震最活跃的俯冲带之一，从浅部到700 km深度左右，表现出强烈的地震活动性，并且由北向南，地震分布不均，有大量的地震空区出现(图 1).前人大量的地震层析成像(Van der Hilst，1995; Hall and Spakman，2002; 鲁人齐等，2013)、各向异性(Chen and Brudzinski，2001，2003; Foley and Long，2011)、地震活动性(Hamburger and Isacks，1987; Okal and Kirby，1998)等研究结果，已经大致描绘出该俯冲带复杂的构造形态及动力学特征.Bonnardot等(2009)通过应力张量反演大致给出了汤加—克马德克俯冲带中深部应力场特征，并根据应力方向不连续特征证实了汤加-克马德克地区21°S以北的板片拆离现象.矿物相变、脱水反应、多重俯冲事件和板片拆离都是造成汤加俯冲带复杂构造形态以及大量深震发生的可能原因，许多动力学过程导致了汤加—克马德克俯冲系统现今复杂的构造格局(Ruellan et al.，2003).
 | 图 1
汤加—克马德克地区构造图与地震分布
图中震源机制球为下半球投影；红色实线框为研究区范围；黑色虚线代表剖面(1)—(9)的中心位置(每个剖面宽度为100 km)；R1—R9为每个剖面所对应的区域；红色箭头代表太平洋板块向澳大利亚板块的俯冲方向. 该地区构造特征据文献(Hall and Spakman，2002).Fig. 1 Regional tectonic setting and seismicity in the Tonga-Kermadec area Focal mechanisms are marked by lower hemisphere projections. The study region is marked by red solid line. The black dashed lines represent central locations of cross-sections of(1)—(9). Every cross-section is 100 km wide. R1—R9 segments are correspond to cross-section(1)—(9). Red arrow represents subducting direction from Pacific Plate to Australian Plate. The regional tectonic map is revised from Hall and Spakman(2002). |
由于该地区构造复杂，虽然已有少量有关该俯冲板片应力场研究结果(Bonnardot et al.，2009)，但应力场反演网格较粗，分辨率较低，研究程度不够深入.本文拟收集汤加—克马德克地区1960—2008年定位精度较高(Mb≥4.7)的全球EHB目录(Engdahl et al.，1998)，拟合出研究区内俯冲带的走向、倾角等定量参数.然后选取1993—2012年Global CMT目录(MW≥4.8)(Dziewonski et al.，1981；Ekström et al.，2012)，将汤加—克马德克俯冲带划分成9个区域(图 1)，求解各个区域不同深度的区域应力场，探讨深部地幔流对主应力方位的影响、复杂构造单元与应力场非均匀性的关系、深部“偏移”板片与主俯冲带的应力状态差异等动力学问题.
2 研究步骤与方法2.1 俯冲带几何定量参数求取首先对汤加—克马德克俯冲带进行平面拟合，旨在得到研究区内俯冲带几何形态的定量参数，为下一步应力反演剖面选取提供依据.本文参考张浪平等(2013，2014)有关全球地震目录选取的经验(地震定位的精度从高到低依次为EHB目录、NEIC目录、CMT目录)，选取定位精度较高的全球EHB地震目录，根据地震的空间展布，将汤加—克马德克俯冲带拟合出一个近似的线性平面.本文假定俯冲平面为到每个地震震源位置距离平方和最小的平面，由此建立求解俯冲平面的数学模型，结合使用模拟退火全局搜索和高斯牛顿局部搜索算法(万永革等，2008)，给出全局搜索最优的俯冲带平面参数(走向、倾角).该方法已经在小震估计大震断层面参数中得到广泛应用，例如，万永革等(2008)首先应用该方法得到了唐山地震发震断层面参数，吕坚等(2013)和盛书中等(2014)分别给出了芦山地震和玉树地震发震断层的定量参数.本文将该方法推广到俯冲板片形态拟合，定量给出俯冲带的走向和倾角.
2.2 应力场反演方法本文利用汤加—克马德克俯冲带研究区内Global CMT(4.8≤MW≤7.8)目录进行区域应力反演，拟得到不同位置与深度的应力场.其中应力张量反演方法采用网格搜索法(Wan et al.，2015; 万永革，2015)，该方法采用滑动方向与剪应力方向最为一致的准则，以更加精细的搜索网格(1°×1°×1°×0.1)反演出三个主应力方向以及R值(Gephart and Forsyth，1984)：
 | 图 2 不同R值所对应的主应力空间分布图 红色代表最大主压应力轴；蓝色代表最大主张应力轴；R=0时最大主张应力确定，最大主压应力不确定； R=1时最大主压应力确定，最大主张应力不确定；R=0.5时最大主压应力、最大主张应力均确定．Fig. 2 Principle stress spatial distribution correspond to different R-value Maximum principle compressive stress axis is shown by red pattern, maximum principle extension stress axis is shown by blue pattern. When R-value equal 0, S3 is certain but S1 is uncertain. When R-value equal 1, S1 is certain but S3 is uncertain. When R-value equal 0.5, both S1 and S3 are certain. |
汤加—克马德克俯冲带范围较大(图 1)，19°S以北板块俯冲形态变化多样(Bonnardot et al.，2009)，28°S以南深部地震活动性较弱，形态复杂.故本文选取18.5°S—28.5°S作为研究区域，该区域构造形态简单，俯冲倾角沿海沟走向变化较缓，有利于进行平面拟合(Giardini and Woodhouse，1984; Jiao et al.，2000; Chen et al.，2004).
由于该俯冲带在深部500～700 km西侧有“板外(outboard)”地震的存在(图 3a灰色椭圆标记区)，Okal和Kirby(1998)将该震群定义为“偏移地震(Offset Event)”，它们不属于主体俯冲带.另外，0～60 km深度区间可能会受到上覆板片以及俯冲带自身弯曲的影响，约400 km深度处，主体俯冲带倾角也发生局部变化(图 3a)，以上这些因素会使拟合平面偏离主体俯冲带，故我们采用Chen等(2004)用于线性拟合所选取该地区地震的深度范围(70～420 km)进行解算.
 | 图 3 拟合俯冲带位置(a)及俯冲带平面两侧的地震频度分布(b) (a) 中黑色小圆点代表EHB目录地震分布，灰色实线为拟合出的俯冲带平面，灰色虚线为俯冲带平面到0 km和700 km的延伸； 灰色虚线椭圆标记为偏移震群； (b) 中DF为参与拟合的地震到俯冲带平面的距离． Fig. 3 The location of fitted subduction zone (a) and histogram of earthquakes sorted by their distances to the subduction plane (b) (a) Earthquakes from EHB catalogs are shown by small black dots. The grey solid line is fitting plane of subduction zone. The grey dashed line represents extending to 0 km and 700 km depth. The offset events are marked by gray dashed ellipse; (b) DF represent distance from earthquakes to subduction plane. |
选取1960—2008年定位精度较高的该地区EHB目录，拟合深度范围内(70～420 km)共1043个地震.通过两个步骤的拟合过程，最终得到研究区俯冲板片的总体几何形态(图 3a中灰色实线)走向为196.0°(标准差为0.3°)，倾角为48.2°(标准差为0.4°)，并且参与拟合的地震位置到俯冲带平面的距离服从正态分布(图 3b).与Chen等(2004)结果相对比(表 1)，该结果很好地表现了研究区俯冲带的整体几何形态.
 | 表 1 汤加—克马德克俯冲带平面拟合结果对比 Table 1 The comparison of Tonga-Kermadec subduction zone of the best linear fitting results |
通过上一步平面拟合得到的俯冲带走向与倾角，我们沿垂直于汤加-克马德克海沟的方向(即俯冲平面倾向)，在研究区内(图 1红色区域)构建9个宽度为100 km的剖面，分别对应R1—R9区域，然后分别反演9个区域内浅部(60～300 km)、中部(300～500 km)、深部(500～700 km)的区域应力场.Bonnardot等(2009)反演该俯冲带应力场时，在28°S—19°S的纬度范围内，将60～300 km 深度区间划分成三个区域，300～700 km深度区间划分成两个区域进行应力场解算，其区域划分太大，应力反演结果平滑程度过高，忽略了局部信息，不能反映细节应力结构.本文沿俯冲带(18.5°S—28.5°S)走向每隔100 km选取一个区域进行浅、中、深部应力张量反演.
本文选取研究区域内Global GMT(4.8≤MW≤7.8)目录，对每个区域(R1—R9)的数据沿N16°E方向进行等面积投影(图 4).根据每个区域深度剖面震源机制解的分布特征，选取主俯冲带地震，确定每个区域内反演所用震源机制解.部分区域(R4—R7)深部(500～700 km)西侧有明显“偏移震群”的存在，为了不影响主俯冲带应力场的反演结果，暂且将其去除(图 4中R4—R7区域灰色虚线椭圆标记部分).
 | 图 4 分区震源机制解剖面图及应力反演区域划分 图中震源机制球为沿俯冲带走向(N16°E)方向里半球投影，灰色空心圆代表EHB目录，灰色实线框包围了每个区域反演所用震源机制解，灰色虚线椭圆区含义见图 3． Fig. 4 The cross-sections of focal mechanisms and stress inversion sub-regions The focal mechanisms are marked by hemisphere projections in N16°E direction. The grey open dots st and for earthquakes from EHB catalogs. The sub-region focal mechanisms used in stress inversion are marked by grey solid line. The meaning of grey dashed ellipse region is the same with Fig. 3. |
沿不同剖面、不同深度划分出22个主俯冲带应力反演区域(去除“偏移”板片的影响).其中每个区域震源机制不少于6个.我们对每个区域内反演所用的震源机制进行P、T轴统计分析(图 5a)，并对汤加—克马德克主俯冲带浅、中、深部每个区域内应力场进行求解(表 2).然后，从每个区域震源机制P、T轴分布(图 5a)、主应力方位分布(图 5c，5d)、主应力倾伏角变化(图 5b)等方面，结合每个剖面地震活动性以及几何形态，对俯冲带不同深度的应力场进行详细描述，为下一步动力学讨论做铺垫.
 | 图 5 汤加—克马德克主俯冲带应力场分布 (a) 区域反演所用震源机制P、T轴统计(下半球投影)； (b) 反演得到的不同深度每个区域最大剪应力所在节面与主应力轴的等面积投影，其中绿色弧线为80%置信度下反演出的最大剪应力所在节面集合, S1和S3周围的封闭曲线为得到的主压应力轴和主张应力轴的80%置信区间； (c)，(d) 分别为不同深度范围S1轴和S3轴最优解水平投影，箭头越长，S1轴或S3轴越水平. 黑 色粗虚线分别代表60 km，300 km，500 km，700 km等深线.黑色细虚线含义同图1. Fig. 5 The stress field of Tonga-Kermadec main subduction zone (a)The lower hemisphere projections illustrate the P and T axes distribution in every stress inversion sub-region;(b)The equal area projections of maximum shear stress planes and principle stress axes of every sub-region by stress inversion in different depth region，the green arcs represent the nodal planes with maximum shear stress in confidence level of 80%. Closed curves surrounding S1 and S3 indicate the confidence interval of the compressive and extensional stress axis respectively，in confidence level of 80%;(c) and (d)represent horizontal projections of S1 and S3 axes with best solutions respectively，the longer arrow is，the more horizontal of the S1 or S3 axis is. Black thick dash lines correspond to depth contours of 60 km，300 km，500 km and 700 km，respectively. The meaning of black thin dash lines is the same with Fig. 1. |
| | 表 2 汤加—克马德克主俯冲带最优应力参数表 Table 2 The best fitting parameters of the Tonga-Kermadec slab stress field |
在60～300 km浅部范围内，俯冲带形态变化较稳定，每个剖面俯冲带倾角较一致，地震活动强烈(图 4)，但也有少量地震空区出现(图 4的R6).由研究区域震源机制解的P、T轴分布(图 5a)可见，P轴分布较集中，主要分布于NNW-EW范围内，少量分布于SW向，且倾伏角偏小；T轴大部分分布于SEE-SE向，少量分布于SW向，倾伏角变化较分散.我们发现R1区P、T轴分布均没有明显的优势方位.Bonnardot等(2009)研究发现在14°S—19°S之间60～300 km深度范围内存在双层地震带，并且该双地震层上层应力结构以沿俯冲带俯冲方向挤压(Down-dip Compression)为主，下层应力结构以沿俯冲方向拉伸(Down-dip Extension)为主，这可能是该区域P、T轴无明显优势方向分布的原因.本文的R1区浅部位于该双地震层的最南端(Bonnardot et al.，2009)，应力结构复杂，故并没有对R1浅部区进行应力张量求解.
反演得到的最优应力方位见图 5c、5d的红色箭头.可以看到R2—R4区域的最大主压应力方向较为一致地趋近EW方向，但R5区的最大主压应力轴发生较大偏转，转向NNE-SSW向，与俯冲带倾向存在较大夹角，这是研究区域浅部应力状态变化的一个转折点.在R6—R9区域又逐渐逆时针偏转到NW-SE的方向.最大主张应力轴在该深度范围内也表现出方向的多变性，R2区张应力轴为N30°W，但R3—R6区张应力轴方位发生逆时针旋转，其方位与俯冲带倾向相同；南部R7、R8、R9区张应力轴方位继续逆时针向南旋转，变为SWW-NEE向.
从反演得到的主应力倾伏角分布特征来看(图 5b)，浅部区域主应力轴倾伏角变化多样： R3、R5、R9区最大压应力轴倾伏角较大，接近垂直；R7、R8区最大压应力轴与最大张应力轴倾伏角均较小，接近水平.还有大量区域(R2、R4、R6)最大压应力轴与最大张应力轴倾伏角均介于30°～50°之间，与本文拟合得到的俯冲带倾角48°接近.
5.2 中部(300～500 km)结果在300～500 km范围内，从EHB目录和Global CMT目录所构建的地震剖面图像中(图 4)可以看出，在该中部深度区内俯冲带几何形态多变，南部区域(R5—R9)地震活动性相对较弱，出现“地震空区”段.故根据震源机制解在该深度范围内的分布特征(图 4)，本研究只选取了4个区域(R1—R4)进行应力张量求解.震源机制解的P、T轴分布特征(图 5a)表明，P轴分布较为集中，R1和R2区域以近NS向为主，R3区域的P轴倾伏角较大，R4区域的P轴接近EW向，倾伏角较小；相对于P轴分布，300～500 km深度范围的T轴分布相对比较分散，大部分倾伏角较小.
图 5c、d展示了该深度范围内反演的最优主应力方位.我们发现R1和R2区域内，最大主压应力轴方位为N5°～25°W，但在与其邻近的R3和R4区域内，主压应力轴方位为N67°—70°W，发生明显偏转，这也与图 5a所示的P轴优势方位相吻合.该深度区内最大主张应力轴方位从R1到R4区也由近EW向逐渐偏转至NNE-SSW向.
从主应力倾伏角分布来看(图 5b)，R1—R3区由北向南最大主压应力轴倾伏角逐渐变大，最大主张应力轴倾伏角逐渐变小；R4区最大主压、张应力轴均接近水平.
5.3 深部(500～700 km)结果在500～700 km范围内，汤加—克马德克俯冲带地震活动相对强烈(图 4)，俯冲带形态也较为复杂.在主俯冲带西侧，有“偏移”震群的影响(图 3a、图 4灰色虚线椭圆标记区).为了不影响主俯冲带应力场，暂且将“偏移”板片的震源机制解去除，只反演主俯冲带的应力状态.观察该深度范围主俯冲带震源机制解的P、T轴分布特征(图 5a)，在R1—R4、R7、R8区P轴分布较集中于NNW-NWW向，倾伏角相对较小；而R5、R6、R9区P轴分布较离散；T轴分布整体较离散，主要分布于SEE-SSE区间内，倾伏角多以近水平为主.
图 5c、5d中的深蓝色箭头给出了该深度范围内的最优主应力方位.可以看到R1—R4区最大主压应力轴方位一致性较好，大体上均沿N42°—54°W向分布，与俯冲带倾向近似一致；而R5—R8区主压应力轴发生较小偏转，大致为N48°—87°W.该深度区内最大主张应力轴方位分布非均匀性更为突出：R2—R4以及R7、R8区主张应力轴方位与俯冲带倾向近似一致，而R1、R5、R6、R9区主张应力轴相对于俯冲带倾向逆时针偏转，为近NE-SW向.
从主应力倾伏角分布来看(图 5b)，对于500～700 km深度范围，最大主压应力轴倾伏角变化较为稳定，介于39°～68°之间，与俯冲带平面拟合结果较为一致.而最大主张应力轴倾伏角变化多样.
6 动力学意义的探讨6.1 应力结构的非均匀性所反映的不同应力环境反演结果显示(图 5b)汤加—克马德克主俯冲带主应力倾伏角多变，不同深度范围局部表现出不同的应力性质，我们认为这可能与该俯冲带受到多重力源作用有关.根据主应力剖面分布图所展示的应力轴倾伏角可以看出(图 6)，浅部(60～300 km)区主压应力轴绝大多数分布趋势为沿俯冲带俯冲方向分布，倾伏角与俯冲带倾角相近，这可能与受到俯冲带上部向下的推挤作用有关(Bonnardot et al.，2009).在该浅部区内，我们明显发现R7、R8应力轴倾伏角均近水平，呈现出走滑型应力结构，说明这两个区域整体处于剪切应力环境中.从R值的分布来看，R值普遍大于0.5(图 6)，说明反演得到的最大主压应力状态较确定，但R2、R4区域R值较小，主张应力状态确定.汤加—克马德克俯冲带浅部应力结构表现出强烈的非均匀特征.
 | 图 6
汤加—克马德克主俯冲带主应力分布剖面图
S1和S3轴最优解为沿俯冲带走向(N16°E)里半球投影，其中蓝色半球投影代表R值为0.1～0.3；红色半球投影代表R值为0.4～0.6；绿色半球投影代表R值为0.7～0.9.主应力空间图样为沿主俯冲带走向视图. 红色和蓝色含义同图 2. 灰色实线含义同图 3. Fig. 6 The cross-sections of principle stress of Tonga-Kermadec main subduction zone The best solutions of S1 and S3 axes are illustrated by hemisphere projections in strike(N16°E)direction. Blue，red and green hemisphere projections represent 0.1～0.3 interval of R-value，0.4～0.6 interval of R-value and 0.7～0.9 interval of R-value，respectively. The pattern of principle stress spatial distribution is rotated to strike of main subduction zone. The meaning of red and blue is the same with Fig. 2. The meaning of gray solid lines is the same with Fig. 3. |
在中部(300~500 km)区应力轴倾伏角出现渐变(图 6)：R1到R3区由北向南最大主压应力轴倾伏角逐渐变大，最大主张应力轴倾伏角逐渐变小，并且R3区表现出明显的正断型应力结构.Richards等(2011)通过构建瓦鲁阿图俯冲拆离板片与俯冲的太平洋板块演化模型，认为该俯冲带北部400～600 km左右深度具有强烈的地震活动性的原因是由于拆离板片与俯冲的太平洋板块发生碰撞的结果，并导致俯冲带的变形.按照上述观点，R1—R2区应处于强烈的挤压环境中而表现出最大主压应力轴水平分布，并且方位应近乎垂直于俯冲带走向.但本文结果显示，中部R1—R2区最大主压应力轴并不是水平分布，方位也未垂直于俯冲带走向(图 6)，说明该拆离板片目前并未对主俯冲带产生影响，这与Chen和Brudzinski(2001)得到的结论相符.R4区显示出主应力轴水平的走滑型应力结构，并且R值为0.5，说明该区域反演得到的三个主应力轴状态确定.这种走滑型应力结构的产生原因可能是400 km左右深度水平向地幔流的动力作用导致局部的剪切应力环境(Giardini and Woodhouse，1984，1986).
在深部500～700 km范围内俯冲带应力结构一致性较好，最大主压应力轴倾伏角介于39°～68°之间，与俯冲带平面倾角拟合结果较为一致(图 6).说明深部主要处于沿俯冲带俯冲方向挤压(Down-dip Compression)的应力环境中.Warren等(2007)通过研究汤加俯冲带深部地震发现存在很多近水平或近垂直的断层面分布，这与我们在该深度范围反演得到的应力节面(图 5b中最大剪应力所在节面)相符合.这种断层面的产生正是由于沿俯冲带俯冲方向挤压的应力分布而导致，并且从反演出的R值中(图 6)，可以看出在深部R值普遍偏大，说明该应力场最大主压应力轴反演结果较确定.产生这种应力状态的动力学原因可以解释为：在地幔转换带随着深度的增加，地幔黏性增加，俯冲的太平洋板块向下运动受到阻力，很难进入到下地幔中，从而导致最大压应力轴沿俯冲方向分布的结果(Isacks and Molnar，1971; Billen et al.，2003; Bonnardot et al.，2009).这种应力结构导致了600 km附近地震活动性剧增，俯冲板片厚度增加(图 4)，地震层析成像的结果证实了这一观点(Van der Hilst，1995; Van der Hilst et al.，1997).
6.2 主应力方位的局部旋转所揭示的动力学特征从本反演得到的主应力方位来看(图 5c、5d)，最大主压应力轴方位在深部分布较均匀，基本上与俯冲带俯冲方向一致，但浅部与中部发生了两次明显偏转.浅部区域R2—R4区最大主压应力轴整体方位与R5—R9区有明显不同，R4到R5区发生近70°偏转，并且R3—R9区由北向南应力形因子均大于0.5(图 6)，压应力轴的方位是可信的.我们推断发生这种偏转现象可能是受到俯冲的太平洋板块洋底构造单元—路易斯维尔海链(LSC)的影响(图 7).路易斯维尔海链是一条抗震带(Giardini and Woodhouse，1986; Ruellan et al.，2003; Bonnardot et al.，2009)，它于5～6百万年前控制了上覆板块劳盆地—阿弗尔槽弧后扩张系统，并与汤加—克马德克海沟斜交(图 7).随着汤加—克马德克俯冲体系不断演化，太平洋板块向西俯冲，海链与汤加—克马德克海沟的交汇点沿海沟走向向西南方不断迁移.由于海链的俯冲，在交汇点处上覆板块弧后区域与海链产生强烈的挤压，导致交汇点处上覆板块(澳大利亚板块)南部与北部不同的区域构造特征(Balance et al.，1989; Ruellan et al.，2003; Bonnardot et al.，2007).我们认为当太平洋板块向下俯冲时，在洋底的路易斯维尔海链对洋脊推力起到一定的阻挡作用.由于海链沿海沟向南迁移，在海链西南侧产生“应力影区”，R5区及其南部附近区域(“应力影区”)受到向下的推挤作用减弱，导致最大主压应力方位发生旋转，并不主要沿俯冲方向分布(图 7).并且根据GPS得到的板块汇聚速率(DeMets et al.，1990)，可以看出在交汇点及南部明显比北部俯冲速率小(图 7)，北部俯冲相对南部要快，这也与北部相对于南部受到向下的推挤作用较大的推论相符.由此处向南的地震活动减弱也支持海链西南侧产生“应力影区”的推断.
 | 图 7
汤加—克马德克俯冲体系示意图，据Smith等(2001)修改
图中太平洋—澳大利亚板块汇聚方向与速率据DeMets等(1990)，推测的路易斯维尔海链的俯冲深度据Bonnardot等(2009). Fig. 7 The sketch of Tonga-Kermadec subduction system(Modified from Smith et al.，2001) The direction and rate of Pacific-Australia plate convergence are from DeMets et al.(1990). The predicted depth of LSC subduction is from Bonnardot et al.(2009). |
在中部(300～500 km)R1、R2区最大压应力轴方位为N5°—25°W，而在R3、R4区逆时针旋转到N67°—70°W，我们推断这可能是向南流动的地幔流对俯冲岩石圈板块产生动力作用，并逐渐减弱的体现(图 7蓝色箭头).Giardini和Woodhouse(1986)通过对该俯冲带400 km左右深度的震源机制研究认为俯冲板块遭受了一个南北向水平地幔流作用.随后根据剪切波分裂所揭示的各向异性结果，Smith等(2001)证实了劳盆地下方由北向南地幔流的存在，这种由北向南流动的地幔流是由于俯冲的太平洋板片北部边界(劳盆地北端)的撕裂作用而进入劳盆地下方的.Conder和Wiens(2007)利用地幔角落流和沿海沟走向的水平地幔流动力学模拟出劳盆地和岛弧下方的各向异性方向为近南北向或沿海沟走向，这也与观测得到的结果相一致(Smith et al.，2001).还有一些研究人员根据萨摩亚地幔柱产生的氦同位素进入劳盆地下方地幔的地球化学特征(图 7)，充分证实了由于板片的撕裂所导致向西流动的地幔物质向南偏转，沿南北向进入劳盆地下方地幔的现象(Turner and Hawkesworth，1998).本文300～500 km深度R1、R2区距离劳盆地北部边界(俯冲板块边界处)较近(图 1，图 7)，由北向南流动的地幔流对俯冲板块的推挤作用在劳盆地北部边界作用较大，可能是导致这两个区域内最大主压应力方位呈近NS向分布的原因.随着向南逐渐远离俯冲板块边界(R3、R4区)，这种地幔流挤压作用也随之减弱，导致300～500 km深度由北向南最大主压应力轴方位发生旋转(图 7).根据反演结果，我们同样发现浅部R2、R3区与深部R1—R4区主压应力方位也表现出相对于俯冲带俯冲方向向南偏转的现象(图 5c，7)，说明该南向地幔流对俯冲带北部浅部与深部区域也产生一定影响，并由北向南逐渐减弱.本文根据应力方位变化提供了南向流动的地幔流对俯冲板块产生推挤作用，并且向南逐渐减弱的证据.
反演结果显示该俯冲带最大主张应力轴方位分布较零散(图 5d)，产生这种现象的原因可能是由于部分R值接近于1(表 2，图 6)：张应力轴可能在垂直于压应力轴的平面任意旋转(万永革，2015)；另一方面从动力学角度分析：与上述导致最大压应力轴方位旋转的原因相同，由北向南流动的地幔流对俯冲带产生影响并向南逐渐减弱，导致300～500 km深度区最大张应力轴方位由北向南发生旋转.
6.3 “偏移”板片与主俯冲带的动力作用关系如图 3a、图 4所示，在汤加—克马德克俯冲带500～650 km深度左右存在极其复杂的构造形态，这种现象反映出处于地幔转换带底部(660 km左右)的俯冲带与下地幔复杂的动力作用关系(Giardini and Woodhouse，1984，1986；Okal and Kirby，1998).主俯冲带最深部(图 3a)表现出强烈的褶皱弯曲变形.在主俯冲带西侧，有一簇不连续的呈带状分布的“板外”震群存在(图 3a、4椭圆标记区)，这些震群主要存在于22°S—27°S俯冲带底部倾角较陡的区域内.Okal和Kirby(1998)将其解释为：早期连续向下俯冲的板片遇到黏性较高的下地幔阻挡，停滞在上地幔底部及转换带内，随着演化过程中俯冲板片后撤(Slab Rollback)，造成早期俯冲板片的断裂，导致偏移板片的产生.Lundgren和Giardini(1994)发现汤加俯冲带深部偏移板片内震源机制的P轴接近垂直，认为这是由于在板片后退过程中，偏移板片逐渐脱离主俯冲带俯冲方向挤压环境的结果.
为了研究偏移板片与主俯冲带的应力场差异，本文将(22°S—27°S)R5—R7区明显的偏移震群(图 4灰色虚线椭圆标记区)进行应力张量反演(R4区偏移地震个数较少，为4个)，并与主体俯冲带应力场对比分析(表 3).结果发现R5—R7区偏移板片最大主压应力轴倾伏角为近垂直或近水平，这与各个区域主俯冲带所呈现的应力结构不同(图 8).偏移板片的最大主压应力轴水平或垂直均反映出不再受到上部俯冲带向下推挤的作用，未产生沿俯冲方向挤压的应力结构，说明主俯冲带与偏移板片主应力方向不连续.该应力反演结果为偏移板片的分离提供了证据.
| | 表 3 深部偏移板片与主俯冲带最优应力参数对比 Table 3 The comparison of the offset slab and main subduction zone best fitting stress parameters |
 | 图 8
偏移板片与主俯冲带应力反演结果比较
灰色实线框包围了R5—R7区偏移板片震源机制解，蓝色实线框包围了主俯冲带深部震源机制解，震源机制不同大小代表不同震级同图 4；主应力空间分布图含义同图 6.(a)—(c)分别为R5—R7区偏移震群P、T轴分布与反演出的主应力轴等面积投影，含义同图 5． Fig. 8 Fig. 8 The inverted stress tensor comparison of the offset slab and main subduction zone The offset focal mechanisms from R5—R7 are marked by gray solid lines，the deep focal mechanisms of main zone are marked by blue solid lines. The meaning of scale of beach balls is the same with Fig. 4. The meaning of stress spatial distribution is the same with Fig. 2. Figure(a)—(c)represent the equal area projections of P，T axes distribution and principle stress axes by inversion of R5~R7 deep offset earthquakes respectively，the meaning of stress plane and color is the same with Fig. 5. |
本文首先采用全球EHB目录对汤加—克马德克俯冲带进行线性拟合，得到俯冲带的定量参数，然后根据震源机制的分布特征选取不同区域对其进行应力张量反演，得到汤加—克马德克俯冲带研究区(18.5°S—28.5°S)高分辨率的应力图像.与Bonnardot等(2009)的研究工作相比，本文应力反演区域划分更加精细，反演结果很好地反映了汤加—克马德克俯冲带研究区内应力场的非均匀特征.我们根据该俯冲带应力场的纵向与横向变化，分析了应力场非均匀性所揭示的动力学意义，得到如下认识：
(1)俯冲带浅部(60～300 km)主应力倾伏角变化多样，并且可能受到路易斯维尔海链(LSC)俯冲的影响，使得24°S附近最大压应力方向发生了强烈旋转.
(2)俯冲带中部(300～500 km)目前未受到北斐济盆地下方拆离板片的动力作用影响，并且最大压、张应力轴由北向南的逐渐旋转表明南北向流动的地幔流对俯冲板块北部产生推挤作用，并且这种推挤作用向南逐渐减弱.
(3)俯冲带深部(500～700 km)由于下地幔的阻挡以及上部板片推挤的作用，导致最大主压应力方向沿俯冲方向分布，支持该俯冲带深部水平及垂直的断层面分布的结论(Warren et al.，2007).
(4)深部主俯冲带西侧偏移板片的最大主压应力轴水平或垂直，与深部主俯冲带沿俯冲方向挤压的应力结构明显不同，证实了偏移板片与主俯冲带是分离的.
致谢 感谢中国科学院大学地球科学学院周元泽教授在本文的撰写过程中提出了宝贵意见. 部分图件采用GMT软件(Wessel and Smith，1995)绘制.| [1] | Arnold R, Townend J. 2007. A Bayesian approach to estimating tectonic stress from seismological data. Geophys. J. Int., 170(3):1336-1356. |
| [2] | Balance P F, Scholl D W, Vallier T L, et al. 1989. Subduction of a late Cretaceous seamount of the Louisville Ridge at the Tonga trench:A model of normal and accelerated tectonic erosion. Tectonics, 8(5):953-962. |
| [3] | Benioff H. 1954. Orogenesis and deep crustal structure:additional evidence from seismology. Geol. Soc. Amer. Bull., 65(5):385-400. |
| [4] | Billen M I, Gurnis M, Simons M. 2003. Multiscale dynamics of the Tonga-Kermadec subduction zone. Geophys. J. Int., 153(2):359-388. |
| [5] | Bonnardot M A, Régnier M, Ruellan E, et al. 2007. Seismicity and state of stress within the overriding plate of the Tonga-Kermadec subduction zone. Tectonics, 26(5), doi:10.1029/2006TC002044. |
| [6] | Bonnardot M A, Régnier M, Christova C, et al. 2009. Seismological evidence for a slab detachment in the Tonga subduction zone. Tectonophysics, 464(1-4):84-99. |
| [7] | Chen P F, Bina C R, Okal E A. 2004. A global survey of stress orientations in subducting slabs as revealed by intermediate-depth earthquakes. Geophys. J. Int., 159(2):721-733. |
| [8] | Chen W P, Brudzinski M R. 2001. Evidence for a large-scale remnant of subducted lithosphere beneath Fiji. Science, 292(5526):2475-2479. |
| [9] | Chen W P, Brudzinski M R. 2003. Seismic anisotropy in the mantle transition zone beneath Fiji-Tonga. Geophys. Res. Lett., 30(13):1682, doi:10.1029/2002GL016330. |
| [10] | Christova C, Scholz C H. 2003. Stresses in the Vanuatu subducting slab:A test of two hypotheses. Geophys. Res. Lett., 30(15):1790, doi:10.1029/2003GL017701. |
| [11] | Christova C. 2004. Stress field in the Ryukyu-Kyushu Wadati-Benioff zone by inversion of earthquake focal mechanisms. Tectonophysics, 384(1-4):175-189. |
| [12] | Christova C, Scholz C H, Kao H. 2004. Stress field in the Vanuatu (New Hebrides) Wadati-Benioff zone inferred by inversion of earthquake focal mechanisms:evidence for systematic lateral and vertical variations of principal stresses. Journal of Geodynamics, 37(2):125-137. |
| [13] | Christova C. 2005. Space distribution of the contemporary stress field in the Izu-Bonin Wadati-Benioff zone by inversion of earthquake focal mechanisms. Journal of Geodynamics, 39(4):413-428. |
| [14] | Christova C V. 2015. Spatial distribution of the contemporary stress field in the Kurile Wadati-Benioff zone by inversion of earthquake focal mechanisms. Journal of Geodynamics, 83:1-17. |
| [15] | Conder J A, Wiens D A. 2007. Rapid mantle flow beneath the Tonga volcanic arc. Earth and Planetary Science Letters, 264(1-2):299-307. |
| [16] | DeMets C, Gordon R G, Argus D F, et al. 1990. Current plate motions. Geophys. J. Int., 101(2):425-478. |
| [17] | Dziewonski A M, Chou T A, Woodhouse J H. 1981. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. J. Geophys. Res., 86(B4):2825-2852. |
| [18] | Ekström G, Nettles M, Dziewoński A M. 2012. The global CMT project 2004-2010:Centroid-moment tensors for 13017 earthquakes. Phys. Earth Planet. Int., 200-201:1-9. |
| [19] | Engdahl E E, Van der Hilst R D, Buland R. 1998. Global teleseismic earthquake relocation with improved travel times and procedures for depth determination. Bull. Seism. Soc. Am., 88(3):722-743. |
| [20] | Foley B J, Long M D. 2011. Upper and mid-mantle anisotropy beneath the Tonga slab. Geophys. Res. Lett., 38(2):L02303, doi:10.1029/2010GL046021. |
| [21] | Fukao Y, Obayashi M. 2013. Subducted slabs stagnant above, penetrating through, and trapped below the 660 km discontinuity. J. Geophys. Res., 118(11):5920-5938. |
| [22] | Gephart J W, Forsyth D W. 1984. An improved method for determining the regional stress tensor using earthquake focal mechanism data:application to the San Fernando earthquake sequence. J. Geophys. Res., 89(B11):9305-9320. |
| [23] | Giardini D, Woodhouse J H. 1984. Deep seismicity and modes of deformation in Tonga subduction zone. Nature, 307(5951):505-509. |
| [24] | Giardini D, Woodhouse J H. 1986. Horizontal shear flow in the mantle beneath the Tonga arc. Nature, 319(6054):551-555. |
| [25] | Hall R, Spakman W. 2002. Subducted slabs beneath the eastern Indonesia-Tonga region:insights from tomography. Earth and Planetary Science Letters, 201(2):321-336. |
| [26] | Hamburger M W, Isacks B L. 1987. Deep earthquakes in the Southwest Pacific:a tectonic interpretation. J. Geophys. Res., 92(B13):13841-13854. |
| [27] | Hardebeck J L, Michael A J. 2006. Damped regional-scale stress inversions:Methodology and examples for Southern California and the Coalinga aftershock sequence. J. Geophys. Res., 111:B11310, doi:10.1029/2005JB004144. |
| [28] | Hayes G P, Wald D J, Johnson R L. 2012. Slab1.0:A three-dimensional model of global subduction zone geometries. J. Geophys. Res., 117:B01302, doi:10.1029/2011JB008524. |
| [29] | Huang J C, Wan Y G. 2015. Present Tectonic stress field in the capital region of China determined from small and strong earthquake focal mechanisms. Earthquake (in Chinese), 35(1):17-27. |
| [30] | Isacks B, Molnar P. 1971. Distribution of stresses in the descending lithosphere from a global survey of focal-mechanism solutions of mantle earthquakes. Rev. Geophys., 9(1):103-174. |
| [31] | Jiao W J, Silver P G, Fei Y W, et al. 2000. Do intermediate- and deep-focus earthquakes occur on preexisting weak zones? An examination of the Tonga subduction zone. J. Geophys. Res., 105(B12):28125-28138. |
| [32] | Lu R Q, Suppe J, He D F, et al. 2013. Deep subducting slab reconstruction and its geometry, kinematics:a case study for the Tonga-Kermadec slab from tomography. Chinese J. Geophys. (in Chinese), 56(11):3837-3845, doi:10.6038/cjg20131125. |
| [33] | Lundgren P, Giardini D. 1994. Isolated deep earthquakes and the fate of subduction in the mantle. J. Geophys. Res., 99(B8):15833-15842. |
| [34] | Lü J, Wang X S, Su J R, et al. 2013. Hypocentral location and source mechanism of the MS7.0 Lushan earthquake sequence. Chinese J. Geophys. (in Chinese), 56(5):1753-1763, doi:10.6038/cjg20130533. |
| [35] | Meighan H E, Pulliam J, Brink U, et al. 2013. Seismic evidence for a slab tear at the Puerto Rico Trench. J. Geophys. Res., 118(6):2915-2923. |
| [36] | Michael A J. 1987. Use of focal mechanisms to determine stress:A control study. J. Geophys. Res., 92(B1):357-368. |
| [37] | Okal E A, Kirby S H. 1998. Deep earthquakes beneath the Fiji Basin, SW Pacific:Earth's most intense deep seismicity in stagnant slabs. Phys. Earth Planet. Int., 109(1-2):25-63. |
| [38] | Richards S, Holm R, Barber G. 2011. When slabs collide:A tectonic assessment of deep earthquakes in the Tonga-Vanuatu region. Geology, 39(8):787-790. |
| [39] | Ruellan E, Delteil J, Wright I, et al. 2003. From rifting to active spreading in the Lau Basin-Havre Trough backarc system (SW Pacific):Locking/unlocking induced by seamount chain subduction. Geochem. Geophys. Geosyst., 4(5):8909, doi:10.1029/2001GC000261. |
| [40] | Sheng S Z, Wan Y G, Wang W L, et al. 2014. The fault plane parameter determination of the 2010 Yushu MS7.1 earthquake. Progress in Geophysics (in Chinese), 29(4):1555-1562, doi:10.6038/pg20140409. |
| [41] | Smith G P, Wiens D A, Fischer K M, et al. 2001. A complex pattern of mantle flow in the Lau backarc. Science, 292(5517):713-716. |
| [42] | Syracuse E M, Abers G A. 2006. Global compilation of variations in slab depth beneath arc volcanoes and implications. Geochem. Geophys. Geosyst., 7(5):Q05017. |
| [43] | Turner S P, Hawkesworth C G. 1998. Using geochemistry to map mantle flow beneath the Lau basin. Geology, 26(11):1019-1022. |
| [44] | Van der Hilst R D. 1995. Complex morphology of subducted lithosphere in the mantle beneath the Tonga trench. Nature, 374(6518):154-157. |
| [45] | Van der Hilst R D, Widiyantoro S, Engdahl E R. 1997. Evidence for deep mantle circulation from global tomography. Nature, 386(6625):578-584. |
| [46] | Wada I, Mazzotti S, Wang K L. 2010. Intraslab stresses in the Cascadia subduction zone from inversion of earthquake focal mechanisms. Bull. Seism. Soc. Am., 100(5A):2002-2013. |
| [47] | Wadati K. 1935. On the activity of deep-focus earthquakes in the Japan islands and neighborhoods. Geophys. Mag., 8:305-325. |
| [48] | Wan Y G, Shen Z K, Diao G L, et al. 2008. An algorithm of fault parameter determination using distribution of small earthquakes and parameters of regional stress field and its application to Tangshan earthquake sequence. Chinese J. Geophys. (in Chinese), 51(3):793-804. |
| [49] | Wan Y G, Sheng S Z, Hsu Y J, et al. 2011. Effect of stress ratio and friction coefficient on composite P wave radiation patterns. Chinese J. Geophys. (in Chinese), 54(4):994-1001, doi:10.3969/j.issn.0001-5733.2011.04.014. |
| [50] | Wan Y G, Sheng S Z, Huang J C, et al. 2015. The grid search algorithm of tectonic stress tensor based on focal mechanism data and its application in the boundary zone of China, Vietnam and Laos. Journal of Earth Science,online,doi:10.1007/s12583-015-0649-1. |
| [51] | Wan Y G. 2015. A grid search method for determination of tectonic stress tensor using qualitative and quantitative data of active faults and its application to the Urumqi area. Chinese J. Geophys. (in Chinese), 58(9):3144-3156, doi:10.6038/cjg20150911. |
| [52] | Warren L M, Hughes A N, Silver P G. 2007. Earthquake mechanics and deformation in the Tonga-Kermadec subduction zone from fault plane orientations of intermediate- and deep-focus earthquakes. J. Geophys. Res., 112(B5):B05314, doi:10.1029/2006JB004677. |
| [53] | Wessel P, Smith W H F. 1995. New version of the generic mapping tools released. Eos Trans. AGU, 76:329. |
| [54] | Zang S X, Wu Z L, Ning J Y. 1992. Stress state and deep deformation of subduction zone in Japan sea and the sea of Okhotsk and the problem of slab penetration into the lower mantle. Chinese J. Geophys. (in Chinese), 35(5):560-573. |
| [55] | Zhang L P, Shao Z G, Ma H S, et al. 2013. The plate contact geometry investigation based on earthquake source parameters at the Burma arc subduction zone. Science China:Earth Sciences, 56(5):806-817. |
| [56] | Zhang L P, Shao Z G, Li Z H. 2014. Dynamic action of mutual subduction between the Indian and the Eurasia plates in Hindu Kush-Pamir region. Chinese J. Geophys. (in Chinese), 57(2):459-471, doi:10.6038/cjg20140212. |
| [57] | 黄骥超, 万永革. 2015. 利用小震与强震震源机制解反演首都圈现今构造应力场. 地震, 35(1):17-27. |
| [58] | 鲁人齐, Suppe J, 何登发等. 2013. 深部俯冲板片三维构造重建及其几何学、运动学研究-以汤加-克马德克地区俯冲板片为例. 地球物理学报, 56(11):3837-3845, doi:10.6038/cjg20131125. |
| [59] | 吕坚, 王晓山, 苏金蓉等. 2013. 芦山7.0级地震序列的震源位置与震源机制解特征. 地球物理学报, 56(5):1753-1763, doi:10.6038/cjg20130533. |
| [60] | 盛书中, 万永革, 王未来等. 2014. 2010年玉树MS7.1地震发震断层面参数的确定. 地球物理学进展, 29(4):1555-1562, doi:10.6038/pg20140409. |
| [61] | 万永革, 沈正康, 刁桂苓等. 2008. 利用小震分布和区域应力场确定大震断层面参数方法及其在唐山地震序列中的应用. 地球物理学报, 51(3):793-804. |
| [62] | 万永革, 盛书中, 许雅儒等. 2011. 不同应力状态和摩擦系数对综合P波辐射花样影响的模拟研究. 地球物理学报, 54(4):994-1001, doi:10.3969/j.issn.0001-5733.2011.04.014. |
| [63] | 万永革. 2015. 联合采用定性和定量断层资料的应力张量反演方法及在乌鲁木齐地区的应用. 地球物理学报, 58(9):3144-3156, doi:10.6038/cjg20150911. |
| [64] | 臧绍先, 吴忠良, 宁杰远. 1992. 日本海-鄂霍次克海下俯冲带的应力状态及其深部形变. 地球物理学报, 35(5):560-573. |
| [65] | 张浪平, 邵志刚, 马宏生等. 2013. 基于地震参数的缅甸弧俯冲带处板块间几何接触方式的研究. 中国科学:地球科学, 43(4):653-664. |
| [66] | 张浪平, 邵志刚, 李志海. 2014. 印度板块与欧亚板块在兴都库什-帕米尔地区相互俯冲的动力作用分析. 地球物理学报, 57(2):459-471, doi:10.6038/cjg20140212. |