随着生物医学的进步，人们对不同疾病的理解逐渐深入，但癌症仍然是尚未攻克的难题^[1]。脑胶质瘤的临床治疗方案一般由治疗过程中肿瘤的分级情况决定。目前，脑胶质瘤分级标准严格基于病理学检查，这种方法有一定的局限性，包括侵入性、采样误差、解释的变化性等^[2]，而基于影像学检查的分级方法^[3]没有侵入性，逐步成为研究的热点。但是，早期医生往往只能通过肉眼来获取图像信息，耗时较长，自动化的分级方法则能有效缩短诊断与治疗时间。磁共振成像常被用于脑胶质瘤检查，它能够提供大量的脑部软组织解剖、功能、代谢等信息，有些与疾病相关的特征，肉眼很难发现，这也需要人们开发自动的工具或算法^[4,5]。因此越来越多的研究者参与到脑胶质瘤磁共振图像的自动化分级中。

       影像组学是近年来兴起的一种图像分析方法，它从医学图像中自动地高通量地提取大量定量特征。假设对于给定模态的医学图像数据，通过自动的或半自动的软件对它们进行定量分析，可能可以从图像中发掘出更多更好的信息^[6]，就能将影像组学与常规磁共振图像相结合，使其在脑胶质瘤分级的应用中发挥巨大作用。本研究用影像组学的方法对基于常规磁共振图像的脑胶质瘤进行了分级。

1　材料与方法

1.1　数据

1.1.1　脑胶质瘤分级

       根据2016年世界卫生组织(World Health Organization，WHO)中枢神经系统肿瘤分类标准，脑胶质瘤分类包括^[7]星形细胞肿瘤、少突胶质细胞肿瘤、室管膜肿瘤、脉络丛肿瘤、神经元及混合神经元神经胶质肿瘤等。脑胶质瘤可以分级为低评级(low grade glioma，LGG)和高评级(high grade glioma，HGG)。

1.1.2　数据集介绍

       本研究使用了2017脑胶质瘤分割竞赛数据集(brain tumor segmentation challenge 2017，BRATS2017)^[8,9]，该数据集提供了专家手工标记的肿瘤区域以及肿瘤的分级情况，其中提供肿瘤分级信息的有285例数据，LGG有75例，HGG有210例，数据模态为T1加权像、增强T1加权像、T2加权像和液体衰减反转恢复序列(fluid attenuated inversion recovery，FLAIR)像。数据来源于19个不同机构，大多采集自3.0 T磁共振成像系统。数据全部经过预处理的，包括配准到相同的解剖模板、插值到相同分辨率等。

       所有图像根据相同的标准由1~4个人进行分割，并由经验丰富的神经放射学家验证。总共分为4个区域，即：增强的核心(enhanced core)、水肿(edema)、非增强的实性核心(non-enhancing solid core)和坏死/囊性核心(necrotic/cystic core)，区域之间没有重叠。数据集中将坏死/囊性核心与非增强的实性核心合并为一个区域，简称为NET区域，以感兴趣区(region of interest，ROI)NET表示。增强的核心区域和水肿区域分别简称为ET区域和ED区域，以ROIET、ROIED表示。图1为BRATS2017数据集中一例患者的胶质瘤示意图，蓝色为T2加权像上显示的整个肿瘤区域，即NET区域、ET区域、ED区域的总和；红色为增强T1加权像上显示的ET区域，绿色为FLAIR像上显示的水肿区域，即ET区域与ED区域的总和。

点击查看大图

图1  BRATS2017胶质瘤示意图。蓝色为肿瘤区域，红色为增强的核心区域，绿色为水肿区域
Fig. 1  The segmentation of BRATS2017 dataset. Blue: the whole tumor, Red: the enhanced core, Green: edema.

1.2　特征提取

       本研究自动提取了117个特征，特征类型包括全局特征、形态特征、灰度特征和Gabor特征。以下简要介绍各类提取的特征。

1.2.1　全局特征

       不同机器、环境、患者扫描得到的磁共振图像具有不同的灰度分布。有些研究^[2,10]利用直方图匹配方法将不同图像匹配到同一灰度分布下，但是这种方法会改变肿瘤区域的灰度分布，进而丢失肿瘤的灰度信息。笔者直接将患者大脑的特征作为肿瘤分级所用的特征的一部分，这些特征称为全局特征。

       全局特征包括脑脊液、白质、灰质在T1、T1增强、T2和FLAIR像上的灰度平均值。要获取各组织的灰度值需要对大脑组织进行分割，现有的最精确的分割方法是基于深度学习的^[11]。另外，Kumar等^[12]比较了不同的图像分割算法，指出K-means算法具有较好的分割表现和较低的计算开销，与其他算法结合能够提升分割效能，考虑到特征提取步骤需要对精确性与效率进行权衡，因此，针对该文的应用场景，采用K-means算法。

       以的形式表示全局特征，其中A为脑脊液(cerebrospinal fluid，CSF)、白质(white matter，WM)和灰质(gray matter，GM)，B为μ表示平均值，C为T1、T2、T1CE、FLAIR。

1.2.2　形态特征

       形态特征包括肿瘤总体积V、肿瘤表面积A、肿瘤似圆度SPH^[13]、肿瘤表面径向分布的熵值RE、ET区域体积VET、ED区域体积VED以及NET区域体积VNET。另外，本研究还提取了VET与V的比值rET、VNET与V的比值rNET以及VED与V的比值rED。

1.2.3　灰度特征

       本研究分别计算了T1、T1增强、T2、Flair图像中ROINET、ROIET、ROIED的灰度平均值、标准差，其中最终的平均值为除以白质的灰度平均值之后得到的平均值。由于肿瘤区域的灰度一般不满足高斯分布，所以笔者还统计了上述3个区域对应图像灰度分布的峰度与偏度。峰度表征某一分布纵向偏离高斯分布的程度，偏度表征某一分布横向偏离高斯分布的程度^[14]。与上节相仿，以的形式表示灰度特征，其中A为NET、ET、ED，分别表示3类区域，B为μ、σ、K和SK，分别表示平均值、标准差、峰度和偏度，C为T1、T2、T1CE、FLAIR。

1.2.4　Gabor特征

       本研究对ROI所处的每一层都用二维盖波滤波器进行卷积，最后对ROI内体素信号进行平均，得到三维图像的Gabor特征^[15]。二维盖波滤波器在位置x，y处的公式为：

       λ为波长，θ为方向角度，γ为纵横比，它决定了卷积核的离心率^[16]，本研究取1，σ为高斯函数的带宽，它决定了像素不同领域范围的权重。上式中的滤波器为复数形式，因此卷积得到的是复图像，笔者取模图来计算Gabor特征。θ取0到之间的5个方向，取，取。

       本研究仅计算T1增强像中肿瘤区域与Flair像中水肿区域的Gabor特征，以的形式表示，其中A为TUMOR、ED，分别表示肿瘤区域和水肿，B取1~5，表示5个方向角度，C为T1CE、FLAIR，D取1~5，表示5个波长。

点击查看大图

1.3　特征选择

       特征选择方法主要包括：过滤(Filter)方法、包装(Wrapper)方法和集成(Ensemble)方法^[17]。本文结合了前两种特征选择方法进行特征选择。具体方法为：(1)使用过滤方法对特征进行排序；(2)添加排名靠前的特征作为新的特征子集，运用包装方法进行特征选择。

1.3.1　过滤方法与t检验

       过滤方法通常从特征集合中剔除与预测变量最不相关的特征，将剩余的变量用于预测。它所用的判据往往是通用的与预测变量的相关性之类的判据，与后续所使用的分类器模型无关。过滤方法得到的特征子集常常会包含冗余的特征，因为它所用的判据忽略了特征间的相关性。由于过滤方法与特定的分类器无关，所以计算开销小，泛化能力强。对于大数据集和高维度特征，过滤方法应用十分广泛。

       本研究使用双总体t检验值的大小作为特征排序规则。该类t检验能够衡量两个分布的差异，因此也能用来检验每个特征的分类能力。假设特征在不同分类下均服从高斯分布，并且它们的方差相等，本研究利用t检验选择在两个分类间分布差异最大的30个特征。

1.3.2　包装方法与遗传算法

       包装方法对特征子集的性能进行评价，它可以发现特征间的关系，选出最优子集。包装方法分类能力很强，但计算开销大，容易过拟合，限制了该方法在高维问题上的应用。本研究采用了遗传算法来搜索特征子集空间。

       遗传算法(genetic algorithm，GA)^[18]是模拟生物在自然环境中的遗传和进化过程而形成的一种自适应全局优化概率搜索算法。对于函数最优化的问题，解的空间一般非常大，在很多情况下几乎求不出解析解。遗传算法能够利用其随机搜索的特性，求出其近似最优解。

       遗传算法中，假设有n维决策向量，X可看作是由n个遗传基因所组成的一个染色体。每个元素的所有取值称为等位基因，本文算法取等位基因为0和1。根据每种特征选择与否就可以得到一个染色体。选择留一交叉验证法下支持向量机(support vector machine，SVM)分类器的准确率作为个体的适应值。

       进行优化时，首先设定初始群体，对群体中的所有个体进行评价，即计算适应值。然后将一系列遗传算子作用到群体上，判断是否达到算法终止条件，若是，将群体中具有最优适应值的个体作为问题的最优解，若否，则继续迭代。遗传算法的典型流程如图2所示。

       遗传算子主要包含选择算子、交叉算子和变异算子^[19]。选择运算将群体中适应度函数值较优的个体以某种规则遗传到下一代中。交叉运算以一定概率交换群体两个个体间的部分染色体，产生新的个体。变异运算以一定概率改变群体中某些个体的一个或者多个等位基因。

       本研究最终选取群体人口为50，最大进化的代数为100，选择算法为锦标赛算法，选择个数为2，即每次选取2个个体中最好的那个，共选取50次。交叉算法为分散交叉，交叉率为0.8，变异算法为均匀变异，变异率为0.1。

点击查看大图

图2  遗传算法步骤示意图
Fig. 2  The procedure of the genetic algorithm.

1.4　分类

       本研究使用过滤特征选择方法，对非线性的SVM^[20]、K近邻法(K nearest neighbour，KNN)^[21]和决策树(decision tree，DTREES)^[16]进行了比较，步骤如下：从最具有区分度的特征开始，逐步添加次一级区分度的特征，找到分类性能最好的特征组合。

       由于数据是不平衡的，需要对数据进行平衡。常用的方法有对数量多的分类样本下采样或对数量小的分类样本上采样，前者不能充分利用数据，后者可能会增加噪声点。根据每种分类算法的特性来分别解决这个问题。支持向量机中一个重要的参数是惩罚因子C。较大的C值会给误分类一个较大的惩罚。本研究还采用带权重的支持向量机^[22]来加大对小样本分类错误的惩罚。在K近邻法中，使用带距离权重的分类判决方法，这使得距离较近的点对分类具有较大贡献。对于决策树，则使用了对数据不平衡不敏感的CART算法^[23]。

1.4.1　SVM

       SVM分类算法的主要思想是找到一个独一无二的超平面，这个超平面能最大化两个分类之间的间隔。给出l组训练数据，x为输入样本，y为标签。

       软间隔的SVM方法需要得到下列优化问题的解：

       约束条件为：

       训练集的特征矢量被映射到高维空间中，C是惩罚因子，用来平衡间隔最大化的要求与分类错误之间的矛盾，取C=0.5。若对不同标签的样本进行不同惩罚，则称为带权重的，这里取正类负类的权重比为1∶1.5。该优化问题的一种简便解法是解它的对偶问题，即：

       约束条件为：

       最后给出一个数据x，得到该数据的分类为，令，则函数K称为核函数，这里取高斯核函数为，，n为特征数量。

点击查看大图

点击查看大图

点击查看大图

点击查看大图

点击查看大图

1.4.2　K近邻法

       KNN算法常用来解决分类问题，它没有显式的学习过程。典型的KNN算法步骤如下：给定一个训练数据集，对新的输入实例，在训练数据集中找到与该实例距离最相近的K个实例，根据分类决策规则将该实例分为某个类。该算法有3个基本要素^[24]：距离度量、k值选择和分类决策规则。常用的分类决策规则有多数表决，即K个实例中较多的那一类为最终分类。这里采用欧氏距离度量，这个距离还能带上权重，取权重为。较小的K值会使模型变复杂容易过拟合，而较大的K值会增大分类误差，这里取K为3。

1.4.3　决策树

       决策树模型是定义在特征空间与类空间上的条件概率分布，主要优点有：模型具有可读性，分类速度快^[24]。决策树算法学习一般包含3个步骤：特征选择、决策树生成和决策树修剪。本文决策树是基于分类回归树(classification and regression tree，CART)算法的，该算法通过递归地构建二叉树来生成决策树。它使用基尼(Gini)指数选择最优特征，同时决定该特征的最优二值切分点。对于给定样本集合D，基尼指数定义为：

       K为标签个数，当K=2时即为二分类问题，标签为k的样本个数。根据特征A的取值分割成D1、D2两部分，则在该条件下基尼指数定义为：

       决策树修剪是通过剪枝算法完成的，算法从完全生长的决策树底端减去一些子树，分为两步^[24]：首先从生成算法的决策树底端开始不断剪枝，直到根节点，形成子数序列，然后通过交叉验证法在独立的验证数据集上测试，选择最优子树。

点击查看大图

点击查看大图

2　结果

2.1　基于Filter方法的分类

       本研究用留一交叉验证法对BRATS2017数据集的包含分级信息的285例数据进行分类。首先，使用基于过滤方法的特征选择方法，特征排序方法为t检验法。选出前30个特征依次添加，并用支持向量机、K近邻法和决策树进行分类。本研究使用了LIBSVM工具箱^[25]来进行SVM算法分类。图1为支持向量机、K近邻法、决策树的分类准确率与特征数量的关系图。

       笔者确定3种方法达到最高准确率时的特征子集，计算对应的敏感性和特异性，见表1。由表1可知，支持向量机算法具有最高的准确率与敏感性，所以选择该算法作为下一步实验的分类器。

点击查看表格

表1  不同方法分类的准确率、敏感性、特异性
Tab.1  Classification accuracy, sensitivity, and specificity of various methods

2.2　基于遗传算法特征选择的分类

       本研究以过滤方法排序得到的前30个特征作为初值，进一步使用基于遗传算法的包装方法进行特征选择。选择上节中选择的SVM作为分类器。图4为每一代适应度函数的平均值、最大值与进化代数的关系图。图5为支持向量机分类的受试者工作特征曲线(receiver operating characteristic，ROC)图，黑线为基于过滤方法的，红线为基于遗传算法的。

       由图4可以看出，算法迭代结束时最小函数值为-93.33，即达到的最大准确率为93.33%，证明遗传算法搜索得到的特征组合效果最好，优于过滤方法的结果。

       由图5可以看出，GA-SVM方法ROC曲线下的面积(area under the curve，AUC)为0.94，大于Filter-SVM方法的0.90，证明具有较好的分类表现，该方法在同等敏感性下有着更低的假阳率，即达到相同效果所付出的代价更低。

点击查看大图

图3  分类准确率与添加特征数量关系图。决策树分类算法随着特征数增加准确率趋于平稳。K近邻分类算法在特征数很少时准确率就达到最大，随后趋于平稳。支持向量机分类算法能达到的准确率最高，但是随着特征数增加，准确率急剧下降
图4  适应度函数平均值及最小值随遗传代数变化的情况
图5  ROC曲线图。Filter特征选择方法的ROC曲线下面积(area under the curve，AUC)=0.90。基于Wrapper与遗传算法的特征选择方法的AUC=0.94
Fig. 3  Classification accuracy vs number of features used. Accuracy of decision tree is quite stable with the increase of the number of features. Accuracy of KNN reaches the maximum when the number of features is small. Accuracy of SVM is the highest among three classifiers, but it rapidly falls off when the number of features increases.
Fig.4  The best or mean fitness function value vs generations.
Fig.5  ROC curve. Filter feature selection method, AUC=0.90. Feature selection method based on Wrapper and genetic algorithm, AUC=0.94.

2.3　特征选择结果

       表2为上两节中4种方法选取得到的特征。由表2可知，基于过滤特征选择方法的支持向量机与决策树算法挑选出了同样的特征，而K近邻法选出了较多特征，但是特征之间有冗余。GA-SVM方法则在这些特征的基础上又选择出了Gabor特征和全局特征。

点击查看表格

表2  不同方法特征的选择结果
Tab.2  The feature selection results of various methods

3　讨论

3.1　脑胶质瘤分级的自动化

       本文利用影像组学方法对脑胶质瘤进行了自动分级。脑胶质瘤的自动分级一般有如下步骤：(1)图像分割；(2)特征提取；(3)特征选择；(4)分类。最常用的特征提取方法有离散小波变换(discrete wavelet transform，DWT)和灰度共生矩阵(gray level co-occurrence matrix，GLCM)；在特征降维方面，主成分分析法(principal component analysis，PCA)和遗传算法被广泛使用^[26]；随机森林(random forests，RF)由于其算法特性，更适合处理高维问题和多分类问题，因此该算法对于脑胶质瘤分级最精确^[8]。本文则提取了多种类型的特征包括全局、形态、灰度、纹理各个方面；并结合两种特征选择方法，进行了多步骤的特征降维，由于是单分类问题，本研究只比较了3种常见的分类算法。笔者首先采用基于t检验的过滤方法进行特征选择，发现支持向量机分类算法能达到的准确率和敏感性最高，但是随着特征数增加，准确率急剧下降。支持向量机算法具有较好性能的原因是：该算法的高斯核函数将特征映射到了高维空间，提高了特征的分类能力。结果表明，支持向量机算法适合作为脑胶质瘤分级问题中的分类器。

       然后，采用基于遗传算法的包装方法进行特征选择，算法得到的准确率和AUC值都超过了过滤方法，结果表明，遗传算法挑选出了特定的特征组合，提升了脑胶质瘤的分级诊断效能。包装方法能够搜索更大的特征空间，但是需要更长的搜索时间，而遗传算法是一种自适应的全局优化概率搜索算法，弥补了该方法的缺点，加快了搜索过程。

3.2　影像组学在常规磁共振图像上的应用

       基于影像组学的方法对脑胶质瘤分级进行的研究相对较少，值得注意的是，邓慧媛等^[27]利用Logistic多元回归相关性分析筛选出了鉴别高低级别脑胶质瘤效能较高的DWI序列特征，高于T1加权增强像、T2加权像的鉴别效能。另外，林坤等^[28]的研究得出多模态MRI检查技术的分级诊断效能高于常规增强MRI及DWI检查。本文则结合公开数据集的大数据量和影像组学发掘了常规磁共振图像在脑胶质瘤分级问题上的巨大潜力。

       在特征选择方面，基于过滤特征选择方法的支持向量机与决策树挑选出了同样的特征，但是它们的分类性能不同；K近邻法选出了较多特征，但是特征之间有冗余。这3种方法选出的特征基本与先验知识相符，即肿瘤越大、肿瘤区域灰度值越高，肿瘤恶性可能性越大。基于遗传算法的特征选择方法则在这些特征的基础上又选择出了Gabor特征和全局特征。这些特征增强了分类性能，而且人的肉眼无法观察到，这凸显了影像组学方法与人相比可能具有的优势。本研究从T1、T1增强和FLAIR像中都提取到了有用的特征，这也表现了常规磁共振图像在脑胶质瘤分级中不可忽视的作用。

3.3　本文工作的不足与局限性

       目前，影像学检查经常使用于脑胶质瘤的分级诊断，常规磁共振图像对于病变的解剖定位意义重大，但是不能同时显示肿瘤的代谢、血液动力学、分子及细胞结构的信息^[29]。这也提示需要克服常规磁共振图像的局限性，将多种MR技术相结合，彼此取长补短，可更好地对脑胶质瘤进行分级诊断。

       另外，本研究所使用的公开数据集中已经提供脑胶质瘤的分割结果。在脑胶质瘤分级的实际应用中，还需要在ROI区域对肿瘤进行分割，在分割结果上进行特征提取。因此，在之后的工作中，对脑胶质瘤的分割方法进行研究也是重点之一。

       总之，本研究对基于常规磁共振图像的脑胶质瘤分级进行了研究，实验结果表明，基于常规磁共振图像，利用影像组学方法，选择合适的特征组合，可以对脑胶质瘤进行自动分级。通过对分类方法的比较，发现非线性支持向量机对脑部肿瘤分级效果最好，在特征筛选方面，笔者提出先利用过滤方法筛选出部分特征，然后在过滤的基础上进一步使用包装方法，用遗传算法选择了最后的特征子集。