CN105006018A - 三维ct岩心图像超分辨率重建方法 - Google Patents

Abstract

本发明公开了一种三维CT岩心图像超分辨率重建方法，包括对CT扫描岩心三维图像进行退化，对退化后的三维图像进行分解、提取特征，得到特征块2；对提取特征应用算法进行字典训练，得到低分辨字典D_l及稀疏系数α；对输入的CT三维图像进行分解、提取特征，得到特征块1；根据特征块1及稀疏系数α求得输入的CT扫描岩心三维图像的高分辨率字典D_h；根据特征块1及低分辨率字典D_l求得新稀疏系数β；根据高分辨率字典D_h及新稀疏系数β，求得高分辨率三维图像块；用插值算法将CT扫描岩心三维图像进行上采样，得到放大的三维图像；将高分辨率三维图像块填充到放大的三维图像中，得到放大后的岩心试样三维模型。本发明解决了三维CT图像分辨率与试样尺寸相矛盾的问题。

Description

三维CT岩心图像超分辨率重建方法

技术领域

本发明涉及CT超分辨率图像重建技术领域，特别涉及一种三维CT岩心图像超分辨率重建方法。

背景技术

在石油地质科学中，具有储存和渗滤流体能力的岩石称为储集岩。储集岩孔隙性的好坏直接决定岩层储存油气的数量，渗透性的好坏控制了储集岩层内所含油气的产能。因此，在孔隙级别上研究岩石的渗流特性，对油气开采具有重要意义。而建立数字岩心则是开展上述研究工作的基础和关键。在开展石油地质分析研究工作中，CT扫描成像技术(Computed Tomography，电子计算机X射线断层扫描技术)得到了越来越广泛的应用，为了分析演示岩石孔隙和颗粒结构，得到岩层连通性、渗透率和形状等微观信息，通常采用对岩石进行CT扫描得到二维序列图，再根据二维序列图进行相应目标(孔隙和颗粒)分割提取，然后重建岩石的三维立体模型，对石油地质进行分析。

采用电子计算机X射线断层扫描成像技术对石油地质进行分析，为了得到准确的石油地质情况，需要获取高分辨率的岩心试样的孔隙结构图像。在采用CT扫描岩心试样时，为了使得到的电子计算机X射线断层扫描成像具有更高的分辨率，通常需要缩小CT扫描视场，甚至在原有岩心上取小样进行CT扫描。但这种以牺牲扫描式样尺寸获取到的高分辨率三维岩心图像，由于其尺寸的限制，只能得到小尺寸岩心试样的CT扫描三维图像，图像所显示的岩心试样孔隙结构，只能是岩心试样局部的孔隙结构，所获得的岩心特性不能表明岩心试样的整体特性，缺乏全局代表性，很难由此对石油地质情况作出准确分析。因此，为了得到岩心试样全面的孔隙结构，准确地分析石油地质情况，需要获取大尺寸岩心试样的整体结构图像。要获取大尺寸岩心试样的整体结构的图像，则需要对大尺寸的岩心试样进行扫描，但基于CT扫描原理和现有技术的图像扫描设备自身的限制，对大尺寸岩心试样扫描的三维CT图像分辨率会相应的降低。因此，现有技术的岩心试样CT三维图像存在图像分辨率与岩心试样尺寸相矛盾的问题，即小尺寸岩心试样的CT扫描三维图像分辨率高，但缺乏全局代表性；大尺寸岩心试样三维CT图像整体性好，但三维图像分辨率低。因此，如何解决这一问题，是所属领域科技工作者面临急待解决的问题。为此，国家设立专门的 基金项目，以求解决该问题。本发明正是受国家自然科学基金项目《岩石微观非均质结构三维图像重建及分辨率提升技术研究》(61372174)资助完成的。

发明内容

针对现有技术的岩心试样CT三维图像重建技术不足，本发明的目的旨在提供一种新的三维CT岩心图像超分辨率重建方法，以解决现有技术的岩心试样CT三维图像重建存在的图像分辨率与岩心试样尺寸相矛盾的问题，提升岩心试样CT三维图像的分辨率。

本发明为了解决分辨率和试样尺寸之间的矛盾，提出了利用超分辨率算法对三维CT图像进行超分辨率重建的方法，在不缩小岩石试样尺寸的情况下，利用现有CT图像进行超分辨率重建，达到分辨率提升的效果，得到更清晰的岩心结构。

本发明的上述目的可通过下述三维CT岩心图像超分辨率重建方法来实现，其技术方案的构成主要包括以下步骤：

(1)对输入计算机的CT扫描岩心试样三维图像进行下采样，并对下采样后的三维图像进行高斯模糊，得到退化后三维图像；

(2)对退化后的三维图像进行分解，提取分解后小区块特征，得到特征块2；

(3)根据稀疏表示原理，对提取到的特征块2应用正交匹配跟踪算法及KSVD算法进行字典训练得到低分辨字典D_l及稀疏系数α；

(4)对输入的CT扫描岩心试样三维图像进行分解，提取分解后的小区块特征，得到特征块1；

(5)根据所述特征块1及所述稀疏系数α，求得输入的CT扫描岩心试样三维图像的高分辨率字典D_h；

(6)根据所述特征块1及所述低分辨率字典D_l，求得新稀疏系数β；

(7)根据所述高分辨率字典D_h及所述新稀疏系数β，求得高分辨率三维图像块；

(8)用插值算法将CT扫描岩心试样的三维图像进行上采样，得到放大的三维图像；

(9)将求得的高分辨率三维图像块填充到放大的三维图像中，以恢复重建结构中缺失的特征信息，得到放大后的岩心试样三维模型。

本发明进一步的三维CT岩心图像超分辨率重建方法，与上述方法的不同之处在于：用小尺寸岩样的高分辨率三维图像来获取超分辨率重建的先验知识，建立相应的高分辨率字典和低分辨率字典，将该高低分辨率字典应用于同类型岩心的大尺寸岩样的低分辨率三维岩心图像，以提高大尺寸岩样低分辨率三维岩心图像的分辨率，其具体技术方案的构成主要包括以下步骤：

(1)对输入计算机的CT扫描小尺寸岩心试样三维图像进行下采样，并对下采样后的三维图像进行高斯模糊，得到退化后三维图像；

(2)对退化后的三维图像进行分解，提取分解后小区块特征，得到特征块2，

(3)根据稀疏表示原理，对步骤(2)提取到的特征应用正交匹配跟踪算法及KSVD算法进行字典训练得到低分辨字典D_l及稀疏系数α；

(4)对输入的小尺寸岩心试样CT扫描三维图像进行分解，提取分解后的小区块特征，得到特征块1；

(5)根据所述特征块1及所述稀疏系数α，求得CT扫描的小尺寸岩心试样三维图像的高分辨率字典D_h；

(6)对输入的与小尺寸岩心试样同类型的大尺寸岩心试样的CT扫描三维图像进行分解，提取分解后的小区块特征，得到大尺寸岩心试样的特征块1；

(7)根据所述大尺寸岩心试样的特征块1及所述低分辨率字典D_l，求得新稀疏系数β；

(8)根据所述高分辨率字典D_h及所述新稀疏系数β，求得大尺寸岩心试样的高分辨率三维图像块；

(9)用插值算法将大尺寸岩心试样的CT扫描三维图像进行上采样，得到其放大的三维图像；

(10)将所述高分辨率三维图像块填充到所谓放大的三维图像中，以恢复重建结构中缺失的特征信息，得到大尺寸岩心试样放大后的三维模型。

在本发明的上述技术方案中，为了取得更好的效果，最好在所述高分辨率三维图像块填充到放大的三维图像中后，使用反向投影算法，对恢复后的重建结构进行灰度修正，得到放大后的三维模型。

在本发明的上述技术方案中，三维图像的分解优先选用按照2x2x2或3x3x3或4x4x4像素大小进行分解。

在本发明的上述技术方案中，提取分解后的小区块特征，包括小区块的边缘信息和高频分量信息等特征。所述特征块2的特征提取，最好采用一阶、二阶梯度滤波、联合拉普拉斯模板和三维拉普拉斯高斯滤波器进行提取。所述特征块1的特征提取，最好先将模型下采样，然后进行三次插值放大得到低频分量模型，再用原始结构信息减去低频分量信息，得到特征块1。

在本发明的上述技术方案中，最好采用三次插值算法将CT扫描岩心试样的三维图像进行下采样和上采样，得到缩小和放大的三维图像。

本发明提出的三维CT岩心图像超分辨率重建方法的原理：

根据通用的成像模型，观测的低分辨率图像序列是理想清晰图像的变形、模糊、下采样和噪声污染等的退化版本。上述过程可以通过如下数学模型进行描述：

Y＝DBMX+n

其中X为理想高分辨率清晰模型，M为形变过程，B为模糊过程，D为下采样过程，n为污染性噪声。由于本发明的方法采用同一输入样本进行下采样进行重建，所以我们认为该模型可简化为：

Y＝DBX

即低分辨率模型可以高分辨率模型经过模糊和下采样得到。同理，通过其反过程，我们可以由低分辨率模型重建高分辨率模型。

稀疏性是信号普遍存在的一种特性。在信号的表示中，信号常常由大量疑似因素中的少量因素决定。一种信号的稀疏性表示与表示手段和度量方式有着密不可分的联系。传统信号的表示建立在正交基上，比如傅里叶变换，小波变换等。如果某种信号能在一组特定的基上稀疏表示，则可以用少量的稀疏表示系数来描述这个信号。

我们给定一个集合D＝{d,r＝1,2,3…K}，它的每一个元素是张成希尔伯特空间H＝R^N的单位矢量，其中K＞＞N，我们把D称为是超完备字典(以后简称字典)，它的每一列素称为一个原子。对于任意给定的信号X，我们希望在字典中自适应的选择一定数量的原子对信号做逼近，即信号表示为字典原子的一组线性组合(如图1)：

X＝Dα

其中，α为稀疏系数矩阵，D为超完备字典，信号X可由超完备字典D和稀疏 系数矩阵α表示。

整个重建过程即转化为求解超完备字典D_h，D_l和稀疏系数的问题。其过程如图2及图3所示。

本发明进一步方案的原理与基本方案基本一致。不同之处在于在图2中训练超完备字典D_h和D_l时，采用由CT扫描的与输入的CT扫描岩心试样三维图像具有相同特征的尺寸较小的高分辨率三维图像进行，在获取训练超完备字典D_h和D_l后，再采用与基本方案同样的步骤进行放大(如图3所示)。由于高分辨率三维图像能对岩心样本的特征进行更精细的描述，所以在图3的重建过程中，在一定程度上能更精确的放大并还原原始岩心样本的高频细节，达到更好的重建效果。

两种方案各有其优点及适用范围：基本方案在只有单一岩心样本的CT扫描岩心试样三维图像时便可对其进行超分辨率重建；而进一步的方案能对即有CT扫描岩心试样三维图像，又有与其同特性的高分辨率三维图像时，对CT扫描岩心试样三维图像进行精确的重建，能达到更好的重建效果。

本发明将二维超分辨率重建技术扩展并应用到三维CT岩心，突破了现有技术的岩心试样CT扫描设备对视域和分辨率的限制，在相同视域情况下，提升了三维模型图像的分辨率，使岩心试样三维结构的孔隙结构特征更为清晰明显，更有利于对岩心模型的后期特性分析。

表1为超分辨率重建PSNR结果对比

实验组	待放大三维图像大小(像素)	三次插值的PSNR	本发明的PSNR
				第一组	246×246×246	27.65	31.82/32.40
第二组	250×250×250	25.50	29.96/30.75
				第三组	296×296×296	27.56	32.02/32.73

附图说明

图1是信号的字典原子表示示意图；

图2是三维CT岩心高、低分辨率字典及稀疏系数求解过程示意图；

图3是三维CT岩心超分辨率重建示意图；

图4-1是本发明实施例中原始的真实岩心CT序列图像之一；

图4-2是本发明实施例中三次插值后重建的CT序列图像之一；

图4-3是本发明实施例中采用方案一重建的CT序列图像之一；

图4-4是本发明实施例中采用方案二重建的CT序列图像之一；

图4-5是本发明实施例中原始的真实岩心CT序列图像之一的切片图；

图4-6是本发明实施例中三次插值后重建的CT序列图像之一的切片图；

图4-7是本发明实施例中采用方案一重建的CT序列图像之一的切片图；

图4-8是本发明实施例中采用方案二重建的CT序列图像之一的切片图；

图5-1是本发明实施例中原始的真实岩心CT序列图像之二；

图5-2是本发明实施例中三次插值后重建的CT序列图像之二；

图5-3是本发明实施例中采用方案一重建的CT序列图像之二；

图5-4是本发明实施例中采用方案二重建的CT序列图像之二；

图5-5是本发明实施例中原始的真实岩心CT序列图像之二的切片图；

图5-6是本发明实施例中三次插值后重建的CT序列图像之二的切片图；

图5-7是本发明实施例中采用方案一重建的CT序列图像之二的切片图；

图5-8是本发明实施例中采用方案二重建的CT序列图像之二的切片图；

图6-1是本发明实施例中原始的真实岩心CT序列图像之三；

图6-2是本发明实施例中三次插值后重建的CT序列图像之三；

图6-3是本发明实施例中采用方案一重建的CT序列图像之三；

图6-4是本发明实施例中采用方案二重建的CT序列图像之三；

图6-5是本发明实施例中原始的真实岩心CT序列图像之三的切片图；

图6-6是本发明实施例中三次插值后重建的CT序列图像之三的切片图；

图6-7是本发明实施例中采用方案一重建的CT序列图像之三的切片图；

图6-8是本发明实施例中采用方案二重建的CT序列图像之三的切片图；

具体实施方式

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

实施例1

为了使本发明所述方法更加便于理解和接近于真实应用，下面对整个流程进详细整体的说明，其中包括本发明的核心方法：

(1)采用三次插值对扫描图像进行下采样，并对下采样后的结果进行高斯模糊，得到退化后的三维图像；

(2)对退化后三维图像按照2x2x2像素大小进行分解

(3)对分解后小区块采用一阶、二阶梯度滤波、联合拉普拉斯模板和三维拉普拉斯高斯滤波器对其进行特征提取，以得到其边缘特征信息和高频分量信息，得到特征块2；

(4)在字典训练过程中，根据稀疏表示原理，应用正交匹配跟踪(OMP)算法对训练对像进行稀疏分解，以特征块2作为初始字典D_l1，将退化后三维图像信号分别投影到字典D_l1中的原子上，并在每次投影分解的过程中，使投影残值与前面每个分量正交，得到矩阵，然后再利用KSVD算法，更新字典D_l1和矩阵，并保证信号的残差ε＝X-D_l1α和系数α的非零项均为最小：

$\{D_{11},\mathrm{\α}\}=\underset{D,\mathrm{\α}}{\min }\mathrm{\λ}{\left|\right|\mathrm{\α}\left|\right|}_p+\left|\right|X-D_{11}\mathrm{\α}|{|}_2^2$

通过以下逼近问题：

$\begin{array}{c}\begin{array}{cc}\{d,g\}=\underset{d,g}{\min }{\left|\right|E-{\mathrm{dg}}^T\left|\right|}_2^2& s.t.{\left|\right|d\left|\right|}_2=1\end{array}\\ E=X_I-\underset{i\≠j}{\mathrm{\Σ}}{d}_j{\mathrm{\α}}_{i,I}\end{array}$

得到最终超完备字典D_l和稀疏系数α，其中E为去掉第j个原子的残差矩阵，

d为字典中升级的原子，g^T为稀疏系数矩阵中新的一行。

(5)对输入的CT扫描的岩心试样三维图像按照2x2x2像素大小进行分解，

(6)提取扫描图像的边缘特征信息和高频分量信息特征，得到特征块1，即 此处提取特征的方法为，先将模型下采样，然后进行三次插值放大得到低频分量模型，再用原始结构信息减去低频分量信息，得到特征块1；

(7)特征块1可以由稀疏系数α和高分辨率字典D_h进行表示：

$F_h{\left\{p_h^i\right\}}_i=D_h\mathrm{\α}$

通过其伪逆计算求解得到高分辨率字典D_h，其中α⁺是α的伪逆：

$D_h=F_h{\left\{p_h^i\right\}}_i{\mathrm{\α}}^{+}=F_h{\left\{p_h^i\right\}}_i{\mathrm{\α}}^T{\left({\mathrm{\α\α}}^T\right)}^{-1}$

(8)根据扫描图像特征及低分辨率退化图像字典D_l，求得新稀疏系数β；

(9)根据高分辨率字典D_h及稀疏系数β，求得高分辨率下的三维结构块；

p_hh＝D_hβ

(10)用三次插值算法将扫描图像进行放大；

(11)将高分辨率的三维结构块填充到放大后的三维模型中，以恢复重建结构中缺失的高频部分信息；

(12)最后使用反向投影算法，对重建后的高分辨率结构进行灰度修正，最终得到放大后的三维模型。

上述步骤中，步骤(3)～(11)为本发明所述重建方法的主要步骤。

本实施例中，我们将真实的岩心CT序列图像人为地进行下采样，模拟低分辨率岩心结构，并对其进行三维超分辨率重建，并重建结果与原始真实岩心CT序列图像进行直观地对比观察，从而验证了本发明的可靠性。

实施例2

(1)使用小尺寸岩样进行高分辨率CT扫描，获得小岩样的高分辨率三维岩心图像I_h；然后对此高分辨率三维岩心图像进行下采样与模糊处理，使用与实施例1中建立高分辨率与低分辨率字典相同的方法，即实施例1步骤(1)-(7)，建立相应的高分辨率字典D_h和低分辨率字典D_l。该高分辨字典D_h和低分辨率字典D_l就可作为用于同类型岩心的低分辨三维图像的超分辨率重建的先验知识。

(2)将同类型岩心的低分辨率三维岩心图像作为待重建图像，对其进行与实施例1中步骤(3)相同的特征提取；然后使用与实施例1中重建高分辨率三维图像相同的重建方法，即是实施例1中的步骤(8)-(11)，最终实现对同类型岩心三维CT低分辨率图像的超分辨率重建。

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

Claims (10)

1.一种三维CT岩心图像超分辨率重建方法，其特征在于主要包括以下步骤：

(1)对输入计算机的CT扫描岩心试样三维图像进行下采样，并对下采样后的三维图像进行高斯模糊，得到退化后三维图像；

(2)对退化后的三维图像进行分解，提取分解后小区块特征，得到特征块2；

(3)根据稀疏表示原理，对提取到的特征块2应用正交匹配跟踪算法及KSVD算法进行字典训练得到低分辨字典D_l及稀疏系数α；

(4)对输入的CT扫描岩心试样三维图像进行分解，提取分解后的小区块特征，得到特征块1；

(5)根据所述特征块1及所述稀疏系数α，求得输入的CT扫描岩心试样三维图像的高分辨率字典D_h；

(6)根据所述特征块1及所述低分辨率字典D_l，求得新稀疏系数β；

(7)根据所述高分辨率字典D_h及所述新稀疏系数β，求得高分辨率三维图像块；

(8)用插值算法将CT扫描岩心试样的三维图像进行上采样，得到放大的三维图像；

(9)将求得的高分辨率三维图像块填充到放大的三维图像中，以恢复重建结构中缺失的特征信息，得到放大后的岩心试样三维模型。

2.一种三维CT岩心图像超分辨率重建方法，其特征在于包括以下步骤：

(1)对输入计算机的CT扫描小尺寸岩心试样三维图像进行下采样，并对下采样后的三维图像进行高斯模糊，得到退化后三维图像；

(2)对退化后的三维图像进行分解，提取分解后小区块特征，得到特征块2，

(3)根据稀疏表示原理，对步骤(2)提取到的特征块2应用正交匹配跟踪算法及KSVD算法进行字典训练得到低分辨字典D_l及稀疏系数α；

(4)对输入的小尺寸岩心试样CT扫描三维图像进行分解，提取分解后的小区块特征，得到特征块1；

(5)根据所述特征块1及所述稀疏系数α，求得CT扫描的小尺寸岩心试样三维图像的高分辨率字典D_h；

(6)对输入的与小尺寸岩心试样同类型的大尺寸岩心试样的CT扫描三维图像进行分解，提取分解后的小区块特征，得到大尺寸岩心试样的特征块1；

(7)根据所述大尺寸岩心试样的特征块1及所述低分辨率字典D_l，求得新稀疏系数β；

(8)根据所述高分辨率字典D_h及所述新稀疏系数β，求得大尺寸岩心试样的高分辨率三维图像块；

(9)用插值算法将大尺寸岩心试样的CT扫描三维图像进行上采样，得到其放大的三维图像；

(10)将所述高分辨率三维图像块填充到所谓放大的三维图像中，以恢复重建结构中缺失的特征信息，得到大尺寸岩心试样放大后的三维模型。

3.根据权利要求1或2所述的三维CT岩心图像超分辨率重建方法，其特征在于，所述高分辨率三维图像块填充到所谓放大的三维图像中后，使用反向投影算法，对恢复重建结构进行灰度修正，得到放大后的三维模型。

4.根据权利要求3所述的三维CT岩心图像超分辨率重建方法，其特征在于，对三维图像按照2x2x2或3x3x3或4x4x4像素大小进行分解。

5.根据权利要求1至4之一所述的三维CT岩心图像超分辨率重建方法，其特征在于，分解后的小区块提取特征包括小区块的边缘信息和高频分量信息。

6.根据权利要求5所述的三维CT岩心图像超分辨率重建方法，其特征在于，所述特征块2的提取方法，采用一阶、二阶梯度滤波、联合拉普拉斯模板和三维拉普拉斯高斯滤波器进行提取。

7.根据权利要求5所述的三维CT岩心图像超分辨率重建方法，其特征在于，所述特征块1的特征提取，先将模型下采样，然后进行三次插值放大得到低频分量模型，再用原始结构信息减去低频分量信息，得到特征块1，即。

8.根据权利要求5所述的三维CT岩心图像超分辨率重建方法，其特征在于，采用三次插值算法将CT扫描岩心试样的三维图像进行下采样和上采样。

9.根据权利要求6所述的三维CT岩心图像超分辨率重建方法，其特征在于，采用三次插值算法将CT扫描岩心试样的三维图像进行下采样和上采样。

10.根据权利要求7所述的三维CT岩心图像超分辨率重建方法，其特征在于，采用三次插值算法将CT扫描岩心试样的三维图像进行下采样和上采样。