一种基于中值滤波和统计学的拉曼光谱spike剔除算法

叶瑞乾，何 浩，郑 鹏，徐梦溪，王 磊*

1.厦门大学航空航天学院，福建 厦门 361101 2.厦门大学化学化工学院，固体表面物理化学国家重点实验室，能源材料化学协同创新中心(iChEM)，福建 厦门 361101

引 言

拉曼光谱能原位地实时探测多物质的化学组分[1]，被广泛应用于细胞生物和电化学[2]、食品安全[3]、临床医学[4-5]等。然而，绝大多数分子的拉曼散射效应极其微弱，因而拉曼光谱中经常包含大量噪声。其中spike是一种在随机位置上出现的宽度极窄、峰高很强的单向尖峰信号[6]。一般认为spike的产生是因为宇宙射线被拉曼光谱仪中广泛使用的电荷耦合器件(CCD)[7]所捕获，宇宙射线可以直接激发CCD像素，产生的spike的强度有时比正常的拉曼峰高出几个数量级，会显著降低光谱的对比度，必须在进一步分析之前被剔除。

目前已经有一些文献介绍不同的spike剔除方法，可以归结于两类，首先是基于硬件的实现。例如，Zhao[8]用直缝扩展狭缝高度以提升信号强度，并分析由此带来的图像弯曲问题，加以校正，将CCD图像分成多个条带，通过条带之间的比较可以识别和剔除spike；Huang[9]等将CCD高度所可以容纳的光纤打包到光纤束中，并以特殊的方法排列，以修复被spike损坏的图像。然而，复杂的仪器操作和额外成本的增加使基于硬件的spike剔除方法的应用较少。第二类剔除方法是基于软件的spike剔除方法。这些方法大致可以分为两种类型。第一种类型是使用滤波器对信号进行处理，包括有鲁棒性的平滑滤波器[10]、加权移动窗口滤波器[11]、小波变换[6，12-13]等。这些算法是基于以下假设，即与真实的拉曼峰相比，spike的带宽更小且强度更高。然而，在实际应用中，一些反映真实分子信息的窄峰会被误认为是spike。此外，使用滤波器可能会让真实的拉曼光谱失真。第二种是基于概率统计原理的剔除方法，包括上限频谱数据矩阵[14]、鲁棒求和方法[15]、基于分布设置阈值[16-17]和差值谱图检测方法[18]。这些方法成立的理论基础是在两个连续光谱中的同一位置发生spike的概率非常小。由于需要采集多次，因此，这些方法需要较长的处理时间来搜索spike，而且参数的选择也很复杂。同时，这两类spike剔除方法的处理对象是基于某条特定的光谱，对于有成千上万条光谱的拉曼成像数据立方体，其应用效率较低且保真性较差。随着拉曼成像技术的发展，针对大量光谱数据的批量spike剔除方法是未来的发展要求。

因此，本文提出一种实用的spike剔除方法，其原理是使用中值滤波器处理原始数据，并对中值滤波后的数据和原始数据的偏差作统计分析来寻找spike，并只替换这些spike值，以保证数据的高保真度，中值滤波的参数设置较于上述方法更加简单且有效。同时该算法有大量数据的处理能力，不仅适用于一维拉曼谱和二维拉曼图像的spike剔除，还适用于拉曼光谱成像得到的三维数据立方体。实验证明，当spike强度超过最高峰值强度20%时，该算法可以准确地识别spike，灵敏度为99.5%。与目前spike识别方法[14，19]相比，该算法具有较高的灵敏度和较小的检测极限，为拉曼光谱分析人员进行后续数据降维及机器学习算法等应用提供了一个重要的光谱预处理工具。

1 实验部分

1.1 通过中值滤波器和统计分析剔除spike

spike可以定义为一维、二维或者三维数据内相邻小区域的异常值。中值滤波器是以当前数据点为中心的小窗口的中值来代替每个数据点，其原理可以用来剔除这些异常值，但是非spike的信号也将随着滤波操作而失真。因此，我们开发了一种分两步进行的spike剔除方法，第一步是定位spike，第二步是只对被定位的spike执行中值滤波。这种策略可以较好地去除spike，并保留原始信号的高保真度。定位spike是其中的关键，可以根据中值滤波后的数据和原始数据作偏差，通过设置阈值来识别spike。将原始数据表示为M，滤波后的数据表示为N，偏差表示为E(E=M-N)。E的分布可以近似表示为高斯函数[式(1)]

(1)

式(1)中，x是随机变量，f(x)是概率密度，μ是期望，σ是标准差。

根据3-σ规则，99.7%的数据在以期望μ为中心, 标准方差σ的3倍区域内。这个特征可以用来识别spike，因为spike通常稀疏地分布在测量数据中，并且具有比正常信号高得多的强度。这些spike可以通过根据E的分布性质设置一个阈值来识别。该阈值将把E中大部分较小的数据点识别为正常值，而当E中的一个数据点的值超过该阈值时，它将被识别为spike。如果原始数据M中没有spike，则E中的大部分值将非常小，拟合分布的期望将接近0。但如果spike存在，使用式(1)直接拟合，拟合曲线将因为异常值的原因偏离它的真实分布。为解决这个问题，本文按升序排序E，并使用中间数据而不是使用E中的所有值来拟合分布。我们使用分位数方法将E切成一个系列连续区域。分位数的函数可以被描述为式(2)

τ=P(X≤xτ)=F(xτ)

(2)

式(2)中，X是随机变量，F(x)是连续分布函数，P是概率，τ是分位数。分位数指满足条件P(X

根据分位数的概念将E中的数据点以概率的方式划分为100个份，通过选择合适的τ，取中间的99个部分作为拟合分布所用的数据，中间99%的数据已剔除强度很大的spike和测量结果为负的CCD坏点，可以表现真实拉曼信号特征，经过实验测试，取99%以上的中间数据会引入个别异常值，取99%以下的数据，对拟合的分布提升不明显，且可能会丢失一些真实信号。原始数据和经过排序后的中间99%数据的拟合分布结果如图1所示。

图1 原始数据和中间99%的数据的分布

图1是使用全部数据及中间99%的数据做拟合的概率密度函数图，其中使用99%的数据的分布更瘦更高，意义为数据更加集中，异常值和正常值的距离更远，且大部分的偏差E在0左右，99%的数据拟合后左右两端比使用全部数据的拟合更窄，意味着正常值的数据更加集中，阈值的设置更加简单。

通过基于概率密度函数(PDF)的原理，为E的真实分布设置一个阈值，可以很容易地识别spike。该原理可以表示为式(3)

(3)

式(3)中，X是随机变量，x是实数，P为概率，f(t)为概率密度函数。

PDF函数可以通过设定一个数来获取数据超过该数的概率，但由这种直接方法设置的阈值将有两个缺陷。一是阈值选择的难度：关于spike和正常值之间的距离尚不清楚，因此没有参考依据。另一个缺陷是缺乏自动化和通用性：针对不同的拉曼数据需要设置不同的阈值，这对分析者来说很不方便且设置的值与使用者的经验有关。为了解决这些问题，我们采用基于逆累积分布函数(ICDF)的反向思想来直接设置概率阈值，通过设定某个概率P(X≥x)，可以根据该分布函数自动确定一个阈值x，数据超过该值为满足这个概率(X≥x)。因此，在不同的拉曼数据中将根据概率自动调整不同的绝对阈值x以剔除spike。为了方便起见，本文以下部分中设置的阈值指的均是概率阈值P。

在通过阈值设置检测到spike之后，下一步是用表示真实拉曼特征的值来替换原始数据M中的spike。最简单、直接的方法是用线性插值来填充空位置。但这种方法引入了伪影，并会因拉曼数据的复杂性而导致错误的结果。应用非线性函数进行插值可以使替换完成的拉曼数据更逼近真实数据。然而，一些非线性函数会引入更多的变量和未知的参数，这将使这个过程变得极其复杂，错误的参数甚至会降低拉曼特征的准确性。一种简单有效的方法是使用中值滤波后的数据N中的值代替原始数据M中的spike。中值将当前像素中的信号与前后像素连接起来，确保该值的可靠性。此外，它仅需要少量的参数来确定，易于操作。

整个算法的基本思想是在偏差分布E中找到超过阈值的数据点，并将其标记为spike。然后去除原始数据M中的这些spike，并填充对应位置的滤波后的N值，具体步骤如表1所示。拉曼数据M可以是一维拉曼光谱、二维拉曼图像和三维拉曼数据立方体，取决于维度D。Median(M,W)表示对数据M应用窗口大小为W的中值滤波器，Fit(E)代表对分布E作曲线拟合，M(I)是数据M在位置I的索引值，I可以是(i)，表示一维数据，(i,j)表示二维数据，(i,j,k)表示三维数据。

表1 算法流程图

1.2 参数优化

参数优化包括中值滤波器的窗口大小的选择和阈值的设定。这些参数应根据拉曼数据的噪声水平、强度和锐度的差异来选择，以在不影响真实拉曼特征的情况下准确地去除spike。本文针对三种维度下不同窗口尺寸和不同阈值进行实验。在此，以一维拉曼谱为例，演示如何选择合适的参数。拉曼光谱中不同参数下的结果如图2所示。

图2a为未经spike剔除算法的一维拉曼光谱数据，使用较大的窗口，有可能无法完全剔除图中的spike(如图2b所示)，或对图中的真实拉曼峰造成影响(如图2d)，使用较小的窗口尺寸的结果如图2e所示，由于spike窄带宽，高强度特征，小窗口识别更成功。因此，我们选择3作为一维拉曼光谱的窗口尺寸。

图2 不同窗口尺寸和阈值的一维spike剔除

在实验中，我们设置阈值1×10-3和1×10-5并进行比较。尽管阈值相差两个数量级但只有一个位于2 000波数附近的spike没有被检测到。引入分位数的方法确保了拟合分布是使用正常值，扩大了正常值和spike的差异，让阈值选择空间足够大。

二维拉曼图像和三维拉曼数据立方体的参数选择与上述操作相似。对于二维图像，将使用二维中值滤波器作用于原始数据。从一维拉曼光谱到二维拉曼图像的扩展不仅仅是从行扩展到列的更改，这涉及到滤波器在不同对象上的操作。在一维数据中，对不同波数下的信号强度进行滤波，是光谱维度上的滤波。在二维中，同一波数上的不同像素信号被滤波，对于采集过程，属于空间上的滤波。原始数据与中值滤波后的数据之间的偏差是二维分布，因此设置的阈值也是二维数据，对应于图像的行和列。为了简化拟合分布，我们首先将二维偏差数据重组为一维偏差数据。对扩展的一维数据的分析与前面描述的一维拉曼光谱数据相同。本文选择3×3作为窗口大小，1×10-6作为设定的阈值用来剔除二维拉曼图像中的spike。

对于三维拉曼数据立方体，数据的中值滤波器和偏差也将是三维的。三维数据的扩展既需要考虑拉曼数据采集的时间，也需要考虑不同的波数，首先是在深度方向上，然后是在宽度方向上，最后是在高度方向上。这将极大增加数据的数量，并使偏差分布更接近真实分布，也更准确。

根据参数的选择，对于三维数据立方体，信号的变化受到多维的限制。因此，一个小的滤波窗口足以去除spike，也可以减少程序的处理时间。由于数据量的扩展，阈值设置为1×10-6。

2 结果与讨论

为了测试该算法在实际应用中的性能，使用真实的实验数据集进行了测试。我们使用Raman 11线扫描快速成像共聚焦拉曼光谱仪器(Nanophoton, Japan)用于扫描生物细胞。该实验总共收集100线的光谱，每线包含400条光谱，光谱深度为1 340个波数。收集到的数据可以以完整光谱的形式进行区分，共获得40 000条光谱，每条光谱包含1 340个波数。如果使用每次采集的光谱数和采集的数量作为区分条件，可以获得1 340张大小为100×400的拉曼图像。此外，我们还可以根据采集深度、采集宽度和采集线数来构建一个大小为1 340×400×100的三维拉曼数据立方体。

为了更好地分析该算法的有效性，我们还在真实的实验数据中引入了不同数量的、不同强度的人工spike来进行验证实验。由于一维拉曼光谱、二维拉曼图像和三维拉曼数据立方体原理的不同，我们在不同的条件下进行了实验。以下分析中列出了具体的实验参数。

2.1 一维spike剔除

一维spike剔除的处理结果如图3所示。Spike出现概率较小，但是强度很大，远超过其他波峰，该算法很好地剔除了spike，同时完整保留了1 000 cm-1处的强拉曼峰(细胞成分：苯丙氨酸)。

图3 一维拉曼光谱spike剔除

作为验证，测试了20个光谱，其中有10个手动添加的人工spike。增加的spike位置是随机的，强度为峰值的0.1～0.4倍。实验结果见图4。

图4 一维拉曼光谱(加人工spike)

用该算法进行了10次独立的验证实验，并将检测结果进行平均，最终结果如表2所示。其中Ratios表示峰值的比例(0.1～0.4倍)，Predicted表示算法所预测出来的spike数目，True为真实spike数量，Percentage是算法的准确度。

表2 一维拉曼光谱spike检测

实验结果所示，在spike强度较低情况下(峰值的比值较低)，spike剔除算法的性能较差，因为人工spike可能由于其随机位置而落入信号的低处，使得叠加的信号太小，无法被检测到。当该比率增加到40%时，检测灵敏度可达到99.5%。参考文献[15]的检测灵敏度为97.5%，参考文献[19]中检测超过峰值50%的spike识别算法，灵敏度可以达到99%，与它们相比，本文介绍的方法可以检测强度更弱的spike，且灵敏度有较大提升。

2.2 二维spike剔除

图5(a)和(b)为波数在2 940 cm-1的真实拉曼图像的二维spike剔除处理结果。其中Y轴是线数，X轴是光谱条数，构成100×400的拉曼图像。颜色代表相对值的大小，蓝色代表小数值，黄色代表高数值。未剔除spike前，高强度的spike使图像其他位置的特征对比度降低。在剔除spike后，严重污染的原始拉曼图像可以显示图像的特征。

图5 二维拉曼图像spike剔除

对于二维验证实验，我们从1 340张图像中随机选择20张，并在每张图像中选择10个随机位置设置spike。然后将每个选定位置的值设定为原本值加上当前图像最大像素值的0.1到0.3倍作为spike。

通过增加不同的比率，我们得到了不同的检测结果。为减少验证实验的随机误差，进行了10个独立验证实验，将检测结果进行平均，最终结果如表3所示。

表3 二维拉曼图像spike检测

对于二维spike剔除方法，当spike强度超过峰值的30%时，可以有效地检测spike。与参考文献[19]中的二维spike检测算法相比，本文介绍的方法可以在spike强度较小的情况下检测和剔除spike。

二维spike检测灵敏度高于一维，通过二维spike剔除方法可以去除具有较低强度的spike。

2.3 三维spike剔除

三维拉曼数据立方体包含三维坐标信息和当前位置的值，它表示四维信息。在图中难以表示四维信息。因此，在使用了三维spike剔除方法后，我们可以简单地使用图6来显示三维拉曼数据立方体spike剔除的效果。

拉曼数据立方体可以表示为1 340张二维拉曼图像，或40 000条一维拉曼光谱。因此，用三维方法剔除spike后，可以直接剔除一维和二维的spike。图6中的立方体是指一个包含1 340个拉曼图像的真实样本数据立方体。图像上的点表示某条拉曼光谱，图的下部显示该光谱的spike剔除结果。图中选取三个数据立方体切片作为效果说明，分别对应波数在1 000，2 130和2 880 cm-1时的拉曼图像。

图6 三维拉曼数据立方体spike剔除

对于三维验证实验，我们随机选择20个波数上的图像，并在随机位置添加10个spike。spike的强度是拉曼数据立方体最高值的0.1～0.2倍。计算结果如表4所示。

表4 三维拉曼数据立方体spike检测

当强度大于峰值的20%时，可剔除所有spike。与一维和二维spike剔除相比，三维方法可以检测更低强度的spike，提高数据处理能力。这是一种新的针对大量拉曼数据的spike剔除方法，在其他文献中很少被提及。

灵敏度随着维度的增加而增加有两个原因。一方面，数据的特征更加明显。另一方面，维数的增加使数据量迅速增加，数据拟合出的分布更准确，spike剔除算法更灵敏。

3 结 论

随机的窄带宽、高强度的spike容易附着在拉曼光谱上导致信号的失真，对后续光谱分析造成困难。本文提出一种spike剔除算法，可以根据原始拉曼数据和中值滤波后数据之间的偏差拟合高斯分布并设置概率阈值来剔除spike。该方法不仅适用于一维拉曼光谱，而且也适用于二维拉曼图像和三维拉曼数据立方体，检测灵敏度将随着维度的增加而提高。当spike比峰值高20%时，该算法可以在保留原始信号的真实性的基础上达到99.5%的灵敏度。该算法阈值设置简单，适用不同维度的拉曼数据，能有效剔除低强度的spike，有望成为拉曼光谱分析人员重要的预处理工具。