局部场地地震动高频衰减系数估计模型

稂子平, 俞瑞芳, 肖亮, 傅磊, 周健

稂子平,俞瑞芳,肖亮,傅磊,周健. 2023. 局部场地地震动高频衰减系数估计模型. 地震学报,45(5):919−928. DOI: 10.11939/jass.20220053
引用本文: 稂子平,俞瑞芳,肖亮,傅磊,周健. 2023. 局部场地地震动高频衰减系数估计模型. 地震学报,45(5):919−928. DOI: 10.11939/jass.20220053
Lang Z P,Yu R F,Xiao L,Fu L,Zhou J. 2023. An estimation model of high frequency attenuation coefficient of ground motion for local site. Acta Seismologica Sinica45(5):919−928. DOI: 10.11939/jass.20220053
Citation: Lang Z P,Yu R F,Xiao L,Fu L,Zhou J. 2023. An estimation model of high frequency attenuation coefficient of ground motion for local site. Acta Seismologica Sinica45(5):919−928. DOI: 10.11939/jass.20220053

局部场地地震动高频衰减系数估计模型

基金项目: 北京市自然科学基金(8212018)、中央级公益性科研院所基本科研业务费专项(DQJB21B36)和国家重点研发计划课题(2017YFC0404901)联合资助
详细信息
    作者简介:

    稂子平,在读硕士研究生,主要从事地震动模拟方面的研究,e-mail:williamlangzp@163.com

    通讯作者:

    俞瑞芳,博士,研究员,主要从事结构抗震理论、地震动特性分析及模拟等方面的研究,e-mail:yrfang126@126.com

  • 中图分类号: P315.9

An estimation model of high frequency attenuation coefficient of ground motion for local site

  • 摘要:

    采用随机有限断层法进行地震动模拟时,选用合理的参数描述特定局部场地近地表高频衰减特征,对评价地震动模拟结果的正确与否具有重要的实践意义。 在工程场址地震动参数预测中,如何快速确定参数的取值,是实际应用中亟需解决的问题。首先对场地高频衰减系数κ0与平均剪切波速vS30的相关性进行了分析;然后,基于国内外学者计算得到的546个κ0系数,采用一定时窗内的κ0均方根值,讨论其随平均剪切波速vS30增加的变化趋势。 结果表明,虽然κ0具有明显的区域差异性,但其均方根值随着vS30的增大呈现出逐渐减小的趋势。 为了得到合理的κ0估计模型,分别采用线性函数、多项式函数、对数线性函数和双对数线性函数对κ0均方根值与vS30的关系进行初步拟合,结果表明,对数线性函数能够较好地描述κ0vS30之间的关系。 最后,基于筛选得到的477个数据,采用最小二乘法对模型参数进行拟合,建立了适合工程应用的κ0-vS30模型。对模型适用性的分析表明,本研究所构建的κ0估计模型能够合理估计地震动的高频衰减影响。

    Abstract:

    When using the stochastic finite fault method for ground motion simulation, how to select reasonable parameters to describe the near-surface high-frequency attenuation characteristics of a specific local site has important practical significance for evaluating the correctness of ground motion simulation results. In the prediction of ground motion parameters of engineering sites, how to quickly determine the value of parameters is an urgent problem to be solved in practical applications. Firstly, we analyzed the correlation between the high-frequency attenuation coefficient κ0 of the site and the average shear wave velocity vS30. Then, based on the 546 κ0 coefficients calculated by domestic and foreign scholars, the root mean square value of κ0 in a certain time window was used to discuss its variation trend with the increase of the average shear wave velocity vS30.The results showed that although κ0 had obvious regional differences, its root mean square value showed a decreasing trend with the increase of vS30. In order to obtain a reasonable κ0 estimation model, the linear function, polynomial function, logarithmic linear function and log-log linear function were used to preliminarily fit the relationship between the root mean square value of κ0 and vS30. The results show that the logarithmic linear function can better describe the relationship between κ0 and vS30. Finally, based on the 477 data obtained from the screening, the model parameters were fitted by the least square method, and a practical model of κ0-vS30 suitable for engineering applications was obtained. The analysis of the applicability of the model shows that the κ0 estimation model constructed in this study can reasonably estimate the high-frequency attenuation of ground motion when predicting engineering site ground motion parameters.

  • 覆盖土层对地震波传播特性具有重要影响(高武平等,2012李平等,2012),通过土层地震反应分析获得工程场地的地面运动特征,为工程结构抗震评估、设计和加固提供依据。许多实际工程场地都可以被简化为水平成层场地模型,即假定覆盖土层的力学性质沿竖向呈水平成层变化,沿横向均匀无限延伸,将地震激励考虑为垂直向上入射的剪切波。按照此模型分析,土层地震反应分析实际即为一维波动问题。由于该模型形式简洁,物理意义明确,至今依然是场地地震反应分析的一种重要途径(Zalachoris,Rathje,2015)。

    常用的一维土层地震反应分析方法主要包括频域和时域两类方法。频域法以等效线性化方法为代表,最早由Idriss和Seed (1968)提出。该方法将不同应变水平下的剪切模量和阻尼比以一个等效剪切模量和阻尼比代替,从而把非线性问题转化为线性问题并在频域内求解。廖振鹏(1989)根据等效线性化方法开发出了土层地震反应分析计算程序LSSRLI-1,在实践中得到了广泛认可。近年来,袁晓铭等(2016)采用直频法动剪模量阻尼比求解技术,提出了新一代土层地震反应分析方法,克服了传统方法低估软弱场地和深厚场地放大效应的缺陷。孙锐和袁晓铭(2021)提出了全局等效剪应变的概念和算法,建立了一种新的等效线性化分析方法。然而,等效线性化方法尚存在一些缺点,比如无法反映土层真实的受力状态(王志良,韩清宇,1981栾茂田,林皋,1992),且忽略了较多的高频成分(齐文浩,薄景山,2007)。为了解决上述问题,发展出了土层地震反应分析的时域方法,可以真实地反映土层在地震作用下的表现。丁海平和周正华(1998)提出了一维土层地震反应分析的时域有限元解法,并指出时域算法能够达到与频域算法相同的计算精度,但其计算耗时远少于频域算法。然而传统有限单元法由于形函数阶次较低,通常需将每个波长内布置6—10个单元(廖振鹏,2002)才能准确刻画出波场特征,在土层复杂或地震动高频成分较丰富的情况下可能会需要较大的计算工作量。对此,邢浩洁等(2017)提出了分析成层场地地震反应的切比雪夫(Chebyshev)谱元法,仅需划分少量单元即可获得较高的计算精度。但其在计算切比雪夫谱单元质量矩阵时沿用了传统谱元法的做法,即利用切比雪夫多项式的性质获得单元质量矩阵的解析解,由此导出的质量矩阵为非对角形式,无法充分发挥显式时间积分算法的高效率优势。

    本文在邢浩洁等(2017)的工作基础上,拟通过节点积分法(Fried,Malkus,1975)建立集中质量切比雪夫谱单元,解决传统切比雪夫谱元法由于需要对质量矩阵求逆而造成计算效率不高的问题,同时避免传统的质量集中方法的随意性,如行和集中法(Zienkiewicz et al,2013)或对角元素放大法(Hinton et al,1976)。在该集中质量切比雪夫谱单元模型中嵌入多次透射人工边界(multi-transmitting formula,缩写为MTF)(Liao,Wong,1984),结合中心差分形成一种求解一维土层波动问题的高阶显式算法.并利用日本Kik-net强震观测台网提供的井下和地面实测记录,以检验本文方法的适用性。

    将水平成层场地简化为如图1所示的计算模型。设坐标原点位于自由地表,建立正方向竖直向下的一维坐标系z。假设有N−1个土层覆盖于基岩半空间上,各覆盖土层的厚度、质量密度和剪切模量分别为hiρiGii = 1,2,…,N-1。在基岩半空间中,从基岩面向下截取一定深度hN,并在底部设置多次透射人工边界模拟无限域对内域波动的影响,以消除外行波场在边界处的反射。从而,N-1个土层和1个基岩层组成的计算模型总厚度为$H {\text{=}} \sum\nolimits_{i {\text{=}} 1}^N {{h_i}}$。入射地震波为从模型底部竖直向上传播的剪切波。采用切比雪夫谱元法对土层模型进行空间离散,每个土层划分为若干个大小相等的一维切比雪夫谱单元,以图1中第i土层为例,将hi厚度的土层划分为两个大小为∆zi的谱单元,而基岩层则被单独划分为一个谱单元。

    图  1  水平成层场地切比雪夫谱元模型
    Figure  1.  Chebyshev spectral element model of horizontal layered soil site

    一维土层分析的切比雪夫谱单元以 [ −1,1 ] 区间内不等间距分布的高斯-洛巴托-切比雪夫(Gauss-Lobatto-Chebyshev,缩写为GLC)节点为参考单元节点,它们是第一类切比雪夫多项式的极值点。对于n阶切比雪夫谱单元,GLC节点的位置为ξi=−cos(iπ/n),i=0,1,···,n

    设单元位移模式为

    $$ u( \xi ) {\text{=}} \sum\limits_{i {\text{=}} 0}^n {{\varphi _i}( \xi )u( {{\xi _i}} )} {\text{,}} $$ (1)

    式中:uξi)为节点水平位移;$\varphi _i $ξ)为建立在单元节点上的拉格朗日形函数;考虑到GLC节点与切比雪夫多项式之间的关系,该形函数亦可表示为一组截断的切比雪夫多项式组合

    $$ {\varphi _i}( \xi ) {\text{=}} \frac{2}{n}\sum\limits_{k {\text{=}} 0}^n {\frac{1}{{{c_i}{c_k}}}{T_k}( {{\xi _i}} ){T_k}( \xi )} {\text{,}}\;\;\;\;i {\text{=}} 0{\text{,}}1{\text{,}}\cdots {\text{,}}n {\text{,}} $$ (2)

    式中,Tkξ)为k阶第一类切比雪夫多项式,其表达式为Tkξ)=cos(kcos−1ξ),cick为多项式系数,当i=0或n时,取值为2,当i=1,···,n-1时,取值为1。

    利用切比雪夫谱单元对水平成层场地模型进行空间离散后,得到离散方程

    $$ {\boldsymbol{M}}\ddot {\boldsymbol{u}} {\text{+}} {\boldsymbol{C}}\dot {\boldsymbol{u}} {\text{+}} {\boldsymbol{Ku}} {\text{=}}{\boldsymbol{0}} {\text{,}} $$ (3)

    式中:$ {\boldsymbol{u}} $$ \dot {\boldsymbol{u}}$$\ddot {\boldsymbol{u}} $分别为位移、速度和加速度;总体质量矩阵M;总体阻尼矩阵C和总体刚度矩阵K分别由单元刚度矩阵Me、单元阻尼矩阵Ce和单元刚度矩阵Ke按照单元编号组装而成。MeCeKe计算如下:

    $$ {M}_{ij}^{\rm e} {\text{=}} \int_{ - 1}^1 {\rho {\varphi _i}( \xi ){\varphi _j}( \xi )\left| {{\boldsymbol{J}}} \right|{\rm d}\xi } {\text{,}} $$ (4)
    $$ {C}_{ij}^{\rm e} {\text{=}} \int_{ - 1}^1 {\eta {{\boldsymbol{J}}^{ - 1}}{{\varphi '}_i}( \xi ){{\boldsymbol{J}}^{ - 1}}{{\varphi '}_j}( \xi )\left| {{\boldsymbol{J}}} \right|{\rm d}\xi } {\text{,}} $$ (5)
    $$ {K}_{ij}^{\rm e} {\text{=}} \int_{ - 1}^1 {G{{\boldsymbol{J}}^{ - 1}}{{\varphi '}_i}( \xi ){{\boldsymbol{J}}^{ - 1}}{{\varphi '}_j}( \xi )\left| {\boldsymbol{J}} \right|{\rm d}\xi } {\text{,}} $$ (6)

    式中ρηG分别表示土体的质量密度、阻尼系数和剪切模量,$ \left| {\boldsymbol{J}} \right| $为雅可比矩阵的行列式。对于一维谱单元,雅可比矩阵仅包含一个元素,即J=∆z/2,∆z为物理单元的长度。

    为了建立集中质量矩阵,本文采用节点积分法求解单元质量矩阵,即以单元节点作为数值积分点,利用高斯-洛巴托(Gauss-Lobatto)数值积分计算单元质量矩阵的积分(式(4))。由于拉格朗日形函数具有克罗内克(Kronecker-δ)性质,由此导出的质量矩阵仅有主对角元素非零,其余非对角元素全部为零。单元GLC节点对应的高斯-洛巴托积分权系数为

    $$ {w_i} = \int_{ - 1}^1 {{\varphi _i}( \xi ){\rm{d}}\xi }{\text{,}} $$ (7)

    式中:wi为与节点ξi相对应的积分权系数;φiξ)为式(2)中的单元形函数。1—5阶切比雪夫谱单元的GLC节点位置及对应的积分权系数,列于表1

    表  1  GLC节点上的高斯−洛巴托积分权系数
    Table  1.  Gauss-Lobatto quadrature weights based on GLC points
    谱单元阶次GLC节点坐标积分权系数
    1±11
    2±1, 00.333 3, 1.333 3
    3±1, ±0.50.111 1, 0.888 9
    4±1, ±0.707 1, 00.066 7, 0.533 3, 0.8
    5±1, ±0.809 0, ±0.309 00.04, 0.360 7, 0.599 3
    下载: 导出CSV 
    | 显示表格

    利用式(7)的数值积分格式计算式(4)单元质量矩阵,则质量矩阵的非对角元素全部为零,主对角元素可进一步计算为

    $$ M_{ii}^{\rm e} {\text{=}} \rho \left| {\boldsymbol{J}} \right|\sum\limits_{k {\text{=}} 0}^n {{w_k}{\varphi _i}( {{\xi _k}} ){\varphi _i}( {{\xi _k}} )} {\text{=}} \rho \left| {\boldsymbol{J}} \right|{w_i} {\text{.}} $$ (8)

    由于切比雪夫正交多项式积分性质表现为

    $$ {B}_{kl}{\text{=}}{\displaystyle {\int }_{-1}^{1}{{T}}_{k}^{\prime }(\xi ){{T}}_{l}^{\prime }(\xi ){\rm d}\xi }{\text{=}}\left\{\begin{array}{l}0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;k{\text{+}}l{\text{为奇数}}{\text{,}} \\ \dfrac{kl}{2}\left[{H}_{\left|{}\frac{{k - l}}{2}\right|}{\text{-}}{H}_{\left|\frac{{k + l}}{2}\right|}\right]\;\;\;\;\;\;\;\;k+l{\text{为偶数}}{\text{,}} \end{array} \right. $$ (9)

    其中:

    $$ {H_i} {\text{=}} \left\{ {\begin{array}{*{20}{l}} {0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i {\text{=}} 0}{\text{,}} \\ { - 4\displaystyle\sum\limits_{j {\text{=}} 1}^i {\frac{1}{{2j {\text{-}} 1}}} {\text{,}}\;\;\;\;\;\;\;\;i {\text{≥}} 1} {\text{,}} \end{array}} \right. $$ (10)

    利用以上性质,CeKe可精确计算为

    $$ C_{ij}^{\rm e} {\text{=}} \eta {\left|{\boldsymbol{ J}} \right|^{ - 1}}\frac{4}{{{n^2}{c_i}{c_j}}}\sum\limits_{k {\text{=}} 0}^n {\sum\limits_{l {\text{=}} 0}^n {\frac{1}{{{c_k}{c_l}}}{T_k}( {{\xi _i}} )} } {T_l}( {{\xi _j}} ){B_{kl}} {\text{,}}$$ (11)
    $$ K_{ij}^{\rm e} {\text{=}} G{\left| {\boldsymbol{J}} \right|^{ - 1}}\frac{4}{{{n^2}{c_i}{c_j}}}\sum\limits_{k {\text{=}} 0}^n {\sum\limits_{l {\text{=}} 0}^n {\frac{1}{{{c_k}{c_l}}}{T_k}( {{\xi _i}} )} } {T_l}( {{\xi _j}} ){B_{kl}} {\text{.}} $$ (12)

    对于土层地震反应分析问题,一般将入射波作为边界处的运动形式施加,因此需要将模型内域节点和边界节点分开处理。在本文建立的水平成层场地谱元模型中,边界节点仅有一个,即模型底部人工边界上的单元节点,其余节点均属于内域节点。将式(3)中的矩阵和向量按照内域和边界进行分块,得

    $$ \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{M}}_{\rm{I}}}}&{\boldsymbol{0}} \\ {\boldsymbol{0}}&{{{\boldsymbol{M}}_{\rm{B}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\ddot {\boldsymbol{u}}}_{\rm{I}}}} \\ {{{\ddot {\boldsymbol{u}}}_{\rm{B}}}} \end{array}} \right] {\text{+}} \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{C}}_{\rm{I}}}}&{{{\boldsymbol{C}}_{{\rm{IB}}}}} \\ {{{\boldsymbol{C}}_{{\rm{BI}}}}}&{{{\boldsymbol{C}}_{\rm{B}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\dot {\boldsymbol{u}}}_{\rm{I}}}} \\ {{{\dot {\boldsymbol{u}}}_{\rm{B}}}} \end{array}} \right] {\text{+}} \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{K}}_{\rm{I}}}}&{{{\boldsymbol{K}}_{{\rm{IB}}}}} \\ {{{\boldsymbol{K}}_{{\rm{BI}}}}}&{{{\boldsymbol{K}}_{\rm{B}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{u}}_{\rm{I}}}} \\ {{{\boldsymbol{u}}_{\rm{B}}}} \end{array}} \right] {\text{=}} {\boldsymbol{0}} {\text{,}} $$ (13)

    式中:${\ddot {\boldsymbol{u}}_{\rm{I}}}$${\dot {\boldsymbol{u}}_{\rm{I}}}$${{\boldsymbol{u}}_{\rm{I}}}$分别为内域节点的加速度、速度和位移向量;${\ddot {\boldsymbol{u}}_{\rm{B}}}$${\dot {\boldsymbol{u}}_{\rm{B}}}$${{\boldsymbol{u}}_{\rm{B}}}$分别为边界节点的加速度、速度和位移;MIMB为总体质量矩阵M按照内域节点和边界节点的节点编号进行分块后的子矩阵;CICIBCBICB以及KIKIBKBIKB分别为总体阻尼矩阵C和总体刚度矩阵K按照相同的方式分块后的子矩阵。将上式中的第一行方程单独写为

    $$ {{\boldsymbol{M}}_{\rm{I}}}{\ddot {\boldsymbol{u}}_{\rm{I}}} {\text{+}} {{\boldsymbol{C}}_{\rm{I}}}{\dot {\boldsymbol{u}}_{\rm{I}}} {\text{+}} {{\boldsymbol{C}}_{{\rm{IB}}}}{\dot {\boldsymbol{u}}_{\rm{B}}} {\text{+}} {{\boldsymbol{K}}_{\rm{I}}}{{\boldsymbol{u}}_{\rm{I}}} {\text{+}} {{\boldsymbol{K}}_{{\rm{IB}}}}{{\boldsymbol{u}}_{\rm{B}}} {\text{=}} {\bf 0} {\text{,}} $$ (14)

    对式(14)进行时间域离散,设时间间隔为∆t,并采用中心差分法进行时步积分。在tp∆t时刻,内域节点的加速度和速度可插值表示为

    $$ \ddot {\boldsymbol{u}}_{\rm{I}}^p {\text{=}} \frac{{{\boldsymbol{u}}_{\rm{I}}^{p {\text{+}} 1} {\text{-}} 2{\boldsymbol{u}}_{\rm{I}}^p {\text{+}} {\boldsymbol{u}}_{\rm{I}}^{p {\text{-}} 1}}}{{\Delta {t^2}}}{\text{,}}\dot {\boldsymbol{u}}_{\rm{I}}^p {\text{=}} \frac{{{\boldsymbol{u}}_{\rm{I}}^{p {\text{+}} 1} {\text{-}} {\boldsymbol{u}}_{\rm{I}}^{p {\text{-}} 1}}}{{2\Delta t}} {\text{,}} $$ (15)

    式中,上标p表示p∆t时刻;边界节点的速度插值表示为

    $$ \dot {\boldsymbol{u}}_{\rm B}^p {\text{=}} \frac{{{\boldsymbol{u}}_{\rm B}^{p {\text{+}} 1} {\text{-}} {\boldsymbol{u}}_{\rm B}^{p {\text{-}} 1}}}{{2\Delta t}} {\text{.}} $$ (16)

    将式(15)和(16)代入式(14),得到内域节点位移的递推格式

    $$ {\boldsymbol{u}}_{ {\rm{I}}}^{p {\text{+}} 1} {\text{=}} {\left( {\frac{1}{{\Delta {t^2}}}{{\boldsymbol{M}}_{\rm{I}}} {\text{+}} \frac{1}{{2\Delta t}}{{\boldsymbol{C}}_{\rm I}}} \right)^{ - 1}}\left[ {\left( {\frac{2}{{\Delta {t^2}}}{{\boldsymbol{M}}_{\rm{I}}} {\text{-}} {{\boldsymbol{K}}_{\rm I}}} \right){\boldsymbol{u}}_{\rm{I}}^p {\text{-}} } {\left( {\frac{1}{{\Delta {t^2}}}{{\boldsymbol{M}}_{\rm I}} {\text{-}} \frac{1}{{2\Delta t}}{{\boldsymbol{C}}_{\rm I}}} \right){\boldsymbol{u}}_{\rm I}^{p {\text{-}} 1} {\text{-}} {{\boldsymbol{C}}_{{\rm{IB}}}}\dot {\boldsymbol{u}}_{\rm{B}}^p {\text{-}} {{\boldsymbol{K}}_{{\rm{IB}}}}{\boldsymbol{u}}_{\rm{B}}^p} \right] {\text{.}} $$ (17)

    一般情况下质量矩阵MI为对角阵,而阻尼矩阵CI通常不会呈对角矩阵形式,故由MICI组合的矩阵亦会是非对角阵。此情形下可考虑构建仅同质量矩阵呈比例的瑞雷阻尼矩阵,或者将阻尼矩阵对角化(Thomson et al,1974)等措施。

    由针对节点的总波场给出的方程式(17)可知,求p+1时刻的内域节点位移时,除已知的p时刻和p-1时刻内域节点位移外,还需要知道p时刻边界点的位移${ { {\boldsymbol{u}}}_{\rm{B}}^p} $。而当采用MTF人工边界时,还需要先将总波场分离为内行波和外行波,再针对外行波应用MTF条件。

    图2所示,人工边界取基岩面以下一定深度,将这一深度范围内的基岩划分为一个谱单元。底部边界节点的总位移场${\boldsymbol u}_{\rm b}^p$可分解为

    图  2  切比雪夫谱单元的一阶MTF插值方案
    Figure  2.  Interpolation scheme for first-order MTF in Chebyshev spectral element
    $$ {\boldsymbol u}_{\rm b}^p {\text{=}} {\boldsymbol u}_{{\rm bi}}^p {\text{+}} {\boldsymbol u}_{{\rm bs}}^p {\text{,}} $$ (18)

    式中,${\boldsymbol{u}}_{{\rm bi}}^p$为已知的p时刻边界上的入射波,${\boldsymbol{u}}_{{\rm bs}}^p$p时刻边界上的外行波。根据一阶MTF条件,边界点的外行波由距离边界caΔt处的外行波递推得到,即

    $$ {\boldsymbol u}_{{\rm {bs}}}^p {\text{=}} {\boldsymbol u}_{{\rm {ms}}}^{p {\text{-}} 1} {\text{,}} $$ (19)

    式中,${\boldsymbol{u}}_{{\rm {ms}}}^{p{{\text{-}}}1}$表示p-1时刻距离边界caΔt处的MTF插值点的外行波场,ca为人工波速,在本文中ca取基岩的剪切波速。由于MTF插值点与单元节点一般不重合,故${{u}}_{{\rm {ms}}}^{p{\text{-}}1}$需要根据p-1时刻整个基岩单元的外行波场插值求出,即

    $$ u_{{\rm {ms}}}^{p {\text{-}} 1}{\text{=}}{a_0}u_{{\rm {bs}}}^{p {\text{-}} 1} {\text{+}} {a_1}u_{{\rm {s1}}}^{p {\text{-}} 1} {\text{+}} \cdots {\text{+}} {a_n}u_{{ {{\rm{s}}n}}}^{p {\text{-}} 1}{\text{,}} $$ (20)

    式中$u_{\rm{s}1}^{p{\text{-}}1}{\text{,}} \cdots {\text{,}}u_{{\rm{s}}n}^{p{\text{-}}1}$为基岩单元中非边界节点在p-1时刻的外行波。根据式(1)所定义的单元位移模式,插值系数a0a1,···,an即为定义在各单元节点上的形函数在MTF插值点处的取值。基岩单元中所有节点在p-1时刻的外行波场则可从已求出的p-1时刻总波场中扣除内行波场得到

    $$ \left[ {\begin{array}{*{20}{c}} {u_{\rm {bs}}^{p-1}} \\ {u_{\rm {s1}}^{p-1}} \\ \vdots \\ {u_{{\rm {s}}n}^{p-1}} \end{array}} \right] {\text{=}} \left[ {\begin{array}{*{20}{c}} {u_{\rm b}^{p-1}} \\ {u_1^{p-1}} \\ \vdots \\ {u_{ n}^{p-1}} \end{array}} \right] {\text{-}} \left[ {\begin{array}{*{20}{c}} {u_{\rm {bi}}^{p-1}} \\ {u_{\rm {i1}}^{p-1}} \\ \vdots \\ {u_{{\rm {i}}n}^{p-1}} \end{array}} \right] {\text{,}} $$ (21)

    式中,$u_{{\rm {bi}}}^{p{\text{-}}1}$p-1时刻边界上的入射波,$u_{\rm {i1}}^{p{{\text{-}}}1}{\text{,}}\cdots {\text{,}}u_{{\rm {i}}{n}}^{p{{\text{-}}}1}$$u_{{\rm {bi}}}^{p{{\text{-}}}1}$ 滞后一定时间间隔,该时间间隔的长度等于各点到边界点的距离与基岩剪切波速之比。至此,求解土层模型所有节点的位移反应得以显式地完成。

    对于土层地震波动的有限元模拟,廖振鹏(2002)建议第i层土层中的有限单元尺寸取为$\Delta {{\textit{z}}_i} {\text{≤}} {{{c_i}{T_{\min }}} \mathord{\left/ {\vphantom {{{c_i}{T_{\min }}} }} \right. } }$αα取6—10),其中ci为第i层土层的剪切波速,Tmin为具有工程意义的输入波的最短周期。而采用本文提出的集中质量切比雪夫谱元法时,经过试算发现一般情况下谱单元阶次n取4即可取得良好的精度,此时谱单元尺寸取$\Delta {{\textit{z}}_i} {\text{≤}} {c_i}{T_{\min }}$即可,单元尺寸明显大于传统有限单元法。

    本文采用的中心差分法为有条件稳定算法,时间步距∆t需要满足稳定性要求

    $$ \Delta t {\text{≤}} {\alpha _{{\rm{CFL}}}}\frac{{{{\rm{d}}{\rm{z}}_{\min }}}}{{{c_{\max }}}}{\text{,}} $$ (22)

    式中:dzmin为模型网格中任意两个节点之间的最小距离,cmax为最大剪切波速;αCFL为无量纲的库朗-弗雷德里奇-李维(Courant-Friedrichs-Lewy,缩写为CFL)条件数,本文将其取为0.75可保证数值结果始终稳定。

    为了检验本文方法的正确性,首先对一个均匀半空间场地模型进行分析,根据行波理论得到其解析解。该模型土体厚度为180 m,土体质量密度为2000 kg/m3,剪切波速为250 m/s。底部设置一阶MTF人工边界,垂直向上输入幅值为1 m、主频为2 Hz的雷克(Ricker)波位移脉冲。将整个模型划分为8个4阶谱单元,则单元尺寸约为输入波最短波长的一半。图3显示了土层底部和地表的位移反应时程,其中底部反应时程的第一个波峰为入射波,第二个波峰为地表反射波。由图可知,入射波从底部传至地表时间约为0.7 s,与解析解吻合,同时,地表位移反应峰值等于入射波峰值的2倍,符合自由地表条件。此外,地表未出现模型底部反射而来的地震波,说明MTF人工边界条件成功实现。本例表明集中质量切比雪夫谱元法能够处理土层地震反应分析问题,而且在网格较为稀疏的情况下依然能够获得较高的计算精度。

    图  3  均匀半空间模型在雷克波入射下的地表和底部位移反应
    Figure  3.  Displacement responses of ground surface and bottom of a homogeneous half-space model under Ricker wave incidence

    利用日本Kik-net强震观测台网(NIED,2021)提供的实测地震动记录和钻孔数据检验本文方法处理实际场地地震反应分析的能力。从Kik-net台网中随机选取4个具有代表性的台站,分别对应中国《GB 50011—2010 建筑抗震设计规范》(中华人民共和国住房和城乡建设部,中华人民共和国国家质量监督检验检疫总局,2010)所定义的Ⅰ1,Ⅱ ,Ⅲ和Ⅳ等四类场地。在每个台站中选择实测地面峰值加速度(peak ground acceleration,缩写为PGA)等于0.05—0.1g,0.1—0.2g和0.2—0.4g的EW向或NS向地震记录各一条,分别对应较弱地震动、中等强度地震动和较强地震动。以井下实测记录作为一维水平成层场地模型底部的基岩输入波,计算地表加速度反应并与地表实测加速度记录进行对比。由于Kik-net数据库中未提供土体质量密度数据,本文利用Boore (2016)提出的公式根据P波波速和剪切波速估算土体密度。各台站的土层剖面剪切波速如图4所示,按照《GB 50011—2010建筑抗震设计规范》计算得到的基本信息列于表2

    图  4  四个Kik-net台站的土层剪切波速图
    Figure  4.  Shear wave velocity profiles at the four Kik-net stations
    表  2  选用Kik-net台站的基本信息
    Table  2.  Basic information of selected Kik-net stations
    台站名称覆盖层厚度/m等效剪切波速/(m·s−1场地类别
    TCGH08 2 170 1
    MYGH10 34 330
    KMMH14 88 200
    IKRH02 >108 112
    下载: 导出CSV 
    | 显示表格

    4个台站在不同强度地震作用下地表的加速度时程反应,如图5所示。可见,在多数情况下,本文方法计算得到的地表加速度时程与实测的地表加速度记录均能较好吻合。与输入的基岩实测加速度时程相比,所有台站的地表计算反应均显示出了土层的放大效应。即便对于IKRH02台站这样覆盖层厚度超过100 m的深厚场地,也依然体现出了土层的放大效应,未出现类似土层时域非线性软件DEEPSOIL严重低估地表反应的问题(袁晓铭等,2016)。

    图  5  不同强度地震动作用下四个台站地表加速度反应时程
    (a) 较弱地震动;(b) 中等强度地震动;(c) 较强地震动
    Figure  5.  Ground acceleration response histories of four stations under ground motions with different intensities
    (a) Weak ground motions;(b) Moderate ground motions; (c) Strong ground motions

    表3给出了按本文方法计算得到的PGA与地震中实际记录PGA的对比。由表中数据可知,对于Ⅱ类场地(MYGH10),本文预测PGA最接近实际值,误差小于2.1%。但对于Ⅲ类场地(KMMH14),本文方法计算得到的PGA大于实测值。总体而言,本文方法计算得到的PGA总是稍小于或偏大于实际值,而未出现严重偏小的情况。表4列出了本文计算得到的PGA放大倍数及实际测得的PGA放大倍数。所选台站的实测放大倍数介于3.28—7.40之间不等,而计算放大倍数在3.67—10.48之间不等。对于Ⅰ 1类场地(TCGH08)、 Ⅱ 类场地(MYGH10)和Ⅳ类场地(IKRH02)在小震和中震情况下本文方法得到的PGA放大倍数均与实测值相差不大,但在大震情况下本文结果偏大。

    表  3  不同强度地震动作用下各台站计算PGA与实测PGA对比
    Table  3.  Comparison of computed PGA and recorded PGA for the stations under ground motions different intensities
    台站名称较弱地震动中等强度地震动较强地震动
    计算PGA/g实测PGA/g计算PGA/g实测PGA/g计算PGA/g实测PGA/g
    TCGH08 0.090 0.065 0.142 0.145 0.226 0.204
    MYGH10 0.066 0.067 0.154 0.156 0.240 0.235
    KMMH14 0.121 0.083 0.227 0.144 0.369 0.228
    IKRH02 0.082 0.084 0.102 0.104 0.329 0.227
    下载: 导出CSV 
    | 显示表格
    表  4  不同强度地震动作用下各台站PGA放大倍数
    Table  4.  Amplification factors of PGA for the stations under ground motions different intensities
    台站名称较弱地震动中等强度地震动较强地震动
    井下实测
    峰值/g
    计算放大
    倍数
    实测放大
    倍数
    井下实测
    峰值/g
    计算放大
    倍数
    实测放大
    倍数
    井下实测
    峰值/g
    计算放大
    倍数
    实测放大
    倍数
    TCGH08 0.015 5.96 4.31 0.037 3.80 3.89 0.028 8.21 7.40
    MYGH10 0.016 4.14 4.22 0.042 3.67 3.71 0.065 3.69 3.62
    KMMH14 0.016 7.36 5.07 0.031 7.36 4.69 0.035 10.48 6.48
    IKRH02 0.018 4.54 4.63 0.023 4.45 4.55 0.069 4.75 3.28
    下载: 导出CSV 
    | 显示表格

    图6为4个台站在不同强度的地震作用下的地表加速度反应谱,阻尼比按5%计算.总体而言,本文计算反应谱在低频段(周期>0.5 s)与实测记录的反应谱较为接近。观察图6可知,在较弱的地震作用下,Ⅰ1和 Ⅱ 类场地的反应谱峰值与实测记录相近,而Ⅲ和Ⅳ类场地的反应谱形状与实测值吻合,其中Ⅳ类场地的计算反应谱曲线与实测反应谱十分接近。在中等强度地震作用下,Ⅰ1和Ⅳ类场地的反应谱峰值与实测记录较为吻合。在较强地震作用下,本文计算反应谱曲线的形状大致与实测记录得到的反应谱一致,但对于Ⅲ和Ⅳ类场地计算得到的反应谱峰值比观测结果偏大.

    图  6  不同强度地震动作用下四个台站的地表加速度反应谱
    (a) 较弱地震动;(b) 中强地震动;(c) 较强地震动
    Figure  6.  Ground acceleration response spectra of four stations under ground motions with different intensities
    (a) Weak ground motions;(b) Moderate ground motion; (c) Strong ground motions

    造成本文计算结果与实际台阵记录存在一定差异的原因可能有:① 本文采用的一维水平成层场地模型是一种高度简化的计算模型,无法反映出实际三维场地的地形特征,因此一维波动的数值计算结果一般无法体现出对实际场地地震动有重要影响的地形放大效应;② 实际场地并不一定仅受剪切波作用,还有可能受到纵波、面波等多种地震波的作用,并且地震波的传播方向也不一定恰好垂直向上的;③ 本文方法未考虑土体的非线性特性。

    本文提出了一种求解水平成层场地地震反应的时域集中质量切比雪夫谱元法。通过节点积分法严格地导出集中质量矩阵,克服了传统切比雪夫谱元法由于质量矩阵为非对角矩阵形式所带来的计算效率不高的问题。采用中心差分法进行时间积分,并嵌入多次透射人工边界,形成了高效的时域显式算法。与传统有限单元法相比,该方法在每个波长内仅需布置一个谱单元即可获得令人满意的计算精度,显著降低了对空间网格尺寸的要求.

    选择日本Kik-net强震台网中分属四种不同场地类型的台站记录检验了本文方法处理实际场地的能力。计算结果表明,本文方法对于Ⅰ 1,Ⅱ和Ⅳ类场地在较弱地震和中等强度地震作用下能够较好地预测地面运动特征,但对于Ⅲ类场地及强震作用下的地表反应计算结果与观测结果相比仍存在一定误差,后续工作中可考虑添加土体的时域非线性本构关系.

  • 图  7   整体数据校验和MTN台站数据校验情况

    Figure  7.   Overall data check and MTN station data check

    图  1   本文研究采用的κ0值随平均剪切波速 vS30的分布

    Figure  1.   The distribution of κ0 value used in this study with mean shear wave velocity vS30

    图  2   κ0值统计数据箱型图

    Figure  2.   κ0 value statistics box plot

    图  3   四种函数对κrms的拟合情况

    Figure  3.   Fitting of four functions to κrms

    图  4   数据量分布情况图

    Figure  4.   Data distribution diagram

    图  5   对数线性函数拟合以及与其它模型对比

    Figure  5.   Log-linear function fitting and comparison with other models

    图  6   对数线性函数拟合残差图

    Figure  6.   Log-linear Function fit residual plot

    表  1   本文研究采用的κ0数目、来源及相应的vS分布范围

    Table  1   The number and source of κ0 used in this study and the corresponding distribution range of vS

    序号数据个数 vS30/(m·s−1地区数据来源
    1 60 106.8—904.2 日本 Cabas等(2017
    2 16 213.2—744.1 日本 Cabas等(2017
    3 50 507.7—1 433.4 日本 van Houtte等(2011
    4 27 1 106.8— 2 394.0 日本 van Houtte等(2011
    5 4 515.7—1 301.3 日本 Laurendeau等(2013
    6 14 170.6—1 428.1 法国 Drouet等(2010
    7 24 192.1—747.1 瑞士 Edwards等(2015
    8 8 1 174.0—1 810.5 瑞士 Edwards等(2011
    9 16 380.1—1 811.5 瑞士 Edwards等(2011
    10 54 160.1—942.8 中国台湾 Huang等(2017
    11 4 233.1—684.8 中国台湾 Lai等(2016
    12 10 167.5—496.4 中国台湾 Lai等(2016
    13 29 191.8—746.9 土耳其 Bora等(2017
    14 16 142.6—1 029.6 意大利 Bora等(2017
    15 5 1 054.7—1 392.5 克罗地亚 Stanko等(2017
    16 4 854.1—953.3 克罗地亚 Stanko等(2017
    17 11 516.6—715.5 中国台湾 van Houtte等(2011
    18 38 299.6—652.9 中国 傅磊和李小军(2017
    19 10 435.7—1 518.3 新西兰 van Houtte等(2018
    20 9 401.0—661.8 亚利桑那 Kishida等(2014
    21 6 550.9—1 000.3 加利福尼亚 van Houtte等(2011
    22 117 170.9—1 531.2 中国台湾 Chang等(2019
    23 7 263.3—660.3 中国 郑旭等(2019
    24 7 531.2—912.1 日本 朱百慧(2016
    下载: 导出CSV

    表  2   κ0值在不同 vS30范围的分组统计

    Table  2   Group statistics of κ0 values in different vS30 ranges

    vS30/(m·s−1κ0
    数量最小值最大值标准差均值
    100—200370.019 60.076 70.014 850.054 48
    200—300650.022 20.072 60.014 580.051 33
    300—400950.009 20.077 20.015 380.043 41
    400—500950.008 80.087 50.016 060.040 43
    500—600730.003 80.067 40.015 160.035 60
    600—700460.003 80.080 90.015 170.029 28
    700—800240.013 00.071 10.014 110.033 48
    800—900190.004 90.073 90.017 840.034 12
    900—1000160.014 30.052 80.010 990.027 61
    1 000—1 10070.015 30.027 00.003 930.022 73
    1 100—1 20070.006 00.025 80.007 520.015 47
    1 200—1 300100.013 60.026 00.004 610.019 67
    1 300—1 40060.010 10.031 50.009 060.020 13
    1 400—1 50090.009 60.026 30.006 150.015 72
    1 500—1 800160.006 90.027 10.006 750.016 06
    1 800—2 100110.002 90.029 30.009 360.013 71
    2 100—2 400100.002 90.025 00.009 300.011 93
    下载: 导出CSV

    表  3   κrms vS30的经验关系

    Table  3   The empirical relationship between кrms and vS30

    拟合函数模型参数
    abSSER2c
    线性函数:κrmsavS30b −1.557×10−5 4.478×10−2 3.61×10−3 7.734×10−1
    多项式函数:κrms=${av^2_{{\rm{S}}30}} $+bvS30c 1.284×10−8 −4.741×10−5 6.559×10−4 9.587×10−1 5.886×10−5
    对数线性函数:κrmsalgvS30b −3.439×10−2 1.286×10−1 1.466×10−3 9.07 8×10−1
    双对数线性函数:lgκrmsalgvS30b −4.488×10−1 −2.72×10−1 2.245×10−3 8.587×10−1
    注:表中abc为模型拟合参数,SSE 表示和方差,R2为拟合优度。
    下载: 导出CSV
  • 傅磊,李小军. 2017. 龙门山地区的kappa(κ0)模型及汶川MS8.0地震的强地震动模拟[J]. 地球物理学报,60(8):2935–2947. doi: 10.6038/cjg20170803

    Fu L,Li X J. 2017. The kappa (κ0) model of the Longmenshan region and its application to simulation of strong ground-motion by the Wenchuan MS8.0 earthquake[J]. Chinese Journal of Geophysics,60(8):2935–2947 (in Chinese).

    郑旭,胡进军,谢礼立,周少辉,鹿子林. 2019. 基于山东省场地的沂沭断裂带地震动高频消减作用[J]. 哈尔滨工业大学学报,51(12):55–62. doi: 10.11918/j.issn.0367-6234.201901170

    Zheng X,Hu J J,Xie L L,Zhou S H,Lu Z L. 2019. High frequency reduction effect of ground motion in Yi-Shu fault zone based on sites of Shandong[J]. Journal of Harbin Institute of Technology,51(12):55–62 (in Chinese).

    朱百慧. 2016. 强震数据Kappa值影响因素的讨论[D]. 哈尔滨: 中国地震局工程力学研究所: 39−41.

    Zhu B H. 2016. Discussion on Influencing Factors of Kappa Value of Strong Earthquake Data[D]. Harbin: Institute of Engineering Mechanics, China Earthquake Administration: 39−41 (in Chinese).

    Anderson J G,Hough S E. 1984. A model for the shape of the Fourier amplitude spectrum of acceleration at high frequencies[J]. Bull Seismol Soc Am,74(5):1969–1993.

    Anderson J G,Humphrey J R. 1991. A least squares method for objective determination of earthquake source parameters[J]. Seismolo Res Lett,62(3/4):201–209.

    Anderson J G,Lee Y,Zeng Y H,Day S. 1996. Control of strong motion by the upper 30 meters[J]. Bull Seismol Soc Am,86(6):1749–1759.

    Atkinson G M. 1984. Attenuation of strong ground motion in Canada from a random vibrations approach[J]. Bull Seismol Soc Am,74(6):2629–2653.

    Boore D M. 1983. Stochastic simulation of high-frequency ground motions based on seismological models of the radiated spectra[J]. Bull Seismol Soc Am,73(6A):1865–1894.

    Boore D M,Thompson E M,Cadet H. 2011. Regional correlations of vS30 and velocities averaged over depths less than and greater than 30 meters[J]. Bull Seismol Soc Am,101(6):3046–3059. doi: 10.1785/0120110071

    Bora S S,Cotton F,Scherbaum F,Edwards B,Traversa P. 2017. Stochastic source,path and site attenuation parameters and associated variabilities for shallow crustal European earthquakes[J]. Bull Earthq Eng,15(11):4531–4561. doi: 10.1007/s10518-017-0167-x

    Cabas A,Rodriguez-Marek A,Bonilla L F. 2017. Estimation of site-specific kappa (κ0)-consistent damping values at KiK-net sites to assess the discrepancy between laboratory-based damping models and observed attenuation (of seismic waves) in the field[J]. Bull Seismol Soc Am,107(5):2258–2271. doi: 10.1785/0120160370

    Castellaro S,Mulargia F,Rossi P L. 2008. VS30:Proxy for seismic amplification?[J]. Seismol Res Lett,79(4):540–543. doi: 10.1785/gssrl.79.4.540

    Chandler A M,Lam N T K,Tsang H H. 2006. Near-surface attenuation modelling based on rock shear-wave velocity profile[J]. Soil Dyn Earthq Eng,26(11):1004–1014. doi: 10.1016/j.soildyn.2006.02.010

    Chang S C,Wen K L,Huang M W,Kuo C H,Lin C M,Chen C T,Huang J Y. 2019. The high-frequency decay parameter (Kappa) in Taiwan[J]. Pure Appl Geophys,176(11):4861–4879. doi: 10.1007/s00024-019-02219-y

    Drouet S,Cotton F,Guéguen P. 2010. VS30,κ,regional attenuation and MW from accelerograms:Application to magnitude 3-5 French earthquakes[J]. Geophys J Int,182(2):880–898. doi: 10.1111/j.1365-246X.2010.04626.x

    Edwards B,Fäh D,Giardini D. 2011. Attenuation of seismic shear wave energy in Switzerland[J]. Geophys J Int,185(2):967–984. doi: 10.1111/j.1365-246X.2011.04987.x

    Edwards B,Ktenidou O J,Cotton F,Abrahamson N,van Houtte C,Fäh D. 2015. Epistemic uncertainty and limitations of the κ0 model for near-surface attenuation at hard rock sites[J]. Geophys J Int,202(3):1627–1645. doi: 10.1093/gji/ggv222

    Hanks T C. 1982. fmax[J]. Bull Seismol Soc Am,72(6A):1867–1879. doi: 10.1785/BSSA07206A1867

    Huang M W,Wen K L,Chang S C,Chang C L,Liu S Y,Chen K P. 2017. The high-cut parameter (Kappa) for the near-surface geology in and around the Taipei basin,Taiwan[J]. Bull Seismol Soc Am,107(3):1254–1264. doi: 10.1785/0120160070

    Kilb D,Biasi G,Anderson J,Brune J,Peng Z,Vernon F L. 2012. A comparison of spectral parameter kappa from small and moderate earthquakes using southern California ANZA seismic network data[J]. Bull Seismol Soc Am,102(1):284–300. doi: 10.1785/0120100309

    Kishida T,Darragh R,Ktenidou O J. 2014. Fourier spectra and kappa (κ0) estimates for rock stations in the NGA-west2 project[J]. Eng Geol,202(4):1–13.

    Lai T S,Mittal H,Chao W A,Wu Y M. 2016. A study on kappa value in Taiwan using borehole and surface seismic array[J]. Bull Seismol Soc Am,106(4):1509–1517. doi: 10.1785/0120160004

    Laurendeau A,Cotton F,Ktenidou O J,Bonilla L F,Hollender F. 2013. Rock and stiff-soil site amplification:Dependency on vS30 and kappa (κ0)[J]. Bull Seismol Soc Am,103(6):3131–3148. doi: 10.1785/0120130020

    Lee V W,Trifunac M D. 2010. Should average shear-wave velocity in the top 30 m of soil be used to describe seismic amplification?[J]. Soil Dyn Earthq Eng,30(11):1250–1258. doi: 10.1016/j.soildyn.2010.05.007

    Silva W J, Wong I G, Darragh R B. 1998. Engineering characterization of earthquake strong ground motions in the Pacific Northwest[G]//Assessing Earthquake Hazards and Reducing Risk in the Pacific Northwest. Reston: U. S. Geological Survey: 313−324.

    Stanko D,Markušić S,Ivančić I,Mario G,Gülerce Z. 2017. Preliminary estimation of kappa parameter in Croatia[J]. IOP Conf Ser:Earth Environ Sci,95(3):032014.

    van Houtte C,Drouet S,Cotton F. 2011. Analysis of the origins of κ (kappa) to compute hard rock to rock adjustment factors for GMPEs[J]. Bull Seismol Soc Am,101(6):2926–2941. doi: 10.1785/0120100345

    van Houtte C,Ktenidou O J,Larkin T,Holden C. 2018. A continuous map of near-surface S-wave attenuation in New Zealand[J]. Geophys J Int,213(1):408–425. doi: 10.1093/gji/ggx559

图(7)  /  表(3)
计量
  • 文章访问数:  465
  • HTML全文浏览量:  209
  • PDF下载量:  115
  • 被引次数: 0
出版历程
  • 收稿日期:  2022-04-10
  • 修回日期:  2022-05-21
  • 网络出版日期:  2022-09-29
  • 刊出日期:  2023-10-29

目录

/

返回文章
返回