2014, Vol. 31文章信息
- 郑慧, 韩明月, 胡炳文, 杨光
- ZHENG Hui, HAN Ming-yue, HU Bing-wen, YANG Guang
- 不同采样模式的固体DQ-SQ实验的压缩感知重建比较
- Comparison of Different Sampling Schemes in Compressed Sensing Reconstruction for DQ-SQ experiments
- 波谱学杂志, 2014, 31(4): 535-547
- Chinese Journal of Magnetic Resonance, 2014, 31(4): 535-547
-
文章历史
收稿日期: 2014-02-21
收修改稿日期: 2014-10-27
对于复杂的化学和生物系统，二维核磁共振谱相比核磁共振一维谱而言，可以提供更多有价值的结构信息.为了保证二维谱有足够的分辨率，我们需要在间接维（t1维）采集更多的点，并且随分辨率的增加而增加[1].二维谱一般在直接维采集N2个t2点（点间距为Δt2），并在间接维方向采集N1个t1点（点间距为Δt1）；若先将时域信号在t2方向做傅里叶变换，再沿t1方向做傅里叶变换，便得到二维频谱，这也就是所谓的傅里叶变换二维谱(FT2D).通常，间接维演化方向的谱宽(SW1=1/Δt1)决定了t1方向的最小采集点数N1，根据关系式N1=SW1/Δv1，可以得到谱的分辨率(Δv1=1/N1Δt1).
减少采样时间最常用的方法是不均匀采样--即降低间接维采样点数N1，但是这样采集到的谱图无法用FT2D直接重建获得好的频率谱.在过去的30多年中提出了很多基于不均匀采集信号的重建算法，如协方差法(Covariance)[2, 3]、最大熵法(Maximum Entropy, MaxEnt)[4, 5]、不均匀傅里叶变换法等[7-9]，这些算法根据不同类型的谱的特点，提出不同的假设，给出适合的解决方案.本文中我们主要讨论压缩感知(Compressed Sensing, CS)重建技术[6, 10, 11]，并根据固体二维(2D)双量子-单量子(DQ-SQ)谱的先验知识，提出本文算法.
压缩感知可以在信号采样率达不到奈奎斯特(Nyquist)采样定理要求的情况下，仍然可能完美地重建出高分辨的谱，同时可以在一定程度上达到去噪的效果[12, 13].压缩感知方法的应用需要满足一些前提：首先，待重建的谱图应具有稀疏性，或在某个变换域是稀疏的，即在某个变换域中，仅有少量的非零值.常用的稀疏变换包括单位矩阵变换、傅里叶变换、离散余弦变换、小波变换、有限差分等[14].二维谱，特别是DQ-SQ谱本身就有很高的稀疏性，可以很好地满足压缩感知对稀疏性的要求.其次，压缩感知需要采用不相干采样获取信号的离散样本，即满足约束等距性(Restricted Isometry Property, RIP)[15]的采样，然后通过非线性重建算法完美的重建信号.不相干采样有很多种方式，例如变密度、径向、均匀螺旋、变密度螺旋、变密度扰动螺旋采样等[16].
基于压缩感知的方法近年来非常流行，在无线通信[17]、太赫兹成像[18]、磁共振成像[13, 19]、磁共振波谱成像[20]、单像素成像[21]等诸多领域中取得了广泛的应用[22].压缩感知的方法可以较大程度地提高信号采集效率和图像质量.
最近压缩感知重建技术被引入到液体核磁共振领域[6, 23]，并取得了良好的效果.屈小波等人[14]指出对于自稀疏的谱，直接使用谱图作为稀疏域比使用小波域的效果更好，其稀疏变换也就是上文提到的单位矩阵变换.由于大部分核磁共振谱在频域都是稀疏的，且数据的维数越大，稀疏性越好[24]，因此在核磁共振波谱中直接利用频域作为稀疏域是可行的.CS技术一般都采用随时间呈指数下降的不均匀采样(NUS)方式[4, 25].NUS采样方式一般在FID前面大量采样，在FID的后面少量采样，采样密度随时间的增加呈现指数型衰减，本文中我们称此类NUS为expNUS方式[4].
expNUS采样得到的时域数据不能用常规的FT2D转换得到频谱图.因此压缩感知算法被应用在二维谱的重建上.本文主要引入了不同于expNUS的新的采样方式，即伪随机采样方式(pseudo-random sampling, PRAS).我们在介绍了本文算法和伪随机采样的实现机制之后，比较了29Si DQ-SQ与13C DQ-SQ的二维谱在相同采样率ξ下不同连续采样区域大小对重建结果的影响，最后分别比较了伪随机采样、全随机采样以及e指数不相干采样这3种采样模式分别应用压缩感知的重建结果，发现对于二维的DQ-SQ核磁共振谱，用伪随机采样模式重建出来的谱图效果最好.
1 基本理论 1.1 压缩感知压缩感知用在二维谱上可以这样描述[13]，假设二维频域数据矩阵用m表示，ψ为反傅里叶变换，Θ为欠采样操作.重建过程是解决以下约束最优化问题：
| $ \hat m = {\rm{min}}{\left\| m \right\|_1}{\rm{ s}}{\rm{.t}}{\rm{. }}{\left\| {\Theta \mathit{\Psi} m-y} \right\|_2} < \varepsilon $ | (1) |
(1) 式中m是重建的二维谱，y是测量到的二维时域数据，ε往往低于噪声水平的期望值.
(1) 式中的目标函数是l1范数，定义为
在本文中，除直接使用频谱进行稀疏性约束外，也使用频谱的全变分(total variation)对频谱图的连续性进行约束.因此(1)式中的目标函数可以改写为：
| $ \hat m = \min \left[{||m|{|_1}{\rm{ + }}\lambda {}_{\rm{1}}T{V_{{f_{\rm{1}}}}}(m) + \lambda {}_{\rm{2}}T{V_{{f_{\rm{2}}}}}(m)} \right]{\rm{ s}}{\rm{.t}}{\rm{. }}||\Theta \mathit{\Psi} m -y|{|_2} < \varepsilon $ | (2) |
式中
本文采用变密度伪随机法在t1方向欠采样，用有限差分(TV)的l1范数作为约束.由于在t1方向欠采样，因此该方向上产生很多震荡信号，导致谱图分辨率低，尤其不易分辨出弱峰.因此我们将TV项的水平方向和竖直方向设置不同的权重，分别取舍，取得更好的效果.
1.2 伪随机采样实现机制伪随机采样模式通常根据特定的概率密度函数，用蒙特卡洛方法产生[13].本文实验所用的概率密度函数如图 1所示，该函数可以分成两个部分：连续采样区域和采样概率衰减区域.连续采样区域内的概率设置为1，表示采集所有信号.由于FID信号的能量集中在演化时间的前期，所以在演化时间的前期连续采集信号，可以保证重建信号的质量.参数r决定了连续采样区域占t1维的比例，例如r=0.2表示对于t1维前20%的点完全采集.采样概率衰减区域中的每个点采样的几率，则随着t1的延长（即随着t1维中采样点的索引k的增大）而降低，本文中采用幂指数函数实现概率的衰减.整个概率密度函数可以表示为：
|
| 图 1 概率密度函数及采样矩阵. (a), (b)中横坐标表示t1维的采样点的索引k，纵坐标表示对应点的采样概率q(k)，范围在0~1之间，最大值为1.ξ为采样率，连续采样区域内概率为1.连续采样区域占t1维总的采样点数的比率r为20%，(a), (b)使用的p值分别为3和5. (c), (d)分别是对应(a), (b)产生的采样矩阵.显然，在间接维采样概率大的采集的越多 Fig. 1 Probability density function and sampling matrix. (a), (b) The x-coordinate represents the index k in t1 dimension, and the y-coordinate gives the corresponding sampling probability q(k) with the range from 0 to 1. The ξ is sampling rate. The uniform sampling region rate r is 20% of the number N1 of t1 dimension. The probability inside the uniform sampling region is one. p values of 3 and 5 are used for (a), (b) respectively. (c), (d) Sampling matrix generated from (a), (b) respectively. Obviously, the bigger the sampling probability is, the more points are sampled |
| $ q(k) = \left\{ \begin{array}{l} 1, {\rm{ }}k \in \left[{1, \left\lfloor {r \times {N_1}} \right\rfloor } \right]\\ {(1 - {d_k})^p}, {\rm{ }}k \in \left[{\left\lfloor {r \times {N_1}} \right\rfloor + 1, {N_1}} \right] \end{array} \right. $ | (3) |
(3) 式中q(k)表示k点对应的采样概率，N1是t1维全采样的总点数，
根据上文所述的概率密度函数的概率确定随机采样点就可以产生伪随机采样模式.因为产生过程的随机性，所以产生的伪随机采样模式的非相干性并不一定好.为了避免产生相干性过高的采样模式，可重复生成伪随机采样模式并测量采样模式的峰值干涉值（即变换到傅立叶域，若有明显的非零值，说明所采的点有特定的周期性，即在一定意义上并不是随机的），选择峰值干涉值最小的结果作为最终的采样模式.
2 实验部分与结果分析 2.1 实验条件图 2(a)所示为29Si DQ-SQ的时域全采的模图数据，采样数据点阵t2×t1=384×136 (68个t1复数点)，完全黑色的部分是Bruker谱仪自动填零的数据.其FT2D频域谱图为图 2(b)所示，显示的是两个维度的实部.该谱采于Bruker AVANCEI 400 MHz谱仪，扫描次数为256，重复时间为3 s.图(c)表示在直接维上第194个t2点对应的间接维切片谱图，位置如图(b)中所画直线.
|
| 图 2 29SiDQ-SQ的时域数据(a)及其FT2D频谱图(b). (c)直接维第194个点对应的间接维的切片，即(b)图中直线所在位置2014-04-11T11:15">图中直线所在位置 Fig. 2 (a) Time domain data of 29Si DQ-SQ; (b) FT2D spectrum. (c) The spectrum of 194 in directive dimension, is that line in (b) |
t1维共68个复数点（实虚共136个点），我们采集24个复数点，即采样率ξ约为35%.在采样率ξ为35%的前提下，比较不同连续采样区域比例r对硅谱重建结果的影响，如r为t1维全采总点数的10%，16%，20%，25%，30%，35%，如图 3所示，图 3(a), (d)为采样矩阵（自上而下为1到N1），图 3(b), (e)为直接傅里叶变换FT2D的结果及直接维上第194个点的间接维切片，图 3(c), (f)为本文算法的重建结果及直接维上第194个点的间接维切片，虚线表示用于显示图 3(e), (f)的等高线.其中λ1为0.1，λ2为0.05，迭代次数为14.从二维谱和一维切片的图示中，可以看出r为20%的时候，谱图质量已达到要求，几乎没有干扰信号；r为35%时，重建效果最好.我们还重建了峰更窄且分辨率更高的13C DQ-SQ谱图，也发现当r等于采样率时，重建效果更佳（见附录）.为了将来讨论的方便，我们将连续采样区域的宽度等于采样率的采样模式称为t1截尾(t1-cutoff, CUO)模式.
|
| 图 3 伪随机采样模式中不同连续采样区域比例r对应谱图的比较，采样率ξ均为35%.图中(a, d)为采样模式（白色表示采样，黑色表示不采样），(b, e)为FT2D变换结果及其直接维上第194个点的间接维切片，(c), (f)为本文算法的重建结果及其直接维上第194个点的间接维切片；从上到下，r为全采总数据的10%，16%，20%，25%，30%，35%.图中虚线表示用于显示的等高线，虚线以下是允许的噪声，虚线以上是信号或干扰信号 Fig. 3 Comparison of different uniform sampling region r in pseudo random sampling schemes. Sampling rate ξ are 35% for all spectra (24 out of total 68 t1 points). (a) Sampling matrixes (lines sampled are indicated by white lines while lines not sampled are indicated by black lines); (b, e) results of FT2D and the spectrum of 194 in directive dimension; (c), (f) results of our algorithm based on CS and the spectrum of 194 in directive dimension. From top to bottom, the r are 10%, 16%, 20%, 25%, 30%, 35% of total data respectively. Dashed lines stand for contour levels for show, and it is noise below the line and signals or interference signals above the line |
在此我们比较伪随机(pseudo-random sampling，PRAS)、t1截尾(t1-cutoff，CUO)模式、全随机(random sampling，RAS)和e指数非均匀采样(exponential-decay non-uniform sampling，expNUS)模式.其中全随机等价于伪随机里使用连续采样区域的比例r=0的模式.采样率ξ为35%的RAS模式和expNUS模式的相关重建数据图见图 4，圈中表示丢失的信号.对比图 3和图 4可知，PRAS和CUO模式优于RAS和expNUS模式，且CUO模式稍微更好.图 5为采样率为22%时的相关数据的比较.
|
| 图 4 35%采样率时e指数非均匀采样(上)与全随机采样(下)模式的比较.(a)采样矩阵；(b)直接傅立叶变换的结果；(c)重建结果.圈中表示丢失的信号或失真的信号 Fig. 4 Comparison between two sampling schemes, namely, exponential-decaying non-uniform sampling (expNUS) (up), and random sampling(RAS) (down), with the sampling factor x of 35%. (a) Sampling matrix; (b) results of FT2D; (c) reconstructed results. Lost signal and distorted signal are in the circle |
|
| 图 5 采样率ξ为22%时伪随机、t1截尾、e指数非均匀、全随机(从上至下)等4种采样模式的比较.(a)采样矩阵；(b)直接傅立叶变换的结果；(c) CS重建的结果.由于采样率降低，伪随机重建结果中存有少量干扰信号 Fig. 5 Comparison between four sampling schemes, namely pseudo-random sampling (PRAS), t1 cutoff sampling (t1-cutoff), exponential-decaying non-uniform sampling (expNUS), random sampling (RAS), with the sampling rate ξ of 22% (from top to bottom). (a) Sampling matrix used; (b) results of FT2D; (c) results of CS reconstruction.There are several interference signals in PRAS result because of low sampling rate |
由图 4和图 5得出，相对于e指数非均匀采样和完全随机采样，伪随机采样应用本文算法能更好地恢复出原始信号.这是因为核磁共振谱图FID开始的信号最强，信噪比最好，且决定着频谱的基本分布，所以演化前期的信号对重建更敏感，多采更有利于谱图重建.
3 结论压缩感知技术之所以应用广泛，是因为它能有效地节省扫描时间，并且能不失真地恢复真实信号.本文将其用于固体DQ-SQ的重建，使用了有限差分约束的l1范数，即有限差分项的水平和竖直方向的差异分别加了不同的权重进行约束.我们对伪随机采样、全随机采样和e指数采样进行了比较，发现伪随机采样表现出更好的重建结果.另外我们发现对于固体DQ-SQ谱图，连续采样区域的比例等于采样率时的CS重建效果最好.我们的发现有助于CS重建在计算机上的应用，即可以使用传统的FT2D采样模式，而不需要在t1维跳跃采样；在这种条件下，我们可以每次采一个t1点，即可以进行CS重建，直到谱图基本稳定为止.
附录图S1(a)为13C DQ-SQ时域全采的模图数据，采样数据点阵t2×t1=320×500 (250个t1复数点).其FT2D频域谱图为图S1(b)所示，显示的是两个维度的模图.该谱采于Bruker AVANCEI 900 MHz谱仪，扫描次数为64，重复时间为1.5 s.
|
| 图 S1 13C DQ-SQ的时域数据(a)及其FT2D频谱图(b) Fig. S1 (a) Time domain data of 13C DQ-SQ; (b) FT2D spectrum |
图S2为关于碳谱的不同连续采样区域比例的相关实验，同样是采样率ξ为35%，半径为全采总点数的6%，15.6%，25.6%，35%时，对应的CS重建结果.同样验证了本文所得的结论.
|
| 图 S2 随机采样模式中不同连续采样区域的比较.采样率ξ均为35%，图中(a)为采样矩阵，(b)为FT2D变换结果，(c)为CS的重建结果；从上到下，连续采样区域比例r为全采总数据的6%，15.6%，25.6%，35% Fig. S2 Comparison of different uniform sampling region rin pseudo random sampling schemes. Sampling factor ξ are 35% for all spectra (24 out of total 68 t1 points). (a) Sampling matrix used; (b) results of FT2D; (c) results of CS (compressed sensing). From top to bottom, the uniform sampling regions occupy 6%, 15.6%, 25.6%, 35% of the total data respectively |
| [1] | Jaravine V I I, Orekhov V Y . Removal of a time barrier for high-resolution multidimensional NMR spectroscopy[J]. Nature Methods , 2006, 3 (8) : 605-607 DOI:10.1038/nmeth900 |
| [2] | Bruschweiler R . Theory of covariance nuclear magnetic resonance spectroscopy[J]. J Chem Phys , 2004, 121 (1) : 409-414 DOI:10.1063/1.1755652 |
| [3] | Hu B W, Zhou P, Noda I et al . An NMR approach applicable to biomolecular structrue characterization[J]. Anal Chem , 2005, 77 (23) : 7534-7538 DOI:10.1021/ac051061o |
| [4] | Barna J C J, Laue E D, Mayger M R et al . Exponential sampling, an alternative method for sampling in two-dimemsional NMR experiments[J]. J Magn Reson , 1987, 73 (1) : 69-77 |
| [5] | Jeffrey C, Hoch A S S . Maximum entropy reconstruction, spectrum analysis and deconvolution in multidimensional nuclear magnetic resonance[J]. Method Enzymol , 2002, 338 : 159-178 DOI:10.1016/S0076-6879(02)38219-3 |
| [6] | Kazimierczuk K, Orekhov V Y . Accelerated NMR spectroscopy by using compressed sensing[J]. Angew Chem Int Ed , 2011, 50 (24) : 5 556-5 559 DOI:10.1002/anie.201100370 |
| [7] | Coggins B E, Zhou P . High resolution 4-D spectroscopy with sparse concentric shell sampling and FFT-CLEAN[J]. J Biomol NMR , 2008, 42 (4) : 225-239 DOI:10.1007/s10858-008-9275-x |
| [8] | Kazimierczuk K, Kozminski W, Zhukov I . Two-dimensional Fourier transform of arbitrarily sampled NMR data sets[J]. J Magn Reson , 2006, 179 (2) : 323-328 DOI:10.1016/j.jmr.2006.02.001 |
| [9] | Kazimierczuk K, Zawadzka A, Kozminski W et al . Random sampling of evolution time space and fourier transform processing[J]. J Biomol NMR , 2006, 36 (3) : 157-168 DOI:10.1007/s10858-006-9077-y |
| [10] | Emmanuel J C, Justin R, Member L et al . Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[J]. IEEE T Inf Theory , 2006, 52 (2) : 489-509 DOI:10.1109/TIT.2005.862083 |
| [11] | David L, Donoho M . Ieee Compressed Sensing[J]. IEEE T Inf Theory , 2006, 52 (4) : 1289-1 306 DOI:10.1109/TIT.2006.871582 |
| [12] | Emmanuel C, Justin R, Tao T . Stable signal recovery from incomplete and inaccurate measurements[J]. Commun Pure Appl Math , 2006, 59 (8) : 1207-1 223 DOI:10.1002/(ISSN)1097-0312 |
| [13] | Lustig M, Donoho D, Pauly J M . Sparse MRI: The application of compressed sensing for rapid MR imaging[J]. Magn Reson Med , 2007, 58 (6) : 1182-1 195 DOI:10.1002/(ISSN)1522-2594 |
| [14] | Qu X B, Guo D, Cao X et al . Reconstruction of self-sparse 2D NMR spectra from undersampled data in the indirect dimension[J]. Sensors (Basel) , 2011, 11 (9) : 8 888-8 909 |
| [15] | Kutyniok G . Theory and applications of compressed sensing[J]. GAMM-Mitteilungen , 2013, 36 (1) : 79-101 DOI:10.1002/gamm.v36.1 |
| [16] | Lustig M, Donoho D L, Santos J M et al . Compressed sensing MRI[J]. IEEE Signal Proc Mag , 2008, 25 (2) : 72-82 DOI:10.1109/MSP.2007.914728 |
| [17] | Liu Y P, Wan Q . Anti-sampling-distortion compressive wideband spectrum sensing for cognitive radio[J]. International Journal of Mobile Communications(IJMC) , 2011, 9 (6) : 604-618 DOI:10.1504/IJMC.2011.042779 |
| [18] | Yu S, Khwaja A Shaharyar, Ma J . Compressed sensing of complex-valued data[J]. Signal Process , 2012, 92 (2) : 357-362 DOI:10.1016/j.sigpro.2011.07.022 |
| [19] | Smith D S, Gore J C, Yankeelov T E et al . Real-time compressive sensing MRI reconstruction using GPU computing and split bregman methods[J]. Int J Biomed Imaging , 2012, 2012 (864827) |
| [20] | Hu S, Lustig M, Chen A P et al . Compressed sensing for resolution enhancement of hyperpolarized 13C flyback 3D-MRSI[J]. J Magn Reson , 2008, 192 (2) : 258-264 DOI:10.1016/j.jmr.2008.03.003 |
| [21] | Duarte M F, Davenport M A, Takhar D et al . Single-pixel imaging via compressive sampling[J]. IEEE Signal Process Mag , 2008, 25 (2) : 83-91 DOI:10.1109/MSP.2007.914730 |
| [22] | Li Shu-Tao(李树涛), Wei Dan(魏丹). A survey on compressive sensing[J]. Acta Automatica Sinica(自动化学报), 2009, 35(11): 1 369-1 377. |
| [23] | Holland D J, Bostock M J, Gladden L F et al . Fast multidimensional NMR spectroscopy using compressed sensing[J]. Angew Chem Int Ed , 2011, 50 (29) : 6 548-6 551 DOI:10.1002/anie.v50.29 |
| [24] | Bostock M J, Holland D J, Nietlispach D . Compressed sensing reconstruction of undersampled 3D NOESY spectra: application to large membrane proteins[J]. J Biomol NMR , 2012, 54 (1) : 15-32 DOI:10.1007/s10858-012-9643-4 |
| [25] | Maciejewski M W, Qui H Z, Rujan I et al . Nonuniform sampling and spectral aliasing[J]. J Magn Reson , 2009, 199 (1) : 88-93 DOI:10.1016/j.jmr.2009.04.006 |