根據對頸動脈斑塊幾何形狀、組織特征及聲學特性的分析，提出基于FIELD II頸動脈斑塊血管的超聲仿真方法，并利用MATLAB中的圖形用戶界面設計開發超聲仿真系統。以基于元球的隱式曲面造型技術對斑塊血管進行三維幾何建模，以空間位置伽馬隨機分布的散射點對其進行組織建模，超聲仿真并建立仿真系統。對系統仿真的超聲圖像進行統計分析，與模型預設值及臨床B超圖像進行比較，結果充分驗證了本文仿真方法及仿真系統的有效性。
引用本文： 姚瑞晗, 張榆鋒, 施繼紅, 胡曉, 楊麗春, 高蓮, 章克信, 李支堯. 頸動脈斑塊的超聲仿真及系統實現. 生物醫學工程學雜志, 2016, 33(6): 1095-1102. doi: 10.7507/1001-5515.20160175 復制
版權信息： ?四川大學華西醫院華西期刊社《生物醫學工程學雜志》版權所有，未經授權不得轉載、改編
引言
心腦血管疾病是世界上三大致死疾病之一，其中頸動脈粥樣硬化是最常見、也是后果嚴重的一種血管病變[1]，其主要特征為頸動脈某些部位的內膜下脂質沉積，并伴有平滑肌細胞和纖維基質成分的增殖，逐步發展形成動脈粥樣硬化性斑塊，造成血管狹窄甚至栓塞。若斑塊破裂，其碎片會隨血流流入腦部，使大腦供血受阻，引發中風等突發性疾病，甚至導致死亡[2]。頸動脈粥樣硬化早期無癥狀，準確診斷疾病的發生發展是及時采取治療措施、防止病情惡化的關鍵，也是后續采取對應治療手段的重要依據。超聲技術因其無損、無輻射、實時、便捷等特點已成為頸動脈粥樣硬化常規檢查與診斷的重要手段[3]。但在實際臨床操作中，B超圖像固有的斑點噪聲，以及成像中存在的環暈偽影、導絲偽影等干擾會影響診斷結果的準確性[4]。隨著計算機信息技術的發展，針對這些干擾因素采取各種消去算法以提高圖像質量，并通過研究特征自動提取和分類技術在一定程度上提高了診斷的可靠性[5]。借助數學建模的方法，研究斑塊形狀、組織特征及聲學特性，開發計算機超聲仿真方法，能夠為這些信號處理和特征提取方法研究提供客觀的數據源和評價依據。此外，對斑塊臨床特征建立數據庫，利用計算機圖形界面建立頸動脈斑塊超聲仿真系統，能夠為醫護人員培訓等應用提供有效手段。
近年來一些研究提出了頸動脈斑塊的計算機超聲仿真方法。Groot等[6]提出了一個仿真血管內組織的超聲模擬器。該方法首先將血管切片以獲得病理組織圖像，并作預處理提取圖像邊緣作為頸動脈斑塊的形狀輪廓。再采用位置隨機分布的散射點建立組織模型，并用射線跟蹤(ray tracing)法模擬聚焦聲束產生掃描線，并以斑塊形狀輪廓排列，插值生成最后圖像。這種方法獲得的斑塊形狀輪廓盡管逼真，但卻缺乏一般性的幾何特征描述。Ramírez等[7]提出的差分背向散射截面法則通過建立一個簡單的物理模型來仿真斑塊形狀，并給定空間分布的單個散射點和背向散射截面來表示斑塊的組織模型，最后模擬換能器獲得回波信號，仿真相應的超聲圖像。但這種方法沒有考慮到斑塊的不同形狀、組織類型以及同一斑塊不同組成部分在超聲圖像上的不同影像特征。國內開展的研究中，汪源源團隊[8]提出了一種基于極坐標圖像生成模型的血管內超聲圖像仿真方法。該方法先在笛卡爾坐標系中確定血管內腔輪廓和斑塊脂質區及鈣化區輪廓，再將確定的各輪廓變換到極坐標下，分別對粥樣硬化斑塊的纖維、脂質和鈣化組成部分設定組織分布，并在組織模型中引入環暈及導絲偽影進行靜態和動態圖像仿真。此外，胡曉等[9]還對頸動脈三層膜管壁結構進行了超聲仿真。以廣義圓柱為基礎，通過控制圓柱的直徑、面積、體積等相關形態學參數得到所需血管結構；根據組織散射點分布的形狀和密度與超聲回聲斑點分布的變化規律，對臨床圖像進行統計分析，確定組織散射點分布，利用計算機產生超聲射頻回波信號并得到B超圖像。但這兩種方法都沒有對斑塊形狀進行系統的分類，也沒有深入開展計算機仿真系統的開發研究。
本文分析總結了臨床頸動脈斑塊的形態、組織結構和回聲斑點特征，提出基于元球造型技術對管狀多分枝組織進行造型的方法建立頸動脈及斑塊的三維幾何模型。采用位置伽馬隨機分布組織散射點模型，研究形狀和密度參數與超聲回聲斑點分布的變化規律。通過對臨床頸動脈斑塊組織特征的分析，選擇合適的參數，利用FIELD II模擬超聲探頭計算超聲射頻回波信號，進而仿真各種類型的頸動脈斑塊B超圖像。在此基礎上，對臨床頸動脈斑塊的幾何結構、組織尺寸和分布特征分類建立數據庫，利用MATLAB中的圖形用戶界面(Graphical User Interface,GUI)設計開發了超聲仿真系統。
1 頸動脈斑塊的分類及特征
對頸動脈斑塊進行超聲計算機建模仿真中，斑塊的形態和組織結構特征的準確描述是關鍵。為此，首先對臨床頸動脈斑塊進行了分析和總結，根據幾何形狀和組織結構對斑塊分類，并歸納了各種類型斑塊的特征。
1.1 斑塊形狀的分類及特征
頸動脈壁上斑塊的形成隨機性較大，不可能有兩個形狀完全相同的斑塊，同時各類型斑塊之間也沒有明確的定義，這就給斑塊形狀的分類造成了困難。Beaumont等[10]為解決這一問題，對俄亥俄州克利夫蘭醫學中心的42例患者的血管內超聲成像數據進行了分析研究，應用K-均值模糊聚類的方法將斑塊形狀分為峰谷(peak-valley)、升序 (ascending)、降序(descending)及彌散(diffuse)四類，如圖 1所示。峰谷是那些表現出顯著突起與峰谷的結構；升序是那些坡度有上升或增加的結構；降序則是坡度下降或減少的結構；彌散即為斜率恒定或具有一定彌漫性的結構。
圖1
斑塊形狀分類
(a)峰谷；(b)升序；(c)降序；(d)彌散
Figure1. Classification of plaque geometry(a) peak-valley; (b) ascending; (c) descending; (d) diffuse
1.2 斑塊組織結構分類及特征
通過動脈內膜切除術取出頸動脈斑塊作為組織標本進行研究[11]，將斑塊分為硬斑、混合斑和軟斑，并描述了三種類型斑塊的構成、形態和超聲回聲表現[12]。為此，我們將頸動脈斑塊組織的構成、幾何尺寸、回聲特征和穩定性表現進行了總結(見表 1)，作為后續斑塊建模、超聲仿真及系統研發的依據。
表1
斑塊類型及特征
Table1.
Plaque types and characteristics
2 頸動脈斑塊的超聲建模及仿真
基于Field II進行頸動脈斑塊的超聲仿真，首先需要建立頸動脈斑塊的三維幾何模型，并用刻畫組織特點的散射點作三維填充，然后調用軟件平臺中的xdc_linear_array函數創建孔徑，xdc_impulse函數和xdc_excitation函數設置孔徑脈沖響應和激勵脈沖，通過xdc_apodization函數建立變跡時間線，xdc_center_focus函數設置聚焦中心點，最后用calc_scat函數計算射頻回波信號，并經包絡提取、插值和對數壓縮最終生成B超圖像。
2.1 頸動脈斑塊的三維幾何模型
以基于元球的隱式曲面造型技術對頸動脈斑塊進行三維幾何建模。該技術采用具有等勢場值的元球點集來定義等勢曲面，生成的曲面光滑，宜于表示可形變的物體，適合對人體、動物器官和液體造型。在一個元球系統中，每個元球可以有不同的勢函數，我們用fi 表示第i個元球Pi的勢函數，通常被定義為一個多項式函數。常用的有分段二次多項式、四次多項式及六次多項式等[13-14]。由于頸動脈斑塊屬于形狀不規則且表面光滑的組織結構，為了使仿真形狀更為逼真，并考慮到計算的復雜度，本文采用四次多項式作為元球勢函數：
| ${{f}_{i}}\left( r \right)=\left\{ \begin{align} &{{\left( 1-{{\left( \frac{r}{{{R}_{i}}} \right)}^{2}} \right)}^{2}},0\le r\le R \\ &0,r\ge {{R}_{i}} \\ \end{align} \right.$ |
其中，Ri是第i個元球的有效半徑，r是空間某一點到第i個元球中心的歐氏距離，定義為：
| $r=\sqrt{{{(x-{{x}_{i}})}^{2}}+{{(y-{{y}_{i}})}^{2}}+{{(z-{{z}_{i}})}^{2}}}$ |
其中，(xi,yi,zi)為第i個元球的球心，(x,y,z)為空間某一點。對于由n個元球所構成的系統，對應的等勢面為滿足以下方程的點集：
| $f\left( x,y,z \right)=\sum\limits_{i=1}^{n}{{}}{{q}_{i}}{{f}_{i}}({{x}_{i}},{{y}_{i}},{{z}_{i}})-T=0$ |
其中fi(xi,yi,zi)為第i個元球的勢函數，T為閾值，即所確定的表面的勢能。qi 為第i個元球的密度值或質量。構造模型時，通過調節元球球心、半徑來產生所需的不同勢函數，再通過式(3)判斷空間點是否滿足所對應的等勢面，并將滿足的點集連接以此來產生與所需的生物體形態相似的等勢面。
以10個元球為例，設元球的密度值為1，通過給定元球勢函數，產生一個管腔的三維結構，再根據狹窄度、纖維帽及血管壁厚度調節位于狹窄部位的球心坐標、半徑及勢能閾值，生成帶有峰谷類型斑塊的管腔三維結構，并根據等勢面獲得頸動脈斑塊的三維幾何模型，如圖 2所示。
圖2
使用元球模型構造頸動脈斑塊
(a)正常頸動脈；(b)有斑塊的頸動脈；(c)頸動脈斑塊三維幾何模型
Figure2. Carotid artery plaque constructed based on the metaball model(a) normal carotid artery; (b) carotid artery with a plaque; (c) three-dimensional geometric model of the carotid artery plaque
2.2 組織散射點模型及其與回聲斑點分布對應關系
頸動脈斑塊組織結構復雜，不同類型斑塊以及同一斑塊的不同組織之間超聲回聲分布特點不同，大致可分為均勻高回聲區、少許透聲分布的高回聲區、透聲與高回聲混合區以及大型透聲區。研究認為，當組織特征為非確定性成分存在，超聲分辨單元內有少量隨機聚集分布的散射點，回聲斑點表現為大型透聲區的K分布型；組織特征為非確定性成分存在，分辨單元內有較多隨機分布的散射點，回聲為透聲區與高回聲區混合存在的瑞利分布；組織特征為確定性成分存在，分辨單元內有充分規律分布的散射點，對應為高回聲區的萊斯分布。Field II仿真中可通過設置散射點不同密度的空間分布模擬產生不同特點的超聲回聲分布。
目前所用的組織散射點模型主要有三種：Neyman-Scott模型、Gibbs-Markov模型和Cramblitt模型[15]。其中，Cramblitt模型因其能夠在不同密度條件下對散射點的空間位置分布從聚集到隨機再到規律連續可調，在組織超聲仿真中更具靈活性與逼真性。一維Cramblitt散射點模型定義為：
| $s\left( x \right)=\Sigma {{a}_{i}}\delta (x-{{X}_{i}})$ |
式中ai為服從正態分布的散射點強度，Xi表示散射點位置，x是一維實線上的位置變量。模型中的散射點間距服從伽馬分布：
| $f\left( d \right)=\frac{{{d}^{\alpha -1}}exp\left( -d/\beta \right)}{\Gamma \left( \alpha \right){{\beta }^{\alpha }}}\alpha ,\beta ,d>0$ |
式中di=|Xi-Xi+1|為散射點間距，α為形狀參數，β為尺度參數，該分布的均值和方差分別為d=αβ和σd2=αβ2=d2/α。因此，散射點分布的特征可通過密度參數(ρ=1/d)和形狀參數(α)連續調整。三維散射點模型是通過Hilbert空間填充曲線[16]對一維散射點進行三維映射實現的。
研究表明[9]，組織散射點模型中不同的密度參數ρ和形狀參數α取值對應回聲斑點不同的分布模型及參數：K分布模型、瑞利分布模型和萊斯分布模型，如下：
K分布：
| $f\left( A \right)=\frac{2b}{\Gamma \left( s \right)}{{\left( \frac{bA}{2} \right)}^{s}}{{K}_{v-1}}\left( bA \right)\text{ }v>0$ |
式中，Kv-1(·)是v-1階的第二類修正貝塞爾函數，Γ(s)是伽瑪函數，A為回聲振幅，b和s為模型參數。
瑞利分布：
| $f\left( A \right)=\frac{A}{{{r}^{2}}}exp(-\frac{{{A}^{2}}}{2{{r}^{2}}})$ |
式中A為回聲振幅，r為模型參數。
萊斯分布：
| $f\left( A \right)=\frac{A}{{{k}^{2}}}exp(-\frac{{{A}^{2}}+{{s}^{2}}}{2{{k}^{2}}}){{I}_{0}}(\frac{As}{{{k}^{2}}})\text{ }s\ge 0$ |
式中，I0(·)為0階第一類修正貝塞爾函數，A為回聲振幅，s和k為模型參數。
2.3 斑塊及管壁組織散射點模型參數的確定
為了獲得逼真的仿真效果，選擇昆明市第三人民醫院90例頸動脈粥樣硬化患者，其中硬斑、軟斑和混合斑各30例，男女比例平均，年齡為(50.5±6.8)歲，均獲得知情同意且符合倫理審查原則。患者取仰臥位，頸后墊枕后仰，頭偏向外側，采用SIEMENS ACUSON S2000彩超進行成像。成像參數：14L5淺表探頭，頻率14 MHz，視野范圍為15 cm×17 cm。對這些臨床圖像中各組織做手工精細分割，分別得到血管壁、纖維帽和脂質核的回聲圖像。對這些圖像作直方圖統計和極大似然擬合，將擬合的分布及參數與2.2節的回聲斑點模型及參數作比較，確定對應的組織散射點模型參數ρ及α。
2.4 基于GUI界面的仿真系統設計
為了給頸動脈斑塊的特征提取、可視化研究和醫護人員培訓提供有效的仿真手段，研究進一步利用MATLAB中的GUI設計開發了超聲仿真系統。MATLAB的GUI是由窗口、菜單、圖標、按鍵、對話框和文本等各種圖形對象組成的用戶與計算機進行信息交流的方式，特別適合圖文對象的交互。本文設計的GUI超聲仿真系統主要包括斑塊類型與幾何描述、組織散射點模型描述以及超聲仿真成像三個模塊。斑塊類型與幾何描述模塊根據選定類型生成纖維帽、血管壁和脂質核三部分對應的尺寸范圍，設定組織的具體幾何尺寸，并生成斑塊的幾何模型；組織散射點模型描述模塊包括選擇確定纖維帽、血管壁和脂質核三部分的組織散射點模型參數。模塊根據選定類型給出各組織對應散射點模型參數的可能取值范圍，用戶可自行搭配設定模型中的形狀和密度參數值；超聲仿真成像模塊基于Field II仿真軟件平臺設置聲學參數，生成射頻回波信號，并經包絡提取、插值和對數壓縮生成B超圖像。
3 實驗結果
根據2.2節給出的組織散射點模型與回聲斑點分布對應關系，表 2列出了ρ和α幾個典型取值對應的回聲分布及參數。
表2
散射點模型中不同ρ和α對應的回聲斑點分布及參數
Table2.
Speckle distributions and parameters corresponding to the different ρ and α in scatterer models
為了準確確定對斑塊各組織仿真的散射點模型參數，應用2.3節描述的方法對臨床圖像進行統計分析。圖 3~圖 5給出了其中一例的結果。圖 3(a)、圖 4(a)及圖 5(a)分別給出了混合斑、硬斑和軟斑臨床B超圖像。子圖(b)、(c)和(d)分別為血管壁、纖維帽和脂質核對應的直方圖、極大似然擬合及模型參數。可以看出，所有斑塊類型的血管壁及斑塊纖維帽的超聲回聲都屬萊斯分布(參數有差異)，而混合斑、硬斑及軟斑的脂質核分別符合瑞利、萊斯和K分布。對各類斑塊的30例B超圖像分別進行手工分割，作統計直方圖和極大似然擬合處理，結果表明同類斑塊各組織回聲斑點的分布都保持一致。但由于臨床上個體之間存在差異，對應的分布參數也有差別。將臨床統計結果進一步與表 2進行對比，得到組織散射點模型參數ρ及α，表 3列出了各參數取值的均值及標準差。仿真時各組織對應散射點模型的形狀和密度參數在表 3列出的可能范圍內隨機取值。
圖3
混合斑塊臨床B超圖及各類組織回聲斑點分布統計結果
(a)混合斑B超圖；(b)血管壁；(c)纖維帽；(d)脂質核
Figure3. Clinical B-mode ultrasound image of a mixed plaque and the statistical distributions for echoed speckles from different kinds of tissue(a) the image; (b) vessel wall; (c) fibrous cap; (d) lipid core
圖4
硬斑臨床B超圖及各類組織回聲斑點分布統計結果
(a)硬斑B超圖；(b)血管壁；(c)纖維帽；(d)脂質核
Figure4. Clinical B-mode ultrasound image of a hard plaque and the statistical distributions for echoed speckles from different kinds of tissue(a) the image; (b) vessel wall; (c) fibrous cap; (d) lipid core
圖5
軟斑臨床B超圖及各類組織回聲斑點分布統計結果
(a)軟斑B超圖；(b)血管壁；(c)纖維帽；(d)脂質核
Figure5. Clinical B-mode ultrasound image of a soft plaque and the statistical distributions for echoed speckles from different kinds of tissue(a) the image; (b) vessel wall; (c) fibrous cap; (d) lipid core
表3
頸動脈斑塊各組織散射點模型相關參數
Table3.
Parameters of scatterer models for different kinds of tissue in plaques
為了驗證仿真方法的有效性，實驗首先對不同類型和幾何形狀的斑塊仿真各30例。仿真主要參數包括：256個陣元的線性陣列探頭，中心頻率12 MHz，采樣頻率100 MHz，聲速為1 540 m/s，波長1.28×10-4 m，激活換能器陣元數64個，陣元寬度和高度分別為0.64×10-4 m及0.005 m，陣元間距為5×10-5 m，固定焦點為[0,0,0.02] m。圖 6分別給出了一例混合斑、硬斑和軟斑不同幾何形狀的B超仿真結果。其中(a)和(b)分別對應升序和彌散兩種形狀的混合斑，(c)對應峰谷形狀的軟斑，(d)對應降序形狀的硬斑。由于回聲的干涉和衍射，圖像中表現出不同程度的斑點分布，但血管壁、纖維帽以及脂質核各組織成分仍清晰可辨，符合超聲成像的規律。從圖像灰度上看，血管壁最亮，回聲強度最高，一致性最好；纖維帽的灰度略低于血管壁，回聲強度及一致性次之；與血管壁和纖維帽相比，脂質核回聲強度最低，其中硬斑脂質核亮度略高于混合斑，但二者都屬于回聲強度較高并含少許透聲分布，而軟斑的脂質核回聲強度較低，含有大量透聲區，亮度也相對較暗。這些回聲強度的特征與臨床超聲圖像表現一致。
圖6
B超仿真圖像
(a)混合斑(升序)；(b)混合斑(彌散)；(c)軟斑(峰谷)；(d)硬斑(降序)
Figure6. Simulated B-mode ulreasound images(a) mixed plaque (ascending); (b) mixed plaque (diffuse); (c) soft plaque (peak-valley); (d) hard plaque (descending)
對仿真的頸動脈斑塊B超圖像進行手工精細分割，分別對各類斑塊的血管壁、纖維帽以及脂質核進行回聲斑點直方圖統計及模型擬合、幾何尺寸測量以及誤差分析。回聲斑點分布的統計結果如表 4所列。與圖 3~圖 5的臨床統計結果比較，仿真回聲斑點的分布和參數值與臨床超聲的總體一致。其中，混合斑、硬斑及軟斑的血管壁和纖維帽都屬于萊斯分布(參數不同)，而混合斑、硬斑及軟斑的脂質核分別屬于瑞利、萊斯和K分布。斑塊中各組織的模型參數雖然存在一些偏差，但與臨床對應的結果大體一致，并處于給定參數的偏差范圍內。
表4
仿真的各類斑塊組織回聲斑點分布及參數
Table4.
Distributions and parameters of echoed speckles for simulated images of different kinds of plaque tissue
幾何尺寸測量及誤差統計結果如表 5所列。與預設值相比，混合斑、硬斑和軟斑的纖維帽相對誤差分別為5%、7.5%及10%，這是由于纖維帽厚度薄，而超聲側向分辨率低，超聲隨機斑點的影響較大，特別是對曲面形狀的纖維帽影響更顯著；而斑塊長度、管腔直徑以及血管壁厚度尺寸較大，超聲隨機斑點的影響較小，最大相對誤差僅為0.01%。混合斑、硬斑和軟斑的最大內膜厚度各不相同，其最大相對誤差分別為0.7%、0.5%及0.6%。
表5
仿真的各類斑塊組織幾何參數預設值及測量結果
Table5.
Preset values and measured results of geometrical parameters for simulated images of different kinds of plaque tissue
圖 7為仿真系統各模塊功能GUI界面。其中(a)為斑塊類型與幾何描述界面。該界面可根據需要選擇斑塊的類型與形狀，設定斑塊的幾何尺寸，生成相應的幾何模型。模塊中還可對生成的幾何模型進行放大、縮小、旋轉及尺寸測量等一系列操作，并可將模型保存下來。(b)為散射點模型描述界面。該界面根據表 3給出了散射點模型中不同的ρ和α取值及對應的回聲斑點分布及參數，用戶可在給定的取值范圍內自行選擇各組織散射點模型進行搭配，建立所需斑塊的散射點模型。(c)為超聲仿真成像界面。該界面基于Field II仿真軟件平臺，生成射頻回波信號，并經包絡提取、插值和對數壓縮生成B超圖像。
圖7
GUI系統的界面
(a)斑塊類型和形狀選擇；(b)散射點模型選擇；(c)生成的B超圖像
Figure7. GUI system interface(a) the selection of plaque types and shapes； (b) the selection of scatterer models； (c) the simulated B-mode ultrasound image
4 結論
本文應用FIELD II軟件對頸動脈斑塊進行超聲仿真，并利用MATLAB中的GUI結合斑塊的臨床特征建立了B超仿真系統。首先對斑塊的形狀尺寸以及特征進行了總結研究，再根據斑塊的回聲斑點特性以及臨床數據分析得到所需組織散射點模型，將模型代入FIELD II軟件平臺中得到超聲射頻回波信號生成超聲圖像，最后利用MATLAB中的GUI界面建立一個B超仿真系統，清晰地對斑塊B超成像過程進行模擬。對仿真圖像進行手工分割，分別對各組織回聲斑點和幾何尺寸分析測量，結果表明仿真圖像表現出的回聲特性能準確刻畫出頸動脈斑塊的組織特征，與臨床成像的結果一致；仿真圖像的幾何尺寸與預設值相比，纖維帽在5%~10%的誤差范圍內，其他組織在0.5%~0.7%的誤差范圍內，均符合超聲成像規律。建立的仿真系統能為頸動脈斑塊的特征提取和可視化研究提供有效的檢測數據源及仿真手段，為醫護人員培訓提供系統的仿真平臺。使用中先根據預設的斑塊類型、幾何尺寸及組織特征生成超聲圖像，受訓人員通過對圖像的分析獲取斑塊的特征，并與預設值相比較，提高醫護人員在臨床上對斑塊性質分析判斷的準確性。但本文方法仍有需要改進的地方，首先，斑塊組織模型只考慮了回聲分布特征，而沒有考慮組織的紋理特性；其次，沒有考慮脈動血壓的動力學因素對斑塊及血管壁的影響。針對以上問題進行研究是下一步工作的重點。
引言
心腦血管疾病是世界上三大致死疾病之一，其中頸動脈粥樣硬化是最常見、也是后果嚴重的一種血管病變[1]，其主要特征為頸動脈某些部位的內膜下脂質沉積，并伴有平滑肌細胞和纖維基質成分的增殖，逐步發展形成動脈粥樣硬化性斑塊，造成血管狹窄甚至栓塞。若斑塊破裂，其碎片會隨血流流入腦部，使大腦供血受阻，引發中風等突發性疾病，甚至導致死亡[2]。頸動脈粥樣硬化早期無癥狀，準確診斷疾病的發生發展是及時采取治療措施、防止病情惡化的關鍵，也是后續采取對應治療手段的重要依據。超聲技術因其無損、無輻射、實時、便捷等特點已成為頸動脈粥樣硬化常規檢查與診斷的重要手段[3]。但在實際臨床操作中，B超圖像固有的斑點噪聲，以及成像中存在的環暈偽影、導絲偽影等干擾會影響診斷結果的準確性[4]。隨著計算機信息技術的發展，針對這些干擾因素采取各種消去算法以提高圖像質量，并通過研究特征自動提取和分類技術在一定程度上提高了診斷的可靠性[5]。借助數學建模的方法，研究斑塊形狀、組織特征及聲學特性，開發計算機超聲仿真方法，能夠為這些信號處理和特征提取方法研究提供客觀的數據源和評價依據。此外，對斑塊臨床特征建立數據庫，利用計算機圖形界面建立頸動脈斑塊超聲仿真系統，能夠為醫護人員培訓等應用提供有效手段。
近年來一些研究提出了頸動脈斑塊的計算機超聲仿真方法。Groot等[6]提出了一個仿真血管內組織的超聲模擬器。該方法首先將血管切片以獲得病理組織圖像，并作預處理提取圖像邊緣作為頸動脈斑塊的形狀輪廓。再采用位置隨機分布的散射點建立組織模型，并用射線跟蹤(ray tracing)法模擬聚焦聲束產生掃描線，并以斑塊形狀輪廓排列，插值生成最后圖像。這種方法獲得的斑塊形狀輪廓盡管逼真，但卻缺乏一般性的幾何特征描述。Ramírez等[7]提出的差分背向散射截面法則通過建立一個簡單的物理模型來仿真斑塊形狀，并給定空間分布的單個散射點和背向散射截面來表示斑塊的組織模型，最后模擬換能器獲得回波信號，仿真相應的超聲圖像。但這種方法沒有考慮到斑塊的不同形狀、組織類型以及同一斑塊不同組成部分在超聲圖像上的不同影像特征。國內開展的研究中，汪源源團隊[8]提出了一種基于極坐標圖像生成模型的血管內超聲圖像仿真方法。該方法先在笛卡爾坐標系中確定血管內腔輪廓和斑塊脂質區及鈣化區輪廓，再將確定的各輪廓變換到極坐標下，分別對粥樣硬化斑塊的纖維、脂質和鈣化組成部分設定組織分布，并在組織模型中引入環暈及導絲偽影進行靜態和動態圖像仿真。此外，胡曉等[9]還對頸動脈三層膜管壁結構進行了超聲仿真。以廣義圓柱為基礎，通過控制圓柱的直徑、面積、體積等相關形態學參數得到所需血管結構；根據組織散射點分布的形狀和密度與超聲回聲斑點分布的變化規律，對臨床圖像進行統計分析，確定組織散射點分布，利用計算機產生超聲射頻回波信號并得到B超圖像。但這兩種方法都沒有對斑塊形狀進行系統的分類，也沒有深入開展計算機仿真系統的開發研究。
本文分析總結了臨床頸動脈斑塊的形態、組織結構和回聲斑點特征，提出基于元球造型技術對管狀多分枝組織進行造型的方法建立頸動脈及斑塊的三維幾何模型。采用位置伽馬隨機分布組織散射點模型，研究形狀和密度參數與超聲回聲斑點分布的變化規律。通過對臨床頸動脈斑塊組織特征的分析，選擇合適的參數，利用FIELD II模擬超聲探頭計算超聲射頻回波信號，進而仿真各種類型的頸動脈斑塊B超圖像。在此基礎上，對臨床頸動脈斑塊的幾何結構、組織尺寸和分布特征分類建立數據庫，利用MATLAB中的圖形用戶界面(Graphical User Interface,GUI)設計開發了超聲仿真系統。
1 頸動脈斑塊的分類及特征
對頸動脈斑塊進行超聲計算機建模仿真中，斑塊的形態和組織結構特征的準確描述是關鍵。為此，首先對臨床頸動脈斑塊進行了分析和總結，根據幾何形狀和組織結構對斑塊分類，并歸納了各種類型斑塊的特征。
1.1 斑塊形狀的分類及特征
頸動脈壁上斑塊的形成隨機性較大，不可能有兩個形狀完全相同的斑塊，同時各類型斑塊之間也沒有明確的定義，這就給斑塊形狀的分類造成了困難。Beaumont等[10]為解決這一問題，對俄亥俄州克利夫蘭醫學中心的42例患者的血管內超聲成像數據進行了分析研究，應用K-均值模糊聚類的方法將斑塊形狀分為峰谷(peak-valley)、升序 (ascending)、降序(descending)及彌散(diffuse)四類，如圖 1所示。峰谷是那些表現出顯著突起與峰谷的結構；升序是那些坡度有上升或增加的結構；降序則是坡度下降或減少的結構；彌散即為斜率恒定或具有一定彌漫性的結構。
圖1
斑塊形狀分類
(a)峰谷；(b)升序；(c)降序；(d)彌散
Figure1. Classification of plaque geometry(a) peak-valley; (b) ascending; (c) descending; (d) diffuse
1.2 斑塊組織結構分類及特征
通過動脈內膜切除術取出頸動脈斑塊作為組織標本進行研究[11]，將斑塊分為硬斑、混合斑和軟斑，并描述了三種類型斑塊的構成、形態和超聲回聲表現[12]。為此，我們將頸動脈斑塊組織的構成、幾何尺寸、回聲特征和穩定性表現進行了總結(見表 1)，作為后續斑塊建模、超聲仿真及系統研發的依據。
表1
斑塊類型及特征
Table1.
Plaque types and characteristics
2 頸動脈斑塊的超聲建模及仿真
基于Field II進行頸動脈斑塊的超聲仿真，首先需要建立頸動脈斑塊的三維幾何模型，并用刻畫組織特點的散射點作三維填充，然后調用軟件平臺中的xdc_linear_array函數創建孔徑，xdc_impulse函數和xdc_excitation函數設置孔徑脈沖響應和激勵脈沖，通過xdc_apodization函數建立變跡時間線，xdc_center_focus函數設置聚焦中心點，最后用calc_scat函數計算射頻回波信號，并經包絡提取、插值和對數壓縮最終生成B超圖像。
2.1 頸動脈斑塊的三維幾何模型
以基于元球的隱式曲面造型技術對頸動脈斑塊進行三維幾何建模。該技術采用具有等勢場值的元球點集來定義等勢曲面，生成的曲面光滑，宜于表示可形變的物體，適合對人體、動物器官和液體造型。在一個元球系統中，每個元球可以有不同的勢函數，我們用fi 表示第i個元球Pi的勢函數，通常被定義為一個多項式函數。常用的有分段二次多項式、四次多項式及六次多項式等[13-14]。由于頸動脈斑塊屬于形狀不規則且表面光滑的組織結構，為了使仿真形狀更為逼真，并考慮到計算的復雜度，本文采用四次多項式作為元球勢函數：
| ${{f}_{i}}\left( r \right)=\left\{ \begin{align} &{{\left( 1-{{\left( \frac{r}{{{R}_{i}}} \right)}^{2}} \right)}^{2}},0\le r\le R \\ &0,r\ge {{R}_{i}} \\ \end{align} \right.$ |
其中，Ri是第i個元球的有效半徑，r是空間某一點到第i個元球中心的歐氏距離，定義為：
| $r=\sqrt{{{(x-{{x}_{i}})}^{2}}+{{(y-{{y}_{i}})}^{2}}+{{(z-{{z}_{i}})}^{2}}}$ |
其中，(xi,yi,zi)為第i個元球的球心，(x,y,z)為空間某一點。對于由n個元球所構成的系統，對應的等勢面為滿足以下方程的點集：
| $f\left( x,y,z \right)=\sum\limits_{i=1}^{n}{{}}{{q}_{i}}{{f}_{i}}({{x}_{i}},{{y}_{i}},{{z}_{i}})-T=0$ |
其中fi(xi,yi,zi)為第i個元球的勢函數，T為閾值，即所確定的表面的勢能。qi 為第i個元球的密度值或質量。構造模型時，通過調節元球球心、半徑來產生所需的不同勢函數，再通過式(3)判斷空間點是否滿足所對應的等勢面，并將滿足的點集連接以此來產生與所需的生物體形態相似的等勢面。
以10個元球為例，設元球的密度值為1，通過給定元球勢函數，產生一個管腔的三維結構，再根據狹窄度、纖維帽及血管壁厚度調節位于狹窄部位的球心坐標、半徑及勢能閾值，生成帶有峰谷類型斑塊的管腔三維結構，并根據等勢面獲得頸動脈斑塊的三維幾何模型，如圖 2所示。
圖2
使用元球模型構造頸動脈斑塊
(a)正常頸動脈；(b)有斑塊的頸動脈；(c)頸動脈斑塊三維幾何模型
Figure2. Carotid artery plaque constructed based on the metaball model(a) normal carotid artery; (b) carotid artery with a plaque; (c) three-dimensional geometric model of the carotid artery plaque
2.2 組織散射點模型及其與回聲斑點分布對應關系
頸動脈斑塊組織結構復雜，不同類型斑塊以及同一斑塊的不同組織之間超聲回聲分布特點不同，大致可分為均勻高回聲區、少許透聲分布的高回聲區、透聲與高回聲混合區以及大型透聲區。研究認為，當組織特征為非確定性成分存在，超聲分辨單元內有少量隨機聚集分布的散射點，回聲斑點表現為大型透聲區的K分布型；組織特征為非確定性成分存在，分辨單元內有較多隨機分布的散射點，回聲為透聲區與高回聲區混合存在的瑞利分布；組織特征為確定性成分存在，分辨單元內有充分規律分布的散射點，對應為高回聲區的萊斯分布。Field II仿真中可通過設置散射點不同密度的空間分布模擬產生不同特點的超聲回聲分布。
目前所用的組織散射點模型主要有三種：Neyman-Scott模型、Gibbs-Markov模型和Cramblitt模型[15]。其中，Cramblitt模型因其能夠在不同密度條件下對散射點的空間位置分布從聚集到隨機再到規律連續可調，在組織超聲仿真中更具靈活性與逼真性。一維Cramblitt散射點模型定義為：
| $s\left( x \right)=\Sigma {{a}_{i}}\delta (x-{{X}_{i}})$ |
式中ai為服從正態分布的散射點強度，Xi表示散射點位置，x是一維實線上的位置變量。模型中的散射點間距服從伽馬分布：
| $f\left( d \right)=\frac{{{d}^{\alpha -1}}exp\left( -d/\beta \right)}{\Gamma \left( \alpha \right){{\beta }^{\alpha }}}\alpha ,\beta ,d>0$ |
式中di=|Xi-Xi+1|為散射點間距，α為形狀參數，β為尺度參數，該分布的均值和方差分別為d=αβ和σd2=αβ2=d2/α。因此，散射點分布的特征可通過密度參數(ρ=1/d)和形狀參數(α)連續調整。三維散射點模型是通過Hilbert空間填充曲線[16]對一維散射點進行三維映射實現的。
研究表明[9]，組織散射點模型中不同的密度參數ρ和形狀參數α取值對應回聲斑點不同的分布模型及參數：K分布模型、瑞利分布模型和萊斯分布模型，如下：
K分布：
| $f\left( A \right)=\frac{2b}{\Gamma \left( s \right)}{{\left( \frac{bA}{2} \right)}^{s}}{{K}_{v-1}}\left( bA \right)\text{ }v>0$ |
式中，Kv-1(·)是v-1階的第二類修正貝塞爾函數，Γ(s)是伽瑪函數，A為回聲振幅，b和s為模型參數。
瑞利分布：
| $f\left( A \right)=\frac{A}{{{r}^{2}}}exp(-\frac{{{A}^{2}}}{2{{r}^{2}}})$ |
式中A為回聲振幅，r為模型參數。
萊斯分布：
| $f\left( A \right)=\frac{A}{{{k}^{2}}}exp(-\frac{{{A}^{2}}+{{s}^{2}}}{2{{k}^{2}}}){{I}_{0}}(\frac{As}{{{k}^{2}}})\text{ }s\ge 0$ |
式中，I0(·)為0階第一類修正貝塞爾函數，A為回聲振幅，s和k為模型參數。
2.3 斑塊及管壁組織散射點模型參數的確定
為了獲得逼真的仿真效果，選擇昆明市第三人民醫院90例頸動脈粥樣硬化患者，其中硬斑、軟斑和混合斑各30例，男女比例平均，年齡為(50.5±6.8)歲，均獲得知情同意且符合倫理審查原則。患者取仰臥位，頸后墊枕后仰，頭偏向外側，采用SIEMENS ACUSON S2000彩超進行成像。成像參數：14L5淺表探頭，頻率14 MHz，視野范圍為15 cm×17 cm。對這些臨床圖像中各組織做手工精細分割，分別得到血管壁、纖維帽和脂質核的回聲圖像。對這些圖像作直方圖統計和極大似然擬合，將擬合的分布及參數與2.2節的回聲斑點模型及參數作比較，確定對應的組織散射點模型參數ρ及α。
2.4 基于GUI界面的仿真系統設計
為了給頸動脈斑塊的特征提取、可視化研究和醫護人員培訓提供有效的仿真手段，研究進一步利用MATLAB中的GUI設計開發了超聲仿真系統。MATLAB的GUI是由窗口、菜單、圖標、按鍵、對話框和文本等各種圖形對象組成的用戶與計算機進行信息交流的方式，特別適合圖文對象的交互。本文設計的GUI超聲仿真系統主要包括斑塊類型與幾何描述、組織散射點模型描述以及超聲仿真成像三個模塊。斑塊類型與幾何描述模塊根據選定類型生成纖維帽、血管壁和脂質核三部分對應的尺寸范圍，設定組織的具體幾何尺寸，并生成斑塊的幾何模型；組織散射點模型描述模塊包括選擇確定纖維帽、血管壁和脂質核三部分的組織散射點模型參數。模塊根據選定類型給出各組織對應散射點模型參數的可能取值范圍，用戶可自行搭配設定模型中的形狀和密度參數值；超聲仿真成像模塊基于Field II仿真軟件平臺設置聲學參數，生成射頻回波信號，并經包絡提取、插值和對數壓縮生成B超圖像。
3 實驗結果
根據2.2節給出的組織散射點模型與回聲斑點分布對應關系，表 2列出了ρ和α幾個典型取值對應的回聲分布及參數。
表2
散射點模型中不同ρ和α對應的回聲斑點分布及參數
Table2.
Speckle distributions and parameters corresponding to the different ρ and α in scatterer models
為了準確確定對斑塊各組織仿真的散射點模型參數，應用2.3節描述的方法對臨床圖像進行統計分析。圖 3~圖 5給出了其中一例的結果。圖 3(a)、圖 4(a)及圖 5(a)分別給出了混合斑、硬斑和軟斑臨床B超圖像。子圖(b)、(c)和(d)分別為血管壁、纖維帽和脂質核對應的直方圖、極大似然擬合及模型參數。可以看出，所有斑塊類型的血管壁及斑塊纖維帽的超聲回聲都屬萊斯分布(參數有差異)，而混合斑、硬斑及軟斑的脂質核分別符合瑞利、萊斯和K分布。對各類斑塊的30例B超圖像分別進行手工分割，作統計直方圖和極大似然擬合處理，結果表明同類斑塊各組織回聲斑點的分布都保持一致。但由于臨床上個體之間存在差異，對應的分布參數也有差別。將臨床統計結果進一步與表 2進行對比，得到組織散射點模型參數ρ及α，表 3列出了各參數取值的均值及標準差。仿真時各組織對應散射點模型的形狀和密度參數在表 3列出的可能范圍內隨機取值。
圖3
混合斑塊臨床B超圖及各類組織回聲斑點分布統計結果
(a)混合斑B超圖；(b)血管壁；(c)纖維帽；(d)脂質核
Figure3. Clinical B-mode ultrasound image of a mixed plaque and the statistical distributions for echoed speckles from different kinds of tissue(a) the image; (b) vessel wall; (c) fibrous cap; (d) lipid core
圖4
硬斑臨床B超圖及各類組織回聲斑點分布統計結果
(a)硬斑B超圖；(b)血管壁；(c)纖維帽；(d)脂質核
Figure4. Clinical B-mode ultrasound image of a hard plaque and the statistical distributions for echoed speckles from different kinds of tissue(a) the image; (b) vessel wall; (c) fibrous cap; (d) lipid core
圖5
軟斑臨床B超圖及各類組織回聲斑點分布統計結果
(a)軟斑B超圖；(b)血管壁；(c)纖維帽；(d)脂質核
Figure5. Clinical B-mode ultrasound image of a soft plaque and the statistical distributions for echoed speckles from different kinds of tissue(a) the image; (b) vessel wall; (c) fibrous cap; (d) lipid core
表3
頸動脈斑塊各組織散射點模型相關參數
Table3.
Parameters of scatterer models for different kinds of tissue in plaques
為了驗證仿真方法的有效性，實驗首先對不同類型和幾何形狀的斑塊仿真各30例。仿真主要參數包括：256個陣元的線性陣列探頭，中心頻率12 MHz，采樣頻率100 MHz，聲速為1 540 m/s，波長1.28×10-4 m，激活換能器陣元數64個，陣元寬度和高度分別為0.64×10-4 m及0.005 m，陣元間距為5×10-5 m，固定焦點為[0,0,0.02] m。圖 6分別給出了一例混合斑、硬斑和軟斑不同幾何形狀的B超仿真結果。其中(a)和(b)分別對應升序和彌散兩種形狀的混合斑，(c)對應峰谷形狀的軟斑，(d)對應降序形狀的硬斑。由于回聲的干涉和衍射，圖像中表現出不同程度的斑點分布，但血管壁、纖維帽以及脂質核各組織成分仍清晰可辨，符合超聲成像的規律。從圖像灰度上看，血管壁最亮，回聲強度最高，一致性最好；纖維帽的灰度略低于血管壁，回聲強度及一致性次之；與血管壁和纖維帽相比，脂質核回聲強度最低，其中硬斑脂質核亮度略高于混合斑，但二者都屬于回聲強度較高并含少許透聲分布，而軟斑的脂質核回聲強度較低，含有大量透聲區，亮度也相對較暗。這些回聲強度的特征與臨床超聲圖像表現一致。
圖6
B超仿真圖像
(a)混合斑(升序)；(b)混合斑(彌散)；(c)軟斑(峰谷)；(d)硬斑(降序)
Figure6. Simulated B-mode ulreasound images(a) mixed plaque (ascending); (b) mixed plaque (diffuse); (c) soft plaque (peak-valley); (d) hard plaque (descending)
對仿真的頸動脈斑塊B超圖像進行手工精細分割，分別對各類斑塊的血管壁、纖維帽以及脂質核進行回聲斑點直方圖統計及模型擬合、幾何尺寸測量以及誤差分析。回聲斑點分布的統計結果如表 4所列。與圖 3~圖 5的臨床統計結果比較，仿真回聲斑點的分布和參數值與臨床超聲的總體一致。其中，混合斑、硬斑及軟斑的血管壁和纖維帽都屬于萊斯分布(參數不同)，而混合斑、硬斑及軟斑的脂質核分別屬于瑞利、萊斯和K分布。斑塊中各組織的模型參數雖然存在一些偏差，但與臨床對應的結果大體一致，并處于給定參數的偏差范圍內。
表4
仿真的各類斑塊組織回聲斑點分布及參數
Table4.
Distributions and parameters of echoed speckles for simulated images of different kinds of plaque tissue
幾何尺寸測量及誤差統計結果如表 5所列。與預設值相比，混合斑、硬斑和軟斑的纖維帽相對誤差分別為5%、7.5%及10%，這是由于纖維帽厚度薄，而超聲側向分辨率低，超聲隨機斑點的影響較大，特別是對曲面形狀的纖維帽影響更顯著；而斑塊長度、管腔直徑以及血管壁厚度尺寸較大，超聲隨機斑點的影響較小，最大相對誤差僅為0.01%。混合斑、硬斑和軟斑的最大內膜厚度各不相同，其最大相對誤差分別為0.7%、0.5%及0.6%。
表5
仿真的各類斑塊組織幾何參數預設值及測量結果
Table5.
Preset values and measured results of geometrical parameters for simulated images of different kinds of plaque tissue
圖 7為仿真系統各模塊功能GUI界面。其中(a)為斑塊類型與幾何描述界面。該界面可根據需要選擇斑塊的類型與形狀，設定斑塊的幾何尺寸，生成相應的幾何模型。模塊中還可對生成的幾何模型進行放大、縮小、旋轉及尺寸測量等一系列操作，并可將模型保存下來。(b)為散射點模型描述界面。該界面根據表 3給出了散射點模型中不同的ρ和α取值及對應的回聲斑點分布及參數，用戶可在給定的取值范圍內自行選擇各組織散射點模型進行搭配，建立所需斑塊的散射點模型。(c)為超聲仿真成像界面。該界面基于Field II仿真軟件平臺，生成射頻回波信號，并經包絡提取、插值和對數壓縮生成B超圖像。
圖7
GUI系統的界面
(a)斑塊類型和形狀選擇；(b)散射點模型選擇；(c)生成的B超圖像
Figure7. GUI system interface(a) the selection of plaque types and shapes； (b) the selection of scatterer models； (c) the simulated B-mode ultrasound image
4 結論
本文應用FIELD II軟件對頸動脈斑塊進行超聲仿真，并利用MATLAB中的GUI結合斑塊的臨床特征建立了B超仿真系統。首先對斑塊的形狀尺寸以及特征進行了總結研究，再根據斑塊的回聲斑點特性以及臨床數據分析得到所需組織散射點模型，將模型代入FIELD II軟件平臺中得到超聲射頻回波信號生成超聲圖像，最后利用MATLAB中的GUI界面建立一個B超仿真系統，清晰地對斑塊B超成像過程進行模擬。對仿真圖像進行手工分割，分別對各組織回聲斑點和幾何尺寸分析測量，結果表明仿真圖像表現出的回聲特性能準確刻畫出頸動脈斑塊的組織特征，與臨床成像的結果一致；仿真圖像的幾何尺寸與預設值相比，纖維帽在5%~10%的誤差范圍內，其他組織在0.5%~0.7%的誤差范圍內，均符合超聲成像規律。建立的仿真系統能為頸動脈斑塊的特征提取和可視化研究提供有效的檢測數據源及仿真手段，為醫護人員培訓提供系統的仿真平臺。使用中先根據預設的斑塊類型、幾何尺寸及組織特征生成超聲圖像，受訓人員通過對圖像的分析獲取斑塊的特征，并與預設值相比較，提高醫護人員在臨床上對斑塊性質分析判斷的準確性。但本文方法仍有需要改進的地方，首先，斑塊組織模型只考慮了回聲分布特征，而沒有考慮組織的紋理特性；其次，沒有考慮脈動血壓的動力學因素對斑塊及血管壁的影響。針對以上問題進行研究是下一步工作的重點。