20世纪70年代Tapponnier等^[1]和Molnar等^[2]提出的青藏高原侧向挤出滑移运动图像及力学机制，对认识亚欧地壳运动及我国西南地区地壳块体的相互运动很有意义。鄢家全等^[3]、邓起东等^[4]阐述了中国构造应力场的总体特征。罗灼礼^[5]分析了构造应力场和地震活动的关系。阚荣举等^[6]讨论了青藏高原东南缘现代构造应力场分区及构造活动特征。许忠淮等^[7]推断青藏高原东南缘平均主应力轴的方向。钟继茂等^[8]根据川滇多个地震区的256次中强地震震源机制解，给出各区应力张量的定量分析结果。

地球物理探测揭示的强震孕育的发震结构和物性条件需要与介质的力学性质相联系^[9]。现代构造应力研究是目前地震学和地球动力学研究的热点，但视野过于局限，只有全面了解区域应力的构造演化过程，才能充分了解地震孕育机制，预判应力发展趋势并分析未来地震危险性。因此，古应力场的研究具有重大的理论和实践意义，也是本文探讨的主要内容。

古构造应力场定量化研究方法很多，但受适用条件等因素限制，效果均不佳。本文采用有限元方法反演应力场。地质体的变形具有变形机制复杂、变形经历时间长和变形大等特点，到目前还没有建立一个能精确描述其复杂变形过程的本构方程，但在一定限制条件(如地壳浅层的短期变形)下, 地质体大体符合弹性(线性或非线性)变形特征。

本文介绍有限元数值模拟方法、数据、建模过程，详细描述后处理中的优化分析过程，给出9种不同的应力边界条件方案。对每种方案进行优化分析后处理，并将相应台站形变模拟值与GPS实测值进行误差分析，选取最佳方案作为青藏高原东南缘古构造应力场。

1 有限元数值模拟

ABAQUS是具有齐全的单元模式、丰富的材料模型以及精确的分析过程的非线性有限元分析软件，特别适用于岩土工程方面的模型分析。

创建有限元模型一般需要建立各个部件，并赋予物理材料性质，然后将其装配形成模型整体。此流程操作时存在一定弊端：2个部件相结合的边界不易吻合，导致后续未能严丝合缝地装配。为此，ABAQUS提供了草图功能，用户可将整个模型内外部所有边界预先画入草图，由于软件不检查模型的自封闭性，后面建立部件时可直接导入草图。

对青藏高原东南缘进行有限元建模前需要先弄清该区的主要断裂和块体分区。

1.1 主要断裂和块体分区

青藏高原东南缘发育多条大型断裂带，围绕菱形块体东侧的是玉树-鲜水河-小江断裂带，包括甘孜-玉树断裂、鲜水河断裂、安宁河断裂、则木河断裂以及小江断裂，属于一级构造边界带；菱形块体西侧边界由澜沧江断裂和红河断裂带构成；同时，内部发育的小金河-丽江断裂将菱形块体一分为二，周围的大断裂则主要有龙门山断裂带和南汀河断裂带。本文将龙门山断裂带北侧的岷江断裂也考虑进建模过程中，因为从速度场结果来看，速度值在岷江断裂附近明显减小，说明两侧地块的物性存在较大差别，应区分对待。

图 1标注了青藏高原东南缘的主要活动断裂带，这些断裂将整个地区划分为5个一级块体：马尔康块体、川滇菱形块体、扬子块体、滇西南块体和拉萨块体。其中，岷江断裂又将马尔康块体划分为松潘-甘孜块体和秦岭块体2个次级地块，小金河-丽江断裂则将菱形块体分为北侧的滇西北块体和南侧的滇中块体^[10]。

对照图 1，选取分区的主要控制点，要求既能反映模型主要特征，又不可过于简略以致与实际不符。连接各控制点完成草图，之后创建部件时导入并逐个建立块体，同时将垂向拉伸长度设置为20 km，即上地壳平均厚度(模型认为断裂带主要分布于上地壳)。模型中除龙门山断裂带以外均设置为直立断层。考虑龙门山断裂带浅陡深缓特征，上地壳浅部倾角约60°~80°，深部则放缓至30°~40°，故设置龙门山断裂倾角为70°^[11]。

1.2 介质参数

建模时选用线弹性模型，主要涉及的物理量为杨氏模量和泊松比。设置断裂带则需要对每个块体进行介质分区，沿断裂带在块体内部作宽度10 km的平行线，将这10 km区域作为弱化带处理，其杨氏模量为所在块体杨氏模量的1/3。各块体介质参数如表 1所示^{[9, 12-13]}，其中滇中块体杨氏模量较其他块体偏大，原因可能为滇中块体与扬子块体本属同一块体。

1.3 网格划分

ABAQUS划分网格有多种策略，首先需要在部件上建立种子来控制网格节点的密度。本文选择全局种子尺寸为30 km，断裂带附近由于精度要求高，故将断裂附近边界种子尺寸适当缩小。选择自由网格划分技术，通过四面体单元来划分网格，共建立61 833个网格、31 822个节点。

1.4 分析步及相互作用

一般来说可建立3个分析步：首先是系统默认的初始步，在该步中可添加模型的初始应力场，以及一直存在的载荷或约束；通常会在第2个分析步建立接触，使模型能够在接触稳定后过渡至后面的步骤，否则可能会出现一些不稳定或是不收敛的解；最后加载边界条件，此时将初始时间增量定为总时间的1/10，从而使载荷平滑地加载到模型，有效减小计算不收敛发生的概率。

定义相互作用时首先选定作用属性为“接触”，接下来添加切向行为，即沿接触面的行为，选用罚函数摩擦模型，并将摩擦系数设置为0.6^[14]。选择模型的主从接触面共9对，勾选“仅调整至消除过盈量”，从而使得主从接触面在分析开始时能够精确接触。本文以4 Ma时瞬时平衡态的位移场和应力场图像作为计算结果进行分析。

1.5 应力边界条件

中国大陆现代构造应力场主要起源于中国大陆周围板块间的相互作用。印度板块NNE向的强烈碰撞作用是中国大陆构造应力场最主要的动力源，它控制了中国大陆构造应力场和构造运动的总体特征。对青藏高原东南缘来说，西北部受到青藏高原隆起及其向东物质挤出的作用；西南部受印度板块的挤压作用；东部受到华南块体的阻挡作用。其中，稳定的四川盆地的阻挡作用导致物质东流的方向发生偏转，从向东转为向东南。青藏高原东南缘的应力场大致存在3个分区：27°N以北的区域应力场方向为近EW向；27°N以南、100°E以西的区域应力场方向为NE-SW向；27°N以南、100°E以东的区域应力场方向为NW-SE向。动力边界条件的具体数据参考前人研究。陈化然等^[15]根据震源机制解给出具体数据为：F₁=105 MPa，F₂=70 MPa，F₃=33 MPa，F₄=66 MPa，F₅=100 MPa。本文据此对5个参数进行调节，确定最佳边界条件作为青藏高原东南缘的古构造应力。模型上表面为自由表面，底部垂直方向固定，水平方向自由。

2 模拟结果后处理

将场输出请求(field output request)设置为我们关心的物理量，创建相应模型的Job并提交给软件进行计算，即可得到模拟结果。本次模拟运行内存为8 G，计算时间约10 min。运行完后，使用Visualization功能即可将计算结果进行可视化处理。

由于本文需要从形变模拟结果中提取出位于研究区域内81个GPS测站所对应的位移模拟值进行匹配，故选用Python语言进行场输出数据的后处理。ABAQUS有限元程序通过集成Python脚本向二次开发用户提供丰富的库函数，可直接操纵ABAQUS内核，实现建模、划分网格、指定材料属性、提交作业、后处理分析等功能。

后处理的第一步是建立模型结果文件odb，然后从odb对象中读取step的最后一步，并从中得到最后一帧(frames[-1])的fieldOutputs场输出结果数据。该数据中包含创建场输出时勾选的所有变量结果，这里选择“U”，即形变结果。最后循环读取每个实例中每个节点的东向形变U₁和北向形变U₂。需要注意的是，ABAQUS是按照块体对节点单独编号的，而不是整体编号，所以需要清楚每个块体有多少个节点，GPS测站所在位置位于哪个实例块体的哪个节点编号上，这样才能准确读取对应台站的模拟值并进行后续匹配。本文模型中，松潘-甘孜块体共10 732个节点，其中有12个对应GPS台站；秦岭块体共5 724个节点，其中有3个对应GPS台站；滇西北块体共4 624个节点，其中有4个对应GPS台站；滇中块体共6 363个节点，其中有11个对应GPS台站；扬子块体共17 437个节点，其中有22个对应GPS台站；滇西南块体共6 789个节点，其中有9个对应GPS台站；拉萨块体共10 164个节点，其中有10个对应GPS台站。利用经纬度换算确定每个GPS台站对应的InstanceName和NodeID，然后使用Python查找这2个条件下对应点的位移值U₁、U₂，保存在文件中，方便后续作误差分析。

3 边界条件误差分析

前面提到，本文对5个参数F₁~F₅进行调节，确定最佳边界条件作为青藏高原东南缘的古构造应力。调节的原则是增大或减小10 MPa，由此得到以下9种应力边界条件选择方案。对于每一种方案，采用上面的方法提取出对应台站的模拟值abaqus_U_1i和abaqus_U_2i(i=1, 2, …70)，对应GPS测站的真实值为gps_U_1i和gps_U_2i(i=1, 2, …70)，由此可以计算形变大小模拟值abaqus_U_i和实测值gps_U_i，形变方向模拟值abaqus_U_azi和实测值gps_U_azi，利用最小二乘原理计算二者之间的误差，公式如下：

$ \delta U = \\ \sqrt {{{\left( {{\rm{abaqus}}\_{U_1}-{\rm{gps}}\_{U_2}} \right)}^2} + {{\left( {{\rm{abaqus}}\_{U_2}-{\rm{gps}}\_{U_2}} \right)}^2} + \cdots }, \\ \;\left( {i = 1, 2, \cdots, 70} \right) $

(1)

$ \delta {U_{{\rm{az}}}} =\\ \sqrt {{{\left( {{\rm{abaqus}}\_{U_{{\rm{a}}{{\rm{z}}_1}}}-{\rm{gps}}\_{U_{{\rm{a}}{{\rm{z}}_2}}}} \right)}^2} + {{\left( {{\rm{abaqus}}\_{U_{{\rm{a}}{{\rm{z}}_2}}}-{\rm{gps}}\_{U_{{\rm{a}}{{\rm{z}}_2}}}} \right)}^2} + \cdots },\\ \;\left( {i = 1, 2, \cdots, 70} \right) $

(2)

4 讨论

由表 2可知，方案4给出的应力边界条件得到的模拟结果与GPS实测值最为匹配，形变大小误差为1.432 8 MPa，形变方向误差为16.372 1°。图 4展示了该边界条件下模拟得到的位移场以及与GPS台站实测值的对比图，图 4(a)中箭头方向代表位移矢量方向，长短以及颜色均代表位移大小；图 4(b)为各基准站模拟结果与GPS实测值的对比，黑色箭头为GPS实测值，红色为模拟值。对比二者可以看出，方案4的数值模拟结果较好地反映了青藏高原东南缘西北部青藏高原推挤、东侧扬子块体阻挡的格局，同时位移方向在菱形块体及拉萨块体内部绕东构造结顺时针旋转的特性。图 4(b)中大部分站点在大小和方向上吻合较好，仅北侧松潘-甘孜块体内部以及南侧滇西南块体内部少数存在一定差异，其原因可能是这些地块内部存在某些特殊构造^[14]。

ABAQUS中绝对主应力最大值的正值代表主张应力，负值代表主压应力。观察青藏高原东南缘最大主应力方向，同样在“Y”字形结构下方呈现高压区；拉萨块体东侧与小金河-丽江断裂带西段交界区域存在一定东西向压力，通过澜沧江断裂带在北侧形成近南北向张力。滇西南块体整体呈现拉张趋势，与拉萨块体交界处的南汀河断裂带西段北侧承受着较大的压力，在通过断裂带挤压过程中转化为张力；红河断裂带中段存在来自滇中块体的南向压力，虽然数值较小，但仍然不容忽视。图 5(b)给出了利用时空应力场反演(spacial and temporal stress inversion, SATSI)程序以及青藏高原东南缘1 131个震源机制解得到的该区现今应力图像。可以看出，青藏高原东南缘西北部受到来自青藏高原NEE向挤压作用；东部则有扬子块体阻挡产生的NWW、NW向应力，该阻挡作用导致来自高原的E向应力在菱形块体及其西侧发生顺时针旋转，菱形块体内部呈现北部SE向、南部SSE向应力分布的情况；另一方面，西南部由于印度板块持续推挤产生NE向应力。

5 结语

分析认为，青藏高原东南缘4 Ma BP的古应力场主要起源于中国大陆周围板块的相互作用，特别是印度板块NNE向强烈碰撞作用，成为中国大陆尤其是西南部青藏高原地区构造应力场最主要的动力来源，控制各个块体相互作用的方式和运动格局。青藏高原东南缘古应力场主要包括以下几个力源：西北部青藏高原侧向挤压造成的WE向应力约105 MPa；西南部直接来自于印度板块的NE-WS向应力约70 MPa；南部NS向作用力33 MPa；东南部扬子块体侧向NW-SE阻挡力56 MPa；北东部受扬子块体强烈EW向阻挡力90 MPa。这些力源共同作用于青藏高原东南缘，经过数百万年的不断演化和构造运动，最终形成现今的复杂应力场。

本文关于应力边界条件的初始设定来自其他文献，仅分为5个区域，未来可参考更多研究成果同时加入程序计算以提高精确度。同时，调整应力作用时，仅考虑了古应力的大小，未调整方向，未来希望在二次开发过程中不断调整其大小和方向，实现结果的自动匹配和优化分析，从而得出更为精确、可靠的古应力场。