随着信息技术的不断发展，通过医学影像信息由计算机辅助诊疗的相关技术也层出不穷。在医疗大数据的背景下，对于肺结节影像的获取也愈来愈容易，如何在海量的肺结节影像中提取有效的特征信息对于肺结节良恶性准确的诊断至关重要。深度学习是根据人脑结构建立的一种分层模型，其在处理海量数据方面优势明显，同时对于肺结节良恶性判断的准确性更高。深度学习可在不依靠手工特征的条件下对海量数据进行并行运算，通过网络训练自动获取数据特征并分析处理。目前深度学习方法主要面向可向量化的特征信息数据，其中在图像处理领域主要包括2维与3维两类目标。

对于面向2维图像信息的深度学习框架，按照其架构可细分为CNN(convolutional neural network)-like的分类/回归训练框架，ELM^[1] (extreme learning machine)-like训练框架与GAN(generative adversarial nets)-like的生成对抗网络框架。按照特征深度化模式可分为CNN-like的卷积特征分类模型与GcForest深度森林^[2]相关的去卷积特征分类模型。CNN-like的分类/回归训练框架通过卷积通道来获得局部特征，其中较为典型的方法包括：区域卷积神经网络(R-CNN)算法^[3]、Faster R-CNN算法^[4]和YOLO算法^[5]等；ELM-like训练框架是一种针对单隐含层前馈神经网络(SLFN)的算法，其能够随机产生输入层与隐含层的连接权值及隐含层神经元的阈值且在训练过程中不需调整，只需要对隐含层神经元数量进行设置就能获得唯一的最优解；GAN是由Goodfellow等人^[6]提出的一种根据对抗过程对生成式模型进行预测的新框架，其与变分自动编码(VAE)^[7]是像素生成中最具有代表性的方法。在实际应用中为满足多分类需求, GAN框架通常延伸采用深度卷积生成对抗网络模型(DCGAN)^[8]进行对抗生成与分类训练，其将原始GAN中的生成模型G和判别模型D替换为深度卷积层，除进行生成图像判断之外还提供生成图像的分类判断。对于VAE而言，最简单的自动编码器网络由3层神经网络构成，第1层是数据输入层，第2层的节点数通常比输入层少且第3层与第1层类似，层与层间互相全连接，由于该网络将输入“编码”变成一个隐藏代码，然后从隐藏表示中“译码”出来。通过简单的测量重构误差和反传网络参数能很好地对该网络进行训练。而VAE与自动编码器的不同点在于其隐藏代码由训练期间学习到的概率分布获得。而GcForest深度森林是一种去卷积运算的图像分类方法，使用随机森林替代卷积运算，通过多粒度扫描增加级联森林的特征，进而进行提取。然而它们具有共同的缺点：输入端要求必须是等长的数据。由于获得的肺结节影像大小不一，因此在输入过程中要对尺度大的影像进行压缩，尺度小的影像进行拉伸，这样无疑会影响特征提取质量, 从而影响最终的分类结果，虽然Lin等人^[9]提出了使用特征金字塔网络融合多层特征，根据不同的特征图给出了特征金字塔网络，由于特征金字塔具有尺度不变性，可有效地解决目标检测中的目标物体大小不一致的问题，但增加了预处理难度且只能应用在2维数据信息上，对于3维数据仍然无法实现。

而在面向3维重建影像的深度学习分类方法中，一般根据较为高剂量的放射线探测方式获得较为严谨的2维图像信息并通过双边插值方法获得3维图像信息。对特征提取方法而言，目前方法均在2维框架上加以扩展。其中，Ji等人^[10]提出了一种3D-CNN的方法提取视频信息，通过在CNN的池化层和卷积层上添加一个参数作为时间维度的信息，能将2维图像处理的方法向视频处理上扩展。2017年，Wu等人^[11]以2维人脸图像及其深度图作为输入端，同时构建两个深度卷积网络模型，将所提取的2维人脸图像及其深度图的高层抽象特征作为神经网络的输入端，并融合2维特征及其深度信息完成3维特征的提取。尽管上述方法在各自场景中均取得了较为突出的结果，但在肺CT中结节图像的分类上，由于肺结节大小不一及其生长位置不确定，对肺CT图像的拍摄角度有着较为严格的要求，而在实际CT拍摄过程中角度因素是不可控的，所以若坚持使用卷积通道特征，需要优先解决海量肺结节角度归一化的问题。由于结节形态各异所以在归一化中仍然存在着归一角度的二义性，因此通过一种球面空间网格构造获得结节表观信息，从而规避角度归一化问题。此外该方法还可以通过取环方法获得不同的球网格比例的纹理与形边特征向量作为分类特征。

根据之前的工作，已将重建后的肺结节特征成功地映射到单位球面上。这部分内容将在基础建模工作中进行简略阐述。而在之后的工作中，本文首先利用网格-LBP的方法计算肺结节表面纹理特征并通过构建环的方式将纹理特征向量化，然后在3维肺结节重建过程中通过重心到肺结节中心点的距离差值构建形状能量损耗，提取肺结节在重建过程中的形状能量损耗并向量化。最后将两个向量化特征输入到改进的深度森林训练框架中，对肺结节良恶性进行训练和测试以验证本文方法。

在之前的工作中，已经使用传统3维参数边界可变形^[12]方法，通过两个概率视觉外观模型控制肺结节在分割过程中的变化将其分割出来，然后利用Fang等人^[13]提出的算法对其完成三角形3维网格化，且可得到网格节点坐标，然后对所有节点坐标进行球坐标变换，最后利用最小二乘法拟合3维肺结节原模型的坐标得到扩展系数，通过扩展系数对模型重建，重建结果及5类肺结节重建相对误差平均值如图 1、图 2所示。

为了达到图 3中肺结节中心点到每个节点都是单位距离的效果，进行如下运算：设$b$代表容斥映射算法的迭代次数，$M$代表总网格节点数，$P_{b, i}$表示第$b$次迭代时网格节点$i(i=1, 2, …, M)$的坐标，设$J$是网格节点相邻的数目，$d_{b, ij}$是第$b$次迭代时$i$点与$j$点之间的欧氏距离(图 4(b))，且$i=1, 2, \cdots, M，j=1, 2, \cdots, J$。设第$b$次迭代时$i$点与$j$点之间的位移是 ${\Delta _{b, ij}} = {P_{b, j}} - {P_{b, i}}$。同时设每个表面中节点之间的位移分别受相容力和排斥力常数$C_{A, 1, }$, $C_{A, 2}$, $C_{R}$控制。

容斥映射算法中的相容阶段通过如下迭代公式对网格中节点$P_{i} (i=1, 2, …, M)$进行位置的调整，使得它作为周围所有节点的中心点

$ \begin{array}{l} {{P'}_{b, i}} = {P_{b, i}} + {C_{A, 1}}\sum\limits_{j = 1;j \ne i}^J {{\Delta _{b, ij}}d_{b, ij}^2 + } \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{C_{A, 2}}\frac{{{\Delta _{b, ij}}}}{{{d_{b, ij}}}} \end{array} $

(1)

随后的排斥阶段将全部的节点向外推动，全部节点都沿着球中心射线的方向向外推动，最终映射到球体上均匀地隔开并且膨胀整个网格。为了避免未移动的节点与已移动的节点间冲突，因此每一个节点$P_{i}(i=1, 2, …, M)$的位置在映射之前进行更新, 即

${{P''}_{b + 1, i}} = {{P'}_{b, i}} + \frac{{{C_R}}}{{2M}}\sum\limits_{j = 1;j \ne i}^M {\left( {\frac{{{\Delta _{b, ij}}}}{{{{\left| {{\Delta _{b, ij}}} \right|}^2}}}} \right)} $

(2)

容斥映射算法如下:

输入：3维肺结节网格表面；

输出：映射完成之后的单位球。

   For $i$=1→$M$

   相容过程：

   选择任何一个网格节点进行处理；

   根据式(1)更新网格节点的位置，并且使其作为相邻所有节点的中心位置。

   排斥过程：

   根据式(2)更新网格节点位置并避免已发生和未发生移动的节点碰撞。

   End

   映射结果如图 5所示。

在3维网格上构建类似于LBP的模式首先需要一种方案来构建中心环周围的小平面环，并以有序方式遍历它们。令$\mathit{\boldsymbol{Z = }}\left\langle {\mathit{\boldsymbol{M}}, \mathit{\boldsymbol{S}}} \right\rangle $表示模型表面的三角网格，其中$\mathit{\boldsymbol{M}}$和$\mathit{\boldsymbol{S}}$分别表示网格的节点和小平面集合。这里设定每个节点周围包括6个三角平面(规则情况)，如图 6(a)(凸轮廓)，将向轮廓外指向的小平面称作$S_{\text{out}}$面。与每个$S_{\text{out}}$面一对一相邻并且位于凸轮廓内部的一组小平面，称作$S_{\text{in}}$面(与相应的$S_{\text{out}}$面共享一个位于凸轮廓上的边)。假设$S_{\text{out}}$面最初是以圆形方式排列在整个轮廓上，给定初始排列，填补每对连续的$S_{\text{out}}$面之间的缝隙，即提取位于两个连续$S_{\text{out}}$面之间共享共同节点的相邻面(称作$\mathit{\boldsymbol{S}}_{\text{gap}}$面)，如图 6(b)。桥接算法是介绍建立$\mathit{\boldsymbol{S}}_{\text{gap}}$面的过程。通过“桥接”迭代填补两个连续$S_{\text{out}}$面之间间隙的过程产生以圆形方式排列的一组面(见图 6(c))。环面序列方式与$S_{\text{out}}$面最初序列方向一致(顺时针或逆时针)。构建环算法的过程通过迭代调用“桥接”过程获得环结构(称作有序环)。而在特殊情况下，$S_{\text{out}}$面包括3个与中心点相邻的小平面，所获得的环由12个有序的小平面(即3个$S_{\text{out}}$小平面，以及连接连续$S_{\text{out}}$小平面之间间隙的9个$\mathit{\boldsymbol{S}}_{\text{gap}}$小平面)组成，如图 6(d)所示。桥接算法如下：

输入：$S_{{\rm{out}}}^i$, $S_{{\rm{out}}}^{i + 1}$(两个共享一个节点的连续$S_{\text{out}}$面)；$S_{{\rm{in}}}^{i }$(与$S_{{\rm{out}}}^i$面共享一个边)。

输出：$\mathit{\boldsymbol{S}}_{{\rm{grp}}}^i$(一组连续的 ${\mathit{\boldsymbol{S}}_{{\rm{grp}}}}$面填补$S_{{\rm{out}}}^i$和$S_{{\rm{out}}}^{i + 1}$面之间的缝隙)。

   procedure桥接($S_{{\rm{out}}}^i$, $S_{{\rm{out}}}^{i + 1}$, $S_{{\rm{in}}}^{i }$)

   $\mathit{\boldsymbol{S}}_{{\rm{grp}}}^i$=[]

   $m$← < $S_{{\rm{out}}}^i$, $S_{{\rm{out}}}^{i + 1}$ > 所共享的节点

   $gs$←与$S_{{\rm{out}}}^i$面相邻且不同于$S_{{\rm{in}}}^i$面同时包含节点$m$的面

   $bef$←$S_{{\rm{out}}}^i$

   While $gs$≠$S_{{\rm{out}}}^{i + 1}$ do

   将$gs$面加入到$\mathit{\boldsymbol{S}}_{{\rm{grp}}}^i$面中

   $new\_{gs}$←与$gs$面相邻且不同于$bef$面同时包含节点$m$的面

   $bef$←$gs$

   $gs$← $new\_{gs}$

   end while

   return $\mathit{\boldsymbol{S}}_{{\rm{grp}}}^i$

   End

构建环算法如下:

输入：$S_{{\rm{out}}}^i$：$S_{{\rm{out}}}^1$, $S_{{\rm{out}}}^2$, …, $S_{{\rm{out}}}^n$(在凸轮廓上); $S_{{\rm{in}}}^i$：$S_{{\rm{in}}}^1$, $S_{{\rm{in}}}^2$, …, $S_{{\rm{in}}}^n$(与 ${S_{{\rm{out}}}}$面一对一相邻，位于凸轮廓分开的区域)。

输出：有序环$\mathit{\boldsymbol{RING}}$。

   procedure构建环($S_{{\rm{out}}}^i$, $S_{{\rm{in}}}^i$)

   $\mathit{\boldsymbol{RING}}$=[]

   For all < $S_{{\rm{out}}}^i$, $S_{{\rm{out}}}^{i\% n + 1}$ > , $i=1, \cdots, n$ do

   将$S_{{\rm{out}}}^i$面加入到$\mathit{\boldsymbol{RING}}$中

   $\mathit{\boldsymbol{S}}_{{\rm{grp}}}^i$←桥接($S_{{\rm{out}}}^i$, $S_{{\rm{out}}}^{i + 1}$, $S_{{\rm{in}}}^i$)

   将$\mathit{\boldsymbol{S}}_{{\rm{grp}}}^i$面加入到$\mathit{\boldsymbol{RING}}$中

   End for

   Return $\mathit{\boldsymbol{RING}}$

   End

令$h\left( f \right):M \to \gamma $是在网格$M$上定义的标量函数^[14]。使用构建环算法获得环的顺序可让我们计算出它的二进制模式，从而与标准的LBP^[15]计算局部二进制的模式相同。通过对构建环算法中，由12个小平面所构成有序环的邻域进行阈值处理来获得中心小平面$f_{\text{c}}$处的网格LBP算子

$ \begin{array}{l} ML({f_{\rm{c}}}) = \sum\limits_{k = 0}^{11} {s(h\left( {{f_k}} \right) - h({f_{\rm{c}}})) \cdot a\left( k \right)} \\ \;\;\;\;\;\;\;\;\;\;\;s\left( x \right) = \left\{ {\begin{array}{*{20}{c}} 1&{x \ge 0}\\ 0&{x < 0} \end{array}} \right. \end{array} $

(3)

式中，$a(k)$是加权函数。二进制模式随着函数$a(k)$定义的不同而改变。本文采用Ojala等人^[16]提出的$a(k)$=$2^k$的传统LBP计算方式。

在映射后的单位球面上任取一个三角网格应用桥接算法和构建环算法，并根据周围12个邻域进行计算网格-LBP值并赋给中心三角网格，将单位球面上所有三角网格均进行相同的操作，最终每个网格都由一个新的数字表示，即肺结节的纹理信息。

通过对肺结节进行三角网格化得到$M$个节点和$F$个网格并且可得网格上每个节点的坐标，将肺结节中心点设置为原点坐标。在重建过程中，当重建阶数$n$=1时，重建模型是一个球体。在1阶重建模型中提取任何一个三角网格，设其节点坐标分别为 ${M_{1, 1}} = \left( {{x_{1, 1}}, {y_{1, 1}}, {z_{1, 1}}} \right)$，${M_{1, 2}} = \left( {{x_{1, 2}}, {y_{1, 2}}, {z_{1, 2}}} \right)$，${M_{1, 3}} = \left( {{x_{1, 3}}, {y_{1, 3}}, {z_{1, 3}}} \right)$(${M_{m, i}}$代表$m$阶第$i$个点的坐标)。求出三角网格重心坐标${G_{1, 1}} = (x_{1, 1}^{\rm{g}}, y_{1, 1}^{\rm{g}}, z_{1, 1}^{\rm{g}})$，即

$ \left\{ \begin{array}{l} x_{1, 1}^{\rm{g}} = \frac{{{x_{1, 1}} + {x_{1, 2}} + {x_{1, 3}}}}{3}\\ y_{1, 1}^{\rm{g}} = \frac{{{y_{1, 1}} + {y_{1, 2}} + {y_{1, 3}}}}{3}\\ z_{1, 1}^{\rm{g}} = \frac{{{z_{1, 1}} + {z_{1, 2}} + {z_{1, 3}}}}{3} \end{array} \right. $

(4)

可得到全部三角网格的重心坐标。

由重心坐标可求出重心到原点之间的距离${R_{1, 1}}$(${R_{m, f}}$代表$m$阶第$f$个网格到中心点的距离)

${R_{1,1}} = \sqrt {{{(x_{1,1}^{\rm{g}})}^2} + {{(y_{1,1}^{\rm{g}})}^2} + {{(z_{1,1}^{\rm{g}})}^2}} $

(5)

经过30阶重建后可得全部三角网格的重心坐标${G_{30, f}} = \left( {x_{30, f}^{\rm{g}}, y_{30, f}^{\rm{g}}, z_{30, f}^{\rm{g}}} \right)$。将肺结节重建过程中1阶重建模型的每个三角网格的重心坐标与其30阶重建模型中每个三角网格的重心坐标做差，可得到肺结节形变能量信息$E$，即形状信息

$ E = \Delta R = \sum\limits_{i = 1}^f {\left( {{R_{1, i}} - {R_{30, i}}} \right)} $

(6)

深度森林模型^[2]是一种去卷积运算的图像分类方法，使用多粒度扫描随机深林替代卷积运算，通过多粒度扫描增加级联森林的特征，级联森林中每一级都是一种决策树的集成方法。

多粒度扫描(multi-grain scanning)^[17]是采用滑动窗口扫描最初输入的样本，输入样本的尺度会随滑动窗口尺度的不同而变化。由于得到的样本均是基于曲面上的网格，因此本文提出了一种基于网格的多粒度扫描方法。将上述构建环的算法向多尺度框架扩展，构成多尺度同心有序环。从第1个环中提取与$\mathit{\boldsymbol{S}}_{\text{gap}}$小平面一对一相邻的小平面序列(图 7(a))，该序列延续了$\mathit{\boldsymbol{S}}_{\text{gap}}$面中的序列方向，通过填充这个序列中每两个面之间的间隙(图 7(b))能够得到一个新的环。迭代这个过程，可以得到多尺度同心有序环(图 7(c))。多尺度同心有序环的构建过程如下:

输入：${\mathit{\boldsymbol{S}}_{{\rm{out\_rt}}}}$(最初有序$S_{\text{out}}$面集合)；${\mathit{\boldsymbol{S}}_{{\rm{in\_rt}}}}$(最初有序$\mathit{\boldsymbol{S}}_{\text{in}}$面集合，与$S_{\text{out}}$面一对一相邻)；$Mr$(围绕$S_{{\rm{in\_rt}}}$面构建环的数目)。

输出：$\mathit{\boldsymbol{RINGS}}$(多尺度同心有序环)。

procedure多尺度同心有序环${\mathit{\boldsymbol{S}}_{{\rm{out\_rt}}}}$, ${\mathit{\boldsymbol{S}}_{{\rm{in\_rt}}}}$, $Mr$)

$\mathit{\boldsymbol{RINGS}}$=[]

${S_{{\rm{out}}}}$←${\mathit{\boldsymbol{S}}_{{\rm{out\_rt}}}}$

${S_{{\rm{in}}}}$←${\mathit{\boldsymbol{S}}_{{\rm{in\_rt}}}}$

For $i$=1，$Mr$ do

($\mathit{\boldsymbol{RING}}$, $New{S_{{\rm{out}}}}$, $\mathit{\boldsymbol{S}}\rm{gap}$)←构建环(${S_{{\rm{out}}}}$, ${S_{{\rm{in}}}}$)

将$\mathit{\boldsymbol{RING}}$加入到$\mathit{\boldsymbol{RINGS}}$中

${S_{{\rm{out}}}}$←$New{S_{{\rm{out}}}}$

${S_{{\rm{in}}}}$←$\mathit{\boldsymbol{S}}_{\rm{gap}}$

End for

Return $\mathit{\boldsymbol{RINGS}}$

End

在本文中环数为1、2、3、4、5这5个等级设为滑动窗口尺度。图 8所示是5环移动实例。由图 8可看出, 1环中有12个网格，2环中有36个网格，3、4、5环中分别有72、120、180个网格，并从中心网格点出发每次移动一个网格点。

级联森林(cascade forest)由多级结构组成^[18]，级联中每级(level)均含有完全随机森林和随机森林两种森林类型。每个完全随机森林中含有1 000棵完全随机树，每棵树随机选取数据特征作为划分节点的父节点。同样，每个随机森林中也含有1 000棵树，但每棵树在划分节点的父节点时，首先通过随机选取$\sqrt D$个特征($ D $是总特征数量)，然后选择一个基尼系数值最大的数作为父节点对节点划分。本文将肺结节分为5个类别，每棵树对输入的$x$(最初分类情况)进行分类时，都会得到一个属于每一个分类的概率。5分类，即产生一个5维向量，对所有树产生的向量取其平均，可得到一个森林中生成的5维向量，如图 9所示。将同一级中的众多森林所产生的向量进行连接，就形成了新$x$，作为下一级的输入。在级联的最后，将众多森林产生的向量取平均值，就产生了一个总的5维向量，哪个概率大，$x$就属于哪一类。图 10是本次实验中的级联结构图，由于本文将纹理和形状信息同时加入到框架中，因此1环总维数变为24，并设置步长和滑动窗口均为12，因此可以得到2个12维向量，再将所有向量经过两种森林进行处理，可产生2个5维(5分类)向量，2~5环计算方式相同。

本文已经利用上述方法将3维肺结节的纹理和形状特征进行向量化建模，现在采用基于网格深度森林的训练框架进行分类。本次实验肺结节数据主要由辽宁中医药大学附属医院的审查委员会(IRB)提供，审批通知编号为：2017111CS(KT)-040-01。数据集中根据肺结节的生长情况分成了独立站位型肺结节、磨玻璃影型肺结节、血管黏连型肺结节和胸腔黏连型肺结节4种(如图 11)。并且又按照肺结节良恶性划分成了5个级别：1——良性，2——疑似良性，3——不明，4——疑似恶性，5——恶性。提取7 326个3维肺结节作为本次实验的数据集，并且按照8：2比例将数据分成训练集和测试集进行验证，即训练集中共有5 861个肺结节，测试集中共有1 465个肺结节。具体数据分布如表 1所示。

本次实验中将肺结节的1、2等级判别为良性，3、4、5等级判定为恶性之后进行实验，这里引入4个衡量指标，分别是准确率(ACC)、特异度(SPE)、灵敏度(SEN)和受试者工作特征曲线下的面积(AUC)值。ACC代表分类模型对正、反例正确分类的比例；真阴性率(SPE, 也称特导度)代表分类模型对反例正确分类的比例；真阳性率(SEN，也称敏感度)代表分类模型对正例正确分类的比例；AUC代表曲线ROC(receiver operating characteristic curve)下方的面积。ROC曲线能使分类结果动态呈现出来，它的纵坐标表示SEN，横坐标表示假阳性率(1-特异度)。为了能够对ROC曲线精确地表达，可以使用AUC值进行表示，AUC取值范围通常在0.51之间，值越大代表分类能力越强。图 12展示了采用不同方法的AUC值比较；表 2展示了本文方法对测试数据集中不同类别的肺结节的分类；表 3展示了不同方法中ACC、SPE和SEN的比较。

对以上数据分析可知，本文方法的SEN达到88.46%，属于次优，但是SPE和ACC两个指标达到了最优，而且AUC的值是0.84，也达到了次优，因此所提方法对肺结节的识别力较高。

本文在基于容斥映射算法将重建肺结节映射到单位球面的基础上，构建了肺结节在球面上的纹理和形状特征并进行向量化建模。提出了一种基于网格向量化特征构造的方法，一方面用网格-LBP计算球表面纹理特征并通过构建环的方法将球面纹理特征进行向量化。另一方面在3维肺结节重建过程中构建形状能量损耗，提取肺结节在重建过程中形状的特征并向量化。最后在基于改进深度森林的训练框架对纹理与形状信息进行训练的过程中，提出了一种构建多尺度同心有序环的扫描方法代替传统的多粒度扫描方法完成了肺结节良恶性判断，结果显示SEN达到了88.46%，SPE和ACC分别是69.46%和76.06%，综合来看对肺结节的良恶性识别效果比较理想, 说明本文方法具有可行性。本文方法不仅考虑了数据的预处理过程，而且利用两个肺结节特征对其良恶性进行检测，相比于传统方法对一个特征进行检测的准确率要高，同时为3维模型中特征的提取及向量化提供了一个有效的方法。尽管如此，目前该方法的准确度受限于医生标定结节与周边组织边界的精确度，一旦在序列中单幅CT图像上存在误差则直接影响球面纹理晶格取值，故在未来的工作中需着重考虑如何使用更为有效的图像分割预处理手段从而增加球面纹理表达的准确性。