2010年2月27日，智利Maule地区近海发生M_W8.8级特大地震，震源深度35 km, 震中位置 (35.85°S, 72.72°W) 位于康塞普西翁 (Concepcion) 东北89 km，距智利首都圣地亚哥西南339 km，此次地震发生在纳兹卡板块与南美板块的边界部.纳兹卡板块俯冲到南美板块之下，相对挤压速率约77 mm·a^-1(Norabuena et al., 1998) 造成了安第斯山脉的隆起、火山的活动以及频繁地震的发生 (图 1).该俯冲带具有明显沿纬度方向火山、地震活动分节现象 (Chinn et al., 1980; 马宗晋, 1983)，Barazangi和Isacks (1976)将该现象归结于南美板块之下的大洋板块具有不同俯冲角度，俯冲角度近水平的板块和上覆板块间只有很少的软流圈物质，无火山喷发物质基础，俯冲角度较陡板块上方则分布着密集火山.本次Maule地震俯冲角度约为18°，破裂从震中开始向南、北两个方向双侧扩展，北侧破裂明显强于南侧，破裂长度约650 km (Tong et al., 2010).20世纪以来, 南美大陆许多7.5级以上地震集中在纳兹卡板块和南美洲板块之间俯冲带上 (Bilek, 2010)，其中最大一次为1960年M_W9.5级地震，数次强震造成的破裂几乎贯穿智利近海两千公里.在秘鲁大陆边缘发生了1960年和1996年“海啸地震”.在该俯冲带还存在一些慢速地震过程，比如1995年M_W8.1 Antofagasta地震的震后余滑 (Pritchard and Simons, 2006).

图 1

Fig. 1

图 1 区域构造图 黑色震源球数据来源美国哈佛大学Global CMT (http://www.globalcmt.org[2015-10-20])；绿色震源球为2010年Maule地震；蓝色正方形表示历史地震，数据来源于 (Bilek, 2010)，红色三角形代表第四季火山活动，数据来源 (http://www.ngdc.noaa.gov/hazard/volcano.shtml[2015-10-20]). Fig. 1 Regional tectonic map Black beach balls are the centroid moment tensor; Green beach ball indicates Maule earthquake; Blue squares are historical earthquakes (Bilek, 2010); Red triangles are Quaternary volcano.

不同观测手段记录到了此次智利特大地震造成的影响.GPS观测结果 (尹志强等, 2011) 表明震中附近产生明显位移，康塞普西翁市整体向西移动达3 m, 智利首都圣地亚哥向西移动约30 cm.地震引发海啸，其中海啸受灾人数为124人 (Fritz et al., 2011)，台湾周边海域的验潮站记录到明显的海啸信号，2月28日海啸最大波高达44 cm (谢燕双等, 2012).Peng等 (2010)在加州中部Coso地热区发现了此次智利地震远程触发的微地震和Tremor，在新西兰也发现了智利地震触发的高频Tremor信号 (Fry et al., 2011).该地震是人类有历史记载以来的第六大地震，那么除了观测手段之外，能否通过计算定量了解这样一个特大地震造成的同震影响，这是本文的研究目的.

目前广泛使用的计算地震同震形变、应力方法中，具有代表性的是Okada (1985, 1992) 弹性半无限空间位错表达式和Wang等 (2003)的EDGRN/EDCMP程序，但它们均采用平面坐标系.石耀霖和朱守彪 (2006)指出在由位移计算应变时，必须使用球坐标，特别当位移量较大或研究区面积大或在高纬度地区时，采用直角坐标带来的误差量级不可忽略.Sun和Okubo (2002)的研究也表明地球曲率对同震垂直位移的影响不可忽略.球型分层弹性地球模型的半解析方法 (Sun，1992) 可以考虑地球曲率，但无法考虑真实地球的介质强横向不均性和地表起伏.在同震形变和应力计算中，地表起伏和介质的不均匀性往往不可忽略，例如：Lin等 (2013)发现青藏高原东缘地形对汶川地震同震水平位移在近断层处最高影响为9%；2001年昆仑山M_S8.1地震同震变形场在断裂带两侧的非对称分布必须考虑断层两侧弹性系数的不均匀 (邵志刚等，2008).真实地球横向不均匀性很大，半解析方法处理该类问题具有局限性，例如付广裕和孙文科 (2012)采用的“微小扰动方法”，该方法只适用于很小的三维构造变化的情况.目前未见公开发表文献在同震形变和应力计算中同时考虑球形地球、介质非均匀性和地表、Moho面地形起伏，本文采用数值方法，在同震应变应力计算中考虑椭球型地球、非均匀介质和界面起伏，更接近真实地球，对地震危险性的定量分析具有重要意义.

2010年智利Maule特大地震发生之后，不少学者计算了该地震引起的同震库仑应力变化 (缪淼和朱守彪，2012；Farías et al., 2011；Lorito et al., 2011)，但这些研究在平面直角坐标系下计算应变和应力，也未能考虑安第斯地区剧烈的地表地形起伏和地壳结构的横向不均匀性.本文采用有限元方法，使用含地表地形起伏和横向不均匀地壳结构椭球地球模型，从力学角度定量计算2010年智利Maule特大地震的同震形变场和应力场，讨论地球曲率、分层结构对同震形变场的影响，并根据库仑应力变化分析周围断层的地震活动性和主震对余震的触发关系.

半无限空间模型和球形地球模型计算的同震位移存在较大差异，Dong等 (2014)发现埋深100 km的走滑型点源位错在震中距约200 km处地球曲率对同震位移影响为5%.从位移计算应变以及应力时也需要在球坐标下进行 (石耀霖和朱守彪，2006)，因此使用平面假设将带来位移、应力、应变计算不可忽略的误差，必须选用球形地球模型.Qu等 (2015)计算2004年苏门答腊地震发现，在断层破裂面所处近场、地形及三维横向不均匀地壳结构对同震水平、垂直位移的影响分别为23%和40%，这样的影响不可忽略.Sun和Okubo (2002)的半解析方法虽然可以考虑地球曲率和分层，却无法处理真实地球的介质强不均匀性和界面起伏，且该方法目前尚未计算同震应力变化，因此要研究真实地球，需要同时考虑地球曲率、介质的不均匀性和界面起伏，必须使用数值方法.

地震位错处理采用张贝等 (2015)提出的“有限元模拟弹性位错的等效体力方法”.该方法对位错处理体现在单元积分中将位错等效为体力添加到平衡方程右端项中.若单元包含位错，则在计算单元积分时，等效体力计算表达式为

(1)

其中，f_m^K为等效体力，Φ_m表示单元K中的第m个形函数，ν_j是位错面的法向矢量，[u (r)_i]为位错，C_{ij pq}为弹性系数，dΣ是单元内部位错面.

计算网格包含整个地球，上至地表，下至核幔边界.自然回避侧边界条件不确定性，只需给出上下边界条件，上边界为自由边界 (正应力、剪应力均为0)，核幔边界为弹簧边界条件，其值由核幔密度差值决定.地表地形数据采用ETOPO1.0(Amante and Eakins, 2009)，地壳采用Crust1.0(Laske et al., 2013) 模型，地幔采用GyPSuM模型 (Simmons et al., 2010)，模型中使用的材料参数按照上述模型波速结构计算得到，地球椭球形状根据WGS-84参考椭球确定.采用根据误差识别自动调整网格疏密度的网格自适应技术，能够合理计算费用，提高复杂问题计算效率，改进结果精度.断层滑动模型使用Shao等 (2010)模型，对断层进行了加密处理，采用带悬点六面体单元，最终单元数目419万，计算网格见图 2.

图 2

Fig. 2

图 2 计算网格地形放大30倍显示. Fig. 2 The computational grids The topography is vertically exaggerated by 30 times to highlight the elevation.

大地震造成断层面及邻近区域的库仑应力变化被广泛用于研究地震对周围断层的应力触发作用，探讨大地震与后续地震关系.一般定义库仑应力变化 (Harris, 1998; 石耀霖和曹建玲, 2010) 为

(2)

式中，ΔCFS为库仑应力变化，Δτ为断层面上剪应力变化量 (与滑动方向同向为正)，Δσ_n为断层面上正应力变化 (拉伸为正)，ΔP为孔隙水压变化，μ为断层面摩擦系数.在实际应用中，常把孔隙水压的作用通过视摩擦系数μ′体现：

(3)

本研究中，取视摩擦系数为典型值μ′=0.4.

分别用有限元方法计算球面分层模型同震位移，用EDGRN/EDCMP2.0程序 (Wang et al., 2003) 计算平面分层模型、平面半空间模型同震形变.这三个模型震源模型一样，介质参数都依据PREM模型，结果差别反映地球曲率和分层结构影响.

图 3a为球面分层模型计算同震位移和GPS观测比较，GPS数据来源于Vigny等 (2011)结果, 考虑Shao等 (2010)断层模型根据地震波反演，此处计算未考虑地形和横向不均匀性，结果虽与观测有一定差别但整体反映了同震位移特征，尤其在远场方向吻合较好，这也反映计算的合理性.由于Shao等 (2010)断层模型在反演中未使用大地测量数据约束，笔者还考察了Tong等 (2010)根据InSAR和近场GPS数据反演的滑动模型，计算发现Tong等模型虽然在近场破裂面附近与GPS吻合更好，在远场与实测位移方向偏离过大 (图 3d).目前公布的滑动模型里，使用远场GPS作为约束的仅有Sladen和Owen (2010)模型，根据Vigny等 (2011)文献，该模型和近场GPS吻合度极差，综合模型计算结果和远、近场GPS吻合程度高，本文选用Shao等 (2010)滑动模型计算Maule大地震同震效应.

图 3

Fig. 3

图 3 同震位移比较 (a) 球面分层模型；(b) 平面分层模型；(c) 平面半空间模型；(d) 平面半空间模型. (a)-(c) 使用Shao等 (2010)模型，(d) 使用Tong等 (2010)模型.黑色、蓝色箭头为GPS观测值，红色、绿色箭头为计算值. Fig. 3 The comparison of co-seismic displacement between the calculation and observation (a) Spherical layered model; (b) Plane layered model; (c) Plane half-space model; (d) Plane half-space model. Model (a)-(c) based on slip distribution of Shao et al (2010), model (d) based on slip model of Tong et al (2010).The black and blue arrows indicate the GPS observation, while the red and green arrows indicate the calculated co-seismic displacement.

根据图 3b和图 3c，平面分层模型和平面半空间模型计算所得水平位移方向和GPS实测位移方向吻合较好，但这2个模型远场位移均比观测值大，图 3a显示球面分层模型远场位移和观测值更为接近.付广裕和孙文科 (2008)计算2004年苏门答腊大地震同震水平位移发现在远场平面半空间模型计算结果约是球位错模型结果1.8倍，笔者结果与其类似.地球曲率、层状结构对远场位移影响不可忽略：球面分层模型远场位移小于平面分层模型远场位移；平面分层模型远场位移大于平面半空间模型远场位移.

地球曲率、层状结构对远场应变影响明显.Maule地震是一个低角度逆冲地震，应变积累/释放主要发生在水平面，考察低、高纬度两种极端情况，采用3种模型计算对比本次地震引起的赤道和南极洲水平位移、应变场.地震断层走向正交方向同震水平位移最大，因此选取位于断层西北方的近赤道区域 (150°W-120°W, 10°S-0°S) 和断层东南方的近南极洲区域 (0°E-30°E, 80°S-70°S) 作为研究区域，研究深度15 km，并采用墨卡托投影.

在赤道区域球面分层模型远场位移、应变量级均小于平面分层模型和平面半空间模型，而平面分层模型和平面半空间模型在赤道的远场位移和应变量级十分接近 (图 4)，因此地球曲率对远场位移影响大于地球层状结构影响.对于2010年Maule地震，在低纬度远场，如果不采用球面分层模型，可能会造成位移和应变量值均被低估.

图 4

Fig. 4

图 4 智利Maule地震引起的赤道区域形变场 (a)-(c) 水平位移场；(d)-(f) 水平主应变场 (红线为主拉应变，黑线为主压应变)；(a)、(d) 球面分层模型结果；(b)、(e) 平面分层模型结果；(c)、(f) 平面半空间模型结果. Fig. 4 Deformation near the Equator (a)-(c) indicate displacement field and (d)-(f) indicate strain field (red lines are principal tensile strain, black lines are principal compression). (a) and (d) based on spherical layered model, (b) and (e) based on plane layered model, (c) and (f) based on plane half-space model.

根据石耀霖和朱守彪 (2006)文献，在高纬度地区采用直角坐标计算应变将带来较大量级误差，因为球坐标应变公式里有一项与余纬的余切相关，纬度越高，该项值越大.在南极洲区域，虽然球面分层模型远场位移值小于平面分层模型和平面半空间模型，但其远场应变值大于平面分层模型和平面半空间模型 (图 5).对于2010年Maule地震，在高纬度远场，如果不采用球面分层模型，可能会造成位移被夸大而应变被低估.

图 5

Fig. 5

图 5 智利Maule地震引起的南极洲区域形变场 (a)-(c) 水平位移场；(d)-(f) 水平主应变场 (红线为主拉应变，黑线为主压应变)；(a)、(d) 球面分层模型结果；(b)、(e) 平面分层模型结果；(c)、(f) 平面半空间模型结果. Fig. 5 Deformation near Antarctica (a)-(c) indicate displacement field and (d)-(f) indicate strain field (red lines are principal tensile strain, black lines are principal compression). (a) and (d) based on spherical layered model, (b) and (e) based on plane layered model, (c) and (f) based on plane half-space model.

2010年Maule地震错动方向近似为西北-东南，相对较大水平位移主要发生在地震断层西北、东南部 (图 6).南美大陆接近一半面积水平位移量级达到0.5 mm，基本上可以被现代GPS观测技术捕捉到，但水平位移较大区域海洋居多，一定程度上影响本次地震同震位移观测.计算近场垂直最大隆升为2.9 m，与Farías等 (2010)根据海岸线上红石灰藻相对位置给出的最大同震垂向位移2.5 m结果较为接近.根据cGPS台站记录到位移 (Vigny et al., 2011)，距离震中约3000 km的南美大陆东北部的法属圭亚那水平位移约为1 mm，南美大陆最南端南巴塔哥尼亚固定台站垂直位移小于1 mm，与计算吻合较好，也反映了本文选取滑动模型在远场具有合理性.总体上，远场垂直位移相对较小，约0.1 mm量级，并表现出正负相间的同心圆样式.

图 6

Fig. 6

图 6 2010年智利地震引起的全球同震位移场等值线图 (单位mm) (a) 水平位移；(b) 垂直位移. Fig. 6 Contour lines of co-seismic displacements caused by the 2010 Maule earthquake (unit:mm) (a) Horizontal displacements; (b) Vertical displacements.

地震发生时应变能释放，产生应力降，根据应力降可以了解区域的应力加/卸载状态.地震引起的15 km深度东西向正应力、南北向正应力和水平剪应力变化对智利和阿根廷中部影响均在10 kPa量级 (图 7a-7c).安第斯地区年应变率在10^-8~10^-7数量级 (Khazaradze and Klotz, 2003)，上地壳的弹性模量估计为70 GPa (Wortel and Cloetingh, 1983), 根据弹性假设东西向压应力年增长量约为0.7~7 kPa，Allmendinger等 (2006)估计阿根廷的应变率约为智利高应变区0.2倍，近似取阿根廷东西向压应力年增长量为1 kPa，本次地震造成阿根廷中部约5~10年构造挤压应力释放、阿根廷北部和南部1~5年构造压应力积累.

图 7

Fig. 7

图 7 同震应力变化 (a)、(d) 东西向正应力变化；(b)、(e) 南北向正应力变化；(c)、(f) 水平剪应力变化. (a)-(c) 考虑横向非均匀介质和地形；(d)-(f) 采用PREM分层介质、不考虑地形；(g)-(i) 为二者差值. Fig. 7 Co-seismic stress change (a), (d) E-W normal stress; (b), (e) S-N normal stress; (c), (f) Horizontal shear stress. Stress change (a)-(c) based on laterally heterogeneous earth model including topography; stress change (d)-(f) based on layered spherical earth model; (g)-(i) stress change difference between 2 models.

安第斯俯冲带地区三维构造复杂，介质横向非均匀性、海陆地形差异对同震应力降的影响不可忽视.计算区域内地壳厚度变化剧烈，大陆地壳厚度最大70 km, 而洋壳厚度小于10 km，计算深度15 km在大洋已经深入到地幔，在安第斯大陆对应于上、中地壳，必然存在横向上介质性质的巨大差异，此外该区海陆地形起伏亦十分剧烈.作为对比，不考虑介质非均匀性和地形起伏，以球面分层模型 (以PREM模型给出分层介质参数) 计算得到相同深度应力降 (图 7d-7f)，结果显示同震应力降被显著低估.在远场同震东西、南北向应力变化相对误差最高均超过100%，考虑到较远处绝对应力降很低的区域，由于绝对量很小，两种计算得到的相对差别可以非常大.在断层破裂面所处的近场，同震东西向应力变化相对误差最高95.4%，南北向应力变化误差最高90.8%，据图 7g-7i可以看出近场两种模型应力降的绝对差值可以达到10 kPa量级.

安第斯山脉西缘为地震多发的板块俯冲带，根据Bird (2003)板块边界数字模型，纳兹卡-南极洲转换断层历史上也发生过许多走滑型地震.Liu等 (2000)根据地质资料认为地壳缩短集中于安第斯山东缘褶皱冲断带，Brooks等 (2011)发现隐伏在玻利维亚安第斯山脉东缘之下一大段断层可能会在一次M_W8.7~8.9地震中破裂.因此，本文选取上述地震多发的板块边界、活动断层及被预言将发生破裂断层作为研究对象.

库仑应力变化计算所需断层走向、倾角和可能滑动方向根据USGS工作报告 (Costa et al., 2000; Lavenu et al., 2000) 和文献 (Allmendinger and González, 2010; Veloza et al., 2012) 确定，对多数具有丰富地震资料的断层按照震源机制解平均确定其滑动角，少数没有地震资料断层根据区域主应力方向 (Assumpcao, 1992) 推测其可能滑动方向，表 1给出接收断层信息.

表 1(Table 1)

表 1 接收断层 Table 1 Receiver faults 编号 断层名称 走向/(°) 倾角/(°) 滑动角/(°) 1 Estern Peru断层 131.6 37.0 98.6 2 Pasco断层 165.3 26.0 65.0 3 Pucallpa断层 158.6 30.0 94.4 4 Manu断层 120.4 30.0 93.3 5 Mande yapecua断层 175.1 60.0 98.3 6 La Dehesa断层 183.4 45.0 94.8 7 Sierra Chica断层 3.2 40.0 94.2 8 Liquiñe-Ofqui断层 14.5 80.0 161.4 9 秘鲁-智利海沟1 326.5 24.4 79.4 10 秘鲁-智利海沟2 315.9 27.1 85.4 11 秘鲁-智利海沟3 3.4 24.0 101.0 12 秘鲁-智利海沟4 11.0 18.0 83.0 13 秘鲁-智利海沟5 17.3 19.2 116.2 14 秘鲁-智利海沟6 12.8 21.4 90.2 15 秘鲁-智利海沟7 6.9 30.0 90.2 16 秘鲁-智利海沟8 6.6 30.0 90.2 17 纳兹卡-南极洲断层1 95.7 80.9 174.0 18 纳兹卡-南极洲断层2 96.3 82.3 175.3 19 纳兹卡-南极洲断层3 110.2 75.0 175.0 20 纳兹卡-南极洲断层4 95.0 84.6 183.7 21 纳兹卡-南极洲断层5 94.7 85.7 173.5

表 1 接收断层 Table 1 Receiver faults

据表 1，将同震应力降投影到断层面滑动方向上得到库仑破裂应力变化 (图 8).虽然构造加载方向并不严格沿每条断层滑动方向，仍可粗略估计断层库仑应力变化的等效构造加载时间，安第斯西缘应变率较高构造应力积累率约7 kPa·a^-1，安第斯东缘应变率较低构造应力积累率约1 kPa·a^-1.Maule地震所在断层 (13号断层)，库仑应力变化约为-100 kPa，相当于释放15年构造积累，板片俯冲造成应力积累在这次地震中得到释放，同一位置短期再发生强震可能性不大.与Maule地震紧邻断层 (12、14号断层) 库仑应力变化最大超过100 kPa，Maule北部Illapel地区 (12号断层) 发生了2015年智利M_W8.3地震，Illapel地区上一次8级以上地震是1943M_W8.2地震，该区8级以上地震复发周期约为80~90年 (Nishenko and Buland, 1987)，Illapel地区库仑应力变化相当于15年构造积累，因此Maule地震显著加速了Illapel地震发生.在Maule南部，单侧破裂长度达1000 km的1960年Valdivia M_W9.5地震破裂区已50余年未再次发生破裂，该破裂区对应14、15、16号断层，尤其14号断层库仑应力变化达到100 kPa，地震危险性增加.阿根廷境内褶皱冲断带La Dehesa断层 (6号) 和Sierra Chica断层 (7号) 库仑应力变化量级在0.4~1 kPa，相当于0.4~1年构造加载, 另一条走滑断层Liquiñe-Ofqui断层 (8号) 滑动方向和震源错动方向几乎正交，因此库仑应力变化的量级较小.转换断层 (17~21) 库仑应力变化以正值为主约0.4 kPa，约为0.1年构造积累.离震源较远的秘鲁、玻利维亚境内的断层库仑应力的加/卸载作用不明显.

图 8

Fig. 8

图 8 断层面上的库仑应力变化 黑线表示断层地表迹线，黑线上数字代表断层编号；绿色五角星代表 2010年Maule地震；黑色五角星代表 2015年Illapel地震. Fig. 8 Coulomb stress on faults Black line indicates surface trace of fault, number indicates fault number; Green star indicates 2010 Maule earthquake; Black star indicates 2015 Illapel earthquake.

在破裂面未知情况下常以最优破裂方向投影库仑应力变化解释余震分布，预测后续地震活动 (King et al., 1994).最优破裂方向计算需知道区域应力场，本文区域应力场按照Luttrell等 (2011)反演的Maule地震破裂区偏应力场确定，即超过静岩压部分最大水平偏应力27 MPa (垂直于海沟方向NE)、最小水平偏应力15 MPa (平行于海沟方向NE)，垂直偏应力0 MPa.为了形象表示计算得到的最优破裂方向，我们采取震源球方式表示，最优破裂方向基本上为逆冲型 (图 9a黑色震源球).最优破裂方向计算依赖于区域应力场和同震应力降叠加，应力降量级低于区域应力场，因此最优破裂方向表现为逆冲型较为合理.

图 9

Fig. 9

图 9 库仑应力变化与余震空间分布 (a) 黑色沙滩球指示最优破裂方向，红色沙滩球指示余震的震源机制; (b) 投影到主震破裂面的库仑应力变化; (c) 投影到最优破裂方向的库仑应力变化. Fig. 9 Coulomb stress change and the spatial distribution of the aftershock (a) Black beach balls indicate the optimal planes, red beach balls indicate aftershocks; (b) Coulomb stress change on the actual fault planes; (c) Coulomb stress change on the optimal planes.

由于后续陆续发生强余震会改变应力分布情况，在此仅讨论主震后一个月内面波震级大于4级的余震，数据来源于USGS (http://neic.usgs.gov/[2015-10-20])，其中181条数据 (图 9a红色震源球) 具有明确震源机制解，且多数为逆冲型，与代表最优破裂方向的黑色震源球吻合较好 (图 9a)，说明本文采用区域应力场有一定合理性.部分余震为正断或走滑，和按照最优破裂方向给出的断层破裂类型不完全一致.原因如下：(1) 有些余震可能发生在隐伏的先存断层上，接收断层不能取为最优破裂方向；(2) 实际的区域应力场远比本文采用的均匀应力场复杂.

余震平均深度约30 km，因此计算深度取为30 km，投影到主震破裂面的库仑应力变化显示约63.4%余震落在库仑应力变化大于库仑应力触发阀值 (Harris, 1998)10 kPa区域 (图 9b)，投影到最优破裂方向的库仑应力变化显示70.9%余震落在库仑应力变化大于10 kPa区域 (图 9c).如果介质参数采用PREM模型且不考虑地形，投影到主震破裂面和最优破裂方向的库仑应变变化显示，分别有52.9%和66.2%余震落在库仑应力变化大于10 kPa区域.然而，真实区域应力场并不是如此简单均匀应力场，余震震源机制也有很大不确定性，不同震源模型库仑应力差别也很大 (石耀霖和曹建玲, 2010)，如何研究主震对余震序列触发还有很多问题需要解决.

2010年3月11日在Pichilemu镇附近区域15 min内接连发生了M_W6.9和M_W7.0两次大的正断型地震，这两个大地震被认为是马乌拉特大地震的余震 (Ruiz et al., 2014)，位于主震以北230 km处.这种大的板内地震在智利俯冲区十分罕见，在Pichilemu地区没有历史地壳浅表地震活动记载.Farías等 (2011)认为该区域发育有新生代伸展构造，存在隐伏的倾向SW的正断层，并根据高异常的V_p/V_s认为马乌拉地震产生的应力变化加强了该区流体流通，最终触发了这两个余震.

大地震发生后的大量余震会产生可以观测到的地表位移场，余震会对周围区域的应力场产生扰动 (万永革等，2005).因此对于发生时间相对较晚的M_W7.0余震事件，在计算中需计入发生在它之前的M_W6.9强余震影响，M_W6.9强余震错动方向按照表 2给出，位错大小根据震级-位错经验公式 (Wells and Coppersmith, 1994) 给出.为了考虑流体活动异常造成的孔隙水压变化，分别计算视摩擦系数μ′=0、μ′=0.4和μ′=0.8情况下余震破裂面库仑应力变化.不同视摩擦系数下，两次余震破裂面上库仑应力变化均超过了库仑应力触发阀值0.01 MPa，由此可见本次俯冲带特大地震的应力释放触发了Pichilemu正断型M_W6.9和M_W7.0余震.

表 2(Table 2)

表 2 两次强余震破裂面上的库仑应力变化 Table 2 Coulomb stress change on 2 strong aftershocks 震中位置 震级 (MW) 深度 /km 走向 /(°) 倾角 /(°) 滑动角 /(°) ΔCFS/MPa 北纬/(°) 东经/(°) μ′=0 μ′=0.4 μ′=0.8 -34.54 -72.11 6.9 12.9 144 55 -90 0.24 0.83 1.28 -34.53 -72.13 7.0 16.3 159 86 -93 1.84 2.30 2.76 注：数据来源美国哈佛大学Global CMT.

表 2 两次强余震破裂面上的库仑应力变化 Table 2 Coulomb stress change on 2 strong aftershocks

本文应用高性能并行有限元方法，建立包含地表地形、Moho面起伏的大规模非均匀椭球地球模型，根据2010年智利Maule特大地震震源破裂分布计算了球坐标系下同震形变场、应力场，以库仑应力变化分析周围断层地震活动性和主震对余震的触发关系.主要得到以下结论.

(1) 地球曲率和层状结构对同震水平位移的影响不可忽略.对于2010年Maule地震，如果不采用球面分层模型，在低纬度远场可能会造成位移和应变量值均被低估，在高纬度远场可能会造成位移被夸大而应变被低估.

(2) 相对较大水平位移主要发生在地震断层的西北和东南部，近场垂直最大隆升为2.9 m，远场垂直位移相对较小，基本上在0.1 mm量级.

(3) 本次地震约造成了阿根廷中部5~10年构造挤压应力释放，阿根廷北部和南部1~5年构造应力积累.如果不考虑介质横向不均匀和地形, 同震应力降被显著低估，断层破裂面近场处同震东西向应力变化误差最高95.4%，南北向应力变化误差最高90.8%.Illapel地区库仑应力变化约10 kPa，相当于15年构造积累，Maule地震加速2015年Illapel地震发生.

(4) 以最优破裂方向投影得到库仑应力变化，经统计70.9%余震落在库仑应力变化大于10 kPa区域.本次特大地震应力释放触发了Pichilemu附近M_W6.9和M_W7.0两次大余震.

俯冲带地区构造加载强烈，构造运动对地震发生起主导作用，但特大地震引起的同震应力变化不容忽视，同震库仑应力变化使近场周边断层提前或延缓滑动时间约为1~10年量级，而俯冲带大地震的复发周期往往不足百年，因此俯冲带地区大地震引起的库仑应力变化仍然值得重视.

地震破裂过程的复杂性以及不同作者所采用数据、方法的不同，带来了滑动模型反演不确定性，造成同震应力应变计算不确定性.目前公开发表滑动模型差异集中如下：主要滑动区 (以5 m滑动量等值线定义的区域) 范围不一样，远震波形数据反演的模型 (Lay et al., 2010) 滑动区延伸到海沟处，大地测量数据反演的模型 (Tong et al., 2010; Vigny et al., 2011) 主要滑动区停止于5~10 km处，大地测量、地震波形数据联合反演的滑动模型主要滑动区停止于10~20 km (Delouis et al., 2010) 处；不同模型的南部滑动区中心和Arauco Peninsula的相对位置也有所区别，Delouis等 (2010)的模型中心在其北部，Tong等 (2010)的模型中心正好在Peninsula，Vigny等 (2011)的模型中心在其西部.Vigny对不同滑动模型的计算值和近场GPS的吻合程度做了对比，这样的对比具有一定片面性，并未比较远场GPS的吻合程度，例如，本文所选用的Shao等 (2010)的滑动模型与近场GPS的吻合程度虽然不如Tong等包含近场GPS约束给出的滑动模型高，但Shao等 (2010)模型与远场GPS的吻合程度较好且远高于Tong等 (2010)模型，这也是本文选用Shao滑动模型的重要原因.目前尚无模型在反演中同时纳入远、近场GPS约束，例如，根据Vigny等 (2011)文献，Sladen和Owen (2010)仅考虑远场GPS的模型与近场GPS吻合较差，据本文结果仅考虑近场GPS的模型 (Tong et al., 2010) 与远场GPS残差也很大，未来需要利用更详实的远、近场大地测量资料和地震波数据反演滑动模型，在此基础上的计算更趋可靠.

本文工作的意义在于同震应力、应变计算的模型更接近真实地球，考虑地表地形起伏、莫霍面起伏、地球介质横向不均匀性以及地球曲率，给出Maule特大地震同震响应.本文还计算了地表地形起伏和介质非均匀性对近场同震水平和垂直位移的影响，其百分比分别为41%和95%，林晓光和孙文科 (2014)计算日本东北大地震引起的同震位移，发现海山地形剧变处，介质非均匀性和地形效应对水平和垂直位移的影响分别为34%和92%，俯冲带地形起伏较大并具有强介质非均匀性，二者对同震形变的影响不能忽略.本文未定量讨论地球曲率和地球分层结构对库仑应力变化影响，未来还需要更多定量计算.值得注意的是，本文所得地球曲率和分层结构对同震位移场的影响是基于2010年智利Maule地震算例的结果，要得到更具普适意义结论需要做更多工作.此次研究采用线弹性本构，未考虑黏弹弹效应，但黏弹效应只在研究时间大于下地壳、地幔应力松弛的的特征时间时更加重要，本文仅探讨震后较短时间内对地震活动性的影响，黏弹效应可以近似忽略.此外，震源模型的不确定性、孔隙流体对视摩擦系数的影响等因素在未来的工作中需要做进一步考虑.