CN108648256B - 一种基于超维的灰度岩心三维重建方法 - Google Patents

Abstract

本发明公开了一种基于超维的灰度岩心三维重建方法。包括以下步骤：提出灰度超维重建的字典训练算法，利用特征提取、模式提取和聚类算法建立二维与三维对应的多层字典，作为重建的先验信息；在重建过程中提出二值约束的邻域匹配算法，对不同位置的图像块采取不同的距离计算方式，选取字典中匹配的块作为重建结果；对岩心灰度图像进行重建实验，将重建结果与训练图像、真实岩心图像的统计特性和形态特征进行比较，以证明算法的有效性。岩心的灰度重建比二值重建能够保留更多的原始信息，本方法在进行匹配块搜索时，增加了权值计算，提高了匹配的鲁棒性。保证了重建结果与原始训练图像在统计特征上的一致性。

Description

一种基于超维的灰度岩心三维重建方法

技术领域

本发明创新性的提出了基于超维的灰度岩心三维重建方法，尤其涉及一种 灰度岩心的三维重建方法，属于图像处理领域。

背景技术

在石油地质研究中，岩心三维微观结构是研究岩心宏观物理特性的基础。 三维数字岩心能够在孔隙尺寸级别上反映孔隙微观结构，计算岩心声学、电学 特性和模拟渗流过程，是分析岩心微观物理特性的有力工具。然而实际工程中 高精度岩心三维图像难以直接获取，利用高精度二维图像进行数字建模，可以 有效重建岩心三维图像。

虽然利用成像设备获取的图像大多是灰度图像，但目前较多的三维重建算 法是针对二值化图像。重建时将灰度图像进行二值化，然后重建二值三维图像。 然而灰度图像比二值图像能反映更多信息，重建灰度图像后可以进行二值化， 但重建的二值图像无法恢复成灰度图像。为了更多保留原始信息，对于灰度图 像的三维重建算法的研究是很有必要的。

目前的灰度岩心重建算法主要有：Tahmasebi等提出的基于互相关函数的模 拟算法(CCSIM，cross-correlation–based simulation)，该方法能够很好继承层与层 之间的连续性，但是不易控制层与层之间的随机变化性。以及纹理合成算法， 其基本思想是以小规模纹理作为样本，合成较大尺寸的结果纹理。

但是上述方法都没有考虑岩心的真实三维结构，超维重建是将已有的真实 岩心图像作为先验信息，利用真实的三维图像信息指导重建过程。本发明提出 了一种基于超维的灰度岩心三维重建方法，保证了重建后三维灰度岩心图像的 统计特征与二维原始岩心图像保持一致，形态特征保持相似性。该研究项目受 国家自然科学基金项目《岩石微观非均质结构三维图像重建及分辨率提升技术 研究》(61372174)资助。

发明内容

本发明的目的在于创新性的提出一种新的灰度岩心三维重建方法，并且保 证重建后三维灰度岩心图像在统计特征上与二维原始岩心图像保持一致。

本发明通过以下技术方案来实现上述目的：

(1)将真实岩心灰度图像I二值化得到I_bw，孔隙相为1，岩石相为0；设 定模板尺寸为n，将真实岩心灰度图像I和二值图像I_bw分成一一对应的h×w×n的 三维子图；

(2)对每一个灰度的三维子图的底面提取特征，得到features；利用模式 提取方法以光栅路径扫描所有三维子图，将features作为灰度二维模式，其对应 的灰度三维图像块作为灰度三维模式，得到灰度模式集Patternsetgray，对二值化 后的三维子图的底面提取二值二维模式，对应的二值的三维图像块作为二值三 维模式，得到二值模式集Patternsetbw；

(3)设定类别数M，使用Kmeans聚类算法对Patternsetgray中的灰度二维模 式进行聚类分析，将模式集分为M类；

(4)将Patternsetgray中的灰度二维模式及其对应的灰度三维模式，Patternsetbw二值二维模式及其二值三维模式按照所属类别对应保存，输出所建立 的字典。

(5)对于待重建的灰度岩心二维图像，每次重叠一行或者一列提取待重建 图像的二维模式{Pattern_2d1,Pattern_2d2,…,Pattern_2di,…,Pattern_2dN}；

(6)对第i个二维模式Pattern_2di计算该模式与训练好的字典中的M个类的 聚类中心的距离D_c＝{d₁,d₂,…d_i,…,d_M}，找到D_c中的最小值D_min；

(7)在D_min对应的类中进行匹配，假设该类中有p个二维与其对应三维的 字典原子，计算Pattern_2di与该类中的所有字典原子的灰度二维模式的距离 D_bg＝{d_bg1,d_bg2,…d_bgi,…,d_bgp}，将Pattern_2di所在的二维图像块二值化，得到Pattern_2dbwi， 计算Pattern_2di与该类中的所有字典原子的二值二维模式的距离 D_bbw＝{d_bbw1,d_bbw2,…d_bbwi,…,d_bbwp}，将D_bg和D_bbw加权得到二维模式距离 D_b＝aD_bg+bD_bbw；

(8)在字典原子中找到满足条件的块，对D_b进行升序排列，距离最小的块 对应的n×n×n灰度三维模式即为匹配的三维图像块，依次完成对整个三维岩心 块的重建。

上述方法的基本原理如下：

超维重建(SDR)是一种全新的三维重建方法，该方法主要借鉴了基于学习 的超分辨率重建算法的思想，将已有的真实岩心三维图像作为先验信息，学习 二维图像块到三维图像块的映射关系，建立对应的字典，在重建时，利用待重 建图像的信息和在先验信息中匹配的部分进行重建，重建出来的三维结构在统 计信息与形态结构与原始图像都比较接近。该算法合理地运用了已有的图像资 源，将成像设备获取三维图像与数学模拟方法结合，解决了利用成像设备获得 大量三维图像价格昂贵和传统三维重建算法在形态上与原始图像差距太大的缺 点。本方法在超维重建的框架下提出了灰度岩心图像超维重建算法，实现了灰 度岩心图像二维到三维的重建，如图1所示。

所述步骤(1)中，采用通过CT扫描获取的真实岩心灰度三维图，并对其 进行二值化得到与之对应的二值岩心图像。

所述步骤(2)中，建立字典时，需要对真实岩心图像进行模式提取，然后 按照一定准则保存二维模式及其对应的三维模式信息。设模板尺寸为n，I(x,y,z) 表示真实岩心三维结构在点(x,y,z)的灰度值，用n×n的模板对I(x,y,z₀)采样， n×n×n的模板对I(x,y,z₀+n)采样，如图2所示为5×5的二维模板和5×5×5的三维 模板在图像上采样，若采样中心为(x₀,y₀,z₀),则得到二维模式Pattern2d及其对应 的三维模式Pattern3d。

Pattern2d＝I(x₀±n/2,y₀±n/2,z₀) (1)

Pattern3d＝I(x₀±n/2,y₀±n/2,z₀±n)

固定模板尺寸对整个真实岩心三维结构按照光栅路径扫描，即可获得灰度 图像二维及其三维结构的模式集Patternset，N表示模板对的个数。为了保证 模式集的连续性与差异性，扫描时中心点位置每次移动n-1，两个块之间只有一 个面相互重叠即相邻块之间重叠一个面，比如扫描时第一个中心点位置为 (x₀,y₀,z₀)，则下一个中心点的位置为(x₀+n-1,y₀,z₀)。

在进行重建时，更加关注结构信息，若直接用图像的灰度值作为模式信息， 匹配条件苛刻，容易出现在字典中无法找到匹配块的情况。岩心CT图像灰度级 反映的是该岩心的成分对X射线的吸收能力，同一块岩心成分相似，可以用待 重建二维图像的灰度信息进行重建，而在灰度岩心重建的过程中更关注结构信 息，因此需要提取二维结构特征来进行字典训练和重建。在岩心图像中，岩石 颗粒和孔隙有着明显的灰度差异，在边缘有明显的灰度跳变，提取梯度信息可 以反映结构信息。基于此本方法提出对二维图像提取一阶梯度和二阶梯度特征， 将特征作为匹配条件。

对于灰度图像f的(x,y)位置处的一阶梯度▽f用向量定义为：

处理的图像是数字量，对于一幅图像的3x3区域如图3(a)所示，z项表 示灰度值，本文使用大小为3x3的模板来近似偏导数，其数学近似有下式给出：

公式(4)可用图3(b)中的两个模板通过滤波整个图像来实现，本文使用 该算子提取一阶梯度。

二阶微分能够体现灰度级变换的方向，拉普拉斯算子是常用的二阶微分算 子，定义为：

其中：

对于一个如图3(a)的3x3的区域：

拉普拉斯算子是无方向性的算子，计算简单，可以用模板与图像进行滤波 计算实现。在本方法中采用如图4的模板对图像提取二阶梯度，同一阶梯度一 起作为灰度二维图像特征。

将features作为灰度二维模式，其对应的灰度三维图像块作为灰度三维模 式，得到灰度模式集Patternsetgray。对二值化后的三维子图的底面提取二值二维 模式，对应的二值的三维图像块作为二值三维模式，得到二值模式集Patternsetbw。

所述步骤(3)中，为了使模式集更加丰富更加完备，将使用较大的真实岩 心图像来训练字典，得到一个非常大的模式集，若直接将获取的模式集作为重 建的解空间，重建时为每一个待重建的二维图像块在整个模式集中搜索一个相 似的二维模式将耗费大量时间，导致重建过程十分漫长，因此本方法提出字典 分类的字典训练方法，提高重建效率。该方法主要利用了经典的Kmeans聚类算 法，对二维模式Pattern2d进行聚类，将模式集划分为不同的类，相同的类的模 式具有相似性，不同类的模式差异较大，将模式集划分到类θ^M中的公式为：

θ^M(Patternset_i)＝ψ(Patternset_i(Pattern2d_i))i＝1,2,…N (8)

式中ψ表示聚类算法，M表示类别数。Kmeans算法计算简单，使用广泛， 可提前预设聚类的类别数，通过计算样本与聚类中心的距离，将样本划分到与 聚类中心距离最近的类中。在灰度超维重建的字典训练中，通过聚类分析，有 助于重建时快速定位待重建的二维模式属于哪一类，然后在该类中去搜索相似 的模式。若聚类个数设为m，则将整个字典按照聚类中心分为m个子空间，将 二维及其对应的三维模式划分到归属的类，作为该子空间的字典原子，字典的 结构示意如图5所示。

所述步骤(4)中，将Patternsetgray中的灰度二维模式及其对应的灰度三维模式，Patternsetbw二值二维模式及其二值三维模式按照所属类别对应保存，得到根 据真实岩心所建立的字典。

所述步骤(5)中，提取待重建图像的二维灰度模式时，为了增加块与块之 间的连续性，能够有效防止重建时产生块效应，所以相邻的模式之间重叠一行 或者一列得到{Pattern_2d1,Pattern_2d2,…,Pattern_2di,…,Pattern_2dN}；

所述步骤(6)中，为了提高重建效率，对于待重建图像的每一个二维灰度 模式，首先在训练好的字典中搜索与该模式距离最近的类。

所述步骤(7)中，对于在步骤(6)中找到的与当前模式距离最近的类中， 如果有p个二维与其对应三维的字典原子，计算Pattern_2di与该类中的所有字典原 子的灰度二维模式的距离D_bg＝{d_bg1,d_bg2,…d_bgi,…,d_bgp}，将Pattern_2di所在的二维图像 块二值化，得到Pattern_2dbwi，计算Pattern_2di与该类中的所有字典原子的二值二维模 式的距离D_bbw＝{d_bbw1,d_bbw2,…d_bbwi,…,d_bbwp}，将D_bg和D_bbw加权得到二维模式距离 D_b＝aD_bg+bD_bbw。在本方法中，我们使用欧式距离来进行不同模式间相似性估算， 两个n维向量x＝(x₁,…,x_n)和y＝(y₁,…,y_n)之间的欧式距离为：

重建时以欧式距离来度量二维模式的相似度，距离越小，表示越相似。

在步骤(8)中，在字典原子中找到满足条件的块，具体的对于不同的位置， 采用如下方法：

i.当Pattern_2di是待重建图像的第一个提取的模式时，对D_b进行升序排列，距 离最小的块对应的n×n×n灰度三维模式即为匹配的三维图像块；

ii.当Pattern_2di是待重建图像的前n行提取的模式时(除去第一个)，对该块 进行重建时，通过计算字典中所有原子的n×n×n三维模式与其左邻域 已重建好的n×n×n三维图像块的距离D_lg＝{d_lg1,d_lg2,…d_lgi,…,d_lgp}，计算字 典原子中二值三维模式与将左邻域的块二值化后的二值三维图像块的 距离D_lbw＝{d_lbw1,d_lbw2,…d_lbwi,…,d_lbwp}，将D_lg和D_lbw加权得到左邻域距离 D_l＝cD_lg+dD_lbw，将左邻域距离与二维模式距离加权得到新的距离D，D＝αD_b+βD_l，α+β＝1，对D进行升序排列，距离最小的块对应的 n×n×n灰度三维模式即为匹配的三维图像块；

iii.当Pattern_2di是待重建图像的前n列提取的模式时(除去第一个)，对该块 进行重建时，通过计算字典中所有原子的n×n×n三维模式与其后邻域 已重建好的n×n×n三维图像块的距离D_u＝{d_u1,d_u2,…d_ui,…,d_up}，计算字典 原子中二值三维模式与将后邻域的块二值化后的二值三维图像块的距 离D_ubw＝{d_ubw1,d_ubw2,…d_ubwi,…,d_ubwp}，将D_ug和D_ubw加权得到后邻域距离 D_u＝eD_ug+fD_ubw，将后邻域距离与二维模式距离加权得到新的距离D， D＝αD_b+γD_u，α+γ＝1，对D进行升序排列，距离最小的块对应的 n×n×n灰度三维模式即为匹配的三维图像块；

iv.当Pattern_2di是待重建图像的其余部分提取的二维模式时，对该块进行重 建时，通过计算字典中所有原子的n×n×n三维模式与其左邻域已重建 好的n×n×n三维图像块的距离D_1g＝d_1g1，d_1g2，…d_1gi，…，d_1gp}，计算字典原子 中二值三维模式与将左邻域的块二值化后的二值三维图像块的距离 D_lbw＝{d_lbw1，d_lbw2，…d_lbwi，…，d_lbwp}，将D_lg和D_lbw加权得到左邻域距离 D_l＝cD_lg+dD_lbw，计算字典中所有原子的n×n×n三维模式与其后邻域已 重建好的n×n×n三维图像块的距离D_u＝{d_u1，d_u2，…d_ui，…，d_up}，计算字典原 子中二值三维模式与将后邻域的块二值化后的二值三维图像块的距离 D_ubw＝{d_ubw1，d_ubw2，…d_ubwi，…，d_ubwp}，将D_ug和D_ubw加权得到后邻域距离 D_u＝eD_ug+fD_ubw，最终与二维模式距离加权得到新的距离D， D＝αD_b+βD_l+γD_u，α+β+γ＝1，对D进行升序排列，距离最小的块对 应的n×n×n灰度三维模式即为匹配的三维图像块。

附图说明

图1是灰度岩心超维重建过程

图2是二维及其三维模式提取过程

图3是图像3x3区域及本文所用模板，其中(a)为图像3x3区域；(b)为本文 所用一阶梯度模板

图4是本文采用的图像二阶梯度模板

图5是字典结构示意图

图6是超维重建实验训练图像

图7-1是二维测试图像的灰度图

图7-2是二维测试图像的二值图

图8是测试图像灰度重建结果图

图9-1是测试图像灰度直方图

图9-2是原始三维结构灰度直方图

图9-3是超维重建结果灰度直方图

图10-1是重建岩心x方向两点相关函数

图10-2是重建岩心y方向两点相关函数

图10-3是重建岩心z方向两点相关函数

图11-1是重建岩心x方向线性路径函数

图11-2是重建岩心y方向线性路径函数

图11-3是重建岩心z方向线性路径函数

具体实施方式

下面结合具体实施例和附图对本发明作进一步说明：

实施例：

(1)在本例中，使用如图6的600×600×600的灰度岩心三维CT图像I作为 训练图像，得到600张无信息确实的真实岩心CT图像并存储在计算机中。并且 将此灰度图像孔隙相为1，岩石相为0，得到二值化后的岩心图像I_bw，设定模板 尺寸为3，将真实岩心灰度图像I和二值图像I_bw分成一一对应的600×600×3的三 维子图；

(2)对每一个灰度的三维子图的底面提取特征，得到features；利用模式 提取方法以光栅路径扫描所有三维子图，将features作为灰度二维模式，其对应 的灰度三维图像块作为灰度三维模式，得到灰度模式集Patternsetgray，对二值化 后的三维子图的底面提取二值二维模式，对应的二值的三维图像块作为二值三 维模式，得到二值模式集Patternsetbw；

(3)设定类别数M，使用Kmeans聚类算法对Patternsetgray中的灰度二维模 式进行聚类分析，将模式集分为M类；

(4)将Patternsetgray中的灰度二维模式及其对应的灰度三维模式，Patternsetbw二值二维模式及其二值三维模式按照所属类别对应保存，输出字典。

(5)使用二维灰度图像进行三维重建，图7-1为待重建灰度图像，图7-2 为待重建二值图像。每次重叠一行或者一列提取待重建图像的二维模式 {Pattern_2d1，Pattern_2d2，…，Pattern_2di，…，Pattern_2dN}；

(6)对第i个二维模式Pattern_2di计算该模式与训练好的字典中的M个类的 聚类中心的距离D_c＝{d₁,d₂,…d_i,…,dM}，找到D_c中的最小值D_min；

(7)在D_min对应的类中进行匹配，假设该类中有p个二维与其对应三维的 字典原子，计算Pattern_2di与该类中的所有字典原子的灰度二维模式的距离 D_bg＝{d_bg1,d_bg2,…d_bgi,…,d_bgp}，将Pattern_2di所在的二维图像块二值化，得到Pattern_2dbwi， 计算Pattern_2di与该类中的所有字典原子的二值二维模式的距离 D_bbw＝{d_bbw1,d_bbw2,…d_bbwi,…,d_bbwp}，将D_bg和D_bbw加权得到二维模式距离 D_b＝aD_bg+bD_bbw；

(8)按照如下准则在字典原子中找到满足条件的块：

i.当Pattern_2di是待重建图像的第一个提取的模式时，对D_b进行升序排列，距 离最小的块对应的3×3×3灰度三维模式即为匹配的三维图像块；

ii.当Pattern_2di是待重建图像的前n行提取的模式时(除去第一个)，对该块 进行重建时，通过计算字典中所有原子的3×3×3三维模式与其左邻域已 重建好的3×3×3三维图像块的距离D_lg＝{d_lg1,d_lg2,…d_lgi,…,d_lgp}，计算字典 原子中二值三维模式与将左邻域的块二值化后的二值三维图像块的距 离D_lbw＝{d_lbw1,d_lbw2,…d_lbwi,…,d_lbwp}，将D_lg和D_lbw加权得到左邻域距离 D_l＝cD_lg+dD_lbw，将左邻域距离与二维模式距离加权得到新的距离D，D＝αD_b+βD_l，α+β＝1，对D进行升序排列，距离最小的块对应的3×3×3 灰度三维模式即为匹配的三维图像块；

iii.当Pattern_2di是待重建图像的前n列提取的模式时(除去第一个)，对该块 进行重建时，通过计算字典中所有原子的3×3×3三维模式与其后邻域 已重建好的n×n×n三维图像块的距离D_u＝{d_u1,d_u2,…d_ui,…,d_up}，计算字典 原子中二值三维模式与将后邻域的块二值化后的二值三维图像块的距 离D_ubw＝{d_ubw1,d_ubw2,…d_ubwi,…,d_ubwp}，将D_ug和D_ubw加权得到后邻域距离 D_u＝eD_ug+fD_ubw，将后邻域距离与二维模式距离加权得到新的距离D， D＝αD_b+γD_u，α+γ＝1，对D进行升序排列，距离最小的块对应的3×3×3 灰度三维模式即为匹配的三维图像块；

iv.当Pattern_2di是待重建图像的其余部分提取的二维模式时，对该块进行重建时，通过计算字典中所有原子的3×3×3三维模式与其左邻域已重建好的 3×3×3三维图像块的距离D_lg＝{d_lg1,d_lg2,…d_lgi,…,d_lgp}，计算字典原子中二 值三维模式与将左邻域的块二值化后的二值三维图像块的距离 D_lbw＝{d_lbw1,d_lbw2,…d_lbwi,…,d_lbwp}，将D_lg和D_lbw加权得到左邻域距离D_l＝cD_lg+dD_lbw，计算字典中所有原子的n×n×n三维模式与其后邻域已 重建好的n×n×n三维图像块的距离D_u＝{d_u1,d_u2,…d_ui,…,d_up}，计算字典原 子中二值三维模式与将后邻域的块二值化后的二值三维图像块的距离 D_ubw＝{d_ubw1,d_ubw2,…d_ubwi,…,d_ubwp}，将D_ug和D_ubw加权得到后邻域距离 D_u＝eD_ug+fD_ubw，最终与二维模式距离加权得到新的距离D， D＝αD_b+βD_l+γD_u，α+β+γ＝1，对D进行升序排列，距离最小的块对 应的3×3×3灰度三维模式即为匹配的三维图像块。

(9)按照步骤(8)的准则，依次完成整幅岩心图像的三维重建，重建结 果如图8所示。

(10)在地质统计学中，通常使用两点相关函数、线性路径函数等统计 特征来描述孔隙空间。因此，在三维重建算法中通常通过对比重建结果与待重 建图像和原始三维结构的统计特征函数来验证重建结果的有效性。在本方法中 我们使用两点相关函数以及先行路径函数来验证本方法的有效性。

从图8的重建结果中可以看出本方法的重建结果在z方向上具有连续性。 如图9-1是测试图像灰度直方图，图9-2是原始三维结构灰度直方图，图9-3是 超维重建结果灰度直方图。从上述图中可以看出该算法能有效重建灰度图像， 重建结果的灰度直方图与待重建图像、原始三维图像的灰度直方图一致。将待 重建图像、原始三维结构和重建的三维结构二值化，然后计算各自的两点相关 函数和线性路径函数，图10-1是重建岩心x方向两点相关函数，图10-2是重建 岩心y方向两点相关函数，图10-3是重建岩心z方向两点相关函数，图11-1是 重建岩心x方向线性路径函数，图11-2是重建岩心y方向线性路径函数，图11-3是重建岩心z方向线性路径函数。从上述图中可以看到本文算法重建结果的两点 概率函数和线性路径函数与原始三维结构、待重建图像相似，证明该算法重建 结果的统计特征与真实灰度岩心三维图像相一致，能取得较好的重建效果，证 明了本方法的有效性。

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

Claims (2)

1.一种基于超维的灰度岩心三维重建方法，其特征在于：包括以下步骤：

(1)将真实岩心灰度图像I二值化得到I_bw，孔隙相为1，岩石相为0；设定模板尺寸为n，将真实岩心灰度图像I和二值图像I_bw分成一一对应的h×w×n的三维子图；

(2)对每一个灰度的三维子图的底面提取特征，得到features；利用模式提取方法以光栅路径扫描所有三维子图，将features作为灰度二维模式，其对应的灰度三维图像块作为灰度三维模式，得到灰度模式集Patternsetgray，对二值化后的三维子图的底面提取二值二维模式，对应的二值的三维图像块作为二值三维模式，得到二值模式集Patternsetbw；

(3)设定类别数M，使用Kmeans聚类算法对Patternsetgray中的灰度二维模式进行聚类分析，将模式集分为M类；

(4)将Patternsetgray中的灰度二维模式及其对应的灰度三维模式，Patternsetbw二值二维模式及其二值三维模式按照所属类别对应保存，输出所建立的字典；

(5)对于待重建的灰度岩心二维图像，每次重叠一行或者一列提取待重建图像的二维模式{Pattern_2d1,Pattern_2d2,…,Pattern_2di,…,Pattern_2dN}；

(6)对第i个二维模式Pattern_2di计算该模式与训练好的字典中的M个类的聚类中心的距离D_c＝{d₁,d₂,…d_i,…,d_M}，找到D_c中的最小值D_min；

(7)在D_min对应的类中进行匹配，假设该类中有p个二维与其对应三维的字典原子，计算Pattern_2di与该类中的所有字典原子的灰度二维模式的距离D_bg＝{d_bg1,d_bg2,…d_bgi,…,d_bgp}，将Pattern_2di所在的二维图像块二值化，得到Pattern_2dbwi，计算Pattern_2di与该类中的所有字典原子的二值二维模式的距离D_bbw＝{d_bbw1,d_bbw2,…d_bbwi,…,d_bbwp}，将D_bg和D_bbw加权得到二维模式距离D_b＝aD_bg+bD_bbw；

(8)在字典原子中找到满足条件的块，对D_b进行升序排列，距离最小的块对应的n×n×n灰度三维模式即为匹配的三维图像块，依次完成对整个三维岩心块的重建。

2.根据权利要求1所述的基于超维的灰度岩心三维重建方法，其特征在于：

所述步骤(1)中，采用通过CT扫描获取的真实岩心灰度三维图，并对其进行二值化得到与之对应的二值岩心图像；

所述步骤(2)中，建立字典时，需要对真实岩心图像进行模式提取，然后按照预定准则保存二维模式及其对应的三维模式信息；设模板尺寸为n，I(x,y,z)表示真实岩心三维结构在点(x,y,z)的灰度值，用n×n的模板对I(x,y,z₀)采样，n×n×n的模板对I(x,y,z₀+n)采样，若采样中心为(x₀,y₀,z₀),则得到二维模式Pattern2d及其对应的三维模式Pattern3d：

Pattern2d＝I(x₀±n/2,y₀±n/2,z₀)

Pattern3d＝I(x₀±n/2,y₀±n/2,z₀±n)

为了保证模式集的连续性与差异性，扫描时中心点位置每次移动n-1，两个块之间只有一个面相互重叠即相邻块之间重叠一个面，扫描时第一个中心点位置为(x₀,y₀,z₀)，则下一个中心点的位置为(x₀+n-1,y₀,z₀)；

Patternset＝{(Pattern2d₁,Pattern3d₁),…(Pattern2d_i,Pattern3d_i),…(Pattern2d_N,Pattern3d_N)}

在岩心图像中，岩石颗粒和孔隙存在明显的灰度差异，在边缘存在灰度跳变，提取梯度信息可以反映结构信息，对二维图像提取一阶梯度和二阶梯度特征，将一阶梯度同二阶梯度一起作为灰度二维图像特征；

将features作为灰度二维模式，其对应的灰度三维图像块作为灰度三维模式，得到灰度模式集Patternsetgray，对二值化后的三维子图的底面提取二值二维模式，对应的二值的三维图像块作为二值三维模式，得到二值模式集Patternsetbw；

所述步骤(3)中，为了提高重建效率，采用字典分类的字典训练方法，利用Kmeans聚类算法，对二维模式Pattern2d进行聚类，将模式集划分为不同的类，相同的类的模式具有相似性，将模式集划分到类θ^M中的公式为：

θ^M(Patternset_i)＝ψ(Patternset_i(Pattern2d_i))i＝1,2,…N

式中ψ表示聚类算法，M表示类别数；通过计算样本与聚类中心的距离，将样本划分到与聚类中心距离最近的类中，在灰度超维重建的字典训练中，通过聚类分析，有助于重建时快速定位待重建的二维模式属于哪一类，然后在该类中去搜索相似的模式，若聚类个数设为m，则将整个字典按照聚类中心分为m个子空间，将二维及其对应的三维模式划分到归属的类，作为该子空间的字典原子；

所述步骤(4)中，将Patternsetgray中的灰度二维模式及其对应的灰度三维模式，Patternsetbw二值二维模式及其二值三维模式按照所属类别对应保存，得到根据真实岩心所建立的字典；

所述步骤(5)中，提取待重建图像的二维灰度模式时，为了增加块与块之间的连续性，防止重建时产生块效应，相邻的模式之间重叠一行或者一列得到{Pattern_2d1,Pattern_2d2,…,Pattern_2di,…,Pattern_2dN}；

所述步骤(6)中，为了提高重建效率，对于待重建图像的每一个二维灰度模式，首先在训练好的字典中搜索与该模式距离最近的类；

所述步骤(7)中，对于在步骤(6)中找到的与当前模式距离最近的类中，如果有p个二维与其对应三维的字典原子，计算Pattern_2di与该类中的所有字典原子的灰度二维模式的距离D_bg＝{d_bg1,d_bg2,…d_bgi,…,d_bgp}，将Pattern_2di所在的二维图像块二值化，得到Pattern_2dbwi，计算Pattern_2di与该类中的所有字典原子的二值二维模式的距离D_bbw＝{d_bbw1,d_bbw2,…d_bbwi,…,d_bbwp}，将D_bg和D_bbw加权得到二维模式距离D_b＝aD_bg+bD_bbw；使用欧式距离来进行不同模式间相似性估算，两个n维向量x＝(x₁,…,x_n)和y＝(y₁,…,y_n)之间的欧式距离为：

重建时以欧式距离来度量二维模式的相似度，距离越小，表示越相似；

所述步骤(8)中，在字典原子中找到满足条件的块，具体的对于不同的位置，采用如下方法：

i.当Pattern_2di是待重建图像的第一个提取的模式时，对D_b进行升序排列，距离最小的块对应的n×n×n灰度三维模式即为匹配的三维图像块；

ii.当Pattern_2di是待重建图像的除去第一个的前n行提取的模式时，对该块进行重建时，通过计算字典中所有原子的n×n×n三维模式与其左邻域已重建好的n×n×n三维图像块的距离D_lg＝{d_lg1,d_lg2,…d_lgi,…,d_lgp}，计算字典原子中二值三维模式与将左邻域的块二值化后的二值三维图像块的距离D_lbw＝{d_lbw1,d_lbw2,…d_lbwi,…,d_lbwp}，将D_lg和D_lbw加权得到左邻域距离D_l＝cD_lg+dD_lbw，将左邻域距离与二维模式距离加权得到新的距离D，D＝αD_b+βD_l，α+β＝1，对D进行升序排列，距离最小的块对应的n×n×n灰度三维模式即为匹配的三维图像块；

iii.当Pattern_2di是待重建图像的除去第一个的前n列提取的模式时，对该块进行重建时，通过计算字典中所有原子的n×n×n三维模式与其后邻域已重建好的n×n×n三维图像块的距离D_u＝{d_u1,d_u2,…d_ui,…,d_up}，计算字典原子中二值三维模式与将后邻域的块二值化后的二值三维图像块的距离D_ubw＝{d_ubw1,d_ubw2,…d_ubwi,…,d_ubwp}，将D_ug和D_ubw加权得到后邻域距离D_u＝eD_ug+fD_ubw，将后邻域距离与二维模式距离加权得到新的距离D，D＝αD_b+γD_u，α+γ＝1，对D进行升序排列，距离最小的块对应的n×n×n灰度三维模式即为匹配的三维图像块；

iv.当Pattern_2di是待重建图像的其余部分提取的二维模式时，对该块进行重建时，通过计算字典中所有原子的n×n×n三维模式与其左邻域已重建好的n×n×n三维图像块的距离D_lg＝{d_lg1,d_lg2,…d_lgi,…,d_lgp}，计算字典原子中二值三维模式与将左邻域的块二值化后的二值三维图像块的距离D_lbw＝{d_lbw1,d_lbw2,…d_lbwi,…,d_lbwp}，将D_lg和D_lbw加权得到左邻域距离D_l＝cD_lg+dD_lbw，计算字典中所有原子的n×n×n三维模式与其后邻域已重建好的n×n×n三维图像块的距离D_u＝{d_u1,d_u2,…d_ui,…,d_up}，计算字典原子中二值三维模式与将后邻域的块二值化后的二值三维图像块的距离D_ubw＝{d_ubw1,d_ubw2,…d_ubwi,…,d_ubwp}，将D_ug和D_ubw加权得到后邻域距离D_u＝eD_ug+fD_ubw，最终与二维模式距离加权得到新的距离D，D＝αD_b+βD_l+γD_u，α+β+γ＝1，对D进行升序排列，距离最小的块对应的n×n×n灰度三维模式即为匹配的三维图像。

Images