吳春燕，李雨澤，丁海艷，陳慧軍

清華大學(xué) 醫(yī)學(xué)院 生物醫(yī)學(xué)工程系，北京 100084

引言

心臟磁共振T1量化成像可以指示心肌組織的一系列病理改變，特別是心臟的炎癥、急（慢）性心肌瘢痕以及纖維化引起的細(xì)胞外間質(zhì)重組[1-2]。臨床上，改良Look-Locker 反轉(zhuǎn)恢復(fù)（Modified Look-Locker Inversion Recovery，MOLLI）序列是標(biāo)準(zhǔn)的心肌T1參數(shù)量化序列，但是掃描過(guò)程中心臟及呼吸運(yùn)動(dòng)將導(dǎo)致圖像質(zhì)量降低，需對(duì)采集的不同幀圖像進(jìn)行配準(zhǔn)以提高圖像質(zhì)量及T1量化精度。

國(guó)際上已有許多研究者對(duì)心肌的T1量化提出了運(yùn)動(dòng)校正方法，Roujol等[3]提出一個(gè)基于光流法的局部非剛性配準(zhǔn)架構(gòu)，同時(shí)估計(jì)運(yùn)動(dòng)場(chǎng)和信號(hào)強(qiáng)度變化；El-Rewaidy等[4]利用一種活動(dòng)形狀模型分割得到的心肌輪廓實(shí)現(xiàn)配準(zhǔn)，提高了各心肌節(jié)段T1參數(shù)量化的準(zhǔn)確度和精度；Tilborghs等[5]提出一種結(jié)合數(shù)據(jù)驅(qū)動(dòng)初始化和基于模型的精細(xì)改良的配準(zhǔn)方法，實(shí)現(xiàn)了更加魯棒的T1參數(shù)量化。但是，這些傳統(tǒng)配準(zhǔn)方法需要多次迭代優(yōu)化，耗時(shí)長(zhǎng)，降低效率[6-7]。Fahmy等[8]提出一種基于深度學(xué)習(xí)的自動(dòng)心肌T1參數(shù)量化方法，然而，這一深度學(xué)習(xí)算法為監(jiān)督式的，需要人工標(biāo)注數(shù)據(jù)，且在自動(dòng)分割效果不佳時(shí)需要額外的后處理操作。

為了克服前述算法的缺點(diǎn)，本研究提出一種基于深度學(xué)習(xí)的心肌T1參數(shù)量化成像運(yùn)動(dòng)校正算法，無(wú)需傳統(tǒng)方法的多次迭代優(yōu)化，可大幅度降低運(yùn)算時(shí)間；同時(shí)，采用自監(jiān)督的訓(xùn)練方法，無(wú)需人工標(biāo)注的金標(biāo)準(zhǔn)圖像，大大降低了醫(yī)生的工作量，為大規(guī)模臨床應(yīng)用提供了可能。

1 材料與方法

1.1 數(shù)據(jù)集

實(shí)驗(yàn)數(shù)據(jù)來(lái)自47例健康志愿者（30例男性，17例女性，年齡30～64歲），納入標(biāo)準(zhǔn)：① 既往無(wú)心源性或系統(tǒng)性疾?。虎?既往CMR掃描無(wú)心臟形態(tài)和功能異常；③ 無(wú)磁共振檢查禁忌證。所有受試者均簽署了知情同意書(shū)。

實(shí)驗(yàn)使用3.0 T磁共振掃描儀（Achieva TX，Philips，荷蘭），32通道的心臟線圈，所有受試者在屏氣情況下進(jìn)行成像。采用MOLLI序列進(jìn)行掃描，每名志愿者采集1～3層圖像，掃描方位為左心室短軸位，每層圖像包含8幀T1加權(quán)圖像。序列的具體掃描參數(shù)如下：重復(fù)時(shí)間2.59 ms；回波時(shí)間1.31 ms；觸發(fā)延遲（Trigger Delay，TD）時(shí)間為舒張末期；完成一次掃描的時(shí)間為8～34 s；翻轉(zhuǎn)角35°；掃描視野30 cm×15 cm；分辨率1.25 mm×1.25 mm；層厚8 mm。

1.2 基于深度學(xué)習(xí)的心肌T1參數(shù)量化成像運(yùn)動(dòng)校正算法

本方法的總體流程主要包括兩個(gè)部分：① 利用基于深度學(xué)習(xí)的自監(jiān)督配準(zhǔn)網(wǎng)絡(luò)將MOLLI圖像中除固定圖像以外的所有圖像配準(zhǔn)到固定圖像上；② 對(duì)配準(zhǔn)后的T1加權(quán)圖像進(jìn)行參數(shù)擬合得到T1量化圖像。

1.2.1 基于自監(jiān)督深度學(xué)習(xí)的非剛性配準(zhǔn)算法

配準(zhǔn)流程圖如圖1所示。圖像配準(zhǔn)通常是將一幅移動(dòng)（或待配準(zhǔn)）圖像配準(zhǔn)到固定（或參考）圖像上。令m、f分別表示定義在2D空間域Ω∩R2上的移動(dòng)、固定圖像，利用卷積神經(jīng)網(wǎng)絡(luò)（Convolutional Neural Network，CNN）對(duì)配準(zhǔn)運(yùn)動(dòng)場(chǎng)u建模，即gθ(f,m)=u，其中g(shù)為神經(jīng)網(wǎng)絡(luò)，θ是可學(xué)習(xí)參數(shù)。網(wǎng)絡(luò)可學(xué)習(xí)每一組輸入對(duì)(m,f)之間的配準(zhǔn)變換；接著，運(yùn)動(dòng)場(chǎng)利用空間變換函數(shù)(o)作用到移動(dòng)圖像m上，可得到運(yùn)動(dòng)校正后的圖像m'。具體地，對(duì)于成像空間中每個(gè)像素p，經(jīng)運(yùn)動(dòng)場(chǎng)u的空間交換(°)得到的m′(p)是其周圍8個(gè)像素點(diǎn)的線性插值，見(jiàn)式(1)。

圖1 基于自監(jiān)督深度學(xué)習(xí)的非剛性配準(zhǔn)算法

其中，Ζ8(p')表示p'的8像素鄰域，p'=p+u(p)。

本方法中的自監(jiān)督體現(xiàn)在，由于固定圖像來(lái)源于直接采集到的圖像，無(wú)需人為標(biāo)注或額外的處理，經(jīng)過(guò)運(yùn)動(dòng)場(chǎng)作用后的圖像m'與固定圖像進(jìn)行損失函數(shù)計(jì)算并反向傳播優(yōu)化網(wǎng)絡(luò)，達(dá)到自監(jiān)督學(xué)習(xí)的目的。

1.2.2 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)

圖2為本方法使用的神經(jīng)網(wǎng)絡(luò)的結(jié)構(gòu)。該網(wǎng)絡(luò)采用Unet結(jié)構(gòu)[9]，主要由編碼器和解碼器兩個(gè)部分組成：編碼器利用卷積層增加CNN的感受野，降低CNN的空間尺寸，同時(shí)學(xué)習(xí)不同尺度下的圖像信息；解碼器采用上采樣恢復(fù)圖像尺寸，并利用跳躍連接融合編碼器和解碼器部分的多尺度信息。編碼和解碼階段都采用尺寸為3×3、步長(zhǎng)為2的卷積和Leaky ReLU激活函數(shù)（小于0部分的斜率為0.2）。

圖2 運(yùn)動(dòng)場(chǎng)估計(jì)的網(wǎng)絡(luò)結(jié)構(gòu)圖

1.2.3 損失函數(shù)

本方法的損失函數(shù)由兩部分構(gòu)成，一部分為自監(jiān)督損失函數(shù)，即圖1中的sim，sim衡量的是校正后圖像與固定圖像的相似度?？紤]到T1加權(quán)圖像不同幀之間對(duì)比度相差大，直接采用基于信號(hào)強(qiáng)度的相似性測(cè)度效果不佳[10]，im選擇對(duì)于對(duì)比度變化不敏感的互相關(guān)（Cross-Correlation，CC）和互信息（Mutual Information，MI）作為相似性測(cè)度，見(jiàn)式(2)～(3)。

P(a)是f中信號(hào)強(qiáng)度等于a的像素比例，P(b)是m'中信號(hào)強(qiáng)度等于b的像素比例，P(a,b)是f和m'信號(hào)強(qiáng)度的聯(lián)合分布。為了使損失函數(shù)可微分以便于神經(jīng)網(wǎng)絡(luò)的反向傳播，使用Parzen窗法近似公式中的概率分布[11]，見(jiàn)式(4)～(6)。

其中，窗函數(shù)W選擇高斯窗。

另外，為了對(duì)比采用基于信號(hào)強(qiáng)度的相似性測(cè)度的運(yùn)動(dòng)校正效果，實(shí)驗(yàn)也驗(yàn)證了sim取均方誤差（Mean Square Error，MSE）的情況，見(jiàn)式(7)。

綜上，網(wǎng)絡(luò)的總損失函數(shù)如式(9)所示。

本研究測(cè)試不同的相似性測(cè)度（MSE、CC、MI）以尋找最優(yōu)的損失函數(shù)。

1.2.4 訓(xùn)練與測(cè)試

網(wǎng)絡(luò)使用TensorFlow 2.1.0框架實(shí)現(xiàn)，優(yōu)化器為Adam，學(xué)習(xí)率設(shè)置為1e-4，beta1=0.9，beta2=0.999；采用的計(jì)算服務(wù)器配置為I7-6950K CPU, NVIDIA TITAN Xp 12 GB GPU及32 GB內(nèi)存。在本研究中，固定圖像為MOLLI序列采集到的T1加權(quán)圖像中的第8幀，因?yàn)榇藭r(shí)心肌-血池對(duì)比度較高，移動(dòng)圖像則為其他幀的圖像。47例受試者數(shù)據(jù)中，36例（2035張圖像）用于網(wǎng)絡(luò)訓(xùn)練，11例（650張圖像）用于測(cè)試。

此外，本文對(duì)比了一種傳統(tǒng)配準(zhǔn)方法：基于互相關(guān)的B樣條自由形變模型配準(zhǔn)算法[12]，由Matlab醫(yī)學(xué)圖像配準(zhǔn)工具箱（Medical Image Registration Toolbox，MIRT）實(shí)現(xiàn)。

1.3 T1參數(shù)擬合

圖像經(jīng)過(guò)配準(zhǔn)后，還需經(jīng)過(guò)T1擬合才可得到T1量化圖像。如圖3所示，本實(shí)驗(yàn)采用的MOLLI序列由兩次ECG觸發(fā)的Look-Locker（LL）實(shí)驗(yàn)（LL1、LL2）組成，分別由三、五次單次激發(fā)讀出。每個(gè)LL實(shí)驗(yàn)中，反轉(zhuǎn)脈沖和TD之后，在遞增的反轉(zhuǎn)恢復(fù)時(shí)間（Inversion Time，TI）心臟舒張末期采集信號(hào)，LL1和LL2之間至少間隔4 s保證完全的信號(hào)恢復(fù)。

圖3 MOLLI序列信號(hào)采集示意圖

本文采用三參數(shù)模型進(jìn)行T1擬合[13]，見(jiàn)式(10)。

對(duì)于MOLLI序列，擬合得到的T1*經(jīng)過(guò)如下的修正后得到T1：T1=T1*(-B/A-1)。本文采用葉猛[14]提出的結(jié)合RD（Reduced Dimension Nonlinear Least Square with Rolarity Restorations）算法和LM（Levenberg-Marquardt）算法的擬合流程，即對(duì)運(yùn)動(dòng)校正后的T1加權(quán)圖像加mask去除背景后，先進(jìn)行RD擬合，再對(duì)擬合殘差過(guò)大的點(diǎn)用LM擬合。

1.4 評(píng)價(jià)指標(biāo)及統(tǒng)計(jì)學(xué)分析

本文評(píng)價(jià)了T1加權(quán)圖像的配準(zhǔn)效果及T1量化誤差。采用Dice相似系數(shù)（Dice Similarity Coefficient，DSC）和平均邊界誤差（Mean Boundary Error，MBE）作為量化評(píng)估配準(zhǔn)誤差的指標(biāo)。這兩個(gè)指標(biāo)的獲得需要在原始加權(quán)圖像上手動(dòng)勾畫出心肌內(nèi)、外輪廓。指標(biāo)具體計(jì)算可參考Tilborghs等[5]的方法。采用T1量化圖像的標(biāo)準(zhǔn)差（Standard Deviation，SD）[15]評(píng)估T1參數(shù)量化的準(zhǔn)確度，可利用T1擬合殘差和測(cè)量噪聲的SD近似得到待估計(jì)參數(shù)的協(xié)方差矩陣，從而得到T1參數(shù)量化的SD。除了SD 量化圖，本文還計(jì)算了心肌部分的平均SD。

本文使用Matlab R2020a軟件進(jìn)行統(tǒng)計(jì)學(xué)分析。DSC、MBE和心肌SD表示為±s；對(duì)不同方法測(cè)量得到的DSC、MBE和心肌SD進(jìn)行Wilcoxon符號(hào)秩檢驗(yàn)，P<0.05表示具有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

2.1 T1加權(quán)圖像配準(zhǔn)結(jié)果

經(jīng)過(guò)本文提出的深度學(xué)習(xí)算法（以CC和MI作為相似性測(cè)度）配準(zhǔn)后，T1加權(quán)圖像的心肌部分得到了較好的運(yùn)動(dòng)校正。T1加權(quán)圖像配準(zhǔn)效果如圖4所示。未配準(zhǔn)時(shí)（圖4a、 e、 i中藍(lán)色箭頭所示），心肌位置相對(duì)于固定圖像發(fā)生了明顯的運(yùn)動(dòng)偏移，采用深度學(xué)習(xí)算法配準(zhǔn)后，相似性測(cè)度選擇CC和MI時(shí)（圖4c、g、k和圖4d、h、l），心肌位置與固定圖像的基本一致。當(dāng)心肌和血池對(duì)比度很低時(shí)（圖4i、k、l），算法也能較好地配準(zhǔn)；相似性測(cè)度選擇MSE時(shí)，配準(zhǔn)后T1加權(quán)圖像出現(xiàn)明顯的變形、失真（圖4b、f、j中紅色箭頭所示）；當(dāng)心肌和血池對(duì)比度很低時(shí)（圖4j中紅色箭頭所示），失真情況更加嚴(yán)重。

圖4 三例未配準(zhǔn)和深度學(xué)習(xí)算法（MSE、CC、 MI）配準(zhǔn)后的T1加權(quán)圖像

定量評(píng)估算法的配準(zhǔn)準(zhǔn)確度結(jié)果如表1所示。深度學(xué)習(xí)（CC）配準(zhǔn)后DSC較原始未配準(zhǔn)顯著提高（P<0.05），而采用傳統(tǒng)算法和深度學(xué)習(xí)（MI）配準(zhǔn)后DSC指標(biāo)沒(méi)有顯著改善（P=0.13和P=0.98）。傳統(tǒng)算法、深度學(xué)習(xí)算法（CC、MI）配準(zhǔn)后MBE均顯著降低（P<0.05），且深度學(xué)習(xí)算法（CC、MI）的MBE顯著低于傳統(tǒng)算法（P<0.05）。而深度學(xué)習(xí)（MSE）配準(zhǔn)后，DSC和MBE均比未配準(zhǔn)時(shí)差。綜合DSC和MBE兩個(gè)指標(biāo)，本文提出的深度學(xué)習(xí)算法（CC）配準(zhǔn)準(zhǔn)確度最高。

表1 未配準(zhǔn)和傳統(tǒng)算法、深度學(xué)習(xí)算法（MSE、CC、MI）配準(zhǔn)后的DSC和MBE

2.2 T1參數(shù)量化結(jié)果

圖5展示了不同配準(zhǔn)算法運(yùn)動(dòng)校正結(jié)果的實(shí)例。從圖5中放大部分可以觀察到，原始T1參數(shù)量化圖心肌內(nèi)膜附近（圖5中箭頭所示）存在運(yùn)動(dòng)偽影，降低了T1量化的準(zhǔn)確性。傳統(tǒng)算法和深度學(xué)習(xí)算法（CC、MI）在不同程度上減少了運(yùn)動(dòng)偽影，其中深度學(xué)習(xí)（CC）的運(yùn)動(dòng)校正效果最佳。觀察對(duì)應(yīng)的SD量化圖，除深度學(xué)習(xí)（MSE）之外，應(yīng)用不同算法配準(zhǔn)后，心肌部分的SD在整體上均減小了，深度學(xué)習(xí)（CC）的心肌SD整體上最小。

圖5 未配準(zhǔn)和傳統(tǒng)算法、深度學(xué)習(xí)算法（MSE、CC、MI）配準(zhǔn)后的T1參數(shù)量化圖和對(duì)應(yīng)的SD量化圖

心肌部分的T1參數(shù)量化誤差如圖6所示。傳統(tǒng)算法、深度學(xué)習(xí)（CC）、深度學(xué)習(xí)（MI）配準(zhǔn)后心肌部分的平均SD分別為（68.74±31.27）、（59.22±29.26）、（60.35±28.47）ms，均小于原始圖像心肌部分的平均SD（74.37±33.29） ms，差異有統(tǒng)計(jì)學(xué)意義（P<0.05），且深度學(xué)習(xí)（CC）和深度學(xué)習(xí)（MI）的平均SD明顯小于傳統(tǒng)算法（P<0.05），深度學(xué)習(xí)（CC）和深度學(xué)習(xí)（MI）的平均SD無(wú)統(tǒng)計(jì)學(xué)差異（P=0.056）。深度學(xué)習(xí)（MSE）配準(zhǔn)后心肌部分的平均SD增加。

3 討論

本研究提出并驗(yàn)證了一種基于自監(jiān)督深度學(xué)習(xí)的心肌T1參數(shù)量化成像運(yùn)動(dòng)校正算法。該方法突出優(yōu)點(diǎn)為基于深度學(xué)習(xí)的配準(zhǔn)算法，提高效率及準(zhǔn)確度；另一個(gè)優(yōu)勢(shì)為采用了自監(jiān)督學(xué)習(xí)框架，無(wú)需人工標(biāo)注數(shù)據(jù)，大大節(jié)省了醫(yī)生的工作量，為心臟T1量化在臨床上大面積使用提供可能。

圖6 未配準(zhǔn)圖像和傳統(tǒng)算法、深度學(xué)習(xí)算法（MSE、CC、MI）配準(zhǔn)后心肌T1參數(shù)量化的平均SD的均值和標(biāo)準(zhǔn)差

首先，深度學(xué)習(xí)是本方法的優(yōu)勢(shì)之一。深度學(xué)習(xí)適合處理大數(shù)據(jù)問(wèn)題，在圖像識(shí)別、分類、分割、配準(zhǔn)等計(jì)算機(jī)視覺(jué)任務(wù)上取得了許多重大研究成果。區(qū)別于傳統(tǒng)配準(zhǔn)算法依賴于迭代優(yōu)化，將深度學(xué)習(xí)應(yīng)用于T1加權(quán)圖像配準(zhǔn)，可達(dá)到一步配準(zhǔn)，大幅提高計(jì)算效率。經(jīng)過(guò)測(cè)算，本研究提出的方法完成一例配準(zhǔn)平均用時(shí)0.4 s，比傳統(tǒng)算法快20倍。在配準(zhǔn)精度及T1量化準(zhǔn)確度方面，本研究在MOLLI序列采集到的47例受試者的左心室短軸切片數(shù)據(jù)上驗(yàn)證了提出的運(yùn)動(dòng)校正算法，與傳統(tǒng)的基于互相關(guān)的B樣條自由形變模型配準(zhǔn)算法相比，深度學(xué)習(xí)（CC）算法的MBE降低了12%，心肌T1量化SD降低了14%，心肌T1加權(quán)圖像的配準(zhǔn)準(zhǔn)確度和T1參數(shù)量化準(zhǔn)確度都有所提高。另外一個(gè)值得討論的地方是，本研究測(cè)試了不同的損失函數(shù)，以尋求最適宜本問(wèn)題的方案。本研究測(cè)試的三種損失函數(shù)中，CC和MI測(cè)度可以有效校正心肌T1參數(shù)量化圖中的運(yùn)動(dòng)偽影，CC測(cè)度的方法也得到了最佳的效果；然而，選擇MSE作為相似性測(cè)度時(shí)，本文提出的算法無(wú)法提供有效的運(yùn)動(dòng)校正，可能的原因?yàn)檎`差無(wú)法解決信號(hào)強(qiáng)度和對(duì)比度變化的圖像配準(zhǔn)問(wèn)題，而CC和MI在衡量圖像相似程度時(shí)摒除了圖像對(duì)比度變化帶來(lái)的影響。

其次，采用自監(jiān)督學(xué)習(xí)框架是本方法的另一個(gè)優(yōu)勢(shì)。深度學(xué)習(xí)算法的一個(gè)挑戰(zhàn)是訓(xùn)練真值的獲得，與Fahmy等[8]提出的方法需要大量預(yù)先勾畫的心肌輪廓作為訓(xùn)練集不同，本文提出的方法是自監(jiān)督的，網(wǎng)絡(luò)訓(xùn)練過(guò)程無(wú)需真值，實(shí)際應(yīng)用中也不需要預(yù)處理和用戶交互。該方法提出的自監(jiān)督配準(zhǔn)網(wǎng)絡(luò)引入了空間變換層，使得在訓(xùn)練階段不需要運(yùn)動(dòng)變換真值就可進(jìn)行損失函數(shù)的計(jì)算和網(wǎng)絡(luò)優(yōu)化，大大降低了醫(yī)生的工作量，可以將本方法自由應(yīng)用在任何相似的圖像處理應(yīng)用當(dāng)中，如造影劑增強(qiáng)的磁共振心臟、肝臟成像等需要運(yùn)動(dòng)配準(zhǔn)的成像方式中。

在后續(xù)的研究中，為了提升算法在運(yùn)動(dòng)偏移幅度過(guò)大情況下的魯棒性和準(zhǔn)確度，將考慮以下改進(jìn)策略：① 將心肌結(jié)構(gòu)相似性損失加入到目前的損失函數(shù)中[16]；② 在目前非剛性配準(zhǔn)網(wǎng)絡(luò)上增加一個(gè)仿射配準(zhǔn)子網(wǎng)絡(luò)[17]以克服運(yùn)動(dòng)幅度過(guò)大的問(wèn)題。

4 結(jié)論

本文提出了一種基于自監(jiān)督深度學(xué)習(xí)的非剛性配準(zhǔn)算法，為心肌T1參數(shù)量化成像提供了快速、有效的運(yùn)動(dòng)校正，增加了T1參數(shù)量化成像大規(guī)模臨床應(yīng)用的可能性。