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

技术特征：

1.岩石三维图像裂缝提取方法，其特征在于经过图像预处理后，算法包括以下步骤：

步骤1：计算预处理得到的三角网格模型A_q中每个三角网格单元的单位法向量与三角面积，并调整单位法向量方向,其中q＝1，2，3…k，k为三角网格模型个数；接着利用三角网格单元的单位法向量及三角面积特征对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中已膨胀的体素点集与原始三维岩石孔隙模型连通分量的体素点集进行逻辑与操作，操作后得到的体素点集即为裂缝。

2.根据权利要求1所述的方法，其特征在于步骤1中所述基于网格单位法向方向及网格三角面积特征的K-均值聚类，计算三角网格单元聚类中心，采用下述计算方法：

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

①确定初始聚类中心；

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

③更新聚类中心；

④反复执行流程②，③直到满足终止条件；

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

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

$\left\{\begin{array}{c}{V}_1=\mathrm{\σ}\\ V_2=\mathrm{\σ}\×\mathrm{\ξ}\\ V_3=\mathrm{\σ}\×(\mathrm{\σ}\×\mathrm{\ξ})\end{array}\right.---\left(1\right)$

其中σ为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个三角网格类中三角网格单元个数；

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

${\mathrm{\θ}}_{i,j}=\left\{\begin{array}{cc}\mathrm{\π}-{\mathrm{\φ}}_{ij}& {\mathrm{\φ}}_{ij}\>\mathrm{\π}/2\\ {\mathrm{\φ}}_{ij}& {\mathrm{\φ}}_{ij}\≤\mathrm{\π}/2\end{array}\right.---\left(2\right)$

其中

${\mathrm{cos\φ}}_{ij}=\frac{{\mathrm{\λ}}_j{V}_i}{\left|{\mathrm{\λ}}_j\right|\left|V_i\right|}---\left(3\right)$

流程③中聚类中心计算公式如下：

$V_i=\frac{{\mathrm{\β}}_n^i}{\left|{\mathrm{\β}}_n^i\right|}--\left(4\right)$

其中

且Add定义为：

$Add(J_1,J_2)=\left\{\begin{array}{cc}{J}_1+J_2& |J_1+J_2|\≥|J_1-J_2|\\ J_1-J_2& |J_1+J_2|\<|J_1-J_2|\end{array}\right.---\left(6\right)$

${\mathrm{\ω}}_j^i=\left\{\begin{array}{cc}5& S_j^i\>S_{max}^i/2\\ 1& S_j^i\≤S_{max}^i/2\end{array}\right.---\left(7\right)$

$S_{max}^i=max\left\{S_j^i\right|j\∈\[1,n\]\}---\left(8\right)$

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

流程④中迭代终止条件为：迭代次数超过50次或者聚类中心变化幅度小于0.001；第K次迭代聚类中心变化幅度值的计算公式如下：

${\∂}^k=\underset{i=1}{\overset{3}{\Σ}}\left(\right|x_{V_i^k}-x_{V_i^{k-1}}|+|y_{V_i^k}-y_{V_i^{k-1}}|+|z_{V_i^k}-z_{V_i^{k-1}}\left|\right)---\left(9\right)$

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

3.根据权利要求1所述的方法，其特征在于步骤2所述三角网格单元分类方法,采用下述计算方法：

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

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

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

4.根据权利要求1所述的方法，其特征在于步骤4所述三角网格连通分量噪声剔除及合并相近三角网格类，采用下述计算方法：

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

$S_{{\mathrm{\η}}_t^i}=\underset{{\mathrm{\Δ}}_j^i\∈{\mathrm{\η}}_t^i}{\Σ}{S}_j^i---\left(12\right)$

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

对中三角网格连通分量若

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

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

${\mathrm{\χ}}_{ij}=\left\{\begin{array}{cc}\mathrm{\π}-{\mathrm{\φ}}_{ij}& {\mathrm{\φ}}_{ij}\>\mathrm{\π}/2\\ {\mathrm{\φ}}_{ij}& {\mathrm{\φ}}_{ij}\≤\mathrm{\π}/2\end{array}\right.---\left(14\right)$

其中

${\mathrm{cos\φ}}_{ij}=\frac{V_i{V}_j}{\left|V_i\right|\left|V_j\right|}---\left(15\right)$

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

5.根据权利要求1所述的方法，其特征在于步骤8所述裂缝的判定,采用以下判定标准：

当下列条件均成立时：

$\left\{\begin{array}{c}{\mathrm{\&Phi;}}_{\&Pi;}\&le;M\\ \mathrm{\&Lambda;}\&GreaterEqual;E\\ {\mathrm{\&Phi;}}_{\mathrm{\&psi;}}\&le;0.05\end{array}\right.---\left(16\right)$

或当以下条件均成立时：

$\left\{\begin{array}{c}{\mathrm{\&Phi;}}_{\&Pi;}\&le;0.04\\ \mathrm{\&Lambda;}\&GreaterEqual;0.17\\ {\mathrm{\&Phi;}}_{\mathrm{\&psi;}}\&le;0.05\end{array}\right.---\left(17\right)$

其中

$M=0.6*\frac{{\log }_2{V}_{\&Pi;}}{20}---\left(18\right)$

$E=\left\{\begin{array}{cc}0.22& M\&lt;0.45\\ 0.17& M\&GreaterEqual;0.45\end{array}\right.---\left(19\right)$

${\mathrm{\&Phi;}}_{\mathrm{\&psi;}}=\frac{36{{\mathrm{\&pi;V}}_{\mathrm{\&psi;}}}^2}{{S_{\mathrm{\&psi;}}}^3}---\left(20\right)$

${\mathrm{\&Phi;}}_{\&Pi;}=\frac{36{{\mathrm{\&pi;V}}_{\&Pi;}}^2}{{S_{\&Pi;}}^3}---\left(21\right)$

$\mathrm{\&Lambda;}=\frac{S_{\mathrm{\&psi;}}}{S_{\&Pi;}}--\left(22\right)$

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

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

6.根据权利要求1所述的方法，其特征在于步骤10所述裂缝的准确提取,使用以下提取方法：

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

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

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