CRF 3D-UNet肺结节分割网络

侯腾璇，赵涓涓+，强 彦，王三虎，王 磐

(1.太原理工大学 信息与计算机学院，山西 晋中 030600；2.吕梁学院 计算机科学与技术系，山西 吕梁 033000)

0 引 言

近年来，随着计算机和胸外科技术的日趋成熟，医学影像在肺癌中起到不可忽视的作用。诊断肺癌的影像学方法包括了X射线检查、正电子发射断层扫描成像和计算机断层扫描成像(computed tomography，CT)等，医生无需对患者的病变组织进行侵入式提取便可观察到病变组织的大小、位置、形态、纹理等特征。肺癌最开始多表现为肺部结节，它们尺寸小、对比度低、形状异质化高，CT图像中肺结节的边缘模糊、灰度不均匀、受到噪声和伪影的影响大。用于肺结节分割的传统方法[1-3]大多数比较局限、结节种类针对性强，分割精度也远远低于医师期望的水平。而且，CT影像数据指数型增长，加重了医生的负担。因此，提出一种对于肺结节种类针对性小甚至无针对性的方法实现肺结节的分割对肺癌的早期诊断十分重要。机器学习的快速发展给肺结节的分割带来了新方向。Olaf Ronneberger等提出UNet[4]网络是一种像素级图像分割的深度卷积网络，比较适用于生物、医学影像领域。以LUNA16为数据集，通过UNet网络训练的迭代优化，分割准确率达到约88.91%，有很大的提升空间。因此，提高肺结节分割的准确度定位在对UNet网络结构的研究改进。

1 相关工作

1.1 UNet分割网络

肺结节的分割不是一个孤立的领域，而是一个从粗略到精细的自然过程。UNet网络结构特点是将分割目标的底层信息和高层信息结合用于目标像素点定位。其中，经过多层的卷积操作和下采样操作得到的低分辨率图像是原始图像对应的底层信息，该信息为分割目标提供在原始图像中的相对位置信息。而高层信息是指编码路径对应到相同高度的解码路径的高分辨率图像，该图像中含有的信息为目标图像的分割提供梯度等精细特征。但是，UNet结构仍存在很多缺陷，使用二维的卷积、池化操作提取肺部CT图像的特征导致肺结节丢失丰富的空间信息；再者，下采样过程中丢失了很多上下文信息，在上采样的过程中并没有完全恢复目标的细节和相应的空间维度，导致上采样的结果比较模糊，对图像中的细节不敏感；UNet网络实现了对肺部图像像素级的分割，但是分割的过程只提取单个像素的特征，完全忽略像素与像素的关联。总结以上分割方法存在的问题：像素级别的分割方法忽视了肺结节的空间一致性[5]，导致上下文信息的丢失，未将像素间的关系考虑在内，这些问题使得UNet分割的结果准确率低。

1.2 全连接条件随机场

Philipp Krahenbuhl等提出了全连接条件随机场[6](fully connected conditional random fields，fully CRF)是一种基于概率的无向图结构，可以用于像素级图像的分割。目标图像中的各个像素都存在对应的类别标签，视目标图像的像素为图的顶点，顶点作为状态特征，连接图中所有像素点作为边，边代表转移特征，求解像素标签时考虑图像中其余像素对该像素的影响，极大地细化了标记和分割，使得边界处分割准确。

1.3 CRF 3D-UNet网络

针对上述存在的问题，本文提出一种将3D-UNet网络与全连接条件随机场结合的模型，用于分割肺结节。本文的数据集使用Lung Nodule Analysis 2016 (LUNA2016)，这个数据集中含有的CT图像有对应的医生标记好的结节，并提供具体位置。第一步对LUNA16数据集中结节样例进行预处理，为3D-UNet网络生成适当的训练集，使用这些训练集来训练有监督分割器。3D-UNet网络分割结果根据输出的各像素是否属于结节的概率，构建CT影像的像素概率地图。通过机器视觉算法[8]，将各区域合并，输出疑似区域集。最后将此疑似区域集经过一系列处理操作然后输入到全连接条件随机场，实现从粗分割到细分割的突破。

2 方法描述

本文提出了CRF 3D-UNet网络结构，3D-UNet网络(如图1所示)对CT序列图像进行操作之后得到了像素概率地图，即粗分割的结果。之后连接全连接条件随机场(如图2所示)，该技术将分割的目标图像中像素间的关系考虑在内，输出具有空间一致性的分割结果图。本文提出的整体模型框架将全连接条件随机场与3D-UNet网络[9,10]技术结合，前端使用3D-UNet网络进行特征粗提取，后端使用全连接条件随机场技术优化前端的输出，端到端网络的构建使得分割图更加精确。

图1 3D-UNet网络结构

图2 全连接条件随机场结构(fully CRF)

2.1 3D-UNet网络结构

由于肺部CT图像是一个三维的断层图，肺结节在体积、形状及许多其它特征如精细度、内部结构、球形度等方面有很大的变化[11]，加之血管、支气管等复杂的构造使肺部CT图像具有丰富的上下文环境。使用二维的卷积、池化、上采样操作会丢失很大一部分空间和上下文信息。对比二维CNN，三维CNN将数据的空间信息利用起来，提炼出图像层与层之间隐含的代表性特征。所以本文采用三维CNN结构，而且为了保留边界处的卷积结果，使用Same填充，处理CT序列影像数据。

本文提到的3D CNN结构指的是3D卷积层、3D池化层、3D上采样层。输入肺结节的CT序列图像，经过3D CNN操作得到不同通道去除冗余之后的显著性特征。

3D卷积层：CRF 3D-UNet网络从3D卷积层开始，输入F和卷积核W之间的操作定义为

(1)

3D最大池化层：3D卷积层之后接着是3D最大池化层，利用不同通道之间的平移不变性对特征进行子采样。3D最大池化操作不做卷积操作，只取当前窗口最大值作为新图的像素值。然后深度上滑动，得到多个特征图。本文的最大池化层包含两种，第一种应用于图1网络结构的第6层和第12层，假设当前第l层是卷积层，而第l+1层是最大池化层，对于经过卷积层的输出结果选择一个立方体领域内的最大体素值进行激活，大小缩为原来的一半；第二种应用于网络结构的第3层和第9层，假设当前第l层是卷积层，而第l+1层是最大池化层，对于经过卷积层的输出结果选择一个立方体领域内的最大体素值进行激活并且将最大值的位置坐标记录下来，之后传递到对应的上采样层，最大池化的操作如图3所示。

图3 3D最大池化+标记坐标

3D上采样层：上采样就是池化的逆过程(索引在上采样过程中发挥作用[8])，在上采样层中可以得到在池化中相对卷积核的位置。经过3D卷积层和3D最大池化层之后得到了最低分辨率的特征地图，很明显经过池化之后，每个卷积核会丢失9个像素值，这些权重是无法复原的，因此，需要通过对特征地图进行3D上采样操作来弥补丢失的信息。如图1网络结构所示，此结构是对称的，第13层和第19层的上采样操作是先将立方体区域的体积翻倍然后对立方体领域进行三线性插值，第16层和第22层是接收第3层和第9层传递过来的最大值的坐标位置，然后将立方体区域的体积加倍，将最大值填入原来的位置，剩余立方体区域中的值补为0。同时这4层融合对应层传递过来的经过卷积之后的深度特征。

综上，网络结构前半部分经过卷积和池化操作损失了分割目标很多的上下文信息和空间信息，而实现分割的关键是像素点的精确定位，所以在池化阶段做出了改进，在最大池化的过程中引入索引功能，融合两种不同的池化操作可以实现像素更精准的定位。

2.2 精分割

后端使用全连接条件随机场技术，考虑输入图像像素间的关联性，实现对肺结节更加精确的分割。从网络结构图1和图2可以看出，前端的输出是经过3D-UNet网络级联浅层特征的结果O1，对结果图进行一系列操作之后输入到全连接条件随机场。这里的操作指的是先将结果图O1与输入图I进行数值型或操作得到图H1，然后在图H1中定位到O1中结节的位置将其像素值置为0，同时将胸腔区域的像素值也置为0，得到图H2；然后综合考虑本文用到的LUNA16数据集和山西某医院的数据集结节的尺寸情况，统计所得结节最大直径大约是27.442 mm。通过改进遗传算法的任意图形最大内接矩形[12]算法找到粗分割图中肺结节的最大内接矩形，然后定位矩形的4个坐标点，将坐标相加除以2找到最大内接矩形的中心点，最后以中心点为中心将H2框到28 mm×28 mm的区域内，将H2作为全连接条件随机场的输入，H2包括有可能成为结节的其它像素点。

对于相同大小的输入图H2，每个像素i具有类别标签Mi，这里的类别标签有两类：肺结节和非肺结节。这样每个像素点作为节点，像素之间的连线作为边，构成了完全无向图。序列M={M1,M2,…Mn}和T={T1,T2,…Tn},标签Y={Y1,Y2}即构成全连接条件随机场(T,M)，T的大小小于输入的肺结节序列图像大小，代表对应序列图像的真实标签，Mj是赋予每个像素点的分类标签。我们通过观测变量T来推测像素i对应的类别标签Mi。条件随机场符合吉布斯[13]分布如式(2)

(2)

其中，g=(v,e)代表序列图M中的节点v和边e，t是图M中的最大团，φt是最大团的势函数，Z(T)是规范化因子[14]，代表一系列最大团之和。吉布斯能量函数E(M|T)如式(3)

E(M|T)=∑t∈Tgφt(Mt|T)

(3)

为了简便，以下省略全局变量T，一元势函数φv和二元势函数φe组成目标能量函数，如式(4)

(4)

其中，φv(Mi)=-logP(Mi)。一元势函数原本的计算只考虑了单个像素点的特征就对像素点进行分类标签，这和前端的输出一致，所以本文中全连接条件随机场直接计算二元势函数即可。二元势函数将像素与像素之间的关系考虑进来，为相似的像素标记相同的标签，而对于差别很大的像素给予不同类型的标签，使像素点标签的分配更符合空间一致性。这样的操作才能使肺结节边界处的分割结果更准确。最终每个像素点对应的标签通过最大后验概率求得如式(5)

X*=argmaxM∈YNP(M|T)

(5)

由于经过3D-UNet网络的每一个像素的一元分类器的输出独立于其它像素的分类器的输出并且通过最大后验概率求得的一元分类的输出存在噪声且不连续，所以对二元势函数的训练很重要。二元势函数的形式如式(6)

(6)

其中，在肺结节的训练中g=2，需要训练两个将特征fi和fj作为自变量的高斯核，每个核赋予一定的权重ω，当Xi≠Xj时，u(Mi,Mj)=1；否则u(Mi,Mj)=0。全连接条件随机场的二元势函数的训练考虑像素点的手工特征如：位置信息和颜色信息，将其编码成向量融入二元势函数如式(7)

(7)

其中，pi代表位置变量，Ti代表颜色向量，θα,θβ,θγ是需要训练的高斯核参数，这些参数从数据中学习得来，控制相似性的程度。在式(7)中第一项表示位置相近而且有相似颜色的像素点成为相同标签的可能性，而第二项表示移除H2图中一些孤立的区域。对于输入图H2，其一元势函数的值已经通过3D-UNet网络得出，二元势函数通过编码手工特征进一步精确定位像素，将一元势函数和二元势函数相加得到每个像素点属于肺结节的概率。

3 实验结果

3.1 数据集和实验环境

本文方法所使用的实验数据是LUNA16数据集和山西某医院的数据集,其中将LUNA16作为训练集，将山西某医院的数据平均分开，一半作为训练集，一半作为测试集。LUNA16数据集包括888个肺部疾病患者，其中肺结节的最大直径是27.442 mm，并且数据集中选取的是至少由3位专家标注的1186个结节，将这些结节作为最后要检测的区域。采集到的山西某影像科数据含有219例病人，经诊断120 例病人是肺结节，这些病人平均含有148张CT数据，17 760张CT图像中含有2800个结节(有可能多个切片表示一个结节)，大直径范围为：1.6 mm-23 mm，肺结节的平均直径为5.1 mm。经一组专家诊断最终确定，2800个结节中良性结节有1180个，恶性结节有1620个。本文网络的运行环境为Python3.4，Keras框架，Theano后端，CentOS7.4，GPU tesla m40,处理器Intel(R) Xeon(R) CPU E5-2630 v4 @ 2.20 GHz。

3.2 网络结构训练过程

3.2.1 3D-UNet网络结构

首先是对图像的预处理过程，读取LUNA16和合作医院作为训练部分的数据的标注得到每个结节序列图对应的包含结节的最大范围的掩膜图；然后提取肺部感兴趣区域生成每个结节对应的肺实质区域以缩小我们的结节搜索范围；二值化图像设置合理的阈值实现能够区分肺和密度更高的组织，之后使用形态学操作填补黑色肺实质区域；将感兴趣区域的掩码作用于CT图像，对感兴趣区域的边界框进行裁剪，之后将得到的图像大小调整为512×512；最终结果是一系列可以作为训练样本的肺部图像。将肺部图像成对读入网络进行训练，随机初始化权重，将最初的学习速率设置为1.0e-5，循环次数设置为20，一次读入的图像数设置为8，使用Adam优化器，从头开始训练经过3个小时训练完毕。损失函数使用Dice[15]系数。

3.2.2 全连接条件随机场推断过程

将经过3D-UNet网络训练得到的2586张粗分割图作为全连接条件随机场的输入，训练二元势函数的参数。全连接条件随机场的推断过程[16,17]主要包括：根据3D-UNet的输出粗分割图初始化一元势函数；信息传递，对于当前像素点i，依据求得不同的高斯核函数来计算其余像素点对当前像素点的影响程度；若两个相近像素点的标签差异很大，添加惩罚项；将二元势函数的结果与一元势函数相加来更新当前像素点的标签。算法流程见表1。

表1 全连接条件随机场的推断过程

3.3 实验结果

3.3.1 评价标准

本文从3个度量标准来对分割结果做出评测，灵敏性、特异性、交叉比[18]。

2016年8月，新的北辰基督教堂建成启用。新教堂建筑面积12000平方米，最高处24米，阶梯式大厅可容纳2500人。

灵敏性是指所有的肺结节得到正确分割的概率；灵敏性=真阳性/(真阳性+假阴性)。特异性是指所有非结节得到正确诊断的概率；特异性=真阴性/(真阴性+假阳性)。交叉比是指本文方法分割出的结果与专家手工标记结果的交集与并集的比值。

将本文方法与之前方法做出比较，结果见表2，我们的方法准确度达到93.25%，目前效果最佳。

3.3.2 结果与分析

我们用山西某医院的219列病人的一半数据作为测试集验证我们方法的准确率。山西某医院的数据包括微小型肺结节、孤立型肺结节、血管粘连型肺结节、胸膜牵拉型肺结节和磨玻璃型肺结节。使用提出的方法对比了在去噪、处理细小组织、边缘模糊目标物体与背景灰度值相近目标物体等方面的分割效果。图4显示了该方法对胸膜牵拉型肺结节的分割的效果，可以明显看出使用本文提出的方法分割出肺结节边缘形状最符合专家手动标记的分割结果，胸膜和结节的连接部分比较细致，灰度值接近，边界区域较难分割，条件随机场技术针对像素点的位置与周围像素点对其产生的影像，本文使用两种算法的融合可以将胸膜和肺结节边界区分开。

表2 不同方法的准确率/%

图4 胸膜牵拉型肺结节的分割

图5显示了该方法对血管粘连型肺结节的分割效果，由于血管的复杂性，将小血管误认为肺结节分割很容易，导致准确分割肺结节困难。可以明显看出图5(c)、图5(e)、图5(f)、图5(g)的分割结果并没有将肺结节和粘连处的血管区分开，本文提出的方法达到了这个效果。

图5 血管粘连型肺结节的分割

图6显示了该方法对磨玻璃型肺结节的分割效果，磨玻璃型结节的特性(形状多变、颜色较淡等)和血管很难区分，很明显看出图6(c)、图6(e)、图6(f)、图6(g)方法对肺结节和血管的分割效果不够理想，而本文提出的方法效果较好，也更贴近专家提供的期望边界。

图6 磨玻璃型肺结节的分割

图7 孤立型肺结节的分割

图8显示对不同类型肺结节的分割结果，可以看出添加了全连接条件随机场的分割效果要整体优于没有添加该技术的方法所得结果，而且在某些毛刺特征明显的结节边缘处分割效果优于人工标记，主要是因为该方法兼顾了灰度对比度与边缘连续性。但是对于胸膜牵拉型的肺结节有时候会出现过分割的现象，主要是由于胸膜牵拉型肺结节的牵拉处灰度值和肺结节的灰度值相比对比度很低。

图8 对不同类型的肺结节进行分割的结果

表3显示了对山西省某医院不同类型肺结节分割的错误率E[18]对比，其错误率定义为式(8)

(8)

其中，R0和Rm分别是本文算法分割的肺结节边界和专家手动标记的肺结节边界。

表3 不同类型肺结节的错分率

3.3.3 与目前成熟分割算法的对比

成熟的分割算法有：模糊C均值聚类算法(FCM)，李越[19]提出将FCM算法作为基础，同时应用小波变换方法针对CT图像展开分解，之后将分解后的低频图的像素点作为FCM算法的基础点,然后采用马氏距离来进一步修正，从而确保更加准确反映医学图像中的信息。但是经实验此算法对于胸膜牵拉型结节和血管粘连型结节的分割效果比较差。Jonathan Long等[17]提出用于图像像素级分割的全卷积网络(fully convolutional networks，FCN)，FCN把原本卷积网络后面接的几个全连接层都换成卷积，这样就可以获得一张2维的feature map，而后接softmax层获得每个像素点的分类信息，从而解决了分割问题。但是此方法存在的问题就是下采样过程丢失的信息并未在反卷积过程弥补完整，所以造成分割结果准确度不高。Simon JégoU等[20]改进DenseNets来处理语义分割问题。但由于此方法主要应用于通用图像，如：分割一幅图像中的建筑物、车、蓝天等像素差别特别大，而且前景和背景没有引起太大的差异。对比肺结节CT图像，正好相反，目标和背景对比差异太大，在一个大的背景下分割小目标，像素对比差异较微小。所以将该方法用于肺部CT图像的分割效果并不佳。表4显示将本文方法与这3种算法用于肺结节分割的结果准确率对比。图9显示了这些算法和本文方法的结果图对比。

表4 不同方法对不同类型肺结节的分割准确率/%

图9 成熟算法与本文算法的结果图对比

4 结束语

由于肺结节CT图像含有丰富的空间信息，结合肺结节的形状、纹理异质性高[21]等特性，本文提出了CRF 3D-UNet网络结构用于肺结节的分割。该网络前端利用了CT图像中丰富的空间信息，将两种最大池化操作结合用于像素的精确定位，上采样的过程中融合浅层的深度特征来弥补下采样过程中丢失的上下文信息，输出肺结节的粗分割图；后端结合全连接条件随机场优化前端的输出，考虑了像素间的关联性，实现从粗分割到精分割的跨越，显著提高了分割的精度。经过对山西某医院数据的测试，最终本文方法获得93.25%的准确率，93.14%的灵敏性，90.21%的特异性。

使用CRF 3D-UNet网络方法分割出CT图像的肺结节之后，后续的计划是对分割出的肺结节进行恶性度分类，分析其后期的生长情况，为肺癌的早期预防提供真正意义的辅助诊断。