本发明涉及一种图像修复增强方法，尤其涉及一种岩心三维图像超分辨方法，属于三维图像超分辨率重建技术领域。
背景技术：
在石油地质分析过程中，通常需要获得三维岩心结构来研究岩石的物理特性，分析岩石储集层中的渗流性质和运移规律。CT(Computed Tomography，计算机断层扫描)是一种三维成像技术，可以提供高质量的物体内部结构图像，被广泛应用于医学图像诊断和地质勘探中。但是在生产生活中CT图像分辨率受到扫描设备和实际情况的制约。通常情况下，想要获得更高精度的CT图像，需要耗费更长的扫描时间，由于扫描时间延长，不可避免地造成CT图像在层间方向上连续性变差和分辨率降低，同时会带来噪声影响。同时，在地质勘探领域中设置高分辨率进行CT扫描，不仅会增加成本，而且随着分辨率升高，会导致可视域FOV(Field ofView)减小。超分辨率是一种常用的提升图像分辨率技术方法，单张图像超分辨率本身是一个病态问题，当高分辨率(High Resolution，HR)图像经过降质得到低分辨率(Low Resolution，LR)图像，其中部分的高频信息已经丢失，属于一个不可逆转的过程。传统的超分辨方法有两种，一种是通过图像内部相关信息完成重建，例如基于邻域嵌入的方法，但是这种方法无法从内部图像中获取更多的缺失信息，因而效果有限；另一种是基于样本学习的方法，通过外部样本学习到LR与HR特征之间的映射关系，在重建过程中通过学习到的先验信息来约束解空间，从而完成重建，此类方法通常具有更好的重建效果。
目前对于三维图像超分辨的研究还较少，大部分还是针对二维图像的超分辨重建。而在三维图像的研究中，主要是针对于医学MRI(Magnetic Resonance Imaging)、CT图像小范围信息缺失的三维图像修复以及视频的修复。神经网络在特征提取、图像分割，目标检测等任务上取得了很好的效果，表现出了很强的泛化能力。Dong C,Chen C L,He K,et al.(Image Super-Resolution Using Deep Convolutional Networks[J].IEEE Transactions on Pattern Analysis&Machine Intelligence 2016,38(2):295-307)将卷积神经网络引入到单张二维图像超分辨重建中，但是该方法只适用于二维图像超分辨率，不能完成对三维卷积操作。而对于岩心三维图像，其在各个方向上的相关性都一样，需同等考虑各个方向的信息，需要在三个维度上提升图像分辨率。利用图像超分辨率技术提升CT图像的分辨率是一种方便有效的方法，可以为下一步的地质研究和图像处理提供更清晰的样本数据，有着重要的研究价值。
技术实现要素：
本发明的目的就是提供一种方便有效的三维图像超分辨重建算法，将深度学习技术引入到三维CT图像超分辨工作中来，提出了一种新的网络结构3DSRCNN(three dimensional convolutional neural network，三维超分辨卷积神经网络)，能够完成针对三维图像的超分辨率重建任务。
本发明通过以下技术方案来实现上述目的：
1、本发明所述岩心三维图像超分辨重建方法包括以下步骤：
(1)将训练集LR图像送入本方法提出的网络，其中，第一层网络负责提取LR图像中的低频特征信息；
(2)第二层至第十一层网络负责学习LR图像中低频特征到HR图像中高频特征之间的映射关系；
(3)第十二层网络利用中间层学习到的映射关系，将原始LR图像中低分辨特征映射为高分辨特征，使用残差的方式最终完成重建，生成HR图像，作为最后一层网络的输出；
(4)HR图像作为基准，使用残差学习的方法计算均方根误差(Mean Square Errors，MSE)损失率，利用动量梯度下降的方法加速训练，不断迭代训练更新各网络层中三维卷积核中的参数；
(5)在训练的过程中，使用自适应学习率，梯度裁剪的方法优化训练过程，使用(1)～(5)步骤中的网络配置不断迭代训练网络参数，直至训练时计算的MSE损失率趋于收敛，停止训练；
(6)使用训练后保存的网络模型完成超分辨重建，对LR图像分块，将所有子块送入网络，经过网络输出后得到HR图像，然后将所有输出的HR子块拼接成原始尺寸大小HR图像。计算重建图像与原始图像PSNR和SSIM值，用于评价重建效果。
上述方法的基本原理如下：
单张图像超分辨率是由于缺乏细节信息而导致的不确定性问题。假设给定真实高分辨率图像为X，降质后的图像为Y，超分辨过程定义为：
X≈F(Y)
超分辨率过程相当于找到一个函数F(Y)，其能够在一定程度上恢复LR图像至HR图像。在三维超分辨率过程中，需要同时对三个维度同时进行处理。由于三维图像数据的计算复杂度和数据量较大，可能会导致训练过程不稳定，因此在数据集预处理时需要对原始三维CT图像进行分块处理，同时采用动量梯度下降优化策略来加速训练。在通常情况下，CT扫描图像不易获取，训练集数量有限，需要尽可能使用较小的训练样本集来进行超分辨网络的训练，且岩石CT图像具有对比度低、纹理单一、孔隙结构复杂等特点，图像中蕴含的高频信息较多，可以使用较深层的网络来更好地提取特征信息来解决这个问题。当网络层数过深会导致训练时出现梯度爆炸的现象，因此使用残差学习和梯度裁剪来抑制梯度爆炸的问题。使用不同缩放倍率样本混合在一起的训练集训练网络，网络能够学习到不同倍率之间的映射关系，从而网络适用于不同倍率的超分辨重建场景。在使用训练好的网络模型完成超分辨重建时，为了解决程序过多占用计算机内存问题，同样需要对较大像素尺寸的图像进行分块送入网络重建，最后拼接网络生成的图像至原始大小像素尺寸。
具体地，所述步骤(1)中，网络共有12层三维卷积神经网络，每一层中使用的均是三维卷积操作，每一层网络均有64个三维卷积核，卷积核尺寸大小为3×3×3，在每次卷积的过程之前先进行填零操作，每一层卷积后通过ReLU(Rectified Linear Unit)激活函数作为当前网络层的输出，ReLU(X)＝max(0,WX+b)，其中X,W分别表示该层网络的输入和网络权重参数，卷积后的结果作为下一层的输入；
所述步骤(2)中，数据经过中间10层网络，每层网络结构相同，如步骤(1)中所述；
所述步骤(4)中，定义输入LR图像为x，HR图像作为标签为y，残差图像r＝y-x，基于均方根误差的损失函数L(Θ)可以写成：其中Θ表示网络参数，f(x)表示为网络的输出，使用动量SGD优化网络参数Θ，批样本图像数量为m,输入LR图像{x1,x2,…,xm},作为标签的HR图像{y1,y2,···,ym},残差图像{r1,r2,···,rm}，梯度更新按照如下规则：
v←ρ*v+g
Θ←Θ-lr*v
其中，ρ表示动量系数，v为速度参数，用来计算累积动量；
所述步骤(5)中，学习率按照如下公式更新：其中lr为训练时的学习率，epoch表示当前训练的次数，step是预设的学习率衰减步长，使用梯度裁剪可以将计算得到的梯度限制在预设的[-θ,θ]范围内，当训练时梯度g超过预设值时调整反复迭代训练网络，直至计算的MSE损失率趋于收敛，停止训练；
所述步骤(6)中，训练好的网络模型实质上为一系列矩阵参数，将LR图像送入网络模型后，网络的输出即为HR图像。在将LR图像送入网络之前，需要对其进行分割预处理，将大像素尺寸的图像分成子块再送入网络，完成所有子块重建后再拼接所有子块重建图像，从而获取完整HR图像，重建结果使用PSNR和SSIM来衡量重建质量，其中其中表示原始HR图像，MAXL为图像灰度级中最高值。给定两张图像x和y，两张图像的结构相似性可SSIM可以表示为：其中μx是x的平均值，μy是y的平均值，是x的方差，是y的方差，σxy是x和y的协方差，c1,c2是用来维持稳定的常数。
本发明的有益效果在于：
本发明可以使用训练好的网络模型对单个三维CT样本在三个方向进行分辨率的提升。相较之前基于学习的超分辨率方法，需要对不同降质倍率的样本进行训练相应的网络从而进行重建。本方法训练好的单个网络模型可以对不同降质倍率的低分辨图像重建，同时，使用GPU加速训练可以在较短的时间内训练出能性能较好地网络，能够节约训练时间，有更强的实用性。训练好的网络内部实质为一系列卷积核权重参数，在重建过程中只进行矩阵运算，解决了以往基于学习的方法在重建过程中迭代求解的方式，大大地减少了计算量，有更快的速度和效率。本发明中提出的方法相较之前的方法在重建质量上有更高的精度。
附图说明
图1是本发明实施例中CT图像构成方式示意图；
图2是本发明实施例中提出的网络结构示意图；
图3是本发明实施例中用于训练网络的原始岩心CT三维图像之一；
图4是本发明实施例中用于训练网络的原始岩心CT三维图像之二；
图5是本发明实施例中用于训练网络的原始岩心CT三维图像之三；
图6是本发明实施例中用于训练网络的原始岩心CT三维图像之四；
图7是本发明实施例中用于训练网络的原始岩心CT三维图像之五；
图8是本发明实施例中用于测试重建的低分辨率CT三维图像之一；
图9是本发明实施例中用于测试重建的低分辨率CT三维图像之二；
图10是本发明实施例中用于测试重建的低分辨率CT三维图像之三；
图11是本发明实施例中用于测试重建的低分辨率CT三维图像之四；
图12是本发明实施例中用于测试重建的低分辨率CT三维图像之五；
图13为经过下采样后通过差值恢复到原尺寸的低分辨率CT三维图像；
图14为图13中任选的某张二维切片图像；
图15为图14中框选位置放大图像；
图16为经过下采样后通过差值恢复到原尺寸的低分辨率CT三维图像；
图17为图16中任选的某张二维切片图像；
图18为图17中框选位置放大图像；
图19为本发明实例中网络训练次数与重建PSNR效果曲线图
具体实施方式
下面结合具体实施例和附图对本发明作进一步说明：
实施例：
为了使本发明所述修复方法更加便于理解和接近于真实应用，下面从原始岩心CT训练样本的预处理开始到CT图像超分辨重建完成为止进行整个流程的整体说明，其中包括本发明的三维超分辨重建网络的使用方法：
(1)使用CT机扫描切割后的小岩心样本，可以得到多张连续的二维图像序列，依次读入并存储生成三维图像。所述的CT三维图像如图1所示。从中抽取连续的400张图片，每单张二维图片裁取其中400*400像素大小的区域，按照这样的规则分别从不同类别的岩石样本中挑选10组，作为训练集和测试集，训练集图像如图3、4、5、6、7所示，测试集图像如图8、9、10、11、12所示。
(2)在实验之前需要对的原始三维CT图像样本进行预处理，每组CT样本集处理方法一致，设原始三维HR图像为{Y}，尺寸为400*400*400，首先对原始图像进行灰度值归一化，然后分别采用×N(N＝2,3,4)倍下采样降质，接着使用双三次插值算法×N倍插值恢复到原始尺寸400*400*400，作为LR图像{X}；对应的LR图像和HR图像相应配对组成为{X；Y}，分别训练阶段网络的输入和计算MSE的标签。从不同缩放倍率的样本集中按照1:1:1的比例选取样本组成混合倍率样本集作为最后的训练集。
(3)本发明提出的三维超分辨重建网络如图2所示，共包含12层三维卷积神经网络，每层网络包含64个三维卷积核，每个三维卷积核尺寸为3×3×3，在每次卷积之前都需要对当前输入的三维块进行填零操作，即对输入三维块各个边界添加一个灰度级为0的像素点。
(4)将原始HR图像Y和LR图像X分别在三个方向上进行分割处理，步长为stride，切分成尺寸为的子块；将X和Y的每一个子块分别对应，随机打乱顺序，制作训练样本集{X，Y}，生成的训练集中样本数量为：其中Iinp表示原始图像块单边尺寸，Isub表示分割之后图像的单边尺寸，此处，Isub设为25。
(5)训练的时候，梯度裁剪范围θ设置为0.4，批训练样本bathch_size为32,每次从训练集中选取32个样本进行训练，输入LR图像为x，HR图像作为标签为y，残差图像r＝y-x，使用计算批次均方根误差，然后初始学习率lr设为0.1，动量因子ρ设为0.9，使用动量梯度下降的方法对各层网络中的三维卷积核权重Θ进行更新，具体规则如下：
v←ρ*v+g
Θ←Θ-lr*v
学习率随着训练次数epoch的增加不断减少，按照如下公式更新：其中step设为10，是预设的学习率衰减步长。训练时使用梯度裁剪来抑制梯度爆炸，设置梯度门限θ＝0.4，在训练时梯度g超过[-0.4,0.4]时，自动调整梯度为
(6)使用(3)～(5)步骤中的网络参数设置，反复迭代训练样本，直至计算的MSE损失率趋于收敛，停止训练。大约需要训练20次左右，网络性能趋于收敛，网络性能随训练变化曲线如图19所示，如果MSE出现波动可以适当增加迭代训练次数。
(7)在使用训练好的网络模型重建之前，需要对待重建的LR图像进行分块处理，原始LR图像像素尺寸为400*400*400，将其分为64个像素值为100*100*100的子块。分别将子块送入训练好的网络，然后将所有子块的拼接成完整像素尺寸为400*400*400的HR图像。
(8)对重建结果与原始图像进行计算PSNR和SSIM，PSNR计算方式为：其中表示真实HR图像，f(x,y,z；Θ)网络重建图像，MAXL为图像灰度级中最高值，此处设为255。两张图像的结构相似性可SSIM可以表示为：其中μx是x的平均值，μy是y的平均值，是x的方差，是y的方差，σxy是x和y的协方差，c1,c2是用来维持稳定的常数。
上述步骤中，步骤(2)～(8)为本发明所述三维图像超分辨重建方法的主要步骤。
本实施例中，我们使用不同倍率下采样真实CT图像得到LR图像，模拟对三维图像造成的信息缺失，作为存在信息缺失的LR三维图像进行修复，并将其修复结果与真实岩心CT序列图像进行PSNR和SSIM的对比，从而验证了本发明的有效性，通过对比低分辨率图13、14、15和重建图16、17、18，可以从人眼视觉上发现重建结果具备超分辨率效果。在同样的实验条件下，通过使用相同的测试图像与之前方法对比，本方法在不同缩放倍率上重建精度更高，具体实验对比结果如表1所示：
表1本方法与其他方法实验测试结果对比
上述实施例只是本发明的较佳实施例，并不是对本发明技术方案的限制，只要是不经过创造性劳动即可在上述实施例的基础上实现的技术方案，均应视为落入本发明专利的权利保护范围内。