徐 樸，郭 莉，馮衍秋，張鑫媛

南方醫(yī)科大學(xué)生物醫(yī)學(xué)工程學(xué)院//廣東省醫(yī)學(xué)圖像處理重點實驗室//廣東省醫(yī)學(xué)成像與診斷技術(shù)工程實驗室//粵港澳大灣區(qū)腦科學(xué)與類腦研究中心，廣東 廣州510515

擴散磁共振成像（dMRI）技術(shù)可通過檢測活體組織內(nèi)水分子的擴散運動狀態(tài)來反映生物組織的微觀結(jié)構(gòu)變化［1］。擴散張量成像（DTI）技術(shù)可通過采集不同擴散編碼方向的擴散加權(quán)（DW）圖像進行張量成像，用于描述水分子擴散運動的三維形式［2］。臨床上常用的DTI量化指標包括各向異性分數(shù)（FA）和平均擴散率（MD），這些參數(shù)可用于檢測急性腦中風(fēng)［3］、多發(fā)性硬化癥［4］、癲癇［5］以及腦腫瘤［6］等腦疾病的細微結(jié)構(gòu)改變。

與常規(guī)MR圖像相比，DW圖像信噪比較低，受噪聲影響嚴重。為提高DW圖像的質(zhì)量，已有多種后處理圖像去噪算法被相繼提出［7-15］。然而磁共振噪聲信號服從Rician［16,17］分布，常用的基于變換域的圖像去噪方法主要是針對加性高斯噪聲，如基于塊匹配的三維濾波（BM3D）［18］等，因此不能直接應(yīng)用于磁共振噪聲圖像上，目前常用的策略是先對圖像進行方差穩(wěn)定性變換（VST）［7,19］，得到噪聲空間分布不變的圖像，然后再利用基于變換域的去噪方法對該圖像去噪，去噪性能依賴于VST算法的效果。針對擴散加權(quán)圖像去噪，有研究提出了一種基于全局指導(dǎo)下的局部高階奇異值分解（GLHOSVD）去噪方法［9］，其去噪性能處于國際領(lǐng)先水平，且能夠降低條形偽影。然而，GL-HOSVD方法是對相似塊組成的高維數(shù)組進行高階奇異值分解（HOSVD）變換，因此在去噪同時不可避免的會引入條形偽影。此外，GL-HOSVD方法只適用于處理加性高斯噪聲圖像，因此也需在去噪之前對DW噪聲圖像進行VST變換。

本文中我們提出了一種聯(lián)合HOSVD稀疏約束與Rician 噪聲校正模型的去噪算法，該方法不需要借助VST變換，可直接處理DW噪聲圖像。此外，該算法不需要構(gòu)造相似塊組進行去噪，因此從根本上解決了條形偽影的問題。為驗證本方法的有效性，我們分別進行了仿真數(shù)據(jù)和真實數(shù)據(jù)實驗研究，將低秩+邊緣約束（LR+Edge）［10］，GL-HOSVD［9］，BM3D［18］和非局部均值（NLM）［20,21］4種去噪算法作為實驗對比方法，結(jié)果表明本方法無論是從定量還是定性方面均優(yōu)于4種實驗對比方法。

1 方法

1.1 高階奇異值分解去噪原理

高階奇異值分解是一種張量分解方法，可看成是對二維矩陣奇異值分解的一種推廣，可處理任意階數(shù)（N≥3）的張量［22］。已知三階張量AI1×I2×I3的HOSVD分解可表示為：

已知U(1)，U(2)，U(3)，核心張量S可由以下公式得到：

假設(shè)一幅無噪圖像具有結(jié)構(gòu)稀疏性，即經(jīng)過HOSVD分解后，可用較少的變換系數(shù)來表示。若A是受加性高斯噪聲污染的張量數(shù)據(jù)，那么HOSVD核心張量S中較小的系數(shù)主要來自于噪聲。因此可直接對S的系數(shù)進行閾值操作得到以達到去噪效果［23-25］。去噪后的張量數(shù)據(jù)可表示為

1.2 Rician噪聲校正模型

磁共振成像中，采集到的原始數(shù)據(jù)是K空間數(shù)據(jù)，經(jīng)傅里葉變換得到空間域的復(fù)數(shù)圖像，再經(jīng)模運算得到幅值圖像，其噪聲不再是簡單的加性高斯噪聲，噪聲信號服從Rician分布，其概率密度函數(shù)表示為［26］：

其中x為無噪信號，y為噪聲信號，σ為噪聲信號在復(fù)數(shù)域中高斯噪聲的標準差，I0為第1類0階修正貝塞爾函數(shù)。在高信噪比區(qū)域，py(y|x,σ)近似均值為x，標準差為σ的高斯噪聲，在低信噪比區(qū)域，py(y|x,σ)會產(chǎn)生偏移，噪聲信號y的均值也會隨之偏移理想信號x，且幅值圖像的噪聲波動大小與圖像信號有關(guān)，是空間變化的。因此基于變換域的去噪方法不能直接用于磁共振圖像上。目前現(xiàn)有方法是先將MR圖像進行VST變換，使得MR幅值圖像的噪聲分布是空間不變的，然后再對其進行去噪。而本文利用了Rician噪聲信號的期望，進行噪聲偏差校正。Rician噪聲信號y的期望可表示為［26］：

1.3 基于HOSVD稀疏約束和Rician噪聲校正模型的去噪算法

本文采用Rician噪聲信號期望進行噪聲偏差校正，將其融合到HOSVD 去噪模型當(dāng)中，從而不需要借助VST變換，可直接處理DW噪聲圖像。具體來說，本文所提去噪模型為：

其中Y為帶有Rician噪聲的DW圖像，數(shù)據(jù)大小為M×N×L，其中M×N為二維DW圖像大小，L為不同編碼方向的DW圖像的總數(shù)，X為待求解的無噪圖像，νNLM（Y）為采用向量形式的非局部均值對Y進行濾波得到的預(yù)處理圖像。

我們采用變量分離法求解以上去噪模型（公式（8）），即將去噪模型分解為以下兩個子問題，采用迭代的形式進行交替求解，最終得到最優(yōu)解，即最終的去噪圖像

對于子問題1，假設(shè)S,U(1),U(2),U(3)均已知，通過l-BFGS優(yōu)化算法得到X。對于子問題2，固定X，對X進行HOSVD變換和稀疏約束，得到U(1),U(2),U(3)。

其中公式（12）為硬閾值操作。本文采用經(jīng)典的自適應(yīng)系數(shù)收縮函數(shù)定義閾值

γ為噪聲估計水平的預(yù)設(shè)參數(shù)，

本文所提去噪算法的具體實施步驟如下：

1.根據(jù)圖像背景計算DW圖像的高斯噪聲σ，公式為Ybg為圖像背景的灰度值；

2.初始化X(1)，公式為

3.針對X(k)，采用滑動窗形式，以步長為Nstep，選取大小為t×t×L的DW三維塊Pi，對每個局部圖像塊Pi，進行HOSVD 分解（公式（11））以及硬閾值操作（公式（12）），得到

5.根據(jù)Z(k)計算vNLM的相似性權(quán)重，根據(jù)更新后的權(quán)重重新計算得到vNLM(Y)；

6.已知Z(k)或是S，U(1)，U(2)，U(3)，根據(jù)公式（9）得到X（k+1）；

7.若滿足迭代停止條件，則退出，X（k+1）為最終的去噪圖像，否則跳到步驟3，繼續(xù)進行迭代。

1.4 實驗數(shù)據(jù)

為了驗證所提方法的有效性，我們分別進行了仿真數(shù)據(jù)和真實數(shù)據(jù)實驗研究。

仿真數(shù)據(jù)實驗中的DW 圖像是由張量模型S0exp（-bgT Dg）產(chǎn)生的，S0為非擴散加權(quán)圖像，g為單位擴散編碼方向，D為擴散張量場，來自于生物醫(yī)學(xué)信息學(xué)研究網(wǎng)絡(luò)數(shù)據(jù)庫中的一個成年大鼠數(shù)據(jù)，該數(shù)據(jù)是由一個b=0圖像以及44個不同擴散編碼方向b=2000 s/mm2的DW圖像構(gòu)成，圖像大小為256×256。然后加入不同噪聲水平的Rician噪聲，公式如下：

其中X表示無噪的DW圖像，N1，N2為標準差為σ的高斯噪聲。經(jīng)過平方-相加-開方運算，幅值噪聲圖像Y服從Rician分布。為了驗證所提方法在不同噪聲水平下是否有效，我們將仿真數(shù)據(jù)最大值歸一化到1，然后仿真不同噪聲水平的DW圖像（σ=0.01，0.03，0.05，0.07和0.09）。

真實數(shù)據(jù)的采集設(shè)備是飛利浦3.0T TX MR掃描儀（Achieva,Philips Medical Systems,Best,The Netherlands），采用8通道腦部接收線圈，采集序列為4次激發(fā)的EPI序列。DW圖像通過multiplexed-SENSE重建得到，沒有進行加速采集。圖像成像參數(shù)為：6個不同擴散編碼方向b=1000 s/mm2的DW 圖像，一個b=0 圖像，TR/TE=3000/82 ms,層內(nèi)分辨率=1 mm×1 mm，層厚=4 mm，矩陣大小=220×220，層數(shù)=14，為了得到高信噪比的參考圖像，重復(fù)采集10次（NEX=10）。

1.5 定量指標

為了定量評價所提算法的有效性，本文采用了DW圖像的峰值信噪比（PSNR），結(jié)構(gòu)相似性測度（SSIM）和FA 圖的均方根誤差（FA-RMSE）分別對去噪后的DW圖像和后續(xù)的量化參數(shù)FA 圖進行定量分析，計算公式如下：

其中X為無噪DW圖像為去噪DW圖像，分別為圖像的均值，分別為圖像的標準差，c1，c2為常數(shù)，分別為無噪DW圖像和去噪后的DW 圖像得到的FA 參數(shù)圖。M1和M2為參與計算PSNR和各向異性分數(shù)均方根誤差（FA-RMSE）的所有像素值的個數(shù)。我們只采用前景區(qū)域的像素進行量化分析，從而排除背景干擾。PSNR值越大，SSIM值越接近于1說明去噪圖像的質(zhì)量越好，RMSE越小說明由去噪圖像得到的量化參數(shù)FA越準確。

1.6 參數(shù)設(shè)置

與BM3D［18］，NLM等［20,21］去噪算法類似，本算法的去噪性能依賴于參數(shù)設(shè)置，主要包括以下5個參數(shù)：圖像局部塊大小t×t，步長Nstep，噪聲標準差重估計尺度γ，基于Rician 噪聲校正模型保真項的權(quán)重λ，以及參與vNLM計算的平滑尺度β。通過優(yōu)化參數(shù)我們發(fā)現(xiàn)，當(dāng)t固定不變時，步長Nstep越小，去噪效果越好，但計算時間越長，通過權(quán)衡計算時間和去噪效果，我們固定，t=60，Nstep=10。通過仿真實驗，我們發(fā)現(xiàn)隨著噪聲水平變大，最優(yōu)值λ也會隨之變大。對于真實數(shù)據(jù)實驗，我們通過目測去噪圖像質(zhì)量來調(diào)節(jié)β和λ。實驗參數(shù)設(shè)置如表1所示。

表1 本文中仿真數(shù)據(jù)和真實數(shù)據(jù)實驗中的具體參數(shù)取值Tab.1 Parameter values in simulation data and real data experiment

2 結(jié)果

2.1 仿真數(shù)據(jù)實驗結(jié)果

圖1為所提方法與4種去噪方法（NLM，BM3D，LR+Edge，GL-HOSVD）在不同噪聲水平下的量化比較結(jié)果。從圖1可知，無論是從去噪圖像的PSNR，SSIM還是后續(xù)參數(shù)FA圖的RMSE來看，所提方法和GL-HOSVD算法的去噪性能均顯著優(yōu)于NLM，BM3D和LR+Edge。此外，在較低的噪聲水平下（1%和3%），所提方法和GL-HOSVD 的量化結(jié)果相似，進一步觀察，我們會發(fā)現(xiàn)所提方法的PSNR 和SSIM 略高于GL-HOSVD，而FA-RMSE 要略低于GL-HOSVD。在噪聲水平為7%和9%時，所提方法的PSNR 和SSIM 要明顯高于GLHOSVD，在噪聲水平為5%和7%時，所提方法的FARMSE要明顯低于GL-HOSVD。圖2顯示了在所提出的算法中不使用vNLM與使用vNLM預(yù)處理的結(jié)果對比。結(jié)果顯示，在所提去噪框架中引入vNLM濾波能夠進一步提高去噪圖像的PSNR以及降低后續(xù)量化參數(shù)FA圖的RMSE，具有更高的去噪性能。圖3給出了噪聲水平為5%時，不同去噪算法得到的去噪圖像、FA圖以及帶有彩色方向編碼信息的FA圖以及相應(yīng)的誤差圖。從圖3中可以看出，5種去噪方法都能有效降低圖像噪聲，恢復(fù)圖像主要的結(jié)構(gòu)信息，且5種方法都能夠顯著提高DW圖像的PSNR以及降低后續(xù)量化參數(shù)FA圖的誤差值。此外，本文方法和GL-HOSVD的去噪效果要明顯優(yōu)于其它3種去噪方法，去噪同時能夠更好的保留圖像的細節(jié)信息，DW圖像和FA圖的邊緣結(jié)構(gòu)信息更加清晰。雖然從整幅圖來看，本文方法和GL-HOSVD得到的去噪圖像以及FA圖的視覺效果相似，但從局部放大圖（圖4）來看，GL-HOSVD與BM3D得到的去噪圖像中含有一些條形偽影（如箭頭所示），并且這些條形偽影結(jié)構(gòu)也會出現(xiàn)在后續(xù)的量化參數(shù)圖中（如箭頭所示），而本文所提方法不但能夠得到高質(zhì)量的去噪圖像以及準確的FA圖，并且不存在這種偽影結(jié)構(gòu)。

圖1 所提方法和4種去噪方法的量化對比（NLM，BM3D，LR+Edge和GL-HOSVD）Fig.1 Comparison of the proposed method with 4 denoising methods(NLM,BM3D,LR+Edge and GL-HOSVD).

圖2 所提方法中不使用vNLM和使用vNLM預(yù)處理的結(jié)果對比Fig.2 Comparison of the results using the proposed method with and without vNLM preprocessing.

圖3 不同去噪方法的仿真實驗結(jié)果對比Fig.3 Comparison of different denoising algorithms on simulation data.

圖4 仿真實驗去噪結(jié)果的局部放大圖Fig.4 Enlarged views of the denoised images for simulation data.

2.2 真實數(shù)據(jù)實驗結(jié)果

圖5為真實數(shù)據(jù)的去噪方法比較結(jié)果。4種去噪算法都能有效降低圖像噪聲，恢復(fù)圖像結(jié)構(gòu)信息。但其中NLM的去噪結(jié)果相比于NEX=10過平滑，模糊了部分細節(jié)，在局部放大圖（圖6）中，NLM存在一些條形偽影。BM3D保留了較多的邊緣結(jié)構(gòu)信息，但同樣也在部分位置引入了條形偽影（圖6箭頭所示）。GL-HOSVD和本文所提方法的去噪圖像與NEX=10的參考圖像相似，甚至比NEX=10 的參考圖像噪聲還要小，F(xiàn)A 圖也是如此。但是從圖像的局部放大圖來看（圖6），GLHOSVD在去噪同時也引入了一些偽影（如箭頭所示），但是本文方法得到的去噪圖像以及FA圖均沒有偽影出現(xiàn)，該實驗結(jié)果與仿真實驗結(jié)果一致。

圖5 不同去噪方法對真實數(shù)據(jù)進行處理的結(jié)果比較Fig.5 Comparison of different denoising algorithms on real data.

圖6 真實數(shù)據(jù)去噪結(jié)果的局部放大圖Fig.6 Enlarged view of the denoised images for real data.arrow:stripe artifact.

3 討論

dMRI技術(shù)可通過采集不同編碼方向的擴散加權(quán)信號，來探測組織內(nèi)水分子擴散的運動狀態(tài)，從而反應(yīng)組織的微觀結(jié)構(gòu)變化，以及進行腦白質(zhì)神經(jīng)纖維束的構(gòu)建［29］。如各向異性分數(shù)和表觀擴散系數(shù)能夠探測多種腦疾病的病理組織改變。通過dMRI技術(shù)得到的大腦神經(jīng)纖維束能夠指導(dǎo)臨床中腦腫瘤手術(shù)方案的制訂，以及進行術(shù)后的療效評估［30］。但與常規(guī)T1和T2圖像相比，DW圖像受噪聲影響嚴重。雖然多次采集平均方法可提高DW圖像的信噪比，但是隨著采集次數(shù)的增加，采集時間也會增加，從而增加了采集過程中運動偽影的可能性。

目前很多研究學(xué)者提出采用后處理去噪技術(shù)提高圖像的信噪比。由于磁共振噪聲服從Rician分布，當(dāng)信噪比較高時，Rician噪聲信號近似高斯噪聲，期望近似等于理想信號，而當(dāng)信噪比較低時，Rician噪聲信號近似瑞利分布，期望有偏于理想信號值。NLM作為圖像域的去噪算法通常需要聯(lián)合Rician 噪聲校正模型，NLM自適應(yīng)差，去噪時容易出現(xiàn)過平滑的問題［8］。LR+Edge可以將Rician噪聲融合到最大后驗中，因此也可以解決Rician噪聲校正的問題，然而LR+Edge需要足夠魯棒的噪聲估計，較差的噪聲估計將導(dǎo)致比較差的去噪效果［10］。MR圖像在實部和虛部信號噪聲服從高斯分布，在相位校正過程中針對復(fù)數(shù)域的實部和虛部去噪，最終可獲得去噪后的幅值圖像，然而在臨床中得到的數(shù)據(jù)通常是重建后的幅值圖像而不是實部和虛部圖像，因此該方法在廣泛應(yīng)用于MR去噪方面有一定的局限性［11］。目前現(xiàn)有的基于變換域的圖像去噪方法不能直接處理磁共振噪聲圖像，尤其是對于這種信噪比較低的DW圖像。據(jù)本文所知，現(xiàn)有基于變換域的去噪方法均是先對MR 圖像進行噪聲方差處理，然后再進行圖像去噪。GL-HOSVD 也是一種基于變換域的去噪方法［9］，GLHOSVD借鑒了BM3D［18］的思想，將結(jié)構(gòu)相似的塊放在一起，構(gòu)成高維數(shù)據(jù)，然后進行HOSVD變換，因此去噪效果極佳，在擴散加權(quán)圖像去噪方面處于國際領(lǐng)先水平。但是GL-HOSVD也需要結(jié)合VST變換使用，去噪效果高度依賴VST算法的性能。此外，通過尋找相似塊構(gòu)造結(jié)構(gòu)更加稀疏的高維數(shù)組這類去噪方法，雖然能夠得到較好的去噪效果，但是在噪聲水平較大的時候，會引入偽影結(jié)構(gòu)，尤其是對這種信噪比較低的DW 圖像。雖然GL-HOSVD通過全局預(yù)濾波來指導(dǎo)后續(xù)的基于塊匹配的HOSVD去噪，可以有效減少去噪過程中引入的偽影，但是并沒有從理論上解決此類問題。

為了解決以上問題，我們提出一種新穎的HOSVD去噪方法，建立了基于HOSVD稀疏約束和Rician噪聲校正模型的聯(lián)合去噪方法，所提去噪模型可以直接處理Rician噪聲，不需要借助VST算法。此外，本文假設(shè)，在噪聲水平較大時，在圖像均勻區(qū)域或是灰度變化比較緩慢的地方，噪聲產(chǎn)生的灰度變化起主導(dǎo)作用，通過尋找相似塊來構(gòu)造高維數(shù)組，會將具有相似噪聲模式的圖像塊聚在一起，在去噪過程中錯把模式相似的噪聲當(dāng)作圖像結(jié)構(gòu)保留下來，形成條形偽影。此外，這些相似塊存在高度重疊區(qū)域，重疊區(qū)域的噪聲值相同，導(dǎo)致某些噪聲值出現(xiàn)的概率變大，因此構(gòu)成的高維數(shù)組中的噪聲并不是嚴格的隨機噪聲，這也是導(dǎo)致偽影問題的一個可能因素。因此，本文所提方法對整幅DW圖像或是每個局部DW圖像單獨進行HOSVD去噪，同一結(jié)構(gòu)塊中，不會重復(fù)出現(xiàn)模式相同的噪聲，因此從根本上解決了偽影問題。此外，我們在Rician噪聲期望保真項中引入了vNLM濾波，從而能夠降低原始信號的噪聲波動，提高模型擬合精度。本方法也可看成是一種雙域的去噪模型，結(jié)合了空間域的非局部均值濾波和變換域的HOSVD去噪，這種雙域的去噪模型能夠優(yōu)勢互補，從而達到更好的去噪效果。

我們進行了仿真實驗和真實實驗驗證，并且與4種典型的DW圖像去噪算法進行了比較，其中GL-HOSVD算法的去噪性能處于國際領(lǐng)先水平。實驗結(jié)果表明本方法能夠有效降低圖像噪聲，保留圖像細節(jié)結(jié)構(gòu)，并且不會引入條形偽影，其去噪性能無論在定量還是定性方面都是最優(yōu)的。與其他優(yōu)化算法類似，本方法的高水平去噪性能依賴于參數(shù)的設(shè)置。由于PSNR/SSIM最高的DW圖像產(chǎn)生的FA圖的RMSE并不是最小的，因此我們同時考慮了DW 圖像的PSNR，SSIM 和FA 圖的RMSE（FA-RMSE）3個指標對所提算法進行參數(shù)優(yōu)化。

本方法是以二維形式進行去噪的，可擴展到三維形式用以處理各向同性的DW數(shù)據(jù)。本方法同樣適用于多通道采集得到的noncentral-χ分布的DW噪聲圖像，只需要將Rician噪聲校正模型換成noncentral-χ分布下的噪聲信號期望即可。目前，該方法只適用于空間不變的噪聲水平，接下來我們可以研究如何對本方法進行改進，從而解決噪聲水平隨空間變化的DW圖像。

綜上所述，本文提出一種新穎的HOSVD 去噪方法，建立了基于HOSVD稀疏約束和Rician噪聲校正的去噪模型。所提去噪方法可以直接處理帶有Rician噪聲的DW圖像，不需要借助VST變換。此外，本方法還從根本上解決了去噪過程中引入偽影的問題。實驗結(jié)果證明了該方法的有效性，其去噪性能比采用了VST變換以及結(jié)構(gòu)相似塊的GL-HOSVD方法更好。