2015, Vol. 51
文章信息
- 任洪娥, 沈雯雯, 白杰云, 官俊
- Ren Hong, Shen Wenwen, Bai Jieyun, Guan Jun
- 基于最小二乘法的落叶松木粉粒径数学建模与分析
- Mathematical Models and Analysis of Particle Size of Coniferous Wood Flour Based on the Least Squares Method
- 林业科学, 2015, 51(4): 164-170
- Scientia Silvae Sinicae, 2015, 51(4): 164-170.
- DOI: 10.11707/j.1001-7488.20150421
-
文章历史
- 收稿日期：2014-04-30
- 修回日期：2014-12-26
-
作者相关文章
2. 黑龙江省林业智能装备工程研究中心 哈尔滨 150040;
3. 哈尔滨工业大学计算机科学与技术学院 哈尔滨 150001
2. Forestry Intelligent Equipment Engineering Research Center, Heilongjiang Province Harbin 150040;
3. School of Computer Science and Technology, Harbin Institute of Technology Harbin 150001
木粉作为一种新型的节能环保原料，用途广泛。在木材工业上，将木粉引入复合材料和高分子材料领域(Kobayashi et al.，2008；朱再胜等，2012；Palaniandy et al.，2009)，木塑复合材料随之产生。木塑复合材料化学稳定性好、强度高，兼具原木材料和高分子材料的优点，因此得以飞速发展(宋永明等，2011；2012)。颗粒的大小和形状是木塑复合材料最重要的物理特性表征量。 生瑜等(2012)研究发现，随木粉粒径的减小，木塑复合材料的拉伸强度、弯曲强度呈现下降趋势，但冲击强度呈上升趋势；宋丽贤等(2013)研究发现，复合材料的拉伸强度和冲击强度随着木粉粒径的减小先升高后降低。这些研究大多集中在粒径大小对于其力学性能的影响上，而很少研究造成力学性能改变的根本原因 。李兰杰等(2006)研究发现，不同粒径的木粉在粉碎过程中所承受的力不同，而不同粒径的木粉具有不同的表面粗糙度和长宽比，大的长宽比有利于材料力学性能的提高(Nourbakhsh et al.，2010)。Guo(2012)研究了425 μm以下木质颗粒平均长宽比与粒径之间的关系，发现长宽比随着颗粒的增大而增大，但没有研究425 μm以上颗粒长宽比的变化规律。系统建立粒径与长宽比的数学模型，可揭示其变化趋势，分析其变化原因，从模型中可探索最佳的木粉粒径，对提高木塑材料的冲击性能具有重要作用。因此本文从长宽比出发，探索了长宽比随粒径的变化趋势，为建立长宽比与力学性能之间的定量关系提供参考。
本文以落叶松(Larix)木粉为研究对象。首先，利用显微图像测算木粉管胞的平均长度、平均宽度以及木粉颗粒的粒径，同时利用数字图像处理方法求得木粉颗粒的长宽比、矩形度：将原始木粉图像由RGB颜色空间转换到Lab颜色空间，对其b分量进行平滑滤波去噪后再进行聚类分割，获得木粉轮廓；利用主轴法计算木粉的最小外接矩形，提取木粉的长宽比和矩形度2个形状特征。然后，分别利用多项式、高斯和傅里叶函数对木粉粒径和木粉长宽比进行数据拟合，通过比较，选用高斯函数得木粉粒径与木粉长宽比的数学模型，根据得出的拟合曲线方程计算出其二阶导数。最后，结合二阶导数与木粉管胞数据对高斯模型进行分析与讨论。
1 材料与方法落叶松木粉图像由光学显微镜Olympus-BX43和Olympus-SZX7配合华通智能的eHAI500工业相机拍摄获得，放大倍率为物镜4×，目镜10×，图像分辨率为2 448×1 920。木粉颗粒的管胞长度、宽度以及木粉粒径大小利用显微镜获取，木粉几何特征利用数字图像处理技术提取。
1.1 管胞长宽度测量针叶材体积的90%以上是由纵行的纤维状厚壁细胞——管胞构成的。显微镜下观察木粉颗粒，管胞两端尖锐，排列整齐，二维图像横断面呈扁平状(图 1)。利用显微镜依次测量每个木粉颗粒中完整管胞的长度和宽度，求取平均值，即为整体目标木粉管胞的平均长度和平均宽度。部分数据见表 1。
 |
图 1 成熟管胞的长度和宽度 Fig. 1 Length and width of the mature tracheid |
|
粒径是单个颗粒几何尺寸的表征量。不同的用途及研究目的采用不同的测量方法，本文中使用的木粉是通过筛选法获得的，其最大粒径为筛孔通过率为100%的最小标准筛孔尺寸，筛分粒径是颗粒可以通过筛网的筛孔尺寸，以1英寸(25.4 mm)宽度的筛网内的筛孔数表示。因此本文中木粉的粒径以其最小外接矩形的长来表示(图 2)。该参数通过光学显微镜测量所得，本研究中木粉的粒径范围为1 100～30 μm。部分数据见表 2。
 |
图 2 木粉颗粒的长和宽 Fig. 2 Length and width of the wood flour particles |
|
|
颗粒形状指一个颗粒的轮廓边界或表面上各点的图像。木粉形状特征通常用二维图像投影面的轮廓曲线来描述。范长胜等(2013)研究发现，木粉形状一般呈片状和纤维状(图 2)。由于木粉此形状特征以及仪器测量的局限性，本文不予考虑木粉厚度，而依据单颗粒木粉的垂直平面投影像的轮廓曲线，主要研究木粉的长宽比(三维意义下又称长径比)和矩形度2个形状特征。
设木粉几何面积为A，最小外接矩形的长为l，宽为b，面积为Sr。
1)长宽比是木粉长轴和短轴的比值，反映木粉的各向异性，用k表示：
| $k=l/b$ | (1) |
2)矩形度是木粉的投影面积与其最小外接矩形面积之比，反映木粉对其外接矩形的充满程度以及形状规则度，用FA表示：
| ${{F}_{a}}=A/{{S}_{\text{r}}}$ | (2) |
考虑到获得的彩色图像受光照影响很大，并且木粉主要表现为黄色，而Lab颜色空间中，L分量表示光照亮度，a分量表示从红色至绿色的范围，b分量表示从蓝色至黄色的范围。因此，本文在Lab颜色空间下，依据b分量进行聚类分析，对木粉进行分割，获得其二值图像(Deng et al.，2013；张冰等，2014；Mok et al.，2012)；利用数学形态学消除二值图像中的噪声块和噪声孔，最后对标记图像进行特征提取(王书涛等，2012；张法全等，2013；卢蓉等，2010)，如图 3。部分数据见表 2。木粉几何特征提取算法步骤如下：
 |
图 3 木粉图像处理过程 Fig. 3 Process of the flour image |
1)输入木粉的原始图像，将RGB颜色空间转为Lab颜色空间，提取b分量。
2)用3×3模板的中值滤波器对b分量图像平滑去噪；用K-means算法进行聚类分析:
① 随机选取k(这里k=2)个初始聚类中心μ1，μ2；
② 定义聚类准则函数：
| $J=\sum\limits_{j=1}^{k}{\sum\limits_{x\in C}{\left\| x-{{\mu }_{j}} \right\|h}}$ |
③ 重复下面过程直到准则函数收敛：对于每一个样例i，计算其应属于类
| ${{C}^{\left( i \right)}}:\arg \underset{j}{\mathop{\min }}\,{{\left\| {{x}^{\left( i \right)}}-{{\mu }_{j}} \right\|}^{2}}$ |
对于每一个类j，重新计算该类质心
| ${{\mu }_{j}}:=\frac{\sum\nolimits_{i=1}^{k}{1\left\{ {{c}^{\left(i \right)}}=j \right\}{{x}^{\left(i \right)}}}}{\sum\nolimits_{i=1}^{k}{\left\{ {{c}^{\left(i \right)}}=j \right\}}}$ |
聚类为2类，得木粉的二值图像A。
3)选取结构元素B，B为5×5的方形结构；运用B对A进行先开启后闭合的数学形态学运算，得图像C：
| $C=\left\{ \left[ \left(A\Theta B \right)\oplus B \right]\oplus B \right\}\oplus B=\left(A\circ B \right)\cdot B$ |
4)利用八连通区域标记算法对图像C中每个独立的木粉颗粒进行标记，得标记图像D。
5)采用目标区域像素点个数统计的方法计算木粉的几何面积：$A=\sum\limits_{\left(x，y \right)\in Q}{1}$。其中，0为背景，1为木粉，Q为图像D中每个木粉的轮廓区域。
6)利用主轴法求取木粉的最小外接矩形：
① 求使二阶中心距从μ11变得最小的旋转角θ：
| $\tan 2\theta =\frac{2{{\mu }_{11}}}{{{\mu }_{20}}-{{\mu }_{02}}}$ |
| $\begin{align} & {{\mu }_{jk}}=\frac{{{M'}_{jk}}}{{{\left({{M'}_{00}} \right)}^{r}}}，\\ & r=\left(\frac{j+k}{2}+1 \right)，\\ & {{M'}_{jk}}=\sum\limits_{x=1}^{N}{\sum\limits_{Y=1}^{M}{{{\left(x-\bar{x} \right)}^{j}}{{\left(y-\bar{y} \right)}^{k}}f\left(x，y \right)}} \\ \end{align}$ |
② 将x，y轴分别旋转θ角得到木粉的主轴x′，y′；
③ 沿x′方向逐行扫描图像，找出x′方向木粉边界点的最小和最大坐标值p1(xmin，y)和p2(xmax，y)；沿y′方向逐列扫描图像，找出y′方向木粉边界点的最小和最大坐标值p3(x，ymin)和p4(x，ymax)；
④ 计算木粉最小外接矩形的长l和宽b:
| $\begin{align} & l={{x}_{\max }}-{{x}_{\min }}，\\ & b={{y}_{\max }}-{{y}_{\min }}，\\ \end{align}$ |
7)输出木粉长宽比k=l/b，木粉矩形度FA=A/Sr。
2 数学建模与分析通过对数据的分析可知，随粒径变化木粉矩形度基本不变，均值在0.6～0.8之间，如图 4所示。矩形度可反映颗粒的形状特征，分析表明木粉的形状规则度不随粒径变化，而木粉长宽比随木粉粒径的变化有着明显的变化趋势。因此，对木粉粒径和长宽比进行数学建模与分析。
 |
图 4 木粉粒径与矩形度散点 Fig. 4 Scatter of relationship between particle size and squareness |
采用最小二乘法建立木粉粒径与长宽比的拟合方程。给定m组木粉原始数据为{(xi，yj)}(i=0，1，…，m)，其中xi表示木粉i的粒径，yi表示木粉i的长宽比。最小二乘法拟合，就是在取定的函数类Ψ中，求p(x)∈Ψ，使算出的函数曲线与观测值之差的平方和E2最小，即使得原始数据yi与拟合数据p(xi)之差的平方和最小(李明阳等，2013)：
| ${{E}^{2}}=\underset{p\left(x \right)\in \psi }{\mathop{\min }}\sum\limits_{i=1}^{m}{{{\left[ p\left({{x}_{i}} \right)-{{y}_{i}} \right]}^{2}}}$ | (3) |
最小二乘法函数拟合的关键就是使得函数p(x)的几何形状与要拟合的函数分布相近，并且形式简单，便于计算，符合实际应用。而对于曲线拟合的评判，采用回归系数R2、残差平方和SSE检验拟合曲线的可信度：
| ${{R}^{2}}=1-\frac{\text{SST}}{\text{SST}}$ | (4) |
| $\text{SSE=}\sum\limits_{i=1}^{m}{{{\left({{y}_{i}}-{{{\hat{y}}}_{i}} \right)}^{2}}，}\text{SST}=\sum\limits_{i=1}^{m}{{{\left({{y}_{i}}-\bar{y} \right)}^{2}}}$ |
在对木粉粒径与长宽比进行曲线拟合的过程中，尝试比较了多种函数拟合的方法，结果发现，多项式拟合、高斯拟合、傅里叶拟合拟合度都较高。多项式拟合随着拟合次数的增加，高次项的拟合系数会越来越小，引起方程的不稳定；高斯函数可以用较少的拟合项数得到较高的拟合精度；傅里叶变换残差值相对稍高。因此，采用高斯拟合对木粉粒径和长宽比进行拟合和分析。图 5是原始数据散点图及3种拟合曲线图。表 3是分别用3种函数类对原始数据拟合的相关系数和残差值。
 |
图 5 木粉粒径与长宽比数据<>散点及3种拟合曲线 Fig. 5 Scatter of relationship between particle size and aspect ratio and three kinds of fitting curves |
|
|
木粉粒径与长宽比高斯函数模型的表达形式为：
| $\psi \left(x \right)=\sum\limits_{i=1}^{m}{a，{{e}^{-{{\left(\frac{x-{{b}_{i}}}{{{e}_{i}}} \right)}^{2}}}}}$ | (5) |
|
|
为分析拟合曲线的几何特性，根据拟合曲线方程(5)，计算出曲线的二阶导数，将二阶导数绘制成图，见图 6。
 |
图 6 木粉粒径与长宽比高斯拟合曲线二阶导数 Fig. 6 Second derivative variation of gaussian fitting curves of particle size and aspect ratio |
木粉粒径在1 100～30 μm范围内，随粒径的减小，木粉长宽比先增大后减小(图 5)，再结合对照二阶导数(图 6)可以发现：粒径在1 100～800 μm时，长宽比从接近1开始增大，二阶导数取值接近0，此时变化速率较缓慢； 800～660 μm时，长宽比增大速率变大，在750～670 μm范围，长宽比变化速率达到最大；而粒径在576 μm时，长宽比达到最大数值4.6；随后，576～330 μm时，长宽比急剧减小；330～50 μm时，长宽比缓慢减小，数值基本无变化，二阶导数再次接近于0；在粒径小于50 μm范围内，长宽比趋近于1。
木粉长宽比随粒径的这种变化与木材细胞本身的性状有关。本文所用目标木粉成熟管胞的平均长度为563.82 μm，平均宽度为46.498 μm，长宽比的变化与管胞破裂密切相关。木粉加工时，在木粉粒径从大于管胞长度到与管胞长度接近的过程中，随粒径减小，长宽比增大，木粉主要是沿纤维方向断裂而呈现纤维状，木粉粒径与管胞长度接近时，长宽比达到了最大；而木粉粒径从与管胞长度接近开始减小的过程中，木粉主要通过横向断裂，使得木粉管胞沿径向破裂而粒径减小，长宽比减小；当木粉粒径与管胞宽度接近时，单个木粉颗粒已为破裂管胞的一部分，其长宽2个方向的承力能力相近，此时木粉呈现薄细片状，其长宽比不再发生明显变化且数值接近1。
本文建立的木粉粒径与长宽比的高斯拟合曲线和木粉粒径与冲击强度的变化趋势相吻合，这种变化趋势是由于木粉结构和几何特征变化引起的，较大的长宽比有利于增强木粉的力学特性，长宽比是影响其力学性能的根本原因之一。
3 结论本文采用光学显微镜获取单颗粒木粉的粒径和木粉管胞的长宽度，利用数字图像处理技术提取木粉几何特征，通过高斯函数模型对木粉粒径与长宽比进行数学建模，为研究木粉加工过程和提高木塑复合材料力学特性提供理论依据。
1)落叶松木粉粒径为1 100～30 μm时，木粉颗粒矩形度为0.6～0.8，此数值可为木粉粉体群特征性能的研究及木粉颗粒的可视化研究提供支持，为木塑复合材料的力学性能提供理论指导。
2)木粉长宽比随粒径的减小先增大后减小，表明了木粉结构在其加工过程中的变化趋势，特别是细胞壁破裂的情况：当木粉粒径大于管胞长度时，主要通过纵向断裂使粒径减小；当木粉粒径小于管胞长度时，主要通过沿管胞径向的横向断裂而使粒径减小。对木粉的断裂分析研究可为管胞破壁和超细木粉的加工提供科学依据。
3)本文发现的木粉长宽比随粒径的变化趋势与前人研究的木塑复合材料冲击强度随粒径的变化趋势是一致的，长宽比是影响材料力学性能本质因素之一；而获得的最大长宽比对应的粒径，可为提高复合材料的力学特性提供理论指导。
4)本文建立了落叶松木粉粒径与长宽比的高斯数学模型，已知木粉粒径时，可求得其长宽比，并为最终建立长宽比与力学性能之间的定量关系提供参考。
| [1] |
范长胜, 杨冬霞, 杨春梅, 等. 2013. 落叶松锯屑在木粉加工过程中的断裂分析. 南京林业大学学报: 自然科学版, 37(6): 127-131. (Fan C S, Yang D X, Yang C M, et al. 2013. Fracture analysis of larch softwood sawdust in the wood flour processing. Journal of Nanjing Forestry University: Natural Sciences Edition, 37(6): 127-131[in Chinese]).(  1)
|
| [2] |
李兰杰, 刘得志, 陈占勋. 2006. 木粉粒径对木塑复合材料性能的影响. 现代塑料加工应用, 17(5): 21-24. (Li L J,Liu D Z,Chen Z X. 2006. Effects of particle diameter of wood flour on the properties of wood-plastics composite.Modern Plastics Procession and Applications, 17(5): 21-24[in Chinese]).( 1)
|
| [3] |
李明阳, 柏鹏, 王徐华, 等. 2013. 一种基于迭代最小二乘法的精确同步方法. 电子与信息学报, 35(4): 832-837. (Li M Y, Bai P, Wang X H, et al. 2013. A precise synchronization method based on iterative least square algorithm. Journal of Electronics & Information Technology, 35(4): 832-837[in Chinese]).( 1)
|
| [4] |
卢蓉,范勇,陈念年, 等. 2010. 一种提取目标图像最小外接矩形的快速算法.计算机工程, 36(21): 178-180. (Lu R, Fan Y, Chen N N, et al. 2010. Fast algorithm for extracting minimum enclosing ectangle of target image. Computer Engineering, 36(21): 178-180[in Chinese]).( 1)
|
| [5] |
生瑜, 朱德钦, 姚巧艺, 等. 2012.TDI 对聚丙烯/木粉复合材料性能的影响——木粉粒径的影响.福建林学院学报, 32(1): 80-83. (Sheng Y, Zhu D Q, Yao Q Y, et al. 2012. Effects of TDI modification and grain size of wood flour on the PP/wood flour composites. Journal of Fujian College of Forestry, 32(1): 80-83[in Chinese]).( 1)
|
| [6] |
宋丽贤, 张平, 姚妮娜, 等. 2013. 木粉粒径和填量对木塑复合材料力学性能影响研究. 功能材料, 44(17): 2451-2454. (Song L X, Zhang P, Yao N N, et al. 2013. Study on effect of particle diameter and filling quantity of wood flour on mechanical properties of wood-plastics composite. Journal of Functional Materials, 44(17): 2451-2454[in Chinese]).( 1)
|
| [7] |
宋永明, 李春桃, 王伟宏, 等. 2011. 硅烷偶联剂对木粉/HDPE 复合材料力学与吸水性能的影响. 林业科学, 47(6): 122-127. (Song Y M, Li C T, Wang W H, et al. 2011. Coupling effects of silane on the mechanical properties and water absorption of wood flour /HDPE composites. Scientia Silvae Sinicae, 47(6): 122-127[in Chinese]).( 1)
|
| [8] |
宋永明, 王清文. 2012. 木塑复合材料流变行为研究进展. 林业科学, 48(8): 143-149. (Song Y M, Wang Q W. 2012. Research progress on rheological behavior of wood-plastic composites. Scientia Silvae Sinicae, 48(8): 143-149[in Chinese]).( 1)
|
| [9] |
孙永厚, 肖威, 黄美发, 等. 2014. 基于移动最小二乘法的渐开线齿轮齿廓曲线拟合方法. 中国机械工程, 25(22): 3017-3021. (Sun Y H, Xiao W, Huang M F, et al. 2014. Fitting method of tooth profile curve for involute cylindrical gears based on MLS. China Mechanical Engineering, 25(22): 3017-3021[in Chinese]).( 1)
|
| [10] |
王书涛, 张金敏, 李圆圆, 等. 2012. 基于数学形态学和模糊聚类的旋转机械故障诊断. 仪器仪表学报, 33(5): 1055-1061. (Wang S T, Zhang J M, Li Y Y, et al. 2012. Rotating machinery fault diagnosis based on mathematical morphology and fuzzy clustering. Chinese Journal of Scientific Instrument, 33(5): 1055-1061[in Chinese]).( 1)
|
| [11] |
张冰, 杨静, 张健沛, 等. 2014. 面向聚类分析的邻域拓扑势熵数据扰动方法. 哈尔滨工程大学学报, 35(9): 1149-1155. (Zhang B, Yang J, Zhang J P, et al. 2014. A neighborhood topological potential entropy data perturbation method for clustering analysis. Journal of Harbin Engineering University, 35(9): 1149-1155[in Chinese]).( 1)
|
| [12] |
张法全, 王国富, 曾庆宁, 等. 2013. 利用重心原理的图像目标最小外接矩形快速算法. 红外与激光工程, 42(5): 1382-1387. (Zhang F Q, Wang G F, Zeng Q N, et al. New algorithm for minimum enclosing rectangle of the object in the image region based on center-of-gravity principle. Infrared and Laser Engineering, 42(5): 1382-1387[in Chinese]).( 1)
|
| [13] |
朱再胜, 盖国胜, 吴成宝, 等. 2012. 不同含水量木质生物质的粉碎特性研究. 中国粉体技术, 18(2): 7-11. (Zhu Z S, Gai G S, Wu C B, et al. 2012. Pulverizing characteristic of woody biomass with different moisture content. China Powder Science and Technology, 18(2): 7-11[in Chinese]).( 1)
|
| [14] |
Allen G I, Grosenick L, Taylor J. 2014. A generalized least-square matrix decomposition. Journal of the American Statistical Association, 109(505): 145-159.(1)
|
| [15] |
Deng M, Liu Q L, Wang J Q, et al. 2013. A general method of spatio-temporal clustering analysis. Science China Information Sciences, 56(10): 1-14.(1)
|
| [16] |
Guo Q, Chen X, Liu H. 2012. Experimental research on shape and size distribution of biomass particle. Fuel, 94: 551-555.(1)
|
| [17] |
Kobayashi N, Guilin P, Kobayashi J, et al. 2008. A new pulverized biomass utilization technology. Powder technology, 180(3): 272-283.(1)
|
| [18] |
Mok P Y, Huang H Q, Kwok Y L, et al. 2012. A robust adaptive clustering analysis method for automatic identification of clusters. Pattern Recognition, 45(8): 3017-3033.(1)
|
| [19] |
Nourbakhsh A, Karegarfard A, Ashori A, et al. 2010. Effects of particle size and coupling agent concentration on mechanical properties of particulate-filled polymer composites. Journal of Thermoplastic Composite Materials, 23(2): 169-174.(1)
|
| [20] |
Palaniandy S, Kadir N A, Jaafar M. 2009. Value adding limestone to filler grade through an ultra-fine grinding process in jet mill for use in plastic industries. Minerals Engineering, 22(7): 695-703.(1)
|