一种基于差商分析的土裂隙分形维数计算方法技术

本发明专利技术涉及一种基于差商分析的土裂隙分形维数计算方法，包括以下步骤：1)对原始彩色裂隙图片进行灰度化和二值化处理，得到裂隙二值图；2)采用多个不同边长r的正方形网格覆盖含有裂隙的网络区域，统计裂隙穿过的正方形的累计个数N；3)对多个边长‑个数对(r,N)取自然对数，得到(ln(1/r)，ln(N))序列，在坐标系下绘制ln(N)‑ln(1/r)的散点；4)基于差商分析法获取同时满足ln(N)‑ln(1/r)一阶差商值在1～2之间缓慢递减和二阶差商绝对值小于1两个条件的区间作为裂隙网络的无标度区；5)在无标度区内采用基于最小二乘法用直线拟合ln(N)‑ln(1/r)散点，得到的直线斜率即为土裂隙的分形盒维数。与现有技术相比，本发明专利技术具有首创性、计算可靠、标准可行等优点。

A Method for Calculating Fractal Dimensions of Soil Fractures Based on Differential Quotient Analysis

The present invention relates to a method for calculating fractal dimension of soil fissures based on difference quotient analysis, which includes the following steps: 1) graying and binarizing the original color fissure images to obtain fissure binary images; 2) covering the network area with fissures with square grids of different side lengths r, counting the cumulative number of square fissures passing through N; 3) counting multiple side lengths (r) pairs (r) (N) Take natural logarithm to get (ln(1/r), ln(N)) sequence and draw the scatter point of ln(N)ln(1/r) in coordinate system; 4) Obtain the interval that satisfies the conditions of slow decline of ln(N)ln(1/r) first-order difference quotient between 1 and 2 and absolute value of second-order difference quotient less than 1 or 2 as scale-free zone of fracture network based on least squares method in scale-free zone. The linear slope of n(N)ln(1/r) scatter points is the fractal box dimension of soil fissures. Compared with the prior art, the present invention has the advantages of originality, reliable calculation and feasible standard.

【技术实现步骤摘要】
一种基于差商分析的土裂隙分形维数计算方法
本专利技术涉及岩土工程
，尤其是涉及一种基于差商分析的土裂隙分形维数计算方法。
技术介绍
土体表面干缩开裂是自然界的一种常见现象。裂隙的开展会影响到土体的压缩性、强度以及渗透性等工程性质。同时，土裂隙网络的形态特征在一定程度上反映了土体失水过程中表面和内部应变场的演化规律。因此，裂隙表面几何形态的定量分析对研究裂隙发育机理具有重要意义，也为土体开裂现象引起的边坡失稳、地基破坏等工程事故的预防和处理提供研究基础。选择合适的评价指标是进行土裂隙定量评价的前提。在众多裂隙评价指标中，分形维数常被用来表征裂隙几何特征的复杂程度。然而，自然界中的分形只是统计意义上的自相似，分形特征只存在于一定尺度范围，即所谓的“无标度区”。在分形研究中，无标度区间的确定是一个难点，至今仍然缺乏一个统一的、公认的标准用于判定无标度区。以往的研究中，土裂隙网络分形维数的计算常常忽略无标度区，这可能会使计算出的分形维数不确定、不可靠。
技术实现思路
本专利技术的目的就是为了克服上述现有技术存在的缺陷而提供一种基于差商分析的土裂隙分形维数计算方法。本专利技术的目的可以通过以下技术方案来实现：一种基于差商分析的土裂隙分形维数计算方法，包括以下步骤：1)对原始彩色裂隙图片进行灰度化和二值化处理，得到裂隙二值图，提取土体干缩裂隙；2)根据数格子法，采用多个不同边长r的正方形网格覆盖含有裂隙的网络区域，统计裂隙穿过的正方形的累计个数N；3)对多个边长-个数对(r,N)取自然对数，得到(ln(1/r)，ln(N))序列，在坐标系下绘制ln(N)-ln(1/r)的散点；4)基于差商分析法获取同时满足ln(N)-ln(1/r)一阶差商值在1～2之间缓慢递减和二阶差商绝对值小于1两个条件的区间作为裂隙网络的无标度区；5)在无标度区内采用基于最小二乘法用直线拟合ln(N)-ln(1/r)散点，得到的直线斜率即为土裂隙的分形盒维数。所述的步骤1)中，在MATLAB软件中采用rgb2gray函数获取像素点的灰度值，得到原始裂隙图片的灰度图。所述的步骤1)中，采用局部阈值分割法进行裂隙图像二值化：先将灰度图分为n2个区域，对每个区域分别用大津法进行阈值分割，提取土裂隙部分，再将各区域的土裂隙进行整合得到新的二值图像。所述的步骤2)中，正方形网格边长r的取值方式为从1像素值开始以1像素值递增直到图片尺寸大小为止，用以保证有足够多的数据。所述的步骤3)中，ln(N)-ln(1/r)一阶差商f1(r1,r2)的表达式为：ln(N)-ln(1/r)二阶差商f2(r1,r2,r3)的表达式为：其中，r1、r2、r3分别为正方形网格不同的边长，N1、N2分别为与边长r1、r2对应的裂隙穿过正方形的累计个数。选择n值的条件为：在局部区域尺寸不小于裂隙网络块区尺寸的前提下，保证每个区域都有裂隙存在，n值越大越好。图像分割后产生的全局分布的杂点，根据实际情况设置一个杂点的面积阈值，在MATLAB中采用bwareaopen函数将面积小于阈值的杂点去掉。在MATLAB中采用imclose函数对裂隙二值图进行桥接用以修复断裂和局部缺失的情况。与现有技术相比，本专利技术具有以下优点：一、本专利技术率先提出了采用差商分析的方法确定土裂隙无标度区，具有首创性。二、本专利技术结合土裂隙(ln(1/r1)，ln(N))散点一、二阶差商的分布特征，提出了差商分析中土裂隙无标度区的确定标准，为以后的无标度区的准确确定提供了可行的方法。三、相比于以往的研究，本专利技术在土裂隙网络分形维数的计算时充分考虑无标度区的影响，提升了分形维数计算结果的可靠性，对分形维数在土裂隙评价中的应用具有重要价值。附图说明图1为土裂隙提取过程图。图2为土裂隙ln(N)-ln(1/r)散点图。图3为土裂隙ln(N)-ln(1/r)一、二阶差商分布图，其中，图(3a)为土裂隙ln(N)-ln(1/r)一阶差商分布图，图(3b)为土裂隙ln(N)-ln(1/r)二阶差商分布图。图4为无标度区内分形维数计算过程。图5为本专利技术方法的流程示意图。具体实施方式下面结合附图和具体实施例对本专利技术进行详细说明。实施例本专利技术提供一种基于差商分析的土裂隙分形维数计算方法，共包括4个阶段：土裂隙提取、“数格子法”获取(ln(1/r)，ln(N))序列、基于差商分析法确定无标度区、计算土裂隙分形维数。第一阶段，土裂隙提取，对原始彩色裂隙图片进行灰度化处理，在MATLAB软件中采用rgb2gray函数对原始图像R、G、B分量进行加权平均，权重分别为0.2989、0.5870、0.1140，计算出像素点的灰度值，得到原始裂隙图片的灰度图；对灰度图进行图像分割，采用局部阈值分割的思想，先将图像分为n2个区域(n为整数，根据实际情况挑选合适的n值)，对每个区域分别用大津法(Ostu法)进行阈值分割，提取裂隙部分，再整合得到新的二值图像；在挑选n值时，在局部区域尺寸不小于裂隙网络块区尺寸的前提下，即保证每个区域都有裂隙存在，n越大越好；图像分割后产生的全局分布的杂点，根据实际情况设置一个杂点的面积阈值，在MATLAB中采用bwareaopen函数将面积小于阈值的杂点去掉；一些细小裂隙由于拍照过程采光不好或者分割、去杂操作的影响，可能出现断裂、局部缺失的情况，可采用闭运算的方法，在MATLAB中使用imclose函数对裂隙二值图进行桥接修复，土裂隙提取过程见附图1，左边图像为原始图像，右边图像为裂隙二值图。第二阶段，“数格子法”获取(ln(1/r1)，ln(N))序列。即采用若干边长(r1,r2,r3,…)的正方形网格去覆盖所研究的裂隙网络区域，相应地统计出裂隙穿过的正方形累计个数N1,N2,N3,…。正方形网格边长r序列取值方式为(1,2,3，…)，直到接近图片尺寸，以保证有尽可能多的数据以备分析。对一系列的(r,N)取自然对数，得到一个(ln(1/r1)，ln(N))序列，在坐标系下绘制ln(N)-ln(1/r)散点，如附图2所示。第三阶段，基于差商分析法确定无标度区，ln(N)-ln(1/r)差商定义如下：式中，f1(r1,r2)是ln(N)关于ln(1/r1)，ln(1/r2)的一阶差商，表示ln(N)在区间[ln(1/r1)，ln(1/r2)]上的平均变化率；f2(r1,r2,r3)是ln(N)关于ln(1/r1)，ln(1/r2)，ln(1/r3)的二阶差商，其值反映了一阶差商的平均变化率。图3表示裂隙图像在区间r∈(1,2,3,…,200)上ln(N)-ln(1/r)的一阶差商和二阶差商分布情况。结合土裂隙的ln(N)-ln(1/r)一、二阶差商的分布特征，找出同时满足ln(N)-ln(1/r)一阶差商值在1～2之间缓慢递减的区间(r∈[1,8])和二阶差商绝对值小于1的区间(r∈[1,12])，并将这两个区间的交集(r∈[1,8])作为裂隙网络的无标度区。第四阶段，在无标度区内，分形盒维数D的计算公式如式(6)所示：在实际应用中只能取有限的r。通常做法是在无标度区内，在ln(N)-ln(1/r)双对数平面内作图，基于最小二乘法用直线拟合ln(N)-ln(1/r)散点，得到的直线斜率即为土裂隙的分形盒维数D，计算得本文档来自技高网...

【技术保护点】
1.一种基于差商分析的土裂隙分形维数计算方法，其特征在于，包括以下步骤：1)对原始彩色裂隙图片进行灰度化和二值化处理，得到裂隙二值图，提取土体干缩裂隙；2)根据数格子法，采用多个不同边长r的正方形网格覆盖含有裂隙的网络区域，统计裂隙穿过的正方形的累计个数N；3)对多个边长‑个数对(r,N)取自然对数，得到(ln(1/r)，ln(N))序列，在坐标系下绘制ln(N)‑ln(1/r)的散点；4)基于差商分析法获取同时满足ln(N)‑ln(1/r)一阶差商值在1～2之间缓慢递减和二阶差商绝对值小于1两个条件的区间作为裂隙网络的无标度区；5)在无标度区内采用基于最小二乘法用直线拟合ln(N)‑ln(1/r)散点，得到的直线斜率即为土裂隙的分形盒维数。

【技术特征摘要】
1.一种基于差商分析的土裂隙分形维数计算方法，其特征在于，包括以下步骤：1)对原始彩色裂隙图片进行灰度化和二值化处理，得到裂隙二值图，提取土体干缩裂隙；2)根据数格子法，采用多个不同边长r的正方形网格覆盖含有裂隙的网络区域，统计裂隙穿过的正方形的累计个数N；3)对多个边长-个数对(r,N)取自然对数，得到(ln(1/r)，ln(N))序列，在坐标系下绘制ln(N)-ln(1/r)的散点；4)基于差商分析法获取同时满足ln(N)-ln(1/r)一阶差商值在1～2之间缓慢递减和二阶差商绝对值小于1两个条件的区间作为裂隙网络的无标度区；5)在无标度区内采用基于最小二乘法用直线拟合ln(N)-ln(1/r)散点，得到的直线斜率即为土裂隙的分形盒维数。2.根据权利要求1所述的一种基于差商分析的土裂隙分形维数计算方法，其特征在于，所述的步骤1)中，在MATLAB软件中采用rgb2gray函数获取像素点的灰度值，得到原始裂隙图片的灰度图。3.根据权利要求2所述的一种基于差商分析的土裂隙分形维数计算方法，其特征在于，所述的步骤1)中，采用局部阈值分割法进行裂隙图像二值化：先将灰度图分为n2个...

【专利技术属性】
技术研发人员：蔡奕，蔡万贤，吴乔枫，石涛，
申请(专利权)人：同济大学，
类型：发明
国别省市：上海,31