朋小秀 張 東

(武漢大學(xué)物理科學(xué)與技術(shù)學(xué)院 武漢 430072)

0 引言

超聲成像技術(shù)由于其操作比較簡單、對人體沒有傷害等優(yōu)點(diǎn)，成為了醫(yī)學(xué)臨床上廣泛使用的診斷工具[1]，從而有大量的超聲圖像需要專家醫(yī)生進(jìn)行解讀。而超聲圖像中存在大量顆粒狀的斑點(diǎn)以及扁平的橢圓狀紋理，即所謂的“散斑噪聲”[2]。散斑噪聲的存在影響了專家對病人病理的判斷，降低了超聲圖像的可靠性[3]。因此，散斑噪聲的相關(guān)問題是醫(yī)學(xué)超聲圖像領(lǐng)域的一大研究重點(diǎn)，其一是對散斑統(tǒng)計(jì)的研究，最常用的模型之一是基于分辨率單元中包含大量散射體的假設(shè)，被稱為完全發(fā)育的斑點(diǎn)噪聲，根據(jù)中心極限定理，振幅被廣泛認(rèn)為滿足瑞利分布[4?5]，但是當(dāng)有孤立的強(qiáng)散射體存在時(shí)，將導(dǎo)致與瑞利統(tǒng)計(jì)的偏差，在這種情況下，后向散射振幅包絡(luò)可以用Rician分布來表征[6?7]。

其二是對散斑噪聲的抑制，由于真實(shí)的超聲圖像不存在無噪聲的原圖，為了評估各種去噪算法[8?11]的優(yōu)劣，在理想的圖像中添加可控噪聲是非常有必要的。關(guān)于超聲圖像的模擬，Perreault等[12]提出一種基于超聲圖像的采集過程，并且加入由Goodman[13]提出、Burckhardt[14]推廣到超聲領(lǐng)域的散斑形成模型，是超聲模擬算法中的一大進(jìn)步，但是這個(gè)算法只能模擬完全發(fā)育的斑點(diǎn)噪聲，即只能產(chǎn)生Rayleigh分布的斑點(diǎn)噪聲。而在真實(shí)的超聲圖像中，分辨率單元中可能存在孤立的強(qiáng)散射體，本文基于這種情況，提出了一種基于Rician 的不完全發(fā)育斑點(diǎn)噪聲的超聲模擬算法，并以合成圖像和腎臟圖像為體模進(jìn)行了模擬實(shí)驗(yàn)，通過在視覺上的對比以及對超聲圖像包絡(luò)進(jìn)行了直方圖統(tǒng)計(jì)并進(jìn)行了分布擬合檢驗(yàn)，在理論和視覺上都證明了該算法的優(yōu)越性。

1 方法

本文算法的主要步驟如圖1 所示，主要部分一共分為4 個(gè)階段，分別為采樣過程、重構(gòu)過程、加噪過程以及插值過程。

圖1 算法流程Fig.1 Algorithm process

1.1 采樣

本文算法的第一步是通過對像素網(wǎng)格進(jìn)行采樣來模擬超聲波束對平面的扇形掃描。圖2 為n條超聲波束的掃描場景。圖3 顯示了采樣網(wǎng)格所需要的參數(shù)，其中Φ表示扇形掃描的角度，w為原圖像Io的寬度，h為原圖像Io的高度，n為超聲波束的條數(shù)，m為每條超聲波束上采樣的像素點(diǎn)數(shù)，h0是換能器距離圖像上邊界的距離，amin和amax分別為軸向掃描的最近和最遠(yuǎn)距離。這兩個(gè)距離可以通過掃描超聲波束與原圖交點(diǎn)得到，但是一般可以根據(jù)實(shí)際超聲圖像的大小以及美觀靈活地調(diào)節(jié)扇形區(qū)域的大小，從而可以得到更加貼近真實(shí)的模擬超聲圖像。用于計(jì)算采樣點(diǎn)數(shù)據(jù)的偽代碼如圖4 所示，圖4描述了對原始圖像Io通過下采樣的過程得到采樣圖Is。

圖2 n 條超聲波束掃描場景Fig.2 n ultrasonic beam scanning scenario

圖3 采樣網(wǎng)格和算法參數(shù)Fig.3 Sampling grid and algorithm parameters

圖4 通過徑向極坐標(biāo)采樣計(jì)算采樣點(diǎn)數(shù)據(jù)的偽代碼Fig.4 Pseudo-code for calculating sampling point data by radial polar sampling

1.2 重構(gòu)

本文算法的第二步是對采樣圖像進(jìn)行校正。這里將由徑向極坐標(biāo)采樣得到的點(diǎn)轉(zhuǎn)換到的網(wǎng)格中，如圖5 所示，其中P點(diǎn)映射到P′點(diǎn)，A區(qū)域映射到A′區(qū)域。用于計(jì)算校正圖像的偽代碼如圖6 所示，將得到校正圖像Ir(i,j)。從某種意義上，圖6 的過程和圖4的過程是等價(jià)的。

圖5 采樣點(diǎn)與校正圖像中像素點(diǎn)的對應(yīng)關(guān)系Fig.5 The correspondence between the sampling points and the pixels in the corrected image

圖6 通過徑向極坐標(biāo)來進(jìn)行采樣的采樣點(diǎn)來計(jì)算校正圖像像素點(diǎn)的偽代碼Fig.6 Pseudocode for calculating corrected image pixel points by sampling points in radial polar coordinates

1.3 加噪

本文算法的第三步是對采樣圖像進(jìn)行加噪，也是核心部分。當(dāng)超聲脈沖的波長小于散射體大小時(shí)，會發(fā)生反射，發(fā)射的回波相位是相干的，經(jīng)疊加后會產(chǎn)生反射分量也稱鏡面分量，而當(dāng)超聲脈沖的波長與散射體大小相當(dāng)時(shí)，會產(chǎn)生散射分量。對于Rayleigh 噪聲的模擬只考慮到散射分量，而當(dāng)分辨率單元中存在孤立的強(qiáng)散射體時(shí)，就會存在一定的鏡面分量[4,6]。

首先單色波方程為

式(1)中，B(x,y,z)表示幅度，t表示時(shí)間，w表示頻率。其中幅度是一個(gè)復(fù)數(shù)值，即有

其中，|B(x,y,z)|為幅度的模，θ(x,y,z)為相位。單色波的幅度可以表示為

其中，f是一個(gè)任意變換(線性、非線性等)。假設(shè)無噪聲圖像Io以及采樣并重構(gòu)后的圖像Ir的幅值(即灰度)為回波振幅的線性變換，按一個(gè)常數(shù)進(jìn)行縮放，不失一般性，則有

根據(jù)物理學(xué)中的知識，復(fù)振幅可以通過回波分量ψi(x,y,z)來進(jìn)行計(jì)算：

其中，M是相量ψi(x,y,z)的總數(shù)，而這些相量又可以分為兩部分，鏡面相量和散射相量，即：

其中，ψsi(x,y,z)是鏡面相量，而ψdi(x,y,z)是散射相量，N為散射相量的數(shù)量，Bs是鏡面相量疊加的結(jié)果，這里稱為B的鏡面分量，Bd是散射相量疊加的結(jié)果，這里稱為B的散射分量。在本文的算法中，鏡面分量來自于降采樣之后的圖像信息，而散射分量來自于散射體的回波在發(fā)生相互干擾之后產(chǎn)生的非相干回波。使用[0,2π]的均勻分布來表示非相干回波的相位獨(dú)立性，不失一般性，假設(shè)每個(gè)相量ψdi(x,y,z)的實(shí)部和虛部都服從均值為0、方差為σ2的高斯分布，則散射分量的復(fù)分布由以下圓高斯分布的概率密度函數(shù)給出：

其中，u、v分別表示散射分量的實(shí)部和虛部，則散射分量的模V=(u2+v2)1/2將滿足Rayleigh分布，概率密度函數(shù)如下：

當(dāng)引入鏡面分量之后，回波復(fù)振幅的模將滿足Rician分布[15]，概率密度函數(shù)如下：

其中，A是主信號幅度的峰值δ2是多徑信號分量的功率，I0()是修正的0階第一類貝塞爾函數(shù)。并且有

其中，K稱作萊斯因子，當(dāng)A →0 時(shí)，萊斯分布轉(zhuǎn)化為瑞利分布。

接下來散斑模擬算法主要分為以下幾個(gè)步驟：

(1) 假設(shè)散射體的實(shí)部和虛部都服從高斯分布，為了方便計(jì)算，這里散射相量的數(shù)量N服從均勻分布U(c,d)，則：

其中，均勻分布的參數(shù)c和d可以隨著像素點(diǎn)(x,y)的變化而變化。

(2) 計(jì)算降采樣理想圖像的每個(gè)像素(x,y)的幅值，每個(gè)像素對應(yīng)的鏡面分量累加到每個(gè)像素的散射特性上，來模擬每個(gè)像素點(diǎn)上復(fù)幅值的實(shí)部。累加散斑特性來模擬每個(gè)像素點(diǎn)上復(fù)幅值的虛部，則：

(3) 通過模擬出的每個(gè)像素點(diǎn)上復(fù)振幅的實(shí)部和虛部來進(jìn)行最終的灰度幅值計(jì)算，來模擬不完全發(fā)育的斑點(diǎn)噪聲：

用于上述計(jì)算的偽代碼如圖7所示。

圖7 對每個(gè)采樣像素點(diǎn)進(jìn)行加散斑噪聲的偽代碼Fig.7 Pseudocode for adding speckle noise to each sampled pixel

1.4 插值

本文算法的最后一步是對重構(gòu)圖像進(jìn)行插值操作，來獲取完整的扇形圖像。為了方便進(jìn)行插值操作，這里使用卷積算子，因此進(jìn)行了算法的第二步，對采樣網(wǎng)格進(jìn)行了矩形校正過程，創(chuàng)建了一個(gè)n×m的臨時(shí)矩形圖像Ir(i,j)。對于采樣扇區(qū)內(nèi)的每個(gè)像素，計(jì)算它的笛卡爾坐標(biāo)，并插值周圍的像素，插值會導(dǎo)致分辨率單元的形狀發(fā)生改變，這與掃描的類型和掃描過程中所使用的插值技術(shù)有關(guān)，此處所采用的插值方法還是為雙線性插值。

插值過程的偽代碼如圖8所示。

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

2.1 評價(jià)指標(biāo)

為了保證模擬的超聲圖像更貼近真實(shí)的超聲圖像，需要對模擬的超聲圖像進(jìn)行直方圖統(tǒng)計(jì)并進(jìn)行擬合，為了驗(yàn)證擬合得到的曲線是否符合超聲圖像的直方圖分布，來定性分析模擬的圖像是否從理論上符合真實(shí)的超聲圖像，需要對擬合結(jié)果進(jìn)行一個(gè)評價(jià)，稱之為擬合優(yōu)度檢驗(yàn)。在這一部分，采用3個(gè)指標(biāo)來評估所提算法的性能。

(1) 判定系數(shù)(Coefficient of determination)[16]：判定系數(shù)在統(tǒng)計(jì)學(xué)中用于衡量因變量的變異中可以由自變量解釋的部分所占的比例，用此來判斷該統(tǒng)計(jì)模型的解釋力。假設(shè)樣本的數(shù)據(jù)集為y1,y2,···,yN，經(jīng)擬合模型計(jì)算的得到的理論值為，可決系數(shù)的定義如下：

(2) KS(Kolmogorov Smirnov)檢驗(yàn)[17?18]：KS檢驗(yàn)是一種非參數(shù)檢驗(yàn)，常用于判斷數(shù)據(jù)樣本與擬合模型給定的分布是否一致，假設(shè)N為樣本數(shù)據(jù)集的數(shù)量，GN(x)為樣本數(shù)據(jù)集的累計(jì)概率函數(shù)，G0(x)為待驗(yàn)證的累計(jì)分布函數(shù)，則KS檢驗(yàn)的統(tǒng)計(jì)量定義為

在本文的計(jì)算中，KS統(tǒng)計(jì)量的簡化計(jì)算公式為

當(dāng)DobsDobs)計(jì)算得出：

(3) 卡方檢驗(yàn)(Chi-Squared)[19]：卡方檢驗(yàn)是統(tǒng)計(jì)樣本數(shù)據(jù)的實(shí)際觀測值和理論期望值之間的偏離程度。假設(shè)樣本數(shù)據(jù)集中有N個(gè)數(shù)據(jù)，將其分到M個(gè)不相交的區(qū)間里面，每個(gè)區(qū)間對應(yīng)的數(shù)據(jù)個(gè)數(shù)為Xi(i=1,2,···,M)，χ2檢驗(yàn)的統(tǒng)計(jì)量為

其中，xi為理論分布對應(yīng)區(qū)間的數(shù)據(jù)個(gè)數(shù)，需要注意的是xi需要大于5，小于5 的區(qū)間應(yīng)該和前面的區(qū)間進(jìn)行合并。卡方統(tǒng)計(jì)量的自由度為M?1?m，m為理論分布中需要估計(jì)的參數(shù)的個(gè)數(shù)，當(dāng)M相對m來說比較大時(shí)，m可以忽略不計(jì)，則自由度近似為M?1。

然后將χ2檢驗(yàn)統(tǒng)計(jì)量與根據(jù)顯著性水平為α、自由度為M?1 時(shí)的χ2分布臨界表進(jìn)行比較。若，則擬合模型是可以被接受的，反之應(yīng)該被拒絕。

卡方分布的累積分布函數(shù)為

其中，n為自由度，另外假設(shè)分布的可信度可以通過FM?1(χ2)得到

若顯著性水平小于P值，不管χ2檢驗(yàn)統(tǒng)計(jì)量大小如何，擬合模型都是可以被接受的。

2.2 實(shí)驗(yàn)材料

本文用的無噪聲圖像即體模是含有不同幾何圖形的合成圖像，根據(jù)前人[8,20]的研究，這種設(shè)計(jì)是合理的。為了更加接近真實(shí)的超聲圖像，來自Field II Simulation Program 的腎臟圖像也作為體模加入到實(shí)驗(yàn)中。腎臟圖像可以在網(wǎng)站(http://field-ii.dk/)上找到。

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

本文采用Microsoft Visual Studio 2019進(jìn)行仿真實(shí)驗(yàn)，實(shí)驗(yàn)分為兩組進(jìn)行，第一組是人工合成的圖像，第二組是腎臟圖像，根據(jù)真實(shí)的超聲成像過程，探頭發(fā)射超聲波的數(shù)量一般為128，即這里扇形的掃描條數(shù)n取128，散射分量的數(shù)量N服從[10,1000]的均勻分布，其中σ分別取0.5、1.0、1.5、2.0。另外這里還把生成的偽超聲圖像與真實(shí)的超聲圖像進(jìn)行了對比。

2.3.1 實(shí)驗(yàn)一：合成圖像

在作為體模的合成圖像中，包括了三角形、圓形、矩形、菱形和箭頭圖案，大小為700×600，整個(gè)算法過程如圖9所示，生成的圖像如圖10所示。

圖9 本文算法以合成圖像為體模的流程圖(以δ=1.0 為例)Fig.9 The flow chart of using the synthetic image as the phantom in the algorithm in this paper (take δ=1.0 as an example)

圖10 本文算法處理結(jié)果Fig.10 The results of the algorithm in this paper

2.3.2 實(shí)驗(yàn)二：腎臟圖像

在這個(gè)部分的實(shí)驗(yàn)，使用腎臟圖像作為模型，大小為522×469，整個(gè)算法過程如圖11所示。生成的圖像如圖12所示。

圖12 本文算法處理結(jié)果Fig.12 The results of the algorithm in this paper

2.3.3 與真實(shí)的超聲圖像的對比

為了從視覺上證明本文算法的合理性，將上述實(shí)驗(yàn)的結(jié)果與真實(shí)的超聲圖像的細(xì)節(jié)進(jìn)行對比，如圖13 所示，其中圖13(a)、圖13(b)中偽超聲圖像的參數(shù)為N服從[10，1000]的均勻分布，其中σ為1.0。本文實(shí)驗(yàn)用到的真實(shí)超聲圖像來自于重慶醫(yī)科大學(xué)附屬醫(yī)院，是高強(qiáng)度聚焦超聲導(dǎo)航的子宮肌瘤圖像。

圖13 本文算法結(jié)果與真實(shí)超聲圖像的對比圖Fig.13 Comparison between the results of the algorithm in this paper and the real ultrasound images

從圖10 和圖12 以及圖13 的實(shí)驗(yàn)結(jié)果可知，該算法模擬的不完全發(fā)育的斑點(diǎn)噪聲在視覺上非常接近真實(shí)的超聲圖像，都有著弧狀的橢圓形紋理噪聲。并且本文算法還可以生成與原圖對應(yīng)的無噪聲圖像，如圖10(b)、圖12(b)所示，這為各種去噪算法的評價(jià)提供了便利條件。

2.4 分布擬合結(jié)果

2.3節(jié)已經(jīng)從視覺上證明了本文算法的優(yōu)越性，為了對實(shí)驗(yàn)結(jié)果進(jìn)行更全面的評價(jià)，本節(jié)將從理論上定性地驗(yàn)證本文算法的正確性。

首先對圖像噪聲分布進(jìn)行相似區(qū)域(即同質(zhì)區(qū))的矩形窗口采樣，再對采樣區(qū)域進(jìn)行直方圖統(tǒng)計(jì)，最后對此進(jìn)行擬合評價(jià)。這里根據(jù)概率論知識里面的各態(tài)歷經(jīng)性，取一定數(shù)量的散射分量的模擬超聲圖進(jìn)行分析(其實(shí)只要保證相同灰度的區(qū)域?qū)?yīng)的散射分量數(shù)量一致即可，這里為了簡便不做個(gè)別區(qū)分)。分別對不同的散射分量數(shù)量的模擬超聲圖的不同位置進(jìn)行分析，不失一般性，這里取散射分量數(shù)量為100、δ為0.15 的情況來進(jìn)行展開分析，這里選的采樣的窗口的大小為50× 50 的正方形區(qū)域，具體采樣區(qū)域如圖14 所示，并通過遺傳算法和列文伯格-馬夸爾特(簡稱LM)算法相結(jié)合的方法對直方圖進(jìn)行Rician 分布的擬合，擬合結(jié)果如圖15 所示。對于上述實(shí)驗(yàn)中的4 組樣本的Rician 分布擬合，所擬合的參數(shù)分別為σ=32.0984、A=124.8432；σ=28.9569、A=120.9332；σ=37.2979、A=125.6149；σ=32.8162、A=115.3947，根據(jù)可決系數(shù)、KS檢驗(yàn)、卡方檢驗(yàn)對擬合的評價(jià)方法，對擬合結(jié)果的評價(jià)具體如表1 所示(其中表2[17,21]為檢驗(yàn)評價(jià)常用的參數(shù)臨界值)，由表1 的結(jié)果可知，KS 和卡方兩種檢驗(yàn)都認(rèn)為所列的顯著性水平下，Rician 分布的假設(shè)是成立的，而可決系數(shù)也是都大于0.98，接近于1，即一致驗(yàn)證了在散射分量數(shù)量為100 的情況下，模擬噪聲圖像在不同的區(qū)域下都是符合Rician分布的。

表1 散射相量數(shù)量為100、δ 為1.5 時(shí)，KS 檢驗(yàn)和卡方檢驗(yàn)在樣本為2500、卡方自由度為47 時(shí)，不同顯著性水平下被接受的情況以及可決系數(shù)的數(shù)值Table 1 When the number of scattering components is 100 and δ is 1.5,the KS test and the chi-square test are accepted at different significance levels and the values of the coefficient of determination when the sample is 2500 and the chi-square degree of freedom is 47

表2 KS 檢驗(yàn)和卡方檢驗(yàn)在樣本為2500、卡方檢驗(yàn)自由度為47 時(shí)，不同顯著性水平下對應(yīng)的臨界值Table 2 Corresponding critical values of KS test and chi-square test at different significance levels when the sample is 2500 and the chi-square degree of freedom is 47

圖14 從模擬的圖像截取不同的區(qū)域用于Rician分布擬合Fig.14 Crop different regions from the simulated image for Rician distribution fitting

圖15 散射相量數(shù)量為100、δ 為1.5 時(shí)的不同區(qū)域的Rician 分布擬合情況Fig.15 Fitting of Rician distribution in different regions when the number of scattering phasors is 100 and δ is 1.5

3 結(jié)論

本文提出了基于Rician 分布的不完全發(fā)育的散斑噪聲超聲圖像模擬算法，考慮到了真實(shí)的分辨率單元中會存在孤立的強(qiáng)散射體的情況，并且結(jié)合了真實(shí)的超聲成像過程。另外本文算法可以生成一張沒有加噪的采樣插值圖像，為后續(xù)的去噪等應(yīng)用的評價(jià)工作創(chuàng)造了條件。同時(shí)利用可決系數(shù)的值以及KS 和卡方兩種檢驗(yàn)方法對模擬噪聲進(jìn)行噪聲分布擬合評價(jià)，證明本文算法模擬的偽超聲圖像從視覺和理論上都接近真實(shí)的超聲圖像。