李 毅,李 晶，陸永萍，黃紅兵

甲狀腺結(jié)節(jié)是一種高發(fā)病，早期發(fā)現(xiàn)并鑒別其良惡性對(duì)臨床治療及手術(shù)選擇有重要意義。甲狀腺結(jié)節(jié)診斷首選的方法是超聲影像檢查[1]，影像醫(yī)生依據(jù)自身經(jīng)驗(yàn)，通過(guò)考察結(jié)節(jié)大小、質(zhì)地、形狀、邊界和彈性等特征來(lái)鑒別其良惡性[2]。 這種特征選擇有強(qiáng)烈的主觀因素，缺乏量化標(biāo)準(zhǔn)，也不能更深入地分析細(xì)微的紋理特征。 因此，本研究使用數(shù)字圖像技術(shù)提取了77 個(gè)甲狀腺結(jié)節(jié)超聲圖像特征。 然而，由于甲狀腺結(jié)節(jié)超聲圖像的特征繁多，各特征間存在冗余性，因此，需要使用特征選擇方法，提取最優(yōu)特征子集，以降低特征維度，節(jié)省分類器的訓(xùn)練時(shí)間，提高分類識(shí)別性能，更好地理解關(guān)鍵特征[3]。本研究引入粒子群優(yōu)化[4]進(jìn)行特征選擇，獲得了理想的實(shí)驗(yàn)結(jié)果，并據(jù)此構(gòu)建了甲狀腺結(jié)節(jié)良惡性風(fēng)險(xiǎn)評(píng)估計(jì)算機(jī)輔助診斷系統(tǒng)。

1 甲狀腺結(jié)節(jié)超聲圖像特征提取

1.1 影像資料 2009 年6 月～2014 年12 月因甲狀腺結(jié)節(jié)性疾病來(lái)云南省第二人民醫(yī)院就診的患者共計(jì)31 654 例，均行超聲成像檢查。 檢查設(shè)備為法 國(guó) 聲 科 影 像 公 司 （SuperSonic Imagine） 的AixplorerR新聲威特種鑒別診斷超聲系統(tǒng)， 探頭為SL15-4，增益及二維條件固定。 在這些病例中，提交病理檢查且有電子存檔的患者共93 例， 其中惡性18 例，男14 例，年齡25～85 歲，平均54 歲，共1334張超聲結(jié)節(jié)影像。 根據(jù)結(jié)節(jié)視野的清晰程度，從中篩選得到245 張圖像，如圖1 示例。

1.2 形態(tài)特征 計(jì)算結(jié)節(jié)的周長(zhǎng)C ，面積A ，質(zhì)心O 到邊緣的距離Ri，i ∈[1,C]，Ri的均值R 、熵H(Ri) (F1)和標(biāo)準(zhǔn)差δ (F2)。 對(duì)結(jié)節(jié)進(jìn)行凸包運(yùn)算，得到凸包的周長(zhǎng)CH和面積AH。進(jìn)而計(jì)算結(jié)節(jié)的緊致度com(F3)、圓形度cir (F4)、平滑度smo (F5)、凹凸性cc(F6)、不規(guī)則度irr(F7)和輪廓溫度HT (F8)：

圖1 超聲結(jié)節(jié)影像示例

計(jì)算結(jié)節(jié)邊緣每個(gè)像素的曲率Ki，i ∈[1,C] ，得到曲率熵H(Ki) (F9)，由公式（7）計(jì)算將直線彎曲成特定曲線形狀所需要的能量，即彎曲能KE (F10)。

1.3 紋理特征 （1）基礎(chǔ)統(tǒng)計(jì)：灰度均值(F11)、標(biāo)準(zhǔn)差(F12)、平均能量(F13)和熵(F14)。 （2）灰度直方圖：二階矩(F15)、三階矩(F15)和四階矩(F16)。 （3）灰度梯度直方圖：對(duì)比度(F17)、熵(F18)和均值(F19)。（4）灰度共生矩陣[5]：均值(F20)、方差(F21)、逆差矩(F22)、對(duì)比度(F23)、熵(F24)、角二階矩（能量）(F25)和相關(guān)性(F26)。 （5）Laws 紋理能量測(cè)度[6]：通過(guò)特定的2D 算子核與原始圖像卷積來(lái)表達(dá)紋理標(biāo)記，即灰度級(jí)L 、邊緣E 、點(diǎn)S、波W 和紋R 五個(gè)算子分別完成一維信號(hào)的灰度平滑、邊緣檢測(cè)、點(diǎn)檢測(cè)、波狀檢測(cè)和紋狀檢測(cè)，從中選取LE 均值(F27)、 EL 均值(F28)、 SL 均值(F29)、 EE 均值(F30)、 LS 均值(F31)、 LE 方差(F32)、 EL 方差(F33)、 SL 方差(F34)、EE 方差(F35)、 LS 方差(F36)。 （6）灰度游程矩陣[7]用于描述紋理基元，定義為具有相同灰度值的最大共線連接的像素集合： 短游程逆矩(F37)、 長(zhǎng)游程矩(F38)、灰度不均勻性(F39)和長(zhǎng)度不均勻性(F40)。（7）分形維度方法基于分形幾何中的自相似性對(duì)紋理的尺度特征經(jīng)行度量，分形維越大，紋理越粗糙。使用盒子維法[8]計(jì)算分形維度量(F41)。（8）局部二值模式是描述圖像局部空間結(jié)構(gòu)的非參數(shù)算子，本研究使用多尺度一致編碼LBP 算子(F42)[9]。 （9）Gabor變換[10]作為一種短時(shí)傅里葉變換，是圖像時(shí)域-頻域分析的基礎(chǔ)，其中心頻率(F43)的差別可反映出紋理所具有的準(zhǔn)周期性特性。 （10） 局部傅里葉變換[11]， 選取極坐標(biāo)下8 個(gè)方向數(shù)量積的均值(F44-F52)和方差(F53-F60)，8 個(gè)方向相角的均值(F61-F68)和方差(F69-F77)。

2 基于粒子群優(yōu)化的特征選擇方法

2.1 特征選擇 特征選擇[3]本質(zhì)上是一個(gè)組合優(yōu)化問(wèn)題[12]，即尋求最優(yōu)化子集以提高分類器的識(shí)別能力和效率。 在搜索特征子集時(shí)，可以使用粒子群優(yōu)化提升搜索效率。

2.2 粒子群優(yōu)化 在鳥(niǎo)群聚和覓食行為的啟發(fā)下，Eberhart 和Kennedy 提出了最早的粒子群算法(particle swarm optimization, PSO )[4]。 假設(shè)在一個(gè)D維的目標(biāo)搜索空間中，有m 個(gè)粒子組成的群落，其中第i 個(gè)粒子表示為一個(gè)D 維的向量xi=(xi1,xi2,...,xiD)T，i =1,2,...,m， 即第i 個(gè)粒子在D 的搜索空間中的位置是xi。 將xi代入目標(biāo)函數(shù)計(jì)算其適應(yīng)度值，其歷史最優(yōu)值記為pi=(pi1,pi2,...,piD)T，群落的歷史最優(yōu)位置為pg=(pg1,pg2,...,pgD)T。 同時(shí)，賦予第i 個(gè)粒子一個(gè)速度，記為vi=(vi1,vi2,...,viD)T。

為了解決組合優(yōu)化問(wèn)題， 離散二進(jìn)制粒子群（binary barticle swarm optimization, BPSO）[13]采用下列公式，在每一代k 中對(duì)粒子的位置和速度進(jìn)行更新：

其中w 為慣性因子，表示粒子當(dāng)前速度對(duì)下一代速度的影響權(quán)重。 c1和c2為加速因子，表示粒子自身和群體運(yùn)動(dòng)的加速度，分別構(gòu)成"認(rèn)知"和"社會(huì)"部分， 是位于[0,1]之間的隨機(jī)數(shù)。

2.3 算法參數(shù) 特征總數(shù)即維度D，粒子可以表示為一個(gè)長(zhǎng)度為D 的二進(jìn)制串， 每一位即代表一個(gè)特征，為1 則該特征被選中，為0 則沒(méi)有被選中。

特征選擇需要考慮特征數(shù)目和分類精度兩個(gè)重要指標(biāo)，采用下列方法作為適應(yīng)度函數(shù)[4]：

其中， A(xi)表示以xi所表征的特征子集進(jìn)行分類的精度，本研究采用留一交叉檢驗(yàn)法[14]評(píng)價(jià)分類精度。 N(xi) 即xi所表征的特征子集的特征數(shù)目。

2.4 算法步驟 （1） 初始化： 設(shè)定最大迭代數(shù)Gmax、適應(yīng)度閾值F(xi)max等參數(shù)值，隨機(jī)生成粒子群， 并隨機(jī)設(shè)置每個(gè)粒子的位置xi0和速度vi0，并初始化各粒子的p0id和p0gd種群的；（2） 根據(jù)適應(yīng)度函數(shù)（公式10）計(jì)算粒子的適應(yīng)度值F(xik) ；（3）根據(jù)（公式8 和9）更新粒子的位置xik和速度vik；（4）更新各粒子的pKid和種群的pkgd；（5）判斷退出條件：達(dá)到最大迭代數(shù)Gmax，或找到超過(guò)適應(yīng)度閾值F(xi)max的特征子集。滿足則輸出最優(yōu)結(jié)果，不滿足則重復(fù)（3）（4）。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 特征選擇方法的比較 本研究對(duì)上述提取的77 個(gè)特征（F1～77）分別施用BPSO 和Relief[15]特征選擇方法，使用支持向量機(jī)[16]分類（表1）。 從表1可以看出，BPSO 獲得了數(shù)目最少的最優(yōu)特征子集和最高的分類準(zhǔn)確率。 在選出的特征集中，緊致度(F3)、平滑度(F5)、不規(guī)則度(F7)和灰度游程矩陣的長(zhǎng)度不均勻性(F40)均被三種算法選中。 也就說(shuō)，在甲狀腺結(jié)節(jié)超聲圖像中，結(jié)節(jié)形態(tài)的緊致度、邊緣的平滑度、形狀的不規(guī)則度于診斷具有重要價(jià)值。

表1 特征選擇方法的比較

3.2 分類器的比較 將BPSO 選出的最優(yōu)特征子集（F3、F5、F7、F14、F23、F24、F40）依次用于訓(xùn)練支持向量機(jī)（SVM）、貝葉斯分類器（Bayes）[17]和決策樹(shù)（DT）[18]分類器，均獲得了較好的分類效果。 其中，支持向量機(jī)仍獲得最高的分類準(zhǔn)確率，決策樹(shù)因其推理的簡(jiǎn)潔性在速度上具有優(yōu)勢(shì)。 見(jiàn)表2。

表2 分類器的比較

3.3 甲狀腺結(jié)節(jié)良惡性風(fēng)險(xiǎn)評(píng)估計(jì)算機(jī)輔助診斷系統(tǒng) 將特征提取、選擇，分類器鑒別集成為計(jì)算機(jī)輔助診斷系統(tǒng)（圖2）。 影像醫(yī)生采集甲狀腺結(jié)節(jié)影像傳入系統(tǒng)，在線標(biāo)記結(jié)節(jié)區(qū)域，系統(tǒng)提取上述77 個(gè)特征，使用分類器進(jìn)行良惡性鑒別，最終輸出結(jié)果。在數(shù)據(jù)積累到迭代條件（每增加100 例患者）時(shí)，將重新施用BPSO 進(jìn)行特征選擇，并生成新一代的最優(yōu)特征集，迭代分類器，以不斷循環(huán)，提升分類性能。

圖2 甲狀腺結(jié)節(jié)良惡性風(fēng)險(xiǎn)評(píng)估計(jì)算機(jī)輔助診斷系統(tǒng)流程圖

4 總結(jié)

臨床中，甲狀腺結(jié)節(jié)一般使用超聲檢查，結(jié)節(jié)圖像存在著豐富的形態(tài)和紋理等圖像特征。為了降低訓(xùn)練成本，提高分類性能，幫助理解關(guān)鍵特征，可以使用特征選擇方法對(duì)原始特征集進(jìn)行處理。本研究引入粒子群優(yōu)化進(jìn)行特征選擇，并施用于提取得到的77 個(gè)結(jié)節(jié)圖像特征， 用支持向量機(jī)訓(xùn)練得到的分類器精度可達(dá)到98.20%， 性能超過(guò)了常規(guī)特征選擇方法得到的同類分類器。 研究結(jié)果表明，結(jié)節(jié)形態(tài)的緊致度、平滑度和不規(guī)則度在甲狀腺結(jié)節(jié)良惡性鑒別中起到關(guān)鍵作用。

[1] 中華醫(yī)學(xué)會(huì)內(nèi)分泌學(xué)會(huì) 《中國(guó)甲狀腺疾病診治指南》 編寫(xiě)組. 中國(guó)甲狀腺疾病診治指南——甲狀腺結(jié)節(jié)[J]. 中華內(nèi)科雜 志, 2008, 47(10)：867-868. DOI： 10. 3321/j. issn： 0578-142 2008. 10.034.

[2] MR Carol, RW Stephanie, JW Charboneau, et al. Diagnostic ultrasound [M]. 3rd ed. Philadelphia： Elsevier Mosby,2005：735-757.

[3] Dash M, Liu H. Feature selection for classification[J]. Intelligent Data Analysis, 1997,1(4)： 131-156.

[4] Kennedy J, Eberhart R. Particle swarm optimization[C]// IEEE International Conference on Neural Networks [C]. Volume 4,1995(4)：1942-1948.

[5] Wouwer GV ,Scheunders P, Dyck DV. Statistical texture characterization from discrete wavelet representations[J]. IEEE Transactions on Image Processing, 1999, 8：592-598.

[6] Shapiro LG,Stockman GC.Computer Vision[M].[S.l.]：Prentice-Hall,Inc., Upper Saddle River, NJ, 2001.

[7] Galloway MM. Texture analysis using gray level run lengths[J].Computer Graphics And Image Processing, 1975, 4(2)：172-179.

[8] Pentland AP. Fractal-based description of natural scenes[J].Pattern Analysis and Machine Intelligence, IEEE Transactions on Medical Imaging, 1984(6)：661-674.

[9] Chu A, Sehgal CM, Greenleaf JF. Use of gray valu distribution of run lengths for texture analysis [J]. Pattern Recognition Letters, 1990, 11(6)：415-419.

[10] Chen CC, Chen CC. Filtering methods for texture discrimination [J]. Pattern Recognition Letters, 1999, 20 (8)：783-790.

[11] Zhou F, Feng JF, Shi Q. Texture feature based on local Fourier transform [C]//Image Processing, 2001. Proceedings. 2001 International Conference on. IEEE, 2001, 2：610-613.

[12] Schrijver A. Combinatorial optimization： polyhedra and efficiency[M]. [S.l.]：Springer, 2003.

[13] Kennedy J, Eberhart RC. A discrete binary version of the particle swarm algorithm[C]//Systems, Man, and Cybernetics,1997. Computational Cybernetics and Simulation, 1997 IEEE International Conference on. IEEE, 1997, 5：4104-4108.

[14] Kearns M, Ron D. Algorithmic stability and sanity-check bounds for leave-one-out cross-validation [J]. Neural Computation, 1999, 11(6)：1427-1453.

[15] Kira K, Rendell LA. The feature selection problem： traditional methods and a new algorithm[C]//AAAI, 1992： 129-134.

[16] Cortes C, Vapnik V. Support-vector networks [J]. Machine Learning, 1995, 20(3)：273-297.

[17] Zhang H. The optimality of naive Bayes[J]. AA, 2004, 1(2)： 3.

[18] Breiman L, Friedman J, Olshen R, et al. CART： classification and regression trees[M]. Wadsworth： Belmont, CA, 1983：156.