CN108920423A - 一种高保真谱重建方法 - Google Patents

Abstract

一种高保真谱重建方法，涉及谱的重建方法。将待恢复谱的时域信号构建为汉克尔矩阵；利用逼近函数来近似计算矩阵的秩；建立一种逼近矩阵秩的低秩重建模型；提出谱的时域信号重建模型的求解算法；数据后处理：对求解得到的时域信号进行傅里叶变换可得到频谱。在磁共振波谱采集中，经常需要获得大量的数据。对大量的数据进行采集需要消耗大量时间，一种方式是通过采集部分信号来加速数据采集。从谱的时域信号的汉克尔矩阵的低秩特性出发来恢复出完整的信号，首先利用逼近函数来近似计算汉克尔矩阵的秩，然后建立谱信号的重建模型，最后通过迭代算法重建信号。重建的谱的精度高，易于操作，可以从少量数据中恢复出高保真的完整谱信号。

Description

一种高保真谱重建方法

技术领域

本发明涉及谱的重建方法，尤其是涉及基于低秩近似的一种高保真谱重建方法。

背景技术

在许多实际应用中，如磁共振波谱和雷达目标定位，我们感兴趣的目标信号可以建模成在频域(相对时域)上若干谱峰的叠加，而采集的数据是时域(相对频域)信号可以表示为一系列指数信号的叠加。在实际的采样过程中，由于受到硬件、物理条件的限制，不得不加快采样速度，因此实际采样得到的数据并不完整或达不到预期的分辨率，需要重建采集到的数据中丢失部分。特别是在高维应用领域，通常数据量非常大，全采样时间过于冗长，往往在测量时采用非均匀欠采样缩短采样时间，通过重建的方法得到完整的数据和预期的分辨率。

以磁共振波谱为例，它在化学分子结构分析领域有着重要的应用，但磁共振实验时间较长，从几分钟到几十天不等。这不但使实验耗费大量谱仪机时，而且增加了不稳定样品的实验难度，从而限制了高维磁共振技术在研究中的应用。其时域信号符合指数函数特征，因此其信号转为的汉克尔矩阵具有低秩特性。为了缩短磁共振实验时间，可以采用欠采样降低间接维采样点数。然而，欠采样容易造成谱峰重叠进而形成伪峰。为了获得高质量频谱，可以通过谱的自稀疏性对欠采样数据进行重建(Xiaobo Qu,Xue Cao,Di Guo,ZhongChen."Compressed sensing for sparse magnetic resonance spectroscopy,"International Society for Magnetic Resonance in Medicine 19th ScientificMeeting,Stockholm,Sweden,pp.3371,2010.)；也可以利用磁共振波谱时间域信号对应汉克尔矩阵的低秩性来实现重建(Xiaobo Qu,Maxim Mayzel,Jian-Feng Cai,Zhong Chen,Vladislav Orekhov."Accelerated NMR spectroscopy with low-rankreconstruction,"AngewandteChemie InternationalEdition,vol.54,no.3,pp.852-854,2015.)，但这些方法在采样率较低时效果不佳，谱峰容易失真。

发明内容

本发明的目的在于提供一种高保真谱重建方法。

本发明包括以下步骤：

1)将待恢复谱的时域信号构建为汉克尔矩阵；

在步骤1)中，所述将待恢复谱的时域信号构建为汉克尔矩阵的具体方法可为：将待恢复谱的时域信号记为x＝[x(1),x(2),…,x(N)]，信号长度为N，其中c_j和z_j均为复数，J为正整数，表示信号x中所包含单指数信号的个数，n为指数的次数；通过线性算子构建x为汉克尔矩阵：

上式算子有两个参数Q和P，分别决定汉克尔矩阵的行数和列数。

2)利用逼近函数来近似计算矩阵的秩；

在步骤2)中，所述利用逼近函数来近似计算矩阵的秩具体方法可为：利用非凸函数可以近似汉克尔矩阵的秩，其中表示汉克尔矩阵第g大的奇异值，φ被定义为：

3)建立一种逼近矩阵秩的低秩重建模型如下：

其中，为欠采样算子，y为采集到的信号，x为待恢复信号，λ是平衡和的正则化参数；

4)提出谱的时域信号重建模型的求解算法；

在步骤4)中，所述提出谱的时域信号重建模型的求解算法的具体方法可为：为求解式(3)中的重建模型，引入中间变量Z，将模型松弛化如下：

其中，β表示正则化参数，与λ共同权衡 和三项的重要性；

当β趋于无穷大时，式(4)的解将趋近式(3)的解，可以利用连续交替方向最小化方法，求解最优化问题式(4)，根据以下式(5)～(7)迭代更新变量：

其中，下标k表示第k次的解，符号“-1”表示求矩阵的逆，上标H为矩阵的共轭转置，对汉克尔矩阵进行奇异值分解可以得到而Z_k+1为引入的中间变量，对比式(3)及式(4)，将函数Θ(Σ_k+1；2a/β,a)定义为：

Θ(Σ_k+1；β,a)＝min{Σ_k+1,max{(Σ_k+1-2a/β)/(1-2a²/β),0}} (7)

其中，max{}表示取元素的最大值，min{}表示取元素的最小值；

函数Θ的作用是将奇异值矩阵Σ_k+1中的奇异值依次进行处理，具体计算过程为：首先将第k+1个奇异值矩阵Σ_k+1中保存的第s个奇异值Σ_s,k+1代入(Σ_s,k+1-2a/β)/(1-2a²/β)，保留集合{(Σ_s,k+1-2a/β)/(1-2a²/β)，0}中两个元素的较大者(如果相等就保留0)作为max{(Σ_s,k+1-2a/β)/(1-2a²/β),0}的结果；然后将max{(Σ_s,k+1-2a/β)/(1-2a²/β),0}与Σ_s,k+1比较，保留集合{Σ_s,k+1,max{(Σ_s,k+1-2a/β)/(1-2a²/β),0}}中两个元素的较小者(如果相等就保留Σ_s,k+1)；最后按照前述2个步骤修改奇异值矩阵Σ_k+1中所有的奇异值，作为Θ(Σ_k+1；β,a)的结果；

式(4)中参数β和λ是正数，当达到迭代停止准则时，迭代停止；迭代停止准则设定为达到最大迭代次数或x在相邻两次迭代中的误差小于设置的阈值η(取值大于0)；当迭代停止时，可根据式(5)得到完整的谱的时域信号。

5)数据后处理：对求解得到的时域信号进行傅里叶变换可得到频谱。

本发明在磁共振波谱采集中，经常需要获得大量的数据。对大量的数据进行采集需要消耗大量时间，一种方式是通过采集部分信号来加速数据采集。本发明从谱的时域信号的汉克尔矩阵的低秩特性出发来恢复出完整的信号，首先利用逼近函数来近似计算汉克尔矩阵的秩，然后建立谱信号的重建模型，最后通过迭代算法重建信号。本发明重建的谱的精度高，易于操作，可以从少量数据中恢复出高保真的完整谱信号。

附图说明

图1为本发明重建后所得的频谱。

图2为全采样的频谱。

具体实施方式

下面通过具体实施例对本发明作进一步的说明，并给出重建结果。本实施例是一个重建二维磁共振波谱的模拟实验，直接维与间接维的大小分别是M＝116和N＝256。根据欠采样模板对二维磁共振波谱时间域信号进行欠采样，采样25％数据，则本实施例中的磁共振波谱数据点为29696点，采样率为25％时得到的总采样数据点数为7424点。具体步骤如下：

1)对待恢复谱的时域信号构建汉克尔矩阵：选定二维磁共振波谱直接维的某一值，沿间接维抽取一条谱的时域信号，由此可以得到116条一维信号。将任意一条信号记为x＝[x(1),x(2),…,x(256)]，信号长度为256，其中c_j和z_j均为复数，J为正整数，表示信号x中所包含单指数信号的个数，n为指数的次数。通过线性算子构建x为汉克尔矩阵其中汉克尔矩阵的行数和列数分别是Q＝128和P＝129。

2)建立一种逼近矩阵秩的低秩重建模型为：

其中，y为采集到的信号，有38个点，x为待恢复信号，为欠采样算子，作用是使需要恢复的完整的信号x变换为欠采样信号y；λ是平衡和的正则化参数，表示汉克尔矩阵第g大的奇异值，φ被定义为：

3)提出谱的时域信号重建模型的求解算法：为求解式(1)中的重建模型，引入中间变量Z，将模型松弛化如下：

其中，β表示正则化参数，与λ共同权衡 和三项的重要性。

当β趋于无穷大时，式(3)的解将趋近式(1)的解。可以利用连续交替方向最小化方法求解最优化问题式(3)，根据以下式(4)～(6)迭代更新变量为：

其中，下标k表示第k次的解，符号“-1”表示求矩阵的逆，上标H为矩阵的共轭转置，对汉克尔矩阵进行奇异值分解可以得到而Z_k+1为引入的中间变量，对比式(1)及式(3)，将因此函数Θ(Σ_k+1；2a/β,a)定义为：

Θ(Σ_k+1；β,a)＝min{Σ_k+1,max{(Σ_k+1-2a/β)/(1-2a²/β),0}} (6)

其中，max{}表示取元素的最大值，min{}表示取元素的最小值。

函数Θ的作用是将奇异值矩阵Σ_k+1中的奇异值依次进行处理，具体计算过程为：

首先将第k+1个奇异值矩阵Σ_k+1中保存的第s个奇异值Σ_s,k+1代入(Σ_s,k+1-2a/β)/(1-2a²/β)，保留集合{(Σ_s,k+1-2a/β)/(1-2a²/β)，0}中两个元素的较大者(如果相等就保留0)作为max{(Σ_s,k+1-2a/β)/(1-2a²/β),0}的结果；

然后将max{(Σ_s,k+1-2a/β)/(1-2a²/β),0}与Σ_s,k+1比较，保留集合{Σ_s,k+1,max{(Σ_s,k+1-2a/β)/(1-2a²/β),0}}中两个元素的较小者(如果相等就保留Σ_s,k+1)；

最后按照前述2个步骤修改奇异值矩阵Σ_k+1中所有的奇异值，作为Θ(Σ_k+1；β,a)的结果。

当a越大，式(3)中的非凸性越强。通过不断更新β，使式(3)为凸函数的a范围越大。当β趋于无穷大时，式(3)的解将趋于式(1)的解。令式(3)中λ＝10⁴，初始值β₁＝1，a₁＝0,内层迭代达到停止准则时，内层迭代停止。内层迭代停止准则设定为达到最大迭代次数500次或x在相邻两次内层迭代中的误差小于设定的阈值10^-4。内层迭代停止后，令β₂＝2×β₁，不断进行内层迭代,内层迭代停止后即更新β_l+1＝2×β_l和直到外层迭代达到停止准则时，求解结束。外层迭代停止准则设定为β达到最大值256或x在相邻两次外层迭代中的误差小于设定的阈值10^-4。最后得到的x为重建的时域信号。

4)数据后处理：对求解得到的信号沿间接维进行傅里叶变换可得到完整的磁共振波谱(如图1所示)。作为参考，将原始的全采样时间信号做二维傅里叶变换得到磁共振波谱(如图2所示)。可以看出，利用采集到的部分数据和本发明的高保真谱重建方法，可以重建得到高质量的磁共振波谱。

Claims (4)

1.一种高保真谱重建方法，其特征在于包括以下步骤：

1)将待恢复谱的时域信号构建为汉克尔矩阵；

2)利用逼近函数来近似计算矩阵的秩；

3)建立一种逼近矩阵秩的低秩重建模型如下：

其中，为欠采样算子，y为采集到的信号，x为待恢复信号，λ是平衡和的正则化参数；

4)提出谱的时域信号重建模型的求解算法；

5)数据后处理：对求解得到的时域信号进行傅里叶变换可得到频谱。

2.如权利要求1所述一种高保真谱重建方法，其特征在于在步骤1)中，所述将待恢复谱的时域信号构建为汉克尔矩阵的具体方法为：将待恢复谱的时域信号记为x＝[x(1),x(2),…,x(N)]，信号长度为N，其中c_j和z_j均为复数，J为正整数，表示信号x中所包含单指数信号的个数，n为指数的次数；通过线性算子构建x为汉克尔矩阵：

上式算子有两个参数Q和P，分别决定汉克尔矩阵的行数和列数。

3.如权利要求1所述一种高保真谱重建方法，其特征在于在步骤2)中，所述利用逼近函数来近似计算矩阵的秩具体方法为：利用非凸函数近似汉克尔矩阵的秩，其中，表示汉克尔矩阵第g大的奇异值，φ被定义为：

4.如权利要求1所述一种高保真谱重建方法，其特征在于在步骤4)中，所述提出谱的时域信号重建模型的求解算法的具体方法为：为求解式(3)中的重建模型，引入中间变量Z，将模型松弛化如下：

其中，β表示正则化参数，与λ共同权衡和三项的重要性；

当β趋于无穷大时，式(4)的解将趋近式(3)的解，利用连续交替方向最小化方法，求解最优化问题式(4)，根据以下式(5)～(7)迭代更新变量：

其中，下标k表示第k次的解，符号“-1”表示求矩阵的逆，上标H为矩阵的共轭转置，对汉克尔矩阵进行奇异值分解得到而Z_k+1为引入的中间变量，对比式(3)及式(4)，将函数Θ(Σ_k+1；2a/β,a)定义为：

Θ(Σ_k+1；β,a)＝min{Σ_k+1,max{(Σ_k+1-2a/β)/(1-2a²/β),0}} (7)

其中，max{}表示取元素的最大值，min{}表示取元素的最小值；

函数Θ的作用是将奇异值矩阵Σ_k+1中的奇异值依次进行处理，具体计算过程为：首先将第k+1个奇异值矩阵Σ_k+1中保存的第s个奇异值Σ_s,k+1代入(Σ_s,k+1-2a/β)/(1-2a²/β)，保留集合{(Σ_s,k+1-2a/β)/(1-2a²/β)，0}中两个元素的较大者作为max{(Σ_s,k+1-2a/β)/(1-2a²/β),0}的结果；然后将max{(Σ_s,k+1-2a/β)/(1-2a²/β),0}与Σ_s,k+1比较，保留集合{Σ_s,k+1,max{(Σ_s,k+1-2a/β)/(1-2a²/β),0}}中两个元素的较小者最后按照前述2个步骤修改奇异值矩阵Σ_k+1中所有的奇异值，作为Θ(Σ_k+1；β,a)的结果；

式(4)中参数β和λ是正数，当达到迭代停止准则时，迭代停止；迭代停止准则设定为达到最大迭代次数或x在相邻两次迭代中的误差小于设置的阈值η，η取值大于0；当迭代停止时，根据式(5)得到完整的谱的时域信号。