基于正交各向异性理论表征断层的变形行为
石富强 朱琳 李玉江 丁晓光
摘要：总结了“介质弱化带”和“摩擦接触”这2类常见断层模型在地壳动力学数值模拟研究中的应用研究进展，再利用2者各自的特点，引入复合材料力学中各向异性理论，结合断层产状给出了针对不同走向、倾角的断层各向异性本构方程。数值模拟试验对比研究显示：基于正交各向异性本构关系的断层模型不但能够避免断层垂直方向上的形变场畸变现象，而且能够有效地刻画断层两盘相对滑动，模拟结果与基于接触摩擦断层模型的模拟结果吻合良好。该模型可以在研究区域构造复杂，接触计算收敛困难的情况下替代摩擦接触模型。
关键词：断层变形；数值模拟；介质弱化模型；接触摩擦模型；正交各向异性理论
中图分类号：P313.5 文献标识码：A 文章编号：1000-0666（2018）01-0064-09
0 引言
地壳形变是地壳在力学、热学等复杂环境作用下表现出来的综合响应，研究地壳形变特征及其演化规律是认识地壳运动及其伴生地质灾害的一个重要途径。随着计算机科学的日益发展，数值模拟技术已经成为地壳动力学机理研究的重要手段（Moreno et al，2010；Hergert，Heidbach，2010；Li et al，2014；Zhu，Zhang，2010；He et al，2013；Lu et al，2011；陈连旺等，1999，2008；邵志刚等，2008；李玉江等，2014）。在复杂地质构造作用下，断层作为岩石圈岩石破裂形成的特殊产物，在现今地壳形变特征以及地震的孕育、发生中有重要的作用。在探索地震发生孕育的过程中，国内外科学家从不同角度阐述了断层物理机制（Sib-son，1977；马胜利，本利彦，1995）。但由于地质构造运动的复杂性，现在还没有一种行之有效的办法可以模拟真实地壳环境中断层的形成、发展变化。因此通过数值模拟解释现今地壳形变特征时，往往需要结合野外考察、深地震反射、地震精定位、震源机制解等多学科方法给出的断层产状及断裂带断层泥等做出一定程度的简化，进而构建断层物理模型。
在长期的实践应用过程中，根据物理模型及具体问题的研究，这些简化的断层模型可以归纳为：劈节点模型、块体模型、介质弱化模型以及摩擦接触模型（和平等，2011）。劈节点模型本质上是在预设断层位错的情况下模拟断层，因而在同震或震后效应研究中具有极大的优势，得到了廣泛应用；块体模型本质基础是Shi（2007）提出的DDA（Discontinuous Deformation Analysis）算法，它将地球动力过程看作为一系列块体相互作用的过程，断层其实就是块体与块体之间的边界。介质弱化模型和摩擦接触模型作为2种原理简单且便于在一般软件中设置的模型，在区域地壳动力学形变特征模拟及断层活动状态模拟中得到大量应用。由于断层的介质力学参数对区域应力场以及形变场均有一定程度的影响（李玉江等，2010），因此在模拟过程中，介质弱化模型是将断层视为有一定厚度的弱介质层，通常通过降低该介质层弹性模量来模拟断层对区域地壳动力过程的影响（Lu et al， 2011，2012；祝爱玉等，2015a；胡勐乾等，2014；杨兴悦等，2013）；摩擦接触模型则是直接在创建好的模型基础上给断层面赋摩擦系数，进而实现断层两盘的相对运动（Hergert，Heidbach，2010； Zhu，Zhang，2010； He et al，2013；李玉江等，2014；祝爱玉等，2015b）。但也存在一些问题：①对于介质弱化模型来说，其本质上是通过降低断层介质弹性模量来近似模拟断层两盘的非连续变形。尽管野外观测到断层带内有糜棱岩或断层泥等弱力学属性材料填充，但这种弱力学属性的填充物已经在长期的构造作用下与周边介质处于力平衡状态。因此，在构建有限元模型时直接赋予该填充物一个较低的弹性模量，以期获得断层两盘的运动差异和断层附近应力状态的同时，往往忽略了降低弹性模量，也使得该填充层内部、垂直于断层走向的力平衡破坏，在断层内部垂直于断层方向产生新的附加力。这可能影响模拟结果的可靠性。此外，实验室关于断层泥的弱力学属性对断层行为的影响也表现在剪切摩擦行为上（马胜利，1986）。在各向同性弹性理论前提下还可能会出现不合理的负泊松比的情况（董培育，石耀霖，2013）。②摩擦接触模型可以很好地刻画断层两盘相对运动的非连续变形，相比介质弱化模型更接近于真实断层状态。但摩擦接触计算是个强非线性问题，涉及到收敛与否以及计算量等问题（李玉江等，2009），特别是在研究区域内断层较多的情况下，如果将每条断层都构建为摩擦基础模型，收敛问题将变得更加紧迫，为此研究中还可以将重点断层考虑为摩擦接触，而将次要断层处理为弱化带模型（Li et al，2015），这一简化处理提高了计算效率，但又增加了弱化介质模型自身的间题。
因此，如何做到既保证计算效率、避免介质弱化模型可能带来的不真实信息，又能快速合理模拟断层两盘相对运动是地壳动力学数值模拟的一个基础性工作。本文在介质弱化带模型基础上引入正交各向异性理论，用正交各向异性本构方程替代传统弱化带中的线弹性本构方程，实现了弹性模量与剪切模量相互独立；同时考虑到实际建模过程中可能遇到不同断层走向、倾角各不相同的情况，通过坐标转换给出了针对任意走向、倾角的断层本构关系。最后通过数值试验对比分析了介质弱化模型、摩擦接触模型以及正交各向异性断层模型。
1 正交各向异性本构关系
假定断层带为夹在块体之间的类似于层合板的一类特殊介质：其在垂直于断层方向力学性质与周边块体相近，不易产生大挤压或拉张变形；而在平行于断层运动方向上的剪切模量较小，能够产生较大的剪切变形。则根据复合材料力学（沈观林，胡更开，2006）基础知识，断层内应力应变关系在其3个主方向上满足如下关系：式中：ε=[e11 e22 e33 e12 e23 e13]T，σ=[s11；s22 s33s12 S23 S13]T分别为3个主方向的应变和应力分量；x为柔度矩阵；E，，E2，E3，v12，V23，v13，G12，G23，G13分别为其3个主方向的工程弹性常数。K可表示为：
由于本文重点考虑断层介质在走向和倾向的剪切行为，其他方向上的介质力学性质与周边介质一致。为此将式（2）进一步简化为：式中：E，v分别为周边介质弹性模量和泊松比；Gd；P为断层在剪切方向上的剪切模量，独立于周边介质参数E，v。
在实际操作过程中，由于区域断层较多且其走向和倾角各不相同，因此需要将其主方向下的应力应变矩阵旋转变换到地理坐标（x，y，z）下。图1为断层面主方向（1，2，3）以及地理坐标系（x，y，z）相对关系，其中：1，2，3分别为断层走向（面向走向，右手在上盘）、倾向（逆冲为正）以及张拉方向（张为正）；x，y，z分别为地理东、北向以及垂直地表向上的方向。α为断层走向，β为断层倾角。令断层面应变和应力分量在地理坐标系（x，y，z）下分别为：
e=[exx eyy ezz exy eyz exz]T
s=[sxx syy szz sxy syz sxz]（4）
由图1b可知，将地理坐标系下的应力张量绕z轴逆时针旋转，然后再绕主轴1逆时针旋转β便可得到主方向下的应力张量。则：式中：T=A*B*结合式（1）和（4），式（5）可以转化为：式中：F为坐标转换矩阵，可表示为：同理可得：将式（6）和式（7）带入式（1），可得：
F·e=K·F·s（8）即：
e=D·s=F-1·K·F·s（9）式中：D为正交各向异性材料在地理坐标系（x、y、z）下的柔度矩阵，且：式中：柔度矩阵D为21个独立参数的对称矩阵，因此在实际模拟计算过程中，只需修改正交各向异性本够关系柔度矩阵K在其相应主方向下对应滑动方向的剪切模量（式（3）中的Gslip），结合坐标转换矩阵F，代入式（10）便可生成断层有限元模型的柔度矩阵，进而可以构建有限元计算的整体刚度矩阵，实现垂直断层方向无大变形，平行断层滑动方向有大变形滑动的断层运动状态。为方便控制，本文令断层走向和倾向剪切模量与周边弹性介质剪切模量之比为λ∈[0，1]，即Gslip=λ
2 数值试验对比研究
2.1 模型构建
为检验正交各向异性本构关系模拟断层的可行性，本文设计1个虚拟理想弹性地壳模型，模型为边长L=6000m的正方形，纵向（z）深2000m且上下地殼各1000m厚，模型中部沿x轴走向有一倾角90°的右旋走滑断裂带，长2000m，宽50m。模型施加狄里克莱边界条件：-y边界约束y方向位移为0，+y边界只有y方向位移，Uy=-0.1m；-x边界上施加一反正切函数的x方向剪切位移，如图2a所示；模型-z底面约束z方向位移为0。假定模型弹性模量E=84GPa，泊松比v=0.25。然后针对常用的介质弱化模型，摩擦接触模型以及本文正交各向异性模型分别设计2组工况：第一组工况（Case 1）假定上地壳断层闭锁，下地壳可以剪切滑动；第二组工况（Case 2）认为整个断层处于非闭锁状态，断层上下均可以相对自由滑动，如图2b所示。
2.2 结果对比
在一个地震轮回过程中，断层在地震发生时由于快速滑动而破裂，使得其力学属性快速降低；在震后经历长时间的愈合，在下一次地震来临前，断层力学属性又增加到相对较强的状态。根据地震轮回的理论模式（Savage，Burford，1973；Scholz，1998），借鉴Hergert和Heidbach（2010）在马尔马拉海的工作思路，本文首先用介质弱化模型、摩擦接触模型以及本文给出的正交各向异性本构模型分别模拟一组震间地壳运动位移场图像。假定上地壳己经经历了长久的愈合，断层力学属性与周边地壳介质基本一致，而深部断层强度低，依然可以自由剪切，断层模型参数设置如图2b（Case 1），模拟结果如图3所示。
基于3种断层模型地壳运动位移场图像的空间分布形态基本一致，但一般文献多用的介质弱化模型在断层附近局部，x方向的位移场存在一定程度的畸变，但并未影响Ux沿y的反正切渐变分布特征。另外在空间分布图像上，基于介质弱化断层模型的Ux和基于摩擦接触模型及正交各向异性模型的Ux结果在定量上有所差异。而基于摩擦接触断层模型和各向异性断层模型的对比结果显示，2者在位移场图像上高度一致，这与本文构建模型的初衷一致，因此在考虑断层较多，摩擦接触建模复杂，计算收敛困难的时候可以考虑用各向异性本构关系来模拟断层，其模型构建与计算量与常用的介质弱化模型完全一致，计算中也不会涉及到非线性收敛困难等情况，不同的是有限元方程中的柔度矩阵（或刚度矩阵）内的具体数值有所不同。
震后短期内，断层没有完全愈合，其介质力学属性依然较低。在分析这些介质力学属性本来偏低的断层时，本文将整个断层的力学参数降低，具体如图2b（Case 2）所示。结果显示基于3种断层模型地壳运动位移场图像空间分布不同程度出现差异，常用的介质弱化模型在水平和垂直位移场上明显与其他2种模拟不一致，在断层附近位移场图像出现巨大的畸变，这与断层弹性模量的降低引起垂直于断层方向发生大幅度的挤压变形有关。实际上由于长期的构造作用，断层两盘的相对运动一般只发生在其走向和倾向滑动方向，并不会在断层上存在这样垂直于断层的大幅挤压变形，因此在应用简单的介质弱化断层模型模拟地壳形变时，随着弹性模量的降低，断层附近的模拟位移场图像将出现不同程度畸变，进而可能会影响到断层附近的模拟应力状态（董培育，石耀霖，2013）等。对比基于摩擦接触模型和本文各向异性模型的模拟结果，尽管2者在断层附近有一定程度的差异，但整体上2者空间分布形态完全相似，可以清晰地给出断层两盘的相对滑动，所不同的是由于摩擦接触模型给断层赋予的是摩擦系数，而本文只是近似给的剪切模量，2者虽然在阻碍变形的机制上相似，但其背后的物理方程不同，因此造成2者在断层两盘相对滑动量的差异。后续的研究需要基于地壳运动位移场图像讨论2者之间的统计关系，这样会进一步提高本文各向异性断层模型在模拟中的应用和理解。
为进一步分析3类断层模型的异同，从图3和图4中切出x=0剖面上的Ux和Uy做对比分析，如图5所示。结果显示当在Case 1（浅部断层愈合，深部断层可以滑动）情况下，总体上这3类断层模型给出的结果基本一致，但从更细致的角度看，基于摩擦接触断层模型和正交各向异性断层模型给出的结果吻合一致。但介质弱化模型给出的结果与其他2种断层模型给出的模拟结果略有差异，并且这种差异在空间上呈现近断层处相对较小，远离断层处增大。这是由于浅部断层与周边介质属性一致，保证了表层形变场的连续性，但基于介质弱化的深部断层模型破坏了垂直于断层方向的平衡状态，產生附加的形变场畸变信息。进一步降低表层断层的介质力学属性（即Case 2），这种畸变信息对表层形变场的影响更加显著，特别是在垂直于断层方向上（图Sd）。但在平行于断层方向上，3类断层模型给出的模拟结果均显示出不同程度的差异，这是因为本文针对3类断层模型的力学参数都是随机赋值，没有探讨3者力学参数之间的等效关系。但3者均可以表征本文预设的断层走滑行为。对于垂直于断层方向，图Sd显示基于正交各向异性理论的断层模型与摩擦接触断层模型的模拟结果一致，表明正交各向异性断层模型相比介质弱化模型，更接近于实际情况。
为进一步分析不同断层模型对断层附近模拟应力状态的影响，考虑到图2设计的模型表现为剪切变形行为，以水平剪应力（τyx）的空间分布为讨论对象。在图2的基础上继续分析了不同1况下，不同断层模型对应的τyx模拟结果，如图6所示。结果显示：较介质弱化模型给出的τyx空间分布（图6a，d），正交各向异性断层模型的模拟结果（图6b，e）与摩擦接触模型的结果（图6c，f）更为相近。在上地壳断层闭锁的情况下（Case 1），水平剪应力τyx云图在空间上相对y=0轴对称（图6b、c），而介质弱化模型给出的结果不同于其他两个模型的模拟结果，水平剪应力呈现出不对称形态（图6a）；在整个断层均贯通非闭锁的情况下（Case 2），尽管正交各向异性断层模型的模拟结果（图6e）与摩擦接触模型（图6f）的模拟结果有一定的差异，但2者基本保持相似的空间分布特征，而基于介质弱化模型给出的模拟结果（图6d）则在空间分布上表现出不规则的分布形态。这些都是由于在数值模拟中，介质弱化模型由于其本构关系的原因，不可避免地在垂直于断层方向上引入了较大的挤压变形（图4d），这种不合理的挤压变形不但影响着位移场的空间分布特征，同时也影响着应力场的空间分布。
3 讨论
近年来，随着GPS地壳速度场观测数据的丰富，基于GPS约束的地壳动力学数值模拟研究日渐深入（Hergert，Heidbach，2010；He et al，2013；Lei et al，2010； Cho，Kuwahara，2013；龙小刚，朱守彪，2015；刘峡等，2010）。一般情况下，研究人员通常都会通过一些途径选定有限元模型的边界条件，模拟计算给出模型内部的速度场模拟结果，然后通过对比模拟结果与观测结果之间的残差来确定模型是否可靠，参数选择是否合理，在大量比较之后给出合理模型的基础上分析区域的动力学机制及应力状态等，这就要求我们首先要保证所采取的模型是合理的。介质弱化模型作为最为简单的断层模型，其建模方便、不涉及非线性问题且不会存在难于收敛的情况，在一些研究中得到大量的应用。因此在应用介质弱化模型的情况下，在对比模拟地壳运动速度场与GPS观测速度场残差的同时，如何评价类似于图4a、d以及图5d的形变场畸变对模拟结果的影响也是不容忽视的问题。在不可避免使用介质弱化断层模型的时候，如能进一步修改其有限元刚度矩阵，使其剪切模量与弹性模量相互独立，对于提高模拟结果的可靠性具有一定的意义。
4 结论
本文在介质弱化带模型基础上，引入正交各向异性理论，同时考虑到实际建模过程中可能遇到不同断层的走向、倾角各不相同的情况，本文通过坐标转换给出了针对任意走向、倾角的正交各向异性断层本构关系。最后通过数值试验对比分析了介质弱化模型、摩擦接触模型以及正交各向异性断层模型在震间地壳速度长图像模拟中的差异，结果显示：（1）无论断层力学性质强弱，传统的介质弱化模型都会对地表地壳运动位移场产生一定程度的影响，断层力学性质（对应于图2不同断层模型中断层的力学参数，如弱化介质模型的弹性模量、各向异性模型的剪切模量以及摩擦接触模型的摩擦系数）越弱，影响越显著，这主要由于是降低断层弹性模量，使得断层内垂直于断层方向的力平衡被打破，进而影响了断层附近位移场图像的空间分布，这在模拟中同时也会对断层附近的应力场和应变场产生一定程度的影响：（2）从位移场图像以及应力场图像的空间分布来看，正交各向异性本构关系与摩擦接触模型给出的模拟结果基本一致，因此在区域构造复杂的情况下，可以考虑使用各向异性本构关系来表征断层的变形机制。
非常感谢陈连旺研究员的悉心指导，同时也感谢中国地震局访问学者计划的支持和地壳应力研究所提供良好的办公环境和软件支持！部分图件采用GMT（Wessel et al，2013）绘制，特此表示感谢！
参考文献：
陈连旺，陆远忠，张杰，等.1999.华北地区三维构造应力场[J].地震学报，21（2）：140-149.
陈连旺，张培震，陆远忠，等.2008.川滇地区强震序列库仑破裂应力加卸载效应的数值模拟[J].地球物理学报，51（5）：1411-1421.
董培育，石耀霖.2013.关于“用单元降刚法探索中国大陆强震远距离跳迁及主体活动区域转移”的讨论——横向各向同性“杀伤单元”才是更好的途径[J].地球物理学报，56（6）：2133-2139.和平，李志雄，陆远忠，等.2011.断层面的有限单元模拟方法综述[J].西北地震学报，33（2）：186-194.
胡勐乾，邓志辉，陆远忠，等.2014.三维数值模拟在华北地区现今构造变形分析中的应用[J].地震地质，36（1）：148-165.
李玉江，陳连旺，刘少峰，等.2014.芦山地震的发生对周围断层影响的数值模拟[J].地球学报，35（5）：627-634.
李玉江，陈连旺，叶际阳，等.2010.岩石物性变化对区域应力场的影响[J].地球物理学进展，25（6）：1941-1946：
李玉江，陈连旺，叶际阳.2009.数值模拟方法在应力场演化及地震科学中的研究进展[J]地球物理学进展，24（2）：418-431.
刘峡，马瑾，傅容珊，等.2010.华北地区现今地壳运动动力学初步研究[J].地球物理学报，53（6）：1418-1427.
龙小刚，朱守彪.2015.台湾碰撞带现今地壳变形场特征及其动力学成因的有限单元法模拟研究[J].地球物理学报，58（7）：2350-2365：
马胜利，本利彦.1995.蒙脱石的脱水作用对断层摩擦本构行为的影响[J].地震地质，17（4）：289-304.
马胜利.1986.模拟断层带摩擦滑动性状与变形特征[J].中国地震，（2）：74-80.
邵志刚，傅容珊，薛霆虓，等，2008.昆仑山MS8.1级地震震后变形场数值模拟与成因机理探讨[J].地球物理学报，51（3）：805-816.
沈观林，胡更开.2006.复合材料力学[M].北京：清华大学出版社.
杨兴悦，陈连旺，杨立明，等.2013.巴颜喀拉块体强震动力学过程数值模拟[J].地震学报，35（3）：304-314.
祝爱玉，张东宁，蒋长胜，等.2015a.川滇地区地壳应变能密度变化率与强震复发间隔的数值模拟[J].地震地质，37（3）：906-927.
祝爱玉，张东宁，蒋长胜.2015b.安宁河—则木河—小江断裂带应力状态分段特征的数值模拟研究[J].中国科学：地球科学，45（12）：1839-1852.
Cho I，Kuwahara Y.2013.Numerical simulation of crustal deformation u-sing a three-dimensional viscoelastic crustal structure model for theJapanese islands under east-west compression[J].Earth，Planetsand Space，65（9）：1041-1046.
He J，Lu S，Wang W.2013.Three-dimensional mechanical modeling ofthe GPS velocity field around the northeastern Tibetan plateau andsurrounding regions[J].Tectonophysics，584；257-266.
Hergert T，Heidbach 0.2010.Slip-rate variability and distributed de-formation in the Marmara Sea fault system[J].Nature Geoscience，3（2）：132-135.
Lei X，Chen Y，Zhao J，et al.2010.Modelling of current crustal tectonicdeformation in the Chinese Tianshan orogenic belt constrained byGPS observations[J].Journal of Geophysics and Engineering，7（4）：431.
Li S，Moreno M，Rosenau M，et al.2014.Splay fault triggering by greatsubduction earthquakes inferred from finite element models[J].Geo-physical Research Letters，41（2）：385-391.
Li Y，Chen L，Liu S，et al.2015.Coseismic Coulomb stress changescaused by the MW6.9 Yutian earthquake in 2014 and its correlationto the 2008 MW7.2 Yutian earthquake[J].Journal of Asian EarthSciences，105：468-475.
Lu Y，Yang S，Chen L，et al.2011.Mechanism of the spatial distributionand migration of the strong earthquakes in China inferred from nu-merical simulation[J].Journal of Asian Earth Sciences，40（4）：990-1001.
Lu Y，Yang S，Chen L，et al.2012.Migration trend of strong earthquakesin North China from numerical simulations[J].Journal of AsianEarth Sciences，50：116-127.
Moreno M，Rosenau M，0ncken O.2010.2010 Maule earthquake slip cor-relates with pre-seismic locking of Andean subduction zone[J].Nature，467（7312）：198-202.
Savage J C，Burford R O.1973.Geodetic determination of relative platemotion in central California[J].Journal of Geophysical Research，78（5）：832-845.
Scholz C H.1998.Earthquakes and friction laws[J].Nature，391（6662）：37-42.
Shi G.2007.Applications of discontinnons deformation analysis（DDA）to rock engineering[M].Berlin，Heidelberg，Computational Mechan-ics springer；136-147.
Sibson R H.1977.Fault rocks and fault mechanisms[J].Journal of theGeological Society，133（3）：191-213.
Wessel P，Smith W H，Scharroo R，et al.2013.Generic mapping tools：Im-proved version released[J].Eos，Transactions American GeophysicalUnion，94（45）：409-410.
Zhu S，Zhang P.2010.Numeric modeling of the strain accumulation andrelease of the 2008 Wenchuan，Sichuan，China，Earthquake[J].Bul-letin of the Seismological Society of America，100（5B）：2825-2839.
摘要：总结了“介质弱化带”和“摩擦接触”这2类常见断层模型在地壳动力学数值模拟研究中的应用研究进展，再利用2者各自的特点，引入复合材料力学中各向异性理论，结合断层产状给出了针对不同走向、倾角的断层各向异性本构方程。数值模拟试验对比研究显示：基于正交各向异性本构关系的断层模型不但能够避免断层垂直方向上的形变场畸变现象，而且能够有效地刻画断层两盘相对滑动，模拟结果与基于接触摩擦断层模型的模拟结果吻合良好。该模型可以在研究区域构造复杂，接触计算收敛困难的情况下替代摩擦接触模型。
关键词：断层变形；数值模拟；介质弱化模型；接触摩擦模型；正交各向异性理论
中图分类号：P313.5 文献标识码：A 文章编号：1000-0666（2018）01-0064-09
0 引言
地壳形变是地壳在力学、热学等复杂环境作用下表现出来的综合响应，研究地壳形变特征及其演化规律是认识地壳运动及其伴生地质灾害的一个重要途径。随着计算机科学的日益发展，数值模拟技术已经成为地壳动力学机理研究的重要手段（Moreno et al，2010；Hergert，Heidbach，2010；Li et al，2014；Zhu，Zhang，2010；He et al，2013；Lu et al，2011；陈连旺等，1999，2008；邵志刚等，2008；李玉江等，2014）。在复杂地质构造作用下，断层作为岩石圈岩石破裂形成的特殊产物，在现今地壳形变特征以及地震的孕育、发生中有重要的作用。在探索地震发生孕育的过程中，国内外科学家从不同角度阐述了断层物理机制（Sib-son，1977；马胜利，本利彦，1995）。但由于地质构造运动的复杂性，现在还没有一种行之有效的办法可以模拟真实地壳环境中断层的形成、发展变化。因此通过数值模拟解释现今地壳形变特征时，往往需要结合野外考察、深地震反射、地震精定位、震源机制解等多学科方法给出的断层产状及断裂带断层泥等做出一定程度的简化，进而构建断层物理模型。
在长期的实践应用过程中，根据物理模型及具体问题的研究，这些简化的断层模型可以归纳为：劈节点模型、块体模型、介质弱化模型以及摩擦接触模型（和平等，2011）。劈节点模型本质上是在预设断层位错的情况下模拟断层，因而在同震或震后效应研究中具有极大的优势，得到了廣泛应用；块体模型本质基础是Shi（2007）提出的DDA（Discontinuous Deformation Analysis）算法，它将地球动力过程看作为一系列块体相互作用的过程，断层其实就是块体与块体之间的边界。介质弱化模型和摩擦接触模型作为2种原理简单且便于在一般软件中设置的模型，在区域地壳动力学形变特征模拟及断层活动状态模拟中得到大量应用。由于断层的介质力学参数对区域应力场以及形变场均有一定程度的影响（李玉江等，2010），因此在模拟过程中，介质弱化模型是将断层视为有一定厚度的弱介质层，通常通过降低该介质层弹性模量来模拟断层对区域地壳动力过程的影响（Lu et al， 2011，2012；祝爱玉等，2015a；胡勐乾等，2014；杨兴悦等，2013）；摩擦接触模型则是直接在创建好的模型基础上给断层面赋摩擦系数，进而实现断层两盘的相对运动（Hergert，Heidbach，2010； Zhu，Zhang，2010； He et al，2013；李玉江等，2014；祝爱玉等，2015b）。但也存在一些问题：①对于介质弱化模型来说，其本质上是通过降低断层介质弹性模量来近似模拟断层两盘的非连续变形。尽管野外观测到断层带内有糜棱岩或断层泥等弱力学属性材料填充，但这种弱力学属性的填充物已经在长期的构造作用下与周边介质处于力平衡状态。因此，在构建有限元模型时直接赋予该填充物一个较低的弹性模量，以期获得断层两盘的运动差异和断层附近应力状态的同时，往往忽略了降低弹性模量，也使得该填充层内部、垂直于断层走向的力平衡破坏，在断层内部垂直于断层方向产生新的附加力。这可能影响模拟结果的可靠性。此外，实验室关于断层泥的弱力学属性对断层行为的影响也表现在剪切摩擦行为上（马胜利，1986）。在各向同性弹性理论前提下还可能会出现不合理的负泊松比的情况（董培育，石耀霖，2013）。②摩擦接触模型可以很好地刻画断层两盘相对运动的非连续变形，相比介质弱化模型更接近于真实断层状态。但摩擦接触计算是个强非线性问题，涉及到收敛与否以及计算量等问题（李玉江等，2009），特别是在研究区域内断层较多的情况下，如果将每条断层都构建为摩擦基础模型，收敛问题将变得更加紧迫，为此研究中还可以将重点断层考虑为摩擦接触，而将次要断层处理为弱化带模型（Li et al，2015），这一简化处理提高了计算效率，但又增加了弱化介质模型自身的间题。
因此，如何做到既保证计算效率、避免介质弱化模型可能带来的不真实信息，又能快速合理模拟断层两盘相对运动是地壳动力学数值模拟的一个基础性工作。本文在介质弱化带模型基础上引入正交各向异性理论，用正交各向异性本构方程替代传统弱化带中的线弹性本构方程，实现了弹性模量与剪切模量相互独立；同时考虑到实际建模过程中可能遇到不同断层走向、倾角各不相同的情况，通过坐标转换给出了针对任意走向、倾角的断层本构关系。最后通过数值试验对比分析了介质弱化模型、摩擦接触模型以及正交各向异性断层模型。
1 正交各向异性本构关系
假定断层带为夹在块体之间的类似于层合板的一类特殊介质：其在垂直于断层方向力学性质与周边块体相近，不易产生大挤压或拉张变形；而在平行于断层运动方向上的剪切模量较小，能够产生较大的剪切变形。则根据复合材料力学（沈观林，胡更开，2006）基础知识，断层内应力应变关系在其3个主方向上满足如下关系：式中：ε=[e11 e22 e33 e12 e23 e13]T，σ=[s11；s22 s33s12 S23 S13]T分别为3个主方向的应变和应力分量；x为柔度矩阵；E，，E2，E3，v12，V23，v13，G12，G23，G13分别为其3个主方向的工程弹性常数。K可表示为：
由于本文重点考虑断层介质在走向和倾向的剪切行为，其他方向上的介质力学性质与周边介质一致。为此将式（2）进一步简化为：式中：E，v分别为周边介质弹性模量和泊松比；Gd；P为断层在剪切方向上的剪切模量，独立于周边介质参数E，v。
在实际操作过程中，由于区域断层较多且其走向和倾角各不相同，因此需要将其主方向下的应力应变矩阵旋转变换到地理坐标（x，y，z）下。图1为断层面主方向（1，2，3）以及地理坐标系（x，y，z）相对关系，其中：1，2，3分别为断层走向（面向走向，右手在上盘）、倾向（逆冲为正）以及张拉方向（张为正）；x，y，z分别为地理东、北向以及垂直地表向上的方向。α为断层走向，β为断层倾角。令断层面应变和应力分量在地理坐标系（x，y，z）下分别为：
e=[exx eyy ezz exy eyz exz]T
s=[sxx syy szz sxy syz sxz]（4）
由图1b可知，将地理坐标系下的应力张量绕z轴逆时针旋转，然后再绕主轴1逆时针旋转β便可得到主方向下的应力张量。则：式中：T=A*B*结合式（1）和（4），式（5）可以转化为：式中：F为坐标转换矩阵，可表示为：同理可得：将式（6）和式（7）带入式（1），可得：
F·e=K·F·s（8）即：
e=D·s=F-1·K·F·s（9）式中：D为正交各向异性材料在地理坐标系（x、y、z）下的柔度矩阵，且：式中：柔度矩阵D为21个独立参数的对称矩阵，因此在实际模拟计算过程中，只需修改正交各向异性本够关系柔度矩阵K在其相应主方向下对应滑动方向的剪切模量（式（3）中的Gslip），结合坐标转换矩阵F，代入式（10）便可生成断层有限元模型的柔度矩阵，进而可以构建有限元计算的整体刚度矩阵，实现垂直断层方向无大变形，平行断层滑动方向有大变形滑动的断层运动状态。为方便控制，本文令断层走向和倾向剪切模量与周边弹性介质剪切模量之比为λ∈[0，1]，即Gslip=λ
2 数值试验对比研究
2.1 模型构建
为检验正交各向异性本构关系模拟断层的可行性，本文设计1个虚拟理想弹性地壳模型，模型为边长L=6000m的正方形，纵向（z）深2000m且上下地殼各1000m厚，模型中部沿x轴走向有一倾角90°的右旋走滑断裂带，长2000m，宽50m。模型施加狄里克莱边界条件：-y边界约束y方向位移为0，+y边界只有y方向位移，Uy=-0.1m；-x边界上施加一反正切函数的x方向剪切位移，如图2a所示；模型-z底面约束z方向位移为0。假定模型弹性模量E=84GPa，泊松比v=0.25。然后针对常用的介质弱化模型，摩擦接触模型以及本文正交各向异性模型分别设计2组工况：第一组工况（Case 1）假定上地壳断层闭锁，下地壳可以剪切滑动；第二组工况（Case 2）认为整个断层处于非闭锁状态，断层上下均可以相对自由滑动，如图2b所示。
2.2 结果对比
在一个地震轮回过程中，断层在地震发生时由于快速滑动而破裂，使得其力学属性快速降低；在震后经历长时间的愈合，在下一次地震来临前，断层力学属性又增加到相对较强的状态。根据地震轮回的理论模式（Savage，Burford，1973；Scholz，1998），借鉴Hergert和Heidbach（2010）在马尔马拉海的工作思路，本文首先用介质弱化模型、摩擦接触模型以及本文给出的正交各向异性本构模型分别模拟一组震间地壳运动位移场图像。假定上地壳己经经历了长久的愈合，断层力学属性与周边地壳介质基本一致，而深部断层强度低，依然可以自由剪切，断层模型参数设置如图2b（Case 1），模拟结果如图3所示。
基于3种断层模型地壳运动位移场图像的空间分布形态基本一致，但一般文献多用的介质弱化模型在断层附近局部，x方向的位移场存在一定程度的畸变，但并未影响Ux沿y的反正切渐变分布特征。另外在空间分布图像上，基于介质弱化断层模型的Ux和基于摩擦接触模型及正交各向异性模型的Ux结果在定量上有所差异。而基于摩擦接触断层模型和各向异性断层模型的对比结果显示，2者在位移场图像上高度一致，这与本文构建模型的初衷一致，因此在考虑断层较多，摩擦接触建模复杂，计算收敛困难的时候可以考虑用各向异性本构关系来模拟断层，其模型构建与计算量与常用的介质弱化模型完全一致，计算中也不会涉及到非线性收敛困难等情况，不同的是有限元方程中的柔度矩阵（或刚度矩阵）内的具体数值有所不同。
震后短期内，断层没有完全愈合，其介质力学属性依然较低。在分析这些介质力学属性本来偏低的断层时，本文将整个断层的力学参数降低，具体如图2b（Case 2）所示。结果显示基于3种断层模型地壳运动位移场图像空间分布不同程度出现差异，常用的介质弱化模型在水平和垂直位移场上明显与其他2种模拟不一致，在断层附近位移场图像出现巨大的畸变，这与断层弹性模量的降低引起垂直于断层方向发生大幅度的挤压变形有关。实际上由于长期的构造作用，断层两盘的相对运动一般只发生在其走向和倾向滑动方向，并不会在断层上存在这样垂直于断层的大幅挤压变形，因此在应用简单的介质弱化断层模型模拟地壳形变时，随着弹性模量的降低，断层附近的模拟位移场图像将出现不同程度畸变，进而可能会影响到断层附近的模拟应力状态（董培育，石耀霖，2013）等。对比基于摩擦接触模型和本文各向异性模型的模拟结果，尽管2者在断层附近有一定程度的差异，但整体上2者空间分布形态完全相似，可以清晰地给出断层两盘的相对滑动，所不同的是由于摩擦接触模型给断层赋予的是摩擦系数，而本文只是近似给的剪切模量，2者虽然在阻碍变形的机制上相似，但其背后的物理方程不同，因此造成2者在断层两盘相对滑动量的差异。后续的研究需要基于地壳运动位移场图像讨论2者之间的统计关系，这样会进一步提高本文各向异性断层模型在模拟中的应用和理解。
为进一步分析3类断层模型的异同，从图3和图4中切出x=0剖面上的Ux和Uy做对比分析，如图5所示。结果显示当在Case 1（浅部断层愈合，深部断层可以滑动）情况下，总体上这3类断层模型给出的结果基本一致，但从更细致的角度看，基于摩擦接触断层模型和正交各向异性断层模型给出的结果吻合一致。但介质弱化模型给出的结果与其他2种断层模型给出的模拟结果略有差异，并且这种差异在空间上呈现近断层处相对较小，远离断层处增大。这是由于浅部断层与周边介质属性一致，保证了表层形变场的连续性，但基于介质弱化的深部断层模型破坏了垂直于断层方向的平衡状态，產生附加的形变场畸变信息。进一步降低表层断层的介质力学属性（即Case 2），这种畸变信息对表层形变场的影响更加显著，特别是在垂直于断层方向上（图Sd）。但在平行于断层方向上，3类断层模型给出的模拟结果均显示出不同程度的差异，这是因为本文针对3类断层模型的力学参数都是随机赋值，没有探讨3者力学参数之间的等效关系。但3者均可以表征本文预设的断层走滑行为。对于垂直于断层方向，图Sd显示基于正交各向异性理论的断层模型与摩擦接触断层模型的模拟结果一致，表明正交各向异性断层模型相比介质弱化模型，更接近于实际情况。
为进一步分析不同断层模型对断层附近模拟应力状态的影响，考虑到图2设计的模型表现为剪切变形行为，以水平剪应力（τyx）的空间分布为讨论对象。在图2的基础上继续分析了不同1况下，不同断层模型对应的τyx模拟结果，如图6所示。结果显示：较介质弱化模型给出的τyx空间分布（图6a，d），正交各向异性断层模型的模拟结果（图6b，e）与摩擦接触模型的结果（图6c，f）更为相近。在上地壳断层闭锁的情况下（Case 1），水平剪应力τyx云图在空间上相对y=0轴对称（图6b、c），而介质弱化模型给出的结果不同于其他两个模型的模拟结果，水平剪应力呈现出不对称形态（图6a）；在整个断层均贯通非闭锁的情况下（Case 2），尽管正交各向异性断层模型的模拟结果（图6e）与摩擦接触模型（图6f）的模拟结果有一定的差异，但2者基本保持相似的空间分布特征，而基于介质弱化模型给出的模拟结果（图6d）则在空间分布上表现出不规则的分布形态。这些都是由于在数值模拟中，介质弱化模型由于其本构关系的原因，不可避免地在垂直于断层方向上引入了较大的挤压变形（图4d），这种不合理的挤压变形不但影响着位移场的空间分布特征，同时也影响着应力场的空间分布。
3 讨论
近年来，随着GPS地壳速度场观测数据的丰富，基于GPS约束的地壳动力学数值模拟研究日渐深入（Hergert，Heidbach，2010；He et al，2013；Lei et al，2010； Cho，Kuwahara，2013；龙小刚，朱守彪，2015；刘峡等，2010）。一般情况下，研究人员通常都会通过一些途径选定有限元模型的边界条件，模拟计算给出模型内部的速度场模拟结果，然后通过对比模拟结果与观测结果之间的残差来确定模型是否可靠，参数选择是否合理，在大量比较之后给出合理模型的基础上分析区域的动力学机制及应力状态等，这就要求我们首先要保证所采取的模型是合理的。介质弱化模型作为最为简单的断层模型，其建模方便、不涉及非线性问题且不会存在难于收敛的情况，在一些研究中得到大量的应用。因此在应用介质弱化模型的情况下，在对比模拟地壳运动速度场与GPS观测速度场残差的同时，如何评价类似于图4a、d以及图5d的形变场畸变对模拟结果的影响也是不容忽视的问题。在不可避免使用介质弱化断层模型的时候，如能进一步修改其有限元刚度矩阵，使其剪切模量与弹性模量相互独立，对于提高模拟结果的可靠性具有一定的意义。
4 结论
本文在介质弱化带模型基础上，引入正交各向异性理论，同时考虑到实际建模过程中可能遇到不同断层的走向、倾角各不相同的情况，本文通过坐标转换给出了针对任意走向、倾角的正交各向异性断层本构关系。最后通过数值试验对比分析了介质弱化模型、摩擦接触模型以及正交各向异性断层模型在震间地壳速度长图像模拟中的差异，结果显示：（1）无论断层力学性质强弱，传统的介质弱化模型都会对地表地壳运动位移场产生一定程度的影响，断层力学性质（对应于图2不同断层模型中断层的力学参数，如弱化介质模型的弹性模量、各向异性模型的剪切模量以及摩擦接触模型的摩擦系数）越弱，影响越显著，这主要由于是降低断层弹性模量，使得断层内垂直于断层方向的力平衡被打破，进而影响了断层附近位移场图像的空间分布，这在模拟中同时也会对断层附近的应力场和应变场产生一定程度的影响：（2）从位移场图像以及应力场图像的空间分布来看，正交各向异性本构关系与摩擦接触模型给出的模拟结果基本一致，因此在区域构造复杂的情况下，可以考虑使用各向异性本构关系来表征断层的变形机制。
非常感谢陈连旺研究员的悉心指导，同时也感谢中国地震局访问学者计划的支持和地壳应力研究所提供良好的办公环境和软件支持！部分图件采用GMT（Wessel et al，2013）绘制，特此表示感谢！
参考文献：
陈连旺，陆远忠，张杰，等.1999.华北地区三维构造应力场[J].地震学报，21（2）：140-149.
陈连旺，张培震，陆远忠，等.2008.川滇地区强震序列库仑破裂应力加卸载效应的数值模拟[J].地球物理学报，51（5）：1411-1421.
董培育，石耀霖.2013.关于“用单元降刚法探索中国大陆强震远距离跳迁及主体活动区域转移”的讨论——横向各向同性“杀伤单元”才是更好的途径[J].地球物理学报，56（6）：2133-2139.和平，李志雄，陆远忠，等.2011.断层面的有限单元模拟方法综述[J].西北地震学报，33（2）：186-194.
胡勐乾，邓志辉，陆远忠，等.2014.三维数值模拟在华北地区现今构造变形分析中的应用[J].地震地质，36（1）：148-165.
李玉江，陳连旺，刘少峰，等.2014.芦山地震的发生对周围断层影响的数值模拟[J].地球学报，35（5）：627-634.
李玉江，陈连旺，叶际阳，等.2010.岩石物性变化对区域应力场的影响[J].地球物理学进展，25（6）：1941-1946：
李玉江，陈连旺，叶际阳.2009.数值模拟方法在应力场演化及地震科学中的研究进展[J]地球物理学进展，24（2）：418-431.
刘峡，马瑾，傅容珊，等.2010.华北地区现今地壳运动动力学初步研究[J].地球物理学报，53（6）：1418-1427.
龙小刚，朱守彪.2015.台湾碰撞带现今地壳变形场特征及其动力学成因的有限单元法模拟研究[J].地球物理学报，58（7）：2350-2365：
马胜利，本利彦.1995.蒙脱石的脱水作用对断层摩擦本构行为的影响[J].地震地质，17（4）：289-304.
马胜利.1986.模拟断层带摩擦滑动性状与变形特征[J].中国地震，（2）：74-80.
邵志刚，傅容珊，薛霆虓，等，2008.昆仑山MS8.1级地震震后变形场数值模拟与成因机理探讨[J].地球物理学报，51（3）：805-816.
沈观林，胡更开.2006.复合材料力学[M].北京：清华大学出版社.
杨兴悦，陈连旺，杨立明，等.2013.巴颜喀拉块体强震动力学过程数值模拟[J].地震学报，35（3）：304-314.
祝爱玉，张东宁，蒋长胜，等.2015a.川滇地区地壳应变能密度变化率与强震复发间隔的数值模拟[J].地震地质，37（3）：906-927.
祝爱玉，张东宁，蒋长胜.2015b.安宁河—则木河—小江断裂带应力状态分段特征的数值模拟研究[J].中国科学：地球科学，45（12）：1839-1852.
Cho I，Kuwahara Y.2013.Numerical simulation of crustal deformation u-sing a three-dimensional viscoelastic crustal structure model for theJapanese islands under east-west compression[J].Earth，Planetsand Space，65（9）：1041-1046.
He J，Lu S，Wang W.2013.Three-dimensional mechanical modeling ofthe GPS velocity field around the northeastern Tibetan plateau andsurrounding regions[J].Tectonophysics，584；257-266.
Hergert T，Heidbach 0.2010.Slip-rate variability and distributed de-formation in the Marmara Sea fault system[J].Nature Geoscience，3（2）：132-135.
Lei X，Chen Y，Zhao J，et al.2010.Modelling of current crustal tectonicdeformation in the Chinese Tianshan orogenic belt constrained byGPS observations[J].Journal of Geophysics and Engineering，7（4）：431.
Li S，Moreno M，Rosenau M，et al.2014.Splay fault triggering by greatsubduction earthquakes inferred from finite element models[J].Geo-physical Research Letters，41（2）：385-391.
Li Y，Chen L，Liu S，et al.2015.Coseismic Coulomb stress changescaused by the MW6.9 Yutian earthquake in 2014 and its correlationto the 2008 MW7.2 Yutian earthquake[J].Journal of Asian EarthSciences，105：468-475.
Lu Y，Yang S，Chen L，et al.2011.Mechanism of the spatial distributionand migration of the strong earthquakes in China inferred from nu-merical simulation[J].Journal of Asian Earth Sciences，40（4）：990-1001.
Lu Y，Yang S，Chen L，et al.2012.Migration trend of strong earthquakesin North China from numerical simulations[J].Journal of AsianEarth Sciences，50：116-127.
Moreno M，Rosenau M，0ncken O.2010.2010 Maule earthquake slip cor-relates with pre-seismic locking of Andean subduction zone[J].Nature，467（7312）：198-202.
Savage J C，Burford R O.1973.Geodetic determination of relative platemotion in central California[J].Journal of Geophysical Research，78（5）：832-845.
Scholz C H.1998.Earthquakes and friction laws[J].Nature，391（6662）：37-42.
Shi G.2007.Applications of discontinnons deformation analysis（DDA）to rock engineering[M].Berlin，Heidelberg，Computational Mechan-ics springer；136-147.
Sibson R H.1977.Fault rocks and fault mechanisms[J].Journal of theGeological Society，133（3）：191-213.
Wessel P，Smith W H，Scharroo R，et al.2013.Generic mapping tools：Im-proved version released[J].Eos，Transactions American GeophysicalUnion，94（45）：409-410.
Zhu S，Zhang P.2010.Numeric modeling of the strain accumulation andrelease of the 2008 Wenchuan，Sichuan，China，Earthquake[J].Bul-letin of the Seismological Society of America，100（5B）：2825-2839.
