楊先鳳，張文俊，彭 博，李小蘭

(西南石油大學(xué)計(jì)算機(jī)科學(xué)學(xué)院，四川 成都 610500)

1 引言

在超聲彈性成像中，可以通過(guò)散斑追蹤技術(shù)來(lái)估計(jì)軟組織運(yùn)動(dòng)位移等信息。峰值跳變問(wèn)題就是超聲散斑運(yùn)動(dòng)追蹤技術(shù)在計(jì)算位移時(shí)受到去相關(guān)因素干擾造成的。當(dāng)回波信號(hào)去相關(guān)并且隨機(jī)相關(guān)峰的幅度超過(guò)與正確散斑運(yùn)動(dòng)相對(duì)應(yīng)的真實(shí)相關(guān)峰的幅度時(shí)，就會(huì)產(chǎn)生峰值跳變點(diǎn)。由于峰值跳變點(diǎn)的存在，超聲散斑追蹤得到的位移圖像不夠精確，對(duì)后期成像和醫(yī)生診斷有較大的影響。

針對(duì)峰值跳變點(diǎn)的研究，目前主要集中在其運(yùn)動(dòng)估計(jì)算法上，如采用先驗(yàn)信息的區(qū)域增長(zhǎng)算法[1]、貝葉斯算法以及采用位移正則化的GLUE算法[2]、提供子采樣位移估計(jì)精度的Coupled[3]算法等。這些算法在去除峰值跳變點(diǎn)時(shí)均只針對(duì)某個(gè)單因素(如壓縮比等)進(jìn)行了去除研究，但未曾研究多種因素影響下產(chǎn)生的任意斑塊半徑峰值跳變點(diǎn)的去除。

經(jīng)過(guò)對(duì)峰值跳變點(diǎn)分析研究，本文發(fā)現(xiàn)峰值跳變點(diǎn)與天文學(xué)星系中的點(diǎn)狀物分布類似，而天文學(xué)中對(duì)點(diǎn)狀物的分析常采用空間統(tǒng)計(jì)方法。因此，本文首次將空間統(tǒng)計(jì)方法-兩點(diǎn)相關(guān)函數(shù)運(yùn)用在峰值跳變點(diǎn)的研究中，找到如斑塊半徑等各因素影響下峰值跳變點(diǎn)規(guī)律；并分析比較多種峰值跳變點(diǎn)去除算法，結(jié)合得到的規(guī)律找到更優(yōu)的去除方案。通過(guò)仿真分析峰值跳變點(diǎn)的空間統(tǒng)計(jì)規(guī)律，提出更優(yōu)的峰值跳變點(diǎn)去除方案，為后續(xù)超聲彈性成像中峰值跳變點(diǎn)的去除研究提供先驗(yàn)信息和參考依據(jù)。

2 峰值跳變點(diǎn)兩點(diǎn)相關(guān)函數(shù)原理

相關(guān)函數(shù)計(jì)算可提供關(guān)于目標(biāo)區(qū)域特定位置不同部分間相互影響的信息，對(duì)目標(biāo)的空間統(tǒng)計(jì)特性進(jìn)行分析。在d維歐幾里得空間Rd中，峰值跳變圖的N點(diǎn)相關(guān)函數(shù)由Torquato[4]定義為：

=〈I(i)(x1)I(i)(x2)I(i)(x3)…I(i)(xn)〉

(1)

上式中，〈 〉代表整體平均值計(jì)算，I(i)(x)為指標(biāo)函數(shù)，定義如下

(2)

(3)

Berryman[5]提出了兩點(diǎn)相關(guān)函數(shù)(2-pcf)的離散化版本，它可以從M×N大小的峰值跳變點(diǎn)矩陣估計(jì)S2(x為水平方向，y為垂直方向)。其計(jì)算公式如下所示

(4)

根據(jù)以上公式，Berryman引入附加等式

(5)

若方程中的值k(固定半徑)和θ(角度)已知，則平均兩點(diǎn)相關(guān)函數(shù)可以寫成：

(6)

這里k是參數(shù)，小于或等于圖像最小尺寸的一半。

但是，采用以上方程計(jì)算高分辨率峰值跳變圖像時(shí)會(huì)出現(xiàn)運(yùn)算時(shí)間過(guò)長(zhǎng)或無(wú)解的情況。因此，Velasquez等人[6]提出一種簡(jiǎn)化的基于蒙特卡羅模擬的兩點(diǎn)相關(guān)函數(shù)計(jì)算，其計(jì)算過(guò)程如圖1所示。

圖1 兩點(diǎn)相關(guān)函數(shù)計(jì)算示意圖

上圖中兩點(diǎn)相關(guān)函數(shù)計(jì)算的主要步驟如下：

1) 隨機(jī)生成點(diǎn)p1，θ角。

2) 以該θ角和固定長(zhǎng)度k生成另一個(gè)隨機(jī)點(diǎn)p2。

3) 分別計(jì)算兩個(gè)隨機(jī)生成的點(diǎn)都在峰值跳變點(diǎn)中的命中數(shù)，記為Nk-hits。

4) 根據(jù)如下所示的公式計(jì)算兩點(diǎn)相關(guān)函數(shù)值為

(7)

其中，Ntotal是每個(gè)長(zhǎng)度k所計(jì)算的總數(shù)。根據(jù)Velasquez等人觀察到Ntotal>10000在計(jì)算中可提供最可靠的結(jié)果，因此在k的每個(gè)值處設(shè)置Ntotal=100000用于計(jì)算。

3 峰值跳變點(diǎn)優(yōu)化分析研究

本文使用了三類數(shù)據(jù)進(jìn)行兩點(diǎn)相關(guān)函數(shù)仿真分析和各去除算法對(duì)比分析，分別為計(jì)算機(jī)模擬數(shù)據(jù)、仿組織體模數(shù)據(jù)及真實(shí)乳腺超聲數(shù)據(jù)。通過(guò)快速歸一化互相關(guān)算法(NCC)[7-9]對(duì)壓縮前后的RF回波信號(hào)進(jìn)行運(yùn)動(dòng)位移估計(jì)，得到帶有峰值跳變點(diǎn)的位移估計(jì)圖。

3.1 峰值跳變分析數(shù)據(jù)獲取

對(duì)于計(jì)算機(jī)模擬數(shù)據(jù)，先在軟件中進(jìn)行組織體模的構(gòu)建，然后利用有限元求解器，獲取局部組織變形的真實(shí)位移信息，再通過(guò)其變形前后的RF信號(hào)來(lái)估計(jì)組織變形的位移信息。

這里列舉各影響因素下計(jì)算機(jī)模擬數(shù)據(jù)在快速歸一化互相關(guān)算法下產(chǎn)生的峰值跳變點(diǎn)圖，如圖2所示。

圖2 模擬數(shù)據(jù)軸向位移圖和峰值跳變點(diǎn)

圖2(a-f)分別展示了在互相關(guān)追蹤窗口尺寸為61*11、彈性模量為6000、簡(jiǎn)單模型和復(fù)雜模型、探頭中心頻率為7MHz、壓縮比為2%時(shí)的軸向位移估計(jì)圖，(a1-f1)表示對(duì)應(yīng)的峰值跳變點(diǎn)，黑色表示估計(jì)誤差小的正確像素點(diǎn)，白色表示估計(jì)誤差大的峰值跳變點(diǎn)。從圖中可以看出同一模型在不同因素下的峰值跳變點(diǎn)具有較大差距，故本文對(duì)各因素做兩點(diǎn)相關(guān)函數(shù)分析。

對(duì)于仿組織體模數(shù)據(jù)，采用凝膠仿組織體模進(jìn)行采集。體模大小為10cm*10cm，其內(nèi)部有一個(gè)半徑為10mm的球形包含物，球形目標(biāo)的彈性模量是背景組織的四倍。該模型的射頻回波信號(hào)使用Siemens SONOLINE AntaresTM采集，并使用40MHz的多行線性陣列換能器(VFX9-4)采樣，組織形變是手動(dòng)壓縮造成的。真實(shí)乳腺超聲數(shù)據(jù)是從人體乳腺病變超聲數(shù)據(jù)庫(kù)中任意選取的一組真實(shí)乳腺病變超聲回波數(shù)據(jù)。

快速歸一化互相關(guān)算法下仿組織體模數(shù)據(jù)和真實(shí)乳腺超聲數(shù)據(jù)產(chǎn)生的峰值跳變點(diǎn)如圖3-4所示，左邊是軸向位移圖，右邊是對(duì)應(yīng)的峰值跳變圖。可以看出仿組織體模中峰值跳變點(diǎn)較少，真實(shí)乳腺數(shù)據(jù)峰值跳變圖中存在大量呈斑塊狀的峰值跳變點(diǎn)。

圖3 仿組織體模軸向位移圖及峰值跳變點(diǎn)

圖4 真實(shí)乳腺數(shù)據(jù)軸向位移圖及峰值跳變點(diǎn)

3.2 兩點(diǎn)相關(guān)函數(shù)分析研究

在得到帶有峰值跳變點(diǎn)的組織位移估計(jì)應(yīng)變圖后，將圖中正確位移估計(jì)點(diǎn)和峰值跳變點(diǎn)看作兩種異質(zhì)材料，帶入兩點(diǎn)相關(guān)函數(shù)進(jìn)行計(jì)算。表1列舉了在圖2中隨機(jī)選取距離為1個(gè)像素(k=1)的兩個(gè)點(diǎn)，經(jīng)過(guò)兩點(diǎn)相關(guān)函數(shù)計(jì)算得到的值。

表1 k=1時(shí)兩點(diǎn)相關(guān)函數(shù)值

結(jié)合圖2和表1發(fā)現(xiàn)，圖中峰值跳變點(diǎn)越密集，對(duì)應(yīng)的兩點(diǎn)相關(guān)函數(shù)值越大。

3.2 峰值跳變點(diǎn)去除方案研究

為了找到在各因素影響下更優(yōu)的峰值跳變點(diǎn)去除算法，本文對(duì)目前已有的五種去除算法進(jìn)行了實(shí)現(xiàn)，包括聯(lián)合互相關(guān)和中值濾波(NCCM)算法、快速混合散斑跟蹤(Hybrid)算法[10]、聯(lián)合區(qū)域增長(zhǎng)貝葉斯(RGBMT)算法[1]、GLUE算法[2]以及Coupled算法[3]。對(duì)每種算法采用客觀評(píng)價(jià)指標(biāo)加權(quán)比噪聲比(簡(jiǎn)稱加權(quán)CNR)評(píng)價(jià)，其計(jì)算公式如下：

(8)

STG和SBK分別是目標(biāo)區(qū)域和背景區(qū)域的平均應(yīng)變，σTG和σBK為目標(biāo)方差和背景方差。CNR由目標(biāo)wTG和背景wBK的相對(duì)總面積的比值貢獻(xiàn)加權(quán)，二者總和為1。對(duì)于二維模型的信噪比，選擇整個(gè)二維包含物作為目標(biāo)，其余的跟蹤ROI(感興趣區(qū)域)為背景。使用加權(quán)CNR值可以描述各算法提高彈性應(yīng)變圖對(duì)比度的效果，值越大，去除峰值跳變點(diǎn)后生成的彈性應(yīng)變圖越精確。

計(jì)算出每個(gè)算法的加權(quán)CNR值后，結(jié)合前面小節(jié)中兩點(diǎn)相關(guān)函數(shù)的分析結(jié)果，找到針對(duì)某種情況下去除峰值跳變點(diǎn)的優(yōu)化方案。

有助于解決“城市熱島效應(yīng)”，維護(hù)城市生態(tài)平衡。透水混凝土能夠減少環(huán)境中的熱量，有助于解決“城市熱島效應(yīng)”。由于淺色和蜂窩狀的結(jié)構(gòu)，使其既不吸收熱量，也不儲(chǔ)存熱量，也就不會(huì)將輻射的熱量返回到大氣中，從而減少環(huán)境中加熱。雨水可通過(guò)透水混凝土迅速滲透下去，還原地下水，使大地恢復(fù)自然的儲(chǔ)水能力，防止地下水枯竭，改善植被的生存環(huán)境，恢復(fù)地表的水循環(huán)系統(tǒng)，營(yíng)造高質(zhì)量的自然生態(tài)環(huán)境，維護(hù)城市生態(tài)平衡。

4 峰值跳變點(diǎn)仿真研究

軟件實(shí)驗(yàn)環(huán)境包括MATLAB R2016b，ANSYS和Field II。為了方便觀察，本文將峰值跳變點(diǎn)通過(guò)一定閾值從位移圖中摘取出來(lái)。

4.1 計(jì)算機(jī)模擬數(shù)據(jù)兩點(diǎn)相關(guān)函數(shù)仿真

通過(guò)ANSYS軟件，對(duì)模型具體參數(shù)(互相關(guān)算法追蹤匹配窗口尺寸、彈性模量、模型復(fù)雜度、探頭中心頻率、壓縮比)進(jìn)行設(shè)置，并計(jì)算兩點(diǎn)相關(guān)函數(shù)值，結(jié)果如圖5所示。

圖5中三角曲線代表該因素影響下的最佳參數(shù)值，進(jìn)行綜合分析可得：五組因素下對(duì)應(yīng)k值增大，峰值跳變點(diǎn)總體呈下降趨勢(shì)；追蹤窗口最優(yōu)參數(shù)為101*11，大量的斑塊半徑(即k)為5，其余尺寸下為20；探頭中心頻率最佳參數(shù)為7MHz，大量斑塊半徑為9；越簡(jiǎn)單病變組織位移估計(jì)后峰值跳變點(diǎn)越少，大量斑塊半徑為16；組織材料因素下，大量斑塊半徑接近10且峰值跳變點(diǎn)較少。

圖5 計(jì)算機(jī)模擬數(shù)據(jù)兩點(diǎn)相關(guān)函數(shù)曲線

4.2 仿組織體模數(shù)據(jù)兩點(diǎn)相關(guān)函數(shù)仿真

仿組織體模數(shù)據(jù)兩點(diǎn)相關(guān)函數(shù)計(jì)算結(jié)果如圖6所示。在k=1處兩點(diǎn)相關(guān)函數(shù)值為0.0143，對(duì)應(yīng)圖3中產(chǎn)生峰值跳變點(diǎn)較少。曲線從k=80時(shí)相關(guān)值逐漸接近于0，說(shuō)明該仿組織體模數(shù)據(jù)產(chǎn)生的峰值跳變點(diǎn)斑塊雖然大，但不連續(xù)。

圖6 仿組織體模數(shù)據(jù)兩點(diǎn)相關(guān)函數(shù)曲線

4.3 真實(shí)乳腺數(shù)據(jù)兩點(diǎn)相關(guān)函數(shù)仿真

真實(shí)乳腺超聲數(shù)據(jù)的兩點(diǎn)相關(guān)函數(shù)值如圖7所示。

圖7 真實(shí)乳腺數(shù)據(jù)兩點(diǎn)相關(guān)函數(shù)曲線

從圖中結(jié)果來(lái)看，真實(shí)乳腺數(shù)據(jù)在k=1處兩點(diǎn)相關(guān)函數(shù)值為0.0926，對(duì)應(yīng)于圖4中峰值跳變點(diǎn)多且連續(xù)復(fù)雜。在k值逼近26個(gè)像素后兩點(diǎn)相關(guān)函數(shù)曲線斜率降低，該數(shù)據(jù)產(chǎn)生的峰值跳變點(diǎn)成塊狀后，大部分斑塊半徑趨近于26個(gè)像素。峰值跳變點(diǎn)斑塊半徑較大，圖片上的斑塊呈連續(xù)狀態(tài)，位移估計(jì)圖質(zhì)量較差。

4.4 峰值跳變點(diǎn)去除算法對(duì)比仿真

峰值跳變點(diǎn)各去除算法結(jié)果應(yīng)變圖的加權(quán)CNR值如表2-7所示，其中表2-5為計(jì)算機(jī)模擬數(shù)據(jù)。加權(quán)CNR值越大，則說(shuō)明去除峰值跳變點(diǎn)后生成的彈性應(yīng)變圖越準(zhǔn)確。

表2 組織材料各算法加權(quán)CNR值

表3 模型復(fù)雜度各算法加權(quán)CNR值

表4 探頭中心頻率各算法加權(quán)CNR值

表5 壓縮百分比各算法加權(quán)CNR值

表6 仿組織體模數(shù)據(jù)各算法加權(quán)CNR值

表7 真實(shí)乳腺數(shù)據(jù)各算法加權(quán)CNR值

經(jīng)過(guò)兩點(diǎn)相關(guān)函數(shù)仿真可得到某種因素下采用對(duì)應(yīng)參數(shù)設(shè)置將會(huì)產(chǎn)生更好的位移估計(jì)，結(jié)合表2-7得出針對(duì)峰值跳變點(diǎn)多因素影響下的去除方案：計(jì)算機(jī)模擬數(shù)據(jù)中，在組織材料因素影響下且斑塊半徑≤10的峰值跳變點(diǎn)使用NCCM算法進(jìn)行去除所獲得的彈性應(yīng)變圖更精確；在模型復(fù)雜度因素影響下且斑塊半徑≥10使用RGBMT算法所獲得的彈性應(yīng)變圖效果最優(yōu)；在探頭因素影響下使用RGBMT算法效果最優(yōu)；壓縮比因素影響下，當(dāng)壓縮比≤1%時(shí)使用NCCM算法效果最優(yōu)，當(dāng)壓縮比≥1%時(shí)使用RGBMT算法效果最優(yōu)。仿組織體模數(shù)據(jù)和真實(shí)乳腺數(shù)據(jù)中，處理半徑較大且不連續(xù)的峰值跳變點(diǎn)斑塊時(shí)，選擇RGBMT算法能夠獲得較好的位移估計(jì)結(jié)果和應(yīng)變圖像。

5 結(jié)論

本文主要對(duì)超聲散斑運(yùn)動(dòng)追蹤中峰值跳變點(diǎn)的分布規(guī)律及其去除算法進(jìn)行研究。首次將天文學(xué)中常用的兩點(diǎn)相關(guān)函數(shù)(2-pcf)運(yùn)用在峰值跳變點(diǎn)影響因素分析上，并對(duì)不同去除算法進(jìn)行仿真比較。最后將比較的結(jié)果和空間分布規(guī)律結(jié)合在一起，提出了一種具有針對(duì)性的峰值跳變點(diǎn)去除方案。這些研究結(jié)果可為今后的峰值跳變點(diǎn)研究及去除算法的改進(jìn)提供參考信息，具有一定的借鑒意義。