1．引言

场地条件对地震地面运动有很大的影响，不同的场地条件对地震波的传播能力不同，其对地震动响应的程度也不同。自从1926年，Wood通过对发生在1906年的美国旧金山大地震的震害分布资料的分析研究，首先认识到了场地条件对地震作用下地面运动的影响(Richter,1950)[1]，岩土地震工程领域的学者们针对场地条件对地震作用下的地面运动及地震反应的影响都进行了大量研究，并逐步认识到了场地条件在工程抗震中的重要性。

一直以来，场地条件对地震地面运动的影响的研究有三种方法，一是间接的近似估计方法，二是经验地震动衰减关系方法，三是直接的理论分析方法[2]。其中，理论分析方法运用的更为广泛，而理论分析中的诸多方法主要分为两大类:确定性方法与不确定性方法。目前，土层地震反应分析方法的研究和工程应用是以确定性分析方法为主，而其中应用最广泛的是频域等效线性化方法与时域非线性方法。

等效线性化方法是一种近似处理土体非线性的办法，使用等效线性方法分析场地响应是一个迭代过程，其中为每个土层提供剪切模量和阻尼的初始估计。频域等效线性法在运动方程的求解方面可以认为是精确的，但是由于在处理土的非线性特性时，采用了等效线性化的方法，不能真实的反应土体的非线性特性,使得等效线性方法在处理地震动较强、土层较软的场地时,其计算结果不合理。为了解决这一问题，则需要进行土层地震非线性反应分析,其核心是发展实用的土体动力本构模型。

2．一维时域非线性方法中土的本构模型发展

非线性方法以土体动力本构模型为核心,力求真实地反映土体在地震等动力荷载作用下的动力特性。随着试验技术的不断发展和土体非线性实验资料积累的日益增加,描述土体应力应变滞回特性的动力本构模型也日益增加,根据建立的形式将其分为两大类[3]:一类是经验函数拟合模型,一类是物理元件组合模型。

2.1．经验函数拟合模型

经验函数拟合模型指的是在经验数据的基础上直接建立应力-应变的经验函数关系，它主要是依据骨架曲线和Masing准则的基础而发展的[3]。在给出初始加载条件下的动应力应变式，再利用Masing的“二倍”法得出卸载再加载条件下的动应力应变关系表达式,从而构成滞回曲线方程。

Masing准则最早由Masing[4]在1926年提出，它主要适用于等幅循环荷载，对于非等幅循环荷载情况，可能出现土体应力超过土体的极限应力的不合理的现象，为此一些学者对Masing准则提出了一些补充[5-6]，其中“上骨架曲线”规则和“上大圈”规则得到广泛的应用。“上骨架曲线”规则是如果卸载和反向加载的应力应变曲线与骨架曲线相交,则后续的加载应力应变曲线将遵循骨架曲线。“上大圈”规则即是如果卸载和反向加载应力应变曲线和先前加载曲线相交,则当前应力应变曲线将遵循先前的应力应变曲线[7-9]。根据Masing准则和补充规则,可以获得任意荷载过程的土体非线性应力应变对应关系,但是当应用于复杂荷载时,这样的做法不便于实现,为此Pyke[9]提出了Masing“n倍”法来对Masing准则进行修正,他提出卸载和再加载曲线以前一次卸载开始点为起点,以土的极限应力水平线为渐进线构造,计算时只需记录前次卸载开始点的坐标,减少了大量记录并且使得土应力不超过土的极限应力。但是由于Pyke方法中采用“n倍关系”使得在等幅循环荷载作用下得到的应力应变曲线不具有对称性，且滞回圈随着加载次数的增加而减小可导致加卸载滞回曲线的顶点不落在骨架曲线上，虽然这与土体强度随加载循环次数的增加而减小的实验结果具有相同趋势，但这是其模型考虑其他因素时带入的，并非由试验结果定量确定，所以不能认为是控制土体强度退化效应的有效措施。Matasovic和Vucetic等人[10-11]在1993年根据饱和砂土往返载荷试验结果，给出了土的初始滞回圈和任意后续滞回圈之间的关系，即假设从第2周期后续滞回圈用衰退骨架曲线和Masing法则来描述，则土的往返衰退特性可以对初始骨架曲线的纵坐标加以折减得到后续骨架曲线的纵坐标的方式来表达，开发出了改进后的压力相关的双曲模型。Hashash和Park[12]在2001年在Matasovic的双曲模型模型中引入了两个额外的参数Beta(β)和s来调整骨架的形状，提出了改进的双曲模型。国内方面，王志良[13]也对Masing准则进行了修正,通过引入了“阻尼比退化系数”使得Masing准则所给出的应力应变关系能符合实验阻尼比结果。李小军[14]提出了“修正骨架曲线”的概念对Masing准则进行修正,使得初始加载过程及应变值超过应变历史最大应变值的加载过程,应力应变曲线遵循骨架曲线,卸载及反向加载应力应变曲线直接指向应力应变历史的最大值点或其反向对称点。并且，李小军[15]在Pyke本构模型的基础上，结合Masing准则，并引入了王志良提出的阻尼比退化系数概念，提出了“动态骨架曲线”的概念，导出了能拟合实验阻尼比的形式简单、应用方便的剪切应力-应变关系的经验函数式。栾茂田[16]把Ramberg-Osgood模型及相应的Masing准则进行了修改,把原模型的定常参数取为剪应变水平的函数,从而得到与实验结果相一致的变参数的Ramberg-Osgood本构模型。栾茂田[17-18]又将定参数的Ramberg-Osgood或双曲线模型等解析经验模型及相应的加卸载准则进行了改进,提出了一种简便可行的半经验半公式化方法,所构造的动态非线性模型能逼近土的非线性与滞回特性,并将这一模型推广至不规则加荷条件。金星，孔戈等[19]将Newmark平均加速度法结合中心差分法，推导出了两者结合的增量形式，之后在透射边界条件下结合Pyke本构模型，提出了一种新的显性有限元方法。陈学良[20]在2006年时，在“阻尼比退化系数”模型中又引入了一个调整参数,从而得到了“广义阻尼比退化系数”模型,该模型用可变化的曲线形状来拟合试验阻尼效应,有很大的灵活性。齐文浩、薄景山[21-22]在Pyke模型和Masing准则的基础上，选用指数函数作为骨架曲线，构造了一个适用于非对称循环荷载的指数形式的本构模型，称之为UE模型，除了Pyke模型的优点外，它更加的简单，便于数值实现。为了更好地拟合各类土的动剪切模量比曲线，Martin[23-24]等人在1978年提出利用3个参数拟合G/Gmax-γ曲线的Davidenkov模型，但是这种做法又存在一个问题,即当剪应变幅值无穷增长时,剪应力也无穷增大，为此陈国兴等人[25]利用上限应变幅值作为分界点,采用分段函数法对骨架曲线进行了修正，解决了这一问题,并推导Davidenkov模型的阻尼比的计算公式。张克绪等[26]基于非Masing准则,在所提出模型中引进了两个模型参数，改进了土弹塑性模型的构成方法,使弹塑性模型更具有一般性。Yee等人[27]在2013年提出使用复合双曲线来控制剪切强度的GQ/H模型，该模型允许定义失效时的剪切强度，同时仍然提供表示小应变土壤行为的灵活性，模型中引入的阻尼减小因子来匹配实验室测量的阻尼曲线，卸载-重载刚度使用非Masing标准。

2.2．物理原件组合模型

物理元件组合模型是用一系列弹性元件、粘性元件和塑性元件串联和并联的不同组合来模拟土体的非线性应力应变关系,这种模型以Iwan模型[28]为代表。早在1930年时就提出过一种假设，一般滞后系统可能是由大量的具有不同屈服水平的理想弹塑性元件构成的，但是在当时这一假设没引起足够的重视,直到1967年Iwan阐述了该类物理模型特性，并将Masing模型的概念引申到空间问题中,才引起人们的注意，故称为Iwan模型[28]。Iwan模型的总应力应变关系以弹塑性元件的应力应变关系为依据，在计算土层反应时可能出现高频振荡以及难以拟合阻尼比随应变变化的问题,一些学者试图利用附加粘性阻尼元件的方法来解决这些问题。例如，李小军[14]以Iwan模型为基础提出一种粘-弹-塑物理模型，该模型不但能很好的拟合试验曲线,而且能很好的满足与加载频率无关的要求。郑大同[29]在Iwan模型的基础上,提出了一个新的物理模型,该模型认为骨架曲线与滞回曲线的两个分支即可以相同,也可不同，而Masing模型只是其中的一种特殊情况。蔡袁强[30]通过动三轴试验得到软黏土的软化指数经验模型,并将软化模型引入了Iwan串联模型中，且考虑循环软化的影响，同时再串联一个理想刚塑性元件，模拟循环过程中产生的累积塑性应变，最后结合Masing准则对软黏土的动应力-应变关系进描述。

3．一维时域非线性土层地震反应程序

在地震工程领域中，一维土层地震反应分析模型是最早被提出的土层反应分析模型，多年的发展使其被国内外工程界广泛应用，其满足了各类常规工程场地地震动分析最基本的要求，也是继续研究二维与三维土层地震反应分析方法的基础。而伴随着计算机技术的急速发展，基于数值计算的数值研究在岩土地震工程学的各个领域起着越来越重要的作用，很多数值分析方法被用来研究场地条件对地震动的影响，发展了各种各样的场地地震反应分析程序，并广泛应用于国际岩土地震工程领域。20世纪70年代以来，国内外学者开始了地震反应分析技术的研究，并取得了一系列研究成果。

3.1.CHARSOIL

Streeter在1974提出的CHARSOIL程序[31]是国际上首个时域非线性土层地震反应分析程序。该程序是基于有限差分法，使用Ramberg-Osgood本构模型来反映土体的动力非线性特征，对土层地震反应过程中土体真实的运动情况做出了模拟，开创了在时域中进行土层地震反应分析的先河。

CHARSOIL程序中使用的边界条件是是刚性边界，刚性边界条件是指边界两侧介质中的有着数值完全相同物理参数，即加速度、速度、应力等参数的数值大小一致。从基岩处输入一条入射波，当该入射波传递至地表的自由边界时，将产生一条反射波，这条反射波与基底输入的入射波相比，有着相同的振动方向，但应力方向相反。之后，这条反射波竖直向下继续传播，当到达土层与基岩的交界处时，由于刚性边界条件的原因，位移变为零，而应力变为原来的二倍，因此产生一条与基岩入射波位移相反，传播方向相同的反射波。该反射波会继续向地表处传播，到达地表处时，又会产生一条新的反射波，如此循环往复，因此导致了在使用原始的CHARSOIL程序进行计算时，会产生与实际不一致的情况[32]。而在Ramberg-Osgood模型中，剪应力与剪应变的关系是单调递增的，随着剪应变的增大而增大，这与实际的情况明显不相符，而且该本构模型不能很好的考虑实际土体的阻尼效应[32]。由于CHARSOIL采用刚性边界条件，若计算模型的边界不为刚性边界时，CHARSOIL计算的得到的结果将会与实际情况不相符，存在较大的差异；而且Ramberg-Osgood本构模型不能准确的拟合土体实际的阻尼比试验曲线，为了较好的拟合试验曲线，还需对CHARSOIL程序中所使用的本构模型进行改进。

目前在国内，刘红帅[33]等人对CHARSOIL程序所用Ramberg-Osgood模型进行改进，修改为李小军提出的“动态骨架曲线模型”,并对修改后的软件可靠性进行了分析，发现改进后的软件计算结果较之前的结果有了相对的提高，但是与实际情况以及DEEPSOIL等新型程序计算结果相比还是有一些差距。

3.2.DESRA-2C

DESRA-2C[24]是由Lee和Finn在1978年开发的一个时域非线性土层地震反应分析程序。该程序是基于有限元法，使用了双曲线应力-应变关系（Duncan和Chang1970）作为土壤本构模型来反映土体的动力非线性特征，对土层地震反应过程中土体真实的运动情况做出了模拟。

DESRA-2C解决了时域中的一维非线性剪切波方程，并结合了粘滞/滞回阻尼[35]。该方法假设一个集总质量土壤系统，并将波动方程应用于集总质量方程组，通过假设土壤根据初始加载的双曲线模型和卸载/重新加载的Masing标准来考虑土壤非线性。同时，DESRA-2C允许进行全部和有效的应力场地响应分析，它还是基于MFS孔隙水压力模型的第一个非线性动态有效应力分析,可以进行孔隙水压力的重新分配和消散。而由于其使用的邓肯-张本构模型，它是以虎克定律为基础，假定模型中的参数是应力状态的函数，与应力路径无关，利用土体常规三轴试验得到的应力-应变曲线建立了模型参数关系，虽然它能较好的反应土体的主要变形特征，考虑了土体非线性变形中的加载模量和卸载模量的不同，但是它不能反映土体的剪胀性及主应力对模量的影响。所以在使用过程中不能模拟出很符合实际的运动情况，需要在后期随着本构模型的发展，对程序所使用的模型进行改进[36]。

3.3.D-MOD

D-MOD[37]是Matasovic在1995年时开发出来的一个时域非线性土层地震反应分析程序，它是DESRA-2软件的后代，它使用了与DESRA-2相同的动态响应模型，即由Lee开发的时域动力学方程的数值解[38]，但是，与DESRA-2不同的是，在D-MOD中的非线性滞回本构模型是MKZ模型。MKZ本构模型[38]是Matasovic和Vucetic创立的，是在KondnerandZelasko(1963)创建的双曲线模型(KZ模型)的改进。

D-MOD2000[39]是D-MOD的最新版本，它被广泛应用于动力非线性的场地的有效应力分析中。在D-MOD2000程序中，通过一个双曲线的本构关系来模拟土层的响应，模量比衰减曲线和阻尼比曲线可以通过计算非线性应力应变关系的MKZ模型得到。同时，DMOD2000程序使用了一维的集中质量系统，这个系统连接了一个剪切弹簧和阻尼器用以模拟成层土层的反应，然后通过Newmark-β求解运动方程，进行完全非线性分析。此外，虽然完全Rayleigh阻尼表达式被用于程序计算中，但是与刚度成正比的阻尼方法也仍然被应用。同时，D-MOD2000所使用的孔隙水压力耗散和再发布模型是DESRA-2C中所使用的模型。

3.4.SUMDES

SUMDES是明海燕、李相崧[40]在1992年开发出来的一种非线性完全耦合有限元地震分析程序，它是按照Biot和Zienkiewicz所描述的完全耦合分析框架建立的，它的基本假设为：场地水平，并且在水平面无限延伸；地震波从底面垂直向上传播；场地表面各种应力为零；场地底面不透水；地下水位以下的土完全饱和；这些假设和饱和土水平自由场的情况基本吻合。虽然程序假设变量在水平方向是不变的，但它能描述土和孔隙流体在空间各方向的运动，因而SUMDES是个一维多方向的水平场地地震响应分析程序，它能模拟压缩波和剪切波同时传播的情况，并计算地震引起的土层水平运动和垂向沉降以及孔隙水压力的累积和消散[41]。只有当选择合适的土体本构模型，振动引起的侧向土压力变化，土的压缩和膨胀，接近液化时土的超低强度等特性也能得到很好的模拟。

而随着土体本构模型的发展，SUMEDS软件也在不断的发展中，SUMEDS2006[42]是深圳大学苏栋在SUMDES的基础上开发基于windows的应用程序。它以SUMDES原程序为基础，提供了用户友好型的图形输入输出界面，并且在目前为止，它已经嵌入了多种本构模型，包括线弹性模型，与平均正应力相关的线弹性模型，以及考虑与状态相关的剪胀性的临界状态砂土模型[43]。

3.5.DEEPSOIL

DEEPSOIL程序[44]是Hashash,DuheePark和Chi-ChinTsai在伊利诺大学UrbanaChampaign分校开发的，它源于D-MOD程序，可以同时进行一维的非线性分析和一维等效线性分析，具有友好的图形用户界面，并且也考虑到孔隙水压的影响[44]。DEEPSOIL运用总应力分析方法，采用修正的双曲线本构模型和MKZ模型，同时程序还引入了两个附加参数β和s来调整骨架曲线的形状。

DEEPSOIL是目前国外使用广泛的一维时域非线性土层地震反应分析程序，它的显著优势在于：（1）软件内采用了多种数值分析方法（频域等效线性和时域非线性方法）和广泛的土体本构模型，可以适应多种需求；（2）界面交互直观高效，用户可以在短时间内建立个各种相对复杂的场地数值模型；（3）目前最新版本的DEEPSOIL软件还引入了Groholski等人在2015年提出的GQ/H模型[45]，通过调整土体的骨干曲线，令大应变阶段土体的剪应力接近剪切强度，实现对强震过程土体剪切模量变化的准确描述，避免了现有大部分模型中基于共振柱试验提出的剪切模量衰减曲线在大应变阶段的不准确性。在目前的使用过程中，DEEPSOIL软件在巨厚场和软土场上的使用还是与实际记录有一定的差距[46-51]，还需要继续对程序使用的分析方法和本构模型进行研究改进。

4．展望

对土体动力本构模型的研究是非常复杂的研究课题，目前，很多动力本构参数及其模型形式等来源于静力本构，如何更合理、更真实的反映动力性态仍需完善。现阶段，动力本构研究主要通过试验实现，而原位动力试验、钻井台阵观测研究还是很有限。因此，提高试验和观测的数量和精度，开展大比例模型振动台试验和离心机试验，加强钻井台阵观测或者原位人工动力实验，进一步检验和修正现有土体非线性本构模型。另外，土体静力固结状态对动力反应影响、以及二维模型向三维模型的推广、模型合理性实验等值得研究。

其次，中国需要加强对自主产权的岩土地震工程软件的研制。岩土工程软件是先进、成熟理论和方法在科研和工程应用中最强有力的工具，是提高岩土地震工程工作效率的必需手段，同时也是实现岩土地震工程智能化、集成化和可视化的保证。国内岩土地震工程软件相对于国外已经落后，但国外软件也并非完美，每一款软件都有其“黑匣子”和不成熟之处，发展“清晰透明”的中国自主产权岩土地震工程软件势在必行。

5．结语

土层地震反应分析方法的研究是目前岩土地震工程研究领域的重要内容，土层地震反应分析技术也是当前工程场地地震安全性评价的核心技术之一。自二十世纪初期以来，人们对土体非线性认识在慢慢的提高，进而提出了越来越多的土体本构模型，使得基于数值分析的一维土层地震反应分析软件也得到了不同的发展。但是对目前使用的最新版时域非线性方法计算分析程序我们还需要进行适用性的分析，对于巨厚场或软土场地等复杂的场地的实际计算结果还需要进行大量的计算比较，并且需要继续开展非线性土层地震反应分析方法研究和随机理论适用于真非线性和复杂场地的研究，这样可以一定程度上促进土层地震反应分析软件的改进。同时，随着目前OpenSees和ABAQUS等有限元软件的发展，可以将其作为平台用于进行基于模型的数值模拟,虽然现在有部分学者开始进行了相关的研究，但是得到的成果还不是很成熟，还需要通过大量的实际地震动数据进行计算、验证、改进。

参考文献：

[3]齐文浩,薄景山,张忠利.土层地震反应分析的研究现状[J].世界工程,2010,26(S1):368-372.

[5]李小军,彭青.不同类别场地地震动参数的计算分析[J].地震工程与工程振动,2001(01):29-36.

[6]李小军,廖振鹏.土应力应变关系的粘-弹-塑模型[J].地震工程与工程振动,1989(03):65-72.

[7]张克绪,谢君斐.土动力学[M].北京:地震出版社,1989.

[8]齐文浩.土层地震反应分析方法的比较研究[D].中国地震局工程力学研究所,2004.

[13]王志良,韩清宇.粘弹塑性土层地震反应的波动分析法[J].地震工程与工程振动,1981(01):117-137.

[14]李小军.土的动力本构关系的一种简单函数表达式[J].岩土工程学报,1992(05):90-94.

[15]李小军,廖振鹏,张克绪.考虑阻尼拟合的动态骨架曲线函数式[J].地震工程与工程振动,1994(01):30-35.

[16]栾茂田.土动力非线性分析中的变参数Ramberg-Osgood本构模型[J].地震工程与工程振动,1992(02):69-78.

[17]栾茂田,林皋.土料非线性滞回本构模型的半解析半离散构造方法[J].大连理工大学学报,1992(06):694-701.

[18]栾茂田,林皋.场地地震反应非线性分析的有效时域算法[J].大连理工大学学报,1994,34(02):228-234.

[19]金星,孔戈,丁海平.水平成层场地地震反应非线性分析[J]地震工程与工程振动,2004,24(3):38-43.

[20]陈学良.土体动力特性、复杂场地非线性地震反应及其方法研究[D].中国地震局工程力学研究所,2006.

[21]齐文浩,王振清,薄景山.土层非线性地震反应分析方法及其检验[J].哈尔滨工程大学学报,2010,31(04):444-450.

[22]齐文浩.土层非线性地震反应分析方法研究[D].中国地震局工程力学研究所,2008.

[25]陈国兴,庄海洋.基于Davidenkov骨架曲线的土体动力本构关系及其参数研究[J].岩土工程学报,2005(08):860-864.

[26]张克绪,李明宰,王治琨.基于非曼辛准则的土动弹塑性模型[J].地震工程与工程振动,1997(02):74-81.

[29]郑大同,王惠昌.循环荷载作用下土的非线性应力应变模型[J].岩土工程学报,1983(01):65-76.

[30]蔡袁强,柳伟,徐长节,黄洪超.基于修正Iwan模型的软黏土动应力-应变关系研究[J].岩土工程学报,2007(09):1314-1319.

[33]韩蓬勃.一维时域非线性场地地震反应分析程序CHARSOIL的改进[D].防灾科技学院,2019.

[42]苏栋,李相菘.离心机自由场地震响应的完全耦合模拟[J].岩土力学,2006,27(S2):1079-1083.

[43]李靖.密砂的力学性质及密砂场地地震响应研究[D].哈尔滨工业大学,2010.

[47]肖遥.场地地震反应分析方法的研究[D].中国地震局工程力学研究所,2011.

[48]李兆焱,袁晓铭,王鸾,王克.巨厚场地三种土层地震反应分析程序对比检验[J].地震工程与工程振动,2017,37(04):42-50.

[50]杨洋,孙锐,杨洪搏,陈卓识.国际上两种典型土层地震反应分析程序对比研究[J].世界地震工程,2017,33(03):17-23.

[51]张海,王震,周泽辉,尤红兵.基于DEEPSOIL的软土场地地震反应研究[J].震灾防御技术,2015,10(02):291-304.

刘冰冰,乔云龙.一维时域非线性土层场地地震反应的研究现状[C].上海来溪会务服务有限公司.2019第四届土木工程国际会议论文集.上海来溪会务服务有限公司:上海来溪会务服务有限公司,2019:53-58.