CN110658484A - 一种磁共振波谱重建方法及系统 - Google Patents
一种磁共振波谱重建方法及系统 Download PDFInfo
- Publication number
- CN110658484A CN110658484A CN201910990393.8A CN201910990393A CN110658484A CN 110658484 A CN110658484 A CN 110658484A CN 201910990393 A CN201910990393 A CN 201910990393A CN 110658484 A CN110658484 A CN 110658484A
- Authority
- CN
- China
- Prior art keywords
- matrix
- magnetic resonance
- rank
- singular
- iteration
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R33/00—Arrangements or instruments for measuring magnetic variables
- G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
- G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
- G01R33/46—NMR spectroscopy
- G01R33/4625—Processing of acquired signals, e.g. elimination of phase errors, baseline fitting, chemometric analysis
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R33/00—Arrangements or instruments for measuring magnetic variables
- G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
- G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
- G01R33/46—NMR spectroscopy
- G01R33/4633—Sequences for multi-dimensional NMR
Landscapes
- Physics & Mathematics (AREA)
- Spectroscopy & Molecular Physics (AREA)
- High Energy & Nuclear Physics (AREA)
- Condensed Matter Physics & Semiconductors (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Signal Processing (AREA)
- Magnetic Resonance Imaging Apparatus (AREA)
Abstract
本发明涉及一种磁共振波谱重建方法及系统,该方法包括:首先获取欠采样矩阵X0;然后对欠采样矩阵X0进行奇异值分解,得到X0=UVT;再将U和V的每一列都转化为汉克尔矩阵,并构建磁共振波谱低秩模型;接着引入中间变量将磁共振波谱低秩模型转化为最小秩优化问题;采用交替方向乘子法求解最小秩优化问题,得到Uk+1和Vk+1,并补全矩阵X:最后当Xk+1为二维磁共振信号时,对Xk+1做二维傅里叶变换得到完整的二维磁共振波谱。本发明利用磁共振信号可以进行范德蒙德分解的性质,对欠采样矩阵X0进行奇异值分解得矩阵U和V,通过采用交替方向乘子法分别求解U和V,完成对U和V的重建,从而实现对欠采样矩阵X0的补全,重建速度快且精度高。
Description
技术领域
本发明涉及属于磁共振波谱重建领域,具体涉及一种磁共振波谱重建方法及系统。
背景技术
核磁共振波谱法是研究原子核对射频辐射的吸收,它是对各种有机和无机物的成分、结构进行定性分析的最强有力的工具之一,有时亦可进行定量分析。核磁共振技术可以提供分子的化学结构和分子动力学的信息,已成为分子结构解析以及物质理化性质表征的常规技术手段,它可以确定几乎所有常见官能团的环境,直观性强,特别是碳谱能直接反映出分子的骨架,谱图解释较为容易,在物理、化学、生物、医药、食品等领域得到广泛应用,在化学中更是常规分析不可少的手段。
在实际应用中,由于较高的实验成本和硬件限制,以及其它不可避免的原因,我们很难对二维或三维磁共振自由衰减信号实现完全采样,因此对磁共振信号进行欠采样并通过欠采样得到的数据重建出完整的磁共振信号显得尤为重要。不同的重建方法在重建精度、最小欠采样率和重建速度上有较大差异,而目前已有的方法都需要较长的重建时间,尤其重建三维磁共振波谱时,所需的时间成比例增加。为此,提供一种重建速度更快且有较高重建精度的新的重建方法成为当前需要解决的问题。
发明内容
(一)要解决的技术问题
本发明提供了一种磁共振波谱重建方法,旨在提供一种重建速度更快且有较高重建精度的新的重建方法。
(二)技术方案
为了解决上述问题,本发明提供一种磁共振波谱重建方法,包括以下步骤:
S1:获取欠采样矩阵X0;
S2:对所述欠采样矩阵X0进行奇异值分解,得到X0=U0ΣV0 T,其中,U0为左奇异向量,V0为右奇异向量,Σ为奇异值;
取所述欠采样矩阵X0的前r个奇异值和与所述r个奇异值对应的左奇异向量和右奇异向量,令Σ=ΣLΣR,U=U0ΣL,VT=ΣRV0 T则X0=UVT;
S3:将所述U和所述V的每一列都转化为汉克尔矩阵,并构建磁共振波谱低秩模型:
s.t.PΩ(X)=X0
其中,r表示矩阵U和V的列数;R为将一维数据转成汉克尔矩阵的操作算子;操作算子Qi用来提取矩阵第i列;rank(·)表示矩阵的秩;表示矩阵的F-范数的平方;β为正则化参数,用于权衡前后两项的重要性,β>0;X表示被补全后的矩阵;集合Ω表示采样点的位置的集合;即PΩ与X的哈达玛积表示按照模板PΩ对所述X进行欠采样,表示为PΩ(X),所述矩阵X对应PΩ数字为1的位置的元素保持不变,对应模板数字0的位置的元素被置为0;
S4:引入中间变量Bi=RQiU、Ci=RQiV,将所述磁共振波谱低秩模型转化为最小秩优化问题:
s.t.Bi=RQiU,Ci=RQiV,PΩ(X)=X0
采用交替方向乘子法求解所述最小秩优化问题,得到所述Uk+1和所述Vk+1,并补全所述矩阵X:
X=Xk+1=PΩ(X0)+PΩc(Uk+1(Vk+1)T)
其中Xk+1为所述矩阵X经过k+1次迭代后的矩阵;Uk+1为所述矩阵U经过k+1次迭代后的矩阵;Vk+1为所述矩阵V经过k+1次迭代后的矩阵;PΩc是PΩ的补充,其值与PΩ相反,即PΩ为1的位置,PΩc在相对应的位置为0;PΩ为0的位置,PΩc在相对应的位置为1;
S5:当所述Xk+1为二维磁共振信号时,对所述Xk+1做二维傅里叶变换得到完整的二维磁共振波谱。
优选地,在所述步骤S2中,所述r在实际的二维磁共振波谱谱峰数6倍之内。
优选地,在所述步骤S4具体为:
S41:将所述最小秩优化问题转化为增广拉格朗日函数:
s.t.PΩ(X)=X0
其中,Di,Mi为拉格朗日乘子,μ为正则化参数,μ>0,<·,·>表示矩阵在希尔伯特空间的内积,<A,B>=trace(AB),trace(·)表示矩阵的迹;
S42:对所述增广拉格朗日函数采用交替方向乘子法交替求解以下子问题得到:
s.t.PΩ(X)=X0
S43:通过求解以下最小二乘问题得到U:
S44:求解以下最小二乘问题得到V:
S45:通过以下公式迭代更新X:
Xk+1=X0+PΩc(Uk+1(Vk+1)T)
通过以下公式迭代更新Bi、Ci:
其中,其中k为迭代次数;R*为操作算子R的伴随算子;Qi*为Qi的伴随算子;为奇异值缩减算子,为所述Bi经过k+1次迭代后的结果,为所述Ci经过k+1次迭代后的结果,μk为所述μ经过k次迭代后的结果。
所述Ci的优化子问题表示为:
其中,σBi表示矩阵Bi的奇异值;σCi表示矩阵Ci的奇异值;||·||0表示矩阵的0范数。
优选地,每次迭代所述μ值逐渐增大,当达到迭代停止准则时,增大所述β并继续迭代,直至所述β达到设定值,迭代结束。
优选地,所述β初值设为20-30,所述μ初值为0.005-0.02,当完成一次迭代后,下一次迭代时所述μ的取值值为上一次迭代时μ值的1.02-1.12倍,当达到迭代停止准则时,下一次迭代时所述β的取值为上次迭代时所述β的值的两倍,继续迭代直至所述β达到最大值222-242,迭代结束;迭代停止准则设定为在两次相邻迭代中所补全的矩阵Xk+1与Xk的误差小于阈值η,η设为10-7-10-5。
优选地,在所述步骤S43具体为:
R*为操作算子R的伴随算子,定义如下:
Qi*为Qi的伴随算子,定义如下:
d为向量,表示将向量d的元素放入矩阵[Qi*d]的p列中;
所以,
其中,C表示复数域,M、N分别为矩阵的行数、列数;w1是一个M×1向量且它的第k个元素为所述矩阵RQiU第k条副对角线上的元素个数,w2是一个N×1向量且它的第k个元素为所述矩阵RQiV第k条副对角线上的元素个数;定义矩阵T1为r列w1、T2为r列w2组合而成;
所以,
并通过以下公式通过更新U的每一行来迭代U:
U(a,:) k+1=Y1(a,:)(μkTa+β(Vk)Tconj(Vk))-1
其中,U(a,:) k+1表示矩阵U经过k+1次迭代后产生的矩阵Uk+1的第a行的元素;Y1(a,:)为矩阵Y1第a行的元素;Ta为对角矩阵,Ta主对角线元素为T1(a,:),T1(a,:)为矩阵T1第a行的元素。
优选地,所述步骤S44具体为:
得到:
并通过以下公式通过更新V的每一行来迭代V:
V(b,:) k+1=Y2(b,:)(μkTb+β(Uk+1)Tconj(Uk+1))-1
其中,V(b,:) k+1表示矩阵V经过k+1次迭代后产生的矩阵VK+1第n行的元素;Y2(b,:)为矩阵Y2第n行的元素;Tb为对角矩阵,Tb主对角线元素为T2(b,:),T2(b,:)为矩阵T2第b行的元素。
优选地,所述步骤S5还包括:
当所述为三维磁共振信号时,对组成所述三维磁共振信号的所有二维磁共振信号分别做二维傅里叶变换,并对所述三维磁共振信号做三维傅里叶变换得到完整的三维磁共振波谱。
优选地,本发明还提供了一种磁共振波谱重建系统,包括:获取模块、奇异值分解模块、低秩模型模块、优化模块和变换模块;
所述获取模块用于获取欠采样矩阵X0;
所述奇异值分解模块用于对欠采样矩阵X0进行奇异值分解,得到X0=U0ΣV0 T,其中,U0为左奇异向量,V0为右奇异向量,Σ为奇异值;
取所述欠采样矩阵X0的前r个奇异值和与所述r个奇异值对应的左奇异向量和右奇异向量,且令Σ=ΣLΣR,U=U0ΣL,VT=ΣRV0 T则X0=UVT;
所述低秩模型模块用于将所述U和所述V的每一列都转化为汉克尔矩阵,并构建磁共振波谱低秩模型:
s.t.PΩ(X)=X0
其中,r表示矩阵U和V的列数;R为将一维数据转成汉克尔矩阵的操作算子;操作算子Qi用来提取矩阵第i列;rank(·)表示矩阵的秩;表示矩阵的F-范数的平方;β为正则化参数,用于权衡前后两项的重要性,β>0;X表示被补全后的矩阵;集合Ω表示采样点的位置的集合;即PΩ与X的哈达玛积表示按照模板PΩ对所述X进行欠采样,表示为PΩ(X),这样所述矩阵X对应PΩ数字为1的位置的元素保持不变,对应模板数字0的位置的元素被置为0;
所述优化模块用于引入中间变量Bi=RQiU、Ci=RQiV,将所述磁共振波谱低秩模型转化为最小秩优化问题:
s.t.Bi=RQiU,Ci=RQiV,PΩ(X)=X0
采用交替方向乘子法求解所述最小秩优化问题,得到所述Uk+1和所述Vk+1,并补全所述矩阵X:
X=Xk+1=PΩ(X0)+PΩc(Uk+1(Vk+1)T)
其中Xk+1为所述矩阵X经过k+1次迭代后的矩阵;Uk+1为所述矩阵U经过k+1次迭代后的矩阵;Vk+1为所述矩阵V经过k+1次迭代后的矩阵;PΩc是PΩ的补充,其值与PΩ相反,即PΩ为1的位置,PΩc在相对应的位置为0;PΩ为0的位置,PΩc在相对应的位置为1;
所述变换模块用于当所述Xk+1为二维磁共振信号时,对所述Xk+1做二维傅里叶变换得到其二维频谱。
(三)有益效果
本发明利用将欠采样矩阵进行奇异值分解,再将分解得到的矩阵的每一列都转化为汉克尔矩阵,并建立磁共振波谱低秩模型,通过采用交替乘子的方法获得重建后的矩阵,最后对重建后的矩阵进行傅里叶变换得到完整的磁共振波谱。这种重建方法重建速度更快且有较高重建精度。
附图说明
图1为本发明一种磁共振波谱重建方法的流程图;
图2为本发明中求解最小秩优化问题的流程图;
图3为本发明一种磁共振波谱重建系统的示意图;
图4为本发明实施例中大小为256×128的全采样二维磁共振波谱;
图5为使用本发明的方法对欠采样矩阵补全重建出的磁共振波谱;
图6为本发明与其它方法在20%采样率下重建结果的对比;
图7为采用不同阈值缩减算子在20%采样率下重建结果的对比。
【附图标记说明】
1:获取模块;2:奇异值分解模块;3:低秩模型模块;4:优化模块;5:变换模块。
具体实施方式
为了更好的解释本发明,以便于理解,下面结合附图,通过具体实施方式,对本发明作详细描述。
需要说明,本发明实施例中所有方向性指示(诸如上、下、左、右、前、后……)仅用于解释在某一特定姿态(如附图所示)下各部件之间的相对位置关系、运动情况等,如果该特定姿态发生改变时,则该方向性指示也相应地随之改变。
另外,在本发明中如涉及“第一”、“第二”等的描述仅用于描述目的,而不能理解为指示或暗示其相对重要性或者隐含指明所指示的技术特征的数量。由此,限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本发明的描述中,“多个”的含义是至少两个,例如两个,三个等,除非另有明确具体的限定。
在本发明中,除非另有明确的规定和限定,术语“连接”、“固定”等应做广义理解,例如,“固定”可以是固定连接,也可以是可拆卸连接,或成一体;可以是机械连接,也可以是电连接;可以是直接相连,也可以通过中间媒介间接相连,可以是两个元件内部的连通或两个元件的相互作用关系,除非另有明确的限定。对于本领域的普通技术人员而言,可以根据具体情况理解上述术语在本发明中的具体含义。
本发明提供一种磁共振波谱重建方法,如图1:本发明一种磁共振波谱重建方法的流程图所示,具体包括以下步骤:
S1:获取欠采样矩阵X0。X0为一个M×N的矩阵,M、N分别表示矩阵X0的行数和列数。
S2:对欠采样矩阵X0进行奇异值分解,得到X0=U0ΣV0 T,其中,U0为左奇异向量,V0为右奇异向量,Σ为奇异值。
取欠采样矩阵X0的前r个奇异值和与r个奇异值对应的左奇异向量和右奇异向量,令Σ=ΣLΣR,U=U0ΣL,VT=ΣRV0 T则X0=UVT。U为M×r矩阵,V为N×r矩阵。
一维磁共振信号可以表示为几个复指数衰减信号的叠加,而二维磁共振信号即是几个两两相乘的复指数衰减信号的叠加,由二维磁共振信号构成的矩阵P中第m行,n列的元素表示为:其中,yi=exp(j2πf1i-τ1i),zi=exp(j2πf2i-τ2i),f表示归一化频率,τ为衰减因子,b表示该二维磁共振波谱有g个谱峰;一个完整的二维磁共振信号矩阵P可以被分解为两个范德蒙德矩阵,如下所示:
P=YDZT,将P分为D=DLDR两个对角矩阵,分别乘入Y与ZT,设U1=YDL,V1 T=DRZT,则X=U1V1 T,U1、V1为范德蒙德矩阵。
欠采样矩阵X0也是对磁共振信号欠采样而得到的,在对矩阵X0进行奇异值分解后得到U、V,如果我们使U、V能够逼近范德蒙德矩阵U1、V1,那么就可以重建得到完整的磁共振信号,本申请利用磁共振信号能够分解为两个范德蒙德矩阵的性质,将补全矩阵X0的问题转为使U、V逼近U1、V1,从而使得补全的方法更加准确快速。
S3:将U和V的每一列都转化为汉克尔矩阵,并构建磁共振波谱低秩模型:
s.t.PΩ(X)=X0
其中,r表示矩阵U和V的列数;R为将一维数据转成汉克尔矩阵的操作算子;操作算子Qi用来提取矩阵第i列;rank(·)表示矩阵的秩;表示矩阵的F-范数的平方;β为正则化参数,用于权衡前后两项的重要性,β>0;X表示被补全后的矩阵;集合Ω表示采样点的位置的集合;即PΩ与X的哈达玛积表示按照模板PΩ对X进行欠采样,表示为PΩ(X),矩阵X对应PΩ数字为1的位置的元素保持不变,对应模板数字0的位置的元素被置为0;
设e为1×A维的数据,则Re为:
其中,k设为0.1A,e(f)表示e中的第f个元素。
S4:引入中间变量Bi=RQiU、Ci=RQiV,将磁共振波谱低秩模型转化为最小秩优化问题:
s.t.Bi=RQiU,Ci=RQiV,PΩ(X)=X0
采用交替方向乘子法求解最小秩优化问题,得到Uk+1和Vk+1,并补全矩阵X:
X=Xk+1=PΩ(X0)+PΩc(Uk+1(Vk+1)T)
其中Xk+1为矩阵X经过k+1次迭代后的矩阵;Uk+1为矩阵U经过k+1次迭代后的矩阵;Vk+1为矩阵V经过k+1次迭代后的矩阵;PΩc是PΩ的补充,其值与PΩ相反,即PΩ为1的位置,PΩc在相对应的位置为0;PΩ为0的位置,PΩc在相对应的位置为1;
S5:当Xk+1为二维磁共振信号时,对Xk+1做二维傅里叶变换得到完整的二维磁共振波谱。
进一步地,在步骤S2中,r在实际的二维磁共振波谱谱峰数6倍之内。在一些实际情况中不确定磁共振信号的谱峰数,故需预设谱峰r,预设数在实际谱峰数6倍之内都可以成功重建,即满足r<6g。
如图2:本发明中求解最小秩优化问题的流程图所示,在优选的实施例中,步骤S4具体为:
S41:将最小秩优化问题转化为增广拉格朗日函数:
s.t.PΩ(X)=X0
其中,Di,Mi为拉格朗日乘子,μ为正则化参数,μ>0,<·,·>表示矩阵在希尔伯特空间的内积,<G,H>=trace(GH),trace(·)表示矩阵的迹;
S42:对增广拉格朗日函数采用交替方向乘子法交替求解以下子问题得到:
s.t.PΩ(X)=X0
S43:通过求解以下最小二乘问题得到U:
S44:求解以下最小二乘问题得到V:
S45:通过以下公式迭代更新X:
Xk+1=X0+PΩc(Uk+1(Vk+1)T)
通过以下公式迭代更新Bi、Ci:
Ci的优化子问题表示为:
其中,σBi表示矩阵Bi的奇异值;σCi表示矩阵Ci的奇异值;||·||0表示矩阵的0范数。
每个RQiU和RQiV都为汉克尔矩阵,且都为秩—1矩阵,对于秩—1,信息大多在第一个奇异值中,采用硬阈值算子迭代可以更好地保留主要信息,去除干扰,在低采样率的情况下它的重建速度与精度效果都更好,说明硬阈值算子更适合本重建问题。矩阵采用0范数模型并用硬阈值算子进行迭代,比其它阈值缩减算子效果更好。
其中,每次迭代μ值逐渐增大,当达到迭代停止准则时,增大β并继续迭代,直至β达到设定值,迭代结束。β初值设为20-30,μ的初值为0.005-0.02,当完成一次迭代后,下一次迭代时μ的取值为上一次迭代时μ值的1.02-1.12倍,当达到迭代停止准则时,下一次迭代时β的取值为上次迭代时β的值的两倍,继续迭代直至β达到最大值222-242,迭代结束;迭代停止准则设定为在两次相邻迭代中所补全的矩阵Xk+1与Xk的误差小于阈值η,η设为10-7-10-5。
进一步地,步骤S43具体为:
R*为操作算子R的伴随算子,定义如下:
Qi*为Qi的伴随算子,定义如下:
其中,C表示复数域,M、N分别为矩阵的行数、列数;w1是一个M×1向量且它的第k个元素为矩阵RQiU第k条副对角线上的元素个数,w2是一个N×1向量且它的第k个元素为矩阵RQiV第k条副对角线上的元素个数;定义矩阵T1为r列w1、T2为r列w2组合而成;
所以,
转化为:
并通过以下公式通过更新U的每一行来迭代U:
U(a,:) k+1=Y1(a,:)(μkTa+β(Vk)Tconj(Vk))-1
其中,U(a,:) k+1表示矩阵U经过k+1次迭代后产生的矩阵Uk+1的第a行的元素;Y1(a,:)为矩阵Y1第a行的元素;Ta为对角矩阵,Ta主对角线元素为T1(a,:),T1(a,:)为矩阵T1第a行的元素。
更进一步地,步骤S44具体为:
得到:
并通过以下公式通过更新V的每一行来迭代V:
V(b,:) k+1=Y2(b,:)(μkTb+β(Uk+1)Tconj(Uk+1))-1
其中,V(b,:) k+1表示矩阵V经过k+1次迭代后产生的矩阵VK+1第n行的元素;Y2(b,:)为矩阵Y2第n行的元素;Tb为对角矩阵,Tb主对角线元素为T2(b,:),T2(b,:)为矩阵T2第b行的元素。
最后,步骤S5还包括:
当为三维磁共振信号时,对组成三维磁共振信号的所有二维磁共振信号分别做二维傅里叶变换,并对三维磁共振信号做三维傅里叶变换得到完整的三维磁共振波谱。
本发明还提供了一种磁共振波谱重建系统,如图3:本发明一种磁共振波谱重建系统的示意图所示,包括:获取模块1、奇异值分解模块2、低秩模型模块3、优化模块4和变换模块5;
获取模块1用于获取欠采样矩阵X0;
奇异值分解模块2用于对欠采样矩阵X0进行奇异值分解,得到X0=U0ΣV0 T,其中,U0为左奇异向量,V0为右奇异向量,Σ为奇异值;
取欠采样矩阵X0的前r个奇异值和与r个奇异值对应的左奇异向量和右奇异向量,且令Σ=ΣLΣR,U=U0ΣL,VT=ΣRV0 T则X0=UVT;
低秩模型模块3用于将U和V的每一列都转化为汉克尔矩阵,并构建磁共振波谱低秩模型:
s.t.PΩ(X)=X0
其中,r表示矩阵U和V的列数;R为将一维数据转成汉克尔矩阵的操作算子;操作算子Qi用来提取矩阵第i列;rank(·)表示矩阵的秩;表示矩阵的F-范数的平方;β为正则化参数,用于权衡前后两项的重要性,β>0;X表示被补全后的矩阵;集合Ω表示采样点的位置的集合;即PΩ与X的哈达玛积表示按照模板PΩ对X进行欠采样,表示为PΩ(X),这样矩阵X对应PΩ数字为1的位置的元素保持不变,对应模板数字0的位置的元素被置为0;
优化模块4用于引入中间变量Bi=RQiU、Ci=RQiV,将磁共振波谱低秩模型转化为最小秩优化问题:
s.t.Bi=RQiU,Ci=RQiV,PΩ(X)=X0
采用交替方向乘子法求解最小秩优化问题,得到Uk+1和Vk+1,并补全矩阵X:
X=Xk+1=PΩ(X0)+PΩc(Uk+1(Vk+1)T)
其中Xk+1为矩阵X经过k+1次迭代后的矩阵;Uk+1为矩阵U经过k+1次迭代后的矩阵;Vk+1为矩阵V经过k+1次迭代后的矩阵;PΩc是PΩ的补充,其值与PΩ相反,即PΩ为1的位置,PΩc在相对应的位置为0;PΩ为0的位置,PΩc在相对应的位置为1;
变换模块5用于当Xk+1为二维磁共振信号时,对Xk+1做二维傅里叶变换得到其二维频谱
现通过具体实施例来说明,欠采样矩阵X0的完整的二维磁共振波谱数据的大小为256×128,包含15个谱峰,此二维磁共振波谱如图4所示。本发明提出的磁共振波谱重建方法重建X0,得到重建数据X,对X做二维傅里叶变换得到其频谱,如图5所示。
实施例1:
如下述表1是本方法与现有的快速块Hankel重建方法、低秩Hankel方法、Hankel矩阵分解重建方法在不同的采样率下重建速度的对比,其中块Hankel方法将二维数据转变为一个更大的块Hankel矩阵再对此块Hankel矩阵进行重建,低秩Hankel方法与Hankel矩阵分解重建方法用于重建一维磁共振信号,应用于二维数据重建时即重建二维矩阵的每一列。Hankel矩阵分解重建方法是对一维数据的Hankel矩阵进行范德蒙德分解,而本方法则直接对二维数据本身进行基于范德蒙德进行奇异值分解,在重建二维或三维数据时,本方法更具优势。
表1:
采样率 | 20% | 25% | 30% | 35% | 40% |
本方法 | 6.5S | 6.9S | 7.5S | 11.2S | 12.0S |
快速块Hankel | 638.7S | 628.2S | 613.3S | 640.2S | 643.0S |
Hankel矩阵分解 | 482.2S | 476.4S | 470.3S | 478.1S | 415.8S |
低秩Hanekl | 6.9S | 7.5S | 9.0S | 11.2S | 13.3S |
可以看出在不同的采样率下,本方法的重建速度都是最快的,且相较快速块Hankel重建方法与Hankel矩阵分解重建方法有很大的提升。
表2是本方法与其他三种方法在不同采样率下重建精度的对比。
表2:
采样率 | 20% | 25% | 30% | 35% | 40% |
本方法 | 0.0601 | 0.0368 | 0.0529 | 0.0483 | 0.0370 |
快速块Hankel | 0.0544 | 0.0532 | 0.0554 | 0.0442 | 0.0341 |
Hankel矩阵分解 | 0.3379 | 0.1764 | 0.0790 | 0.0735 | 0.0526 |
低秩Hanekl | 0.4498 | 0.3213 | 0.2041 | 0.1349 | 0.0668 |
重建精度为重建结果与真实数据的相对误差。可以看出在不同的采样率下,本方法的重建精度与快速块Hankel方法相差很小,Hankel矩阵分解方法在低采样率时精度较差,而低秩Hankel方法需要更高的采样率。
实验结果如图6所示,是在20%的采样率下四种方法重建结果的对比,可以看出只有本方法与快速块Hankel方法有较高的重建质量。综上,在相同质量的重建结果下,本方法所需的重建时间较现有方法有很大提升,且重建要求的采样率也较低。
实施例2:
本实施例基于本申请的方法,采用不同的阈值缩减算子迭代重建,并比较它们的重建速度与精度,如下述表3:是不同的采样率下采用硬阈值算子与软阈值算子、拟软阈值算子在求解该问题时迭代次数的对比。软阈值算子对所有奇异值进行同等的缩减,硬阈值算子只缩减小于阈值的奇异值,而拟软阈值是软硬阈值结合的阈值缩减算子。
表3:
采样率 | 20% | 25% | 30% | 35% | 40% |
硬阈值 | 358 | 356 | 355 | 502 | 537 |
软阈值 | 729 | 727 | 723 | 715 | 710 |
拟软阈值 | 369 | 356 | 357 | 500 | 501 |
可以看出,硬阈值算子与拟软阈值算子的迭代次数都比软阈值少。
表4是在不同的采样率下,不同的阈值缩减算子重建精度的对比。
表4:
采样率 | 20% | 25% | 30% | 35% | 40% |
硬阈值 | 0.0579 | 0.0368 | 0.0529 | 0.0483 | 0.0370 |
软阈值 | 0.3242 | 0.2879 | 0.1978 | 0.1290 | 0.1003 |
拟软阈值 | 0.2350 | 0.1561 | 0.0540 | 0.0487 | 0.0446 |
可以看出,硬阈值算子在不同的采样率下重建精度都在0.1以下,而拟软阈值在较低采样率时精度不高,软阈值算子重建精度都较低。
实验结果如图7所示,在20%采样率下采用不同阈值算子的重建结果对比,可以看出硬阈值算子更加接近真实结果。综上,由于每个Hankel矩阵都应是秩-1矩阵,信息大多在第一个奇异值中,采用硬阈值算子迭代可以更好地保留主要信息,去除干扰,在低采样率的情况下它的重建速度与精度效果都更好,说明硬阈值算子更适合本重建问题。
需要理解的是,以上对本发明的具体实施例进行的描述只是为了说明本发明的技术路线和特点,其目的在于让本领域内的技术人员能够了解本发明的内容并据以实施,但本发明并不限于上述特定实施方式。凡是在本发明权利要求的范围内做出的各种变化或修饰,都应涵盖在本发明的保护范围内。
Claims (10)
1.一种磁共振波谱重建方法,其特征在于,所述磁共振波谱重建方法具体包括以下步骤:
S1:获取欠采样矩阵X0;
S2:对所述欠采样矩阵X0进行奇异值分解,得到X0=U0ΣV0 T,其中,U0为左奇异向量,V0为右奇异向量,Σ为奇异值;
取所述欠采样矩阵X0的前r个奇异值和与所述r个奇异值对应的左奇异向量和右奇异向量,令Σ=ΣLΣR,U=U0ΣL,VT=ΣRV0 T则X0=UVT;
S3:将所述U和所述V的每一列都转化为汉克尔矩阵,并构建磁共振波谱低秩模型:
s.t.PΩ(X)=X0
其中,r表示矩阵U和V的列数;R为将一维数据转成汉克尔矩阵的操作算子;操作算子Qi用来提取矩阵第i列;rank(·)表示矩阵的秩;表示矩阵的F-范数的平方;β为正则化参数,用于权衡前后两项的重要性,β>0;X表示被补全后的矩阵;集合Ω表示采样点的位置的集合;即PΩ与X的哈达玛积表示按照模板PΩ对所述X进行欠采样,表示为PΩ(X),所述矩阵X对应PΩ数字为1的位置的元素保持不变,对应模板数字0的位置的元素被置为0;
S4:引入中间变量Bi=RQiU、Ci=RQiV,将所述磁共振波谱低秩模型转化为最小秩优化问题:
s.t.Bi=RQiU,Ci=RQiV,PΩ(X)=X0
采用交替方向乘子法求解所述最小秩优化问题,得到所述Uk+1和所述Vk+1,并补全所述矩阵X:
X=Xk+1=PΩ(X0)+PΩc(Uk+1(Vk+1)T)
其中Xk+1为所述矩阵X经过k+1次迭代后的矩阵;Uk+1为所述矩阵U经过k+1次迭代后的矩阵;Vk+1为所述矩阵V经过k+1次迭代后的矩阵;PΩc是PΩ的补充,其值与PΩ相反,即PΩ为1的位置,PΩc在相对应的位置为0;PΩ为0的位置,PΩc在相对应的位置为1;
S5:当所述Xk+1为二维磁共振信号时,对所述Xk+1做二维傅里叶变换得到完整的二维磁共振波谱。
2.如权利要求1所述的磁共振波谱重建方法,其特征在于,在所述步骤S2中,所述r在实际的二维磁共振波谱谱峰数6倍之内。
3.如权利要求1所述的磁共振波谱重建方法,其特征在于,在所述步骤S4具体为:
S41:将所述最小秩优化问题转化为增广拉格朗日函数:
s.t.PΩ(X)=X0
其中,Di,Mi为拉格朗日乘子,μ为正则化参数,μ>0,<·,·>表示矩阵在希尔伯特空间的内积,<A,B>=trace(AB),trace(·)表示矩阵的迹;
S42:对所述增广拉格朗日函数采用交替方向乘子法交替求解以下子问题得到:
s.t.PΩ(X)=X0
S43:通过求解以下最小二乘问题得到U:
S44:求解以下最小二乘问题得到V:
S45:通过以下公式迭代更新X:
Xk+1=X0+PΩc(Uk+1(Vk+1)T)
通过以下公式迭代更新Bi、Ci:
5.如权利要求3所述的磁共振波谱重建方法,其特征在于,每次迭代所述μ值逐渐增大,当达到迭代停止准则时,增大所述β并继续迭代,直至所述β达到设定值,迭代结束。
6.如权利要求5所述的磁共振波谱重建方法,其特征在于,所述β初值设为20-30,所述μ初值为0.005-0.02,当完成一次迭代后,下一次迭代时所述μ的取值值为上一次迭代时μ值的1.02-1.12倍,当达到迭代停止准则时,下一次迭代时所述β的取值为上次迭代时所述β的值的两倍,继续迭代直至所述β达到最大值222-242,迭代结束;迭代停止准则设定为在两次相邻迭代中所补全的矩阵Xk+1与Xk的误差小于阈值η,η设为10-7-10-5。
7.如权利要求3所述的磁共振波谱重建方法,其特征在于,在所述步骤S43具体为:
R*为操作算子R的伴随算子,定义如下:
Qi*为Qi的伴随算子,定义如下:
其中,C表示复数域,M、N分别为矩阵的行数、列数;w1是一个M×1向量且它的第k个元素为所述矩阵RQiU第k条副对角线上的元素个数,w2是一个N×1向量且它的第k个元素为所述矩阵RQiV第k条副对角线上的元素个数;定义矩阵T1为r列w1、T2为r列w2组合而成;
所以,
并通过以下公式通过更新U的每一行来迭代U:
U(a,:) k+1=Y1(a,:)(μkTa+β(Vk)Tconj(Vk))-1
其中,U(a,:) k+1表示矩阵U经过k+1次迭代后产生的矩阵Uk+1的第a行的元素;Y1(a,:)为矩阵Y1第a行的元素;Ta为对角矩阵,Ta主对角线元素为T1(a,:),T1(a,:)为矩阵T1第a行的元素。
9.如权利要求1所述的磁共振波谱重建方法,其特征在于,所述步骤S5还包括:
当所述为三维磁共振信号时,对组成所述三维磁共振信号的所有二维磁共振信号分别做二维傅里叶变换,并对所述三维磁共振信号做三维傅里叶变换得到完整的三维磁共振波谱。
10.一种磁共振波谱重建系统,其特征在于,所述二维磁共振波谱重建系统包括:获取模块、奇异值分解模块、低秩模型模块、优化模块和变换模块;
所述获取模块用于获取欠采样矩阵X0;
所述奇异值分解模块用于对欠采样矩阵X0进行奇异值分解,得到X0=U0ΣV0 T,其中,U0为左奇异向量,V0为右奇异向量,Σ为奇异值;
取所述欠采样矩阵X0的前r个奇异值和与所述r个奇异值对应的左奇异向量和右奇异向量,且令Σ=ΣLΣR,U=U0ΣL,VT=ΣRV0 T则X0=UVT;
所述低秩模型模块用于将所述U和所述V的每一列都转化为汉克尔矩阵,并构建磁共振波谱低秩模型:
s.t.PΩ(X)=X0
其中,r表示矩阵U和V的列数;R为将一维数据转成汉克尔矩阵的操作算子;操作算子Qi用来提取矩阵第i列;rank(·)表示矩阵的秩;表示矩阵的F-范数的平方;β为正则化参数,用于权衡前后两项的重要性,β>0;X表示被补全后的矩阵;集合Ω表示采样点的位置的集合;即PΩ与X的哈达玛积表示按照模板PΩ对所述X进行欠采样,表示为PΩ(X),这样所述矩阵X对应PΩ数字为1的位置的元素保持不变,对应模板数字0的位置的元素被置为0;
所述优化模块用于引入中间变量Bi=RQiU、Ci=RQiV,将所述磁共振波谱低秩模型转化为最小秩优化问题:
s.t.Bi=RQiU,Ci=RQiV,PΩ(X)=X0
采用交替方向乘子法求解所述最小秩优化问题,得到所述Uk+1和所述Vk+1,并补全所述矩阵X:
X=Xk+1=PΩ(X0)+PΩc(Uk+1(Vk+1)T)
其中Xk+1为所述矩阵X经过k+1次迭代后的矩阵;Uk+1为所述矩阵U经过k+1次迭代后的矩阵;Vk+1为所述矩阵V经过k+1次迭代后的矩阵;PΩc是PΩ的补充,其值与PΩ相反,即PΩ为1的位置,PΩc在相对应的位置为0;PΩ为0的位置,PΩc在相对应的位置为1;
所述变换模块用于当所述Xk+1为二维磁共振信号时,对所述Xk+1做二维傅里叶变换得到其二维频谱。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910990393.8A CN110658484B (zh) | 2019-10-17 | 2019-10-17 | 一种磁共振波谱重建方法及系统 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201910990393.8A CN110658484B (zh) | 2019-10-17 | 2019-10-17 | 一种磁共振波谱重建方法及系统 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN110658484A true CN110658484A (zh) | 2020-01-07 |
CN110658484B CN110658484B (zh) | 2021-07-20 |
Family
ID=69041146
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201910990393.8A Active CN110658484B (zh) | 2019-10-17 | 2019-10-17 | 一种磁共振波谱重建方法及系统 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN110658484B (zh) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111538944A (zh) * | 2020-04-13 | 2020-08-14 | 厦门理工学院 | 一种基于子空间的磁共振波谱快速重建方法 |
CN113180636A (zh) * | 2021-04-29 | 2021-07-30 | 杭州微影医疗科技有限公司 | 干扰消除方法、介质及设备 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103654789A (zh) * | 2013-12-10 | 2014-03-26 | 深圳先进技术研究院 | 磁共振快速参数成像方法和系统 |
CN105808869A (zh) * | 2016-03-16 | 2016-07-27 | 厦门理工学院 | 一种基于块Hankel矩阵的磁共振波谱重建方法 |
CN106649201A (zh) * | 2016-09-27 | 2017-05-10 | 厦门大学 | 一种基于指数信号的范德蒙分解的数据补全方法 |
CN106646303A (zh) * | 2016-11-17 | 2017-05-10 | 厦门理工学院 | 一种欠采样磁共振波谱的快速重建方法 |
CN109165432A (zh) * | 2018-08-09 | 2019-01-08 | 厦门理工学院 | 一种基于部分奇异值和的磁共振波谱重建方法 |
CN109191540A (zh) * | 2018-07-24 | 2019-01-11 | 厦门理工学院 | 一种基于截断核范数的磁共振波谱重建方法 |
CN109903259A (zh) * | 2019-01-25 | 2019-06-18 | 厦门大学 | 一种基于深度学习的磁共振波谱重建方法 |
-
2019
- 2019-10-17 CN CN201910990393.8A patent/CN110658484B/zh active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103654789A (zh) * | 2013-12-10 | 2014-03-26 | 深圳先进技术研究院 | 磁共振快速参数成像方法和系统 |
CN105808869A (zh) * | 2016-03-16 | 2016-07-27 | 厦门理工学院 | 一种基于块Hankel矩阵的磁共振波谱重建方法 |
CN106649201A (zh) * | 2016-09-27 | 2017-05-10 | 厦门大学 | 一种基于指数信号的范德蒙分解的数据补全方法 |
CN106646303A (zh) * | 2016-11-17 | 2017-05-10 | 厦门理工学院 | 一种欠采样磁共振波谱的快速重建方法 |
CN109191540A (zh) * | 2018-07-24 | 2019-01-11 | 厦门理工学院 | 一种基于截断核范数的磁共振波谱重建方法 |
CN109165432A (zh) * | 2018-08-09 | 2019-01-08 | 厦门理工学院 | 一种基于部分奇异值和的磁共振波谱重建方法 |
CN109903259A (zh) * | 2019-01-25 | 2019-06-18 | 厦门大学 | 一种基于深度学习的磁共振波谱重建方法 |
Non-Patent Citations (1)
Title |
---|
GUO, DI 等: "Improved Reconstruction of Low Intensity Magnetic Resonance Spectroscopy With Weighted Low Rank Hankel Matrix Completion", 《IEEE ACCESS》 * |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111538944A (zh) * | 2020-04-13 | 2020-08-14 | 厦门理工学院 | 一种基于子空间的磁共振波谱快速重建方法 |
CN111538944B (zh) * | 2020-04-13 | 2022-06-24 | 厦门理工学院 | 一种基于子空间的磁共振波谱快速重建方法 |
CN113180636A (zh) * | 2021-04-29 | 2021-07-30 | 杭州微影医疗科技有限公司 | 干扰消除方法、介质及设备 |
Also Published As
Publication number | Publication date |
---|---|
CN110658484B (zh) | 2021-07-20 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
WO2020151355A1 (zh) | 一种基于深度学习的磁共振波谱重建方法 | |
CN106646303B (zh) | 一种欠采样磁共振波谱的快速重建方法 | |
CN110658484B (zh) | 一种磁共振波谱重建方法及系统 | |
CN105808869A (zh) | 一种基于块Hankel矩阵的磁共振波谱重建方法 | |
CN109615675B (zh) | 一种多通道磁共振成像的图像重建方法 | |
CN111324861B (zh) | 一种基于矩阵分解的深度学习磁共振波谱重建方法 | |
CN104739410B (zh) | 一种磁共振图像的迭代重建方法 | |
CN107423543B (zh) | 一种超复数磁共振波谱的快速重建方法 | |
Lyashenko et al. | Ideology of Image Processing in Infocommunication Systems | |
CN108828482B (zh) | 结合稀疏和低秩特性的欠采样磁共振扩散谱的重建方法 | |
CN111783631B (zh) | 一种基于稀疏表示的深度学习磁共振波谱重建方法 | |
CN114820849A (zh) | 基于深度学习的磁共振cest图像重建方法、装置及设备 | |
Wang et al. | Efficient approximation of Jacobian matrices involving a non-uniform fast Fourier transform (NUFFT) | |
CN110598579B (zh) | 一种基于深度学习的超复数磁共振波谱重建方法 | |
CN109165432A (zh) | 一种基于部分奇异值和的磁共振波谱重建方法 | |
CN108537738A (zh) | 一种矩阵补全方法 | |
CN111538944B (zh) | 一种基于子空间的磁共振波谱快速重建方法 | |
CN109191540B (zh) | 一种基于截断核范数的磁共振波谱重建方法 | |
CN113143243B (zh) | 一种基于谱分解的深度学习磁共振波谱重建方法 | |
CN115471580A (zh) | 一种物理智能高清磁共振扩散成像方法 | |
CN108920423B (zh) | 一种高保真谱重建方法 | |
CN109033980B (zh) | 基于增量局部残差最小二乘法的高光谱图像Gabor特征分类方法 | |
CN113034639B (zh) | 一种基于可分离汉克尔矩阵的磁共振成像的图像重建方法 | |
CN114972562B (zh) | 联合线圈灵敏度估计与图像重建的快速磁共振成像方法 | |
CN117078785B (zh) | 一种快速非笛卡尔磁共振智能成像方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |