聚类算法在竹种造纸适应性判定中的应用1)

丁梦雅 余雁王 汉坤 郑文瑞 韩亚鲁 金秀 李绍稳

(安徽农业大学，合肥，230036)(福建农林大学)(国际竹藤中心)(安徽农业大学)

目前世界范围内造纸原料主要是木材纤维原料，木浆比占90%以上。随着造纸需求增加和木材供应不足，竹材成为造纸纤维原料的第二供应者。我国有丰富的竹林资源，分布区域广泛。人工经营的竹林面积近600万hm2，其中毛竹占420万hm2。据统计资料显示，可用于纸浆造纸的竹林约140万hm2，其中集约经营的梁山慈竹、慈竹、硬头黄竹等约20万hm2，主要集中在四川省泸州市、宜宾市和贵州省赤水市，年产竹材约400万t[1]。全球有三大竹区，分别是亚太竹区、美洲竹区、非洲竹区。亚太竹区是世界上最大的竹区，面积约1 400万hm2，其中中国、印度、日本、缅甸、泰国、菲律宾、孟加拉等国是竹子资源较丰富的国家。截至2018年，印度竹浆年产量近50万t；缅甸、泰国、菲律宾等国，单个竹材制浆厂的年产量不超过10万t。美洲竹区的面积约160万hm2，其中以丛生竹为主，且多为小型低矮竹种，经济价值低。非洲竹区面积约150万hm2，其中乡土竹种很少，而东非的马达加斯加岛的竹种较为丰富，有11属40种[2]。然而，全球可用于纸浆造纸的竹材只局限于极少数竹种，大多数竹种的造纸性能并不知晓，也不为造纸工业所用，造成了资源浪费。并且随着可用于竹林的土地资源越来越紧张，靠扩大竹林种植面积增加竹材产量不太现实。只能通过选育优良竹种提高增益，扩大竹材产量。

近年来，在制浆造纸技术方面，国内外已有较多新成果、新技术[3-6]。但是，对竹种制浆造纸性能筛选技术和方法方面的研究较少，使得选育优良竹种陷入了困境，成为亟待解决的难题。为此，本研究依据“十二五”国家科技支撑计划课题“竹藤种质资源基础数据库构建”，以“竹藤种质资源数据库”中126个竹种为基础数据，应用K均值聚类算法(K-means算法)、小批量K均值聚类算法(Mini Batch K means算法)、凝聚层次聚类算法，依据竹种酸不溶木质素、苯醇提取物、综纤维素、基本密度、纤维长度、纤维双壁厚、纤维壁腔比7个特征对竹种筛选分级；评估不同聚类算法的有效性以及不同竹种的造纸适应性，选出一批优良造纸竹种，旨在为合理利用竹林资源提供参考。

1 研究方法

1.1 数据来源

本研究数据来源于“十二五”国家科技支撑计划课题“竹藤种质资源基础数据库构建”项目，由国际竹藤中心提供竹种材性数据，由安徽农业大学开发竹藤种质资源数据管理平台(以下简称“管理平台”)并构建竹藤种质资源数据库(以下简称“数据库”)。竹种材性数据采集地点，来自安徽太平、中国科学院西双版纳热带植物园、云南金平、四川都江堰、广西南宁等地；竹种材性指标由国际竹藤中心的专家核定。

管理平台，目前共录入43个竹属，551个竹种，由于跨科属种类别大，竹种间材性差异较大。竹藤专家参照GB/T 15780—1995《竹材物理力学性质实验方法》标准，采用物理化学等实验与测量方法，获得了大量数据，这些数据是本文建模工作的主要依据。数据库中，由竹种的化学性质、物理性质、纤维形态的部分材性指标的数据量统计(见表1)可见，除了酸溶木质素只有54条记录，其他几项指标数据量都在100条以上。木质素由酸不溶木质素和酸溶木质素组成，由于酸溶木质素数据量过少，本研究采用酸不溶木质素替代木质素研究不同竹种之间木质素的差异性。

表1 竹种部分材性数据统计

材性数据分布在数据库中化学成分、物理性质、纤维形态3个表中，由于每个表数据量不一致，合并为1个表后会产生部分缺失值。为保证数据中缺失值尽可能少，本研究将合并后超过3条材性指标缺失的竹种删除，最终保留了126个竹种。此时，基本密度、纤维双壁厚、纤维长度、纤维壁腔比分别缺失33、48、8、22个数据；酸不溶木质素、苯醇提取物、综纤维素没有缺失值。竹种数据进行分析之前，用各指标值的中位数填充各自缺失值，将数据补充完整。

1.2 聚类算法确定及评价方法

根据竹种的纤维形态、化学性质、物理性质将竹种分类，结合造纸适应性评价指标选出优良造纸竹种。聚类分析法可以根据多条性质的相似性程度对样本进行分类。最终目的是将相似性大的样本聚到一个簇，相似性小的样本聚到不同的簇。本研究采用了6种聚类算法：K均值聚类算法、小批量K均值聚类算法、利用层次方法的平衡迭代规约和聚类算法(BIRCH聚类算法)、谱聚类算法、凝聚层次聚类、高斯混合模型聚类算法。

——K均值聚类算法。通过不断更新聚类中心和分配样本点的过程，最终将样本分为几个不同的簇，实现簇内样本点紧凑，簇间样本点分离[7]。本研究采用K-means++算法初始化聚类中心，以解决传统K均值聚类算法收敛到局部最优的问题。

——小批量K均值聚类算法。是K均值聚类算法的一种改进算法，其在分配样本点给聚类中心时，随机抽取少批量的样本点，更新聚类中心，因此聚类速度比K均值聚类算法快，适合大规模数据集的聚类。

——利用层次方法的平衡迭代规约和聚类算法。是适合大规模数据集的聚类。与小批量K均值聚类算法不同，小批量K均值聚类算法适合于样本量大、数据特征维度高、类别数小的情况，而利用层次方法的平衡迭代规约和聚类算法适合于样本量大、数据特征维度低、类别数大的情况。

——谱聚类。以无向权重图的思想聚类数据集中的样本点，距离较远的两点之间边的权重较低，距离较近的两点之间边的权重较高，通过切图让子图内边的权重之和尽可能高，不同子图间边的权重之和尽可能低，从而达到聚类的目的。

——凝聚层次聚类算法。通过不断迭代循环地找距离最近的两个簇进行合并，最终形成一棵聚类树，簇与簇之间的相似度也在树的生成过程中得到。该算法可以不用预先给定聚类个数，避免了聚类个数调参的问题。

——高斯混合模型聚类。通过高斯模型估计样本数据的概率密度分布，最后得到样本数据属于某个簇的概率，概率最大的簇即为该样本所属的簇。

聚类有效性指标是用来评价聚类算法好坏的一种适应性指标。本研究在聚类前不知道竹种的分类情况，即数据集的分类标签未知，故选用内部有效性指标评价聚类算法的有效性。聚类算法在有效性指标达到最优时取得的聚类个数称为“最佳聚类数”[8]，此时聚类算法的性能最好。内部有效性指标由紧凑度、分离度两个指标构成，紧凑度衡量一个簇中不同样本之间的距离是否足够近、分离度衡量不同簇之间的距离是否足够远。本研究重点研究竹种数据不同簇之间的分离度，因此采用了仅考量分离度的指标(决定系数(R2))以及综合考量紧凑度和分离度的指标(包括卡林斯其-哈瑞本斯(CH)指标、戴维森堡丁指数(DBI)、轮廓系数)。

——决定系数。由簇之间距离平方和与所有数据距离平方和的比值得到，该指标在“肘点”[9]处取得最佳聚类数。

——卡林斯其-哈瑞本斯指标。通过计算簇中样本点与簇中心的距离平方和，度量簇内的紧凑度；通过计算子簇中心点与数据集中心点距离平方和，度量数据集的分离度；通过分离度与紧凑度的比值得到指标值，并且在指标最大值处取得最佳聚类数。

——戴维森堡丁指数。通过计算一个簇C与其他簇的相似度，将最大值作为簇C的相似度，最后平均所有簇的相似度即为指标值。该指标在最小值处取得最佳聚类数。

——轮廓系数。计算方法为S(i)=[b(i)-a(i)]/max{a(i),b(i)}，式中的a(i)为样本i到同簇其他样本的平均距离、b(i)为样本i到其他簇的样本的平均距离。轮廓系数指标最大值对应的聚类数为最佳聚类数。

1.3 纸浆造纸适应性的评价指标设计

苏文会等[10]认为，综纤维素质量分数较高而木质素质量分数较少或中等的造纸材，其纤维得率高、易打浆、化学药品消耗少。范友华等[11]在对湖南益阳毛竹的研究中发现，综纤维素质量分数高于70%、酸不溶木质素质量分数低于26%的竹种，适合做造纸优良竹材。赵燕等[12]认为，苯醇提取物质量分数过多，在制浆时会增加化学药品的消耗。杜复元等[13]认为，基本密度是影响纸浆得率的重要因素，一般认为基本密度介于0.4～0.6 g/cm3之间的竹种适合做优良造纸原料。马灵飞等[14]认为，细长的纤维能增加纸张的耐折度、耐破度、撕裂度。国际木材解剖协会将纤维按长度大于1.6 mm的纤维称为长纤维[15]。王菊华[16]认为，纤维壁薄且壁腔比小的原料，纤维结合强度高，抄制的纸张结合强度大，可压扁性好，通常认为壁腔比小于1的纤维适合造纸；而陈富枢[17]在实验中表明，竹种纤维壁腔比小于2.24，同样可以用作纸浆造纸。参考上述文献，本研究设计了竹种造纸的适应性评价指标(见表2)。本研究选定的126个竹种的苯醇提取物均值是3.669%，故将苯醇提取物质量分数低于等于3.669%的竹种视为适合造纸的优良竹种，高于3.669%的竹种视为不适合造纸的劣等竹种。

综合上述分析，本研究选用数据库中126个竹种，每个竹种包含酸不溶木质素、苯醇提取物、综纤维素、基本密度、纤维长度、纤维双壁厚、纤维壁腔比7个特征，进行聚类分析，结合表2，筛选出优良造纸竹种。

表2 造纸适应性评价指标值

2 结果与分析

2.1 应用聚类算法对竹种的聚类结果

对竹种数据进行缺失数据填补、异常数据处理等数据整理后，选定聚类个数在2～7之间，分别进行K均值聚类、小批量K均值聚类、利用层次方法的平衡迭代规约和聚类算法、谱聚类、凝聚层次聚类、高斯混合模型聚类，每个算法通过决定系数、卡林斯其-哈瑞本斯指标、戴维森堡丁指数、轮廓系数有效性指标评价聚类效果(见图1)。

由图1(a)可见：K均值聚类算法、小批量K均值聚类算法、利用层次方法的平衡迭代规约和聚类算法、凝聚层次聚类算法在聚类个数为6时，决定系数指标取得“肘点”，即竹种数据聚为6类时4个算法的聚类有效性最大，此时最佳聚类数为6。综合图1中4个子图可见：小批量K均值聚类算法在聚类个数为6时，决定系数指标取得“肘点”、卡林斯其-哈瑞本斯指标取得最大值、戴维森堡丁指数取得最小值、轮廓系数取得最大值，4个有效性指标均达到最优，表明小批量K均值聚类算法依据4个聚类有效性指标的最佳聚类数均为6。在图1(b)中，K均值聚类算法、小批量K均值聚类算法、凝聚层次聚类算法、利用层次方法的平衡迭代规约和聚类算法在聚类个数为6时，指标达到最大值，且K均值聚类算法的卡林斯其-哈瑞本斯指标值最高，达到63.568，聚类有效性最好。而K均值聚类算法、小批量K均值聚类算法、凝聚层次聚类算法的卡林斯其-哈瑞本斯指标值很相近，远高于谱聚类算法、高斯混合模型算法。

因此，本研究主要讨论K均值聚类算法、小批量K均值聚类算法、凝聚层次聚类算法对竹种数据的聚类分析。上述3种聚类算法依据4个聚类有效性指标得到的最佳聚类数见表3。

表3 3种聚类算法依据4个有效性指标得到的最佳聚类数

将K均值聚类算法、小批量K均值聚类算法、凝聚层次聚类算法分别应用到竹种数据进行聚类分析，得到聚类结果。由表4可见：3种聚类算法得到的竹种聚类情况大致相似，大多数竹种在3个聚类算法的分析中被聚为同一类，如甜笋竹、金丝毛竹、黄槽毛竹等11个竹种始终属于同一类，糯竹、空竹、泰竹等51个竹种属于一类，油簕竹、酸竹等19个竹种属于一类，斑竹自成一类，泡竹、茶竿竹、版纳甜龙竹属于一类，毛竹、思劳竹等28个竹种属于一类。说明这些竹种间的纤维形态、理化性质具有很多相似之处。

表4 3种聚类算法的聚类结果

续(表4)

A1为K均值聚类算法、A2为小批量K均值聚类算法、A3为利用层次方法的平衡迭代规约和聚类算法、A4为谱聚类算法、A5为高斯混合模型算法、A6为凝聚层次聚类算法。

2.2 竹种的造纸适应性判定

竹种数据经过3种聚类算法分析后均聚为6类。为使竹种纸浆造纸适应性评价较为客观，以各类材性指标的平均值为依据，综合判定各类竹种的造纸适应性等级(见表5)。

表5 3种聚类方法聚类后各类指标的均值

由表5可见：K均值算法聚类结果中，第4类只有1个竹种(是斑竹)，综纤维素质量分数为18.13%，远低于竹种造纸纤维标准70%。第5类有3个竹种(是泡竹、茶竿竹、版纳甜龙竹)，纤维双壁厚均值为32.17 μm，纤维壁腔比均值为4.76，纤维壁腔比值远高于2.24。故认为第4类的竹种、第5类的竹种不适合造纸，在造纸适应性判定里不考虑这两类竹种。

结合表2中竹种造纸适应性评价指标，判定K均值算法聚类结果的第1、2、3、6类竹种的造纸适应性。然后分别进行统计，6类材性指标均达到设定值或者有5类达标且纤维细长的视为Ⅰ级造纸竹种，有4类或3类达标的视为Ⅱ级造纸竹种，只有2类或1类达标的视为Ⅲ级造纸竹种(见表6)。同理可得出，应用小批量K均值聚类算法、凝聚层次聚类2种算法的竹种造纸适用性的评定情况(见表6)。

表6 3种聚类方法的竹种造纸适应性评定结果

由表6可见：K均值聚类算法评为Ⅰ级的造纸竹种为第2类(糯竹、空竹等)，小批量K均值聚类算法评为Ⅰ级的为第1类(毛竹、思劳竹等)、第5类(糯竹、空竹等)，凝聚层次聚类评为Ⅰ级的为第2类(油簕竹、酸竹等)、第5类(糯竹、空竹等)。K均值聚类算法评为Ⅱ级的为第3类(油簕竹、酸竹等)、第6类(毛竹、思劳竹等)，小批量K均值聚类算法评为Ⅱ级的为第3类(油簕竹、酸竹等)，凝聚层次聚类评为Ⅱ级的为第3类(毛竹、思劳竹等)。K均值聚类算法、小批量K均值聚类算法、凝聚层次聚类评为Ⅲ级的均为甜笋竹、金丝毛竹等。由于纤维双壁厚没有确切的划分范围，且一般纤维壁腔比小的竹种纤维双壁厚小，故表6中不含纤维双壁厚指标的评定。

虽然3种算法的聚类结果不完全一致，但是大部分竹种的分类情况是一致的。综合竹种的造纸适应性评定，3种算法均将糯竹、空竹等51个竹种评为Ⅰ级造纸竹，而甜笋竹、金丝毛竹等11个竹种为Ⅲ级造纸竹。云南甜龙竹、短穗竹，在K均值聚类算法、小批量K均值聚类算法聚类时评为Ⅲ级，在凝聚层次聚类时评为Ⅱ级，故本研究将这2个竹种划分为Ⅲ级造纸竹。其余62个竹种，依据3种聚类算法评为Ⅱ级造纸竹2次、Ⅰ级造纸竹1次，故本研究将这些竹种划分为Ⅱ级造纸竹。为此，本研究给出造纸竹种的综合评定结果(见表7)。

表7 造纸竹种的综合评定结果

3 讨论

由竹种造纸适应性的评定可知，应用K均值聚类算法评定的，Ⅰ级造纸竹种的6个材性指标均达到设定值、Ⅱ级造纸竹种达标了4个、Ⅲ级造纸竹种仅达标1个。应用小批量K均值聚类算法评定的，Ⅰ级造纸竹种达标5个、Ⅱ级造纸竹种达标了4个、Ⅲ级造纸竹种达标1个。应用凝聚层次聚类算法评定的，Ⅰ级造纸竹种达标5个或6个、Ⅱ级造纸竹种达标了3个、Ⅲ级造纸竹种达标2个。

可以看到，小批量K均值聚类算法的Ⅰ级造纸竹种存在苯醇提取物略高、基本密度略低2个缺陷，包含部分化学处理困难、得浆率低的竹种。凝聚层次聚类算法同样是在Ⅰ级造纸竹种的选择上有不足，其第2类和第5类均评为Ⅰ级，而第2类综纤维素略低，包含部分纤维得率低、不易打浆的竹种。K均值聚类算法分级最清晰，其Ⅰ级造纸竹种明显比Ⅱ级竹种造纸适应性高，Ⅲ级竹种明显不适合造纸。上述分析表明，K均值聚类算法聚类效果最好，达到了类内相似度大、类间差异性大的目的。这与前面表述的“K均值聚类算法的卡林斯其-哈瑞本斯指标得分最高，算法有效性最高”的结果一致。

由张喜[18]对贵州竹种的纤维形态和造纸性能的研究表明，苦竹、慈竹、吊丝单、早园竹都是优良的造纸竹种，这4个竹种均在本研究的Ⅰ级造纸竹种中。陈其兵等[19]通过实验分析认为，麻竹、绿竹、绵竹、料慈竹是优良造纸竹种，同样存在于本研究的Ⅰ级造纸竹种中。刘晓波等[20]研究的适合做造纸竹种的巨龙竹，在本研究的Ⅱ级造纸竹种可以找到。赵燕等[12]认为，撑绿竹适合纸浆造纸，其在本研究的Ⅱ级造纸竹种可以找到。苏文会等[10]研究的造纸竹种大木竹，在本研究的Ⅱ级造纸竹种可以找到。上述分析表明，本研究采用K均值聚类算法对竹种的造纸适应性评价与多数学者的研究相吻合，说明K均值聚类算法对探索未知造纸性能的竹种是一个切实可行的方法。通过酸不溶木质素、苯醇提取物、综纤维素、基本密度、纤维长度、纤维双壁厚、纤维壁腔比7个造纸适应性评价指标，选出了较为可靠的优良造纸竹种，表明竹种这7个指标是评价其造纸特性的关键因素。

4 结论

本研究重点研究了K均值聚类算法、小批量K均值聚类算法、凝聚层次聚类算法对竹种数据聚类的有效性，结合不同竹种纤维形态、化学性质、物理性质适用于造纸的指标，对竹种用于纸浆造纸的适应性分级，选出优良造纸竹种。最后将126个竹种评定出Ⅰ级、Ⅱ级、Ⅲ级，Ⅰ级竹种的造纸性能最好、Ⅱ级竹种的造纸性能次之、Ⅲ级竹种的造纸性能最差不适合做造纸纤维原料。

根据本研究，造纸工业可以优先发展Ⅰ级竹种作为制浆造纸原料，其具有酸不溶木质素含量低、综纤维素含量高、苯醇提取物含量低、基本密度适中、纤维长度长、纤维壁腔比小的特点，适合制造高撕裂度、高耐破度、高耐折度的优质纸张。Ⅱ级竹种的苯醇提取物和纤维壁腔比略高，可以在具备良好排污能力的造纸厂作为备用造纸原料。Ⅲ级竹种的苯醇提取物高、酸不溶木质素高、基本密度大、纤维长度小、纤维壁厚度大，不适合造纸。