一种基于图像纹理的岩体结构均质区自动分区方法与流程

本发明涉及工程现场岩体结构均质区的自动分区方法。属于岩体工程的现场地质调查技术领域。适用于岩体工程全域范围内岩体地质分区工作。

背景技术：

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最大灰度值。

g的控制参数：(1)网格的坡向，0°～360°；(2)光源方位角，0°～360°；

(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)；

fork＝1：m

forkk＝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中的样本的数量。

步骤5：令为sj与相关联的聚类中心zj之间的平均距离

步骤6：是这些距离的总体平均值。

步骤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来更新迭代。