CN106646303A - 一种欠采样磁共振波谱的快速重建方法 - Google Patents

Abstract

一种欠采样磁共振波谱的快速重建方法，涉及磁共振波谱。首先对磁共振波谱的时间域信号进行欠采样，将采样数据填充在一个汉克尔矩阵当；然后利用矩阵的低秩特性重建出完整的时间域信号，再经傅里叶变换获得高质量波谱；在此过程中，利用矩阵因子分解方法，通过矩阵弗罗贝尼乌斯范数项的最小化来逼近矩阵核范数项的最小化，避免在求解核范数项时使用奇异值分解，从而达到加速重建的效果。效果优良，易于操作，可以应用在一维、二维以及更高维的欠采样的磁共振波谱重建中。

Description

一种欠采样磁共振波谱的快速重建方法

技术领域

本发明涉及磁共振波谱，尤其是涉及基于时间域信号低秩汉克尔矩阵的一种欠采样磁共振波谱的快速重建方法。

背景技术

磁共振波谱被广泛地应用在临床医学和生物技术领域，其中高维磁共振波谱因能提供丰富信息而备受青睐。但对于高维核磁共振波谱来说，采样时间长是亟待解决的问题之一。用欠采样(相对于满足奈奎斯特准则下的全采样)来降低间接维的采样点数是一种缩短核磁共振采样时间的方法。但欠采样因为不满足奈奎斯特采样条件，容易造成谱峰重叠进而形成伪峰。为了获得高质量的谱，可以通过谱的自稀疏性从欠采样数据来重建这些谱(Xiaobo Qu,Xue Cao,Di Guo,Zhong Chen,"Compressed sensing for sparse magneticresonance spectroscopy,"International Society for Magnetic Resonance inMedicine 18th Scientific Meeting.Stockholm,Sweden,pp.3371,2010.)(Xiaobo Qu,DiGuo,Xue Cao,Shuhui Cai,Zhong Chen,"Reconstruction of self-sparse 2D NMRspectra from undersampled data in indirect dimension,"Sensors,vol.11,no.9,pp.8888-8909,2011.)。然而，谱峰宽度的增加将降低谱的自稀疏性，容易导致重建的宽谱峰发生畸变甚至完全丢失。屈小波等(Xiaobo Qu,Maxim Mayzel,Jian-Feng Cai,ZhongChen,Vladislav Orekhov."Accelerated NMR spectroscopy with low-rankreconstruction,"Angewandte Chemie International Edition,vol.54,no.3,pp.852-854,2015.)提出了一种基于低秩汉克尔矩阵的谱重建方法，该方法能同时重建不同宽度的谱峰，能够使用更少的欠采样数据点重建出高质量的谱。但该方法在迭代求解中最小核范数项时需要进行复杂度较高的奇异值分解，导致谱重建时间较长。

发明内容

本发明的目的在于提供效果优良、易于操作、重建时间短的欠采样磁共振波谱重建方法。

本发明包括以下步骤：

1)欠采样磁共振波谱的时间域信号，具体方法可为：

对磁共振信号进行欠采样，给定一组与预期的全采情况下信号的维度相同的且取值为0或1的数字，用于表示全采样磁共振波谱的特定时间点的信号是否采集，其中0表示不采样，1表示采样，根据0和1构成欠采样算子U；

2)构建时间域信号的汉克尔矩阵，具体方法可为：

假设x是待重建的磁共振波谱时间域信号，其维度与全采情况下数据的维度相同，然后通过线性算子R构建汉克尔矩阵，一维的磁共振波谱信号通过线性算子R构建成汉克尔矩阵，二维的磁共振波谱信号通过线性算子R构建成分块汉克尔矩阵，更高维的磁共振信号通过线性算子R构建成多维汉克尔矩阵；

3)设基于矩阵因子分解的低秩汉克尔矩阵的重建模型：

式(1)中，x为待重建磁共振波谱时间域信号，U为欠采样算子，y为欠采的磁共振波谱时间域信号，R是构建汉克尔矩阵的线性算子；||·||₂为向量的二范数，λ是一个正则化参数，用于权衡和两项的重要性，符号“H”表示复共轭转置，||·||_F表示矩阵的弗罗贝尼乌斯范数，P,Q是满足约束条件的任意矩阵；

4)采用基于矩阵因子分解的低秩汉克尔矩阵的重建算法：

式(1)的增广拉格朗日形式可以写成式(2)：

其中D是拉格朗日乘子，其维度与Rx相同，β为大于零的参数，<·,·>表示内积；

然后通过式(3)进行迭代求解：

其中，下标k表示第k次的解，符号“-1”表示求矩阵的逆，τ是一个大于0的参数。当k＝1时，P_k和Q_k初始为随机矩阵，D_k是一个全为1的矩阵；当达到迭代停止准则时，迭代停止；

在步骤4)中，所述迭代停止的准则设定为达到最大迭代次数10³或x在相邻两次迭代中的误差小于设定的阈值η，η设为10^-6。

5)对补全的时间信号进行傅立叶变换即可得到谱图。

本发明利用弗罗贝尼乌斯范数项的最小化来逼近核范数项的最小化的性质(Srebro,Nathan."Learning with matrix factorizations,"Massachusetts Instituteof Technology,2004.)，利用矩阵的因子分解来避免迭代过程中的奇异值分解，通过对磁共振波谱的时间信号进行汉克尔矩阵的低秩重建，并设计对应的迭代重建算法，实现包含多维磁共振波谱在内的欠采样信号的快速重建，进而得到高质量的磁共振谱。本发明首先对磁共振波谱的时间域信号进行欠采样，将采样数据填充在一个汉克尔矩阵当；然后利用矩阵的低秩特性重建出完整的时间域信号，再经傅里叶变换获得高质量波谱；在此过程中，利用矩阵因子分解方法，通过矩阵弗罗贝尼乌斯范数项的最小化来逼近矩阵核范数项的最小化，避免在求解核范数项时使用奇异值分解，从而达到加速重建的效果。效果优良，易于操作，可以应用在一维、二维以及更高维的欠采样的磁共振波谱重建中。

附图说明

图1是实施例中的采样模板。

图2是全采样的谱。

图3是本发明方法重建的谱。

具体实施方式

本实施例将重建二维磁共振波谱，直接维与间接维大小分别为N₁＝256和N₂＝116，包含87个谱峰。根据欠采样模板对磁共振波谱时间域信号进行欠采样，采集30％的数据，则本实施例中总的磁共振波谱数据点数为29696点，欠采样率为30％时得到的总采样数据点数为8909点。正则化参数λ＝10⁶，β＝1。参见图1，本实施例的具体步骤如下：

1)欠采样磁共振波谱的时间域信号：对二维磁共振信号进行欠采样，给定一组与预期的全采情况下信号的维度相同的且取值为0或1的数字(如图2所示，黑色表示0，白色表示1)，用于表示对全采样磁共振波谱的特定时间点的信号是否采集，其中0表示不采样，1表示采样。根据0和1构成欠采样算子U。

2)构建时间域信号的分块汉克尔矩阵：假设x是待重建的磁共振波谱时间域信号，其维度为256×116，然后通过线性算子R构建分块汉克尔矩阵，得到维度为7424×7611的分块汉克尔矩阵。

3)设基于矩阵因子分解的低秩汉克尔矩阵的重建模型：

式(1)中，x为待重建磁共振波谱时间域信号，U为欠采样算子，维度为256×116，y为欠采的磁共振波谱时间域信号，此实施例中为欠采样信号中丢失数据点进行填零后的欠采样数据，其维度为256×116，R是构建汉克尔矩阵的线性算子；||·||₂为向量的二范数，正则化参数λ＝10⁶，符号“H”表示复共轭转置，||·||_F表示矩阵的弗罗贝尼乌斯范数，P,Q是满足约束条件的矩阵，此实施例中P的维度为7424×400，Q的维度为7611×400。

4)基于矩阵因子分解的低秩汉克尔矩阵的重建算法：

式(1)的增广拉格朗日形式可以写成式(2)：

其中D是拉格朗日乘子，其维度是7424×7611；β为大于零的参数，取值为1；<·,·>表示内积。

然后通过式(3)进行迭代求解如下：

其中，下标k表示第k次的解，符号“-1”表示求矩阵的逆，τ的取值为1。初值化算法(也就是k＝1时)中，P_k和Q_k初始为随机矩阵，D_k是一个全为1的矩阵。当达到迭代停止准则时，迭代停止。迭代停止准则设定为达到最大迭代次数10³或x在相邻两次迭代中的误差小于设定的阈值η，η设为10^-6。

5)数据后处理：对重建的时间信号进行傅立叶变换即可得到谱图(如图3所示)。作为参考，我们将原始的全采样时间信号做傅立叶变换得到谱图。可以看出，利用采集到的部分数据和本发明的数据补全方法，可以重建得到高质量的谱。

本发明的谱重建耗时是采用奇异值分解的谱重建(Xiaobo Qu,Maxim Mayzel,Jian-Feng Cai,Zhong Chen,Vladislav Orekhov."Accelerated NMR spectroscopy withlow-rank reconstruction,"Angewandte Chemie International Edition,vol.54,no.3,pp.852-854,2015.)耗时的7.6％，奇异值分解的谱重建时间与按照本发明方法的谱重建时间如表1所示。

表1

磁共振波谱重建方法	重建耗时(单位：s)
		基于奇异值分解的磁共振波谱重建	13758.4
本发明	1052.2

由表1可见,本发明可以达到快速重建波谱的目的。

Claims (4)

1.一种欠采样磁共振波谱的快速重建方法，其特征在于包括以下步骤：

1)欠采样磁共振波谱的时间域信号；

2)构建时间域信号的汉克尔矩阵；

3)设基于矩阵因子分解的低秩汉克尔矩阵的重建模型：

$\begin{array}{cccc}\underset{x,P,Q}{min}& \frac{1}{2}\left(\right|\left|P\right|{|}_F^2+\left|\right|Q\left|{|}_F^2\right)+\frac{\mathrm{\&lambda;}}{2}\left|\right|y-Ux|{|}_2^2& s.t.& Rx={\mathrm{PQ}}^H\end{array}---\left(1\right)$

式(1)中，x为待重建磁共振波谱时间域信号，U为欠采样算子，y为欠采的磁共振波谱时间域信号，R是构建汉克尔矩阵的线性算子；||·||₂为向量的二范数，λ是一个正则化参数，用于权衡和两项的重要性，符号“H”表示复共轭转置，||·||_F表示矩阵的弗罗贝尼乌斯范数，P,Q是满足约束条件的任意矩阵；

4)采用基于矩阵因子分解的低秩汉克尔矩阵的重建算法：

式(1)的增广拉格朗日形式写成式(2)：

$\underset{x,P,Q,D}{min}\frac{1}{2}\left|\right|P|{|}_F^2+\frac{1}{2}|\left|Q\right|{|}_F^2+\frac{\mathrm{\&lambda;}}{2}\left|\right|y-Ux|{|}_2^2+<D,Rx-{\mathrm{PQ}}^H>+\frac{\mathrm{\&beta;}}{2}||Rx-{\mathrm{PQ}}^H|{|}_F^2---\left(2\right)$

其中D是拉格朗日乘子，其维度与Rx相同，β为大于零的参数，<·,·>表示内积；

然后通过式(3)进行迭代求解：

$\left\{\begin{array}{c}{x}_{k+1}={({\mathrm{\&lambda;U}}^{T}U+{\mathrm{\&beta;R}}^{H}R)}^{-1}\&lsqb;{\mathrm{\&lambda;U}}^{H}y+{\mathrm{\&beta;R}}^H(P_k{Q}_k^H-D_k)\&rsqb;\\ P_{k+1}=\mathrm{\&beta;}({\mathrm{Rx}}_{k+1}+D_k)Q_k{({\mathrm{\&beta;Q}}_k^H{Q}_k+I)}^{-1}\\ Q_{k+1}=\mathrm{\&beta;}{({\mathrm{Rx}}_{k+1}+D_k)}^H{P}_{k+1}{({\mathrm{\&beta;P}}_k^H{P}_k+I)}^{-1}\\ D_{k+1}=D_k+{\mathrm{\&tau;}}_k(Rx-{\mathrm{PQ}}^H)\end{array}\right.---\left(3\right)$

其中，下标k表示第k次的解，符号“-1”表示求矩阵的逆，τ是一个大于0的参数；当k＝1时，P_k和Q_k初始为随机矩阵，D_k是一个全为1的矩阵；当达到迭代停止准则时，迭代停止；

5)对补全的时间信号进行傅立叶变换即得到谱图。

2.如权利要求1所述一种欠采样磁共振波谱的快速重建方法，其特征在于在步骤1)中，所述对欠采样磁共振波谱的时间域信号的具体方法为：磁共振信号进行欠采样，给定一组与预期的全采情况下信号的维度相同的且取值为0或1的数字，用于表示全采样磁共振波谱的特定时间点的信号是否采集，其中0表示不采样，1表示采样，根据0和1构成欠采样算子U。

3.如权利要求1所述一种欠采样磁共振波谱的快速重建方法，其特征在于在步骤2)中，所述构建时间域信号的汉克尔矩阵的具体方法为：假设x是待重建的磁共振波谱时间域信号，其维度与全采情况下数据的维度相同，然后通过线性算子R构建汉克尔矩阵，一维的磁共振波谱信号通过线性算子R构建成汉克尔矩阵，二维的磁共振波谱信号通过线性算子R构建成分块汉克尔矩阵，更高维的磁共振信号通过线性算子R构建成多维汉克尔矩阵。

4.如权利要求1所述一种欠采样磁共振波谱的快速重建方法，其特征在于在步骤4)中，所述迭代停止的准则设定为达到最大迭代次数10³或x在相邻两次迭代中的误差小于设定的阈值η，η设为10^-6。