Study on Seismic Response of Underground Structure in Saturated Soil Deposit Based on OpenSees
-
摘要: 应用OpenSees有限元计算程序,选取某地铁车站结构为主要分析对象,构建了饱和土-地下结构体系地震反应运算模型,并利用该模型进行系统的地震反应计算。将现场饱和土动力反应视为饱和两相介质近场波动问题,选取时域显式数值算法进行计算,同时考虑了土体的弹塑性。研究结果显示:(1)因选取弹塑性土体本构,土体-地下结构的位移反应时程和输入地震动的位移时程体现出了显著性差异。(2)对于两层三跨地下结构,在以剪切波形式输入的地震动作用下,顶板的峰值加速度和侧向位移最大,中板次之,底板最小,场地土体对地震波具有放大效应,顶层的层间位移大于底层。(3)在地震动作用下,地下结构不同区域应力时程的改变规律存在非常显著的差异。中柱底部与底板节点处的应力峰值最大。地震动输入结束时,结构存在残余应力。Abstract: Using the finite element calculation software OpenSees, this study focuses on a subway station structure as the primary analysis object while establishing an operational seismic response model for a saturated soil-underground structural system. This model enables a systematic calculation of the seismic response. The study considers the dynamic response of the saturated soil on site as a near-field fluctuation problem of a saturated two-phase medium. A time-domain explicit numerical algorithm is employed for calculations, incorporating the elastic-plastic dynamic response characteristics of the soil. The results of this study indicate several key findings: (1) The elastic-plastic dynamic response characteristics of the site soil are accounted for during the simulation analysis. This consideration reveals a significant difference between the displacement response time histories of the soil and underground structures compared to the displacement time histories of the input ground motion. (2) In the case of a two-story, three-span underground structure subjected to ground motion in the form of shear waves, the roof exhibits the highest peak acceleration and lateral displacement, followed by the middle plate, while the bottom plate shows the smallest values. Additionally, the inter-story displacement of the top layer is greater than that of the bottom layer. (3) Under earthquake loading, there are pronounced differences in the stress time history patterns across various regions of the underground structure. The stress peak at the joint between the bottom of the center column and the bottom plate is the highest observed. Furthermore, when ground motion input ceases, residual stresses remain in the structure.
-
引言
我国是世界上地震活动强烈和地震灾害严重的国家之一,地震带来的人员伤亡、经济损失巨大,灾害发生后的应急救援工作十分重要。目前,受灾区范围大、破坏严重、地形复杂、天气恶劣等条件的限制,传统地面现场调研方法不能满足震后应急需求,因此各种遥感影像广泛应用于震后救援与应急中(宿渊源等,2012;杨乐等,2012),使遥感技术成为震害提取和损失评估的重要手段(陈鑫连,1995;王晓青等,2008;姜立新等,2012)。大部分光学影像不能提供高程信息,个别可提供高程信息,但精度不高,制约了对震害信息的提取,而机载激光雷达技术则可以弥补此不足。机载激光雷达扫描(Light Detection and Ranging,下文简称LiDAR)技术能够提供高精度的三维空间信息,能增强对受损建筑物的识别。同时,与光学遥感测量相比,机载激光雷达测量受日照和云雾天气条件的影响较小,基本可实现全天候观测。在2008年汶川8.0级地震(马洪超等,2008)、2010年海地7.0级地震(Margesson等,2010)、2011年3月11日日本地震海啸(Fritz等,2011)等地震的应急响应阶段,应急人员均利用LiDAR技术快速获取了灾区的高精度空间数据。特大地震灾害应急响应实践证实,机载LiDAR技术在震后交通不畅、地形复杂、气象条件差的地区,具备灾情数据获取能力。因此,基于震后机载激光雷达点云进行地物识别,尤其是对受损房屋建筑进行提取具有重要意义。
目前,机载LiDAR技术在地震应急领域的应用日益广泛,主要研究方法可分为基于不同时相的激光雷达点云信息提取、LiDAR点云与其他遥感影像相结合的信息处理和单时相LiDAR点云信息提取。当前大部分地区的震前LiDAR数据匮乏,基于不同时相的点云信息提取在地震应急过程中较难实现。而LiDAR点云与高分辨率卫星影像相结合(程亮等,2008;窦爱霞等,2013),或者与航空影像相结合(张志友,2008;邓非等,2010;Trinder等,2012),综合运用遥感影像的光谱信息和LiDAR数据的三维信息,精确提取地物目标的方法已被广泛应用,但存在不同数据融合和匹配困难、易引入分辨率差异等缺点。利用震后单时相机载LiDAR数据则可避免多时相LiDAR数据的数据源匮乏、多源数据的匹配融合等问题,较适合应用于地震应急救援。Dorninger等(2008)提出了基于机载LiDAR点云,通过建筑物轮廓、面片的提取实现城市建筑物三维自动建模的方法;Labiak等(2011)利用移动窗口的方法对点云多个特征因子进行统计,得到栅格影像,以尽可能识别出树木,但仅用坡度特征因子对受损建筑物进行提取和定量化,忽视了其他特征对受损建筑物的识别能力。目前,基于震后单时相机载LiDAR数据对建筑物震害识别开展的研究较少,且研究方法多面向单个点的特征。针对单个点特征的研究可以避免点云划分,即同一地物对象点云聚合,但没有考虑地物对象整体点云的特征关系。
本文在前人研究的基础上,提出反映地物对象点云整体特征的回波次数特征指标,并结合回波强度、法向量、坡角及高差等已有的机载LiDAR点云特征,以海地地震部分灾区为研究区,分析倒塌、未倒塌建筑物以及树木等的LiDAR点云的各特征均值和标准差分布情况,挑选有效的分类敏感特征,再利用K-最近邻分类方法对测试样本进行建筑物震害程度的分类提取与验证,以探讨利用LiDAR数据对未倒塌建筑物、倒塌建筑物和树木进行识别的效果。
1. 方法
1.1 分类特征
通常机载LiDAR点云信息除点的三维坐标信息外,还包括点的回波强度、回波次数信息。同时,同一地物对象包含一系列的点云,其中某一点邻域内的点集,将从空间形态、回波强度和回波次数的统计分布上反映所在地物对象的特征。因此,本文提出了回波次数比(R)点云特征量,同时,引入回波强度(I)、归一化强度(IN)、最邻近点高差(H)、法向量与天顶方向夹角(A)和X向坡角(SX)、Y方向坡角(SY)等,作为建筑物震害LiDAR点云分类特征的参数。
(1)回波次数比(R)
设研究区域点云所在的投影坐标系为O-XYZ,则LiDAR点云中某一点i的Ri为i点在水平投影面O-XY上,n个最近邻的投影点(如图 1)中,回波次数大于1的点所占的比例:
$$ {R_i} = \frac{{{m_i}}}{n} $$ (1) 式中,mi为回波次数大于1的投影点的数量。
(2)回波强度(I)与归一化回波强度(IN)
点云回波强度(I)取决于介质表面的反射系数,而激光的波长、介质材料及介质表面的明暗黑白程度等都是影响反射系数的因素。反射介质表面越亮,反射率越高(张小红,2007)。理论上未倒塌建筑物的屋顶表面多平滑,粗糙度小,相对于倒塌建筑物和树木的离散组成结构和粗糙表面而言,同等激光波长下其反射系数更大,回波强度更大。
为表征不同地物回波强度的分布特征,本文对回波强度采用归一化处理(Labiak等,2011),即某地物对象g的点云集合中,某一点i的归一化强度(INi)等于该点的回波强度(Ii)与该地物的最小强度值之差与地物的回波最大、最小强度值之差的比值:
$$ {I_{Ni}} = \frac{{{I_i} - {I_s}}}{{{I_l} - {I_s}}} $$ (2) 式中,Is、Il分别为地物对象g点云中的回波强度最小值和最大值。
(3)最邻近点高差(H)
最邻近点高差(H)指目标点高程与其水平投影面上距离其最近的点的高程差值。显然未倒塌建筑物单体点云高差变化极小,理想状态下结构简单的平面屋顶或人字形屋顶的点云高差为定值,对于屋顶结构复杂的未倒塌建筑物,其高差会有所变化,但仍呈现较强的规律性;而倒塌建筑物单体高差则出现多值性,规律性也较差。因此,高差对地物的识别有一定作用。
(4)法向量夹角(A)
法向量夹角(A)是指法向量与天顶方向的夹角(黄树松等,2016),其范围是[0°,180°],其中法向量通过Hoppe等(1992)提出的基于局部表面拟合的方法来确定。本文通过搜索三维空间中距目标点最近的20个点构成最小二乘平面,再求取该平面在目标点处的法向量,从而得到法向量与天顶方向的夹角。法向量夹角对建筑物的倒塌部分十分敏感(Dou等,2016)。
(5)X方向坡角(SX)与Y方向坡角(SY)
在机载LiDAR三维点云所在的地理坐标系统(O-XYZ)中将点云投影在O-XY平面上,求取X方向上与目标投影点最邻近的投影点,目标点与该投影点对应的LiDAR点所构成的向量与O-XY平面的夹角称为目标点的X方向坡角(SX)。同理可以确定Y方向坡角(SY),如图 2。为便于统计,本文规定坡角范围是[-90°,90°]。通常,未倒塌建筑物单体点云的X方向坡角和Y方向坡角的大小均一,分布规律;倒塌建筑物单体坡角值则大小起伏,倒塌与未倒塌部分交界处有明显变化,一般倒塌建筑物的坡角比未倒塌建筑物的坡角大,呈现局部规则分布;而树木坡角值则大小不一,分布无规律(黄树松等,2016)。
1.2 K-最近邻分类法
K-最近邻分类法是惰性学习法,基于类比学习,其基本原理是将测试元组与已知类别的测试元组进行比较学习。如图 3所示(以二维特征为例),对具有两个特征的数据组利用K-最近邻分类时,我们计算待分类点Xu距离K个邻近点的距离,则距离最小的点所属的类别即为该待分类点的类别(ω1、ω2或ω3)。K-最近邻分类法最大的特点是不需要使用训练样本设计分类器,可直接用训练样本对测试样本进行分类,原理简单,便于实现,应用广泛(滕敏等,2015)。目前,对于分类速度的改进主要有浓缩训练样本集算法和未知样本预分类算法,其基本思想是结合聚类算法或其他分类器,先对原始样本进行预分类,选择预分类后的代表子集参与K-最近邻分类,实现样本压缩,提高分类速度;对于分类准确度的改进主要通过选用合适的距离机制等方法实现(王娜等,2009)。
该分类算法的数据类型可以是数值型和分类型,数值型可以用于计算距离,但分类型则不适合,因此针对分类型数据,通常对其进行数字化再比较两个样本中对应属性值的差值。在本研究中,所有样本的属性值均为数值型,可直接进行距离计算,对输出类别进行数字化。
K-最近邻分类法中度量距离的选择也很重要,最常用的是欧几里得距离,不同情况下也可选用Manhattan距离或其他距离度量(Han等,2007)。本研究使用欧几里得距离,假设某测试验本为Y(y1,y2,…,yn),训练样本为X(x1,x2,…,xn),计算方法如下:
$$ d = \sqrt {\sum\limits_{i = 1}^n {{{\left({{y_i} - {x_i}} \right)}^2}} } $$ (3) 其中,d为两样本之间的欧几里得距离。
如此计算该测试样本距离邻近K个训练样本的距离,找到其中距离的最小值及其对应的训练样本,则测试样本与该训练样本类别相同。依此对所有测试样本进行类别判断,完成分类。
2. 研究数据及其预处理
2.1 研究数据
2010年1月12日海地首府太子港附近发生7.0级地震,震源深度约10km,后发生多次余震,余震最大达6.3级,造成包括海地总统府、国会大厦在内的大多数建筑物损毁和几十万人员伤亡。世界银行全球减灾和灾后恢复机构(Globe Facility for Disaster Deduction and Recovery,简称GFDDR)于震后获取了灾区机载LiDAR点云,数据获取时间为2010年1月21日,点云密度为2—5点/m2。本文选取位于海地首都太子港城区中约0.8km2的区域作为研究区。辅助使用的光学遥感影像是同机相机获取的光学航空影像。该区域房屋建筑类型多样,有大型多层建筑区,如公寓、宾馆、教堂等,建筑面积较大,分布较稀疏,其间夹杂少量树木或小范围的树林;亦有低矮简易房屋,建筑物单体面积较小,分布较密集,在卫星影像上成片分布,树木较少,无植被密集区。该研究区域基本涵盖了城区所有房屋的建筑类型,具有代表性。利用光学影像在研究区选取了210个典型样本,作为本文特征选择、分类及验证的样本。研究中随机选取了未倒塌建筑物、倒塌建筑物和树木各50个样本作为训练样本,分布如图 4(a)和4(b);各20个样本作为测试样本,分布如图 4(c)所示。
2.2 数据预处理
获取点云数据的过程中,由于扫描仪本身的缺陷、环境干扰(如空中飞鸟)等因素,点云中常包含少量噪声点,剔除这些噪声点的过程称为点云去噪。本文基于点云魔方软件,采用基于点云局部空间分布统计的去噪算法,将局部点密度与整体点密度相差较大的点标记为噪声点并去除。为获得非地面点,本文采用CSF(Cloth Simulation Filter)滤波方法(Zhang等,2016),对去噪后的点云进行滤波,以便从激光脚点数据点云中提取数字地面模型(DTM)和数字高程模型(DEM)(张小红,2007),即提取地面点。CSF采用布料模拟原理,将原始点云倒置,再将具有一定硬度的布料覆盖其上,得到近似的地形估计,最后依据估计的地形,分离地面点与非地面点。以图 5(a)研究区1的点云为例进行实验分析,滤波所得结果如图 5(b)所示,可见该方法基本分离了地面点与非地面点,以CSF滤波法后所得的非地面点作为后续研究的基础。
3. 基于机载LiDAR点云参数的震害分类特征分析
3.1 震害点云参数特征分析
本文选取了回波次数比(R)、回波强度(I)、归一化回波强度(IN)、最邻近点高差(H)、法向量夹角(A)、X方向坡角(SX)和Y方向坡角(SY)等7个因子,对未倒塌建筑物、倒塌建筑物和树木的训练样本的点云数据进行了统计分析,结果如图 6所示。可见,不同类型的地物,其各因子的分布有明显差异。其中未倒塌建筑物的回波次数比的值基本为0;倒塌建筑物多大于0,部分点的R值在0.1—0.2范围内;树木R值多大于0.1,最高可达0.6。关于回波强度,未倒塌建筑物屋顶表面多平滑,粗糙度小,相对于倒塌建筑物和树木的离散组成结构和粗糙表面而言,同等激光波长下其反射系数更大,回波强度更大。从图 6中可见,未倒塌建筑物的最邻近点高差值H基本小于0.1m;倒塌建筑物的倒塌部分高差发生明显变化,H值可达5m;而树木的H值为0—16m,变化无规则。3类地物的法向量夹角(A)和X、Y方向坡角变化也比较明显,未倒塌建筑物角度呈现规则变化,倒塌建筑物的角度分布范围较大,但有明显的变化边界,树木则大小不一,分布无规律。
3.2 震害点云特征选择
为实现不同建筑物震害的识别,对选取的150个训练样本和60个测试样本分别计算了各样本点云的上述7个特征因子值及其均值和标准差,以进一步从统计上研究不同类型地物上述各特征因子的分布情况,选择适宜的分类特征。未倒塌建筑物、倒塌建筑物和树木样本的各个因子均值和标准差的概率密度曲线分布如图 7所示,其中各地物类别的回波强度、归一化回波强度的均值和标准差、X方向和Y方向坡角均值的曲线交叉严重,重叠部分占比较大,
说明这些特征值对地物的可分性较差,各类别其他8个特征值的曲线交叉部分相对较少,即有较强的可分性。因此,本文选取回波次数比均值和标准差、最邻近点高差均值和标准差、法向量夹角均值和标准差、X方向和Y方向坡角标准差作为分类特征。
4. 分类结果与讨论
基于上节选取的回波次数比均值、标准差等8个分类特征,利用150个训练样本,60个测试样本(其中未倒塌建筑物、倒塌建筑物、树木的训练样本各50个,测试样本各20个,类别编码依次为1、2、3),通过K-最近邻分类法进行地物自动分类。K-最近邻分类算法基于Matlab平台实现,流程如图 8所示。
为了确定最优的K值,分别测试了K为1—100的分类正确率,分布如图 9所示,当K=3时,分类准确率为98.3%,达到最高。所有分类准确率均在85%以上,这表明回波次数比均值和标准差、最邻近点高差均值和标准差、法向量夹角均值和标准差、X方向和Y方向坡角标准差这8个特征是区分未倒塌建筑物、倒塌建筑物和树木的有效分类特征。
统计K取1—100的分类结果,发现错误最多的为未倒塌建筑物和倒塌建筑物的分类,如图 10中黄色圈内所示,有2个未倒塌样本错分为倒塌样本,一个倒塌建筑物错分为树木,2个倒塌建筑物错分为未倒塌建筑物。树木样本在每次分类中全部分类正确,这主要是由样本选取精度及样本本身质量造成的。
本文中未倒塌建筑物、倒塌建筑物和树木样本全部为手动选取,因此分类方法的自动性有待探究,如何实现滤波后非地面点云的自动分割是未来研究的重要内容。本文中对建筑物仅划分了未倒塌建筑和倒塌建筑物,若能对倒塌建筑物进行更细致的分类,将进一步提高震害分类精度,为灾后应急救援工作提供帮助。此外,各类地物的测试样本数量有限,因此该分类方法的普适性仍有待提高。研究中所使用的K-最近邻分类法方法在分类过程中,最邻近样本数K的选取直接影响分类效果,因此对不同的K值进行试验也是研究的重要工作。
5. 结论
本文分析了未倒塌建筑物、倒塌建筑物和树木3类地物的机载LiDAR点云特征,选取了回波强度、回波次数比、最邻近点高差、法向量夹角及坡角特征,计算其均值和标准差并统计各样本特征值分布,选择其中分割效果最好的8个分类特征,在此基础上利用K-最近邻分类法进行分类试验。结果表明,基于回波次数比均值和标准差、最邻近点高差均值和标准差、法向量夹角均值和标准差、X方向和Y方向坡角标准差这几种分类特征,利用K-最近邻分类法所得结果与真实分类具有非常高的一致性,证明了本研究中所选分类特征的有效性,同时利用该方法可以有效实现地物点云分类,能够为震后应急救援工作提供有效指导,为灾后经济损失评估提供依据。
-
表 1 场地土体物理力学参数
Table 1. Site soil physical and mechanical parameters
土层 重度/(kN·m−3) 内摩擦角/(°) 弹性模量/MPa 泊松比 孔隙率 素填土 19.0 16 1.0 0.4 — 细砂 19.0 30 7.5 0.3 0.474 黏土 19.3 16 13.2 0.42 — -
陈少林,朱学江,赵宇昕等,2019. 考虑土骨架非线性的饱和土-结构相互作用分析. 地震工程与工程振动,39(1):114−127.Chen S. L., Zhu X. J., Zhao Y. X., et al., 2019. Analysis of saturated soil-structure interaction considering soil skeleton nonlinearity. Earthquake Engineering and Engineering Dynamics, 39(1): 114−127. (in Chinese) 程学磊,李文东,海然等,2022. 地震作用下饱和软土场地地下结构动力参数敏感性模拟分析. 广西大学学报(自然科学版),47(1):92−102.Cheng X. L., Li W. D., Hai R., et al., 2022. Parametric sensitivity simulation analyses of subway station structure surrounded by saturated soft soil foundation under earthquake. Journal of Guangxi University (Natural Science Edition), 47(1): 92−102. (in Chinese) 崔智谋,2012. 基于流固耦合动力模型的饱和两相介质-地下结构动力相互作用研究. 北京:北京工业大学.Cui Z. M., 2012. Study of dynamic interaction between fluid-saturated porous media of fluid-soild coupling model and underground structure. Beijing:Beijing University of Technology. (in Chinese) 丁伯阳,宋宥整,2019. 饱和土地下源u-P形式解答动力响应计算. 岩土力学,40(2):474−480.Ding B. Y., Song Y. Z., 2019. Dynamic response calculation for u-P solution in saturated soil subjected to an underground point source. Rock and Soil Mechanics, 40(2): 474−480. (in Chinese) 丁海滨,管凌霄,童立红等,2023. 基于非局部Biot理论的循环荷载下饱和土地基动力特性研究. 工程力学,40(3):141−152.Ding H. B., Guan L. X., Tong L. H., et al., 2023. On investigating the dynamic characteristics of saturated soil foundation subjected to cyclic load based on nonlocal Biot theory. Engineering Mechanics, 40(3): 141−152. (in Chinese) 谷音,庄舒曼,卓卫东等,2015. 考虑饱和土的地铁车站结构非线性地震反应研究. 岩土力学,36(11):3243−3251.Gu Y., Zhuang S. M., Zhuo W. D., et al., 2015. Analysis of nonlinear seismic response of subway station considering saturated soil. Rock and Soil Mechanics, 36(11): 3243−3251. (in Chinese) 李立云,2007. (准)饱和土与地下结构非线性动力相互作用问题研究. 北京:北京工业大学. 刘光磊,宋二祥,刘华北,2007. 可液化地层中地铁隧道地震响应数值模拟及其试验验证. 岩土工程学报,29(12):1815−1822. doi: 10.3321/j.issn:1000-4548.2007.12.012Liu G. L., Song E. X., Liu H. B., 2007. Numerical modeling of subway tunnels in liquefiable soil under earthquakes and verification by centrifuge tests. Chinese Journal of Geotechnical Engineering, 29(12): 1815−1822. (in Chinese) doi: 10.3321/j.issn:1000-4548.2007.12.012 刘华北,宋二祥,2005. 可液化土中地铁结构的地震响应. 岩土力学,26(3):381−386,391.Liu H. B., Song E. X., 2005. Earthquake induced liquefaction response of subway structure in liquefiable soil. Rock and Soil Mechanics, 26(3): 381−386,391. (in Chinese) 宋佳,2017. 饱和土场地波动数值模拟方法及其工程应用. 北京:北京工业大学.Song J., 2017. Wave numerical method of saturated site soil and its engineering application. Beijing:Beijing University of Technology. (in Chinese) 王相宝,2014. 基于流固耦合两相介质动力模型的饱和土体-地下结构体系近场波动问题研究. 北京:北京工业大学. 王子辉,2008. 饱和两相与单相土互层场地中地铁车站地震反应分析. 北京:北京交通大学.Wang Z. H., 2008. Seismic analysis of subway station in a site interbedded by saturated two-phase and single-phase soil. Beijing:Beijing Jiaotong University. (in Chinese) 许民泽,崔春义,李静波等,2021. 饱和砂土场地中地铁车站结构地震易损性分析. 工程力学,38(S1):251−258.Xu M. Z., Cui C. Y., Li J. B., et al., 2021. Seismic vulnerability analysis of subway station embedded in saturated sand layers. Engineering Mechanics, 38(S1): 251−258. (in Chinese) 杨军,宋二祥,陈肇元,2003. 饱和土一维简谐响应解析解的求解和应用:Ⅱ应用. 岩土力学,24(5):710−714. doi: 10.3969/j.issn.1000-7598.2003.05.008Yang J., Song E. X., Chen Z. Y., 2003. Application of analytical solution of 1-D harmonic response in saturated soil. Rock and Soil Mechanics, 24(5): 710−714. (in Chinese) doi: 10.3969/j.issn.1000-7598.2003.05.008 禹海涛,王治坤,刘中宪,2022. SV波入射下均匀饱和地层渗透系数对深埋隧道的影响机制. 岩土工程学报,44(2):201−211.Yu H. T., Wang Z. K., Liu Z. X., 2022. Influence mechanism of permeability coefficient in homogeneously saturated strata on responses of deep tunnels under incidence of SV waves. Chinese Journal of Geotechnical Engineering, 44(2): 201−211. (in Chinese) 袁宗浩,蔡袁强,曾晨,2015. 地铁列车荷载作用下轨道系统及饱和土体动力响应分析. 岩石力学与工程学报,34(7):1470−1479.Yuan Z. H., Cai Y. Q., Zeng C., 2015. Dynamic response of track system and underground railway tunnel in saturated soil subjected to moving train loads. Chinese Journal of Rock Mechanics and Engineering, 34(7): 1470−1479. (in Chinese) 周巧芳,2008. 饱和土中地铁车站的地震响应分析. 北京:北京交通大学. Biot M. A., 1956. Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Low‐frequency range. The Journal of the Acoustical Society of America, 28(2): 168−178. doi: 10.1121/1.1908239 Li L., Jiao H. Y., Du X. L., et al., 2020. Fully fluid-solid coupling dynamic model for seismic response of underground structures in saturated soils. Earthquake Engineering and Engineering Vibration, 19(2): 257−268. doi: 10.1007/s11803-020-0560-3 Yang Z. H. , Lu J. C. , Elgamal A. , 2008. OpenSees soil models and solid-fluid fully coupled elements user’s manual. San Diego: University of California. Zhu J., Li X. J., Liang J. W., 2020. 3D seismic responses of a long lined tunnel in layered poro-viscoelastic half-space by a hybrid FE-BE method. Engineering Analysis with Boundary Elements, 114: 94−113. doi: 10.1016/j.enganabound.2020.02.007 Zienkiewicz O. C., 1982. Basic formulation of static and dynamic behaviours of soil and other porous media. Applied Mathematics and Mechanics, 3(4): 457−468. doi: 10.1007/BF01908222 -