基于Rician分布散斑噪聲的超聲圖像模擬算法?

朋小秀 張 東

(武漢大學物理科學與技術學院 武漢 430072)

0 引言

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

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

1 方法

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

圖1 算法流程Fig.1 Algorithm process

1.1 采樣

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

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

圖3 采樣網格和算法參數Fig.3 Sampling grid and algorithm parameters

圖4 通過徑向極坐標采樣計算采樣點數據的偽代碼Fig.4 Pseudo-code for calculating sampling point data by radial polar sampling

1.2 重構

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

圖5 采樣點與校正圖像中像素點的對應關系Fig.5 The correspondence between the sampling points and the pixels in the corrected image

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

1.3 加噪

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

首先單色波方程為

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

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

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

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

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

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

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

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

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

其中，K稱作萊斯因子，當A →0 時，萊斯分布轉化為瑞利分布。

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

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

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

(2) 計算降采樣理想圖像的每個像素(x,y)的幅值，每個像素對應的鏡面分量累加到每個像素的散射特性上，來模擬每個像素點上復幅值的實部。累加散斑特性來模擬每個像素點上復幅值的虛部，則：

(3) 通過模擬出的每個像素點上復振幅的實部和虛部來進行最終的灰度幅值計算，來模擬不完全發育的斑點噪聲：

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

圖7 對每個采樣像素點進行加散斑噪聲的偽代碼Fig.7 Pseudocode for adding speckle noise to each sampled pixel

1.4 插值

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

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

2 實驗及結果分析

2.1 評價指標

為了保證模擬的超聲圖像更貼近真實的超聲圖像，需要對模擬的超聲圖像進行直方圖統計并進行擬合，為了驗證擬合得到的曲線是否符合超聲圖像的直方圖分布，來定性分析模擬的圖像是否從理論上符合真實的超聲圖像，需要對擬合結果進行一個評價，稱之為擬合優度檢驗。在這一部分，采用3個指標來評估所提算法的性能。

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

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

在本文的計算中，KS統計量的簡化計算公式為

當DobsDobs)計算得出：

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

其中，xi為理論分布對應區間的數據個數，需要注意的是xi需要大于5，小于5 的區間應該和前面的區間進行合并。卡方統計量的自由度為M?1?m，m為理論分布中需要估計的參數的個數，當M相對m來說比較大時，m可以忽略不計，則自由度近似為M?1。

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

卡方分布的累積分布函數為

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

若顯著性水平小于P值，不管χ2檢驗統計量大小如何，擬合模型都是可以被接受的。

2.2 實驗材料

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

2.3 實驗結果

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

2.3.1 實驗一：合成圖像

在作為體模的合成圖像中，包括了三角形、圓形、矩形、菱形和箭頭圖案，大小為700×600，整個算法過程如圖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 本文算法處理結果Fig.10 The results of the algorithm in this paper

2.3.2 實驗二：腎臟圖像

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

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

2.3.3 與真實的超聲圖像的對比

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

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

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

2.4 分布擬合結果

2.3節已經從視覺上證明了本文算法的優越性，為了對實驗結果進行更全面的評價，本節將從理論上定性地驗證本文算法的正確性。

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

表1 散射相量數量為100、δ 為1.5 時，KS 檢驗和卡方檢驗在樣本為2500、卡方自由度為47 時，不同顯著性水平下被接受的情況以及可決系數的數值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 檢驗和卡方檢驗在樣本為2500、卡方檢驗自由度為47 時，不同顯著性水平下對應的臨界值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 從模擬的圖像截取不同的區域用于Rician分布擬合Fig.14 Crop different regions from the simulated image for Rician distribution fitting

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

3 結論

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