GPS约束下川滇地区下地壳拖曳作用及断裂活动性有限元模拟

刘昌伟 常祖峰 李春光 曾平 马琳

摘要：提取川滇地区GPS速度场边界值作为约束条件，依据地质资料，利用ANSYS有限元模拟软件建立川滇地区三维有限元模型，将川滇地区划分为4个块体，利用静力学分析模块模拟下地壳拖曳作用及2条主要活动断裂的走滑速率。结果显示：模型四周仅加载GPS速度场边界值，川滇菱形块体位移速度较实测速度偏差明显，不满足此地区地球动力学特征，增加下地壳拖曳作用荷载后速度场模拟结果得到优化，川滇菱形块体位移速度差值减小;研究模拟下地壳拖曳作用产生的剪切力发现，拖曳作用在小金河断裂附近由南东转向正南时才能得到最优模拟结果，暗示其对川滇菱形块体拖曳方向发生了偏转。鲜水河、小江断裂模拟走滑速率分别为85 mm/a，6 mm/a，与实际走滑速率具有较好一致性，优化后可对断裂带闭锁研究提供参考。

关键词：川滇地区;有限元模拟;拖曳作用;断层走滑速率

中图分类号：P315727   文献标识码：A   文章编号：1000-0666（2019）03-0385-08

0 引言

川滇地区位于青藏高原东南缘，处于羌塘、巴彦喀拉以及华南块体交接部位，在印度板块与欧亚板块相互碰撞挤压下其地质构造错综复杂。曾融生和孙为国（1992）研究青藏高原下地壳及地幔时发现，东邻高原的青海、川滇地区上地幔存在低速带，分析认为下地壳物质向东南向流动并堆积于此;相关数值模拟、重力异常以及横波分裂等数据均显示青藏高原东南部下地壳较软，更易于流动，导致上地壳与地幔存在解耦（Royden et al，1997;熊熊等，2001;Flesch et al，2005）;受到NE向扬子板块的阻挡，青藏高原物质向ES方向流动，尽管前人提出的下地壳通道流模型如 Couette流以及Poiseuille流等还具有争议，但川滇地区下地壳作为青藏高原物质EN向流动的通道已基本不是争论的焦点（Klemperer，2006，Beaumont et al，2004）;部分研究者通过给定地壳粘滞系数模拟得出了较为可信的川滇地区下地壳流动速度，认为其较上地壳运动快约10 mm/a（Wang，2007;曹建玲等，2013），朱守彪和石耀霖（2004）、王辉等（2007）等在川滇地区上、下地壳相互作用的研究过程中提出下地壳产生了拖曳作用。前人针对川滇地区数值模拟取得了丰富的研究成果，主要思路分为2类：一类将川滇地区四周GPS位移作为荷载，断裂带设置为软弱带，模拟给出应力分布、断层相互作用等，并与发震规律进行对比分析（王凯英，马瑾，2004;陈连旺等，2008），再根据模拟应力场分布划分出了地震危险区（廖思佩等，2016）;另一类将川滇地区划分为多个块体，分析各个块体内部应变状态与相互作用（蒋锋云等，2013），或将断裂带模拟走滑速率与实际走滑速率进行对比，分析断裂带活动性或闭锁状态（王阎昭等，2008;申重阳等，2002;赵静等，2015;李长军等，2018）。本文结合以上研究思路，利用ANSYS有限元模拟软件与地质资料建立川滇地区三维模型，以近期GPS观测数据作为边界条件，模拟对比下地壳拖曳荷载对模型产生的影响;断裂设置为软弱带，与各块体间摩擦接触，将模拟与实际断裂带走滑距离进行对比，探讨摩擦接触定义下断裂带有限元模拟与真实断裂的可比性。

1 川滇地区构造特征及模型建立

11 主要断裂与块体划分

李玶和汪良谋（1975），阚荣举等（1977）首先提出了川滇菱形块体的概念，以丽江—小金河断裂为界将川滇菱形块体划分为川西北与滇中2个次级块体（王铠元等，1987;李铁明等，2003）。巴彦喀拉块体以东昆仑断裂、龙门山断裂和甘孜—玉树断裂为界，早期也被称为川青块体（韩渭宾，蒋国芳，2003），一般根据龙日坝断裂带和虎牙断裂带将其划分为阿坝次级块体、马尔康次级块体和龙门山次级块体（徐锡伟等，2008，陈长云等，2013）。滇东块体活动性断裂分布较为有限，以大凉山断裂为界划分出大凉山次级块体和川中次级块体（宋方敏等，2002）。滇西南块体划分观点较多，主流的划分方式是以澜沧江—红河断裂为界划分出印支板块和滇缅泰板块（苏有锦，秦嘉政，2001）。综上，本文将川滇地区划分成川滇菱形块体、巴颜喀拉块体、滇东块体和滇西南块体4个一级块体（图1）。将连接一级块体的断裂定义为一级断裂，有甘孜—玉树—鲜水河断裂、龙门山断裂、安宁河—则木河—小江断裂、怒江断裂及红河断裂。次级块体间接触的为次级断裂，有东昆仑断裂、里塘断裂、金沙江断裂、丽江小金河断裂、南华—楚雄—建水断裂、南汀河断裂以及怒江断裂南段。

12 介质参数与网格划分

GPS 观测结果显示，川滇地区上地壳形变满足弹性形变的特征（Shen et al，2005）。本文ANSYS模拟采用线弹性静力学分析模块，取川滇地区平均上地壳厚度为30 km，参考前人对川滇地区地壳泊松比研究，将上述4个一级块体赋予不同的密度、弹性模量以及泊松比，断裂设置为软弱带，其力学性质参考可塑状态粘土（含水），各块体介质参数如表1所示。川滇地区断裂以走滑为主，龙门山推覆构造带产状具有浅部陡深部缓的特点（张竹琪等，2010），因此本文将龙门山断裂倾角简化设置为65°，其余断裂带倾角设置为直立，断裂带倾角与块体为摩擦接触，当滑动速率<1 mm/s且走滑距离小于1 mm时质密石英岩摩擦系数为06（Dieterich，1978;Dieterich，Kilgore，1996），在不影响计算结果的前提下断裂带宽度均简化为2 km，模型中共计12条断裂带、15个摩擦面。本文模型均采用四面体网格划分方法，块体网格尺寸为40 km，断裂带网格尺寸为5 km，断裂带与块体接触位置进行局部网格加密，共划分出3 458个网格，28 567个节点（图2）。

13 边界条件

本文GPS数据来自中国地震局GNSS服务平台提供的基础数据产品 http：//wwwcgpsaccn。中国大陆构造环境监测网络共建设有260个基准站和2 000多个区域站，其中66个基准站位于研究区（96°～107°E，21°～35°N）内。本文采用ANSYS静力学分析模块进行求解，时间在此求解模块中无实际物理意义，故荷载均设置为1年产生的位移量。建模过程中对地形起伏与部分断裂构造进行了简化处理，断裂模型均采用了非线性摩擦接触分析，由于具有大量非线性分析，荷载步设置为5，每次加载20%，使其荷载平滑加载以避免难以收敛的问题。为得到最理想的模拟结果，共设计3个模型。模型Ⅰ在不考虑位移随深度发生变化的基础上，将川滇地区66个GPS基准站以及142个流动站（共计208个）实测的年均速度值（2013—2015年、2015—2017年2期）作为控制点，利用Matlab插值功能求出模型四周边界处的位移速度，依次将位移荷载加载在模型四周（图1b），将模型Ⅰ作为初始参考模型。模型Ⅱ在模型Ⅰ的基础上使川滇菱形块体底部受到下地壳拖曳产生的剪切力作用（图3a），位移方向与川滇菱形块体长轴方向一致（约140°方向），位移大小通过多次模拟试验确定理想值。模型Ⅲ在模型Ⅱ的基础上考虑川滇菱形块体底部受到的剪切力在小金河断裂附近发生由东向南偏转， 偏转角度与地表实测GPS偏转角平均值一致约为30°（图3b），位移值大小同样通过多次模拟實验确定。

2 模拟结果

使用GPS位移作为边界条件加载在模型Ⅰ（图4a）四周，为了进一步讨论其是否满足川滇地区地球动力学特征，需与GPS实测点位移进行对比验证。图4a为GPS实测值（基准站点）与同坐标模拟值对比结果，图中箭头指向表示位移矢量方向。从图中可以看出，巴颜喀拉块体、滇东块体以及滇西南块体模拟值能够反映出西北部青藏高原的挤压以及东部四川盆地阻挡的构造格局，与GPS实测值较为接近，而川滇菱形块体模拟值却与实测值存在较大偏差，其模拟—实测速度差分布显示，滇中次级块体（菱块南部）差值尤为明显。

模型Ⅱ（图4b）在川滇菱形块体底部添加一个与剪切力等效的位移荷载，加载方向为川滇菱形块体长轴方向（140°），荷载大小以川滇菱形块体GPS实测位移的平均值148 mm/a为基数进行调整，经过多次模拟实验，确定最佳位移量为9 mm/a，但滇中次级块体残差并未发生明显改善。模型Ⅲ（图4c）保持底部最佳位移量不变，在小金河断裂附近对荷载进行偏转调整后，川滇菱形块体位移速度差分布如图4左侧橙色点所示，较GPS实测东西向、南北向位移速度差在±5 mm/a以内，无系统性误差并呈现呈正态分布，因此模型Ⅲ为最佳拟合结果。

川滇地区平均震源深度为15 km，截取模型Ⅲ在该深度的等效应力并去除边界效应影响（图5）。由图5可见，巴颜喀拉块体应力主要集中于北部以及东昆仑断裂两端;川滇菱形块体应力分布由WN向ES递减，其中小金河以及鲜水河断裂交汇处出现局部应力集中，滇中次级块体内基本无应力分布;滇西南块体应力主要集中于怒江断裂、南汀河断裂两端，块体南缘速度场分布复杂，局部地区出现应力集中。断裂带走滑速率如图6所示，鲜水河断裂、则木河断裂、小江断裂、龙门山断裂以及红河断裂走滑速率明显。其中鲜水河断裂走滑速率由WN至ES表现出先增大后减小的特点，走滑速率变化范围较大，为11～104 mm/a;则木河断裂、小江断裂走滑速率分布较为均匀，分别为26～56 mm/a，52～65 mm/a;龙门山断裂走滑速度整体偏小并集中于中段，为11～25 mm/a;红河断裂走滑速率为32～55 mm/a，集中于与小江断裂交汇处，其余分段未发生明显走滑。

3 讨论

31 川滇菱形块体下地壳拖曳作用

模型Ⅰ在排除了加载方式与介质参数设置不合理的情况下，偏差依然存在，笔者认为是川滇菱形块体的模拟中缺少位移约束所致。下地壳流动对川滇菱形块体产生的“拖曳作用”可能是缺少的主要约束之一。巴颜喀拉块体、滇东块体以及滇西南块体的模拟值与实测值吻合度较高，“拖曳作用”可能并未在这些块体产生显著影响，其作用范围仅在川滇菱形块体内，支持川滇菱形块体下地壳作为青藏高原东南方向流动通道的观点。“拖曳作用”产生的剪切力方向发生了偏转，如图4c所示，川滇菱形块体的2个次级块体（川西北块体、滇中块体）GPS实测位移方向由SE向偏转至近S向，偏转处大致沿小金河断裂分布。前人研究川滇地区地壳结构时发现鲜水河断裂与小金河断裂的交叉处（深紫色区域）地壳泊松比大于03（图7），认为下地壳存在低速层甚至局部熔融（Xu et al，2007;Wang et al，2010）。在此区域出现2种异常并非巧合，笔者推测软弱的下地壳物质由青藏高原沿ES向挤入川滇菱形块体（红色箭头），受到四川盆地刚性块体的“阻挠”，局部因高温挤压形成低速层或熔融，剩余物质转向继续向南流动（紫色箭头），进而下地壳牵引作用对川滇菱形块体产生的剪切力也随之发生改变（图7）。与模型Ⅰ对比显示，模型Ⅲ在川滇菱形块体底部添加剪切力后，其模拟速度场得到优化，EW向与NS向位移差明显改善，十分接近GPS实测值。其次，模型Ⅲ川滇菱块底部剪切力方向由NE向正S偏转30°后才与GPS实测值拟合最佳，拖曳作用随下地壳流向改变而发生偏转。

32 断层走滑速率对比分析

大地构造应力场与模拟得出的理想应力场主要区别在于断裂带局部形成的“閉锁”，闭锁位置、程度等均不可测，而闭锁势必产生局部应力集中，直接利用模拟得出的应力云图去讨论似乎并不合理。本文以鲜水河断裂、小江断裂为例，对断裂带模型摩擦滑动速率与实际活动断层走滑速率进行对比研究。鲜水河断裂是川滇地区现今强烈活动的一条大规模左旋走滑断裂带，附近地震活动频率高，以惠远寺为界分为NW，SE两段。前人对鲜水河断裂进行了大量研究（表2），其整体走滑速率呈中间快两端慢的特点，为888～973 mm/a;本文模拟结果如图6所示，鲜水河断裂SE段与龙门山断裂交汇处走滑速率为104 mm/a，SE段均值为85 mm/a。小江断裂是云南境内地震活动最为频繁的断裂带之一，主要表现为左旋走滑断裂，毛玉平和韩新民（2003）、施发奇等（2012）以及何宏林等（2002）计算出小江断裂走滑速率约7 mm/a，本文模拟走滑速率均值为6 mm/a，与前人结果对比较为接近。值得注意的是，在对龙门山断裂及红河断裂等多条断裂带的对比中发现大部分模拟走滑速率高于实际走滑速率（未发表数据），笔者认为原因除模型简化造成的误差之外，断裂带局部闭锁也是造成实际与模拟断裂走滑速率不一致的重要因素：若能精细化断裂带模型的介质参数与边界条件，有限元法可以得出更真实有效的模拟结果，为研究断裂带闭锁提供更多手段。

33 地质体仿真模拟的局限

尽管本文利用ANSYS模拟得出了部分结果，但地质体仿真模拟存在的局限主要有3点：①深部地应力数据匮乏。地壳具有漫长的构造演化史，后期的改造均叠加在前期地壳固有的地应力之上，但我们对这部分应力分布及大小知之甚少。岩土、采矿领域在地下1 km内取得了部分地应力成果，但其研究过程已十分困难（赵德安等，2007;汪斌等，2012）。②断裂带延伸情况或相互穿插造成的影响，真实地质环境中断裂往往表现为雁列式、共轭式等不规则分布，局部次生断裂的加密很可能对模型整体力学性质产生影响，大量隐伏断裂的存在使得断裂的延伸情况变得更加复杂，单一摩擦模型很难还原真实断裂构造过程。③由于运算速度的限制，若在模型内考虑下地壳粘性、塑形、蠕变等非线性分析以及地幔对流、壳幔作用等热力学条件，运算时间会成倍增加，在不影响地球动力学特征的前提下恰当优化模型变得必要。

4 结论

本文以川滇地区GPS速度场边界值为约束条件，依据地质资料建立川滇地区三维有限元模型，使用ANSYS静力学分析模块进行模拟，主要结论如下：

（1）川滇地区下地壳物质流动对川滇菱形块体产生“拖曳作用”，其产生的剪切力荷载随下地壳流动发生了偏转，仿真模拟考虑“拖曳作用”的影响后能在一定程度上优化川滇地区位移速度场模拟结果。

（2）鲜水河、小江断裂模拟走滑速率分别为85 mm/a和6 mm/a，与实测走滑速率具有较好一致性，对断裂优化后可为“闭锁”研究提供更多手段。

（3）地质体仿真模拟具有局限性，应该从壳内地应力、断裂带精细分布以及模型效率方面进行改进。

参考文献：

曹建玲，王辉，张晶2013青藏高原下地壳流动方式的数值模拟研究[J].地震，33（4）：55-63

陈长云，任金卫，孟国杰，等2013巴颜喀拉块体东部活动块体的划分、形变特征及构造意义[J].地球物理学报，56（12）：4125-4141

陈连旺，李红，陆远忠2008汶川MS80地震同震效应的三维非线性数值模拟分析[J].国际地震动态，（11）：10

韩渭宾，蒋国芳2003川青块体及其向南东方向运动的新证据[J].地震工程学报，25（2）：175-178

何宏林，池田安隆，宋方敏，等2002小江断裂带第四纪晚期左旋走滑速率及其构造意义[J].地震地质，24（1）：14-26

蒋锋云，朱良玉，王双绪2013川滇地区地壳块体运动特征研究[J].地震研究，36（3）：263-268

阚荣举，张四昌，晏凤桐，等1977我国西南地区现代构造应力场与现代构造活动特征的探讨[J].地球物理学报，20（2）：96-109

李长军，甘卫军，秦姗兰，等2018滇中主要活断层现今活动性研究[J]地震研究，41（3）：381-389

李玶，汪良谋1975云南川西地区地震地质基本特征的探讨[J].地质科学，10（4）：308-326

李铁明，邓志辉，吕弋培2003川滇地区现今地壳形变及其与强震时空分布的相关性研究[J].中国地震，19（2）：132-147

李延兴，杨国华，李智，等2003中国大陆活动地块的运动与应变状态[J].中国科学：地球科学，33（增刊1）：65-81

廖思佩，侯强，杜永超2016基于GPS形变资料的川滇地区应力场数值模拟研究[J].大地测量与地球动力学，36（7）：645-649

毛玉平，韩新民2003云南地区强震（M≥6）研究[M].昆明：云南科技出版社

申重阳，王琪，吴云，等2002川滇菱形块体主要边界运动模型的GPS数据反演分析[J].地球物理学报，45（3）：352-361

施发奇，尤伟，付云文2012GPS资料揭示的小江断裂近期运动特征[J].地震研究，35（2）：207-212

宋方敏，李如成，徐锡伟2002四川大凉山断裂带古地震研究初步结果[J].地震地质，24（1）：27-34

苏有锦，秦嘉政2001川滇地区强地震活动与区域新构造运动的关系[J].中国地震，17（1）：24-34

唐文清，陈智梁，刘宇平，等2005青藏高原东缘鲜水河断裂与龙门山断裂交会区现今的构造活动[J].地质通报，24（12）：1169-1172

唐文清，刘宇平，陈智梁，等2007基于GPS技术的活动断裂监测——以鲜水河、龙门山断裂为例[J].山地学报，25（1）：103-107

汪斌，李新平，周桂龙2012我国大陆实测深部地应力分布规律研究[J].岩石力学与工程学报，31（1）：2875-2880

王辉，曹建玲，张怀，等2007川滇地区下地壳流动对上地壳运动变形影响的数值模拟[J].地震学报，29（6）：581-591

王凯英，马瑾2004川滇地区断层相互作用的地震活动证据及有限元模拟[J].地震地质，26（2）：259-272

王铠元，孙克祥，薛啸锋1987川滇交界金河—洱海深断裂带的形成机理和构造属性探讨[J].地震研究，10（3）：92-105

王敏，沈正康，牛之俊，等2004现今中国大陆地壳运动与活动块体模型[J].国际地震动态，（增刊1）：21-32

王阎昭，王恩宁，沈正康2008基于GPS资料约束反演川滇地区主要断裂现今活动速率[J].中国科学：地球科学，38（5）：582-597

熊探宇，姚鑫，张永双2010鲜水河断裂带全新世活动性研究进展综述[J].地质力学学报，16（2）：176-188

熊熊，許厚泽，滕吉文2001青藏高原物质东流的岩石层力学背景探讨[J].大地测量与地球动力学，21（2）：1-6

徐锡伟，闻学泽，陈桂华，等2008巴颜喀拉地块东部龙日坝断裂带的发现及其大地构造意义[J].中国科学：地球科学，38（5）：529-542

曾融生，孙为国1992青藏高原及其邻区的地震活动性和震源机制以及高原物质东流的讨论[J].地震学报，（增刊1）：534-564

张竹琪，张培震，王庆良2010龙门山高倾角逆断层结构与孕震机制[J].地球物理学报，53（9）：2068-2082

赵德安，陈志敏，蔡小林，等2007中国地应力场分布规律统计分析[J].岩石力学与工程学报，26（6）：1265-1271

赵静，江在森，牛安福，等2015川滇菱形块体东边界断层闭锁程度与滑动亏损动态特征研究[J].地球物理学报，58（3）：872-885

朱守彪，石耀霖2004用遗传有限单元法反演川滇下地壳流动对上地壳的拖曳作用[J].地球物理学报，47（2）：232-239