核磁共振圖像重構(gòu)的即插即用算法

李金城，謝朝陽(yáng)，李金蘭，張 濤，鄒 健

長(zhǎng)江大學(xué) 信息與數(shù)學(xué)學(xué)院，湖北 荊州 434023

自20世紀(jì)70年代早期開始，核磁共振成像（magnetic resonance imaging，MRI）因其無創(chuàng)性且能有效地描述軟組織病變等優(yōu)點(diǎn)，徹底改變了放射學(xué)和醫(yī)學(xué)，成為廣泛應(yīng)用于臨床診斷的方法。限制MRI應(yīng)用的主要因素是其數(shù)據(jù)采集過程比其他技術(shù)慢。在不降低重構(gòu)質(zhì)量的前提下，減少數(shù)據(jù)采集量一直是MRI成像研究的重點(diǎn)[1-3]。在過去的幾十年里，基于稀疏先驗(yàn)的壓縮感知（compressed sensing，CS）利用圖像在某些變換域（如傅里葉變換、小波變換、梯度變換）中的稀疏性，精確地從少量隨機(jī)測(cè)量數(shù)據(jù)中重構(gòu)稀疏信號(hào)。壓縮感知核磁共振成像（CS-MRI）使得從高度欠采樣的K空間數(shù)據(jù)重構(gòu)核磁共振圖像成為可能，并成為多種重構(gòu)方法的基礎(chǔ)[4-5]。在MRI模型中，獲取的K空間數(shù)據(jù)y∈Rm通常被表示為：

其中，x∈Rn為待重構(gòu)圖像；ηe∈Rm為測(cè)量噪聲或干擾；A∈Rm×n(m＜n)是壓縮感知采樣矩陣。通常A=P×F,P∈Rm×n為欠采樣模板，如隨機(jī)采樣模板、徑向采樣模板、笛卡爾采樣模板等；F∈Rn×n為稀疏變換算子，如小波變換、傅里葉變換等。此時(shí)獲取的K空間數(shù)據(jù)y要比觀測(cè)核磁共振圖像x的規(guī)模小得多。從低維y中重構(gòu)唯一的高維圖像x可以通過求解如下優(yōu)化問題得到：

MRI中最常用的正則函數(shù)是L1正則，即φ(x)=L1范數(shù)的凸性可以保證問題（2）得到高效的求解。但是L1正則得到的解是真實(shí)信號(hào)的一個(gè)有偏估計(jì)，往往會(huì)低估x的高振幅分量[6]。非凸正則通常比L1正則更能有效地促進(jìn)稀疏和更準(zhǔn)確地估計(jì)x的高振幅分量。常用的非凸正則包括：光滑剪切絕對(duì)偏差（smoothly clipped absolute deviation，SCAD）正則，Lq范數(shù)φ(x)=‖x ‖q(0＜q＜1)等[7]。然而這些非凸正則，會(huì)使得問題（2）的目標(biāo)函數(shù)是非凸的，使其陷入次優(yōu)的局部最小值，增加了求解難度。Selesnick等人通過引入L1范數(shù)的Moreau包絡(luò)和最小最大凹罰函數(shù)，構(gòu)造了凸-非凸（convex-nonconvex，CNC）稀疏正則[6，8-9]，其表達(dá)式如下：

MRI重構(gòu)的另一個(gè)重要問題是重構(gòu)算法。近似點(diǎn)梯度下降（proximal gradient descent，PGD）是一種高效的求解非光滑問題（2）的一階算法，其基本原理是利用不同策略來處理目標(biāo)函數(shù)的光滑和非光滑部分，對(duì)光滑部分使用梯度迭代，對(duì)非光滑部分利用近似點(diǎn)算子求解[11-12]。近年來，隨著深度學(xué)習(xí)的興起，將傳統(tǒng)優(yōu)化算法與深度神經(jīng)網(wǎng)絡(luò)結(jié)合的數(shù)據(jù)驅(qū)動(dòng)算法得到快速發(fā)展。即插即用（plug-and-play，PnP）算法將用合適的去噪器來替代傳統(tǒng)稀疏優(yōu)化算法中的近似點(diǎn)算子[13-14]，其優(yōu)點(diǎn)是當(dāng)沒有足夠的數(shù)據(jù)進(jìn)行端到端訓(xùn)練時(shí)，傳統(tǒng)或者預(yù)先訓(xùn)練好的去噪器可以作為即插即用框架下的模塊化部件靈活插入使用。常用的去噪器包括傳統(tǒng)的手工設(shè)計(jì)去噪器，如全變差（total variation，TV）、三維塊匹配協(xié)同濾波（block-matching 3D，BM3D）[15]等。近年來，深度神經(jīng)網(wǎng)絡(luò)去噪器可以從大量自然圖像中學(xué)習(xí)到更為豐富的圖像先驗(yàn)，顯示了優(yōu)異的去噪性能[13，16]。基于深度學(xué)習(xí)的即插即用算法在圖像重構(gòu)、圖像去噪、語(yǔ)音去噪和運(yùn)動(dòng)設(shè)備等方向得到了廣泛的應(yīng)用[17-19]。

本文利用即插即用近似點(diǎn)梯度下降算法來求解帶凸-非凸正則的核磁共振重構(gòu)模型，其主要貢獻(xiàn)如下：

（1）給出了凸-非凸正則的近似點(diǎn)算子的迭代解形式，解決了采用凸-非凸正則的優(yōu)化問題求解過程中較為困難的部分，從而提高計(jì)算效率。

（2）提出求解基于凸-非凸正則的核磁共振重構(gòu)模型的近似點(diǎn)梯度下降算法，進(jìn)而利用去噪器替換近似點(diǎn)算子，得到求解基于凸-非凸正則的核磁共振重構(gòu)模型的即插即用的近似點(diǎn)梯度下降算法。其中即插即用算法可以靈活地引入更多先進(jìn)的去噪器，可以有效提升去噪效果。另外在MRI重構(gòu)中，采用凸-非凸正則能夠更精確地重構(gòu)圖像在頻域中的高頻分量，進(jìn)而能減少重構(gòu)圖像邊緣輪廓和細(xì)小紋理等細(xì)節(jié)信息的丟失。

1 相關(guān)工作

定義1對(duì)于任意的α＞0，對(duì)于給定的適當(dāng)、閉、凸函數(shù)g(x)，其近似點(diǎn)算子定義為：

其中，xi為x的第i個(gè)元素。從軟閾值的定義可以看出，當(dāng)xi的絕對(duì)值小于某一閾值α?xí)r，軟閾值函數(shù)會(huì)將xi的值變?yōu)?，因此軟閾值函數(shù)可理解為一個(gè)去噪器，參數(shù)α控制其去噪強(qiáng)度。

由近似點(diǎn)算子的定義，可得到近似點(diǎn)梯度下降算法的一般形式?？紤]復(fù)合優(yōu)化問題：

其中，f(x)凸且可微的，g(x)凸但不可微，且g(x)的近似點(diǎn)算子容易求得。

近似點(diǎn)梯度下降算法的基本迭代形式如下：

其中，?f(x)是f(x)的梯度函數(shù)，α是步長(zhǎng)，g的近似點(diǎn)算子proxαg如式（4）所定義。

2 凸-非凸稀疏正則的近似點(diǎn)算子

凸-非凸稀疏正則φb(x)的近似點(diǎn)算子定義如下：

由于φb(x)形式較為復(fù)雜，此時(shí)proxλφb(y)不再具有封閉解，本文將給出其迭代解。

證明 由近似點(diǎn)算子的定義可得Sb(x)的最優(yōu)值點(diǎn)為：

代入Sb(x)中，則有：

容易求得Sb(x)的梯度如下：

命題1得證。

命題2 φb(x)的近似點(diǎn)算子為：

證明 將式（8）中的目標(biāo)函數(shù)記為：

由近似點(diǎn)梯度下降算法，可得：

取步長(zhǎng)μ=1，命題2得證。

3 即插即用近似點(diǎn)梯度下降算法

3.1 近似點(diǎn)梯度下降算法

本節(jié)將考慮運(yùn)用近似點(diǎn)梯度下降算法求解基于凸-非凸稀疏正則的核磁共振圖像重構(gòu)模型（2）。令式（2）中,結(jié)合命題2，求解基于凸-非凸稀疏正則的近似點(diǎn)梯度下降（proximal gradient descent algorithm for convex-nonconvex sparse regularization，PGD-CNC）算法的主要步驟如下：

步驟1

步驟2

3.2 即插即用近似點(diǎn)梯度下降算法

此時(shí)若選擇兩個(gè)適當(dāng)?shù)娜ピ肫鱄σ1、Hσ2分別替換近似點(diǎn)算子、proxαλ‖·‖1，其中σ1、σ2為去噪?yún)?shù)，值越大去噪強(qiáng)度越大，則式（19）可改寫為：

綜合式（18）、（20），可得到求解基于凸-非凸稀疏正則的即插即用近似點(diǎn)梯度下降（plug-and-play proximal gradient descent algorithm for convex-nonconvex sparse regularization，PnP-PGD-CNC）算法如算法1所示。

算法1 MRI圖像重構(gòu)的PnP-PGD-CNC算法

即插即用的近似點(diǎn)梯度下降算法與近似點(diǎn)梯度下降算法的區(qū)別：前者可將預(yù)先訓(xùn)練好的去噪器直接插入到傳統(tǒng)近似點(diǎn)梯度下降算法中的一個(gè)迭代步中。當(dāng)去噪器滿足一些簡(jiǎn)單條件，如非擴(kuò)張、收縮時(shí)，可以保證算法的收斂性[17，19]。

算法1中計(jì)算復(fù)雜度較高的是z(k+1)的迭代，涉及到多個(gè)矩陣相乘，但在實(shí)際實(shí)驗(yàn)中都可以進(jìn)行快速計(jì)算。由于P是欠采樣矩陣，則PTP為矩陣元素間的乘積。若F為傅里葉變換，則有FT=F-1,F和FT在實(shí)驗(yàn)中可通過快速傅里葉變換（fast Fourier transform，F(xiàn)FT）和逆快速傅里葉變換（inverse fast Fourier transform，IFFT）實(shí)現(xiàn)。

4 數(shù)值實(shí)驗(yàn)與結(jié)果分析

4.1 實(shí)驗(yàn)設(shè)置

壓縮感知核磁共振成像通過下采樣獲取更少的數(shù)據(jù)來加速核磁共振成像。本實(shí)驗(yàn)選取了三種不同的核磁共振圖像：大腦、腦部血管和上半身。所有圖像灰度值取值范圍為0～1，尺寸設(shè)置為256×256像素；ηe為測(cè)量噪聲且滿足ηe～N(0,σeIk)，噪聲級(jí)別為σe=15/255；P是欠采樣模板，本實(shí)驗(yàn)應(yīng)用了采樣率均為30%的欠采樣隨機(jī)采樣模板、徑向采樣模板和笛卡爾采樣模板；F為傅里葉算子。待重構(gòu)圖像和采樣模板如圖1所示。

為了驗(yàn)證PnP-PGD-CNC算法的有效性，實(shí)驗(yàn)將ADMM-Net[20]、BM3D-MRI[21]、基于L1正則的近似點(diǎn)梯度下降算法（PGD-L1）[17]、基于L1正則的即插即用近似點(diǎn)梯度下降算法（PnP-PGD-L1）[17]以及本文提出的基于凸-非凸稀疏正則的近似點(diǎn)梯度下降算法（PGD-CNC）、基于凸-非凸稀疏正則的即插即用的近似點(diǎn)梯度下降算法（PnP-PGD-CNC）的核磁共振重構(gòu)圖像進(jìn)行對(duì)比。其中ADMM-Net[20]是將傳統(tǒng)的交替方向乘子法（alternating direction method of multipliers，ADMM）展開到深度神經(jīng)網(wǎng)絡(luò)中；BM3D-MRI[21]是使用非局部的BM3D圖像模型進(jìn)行MRI重構(gòu)；其他四種算法的主要迭代步驟如表1所示。

表1中PGD-L1算法為正則項(xiàng)φ(x)=‖x‖1時(shí)，使用PGD算法得出的迭代解，其解由軟閾值函數(shù)給出[17]；PnP-PGD-L1算法為文獻(xiàn)[17]中提到的，由上述PGD算法與PnP算法相結(jié)合得出的迭代解；PGD-CNC算法即由式（18）、（19）給出；PnP-PGD-CNC算法即本文所提出的算法（算法1）。

表1 算法的迭代式Table 1 Iterative formula of algorithms

在PnP-PGD-L1和PnP-PGD-CNC算法中選取BM3D[15]和五種神經(jīng)網(wǎng)絡(luò)去噪器作為PnP框架下不同的去噪器。表2簡(jiǎn)要羅列了各種去噪器的特點(diǎn)。

表2 PnP框架下的去噪器Table 2 Denoiser under PnP framework

特別地，文獻(xiàn)[25]中提到DnCNN需要為每個(gè)噪聲水平分別學(xué)習(xí)一個(gè)模型，其他的去噪器通過單個(gè)模型可以處理各個(gè)噪聲水平的噪聲圖片。因此，在PnP-PGDCNC算法中選用DnCNN去噪器時(shí)，去噪?yún)?shù)σ1=5σe/3=25/255,σ2=σe/3=5/255,在其余算法中σ1=σ2=σe=15/255。本實(shí)驗(yàn)中，參數(shù)設(shè)置均是在大量的實(shí)驗(yàn)下選擇圖像恢復(fù)效果較好時(shí)的參數(shù)，迭代次數(shù)均設(shè)置為各算法在50次內(nèi)效果最好時(shí)的迭代次數(shù)。五種神經(jīng)網(wǎng)絡(luò)去噪器由多種不同的自然數(shù)據(jù)集訓(xùn)練得到[22-25]。實(shí)驗(yàn)操作都是在Intel?CoreTMi7-8750HM、2.20 GHz CPU、16 GB內(nèi)存，NVIDIA GeForce GTX 1060 6 GB顯卡的PC上進(jìn)行，編程環(huán)境使用Python3.8，深度學(xué)習(xí)框架使用Pytorch。本文的主要實(shí)驗(yàn)代碼、參數(shù)設(shè)置和結(jié)果可以從網(wǎng)站鏈接https：//github.com/zj15001/PNP-PGD-CNC下載。

本文實(shí)驗(yàn)結(jié)果由兩方面進(jìn)行呈現(xiàn)，第一方面是視覺效果圖像對(duì)比，第二方面是數(shù)值結(jié)果比較。

為了從視覺效果上直觀的評(píng)價(jià)重構(gòu)效果，本文分別給出了表1中不同算法的重構(gòu)效果圖，以及重構(gòu)后效果圖與待重構(gòu)核磁共振圖像的差值圖像。為方便更好地顯示差異，選取圖像的兩個(gè)局部細(xì)節(jié)進(jìn)行放大顯示，并將差值圖像的灰度值取值范圍設(shè)置為0～0.2。

數(shù)值結(jié)果比較采用峰值信噪比(PSNR)的數(shù)值大小來作為評(píng)價(jià)核磁共振圖像重構(gòu)效果的重要指標(biāo)，其定義如下：

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

首先，采用隨機(jī)采樣模板對(duì)大腦圖像進(jìn)行采樣后加噪聲，再運(yùn)用表1中的四種算法進(jìn)行圖像重構(gòu)，具體重構(gòu)圖像和差值圖像如圖2所示。其中，每張子圖的第一排為重構(gòu)圖像，第二排為重構(gòu)圖像與待重構(gòu)圖像的差值圖像，每排下方的兩個(gè)小圖則是將兩個(gè)局部細(xì)節(jié)進(jìn)行放大顯示；PnP-PGD-L1和PnP-PGD-CNC算法的重構(gòu)結(jié)果中顯示了部分去噪器（BM3D、FFDNet、IRCNN、DRUNet）的重構(gòu)圖像和差值圖像。

由圖2中各算法的重構(gòu)圖像及其差值圖像的對(duì)比可以看出，PGD-L1算法的重構(gòu)效果最差，圖像分辨率較低，保留了較少的邊緣和細(xì)節(jié)信息。PnP-PGD-L1和PGD-CNC算法相較PGD-L1算法都有不同程度的性能提升。本文提出的PnP-PGD-CNC算法整體要比其他算法的重構(gòu)效果更好，其中效果最好的是采用了DRUNet去噪器的圖2（j）。其左上角方框中的邊緣線條和紋路更加清晰，右下角方框中的棱角部分更加分明，保留了更好的邊緣輪廓信息，且差值圖的顏色明顯變深，表明了重構(gòu)圖像與原始圖像的一致性較高。從數(shù)值結(jié)果來看，PGD-L1算法的重構(gòu)圖像（圖2（a））峰值信噪比為23.78 dB，而在PnP-PGD-CNC框架下使用DRUNet去噪器的重構(gòu)圖像（圖2（j））峰值信噪比為31.60 dB，顯然PnP-PGD-CNC的視覺效果和數(shù)值結(jié)果更好。

最后，表3給出對(duì)于不同圖像使用不同算法進(jìn)行圖像重構(gòu)的數(shù)值比較結(jié)果。實(shí)驗(yàn)數(shù)據(jù)表明，對(duì)于本文給出的三種不同采樣模板，PnP-PGD-L1算法和PGD-CNC算法的重構(gòu)效果都要明顯高于PGD-L1算法，分別突出了即插即用算法和凸-非凸稀疏正則的有效性。同時(shí)對(duì)比PnP-PGD-L1算法和PnP-PGD-CNC算法中使用不同去噪器時(shí)，神經(jīng)網(wǎng)絡(luò)去噪器的重構(gòu)效果要優(yōu)于其他去噪器，能夠得到更高的峰值信噪比值，其中峰值信噪比值最高的是選用DRUNet去噪器。當(dāng)PnP-PGD-CNC算法選用DRUNet去噪器時(shí)，較PGD-L1算法峰值信噪比值平均提升了6.26 dB。同時(shí)從PSNR值的大小也可以看出PnP-PGD-CNC算法的重構(gòu)效果明顯優(yōu)于ADMMNet和BM3D-MRI算法。

表3 圖像重構(gòu)的峰值信噪比值Table 3 Peak signal-to-noise ratio of image reconstruction 單位：dB

5 結(jié)束語(yǔ)

本文提出了一種基于凸-非凸稀疏正則和即插即用近似點(diǎn)梯度下降的核磁共振圖像重構(gòu)算法，主要有以下特點(diǎn)：

（1）該算法是基于凸-非凸稀疏正則的，可以保證目標(biāo)函數(shù)的整體凸性，并且能夠準(zhǔn)確地估計(jì)高振幅分量。本文給出了該正則項(xiàng)的近似點(diǎn)算子的迭代解形式，這為之后基于凸-非凸稀疏正則的優(yōu)化問題的計(jì)算提供了便利。

（2）該算法中融入的即插即用算法，是將傳統(tǒng)優(yōu)化算法與深度學(xué)習(xí)相結(jié)合，與傳統(tǒng)圖像重構(gòu)方法不同，它可以通過去噪器隱式定義即插即用先驗(yàn)。該框架具有足夠的靈活性，可以與各種最先進(jìn)的去噪器兼容。同時(shí)可以從理論上保證算法的收斂性，體現(xiàn)出其便捷、高效等優(yōu)點(diǎn)。

為驗(yàn)證本文所提算法的有效性，與另外幾種算法的重構(gòu)結(jié)果進(jìn)行比較。從視覺效果圖像和數(shù)值(PSNR)結(jié)果比較都可以得出：本文所提算法對(duì)邊緣輪廓和細(xì)小紋理的處理具有較為顯著的性能提升，并且可以得到更高的峰值信噪比。

進(jìn)一步的工作研究主要分為三方面：一是將繼續(xù)探究本文所提出算法的收斂性理論證明；二是考慮設(shè)計(jì)更為高效的基于深度學(xué)習(xí)的去噪器，以達(dá)到更好的去噪效果；三是將凸-非凸稀疏正則和即插即用算法融入到其他的凸優(yōu)化算法進(jìn)行圖像重構(gòu)的研究。