基于Nakagami成像紋理特征檢測微波熱消融凝固區(qū)的方法研究

王 赫，周著黃，吳水才

（北京工業(yè)大學(xué)環(huán)境與生命學(xué)部，北京 100124）

0 引言

肝細(xì)胞癌（hepatocellular carcinoma，HCC）是最常見的原發(fā)性肝臟惡性腫瘤，是癌癥死亡的主要原因之一[1]。微波熱消融術(shù)作為20世紀(jì)90年代的一種新興技術(shù)，目前已成為治療肝腫瘤的一種重要手段。微波熱消融是一種使用頻率≥900 MHz的設(shè)備誘導(dǎo)腫瘤死亡的電磁方法，其原理是通過水分子劇烈運(yùn)動產(chǎn)生的熱量造成腫瘤不可逆的凝固性壞死[2]。微波熱消融術(shù)的理想狀態(tài)是使壞死的凝固區(qū)范圍完全覆蓋腫瘤并形成5~10 mm的安全邊界，從而達(dá)到“適形消融”，這樣既可避免消融不足導(dǎo)致的腫瘤復(fù)發(fā)，又可防止過度消融對器官造成損害。因此，消融過程中的凝固區(qū)監(jiān)測是國內(nèi)外學(xué)者的研究重點(diǎn)。

為了彌補(bǔ)傳統(tǒng)B超成像在監(jiān)測過程中的不足，定量超聲通過分析背散射回波射頻信號來檢測消融凝固區(qū)，其中包括散射子濃度、間距及直徑等信息。目前，應(yīng)用在肝腫瘤消融領(lǐng)域的定量超聲主要有Nakagami統(tǒng)計(jì)分布模型法、零差K分布模型法、超聲衰減系數(shù)法、散射子平均間距法、回波信號去相關(guān)法等[3]。Nakagami統(tǒng)計(jì)分布模型法應(yīng)用比較廣泛，該方法通過定量映射構(gòu)建Nakagami參數(shù)分布圖描述散射介質(zhì)中超聲散射子排列和濃度。由于消融過程中消融區(qū)的氣泡作為強(qiáng)散射子導(dǎo)致該處散射子濃度增加，通過Nakagami參數(shù)可定量描述凝固區(qū)變化情況[4-6]。

圖像紋理特征是一種人類視覺認(rèn)知的關(guān)鍵組成部分，近年來常被應(yīng)用于計(jì)算機(jī)視覺系統(tǒng)中，灰度共生矩陣（gray-level co-occurrence matrix，GLCM）是其中的代表性方法之一[7]，常用于肝臟、乳腺腫瘤、腦腫瘤的醫(yī)學(xué)影像分析[7-10]，但目前還未見將GLCM應(yīng)用于肝腫瘤微波熱消融定量超聲圖像分析的報(bào)道。支持向量機(jī)（support vector machine，SVM）是由Vapnik[11]提出的數(shù)據(jù)分類方法，通過訓(xùn)練數(shù)據(jù)生成SVM模型并基于該模型實(shí)現(xiàn)分類。

本文提出一種基于GLCM和SVM的微波熱消融Nakagami參數(shù)圖檢測凝固區(qū)的新方法。通過采集Nakagami參數(shù)圖中凝固區(qū)與正常組織部分的圖像獲取SVM訓(xùn)練樣本，并基于訓(xùn)練樣本獲取SVM凝固區(qū)識別模型；對整幅Nakagami參數(shù)圖逐點(diǎn)滑動識別得到等大小的二值識別結(jié)果圖；基于識別結(jié)果圖使用多項(xiàng)式擬合（polynomial approximation，PAX）技術(shù)得到凝固區(qū)最終識別結(jié)果。

1 原理與方法

1.1 實(shí)驗(yàn)儀器與數(shù)據(jù)采集

1.1.1 實(shí)驗(yàn)儀器

本實(shí)驗(yàn)采用離體豬肝進(jìn)行微波熱消融實(shí)驗(yàn)。實(shí)驗(yàn)裝置平臺如圖1所示，主要實(shí)驗(yàn)儀器設(shè)備包括南京康友公司生產(chǎn)的KY-2000微波消融儀、美國Terason公司生產(chǎn)的T3000型超聲儀（用于采集超聲射頻回波信號）。

圖1 實(shí)驗(yàn)裝置平臺

1.1.2 數(shù)據(jù)采集

用80 W微波熱消融功率對離體豬肝進(jìn)行60 s的消融實(shí)驗(yàn)，通過超聲儀實(shí)時(shí)成像并以2幀/s的速度采集超聲背散射信號（數(shù)據(jù)采樣點(diǎn)為1 558，掃描線為256），共采集20組數(shù)據(jù)，獲取其超聲背散射信號，并沿消融針方向切開離體豬肝，測得其凝固區(qū)大小作為金標(biāo)準(zhǔn)。

1.2 Nakagami統(tǒng)計(jì)模型和成像

1.2.1 Nakagami統(tǒng)計(jì)模型

Nakagami統(tǒng)計(jì)模型是一種描述超聲背向散射包絡(luò)統(tǒng)計(jì)的通用模型[12]。在微波熱消融過程中，組織由于受熱產(chǎn)生氣泡并發(fā)生微小的結(jié)構(gòu)變化，在聲學(xué)模型中可將這種變化視為超聲散射介質(zhì)中的散射子濃度或排列發(fā)生變化。在Nakagami統(tǒng)計(jì)模型中，超聲背散射包絡(luò)振幅r的概率密度函數(shù)公式為

式中，Γ（·）表示伽馬函數(shù)；U（·）表示單位階躍函數(shù)；Ω表示尺度參數(shù)；m表示Nakagami參數(shù)。

假定E（·）表示統(tǒng)計(jì)均值，則可由公式（2）、（3）計(jì)算尺度參數(shù)Ω和Nakagami參數(shù)m：

式中，R為超聲背向散射包絡(luò)。

不同的m參數(shù)可定量表征由組織中散射體變化造成的不同包絡(luò)統(tǒng)計(jì)分布。當(dāng)01時(shí)，包絡(luò)統(tǒng)計(jì)分布服從后瑞利分布[4]。

1.2.2 Nakagami成像

Nakagami成像是基于Nakagami參數(shù)m映射的一種定量超聲圖像。通過滑動窗口采集包絡(luò)成像中的局部背散射包絡(luò)信號并估計(jì)局部m參數(shù)值，將其作為新像素賦值給當(dāng)前窗口；將上述滑動窗口以特定的窗口重疊率遍歷整幅包絡(luò)成像可獲得m參數(shù)矩陣?；瑒哟翱诘膶挾葘⑼瑫r(shí)影響Nakagami成像分辨力和m參數(shù)估算準(zhǔn)確度，為了更好地均衡二者關(guān)系，將窗口寬度設(shè)定為換能器脈沖長度的3倍[13]。圖2為離體豬肝微波熱消融B模式超聲成像、Nakagami成像和病理大體圖，其中圖2（a）為B模式超聲成像圖，將m參數(shù)矩陣加以顏色映射并成像可得到Nakagami成像，如圖2（b）所示。本研究選用的圖像紋理特征GLCM要求計(jì)算圖像為灰度圖像，因此將m參數(shù)矩陣進(jìn)行256灰度級數(shù)映射可得到如圖2（c）所示的Nakagami灰度成像。微波熱消融實(shí)驗(yàn)可得到消融針坐標(biāo)及凝固區(qū)長短軸大小，根據(jù)該長短軸計(jì)算橢圓偏心率并以消融針坐標(biāo)為中心勾畫凝固區(qū)范圍。圖2（d）為離體豬肝剖面大體病理圖（視為參考標(biāo)準(zhǔn)），以離體豬肝最大剖面的測量結(jié)果為金標(biāo)準(zhǔn)，檢驗(yàn)本文提出的算法的準(zhǔn)確度。

圖2 離體豬肝微波熱消融B模式超聲成像、Nakagami成像和病理大體圖

1.3 GLCM樣本采集及特征提取

1.3.1 GLCM樣本采集

Nakagami成像既包括凝固區(qū)范圍又包括正常組織部分，因此需要采集不同部分的樣本用于計(jì)算GLCM特征。本研究以像素大小為33×33的感興趣區(qū)域在每一幅Nakagami成像中分別采集4例凝固區(qū)圖像樣本和4例正常組織圖像樣本，凝固區(qū)樣本采集以靠近消融針位置處為宜。最終共在20例微波熱消融數(shù)據(jù)中采集樣本160例，其中凝固區(qū)圖像樣本80例，正常組織圖像樣本80例。圖3為樣本采集示意圖。

1.3.2 特征提取

GLCM是一種提取灰度圖像二階統(tǒng)計(jì)量的方法，其表現(xiàn)形式是一個(gè)以圖像灰度階數(shù)Ng為邊長的正方形矩陣。矩陣的元素p（i，j）為特定距離d和角度θ下的二階統(tǒng)計(jì)量概率值[14-15]。簡單地說，假定圖像分別有Nx行和Ny列，且灰度階數(shù)為Ng，那么圖像矩陣的行列可表示為Lx={1，2，…，Nx}及Ly={1，2，…，Ny}，圖像的灰度范圍可表示為Gx={0，1，…，Ng-1}。將灰度級數(shù)分配至每個(gè)Lx×Ly范圍內(nèi)的像素點(diǎn)，由相對頻率矩陣C（i，j）表示圖像的紋理特征關(guān)系。C（i，j）為在特定的（d，θ）窗口中灰度級別i和j出現(xiàn)的次數(shù)，概率密度可表示為

GLCM矩陣元素p（i，j）定義如下：

本實(shí)驗(yàn)基于GLCM矩陣計(jì)算得到4個(gè)特征值：能量（Energy）、對比度（Contrast）、同質(zhì)性（Homogeneity）、相關(guān)性（Correlation）。上述特征值的計(jì)算公式分別如（6）~（9）所示：

式中，μ、σ分別為均值和標(biāo)準(zhǔn)差。

能量是一種度量圖像灰度均勻性的特征，均勻的圖像通常僅由幾個(gè)灰度階數(shù)構(gòu)成，因此圖像越均勻其能量值越大；對比度反映了圖像中的局部灰度階數(shù)變化，若圖像中存在大量的灰度變化，則其具有高對比度值；同質(zhì)性是圖像中不同于能量的另一種灰度均勻性的度量特征值，圖像灰度越均勻，對應(yīng)的同質(zhì)性越大；相關(guān)性用于度量相鄰像素或特定位置像素之間灰度的線性相關(guān)性，它反映了灰度圖像的局部灰度依賴性，若灰度值沿某一方向延伸越長，則其相關(guān)性越大[14]。

GLCM的大小為灰度圖像的灰度階數(shù)Ng的二次方，過大的Ng會帶來大量的計(jì)算量進(jìn)而影響計(jì)算效率。為方便GLCM計(jì)算，將輸入圖像的Ng縮放至16位。設(shè)定步長d為1，方向（角度）θ為0°、45°、90°、135°。對上述4個(gè)特征值分別以設(shè)定的d和θ計(jì)算特征值，并以不同方向上的同一特征值的均值作為最終結(jié)果。

1.4 SVM模型訓(xùn)練

SVM是一種基于統(tǒng)計(jì)學(xué)習(xí)理論的機(jī)器學(xué)習(xí)方法，也是廣泛應(yīng)用于模式識別領(lǐng)域的常用工具。利用SVM可以通過不同的核函數(shù)組合多個(gè)自變量實(shí)現(xiàn)對因變量的估計(jì)。假定D為訓(xùn)練樣本集合，則有D={（X1，Y2），（X2，Y2），…，（Xn，Yn）}，其中第i個(gè)特征量Xi∈Rn，標(biāo)簽Yi∈{0，1}。SVM的實(shí)質(zhì)是建立一個(gè)具有最大間隔的劃分超平面，實(shí)現(xiàn)對訓(xùn)練樣本的劃分，也就是對公式（10）的求解：

式中，ω為決定超平面方向的法向量；b為位移項(xiàng)，決定超平面與原點(diǎn)的距離。由于并非所有樣本都可以在二維空間被線性可分，因此需要引用核函數(shù)將樣本映射到更高的維度，從而尋找超平面進(jìn)行劃分。常用的幾種核函數(shù)為線性核、多項(xiàng)式核、徑向基函數(shù)（radial basis function，RBF）核、拉普拉斯核、Sigmoid核。為了避免模型分類時(shí)出現(xiàn)過擬合或欠擬合，需要確定最優(yōu)懲罰參數(shù)c及核函數(shù)參數(shù)g。

本文選用LIBSVM工具箱進(jìn)行SVM模型訓(xùn)練[16]，確定核函數(shù)為RBF核。選用十折交叉驗(yàn)證法選擇最優(yōu)參數(shù)c和g，最終確定c和g分別為2.639和0.5。

1.5 凝固區(qū)識別

在獲取SVM模型后，對Nakagami參數(shù)圖以33×33大小的像素矩陣進(jìn)行逐點(diǎn)滑動識別。具體而言，首先對Nakagami參數(shù)圖以0進(jìn)行邊界填充，將填充后的圖像以33×33大小的像素塊進(jìn)行逐點(diǎn)塊分割處理。原始Nakagami參數(shù)圖大小為360×360，邊界填充并塊處理后共分割為129 600個(gè)子圖像塊。對每一個(gè)圖像塊分別進(jìn)行GLCM特征提取并進(jìn)行歸一化，由SVM進(jìn)行分類處理，凝固區(qū)標(biāo)簽為1，正常組織標(biāo)簽為0，最終重建可得到凝固區(qū)識別圖（二值圖像）。對凝固區(qū)識別圖進(jìn)行刪除小面積分量和孔洞填充2種形態(tài)學(xué)圖像處理，并結(jié)合PAX獲取最終凝固區(qū)識別結(jié)果[17-18]。圖4為本文算法的流程圖。

圖4 本文算法流程圖

2 實(shí)驗(yàn)結(jié)果

2.1 GLCM特征對比

箱線圖能夠直觀反映數(shù)據(jù)的分布關(guān)系，GLCM特征值箱線圖如圖5所示。皮爾森相關(guān)系數(shù)是一種描述2個(gè)變量之間線性關(guān)系強(qiáng)弱的統(tǒng)計(jì)量，2個(gè)變量的皮爾森相關(guān)系數(shù)的絕對值越接近0，則表明二者線性相關(guān)性越弱；越接近1，表明二者線性相關(guān)性越強(qiáng)[18]。分別計(jì)算上述4個(gè)GLCM特征值中凝固區(qū)樣本和正常組織樣本的皮爾森相關(guān)系數(shù)，結(jié)果顯示能量、對比度、同質(zhì)性、相關(guān)性4種GLCM參數(shù)的皮爾森相關(guān)系數(shù)的絕對值分別為0.104 7、0.091 8、0.048 6和0.109 3，4種GLCM特征值的樣本間皮爾森相關(guān)系數(shù)都接近0。結(jié)果表明，凝固區(qū)樣本與正常組織樣本的4種GLCM特征線性無關(guān)，可以作為訓(xùn)練樣本進(jìn)行SVM模型訓(xùn)練。

圖5 GLCM特征值箱線圖

2.2 SVM訓(xùn)練結(jié)果

使用分類準(zhǔn)確率（accuracy，ACC）、平均平方誤差（mean squared error，MSE）、相關(guān)系數(shù)r2作為衡量SVM識別能力的量化指標(biāo)。公式（11）~（13）列出了上述3種指標(biāo)的計(jì)算方法：

式中，TP為真正例；TN為真負(fù)例；FP為假正例；FN為假負(fù)例。將160例樣本數(shù)據(jù)按3∶1隨機(jī)劃分為120例用作訓(xùn)練集，40例用作驗(yàn)證集。實(shí)驗(yàn)結(jié)果顯示，訓(xùn)練后的SVM模型在驗(yàn)證集上的ACC為90%，MSE為0.1，r2為0.646 5。

2.3 凝固區(qū)識別結(jié)果

凝固區(qū)識別結(jié)果如圖6所示。圖6（a）為SVM凝固區(qū)識別模型初步識別結(jié)果，可以看到圖像下方有一些多余的識別結(jié)果。這些多余識別結(jié)果可能是圖像中存在一些高亮度區(qū)域噪聲，因而被SVM模型錯誤識別為凝固區(qū)。經(jīng)過圖像形態(tài)學(xué)處理（孔洞填充、小面積分量刪除）后得到圖6（b）。經(jīng)過實(shí)驗(yàn)測試得到，當(dāng)PAX的參數(shù)為4階、-6 dB等高線時(shí)候擬合效果最佳。圖6（c）為PAX成像，圖6（d）為PAX成像疊加在原始Nakagami成像后的結(jié)果圖。

圖6 凝固區(qū)識別結(jié)果圖

2.4 凝固區(qū)面積檢測準(zhǔn)確度

通常凝固區(qū)可近似看作橢圓形，因此利用公式（14）可計(jì)算凝固區(qū)面積S，根據(jù)公式（15）可得到PAX成像識別凝固區(qū)準(zhǔn)確度。

式中，SACC為凝固區(qū)面積檢測準(zhǔn)確度；SPAX為PAX成像計(jì)算面積；Sgold為豬肝消融后測量得到的面積（金標(biāo)準(zhǔn)）。本實(shí)驗(yàn)總計(jì)驗(yàn)證20例豬肝組織微波熱消融數(shù)據(jù)。傳統(tǒng)Nakagami方式在得到Nakagami參數(shù)成像后直接進(jìn)行PAX得到近似橢圓形凝固區(qū)，也可通過公式（14）進(jìn)行面積計(jì)算。結(jié)果顯示，經(jīng)過SVM模型識別后的凝固區(qū)與傳統(tǒng)Nakagami成像準(zhǔn)確度分別為（87.75±7.74）%和（84.38±13.52）%。

Bland-Altman是在分析化學(xué)及生物醫(yī)學(xué)中常見的相關(guān)性統(tǒng)計(jì)法[19-20]。20例數(shù)據(jù)經(jīng)SVM模型識別后的凝固區(qū)面積SSVM與金標(biāo)準(zhǔn)的Bland-Altman回歸分析圖如圖7所示。從圖中可以看出，20例數(shù)據(jù)中有18例（90%）處于95%一致性界限區(qū)間內(nèi)，95%的一致性界限區(qū)間的計(jì)算方法為Mean±1.96Std（Mean和Std分別為平均值和標(biāo)準(zhǔn)差）。結(jié)果表明，本實(shí)驗(yàn)提出的SVM模型識別的凝固區(qū)面積與金標(biāo)準(zhǔn)具有很高的一致性。

圖7 Bland-Altman回歸分析圖

3 討論

微波熱消融憑借其微創(chuàng)性、易恢復(fù)、可操作性強(qiáng)等優(yōu)勢在肝腫瘤治療領(lǐng)域的地位不斷提高，但消融過程中的凝固區(qū)檢測一直是亟須解決的問題。超聲作為常見的醫(yī)學(xué)成像技術(shù)被廣泛應(yīng)用于微波熱消融術(shù)中。常見的超聲監(jiān)測方法包括超聲彈性成像[21-22]和定量超聲成像[4，17，23-26]。

目前，Nakagami參數(shù)成像作為定量超聲監(jiān)測腫瘤微波熱消融熱凝固區(qū)變化的常用方法之一，已被證實(shí)可有效監(jiān)測凝固區(qū)面積。Zhang等[25-26]使用Nakagami參數(shù)成像監(jiān)測高強(qiáng)度聚焦超聲消融和微波熱消融，對Nakagami參數(shù)變化加以顏色映射，并將m參數(shù)發(fā)生顯著變化的區(qū)域視為凝固區(qū)。Zhou等[4]對Nakagami參數(shù)成像圖應(yīng)用PAX勾畫出等高線，從而確定凝固區(qū)范圍。本文通過提取Nakagami參數(shù)圖像的紋理特征對凝固區(qū)進(jìn)行識別，圖像紋理特征可表述肉眼無法觀察到的圖像局部高階信息從而對圖像進(jìn)行識別分割[27]。傳統(tǒng)的Nakagami參數(shù)成像僅利用m參數(shù)的大小變化單一地判斷凝固區(qū)范圍，本文提出的方法通過分析更深層的灰度階數(shù)變化關(guān)系來判斷凝固區(qū)范圍。這種灰度階數(shù)的變化通過計(jì)算特定距離、角度的像素對的概率統(tǒng)計(jì)對選定區(qū)域內(nèi)的紋理特征進(jìn)行定量描述，有效彌補(bǔ)了傳統(tǒng)方法簡單依靠m參數(shù)數(shù)值大小關(guān)系判定凝固區(qū)的不足。

本文選擇能量、同質(zhì)性、相關(guān)性、對比度共4個(gè)GLCM特征參數(shù)。能量和同質(zhì)性反映了圖像灰度階數(shù)的均勻程度。由于正常肝臟組織的肝細(xì)胞直徑為20~30μm，因此可被視為回波散射性低的弱散射子；而凝固區(qū)由于肝臟細(xì)胞發(fā)生微小結(jié)構(gòu)變化并且上方出現(xiàn)氣泡，因此可被視為回波散射性高的強(qiáng)散射子。上述關(guān)系在Nakagami-m參數(shù)上表現(xiàn)為m參數(shù)在正常組織處變化相較于凝固區(qū)處更微小，且這一特性在計(jì)算GLCM時(shí)的技術(shù)壓縮后更明顯。因此，正常組織樣本的能量及同質(zhì)性大于凝固區(qū)，在圖像上反映出的紋理特征更加均勻。對比度反映了圖像的明暗變化關(guān)系。由于熱消融過程中的溫度以消融針為中心向外逐漸減弱，因此對組織的損傷程度也逐步減弱，在Nakagami-m參數(shù)中表現(xiàn)為m值由凝固區(qū)邊界至消融中心逐步增大，因此凝固區(qū)內(nèi)的灰度階數(shù)變化相較于正常組織更大，在圖像紋理特征中表現(xiàn)出的明暗關(guān)系更復(fù)雜。同理，由于m參數(shù)由外向內(nèi)逐漸增大，在圖像紋理特征上表現(xiàn)為灰度值沿某一方向延伸，因此凝固區(qū)相關(guān)性高于正常組織部分。

將本文提出的方法與傳統(tǒng)Nakagami參數(shù)成像進(jìn)行對比可以看出，凝固區(qū)在傳統(tǒng)Nakagami參數(shù)成像上表現(xiàn)為高亮度區(qū)域但邊界并不明確，而基于圖像紋理特征的SVM識別模型可以識別出完整輪廓，結(jié)合PAX可更精準(zhǔn)地?cái)M合出凝固區(qū)邊界。通過對20例離體豬肝微波熱消融數(shù)據(jù)進(jìn)行驗(yàn)證后發(fā)現(xiàn)，本文所提的方法準(zhǔn)確度較原始方法提高約3.37%。經(jīng)過Bland-Altman回歸分析后發(fā)現(xiàn)，基于圖像紋理特征的識別結(jié)果與金標(biāo)準(zhǔn)高度一致。

本文的研究在未來應(yīng)在以下幾個(gè)方向進(jìn)行改進(jìn)。首先，通過改變不同的微波功率及消融時(shí)間，增加樣本數(shù)量并使樣本更加多元；其次，Nakagami成像本身存在的一些噪聲會影響SVM模型識別準(zhǔn)確度，因此需探究一種噪聲抑制方法改進(jìn)Nakagami成像；最后，可嘗試通過融合不同定量超聲圖像的GLCM特征值及定量參數(shù)訓(xùn)練SVM識別模型，提高凝固區(qū)的識別準(zhǔn)確度。

4 結(jié)論

本文提出一種基于GLCM和SVM識別Nakagami參數(shù)成像凝固區(qū)范圍的新方法。結(jié)果顯示，本方法檢測凝固區(qū)面積的準(zhǔn)確度達(dá)到（87.75±7.74）%，Bland-Altman回歸分析證實(shí)本方法估算面積與金標(biāo)準(zhǔn)具有高度一致性，可作為一種新方法對凝固區(qū)進(jìn)行識別檢測。