CN104966074A - 基于变换域广义伽玛分布的煤岩识别方法 - Google Patents

Abstract

本发明公开了基于变换域广义伽玛分布的煤岩识别方法。在样本训练阶段，获取满足一定条件的多幅煤炭训练样本图像和多幅岩石训练样本图像，对每幅图像进行两路小波变换，经过相关处理得到按一定排序规则排列的均方根小波系数子带，计算这些均方根小波系数子带在服从广义伽玛分布条件下的各参数，并利用这些参数构造训练样本特征列向量；在煤岩识别阶段，采集满足一定条件的未知类别测试样本图像，与训练样本图像的处理方式类似，对测试样本图像进行一系列处理，最终得到测试样本特征列向量。通过比较测试样本特征列向量与训练样本特征列向量之间的相似度判断测试样本所属的煤岩类型。本发明可靠性好，正确识别率高，软硬件实现容易。

Description

基于变换域广义伽玛分布的煤岩识别方法

技术领域

本发明涉及基于变换域广义伽玛分布的煤岩识别方法，属于图像识别技术领域。

背景技术

煤岩识别是指通过各种技术手段自动判别出煤炭与岩石。在煤炭资源开采及运输过程中，存在许多生产环节需要判别区分煤炭与岩石，如采煤机滚筒高度调节、综采放顶煤过程控制、选煤厂原煤选矸等。从20世纪50年代开始，南非、澳大利亚、德国、美国、中国等世界主要产煤国家对煤岩识别方法展开了一系列研究，相继产生了一些代表性的研究成果，如自然γ射线探测法、雷达探测法、红外探测法、有功功率检测法、振动信号检测法、声音信号检测法等。然而这些方法均存在以下共性问题：(1)需要在现有设备上安装部署各种传感器、导致装置结构复杂，制造成本高；(2)采煤机、掘进机等机械设备在煤炭生产过程中受力复杂、振动剧烈、磨损严重，传感器部署相对比较困难，其电子线路也容易受到损坏，装置可靠性差；(3)针对不同类型的机械载体设备，传感器的选型和安装位置的选择存在较大区别，这就需要进行个性化定制，因此其普适性不佳。

通过对块状的煤炭、岩石样本的观察，发现煤炭与岩石在颜色、光泽、纹理等方面存在较大的差异。当通过现有的数字摄像设备对煤炭或岩石进行成像时，煤炭与岩石的视觉差异信息就隐藏在采集得到的数字图像中了，因此提出通过挖掘煤岩数字图像中的视觉信息来区分煤炭与岩石。现有的基于图像处理的煤岩识别方法在鲁棒性、识别率等方面还存在着较大的提升空间。

发明内容

为了克服现有煤岩识别方法存在的不足，本发明提出基于变换域广义伽玛分布的煤岩识别方法，该方法具有实时性强、识别率高、稳健性好等优点，能够为现代化煤矿安全高效生产提供有力的技术保障。

本发明所述的煤岩识别方法采用如下技术方案实现，包括样本训练阶段和煤岩识别阶段，具体步骤如下：

A.在样本训练阶段，获取相同光照条件下拍摄的m幅已知煤炭样本图像和m幅已知岩石样本图像，截取宽度和高度均为偶数个像素点并且不含非煤岩背景的子图像，煤炭样本子图和岩石样本子图分别记为c₁，c₂，…，c_m和r₁，r₂，…，r_m；

B.分别对样本子图c₁，c₂，…，c_m和r₁，r₂，…，r_m进行N级两路小波变换，每幅样本子图得到12×N个高频部分小波系数子带，每个高频部分小波系数子带记为 表示样本子图在经过第i级两路小波变换以后的第d个方向的第j路高频部分小波系数子带，其中i表示两路小波变换的级数序号，i的取值为1，2，…，N，d表示两路小波变换的方向序号，d的取值为1，2，…，6，j表示两路小波变换的路序号，j的取值为1，2；

C.对步骤B所述的每幅样本子图的级数序号相同并且方向序号也相同的两个高频部分小波系数子带进行求均方根操作，得到6×N个均方根小波系数子带，每一个均方根小波系数子带记为RmsCoef_i-d，RmsCoef_i-d表示样本子图在经过第i级两路小波变换以后的第d个方向的均方根小波系数子带；

D.分别计算步骤C所述的每一个均方根小波系数子带RmsCoef_i-d的均值μ_i-d和方差

E.根据步骤D所计算的每一个均方根小波系数子带的均值和方差，对步骤C所述的每一级两路小波变换以后的6个方向的均方根小波系数子带按照均值和方差的乘积值从大到小的顺序排列，每一个经过排序以后的均方根小波系数子带记为RmsCoefS_i-n，RmsCoefS_i-n表示样本子图在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带，n表示均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排列的排列序号，n的取值为1，2，…，6；

F.分别计算步骤E所述的经过排序以后的每一个均方根小波系数子带在服从概率密度函数为的广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ，其中x表示均方根小波系数子带中的元素，也就是均方根小波系数，函数在Γ(α)函数表达式中，α为自变量，t为积分变量，积分区间为[0，+∞)，e表示自然常数；

G.利用每一幅煤炭样本子图c_k或岩石样本子图r_k经过N级两路小波变换、高频部分小波系数子带求均方根操作并且按上述步骤E排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ，构造出一个6×N×3维的特征列向量或其中下标k表示样本序号，k的取值为1，2，…，m，然后把这2×m个特征列向量保存到数据存储设备中，用于后续的煤岩识别阶段；

H.在煤岩识别阶段，获取相同光照条件下拍摄的未知类别煤岩图像，截取宽度和高度均为偶数个像素点并且不含非煤岩背景的待识别子图u_x；

I.与上述步骤B类似，对u_x进行N级两路小波变换，从而得到12×N个高频部分小波系数子带，每个高频部分小波系数子带记为 表示待识别子图u_x在经过第i级两路小波变换以后的第d个方向的第j路高频部分小波系数子带；

J.与上述步骤C类似，对步骤I所述的待识别子图u_x的级数序号相同并且方向序号也相同的两个高频部分小波系数子带进行求均方根操作，得到6×N个均方根小波系数子带，每一个均方根小波系数子带记为RmsUxCf_i-d，RmsUxCf_i-d表示待识别子图u_x在经过第i级两路小波变换以后的第d个方向的均方根小波系数子带；

K.与上述步骤D类似，分别计算步骤J所述的每一个均方根小波系数子带RmsUxCf_i-d的均值和方差

L.与上述步骤E类似，根据步骤K所计算的每一个均方根小波系数子带的均值和方差，对步骤J所述的每一级两路小波变换以后的6个方向的均方根小波系数子带按照均值和方差的乘积值从大到小的顺序排列，每一个经过排序以后的均方根小波系数子带记为RmsUxCfS_i-n，RmsUxCfS_i-n表示待识别子图u_x在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带；

M.与上述步骤F类似，分别计算步骤L所述的经过排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ；

N.与上述步骤G类似，利用待识别子图u_x经过N级两路小波变换、高频部分小波系数子带求均方根操作并且按上述步骤L排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ，构造出一个6×N×3维的特征列向量

O.比较特征列向量与特征列向量之间的相似度，其中下标k表示样本序号，k的取值为1，2，…，m，从而判定待识别子图u_x所属的煤岩类别。

所述N级两路小波变换包括以下步骤：

(1)初始化设置两路小波变换的级数序号i＝1；

(2)彩色图像灰度化处理，即把24位煤岩彩色图像转换成8位煤岩灰度图像，采用的转换公式为Y＝0.299R+0.587G+0.114B，其中R，G和B分别表示转换前彩色图像的红色，绿色和蓝色分量，Y表示转换后灰度图像的像素值；

(3)用一维的低通滤波器h₀对煤岩灰度图像进行逐列滤波，从而得到低频部分小波系数子带Coef_Lo；

(4)用一维的高通滤波器h₁对煤岩灰度图像进行逐列滤波，从而得到高频部分小波系数子带Coef_Hi；

(5)用步骤(3)所述的低通滤波器h₀对步骤(4)所述的高频部分小波系数子带Coef_Hi进行逐行滤波，从而得到高频部分小波系数子带Coef_HL；

(6)用步骤(4)所述的高通滤波器h₁对步骤(3)所述的低频部分小波系数子带Coef_Lo进行逐行滤波，从而得到高频部分小波系数子带Coef_LH；

(7)用步骤(4)所述的高通滤波器h₁对步骤(4)所述的高频部分小波系数子带Coef_Hi进行逐行滤波，从而得到高频部分小波系数子带Coef_HH；

(8)用步骤(3)所述的低通滤波器h₀对步骤(3)所述的低频部分小波系数子带Coef_Lo进行逐行滤波，从而得到低频部分小波系数子带Coef_LL；

(9)利用步骤(5)所述的高频部分小波系数子带Coef_HL，分别构造煤岩图像在经过第1级两路小波变换以后的第1个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第1个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第6个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第6个方向的第2路高频部分小波系数子带

(10)利用步骤(6)所述的高频部分小波系数子带Coef_LH，分别构造煤岩图像在经过第1级两路小波变换以后的第3个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第3个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第4个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第4个方向的第2路高频部分小波系数子带

(11)利用步骤(7)所述的高频部分小波系数子带Coef_HH，分别构造煤岩图像在经过第1级两路小波变换以后的第2个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第2个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第5个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第5个方向的第2路高频部分小波系数子带

(12)两路小波变换的级数序号自增1，即i＝i+1；

(13)动态申请与Coef_LL数据结构完全一致的临时数据存储区TempCoef，并将其数据内容初始化为Coef_LL的数据内容；

(14)判断i≤N条件是否成立，如果是，则进入以下步骤(15)-(28)的迭代循环，如果否，则转到步骤(29)；

(15)如果TempCoef数据的总行数不是4的倍数，那么对TempCoef数据进行修正，在TempCoef数据的第一行Row_first前面追加一行Row′_first并且Row′_first的数据内容用Row_first的数据内容进行填充，在TempCoef数据的最后一行Row_last后面追加一行Row′_last并且Row′_last的数据内容用Row_last的数据内容进行填充；

(16)如果TempCoef数据的总列数不是4的倍数，那么对TempCoef数据进行修正，在TempCoef数据的第一列Col_first前面追加一列Col′_first并且Col′_first的数据内容用Col_first的数据内容进行填充，在TempCoef数据的最后一列Col_last后面追加一列Col′_last并且Col′_last的数据内容用Col_last的数据内容进行填充；

(17)用两组一维的低通滤波器h₀₀、h₀₁对TempCoef数据进行双重列滤波，从而得到低频部分小波系数子带Coef_L2-i，其中下标i表示两路小波变换的级数序号；

(18)用两组一维的高通滤波器h₁₀、h₁₁对TempCoef数据进行双重列滤波，从而得到高频部分小波系数子带Coef_H2-i，其中下标i表示两路小波变换的级数序号；

(19)用步骤(17)所述的低通滤波器h₀₀、h₀₁对步骤(18)所述的高频部分小波系数子带Coef_H2-i进行双重行滤波，从而得到高频部分小波系数子带Coef_HL2-i，其中下标i表示两路小波变换的级数序号；

(20)用步骤(18)所述的高通滤波器h₁₀、h₁₁对步骤(17)所述的低频部分小波系数子带Coef_L2-i进行双重行滤波，从而得到高频部分小波系数子带Coef_LH2-i，其中下标i表示两路小波变换的级数序号；

(21)用步骤(18)所述的高通滤波器h₁₀、h₁₁对步骤(18)所述的高频部分小波系数子带Coef_H2-i进行双重行滤波，从而得到高频部分小波系数子带Coef_HH2-i，其中下标i表示两路小波变换的级数序号；

(22)用步骤(17)所述的低通滤波器h₀₀、h₀₁对步骤(17)所述的低频部分小波系数子带Coef_L2-i进行双重行滤波，从而得到低频部分小波系数子带Coef_LL2-i，其中下标i表示两路小波变换的级数序号；

(23)利用步骤(19)所述的高频部分小波系数子带Coef_HL2-i，分别构造煤岩图像在经过第i级两路小波变换以后的第1个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第1个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第6个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第6个方向的第2路高频部分小波系数子带

(24)利用步骤(20)所述的高频部分小波系数子带Coef_LH2-i，分别构造煤岩图像在经过第i级两路小波变换以后的第3个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第3个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第4个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第4个方向的第2路高频部分小波系数子带

(25)利用步骤(21)所述的高频部分小波系数子带Coef_HH2-i，分别构造煤岩图像在经过第i级两路小波变换以后的第2个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第2个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第5个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第5个方向的第2路高频部分小波系数子带

(26)释放临时数据存储区TempCoef所占用的内存空间，重新动态申请与步骤(22)所述Coef_LL2-i的数据结构完全一致的临时数据存储区TempCoef，并将其数据内容初始化为Coef_LL2-i的数据内容；

(27)两路小波变换的级数序号自增1，即i＝i+1；

(28)转到步骤(14)，进行下一级两路小波变换；

(29)释放临时数据存储区TempCoef所占用的内存空间，完成N级两路小波变换。

所述一维的低通滤波器h₀和高通滤波器h₁的序列长度均为奇数，并且h₀序列和h₁序列均为中心对称序列。

所述一维的低通滤波器h₀₀、h₀₁和高通滤波器h₁₀、h₁₁的序列长度相等并且均为偶数，h₀₁序列为h₀₀序列的反转序列，h₁₀序列通过对h₀₀序列进行反转操作然后对偶数序号采样点的序列元素进行符号取反操作而得到，h₁₁序列为h₁₀序列的反转序列。

在给定两组序列长度相等的一维滤波器h_f、h_s和滤波对象X的前提下，所述双重列滤波包括以下步骤：

STC1.对X进行行对称延拓得到X_ex，如果X的总行数和总列数分别为RowsX和ColsX，h_f和h_s的序列长度为len，那么X_ex的总行数和总列数分别为(RowsX+2×len)和ColsX，并且X_ex的第1，2，…，len行的数据内容分别用X的第len，len-1，…，1行的数据内容进行填充，X_ex的第len+1，len+2，…，len+RowsX行的数据内容分别用X的第1，2，…，RowsX行的数据内容进行填充，X_ex的第len+RowsX+1，len+RowsX+2，…，len+RowsX+len行的数据内容分别用X的第RowsX，RowsX-1，…，RowsX-(len-1)行的数据内容进行填充；

STC2.抽取h_f中序号为奇数的滤波系数构成子滤波器h_fo，抽取h_f中序号为偶数的滤波系数构成子滤波器h_fe，抽取h_s中序号为奇数的滤波系数构成子滤波器h_so，抽取h_s中序号为偶数的滤波系数构成子滤波器h_se；

STC3.对步骤STC1所述的X_ex进行修正，即舍弃X_ex的起始两行和末尾两行，从而得到X′_ex；

STC4.抽取步骤STC3所述X′_ex中行号为奇数的数据构成X′_exO，抽取步骤STC3所述X′_ex中行号为偶数的数据构成X′_exE；

STC5.抽取步骤STC4所述X′_exE中行号为偶数的数据构成X′_exA，抽取步骤STC4所述X′_exO中行号为偶数的数据构成X′_exB，抽取步骤STC4所述X′_exE中行号为奇数的数据构成X′_exC，抽取步骤STC4所述X′_exO中行号为奇数的数据构成X′_exD；

STC6.用步骤STC2所述的子滤波器h_fo对步骤STC5所述的X′_exB进行逐列滤波，从而得到ResB；

STC7.用步骤STC2所述的子滤波器h_fe对步骤STC5所述的X′_exD进行逐列滤波，从而得到ResD；

STC8.用步骤STC2所述的子滤波器h_so对步骤STC5所述的X′_exA进行逐列滤波，从而得到ResA；

STC9.用步骤STC2所述的子滤波器h_se对步骤STC5所述的X′_exC进行逐列滤波，从而得到ResC；

STC10.把步骤STC6所述的ResB和步骤STC7所述的ResD中行号相同且列号也相同的元素相加得到ResBD，即ResBD＝ResB+ResD；

STC11.把步骤STC8所述的ResA和步骤STC9所述的ResC中行号相同且列号也相同的元素相加得到ResAC，即ResAC＝ResA+ResC；

STC12.申请一块总行数为(RowsAC+RowsBD)且总列数为ColsAC或ColsBD的数据存储区FinalRes，其中RowsAC表示步骤STC11所述ResAC的总行数，RowsBD表示步骤STC10所述ResBD的总行数，ColsAC表示步骤STC11所述ResAC的总列数，ColsBD表示步骤STC10所述ResBD的总列数，并且满足RowsAC＝RowsBD和ColsAC＝ColsBD的条件；

STC13.把一维滤波器h_f和h_s视为两个向量，如果这两个向量的内积大于0，那么步骤STC12所述FinalRes的第1，3，5，7，…，(RowsAC+RowsBD-1)行的数据内容分别用步骤STC10所述ResBD的第1，2，3，4，…，RowsBD行的数据内容进行填充，步骤STC12所述FinalRes的第2，4，6，8，…，(RowsAC+RowsBD)行的数据内容分别用步骤STC11所述ResAC的第1，2，3，4，…，RowsAC行的数据内容进行填充，如果这两个向量的内积小于或等于0，那么步骤STC12所述FinalRes的第1，3，5，7，…，(RowsAC+RowsBD-1)行的数据内容分别用步骤STC11所述ResAC的第1，2，3，4，…，RowsAC行的数据内容进行填充，步骤STC12所述FinalRes的第2，4，6，8，…，(RowsAC+RowsBD)行的数据内容分别用步骤STC10所述ResBD的第1，2，3，4，…，RowsBD行的数据内容进行填充，其中RowsAC表示步骤STC11所述ResAC的总行数，RowsBD表示步骤STC10所述ResBD的总行数；

STC14.完成双重列滤波操作，其最终结果输出为步骤STC12所述的FinalRes。

在给定两组序列长度相等的一维滤波器h_f、h_s和滤波对象X的前提下，所述双重行滤波包括以下步骤：

STR1.对X进行列对称延拓得到X_ex，如果X的总行数和总列数分别为RowsX和ColsX，h_f和h_s的序列长度为len，那么X_ex的总行数和总列数分别为RowsX和(ColsX+2×len)，并且X_ex的第1，2，…，len列的数据内容分别用X的第len，len-1，…，1列的数据内容进行填充，X_ex的第len+1，len+2，…，len+ColsX列的数据内容分别用X的第1，2，…，ColsX列的数据内容进行填充，X_ex的第len+ColsX+1，len+ColsX+2，…，len+ColsX+len列的数据内容分别用X的第ColsX，ColsX-1，…，ColsX-(len-1)列的数据内容进行填充；

STR2.抽取h_f中序号为奇数的滤波系数构成子滤波器h_fo，抽取h_f中序号为偶数的滤波系数构成子滤波器h_fe，抽取h_s中序号为奇数的滤波系数构成子滤波器h_so，抽取h_s中序号为偶数的滤波系数构成子滤波器h_se；

STR3.对步骤STR1所述的X_ex进行修正，即舍弃X_ex的起始两列和末尾两列，从而得到X′_ex；

STR4.抽取步骤STR3所述X′_ex中列号为奇数的数据构成X′_exO，抽取步骤STR3所述X′_ex中列号为偶数的数据构成X′_exE；

STR5.抽取步骤STR4所述X′_exE中列号为偶数的数据构成X′_exA，抽取步骤STR4所述X′_exO中列号为偶数的数据构成X′_exB，抽取步骤STR4所述X′_exE中列号为奇数的数据构成X′_exC，抽取步骤STR4所述X′_exO中列号为奇数的数据构成X′_exD；

STR6.用步骤STR2所述的子滤波器h_fo对步骤STR5所述的X′_exB进行逐行滤波，从而得到ResB；

STR7.用步骤STR2所述的子滤波器h_fe对步骤STR5所述的X′_exD进行逐行滤波，从而得到ResD；

STR8.用步骤STR2所述的子滤波器h_so对步骤STR5所述的X′_exA进行逐行滤波，从而得到ResA；

STR9.用步骤STR2所述的子滤波器h_se对步骤STR5所述的X′_exC进行逐行滤波，从而得到ResC；

STR10.把步骤STR6所述的ResB和步骤STR7所述的ResD中行号相同且列号也相同的元素相加得到ResBD，即ResBD＝ResB+ResD；

STR11.把步骤STR8所述的ResA和步骤STR9所述的ResC中行号相同且列号也相同的元素相加得到ResAC，即ResAC＝ResA+ResC；

STR12.申请一块总行数为RowsAC或RowsBD且总列数为(ColsAC+ColsBD)的数据存储区FinalRes，其中RowsAC表示步骤STR11所述ResAC的总行数，RowsBD表示步骤STR10所述ResBD的总行数，ColsAC表示步骤STR11所述ResAC的总列数，ColsBD表示步骤STR10所述ResBD的总列数，并且满足RowsAC＝RowsBD和ColsAC＝ColsBD的条件；

STR13.把一维滤波器h_f和h_s视为两个向量，如果这两个向量的内积大于0，那么步骤STR12所述FinalRes的第1，3，5，7，…，(ColsAC+ColsBD-1)列的数据内容分别用步骤STR10所述ResBD的第1，2，3，4，…，ColsBD列的数据内容进行填充，步骤STR12所述FinalRes的第2，4，6，8，…，(ColsAC+ColsBD)列的数据内容分别用步骤STR11所述ResAC的第1，2，3，4，…，ColsAC列的数据内容进行填充，如果这两个向量的内积小于或等于0，那么步骤STR12所述FinalRes的第1，3，5，7，…，(ColsAC+ColsBD-1)列的数据内容分别用步骤STR11所述ResAC的第1，2，3，4，…，ColsAC列的数据内容进行填充，步骤STR12所述FinalRes的第2，4，6，8，…，(ColsAC+ColsBD)列的数据内容分别用步骤STR10所述ResBD的第1，2，3，4，…，ColsBD列的数据内容进行填充，其中ColsAC表示步骤STR11所述ResAC的总列数，ColsBD表示步骤STR10所述ResBD的总列数；

STR14.完成双重行滤波操作，其最终结果输出为步骤STR12所述的FinalRes。

所述计算均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ包括以下步骤：

STEP1.把τ₀＝0.05作为τ参数的初始迭代值，进入以下步骤STEP2-STEP4的循环迭代过程；

STEP2.把τ作为自变量代入并计算

$f\left(\mathrm{\τ}\right)=\ln \left(\frac{1}{\mathrm{Num}}\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{2\mathrm{\τ}}\right)-2\ln \left(\frac{1}{\mathrm{Num}}\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)-\ln [1+\mathrm{\τ}(\frac{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\ln x_i}{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}}-\frac{1}{\mathrm{Num}}\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}\ln x_i)],$

其中x_i表示均方根小波系数子带中的第i个均方根小波系数，Num表示均方根小波系数子带中的均方根小波系数总数，ln表示以自然常数e为底数的对数函数，即自然对数；

STEP3.把τ作为自变量代入并计算

$\begin{array}{c}g\left(\mathrm{\τ}\right)=\frac{2\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{2\mathrm{\τ}}\ln x_i}{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{2\mathrm{\τ}}}-\frac{2\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\ln x_i}{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}}\\ -\frac{\mathrm{Num}\·\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)\·\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\ln x_i\right)-\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}\ln x_i\right)\·{\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)}^2+\mathrm{Num}\·\mathrm{\τ}[\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)\·\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}{\ln }^2{x}_i\right)-{\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\ln x_i\right)}^2]}{\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)\·\{\mathrm{Num}\·\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)+\mathrm{\τ}[\mathrm{Num}\·\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\ln x_i\right)-\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}\ln x_i\right)\·\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)\left]\right\}}\end{array}$

其中x_i表示均方根小波系数子带中的第i个均方根小波系数，Num表示均方根小波系数子带中的均方根小波系数总数，ln表示以自然常数e为底数的对数函数，即自然对数；

STEP4.计算若|Δτ|大于计算精度ε，则通过公式τ＝τ-|Δτ|来更新τ，然后执行步骤STEP2-STEP4的循环迭代过程，若|Δτ|小于或等于计算精度ε，则把当前τ值作为最终τ值并进入步骤STEP5；

STEP5.根据步骤STEP4获得的最终τ值，计算得到最终的α值和λ值，最终的α值，λ值的计算公式分别为

$\mathrm{\α}=\frac{1}{\mathrm{\τ}\·(\frac{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\ln x_i}{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}}-\frac{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}\ln x_i}{\mathrm{Num}})},$

$\mathrm{\λ}={[\frac{\mathrm{\τ}}{\mathrm{Num}}\·\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)\·(\frac{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\ln x_i}{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}}-\frac{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}\ln x_i}{\mathrm{Num}})]}^{\frac{1}{\mathrm{\τ}}},$

其中x_i表示均方根小波系数子带中的第i个均方根小波系数，Num表示均方根小波系数子带中的均方根小波系数总数，ln表示以自然常数e为底数的对数函数，即自然对数。

本发明所述6×N×3维的特征列向量的数据结构采用[α_1-1，τ_1-1，λ_1-1，α_1-2，τ_1-2，λ_1-2，α_1-3，τ_1-3，λ_1-3，α_1-4，τ_1-4，λ_1-4，α_1-5，τ_1-5，λ_1-5，α_1-6，τ_1-6，λ_1-6，α_2-1，τ_2-1，λ_2-1，α_2-2，τ_2-2，λ_2-2，α_2-3，τ_2-3，λ_2-3，α_2-4，τ_2-4，λ_2-4，α_2-5，τ_2-5，λ_2-5，α_2-6，τ_2-6，λ_2-6，…，α_N-1，τ_N-1，λ_N-1，α_N-2，τ_N-2，λ_N-2，α_N-3，τ_N-3，λ_N-3，α_N-4，τ_N-4，λ_N-4，α_N-5，τ_N-5，λ_N-5，α_N-6，τ_N-6，λ_N-6]^T的格式，特征列向量中的每一个α元素、τ元素、λ元素可以分别表示为α_i-n、τ_i-n、λ_i-n，其中下标i表示两路小波变换的级数序号，i的取值为1，2，…，N，n表示均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排列的排列序号，n的取值为1，2，…，6，α_i-n、τ_i-n和λ_i-n分别表示煤岩样本子图在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ。

所述比较特征列向量与特征列向量之间的相似度，从而判定待识别子图u_x所属煤岩类别的过程包括以下步骤：

S1.计算特征列向量与特征列向量之间的相似度，特征列向量与特征列向量之间的相似度计算公式，特征列向量与特征列向量之间的相似度计算公式分别为

$\begin{array}{c}\mathrm{Sm}(V_{u_x},V_{c_k})=\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}({\mathrm{ux}}_{3i-1}\·{\mathrm{ux}}_{3i-2}-{\mathrm{ck}}_{3i-1}\·{\mathrm{ck}}_{3i-2})[\ln \left(\frac{{\mathrm{ux}}_{3i}}{{\mathrm{ck}}_{3i}}\right)+\frac{\mathrm{\ψ}\left({\mathrm{ux}}_{3i-2}\right)}{{\mathrm{ux}}_{3i-1}}-\frac{\mathrm{\ψ}\left({\mathrm{ck}}_{3i-2}\right)}{{\mathrm{ck}}_{3i-1}}]\\ +\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}{\left(\frac{{\mathrm{ux}}_{3i}}{{\mathrm{ck}}_{3i}}\right)}^{{\mathrm{ck}}_{3i-1}}\·\frac{\mathrm{\Γ}({\mathrm{ux}}_{3i-2}+\frac{{\mathrm{ck}}_{3i-1}}{{\mathrm{ux}}_{3i-1}})}{\mathrm{\Γ}\left({\mathrm{ux}}_{3i-2}\right)}+\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}{\left(\frac{{\mathrm{ck}}_{3i}}{{\mathrm{ux}}_{3i}}\right)}^{{\mathrm{ux}}_{3i-1}}\·\frac{\mathrm{\Γ}({\mathrm{ck}}_{3i-2}+\frac{{\mathrm{ux}}_{3i-1}}{{\mathrm{ck}}_{3i-1}})}{\mathrm{\Γ}\left({\mathrm{ck}}_{3i-2}\right)}-\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}({\mathrm{ux}}_{3i-2}+{\mathrm{ck}}_{3i-2}),\end{array}$

$\begin{array}{c}\mathrm{Sm}(V_{u_x},V_{r_k})=\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}({\mathrm{ux}}_{3i-1}\·{\mathrm{ux}}_{3i-2}-{\mathrm{rk}}_{3i-1}\·{\mathrm{rk}}_{3i-2})[\ln \left(\frac{{\mathrm{ux}}_{3i}}{{\mathrm{rk}}_{3i}}\right)+\frac{\mathrm{\ψ}\left({\mathrm{ux}}_{3i-2}\right)}{{\mathrm{ux}}_{3i-1}}-\frac{\mathrm{\ψ}\left({\mathrm{rk}}_{3i-2}\right)}{{\mathrm{rk}}_{3i-1}}]\\ +\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}{\left(\frac{{\mathrm{ux}}_{3i}}{{\mathrm{rk}}_{3i}}\right)}^{{\mathrm{rk}}_{3i-1}}\·\frac{\mathrm{\Γ}({\mathrm{ux}}_{3i-2}+\frac{{\mathrm{rk}}_{3i-1}}{{\mathrm{ux}}_{3i-1}})}{\mathrm{\Γ}\left({\mathrm{ux}}_{3i-2}\right)}+\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}{\left(\frac{{\mathrm{rk}}_{3i}}{{\mathrm{ux}}_{3i}}\right)}^{{\mathrm{ux}}_{3i-1}}\·\frac{\mathrm{\Γ}({\mathrm{rk}}_{3i-2}+\frac{{\mathrm{ux}}_{3i-1}}{{\mathrm{rk}}_{3i-1}})}{\mathrm{\Γ}\left({\mathrm{rk}}_{3i-2}\right)}-\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}({\mathrm{ux}}_{3i-2}+{\mathrm{rk}}_{3i-2}),\end{array}$

其中表示特征列向量与特征列向量之间的相似度，表示特征列向量与特征列向量之间的相似度，特征列向量和中的下标k分别表示煤炭样本序号和岩石样本序号，k的取值为1，2，…，m，N表示两路小波变换的总级数，和中的每一个元素可以分别表示为ux_p、ck_p和rk_p，p的取值为1，2，…，(6×N×3)，ux_p表示的第p个元素，ck_p表示的第p个元素，rk_p表示的第p个元素，函数在Γ(z)函数表达式中，z为自变量，t为积分变量，积分区间为[0，+∞)，e表示自然常数，函数在ψ(z)函数表达式中，z为自变量，Γ′(z)为Γ(z)的导函数；

S2.进行相似度比较，如果那么判定待识别子图u_x属于煤炭，如果那么判定待识别子图u_x属于岩石。

附图说明

图1是基于变换域广义伽玛分布的煤岩识别方法的基本流程图；

图2是本发明所述N级两路小波变换的基本流程图；

图3是本发明所述双重列滤波的基本流程图；

图4是本发明所述双重行滤波的基本流程图；

图5是本发明所述计算均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ的基本流程图。

具体实施方式

在对我国山西、河南等地主要煤种和岩种的图像进行实验分析的基础上，本发明提出了基于变换域广义伽玛分布的煤岩识别方法，该方法可以有效判别煤岩与岩石。

下面结合附图和具体实施方式对本发明作进一步的详细描述。

参照图1，基于变换域广义伽玛分布的煤岩识别方法的具体步骤如下：

A.在样本训练阶段，获取相同光照条件下拍摄的m幅已知煤炭样本图像和m幅已知岩石样本图像，截取宽度和高度均为偶数个像素点并且不含非煤岩背景的子图像，煤炭样本子图和岩石样本子图分别记为c₁，c₂，…，c_m和r₁，r₂，…，r_m；

B.分别对样本子图c₁，c₂，…，c_m和r₁，r₂，…，r_m进行N级两路小波变换，每幅样本子图得到12×N个高频部分小波系数子带，每个高频部分小波系数子带记为表示样本子图在经过第i级两路小波变换以后的第d个方向的第j路高频部分小波系数子带，其中i表示两路小波变换的级数序号，i的取值为1，2，…，N，d表示两路小波变换的方向序号，d的取值为1，2，…，6，j表示两路小波变换的路序号，j的取值为1，2；

C.对步骤B所述的每幅样本子图的级数序号相同并且方向序号也相同的两个高频部分小波系数子带进行求均方根操作，得到6×N个均方根小波系数子带，每一个均方根小波系数子带记为RmsCoef_i-d，RmsCoef_i-d表示样本子图在经过第i级两路小波变换以后的第d个方向的均方根小波系数子带，RmsCoef_i-d中每一个元素的计算公式为

$\mathrm{Rms}{\mathrm{Coef}}_{i-d}\left(u\right)\left(v\right)=\sqrt{\frac{{\left[{\mathrm{Coef}}_{i-d}^1\left(u\right)\left(v\right)\right]}^2+{\left[{\mathrm{Coef}}_{i-d}^2\left(u\right)\left(v\right)\right]}^2}{2}},$

式中u和v分别表示RmsCoef_i-d、的行号和列号，RmsCoef_i-d(u)(v)表示RmsCoef_i-d中行号为u、列号为v的元素，表示中行号为u、列号为v的元素，表示中行号为u、列号为v的元素；

D.分别计算步骤C所述的每一个均方根小波系数子带RmsCoef_i-d的均值μ_i-d和方差μ_i-d的计算公式为其中U和V分别表示RmsCoef_i-d的总行数和总列数，的计算公式为 ${\mathrm{\σ}}_{i-d}^2=\frac{1}{\mathrm{UV}-1}\underset{u=1}{\overset{U}{\mathrm{\Σ}}}\left\{\underset{r=1}{\overset{V}{\mathrm{\Σ}}}{[\mathrm{Rms}{\mathrm{Coef}}_{i-d}\left(u\right)\left(v\right)-{\mathrm{\μ}}_{i-d}]}^2\right\},$ 其中U和V分别表示RmsCoef_i-d的总行数和总列数；

E.根据步骤D所计算的每一个均方根小波系数子带的均值和方差，对步骤C所述的每一级两路小波变换以后的6个方向的均方根小波系数子带按照均值和方差的乘积值从大到小的顺序排列，每一个经过排序以后的均方根小波系数子带记为RmsCoefS_i-n，RmsCoefS_i-n表示样本子图在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带，n表示均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排列的排列序号，n的取值为1，2，…，6；

F.分别计算步骤E所述的经过排序以后的每一个均方根小波系数子带在服从概率密度函数为的广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ，其中x表示均方根小波系数子带中的元素，也就是均方根小波系数，函数在Г(α)函数表达式中，α为自变量，t为积分变量，积分区间为[0，+∞)，e表示自然常数；

G.利用每一幅煤炭样本子图c_k或岩石样本子图r_k经过N级两路小波变换、高频部分小波系数子带求均方根操作并且按上述步骤E排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ，构造出一个6×N×3维的特征列向量或其中下标k表示样本序号，k的取值为1，2，…，m，然后把这2×m个特征列向量保存到数据存储设备中，用于后续的煤岩识别阶段；

H.在煤岩识别阶段，获取相同光照条件下拍摄的未知类别煤岩图像，截取宽度和高度均为偶数个像素点并且不含非煤岩背景的待识别子图u_x；

I.与上述步骤B类似，对u_x进行N级两路小波变换，从而得到12×N个高频部分小波系数子带，每个高频部分小波系数子带记为表示待识别子图u_x在经过第i级两路小波变换以后的第d个方向的第j路高频部分小波系数子带；

J.与上述步骤C类似，对步骤I所述的待识别子图u_x的级数序号相同并且方向序号也相同的两个高频部分小波系数子带进行求均方根操作，得到6×N个均方根小波系数子带，每一个均方根小波系数子带记为RmsUxCf_i-d，RmsUxCf_i-d表示待识别子图u_x在经过第i级两路小波变换以后的第d个方向的均方根小波系数子带，RmsUxCf_i-d中每一个元素的计算公式为

$\mathrm{Rms}{\mathrm{UxCf}}_{i-d}\left(p\right)\left(q\right)=\sqrt{\frac{{\left[{\mathrm{UxCf}}_{i-d}^1\left(p\right)\left(q\right)\right]}^2+{\left[{\mathrm{UxCf}}_{i-d}^2\left(p\right)\left(q\right)\right]}^2}{2}},$

式中p和q分别表示RmsUxCf_i-d、的行号和列号，RmsUxCf_i-d(p)(q)表示RmsUxCf_i-d中行号为p、列号为q的元素，表示中行号为p、列号为q的元素，表示中行号为p、列号为q的元素；

K.与上述步骤D类似，分别计算步骤J所述的每一个均方根小波系数子带RmsUxCf_i-d的均值和方差的计算公式为 ${\widehat{\mathrm{\μ}}}_{i-d}=\frac{1}{\mathrm{PQ}}\underset{p=1}{\overset{P}{\mathrm{\Σ}}}\left[\underset{q=1}{\overset{Q}{\mathrm{\Σ}}}\mathrm{RmsUxC}{f}_{i-d}\left(p\right)\left(q\right)\right],$ 其中P和Q分别表示均方根小波系数子带RmsUxCf_i-d的总行数和总列数，的计算公式为 ${\widehat{\mathrm{\σ}}}_{i-d}^2=\frac{1}{\mathrm{PQ}}\underset{p=1}{\overset{P}{\mathrm{\Σ}}}\left\{\right[\underset{q=1}{\overset{Q}{\mathrm{\Σ}}}\mathrm{RmsUxC}{f}_{i-d}\left(p\right)\left(q\right)-{\widehat{\mathrm{\μ}}}_{i-d}{]}^2\},$ 其中P和Q分别表示均方根小波系数子带RmsUxCf_i-d的总行数和总列数；

L.与上述步骤E类似，根据步骤K所计算的每一个均方根小波系数子带的均值和方差，对步骤J所述的每一级两路小波变换以后的6个方向的均方根小波系数子带按照均值和方差的乘积值从大到小的顺序排列，每一个经过排序以后的均方根小波系数子带记为RmsUxCfS_i-n，RmsUxCfS_i-n表示待识别子图u_x在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带；

M.与上述步骤F类似，分别计算步骤L所述的经过排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ；

N.与上述步骤G类似，利用待识别子图u_x经过N级两路小波变换、高频部分小波系数子带求均方根操作并且按上述步骤L排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ，构造出一个6×N×3维的特征列向量

O.比较特征列向量与特征列向量之间的相似度，其中下标k表示样本序号，k的取值为1，2，…，m，从而判定待识别子图u_x所属的煤岩类别。

参照图2，N级两路小波变换的具体步骤如下：

(1)初始化设置两路小波变换的级数序号i＝1；

(2)彩色图像灰度化处理，即把24位煤岩彩色图像转换成8位煤岩灰度图像，采用的转换公式为Y＝0.299R+0.587G+0.114B，其中R，G和B分别表示转换前彩色图像的红色，绿色和蓝色分量，Y表示转换后灰度图像的像素值；

(3)用一维的低通滤波器h₀对煤岩灰度图像进行逐列滤波，从而得到低频部分小波系数子带Coef_Lo；

(4)用一维的高通滤波器h₁对煤岩灰度图像进行逐列滤波，从而得到高频部分小波系数子带Coef_Hi；

(5)用步骤(3)所述的低通滤波器h₀对步骤(4)所述的高频部分小波系数子带Coef_Hi进行逐行滤波，从而得到高频部分小波系数子带Coef_HL；

(6)用步骤(4)所述的高通滤波器h₁对步骤(3)所述的低频部分小波系数子带Coef_Lo进行逐行滤波，从而得到高频部分小波系数子带Coef_LH；

(7)用步骤(4)所述的高通滤波器h₁对步骤(4)所述的高频部分小波系数子带Coef_Hi进行逐行滤波，从而得到高频部分小波系数子带Coef_HH；

(8)用步骤(3)所述的低通滤波器h₀对步骤(3)所述的低频部分小波系数子带Coef_Lo进行逐行滤波，从而得到低频部分小波系数子带Coef_LL；

(9)利用步骤(5)所述的高频部分小波系数子带Coef_HL，分别构造煤岩图像在经过第1级两路小波变换以后的第1个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第1个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第6个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第6个方向的第2路高频部分小波系数子带它们中各元素的计算公式分别为

${\mathrm{Coef}}_{1-1}^1\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HL}}(2u-1)(2v-1)-{\mathrm{Coef}}_{\mathrm{HL}}\left(2u\right)\left(2v\right),$

${\mathrm{Coef}}_{1-1}^2\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HL}}(2u-1)\left(2v\right)+{\mathrm{Coef}}_{\mathrm{HL}}\left(2u\right)(2v-1),$

${\mathrm{Coef}}_{1-6}^1\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HL}}(2u-1)(2v-1)+{\mathrm{Coef}}_{\mathrm{HL}}\left(2u\right)\left(2v\right),$

${\mathrm{Coef}}_{1-6}^2\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HL}}(2u-1)\left(2v\right)-{\mathrm{Coef}}_{\mathrm{HL}}\left(2u\right)(2v-1),$

上述各式中u和v分别表示的行号和列号， 和分别表示和中行号为u、列号为v的元素，Coef_HL(2u-1)(2v-1)表示Coef_HL中行号为2u-1、列号为2v-1的元素，Coef_HL(2u)(2v)表示Coef_HL中行号为2u、列号为2v的元素，Coef_HL(2u-1)(2v)表示Coef_HL中行号为2u-1、列号为2v的元素，Coef_HL(2u)(2v-1)表示Coef_HL中行号为2u、列号为2v-1的元素，u和v的最小取值为1，u的最大取值为Coef_HL总行数的一半，v的最大取值为Coef_HL总列数的一半；

(10)利用步骤(6)所述的高频部分小波系数子带Coef_LH，分别构造煤岩图像在经过第1级两路小波变换以后的第3个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第3个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第4个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第4个方向的第2路高频部分小波系数子带它们中各元素的计算公式分别为

${\mathrm{Coef}}_{1-3}^1\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{LH}}(2u-1)(2v-1)-{\mathrm{Coef}}_{\mathrm{LH}}\left(2u\right)\left(2v\right),$

${\mathrm{Coef}}_{1-3}^2\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{LH}}(2u-1)\left(2v\right)+{\mathrm{Coef}}_{\mathrm{LH}}\left(2u\right)(2v-1),$

${\mathrm{Coef}}_{1-4}^1\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{LH}}(2u-1)(2v-1)+{\mathrm{Coef}}_{\mathrm{LH}}\left(2u\right)\left(2v\right),$

${\mathrm{Coef}}_{1-4}^2\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{LH}}(2u-1)\left(2v\right)-{\mathrm{Coef}}_{\mathrm{LH}}\left(2u\right)(2v-1),$

上述各式中u和v分别表示的行号和列号， 和分别表示和中行号为u、列号为v的元素，Coef_LH(2u-1)(2v-1)表示Coef_LH中行号为2u-1、列号为2v-1的元素，Coef_LH(2u)(2v)表示Coef_LH中行号为2u、列号为2v的元素，Coef_LH(2u-1)(2v)表示Coef_LH中行号为2u-1、列号为2v的元素，Coef_LH(2u)(2v-1)表示Coef_LH中行号为2u、列号为2v-1的元素，u和v的最小取值为1，u的最大取值为Coef_LH总行数的一半，v的最大取值为Coef_LH总列数的一半；

(11)利用步骤(7)所述的高频部分小波系数子带Coef_HH，分别构造煤岩图像在经过第1级两路小波变换以后的第2个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第2个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第5个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第5个方向的第2路高频部分小波系数子带它们中各元素的计算公式分别为

${\mathrm{Coef}}_{1-2}^1\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HH}}(2u-1)(2v-1)-{\mathrm{Coef}}_{\mathrm{HH}}\left(2u\right)\left(2v\right),$

${\mathrm{Coef}}_{1-2}^2\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HH}}(2u-1)\left(2v\right)+{\mathrm{Coef}}_{\mathrm{HH}}\left(2u\right)(2v-1),$

${\mathrm{Coef}}_{1-5}^1\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HH}}(2u-1)(2v-1)+{\mathrm{Coef}}_{\mathrm{HH}}\left(2u\right)\left(2v\right),$

${\mathrm{Coef}}_{1-5}^2\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HH}}(2u-1)\left(2v\right)-{\mathrm{Coef}}_{\mathrm{HH}}\left(2u\right)(2v-1),$

上述各式中u和v分别表示的行号和列号， 和分别表示和中行号为u、列号为v的元素，Coef_HH(2u-1)(2v-1)表示Coef_HH中行号为2u-1、列号为2v-1的元素，Coef_HH(2u)(2v)表示Coef_HH中行号为2u、列号为2v的元素，Coef_HH(2u-1)(2v)表示Coef_HH中行号为2u-1、列号为2v的元素，Coef_HH(2u)(2v-1)表示Coef_HH中行号为2u、列号为2v-1的元素，u和v的最小取值为1，u的最大取值为Coef_HH总行数的一半，v的最大取值为Coef_HH总列数的一半；

(12)两路小波变换的级数序号自增1，即i＝i+1；

(13)动态申请与Coef_LL数据结构完全一致的临时数据存储区TempCoef，并将其数据内容初始化为Coef_LL的数据内容；

(14)判断i≤N条件是否成立，如果是，则进入以下步骤(15)-(28)的迭代循环，如果否，则转到步骤(29)；

(15)如果TempCoef数据的总行数不是4的倍数，那么对TempCoef数据进行修正，在TempCoef数据的第一行Row_first前面追加一行Row′_first并且Row′_first的数据内容用Row_first的数据内容进行填充，在TempCoef数据的最后一行Row_last后面追加一行Row′_last并且Row′_last的数据内容用Row_last的数据内容进行填充；

(16)如果TempCoef数据的总列数不是4的倍数，那么对TempCoef数据进行修正，在TempCoef数据的第一列Col_first前面追加一列Col′_first并且Col′_first的数据内容用Col_first的数据内容进行填充，在TempCoef数据的最后一列Col_last后面追加一列Col′_last并且Col′_last的数据内容用Col_last的数据内容进行填充；

(17)用两组一维的低通滤波器h₀₀、h₀₁对TempCoef数据进行双重列滤波，从而得到低频部分小波系数子带Coef_L2-i，其中下标i表示两路小波变换的级数序号；

(18)用两组一维的高通滤波器h₁₀、h₁₁对TempCoef数据进行双重列滤波，从而得到高频部分小波系数子带Coef_H2-i，其中下标i表示两路小波变换的级数序号；

(19)用步骤(17)所述的低通滤波器h₀₀、h₀₁对步骤(18)所述的高频部分小波系数子带Coef_H2-i进行双重行滤波，从而得到高频部分小波系数子带Coef_HL2-i，其中下标i表示两路小波变换的级数序号；

(20)用步骤(18)所述的高通滤波器h₁₀、h₁₁对步骤(17)所述的低频部分小波系数子带Coef_L2-i进行双重行滤波，从而得到高频部分小波系数子带Coef_LH2-i，其中下标i表示两路小波变换的级数序号；

(21)用步骤(18)所述的高通滤波器h₁₀、h₁₁对步骤(18)所述的高频部分小波系数子带Coef_H2-i进行双重行滤波，从而得到高频部分小波系数子带Coef_HH2-i，其中下标i表示两路小波变换的级数序号；

(22)用步骤(17)所述的低通滤波器h₀₀、h₀₁对步骤(17)所述的低频部分小波系数子带Coef_L2-i进行双重行滤波，从而得到低频部分小波系数子带Coef_LL2-i，其中下标i表示两路小波变换的级数序号；

(23)利用步骤(19)所述的高频部分小波系数子带Coef_HL2-i，分别构造煤岩图像在经过第i级两路小波变换以后的第1个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第1个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第6个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第6个方向的第2路高频部分小波系数子带它们中各元素的计算公式分别为

${\mathrm{Coef}}_{i-1}^1\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HL}2-i}(2u-1)(2v-1)-{\mathrm{Coef}}_{\mathrm{HL}2-i}\left(2u\right)\left(2v\right),$

${\mathrm{Coef}}_{i-1}^2\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HL}2-i}(2u-1)\left(2v\right)+{\mathrm{Coef}}_{\mathrm{HL}2-i}\left(2u\right)(2v-1),$

${\mathrm{Coef}}_{i-6}^1\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HL}2-i}(2u-1)(2v-1)+{\mathrm{Coef}}_{\mathrm{HL}2-i}\left(2u\right)\left(2v\right),$

${\mathrm{Coef}}_{i-6}^2\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HL}2-i}(2u-1)\left(2v\right)+{\mathrm{Coef}}_{\mathrm{HL}2-i}\left(2u\right)(2v-1),$

上述各式中u和v分别表示的行号和列号， 和分别表示和中行号为u、列号为v的元素，Coef_HL2-i(2u-1)(2v-1)表示Coef_HL2-i中行号为2u-1、列号为2v-1的元素，Coef_HL2-i(2u)(2v)表示Coef_HL2-i中行号为2u、列号为2v的元素，Coef_HL2-i(2u-1)(2v)表示Coef_HL2-i中行号为2u-1、列号为2v的元素，Coef_HL2-i(2u)(2v-1)表示Coef_HL2-i中行号为2u、列号为2v-1的元素，u和v的最小取值为1，u的最大取值为Coef_HL2-i总行数的一半，v的最大取值为Coef_HL2-i总列数的一半；

(24)利用步骤(20)所述的高频部分小波系数子带Coef_LH2-i，分别构造煤岩图像在经过第i级两路小波变换以后的第3个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第3个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第4个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第4个方向的第2路高频部分小波系数子带它们中各元素的计算公式分别为

${\mathrm{Coef}}_{i-3}^1\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{LH}2-i}(2u-1)(2v-1)-{\mathrm{Coef}}_{\mathrm{LH}2-i}\left(2u\right)\left(2v\right),$

${\mathrm{Coef}}_{i-3}^2\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{LH}2-i}(2u-1)\left(2v\right)+{\mathrm{Coef}}_{\mathrm{LH}2-i}\left(2u\right)(2v-1),$

${\mathrm{Coef}}_{i-4}^1\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{LH}2-i}(2u-1)(2v-1)+{\mathrm{Coef}}_{\mathrm{LH}2-i}\left(2u\right)\left(2v\right),$

${\mathrm{Coef}}_{i-4}^2\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{LH}2-i}(2u-1)\left(2v\right)+{\mathrm{Coef}}_{\mathrm{LH}2-i}\left(2u\right)(2v-1),$

上述各式中u和v分别表示的行号和列号， 和分别表示和中行号为u、列号为v的元素，Coef_LH2-i(2u-1)(2v-1)表示Coef_LH2-i中行号为2u-1、列号为2v-1的元素，Coef_LH2-i(2u)(2v)表示Coef_LH2-i中行号为2u、列号为2v的元素，Coef_LH2-i(2u-1)(2v)表示Coef_LH2-i中行号为2u-1、列号为2v的元素，Coef_LH2-i(2u)(2v-1)表示Coef_LH2-i中行号为2u、列号为2v-1的元素，u和v的最小取值为1，u的最大取值为Coef_LH2-i总行数的一半，v的最大取值为Coef_LH2-i总列数的一半；

(25)利用步骤(21)所述的高频部分小波系数子带Coef_HH2-i，分别构造煤岩图像在经过第i级两路小波变换以后的第2个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第2个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第5个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第5个方向的第2路高频部分小波系数子带它们中各元素的计算公式分别为

${\mathrm{Coef}}_{i-2}^1\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HH}2-i}(2u-1)(2v-1)-{\mathrm{Coef}}_{\mathrm{HH}2-i}\left(2u\right)\left(2v\right),$

${\mathrm{Coef}}_{i-2}^2\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HH}2-i}(2u-1)\left(2v\right)+{\mathrm{Coef}}_{\mathrm{HH}2-i}\left(2u\right)(2v-1),$

${\mathrm{Coef}}_{i-5}^1\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HH}2-i}(2u-1)(2v-1)+{\mathrm{Coef}}_{\mathrm{HH}2-i}\left(2u\right)\left(2v\right),$

${\mathrm{Coef}}_{i-5}^2\left(u\right)\left(v\right)={\mathrm{Coef}}_{\mathrm{HH}2-i}(2u-1)\left(2v\right)+{\mathrm{Coef}}_{\mathrm{HH}2-i}\left(2u\right)(2v-1),$

上述各式中u和v分别表示的行号和列号， 和分别表示和中行号为u、列号为v的元素，Coef_HH2-i(2u-1)(2v-1)表示Coef_HH2-i中行号为2u-1、列号为2v-1的元素，Coef_HH2-i(2u)(2v)表示Coef_HH2-i中行号为2u、列号为2v的元素，Coef_HH2-i(2u-1)(2v)表示Coef_HH2-i中行号为2u-1、列号为2v的元素，Coef_HH2-i(2u)(2v-1)表示Coef_HH2-i中行号为2u、列号为2v-1的元素，u和v的最小取值为1，u的最大取值为Coef_HH2-i总行数的一半，v的最大取值为Coef_HH2-i总列数的一半；

(26)释放临时数据存储区TempCoef所占用的内存空间，重新动态申请与步骤(22)所述Coef_LL2-i的数据结构完全一致的临时数据存储区TempCoef，并将其数据内容初始化为Coef_LL2-i的数据内容；

(27)两路小波变换的级数序号自增1，即i＝i+1；

(28)转到步骤(14)，进行下一级两路小波变换；

(29)释放临时数据存储区TempCoef所占用的内存空间，完成N级两路小波变换。

所述一维的低通滤波器h₀和高通滤波器h₁的序列长度均为奇数，并且h₀序列和h₁序列均为中心对称序列。

所述一维的低通滤波器h₀₀、h₀₁和高通滤波器h₁₀、h₁₁的序列长度相等并且均为偶数，h₀₁序列为h₀₀序列的反转序列，h₁₀序列通过对h₀₀序列进行反转操作然后对偶数序号采样点的序列元素进行符号取反操作而得到，h₁₁序列为h₁₀序列的反转序列。

如图3所示，在给定两组序列长度相等的一维滤波器h_f、h_s和滤波对象X的前提下，所述双重列滤波包括以下步骤：

STC1.对X进行行对称延拓得到X_ex，如果X的总行数和总列数分别为RowsX和ColsX，h_f和h_s的序列长度为len，那么X_ex的总行数和总列数分别为(RowsX+2×len)和ColsX，并且X_ex的第1，2，…，len行的数据内容分别用X的第len，len-1，…，1行的数据内容进行填充，X_ex的第len+1，len+2，…，len+RowsX行的数据内容分别用X的第1，2，…，RowsX行的数据内容进行填充，X_ex的第len+RowsX+1，len+RowsX+2，…，len+RowsX+len行的数据内容分别用X的第RowsX，RowsX-1，…，RowsX-(len-1)行的数据内容进行填充；

STC2.抽取h_f中序号为奇数的滤波系数构成子滤波器h_fo，抽取h_f中序号为偶数的滤波系数构成子滤波器h_fe，抽取h_s中序号为奇数的滤波系数构成子滤波器h_so，抽取h_s中序号为偶数的滤波系数构成子滤波器h_se；

STC3.对步骤STC1所述的X_ex进行修正，即舍弃X_ex的起始两行和末尾两行，从而得到X′_ex；

STC4.抽取步骤STC3所述X′_ex中行号为奇数的数据构成X′_exO，抽取步骤STC3所述X′_ex中行号为偶数的数据构成X′_exE；

STC5.抽取步骤STC4所述X′_exE中行号为偶数的数据构成X′_exA，抽取步骤STC4所述X′_exO中行号为偶数的数据构成X′_exB，抽取步骤STC4所述X′_exE中行号为奇数的数据构成X′_exC，抽取步骤STC4所述X′_exO中行号为奇数的数据构成X′_exD；

STC6.用步骤STC2所述的子滤波器h_fo对步骤STC5所述的X′_exB进行逐列滤波，从而得到ResB；

STC7.用步骤STC2所述的子滤波器h_fe对步骤STC5所述的X′_exD进行逐列滤波，从而得到ResD；

STC8.用步骤STC2所述的子滤波器h_so对步骤STC5所述的X′_exA进行逐列滤波，从而得到ResA；

STC9.用步骤STC2所述的子滤波器h_se对步骤STC5所述的X′_exC进行逐列滤波，从而得到ResC；

STC10.把步骤STC6所述的ResB和步骤STC7所述的ResD中行号相同且列号也相同的元素相加得到ResBD，即ResBD＝ResB+ResD；

STC11.把步骤STC8所述的ResA和步骤STC9所述的ResC中行号相同且列号也相同的元素相加得到ResAC，即ResAC＝ResA+ResC；

STC12.申请一块总行数为(RowsAC+RowsBD)且总列数为ColsAC或ColsBD的数据存储区FinalRes，其中RowsAC表示步骤STC11所述ResAC的总行数，RowsBD表示步骤STC10所述ResBD的总行数，ColsAC表示步骤STC11所述ResAC的总列数，ColsBD表示步骤STC10所述ResBD的总列数，并且满足RowsAC＝RowsBD和ColsAC＝ColsBD的条件；

STC13.把一维滤波器h_f和h_s视为两个向量，如果这两个向量的内积大于0，那么步骤STC12所述FinalRes的第1，3，5，7，…，(RowsAC+RowsBD-1)行的数据内容分别用步骤STC10所述ResBD的第1，2，3，4，…，RowsBD行的数据内容进行填充，步骤STC12所述FinalRes的第2，4，6，8，…，(RowsAC+RowsBD)行的数据内容分别用步骤STC11所述ResAC的第1，2，3，4，…，RowsAC行的数据内容进行填充，如果这两个向量的内积小于或等于0，那么步骤STC12所述FinalRes的第1，3，5，7，…，(RowsAC+RowsBD-1)行的数据内容分别用步骤STC11所述ResAC的第1，2，3，4，…，RowsAC行的数据内容进行填充，步骤STC12所述FinalRes的第2，4，6，8，…，(RowsAC+RowsBD)行的数据内容分别用步骤STC10所述ResBD的第1，2，3，4，…，RowsBD行的数据内容进行填充，其中RowsAC表示步骤STC11所述ResAC的总行数，RowsBD表示步骤STC10所述ResBD的总行数；

STC14.完成双重列滤波操作，其最终结果输出为步骤STC12所述的FinalRes。

如图4所示，在给定两组序列长度相等的一维滤波器h_f、h_s和滤波对象X的前提下，所述双重行滤波包括以下步骤：

STR1.对X进行列对称延拓得到X_ex，如果X的总行数和总列数分别为RowsX和ColsX，h_f和h_s的序列长度为len，那么X_ex的总行数和总列数分别为RowsX和(ColsX+2×len)，并且X_ex的第1，2，…，len列的数据内容分别用X的第len，len-1，…，1列的数据内容进行填充，X_ex的第len+1，len+2，…，len+ColsX列的数据内容分别用X的第1，2，…，ColsX列的数据内容进行填充，X_ex的第len+ColsX+1，len+ColsX+2，…，len+ColsX+len列的数据内容分别用X的第ColsX，ColsX-1，…，ColsX-(len-1)列的数据内容进行填充；

STR2.抽取h_f中序号为奇数的滤波系数构成子滤波器h_fo，抽取h_f中序号为偶数的滤波系数构成子滤波器h_fe，抽取h_s中序号为奇数的滤波系数构成子滤波器h_so，抽取h_s中序号为偶数的滤波系数构成子滤波器h_se；

STR3.对步骤STR1所述的X_ex进行修正，即舍弃X_ex的起始两列和末尾两列，从而得到X′_ex；

STR4.抽取步骤STR3所述X′_ex中列号为奇数的数据构成X′_exO，抽取步骤STR3所述X′_ex中列号为偶数的数据构成X′_exE；

STR5.抽取步骤STR4所述X′_exE中列号为偶数的数据构成X′_exA，抽取步骤STR4所述X′_exO中列号为偶数的数据构成X′_exB，抽取步骤STR4所述X′_exE中列号为奇数的数据构成X′_exC，抽取步骤STR4所述X′_exO中列号为奇数的数据构成X′_exD；

STR6.用步骤STR2所述的子滤波器h_fo对步骤STR5所述的X′_exB进行逐行滤波，从而得到ResB；

STR7.用步骤STR2所述的子滤波器h_fe对步骤STR5所述的X′_exD进行逐行滤波，从而得到ResD；

STR8.用步骤STR2所述的子滤波器h_so对步骤STR5所述的X′_exA进行逐行滤波，从而得到ResA；

STR9.用步骤STR2所述的子滤波器h_se对步骤STR5所述的X′_exC进行逐行滤波，从而得到ResC；

STR10.把步骤STR6所述的ResB和步骤STR7所述的ResD中行号相同且列号也相同的元素相加得到ResBD，即ResBD＝ResB+ResD；

STR11.把步骤STR8所述的ResA和步骤STR9所述的ResC中行号相同且列号也相同的元素相加得到ResAC，即ResAC＝ResA+ResC；

STR12.申请一块总行数为RowsAC或RowsBD且总列数为(ColsAC+ColsBD)的数据存储区FinalRes，其中RowsAC表示步骤STR11所述ResAC的总行数，RowsBD表示步骤STR10所述ResBD的总行数，ColsAC表示步骤STR11所述ResAC的总列数，ColsBD表示步骤STR10所述ResBD的总列数，并且满足RowsAC＝RowsBD和ColsAC＝ColsBD的条件；

STR13.把一维滤波器h_f和h_s视为两个向量，如果这两个向量的内积大于0，那么步骤STR12所述FinalRes的第1，3，5，7，…，(ColsAC+ColsBD-1)列的数据内容分别用步骤STR10所述ResBD的第1，2，3，4，…，ColsBD列的数据内容进行填充，步骤STR12所述FinalRes的第2，4，6，8，…，(ColsAC+ColsBD)列的数据内容分别用步骤STR11所述ResAC的第1，2，3，4，…，ColsAC列的数据内容进行填充，如果这两个向量的内积小于或等于0，那么步骤STR12所述FinalRes的第1，3，5，7，…，(ColsAC+ColsBD-1)列的数据内容分别用步骤STR11所述ResAC的第1，2，3，4，…，ColsAC列的数据内容进行填充，步骤STR12所述FinalRes的第2，4，6，8，…，(ColsAC+ColsBD)列的数据内容分别用步骤STR10所述ResBD的第1，2，3，4，…，ColsBD列的数据内容进行填充，其中ColsAC表示步骤STR11所述ResAC的总列数，ColsBD表示步骤STR10所述ResBD的总列数；

STR14.完成双重行滤波操作，其最终结果输出为步骤STR12所述的FinalRes。

如图5所示，所述计算均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ包括以下步骤：

STEP1.把τ₀＝0.05作为τ参数的初始迭代值，进入以下步骤STEP2-STEP4的循环迭代过程；

STEP2.把τ作为自变量代入并计算

$f\left(\mathrm{\τ}\right)=\ln \left(\frac{1}{\mathrm{Num}}\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{2\mathrm{\τ}}\right)-2\ln \left(\frac{1}{\mathrm{Num}}\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)-\ln [1+\mathrm{\τ}(\frac{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\ln x_i}{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}}-\frac{1}{\mathrm{Num}}\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}\ln x_i)],$

其中x_i表示均方根小波系数子带中的第i个均方根小波系数，Num表示均方根小波系数子带中的均方根小波系数总数，ln表示以自然常数e为底数的对数函数，即自然对数；

STEP3.把τ作为自变量代入并计算

$g\left(\mathrm{\τ}\right)=\frac{2\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{2\mathrm{\τ}}\ln x_i}{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{2\mathrm{\τ}}}-\frac{2\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\ln x_i}{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}}$

$-\frac{\mathrm{Num}\·\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)\·\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\ln x_i\right)-\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}\ln x_i\right){\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)}^2+\mathrm{Num}\·\mathrm{\τ}[\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)\·\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}{\ln }^2{x}_i\right)-{\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\ln x_i\right)}^2]}{\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)\·\{\mathrm{Num}\·\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)+\mathrm{\τ}[\mathrm{Num}\·\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\ln x_i\right)-\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}\ln x_i\right)\·\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)\left]\right\}}$

其中x_i表示均方根小波系数子带中的第i个均方根小波系数，Num表示均方根小波系数子带中的均方根小波系数总数，ln表示以自然常数e为底数的对数函数，即自然对数；

STEP4.计算若|Δτ|大于计算精度ε，则通过公式τ＝τ-|Δτ|来更新τ，然后执行步骤STEP2-STEP4的循环迭代过程，若|Δτ|小于或等于计算精度ε，则把当前τ值作为最终τ值并进入步骤STEP5；

STEP5.根据步骤STEP4获得的最终τ值，计算得到最终的α值和λ值，最终的α值，λ值的计算公式分别为

$\mathrm{\α}=\frac{1}{\mathrm{\τ}\·(\frac{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\ln x_i}{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}}-\frac{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}\ln x_i}{\mathrm{Num}})},$

$\mathrm{\λ}={[\frac{\mathrm{\τ}}{\mathrm{Num}}\·\left(\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\right)\·(\frac{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}\ln x_i}{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}{x}_i^{\mathrm{\τ}}}-\frac{\underset{i=1}{\overset{\mathrm{Num}}{\mathrm{\Σ}}}\ln x_i}{\mathrm{Num}})]}^{\frac{1}{\mathrm{\τ}}},$

其中x_i表示均方根小波系数子带中的第i个均方根小波系数，Num表示均方根小波系数子带中的均方根小波系数总数，ln表示以自然常数e为底数的对数函数，即自然对数。

本发明所述6×N×3维的特征列向量的数据结构采用[α_1-1，τ_1-1，λ_1-1，α_1-2，τ_1-2，λ_1-2，α_1-3，τ_1-3，λ_1-3，α_1-4，τ_1-4，λ_1-4，α_1-5，τ_1-5，λ_1-5，α_1-6，τ_1-6，λ_1-6，α_2-1，τ_2-1，λ_2-1，α_2-2，τ_2-2，λ_2-2，α_2-3，τ_2-3，λ_2-3，α_2-4，τ_2-4，λ_2-4，α_2-5，τ_2-5，λ_2-5，α_2-6，τ_2-6，λ_2-6，…，α_N-1，τ_N-1，λ_N-1，α_N-2，τ_N-2，λ_N-2，α_N-3，τ_N-3，λ_N-3，α_N-4，τ_N-4，λ_N-4，α_N-5，τ_N-5，λ_N-5，α_N-6，τ_N-6，λ_N-6]^T的格式，特征列向量中的每一个α元素、τ元素、λ元素可以分别表示为α_i-n、τ_i-n、λ_i-n，其中下标i表示两路小波变换的级数序号，i的取值为1，2，…，N，n表示均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排列的排列序号，n的取值为1，2，…，6，α_i-n、τ_i-n和λ_i-n分别表示煤岩样本子图在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ。

所述比较特征列向量与特征列向量之间的相似度，从而判定待识别子图u_x所属煤岩类别的过程包括以下步骤：

S1.计算特征列向量与特征列向量之间的相似度，特征列向量与特征列向量之间的相似度计算公式，特征列向量与特征列向量之间的相似度计算公式分别为

$\begin{array}{c}\mathrm{Sm}(V_{u_x},V_{c_k})=\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}({\mathrm{ux}}_{3i-1}\·{\mathrm{ux}}_{3i-2}-{\mathrm{ck}}_{3i-1}\·{\mathrm{ck}}_{3i-2})\left[\ln \left(\frac{{\mathrm{ux}}_{3i}}{{\mathrm{ck}}_{3i}}\right)+\frac{\mathrm{\ψ}\left({\mathrm{ux}}_{3i-2}\right)}{{\mathrm{ux}}_{3i-1}}-\frac{\mathrm{\ψ}\left({\mathrm{ck}}_{3i-2}\right)}{{\mathrm{ck}}_{3i-1}}\right]\\ +\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}{\left(\frac{{\mathrm{uk}}_{3i}}{{\mathrm{ck}}_{3i}}\right)}^{{\mathrm{ck}}_{3i-1}}\·\frac{\mathrm{\Γ}({\mathrm{ux}}_{3i-2}+\frac{{\mathrm{ck}}_{3i-1}}{{\mathrm{ux}}_{3i-1}})}{\mathrm{\Γ}\left({\mathrm{uk}}_{3i-2}\right)}+\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}{\left(\frac{{\mathrm{ck}}_{3i}}{{\mathrm{ux}}_{3i}}\right)}^{{\mathrm{ux}}_{3i-1}}\·\frac{\mathrm{\Γ}({\mathrm{ck}}_{3i-2}+\frac{{\mathrm{uk}}_{3i-1}}{{\mathrm{ck}}_{3i-1}})}{\mathrm{\Γ}\left({\mathrm{ck}}_{3i-2}\right)}-\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}({\mathrm{uk}}_{3i-2}+{\mathrm{ck}}_{3i-2}),\end{array}$

$\begin{array}{c}\mathrm{Sm}(V_{u_x},V_{r_k})=\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}({\mathrm{ux}}_{3i-1}\·{\mathrm{ux}}_{3i-2}-{\mathrm{rk}}_{3i-1}\·{\mathrm{rk}}_{3i-2})\left[\ln \left(\frac{{\mathrm{ux}}_{3i}}{{\mathrm{rk}}_{3i}}\right)+\frac{\mathrm{\ψ}\left({\mathrm{ux}}_{3i-2}\right)}{{\mathrm{ux}}_{3i-1}}-\frac{\mathrm{\ψ}\left({\mathrm{rk}}_{3i-2}\right)}{{\mathrm{rk}}_{3i-1}}\right]\\ +\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}{\left(\frac{{\mathrm{uk}}_{3i}}{{\mathrm{rk}}_{3i}}\right)}^{{\mathrm{rk}}_{3i-1}}\·\frac{\mathrm{\Γ}({\mathrm{ux}}_{3i-2}+\frac{{\mathrm{rk}}_{3i-1}}{{\mathrm{ux}}_{3i-1}})}{\mathrm{\Γ}\left({\mathrm{ux}}_{3i-2}\right)}+\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}{\left(\frac{{\mathrm{rk}}_{3i}}{{\mathrm{ux}}_{3i}}\right)}^{{\mathrm{ux}}_{3i-1}}\·\frac{\mathrm{\Γ}({\mathrm{rk}}_{3i-2}+\frac{{\mathrm{uk}}_{3i-1}}{{\mathrm{rk}}_{3i-1}})}{\mathrm{\Γ}\left({\mathrm{rk}}_{3i-2}\right)}-\underset{i=1}{\overset{6\×N}{\mathrm{\Σ}}}({\mathrm{ux}}_{3i-2}+{\mathrm{rk}}_{3i-2}),\end{array}$

其中表示特征列向量与特征列向量之间的相似度，表示特征列向量与特征列向量之间的相似度，特征列向量和中的下标k分别表示煤炭样本序号和岩石样本序号，k的取值为1，2，…，m，N表示两路小波变换的总级数，和中的每一个元素可以分别表示为ux_p、ck_p和rk_p，p的取值为1，2，…，(6×N×3)，ux_p表示的第p个元素，ck_p表示的第p个元素，rk_p表示的第p个元素，函数在Γ(z)函数表达式中，z为自变量，t为积分变量，积分区间为[0，+∞)，e表示自然常数，函数在ψ(z)函数表达式中，z为自变量，Γ′(z)为Γ(z)的导函数；

S2.进行相似度比较，如果那么判定待识别子图u_x属于煤炭，如果那么判定待识别子图u_x属于岩石。

需要指出的是，以上所述实施实例用于进一步说明本发明，实施实例不应被视为限制本发明的范围。

Claims (9)

1.基于变换域广义伽玛分布的煤岩识别方法，其特征在于，包括以下步骤：

A.在样本训练阶段，获取相同光照条件下拍摄的m幅已知煤炭样本图像和m幅已知岩石样本图像，截取宽度和高度均为偶数个像素点并且不含非煤岩背景的子图像，煤炭样本子图和岩石样本子图分别记为c₁，c₂，…，c_m和r₁，r₂，…，r_m；

B.分别对样本子图c₁，c₂，…，c_m和r₁，r₂，…，r_m进行N级两路小波变换，每幅样本子图得到12×N个高频部分小波系数子带

表示样本子图在经过第i级两路小波变换以后的第d个方向的第j路高频部分小波系数子带，其中i表示两路小波变换的级数序号，i的取值为1，2，…，N，d表示两路小波变换的方向序号，d的取值为1，2，…，6，j表示两路小波变换的路序号，j的取值为1，2；

C.对步骤B所述的每幅样本子图的级数序号相同并且方向序号也相同的两个高频部分小波系数子带进行求均方根操作，得到6×N个均方根小波系数子带RmsCoef_i-d，RmsCoef_i-d表示样本子图在经过第i级两路小波变换以后的第d个方向的均方根小波系数子带；

D.分别计算步骤C所述的每一个均方根小波系数子带RmsCoef_i-d的均值μ_i-d和方差 

E.根据步骤D所计算的每一个均方根小波系数子带的均值和方差，对步骤C所述的每一级两路小波变换以后的6个方向的均方根小波系数子带按照均值和方差的乘积值从大到小的顺序排列，每一个经过排序以后的均方根小波系数子带记为RmsCoefS_i-n；

RmsCoefS_i-n表示样本子图在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带，n表示均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排列的排列序号，n的取值为1，2，…，6；

F.分别计算步骤E所述的经过排序以后的每一个均方根小波系数子带在服从概率密度函数为的广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ，其中x表示均方根小波系数子带中的元素，也就是均方根小波系数，函数 在Г(α)函数表达式中，α为自变量，t为积分变量，积分区间为[0，+∞)，e表示自然常数；

G.利用每一幅煤炭样本子图c_k或岩石样本子图r_k经过N级两路小波变换、高频部分小波 系数子带求均方根操作并且按上述步骤E排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ，构造出一个6×N×3维的特征列向量或其中下标k表示样本序号，k的取值为1，2，…，m，然后把这2×m个特征列向量保存到数据存储设备中，用于后续的煤岩识别阶段；

H.在煤岩识别阶段，获取相同光照条件下拍摄的未知类别煤岩图像，截取宽度和高度均为偶数个像素点并且不含非煤岩背景的待识别子图u_x；

I.与上述步骤B类似，对u_x进行N级两路小波变换，从而得到12×N个高频部分小波系数子带，每个高频部分小波系数子带记为表示待识别子图u_x在经过第i级两路小波变换以后的第d个方向的第j路高频部分小波系数子带；

J.与上述步骤C类似，对步骤I所述的待识别子图u_x的级数序号相同并且方向序号也相同的两个高频部分小波系数子带进行求均方根操作，得到6×N个均方根小波系数子带，每一个均方根小波系数子带记为RmsUxCf_i-d，RmsUxCf_i-d表示待识别子图u_x在经过第i级两路小波变换以后的第d个方向的均方根小波系数子带；

K.与上述步骤D类似，分别计算步骤J所述的每一个均方根小波系数子带RmsUxCf_i-d的均值和方差

L.与上述步骤E类似，根据步骤K所计算的每一个均方根小波系数子带的均值和方差，对步骤J所述的每一级两路小波变换以后的6个方向的均方根小波系数子带按照均值和方差的乘积值从大到小的顺序排列，每一个经过排序以后的均方根小波系数子带记为RmsUxCfS_i-n，RmsUxCfS_i-n表示待识别子图u_x在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带；

M.与上述步骤F类似，分别计算步骤L所述的经过排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ；

N.与上述步骤G类似，利用待识别子图u_x经过N级两路小波变换、高频部分小波系数子带求均方根操作并且按上述步骤L排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ，构造出一个6×N×3维的特征列向量

O.比较特征列向量与特征列向量之间的相似度，如果 那么判定待识别子图u_x属于煤炭，如果 那么判定待识别子图u_x属于岩石；

其中，表示特征列向量与特征列向量之间的相似度，表示特征列向量与特征列向量之间的相似度，其中下标k表示样本序号，k的取值为1，2，…，m。

2.根据权利要求1所述的基于变换域广义伽玛分布的煤岩识别方法，其特征在于，所述N级两路小波变换包括以下步骤：

(1)初始化设置两路小波变换的级数序号i＝1；

(2)彩色图像灰度化处理，即把24位煤岩彩色图像转换成8位煤岩灰度图像，采用的转换公式为Y＝0.299R+0.587G+0.114B，其中R，G和B分别表示转换前彩色图像的红色，绿色和蓝色分量，Y表示转换后灰度图像的像素值；

(3)用一维的低通滤波器h₀对煤岩灰度图像进行逐列滤波，从而得到低频部分小波系数子带Coef_Lo；

(4)用一维的高通滤波器h₁对煤岩灰度图像进行逐列滤波，从而得到高频部分小波系数子带Coef_Hi；

(5)用步骤(3)所述的低通滤波器h₀对步骤(4)所述的高频部分小波系数子带Coef_Hi进行逐行滤波，从而得到高频部分小波系数子带Coef_HL；

(6)用步骤(4)所述的高通滤波器h₁对步骤(3)所述的低频部分小波系数子带Coef_Lo进行逐行滤波，从而得到高频部分小波系数子带Coef_LH；

(7)用步骤(4)所述的高通滤波器h₁对步骤(4)所述的高频部分小波系数子带Coef_Hi进行逐行滤波，从而得到高频部分小波系数子带Coef_HH；

(8)用步骤(3)所述的低通滤波器h₀对步骤(3)所述的低频部分小波系数子带Coef_Lo进行逐行滤波，从而得到低频部分小波系数子带Coef_LL；

(9)利用步骤(5)所述的高频部分小波系数子带Coef_HL，分别构造煤岩图像在经过第1级两路小波变换以后的第1个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第1个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第6个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第6个方向的第2路高频部分小波系数子带

(10)利用步骤(6)所述的高频部分小波系数子带Coef_LH，分别构造煤岩图像在经过第 1级两路小波变换以后的第3个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第3个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第4个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第4个方向的第2路高频部分小波系数子带

(11)利用步骤(7)所述的高频部分小波系数子带Coef_HH，分别构造煤岩图像在经过第1级两路小波变换以后的第2个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第2个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第5个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第5个方向的第2路高频部分小波系数子带

(12)两路小波变换的级数序号自增1，即i＝i+1；

(13)动态申请与Coef_LL数据结构完全一致的临时数据存储区TempCoef，并将其数据内容初始化为Coef_LL的数据内容；

(14)判断i≤N条件是否成立，如果是，则进入以下步骤(15)-(28)的迭代循环，如果否，则转到步骤(29)；

(15)如果TempCoef数据的总行数不是4的倍数，那么对TempCoef数据进行修正，在TempCoef数据的第一行Row_first前面追加一行Row′_first并且Row′_first的数据内容用Row_first的数据内容进行填充，在TempCoef数据的最后一行Row_last后面追加一行Row′_last并且Row′_last的数据内容用Row_last的数据内容进行填充；

(16)如果TempCoef数据的总列数不是4的倍数，那么对TempCoef数据进行修正，在TempCoef数据的第一列Col_first前面追加一列Col′_first并且Col′_first的数据内容用Col_first的数据内容进行填充，在TempCoef数据的最后一列Col_last后面追加一列Col′_last并且Col′_last的数据内容用Col_last的数据内容进行填充；

(17)用两组一维的低通滤波器h₀₀、h₀₁对TempCoef数据进行双重列滤波，从而得到低频部分小波系数子带Coef_L2-i，其中下标i表示两路小波变换的级数序号；

(18)用两组一维的高通滤波器h₁₀、h₁₁对TempCoef数据进行双重列滤波，从而得到高频部分小波系数子带Coef_H2-i，其中下标i表示两路小波变换的级数序号；

(19)用步骤(17)所述的低通滤波器h₀₀、h₀₁对步骤(18)所述的高频部分小波系数子带Coef_H2-i进行双重行滤波，从而得到高频部分小波系数子带Coef_HL2-i，其中下标i表示两 路小波变换的级数序号；

(20)用步骤(18)所述的高通滤波器h₁₀、h₁₁对步骤(17)所述的低频部分小波系数子带Coef_L2-i进行双重行滤波，从而得到高频部分小波系数子带Coef_LH2-i，其中下标i表示两路小波变换的级数序号；

(21)用步骤(18)所述的高通滤波器h₁₀、h₁₁对步骤(18)所述的高频部分小波系数子带Coef_H2-i进行双重行滤波，从而得到高频部分小波系数子带Coef_HH2-i，其中下标i表示两路小波变换的级数序号；

(22)用步骤(17)所述的低通滤波器h₀₀、h₀₁对步骤(17)所述的低频部分小波系数子带Coef_L2-i进行双重行滤波，从而得到低频部分小波系数子带Coef_LL2-i，其中下标i表示两路小波变换的级数序号；

(23)利用步骤(19)所述的高频部分小波系数子带Coef_HL2-i，分别构造煤岩图像在经过第i级两路小波变换以后的第1个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第1个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第6个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第6个方向的第2路高频部分小波系数子带

(24)利用步骤(20)所述的高频部分小波系数子带Coef_LH2-i，分别构造煤岩图像在经过第i级两路小波变换以后的第3个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第3个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第4个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第4个方向的第2路高频部分小波系数子带

(25)利用步骤(21)所述的高频部分小波系数子带Coef_HH2-i，分别构造煤岩图像在经过第i级两路小波变换以后的第2个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第2个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第5个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第5个方向的第2路高频部分小波系数子带

(26)释放临时数据存储区TempCoef所占用的内存空间，重新动态申请与步骤(22)所述Coef_LL2-i的数据结构完全一致的临时数据存储区TempCoef，并将其数据内容初始化为Coef_LL2-i的数据内容；

(27)两路小波变换的级数序号自增1，即i＝i+1；

(28)转到步骤(14)，进行下一级两路小波变换；

(29)释放临时数据存储区TempCoef所占用的内存空间，完成N级两路小波变换。

3.根据权利要求1所述的基于变换域广义伽玛分布的煤岩识别方法，其特征在于，所述计算均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ包括以下步骤：

STEP1.把τ₀＝0.05作为τ参数的初始迭代值，进入以下步骤STEP2-STEP4的循环迭代过程；

STEP2.把τ作为自变量代入并计算

其中x_i表示均方根小波系数子带中的第i个均方根小波系数，Num表示均方根小波系数子带中的均方根小波系数总数，ln表示以自然常数e为底数的对数函数，即自然对数；

STEP3.把τ作为自变量代入并计算

其中x_i表示均方根小波系数子带中的第i个均方根小波系数，Num表示均方根小波系数子带中的均方根小波系数总数，ln表示以自然常数e为底数的对数函数，即自然对数；

STEP4.计算若|Δτ|大于计算精度ε，则通过公式τ＝τ-|Δτ|来更新τ，然后执行步骤STEP2-STEP4的循环迭代过程，若|Δτ|小于或等于计算精度ε，则把当前τ值作为最终τ值并进入步骤STEP5；

STEP5.根据步骤STEP4获得的最终τ值，计算得到最终的α值和λ值，最终的α值，λ值的计算公式分别为

其中x_i表示均方根小波系数子带中的第i个均方根小波系数，Num表示均方根小波系数子带中的均方根小波系数总数，ln表示以自然常数e为底数的对数函数，即自然对数。

4.根据权利要求1所述的基于变换域广义伽玛分布的煤岩识别方法，其特征在于，所述6×N×3维的特征列向量的数据结构采用[α_1-1，τ_1-1，λ_1-1，α_1-2，τ_1-2，λ_1-2，α_1-3，τ_1-3，λ_1-3，α_1-4，τ_1-4，λ_1-4，α_1-5，τ_1-5，λ_1-5，α_1-6，τ_1-6，λ_1-6，α_2-1，τ_2-1，λ_2-1，α_2-2，τ_2-2，λ_2-2，α_2-3，τ_2-3，λ_2-3，α_2-4，τ_2-4，λ_2-4，α_2-5，τ_2-5，λ_2-5，α_2-6，τ_2-6，λ_2-6，…，α_N-1，τ_N-1，λ_N-1，α_N-2，τ_N-2，λ_N-2，α_N-3，τ_N-3，λ_N-3，α_N-4，τ_N-4，λ_N-4，α_N-5，τ_N-5，λ_N-5，α_N-6，τ_N-6，λ_N-6]^T的格式，特征列向量中的每一个α元素、τ元素、λ元素可以分别表示为α_i-n、τ_i-n、λ_i-n，其中下标i表示两路小波变换的级数序号，i的取值为1，2，…，N，n表示均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排列的排列序号，n的取值为1，2，…，6，α_i-n、τ_i-n和λ_i-n分别表示煤岩样本子图在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ。

5.根据权利要求1所述的基于变换域广义伽玛分布的煤岩识别方法，其特征在于，所述特征列向量与特征列向量之间的相似度，以及特征列向量与特征列向量之间的相似度计算公式分别为

N表示两路小波变换的总级数，和中的每一个元素可以分别表示为ux_p、ck_p和rk_p，p的取值为1，2，…，(6×N×3)，ux_p表示的第p个元素，ck_p表示的第p个元素，rk_p表示的第p个元素，函数在Г(z)函数表达式中，z为自变量，t为积分变量，积分区间为[0，+∞)，e表示自然常数，函数在ψ(z)函数表达式中，z为自变量，Г′(z)为Г(z)的导函数。

6.根据权利要求2所述的基于变换域广义伽玛分布的煤岩识别方法，其特征在于，所述一维的低通滤波器h₀和高通滤波器h₁的序列长度均为奇数，并且h₀序列和h₁序列均为中心对称序列。

7.根据权利要求2所述的基于变换域广义伽玛分布的煤岩识别方法，其特征在于，所述一维的低通滤波器h₀₀、h₀₁和高通滤波器h₁₀、h₁₁的序列长度相等并且均为偶数，h₀₁序列为h₀₀序列的反转序列，h₁₀序列通过对h₀₀序列进行反转操作然后对偶数序号采样点的序列元素进行符号取反操作而得到，h₁₁序列为h₁₀序列的反转序列。

8.根据权利要求2所述的基于变换域广义伽玛分布的煤岩识别方法，其特征在于，在给定两组序列长度相等的一维滤波器h_f、h_s和滤波对象X的前提下，所述双重列滤波包括以下步骤：

STC1.对X进行行对称延拓得到X_ex，如果X的总行数和总列数分别为RowsX和ColsX，h_f和h_s的序列长度为len，那么X_ex的总行数和总列数分别为(RowsX+2×len)和ColsX，并且X_ex的第1，2，…，len行的数据内容分别用X的第len，len-1，…，1行的数据内容进行填充，X_ex的第len+1，len+2，…，len+RowsX行的数据内容分别用X的第1，2，…，RowsX行的数据内容进行填充，Xex的第len+RowsX+1，len+RowsX+2，…，len+RowsX+len行的数据内容分别用X的第RowsX，RowsX-1，…，RowsX-(len-1)行的数据内容进行填充；

STC2.抽取h_f中序号为奇数的滤波系数构成子滤波器h_fo，抽取h_f中序号为偶数的滤波系数构成子滤波器h_fe，抽取h_s中序号为奇数的滤波系数构成子滤波器h_so，抽取h_s中序号为 偶数的滤波系数构成子滤波器h_se；

STC3.对步骤STC1所述的X_ex进行修正，即舍弃X_ex的起始两行和末尾两行，从而得到X′_ex；

STC4.抽取步骤STC3所述X′_ex中行号为奇数的数据构成X′_exO，抽取步骤STC3所述X′_ex中行号为偶数的数据构成X′_exE；

STC5.抽取步骤STC4所述X′_exE中行号为偶数的数据构成X′_exA，抽取步骤STC4所述X′_exO中行号为偶数的数据构成X′_exB，抽取步骤STC4所述X′_exE中行号为奇数的数据构成X′_exC，抽取步骤STC4所述X′_exO中行号为奇数的数据构成X′_exD；

STC6.用步骤STC2所述的子滤波器h_fo对步骤STC5所述的X′_exB进行逐列滤波，从而得到ResB；

STC7.用步骤STC2所述的子滤波器h_fe对步骤STC5所述的X′_exD进行逐列滤波，从而得到ResD；

STC8.用步骤STC2所述的子滤波器h_so对步骤STC5所述的X′_exA进行逐列滤波，从而得到ResA；

STC9.用步骤STC2所述的子滤波器h_se对步骤STC5所述的X′_exC进行逐列滤波，从而得到ResC；

STC10.把步骤STC6所述的ResB和步骤STC7所述的ResD中行号相同且列号也相同的元素相加得到ResBD，即ResBD＝ResB+ResD；

STC11.把步骤STC8所述的ResA和步骤STC9所述的ResC中行号相同且列号也相同的元素相加得到ResAC，即ResAC＝ResA+ResC；

STC12.申请一块总行数为(RowsAC+RowsBD)且总列数为ColsAC或ColsBD的数据存储区FinalRes，其中RowsAC表示步骤STC11所述ResAC的总行数，RowsBD表示步骤STC10所述ResBD的总行数，ColsAC表示步骤STC11所述ResAC的总列数，ColsBD表示步骤STC10所述ResBD的总列数，并且满足RowsAC＝RowsBD和ColsAC＝ColsBD的条件；

STC13.把一维滤波器h_f和h_s视为两个向量，如果这两个向量的内积大于0，那么步骤STC12所述FinalRes的第1，3，5，7，…，(RowsAC+RowsBD-1)行的数据内容分别用步骤STC10所述ResBD的第1，2，3，4，…，RowsBD行的数据内容进行填充，步骤STC12所述FinalRes的第2，4，6，8，…，(RowsAC+RowsBD)行的数据内容分别用步骤STC11所述ResAC的第1，2，3，4，…，RowsAC行的数据内容进行填充，如果这两个向量的内积小于 或等于0，那么步骤STC12所述FinalRes的第1，3，5，7，…，(RowsAC+RowsBD-1)行的数据内容分别用步骤STC11所述ResAC的第1，2，3，4，…，RowsAC行的数据内容进行填充，步骤STC12所述FinalRes的第2，4，6，8，…，(RowsAC+RowsBD)行的数据内容分别用步骤STC10所述ResBD的第1，2，3，4，…，RowsBD行的数据内容进行填充，其中RowsAC表示步骤STC11所述ResAC的总行数，RowsBD表示步骤STC10所述ResBD的总行数；

STC14.完成双重列滤波操作，其最终结果输出为步骤STC12所述的FinalRes。

9.根据权利要求2所述的基于变换域广义伽玛分布的煤岩识别方法，其特征在于，在给定两组序列长度相等的一维滤波器h_f、h_s和滤波对象X的前提下，所述双重行滤波包括以下步骤：

STR1.对X进行列对称延拓得到X_ex，如果X的总行数和总列数分别为RowsX和ColsX，h_f和h_s的序列长度为len，那么X_ex的总行数和总列数分别为RowsX和(ColsX+2×len)，并且X_ex的第1，2，…，len列的数据内容分别用X的第len，len-1，…，l列的数据内容进行填充，X_ex的第len+1，len+2，…，len+ColsX列的数据内容分别用X的第1，2，…，ColsX列的数据内容进行填充，X_ex的第len+ColsX+1，len+ColsX+2，…，len+ColsX+len列的数据内容分别用X的第ColsX，ColsX-1，…，ColsX-(len-1)列的数据内容进行填充；

STR2.抽取h_f中序号为奇数的滤波系数构成子滤波器h_fo，抽取h_f中序号为偶数的滤波系数构成子滤波器h_fe，抽取h_s中序号为奇数的滤波系数构成子滤波器h_so，抽取h_s中序号为偶数的滤波系数构成子滤波器h_se；

STR3.对步骤STR1所述的X_ex进行修正，即舍弃X_ex的起始两列和末尾两列，从而得到X′_ex；

STR4.抽取步骤STR3所述X′_ex中列号为奇数的数据构成X′_exO，抽取步骤STR3所述X′_ex中列号为偶数的数据构成X′_exE；

STR5.抽取步骤STR4所述X′_exE中列号为偶数的数据构成X′_exA，抽取步骤STR4所述X′_exO中列号为偶数的数据构成X′_exB，抽取步骤STR4所述X′_exE中列号为奇数的数据构成X′_exC，抽取步骤STR4所述X′_exO中列号为奇数的数据构成X′_exD；

STR6.用步骤STR2所述的子滤波器h_fo对步骤STR5所述的X′_exB进行逐行滤波，从而得到ResB；

STR7.用步骤STR2所述的子滤波器h_fe对步骤STR5所述的X′_exD进行逐行滤波，从而得 到ResD；

STR8.用步骤STR2所述的子滤波器h_so对步骤STR5所述的X′_exA进行逐行滤波，从而得到ResA；

STR9.用步骤STR2所述的子滤波器h_se对步骤STR5所述的X′_exC进行逐行滤波，从而得到ResC；

STR10.把步骤STR6所述的ResB和步骤STR7所述的ResD中行号相同且列号也相同的元素相加得到ResBD，即ResBD＝ResB+ResD；

STR11.把步骤STR8所述的ResA和步骤STR9所述的ResC中行号相同且列号也相同的元素相加得到ResAC，即ResAC＝ResA+ResC；

STR12.申请一块总行数为RowsAC或RowsBD且总列数为(ColsAC+ColsBD)的数据存储区FinalRes，其中RowsAC表示步骤STR11所述ResAC的总行数，RowsBD表示步骤STR10所述ResBD的总行数，ColsAC表示步骤STR11所述ResAC的总列数，ColsBD表示步骤STR10所述ResBD的总列数，并且满足RowsAC＝RowsBD和ColsAC＝ColsBD的条件；

STR13.把一维滤波器h_f和h_s视为两个向量，如果这两个向量的内积大于0，那么步骤STR12所述FinalRes的第1，3，5，7，…，(ColsAC+ColsBD-1)列的数据内容分别用步骤STR10所述ResBD的第1，2，3，4，…，ColsBD列的数据内容进行填充，步骤STR12所述FinalRes的第2，4，6，8，…，(ColsAC+ColsBD)列的数据内容分别用步骤STR11所述ResAC的第1，2，3，4，…，ColsAC列的数据内容进行填充，如果这两个向量的内积小于或等于0，那么步骤STR12所述FinalRes的第1，3，5，7，…，(ColsAC+ColsBD-1)列的数据内容分别用步骤STR11所述ResAC的第1，2，3，4，…，ColsAC列的数据内容进行填充，步骤STR12所述FinalRes的第2，4，6，8，…，(ColsAC+ColsBD)列的数据内容分别用步骤STR10所述ResBD的第1，2，3，4，…，ColsBD列的数据内容进行填充，其中ColsAC表示步骤STR11所述ResAC的总列数，ColsBD表示步骤STR10所述ResBD的总列数；

STR14.完成双重行滤波操作，其最终结果输出为步骤STR12所述的FinalRes。