CN105931225A - 一种基于实时图像检测技术分析晶体生长形状和尺寸分布的方法 - Google Patents

Abstract

本发明提供一种基于实时图像检测技术的在线分析反应釜结晶过程晶体生长形状和尺寸分布方法。首先，针对图像尺寸较大问题，对实时采集的图像进行压缩，并滤波去噪；针对图像光照不均问题，采用图像增强技术，便于图像分割。其次，对分割后的图像的进行尺寸、形状和纹理特征提取，引入内距描述因子区分不同晶体的基本形状。利用聚类算法筛除图像中的颗粒碎片和一些无统计意义的颗粒。然后，采用基于亚像素边缘检测和圆拟合的尺寸标定算法来计算晶体颗粒的尺寸，得到尺寸分布。最后，利用特征降维和支持向量机识别晶体的形状。本发明能够在线提取多幅图像中的晶体颗粒，在线获取晶体生长形状与尺寸分布信息，有助于实时分析反应釜内结晶工况和调控优化。

Description

一种基于实时图像检测技术分析晶体生长形状和尺寸分布的 方法

技术领域

本发明属于工业结晶过程控制技术领域，涉及反应釜结晶过程的实时图像检测技术，具体涉及一种利用高速高分辨率摄像设备实时原位监测结晶反应釜内晶体生长形状与尺寸分布的方法。

背景技术

在化工和制药工程等行业中，需要采用结晶工艺提取产物和提纯产品，结晶过程中各操作参数的控制优劣，对结晶工艺得到的晶体产品质量产生直接影响。实际生产中，通常是对晶体颗粒的生长形状和尺寸分布进行调控优化来保证产品质量和收率。因此，实时获得晶体的形状和尺寸分布是具有重要意义的。目前，实时图像技术是获得晶体形状和二维尺寸的主要检测手段。然而国内外很少有文献给出能应用于工程实践的基于实时图像在线分析反应釜结晶过程晶体生长形状和尺寸分布方法和应用技术。

目前结晶工程实践中，已愈来愈多地采用图像检测技术来分析结晶过程和晶体产品，但是大多数是以离线检测为主，主要是提取一定量的晶体溶液，经过后期处理，利用显微镜观测晶体的形状等信息。这种离线检测方法的主要缺点是不能实时掌握晶体的生长状态，不便于实时调节结晶过程控制参数。因此，如何对反应釜内的晶体生长状态进行实时监测是目前的研究和应用难题。如英国利兹大学化工学院教授王学重在最近研究论文“Stereo imaging of crystal growth”(简译：晶体生长的立体成像，发表在化工领域国际重要刊物AIChE Journal,2016,62(1),18-25)中指出，使用图像可以有效监测结晶生长过程，目前实时图像监测技术在实时性和抗噪辨识方面仍有待提高；国际结晶工程PAT领域专家Z.K.Nagy在近期综述论文“Recent advances in the monitoring,modeling andcontrol of crystallization systems”(简译：结晶系统监测、建模和控制的近期进展，发表在化工领域国际重要刊物Chemical Engineering Research&Design,2013,91(10),1903–1922)中指出，虽然实时图像监测已被广泛认为是对结晶过程控制具有重要意义，然而如何发展和应用这样的技术在线分析反应釜结晶过程晶体生长形状和尺寸分布，以便于实时调控优化结晶过程，存在诸多技术问题和挑战有待于解决。

发明内容

本发明要解决的技术问题是如何利用实时检测的结晶反应釜溶液图像信息来在线分析结晶过程中晶体生长形状与尺寸分布问题。为解决上述这一问题，本发明提出对结晶反应釜内晶体形状与尺寸进行原位测量与统计的技术方法，以实现对结晶过程晶体生长形状和尺寸分布在线检测的效果。

本发明利用非接触式高速高分辨率的图像采集设备，实时采集反应釜内晶体图像，从而分析结晶过程晶体的生长状态。首先，对实时采集的图像进行压缩，并做图像滤波、增强处理，对晶体图像进行分割。其次，对分割后的晶体颗粒图像进行特征提取，把晶体特征分为尺寸特征、形状特征和纹理特征，通过定义内距描述因子(IDD)来区分不同晶体的基本形状。然后，利用聚类算法筛除实时检测图像中的颗粒碎片和粘连颗粒等一些无统计意义的颗粒图像信息，并采用一种基于亚像素边缘检测和圆拟合的尺寸标定算法来计算晶体颗粒的尺寸以及分布。最后，利用特征降维和支持向量机方法来识别晶体的形状。

本发明的技术方案如下：

一种基于实时图像检测技术分析晶体生长形状和尺寸分布的方法，主要针对L-谷氨酸结晶过程，具体包括以下步骤：

(1)处理采集的图像

由于采集的图像具有较大的尺寸，在实时分析中造成时间滞后，因此采用小波变换的图像压缩技术降低图像尺寸。采用导向滤波和冲击滤波结合的方式进行图像滤波。在实时结晶工况下，采集的图像会存在光照不均，背景强度分布不一等现象，这样容易影响后续的分割和提取效果，所以采用图像增强技术，对感兴趣的目标区域进行基于多尺度Retinex算法的增强处理。最后，采用阈值图像分割来提取晶体轮廓。

(2)特征提取

为了计算晶体尺寸和分析晶习，结合晶体控制因素，对晶体的尺寸特征、形状特征和纹理特征进行特征提取，计算晶体尺寸和分析晶习和分析。所述的尺寸特征包括长度、宽度、面积、周长；所述的形状特征包括圆形度、矩形度、长宽比、离心率、凹凸度、傅里叶描述因子、几何矩、内距描述因子；所述的纹理特征包括对比度、能量、熵、分形维数。其中，内距描述因子IDD区分不同晶体的基本形状。

(3)筛选颗粒

晶体结晶过程中，由于晶体运动和过饱和度等因素影响，晶体会出现断裂、破碎、粘连等现象，产生碎片和粘连等无意义的晶体颗粒。人工设定面积阈值方法去除小碎片，尤其在晶体生长初期时容易误删有意义的晶体，按照颗粒的面积进行K-means聚类，自动筛除相对较小的颗粒，这样就避免人工阈值法误删有效晶体的缺陷。而且，按照凸凹度的标准筛除粘连严重的颗粒，根据经验，把凸凹度较小的颗粒划分为粘连晶体。

(4)测量晶体尺寸

在晶体尺寸测量之前，利用圆形的微米尺对图像进行像素当量标定。首先，将微米尺放置在成像点处。其次，为了提高边缘检测定位的精度，采用亚像素边缘检测，对边缘点进行圆拟合，求出圆的半径。最后，求出本成像系统的像素当量。

测量晶体的二维尺寸，采用最小最合适矩形拟合颗粒轮廓，求出该矩形的长轴和短轴，即为晶体的长度和宽度。

(5)识别晶体形状

分为两步：晶体特征降维和晶体形状分类。首先，由于特征集具有非线性，所以采用基于谱回归核判别法(SRKDA)建立晶体特征降维模型对特征集进行降维；其次，采用支持向量机(SVM)对降维后的晶体特征集进行分类。其中，在SVM训练模型参数时，由于SVM的C、g两个参数需预先给出，采用交叉验证法选取这两个参数，选取原则是在识别率最高的情况下，C取最小值时对应g的取值。

本发明的有益效果为：本发明能够在线提取多幅图像中的晶体颗粒，在线获取晶体生长形状与尺寸分布信息，有助于实时分析反应釜内结晶工况和调控优化；能够实现对反应釜内实时结晶图像的快速处理，较为准确地分析结晶过程晶体的生长状态。该方法可操作性强，对经验技术要求较低，能够达到自动快速检测结晶状态效果，便于实际工业应用和推广。

附图说明

图1为本发明的设备装置示意图；

图2为本发明的实时晶体检测过程实现图；

图3为本发明中的图像处理流程图；

图4(a)为原图；

图4(b)为增强图；

图4(c)为分割结果图。

图5(a)为三角形内距描述因子(IDD)说明图；

图5(b)为方形内距描述因子(IDD)说明图；

图5(c)为圆形内距描述因子(IDD)说明图；

图5(d)为六边形内距描述因子(IDD)说明图；

图6为晶体测量结果图；

图7为本发明中的晶体形状识别的过程图；

图中：1温度探头；2循环温控装置；3摄像镜头；4ATR-FTIR探头；5ATR-FTIR计算机；6图像计算机。

具体实施方式

实施例采用4升结晶玻璃反应釜，内置4叶搅拌桨，反应釜内注入了2升的L-谷氨酸水溶液。在结晶反应釜外配置一个非接触式图像采集设备，包括两个高速高分辨率的摄像镜头3。同时，配备了循环温控装置2、温度探头1、ATR-FTIR探头4、ATR-FTIR计算机5和图像计算机6等，该装置详见图1。在结晶过程中，先将结晶溶液温度升到75度，使晶体颗粒充分溶解，再以一定速度降温进行冷却结晶。采集软件配置采用定时采集方式，在采集时间区间每秒获取一幅图像。设定L-谷氨酸的晶体形状分为“α型”、“β型”和“其它”三种形状。在这个冷却结晶过程中，对实时采集的晶体图像进行在线分析和统计结果，全过程如图2所示。

采用本发明的具体实施步骤如下：

第一步，图像处理

晶体图像处理分为图像压缩，图像滤波，图像增强和图像分割，其处理流程如图3所示，其处理结果如图4a-图4c所示。具体过程实现如下：

1.1)图像压缩

设原始晶体图像为I(x,y)，其尺寸为M×N，使用二维离散小波变换进行图像压缩。设行为m，列为n，尺度为j。I(x,y)的小波变换为：

其中，

其中，i＝{H,V,D}分别代表横向、纵向和对角线。那么，原始图像I(x,y)分解成四个部分：低频部分和3个高频部分低频部分保留作为压缩图像。

1.2)图像滤波

采用导向滤波和冲击滤波结合的方式来实现图像滤波过程。首先，使用导向滤波来消除噪声，保持边缘信息。在一个序号为k的窗口ω_k中，输入图像中的点p_i与输出图像中的点q_i关系表示为：

q_i＝a_kp_i+b_k (3)

其中，最小化窗口ω_k的代价函数定义为：

$(a_k,b_k)=\arg \underset{a_k,b_k}{min}\underset{i\∈{\mathrm{\ω}}_k}{\Σ}({\left(a_k{p}_i+b_k-p_i\right)}^2+{\mathrm{\δa}}_k^2)---\left(4\right)$

从而确定(a_k,b_k)的值为

$a_k=\frac{\frac{1}{\left|\mathrm{\ω}\right|}{\mathrm{\Σ}}_{i\∈{\mathrm{\ω}}_k}{p_i}^2-{\mathrm{\μ}}_k{\stackrel{\‾}{p}}_k}{{\mathrm{\σ}}_k^2+\mathrm{\δ}}---\left(5\right)$

和

$b_k={\stackrel{\‾}{p}}_k-a_k{\mathrm{\μ}}_k---\left(6\right)$

其中，

导向滤波输出的结果为：

$q_i=\frac{1}{\left|\mathrm{\ω}\right|}{\mathrm{\Σ}}_{i\∈{\mathrm{\ω}}_k}{a}_k{p}_i+b_k---\left(7\right)$

其次，采用冲击滤波来强化边缘和纹理。设输入图像为q(x,y)，输出图像s(x,y)由下式得出，

$s(x,y)=-\mathrm{sgn}\left(q_{\mathrm{\η}\mathrm{\η}}\right(x,y\left)\right)|\▿q(x,y)|---\left(8\right)$

其中，q_ηη(x,y)为η的二阶导数，为输入图像的梯度图，sgn是符号函数。

1.3)图像增强：

采用图像增强技术解决晶体图像的光照不均问题，其增强图像的结果为r(x,y)，其表达式如下：

$r(x,y)=\underset{k=1}{\overset{3}{\Σ}}{W}_k\{\log s(x,y)-log\[F_k(x,y)*s(x,y)\]\}---\left(9\right)$

其中，*为卷积，s(x,y)是输入图像，k是尺度值，W_k为权重，W₁＝W₂＝W₃＝1/3，k＝1,2,3.F_k(x,y)定义为：

$F_k(x,y)={\mathrm{\λ}}_k{e}^{-(x^2+y^2)/{c_k}^2}---\left(10\right)$

其中，c_k分别为10,80和150，λ_k满足下式

∫∫F_k(x,y)dxdy＝1 (11)

1.4)图像分割：

采用最小交叉熵分割方法来实现晶体图像的颗粒提取。图像一维直方图为h(i)(i＝1,2,...,L)，交叉熵判别函数定义为：

$D\left(t\right)=\underset{i=1}{\overset{t-1}{\Σ}}ih\left(i\right)\ln i/n(1,t)+\underset{i=t}{\overset{L}{\Σ}}ih\left(i\right)\ln i/n(t,L+1)---\left(12\right)$

其中，i是灰度值；t是阈值化时的阈值；n(1,t)和n(t,L+1)是类内均值，分别代表分割后得到的分割图中目标和背景的灰度。

最佳阈值为：

$t^{*}=\arg \underset{t}{min}\left\{D\left(t\right)\right\}---\left(13\right)$

最后，得到二值图像表示为：

$b(x,y)=\left\{\begin{array}{cc}0,& f(x,y)<t^{*}\\ 1,& f(x,y)\&GreaterEqual;t^{*}\end{array}\right.---\left(14\right)$

图像内的颗粒分割后，使用形态学方法填充分割图中颗粒内部的孔洞，同时去除与图像边界相连的不完整颗粒图像。

第二步，特征提取

2.1)晶体尺寸特征

晶体的基本特征直观反映晶体的尺寸。在特征分析中，对晶体成像投影按有效数量平均来计算尺寸参数，采取的特征如下：

长度：晶体投影图像的最佳拟合矩形的长轴像素个数。

宽度：晶体投影图像的最佳拟合矩形的短轴像素个数。

面积：计算图像中颗粒所占像素点的个数。

周长：计算颗粒图像边缘的像素个数。

2.2)晶体形状特征

晶体形状(晶习)是可以通过外观形状直观表现出来，形状特征应该具有旋转、平移和尺度不变性，采用的特征如下：

圆形度e：反映颗粒接近圆的程度。计算公式如下：

e＝4πS/L² (15)

其中，S为颗粒面积，L为颗粒周长。

矩形度E_x：反映颗粒与矩形的相似程度。同时在区域和其最小外接矩形中的像素比例。

E_x＝S/S_e (16)

其中，S_e为最小外接矩形的面积。

长宽比E_r：反应颗粒的伸长情况，是长度P_l与宽度P_m的比值。

E_r＝P_l/P_m (17)

凸凹度S_o：反映晶体形状凸凹程度的一个重要度量。最小凸包采用Graham扫描法来获得。

S_o＝S/S_c (18)

其中，S_c为最小外接多边形的面积。

傅里叶描述因子F：它是物体形状边界曲线的傅里叶变换系数，也即物体边界曲线信号频域分析的结果。傅里叶变换的低频分量对应晶体总体形状，高频分量对应一些晶体形状细节，仅用低频分量可近似描述晶体轮廓形状。

几何矩：由七个Hu不变矩的参数组成。

内距描述因子IDD的计算过程：

定义晶体颗粒中心到边缘的距离为内距。其表示曲线为：

$\left\{\begin{array}{c}x=x\left(\mathrm{\&theta;}\right)\\ y=y\left(\mathrm{\&theta;}\right)\end{array}\right.,-\mathrm{\&pi;}\&le;\mathrm{\&theta;}\&le;\mathrm{\&pi;}---\left(19\right)$

记边缘点为(x_n,y_n)，n＝1,2,…,N，中心(x_c,y_c)可表示为：

$\left\{\begin{array}{c}{x}_c=\frac{1}{N}\underset{n=0}{\overset{N-1}{\&Sigma;}}{x}_n\\ y_c=\frac{1}{N}\underset{n=0}{\overset{N-1}{\&Sigma;}}{y}_n\end{array}\right.---\left(20\right)$

内距计算公式为：

$d_n=\sqrt{{(x_c-x_n)}^2+{(y_c-y_n)}^2}---\left(21\right)$

内距方差表示为：

${\mathrm{\&rho;}}_n=d_n-\stackrel{\&OverBar;}{d}---\left(22\right)$

其中，为均值，那么，

$\stackrel{\&OverBar;}{d}=\frac{1}{N}\underset{n=1}{\overset{N}{\&Sigma;}}{d}_n---\left(23\right)$

内距描述因子(IDD)可算为：

$Bd=S_d/\stackrel{\&OverBar;}{d}---\left(24\right)$

其中，S_d为ρ_n的方差。

使用内距描述因子，对晶体的一些基本形状进行量化标注，如图5a-图5d所示。

2.3)晶体纹理特征

晶体纹理特征是针对晶体表面直观成像的综合描述。纹理特征不是基于像素点的特征，它需要在包含多个像素点的区域中进行统计计算。因此采用灰度共生矩阵(GLCM)描述内部纹理特征。基于灰度共生矩阵的三个纹理描述量：对比度、能量、熵。同时采用分形维参数对图像纹理进行区分。

第三步，颗粒筛选

首先,根据晶体颗粒的面积来筛颗粒碎片。将颗粒面积分为4组，分别为大、中、小和极小。记面积集合为X＝{x(i)|i＝1,2,…,m}，聚类中心为μ(j)，j＝1,…,k。基于K均值聚类的筛选有如下两个步骤：

3.1)把x(i)分配到类别c(i)中，其表达式如下：

$c\left(i\right)=\arg \underset{j}{min}\left|\right|x\left(i\right)-\mathrm{\&mu;}\left(j\right)|{|}^2---\left(25\right)$

3.2)更新X的类中心μ(j)，有

$\mathrm{\&mu;}\left(j\right)=\frac{\underset{i=1}{\overset{m}{\&Sigma;}}1\{c\left(i\right)=j\}x\left(i\right)}{\underset{i=1}{\overset{m}{\&Sigma;}}1\{c\left(i\right)=j\}}---\left(26\right)$

聚类后,对应最小μ(j)的j^*可由下式得出，

$j^{*}=\arg \underset{j}{min}\mathrm{\&mu;}\left(j\right),j=1,...,k---\left(27\right)$

因此，所有在j^*类内的颗粒作为极小颗粒被删除。

同理，计算颗粒的凸凹度，将其分为2类，最小类中心的颗粒划分为粘连颗粒，将被删除。这样，留下的颗粒将用作后续的测量和识别。

第四步，晶体尺寸测量

采用亚像素边缘检测来提取圆的边界，使用最小二乘法拟合圆。实现方法如下：

输入：微米尺图像。

4.1)使用阈值法二值化待检测图像；

4.2)去掉连接边界的连通区域；

4.3)使用灰度矩亚像素边缘方法检测边缘；

4.4)得到边缘点集，用最小二乘法拟合圆；

4.5)求出圆半径的亚像素级值；

4.6)计算亚像素当量。

输出：输出像素当量。

像素当量计算公式为：

P_e＝P_a/P_v (28)

其中，P_a为标定尺的实际尺寸(半径)，P_v为标定尺的像素级尺寸(半径)。

晶体的实际长度L_p和实际宽度W_p计算公式，如下：

$\left\{\begin{array}{c}{L}_p=P_l\&times;P_e\\ W_p=P_m\&times;P_e\end{array}\right.---\left(29\right)$

其中，P_l是像素级长度，P_m是像素级宽度。这样，可以求出颗粒的长度与宽度实际长度，示例测量结果如图6所示。

第五步，晶体形状分类

晶体形状分类的过程如图7所示，包括特征降维和形状分类。经过图像处理和特征提取，随机选择训练样本50个，测试样本100个。

5.1)特征降维

首先，选择傅里叶描述因子个数，对于采集的图像中晶体最小周长小于22的情况，选择最小周长的整数作为描述因子个数，对于最小周长大于等于22的情况，选择描述因子个数为22。其次，采用谱回归核鉴别法建立晶体特征降维模型(SRKDA)，记晶体特征为X＝[x₁,x₂,…,x_l]^T，其核空间映射为φ(x)＝[φ(x₁),φ(x₂),…,φ(x_l)]^T，分类数为c，则

其中，k＝1,…,c and y₀＝[1,1,…,1]^T.

$\left\{\begin{array}{c}Vy=\mathrm{\&lambda;}y\\ \mathrm{\&kappa;}\mathrm{\&alpha;}=y\end{array}\right.---\left(31\right)$

其中，κ(x_i,x_j)＝φ(x_i)^Tφ(x_j)，V为

$V_{ij}=\left\{\begin{array}{cc}1/l_k,& x_{k}and{x}_j\&Subset;classk\\ 0,& otherwise.\end{array}\right.---\left(32\right)$

线性关系为：

(κ+δI)α＝y (33)

其中，δ为调节参数，I为标准矩阵。

α^*＝(κ+δI)^-1y (34)

则

$\mathrm{\&alpha;}=\arg \underset{\mathrm{\&alpha;}}{min}(\underset{i=1}{\overset{N}{\&Sigma;}}{\left(\mathrm{\&kappa;}{\left(x,x_i\right)}^T\mathrm{\&alpha;}-y_i\right)}^2+\mathrm{\&delta;}\underset{i=1}{\overset{N}{\&Sigma;}}|{\mathrm{\&alpha;}}_i\left|\right)---\left(35\right)$

其中，使用线性最小距离法来训练核函数RBF的最优参数。

5.2)形状分类

采用一对一的支持向量机模型来分类，构造3个分类器。在第i类和第j类的区分中，对于训练向量x_t,t＝1,2,…,50，最优分类函数为：

$\begin{array}{c}\underset{w^{ij},{\mathrm{\&xi;}}^{ij},b^{ij}}{\min }\frac{1}{2}{\left(w^{ij}\right)}^T{w}^{ij}+C_p\left(\underset{t}{\&Sigma;}{\left({\mathrm{\&xi;}}^{ij}\right)}_t\right)\\ \left({\left(w^{ij}\right)}^T\mathrm{\&phi;}\right(x_t\left)\right)+b^{ij})\&GreaterEqual;1-{\mathrm{\&xi;}}_t^{ij},if{x}_{t}istheithclass,\\ \left({\left(w^{ij}\right)}^T\mathrm{\&phi;}\right(x_t\left)\right)+b^{ij})\&le;-1+{\mathrm{\&xi;}}_t^{ij},if{x}_{t}isthejthclass,\\ {\mathrm{\&xi;}}_t^{ij}\&GreaterEqual;0.\end{array}---\left(36\right)$

其中，w^ij为权值向量，b^ij为偏项，C_p为惩罚参数。支持向量机的参数由交叉验证法训练得出。针对L-谷氨酸的形状识别，测试样本为100个，其准确识别数为96个，识别率为96％。

Claims (1)

1.一种基于实时图像检测技术分析晶体生长形状和尺寸分布方法，其特征在于，包括以下步骤：

第一步，图像处理

采用小波变换的图像压缩技术降低图像尺寸，防止实时分析中图像较大造成时间滞后；采用导向滤波和冲击滤波结合的方式进行图像滤波；采用图像增强技术，对目标区域进行基于多尺度Retinex算法的增强处理，解决晶体图像光照不均的问题；采用最小交叉熵分割方法提取晶体轮廓；

第二步，结合晶体控制因素，对晶体的尺寸特征、形状特征和纹理特征进行特征提取，计算晶体尺寸和分析晶习；所述的晶体的尺寸特征包括长度、宽度、面积和周长；所述的晶体的形状特征包括圆形度、矩形度、长宽比、离心率、凹凸度、傅里叶描述子、几何矩和内距描述子IDD；所述的晶体的纹理特征包括对比度、能量、熵和分形维数；

其中，内距描述因子IDD区分不同晶体的基本形状，提取晶体的内距描述子IDD过程为：

晶体颗粒中心到边缘的距离称为内距，其表示曲线为：

$\left\{\begin{array}{c}x=x\left(\mathrm{\&theta;}\right)\\ y=y\left(\mathrm{\&theta;}\right)\end{array}\right.,-\mathrm{\&pi;}\&le;\mathrm{\&theta;}\&le;\mathrm{\&pi;}---\left(19\right)$

边缘点为(x_n,y_n)，n＝1,2,…,N，中心(x_c,y_c)为：

$\left\{\begin{array}{c}{x}_c=\frac{1}{N}\underset{n=0}{\overset{N-1}{\&Sigma;}}{x}_n\\ y_c=\frac{1}{N}\underset{n=0}{\overset{N-1}{\&Sigma;}}{y}_n\end{array}\right.---\left(20\right)$

内距计算公式为：

$d_n=\sqrt{{(x_c-x_n)}^2+{(y_c-y_n)}^2}---\left(21\right)$

内距方差为：

${\mathrm{\&rho;}}_n=d_n-\stackrel{\&OverBar;}{d}---\left(22\right)$

其中，为均值，

内距描述因子IDD为：

$Bd=S_d/\stackrel{\&OverBar;}{d}---\left(24\right)$

其中，S_d为ρ_n的方差；

第三步，筛选颗粒

按照颗粒的面积进行K-means聚类，自动筛除相对较小的颗粒；按照凸凹度的标准筛除粘连严重的颗粒，凸凹度较小的颗粒划分为粘连晶体；具体过程为：

将颗粒面积分为大、中、小和极小4组；设面积集合为X＝{x(i)|i＝1,2,…,m}，聚类中心为μ(j),j＝1,…,k，基于K均值聚类的筛选为：

3.1)把x(i)分配到c(i)类别中，表达式为：

$c\left(i\right)=\arg \underset{j}{min}\left|\right|x\left(i\right)-\mathrm{\&mu;}\left(j\right)|{|}^2---\left(25\right)$

3.2)更新X的类中心μ(j)，有

$\mathrm{\&mu;}\left(j\right)=\frac{\underset{i=1}{\overset{m}{\&Sigma;}}1\{c\left(i\right)=j\}x\left(i\right)}{\underset{i=1}{\overset{m}{\&Sigma;}}1\{c\left(i\right)=j\}}---\left(26\right)$

聚类后,对应最小μ(j)的j^*由下式得出，

$j^{*}=\arg \underset{j}{min}\mathrm{\&mu;}\left(j\right),j=1,...,k---\left(27\right)$

所有在j^*类内的颗粒作为极小颗粒被删除；

同理，计算颗粒的凸凹度，将其分为2类，最小类中心的颗粒为粘连颗粒，将被删除；留下的颗粒将用作后续的测量和识别；

第四步，测量晶体尺寸

在测量晶体尺寸之前，用圆形微米尺对图像进行像素当量标定：将微米尺放置在成像点处；为了提高边缘检测定位的精度，采用亚像素边缘检测，对边缘点进行圆拟合，求出圆的半径；求出本成像系统的像素当量；

测量晶体的二维尺寸，采用最小最适矩形拟合颗粒轮廓，求出该矩形的长轴和短轴，即为该颗粒的长度和宽度；

第五步，识别晶体形状

5.1)特征集为非线性，采用基于谱回归核判别法SRKDA的特征降维模型对特征集进行降维；

5.2)采用支持向量机SVM对降维后的特征集进行分类；其中，在SVM训练模型参数时，由于SVM的C、g两个参数需预先给出，采用交叉验证法选取C、g两个参数；选取原则为在模型识别率最高的情况下，C取最小值时对应g的值。