0.   引言

作为一种典型复杂场地，沉积盆地对地震动的显著放大作用已在多次震害调查和强震观测中得以证实^[1]。并且由于众多大城市建设于沉积盆地内部或边缘（洛杉矶、大阪、台北和北京等），这些城市人口财富聚集，一旦破坏性地震发生，将产生灾难性后果^[2]。由于现有强震观测台站还不足以覆盖所有城镇区域，观测数据还远不够预测工程中需要的地表反应特征。因此，通过数值模拟方法来模拟震源破裂、地震波的传播、局部场地效应的全过程物理模拟逐渐成为工程地震学和理论地震学的一个热点课题，并基于此来建立复杂场地的强地震动场，可为现有的经验地震动预测方法，特别是对于观测记录稀疏的近断层区域提供了有价值的补充。

针对实际三维沉积盆地区域，许多学者采用不同数值模拟方法：有限差分法^[3]（finite difference method, FDM）、有限元法^[4]（finite element method, FEM）、谱元法^[5]（spectral element method, SEM）等，进行设定地震的全过程地震动模拟。但值得注意的是，受限于复杂场地特性探测数据的匮乏及超大规模数值模拟所需的计算机容量和计算效率，目前相关研究主要集中于弹性地震动场的建立，非线性范围下地震动场的研究则进展有限。然而，实际地形地质结构及软土层的强非线性和夹层的弱非线性特征使得场地的地表地震动具有很强的空间变异性。在中强地震动下，实际观测发现土体的非线性滞回耗能特性对地震动有明显的削弱作用，影响地形和场地的耦合响应^[6]。因此，迫切需要在物理模拟中纳入土体非线性，进一步开展大尺度非线性地震动模拟，给出更加合理的地面运动和地震动场分布。

目前，已有部分学者注意到上述问题，尝试将土体非线性纳入大尺度复杂场地的三维地震波场模拟之中。Xu等^[7]基于FEM，采用Drucker-Prager弹塑性模型研究了土体非线性特性对地震动幅值和主频的影响。Taborda等^[8]基于FEM，对实际三维盆地进行基于完美弹塑性模型的非线性地震动分析。Roten等^[9]在基于FDM的AWP-ODC程序基础上开发了的Drucker-Prager弹塑性模型，模拟了频率高达4 Hz的San Andreas断层下的M7.7设定地震。Dupros等^[10]同样基于FDM的AWP-ODC程序，采用Mohr-Coulomb模型，针对法国Riviera地区开展了考虑土体非线性的大规模并行模拟。Fu等^[11]优化了基于FDM的AWP-ODC和CG-FDM程序，利用TaihuLight超算对唐山地区进行了频率上限到18 Hz的非线性地震动模拟。Esmaeilzadeh等^[12]基于FDM开展了加拿大西北-东南长8 km的Kinburn沉积盆地的0.1~1 Hz的非线性地震动模拟。Chen等^[13]基于SEM，采用一维等效线性方法考虑了近地表土体非线性反应，针对2016年M_w7.0熊本地震开展了地震动模拟。

上述针对大尺度沉积盆地的非线性地震反应分析的研究仍有所欠缺。一维等效线性方法在处理软土大变形方面的缺陷尤为突出，无法实时反映土体参数的演变和模型的三维效应。弹塑性分析方法发展出了不同应变硬化机制和流动法则下的模型，涉及许多没有明确物理意义的中间变量，缺乏针对不同土类可靠的动力反应分析方法，并且在大尺度场地地震动分析中计算效率低。考虑上述问题，本文借助谱元法开源软件SPECFEM3D在模拟实际高度复杂异质体和易于并行化等优势，将Martin-Seed-Davidenkov本构模型开发到程序中，建立了从震源到场地（包含复杂地层、沉积盆地和场地非线性地）的全过程地震动数值模拟方法，实现了完整的全过程强地面运动模拟。上述研究可为大尺度沉积盆地地区提供更可靠具体的工程抗震设计地震动参数。

1.   土体非线性本构在谱元法中的开发

在Masing类非线性动力本构模型中，双曲线类模型由于其概念清晰和运用方便在工程中应用较为广泛^[14-16]。本文即选取Martin等在此基础上提出了三参数的Davidenkov模型来开展在SPECFEM3D中的非线性本构开发工作。本节首先对Davidenkov模型相关理论公式进行介绍，其次给出了SPECFEM3D中非线性本构的具体开发步骤，最后通过与一维非线性动力分析软件DEEPSOIL结果进行对比说明开发程序的正确性。

1.1   Davidenkov模型基本理论

Davidenkov模型的骨架曲线^[17]可表示为

式中：${G_{\max }}$根据共振柱试验或经验公式确定；A，B，${\gamma _0}$根据不同土性的试验剪切模量比曲线拟合确定。对于三维问题，需选取一个广义剪应变以在一维应力-应变空间中进行滞回曲线分析，通常选取八面体剪应变或偏应变的第二不变量来表示，本文采用后者方法将广义剪应变表示为

动应力-应变滞回曲线的推导则需要借助Masing法则，在传统的再加载和卸载曲线的坐标比例尺为初始加载曲线二倍的假设下，再加载和卸载曲线可以表示为

Pyke^[18]以n倍法修正了传统的二倍法，有效控制了滞回曲线的发展方向处于骨架曲线渐近线内。赵丁凤等^[19]进一步对滞回曲线进行了修正，使得加卸载曲线沿最新拐点指向历史最值点方向前行并始终能回到骨架曲线上。

Davidenkov骨架曲线和滞回曲线的剪切模量可分别通过骨架曲线和滞回曲线对剪应变求导(${\text{d}}\tau$/${\text{d}}\gamma$和d($\tau - {\tau _{\text{c}}}$)/d($\gamma - {\gamma _{\text{c}}}$))得到。本文基于上述Davidenkov非线性本构模型在SPECFEM3D中进行二次开发，以实现地震中土体非线性特性的模拟。

1.2   SPECFEM3D开发过程

SPECFEM3D中内置程序为FORTRAN语言编写，本次开发主要任务是根据当前时间步GLL节点的应变增量，结合基于不规则加卸载准则的Davidenkov本构模型，给出相应的应力增量；更新并传递计算过程中的状态变量。图 1为谱元法SPECFEM3D软件中的开发流程，其中主要包含4个主要步骤。

图  1  谱元法中非线性本构开发流程

Figure  1.  Development process of nonlinear constitutive in spectral element method

下载: 全尺寸图片 幻灯片

（1）标定参数和初始化全局变量。根据不同土类的G-$\gamma$曲线拟合骨架曲线中的参数A、B和${\gamma _0}$，同时初始化全局变量，包括上一步应变张量ε^t和应力张量σ^t、上一步拐点处的应变张量ε_{ij, c}、拐点处等效应变${\gamma _{\text{c}}}$和应力${\tau _{\text{c}}}$、上一步等效应变${\gamma _{{\text{eq}}}}$和广义应变${\gamma _{{\text{gen}}}}$、上一步广义应力τ_eq、上一步等效应变和应力极值${\gamma _{{\text{ex}}}}$和${\tau _{{\text{ex}}}}$以及加卸载标识sign=1（1表示加载，-1表示卸载）。

（2）更新各GLL节点处应变并计算等效应变增量。以拐点处的应变张量更新GLL节点t+1时刻的应变$\varepsilon _{ij}^{t + \Delta t} = \varepsilon _{ij}^{t + \Delta t} - {\varepsilon _{ij,c}}$并按式（2）计算等效应变$\gamma _{{\text{gen}}}^{t + \Delta t}(e_{ij}^{t + \Delta t})$，计算广义应变增量$\Delta \gamma _{{\text{incre}}}^{t + \Delta t}(e_{ij}^{}) =$$\gamma _{{\text{gen}}}^{t + \Delta t}(e_{ij}^{t + \Delta t}) - \gamma _{{\text{gen}}}^t(e_{ij}^t)$。更新广义等效应变$\gamma _{{\text{gen}}}^t = \gamma _{{\text{gen}}}^{t + \Delta t}$。

（3）判断应变转向点是最为关键的步骤。①若等效应变增量小于零，则判定为当前时间步出现加卸载拐点。令加卸载标识取为其相反数sign=-sign，记录上一步应变张量ε^t并将最新广义应变清零，即令$\gamma _{{\text{gen}}}^t = 0$。更新拐点处的广义应变和应力值，即$\gamma _{\text{c}}^{t + \Delta t} =$ $\gamma _{{\text{eq}}}^t$和$\tau _{\text{c}}^{t + \Delta t} = \tau _{{\text{eq}}}^t$。②否则，以广义应变增量更新等效应变$\gamma _{{\text{eq}}}^{t + \Delta t} = \gamma _{{\text{eq}}}^t + {\text{sign}}|\Delta \gamma _{{\text{incre}}}^{t + \Delta t}|$。判断其与等效应变的关系，若超过则更新等效应变极值$\gamma _{{\text{ex}}}^{} = \gamma _{{\text{eq}}}^{t + \Delta t}$并按式（1）骨架曲线公式计算等效应力$\tau _{{\text{eq}}}^{t + \Delta t}$；若未超过则根据式（3）计算$\tau _{{\text{eq}}}^{t + \Delta t}$。

（4）根据初始加载或后继加卸载状态计算剪切模量G^t+Δt。根据最新剪切模量和当前传来应变增量确定应力值，即σ^t+Δt=σ^t+λ^t+ΔtΔε_kk+G^t+ΔtΔε_ij。图 1给出了开发的主要流程图。

2.   非线性谱元法正确性检验

通过与一维非线性软件DEEPSOIL结果对比来说明在SPECFEM3D中开发非线性本构的正确性。DEEPSOIL是美国Illinois大学开发的一款场地非线性分析软件，采用Newmark数值方法求解一维系统多自由度动态平衡方程。采用DEEPSOIL中非线性MKZ本构模型，不考虑围压影响，骨架曲线表示为

式中：${G_{{\text{m0}}}}$为初始剪切模量；β，s，${\gamma _{\text{r}}}$为对应特定土体试验数据的拟合参数。黏性阻尼比则等效为正比于滞回圈面积的变量，由此可由滞回曲线求得特定土体的动剪切模量比和阻尼比曲线。DEEPSOIL中对阻尼项的考虑采用扩展Rayleigh阻尼公式，可考虑土体的N阶特征频率，相应的阻尼矩阵可表示为

式中：${a_b}$为与阻尼比和特征频率有关的参数，

本节采用的验证模型及材料参数如图 2所示，拟合得到确定骨架曲线的相关系数也标于图 2（a）中。模型底部入射垂直SV波，波形为图 2（b）所示的Ricker子波形。为将结果退化至一维与DEEPSOIL结果进行对比，在厚8 m的土层下方设置2 m刚度很大近似为刚体的土层，限制模型在y和z方向的运动，且在模型底部设置吸收边界。这里在平面波垂直入射情况下，采用Stacey底部吸收边界即可满足精度要求。传播：

图  2  与DEEPSOIL对比验证的SPECFEM3D三维模型和入射Ricker波形

Figure  2.  SPECFEM3D 3D model verified by comparison with DEEPSOIL (a) and incident Ricker waveform (b)

下载: 全尺寸图片 幻灯片

式中：n为边界上正法向量；s为切向向量；T_total和u_total为SEM求解得到的边界上的总牵引力和位移场值，${\boldsymbol{T}_{{\text{FK}}}}$和${\boldsymbol{u}_{{\text{FK}}}}$为FK方法求解得到的内外域边界上的牵引力和位移场值^[20]。

首先针对滞回曲线进行验证，图 3给出了两种不同应变路径条件下分别根据原始Masing准则的DEEPSOIL和本文基于不同加卸载准则计算得到的土体应力-应变滞回曲线结果。从图 3中可以看出，扩展Masing准则构造滞回曲线需要存储全部拐点信息并通过数值方法寻找后续波与先行波滞回曲线的交点，而本文采用的文献[19]中提出的不规则荷载加卸载准则，加卸载曲线以最近拐点为起点并指向历史最值点，使得编程时仅需记录最新拐点信息，且该规则更符合试验结果。

图  3  给定应变路径下基于Masing准则和本文不规则加卸载准则的滞回曲线对比

Figure  3.  Comparison of hysteresis curves based on Masing criterion and irregular loading and unloading criterion under given strain path

下载: 全尺寸图片 幻灯片

其次对DEEPSOIL和SPECFEM3D求解地震响应结果进行比较，在土层底部入射垂直SV波，输入波形选取中心频率f₀=2.0 Hz的Ricker波，峰值加速度分别为调整为0.2g，0.6g，1.0g，1.4g。为保证加卸载过程中拐点正确识别，SPECFEM3D中时间步长取dt=2.0×10^-5 s，由于非线性计算中剪切模量的变化，结果含有复杂的高频成分，因此采用三阶Butterworth低通滤波器将结果滤波至输入地震动的频带上限（3.27 Hz）。图 4给出了在输入不同强度地震动情况下SPECFEM3D与DEEPSOIL线性和非线性地表加速度结果的对比，可以看出二者的线性计算结果完全吻合，说明三维模型退化至一维施加的边界条件正确，二者非线性结果加速度时程及其峰值基本吻合。

图  4  不同强度地震动下SPECFE3D和DEEPSOIL的地表原点处水平加速度对比

Figure  4.  Comparison of horizontal acceleration of SPECFE3D and DEEPSOIL under different intensity ground motions

下载: 全尺寸图片 幻灯片

3.   施甸盆地非线性地震动模拟

施甸盆地位于中国地震活动性强烈的滇西南地区，位于云南省西部边陲，怒江东岸，保山市南部，地处东经98˚54—99˚21′，北纬24˚16′—25˚00′之间，与昆明市相距约360 km，施甸县城坐落在该盆地的南部。1976年5月29日，在施甸县城西南约50 km的龙陵先后发生的7.3和7.4级地震，施甸盆地均出现了Ⅶ~Ⅷ度的高烈度异常区，相比周围地区的地震烈度高了1~2度。1988年11月6日，在施甸东南150~200 km左右的澜沧－耿马发生7.6和7.2级地震时，施甸县城再次遭到破坏，施甸盆地内出现了Ⅵ度高烈度异常区，而周围地区地震烈度在Ⅴ度左右，无明显震害现象。2001年施甸发生5.9级地震中，施甸盆地的震害也较周围严重。施甸盆地在历史强震中，重复出现高烈度异常的现象绝不是偶然的，与该盆地的地质背景和工程场地条件有关。针对施甸盆地的高烈度异常的原因、产生机理的研究，对于研究沉积盆地对地震动作用机理，以及对施甸地区的防震减灾具有重要意义。本节基于谱元法中开发的非线性本构模型针对施甸盆地进行设定地震地震动模拟，研究沉积盆地土体非线性对地震动产生的影响。

选择研究计算区域如图 5所示（东经98˚54′— 99˚21′，北纬24˚16′—25˚00′），其中黑色框内为沉积盆地位置。施甸盆地为第四纪较为发育的向斜断陷盆地，盆地大致沿南北走向发育，东西宽约为3 km，南北长约21 km。通过国家地理空间数据云获取该区域数字高程信息并结合文献[21]中建立的施甸盆地的三维速度结构模型，利用剖分软件Trelis建立包含起伏地形的三维速度结构模型。为简化计算，盆地介质区域内P波和S波的速度取为1100，500 m/s，盆地外参考川滇地区的一维地壳速度结构模型以及CRUST1.0提供的速度结构，介质参数如表 1中所示。其中，本文采用粗粒度（coarse-grain）方法^[22]根据地壳波速结构计算品质因子${Q_{\text{P}}}$和${Q_{\text{S}}}$，${Q_{\text{S}}}$＝0.1${v_{\text{S}}}$，${Q_{\text{P}}}$＝2${Q_{\text{S}}}$。

图  5  施甸盆地模拟研究区域

Figure  5.  Simulation study area of Shidian Basin

下载: 全尺寸图片 幻灯片

表  1  施甸盆地覆盖层及地壳层材料参数

Table  1.  Wave velocity structure of Shidian Basin and crustal layers

深度/km	${v_{\text{P}}}$/(m·s-1)	${v_{\text{s}}}$/(m·s-1)	密度/(kg·m-3)	${Q_{\text{P}}}$	${Q_{\text{S}}}$
沉积	1100	500	1800	100	50
3	4570	2670	2300	534	267
5	4626	2868	2572	574	287
8	5577	3278	2600	654	327
10	5898	3500	2676	700	350
15	6013	3509	2750	716	358
20	6288	3654	2803	732	366

下载: 导出CSV 

| 显示表格

盆地外部不考虑土体非线性，内部沉积采用本文开发的非线性本构。经查阅相关地质资料，施甸盆地内覆盖层中淤泥质和泥炭质黏土发育。盆地基地地层为泥盆系和石炭系地层，以白云岩、灰岩等碳酸岩类构成，盆地边缘为奥陶系和志留系地层，以页岩、砂岩、白云岩和灰岩为主。采用Darendelim等^[23]提出的G-$\gamma$曲线拟合骨架曲线相关参数，其中，假设土体为正常固结的非塑性土（I_P=0，OCR=1），N_c和F_req取建议值10.0和1.0，K₀通过土体的有效内摩擦角确定K₀=1-sinφ，这里取φ=40°。通过土体塑性指数（I_P）、超固结比（OCR）、静止侧压力系数（K₀）、加载循环数（N_c）和频率（F_req）5个参数可确定G-$\gamma$曲线（图 6），与彭胜恩等^[24]基于岩石试验数据得到的动剪切模型比和阻尼比拟合曲线对比，验证了模拟中场地土选用的动剪切模量比与阻尼比曲线的合理性。为保证模拟的频率能达到2 Hz，盆地内的网格尺寸设置为水平向50 m左右，垂直向20 m左右。为减少计算量，盆地外部区域的网格大小从500~2000 m。计算模型共包含2638125个谱单元。在国家超算中心“天河一号”平台上进行计算。调用400个进程的并行计算系统上完成求解，时间步距取为0.002 s。模拟30 s内的地震波传播，非线性计算时间约4.5 h。

图  6  场地土的动剪切模量比与阻尼比曲线

Figure  6.  Curves of shear modulus ratio and damping ratio of site soils

下载: 全尺寸图片 幻灯片

由于目前没有学者给出2001年施甸5.9级地震的震源反演滑动分布结果，故本文采用能同时考虑破裂面上低波数确定性错动和高波数随机性错动的GP14.3运动学混合震源模型建立此次地震的震源滑动分布^[25]，建立过程如下：①根据云南数字强震台网P波初动资料计算的M_s5.9级地震的震源机制解以及定标律公式^[26]，建立低波数确定性的凹凸体震源模型，具体参数详见表 2，3；②在低波数确定性的凹凸体震源模型中采用波数衰减满足von Karman自相关函数的波数谱并引入随机数ϕ = [-π, π]表达断层破裂面上小尺度的随机变量，引入高波数随机性错动，确保合成宽频地面运动的带宽有效性，构建GP14.3运动学混合震源模型，断层面的时空分布见图 7所示。

表  2  施甸5.9级地震全局震源参数

Table  2.  Global source parameters of Shidian M5.9 earthquake

走向/(°)	倾角/(°)	滑动角/(°)	断层面长度/km	断层面宽度/km	震中深度/km	破裂方式
161	65	-150	8.1	7.9	7.5	中心破裂

下载: 导出CSV 

| 显示表格

表  3  施甸5.9级地震局部震源参数

Table  3.  Local source parameters of Shidian M5.9 earthquake

局部参数	单位	定标律	参数值
面积Sm	km2	lgSm=lgS–0.80	9.10
平均错动量Dm	cm	lgDm=lgD+0.39	46.66
长度Lm	km	lgLm=lgL–0.45	2.68
宽度Wm	km	Wm=Sm/Lm	3.40
沿走向中心Xm	km	lgXm=lgL–0.31	3.79
沿倾向中心Ym	km	lgYm=lgW–0.35	3.79

下载: 导出CSV 

| 显示表格

图  7  2001年施甸5.9级地震震源时空分布

Figure  7.  Spatial-temporal distribution of 2001 Shidian 5.9 earthquake

下载: 全尺寸图片 幻灯片

图 8分别对比了本文计算的2个基岩台站R1（震中距13.52 km）和R2（震中距21.56km）的三分量加速度时程和强震记录，黑线代表强震记录，红线代表本文模拟结果，每条加速度时程的时长均为20 s，PGA均标注于曲线的右上方，单位为cm/s²。强震记录及本文计算结果均经4阶Butterworth滤波至2 Hz。对比模拟结果和观测记录可知：模拟得到的地震动三分量时程结果与观测记录的幅值和波形特征基本上一致，模拟结果反应谱的各周期成分的PGA与强震记录基本吻合。根据上述地震动时程和反应谱的比较分析，验证了本文提出方法的可靠性，总体上说明本文采用的SEM和混合震源模型适用于2001年施甸5.9级地震动模拟。

图  8  地震动时程、反应谱模拟结果与观测记录对比

Figure  8.  Comparison among simulated ground motion time history, Fourier amplitude spectra and ground motion records

下载: 全尺寸图片 幻灯片

图 9给出了施甸盆地弹性和非线性模拟的EW方向速度地震波场快照。图 9中红色表示速度为正值、蓝色表示速度为负值，黑色框位置为施甸盆地位置，五角星位置代表震源位置。波场快照图均从第2.2 s开始，并每隔0.4 s输出相应的波场快照结果。由图 9中可看出：①t=2.2 s时震源激发的地震波已经传播到地表位置处；②当t=2.6 s时震源激发的地震波传播至盆地边缘，此时可见由于震中受局部地形的影响，地震波大致呈近似的圆形向外传播，此时地震波已经引起了沉积盆地内部较小的地震响应，但是还未见到有明显的面波出现；③当t=3.0 s时可以看出地震波已经进入沉积盆地内部，盆地内波前的到时滞后于盆地外，波前发生扭曲，并且由于盆地内外较高阻抗比，可以很明显的看出弹性介质下沉积盆地内部的地震响应明显高于沉积盆地周边地区，但观察同时刻非线性计算的沉积盆地的地震响应小于弹性计算的地震响应；④当t=3.4，3.6 s时，发现地震波的能量大部分被捕捉在弹性沉积盆地之中，形成了振幅较大的面波，面波不断在沉积内部反射叠加，导致盆地内部地震动明显放大。

图  9  地表EW水平方向速度波场快照（每时刻第一和第二行分别对应线性和非线性结果）

Figure  9.  Snapshot of surface EW velocity wave field (first and second lines correspond to linear and nonlinear results respectively at each moment)

下载: 全尺寸图片 幻灯片

图 10，11分别给出了弹性和非线性模拟结果的水平加速度峰值（PGA）分布图和速度峰值(PGV) 分布图，水平加速度和速度取为WE和NS分量的算术平方根。结果显示弹性介质和非线性介质在施甸盆地出现明显差异：①针对弹性介质的施甸盆地，水平峰值速度达到50 cm/s，水平峰值加速度达到0.3g左右，施甸盆地出现地震动的聚集放大效应在PGA和PGV分布图上都有比较清晰的显示，该情况进一步表明了由于沉积盆地与基岩的阻抗比较高，导致地震反射波强烈，地震动的放大效应明显。同时由于是三维沉积盆地，地下构造的不均匀性进一步加强了地震反应，这一规律与所认知的沉积盆地的聚集放大效应规律保持一致。②当考虑土体非线性时施甸盆地水平峰值速度接近35 cm/s，水平峰值加速度达到0.23g左右，较线性计算的幅值降低约30%，并且由于从地震波频率成分而言，PGA一般与高频成分相关，而PGV则与地震波中低频成分相关，由图 11中峰值分布图中可以看出，非线性土体对沉积盆地PGV的影响更加明显。考虑非线性结果对速度结果影响更为明显，为进一步对比沉积盆地不同位置处土体非线性对地震动时程的影响，图 12，13分别给出了近断层地震动下盆地内地表观测点1~6的三分量速度时程和速度反应谱。

图  10  模拟区域以及盆地范围内的地表PGA结果（左图为线性结果，右图为非线性结果）

Figure  10.  PGA results in simulated area and basin (left figure is a linear result, and right figure is a nonlinear result)

下载: 全尺寸图片 幻灯片

图  11  模拟区域以及盆地范围内的地表PGV结果（左图为线性结果，右图为非线性结果）

Figure  11.  PGV results in simulated area and basin (left figure is a linear result, and right figure is a nonlinear result)

下载: 全尺寸图片 幻灯片

图  12  盆地内地表点^#1~^#6处速度三分量时程（图 5中对应盆地内地表观测点^#1~^#6）

Figure  12.  Three-component time histories of velocity at surface points ^#1~^#6 in basin (Figure 5 corresponds to surface receivers ^#1~^#6 in basin)

下载: 全尺寸图片 幻灯片

图  13  盆地内地表点^#1~^#6处速度反应谱（图 5中对应盆地内地表观测点^#1~^#6）

Figure  13.  Response spectral velocities of surface points ^#1~^#6 in basin (Figure 5 corresponds to surface receivers ^#1~^#6 in the basin)

下载: 全尺寸图片 幻灯片

由图 12所示，考虑土体非线性时由结果可看出沉积非线性整体上显著削弱了水平速度幅值且峰值时刻略微时间延后，竖直速度分量幅值有所增加。此外，不同分量上沉积非线性影响存在差异：对于NS方向分量，在^#1，^#2，^#3，^#4，^#5和^#6点处线性模拟速度峰值约为其非线性模拟结果的1.6倍、1.7倍、1.2倍、1.2倍、1.3倍和1.7倍；而对于EW分量，在^#1，^#2，^#3，^#4，^#5和^#6点处线性模拟速度峰值约为其非线性模拟结果的90%、1.2倍、1.7倍、1.2倍、1.5倍和1.7倍，整体上由于虚拟观测点^#1和^#6位于沉积盆地边缘，所处复杂地形，强地震动导致沉积剪切模量的显著衰减，非线性模拟的地震动幅值削减最为显著，盆地中心位置处非线性土体对地震动的影响相对较小。对于竖向UD分量，尽管速度峰值都较大，但非线性结果与线性结果差异很小，这可能是由于沉积盆地内部覆盖层厚度较小，由此产生的非线性动力特性在深度方向变化迟缓导致的。

由图 13所示的速度反应谱可看出考虑盆地土体非线性整体上较线性结果降低了水平方向的反应谱峰值，竖向的反应谱峰值结果差异同样不明显。①接收点^#1，^#2，^#3，^#4，^#5和^#6点处考虑非线性的水平速度反应谱幅值较线性均有所降低：NS分量分别降低了30.25%，18.67%，18.18%，8.73%，29.16%和37.75%，EW分量分别降低了5.63%，9.76%，31.76%，18.74%，25.71%和37.63%。对于竖向的速度反应峰值，其整体反应谱曲线较低且非线性与线性的差异不大，原因同对速度时程的分析。并且从图中观察，除观测点^#6之外，由于土体的刚度弱化，考虑土体非线性后的峰值周期有向长周期移动的趋势。在点^#6对应盆地底部地形较复杂处，考虑土体非线性计算的速度分量约为线性结果的一半，而二者PSV峰值对应特征频率差别不大。上述结果表明，近断层地震动导致沉积介质剪切模量的显著衰减，而近地表土体非线性响应依赖于各分量上的应变发展，由此造成了线性结果与非线性结果的差异，以及各分量上地震反应的差异。场地中介质非线性特性严重影响对风险水平的判断，需特别重视。

4.   结论

本文基于谱元法SPECFEM3D软件完成了土体非线性本构的开发工作，并针对云南施甸盆地开展了设定地震下非线性地震动响应模拟，得到以下3点结论。

（1）基于Davidenkov非线性本构模型在谱元程序SPECFEM3D中进行二次开发，通过与一维非线性软件DEEPSOIL结果对比证明SPECFEM3D中开发非线性本构的正确性，可进一步应用于三维大尺度复杂场地非线性地震动模拟中。

（2）由于沉积盆地与基岩的阻抗比较高，导致地震反射波强烈，地震动的放大效应明显，考虑沉积非线性的水平PGA和PGV较线性结果均降低，并且由于从地震波频率成分而言，PGA一般与高频成分相关，而PGV则与地震波中低频成分相关，非线性土体对沉积盆地PGV的影响更加明显，较线性计算的幅值降低约30%。

（3）考虑沉积非线性使沉积内部观测点整体上速度及反应谱幅值降低，且沉积非线性整体上使特征频率向长周期方向移动，部分地形较复杂尤其是靠近盆地边缘区观测点处竖向模拟峰值结果出现放大情况。