[發(fā)明專利]基于奇異譜分析的SPEED快速磁共振成像方法

在審

• 申請?zhí)枺?/td>: 201711337327.8
• 申請日: 2017-12-14
• 公開（公告）號: CN108051766A
• 公開（公告）日: 2018-05-18
• 發(fā)明（設(shè)計）人: 金朝陽;向清三
• 申請（專利權(quán)）人: 杭州電子科技大學
• 主分類號: G01R33/48
• 分類號: G01R33/48;G01R33/56;G01R33/561
• 代理公司: 杭州君度專利代理事務(wù)所(特殊普通合伙) 33240
• 代理人: 朱月芬

1.基于奇異譜分析的SPEED快速磁共振成像方法,其特征在于，該方法包括以下步驟：

步驟1：k空間數(shù)據(jù)采集

在k空間的相位編碼方向，即PE方向，每隔N行采集一行數(shù)據(jù)，共采集兩組，分別用S1和S2表示，用kx表示在k空間頻率編碼FE方向的位置，ky表示在PE方向的位置，則S1(kx,ky)和S2(kx,ky)分別表示S1和S2中的一個數(shù)據(jù)點的值，用d1,d2表示每組欠采樣數(shù)據(jù)在PE方向上的偏移量，采樣方式用N(d1,d2)表示；根據(jù)圖像大小，在PE方向的k空間中心區(qū)域采集16至64行數(shù)據(jù)，用Sc表示，Sc(kx,ky)表示Sc中的一個數(shù)據(jù)點；在FE方向的數(shù)據(jù)都為全采集；

步驟2：填零重建

填零重建包含3個步驟：兩組欠采樣數(shù)據(jù)的填零重建、差分變換、k空間中心數(shù)據(jù)填零重建成低分辨率圖像；

步驟2-1：兩組欠采樣數(shù)據(jù)的填零重建

對于兩組欠采樣數(shù)據(jù)S1和S2，其對應(yīng)k空間中沒有進行數(shù)據(jù)采集的點用0表示，進行常規(guī)的填零傅立葉重建，重建后圖像分別用I1和I2表示，用(x,y)表示圖像空間中每個像素點的坐標值，則I1(x,y)和I2(x,y)表示I1和I2中的一個像素點的值；由于k空間中每隔N行采集一行數(shù)據(jù)，因此每組數(shù)據(jù)對應(yīng)的填零傅立葉重建圖像中有N層重疊的鬼影，每個像素點上最多有N層重疊的鬼影；

步驟2-2：差分變換

對步驟2-1得到的圖像I1和I2分別進行差分變換，得到稀疏的邊緣鬼影圖像E1和E2；稀疏的鬼影圖像E1和E2中的每個像素點E1(x,y)和E2(x,y)分別表示為：

E1(x,y)＝I1(x,y)-I1(x,y-1) [1]

E2(x,y)＝I2(x,y)-I2(x,y-1) [2]

步驟2-3：k空間中心數(shù)據(jù)填零重建成低分辨率圖像

采集到的k空間中心部分數(shù)據(jù)Sc，其對應(yīng)k空間中沒有進行數(shù)據(jù)采集的點用0表示，進行常規(guī)的填零傅立葉重建，形成一個低分辨率的重建圖像Ilow，其對應(yīng)的每一個像素點用Ilow(x,y)表示；

步驟3：奇異譜分析重建

奇異譜分析重建包含8個步驟：差分變換、確立采集矩陣、尋找奇異點、層析法更新差分圖像、構(gòu)建奇異譜函數(shù)、計算奇異度、重建出k空間缺失數(shù)據(jù)、基于離散傅立葉逆變換的圖像重建；

步驟3-1：差分變換

對填零重建圖像Ilow在y方向上進行差分變換，得到稀疏的差分圖像ΔIlow，對應(yīng)的每個像素點表示為ΔIlow(x,y)：

ΔIlow(x,y)＝Ilow(x,y)-Ilow(x,y-1) [3]

閾值T設(shè)為差分圖像ΔIlow中幅值最大值的十分之一；

步驟3-2：確立采集矩陣

采集矩陣Mask中的每一點Mask(kx,ky)表示為：

采集矩陣對應(yīng)的頻域點表示為：

i(x,y)＝IDFT(Mask(kx,ky)) [5]

步驟3-3：尋找奇異點

取ΔIlow中幅值最大的點為奇異點j(xj,yj)，加入奇異點集SP中；把i(x,y)矩陣的原點移動至奇異點j(xj,yj)處，更新采集矩陣的頻域表示為i(x-xj,y-yj)；

步驟3-4：層析法更新差分圖像

更新后的差分圖像ΔIlow(x,y)new表示為，

ΔIlow(x,y)new＝ΔIlow(x,y)-α·i(x-xj,y-yj) [6]

α＝ΔIlow(x,y)||Mask||/(M×N) [7]

其中M和N表示x和y方向的像素點數(shù)，||Mask||表示采集矩陣中非零點的總數(shù)；

如果ΔIlow(x,y)new≥T，則令ΔIlow(x,y)＝ΔIlow(x,y)new，轉(zhuǎn)到步驟3-3繼續(xù)尋找奇異點；如果ΔIlow(x,y)new＜T，則停止奇異點的尋找，奇異點集SP的奇異點的總個數(shù)為Q，繼續(xù)執(zhí)行步驟3-5；

步驟3-5：構(gòu)建奇異譜函數(shù)

對奇異點集SP中的奇異點分別構(gòu)建奇異譜函數(shù)Wj(kx,ky)，用公式表示為：

Wj(kx,ky)＝IDFT(δ(y-yj)u(x-xj)) [8]

其中，

IDFT表示離散傅立葉逆變換；

步驟3-6：計算奇異度

利用已采集到的k空間數(shù)據(jù)Sc通過求解方程[11]來得到各個奇異點對應(yīng)的奇異度aj，

步驟3-7：重建出k空間缺失數(shù)據(jù)

用奇異譜函數(shù)和奇異度來估算出k空間中沒有進行過數(shù)據(jù)采集處的數(shù)據(jù)Sr(kx,ky)：

步驟3-8：基于離散傅立葉逆變換的圖像重建

對于k空間中進行過數(shù)據(jù)采集的點用實際采集到的數(shù)據(jù)Sc(kx,ky)，沒有進行過采集的點用步驟3.7估算得到的Sr(kx,ky)，然后利用離散傅立葉逆變換進行圖像重建，得到奇異譜重建圖像ISFA，其對應(yīng)的像素點表示為，

步驟4：鬼影定位

鬼影定位包含4個步驟：合成定位圖像、差分變換、建立重疊鬼影圖和確立鬼影階數(shù)；

步驟4-1：合成定位圖像

將步驟2-3得到的低分辨率圖像Ilow與步驟3-8得到的重建圖像ISFA進行對應(yīng)像素點之間的相乘操作，合成用于鬼影定位的定位圖像Ic，對應(yīng)的像素點表示為：

Ic(x,y)＝Ilow(x,y)×ISFA(x,y) [14]

步驟4-2：差分變換

對Ic進行差分變換，得到稀疏的邊緣鬼影圖像Ec；

步驟4-3：建立重疊鬼影圖

在相位編碼方向?qū)c分別進行長度為Ny×n/N的平移，其中Ny表示沿PE方向的數(shù)據(jù)矩陣的大小，n表示邊緣鬼影的階數(shù)，n＝0,1,2,…,N-1；這n個稀疏的邊緣鬼影相加后形成一個重疊的鬼影映射圖Ec,n；

步驟4-4：確立鬼影階數(shù)

在鬼影映射圖Ec,n中，為每個像素點找出兩個最強的鬼影，并記錄下它們對應(yīng)的鬼影階數(shù)(n1,n2)；

步驟5：基于雙層鬼影模型的SPEED圖像重建

基于雙層鬼影模型的SPEED圖像重建包含4個步驟：雙層鬼影模型求解、重疊鬼影的分離、多個分離鬼影圖像的配準求和、逆濾波重建；

步驟5-1：雙層鬼影模型求解

稀疏邊緣鬼影圖像E1和E2中，采用雙層稀疏邊緣鬼影模型來描述，表示為：

公式[15]中為相位因子，Gn1和Gn2分別為每個像素點上需要確定的不同階的鬼影，n1和n2分別表示不同的鬼影階數(shù)；定義為：

公式[16]中d表示每組欠采樣數(shù)據(jù)在PE方向上的偏移量d1和d2，n為鬼影階數(shù)；

在公式[15]中，由于E1、E2、d和N已知，步驟3-3得到了鬼影階數(shù)(n1,n2)，即公式[15]中的兩個方程，僅有兩個未知數(shù)Gn1和Gn2，因此直接解出公式[15]中的兩個重疊的鬼影Gn1和Gn2；

步驟5-2：重疊鬼影的分離

對步驟5-1得到的Gn1和Gn2中的像素點，按不同的鬼影階數(shù)n進行分類，產(chǎn)生N個分離的鬼影映射圖Gn,其中n＝0,1,…,N-1；

步驟5-3：多個分離鬼影圖像的配準求和

步驟5-2得到的N個鬼影映射圖Gn，各自對應(yīng)的鬼影位置不同，通過移位和對齊來配準；配準后各鬼影圖對應(yīng)的像素點求和后，得到?jīng)]有重疊鬼影的邊緣映射圖像E0；

步驟5-4：逆濾波重建

步驟5-3得到的邊緣映射圖像E0經(jīng)離散傅立葉變換到k空間，其對應(yīng)k空間中實際進行數(shù)據(jù)采集的點的值用實際采集的數(shù)據(jù)替代，得到k空間數(shù)據(jù)R0；基于逆濾波公式[17]重建出最終的SPEED圖像I0；

公式[17]中ky表示沿PE方向的k空間位置。