一種基于多尺度低秩模型的心臟磁共振成像方法
【專利摘要】本發(fā)明公開了一種基于多尺度低秩模型的心臟磁共振成像方法，為心臟磁共振成像研究探索出新的方法；其采用徑向采樣軌跡實現心臟K空間數據的欠采樣，可實現全心臟數據徑向欠采樣，加快磁共振數據采集速度，降低磁共振設備的掃描時間；其基于多尺度低秩分解的心臟磁共振數據的最小化表示，能夠提高磁共振成像的精度；其通過交替方向乘子方法求解重構圖像的凸優(yōu)化問題，能夠提高磁共振成像重構的速度，且重構的圖像結果紋理清晰、邊緣光滑。
【專利說明】
一種基于多尺度低秩模型的心臟磁共振成像方法
技術領域
[0001] 本發(fā)明屬于磁共振成像技術領域，具體涉及一種基于多尺度低秩模型的心臟磁共 振成像方法。
【背景技術】
[0002] 心血管疾病是導致心臟突然停止跳動的主要原因。當前心血管病發(fā)病和死亡率居 高不下，使心血管病防治負擔加重，成為重要公共衛(wèi)生問題，加強心血管病防治刻不容緩。 心臟磁共振成像（Cardiac Magnetic Resonance Imaging，CMR)是利用核磁共振原理進行 人體心臟斷層成像的技術，能準確地反映心臟的解剖結構、形態(tài)功能、血流特性和心肌活 性，已迅速發(fā)展成為心臟疾病診斷中的主要工具。心臟磁共振成像在成像過程中具有較好 的軟組織對比度，沒有任何放射性污染，分辨率高，可任意層面成像;而且由于參與磁共振 成像的因素較多，得到的圖像信息量大，優(yōu)于現有的其它各種影像學成像技術，在心臟疾病 診斷中有很大的優(yōu)越性和應用潛力。
[0003] 在早中期的磁共振成像中，設備掃描需時較長，如一次心臟掃描需要大約數小時 左右，甚至更長時間，這限制了心臟磁共振成像的應用范圍。磁共振成像大多被用在靜態(tài)部 位的成像中，通常不適用于動態(tài)成像，這是因為在動態(tài)成像中，數據量的增加會大幅度延長 數據獲取時間，需時會更長，病人常常會因為長時間的等待數據掃描感覺不適;或在數據獲 取過程中發(fā)生自主和非自主的運動，導致圖像中出現各種偽影。提高磁共振成像速度，縮短 磁共振成像掃描時間的意義不僅僅在于提高磁共振設備的工作效率，減輕病人痛苦，更重 要的是它有利于消除心臟運動以及呼吸等造成運動偽影的影響。
[0004] 目前所研究的磁共振圖像重構方法主要兩種：一種是多線圈并行成像技術，主要 是利用相控陣線圈中單個接收線圈的空間敏感度差異來編碼空間信息，降低成像所必需的 梯度編碼步數，其采用多線圈陣列同時采集信號，允許對K空間進行欠采樣以減少相位編碼 步數，在保持圖像空間分辨率不變的同時，能大幅度縮短掃描時間，提高成像速度;但多線 圈并行成像技術涉及多線圈敏感度分布的估計，需要增加計算量。
[0005] 另一種基于壓縮感知理論的磁共振成像重構方法，由于磁共振影像具有低秩性和 稀疏性，可以采用壓縮感知理論從隨機欠采樣的k空間數據進行圖像重構，減少采樣數據， 提高成像速度，目前常用的有SLR(k_t Sparsity and Low_Rank)、L+S(Low Rank plus Sparse)等低秩變換方法，這些方法雖然考慮到了低秩特性，但沒有考慮磁共振圖像的多尺 度低秩特性。

【發(fā)明內容】

[0006] 針對現有技術所存在的上述技術問題，本發(fā)明提供了一種基于多尺度低秩模型的 心臟磁共振成像方法，該方法將矩陣分解成多尺度的塊低秩矩陣之和，并將多尺度塊低秩 矩陣之和的最小化作為約束條件用于磁共振成像，可提高磁共振成像精度。
[0007] -種基于多尺度低秩模型的心臟磁共振成像方法，包括如下步驟：
[0008] (1)利用徑向采樣軌跡模式對人體心臟的磁共振K空間數據進行欠采樣，得到心臟 的Κ空間欠采樣數據；
[0009] (2)對所述的Κ空間欠采樣數據進行傅里葉逆變換，得到初始心臟磁共振圖像；
[0010] (3)通過將心臟磁共振圖像分解成對應多尺度的若干個塊低秩矩陣，并建立重構 心臟磁共振圖像的目標函數如下：
[0012] 其中：F為系統(tǒng)矩陣（該矩陣由傅立葉變換矩陣和磁共振Κ空間數據欠采樣掩模構 成），f為κ空間欠采樣數據，II ||2為2范數，II II*為核范數，Xi為第i個塊低秩矩陣，i為自 然數且l<i<L，L為塊低秩矩陣的個數且為大于1的自然數，Μ為對應&的正則化系數；
[0013] (4)根據初始心臟磁共振圖像對上述目標函數進行最小化求解，得到L個塊低秩矩 陣心~1，進而對這些塊低秩矩陣進行求和即重建得到心臟磁共振圖像。
[0014] 所述重建得到的心臟磁共振圖像、Κ空間欠采樣數據、初始心臟磁共振圖像以及塊 低秩矩陣Xi的尺度大小均為Μ X Ν，Μ和Ν均為大于1的自然數。
[0015] 所述的塊低秩矩陣X,被分割由多個塊矩陣組成，±夬低秩矩陣中的塊矩陣對 應為不同尺度大小的方陣(如1X1、4X 4、16X16或64X64)。
[0016] 所述的步驟（4 )中米用ADMM算法（A11ernat ing direc t i on me thod of mu 11 i p 1 i er，交替方向乘子算法)對目標函數進行最小化求解。
[0017] 所述的ADMM算法基于以下迭代方程：
[0020] 其中：ΧΓ為塊低秩矩陣Xi對應第k+Ι次迭代的結果，Y〇為初始心臟磁共振圖像，zf 和Zf1分別為伴隨矩陣冗,對應第k次迭代和第k+1次迭代的結果，uf和uf1分別為伴隨矩陣 Ui對應第k次迭代和第k+Ι次迭代的結果，初始化zf和Uf均為與塊低秩矩陣Xi尺度大小相同 的零矩陣；表示以λ,/ρ為軟閾值因子對應括號內目標矩陣為 + 的塊形變SVT(奇異值閾值)函數，Ρ為預設的收斂系數，k為迭代次數。
[0021] 所述塊形變SVT函數&^51^^< +1 +〇的具體計算過程如下：
[0022] 〈1>從目標矩陣X^+Uf1中提取任一塊矩陣，通過塊形變因子將其形變成指定大 小為nu Xm的小塊b;目標矩陣+U"中的塊矩陣與塊低秩矩陣X,中的塊矩陣尺度大小 一致；
[0023] 〈2>對小塊b進行奇異值閾值處理后得到小塊fe，將小塊通過塊形變因子逆過程 形變回塊矩陣的尺度大小并嵌入目標矩陣Χ? +1 + Ι?+1中的原有位置，進而將目標矩陣 XT1 +1^1中其余元素值均置為ο，對應得到一個數據矩陣；
[0024] 〈3>根據步驟〈1>和〈2>遍歷目標矩陣Xf +11^中所有塊矩陣，得到對應的多個數 據矩陣，進而將這些數據矩陣進行累加即得到塊形變SVT函數
[0025] 所述的步驟〈2>中對小塊b進行奇異值閾值處理，具體過程如下：
[0026]首先，根據公式b = USVT對小塊b進行奇異值分解，得到兩個正交矩陣U和V以及一 個對角矩陣S，T表示轉置；
[0027] 然后，使對角矩陣S軟閾值化：即將對角矩陣S中的每個對角線元素均減去λ^ρ，對 于減去Μ/ρ后的對角線元素，若其小于0則將其置為0;
[0028] 最后，根據公式= U§VT計算出小塊?)，§為對角矩陣S軟閾值化后的結果。 [0029]所述的正則化系數\由以下算式確定：
[0031] 本發(fā)明基于多尺度低秩模型的心臟磁共振成像方法為心臟磁共振成像研究探索 出新的方法，其基于徑向采樣模式實現心臟Κ空間數據的欠采樣，可實現全心臟數據徑向欠 采樣，加快磁共振數據采集速度，降低磁共振設備的掃描時間；其基于多尺度低秩分解的心 臟磁共振數據的最小化表示，提高磁共振成像的精度;其基于ADMM算法的成像研究，提高磁 共振成像重構的速度。
【附圖說明】
[0032] 圖1為本發(fā)明心臟磁共振成像方法的流程示意圖。
[0033] 圖2為二維徑向采樣軌跡模式的示意圖。
[0034]圖3為塊形變因子的轉換示意圖。
[0035] 圖4(a)為采用數據集Α對應生成的心臟磁共振圖像。
[0036] 圖4(b)為采用數據集B對應生成的心臟磁共振圖像。
[0037]圖5為加速因子R=3情況下本發(fā)明MSL算法與現有k-t SLR和L+S方法對數據集A進 行心臟磁共振成像關于信號誤差比SER的對比示意圖。
[0038]圖6為加速因子R=8情況下本發(fā)明MSL算法與現有k-t SLR和L+S方法對數據集B進 行心臟磁共振成像關于信號誤差比SER的對比示意圖。
【具體實施方式】
[0039] 為了更為具體地描述本發(fā)明，下面結合附圖及【具體實施方式】對本發(fā)明的技術方案 進行詳細說明。
[0040] 如圖1所示，本發(fā)明基于多尺度低秩模型的心臟磁共振成像方法，具體實施步驟如 下：
[0041] (1)利用徑向采樣軌跡實現K空間數據的欠采樣。一個采樣層面的徑向采樣軌跡包 含Μ條投影線，每條投影線上包含有N個樣本點。選定某個采樣點之后，如圖2所示，會形成二 維的徑向采樣軌跡圖。
[0042] (2)利用塊形變因子Rb()，從目標矩陣中獲取塊大小b(大小為1 X 1、4 X 4、16 X 16 或64X64)并將其形變成指定大小后，而其伴隨矩陣/?,!()將形變后的矩陣轉換成目標矩陣 中的塊矩陣大小，將其嵌入全零矩陣中，如圖3所示。由此建立一個由塊低秩矩陣Xi組成的 多尺度低秩模型如下：
[0045] 其中：Ub，Sb和Vb是大小分別為nu Xrb，rb X rb和m X rb的矩陣，三者形成Rb(Xi)秩為 η的簡化奇異值分解。
[0046] (3)心臟磁共振重構公式表示如下：
[0048] 其中：F為系統(tǒng)矩陣，f為Κ空間欠采樣數據，I I I |2為2范數，I I I I*為核范數，Xi為 第i個塊低秩矩陣，i為自然數且l<i<L，L為塊低秩矩陣的個數且為大于1的自然數，λ?為 對應Xi的正則化系數。
[0049] 對上述目標函數優(yōu)化求解的詳細過程如下：
[0050] 首先，設定初始狀態(tài);將欠采樣的K空間數據進行傅里葉逆變換得到初始圖像Yo。
[0051] 然后，對形成的凸優(yōu)化方程進行約束，我們選擇使用正則化系數λ,，表示對i尺度 矩陣的正則化參數，定義如下：
[0053]其中，MXN表示初始心臟磁共振圖像Yo的大小，nu Xm表示從X矩陣中將b塊形變的 大小，隨后將nu Xm大小的矩陣嵌入到全零矩陣中。
[0054]進而根據交替方向乘子法（ADMM，Alternating direction method of 1111111^?1161')重構上述的基于多尺度低秩模型的磁共振圖像，所用的4011算法如下：
[0057]其中為塊低秩矩陣Xi對應第k+1次迭代的結果，YQ為初始心臟磁共振圖像，Zj1 和Ζ廣分別為伴隨矩陣Ζ,對應第k次迭代和第k+1次迭代的結果，Df和Uf+1分別為伴隨矩陣 山對應第k次迭代和第k+Ι次迭代的結果，初始化Zf和Uf均為與塊低秩矩陣&尺度大小相同 的零矩陣；
表示以λ,/ρ為軟閾值因子對應括號內目標矩陣為 xp+1J"的塊形變SVT(奇異值閾值)函數，p為預設的收斂系數，設定 p=20，k為迭代次數。
[0058] SVI\(X)和 (X)的函數表達式為：
[0059] SVTa(X)=U max(E-A,0)VT
[0061]其中，U和V表示正交陣，而Σ表示對角陣，λ為奇異值軟閾值因子，Rb〇為塊形變因 子，其作用是從目標矩陣中獲取塊大小b并將其形變成指定大小，表示伴隨矩陣，Pi表示 的是Xi所對應的分塊低秩部分，塊大小b隨著Pi的不同而變化。
[0062]最后，判斷迭代收斂，或直至達到最大迭代次數k = 700為止。
[0063] 為了測試多尺度低秩模型的磁共振圖像重構方法的重建性能，利用兩種不同的MR 原始數據集進行研究，即心臟灌注和心臟電影。數據集A:心臟灌注MR數據集，是在Utah大學 采用3T西門子掃描儀飽和恢復序列（磁共振掃描參數TR/TE為=2.5/1毫秒，飽和度恢復時 間為1〇〇毫秒），它含有70幀的時間序列數據，采集每幀的笛卡爾網格具有90 X190(相位編 碼X頻域編碼）的K空間矩陣。數據集B:心臟電影MR數據，在Yonsei大學醫(yī)療中心采用3T飛 利浦磁共振掃描儀采集得到，它是由具有256 X 256面內的空間分辨率的25幀時間序列數據 組成(本實驗中只取其一部分:64X64的25幀時間序列數據）。獲取數據所用成像序列是一 個具有60°翻轉角度以及磁共振掃描參數TR為3.45ms的穩(wěn)態(tài)自由進動(SSFP，Steady State Free Precession)序列，視野（F0V，Field of View)為345mmX6270mm，切片厚度為10mm。
[0064] 兩組數據集采用本發(fā)明方法對應生成的心臟磁共振圖像如圖4所示。
[0065]以下將本發(fā)明多尺度低秩磁共振成像方法與k-t SLR、L+S方法進行比較，從多個 方面評價磁共振重構圖像的質量。首先利用信差比、相對誤差評價重構算法的性能;其次對 重構圖像及誤差圖，從視覺角度評價重構質量。其中信號誤差比（SER, signal to error ratio)和相對誤差(Relative Error，RE)定義如下：
[0067] 表1是兩組數據六種不同的加速因子所對應的三種方法的SER值，可以發(fā)現，本發(fā) 明基于MSL (多尺度低秩)算法的SER值最大。
[0068] 表 1
[0070]本發(fā)明MSL算法與另外兩種方法在重構的心臟磁共振圖像結果及其性能上分別如 圖5和圖6所示(其中數據集A的加速因子R = 3,數據B的加速因子R = 8);從實驗結果可以看 出本發(fā)明MSL算法相對現有的k-t SLR、L+S方法，在磁共振成像的精度上有顯著的提高。 [0071]上述的對實施例的描述是為便于本技術領域的普通技術人員能理解和應用本發(fā) 明。熟悉本領域技術的人員顯然可以容易地對上述實施例做出各種修改，并把在此說明的 一般原理應用到其他實施例中而不必經過創(chuàng)造性的勞動。因此，本發(fā)明不限于上述實施例， 本領域技術人員根據本發(fā)明的揭示，對于本發(fā)明做出的改進和修改都應該在本發(fā)明的保護 范圍之內。
【主權項】
1. 一種基于多尺度低秩模型的屯、臟磁共振成像方法，包括如下步驟： (1) 利用徑向采樣軌跡模式對人體屯、臟的磁共振K空間數據進行欠采樣，得到屯、臟的K 空間欠采樣數據； (2) 對所述的K空間欠采樣數據進行傅里葉逆變換，得到初始屯、臟磁共振圖像； (3) 通過將屯、臟磁共振圖像分解成對應多尺度的若干個塊低秩矩陣，并建立重構屯、臟 磁共振圖像的目標函數如下： 、.一 其中：F為系統(tǒng)矩陣，f為K空間欠采樣數據，II IU為2范數，II II*為核范數，Xi為第i個 塊低秩矩陣，i為自然數且1《1《^1^為塊低秩矩陣的個數且為大于1的自然數，、為對應乂1 的正則化系數； (4) 根據初始屯、臟磁共振圖像對上述目標函數進行最小化求解，得到L個塊低秩矩陣Xi ~Xl，進而對運些塊低秩矩陣進行求和即重建得到屯、臟磁共振圖像。2. 根據權利要求1所述的屯、臟磁共振成像方法，其特征在于:所述重建得到的屯、臟磁共 振圖像、K空間欠采樣數據、初始屯、臟磁共振圖像W及塊低秩矩陣Xi的尺度大小均為ΜΧΝ,Μ 和Ν均為大于1的自然數。3. 根據權利要求4所述的屯、臟磁共振成像方法，其特征在于:所述的塊低秩矩陣Xi被分 割由多個塊矩陣組成，塊低秩矩陣Xi~Xl中的塊矩陣對應為不同尺度大小的方陣。4. 根據權利要求3所述的屯、臟磁共振成像方法，其特征在于：所述的步驟（4)中采用 ADMM算法對目標函數進行最小化求解。5. 根據權利要求4所述的屯、臟磁共振成像方法，其特征在于:所述的ADMM算法基于W下 迭代方程：其中為塊低秩矩陣Xl對應第k+1次迭代的結果，1^0為初始屯、臟磁共振圖像，巧和 巧"分別為伴隨矩陣Zi對應第k次迭代和第k+1次迭代的結果，Uf和Uf"分別為伴隨矩陣化 對應第k次迭代和第k+1次迭代的結果，初始化Zf和Uf均為與塊低秩矩陣Xl尺度大小相同的 零矩陣：表示Κλι/ρ為軟闊值因子對應括號內目標矩陣為 XfTi+Uf的塊形變SVT函數，Ρ為預設的收斂系數，k為迭代次數。6. 根據權利要求5所述的屯、臟磁共振成像方法，其特征在于：所述塊形變SVT函數的具體計算過程如下： <1〉從目標矩陣Xf+i+llf中提取任一塊矩陣，通過塊形變因子將其形變成指定大小為 mi Xm的小塊b;目標矩陣餐+1 +llfi中的塊矩陣與塊低秩矩陣Xi中的塊矩陣尺度大小一致； <2〉對小塊b進行奇異值闊值處理后得到小塊b，將小塊Ι?通過塊形變因子逆過程形變回 塊矩陣的尺度大小并嵌入目標矩陣Xf+i +Uf^i中的原有位置，進而將目標矩陣+Uf 1中 其余元素值均置為0，對應得到一個數據矩陣； <3〉根據步驟<1〉和<2〉遍歷目標矩陣Xfi+Uf"中所有塊矩陣，得到對應的多個數據矩 陣，進而將運些數據矩陣進行累加即得到塊形變SVT函數7. 根據權利要求6所述的屯、臟磁共振成像方法，其特征在于:所述的步驟<2〉中對小塊b 進行奇異值闊值處理，具體過程如下： 首先，根據公式b=usyT對小塊b進行奇異值分解，得到兩個正交矩陣U和VW及一個對角 矩陣S，T表示轉置； 然后，使對角矩陣S軟闊值化：即將對角矩陣S中的每個對角線元素均減去λι/ρ，對于減 去λι/ρ后的對角線元素，若其小于0則將其置為0; 最后，根據公式b二U§vT計算出小塊I)，S為對角矩陣S軟闊值化后的結果。8. 根據權利要求7所述的屯、臟磁共振成像方法，其特征在于:所述的正則化系數λι由W 下算式確定：
【文檔編號】G01R33/20GK106093814SQ201610387992
【公開日】2016年11月9日
【申請日】2016年6月2日 公開號201610387992.7, CN 106093814 A, CN 106093814A, CN 201610387992, CN-A-106093814, CN106093814 A, CN106093814A, CN201610387992, CN201610387992.7
【發(fā)明人】蔣明峰, 陸雨, 黃文清, 馮杰, 鄭俊褒
【申請人】浙江理工大學