CN110135515B - 一种基于图像纹理的岩体结构均质区自动分区方法 - Google Patents
一种基于图像纹理的岩体结构均质区自动分区方法 Download PDFInfo
- Publication number
- CN110135515B CN110135515B CN201910437094.1A CN201910437094A CN110135515B CN 110135515 B CN110135515 B CN 110135515B CN 201910437094 A CN201910437094 A CN 201910437094A CN 110135515 B CN110135515 B CN 110135515B
- Authority
- CN
- China
- Prior art keywords
- texture
- image
- rock mass
- mass structure
- window
- 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.)
- Active
Links
Images
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/2135—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on approximation criteria, e.g. principal component analysis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/23—Clustering techniques
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/40—Analysis of texture
- G06T7/41—Analysis of texture based on statistical description of texture
- G06T7/45—Analysis of texture based on statistical description of texture using co-occurrence matrix computation
Landscapes
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- General Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Artificial Intelligence (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- General Engineering & Computer Science (AREA)
- Probability & Statistics with Applications (AREA)
- Computing Systems (AREA)
- Mathematical Physics (AREA)
- Image Analysis (AREA)
Abstract
本发明是一种基于图像纹理的岩体结构均质区自动分区方法。利用岩体结构面在露头区赋存几何形态的纹理特征,在摄影测量方法获取具有实际尺寸的岩体结构面赋存影像的基础上,采用灰度共生矩阵的10个纹理参数来描述岩体结构,通过点对距离、点对方向角与像素尺度的敏感性检测与主成分分析方法筛选出主控纹理参数,以主控纹理参数为聚类指标,提出基于ISODATA聚类方法的岩体结构均质区自动分区方法。包括影像数据采集、灰度图像转换、纹理信息增强、纹理特征提取、纹理参数确定、聚类窗口确定、纹理图像分块、子域纹理参数、矢量数据转换、ISODATA聚类10步骤。规避了现场地质测绘与统计工作量,数据处理过程可实现自动化,可应用于工程全域岩体结构均质区分区。
Description
技术领域
本发明涉及工程现场岩体结构均质区的自动分区方法。属于岩体工程的现场地质调查技术领域。适用于岩体工程全域范围内岩体地质分区工作。
背景技术
1978年国际岩石力学学会实验室和野外试验标准化专门委员会提出了“对岩体中结构面定量描述的推荐方法”,其中规定了对结构面的10个描述指标,包括:产状(Orientation)、间距或密集程度(Spacing)、延续性(Persistence)、粗糙程度(Roughness)、隙宽(Aperture)、岩壁抗压强度、充填情况、渗流、组数和块度等[2]。自上世纪70年代,Hudson等人用概率统计方法对结构面的几何特征进行定量数学描述以来,许多学者对结构面几何特征参量的数学描述研究开展了系统研究[73-83],包括结构面产状(Shanley和Mahtab,1976;Hammah,;Klose,;蔡美峰,;陈剑平,;周志芳,;唐辉明,;);结构面迹长(Hudson、Cruden、Kulatilake、Einstein、Mauldon、黄润秋、陈剑平、范留明、唐辉明);结构面间距或密度(Hudson,1976;King,1980;伍法权,2008;Cruden,;Priest和Hudson,;Kayzulovic和Goodman,;Kulatilake和Wu,;王思敬与范建军,;);结构面连通率(绪方正虔,1978;Laslett,1982;Kulatilake,1986;黄润秋和范留明,2003;陈剑平,2005;)。岩体结构面描述研究最终目标是通过这些参数来综合评估岩体结构,而岩体结构表征则是工程岩体分级中的重要指标。目前工程岩体分级主要有RMR、Q、BQ、Heok-Brown法,其中岩体结构表征主要有RQD(钻孔)、Kv(波速测试)、GSI(综合估计)、节理密度(测区)、单位体积节理数(三维模拟)等参数指标,这些指标定义与岩体结构测量方法密切相关,当然这些指标之间也存在一定的关联关系,黄润秋、陈剑平、聂德新、胡卸文、刘长武等[84-89]学者试图构建结构面密度与RQD、岩体块度指数等岩体结构表征参数之间关系,使得岩体质量指标评判指标更为合理。
1983年美国怀俄明大学的Miller采用概率论中的关联表和施密特等面积投影网结合的办法(以下简称Miller法),成功地用地质结构面产状为指标进行了岩体结构均质区划分;后来Kulatilake等对此法进行了改进,并成功运用于三峡永久船闸附近隧洞的均质区划分中;瑞典的SKB公司在进行核废料处置库场址评价时也采用了改进的Miller法;1984年美国哥伦比亚大学的Mahtab等提出了适用于地质结构面产状变化复杂的显著性簇产状相似法;Martin在划分某钻石矿场均质区时提出了相关系数法;另外,Kulatilake等及陈剑平等近年来逐步将分形维数引入均质区划分中,成功计算了施密特投影图的分形维数。上述各法所用的划分指标不尽相同,有的以结构面产状、密度、迹长等单参数作为划分依据(如Miller、Mahtab法);有些则考虑多参数的综合效应(如Kulatilake的分形维数法)。可见即便同一区域,划分法的选取不同,划分结果亦或不同。因此,针对研究区地质结构面分布特征开展均质区划分法的适宜性研究,有助于区分各法的判别力,进而评价划分结果的可信度,提升所建结构面网络模型的可靠度。对场址评价和区域稳定性分析具有重要意义。
然而,上述传统方法在概率统计岩体结构面的密度、产状、间距与迹长等参量的基础上,采用关联表、相关系数法和施密特等面积投影网相结合开展岩体结构均质区分区。不仅现场地质测绘与统计工作量巨大,而且后续数据处理流程多,其中一些步骤难以实现自动化,方法只能应用于局部小范围区域岩体结构均质区分区。
发明内容
技术问题:本发明的目的是提供一种基于图像纹理的岩体结构均质区自动分区方法,绕开传统地质统计参量,避免现场地质测量与统计的繁杂工作量,解决后续数据的自动化,从岩体结构面在露头区赋存的纹理特征出发,利用具有实际尺度特征数字图像分析技术,运用纹理参数表征岩体结构的均质区特征,采用聚类方法对子域纹理参数进行分区聚类,实现基于灰度共生矩阵与ISODATA聚类的岩体结构均质区自动分区方法。
技术方案:本发明的一种基于图像纹理的岩体结构均质区自动分区方法,采用具有实际尺寸的岩体结构纹理图像作为分区工作的数据来源,以图像中的岩体结构纹理特征为对象,采用图像纹理分析方法描述与表征岩体结构均质区,通过聚类方法实现自动化分区;规避了现场地质测绘与统计繁重工作量。
该方法包括以下步骤:
第一步:影像数据采集:借助多幅定焦摄影测量技术,利用Pix4DMapper或PhotoScan软件获得描述目标体几何特征与影像特征的三维点云数据和图像数据;针对点云数据的真实坐标尺寸和图像数据的像素影像的对应关系,建立图像数据的像素尺度;
第二步:灰度图像转换:采用特征值线性加权子空间投影法将彩色图像转换为灰度图像;
第三步:纹理信息增强:采用Hill-shading算法对灰度图像中的纹理信息进行图像增强;
第四步:纹理特征提取:采用多角度组合剖面法对灰度图像中提取出纹理特征,获得岩体结构的纹理图像;
第五步:纹理参数确定:在目测岩体结构均质区内预设纹理图像窗口,构建预设窗口的灰度共生矩阵,计算窗口内纹理图像的灰度共生矩阵10个纹理参数,通过针对灰度共生矩阵的点对距离与点对角度以及像素尺度的参数鲁棒性检验,找出鲁棒性强的纹理参数,最后采用主成分分析法在鲁棒性强的纹理参数中筛选出相对独立的三个主控纹理参数;
第六步:聚类窗口确定:在目测岩体结构均质区内预设小尺度窗口,构建预设窗口的灰度共生矩阵,计算窗口内纹理图像的三个主控纹理参数值,然后不断扩大窗口尺度,计算对应尺度窗口内纹理图像的三个主控纹理参数值,比较扩大窗口前后的主控纹理参数值,当两者方差小于3%,扩大前计算窗口尺度作为聚类窗口尺度;
第七步:纹理图像分块:按照第六步确定的聚类窗口尺度,将第一步采集的图像数据进行图像分块,获得若干图像子域;
第八步:子域纹理参数:压缩灰度图像中的灰度值为16级,构建子域的灰度共生矩阵,并归一化处理,计算所有子域的主控纹理参数;
第九步:矢量数据转换:借助ArcGIS的Fishnet工具对具备纹理参数的栅格数据转换为矢量数据;
第十步:ISODATA聚类:认定三个主控纹理参数一致或相近时同为一类岩体结构均质区,采用ISODATA聚类方法对具备三个主控纹理参数共同特征的区域进行识别、拆分与合并,实现岩体结构均质区的自动分区。
其中,
所述第一步中针对点云数据的真实坐标尺寸和图像数据的像素影像的对应关系,建立图像数据的像素尺度,具体方法为:采用地面激光扫描仪或摄影测量方法,获取目标区域的图像数据与点云数据,将具有实际坐标信息的点云数据与具有丰富影像的图像数据进行特征点位置配准,使得图像数据也具有了点云数据在图像数据平面上的投影坐标,继而图像数据实际面积除以图像数据像素数量,便可获知像素尺度。
所述第五步中计算窗口内纹理图像的灰度共生矩阵10个纹理参数包括:二阶角矩、对比度、相关度、均值和、方差、方差和、逆矩差、熵、差熵、差的方差。
所述第五步中计算的10个纹理参数,需要通过鲁棒性检验和主成分分析,筛选出三个主控纹理参数。
有益效果:传统方法不仅现场地质测绘与统计作业量大,而且数据处理流程多,尤其难以实现方法自动化,只能应用于小尺度域岩体结构均质区分区。本发明绕开传统地质统计参量,避免现场地质测量与统计的繁杂工作量,解决后续数据的自动化,从岩体结构面在露头区赋存的纹理特征出发,利用具有实际尺度特征数字图像分析技术,运用纹理参数表征岩体结构的均质区特征,采用聚类方法对子域纹理参数进行分区聚类,实现基于灰度共生矩阵与ISODATA聚类的岩体结构均质区自动分区方法。上述方法规避了现场地质测绘与统计工作量,数据处理过程可实现自动化,可应用于工程全域岩体结构均质区分区。
附图说明
图1为岩体结构纹理特征组合剖面及选点规则示意图,
图2为4×4图像窗口。
具体实施方式
下面举例对本发明的技术方案进行详细说明:
本发明的一种基于图像纹理的岩体结构均质区自动分区方法采用具有实际尺寸的岩体结构纹理图像作为分区工作的数据来源,以图像中的岩体结构纹理特征为对象,采用图像纹理分析方法描述与表征岩体结构均质区,通过聚类方法实现自动化分区;规避了现场地质测绘与统计繁重工作量。
具体方法如下:
第一步:影像数据采集,借助多幅定焦摄影测量技术,利用Pix4DMapper或Photoscan软件获得描述目标体几何特征与影像特征的三维点云数据和图像数据;针对点云数据的真实坐标尺寸和图像数据的像素影像的对应关系,建立图像数据的像素尺度。
第二步:灰度图像转换,采用特征值线性加权子空间投影法将彩色图像转换为灰度图像;
该算法步骤如下:
步骤1:输入影像图片,转YCbCr空间,将均值减去,转换公式如下:
Iycc=f(Irgb),Iycc=Iycc-Iycc,avg
其中,Iycc,avg是Iycc的均值。
步骤2:对Iycc进行主成分分析并归一化处理,λ1,λ2,λ3是其特征值,v1,v2,v3是其对应的特征向量。
步骤3:通过线性加权平均算法得到初始灰度,并将灰度值控制在值域[0,255]。
达成条件|Iy-Igray|>|Iy+Igray-255|时,满足Igray=255-Igray,转换的最终图片即为灰度图片。
第三步:纹理信息增强,采用Hill-shading算法对灰度图像中的纹理信息进行图像增强;
光照模拟常用的计算方法通过计算每个栅格网格的相对辐射值,向灰度值转换,其公式如下:
公式中分别包含相对辐射与入射辐射两种模型。G的取值区间是[0,255],0~255分别对应最暗到最亮。Gmax最大灰度值。
(3)Hf:网格坡度,0°~90°;(4)Hs:光源高度角,0°~90°;
第四步:纹理特征提取,采用多角度组合剖面法对灰度图像中提取出纹理特征,获得岩体结构的纹理图像;
岩体结构灰度图像的纹理特征借助局部分析窗口计算,单一灰度剖面很难顾及剖面线上栅格点周围邻域范围内的岩体结构纹理特征所反映出来灰度值起伏特征。如图1,箭头方向为剖面线提取方向,以平行于岩体结构纹理特征、左右以一个图像栅格为步长,依次提取N对(每对2条,于纹理剖面线左右对称)辅助剖面线,N的选择和计算分析窗口r一致,其公式为:
N=(r-1)/2
注:分析窗口为3×3时,N=1;分析窗口为5×5时,N=2。
如图1,
按照上述计算方式提取结构面剖面线及辅助剖面线,可以检测结构面局部区域内表达岩体结构纹理特征的灰度值变化的曲率、不平整程度。当剖面线与网格方向一致时,分析窗口内的数据均参与计算;当其方向随机时,参与计算的点为分析窗口与剖面线的交集。基于以上规则,定量计算其灰度值的变化,当灰度标准差较小时,岩体结构纹理特征不清晰,当灰度标准差较大时,岩体结构纹理特征比较清晰。
第五步:纹理参数确定,在目测岩体结构均质区内预设纹理图像窗口,构建预设窗口的灰度共生矩阵,计算窗口内纹理图像的灰度共生矩阵10个纹理参数,包括二阶角矩、对比度、相关度、均值和、方差、方差和、逆矩差、熵、差熵、差的方差。通过针对灰度共生矩阵的点对距离与点对角度以及像素尺度的参数鲁棒性检验,找出鲁棒性强的纹理参数,最后采用主成分分析法在鲁棒性强的纹理参数中筛选出相对独立的三个主控纹理参数;
灰度共生矩阵依据灰度图像中的灰度值,对不同灰度值出现的概率分布进行统计分析。灰度共生矩阵的构造离不开两个重要的元素,一是图像栅格点对间的方向角,二是图像栅格点对间的平面距离。通过对相同灰度的图像栅格统计,用二阶联合条件概率密度,构建灰度共生矩阵。其算法描述如下:
记坐标轴x轴的图像栅格数量为Nx,y轴的图像栅格数量为Ny,统计图像中相同位置图像栅格的灰度值,G表示图像栅格聚合后的总数量,则Ng定义为最高灰度值,公式如下所示:
Lx={1,2,......,Nx}
Ly={1,2,......,Ny}
G={1,2,......,Ng}
Lx*Ly中的每一点代表一个图像栅格,对应这个图像栅格的灰度值G,则图像中灰度值由Lx*Ly→G转化。
定义灰度共生矩阵中图像栅格点对间的方向角为θ,图像栅格点对间的平面距离为d:
Pc=p(i,j,d,θ)
其中,i表示矩阵中的行数,j表示矩阵中的列数,i、j位置的图像栅格由Pc=p(i,j,d,θ)表示,(i,j)∈G*G,i,j图像栅格点对间的平面距离为d,θ=0°,45°,90°,135°。当距离较小时,纹理图像的灰度变化较为平缓,其灰度共生矩阵对角线上的数值越大,图像栅格分布主要集中在作对角线的位置;若纹理图像灰度变化差异较大,则像元在对角线上的两侧呈均匀分布,而对角线上的数值较小。
设逆时针的方向为正,角度为θ时的图像栅格定义如下:
p(i,j,d,0°)=#{(k,l)(m,n)∈(Lx*Ly)*(Lx*Ly)|k-m=0,|l-n|=d;
f(k,l)=i,f(m,n)=j}
p(i,j,d,45°)=#{(k,l)(m,n)∈(Lx*Ly)*(Lx*Ly)|(k-m=d,l-n=d)或
(k-m=-d,l-n=-d);f(k,l)=i,f(m,n)=j}
p(i,j,d,90°)=#{(k,l)(m,n)∈(Lx*Ly)*(Lx*Ly)||k-m|=d,l-n=0;
f(k,l)=i,f(m,n)=j}
p(i,j,d,135°)=#{(k,l)(m,n)∈(Lx*Ly)*(Lx*Ly)|(k-m=d,l-n=-d)}
或(k-m=-d,l-n=d);f(k,l)=i,f(m,n)=j
其中,#{x}表示x的几何数,公式表述了矩阵里θ方向所有距离为d的栅格。
Pc=p(i,j,d,θ)除以一个常系数R,即p(i,j)/R→p(i,j),这个常数叫正规化常数,提高纹理的分辨率。其中d=1,θ=0°,R=2Ny(Nx-1),当d=1,θ=45°,R=2(Ny-1)(Nx-1)。当θ=90°,135°时,其相邻个数满足的以下公式:
设一个4×4的图像窗口如图2所示,构建灰度共生矩阵统计各个方向的灰度值出
现的次数,则0°,45°,90°,135°方向的矩阵如下:0°,45°,90°,135°灰度矩阵,
第六步:聚类窗口确定,在目测岩体结构均质区内预设小尺度窗口,构建预设窗口的灰度共生矩阵,计算窗口内纹理图像的三个主控纹理参数值,然后不断扩大窗口尺度,计算对应尺度窗口内纹理图像的三个主控纹理参数值,比较扩大窗口前后的主控纹理参数值,当两者方差小于3%,扩大前计算窗口尺度作为聚类窗口尺度。
第七步:纹理图像分块,按照第六步确定的聚类窗口尺度,将第一步采集的图像数据进行图像分块,获得若干图像子域。
基于Matlab实现的主要代码如下:
I=imread(′3-3.jpg′);%%实现数据载入,读取图片各像素的灰度值,构建double矩阵。
rs=size(I,1);cs=size(I,2);%%设置图像分块的行与列
sz=75;%%设置分块后的图像像素。
numr=rs/sz;
numc=cs/sz;
ch=sz;cw=sz;
t1=(0:numr-1)*ch+1;t2=(1:numr)*ch;%%计算分块图像的初始行像素值。
t3=(0:numc-1)*cw+1;t4=(1:numc)*cw;%%计算分块图像的初始列像素值。
第八步:子域纹理参数,压缩灰度图像中的灰度值为16级,构建子域的灰度共生矩阵,并归一化处理。计算所有子域的主控纹理参数。
纹理特征参数值的主要计算代码如下:
Gray=Gray/16;%%灰度级压缩为16级。
Gray(i,j)==m-1&&Gray(i,j+3)==n-1;%%点对距离为3,0°方向的GLCM矩阵。
Gray(i,j)==m-1&&Gray(i-3,j+3)==n-1;%%点对距离为3,45°方向的GLCM矩阵。
Gray(i,j)==m-1&&Gray(i+3,j)==n-1;%%点对距离为3,90°方向的GLCM矩阵。
Gray(i,j)==m-1&&Gray(i+3,j+3)==n-1;%%点对距离为3,135°方向的GLCM矩阵。
P(:,:,n)=P(:,:,n)/sum(sum(P(:,:,n)));%%归一化处理
主要函数代码:
H(n)=-P(i,j,n)*log(P(i,j,n))+H(n);
a2=mean(H);%%GLCM矩阵中4个方向上熵的均值。
I(n)=P(i,j,n)/(1+(i-j)^2)+I(n);
a3=mean(I);%%GLCM矩阵中4个方向上逆差矩的均值。
Ux(n)=i*P(i,j,n)+Ux(n);%相关度中参数μx
Uy(n)=j*P(i,j,n)+Uy(n);%相关度中参数μy
U(n)=(i+j)*P(i,j,n)+U(n);
a4=mean(C);%%GLCM矩阵中4个方向上相关度的均值。
结合图像分块后的像素大小,采用滑动窗口法读取每个分块的特征值,其滑窗法读取源代码如下:
A=imread(′3-4.jpg′);
[m,n]=size(A);
for k=1:m
for kk=1:n
B=A(k:k+75,kk:kk+75)
mwrite(B,[′3-4jpg′]);
第九步:矢量数据转换,借助ArcGIS的Fishnet工具对具备纹理参数的栅格数据转换为矢量数据。
第十步:ISODATA聚类,认定三个主控纹理参数一致或相近时同为一类岩体结构均质区,采用ISODATA聚类方法对具备三个主控纹理参数共同特征的区域进行识别、拆分与合并,实现岩体结构均质区的自动分区。
ISODATA算法的详细步骤如下所示:
步骤1:输入N个点的数据集合{xi,i=1,2,…,N};在这项研究中,输入的数据集合{xi,i=1,2,…,N}被分别设置为三个数据层,包括栅格网格的相关度,逆差矩和熵。N是网格的数量。设置Nc为初始聚类中心没有必要使Nc的值等于聚类的数量,中心的初始值可以随机选择。K:聚类中心的预期数量;θn:可以形成聚类的最少点数。如果一个聚类中的点数少于θn个,它将与另一个聚类相结合。
θs:沿着每个轴线的聚类中心点的最大标准偏差;
θc:两个聚类中心之间的最小距离。如果两个聚类的距离小于θc将合并;
L:在一次迭代中可以合并的最大数量的聚类对数;
I:最多迭代计算的次数;
步骤2:将数据集分配给聚类。如果Di=min{||xi-yi||,j=1,2,…,Nc},则xi∈Sj。
步骤3:如果点数Sj小于θn,省略样本的子集;同时Nc减1。
步骤4:重新计算聚类中心。
其中Nj是Sj中的样本的数量。
步骤7:确定是否需要拆分或合并。如果迭代计算次数已经达到I,则设,θc=0进入步骤11。如果Nc≤K/2,则进入步骤8,并拆分聚类。
步骤8:对于每个聚类Sj,计算向量σj。每个σij代表坐标i点的σj直接从zj指向聚类Sj中的每个点的标准差。
σj=(σ1j,σ2j,…σnj)T
其中i(i=1,2,…,n)是样本特征向量的维数,j(j=1,2,…,Nc)是聚类的数量,Nj和Sj是样本的数量。
步骤9:令σjmax=max(σ1,σ2,σ3,…,σj)。
步骤10:对于每一个聚类Sj,如果满足σjmax>θs,并且满足(和Nj>2(θN+1))或者Nc≤K/2,就根据σjmax计算两个新的聚类,用这两个新的聚类替换zj,增加k,并将Sj分成两个聚类。如果激活了分割操作,进入步骤2。
步骤11:计算所有不同的聚类中心之间成对的聚类间距离。
Dij=||zi-zj||,i=1,2,…,Nc-1;j=i+1,…,Nc
步骤12:比较Dij和θc,并选择全部Dij小于θc作为一个新的子集。然后按升序对子集的元素进行排序:
{Di1j1,Di2j2,…,DiLjL}
其中Di1j1<Di2j2<…<DiLjL
步骤13:结合两个聚类中心zik和zjk,然后计算新的聚类中心:
步骤14:如果这是最后一次迭代(即I次),则算法结束。否则,如果需要更改输入参数,则返回步骤1;否则返回步骤2,通过加1来更新迭代。
Claims (4)
1.一种基于图像纹理的岩体结构均质区自动分区方法,其特征在于该方法采用具有实际尺寸的岩体结构纹理图像作为分区工作的数据来源,以图像中的岩体结构纹理特征为对象,采用图像纹理分析方法描述与表征岩体结构均质区,通过聚类方法实现自动化分区;规避了现场地质测绘与统计繁重工作量;
该方法包括以下步骤:
第一步:影像数据采集:借助多幅定焦摄影测量技术,利用Pix4DMapper或PhotoScan软件获得描述目标体几何特征与影像特征的三维点云数据和图像数据;针对点云数据的真实坐标尺寸和图像数据的像素影像的对应关系,建立图像数据的像素尺度;
第二步:灰度图像转换:采用特征值线性加权子空间投影法将彩色图像转换为灰度图像;
第三步:纹理信息增强:采用Hill-shading算法对灰度图像中的纹理信息进行图像增强;
第四步:纹理特征提取:采用多角度组合剖面法对灰度图像中提取出纹理特征,获得岩体结构的纹理图像;
第五步:纹理参数确定:在目测岩体结构均质区内预设纹理图像窗口,构建预设窗口的灰度共生矩阵,计算窗口内纹理图像的灰度共生矩阵10个纹理参数,通过针对灰度共生矩阵的点对距离与点对角度以及像素尺度的参数鲁棒性检验,找出鲁棒性强的纹理参数,最后采用主成分分析法在鲁棒性强的纹理参数中筛选出相对独立的三个主控纹理参数;
第六步:聚类窗口确定:在目测岩体结构均质区内预设小尺度窗口,构建预设窗口的灰度共生矩阵,计算窗口内纹理图像的三个主控纹理参数值,然后不断扩大窗口尺度,计算对应尺度窗口内纹理图像的三个主控纹理参数值,比较扩大窗口前后的主控纹理参数值,当两者方差小于3%,扩大前计算窗口尺度作为聚类窗口尺度;
第七步:纹理图像分块:按照第六步确定的聚类窗口尺度,将第一步采集的图像数据进行图像分块,获得若干图像子域;
第八步:子域纹理参数:压缩灰度图像中的灰度值为16级,构建子域的灰度共生矩阵,并归一化处理,计算所有子域的主控纹理参数;
第九步:矢量数据转换:借助ArcGIS的Fishnet工具对具备纹理参数的栅格数据转换为矢量数据;
第十步:ISODATA聚类:认定三个主控纹理参数一致或相近时同为一类岩体结构均质区,采用ISODATA聚类方法对具备三个主控纹理参数共同特征的区域进行识别、拆分与合并,实现岩体结构均质区的自动分区。
2.根据权利要求1所述的一种基于图像纹理的岩体结构均质区自动分区方法,其特征在于所述第一步中针对点云数据的真实坐标尺寸和图像数据的像素影像的对应关系,建立图像数据的像素尺度,具体方法为:采用地面激光扫描仪或摄影测量方法,获取目标区域的图像数据与点云数据,将具有实际坐标信息的点云数据与具有丰富影像的图像数据进行特征点位置配准,使得图像数据也具有了点云数据在图像数据平面上的投影坐标,继而图像数据实际面积除以图像数据像素数量,便可获知像素尺度。
3.根据权利要求1所述的一种基于图像纹理的岩体结构均质区自动分区方法,其特征在于所述第五步中计算窗口内纹理图像的灰度共生矩阵10个纹理参数包括:二阶角矩、对比度、相关度、均值和、方差、方差和、逆矩差、熵、差熵、差的方差。
4.根据权利要求1所述的一种基于图像纹理的岩体结构均质区自动分区方法,其特征在于所述第五步中计算的10个纹理参数,需要通过鲁棒性检验和主成分分析,筛选出三个主控纹理参数。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910437094.1A CN110135515B (zh) | 2019-05-23 | 2019-05-23 | 一种基于图像纹理的岩体结构均质区自动分区方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910437094.1A CN110135515B (zh) | 2019-05-23 | 2019-05-23 | 一种基于图像纹理的岩体结构均质区自动分区方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110135515A CN110135515A (zh) | 2019-08-16 |
CN110135515B true CN110135515B (zh) | 2021-04-27 |
Family
ID=67572897
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910437094.1A Active CN110135515B (zh) | 2019-05-23 | 2019-05-23 | 一种基于图像纹理的岩体结构均质区自动分区方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110135515B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111553292B (zh) * | 2020-04-30 | 2023-05-02 | 南京工业大学 | 一种基于点云数据的岩体结构面识别与产状分类方法 |
CN111797679A (zh) * | 2020-05-19 | 2020-10-20 | 中国地质大学(武汉) | 一种遥感纹理信息处理方法、装置、终端及存储介质 |
CN114882038B (zh) * | 2022-07-12 | 2022-09-30 | 济宁鸿启建设工程检测有限公司 | 一种建筑外墙保温类材料检测方法及检测设备 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2013163467A1 (en) * | 2012-04-26 | 2013-10-31 | Arizona Chemical Company, Llc | Rejuvenation of reclaimed asphalt |
CN103529455A (zh) * | 2013-10-21 | 2014-01-22 | 中铁第四勘察设计院集团有限公司 | 一种基于机载激光雷达三维的危岩落石调查方法 |
CN105180890A (zh) * | 2015-07-28 | 2015-12-23 | 南京工业大学 | 融合激光点云和数字影像的岩体结构面产状测量方法 |
Family Cites Families (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN102081791B (zh) * | 2010-11-25 | 2012-07-04 | 西北工业大学 | 一种基于多尺度特征融合的sar图像分割方法 |
CN104268836A (zh) * | 2014-09-24 | 2015-01-07 | 江西理工大学 | 一种基于局域均质指标的分水岭分割标记点提取方法 |
CN105893723A (zh) * | 2014-10-15 | 2016-08-24 | 长沙矿山研究院有限责任公司 | 一种基于微震事件簇群pca法岩体断层滑移面产状计算方法 |
CN106225770B (zh) * | 2016-08-26 | 2018-12-25 | 招商局重庆交通科研设计院有限公司 | 隧道掌子面地质多维数字化记录识别方法及系统 |
CN106447776A (zh) * | 2016-09-22 | 2017-02-22 | 北京科技大学 | 基于3d打印制作的复杂裂隙岩体物理模型及建模方法 |
TWI626622B (zh) * | 2017-07-04 | 2018-06-11 | System and method for stereoscopic imaging of underground rock formation characteristics | |
CN107703560B (zh) * | 2017-09-29 | 2019-12-13 | 西南石油大学 | 一种基于三重信息的泥页岩岩相精细识别方法 |
-
2019
- 2019-05-23 CN CN201910437094.1A patent/CN110135515B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2013163467A1 (en) * | 2012-04-26 | 2013-10-31 | Arizona Chemical Company, Llc | Rejuvenation of reclaimed asphalt |
CN103529455A (zh) * | 2013-10-21 | 2014-01-22 | 中铁第四勘察设计院集团有限公司 | 一种基于机载激光雷达三维的危岩落石调查方法 |
CN105180890A (zh) * | 2015-07-28 | 2015-12-23 | 南京工业大学 | 融合激光点云和数字影像的岩体结构面产状测量方法 |
Also Published As
Publication number | Publication date |
---|---|
CN110135515A (zh) | 2019-08-16 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Chen et al. | Machine learning-based classification of rock discontinuity trace: SMOTE oversampling integrated with GBT ensemble learning | |
Guo et al. | Towards semi-automatic rock mass discontinuity orientation and set analysis from 3D point clouds | |
Xu et al. | Extraction and statistics of discontinuity orientation and trace length from typical fractured rock mass: A case study of the Xinchang underground research laboratory site, China | |
CN110135515B (zh) | 一种基于图像纹理的岩体结构均质区自动分区方法 | |
Chen et al. | Towards semi-automatic discontinuity characterization in rock tunnel faces using 3D point clouds | |
CN103345566B (zh) | 基于地质内涵的化探异常识别与评价方法 | |
US20170018073A1 (en) | Digital Rock Physics-Based Trend Determination and Usage for Upscaling | |
Mohebbi et al. | Rock mass structural data analysis using image processing techniques (Case study: Choghart iron ore mine northern slopes) | |
CN108038081B (zh) | 基于特征函数空间滤值的滑坡灾害logistic回归分析方法 | |
CN110363299B (zh) | 面向露头岩层分层的空间案例推理方法 | |
Ge et al. | Rock discontinuities identification from 3D point clouds using artificial neural network | |
CN108830317B (zh) | 基于数字摄影测量的露天矿山边坡岩体节理产状快速精细取值方法 | |
Ge et al. | Estimation of the appropriate sampling interval for rock joints roughness using laser scanning | |
Jung et al. | A line-based progressive refinement of 3D rooftop models using airborne LiDAR data with single view imagery | |
Chen et al. | A novel image-based approach for interactive characterization of rock fracture spacing in a tunnel face | |
Chen et al. | A critical review of automated extraction of rock mass parameters using 3D point cloud data | |
Yi et al. | An efficient method for extracting and clustering rock mass discontinuities from 3D point clouds | |
Zahs et al. | Classification of structural building damage grades from multi-temporal photogrammetric point clouds using a machine learning model trained on virtual laser scanning data | |
CN112200417A (zh) | 基于摄影测量、BQ、RQDt、地应力的改进Mathews稳定图评价方法 | |
Michalak et al. | Using Delaunay triangulation and cluster analysis to determine the orientation of a sub-horizontal and noise including contact in Kraków-Silesian Homocline, Poland | |
Ge et al. | Rock joint detection from 3D point clouds based on colour space | |
Van Knapen et al. | Identification and characterisation of rock mass discontinuity sets using 3D laser scanning | |
Wijaya et al. | Building crack due to lombok earthquake classification based on glcm features and svm classifier | |
Liu et al. | LOSN: lightweight ore sorting networks for edge device environment | |
CN115578212A (zh) | 一种地质勘查内外业一体化的信息化系统 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |