岩石三维图像裂缝提取方法与流程

本发明涉及一种图像分割与识别技术，尤其涉及一种岩石三维图像裂缝提取方法，属于图像分割及识别技术领域。

技术背景

岩石在应力作用下产生机械性破坏，无明显位移的断裂构造叫裂缝。裂缝是一种特殊的孔隙。在石油地质领域，裂缝识别是油藏和储存非均质研究的重要内容，对油气开发和原油采收率的提高具有重要意义。

本发明针对的是岩石CT序列图像的裂缝识别。本发明需先对岩石CT序列图中的每幅图像进行二值化，再利用二值化的序列图构建出三维岩石孔隙模型，并在三维空间识别出岩石裂缝。经由CT获取的岩石图像是一组连续的二维序列图，利用传统二维图像分割算法，可依次对序列图中每一幅CT图像作裂缝提取，进而便可得到裂缝的三维图像。但上述方法仅利用了二维图像信息，忽略了岩石裂缝在多张连续岩石CT图像之间的关联信息。这将导致CT图像呈狭长线条状而三维空间呈非裂缝形态的普通孔隙易被误判为裂缝，或CT图像呈非狭长线条状而三维空间呈裂缝形态的裂缝易被误判为非裂缝。

目前已有的基于三维空间的裂缝识别方法是依据每个孤立三维孔隙的最小外接球半径与等效球半径的比值、形状因子、近似最小外接长方体最长边与最短边的比值等参数来识别岩石裂缝。但由于裂缝通常不是独立存在的，裂缝与普通孔隙之间关联紧密，故此方法对裂缝的误识率较高，如将与普通孔隙相连的裂缝判为非裂缝，或将与裂缝相连的普通孔隙判为裂缝。

近年来，研究人员在不断地研究与探索三维岩石图像裂缝识别问题，但目前还未提出一种可对三维岩石图像裂缝进行高效准确识别的方法。本发明正是基于这种研究现状，提出一种可高效准确识别三维岩石图像裂缝的方法。

技术实现要素：

本发明的目的在于克服现有技术存在的缺陷和不足，提供一种三维岩石图像裂缝的识别方法。

本发明提供的三维岩石图像裂缝识别方法基于下述思想：裂缝是一种特殊的孔隙，常与普通孔隙相连接，故不能直接判定某个孔隙是否为裂缝。本发明对三维岩石孔隙模型的每个连通分量执行表面重建、拉普拉斯网格平滑、网格简化、网格单元单位法向量计算等一系列操作。由于裂缝存在面状特征，位于裂缝表面的三角网格单元被简化的概率高。三角网格单元被简化后网格三角面积将变大，且网格单位法向量方向趋于统一。根据网格三角面积和网格单位法向量方向特征，对三角网格单元执行K-均值聚类，根据聚类后得到的聚类中心将三角网格类重新划分，位于同一裂缝表面的三角网格单元将聚于同一类中。进而通过形状因子等参数即可判定每个三角网格类构成的三维空间结构是否具有裂缝特征。为了提高裂缝识别的精确度，对具有裂缝特征的三维空间结构执行三维图像形态学膨胀操作，之后将膨胀得到的三维空间结构的体素点集与原始三维岩石孔隙模型连通分量的体素点集进行逻辑与操作，与操作的结果即为岩石裂缝。

本发明提供的三维岩石图像裂缝的识别方法，包括以下操作步骤：

预处理：对一组待识别的岩石二值序列图构建其相应的三维岩石孔隙模型，并提取模型所有连通分量；剔除其中体素点过少的连通分量；使用移动立方体等值面提取(Marching Cube)算法对每个连通分量进行表面重建，得到每个连通分量相应的三角网格模型A_q，其中q＝1，2，3…k，k为三角网格模型个数；对得到的三角网格模型作拉普拉斯网格平滑处理，并使用网格简化算法对三角网格模型进行网格简化；接着对每一个A_q独立执行下述步骤1到步骤10的操作；

步骤1：计算预处理得到的三角网格模型A_q中每个三角网格单元的单位法向量与三角面积，并调整单位法向量方向；接着利用三角网格单元的单位法向量及三角面积特征对A_q的三角网格单元进行K-均值聚类；

步骤2：根据步骤1得到的聚类中心将A_q的三角网格单元重新分类，分类得到多个三角网格类；

步骤3：提取出步骤2中每个三角网格类的三角网格连通分量；

步骤4：剔除步骤3得到的三角网格连通分量中的噪声成分，重新计算三角网格类的聚类中心，并将重新计算得到的聚类中心相邻的所有三角网格类合并为一类；接着对每个三角网格类分别独立执行步骤5到步骤10的操作；

步骤5：步骤4得到的每个三角网格类中的所有三角网格单元组成一个三维面结构，对这个三维面结构执行三维模型孔洞修补操作，操作后可得到一个或多个封闭的三维面结构，这些三维面结构组成一个封闭面结构集合；

步骤6：计算步骤5获得的封闭面结构集合的最小外接长方体；

步骤7：计算步骤5获得的封闭面结构集合的形状因子；计算步骤6得到的最小外接长方体的形状因子；计算面积比参数；

步骤8：利用步骤7得到的两个形状因子及面积比等参数，判定当前封闭面结构集合是否为裂缝。若判定结果为真，则继续执行后续操作；若判定结果为假，则当前封闭面结构集合不是裂缝；

步骤9：提取步骤8中判定结果为真的封闭面结构集合表面及内部的体素点，对提取的体素点集作三维图像形态学膨胀处理；

步骤10：将步骤9中已膨胀的体素点集与原始三维岩石孔隙模型连通分量的体素点集进行逻辑与操作，操作后得到的体素点集即为裂缝。

上述技术方案中，预处理中所述剔除体素点过少的连通分量，优先采用的体素点个数阈值为100。

上述技术方案中，预处理中所述对三角网格模型执行的拉普拉斯网格平滑处理，优先采用以下参数设置：拉普拉斯网格平滑指数为0.2，平滑次数为100。

上述技术方案中，预处理中所述三角网格简化处理，优先采用网格抽取(Decimation of Triangle Meshes)算法，网格抽取比例控制在0.9，即移除90％的三角网格单元。

上述技术方案中，步骤1中所述对三角网格单元单位法向量方向的调整，优先采用以下调整方式：

设计算得到单位法向量为(x,y,z)，若z＜0时，则调整单位法向量为(-x,-y,-z)。

上述技术方案中，步骤1中所述对三角网格单元的K-均值聚类，计算聚类中心，优先采用下述方式进行计算：

K-均值聚类算法流程如下：

(1)确定初始聚类中心；

(2)将所有三角网格单元分配到距其最近的聚类中心所属三角网格类；

(3)更新聚类中心；

(4)反复执行(2)(3)步直到满足终止条件。

K-均值聚类算法中相关计算公式如下所示：

流程(1)中初始聚类中心V_i的计算公式如下所示，其中i＝1,2,3：

其中σ为A_q中三角面积最大的三角网格单元的单位法向量；ξ为A_q中三角面积最小的三角网格单元的单位法向量。

令Δ_j为A_q中第j个三角网格单元，λ_j为Δ_j单位法向量，其中j∈[1,m]，m为A_q含有的三角网格单元个数；为A_q中第i个三角网格类，其聚类中心为V_i，其中i＝1,2,3；为第i个三角网格类中第j个三角网格单元，为单位法向量，为三角面积，其中j∈[1,n]，n为第i个三角网格类中三角网格单元个数。

流程(2)中三角网格单元Δ_j与各聚类中心V_i距离θ_ij的计算公式如下：

其中

流程(3)中聚类中心计算公式如下：

其中

Add定义为：

其中n为第i个三角网格类中三角网格单元个数。

流程(4)中迭代终止条件为：迭代次数超过50次或聚类中心变化幅度小于0.001；

第K次迭代聚类中心变化幅度值的计算公式如下：

其中为第k次迭代结束时聚类中心的三个坐标值。

上述技术方案中，步骤2中所述三角网格单元的重新分类，优先采用下述方法进行分类：

三角网格单元Δ_j与聚类中心V_i的距离为θ_ij，且

θ_kj＝min{θ_ij|i∈[1,2,3]} (10)

若θ_kj＜0.59359877，则Δ_j属于第k三角网格类否则Δ_j不属于任何一类。

上述技术方案中，步骤3中所述提取三角网格类的三角网格连通分量，优先采用下述方法：

若三角网格单元Δ_i的三角顶点与三角网格单元Δ_j的三角顶点有至少一个顶点相同，则Δ_i与Δ_j空间相邻。每个三角网格类中，按空间相邻关系提取三角网格连通分量其中t∈[1,k]，k是三角网格类的三角网格连通分量个数。三角网格类与三角网格连通分量的关系如下：

上述技术方案中，步骤4中所述对三角网格连通分量噪声剔除的方法，优先采用下述方法：

三角网格类的三角网格单元三角面积之和为三角网格连通分量的三角网格单元三角面积之和为其中t∈[1,k]，k是三角网格类的三角网格连通分量个数。计算公式如下：

其中为三角网格单元三角面积。

对中所有三角网格连通分量若

则将三角网格连通分量视为噪声干扰，并从三角网格类中剔除。

上述技术方案中，步骤4中所述合并相近三角网格类的方法，优先采用下述方法进行合并：

设聚类中心V_i与聚类中心V_j的距离为χ_ij，其计算公式如下：

其中

若χ_ij<0.7853981，则视两个聚类中心所属三角网格类相邻。将三个三角网格类中所有相邻的类合并为一个类。合并后的三角网格类仍记为其中1≤i≤3。

上述技术方案中，步骤6中所述封闭面结构集合的最小外接长方体计算，优先采用以下文献中所述的方法计算：三维岩心图像裂缝自动识别[期刊论文]邓知秋,滕奇志-《计算机与数字工程》2013年1期。

上述技术方案中，步骤7中所述形状因子、面积比的计算，优先采用下述方法：

①封闭面结构集合ψ的形状因子Φ_ψ计算公式如下：

②最小外接长方体形Π的形状因子Φ_Π计算公式如下：

③面积比Λ计算公式如下：

其中V_ψ为封闭面结构集合ψ的体积；S_ψ为ψ孔洞修补前所有三角网格单元三角面积之和；V_Π为封闭面结构集合最小外接长方体Π的体积；S_Π为Π的表面积。

上述技术方案中，步骤8中所述裂缝的判定标准，优先使用以下判定标准：

当下列条件均成立时：

或当以下条件均成立时：

其中

则判定当前封闭面结构集合ψ表面及内部的体素点集δ属于裂缝的一部分，并执行后续步骤；否则判定为非裂缝。

上述技术方案中，步骤9中所述三维图像形态学膨胀操作，优先采用以下方法：

其中B为中心位于(0,0,0)、尺寸为3*3的立方体；z为空间任意点(z₁,z₂,z₃)

(B)_z＝{c|c＝b+z,b∈B} (25)

上述技术方案中，步骤10中所述逻辑与操作，优先采用以下方法：

裂缝体素点集ζ计算公式为：

ζ＝{p|p∈H_q∩δ'} (26)

其中H_q为A_q对应的原始三维岩石孔隙模型连通分量的体素点集，p为任意体素点。δ'为步骤9中形态学膨胀处理后得到的体素点集。

本发明与现有技术相比，所具有的优点及有益的技术效果如下：

本发明所述的三维岩石图像裂缝识别方法，以实际的岩石二值序列图构建三维岩石孔隙模型，使用本发明所述的方法对其裂缝进行识别，并对识别结果进行观察，从而验证了本发明所用方法的可靠性与实用性。本发明解决了传统方法中裂缝识别不正确以及裂缝提取不准确的问题，本发明可以正确识别并精确提取岩石裂缝。

本发明受国家自然科学基金“岩石微观非均质结构三维图像重建及分辨率提升技术研究(61372174)”资助。

附图说明

图1-1是本发明实施例中岩石二值序列图之一；

图1-2是本发明实施例中由岩石二值序列图构造的三维岩石孔隙模型；

图2-1是本发明实施例中对图1-2中模型去噪后的三维岩石孔隙模型；

图2-2是本发明实施例中对图2-1中模型进行表面重建的结果；

图3-1是本发明实施例中对图2-2中模型执行拉普拉斯平滑后的局部放大图片；

图3-2是本发明实施例中对图2-2中模型执行拉普拉斯平滑后的局部放大图片；

图4-1是本发明实施例中对图3-1中模型执行网格抽取后的局部放大图片；

图4-2是本发明实施例中对图3-1中模型执行网格抽取后的局部放大图片；

图5是本发明实施例中对图4-1中模型计算网格单位法向量的局部示意图；

图6-1是本发明实施例中对图5中模型的三角网格单元K-均值聚类后重新分类的结果图；

图6-2是本发明实施例中对图5中模型的三角网格单元K-均值聚类后重新分类的结果图；

图7-1是本发明实施例中对图6-1中分类后的三角网格单元执行连通分量标记的结果图；

图7-2是本发明实施例中对图6-1中分类后的三角网格单元执行连通分量标记的结果图；

图7-3是本发明实施例中对图6-1中分类后的三角网格单元执行连通分量标记的结果图；

图8-1是本发明实施例中三维岩石孔隙模型裂缝识别的效果图；

图8-2是本发明实施例中图8-1中第一条裂缝的最小外接长方体示意图；

图8-3是本发明实施例中图8-1中第二条裂缝的最小外接长方体示意图；

图8-4是本发明实施例中图8-1中第三条裂缝的最小外接长方体示意图；

图8-5是本发明实施例中三维岩石孔隙模型裂缝识别的效果图；

图8-6是本发明实施例中图8-5中裂缝与原始孔隙同时显示的效果图；

图8-7是本发明实施例中三维岩石孔隙模型裂缝识别的效果图；

图8-8是本发明实施例中图8-7中裂缝与原始孔隙同时显示的效果图；

图8-9是本发明实施例中三维岩石孔隙模型裂缝识别的效果图；

图8-10是本发明实施例中图8-9中裂缝与原始孔隙同时显示的效果图；

图9是本发明所述的方法流程图。

具体实施方式

下面通过具体实施例并结合附图对本发明作进一步详细说明，有必要在此指出的是所述实施例只是用于对本发明的进一步描述，但并不意味着是对本发明保护范围的任何限定。

实施例：为了使本发明所述识别方法更加便于理解和接近于真实应用，下面从岩石二值序列图像开始到三维岩石孔隙模型裂缝识别完成为止，对整个操作流程进行说明。具体操作步骤如下：

(1)准备一组待识别的岩石二值序列图，图片数量为256，格式为BMP，尺寸256×256，其中代表如图1-1所示；由岩石二值序列图构建的三维岩石孔隙模型，如图1-2所示,图中灰色部分为岩石孔隙。

(2)提取三维岩石孔隙模型的所有连通分量，剔除其中体素点个数小于100的连通分量，去噪结果如图2-1所示。用移动立方体等值面提取(Marching Cube)算法对每个连通分量进行表面重建，得到每个连通分量相应的三角网格模型A_q，其中q＝1，2，3……k，k为三角网格模型个数。表面重建结果如图2-2所示。

(3)对每个表面重建结果A_q执行平滑指数0.2，平滑次数100的拉普拉斯网格平滑处理；网格平滑效果如模型局部放大图3-1和3-2所示。

(4)使用三角网格抽取算法对平滑后的三角网格模型执行网格简化处理，网格抽取比例控制在0.9。网格简化效果如模型局部放大图4-1和4-2所示，图中直线段为三角网格的边界。

(5)计算所有三角网格单元的单位法向量，设计算得到单位法向量为(x,y,z)，若z＜0时，则调整单位法向量为(-x,-y,-z)。单位法向量计算结果如图5所示，图中三角网格单元上的箭头为其单位法向量。

(6)利用三角网格单元的单位法向量及三角网格单元三角面积对A_q中的三角网格单元进行K-均值聚类，得到三个聚类中心V₁、V₂、V₃。

K-均值聚类算法流程如下：

1)求A_q中三角面积最大的三角网格单元的单位法向量σ及A_q中三角面积最小的三角网格单元的单位法向量ξ。则令σ、σ×ξ、σ×(σ×ξ)分别为初始的三个聚类中心；

2)计算三角网格单元Δ_j与各聚类中心V_i距离θ_ij，计算公式为：

其中

其中λ_j为三角网格Δ_j的单位法向量。

将三角网格单元分配到距离其最近的聚类中心所属三角网格类。

3)令Δ_j为A_q中第j个三角网格单元，λ_j为Δ_j单位法向量，其中j∈[1,m]，m为A_q含有的三角网格单元个数；为A_q中第i个三角网格类，其聚类中心为V_i，其中i＝1,2,3；为第i个三角网格类中第j个三角网格单元，为单位法向量，为三角面积，其中j∈[1,n]，n为第i个三角网格类中三角网格单元个数；

重新计算聚类中心V_i，计算公式如下：

其中

Add定义为：

其中n为第i个三角网格类中三角网格单元个数。

4)计算第k次迭代结束时，聚类中心变化幅度值

其中为第k次迭代结束时聚类中心的三个坐标值。

若迭代次数k超过50或者聚类中心变化幅度小于0.001则结束迭代，否则返回流程2)。

(7)计算每个三角网格单元Δ_j与聚类中心V_i的距离θ_ij，且

θ_kj＝min{θ_ij|i∈[1,2,3]} (9)

若θ_kj＜0.59359877则将Δ_j归为第k三角网格类否则Δ_j不属于任何一类。

(8)提取三角网格类的三角网格连通分量：若三角网格单元Δ_i的三角顶点与Δ_j的三角顶点至少有一个顶点相同，则Δ_i与Δ_j空间相邻。每个三角网格类中，按空间相邻关系提取三角网格连通分量其中t∈[1,k]，k是三角网格类的三角网格连通分量个数。三角网格连通分量提取结果如图7-1、7-2和7-3所示。

(9)计算每个三角网格类中所有三角网格单元的三角面积和计算三角网格连通分量中所有三角网格单元的三角面积和其中t∈[1,k]，k是三角网格类的三角网格连通分量个数。计算公式分别为：

其中为三角网格单元的三角面积。

对中所有三角网格连通分量若

则将此三角网格连通分量视为噪声干扰，并从三角网格类中剔除。

剔除噪声后，重新计算每个类的聚类中心，并合并相近的类。计算聚类中心V_i与聚类中心V_j的距离χ_ij，计算公式如下：

其中

若χ_ij<0.7853981，则视两个聚类中心所属三角网格类相邻。将三个三角网格类中所有相邻的类合并为一个类。合并后的三角网格类仍记为其中1≤i≤3。相邻类的合并结果如图6-1、6-2所示，图中不同灰度分别表示一个类，白色三角网格不属于任何一类。

(10)三维岩石孔隙连通分量A_q中的每个三角网格类中的所有三角网格单元组成一个三维面结构。对这个三维面结构执行三维模型孔洞修补操作，操作后得到一个或多个封闭的三维面结构，这些三维面结构组成一个封闭面结构集合。

(11)按如下文献计算封闭面结构集合的最小外接长方体：三维图像裂缝自动识别[期刊论文]邓知秋,滕奇志-《计算机与数字工程》2013年1期。

(12)计算封闭面结构集合ψ的形状因子Φ_ψ：

计算最小外接长方体形Π的形状因子Φ_Π：

计算面积比Λ：

其中V_ψ为封闭面结构集合ψ的体积；S_ψ为ψ孔洞修补前的所有三角网格单元三角面积之和；V_Π为封闭面结构集合最小外接长方体Π的体积；S_Π为Π的表面积。

(13)当下列条件均成立时：

或者当以下条件均成立时：

其中

则判定当前封闭面结构集合ψ表面及内部包含的体素点集属于裂缝的一部分，并执行下一步骤；否则判定为非裂缝。

(14)提取(13)中判定为真的封闭面结构集合的表面以及内部的所有体素点，对提取得到的体素点集进行三维图像形态学膨胀处理。

(15)：将(14)中膨胀得到的体素点集与原始三维岩石孔隙模型连通分量的体素点集进行逻辑与操作，操作后得到的体素点集即为裂缝。与操作结果如图8-1至8-10所示：

图8-1为提取的裂缝示意图，图中共有3个裂缝；

图8-2至8-4为每个裂缝的最小外接长方体示意图；

图8-5至8-10为裂缝不同方位显示图，以及将其与原始孔隙同时显示的效果图。

本实施例中，以实际的岩石二值序列图构建三维模型并使用本发明的方法对其进行裂缝识别。对识别结果进行观察，从而验证了本发明识别方法的可靠性与准确性。

本发明上述实施例只是本发明的优选实施例，并不是对本发明所述技术方案的任何限制，只要是不经过创造性劳动即在上述实施例的基础上实现的技术方案，均应视为落入本发明所保护的范围内。