基于交替方向乘子法(ADMM)直接重建心臟11C-acetate PET動(dòng)力學(xué)參數(shù)圖

何升級(jí)，劉 帥，張 輝，霍 力，尚 斐*

(1.北京理工大學(xué)生命學(xué)院生物醫(yī)學(xué)工程系，北京 100081；2.清華大學(xué)醫(yī)學(xué)院生物醫(yī)學(xué)工程系，北京 100084； 3.中國(guó)醫(yī)學(xué)科學(xué)院 北京協(xié)和醫(yī)學(xué)院 北京協(xié)和醫(yī)院核醫(yī)學(xué)科，北京 100730)

11C-acetate是臨床非侵入式定量評(píng)估心肌代謝和血流灌注的重要示蹤劑之一[1-2]。目前多采用單組織房室模型(即單室模型)計(jì)算心臟11C-acetate PET動(dòng)態(tài)成像動(dòng)力學(xué)參數(shù)[3-4]，一般分為兩步：①基于PET投影數(shù)據(jù)重建PET動(dòng)態(tài)圖像；②采用動(dòng)力學(xué)模型逐像素?cái)M合時(shí)間-活度曲線，得到動(dòng)力學(xué)參數(shù)圖。為提高動(dòng)力學(xué)參數(shù)圖像質(zhì)量，有學(xué)者[5]將圖像重建與動(dòng)力學(xué)模型相結(jié)合而提出直接重建算法，即根據(jù)投影數(shù)據(jù)直接重建動(dòng)力學(xué)參數(shù)圖，以降低噪聲對(duì)參數(shù)擬合的影響。現(xiàn)有相關(guān)研究[6-8]多采用圖模型(如Patlak和Logan模型)重建動(dòng)力學(xué)參數(shù)圖。本研究觀察基于交替方向乘子法(alternating direction method of multipliers, ADMM)直接重建心臟11C-acetate PET動(dòng)力學(xué)參數(shù)圖的應(yīng)用價(jià)值。

1 資料與方法

1.1 數(shù)據(jù)來(lái)源 PET/CT數(shù)據(jù)來(lái)自1例有30年飲酒史、無(wú)其他心臟疾病，于2017年7月在北京協(xié)和醫(yī)院接受PET/CT檢查的51歲男性早期酒精性心肌病患者。本研究獲醫(yī)院倫理委員會(huì)批準(zhǔn)(倫理審批號(hào)：20160012)，檢查前患者簽署知情同意書。

1.2 儀器與方法 采用PoleStar m660 PET/CT掃描儀。囑患者仰臥，先行低劑量肺部CT掃描，管電壓120 kV，管電流140 mA；之后經(jīng)靜脈注射11C-acetate 740 MBq，行40 min動(dòng)態(tài)PET心臟掃描；基于CT圖像進(jìn)行衰減校正后，采用有序子集最大期望值法(ordered subsets expectation maximization, OSEM)及飛行時(shí)間技術(shù)重建PET圖像，圖像大小192×192×117，分辨率3.147 mm×3.147 mm×1.866 mm，共得到53幀重建圖像，包括15幀10 s、15幀30 s、16幀60 s及7幀120 s圖像。

1.3 構(gòu)建模型

1.3.1 PET動(dòng)態(tài)圖像重建模型 將第n幀PET圖像記為xn∈N×1(N表示圖像的像素?cái)?shù)目，表示實(shí)數(shù)域)，將測(cè)得的第n幀投影數(shù)據(jù)記為yn∈M×1，其期望為M×1(M表示響應(yīng)線數(shù)目)；則與xn的關(guān)系為：

(1)

式中，E為期望符號(hào)，P∈M×N為系統(tǒng)矩陣，rn∈M×1表示第n幀隨機(jī)和散射事件的期望；則動(dòng)態(tài)圖像重建似然方程L為：

(2)

式中，y為動(dòng)態(tài)投影數(shù)據(jù)，y=[y1,y2,…,yT]；x為動(dòng)態(tài)圖像，x=[x1,x2,…,xT]；T為圖像重建幀數(shù)；i表示投影數(shù)據(jù)射線的索引。

1.3.2 心臟11C-acetate單室模型 見圖1。

圖1 心臟11C-acetate單室模型

其中，CP(t)和CT(t)分別表示血液和心肌中的11C-acetate濃度隨時(shí)間t的變化值；K1和k2表示11C-acetate在血液和心肌之間交換的速率常數(shù)，是臨床量化評(píng)估心肌的重要指標(biāo)，K1與心肌血流量相關(guān)，k2與心肌耗氧量相關(guān)[1,4]。該模型用微分方程表示為：

(3)

對(duì)公式(3)求解，可得：

CT(t)=K1CP(t)?e-k2t

(4)

式中，?代表卷積運(yùn)算；由于PET成像存在容積效應(yīng)，需要對(duì)心肌區(qū)域進(jìn)行溢出及恢復(fù)校正[9]，校正公式為：

CS(t)=(1-vl)×CT(t)+vl×CP(t)

(5)

式中，CS(t)為PET圖像中實(shí)際測(cè)得的心肌區(qū)域示蹤劑濃度，vl∈[0,1]為血液容積的比例因子。

圖2 直接重建算法流程圖

對(duì)公式(4)和(5)進(jìn)行擬合，可得到參數(shù)K1、k2和vl，參數(shù)集合為θ=[K1,k2,vl]。

1.3.3 直接重建流程 將重建PET動(dòng)態(tài)圖像和單室模型參數(shù)求解兩個(gè)獨(dú)立步驟合二為一，以實(shí)現(xiàn)直接重建參數(shù)圖像。動(dòng)力學(xué)模型中的參數(shù)θ與動(dòng)態(tài)圖像x關(guān)系為x=K(θ)，其中K(θ)表示上述心臟11C-acetate單室模型；結(jié)合公式(1)可知，參數(shù)θ與投影數(shù)據(jù)的關(guān)系為：

(6)

(7)

對(duì)θ則可通過如下公式進(jìn)行優(yōu)化并求解：

(8)

(9)

再通過增廣拉格朗日法將公式(9)轉(zhuǎn)化為無(wú)約束優(yōu)化問題，對(duì)應(yīng)優(yōu)化目標(biāo)函數(shù)為：

(10)

最后通過ADMM框架進(jìn)行迭代求解：

(11)

(12)

μn+1=μn+xn+1-K(θn+1)

(13)

式中，ρ為懲罰因子，μ為對(duì)偶變量，n表示迭代次數(shù)。求解流程見圖2。

1.4 動(dòng)力學(xué)參數(shù)重建

1.4.1 參考標(biāo)準(zhǔn) 對(duì)本例PET圖像所示心臟區(qū)域進(jìn)行裁剪，裁剪后圖像大小為64×64×48；再由1名具有10年工作經(jīng)驗(yàn)的放射科主治醫(yī)師勾畫左心室心外膜和心內(nèi)膜，獲得左心室血池和心肌區(qū)域，并基于左心室血池獲取CP(t)，根據(jù)公式(2)和(3)計(jì)算得到K1、k2和vl參數(shù)圖；最后根據(jù)美國(guó)心臟協(xié)會(huì)提出的17節(jié)段模型[10]將左心室心肌區(qū)域參數(shù)圖映射生成參數(shù)圖，并以之作為標(biāo)準(zhǔn)，間接重建采用OSEM方法重建動(dòng)態(tài)圖像。

1.4.2 直接重建 采用MATLAB R2017a軟件，操作系統(tǒng)為64位Windows 10，CPU為英特爾i7-9700F，16 GB內(nèi)存?；趩问夷Ｐ陀蓞?shù)圖獲得53幀理想動(dòng)態(tài)圖像，并對(duì)其進(jìn)行120個(gè)角度(0～180°平均劃分)投影后，加入泊松噪聲；采用開源NiftyRec工具箱[11]實(shí)現(xiàn)基于ADMM的直接重建算法，總迭代次數(shù)15，懲罰項(xiàng)ρ為2.0×10-7，子迭代次數(shù)2，初始μ0為0；并通過濾波反投影法重建得到初始圖像x0，通過11C-acetate單室模型對(duì)x0進(jìn)行逐像素曲線擬合，得到初始參數(shù)圖θ0。見圖3。

1.5 圖像評(píng)估 通過偏差(bias)評(píng)估圖像重建結(jié)果，其定義式為：

(14)

式中，Iest為通過直接或間接重建方法計(jì)算得到的值，Igs為參考標(biāo)準(zhǔn)。

1.6 統(tǒng)計(jì)學(xué)分析 采用線性回歸分析參考標(biāo)準(zhǔn)與直接、間接重建圖像的相關(guān)性，并計(jì)算Pearson相關(guān)系數(shù)；以0<│r│<0.3為弱相關(guān)，0.3≤│r│<0.5為低度相關(guān)，0.5≤│r│<0.8為中度相關(guān)，│r│≥0.8為高度相關(guān)。P<0.05為差異有統(tǒng)計(jì)學(xué)意義。

2 結(jié)果

間接重建耗法耗時(shí)23 min，直接重建法耗時(shí)140 min。

圖3 心臟11C-acetate PET動(dòng)力學(xué)參數(shù)圖直接重建流程圖

圖4 參考標(biāo)準(zhǔn)、間接重建及直接重建參數(shù)圖及其偏差

參考標(biāo)準(zhǔn)與間接重建得到的K1、k2和vl值的平均偏差分別為2.97%、1.61%和10.19%，與通過直接重建得到的K1、k2和vl值的平均偏差分別為1.16%、1.06%和4.30%。見圖4。

參考標(biāo)準(zhǔn)與通過間接重建得到的K1(r=0.94,P<0.001)、k2(r=0.94,P<0.001)和vl值(r=0.91,P<0.001)及與通過直接重建得到的K1(r=0.99,P<0.001)、k2(r=0.97,P<0.001)和vl值(r=0.99,P<0.001)均高度相關(guān)，且參考標(biāo)準(zhǔn)與直接重建法的相關(guān)系數(shù)略大。見圖5。

3 討論

心臟11C-acetate PET動(dòng)態(tài)成像有助于評(píng)估心肌代謝和血流灌注，目前多采用圖模型(如Patlak和Logan模型)重建其動(dòng)力學(xué)參數(shù)圖[6-8]。GONG等[7]結(jié)合18F-FDG PET數(shù)據(jù)及MRI所示結(jié)構(gòu)信息，以基于ADMM的非監(jiān)督深度學(xué)習(xí)法直接重建，得到腦部Patlak動(dòng)力學(xué)參數(shù)圖。GALLEZOT等[8]對(duì)傳統(tǒng)Logan模型進(jìn)行改良，采用最大期望值法(expectation maximization, EM)直接重建，得到腦部11C-PBR28動(dòng)力學(xué)參數(shù)圖。上述研究結(jié)果均表明，通過圖模型直接重建法得到的參數(shù)圖質(zhì)量較高。

圖5 參考標(biāo)準(zhǔn)與間接重建K1(A)、k2(B)和vl(C)參數(shù)圖及與直接重建K1(D)、k2(E)和vl(F)參數(shù)圖的相關(guān)性分析圖

SHI等[12]發(fā)現(xiàn)，計(jì)算11C-acetate動(dòng)力學(xué)參數(shù)時(shí)，基于OSEM重建圖像的穩(wěn)定性和魯棒性更好。本研究采用OSEM間接重建心臟11C-acetate PET動(dòng)力學(xué)參數(shù)圖，并基于ADMM實(shí)現(xiàn)直接重建，結(jié)果表明直接及間接重建得到的K1、k2和vl均與參考標(biāo)準(zhǔn)高度相關(guān)；直接重建法與參考標(biāo)準(zhǔn)的相關(guān)系數(shù)略大于間接重建，且直接重建得到的K1、k2和vl與參考標(biāo)準(zhǔn)間的偏差更小。分析原因，主要在于間接重建過程中動(dòng)態(tài)圖像易受噪聲干擾而在擬合時(shí)出現(xiàn)奇異值，降低計(jì)算動(dòng)力學(xué)參數(shù)的準(zhǔn)確性，同時(shí)由于參數(shù)擬合時(shí)忽略了動(dòng)態(tài)圖像像素間的相關(guān)性，導(dǎo)致所獲參數(shù)圖常具有較大噪聲；而在直接重建迭代過程中，每一次迭代均包含圖像重建和動(dòng)力學(xué)參數(shù)擬合兩個(gè)步驟，且優(yōu)化過程中的θ與x相互約束，加之重建圖像時(shí)考慮了空間相關(guān)性，使得參數(shù)圖的擬合誤差較小，抗噪性能則有所提高。

本研究間接重建法耗時(shí)23 min，直接重建法耗時(shí)140 min，主要原因在于采用單室模型進(jìn)行直接重建，其動(dòng)力學(xué)參數(shù)與PET數(shù)據(jù)的關(guān)系為非線性，使算法較為復(fù)雜，增加了時(shí)間成本；但非線性單室模型是以示蹤劑在人體內(nèi)的代謝過程為基礎(chǔ)而建立的動(dòng)力學(xué)數(shù)學(xué)模型，與人體內(nèi)生化過程緊密相關(guān)，有助于得到更具生理意義指標(biāo)，以定量評(píng)估靶器官功能。

目前已有學(xué)者嘗試將正則化方法，如全變分正則化[13]、核方法[14]等，加入直接重建框架中，以提高重建圖像質(zhì)量，且已知采用深度神經(jīng)網(wǎng)絡(luò)重建PET圖像可提高成像質(zhì)量[15]。后續(xù)研究可將深度神經(jīng)網(wǎng)絡(luò)嵌入直接重建PET參數(shù)圖中，以進(jìn)一步提高其動(dòng)力學(xué)參數(shù)圖的圖像質(zhì)量。

綜上所述，相比間接重建方法，基于ADMM直接重建心臟11C-acetate PET動(dòng)力學(xué)參數(shù)圖的準(zhǔn)確性更高。