作者简介:王杏杏(1992—),女,天津城建大学硕士研究生.研究方向:土的本构关系.E-mail:964479377@qq.com
中文责编:坪 梓; 英文责编:之 聿
1)天津城建大学土木工程学院, 天津 300384; 2)新乡学院土木工程与建筑学院, 河南新乡 453003; 3)大连民族大学土木建筑工程学院, 辽宁大连 116600
Wang Xingxing1, Pan Lin2, Gao Lingxia3, Xia Jinhong2, and Li Shunqun11)School of Civil Engineering, Tianjin Chengjian University, Tianjin 300384, P.R.China2)School of Civil Engineering and Architecture, Xinxiang University, Xinxiang 453003, Henan Province, P.R.China3)School of Civil and Architecture Engineering, Dalian Nationalities University, Dalian 116600, Liaoning Province, P.R.China
geotechnical engineering; pedigree clustering analysis; inter-class distance; loess; microstructure; collapsible; geometric properties
DOI: 10.3724/SP.J.1249.2016.04394
为表述黄土湿陷前后的微结构特征,引入谱系聚类分析法,对采集的21张黄土扫描电镜图像的微结构参数进行研究.根据聚类原理,定义了最长距离法、最短距离法以及Ward法的类间距离,给出了参数间具体的聚类过程.研究表明,对黄土湿陷前后微结构进行评估时,最长距离法和Ward法具有较好的稳定性和可重复性,最短距离法无法客观评价其特征.采用的谱系聚类法,可以综合评估域内颗粒体几何属性方面的差别.
In order to visually describe the microstructure characteristics of loess wet depression, the microstructure parameters of the twenty-one SEM images of loess were studied by clustering analysis method. According to the clustering principle, the three types of inter-class distance in the longest distance method, the shortest distance method and the Ward method were defined respectively, and the specific clustering process was given. The research shows that, when the microstructure characteristics of loess wet depression were evaluated, the longest distance method and the Ward method have good stability and repeatability, but the shortest distance method cannot objectively evaluate its characteristics. The pedigree clustering method can be used to evaluate the difference of the geometrical properties of the particles in the region.
从工程角度出发,研究土的微结构是为了解释某些工程力学现象,以及建立与宏观特性的定量联系[1].在利用微结构解释工程力学现象方面国内外已有了许多成果[2- 4].李兰[5]利用扫描电镜(scanning electron microscope,SEM)技术研究黄土微结构,发现了黄土液化与黄土微结构间的内在联系.高国瑞[6]阐述了黄土湿陷性的本质和失陷变形的机制.宋章等[7]对黄土微结构做了研究,同时分析了对黄土湿陷性状的影响.孙强等[8]讨论了湿陷性黄土的特征和微结构失稳之间的关系. Hu等[9]研究了强夯黄土的力学性能和微观结构变化.谷天峰等[10]利用技术处理后的扫描电镜图像,研究动力作用对Q3黄土微结构参数的影响.在微结构和宏观特性间建立定量联系方面,近几年来展开了一些研究[11-13].田堪良等[14]为了黄土结构性的定量描述,提出了一系列结构性参数.王兰民等[15]结合扫描电镜的相关数据,定量化分析了黄土的震陷性.这些研究揭示了黄土的湿陷性与其微结构间的密切联系.但是,由于研究微结构的方法众多,各方法得到的微结构参数量大,参数携带的信息也多有重复.以致利用微结构参数描述黄土的宏观特性(如湿陷性、结构性等)时存在局限性.
统计学中常用谱系聚类分析法将样品按照亲疏程度的不同进行分类[16].为了利用较少的微结构参数直观表述黄土湿陷前后微结构特征,避免信息重复.本研究引入谱系聚类分析法对采集的微结构参数进行研究.
黄土样品取自陕西省西安市灞桥镇灞河北岸,电镜扫描下的黄土骨架颗粒结构多为支架结构和镶嵌结构,颗粒体间的连接方式多为焊接连接,如图1.
图1 黄土骨架的颗粒结构及连接方式(2 000×)
Fig.1(Color online)Structure and connection of the loess skeleton particles(2 000×)
Fig.1(Color online)Structure and connection of the loess skeleton particles(2 000×)
Leica QWin图像处理系统可以提取域(field)微结构参数和特征(feature)微结构参数.前者提取的是整个分析域的微结构参数,而后者提取的是每个颗粒单元的微结构参数.本研究在Leica QWin图像处理系统提取的域微结构参数的基础上,给出21张黄土SEM照片有关数据的谱系聚类分析.提取的域微结构参数反映的为整个分析域,共有12个,包含9个原始参数和3个导出参数.本研究主要研究9个原始参数.提取所有照片同一分析域的信息,使所得数据具有一般性和可比性.分析域具体大小和位置如图2.图2中截取的分析域为800 μm×600 μm的矩形,其中, pi表示第i个颗粒体; εi表示垂直长度; νi表示水平长度.长度单位为像素,1 μm=20像素.
图2 分析域及参数定义
Fig.2(Color online)Analysis field and parameter definition
Fig.2(Color online)Analysis field and parameter definition
每个颗粒体提取的9个原始参数,组成1个9维向量,9维向量中每个元素指代的含义如下.水平截距为x1, 表示颗粒体水平投影长度之和; 垂直截距为x2, 表示颗粒体垂直投影长度之和; 颗粒体的周长之和为x3; 颗粒总数为x4; 总孔隙面积百分比为x5; 各向异性系数为x6; 填充比为x7; 平均弦为x8; 256灰度自适应阈值为x9[17].
聚类的目的是将数据归类,减少研究对象数目并方便获取对象的相关关系.目前,已有学者将聚类分析方法引入土体研究[18-19].蔡国军等[20]证明在划分土类时,结合聚类分析原理,能更可靠有效地划分地质土层.本研究引入统计学中谱系聚类分析法的最长距离法、最短距离法和Ward法,定义了3种方法的类间距离及其递推公式.
本研究的每个颗粒体都是1个9维向量,且这9个原始参数物理意义和量纲各不相同.为了比较数据间的差别,在系统聚类前,需将每一参数进行标准化.此次谱系聚类分析对象可看作9维空间中的21个样品点,故标准化后所得的参数矩阵应为9×21矩阵,即
a=(a1,a2,a3,a4,a5,a6,a7,a8,a9)T(1)
标准化过程举例说明.每个样品点对应1个9维的列向量(即每个样品点观察9个指标),假设4个样品点数据为
n*1=(1,2,3,4,5,6,7,8,9)T(2)
n*2=(0.8,2.1,3.2,4.2,4.8,6.1,7.2,
8.0,9.3)T(3)
n*3=(1.2,2.2,2.8,4.2,5.3,6.2,6.9,
8.4,8.9)T(4)
n*4=(1.4,2.2,2.7,3.9,5.0,6.5,7.0,
8.2,9.4)T(5)
可得9×4矩阵,按照以下关系
x'ij=(xij-xj)/(Sj)
i=1,2,…,n; j=1,2,…,p(6)
xj=1/n∑ni=1xij, j=1,2,…,p(7)
Sj=[1/(n-1)∑ni=1(xij-xj)2]1/2(8)
其中,xij为矩阵中第i行第j列的元素; x'ij为标准化后矩阵中第i行第j列的元素; Sj为标准差; xj为第j列的元素; n为总行数; p为总列数.将4个样品点中的每个指标标准化,即矩阵中的每行数据标准化,得
n1=(-0.39,-0.78,0.34,-0.50,-0.12,
-0.93,-0.20,-0.78,-0.63)T(9)
n2=(-1.16,0.26,1.24,0.83,-1.09,
-0.46,1.39,-0.78,0.63)T(10)
n3=(0.39,-0.78,-0.56,0.83,1.33,
0,-0.99,1.30,-1.05)T(11)
n4=(1.16,1.30,-1.02,-1.16,-0.12,
1.39,0.20,0.26,1.05)T(12)
基于欧氏距离,定义元素ai与aj之间的距离为lij=d(ai, aj), 即
lij=(∑9k=1(aki-akj)2)1/2(13)
计算4个样品点 n1、 n2、 n3和 n4相互间的距离如表1,组成距离矩阵N(0).
N(0)=[0 4.77 5.64 9.23
4.77 0 12.80 11.28
5.64 12.80 0 10.02
9.23 11.28 10.02 0]
系统聚类方法有很多,选用不同的方法对应不同的距离,现定义最长距离法、最短距离法以及Ward法的类间距离.
1)最长距离法
Dpq=maxi∈Gp,j∈Gqd(ai,aj)(14)
2)最短距离法
Dpq=mini∈Gp,j∈Gqd(ai,aj)(15)
3)Ward法中的类间距离
Dpq=((npnq)/(np+nq)d2<sub>(-overx)p(-overx)q)1/2(16)
4)类的直径
DG=maxi,j∈G dij(17)
其中, Dpq为两类间的距离; Gp和Gq为两个类; np和nq分别为两类中的样品数;(-overx)p和(-overx)q分别为两类的重心; dij为某类中两元素间的距离.
已知类Gp与类Gk之间的距离为Dpk, 类Gq与类Gk之间的距离为Dqk, 则类Gp和Gq合并为新类Gr后,与类Gk间的距离Drk为
1)最长距离
Drk=max{Dpk,Dqk}(18)
2)最短距离
Drk=min{Dpk,Dqk}(19)
3)Ward法中的类间距离
Drk=((np+nk)/(nr+nk)D2pk+(nq+nk)/(nr+nk)D2qk-(nk)/(nr+nk)D2pq)1/2(20)
4)直径
Drk=maxi,j∈Gp∪Gqdij(21)
研究选用最长距离法、最短距离法和Ward法,对300 kPa压力作用下黄土湿陷前后的微结构参数进行评估,具体步骤如下.
1)本次谱系聚类分析针对21个样品,初始时将每一个样品的微结构特征当作独立的一个类,即共有21个类.利用式(13),计算21个类相互间的距离,得到一个21×21的距离方阵,为
D(0)=[0 l1,2 l1,3 … l1,21
l2,1 0 l2,3 … l2,21
l3,1 l3,2 0 … l3,21
l21,1 l21,2 l21,3 … 0](22)
2)从21阶方阵D(0)中选择最小非零距离,将其对应的类Gp和Gq合并为一个新类Gr, 同时用新类Gr替换D(0)中的类Gp和Gq, 新类与剩下的类之间的距离组成新的行列,从而得到一个新的20阶方阵D(1).
3)从D(1)出发,按照所选用的系统聚类方法,结合步骤2),将20阶方阵更新为19阶距离方阵D(2). 按照上述方法递推,当所有的类全部聚为一类时,聚类结束.
4)由具体聚类情况,按顺序为合并的类编号,最后绘制聚类谱系图.
基于上述4个样品n1、 n2、 n3和n4的数据,描述最长距离法、最短距离法和Ward法的聚类过程.
将 n1、 n2、 n3和 n4每个样品当作独立的一类,由表1可知所有距离中l12=4.77最小.由于 n1和 n2距离最短,将其合并为新类,按顺序定为第5类.
按式(14)和(18)计算类5与当前类的距离
d53=max{d13,d23}=max{5.64,12.8}=12.8
d54=max{d14,d24}=
max{9.23,11.28}=11.28
得到新距离矩阵为
N(1)=[0 10.02 12.80
10.02 0 11.28
12.80 11.28 0](23)
继续合并距离最短的两类n3和n4为新类,定为第6类,计算新类与当前类的距离
d65=max{d35,d45}=max{12.8,11.28}=12.8
此时所有元素已聚为一个大类,聚类结束.
与最长距离法类似,将距离最短的两类,即 n1和 n2, 合并为新类,按顺序定为第5类.按式(15)和(19)计算新类5与当前类的距离,如
d53=min{d13,d23}=min{5.64,12.8}=5.64
d54=min{d14,d24}=min{9.23,11.28}=9.23
得到新距离矩阵为
N(1)=[0 10.02 5.64
10.02 0 9.23
5.64 9.23 0](24)
继续合并距离最近的两类为新类.因d35=5.64最小,合并为新类第6类,其与当前类的距离如下
d65=min{d35,d45}=min{5.64,9.23}=5.64
此时所有元素已聚为一个大类,聚类结束.
结合表1,利用式(16)可求得Ward法中的d2ij, 表示为距离矩阵
N(0)=[0 2.38 2.82 4.62
2.38 0 6.40 5.64
2.82 6.40 0 5.01
4.62 5.64 5.01 0](25)
由于d212=2.38最小,合并为新类,按顺序定为第5类,其与当前类的距离为
d253=2/3d213+2/3d223-1/3d212=2/3×2.82+2/3×6.40-1/3×2.38=5.35
d254=2/3d214+2/3d224-1/3d212=2/3×4.62+2/3×5.64-1/3×2.38=6.05
得到新距离矩阵
N(1)=[0 5.01 5.35
5.01 0 6.05
5.35 6.05 0](26)
由于d234=5.01最小,合并为新类第6类,其与当前类的距离为
d265=3/4d235+3/4d245-2/4d234=3/4×5.35+3/4×6.05-2/4×5.01=6.05
此时所有元素已聚为一个大类,聚类结束.
对黄土微结构参数进行分析时,选取不同的谱系聚类方法,得到的聚类结果不同.本研究选取3种谱系聚类方法分析黄土微结构参数.
采用最长距离法对21个样品微结构参数进行分析时,具体步骤类比4个样品示例.首先利用式(13),计算21个类相互间的距离,得到一个21×21的距离方阵D(0). 从D(0)中可以发现类B6和B7对应的距离最短,将它们聚为新类,命名为CL20.此时原有的21阶方阵D(0)更新为20阶方阵D(1). 从D(1)中观察到类B8和B9对应距离最短,聚为新类,命名为CL19.聚类完成后方阵D(1)更新为19阶方阵D(2).
比较D(2)中的元素,发现类CL19和B10间的最长距离最小,将其聚为新类,命名为CL18.聚类完成后方阵D(2)更新为18阶方阵D(3). 观察D(3), 发现类B4和B5对应最长距离最小,聚为新类,命名为CL17.按照此种方法类推,直到所有的类全部聚为一类为止.聚类完成后,根据具体的聚类情况,绘制聚类谱系图,如图3.
图3所示的聚类结果中,A1至A10表示的是黄土湿陷前的样品微结构参数,可以聚为一类(CL2); B1至B11为黄土湿陷后的样品微结构参数,可以聚为一类(CL3).从图3可发现黄土湿陷前后的微结构特征聚类效果明显.也就是说,基于最长距离法的谱系聚类分析可以反映湿陷这一特征,该类间距离可以作为合成微结构参数评价黄土湿陷性.为使图形清晰,图中略去任意两类间的最长距离,仅给出SAB表示湿陷前后两大类间的最长距离.
结合4个样品的示例以及类比最长距离的聚类过程,利用最短距离法分析21个样品的微结构参数,聚类谱系图如图4.
观察图4,发现聚类结果没有规律,湿陷前的样品微结构参数和湿陷后的样品微结构参数不能各聚为一类,聚类结果与实际差距较大.主要原因是土颗粒体微结构间以及不同类间的距离无法依据最短距离来评价.因此最短距离法无法客观评估黄土湿陷前后的微结构特征.
本研究利用谱系聚类分析法表述黄土微结构特征,根据聚类原理,定义了最长距离法、最短距离法以及Ward法的类间距离,并给出了参数的聚类过程.由聚类结果可得:基于最长距离法和Ward法的聚类方法具有较好的稳定性和可重复性,而最短距离法无法客观评价微结构特征.本研究给出的谱系聚类法,可以用于综合评估域内颗粒体在几何属性方面的差别.
深圳大学学报理工版
JOURNAL OF SHENZHEN UNIVERSITY SCIENCE AND ENGINEERING
(1984年创刊 双月刊)
主 管 深圳大学
主 办 深圳大学
编辑出版 深圳大学学报理工版编辑部
主 编 阮双琛
国内发行 深圳市邮电局
国外发行 中国国际图书贸易集团有限公司(北京399信箱)
地 址 北京东黄城根北街16号
邮 编 100717
电 话 0755-26732266
0755-26538306
Email journal@szu.edu.cn
标准刊号 ISSN 1000-2618
CN 44-1401/N