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

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

1 谱元法计算同震形变 1.1 控制方程

在图 1所示设置中求解控制同震形变的准静态方程。计算区域Ω被嵌入在准半空间中，包括自由表面Γ₀和人工边界Γ，半空间包括自由表面Γ₀和Γ_∞。$ \hat{{\boldsymbol{n}}}$为计算区域边界单位法向量，Ω包含单位法向量为$ {\hat{{\boldsymbol{\upsilon}} }}$的断层面Σ。

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

$ \nabla \cdot T+f=0 $

(1)

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

在自由表面Γ₀上满足边界条件：

$ \hat{\boldsymbol{n}} \cdot T=0 $

(2)

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

1.2 分裂节点法

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

$ \Delta \boldsymbol{s}=\boldsymbol{s}^{+}-\boldsymbol{s}^{-}=[\boldsymbol{s}]_{-}^{+} $

(3)

式中，[·]_－⁺表示捕获断层从正侧到负侧封闭量的跃变。对于滑动断层，位移的法向分量在断层上必须连续，即$ \hat{\boldsymbol{\upsilon}} \cdot \Delta \boldsymbol{s}=0$。还可以在断层上引入不连续张力$ \Delta t=[\hat{\boldsymbol{\upsilon}} \cdot T]_{-}^{+}$。然而，对于自发破裂，由于无外部施加的作用力，张力必须连续，即Δt=0。

本文采用分裂节点法对模型域进行网格划分，使元素边界满足断层面^[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边界条件。

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%。

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

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

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

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

3 地形因素对2013年芦山M_S7.0地震同震形变的影响 3.1 数据

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

将芦山地区真实地表-地壳结构纳入模型中，地形网格文件使用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个节点。

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

2013年芦山M_S7.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)，也有类似的情况。

本文利用林晓光等^[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可知，随着断层深度变浅，地表同震位移逐渐变大。这是因为浅层断层距离地表更近，断层运动能更直接地影响地表，从而使地表同震形变更为明显。

在原始断层几何参数基础上，将断层走滑分量和倾滑分量分别设置为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芦山M_S7.0地震，本文分别构建芦山区域地形模型和无地形模型，探讨地形效应对同震形变的影响。结果表明，地形因素不仅会改变地表位移场分布形态，而且量值在不同方向上也有所差异。对于水平同震形变分量而言，沿断层走向和垂直断层走向的影响最大分别约为10.3%、12.4%；对于垂直同震形变分量而言，最大影响可达约11.9%。反演此次地震断层滑动速率时，若不考虑地形因素的影响，计算结果会低估断层走滑分量和高估倾滑分量。

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

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