胡偉鵬,鄒 孝,劉 備,趙新民,錢盛友

(湖南師范大學物理與電子科學學院,長沙 410081)

高強度聚焦超聲治療就是將超聲波的能量聚焦于靶區(qū),使靶區(qū)內(nèi)產(chǎn)生高溫,利用蛋白質(zhì)在一定的溫度下會產(chǎn)生變性使組織細胞壞死的原理,達到治療或組織切除目的[1]。在治療過程中,隨著超聲波聚焦于靶區(qū),靶區(qū)的細胞逐漸失去活性[2],位于超聲波傳播路徑上的正常組織也會吸收超聲能量,且越接近靶區(qū)的組織吸收的能量越多,因此使用超聲治療中需要實時監(jiān)測治療區(qū)域情況。迄今為止,超聲領域的研究人員已經(jīng)從超聲衰減系數(shù),聲速和熵等參數(shù)對采集的超聲信號進行研究,希望找到能準確反映生物組織特性的參數(shù)[3-6]。但采集的超聲背散射回波信號含有大量的噪聲,對判斷生物組織是否變性造成了很大的困難。減少超聲背散射信號中的噪聲,對準確判斷生物組織是否發(fā)生變性非常重要。

經(jīng)驗模態(tài)分解(Empirical Mode Decomposition,EMD)與變分模態(tài)分解(Variational Mode Decomposition,VMD)憑借其強大的分析和去噪能力廣泛應用于各個領域[7-11]。但是EMD得到的本征模函數(shù)(Intrinsic Mode Function,IMF)存在模態(tài)混疊,微弱信號難以提取,容易產(chǎn)生虛假信號[12-13],對于頻率相近的分量難以正確分離。與EMD相比VMD的魯棒性較好[14],在分析頻率相近的分量時更加容易分離,但是VMD分解的模態(tài)數(shù)是一個定值,在分解過程中部分噪聲信號會與有用信號分在同一個IMF中。本文基于VMD提出了多迭代變分模態(tài)分解(Multi-Iterative Variational Mode Decomposition,MIVMD),MIVMD在計算過程中,計算了各個IMF的能量比,將低能量比的IMF濾除,并根據(jù)重構信號中有用頻率的能量比進行多次分解重構。

由于超聲波在生物組織內(nèi)的非線性傳播特性,實驗采集的超聲背散射信號是一種非線性信號[15]。信號的復雜程度會根據(jù)生物組織變性程度的改變而改變。目前為止,我們通常通過計算信號的熵值來評估信號的復雜程度[16-18]。散布熵(Dispersion Entropy,DE)是一種新的評估時間序列動態(tài)特征的參數(shù),DE不會產(chǎn)生未定義的熵值,且時間序列的微小改變不會對DE造成很大影響[19]。復合多尺度散布熵(Composite Multi-scale Dispersion Entropy,CMDE)是DE的改進,進行多尺度處理,能從更多尺度下對時間序列的復雜性進行分析[20]。復合化能克服多尺度散布熵的缺點,能結合多個粗?；蛄械男畔?降低熵值的標準偏差使熵值更可靠。針對生物組織在一定溫度下逐漸產(chǎn)生變性的特性,GK模糊聚類可以根據(jù)輸入特征產(chǎn)生一個聚類中心,根據(jù)特征對聚類中心的隸屬度判斷輸入特征屬于哪一類,因此可以判斷生物組織是否變性。

通過上述分析,本文提出了MIVMD-CMDE-GK模糊聚類的生物組織變性識別方法。將采集的超聲背散射信號進行MIVMD重構,計算重構信號的CMDE,并將得到的所有CMDE作為特征參數(shù)通過GK模糊聚類得到兩個聚類中心,計算熵值與聚類中心的貼近度,根據(jù)貼近度的大小判斷生物組織是否變性。實驗結果表明,本文所提的MIVMD-CMDE-GK方法有效提高了超聲背散射信號對生物組織變性識別的能力,為超聲治療診斷提供了參考。

1 原理與方法

1.1 模態(tài)分解

1.1.1 變分模態(tài)分解

VMD算法中,原始信號分解得到的IMF為調(diào)頻-調(diào)幅信號,表達式為:

uk(t)=Ak(t)cos[φk(t)]

(1)

(2)

式中:uk(t)為VMD分解得到的第k個IMF,Ak(t)為瞬時幅值,φk(t)為瞬時相位,ω′k(t)為瞬時頻率。

VMD計算過程中,不斷更新IMF的中心頻率和帶寬,在各個IMF的帶寬和最小時停止計算,且各IMF之和為原始輸入信號。約束變分問題函數(shù)如下:

(3)

式中:{ωk}={ω1,ω2,…,ωk}為IMF的中心頻率。通過使用增廣拉格朗日函數(shù)對約束變分問題進行求解,如式(4)所示:

(4)

式中:α為二次懲罰參數(shù),λ(t)為拉格朗日懲罰算子。使用交替方向乘數(shù)法求解式(4)的鞍點,同時更新uk,ωk,λ三種參數(shù)。將原始信號x(t),通過VMD分解成k個IMF的方法如下[21]:

②更新uk,ωk,λ:

(5)

(6)

(7)

③設置精度ω,滿足當式(8)時,輸出結果uk。否則重復②步。

(8)

1.1.2 多迭代變分模態(tài)分解

由于噪聲頻率分布范圍廣,而VMD分解重構得到的模態(tài)數(shù)有限,因此得到每個IMF中都含有噪聲。為了更好的將信號中的噪聲濾除,提出了MIVMD,通過多次選擇有用信號能量占比高的分量重構,達到濾除噪聲的目的。MIVMD計算方法入下:

①對輸入信號進行VMD分解,IMF數(shù)為k。

②計算各IMF中的主頻能量占此IMF能量的比,能量比越大說明噪聲信號越少。

③將能量比小的IMF分量濾除,使用剩余分量重構信號。

④計算重構信號的頻譜能量,計算最大的前k-1個頻譜能量占重構信號總能量的比。

⑤設定占比閾值,當能量占比小于閾值時,將重構信號視為輸入信號并重復步驟(1)～步驟(4)。當能量占比大于閾值時,停止迭代,輸出重構信號。

1.2 復合多尺度散布熵

對長度為N采集信號x(t),計算信號的DE如下[22]:

①通過正太分布函數(shù),將采集信號x(t)映射成y(t)

(9)

式中:σ為標準差,mean為采集信號的平均值。

②將y(t)映通過式(10)映射成到1至c的集合z(t),c為類別數(shù)。

zc(t)=round[cy(t)+0.5]

(10)

(11)

④每種分散模式的概率為:

(12)

⑤計算DE:

(13)

CMDE是對DE的改進,先對原始輸入信號x(t)進行粗粒化尺度因子為τ的處理,得到第K個粗?；蛄?

(14)

然后對得到的τ個粗?；蛄猩⒉检厍缶?得到該粗?；叨纫蜃酉碌腃MDE為:

(15)

1.3 GK模糊聚類

GK模糊聚類是一種距離協(xié)方差矩陣能自適應動態(tài)度量的模糊聚類算法,假設輸入的特征為B={b1,b2,…,bN},GK模糊聚類的目標函數(shù)為[23]:

(16)

可以得到e個聚類中心Oi(i=1,2,…,e);f為模糊指數(shù),模糊指數(shù)會影響聚類效果,f太大會導致各類之間相互重疊;U為隸屬度矩陣U=[μij]e×N;μij表示第j個元素隸屬于第i類的隸屬度且滿足:

(17)

Dij表示第j個元素與聚類中心Oi的距離泛數(shù):

(18)

(19)

常用拉格朗日乘法對GK模糊聚類的目標函數(shù)進行優(yōu)化,得到最小值點(U,V),其必要條件為:

(20)

(21)

1.4 評價參數(shù)

為了便于比較不同方法的去噪效果,使用信噪比和均方根誤差來進行評價。

信噪比:

圖1 仿真信號與重構信號的頻譜對比

(22)

Ps,Pn分別代表有用信號和噪聲的有效功率。信噪比越大說明去噪效果越好。

均方根誤差:

(23)

signal′,signal,l分別為重構信號、未加噪信號、信號長度。均方根誤差越小說明去噪效果越好。

本文采用劃分系數(shù)PC、劃分熵系數(shù)PE和Xie-beni指數(shù)評估聚類效果[24]。PC,PE均與隸屬度矩陣有關,當PC越接近1與PE越大接近0時,聚類的效果越好。XB是用來評估類間距離的一種指數(shù),值越小,聚類類間距離越大,聚類效果越好。

劃分系數(shù)PC:

(24)

劃分熵系數(shù)PE:

(25)

Xie-beni指數(shù)表達式為:

(26)

式中:θ為類的平均方差,Jmin是類間最短模糊距離。

2 實驗分析

2.1 信號去噪

使用VMD、MIVMD、EMD,分別對加噪信號進行重構分解。未加噪仿真信號為x(t)=3sin(40πt)+2sin(200πt),加噪仿真信號為Q(t)=x(t)+η,η為高斯白噪聲。對加噪仿真信號分別進行VMD分解重構、MIVMD分解重構并與EMD分解重構進行比較。圖1(a)為信噪比為2.31 dB時加噪仿真信號頻譜圖,圖1(b)～1(d)分別為VMD,MIVMD,EMD三種方法重構信號的頻譜圖。從圖1可以發(fā)現(xiàn),EMD和VMD均能抑制中高頻噪聲,但是對兩個有用頻率之間的噪聲抑制效果較差。而本文提出的MIVMD不僅能抑制中高頻噪聲,還能有效地抑制兩個有用頻率之間的噪聲。

為了驗證MIVMD對低信噪比信號有更好的濾波能力,使用EMD、VMD和MIVMD對信噪比不同的仿真信號進行重構,并計算重構信號的信噪比和均方根誤差如圖2所示。由圖2可知,MIVMD重構信號的信噪比高于EMD與VMD重構信號的信噪比,MIVMD重構信號的均方根誤差低于EMD與VMD重構信號的均方根誤差。同時可以發(fā)現(xiàn),VMD與EMD重構信號的信噪比隨輸入信號信噪比增加而增加,而MIVMD在計算過程中,對各個IMF分量的能量進行了分析,因此MIVMD重構信號的信噪比在某一定值附近浮動。

圖3(a)、3(b)分別為通過光纖水聽器采集的正常與變性組織的超聲背散射信號,使用MIVMD對采集的超聲背散射信號進行分解重構,模態(tài)數(shù)設定為5,結果分別如圖3(c)、3(d)所示。可以觀察到通過MIVMD重構獲得的回波信號有明顯的脈沖波形,且波形震蕩衰減。

圖2 不同輸入信噪比下不同方法的去噪效果比較

圖3 超聲背散射回波信號

2.2 變性識別

為了驗證本文提出的生物組織變性識別方法的有效性,本文將該方法的判定結果與實際切片判斷結果進行比較。實驗使用HIFU對新鮮離體豬肉組織輻照來改變生物組織特性,使用熱敏電阻測量聲焦域處的生物組織溫度,最大溫度在90 ℃以內(nèi);考慮到HIFU環(huán)境中強超聲的影響,使用光纖水聽器(FOPH2000;RP acoustics)獲取監(jiān)控超聲的回波信號,并經(jīng)數(shù)字示波器(Model MDO3032;Tektronix)轉(zhuǎn)化為數(shù)字信號后進行保存。采集了16組樣本共288例實驗數(shù)據(jù)。為找到CMDE的最佳尺度,研究了3到20尺度下,第一組樣本的超聲背散射信號。為驗證MIVMD處理能使提取的特征參數(shù)更加有效,使用CMDE對經(jīng)VMD和MIVMD處理的重構信號進行分析。在使用CMDE進行分析時,為了避免將所有數(shù)據(jù)視為一種散布模式,需要將類別數(shù)設定大于1,本文將變性與未變性組織的信號視為兩種類別,因此設定類別數(shù)為2。計算CMDE時散布模式數(shù)應小于信號長度,而信號進行粗粒化處理將會減少信號長度,為了在分析不同尺度下CMDE時不受維數(shù)影響,將維數(shù)設定為2。設定延時是為了從不同時間分辨率下對信號進行分析,但延時過大時,信號的有用信息會被丟失且產(chǎn)生混疊,通過比較將CMDE的延時設定為2。

由圖4(a)可以發(fā)現(xiàn),使用VMD-CMDE方法得到的變性與未變性組織熵值之間存在差值,但是在尺度為6、16、19時變性與未變性組織熵值之間的差值小。通過圖4(b)可以發(fā)現(xiàn),由MIVMD-CMDE得到的變性與未變性組織熵值差值比VMD-CMDE方法得到的熵值差值大,在尺度8到11時,變性組織與未變性組織的熵值差較大,其中在尺度為9時差值最大。經(jīng)過上述分析,經(jīng)過MIVMD重構的實際超聲背散射信號的CMDE值能用于識別生物組織變性識別。

為了證明尺度為9時CMDE能較好的區(qū)分變性組織與未變性組織,分別計算了尺度為9時,經(jīng)MIVMD處理的16組實驗的變性組織與未變性組織的超聲回波信號的CMDE均值。同時與文獻[16-18]中的近似熵(Approximate entropy,ApEn)、樣本熵(Sample Entropy,SE)、模糊熵(FuzzyEntropy,FE)進行比較。如圖5所示,可以發(fā)現(xiàn),部分樣本組的SE與ApEn中出現(xiàn)變性時的熵值比未變性時的熵值小。所有樣本組的變性與未變性的FE不存在交疊,但是部分樣本組的變性與未變性的熵值差較小。所有樣本組的變性與未變性的CMDE值不存在交疊;相較于FE,變性與未變性的CEMD差值較大。變性組織的超聲回波信號的CMDE均值比未變性組織的超聲回波信號的CMDE均值高出0.110 8,約為8.12%。

圖4 不同信號處理方式時尺度對CMDE值的影響

圖5 變性與未變性組織的回波信號的不同熵

將所有特征參量視為未知參量進行GK聚類得到未變性與變性兩個聚類中心,計算CMDE與兩個聚類中心的貼近度,根據(jù)貼近度判斷樣本是否變性,并分別計算MIVMD-SE-GK、MIVMD-ApEn-GK、MIVMD-FE-GK和MIVMD-CMDE-GK四種方法的PC、PE、XB及變性識別率,結果如表1所示。

表1 特征參數(shù)對聚類效果的影響

通過對MIVMD-SE-GK、MIVMD-ApEn-GK、MIVMD-FE-GK和MIVMD-CMDE-GK方法分析可以發(fā)現(xiàn),MIVMD-ApEn-GK方法的PC指數(shù)最小且PE指數(shù)最大,聚類劃分在四種方法中最模糊,XB指數(shù)最大說明MIVMD-ApEn-GK方法的類間分離最小,聚類離散程度最大。與MIVMD-ApEn-GK方法相比,MIVMD-SE-GK和MIVMD-FE-GK方法的PC、PE、XB評價指數(shù)都有所改善。在四種方式中,MIVMD-CMDE-GK方法的PC指數(shù)最大且PE指數(shù)最小,聚類劃分在四種方法中最清晰,XB指數(shù)最小說明MIVMD-CMDE-GK方法的類間分離最大,聚類離散程度最小。在四種方法中,MIVMD-CMDE-GK方法的識別率最高,說明CMDE對變性生物組織的識別能力比ApEn、SE、FE要強。

3 結論

超聲治療中采集的超聲背散射信號含有大量的噪聲,且噪聲頻率范圍廣,對判斷生物組織是否變性造成了影響。針對上述問題提出了MIVMD,同時采用了CMDE算法對超聲背散射回波信號進行特征提取并結合GK模糊聚類進行生物組織變性識別。通過仿真和實例驗證可以得到一下結論:①MIVMD分解重構的去噪能力比EMD和VMD更好。②利用CMDE對處理過后的超聲背散射信號進行了特征提取,結果表明:當延時、類別、嵌入維數(shù)均設定為2時,變性組織超聲背散射信號的熵值要高于未變性組織超聲背散射信號的熵值。經(jīng)MIVMD處理的超聲背散射信號的CMDE值在尺度為9時兩者差別最大。③對MIVMD分解重構處理的超聲背散射信號進行SE、ApEn、FE和CMDE特征提取并進行GK聚類發(fā)現(xiàn),CMDE有更好地表征生物組織是否變性的能力。