一种用于岩体节理裂隙模型检验的迹线图相似度判别方法与流程

本发明涉及岩土工程、地质三维建模分析以及计算机图形学领域，特别涉及一种用于岩体节理裂隙模型检验的迹线图相似度判别方法。

背景技术：

岩体结构面反映了岩体不连续、非均匀、各向异性的本质，并对岩体的力学特性、水力学特性和工程稳定性起着关键控制作用。然而，岩体内部的节理裂隙难以直接被观测到，且数量众多，这给岩体稳定性分析带来了极大的挑战。同时，随着计算机技术的发展，地质三维建模技术也越来越成熟，而在地层中加入离散裂隙网络的岩体精细化建模也逐渐成为一种趋势，且已被证明是解决上述问题的一种有效的方法，如熊自明等人公开了一种基于分形算法的区域地质体岩体裂隙模拟的方法；吴云公开了一种适用于岩体的建模方法，并能够对更复杂的破坏模式进行力学的稳定性分析。

然而，对于目前的精细化建模过程，其中的节理裂隙的各种几何参数都是通过随机模拟得到的，模拟结果是否有效，构建的模型是否能够代表岩体内部真实的裂隙，是裂隙模拟成功与否的关键。对于模型的有效性检验，目前主要存在两种方法：数值检验和图形检验。其中数值检验较有依据性，其主要的内容是运用统计学和概率论的方法将实际节理裂隙的各种特征与模拟生成的节理裂隙进行对比分析。一方面，以差分比的方式对模型数据与实际数据对进行均值误差分析；另一方面，判断模型数据的离散性是否达到要求。而对于图形检验，则存在较大的主观性。其方法就是简单地将实际工程勘察中的迹线素描图与模型中对应位置生成的迹线图进行对比，通过肉眼直观地对其进行相似与否的评价。这种方法主观性太强，得出的结果因人而异，难以达成一个统一的标准；另一方面，图形检验又拥有的数值检验所不能替代的特性，如通过对图形进行分析，能从宏观到微观综合反映出节理裂隙的密度、长度、间距甚至相对位置等信息。

技术实现要素：

本发明的目的是克服现有技术中的不足，提供一种用于岩体节理裂隙模型检验的迹线图相似度判别方法，将计算机图像处理技术、三维地质建模技术与岩土工程专业特点相结合，克服了以往三维地形建模有效性检验过程中图形检验的主观性，为地质建模的有效性检验提供了定量的标准和客观的依据。

本发明所采用的技术方案是：一种用于岩体节理裂隙模型检验的迹线图相似度判别方法，包括以下步骤：

步骤a：列出原始的实际迹线图和模拟迹线图；

步骤b：确定相似阈值，作为下述步骤中判断相似性的基本标准；

步骤c：原始的实际迹线图和模拟迹线图的预处理；

步骤d：将预处理后的实际迹线图和模拟迹线图的总体灰度值相似度与相似阈值进行比较，判断迹线总量相似度；

步骤e：如果步骤d中判别为相似，将两幅迹线图的灰度级配曲线相似度与相似阈值进行比较，判断迹线图局部特点相似度；

步骤f：如果步骤e中判别为相似，对两幅迹线图进行特征方向相似度匹配；

步骤g：如果步骤f中判别为相似，将两幅迹线图的灰度密度分布曲线相似度与相似阈值进行比较，完成用于岩体节理裂隙模型检验的迹线图相似度判别。

步骤b中的相似阈值根据实际工程精度要求设置在[0.7,0.9]的范围。

步骤c进一步包括：

步骤c1：如果实际迹线图或模拟迹线图为彩色图像，将其转换为灰度图；

步骤c2：强化迹线清晰度，增强迹线可识别性；

步骤c3：对两幅迹线图进行去边框处理，以避免算法将边框识别为迹线的一部分；

步骤c4：以裁剪的方式对两幅迹线图进行微调，统一迹线图尺寸。

步骤d的过程为：分别求出预处理后的实际迹线图与模拟迹线图的总体灰度值，计算其相似度并与相似阈值进行比较，如果大于等于相似阈值，则进行下一步分析，如果小于相似阈值，就认为这两幅迹线图是不相似的。

步骤e进一步包括：

步骤e1：对两幅预处理后的迹线图进行网格划分，求得每个网格的平均灰度值；

步骤e2：制作出不同孔径的筛网，筛网按孔径由小到大排列；

步骤e3：以平均灰度值为标准，将每个网格类比成直径不同的“土质颗粒”，通过筛网对这堆“土质颗粒”进行筛选，分别绘制两幅迹线图的灰度级配曲线；

步骤e4：计算预处理后的实际迹线图的灰度级配曲线和模拟迹线图的灰度级配曲线的相似度，比将其与相似阈值进行比较，如果大于等于相似阈值，则进行下一步分析，如果小于相似阈值，就认为这两幅迹线图是不相似的。

步骤f进一步包括：

步骤f1：分别对两幅预处理后的迹线图进行拉东变换，变换的强度为15°一次；

步骤f2：在拉东变换过程中，将灰度的线性积分结果进一步转变为灰度密度分布曲线；

步骤f3：对拉东变换过程中得到的所有的灰度密度分布曲线进行离散性分析，找出离散性最大的前四个曲线，并规定其对应的测线方向为特征方向；

步骤f4：比较两幅图特征方向的相似性，如果有至少两个特征方向是一致的，则进行下一步分析，如果不是，就认为这两幅迹线图是不相似的。

步骤g进一步包括：

步骤g1：以步骤f中得到的实际迹线图的四个特征方向为标准，分别提取实际迹线图和模拟迹线图的在这四个特征方向上的灰度密度分布曲线，并一一对应，得到四对灰度密度分布曲线；

步骤g2：对于每一对灰度密度分布曲线，进行回环滚动余弦相似度匹配，并提取出匹配过程中得到的最大值；

步骤g3：对步骤g2中得到了的四个最大匹配值按照加权平均的计算方法计算相似度，将其与相似阈值进行比较，如果大于等于相似阈值，就认为这两幅迹线图是相似的，如果小于相似阈值，就认为这两幅迹线图是不相似的。

本发明的有益效果是：本发明结合了计算机图形学和数学地质学，克服了以往三维地形建模有效性检验过程中图形检验主观性太强的缺点，将人对图像直观的视觉分析过程转化为计算机算法，为模型的有效性检验提供了合理的标准和充分的依据。

附图说明

图1：本发明一种用于岩体节理裂隙模型检验的迹线图相似度判别方法流程图；

图2：本发明工程实例实际迹线图；

图3：本发明工程实例模拟迹线图；

图4：本发明工程实例预实际迹线图预处理；

图5：本发明工程实例预模拟迹线图预处理；

图6：本发明工程实例实际迹线图总体灰度值统计图；

图7：本发明工程实例模拟迹线图总体灰度值统计图；

图8：本发明工程实例灰度级配曲线图；

图9：本发明工程实例特征方向对比图；

图10：本发明工程实例特征方向相似度计算示意图；

图11：本发明工程实例滚动余弦相似度。

具体实施方式

下面结合附图以及具体公式对本发明作进一步的描述。

如附图1所示，一种用于岩体节理裂隙模型检验的迹线图相似度判别方法，包括以下步骤：

步骤a：列出原始的实际迹线图和模拟迹线图；

步骤b：确定相似阈值r，作为下述步骤中判断相似性的基本标准，所述相似阈值r根据实际工程精度要求设置在[0.7,0.9]的范围，如设定r＝0.85，则代表在下述步骤中产生的相关指标均达到了0.85，才能将两幅迹线图判定为相似。

步骤c：原始的实际迹线图和模拟迹线图的预处理。具体包括：

步骤c1：如果实际迹线图或模拟迹线图为彩色图像，将其转换为灰度图。

步骤c2：强化迹线清晰度，增强迹线可识别性：将灰度值低于或等于230(230≈255×0.9)的像素点统一设置为0，将灰度值高于230的像素点统一设置为255；其中，在处理图像的时候，255是白和0是黑，0～255之间是灰色；0.9是用户自已定义的一个界限。

步骤c3：对两幅迹线图进行去边框处理，以避免算法将边框识别为迹线的一部分：从上到下逐行搜索灰度矩阵，当搜索到某一行的灰度值骤减，记录该行到上边缘的距离，此距离大小则为边框的厚度，并以此厚度为标准，去掉四个边框线段。

步骤c4：以裁剪的方式对两幅迹线图进行微调，统一迹线图尺寸。通常情况下，实际迹线图与模拟迹线图大小基本一致，但可能有很细微的差别；而在算法处理的过程对两幅迹线图的要求是尺寸严格一致。具体处理的方法是：将两张图像重叠，中心点重合，且保持横向中心轴与纵向中心轴皆重合，将各自存在的不重叠的部分裁剪掉。

步骤d：将预处理后的实际迹线图和模拟迹线图的总体灰度值相似度η与相似阈值r进行比较，判断迹线总量相似度。具体过程为：

经预处理后的迹线图为白底黑线，其中黑线部分代表全部的节理裂隙，分别统计两幅迹线图的总体灰度值sa和sb，并通过z分数计算它们的相似度η：

式中，sa为实际迹线图的总体灰度值，sb为模拟迹线图的总体灰度值。

将相似度η与相似阈值r进行比较：如果η≥r，就认为这两幅迹线图是相似的，则继续进行下述步骤；如果η<r，就认为这两幅迹线图是不相似的。

步骤e：如果步骤d中判别为相似，将两幅迹线图的灰度级配曲线相似度与相似阈值进行比较，判断迹线图局部特点相似度。具体包括：

步骤e1：对两幅预处理后的迹线图进行网格划分，求得每个网格的平均灰度值：参考其它图像处理算法中图像分级处理的方法，分别将两幅图划分成2ⁿ×2ⁿ的均等网格(考虑到网格数量不能太少，同时又不能太大以致于使计算量太大，优选的，将迹线图划分成16×16的均等网格)，计算每个网格的平均灰度值dij。

步骤e2：制作出不同孔径的筛网，筛网按孔径由小到大排列：在两幅迹线图像的网格中，找到平均灰度值dij的最大值dmax作为筛网的最大筛网孔径，在[0,dmax]范围内，优选的以0.005为增量划分出一系列孔径不同的筛网，并由小到大排列。

步骤e3：以平均灰度值dij为标准，将每个网格类比成直径不同的“土质颗粒”，通过筛网对这堆“土质颗粒”进行筛选，分别绘制两幅迹线图的灰度级配曲线[a1,a2,...,an]和[b1,b2,...,bn]。

步骤e4：计算实际迹线图的灰度级配曲线[a1,a2,...,an]和模拟迹线图的灰度级配曲线[b1,b2,...,bn]的相似度ξ：

将相似度ξ与相似阈值r进行比较：如果ξ≥r，就认为这两幅迹线图是相似的，则继续进行下述步骤；如果ξ<r，就认为这两幅迹线图是不相似的。

步骤f：如果步骤e中判别为相似，对两幅迹线图进行特征方向相似度匹配。具体包括：

步骤f1：分别对两幅预处理后的迹线图进行拉东变换(radon变换)，变换的强度为15°一次。

radon变换可定义在做任意维变量空间域，对于图像面言，考虑二维欧式空间的变换形式：

式中：d为整个图像平面；f(x,y)为图像上某一点(x,y)的像素点灰度值；特征函数δ为狄拉克函数；ρ为(x,y)平面内直线到原点的距离；θ为原点到直线的垂线与x轴的夹角，特征函数δ使积分沿着直线ρ＝xcosθ-ysinθ进行。此处，令θ按逆时针每15°变换一次，ρ每3个单位变换一次。

步骤f2：在拉东变换过程中，将灰度的线性积分结果进一步转变为灰度密度分布曲线。即：将原有的radon变换公式修改为如下：

式中各项含意同步骤f1。

步骤f3：对拉东变换过程中得到的所有的灰度密度分布曲线进行离散性分析，找出离散性最大的前四个曲线，并规定其对应的测线方向为图像的特征方向。即：对于每一个θ，生成由不同的ρ产生的曲线，并求得其方差，其中方差最大的前四个灰度密度分布曲线所对应的θ被定义为该迹线图的特征方向。

步骤f4：比较两幅迹线图特征方向的相似性λ，包括重现率与顺序相似性判断。其中，重现率指的是模拟图的四个特征方向中有多少个是同时存在实际迹线图的特征方向中，用符号λ1表示，λ1的取值为[0,0.25,0.5,0.75,1]。顺序相似度是表征了当有几个特征方向既属于模拟迹线图又属于实际迹线图时，它们在模拟迹线图的排列顺序是否与实际迹线图相似，用λ2表示。具体的计算方法为：

(1)设[rd1,rd2,rd3,rd4]为实际迹线图的特征方向，[sd1,sd2,sd3,sd4]为模拟迹线图的特征方向。然后建立矩阵m：

(2)λ1可按下式计算：

(3)找到m矩阵中所有值为1的元素的坐标，并记录坐标，如下式所示：

p＝{(i,j)|i,j∈z,0≤i,j≤4,mij＝1}(7)

在平面坐标中，将p中所有的点两两相连，得到一组向量v(v＝v1,v2,…)。数出所有斜率为正值的vi的个数，并以该数比上v中向量的总数，得到λ2。

对于得到的λ1和λ2，采用“二八法则”，按下式计算出综合的值：

λ＝0.8λ1+0.2λ2(8)

当λ≥0.6时，则可以继续进行下述步骤；如果不是，就直接认为这两幅迹线图不相似。

步骤g：如果步骤f中判别为相似，将两幅迹线图的灰度密度分布曲线相似度与相似阈值进行比较，完成用于岩体节理裂隙模型检验的迹线图相似度判别。具体包括：

步骤g1：对于实际迹线图，提取步骤f3过程中产生的四个特征方向对应的四个灰度密度分布曲线；以实际迹线图的四个特征方向为标准，提取模拟迹线图的在这四个特征方向上的灰度密度分布曲线，并一一对应，得到四对灰度密度分布曲线。

步骤g2：对于每一对灰度密度分布曲线，进行回环滚动余弦相似度匹配，并提取出匹配过程中得到的最大值。以一对曲线[x1,x2,...,xn]和[y1,y2,...,yn]为例，余弦相似度公式如下：

式中，xi为实际迹线图灰度密度分布曲线的第i个值，yi为模拟迹线图灰度密度分布曲线的第i个值。其中“回环滚动”的方法为，分别将两条灰度密度分布曲线各自首尾相连，围成两个曲线环，重合两个环的中心轴并旋转错动，同时边错动边用余弦相似度去匹配，最终滚动匹配过程中产生最大余弦相似度的值，即最大匹配值。

步骤g3：对步骤g2中得到了的四个最大匹配值从小到大排列为[α1,α2,α3,α4]，按照加权平均的计算方法计算相似度ζ：

ζ＝0.1α1+0.2α2+0.3α3+0.4α4(10)

将相似度ζ与相似阈值r进行比较，如果ζ≥r，就认为这两幅迹线图是相似的，如果ζ＜r，就认为这两幅迹线图是不相似的，完成用于岩体节理裂隙模型检验的迹线图相似度判别。

下面我们将借助一个工程实例对本方法进行验证，整个算法通过matlab实现：

图2为某工程地质勘察中平硐勘察的部分实际迹线图，图3为三维地质模型中对应位置的模拟迹线图，其中不同颜色代表不同分组的优势节理裂隙。现对模拟结果的有效性进行图形检验：

(1)取阈值r＝0.85。

(2)对两幅迹线图进行预处理后得到的结果如图4、图5所示。

(3)对预处理后的两幅迹线图进行总体灰度值统计，如图6、图7所示，其中实际迹线图与模拟迹线图总体灰度值相似值η＝0.9148>r。

(4)两幅迹线图的灰度级配曲线分析如图8所示，其中dmax＝0.375，ξ＝1>r。

(5)特征方向分析结果如图9所示，实际迹线图的特征方向为[105°,120°,135°,150°]，模拟迹线图的特征方向为[120°,135°,105°,15°]。以此建立的m矩阵及重现率和顺序相似度的计算结果如图10所示。

重现率为：λ1＝3/4＝0.75；

顺序相似度为：λ2＝1/3＝0.33

特征方向相似度为：λ＝0.8×0.75+0.2×0.33＝0.67>0.6。

(6)回环滚动余弦相似度的结果如图11所示，ζ＝0.8506＞r。

综上，模拟迹线图与实际迹线相似，图形检验合格。