引言

隐伏活动性断层在中国广泛分布(Xu等,2016)。活断层不仅是产生地震的根源,地震时沿断层线两侧地区的破坏也最为严重,因而,研究隐伏断层错动时断层带两侧上覆土体的破裂变形特征及影响因素对于地震灾害的预防具有重要意义。汶川地震中断层错动引起的地表破裂如图 1所示(Li等,2015)。


图 1 汶川地震诱发的地表破裂 Fig. 1 The surface rupture caused by Wenchuan Earthquake

对活动断层错动导致地表破裂的研究,目前主要采用野外调查、构造物理模型试验和数值模拟等方法。在地震断层突然错动导致的地表位移场计算方面,主要采用2种方法:一种是解析方法,另一种是有限元法(徐锡伟等,2011)。 Okada(1992)针对断层错动过程中由于拉伸和剪切造成的地表土体的内部位移和变形,提出了一套半无限空间弹性体的完整解析表达式,并应用表达式进行了数值计算。韩竹军等(2002)采用Okada地表位移的计算方法,结合《建筑抗震设计规范》(GB 50011—2001)(中华人民共和国建设部等,2001)给出了隐伏断层错动引发地表破裂的临界值,讨论了地表破裂带宽度随埋深、倾角和断面错动量的关系。赵雷等(2007)对含有软土夹层的上覆土层在断层错动过程中的破裂特征进行了数值模拟研究,发现软弱夹层所处的位置及薄厚程度对破裂进程有一定影响。Loukidis等(2009)对倾滑断层错动过程中覆盖土层的破裂过程开展了二维数值模拟研究,结果显示:从基岩到地表的破裂扩展偏离基岩断层线在地表的投影,对于非黏性土和小断层倾角可能形成地堑;对于正断层,破裂扩展至地表需要的断层位错量为覆盖土层厚度的0.2%—0.4%。赵颖等(2013)采用拟静力弹塑性有限元方法,对走滑断层作用下粉质黏土和黏土覆盖层的地表破裂进行了数值模拟分析,结果表明走滑断层引发的地表破裂不仅与震级有关,还与土层厚度和土层性质有关。崔光耀等(2018)对位于黏滑逆断层上方隧道的减小错动效果进行了分析,研究表明断层的黏滑错动对位于断层上盘的隧道影响远大于下盘。

以往的研究对不同错动方式下断层错动的地表破裂响应缺乏必要的分析,覆盖土层的破裂标准没有明确的指标或处于初步讨论阶段。为深入理解正断层错动导致的上覆土体破裂特征,本文利用大型通用有限元软件ADINA建立了三维数值模型,结合覆盖土层破坏时的力学机制和破坏标准,对不同错动方式的正断层错动致地表破裂效应进行研究,并讨论了地震断层地表破裂与强变形带的分布特征及影响因素。

1 有限元模型建立
1.1 正断层错动方式

野外调查表明,正断层的倾角一般大于45°,沿断层走向存在均匀错动、倾斜错动和翘倾错动等复杂的错动方式(图 2)。均匀错动是断层上盘块体沿断层倾角向下整体错动(图 2(a));倾斜错动是指在上盘块体上施加的位移大小沿断层走向呈线性变化,但沿断层倾向的位移大小不变(图 2(b));翘倾错动是指在上盘块体上所施加的位移大小沿断层倾向呈线性变化,但沿断层走向的位移大小不变(图 2(c))。考虑到7级以上地震断层的错动量往往大于2m,因此在数值模拟中所施加的最大垂向位移量为3m。图 2z为施加的垂向位移,h为施加的水平位移,β为断层倾角。


图 2 正断层错动的加载方式 Fig. 2 The loading patterns of normal fault
1.2 有限元计算模型

考虑到活断层的分段性(丁国瑜,1995)和模型单元的数量,选取模型长为5000m,宽为4000m,基岩厚度400m,上覆土体的厚度分别为20m、50m、100m和150m进行模拟计算,隐伏断层走向沿模型长5000m的方向(图 3)。根据钻孔探测结果,在上覆松散土层中,波速变化比较均匀,物理力学性质较均匀,可以将覆盖土体看作1层(徐锡伟等,2011),因而在建立几何模型时,将模型的地层划分为下伏基岩和上覆土体2部分。上覆土体为均一土层,采用摩尔库伦岩土材料模型进行模拟,下伏基岩中含有断层带,下伏基岩和断层带均采用线弹性材料模型,断层带的模拟参考Hernandez等(2012)所使用的模拟方法。研究表明,在隐伏活动断层发生突然错动时,介质的物理力学参数所能引起的地表位移量的差值不超过1‰(徐锡伟等,2011),因而可以采用一组岩土参数来研究地表变形,所选取的模型参数见表 1


图 3 有限元模型 Fig. 3 Finite element model
表 1 模型参数 Table 1 Model parameters

为便于计算分析,选择垂直断层走向的剖面线AA'BB'CC',它们距断层一端的位置分别为100m、2500m和4900m。模型共划分单元68262个,节点54332个,上覆土层和下伏基岩采用粘结网格法联接在一起。为提高断层带附近上覆土体的计算精度,在该位置采用网格加密技术。为消除边界影响,上覆土体的地表设置为无约束自由面,模型周边以及底部施加无限远边界,无限远边界有粘结网格法和弹簧单元法(马野等,2011),本文采用弹簧单元法定义粘弹性的接地弹簧边界(图 3),通过在断层上盘块体施加位移荷载实现断层的错动。

2 地表土体的应力路径变化特征

土体的变形与强度不仅与受力大小和应力历史有关,还与土体应力路径有关。土体的应力路径可以模拟土体的实际应力变化,有助于全面研究应力变化过程中土体的力学性质。通过对不同上覆土体厚度(20m、50m、100m和150m)、不同断层倾角(60°、75°和85°)、不同错动方式(均匀、倾斜和翘倾)的有限元模型进行计算,获得了地表各位置的平均应力p和偏应力q,以此描述地表土体的应力路径变化特征。为便于描述,以断层上断点在地表的投影为界限,划分为断层上盘一侧和断层下盘一侧的土体。限于篇幅,以上覆土体厚度为100m的模型计算结果为例进行展示。

为便于分析地表强变形带的应力路径特征,沿剖面BB′选取地表不同位置的节点,各节点位置以距断层带的距离来表示,如表 2所示。表中距离所对应的正值表示节点位于断层上盘一侧,负值表示节点位于断层下盘一侧。选取节点号1、8、16、24、32和40描述断层带两侧上覆土体的应力变化特征,其中节点24为断层带对应的上覆土体位置。各节点的应力由该节点所在单元的高斯点应力采用RST样条插值法获得。

表 2 节点的位置及编号 Table 2 The location and number of the nodes

随着断层上盘块体的错动,地表土体的应力路径变化过程如图 4所示。由图 4(a)可知,在均匀错动加载方式下,断层下盘一侧地表土体(节点1、8和16)的应力路径变化近似于三轴拉伸状态,应力路径斜向左上方发展,直至达到临界状态线。位于断层带上部的地表土体(节点24)的应力路径开始为斜向右上方发展,类似于三轴压缩状态,但发展到一定程度时,其应力路径转折向临界状态线;对于断层倾角85°的情况,地表土体应力路径(节点24)未达到临界状态线,而对于断层倾角75°、60°的情况,地表土体的应力路径达到了破坏标准。这说明倾角小的断层,其断层带上部的地表土体更容易破坏。断层上盘一侧地表土体(节点32和40)的应力路径先向右下方呈卸载状态,然后出现转折,向右上方发展,直至达到破坏。所不同的是,对较小倾角(如60°)的断层,距离断层较远的地表土体(节点40)未达到临近状态线。这说明断层倾角越小,上盘一侧地表破坏土体的范围向断层带收缩。


图 4 地表不同位置土体应力路径的变化过程 Fig. 4 The changes of stress path of soil mass on the ground surface

图 4(a)图 4(b)可知,在倾斜错动与均匀错动方式下,地表土体的应力路径变化总体趋势基本一致,但对于断层带上部的地表土体(节点24)的应力路径有差别。在倾斜错动方式下,应力路径开始也向斜右上方发展,类似于三轴压缩状态,但随后转折的方向垂直于横坐标,类似于纯剪状态。由图 4(c)可见,与上述2种错动方式相比,应力路径的变化总体趋势一致,但当断层倾角较小时(如60°),上盘一侧距离断层带较远的地表土体(节点40)达到了临界状态线。

上述研究表明,在断层错动过程中,断层下盘一侧地表土体的应力路径呈三轴拉伸变化状态,最先达到破坏;断层上盘一侧地表土体的应力路径变化为:先卸载,再转折为三轴压缩状态,随后向临界状态线靠近,直至达到破坏;断层带在地表投影位置附近土体的应力路径变化比较复杂,既受错动方式的影响,又受断层倾角大小的影响。

不同错动方式下地表土体产生塑性破坏的形式见图 5,图中深色部分为上覆土体出现的塑性破坏区,3种错动方式下地表土体的塑性区均平行于断层走向,然而均匀错动和翘倾错动方式下的地表塑性区长度与隐伏断层长度相同,而倾斜错动方式下地表土体的塑性区长度短于断层长度,且破坏宽度有变化。


图 5 地表土体的塑性破坏形式 Fig. 5 The plastic failure features of the soil mass on the ground surface
3 地表破裂与强变形带
3.1 变形特征

地震灾害具有沿发震断层带呈狭窄带状分布的特征,因而研究断层带(面)两侧上覆土体的同震地表错动和强变形带对于抗震防灾具有意义。韩竹军等(2002)利用位移矢量差计算了地表的变形,计算方法如下:

水平位移差:

$ \Delta {{h}_{i}}={{h}_{i+1}}-{{h}_{i}} $ (1)

垂直位移差:

$ \Delta {{z}_{i}}={{z}_{i+1}}-{{z}_{i}} $ (2)

总位移差:

$ \Delta {{d}_{i}}=\sqrt{\Delta h_{i}^{2}+\Delta z_{i}^{2}} $ (3)

Δdi的空间位置:

$ {{d}_{j}}={({{d}_{i}}+{{d}_{i+1}})}/{2}\; $ (4)

其中,hihi+1为空间相邻2点的水平位移;zizi+1为空间相邻2点的垂直位移;didi+1为空间相邻2点的位置。如果将相邻2点的距离作为临近建筑物的跨度L,则可以预测该跨度下建筑物可能遭遇的地震错动变形量,为抗震设防提供依据。由于地表变形依赖于两相邻点的跨度,因而可以用垂直位移比Δzi/L、总位移比Δdi/L来评价地表变形的程度。

地表变形的强烈地方则建筑物受损严重,然而由于隐伏断层错动过程中上覆土体的应力发生调整,地表变形强烈的地方,其地表土体未必发生塑性破坏。本文采用摩尔库伦土体模型,以是否出现塑性变形来表征地表土体是否发生破坏,常用的塑性变形为等效塑性应变(张学言等,2004),其表达式如下:

$ \varepsilon _{\text{ep}}^{p}=\sqrt{\frac{2}{9}\left[ {{(\varepsilon _{xx}^{p}-\varepsilon _{yy}^{p})}^{2}}+{{(\varepsilon _{xx}^{p}-\varepsilon _{zz}^{p})}^{2}}+{{(\varepsilon _{yy}^{p}-\varepsilon _{zz}^{p})}^{2}}+6({{(\varepsilon _{xy}^{p})}^{2}}+{{(\varepsilon _{xz}^{p})}^{2}}+{{(\varepsilon _{yz}^{p})}^{2}}) \right]} $ (5)

其中,εxxpεyypεzzp分别为xyz方向的塑性应变;εxypεxzpεyzp分别为xyxzyz面上的塑性剪应变。沿剖面BB′(图 3)选取地表不同位置的44个节点(表 2),通过计算获得各位置的垂直位移比、总位移比和等效塑性应变。为方便展示,在断层倾角一定的情况下,选取上盘块体垂直错动分别为0.6m、1.8m和3.0m时的垂直位移比Δzi/L、总位移比Δdi/L和等效塑性应变εep在断层带两侧的变化情况,见图 6。由图可知,在不同断层错动方式下,垂直位移比Δzi/L和总位移比Δdi/L的曲线呈中部高两侧低的形态,即断层带附近的数值大,向两侧逐步减小。然而,等效塑性应变在断层带附近为0,两侧的数值大,说明地表强变形带与地表破裂带并不一致,因为等效塑性应变大于0,土体才产生塑性破坏,然而在地表变形最大的部位,等效塑性应变为0,从塑性力学的角度并没有产生破坏或破裂。出现上述不一致的原因可解释为:由于断层错动,上覆土体在断层带附近形成1个陡坡,坡体的顶部在下盘一侧,容易产生拉伸破坏,坡体的底部在上盘一侧,剪应力集中,容易产生剪切破坏,坡体的中部也是断层带附近的土体,相对较为稳定,但此处的总位移比Δdi/L最大,地表变形程度大,对建筑物安全的影响强。图 6还表明,在断层带两侧较远位置(大于75m)的地表土体,等效塑性应变与总位移比的变化趋势基本一致,均为距断层带越近,其值越大,因而可以综合考虑等效塑性应变和总位移比2个指标,联合确定隐伏断层错动对地表以及建筑物的影响。


图 6 断层错动过程中沿断层带上部地表土体的变形情况 Fig. 6 The deformations at the ground surface above the fault zone during fault movement

上述研究表明,从建筑物受损的角度考虑,地表强变形带和地表破裂带均具有重要的意义,但是两者并非同1个问题。断层带在地表投影一定范围内的土体,地表变形量最大,但此处土体未发生塑性破坏;在断层带两侧较远位置(大于75m)的土体,地表变形程度和塑性应变的变化趋势一致。

3.2 判别方法

地震地表破裂带和强变形带对于工程建设的抗震设防至关重要。徐锡伟等(2002)通过统计学方法研究了地震地表破裂带,将地震地表破裂带作为确定活断层“避让带”的重要依据。然而对于隐伏活断层,上覆松散土层的破裂是1个复杂的力学问题,土体变形既有可恢复弹性变形,又有不可恢复的塑性变形(永久变形),庄妍等(2016)在对路面材料分析时采用等效塑性应变达到0.2%时的σ值作为材料的屈服强度,如果按照此标准判别地表土体的破裂,则与徐锡伟等(2002)提出的地震地表破裂带的位置不一致,因为从力学角度判别地表土体的破裂和从工程建设受损程度考虑“避让带”是不一致的。从力学角度,地表土体的破坏可能会影响建筑物的安全运行,但并不一定需要避让。然而,地表变形强烈的部位不一定产生塑性破坏(图 6),却对工程建设破坏力极大。因此,“避让带”应该是1个上覆土体发生强烈变形且在某些部位产生破裂的地带。

韩竹军等(2002)将上覆土层作为弹性体,提出了地表破裂带的判别方法,获得了隐伏活断层突然错动产生地表破裂带的临界值,即总位移比Δdi/L为0.02。采用相同的模型(图 3)将上覆土体改为弹性体,密度、变形模量和泊松比与表 1土层的参数相同,并与摩尔库伦弹塑性体的计算结果进行对比分析,地表变形的分布情况见图 7。由图可知,当断层错动位移达到3.0m时,距断层带两侧100m范围内,地表变形较大,峰值出现在上盘一侧;随着土层厚度增大,地表变形逐步变小。当土层为弹性体时,地表强变形带曲线呈“帽”状分布,地表变形的峰值位置基本不随土层厚度而变化,位于距断层带约35m的上盘一侧;而当土层为摩尔库伦体时,地表变形的幅度增大,曲线形状不规则,但地表变形的峰值位置向断层上盘一侧方向迁移。当土层厚度达到150 m时,土层作为弹性或摩尔库伦体对于地表变形的影响不大。上述研究表明,对于厚度较小的土层,考虑土层的塑性,强震将造成距离断层带较远的上盘一侧建筑物出现更大的破坏性。在实际强震调查中也发现,隐伏正断层上盘一侧的建筑物出现的破坏更大。


图 7 不同厚度、不同性质土层的地表变形特征 Fig. 7 The deformation features of the surface soil layer with different thickness and properties

除了上述地震地表破裂带以外,工程建设特别关注的地带,即建筑地基出现较大地表变形或出现小裂缝的地带,对于工程建设的安全运行较为重要。根据《建筑地基基础设计规范》(GB 50007—2011)(中华人民共和国住房和城乡建设部,2011),工业与民用建筑相邻柱基的沉降差与相邻柱基的中心距离的比值不能超过0.005,可以考虑将垂直位移比Δzi/L≥0.005作为工程建设特别”关注带”的标准。地表出现小裂缝的地带,即是地表土体达到了破坏且出现了一定量的塑性变形,根据庄妍等(2016)的研究,可以采用等效塑性应变达到0.2%作为工程建设特别关注带的标准。综合考虑地表变形大小和小裂缝,工程建设特别关注带(后简称“关注带”)的范围包含2个部分,即垂直位移比≥0.005的地段和等效塑性应变≥0.2%的地段。

4 “避让带”与“关注带”的宽度及起始位置的比较分析

在不同土层厚度、不同倾角、不同土层性质和不同错动方式条件下,对隐伏正断层错动导致地表变形的计算结果进行了分析,对工程建设“避让带”、工程建设特别“关注带”的宽度、起始位置进行了研究,得到了“避让带”、“关注带”的变化特征及影响因素。

4.1 “避让带”变化特征及影响因素

当断层错动位移达到3.0m时,计算获得沿剖面BB'的“避让带”宽度(总位移比≥0.02)、“避让带”起始位置的变化曲线,见图 8。由图 8(a)可见,断层错动导致地表破裂的宽度在10—90m范围内变化,与上覆土体厚度的关系最为密切,总体趋势随上覆土体厚度增大、“避让带”宽度减小,当厚度达到100m时“避让带”的宽度已经很小。总体上看,在其它条件相同的情况下,断层倾角的减小使“避让带”宽度呈现增大的趋势,摩尔库伦模型相比弹性模型计算得出的上覆土体“避让带”宽度大,均匀错动和翘倾错动对“避让带”宽度的影响较小。由图 8(b)可见,“避让带”的起始位置受断层倾角的影响大,随着断层倾角、上覆土体厚度的增大,“避让带”的起始位置总体上均由断层下盘一侧向上盘一侧移动。在其它条件相同的情况下,当上覆土体为摩尔库伦体时,与弹性体相比,“避让带”起始位置向下盘一侧方向移动。


图 8 “避让带”的宽度和起始位置变化特征(1表示弹性体,2表示摩尔库伦弹塑性体) Fig. 8 The features of width and starting position the avoiding zone

上述研究表明,上覆土体的厚度、隐伏断层倾角对工程建设“避让带”宽度和起始位置的影响最大,其次是上覆土体的弹、弹塑性质。

4.2 “关注带”变化特征及影响因素

当断层错动位移达到3.0m时,沿剖面BB'计算得到的“关注带”(垂直位移比≥0.005或等效塑性应变≥0.2%的地段)宽度、起始位置的变化特征,见图 9。由图 9(a)可知,“关注带”的宽度受上覆土体性质的影响大,当上覆土体作为弹性体时,工程建设“关注带”的宽度较小,在150—200m范围内变化,随上覆土体厚度的增大而增大;当上覆土体为摩尔库伦弹塑性体时,“关注带”的宽度明显增大,在350—400m范围内变化。由图 9(b)可见,“关注带”的起始位置均位于下盘一侧,向上盘一侧发展;“关注带”起始位置受土体性质的影响较大,当上覆土体作为弹性体时,起始位置位于下盘一侧80—120m范围内。而作为摩尔库伦体时,起始位置则向下盘一侧移动到240—280m范围内。


图 9 工程建设“关注带”的宽度和起始位置变化特征(1表示弹性体,2表示摩尔库伦弹塑性体) Fig. 9 The features of width and starting positions of attention zone for engineering construction

上述研究表明,工程建设“关注带”的宽度和起始位置受上覆土体的弹、弹塑性质影响最大,其次为上覆土体的厚度。

5 结论

通过分析不同上覆土体厚度、不同断层倾角、不同错动方式、不同土体性质的断层带两侧地表土体的应力路径、地表变形参数的变化特征及影响因素,得到以下结论:

(1)在隐伏正断层错动过程中,断层两侧地表土体的应力路径变化特征为:位于断层下盘一侧的地表土体,其应力路径呈三轴拉伸变化状态,最先达到破坏。位于断层上盘一侧的地表土体,其应力路径变化为:先卸载,再转折为三轴压缩状态,并向临界状态线靠近,直至破坏。断层带在地表投影土体的应力路径变化比较复杂,既受错动方式的影响,也受断层倾角大小的影响。

(2)从建筑物受损的角度考虑,地表强变形带和地表破裂带均具有重要意义。从力学角度判别地表土体的破裂带和从工程建设受损程度考虑的“避让带”是不一致的,需要综合考虑等效塑性应变和总位移比2个指标来评价同震地表错动对建筑物的影响。

(3)在断层错动过程中,位于下盘一侧的地表土体先出现塑性破坏,随着断层错动位移量的增大,上盘一侧地表土体才出现塑性破坏,地表土体变形较大的部位也逐步移动到上盘一侧。因而,在断层错动位移量较大的情况下,上盘一侧一定范围内地表土体的变形强烈,对建筑物破坏最大,与实际情况相符。

(4)根据建筑行业的规范,并考虑地表产生小裂缝的条件,提出了工程建设“关注带”的确定方法,即地表垂直位移比≥0.005或等效塑性应变≥0.2%的地段。工程建设“关注带”的宽度和起始位置受上覆土体性质的影响最大。当上覆土体为弹性体时,“关注带”的宽度在150—200m范围内变化,起始位置位于断层下盘一侧,距断层带的水平距离较小;当上覆土体为摩尔库伦弹塑性体时,“关注带”的宽度在350—400m范围内变化,起始位置位于断层下盘一侧,但距断层带的水平距离较大。

(5)根据工程建设“避让带”和地表破裂临界值,计算获得了“避让带”的宽度和起始位置。当断层垂直错动位移达到3.0m时,工程建设“避让带”宽度在10—90m范围内变化,受上覆土体厚度、断层倾角的影响最大,总体趋势随着上覆土体厚度和断层倾角的增大、“避让带”宽度逐渐变小;“避让带”的起始位置受断层倾角的影响大,总体趋势随断层倾角的增大,从断层下盘一侧向上盘一侧地表土体移动。

参考文献
崔光耀, 伍修刚, 王明年, 等, 2018. 黏滑断层隧道减错措施参数对减错效果的影响分析[J]. 震灾防御技术, 13(3): 502-511.
丁国瑜, 1995. 活断层的分段模型[J]. 地学前缘, 2(1-2): 195-202.
韩竹军, 冉勇康, 徐锡伟, 2002. 隐伏活断层未来地表破裂带宽度与位错量初步研究[J]. 地震地质, 24(4): 484-494. DOI:10.3969/j.issn.0253-4967.2002.04.002
马野, 袁志丹, 曹金凤, 2011. ADINA有限元经典实例分析[M]. 北京: 机械工业出版社, 118-120.
徐锡伟, 于贵华, 马文涛, 等, 2002. 活断层地震地表破裂"避让带"宽度确定的依据与方法[J]. 地震地质, 24(4): 470-483. DOI:10.3969/j.issn.0253-4967.2002.04.001
徐锡伟, 赵伯明, 马胜利, 等, 2011. 活动断层地震灾害预测方法与应用[M]. 北京: 科学出版社, 132-165.
张学言, 闫澍旺, 2004. 岩土塑性力学基础[M]. 天津: 天津大学出版社, 35-36.
赵雷, 李小军, 霍达, 2007. 数值模拟软夹层对基岩上覆土层破裂的影响[J]. 北京工业大学学报, 33(3): 289-292. DOI:10.3969/j.issn.0254-0037.2007.03.013
赵颖, 郭恩栋, 王琼, 等, 2013. 走滑断层地震地表断裂位错估计方法研究[J]. 岩土力学, 34(5): 1403-1408.
中华人民共和国建设部, 国家质量监督检验检疫总局, 2001. 建筑抗震设计规范(GB 50011-2001)[M]. 北京: 中国建筑工业出版社.
中华人民共和国住房和城乡建设部, 中华人民共和国国家质量监督检验检疫总局, 2011, 建筑地基基础设计规范(GB 50007-2011).北京: 中国建筑工业出版社.
庄妍, 王康宇, 2016. 基于Von-Mises屈服准则的结构安定性研究[J]. 地下空间与工程学报, 12(S1): 170-174, 191.
Hernandez-Marin M., BurbeyT. J., 2012. Fault-controlled deformation and stress from pumping-induced groundwater flow[J]. Journal of Hydrology, 428-429: 80-93. DOI:10.1016/j.jhydrol.2012.01.025
Li B., Su J. Y., Ma D. H., et al, 2015. The study on forecasting the surface rupture width under strong earthquake based on information diffusion methodology[J]. Natural Hazards, 75(2): 1871-1882. DOI:10.1007/s11069-014-1403-1
Loukidis D., Bouckovalas G. D., Papadimitriou A. G., et al, 2009. Analysis of fault rupture propagation through uniform soil cover[J]. Soil Dynamics and Earthquake Engineering, 29(11-12): 1389-1404. DOI:10.1016/j.soildyn.2009.04.003
Okada B. Y., 1992. Internal Deformation due to shear and tensile faults in a half-space[J]. Bulletin of the Seismological Society of America, 82(2): 1018-1040.
Xu L. Q., Li S. Z., Cao X. Z., et al, 2016. Holocene intracontinental deformation of the northern North China Plain:Evidence of tectonic ground fissures[J]. Journal of Asian Earth Sciences, 119: 49-64. DOI:10.1016/j.jseaes.2016.01.003