CN104966074A - 基于变换域广义伽玛分布的煤岩识别方法 - Google Patents
基于变换域广义伽玛分布的煤岩识别方法 Download PDFInfo
- Publication number
- CN104966074A CN104966074A CN201510418701.1A CN201510418701A CN104966074A CN 104966074 A CN104966074 A CN 104966074A CN 201510418701 A CN201510418701 A CN 201510418701A CN 104966074 A CN104966074 A CN 104966074A
- Authority
- CN
- China
- Prior art keywords
- wavelet coefficient
- coefficient subband
- row
- hfs
- mean square
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 239000003245 coal Substances 0.000 title claims abstract description 144
- 239000011435 rock Substances 0.000 title claims abstract description 56
- 238000000034 method Methods 0.000 title claims abstract description 51
- 238000009826 distribution Methods 0.000 title claims abstract description 37
- 230000009466 transformation Effects 0.000 claims abstract description 164
- 239000013598 vector Substances 0.000 claims abstract description 65
- 238000012549 training Methods 0.000 claims abstract description 10
- 238000012545 processing Methods 0.000 claims abstract description 6
- 238000001914 filtration Methods 0.000 claims description 92
- 238000000605 extraction Methods 0.000 claims description 36
- 230000006870 function Effects 0.000 claims description 33
- 230000009977 dual effect Effects 0.000 claims description 29
- 230000015572 biosynthetic process Effects 0.000 claims description 24
- 239000000284 extract Substances 0.000 claims description 18
- 102100030510 Stanniocalcin-2 Human genes 0.000 claims description 15
- 238000004364 calculation method Methods 0.000 claims description 12
- 230000008569 process Effects 0.000 claims description 11
- 238000006243 chemical reaction Methods 0.000 claims description 9
- 238000013500 data storage Methods 0.000 claims description 9
- 102100030511 Stanniocalcin-1 Human genes 0.000 claims description 6
- 238000005286 illumination Methods 0.000 claims description 6
- 230000010354 integration Effects 0.000 claims description 6
- 230000000052 comparative effect Effects 0.000 claims description 5
- 238000009434 installation Methods 0.000 claims description 5
- 101100279959 Gibberella fujikuroi (strain CBS 195.34 / IMI 58289 / NRRL A-6831) STC3 gene Proteins 0.000 claims description 3
- 101100393982 Gibberella fujikuroi (strain CBS 195.34 / IMI 58289 / NRRL A-6831) STC5 gene Proteins 0.000 claims description 3
- 101000701440 Homo sapiens Stanniocalcin-1 Proteins 0.000 claims description 3
- 101000701446 Homo sapiens Stanniocalcin-2 Proteins 0.000 claims description 3
- 101150108015 STR6 gene Proteins 0.000 claims description 3
- 101100386054 Saccharomyces cerevisiae (strain ATCC 204508 / S288c) CYS3 gene Proteins 0.000 claims description 3
- 101150035983 str1 gene Proteins 0.000 claims description 3
- 238000012360 testing method Methods 0.000 abstract 5
- 239000000523 sample Substances 0.000 description 30
- 238000001514 detection method Methods 0.000 description 4
- 238000004519 manufacturing process Methods 0.000 description 4
- 238000005065 mining Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 230000004069 differentiation Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000005251 gamma ray Effects 0.000 description 1
- 238000002309 gasification Methods 0.000 description 1
- 238000003384 imaging method Methods 0.000 description 1
- 238000004886 process control Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
- G06F18/2134—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on separation criteria, e.g. independent component analysis
- G06F18/21347—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on separation criteria, e.g. independent component analysis using domain transformations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/0002—Inspection of images, e.g. flaw detection
- G06T7/0004—Industrial image inspection
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
- G06F2218/02—Preprocessing
- G06F2218/04—Denoising
- G06F2218/06—Denoising by applying a scale-space analysis, e.g. using wavelet analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10004—Still image; Photographic image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20048—Transform domain processing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20048—Transform domain processing
- G06T2207/20064—Wavelet transform [DWT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30108—Industrial image inspection
- G06T2207/30132—Masonry; Concrete
Landscapes
- Engineering & Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Artificial Intelligence (AREA)
- Life Sciences & Earth Sciences (AREA)
- Quality & Reliability (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Compression Or Coding Systems Of Tv Signals (AREA)
- Image Processing (AREA)
Abstract
本发明公开了基于变换域广义伽玛分布的煤岩识别方法。在样本训练阶段,获取满足一定条件的多幅煤炭训练样本图像和多幅岩石训练样本图像,对每幅图像进行两路小波变换,经过相关处理得到按一定排序规则排列的均方根小波系数子带,计算这些均方根小波系数子带在服从广义伽玛分布条件下的各参数,并利用这些参数构造训练样本特征列向量;在煤岩识别阶段,采集满足一定条件的未知类别测试样本图像,与训练样本图像的处理方式类似,对测试样本图像进行一系列处理,最终得到测试样本特征列向量。通过比较测试样本特征列向量与训练样本特征列向量之间的相似度判断测试样本所属的煤岩类型。本发明可靠性好,正确识别率高,软硬件实现容易。
Description
技术领域
本发明涉及基于变换域广义伽玛分布的煤岩识别方法,属于图像识别技术领域。
背景技术
煤岩识别是指通过各种技术手段自动判别出煤炭与岩石。在煤炭资源开采及运输过程中,存在许多生产环节需要判别区分煤炭与岩石,如采煤机滚筒高度调节、综采放顶煤过程控制、选煤厂原煤选矸等。从20世纪50年代开始,南非、澳大利亚、德国、美国、中国等世界主要产煤国家对煤岩识别方法展开了一系列研究,相继产生了一些代表性的研究成果,如自然γ射线探测法、雷达探测法、红外探测法、有功功率检测法、振动信号检测法、声音信号检测法等。然而这些方法均存在以下共性问题:(1)需要在现有设备上安装部署各种传感器、导致装置结构复杂,制造成本高;(2)采煤机、掘进机等机械设备在煤炭生产过程中受力复杂、振动剧烈、磨损严重,传感器部署相对比较困难,其电子线路也容易受到损坏,装置可靠性差;(3)针对不同类型的机械载体设备,传感器的选型和安装位置的选择存在较大区别,这就需要进行个性化定制,因此其普适性不佳。
通过对块状的煤炭、岩石样本的观察,发现煤炭与岩石在颜色、光泽、纹理等方面存在较大的差异。当通过现有的数字摄像设备对煤炭或岩石进行成像时,煤炭与岩石的视觉差异信息就隐藏在采集得到的数字图像中了,因此提出通过挖掘煤岩数字图像中的视觉信息来区分煤炭与岩石。现有的基于图像处理的煤岩识别方法在鲁棒性、识别率等方面还存在着较大的提升空间。
发明内容
为了克服现有煤岩识别方法存在的不足,本发明提出基于变换域广义伽玛分布的煤岩识别方法,该方法具有实时性强、识别率高、稳健性好等优点,能够为现代化煤矿安全高效生产提供有力的技术保障。
本发明所述的煤岩识别方法采用如下技术方案实现,包括样本训练阶段和煤岩识别阶段,具体步骤如下:
A.在样本训练阶段,获取相同光照条件下拍摄的m幅已知煤炭样本图像和m幅已知岩石样本图像,截取宽度和高度均为偶数个像素点并且不含非煤岩背景的子图像,煤炭样本子图和岩石样本子图分别记为c1,c2,…,cm和r1,r2,…,rm;
B.分别对样本子图c1,c2,…,cm和r1,r2,…,rm进行N级两路小波变换,每幅样本子图得到12×N个高频部分小波系数子带,每个高频部分小波系数子带记为 表示样本子图在经过第i级两路小波变换以后的第d个方向的第j路高频部分小波系数子带,其中i表示两路小波变换的级数序号,i的取值为1,2,…,N,d表示两路小波变换的方向序号,d的取值为1,2,…,6,j表示两路小波变换的路序号,j的取值为1,2;
C.对步骤B所述的每幅样本子图的级数序号相同并且方向序号也相同的两个高频部分小波系数子带进行求均方根操作,得到6×N个均方根小波系数子带,每一个均方根小波系数子带记为RmsCoefi-d,RmsCoefi-d表示样本子图在经过第i级两路小波变换以后的第d个方向的均方根小波系数子带;
D.分别计算步骤C所述的每一个均方根小波系数子带RmsCoefi-d的均值μi-d和方差
E.根据步骤D所计算的每一个均方根小波系数子带的均值和方差,对步骤C所述的每一级两路小波变换以后的6个方向的均方根小波系数子带按照均值和方差的乘积值从大到小的顺序排列,每一个经过排序以后的均方根小波系数子带记为RmsCoefSi-n,RmsCoefSi-n表示样本子图在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带,n表示均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排列的排列序号,n的取值为1,2,…,6;
F.分别计算步骤E所述的经过排序以后的每一个均方根小波系数子带在服从概率密度函数为的广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ,其中x表示均方根小波系数子带中的元素,也就是均方根小波系数,函数在Γ(α)函数表达式中,α为自变量,t为积分变量,积分区间为[0,+∞),e表示自然常数;
G.利用每一幅煤炭样本子图ck或岩石样本子图rk经过N级两路小波变换、高频部分小波系数子带求均方根操作并且按上述步骤E排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ,构造出一个6×N×3维的特征列向量或其中下标k表示样本序号,k的取值为1,2,…,m,然后把这2×m个特征列向量保存到数据存储设备中,用于后续的煤岩识别阶段;
H.在煤岩识别阶段,获取相同光照条件下拍摄的未知类别煤岩图像,截取宽度和高度均为偶数个像素点并且不含非煤岩背景的待识别子图ux;
I.与上述步骤B类似,对ux进行N级两路小波变换,从而得到12×N个高频部分小波系数子带,每个高频部分小波系数子带记为 表示待识别子图ux在经过第i级两路小波变换以后的第d个方向的第j路高频部分小波系数子带;
J.与上述步骤C类似,对步骤I所述的待识别子图ux的级数序号相同并且方向序号也相同的两个高频部分小波系数子带进行求均方根操作,得到6×N个均方根小波系数子带,每一个均方根小波系数子带记为RmsUxCfi-d,RmsUxCfi-d表示待识别子图ux在经过第i级两路小波变换以后的第d个方向的均方根小波系数子带;
K.与上述步骤D类似,分别计算步骤J所述的每一个均方根小波系数子带RmsUxCfi-d的均值和方差
L.与上述步骤E类似,根据步骤K所计算的每一个均方根小波系数子带的均值和方差,对步骤J所述的每一级两路小波变换以后的6个方向的均方根小波系数子带按照均值和方差的乘积值从大到小的顺序排列,每一个经过排序以后的均方根小波系数子带记为RmsUxCfSi-n,RmsUxCfSi-n表示待识别子图ux在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带;
M.与上述步骤F类似,分别计算步骤L所述的经过排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ;
N.与上述步骤G类似,利用待识别子图ux经过N级两路小波变换、高频部分小波系数子带求均方根操作并且按上述步骤L排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ,构造出一个6×N×3维的特征列向量
O.比较特征列向量与特征列向量之间的相似度,其中下标k表示样本序号,k的取值为1,2,…,m,从而判定待识别子图ux所属的煤岩类别。
所述N级两路小波变换包括以下步骤:
(1)初始化设置两路小波变换的级数序号i=1;
(2)彩色图像灰度化处理,即把24位煤岩彩色图像转换成8位煤岩灰度图像,采用的转换公式为Y=0.299R+0.587G+0.114B,其中R,G和B分别表示转换前彩色图像的红色,绿色和蓝色分量,Y表示转换后灰度图像的像素值;
(3)用一维的低通滤波器h0对煤岩灰度图像进行逐列滤波,从而得到低频部分小波系数子带CoefLo;
(4)用一维的高通滤波器h1对煤岩灰度图像进行逐列滤波,从而得到高频部分小波系数子带CoefHi;
(5)用步骤(3)所述的低通滤波器h0对步骤(4)所述的高频部分小波系数子带CoefHi进行逐行滤波,从而得到高频部分小波系数子带CoefHL;
(6)用步骤(4)所述的高通滤波器h1对步骤(3)所述的低频部分小波系数子带CoefLo进行逐行滤波,从而得到高频部分小波系数子带CoefLH;
(7)用步骤(4)所述的高通滤波器h1对步骤(4)所述的高频部分小波系数子带CoefHi进行逐行滤波,从而得到高频部分小波系数子带CoefHH;
(8)用步骤(3)所述的低通滤波器h0对步骤(3)所述的低频部分小波系数子带CoefLo进行逐行滤波,从而得到低频部分小波系数子带CoefLL;
(9)利用步骤(5)所述的高频部分小波系数子带CoefHL,分别构造煤岩图像在经过第1级两路小波变换以后的第1个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第1个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第6个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第6个方向的第2路高频部分小波系数子带
(10)利用步骤(6)所述的高频部分小波系数子带CoefLH,分别构造煤岩图像在经过第1级两路小波变换以后的第3个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第3个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第4个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第4个方向的第2路高频部分小波系数子带
(11)利用步骤(7)所述的高频部分小波系数子带CoefHH,分别构造煤岩图像在经过第1级两路小波变换以后的第2个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第2个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第5个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第5个方向的第2路高频部分小波系数子带
(12)两路小波变换的级数序号自增1,即i=i+1;
(13)动态申请与CoefLL数据结构完全一致的临时数据存储区TempCoef,并将其数据内容初始化为CoefLL的数据内容;
(14)判断i≤N条件是否成立,如果是,则进入以下步骤(15)-(28)的迭代循环,如果否,则转到步骤(29);
(15)如果TempCoef数据的总行数不是4的倍数,那么对TempCoef数据进行修正,在TempCoef数据的第一行Rowfirst前面追加一行Row′first并且Row′first的数据内容用Rowfirst的数据内容进行填充,在TempCoef数据的最后一行Rowlast后面追加一行Row′last并且Row′last的数据内容用Rowlast的数据内容进行填充;
(16)如果TempCoef数据的总列数不是4的倍数,那么对TempCoef数据进行修正,在TempCoef数据的第一列Colfirst前面追加一列Col′first并且Col′first的数据内容用Colfirst的数据内容进行填充,在TempCoef数据的最后一列Collast后面追加一列Col′last并且Col′last的数据内容用Collast的数据内容进行填充;
(17)用两组一维的低通滤波器h00、h01对TempCoef数据进行双重列滤波,从而得到低频部分小波系数子带CoefL2-i,其中下标i表示两路小波变换的级数序号;
(18)用两组一维的高通滤波器h10、h11对TempCoef数据进行双重列滤波,从而得到高频部分小波系数子带CoefH2-i,其中下标i表示两路小波变换的级数序号;
(19)用步骤(17)所述的低通滤波器h00、h01对步骤(18)所述的高频部分小波系数子带CoefH2-i进行双重行滤波,从而得到高频部分小波系数子带CoefHL2-i,其中下标i表示两路小波变换的级数序号;
(20)用步骤(18)所述的高通滤波器h10、h11对步骤(17)所述的低频部分小波系数子带CoefL2-i进行双重行滤波,从而得到高频部分小波系数子带CoefLH2-i,其中下标i表示两路小波变换的级数序号;
(21)用步骤(18)所述的高通滤波器h10、h11对步骤(18)所述的高频部分小波系数子带CoefH2-i进行双重行滤波,从而得到高频部分小波系数子带CoefHH2-i,其中下标i表示两路小波变换的级数序号;
(22)用步骤(17)所述的低通滤波器h00、h01对步骤(17)所述的低频部分小波系数子带CoefL2-i进行双重行滤波,从而得到低频部分小波系数子带CoefLL2-i,其中下标i表示两路小波变换的级数序号;
(23)利用步骤(19)所述的高频部分小波系数子带CoefHL2-i,分别构造煤岩图像在经过第i级两路小波变换以后的第1个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第1个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第6个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第6个方向的第2路高频部分小波系数子带
(24)利用步骤(20)所述的高频部分小波系数子带CoefLH2-i,分别构造煤岩图像在经过第i级两路小波变换以后的第3个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第3个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第4个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第4个方向的第2路高频部分小波系数子带
(25)利用步骤(21)所述的高频部分小波系数子带CoefHH2-i,分别构造煤岩图像在经过第i级两路小波变换以后的第2个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第2个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第5个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第5个方向的第2路高频部分小波系数子带
(26)释放临时数据存储区TempCoef所占用的内存空间,重新动态申请与步骤(22)所述CoefLL2-i的数据结构完全一致的临时数据存储区TempCoef,并将其数据内容初始化为CoefLL2-i的数据内容;
(27)两路小波变换的级数序号自增1,即i=i+1;
(28)转到步骤(14),进行下一级两路小波变换;
(29)释放临时数据存储区TempCoef所占用的内存空间,完成N级两路小波变换。
所述一维的低通滤波器h0和高通滤波器h1的序列长度均为奇数,并且h0序列和h1序列均为中心对称序列。
所述一维的低通滤波器h00、h01和高通滤波器h10、h11的序列长度相等并且均为偶数,h01序列为h00序列的反转序列,h10序列通过对h00序列进行反转操作然后对偶数序号采样点的序列元素进行符号取反操作而得到,h11序列为h10序列的反转序列。
在给定两组序列长度相等的一维滤波器hf、hs和滤波对象X的前提下,所述双重列滤波包括以下步骤:
STC1.对X进行行对称延拓得到Xex,如果X的总行数和总列数分别为RowsX和ColsX,hf和hs的序列长度为len,那么Xex的总行数和总列数分别为(RowsX+2×len)和ColsX,并且Xex的第1,2,…,len行的数据内容分别用X的第len,len-1,…,1行的数据内容进行填充,Xex的第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.抽取hf中序号为奇数的滤波系数构成子滤波器hfo,抽取hf中序号为偶数的滤波系数构成子滤波器hfe,抽取hs中序号为奇数的滤波系数构成子滤波器hso,抽取hs中序号为偶数的滤波系数构成子滤波器hse;
STC3.对步骤STC1所述的Xex进行修正,即舍弃Xex的起始两行和末尾两行,从而得到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所述的子滤波器hfo对步骤STC5所述的X′exB进行逐列滤波,从而得到ResB;
STC7.用步骤STC2所述的子滤波器hfe对步骤STC5所述的X′exD进行逐列滤波,从而得到ResD;
STC8.用步骤STC2所述的子滤波器hso对步骤STC5所述的X′exA进行逐列滤波,从而得到ResA;
STC9.用步骤STC2所述的子滤波器hse对步骤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.把一维滤波器hf和hs视为两个向量,如果这两个向量的内积大于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。
在给定两组序列长度相等的一维滤波器hf、hs和滤波对象X的前提下,所述双重行滤波包括以下步骤:
STR1.对X进行列对称延拓得到Xex,如果X的总行数和总列数分别为RowsX和ColsX,hf和hs的序列长度为len,那么Xex的总行数和总列数分别为RowsX和(ColsX+2×len),并且Xex的第1,2,…,len列的数据内容分别用X的第len,len-1,…,1列的数据内容进行填充,Xex的第len+1,len+2,…,len+ColsX列的数据内容分别用X的第1,2,…,ColsX列的数据内容进行填充,Xex的第len+ColsX+1,len+ColsX+2,…,len+ColsX+len列的数据内容分别用X的第ColsX,ColsX-1,…,ColsX-(len-1)列的数据内容进行填充;
STR2.抽取hf中序号为奇数的滤波系数构成子滤波器hfo,抽取hf中序号为偶数的滤波系数构成子滤波器hfe,抽取hs中序号为奇数的滤波系数构成子滤波器hso,抽取hs中序号为偶数的滤波系数构成子滤波器hse;
STR3.对步骤STR1所述的Xex进行修正,即舍弃Xex的起始两列和末尾两列,从而得到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所述的子滤波器hfo对步骤STR5所述的X′exB进行逐行滤波,从而得到ResB;
STR7.用步骤STR2所述的子滤波器hfe对步骤STR5所述的X′exD进行逐行滤波,从而得到ResD;
STR8.用步骤STR2所述的子滤波器hso对步骤STR5所述的X′exA进行逐行滤波,从而得到ResA;
STR9.用步骤STR2所述的子滤波器hse对步骤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.把一维滤波器hf和hs视为两个向量,如果这两个向量的内积大于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=0.05作为τ参数的初始迭代值,进入以下步骤STEP2-STEP4的循环迭代过程;
STEP2.把τ作为自变量代入并计算
其中xi表示均方根小波系数子带中的第i个均方根小波系数,Num表示均方根小波系数子带中的均方根小波系数总数,ln表示以自然常数e为底数的对数函数,即自然对数;
STEP3.把τ作为自变量代入并计算
其中xi表示均方根小波系数子带中的第i个均方根小波系数,Num表示均方根小波系数子带中的均方根小波系数总数,ln表示以自然常数e为底数的对数函数,即自然对数;
STEP4.计算若|Δτ|大于计算精度ε,则通过公式τ=τ-|Δτ|来更新τ,然后执行步骤STEP2-STEP4的循环迭代过程,若|Δτ|小于或等于计算精度ε,则把当前τ值作为最终τ值并进入步骤STEP5;
STEP5.根据步骤STEP4获得的最终τ值,计算得到最终的α值和λ值,最终的α值,λ值的计算公式分别为
其中xi表示均方根小波系数子带中的第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个的均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ。
所述比较特征列向量与特征列向量之间的相似度,从而判定待识别子图ux所属煤岩类别的过程包括以下步骤:
S1.计算特征列向量与特征列向量之间的相似度,特征列向量与特征列向量之间的相似度计算公式,特征列向量与特征列向量之间的相似度计算公式分别为
其中表示特征列向量与特征列向量之间的相似度,表示特征列向量与特征列向量之间的相似度,特征列向量和中的下标k分别表示煤炭样本序号和岩石样本序号,k的取值为1,2,…,m,N表示两路小波变换的总级数,和中的每一个元素可以分别表示为uxp、ckp和rkp,p的取值为1,2,…,(6×N×3),uxp表示的第p个元素,ckp表示的第p个元素,rkp表示的第p个元素,函数在Γ(z)函数表达式中,z为自变量,t为积分变量,积分区间为[0,+∞),e表示自然常数,函数在ψ(z)函数表达式中,z为自变量,Γ′(z)为Γ(z)的导函数;
S2.进行相似度比较,如果那么判定待识别子图ux属于煤炭,如果那么判定待识别子图ux属于岩石。
附图说明
图1是基于变换域广义伽玛分布的煤岩识别方法的基本流程图;
图2是本发明所述N级两路小波变换的基本流程图;
图3是本发明所述双重列滤波的基本流程图;
图4是本发明所述双重行滤波的基本流程图;
图5是本发明所述计算均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ的基本流程图。
具体实施方式
在对我国山西、河南等地主要煤种和岩种的图像进行实验分析的基础上,本发明提出了基于变换域广义伽玛分布的煤岩识别方法,该方法可以有效判别煤岩与岩石。
下面结合附图和具体实施方式对本发明作进一步的详细描述。
参照图1,基于变换域广义伽玛分布的煤岩识别方法的具体步骤如下:
A.在样本训练阶段,获取相同光照条件下拍摄的m幅已知煤炭样本图像和m幅已知岩石样本图像,截取宽度和高度均为偶数个像素点并且不含非煤岩背景的子图像,煤炭样本子图和岩石样本子图分别记为c1,c2,…,cm和r1,r2,…,rm;
B.分别对样本子图c1,c2,…,cm和r1,r2,…,rm进行N级两路小波变换,每幅样本子图得到12×N个高频部分小波系数子带,每个高频部分小波系数子带记为表示样本子图在经过第i级两路小波变换以后的第d个方向的第j路高频部分小波系数子带,其中i表示两路小波变换的级数序号,i的取值为1,2,…,N,d表示两路小波变换的方向序号,d的取值为1,2,…,6,j表示两路小波变换的路序号,j的取值为1,2;
C.对步骤B所述的每幅样本子图的级数序号相同并且方向序号也相同的两个高频部分小波系数子带进行求均方根操作,得到6×N个均方根小波系数子带,每一个均方根小波系数子带记为RmsCoefi-d,RmsCoefi-d表示样本子图在经过第i级两路小波变换以后的第d个方向的均方根小波系数子带,RmsCoefi-d中每一个元素的计算公式为
式中u和v分别表示RmsCoefi-d、的行号和列号,RmsCoefi-d(u)(v)表示RmsCoefi-d中行号为u、列号为v的元素,表示中行号为u、列号为v的元素,表示中行号为u、列号为v的元素;
D.分别计算步骤C所述的每一个均方根小波系数子带RmsCoefi-d的均值μi-d和方差μi-d的计算公式为其中U和V分别表示RmsCoefi-d的总行数和总列数,的计算公式为 其中U和V分别表示RmsCoefi-d的总行数和总列数;
E.根据步骤D所计算的每一个均方根小波系数子带的均值和方差,对步骤C所述的每一级两路小波变换以后的6个方向的均方根小波系数子带按照均值和方差的乘积值从大到小的顺序排列,每一个经过排序以后的均方根小波系数子带记为RmsCoefSi-n,RmsCoefSi-n表示样本子图在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带,n表示均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排列的排列序号,n的取值为1,2,…,6;
F.分别计算步骤E所述的经过排序以后的每一个均方根小波系数子带在服从概率密度函数为的广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ,其中x表示均方根小波系数子带中的元素,也就是均方根小波系数,函数在Г(α)函数表达式中,α为自变量,t为积分变量,积分区间为[0,+∞),e表示自然常数;
G.利用每一幅煤炭样本子图ck或岩石样本子图rk经过N级两路小波变换、高频部分小波系数子带求均方根操作并且按上述步骤E排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ,构造出一个6×N×3维的特征列向量或其中下标k表示样本序号,k的取值为1,2,…,m,然后把这2×m个特征列向量保存到数据存储设备中,用于后续的煤岩识别阶段;
H.在煤岩识别阶段,获取相同光照条件下拍摄的未知类别煤岩图像,截取宽度和高度均为偶数个像素点并且不含非煤岩背景的待识别子图ux;
I.与上述步骤B类似,对ux进行N级两路小波变换,从而得到12×N个高频部分小波系数子带,每个高频部分小波系数子带记为表示待识别子图ux在经过第i级两路小波变换以后的第d个方向的第j路高频部分小波系数子带;
J.与上述步骤C类似,对步骤I所述的待识别子图ux的级数序号相同并且方向序号也相同的两个高频部分小波系数子带进行求均方根操作,得到6×N个均方根小波系数子带,每一个均方根小波系数子带记为RmsUxCfi-d,RmsUxCfi-d表示待识别子图ux在经过第i级两路小波变换以后的第d个方向的均方根小波系数子带,RmsUxCfi-d中每一个元素的计算公式为
式中p和q分别表示RmsUxCfi-d、的行号和列号,RmsUxCfi-d(p)(q)表示RmsUxCfi-d中行号为p、列号为q的元素,表示中行号为p、列号为q的元素,表示中行号为p、列号为q的元素;
K.与上述步骤D类似,分别计算步骤J所述的每一个均方根小波系数子带RmsUxCfi-d的均值和方差的计算公式为 其中P和Q分别表示均方根小波系数子带RmsUxCfi-d的总行数和总列数,的计算公式为 其中P和Q分别表示均方根小波系数子带RmsUxCfi-d的总行数和总列数;
L.与上述步骤E类似,根据步骤K所计算的每一个均方根小波系数子带的均值和方差,对步骤J所述的每一级两路小波变换以后的6个方向的均方根小波系数子带按照均值和方差的乘积值从大到小的顺序排列,每一个经过排序以后的均方根小波系数子带记为RmsUxCfSi-n,RmsUxCfSi-n表示待识别子图ux在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带;
M.与上述步骤F类似,分别计算步骤L所述的经过排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ;
N.与上述步骤G类似,利用待识别子图ux经过N级两路小波变换、高频部分小波系数子带求均方根操作并且按上述步骤L排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ,构造出一个6×N×3维的特征列向量
O.比较特征列向量与特征列向量之间的相似度,其中下标k表示样本序号,k的取值为1,2,…,m,从而判定待识别子图ux所属的煤岩类别。
参照图2,N级两路小波变换的具体步骤如下:
(1)初始化设置两路小波变换的级数序号i=1;
(2)彩色图像灰度化处理,即把24位煤岩彩色图像转换成8位煤岩灰度图像,采用的转换公式为Y=0.299R+0.587G+0.114B,其中R,G和B分别表示转换前彩色图像的红色,绿色和蓝色分量,Y表示转换后灰度图像的像素值;
(3)用一维的低通滤波器h0对煤岩灰度图像进行逐列滤波,从而得到低频部分小波系数子带CoefLo;
(4)用一维的高通滤波器h1对煤岩灰度图像进行逐列滤波,从而得到高频部分小波系数子带CoefHi;
(5)用步骤(3)所述的低通滤波器h0对步骤(4)所述的高频部分小波系数子带CoefHi进行逐行滤波,从而得到高频部分小波系数子带CoefHL;
(6)用步骤(4)所述的高通滤波器h1对步骤(3)所述的低频部分小波系数子带CoefLo进行逐行滤波,从而得到高频部分小波系数子带CoefLH;
(7)用步骤(4)所述的高通滤波器h1对步骤(4)所述的高频部分小波系数子带CoefHi进行逐行滤波,从而得到高频部分小波系数子带CoefHH;
(8)用步骤(3)所述的低通滤波器h0对步骤(3)所述的低频部分小波系数子带CoefLo进行逐行滤波,从而得到低频部分小波系数子带CoefLL;
(9)利用步骤(5)所述的高频部分小波系数子带CoefHL,分别构造煤岩图像在经过第1级两路小波变换以后的第1个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第1个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第6个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第6个方向的第2路高频部分小波系数子带它们中各元素的计算公式分别为
上述各式中u和v分别表示的行号和列号, 和分别表示和中行号为u、列号为v的元素,CoefHL(2u-1)(2v-1)表示CoefHL中行号为2u-1、列号为2v-1的元素,CoefHL(2u)(2v)表示CoefHL中行号为2u、列号为2v的元素,CoefHL(2u-1)(2v)表示CoefHL中行号为2u-1、列号为2v的元素,CoefHL(2u)(2v-1)表示CoefHL中行号为2u、列号为2v-1的元素,u和v的最小取值为1,u的最大取值为CoefHL总行数的一半,v的最大取值为CoefHL总列数的一半;
(10)利用步骤(6)所述的高频部分小波系数子带CoefLH,分别构造煤岩图像在经过第1级两路小波变换以后的第3个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第3个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第4个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第4个方向的第2路高频部分小波系数子带它们中各元素的计算公式分别为
上述各式中u和v分别表示的行号和列号, 和分别表示和中行号为u、列号为v的元素,CoefLH(2u-1)(2v-1)表示CoefLH中行号为2u-1、列号为2v-1的元素,CoefLH(2u)(2v)表示CoefLH中行号为2u、列号为2v的元素,CoefLH(2u-1)(2v)表示CoefLH中行号为2u-1、列号为2v的元素,CoefLH(2u)(2v-1)表示CoefLH中行号为2u、列号为2v-1的元素,u和v的最小取值为1,u的最大取值为CoefLH总行数的一半,v的最大取值为CoefLH总列数的一半;
(11)利用步骤(7)所述的高频部分小波系数子带CoefHH,分别构造煤岩图像在经过第1级两路小波变换以后的第2个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第2个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第5个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第5个方向的第2路高频部分小波系数子带它们中各元素的计算公式分别为
上述各式中u和v分别表示的行号和列号, 和分别表示和中行号为u、列号为v的元素,CoefHH(2u-1)(2v-1)表示CoefHH中行号为2u-1、列号为2v-1的元素,CoefHH(2u)(2v)表示CoefHH中行号为2u、列号为2v的元素,CoefHH(2u-1)(2v)表示CoefHH中行号为2u-1、列号为2v的元素,CoefHH(2u)(2v-1)表示CoefHH中行号为2u、列号为2v-1的元素,u和v的最小取值为1,u的最大取值为CoefHH总行数的一半,v的最大取值为CoefHH总列数的一半;
(12)两路小波变换的级数序号自增1,即i=i+1;
(13)动态申请与CoefLL数据结构完全一致的临时数据存储区TempCoef,并将其数据内容初始化为CoefLL的数据内容;
(14)判断i≤N条件是否成立,如果是,则进入以下步骤(15)-(28)的迭代循环,如果否,则转到步骤(29);
(15)如果TempCoef数据的总行数不是4的倍数,那么对TempCoef数据进行修正,在TempCoef数据的第一行Rowfirst前面追加一行Row′first并且Row′first的数据内容用Rowfirst的数据内容进行填充,在TempCoef数据的最后一行Rowlast后面追加一行Row′last并且Row′last的数据内容用Rowlast的数据内容进行填充;
(16)如果TempCoef数据的总列数不是4的倍数,那么对TempCoef数据进行修正,在TempCoef数据的第一列Colfirst前面追加一列Col′first并且Col′first的数据内容用Colfirst的数据内容进行填充,在TempCoef数据的最后一列Collast后面追加一列Col′last并且Col′last的数据内容用Collast的数据内容进行填充;
(17)用两组一维的低通滤波器h00、h01对TempCoef数据进行双重列滤波,从而得到低频部分小波系数子带CoefL2-i,其中下标i表示两路小波变换的级数序号;
(18)用两组一维的高通滤波器h10、h11对TempCoef数据进行双重列滤波,从而得到高频部分小波系数子带CoefH2-i,其中下标i表示两路小波变换的级数序号;
(19)用步骤(17)所述的低通滤波器h00、h01对步骤(18)所述的高频部分小波系数子带CoefH2-i进行双重行滤波,从而得到高频部分小波系数子带CoefHL2-i,其中下标i表示两路小波变换的级数序号;
(20)用步骤(18)所述的高通滤波器h10、h11对步骤(17)所述的低频部分小波系数子带CoefL2-i进行双重行滤波,从而得到高频部分小波系数子带CoefLH2-i,其中下标i表示两路小波变换的级数序号;
(21)用步骤(18)所述的高通滤波器h10、h11对步骤(18)所述的高频部分小波系数子带CoefH2-i进行双重行滤波,从而得到高频部分小波系数子带CoefHH2-i,其中下标i表示两路小波变换的级数序号;
(22)用步骤(17)所述的低通滤波器h00、h01对步骤(17)所述的低频部分小波系数子带CoefL2-i进行双重行滤波,从而得到低频部分小波系数子带CoefLL2-i,其中下标i表示两路小波变换的级数序号;
(23)利用步骤(19)所述的高频部分小波系数子带CoefHL2-i,分别构造煤岩图像在经过第i级两路小波变换以后的第1个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第1个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第6个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第6个方向的第2路高频部分小波系数子带它们中各元素的计算公式分别为
上述各式中u和v分别表示的行号和列号, 和分别表示和中行号为u、列号为v的元素,CoefHL2-i(2u-1)(2v-1)表示CoefHL2-i中行号为2u-1、列号为2v-1的元素,CoefHL2-i(2u)(2v)表示CoefHL2-i中行号为2u、列号为2v的元素,CoefHL2-i(2u-1)(2v)表示CoefHL2-i中行号为2u-1、列号为2v的元素,CoefHL2-i(2u)(2v-1)表示CoefHL2-i中行号为2u、列号为2v-1的元素,u和v的最小取值为1,u的最大取值为CoefHL2-i总行数的一半,v的最大取值为CoefHL2-i总列数的一半;
(24)利用步骤(20)所述的高频部分小波系数子带CoefLH2-i,分别构造煤岩图像在经过第i级两路小波变换以后的第3个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第3个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第4个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第4个方向的第2路高频部分小波系数子带它们中各元素的计算公式分别为
上述各式中u和v分别表示的行号和列号, 和分别表示和中行号为u、列号为v的元素,CoefLH2-i(2u-1)(2v-1)表示CoefLH2-i中行号为2u-1、列号为2v-1的元素,CoefLH2-i(2u)(2v)表示CoefLH2-i中行号为2u、列号为2v的元素,CoefLH2-i(2u-1)(2v)表示CoefLH2-i中行号为2u-1、列号为2v的元素,CoefLH2-i(2u)(2v-1)表示CoefLH2-i中行号为2u、列号为2v-1的元素,u和v的最小取值为1,u的最大取值为CoefLH2-i总行数的一半,v的最大取值为CoefLH2-i总列数的一半;
(25)利用步骤(21)所述的高频部分小波系数子带CoefHH2-i,分别构造煤岩图像在经过第i级两路小波变换以后的第2个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第2个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第5个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第5个方向的第2路高频部分小波系数子带它们中各元素的计算公式分别为
上述各式中u和v分别表示的行号和列号, 和分别表示和中行号为u、列号为v的元素,CoefHH2-i(2u-1)(2v-1)表示CoefHH2-i中行号为2u-1、列号为2v-1的元素,CoefHH2-i(2u)(2v)表示CoefHH2-i中行号为2u、列号为2v的元素,CoefHH2-i(2u-1)(2v)表示CoefHH2-i中行号为2u-1、列号为2v的元素,CoefHH2-i(2u)(2v-1)表示CoefHH2-i中行号为2u、列号为2v-1的元素,u和v的最小取值为1,u的最大取值为CoefHH2-i总行数的一半,v的最大取值为CoefHH2-i总列数的一半;
(26)释放临时数据存储区TempCoef所占用的内存空间,重新动态申请与步骤(22)所述CoefLL2-i的数据结构完全一致的临时数据存储区TempCoef,并将其数据内容初始化为CoefLL2-i的数据内容;
(27)两路小波变换的级数序号自增1,即i=i+1;
(28)转到步骤(14),进行下一级两路小波变换;
(29)释放临时数据存储区TempCoef所占用的内存空间,完成N级两路小波变换。
所述一维的低通滤波器h0和高通滤波器h1的序列长度均为奇数,并且h0序列和h1序列均为中心对称序列。
所述一维的低通滤波器h00、h01和高通滤波器h10、h11的序列长度相等并且均为偶数,h01序列为h00序列的反转序列,h10序列通过对h00序列进行反转操作然后对偶数序号采样点的序列元素进行符号取反操作而得到,h11序列为h10序列的反转序列。
如图3所示,在给定两组序列长度相等的一维滤波器hf、hs和滤波对象X的前提下,所述双重列滤波包括以下步骤:
STC1.对X进行行对称延拓得到Xex,如果X的总行数和总列数分别为RowsX和ColsX,hf和hs的序列长度为len,那么Xex的总行数和总列数分别为(RowsX+2×len)和ColsX,并且Xex的第1,2,…,len行的数据内容分别用X的第len,len-1,…,1行的数据内容进行填充,Xex的第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.抽取hf中序号为奇数的滤波系数构成子滤波器hfo,抽取hf中序号为偶数的滤波系数构成子滤波器hfe,抽取hs中序号为奇数的滤波系数构成子滤波器hso,抽取hs中序号为偶数的滤波系数构成子滤波器hse;
STC3.对步骤STC1所述的Xex进行修正,即舍弃Xex的起始两行和末尾两行,从而得到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所述的子滤波器hfo对步骤STC5所述的X′exB进行逐列滤波,从而得到ResB;
STC7.用步骤STC2所述的子滤波器hfe对步骤STC5所述的X′exD进行逐列滤波,从而得到ResD;
STC8.用步骤STC2所述的子滤波器hso对步骤STC5所述的X′exA进行逐列滤波,从而得到ResA;
STC9.用步骤STC2所述的子滤波器hse对步骤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.把一维滤波器hf和hs视为两个向量,如果这两个向量的内积大于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所示,在给定两组序列长度相等的一维滤波器hf、hs和滤波对象X的前提下,所述双重行滤波包括以下步骤:
STR1.对X进行列对称延拓得到Xex,如果X的总行数和总列数分别为RowsX和ColsX,hf和hs的序列长度为len,那么Xex的总行数和总列数分别为RowsX和(ColsX+2×len),并且Xex的第1,2,…,len列的数据内容分别用X的第len,len-1,…,1列的数据内容进行填充,Xex的第len+1,len+2,…,len+ColsX列的数据内容分别用X的第1,2,…,ColsX列的数据内容进行填充,Xex的第len+ColsX+1,len+ColsX+2,…,len+ColsX+len列的数据内容分别用X的第ColsX,ColsX-1,…,ColsX-(len-1)列的数据内容进行填充;
STR2.抽取hf中序号为奇数的滤波系数构成子滤波器hfo,抽取hf中序号为偶数的滤波系数构成子滤波器hfe,抽取hs中序号为奇数的滤波系数构成子滤波器hso,抽取hs中序号为偶数的滤波系数构成子滤波器hse;
STR3.对步骤STR1所述的Xex进行修正,即舍弃Xex的起始两列和末尾两列,从而得到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所述的子滤波器hfo对步骤STR5所述的X′exB进行逐行滤波,从而得到ResB;
STR7.用步骤STR2所述的子滤波器hfe对步骤STR5所述的X′exD进行逐行滤波,从而得到ResD;
STR8.用步骤STR2所述的子滤波器hso对步骤STR5所述的X′exA进行逐行滤波,从而得到ResA;
STR9.用步骤STR2所述的子滤波器hse对步骤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.把一维滤波器hf和hs视为两个向量,如果这两个向量的内积大于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=0.05作为τ参数的初始迭代值,进入以下步骤STEP2-STEP4的循环迭代过程;
STEP2.把τ作为自变量代入并计算
其中xi表示均方根小波系数子带中的第i个均方根小波系数,Num表示均方根小波系数子带中的均方根小波系数总数,ln表示以自然常数e为底数的对数函数,即自然对数;
STEP3.把τ作为自变量代入并计算
其中xi表示均方根小波系数子带中的第i个均方根小波系数,Num表示均方根小波系数子带中的均方根小波系数总数,ln表示以自然常数e为底数的对数函数,即自然对数;
STEP4.计算若|Δτ|大于计算精度ε,则通过公式τ=τ-|Δτ|来更新τ,然后执行步骤STEP2-STEP4的循环迭代过程,若|Δτ|小于或等于计算精度ε,则把当前τ值作为最终τ值并进入步骤STEP5;
STEP5.根据步骤STEP4获得的最终τ值,计算得到最终的α值和λ值,最终的α值,λ值的计算公式分别为
其中xi表示均方根小波系数子带中的第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个的均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ。
所述比较特征列向量与特征列向量之间的相似度,从而判定待识别子图ux所属煤岩类别的过程包括以下步骤:
S1.计算特征列向量与特征列向量之间的相似度,特征列向量与特征列向量之间的相似度计算公式,特征列向量与特征列向量之间的相似度计算公式分别为
其中表示特征列向量与特征列向量之间的相似度,表示特征列向量与特征列向量之间的相似度,特征列向量和中的下标k分别表示煤炭样本序号和岩石样本序号,k的取值为1,2,…,m,N表示两路小波变换的总级数,和中的每一个元素可以分别表示为uxp、ckp和rkp,p的取值为1,2,…,(6×N×3),uxp表示的第p个元素,ckp表示的第p个元素,rkp表示的第p个元素,函数在Γ(z)函数表达式中,z为自变量,t为积分变量,积分区间为[0,+∞),e表示自然常数,函数在ψ(z)函数表达式中,z为自变量,Γ′(z)为Γ(z)的导函数;
S2.进行相似度比较,如果那么判定待识别子图ux属于煤炭,如果那么判定待识别子图ux属于岩石。
需要指出的是,以上所述实施实例用于进一步说明本发明,实施实例不应被视为限制本发明的范围。
Claims (9)
1.基于变换域广义伽玛分布的煤岩识别方法,其特征在于,包括以下步骤:
A.在样本训练阶段,获取相同光照条件下拍摄的m幅已知煤炭样本图像和m幅已知岩石样本图像,截取宽度和高度均为偶数个像素点并且不含非煤岩背景的子图像,煤炭样本子图和岩石样本子图分别记为c1,c2,…,cm和r1,r2,…,rm;
B.分别对样本子图c1,c2,…,cm和r1,r2,…,rm进行N级两路小波变换,每幅样本子图得到12×N个高频部分小波系数子带
表示样本子图在经过第i级两路小波变换以后的第d个方向的第j路高频部分小波系数子带,其中i表示两路小波变换的级数序号,i的取值为1,2,…,N,d表示两路小波变换的方向序号,d的取值为1,2,…,6,j表示两路小波变换的路序号,j的取值为1,2;
C.对步骤B所述的每幅样本子图的级数序号相同并且方向序号也相同的两个高频部分小波系数子带进行求均方根操作,得到6×N个均方根小波系数子带RmsCoefi-d,RmsCoefi-d表示样本子图在经过第i级两路小波变换以后的第d个方向的均方根小波系数子带;
D.分别计算步骤C所述的每一个均方根小波系数子带RmsCoefi-d的均值μi-d和方差
E.根据步骤D所计算的每一个均方根小波系数子带的均值和方差,对步骤C所述的每一级两路小波变换以后的6个方向的均方根小波系数子带按照均值和方差的乘积值从大到小的顺序排列,每一个经过排序以后的均方根小波系数子带记为RmsCoefSi-n;
RmsCoefSi-n表示样本子图在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带,n表示均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排列的排列序号,n的取值为1,2,…,6;
F.分别计算步骤E所述的经过排序以后的每一个均方根小波系数子带在服从概率密度函数为的广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ,其中x表示均方根小波系数子带中的元素,也就是均方根小波系数,函数 在Г(α)函数表达式中,α为自变量,t为积分变量,积分区间为[0,+∞),e表示自然常数;
G.利用每一幅煤炭样本子图ck或岩石样本子图rk经过N级两路小波变换、高频部分小波 系数子带求均方根操作并且按上述步骤E排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ,构造出一个6×N×3维的特征列向量或其中下标k表示样本序号,k的取值为1,2,…,m,然后把这2×m个特征列向量保存到数据存储设备中,用于后续的煤岩识别阶段;
H.在煤岩识别阶段,获取相同光照条件下拍摄的未知类别煤岩图像,截取宽度和高度均为偶数个像素点并且不含非煤岩背景的待识别子图ux;
I.与上述步骤B类似,对ux进行N级两路小波变换,从而得到12×N个高频部分小波系数子带,每个高频部分小波系数子带记为表示待识别子图ux在经过第i级两路小波变换以后的第d个方向的第j路高频部分小波系数子带;
J.与上述步骤C类似,对步骤I所述的待识别子图ux的级数序号相同并且方向序号也相同的两个高频部分小波系数子带进行求均方根操作,得到6×N个均方根小波系数子带,每一个均方根小波系数子带记为RmsUxCfi-d,RmsUxCfi-d表示待识别子图ux在经过第i级两路小波变换以后的第d个方向的均方根小波系数子带;
K.与上述步骤D类似,分别计算步骤J所述的每一个均方根小波系数子带RmsUxCfi-d的均值和方差
L.与上述步骤E类似,根据步骤K所计算的每一个均方根小波系数子带的均值和方差,对步骤J所述的每一级两路小波变换以后的6个方向的均方根小波系数子带按照均值和方差的乘积值从大到小的顺序排列,每一个经过排序以后的均方根小波系数子带记为RmsUxCfSi-n,RmsUxCfSi-n表示待识别子图ux在经过第i级两路小波变换以后的6个均方根小波系数子带中均方根小波系数子带均值和均方根小波系数子带方差的乘积值由大到小顺序排第n个的均方根小波系数子带;
M.与上述步骤F类似,分别计算步骤L所述的经过排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ;
N.与上述步骤G类似,利用待识别子图ux经过N级两路小波变换、高频部分小波系数子带求均方根操作并且按上述步骤L排序以后的每一个均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ,构造出一个6×N×3维的特征列向量
O.比较特征列向量与特征列向量之间的相似度,如果 那么判定待识别子图ux属于煤炭,如果 那么判定待识别子图ux属于岩石;
其中,表示特征列向量与特征列向量之间的相似度,表示特征列向量与特征列向量之间的相似度,其中下标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)用一维的低通滤波器h0对煤岩灰度图像进行逐列滤波,从而得到低频部分小波系数子带CoefLo;
(4)用一维的高通滤波器h1对煤岩灰度图像进行逐列滤波,从而得到高频部分小波系数子带CoefHi;
(5)用步骤(3)所述的低通滤波器h0对步骤(4)所述的高频部分小波系数子带CoefHi进行逐行滤波,从而得到高频部分小波系数子带CoefHL;
(6)用步骤(4)所述的高通滤波器h1对步骤(3)所述的低频部分小波系数子带CoefLo进行逐行滤波,从而得到高频部分小波系数子带CoefLH;
(7)用步骤(4)所述的高通滤波器h1对步骤(4)所述的高频部分小波系数子带CoefHi进行逐行滤波,从而得到高频部分小波系数子带CoefHH;
(8)用步骤(3)所述的低通滤波器h0对步骤(3)所述的低频部分小波系数子带CoefLo进行逐行滤波,从而得到低频部分小波系数子带CoefLL;
(9)利用步骤(5)所述的高频部分小波系数子带CoefHL,分别构造煤岩图像在经过第1级两路小波变换以后的第1个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第1个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第6个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第6个方向的第2路高频部分小波系数子带
(10)利用步骤(6)所述的高频部分小波系数子带CoefLH,分别构造煤岩图像在经过第 1级两路小波变换以后的第3个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第3个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第4个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第4个方向的第2路高频部分小波系数子带
(11)利用步骤(7)所述的高频部分小波系数子带CoefHH,分别构造煤岩图像在经过第1级两路小波变换以后的第2个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第2个方向的第2路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第5个方向的第1路高频部分小波系数子带煤岩图像在经过第1级两路小波变换以后的第5个方向的第2路高频部分小波系数子带
(12)两路小波变换的级数序号自增1,即i=i+1;
(13)动态申请与CoefLL数据结构完全一致的临时数据存储区TempCoef,并将其数据内容初始化为CoefLL的数据内容;
(14)判断i≤N条件是否成立,如果是,则进入以下步骤(15)-(28)的迭代循环,如果否,则转到步骤(29);
(15)如果TempCoef数据的总行数不是4的倍数,那么对TempCoef数据进行修正,在TempCoef数据的第一行Rowfirst前面追加一行Row′first并且Row′first的数据内容用Rowfirst的数据内容进行填充,在TempCoef数据的最后一行Rowlast后面追加一行Row′last并且Row′last的数据内容用Rowlast的数据内容进行填充;
(16)如果TempCoef数据的总列数不是4的倍数,那么对TempCoef数据进行修正,在TempCoef数据的第一列Colfirst前面追加一列Col′first并且Col′first的数据内容用Colfirst的数据内容进行填充,在TempCoef数据的最后一列Collast后面追加一列Col′last并且Col′last的数据内容用Collast的数据内容进行填充;
(17)用两组一维的低通滤波器h00、h01对TempCoef数据进行双重列滤波,从而得到低频部分小波系数子带CoefL2-i,其中下标i表示两路小波变换的级数序号;
(18)用两组一维的高通滤波器h10、h11对TempCoef数据进行双重列滤波,从而得到高频部分小波系数子带CoefH2-i,其中下标i表示两路小波变换的级数序号;
(19)用步骤(17)所述的低通滤波器h00、h01对步骤(18)所述的高频部分小波系数子带CoefH2-i进行双重行滤波,从而得到高频部分小波系数子带CoefHL2-i,其中下标i表示两 路小波变换的级数序号;
(20)用步骤(18)所述的高通滤波器h10、h11对步骤(17)所述的低频部分小波系数子带CoefL2-i进行双重行滤波,从而得到高频部分小波系数子带CoefLH2-i,其中下标i表示两路小波变换的级数序号;
(21)用步骤(18)所述的高通滤波器h10、h11对步骤(18)所述的高频部分小波系数子带CoefH2-i进行双重行滤波,从而得到高频部分小波系数子带CoefHH2-i,其中下标i表示两路小波变换的级数序号;
(22)用步骤(17)所述的低通滤波器h00、h01对步骤(17)所述的低频部分小波系数子带CoefL2-i进行双重行滤波,从而得到低频部分小波系数子带CoefLL2-i,其中下标i表示两路小波变换的级数序号;
(23)利用步骤(19)所述的高频部分小波系数子带CoefHL2-i,分别构造煤岩图像在经过第i级两路小波变换以后的第1个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第1个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第6个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第6个方向的第2路高频部分小波系数子带
(24)利用步骤(20)所述的高频部分小波系数子带CoefLH2-i,分别构造煤岩图像在经过第i级两路小波变换以后的第3个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第3个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第4个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第4个方向的第2路高频部分小波系数子带
(25)利用步骤(21)所述的高频部分小波系数子带CoefHH2-i,分别构造煤岩图像在经过第i级两路小波变换以后的第2个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第2个方向的第2路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第5个方向的第1路高频部分小波系数子带煤岩图像在经过第i级两路小波变换以后的第5个方向的第2路高频部分小波系数子带
(26)释放临时数据存储区TempCoef所占用的内存空间,重新动态申请与步骤(22)所述CoefLL2-i的数据结构完全一致的临时数据存储区TempCoef,并将其数据内容初始化为CoefLL2-i的数据内容;
(27)两路小波变换的级数序号自增1,即i=i+1;
(28)转到步骤(14),进行下一级两路小波变换;
(29)释放临时数据存储区TempCoef所占用的内存空间,完成N级两路小波变换。
3.根据权利要求1所述的基于变换域广义伽玛分布的煤岩识别方法,其特征在于,所述计算均方根小波系数子带在服从广义伽玛分布条件下的指数形状参数α、形状参数τ和尺度参数λ包括以下步骤:
STEP1.把τ0=0.05作为τ参数的初始迭代值,进入以下步骤STEP2-STEP4的循环迭代过程;
STEP2.把τ作为自变量代入并计算
其中xi表示均方根小波系数子带中的第i个均方根小波系数,Num表示均方根小波系数子带中的均方根小波系数总数,ln表示以自然常数e为底数的对数函数,即自然对数;
STEP3.把τ作为自变量代入并计算
其中xi表示均方根小波系数子带中的第i个均方根小波系数,Num表示均方根小波系数子带中的均方根小波系数总数,ln表示以自然常数e为底数的对数函数,即自然对数;
STEP4.计算若|Δτ|大于计算精度ε,则通过公式τ=τ-|Δτ|来更新τ,然后执行步骤STEP2-STEP4的循环迭代过程,若|Δτ|小于或等于计算精度ε,则把当前τ值作为最终τ值并进入步骤STEP5;
STEP5.根据步骤STEP4获得的最终τ值,计算得到最终的α值和λ值,最终的α值,λ值的计算公式分别为
其中xi表示均方根小波系数子带中的第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表示两路小波变换的总级数,和中的每一个元素可以分别表示为uxp、ckp和rkp,p的取值为1,2,…,(6×N×3),uxp表示的第p个元素,ckp表示的第p个元素,rkp表示的第p个元素,函数在Г(z)函数表达式中,z为自变量,t为积分变量,积分区间为[0,+∞),e表示自然常数,函数在ψ(z)函数表达式中,z为自变量,Г′(z)为Г(z)的导函数。
6.根据权利要求2所述的基于变换域广义伽玛分布的煤岩识别方法,其特征在于,所述一维的低通滤波器h0和高通滤波器h1的序列长度均为奇数,并且h0序列和h1序列均为中心对称序列。
7.根据权利要求2所述的基于变换域广义伽玛分布的煤岩识别方法,其特征在于,所述一维的低通滤波器h00、h01和高通滤波器h10、h11的序列长度相等并且均为偶数,h01序列为h00序列的反转序列,h10序列通过对h00序列进行反转操作然后对偶数序号采样点的序列元素进行符号取反操作而得到,h11序列为h10序列的反转序列。
8.根据权利要求2所述的基于变换域广义伽玛分布的煤岩识别方法,其特征在于,在给定两组序列长度相等的一维滤波器hf、hs和滤波对象X的前提下,所述双重列滤波包括以下步骤:
STC1.对X进行行对称延拓得到Xex,如果X的总行数和总列数分别为RowsX和ColsX,hf和hs的序列长度为len,那么Xex的总行数和总列数分别为(RowsX+2×len)和ColsX,并且Xex的第1,2,…,len行的数据内容分别用X的第len,len-1,…,1行的数据内容进行填充,Xex的第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.抽取hf中序号为奇数的滤波系数构成子滤波器hfo,抽取hf中序号为偶数的滤波系数构成子滤波器hfe,抽取hs中序号为奇数的滤波系数构成子滤波器hso,抽取hs中序号为 偶数的滤波系数构成子滤波器hse;
STC3.对步骤STC1所述的Xex进行修正,即舍弃Xex的起始两行和末尾两行,从而得到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所述的子滤波器hfo对步骤STC5所述的X′exB进行逐列滤波,从而得到ResB;
STC7.用步骤STC2所述的子滤波器hfe对步骤STC5所述的X′exD进行逐列滤波,从而得到ResD;
STC8.用步骤STC2所述的子滤波器hso对步骤STC5所述的X′exA进行逐列滤波,从而得到ResA;
STC9.用步骤STC2所述的子滤波器hse对步骤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.把一维滤波器hf和hs视为两个向量,如果这两个向量的内积大于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所述的基于变换域广义伽玛分布的煤岩识别方法,其特征在于,在给定两组序列长度相等的一维滤波器hf、hs和滤波对象X的前提下,所述双重行滤波包括以下步骤:
STR1.对X进行列对称延拓得到Xex,如果X的总行数和总列数分别为RowsX和ColsX,hf和hs的序列长度为len,那么Xex的总行数和总列数分别为RowsX和(ColsX+2×len),并且Xex的第1,2,…,len列的数据内容分别用X的第len,len-1,…,l列的数据内容进行填充,Xex的第len+1,len+2,…,len+ColsX列的数据内容分别用X的第1,2,…,ColsX列的数据内容进行填充,Xex的第len+ColsX+1,len+ColsX+2,…,len+ColsX+len列的数据内容分别用X的第ColsX,ColsX-1,…,ColsX-(len-1)列的数据内容进行填充;
STR2.抽取hf中序号为奇数的滤波系数构成子滤波器hfo,抽取hf中序号为偶数的滤波系数构成子滤波器hfe,抽取hs中序号为奇数的滤波系数构成子滤波器hso,抽取hs中序号为偶数的滤波系数构成子滤波器hse;
STR3.对步骤STR1所述的Xex进行修正,即舍弃Xex的起始两列和末尾两列,从而得到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所述的子滤波器hfo对步骤STR5所述的X′exB进行逐行滤波,从而得到ResB;
STR7.用步骤STR2所述的子滤波器hfe对步骤STR5所述的X′exD进行逐行滤波,从而得 到ResD;
STR8.用步骤STR2所述的子滤波器hso对步骤STR5所述的X′exA进行逐行滤波,从而得到ResA;
STR9.用步骤STR2所述的子滤波器hse对步骤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.把一维滤波器hf和hs视为两个向量,如果这两个向量的内积大于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。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510418701.1A CN104966074B (zh) | 2015-07-17 | 2015-07-17 | 基于变换域广义伽玛分布的煤岩识别方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510418701.1A CN104966074B (zh) | 2015-07-17 | 2015-07-17 | 基于变换域广义伽玛分布的煤岩识别方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104966074A true CN104966074A (zh) | 2015-10-07 |
CN104966074B CN104966074B (zh) | 2018-06-19 |
Family
ID=54220110
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510418701.1A Active CN104966074B (zh) | 2015-07-17 | 2015-07-17 | 基于变换域广义伽玛分布的煤岩识别方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104966074B (zh) |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106845560A (zh) * | 2017-03-01 | 2017-06-13 | 中国矿业大学(北京) | 基于可控塔式分解和字典学习的煤岩识别方法 |
CN107169524A (zh) * | 2017-05-31 | 2017-09-15 | 中国矿业大学(北京) | 基于完备局部二值模式重构残差的煤岩识别方法 |
CN111779524A (zh) * | 2020-06-30 | 2020-10-16 | 中国矿业大学 | 一种综放工作面液压支架智能放煤方法 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20030086593A1 (en) * | 2001-05-31 | 2003-05-08 | Chengjun Liu | Feature based classification |
CN102880858A (zh) * | 2012-08-30 | 2013-01-16 | 中国矿业大学(北京) | 一种煤岩图像自动识别方法 |
CN102930253A (zh) * | 2012-10-31 | 2013-02-13 | 中国矿业大学(北京) | 一种基于图像离散多小波变换的煤岩识别方法 |
CN103207999A (zh) * | 2012-11-07 | 2013-07-17 | 中国矿业大学(北京) | 一种基于煤岩图像特征抽取以及分类识别的煤岩分界方法和系统 |
CN103530621A (zh) * | 2013-11-04 | 2014-01-22 | 中国矿业大学(北京) | 一种基于bp神经网络的煤岩图像识别方法 |
CN104732239A (zh) * | 2015-04-08 | 2015-06-24 | 中国矿业大学(北京) | 基于小波域非对称广义高斯模型的煤岩分类方法 |
-
2015
- 2015-07-17 CN CN201510418701.1A patent/CN104966074B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20030086593A1 (en) * | 2001-05-31 | 2003-05-08 | Chengjun Liu | Feature based classification |
CN102880858A (zh) * | 2012-08-30 | 2013-01-16 | 中国矿业大学(北京) | 一种煤岩图像自动识别方法 |
CN102930253A (zh) * | 2012-10-31 | 2013-02-13 | 中国矿业大学(北京) | 一种基于图像离散多小波变换的煤岩识别方法 |
CN103207999A (zh) * | 2012-11-07 | 2013-07-17 | 中国矿业大学(北京) | 一种基于煤岩图像特征抽取以及分类识别的煤岩分界方法和系统 |
CN103530621A (zh) * | 2013-11-04 | 2014-01-22 | 中国矿业大学(北京) | 一种基于bp神经网络的煤岩图像识别方法 |
CN104732239A (zh) * | 2015-04-08 | 2015-06-24 | 中国矿业大学(北京) | 基于小波域非对称广义高斯模型的煤岩分类方法 |
Non-Patent Citations (3)
Title |
---|
孙继平等: "基于小波的煤岩图像特征抽取与识别", 《煤炭学报》 * |
孙继平等: "基于支持向量机的煤岩图像特征抽取与分类识别", 《煤炭学报》 * |
秦先祥等: "基于SISE方程的广义gamma分布参数估计方法", 《电子与信息学报》 * |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106845560A (zh) * | 2017-03-01 | 2017-06-13 | 中国矿业大学(北京) | 基于可控塔式分解和字典学习的煤岩识别方法 |
CN106845560B (zh) * | 2017-03-01 | 2020-04-28 | 中国矿业大学(北京) | 基于可控塔式分解和字典学习的煤岩识别方法 |
CN107169524A (zh) * | 2017-05-31 | 2017-09-15 | 中国矿业大学(北京) | 基于完备局部二值模式重构残差的煤岩识别方法 |
CN111779524A (zh) * | 2020-06-30 | 2020-10-16 | 中国矿业大学 | 一种综放工作面液压支架智能放煤方法 |
Also Published As
Publication number | Publication date |
---|---|
CN104966074B (zh) | 2018-06-19 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN110321963B (zh) | 基于融合多尺度多维空谱特征的高光谱图像分类方法 | |
CN101299237B (zh) | 一种基于信息量维数序列的高光谱数据监督分类方法 | |
CN110929607A (zh) | 一种城市建筑物施工进度的遥感识别方法和系统 | |
CN110232419A (zh) | 一种边坡岩石类别自动识别的方法 | |
CN103927514A (zh) | 一种基于随机局部图像特征的煤岩识别方法 | |
CN102297822A (zh) | 一种通过图像分析预测煤粒灰分的方法 | |
CN104850854A (zh) | 一种滑石矿品分选处理方法及滑石矿品分选系统 | |
CN104966074A (zh) | 基于变换域广义伽玛分布的煤岩识别方法 | |
CN113685188A (zh) | 一种基于岩渣物理特征的tbm掘进优化方法 | |
CN106557758A (zh) | 一种沙粒显微图像的多目标自动鉴别方法 | |
CN103927528A (zh) | 基于紧邻域像素灰度联合分布特征的煤岩识别方法 | |
CN104134069A (zh) | 一种页岩显微薄片自动鉴别方法 | |
CN104835142B (zh) | 一种基于纹理特征的车辆排队长度检测方法 | |
CN114898109A (zh) | 一种基于深度学习的斑岩浅成低温热液型矿产预测方法及系统 | |
Wang et al. | Fusion of hyperspectral and lidar data based on dual-branch convolutional neural network | |
CN115272826A (zh) | 一种基于卷积神经网络的图像识别方法、装置及系统 | |
CN112488043B (zh) | 一种基于边缘智能的无人机车辆目标检测方法 | |
CN105243401A (zh) | 基于煤岩结构元学习的煤岩识别方法 | |
CN103942576A (zh) | 一种用空域多尺度随机特征识别煤岩的方法 | |
CN112818777B (zh) | 一种基于密集连接与特征增强的遥感图像目标检测方法 | |
CN104732239B (zh) | 基于小波域非对称广义高斯模型的煤岩分类方法 | |
CN105350963B (zh) | 一种基于相关性度量学习的煤岩识别方法 | |
CN113158770A (zh) | 一种改进的全卷积孪生神经网络的矿区变化检测方法 | |
Ghasemloo et al. | Road and tunnel extraction from SPOT satellite images using neural networks | |
CN115497006B (zh) | 基于动态混合策略的城区遥感图像变化深度监测方法及系统 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |