基于級聯(lián)卷積神經(jīng)網(wǎng)絡(luò)的非笛卡爾磁共振重建

覃 芹，張 利，2，張達(dá)敏，蔡宇佳，冉文兵

（1 貴州大學(xué) 大數(shù)據(jù)與信息工程學(xué)院，貴陽 550025；2 貴州大學(xué) 省部共建公共大數(shù)據(jù)國家重點(diǎn)實(shí)驗(yàn)室，貴陽 550025）

0 引言

磁共振成像（Magnetic Resonance Imaging，MRI）利用核磁共振原理進(jìn)行圖像重建，具有無放射性損害、對比度高等特點(diǎn)，現(xiàn)已廣泛應(yīng)用在醫(yī)學(xué)成像上，是臨床診斷中不可或缺的技術(shù)工具［1］。MRI 過程由數(shù)據(jù)采集過程和重建過程組成，需用到的數(shù)據(jù)采集方式有笛卡爾采樣和非笛卡爾采樣。其中，非笛卡爾采樣是加速磁共振成像的一種重要方法，能夠極大地減少數(shù)據(jù)采集時(shí)間。

傳統(tǒng)的非笛卡爾重建，主要有NUFFT（Non-Uniform Fast Fourier Transform）算法和網(wǎng)格化（Gridding）算法。Gridding 算法是1975年由Brouw［2］最先提出并應(yīng)用在天文數(shù)據(jù)處理上，后經(jīng)發(fā)展被廣泛應(yīng)用在核磁共振成像和計(jì)算機(jī)層析等成像中。1993 年，Dutt 等人［3］全面研究了NDFT 問題的4 種格式，給出了快速算法使用插值公式將數(shù)據(jù)值從非均勻坐標(biāo)點(diǎn)轉(zhuǎn)至均勻的坐標(biāo)點(diǎn)處，并使用標(biāo)準(zhǔn)的傅里葉變換達(dá)到快速計(jì)算的目的，這是NUFFT首次被提出［3-4］，此后發(fā)展至nufftt 可由tfkbnufft 和torchnufft 實(shí)現(xiàn)。在實(shí)現(xiàn)上，網(wǎng)格化算法和nufft 算法都是卷積插值的過程，即利用合適的卷積函數(shù)與非笛卡爾采集到的磁共振數(shù)據(jù)進(jìn)行卷積，將其轉(zhuǎn)換到均勻分布的笛卡爾網(wǎng)格上。該過程存在大量冗余的計(jì)算量，同時(shí)也由于卷積函數(shù)選擇的問題還存在較大誤差，冗余計(jì)算和誤差是導(dǎo)致非笛卡爾成像速度慢和質(zhì)量差的主要原因。

近年來，不少學(xué)者將深度學(xué)習(xí)應(yīng)用在欠采樣的MRI 重建領(lǐng)域中，與傳統(tǒng)的基于模型的方法有所不同，基于深度學(xué)習(xí)的方法運(yùn)用數(shù)據(jù)驅(qū)動技術(shù)，只要網(wǎng)絡(luò)設(shè)計(jì)得當(dāng)就可以通過大量數(shù)據(jù)訓(xùn)練實(shí)現(xiàn)非笛卡爾重建。Schlemper 等人［5］提出了一個(gè)深度級聯(lián)神經(jīng)網(wǎng)絡(luò)，并使用數(shù)據(jù)一致性模塊（Data Consistency，DC）來保障數(shù)據(jù)的保真度。隨后榮楚譽(yù)［6］在級聯(lián)網(wǎng)絡(luò)的基礎(chǔ)上，通過添加非均勻變換層，將級聯(lián)網(wǎng)絡(luò)推廣至非笛卡爾情況。2018 年，Zhu 等人［7］提出了一種基于流行逼近的自動轉(zhuǎn)換（Automated Transform by Mainfold Approximation，AUTOMAP）的成像方法，提出一種統(tǒng)一的圖像重建框架，該方法難以找到契合的流形估計(jì)，導(dǎo)致模型難以收斂，可能存在過擬合。

為了提高非笛卡爾磁共振重建質(zhì)量，本文提出了一種深度神經(jīng)網(wǎng)絡(luò)—級聯(lián)卷積神經(jīng)網(wǎng)絡(luò)［8-9］，利用密度補(bǔ)償來糾正非笛卡爾k 采樣空間的不均勻加權(quán)，利用forward nufft 運(yùn)算子和adjoint nufft 運(yùn)算子交替實(shí)現(xiàn)測量空間和圖像空間的交替修正問題。其中，圖像空間由3 個(gè)卷積層和1 個(gè)卷積去噪自編碼器組成，本方法的測量空間能夠充分利用原始k空間數(shù)據(jù)，圖像空間利用卷積去噪自編碼器能夠?qū)Ｗ⒂诟哳l信息的學(xué)習(xí)。

1 相關(guān)工作

磁共振是一項(xiàng)生物磁學(xué)核自旋成像技術(shù)，通過施加一個(gè)射頻電磁波來輻射被檢測物體，以此產(chǎn)生磁共振信號，該磁共振信號就是k空間數(shù)據(jù)。研究中將利用計(jì)算機(jī)來分析磁共振信號，進(jìn)而得到磁共振圖像。因此對于磁共振成像主要分為2 個(gè)階段：數(shù)據(jù)采集階段和重建階段。其中，數(shù)據(jù)采集階段就是k空間數(shù)據(jù)的填充過程，重建階段是實(shí)現(xiàn)從頻率域數(shù)據(jù)到圖像的過程。

1.1 問題建模

數(shù)據(jù)采集階段由采集設(shè)備決定，在笛卡爾情況下，得到的完全采樣的k空間，將其IFFT 即可轉(zhuǎn)換為磁共振圖像，在非笛卡爾情況下，k空間數(shù)據(jù)不是規(guī)則分布，需要利用數(shù)據(jù)重采樣技術(shù)將其轉(zhuǎn)換到笛卡爾網(wǎng)格上后再使用IFFT 進(jìn)行成像。完全采樣的k空間數(shù)據(jù)和磁共振圖像之間的關(guān)系由式（1）來描述：

其中，x（i，j）是需要重建的目標(biāo)圖像，yt(ki，kj) 是檢測到的MRI 信號，是頻率空間的表示。從式（1）可知，重建的圖像和檢測到的信號是一對傅里葉關(guān)系對。因此當(dāng)MRI 信號已知時(shí)，圖像重建的過程可由式（2）來描述：

由于MRI 信號數(shù)據(jù)是復(fù)數(shù)形式，無法對其進(jìn)行直接展示，所以在數(shù)據(jù)處理中一般是提取復(fù)數(shù)的幅值和相位數(shù)據(jù)，再從幅值和相位數(shù)據(jù)中重構(gòu)出空間圖像。對于k空間數(shù)據(jù)，其中心部分對應(yīng)于圖像輪廓，高頻部分對應(yīng)圖像的細(xì)節(jié)部分，因此目前的圖像重建挑戰(zhàn)即是如何充分利用高頻信息得到高質(zhì)量的圖像。

為了逼近更真實(shí)的磁共振重建場景，提高數(shù)據(jù)重建速度，提出非笛卡爾采樣條件，但是在該采樣條件下不能直接使用IFFT，需要利用數(shù)據(jù)重采樣將MRI 信號轉(zhuǎn)換到笛卡爾網(wǎng)格上。數(shù)據(jù)重采樣就是卷積插值的過程，其目的是將非笛卡爾采樣的MRI 信號卷積插值到均勻分布的網(wǎng)格上，該算法依賴于卷積函數(shù)的選取，其數(shù)學(xué)實(shí)現(xiàn)可由式（3）來描述：

其中，M*表示卷積插值后的數(shù)據(jù)；(ki，kj) 是測量的非笛卡爾采樣數(shù)據(jù)；yt(ki，kj) 為笛卡爾采樣數(shù)據(jù)；S（ki，kj）為非笛卡爾采樣函數(shù)，本文中探討的非笛卡爾采樣是徑向采樣，因此S（ki，kj）為徑向采樣函數(shù)；C(ki，kj) 是網(wǎng)格卷積函數(shù)。

1.2 相關(guān)算法

非笛卡爾磁共振重建常用的算法有網(wǎng)格化（Gridding）和非均勻傅里葉變換（NUFFT）。一般情況下，Gridding 在射電天文成像中應(yīng)用較為廣泛，NUFFT 算法則更多應(yīng)用在磁共振圖像重建中。而NUFFT 算法主要分為2 個(gè)運(yùn)算子。一個(gè)是前向（Forward）NUFFT 運(yùn)算子，前向運(yùn)算的目的是模擬成像過程中的非均勻頻率域數(shù)據(jù)。Forward NUFFT運(yùn)算子的操作流程如圖1 所示。另一個(gè)是伴隨（Adjoint）NUFFT 運(yùn)算子，這是一個(gè)成像過程，是將非均勻的數(shù)據(jù)變換到規(guī)則的網(wǎng)格上后、再進(jìn)行成像。

圖1 Forward NUFFT 運(yùn)算子的操作流程Fig. 1 Operation flow of Forward NUFFT operator

根據(jù)迭代重建的原理，可將非笛卡爾欠采樣磁共振由式（7）描述：

其中，x是待重建的目標(biāo)圖像；y是原始k空間數(shù)據(jù)；Fs表示系統(tǒng)矩陣；Fsx是非笛卡爾欠采樣數(shù)據(jù)；C是網(wǎng)格卷積函數(shù)；‖F(xiàn)sx*C -y‖為數(shù)據(jù)保真項(xiàng)；?（x）為正則化函數(shù)，對于傳統(tǒng)的方法來說，?（x）通常是x特定變換域的l0或l1泛式。

在選取網(wǎng)格卷積函數(shù)時(shí)，要滿足函數(shù)時(shí)域和頻域都具有良好集中性的要求，常用的網(wǎng)格卷積函數(shù)有Gaussian 函數(shù)、Kaiser-Bessel 函數(shù)和長波球函數(shù)PSWF 等。目前在磁共振非笛卡爾重建中，常用的是具有較好靈活性的Kaiser-Bessel 函數(shù)，研究時(shí)用到的時(shí)域表達(dá)式和頻域表達(dá)式分別見式（6）、式（7）：

在執(zhí)行插值運(yùn)算時(shí)，要制定過采樣系數(shù)K和插值寬度W，每一個(gè)采樣點(diǎn)都有W2大小的窗口，從采樣點(diǎn)到插值窗口內(nèi)的每個(gè)均勻點(diǎn)的距離決定插值權(quán)值。過采樣系數(shù)決定最終成像的分辨率。二者的選擇直接影響算法的速度，并且隨著過采樣系數(shù)K和插值寬度W的增大，算法的計(jì)算復(fù)雜度也會變大，由此即會帶來硬件計(jì)算問題，但這樣成像的誤差就越小，圖像質(zhì)量也越高。綜上所述可知，隨著其值的增大，成像質(zhì)量有所提高，但是成像速度較慢，因此在研究數(shù)據(jù)重采樣算法時(shí)需要平衡成像速度和成像質(zhì)量，就要選擇合適的過采樣系數(shù)和插值寬度。

2 本文方法

非笛卡爾磁共振重建質(zhì)量由重建算法決定，目前的重建算法存在冗余計(jì)算和計(jì)算誤差大的問題，導(dǎo)致重建結(jié)果較差，同時(shí)現(xiàn)有的神經(jīng)網(wǎng)絡(luò)算法不能很好地提取高頻信息。為了實(shí)現(xiàn)高質(zhì)量的非笛卡爾磁共振重建，本文提出一種基于密度補(bǔ)償?shù)募壜?lián)卷積神經(jīng)網(wǎng)絡(luò)。

2.1 方法理論

本文提出的非笛卡爾欠采樣磁共振圖像重建方法是一種迭代重建方法，該網(wǎng)絡(luò)是基于密度補(bǔ)償?shù)募壜?lián)卷積神經(jīng)網(wǎng)絡(luò)。方法的整體流程可做闡述如下：

（1）將采集的k空間數(shù)據(jù)進(jìn)行非笛卡爾徑向欠采樣得到模型的輸入數(shù)據(jù)ync。

（2）將全采樣k空間數(shù)據(jù)進(jìn)行逆傅里葉變換，得到真實(shí)的重建圖像xGT。

（3）將ync和xGT作為模型的訓(xùn)練數(shù)據(jù)集對，經(jīng)過網(wǎng)絡(luò)I次迭代后得到網(wǎng)絡(luò)重建結(jié)果xncpdnet。推得的數(shù)學(xué)模型公式可寫為：

其中，i為網(wǎng)絡(luò)當(dāng)前迭代次數(shù)，i∈［1，I］；I為總的迭代次數(shù)，I ＝10；θi（i∈［1，I］）是每次迭代過程中的網(wǎng)絡(luò)參數(shù)?；诿芏妊a(bǔ)償?shù)纳窠?jīng)網(wǎng)絡(luò)依賴于密度補(bǔ)償來糾正k空間的不均勻加權(quán)，利用非均勻變換層來連接測量空間和圖像空間，其中測量空間和圖像空間交替修正以補(bǔ)償更多的細(xì)節(jié)信息，實(shí)現(xiàn)高質(zhì)量重建。這里給出網(wǎng)絡(luò)結(jié)構(gòu)設(shè)計(jì)如圖2 所示。對此擬展開研究分述如下。

圖2 網(wǎng)絡(luò)結(jié)構(gòu)圖Fig. 2 The structure diagram of the proposed network

（1）測量空間網(wǎng)絡(luò)。測量空間由數(shù)據(jù)一致性模塊（DC）組成，在每個(gè)步驟中引入原始測量數(shù)據(jù)，保證數(shù)據(jù)的一致性。DC 最早是由Schlemper 等人［5］應(yīng)用在欠采樣的MRI 中。測量空間可由式（9）進(jìn)行描述：

其中，y為原始全采樣k空間；yi為第i次迭代的全采樣k空間；Ω為采樣軌跡；MΩ為非笛卡爾徑向欠采樣。DC 模塊可以作為先驗(yàn)知識來指導(dǎo)網(wǎng)絡(luò)的學(xué)習(xí)。

（2）圖像空間網(wǎng)絡(luò)。該部分網(wǎng)絡(luò)由3 個(gè)3×3conv2D組成，由于卷積是對實(shí)值進(jìn)行操作，故而在輸入數(shù)據(jù)傳送至該模塊前將復(fù)數(shù)表現(xiàn)形式的k空間提取實(shí)部和虛部，再將其傳入新的通道數(shù)中。這就使得傳入圖像網(wǎng)絡(luò)空間數(shù)據(jù)的通道數(shù)比此前要多2 個(gè)通道數(shù)。該部分用于提取從測量空間到圖像空間映射的特征。

（3）卷積去噪自編碼器。自編碼器由編碼層和解碼層組成。編碼器旨在提取特征，解碼器是將特征通過與編碼層相反的網(wǎng)絡(luò)結(jié)構(gòu)還原為原始數(shù)據(jù)。通過借助稀疏編碼的思想，使用稀疏的高階特征重新組合來重構(gòu)自身，達(dá)到輸出與原始數(shù)據(jù)盡可能相似的目的。但是在實(shí)際圖像重建中存在重構(gòu)誤差，因此在自編碼器的基礎(chǔ)上，可以添加噪聲來約束自編碼器重構(gòu)。這可使得編碼器學(xué)習(xí)提取最重要的特征并學(xué)習(xí)輸入數(shù)據(jù)中更加魯棒的表示。

2.2 損失函數(shù)

磁共振重建是一個(gè)迭代重建的過程，從數(shù)學(xué)角度看是一個(gè)不斷最小化損失函數(shù)的過程。這一過程可以被描述為式（11）～式（13）：

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

3.1 實(shí)驗(yàn)數(shù)據(jù)

本文數(shù)據(jù)使用的是fastMRI 數(shù)據(jù)集（https：／／fastmri.med.nyu.edu／），該數(shù)據(jù)集是Facebook AI Research（FAIR）和NYU Langone Health 之間的合作研究項(xiàng)目，目的是研究使用AI 提高M(jìn)RI 成像速度，并且從采樣不足的k空間數(shù)據(jù)中高效恢復(fù)出圖像，AI 進(jìn)行MRI 圖像重建可以使MRI 更廣泛應(yīng)用在醫(yī)療成像領(lǐng)域上，且符合更逼真的重建場景。fastMRI數(shù)據(jù)集包括2 種類型：膝蓋磁共振圖像和腦部磁共振圖像。數(shù)據(jù)集包括訓(xùn)練集、驗(yàn)證集以及masked 測試數(shù)據(jù)集。本實(shí)驗(yàn)采用的是單線圈膝蓋數(shù)據(jù)集，圖像尺寸大小為320×320，平面分辨率為0.5 mm×0.5 mm，切片厚度為3 mm，切片之間無縫隙。訓(xùn)練集總共有973 卷（973 個(gè)h5 文件，34 742 張圖片），驗(yàn)證集共有199 卷（199 個(gè)h5 文件，135 張圖片）。

3.2 欠采樣模式

為了適應(yīng)更真實(shí)的磁共振成像場景，提高數(shù)據(jù)的掃描速度，本文的欠采樣模式為非笛卡爾徑向欠采樣。圖3 顯示的是tfkbnufft 不同采樣模式下圖像重建的結(jié)果。圖3（a）～（c）中，左邊圖像是k-空間可視化，右邊圖像是在圖像域上的可視化。

圖3 徑向欠采樣模式下k 空間與重建圖像的可視化Fig. 3 Visualization of k-space and reconstructed images in radial undersampling mode

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

為了驗(yàn)證本方法，采用每張圖像重建時(shí)間T和重建質(zhì)量作為評價(jià)指標(biāo)。重建質(zhì)量指標(biāo)為峰值信噪比和結(jié)構(gòu)相似性。其中，峰值信噪比（PSNR）的定義是重建圖像與原始圖像之間對比質(zhì)量的好壞，PSNR越高、表明重建質(zhì)量越好；結(jié)構(gòu)相似性（SSIM）的定義是比較重建圖像與原始圖像之間的相似性，主要從亮度、對比度和結(jié)構(gòu)三個(gè)方面來進(jìn)行比較，其值范圍為［0，1］，即值越接近1、就表示重建圖像越接近于原始圖像。

3.4 實(shí)驗(yàn)環(huán)境與參數(shù)配置

本文基于密度補(bǔ)償?shù)纳窠?jīng)網(wǎng)絡(luò)模型是基于Tensorflow 2.7.0 版本實(shí)現(xiàn)的，顯卡型號為GPU RTX 3090 * 3，其中顯存大小為24 GB。CPU 型號為42核Intel（R）Xeon（R）Gold 6330 CPU ＠ 2.00 GHz，內(nèi)存大小為180 GB。

網(wǎng)絡(luò)迭代次數(shù)I ＝10，測量空間網(wǎng)絡(luò)由數(shù)據(jù)一致塊構(gòu)成，以殘差的方式不斷修正測量k空間和原始k空間數(shù)據(jù)之間的差距，圖像空間網(wǎng)絡(luò)由3 個(gè)卷積層組成，卷積大小都為3×3，由于卷積層只能處理實(shí)數(shù)，就需要k空間數(shù)據(jù)實(shí)部和虛部分別放入2 個(gè)通道，因此圖像網(wǎng)絡(luò)的輸入通道數(shù)為12，輸出通道數(shù)為10，核的初始化方式為glorot_uniform 方式。卷積去噪自編碼器用于提取高頻信息，該模塊由2 個(gè)卷積層和1 個(gè)反卷積層組成，其中卷積層卷積核大小為5×5。在網(wǎng)絡(luò)訓(xùn)練時(shí)，batch_size設(shè)置為一卷數(shù)據(jù)集、即為一個(gè)h5 文件，損失函數(shù)為mse，優(yōu)化器為Adamw，初始學(xué)習(xí)率為lr ＝10-5。訓(xùn)練輪次為200 輪。

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

為了驗(yàn)證方法在非笛卡爾欠采樣磁共振圖像重建中的有效性，以及卷積去噪自編碼器在高頻信息提取中的高效性，本文將比較3 種方法：

（1）傳統(tǒng)tfkbnufft 重建算法，是torchkbnufft 的進(jìn)一步發(fā)展，在重建速度上優(yōu)于torchkbnufft。

（2）基于密度補(bǔ)償?shù)纳窠?jīng)網(wǎng)絡(luò)（ncpdnet），通過在網(wǎng)絡(luò)中添加非均勻變換層重建非笛卡爾數(shù)據(jù)。

（3）基于密度補(bǔ)償?shù)募壜?lián)卷積神經(jīng)網(wǎng)絡(luò)，在ncpdnet 網(wǎng)絡(luò)的基礎(chǔ)上，為了提取高頻信息，使網(wǎng)絡(luò)具有更好的魯棒性，該網(wǎng)絡(luò)利用卷積去噪自編碼器學(xué)習(xí)更多細(xì)節(jié)信息，以下稱為ours。

運(yùn)行后，4 倍非笛卡爾徑向欠采樣模式下重建質(zhì)量和效率對比結(jié)果見表1，6 倍非笛卡爾徑向欠采樣模式下重建質(zhì)量和效率對比結(jié)果見表2，4 倍欠采樣模式下不同迭代次數(shù)下ours 重建質(zhì)量與效率對比結(jié)果見表3。圖4 展示的是在4 倍非笛卡爾欠采樣條件下，3 種方法的重建結(jié)果對比圖。

表1 4 倍非笛卡爾徑向欠采樣模式下重建質(zhì)量和效率對比Tab.1 Comparison of reconstruction quality and efficiency in the 4x non-cartesian radial undersampling mode

表2 6 倍非笛卡爾徑向欠采樣模式下重建質(zhì)量和效率對比Tab.2 Comparison of reconstruction quality and efficiency in the 6x non-cartesian radial undersampling mode

表3 4 倍欠采樣模式下不同迭代次數(shù)下ours 重建質(zhì)量與效率對比Tab.3 Comparison of reconstruction quality and efficiency of ours under different iterations in 4x undersampling mode

圖4 在4 倍欠采樣下3 種方法的重建結(jié)果圖Fig. 4 Reconstruction maps of the three methods under 4x undersampling mode

實(shí)驗(yàn)結(jié)果表明，基于神經(jīng)網(wǎng)絡(luò)的非笛卡爾欠采樣重建速度比傳統(tǒng)方法更快，主要是神經(jīng)網(wǎng)絡(luò)的方法基于數(shù)據(jù)驅(qū)動，不需要進(jìn)行插值運(yùn)算；ours 相比ncpdnet 方法重建質(zhì)量具有較大提升，主要是由于卷積去噪自編碼器能夠在密度補(bǔ)償展開網(wǎng)絡(luò)中發(fā)揮高頻信息提取作用，以重建出更多的圖像細(xì)節(jié)信息。同時(shí)，為了驗(yàn)證I ＝10 是最佳迭代次數(shù)，對比了不同迭代次數(shù)下ours 重建性能。雖然本文提出的方法相比于ncpdnet 重建時(shí)間略長，這是由于模型參數(shù)量比ncpdnet 大造成的，不過這個(gè)細(xì)微的重建時(shí)間差別并不影響本文方法的性能。上述結(jié)果均是在測試集表現(xiàn)的結(jié)果。

4 結(jié)束語

本文提出了一種基于密度補(bǔ)償?shù)募壜?lián)卷積神經(jīng)網(wǎng)絡(luò)，該網(wǎng)絡(luò)用于非笛卡爾欠采樣磁共振圖像重建，在本方法中，將MRI 圖像重建過程視為從非笛卡爾數(shù)據(jù)中迭代重建的過程，為了提取重建過程中的高頻信息部分，將卷積去噪自編碼器用于該重建任務(wù)中。同時(shí)對比了不同加速因子欠采樣模式下重建質(zhì)量和重建效率，通過實(shí)驗(yàn)證明了本文方法具有更強(qiáng)的細(xì)節(jié)信息提取能力，在4 倍欠采樣和6 倍欠采樣模式下，本文的方法能夠有效提高重建圖像的峰值信噪比和結(jié)構(gòu)相似性。本實(shí)驗(yàn)建立在單線圈成像場景上，為了驗(yàn)證卷積去噪自編碼器在密度補(bǔ)償神經(jīng)網(wǎng)絡(luò)中提取高頻信息的能力，未來的工作可以將其應(yīng)用到多線圈膝蓋以及腦部圖重建中。