CN104793159A - 一种高维核磁共振时域信号补全方法 - Google Patents

Abstract

一种高维核磁共振时域信号补全方法，涉及核磁共振高维谱信号处理。先根据采集得到的数据和给定的核磁共振谱频谱宽度和分辨率，确定需要补全的时域信号的位置并设计模板；再利用提出的高维核磁共振时域信号补全方法来构建重建模型，然后通过最优化算法求解出完整的高维核磁共振时域信号；最后对补全后的时域信号做傅立叶变换得到核磁共振谱。实现了在高维核磁共振实验中对有丢失的核磁共振时域信号进行信号补全，进而得到完整的核磁共振时域信号。由于利用了高维核磁共振信号自身的特征，可以补全任意高维核磁共振时域信号。可以达到降低采样时间、提高信噪比、达到给定核磁共振谱频谱宽度和分辨率的目的。

Description

一种高维核磁共振时域信号补全方法

技术领域

本发明涉及核磁共振高维谱信号处理，特别是涉及一种高维核磁共振时域信号补全方法。

背景技术

核磁共振波谱(NMR)，简称为核磁共振谱，已经成为物理学、化学以及生命科学等诸多学科研究物质组成的成分、组织形态及其变化、脑功能、分子结构和动力学强有力的手段。特别在化学结构分析领域，核磁共振谱有着重要的应用。

随着检测对象和设备的不断发展，高维NMR(≥2维)在应用领域发挥着越来越重要的作用。例如，2维、3维和4维NMR广泛用于解析复杂蛋白结构(D.Sakakibara,A.Sasaki,T.Ikeya,J.Hamatsu,T.Hanashima,M.Mishima,et al.,"Protein structure determination in living cells byin-cell NMR spectroscopy,"Nature,vol.458,pp.102-105,2009)。但数据采集时间随着维数的增加而指数级增长。2维谱到4维谱的采集时间从几分钟不断上升到几十天(M.Mobli and J.C.Hoch,"Nonuniform sampling and non-Fourier signal processing methods in multidimensionalNMR,"Progress in Nuclear Magnetic Resonance Spectroscopy,vol.83,pp.21-41,2014)。冗长的高维蛋白质NMR实验时间，不但使得实验必须耗费大量的谱仪机时，而且提高了不稳定蛋白质样品的实验难度，从而限制了高维NMR技术的在蛋白质研究中的应用。为了缩短高维NMR实验时间，一个常见的做法是采用非均匀采样的方式来减少采样的点数，但同时也造成了信号的丢失，需要信号补全。

在NMR实验中，采样得到的时域信号不可避免会受到噪声的影响。由于核磁共振信号末尾处的信号强度小造成这些数据点的信噪比低，截掉这些数据点有利于提高核磁共振谱的信噪比，但这也会造成信号丢失，需要信号补全。

频谱宽度和分辨率是评价核磁共振谱质量的两个重要指标。当采样所得的时间信号未达到预期的频谱宽度和分辨率时，需要对时域信号进行信号插值。信号插值也是一种信号补全。

因此信号补全技术在NMR实验中应用很广泛，设计有效的信号补全技术意义重大。

发明内容

本发明的目的在于针对高维核磁共振波谱实验存在的上述问题，提供一种高维核磁共振时域信号补全方法。

本发明包括以下步骤：

1)设计模板确定信号丢失位置：根据实验中的采样方式和给定的核磁共振谱频谱宽度和分辨率，确定需要补全的信号位置并设计成模板，模板大小与预期完整的信号大小相同；模板由数字0和1构成，数字0表示该位置对应的信号需要补全，数字1表示该位置对应的信号未丢失；

2)构建信号重建模型：在m维核磁共振实验中，预期完整的信号X是一个m维矩阵，其中第一维用j₁表示，这一维的信号长度是J₁，第二维用j₂表示，这一维的信号长度是J₂，以此类推；第一维是直接维，其余是间接维，当给定除j_k以外的其它维度时，可以得到一个一维数据，称为一维阵列，用表示；根据核磁共振信号具有指数函数形式的特点，信号中沿任意维度的一个一维阵列排列成一个汉克尔矩阵，该矩阵的秩比行数小很多，这种矩阵称为低秩矩阵；构造以下的模型来求解完整的核磁共振信号X：

$\begin{array}{cc}\underset{x}{\min }\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{D}_k& s.t.{\left|\right|y-\mathrm{AX}\left|\right|}_2^2\≤\mathrm{\ϵ}---\left(1\right)\end{array}$

其中，

$D_k=\underset{j_m=1}{\overset{J_m}{\mathrm{\Σ}}}\·\·\·\underset{j_{k+1}=1}{\overset{J_{k+1}}{\mathrm{\Σ}}}\underset{j_{k-1}=1}{\overset{J_{k-1}}{\mathrm{\Σ}}}\·\·\·\underset{j_1=1}{\overset{J_1}{\mathrm{\Σ}}}{\left|\right|H_{j_1{j}_2\·\·\·j_{k-1}{j}_{k+1}\·\·\·j_m}\left|\right|}_{*}---\left(2\right)$

$H_{j_1{j}_2\·\·\·j_{k-1}{j}_{k+1}\·\·\·j_m}={\mathrm{RX}}_{j_1{j}_2\·\·\·j_{k-1}:j_{k+1}\·\·\·j_m}---\left(3\right)$

在公式(1)中，算子A根据步骤1)中设计的模板确定X中信号丢失的位置，并把未丢失的信号排列成一维向量；y表示采集到并且没有丢失的数据；表示(y-AX)中的元素的平方和；表示将一维阵列排列成一个汉克尔矩阵 表示矩阵的核范数，即对矩阵的奇异值求和；参数ε与实际存在的噪声大小有关，ε＞0；

求解公式(1)中的模型时，可以把约束优化问题转化为以下的非约束优化问题，

$\begin{array}{c}\underset{x}{\min }\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{D}_k+\frac{\mathrm{\λ}}{2}{\left|\right|y-\mathrm{AX}\left|\right|}_2^2---\left(4\right)\end{array}$

其中λ是正则化参数(λ>0)，用于权衡和两项的重要性；

3)求解最优化问题：求解最优化问题即公式(4)，即可得到完整的m维核磁共振信号。

在步骤3)中，所述求解最优化问题即公式(4)是一个核范数最优化问题，可以采用奇异值阈值化法(J.F.Cai,E.J.Candes,and Z.W.Shen,"A singular value thresholding algorithm formatrix completion,"SIAM Journal on Optimization,vol.20,pp.1956-1982,2010.)(张贤达，矩阵分析与应用[M].第二版，北京：清华大学出版社，2013.)和交替方向乘子法(X.Qu,M.Mayzel,J.-F.Cai,Z.Chen,and V.Y.Orekhov,Accelerated NMR spectroscopy with low-rankreconstruction,Angewandte Chemie International Edition,vol.54,pp.852-854,2015.)等数值计算方法求解，得到完整的m维核磁共振时间信号X。然后对X进行傅立叶变换得到核磁共振谱。

本发明的突出效果如下：由于本发明利用了高维核磁共振信号自身的特征，可以补全任意高维核磁共振时域信号。通过采用本发明提出的一种高维核磁共振时域信号补全方法，可以达到降低采样时间、提高信噪比、达到给定核磁共振谱频谱宽度和分辨率的目的。

附图说明

图1为信号丢失的示意图。

图2为信号丢失位置模板的示意图。

图3为实施例中设计的信号丢失位置模板。

图4为信号补全后所得的核磁共振谱。

图5为全采样所得的核磁共振谱。

具体实施方式

下面通过具体实施例对本发明作进一步的描述，并给出补全的结果。

本实施例是一个二维核磁共振实验。

第一步：设计模板确定信号丢失位置

本实施例中，预期完整的核磁共振信号X是一个二维矩阵，直接维有256点，间接维有128点。在实验中，对直接维进行均匀采样，采样点数是256；对间接维进行非均匀采样，采样点数是32，从而得到二维的时域信号。为了提高信噪比，截掉直接维尾部64个点。为了达到给定的频率分辨率，需要对间接维的尾部做信号插值。图1是信号丢失的示意图，图2是信号丢失位置模板的示意图，“1”代表信号未丢失，“0”代表信号丢失，需要补全。图3是实施例中设计的信号丢失位置模板。模板中白色的点代表“1”，表示该位置对应的信号未丢失，黑色的点代表“0”，表示该位置对应的信号丢失，需要补全，模板大小为256×128。第二步：构建信号重建模型

本实施例中，预期完整的核磁共振信号用X表示，j₁表示X的直接维，直接维有256点，j₂表示间接维，间接维有128点。当给定j₁时，得到一个一维阵列当给定j₂时，得到一个一维阵列构造如下优化模型。

$\begin{array}{cc}\underset{x}{\min }\underset{j_1=1}{\overset{256}{\mathrm{\Σ}}}{\left|\right|{\mathrm{RX}}_{j_1:}\left|\right|}_{*}+\underset{j_2=1}{\overset{128}{\mathrm{\Σ}}}{\left|\right|{\mathrm{RX}}_{:j_2}\left|\right|}_{*}& s.t.{\left|\right|y-\mathrm{AX}\left|\right|}_2^2\≤\mathrm{\ϵ}---\left(5\right)\end{array}$

在公式(5)中，X表示预期完整的核磁共振信号，算子A根据步骤1)中设计的模板确定X中信号丢失的位置，并把未丢失的信号排列成一维向量；y表示采集到并且没有丢失的数据；表示将(y-AX)中的每个元素求平方，再求和；表示将一维阵列排列成一个汉克尔矩阵表示将一维阵列排列成一个汉克尔矩阵表示求矩阵的核范数，即对矩阵的奇异值求和；表示求矩阵的核范数；参数ε与实际存在的噪声大小有关，本实施例中取ε＝0.1。

在求解最优化问题(5)时，可以把约束优化问题转化为如下的非约束优化问题

$\underset{x}{\max }\underset{j_1=1}{\overset{256}{\mathrm{\Σ}}}{\left|\right|{\mathrm{RX}}_{j_1:}\left|\right|}_{*}+\underset{j_2=1}{\overset{128}{\mathrm{\Σ}}}{\left|\right|{\mathrm{RX}}_{{:j}_2}\left|\right|}_{*}+\frac{\mathrm{\λ}}{2}{\left|\right|y-\mathrm{AX}\left|\right|}_2^2---\left(6\right)$

其中λ是正则化参数，本实施例中取λ＝10⁴。

第三步：求解最优化问题

求解最优化问题(6)即可得到完整的m维核磁共振信号。公式(6)是一个核范数最优化问题，本实施例采用交替方向乘子法(X.Qu,M.Mayzel,J.-F.Cai,Z.Chen,and V.Y.Orekhov,Accelerated NMR spectroscopy with low-rank reconstruction,Angewandte Chemie InternationalEdition,vol.54,pp.852-854,2015.)求解，得到完整的m维核磁共振时间信号X。然后对X进行傅立叶变换得到核磁共振谱。

在核磁共振实验中，由于截掉时域信号末尾的数据可以提高信噪比，欠采样信号能够缩短采样时间等原因，实际采集到的核磁共振时域信号经常存在信号丢失，因此需要信号补全。为了达到给定的核磁共振谱频谱宽度和分辨率，需要对时域信号插值，信号插值也是一种信号补全方法。本发明的目的在于提供了一种高维核磁共振时域信号补全方法。本发明先根据采集得到的数据和给定的核磁共振谱频谱宽度和分辨率，确定需要补全的时域信号的位置并设计模板；再利用提出的高维核磁共振时域信号补全方法来构建重建模型，然后通过最优化算法求解出完整的高维核磁共振时域信号。最后对补全后的时域信号做傅立叶变换得到核磁共振谱。本发明实现了在高维核磁共振实验中对有丢失的核磁共振时域信号进行信号补全，进而得到完整的核磁共振时域信号。

Claims (3)

1.一种高维核磁共振时域信号补全方法，其特征在于包括以下步骤：

1)设计模板确定信号丢失位置：根据实验中的采样方式和给定的核磁共振谱频谱宽度和分辨率，确定需要补全的信号位置并设计成模板，模板大小与预期完整的信号大小相同；模板由数字0和1构成，数字0表示该位置对应的信号需要补全，数字1表示该位置对应的信号未丢失；

2)构建信号重建模型：在m维核磁共振实验中，预期完整的信号X是一个m维矩阵，其中第一维用j₁表示，这一维的信号长度是J₁，第二维用j₂表示，这一维的信号长度是J₂，以此类推；第一维是直接维，其余是间接维，当给定除j_k以外的其它维度时，可以得到一个一维数据，称为一维阵列，用表示；根据核磁共振信号具有指数函数形式的特点，信号中沿任意维度的一个一维阵列排列成一个汉克尔矩阵，该矩阵的秩比行数小很多，这种矩阵称为低秩矩阵；构造以下的模型来求解完整的核磁共振信号X：

$\begin{array}{cc}\underset{x}{\min }\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{D}_k& s.t.{\left|\right|y-\mathrm{AX}\left|\right|}_2^2\≤\mathrm{\ϵ}---\left(1\right)\end{array}$

其中，

$D_k=\underset{j_m=1}{\overset{J_m}{\mathrm{\Σ}}}...\underset{j_{k+1}=1}{\overset{J_{k+1}}{\mathrm{\Σ}}}\underset{j_{k-1}=1}{\overset{J_{k-1}}{\mathrm{\Σ}}}...\underset{J_1=1}{\overset{J_1}{\mathrm{\Σ}}}{\left|\right|H_{j_1{j}_2...j_{k-1}{j}_{k+1}...j_m}\left|\right|}_{*}---\left(2\right)$

$H_{j_1{j}_2...j_{k-1}{j}_{k+1}...j_m}={\mathrm{RX}}_{j_1{j}_2...j_{k-1}:j_{k+1}...j_m}---\left(3\right)$

在公式(1)中，算子A根据步骤1)中设计的模板确定X中信号丢失的位置，并把未丢失的信号排列成一维向量；y表示采集到并且没有丢失的数据；表示(y-AX)中的元素的平方和；表示将一维阵列排列成一个汉克尔矩阵 表示矩阵的核范数，即对矩阵的奇异值求和；参数ε与实际存在的噪声大小有关，ε＞0；

求解公式(1)中的模型时，可以把约束优化问题转化为以下的非约束优化问题，

$\underset{X}{\min }\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{D}_k+\frac{\mathrm{\λ}}{2}{\left|\right|y-\mathrm{AX}\left|\right|}_2^2---\left(4\right)$

其中λ是正则化参数(λ>0)，用于权衡和两项的重要性；

3)求解最优化问题：求解最优化问题即公式(4)，即可得到完整的m维核磁共振信号。

2.如权利要求1所述一种高维核磁共振时域信号补全方法，其特征在于在步骤3)中，所述求解最优化问题即公式(4)是一个核范数最优化问题，采用数值计算方法求解，得到完整的m维核磁共振时间信号X；然后对X进行傅立叶变换得到核磁共振谱。

3.如权利要求2所述一种高维核磁共振时域信号补全方法，其特征在于所述数值计算方法包括但不限于奇异值阈值化法或交替方向乘子法。