CN110658484B - 一种磁共振波谱重建方法及系统 - Google Patents

Abstract

本发明涉及一种磁共振波谱重建方法及系统，该方法包括：首先获取欠采样矩阵X₀；然后对欠采样矩阵X₀进行奇异值分解，得到X₀＝UV^T，再将U和V的每一列都转化为汉克尔矩阵，并构建磁共振波谱低秩模型；接着引入中间变量将磁共振波谱低秩模型转化为最小秩优化问题；采用交替方向乘子法求解最小秩优化问题，得到U^k+1和V^k+1,并补全矩阵X：最后当X^k+1为二维磁共振信号时，对X^k+1做二维傅里叶变换得到完整的二维磁共振波谱。本发明利用磁共振信号可以进行范德蒙德分解的性质，对欠采样矩阵X₀进行奇异值分解得矩阵U和V，通过采用交替方向乘子法分别求解U和V，完成对U和V的重建，从而实现对欠采样矩阵X₀的补全，重建速度快且精度高。

Description

一种磁共振波谱重建方法及系统

技术领域

本发明涉及属于磁共振波谱重建领域，具体涉及一种磁共振波谱重建方法及系统。

背景技术

核磁共振波谱法是研究原子核对射频辐射的吸收，它是对各种有机和无机物的成分、结构进行定性分析的最强有力的工具之一，有时亦可进行定量分析。核磁共振技术可以提供分子的化学结构和分子动力学的信息，已成为分子结构解析以及物质理化性质表征的常规技术手段，它可以确定几乎所有常见官能团的环境，直观性强，特别是碳谱能直接反映出分子的骨架，谱图解释较为容易，在物理、化学、生物、医药、食品等领域得到广泛应用，在化学中更是常规分析不可少的手段。

在实际应用中，由于较高的实验成本和硬件限制，以及其它不可避免的原因，我们很难对二维或三维磁共振自由衰减信号实现完全采样，因此对磁共振信号进行欠采样并通过欠采样得到的数据重建出完整的磁共振信号显得尤为重要。不同的重建方法在重建精度、最小欠采样率和重建速度上有较大差异，而目前已有的方法都需要较长的重建时间，尤其重建三维磁共振波谱时，所需的时间成比例增加。为此，提供一种重建速度更快且有较高重建精度的新的重建方法成为当前需要解决的问题。

发明内容

(一)要解决的技术问题

本发明提供了一种磁共振波谱重建方法，旨在提供一种重建速度更快且有较高重建精度的新的重建方法。

(二)技术方案

为了解决上述问题，本发明提供一种磁共振波谱重建方法，包括以下步骤：

S1：获取欠采样矩阵X₀；

S2：对所述欠采样矩阵X₀进行奇异值分解，得到X₀＝U₀ΣV₀ ^T，其中，U₀为左奇异向量，V₀为右奇异向量，Σ为奇异值；

取所述欠采样矩阵X₀的前r个奇异值和与所述r个奇异值对应的左奇异向量和右奇异向量，令Σ＝Σ_LΣ_R，U＝U₀Σ_L，V^T＝Σ_RV₀ ^T则X₀＝UV^T；

S3：将所述U和所述V的每一列都转化为汉克尔矩阵，并构建磁共振波谱低秩模型：

s.t.P_Ω(X)＝X₀

其中，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)＝X₀

采用交替方向乘子法求解所述最小秩优化问题，得到所述U^k+1和所述V^k+1,并补全所述矩阵X：

X＝X^k+1＝P_Ω(X₀)+P_Ωc(U^k+1(V^k+1)^T)

其中X^k+1为所述矩阵X经过k+1次迭代后的矩阵；U^k+1为所述矩阵U经过k+1次迭代后的矩阵；V^k+1为所述矩阵V经过k+1次迭代后的矩阵；P_Ωc是P_Ω的补充，其值与P_Ω相反，即P_Ω为1的位置，P_Ωc在相对应的位置为0；P_Ω为0的位置，P_Ωc在相对应的位置为1；

S5：当所述X^k+1为二维磁共振信号时，对所述X^k+1做二维傅里叶变换得到完整的二维磁共振波谱。

优选地，在所述步骤S2中，所述r在实际的二维磁共振波谱谱峰数6倍之内。

优选地，在所述步骤S4具体为：

S41：将所述最小秩优化问题转化为增广拉格朗日函数：

s.t.P_Ω(X)＝X₀

其中，Di,Mi为拉格朗日乘子，μ为正则化参数，μ>0,<·，·>表示矩阵在希尔伯特空间的内积，<A,B>＝trace(AB),trace(·)表示矩阵的迹；

S42：对所述增广拉格朗日函数采用交替方向乘子法交替求解以下子问题得到：

s.t.P_Ω(X)＝X₀

S43：通过求解以下最小二乘问题得到U：

S44：求解以下最小二乘问题得到V：

S45：通过以下公式迭代更新X：

X^k+1＝X₀+P_Ωc(U^k+1(V^k+1)^T)

通过以下公式迭代更新B_i、C_i：

其中，其中k为迭代次数；R*为操作算子R的伴随算子；Qi*为Qi的伴随算子；

为奇异值缩减算子，

为所述B_i经过k+1次迭代后的结果，

为所述C_i经过k+1次迭代后的结果,μ^k为所述μ经过k次迭代后的结果。

优选地，所述奇异值缩减算子

为硬阈值算子，矩阵的秩为奇异值的0范数，在所述步骤S42中，B_i的优化子问题表示为：

所述C_i的优化子问题表示为：

其中，σ_Bi表示矩阵Bi的奇异值；σ_Ci表示矩阵Ci的奇异值；||·||₀表示矩阵的0范数。

优选地，每次迭代所述μ值逐渐增大，当达到迭代停止准则时，增大所述β并继续迭代，直至所述β达到设定值，迭代结束。

优选地，所述β初值设为20-30，所述μ初值为0.005-0.02，当完成一次迭代后，下一次迭代时所述μ的取值值为上一次迭代时μ值的1.02-1.12倍，当达到迭代停止准则时，下一次迭代时所述β的取值为上次迭代时所述β的值的两倍，继续迭代直至所述β达到最大值2²²-2⁴²，迭代结束；迭代停止准则设定为在两次相邻迭代中所补全的矩阵X^k+1与X^k的误差小于阈值η，η设为10^-7-10^-5。

优选地，在所述步骤S43具体为：

R^*为操作算子R的伴随算子，定义如下：

R^*D将一个矩阵D转化为一个向量，

表示该向量的第l个元素为矩阵D的第l条副对角线上的元素之和；

Qi*为Qi的伴随算子，定义如下：

d为向量，

表示将向量d的元素放入矩阵[Qi^*d]的p列中，x为向量d的元素；

所以，

其中，C表示复数域，M、N分别为矩阵的行数、列数；w₁是一个M×1向量且它的第k个元素为所述矩阵RQiU第k条副对角线上的元素个数，w₂是一个N×1向量且它的第k个元素为所述矩阵RQiV第k条副对角线上的元素个数；定义矩阵T1为r列w₁、T2为r列w₂组合而成；

令

所以，

转化为：

并通过以下公式通过更新U的每一行来迭代U：

U_(a,:) ^k+1＝Y1_(a,:)(μ^kT_a+β(V^k)^Tconj(V^k))^-1

其中，U_(a,:) ^k+1表示矩阵U经过k+1次迭代后产生的矩阵U^k+1的第a行的元素；Y1_(a,:)为矩阵Y1第a行的元素；T_a为对角矩阵，T_a主对角线元素为T1_(a,:)，T1_(a,:)为矩阵T1第a行的元素。

优选地，所述步骤S44具体为：

令

得到：

转化为

并通过以下公式通过更新V的每一行来迭代V：

V_(b,:) ^k+1＝Y2_(b,:)(μ^kT_b+β(U^k+1)^Tconj(U^k+1))^-1

其中，V_(b,:) ^k+1表示矩阵V经过k+1次迭代后产生的矩阵V^k+1第b行的元素；Y2_(b,:)为矩阵Y2第b行的元素；T_b为对角矩阵，T_b主对角线元素为T2_(b,:)，T2_(b,:)为矩阵T2第b行的元素。

优选地，所述步骤S5还包括：

当所述X^k+1为三维磁共振信号时，对组成所述三维磁共振信号的所有二维磁共振信号分别做二维傅里叶变换，并对所述三维磁共振信号做三维傅里叶变换得到完整的三维磁共振波谱。

优选地，本发明还提供了一种磁共振波谱重建系统，包括：获取模块、奇异值分解模块、低秩模型模块、优化模块和变换模块；

所述获取模块用于获取欠采样矩阵X₀；

所述奇异值分解模块用于对欠采样矩阵X₀进行奇异值分解，得到X₀＝U₀ΣV₀ ^T，其中，U₀为左奇异向量，V₀为右奇异向量，Σ为奇异值；

取所述欠采样矩阵X₀的前r个奇异值和与所述r个奇异值对应的左奇异向量和右奇异向量，且令Σ＝Σ_LΣ_R，U＝U₀Σ_L，V^T＝Σ_RV₀ ^T则X₀＝UV^T；

所述低秩模型模块用于将所述U和所述V的每一列都转化为汉克尔矩阵，并构建磁共振波谱低秩模型：

s.t.P_Ω(X)＝X₀

其中，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)＝X₀

采用交替方向乘子法求解所述最小秩优化问题，得到所述U^k+1和所述V^k+1,并补全所述矩阵X：

X＝X^k+1＝P_Ω(X₀)+P_Ωc(U^k+1(V^k+1)^T)

其中X^k+1为所述矩阵X经过k+1次迭代后的矩阵；U^k+1为所述矩阵U经过k+1次迭代后的矩阵；V^k+1为所述矩阵V经过k+1次迭代后的矩阵；P_Ωc是P_Ω的补充，其值与P_Ω相反，即P_Ω为1的位置，P_Ωc在相对应的位置为0；P_Ω为0的位置，P_Ωc在相对应的位置为1；

所述变换模块用于当所述X^k+1为二维磁共振信号时，对所述X^k+1做二维傅里叶变换得到其二维频谱。

(三)有益效果

本发明利用将欠采样矩阵进行奇异值分解，再将分解得到的矩阵的每一列都转化为汉克尔矩阵，并建立磁共振波谱低秩模型，通过采用交替乘子的方法获得重建后的矩阵，最后对重建后的矩阵进行傅里叶变换得到完整的磁共振波谱。这种重建方法重建速度更快且有较高重建精度。

附图说明

图1为本发明一种磁共振波谱重建方法的流程图；

图2为本发明中求解最小秩优化问题的流程图；

图3为本发明一种磁共振波谱重建系统的示意图；

图4为本发明实施例中大小为256×128的全采样二维磁共振波谱；

图5为使用本发明的方法对欠采样矩阵补全重建出的磁共振波谱；

图6为本发明与其它方法在20％采样率下重建结果的对比；

图7为采用不同阈值缩减算子在20％采样率下重建结果的对比。

【附图标记说明】

1：获取模块；2：奇异值分解模块；3：低秩模型模块；4：优化模块；5：变换模块。

具体实施方式

为了更好的解释本发明，以便于理解，下面结合附图，通过具体实施方式，对本发明作详细描述。

需要说明，本发明实施例中所有方向性指示(诸如上、下、左、右、前、后……)仅用于解释在某一特定姿态(如附图所示)下各部件之间的相对位置关系、运动情况等，如果该特定姿态发生改变时，则该方向性指示也相应地随之改变。

另外，在本发明中如涉及“第一”、“第二”等的描述仅用于描述目的，而不能理解为指示或暗示其相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本发明的描述中，“多个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。

在本发明中，除非另有明确的规定和限定，术语“连接”、“固定”等应做广义理解，例如，“固定”可以是固定连接，也可以是可拆卸连接，或成一体；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通或两个元件的相互作用关系，除非另有明确的限定。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。

本发明提供一种磁共振波谱重建方法，如图1：本发明一种磁共振波谱重建方法的流程图所示，具体包括以下步骤：

S1：获取欠采样矩阵X₀。X₀为一个M×N的矩阵，M、N分别表示矩阵X₀的行数和列数。

S2：对欠采样矩阵X₀进行奇异值分解，得到X₀＝U₀ΣV₀ ^T，其中，U₀为左奇异向量，V₀为右奇异向量，Σ为奇异值。

取欠采样矩阵X₀的前r个奇异值和与r个奇异值对应的左奇异向量和右奇异向量，令Σ＝Σ_LΣ_R，U＝U₀Σ_L，V^T＝Σ_RV₀ ^T则X₀＝UV^T。U为M×r矩阵，V为N×r矩阵。

一维磁共振信号可以表示为几个复指数衰减信号的叠加，而二维磁共振信号即是几个两两相乘的复指数衰减信号的叠加，由二维磁共振信号构成的矩阵P中第m行，n列的元素表示为：

其中，y_i＝exp(j2πf_1i-τ_1i),z_i＝exp(j2πf_2i-τ_2i)，f表示归一化频率，τ为衰减因子，b表示该二维磁共振波谱有g个谱峰；一个完整的二维磁共振信号矩阵P可以被分解为两个范德蒙德矩阵,如下所示：

D＝diag[d₁,d₂，…d_g]

P＝YDZ^T，将P分为D＝D_LD_R两个对角矩阵，分别乘入Y与Z^T,设U₁＝YD_L,V₁ ^T＝D_RZ^T,则X＝U₁V₁ ^T,U₁、V₁为范德蒙德矩阵。

欠采样矩阵X₀也是对磁共振信号欠采样而得到的，在对矩阵X₀进行奇异值分解后得到U、V，如果我们使U、V能够逼近范德蒙德矩阵U₁、V₁，那么就可以重建得到完整的磁共振信号，本申请利用磁共振信号能够分解为两个范德蒙德矩阵的性质，将补全矩阵X₀的问题转为使U、V逼近U₁、V₁，从而使得补全的方法更加准确快速。

S3：将U和V的每一列都转化为汉克尔矩阵，并构建磁共振波谱低秩模型：

s.t.P_Ω(X)＝X₀

其中，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)＝X₀

采用交替方向乘子法求解最小秩优化问题，得到U^k+1和V^k+1,并补全矩阵X：

X＝X^k+1＝P_Ω(X₀)+P_Ωc(U^k+1(V^k+1)^T)

其中X^k+1为矩阵X经过k+1次迭代后的矩阵；U^k+1为矩阵U经过k+1次迭代后的矩阵；V^k+1为矩阵V经过k+1次迭代后的矩阵；P_Ωc是P_Ω的补充，其值与P_Ω相反，即P_Ω为1的位置，P_Ωc在相对应的位置为0；P_Ω为0的位置，P_Ωc在相对应的位置为1；

S5：当X^k+1为二维磁共振信号时，对X^k+1做二维傅里叶变换得到完整的二维磁共振波谱。

进一步地，在步骤S2中，r在实际的二维磁共振波谱谱峰数6倍之内。在一些实际情况中不确定磁共振信号的谱峰数，故需预设谱峰r，预设数在实际谱峰数6倍之内都可以成功重建，即满足r<6g。

如图2：本发明中求解最小秩优化问题的流程图所示，在优选的实施例中，步骤S4具体为：

S41：将最小秩优化问题转化为增广拉格朗日函数：

s.t.P_Ω(X)＝X₀

其中，Di,Mi为拉格朗日乘子，μ为正则化参数，μ>0,<·，·>表示矩阵在希尔伯特空间的内积，<G,H>＝trace(GH),trace(·)表示矩阵的迹；

S42：对增广拉格朗日函数采用交替方向乘子法交替求解以下子问题得到：

s.t.P_Ω(X)＝X₀

S43：通过求解以下最小二乘问题得到U：

S44：求解以下最小二乘问题得到V：

S45：通过以下公式迭代更新X：

X^k+1＝X₀+P_Ωc(U^k+1(V^k+1)^T)

通过以下公式迭代更新B_i、C_i：

其中，其中k为迭代次数；R*为操作算子R的伴随算子；Qi*为Qi的伴随算子；

为奇异值缩减算子，

为B_i经过k+1次迭代后的结果，

为C_i经过k+1次迭代后的结果,μ^k为μ经过k次迭代后的结果。

在更优选的实施方式中，奇异值缩减算子

为硬阈值算子，矩阵的秩为奇异值的0范数，在步骤S42中，B_i的优化子问题表示为：

C_i的优化子问题表示为：

其中，σ_Bi表示矩阵Bi的奇异值；σ_Ci表示矩阵Ci的奇异值；||·||₀表示矩阵的0范数。

每个RQiU和RQiV都为汉克尔矩阵，且都为秩—1矩阵，对于秩—1，信息大多在第一个奇异值中，采用硬阈值算子迭代可以更好地保留主要信息，去除干扰，在低采样率的情况下它的重建速度与精度效果都更好，说明硬阈值算子更适合本重建问题。矩阵采用0范数模型并用硬阈值算子进行迭代，比其它阈值缩减算子效果更好。

其中，每次迭代μ值逐渐增大，当达到迭代停止准则时，增大β并继续迭代，直至β达到设定值，迭代结束。β初值设为20-30，μ的初值为0.005-0.02，当完成一次迭代后，下一次迭代时μ的取值为上一次迭代时μ值的1.02-1.12倍，当达到迭代停止准则时，下一次迭代时β的取值为上次迭代时β的值的两倍，继续迭代直至β达到最大值2²²-2⁴²，迭代结束；迭代停止准则设定为在两次相邻迭代中所补全的矩阵X^k+1与X^k的误差小于阈值η，η设为10^-7-10^-5。

进一步地，步骤S43具体为：

R^*为操作算子R的伴随算子，定义如下：

R^*D将一个矩阵D转化为一个向量，

表示该向量的第l个元素为矩阵D的第l条副对角线上的元素之和；

Qi*为Qi的伴随算子，定义如下：

d为向量，

表示将向量d的元素放入矩阵[Qi^*d]的p列中，x为向量d的元素；

所以，

其中，C表示复数域，M、N分别为矩阵的行数、列数；w₁是一个M×1向量且它的第k个元素为矩阵RQiU第k条副对角线上的元素个数，w₂是一个N×1向量且它的第k个元素为矩阵RQiV第k条副对角线上的元素个数；定义矩阵T1为r列w₁、T2为r列w₂组合而成；

令

所以，

转化为：

并通过以下公式通过更新U的每一行来迭代U：

U_(a,:) ^k+1＝Y1_(a,:)(μ^kT_a+β(V^k)^Tconj(V^k))^-1

其中，U_(a,:) ^k+1表示矩阵U经过k+1次迭代后产生的矩阵U^k+1的第a行的元素；Y1_(a,:)为矩阵Y1第a行的元素；T_a为对角矩阵，T_a主对角线元素为T1_(a,:)，T1_(a,:)为矩阵T1第a行的元素。

更进一步地，步骤S44具体为：

令

得到：

转化为

并通过以下公式通过更新V的每一行来迭代V：

V_(b,:) ^k+1＝Y2_(b,:)(μ^kT_b+β(U^k+1)^Tconj(U^k+1))^-1

其中，V_(b,:) ^k+1表示矩阵V经过k+1次迭代后产生的矩阵V^k+1第b行的元素；Y2_(b,:)为矩阵Y2第b行的元素；T_b为对角矩阵，T_b主对角线元素为T2_(b,:)，T2_(b,:)为矩阵T2第b行的元素。

最后，步骤S5还包括：

当X^k+1为三维磁共振信号时，对组成三维磁共振信号的所有二维磁共振信号分别做二维傅里叶变换，并对三维磁共振信号做三维傅里叶变换得到完整的三维磁共振波谱。

本发明还提供了一种磁共振波谱重建系统，如图3：本发明一种磁共振波谱重建系统的示意图所示，包括：获取模块1、奇异值分解模块2、低秩模型模块3、优化模块4和变换模块5；

获取模块1用于获取欠采样矩阵X₀；

奇异值分解模块2用于对欠采样矩阵X₀进行奇异值分解，得到X₀＝U₀ΣV₀ ^T，其中，U₀为左奇异向量，V₀为右奇异向量，Σ为奇异值；

取欠采样矩阵X₀的前r个奇异值和与r个奇异值对应的左奇异向量和右奇异向量，且令Σ＝Σ_LΣ_R，U＝U₀Σ_L，V^T＝Σ_RV₀ ^T则X₀＝UV^T；

低秩模型模块3用于将U和V的每一列都转化为汉克尔矩阵，并构建磁共振波谱低秩模型：

s.t.P_Ω(X)＝X₀

其中，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)＝X₀

采用交替方向乘子法求解最小秩优化问题，得到U^k+1和V^k+1,并补全矩阵X：

X＝X^k+1＝P_Ω(X₀)+P_Ωc(U^k+1(V^k+1)^T)

其中X^k+1为矩阵X经过k+1次迭代后的矩阵；U^k+1为矩阵U经过k+1次迭代后的矩阵；V^k+1为矩阵V经过k+1次迭代后的矩阵；P_Ωc是P_Ω的补充，其值与P_Ω相反，即P_Ω为1的位置，P_Ωc在相对应的位置为0；P_Ω为0的位置，P_Ωc在相对应的位置为1；

变换模块5用于当X_k+1为二维磁共振信号时，对X_k+1做二维傅里叶变换得到完整的二维磁共振波谱。

现通过具体实施例来说明，欠采样矩阵X₀的完整的二维磁共振波谱数据的大小为256×128，包含15个谱峰，此二维磁共振波谱如图4所示。本发明提出的磁共振波谱重建方法重建X₀，得到重建数据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 (8)

1.一种磁共振波谱重建方法，其特征在于，所述磁共振波谱重建方法具体包括以下步骤：

S1：获取欠采样矩阵X₀；

S2：对所述欠采样矩阵X₀进行奇异值分解，得到X₀＝U₀ΣV₀ ^T，其中，U₀为左奇异向量，V₀为右奇异向量，Σ为奇异值；

取所述欠采样矩阵X₀的前r个奇异值和与所述r个奇异值对应的左奇异向量和右奇异向量，令Σ＝Σ_LΣ_R，U＝U₀Σ_L，V^T＝Σ_RV₀ ^T则X₀＝UV^T；

S3：将所述U和所述V的每一列都转化为汉克尔矩阵，并构建磁共振波谱低秩模型：

s.t.P_Ω(X)＝X₀

其中，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)＝X₀

采用交替方向乘子法求解所述最小秩优化问题，得到U^k+1和V^k+1,并补全所述矩阵X：

X＝X^k+1＝P_Ω(X₀)+P_Ωc(U^k+1(V^k+1)^T)

其中X^k+1为所述矩阵X经过k+1次迭代后的矩阵；U^k+1为所述矩阵U经过k+1次迭代后的矩阵；V^k+1为所述矩阵V经过k+1次迭代后的矩阵；P_Ωc是P_Ω的补充，其值与P_Ω相反，即P_Ω为1的位置，P_Ωc在相对应的位置为0；P_Ω为0的位置，P_Ωc在相对应的位置为1；

所述步骤S4具体为：

S41：将所述最小秩优化问题转化为增广拉格朗日函数：

s.t.P_Ω(X)＝X₀

其中，Di,Mi为拉格朗日乘子，μ为正则化参数，μ>0,<·，·>表示矩阵在希尔伯特空间的内积，<A,B>＝trace(AB),trace(·)表示矩阵的迹；

S42：对所述增广拉格朗日函数采用交替方向乘子法交替求解以下子问题得到：

s.t.P_Ω(X)＝X₀

S43：通过求解以下最小二乘问题得到U：

S44：求解以下最小二乘问题得到V：

S45：通过以下公式迭代更新X：

X^k+1＝X₀+P_Ωc(U^k+1(V^k+1)^T)

通过以下公式迭代更新B_i、C_i：

其中，其中k为迭代次数；R*为操作算子R的伴随算子；Qi*为Qi的伴随算子；

为奇异值缩减算子，

为所述B_i经过k+1次迭代后的结果，

为所述C_i经过k+1次迭代后的结果,μ^k为所述μ经过k次迭代后的结果；

所述奇异值缩减算子

为硬阈值算子，矩阵的秩为奇异值的0范数，在所述步骤S42中，B_i的优化子问题表示为：

所述C_i的优化子问题表示为：

其中，σ_Bi表示矩阵Bi的奇异值；σ_Ci表示矩阵Ci的奇异值；||·||₀表示矩阵的0范数；

S5：当所述X^k+1为二维磁共振信号时，对所述X^k+1做二维傅里叶变换得到完整的二维磁共振波谱。

2.如权利要求1所述的磁共振波谱重建方法，其特征在于，在所述步骤S2中，所述r在实际的二维磁共振波谱谱峰数6倍之内。

3.如权利要求1所述的磁共振波谱重建方法，其特征在于，每次迭代所述μ值逐渐增大，当达到迭代停止准则时，增大所述β并继续迭代，直至所述β达到设定值，迭代结束。

4.如权利要求3所述的磁共振波谱重建方法，其特征在于，所述β初值设为20-30，所述μ初值为0.005-0.02，当完成一次迭代后，下一次迭代时所述μ的取值为上一次迭代时μ值的1.02-1.12倍，当达到迭代停止准则时，下一次迭代时所述β的取值为上次迭代时所述β的值的两倍，继续迭代直至所述β达到最大值2²²-2⁴²，迭代结束；迭代停止准则设定为在两次相邻迭代中所补全的矩阵X^k+1与X^k的误差小于阈值η，η设为10^-7-10^-5。

5.如权利要求1所述的磁共振波谱重建方法，其特征在于，所述步骤S43具体为：

R^*为操作算子R的伴随算子，定义如下：

R^*D将一个矩阵D转化为一个向量，

表示该向量的第l个元素为矩阵D的第l条副对角线上的元素之和；

Qi*为Qi的伴随算子，定义如下：

d为向量，

表示将向量d的元素放入矩阵[Qi^*d]的p列中，x为向量d的元素；

所以，

其中，C表示复数域，M、N分别为矩阵的行数、列数；w₁是一个M×1向量且它的第k个元素为所述矩阵RQiU第k条副对角线上的元素个数，w₂是一个N×1向量且它的第k个元素为所述矩阵RQiV第k条副对角线上的元素个数；定义矩阵T1为r列w₁、T2为r列w₂组合而成；

令

所以，

转化为：

并通过以下公式通过更新U的每一行来迭代U：

U_(a,:) ^k+1＝Y1_(a,:)(μ^kT_a+β(V^k)^Tconj(V^k))^-1

其中，U_(a,:) ^k+1表示矩阵U经过k+1次迭代后产生的矩阵U^k+1的第a行的元素；Y1_(a,:)为矩阵Y1第a行的元素；T_a为对角矩阵，T_a主对角线元素为T1_(a,:)，T1_(a,:)为矩阵T1第a行的元素。

6.如权利要求5所述的磁共振波谱重建方法，其特征在于，所述步骤S44具体为：

令

得到：

转化为

并通过以下公式通过更新V的每一行来迭代V：

V_(b,:) ^k+1＝Y2_(b,:)(μ^kT_b+β(U^k+1)^Tconj(U^k+1))^-1

其中，V_(b,:) ^k+1表示矩阵V经过k+1次迭代后产生的矩阵V^k+1第b行的元素；Y2_(b,:)为矩阵Y2第b行的元素；T_b为对角矩阵，T_b主对角线元素为T2_(b,:)，T2_(b,:)为矩阵T2第b行的元素。

7.如权利要求1所述的磁共振波谱重建方法，其特征在于，所述步骤S5还包括：

当所述X^k+1为三维磁共振信号时，对组成所述三维磁共振信号的所有二维磁共振信号分别做二维傅里叶变换，并对所述三维磁共振信号做三维傅里叶变换得到完整的三维磁共振波谱。

8.一种磁共振波谱重建系统，其特征在于，所述磁共振波谱重建系统包括：获取模块、奇异值分解模块、低秩模型模块、优化模块和变换模块；

所述获取模块用于获取欠采样矩阵X₀；

所述奇异值分解模块用于对欠采样矩阵X₀进行奇异值分解，得到X₀＝U₀ΣV₀ ^T，其中，U₀为左奇异向量，V₀为右奇异向量，Σ为奇异值；

取所述欠采样矩阵X₀的前r个奇异值和与所述r个奇异值对应的左奇异向量和右奇异向量，且令Σ＝Σ_LΣ_R，U＝U₀Σ_L，V^T＝Σ_RV₀ ^T则X₀＝UV^T；

所述低秩模型模块用于将所述U和所述V的每一列都转化为汉克尔矩阵，并构建磁共振波谱低秩模型：

s.t.P_Ω(X)＝X₀

其中，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)＝X₀

采用交替方向乘子法求解所述最小秩优化问题，得到U^k+1和V^k+1,并补全所述矩阵X：

X＝X^k+1＝P_Ω(X₀)+P_Ωc(U^k+1(V^k+1)^T)

其中X^k+1为所述矩阵X经过k+1次迭代后的矩阵；U^k+1为所述矩阵U经过k+1次迭代后的矩阵；V^k+1为所述矩阵V经过k+1次迭代后的矩阵；P_Ωc是P_Ω的补充，其值与P_Ω相反，即P_Ω为1的位置，P_Ωc在相对应的位置为0；P_Ω为0的位置，P_Ωc在相对应的位置为1；

所述优化模块具体用于：

将所述最小秩优化问题转化为增广拉格朗日函数：

s.t.P_Ω(X)＝X₀

其中，Di,Mi为拉格朗日乘子，μ为正则化参数，μ>0,<·，·>表示矩阵在希尔伯特空间的内积，<A,B>＝trace(AB),trace(·)表示矩阵的迹；

对所述增广拉格朗日函数采用交替方向乘子法交替求解以下子问题得到：

s.t.P_Ω(X)＝X₀

通过求解以下最小二乘问题得到U：

求解以下最小二乘问题得到V：

通过以下公式迭代更新X：

X^k+1＝X₀+P_Ωc(U^k+1(V^k+1)^T)

通过以下公式迭代更新B_i、C_i：

其中，其中k为迭代次数；R*为操作算子R的伴随算子；Qi*为Qi的伴随算子；

为奇异值缩减算子，

为所述B_i经过k+1次迭代后的结果，

为所述C_i经过k+1次迭代后的结果,μ^k为所述μ经过k次迭代后的结果；

所述奇异值缩减算子

为硬阈值算子，矩阵的秩为奇异值的0范数，B_i的优化子问题表示为：

所述C_i的优化子问题表示为：

其中，σ_Bi表示矩阵Bi的奇异值；σ_Ci表示矩阵Ci的奇异值；||·||₀表示矩阵的0范数；

所述变换模块用于当所述X^k+1为二维磁共振信号时，对所述X^k+1做二维傅里叶变换得到完整的二维磁共振波谱。

Images