内、中膜厚度是临床上用于评价动脉粥样硬化发展程度的主要指标。目前，基于 B 超图像测量内、中膜厚度通常由专业医生手动标记内、中膜边界来实现，过程繁琐耗时，人为影响因素多。本文提出了一种基于高斯混合模型（GMM）的聚类灰度阈值法，以检测 B 超图像中颈动脉内、中膜厚度。首先基于 GMM 对颈动脉图像灰度聚类，然后用灰度阈值法检测血管壁内、中膜的分界，最后测量二者的厚度。与直接使用灰度阈值法的测量技术相比，颈动脉 B 超图像的聚类解决了内、中膜灰度边界模糊的问题，从而提高了灰度阈值法的稳定性与检测精度。本研究选取 120 例健康颈动脉临床试验数据，以两名专家分别手动精细测量 4 次的内、中膜厚度的均值作为参考值，最终研究结果显示，经 GMM 聚类后估计的内、中膜厚度的归一化均方根误差（NRMSE）分别为 0.104 7 ± 0.076 2 和 0.097 4 ± 0.068 3；与直接进行灰度阈值估计的结果相比，NRMSE 的均值分别减小 19.6% 和 22.4%，表明本文所提方法测量精度有所提高；标准差分别减小 17.0% 和 21.7%，表明所提方法稳定性增加。综上，本文方法有助于动脉粥样硬化等血管疾病的早期诊断和病程监测。
引用本文： 戚贵玲, 何冰冰, 张榆锋, 李支尧, 莫鸿, 程洁. 基于高斯混合模型聚类的 B 超图像颈动脉内膜和中膜厚度检测. 生物医学工程学杂志, 2020, 37(6): 1080-1088, 1094. doi: 10.7507/1001-5515.201906027 复制
引言
心脑血管动脉硬化患病率高，由此引发的脑中风、冠心病等心脑血管病具有高致残率及高死亡率的特点，严重威胁人类健康。世界卫生组织调查显示，全球每年死于心血管疾病的人数高达 1 700 万，占世界人口死亡率的 31%，居各种死因之首[1]。研究表明，颈动脉粥样硬化导致血管壁增厚和粥样斑块形成，最早累及内膜，是引起冠心病、脑卒中等心血管疾病的首要原因[2]。目前，颈动脉内、中膜厚度（intima-media thickness，IMT）是临床上用于定量评价动脉粥样硬化发展程度的主要指标[3-4]。准确测量 IMT 对于心血管疾病的早期诊断以及疾病监测具有重要意义。超声成像技术能够完整呈现颈动脉的内、中膜结构，具有无创、低成本和可重复等优点，为测量 IMT 提供了可靠途径[5]。基于 B 超图像测量 IMT 通常由专业医生手动标记颈动脉内、中膜边界来实现，但过程繁琐、耗时且高度依赖医生的从业经验及主观评价。
随着高分辨、多功能超声诊断仪和计算机智能化信息处理的不断发展，大量非人工测量的新技术相继被提出，并应用于检测 B 超图像中的 IMT，其中包括动态规划、活动轮廓模型、统计建模以及深度学习等[6-16]。刘一学等[13]提出使用支持向量机测量 IMT，首先根据训练样本训练支持向量机得出分类模型，将感兴趣区的像素分为边界点和非边界点，再对边界点分类为管腔—内膜边界点和中膜—外膜边界点，采用启发式搜索对分类结果进行甄别，去除错分类的像素点，最终测量 IMT。孙萍等[14]基于卷积神经网络（convolutional neural networks，CNN）提出了一种全自动分割算法，使用 CNN 从图像中识别包含内、中膜信息的图像感兴趣区，并采用堆栈式自编码器对感兴趣区像素点进行分类，最后通过曲线拟合完成厚度测量。Savaş等[15]开发了一种具有多个隐藏层的深度学习策略，使用 CNN 对 IMT 进行分类。上述人工智能的方法显著提高了 IMT 的检测精度，同时具备良好的灵敏度和特异性，解决了专家手工测量 IMT 效率低下的问题，也避免了人工测量的主观性和可变性。然而，上述方法均需使用大量手工测量结果作为训练样本，训练样本的代表性和准确度与算法性能密切相关。Xiao 等[16]提出了一种基于马尔科夫随机场模型自动检测 IMT 的方法，该方法从 B 超图像中提取到含有内、中膜感兴趣区的子图像，经双边滤波后转换成一系列的同样宽度、非重叠的二进制图像片段，使用霍夫变换和迭代条件模式算法来分割出管腔—内膜的边界和中膜—外膜的边界，计算两条边界绝对距离的均值得到 IMT。结果表明，该方法能够准确提取到 IMT，是一种寻求最优解的高效计算方法，但该方法在 B 超图像有噪声或干扰过多的情况下检测效果不理想。上述算法应用于检测内膜与中膜的整体厚度，而最新研究表明，除了内膜与中膜的整体厚度，血管壁内膜和中膜各自的厚度能够更精确地表征和评估冠心病、心肌梗塞、中风等动脉粥样硬化程度[17-18]。
近年来，上述心血管疾病发展引起的颈动脉内膜和中膜各自厚度的变化吸引了人们的关注。与整体厚度的检测方法相比，检测内膜和中膜各自厚度的研究难点在于提取内、中膜边界，目前相关提取方法主要有边缘提取和深度学习两类技术。Bae 等[19]提出根据 B 超图像中颈动脉内膜和中膜超声回声强度不同导致所呈现的亮暗回波线变化这一情况，可用电子卡尺测量每个像素并估计内膜和中膜的厚度，改进了传统边缘检测方法。为了研究超声图像对比度对测量结果的影响，Macioch 等[20]对加入造影剂前后的颈动脉内膜和中膜厚度进行测量，得出管壁近端的内膜和中膜厚度在低对比度时低于高对比度以及手动测量的结果，表明 B 超图像对比度增强有利于内膜和中膜厚度的检测。Loizou 等[21]改进了活动轮廓模型，使用贪吃蛇（snake）算法自动检测超声图像颈动脉管壁的内膜和中膜轮廓以估计其厚度。该方法假设内膜和中膜两个成像区域灰度不变，将求解两个同质区域问题转化为求解一个能量方程的最小值问题，但该方法仅适用于 B 超图像中动脉壁轮廓边缘较明显的情况。上述边缘提取测量方法的性能依赖于在 B 超图像中内、中膜边界的灰度梯度信息，但易受斑点噪声影响。在理想的超声图像中，内、中膜边界显示清晰，灰度级呈现垂直跃迁变化。而临床上，实际的 B 超图像受超声斑点、混响伪影、低分辨率等不利因素的影响，加之内、中膜自身厚度薄（正常内膜的厚度约为 0.1～0.3 mm，占管壁厚度的 1/6），使内、中膜边界存在模糊、灰度不均匀等问题。袁绍锋等[22]利用 CNN 方法，学习原始血管内超声图像与所对应手动分割图像间的映射，预测内膜、中膜及外膜的概率图，实现医学图像语义分割，进而采用形态学闭、开操作平滑内膜和中外膜边界。CNN 方法对噪声干扰容忍度高，有效改善了内、中膜边界模糊的问题，显著提高了内膜和中膜厚度的测量精度。但该方法使用大量的手动分割结果作为训练样本，训练样本的数量及准确性与算法性能密切相关。综上所述，由于超声斑点噪声对图像的影响，手动分割内膜和中膜边界是一项难度较高的工作。
无需训练的图像聚类技术有望改善 B 超图像质量，解决内、中膜边界模糊的问题。图像聚类是运用一些准则函数将图像中灰度、色彩、纹理、结构等按照类内属性相同或相似进行分类，类间属性区别最大的区域合并，使得变化不连续的图像边界呈现为特征目标一致性较好的连续轮廓，从而获得目标边缘明显的图像[23-24]。由于针对整个特征空间，图像聚类更容易把握全局信息，受局部噪声的影响较小，改善了常规边缘检测方法的缺陷，目前在医学图像分割、配准等方面已成为一大研究热点。选择合适的概率模型描述颈动脉超声图像的灰度信息对图像聚类的效果至关重要。不同组织超声回声分布的研究结果表明，瑞利分布、莱斯分布等单一分布只能描述一种同质均匀组织，不适用于包括血液、三层膜以及周围组织等复杂组织结构的颈动脉 B 超图像，也难以描述三层膜随着动脉粥样硬化病程发展的病变情况[25-27]。对于结构复杂、类型多样的生物组织，混合分布通用性更强，泛化性更好。为了对比各类混合分布描述血管超声散斑分布的适用性，Mohana[28]使用伽马混合模型、瑞利混合模型和 Nakagami 混合模型估计组织的概率分布，结果表明伽马混合模型优于其他两种混合模型。在此研究的基础上，柴五一等[29]进一步对比了高斯混合模型（gaussian mixture model，GMM）和伽马混合模型，结果表明 GMM 更符合钙化斑块和正常血管区域的超声斑点概率分布。
为此，本文提出一种基于 GMM 的聚类灰度阈值法，以检测颈动脉内膜和中膜各自的厚度。假设将颈动脉 B 超图像中具有相同像素点特征值的每个区域看作单高斯分布，整体图像即可用 GMM 分布来描述。首先对内膜和中膜边界灰度变化不明显的 B 超图像采用 GMM 聚类处理以提高图像对比度，选择各像素点的灰度为其特征，用最大期望（expectation-maximization，EM）算法估计内膜和中膜特征参数，得到颈动脉血管内膜和中膜边界清晰的灰度图像，然后使用灰度阈值法自动检测血管壁内膜和中膜的分界位置，估计出内膜和中膜厚度。最后，为验证本文方法的可靠性与有效性，以两名专家手动精细测量的结果为参考值，分别使用直接阈值法和 GMM 聚类阈值法测量 120 例健康受试者颈动脉的内膜和中膜的厚度，并对测量结果进行误差分析。本文研究使用 GMM 对 B 超图像进行聚类可改善血管壁内、中膜边界模糊的问题，能够提高灰度阈值法检测颈动脉内、中膜厚度的估计精度，今后将有助于动脉粥样硬化等血管疾病的早期诊断和病程监测。
1 方法
1.1 GMM 及其参数估计
GMM 分布由 K 个高斯成分线性组合而成。对于观测数据集 
中的单个采样 xi，设第 k 个高斯分量的概率密度函数
如式（1）所示：
 |
式中，参数集 
，
为均值，
为标准差。则高斯混合分布的概率密度函数 
，如式（2）所示：
 |
式中，
是第 k 个高斯分量的权重，表示各混合成分的先验概率；
为各混合成分的参数矢量，属于某个参数空间。定义整个观测集数据的对数似然函数 
，如式（3）所示：
 |
其中 n 为样本总数。则混合模型参数 
，如式（4）所示：
 |
对式（3）的估计通常用 EM 算法，该算法是一种迭代算法，每次迭代由两步组成，第一步求期望，第二步求极大值[30]。在迭代开始之前，需定义隐随机变量 Z = {Zi}，使得 Zi = k 时，xi 属于第 k 个分布，并初始化参数 
、
和 
，如式（5）所示：
 |
其中，nk表示第k类中的采样数量。
第一步求期望的过程中，利用贝叶斯理论结合隐随机变量 Z，计算每个样本分别来源于高斯分布类别 k 的先验概率，如式（6）所示：
 |
第二步求极大值的过程中，对 U 函数进行最大化求解，估计下一次迭代的各分布模型参数 
、
和 
，如式（7）所示：
 |
重复（6）、（7）直到 
（
为收敛误差）时，迭代结束。GMM 模型参数，如式（8）至（11）所示：
 |
 |
 |
 |
其中，
、
和 
分别表示第 (t + 1) 次迭代中第 k 类的均值、方差以及所占权重。
1.2 基于 GMM 聚类的颈动脉 B 超图像处理
GMM 图像聚类是利用图像灰度特征集合的相似性准则来对图像进行分组，进而把图像不同目标划分为不相交的区域。人体颈动脉的管腔、管壁内膜、中膜以及周围组织等不同对象的声阻抗有显著差异，因此在 B 超图像中这些不同对象的超声斑点分布也不同。将每个对象的斑点分布各看作一个高斯分布，整个 B 超图像即为一种各类数据集中的 GMM。对 B 超图像进行 GMM 聚类处理会有助于准确测量颈动脉内膜及中膜各自的厚度，具体流程如图 1 所示。
图1
GMM 聚类阈值法测量颈动脉内、中膜各自厚度的流 程图
Figure1.
Flow chart of measuring the thickness of intima and media carotid artery by GMM clustering threshold method
给定一幅颈动脉血管 B 超图像，灰度空间像素值作为 GMM 输入样本；然后根据图像中对象个数设定分类数 k；对于每个样本 xi，如式（5）所示，初始化分类参数；如式（6）所示，计算样本最可能来自于类别 k 的概率；根据设定的最大迭代次数及收敛条件 
，如式（6）和（7）所示，完成每一类的特征参数估计。最终找到一组最优的参数，以产生 B 超图像灰度值的样本数据。最后将聚类结果转换到图像空间，得到聚类处理的颈动脉 B 超图像。
1.3 算法时间复杂度分析
在本文提出的 GMM 聚类阈值法中，对颈动脉超声图像进行 GMM 聚类处理后，使用灰度阈值测量内膜和中膜厚度再计算复杂度与直接灰度阈值法相同。因此，本文算法的运行时间主要受 EM 算法的影响。EM 算法的时间复杂度由求期望与求极大值两步的时间复杂度，以及收敛所需迭代次数共同决定[31]。设 L 是 EM 算法的迭代次数，k 是聚类个数，n 是超声图像中的像素总数，算法时间复杂度为 O(Lkn)。研究表明使用 EM 算法估计 GMM 模型参数时，其迭代次数L是样本总数 n 与聚类个数 k 的线性函数，则 EM 算法的时间复杂度是 O(k2n2)[32]。在本文中，聚类个数k预设为 4，因此，本文算法的时间复杂度为 O(n2)。关于空间复杂度，算法在迭代过程中需要一定空间存储中间变量，除此之外，无其他内存消耗。
1.4 颈动脉管壁灰度阈值三层膜检测
基于边缘的灰度阈值检测常用于图像目标测量、特征提取和分析，也用于测量颈动脉三层膜厚度[33-35]。在 B 超图像中，颈动脉管腔及管壁内、中以及外膜超声回声分布依次呈现为暗—亮—暗—亮的变化。如图 2 所示，右图为左侧原始 B 超图像中黄色虚线扫描线上红色箭头指示范围的灰度变化曲线，曲线上的最大峰值点和次大峰值点分别代表了血管壁的外膜和内膜中心。为了确定血管壁内、中膜的边界，根据内膜的灰度峰值 Tm，设定一个百分比例作为测量内、中膜边界的灰度阈值 T（图 2 右图中黑色虚线）。在血管壁灰度变化曲线中，灰度值与 T 相等的深度位置分别记为 y1、y2、y3，则内膜厚度为 Vi = y2 − y1，中膜厚度为 Vm = y3 − y2。基于灰度阈值法检测颈动脉三层膜精度对阈值的选取有较强的依赖性，目前对合适阈值的选择还没有一个客观的共识。
图2
灰度阈值法检测颈动脉血管壁三层膜示意图
Figure2.
Schematic diagram of three layers of carotid blood vessel wall detected by gray threshold method
2 数据来源与性能评价
本研究中所使用的 120 例颈动脉临床试验数据来源于 120 名健康人，其中 73 名男性、47 名女性，年龄范围：19～49 岁，由昆明医科大学第三附属医院提供并授权使用。数据采集试验经昆明医科大学第三附属医院伦理委员会批准，试验前所有受试对象都已被告知试验相关信息并自愿签署知情同意书。使用超声系统（Aplio 500，Toshiba Inc.，日本）及线阵探头（PLT-805AT，Toshiba Inc.，日本），中心频率为 5 MHz，频带范围为 5～12 MHz，室温 20～25 ℃，整个数据采集环境保持安静；受试者仰卧静息 10～15 min，头稍偏向一侧。利用脉冲超声以 30 Hz 的帧频扫描颈总动脉，声束与血管垂直。
为了评估 GMM 对颈动脉 B 超图像聚类处理的性能，本研究对比了 GMM 聚类前后的 B 超图像及血管壁三层膜的灰度变化曲线，分析血管壁边缘对比度及清晰度。本文分别使用直接阈值法和 GMM 聚类阈值法检测不同阈值下的内、中膜厚度，以分析 GMM 聚类处理对阈值法稳定性的改进。另外，考虑到专家手动精细分割测量是目前临床上评估内、中膜厚度测量精度的金标准，本文以两名专家手动精细分割 4 次结果的均值作为参考值，计算直接阈值法和 GMM 聚类阈值法对内、中膜厚度检测结果的归一化均方根误差，定量分析 GMM 聚类处理对阈值法检测精度的影响。试验在 4 GB 内存及中央处理器（i5-5200U，Intel Inc.，美国）上利用数据分析软件 MATLAB R2016a（MathWorks.，美国）编程实现。
3 实验结果与分析
3.1 GMM 聚类前后的 B 超图像对比
利用 GMM 算法对 120 幅 B 超图像进行聚类，算法参数设置为：类别k = 4，收敛误差 
；算法的极限终止条件设置为：最大类别为 10，最大迭代次数为 800。随机选取一位 26 岁男性颈动脉中段和一位 32 岁女性颈动脉近分叉处为例，如图 3～图 6 所示，对比原始 B 超图像和 GMM 聚类 B 超图像，并在从左到右 3 个不同位置（图 3 和图 5 中的黄色虚线所示）的血管壁成像区域（图 3 和图 5 中的红线所示）分析 GMM 聚类引起的灰度变化。
图3
26 岁男性颈动脉 GMM 处理前后 B 超图像
Figure3.
Ultrasound B-mode images before and after GMM clustering in the carotid artery of a 26-year-old male
图4
26 岁男性颈动脉 GMM 处理前后不同位置灰度变化曲线
Figure4.
Gray-scale curves at different positions before and after GMM clustering in the carotid artery of a 26-year-old male
图5
32 岁女性颈动脉 GMM 处理前后 B 超图像
Figure5.
Ultrasound B-mode images before and after GMM clustering in the carotid artery of a 32-year-old female
图6
32 岁女性颈动脉 GMM 处理前后不同位置灰度变化曲线
Figure6.
Gray-scale curves at different positions before and after GMM clustering in the carotid artery of a 32-year-old female
如图 3 所示，相比原始 B 超图像，GMM 聚类处理后的图像整体变化趋势更显著，各组织间超声回声分布差异明显，颈动脉内、中、外膜边缘划分更清晰。如图 4 所示，从血管壁三层膜灰度变化曲线看，不同位置的曲线呈现一致性的变化，表明 GMM 聚类处理一致性较好。对比同一位置 GMM 处理前后的灰度变化曲线，GMM 聚类结果在内膜和中膜、中膜和外膜灰度变化没有过渡状态，呈现垂直跃迁性突变；而在内膜、中膜各自的区域内灰度分布更均匀，变化显著。上述特征能够为后续灰度阈值法检测内、中膜厚度奠定良好的基础。
如图 5 和图 6 所示为另一例 32 岁健康女性颈动脉近分叉处 B 超原始图像及 3 个不同位置（黄色虚线）的血管壁（红线）灰度变化曲线，由图 5、图 6 可见，经 GMM 聚类处理后的结果整体呈现与图 3、图 4 所示 26 岁男性案例一致的变化趋势和特征。
3.2 基于阈值法的内、中膜测量结果
对图 3 和图 5 所示的两例 GMM 聚类处理前后的 B 超图像各随机选取 20 个位置，采用不同阈值测量颈动脉管壁内、中膜厚度。
如图 7 所示，是上述随机选取的 26 岁男性的测量结果，使用其原始图像检测的内、中膜厚度均值在 0.24～0.11 mm 之间随阈值的增大而减小，标准差呈随机变化；而 GMM 聚类图像的测量结果中，当阈值大于 20，内、中膜厚度的测量均值分别为（0.172
0.016）mm 和（0.221
0.019）mm。与原始图像相比，GMM 聚类阈值法的测量标准差减小。
图7
26 岁男性颈动脉内膜和中膜 GMM 聚类前后不同阈值测量结果
Figure7.
Measurement results of different thresholds before and after GMM clustering in the carotid intima and media of the 26-year-old male
如图 8 所示，32 岁女性颈动脉内膜和中膜 GMM 聚类前后不同阈值测量结果也呈现相似特征，表明 GMM 对 B 超图像的聚类处理避免了预设阈值对测量结果造成的不良影响，提高了阈值法检测内、中膜厚度的稳定性。
图8
32 岁女性颈动脉内膜和中膜 GMM 聚类前后不同阈值测量结果比较
Figure8.
Measurement results of different thresholds before and after GMM clustering in the carotid intima and media of the 32-year-old female
针对 120 幅原始 B 超图像，两名专家对原始 B 超图像手动精细测量内膜和中膜厚度 4 次，对 4 次测量结果的均值做线性回归分析。如图 9 所示，内、中膜测量结果拟合的直线斜率
均近似于 1，截距
较小，说明两名专家手动测量的结果一致性较好。以两名专家手动测量结果的均值作为参考值，如表 1 所示给出了直接阈值法和 GMM 聚类阈值法检测内、中膜厚度的归一化均方根误差（normalized root mean square error，NRMSE）。GMM 聚类阈值法估计的内、中膜厚度的 NRMSE 分别为 0.104 7 ± 0.076 2 和 0.097 4 ± 0.068 3，与直接阈值法估计的结果相比，NRMSE 的均值分别减小 19.6% 和 22.4%，标准差分别减小 17.0% 和 21.7%，这表明 GMM 聚类处理提高了灰度阈值法的检测精度（P < 0.05）。
表1
GMM 聚类前后阈值法检测内、中膜厚度的 NRMSE
Table1.
NRMSEs of intima and media thicknesses measure ments using the threshold method before and after the clustering of GMM.
图9
专家手动精细分割测量内膜和中膜厚度的线性回归分析
Figure9.
Linear regression analysis of intima and media thickness measured by expert manual fine segmentation
基于数据分析软件 MATLAB R2016a（MathWorks., 美国）的软件平台，以及 4 GB 内存、中央处理器（i5-5200U, Intel Inc., 美国）的硬件条件，对于 120 幅像素为 960 × 720 的 B 超图像，GMM 聚类处理的平均运行时间是 6.76 s。如果使用 C 语言实现并优化算法程序，可以进一步缩短算法运行时间，提高本文 GMM 聚类阈值法的临床实用性。
4 结论与讨论
本文提出一种 GMM 聚类阈值法，基于 GMM 对 B 超图像进行聚类处理以提高灰度阈值法检测颈动脉内、中膜厚度的估计精度和稳定性。通过对临床 B 超图像不同组织基于 GMM 进行聚类处理，可改善血管壁内、中膜边界模糊的问题，从而精确分割颈动脉内、中膜，有助于灰度阈值法检测内、中膜厚度。120 幅颈动脉 B 超图像的检测结果表明，与直接阈值法相比，内、中膜测量结果的 NRMSE 的均值分别减小了 19.6% 和 22.4%，标准差分别减小 17.0% 和 21.7%，因此，本文方法获得的结果精度更高。同时，在不同阈值下获得的结果具有较好的一致性，有助于降低预设阈值引入的测量误差，提高临床医生对相关疾病的诊断和治疗水平。然而本文方法也存在一定不足，没有系统研究图像噪声对检测结果的影响，研究主要针对正常血管，对于粥样硬化的病变血管无涉及，且局限于对超声图像斑点进行分类处理的情况，因此，通过拓展模型类别开展检测结果的对比分析将是下一步研究的重点。
利益冲突声明：本文全体作者均声明不存在利益冲突。
引言
心脑血管动脉硬化患病率高，由此引发的脑中风、冠心病等心脑血管病具有高致残率及高死亡率的特点，严重威胁人类健康。世界卫生组织调查显示，全球每年死于心血管疾病的人数高达 1 700 万，占世界人口死亡率的 31%，居各种死因之首[1]。研究表明，颈动脉粥样硬化导致血管壁增厚和粥样斑块形成，最早累及内膜，是引起冠心病、脑卒中等心血管疾病的首要原因[2]。目前，颈动脉内、中膜厚度（intima-media thickness，IMT）是临床上用于定量评价动脉粥样硬化发展程度的主要指标[3-4]。准确测量 IMT 对于心血管疾病的早期诊断以及疾病监测具有重要意义。超声成像技术能够完整呈现颈动脉的内、中膜结构，具有无创、低成本和可重复等优点，为测量 IMT 提供了可靠途径[5]。基于 B 超图像测量 IMT 通常由专业医生手动标记颈动脉内、中膜边界来实现，但过程繁琐、耗时且高度依赖医生的从业经验及主观评价。
随着高分辨、多功能超声诊断仪和计算机智能化信息处理的不断发展，大量非人工测量的新技术相继被提出，并应用于检测 B 超图像中的 IMT，其中包括动态规划、活动轮廓模型、统计建模以及深度学习等[6-16]。刘一学等[13]提出使用支持向量机测量 IMT，首先根据训练样本训练支持向量机得出分类模型，将感兴趣区的像素分为边界点和非边界点，再对边界点分类为管腔—内膜边界点和中膜—外膜边界点，采用启发式搜索对分类结果进行甄别，去除错分类的像素点，最终测量 IMT。孙萍等[14]基于卷积神经网络（convolutional neural networks，CNN）提出了一种全自动分割算法，使用 CNN 从图像中识别包含内、中膜信息的图像感兴趣区，并采用堆栈式自编码器对感兴趣区像素点进行分类，最后通过曲线拟合完成厚度测量。Savaş等[15]开发了一种具有多个隐藏层的深度学习策略，使用 CNN 对 IMT 进行分类。上述人工智能的方法显著提高了 IMT 的检测精度，同时具备良好的灵敏度和特异性，解决了专家手工测量 IMT 效率低下的问题，也避免了人工测量的主观性和可变性。然而，上述方法均需使用大量手工测量结果作为训练样本，训练样本的代表性和准确度与算法性能密切相关。Xiao 等[16]提出了一种基于马尔科夫随机场模型自动检测 IMT 的方法，该方法从 B 超图像中提取到含有内、中膜感兴趣区的子图像，经双边滤波后转换成一系列的同样宽度、非重叠的二进制图像片段，使用霍夫变换和迭代条件模式算法来分割出管腔—内膜的边界和中膜—外膜的边界，计算两条边界绝对距离的均值得到 IMT。结果表明，该方法能够准确提取到 IMT，是一种寻求最优解的高效计算方法，但该方法在 B 超图像有噪声或干扰过多的情况下检测效果不理想。上述算法应用于检测内膜与中膜的整体厚度，而最新研究表明，除了内膜与中膜的整体厚度，血管壁内膜和中膜各自的厚度能够更精确地表征和评估冠心病、心肌梗塞、中风等动脉粥样硬化程度[17-18]。
近年来，上述心血管疾病发展引起的颈动脉内膜和中膜各自厚度的变化吸引了人们的关注。与整体厚度的检测方法相比，检测内膜和中膜各自厚度的研究难点在于提取内、中膜边界，目前相关提取方法主要有边缘提取和深度学习两类技术。Bae 等[19]提出根据 B 超图像中颈动脉内膜和中膜超声回声强度不同导致所呈现的亮暗回波线变化这一情况，可用电子卡尺测量每个像素并估计内膜和中膜的厚度，改进了传统边缘检测方法。为了研究超声图像对比度对测量结果的影响，Macioch 等[20]对加入造影剂前后的颈动脉内膜和中膜厚度进行测量，得出管壁近端的内膜和中膜厚度在低对比度时低于高对比度以及手动测量的结果，表明 B 超图像对比度增强有利于内膜和中膜厚度的检测。Loizou 等[21]改进了活动轮廓模型，使用贪吃蛇（snake）算法自动检测超声图像颈动脉管壁的内膜和中膜轮廓以估计其厚度。该方法假设内膜和中膜两个成像区域灰度不变，将求解两个同质区域问题转化为求解一个能量方程的最小值问题，但该方法仅适用于 B 超图像中动脉壁轮廓边缘较明显的情况。上述边缘提取测量方法的性能依赖于在 B 超图像中内、中膜边界的灰度梯度信息，但易受斑点噪声影响。在理想的超声图像中，内、中膜边界显示清晰，灰度级呈现垂直跃迁变化。而临床上，实际的 B 超图像受超声斑点、混响伪影、低分辨率等不利因素的影响，加之内、中膜自身厚度薄（正常内膜的厚度约为 0.1～0.3 mm，占管壁厚度的 1/6），使内、中膜边界存在模糊、灰度不均匀等问题。袁绍锋等[22]利用 CNN 方法，学习原始血管内超声图像与所对应手动分割图像间的映射，预测内膜、中膜及外膜的概率图，实现医学图像语义分割，进而采用形态学闭、开操作平滑内膜和中外膜边界。CNN 方法对噪声干扰容忍度高，有效改善了内、中膜边界模糊的问题，显著提高了内膜和中膜厚度的测量精度。但该方法使用大量的手动分割结果作为训练样本，训练样本的数量及准确性与算法性能密切相关。综上所述，由于超声斑点噪声对图像的影响，手动分割内膜和中膜边界是一项难度较高的工作。
无需训练的图像聚类技术有望改善 B 超图像质量，解决内、中膜边界模糊的问题。图像聚类是运用一些准则函数将图像中灰度、色彩、纹理、结构等按照类内属性相同或相似进行分类，类间属性区别最大的区域合并，使得变化不连续的图像边界呈现为特征目标一致性较好的连续轮廓，从而获得目标边缘明显的图像[23-24]。由于针对整个特征空间，图像聚类更容易把握全局信息，受局部噪声的影响较小，改善了常规边缘检测方法的缺陷，目前在医学图像分割、配准等方面已成为一大研究热点。选择合适的概率模型描述颈动脉超声图像的灰度信息对图像聚类的效果至关重要。不同组织超声回声分布的研究结果表明，瑞利分布、莱斯分布等单一分布只能描述一种同质均匀组织，不适用于包括血液、三层膜以及周围组织等复杂组织结构的颈动脉 B 超图像，也难以描述三层膜随着动脉粥样硬化病程发展的病变情况[25-27]。对于结构复杂、类型多样的生物组织，混合分布通用性更强，泛化性更好。为了对比各类混合分布描述血管超声散斑分布的适用性，Mohana[28]使用伽马混合模型、瑞利混合模型和 Nakagami 混合模型估计组织的概率分布，结果表明伽马混合模型优于其他两种混合模型。在此研究的基础上，柴五一等[29]进一步对比了高斯混合模型（gaussian mixture model，GMM）和伽马混合模型，结果表明 GMM 更符合钙化斑块和正常血管区域的超声斑点概率分布。
为此，本文提出一种基于 GMM 的聚类灰度阈值法，以检测颈动脉内膜和中膜各自的厚度。假设将颈动脉 B 超图像中具有相同像素点特征值的每个区域看作单高斯分布，整体图像即可用 GMM 分布来描述。首先对内膜和中膜边界灰度变化不明显的 B 超图像采用 GMM 聚类处理以提高图像对比度，选择各像素点的灰度为其特征，用最大期望（expectation-maximization，EM）算法估计内膜和中膜特征参数，得到颈动脉血管内膜和中膜边界清晰的灰度图像，然后使用灰度阈值法自动检测血管壁内膜和中膜的分界位置，估计出内膜和中膜厚度。最后，为验证本文方法的可靠性与有效性，以两名专家手动精细测量的结果为参考值，分别使用直接阈值法和 GMM 聚类阈值法测量 120 例健康受试者颈动脉的内膜和中膜的厚度，并对测量结果进行误差分析。本文研究使用 GMM 对 B 超图像进行聚类可改善血管壁内、中膜边界模糊的问题，能够提高灰度阈值法检测颈动脉内、中膜厚度的估计精度，今后将有助于动脉粥样硬化等血管疾病的早期诊断和病程监测。
1 方法
1.1 GMM 及其参数估计
GMM 分布由 K 个高斯成分线性组合而成。对于观测数据集 中的单个采样 xi，设第 k 个高斯分量的概率密度函数如式（1）所示：
|
式中，参数集 ， 为均值， 为标准差。则高斯混合分布的概率密度函数 ，如式（2）所示：
|
式中， 是第 k 个高斯分量的权重，表示各混合成分的先验概率； 为各混合成分的参数矢量，属于某个参数空间。定义整个观测集数据的对数似然函数 ，如式（3）所示：
|
其中 n 为样本总数。则混合模型参数 ，如式（4）所示：
|
对式（3）的估计通常用 EM 算法，该算法是一种迭代算法，每次迭代由两步组成，第一步求期望，第二步求极大值[30]。在迭代开始之前，需定义隐随机变量 Z = {Zi}，使得 Zi = k 时，xi 属于第 k 个分布，并初始化参数 、 和 ，如式（5）所示：
|
其中，nk表示第k类中的采样数量。
第一步求期望的过程中，利用贝叶斯理论结合隐随机变量 Z，计算每个样本分别来源于高斯分布类别 k 的先验概率，如式（6）所示：
|
第二步求极大值的过程中，对 U 函数进行最大化求解，估计下一次迭代的各分布模型参数 、 和 ，如式（7）所示：
|
重复（6）、（7）直到 （为收敛误差）时，迭代结束。GMM 模型参数，如式（8）至（11）所示：
|
|
|
|
其中，、 和 分别表示第 (t + 1) 次迭代中第 k 类的均值、方差以及所占权重。
1.2 基于 GMM 聚类的颈动脉 B 超图像处理
GMM 图像聚类是利用图像灰度特征集合的相似性准则来对图像进行分组，进而把图像不同目标划分为不相交的区域。人体颈动脉的管腔、管壁内膜、中膜以及周围组织等不同对象的声阻抗有显著差异，因此在 B 超图像中这些不同对象的超声斑点分布也不同。将每个对象的斑点分布各看作一个高斯分布，整个 B 超图像即为一种各类数据集中的 GMM。对 B 超图像进行 GMM 聚类处理会有助于准确测量颈动脉内膜及中膜各自的厚度，具体流程如图 1 所示。
图1
GMM 聚类阈值法测量颈动脉内、中膜各自厚度的流 程图
Figure1.
Flow chart of measuring the thickness of intima and media carotid artery by GMM clustering threshold method
给定一幅颈动脉血管 B 超图像，灰度空间像素值作为 GMM 输入样本；然后根据图像中对象个数设定分类数 k；对于每个样本 xi，如式（5）所示，初始化分类参数；如式（6）所示，计算样本最可能来自于类别 k 的概率；根据设定的最大迭代次数及收敛条件 ，如式（6）和（7）所示，完成每一类的特征参数估计。最终找到一组最优的参数，以产生 B 超图像灰度值的样本数据。最后将聚类结果转换到图像空间，得到聚类处理的颈动脉 B 超图像。
1.3 算法时间复杂度分析
在本文提出的 GMM 聚类阈值法中，对颈动脉超声图像进行 GMM 聚类处理后，使用灰度阈值测量内膜和中膜厚度再计算复杂度与直接灰度阈值法相同。因此，本文算法的运行时间主要受 EM 算法的影响。EM 算法的时间复杂度由求期望与求极大值两步的时间复杂度，以及收敛所需迭代次数共同决定[31]。设 L 是 EM 算法的迭代次数，k 是聚类个数，n 是超声图像中的像素总数，算法时间复杂度为 O(Lkn)。研究表明使用 EM 算法估计 GMM 模型参数时，其迭代次数L是样本总数 n 与聚类个数 k 的线性函数，则 EM 算法的时间复杂度是 O(k2n2)[32]。在本文中，聚类个数k预设为 4，因此，本文算法的时间复杂度为 O(n2)。关于空间复杂度，算法在迭代过程中需要一定空间存储中间变量，除此之外，无其他内存消耗。
1.4 颈动脉管壁灰度阈值三层膜检测
基于边缘的灰度阈值检测常用于图像目标测量、特征提取和分析，也用于测量颈动脉三层膜厚度[33-35]。在 B 超图像中，颈动脉管腔及管壁内、中以及外膜超声回声分布依次呈现为暗—亮—暗—亮的变化。如图 2 所示，右图为左侧原始 B 超图像中黄色虚线扫描线上红色箭头指示范围的灰度变化曲线，曲线上的最大峰值点和次大峰值点分别代表了血管壁的外膜和内膜中心。为了确定血管壁内、中膜的边界，根据内膜的灰度峰值 Tm，设定一个百分比例作为测量内、中膜边界的灰度阈值 T（图 2 右图中黑色虚线）。在血管壁灰度变化曲线中，灰度值与 T 相等的深度位置分别记为 y1、y2、y3，则内膜厚度为 Vi = y2 − y1，中膜厚度为 Vm = y3 − y2。基于灰度阈值法检测颈动脉三层膜精度对阈值的选取有较强的依赖性，目前对合适阈值的选择还没有一个客观的共识。
图2
灰度阈值法检测颈动脉血管壁三层膜示意图
Figure2.
Schematic diagram of three layers of carotid blood vessel wall detected by gray threshold method
2 数据来源与性能评价
本研究中所使用的 120 例颈动脉临床试验数据来源于 120 名健康人，其中 73 名男性、47 名女性，年龄范围：19～49 岁，由昆明医科大学第三附属医院提供并授权使用。数据采集试验经昆明医科大学第三附属医院伦理委员会批准，试验前所有受试对象都已被告知试验相关信息并自愿签署知情同意书。使用超声系统（Aplio 500，Toshiba Inc.，日本）及线阵探头（PLT-805AT，Toshiba Inc.，日本），中心频率为 5 MHz，频带范围为 5～12 MHz，室温 20～25 ℃，整个数据采集环境保持安静；受试者仰卧静息 10～15 min，头稍偏向一侧。利用脉冲超声以 30 Hz 的帧频扫描颈总动脉，声束与血管垂直。
为了评估 GMM 对颈动脉 B 超图像聚类处理的性能，本研究对比了 GMM 聚类前后的 B 超图像及血管壁三层膜的灰度变化曲线，分析血管壁边缘对比度及清晰度。本文分别使用直接阈值法和 GMM 聚类阈值法检测不同阈值下的内、中膜厚度，以分析 GMM 聚类处理对阈值法稳定性的改进。另外，考虑到专家手动精细分割测量是目前临床上评估内、中膜厚度测量精度的金标准，本文以两名专家手动精细分割 4 次结果的均值作为参考值，计算直接阈值法和 GMM 聚类阈值法对内、中膜厚度检测结果的归一化均方根误差，定量分析 GMM 聚类处理对阈值法检测精度的影响。试验在 4 GB 内存及中央处理器（i5-5200U，Intel Inc.，美国）上利用数据分析软件 MATLAB R2016a（MathWorks.，美国）编程实现。
3 实验结果与分析
3.1 GMM 聚类前后的 B 超图像对比
利用 GMM 算法对 120 幅 B 超图像进行聚类，算法参数设置为：类别k = 4，收敛误差 ；算法的极限终止条件设置为：最大类别为 10，最大迭代次数为 800。随机选取一位 26 岁男性颈动脉中段和一位 32 岁女性颈动脉近分叉处为例，如图 3～图 6 所示，对比原始 B 超图像和 GMM 聚类 B 超图像，并在从左到右 3 个不同位置（图 3 和图 5 中的黄色虚线所示）的血管壁成像区域（图 3 和图 5 中的红线所示）分析 GMM 聚类引起的灰度变化。
图3
26 岁男性颈动脉 GMM 处理前后 B 超图像
Figure3.
Ultrasound B-mode images before and after GMM clustering in the carotid artery of a 26-year-old male
图4
26 岁男性颈动脉 GMM 处理前后不同位置灰度变化曲线
Figure4.
Gray-scale curves at different positions before and after GMM clustering in the carotid artery of a 26-year-old male
图5
32 岁女性颈动脉 GMM 处理前后 B 超图像
Figure5.
Ultrasound B-mode images before and after GMM clustering in the carotid artery of a 32-year-old female
图6
32 岁女性颈动脉 GMM 处理前后不同位置灰度变化曲线
Figure6.
Gray-scale curves at different positions before and after GMM clustering in the carotid artery of a 32-year-old female
如图 3 所示，相比原始 B 超图像，GMM 聚类处理后的图像整体变化趋势更显著，各组织间超声回声分布差异明显，颈动脉内、中、外膜边缘划分更清晰。如图 4 所示，从血管壁三层膜灰度变化曲线看，不同位置的曲线呈现一致性的变化，表明 GMM 聚类处理一致性较好。对比同一位置 GMM 处理前后的灰度变化曲线，GMM 聚类结果在内膜和中膜、中膜和外膜灰度变化没有过渡状态，呈现垂直跃迁性突变；而在内膜、中膜各自的区域内灰度分布更均匀，变化显著。上述特征能够为后续灰度阈值法检测内、中膜厚度奠定良好的基础。
如图 5 和图 6 所示为另一例 32 岁健康女性颈动脉近分叉处 B 超原始图像及 3 个不同位置（黄色虚线）的血管壁（红线）灰度变化曲线，由图 5、图 6 可见，经 GMM 聚类处理后的结果整体呈现与图 3、图 4 所示 26 岁男性案例一致的变化趋势和特征。
3.2 基于阈值法的内、中膜测量结果
对图 3 和图 5 所示的两例 GMM 聚类处理前后的 B 超图像各随机选取 20 个位置，采用不同阈值测量颈动脉管壁内、中膜厚度。
如图 7 所示，是上述随机选取的 26 岁男性的测量结果，使用其原始图像检测的内、中膜厚度均值在 0.24～0.11 mm 之间随阈值的增大而减小，标准差呈随机变化；而 GMM 聚类图像的测量结果中，当阈值大于 20，内、中膜厚度的测量均值分别为（0.1720.016）mm 和（0.2210.019）mm。与原始图像相比，GMM 聚类阈值法的测量标准差减小。
图7
26 岁男性颈动脉内膜和中膜 GMM 聚类前后不同阈值测量结果
Figure7.
Measurement results of different thresholds before and after GMM clustering in the carotid intima and media of the 26-year-old male
如图 8 所示，32 岁女性颈动脉内膜和中膜 GMM 聚类前后不同阈值测量结果也呈现相似特征，表明 GMM 对 B 超图像的聚类处理避免了预设阈值对测量结果造成的不良影响，提高了阈值法检测内、中膜厚度的稳定性。
图8
32 岁女性颈动脉内膜和中膜 GMM 聚类前后不同阈值测量结果比较
Figure8.
Measurement results of different thresholds before and after GMM clustering in the carotid intima and media of the 32-year-old female
针对 120 幅原始 B 超图像，两名专家对原始 B 超图像手动精细测量内膜和中膜厚度 4 次，对 4 次测量结果的均值做线性回归分析。如图 9 所示，内、中膜测量结果拟合的直线斜率均近似于 1，截距较小，说明两名专家手动测量的结果一致性较好。以两名专家手动测量结果的均值作为参考值，如表 1 所示给出了直接阈值法和 GMM 聚类阈值法检测内、中膜厚度的归一化均方根误差（normalized root mean square error，NRMSE）。GMM 聚类阈值法估计的内、中膜厚度的 NRMSE 分别为 0.104 7 ± 0.076 2 和 0.097 4 ± 0.068 3，与直接阈值法估计的结果相比，NRMSE 的均值分别减小 19.6% 和 22.4%，标准差分别减小 17.0% 和 21.7%，这表明 GMM 聚类处理提高了灰度阈值法的检测精度（P < 0.05）。
表1
GMM 聚类前后阈值法检测内、中膜厚度的 NRMSE
Table1.
NRMSEs of intima and media thicknesses measure ments using the threshold method before and after the clustering of GMM.
图9
专家手动精细分割测量内膜和中膜厚度的线性回归分析
Figure9.
Linear regression analysis of intima and media thickness measured by expert manual fine segmentation
基于数据分析软件 MATLAB R2016a（MathWorks., 美国）的软件平台，以及 4 GB 内存、中央处理器（i5-5200U, Intel Inc., 美国）的硬件条件，对于 120 幅像素为 960 × 720 的 B 超图像，GMM 聚类处理的平均运行时间是 6.76 s。如果使用 C 语言实现并优化算法程序，可以进一步缩短算法运行时间，提高本文 GMM 聚类阈值法的临床实用性。
4 结论与讨论
本文提出一种 GMM 聚类阈值法，基于 GMM 对 B 超图像进行聚类处理以提高灰度阈值法检测颈动脉内、中膜厚度的估计精度和稳定性。通过对临床 B 超图像不同组织基于 GMM 进行聚类处理，可改善血管壁内、中膜边界模糊的问题，从而精确分割颈动脉内、中膜，有助于灰度阈值法检测内、中膜厚度。120 幅颈动脉 B 超图像的检测结果表明，与直接阈值法相比，内、中膜测量结果的 NRMSE 的均值分别减小了 19.6% 和 22.4%，标准差分别减小 17.0% 和 21.7%，因此，本文方法获得的结果精度更高。同时，在不同阈值下获得的结果具有较好的一致性，有助于降低预设阈值引入的测量误差，提高临床医生对相关疾病的诊断和治疗水平。然而本文方法也存在一定不足，没有系统研究图像噪声对检测结果的影响，研究主要针对正常血管，对于粥样硬化的病变血管无涉及，且局限于对超声图像斑点进行分类处理的情况，因此，通过拓展模型类别开展检测结果的对比分析将是下一步研究的重点。
利益冲突声明：本文全体作者均声明不存在利益冲突。
