地应力是地震活动、造山事件、沉积盆地形成以及地下流体运动等地质过程的直接控制因素，也是采矿、水利、桥梁等工程建设中必须考虑的地质条件.中国大陆属于欧亚板块的一部分，被夹持在印度板块、菲律宾海板块、太平洋板块和西伯利亚-蒙古块体之间，受印度板块与欧亚板块的碰撞和菲律宾板块的俯冲的影响，其内部构造运动极为活跃，板块构造环境所控制的地应力场是活动构造的直接动力因素.通过活动断裂、原地应力测量和地震活动的震源机制等研究能够反映局部的地应力场状态，并在一定程度上反映区域地应力场的分布趋势特征^[1]，尽管中国大陆地应力场在区域上受控于周缘的板块构造环境，而在局部地区则显示了活动断裂构造和块体构造单元地壳结构的影响，地应力场在不同的空间尺度具有不完全一致的控制因素.活动构造与地应力场的关系相互影响、密不可分，一方面，活动构造的活动性体现了局部应力场的特征；另一方面，活动构造的存在会造成局部地应力场的调整，从而影响地应力场在局部地区的分布状态^[2-5].地形起伏和地壳介质密度分布不均匀引起的重力势所产生的地应力场也是地应力分析中重要的因素之一，并引起许多学者的关注和研究^[6-8].中国大陆在不同区域地壳厚度分布不均匀，内部结构也存在很大的差异^[9-16]，青藏高原的地壳厚达70 多公里，而东部华北地区只有30km 左右；地形地貌起伏很大，青藏高原平均高程达4km.地应力场在不同区域、不同深度其主要控制因素也有所区别.众多学者^[17-19]利用数值模拟方法从总体上分析了中国大陆构造应力场的特征和成因，模拟过程中主要有以下几点考虑不够：（1）模拟过程中都采用平面直角坐标系，而对于中国大陆整体上经纬度跨度都比较大，选用球壳模型更能够接近实际；（2）在模拟过程中，边界条件和模拟结果约束或者只选用GPS，或者只选用实际地应力场观测数据，没有把二者给予统一考虑；（3）重力作为影响地应力场分布的重要因素之一，缺乏对其定量的分析.完整地理解现今地应力场的状态，不仅需要现今的构造环境和岩石圈展布特征，而且需要了解岩石圈在最近地质历史时期内的运动状态和动力学演化过程，特别是活动构造所反映的地应力场状态和演化.本文以中国大陆主要活动断裂所控制的块体模型为基础，采用黏弹性壳体有限单元，综合考虑GPS、地应力实测结果以及重力等影响，分析中国大陆地应力场的特征和控制因素．

从岩石力学上来说，尽管整个地壳的岩石主要表现为弹性，但是目前在中地壳普遍存在一个岩石强度相对弱的区域.地壳上部的构造块体单元可以看作是完全弹性介质，活动断裂作为块体边界具有与构造块体本身完全不同的性质，活动断裂运动除了伴随地震发生瞬间明显的断层滑动外，更多的运动则表现为大量的无震滑动和长期蠕动等形式^[20]，活动断裂带可以视为调整块体弹性应力的软弱带.

现今中国大陆地应力场是在周缘板块构造环境条件下，地应力场长期演化过程中的一个瞬间状态.地壳岩石在大时间尺度上通常可用标准线性体来近似，Maxwell体为标准线性体的一个特例（图 1）^[21-24]，作为黏弹性的岩石具有应力松弛和应变蠕变的特性.其本构关系在Laplace域可表示为

图 1

Fig. 1

图 1 Maxwell标准线性黏弹体（a）及其蠕变（b）和松弛曲线（c) Fig. 1 (a) The axwell form of standard linear solid； (b) The creep compliance;（c）The relaxation modulus

(1)

其中

(2)

τ 为弛豫时间：

(3)

当应变为常数时，即，意味着边界条件近似为恒定位移，应力σ 的变化应符合松弛曲线；当应力为常数时，即意味着边界条件为给定一个持续恒定的应力，应变的变化应与蠕变曲线一致.

中国大陆块体通过板块边界受到相邻板块的作用力，这些作用力一方面通过完全弹性变形作用于中国大陆内部，另一方面，相邻板块持续的作用通过岩石的黏滞性对中国大陆内部作用，这种作用不同于完全弹性，在不同的时间不同块体具有不同的状态.由于受到断层库仑摩擦破裂强度的限制，块体内能够维持的地应力通常符合Byerlee准则^[25].当板块边界作用于大陆地壳内部时，板块边界的作用力不会超过Byerlee准则所确定最大应力，所以边界作用力在一定的地质历史时间内，当板块构造环境没有发生变化时，应为一个持续恒定的应力，这种边界作用力除了瞬时的弹性应力作用于块体内部，还会由于不同块体具有不同岩石黏滞性，导致在不同的时段具有不同的应力状态，并最终达到平衡，大陆地壳整体都处于一个临界的断裂平衡的状态，全球广泛存在有水库蓄水和流体注入引起的地震活动、一个地震的发生会触发另一个地震的发生，原地深孔地应力测量结果符合Byerlee准则，这些现象恰好说明了目前大陆内部，特别是活动构造断裂附近处于临界断裂平衡状态^[26-27].

地应力场通常为构造应力场与重力势引起的应力场的叠加，重力势引起的应力场主要是由地形起伏和岩石层密度不均匀造成的，当地势平缓且地下岩石分布均匀时，重力对地应力场影响可以忽略不计，构造应力场通常由于受到岩石强度的限制，在不同的构造环境具有不同的变化规律.

中国大陆是由不同活动块体组成，块体之间以大型的活动断裂带为分界，板内块体具有不同活动程度的次级断裂^[28-29].中国大陆及其邻区的活动地块主要有拉萨、羌塘、巴颜喀拉、柴达木、祁连、川滇、塔里木、天山、准噶尔、华北等（图 2）.随着现代GPS观测资料的不断积累，初步构建了中国地壳活动的运动模型^[30-33]，并在此基础上分析了其应变场的特征^[34-35]，中国大陆的现今构造变形主要以刚性地块运动为主，如塔里木、鄂尔多斯、华南等地块；只有少数块体内部存在比较明显的连续变形，如青藏高原和天山.块体分界的大型活动断裂区是变形的主要区域，边界断裂区的应变率远大于块体内部的应变率，通常边界断裂带都具有数十公里的宽度，发育了中国大陆的主要地震活动^{[29， 36]}.周缘的印度板块、菲律宾板块和太平洋板块的运动是中国大陆板块内部构造活动的主要动力来源.印度板块与欧亚板块为陆陆碰撞，形成挤压环境；菲律宾板块和太平洋板块与欧亚板块为洋陆碰撞，太平洋洋壳和菲律宾的洋壳向欧亚板块之下俯冲，形成了日本海和冲绳海槽，形成拉张的边界环境.

图 2

Fig. 2

图 2 中国大陆构造块体分区及板块构造环境 Fig. 2 China mainland tectonic blocks and its surrounding platete ctonics

中国及邻区地壳厚度分布的特点为从东到西逐渐增厚，在构造块体内部地壳厚度通常变化较缓，如蒙古高原、华北、华中、华南、青藏高原内部等地区（图 2）.西部青藏高原的地壳厚度平均为60～70km，局部达到74km，而东部地区的厚度平均为30～36km，沿海地区仅有28～30km；地壳局部厚度变化与构造环境密切相关，大型地块的边界带处的地壳厚度梯度也比较大，南北构造带以东地壳厚度为40km，而在南北构造带以西达到60km，在很短的距离上地壳厚度发生了很大的变化.中国大陆的地壳在沉降带或盆地发育地方的地壳较周围薄一些，如塔里木、准格尔、柴达木、鄂尔多斯等盆地，而抬升区或造山带处的地壳较周围厚一些，如青藏高原、天山、太行山等地区^[37-39].稳定的克拉通塔里木、华南、鄂尔多斯地块处，其地壳内低速层不发育，而其它地区都普遍有地壳低速层发育，特别在青藏高原的中地壳低速层整体上普遍存在^[9].而细致地壳的流变结构成果^[40]为本文进行地应力场的数值分析提供了直接的参数依据.

在模拟过程中，以中国大陆块体模型为基础，采用三维球壳黏弹线性体单元，对块体模型进行无结构网格剖分（图 3），模型的地形分布参照GTOPO₃₀数字高程模型，按照区域板块构造环境，设定边界条件和加载.块体边界的大规模活动断裂带的网格单元划分相对细密，从而保证边界断裂带的数值计算的稳定性.各个块体和断裂带分层岩石物性参数见表 1，各个块体分为上地壳、中地壳和下地壳三层，表中的杨氏模量和泊松比为瞬态泊松比，剪切模量系数是指岩石剪切模量在经过足够长的时间后，其长期剪切松弛模量与瞬时剪切模量的近似比值，长时间的泊松比近似为0.5，岩石密度上地壳取为2750kg/m³，中下地壳为3300kg/m³.

图 3

Fig. 3

图 3 有限元网格模型及边界条件 Fig. 3 FEM mesh model and boundary condition

表 1(Table 1)

表 1 模型岩石物性参数 Table 1 Physical property parameters of rocks in model

表 1 模型岩石物性参数 Table 1 Physical property parameters of rocks in model

模型北边界中段为南北向和东西向均为约束边界；北边界东西两端和南边界东段均为南北向约束，东西向自由；西边界北段和东边界南端均为东西向约束，南北向自由；其他边界均按照构造应力环境加载应力.

根据Anderson断层理论可以确定不同性质的断层所在区域的构造应力特征^{[27， 41]}.在正断层、走滑断层和逆冲断层活动的不同区域，其应力特征分别为：

(4)

(5)

(6)

其中σ₁ 为最大主应力，σ₃ 为最小主应力，S_Hmax 为水平最大主应力，S_hmin 为水平最小主应力，S_V 为垂向应力，P_f 为孔隙流体压力，μ 为断层临界摩擦系数，根据Byerlee 实验数据^[25]和深部地应力测量结果^[27]，其取值范围在0.6～1.0之间.由此可得出正断层区域，垂直断层的方向为最小主应力方向，也是水平最小主应力，其大小约为

(7)

在逆冲断层区域，垂直断层的方向为最大主应力方向，也是水平最大主应力，其大小约为

(8)

通常岩石流体密度比ρ_rock/ρ_liquid约为3.0，

当μ=0.6时，S_hmin ≥0.55ρ_rockgh，S_Hmax ≤2.41ρ_rockgh，当μ=1.0时，S_hmin ≥0.45ρ_rockgh，S_Hmax ≤4.22ρ_rockgh.

由于板块边界的最大主应力方向与板块边界近似垂直，所以按照板块构造环境对边界加载最大主应力或最小主应力的应力载荷，除去重力作用的影响，边界条件见表 2.

表 2(Table 2)

表 2 边界加载应力表 Table 2 The load on the boundary

表 2 边界加载应力表 Table 2 The load on the boundary

图 4为在边界加载压力约10万年后的近地表的速度分布，在这一时间点上，模拟速度整体上能够很好地与现今实测的GPS 速度场吻合，因此，现今构造应力场的分布特征可参考该时点的应力状态模拟结果.图 5为在只加载边界应力，而不考虑重力作用，约10万年后的距地表1000m 深度的构造应力场模拟结果，总体上构造应力场在西部为挤压，东部为拉张，构造应力场最大主应力的方向受到板块构造环境和活动构造分布的控制，印度板块与欧亚板块的持续碰撞形成了青藏高原及其周缘的挤压性质的构造应力场，而东部菲律宾板块与太平洋板块的俯冲形成了黄海、东海和环渤海区域的拉张性质构造应力场，中间为拉张环境和挤压环境的过渡区域.

图 4

Fig. 4

图 4 模拟速度与现今实测GPS速度对比图 图中白色的箭头为实测的GPS速度场，参照系为稳定欧亚大陆板块（据王琪等[42]和中国地壳运动观测网络）,黑色箭头为数值模拟结果. Fig. 4 Comparison between the modeling velocity and the measured velocity with GPS The white arrows are GPS measured velocity，the reference frame is stable Europe-Asia continental plate (reference Wang[42]

图 5

Fig. 5

图 5 中国大陆构造应力场水平最大主应力模拟结果分布图（黏弹性模拟，1000 m深度） 其中红色箭头为水平最大主应力的方向，云图为水平最大 主应力值，h为深度，红色为张应力，蓝色为压应力 Fig. 5 Modeling tectonic horizontal maximum principal stress of China mainland with viscoelastic model at 1000 meters depth Red line arrows are the direction of maximum principal stress, cloud image is the valueof maximum principal stress,h is depth, red zone is extensive, blue zone is compressive.

图 6为中国大陆在1000m深度重力势引起的偏应力模拟结果，由于中国大陆地势西高东低，整体上形成三级阶梯，第一级最高的阶梯为青藏高原，平均海拔在4000 m 以上；第二级阶梯位于青藏高原与大兴安岭─太行山─ 巫山─ 雪峰山之间，包括内蒙古高原、黄土高原、云贵高原和塔里木盆地、准噶尔盆地、四川盆地等地区，海拔一般为1000～2000m，惟四川盆地较低，海拔在500 m 以下；第三级最低阶梯，在大兴安岭─太行山─巫山─雪峰山以东，自北而南，有海拔200 m 以下的东北平原、华北平原和长江中下游平原.重力势引起的应力场与构造应力场的结果大致相反，在青藏高原的藏东南、四川西部和云南西北地区为拉张性质的应力场，而在塔里木盆地和东部的平原区为挤压性质的应力场，应力场的方向受地势展布特征控制.图 7 为1000 m 深度位置的中国大陆地应力场，是由构造应力的模拟结果与重力作用模拟结果进行张量求和合成，可以作为浅部地应力场的定量参考模型，地壳深部的地应力场应与构造应力场相近.模拟地应力场与实测的地应力分布特征^[43]（图 8）也比较一致，水平最大主应力的方向主要受到板块构造环境和活动断裂带的控制，总体上以南北构造带为界，以西以较强的压性构造环境为主，以东为较弱的压性地应力环境，由于受到重力的影响，藏南和滇西北局部地区的浅部存在有张性地应力环境，而在塔里木区域，由于重力势引起的应力场与构造应力场都为挤压性质，因此挤压强度得以增加.图 9 为初期瞬时的模拟结果，可作为完全弹性结果，完全弹性模型的构造应力场模拟结果与实测的地应力场存在一定的偏差，如在龙门山断裂带的北段地区，水平最大主应力的方向总体为北西-南东向，而完全弹性模型的结果为北东-南西向.青藏高原南缘在重力作用下形成了张性应力环境，但是由于重力作用影响有限，在深部的张性环境应与深部物质上涌等局部构造活动有更为密切的关系.

图 6

Fig. 6

图 6 中国大陆重力引起的水平最大主应力模拟结果分布图（1000m 深度） 其中红色箭头为水平最大主应力的方向，云图为水平最大主应力值，h为深度，红色为张应力，蓝色为压应力。 Fig. 6 Modeling maximum principal stress of China mainland induced by gravity（at depth 1000m） Red line arrows are the direction of maximum principal stress, cloud image is the value of maximum principal stress,h is depth, red zone is extensive,blue zone is compressive.

图 7

Fig. 7

图 7 中国大陆地应力场的水平主应力模拟结果分布图（1000m 深度） （a）水平最大主应力，其中红色箭头为水平最大主应力的方向，云图为水平最大主应力值，h为深度，红色为张应力，蓝色为压应力；（b）水平最小主应力，其中红色箭头为水平最小主应力的方向，云图为水平最小主应力值，h为深度，红色为张应力，蓝色为压应力；（c）为剪应力，其中红色箭头为水平最大主应力的方向，云图为剪应力值，h为深度。 Fig. 7 Modeling crustal horizontal stress in China mainland (at depth 1000 m) (a) Horizontal maximum principal stress， Red line arrows are the direction of maximum principal stress，cloud image is the value of maximum principal stress，h is depth，red zone is extensive，blue zone s compressive； (b) Horizontal mimm principal stress，Red line arrows are the direction of minim principal stress,cloud image is the value of minim principal stress,h is depth,red zone is extensive,blue zone is compressive;（c）Shear stress,Red line arrows are the direction of maximum principal stress,cloud image is the value of shear stress,his depth

图 8

Fig. 8

图 8 中国现代构造应力场图（谢富仁等[43],2003) Fig. 8 The measured crust stress of China（reference the Crust Stress Databese of Institute of Crustal Dynamics[43]）

图 9

Fig. 9

图 9 中国大陆构造应力场水平最大主应力模拟结果分布图（完全弹性模拟） 其中红色箭头为水平最大主应力的方向，云图为水平最大主应力值，h为深度，红色为张应力，蓝色为压应力． Fig. 9 Modeling tectonic horizontal maximum principal stress in China mainland with pure elastic model Red line arrows are the direction of maximum principal stress， cloud image is the value of maximum principal stress， h is depth，red zone is extensive，blue zone is compressive.

中国大陆地应力场和现今活动构造的性质总体上主要受到印度洋板块的碰撞挤压和太平洋板块与菲律宾板块俯冲的影响.由于受到华南板块、塔里木克拉通以及鄂尔多斯块体的阻挡，使得印度洋板块碰撞作用集中在青藏高原及其周边，青藏高原主体多为活动断裂的走滑，块体内部的连续变形较主要断裂活动要小很多.华北块体、山西地堑除了受到印度洋板块碰撞的影响之外，菲律宾和太平洋俯冲板块所形成的弧后张性环境的影响不容忽视.

岩石蠕变曲线表明，在给定的应力条件下，岩石蠕变逐渐趋于一个常数，这就意味着随时间的演化，岩石块体整体可能会运动，但是由于其应变率逐渐趋于0，所以岩石块体内部处于一个稳定的无变形状态.一方面，上地壳的变形受控于自身的构造块体结构，在板块构造环境控制下，如果不考虑其底部中下地壳的作用，在仅受板块构造环境的作用，其应力状态和变形特征应接近于完全弹性模型的数值分析结果；另一方面，当柔性较强的中下地壳受到与上地壳同样的板块构造环境的作用时，中下地壳会表现出更强的运动性，从而对上地壳的变形和应力分布产生影响.例如在青藏高原，由于中下地壳向两侧挤出，导致上地壳被中下地壳拖曳，使得上地壳在活动断裂的控制下随中下地壳运动，现今中下地壳对上地壳的变形模式和应力分布特征的影响作用不容忽视.

青藏高原东南部主应力方向出现转折，主要有以下几个因素控制：（1）由于印度板块持续的碰撞，使得青藏高原整体受到塔里木克拉通、华南板块和鄂尔多斯块体的阻挡，使得青藏高原随着活动断裂的滑动向两侧挤出^[44]，在挤出过程中，刚性的上地壳块体构造的运动受边界断裂的控制；（2）由于中下地壳较上地壳有更强的柔性，所以挤出时上地壳的运动受到中下地壳拖曳的影响；（3）印度板块向北运动过程中，通过缅甸的实皆断裂（SagaingFault）剪切拉伸作用于欧亚板块，使得上地壳中青藏高原东南部应力和地表运动方向都发生了转向.

本文基于线性黏弹体模型对现今中国大陆的构造应力场进行了数值模拟分析，其前提假设是中国大陆周边的板块构造环境相对稳定，且维持了足够长的时间，从而能够保证模拟分析过程中边界条件的一致性；大型活动断裂所决定的块体结构保持稳定，从而保证在模拟分析过程中地质模型的物性参数保持一致.该结果主要考虑了板块构造作用和区域性的地貌和地壳结构所产生的重力的影响，在进一步的研究中需要考虑底部地幔上涌和岩石圈热演化等影响因素，深入讨论地应力场控制因素.

整体上，中国大陆构造应力场的分布特征与活动构造的展布密切相关，是在板块构造环境控制下，整个地壳岩石黏弹特性长期演化和活动断裂无震蠕动的结果.而伴随地震发生的断层滑动引起的地应力场的瞬时或短期调整，可参考完全弹性分析的结果.陈群策等（私人通信）最近在广元附近进行的汶川地震后的地应力测量，所测的地应力水平最大主应力方向为北东-南西向，而该地区通常的区域应力场的分布特征为北西-南东向.应力场作为动态演化的过程，特别是断层附近的地应力场，由于受到断层蠕动和滑动的影响，地应力的变化会更为突出.郭启良^[45]和廖椿庭^[46]分别在汶川地震前后和昆仑山地震前后的地应力测量所反映的地应力的变化，也同样说明，活动断裂附近的地应力场受到局部构造的影响，从而体现了不稳定性.随着地应力监测网络的建立，可以有效地利用地应力动态的信息，更好地认识区域构造应力与局部地质构造和岩石物性之间的关系.

感谢匿名评审人提出的宝贵意见！