Determination of the surface-wave magnitude of the Wenchuan earthquake and its seismic Doppler effect
-
摘要: 本文介绍了全球主要地震机构对2008年5月12日汶川地震参数的速报和修订情况,分析了美国地质调查局国家地震信息中心测定的面波震级。通过对比198个全球地震台站测定的面波震级和面波周期,得出如下结论:测定面波震级偏大的台站主要分布在震中的东北方向,测定面波震级偏小的台站主要分布在震中的西南和东南方向,面波周期偏小的台站主要分布在震中东北方向。由于此次地震破裂方向是以北东向单侧破裂为主,且地震多普勒效应导致震中东北方向振动加强,因此该方向上的面波震级偏大,地震烈度衰减慢;而震中西南方向的振动减弱,此方向面波震级偏小,地震烈度衰减快。从而造成地震烈度沿中央断裂带的北东方向衰减慢,而南西方向衰减快的特征分布。Abstract: The fast report and revision of the parameters of the Wenchuan earthquake on 12 May 2008 by some major international seismological institutions are introduced, and the surface-wave magnitude measured by National Earthquake Information Center of United States Geological Survey is analyzed in this paper. The comparison of the measured surface-wave magnitudes and periods of 198 global seismic stations reveals that, those stations with larger magnitudes are mainly in NE direction of the epicenter while those with smaller magnitudes are in SW and SE direction, and stations with smaller periods are mainly in NE direction. The seismic intensity map of the Wenchuan earthquake shows that the intensity attenuated slowly in the NE direction of the central fracture zone while the intensity attenuated fast in the SW direction. The reason is that Wenchuan earthquake is a unilateral rupture mainly in NE direction, due to the seismic Doppler effect, the stronger shaking leads to the larger surface-wave magnitudes and slower intensity attenuation in the NE direction of the epicenter, vice versa the weaker shaking leads to the smaller surface-wave magnitudes and faster intensity attenuation in the SW direction.
-
Keywords:
- Wenchuan earthquake /
- surface-wave magnitude /
- Doppler effect /
- seismic intensity
-
引言
地震波正演模拟有助于研究地震波传播规律,是研究地球内部结构的基础,同时在资源勘探和地震灾害评估中发挥着重要作用(Tromp et al,2005;Tong et al,2014;Liu et al,2015a;Shen et al,2022;祝贺君等,2023).就地震波运动方程的直接求解而言,研究人员已发展出多种数值模拟方法,例如:有限差分法(Zhang,Yao,2013)、有限元法(finite difference method,缩写为FDM)(Marfurt et al,1984)、谱元法(spectral-element method,缩写为SEM)(Komatitsch,Vilotte,1998;刘少林等,2021)等,每种方法各有其优缺点,其中有限差分法和谱元法是目前关于地震波运动方程的数值求解使用最广泛的方法(桑莹泉等,2021).
有限差分法(FDM)采用泰勒(Taylor)级数近似原理,利用网格点上的波场值直接对方程包含的空间微分项离散(Zhou et al,2021),该方法具有数值格式简单、程序易于实现、计算速度快等优点(Zhang,Yao,2013),因此,FDM在不同尺度的地震波场模拟问题中得到了广泛应用(Moczo et al,2007;Zhang et al,2019).然而,FDM在模拟地震波时存在一些不足,其不足表现为网格划分不灵活、自由边界条件处理困难以及数值频散现象明显等(Liu,2014;Zhang,Yao,2013).前人提出了多种策略克服FDM存在的不足,例如:曲线网格方法(Rao,Wang,2013)、优化系数方法(Zhang,Yao,2013)、自由地表边界条件处理方法(Lan,Zhang,2011).虽然这些改进策略对FDM某方面的不足具有一定的改进效果,但其它方面的不足依然存在。例如:曲线网格的FDM难以以高阶差分格式应用于地震波数值模拟;FDM对自由地表边界条件的近似精度低(苏波等,2019)。
谱元法(SEM)基于变分原理求解弱形式的地震波运动方程,其结合了谱方法的高精度性和有限元法的网格灵活性(Komatitsch,Vilotte,1998),且能自然满足自由地表边界条件(刘少林等,2022)。当采用勒让德(Legendre)型正交插值多项式和高斯-勒让德-洛巴托(Gauss-Legendre-Lobatto,缩写为GLL)数值求积公式时,由于数值积分点与插值点重合的特点(Komatitsch,Tromp,2002;刘少林等,2021),SEM得到的质量矩阵为对角矩阵,而对角矩阵的求逆运算计算量较小,因此SEM用于地震波场模拟时具有较高的计算效率(刘少林等,2021)。在波形伴随层析成像之中,越来越多的研究利用SEM计算地震波形(Chen et al,2015;Lei et al,2020)。
虽然谱元法得到了较为广泛的应用,但该方法存在两方面不足。其一,相比于规则网格的FDM,SEM的计算量明显要大(刘少林等,2021)。在物理坐标向局部坐标转换的过程中,任意方向的空间微分算子要转换为局部坐标所有方向微分算子的加权(Komatitsch,Tromp,2002),该过程增加了计算量;其二,当模型内界面复杂时,网格剖分存在一定困难。谱元法采用六面体单元(三维情形)刻画地下介质结构,当模型内部界面复杂时,难以简单快捷构建高质量网格(Shen et al,2022)。
为了实现复杂含起伏地表模型中的地震波高效数值模拟,前人发展了有限元-有限差分混合方法(Galis et al,2008)。在该混合方法中,利用有限元法模拟起伏地表附近地震波场,采用FDM模拟其他区域中的地震波场。由于有限元法计算量较大,为了提高计算效率,有限元法所计算的区域仅占模型的一小部分(Galis et al,2008)。特别需要指出的是,在该混合方法中FEM难以使用高阶格式得到高精度地震波场。
为了实现复杂模型中地震波场的高精度、高效数值模拟,本文发展了SEM-FDM混合方法。该方法在地表附近采用SEM模拟地震波传播,在远离地表区域利用FDM计算地震波场。混合方法吸收了SEM自由地表边界条件自然满足的优点,具有对复杂地形刻画能力;同时,采用针对任意非均匀介质的有限差分格式,模拟地震波在远离地表区域中的传播,在保证算法计算效率的同时避免了生成复杂网格的难点问题。通过模型试算,并与SEM和FDM等方法进行对比,以证实SEM-FDM混合方法是一种高效的地震波模拟方法。
1. SEM-FDM混合方法原理
以2D非均匀各向同性介质中地震波运动方程为例介绍SEM-FDM混合方法,对于3D或各向异性介质的SEM-FDM混合方法与之类似。首先简要介绍SEM,更详细的推导请见刘少林等(2021)。地震波运动方程可写为:
$$ \begin{array}{l}\left\{\begin{array}{l}\rho \ddot{\boldsymbol{u}}=\nabla ·{\boldsymbol{T}} + {\boldsymbol{f}} \text{,} \\ {\boldsymbol{T}}={\boldsymbol{C}}:\nabla {\boldsymbol{u}} \text{,} \end{array}\right.\end{array} $$ (1) 式中$ \boldsymbol{u}=\left(u\text{,} v\right) $为位移向量,T为应力张量,$ \boldsymbol{f}=\left({f}_{x}\text{,} {f}_{{\textit{z}}}\right) $,$ \nabla =\left(\frac{ {\text{∂}} }{ {\text{∂}} x}\text{,} \frac{ {\text{∂}} }{ {\text{∂}} {\textit{z}}}\right) $,ρ为介质质量密度,$ \boldsymbol{C} $为弹性张量,$ \boldsymbol{u} $上两点表示时间二阶微分。利用任意测试向量w乘以式(1),并利用分部积分和格林公式,式(1)可变为:
$$ \begin{array}{c}{\int }_{\mathrm{\Omega }}\mathrm{\rho }\boldsymbol{w}· \ddot{\boldsymbol{u}}{\mathrm{d}}^{2}{\boldsymbol{x}}=-{\int }_{\mathrm{\Omega }}\nabla \boldsymbol{w}:\boldsymbol{T}{\mathrm{d}}^{2}{\boldsymbol{x}} + {\int }_{ {\text{∂}} \mathrm{\Omega }}\boldsymbol{w}· \boldsymbol{T}· \mathbf{n}\mathrm{d}{\mathrm{s}} + {\int }_{\mathrm{\Omega }}\boldsymbol{w}· \boldsymbol{f}{\mathrm{d}}^{2}{\boldsymbol{x}} \text{,} \end{array} $$ (2) 式中Ω为研究区域,n为边界外法向。根据自由地表边界条件和自然边界条件,式(2)中关于边界的积分项为零。将模型剖分为Ne个非重叠单元,即有$ \mathrm{\mathit{\Omega}}=\bigcup_{\mathrm{i}=1}^{N_{\mathrm{e}}}\mathrm{\mathit{\Omega}}_{\mathrm{i}} $。在每个单元上将物理坐标(x,z)变换成局部坐标(ξ,η),其中ξ和η的取值范围为$ [ $−1,1$ ] $。对于每个局部坐标,根据方程$ ( 1-{\mathrm{\alpha }}^{2} ) {\mathrm{L}}_{\mathrm{N}}' ( \mathrm{\alpha } ) =0 $(其中LN(α)为$ N $阶勒让德多项式)的零点确定插值点。确定插值点以后,每个单元上利用插值多项式近似$ \boldsymbol{u} $和$ \boldsymbol{w} $,其近似公式为:
$$ \begin{array}{c}\left\{\begin{array}{c}{\boldsymbol{u}}^{\mathrm{e}} ( \xi \text{,} \eta ) {\text{≈}} \sum\limits _{i=0}^{N}\sum\limits _{j=0}^{N}{l}_{i} ( \xi ) {l}_{j} ( \eta ) {\boldsymbol{u}}_{ij}^{\mathrm{e}} \text{,} \\ {\boldsymbol{w}}^{\mathrm{e}} ( \xi \text{,} \eta ) {\text{≈}} \sum\limits _{i=0}^{N}\sum\limits _{j=0}^{N}{l}_{i} ( \xi ) {l}_{j} ( \eta ) {\boldsymbol{w}}_{ij}^{\mathrm{e}} \text{,} \end{array}\right.\end{array} $$ (3) 式中l为插值多项式,上标e表示单元编号,下标i和j表示插值节点索引。将式(3)代入式(2),在每个单元上采用GLL数值求积公式,式(2)对应的离散形式为:
$$ \begin{array}{c}\sum\limits_{\mathit{\mathit{\mathrm{\mathit{e}}}}=1}^{N_{\mathrm{e}}}\mathbf{\mathit{\boldsymbol{M}}}^{\mathrm{e}}\ddot{\mathbf{\mathit{\boldsymbol{U}}}}^{\mathrm{e}}=-\sum\limits_{\mathrm{e}=1}^{N_e}\mathbf{\mathit{\boldsymbol{K}}}^{\mathrm{e}}\mathbf{\mathit{\mathbf{U}}}^{\mathrm{e}}+\sum\limits_{\mathrm{e}=1}^{N_e}\mathbf{\boldsymbol{F}}^{\mathrm{e}}\text{,} \end{array} $$ (4) 式中Ue为离散位移向量,Me为单位质量矩阵,Ke为单元刚度矩阵,Fe为单元震源向量。
为了便于介绍FDM,将式(1)写成如下等价形式:
$$ \begin{array}{c}\left\{\begin{array}{c}\rho \dfrac{{ {\text{∂}} }^{2}u}{ {\text{∂}} {t}^{2}}=\dfrac{ {\text{∂}} }{ {\text{∂}} x}\left[ ( \lambda + 2\mu ) \dfrac{ {\text{∂}} u}{ {\text{∂}} x} + \lambda \dfrac{ {\text{∂}} v}{ {\text{∂}} {\textit{z}}}\right] + \dfrac{ {\text{∂}} }{ {\text{∂}} {\textit{z}}}\left[\mu \left(\dfrac{ {\text{∂}} v}{ {\text{∂}} x} + \dfrac{ {\text{∂}} u}{ {\text{∂}} {\textit{z}}}\right)\right] + {f}_{x} \text{,} \\ \rho \dfrac{{ {\text{∂}} }^{2}v}{ {\text{∂}} {t}^{2}}=\dfrac{ {\text{∂}} }{ {\text{∂}} {\textit{z}}}\left[\lambda \dfrac{ {\text{∂}} u}{ {\text{∂}} x} + ( \lambda + 2\mu ) \dfrac{ {\text{∂}} v}{ {\text{∂}} {\textit{z}}}\right] + \dfrac{ {\text{∂}} }{ {\text{∂}} x}\left[\mu \left(\dfrac{ {\text{∂}} v}{ {\text{∂}} x} + \dfrac{ {\text{∂}} u}{ {\text{∂}} {\textit{z}}}\right)\right] + {f}_{{\textit{z}}} \text{,} \end{array}\right.\end{array} $$ (5) 式中λ和μ为各项同性介质的拉梅常数。在式(5)的右端项中,对于只包含关于x或z的偏导数项的FDM离散格式,其离散格式具有相似性。以$ \frac{ {\text{∂}} }{ {\text{∂}} x}\left[\left(\lambda + 2\mu \right)\frac{ {\text{∂}} u}{ {\text{∂}} x}\right] $为例介绍,2N阶有限差分近似格式为:
$$ \begin{split} \frac{ {\text{∂}} }{ {\text{∂}} x}{\left[ ( \lambda + 2\mu ) \frac{ {\text{∂}} u}{ {\text{∂}} x}\right]}_{{i}_{0}\text{,} {j}_{0}} {\text{≈}} & \frac{1}{\Delta x}\sum _{i=0}^{N-1}{b}_{i}\left[{ ( \lambda + 2\mu ) }_{{i}_{0} + i + \frac{1}{2}\text{,} {j}_{0}}{\left(\frac{ {\text{∂}} u}{ {\text{∂}} x}\right)}_{{i}_{0} + i + \frac{1}{2}\text{,} {j}_{0}}-{ ( \lambda + 2\mu ) }_{{i}_{0}-i-\frac{1}{2}\text{,} {j}_{0}}{\left(\frac{ {\text{∂}} u}{ {\text{∂}} x}\right)}_{{i}_{0}-i-\frac{1}{2}\text{,} {j}_{0}}\right] =\\ & \frac{1}{{\Delta x}^{2}}\sum _{i=0}^{N-1}{b}_{i}\left[{ ( \lambda + 2\mu ) }_{{i}_{0} + i + \frac{1}{2}\text{,} {j}_{0}}\sum _{j=0}^{N-1}{b}_{j}\left({u}_{{i}_{0} + i + j + 1\text{,} {j}_{0}}-{u}_{{i}_{0} + i-j\text{,} {j}_{0}}\right) \right. -\\ & \left. { ( \lambda + 2\mu ) }_{{i}_{0}-i-\frac{1}{2}\text{,} {j}_{0}}\sum _{j=0}^{N-1}{b}_{j}\left({u}_{{i}_{0}-i + j\text{,} {j}_{0}}-{u}_{{i}_{0}-i + j-1\text{,} {j}_{0}}\right)\right]\text{\text{,}} \end{split} $$ (6) 式中下标i0和j0分别为x和z方向上的有限差分网格点编号,bi为有限差分系数(Liu,2014),Δx为x方向的网格间距。应该说明的是,式(6)涉及半网格点上的拉梅系数值,可利用整数网格点上的值和调和平均公式(Ma et al,2018)计算。以$ \frac{ {\text{∂}} }{ {\text{∂}} {\textit{z}}}\left[\lambda \frac{ {\text{∂}} v}{ {\text{∂}} x}\right] $为例介绍式(5)中的混合导数有限差分离散格式,$ \frac{ {\text{∂}} }{ {\text{∂}} {\textit{z}}}\left[\lambda \frac{ {\text{∂}} v}{ {\text{∂}} x}\right] $的近似公式如下(Ma et al,2018):
$$ \begin{split} {\left[\frac{ {\text{∂}} }{ {\text{∂}} {\textit{z}}}\left(\lambda \frac{ {\text{∂}} v}{ {\text{∂}} x}\right)\right]}_{{i}_{0} \text{,} {j}_{0}}=& \sum _{j=1}^{N}{\omega }_{j}\frac{\left[{\lambda }_{{i}_{0} \text{,} {j}_{0} + j}{\left(\frac{ {\text{∂}} v}{ {\text{∂}} x}\right)}_{{i}_{0} \text{,} {j}_{0} + j}-{\lambda }_{{i}_{0} \text{,} {j}_{0}-j}{\left(\frac{ {\text{∂}} v}{ {\text{∂}} x}\right)}_{{i}_{0} \text{,} {j}_{0-j}}\right]}{\mathrm{\Delta }{\textit{z}}}=\\& \sum _{j=1}^{N}\frac{{\omega }_{j}}{\mathrm{\Delta }{\textit{z}}}\left[{\lambda }_{{i}_{0} \text{,} {j}_{0} + j}\sum _{{i}_{0}=1}^{N}{\omega }_{{i}_{0}}\frac{\left[{v}_{{i}_{0} + i \text{,} {j}_{0} + j}-{v}_{{i}_{0}-i \text{,} {j}_{0} + j}\right]}{\mathrm{\Delta }x}-{\lambda }_{{i}_{0} \text{,} {j}_{0}-j}\sum _{{i}_{0}=1}^{N}{\omega }_{{i}_{0}}\frac{\left[{v}_{{i}_{0} + i \text{,} {j}_{0}-j}-{v}_{{i}_{0}-i \text{,} {j}_{0}-j}\right]}{\mathrm{\Delta }x}\right]=\\&\frac{1}{\mathrm{\Delta }x\mathrm{\Delta }{\textit{z}}}\sum _{{j}_{0}=1}^{N}\frac{{\omega }_{{j}_{0}}}{2{j}_{0}}\left[{\lambda }_{{i}_{0} \text{,} {j}_{0} + j}\sum _{{i}_{0}=1}^{N}\frac{{\omega }_{{i}_{0}}}{2{i}_{0}}\left({v}_{i + {i}_{0} \text{,} j + {j}_{0}}-{v}_{i-{i}_{0} \text{,} j + {j}_{0}}\right)\right.-\\&\left. {\lambda }_{{i}_{0} \text{,} {j}_{0}-j}\sum _{{i}_{0}=1}^{N}\frac{{\omega }_{{i}_{0}}}{2{i}_{0}}\left({v}_{i + {i}_{0} \text{,} j-{j}_{0}}-{v}_{i-{i}_{0} \text{,} j-{j}_{0}}\right)\right]\text{\text{,}} \end{split} $$ (7) 式中Δz为z方向网格间距,ωi为有限差分系数。应该指出的是,虽然bi和ωi都为有限差分系数,但其值不同。
使用SEM模拟起伏地表附近地震波传播,利用FDM模拟远离地表处的地震波场,即为SEM-FDM混合方法,该方法的优势在于网格建模只需关注起伏地表地形,而无需考虑远离地表处的非规则结构,这种处理方式在保持SEM对起伏地表模型处理能力的同时,极大地简化了SEM网格建模(图1)。SEM-FDM混合方法的核心在于在SEM网格的底边界和FDM的上边界设置数据交换区,设置数据交换区是为了在SEM网格的底边界和FDM网格的上边界满足位移连续性边界条件。将计算区域Ω分解为Ω1和Ω2,即Ω=Ω1+Ω2,Ω1为模型浅部区域,包含起伏地表,Ω2为模型深部区域,地震波运动方程的求解可表示为(假设震源位于Ω1):
$$ \left\{\begin{array}{l}\ddot{\boldsymbol{U}}=\left(\genfrac{}{}{0pt}{}{{\boldsymbol{U}}_{1}}{{\boldsymbol{U}}_{2}}\right)\text{,}\boldsymbol{x}\in {\varOmega } \text{,} \\ {\ddot{\boldsymbol{U}}}_{1}={K}^{\mathrm{S}\mathrm{E}\mathrm{M}}\left({\boldsymbol{U}}_{1} + \boldsymbol{f}\right)\text{,}\boldsymbol{x}\in {{\varOmega }}_{1} \text{,} \\ {\ddot{\boldsymbol{U}}}_{2}={K}^{\mathrm{F}\mathrm{D}\mathrm{M}}\left({\boldsymbol{U}}_{2}\right)\text{,}\boldsymbol{x}\in {{\varOmega }}_{2} \text{,} \end{array}\right. $$ (8) 式中U为离散网格点上的位移向量,其下标1和2分别代表Ω1和Ω2研究区域,KSEM表示SEM的空间离散算子(Liu et al,2017),KFDM表示FDM的空间离散算子(Liu et al,2015a)。当有限差分的阶数为2N时,数据交换区的厚度为$N · \Delta {\textit{z}} $,即Ω1和Ω2重叠区域的厚度为$N · \Delta {\textit{z}} $。此时通过以下三步操作保证SEM网格的底边界和FDM网格的上边界满足位移连续边界条件。第一,在t时刻,由式(7)计算Ω1区域GLL节点上的加速度值,以及Ω2区域有限差分节点上的加速度值,并根据时间中心差分公式(刘少林等,2022),得到t+Δt (其中Δt为时间离散步长)时刻GLL节点和有限差分节点上的波场值。第二,在t+Δt时刻,数据交换区域第一至第N层有限差分网格点上的波场可通过邻近GLL点上的波场值插值近似(图2a),插值近似多项式为:
$$ \left\{\genfrac{}{}{0pt}{}{{{\boldsymbol{U}}}_{\mathrm{F}\mathrm{D}}\left(x,{\textit{z}}\right) {\text{≈}} \sum _{i=0}^{M}\sum _{j=0}^{M}{L}_{ij}{{\boldsymbol{U}}}_{\mathrm{G}\mathrm{L}\mathrm{L}}\left({x}_{i},{{\textit{z}}}_{j}\right),}{{L}_{ij}\left(x,{\textit{z}}\right)=\prod _{k=0 \text{,} k\ne i}^{M}\dfrac{x-{x}_{k}}{{x}_{i}-{x}_{k}}\prod _{k=0 \text{,} k\ne j}^{M}\dfrac{{\textit{z}}-{{\textit{z}}}_{k}}{{{\textit{z}}}_{j}-{{\textit{z}}}_{k}},}\right. $$ (9) 式中下标FD表示有限差分网格节点,下标GLL表示SEM的GLL网格节点,L为拉格朗日插值多项式,M为插值多项式阶数。实际计算表明,插值阶数为3时,可保证数据交换的精度较高;进一步提高插值阶数,对于数据交换精度提升贡献较小。第三,在t+Δt时刻,数据交换区底边界GLL点的位移值可通过临近的有限差分节点上的波场值插值近似(图2b),近似公式如下:
$$ \boldsymbol{U}_{\mathrm{G}\mathrm{L}\mathrm{L}}\left(x,\text{z}\right)\text{≈}\sum_{i=0}^{\mathrm{\mathit{M}}}\sum_{j=0}^{\mathrm{\mathit{M}}}L_{ij}\boldsymbol{U}_{\mathrm{F}\mathrm{D}}\left(x_i,\text{z}_j\right). $$ (10) SEM-FDM混合方法模拟地震波场的流程可由图3描述。相比于SEM和FDM,SEM-FDM混合方法通过设置数据交换区的方式满足研究区内部位移连续性边界条件。应该指出的是,在t+Δt时刻,在近似数据交换区第一至第N层有限差分网格点上的位移时,避免使用数据交换区底边界上GLL点上的位移值用于插值近似。此外,在数据交换区,除了底边界GLL点上位移值需要插值近似得到t+Δt时刻值,其它GLL点上的位移值已经满足连续性边界条件而无需再插值近似。
数值算法精度主要由网格尺度和算法数值阶数来决定。研究表明网格尺寸(或一个波长内的采样点数)是影响数值算法的最主要因素(Liu et al,2015b,2017)。在使用SEM-FDM混合方法模拟地震波场时,为了避免因SEM和FDM精度不匹配而在数据交换区产生虚假散射,在网格设计时,SEM的GLL点平均间距应与FDM的网格间距相近。SEM和FDM的高阶格式能有效较弱数值频散(De Basabe,Sen,2007;Liu et al,2015a),但数值格式阶数过高会导致计算效率下降,SEM的常用阶数为4—5 (Komatitsch,Tromp,2002),FDM的常用差分阶数为4—6 (Tan,Huang,2014)。为了在数值精度与计算效率之间达到平衡,建议一个波长内采样点数为8—10,SEM-FDM混合方法中SEM的阶数选择5,FDM的阶数选择8—10。
2. 稳定性分析
根据Ma等(2014)和Liu等(2015b,2017)的研究,地震波运动方程的稳定性条件可表示为:
$$ \frac{\Delta t\mathrm{\mathit{v}}_{\mathrm{P}}}{h}\sqrt{1+\frac{a}{r^2}}\text{≤}b\text{,} $$ (11) 式中vP为P波速度,γ=vP/vS为P波与S波速度比,h为网格间距(或GLL点的平均间距),a和b为待定常数。De Basabe等(2007)和Liu等(2017)指出,当SEM的空间插值阶数为1时,SEM与标准网格2阶的FDM等价(Liu et al,2017),此时a和b的值都为1;当SEM的空间插值阶数大于1时,a为0,b值列于表1。在对式(5)中只包含x或z方向的导数项近似时,由于利用了半网格点处的空间一阶导数,以至于本文SEM-FDM混合方法中的FDM差分格式与标准网格FDM有一定差别,且稳定性条件与标准网格的FDM之间具有一定差别。根据平面波分析原理(Liu et al,2015b)可得到SEM-FDM混合方法中FDM的稳定性条件,a和b的取值见表2。
表 1 空间不同插值阶数SEM的稳定性参数a和b值Table 1. The stability parameters a and b for the SEMs with different interpolation orders空间插值阶数 a b 1 1 1 2 0 0.66 3 0 0.56 4 0 0.48 5 0 0.41 6 0 0.35 7 0 0.31 8 0 0.28 9 0 0.25 表 2 空间不同差分阶数FDM的稳定性参数a和b值Table 2. The stability parameters a and b for the different order FDMs空间插值阶数 a b 2 1.08 1.00 4 1.07 0.80 6 0.93 0.73 8 0.73 0.69 10 0.70 0.67 12 0.55 0.65 14 0.36 0.63 16 0.31 0.62 18 0.21 0.61 在获得了SEM和FDM各自的稳定性条件以后,可根据SEM和FDM稳定性的下限得到SEM-FDM混合方法的稳定性条件,例如SEM-FDM混合方法中SEM采用5阶插值和FDM采用8阶空间差分,其稳定性条件为:
$$ \begin{split} \\[-8pt] \frac{{\Delta t\mathrm{V}}_{\mathrm{P}}}{h} {\text{≤}} 0.41. \end{split} $$ (12) 3. 数值算例
本节将利用四个数值算例检验SEM-FDM混合方法的计算精度与效率。数值算例中SEM-FDM混合方法的SEM阶数选为5,FDM空间差分阶数选为8。为了减弱边界截断而引起的虚假反射,采用二阶位移形式的完美匹配层(perfectly matched layer,缩写为PML)吸收边界条件(刘有山等,2012)吸收虚假反射。采有较小的时间步长,以减弱时间离散误差对结果的干扰(衡量时间离散大小的标准是式(10)中的库朗数),3.1—3.3节中的时间离散间隔为0.1 ms,3.4节中时间离散间隔为1 ms。
3.1 Lamb问题
均匀介质模型如图4所示,P波速度为3.2 km/s,S波速度为1.8 km/s,密度ρ为2.6×103 kg/m3,一个单位集中力位于(
1000 m,100 m)处,震源时间函数是频率为15 Hz的雷克(Ricker)子波。两个接收器R1和R2分别位于(1500 m,800 m)和(1700 m,0 m)。模型0—1 km深度部分用SEM计算,正方形单元的边长为20 m;模型1—2 km深度部分用FDM计算,网格间隔为5 m。数据交换区位于z=1000 m附近。图 4 用于检验SEM-FDM混合方法计算效率的均匀介质模型模型0—1 km深度部分用SEM计算;1—2 km部分用FDM计算。黑色五角星代表震源;两个黑色三角形对应两个接收器Figure 4. A homogeneous model used to verify the efficiency of the SEM-FDM hybrid methodThe wavefields in upper (0−1 km depths) and lower (1−2 km depths) parts of the model are calculated by the SEM and FDM,respectively. A star represents a source;two triangles denote two receivers R1 and R2上述模型的波场模拟结果如图5所示。为了对比需要,图5也展示了5阶SEM和8阶FDM的波场模拟结果。从t=0.5 s的快照可看出SEM-FDM混合方法、SEM和FDM波场模拟结果大致相同(图5)。波场快照显示除了明显的直达波,还可看到地表的反射波和转换波;由于震源离地表较近,可在地表附近观测到明显的瑞雷面波。从图5a和5b可看到地震波场具有较好的连续性,在z=1 000 m附近的数据交换区未见明显的虚假数值散射现象。
为了说明SEM-FDM混合方法的波场模拟精度,图6展示了合成波形模拟结果,并将数值模拟结果与解析解(De Hoop,1959)对比。可以看出:在模型内部的接收器R1处SEM-FDM混合方法、SEM和FDM的模拟结果与解析解均较接近(图6a和6b),最大相对误差约为3%(图6c和6d);在地表接收器R2处SEM-FDM混合方法和SEM模拟精度较高(图6e和6f),FDM模拟结果存在明显误差(图6e和6f),最大相对误差约15%,其误差明显大于SEM-FDM混合方法和SEM,这是缘于FDM对自由地表边界条件近似精度减低(苏波等,2019)。由于SEM-FDM混合方法能自动满足自由地表边界条件,因此能高精度模拟瑞雷面波(图6e和6f)。
图 6 由解析解得到的合成地震记录及其与三种数值方法得到的地震记录对比(a,b) R1处位移水平和垂直分量;(c,d) R1处合成波形水平和垂直分量的误差;(e,f) R2处位移水平和垂直分量;(g,h) R2处合成波形水平和垂直分量的误差. 黑线、蓝线、红线和青色为解析解、SEM-FDM、SEM和FDMFigure 6. Comparison of synthetic seismograms generated by analytical solution and three numerical methods(a,b) The horizontal and vertical components of displacements recorded at the R1;(c,d) the errors of the horizontal and vertical components of the displacements at the R1;(e,f) the horizontal and vertical components of displacements at the R2,respectively;(g,h) the errors of the horizontal and vertical components of the displacements at the R2表3对比了SEM-FDM混合方法、SEM和FDM三种方法的计算内存和时间。从表3可看出SEM-FDM混合方法的计算内存和时间介于SEM与FDM之间。SEM除了需要存储GLL点上的波场值以外,还需要存储GLL节点上雅克比(Jacobi)矩阵行列式的值(刘少林等,2021,2022),因此SEM的内存开销要大于FDM,但SEM和FDM的内存消耗仍处于同一数量级(表3;刘少林等,2021)。由于SEM-FDM混合方法为SEM和FDM两种方法的综合,因此其内存消耗介于SEM与FDM之间。当使用SEM方法计算空间任意方向的微分时,需要将空间坐标转化为三个参考坐标的加权(刘少林等,2021),而规则网格的FDM方法不存在坐标转化过程(Liu,2014),因此SEM的计算量略大于FDM,但SEM和FDM的计算量仍处于统一数量级。因此对于SEM-FDM混合方法而言,其计算量介于SEM与FDM之间。
表 3 三种地震波正演模拟方法运行内存和计算时间对比Table 3. Computer memory and computation time of three forward modeling methods正演模拟方法 运行内存/kb 计算时间/s SEM-FDM 68157 912.12 SEM 80740 961.31 FDM 28835 821.25 3.2 层状模型
为了说明SEM-FDM混合方法对模型内界面的处理能力,设计如图7所示模型,该模型具有分层结构,且内部含倾斜界面。在SEM-FDM混合方法中,深度500 m之上采用SEM计算,SEM的网格单元平均尺寸为20 m,深度500 m之下利用FDM计算,FDM网格间隔为5 m,数据交换区位于z=500 m附近。模型上层介质P波速度为2.9 km/s,S波速度为1.7 km/s,密度为2.3×103 kg/m3;下层介质P波速度为3.2 km/s,S波速度为1.8 km/s,密度为2.6×103 kg/m3。一个爆炸源位于(1 000 m,100 m),其时间函数为15 Hz的雷克子波。两个接收器R1和R2分别位于(500 m,800 m)和(500 m,
1400 m)。SEM-FDM混合方法地震波场模拟结果如图8所示。为了对比,也给出了5阶SEM波场模拟结果。需要指出的是,SEM网格考虑了模型内部界面,即网格剖分时沿着倾斜界面剖分。图8显示SEM-FDM混合方法与SEM方法所得的波场快照模拟结果具有很好的一致性,从波场快照可看到穿过倾斜界面的透射P波和pP波,以及倾斜界面的反射Pp波、Ps波和pPp波。接收器R1和R2处的合成地震记录如图9所示,可以看出SEM-FDM混合方法模拟结果与参考波形(6阶SEM计算得到)高度一致,相对误差在1%左右,这说明SEM-FDM混合方法能较好处理含速度间断面的模型。
图 8 层状模型中t=0.5 s时刻的波场水平分量(左)和垂直分量(右)快照(a,b) SEM-FDM计算所得;(c,d) SEM计算所得.Figure 8. Snapshots of displacement wavefields of horizontal component (left) and horizontal component (right) in a layered model at t=0.5 s generated by the SEM-FDM (a,b) and SEM (c,d)The left plots show the displacement wavefields of horizontal component;left plots display the displacement wavefields of vertical component. The upper and lower plots are generated by the SEM-FDM and SEM,respectively图 9 混合方法得到的合成地震记录与参考解对比(a,b) R1处位移水平和垂直分量;(c,d) R1处合成波形水平和垂直分量的误差;(e,f) R2处位移水平和垂直分量;(g,h) R2处合成波形水平和垂直分量的误差Figure 9. Comparison of synthetic seismograms generated by the SEM-FDM hybrid method and reference solution(a,b) The horizontal and vertical components of displacements recorded at the R1,respectively;(c,d) Errors of the horizontal and vertical components of the displacements at the R1;(e,f) the horizontal and vertical components of displacements at the R2;(g,h) the errors of the horizontal and vertical components of the displacements at the R23.3 起伏地表模型
为了检验SEM-FDM混合方法对起伏地表模型的处理能力,设计了如图10所示的起伏地表模型。模型分为两层,每层介质参数与3.2节中模型相同。一个爆炸震源位于(1 000 m,100 m)处,震源时间函数是频率为15 Hz的雷克子波。200个接收器均匀放置于地表,其平均水平间隔为16 m。
图 10 含起伏地表的层状模型中地震波场模拟地表曲线代表起伏地形; 模型内部实线表示倾斜界面. 五角星表示震源,其位置为(1 000 m,100 m). 三角形表示接收器,其平均水平间隔为16 m,为了显示需要,接收器每隔四个显示一个Figure 10. Simulation of seismic wavefields in a layered model with topographyThe topography is marked by the black line. The star located at (1 000 m,100 m) represents a source. The triangles denote receivers,the average horizontal distance of which is about 16 m. For visualization purposes,only one-fifth stations are displayed单炮地震记录模拟结果如图11所示。图11a和11d由SEM-FDM混合方法计算得到;为了对比需要,图11b和11e展示了5阶SEM计算结果,可见SEM-FDM混合方法和SEM计算结果具有很好的一致性;图11c和11f展示了两种方法模拟结果的相对差,相对差在0.8%以内。通过对比可见,SEM-FDM混合方法保持了SEM对起伏地表的处理能力,能高精度模拟起伏地表模型中的地震波传播。
图 11 SEM-FDM混合方法计算的单炮地震记录与SEM计算结果对比图(a)和(b)分别为SEM-FDM和SEM方法计算的位移水平分量,图(c)为其相对差;图(d)和(e)分别为SEM-FDM和SEM方法计算的位移垂直分量,图(f)为其相对差Figure 11. A comparison of single shot profiles generated by the SEM-FDM hybrid method and SEMPlots (a) and (b) show displacement wavefields of horizontal component by the SEM-FDM and SEM,respectively,and plot (c) shows the relative differences between the wavefields shown in plots (a) and (b). Plots (d) and (e) display displacement wavefields of vertical component by the SEM-FDM and SEM,respectively,and plot (f) displays the relative differences between the wavefields shown in plots (d) and (e)3.4 龙门山断裂带地壳模型
为了检验SEM-FDM混合方法在实际地质模型中的地震波模拟能力,选取龙门山断裂带地壳模型(图12)。图12显示龙门山断裂带东西两侧的地形相差较大,东侧的高程较低,西侧的高程较高。模型P波(图12b)和S波(图12c)速度结构来源于川滇地区速度模型2.0 (Liu等,2023)。选取二维垂直剖面(图12a中红色线表示其位置)用于正演合成测试。一个爆炸震源位于(110 km,13 km),震源时间函数是主频为1 Hz的雷克子波。两个接收器R1和R2分别位于(77.4 km,−3.4 km)和(136.4 km,−0.6 km)(其中z坐标的负值表示海平面之上)。SEM-FDM混合方法在深度z=5 km之上利用SEM计算,网格单元的平均尺寸为800 m,在z=5 km之下采用FDM计算,网格间隔为320 m,数据交换区位于z=5 km附近。
图 12 龙门山断裂带地壳模型F1:龙门山断裂带(LMSF);F2:岷江断裂(MJF);F3:虎牙断裂(HYF);F4:龙泉山断裂(LQSF);F5:安宁河断裂(ANHF);SGB: 松潘-甘孜块体; SCB: 四川盆地. (a) 地质构造图. 红色线为AB垂直剖面所在位置,用于二维数值测试;(b) P波速度结构;(c) S波速度结构Figure 12. The crustal model of the Longmenshan fault zoneF1:Longmenshan fault zone (LMSF);F2:Minjiang fault (MJF);F3:Huya fault (HYF);F4:Longquanshan fault (LQSF);F5:Anninghe fault (ANHF). (a) The topography around the Longmenshan fault zone. The black lines denote faults. The red line AB represent a vertical profile,that is used for a 2D numerical test.(b) P-wave velocity model;(c) S-wave velocity modelSEM-FDM混合方法的地震波数值模拟结果如图13和14所示。为了对比需要,图13和14也展示了5阶SEM数值模拟结果。可见:SEM-FDM混合方法得到的t=13 s时刻波场快照显示了清晰的P,pP和pS震相,波场模拟结果(图13a和13c)与SEM模拟结果一致(图13b和图13d)。R1和R2处的合成波形模拟结果(图14)表明SEM-FDM混合方法与SEM得到的合成波形具有很好的一致性。以上模拟结果说明SEM-FDM能模拟复杂地壳模型中的地震波场。
图 13 地壳模型中t=13 s时刻的波场快照. 左图为位移水平分量;右图为位移垂直分量(a)和(b)为SEM-FDM混合方法计算得到;(c)和(d)为SEM计算得到Figure 13. Snapshots of wavefields in the crustal model at t=13 s. The left and right plots show horizontal and vertical components of displacement wavefieldsPlots (a) and (b) are generated by the SEM-FDM hybrid method;plots (c) and (d) are computed by the SEM图 14 地壳模型中合成地震记录左图为位移水平分量;右图为位移垂直分量。图(a)和(b)为R1处合成地震记录;图(c)和(d)为R2处合成地震记录. 黑线为参考地震图; 蓝线为SEM-FDM混合方法计算得到Figure 14. Synthetic seismograms in a crustal modelLeft and right plots show synthetic seismograms of horizontal and vertical components,respectively;Plots (a) and (b) present seismograms recorded by the R1;plots (a) and (b) display seismograms recorded by the R2. The black lines are the reference seismograms. The blue lines are the synthetic seismograms computed by the SEM-FDM hybrid method4. 讨论与结论
本文提出了模拟地震波场的SEM-FDM混合方法,该方法在接近地表区域利用SEM模拟地震波传播,在远离地表区域采用FDM计算地震波场,通过设置数据交换区将SEM和FDM耦合,使得模型内部的位移连续。发展SEM-FDM混合方法是为了简化网格建模,即在网格建模时只用考虑起伏地表附近模型的非规则性,在远离地表区域利用规则有限差分网格刻画模型内部界面。SEM-FDM混合方法综合了SEM和FDM的各自优势,可简单、便捷地模拟复杂模型中的地震波场。通过四组数值算例证实了SEM-FDM混合方法模拟地震波场的有效性。本文设计的地下介质模型相对简单,主要目的是为了简化参考、对比对象SEM的网格建模,便于SEM-FDM混合方法计算结果与SEM对比。SEM-FDM混合方法可适用于任意模型,并不受模型的限制。
在利用SEM-FDM混合方法模拟地震波场过程中,利用商业软件Hypermesh (https://www.altair.com)建立地表附近SEM网格模型,为了进一步提高网格建模效率,可构造高阶多项式拟合地表地形,在此基础之上自动生成SEM网格,作者将另文介绍SEM网格自动生成技术。基于网格自动生成技术,SEM-FDM混合方法将成为一种高效、便捷的正演模拟工具,可模拟任意复杂介质中的地震波场。
在波形伴随成像过程中,随着模型迭代更新,反演得到的速度模型与初始模型之间可能存在较大差别。目前,基于SEM的波形伴随成像在模型迭代过程中难以根据反演模型重新剖分网格。由于SEM网格未能有效地刻画反演后的速度结构异常体,这必然造成合成波形精度较低,进而影响成像结果。本文的SEM-FDM混合方法利用较密的有限差分网格刻画速度结构的高波数成分,因此无须在反演迭代过程中再重新划分网格。作者将另文介绍基于SEM-FDM混合方法的波形伴随成像方法与软件。
-
表 1 国际地震机构速报的汶川地震的参数
Table 1 Parameters of the Wenchuan earthquake quickly reported by major international seismological institutions
序号 发震时刻 (UTC)
时:分:秒震中位置 震源
深度/kmMS 测定机构 北纬/° 东经/° 中文名称 代码 1 06:28:04.1 30.95 103.40 14 8.0 中国地震台网中心 CENC 2 06:28:00.9 31.10 103.30 10 7.8 美国国家地震信息中心 NEIC 3 06:27:59.3 31.10 103.30 10 8.0 俄罗斯科学院 RAS 4 06:27:58.9 31.10 103.20 10 7.5 欧洲地中海地震中心 EMSC 5 06:28:00.8 30.80 103.40 10 6.8 (mb) 罗马尼亚地球物理研究所 RNIEP 注:表中序号1数据来自中国地震台网中心 (2008a),其它数据来自Swiss Seismological Service (2008)。 表 2 国际地震机构测定的汶川大地震的参数
Table 2 Parameters of the Wenchuan earthquake observed by major international seismological institutions
序号 发震时刻 (UTC)
时:分:秒震中位置 震源深
度/km震级类型 测定机构 数据来源 北纬/° 东经/° mb MS MW 名称 代码 1 06:27:59.5 31.01 103.42 14 6.4 8.2 中国地震台网中心 CENC 中国地震台网中心 (2008a) 2 06:28:01.8 31.00 103.32 19 6.9 8.1 7.9 美国国家地震信息中心 NEIC NEIC (2008b) 3 06:28:41.4 31.49 104.11 12 7.8 7.9 全球矩心矩张量项目
数据中心GCMT GCMT (2008) 4 06:27:59.0 31.10 103.20 10 7.9 欧洲地中海地震中心 EMSC EMSC (2008) 5 06:28:03.7 31.60 103.70 33 7.2 8.4 德国格拉芬堡地震
观测中心SZGRF SZGRF (2008) 注:全球矩心矩张量项目数据中心测定的是“矩心”(即所释放的地震矩的“时-空几何中心”)的位置,其物理意义与传统的震源位置及发震时刻(地震初始破裂的位置与时刻)不相同,不具有简单的可比性。 表 3 《地震数据报告》中列出的汶川地震的209个台站的参数和MSZ (NEIC,2008a)
Table 3 The related data of 209 stations and MSZ of Wenchuan earthquake in Earthquake Data Report(NEIC,2008a)
序号 台站代码 震中距/° 方位角/° 周期/s MSZ 序号 台站代码 震中距/° 方位角/° 周期/s MSZ 1 INCN 20.30 65.1 19.0 8.4 37 KONO 65.59 326.2 20.0 8.0 2 HIA 22.06 29.5 20.0 7.7 38 CTAO 65.42 135.3 21.0 7.7 3 MDJ 24.68 49.2 20.0 7.9 39 TIR 65.73 304.9 22.0 7.9 4 TKM2 25.03 306.1 22.0 7.8 40 AKUT 66.02 40.1 18.0 8.1 5 KULM 25.70 186.1 20.0 7.9 41 FOO 66.98 329.1 18.9 8.3 6 AAK 25.71 304.9 22.0 7.9 42 BER 67.32 327.8 18.0 8.3 7 AML 26.14 303.4 19.0 8.0 43 TRI 68.27 311.0 19.0 7.9 8 EKS2 26.22 304.6 22.0 7.8 44 MIDW 67.56 69.7 19.0 7.9 9 KURK 26.93 324.1 19.0 8.1 45 GRF 68.33 315.7 20.0 8.4 10 KKM 27.62 151.4 20.0 7.6 46 GRA1 68.33 315.7 20.0 8.4 11 MAJO 29.41 69.7 19.0 8.0 47 TIP 68.75 303.6 22.0 7.9 12 KSM 30.10 166.0 20.0 7.7 48 SDPT 68.43 37.6 20.0 8.1 13 BRVK 32.53 322.3 20.0 7.9 49 RER 69.10 227.5 22.0 7.7 14 ERM 33.58 59.7 22.0 7.8 50 COLA 69.81 25.5 20.0 8.1 15 YSS 34.14 50.8 19.0 8.3 51 CEL 69.73 303.0 21.0 7.7 16 YAK 35.44 21.2 19.0 8.0 52 AQU 69.82 307.8 20.0 8.0 17 ARU 40.10 322.6 22.0 7.9 53 OHAK 71.22 34.0 21.0 8.3 18 GUMO 41.88 105.0 20.0 7.6 54 VLC 71.09 310.6 22.0 8.1 19 COCO 43.39 189.3 19.0 7.9 55 KDAK 71.29 33.3 19.0 8.3 20 PET 45.29 44.4 21.0 8.3 56 ECH 71.38 315.4 22.0 7.9 21 GNI 47.76 297.9 20.0 7.7 57 WLF 71.34 317.1 21.0 8.0 22 DGAR 48.30 223.0 19.0 7.6 58 WDD 71.72 301.2 20.0 7.6 23 KIV 48.81 303.2 19.0 8.2 59 EGAK 72.27 24.0 19.0 8.3 24 BILL 51.55 25.3 19.0 8.2 60 TARA 71.87 99.6 21.0 7.6 25 MBWA 54.20 161.0 22.0 7.7 61 BNI 73.08 312.6 21.0 8.3 26 SMY 54.64 44.8 19.0 8.2 62 ABPO 73.31 235.1 19.0 7.8 27 BR13 56.21 299.5 21.0 7.6 63 MID 73.27 30.0 18.0 8.3 28 KEV 56.55 336.1 20.0 8.2 64 ESK 73.69 325.2 20.0 8.2 29 MSEY 57.80 240.8 21.0 7.8 65 LOR 73.79 315.6 18.5 8.4 30 WAKE 57.84 85.6 20.0 7.8 66 VSL 73.80 306.5 22.0 7.7 31 ISP 59.07 298.1 22.0 7.9 67 SSB 74.31 313.6 20.0 8.1 32 KBS 60.08 347.1 20.0 8.1 68 FLN 75.60 318.4 20.0 8.3 33 ADK 60.35 44.6 22.0 8.1 69 BORG 75.31 338.6 19.0 8.0 34 PSZ 63.55 311.7 21.0 7.8 70 RJF 76.16 314.6 19.0 8.2 35 NAO0 64.60 327.2 20.0 7.9 71 SKAG 77.73 26.5 18.0 8.3 36 NWAO 64.93 167.1 22.0 8.1 72 CAN 78.72 143.4 20.0 7.9 73 UCH 25.58 304.0 21.0 7.8 109 SNZO 97.50 133.9 20.0 7.8 74 SIT 79.33 28.3 21.0 8.2 110 LAO 98.19 20.1 20.0 8.2 75 SFJD 80.32 349.9 21.0 8.3 111 HOPS 98.07 34.7 21.0 8.1 76 CRAG 81.31 28.7 18.0 8.2 112 LKWY 98.75 23.7 22.0 8.4 77 PAB 82.91 312.1 19.0 8.1 113 RLMT 98.64 22.7 21.0 8.2 78 ESLA 82.59 312.0 19.0 8.8 114 AGMN 99.21 12.8 19.0 8.2 79 FUNA 82.43 104.6 20.0 7.6 115 BMN 99.61 30.2 22.0 8.1 80 TAU 83.91 149.1 22.0 7.8 116 EYMN 100.29 10.1 20.0 8.3 81 MTE 84.41 314.2 20.0 8.0 117 CMB 100.12 33.7 22.0 8.0 82 PAF 85.26 201.0 21.0 7.7 118 AHID 100.07 25.0 19.0 8.2 83 LSZ 85.50 249.3 22.0 7.8 119 ELK 100.19 28.8 21.0 8.2 84 SFS 85.88 310.4 20.0 8.5 120 BW06 100.62 24.0 20.0 8.2 85 KIP 86.24 67.4 22.0 7.6 121 HWUT 100.91 25.9 21.0 8.2 86 RTC 87.76 308.7 20.0 7.8 122 RSSD 101.17 19.7 19.0 7.6 87 POHA 89.09 67.5 20.0 7.8 123 DUG 101.64 27.5 19.0 7.9 88 NLWA 91.17 29.9 21.0 8.2 124 COWI 102.43 8.8 21.0 8.4 89 FFC 91.88 14.4 21.0 8.0 125 DBIC 101.94 285.3 20.0 8.0 90 LBTB 92.86 242.6 20.0 7.8 126 LIC 102.33 285.0 21.0 7.7 91 NEW 93.17 25.7 21.0 8.2 127 MAW 102.89 194.8 20.0 8.1 92 HAWA 93.77 28.1 20.0 8.1 128 ECSD 103.42 14.7 19.0 8.3 93 SCHQ 94.10 354.3 20.0 8.2 129 OGNE 104.67 19.6 20.0 8.2 94 MSO 95.54 24.7 22.0 8.0 130 ISCO 104.61 22.7 19.0 8.2 95 EGMT 95.92 21.6 21.0 8.1 131 LONY 104.72 358.5 20.0 8.1 96 HUMO 95.30 32.4 22.0 8.1 132 JFWS 105.40 10.3 19.0 8.2 97 BMO 95.93 27.8 20.0 8.0 133 NCB 105.36 358.2 20.0 8.3 98 YBH 96.10 32.8 19.0 7.9 134 SCIA 105.87 12.8 21.0 8.3 99 TSUM 96.16 251.4 21.0 7.9 135 MVCO 105.91 26.0 19.0 8.1 100 CASY 97.13 177.1 21.0 8.0 136 WUAZ 106.23 28.9 20.0 8.0 101 ULM 97.33 12.3 20.0 8.3 137 SDCO 106.50 23.5 21.0 8.1 102 DGMT 97.14 18.1 22.0 8.3 138 BINY 107.16 359.5 20.0 8.2 103 WDC 97.08 33.4 22.0 8.0 139 KSU1 107.83 16.2 20.0 8.3 104 BOZ 97.38 23.9 20.0 8.1 140 HDIL 107.85 10.1 20.0 8.2 105 MOD 97.07 31.3 19.0 8.1 141 ANMO 108.66 25.5 22.0 8.0 106 XMAS 96.86 83.5 19.0 7.7 142 ACSO 108.89 5.1 20.0 8.2 107 WVOR 97.36 29.9 20.0 8.0 143 CCM 109.97 12.2 19.0 7.8 108 HLID 98.15 26.7 20.0 8.1 144 TUC 109.23 30.2 21.0 8.1 145 WVT 112.41 9.8 19.0 8.1 178 LVZ 53.66 334.0 21.0 8.5 146 WCI 110.54 8.1 19.0 7.4 179 JOHN 79.19 77.0 21.0 7.8 147 AMTX 110.36 21.8 21.0 8.1 180 MCCM 98.82 35.2 20.0 8.0 148 MNTX 111.96 26.2 20.0 8.0 181 SAO 100.61 35.1 20.0 7.9 149 WMOK 111.38 19.4 20.0 8.4 182 TPH 101.55 31.6 20.0 7.9 150 SBA 114.81 168.0 19.0 8.0 183 DAC 102.85 32.9 20.0 8.0 151 MIAR 112.89 15.1 21.0 8.2 184 MVU 103.33 27.9 22.0 8.1 152 NATX 115.29 16.9 19.0 8.3 185 GLMI 104.19 5.8 20.0 8.2 153 GOGA 115.61 6.3 22.0 8.3 186 LBNH 105.00 356.5 19.0 8.1 154 VBMS 115.77 13.0 19.0 8.5 187 AAM 106.78 5.4 19.0 8.2 155 LRAL 115.57 9.6 19.0 8.3 188 CBKS 107.30 18.7 21.0 8.2 156 HKT 116.71 18.6 20.0 7.6 189 MCWV 109.64 2.6 20.0 8.4 157 BRAL 117.54 10.0 19.0 8.2 190 BLA 112.05 3.2 20.0 8.4 158 KVTX 118.36 21.4 19.0 8.3 191 CNNC 114.10 1.1 22.0 8.2 159 DWPF 121.04 4.9 20.0 7.9 192 SHEL 114.20 265.4 20.0 7.6 160 GRTK 127.52 353.5 20.0 8.0 193 BBSR 115.96 348.7 22.0 8.1 161 TGUH 134.06 14.4 19.0 8.0 194 NHSC 116.15 3.3 20.0 8.3 162 ANWB 129.54 341.5 20.0 8.0 195 MOL 65.53 329.6 23.1 8.2 163 SDDR 130.03 353.3 21.0 8.2 196 KMBO 70.23 256.2 20.0 7.5 164 SJG 130.13 346.9 20.0 8.4 197 COR 93.6 31.4 21.0 7.3 165 GTBY 129.34 358.1 22.0 8.2 198 PFO 105.57 12.8 21.0 7.5 166 MTDJ 131.05 1.1 20.0 8.1 199 BJO1 58.90 341.8 17.7 8.5 167 FDF 132.14 339.6 21.0 8.1 200 TRO 59.37 336.1 13.8 8.5 168 JTS 138.18 12.3 21.0 8.1 201 NSS 63.04 331.0 17.2 8.3 169 GRGR 134.75 339.1 22.0 8.1 202 PMR 70.95 28.8 17.0 8.2 170 RCBR 134.32 294.7 19.0 8.0 203 CLL 66.67 316.9 22.1 8.1 171 BCIP 139.95 4.9 22.0 7.9 204 CPUP 162.4 280.5 21.0 7.2 172 PMSA 145.30 189.5 21.0 7.9 205 NNA 161.08 0.5 21.0 7.8 173 PAYG 147.19 25.7 20.0 7.9 206 LPAZ 163.4 330.0 22.0 8.1 174 OTAV 148.89 3.4 20.0 7.9 207 TRQA 166.0 235.9 21.0 8.2 175 RPN 151.25 90.6 21.0 7.8 208 PLCA 169.1 205.4 19.0 7.9 176 SPB 153.04 278.7 20.0 8.2 209 LCO 174.5 289.4 20.0 8.2 177 EFI 155.20 208.3 21.0 8.0 平均 8.1 注:表中平均值计算未包括增加的11个台站资料。 -
陈运泰. 2008. 汶川地震的成因断层、破裂过程与成灾机理[C]//中国科学院第十四次院士大会学部学术报告汇编. 北京: 中国科学院: 38–39. Chen Y T. 2008. Generating faults, rupture process and disaster mechanism of the Wenchuan earthquake[C]//Collection of the Academician Reports on the 14th General Assembly, Chinese Academy of Sciences. Beijing: Chinese Academy of Sciences: 38–39 (in Chinese).
陈运泰, 许力生, 张勇, 杜海林, 冯万鹏, 刘超, 李春来, 张红霞. 2009. 汶川特大地震震源特性分析报告[G]//汶川大地震工程震害调查分析与研究. 北京: 中国岩石力学与工程学会: 6–17. Chen Y T, Xu L S, Zhang Y, Du H L Feng W P, Liu C, Li C L , Zhang H X. 2009. Analysis report of the source characters of the great Wenchuan earthquake[G]//Investigation, Analysis and Research on the Engineering Earthquake Damages of the Great Wenchuan Earthquake. Beijing: Chinese Society for Rock Mechanics and Engineering: 6–17 (in Chinese).
郭履灿, 庞明虎. 1981. 面波震级和它的台基校正值[J]. 地震学报, 3 (3): 312-320. Guo L C, Pang M H. 1981. Surface-wave magnitude and its station correction [J]. Acta Seismologica Sinica, 3 (3): 312-320 (in Chinese).
胡聿贤. 2006. 地震工程学[M]. 北京: 地震出版社: 54–56. Hu Y X. 2006. Earthquake Engineering [M]. Beijing: Seismological Press: 54–56 (in Chinese).
刘瑞丰, 陈运泰, 任枭, 徐志国, 王晓欣, 邹立晔, 张立文. 2015. 震级的测定[M]. 北京: 地震出版社: 1–10. Liu R F, Chen Y T, Ren X, Xu Z G, Wang X X, Zou L Y, Zhang L W. 2015. Determination of Magnitude [M]. Beijing: Seismological Press: 1–10 (in Chinese).
刘瑞丰, 陈运泰, Peter Bormann, 任枭, 侯建民, 邹立晔. 2006. 中国地震台网与美国地震台网测定震级的对比(Ⅱ): 面波震级[J]. 地震学报, 28(1): 1-7. Liu R F, Chen Y T, Bormann P, Ren X, Hou J M, Zou L Y. 2006. Comparison between earthquake magnitudes determined by China seismograph network and U.S. seismograph network (Ⅱ): Surface wave magnitude [J]. Acta Seismologica Sinica, 28 (1): 1-7(in Chinese).
许绍燮. 1999. 中华人民共和国国家标准《地震震级的规定》(GB 17740—1999)宣贯教材[M]. 北京: 中国标准出版社: 1–10. Xu S X. 1999. Teaching Book for Propagandizing and Implementing the National Standard of the People’s Republic of China: General Ruler for Earthquake Magnitude (GB 17740–1999) [M]. Beijing: Standards Press of China: 1–10 (in Chinese).
袁一凡. 2008. 四川汶川8.0级地震灾害损失评估[J]. 地震工程与工程振动, 28(5): 10-19. Yuan Y F. 2008. The earthquake disaster loss evaluation[J]. Earthquake Engineering and Engineering Vibration, 28(5): 10-19(in Chinese).
张勇, 许力生, 陈运泰. 2009. 2008年汶川大地震震源机制的时空变化[J]. 地球物理学报, 52 (2): 379-389. Zhang Y, Xu L S, Chen Y T. 2009. Spatio-temporal variation of the source mechanism of the 2008 great Wenchuan earthquake [J]. Chinese Journal of Geophysicis, 52 (2): 379-389(in Chinese).
中国地震台网中心. 2008a. 地震信息[DB/OL]. [2008–05–12]. http://www.cenc.ac.cn. China Earthquake Networks Center. 2008a. Earthquake information [DB/OL]. [2008–05–12]. http://www.cenc.ac.cn (in Chinese).
中国地震台网中心. 2008b. 中国地震台站观测报告[R]. 北京: 中国地震台网中心: 1–10. China Earthquake Networks Center. 2008b. Seismological Report of the Chinese Seismic Stations[R]. Beijing: China Earthquake Networks Center: 1–10 (in Chinese).
中国地震台网中心. 2008c. 中国数字地震台网观测报告[R]. 北京: 中国地震台网中心: 1–10. China Earthquake Networks Center. 2008c. Seismological Report of the Chinese Digital Seismic Network[R]. Beijing: China Earthquake Networks Center: 1–10 (in Chinese).
周庆, 张春山, 陈献程. 2011. 汶川MS8.0地震灾害的非对称分布与成因分析[J]. 地震学报, 33(4): 492-504. Zhou Q, Zhang C S , Chen X C. 2011. Asymmetric disaster distribution and its cause analysis of the MS8.0 Wenchuan earthquake. Acta Seismologica Sinica, 33 (4): 492-504(in Chinese).
EMSC. 2008. M7.9 − Eastern Sichuan, China − 2008-05-12 06:28:00 UTC[EB/OL]. [2017–10–22]. https://www.emsc-csem.org/Earthquake/earthquake.php?id=85969.
GCMT. 2008. Monthly CMT solutions[EB/OL]. [2017–10–22]. http://www.ldeo.columbia.edu/~gcmt/projects/CMT/catalog/NEW_MONTHLY/2008/mar08.ndk.
Gutenberg B. 1945.Amplitudes of surface waves and magnitudes of shallow earthquakes[J].Bull. Seism. Soc. Am. 35: 3-12.
NEIC. 2008a. The preliminary determination of epicenters (PDE) bulletin[EB/OL]. [2014–10–01]. ftp://hazards.cr.usgs.gov/NEICPDE/olderPDEdata/mchedr/.
NEIC. 2008b. M7.9−eastern Sichuan, China [EB/OL]. [2017–10–22]. https://earthquake.usgs.gov/earthquakes/eventpage/usp000g650
Somerville P G, Smith N F, Graves R W, Abrahamson N A. 1997.Modification of empirical strong ground motion attenuation relations to include the amplitude and duration effects of rupture direction[J]. Seis Res Lett, 68(1): 199-222.
Swiss Seismological Service. 2008. Recent earthquakes Switzerland[DB/OL]. [2008–05–12]. http://www.seismo.ethz.ch.
SZGRF. 2008. Seismological center observatory[DB/OL]. [2008–07–01]. https://www.szgrf.bgr.de/bulletins.html.