引言

土的动剪切模量比和阻尼比是土动力学特征的重要参数（Hardin等，1968），在工程场地土层地震反应分析和地震安全性评价工作中不可缺少（郭婷婷等，2016）。通过描述土在动力荷载作用下的动应力应变关系（动本构关系）模型，将动剪切模量和阻尼比与动剪应变幅值的函数关系具体化。Seed等（1970）首先给出了砂土和黏性土的动剪切模量和阻尼比与动剪应变幅值的关系曲线。Hardin等（1972）给出了动剪切模量与动剪应变的关系（即Hardin-Drnevich双曲线模型）。Martin等（1982）在Hardin-Drnevich模型的基础上进行了改进，提出了三参数的Davidenkov模型和具有幂次形式的阻尼比拟合公式。陈国兴等（2005）基于Davidenkov模型，采用上限剪应变幅值作为分界点，推导出了修正的Davidenkov模型阻尼比计算公式。在所有模型中，Hardin-Drnevich模型需要求解的参数最少，只需要2个参数，因而较为简便，目前应用较为广泛。下文中讨论动剪切模量和阻尼比与动剪应变幅值关系就基于此模型。

土的动剪切模量和阻尼特性测定是地震等动荷载作用下土工建筑物设计和计算的基本依据（赵红芬等，2009）。通常，为获得土体的动剪切模量和阻尼比，最直接的方法是进行土动力学试验。在土的动力学特性与岩土地震工程试验研究中最常用的方法主要是动三轴和自振柱试验两种方法。其中，动三轴试验主要用于在室内测定粘性土的动强度、各类土的动模量和阻尼比、砂土抗液化强度等相关参数（吴世明等，2001；张涛，2004；申权等，2013；郭可骍，2015）。基于前人对土动力特性的分析，同时考虑到试验测试及数据采集系统的特征，试验荷载一般采用正弦波形循环荷载（南京水利科学研究院，2003）。在动三轴试验中，数据处理是保证试验成果有效产出的关键，而应力应变滞回曲线椭圆拟合则是判断试验成果准确性的决定因素。

1.   循环应力动三轴试验

在循环应力动三轴试验中，土试样直径选取39.1mm，高取为80mm，对试样施加相等的固结压应力σ_1c=σ_3c=σ₀，此时固结应力比K_c，为等压固结（Hardin等，1972；张茹等，2006）。动力特性试验荷载激振频率定为1Hz，逐级施加动荷，每级振次为10次，在排水条件下进行振动试验。在动三轴试验中，测定动模量和阻尼比的方法是将轴向动应力作用于试样，应力幅值由小逐级增大进行试验，当试样变形超过允许值时，试验终止。将试验过程中采集到的各级荷载下的动应力应变滞回曲线进行整理，用以确定不同动应变时的动模量和阻尼比，如图 1所示。动弹性模量E_d为椭圆长轴斜率，阻尼比λ_d定义为循环消耗的能量与施加的最大应变能之比的1/4π倍（李松林，1990），根据应变能定义，一个循环内损耗的能量为滞回圈椭圆的面积πab，其中a、b为椭圆的长短轴；施加的最大应变能为P_dε_d/2，即图中三角形AOB的面积，P_d和ε_d为椭圆与长轴交点中坐标值较大的一点坐标。

图 1  确定动弹性模量和阻尼比示意图

Figure 1.  Schematic diagram of calculating dynamic elastic modulus and damping ratio

下载: 全尺寸图片 幻灯片

2.   拟合滞回曲线椭圆的方法分析

试验采用正弦波循环荷载，即动应力可以表示为：

式中，P_a表示动应力，P_max为动应力最大变化幅度，P₀为静态基准应力，ω为荷载振动的角频率，t为时间，${{\varphi }_{1}}$为动应力的初相。由于存在阻尼，动应变与动应力不再满足简单的线性关系，而是存在一个相位差。

对应的动应变可以表示为：

式中，ε_a为动应变，ε_max为动应变最大变化幅度，ε₀为静态基准应变。公式（1）和（2）中，角频率ω与时间t是相同的，但因为相位差的存在使得${{\varphi }_{1}}\ne {{\varphi }_{2}}$，因此，动应力-应变关系不再是直线，而是在直角坐标系中的长短轴均与坐标轴斜交的椭圆。用于计算阻尼比λ_d的P_d和ε_d可直接从应力-应变图上读取（图 4（b）），但椭圆的长短轴必须通过拟合得到，因此椭圆拟合成为获得阻尼比的关键。

图 4  椭圆拟合迭代前（a）后（b）变化示意图

Figure 4.  Schematic diagrams of ellipse fitting iteration, before(a) and after(b)

下载: 全尺寸图片 幻灯片

理论上，试验过程中应力及应变响应均为时间的函数，是一条光滑的单频谐波曲线，但由于误差的存在，实际得到的曲线除了设定频率的主频外，通常含有其它频率成份的噪声，如图 2（a）和（c）所示。此时，可以通过设置滤波器来滤掉主频以外的成份。为了保证不发生相位的移动，应该采用无相移的滤波器；同时考虑振幅失真最小，我们选用通带内具有最大平坦特性的巴特沃斯滤波器。图 2（b）和2（d）是采用4阶无相移巴特沃斯2Hz低通滤波器滤波后的结果，可见经滤波后曲线光滑了很多，振幅与滤波前变化不大。原则上，可以在图 2（b）和（d）中直接读取初相${{\varphi }_{1}}$和${{\varphi }_{2}}$，然后通过三角函数和差化积相关公式求得待拟合椭圆的相关参数，但是这样做误差较大，因此多数情况下仍采用椭圆拟合的方法来得到椭圆的相关参数。

图 2  应力、应变初始（（a）、（b））及滤波后（（c）、（d））的时间曲线

Figure 2.  Initial time series of stress (a) and strain (b), filtered time series of stress (c) and strain (d)

下载: 全尺寸图片 幻灯片

椭圆拟合的方法有很多种，如蒙特卡罗法、非线性最小二乘法（许正文等，2008；闫蓓等，2008）、参数拟合法（陈基伟，2007a；2007b）和几何法（彭青玉，2003）等。蒙特卡罗法也称统计模拟法，需要借助计算机生成大量样本并对其特征进行统计分析，通常较为费时；非线性最小二乘法直观明了，但是需要求解5个参数，同时还要求二元二次函数满足椭圆的约束条件，事实上是一个二次规划问题，容易因陷入局部极小而不能很好地拟合观测数据；参数拟合法处理已知初相的问题较为方便，但对本文涉及的问题无法线性化，求解存在困难。彭青玉（2003）根据椭圆到两定点距离之和为常数的几何定义出发，提出了椭圆拟合的几何方法。与传统方法相比，该方法直接基于几何定义，因此相对稳定，并且很容易线性化，可操作性强，下文中将对该方法进行详细介绍。

椭圆拟合的几何方法中，需要给定半长、短轴以及长轴与横坐标的夹角。我们的做法是首先将滤波后的应变-应力曲线画在一张图上，然后手动删去明显异常的点，如图 3。

图 3  手动删除异常点前（a）后（b）应力应变曲线

Figure 3.  Stress-strain curves before (a) and after (b) deleted abnormity point manually

下载: 全尺寸图片 幻灯片

由于椭圆拟合问题本质上是非线性问题，初始值越接近于拟合值越容易收敛，计算也越快。我们首先通过二维主成份分析确定长轴的方向。主成份分析通过正交变换将原坐标系进行旋转得到一个新的坐标系，使变量在新坐标系下方差最大，通过确定旋转角度θ就可以得到椭圆长短轴的方向。关于主成份分析理论及计算过程可参考文献（杨淑莹等，2015）。尽管主成份分析的数学推导过程较为复杂，但是利用Matlab仅需要一条命令就可以实现。椭圆中心初始值取所有点的质心（ε_c，P_c），然后可以确定椭圆长轴所在的直线方程。由于长短轴是正交的，且都过椭圆的中心，因此椭圆的短轴方向也随之确定了。椭圆半长轴长度取所有点到短轴距离的最大值a，半短轴长度取所有点到长轴的距离最大值b。椭圆函数是二次函数，可以从实际椭圆的外部或内部逼近拟合，经尝试后发现从内部逼近较容易收敛。为此，取初始半长轴为所有点到短轴距离最大值的0.8倍。初始焦距长度为$c=\sqrt{{{a}^{2}}-{{b}^{2}}} $，两个焦点（ε₁，P₁）和（ε₂，P₂）的位置分别为：

理想情况下椭圆上的点到两焦点距离之和为常数2a，实际中由于误差的存在，椭圆拟合方程表示为：

式中υ为残差。设：

l₁₀对ε₁，P₁的偏导数为：

l₂₀对ε₂，P₂的偏导数具有与公式（6）相同的形式。将（4）式在初值处进行泰勒展开可得：

其中，Δε₁、ΔP₁，Δε₂、ΔP₂分别为两个焦点位置的改正量，Δa为半长轴的改正量。对于（7）式，每一个已知点对应于一个方程，若有N个试验值就可得到一个Gm=d形式的线性方程组，其中G为N×5的系数矩阵，m为焦点位置与长轴改正量这5个系数的列矢量，d为由初始值构成的N元素列矢量，采用最小二乘可以求得5个改正量。实际中，将上一次的初始值加上在该初始值下得到的改正量，作为下一次计算的初始值进行新一轮的计算。通过迭代，不断更新焦点位置和长轴的长度，直至相邻两次的值之差或残差小于一定的阈值，一般5—10步就可以收敛。实际中，为了保证系数矩阵的非奇异，可以加一个微小的阻尼，即采用阻尼最小二乘法。图 4（b）中的拟合椭圆是以图 4（a）的初始椭圆迭代6次的结果。

3.   试验结果

采用残余变形（孔压）较小的若干级循环，在每级循环中选取有代表性的滞回曲线，用本文的方法拟合椭圆，并计算动剪切模量和阻尼比（图 4（b））。然后，采用Hardin-Drnevich双曲线模型（Hardin等，1972）对这一系列的值进行拟合得到动剪切模量G与动剪应变γ_d、阻尼比λ_d的关系，最终将某一粘土样的试验结果绘制成G/G_max-γ_d及λ_d-γ_d的关系曲线，如图 5所示（G_max为最大动剪切模量），不同动剪应变γ_d对应的G/G_max值、λ_d值如表 1所示。

图 5  不同动剪应变下的G/G_max及动阻尼比λ_d曲线

Figure 5.  Curves of G/G_max and λ_d at different dynamic shear strain levels

下载: 全尺寸图片 幻灯片

表 1  试样动三轴试验结果，不同动剪应变下的G/G_max及λ_d

Table 1.  The dynamical tri-axial test results of G/G_max and λ_d at different shear strain levels

施加围压/kPa	土动力特性参数	动剪应变
		5×10－6	1×10－5	5×10－5	1×10－4	5×10－4	1×10－3	5×10－3	1×10－2
200	G/Gmax	0.9937	0.9895	0.9505	0.9058	0.6592	0.4930	0.1633	0.0894
	λd	0.0028	0.0049	0.0210	0.0372	0.0989	0.1246	0.1582	0.1641

下载: 导出CSV 

| 显示表格

从表中可以看出，不同剪应变下，G/G_max及动阻尼比λ_d值在粘性土合理值范围内（孙静等，2004），满足工程项目要求。

从图 5可以看出，本文中椭圆拟合的几何方法能够很好地拟合曲线，是能够应用于工程实践的实用方法。

4.   结论

通过上述试验和分析，得到了不同动剪应变下该粘土样的动剪切模量和阻尼比随剪应变变化的关系曲线。可以看出，本文进行椭圆拟合的几何方法是一种相对稳定、可操作性强的方法。在实际工程中，采用人机交互可有效地控制数据质量。将该方法应用于同一批其他试验土样的滞回曲线拟合，均得到了较好的结果，进一步证明了该方法可以确保土动三轴试验数据结果的可靠性以及各土动力学特性参数取值的准确性。