本發(fā)明涉及磁共振成像技術(shù)領(lǐng)域，具體涉及一種壓縮感知理論下的基于高階全變分正則化模型的增廣拉格朗日快速迭代磁共振圖像重建方法。

背景技術(shù)：

磁共振成像因其非創(chuàng)傷性、無電離輻射等優(yōu)勢，在醫(yī)學(xué)檢測中發(fā)揮著重要作用。然而，過長的掃描時間等問題阻礙著磁共振成像的發(fā)展。壓縮感知理論利用圖像的稀疏性，在采集部分k空間樣本的情況下，通過非線性算法重建磁共振圖像，實現(xiàn)加速磁共振成像這一目標。

在基于壓縮感知理論的磁共振圖像重建算法中，常用方法是正則化方法，即將圖像重建問題歸納為一個求解代價方程的問題，代價方程中包含數(shù)據(jù)一致性項和對所求解問題進行約束的正則項；再通過求解此代價方程獲得重建圖像。

目前使用廣泛的一種正則化方法為全變分方法，這一方法因其能在平滑圖像的同時有效地保留圖像的紋理和邊緣等細節(jié)信息而被廣泛應(yīng)用。然而，這一方法仍然存在著許多問題和局限性，例如全變分正則化比較適合具有分塊常量特性的圖像處理問題，而在重建具有較多細節(jié)的磁共振圖像時常會產(chǎn)生片狀偽影等。這主要是由于全變分方法利用圖像的一階導(dǎo)數(shù)，因此對于平滑區(qū)域、細節(jié)區(qū)域等重建效果不佳。為了改善這一狀況，研究人員提出了許多擴展方法，相比于全變分方法來說，這些方法雖然重建效果有所提高，但大多以增加計算復(fù)雜度為代價，降低了計算效率。因此，設(shè)計一種計算效率高、耗時短、且能克服全變分產(chǎn)生的偽影問題的磁共振圖像重建方法是目前該領(lǐng)域中的一個重要目標。

技術(shù)實現(xiàn)要素：

本發(fā)明的目的是針對目前磁共振圖像重建方法的不足，提供一種可以同時提高重建圖像質(zhì)量及計算效率的基于高階全變分正則化的快速迭代磁共振圖像重建方法。

本發(fā)明為解決上述技術(shù)問題采取的技術(shù)方案是：

(1)利用預(yù)先設(shè)置的欠采樣模板獲取磁共振圖像的部分k空間數(shù)據(jù)；

(2)利用獲取的部分k空間數(shù)據(jù)和高階全變分正則化法建立磁共振圖像重建模型；

(3)對部分k空間數(shù)據(jù)直接進行傅里葉逆變換，得到空間域預(yù)估磁共振圖像作為初始重建圖像；

(4)引入輔助變量和增廣拉格朗日乘子以便重建模型的快速迭代求解；

(5)計算初始重建圖像中每個像素的高階方向?qū)?shù)，利用軟閾值收縮方法更新高階方向?qū)?shù)，然后用更新的高階方向?qū)?shù)獲得本次迭代得到的磁共振重建圖像；

(6)判斷當前重建圖像結(jié)果是否滿足收斂條件，若是則獲得最終重建磁共振圖像，否則進入步驟(7)；

(7)增加迭代參數(shù)取值，以當前迭代步驟中更新的磁共振圖像為初始重建圖像，返回步驟(5)繼續(xù)進行循環(huán)迭代操作。

上述步驟(2)中重建模型的建立，具體采用如下方式進行：

若將D_u，n定義為n階微分算子，f為目標重建圖像，u∈S表示單位圓S上任一方向的向量，則圖像在某一方向上的n階導(dǎo)數(shù)可表示為如下向量形式：

其中，D_u，nf表示圖像在單位向量u＝(u_x，u_y)方向上的n階方向?qū)?shù)，s_n(u)為n階導(dǎo)數(shù)下的單位向量分量，為圖像的n階偏導(dǎo)。實驗發(fā)現(xiàn)，方向?qū)?shù)階數(shù)過高會導(dǎo)致過度平滑現(xiàn)象，因此本發(fā)明所使用的主要為二階方向?qū)?shù)，即：

其中，

因此，可以根據(jù)n階方向?qū)?shù)，建立高階全變分正則項即基于高階方向?qū)?shù)的圖像L₁-L₁范數(shù)：

利用高階全變分正則項和欠采樣獲取的k空間數(shù)據(jù)，建立基于高階全變分正則化的壓縮感知圖像重建模型：

其中A為磁共振圖像傅里葉變換后k空間欠采樣操作算子，b為獲得的k空間欠采樣數(shù)據(jù)，f為待重建圖像，D_u，nf為圖像的每一像素點在單位圓上任意方向的n階方向?qū)?shù)，表示L₁范數(shù)，λ為用于平衡數(shù)據(jù)一致性項||Af-b||²和正則項的正則化參數(shù)，C(f)為需要求解的代價方程。式(5)的目的即為通過最小化代價方程C(f)，獲得重建出的磁共振結(jié)果圖像

上述步驟(4)具體包括：

由于圖像重建模型優(yōu)化問題(5)中D_u，nf的L₁范數(shù)不是處處可微，為了便于求解，將L₁范數(shù)用Huber方程進行估計，Huber方程的形式如下：

利用Huber方程估計后的代價方程形式變?yōu)椋?/p>

上式中，β為迭代參數(shù)，在每次迭代中，逐漸增加β的取值，當β→∞時，C_β(f)＝C(f)。

上述步驟(4)具體包括：

對于式(7)的求解，我們引入輔助變量z_u和可進一步加快求解過程收斂速度的拉格朗日乘子Λ_u，得到的C_β(f)的逼近形式如下：

至此，原始圖像重建優(yōu)化問題(5)便轉(zhuǎn)換成為一個較易求解的優(yōu)化問題，即：

上述步驟(5)具體包括：

分別固定方程中的z_u和D_u，nf，可利用交替最小化求解方程，獲得重建圖像。

首先，假定D_u，nf為固定常數(shù)，以z_u為未知數(shù)求解(9)方程，具體可表示為：

利用軟閾值方法求解上式可得：

其次，假定z_u為固定常數(shù)，以D_u，nf為未知數(shù)求解(9)方程，忽略式中與z_u相關(guān)項可得如下方程：

最小化上式可通過直接求一階導(dǎo)數(shù)來得到結(jié)果，因此我們可以得出關(guān)于f的二次方程(12)最小時滿足下式：

利用導(dǎo)數(shù)算子的旋轉(zhuǎn)可操控性質(zhì)，可以得到和的簡單表示形式：

其中，q_n＝∫_Ss_n(u)z_udu，r_n＝∫_Ss_n(u)Λ_udu，E_n為n階偏導(dǎo)算子，即本發(fā)明中的二階偏導(dǎo)算子即為E₂f＝[f_xx，f_xy，f_yy]^T；

在求解z_u的過程中，為了能夠獲得較好的圖像重建結(jié)果，需在所有的方向上對z_u進行考慮，這使得當z_u在某一方向上數(shù)量級比較大時所需的計算存儲容量會變得很大。但在f的求解中可發(fā)現(xiàn)，實際只需要計算z_u在s(u)上的投影q_n(u)＝∫_Ss(u)z_udu，這樣的簡化便使得計算所需的存儲代價減少了很多，因此可以大幅提高計算效率。

根據(jù)上述推導(dǎo)，式(13)等價于：

在離散傅里葉域，二階算子能夠簡單地表達為如下形式：

其中，為離散傅里葉變換算子。當對磁共振圖像的欠采樣在笛卡爾坐標系下進行時，可將A算子表示為其中S為頻域采樣算子，表示對圖像的傅里葉變換和欠采樣兩個步驟。因此，式(17)可以寫成：

由此，可以得到f的解析表達式，即本次迭代中的重建圖像結(jié)果：

當采樣不在非笛卡爾坐標系下進行，也可簡單地利用共軛梯度法進行最小化計算，求解(17)。

在每一次迭代計算后需要將引入的拉格朗日乘子Λ_u進行如下更新：

其中為本次迭代所使用的拉格朗日乘子，為用于下一次迭代中的拉格朗日乘子，D_u，nf和z_u分別為當前迭代所使用初始圖像的高階方向?qū)?shù)和z_u的解。

上述步驟(6)具體包括：

判斷當前重建圖像是否滿足收斂條件的方法為計算相鄰兩次重建圖像的L₂范數(shù)誤差，即當誤差小于預(yù)定閾值時，停止迭代，獲得最終重建磁共振圖像，否則繼續(xù)循環(huán)迭代。

上述步驟(7)具體包括：

每次迭代中增加迭代參數(shù)β的取值：

β＝β₀·β_inc (22)

其中β₀為本次迭代中所使用的參數(shù)值，β_inc為每次循環(huán)所增加的倍數(shù)。

本發(fā)明的有益效果是：

本發(fā)明所提出的算法計算效率較高，CPU運行時間較短，是一種高效的算法；同時利用本發(fā)明提出的算法，可以較好地重建磁共振圖像，較完整地保留圖像中的細節(jié)區(qū)域。實驗表明，與全變分方法、圖像高階導(dǎo)數(shù)Laplacian方法、小波方法等方法相比，本發(fā)明能夠在較高加速系數(shù)下獲得較高質(zhì)量的重建圖像，同時提高了重建速度，從而實現(xiàn)了加速磁共振圖像重建的目的

附圖說明

圖1為本發(fā)明方法流程圖(高階全變分正則化磁共振圖像重建算法流程圖)。

圖2為仿真實驗所用參考磁共振圖像；圖中：(a)為側(cè)向腦部圖像，(b)為軸向腦部圖像，(c)為血管造影圖像，(d)為腕部圖像。

圖3為利用不同方法對圖2-(a)側(cè)向腦部磁共振圖像在不同欠采樣參數(shù)和噪聲水平下重建圖像對比(磁共振圖像側(cè)向腦部圖像重建結(jié)果)；圖中：

(a)為側(cè)向腦部圖像放大圖(原磁共振圖像)，(b)為三階變分重建結(jié)果圖(4.35倍欠采樣，添加40dB噪聲；信噪比SNR＝21.01dB)，(c)為曲波變換重建結(jié)果圖(4.35倍欠采樣，添加40dB噪聲；SNR＝17.25dB)，(d)為拉普拉斯重建結(jié)果圖(4.35倍欠采樣，添加40dB噪聲；SNR＝16.19dB)，

(e)為全變分重建結(jié)果圖(SNR＝23.55dB；4.35倍欠采樣，添加40dB噪聲)，(f)為二階變分重建結(jié)果圖(SNR＝24.35dB；4.35倍欠采樣，添加40dB噪聲)，(g)為各項同性二階變分圖(SNR＝16.19dB；4.35倍欠采樣，添加40dB噪聲)，(h)為Lysaker重建結(jié)果圖(SNR＝16.19dB；4.35倍欠采樣，添加40dB噪聲)，

(i)為全變分重建結(jié)果圖(SNR＝28.40dB；2倍欠采樣，添加40dB噪聲)，(j)為二階變分重建結(jié)果圖(SNR＝30.15dB；2倍欠采樣，添加40dB噪聲)，(k)為各項同性二階變分圖(SNR＝29.90dB；2倍欠采樣，添加40dB噪聲)，(l)為Lysaker重建結(jié)果圖(SNR＝26.89dB；2倍欠采樣，添加40dB噪聲)。

圖4為利用不同方法對圖2-(d)腕部磁共振圖像在不同欠采樣參數(shù)和噪聲水平下重建圖像對比(欠采樣和添加噪聲下的腕部圖像重建結(jié)果圖)；圖中；

(a)為原圖像，(b)為全變分正則化重建結(jié)果圖(3倍欠采樣，添加噪聲35dB，重建結(jié)果SNR＝25.34dB)，(c)全變分正則化誤差圖(3倍欠采樣，添加噪聲35dB，重建結(jié)SNR＝25.34dB)，

(d)為原圖像放大圖，(e)為二階變分正則化結(jié)果圖(3倍欠采樣，添加噪聲35dB，重建結(jié)果SNR＝25.65dB)，(f)為二階變分正則化誤差圖(3倍欠采樣，添加噪聲35dB，重建結(jié)果SNR＝25.65dB)。

圖5為利用不同方法對圖2-(c)血管造影磁共振圖像在不同欠采樣參數(shù)和噪聲水平下重建圖像對比(欠采樣和添加噪聲下的血管造影圖像重建結(jié)果)；圖中：

(a)為原圖像，(b)為全變分方法重建結(jié)果圖(2.86倍欠采樣，添加噪聲30dB，重建結(jié)果SNR＝29.24dB)，(c)為全變分方法誤差圖(2.86倍欠采樣，添加噪聲30dB，重建結(jié)果SNR＝29.24dB)，

(d)為原圖像放大圖，(e)為二階全變分方法重建結(jié)果圖(2.86倍欠采樣，添加噪聲30dB，重建結(jié)果SNR＝31.30dB)，(f)為二階全變分方法誤差圖(2.86倍欠采樣，添加噪聲30dB，重建結(jié)果SNR＝31.30dB)。

圖6為本發(fā)明方法的收斂性能及運算效率分析對比圖；圖中：(a)為不同β值所對應(yīng)的收斂趨勢曲線圖，(b)為本發(fā)明所用算法的CPU計算時間(圖中，短的淺色曲線為對尺寸為256×256的圖像進行重建；長的深色曲線為對尺寸為512×512的圖像進行重建)。

具體實施方式

下面結(jié)合附圖和實施例對本發(fā)明進行詳細說明。

如圖1所示，本發(fā)明的具體實施步驟如下：

(1)利用預(yù)先設(shè)置的欠采樣模板獲取部分k空間數(shù)據(jù)，為了驗證本發(fā)明的效果，采用了四組參考磁共振圖像，如圖2所示，分別為側(cè)向腦部磁共振圖像(a)，軸向腦部磁共振圖像(b)，血管造影磁共振圖像(c)，腕部磁共振圖像(d)，對參考圖像進行傅里葉變換，采集原始k空間數(shù)據(jù)，采集到的欠采樣k空間數(shù)據(jù)表示為b＝Af+n，其中A為對磁共振圖像進行傅里葉變換后在k空間進行欠采樣操作算子，n為實際采樣中可能存在的加性噪聲，b為獲得的k空間欠采樣數(shù)據(jù)，f為待重建圖像；

對測量數(shù)據(jù)b直接進行補零傅里葉逆變換得到空間域初始重建圖像其中表示對測量數(shù)據(jù)做傅里葉逆變換；

(2)利用獲取的部分k空間數(shù)據(jù)和高階全變分正則化法建立磁共振圖像重建模型；

(3)對部分k空間數(shù)據(jù)直接進行傅里葉逆變換，得到空間域預(yù)估磁共振圖像作為初始重建圖像；

(4)引入輔助變量和增廣拉格朗日乘子以便重建模型的快速迭代求解；

(5)計算初始重建圖像中每個像素的高階方向?qū)?shù)，利用軟閾值收縮方法更新高階方向?qū)?shù)，然后用更新的高階方向?qū)?shù)獲得本次迭代得到的磁共振重建圖像；

(6)判斷當前重建圖像結(jié)果是否滿足收斂條件，若是則獲得最終重建磁共振圖像，否則進入步驟(7)；

(7)增加迭代參數(shù)取值，以當前迭代步驟中更新的磁共振圖像為初始重建圖像，返回步驟(5)繼續(xù)進行循環(huán)迭代操作。

上述步驟(2)中重建模型的建立，具體采用如下方式進行：

將圖像在某一方向上的n階方向?qū)?shù)表示為向量形式：

其中f為目標圖像，D_u，n為n階微分算子，D_u，nf表示圖像在單位向量u＝(u_x，u_y)方向上的n階方向?qū)?shù)，s_n(u)為n階導(dǎo)數(shù)下的單位向量分量，為n階圖像偏導(dǎo)。實驗發(fā)現(xiàn)方向?qū)?shù)階數(shù)過高會導(dǎo)致過度平滑現(xiàn)象，因此本發(fā)明所使用的主要為二階方向?qū)?shù)，即：

其中，

因此，可以根據(jù)n階方向?qū)?shù)，建立高階全變分正則項即基于高階方向?qū)?shù)的圖像L₁-L₁范數(shù)：

本發(fā)明中所使用的是二階全變分正則項，即：

利用高階全變分正則項和欠采樣獲取的k空間數(shù)據(jù)，建立基于高階全變分正則化的壓縮感知圖像重建模型：

其中A為磁共振圖像傅里葉變換后k空間欠采樣操作算子，b為獲得的k空間欠采樣數(shù)據(jù)，f為待重建圖像，D_u，nf為圖像的每一像素點在單位圓上任意方向的n階方向?qū)?shù)，λ為用于平衡數(shù)據(jù)一致性項||Af-b||²和正則項的正則化參數(shù)，C(f)為需要求解的代價方程。式(6)的目的即為通過最小化代價方程C(f)，獲得重建出的磁共振圖像

上述步驟(4)具體包括：

由于圖像重建模型優(yōu)化問題(6)中L₁范數(shù)方程不是處處可微，為了便于求解，將L₁范數(shù)用Huber方程進行估計，Huber方程的形式如下：

利用Huber方程估計后的代價方程形式變?yōu)椋?/p>

上式中，當β→∞時，C_β(f)＝C(f)。

對于式(8)的求解，我們引入輔助變量z_u和可進一步加快求解過程收斂速度的拉格朗日乘子Λ_u，得到的C_β(f)的逼近形式如下：

至此，原始圖像重建優(yōu)化問題(6)便轉(zhuǎn)換成為一個較易求解的優(yōu)化問題，即：

上述步驟(5)具體包括：

分別固定方程中的z_u和D_u，nf，可利用交替最小化求解方程，獲得重建圖像。

首先，假定D_u，nf為固定常數(shù)，以z_u為未知數(shù)求解(10)方程，具體可表示為：

利用軟閾值方法求解上式可得：

其次，假定z_u為固定常數(shù)，以D_u，nf為未知數(shù)求解(10)方程，忽略式中與z_u相關(guān)項可得如下方程：

最小化上式可通過直接求一階導(dǎo)數(shù)來得到結(jié)果，因此我們可以得出關(guān)于f的二次方程(13)最小時滿足下式：

利用導(dǎo)數(shù)算子的旋轉(zhuǎn)可操控性質(zhì)，可以得到和的簡單表示形式：

其中，q_n＝∫_Ss_n(u)z_udu，r_n＝∫_Ss_n(u)Λ_udu，E_n為n階偏導(dǎo)算子，即本發(fā)明中的二階偏導(dǎo)算子即為E₂f＝[f_xx，f_xy，f_yy]^T；

在求解z_u的過程中，為了能夠獲得較好的圖像重建結(jié)果，需在所有的方向上對z_u進行考慮，這使得當z_u在某一方向上數(shù)量級比較大時所需的計算存儲容量會變得很大。但在f的求解中可發(fā)現(xiàn)，實際只需要計算z_u在s(u)上的投影q_n(u)＝∫_Ss(u)z_udu，這樣的簡化便使得計算所需的存儲代價減少了很多，因此可以大幅提高計算效率。

根據(jù)上述推導(dǎo)，式(14)等價于：

在離散傅里葉域，二階算子能夠簡單地表達為如下形式：

其中，為離散傅里葉變換算子。當對磁共振圖像的欠采樣在笛卡爾坐標系下進行時，可將A算子表示為其中S為頻域采樣算子，表示對圖像的傅里葉變換和欠采樣兩個步驟。因此，式(18)可以寫成：

因此，可以得到f的解析表達式，即本次迭代中的重建圖像結(jié)果：

當采樣不在非笛卡爾坐標系下進行，也可簡單地利用共軛梯度法進行最小化計算，求解(18)。

在每一次迭代計算后需要將引入的拉格朗日乘子Λ_u進行如下更新：

其中為本次迭代所使用的拉格朗日乘子，為用于下一次迭代中的拉格朗日乘子，D_u，nf和z_u分別為當前迭代所使用初始圖像的高階方向?qū)?shù)和z_u的解。

上述步驟(6)具體包括：

判斷當前重建圖像是否滿足收斂條件的方法為計算相鄰兩次重建圖像的L₂范數(shù)誤差，即當誤差小于預(yù)定閾值時，停止迭代，獲得最終重建磁共振圖像，否則繼續(xù)循環(huán)迭代。

上述步驟(7)具體包括：

每次迭代中增加迭代參數(shù)β的取值：

β＝β₀·β_inc (23)

其中β₀為本次迭代中所使用的參數(shù)值，β_inc為每次循環(huán)所增加的倍數(shù)。

為了定量分析本發(fā)明中所進行實驗的結(jié)果，本發(fā)明采用SNR和CPU運行時間兩個指標對結(jié)果進行分析：

其中f_ref為參考圖像，為重建出的磁共振結(jié)果圖像。

圖3為利用圖2(a)側(cè)向腦部磁共振圖像，在不同欠采樣參數(shù)下基于不同方法重建出的結(jié)果圖像對比。其中，(b)-(h)為在4.35倍欠采樣參數(shù)下，分別使用三階全變分、曲波變換、拉普拉斯高階方法、全變分、二階全變分、各向同性二階全變分和Lysaker高階方法這七種方法對其進行重建的效果圖。可以看出，三階全變分由于圖像導(dǎo)數(shù)階數(shù)較高，因此造成了圖像過度模糊，信噪比較低；全變分導(dǎo)致了重建結(jié)果出現(xiàn)片狀偽影；曲波變換和拉普拉斯高階方法在重建圖像中產(chǎn)生曲線狀和圓形偽影，因此重建效果不佳。綜合考慮，二階全變分對圖像重建效果最好。(i)-(l)為2倍欠采樣參數(shù)下分別用全變分、二階全變分、各向同性二階全變分和Lysaker高階方法這四種方法的重建結(jié)果圖，結(jié)合箭頭所指位置進行比對，二階全變分優(yōu)于其他方法。

圖4和圖5分別為對圖2(d)血管造影磁共振圖像和2(c)腕部磁共振圖像進行3倍欠采樣同時添加35dB噪聲，對圖2-(c)進行2.86倍欠采樣同時添加30dB噪聲后全變分和二階全變分方法重建結(jié)果比對。圖4(c)、4(f)、5(c)、5(f)是對原圖像和重建結(jié)果圖像做差值得出的誤差圖，從中可以看出，本文所提出的二階全變分方法無論在定性方面或者定量方面都優(yōu)于全變分正則化方法。

在圖6(a)中，通過在迭代過程中改變β參數(shù)來觀察其對收斂速度的影響，兩紅顏色折線代表的是β初始迭代值為15，在每一次迭代中分別增大2倍和5倍的誤差收斂趨勢；兩藍顏色折線是初始迭代值為50，在每一次迭代中分別增大2倍和5倍的誤差收斂趨勢。從此圖中可看出，β的初始值以及每次迭代的增量大小并不會對收斂速度產(chǎn)生特別明顯的影響。因此，只要滿足β值隨每次迭代而增加，都能在一定迭代次數(shù)內(nèi)獲得較好的重建結(jié)果。圖6(b)為本發(fā)明所使用算法的CPU運行時間。其中淺顏色與深顏色曲線分別為使用本發(fā)明方法對尺寸為256×256和512×512的磁共振圖像進行重建的信噪比改善趨勢，從對應(yīng)坐標中可看出，本發(fā)明所使用的算法CPU運行時間非常短，是一種高效的算法。

表1為部分實驗數(shù)據(jù)，若以結(jié)果圖的SNR為衡量重建質(zhì)量的標準，可看出高階變分方法較其他方法效果更好。表1為利用不同方法對圖2中四幅參考磁共振圖像在不同欠采樣在不同欠采樣參數(shù)和噪聲水平下重建圖像的信噪比值列表。

表1在不同采樣參數(shù)和噪聲水平下利用不同方法重建磁共振圖像的SNR列表