基于自一致性的含聯(lián)合全變分的并行磁共振成像高質(zhì)量重構(gòu)方法
【技術(shù)領(lǐng)域】
[0001] 本發(fā)明涉及一種基于自一致性的含聯(lián)合全變分的并行磁共振成像高質(zhì)量重構(gòu)方 法。
【背景技術(shù)】
[0002] 磁共振成像（Magnetic Resonance Imaging, MRI)能夠提供良好的人體軟組織對 比度，并且沒有放射性，因此逐漸成為現(xiàn)代臨床醫(yī)學(xué)中不可缺少的成像工具。然而，受物理 和生理因素的限制，采集磁共振信號的速度很慢。并行成像是一種提高采集速度的常用技 術(shù)。SMASH和SENSE的提出，標志著并行成像成為切實可行的技術(shù)。并行成像使用對磁共振 信號靈敏度不同的多個線圈同時采集磁共振信號。靈敏度信息的使用減少了用于重構(gòu)的數(shù) 據(jù)個數(shù)，從而提高了采集的速度。
[0003] 過去十多年，產(chǎn)生了很多種并行成像重構(gòu)算法。這些算法的不同之處在于：1)顯 式還是隱式使用靈敏度信息；2)恢復(fù)單一圖像還是逐線圈恢復(fù)多幅圖像。SMASH，SENSE， SPACE-RIP，kSPA屬于采用顯式靈敏度信息恢復(fù)單一圖像的方法。PILS，RAPPA，SPIRiT屬 于采用隱式靈敏度信息進行逐線圈恢復(fù)的方法。而AUTO-SMASH采用隱式靈敏度信息恢復(fù) 單一圖像，PARS則采用顯式靈敏度信息進行逐線圈恢復(fù)。
[0004] SENSE是各種算法中應(yīng)用最普遍的，它易于集成圖像的先驗知識。近年來，有研究 人員將壓縮感知與SENSE結(jié)合，增加正則項，取得了較好重構(gòu)效果。當(dāng)靈敏度已知時，SENSE 能獲得最優(yōu)解。然而準確估計線圈的靈敏度往往非常困難。因此，不需顯式使用線圈靈敏 度信息的自校準方法就體現(xiàn)出了優(yōu)勢。
[0005] GRAPPA是一種應(yīng)用廣泛的逐線圈的自校準重構(gòu)方法，由于不需要顯式的使用靈敏 度信息，從而避免了棘手的靈敏度估計。Lustig等在GRAPPA的基礎(chǔ)上提出了一種更為有效 的并行成像技術(shù)的理論框架--SPIRiT。和SENSE類似，SPIRiT將重構(gòu)問題轉(zhuǎn)化成一個逆 問題來求解，易于集成圖像的先驗知識，而且能夠重構(gòu)任意k空間采樣的數(shù)據(jù)。Murphy和 Vasanawala等提出Ll-SPIRiT模型，引入Joint Ll正則項，并使用POCS算法來求解。

【發(fā)明內(nèi)容】

[0006] 本發(fā)明的目的在于克服現(xiàn)有技術(shù)的不足，提供一種基于自一致性的含聯(lián)合全變分 的并行磁共振成像高質(zhì)量重構(gòu)方法。
[0007] 本發(fā)明的目的是通過以下技術(shù)方案來實現(xiàn)的：基于自一致性的含聯(lián)合全變分的并 行磁共振成像高質(zhì)量重構(gòu)方法，其特征在于：它包括以下步驟：
[0008] SO :初始化，令# =0，u〇= D Ty，z〇= 0, t 1= 1，j = 1 ;
[0009] 式中，x為全部線圈的頻域數(shù)據(jù)，r表示空間位置索引；《.=聲1，，^表示逐線圈 的傅里葉變換，々=/_)/.?+ D'>'，D和D。分別表示選擇原頻率采樣點和未采樣點，D τ表示選 擇原采樣頻率點并將元采樣頻率點放回頻率域的原位置，巧表示選擇未被采樣的頻率點 并將未被采樣的頻率點放回頻率域的原位置，$表示未被采樣的頻率點，y表示采集到的頻 率點，Z表示中間變量，t表示與算法加速相關(guān)的中間變量，j表示循環(huán)變量；
[0010] Sl :計算未采樣的重構(gòu)數(shù)據(jù)Xg，計算公式如下：
[0012] 式中，Xg表示未采樣的重構(gòu)數(shù)據(jù)，，b = -(G_I)DTy，G為頻域內(nèi)插操 作子；L為
I的梯度的Lipschitz常數(shù)；
[0013] S2 :將未采樣的重構(gòu)數(shù)據(jù)Xg轉(zhuǎn)換成步驟S3需要的含噪圖像數(shù)據(jù)V，計算公式如 下：
[0016] S3 :計算出去噪后的圖像數(shù)據(jù)包括以下子步驟：
[0017] S30:初始化，令
[0018] 式中，.τ, 表示多線圈圖像變量，N = mXn，C為多線圈個數(shù)，m和η分別為單 個線圈二維圖像的行數(shù)和列數(shù)；d = DhvX1, Z1= Ψχ ρ d和Z1為任意變量；b d和氣,.為對偶變 量；D4,,=(代,A' Γ e K_'VxW，I), = A ? /," e Rjvxjv ' Dv = 4 €) Z)," e ITxjv，其中隊和 D "分別 表示nXn和mXm的循環(huán)矩陣，所述的循環(huán)矩陣結(jié)構(gòu)如下：
[0020] S31 :令k的值為0, k表示循環(huán)變量；
[0021] S32 :計算
[0022] 式中，β 1表示懲罰參數(shù)，shrink2J〇表示聯(lián)合二維收縮算子，計算公式如下：
CN 10bl847bb A ^ "n \J ：V1U 貝
[0025] S32 :計算
[0026] 式中，β 2表示懲罰參數(shù)，shrinkjO表示聯(lián)合一維收縮算子，計算公式如下：
[0028] 式中：
Ψ表示逐線圈的小波變換；
[0029] S33 :令循環(huán)變量c (只在步驟S33~S35中應(yīng)用）的值為1 ;
[0030] S34 :計算
[0031] 式中，a JP α 2為進行rescale處理的正則化參數(shù)；
[0032] S35 :判斷循環(huán)變量c的值是否大于C :若不大于則對循環(huán)變量c做加一操作之后 返回步驟S34,否則進入步驟S36 ;
[0033] S36 :計算
[0034] S37 :計算
[0035] S38 :判斷k值的大小是否大于K :若不大于K則對k做加一操作之后返回步驟 S32,否則進入步驟S39 ;K為循環(huán)次數(shù)；
[0036] S4 :將步驟S34中得到的所有線圈的去噪后的圖像數(shù)據(jù)的集合u]轉(zhuǎn)換成步驟S5的 迭代算法需要的未采樣重構(gòu)數(shù)據(jù)ζ]，計算公式如下：
[0038] S5 :采用來自于FISTA的迭代算法進行差值加速，包括以下子步驟：
[0039] S51 :更新t·"1，計算公式如下：
[0043] S6 :判斷是否滿足條件，若滿足條件則進入步驟S7,否則返回步驟Sl ;所述的條 件為迭代達到最大迭代次數(shù)或者當(dāng)重構(gòu)圖像與前次迭代重構(gòu)出的圖像的相對誤差小于某 值；
[0044] S7 :計算最終重構(gòu)的多線圈頻域數(shù)據(jù)，計算公式如下：
[0045] CN 10bl847bb A ^ "n \J VlU 貝
[0046] S8 :對步驟S7得到的各個線圈的頻域數(shù)據(jù)進行傅里葉反變換，采用SRSOS對各個 線圈圖像進行聯(lián)合，得到最終的單幅重構(gòu)圖像，計算公式如下：
[0048] 本發(fā)明的有益效果是：本發(fā)明基于SPIRiT框架，針對含有JTV和JLl復(fù)合正則項 的并行成像的重構(gòu)問題，提出了一種高質(zhì)量的重構(gòu)算法。首先將有約束的重構(gòu)問題，轉(zhuǎn)化為 無約束的最優(yōu)化問題，再對數(shù)據(jù)保真項和自校正項進行簡化，再采用算子分裂技術(shù)將簡化 后的重構(gòu)問題轉(zhuǎn)化成一個梯度計算問題和一個含有JTV和幾1復(fù)合正則項的去噪問題，復(fù) 合正則項的去噪問題通過全新設(shè)計的基于Split Bregman技術(shù)的算法求解。最后再通過 FISTA進行加速。本發(fā)明設(shè)計實驗，比較了新算法與其他常用算法的重構(gòu)性能。實驗仿真表 明，新算法的收斂速度與POCS算法相當(dāng)，而重構(gòu)圖像的SNR有較大提升。
【附圖說明】
[0049] 圖1為本發(fā)明方法流程圖；
[0050] 圖2為本發(fā)明步驟S3流程圖；
[0051] 圖3為Tl加權(quán)的使用SPGR序列獲取的全采樣的高分辨率大腦成像數(shù)據(jù)（即 datal)圖；
[0052] 圖4為使用2維旋轉(zhuǎn)回城序列獲取的全采樣的軸向腦成像數(shù)據(jù)（即data2)圖；
[0053] 圖5為加速比約為6的泊松亞采樣示意圖；
[0054] 圖6為亞采樣率為6的采用datal測試序列時本發(fā)明與POCS性能比較；
[0055] 圖7為亞采樣率為6的采用data2測試序列時本發(fā)明與POCS性能比較；
[0056] 圖8為測試序列為datal時POCS算法重構(gòu)圖像與原始圖像的差值圖像；
[0057] 圖9為測試序列為datal時本發(fā)明方法重構(gòu)圖像與原始圖像的差值圖像；
[0058] 圖10為測試序列為data2時POCS算法重構(gòu)圖像與原始圖像的差值圖像；
[0059] 圖11為測試序列為data2時本發(fā)明方法重構(gòu)圖像與原始圖像的差值圖像。
【具體實施方式】
[0060] 下面結(jié)合附圖進一步詳細描述本發(fā)明的技術(shù)方案：
[0061] 本發(fā)明是基于SPITiT框架提出的一種高效的重構(gòu)方法。
[0062] SPIRiT是一個基于自校準的并行成像重構(gòu)方法，該方法逐線圈內(nèi)插出亞采樣時缺 失的頻率點，再將多線圈圖像合并成一幅圖像。在SPIRiT中，一個內(nèi)插核gl]通過對頻域中 心全采樣的數(shù)據(jù)（通常稱為自校準信號）進行校準得到。若以X1表示第i個線圈的整個 頻域數(shù)據(jù)，則基于校準的一致性準則可寫成：
[0064] 式中，N。表示線圈數(shù)，表示卷積操作。卷積核g u被稱作SPIRiT核。
[0065] 全部線圈的一致性準則可簡化成矩陣形式：
[0066] X = Gx ； (2)
[0067] 式中，X為全部線圈的頻域數(shù)據(jù)。G為頻域內(nèi)插操作子，使用整個頻域數(shù)據(jù)與gl]相 卷積來完成內(nèi)插。
[0068] 當(dāng)然，除了基于校準的一致性，重構(gòu)還需要與原采樣的數(shù)據(jù)保持一致性，表示成矩 陣形式：
[0069] y = Dx ； (3)
[0070] 式中，D為選擇操作子，從整個頻域中選出采樣的頻率點，y為采集到的頻率點。
[0071] 由于不同線圈的圖像在小波域具有相似性，因此引入聯(lián)合稀疏性模型（Group Sparse)：
[0073] 式中，c為線圈索引，r為空間位置索引。
[0074] 根據(jù)公式（2)和公式（3)的兩個限制，結(jié)合聯(lián)合稀疏性模型公式（4)，并行成像的 重構(gòu)問題可轉(zhuǎn)化成求解帶聯(lián)合稀疏性正則項（幾1)的優(yōu)化問題：
[0076] 式中，其中Ψ為逐線圈的小波變換，.F為逐線圈的傅里葉變換。Lustig等采用 Projection Over Convex Sets(POCS)來求解公式（5)〇
[0077] 另外，Murphy和Vasanawala使用圖像域校正，將重構(gòu)問題寫成：
[0079] 式中，G為圖像域內(nèi)插操作子。Murphy和Vasanawala都是用POCS從圖像域直接 重構(gòu)各個線圈圖像。頻域和圖像域校正的POCS算法在本質(zhì)上類似，因此具有非常相似的重 構(gòu)性能。
[0080] 在考慮到多線圈圖像梯度的相似性，我們在重構(gòu)時引入Joint Total Variation(JTV)正則項。引入JTV之后，基于頻域校正的并行成像的重構(gòu)問題可轉(zhuǎn)化為以 下的優(yōu)化問題：
[0082] 式中，Ψ表示逐線圈的小波變換，，' 表示逐線圈的傅里葉變換，G為頻域內(nèi)插操 作子。
[0083] 利用二次罰函數(shù)技術(shù)，將公式（5)變換為：
[0085] 在很多情況下，保持已采樣的頻率系數(shù)不變是很有必要的。為此，我們提出了一種 簡化公式（8)的方法。
[0086] 假設(shè)i表示未被采樣的頻率點，D和D。分別表示選擇原頻率采樣點和未采樣點，D τ 表示選擇原采樣頻率點并將其放回頻率域的原位置，D/表示選擇未被采樣的頻率點并將其 放回頻率域的原位置，則有X可以表示成j+D/.h因此，式⑶可以簡化成：
[0088] 式中，.才=(?-b = _(G-I)DTy，