芦山MS7.0地震同震形变的地形效应

王金驰 段虎荣,2 张成浩 梁文康 刘 鹏

1 西安科技大学测绘科学与技术学院,西安市雁塔中路58号,710054 2 中国科学院精密测量科学与技术创新研究院大地测量与地球动力学国家重点实验室,武汉市徐东大街340号,430077

前人根据地震同震形变反演断层参数时,为了简化模型,常常会忽略真实地形因素。在计算同震形变时地表真实地形会产生多大影响至今仍然不明确,部分学者认为地形效应对同震形变影响较大[1-3],但也有学者持相反意见[4-6]。

2013-04-20四川芦山发生MS7.0地震,震中位于四川盆地西缘。四川盆地平均海拔为250～750 m,而青藏高原平均海拔为3 000～5 000 m,地形梯度变化较大,因此此次地震是研究地形对地震同震形变影响的较好案例。谱元法(spectral element method,SEM)可兼顾有限元法网格的灵活性和伪谱法的高精度、快速收敛特性,在复杂地震动数值模拟方面得到广泛应用[7-9]。本文利用谱元法模拟2013年芦山MS7.0地震同震形变,探讨地形对地表同震形变的影响。

1 谱元法计算同震形变

1.1 控制方程

图1 准半空间中断层示意图

区域Ω的控制方程可表示为:

∇·T+f=0

(1)

式中,T为应力张量,f为外力。

在自由表面Γ0上满足边界条件:

(2)

同震形变由弹性本构关系T=c:ε控制,c为4阶弹性张量,ε为应变张量。

1.2 分裂节点法

地震震源由断层面Σ上的滑动矢量Δs定义。如图1所示,当从断层Σ的正侧向负侧移动时,滑动矢量捕获到位移场s的不连续性[10]:

(3)

本文采用分裂节点法对模型域进行网格划分,使元素边界满足断层面[11](图2)。在分裂节点法中,通过在断层面上规定一个不连续的位移场来引入载荷,因此控制方程(1)中f为0。

图中使用谱元素(浅灰色单元)离散化研究区域Ω,在区域外添加一层无限元素(深灰色单元);黑色粗实线表示断层

2 数值模拟

通过模拟基于弹性半空间各向同性介质的垂直断层走滑和倾滑,验证谱元法在地震同震形变模拟中的可靠性。采用图3所示模型,尺寸为100 km×80 km×30 km,泊松比为0.25,杨氏弹性模量为5.68 GPa。图中深灰色矩形表示断层,其顶部距离地表4 km,底部距离地表14 km,中心距离模型左右两侧边界50 km。使用近似为2 km×2 km×2 km的六面体单元对模型进行网格划分,共包含30 000个单元和253 611个节点,节点自由度为725 370。模型顶面施加应力自由边界条件,侧边和底面施加Dirichlet边界条件。

表1 不同倾角下计算的同震位移

图3 垂直断层模型

2.1 走滑断层引起的同震形变

设定断层右旋纯走滑量为4 m,并使用谱元法计算地表三维同震形变,将数值模拟结果与Okada解析法结果进行对比(图4)。从图4(a)和4(b)可以看出,采用2种方法计算得到的水平位移场沿断层两侧呈反对称分布,靠近断层位置的水平位移较大且与断层走向平行。从图4(c)和4(d)可以看出,采用2种方法计算得到的垂直位移场呈四象限对称分布。水平和垂直位移场的分布均符合右旋走滑断层特征,并且谱元法与Okada解析法结果在数值上具有很好的一致性。

粗黑直线代表断层所在位置,红色箭头代表断层走向

为详细研究谱元法与Okada解析法的差异,选择位移场剖面X=-10 km位置,将2种方法的模拟结果进行对比。从图5可以看出,对于水平位移x分量,除边界区域外,两者的振幅和相位基本一致。对于水平位移y分量,在断层附近,左侧Okada解析法结果大于谱元法,右侧谱元法结果大于Okada解析法;而在超过7 km之后,Okada解析法结果大于谱元法。对于垂直位移z分量,在断层附近,谱元法结果大于Okada解析法。总体而言,两者差异最大不超过10%。

图5 走滑断层在X=-10 km处地表位移剖面

2.2 倾滑断层引起的同震形变

设定断层纯倾滑量为4 m,同样采用上节中所构建的模型区域和介质属性。对比谱元法结果与Okada解析法结果(图6)可知,水平和垂直位移场的分布均符合断层倾滑特征,并且谱元法与Okada解析法结果在数值上具有很好的一致性。

粗黑实线代表断层所在位置,红色圆圈和叉代表断层倾向

选取X=10 km地表位移剖面,对比谱元法数值模拟结果与Okada解析法结果。由图7可知,x和z分量呈中心对称分布,y分量呈轴对称分布。靠近断层时,2种方法所有分量的振幅和相位吻合较好;远离断层时,2种方法所有分量均存在一定差异,这可能是由建模误差或数值方法的近似性所致。在下一步研究中,将考虑添加无限元素以更好地探究这些差异的原因。

图7 倾滑断层在X=10 km处地表位移剖面

综合走滑运动和倾滑运动数值模拟结果可知,谱元法与Okada解析法计算结果一致性良好。因此,认为使用谱元法模拟地震地表同震形变具有可行性。

3 地形因素对2013年芦山MS7.0地震同震形变的影响

3.1 数 据

本文使用的GPS同震位移资料来自于33个GPS连续站,GPS观测数据最大水平误差小于3 mm,最大垂直误差小于12.4 mm[12]。图8为本文研究区,震源中心位于30.3°N、103.0°E。使用单一平面断层生成格林函数,断层参数参照Duan等[13]的研究成果。

图中沙滩球为CENC CMT震源机制解;红色五角星为震中位置;U-V为谱元法计算时所采用的坐标系,U为沿断层走向方向,V为垂直断层走向方向

将芦山地区真实地表-地壳结构纳入模型中,地形网格文件使用NOAA发布的ETOPO1数据[14],分辨率为1′,包含陆地地形和海洋水深数据。模型介质参数参考谈洪波等[15]的研究成果。

3.2 模型构建

为研究地形因素对同震形变的影响,构建图9中2种模型,其中A为地形模型(图9(a)),B为无地形模型(图9(b))。本文模型基于均匀弹性介质,泊松比为0.25,杨氏弹性模量为79.8 GPa,水平尺寸为116 km×133 km,深度为40 km(不包括地形)。离散模型使用六面体单元,所有模型平均网格间距设置为2 km,模型A共包含79 750个单元和675 608个节点,模型B共包含81 200个单元和686 868个节点。

粉红色矩形为断层位置;U为沿断层走向方向,V为沿垂直断层走向方向

3.3 基于分裂节点法计算地形对芦山地震同震形变的影响

2013年芦山MS7.0地震是发生在龙门山断裂带西南段的一次浅源逆冲型地震,破裂长度20 km[16],地震引起明显的地表水平位移及垂直形变,可为研究地形效应对同震形变的影响提供理想的条件。

本文分别基于模型A和模型B,利用谱元法计算地表同震位移,并探讨地形对地表同震位移的影响。图10为2种模型生成的地表同震位移场,可以看出,模型A计算的沿断层走向方向位移为-71～0 mm,垂直断层走向方向位移为-18～28 mm,垂直方向位移为-42～67 mm;模型B则分别为-80～2.5 mm、-19～31 mm和-47～75 mm。由图9可知,同震形变场分布形态符合逆冲型地震特征,距离震中越近,同震形变越大。地形的存在不仅会改变地表同震形变场分布形态,也会改变同震位移的大小。从图10(a)和10(d)可以看出,在震中上方区域,模型A形变场分布不如模型B平滑,且模型A的最大值比模型B少9 mm,对比图10(b)与10(e)、10(c)与10(f),也有类似的情况。

第1、2行分别为基于模型A、B计算的同震位移,第3行为模型A与模型B之差;第1、2、3列分别代表沿断层走向的水平位移、垂直断层走向的水平位移、垂直位移

本文利用林晓光等[3]的方法计算地形效应的影响。结果表明,沿断层走向方向、垂直断层走向方向和垂直方向,地形对同震形变的影响最大可达约16.8%、19.2%和18.3%。

4 讨 论

本文基于平面断层模型研究地形对同震形变的影响,断层参数在一定程度上会影响地表同震位移,这里主要讨论断层倾角和深度对地表同震位移的影响。

在Duan等[13]的研究结果基础上,构建6种模型,采用控制变量法研究断层倾角和深度对地表同震位移的影响,所有模型使用与§3.2相同的介质参数。数值模拟结果如表1、2所示。

根据所选研究区域,使用8个GPS站点观测值,为了与谱元法结果进行对比,将这些GPS观测值转换到U-V坐标系下(谱元法计算时所采用的坐标系),沿断层走向方向U分量变化范围为-63.8～2.4 mm,沿垂直断层走向方向V分量变化范围为-17.6～22.1 mm,垂直方向Z分量变化范围为-4.9～83.6 mm。

通过表1可以看出,随着断层倾角增大,沿断层走向方向分量、垂直断层走向方向分量和垂直方向分量位移均减小。

由表2可知,随着断层深度变浅,地表同震位移逐渐变大。这是因为浅层断层距离地表更近,断层运动能更直接地影响地表,从而使地表同震形变更为明显。

表2 不同深度下计算的同震位移

在原始断层几何参数基础上,将断层走滑分量和倾滑分量分别设置为0.15 m和0.70 m,计算得到地形效应对沿断层走向方向同震形变的影响最大可达约10.3%,对垂直断层走向方向同震形变的影响最大可达约12.4%。这与Lin等[2]计算汶川地震地形效应对同震位移的影响(9.05%)不同,可能是由于其使用的是二维有限元模型,而本文采用的是三维谱元法模型。Langer等[1]计算智利地震地形效应对同震位移的影响为30%,远大于本文计算值,这可能是由于秘鲁-智利海沟地形急剧变化,相对高差过大,而且断层两侧介质不同,一侧为海域,而另一侧为陆地。

值得一提的是,GPS结果是某一确定点的形变数据,而谱元法结果是整个网格区域的形变情况。为了方便二者对比,谱元法结果使用与GPS观测点邻近区域的位移值。本文选取同震区域位移较大的LS05、LS07、SCTQ测站进行对比。LS05测站东、北分量观测值分别为-9.9 mm、-66.8 mm,谱元法模拟值分别为-10.4 mm、-53.3 mm;LS07测站东、北分量观测值分别为24.1 mm、-17.8 mm,谱元法模拟值分别为17.1 mm、-17.0 mm;SCTQ测站东、北分量观测值分别为-9.8 mm、-18.7 mm,谱元法模拟值分别为-7.0 mm、-20.1 mm。除LS05测站北分量差异较大外(13.5 mm),其他结果二者较为接近。

5 结 语

在不考虑地形条件下,采用谱元法与Okada解析法分别计算垂直断层走滑、倾滑运动引起的地表同震位移。结果表明,无论对于走滑断层还是倾滑断层,2种方法结果的一致性均较好,表明谱元法模拟同震形变具有可靠性。

针对2013-04-20芦山MS7.0地震,本文分别构建芦山区域地形模型和无地形模型,探讨地形效应对同震形变的影响。结果表明,地形因素不仅会改变地表位移场分布形态,而且量值在不同方向上也有所差异。对于水平同震形变分量而言,沿断层走向和垂直断层走向的影响最大分别约为10.3%、12.4%;对于垂直同震形变分量而言,最大影响可达约11.9%。反演此次地震断层滑动速率时,若不考虑地形因素的影响,计算结果会低估断层走滑分量和高估倾滑分量。

本文基于平面断层几何模型探究地形对地震同震形变的影响,然而芦山地区断层几何复杂,为了更全面地理解同震形变的影响,后续将进一步收集同震DInSAR数据,并考虑更精细的断层几何模型,以便更准确地研究地形对同震形变的影响。另外,对于复杂的断层几何模型,不同构建方法所得结果的差异也值得进一步探索。

致谢：感谢Gharti教授提供SPECFEM-X软件。