地球物理学报  2015, Vol. 58 Issue (7): 2350-2365   PDF    
台湾碰撞带现今地壳变形场特征及其动力学成因的有限单元法模拟研究
龙小刚1, 朱守彪1,2    
1. 中国地震局地壳应力研究所(地壳动力学重点实验室), 北京 100085;
2. 中国科学院计算地球动力学重点实验室, 北京 100049
摘要:台湾是菲律宾海板块与欧亚板块会聚、碰撞的产物,区内地质构造复杂,地震频发,变形强烈.为定量研究台湾地区变形特征及其动力学成因,文中运用二维不连续变形体弹性力学的有限单元计算方法,利用台湾地区1995-2005年GPS的观测结果作为边界约束,对台湾地区地壳变形场进行了模拟.计算结果显示,台湾中部地区主要呈压缩状态,但在台湾的东北部及南部地区出现了拉张的变形环境.变形程度最大的区域位于台湾东部海岸山脉及其附近海域,同时,计算给出台湾纵谷断裂滑移速率约为13.8~23.5 mm/yr,由于纵谷断裂对变形的吸收,因此变形在纵谷断裂以西及西北地区迅速衰减.此外,计算结果还发现,计算给出的台湾岛上的速度值与GPS观测结果吻合得较好;计算给出的主应力方位与地应力观测及震源机制解结果也颇为一致.由此说明文中的有限元模型是合理的.此外,计算结果还表明,菲律宾海板块向北西方向的推挤,板块边界形状,观音、北康高地的抵阻, 冲绳海槽扩张及琉球海沟的向南后撤以及断裂作用等,共同造成了台湾地区现今变形场的主要格局.
关键词GPS观测     变形场     动力学成因     有限单元     台湾    
Present-day crustal strain field in the Taiwan active collision zone and its geodynamic mechanisms: Insight from FEM simulations
LONG Xiao-Gang1, ZHU Shou-Biao1,2    
1. Key Laboratory of Crustal Dynamics, Institute of Crustal Dynamics, China Earthquake Administration, Beijing 100085, China;
2. Key Laboratory of Computational Geodynamics, Chinese Academy of Science, Beijing 100049, China
Abstract: The Taiwan Island is the product of convergence and collision between the Eurasian plate and the Philippine Sea plate, where the geological structure is complex, seismicity is high and deformation is strong. As ongoing arc-continent collision around Taiwan is rare in the world, Taiwan has been described as an ideal natural laboratory for plate collision studies. Taiwan has established its GPS Network since the 1980s, which has cumulated mass data since then and makes it possible to study the crustal deformation in and around it.
In order to quantitatively study the characteristics of crustal deformation in Taiwan and to understand their geodynamic mechanisms, we calculated the strain rate field in and around Taiwan by using finite element method (FEM), utilizing GPS data from 1995 to 2005 as boundary constraints in simulation. The heterogeneity of the medium, slippage of the faults and some other factors are taken into consideration to quantitatively study the interseismic deformation in Taiwan and we try to find out the dominant factors that control the deformation. To acquire more ideal computed results, we have taken advantage of the optimal design function of the ANSYS(Design Opt). We take the fitted difference between GPS observations and their corresponding finite element calculated results as objective function, take the material parameters of each subregion and friction coefficients of the faults as design parameters to obtain the optimal finite element model.
The computed results of the optimal model show that the largest calculated maximum principal strain of Taiwan is located in the Coastal Range and adjacent waters in eastern Taiwan, the superiority trend of maximum principal compressive strain is NW-SE, which agrees with the orientation of arc-continent collision between the Eurasian plate and the Philippine Sea plate and converges obliquely to Taiwan Island. The computed results also show that contractions are predominant in the central part of Taiwan, while extensions exist in the northeastern and southern parts, respectively. The principal tensile strain rate is larger than the principal compressive strain rate in I-lan Plain, northeast Taiwan, indicating that the region is in the tensile stress environment which may be closely related to the extension of the Okinawa trough. Similarly, tensile stress environment also appears in southern Taiwan. In addition, the azimuths of calculated principal tensile strain are mainly consistent with the directions of the world stress map results and the P and T axes of focal mechanism solutions.The computed velocity fields in Taiwan show a fan shape. The displacement vectors in central and eastern part are in the direction of NW, the direction of velocity changes clockwise towards north and anticlockwise towards south, which reflect the characteristics of Taiwanese materials escaping to northeast and southwest. The Philippine Sea Plate squeezes Taiwan in the direction of NW, but the movement to northwest is blocked by the Kuanyin High(KH) and Peikang High(PH). This may result in the escape phenomenon. In general, the magnitudes of velocities gradually decrease northwestward from about 81.5 mm/yr at Lanyu and Lutao, SE offshore of Taiwan, to nearly no deformation in NW Taiwan. Besides, the simulation results show that most faults are locked(or with large friction coefficients) during interseismic deformation, while the calculated friction coefficient of the Longitudinal Valley Fault (LVF) is 0.5, and the slip rate is about 13.8~23.5 mm/yr. As part of the convergence is absorbed by LVF, deformation west of LVF decays rapidly westwards and northwestwards.
Because of the special geology structure in Taiwan region, the crustal strain field is complicated. The calculated results imply that the general framework of present-day deformation in Taiwan results from interactions of many factors such as plate collision between the Eurasian plate and the Philippine Sea Plate, geometry of the plate boundaries, the obstruction of Kuanyin and Peikang high, faulting and rifting, opening of the Okinawa Trough and retreat of the Ryukyu Trench.
Key words: GPS observation     Deformation field     Geodynamic mechanism     Finite element method     Taiwan    
1 引言

台湾位于欧亚板块与菲律宾海板块的会聚边界上(见图 1),是由吕宋岛弧与欧亚大陆边缘会聚、碰 撞形成的(Ho,1986; Teng,198719901996; Huang et al.,2006). 这种碰撞至少开始于8百万年前,并且现在仍在进行之中(Ho,1988; Teng,1990; Kao et al.,1998; Lallemand et al.,2001; Bos et al.,2003).台湾及其周边地区的地质构造极其复杂:在其东北方向,菲律宾海板块沿琉球海沟向西北方向俯冲到欧亚板块之下;同时,在西南方向,欧亚板块沿着马尼拉海沟向东南俯冲到菲律宾海板块之下(Wu,1978; Tsai,1986; Pezzopane and Wesnousky,1989). 台湾处于现代沟-弧-盆体系内,是世界上少见的“弧陆碰撞”正在进行的区域,是研究板块碰撞的天然实验场(Kao and Chen,2000; Chiang et al.,2008).由于正在进行的板块俯冲和碰撞作用,有记录以来台湾地区发生了许多破坏性地震(Bonilla,19751977; Hsu,1980; Chi et al.,1981; Chen and Tsai,2008).例如:1994年9月16日发生在台湾西边的台湾海峡地震(MW=6.8),1999年9月21日发生在台湾中部的集集地震(MW=7.6)以及2006年12月26日的屏东地震(MW=7.0)(Kao and Wu,1996; Chang et al.,2000; Wu et al.,2009; 孙玉军等,2011).图 2展示了1900年1月1日-2014年12月31日期间,台湾地区震源深度小 于40 km,震级为5≤M≤8的 地震震中分布(图中 数据来源于http://earthquake.usgs.gov/earthquakes/search/,最后访问时间2015年1月1日).研究台湾地区的地壳形变特征、活动断层及其滑动速率、动力学机制对于研究台湾地区地震成因及分析未来潜在地震危险性有重要的科学意义.

图 1 台湾地区的地质背景 Fig. 1 Geological Setting of the area around Taiwan

图 2 台湾地区1900-2014年震中分布(震级5≤M≤8,震源深度小于40 km)和1995-2005年GPS速度场分布图.其中A-E是台湾地区主要活动构造,分别为变 形前缘,车笼铺断层,屈尺断层,梨山断层和纵谷断层(Ching et al.,2011) Fig. 2 The epicenter distribution in and around Taiwan from 1900 to 2014(5≤M≤8,with focal depth less than 40 km) and the GPS horizontal velocity field in Taiwan from 1995 to 2005 with respect to station S01R(black solid triangle)in the continental margin of the Penghu Isl and at the Taiwan Strait. A-E indicate the main active structures in Taiwan,the deformation front,the Chelungpu fault,the Chihtsu fault,the Lishan fault,and the Longitudinal Valley fault

台湾地区自20世纪80年代以来建立GPS观测系统,陆续累积记录了大量关于地壳运动的观测资料,这些数据为研究整个台湾地区的地壳变形提供了可能.Yu等(1997)利用台湾GPS观测网络记录结果,计算给出了1990-1995年台湾地区地壳运动的速度场.利用台湾的这些GPS观测资料,许多学者研究了台湾地区的变形场,如:何玉梅和姚振兴(2002)研究了台湾南部地区的地壳变形,给出了菲律宾海板块与欧亚板块在台湾地区的会聚过程中,有一半左右的变形量被台东纵谷吸收,在中央山脉以西,块体运动大致呈扇形分布,并与应力方向一致;Bos等(2003)利用GPS资料综合反演了变形及断层滑移速度,给出了整个台湾的应变场;郝金来等(2009)利用台湾地区GPS观测速度反演得到了该地区的应变率场及块体的旋转速率;Lin等(2010)利用2003-2005年的连续GPS数据对台湾地表变形场进行计算;朱守彪等(2011)利用集集地震震前(1990-1995年)及震后(2003-2005年)的GPS观测资料分别对台湾及附近地区的地应变率场进行了计算,给出了集集地震前后的变形场特征;Ching等(2011)利用1995-2005年的GPS数据(见图 2)计算构造块体的运动和断层滑动速率进而分析评估了潜在大地震发生的可能性.

在此基础上,为定量研究台湾地区地壳变形的力学成因,Huchon等(1986)利用有限元方法计算了受板块推挤碰撞引起的台湾地区的应力应变场,给出了主应力迹线的分布格局;Hu等(199619972001)利用有限元和离散元方法模拟了典型地质构造对台湾南部变形场的影响及台湾地区速度场,给出了台湾地区速度场的分布及其主要的影响因素,并分析了台湾地区挤压、形变与应力分布之间的关系.

但是,尽管前人对台湾地区的变形场进行了许多研究,给出了很多有意义的结果,但利用最新、更为全面的GPS观测资料,考虑断层的解耦等更多因素来定量研究台湾地区的变形特征,近年来还未见到.所以,本文利用有限元方法,综合利用台湾地区1995-2005年的GPS观测资料(图 2),考虑介质的不均匀性及断层的滑移等多种复杂因素,定量模拟台湾地区的震间变形,并力图给出决定变形的主要控制因素.

2 地质背景

台湾造山带及其频繁的地震活动是菲律宾海板 块上的吕宋火山岛弧与欧亚板块的大陆边缘之间倾斜碰撞的结果(Suppe,1984; Yu et al.,1997).欧亚板块与菲律宾海板块之间的会聚速率可达81.5 mm/yr(Yu et al.,19971999),台湾地区的运动框架及震中分布特征与两个俯冲系统有关,一个是台湾以东的琉球沟弧系统,另一个是台湾以南的吕宋沟弧系统(Hsu,1990; Rau and Wu,1995; Ching et al.,2011)(见图 1).在琉球沟弧系统中,菲律宾海板块向北沿琉球岛弧俯冲至欧亚板块之下,而在吕宋沟弧系统中,菲律宾海板块向西沿吕宋岛弧覆盖于欧 亚板块之上(Wu,1978; Tsai,1986; Pezzopane and Wesnousky,1989). 琉球沟弧系统从琉球南部一直延伸到台湾东部的花莲附近,与之相关的弧后扩张形成的冲绳海槽向西延伸到台湾的东北部宜兰平原(Ho,1986; Tsai,1986; Hu et al.,1996).吕宋沟弧系统从菲律宾延伸到台湾西南部高雄附近,在这组成了台湾山脉的一部分(Hu et al.,1996).位于琉球岛弧与吕宋岛弧之间的纵谷断层(Longitudinal Valley Fault)通常被认为是吕宋岛弧与欧亚大陆边缘正在进行的碰撞的缝合带,也是台湾东部主要的 活动构造(Biq,1972; Hsu,1976; Wu,1978; Barrier and Angelier,1986; Ho,1986; Angelier et al.,2000; Lee et al.,2001). 纵谷断层是一个左旋逆冲断层,左旋与逆冲分量的比值为1 ∶ 3(Yu and Liu,1989; Lee and Angelier,1993a),它吸收了菲律宾海板块向西北挤压速率的40%,使海岸山脉和中央山脉之间的相对运动速率达到28~31 mm/yr(Yu and Chen,1994).

在台湾地区断层形成的过程中,台湾西边的北康高地和观音高地的作用十分重要(图 3),它们代表前白垩纪中国大陆浅层基底,构造比较稳定(Bos et al.,2003),由于北康高地相对于台湾地区来说是个坚硬的基底,阻止相邻的构造带继续往西运动,因此处在马尼拉消减带和北康高地之间的台湾西南 地区的物质被迫向西南方向运动(Lu and Malavieille,1994; Lu et al.,1998; Lin and Watts,2002). 北康高地和观音高地的阻抵及菲律宾海板块向西北的推挤,使得台湾地区的山脉基本呈北北东走向拉长的“S”型构造(Biq,1972; Lu and Malavieille,1994; Lu et al.,1998; Chang et al.,2003;Ching et al.,2011).

图 3 有限元模型结构 其中Ⅰ为一组断层(黑色粗线),包括:纵谷断层和吕宋岛弧北缘.1-13为模型中使用的断层模型:1-美仑断层;2-纵谷断层;3-吕宋岛弧北沿;4-梨山断层;5-潮州-旗山断层;6-恒春断层;7-屈尺断层;8-右昌断层;9-左镇断层;10-大尖山-触口断层;11-车笼辅断层;12-台湾西北部断层;13-变形前缘.模型被断层(组)划分成的子区域见图 4a. 由空 心箭头表示菲律宾海板块相对欧亚大陆板块的运动速率.蓝色实心箭头表示物质逃逸方向(Angelier et al.,2009). Fig. 3 Map of finite elemtne structures Ⅰ(the bold black line)in the figure is a fault group,in which LVF and north margin of the Luzon arc are included. 1-13 are fault models implied in simulation: 1 Meinung Fault; 2 LVF; 3 North margin of the Luzon arc; 4 Lishan Fault; 5 Chaochou-Chishan Fault; 6 Hengchun Fault; 7 Chuchih Fault; 8 Youchang Fault;9 Cuochen Fault; 10 Choukou Fault; 11 Chelungpu Fault; 12 Faults at Northern Taiwan; 13 Deformation Front. The study area is divided into several subregions by those faults above,see Figure 4a. Large open arrow indicating the direction of plate convergence of Philippine Sea plate relative to Eurasia,the number inside it is the convergence rate. The blue solid arrows indicate extrusion and tectonic escape to Ilan Plain and Pingtung Plain areas.

台湾岛通常被主要的断层及逆冲断裂划分为几个主要的地质单元,各单元基本呈北北东走向,由西向东依次是彭湖列岛、海岸平原、西部山麓、雪山山脉、中央山脉、纵谷断层和海岸山脉(Ho,1986;Yu et al.,1997)(由图 3所示).纵谷断层是海岸山脉与中央山脉之间的分界线,潮州-旗山断层是中央山脉南部与西部山麓之间的分界,梨山断层是中央山脉与雪山山脉之间的界限,海岸平原与西部山麓之间的分界是台湾地区变形前缘(图 3)(Bos et al.,2003).彭湖列岛位于台湾海峡,主要由更新世泛布玄武岩构成;海岸平原主要成分为来自中央山脉和西部山麓的第四纪河流沉积物;西部山麓主要由从晚渐新世、中新世到早更新世以来的浅海沉积物厚层经过推覆构造断层活动挤压变形形成;雪山山脉的主要成分是已经变质为绿片岩相的始新世和渐新世大陆边缘海洋沉积物;中央山脉西部是新生代变质硅质黏土岩,东部由前第三纪基底杂岩、晚第三纪绿片岩和中生代到新生代多相变质岩组成;纵谷断层是台湾地区重要的剪切变形带,由第四纪河流沉积物填充;海岸山脉以及绿岛和兰屿两个小岛是北吕宋岛弧的一部分,由晚第三纪安山岩组成(Yu et al.,1997; Hsu et al.,2009; Ching et al.,2011).

台湾北部是一个转换带,在这里菲律宾海板块由覆盖在向东消减的欧亚板块之上转换为向北沿琉球岛弧俯冲到欧亚板块之下,因此台湾北部已经运动到琉球沟弧系统之中而不再处在碰撞边界上了(Suppe,1984; Teng,1996; Clift et al.,2003; Rau et al.,2008; Ching et al.,2011). 东北部的宜兰平原在冲绳海槽的最西南端(Ho,1986; Tsai,1986; Hu et al.,1996),其北边是雪山山脉,南边是中央山脉(图 3),呈三角形状,由于冲绳海槽的扩展(opening),宜兰平原处在拉张的变形环境中,因此平原附近多发生正断层型地震(Tsai,1986).

台湾西南部的屏东平原呈南北走向,其东边是中央山脉,北边和西边是西部山麓,而南边是马尼拉 增生楔(accretionary wedge)(Angelier et al.,2009). 台湾南部主要是马尼拉增生楔在岛上的延伸部分,相对于马尼拉消减带来说是一个软弱带,马尼拉增生楔的北部边界为旗山断层,西部边界为变形前缘,通常认为由菲律宾海板块和欧亚板块碰撞引起的变形都分布在变形前缘以东(Bos et al.,2003).

3 有限元模型

为定量研究台湾地区的变形场,并探索断层等因素对变形场的影响,文中用数值模拟的方法来进行定量分析.为抓住主要矛盾,揭示变形场的主要特征,按照朱守彪和石耀霖(200420052007)以及石耀霖和朱守彪(2006)的作法,取地下5 km的平面 作为研究平面来建立二维有限元模型.根据台湾地区的地质构造(图 1),文中选定研究区域为119.5°E-122.3°E,21.8°N-25.5°N,这样整个台湾岛处在研究区域内.由于台湾地区构造复杂,断层发育,为了更好地模拟变形的主体特征,研究中考虑台湾地区断层形成的运动学界面对变形场的影响.根据台湾地区断层的几何特征及前人的研究(Ho,1986;郑世楠,1995; Hu et al.,2001; 何玉梅和姚振兴,2002; Bos et al.,2003; Ching et al.,2007; 郝金来等,2009),在有限元模型中将断层画成直线或简单曲线来近似,具体断层及有限元空间子区域的分块情况见图 3图 4a

图 4 (a)有限元模型子区域划分,图中子区域19,20,21代表纵谷断层内的软弱区域;其他的子区域均由台湾地区各断层几何(见图 3)与台湾岛轮廓划分而得.(b)模型的几何、边界条件及有限元网格.模型采用自由网格划分,最大边长限定为2 km.图(b)中边界上的黑色三角形表示固定边界,空心箭头为边界上所施加的位移的示意图,其中DE-EF段施加的位移边界为81.5 mm/yr,306° Fig. 4 (a)Subregions of the finite element model. Subregions 19-21 represent the weak zone of LVF. Other subregions are devided by faults and Taiwan Isl and ′s outlines.(b)Model geometry,boundary conditions and finite element meshes. The model is free meshed and the length of the elements is no longer than 2 km. The solid triangles on the boundaries indicate fixed boundary condition,the open arrows show the displacement boundary conditions. The displacement on southeast boundaries DE and EF is 81.5 mm/yr with azimuth 306°

图 3所示,以台湾地区断层作为连续或不连续的界面确立了模型结构之后,再加入台湾岛的轮廓,把整个研究区域划分为如下25个小块体(图 4a),其中区域19-21是考虑到断层区上物质属性与周边不同来划分的,其他区域由断层和台湾岛的轮廓划分.如上文所述,台东纵谷断层对台湾的变形场起着重要作用,但其几何目前还不是很清楚,Lee等(2003)Huang等(2010)认为纵谷北部倾角很大,Ching等(2011)取纵谷倾角为70°~75°,深度为20 km.由于本文为二维模型,为更好地模拟断层的作用,取纵谷断层的宽度为4 km,其中被较柔软的物质所充填.

模型中介质的物质属性是非常重要的物理参数,由于本文模拟震间较短的时间尺度,并且研究平面位于上地壳,所以模型中介质视为线弹性.参考Hu等(1996)对台湾及周边地区的模拟、地震层析成像以及前人反演的台湾地区物质属性结果(Lin et al.,1998; Hu et al.,2001;朱守彪和蔡永恩,2006; Wu et al.,2007; 朱守彪和蔡永恩,2009; Kuo-Chen et al.,2012),模型中杨氏模量的取值范围为1~70 GPa,泊松比为0.20~0.30,摩擦系数为0.05~1.20(断层组Ⅰ上).

图 4b是研究中选用的模型几何及有限元网格剖分.模型的几何尺度为310.8 km×410.7 km,有限元模型全部采用三角形单元,节点数为53831,单元总数为106874个.

计算中,模型中的断层利用有限元中的接触单元技术来模拟,对各界面在计算过程中接触行为的判定,使用了库仑摩擦模型(Zhu and Zhang,20102013戴黎明等,2013),其形式如下:

上式中,τ 表示接触面上的剪应力,p表示接触面间的正应力,μ 表示接触面上的摩擦系数.在模拟过程中,接触面之间的作用力可以分解为法向和切向,法向力垂直于接触面,而切向力平行于接触面,公式(1)中的剪应力的相对大小则控制断层面的接触状态:当切向力小于剪切力时,断层面表现为接触粘合状态;当切向力大于剪切力时,断层面表现为接触滑动状态.一对接触面由接触面和目标面组成,为了防止接触面和目标面之间相互穿透,本文利用罚函数法来控制法向约束,此时法向接触力为

式中f为接触力,n为接触面外法线方向,En为罚因子,gn为法向穿透距离,且 gn=n·(xs-xc),其中xsxc分别为节点在接触面上的位置和其在目标面上的投影点的位置.为了防止穿透,通常把En设置成一个很大的值(朱守彪 et al.,2008).

模拟中,边界条件决定着模拟结果.本文的边界条件主要是根据GPS观测及前人的研究来选定(Hu et al.,19972001).台湾处在中国大陆边缘,其西北是稳定的大陆,西边是观音高地和北康高地,相对比较稳定,而且台湾地区GPS数据也是以位于台湾海峡的白沙台为基准,因此模型中北部和西部边界的大部分区段(图 4b中AB和HA段)视为固定边界,而处在菲律宾海板块部分的边界则为均匀边界,速率大小为81.5 mm/yr,方位角为306°,这与 Hu等(2001)模拟中的处理方法相同.模型北部边界东段(BC段),东部边界北段(CD段),南部边界西段(FG段)和西部边界南段(GH段)约束条件是通过GPS观测数据插值来获得的.为防止利用台湾岛内GPS数据外插产生的边界所带来的差错,文中在利用GPS速度场选定模型边界时,还吸纳了研究区域外围的吕宋岛弧和琉球岛弧上另外10个台站的GPS数据(Yu et al.,1999).经处理后这10个台站以白沙台为基准,由于这10个台站均处于研究区域以外,图 2中无法显示,详见表 1.

表 1 附加台站GPS速度观测结果 Table 1 Additional GPS station velocities

本文的模拟计算是利用大型有限元商业软件ANSYS进行的,模型设计及网格剖分等都是利用ANSYS中APDL脚本编写指令来实现.最后,将后处理结果的数据导出,利用GMT等软件来进行绘图展示.

4 模拟结果

模拟中,为获得理想的计算结果,本文利用ANSYS中的优化设计功能(Design Opt)(类似于地球物理中的反演技术).优化设计本身是一种确定结构最优设计方案的技术,目的是让方案满足所有的设计要求的同时,又让方案的某些方面(如模型结构、材料参数等)最为合理.为满足设计要求,需设置设计变量和状态变量要满足的条件,让方案的某些方面最合理,即让目标函数最小(张朝晖,2008).

设计变量(DVs)为自变量,优化结果的取得就是通过改变设计变量的数值来实现的.每个设计变量都有上下限,即要给定每个设计变量的变化范围.ANSYS优化程序允许用户定义60个设计变量.状态变量(SVs)是约束设计变量的数值,是设计变量的函数,是因变量.目标函数(Objective)是要尽量减小的数值,它是设计变量的函数.ANSYS中设计变量和目标函数是必须要有的,而状态变量不是必须的(张朝晖,2008).

为了尽量减少设计变量的个数,参考Hu等(2001)的处理办法,将研究区域所分块体1、2、24和25(见图 4a)的杨氏模量和泊松比分别设置为60GPa和0.25.另外,根据研究区域内的地质背景,把块体划分为若干个组,每组内块体的杨氏模量和泊松比相同(表 2),把断层分为两组,组内各断层摩擦系数相同,而且根据Bos等(2003)的结果及Hu等(2001)的处理方法,把除断层组Ⅰ外的其他的断层上的摩擦系数设置为100.因此,在确定了模型结构、边界条件之后,我们以模型中各块体组的杨氏模量和泊松比以及各断层组上摩擦系数为设计变量,共设置了15个设计变量(7个块体组的杨氏模量及泊松比共14个设计变量,断层组有1个,即摩擦系数有1个设计变量),分别以前文所述各子区域的杨氏模量、泊松比以及断层上的摩擦系数的范围为相应设计变量的上下限,以速度计算值和GPS观测值之间的拟合差为目标函数,进行优化设计.经过大量数值模拟实验,最后获得了拟合差最小的最优模型,也得到了所有子区域材料参数(见表 2)及断层组Ⅰ上摩擦系数值0.5.

表 2 ANSYS优化获得的各子区域材料参数 Table 2 The ANSYS optimized material parameters of all subregions

从优化结果可以看出,代表断层软弱带的区域19-21,杨氏模量是最小的,与实际情况相吻合.断层组Ⅰ上的摩擦系数为0.5,可见该组断层是能够滑移的.详细分析见下文.

4.1 计算给出的台湾地区速度场及纵谷断层滑移特征

上述的ANSYS优化不仅能给出模型参数,同时也给出了模型中各节点的位移结果(或称为速度).图 5是计算给出的均匀节点处的速度值.由图可见,速度场的方向由北向南基本呈扇形分布,中东部速度场的方向为北西向,往北速度场方向发生了顺时针旋转,向南则为逆时针旋转.这是由于菲律宾海板块向西北方向挤压台湾,而台湾西边是凸出的北康高地和观音高地(图 3),向西的运动受到北康及观音高地的阻挡,造成台湾地区物质向东北和西南两个方向逃逸的结果(Lu and Malavieille,1994; Lu et al.,1998; Lin and Watts,2002).总体来说,速度场大小由东南向西北逐渐减小,速度从东南部绿岛和兰屿附近的81.5 mm/yr一直减小到西北部台北盆地的接近于0(见图 5).在海岸山脉地区,速度场大小为37.73~58.86 mm/yr,方位角为303°~323°.从海岸山脉往西跨过纵谷断层,中央山脉西部的速度场大小降为19.62~48.30 mm/yr,方位角为285°~310°,这与Lin等(2010)给出的结果是一致的.跨过纵谷断层速度场大小有较大的衰减,这可能是因为海岸山脉中断层的逆冲运动吸收了一部分水平向运动导致的.西部山麓速度场从北往南递增,大 小为15.09~31.69 mm/yr,方位角为275°~318°.西部海岸平原地区,速度场大小为9.06~13.58 mm/yr,方位角在260°~322°之间.

图 5 模型速度场计算结果和断层Ⅰ两侧速度矢量对比.其中蓝色、红色和黄色箭头均为速度场计算结果,红色和黄色箭头分别表示的是断层Ⅰ上西侧和东侧对应 点处的速度矢量,黑色箭头表示的是断层Ⅰ上的位错 Fig. 5 The calculated results of velocity field and the velocity vectors comparison between two sides of fault group Ⅰ. Computed velocity vectors as blue,red and yellow arrows,displacements on fault group Ⅰ as black arrows. The red and yellow arrows indicate the velocities at the corresponding points on both sides of fault group Ⅰ

在台湾南部,中央山脉南部和恒春半岛地区,速度场大小为49.81~58.86 mm/yr,方位角为275°~280°,变化不是很大.在高雄-屏东平原,速度为43.77~49.81 mm/yr,方向由西逆时针转向西南.这一地区的速度场逆时针旋转可能是由于台湾受到菲律宾海板块的推挤,向西北方向运动,但是受到西边北康高地的阻挡,不能继续往西北运动,于是台湾南部物质被迫向其西南边的较软弱的马尼拉增生楔逃逸,再加上吕宋岛弧往西扩张,增强了物质向西南的逃逸,所以台湾西部速度方向由北西转向正西,台湾西南部则由西转向西南.

丰滨、花莲以北的台湾北部地区的速度场以梨山断层为界可以划分为东西两个主要部分.梨山断层以西的雪山山脉等地区,速度场方向基本与菲律宾海板块的推挤方向一致,方位角大于315°,在雪山山脉北部方位角接近正北,在台北盆地顺时针旋转到北东方向,大小为0.2~7.55 mm/yr,台北盆地西边和北边,速度场都很小,这与Rau等(2008)得出的相同地区速度场大小0.2~11.5 mm/yr的结论是接近的.这可能是因为台湾北部的西边是观音高地,高地的阻挡作用导致物质无法继续向西北运动,而北边是稳定的中国大陆,因此也没法向北运动,所以梨山断层以西的雪山山脉和台北盆地等地区速度场都比较小.梨山断层以东的中央山脉北部和宜兰平原地区,速度场方向发生了一个巨大的转变,方位角 从320°顺时针转向正北,再继续从正北旋转到北东、 东和东南,在宜兰平原附近达到~140°,方位角顺时针旋转了~180°,速度大小则先是从 ~15 mm/yr向北衰减为0 mm/yr,然后再往东增加到 16.86 mm/yr,在宜兰平原以东的海洋中达到34.71 mm/yr. 这一地区速度场方向的巨大转变则明显是受冲绳海槽的弧后扩张的影响,该地区西边是观音高地,而北边是稳定的中国大陆,所以物质有向东侧逃逸的趋势,加上宜兰平原以东的琉球沟弧 向南后撤,于是使得这一区域的物质转而向东南运动.

台湾西部海岸平原和台湾海峡大部分处于变形前缘以西,靠近澎湖群岛,属于中国大陆边缘,因此相对白沙站点运动的速度都比较小,方向基本为北西向,与菲律宾海板块向欧亚大陆板块运动的方向一致.

从绿岛、兰屿和海岸山脉跨过断层Ⅰ,速度场有很大变化.如图 5所示,红色和黄色箭头分别表示的是断层Ⅰ上西侧和东侧对应点处的速度矢量,黑色箭头表示的是断层Ⅰ上的位错(两侧位移之差).由图中可以看出:跨过断层Ⅰ速度矢量大小和方向均发生了变化,变化量由北向南增加,在断层面南部比较突出,表现为由东向西,速度矢量逆时针旋转,方位角减小了10°~20°,速度大小减小了9.67~17.96 mm/yr,可见断层Ⅰ对台湾地区地壳变形场有着重要的影响.计算得到的断层上的滑移速率为13.81~23.48 mm/yr,由北向南递增,方向与断层 Ⅰ走向一致,与Bos等(2003)计算得到的断层上的位错7.8~25 mm/yr是比较接近的.

4.2 台湾地区应变率场

图 6a是主应变率分布图.图中显示,主压应变最大的地区位于海岸山脉和台湾东部海域,最大主压应变优势方位呈北西-南东向;其方向与菲律宾海板块与欧亚板块的碰撞方向一致,与台湾岛倾斜相交.在台湾的东北部的宜兰平原,主应变率中主张应变值大于主压应变,说明该区域处在拉张的应力环境,这与冲绳海槽的张开密切相关.同样,在台湾南部地区,也出现了拉张的应力环境.具体地质环境,上面已经作了分析,此处不再赘述.

图 6 (a)模型计算的最大主应变率与震源机制结果的对比,(b)模型主应力计算结果与世界地应力结果对比.图中黑色十字箭头为计算所得主应力,箭头向外表示张应变,向内则为压应变,红色线段为世界地应力数据 Fig. 6 (a)Comparison between computed maximum principal strain rate and focal mechanism.(b)Comparison between computed maximum principal strain rate and Word Stress Map. Black cross arrows are computed principal stress,outward arrows indicate compressive strain and inward arrows indicate tensile strain. Red lines are World Stress data

为验证计算所得的主应变率结果的合理性,文中将计算的主应变方位与世界地应力图以及震源机制解给出的P、T轴方位结果进行比较.

文中选取了发生在1976年1月1日-2014年12月31日间,震源深度小于40 km,并且震级为 4≤MW <10的地震(数据来源于http://www.globalcmt.org/CMTsearch.html,最后访问时间2015年1月1日)(Dziewonski et al.,1981; Ekstrm et al.,2012).由于地震事件的空间不均匀性,许多地震震中在空间重叠,无法辨别.因此,为便于观察比较,将研究区域划分为20×26个网格,若在某个小网格里发生了2个以上的地震,则将这些地震的地震矩进行平均,再计算其震源机制,然后在其网格中心绘出其结果.图 6a展示了台湾地区平均后的震源机制解的结果.图中显示,计算给出的台湾地区主应变场与地震的震源机制比较吻合,如在海岸山脉和西部山麓地区,主压应变率大于主张应变率,应变状态属于挤压,此时震源机制也显示这些地区的地震多为逆冲型地震;在中央山脉南部、北部和宜兰平原及附近地区,主张应变率大于主压应变率,表示该地区处于拉张型应变状态,而震源机制解结果也表明在这几个地区,正断层型地震居多(图 6a).在中央山脉北部及南部,主张应变大于主压应变,说明该地区的变形性质与海岸山脉的挤压不同,属于引张变形,这与Hu等(1996)Bos等(2003)朱守彪等(2011)根据GPS数据计算得到的结果是一致的.这可能是因为菲律宾海板块向西北推挤碰撞欧亚板块边缘,导致台湾地区物质向西南和东北两个方向侧向挤出,而马尼拉增生楔的隆起引张作用可能也对该地区的变形场产生影响.在西部山麓和西部海岸平原地区,主压应变大于主张应变占据主导地位,越往西北,主压应变率越小,这可能是因为这些地区离碰撞中心较远,应变衰减了的缘故;实际上该地区地震活动也少.台湾西南部的屏东平原地区主压应变比主张应力大,呈挤压状态,与Hu等(1996)Bos等(2003)朱守彪等(2011)得出的结论相似;同时震源机制解结果也显示,屏东平原地区多为逆断层地震.图中还显示,该地区的主压应变轴发生了逆时针旋转,这可能与台湾西南部地区物质向西南逃逸,马尼拉海沟向西扩张和西边北康高地的阻挡三者有关.在东北部宜兰平原,变形状态比较复杂,总体来说处于拉张的状态,但是其南部又处于挤压状态,这与震源机制解结果表现也是一致的.这种应变场特征可能是由台湾地区物质向东北方向逃逸和冲绳海槽扩展以及琉球沟弧系统向南后撤造成的.

此外,下面还将计算结果与世界地应力图(World Stress Map,简称WSM)给出的台湾地区应力测量结果进行比较(图 6b). WSM是全球当前地壳应力场图,地应力数据有两种形式,一种是原始数据,一种是平滑过的数据,由于原始数据空间分布不均匀,各地区数据密度也相差较大,为了便于分析比较,我们选取并下载了台湾地区的平滑过的地应力数据(数据来源于http://dc-app3-14.gfz-potsdam. de/pub/introduction/introduction_frame.html,http: //dc-app3-14.gfz-potsdam.de/pub/stress_data/stress_data_frame.htm,最后访问时间2015年1月1日)(Heidbach et al.,2010).从图 6b可以看出,除了台湾北部局部地区之外,在整个研究区域,计算所得的主应力方位与WSM结果比较吻合,说明模型的计算结果是合理的.

通过将计算结果与震源机制结果及WSM数据比较,可知本文所采用的模型基本上可以勾画台湾地区复杂的地壳变形情况.

图 7为计算得到的台湾地区最大剪应变率场分布图,由于我们主要关心台湾岛上的情况,因此只给出了岛内的应变场.从图中可以看出,总体上台湾南部的最大剪应力比北部的大.最大值出现在台湾西南部的屏东平原地区,次大值分布在台湾东部沿海从花莲到台东即海岸山脉地区.这与Bos等(2003)朱守彪等(2011)计算的结果有一些不同,可能是由于对模型中西南方向的边界处理不同造成的.台湾北部和西部地区的最大剪应变率比较小,这可能是因为菲律宾海板块与欧亚板块碰撞产生的变形一部分被海岸山脉和纵谷断裂所吸收,因此海岸山脉地区的最大剪应变率是比较大的.

图 7 计算得到的台湾地区最大剪应变率场 Fig. 7 Computed maximum shear strain rates field
5 问题讨论 5.1 计算的速度场与GPS观测结果的对比

为考察计算结果的合理性,下面将计算的速度矢量与GPS观测结果进行对比.图 8是计算得到的在GPS测点处的速度矢量与原GPS观测结果的比较,图中红色和蓝色箭头分别为观测和计算结果.需要指出的是,图中GPS的观测数据来自Ching等(2011).尽管台湾及周边地区在GPS观测时段里发生了很多地震,如1999年的集集地震等,但Ching等(2011)删除了受大震同震影响的GPS数据,因此采用的GPS数据反映的是震间地壳运动.

图 8 模型中计算的速度矢量与GPS观测结果的对 比. 图中红色箭头表示GPS观测结果,蓝色箭头表示在GPS测点处计算的速度矢量 Fig. 8 Comparison between computed velocity vectors and GPS observation results. Computed velocity vectors as blue arrows and GPS observation results as red arrows

图 8显示,总体来说模型计算的速度矢量与观测结果比较符合.由于菲律宾海板块向亚欧大陆板块挤压的方位角是306°,即N54°W,而纵谷断层的走向大体为N20°E,这之间有~74°的夹角,所以存在一个平行于纵谷断层的剪切应力分量,也就导致在纵谷断层两侧产生剪切变形.因此,模型中纵谷断层两侧位移计算结果的方向明显由北-西向转向北西西-南东东方向,这一点在断层Ⅰ的南部表现得最为明显(上文已作详细讨论).断层Ⅰ的存在,正好解释了为什么由东向西跨过纵谷断层,速度矢量发生变化.台湾西南部屏东平原,速度矢量发生了逆时针旋转,与GPS观测表现一致;在中央山脉北部和台北盆地,计算结果展示了速度发生顺时针旋转,与观测值具有较高的一致性;并且这一区域速度矢量很小,与GPS观测也一致;宜兰平原和台湾东北部沿海,计算所得速度矢量的方向由东-西向转向南-东向,方向变化趋势与观测结果比较吻合,速度大小也与观测较为符合.

由于纵谷断层不是一个纯左旋走滑断层,而是 左旋兼逆冲断层(Yu and Liu,1989; Lee and Angelier,1993b),而且逆冲成分很大.纵谷断层北部,逆冲成分比较小,断层两侧相对运动主要是左旋走滑,速度 为6~16 mm/yr;中部断层两侧的相对运动左旋走滑部分为19~24 mm/yr,而逆冲部分为16~22 mm/yr,接近1 ∶ 1;南部断层两侧的相对运动左旋走滑和逆冲部分均为29~34 mm/yr(Yu et al.,1997).而本文模型为二维平面应力模型,无法模拟地表垂直运动,尽管模型中的断层具有一定的宽度并进行了处理(其内充填软弱物质),但是还是没法模拟出断层的逆冲运动特征.因此,在纵谷断层的中部和北部即美仑断层附近,计算结果与观测值比较符合,但在纵谷断层的南部逆冲分量比较高的地区,计算结果与观测值偏差较大.另外宜兰平原南部计算所得速度比观测值偏小,则可能是由该地区东边GPS数据较少,难以进行很好的约束造成的.

总之,模型计算的台湾地区震间的速度矢量与GPS实际观测结果有较好的一致性,由此充分说明文中的有限元模型是合理的,计算结果是可靠的.

5.2 摩擦系数对速度场的影响

由ANSYS优化结果可知除断层Ⅰ外,其他断 层上的摩擦系数均很大,对台湾地区变形场的影响比较小(Bos et al.,2003).为了抓住主要的因素进行分析,下面只对断层Ⅰ上的摩擦系数对台湾地区变形场的影响进行定量分析.上文中ANSYS优化给出的断层Ⅰ上的摩擦系数为μ=0.5,为了研究摩擦系数对变形场的影响,本文还建立了另外两个模型,分别取断层Ⅰ上摩擦系数μ1=0.1和μ2=1.0(其他断层上的摩擦系数均保持不变),然后把计算结果与μ=0.5的模型结果进行对比.图 9a和9b分别是模型在摩擦系数为μ1=0.1和μ2=1.0时的速度场计算结果.比较图 5图 8图 9a9b可以看出,三种情况下速度场整体变化都不大,但是断层Ⅰ两侧速度矢量对比随断层摩擦系数变化却很明显.当断层面上摩擦系数为0.1时(图 9a),断层两侧速度差明显增大,相反增大断层面上的摩擦系数(图 9b),断层两侧速度差异变小.这一结果说明台湾地区速度场对断层上的摩擦系数变化是敏感的,也说明断层Ⅰ是台湾地区位移的不连续面,断层Ⅰ的存在对台湾地区地壳变形有着重要的影响,这与Yu等(1990)Yu等(1999)Hu等(2001)Bos等(2003)的观点都是一致的.比较图 9a图 5可以看出,图 9a中计算值与观测值有较大偏离.虽然图 9b位移 场看起来与图 5相差不大,但是图 9b断层两侧速度相差太小,偏离Yu等(1990)Yu等(1999)的结论30 mm/yr较多(Yu et al.,19901999),断层上滑移速率较小,与Bos等(2003)计算得到的断层上的位错7.8~25 mm/yr相差较大(Bos et al.,2003),因此本文没有采纳μ2=1.0的模型.

图 9 断层面摩擦系数不同时的速度场计算结果和断层Ⅰ两侧速度矢量对比.(a)μ1=0.1,(b)μ2=1.0.插图中黑色和黄色箭头分别表示的是断层Ⅰ上同一点处西侧和东侧的速度矢量 Fig. 9 Computed velocity field results and the velocity comparison between two sides of fault group Ⅰ with different fault Friction coefficients.(a)Friction coefficient of fault group Ⅰ,μ1=0.1.(b)Friction coefficient of fault group Ⅰ,μ2=1.0. The black and yellow arrows in the insert indicate the velocities at the corresponding points on both sides of fault group Ⅰ

由于本文采用的二维平面应力模型是对实际三维情况的简化,而且断层上不同区段的摩擦系数也有可能会不同,所以本文给出的断层上的摩擦系数 0.5并不一定能准确反映相应断层上的真实摩擦系数.

由于台湾地区特殊的地质构造,变形场十分复杂,文中利用二维有限元模型来进行计算,尽管对震间的地壳变形、地面运动速度等进行了模拟,取得了一定的认识,但其结果只能说是一级近似.要得到符合实际地质构造的理想结果,还需考虑更多因素,进行三维的有限元模拟.

6 结论

利用台湾地区GPS观测结果作为有限元模拟的边界约束,计算获得了台湾地区地壳运动的速度场、变形场以及应力场,同时给出了影响台湾地区变形的主要因素,现总结如下.

本文通过计算得到的台湾地区的速度场与GPS观测结果基本吻合,计算给出的主应力与WSM结果一致性很好,以及计算的应力结果与震源机制有可比性.这些都说明文中所建立的模型是 合理的.因此,模拟给出的台湾地区变形场是可靠的.

计算所得的台湾地区速度场方向由北向南基本呈扇形分布,总体来说速度场大小从东南往西北逐渐减小,中东部速度场方向为西北向,往北速度场方向顺时针旋转,往南则逆时针旋转.这可能是由于菲律宾海板块向西北方向挤压台湾,而台湾西边是凸出的北康高地和观音高地,向西的运动受到强硬高地的阻挡,造成台湾地区物质向东北和西南两个方 向侧向挤出造成的.跨过纵谷断层Ⅰ,速度场大小有较大的衰减,纵谷断裂滑移速率约为13.8~23.5 mm/yr,说明一部分的变形被纵谷断裂所吸收.

主压应变率最大的地区是海岸山脉和台湾东部海域,主应变率场呈北西-南东向分布,方向与菲律宾海板块碰撞台湾的方向一致.由海岸山脉跨过纵谷断层往西和西北,主压应变率逐渐减小.中央山脉北部及南部地区,主张应变大于主压应变,相应的该地区的震源机制解也表现为正断层地震.这可能是因为菲律宾海板块向西北推挤碰撞欧亚板块边缘,导致台湾地区物质向西南和东北两个方向侧向挤出造成的.在西部山麓和西部海岸平原地区,主压应变大于主张应变,越往西北,主压应变率越小,这可能是因为这些地区远离碰撞中心、应变衰减了的缘故.台湾西南部的屏东平原地区呈挤压状态,震源机制解结果也显示屏东平原地区多逆冲型地震.该地区的主压应变轴发生逆时针方向旋转,这可能与台湾西南部地区物质向西南逃逸,马尼拉海沟向西扩张和西边北康高地的阻挡有关.在台湾东北部宜兰平原,变形状态也较复杂,总体来说处于拉张的应力状态,但是其南部又处于挤压状态,这与震源机制解结果表现也是一致的.这种应变场特征可能是由受碰撞挤压物质向东北方向逃逸以及冲绳海槽的扩展和琉球海沟向南后撤造成的.

从文中的计算和分析还可以看出,台湾地区的变形场受多种因素控制.其中菲律宾海板块向北西方向的推挤,台湾西部和西北部稳定的北康、观音高地的抵阻,台湾东北部冲绳海槽的扩张及琉球海沟的向南后撤,以及台东纵谷断裂的无震错动等诸多影响,共同造成了台湾地区现今应变场的主要格局.

致谢 十分感谢审稿专家提出的宝贵意见.

参考文献
[1] Angelier J, Chu H T, Lee J C, et al. 2000. Active faulting and earthquake hazard: The case study of the Chihshang Fault, Taiwan. Journal of Geodynamics, 29(3-5): 151-185.
[2] Angelier J, Chang T Y, Hu J C, et al. 2009. Does extrusion occur at both tips of the Taiwan collision belt? Insights from active deformation studies in the Ilan Plain and Pingtung Plain regions. Tectonophysics, 466(3-4): 356-376.
[3] Barrier E, Angelier J. 1986. Active collision in eastern Taiwan: the Coastal Range. Tectonophysics, 125(1-3): 39-72.
[4] Biq C C. 1972. Transcurrent buckling, transform faulting and transpression, their relevance in eastern Taiwan kinematics. Petrol. Geol. Taiwan, 20: 1-39.
[5] Bonilla M G. 1975. A review of recently active faults in Taiwan. US: US Geological Survey.
[6] Bonilla M G. 1977. Summary of Quaternary faulting and elevation changes in Taiwan. Mem. Geol. Soc. China, 2: 43-55.
[7] Bos A G, Spakman W, Nyst M C J. 2003. Surface deformation and tectonic setting of Taiwan inferred from a GPS velocity field. Journal of Geophysical Research: Solid Earth(1978—2012), 108(B10):2458.
[8] Chang C H, Wu Y M, Shin T C, et al. 2000. Relocation of the 1999 Chi-Chi earthquake in Taiwan. Terrestrial Atmospheric and Oceanic Sciences, 11(3): 581-590.
[9] Chang C P, Chang T Y, Angelier J, et al. 2003. Strain and stress field in Taiwan oblique convergent system: constraints from GPS observation and tectonic data. Earth and Planetary Science Letters, 214(1-2): 115-127.
[10] Chen K P, Tsai Y B. 2008. A catalog of Taiwan earthquakes (1900—2006) with homogenized MW magnitudes. Bulletin of the Seismological Society of America, 98(1): 483-489.
[11] Cheng S N. 1995. The study of stress distribution in and around Taiwan[Ph. D. thesis](in Chinese). Taiwan: National Central University.
[12] Chi W R, Namson J, Suppe J. 1981. Stratigraphic record of plate interactions in the Coastal Range of eastern Taiwan. Mem. Geol. Soc.China, 4: 155-194.
[13] Chiang C W, Chen C C, Chen C S, et al. 2008. Deep electrical structure of the southern Taiwan orogeny and its tectonic implications by MT data. // Abstracts of the 19th International Workshop on Electromagnetic Induction in the Earth (Volume 1).
[14] Ching K E, Rau R J, Lee J C, et al. 2007. Contemporary deformation of tectonic escape in SW Taiwan from GPS observations, 1995—2005. Earth and Planetary Science Letters, 262(3-4): 601-619.
[15] Ching K E, Rau R J, Johnson K M, et al. 2011. Present-day kinematics of active mountain building in Taiwan from GPS observations during 1995—2005. Journal of Geophysical Research: Solid Earth, 116(B9): B09405, doi: 10.1029/2010jb008058.
[16] Clift P D, Schouten H, Draut A E. 2003. A general model of arc-continent collision and subduction polarity reversal from Taiwan and the Irish Caledonides. Geological Society, London, Special Publications, 219(1): 81-98.
[17] Dai L M, Li S Z, Lou D, et al. 2013. Numerical modeling on the stress filed in the Huanghua depression, Bohai Bay basin. Chinese J. Geophys. (in Chinese), 56(3): 929-942, doi: 10.6038/cjg20130321.
[18] 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. Journal of Geophysical Research: Solid Earth (1978—2012), 86(B4): 2825-2852.
[19] Ekström G, Nettles M, Dziewoński A M. 2012. The global CMT project 2004—2010: centroid-moment tensors for 13017 earthquakes. Physics of the Earth and Planetary Interiors, 200-201: 1-9.
[20] Hao J L, Wang W M, Wang J, et al. 2009. A dislocation model of elastic block for aseismic crustal deformation in Taiwan. Chinese J. Geophys. (in Chinese), 52(5): 1223-1232, DOI: 10.3969/j.issn.0001-5733.2009.05.011.
[21] He Y M, Yao Z X. 2002. Dislocation model for current crustal deformation in and around southern Taiwan, China. Chinese J. Geophys.(in Chinese), 45(5): 638-645.
[22] Heidbach O, Tingay M, Barth A, et al. 2010. Global crustal stress pattern based on the World Stress Map database release 2008. Tectonophysics, 482(1-4): 3-15.
[23] Ho C S. 1986. A synthesis of the geologic evolution of Taiwan. Tectonophysics, 125(1-3): 1-16.
[24] Ho C S. 1988. An introduction to the geology of Taiwan: explanatory text of the geologic map of Taiwan. Taipei, Taiwan: Central Geological Survey, Ministry of Economic Affairs.
[25] Hsu M T. 1980. Earthquake Catalogues in Taiwan (from 1644 to 1979)[Ph. D. thesis]. Taipei: National Taiwan University, 77.
[26] Hsu T L. 1976. Neotectonics of the Longitudinal Valley, eastern Taiwan. Bull. Geol. Surv. Taiwan, 25: 53-62.
[27] Hsu V. 1990. Seismicity and tectonics of a continent-island arc collision zone at the Island of Taiwan. Journal of Geophysical Research: Solid Earth (1978—2012), 95(B4): 4725-4734.
[28] Hsu Y J, Yu S B, Simons M, et al. 2009. Interseismic crustal deformation in the Taiwan plate boundary zone revealed by GPS observations, seismicity, and earthquake focal mechanisms. Tectonophysics, 479(1-2): 4-18.
[29] Hu J C, Angelier J, Lee J C, et al. 1996. Kinematics of convergence, deformation and stress distribution in the Taiwan collision area: 2-D finite-element numerical modelling. Tectonophysics, 255(3-4): 243-268.
[30] Hu J C, Angelier J, Yu S B. 1997. An interpretation of the active deformation of southern Taiwan based on numerical simulation and GPS studies. Tectonophysics, 274(1): 145-169.
[31] Hu J C, Yu S B, Angelier J, et al. 2001. Active deformation of Taiwan from GPS measurements and numerical simulations. Journal of Geophysical Research: Solid Earth (1978—2012), 106(B2): 2265-2280.
[32] Huang M H, Hu J C, Hsieh C S, et al. 2006. A growing structure near the deformation front in SW Taiwan as deduced from SAR interferometry and geodetic observation. Geophysical Research Letters, 33(12): L12305, doi: 10.1029/2005gl025613.
[33] Huang W J, Johnson K M, Fukuda J, et al. 2010. Insights into active tectonics of eastern Taiwan from analyses of geodetic and geologic data. Journal of Geophysical Research: Solid Earth (1978—2012), 115(B3): B03413.
[34] Huchon P, Barrier E, de Bremaecker J C, et al. 1986. Collision and stress trajectories in Taiwan: a finite element model. Tectonophysics, 125(1-3): 179-191.
[35] Kao H, Wu F T. 1996. The 16 September 1994 earthquake (mb=6.5) in the Taiwan strait and its tectonic implications. Terrestrial, Atmosphere and Oceanic Science, 7(1): 13-29.
[36] Kao H, Shen S S J, Ma K F. 1998. Transition from oblique subduction to collision: Earthquakes in the southernmost Ryukyu arc-Taiwan region. Journal of Geophysical Research: Solid Earth (1978—2012), 103(B4): 7211-7229.
[37] Kao H, Chen W P. 2000. The Chi-Chi earthquake sequence: active, out-of-sequence thrust faulting in Taiwan. Science, 288(5475): 2346-2349.
[38] Kuo-Chen H, Wu F T, Jenkins D M, et al. 2012. Seismic evidence for the α-β quartz transition beneath Taiwan from Vp/Vs tomography. Geophysical Research Letters, 39(22): L22302.
[39] Lallemand S, Font Y, Bijwaard H, et al. 2001. New insights on 3-D plates interaction near Taiwan from tomography and tectonic implications. Tectonophysics, 335(3-4): 229-253.
[40] Lee J C, Angelier J. 1993a. Localisation des déformations actives et traitement des données géodésiques: l'exemple de la faille de la Vallée Longitudinale, Tawan. Bull. Soc. Géol. France, 164(4): 533-570.
[41] Lee J C, Angelier J. 1993b. Location of active deformations and geodetic data analyses: An example of the longitudinal valley fault, taiwan. Bull. Soc. Géol. France, 164(4): 533-540.
[42] Lee J C, Angelier J, Chu H T, et al. 2001. Continuous monitoring of an active fault in a plate suture zone: a creepmeter study of the Chihshang Fault, eastern Taiwan. Tectonophysics, 333(1-2): 219-240.
[43] Lee J C, Angelier J, Chu H T, et al. 2003. Active fault creep variations at Chihshang, Taiwan, revealed by creep meter monitoring, 1998—2001. Journal of Geophysical Research: Solid Earth (1978—2012), 108(B11): 2528.
[44] Lin A T, Watts A B. 2002. Origin of the West Taiwan basin by orogenic loading and flexure of a rifted continental margin. Journal of Geophysical Research: Solid Earth (1978—2012), 107(B9): 2185.
[45] Lin C H, Yeh Y H, Yen H Y, et al. 1998. Three-dimensional elastic wave velocity structure of the Hualien region of Taiwan: Evidence of active crustal exhumation. Tectonics, 17(1): 89-103.
[46] Lin K C, Hu J C, Ching K E, et al. 2010. GPS crustal deformation, strain rate, and seismic activity after the 1999 Chi-Chi earthquake in Taiwan. Journal of Geophysical Research: Solid Earth, 115(B7): B07404, doi: 10.1029/2009JB0064-17.
[47] Lu C Y, Jeng F S, Chang K J, et al. 1998. Impact of basement high on the structure and kinematics of the western Taiwan thrust wedge: insights from sandbox models. Terrestrial, Atmospheric and Oceanic Sciences, 9(3): 533-550.
[48] Lu C Y, Malavieille J. 1994. Oblique convergence, indentation and rotation tectonics in the Taiwan Mountain Belt: Insights from experimental modelling. Earth and Planetary Science Letters, 121(3-4): 477-494.
[49] Pezzopane S K, Wesnousky S G. 1989. Large earthquakes and crustal deformation near Taiwan. Journal of Geophysical Research: Solid Earth (1978—2012), 94(B6): 7250-7264.
[50] Rau R J, Wu F T. 1995. Tomographic imaging of lithospheric structures under Taiwan. Earth and Planetary Science Letters, 133(3-4): 517-532.
[51] Rau R J, Ching K E, Hu J C, et al. 2008. Crustal deformation and block kinematics in transition from collision to subduction: Global positioning system measurements in northern Taiwan, 1995—2005. Journal of Geophysical Research: Solid Earth (1978—2012), 113(B9): B09404.
[52] Shi Y L, Zhu S B. 2006. Discussion on method of calculating strain with GPS displacement data. Journal of Geodesy and Geodynamics (in Chinese), 26(1): 1-8.
[53] Sun Y J, Zhang H, Shi Y L, et al. 2011. Numerical investigation on the geodynamical mechanism of the first major shock of 2006 Pingtung MW7.0 earthquake. Sci. China Earth Sci. (in Chinese), 54(5): 631-639, doi: 10.1007/s11430-010-4114-9.
[54] Suppe J. 1984. Kinematics of arc-continent collision, flipping of subduction, and back-arc spreading near Taiwan. Mem. Geol. Soc.China, 6: 21-33.
[55] Teng L S. 1987. Stratigraphic records of the late Cenozoic Penglai orogeny of Taiwan. Acta Geologica Taiwanica, (25): 205-224.
[56] Teng L S. 1990. Geotectonic evolution of late Cenozoic arc-continent collision in Taiwan. Tectonophysics, 183(1-4): 57-76, doi: http://dx.doi.org/10.1016/0040-1951(90)90188-E.
[57] Teng L S. 1996. Extensional collapse of the northern Taiwan mountain belt. Geology, 24(10): 949-952.
[58] Tsai Y B. 1986. Seismotectonics of Taiwan. Tectonophysics, 125(1-3): 17-28, 31-37.
[59] Wu F T. 1978. Recent tectonics of Taiwan. Journal of Physics of the Earth, 26(Suppl): S265-S299.
[60] Wu Y M, Chang C H, Zhao L, et al. 2007. Seismic tomography of Taiwan: Improved constraints from a dense network of strong motion stations. Journal of Geophysical Research: Solid Earth (1978—2012), 112(B8): B08312.
[61] Wu Y M, Zhao L, Chang C H, et al. 2009. Relocation of the 2006 Pingtung Earthquake sequence and seismotectonics in Southern Taiwan. Tectonophysics, 479(1-2): 19-27.
[62] Yu S B, Liu C C. 1989. Fault creep on the central segment of the Longitudinal Valley fault, eastern Taiwan. Proc. Geol. Soc. China, 32(3): 209-231.
[63] Yu S B, Jackson D D, Yu G K, et al. 1990. Dislocation model for crustal deformation in the Longitudinal Valley area, Eastern Taiwan. Tectonophysics, 183(1-4): 97-109.
[64] Yu S B, Chen H Y. 1994. Global positioning system measurements of crustal deformation in the Taiwan arc-continent collision zone. Terrestrial, Atmospheric and Oceanic Sciences, 5(4): 477-498.
[65] Yu S B, Chen H Y, Kuo L C. 1997. Velocity field of GPS stations in the Taiwan area. Tectonophysics, 274(1): 41-59.
[66] Yu S B, Kuo L C, Punongbayan R S, et al. 1999. GPS observation of crustal deformation in the Taiwan-Luzon Region. Geophysical Research Letters, 26(7): 923-926, doi: 10.1029/1999gl900148.
[67] Zhang C H. 2008. ANSYS 11.0 Engineering Application Examples of Analytic Structure (in Chinese). Beijing: China Machine Press.
[68] Zhu S B, Shi Y L. 2004. Genetic algorithm-finite element inversion of drag forces exerted by the lower crust on the upper crust in the Sichuan-Yunnan area. Chinese J. Geophys. (in Chinese), 47(2): 232-239.
[69] Zhu S B, Shi Y L. 2005. Genetic algorithm-finite element inversion of topographic spreading forces and drag forces of lower crust to upper crust in Tibetan Plateau. Acta Scientiarum Naturalium Universitatis Pekinensis (in Chinese), 41(2): 225-234.
[70] Zhu S B, Cai Y E. 2006. Inversion of viscous properties of crust and mantle from the GPS temporal series measurements. Chinese J. Geophys.(in Chinese), 49(3): 771-777.
[71] Zhu S B, Shi Y L. 2007. Origin of tectonic stresses in the Chinese continent and adjacent areas. Science in China Series D: Earth Sciences (in Chinese), 50(1): 67-74.
[72] Zhu S B, Xing H L, Xie F R, et al. 2008. Simulation of earthquake processes by finite element method: The case of megathrust earthquakes on the Sumatra subduction zone. Chinese J. Geophys. (in Chinese), 51(2): 460-468.
[73] Zhu S B, Cai Y E. 2009. Dynamic mechanisms of the post-seismic deformation following large events: Case study of the 1999 Chi-Chi earthquake in Taiwan of China. Science in China Series D: Earth Sciences (in Chinese), 52(11): 1813-1824.
[74] Zhu S B, Zhang P Z. 2010. Numeric modeling of the strain accumulation and release of the 2008 Wenchuan, Sichuan, China, earthquake. Bull. Seismol. Soc. Am., 100(5B): 825-2839.
[75] Zhu S B, Zhao X Y, Liu Y, et al. 2011. Distribution of strain rates in Taiwan before and after 1999 Chi-Chi earthquake and their geodynamic mechanisms. Journal of Geodesy and Geodynamics (in Chinese), 31(2): 38-43.
[76] Zhu S B, Zhang P Z. 2013. FEM simulation of interseismic and coseismic deformation associated with the 2008 Wenchuan earthquake. Tectonophysics, 584: 64-80, doi: 10.1016/j.tecto.2012.06.024.
[77] 郑世楠. 1995. 台湾及其邻近地区大地应力分布的研究 [博士论文]. 台湾: 国立“中央”大学地球物理研究所.
[78] 戴黎明, 李三忠, 楼达等. 2013. 渤海湾盆地黄骅坳陷应力场的三维数值模拟分析. 地球物理学报, 56(3): 929-942.
[79] 郝金来, 王卫民, 王建等. 2009. 台湾地区地壳形变的弹性块体位错模型. 地球物理学报, 52(5): 1223-1232.
[80] 何玉梅, 姚振兴. 2002. 中国台湾南部及其周边岛屿现今地壳形变的位错模型. 地球物理学报, 45(5): 638-645.
[81] 石耀霖, 朱守彪. 2006. 用GPS位移资料计算应变方法的讨论. 大地测量与地球动力学, 26(1): 1-8.
[82] 孙玉军, 张怀, 石耀霖. 2011. 2006年12月26日台湾屏东地震机制的地球动力学成因探讨. 中国科学: 地球科学, 40(10): 1301-1309.
[83] 张朝晖. 2008. ANSYS 11.0结构分析工程应用实例解析. 北京: 机械工业出版社.
[84] 朱守彪, 石耀霖. 2004. 用遗传有限单元法反演川滇下地壳流动对上地壳的拖曳作用. 地球物理学报, 47(2): 232-239, doi: 10.3321/j.issn:0001-5733.2004.02.009.
[85] 朱守彪, 石耀霖. 2005. 青藏高原地形扩展力以及下地壳对上地壳的拖曳力的遗传有限单元法反演. 北京大学学报(自然科学版), 41(2): 225-234.
[86] 朱守彪, 蔡永恩. 2006. 利用GPS观测的时间序列资料反演地壳地幔黏性结构. 地球物理学报, 49(3): 771-777.
[87] 朱守彪, 石耀霖. 2007. 中国大陆及邻区构造应力场成因的研究. 中国科学D 辑 地球科学, 36(12): 1077-1083.
[88] 朱守彪, 邢会林, 谢富仁等. 2008. 地震发生过程的有限单元法模拟——以苏门答腊俯冲带上的大地震为例. 地球物理学报, 51(2): 460-468, doi: 10.3321/j.issn:0001-5733.2008.02.018.
[89] 朱守彪, 蔡永恩. 2009. 强震后地表变形的动力学机制研究—以1999年台湾集集地震为例. 中国科学 D辑: 地球科学, 39(9): 1209-1219.
[90] 朱守彪, 赵晓燕, 刘杨等. 2011. 1999 年集集地震前后台湾地区地应变率场的分布及其动力学成因. 大地测量与地球动力学, 31(2): 38-43.