本发明属于机器故障诊断领域的铁谱分析技术，具体涉及离线铁谱图像数字化自动识别技术，特别涉及一种多纹理特征融合的磨粒类型自动识别方法。
背景技术：
机械设备零件失效的主要是由磨损、疲劳和腐蚀造成的，而磨损是造成机械故障和失效的主要原因，80％以上的设备故障是由于磨损所致。磨粒则是摩擦副在发生磨损过程中其接触表面发生相对运动，并与界面之间的介质发生作用的表面磨损的直接产物，其携带了丰富的磨损信息，通过对油液中的磨粒进行分析研究来提取磨损颗粒的数量、形貌、浓度、成分等信息，进而对设备的磨损状态进行定量或定性的分析，找出磨损故障发生的诱导因素，从而预测机械设备的健康状态发展，提高设备运行的可靠性、安全性和可维护性。因此，磨粒分析(wda)技术已经成为设备运行状态监测和故障诊断最有效的方法。
在众多的磨粒分析(wda)技术中，基于磨粒图像分析的铁谱技术由于其分析精度高、准确性好，能够准确获得磨粒的机理与磨损类型，从而提供深层次磨损致因，已经成为了机械系统磨损分析的主要技术依据。在机械设备磨损过程中常见的磨损类型主要有正常磨损、粘着磨损、磨料磨损、疲劳磨损、腐蚀磨损等，与其对应的铁谱分析中的磨粒类型主要为正常磨粒、严重滑动磨粒、切削磨粒、疲劳磨粒、球形磨粒等。不同类型的磨粒具有不同的形态特征，离线铁谱分析技术的核心就在于通过人脑对磨粒的形貌、纹理、颜色、尺寸等特征进行分析，其主要依靠领域专家或有丰富经验的分析员,这要求操作者需要具有很高的专业知识,难以推广应用。由于这种定性分析的结果往往过分依赖于分析员的主观判断及经验,可能会造成很大的误差,而且人工分析效率低,这就会在一定程度上导致人力资源的浪费。因此磨粒的自动识别技术则应运而生，针对离线铁谱图像，上海海事大学提出了一种铁谱磨粒纹理特征提取和模式识别的方法(公开号：cn104484675a)，中国人民解放军总参谋部陆航研究所提出了一种交互式磨粒图像标注方法(公开号：cn102768730a)。随着计算机技术与模式识别方法的不断发展，其不断推动离线铁谱磨粒图像分析朝着自动化、数字化方向发展。
目前在铁谱图像磨粒分割领域应用最广的是灰度阈值法,其简单易行,但由于磨粒一般包括深色、浅色及混合色3大类,一幅铁谱图像中可能同时出现亮度高于背景灰度的亮磨粒和亮度低于背景灰度的暗磨粒,若采用单一固定的阈值则会造成分割不完全,进而影响接下来的磨粒特征提取与识别。清华大学的于世强等在“基于背景色彩识别的磨粒图像分割方法”中首先对彩色磨粒图像采用灰度阈值法分两次对深色磨粒和浅色磨粒进行提取，再用原图像分别减去两次提取的结果得到预处理后的背景图像；然后提取背景图像的红、绿、蓝三套灰度信息所对应的灰度平均值以及标准差，再结合多次试验各颜色分量的误差范围对预处理后的背景图像进行提取，从而得到磨粒图像。该方法需要多次进行阈值分割操作效率低，且其经验应用对象较为固定，针对不同的研究对象以及不同的铁谱显微镜都需要重新进行多次试验来求取误差经验值。海军航空工程学院的郭恒光等在“基于分型特征的磨粒图像分割”中利用磨粒图像的分型维数，多重分型维数，图像的灰度信息，结合自组织特征映射神经网络进行磨粒图像分割。该方法需要提取大量的磨粒图像分型维数特征进行训练神经网络，时间久，效率低，且容易出现孔洞缺陷等问题。此外疲劳磨粒与严重滑动磨粒的尺寸、形状方面特征区分性较小，而表面纹理又具有一定的相似性，如何选取合适的纹理特征则是提高疲劳磨粒与严重滑动磨粒识别准确率的关键问题。南京航空航天大学的吴振锋等在“磨损微粒显微镜形态学特征量化描述体系”中研究了超过100个磨粒参数，包括磨粒形状参数、颜色参数(rgb)、纹理参数和分形参数等。其中基于灰度共生矩阵纹理特征参数在实际应用时针对严重滑动磨粒以及疲劳磨粒区分度较小，导致该两种磨粒自动辨识精度较低。
综上所述，磨粒自动识别存在目标磨粒提取以及磨粒特征提取等困难，从而导致自动识别准确率较低，不能满足工业实际应用需求。因此,将现代计算机技术与铁谱诊断技术相结合,用一种快速、有效的定量分析方法来实现对磨粒的自动识别来满足工业实际应用需求已是刻不容缓。
技术实现要素：
为了解决现有磨粒自动分析的技术瓶颈，本发明的目的在于提供一种多纹理特征融合的磨粒类型自动识别方法，针对磨粒图像中出现的深色、浅色以及混合色的磨粒需要采用多种图像处理技术方法对磨粒进行分割，将磨粒与背景进行分离，提取出目标磨粒。针对疲劳磨粒与严重滑动磨粒的相似性采用基于灰度共生矩阵以及tamura纹理特征参数相融合的方法获取磨粒的纹理特征，最后结合机器学习方法实现磨粒的分类识别，进而为重大装备磨损状态和磨损机理的研究提供参考依据。
为实现上述目的，本发明采用如下技术方案：
一种多纹理特征融合的磨粒类型自动识别方法，包括以下步骤：
步骤一：从标准磨粒图像库中选取出典型的各类数量相同的磨粒图片，作为训练样本；
步骤二：对各类磨粒图像进行处理获得磨粒图像rgb值，再对rgb图像进行灰度化处理，且将原rgb图像转换至hsv颜色空间下的s、v分量所对应的磨粒图像；
步骤三：分别对rgb图像灰度化后的图像以及s、v分量所对应的磨粒图像进行自适应阈值分割以及形态学处理，利用区域标记联通算法求出rgb图像灰度化后的图像、s分量磨粒图像、v分量磨粒图像三幅磨粒图像的磨粒面积，选出其中磨粒面积最大的磨粒图像；
具体步骤如下：
s1，对三幅磨粒灰度图像进行阈值分割，具体步骤如下：
1)记t为磨粒图像前景与背景的分割阈值，其取值范围为0～255,对t依次取0～255值进行计算；
2)根据当前阈值t将图像分割为前景和背景，分别求出前景点数占图像比例为w0，平均灰度为u0；背景点数占图像比例为w1，平均灰度为u1，图像的总平均灰度为u；
根据公式(4)计算前景和背景图象的方差：
公式(4)：δ2(t)＝w0×(u0-u)2+w1×(u1-u)2
3)分别求出t取0～255时所对应的方差δ2(t)，对比256个方差值，选择方差最小时所对应的t值作为磨粒图像分割的阈值；
4)再将灰度化图像按照公式(5)进行二值化运算，得到磨粒为黑色，背景为白色的mask图：
公式(5)：
公式(5)中f(x,y)表示磨粒图像的灰度图；
s2，对得到的灰度图像进行形态学处理，即进行“反向”处理，“闭运算”处理，和“孔洞填充”操作；其中“反向处理”为“闭运算”操作做准备，“闭运算”操作能够连接平滑磨粒图像边缘像数点，“孔洞填充”则针对封闭的磨粒图像内部出现的孔洞进行填充，保证分割出的磨粒图像的完整性；
s3，针对磨粒图像中存在多个磨粒，利用“八连通域算法”找出磨粒面积最大的保留作为研究观察对象，剔除其他较小的磨粒，其算子如图所示，中间点则为当前像素点；具体步骤如下：
1)判断此点八邻域中的最左，左上，最上，上右点的情况；如果都没有点，则表示一个新的区域的开始；
2)如果此点八邻域中的最左有点，上右都有点，则标记此点为这两个中的最小的标记点，并修改大标记为小标记；
3)如果此点八邻域中的左上有点，上右都有点，则标记此点为这两个中的最小的标记点，并修改大标记为小标记；
4)否则按照最左，左上，最上，上右的顺序，标记此点为四个中的一个；
5)标记图像中所有磨粒，并统计各个磨粒面积，保留图像中最大的磨粒；
s4，对比三幅磨粒图像最终所得到保留有最大磨粒的的磨粒图像，选取分割出磨粒面积最大的磨粒图像作为进一步磨粒分析研究对象；
步骤四：对所选出的磨粒图像进行基于概率统计的灰度共生矩阵纹理特征提取，获得与磨粒纹理特征密切相关的4个特征参数，分别为能量、惯性矩、相关性、熵；
步骤五：对所选出的磨粒图像进行基于人类视觉感官的tamura纹理特征提取，获得与磨粒纹理特征密切相关的三个特征参数，分别为粗糙度、对比度、方向度；
步骤六：选取灰度共生矩阵纹理特征及tamura纹理特征利用主成分分析法(pca)求得综合纹理特征参数作为神经网络的输入特征进行训练，构建磨粒的神经网络分类器；
步骤七：选取出待识别的磨粒图像重复步骤二至步骤五，将提取出的待识别的磨粒特征作为训练完毕的神经网络分类器的输入特征向量，实现磨粒的自动分类识别。
所述的步骤二，具体步骤如下：
s1，将磨粒图像中的每个像素点(x,y)的rgb值按照公式(1)进行灰度化处理，得到灰度图像，
公式(1)：f(x,y)＝0.299r(x,y)+0.587g(x,y)+0.114b(x,y)公式(1)中r(x,y)、g(x,y)、b(x,y)代表图像中坐标(x,y)处像素点的r、g、b分量值，下述步骤公式中r(x,y)、g(x,y)、b(x,y)表示意义与公式(1)相同；
s2，将磨粒图像中的每个像素点(x,y)的rgb值按照公式(2)进行颜色空间坐标转换，得到v分量下的磨粒图像，
公式(2)：
s3，将磨粒图像中的每个像素点(x,y)的rgb值按照公式(3)进行颜色空间坐标转换，得到s分量下的磨粒图像，
公式(3)：
所述的步骤六，具体步骤如下：
s1，基于所选磨粒纹理特征参数，构建纹理特征参数矩阵，如公式(20)所示：
公式(19)：x＝[xi1xi2…xij]，(i＝1,2…,n；j＝1,2,…,p)
其中xij是纹理特征参数，i是样本磨粒个数，j是纹理特征参数个数；
s2，计算纹理特征参数相关系数矩阵，计算如公式(20)所示：
公式(20)：
s3，计算纹理特征矩阵的特征值与特征向量，如公式(21)所示：
公式(21)：|λi-r|u＝0
式中i是单位向量，λi(i＝1,2,…,n)，λ1≥λ2…≥λn≥0，u＝[ui1ui2…uin]为特征值λi的特征向量；
s4，计算表面纹理特征参数的贡献率，如公式(22)所示：
公式(22)：
s5，计算磨粒表面纹理综合参数，如公式(23)所示：
公式(23)：zi×k＝xut，(i＝1,2…,n；k＝1,2,…,n)
其中zi×k是与xi×j相对应的主成分特征矩阵，z＝[zi1zi2…zik]即为最终得到的磨粒综合纹理特征向量；
s6，基于所得到的综合纹理特征参数构建磨粒的神经网络分类器其中输入层为z＝[zi1zi2…zik]，在训练前需将z中各个纹理特征参数归一化到[-1,1]，输出层为严重滑动磨粒和疲劳磨粒，输入层与隐层采用tanh激活函数，训练采用动量梯度下降法，其中动量、学习率，学习率增长速度、学习率减少速率分别为0.9、0.01、.05、0.7，训练的终止条件为预测值与期望值的误差小于0.001。
有益效果：
本发明应用于机械设备磨损状态监测领域，具有以下有益效果：
1、本发明将铁谱分析技术与图像处理技术相结合，利用铁谱显微镜获取具有多个磨粒的谱片图像，通过将原rgb图像分别转到hsv颜色空间坐标下的单个分量下的三幅磨粒图像，对所得图像进行分割，自动获取其中分割效果最好的分割后的磨粒图像，实现具有多种情况的磨粒图像分割，适用于铁谱分析领域所有磨粒图像的分割处理；
2、本发明基于图像纹理特征，针对疲劳磨粒和严重滑动磨粒表面的纹理特性，提取了合适的纹理特征参数，为磨粒自动识别提供了准确的表征参数；
3、本文基于机器学习算法，利用主成分分析法(pca)获取了磨粒表面的综合纹理特征，结合神经网络模型进行磨粒的分类识别，提高了疲劳磨粒和严重滑动磨粒的识别精度。
附图说明
图1为一种基于离线铁谱图像磨粒分割和多纹理特征融合的磨粒类型自动识别总体流程图。
图2为部分训练样本磨粒图。
图3(a)为同时具有深色和浅色磨粒的图像。
图3(b)只分割出深色磨粒的图像。
图3(c)深色和浅色磨粒都分割出来的图像。
图3(d)只分割出浅色磨粒的图像。
图4(a)-(d)分别为严重滑动磨粒和疲劳磨粒所提取的相关性、惯性矩、熵、能量特征参数箱体图。
图4(e)-(g)分别为严重滑动磨粒和疲劳磨粒所提取的方向度、对比度、粗糙度特征参数箱体图。
图5(a)-(b)为图像的3×3掩模。
图6(a)为测试磨粒为疲劳磨粒的分类结果示意图。
图6(b)为测试磨粒为严重滑动磨粒的分类结果示意图。
具体实施方式
下面结合附图对本方法进行说明。
参照图1，一种多纹理特征融合的磨粒类型自动识别方法，包括以下步骤：
步骤一、从标准磨粒图像库中选取出典型的各类数量相同的磨粒图片(严重滑动磨粒和疲劳磨粒)，作为训练样本，并标注不同磨粒图片的放大倍数，部分样本如图2所示。
步骤二、由于磨粒图像存在深色、浅色以及混合色三种类型的磨粒，为了提高各种磨粒图像的分割效果，对磨粒图像进行处理获得磨粒图像rgb值，再对rgb图像进行灰度化处理，且将原rgb图像转换至hsv颜色空间下的s、v分量所对应的磨粒图像，将原来一幅彩色磨粒图像转化成三幅灰度化后的磨粒图像，具体步骤如下：
s1，将磨粒图像中的每个像素点(x,y)的rgb值按照公式(1)进行灰度化处理，得到灰度图像，
公式(1)：f(x,y)＝0.299r(x,y)+0.587g(x,y)+0.114b(x,y)公式(1)中r(x,y)、g(x,y)、b(x,y)代表图像中坐标(x,y)处像素点的r、g、b分量值，下述步骤公式中r(x,y)、g(x,y)、b(x,y)表示意义与公式(1)相同；
s2，将磨粒图像中的每个像素点(x,y)的rgb值按照公式(2)进行颜色空间坐标转换，得到v分量下的磨粒图像，
公式(2)：
s3，将磨粒图像中的每个像素点(x,y)的rgb值按照公式(3)进行颜色空间坐标转换，得到s分量下的磨粒图像，
公式(3)：
步骤三：通过颜色空间转化得到s、v分量所对应的磨粒图像，结合rgb图像灰度化后的图像，利用最大类间方差法(otsu)分别对三幅图像进行自适应阈值图像分割，再利用区域标记联通算法分别求出分割后的三幅磨粒图像的磨粒面积，选出其中磨粒面积最大的磨粒图像作为研究对象,分割效果如图3所示，具体步骤如下：
s1，对三幅磨粒灰度图像进行阈值分割，具体步骤如下：
1)记t为磨粒图像前景与背景的分割阈值，其取值范围为0～255,对t依次取0～255值进行计算。
2)根据当前阈值t将图像分割为前景和背景，分别求出前景点数占图像比例为w0，平均灰度为u0；背景点数占图像比例为w1，平均灰度为u1，图像的总平均灰度为u。
3)根据公式(4)计算前景和背景图象的方差：
公式(4)：δ2(t)＝w0×(u0-u)2+w1×(u1-u)2
4)分别求出t取0～255时所对应的方差δ2(t)，对比256个方差值，选择方差最小时所对应的t值作为磨粒图像分割的阈值。
5)再将灰度化图像按照公式(5)进行二值化运算，得到磨粒为黑色，背景为白色的mask图：
公式(5)：
公式(5)中f(x,y)表示磨粒图像的灰度图。
s2，对得到的灰度图像进行形态学处理，即进行“反向”处理，“闭运算”处理，和“孔洞填充”操作。其中“反向处理”为“闭运算”操作做准备，“闭运算”操作能够连接平滑磨粒图像边缘像数点，“孔洞填充”则针对封闭的磨粒图像内部出现的孔洞进行填充，保证分割出的磨粒图像的完整性。
s3，针对磨粒图像中存在多个磨粒，利用“八连通域算法”找出磨粒面积最大的保留作为研究观察对象，剔除其他较小的磨粒，其算子如图所示，中间点则为当前像素点。具体步骤如下：
1)判断此点八邻域中的最左，左上，最上，上右点的情况。如果都没有点，则表示一个新的区域的开始。
2)如果此点八邻域中的最左有点，上右都有点，则标记此点为这两个中的最小的标记点，并修改大标记为小标记。
3)如果此点八邻域中的左上有点，上右都有点，则标记此点为这两个中的最小的标记点，并修改大标记为小标记。
4)否则按照最左，左上，最上，上右的顺序，标记此点为四个中的一个。
5)标记图像中所有磨粒，并统计各个磨粒面积，保留图像中最大的磨粒。
s4，对比三幅磨粒图像最终所得到保留有最大磨粒的的磨粒图像，选取分割出磨粒面积最大的磨粒图像作为进一步磨粒分析研究对象。
步骤四：严重滑动磨粒与疲劳磨粒在形状、尺寸特征参数较为接近，其主要区别在于表面纹理特征，因此对所选出的磨粒图像进行基于概率统计的灰度共生矩阵纹理特征提取，获得与磨粒纹理特征密切相关的4个特征参数，分别为能量、惯性矩、相关性、熵。所提取的各个参数值如图4(a)～(d)所示，各个参数的求取公式分别如下公式(6)～(9)所示：
公式(6)：
其中p(i,j,d,θ)——灰度共生矩阵元素值；(i,j)——灰度共矩阵中元素位置；d为像元距离；θ——纹理方向，一般为0°，45°，90°或135°。
公式(7)：
公式(8)：
式中
公式(9)：
步骤五：tamura纹理特征主要以人类的主观心理度量作为标准，其很好的对应了人类视觉感知，其能够将磨粒表面纹理直观的以特征参数的形式表达出来，有利于疲劳磨粒和严重滑动磨粒表面纹理的区分。本研究主要提取了磨粒表面的粗糙度、对比度和方向度，训练样本所提取的各个参数值如图4(e)～(g)所示，具体参数求取步骤如下：
s1，粗糙度对具有不同结构的纹理模式来说，基元尺寸越大或基元重复次数越小，给人的感觉越粗糙。其计算方法如下：
1)计算图像中大小为2k×2k矩形区域的活动窗口中每个象数点的亮度均值如式(10)所示。其中(x，y)代表所选定区域中的第(i，j)点的亮度值，通过k来确定象数的范围，例如1×1，2×2，3×3...,32×32。
公式(10)：
2)对于每个象数点分别计算它在水平和垂直方向上互不重叠的活动窗口之间的平均强度差。如公式(11)、(12)所示：
公式(11)：ek,h＝|ak(x+2k-1,y)-ak(x-2k-1,y)|
公式(12)：ek,v＝|ak(x,y+2k-1)-ak(x,y-2k-1)|
其中ek,h代表该象数点水平方向差值，ek,v代表该象数点垂直方向差值。根据每个象数点，找到能使e值达到最大的最佳尺寸sbest如式(13)所示，其中，无论任何方向能使e值达到最大的k，如式(14)所示。
公式(13)：sbest(x,y)＝2k
公式(14)：ek＝emax＝max(e1,e2,...,eh)
3)通过计算整幅图像中sbest的平均值得到粗糙度fcrs，其中m和n为图像的宽度和高度，如式(15)所示：
公式(15)：
s2，对比度是指一幅图像中明暗区域最亮的白和最暗的黑之间的亮度层次，象数的差异范围越大代表对比度越大，反之亦然。对每个象数的邻域都分别计算其均值、方差、峰态等统计特征，以衡量整个图像或者区域中对比度的全局变量。其计算如公式(16)所示：
公式(16)：
其中σ是图像灰度的标准方差，σ4表示图像灰度值得峰态，通过α4＝u4/σ4定义；u4是四阶矩阵均值，σ2表示图像灰度值的方差。
s3，方向度描述纹理是如何沿某些方向散步或者集中的，它与纹理基元的形状和排列的规则有关。其计算方法如下：
1)计算每个象素处的梯度向量的模和局部边缘方向。其计算如式(17)，式(18)所示：
公式(17)：|δg|＝(|δh|+|δv|)/2
公式(18)：
这两个公式可以通过图像象素点周围3×3的矩形区域卷积图5中的(a)和(b)的两个3×3掩模实现。
2)把0-π区域划分成16等分，取每个区间的最大值
公式(19)：
其中p代表某个峰值，wp代表谷之间该峰值得范围。该特征代表了整个图像的方向性排列一致度。
步骤六：针对所选的灰度共生矩阵纹理特征及tamura纹理特征利用主成分分析法(pca)求得综合纹理特征参数作为神经网络的输入特征进行训练，构建磨粒的神经网络分类器，具体步骤如下：
s1，基于所选磨粒纹理特征参数，构建纹理特征参数矩阵，如公式(20)所示：
公式(19)：x＝[xi1xi2…xij]，(i＝1,2…,n；j＝1,2,…,p)
其中xij是纹理特征参数，i是样本磨粒个数，j是纹理特征参数个数；
s2，计算纹理特征参数相关系数矩阵，计算如公式(20)所示：
公式(20)：
s3，计算纹理特征矩阵的特征值与特征向量，如公式(21)所示：
公式(21)：|λi-r|u＝0
式中i是单位向量，λi(i＝1,2,…,n)，λ1≥λ2…≥λn≥0，u＝[ui1ui2…uin]为特征值λi的特征向量。
s4，计算表面纹理特征参数的贡献率，如公式(22)所示：
公式(22)：
s5，计算磨粒表面纹理综合参数，如公式(23)所示：
公式(23)：zi×k＝xut，(i＝1,2…,n；k＝1,2,…,n)
其中zi×k是与xi×j相对应的主成分特征矩阵，z＝[zi1zi2…zik]即为最终得到的磨粒综合纹理特征向量。
s6，基于所得到的综合纹理特征参数构建磨粒的神经网络分类器其中输入层为z＝[zi1zi2…zik]，在训练前需将z中各个纹理特征参数归一化到[-1,1]，输出层为严重滑动磨粒和疲劳磨粒，输入层与隐层采用tanh激活函数，训练采用动量梯度下降法，其中动量、学习率，学习率增长速度、学习率减少速率分别为0.9、0.01、.05、0.7，训练的终止条件为预测值与期望值的误差小于0.001。
步骤七：选取出待识别的磨粒图像重复步骤二至步骤五，将提取出的待识别的磨粒特征作为训练完毕的神经网络分类器的输入特征向量，实现磨粒的自动分类识别，其中测试磨粒分类结果示意图如图6所示。