冯 宝, 李昌林, 李 智, 刘壮盛

(1. 中山大学 生物医学工程学院, 广东 广州 511400; 2. 桂林航天工业学院 电子信息与自动化学院, 广西 桂林 541004; 3. 中山大学附属江门市中心医院 放射科, 广东 江门 529000)

乳腺癌是发病率最高的女性恶性肿瘤[1-3].肿瘤对周围脉管(包括淋巴管和血管)的侵犯状态(lymphovascular invasion,LVI)是预测乳腺癌患者局部复发和转移的重要指标[4-6].目前,对肿瘤切除术后的标本行病理检查是诊断乳腺癌LVI状态的金标准和常规方法,但该方法的LVI信息是在肿瘤切除术后获得,这对乳腺癌病人术前新辅助化疗和手术方式优化无参考价值.医学影像技术能术前无创提供整个肿瘤信息,在指导临床治疗上潜力巨大[7-8].在乳腺癌患者的影像诊断中,动态对比增强磁共振成像(dynamic contrast enhanced magnetic resonance imaging,DCE-MRI)凭借高敏感性、高空间分辨率等特点成为乳腺癌影像检查的主要工具.因此,开展基于DCE-MRI的乳腺癌患者术前LVI状态预测,对指导和优化患者诊疗方案意义重大.

影像组学是当前医工交叉学科的研究热点之一[9].它利用临床、基因和影像等信息,通过模式识别方法来提取病灶区的有效特征,构建描述肿瘤表型的定量辅助诊断模型,再进行计算机辅助诊断.在影像组学的研究中,肿瘤病灶的准确分割和定量辅助诊断模型的构建是两个关键技术.其中病灶的准确分割是一个重要的预处理步骤,病灶区分割的准确性将直接影响后续影像特征提取的准确性和稳定性,进而决定模型构建的有效性；而模型构建的有效性可以作为肿瘤病灶分割的评价标准之一.

近年来,基于影像组学的乳腺癌影像分析受到广泛关注,已取得较多成果.刘桐桐等[10]利用动态活动轮廓模型对乳腺癌病灶分割,然后提取形态、纹理等特征,利用支持向量机建立乳腺癌雌激素受体表达情况的预测模型.该算法获得的AUC值为0.738,准确率为0.712,敏感度为0.732,特异度为0.758.吴佩琪等[11]提出自动阈值方法获取病灶初始轮廓,同时结合医生手动分割调整来获取乳腺癌DCE-MRI图像的病灶,提取了灰度共生矩阵(GLCM)等纹理特征,利用逻辑回归方法建立乳腺导管癌病理分级模型.该方法获得的AUC值为0.691,准确率为0.676,敏感度为0.905,特异度为0.345.Huang等[12]提出基于模板匹配算法对乳腺癌DCE-MRI病灶图像分割,提取了基于GLCM的纹理特征,采用模糊C均值聚类方法对乳腺肿瘤良恶性分类.Fan等[13]利用医生手动分割方法来分割乳腺癌病灶,通过逻辑回归方法建立亚型分类模型(Luminal A型和 non-Luminal A型),获得AUC值为0.718.此外,深度学习方法也受到关注.孙利雷等[14]提出多尺度卷积神经网络对乳腺肿块良恶性分类,得到AUC值为0.713.

上述方法一定程度上提高了乳腺癌计算机辅助诊断效率,但病灶区提取还是以半自动分割、医生手动分割为主,耗时耗力,重复性欠佳,不适合大数据集上的影像分析.此外,大部分方法采用单分类器,应用在具有多样性的乳腺癌影像分析中存在一定局限性,临床应用很难达到满意的效果.

针对乳腺癌患者术前LVI状态预测问题,本文提出活动轮廓模型和影像组学结合的分析方法.首先对乳腺癌病灶区在DCE-MRI图像中存在的边界模糊、内部亮度不均匀等特点,提出一种基于后验概率和模糊速度函数的活动轮廓模型方法来分割病灶.通过在小波域下构建基于后验概率的活动轮廓模型的区域项,同时利用模糊速度函数改进活动轮廓模型的边界项,能提高分割的准确性.在构建定量辅助诊疗模型时,提取形态、灰度、纹理等影像特征来表征乳腺癌病灶,再用最小收缩和选择算子(least absolute shrinkage and selection operator,LASSO)方法选出有效特征.最终利用集成分类随机森林方法构建分类器.通过集成多棵决策树提高分类准确率,完成乳腺癌患者术前LVI状态预测任务.

1 LVI状态预测方法

乳腺癌患者术前LVI状态预测问题有两个关键技术点:病灶的准确分割和分类预测模型的有效构建.病灶分割是重要的预处理步骤,分割准确度直接影响特征提取的准确性和稳定性,从而影响分类模型构建,因此病灶的准确分割是本文的关键;而分类模型的有效构建可以作为病灶分割的评价指标之一,更好地对乳腺癌患者进行术前LVI状态预测.

1.1 基于后验概率和模糊速度函数的乳腺癌灶分割

由于部分乳腺癌病灶沿腺体分布,导致其形态极不规则,边界模糊,对比度低,准确分割困难极大[15].针对上述问题,提出一种基于后验概率和模糊速度相结合的活动轮廓模型方法来分割病灶:①利用小波变换将图像灰度信息转变为小波能量,增强目标区域与背景差异,再在能量图中计算每个像素点的后验概率,建立基于后验概率的活动轮廓模型区域项;②结合Gabor纹理特征和灰度特征构建模糊速度函数,并引入到活动轮廓模型中作为边界检测项;③计算构建的活动轮廓模型的梯度下降流,驱动轮廓曲线演变,完成病灶的分割.

1.1.1 活动轮廓模型的能量函数定义

假定C是一条轮廓曲线,该曲线处于封闭状态并把曲线区域划分为内部和外部;因此,定义活动轮廓模型的能量函数为

E(C,φ)=μERegion(φ)+βEContour(C) .

(1)

式中:ERegion(φ)和EContour(C)分别是控制轮廓曲线演变的区域项和边界检测项;μ和β均是常数(文中μ和β均取1)[16];φ是水平集函数,在C内部取正数,在外部取负数.

1.1.2 小波能量引导下基于后验概率的活动轮廓模型区域项

图像任意像素点(x,y)的小波能量定义为

W(x,y)=EA(x,y)+EH(x,y)+

EV(x,y)+ED(x,y) .

(2)

式中,EA(x,y),EH(x,y),EV(x,y)和ED(x,y)分别是图像的近似细节、水平细节、垂直细节和对角细节小波能量.

将原始DCE-MRI图像分成病灶区和背景区,假定这两个区域中像素的小波能量值都满足均值为ur、方差为σr(r=1,2)的高斯分布,采用一个混合高斯模型[17]来表达整张图像.将图像像素定义为{x0,x1,…,xN}(N是像素个数),xi是第i个像素点小波能量值,当区域r有xi值时,混合高斯分布定义为

P(xi,r)=P(xi|r)Pr.

(3)

(4)

通过最大期望算法[18]求均值ur、方差σr和Pr,得到区域r中有xi值的后验概率P(r|Wxy,θ(n)),θ(n)为通过最大期望算法求解最大似然函数获得的模型参数,Wxy表示像素点(x,y)处的小波能量.故活动轮廓模型区域项为

(5)

式中：H(·)是Heaviside函数；Ω为图像区域.

图1是乳腺癌病灶小波能量分布和后验概率分布图.由图1b可看出病灶和正常腺体灰度值相近,难以区分.由图1c可看到正常腺体的小波能量没有病灶区强,可据此区分病灶和正常腺体,此外还可看到病灶区内部亮度不均匀.由图1d可看出病灶区得到明显增强.

1.1.3 构建基于模糊速度的边界检测项

利用灰度特征和Gabor纹理特征相结合的模糊聚类算法,计算模糊隶属度,构建模糊速度函数,并将其引入活动轮廓模型中作为边界停止项.

模糊聚类算法需要多次迭代寻找最优的目标函数,目标函数[19]被定义为

(6)

式中:xl是二维特征向量;U=(ulh)是病灶区和背景区的模糊隶属度矩阵,ulh是xl属于h类(h=1表示背景,h=2表示病灶)的隶属度大小;v=[v1,v2]是聚类中心向量;m是模糊控制常数;vh是距离中心,d2(xl,vh)是xl和vh的距离.重复迭代计算ulh和vh,所有的二维向量被正确分类时,J(U,v)目标函数达到最小,隶属病灶的模糊隶属度矩阵Z由ulh得到.

为了更好分割病灶,需设定较好的模糊速度函数:离病灶边界处越远,模糊速度越大;边界处,模糊速度为零.

将图像像素值(x,y)的模糊速度函数定义为

V(x,y)=eα[Z(x,y)-0.5]2-1 .

(7)

式中：α是常数；Z(x,y)是(x,y)隶属病灶区的模糊隶属度矩阵.

把V(x,y)加到活动轮廓模型中,边界检测项E2(C)定义为

(8)

1.1.4 活动轮廓模型能量泛函演化

将式(5)区域项和式(8)边界检测项代入式(1)能量函数,可得改进后的活动轮廓模型能量泛函:

(9)

病灶分割是最小化能量函数E(φ).通过构建E(φ)关于水平集函数φ的Euler-Lagrange方程,则

(10)

根据变分原理,求出φ的梯度下降流为

(11)

1.2 特征提取和特征选择

分割完病灶后,提取其形态特征共4个[20](直径、体积、偏心率和密度),灰度特征共3个[20](方差、偏度和峰度).为提高特征的鲁棒性和可重复性,对三维重建的病灶进行各向同性重采样[21]和灰度离散化[22]处理,获得不同参数下的特征,提取了基于灰度共生矩阵(GLCM)在不同重采样系数(重采样系数=[0.8,1,1.2,1.5,2])和不同灰度阶(灰度阶=[8,16,32,64])下的纹理特征(纹理能量、对比度、相关性、同质性、纹理方差、和均值、纹理熵、差异性、自相关)[23-24]，共5×4×9=180个.

由于特征间的相关性和冗余性会降低分类准确率,过多无关特征易增加分类器复杂性,出现过拟合;故需对提取的特征进行优化选择[25].本文选择“最小收缩和选择算子(LASSO)”方法进行特征选择[26].这是一种能稀疏特征的方法,本质是在回归模型中引入L1范数约束项,将一些不显著特征对应的回归系数置0,实现特征的自动选择.

1.3 基于集成分类器随机森林的LVI状态预测模型

特征选择后,采用集成分类器随机森林算法构建预测模型:通过集成学习的思想将多棵决策树进行集成[27].集成分类器随机森林算法能处理高维特征的输入样本,训练速度快,可有效避免过拟合.

集成分类器随机森林模型构建步骤如下:

①选择乳腺癌LVI的训练集数据,采用Bootstrap方法随机,筛选出k个样本,生成m个训练样本子集.

②每个样本子集训练一棵分类决策树,分别训练m个决策树模型.每次训练过程中,随机选取n个特征作为决策树每个节点的候选分裂特征,然后通过计算基尼系数[24]来选取最好的特征.

③利用最优特征对决策树每个节点进行分裂,划分决策树的左、右子数,直到节点的所有训练样本属于同一类别.决策树尽可能在不修剪的前提下生长,直到达到最大迭代次数或者满足停止标准,则停止该枝生长.多个训练子集形成多棵决策树.

④将所有训练得到的决策树组成随机森林,形成随机森林分类器.

⑤用训练好的随机森林分类器对测试集数据进行预测.再根据多数投票准则,将票数最多的分类结果作为最终预测结果.

1.4 分割和分类算法性能评价指标

1.4.1 分割性能评价指标

乳腺癌DCE-MRI的病灶分割以影像科医生的手动分割结果为金标准.医生首先判定病灶区,然后进行手动精准分割,最后将分割的病灶区作为金标准.分割指标采用阳性率(true positive ratio,TPR)、假阳性率(false positive ratio,FPR)和相似度(similarity degree,SD)来评判:

(12)

式中:Am是金标准区域;Aa是算法分割区域.TPR是Aa对Am的覆盖率,TPR越大,分割效果越好.FPR是Aa的背景区域对Am的比值,FPR值越小,Aa包含的背景区域越小,分割效果越好.SD值是Aa和Am的相似比值,SD越大表明分割效果越好.

1.4.2 分类性能评价指标

为了评估分类算法性能,采用3个指标:准确率(Acc)、敏感度(Sen)、特异度(Spe),计算公式为

(13)

式中:TP和FP分别表示将LVI状态正确和错误地划分为阳性的病例数;TN和FN分别表示将LVI状态正确和错误地划分为阴性的病例数.

2 实验结果与分析

2.1 实验数据描述

本文选择中山大学附属江门市中心医院的149例乳腺癌患者DCE-MRI数据,90例为训练集,59例为测试集.根据术后病理结果,将患者LVI状态分为阳性(肿瘤侵犯周围脉管)和阴性(肿瘤没有侵犯周围脉管).训练集包括22例LVI阳性和68例LVI阴性,测试集包括22例LVI阳性和37例LVI阴性.

DCE-MRI图像数据的采集设备为西门子 Magnetom Espree 1.5T MR扫描仪,配备4通道乳腺专用相控阵表面线圈.全部患者为俯卧位.选择轴向三维T1梯度回波法,利用脂肪饱和序列,最终采集DCE-MRI图像(TR:4.6 ms;TE:1.7 ms;翻转角:7°;像素大小:1.0 mm×1.0 mm;FOV:280 mm×340 mm;矩阵大小:280×340;截面厚度:1.0 mm;差距:0 mm;采集时间:75 s/相).所有患者均顺利完成MRI检查.对于多发乳腺病变的病人,仅选取最大的癌变病灶.

2.2 不同分割算法的分割性能比较

为了检验本文分割算法性能,选择4个形态不规则、边界模糊和对比度低的数据进行验证.将本文分割算法和基于边界的ACM方法[28]、基于区域的ACM方法[29]对比分析,图2是对比结果.

从图2b知,病灶区和周围正常腺体对比度较低,采用基于边界的ACM方法分割,发现正常腺体较多,部分边界易泄漏.图2c采用基于区域的ACM方法,分割效果也不佳,出现误分割.比较图2a和图2d知,本文算法接近医生分割结果.

基于区域的ACM方法和基于边界的ACM方法无法精准分割病灶的原因是:病灶和正常腺体灰度值相近,在轮廓曲线演变时只用了图像梯度信息;而图像对比度低时,无法检测低对比度边界是否为病灶边界,导致边界泄漏,无法精准分割.

为了验证算法分割性能,选择20个临床样本.采用TPR,FPR,SD的平均值作为分割综合评价指标.从表1看出,基于区域的ACM方法和基于边界的ACM方法能获得较好的TPR值,但FPR值太大,SD值较小,表明这两种分割方法尽管能把病灶分割出来,但出现大量背景,相似度低;而本文算法分割背景少,相似度高,分割性能更优.

表1 20个临床样本分割结果
Table 1 20 clinical sample segmentation results

算法TPRFPRSD基于边界的ACM方法0.96391.97670.4041基于区域的ACM方法0.94211.17420.4994本文分割算法0.91790.28070.7280

2.3 基于LASSO的特征选择

当全部特征经LASSO方法选择后,最终选出两个任务相关有效特征,分别是方差(灰度特征)和熵(纹理特征).图3为LASSO特征选择可视化图.图3a是10重交叉验证结果,虚线①定义调整参数λ最优值,此时模型提供对数据的最佳拟合;图3b是全部特征LASSO回归系数分布图,虚线②表明当λ达到最优时非零回归系数对应的两个特征.

2.4 基于集成分类器随机森林的LVI状态预测

在训练集成分类器随机森林模型时,决策树的数量N1对随机森林预测的准确性起到重要作用.表2列举了在选择不同决策树数量时得到的分类性能.从表中看出,N1大于100或小于100时,敏感度、特异性和准确率均不如N1为100时的好,只有当N1=100时分类性能最好.

表2N1取不同值对应的分类性能比较
Table 2Comparison of classification performancecorresponding to different values ofN1

N1TPFPFNTNSenSpeAcc104183340.57140.65380.6441203193340.50000.64150.6271302202350.50000.63640.6271402202350.50000.63640.6271502203340.40000.62960.6102602203340.40000.62960.6102702202350.50000.63640.6271802202350.50000.63640.6271902202350.50000.63640.62711004182350.66670.66040.66102003192350.60000.64810.64413003192350.60000.64810.64414002202350.50000.63640.62715002202350.50000.63640.62716002202350.50000.63640.6271

2.5 不同分割算法在随机森林模型的分类比较

采用三种分割算法进行乳腺癌病灶分割,然后构建不同分割算法下的集成分类器随机森林模型,再进行乳腺癌患者LVI状态预测.为了公平起见,三种分割算法均采用LASSO选出的2个特征.从表3可以看出,本文提出的分割算法敏感度、特异度和准确率均比基于边界的ACM方法和基于区域的ACM方法好,表明本文提出的分割算法能够提高分类器的分类性能.

表3 不同分割算法下随机森林预测模型的分类性能Table 3 Classification performance of random forest prediction models under different segmentation algorithms

2.6 不同分类器对测试集的分类性能比较

采用6种不同的分类器对测试集进行分类性能对比.从表4看出,集成分类器随机森林比其他分类器性能好,尽管贝叶斯线性回归的特异度比集成分类器随机森林高,但敏感度和准确率均较低.ResNet18的敏感度高,但特异度和准确率低.

表4 不同分类器对测试集的分类性能Table 4 Classification performance of test sets by different classifiers

3 结 语

本文结合活动轮廓模型和影像组学方法,完成乳腺癌患者术前LVI状态预测.针对乳腺癌病灶在DCE-MRI图像中边界模糊、亮度不均匀的特点,提出基于后验概率和模糊速度函数的活动轮廓模型方法,完成病灶的分割,然后进行影像学特征提取,并利用LASSO方法选择有效特征,最后利用集成分类器随机森林模型进行测试.实验结果表明，本文提出的分割算法可以提高乳腺癌病灶分割的准确性,所构建的集成分类器随机森林模型能够对乳腺癌患者的术前LVI状态进行较准确的预测.