基于變分模態(tài)分解與多尺度排列熵的生物組織變性識別*

劉備 胡偉鵬 鄒孝 丁亞軍 錢盛友

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

(2018 年9 月26日收到; 2018 年10 月26日收到修改稿)

根據(jù)高強度聚焦超聲(HIFU)治療中超聲散射回波信號的特點, 本文利用變分模態(tài)分解(VMD)與多尺度排列熵(MPE)對生物組織變性識別進行了研究. 首先對生物組織中的超聲散射回波信號進行變分模態(tài)分解, 根據(jù)各階模態(tài)的功率譜信息熵值分離出噪聲分量和有用分量; 對分離出的有用信號進行重構(gòu)并提取其多尺度排列熵; 然后通過Gustafson-Kessel (GK)模糊聚類確定聚類中心, 采用歐氏貼近度與擇近原則對生物組織進行變性識別. 將所提方法應(yīng)用于HIFU治療中超聲散射回波信號實驗數(shù)據(jù), 用遺傳算法對多尺度排列熵的參數(shù)優(yōu)化后, 對293例未變性組織和變性組織的超聲散射回波信號數(shù)據(jù)進行了多尺度排列熵分析, 發(fā)現(xiàn)變性組織的超聲散射回波信號的多尺度排列熵值要高于未變性組織; 多尺度排列熵可以較好地識別生物組織是否變性. 相對于EMD-MPE-GK模糊聚類以及VMD-小波熵(WE)-GK模糊聚類變性識別方法, 本文所提方法中變性與未變性組織特征交疊區(qū)域數(shù)據(jù)點更少, 聚類效果和分類性能更好; 本實驗環(huán)境下生物組織變性識別結(jié)果表明, 該方法的識別率更高, 高達(dá)93.81%.

1 引 言

高強度聚焦超聲(high intensity focused ultrasound, HIFU)治療是一種新的無創(chuàng)腫瘤治療技術(shù). 它通過聚焦方式將聲能聚集于治療靶區(qū), 使靶區(qū)產(chǎn)生高溫, 從而使病變組織細(xì)胞內(nèi)的蛋白質(zhì)發(fā)生固化、變性和壞死, 同時又不損傷靶區(qū)之外的正常組織. 通過監(jiān)測生物組織的變性情況, 能掌握HIFU的治療效果, 對確保HIFU治療的安全高效有重要意義[1-3]. 迄今為止, 超聲領(lǐng)域的研究人員從多個方面對所取得的超聲信號進行研究, 期望提取出能準(zhǔn)確反映組織損傷變性的特征參數(shù), 如回波能量、聲衰減系數(shù)、頻率偏移、聲速和熵等[4-8]. 當(dāng)組織溫度升高到65 ℃以上時, 背散射信號能量會發(fā)生顯著變化[4-9]. Seip和Ebbini[10]采用回波信號的能量特征來檢測組織損傷變性, 實驗結(jié)果表明,信號能量識別準(zhǔn)確率達(dá)到82%. 在文獻(xiàn)[11]中, 生物組織的超聲衰減系數(shù)特征被用于估計組織損傷區(qū)域的溫度. 盛磊等[6]利用頻率偏移等特征參數(shù)檢測HIFU導(dǎo)致的組織凝固性壞死. 聲速可以非侵入式地估計組織溫度, 但是聲速受非線性以及熱膨脹的影響, 導(dǎo)致通過聲速變化測量體內(nèi)溫度的準(zhǔn)確性大大降低[7-12]. 明文等[8]采用超聲散射回波信號的小波熵(wavelet entropy, WE)特征對HUFU治療過程中的組織損傷變性進行評價研究, 同時對組織樣本進行損傷級歸類. 但到目前為止, 將超聲回波信號的多尺度排列熵(multi-scale permutation entropy, MPE)特征用于生物組織變性識別是否有效的研究還未見報道.

經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition,EMD)是將一個時間序列信號分解成一組有限的本征模態(tài)函數(shù)[13], 并且文獻(xiàn)[14]已經(jīng)將其應(yīng)用至從人體獲得的超聲回波信號, 然而, EMD存在端點效應(yīng)、模態(tài)混疊問題[15]. 為了解決上述問題,Dragomiretskiy和Zosso[15]提出一種新的變分模態(tài)分解(variational mode decomposition, VMD)方法, 能較好地抑制信號EMD過程中的端點效應(yīng)和模態(tài)混疊現(xiàn)象, 對噪聲具有更強的魯棒性.

從非線性角度分析生物組織中獲取的信號, 生物組織的損傷變性會導(dǎo)致信號的復(fù)雜度不同. 目前非線性分析算法可以通過其熵值有效地分析信號的復(fù)雜度, 例如尺度熵、近似熵、多尺度熵等[16-18].排列熵作為一種非線性分析算法, 具有計算簡單、抗噪能力強、魯棒性強的優(yōu)點, 因此被廣泛應(yīng)用于時間序列的復(fù)雜性分析中[19]. 多尺度排列熵將多尺度和排列熵結(jié)合起來, 可以更加有效地分析序列信息[20,21]. 本文從非線性角度分析生物組織中的超聲散射回波信號, 利用多尺度排列熵算法研究未變性生物組織和變性生物組織超聲散射回波信號的差異, 為HIFU治療診斷提供有效的依據(jù)和幫助. 然而, 若MPE算法中嵌入維數(shù)和尺度因子參數(shù)設(shè)置不合理, 將導(dǎo)致算法無法達(dá)到最佳的處理效果. 同時, 由于生物組織從未變性狀態(tài)到變性狀態(tài)的演變是一個漸變的過程[22], 因此提取的組織變性特征具有一定的模糊性, GK模糊聚類算法為解決這類問題提供了一條有效途徑[23].

本文結(jié)合VMD與MPE算法的優(yōu)點, 對生物組織變性的特征進行提取, 同時采用遺傳算法確定多尺度排列熵尺度因子參數(shù), 討論了本實驗環(huán)境下多尺度排列熵的最佳嵌入維數(shù); 并用GK模糊聚類算法實現(xiàn)了生物組織變性的識別分類.

2 原理與方法

2.1 VMD理論

VMD是一種基于維納濾波, 希爾伯特變換和外差解調(diào)的新型多分量信號分解算法. VMD可以將一個實際信號f(t) 分解成K個離散的模態(tài)μk(k=1,2,3,···,K), 與EMD算法不同,μk在頻域中的帶寬都具有特定的稀疏屬性, 是調(diào)幅調(diào)頻(AM-FM)信號, 可以有效地抑制EMD中出現(xiàn)的模態(tài)混疊現(xiàn)象. 每個μk緊湊圍繞相應(yīng)的中心頻率,并且其帶寬通過高斯平滑解調(diào)獲得. VMD處理過程中的約束變分問題為

VMD將信號分解為K個模態(tài)分量的具體步驟如下:

1) 初始化

2)μk和ωk分別由(3)式和(4)式迭代更新:

3) 根據(jù)(5)式更新λ:

4) 重復(fù)步驟 2)和 3), 直至滿足迭代終止條件

式中ε為判別精度,ε＞0 ;

5) 輸出結(jié)果, 得到K個模態(tài)分量.

2.2 功率譜信息熵

基于功率譜信息熵判斷各模態(tài)分量為有效超聲散射回波信號還是噪聲以及無效超聲散射回波信號的方法流程如下.

1) 計算各模態(tài)的功率譜, 并對幅值進行歸一化.

2) 根據(jù)歸一化后功率譜的最大值和最小值,得出功率譜的變化范圍. 對功率譜的變化范圍進行分段. 例如, 若信號歸一化后功率譜幅值的變化范圍是0—1, 分成10段, 那么幅值范圍為0—0.1為第一段, 0.1—0.2為第二段, 依次類推.

3) 統(tǒng)計每一段在功率譜幅值中出現(xiàn)的次數(shù),并求出該段在整個數(shù)據(jù)點中出現(xiàn)的概率.

4) 按照(7)式求解信息熵, 并根據(jù)熵值判斷該模態(tài)是否為噪聲分量或無效超聲散射信號:

H(x)為功率譜信息熵, 其中pi表示功率譜分段后第i段所對應(yīng)的概率,M表示總共的分段數(shù).

2.3 多尺度排列熵

多尺度排列熵是在排列熵基礎(chǔ)上的改進, 基本思想是將時間序列進行多尺度粗?；?然后計算其排列熵. 計算具體步驟如下.

1) 對序列長度為N的時間序列X={xi,i=1,2,···,N}進行粗?；幚? 得到粗?；蛄衴(s)j

得到粗?；蛄?/p>

式中s為尺度因子; [N/s]表示對N/s取整.

2) 對y(s)進行時間重構(gòu)得到

其中m為嵌入維數(shù),τ為延遲時間;l為第l個重構(gòu)分量,l=1,2,···,N-(m-1)τ.

4) 根據(jù)(10)式計算每個粗?；蛄械呐帕徐? 由此得到時間序列在多尺度下的排列熵.

當(dāng)Pr=1/m! 時 ,Hp(m) 達(dá) 到 最 大 值 l n(m!) ; 通 常將多尺度排列熵值Hp(m) 進行歸一化處理, 即

2.4 GK聚類與變性識別

GK模糊聚類算法是利用協(xié)方差矩陣能自適應(yīng)動態(tài)度量的模糊聚類算法, 假設(shè)輸入數(shù)據(jù)為GK模糊聚類的目標(biāo)函數(shù)為

式中θ≥1 為模糊指數(shù), 模糊指數(shù)會影響聚類效果,其值太大會導(dǎo)致各類之間相互重疊;Dij為第j個樣本與第i類聚類中心的馬氏距離,Di2j是一個平方內(nèi)積范數(shù).

利用拉格朗日乘法對目標(biāo)函數(shù)進行優(yōu)化, 使(14)式取得極小值, 其必要條件為

GK模糊聚類算法的具體步驟如下:

1) 確定聚類中心數(shù)目c以及模糊指數(shù)θ, 根據(jù)(16)式對隸屬度矩陣U進行初始化, 通過(18)式計算聚類中心νi;

2) 計算協(xié)方差矩陣Fi

由協(xié)方差矩陣Fi求出正定對稱矩陣Zi, 之后根據(jù)(15)式計算出平方內(nèi)積范數(shù)利用得到的平方內(nèi)積范數(shù)根據(jù)(17)式更新隸屬度矩陣U,若不滿足則增加迭代次數(shù),若滿足判斷式, 則終止迭代運算.

GK模糊聚類的聚類效果可用劃分系數(shù)和Xie-Beni (XB)指數(shù)進行檢驗.

劃分系數(shù)為

Xie-Beni指數(shù)為

式中δ為類的平均方差,dmin是類間最短模糊距離;對于聚類效果評價, 其中PC越接近1, 代表劃分越清晰, 反之劃分越模糊; XB越小, 代表類間分離的聚類就越好.

本文采用歐氏貼進度和擇近原則實現(xiàn)組織變性模式識別. 貼近度越大, 代表兩個模糊子集的相近程度越大. 對于待識別樣本模糊子集A、標(biāo)準(zhǔn)聚類中心模糊子集B, 其計算如下:

3 模擬仿真及應(yīng)用實例

3.1 信號分解與去噪

為研究EMD與VMD算法對噪聲的魯棒性,對含噪仿真信號進行分析, 所采用的仿真信號表達(dá)式為cos(2π·288t)+η, 其中η為高斯白噪聲, 標(biāo)準(zhǔn)偏差取0.1. 仿真信號及各對應(yīng)成分如圖1所示.

圖1 仿真信號及其對應(yīng)非噪聲成分的時域圖 (a)加入噪聲的仿真信號時域圖; (b)對應(yīng)非噪聲成分的時域圖Fig.1. Time-domain diagram of the simulated signal and its non-noise components: (a) Time-domain diagram of the simulated signal with added noise; (b) time-domain diagram corresponding to non-noise components.

圖2為加入噪聲的仿真信號經(jīng)EMD和VMD處理后的結(jié)果. 由圖2(b)可以明顯看出, 隨著EMD分解層數(shù)增多, 分解結(jié)果出現(xiàn)嚴(yán)重失真, 產(chǎn)生虛假分量, 即模態(tài)混疊現(xiàn)象. 而根據(jù)圖2(a)的VMD結(jié)果, 發(fā)現(xiàn)分解得到的前三個單頻模態(tài)分量與原信號對應(yīng)各成分比較一致, 說明VMD可以較好地解決信號分解過程中的模態(tài)混疊現(xiàn)象, VMD對噪聲具有較強的魯棒性.

圖3(a)為實驗獲得的超聲回波信號. 實驗中用HIFU輻照新鮮離體豬肉組織來改變其特性, 通過B超監(jiān)控來獲取超聲回波信號, 并經(jīng)數(shù)字示波器(Model MDO3032; Tektronix)轉(zhuǎn)化為數(shù)字信號后進行保存. B超探頭的中心頻率為3.5 MHz. 預(yù)先設(shè)定VMD算法中的分解模態(tài)的個數(shù)K= 5, 超聲回波信號經(jīng)VMD算法分解后, 不同K值下各模態(tài)分量如圖3(b)所示.

圖3 實際超聲回波信號與VMD結(jié)果 (a)實際超聲回波信號; (b)超聲回波信號VMD結(jié)果Fig.3. Actual ultrasonic echo signals and their VMD results: (a) Actual ultrasonic echo signals; (b) VMD results of ultrasonic echo signals.

對于每個模態(tài)分量, 其對應(yīng)的功率譜信息熵值越小, 就表示是噪聲分量或無效超聲散射回波分量的概率越大, 反之就認(rèn)為其中含有用超聲散射回波分量概率越大. 在本實驗環(huán)境下得到的超聲散射回波信號經(jīng)過VMD算法分解得到各階模態(tài)分量, 計算各模態(tài)的功率譜, 利用2.2節(jié)的方法得出各模態(tài)的功率譜信息熵. 實驗發(fā)現(xiàn)各模態(tài)分量的功率譜幅值分段數(shù)為16—21時, 噪聲分量、無效超聲散射回波模態(tài)分量與有用超聲散射回波模態(tài)分量區(qū)分效果最佳. 對于本實驗環(huán)境下的293例超聲回波數(shù)據(jù)經(jīng)過VMD得到的各階模態(tài)分量, 其中噪聲分量和無效超聲散射回波信號模態(tài)分量的功率譜信息熵在0.028—0.055之間; 如果模態(tài)分量中含有明顯的散射信號波形, 其功率譜信息熵在0.13—0.35之間. 上述兩類模態(tài)分量功率譜信息熵值的差異, 可以較好地區(qū)分噪聲分量、無效超聲散射回波模態(tài)分量和有用超聲散射回波模態(tài)分量.

3.2 多尺度排列熵的計算

通過篩選功率譜信息熵值在0.13—0.35的模態(tài)分量對信號進行重構(gòu), 對重構(gòu)后的信號進行參數(shù)優(yōu)化后的MPE分析, 利用遺傳算法來優(yōu)化尺度因子s. 參數(shù)設(shè)置為: 最大進化代數(shù)為100; 種群最大數(shù)量為20; 交叉概率為0.4; 變異概率為0.09; 超聲回波信號的截取點數(shù)為675, 一般來說, 嵌入維數(shù)m通常取3—8, 延遲時間τ=2 . 對重構(gòu)后的超聲散射回波信號進行尺度因子參數(shù)優(yōu)化, 取整后得到多尺度排列熵算法的尺度因子優(yōu)化參數(shù). 三組樣本的MPE參數(shù)優(yōu)化結(jié)果: 第一組和第二組樣本的MPE分析中, 尺度因子參數(shù)為12; 第三組樣本的MPE分析中, 尺度因子參數(shù)為13.

每組樣本包括總共15例超聲回波信號數(shù)據(jù),其中未變性狀態(tài)8例, 變性狀態(tài)7例. 圖4為嵌入維數(shù)m分別取3, 5和7時樣本1的MPE隨尺度因子的變化情況.

從圖4可以明顯看出, 嵌入維數(shù)m= 3和5時, 變性組織與未變性組織各尺度因子下的MPE值有較多交疊, 此時MPE不能較好地區(qū)分變性組織和未變性組織; 而嵌入維數(shù)m= 7時, 變性組織和未變性組織的MPE值有明顯區(qū)別, 變性組織的多尺度排列熵值要明顯高于未變性組織的熵值. 表1為嵌入維數(shù)為7時, 三組樣本中變性與未變性組織的各尺度排列熵平均值.

圖4 不同嵌入維數(shù)時變性與未變性情況的MPE分布 (a) m= 3; (b) m= 5; (c) m= 7Fig.4. MPE distribution of denatured and undenatured cases with different embedding dimension: (a) m= 3; (b) m= 5; (c) m= 7.

表1 三組樣本在嵌入維數(shù)m= 7時各尺度下排列熵平均值Table 1. The average entropy of the three samples at each scale when the embedding dimension m= 7.

由表1可見, 嵌入維數(shù)為7時, 尺度因子改變,同一樣本的MPE值發(fā)生變化; 尺度因子越大,MPE值越小, 且在不同尺度因子下, 變性組織的MPE值均高于未變性組織的MPE值. 另外, 尺度因子為12時, 第一組樣本與第二組樣本變性狀態(tài)與未變性狀態(tài)的MPE平均熵值差別較大, 與本文尺度因子參數(shù)優(yōu)化結(jié)果相符合; 三組樣本的MPE值隨尺度因子的變化趨勢類似. 后續(xù)研究中,延遲時間τ為2, 嵌入維數(shù)m取為7, 綜合考慮優(yōu)化結(jié)果選取尺度因子s為12. 通過以上分析可以發(fā)現(xiàn), MPE特征可以作為識別生物組織是否變性的一個重要特征參數(shù).

3.3 GK聚類與生物組織變性識別

圖5 不同聚類方法對未變性與變性生物組織超聲散射回波信號的聚類效果 (a) EMD-MPE-GK; (b) VMD-WE-GK;(c) VMD-MPE-GKFig.5. Clustering effect on ultrasonic scattering echo signals of undenatured and denatured biological tissues through different clustering methods: (a) EMD-MPE-GK; (b) VMD-WE-GK; (c) VMD-MPE-GK.

表2 變性與未變性組織識別結(jié)果Table 2. Recognition results of denatured and undenatured tissues.

取超聲散射回波信號的重構(gòu)后信號293例數(shù)據(jù), 構(gòu)成20組樣本. 隨機選取其中150例數(shù)據(jù)(共10組樣本)作為已知樣本, 后143例數(shù)據(jù)(10組樣本)作為待識別樣本, 選取各自的參數(shù)優(yōu)化后的MPE值作為特征參量. 提取的特征向量經(jīng)過GK模糊聚類算法處理后, 結(jié)合擇近原則計算歐氏貼近度, 從而實現(xiàn)生物組織變性的模式識別. 為進一步驗證所提方法的有效性, 分別采用EMD-MPE方法以及文獻(xiàn)[8]報道的WE分析(VMD-WE)方法提取特征參量, 最后采用GK模糊聚類對生物組織變性分類識別, 圖5分別為利用EMD-MPE,VMD-WE和VMD-MPE提取組織超聲散射回波信號多尺度排列熵的特征向量、經(jīng)GK模糊聚類后的二維空間分布以及二維等高線. 計算可得出各種方法的聚類PC, XB以及變性識別率, 識別結(jié)果見表2. 由圖5可以看出, 帶有組織變性與未變性特征的超聲散射回波信號經(jīng)本文方法處理后, 多尺度排列熵特征參量按照變性與未變性基本分布在2個聚類中心的周圍, 并且不同類型之間區(qū)分較為明顯, 相較于EMD-MPE-GK聚類和VMD-WEGK聚類方法, 本文方法變性與未變性特征二維空間交疊區(qū)域數(shù)據(jù)點更少, 而且隸屬度等高線區(qū)分變性組織與未變性組織的誤識別數(shù)據(jù)點更少, 類內(nèi)變性MPE特征更緊密, 分類效果更理想, 這說明基于VMD, MPE和GK模糊聚類的算法對生物組織是否變性具有較好的分類識別效果.

由表2可知, 對比VMD-MPE-GK聚類方法與EMD-MPE-GK聚類方法的識別結(jié)果, VMD分解重構(gòu)信號的XB更小, PC與變性識別率高于EMD分解重構(gòu)信號; 對比VMD-MPE-GK聚類方法與VMD-WE-GK聚類方法的識別結(jié)果, MPE的聚類效果優(yōu)于WE, 且變性識別率更高. 基于VMD, MPE和GK模糊聚類的生物組織變性識別方法能夠較準(zhǔn)確地識別生物組織是否變性, 變性識別率高達(dá)93.81%, 且本文所提方法相較于EMDMPE-GK聚類和VMD-WE-GK聚類方法, PC系數(shù)更接近1, XB指數(shù)更小, 從而證明了所提方法對生物組織變性識別分類的優(yōu)越性.

4 結(jié) 論

本文針對HIFU治療中超聲散射回波信號的特點, 利用VMD與MPE對生物組織變性識別進行了研究, 得出以下結(jié)論.

1) VMD能夠有效分解出包含生物組織變性信息的超聲散射回波信號模態(tài)分量; 利用參數(shù)優(yōu)化后的多尺度排列熵對重構(gòu)后的超聲散射回波信號進行特征提取, 結(jié)果表明, 選取延遲時間為2, 嵌入維數(shù)為7, 尺度因子為12時, 變性生物組織超聲散射回波信號的多尺度排列熵值要高于未變性生物組織的熵值, 多尺度排列熵可以較好地區(qū)分未變性組織和變性組織.

2) 利用多尺度排列熵對VMD去噪后超聲信號進行特征提取, 并結(jié)合GK模糊聚類和歐氏貼近度進行變性識別; 與EMD-MPE-GK聚類和VMDWE-GK聚類方法相比, 本文方法的聚類效果和變性識別率結(jié)果證明了所提方法的分類性能更好, 變性識別率更高, 可用于HIFU治療中生物組織變性識別.