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

范贤光， 吴腾达， 支瑜亮， 王 昕*

1. 厦门大学航空航天学院仪器与电气系， 福建 厦门 361005 2. 传感技术福建省高等学校重点实验室， 福建 厦门 361005 3. 厦门市光电传感技术重点实验室， 福建 厦门 361005

引 言

拉曼成像是一种建立在拉曼光谱之上的成像技术， 由于其能够良好地描述物质的空间分布， 兼具拉曼光谱的快速、 无创、 无损等特性， 在医学、 生物化学、 食品安全等领域得到了广泛的应用[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]， 存在正交矩阵U∈Rm×m和V∈Rn×n， 使得成像矩阵

D=DΣVT

(2)

σ1≥σ2≥…≥σr>0，r=rank(D)

(3)

排列。

由于U和V均为正交矩阵， 满足UTU=Im和VTV=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-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.482 6×Mi(|σi-Mj(σj)|)

(5)

其中，σi，σj∈{σ1，σ2， …，σr}，Mi和Mj表示对序列取中位数。

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

(6)

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

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

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

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

(2)对D进行奇异值分解， 得到U，Σ和V；

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

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

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

2 实验部分

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

3 结果与讨论

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

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

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

图1 Hela癌细胞的平均拉曼信号蓝色虚线表示用于成像的波数位置

The blue dotted line indicates the selected wave number for imaging

表1 本文算法获得的前10个奇异值及对应的离群值指示

为了进一步验证通过中位数绝对偏差法确定前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 withk=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 结 论

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