0 引言
作为自动化开采关键技术之一，煤岩识别技术能够有效实现采煤机滚筒截割高度的自动调节，提高煤炭采出率，降低原煤含矸量，同时可避免在截割顶底板岩过程中引发瓦斯爆炸。目前该技术已成为国内外研究热点，已有γ射线探测[1]、雷达探测[2]、应力截齿探测[3]、红外探测[4]、振动检测[5]、声音检测[6]等识别方法，但这些方法普适性较差，无法满足矿井开采实际需求，且精度与可靠性均较低。
根据现有识别技术，结合煤岩图像的纹理方向、周期等特性，使用纹理特征能够较好实现煤、岩图像分类，相关算法的应用正逐渐成为研究热点。目前，在煤岩图像纹理识别方面，局部二值模式(Local Binary Patterns，LBP)[7]以其特征维数低、具有旋转不变性等良好性能正得到广泛研究与关注，发展出旋转不变局部二值模式(Rotation Invariant Local Binary Patterns，RILBP)[8]、统一局部二值模式(Uniform Local Binary Patterns，ULBP)[9]、判别局部二值模式(Discriminative Local Binary Patterns，DisLBP)[10]等。文献[11]提出利用传统LBP算法实现煤岩纹理识别，但该应用只完成了简单的基本算法实现，运行时间及算法效果未依据待提取煤岩图像进行调整，分类准确率较低。文献[12]采用RILBP对煤岩的显微组分图像纹理特征进行识别，算法实现了煤岩纹理的旋转不变识别，但同时也增加了运算复杂度，且分类效果仍较差，无法在实际应用中推广。文献[13]提出了基于CLBP和支持向量诱导字典学习的煤岩识别方法，在一定程度上摆脱了训练样本数量的制约，但识别效果容易受到图像质量的影响，鲁棒性欠佳。同时，以上各类LBP算法均未考虑局部纹理与整体之间的关系，而煤岩图像局部纹理的整体空间分布信息对煤岩分类过程至关重要。
基于煤岩识别应用需求与现有算法存在的缺陷，本文提出了基于变差函数和局部方差图(Variogram-Local Variance Image, V-LVI)的煤岩图像纹理特征提取算法，其在传统LBP算法基础上提出了局部方差图的概念，从局部方差空间分布关系中提取煤岩纹理特征，引入变差函数对传统LBP丢失的纹理信息进行再利用，使得提取的特征更加准确，在保证算法效率的同时，对旋转纹理识别具有很强的鲁棒性,煤岩分类准确率较高。
1 煤岩图像纹理识别基础
煤岩图像识别系统主要依靠煤与岩在图像特征上的差异来对煤与岩进行区分，这些差异是由煤、岩的物理特性不同决定的，包括颜色、光泽、端口、截面纹理等。其中，煤和顶底板岩的纹理特征存在较大差异，是进行煤、岩图像分类的最佳属性。
目前煤炭主要有烟煤、无烟煤，顶底板岩主要有砂岩与泥岩，部分煤岩样本纹理图像如图1所示。
图1 部分煤岩样本纹理图像
Fig.1 Part of texture images of coal-rock samples
作为自然纹理之一，煤岩纹理呈现不规则性，且随机性强，其提取目标应满足特征维数适中、鉴别能力强、稳健性好、计算量小的要求，综合考虑以上因素，本文根据统计法中的LBP思路，提出了基于V-LVI的煤岩图像纹理特征提取算法。
2 基于V-LVI的煤岩图像纹理特征提取算法
V-LVI以LBP为基础，首先在局部二值模式理论框架中逐像素计算局部方差,得到局部方差图，然后使用变差函数从局部方差图中提取纹理特征，最后将所提取的特征与LBP特征进行融合并分类,用于实现煤岩图像纹理识别。
2.1 局部方差图
在LBP的计算过程中，像素LBP值的计算是通过比较邻域点与中心像素灰度值相对大小并赋值“0/1”来实现的，然而，该过程使得LBP算子丢失了描述邻域点与中心像素灰度值相对大小程度的能力。
局部方差是描述数值信息的有效参数，一个像素的局部方差被定义为其LBP计算过程中多个邻域点与中心像素灰度差值的方差，文献[14]在局部二值模式方差(Local Binary Patterns Variance，LBPV)算法中使用了局部方差，然而LBPV存在很多不足，尤其是没有考虑局部方差的空间分布关系。实验发现，不同类型纹理的局部方差空间分布是不同的，从局部方差空间分布关系中提取的特征可以用于纹理识别。
对于以某一像素为中心的圆形区域，设其半径为R，中心像素灰度值为gc，均匀分布在圆周上的n个邻域点的灰度值分别为g0,g1,…,gn-1，该区域纹理T可表示为
T=t(g0-gc,g1-gc,…,gn-1-gc)
(1)
当采用LBP算子描述时，T可通过符号函数进一步表示，如式(2)所示，其中s(x)为符号函数。
T=t(s(g0-gc),s(g1-gc),…,s(gn-1-gc))
(2)
(3)
当使用方差来描述局部纹理时，T可表示为
T=V(K)=V(g0-gc,g1-gc,…,gn-1-gc)
(4)
式中V(K)表示计算集合K中元素的方差。
将V(g0-gc,g1-gc,…,gn-1-gc)赋值给中心像素，那么中心像素值就是局部方差值，它可以用来表示局部纹理信息，像素k处的局部方差表示为var(k)。
对全图像逐像素计算局部方差可以得到局部方差图，因为局部方差图是在LBP的理论框架下计算得到的，所以，也称之为LBP方差图。不同种类煤、岩的纹理图像局部方差图如图2所示。
在局部方差图中，像素值越大，说明在以该像素为中心的局部区域中像素之间的灰度值对比越明显，通常来说，边界的局部方差会很大，而灰度变化较为平缓的区域局部方差很小，局部方差的空间分布反映了纹理的空间分布，这也是V-LVI算法提取纹理特征的理论基础。从图2可看出，煤、岩的局部方差图纹理结构有较大的不同。
2.2 基于变差函数的局部方差空间分布关系描述
局部方差图包含了一定的纹理信息，使用变差函数提取这些纹理信息,并对其进行处理。
变差函数又称为半方差图[15]。变差函数是一个以步长为自变量的函数，在灰度图像某一方向上，一定步长下的像素对灰度之差组成一个集合，这个集合中元素值方差的一半就是该步长下的变差函数值，通过改变步长可以得到多个与步长对应的函数值，通过近似连续的步长计算可以得到变差函数曲线，变差函数的计算如图3所示，图中分别用实线和虚线标示出了步长h为1和2时的像素对组合。
图2 不同煤岩纹理图像及其局部方差图
Fig.2 Coal-rock images with different texture and its local variance images
图3 变差函数计算
Fig.3 Variogram calculation process
设K表示灰度图像沿某一方向顺序分布的像素点集，点集K中像素的个数为q，像素ki表示集合K中的第i个像素，其灰度值的大小为g(ki)，步长为h，那么变差函数r(h)的计算公式为
(5)
式中N(h)表示集合K中相距步长h的像素对个数，N(h)=q-h。
由式(5)可看出，变差函数值为方差的一半。
将式(5)中的灰度值g(ki)替换为局部方差var(ki)，可以得到局部方差图中变差函数的计算公式为
(6)
像素灰度值在空间中的分布具有连续性，即像素灰度具有空间自相关性。研究表明,像素灰度的自相关性随着距离的增大而递减，相近位置像素的灰度值相关性强，而随着距离的增大，这种相关性会变得很弱，以至于可以忽略。LBP算法的邻域半径及统计距离都不应取较大的值，正是受这种空间自相关性的影响。同样地，为了保证纹理信息提取的有效性，变差函数的步长一般不应该超过对应方向像素长度的一半。
在对局部方差图计算变差函数时，方向的选择有无数种，不可能对每个方向计算变差函数。文献[16]认为，纹理图像包含很多显性邻域结构(Domindant Neighborhood Structure，DNS)，利用DNS的差异可以对纹理进行分类，并以此为基础提出了DNS纹理特征提取算法。在对局部方差图计算变差函数时，可以参考DNS算法的思想。为了使最终特征更具有代表性，首先将局部方差图均匀分割成大小为a×a的互不重叠的子图像，然后在每张子图像中选择经过子图像中心位置且在[0°，180°]范围内均匀分布的多个方向计算变差函数。在选择子图像尺寸时，应保证a×a大小的窗口足以涵盖目标纹理的基本信息，DNS算法的相关研究可以作为选择参数a的参考，一般应保证a不小于50像素。
当选择的方向不是轴向时，首先以子图像中心像素为参考,确定该方向对应亚像素及其邻域点位置，然后通过双线性插值计算邻域点的灰度值，进而求得对应亚像素的局部方差，最后在该方向上的亚像素局部方差集合中计算变差函数。设P1,P2,P3,P4为空间中呈正方形角点关系分布的4个相邻像素，P0为正方形区域内或边界上的任一亚像素。P0与正方形边界的距离分别为m1,m2,w1,w2，其在边界上的对应点分别为P0A,P0B，P0C，P0D，如图4所示。
图4 双线性插值
Fig.4 Bilinear interpolation
已知P1，P2，P3，P4的局部方差，那么亚像素P0的局部方差可通过式(7)计算得到。
(7)
对多个子图像对应方向的变差函数值取平均，即可得到目标纹理的变差函数向量，多个方向的变差函数向量组成变差函数矩阵。通过变差函数矩阵计算目标纹理特征的方法有2种：一种是通过变差函数矩阵拟合计算出不同方向变差函数参数，并将其作为纹理特征，这种方法分析效果有限且计算量较大；另一种方法是将变差函数矩阵中的各方向变差函数向量顺序连接得到最终纹理特征，相对于第1种方法，这种方法更为简单有效，因此，本文采用该方法。通过变差函数在局部方差图中提取纹理特征的过程如图5所示，图中选择了0,45,90,135°四个方向计算变差函数，此处的角度指的是其对应方向与x轴正向的夹角，下同。
图5 基于变差函数的纹理特征提取过程
Fig.5 Texture feature extraction process based on variogram
2.3 特征融合及旋转不变性
V-LVI算法提取的煤岩纹理特征本身是对LBP丢失信息的再利用，为了使其在煤岩纹理图像分类时得到更好的效果，需要将该算法提取的特征与LBP相关算法提取的特征进行融合。
特征融合在特征相似度计算阶段进行，首先对2种算法提取的特征分别计算相似度，然后通过融合系数对相似度加权得到最终相似度。LBP相关算法得到的特征为直方图特征，其相似度一般通过直方图距离来描述，比如直方图相交[17]、巴氏距离[18]和卡方距离[19]等，其中卡方距离的描述效果好且应用广泛，其计算公式为
(8)
V-LVI算法提取的特征是变差函数向量，其各步长对应的函数值大小直接反映了纹理信息，采用欧氏距离来描述其相似度。设X和Y分别表示一个训练样本和一个测试样本，X1和Y1表示用LBP算法计算得到的特征向量，X2和Y2表示用V-LVI算法得到的特征向量，那么，特征融合后的最终相似度可以表示为
D=αD1(X1,Y1)+(1-α)D2(X2,Y2)
(9)
式中：D反映了X和Y两种样本的不相似性，D越大，相似度越小；D1(x,y)和D2(x,y)分别为卡方距离和欧氏距离；α为融合系数,0<α<1，通过改变α的大小可以提高分类效果。
在2.2节中，将变差函数矩阵中的变差函数向量顺序连接得到纹理特征，但该纹理特征并不具有旋转不变性，特别是对于各方向变化不均匀的纹理。将变差函数向量依次顺序连接并分别计算相似度可以很好地解决这一问题，为了更清晰地描述，以图6为例进行说明。对于训练样本X和测试样本Y，分别计算出它们4个方向的变差函数向量aX,bX,cX,dX和aY,bY,cY,dY，各向量的方向关系如图6所示，对于训练样本X，从aX开始顺序连接各变差函数向量得到特征向量X1，对于测试样本Y，分别从aY,bY,cY,dY开始顺序连接各变差函数向量得到4个特征向量Y1，Y2，Y3，Y4，这4个特征向量组成旋转特征向量组。分别计算旋转特征向量组中各向量与X1的相似度，并取最小值作为样本Y与X的相似度，即
D(X,Y)=min[(X1,Y1),(X1,Y2),(X1,Y3),(X1,Y4)]
(10)
图6 变差函数向量的方向关系
Fig.6 Directional relationship of variogram vectors
变差函数的计算方向越多，V-LVI算法对旋转纹理的描述能力越强，但这会相应地增加计算量。
2.4 算法流程
对V-LVI算法提取的特征与LBP特征进行分类。V-LVI算法流程如图7所示。分类的过程实现了特征融合，并使纹理特征具有旋转不变性。
图7 V-LVI算法流程
Fig.7 Flow of V-LVI algorithm
3 实验与结果分析
为验证V-LVI算法在煤岩纹理识别中应用的有效性，同时为使算法研究具有代表性，本文选择烟煤与无烟煤作为煤样本，砂岩与泥岩作为岩样本，并由此组成煤岩纹理图像样本数据库，进行纹理分类测试。
在进行特征提取前，对样本进行随机角度的旋转，以验证旋转不变性。对5种经典LBP相关纹理特征提取算法进行对比实验。参与对比的算法包括RILBP、ULBP、LBPV、DisLBP、多尺度局部二值模式共生矩阵(Multi-Scale Local Pattern Co-occurrence Matrix,MS-LPCM)[20]。其中MS-LPCM算法是一种融合LBP算法和共生矩阵思想的纹理特征提取算法。
以上5种算法提取的特征为直方图特征，本文选择卡方距离作为所提取特征的相似性度量。K近邻(K-Nearest Neighbor，KNN)分类器具有简单有效、实用性强的特点，选其作为本实验的分类器。各种纹理特征提取算法的分类效果通过分类准确率来反映。
3.1 煤岩纹理图像样本数据集预处理
实验采集的原始样本包含112张纹理图像，其中样本纹理类型为无烟煤纹理、长焰煤纹理、砂岩纹理和泥岩纹理。
在样本数据集中随机选取50张图像，对每张图像做4次随机角度的旋转，得到4张旋转图像，使用大小为160×160的窗口，从原图像及其旋转图像中截取子图像，保证来自同一张图像的任意2张子图像互不重叠。每个子图像作为一个样本，每张纹理图像可以得到36个样本，最终得到1 800个样本，这1 800个样本组成了本实验的样本集合。从该样本集合的每个类别中随机选取100个样本加入训练集，剩下的样本加入测试集，最终得到的训练集包含400个样本，测试集包含1 400个样本。
3.2 算法参数设置与实验结果分析
计算得到LBP的半径和邻域点个数组合(R,n)分别为(3,8)、(4,8)，DisLBP算法中的模式频数阈值(Pattern Occurrence)设置为90%，MS-LPCM算法的模式对统计距离d=R，V-LVI算法的窗口大小为51×51，变差函数计算方向包括0，45，90，135°四个方向，对V-LVI特征与RILBP特征进行融合。对于不同算法提取的特征使用KNN分类器进行分类，重复做10次实验，取10次分类准确率的平均值作为最终分类准确率，分类结果见表1。图8为R=3,n=8时，不同算法的煤岩图像纹理分类准确率柱状图，图9为R=4,n=8时，不同算法的煤岩图像纹理分类准确率柱状图。
表1 不同算法的煤岩图像纹理分类准确率
Table 1 Classification precision of texture of coal-rock images using different algorithms %
从分类结果可以看出，V-LVI算法的分类准确率最高，达到了86.02%，这是因为V-LVI算法对LBP特征信息和局部方差图中的方差分布信息进行了融合。在其他算法中，MS-LPCM算法比较充分地利用了局部信息和模式分布信息，其分类准确率也较高，但是MS-LPCM算法对距离近似为R的模式对也进行了统计，因此,所提取的特征不够准确；ULBP算法的分类准确率要低于RILBP算法,说明在该数据集中非统一模式同样包含了对分类有用的信息；LBPV算法将反映局部纹理信息的方差融入到LBP统计中，所以,其分类准确率较传统LBP算法有所提高；虽然DisLBP算法中融合了费舍尔判别思想，但其对类别信息的利用只是停留在局部模式数量的层次，这种对类别信息的利用方式并不理想，有时甚至会丢失关键信息，因此,相对于其他算法其分类准确率并不高。
图8R=3,n=8时不同算法的煤岩图像纹理分类准确率柱状图
Fig.8 The histogram of classification precision of texture of coal-rock images using different algorithms under condition ofR=3,n=8
图9R=4,n=8时不同算法的煤岩图像纹理分类准确率柱状图
Fig.9 The histogram of classification precision of texture of coal-rock images using different algorithms under condition ofR=4,n=8
4 结论
基于变差函数和局部方差图的煤岩纹理特征提取算法依据煤岩纹理识别需求，结合变差函数和局部方差图获得煤岩纹理特征，解决了现有LBP算法的局部二值化纹理丢失问题，实现了对LBP丢失信息的再利用，增强了煤岩图像纹理描述能力，实现了煤岩图像纹理识别的高鲁棒旋转不变性，保证了算法的高效率与实时性。
实验结果表明，该算法能够稳定准确地提取纹理特征，其对纹理的提取结果优于RILBP，ULBP，LBPV，DisLBP，MS-LPCM等相关算法，具有较高的煤岩分类准确率，达到了86.02%。
参考文献(References)：
[1] 张宁波.综放开采煤矸自然射线辐射规律及识别研究[D].徐州：中国矿业大学,2015.
[2] 李亮，王昕，胡克想，等.探地雷达探测煤岩界面的方法与试验[J].工矿自动化,2015,41(9):8-11.
LI Liang,WANG Xin,HU Kexiang,et al.Coal-rock interface detection method using ground penetrating radar and its experiment[J].Industry and Mine Automation,2015,41(9):8-11.
[3] 李文华，刘娇，柴博.采煤机截齿与牵引速度和滚筒转速之间关系[J].测控技术，2017，36(1):133-136.
LI Wenhua,LIU Jiao,CHAI Bo.Relationship between coal winning machine cutting teeth and traction speed and drum rotational speed[J].Measurement & Control Technology,2017,36(1):133-136.
[4] 张强,王海舰,王兆,等.基于红外热像检测的截齿煤岩截割特性与闪温分析[J].传感技术学报,2016,29(5):686-692.
ZHANG Qiang,WANG Haijian,WANG Zhao,et al.Analysis of coal-rock's cutting characteristics and flash temperature for peak based on infrared thermal image testing[J].Chinese Journal of Sensors and Actuators,2016,29(5):686-692.
[5] 王水生.基于振动特性分析的采煤机煤岩识别控制系统[J].工矿自动化,2015,41(5):83-87.
WANG Shuisheng.Coal-rock recognition and control system of shearer based on vibration characteristics analysis[J].Industry and Mine Automation,2015,41(5):83-87.
[6] 李力，欧阳春平.基于超声相控阵的煤岩界面识别研究[J].中国矿业大学学报，2017,46(3):485-492.
LI Li,OUYANG Chunping.Research on coal-rock interface recognition based on ultrasonic phased array[J].Journal of China University of Mining & Technology,2017,46(3):485-492.
[7] OJALA T,PIETIKAINEN M,HARWOOD D.A comparative study of texture measures with classification based on featured distributions[J].Pattern Recognition,1996,29(1):51-59.
[8] OJALA T,PIETIKAINEN M,MAENPAA T.Multiresolution gray-scale and rotation invariant texture classification with local binary patterns[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2002,24(7):971-987.
[9] HAFIANE A,SEETHARAMAN G,ZAVIDOVIQUE B.Median binary pattern for textures classification[C]// International Conference on Image Analysis and Recognition,2007:387-398.
[10] GUO Yimo,ZHAO Guoying,PIETIKAINEN M.Discriminative features for texture description[J].Pattern Recognition,2012,45(10):3834-3843.
[11] 孙继平,杨坤.一种基于图像LBP的煤岩识别方法：CN105184300A[P].2015-12-23.
[12] 王素婷，朱宪坤，吕青.基于RILBP-GLCM算法的煤岩显微组分识别[J].煤炭技术,2017,36(3):142-144.
WANG Suting,ZHU Xiankun,LYU Qing.Coal rock macerals recognition based on RILBP-GLCM algorithm[J].Coal Technology,2017,36(3):142-144.
[13] 孙继平，陈滨.基于CLBP和支持向量诱导字典学习的煤岩识别方法[J].煤炭学报,2017,42(12):3338-3348.
SUN Jiping,CHEN Bin.Coal-rock recognition approach based on CLBP and support vector guided dictionary learning[J].Journal of China Coal Society,2017,42(12):3338-3348.
[14] GUO Zhenhua,ZHANG Lei,ZHANG David.Rotation invariant texture classification using LBP variance (LBPV) with global matching[J].Pattern Recognition,2010,43(3):706-719.
[15] MATHERON G.Principles of geostatistics[J].Economic Geology,1963,58(8):1246-1266.
[16] KHELLAH F M.Texture classification using dominant neighborhood structure[J].IEEE Transactions on Image Processing,2011,20(11):3270-3279.
[17] 李盼.融入直方图相交核的局部稀疏编码图像分类算法研究[D].哈尔滨：哈尔滨工业大学,2014.
[18] 宣国荣，柴佩琪.基于巴氏距离的特征选择[J].模式识别与人工智能,1996(4):324-329.
XUAN Guorong,CHAI Peiqi.Feature Selection based on Bhattacharyya distance[J].Pattern Recognition and Artificial Intelligence,1996(4):324-329.
[19] 谢红，赵洪野.基于卡方距离度量的改进KNN算法[J].应用科技,2015(1):10-14.
XIE Hong,ZHAO Hongye.An improved KNN algorithm based on Chi-square distance measure[J].Applied Science and Technology,2015(1):10-14.
[20] SUN Xiangping,WANG Jin,CHEN Ronghua,et al.Multi-scale local pattern co-occurrence matrix for textural image classification[C]//The 2012 IEEE International Joint Conference on Neural Networks,2012:1-7.