地震学报  2012, Vol. 34 Issue (6): 727-740
引用本文
刘 鎏, 魏东平. 2012. 中国大陆及邻区板内应力场的数值模拟及动力机制探讨. 地震学报, 34 (6): 727-740.
Liu, Liu, Wei, Dongping. 2012. Numerical simulation of the intraplate stress field in China mainland and adjacent areas and its dynamic implications. Acta Seismologica Sinica , 34 (6): 727-740.
中国大陆及邻区板内应力场的数值模拟及动力机制探讨
刘 鎏1,2, 魏东平1,2     
1. 中国北京100049中国科学院计算地球动力学重点实验室;
2. 中国北京100049中国科学院大学地球科学学院
收稿日期 2011-11-22收到初稿， 2012-01-16决定采用修改稿.
基金项目： 国家自然科学基金项目(40874047, 40628004)资助.
通讯联系人,E-mail:dongping@ucas.ac.cn
摘要：中国大陆位于欧亚板块的东南部， 受到印度板块、 太平洋板块和菲律宾海板块的碰撞挤压与俯冲作用， 其构造应力场形态和动力学机制相当复杂. 本文采用伪三维有限元方法， 以世界应力图 2008年版本数据(WSM2008)的应力方向和应力型两类指标作为主要约束， 对中国大陆及邻区的动力驱动机制进行数值模拟， 给出了中国大陆周边地区板块边界力的大小和方向估计. 同时对3个典型情况的数值模型进行了分析. 结果显示, 软流层静压推力对该区域构造应力场影响相对较小， 板块边界力作用则起主导作用； 印度板块在喜马拉雅造山带对欧亚板块的碰撞控制了中国大陆地区应力场的基本形态, 是形成川滇地区走滑型地震为主的重要原因； 琉球海沟—南海海槽俯冲带边界力显示了挤压-张性的分段特性, 贝加尔裂谷表现为拉张作用. 进一步的分析表明， 中国大陆大部分区域内最大水平剪应力分布图像与该地区地震辐射能量密度的分布存在较好的空间正相关性.
关键词： 中国大陆     构造应力场     WSM2008     数值模拟     软流层静压推力    
Numerical simulation of the intraplate stress field in China mainland and adjacent areas and its dynamic implications
Liu1,2, Liu1,2     
1. Key Laboratory of Computational Geodynamics, Chinese Academy of Sciences, Beijing 100049, China;
2. College of Earth Science, University of Chinese Academy of Sciences, Beijing 100049, China
Abstract: Located in the southeast of Eurasian plate, China mainland is subject to collision, extrusion and subduction of Indian plate, Pacific plate and Philippine Sea plate, so the tectonic stress field there and its dynamical mechanism are very complex. Using a pseudo 3D finite element method and with the orientation and regime indicators from World Stress Map 2008 (WSM2008) database as principal constraints, this paper numerically simulated driving forces on China mainland and adjacent areas, and estimated the magnitude and direction of plate boundary forces. We analyzed three typical models. The results show that the asthenospheric static push forces have small influence on the tectonic stress field, while the plate boundary forces play a significant role. The collision of Indian plate with Eurasian plate in the Himalayan orogenic belt controls the basic pattern of stress field in China mainland, and is an important cause of lots of strike-slip earthquakes in Sichuan-Yunnan region. The boundary forces of Ryukyu Trench-Nankai Trough subduction zone display push-to-pull segmental feature. Baikal rift is under the action of pull force. Further analysis reveals that in most areas of China mainland the distribution of maximum horizontal shear stress shows a good positive spatial correlation with the radiated seismic energy density in this region.
Key words: China mainland     tectonic stress field     WSM2008     numerical simulation     asthenospheric static push force    

引言

中国大陆及邻区位于欧亚板块东南部， 东邻环太平洋地震带， 南接欧亚地震带， 受到印度板块、 太平洋板块和菲律宾海板块碰撞、 挤压和俯冲的影响． 各种观测资料显示， 该地区表现出很强的地震活动性， 表明在全球性板块作用力及区域构造应力作用下， 应变积累到一定程度， 地球内部岩层会发生破裂与错动， 导致地震发生. 因此研究构造应力场的分布和演化规律是帮助我们了解地震成因机理的重要手段之一．

Richardson等(1979)分析全球最大水平主应力分布情况， 给出了与板块运动相关的驱动力； Bird等(2008)采用最新板块模型对板块驱动力加以详细说明． Zoback等(1989)指出印度板块对欧亚板块的碰撞强烈影响东亚地区的应力图像， 板块驱动力控制了板块内部的应力场分布； Zoback和Magee(1991), Coblentz和Richardson (1996)得到了类似结论． Zoback(1992)指出构造应力场包含一阶及二阶分量， 一阶分量主要来源于板块边界； Flesch等(2001)， Reynolds等(2002)证实了这个结论．

近30年来， 我国学者从不同角度, 采用不同方法对中国大陆及邻区开展了构造应力场的研究工作． 邓起东等(1979)认为, 中国大陆构造应力场是受印度板块和太平洋板块联合作用控制的， 有明显的分区特点． 汪素云和陈培善(1980)、 汪素云等(1996)及许忠淮等(1992)反演了中国大陆周边板块作用力的大小， 表明印度板块作用力最大， 其次是太平洋板块， 菲律宾海板块作用力最小． 张东宁和许忠淮(1999)研究表明, 青藏高原物质的侧向扩展是中国大陆构造特征的主要动力， 菲律宾海板块对欧亚板块的挤压作用很微弱． 谢富仁等(2004)总结了中国大陆及邻区构造应力场的基本特征， 给出了划分应力场分区的原则． 朱守彪和石耀霖(2006)反演了中国大陆及邻区构造应力场的成因， 表明板块边界力起主导作用． 徐纪人和赵志新(2006)及徐纪人等(2008)认为， 印度板块对欧亚板块的碰撞挤压主导了中国中西部地区应力场， 南北地震带是印度板块、 菲律宾海板块在中国大陆影响控制范围的分界线． 梁海华等(1999)对中国大陆构造应力场与大震周期进行了模拟， 表明可用应力衰减图像解释不同地区大震周期的差异性； 焦明若等(1999)及章纯(2007)研究认为， 剪应力、 剪应变能等与地震活动性有关．

Wei (1997)对欧亚板块构造应力场进行模拟计算， 主要考虑了一阶作用力， 同时以世界应力图计划的观测应力取向和应力型两方面指标作为约束， 推测出中国大陆及邻区的边界作用力在1013 N/m量级． 本文在Wei (1997)研究方法的基础上， 以世界应力图计划最新观测数据库WSM2008的最大水平主应力(SH max )方向和应力型作为主要约束条件， 进行数值模拟计算， 分析了中国大陆及邻区边界力对板内应力场的影响， 并利用模拟结果得到了最大剪应力分布， 探讨了其与该地区地震辐射能量密度分布的关系．

1 模型数据及计算方法 1.1 岩石层模型

以PB2002板块模型(Bird, 2003)和MORVEL板块模型(DeMets et al, 2010)为参考， 建立了中国大陆及邻区岩石层三维弹性模型(命名为模型CHN)， 区域包括阿穆尔板块(Amurian Plate, 简写为AM)、 扬子板块(Yangtze Plate, 简写为YA)、 冲绳板块(Okinawa Plate, 简写为ON)、 波斯—青藏高原—缅甸(Persia-Tibet-Burma, 简写为P-T-B)造山带的中国大陆部分(图 1a)． 对该区域采用兰伯特等积方位投影(Lambert azimuthal equal-area projection)建立平面有限元网格， 模型共包含1 841个节点， 3 505个三角形单元(图 1b)， 投影中心为107°E， 36°N．

图 1 模型CHN研究区域分布. (a) PB2002板块模型(灰色区域为P-T-B造山带); (b) 本文研究模型区域和网格划分AM: 阿穆尔板块; EU: 欧亚板块; IN: 印度板块; OK: 鄂霍次克海板块; ON: 冲绳板块; PA: 太平洋板块; PS: 菲律宾海板块; YA: 扬子板块Fig. 1 Regional distribution of the research model CHN. (a) PB2002 Plate model (the gray area is P-T-B orogen); (b) Research model areas and meshing AM—Amurian, EU—Eurasian, IN—Indian, OK—Okhotsk, ON—Okinawa, PA—Pacific, PS—Phillippine Sea and YA—Yangtze plates

依据板块性质不同， 岩石层板块包括大陆型和海洋型两部分. 大陆型板型由上地壳、 下地壳、 固体上地幔3层组成， 地壳厚度数据来源于CRUST2.0模型(Bassin et al, 2000)， 将3—5层归为上地壳， 6—7层归为下地壳， 固体上地幔厚度统一取为100 km； 海洋型板块岩石层厚度与海洋板块年龄有关， 其厚度依据板块冷却模式(McKenzie, 1967)给定， 其中板块年龄数据来自全球海洋等时线年龄(Müller et al, 2008)， 板块的渐进厚度取为95 km (Stein,Stein, 1992)．

给定岩石层板块的地表温度， 利用全球大地热流最新数据(http://www. und. edu/org/ihfc)进行约束， 计算岩石层的三维温度分布． 对大陆节点， 采用一维稳态热传导模型； 对大洋节点， 采用板块冷却模式， 利用公式(1)计算岩石层的三维密度分布.

式中， ρ0是0℃时岩石的密度， α为热膨胀系数． 模型主要参数的具体数值见表 1．
表 1 建立模型的主要参数及其参考值 Table 1 Principal parameters and their reference values in modeling

岩石层板块几何结构与软流层静压推力和应力场垂直分量直接相关(魏东平, 2000)． 本研究中假定整个模型处于重力均衡状态， 以ETOPO5数据作为地形约束， 调整差别较大区域的参数(An, Shi, 2006; Conrad, Lithgow, 2006; Li et al, 2006; 胥颐等, 2006; 李廷栋, 2010)， 得到与实际接近的岩石层模型．

1.2 计算方法

陈颙(1988)与Saito(1995)的研究表明， 可以假定应力场的两个水平主应力分量(最大水平主应力SHmin、 最小水平主应力SHmin)与垂直主应力分量(SV)互相解耦， 因此可利用平面有限元方法计算应力场的水平分量， 垂直方向计算岩石层自重造成的应力场， 同时考虑不同边界类型作用力、 板块基底拖曳力等共同作用， 使得整个模型处于动力平衡状态．

采用WSM2008质量为AD级的观测数据作为主要约束条件， 模型内部有1669个数据， 分布较均匀， 因此可以在大尺度上反映板内应力场的形态． 采用Wei (1997)、 魏东平, 2001的方法计算单元内平均应力型， 共有620个单元可用. 其中214个为挤压型单元， 278个为走滑型单元， 128个为扩张型单元； 包含应力方向的单元共688个可用， 观测应力场分布见图 2． 对于应力型数据， 如果计算值与统计值相同， 则认为单元通过检验； 对于应力方向数据， 如果计算值与统计值偏差不大于20°， 则认为单元通过检验(Zoback, 1992; Wei,1997)．

图 2 观测应力场应力方向指标分布(a)和观测应力场应力型指标分布(b)Fig. 2 Distribution of observed stress orientation indicators (a) and observed stress regime indicators (b)

研究表明， 影响板块运动和板内应力场的作用力主要为洋脊推力、 碰撞带推力、 俯冲带作用力及地幔拖曳力等(Forsyth, Uyeda, 1975; Cloetingh, Wortel, 1986； Coblentz,Richardson, 1996; Pacanovsky et al， 1999)． 本研究中， 以软流层静压推力取代传统的洋 脊推力. 在边界上施加不同类型的作用力， 最后用地幔拖曳力平衡其它各力的合力(包括力平衡与力矩平衡， 见公式(2))， 保证岩石层系统处于动力平衡状态(魏东平, 2000)． 板内应力场与驱动板块运动的作用力有很好的对应关系(Cloetingh，Wortel, 1986; Zoback, 1992)， 因此参考板块相对运动方向， 在一定方向区间内寻找力的最佳作用方向， 同时不断调整边界力的大小， 进行模拟计算， 综合考虑应力型和应力方向两方面条件约束， 得到效果较好的模型(冯向东等， 2006； 皇甫海泉,魏东平， 2010)．

式中， l1为碰撞带长度， l2为俯冲带长度， r为相应点的矢径，Fc为单位长度碰撞力大小， F s为单位长度板块拉力或吸力大小， M d为地幔拖曳力矩．

1.3 模型边界类型

参照Bird(2003)给出的边界类型， 本研究模型边界类型见图 3a． 边界带Ⅰ ： 喜马拉雅山带， 为大陆汇聚边界(continental convergent boundary, 简写为CCB)； 边界带Ⅱ-1、 Ⅱ-2： 印度板块与欧亚板块边界， 前者为大陆转换边界(continental transform fault, 简写为CTF)， 后者为大陆汇聚边界； Ⅱ-3： 巽达板块与欧亚板块边界带， 大陆汇聚边界； 边界带Ⅲ-1： 巽达板块与扬子板块边界， 大陆转换边界； 边界带Ⅲ-2： 台湾—耶亚玛海脊， 大陆汇聚边界； 边界带Ⅲ-3： 琉球海沟—南海海槽边界， 俯冲边界(subduction zone, 简写为SUB)； 边界带Ⅳ： 阿穆尔板块东边界， Ⅳ-1、 Ⅳ-3为大陆汇聚边界， Ⅳ-2为海洋汇聚边界(oceanic convergent boundary, 简写为OCB)； 边界带Ⅴ-1： 斯塔诺夫山脉西段， 大陆转换边界， Ⅴ-2： 贝加尔裂谷， 大陆扩张边界(continental rift boundary, 简写为CRB)； P-T-B造山带和模型西边界位于欧亚板块内部， 不是板块边界， 故设定为自由边界(Free Bou)．

2 模拟结果 2.1 各边界带作用力初步估计

单独在某一边界段施加作用力， 其它边界段设置为自由边界， 即施加作用力为0， 对边界带的作用力进行初步估计(图 3b)．

图 3 中国大陆及邻区模型边界类型(a)和初步估计得到的模型边界力大小范围(b)CCB: 大陆汇聚边界; CTF: 大陆转换边界; CRB: 大陆扩张边界; OCB: 海洋汇聚边界; SUB: 俯冲边界; Free Bou: 自由边界Fig. 3 Boundary types of the model of China mainland and adjacent areas (a) and estimated magnitude of boundary forces in the model (b)CCB: Continental convergent boundary; CTF: Continental transform fault; CRB: Continental rift boundary; OCB: Oceanic convergent boundary; SUB: Subduction zone; Free Bou: Free boundary

边界带Ⅰ ： 喜马拉雅山带， 施加作用力1.0×1013—14.0×1013 N/m， 作用方向NE0°—30°． 当作用力大于5.0×1013 N/m时， 应力型在喜马拉雅山脉地区由走滑型变为挤压型， 在川滇地区转变为走滑型， 与实际接近． 随着作用力增大， 应力方向单元检验通过量在青藏高原地区逐渐增加， 初步推断该边界带的作用力为6.0×1013—9.0×1013 N/m(对于170 km厚的岩石层， 应力水平约相当于350—500 MPa, 下同). 这与Cloetingh和Wortel(1986)， Kincaid和Silver (1996)， 以及Wei (1997)给出的几百MPa量级类似， 但比Coblentz等(1998)和Copley等(2010)给出的几十MPa量级要大． 印度板块对欧亚板块的碰撞挤压影响到华北甚至东北及其以北地区的应力场， 主导了中国大陆构造应力场， 这与诸多研究结果类似(张东宁， 1994； 傅容珊， 2000； 朱守彪，石耀霖， 2006)．

边界带Ⅱ-2、 Ⅱ-3： 施加作用力0.5×1013—3.0×1013 N/m， 当作用力增大， 应力型通过量近乎稳定， 但应力方向通过量逐渐下降， 初步推断该边界带作用力不大于2.0×1013 N/m(应力水平150 MPa)．

边界带Ⅱ-2、 Ⅲ-1： 该边界段以大陆转换边界为主， 暂不考虑施加作用力．

边界带Ⅲ-2： 台湾—耶亚玛海脊， 施加作用力1.0×1013—6.0×1013 N/m， 方向NW30°—60°． 当作用力达到3.5×1013 N/m后， 应力型和应力方向通过量趋于稳定． 初步推断该段作用力大小为3.5×1013—5.0×1013 N/m(应力水平350—500 MPa)．

边界带Ⅲ-3： 该边界区域构造环境复杂， 应力方向和应力型数据在不同区域内均差别较大， 因此不同边界段作用力类型可能不同． 施加作用力-4.0×1013—2.0×1013 N/m， 方向NW30°—60°. 当作用力为正值时， 应力型和应力方向通过量均下降， 与实际不符； 考虑作用力为负值， 即为拉张作用(Wei, 1997; 许忠淮,吴少武, 1997； 臧绍先,宁杰远, 2002)． 当作用力为-2.0×1013 N/m时， 应力方向通过量达到最大， 华北和华东地区应力型由挤压型转变为走滑型， 与观测结果相同． 初步估计琉球海沟南段作用力为正， 大小范围为0.5×1013—1.5×1013 N/m(应力水平50—150 MPa). 其它边界段作用力范围为-1.0×1013—-3.0×1013 N/m(应力水平200—300 MPa)．

边界带Ⅳ ： 该边界带为碰撞带， 施加作用力1.0×1013—5.0×1013 N/m， 作用方向NW75°—90°． 当作用力大于3.0×1013 N/m时， 应力方向在日本岛地区通过量提高， 东北及华北地区几乎全部通过； 应力型通过量近似稳定． 初步估计日本岛边界段作用力为3.0×1013—5.0×1013 N/m(应力水平300—500 MPa)， 其它边界段为1.0×1013—3.0×1013 N/m(应力水平75—250 MPa)． 该结果表明， 鄂霍次克海板块与阿穆尔板块的碰撞影响到中国大陆东北及华北地区的应力场．

边界带Ⅴ ： 两个边界力同时作用， 作用力-1.0×1013—-4.0×1013 N/m， 作用方向NW30°—60°. 当作用力在Ⅴ-1段为-2.0×1013 N/m(水平应力140 MPa)， 在Ⅴ-2段为-3.0×1013—-4.0×1013 N/m(应力水平200—300 MPa)时， 东北地区出现走滑型单元， 应力型和应力方向通过量均较高且近似稳定， 表明贝加尔裂谷影响着阿穆尔板块的应力场分布．

板内应力场是多种力共同作用的结果， 我们对此进行了大量的模拟计算， 选取了3个有代表性的模型， 来研究不同类型边界对模型内部应力场的影响．

2.2 3个典型数值模型分析 2.2.1 模型CHN1： 仅考虑软流层静压推力

仅考虑软流层静压推力作用， 在这种情况下， 软流层静压推力本身达到力平衡， 仅存在很微小的板块基底拖曳力来平衡软流层静压推力导致的力矩作用． 通过应力方向检验的单元数为137， 通过率19.9%， 通过单元分布不均匀(图 4a)． 帕米尔高原、 青藏高原东边界、 琉球群岛、 中国大陆东北部等地区有少量单元通过； 青藏高原地区以近SE60°—80°为主， 与实际的NNE向为主差别较大； 中国大陆中东部地区以近南北向为主， 与实际的近东西向差别甚远； 台湾及琉球群岛北部以近南北向为主， 与实际差别60°左右， 日本岛地区以NW45°为主， 与观测的近东西向差别45°左右．

图 4 模型CHN1的检验结果(只考虑软流层静压推力)(a) 应力取向检验结果； (b) 应力型检验结果Fig. 4 Simulated results of model CHN1 (only the asthenospheric static push is considered) (a) Test result of SH max orientation; (b) Test result of stress regime

通过应力型检验的单元数为222， 通过率为35.8%. 通过检验的单元大部分位于105°E 以西地区和琉球群岛(图 4b)． 帕米尔高原和塔里木盆地以北地区以走滑型和挤压型为主， 通过率超过70%； 青藏高原中西部地区以走滑型为主， 东部以挤压型为主， 通过量超过一半， 但喜马拉雅山带挤压型过少， 与实际不符； 琉球群岛北部及日本岛地区以挤压型为主， 与实际走滑型和扩张型为主不符， 通过率约30%， 其它地区通过量较低．

模型CHN1结果显示， 仅考虑软流层静压推力时， 通过检验的单元数较少， 模拟的应力场与观测应力场差别较大， 表明软流层静压推力对该区域构造应力场影响较小．

2.2.2 模型CHN2： 考虑软流层静压推力、 边界带阻力及板块基底拖曳力

在模型CHN1的基础上加入边界带作用力进行模拟计算， 考虑了除扩张和走滑型边界外的所有边界带作用力．

根据2.1小节对模型边界力大小和方向的初步估计， 在边界带Ⅰ施加8.0×1013 N/m(应力水平450 MPa)的作用力， 作用方向NE10°—20°； 在边界带Ⅱ-1与Ⅱ-2施加2.0×1013 N/m(应力水平140 MPa)的平均作用力， 作用方向NE45°； 在边界带Ⅱ-3施加0.5× 1013 N/m(应力水平50 MPa)的作用力， 作用方向NE15°； 在边界带Ⅲ-2施加4.0×1013 N/m (应力水平400 MPa)的作用力， 作用方向NW45°； 在边界带Ⅲ-3施加大小不同的作用力， 从1.0×1013 N/m(应力水平100 MPa)到-2.0×1013 N/m(应力水平200 MPa)， 作用方向NW50°； 在边界带Ⅳ-1施加4.0×1013 N/m(应力水平400 MPa)的作用力， 作用方向EW；在边界带Ⅳ-2施加2.0×1013 N/m(应力水平150 MPa)的作用力， 作用方向NW75°； 在边界带Ⅳ-3施加1.0×1013 N/m(应力水平75 MPa)的作用力， 见图 5a所示．

图 5 模型CHN2的边界带作用力和模拟结果(a) 模型CHN2 施加的边界力； (b) 应力方向检验结果； (c) 应力型检验结果Fig. 5 Boundary forces of model CHN2 and simulated results(a) Boundary forces of model CHN2; (b) Test result of SHmax orientation;(c) Test result of stress regime

研究区域内通过应力方向检验的单元数为275， 通过率为40.0%， 比模型CHN1提高了近1倍(图 5b)． 中国青藏高原、 华北和东北南部及台湾地区， 日本岛南部， 蒙古地区通过量增加较多； 帕米尔高原及贝加尔裂谷区域通过量增加较少； 中国华中华东地区为近南北向， 与观测的近东西向差别近90°； 中国川滇地区通过量较少， 与观测方向差别约60°—90°．

通过应力型检验的单元数为248， 通过率40.0%， 与模型CHN1相比提高较少(图 5c)． 帕米尔高原地区通过量增加； 喜马拉雅山带完全变成挤压型， 与观测一致， 青藏高原北部地区部分单元由走滑型转变为挤压型， 通过量降低； 川滇地区由挤压型转变为走滑型， 通过量增加较多； 蒙古西部地区由走滑型变成挤压型， 通过量未明显变化； 琉球群岛中部及韩朝地区部分单元由挤压型变成走滑型， 与观测相同， 其它地区变化较小．

由模型CHN2分析可知， 板块边界带作用力对中国大陆及邻区板内应力场影响巨大．

2.2.3 模型CHN3： 考虑软流层静压推力、 所有边界带作用力、 板块基底拖曳力和局部调整的共同作用

在模型CHN2的基础上， 考虑了冲绳海槽和贝加尔裂谷的扩张作用， 边界带Ⅱ-1、 Ⅴ-1 的等效作用力， 施加作用力如图 6a所示(空心箭头为等效作用力)． 通过应力型检验的单元数为260， 通过率42.0%(图 6b)； 通过应力方向检验的单元数为293， 通过率为42.6%(图 6c). 两者的检验通过率几乎相同. 与模型CHN2相比， 两者通过量均有少量提高． 东北和东海地区单元由挤压型变成走滑型， 与实际吻合. 但华北华东地区应力型仍为挤压型， 与实际走滑型为主不吻合． 川滇地区的应力型与实际相同， 但应力方向与实际差别约60°， 这与该地区复杂的地质构造和大量的断裂带有关．

图 6 模型CHN3的边界带作用力和模拟结果(a) 模型CHN3施加的作用力； (b) 应力方向检验结果； (c) 应力型检验结果；(d) 地幔拖曳力方向与板块绝对运动方向关系图Fig. 6 Boundary forces of model CHN3 and simulated results(a) Boundary forces of model CHN3; (b) Test result of SHmax orientation; (c) Test result of stress regime; (d) Relationship between mantle drag force direction and absolute plate motion direction

地幔拖曳力平衡板块的运动， 在一些板块内与板块绝对运动方向有较好的相关性(Richardson et al, 1979; Coblentz,Richardson, 1996; Bird et al, 2008)． 模型CHN3计算得到的地幔拖曳力以近南北向为主， 并向东偏转， 与Lithgow-Bertelloni和Guyyn(2004)给出的结果类似； 板块绝对运动模型参考HS2-NUVEL 1A(Gripp,Gordon, 1990)及HS3-NUVEL 1A(Gripp,Gordon, 2002)， 两者计算得到的模型绝对运动方向类似， 从西到东模型的绝对运动方向发生北向偏转(图 6d)． 地幔拖曳力方向与板块绝对运动方向斜交， 相关性很差. 主要因为中国大陆及邻区构造应力场受洋脊推力影响较小， 板块边界力起主导作用， 尤其是喜马拉雅山带边界的作用力较其它边界带作用力大得多， 地幔拖曳力主要用来平衡此边界带作用力； 其次模型的绝对运动速度较小， 且存在很大的不确定性， 计算的运动方向与真实情况存在差别(Zoback et al, 1989; Wei, 1997)．

2.3 模拟结果与地震相关性分析

地震发生时释放出巨大能量， 其中一部分能量以地震波形式传播出去， 称为地震辐射 能量． 计算面波震级辐射能量采用Gutenberg-Richter公式lgEs=1.5MS+4.8(Gutenberg,Richter, 1956)， 矩震级采用公式lgEs=1.5MW+4.25(Hanks,Kanamori, 1979; Choy,Boatwright, 1995)． 地震目录数据取自USGS， 震级大于3.0， 时间范围1973—2011年． 定义地震辐射能量密度为单位面积内的地震辐射总能量， 以此表明不同地区地震活动性的强度．

实践表明， 岩石受到的剪应力超过其强度时就会发生破裂， 导致地震发生． 选取了模型CHN3的数据结果对中国大陆地区进行分析， 计算其最大剪应力值． 从图 7对比可看出， 中国大陆102°E以西和华北、 东北局部地区等地震辐射能量密度高值区， 其对应的最大水平剪应力也为高值区， 华中、 华东、 鄂尔多斯盆地等地震能量密度低值区对应的最大剪切应力也为低值区； 塔里木、 准格尔等盆地和东北部分地区处于最大剪应力相对高值区， 但由于其位于稳定的古板块构造中， 地震发生较少， 地震能量密度较低．

图 7 USGS地震辐射能量密度分布(a)和模型CHN3最大水平剪应力分布(b)Fig. 7 Distribution of radiated seismic energy density (a) and maximum horizontal shear stress of model CHN3 (b)

3 讨论与结论

以上模拟结果表明， 在软流层静压推力、 板块边界带作用力及板块基底拖曳力共同作 用下， 造就了中国大陆及邻区的构造应力场． 各主要边界带的作用力类型和应力水平见表 2．

表 2 模拟得到的主要边界段的特征作用力及相应的应力水平Table 2 Simulated characteristic forces acting on main boundaries and related stress levels

1) 软流层静压推力是广义的洋脊推力， 实质是由板块自重、 岩石层底部深度和地形等因素造成的作用力， 在海洋与大陆板块系统中处于动力平衡状态， 不驱动板块运动(魏东平， 2000)． 只考虑软流层静压推力的模拟结果与实际差别较大， 表明软流层静压推力并非影响中国大陆及邻区构造应力场的控制性因素．

2) 板块边界力对中国大陆及邻区应力场起主导作用， 印度板块对欧亚板块的碰撞挤压主导了中国大陆地区应力场的基本形态； 碰撞造成了川滇地区出现大量走滑型地震， 这与应力受到四川盆地的阻挡发生向南偏转有很大关系； 碰撞甚至影响到阿穆尔板块地区的应力场．

3) 琉球海沟—南海海槽俯冲带自南向北， 边界力发生了挤压-张性的变化， 以张性为主， 显示明显的分段特性. 该边界带作用力主要影响华北和华东地区应力场， 对华南地区影响很小． 阿穆尔板块东边界作用力由南向北逐渐递减， 主要由于其受到太平洋板块俯冲挤压作用减小所致， 影响了东北和华北地区的应力场．

4) 扩张边界的演化常常会改变岩石层的三维结构， 造成岩石层减薄、 地幔流上升等现象(江为为等, 2002； Kubo,Fukuyama, 2003). 而一些板块动力学模拟计算中(Wei, 1997; Zhao, Müller, 2003)将扩张边界视为自由边界， 则过于简化， 不能较好反映真实的作用力． 本研究考虑了贝加尔裂谷、 冲绳海槽的拉张作用， 取得了较好的效果. 结果表明， 贝加尔裂谷影响了内蒙古东部和东北地区的应力场．

5) 中国大陆剪应力分布与该地区地震辐射能量密度分布在大部分区域内有较好的空间正相关性. 但华南和塔里木等盆地为异常区域， 主要由于这些地区位于稳定的古板块构造中， 地震发生较少．

本文研究模型在一定程度上对板块边界力做出了合理估计， 并对中国大陆地区的构造应力场， 尤其是一阶分量给出解释， 但仍存在一些不足． 首先模型的岩石层性质只考虑为弹性， 未考虑模型内部断裂带的影响， 这与中国大陆实际复杂的环境不符合. 同时岩石层厚度设定过于单一， 无法反映局部小尺度的应力场特征， 这也就导致了川滇地区应力型与实际符合， 而应力方向却与实际差别较大以及华北和华中地区应力型与实际不符的结果． 冯向东等(2005)单独分析了大华北地区， 并考虑主要断层带影响， 取得了较好效果． 若考虑以上因素影响， 需要进一步的研究.

参考文献
[1] 陈颙.1988. 地壳岩石的力学性能: 理论基础与试验方法[M]. 北京: 地震出版社: 14-16.(1)
[2] 邓起东, 张裕明, 徐桂林, 范福田. 1979. 中国构造应力场特征及其与板块运动的关系[J]. 地震地质, 1(1): 11-22.(1)
[3] 冯向东, 魏东平, 陈棋福. 2005. 基于观测应力场的大华北地区动力学机制探讨[J]. 地震学报, 27(1): 110-119.(1)
[4] 冯向东, 刘卫华, 魏东平. 2006. 印度—澳大利亚板块应力场数值模拟及其动力学分析[J]. 地质力学学报, 12(4): 429-440.(1)
[5] 傅容珊, 黄建华, 徐耀民, 李力刚, 常筱华. 2000. 印度与欧亚板块碰撞的数值模拟和现代中国大陆形变[J]. 地震学报, 22(1): 1-7.(1)
[6] 皇甫海泉, 魏东平. 2010. 冲绳板块应力场数值模拟及其动力学分析[J]. 地震学报, 32(5): 529-543.(1)
[7] 江为为, 刘少华, 郝天珧, 宋海斌. 2002. 应用重力资料估算东海冲绳海槽地壳厚度[J]. 地球物理学进展, 17(1): 35-41.(1)
[8] 焦明若, 张国民, 车时, 刘杰. 1999. 中国大陆及其周边地区构造应力场的数值计算及其在地震活动性解释上的应用[J]. 地震学报, 21(2): 123-132.(1)
[9] 李廷栋. 2010. 中国岩石圈的基本特征[J]. 地学前缘, 17(3): 1-13. (1)
[10] 梁海华, 侯建军, 刘树文, 常祖峰, 阎志梅. 1999. 中国构造应力场与大震复发周期关系的数值模拟[J]. 地震地质, 21(1): 51-57.(1)
[11] 汪素云, 陈培善. 1980. 中国及邻区现代构造应力场的数值模拟[J]. 地球物理学报, 23(1): 35-45.(1)
[12] 汪素云, 许忠淮, 俞言样, 张琳. 1996. 中国及其邻区周围板块作用力的研究[J]. 地球物理学报, 39(6): 764-771.(1)
[13] 魏东平. 2000. 软流层静压推力及其板块动力学意义[J]. 地质力学学报, 6(1): 4-14.(3)
[14] 魏东平, 左如斌, 濑野徹三. 2001. 全球观测应力场应力取向数据的加权统计分析[J]. 中国科学技术大学学报, 31(1): 50-56.(1)
[15] 谢富仁, 崔效锋, 赵建涛, 陈群策, 李宏. 2004. 中国大陆及邻区现代构造应力场分区[J]. 地球物理学报, 47(4): 654-662.(1)
[16] 徐纪人, 赵志新. 2006. 中国岩石圈应力场与构造运动区域特征[J]. 中国地质, 33(4): 782-792.(1)
[17] 徐纪人, 赵志新, 石川有三. 2008. 中国大陆地壳应力场与构造运动区域特征研究[J]. 地球物理学报, 51(3): 770-781.(1)
[18] 许忠淮, 汪素云, 俞言祥. 1992. 根据观测的应力方向利用有限单元方法反演板块边界作用力[J]. 地震学报, 14(4): 446-455.(1)
[19] 许忠淮, 吴少武. 1997. 南黄海和东海地区现代构造应力场特征的研究[J]. 地球物理学报, 40(6): 773-781.(1)
[20] 胥颐, 刘建华, 郝天珧, 刘劲松, 李志伟. 2006. 中国东部海域及邻区岩石层地幔的P波速度结构与构造分析[J]. 地球物理学报, 49(4): 1053-1061.(1)
[21] 臧绍先, 宁杰远. 2002. 菲律宾海板块与欧亚板块的相互作用及其对东亚构造运动的影响[J]. 地球物理学报, 45(2): 188-197.(1)
[22] 章纯. 2007. 中国东部地区地震活动与构造应力场关系的有限元数值模拟[J]. 西北地震学报, 29(3): 230-235.(1)
[23] 张东宁, 许忠淮. 1994. 青藏高原现代构造应力状态及构造运动的三维弹粘性数值模拟[J]. 中国地震, 10(2): 136-143.(1)
[24] 张东宁, 许忠淮. 1999. 中国大陆岩石层动力学数值模型的边界条件[J]. 地震学报, 21(2): 133-139.(1)
[25] 朱守彪, 石耀霖. 2006. 中国大陆及邻区构造应力场成因的研究[J]. 中国科学: D辑, 36(12): 1077-1083.(2)
[26] An M, Shi Y. 2006. Lithospheric thickness of the Chinese continent[J]. Phys Earth Planet Interi, 159(3-4): 257-266.(1)
[27] Bassin C, Laske G, Masters G. 2000. The current limits of resolution for surface wave tomography in North America[J]. EOS Trans AGU, 81: F897.(1)
[28] Bird P. 2003. An updated digital model of plate boundaries[J]. Geochem Geophys Geosyst, 4(3): 1027, doi:10.1029/2001GC000252.(2)
[29] Bird P, Liu Z, Rucker W K. 2008. Stresses that drive the plates from below: Definitions, computational path, model optimization, and error analysis[J]. J Geophys Res, 113(B11406): 32, doi:10.1029/2007JB005460.(2)
[30] Choy G L, Boatwright J L. 1995. Global patterns of radiated seismic energy and apparent stress[J]. J Geophys Res, 100(B9): 18205-18218, 18228.(1)
[31] Cloetingh S, Wortel R. 1986. Stress in the Indo-Australian plate[J]. Tectonophysics, 132(1-3): 49-67.(3)
[32] Coblentz D D, Richardson R M. 1995. Statistical trends in the intraplate stress field[J]. J Geophys Res, 100(B10): 20245-20220, 20255.(0)
[33] Coblentz D D, Richardson R M. 1996. Analysis of the South American intraplate stress field[J]. J Geophys Res, 101(B4): 8643-8657.(3)
[34] Coblentz D D, Zhou S, Hills R R, Richardson R M, Sandiford M. 1998. Topography, boundary forces and the Indo-Australian intraplate stress field[J]. J Geophys Res, 103( B1): 919-932.(1)
[35] Conrad C, Lithgow-Bertelloni C. 2006. Influence of continental roots and asthenosphere on plate-mantle coupling[J]. Geophys Res Lett, 33(5).(1)
[36] Copley A, Avouac J P, Royer J Y. 2010. India-Asia collision and the Cenozoic slowdown of the Indian plate: Implications for the forces driving plate motions[J]. J Geophys Res B: Solid Earth, 115(B03410): 14, doi.10.1029/2009JB006634.(1)
[37] DeMets C, Gordon R G, Argus D F. 2010. Geologically current plate motions[J]. Geophys J Int, 181(1): 1-80.(1)
[38] Flesch L M, Haines A J, William E H. 2001. Dynamics of the India-Eurasia collision zone[J]. J Geophys Res, 106(B8): 16435-16460.(1)
[39] Forsyth D, Uyeda S. 1975. On the relative importance of the driving forces of plate motion[J]. Geophys J R astr Soc, 43(1): 163-200.(1)
[40] Gripp A E, Gordon R G. 1990. Current plate velocities relative to the hotspots incorporating the NUVEL 1 global plate motion model[J]. Geophys Res Lett, 17(8): 1109-1112.(1)
[41] Gripp A E, Gordon R G. 2002. Young tracks of hotspots and current plate velocities[J]. Geophys J Int, 150(2): 321-361.(1)
[42] Gutenberg B, Richter C F. 1956. Earthquake magnitude, intensity, energy, and acceleration (Second paper) [J]. Bull Seism Soc Amer, 46(2): 105-145.(1)
[43] Hanks T C, Kanamori H. 1979. A moment magnitude scale[J]. J Geophys Res, 84(B5): 2348-2350.(1)
[44] Kincaid C, Silver P. 1996. The role of viscous dissipation in the orogenic process[J]. Earth Planet Sci Lett, 142(3-4): 271-288.(1)
[45] Kubo A, Fukuyama E. 2003. Stress field along the Ryukyu Arc and the Okinawa Trough inferred from moment tensors of shallow earthquakes[J]. Earth Planet Sci Lett, 210: 305-316.(1)
[46] Li S, Mooney W D, Fan J. 2006. Crustal structure of mainland China from deep seismic sounding data[J]. Tectonophysics, 420(1-2): 239-252.(1)
[47] Lithgow-Bertelloni C, Guynn J H. 2004. Origin of the lithospheric stress field[J]. J Geophys Res, 109(B01408): 32, doi:10.1029/2003JB0062467.(1)
[48] McKenzie D. 1967. Some remarks on heat flow and gravity anomalies[J]. J Geopys Res, 72(24): 6261-6273.(1)
[49] Müller R D, Sdrolias M, Gaina C, Roest W R. 2008. Age, spreading rates, and spreading asymmetry of the world’s ocean crust[J]. Geochem Geophys Geosyst, 9(Q04006), doi:10.1029/2007GC001743.(1)
[50] Pacanovsky K M, Davis D M, Richardson R M, Coblentz D. 1999. Intraplate stresses and plate-driving forces in the Philippine Sea Plate[J]. J Geophys Res, 104(B1): 1095-1110.(1)
[51] Reynolds S D, Coblentz D D, Richard R H. 2002. Tectonic forces controlling the regional intraplate stress field in continental Australia: Results from new finite element modeling[J]. J Geophys Res, 107(2131): 15, doi:10.1029/2001JB000408.(1)
[52] Richardson R M, Solomon S C, Sleep N H. 1979. Tectonic stress in the plates[J]. Rev Geophys, 17(5): 981-1019.(2)
[53] Saito A. 1995. Differential Stress Fields in the Nazca and South American Plates: Two Dimensional Modeling Taking into Account the Plate and Crustal Structure Variations and Plate Driving Forces[D]. Tokyo: Earthquake Research Institute of University of Tokyo: 1-52.(1)
[54] Stein C A, Stein S. 1992. A model for the global variation in oceanic depth and heat flow with lithospheric age[J]. Nature, 359: 123-129.(1)
[55] Wei D P. 1997. Pseudo-3-D Spherical Modeling of the Intraplate Sresses of the Eurasian Plate: Implications to Plate Dynamics[D]. Tokyo: Earthquake Research Institute of University of Tokyo: 1-45.(8)
[56] Zhao S, Müller R D. 2003. Three-dimensional finite-element modeling of the tectonic stress field in continental Australia[J]. Geol Soc of Aus Spec Publ, 372: 71-89. (1)
[57] Zoback M L, Zoback M D, Adams J, Assumpcao M, Bell S, Bergman E A, Blumling P, Brereton N R, Denham D, Ding J, Fuchs K, Gay N, Gregersen S, Gupta H K, Gvishiani A, Jacob K, Klein R, Knoll P, Magee M, Mercier J L, Muller B C, Paquin C, Rajendran K, Stephansson O, Suarez G, Suter M, Udias A, Xu Z H, Zhizhin M. 1989. Global patterns of tectonic stress[J]. Nature, 341: 291-298.(2)
[58] Zoback M L, Magee M. 1991. Stress magnitudes in the crust: Constraints from stress orientation and relative magnitude data[J]. Phil Trans R Soc Lond A: Physical and Engineering Sciences, 337(1645): 181-194.(1)
[59] Zoback M L. 1992. First- and second-order patterns of stress in the lithosphere: The World Stress Map Project[J]. J Geophys Res, 97(B8): 11703-11728.(3)