基于奇异值分解和中位数绝对偏差的拉曼成像数据去噪方法
范贤光1,2,3, 吴腾达1, 支瑜亮1, 王昕1,2,3,*
1. 厦门大学航空航天学院仪器与电气系, 福建 厦门 361005
2. 传感技术福建省高等学校重点实验室, 福建 厦门 361005
3. 厦门市光电传感技术重点实验室, 福建 厦门 361005
*通讯联系人 e-mail: xinwang@xmu.edu.cn

作者简介: 范贤光, 1980年生, 厦门大学航空航天学院仪器与电气系教授 e-mail: fanxg@xmu.edu.cn

摘要

拉曼成像是一种无损伤、 无需标记的光谱成像技术, 它可以提供样品的不同组分的分子指纹信息以及空间分布特征, 相比其他成像技术有着更重要的应用。 但是拉曼散射的截面积小, 灵敏度低, 加上在很多实验中为了观察某些组分的动态分布而缩短扫描时间, 导致最终得到的成像数据被噪声干扰, 因此往往需要对信号进行去噪处理。 常规的算法一般都是基于一个给定的数学模型对光谱进行处理, 容易造成过滤波, 使得信号失真; 另外, 在处理拉曼成像数据时, 常规算法往往是对数据进行逐条光谱去噪, 从而忽略了多条光谱之间的相互关系, 导致最终的拉曼图像仍然受许多噪点干扰。 因此, 提出了一种基于奇异值分解和中位数绝对偏差的拉曼成像的信号处理方法, 用于拉曼成像数据的去噪处理。 该方法首先对拉曼成像数据进行奇异值分解, 获得一个奇异值矩阵与两个正交矩阵; 然后通过中位数绝对偏差法对奇异值矩阵中的各奇异值进行离群值检测, 选取前 k个被连续标记的离群值作为要保留的奇异值, 并将其余的奇异值赋值为零, 得到新的奇异值矩阵; 最后用新的奇异值矩阵与两个正交矩阵重新求解得到去噪后的拉曼成像数据。 实验中, 首先验证了中位数绝对偏差法确定前 k个奇异值的正确性, 其次分别从处理后的图像质量和信号波形两方面对比了该算法与常规算法的去噪效果。 结果证明, 中位数绝对偏差法可以快速地确定出合理的 k值大小, 而且, 依据该 k值处理后的成像数据不仅在成像质量上消除了大量的噪点, 使得组分的空间分布特征清晰可见, 也在信号波形上较完美地保留了微小谱峰, 并恢复光谱信号。 该算法不同于常规算法, 能同时对整个拉曼成像数据进行处理, 并保留光谱之间的统计特征, 是一种更加有效的拉曼成像数据的去噪方法。

关键词: 去噪; 拉曼成像; 奇异值分解; 中位数绝对偏差
中图分类号:O657.3 文献标志码:A
Denoising Method for Raman Imaging Data Based on Singular Value Decomposition and Median Absolute Deviation
FAN Xian-guang1,2,3, WU Teng-da1, ZHI Yu-liang1, WANG Xin1,2,3,*
1. Department of Instrumental and Electrical Engineering, School of Aerospace Engineering, Xiamen University, Xiamen 361005, China
2. Fujian Key Laboratory of Universities and Colleges for Transducer Technology, Xiamen 361005, China
3. Xiamen Key Laboratory of Optoelectronic Transducer Technology, Xiamen 361005, China
*Corresponding author
Abstract

Raman imaging is a noninvasive, marker-free spectral imaging technique that provides molecular fingerprinting and spatial distribution of different components of a sample, and is more important than other imaging techniques. However, the Raman scattering has a small cross-sectional area and low sensitivity. In addition, in many experiments, in order to observe the dynamic distribution of certain components, the scanning time is shortened, and the resulting imaging data are disturbed by noise, so it is often necessary to denoise the signal. Conventional algorithms generally process the spectrum based on a given mathematical model, which is likely to cause excessive filtering and distortion of the signal. In addition, when processing Raman imaging data, conventional algorithms tend to denoise the data one by one. This neglects the relationship between multiple spectra, resulting in the final Raman image still being disturbed by many noises. Therefore, a signal processing method based on singular value decomposition (SVD) and median absolute deviation (MAD) is proposed for denoising Raman imaging data. Firstly, the singular value decomposition is performed on the Raman imaging data to obtain a singular value matrix and two orthogonal matrices. Then, all singular values in the singular value matrix are detected by the median absolute deviation method. The consecutively labeled outliers are used as singular values to be preserved, and the remaining singular values are assigned to zero to obtain a new singular value matrix. Finally, the new singular value matrix and two orthogonal matrices are solved again to obtain a denoised Raman imaging data. In the experiment, we first verify the correctness of the median absolute deviation method in determining the k value, and then the proposed algorithm is compared with the conventional algorithm from the aspects of image quality and signal waveform. The results show that the median absolute deviation method can quickly determine a reasonable value, and the imaging data processed according to this value not only eliminate a lot of noise in the imaging quality, but also make the spatial distribution characteristics of the components clearly visible. The tiny peaks are also perfectly preserved on the signal waveform and the spectral signal is recovered. This algorithm is different from the conventional algorithm in that it can process the entire Raman imaging data at the same time and preserve the statistical features between the spectra. It is a more effective denoising method for Raman imaging data.

Keyword: Denoising; Raman imaging; Singular value decomposition; Median absolute deviation
引言

拉曼成像是一种建立在拉曼光谱之上的成像技术, 由于其能够良好地描述物质的空间分布, 兼具拉曼光谱的快速、 无创、 无损等特性, 在医学、 生物化学、 食品安全等领域得到了广泛的应用[1]。 高信噪比的拉曼成像数据是研究物质组分特征的前提。 由于拉曼散射的截面积小, 灵敏度低, 容易被噪声干扰, 而在很多实验中, 我们为了观察某些动态过程, 往往缩短扫描时间, 这使得获取到的拉曼成像数据信噪比不高, 导致最终的成像质量下降。 近年来, 国内外学者提出了各种算法用于拉曼光谱的去噪处理, 如: S-G滤波、 小波变换等。 S-G滤波是一种基于滑动窗口的平滑算法, 其基本思想是定义奇数大小的滑动窗口沿序列平移, 采用多项式拟合滑动窗口内的数据点, 从而达到去噪目的[2]。 小波变换则是一种时频分析算法, 其将带噪信号分解为低频部分和高频部分, 对高频部分采用阈值法扣除噪声, 最后经小波逆变换得到去噪后的信号[3, 4]

这些算法都是对单条光谱进行处理, 而在拉曼成像数据中, 多条光谱之间往往存在相互联系, 故现有的算法对这样的信号处理效果并不理想。 因此, 本文提出了一种适用于拉曼成像数据上的去噪方法。 该方法首先对拉曼成像数据进行奇异值分解(singular value decomposition, SVD), 得到一系列奇异值。 其次, 通过中位数绝对偏差法(median absolute deviation, MAD)对奇异值进行离群值检测, 确定出需要保留的前k个奇异值, 并扣除剩余的奇异值。 最后再用新的奇异值矩阵重新求解得到去噪后的拉曼成像数据。 该算法不仅去除了数据中的绝大部分噪声, 同时也大幅改善了成像质量。 实验中, 我们分别对比了该算法和S-G滤波、 小波变换等算法的去噪和成像效果。

1 理 论
1.1 奇异值分解

奇异值分解作为数值线性代数中最有效的工具之一, 已被广泛应用于统计分析、 物理和应用科学等领域[5]

在拉曼成像领域中, 拉曼成像矩阵D∈ ▯m× n可以描述为

D=S+E(1)

其中, S∈ ▯m× n为光谱信号矩阵, E∈ ▯m× n为误差矩阵。

根据奇异值分解的定义[6], 存在正交矩阵URm× mVRn× n, 使得成像矩阵

D=VT(2)

其中, Σ = Σ1OOORm× n, 且Σ 1=diag(σ 1, σ 2, …, σ r), 其对角元素按照顺序

σ1σ2σr> 0,  r=rank(D)(3)

排列。

由于UV均为正交矩阵, 满足UTU=ImVTV=In, 所以分别对式(2)左乘UT和右乘V, 并结合式(1), 不难得到

Σ=ΣS+ΣE(4)

其中, Σ S=UTSV, Σ E=UTEV

由式(4)可以看出, 对原始成像矩阵D的奇异值分解结果等价于分别对信号矩阵S与误差矩阵E的奇异值分解结果之和。 误差矩阵E通常是由噪声引起的, 噪声具有随机性, 由于奇异值分解本质上是一种线性变换, 因此对误差矩阵E奇异值分解的结果Σ E, 其对角线上的奇异值仍然是一系列噪声。 当成像矩阵拥有一定的信噪比时, 误差矩阵E的奇异值都很小, 并且都很接近, 而信号主要集中在大的奇异值上, 故可以通过保留Σ 中大奇异值, 扣除小奇异值, 再对成像矩阵进行重新求解来达到去噪的目的。

1.2 中位数绝对偏差

定义奇异值矩阵Σ 中的前k个对角元素σ 1, σ 2, …, σ k为要保留的奇异值, 而其余对角元素σ k+1, σ k+2, …, σ r为要扣除的奇异值。 合理的k值选择决定了去噪后的图像质量, 比较常见的方法是通过直接观察奇异值的大小来确定[7, 8, 9]。 该方法依靠人来配合, 并不能满足生产需求。 实际上, 在选择合理k值的过程中, 通过多次观察奇异值的分布会发现, σ k+1, σ k+2, …, σ r往往数量居多, 整体呈线性分布; 而σ 1, σ 2, …, σ k数量少, 与σ k+1, σ k+2, …, σ r的分布相距较远, 故在某种意义上可以认为要保留的奇异值σ 1, σ 2, …, σ k为拥有离群特征。 基于这一经验假设, 提出采用离群值检测的方法来发现奇异值中的离群值, 进而确定k值。

中位数绝对偏差法是一种非常简单精妙的离群值检测方法, 在统计学等领域有着重要的应用[10]。 已知奇异值矩阵Σ 中的奇异值序列σ 1, σ 2, …, σ r, 首先确定中位数绝对偏差[11]

MAD=1.4826×Mi(|σi-Mj(σj)|)(5)

其中, σ i, σ j∈ {σ 1, σ 2, …, σ r}, MiMj表示对序列取中位数。

分析序列σ 1, σ 2, …, σ r容易得到该序列的中位数σ M, 则离群值指示序列o1, o2, …, or可由式(6)的决策法则来确定

oi=1,  σi-σMMAD> |±3|0, σi-σMMAD|±3|(6)

我们只对前几个大的奇异值感兴趣, 故在离群值指示序列o1, o2, …, or中取第一条连续为1的子序列, 该子序列的长度即为k

1.3 基于奇异值分解和中位数绝对偏差的拉曼成像数据去噪方法

本文利用奇异值分解得到原始拉曼成像矩阵的一系列奇异值, 通过中位数绝对偏差法从中选取前k个大的奇异值进行重新求解即得到去噪后的成像矩阵, 提高成像质量。 具体算法步骤描述如下:

(1)扫描样品获得样品的拉曼成像矩阵D;

(2)对D进行奇异值分解, 得到U, Σ V;

(3)对Σ 中的奇异值使用中位数绝对偏差法作离群值检测, 得到离群值指示序列o1, o2, …, or;

(4)从左往右扫描o1, o2, …, or, 统计连续为1的个数, 当碰到0时终止扫描, 得到k;

(5)调整Σ , 保留前k个大奇异值, 其余赋值为0, 得到Σ * ;

(6)计算D* =* V, 最终得到去噪后的成像矩阵D*

2 实验部分

本文实验仪器为nanophoton公司生产的第三代显微拉曼成像系统Raman-11。 实验样品为Hela癌细胞, 采用线激光作为激发光, 每次扫描曝光时间为10 s。

3 结果与讨论

由于激发光能量低、 拉曼散射的截面积小等原因, 实验得到的拉曼成像数据往往带有很大的噪声, 这直接影响到成像质量。 为了解决该问题, 本文分别采用了提出的算法和常规算法对实测Hela癌细胞的拉曼光谱数据进行去噪处理。

图1所示为Hela癌细胞的平均光谱信号。 我们取苯丙氨酸的特征峰, 即波数1 002 cm-1处所对应的峰高进行成像, 来对比本文算法和常规算法的性能。

图1 Hela癌细胞的平均拉曼信号
蓝色虚线表示用于成像的波数位置
Fig.1 Average Raman signal of Hela cells
The blue dotted line indicates the selected wave number for imaging

表1所示为本文算法对Hela癌细胞数据进行处理后得到的前10个奇异值以及对应的离群值指示。 可以看到, 前4个奇异值被连续地标记为离群值, 因此保留前4个奇异值, 并将剩余的奇异值置为0, 来扣除成像数据的噪声。

表1 本文算法获得的前10个奇异值及对应的离群值指示 Table 1 Singular values and outlier indicators obtained by the proposed algorithm

为了进一步验证通过中位数绝对偏差法确定前k个奇异值的正确性, 分别取值为3, 4和5, 对比去噪后的成像数据中第一条拉曼光谱信号的保真情况, 如图2所示。 通过图2的处理结果可知, 当k为3时, 去噪后的光谱信号强度在部分区域(比如蓝色虚线矩形圈出的区域)发生了偏移; 而当k分别为4和5时, 二者去噪后的光谱信号基本吻合, 并且信号强度与原始光谱信号相比大体上得到了保留。 由于k为5及以上的奇异值比较小, 对成像数据的贡献不大, 因此认为k取4是合理的。

图2 原始成像数据与分别取k=3, 4, 5去噪后成像数据中的第一条拉曼光谱信号Fig.2 First Raman spectral signal in the original imaging data and the denoised imaging data with k=3, 4, and 5 respectively

图3给出了各算法处理后的成像结果及拉曼光谱。 图3(a)为Hela癌细胞的原始拉曼图像, 可以看到图中带有很多噪点。 图3(c), (e)和(g)分别为采用S-G滤波、 小波变换、 本文算法对(a)中数据进行去噪处理后的结果。 从成像结果中可以看出, 相比于前两种常规算法, 经本文算法处理得到的图像几乎不带任何噪点, 信号特征损失很少, 成像清晰度也得到了极大的提高。 继续观察经各算法处理后的拉曼图像中不同位置的拉曼光谱去噪情况。 图3(b)为(a)中数据在A和B两点处的拉曼光谱信号, 可以看到信号中带有一定的噪声。 图3(d), (f)和(h)分别为S-G滤波、 小波变换和本文算法对(b)中数据处理后的拉曼光谱信号。 可以看到, S-G滤波只抑制了部分噪声, 小波变换在去噪的同时也使一些原本平整的光谱信号产生波动, 而本文算法能够较好地把噪声扣除。

为了客观验证算法的去噪能力, 引入基于Kaiser窗口的信噪比估计算法[12], 并计算了A和B两点处原始光谱与各算法处理后光谱的信噪比, 如图4所示。 相比于原始光谱, 各算法处理后的光谱信噪比均有所提升, 其中, 本文算法对B点光谱处理得到的信噪比最大, 说明本文算法对弱信号的处理能力最好, 该结论从图3处理后的拉曼成像结果也能得到, 比如由本文算法处理的拉曼图像中蓝色区域噪点很少, 明显优于常规算法。 而在A点处, 本文算法处理得到的光谱信噪比介于S-G和小波变换之间, 直观的推论是本文算法在处理强信号的能力上有所欠缺, 然而分析图3拉曼成像结果会发现, 提出的算法对强信号的处理仍然很细腻, 并且在波形上保留了一些微小信号, 对照图1的平均拉曼光谱来看, 这些信号很有可能是我们所需要的。 下面我们将通过单独对比A点处各算法的处理结果来说明这一点。

图3 (a)1 002 cm-1处Hela细胞的原始拉曼图像; (b)图(a)中A和B两点的拉曼光谱; (c), (e)和(g)分别为由S-G滤波、 小波变换和提出的算法对图(a)处理后的拉曼图像; (d), (f)和(h)分别为图(c), (e)和(g)中A和B两点的拉曼光谱Fig.3 (a) Original Raman image of Hela cells at 1 002 cm-1; (b) Raman spectra at point A and B in (a); (c), (e) and (g) Denoised Raman image of (a) using S-G filter, wavelet and the proposed SVD algorithm, respectively; (d), (f) and (h) Raman spectra at point A and B of (c), (e) and (g), respectively

图4 原始光谱与分别由S-G滤波、 小波变换以及本文提出的算法处理后的光谱在A和B两点处的信噪比Fig.4 Signal-to-noise ratio of the original spectra and the spectra processed by S-G filtering, wavelet transform and the algorithm proposed in this paper at two points A and B respectively

图5所示为图3(a)中数据在A点处的原始拉曼光谱与各种算法处理后的结果对比。 注意到原始拉曼光谱在2 326 cm-1处(紫色虚线矩形圈出的区域)存在一个弱小谱峰。 该谱峰区域各算法的处理情况, 如图5所示。 通过图5紫色虚线矩形圈出的区域可以很容易看出, 由S-G滤波处理后的该区域内的谱峰隐约可见, 但仍受噪声影响; 小波变换直接丢失了该谱峰信号; 而本文算法不仅去除了噪声, 同时也较完整地保留了该区域的谱峰特征, 这也从侧面说明本文算法在处理弱小信号特征时优于前两种常规算法。

图5 点A处原始拉曼光谱与分别由S-G滤波、 小波变换以及本文提出的算法处理后的结果对比Fig.5 Comparison of original Raman spectra at point A with the results of S-G filtering, wavelet transform and the proposed algorithm

4 结 论

提出了一种基于奇异值分解和中位数绝对偏差的拉曼成像的信号处理方法, 用于拉曼成像数据的去噪处理, 提高了成像质量。 本文通过实验对该算法的性能进行了验证, 并与常规算法进行了比较。 事实证明, 使用基于奇异值分解和中位数绝对偏差的去噪方法处理得到的拉曼成像数据不仅在成像质量上消除了大量的噪点, 提高了清晰度, 同时在波形上也较完美地吻合了原始光谱信号。 因此, 该算法为处理带噪声的拉曼成像数据提供了一个强有力的工具。

参考文献
[1] HU Xiao-hong, ZHOU Jin-chi(胡晓红, 周金池). Analytical Instrumentation(分析仪器), 2011, 6: 1. [本文引用:1]
[2] Schafer R W. IEEE Signal Processing Magazine, 2011, 28(4): 111. [本文引用:1]
[3] ZHAO Ming-fu, TANG Ping, TANG Bin, et al(赵明富, 唐平, 汤斌, ). Spectroscopy and Spectral Analysis(光谱学与光谱分析), 2018, 38(3): 844. [本文引用:1]
[4] Srivastava M, Anderson C L, Freed J H. IEEE Access, 2016, 4: 3862. [本文引用:1]
[5] Ghaderyan P, Abbasi A, Saber S, et al. Expert Systems with Applications, 2018, 114: 428. [本文引用:1]
[6] ZHANG Xian-da(张贤达). Matrix Analysis and Applications(矩阵分析与应用). Beijing: Tsinghua University Press(北京: 清华大学出版社), 2013. 288. [本文引用:1]
[7] Ahmad F M, Shen R, Zaheer A B, et al. Meteorology and Atmospheric Physics, 2018, 130(6): 689. [本文引用:1]
[8] Feng Liang, Zhou Cangqi, Zhao Qianchuan. Physica A-Statistical Mechanics and Its Applications, 2019, 513: 424. [本文引用:1]
[9] Prakash N, Ramachand ran A, Varma R, et al. Analyst, 2018, 143(14): 3284. [本文引用:1]
[10] Aydin D. Wind and Structures, 2018, 26(6): 383. [本文引用:1]
[11] Leys C, Ley C, Klein O, et al. Journal of Experimental Social Psychology, 2013, 49(4): 764. [本文引用:1]
[12] Kumar K S, Yazdanpanah B, Kumar P R. International Conference on Communications and Signal Processing, 2015, 4: 157. [本文引用:1]