基于低秩Hankel重建的NMR波谱自动去噪法 (中文，English)

邱天予¹, Wenjing Liao², 黄奕晖¹, 吴瑾瑜¹, 郭迪³, 刘东宝¹, 王新¹, Jian-Feng Cai⁴, 胡炳文⁵, 屈小波^1,*

¹ 厦门大学，电子科学系，福建省等离子体与磁共振研究重点实验室，生物医学智能云研发中心，中国，厦门，361005.
² School of Mathematics, Georgia Institute of Technology, Atlanta, GA 30332 USA.
³ 厦门理工学院，计算机与信息工程学院，中国，厦门，361024.
⁴ Department of Mathematics, The Hong Kong University of Science and Technology, Hong Kong, China.
⁵ 华东师范大学，物理学院，中国，上海，200062.
^* Emails: quxiaobo <at> xmu.edu.cn or quxiaobo2009 <at> gmail.com

引用

Tianyu Qiu, Wenjing Liao, Yihui Huang, Jinyu Wu, Di Guo, Dongbao Liu, Xin Wang, Jian-Feng Cai, Bingwen Hu, Xiaobo Qu, An Automatic Denoising Method for NMR Spectroscopy Based on Low-Rank Hankel Model, IEEE Transactions on Instrumentation & Measurement, 70: 1-12, 2021.

摘要

一直以来，由于磁共振波谱的低灵敏度，磁共振去噪都是一个重要的问题。由于无噪FID信号的指数特性，去噪问题可以被认为是一个基于低秩Hankel矩阵的含衰减因子的指数信号正则化重建问题。在这一问题中，去噪结果主要取决于正则化参数的选择。本文基于利用核磁共振时域信号低秩特性的凸优化去噪方法，提供了参数界的理论估计，实现去噪参数的自动设置。

方法

1. 背景

在生物医学研究中，核磁共振(NMR) 波谱是一项重要的研究手段 [1-3]。然而，在实际应用中，NMR波谱往往被高斯噪声所污染。重复采样求平均是最常见的去噪方法，但这一方案使得采样时间大幅增加。
指数信号的低秩Hankel特性被广泛应用于NMR波谱去噪中，例如低秩Hankel重建法(LRHM) [4]。然而，如何选取一个合适的正则化参数仍然主要依靠个人经验。
本文基于LRHM，提出了一种凸优化去噪方法，通过对加权Hankel矩阵谱范数界的准确估计，实现参数的自动设置。该方法仅需要对噪声进行估计，便可以有效地计算出合适的参数值，并获得可信的去噪结果。

2. 基于LRHM的自适应正则化参数去噪(CHORD)

NMR波谱的时域信号 (FID)可以建模为R个衰减指数分量之和[4][5][10-14]

$$ x_{0}(t_{n})=\sum_{r=1}^{R}a_{r}e^{\left({j2{\pi}f_{r}-{\tau}_{r}}\right)t_{n}}, \ n=0,\ldots,2N, (1) $$ 这里$a_r$、$f_r$、$\tau_r$分别为信号幅度、中心频率和衰减因子。
本文所提方法 (CHORD)可以表示为对以下问题的求解 $$ \mathbf{\hat{x}}=\arg\min_{\mathbf{x} \in \mathbb{C}^{2N+1}} \left\|{\mathcal{R}\mathbf{x}}\right\|_{*}+\frac{\lambda}{2}\left\|{\mathbf{y}-\mathbf{x}}\right\|_{2}^{2}, (2) $$ 其中$\mathbf{y}{\in}\mathbb{C}^{(2N+1){\times}1}$表示包含高斯噪声z的观测信号；$\mathcal{R}$是将向量排成Hankel矩阵的算子；λ是正则化参数。
根据核范数的次梯度[6]以及近似理论，参数λ值应取到
$$ \lambda \leq {\frac{1}{\left(\left|\left\|\mathbf{Z}\right\|_{2}-\left\|\mathbf{\tilde{X}}\right\|_{2}\right|\right)}}, (3) $$ 其中$\mathbf{Z}$表示由$\mathbf{z}$转化所得的加权Hankel矩阵 $$ \mathbf{Z}= \left( \begin{array}{cccc} z_1 & \frac{z_2}{2} & \cdots & \frac{z_{N+1}}{N+1} \\ \frac{z_2}{2} & \frac{z_3}{3} & \cdots & \frac{z_{N+2}}{N} \\ \vdots & \vdots & \cdots & \vdots \\ \frac{z_{N+1}}{N+1} & \frac{z_{N+2}}{N} & \cdots & z_{2N+1} \end{array} \right), (4) $$ 类似地，$\mathbf{\tilde{X}}$是由$\mathbf{x}_{0}-\mathbf{\hat{x}}$转化所得的加权Hankel矩阵，其中$\mathbf{\hat{x}}$是问题(2)的解。
通过观察加权范数谱范数与噪声强度和矩阵规模之间的关系，我们发现，$\left\|{\mathbf{Z}}\right\|_2$与$\left\|{\mathbf{\tilde{X}}}\right\|_2$的均值基本与N无关 (见图1)。并且，该均值与噪声标准差(σ)呈线性关系。
在实际应用中，用户倾向于保留尽可能多的信号细节，因此，参数λ取值如下 $$ {\lambda}^{*}=\frac {1} {\left|\mathbb{E}\left\|\mathbf{Z}\right\|_2-\mathbb{E}\left\|\mathbf{\tilde{X}}\right\|_2\right|}. (5) $$

图1. 谱范数与噪声标准差(σ)的关系。(a) 100次蒙特卡洛实验中$\left\|{\mathbf{Z}}\right\|_2$与σ的关系；(b) 50次蒙特卡洛实验中$\left\|{\mathbf{\tilde{X}}}\right\|_2$与σ的关系。矩阵规模为(N+1)×(N+1) ，其中(N+1)分别为64、128、256、和512。曲线表示谱范数均值与标准差σ的关系，垂直误差线表示谱范数标准差。注：$\mathbf{x}_0$是由随机参数$a_r$、$f_r$、和$\tau_r$构成的衰减指数信号。$\mathbf{\hat{x}}$是CHORD所得去噪信号。参数λ是使得误差(NRMSE)取到最小的值。

由证明可得，$\mathbb{E}\left\|{\mathbf{Z}}\right\|_2$的上下界如下 $$ {\sigma}\frac{C(N+1)}{2N+1}\sqrt{R_{N}^{2}\left({1+\log\frac{R_{N}^{4}}{Q_{N}^{4}}}\right)}{\leq}\mathbb{E}\left\|{\mathbf{Z}}\right\|_{2}{\leq}{\sigma}\sqrt{2C_{\mathbf{w}}\log\left({2N+2}\right)}, (6) $$ 其中$C_{\mathbf{w}} = \max(\sum_{k=0}^{N}{w_{k}^{-2}},\sum_{k=1}^{N+1}{w_{k}^{-2}},\ldots,\sum_{k=N}^{2N}{w_{k}^{-2}})$；(4)种的权重$\mathbf{w}\in\mathbb{R}^{(2N+1)\times{1}}$定义为 $\mathbf{w}=\left[{ \begin{array}{ccccccc} 1 & 2 & \cdots & N+1 & \cdots & 2 & 1 \end{array} }\right]^{T}$。序列$R_{N}$和$Q_{N}$分别定义为$R_{N}^{2}=\sum_{k=0}^{2N}\left|{d_k}\right|^{2}$和$Q_{N}^{4}=\sum_{k=0}^{2N}\left|{d_k}\right|^{4}$，其中$ d_{k}= \left\{ \begin{array}{ll} \frac{2}{(k+1)(k+2)}\sum_{m=0}^{k}\frac{1}{m+1}, 0\leq k\leq N \\ \frac{2}{(2N-k+1)(k+2)}\sum_{m=k}^{2N}\frac{1}{m-N+1}, N< k\leq 2N \\ \end{array} \right. $。当N足够大时，上下界仅相差$\sqrt{\log N}$倍。
选择$\mathbb{E}\left\|\mathbf{Z}\right\|_{2}$的下界可以获得一个更大的$\lambda$值，这有益于保留更多的信号细节。因此在本文中，我们建议选择下界。对于其中的常数C，我们建议取$C=2.9$用于去噪。
本文采用不同的N、σ、λ、信号和噪声，通过蒙特卡罗实验，经验性估计$\left\|{\mathbf{\tilde{X}}}\right\|_2$的值。实验结果说明，σ是确定$\mathbb{E}\left\|{\mathbf{\tilde{X}}}\right\|_2$的主要因素。此外，$\left\|{\mathbf{\tilde{X}}}\right\|_2$的均值与σ呈线性关系。在去噪问题中，我们建议$\mathbb{E}\left\|\mathbf{\tilde{X}}\right\|_2=1.94σ$。

3. 主要结果

这里我们将CHORD用于一维NMR波谱的去噪，并将其结果与典型方法Cadzow [7], [8]与主流方法rQRd [7]相比较。
该实测一维$^1$H NMR波谱样品为肌酸、胆碱、柠檬酸镁和柠檬酸钙混合物，其浓度分别为0.03、0.03、0.06、0.06 g/mL（采样参数细节见原文）。
本文通过截取带噪FID信号的末端用以估计噪声强度。截断长度由Kolmogorov–Smirnov (KS) 估计 [9] 进行判断。
去噪结果展示在图2中。在相对强的噪声下(σ = 0.035)，Cadzow平滑了波谱，一方面，有效地抑制了噪声；另一方面，这也导致部分谱峰的丢失(例如6.8 ppm处的峰)。rQRd产生了一个包含明显噪声残留的去噪谱[橙线]，并且伴随着低强度谱峰的削弱(例如6.8 ppm处的峰)。CHORD能有效去除噪声，并保存更多的谱峰细节。

图3. 代谢谱去噪结果(σ = 0.035)。绿线表示高SNR参考谱；黑线表示带噪谱；蓝、橙和红线分别为Cadzow、rQRd和CHORD的去噪谱。注： Cadzow和rQRd的去噪结果对应最低NRMSE。棕色虚线代表去噪谱的MAE值。

下载

所提方法的MATLAB代码可以从这里下载。
本文也可以在这里下载。

致谢

这项工作得到了国家自然科学基金（6212200447、61971361、61871341、61811530021），国家重点研发计划（2017YFC0108700），厦门市科学技术计划（3502Z20183053），福建省卫生教育联合攻关计划项目（2019-WJ-31），和厦门大学南强拔尖人才计划的资助。

参考文献:

[1] A. Wegemann, et al., "A portable NMR spectrometer with a probe head combining RF and DC capabilities to generate pulsed-field gradients," IEEE Trans. Instrum. Meas., vol. 69, no. 10, pp. 8628–8636, Oct. 2020.
[2] O. Beckonert, et al., "High-resolution magic-angle-spinning NMR spectroscopy for metabolic profiling of intact tissues," Nature Protocols, vol. 5, no. 6, pp. 1019–1032, Jun. 2010.
[3] M. C. Preul, et al., "Accurate, noninvasive diagnosis of human brain tumors by using proton magnetic resonance spectroscopy," Nature Med., vol. 2, no. 3, pp. 323–325, Mar. 1996.
[4] X. Qu, et al., "Accelerated NMR spectroscopy with low-rank reconstruction," Angew. Chem. Int. Edit., vol. 54, no. 3, pp. 852–854, Jan. 2015.
[5] T. Qiu, et al., "Review and prospect: NMR spectroscopy denoising and reconstruction with low-rank Hankel matrices and tensors," Magn. Reson. Chem., vol. 59, pp. 324–345, Mar. 2021.
[6] J.-F. Cai, et al., "A singular value thresholding algorithm for matrix completion," SIAM J. Optim., vol. 20, no. 4, pp. 1956–1982, 2010.
[7] L. Chiron, et al., "Efficient denoising algorithms for large experimental datasets and their applications in Fourier transform ion cyclotron resonance mass spectrometry," Proc. Nat. Acad. Sci. USA, vol. 111, no. 4, pp. 1385–1390, Jan. 2014.
[8] J. Gillard, "Cadzow’s basic algorithm, alternating projections and singular spectrum analysis," Statist. Interface, vol. 3, no. 3, pp. 335–343, 2010.
[9] F. J. Massey Jr., "The Kolmogorov-Smirnov test for goodness of fit," J. Amer. Statist. Assoc., vol. 46, no. 253, pp. 68–78, 1951.
[10] J. Ying, et al., "Hankel matrix nuclear norm regularized tensor completion for N-dimensional exponential signals," IEEE Trans. Signal Proces., vol. 65, no. 14, pp. 5520-5533, 2017.
[11] J. Ying, et al., "Vandermonde factorization of Hankel matrix for complex exponential signal recovery—Application in fast NMR spectroscopy," IEEE Trans. Signal Proces., vol. 66, no. 21, pp. 5520–5533, 2018.
[12] H. Lu, et al., "Low rank enhanced matrix recovery of hybrid time and frequency data in fast magnetic resonance spectroscopy," IEEE Trans. Biomed. Eng., vol. 65, no. 4, pp. 809–820, 2017.
[13] D. Guo, et al., "A fast low rank Hankel matrix factorization reconstruction method for non-uniformly sampled magnetic resonance spectroscopy," IEEE Access, vol. 5, pp. 16033-16039, 2017.
[14] D. Chen, et. al., "Review and prospect: Deep learning in nuclear magnetic resonance spectroscopy," Chem. Eur. J., vol. 26, no. 46, p. 10, Aug. 2020.