一种三维空间裂缝分离识别与表征方法

摘要

本发明公开一种三维空间裂缝分离识别与表征方法，对采集的图像进行如下处理，以实现对三维空间裂缝进行分离识别与表征：1)图像数据预处理；2)统计分析图像数据基本信息：图像数据的基本信息包括孔隙度、各个孔隙的连通性、孔隙大小统计、每个孔隙的位置、大小、方向以及各向异性；3)过滤：去除图像数据中非裂隙结构；4)平滑：对图像数据进行平滑和修补；5)减薄：使空隙结构在三维中最窄方向上减薄到厚度d；6)分离：以断开连接的方式分离裂缝网络中互相交错的裂缝；7)合并：合并在前一步中被断开的延伸较长的裂缝，整合分离过程形成的微小结构并恢复至减薄前状态，最后给出裂缝的表征。

权利要求

1.一种三维空间裂缝分离识别与表征方法，其特征在于，对采集的图像进行如下处理，以实现对三维空间裂缝进行分离识别与表征：

1)图像数据预处理；

2)统计分析图像数据基本信息：图像数据的基本信息包括孔隙度、各个空隙结构的连通性、位置、大小、方向以及各向异性；

3)过滤：去除图像数据中非裂隙结构；

4)平滑：对图像数据进行平滑和修补；

5)减薄：使空隙结构在三维中最窄方向上减薄到厚度d；

6)分离：以断开连接的方式分离裂缝网络中互相交错的裂缝；

7)合并：合并在前一步中被断开的空隙结构，整合分离过程形成的微小结构，最后分析各个空隙结构的孔隙度、连通性、位置、大小、方向以及各向异性，完成对三维空间裂缝的分离识别与表征；

所述步骤2)结构的各向异性的计算方式为：

对于一个包含了n个体像素的空隙结构，每一个体像素i都表达为从团簇中心点到当前位置的向量，空隙结构中心点坐标通过计算孔隙中所有点坐标的平均值算出；那么此空隙结构的整体各向异性能够用方向矩阵R来表示：该矩阵有三个特征值τ1<τ2<τ3，以及对应的特征向量μ1，μ2，μ3；各向同性指标I＝τ1/τ3和延伸指数E＝1-τ1/τ3用来定义结构的各向异性；当E→0且I→1时就表明该结构是各向同性的；

其中团簇是在二值数据中由具有共同面的相同材料体像素构成的独立结构，如果一个单独的体像素标识为1，没有与周围其它标示为1的体像素共面，则这一个体像素是一个团簇；多个相同标识的体像素能够通过共面方式连接成一个很大且很复杂的团簇，多条裂隙相互交错也构成一个团簇；

团簇是内部相互连通，但是与外界不连通的一个结构；

步骤3)所述过滤的过程为：当团簇满足过滤条件时，修改团簇内的所有体像素的材料标识，即将原空隙标识1修改为0，表示这一个小的空隙结构转变为岩石的一部分，在后续的分析中不存在，标识1表示为空隙结构，标识0表示为岩石的一部分；

步骤4)所述平滑的过程为：对过滤后的图像数据进行平滑和修补操作，具体是使用数学形态学中的闭运算实现图像的平滑和修补，形态学闭运算包括一个膨胀运算和一个紧接着的腐蚀运算；

膨胀运算是形态学中的基本运算，是一个求局部最大值的操作，首先用一个事先定义的结构元素B与图像A进行卷积，也即先计算结构元素B所覆盖的图像A的区域中的像素点灰度的最大值，并把这个最大值赋值给结构元素B的原点指定的像素点；这会使得灰度图像中高亮区域逐渐扩张，对于二值图像来说，也相当于将结构元素B与图像A进行或运算，运算结果决定原点指定的像素点的值；

其中结构元素是指一个原点位于中心的形状，该形状能为任意形状和大小；

腐蚀运算是形态学中的基本运算，与膨胀互为对偶，腐蚀本质上是一个求局部最小值的操作，操作方法与膨胀类似，最终用局部最小值赋值给结构元素B原点指定的像素点，此操作会使得灰度图像中的高亮区域逐渐缩减；对于二值图像来说，也相当于将结构元素B与图像A进行与运算，运算结果决定原点指定的像素点的值；由于腐蚀与膨胀互为对偶，故图像A被结构元素B腐蚀的补集等于A的补集被B膨胀；

形态学闭运算是对过滤后的图像先进行一次膨胀操作，再做一次腐蚀操作；

步骤5)所述减薄的过程为：判断当前点是否是需要处理的空隙点，对于空隙点，计算它在三个方向的延伸情况，即在x、y、z三个方向上逐点进行检查，如果在其中一个方向上遇到一个非空隙点，则停止该方向上的检查，并记录延伸长度；三个方向都检查完毕后，选取其中延伸长度最小的方向作为待减薄方向；之后将该方向减薄至3个体像素厚度，除此三个体像素，此次处理过程中该方向上的其他体像素全部改变赋值，即标识为0，使其成为岩石部分，并且记录相关信息；处理完成之后继续下一个点的计算；

步骤6)所述分离是对减薄后的数据进行分离，使得原本相互交错的结构分开，成为独立的结构，分离的具体过程为：(1)利用局部网格从三维上进行逐点处理；以当前处理像素为中心，分别以给定的一大一小两个分析半径，在三个方向上分别生成两个二维正方形局部网格，每次逐点处理只利用网格里的图像信息；将当前处理像素视为原点，两个局部网格都被原点划分为四个象限，对于外部半径大的局部网格来说，若存在任意两个象限内的图像是线性的，且它们之间的夹角处于所限定范围之内，同时内部半径较小网格内的相应象限中含有的空隙体像素数量不少于预设的数量，则它们被视为两条相交的方向不同的裂缝在当前截面上的投影，所以改变当前原点处像素的值，使其成为岩石部分；

(2)对于三维图像，处理步骤是将依次从三个方向上进行，直到确定当前像素是否应当被归为岩石。

2.根据权利要求1所述的三维空间裂缝分离识别与表征方法，其特征在于，步骤7)所述合并的过程为：首先判断每一个空隙结构是否应当与另外一个空隙结构合并，判断的条件包括：i)两个空隙结构的法线方向是否近似平行；ii)在分离步骤之前它们是否原本就属于同一空隙结构；iii)它们的法线是否都与它们的中心连线垂直；满足以上三个条件，则被判定为应当合并。

3.根据权利要求2所述的三维空间裂缝分离识别与表征方法，其特征在于，该方法还包括基于减薄步骤中记录的信息，将各个孔隙的形状恢复至减薄之前；最后分析各个空隙结构的孔隙度、连通性、位置、大小、方向以及各向异性，完成对三维空间裂缝的分离识别与表征。

说明书

一种三维空间裂缝分离识别与表征方法

技术领域

[0001] 本发明涉及裂缝识别领域，更具体地，涉及一种三维空间裂缝分离识别与表征方法。

背景技术

[0002] 岩石(和其它脆性材料)中发育裂缝十分常见。裂缝本身是材料受力破坏的产物，对其后续强度产生严重影响。同时，裂隙对流体运移(如地下水和油气)起着十分重要的作用。裂缝的表征(Characterization)或描述是所有与裂缝相关的力学研究的基础。裂缝的定量表征一般包括裂缝的长度、宽度、方向及分布密度等主要参数。对于一个由裂缝相互交错构成的裂缝网络，裂缝的表征首先需要将单条裂缝逐一从网络中分离出来。

[0003] 已有裂缝研究所采用的观测手段，大多以野外露头和岩心观测为主，也可以通过岩石薄片的显微镜观测获得。这些手段都只能提取观测面上裂缝分布的二维特征。基于二维图像或数据的裂缝分离识别，由于图像和数据的直观而相对很容易实现。如依据图1所示平面裂缝网络，人工进行裂缝分离识别显然是十分简易可行的。一个相互交错连接的二维裂缝网络，如图1(a)中各条裂缝可以被分离识别成单条裂缝的组合，如图1(b)，表示分离识别后的单条裂缝。

[0004] 当问题扩展到三维时，由于裂缝在三维空间中的分布和交互关系可能极其复杂，裂缝表征首先要实现的单条裂缝的逐一分离识别具有非常的技术难度。从图2所示的复杂三维裂缝网络可以看出，将其中各个单一裂缝分离识别出来将是巨大的挑战。

[0005] 三维裂缝的分离识别

[0006] 国内外都有学者对二维岩石裂缝的识别与表征进行过初步研究。Dare等(2002)曾提出了一种根据混凝土表面数字图像计算其中单条裂缝宽度的方法。Lee等(2013)提出了同样基于混凝土表面数字图像的方法，该方法可以自动提取裂缝的宽度、长度、方向等特征，并且使用了神经网络来识别裂缝模式。Jahanshahi和Masri(2013)也提出了利用图像细化和距离变换来完成裂缝图像特征提取的方法。Arena等(2014)通过识别交点实现了二维图像中多条相交裂缝的识别、分离与量化分析。国内学者中，杨松等人(2012)曾根据裂缝与其他元素在分形特征上的区别，提出过一种基于骨架和分形的混凝土裂缝图像识别算法。王平让等(2012)提出一种基于图像局部网格特征的裂缝识别技术，并在识别过程中自动计算裂缝的走向、长度和宽度等特征。

[0007] 需要强调指出的是，以上研究全部基于二维图像，三维的裂缝识别表征研究仅见一例。DellePiane等(2015)首先利用 软件的颗粒分离功能，将具有较完整颗粒形态的受热大理岩样品中的颗粒进行分离；继而将三个(或以上)颗粒的相交点确定为裂缝相交点，以这些相交点分割构成网络的裂缝。裂缝分离的这一尝试性成功，十分依赖于裂缝的形态。注意到该方法中裂缝的分割点是由三个(或以上)颗粒的相交点决定的，因此，该方法仅适用于构成封闭面的裂缝系统。对于一般情况的裂缝网络(如一种极简单情形：两条相交而没有达到模型边界的裂缝)并不适用。另外该方法需借助 软件的颗粒分割功能进行操作，一方面， 软件的颗粒分离获得的分离边界构成的裂缝网络可能与原有图像中可识别的裂缝网络具有较大差异，同时 软件所能处理的体积大小较有限(Liu et al.2016)。

[0008] 裂缝的表征

[0009] 微观层析成像方法发展以来，相关数据分析更多地关注孔隙介质，并在此基础上较好地发展了孔隙结构的精细表征技术和方法(Ketcham&Iturrino,2005；Nakashima&Kamiya,2007,Liu et al.,2009,2014)。

[0010] 裂缝结构的精细表征由于裂缝分离识别技术没有实现而难以发展。在裂缝分离识别可实现的情况下，获得裂缝的长度、宽度、方向、密度等主要参数均不存在任何难度。DellePiane等(2015)在实现其样品中特定结构裂缝的分离之后，利用 软件的统计分析功能，给出了相应的裂缝表征和统计结果。

[0011] 在裂缝分离识别实现后，每一条裂缝被定义为一个“团簇”(可以理解为一个独立的结构)。这时，微观层析成像数据的结构表征分析程序ctsta10.f90可用于裂缝表征。该程序原来主要用于孔隙结构的精细表征，输出参数包括：孔隙度、比表面积、各单个团簇(结构)的大小、位置、形态、方向等，也可以统计颗粒和孔径的大小分布。

[0012] 鉴于此，实现对三维空间裂缝的表征，关键技术在于裂缝的分离识别。

发明内容

[0013] 为了克服上述现有技术的不足，本发明提出一种三维空间裂缝分离识别与表征方法。

[0014] 为了解决上述技术问题，本发明的技术方案为：

[0015] 一种三维空间裂缝分离识别与表征方法，对采集的图像进行如下处理，以对三维空间裂缝进行分离识别与表征：

[0016] 1)图像数据预处理；

[0017] 2)统计分析图像数据基本信息：图像数据的基本信息包括孔隙度、各个空隙结构的连通性、位置、大小、方向以及各向异性；

[0018] 3)过滤：去除图像数据中非裂隙结构；

[0019] 4)平滑：对图像数据进行平滑和修补；

[0020] 5)减薄：使空隙结构在三维中最窄方向上减薄到厚度d；

[0021] 6)分离：以断开连接的方式分离裂缝网络中互相交错的裂缝；

[0022] 7)合并：合并在前一步中被断开的空隙结构，整合分离过程形成的微小结构，最后分析各个空隙结构的孔隙度、连通性、位置、大小、方向以及各向异性，完成对三维空间裂缝的分离识别与表征。

[0023] 优选的，所述步骤2)孔隙的各向异性的计算方式为：

[0024] 对于一个包含了n个体像素的空隙结构，每一个体像素i都表达为从团簇中心点到当前位置的向量ai＝(axi,ayi,azi)T，空隙中心点坐标通过计算孔隙中所有点坐标的平均值算出；那么此空隙结构的整体各向异性能够用方向矩阵R表示：

[0025]

[0026] 该矩阵有三个特征值τ1＜τ2＜τ3，以及对应的特征向量μ1，μ2，μ3；各向同性指标I＝τ1/τ3和延伸指数E＝1-τ1/τ3用来定义空隙结构的各向异性；当E→0且I→1时就表明该空隙结构是各向同性的；

[0027] 其中团簇是在二值数据中由具有共同面的相同材料体像素构成的独立结构，如果一个单独的体像素标识为1，没有与周围其它标示为1的体像素共面，则这一个体像素是一个团簇；多个相同标识的体像素能够通过共面方式连接成一个很大且很复杂的团簇，多条裂隙相互交错也构成一个团簇；

[0028] 团簇是内部相互连通，但是与外界不连通的一个结构。

[0029] 优选的，步骤3)所述过滤的过程为：当团簇满足过滤条件时，修改团簇内的所有体像素的材料标识，即将原孔隙标识1修改为0，表示这一个小的空隙结构转变为岩石的一部分，在后续的分析中不存在，标识1表示为空隙结构，标识0表示为岩石的一部分。

[0030] 优选的，步骤4)所述平滑的过程为：对过滤后的图像数据进行平滑和修补操作，具体是使用数学形态学中的闭运算实现图像的平滑和修补，形态学闭运算包括一个膨胀运算和一个紧接着的腐蚀运算；

[0031] 膨胀运算是形态学中的基本运算，是一个求局部最大值的操作，首先用一个事先定义的结构元素B与图像A进行卷积，也即先计算结构元素B所覆盖的图像A的区域中的像素点灰度的最大值，并把这个最大值赋值给结构元素B的原点指定的像素点；这会使得灰度图像中高亮区域逐渐扩张，对于二值图像来说，也相当于将结构元素B与图像A进行或运算，运算结果决定原点指定的像素点的值；

[0032] 其中结构元素是指一个原点位于中心的形状，该形状能为任意形状和大小；

[0033] 腐蚀运算是形态学中的基本运算，与膨胀互为对偶，腐蚀本质上是一个求局部最小值的操作，操作方法与膨胀类似，最终用局部最小值赋值给结构元素B原点指定的像素点，此操作会使得灰度图像中的高亮区域逐渐缩减；对于二值图像来说，也相当于将结构元素B与图像A进行与运算，运算结果决定原点指定的像素点的值；由于腐蚀与膨胀互为对偶，故图像A被结构元素B腐蚀的补集等于A的补集被B膨胀；

[0034] 形态学闭运算是对过滤后的图像先进行一次膨胀操作，再做一次腐蚀操作。

[0035] 优选的，步骤5)所述减薄的过程为：判断当前点是否是需要处理的空隙点，对于空隙点，计算它在三个方向的延伸情况，即在x、y、z三个方向上逐点进行检查，如果在其中一个方向上遇到一个非空隙点，则停止该方向上的检查，并记录延伸长度；三个方向都检查完毕后，选取其中延伸长度最小的方向作为待减薄方向；之后将该方向减薄至3个体像素厚度，除此三个体像素，此次处理过程中该方向上的其他体像素全部改变赋值，即标识为0，使其成为岩石部分，并且记录相关信息；处理完成之后继续下一个点的计算。

[0036] 优选的，步骤6)所述分离是对减薄后的数据进行分离，使得原本相互交错的结构分开，成为独立的结构，分离的具体过程为：

[0037] (1)利用局部网格从三维上进行逐点处理；以当前处理像素为中心，分别以给定的一大一小两个分析半径，在三个方向上分别生成两个二维正方形局部网格，每次逐点处理只利用网格里的图像信息；将当前处理像素视为原点，两个局部网格都被原点划分为四个象限，对于外部半径大的局部网格来说，若存在任意两个象限内的图像是线性的，且它们之间的夹角处于所限定范围之内，同时内部半径较小网格内的相应象限中含有的空隙体像素数量不少于预设的数量(本发明将这个数量设置为一个可调节的输入参数，以便于不同裂缝条件的具体操作)，则它们被视为两条相交的方向不同的裂缝在当前截面上的投影，所以改变当前原点处像素的值，使其成为岩石部分；

[0038] (2)对于三维图像，处理步骤是将依次从三个方向上进行，直到确定当前像素是否应当被归为岩石。

[0039] 优选的，步骤7)所述合并的过程为：首先判断每一个空隙结构是否应当与另外一个空隙结构合并，判断的条件包括：i)两个空隙结构的法线方向是否近似平行；ii)在分离步骤之前它们是否原本就属于同一空隙；iii)它们的法线是否都与它们的中心连线垂直；满足以上三个条件，则被判定为应当合并。合并过程就是将分离步骤之后分别属于不同团簇标号的空隙结构在判定为属于同一条裂缝的前提下赋予相同的团簇标号。

[0040] 优选的，该方法还包括基于减薄步骤中记录的信息，将各个空隙的形状恢复至减薄之前；即得到最终的处理结果。最终结果中的每一个团簇代表一条裂缝，并且是恢复至减薄之前的状态，每一个团簇的位置、大小、表面积、方向等参数都可以准确计算和表述，从而实现裂缝的精确表征。

[0041] 与现有技术相比，本发明的有益效果为：本发明从三维复杂裂缝网络中逐一识别单条裂缝，进而实现对裂缝的表征，该方法为开展裂缝形态、渗透率及力学响应、尺度放大等研究提供重要技术支持。

附图说明

[0042] 图1为相互交错连接的二维裂缝网络示意图，(a)为识别前的示意图，(b)为各条裂缝被分离识别成单条裂缝的组合示意图。

[0043] 图2为岩石样品三维裂缝结构示例图。

[0044] 图3为裂缝分离识别与表征方法流程图。

[0045] 图4为从二维切片看二值化效果示意图。

[0046] 图5为岩石过滤前后对比图。

[0047] 图6为膨胀算法实现过程图。

[0048] 图7为腐蚀算法实现过程图。

[0049] 图8为闭运算实现过程图。

[0050] 图9为平滑效果示意图。

[0051] 图10为减薄前后对比图。

[0052] 图11为分离过程示意图，(a)断开前状态；(b)为裂缝交叉点集被断开。

[0053] 图12为分离前后对比图。

[0054] 图13分离后的各个结构图。

[0055] 图14为T型相交示意图。

[0056] 图15为拐弯型相交示意图。

[0057] 图16为十字型相交示意图。

[0058] 图17为合并过程图。

[0059] 图18为典型情况1的示意图。

[0060] 图19为典型情况2的示意图。

[0061] 图20为合并前后对比图(二维平面显示)。

[0062] 图21为合并之后的各个结构的依次示意图。

具体实施方式

[0063] 下面结合附图对本发明做进一步的描述，但本发明的实施方式并不限于此。

[0064] 一般采用微观CT图像需要进行一些常规预处理，包括平滑、去噪、边缘切割、图像分割(识别裂缝或孔隙)。因此实际工作需要7个步骤，图像分割之后的主要步骤为6个，如图3所示，即：

[0065] 0)图像数据预处理

[0066] 1)基本信息统计分析：图像数据的基本信息，如孔隙度、连通性、结构各向异性等[0067] 2)过滤：去除数据中非裂隙结构，如很小的或各向同性强的结构

[0068] 3)平滑：对图像数据进行平滑和修补

[0069] 4)减薄：使空隙结构在三维中最窄方向上减薄到一定厚度

[0070] 5)分离：以断开连接的方式分离裂缝网络中互相交错的裂缝

[0071] 6)合并：合并在前一步中被强行断开的空隙，整合分离过程形成的微小结构，并将各个空隙的形状恢复至减薄之前，最后分析各个空隙结构的孔隙度、连通性、位置、大小、方向以及各向异性，完成对三维空间裂缝的分离识别与表征。

[0072] 下面先简略介绍数据预处理，再分别介绍6个主要实施步骤。

[0073] 微观CT图像的预处理可以由多种图像处理软件实现，其中 软件最为通用。首先，一组岩石样品的三维CT数据被载入到图像软件中，原始图像一般为灰度图像，在图像质量不够好的情况下，可能需要进行图像的平滑、去噪。其次，需要从原始数据中切割出需要的部分，或者说，去除图像中非样品部分。图像切割一般是切割为平行六面体，也可以切割为圆柱体。第三，需要将切割后的图像数据转化为二值图像，这一步称为图像分割，也称为二值化，是将感兴趣的结构或材料识别出来，其余部分理解为基质。分割之后的图像只有两个值，在计算机中用0或1代替，图像显示则转变为黑白图像，如图4。二值化之后的三维图像可参见图2。图像分割的方法有多种，这里不做介绍。

[0074] 基本信息统计分析

[0075] 这一步目的在于获得图像数据的主要基本信息，包括：样品孔隙度、各个团簇的连通性、位置、大小、方向以及各向异性。这些信息可用于确定后续分析步骤中需要用到的各个参数的取值。实施程序为ctsta10.f90，为专利发明人刘洁所编写。

[0076] 在基本信息分析处理中，涉及到“团簇(cluster)”这个概念。一个团簇是在二值数据中由具有共同面的相同材料(如孔隙，材料标识为1)体像素构成的独立结构。如果一个单独的体像素标识为1(表示孔隙)，没有与周围其它标示为1的体像素共面，则这一个体像素是一个团簇。多个相同标识的体像素可以通过共面方式连接成一个很大且很复杂的团簇，如砂岩中的孔隙网络。多条裂隙相互交错也构成一个团簇。团簇也可以理解为内部相互连通，但是与外界不连通的一个结构。

[0077] 众多参数中，孔隙的各向异性为重点，通过以下方法计算各向异性：

[0078] 对于一个包含了n个体像素的团簇，每一个体像素i都可以表达为从团簇中心点到当前位置的向量ai＝(axi,ayi,azi)T，团簇中心点坐标可通过计算团簇中所有点坐标的平均值算出。那么，此团簇结构的整体各向异性可以用方向矩阵R来表示：

[0079]

[0080] 该矩阵有三个特征值τ1＜τ2＜τ3，以及对应的特征向量μ1，μ2，μ3；各向同性指标I＝τ1/τ3和延伸指数E＝1-τ1/τ3用来定义空隙结构的各向异性；当E→0且I→1时就表明该空隙结构是各向同性的。

[0081] 基本信息统计分析有7个输出文件，每个文件记录了图像数据的不同信息。

[0082] 1)基本信息：给出了图像数据的基本信息，包括孔隙度、比表面积、团簇数目(此处团簇即样品内部的空隙结构)、连通性以及其他统计信息。

[0083] 2)颗粒/孔喉分布。

[0084] 3)连通性结果：简单列出了当前图像数据中所有的连通团簇信息，包括它们的编号以及连通状态。在x，y，z三个方向上，0代表不连通，1代表连通。

[0085] 4)团簇大小统计：给出了包括团簇编号、包含体像素数目、等效半径平方、团簇累计包含体像素数目以及占总空隙体像素的累计百分比。

[0086] 5)团簇形状：给出了包括团簇大小、表面积、主方向延伸长度、各向同性指数、延伸指数的信息。

[0087] 6)团簇方向矩阵：给出了每个团簇的中心坐标以及方向矩阵。

[0088] 7)归一化的团簇方向矩阵：给出了每个团簇的归一化方向矩阵以及每个特征向量的方向角。

[0089] 过滤

[0090] 这一步的目的在于去除模型中不重要的团簇，使得后续计算分析更方便。可以选择的过滤参数包括团簇的大小、位置、是否连通、各向异性等。

[0091] 例如，一般情况下，自然样品中会含有很多细小的结构，代表样品中存在的小孔隙和图像数据获取时产生的噪声，这些小孔隙在裂缝分析中可以忽略，因此将这些过小的结构过滤掉有助于后续分析。

[0092] 基本信息分析处理的结果中得到的结构参数为确定过滤阀值提供了依据。

[0093] 当某一个团簇(结构)满足过滤条件时，该团簇内的所有体像素的材料标识将被修改，即由原来的空隙标识1修改为0，表示这一个小的空隙结构转变为岩石的一部分，在后续的分析中不存在。

[0094] 除了岩石中包含的微小孔隙的情形以外，也存在空隙结构中包含标识为0的固体岩石的情形。对此可以通过调换数据文件中标号的方式(如labmat＝1-labmat)，生成一个相反标识的文件，再次运行过滤。处理完成之后需要将数据再次反转并存储。图5显示过滤后一些过小的空隙和图像噪声得到了显著压制。

[0095] 平滑

[0096] 过滤后的数据作为这一步的输入数据。由于自然样品结构的复杂性，仅经过过滤的图像数据仍然难以处理，所以在这一步中，对图像数据进行了平滑与修补操作，使图像形状尽量简洁，便于处理。使用数学形态学中的闭运算实现图像的平滑和修补。形态学闭运算包括一个膨胀运算和一个紧接着的腐蚀运算。

[0097] 结构元素：可以为任意形状和大小，并拥有一个参考点，称之为原点。一般情况下，结构元素都指一个原点位于中心的正方形或圆盘，但也可以使用其他形状。

[0098] 膨胀运算：形态学中的基本运算。本质上是一个求局部最大值的操作，我们用一个事先定义的结构元素B与图像A进行卷积，也即先计算结构元素B所覆盖的A的区域中的像素点灰度的最大值，并把这个最大值赋值给结构元素B的原点指定的像素点。这会使得灰度图像中高亮区域逐渐扩张。对于二值图像来说，也相当于将结构元素B与图像A进行或运算，运算结果决定原点指定的像素点的值(图6)。

[0099] 腐蚀运算：形态学中的基本运算。与膨胀互为对偶，腐蚀本质上是一个求局部最小值的操作，操作方法与膨胀类似，最终用局部最小值赋值给结构元素B原点指定的像素点。此操作会使得灰度图像中的高亮区域逐渐缩减。对于二值图像来说，也相当于将结构元素B与图像A进行与运算，运算结果决定原点指定的像素点的值(图7)。由于腐蚀与膨胀互为对偶，故图像A被结构元素B腐蚀的补集等于A的补集被B膨胀。

[0100] 闭运算：对图像先进行一次膨胀操作，再做一次腐蚀操作。一般来说，闭运算可以填平小坑洞，弥合小裂缝，图像总体的位置和形状保持不变(图8)。

[0101] 在三维数据的处理中，结构元素相应地由正方形和圆盘改变为立方体和球体。不同形状和大小的结构元素会产生不同的处理结果，结构元素应当比所有要去除的噪声块的尺寸都要大。也可以适当增加膨胀和腐蚀的次数，如多次膨胀之后多次腐蚀，来使平滑效果更好。但应注意，过大的结构元素和过多次膨胀之后腐蚀有可能会去除原来图像数据的某些特征，所以应当多次尝试，根据处理效果选取适当的结构元素种类、大小以及操作次数(图9)。

[0102] 减薄

[0103] 对于裂缝分析，需要对图像采用减薄算法进行处理，得到与原来物体区域形状近似的、由较为简单的曲面组成的图形。

[0104] 本实施例提出的减薄算法与二维图像处理中的形态学细化算法有些类似，但并不完全相同。应该注意到，本发明处理的是三维数据，如果仅从二维的角度去处理三维图像，需要将三维图像看作是一组二维图像的集合，并对每张二维图像采用二维中方法。而由于三维图像在二维截面中的表现并不稳定，即使是相邻的截面，也可能存在差异较大的情况，这样如果对两个二维截面单独使用细化算法，处理后的两个截面中的图像会错开，以至于使得三维中原来的一个结构被误分为两个独立的结构。另外，仅从二维角度去处理三维图像，就意味着需要同时对三个方向上的截面都进行处理，而三个方向的截面处理后的结果势必存在着不一致性，如何综合这三个方向的结果也是一个非常棘手的问题。故二维图像中常用的图像细化方法在此处并不适用。

[0105] 本实施例提出的减薄算法则是从三维的角度进行计算，首先判断当前点是否是需要处理的空隙点，对于空隙点，计算它在三个方向的延伸情况，即在x、y、z三个方向上逐点进行检查，如果在其中一个方向上遇到一个非空隙点，则停止该方向上的检查，并记录延伸长度。三个方向都检查完毕后，选取其中延伸长度最小的方向作为待减薄方向。之后将该方向减薄至3个体像素厚度，除此三个体像素，此次处理过程中该方向上的其他体像素全部改变赋值，使其成为岩石部分，并且记录相关信息(图10)。处理完成之后继续下一个点的计算。需要注意的是，减薄是便于处理图像数据的一个手段，并非所需的图像的形态，故在后面的步骤中，还会利用这里记录的信息来恢复被减薄的孔隙。

[0106] 分离

[0107] 这是本发明的关键所在。这一步以上一步经过减薄的数据为输入，根据一系列人工设定的条件，对图像中的空隙结构进行分离处理，使得原本互相交错的结构分开，成为独立的结构。

[0108] 在这一步中，首先给各个空隙结构按照包含体像素的多少编号，最大的结构为1号。然后单独计算每个结构在三维空间中占据的空间尺寸，即x、y、z三个方向上该结构的最小坐标与最大坐标。为避免受到其他结构影响，接下来将单独对每个孔隙结构进行处理，而处理的空间也就是三维方向上最小坐标与最大坐标构成的长方体空间，这样在处理过程中程序的搜索空间小于原数据的空间，会大大节省程序的运行时间。

[0109] 具体的处理步骤分为两部分。

[0110] 第一部分利用局部网格从三维上进行逐点处理。以当前处理像素为中心，分别以给定的一大一小两个分析半径，生成两个二维正方形局部网格，每次逐点处理只利用网格里的图像信息。可以将当前处理像素视为原点，两个局部网格都被原点划分为四个象限，对于外部半径较大的局部网格来说，若存在任意两个象限内的图像是线性的，且它们之间的夹角处于限定范围之内，同时内部网格的相应象限中含有的空隙体像素数量不能过少，则它们被视为两条相交的方向不同的裂缝在当前截面上的投影，所以改变当前原点处像素的值，使其成为岩石部分(图11)。使用内外两个半径作为判断条件的原因在于，小半径因所含体像素数量较少，用它计算得出的图像方向会误差较大，因此我们利用了较大半径的局部网格判断图像的方向，同时利用较小半径的网格来避免大半径导致的去除过多的现象。

[0111] 由于处理的图像是三维的，故以上处理步骤将依次从三个方向上进行，直到确定当前像素是否应当被归为岩石。经过这一步，大部分交错的空隙结构都能被断开成两个或多个独立的空隙结构。断开的效果与参数设置有关，所以应尝试多次并根据实验结果适当选取，例如外部半径过大时，仍然会导致断开处宽度过大，对图像造成破坏；太小的半径所包含的空隙体像素个数过少，计算出的图像各向异性非常不准确。夹角的范围同样需要人工设置，在本实施例两个结构夹角在0°～30°以及150°～180°时可以视为近似平行，而夹角在30°～150°之间的结构则视为互相交错的不同裂缝结构。

[0112] 在第一部分中，如果当前局部网格包含的空隙体像素个数过少，就可能造成图像并非是线性结构，从而使得这些点不会在第一部分被去除。因此存在这样一种情况：一些较大的空隙仅通过几个体像素连接，成为了一个团簇。这种情况在第一部分就不能得到很好的处理，此时需要继续进行第二部分处理。

[0113] 第二部分同样是从三个方向上进行处理，但不再为逐点处理，而是逐层处理。将当前分析的团簇所占空间划分为若干个不重叠的正方体空间，然后对每一个正方体空间进行逐层分析。在小正方体内，首先沿Z轴方向分析三维局部网格的第一层图像的各向异性，并判断图像是否为线性，随后计算下一层图像的各向异性并判断是否为线性。如果相邻两层的图像都为线性，且各向异性相差较大，即相邻层内结构非近似平行，则它们被视为不同的孔隙结构，第二层图像中所有原属于孔隙的点被标识为岩石类别，使得不同的孔隙结构不再接触。接着程序继续使用相同方法处理下一个相邻的两层图像，直到计算至最后一层。再分别从Y轴方向和X轴方向重复此操作。全部完成后再进入下一个三维局部网格的计算。

[0114] 经历这两部分分离操作后的结果如果课见于图12以及图13。

[0115] 这一步骤中用来计算局部网格内图像方向的算法与三维方向矩阵R的计算(见基本信息统计分析)原理相同，在二维情况下相应地有：

[0116]

[0117] 该矩阵有两个特征值τ1＜τ2，以及对应的特征向量μ1，μ2。各向同性指标I＝τ1/τ2。使用R的值来判断该空隙在二维上的投影是否为线性结构，利用Fortran数学库函数中ssyev函数可以计算得到该矩阵的两个特征向量，两个特征向量可以分别代表该象限中图像的短轴和长轴与x轴和y轴夹角的余弦值，即可以判断图像方向。

[0118] 以下为实际处理的图像数据中常出现的三种典型情况：

[0119] 1)如图14，裂缝为T型相交的两条时，以中心蓝色部分为原点，整个团簇被视为由三个部分组成：第一象限的A，第二象限的B以及第四象限的C。程序将分别计算A、B、C的各向异性(或方向)，在原点周围的空隙体像素数目足够多时，会继续计算A、B、C之间的夹角，若其中存在夹角的余弦值的绝对值小于给定阈值，则意味着A、B、C中存在着夹角较大的相交情况，需要将原点处的空隙体像素改为岩石体像素。

[0120] 2)裂缝为图15所示时，裂缝出现大角度转折，也可以认为是两条裂缝相连，需要断开。以中心蓝色部分为原点，整个团簇被视为由两个部分组成：第一象限的A和第四象限的B。程序将分别计算A、B的各向异性(或方向)，在原点周围的空隙体像素数目足够多时，会继续计算A、B之间的夹角，若A、B存在夹角的余弦值的绝对值小于给定阈值，则意味着A、B属于夹角较大的相交情况，需要将原点处的空隙体像素改为岩石体像素。

[0121] 3)裂缝为图16所示的十字型相交时，以中心蓝色部分为原点，整个团簇被视为由四个部分组成：第一象限的A、第二象限的B、第三象限的C和第四象限的D。程序将分别计算A、B、C、D的各向异性(或方向)，在原点周围的空隙体像素数目足够多时，会继续计算A、B、C、D之间的夹角，若存在任意两者之间夹角的余弦值的绝对值小于给定阈值，则意味着它们属于夹角较大的相交情况，需要将原点处的空隙体像素改为岩石体像素。

[0122] 合并

[0123] 这一步以上一步的结果为输入数据，是之前步骤的修正，即对于实际上连通的裂缝，在分离过程中被段开，通过这一步操作使之归并为一个团簇，也包括对一些分离过程形成的小团簇进行处理，并且将裂缝恢复至减薄之前的状态。

[0124] 首先，对图像数据进行标号，因为经过了分离步骤，所以这次的图像中包含的团簇数会比之前更多。完成之后，开始一一判断每一个空隙结构是否应当与另外一个空隙结构合并，判断的条件包括：i)两个空隙结构的法线方向是否近似平行；ii)在分离步骤之前它们是否原本就属于同一空隙；iii)它们的法线是否都与它们的中心连线垂直。满足以上三个条件，才能被判定为应当合并。

[0125] 合并操作包括改写供给方团簇编号为接收方团簇编号，改写两个团簇含有的体像素总数，以及记录接收方团簇的编号，这样供给方团簇就被并入了接收方团簇。程序在运行时很可能会遇到这样的情况：作为接收方的团簇之前也作为供给方，并且已经被并入其他团簇。此时需要用到之前合并时记录的接收方编号，找到接收方团簇，并将此孔隙作为新的接收方接受并入(图17)。

[0126] 再如图18所示，由于分离算法，原本属于同一个团簇的裂缝现在被分为了三条，中心部分已经更改为岩石，即黑色的体像素。对于这三条独立的裂缝，它们的法线分别为L1、L2和L3，L4为其中两条裂缝的中心点连线。由于法线L1并不与L2和L3近似平行，所以它所在的裂缝首先被判定不与其他两条裂缝合并。而L2与L3近似平行，则继续判断这两个裂缝的中心点连线是否与它们的法线垂直。由图中可以看出，L4与L1和L2都近似垂直，加上这两条裂缝原本属于同一个团簇，则它们被判定可以合并。

[0127] 而如图19的情况中，两条裂缝的法线L1、L2平行，但它们的中心连线L3与L1、L2并不能构成近似垂直，所以这张图中的两条裂缝不能合并。

[0128] 由于分离过程中可能产生一些较小的碎片式的空隙(如图15最后一幅图)，这些空隙因为尺寸较小，往往不具有明确的取向，也就无法计算其法线方向，因此在上面的合并操作中很难得到妥善处理。于是在这里，采用计算类间距离的方法将他们并入较大的空隙结构。取这些小空隙的各自的中心点代表它们的坐标，依次计算它们到其他较大孔隙的体像素的最小距离，选择距离最近的较大空隙，将小空隙并入。

[0129] 最后，利用减薄步骤中记录的信息，将各个孔隙的形状恢复至减薄之前。即得到最终的处理结果(图20、21)。

[0130] 裂缝表征

[0131] 完成全部处理之后，再对分离后的图像数据进行分析表征。这里仅以图13所示的实验数据为例，将一些主要结果以列表形式给出。表1显示，经过处理的数据中含有的空隙体像素数目减少了，这是因为裂缝分离时将一些体像素的标识值修改了，表面积的增加也同样是因为相交裂缝断开造成。分离前仅有3个团簇，分离后为8个。团簇的大小和形态参数变化十分明显。表3数据可以给出模型中裂缝的长度、宽度的统计规律。更多信息可以从ctsta10.f90的7个输出文件中提取。

[0132] 表1裂缝分离识别前后主要参数对比

[0133]参数 分离前 分离后
孔隙度 3.18 2.847％
表面积 177856 177346

[0134] 表2裂缝分离识别前后团簇大小对比

[0135]

[0136] 表3裂缝分离识别前后团簇形态对比

[0137]

[0138] 以体像素为基本单位的二值图像数据，通过基本特征分析、过滤、平滑、减薄、分离、合并6个步骤，可以完成裂缝的分离识别和表征。这一技术将为开展裂缝形态、渗透率及力学响应、尺度放大等研究提供重要技术支持。

[0139] 在上述的文件中，空隙主要是强调图像中表示出来的“空”隙而非岩石固体。在本发明申请中，主要针对空隙呈裂缝状态显现的情况，也不排除以小孔隙状态显现的情况。典型的孔隙型结构不在本发明申请的保护范围内。

[0140] 以上所述的本发明的实施方式，并不构成对本发明保护范围的限定。任何在本发明的精神原则之内所作出的修改、等同替换和改进等，均应包含在本发明的权利要求保护范围之内。