拉曼成像光谱数据中宇宙射线Spike的识别及修正方法

摘要

本发明公开一种拉曼成像光谱数据中宇宙射线Spike的识别及修正方法，不再直接用两个光谱的线性近似谱进行比较，而是利用可能包含Spike和小量尖峰信息的目标光谱ss进行比较，减少连续光谱相似度低时为判别Spike带来的干扰，提高了识别的正确率，能有效在拉曼成像数据中集中识别出含有Spike的拉曼光谱并修正。

说明书

技术领域

本发明涉及一种光谱成像处理方法，尤其是涉及一种拉曼成像光谱数据中宇宙射线Spike的识别及修正方法。

背景技术

拉曼成像技术能提供样品空间上光谱上的信息，有效地提取出样品成分的空间分布信息，被广泛应用生物组织、材料成分分析等领域。耦合电荷器件(CCD)是目前拉曼光谱仪的主要传感器件，具有高灵敏度和低噪声等优点，但却容易受到宇宙射线的干扰，在拉曼光谱上产生狭窄而尖锐的峰，称为Spike。由于宇宙射线的干扰是随机的，所以区别于一般的拉曼特征峰，Spike是随机出现的。

拉曼成像技术一次能扫描样品上数千个点，每个样品点都生成一条拉曼光谱，故拉曼成像的数据集包含了样品的数千条拉曼光谱，这些光谱全部拥有相同的波数。在这些光谱中，有一些光谱受到宇宙射线干扰，产生Spike，这些Spike通常强度值很高，若不去除，会对后续的数据处理造成影响。但由于拉曼成像的数据集中包含了太多的拉曼光谱，无法用人眼一一鉴别是否包含Spike，因此需要一种技术在拉曼成像数据集中自动识别包含Spike的光谱，并剔除这些Spike。

在拉曼光谱中去除Spike的现有技术可以分为两种，一种是只利用单一光谱的信息去除Spike的方法。另一种是需要利用多个光谱，通过比较光谱之间差异来识别、剔除Spike的方法。

单一拉曼光谱修正方法主要依靠Spike峰宽窄，强度高的特点，比较典型的有中值滤波，小波转换。这些方法基于Spike峰宽远小于正常拉曼峰宽的假设，而现有文献报道某些Spike的峰宽高达10个波数；另外中值滤波，小波转换这两种方法会对正常拉曼峰造成一定的扭曲。

如上文提及，Spike最重要的特征是其随机性，不可重现这一点，单一拉曼光谱修正Spike方法仅依靠峰宽和峰值强度作为判断Spike的依据，有将正常拉曼峰误判为Spike或遗漏识别Spike的可能——因为有可能出现峰宽狭窄的拉曼峰或峰值强度很小的Spike。这是单一拉曼光谱修正Spike方法的先天不足之处，单条光谱不能提供识别Spike的充分信息。

多拉曼光谱修正Spike方法，考虑到Spike出现的随机性，克服单一拉曼光谱修正方法缺点是比较两个相似的光谱，比如使用于连续光谱扫描的最近邻比较法，需要提供额外目标光谱的上界光谱法，一阶、二阶差分，利用矩阵理论，如PCA(principal componentanalysis)进行分析的方法。

其中最近邻比较法的主要问题是不能保证两个连续的光谱在拉曼峰位置足够相似，峰值过大的差异会造成误判，上界光谱法需要提供在相同条件下测量的目标光谱，显然不是一种实际的办法；基于PCA的方案则要保证光谱矩阵的主要信息不被丢失，但实际上各个光谱主成分的关系复杂，在去除噪声和Spike时，难免会损失一部分有用的信息。

发明内容

本发明的第一目的在于提供一种拉曼成像的光谱数据中宇宙射线Spike的识别方法，能提高识别正确率。

本发明的第二个目的在于提供一种从拉曼成像的光谱数据中宇宙射线Spike的修正方法。可以正确识别出整个拉曼成像光谱中含有的全部宇宙射线Spike并去除。

本发明第一目的可通过以下技术措施来实现：一种拉曼成像的光谱数据中宇宙射线Spike的识别方法，包括以下步骤：

1.拉曼成像中生成了m个原始拉曼光谱s，每个光谱有n个波数，对其中一个原始拉曼光谱s执行半窗口值为w的开运算，得到开运算光谱so，半窗口值w取值应不小于拉曼光谱仪产生的最大的Spike峰宽的一半；

2.用原始拉曼光谱s每个波数的强度值减去开运算光谱so对应波数的强度值，得到比较光谱sc，即sc＝s-so，再对比较光谱sc执行半窗口值为6的膨胀运算，得光谱sc'，再用光谱sc'中每个波数的强度值减去该光谱sc'所有强度值的最小值，得到目标光谱ss；

3.重复步骤1～2，得到与m个原始拉曼光谱s对应的m个开运算光谱so，m个比较光谱sc和m个目标光谱ss；

4.求第i个和第i-1个目标光谱ss的夹角的余弦值cos＜ss_i-1,ss_i＞，计算公式为为：

式中<ss_i-1,ss_i＞表示ss_i-1和ss_i的内积，||ss_i-1||表示ss_i-1的二范数；

5.若cos＜ss_i-1,ss_i＞小于第一阈值th₁，则第i个拉曼光谱s_i疑似含有Spike；

6.对步骤(5)中第i个拉曼光谱s_i进一步确认Spike的存在，将比较光谱sc_i第j个波数处的强度值记为sc_i(j)，遍历全部n个波数，找到使得第i个比较光谱sc_i和第i-1个比较光谱sc_i-1强度值之差取得最大值的第k个波数，即：

7.若sc_i(k)的强度值大于第二阈值th₂，则判定原始光谱s_i在第k个波数附近出现了Spike，其中th₂设定为sc_i的平均值加上4倍的sc_i标准差，公式为：

本发明中第一阈值th₁的取值范围为0.8～0.3，th₁取值偏大有助于识别强度微弱的Spike但会增加误判的可能，而th₁取值偏小可能遗漏Spike的识别，若拉曼光谱之间的连续性好，即连续得到的两个拉曼光谱形态非常相似，可采用偏大的th₁，反之应采用偏小的th₁。

本发明第二目的可通过以下技术措施来实现：一种拉曼成像数据中宇宙射线Spike的修正方法，按上述识别方法，取i值为2，即从第2个原始光谱S₂开始判断是否有Spike，如没有，则继续按序识别下一个原始光谱S₃；如有，则去除当前宇宙射线Spike进行修正，并用修正后的原始光谱S₂重新计算比较光谱sc₂和目标光谱ss₂，用作第3个原始光谱S₃的识别；依次重复，除第一个原始光谱外，其它m-1个原始光谱中集中判定出所含有Spike的光谱。而第一个原始光谱，不对其进行Spike的识别，或设定不存在Spike。

上述去除当前宇宙射线Spike的具体过程：将上述判定出现了Spike的原始光谱s_i第k个波数前后w处的强度值：s_i(j+r),r＝-w,...,w；一共2w+1个强度值替换为开运算光谱so_i对应的强度值；若k处于光谱的开头或末尾，不足2w+1个强度值，忽视那些超过开头或末尾的部分，完成强度值的替换后获得去除当前Spike的修正后的拉曼光谱。

本发明具有以下技术效果：

本发明技术方案不再直接用两个光谱的线性近似谱进行比较，而是利用可能包含Spike和小量尖峰信息的目标光谱ss进行比较，减少连续光谱相似度低时为判别Spike带来的干扰，提高了识别的正确率，能有效在拉曼成像数据中集中识别出含有的全部Spike，并对拉曼光谱进行修正。

附图说明

图1异常光谱的位置信息，图(a)～(d)分别对应第2748～2751条光谱的cos＜ss_i-1,ss_i＞值，图中的虚线为y＝0.75，图(a)～(d)中标注的四个点均在虚线以下，有存在Spike的可能。

图2比较光谱sc₂₇₄₇和sc₂₇₄₈，其强度值在1024cm^-1差异明显。

图3(a)第2748～2751条原始拉曼光谱，这四条连续的光谱均在波数1042cm-1的位置出现Spike；(b)修复后的第2748～2751条拉曼光谱。

图4连续得到两条信号较弱的拉曼光谱s₁和s₂，光谱s₂在波数1238cm^-1的位置出现低强度的Spike；

图5目标光谱ss₂和ss₁，目标光谱ss₂保留了Spike的特征。

图6连续得到两条拉曼光谱s_a和s_b，在拉曼峰的位置，强度出现巨大差异。

图7目标光谱ss_a和ss_b，拉曼峰位置的剧烈变化在目标光谱ss中得到减弱，因此不会出现误判Spike的情况。

图8光谱s_b和其拟合光谱s'_b在拉曼峰的位置仍有较大拟合误差。

图9差谱s_b-s'_b在拉曼峰的1320cm^-1位置出现高强度值318，显示出光谱s_b和其拟合光谱s'_b在1320cm^-1的位置仍有较大拟合误差容易引起Spike的误判。

具体实施方式

本发明实施例中采用的常用的运算方式如下：

腐蚀运算(Erosion)：

E(s(j))＝min(s(j+r)),r＝-w,...,w (1)

膨胀运算(Dilation)：

D(s(j))＝max(s(j+r)),r＝-w,...,w (2)

开运算(Open):

O(s(j))＝D(E(s(j+r))),r＝-w,...,w (3)

式中s(j)是光谱s在第j个波数处的强度值，w是半窗口值，腐蚀运算是把光谱某点的强度值替换为窗口内共2*w+1个强度值的最小值，膨胀运算则相反，若第j个波数处于光谱开头或末尾，使窗口内不能凑齐2*w+1个强度值，则忽视那些超过开头或末尾的部分，开运算是对光谱先执行腐蚀运算再执行膨胀运算。

实例一

实验仪器为雷尼绍公司的共聚焦拉曼显微镜(Renishaw in Via Confocal RamanMicroscope)，对纳米金孵育后的HepG2细胞进行拉曼成像，共计得到3640条拉曼光谱。

对这些光谱中的Spike进行识别和去除，本例中半窗口值w取值4，第一阈值th₁取值0.75。

1.从第1个光谱s₁开始，执行半窗口值w为4的开运算，求得开运算光谱so₁。

2.用s₁减去so₁得到比较光谱sc₁，对sc₁执行半窗口值为6的膨胀运算，得光谱sc'₁。

3.光谱sc'₁的每一个波数的强度值减去光谱sc'₁强度的最小值min(sc'₁)，得到目标光谱ss₁。

4.按上述步骤求得1～3640个原始拉曼光谱对应的开运算光谱so，比较光谱sc和目标光谱ss。

5.令i遍历2到3690，计算第i-1和i个目标光谱ss的夹角余弦值cos＜ss_i-1,ss_i＞，若cos＜ss_i-1,ss_i＞不小于第一阈值th₁，则继续计算第i和i+1个目标光谱ss的夹角余弦值cos＜ss_i,ss_i+1＞，若cos＜ss_i-1,ss_i＞小于第一阈值th₁，则如步骤6～10进一步分析；共计算得3639个cos＜ss_i-1,ss_i＞值，以cos＜ss_i-1,ss_i＞为纵轴，以光谱序号为横轴，作图1。

6.以第2748～2751连续的4条拉曼光谱为例，若cos＜ss_i-1,ss_i＞小于th₁后的步骤；第2748的cos＜ss_i-1,ss_i＞值为0.5657，如图1(a)，小于th₁＝0.75，初步怀疑该光谱含有Spike。

7.进一步判断第2748条光谱是否含有Spike，找到第2748条光谱的比较光谱sc₂₇₄₈，及前一条不含spike的光谱或已经修正过的第2747条光谱的比较光谱sc₂₇₄₇，如图2。

8.满足的k值为365，而第k＝365个波数为1042cm^-1，即两个比较光谱sc₂₇₄₇和sc₂₇₄₈强度相差最大的位置是1042cm^-1，而sc₂₇₄₈在1042cm^-1的强度值为1951，sc₂₇₄₈的均值和标准差为别为24.3376和72.1299，算得th₂＝24.3376+4*72.1299＝312.8573，sc₂₇₄₈在1042cm^-1的强度值远大于th₂，到此确定第2748条光谱在1024cm^-1的位置出现Spike。

9.将第2748条光谱s₂₇₄₈第(k-w)～(k+w)个强度值，即第361～369个，共9个强度值换成光谱so₂₇₄₈对应波数的强度值，执行后原始光谱剔除了Spike，并以修正的光谱重新计算其比较光谱sc₂₇₄₈目标光谱ss₂₇₄₈，

10.第2749～2751条光谱，在前一条完成修正光谱的基础上，以同样的方法可以算得2749～2751条光谱的cos＜ss_i-1,ss_i＞值分别为0.3136，0.3211，0.6491，如图1(b)～(d)，均小于0.75，初步怀疑这些光谱含有Spike。

11.依次重复7～9的步骤，完成Spike的识别及修正。第2748～2751条原始光谱如图3(a)；执行后剔除了Spike的第2748～2751条光谱如图3(b)。

12.本例对第一条光谱外3639条光谱进行Spike的识别，正确识别及修正9处Spike。第一条原始光谱，不对其进行Spike的识别，或设定不存在Spike。

实例二

实例二显示了本发明在低信噪比下识别Spike的能力。如图4，展示了两条在拉曼成像中连续得到的拉曼光谱和s₂，其中s₂在波数1238cm^-1的位置出现Spike，

s1

但强度非常低，为161.805。本例中将半窗口值w取值为4，第一阈值th₁取值为0.75。

1.对拉曼光谱s₁和s₂执行半窗口值w为4的开运算，求得开运算光谱so₁和so₂。

2.用s_i减去so_i得到比较光谱sc_i，对sc_i执行半窗口值为6的膨胀运算，得光谱sc'_i，光谱sc'_i的每一个波数的强度值减去光谱sc'_i强度的最小值min(sc'_i)，得到目标光谱ss_i，如图5，其中i＝1，2。

3.计算cos＜ss₁,ss₂＞值为0.6849，小于第一阈值th₁＝0.75，进一步进行判断。

4.满足的k值为543，而第k＝543个波数为1238cm^-1，光谱sc₂在1238cm^-1强度值为145.6141，算得th₂＝45.286，sc₂在1238cm^-1强度值大于th₂，正确判断第2条光谱在1238cm^-1处出现Spike。

实例三

实例三显示本发明在信号较强，光谱变化剧烈的场合不会错误识别Spike的能力。两条在拉曼成像中连续得到的拉曼光谱s_a、s_b，由于成像过程中空间位置的变化，引起扫描样品点成分的变化，这种变化可能造成在拉曼峰的位置，强度出现巨大差异，如图6，两条连续光谱在波数1320cm^-1的位置，强度差高达500以上，其值几乎是例二中Spike强度的三倍！若使用一般的近邻光谱的比较方法，很可能出现误判为Spike的情况。本例中同样将半窗口值w取值为4，第一阈值th₁取值为0.75。

1.对拉曼光谱s_a和s_b执行半窗口值w为4的开运算，求得开运算光谱so_a和so_b。

2.用s_i减去so_i得到比较光谱sc_i，对sc_i执行半窗口值为6的膨胀运算，得光谱sc'_i，光谱sc'_i的每一个波数的强度值减去光谱sc'_i强度的最小值min(sc'_i)，得到目标光谱ss_i，如图7，其中i＝a，b。

3.计算cos＜ss_a,ss_b＞值为0.9362，大于第一阈值th₁＝0.75，s_b判定为正常光谱，没有出现误判。

本发明能有效在拉曼成像数据集中识别出含有Spike的拉曼光谱，比起现有方案，识别的正确率更高。对比分析如下：

1.相比单单依靠峰值强度或峰宽开辨识Spike的方法更加灵活可靠。如专利CN103674251B中提及，判断为Spike之依据：峰值与峰宽之比大于1000；而本发明实例二中出现的Spike峰值强度不足200，显然由上述方法无法正确判断Spike的存在。

2.相比现有的最近邻比较法，能在近邻光谱相差较大时减小误判的可能。如本发明实例三，按照传统最近邻比较法，使用最近邻光谱的一个线性组合作为近似，如式(4)。其中p、q是拟合系数，均为标量，

s'_b＝p*s_a+q>

但其拟合效果有时还是不能令人满意，如图8。进一步作出光谱s_b和其拟合光谱s'_b之差s_b-s'_b，如图9。不少文献将判断的阈值设置为差谱s_b-s'_b的标准差的3.5倍，若差谱的强度值超过此阈值，则判断该处出现Spike。而图9中强度值的标准差值为47.798，最大值为318，最大强度是标准差的6倍以上，显然会造成Spike的误判。

本发明的上述实施例并不是对本发明保护范围的限定，本发明的实施方式不限于此，凡此种种根据本发明的上述内容，按照本领域的普通技术知识和惯用手段，在不脱离本发明上述基本技术思想前提下，对本发明上述结构做出的其它多种形式的修改、替换或变更，均应落在本发明的保护范围之内。