內、中膜厚度是臨床上用于評價動脈粥樣硬化發展程度的主要指標。目前，基于 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%，因此，本文方法獲得的結果精度更高。同時，在不同閾值下獲得的結果具有較好的一致性，有助于降低預設閾值引入的測量誤差，提高臨床醫生對相關疾病的診斷和治療水平。然而本文方法也存在一定不足，沒有系統研究圖像噪聲對檢測結果的影響，研究主要針對正常血管，對于粥樣硬化的病變血管無涉及，且局限于對超聲圖像斑點進行分類處理的情況，因此，通過拓展模型類別開展檢測結果的對比分析將是下一步研究的重點。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
