本发明创新性的提出了基于超维的灰度岩心三维重建方法，尤其涉及一种灰度岩心的三维重建方法，属于图像处理领域。
背景技术：
在石油地质研究中，岩心三维微观结构是研究岩心宏观物理特性的基础。三维数字岩心能够在孔隙尺寸级别上反映孔隙微观结构，计算岩心声学、电学特性和模拟渗流过程，是分析岩心微观物理特性的有力工具。然而实际工程中高精度岩心三维图像难以直接获取，利用高精度二维图像进行数字建模，可以有效重建岩心三维图像。
虽然利用成像设备获取的图像大多是灰度图像，但目前较多的三维重建算法是针对二值化图像。重建时将灰度图像进行二值化，然后重建二值三维图像。然而灰度图像比二值图像能反映更多信息，重建灰度图像后可以进行二值化，但重建的二值图像无法恢复成灰度图像。为了更多保留原始信息，对于灰度图像的三维重建算法的研究是很有必要的。
目前的灰度岩心重建算法主要有：tahmasebi等提出的基于互相关函数的模拟算法(ccsim，cross-correlation–basedsimulation)，该方法能够很好继承层与层之间的连续性，但是不易控制层与层之间的随机变化性。以及纹理合成算法，其基本思想是以小规模纹理作为样本，合成较大尺寸的结果纹理。
但是上述方法都没有考虑岩心的真实三维结构，超维重建是将已有的真实岩心图像作为先验信息，利用真实的三维图像信息指导重建过程。本发明提出了一种基于超维的灰度岩心三维重建方法，保证了重建后三维灰度岩心图像的统计特征与二维原始岩心图像保持一致，形态特征保持相似性。该研究项目受国家自然科学基金项目《岩石微观非均质结构三维图像重建及分辨率提升技术研究》(61372174)资助。
技术实现要素：
本发明的目的在于创新性的提出一种新的灰度岩心三维重建方法，并且保证重建后三维灰度岩心图像在统计特征上与二维原始岩心图像保持一致。
本发明通过以下技术方案来实现上述目的：
(1)将真实岩心灰度图像i二值化得到ibw，孔隙相为1，岩石相为0；设定模板尺寸为n，将真实岩心灰度图像i和二值图像ibw分成一一对应的h×w×n的三维子图；
(2)对每一个灰度的三维子图的底面提取特征，得到features；利用模式提取方法以光栅路径扫描所有三维子图，将features作为灰度二维模式，其对应的灰度三维图像块作为灰度三维模式，得到灰度模式集patternsetgray，对二值化后的三维子图的底面提取二值二维模式，对应的二值的三维图像块作为二值三维模式，得到二值模式集patternsetbw；
(3)设定类别数m，使用kmeans聚类算法对patternsetgray中的灰度二维模式进行聚类分析，将模式集分为m类；
(4)将patternsetgray中的灰度二维模式及其对应的灰度三维模式，patternsetbw二值二维模式及其二值三维模式按照所属类别对应保存，输出所建立的字典。
(5)对于待重建的灰度岩心二维图像，每次重叠一行或者一列提取待重建图像的二维模式{pattern2d1,pattern2d2,…,pattern2di,…,pattern2dn}；
(6)对第i个二维模式pattern2di计算该模式与训练好的字典中的m个类的聚类中心的距离dc＝{d1,d2,…di,…,dm}，找到dc中的最小值dmin；
(7)在dmin对应的类中进行匹配，假设该类中有p个二维与其对应三维的字典原子，计算pattern2di与该类中的所有字典原子的灰度二维模式的距离dbg＝{dbg1,dbg2,…dbgi,…,dbgp}，将pattern2di所在的二维图像块二值化，得到pattern2dbwi，计算pattern2di与该类中的所有字典原子的二值二维模式的距离dbbw＝{dbbw1,dbbw2,…dbbwi,…,dbbwp}，将dbg和dbbw加权得到二维模式距离db＝adbg+bdbbw；
(8)在字典原子中找到满足条件的块，对db进行升序排列，距离最小的块对应的n×n×n灰度三维模式即为匹配的三维图像块，依次完成对整个三维岩心块的重建。
上述方法的基本原理如下：
超维重建(sdr)是一种全新的三维重建方法，该方法主要借鉴了基于学习的超分辨率重建算法的思想，将已有的真实岩心三维图像作为先验信息，学习二维图像块到三维图像块的映射关系，建立对应的字典，在重建时，利用待重建图像的信息和在先验信息中匹配的部分进行重建，重建出来的三维结构在统计信息与形态结构与原始图像都比较接近。该算法合理地运用了已有的图像资源，将成像设备获取三维图像与数学模拟方法结合，解决了利用成像设备获得大量三维图像价格昂贵和传统三维重建算法在形态上与原始图像差距太大的缺点。本方法在超维重建的框架下提出了灰度岩心图像超维重建算法，实现了灰度岩心图像二维到三维的重建，如图1所示。
所述步骤(1)中，采用通过ct扫描获取的真实岩心灰度三维图，并对其进行二值化得到与之对应的二值岩心图像。
所述步骤(2)中，建立字典时，需要对真实岩心图像进行模式提取，然后按照一定准则保存二维模式及其对应的三维模式信息。设模板尺寸为n，i(x,y,z)表示真实岩心三维结构在点(x,y,z)的灰度值，用n×n的模板对i(x,y,z0)采样，n×n×n的模板对i(x,y,z0+n)采样，如图2所示为5×5的二维模板和5×5×5的三维模板在图像上采样，若采样中心为(x0,y0,z0),则得到二维模式pattern2d及其对应的三维模式pattern3d。
pattern2d＝i(x0±n/2,y0±n/2,z0)(1)
pattern3d＝i(x0±n/2,y0±n/2,z0±n)
固定模板尺寸对整个真实岩心三维结构按照光栅路径扫描，即可获得灰度图像二维及其三维结构的模式集patternset，n表示模板对的个数。为了保证模式集的连续性与差异性，扫描时中心点位置每次移动n-1，两个块之间只有一个面相互重叠即相邻块之间重叠一个面，比如扫描时第一个中心点位置为(x0,y0,z0)，则下一个中心点的位置为(x0+n-1,y0,z0)。
在进行重建时，更加关注结构信息，若直接用图像的灰度值作为模式信息，匹配条件苛刻，容易出现在字典中无法找到匹配块的情况。岩心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(patternseti)＝ψ(patternseti(pattern2di))i＝1,2,…n…n(8)
式中ψ表示聚类算法，m表示类别数。kmeans算法计算简单，使用广泛，可提前预设聚类的类别数，通过计算样本与聚类中心的距离，将样本划分到与聚类中心距离最近的类中。在灰度超维重建的字典训练中，通过聚类分析，有助于重建时快速定位待重建的二维模式属于哪一类，然后在该类中去搜索相似的模式。若聚类个数设为m，则将整个字典按照聚类中心分为m个子空间，将二维及其对应的三维模式划分到归属的类，作为该子空间的字典原子，字典的结构示意如图5所示。
所述步骤(4)中，将patternsetgray中的灰度二维模式及其对应的灰度三维模式，patternsetbw二值二维模式及其二值三维模式按照所属类别对应保存，得到根据真实岩心所建立的字典。
所述步骤(5)中，提取待重建图像的二维灰度模式时，为了增加块与块之间的连续性，能够有效防止重建时产生块效应，所以相邻的模式之间重叠一行或者一列得到{pattern2d1,pattern2d2,…,pattern2di,…,pattern2dn}；
所述步骤(6)中，为了提高重建效率，对于待重建图像的每一个二维灰度模式，首先在训练好的字典中搜索与该模式距离最近的类。
所述步骤(7)中，对于在步骤(6)中找到的与当前模式距离最近的类中，如果有p个二维与其对应三维的字典原子，计算pattern2di与该类中的所有字典原子的灰度二维模式的距离dbg＝{dbg1,dbg2,…dbgi,…,dbgp}，将pattern2di所在的二维图像块二值化，得到pattern2dbwi，计算pattern2di与该类中的所有字典原子的二值二维模式的距离dbbw＝{dbbw1,dbbw2,…dbbwi,…,dbbwp}，将dbg和dbbw加权得到二维模式距离db＝adbg+bdbbw。在本方法中，我们使用欧式距离来进行不同模式间相似性估算，两个n维向量x＝(x1,…,xn)和y＝(y1,…,yn)之间的欧式距离为：
重建时以欧式距离来度量二维模式的相似度，距离越小，表示越相似。
在步骤(8)中，在字典原子中找到满足条件的块，具体的对于不同的位置，采用如下方法：
i.当pattern2di是待重建图像的第一个提取的模式时，对db进行升序排列，距离最小的块对应的n×n×n灰度三维模式即为匹配的三维图像块；
ii.当pattern2di是待重建图像的前n行提取的模式时(除去第一个)，对该块进行重建时，通过计算字典中所有原子的n×n×n三维模式与其左邻域已重建好的n×n×n三维图像块的距离dlg＝{dlg1,dlg2,…dlgi,…,dlgp}，计算字典原子中二值三维模式与将左邻域的块二值化后的二值三维图像块的距离dlbw＝{dlbw1,dlbw2,…dlbwi,…,dlbwp}，将dlg和dlbw加权得到左邻域距离dl＝cdlg+ddlbw，将左邻域距离与二维模式距离加权得到新的距离d，d＝αdb+βdl，α+β＝1，对d进行升序排列，距离最小的块对应的n×n×n灰度三维模式即为匹配的三维图像块；
iii.当pattern2di是待重建图像的前n列提取的模式时(除去第一个)，对该块进行重建时，通过计算字典中所有原子的n×n×n三维模式与其后邻域已重建好的n×n×n三维图像块的距离du＝{du1,du2,…dui,…,dup}，计算字典原子中二值三维模式与将后邻域的块二值化后的二值三维图像块的距离dubw＝{dubw1,dubw2,…dubwi,…,dubwp}，将dug和dubw加权得到后邻域距离du＝edug+fdubw，将后邻域距离与二维模式距离加权得到新的距离d，d＝αdb+γdu，α+γ＝1，对d进行升序排列，距离最小的块对应的n×n×n灰度三维模式即为匹配的三维图像块；
iv.当pattern2di是待重建图像的其余部分提取的二维模式时，对该块进行重建时，通过计算字典中所有原子的n×n×n三维模式与其左邻域已重建好的n×n×n三维图像块的距离d1g＝d1g1，d1g2，…d1gi，…，d1gp}，计算字典原子中二值三维模式与将左邻域的块二值化后的二值三维图像块的距离dlbw＝{dlbw1，dlbw2，…dlbwi，…，dlbwp}，将dlg和dlbw加权得到左邻域距离dl＝cdlg+ddlbw，计算字典中所有原子的n×n×n三维模式与其后邻域已重建好的n×n×n三维图像块的距离du＝{du1，du2，…dui，…，dup}，计算字典原子中二值三维模式与将后邻域的块二值化后的二值三维图像块的距离dubw＝{dubw1，dubw2，…dubwi，…，dubwp}，将dug和dubw加权得到后邻域距离du＝edug+fdubw，最终与二维模式距离加权得到新的距离d，d＝αdb+βdl+γdu，α+β+γ＝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，得到二值化后的岩心图像ibw，设定模板尺寸为3，将真实岩心灰度图像i和二值图像ibw分成一一对应的600×600×3的三维子图；
(2)对每一个灰度的三维子图的底面提取特征，得到features；利用模式提取方法以光栅路径扫描所有三维子图，将features作为灰度二维模式，其对应的灰度三维图像块作为灰度三维模式，得到灰度模式集patternsetgray，对二值化后的三维子图的底面提取二值二维模式，对应的二值的三维图像块作为二值三维模式，得到二值模式集patternsetbw；
(3)设定类别数m，使用kmeans聚类算法对patternsetgray中的灰度二维模式进行聚类分析，将模式集分为m类；
(4)将patternsetgray中的灰度二维模式及其对应的灰度三维模式，patternsetbw二值二维模式及其二值三维模式按照所属类别对应保存，输出字典。
(5)使用二维灰度图像进行三维重建，图7-1为待重建灰度图像，图7-2为待重建二值图像。每次重叠一行或者一列提取待重建图像的二维模式{pattern2d1，pattern2d2，…，pattern2di，…，pattern2dn}；
(6)对第i个二维模式pattern2di计算该模式与训练好的字典中的m个类的聚类中心的距离dc＝{d1,d2,…di,…,dm}，找到dc中的最小值dmin；
(7)在dmin对应的类中进行匹配，假设该类中有p个二维与其对应三维的字典原子，计算pattern2di与该类中的所有字典原子的灰度二维模式的距离dbg＝{dbg1,dbg2,…dbgi,…,dbgp}，将pattern2di所在的二维图像块二值化，得到pattern2dbwi，计算pattern2di与该类中的所有字典原子的二值二维模式的距离dbbw＝{dbbw1,dbbw2,…dbbwi,…,dbbwp}，将dbg和dbbw加权得到二维模式距离db＝adbg+bdbbw；
(8)按照如下准则在字典原子中找到满足条件的块：
i.当pattern2di是待重建图像的第一个提取的模式时，对db进行升序排列，距离最小的块对应的3×3×3灰度三维模式即为匹配的三维图像块；
ii.当pattern2di是待重建图像的前n行提取的模式时(除去第一个)，对该块进行重建时，通过计算字典中所有原子的3×3×3三维模式与其左邻域已重建好的3×3×3三维图像块的距离dlg＝{dlg1,dlg2,…dlgi,…,dlgp}，计算字典原子中二值三维模式与将左邻域的块二值化后的二值三维图像块的距离dlbw＝{dlbw1,dlbw2,…dlbwi,…,dlbwp}，将dlg和dlbw加权得到左邻域距离dl＝cdlg+ddlbw，将左邻域距离与二维模式距离加权得到新的距离d，d＝αdb+βdl，α+β＝1，对d进行升序排列，距离最小的块对应的3×3×3灰度三维模式即为匹配的三维图像块；
iii.当pattern2di是待重建图像的前n列提取的模式时(除去第一个)，对该块进行重建时，通过计算字典中所有原子的3×3×3三维模式与其后邻域已重建好的n×n×n三维图像块的距离du＝{du1,du2,…dui,…,dup}，计算字典原子中二值三维模式与将后邻域的块二值化后的二值三维图像块的距离dubw＝{dubw1,dubw2,…dubwi,…,dubwp}，将dug和dubw加权得到后邻域距离du＝edug+fdubw，将后邻域距离与二维模式距离加权得到新的距离d，d＝αdb+γdu，α+γ＝1，对d进行升序排列，距离最小的块对应的3×3×3灰度三维模式即为匹配的三维图像块；
iv.当pattern2di是待重建图像的其余部分提取的二维模式时，对该块进行重建时，通过计算字典中所有原子的3×3×3三维模式与其左邻域已重建好的3×3×3三维图像块的距离dlg＝{dlg1,dlg2,…dlgi,…,dlgp}，计算字典原子中二值三维模式与将左邻域的块二值化后的二值三维图像块的距离dlbw＝{dlbw1,dlbw2,…dlbwi,…,dlbwp}，将dlg和dlbw加权得到左邻域距离dl＝cdlg+ddlbw，计算字典中所有原子的n×n×n三维模式与其后邻域已重建好的n×n×n三维图像块的距离du＝{du1,du2,…dui,…,dup}，计算字典原子中二值三维模式与将后邻域的块二值化后的二值三维图像块的距离dubw＝{dubw1,dubw2,…dubwi,…,dubwp}，将dug和dubw加权得到后邻域距离du＝edug+fdubw，最终与二维模式距离加权得到新的距离d，d＝αdb+βdl+γdu，α+β+γ＝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方向线性路径函数。从上述图中可以看到本文算法重建结果的两点概率函数和线性路径函数与原始三维结构、待重建图像相似，证明该算法重建结果的统计特征与真实灰度岩心三维图像相一致，能取得较好的重建效果，证明了本方法的有效性。
上述实施例只是本发明的较佳实施例，并不是对本发明技术方案的限制，只要是不经过创造性劳动即可在上述实施例的基础上实现的技术方案，均应视为落入本发明专利的权利保护范围内。