鼻咽癌（nasopharyngeal carcinoma，NPC）是人体鼻咽部最常见的恶性肿瘤. 据世界卫生组织报告显示，全世界约80%鼻咽癌患者集中出现于中国，其余患者大部分常见于东南亚、中东及北非地区^[1-2]. 据统计，广州市的NPC发病率为17.8人/10万人，且呈现发病率上升、发病年轻化的趋势^[3]. 大部分患者在确诊时已错过最佳治疗时期，因此NPC的早期诊断对抓住最佳治疗时机至关重要. 目前，NPC的治疗方式主要是放疗，而NPC病灶的准确定位，是制定放疗计划与实施放疗的重要依据.

MR图像由于具有分辨率高、安全性强等优点，常被应用于NPC的诊断和病灶定位. 在制定治疗计划前，临床医生通常需要检查每张MR图像中是否存在病变，并手动勾画出肿瘤边界，由于不同NPC患者间肿瘤形状、大小、位置复杂多变，因此该过程存在工作负担重、对医生经验要求高、受主观影响大等缺点. 因此，研究者纷纷寻求NPC病灶区域的自动勾画方法.

NPC肿瘤分割方法主要包括阈值法、区域增长法、统计理论等传统图像分割方法^[4-5]，以及支持向量机、人工神经网络等机器学习方法^[6-7]. 此类方法通常需要特征提取、降维等人工干预过程，存在模型鲁棒性差、对噪声敏感等缺点，因此，使用此类方法难以实现NPC的快速自动分割.

近年来，卷积神经网络（convolutional neural network，CNN）等深度学习（deep learning，DL）方法由于能够直接自动提取图像中从低阶到高阶、从抽象到具体的特征，被广泛应用于医学图像分类、分割以及配准中. 例如，在NPC分割方面，Men等^[8]采用编码-解码的全卷积神经网络对NPC患者的CT图像进行分割，并与VGG网络^[9]进行比较。结果表明：该网络能够有效提高NPC的分割效果. 同样，Li等^[10]采用编码−解码的全卷积神经网络对27例NPC患者的MR图像进行分割，采用留一交叉验证法实现了NPC的自动有效分割，分割Dice值达到0.89. 另外，Ma等^[11]采用卷积神经网络对30例NPC患者的T1W模态MR图像病灶区域进行分割，再利用3D图割算法对分割结果进行改进，得到Dice值为0.851. 以上研究所采用的模型均为2D模型，并未考虑3D空间上图像特征的关联，网络结构较简单，且实验数据较少、仅采用了单一模态的图像进行分割，而单种模态影像数据所反映的NPC肿瘤局部信息有限，故模型鲁棒性有待提高.

虽然2D网络具有参数少、模型拟合速度快的优点，但因其没有有效利用层间拓扑结构信息，易出现分割结果准确性较差、边界不连续的情况. 3D网络能克服此不足，然而，3D网络具有参数量较大、模型拟合慢甚至难拟合等缺点. 为此，本研究在Li等^[12]提出的H-DenseUnet模型的基础上，将2D与3D信息融合，即以2D网络的快速分割结果来指导3D模型的学习并实施分割，采用T1W、T2W、T1C三种模态的MR图像，提出一种新的多模态多维信息融合的深度分割方法，建立多模态2D-ResUNet与3D-ResUNet多维特征融合模型，以实现NPC病灶的自动准确分割.

深度卷积神经网络结构如图1所示，主要包括多模态2D-ResUNet结构、3D-ResUNet结构以及2D+3D融合层. 整个模型的输入为3D图像 $I \in {{\bf{R}}^{1 \times 384 \times 384 \times {b} \times 3}}$，尺寸为1×384×384×b×3，分别表示输入网络的批量大小（batch size）、图像高度（h）、图像宽度（w）、图像深度（b）以及图像通道数（c），其中，图像通道数c=3表示3种图像模态T1W、T2W以及T1C.

首先，假设函数 ${T}$表示将3D图像转化为2D的过程， ${{T}^{ - 1}}$表示该变换的逆过程，则 ${{I}_{{\rm{2d}}}} = {T}\left( {I} \right)$，得到3种模态2D图像 ${I_{{\rm{2d - T1W}}}}$、 ${I_{{\rm{2d - T2W}}}}$以及 ${I_{{\rm{2d - T1C}}}}$；其次，将2D网络表示为 ${{f}_{{\rm{2d}}}}$，3D网络表示为 ${{f}_{{\rm{3d}}}}$，则多模态2D图像经多模态2D-ResUNet后的特征图和概率图分别为

(1) $ \begin{array}{*{20}{c}} {{{F}_{\rm{2}}}_{\rm{d}} = {{f}_{\rm{2}}}_{\rm{d}}({{I}_{{\rm{2d - T1W}}}},{{I}_{{\rm{2d - T2W}}}},{{I}_{{\rm{2d - T1C}}}};{\theta _{\rm{2}}}_{\rm{d}}),}\\ {{{F}_{\rm{2}}}_{\rm{d}} \in {{\rm{R}}^{b \times {\rm{384}} \times {\rm{384}} \times {\rm{16}}}};}\\ {{{y}_{\rm{2}}}_{\rm{d}} = {{f}_{\rm{2}}}_{{\rm{dcls}}}({{F}_{\rm{2}}}_{\rm{d}};{\rm{ }}{\theta _{\rm{2}}}_{{\rm{dcls}}}) ,\;\;{{y}_{\rm{2}}}_{\rm{d}} \in {{\rm{R}}^{b \times 384 \times 384 \times {\rm{2}}}}}. \end{array} $

式中： ${\theta _{{\rm{2d}}}}$、 ${\theta _{{\rm{2dcls}}}}$分别为卷积网络和分类网络的参数. 为了将2D网络所得结果与3D网络结合， ${{F}_{{\rm{2d}}}}$与 ${y_{{\rm{2d}}}}$需要经过以下反变换，得到对应的3D特征图：

(2) $ \begin{array}{*{20}{c}} {{{\widehat {F}}_{{\rm{2}}{\rm{d}}}}{\rm{ = }}{{T}^{{\rm{ - 1}}}}{\rm{(}}{{F}_{{\rm{2}}{\rm{d}}}}{\rm{)}},\;\;{{\widehat {F}}_{{\rm{2}}{\rm{d}}}} \in {{\bf{R}}^{{\rm{1}} \times {\rm{384}} \times {\rm{384}} \times b \times {\rm{16}}}};}\\ {{{\widehat {y}}_{{\rm{2}}{\rm{d}}}}= {{T}^{{\rm{ - 1}}}}{\rm{(}}{{\rm{y}}_{{\rm{2}}{\rm{d}}}}{\rm{)}},\;\;{{\widehat {y}}_{{\rm{2}}{\rm{d}}}} \in {{\bf{R}} ^{{\rm{1}} \times {\rm{384}} \times {\rm{384}} \times b \times {\rm{2}}}}}. \end{array} $

将 ${\widehat {y}_{{\rm{2d}}}}$与 $I$合并，一同输入3D-ResUNet中，可得到3D网络的特征图：

(3) $ {{F}_{{\rm{3}}{\rm{d}}}} = {{f}_{{\rm{3}}{\rm{d}}}}({I},{\widehat {y}_{{\rm{2}}{\rm{d}}}};{\rm{ }}{\theta _{{\rm{3}}{\rm{d}}}}),\;\;{{F}_{{\rm{3}}{\rm{d}}}} \in {{\bf{R}}^{{\rm{1}} \times {\rm{384}} \times {\rm{384}} \times {b} \times {\rm{16}}}}. $

式中： ${\theta _{{\rm{3d}}}}$为3D网络的参数. 将 ${\widehat {F}_{{\rm{2d}}}}$与 ${{F}_{{\rm{3d}}}}$求和后得到Z，再输入2D+3D融合层 ${{f}_{{\rm{HF}}}}$，进行卷积计算得到H，再经过分类层 ${{f}_{{\rm{HFcls}}}}$，可得到3D的分割结果 ${{y}_{\rm{H}}}$.

(4) $ \begin{array}{*{20}{c}} {{Z = }{\widehat {F}_{{\rm{2d}}}}{\rm{ + }}{{F}_{{\rm{3d}}}},}\\ {{H = }{{f}_{{\rm{HF}}}}({Z};{\theta _{{\rm{HF}}}}),}\\ {{{y}_{\rm{H}}}{\rm{ = }}{f}_{{\rm{HFcls}}}}({{H};{\theta _{{\rm{HCFcls}}}})}. \end{array} $

式中： ${\theta _{{\rm{HF}}}}$与 ${\theta _{{\rm{HFcls}}}}$分别为融合层的卷积层 ${{f}_{\rm{HF}}}$的参数与分类层 ${{f}_{\rm{HFcls}}}$的参数.

2D-ResUNet与3D-ResUNet的网络层次结构及相应参数如表1所示，采用带残差结构的ResUNet构建模型是因为残差结构能够有效解决梯度消失等问题，并且在计算机视觉任务中已被证明性能优于简单进行卷积层堆叠的网络.该网络结合了2D网络拟合速度快和3D网络空间信息利用充分的优点，利用2D网络的分割结果去指导3D模型的拟合，从而有效地实现模型的训练和测试.

多模态2D-ResUNet结构如图2所示，主要包括3种模态图像的编码器、多模态卷积结构、解码器以及编码器-解码器间的跳跃连接结构. 编码器由带残差结构的卷积块与最大池化层组成，解码器由带残差结构的卷积块和反卷积层组成，将3种模态图像经编码器中的每一池化层进行池化，再将经池化后的3种特征图进行卷积运算，从而实现多模态2D特征的融合. 具体实现过程如下：3种模态图像经过相同层次的卷积和池化后的特征图具有相同的尺寸，通过在图像深度的方向上对3种特征图进行合并，可生成深度为3的3D特征图，再通过大小为（1，1，3）的卷积核以（1，1，3）为步长对该特征图进行卷积，从而将特征图深度转化为1，实现3种模态特征的深度融合. 经多模态卷积后的特征图，通过跳跃连接结构与解码器相同层次的特征图相加，从而减少反卷积过程中的信息丢失.

如图1所示，多模态3D-ResUNet的输入为多模态2D-ResUNet的分割概率图 ${{y}_{2{\rm{d}}}}$经变换 ${{T}^{ - 1}}$后所得概率图 ${{\widehat {y}}_{\rm{2d}}}$与模型的原始输入 $I$合并后的图像，因此输入尺寸为1×384×384×b×5，其中“5”表示网络输入为5通道，分别为3种模态图像以及多模态2D-ResUNet网络所得肿瘤与背景2种区域分割概率图. 将该图像输入到带有跳跃连接以及残差结构的3D-ResUNet网络中，即可实现3D分割. 因此，该网络将3种模态图像作为输入图像的三通道，将多模态2D-ResUNet分割概率图作为输入图像的二通道，从而将2D网络快速分割结果用于指导3D模型分割.

本研究的鼻咽癌患者MR图像中，肿瘤区域相对于整幅图像而言，占比很小，即肿瘤区域面积远小于非肿瘤区域面积，因此采用Dice损失^[13]作为基本损失函数，Dice定义如下：

(5) $ {\rm{Dice}} = \frac{{2|P \! \cap G|}}{{|P| + |G|}}. $

式中：P为模型的输出结果，G为真实标签，即人工勾画的肿瘤轮廓.

Dice损失定义为

(6) $ L = 1{\rm{ - Dice}}. $

模型总损失包括两部分：2D-ResUNet的损失L_2D以及2D+3D融合层损失L_F，模型总损失为这2个损失的加权和：

(7) $ {L_{\rm{Total}}} = \alpha {L_{\rm{2D}}} + {L_{\rm{F}}}. $

式中： $\alpha $为2D-ResUNet损失的权重，本研究中设为0.5，表示模型更关注最终输出的损失.

共收集421名NPC患者T1W、T2W、T1C三种脑结构MR图像，采用1.5T GE MEDICAL SYSTEMS水平位扫描。所有患者同种序列图像的成像参数相同，其中，T1W图像的参数为TR=666 ms，TE=9.17 ms，DFOV=20 0 mm，层厚为5 mm，共32层；T1C图像的参数为TR=650 ms，TE=9.17 ms，DFOV=200 mm，层厚为5 mm，共32层；T2W图像的参数为TR=3 100 ms，TE=82.51 ms，DFOV=200 mm，层厚为5 mm，共32层，所有3D图像水平位切片分辨率均为512×512. 由2名经验丰富的临床医生对T2W模态图像肿瘤区域进行手工勾画，确定分割标签. 从所有421名被试人中随机选择346名患者的图像作为训练集，剩余75名患者的图像作为独立测试集，训练集和测试集的信息如表2所示.

采用Dice系数、豪斯多夫距离（Hausdorff distance，HD）以及面积差占比（percentage of area difference，PAD）作为模型效果的评估指标. HD以及PAD定义如下：

(8) ${\rm{HD}}(P,G) = \mathop {\max }\limits_{a \in P}\;\{ \mathop {\max }\limits_{b \in G}\; [d(a,b)]\}. $

式中： $d(a,b)$为欧氏距离.

(9) ${\rm{PAD = }}\frac{{{|P -{ G|}}}}{G}.$

HD与PAD的值越小，表明模型分割结果与人工勾画结果越接近，说明模型性能越好.

首先利用单模态2D-ResUNet和3D-ResUNet以及2D+3D融合层对T1W、T2W和T1C三种模态数据分别构建单模态多维融合模型，分别为T1W-MDF、T2W-MDF、T1C-MDF；然后，将T1W、T2W和T1C两两分别组合，即将T1W+T2W、T1W+T1C和T2W+T1C共3种双模态图像分别作为3D图像的2个通道，对3D图像进行重组后，将2种模态图像作为多模态2D-ResUNet各编码器的输入，将2D-ResUNet输出概率图与2种模态3D输入图像合并后输入到3D-ResUNet中，再通过2D+3D融合层构建两模态多维信息融合分割模型，分别为T1W+T2W-MDF、T1W+T1C-MDF、T2W+T1C-MDF；最后，将T1W、T2W和T1C三种模态图像分别作为3D图像的三通道，对3D图像进行重组后，将各模态图像作为多模态2D-ResUNet各编码器的输入，而将多模态2D-ResUNet输出概率图与多模态3D输入图像合并后输入到3D-ResUNet中，再通过2D+3D融合层构建多模态多维信息融合（multi-modality multi-dimension fusion, MMMDF）分割模型.

采用Tensorflow^[14]软件库构建模型，训练阶段的初始学习率为0.001，每4轮乘以0.9进行衰减，采用Adam优化器进行优化. 显卡为NIVIDIA Titan XP GPU，单模态多维融合模型、双模态多维融合模型与多模态多维融合模型的训练时间分别为23 h 37 min、30 h 24 min以及34 h 47 min. 测试阶段3类模型中每个患者的测试时间分别约为13、18以及22 s.

使用相同的训练集和测试集，对文献[8]、[10]、[15]中的算法进行训练与测试，将结果与本文提出的MMMDF方法进行对比，如表3所示. 本文不同模型结果的箱型图如图3所示.

由表2及图3可知，MMMDF模型的分割性能优于以T1W、T2W和T1C任一单模态作为输入的模型，并优于任一两模态多维融合分割模型. 将MMMDF分割结果与3种两模态多维融合分割结果进行统计检验，结果显示，dice系数、面积差占比以及豪斯多夫距离的P值均小于0.05，表明MMMDF模型的分割性能比任一两模态多维融合模型性能更优，且具有统计性差异. 采用多模态作为2D-ResUNet的输入，在网络卷积过程中，将不同模态的信息进行融合：另外，采用多模态图像作为3D-ResUNet的输入，将不同模态间的信息再次进行融合，因此融合模型的效果优于单模态输入的效果. 此外，由于采用多维度模型串联的方式，同时结合了2D模型拟合速度快和3D模型特征学习充分的优点，最终在测试集上得到了较优的分割效果.

由表2可知，相比于其他文献中提到的方法，本文提出的MMMDF分割模型在NPC分割性能上有显著的提高. 这是由于这些文献中的研究仅采用简单的2D编码-解码网络，或基于此对网络结构进行微调，同时仅采用单种模态图像.

如图4所示为T1W-MDF、T2W-MDF、T1C-MDF以及MMMDF模型得到的部分分割结果. 如图5所示为相应的的三维表面图. 由图4和5可知，MMMDF模型对于不同形态、体积、不同区域上的肿瘤均表现出优良的分割效果. 在图4和5中，每一行表示不同患者，第1列表示输入图像；第2~8列分别表示T1W-MDF、T2W-MDF、T1C-MDF、文献[8]、文献[10]、文献[15]的方法以及MMMDF的结果. 图4的第1列为原图像，第2~8列为第1列窗口区域放大图，2条线分别表示人工勾画肿瘤区域（金标准）和模型分割结果.

为了实现鼻咽癌患者病灶的准确分割，本研究收集了NPC患者头颈部T1W、T2W和T1C三种模态MR图像，建立了一种新的多模态多维信息融合深度学习分割模型. 实验结果表明，与单种模态多维融合模型及公开发表的其他鼻咽癌病灶分割方法相比，多模态多维信息融合模型能更准确地定位病灶，显著提高分割效果. 本文提出的方法能够有效地实现NPC肿瘤的自动准确定位，为NPC的诊断与治疗提供了客观依据，提高了诊治效率与水平.

本研究主要有以下4点不足：1）采集了421名患者的3种模态MR图像，患者样本数量少，继续增加样本量，会有利于提高分割模型的泛化能力；2）MR图像的分辨率较低，层厚为5mm，空间结构信息不连续或部分缺失，一方面，给训练样本的肿瘤区域勾画带来困难，另一方面，不利于3D分割模型有效利用图像的层面拓扑结构信息；3）部分面积非常小的肿瘤区域，模型的分割性能较低，需要重点研究如何利用相邻层的肿瘤区域信息，改山分割结果；4）不同模态图像，部分出现位置差异。保证空间位置一致将有助于提高模型的可靠性与适应性。