Discussions on the response spectral solution and fitting of spectrum-compatible artificial seismic waves
-
摘要: 为更高效合理地生成与设计反应谱相符的人造地震动信号，本文基于自主开发的地震动信号处理软件EQSignal对人造地震动生成过程中涉及的几个问题进行了探讨，并通过对单自由度简谐共振体系的响应分析，提出了具体的解决建议：① 反应谱求解时应该对高频段和低频段分别采用频域传递函数法和逐步积分法求解；② 反应谱周期控制点的分布模式也应分段处理，短周期段宜采用对数平均分布，长周期段宜采用线性平均分布，反应谱总周期控制点不宜少于120个；③ 人造地震动反应谱与设计反应谱拟合的过程中，使用频域方法与时域方法相结合可兼顾效率与收敛性。Abstract: In order to generate design-spectrum-compatible artificial seismic waves more efficiently and rationally, three problems that arise during the procedure of artificial seismic signals generation are studied and discussed in this paper by the aid of an earthquake signal processing software named EQSignal which is developed by the authors, and through the response analysis of harmonic resonance system with single degree of freedom, the corresponding suggestions to solve these questions are also provided. It is suggested that transfer function method should be adopted within the high-frequency segment of a response spectrum, and step-by-step integration method can be used within the low-frequency segment. There should be more than 120 period points in a seismic response spectrum, and the period points should be log-uniformly distributed at short-period segment and uniformly distri-buted at long-period segment. During seismic spectrum fitting, both the two iterative methods in frequency domain and time domain are recommended to use in order to give consideration to both efficiency and convergence.
-
引言
覆盖土层对地震波传播特性具有重要影响（高武平等，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. 水平成层场地的集中质量切比雪夫谱元模型
1.1 水平成层场地简化模型
将水平成层场地简化为如图1所示的计算模型。设坐标原点位于自由地表，建立正方向竖直向下的一维坐标系z。假设有N−1个土层覆盖于基岩半空间上，各覆盖土层的厚度、质量密度和剪切模量分别为hi，ρi和Gi，i = 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 site1.2 切比雪夫谱单元位移模式
一维土层分析的切比雪夫谱单元以 ［ −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ξ），ci和ck为多项式系数，当i＝0或n时，取值为2，当i＝1，···，n－1时，取值为1。
1.3 切比雪夫谱单元特性矩阵
利用切比雪夫谱单元对水平成层场地模型进行空间离散后，得到离散方程
$$ {\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按照单元编号组装而成。Me，Ce和Ke计算如下：$$ {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 ±1 1 2 ±1, 0 0.333 3, 1.333 3 3 ±1, ±0.5 0.111 1, 0.888 9 4 ±1, ±0.707 1, 0 0.066 7, 0.533 3, 0.8 5 ±1, ±0.809 0, ±0.309 0 0.04, 0.360 7, 0.599 3 利用式（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) 利用以上性质，Ce和Ke可精确计算为
$$ 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) 2. 时域逐步积分
对于土层地震反应分析问题，一般将入射波作为边界处的运动形式施加，因此需要将模型内域节点和边界节点分开处理。在本文建立的水平成层场地谱元模型中，边界节点仅有一个，即模型底部人工边界上的单元节点，其余节点均属于内域节点。将式（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}}}$ 分别为边界节点的加速度、速度和位移；MI和MB为总体质量矩阵M按照内域节点和边界节点的节点编号进行分块后的子矩阵；CI，CIB，CBI和CB以及KI，KIB，KBI和KB分别为总体阻尼矩阵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，并采用中心差分法进行时步积分。在t＝p∆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通常不会呈对角矩阵形式，故由MI与CI组合的矩阵亦会是非对角阵。此情形下可考虑构建仅同质量矩阵呈比例的瑞雷阻尼矩阵，或者将阻尼矩阵对角化（Thomson et al，1974）等措施。
3. 基于透射边界的地震动输入
由针对节点的总波场给出的方程式（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）所定义的单元位移模式，插值系数a0，a1，···，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可保证数值结果始终稳定。
4. 算例分析
4.1 均匀半空间土体地震反应分析
为了检验本文方法的正确性，首先对一个均匀半空间场地模型进行分析，根据行波理论得到其解析解。该模型土体厚度为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 incidence4.2 Kik-net台站记录模拟
利用日本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 Ⅳ 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 表 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 图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造成本文计算结果与实际台阵记录存在一定差异的原因可能有：① 本文采用的一维水平成层场地模型是一种高度简化的计算模型，无法反映出实际三维场地的地形特征，因此一维波动的数值计算结果一般无法体现出对实际场地地震动有重要影响的地形放大效应；② 实际场地并不一定仅受剪切波作用，还有可能受到纵波、面波等多种地震波的作用，并且地震波的传播方向也不一定恰好垂直向上的；③ 本文方法未考虑土体的非线性特性。
5. 结论
本文提出了一种求解水平成层场地地震反应的时域集中质量切比雪夫谱元法。通过节点积分法严格地导出集中质量矩阵，克服了传统切比雪夫谱元法由于质量矩阵为非对角矩阵形式所带来的计算效率不高的问题。采用中心差分法进行时间积分，并嵌入多次透射人工边界，形成了高效的时域显式算法。与传统有限单元法相比，该方法在每个波长内仅需布置一个谱单元即可获得令人满意的计算精度，显著降低了对空间网格尺寸的要求.
选择日本Kik-net强震台网中分属四种不同场地类型的台站记录检验了本文方法处理实际场地的能力。计算结果表明，本文方法对于Ⅰ 1，Ⅱ和Ⅳ类场地在较弱地震和中等强度地震作用下能够较好地预测地面运动特征，但对于Ⅲ类场地及强震作用下的地表反应计算结果与观测结果相比仍存在一定误差，后续工作中可考虑添加土体的时域非线性本构关系.
-
![]()
图 1 不同算法求解单自由度体系共振响应的加速度时程响应精度对比
（a） 数值解与精确解的相关系数C；（b） 数值解与精确解的峰值比P
Figure 1. Accuracy comparison of acceleration history by different algorithms for solution of resonant response of a single-degree-of-freedom system
（a） Correlation coefficient between numerical and exact solutions；（b） Peak ratio to numerical and exact solutions
![]()
图 2 单自由度体系的共振响应时程曲线
Figure 2. Time history curves of resonant response of a single-degree-of-freedom system
![]()
图 3 不同反应谱周期控制点分布模式下加速度反应谱的对比
（a） x轴刻度为线性刻度；（b） x轴刻度为对数刻度
Figure 3. Comparison of acceleration response spectrum on the condition of distribution patterns of different period controling points
（a） x axis is linear scale；（b） x axis is logarithmic scale
![]()
图 4 周期控制点数量对CHICHI04-CHY039-E （a）和CHICHI-ILA004-N （b）时程的反应谱曲线的影响
Figure 4. Influence of period contorling point counts on the shape of response spectrum curve of the histories CHICHI04-CHY039-E （a） and CHICHI-ILA004-N （b）
![]()
图 5 周期控制点加密前后反应谱对照图
（a） 20个拟合控制点；（b） 60个拟合控制点；（c） 120个拟合控制点
Figure 5. Comparison of response spectra before and after increasing the number of period controling points
（a） Period point counts 20； （b） Period point counts 60； （c） Period point counts 120
![]()
图 6 反应谱拟合所用周期控制点数对拟合结果可靠性的影响
（a） 平均误差；（b） 最大误差
Figure 6. Influence of period controling point count used in response spectral fitting on the reliability of fitting results
（a） Average error；（b） Maximum error
![]()
图 7 频域反应谱拟合过程中的时程波形变化
（a） 原始信号；（b） 5次迭代后的波形；（c） 20次迭代后的波形
Figure 7. Variation of time-history waveform during frequency-domain response spectral fitting
（a） Original signal；（b） Waveform after five iterations；（c） Waveform after 20 iterations
![]()
图 8 时域反应谱拟合过程中的时程波形变化
（a） 原始信号；（b） 5次迭代后波形；（c） 20次迭代后波形
Figure 8. Variation of time-history waveform during time-domain response spectral fitting
（a） Original signal；（b） Waveform after five iterations；（c） Waveform after 20 iterations
![]()
图 9 反应谱拟合前后的傅里叶幅值谱（左）和相位差谱（右）
（a） 原始信号；（b） 频域拟合结果；（c） 时域拟合结果
Figure 9. Fourier amplitude spectrum （left） and phase difference spectra （right） before and after response spectral fitting
（a） Original signal；（b） Frequency-domain fitting result；（c） Time-domain fitting result
表 1 反应谱拟合所用周期控制点数对拟合结果的影响
Table 1 Influence of period controling points number on fitting results of the response spectrum
N 反应谱拟合平均误差 反应谱拟合最大误差 M=30 M=60 M=90 M=120 M=150 M=30 M=60 M=90 M=120 M=150 30 1.7% 6.2% 6.9% 6.7% 6.6% 5.7% 21.5% 23.5% 24.1% 20.1% 60 2.7% 2.3% 5.0% 4.9% 4.5% 7.8% 5.7% 20.8% 21.5% 17.0% 90 2.2% 3.8% 2.0% 3.6% 3.6% 6.3% 14.4% 5.2% 14.4% 11.6% 120 1.9% 2.4% 2.2% 1.9% 2.4% 5.0% 8.5% 7.3% 5.5% 9.9% 150 1.9% 2.4% 2.2% 2.3% 1.8% 4.9% 6.6% 6.8% 7.6% 5.8% 注：M表示误差分析周期控制点数，N表示反应谱拟合周期控制点数。 
下载: 导出CSV
表 2 地震动信号利用频域法（F）和时域法（T）进行反应谱拟合的实例对比
Table 2 Response spectral fitting examples of seismic waves by using frequency-domain method（F）and time-domain method（T）
地震动信号 迭代次数 运算耗时/s 最大误差 平均误差 是否收敛 F T F T 拟合前 拟合后 拟合前 拟合后 F T F T F T CHICHI04-CHY039-E 30 23 0.52 49.4 90.6% 5.9% 4.4% 47.2% 1.8% 0.7% 否 是 CHICHI04-CHY039-N 30 30 0.50 64.1 122.4% 8.6% 7.4% 49.9% 2.5% 1.8% 否 否 LOMAP-A10090 30 24 0.31 26.4 94.5% 12.4% 4.5% 49.0% 3.2% 1.2% 否 是 LOMAP-A10000 30 15 0.29 16.5 72.0% 8.4% 3.9% 42.8% 2.6% 1.7% 否 是 MORGAN-A01310 30 29 0.29 30.8 98.1% 11.4% 4.7% 68.3% 2.8% 1.6% 否 是 MORGAN-A01040 30 23 0.30 24.9 97.1% 24.7% 4.8% 67.3% 6.6% 1.4% 否 是 YOUNTVL-0436c090 30 23 0.46 46.7 111.7% 16.3% 3.9% 73.2% 3.7% 1.5% 否 是 YOUNTVL-0436a360 30 11 0.45 22.2 97.0% 15.5% 4.5% 69.9% 3.8% 1.3% 否 是 CHICHI-ILA004-W 30 19 1.09 118.1 431.7% 9.4% 4.2% 128.4% 1.5% 1.6% 否 是 CHICHI-ILA004-N 30 23 1.10 142.0 485.6% 7.0% 4.2% 136.6% 1.5% 1.4% 否 是 LOMAP-A02133 30 30 0.29 32.3 91.6% 13.6% 5.3% 65.4% 4.6% 1.3% 否 否 LOMAP-A02043 30 30 0.30 33.1 105.5% 7.3% 12.8% 63.3% 2.5% 3.1% 否 否 LOMAP-MEN360 30 26 0.27 26.9 91.3% 15.9% 4.7% 58.8% 4.8% 1.4% 否 是 LOMAP-MEN270 30 30 0.28 31.6 82.1% 5.6% 8.9% 46.6% 1.9% 2.9% 否 否
下载: 导出CSV
-
郭子雄, 王妙芳. 2006. 人造地震动合成的研究现状及展望[J]. 华侨大学学报: 自然科学版, 27（1）: 7-11. Guo Z X, Wang M F. 2006. State of the art and prospect of artificial ground motion[J]. Journal of Huaqiao University: Natural Science, 27（1）: 7-11 （in Chinese）.
胡聿贤. 2006. 地震工程学[M]. 北京: 地震出版社: 8–12. Hu Y X. 2006. Earthquake Engineering[M]. Beijing: Seismological Press: 8–12 （in Chinese）.
胡聿贤, 何训. 1986. 考虑相位谱的人造地震动反应谱拟合[J]. 地震工程与工程振动, 6（2）: 37-51. Hu Y X, He X. 1986. Phase angle consideration in generating response spectrum-compatible ground motion[J]. Earthquake Engineering and Engineering Vibration, 6（2）: 37-51 （in Chinese）.
全伟, 李宏男. 2008. 调整已有地震动拟和规范反应谱人造地震动方法比较[J]. 防灾减灾工程学报, 28（1）: 91-97. Quan W, Li H N. 2008. Comparison of methods of generation of spectrum-compatible artificial earthquake through the modification of available records[J]. Journal of Disaster Prevention and Mitigation Engineering, 28（1）: 91-97 （in Chinese）.
徐科军. 2001. 信号处理技术[M]. 武汉: 武汉理工大学出版社: 5–127. Xu K J. 2001. Signal Processing Technology[M]. Wuhan: Wuhan University of Technology Press: 5–127 （in Chinese）.
中华人民共和国住房和城乡建设部. 2010. GB 50011—2010建筑结构抗震设计规范[S]. 北京: 中国建筑工业出版社: 33–35. Ministry of Housing and Urban-Rural Construction of the People’s Republic of China. 2010. GB 50011—2010 Code for Seismic Design of Buildings[S]. Beijing: China Building Industry Press: 33–35 (in Chinese).
大崎順彦. 1980. 新•地震动の ス ぺタ ト ル解析入门[M]. 东京都: 鹿岛出版会. Ozaki N. 1980. Introduction to Spectral Analysis of Seismic Ground Motion[M]. Tokyo: Kashima Press （in Janpanese）.
Al Atik L, Abrahamson N. 2010. An improved method for nonstationary spectral matching[J]. Earthquake Spectra, 26（3）: 601-617.
Chopra A K. 2007. Dynamics of Structures: Theory and Applications to Earthquake Engineering[M]. 3rd edition. New Jersey: Prentice Hall: 1–812.
Clough R W, Penzien J, Griffin D S. 1990. Dynamics of Structures[M]. New Jersey: Prentice Hall.
Gao Y, Wu Y, Li D.2013. An improved method for the generating of spectrum compatible time series using wavelets[J]. Earthquake Spectra, 30（4）: 1467-1485.
Hancock J, Watson-Lamprey J, Abrahamson N A, Bommer J J, Markatis A, Mccoy E, Mendis R. 2006. An improved method of matching response spectra of recorded earthquake ground motion using wavelets[J]. J Earthquake Eng, 10（S1）: 67-89.
Housner G W, Martel R R, Alford J L. 1953. Spectrum analysis of strong-motion earthquakes[J]. Bull Seismol Soc Am, 43（2）: 97-119.
Pan C, Zhang R F. 2016. EQSignal: A useful tool to process and generate earthquake signals[EB/OL]. [2016-08-01]. www.eqsignal.com.
Pan C, Zhang R F, Luo H, Shen H. 2017. Target-based algorithm for baseline correction of inconsistent vibration signals[J]. Journal of Vibration and Control, 24（12）: 2562-2575.
Rizzo P C, Shaw D E, Jarecki S J. 1975. Development of real/synthetic time histories to match smooth design spectra[J]. Nucl Eng Des, 32（1）: 148-155.
Tsai N C. 1972. Spectrum compatible motions for design purposes[J]. J Eng Mechans Divis, 98（2）: 345-356.
-
期刊类型引用(7)
1. 陈珍，郝冰，李远东，周正华，卞祝，韩轶. 含软弱土层场地地震动加速度反应谱特征周期调整方法. 地震学报. 2024(04): 734-750 . 
本站查看
2. 杨兰兰，傅梓岳，王登峰，XIE Weichau. 基于数字滤波技术拟合规范反应谱的地震动研究. 振动与冲击. 2023(09): 57-67+105 . 百度学术
3. 李晓莉，赵月敏，邹积娜，王东升. 基于小波变换的长周期桥梁选波方法研究. 中国地震. 2023(04): 795-809 . 百度学术
4. 禹海涛，李晶，王祺. 软土隧道基于地震响应的输入地震动排序. 地震学报. 2022(01): 123-131 . 本站查看
5. 刘帅，夏舟，张悦超. 软土场地大型LNG储罐组合滑移隔震. 吉林大学学报(工学版). 2022(04): 856-864 . 百度学术
6. 谢皓宇，郑万山，仉文岗，高文军. 考虑迭代相关及相位谱的人工地震波反应谱拟合. 地震学报. 2020(03): 341-348+378 . 本站查看
7. 刘帅，潘超，周志光. 耗能联肢墙体系的减震性能及参数影响. 浙江大学学报(工学版). 2019(03): 492-502 . 百度学术
其他类型引用(5)
计量
- 文章访问数: 1647
- HTML全文浏览量: 725
- PDF下载量: 82
- 被引次数: 12
