基于有限元方法的冕宁跨断层水准大幅变化机理

岳冲 赵静 牛安福 吉平 李晓帆 苏琴 杨凡

摘要：分析计算安宁河断裂带上冕宁跨断层水准场地原始观测数据及年均速率结果，结合1-1C测段加密观测数据的差异性特征，重点对滑坡、抽水2个主要影响因素，利用有限元方法构建该场地三维Mohr-Coulomb弹塑性滑坡体模型及二维抽水影响的单向流固耦合模型，模拟获得2种主要因素的影响模式。结果显示：冕宁水准1-1C测线中距离1C点约300 m范围为古滑坡体所处边坡的坡脚区域，古滑坡体对坡脚区域长期的累积作用使得该区域存在隆升变形特征，进而使1-1C测段在抽水造成的大幅下降过程中出现沉降量逐渐减小后再次增大的特征;地下水抽取造成沿水井形成了较为明显的抽水漏斗沉降变形特征，在含水系统平衡被打破后，1号点所处土体边坡更易受到抽水带来的沉降变形影响，出现与井口处相当甚至比井口附近更大的沉降变形。

关键词：冕宁跨断层水准;有限元法;Mohr-Coulomb;流固耦合模型

中图分类号：P315.725文献标识码：A文章编号：1000-0666（2020）04-0666-08

0引言

断层作为地壳形变的薄弱地带，强地震的孕育及发生往往伴随着显著的断层变化（王永安等，2003，2004;岳冲等，2019）。定点跨断层测量作为断层形变监测的有效手段，能够直接反映测点处断层段落的运动变化，但受布设场地等环境因素的影响，并不能保证所有的观测桩均布设在基岩上，特别是位于川滇地区的跨断层场地，由于地形起伏变化较大，部分观测桩布设于山体底部土体内，观测桩稳定性的差异会影响测项异常变化的可信程度。在实际的异常分析过程中，除了分析跨断层场地所跨断层本身调整运动带来的变化以外，更多需要结合跨断层水准场地条件分析观测桩周边地下水位升降、人为载荷干扰、坡体稳定性变化等因素的影响（马伶俐等，2019），只有通对比分析不同影响模式下变形特征，才能结合实际异常测项变化，进一步对异常信度等做出评判。

位于安宁河断裂带上的冕宁跨断层水准场地3-1测段水准自2016年出现大幅下降变化，国内学者进行了大量研究工作，赵强和苏琴（2018）利用BABS-InSAR技术获取冕宁水准场地地面沉降结果，提出该场地受温泉开采影响的结论，但受测点所处山坡植被的影响，未能给出测线变化较大的1号点及所处边坡的位移变化，其计算结果显示沉降漏斗中心最大的沉降速率为27 mm/a，与1号点附近过渡点沉降量值相当，与沉降漏斗沉降量向边缘递减的认识并不一致;马伶俐等（2019）针对该测项周边滑坡体及抽水井等影响因素的现场落实以及数据的进一步跟踪分析提出了场地受温泉单一因素影响的质疑。因此针对冕宁跨断层水准场地大幅变化的成因及机理，本文重点讨论滑坡及抽水2个影响因素，通过有限元方法构建冕宁场地三维Mohr-Coulomb弹塑性滑坡体模型及二维抽水影响的单向流固耦合模型，获得滑坡及抽水2种因素的影响模式，进而对跨断层水准异常变化的成因机理进行研究。

1研究区概况及异常分析

1.1研究区概况

本文主要研究位于安宁河断裂带上的冕宁跨断层水准场地（28.56°N，102.21°E），据卫星影像和实地观察，安宁河断裂带在该场地主要由安宁河东支断裂（主断裂）和安宁河西支断裂组成（图1），安宁河西支断裂在冶勒盆地大坝子切割全新世洪积扇，形成边坡脊、断塞塘等显著的断错地貌，呈近NS向延伸;安宁河东支断裂从紫马垮向南经彝海、小盐井至马尿河边，总体走向近NS，沿断裂走向分布有串珠状断塞塘、断层崖、反向陡坎等断错地貌，表明其具有强烈的全新世活动性。东支断裂具有规模较大的破碎带，冕宁水准场地3号位于基岩破碎带西侧，其北侧存在2～3条分支断裂，向南在马尿河南岸至秧柴沟一线，形成显著的槽谷地貌。地震研究43卷第4期岳冲等：基于有限元方法的冕宁跨断层水准大幅变化机理

1.2冕宁跨断层水准异常分析

冕宁跨断层水准测线由3个观测端点构成（图1c）（赵强，苏琴，2018;马伶俐等，2019），其中1号点位于断层下盘，其台基位于山坡底部，主要由松散的冲洪积块碎石层组成，充填物以砂为主，含少量泥质成分;2号点位于河流西侧，台基为昔格达组（Q1x）深灰色半固结的泥岩、粉砂岩;3号点位于断层上盘，其台基处于马尿河右岸阶地上，阶地堆积物由大块的漂石、卵石及砂组成，厚度约2 m，阶地基座为花岗岩。整个场地长约3 km（其中1-2测段测线长约1.67 km，2-3测段测线长约1.33 km），为单线水准场地。整个观测测线建于1980年11月，于1981年9月开始跨断层短水准观测，观测之初，每年观测3～4个周期;1996—2012年，每年观测6个周期;2013年4月以后，每年观测12期，于2017年12月停测。

由图2a原始曲线可以看出：自20世纪80年代开始观测以来，3-2测段数据变化较为稳定，证明2号、3号点稳定性较强，观测期间未发生显著相对高差变化;相较之下，3-1测段于2006年之后整体运动速率有所增强，证明1号点自2006年以来出现下沉变化，自2016年6月开始，出现显著压缩变化，达19 mm，而3-2测段未出现同步异常变化，证明该测段变化可能主要由1号点引起。

为进一步定量分析跨断层水准变化，排除降雨、气温等年度影响因素，本文通过计算断层年均速率对比2个测段的活动强度变化（江在森等，2001;岳冲等，2017）。其中，断层的年均速率计算公式为：vt=1n∑ni=1（hti-ht-1i）（i=1，2，…，n）（1）式中：n为一年中观测期数;hti为第t年的第i个观测值。

图2b计算结果显示：2016年前，3-1，3-2测段年均速率分别为-0.27，0.003 mm/a，3-1测段年均速率高于3-2测段，也证明1号测点的活动性较强，但2个测段年均速率均在0值附近變化，也反映冕宁跨断层水准场地所在区域较弱的断层垂向活动特征。

为进一步获得1号点周围区域内位移下降变化，收集2016年1月以来1-1C测段（图1c ）内观测过渡点结果，测量结果（以2016年1月测量结果为对比数据，其它月份均与此月结果对比）如图3所示。通过差分数据可以看出：2016年1—6月（图3a），距离1号点800 m范围内，点位变化相对稳定，而在距离1号点约800～1 300 m范围内出现小幅隆升-下降变化;自2016年7月起（图3b），距离1号点800 m范围内陆续出现大幅沉降变化，沉降变化时段主要集中于2016年7—10月，2016年11月后整个测段沉降速率降低，其中最大的沉降点位于距离1号点64.5 m区域，最大沉降幅度达到27.5 mm;而在距离1号点约800～1 300 m范围出现先隆升后下降变化，也证明在该区域内可能存在其它影响因素，使得其变化与0～800 m范围内存在较大差异。

2计算方法及模型设计

2.1强度折减法

针对边坡稳定性的分析方法主要分为极限平衡法和有限元（或有限差分）法，而有限元法中强度折减法通过分析獲得安全系数，既能保持有限元模拟复杂问题的优势，又能给出明确的安全系数，在边坡稳定性模拟计算中得到了广泛的应用（郑颖人等，2005;黄盛铨等，2008;吕庆等，2008;江权等，2009;张强等，2011;Griffiths，Lane，1999）。

强度折减法最早由Zienkiewic等（1975）提出，抗剪强度折减系数定义为：在外载荷保持不变的情况下，边坡内土体所能提供的最大抗剪强度与外载荷在边坡内所产生的实际剪应力之比。折减后的抗剪强度参数为：c′=c/Fcφ′=arctan（tanφ/Fφ）（2）式中：c，φ分别为材料的峰值粘聚力和内摩擦角;c′，φ′分别为临界破坏状态下材料的粘聚力和内摩擦角;Fc和Fφ为粘聚力和抗滑动折减系数，通过不断迭代模拟计算，逐渐增大折减系数，直至边坡达到发生失稳破裂的临界状态（计算不收敛的边界）。

2.2滑坡体模拟

本文以1号点所处边坡为主要研究区域设计三维滑坡体模型，模型尺寸为800 m×1 000 m，整个边坡高度参考场地等高线分布为300 m（马伶俐等，2019），边坡角度为30°，模型示意图见图4a。本文使用有限元软件ABAQUS完成整个模型构建及模拟计算工作，其中采用线性六面体单位（C3D8R）进行网格剖分，网格间距为10 m，整个模型共获得159 200个单元，171 296个节点。

1号点台基主要由松散的冲积碎石组成，内含砂、泥质等（马伶俐等，2019），土体采用理想线弹塑性Mohr-Coulomb模型，模型的基本参数设置如下：密度为1 500 kg/m3，弹性模量为10 MPa，泊松比为0.35，粘聚力为20 kPa（陈正汉，郭楠，2019），内摩擦角为30°（李宁等，2018;Cai，Ugal，2004）。因沉积土深度的不确定性，分别考虑沉积层厚度分别为20，40，60，80 m和全尺寸沉积土共5个模型，以探索不同厚度沉积层对地表位移的影响，其中基岩考虑为弹性体，密度为2 630 kg/m3，弹性模量为300 MPa，泊松比0.36（陈卫兵等，2007;吕庆等，2008;许瑞宁，2016;陈正汉，郭楠，2019）。

在整个模拟计算过程中，模型的边界条件设置如下：模型底部固定，其西侧和东侧边界沿X方向固定，其北侧和南侧边界沿Z方向固定，模型表面保持自由（边界条件按照不同坐标轴颜色进行区分）。计算时考虑重力作用影响，重力加速度为9.8 m/s2。计算步骤如下：首先对模型施加重力载荷作用，建立模型应力平衡状态;在应力平衡状态下逐渐增加强度折减系数（折减系数的峰值为2），迭代计算，直至数值不收敛，以此作为土坡达到临界破坏的标准。

2.3抽水影响模拟

本文结合冕宁跨断层水准场地周边遥感影像（图1a）、等高线分布及异常落实结果（马伶俐等，2019），构建跨1号点至水井的二维剖面模型（图4b）。根据InSAR等监测结果显示，沿短轴方向沉降范围约0.4 km（赵强，苏琴，2018），因此，整个模型设计长600 m（水井中心距左侧边界400 m），井口处高200 m，1号点距井口215 m，1号点所处边坡比井口处高50 m，水井深度为70 m。根据场地周边大于10 m的沉积层结果（马伶俐等，2019），井口以下20 m设置为沉积层（弹性模量5 MPa，泊松比0.4）（金小荣等，2005;曾超峰，2014）;马伶俐等（2019）认为该水井可能将55 m处含水细层打通，因此，井下20～55 m处设置为砂岩（弹性模型10 MPa，泊松比0.32），井口55 m以下为昔格达组地层（杨氏模量80 MPa，泊松比0.4）（孙长升，2015;许瑞宁，2016;张伟松，卢海峰，2017;杨宝，2018;杜宇翔等，2020）。整个模型两侧沿水平方向固定，底部沿垂直方向固定，其上侧地表均匀施加抽水引起的地表压强变化，水井井底和边坡表面设置渗透系数为0，主要模拟沉降变化最大时段2016年7—10月共计4个月抽水量引起的地表变形。

3模拟结果及分析

3.1累积位移变化分析

本文重点研究边坡失稳对跨断层水准的影响，共模拟5组不同厚度沉积层模型失稳破坏前累积位移结果。通过在模型地表构建剖面线（图5a中黑线所示，从模型左侧至右侧点号依次为1～86号，其中，坡脚点为34号点，坡顶点为75号点），模拟结果显示折减系数最大达1.16时，模型达到失稳破坏的临界状态，我们获得接近失稳前折减系数为1.1时，5组模型地表垂向累积位移结果如图5b所示。从模拟结果可以看出：松散沉积层厚度由20 m增加至全尺寸的沉积土模型，在滑坡体边坡处的变形区（位移绝对值大于0 m）逐渐向坡脚处（34号点）扩展，在坡顶处变形区逐步向模型右侧（86号点）扩展，即滑坡失稳造成的变形范围以及沉降位移量均显著增大，最大位移量变化主要集中在坡顶点附近区域;而随着沉积层厚度的不断增加，失稳带来的坡体运动逐步在坡脚区域累积，进而造成边坡坡脚区域地表的隆升变形。

为了进一步展示边坡失稳初期模型内部位移特征，图5a为强度折减系数为1.01时的全尺寸沉积土模型位移结果。图5a显示，随着土边坡整体强度不断降低，坡顶区域土体逐渐失稳并向边坡中部区域迁移运动，当强度折减系数达到1.16时，在边坡中部区域形成较大垮塌位移，进而造成滑坡运动，如图1a中紫色框即为在山坡中上部区域发生的古滑坡迹象，而滑坡发生后整个坡顶区域土体向坡体中部及坡脚区域堆积，进一步带来坡脚区域的隆升。而本文所研究的跨断层水准1-1C测线距离1C点约300 m范围即为该古滑坡体所处边坡的坡脚区域，古滑坡体对坡脚区域的长期累积作用使得该区域存在隆升变形特征，进而使1-1C测段在大幅沉降过程中出现沉降量先减小后增大的迹象。而冕宁水准2号点跨过了该滑坡区域，且2号点布设于半固结的昔格达组地层上，相对较为稳定，因此滑坡体并未影响2号点变化，冕宁3-2测线变化相对较为稳定。3.2抽水影响模拟结果

施加4个月抽水量模拟地表沉降的结果如图6所示。为对比抽水量变化对地表沉降量的影响，模型模拟抽水量由10 m2/h（图6b中模型1）增加至20 m2/h（图6b中模型2），图6a展示抽水量为20 m2/h时的地表沉降位移结果。从图6a可以看出：伴随着不断抽水，原来平衡状态下的含水系统被打破，孔隙水压力减小，导致有效应力变化，进而造成土体变形，沿水井形成了较为明显的抽水漏斗沉降变形，而整个左侧边坡在抽水过程中作为自然排水面，加之较为松散的沉积层结构，存在放大沉降位移量现象;从具体沉降量（图6b）可以看出，模型1中井口处沉降量达到13.3 mm，井口左侧边坡区域作为自由排水面，在整个抽水过程中坡脚区域随水体逐步渗透排出、沉积土体固结并向坡脚区域累积，沉降量出现先减小后逐渐增加趋势，最大达到15.8 mm后逐步稳定;而当抽水量达到20 m2/h时，整个区域沉降量明显增加，同时整个边坡沉降变化速率也显著增大，其中井口处沉降量达到26.6 mm，坡脚最小沉降量为24.7 mm，1号点区域沉降量达到31.6 mm。

因此，作为土体边坡上的1号点，因受边坡结构本身以及沉积层的影响，在整个抽水过程中，更易受到抽水带来的沉降变形，出现与井口处相当甚至比井口附近更大的沉降变形，但后期随着抽水量的稳定以及地表水、地下水的不断补给，地下含水系统将再次达到平衡状态，地下水水位下降一定幅度后将趋于稳定，进而地表沉降位移也将趋于稳定。

4结论

本文针对冕宁跨断层水准3-1测线出现的大幅沉降变化异常，计算并分析了冕宁场地年均变化速率、测段加密数据，结合异常落实过程中的滑坡及抽水2个主要因素，利用有限元方法构建冕宁场地Mohr-Coulomb弹塑性滑坡体模型及抽水影响的单向流固耦合模型，通过模拟获得2种主要因素的影响模式，结果显示：

（1）断层年均速率结果显示3-1，3-2测段年均速率分别为-0.27和0.003 mm/a，3-1测段运动速率高于3-2测段，也证明1号点活动性较强;1-1C测段加密观测结果显示沉降变形主要集中于2016年7—10月，最大的沉降点位于距离冕宁1号点64.5 m点位区域，最大沉降幅度达到27.5 mm，1-1C測段沉降量呈“增加—减小—增加—减小”的变化特征。

（2）模拟5组不同厚度沉积层模型失稳破坏前累积位移结果显示：随松散沉积层厚度的不断增加，失稳变形的范围逐步扩大，边坡沉降位移也逐渐增大，最大位移变化主要集中在坡顶点附近区域，失稳带来的坡体运动逐步在坡脚区域形成累积，进而造成边坡坡脚区域地表的隆升变形。冕宁水准1-1C测线中距离1C点约300 m范围为古滑坡体所处边坡的坡脚区域，古滑坡体对坡脚区域长期的累积作用使得该区域存在隆升变形特征，进而使1-1C测段在抽水造成的大幅下降过程中出现沉降量逐渐减小后再次增大的迹象。

（3）模拟4个月抽水造成地表沉降的结果显示：沿水井形成了较为明显的抽水漏斗沉降变形特征，1号点所处土体边坡，因本身自由排水面以及松散沉积层的影响，在含水系统平衡被打破后，更易受到抽水带来的沉降变形影响，出现与井口处相当甚至比井口附近更大的沉降变形。

受该地区地下水文资料、水井钻孔资料等限制，本文对地下介质分布、物性参数等进行必要的简化，模拟得到的造成1号点沉降变化的抽水量介于10～20 m2/h，与实地异常落实的10 m2/h结果（马伶俐等，2019）基本一致，通过定量模拟滑坡与抽水两种影响因素的影响模式，为跨断层测量出现大幅变化后的判定提供定性分析的依据，未来随着该区域研究资料的进一步完善，作者也将进一步跟踪补充该区域的研究结果。

两位审稿人对本文进行认真的审阅并提出了宝贵建议，四川省地震局测绘工程院提供高质量的测量数据，在此表示衷心感谢。

参考文献：

陈卫兵，郑颖人，雷文杰.2007.沉埋桩加固滑坡体模型试验的三维有限元模拟[J].岩土力学，28（增刊1）：273-276.

陈正汉，郭楠.2019.非饱和土与特殊土力学及工程应用研究的新进展[J].岩土力学，40（1）：1-54.

杜宇翔，盛谦，王帅，等.2020.昔格达组半成岩微观结构与力学性质研究[J].岩土力学，41（4）：1-13.

黄盛铨，刘君，孔宪京.2008.强度折减DDA法及其在边坡稳定分析中的应用[J].岩石力学与工程学报，（增刊1）：2799-2806.

江权，冯夏庭，向天兵.2009.基于强度折减原理的地下洞室群整体安全系数计算方法探讨[J].岩土力学，30（8）：2483-2488.

江在森，丁平，王双绪，等.2001.中国西部大地形变监测与地震预报[M].北京：地震出版社.

金小荣，俞建霖，祝哨晨，等.2005.基坑降水引起周围土体沉降性状分析[J].岩土力学，26（10）：54-60.