基于独立分量分析的黏土欧拉数与物理力学指标关系研究

单 浩 张思卿 赵 晔 曹永刚 徐寿政 刘 鑫

(1.徐州工程学院 土木工程学院,江苏 徐州 221018;2.河海大学 岩土力学与堤坝工程教育部重点实验室,南京 210024;3.中交三航局第三工程有限公司,南京 210011)

在工程建设和运行过程中,由于土体强度的不足,公路路基失稳滑坡、基坑塌方等工程事故时有发生,通常可在研究土体的物理力学性质的基础上分析事故发生的机理.近年来,随着科技的进步,人们不再局限于宏观的土工试验方法分析黏土的强度特性,开始采用细观试验方法研究黏土受力变形机理[1-2].一般情况下,土的宏观试验可以测定物理力学指标以指导工程设计和施工,而黏土的抗剪强度受到土体细观结构中颗粒间的联结形态的影响,一般以欧拉数反映联结情况的好坏,通过分析欧拉数的变化有助于揭示黏土受力变形机理[3].蒋明镜[4]在岩土破损力学理论的基础上,通过考虑土体结构性破损的微观机理,建立了结构性土的本构模型;张先伟等[5]对湛江原状黏土与重塑土开展SEM-EDS试验研究,建立了微观结构-化学物质-宏观力学行为的相互关系;刘乐青等[6]对黄土进行了冻融循环试验、无侧限抗压强度试验和电镜扫描试验,分析冻融循环对黄土无侧限抗压强度和微观结构的影响;王静等[7]基于灰色关联理论,分析了微观结构特性与粉砂土的内摩擦角、黏聚力以及粉砂土路基抗剪强度等宏观力学行为的关联性.可以发现,现有研究在黏土细观特性以及与宏观物理力学指标之间关系方面取得丰富的成果,然而鲜见关于欧拉数与黏土物理力学指标间建立直接关系的研究,目前欧拉数仍主要通过开展细观结构试验进行测定,在工程应用中对技术人员的要求较高,不利于工程早期分析黏土受力变形机理.

黏土欧拉数反映了土颗粒联结情况,而物理力学性质指标中,孔隙比反映了土颗粒和孔隙宏观上的体积关系,土颗粒相互联结、相互作用情况直接影响了土体强度指标黏聚力和内摩擦角,因此本文首先开展室内土工试验测定不同黏土的物理力学指标孔隙比e、黏聚力c和内摩擦角φ,以及细观结构试验测定欧拉数E0,然后利用独立分量分析方法建立欧拉数E0与黏土物理力学指标间的关系式,最后对关系式的可靠性和稳定性进行验证,为工程中早期开展黏土受力变形机理分析奠定基础.

1 黏土宏微观试验

路堤填料黏土试样取自广东河惠莞高速公路(土样1)、广佛肇高速公路A 标(土样2)、佛江高速(土样3)、梅平高速(土样4)、江苏南京横江大道(土样5)和河北京秦高速遵化至秦皇岛段(土样6)等高等级公路项目,共获取6组.

1.1 土体孔隙比

明显影响土体细观结构的物理参数主要为孔隙比,因此,首先开展室内土工试验,测定黏土的孔隙比.

针对获取的6组黏土试样分别开展室内土工试验,即烘干法测定含水率w、环刀法测定密度ρ和比重瓶法测定土粒相对密度Gs,并根据式(1)计算土体孔隙比e.

式中:ρw为水在4℃时的密度.

1.2 三轴剪切试验

1)试验仪器

为加快试验进度,采用TSZ-2 全自动三轴仪进行三轴剪切试验.

2)破坏取值标准的选择

根据试验规程,黏土CU 试验应变控制在轴向应变达到20%时停止,当主应力差和轴向应变关系曲线存在峰值时,以峰值点为土体破坏点,无峰值时,取轴向应变的15%对应点作为土体破坏点[8].

3)试验过程

对不同土样分别用围压σ3为100、200、300、400 k Pa进行等压固结,直至孔隙水压力读数接近于零,待固结完成后,进行不排水剪切试验,剪切速率控制为轴向应变0.062 5%/min(即0.05 mm/min).通过试验记录剪切过程的偏应力与轴向应变的关系,得到土体的黏聚力c和内摩擦角φ.

1.3 细观结构试验

土体结构特性由颗粒联结形态决定,而颗粒间的联结带通过常规试验无法获取,因此可以采用具有拓扑性质的欧拉数对联结带的几何形态进行描述,并通过科学合理的试验手段观测细观结构,测定欧拉数.

1.3.1 试验仪器

土体微细结构测试系统由图像采集系统、三维位移控制系统、图片处理及数据提取系统等部分组成,如图1所示.试验中配合使用冷光源作为补光来源,当被观测黏土体受到光线的照射,其表面的反射光经CCD摄像头接受和图像采集卡成像后,通过相配套的Ueye软件在计算机中显像和保存[9-10].

图1 土体微细结构测试系统

1.3.2 测试方法

黏土细观结构测试方法主要分为选择测试土块、确定放大倍率、图像处理和细观特征参数提取等4个部分.

1)选择测试土块:为便于观察土的细观结构,选取土体直径为2 cm 左右.由于用刀切割土块将导致土颗粒和孔隙的破坏,因此采用人工掰开土块后以新鲜凹凸剖面作为观测对象,此时观测剖面具有随机性和代表性.将待观测土块放置在观测平台上,观测过程中尽量避开土体中石英等砾状杂质,选择质地比较均匀的区域作为观测区域.

2)确定放大倍率:通过对比分析不同放大倍数的图像,确定当放大140倍时所采集图像能比较明显反映颗粒与孔隙分布情况.

3)图像处理:通过微调被观测物体与镜头之间的距离来调节聚焦的区域,得到同一视野不同聚焦区域的图片,然后对该图片进行融合,得到全景较清晰的图片,图2给出了土样1融合后的细观结构图片.

图2 土样1细观结构图

4)细观特征参数提取:通过细观特征参数采集软件Geoimage进行图像参数提取,选取最大方差自动取阈的方式对图片进行二值化处理,并提取图片中包含的颗粒联结形态信息.

1.3.3 欧拉数提取

欧拉数在图像变换下可以保持不变,且与几何形状无关.土体中联结形态在细观结构图像中的欧拉数Eo定义为联结部分的数目Ω与孔隙数Nh之差,则计算公式为:

由于根据定义的计算方法仅能进行全局测量计算,计算过程相对繁琐,因此常通过局部性质进行测量计算确定欧拉数,许多学者提出了有意义的计算方法,目前常用的有根据相邻像素值进行计算的4-连通算法和8-连通算法.采用林小竹等[11]基于4-连通算法和8-连通算法提出的一种更为统一的描述欧拉数与图段相邻数关系的局部计算方法:

式中:I表示图像的总行数;N(i)为图像第i行的图段数;V in为第i行第n个图段的相邻数.当n=0时,取V in=1.

2 试验结果及分析

2.1 试验结果

2.1.1 土样物理参数

6组黏土试样含水率、密度、土粒相对密度和孔隙比试验及计算结果见表1.

表1 土样物理参数

2.1.2 三轴试验结果

根据上述CU 试验方法,分别针对6组不同区域的黏土测定抗剪强度指标,见表2.

2.1.3 细观结构试验结果

试验中,取不同的黏土测试土块各3块用于微细观结构测试,同时各测试土块取3个测试点,则可以获取9张同一黏土微细观结构图像,并设定为1组;然后将每组9张微细观图像分别采用Geoimage软件处理,计算每个测点的欧拉数,最后以平均值作为最终的试验结果.黏土欧拉数试验结果见表3.

2.2 试验结果分析

从表1～3可以发现,欧拉数越大,联结带数目越多,颗粒接触较为紧密,凝结趋势较强,宏观上表现为孔隙比越小,强度指标越大;反之,欧拉数越小,联结带数目越少,土体结构较疏松,孔隙比越大,宏观上表现为强度指标越小.表明欧拉数与孔隙比和强度指标有明显的相关性,可以通过数学方法建立4个参数间的相关表达式.

3 欧拉数与黏土物理力学指标关系式

为建立欧拉数与黏土物理力学指标的关系模型,传统方法是将欧拉数作为因变量,物理力学指标作为自变量,建立多元线性回归方程,然而自变量之间容易出现多重共线性问题.为克服传统建模方法的不足,引入独立分量分析(ICA)这一信号处理的新方法[12].通过分析ICA 的基本原理和主要算法,然后根据研究目的采用基于独立分量分析的回归方法建立欧拉数与黏土物理力学指标关系式.

3.1 独立分量分析原理和算法

独立分量分析是在盲源分离技术(BSS)的基础上发展起来的处理多维数据的一种新方法,近年来在信号数据处理方面受到广泛关注.ICA 是针对试验数据,通过建立目标函数运用优化算法分解试验数据,得到相互独立的源数据,具有较强的物理意义[13].

ICA 的一般表示是:设n个独立的影响因素源数据s i(t)组成一个n维矢量,i=1,…,n,S(t)=[s1(t),…,s i(t),…,s n(t)].其经混合矩阵A线性组合成m维的试验数据矢量X(t),A是m×n阵,j=1,…,m,X(t)=[x1(t),…,x j(t),…,x m(t)],则:

除已知s i(t)相互独立外,没有关于S(t)和A的其他先验知识,即无法用具体的数学物理关系进行求解.此时求解一个n×m的解混矩阵W,使得X(t)经过W混合得到的数据是独立源数据的最佳逼近,即Y(t)是S(t)的最佳逼近,Y(t)=[y1(t),…,y n(t)],则

ICA 原理如图3所示.为了求解混矩阵和独立分量,由芬兰学者Aapo Hyvarinen[14-16]在改进和简化的固定点算法得到较为成熟的应用.该算法是基于负熵最大化作为非高斯性度量判据确定的目标函数,采用牛顿迭代原理进行计算,具有收敛速度快的优点,因此又被称为快速ICA(FastICA)算法.

图3 独立分量分析原理图

3.2 独立分量回归

数据分析的主要目的是希望通过合理的数学统计方法预测指标的发展规律,基于ICA 的回归方法称为独立分量回归(ICR),即首先获取黏土物理力学指标试验数据的独立源数据,由于独立源数据去除了向量间的相关性和共线性,从而可以进一步建立多元线性回归模型来达到预测欧拉数的目的.

3.2.1 ICA 处理

用FastICA 算法求X的解混信号Y:

式中:E(X)表示X均值;Q为X-E(X)的白化矩阵;w ik为解混矩阵WQ的第i行第k列元素.i=1,…,n;k=1,…,m.

3.2.2 确定回归系数

建立因变量C y与自变量的解混信号Y的n个分量间的多元线性回归模型,见式(7);利用最小二乘法进行求解系数矩阵D,见式(8).

设自变量与因变量存在多元线性回归模型,即

式中:b0,…,b m为待定系数.则联立式(6)～(9),得如式(10)所示的系数矩阵:

式中:x k为X各分量的均值.

3.3 欧拉数与物理力学指标关系式

取上述第6组土样为对照组,以另外5组的黏聚力c、内摩擦角φ和孔隙比e作为ICA 处理的试验观察数据,设

从而根据FastICA 算法,运用Matlab 编制的FastICA 计算程序,分析X的解混矩阵W和解混信号Y.

X数据各分量的均值分别为c=22.58=9.96=0.77,然后求得中心化后的观察信号的协方差矩阵:

则白化矩阵

敏感指标的白化数据如图4所示.

图4 物理力学指标白化数据

解混后的独立数据与解混后数据维度相等,即解混后仍为3个数据维度.根据FastICA 算法依次进行迭代、正交化和归一化处理,从而求得解混矩阵.

根据式(6)计算独立数据Y,则Y的解混信号,如图5所示.

图5 解混信号图

分别以5组不同黏土的细观结构试验获取的欧拉数为因变量,黏聚力c、内摩擦角φ和孔隙比e的解混信号的3个分量为自变量,以欧拉数和解混信号Y各分量采用多元线性模型建立回归方程,可求得:

由式(13)和式(14)得

设欧拉数与黏聚力c、内摩擦角φ和孔隙比e存在如下关系:

式中:β0～β3为待定系数.

则由式(10),可得欧拉数与黏聚力c、内摩擦角φ和孔隙比e的关系式

4 关系式的验证

基于ICA 方法建立的欧拉数与黏土物理力学指标关系式,可以从两个方面对其可靠性进行验证,一是回归方程的检验,即回归方程和回归方程系数进行相关性和显著性检验;二是以对照组的欧拉数与物理力学指标试验值与模型计算值进行对比分析验证.

4.1 数值稳定性和相关性

回归方程的检验中,首先通过分析在基于最小二乘法确定系数时的解混矩阵的条件数,以判断所求待定系数的数值稳定性,从R检验对ICA 方法建立的回归方程进行可靠性和显著性检验.

在基于ICA 方法建立的关系式时,均采用所获取的独立源信号作为自变量,在确定其待定系数时采用最小二乘法得到以待定系数为未知数的线性方程组.根据数值分析知识,由于试验数据存在一定的变动范围,当未知数的系数矩阵条件数越大时,试验数据的变动为待定系数的确定带来更大的误差.在ICA过程中,所获取的解混信号经过了正交化和归一化处理,因此以待定系数为未知数的线性方程组的系数矩阵的条件数恒等于1,即=1,表明此时获取的待定系数数值最稳定.若直接建立欧拉数与物理力学性质指标的多元线性回归方程,其以待定系数为未知数的线性方程组的系数矩阵条件数为2.67×107,表明方程病态性严重,所求系数的数值稳定性较差,分析原因是物理力学性质指标间存在较为明显的多重共线性,不适合直接建立欧拉数与物理力学指标的多元线性回归方程,这也表明了ICA方法的优越性.

为检验回归关系式的拟合优度,采用R检验相关系数进行判断.本文采用Matlab计算关系式的回归系数,直接获取相关系数R2=0.978 0,计算结果表明回归方程相关性较接近1,拟合效果较好.

4.2 试验值与计算值对比分析

根据试验获取土样1～5的物理力学性质指标,代入式(18),可得由ICA 方法获取表达式的计算结果,从而对比分析关系式计算值与欧拉数试验值的误差,如图6所示.

图6 关系式计算值与欧拉数试验值对比图

从图6可以看出,基于ICA 建立的关系式欧拉数计算值与试验值最大误差发生在土样3,绝对误差为7,相对误差为14.6%,最小误差为0,发生在土样5,表明计算值与试验值误差较小.此外,以本文中第6组黏土对关系式进行验证,根据物理力学性质指标试验结果,代入式(18),模型计算结果为74.5,试验值为74,误差仅为0.67%,表明本文基于ICA 建立的黏土欧拉数与物理力学指标关系式具有较高的准确性,有助于为准确预测欧拉数和分析细观结构特征提供条件.

5 结论

1)欧拉数越大,联结带数目越多,颗粒接触较为紧密,凝结趋势较强,宏观上表现为孔隙比越小,强度指标越大,表明欧拉数与孔隙比和强度指标有明显的相关性,可以通过数学方法建立4个参数间的相关表达式.

2)黏土物理力学指标间存在较为明显的多重共线性,不适合直接建立欧拉数与物理力学指标的多元线性回归方程,基于ICA 建立关系式过程中,对所获取的解混信号经过了正交化和归一化处理,关系式系数矩阵的条件数恒等于1,获取的待定系数值最稳定.

3)基于ICA 建立的黏土欧拉数与孔隙比、黏聚力和内摩擦角关系式具有较高稳定性、相关性和准确性,有助于为准确预测欧拉数和分析细观结构特征提供条件.