CN114720433A - 荧光光谱无用散射自动去除方法 - Google Patents

Abstract

本发明属于分析化学领域，具体涉及一种自动识别散射区域，去除散射数据并填补去除部分的计算方法。包括以下步骤：首先，对于二维的光谱数据，通过拉曼散射和瑞利散射机理以及其对应的能量判断散射位置；然后，在一维光谱数据的基础上，通过散射位置和散射强度判断散射能量的覆盖范围，并将无用能量进行剔除；最后，对于去除散射数据的区域进行插值处理，达到去无用除散射并补齐的目的。本发明可以快速的自动识别出无用散射的位置，不需要人为设定散射能量的宽度，在可接受的误差范围内，最大限度的保留数据，这大大的提高了去除散射过程中的准确性和实现效率。这为实现大量含散射的光谱数据的批量处理提供可能。

Description

荧光光谱无用散射自动去除方法

技术领域

本发明利用数值计算方法去除三维荧光光谱中的无用散射数据，属于分析化学领域，具体涉及一种自动识别散射区域，去除散射数据并填补去除部分的计算方法。

背景技术

荧光激发发射矩阵测量在食品科学、分析化学、生物化学和环境科学等领域都有广泛的应用。但由于荧光光谱的生成机制所致，存在瑞利散射和拉曼散射，对其进行分析和去除往往是复杂的。现阶段存在若干种方法来处理无用散射的影响。然而，现有的方法大部分都存在相应的不足和待改进的地方，往往存在只能单次处理样本数据，同时在处理数据的过程中不能够自动识别散射位置只能主观判断散射的位置的问题。

本发明是基于现有的去散射的方法的基础上做出的改良，原有的方法依赖人为设定参数，需要光谱的使用人员主观判断无用散射能量所覆盖的区域，而本发明是基于光谱数据本身，通过对图像数据进行处理的方式自动判断散射的位置和范围，达到去散射并补齐数据的目的。

发明内容

本发明目的是提供一种本发明提出了高效自动去除散射的计算方法。

本发明为实现上述目的所采用的技术方案是：

荧光光谱无用散射自动去除方法，包括以下步骤：

1)通过荧光光谱仪获取有机物中的荧光光谱数据，并对荧光光谱数据进行预处理；

2)取出预处理后的荧光光谱数据中荧光强度大于阈值的若干个数据，并对其进行拟合，得到拟合参数；

3)根据拟合参数，构建荧光光谱数据的表达式：y＝kx+b，其中，y为荧光光谱的激发波长，x为荧光光谱的发射波长，k、b为拟合参数；

4)取x为某一值时，y所对应的荧光光谱范围中，荧光强度的最大值；

5)对步骤4)中的荧光强度最大值进行处理，去除该发射波长下的散射；

6)重复步骤4)～步骤5)，直至对全部荧光光谱完成处理，对每一个发射波长下处理后的荧光光谱进行整合，得到完整的荧光光谱；

7)判断处理后的荧光光谱是否还存在散射，若不存在，则输出荧光光谱至上位机，使其可视化；若存在，则重复步骤2)～步骤6)继续对其处理，并将处理后的荧光光谱输出至上位机，使其可视化。

所述散射包括：瑞利散射以及拉曼散射，其中，瑞利散射与拉曼散射的拟合参数不同，所述步骤7)中若还存在散射，则存在的散射为拉曼散射。

采用最小二乘法对预处理后的荧光光谱数据进行拟合。

所述对步骤5)中的荧光强度最大值进行处理，具体为：

对荧光强度曲线进行平滑处理，并找到该曲线中与波峰左右相邻的两个波谷；

去除两个波谷之间的数值；

对去除数据的曲线进行插值处理，得到去除散射的荧光强度曲线。

使用均值滤波器对荧光强度曲线进行平滑处理。

本发明具有以下有益效果及优点：

1.本发明可以准确的自动的识别到散射位置和范围，避免了人为输入散射范围主观性，增加了算法的普棒性，所得到的数据结果更加的快速，耗时更短。

2.本发明采用了更加精确的方法识散射区间，比较原有算法人为输入散射区间，更加的精准，大大减小了误差，更符合实际。

3.本发明用插值法拟合数据，体现非散射部分光谱数据的性质。

附图说明

图1是三维光谱自动去散射流程图；

图2是待处理拉曼散射和瑞利散射的三维光谱图；

图3是保留能量最高点流程图；

图4a是激发波长与荧光强度原始数据示意图；

图4b是平滑处理原始数据示意图；

图4c是处理后的数据一阶求导示意图；

图4d是处理后的数据二阶求导示意图；

图5是未去散射前的数据图一；

图6是去散射能量值的数据图一；

图7是去散射后插入数据的数据图一；

图8是未去散射前的数据图二；

图9是去散射能量值的数据图二；

图10是去散射后插入数据的数据图二；

图11是去除拉曼散射和瑞利散射的三维光谱图；

图12a是待处理拉曼散射和瑞利散射的三维光谱数据示意图一；

图12b是去除拉曼散射和瑞利散射的三维光谱数据示意图二；

图13a是待处理拉曼散射和瑞利散射的三维光谱数据示意图一；

图13b是去除拉曼散射和瑞利散射的三维光谱数据示意图二；

图14a是待处理拉曼散射和瑞利散射的三维光谱数据示意图一；

图14b是去除拉曼散射和瑞利散射的三维光谱数据示意图二；

图15a是待处理拉曼散射和瑞利散射的三维光谱数据示意图一；

图15b是去除拉曼散射和瑞利散射的三维光谱数据示意图二；

图16a是待处理拉曼散射和瑞利散射的三维光谱数据示意图一；

图16b是去除拉曼散射和瑞利散射的三维光谱数据示意图二；

图17a是待处理拉曼散射和瑞利散射的三维光谱数据示意图一；

图17b是去除拉曼散射和瑞利散射的三维光谱数据示意图二；

图18a是待处理拉曼散射和瑞利散射的三维光谱数据示意图一；图18b是去除拉曼散射和瑞利散射的三维光谱数据示意图二。

具体实施方式

下面结合附图及实施例对本发明做进一步的详细说明。

三维荧光光谱是近些年广泛应用而发展起来的一种荧光分析技术，它在测量在食品科学、分析化学、生物化学和环境科学等领域都有广泛的应用。三维荧光光谱是由激发波长(Emission，Em),发射波长(Excitation，Ex),荧光强度(荧光强度数值用Z表示)的三维坐标所表征的矩阵光谱，是描述荧光强度及同时随激发波长和发射波长变化的关系图谱。三维荧光光谱不仅能够获得激发波长与发射波长，同时能够获取变化时的荧光强度信息。

目前存在若干种处理散射效应的方法，但受荧光光谱的生成机制所影响，在荧光光谱中往往存在瑞利散射和拉曼散射对应的能量，这些能量的幅值相对较大，往往会对目标物的分析产生巨大的影响。所以，在分析目标物能量之前，有必要将其去除。本发明是基于原有的三维荧光光谱去散射算法进行的改进，三维荧光光谱是由若干张二维荧光光谱组成，本发明中所有的说明均是基于二维荧光光谱进行的。如图1所示，包括以下步骤：

首先，对于二维的光谱数据，通过拉曼散射和瑞利散射机理以及其对应的能量判断散射位置；然后，在一维光谱数据的基础上，通过散射位置和散射强度判断散射能量的覆盖范围，并将无用能量进行剔除；最后，对于去除散射数据的区域进行插值处理，达到去无用除散射并补齐的目的。

迭代选择高能量点，拟合无用散射位置表达式，根据无用散射位置直线表达式，确定无用散射最高能量点的位置，包括以下步骤：

将受散射影响的三维光谱荧光数据导入，一般拉曼散射或瑞利散射的能量值，要高于目标物的散射能量值，且在荧光光谱上呈现为近似直线区域。取拉曼散射或瑞利散射部分最高能量点，通过迭代，不断判断拟合出的散射位置是否合理，在处理二维荧光数据时，通过直线表达式判断无用散射能量最高点的位置。

平滑数据，根据无用散射最高能量点，确定散射两侧波谷，去除无用散射区域，包括以下步骤：

采用移动均值滤波器对荧光数据列向量进行平滑处理，通过之前步骤中确定的最高能量点，搜索到其对应的两侧的波谷，确定无用散射的区间位置，去除无用散射区间内的散射数据。

本发明的内容在于提供了一种基于三维荧光光谱去无用散射的计算方法。方法的具体步骤：

(1)拟合出瑞利散射的直线表达式

为了分析本算法的能力，模拟生成了一组三维荧光光谱数据，该数据集存在拉曼散射，瑞利散射和目标物散射的荧光数据，以等高线的形式表示，具体的如图2所示。以此数据集为例进行说明。

第一步 荧光光谱数据的导入说明。

通过荧光光谱仪，选择合适的激发波长和发射波长，调整参数扫描样本，记录荧光强度的变化，得到三维光谱图，每个样本的三维荧光光谱生成一个代表了样本中有机物荧光特征图谱的荧光数据。

荧光数据是由激发波长em(对应图中的y轴),发射波长ex(对应图中的x轴),荧光强度(对应图中的z轴)三维坐标所表征的矩阵光谱。实验数据集用纳米间隔测量，激发波长在200-250纳米之间，发射波长范围为200-250纳米。

数据的存储形式为dat.xls或者为dat.txt，将待去散射的三维荧光光谱的数据以矩阵的形式导入。该矩阵的横坐标为矩阵的列序号(从1开始计数)，该矩阵的纵坐标为矩阵的行序号，具体的行列对应的数值为该点的荧光强度数据(z_(l,r)为具体荧光强度值，l与r为所在位置的坐标)。例如，数据矩阵为

以z_(50,1)为例，其横坐标为50，纵坐标为1，z轴坐标为z_(50,1)。表示为发射波长为50nm，激发波长为1nm时，荧光强度为z_(50,1)。

如图2所示，通过等高线的形式将三维荧光光谱强度数据进行可视化处理。其中横轴是发射波长ex；纵轴是激发波长em；等高线的值表示荧光强度的大小，从0.2开始，间隔为0.2，至1.4结束，等高线的颜色越深代表荧光强度越大。通过数据集生成的荧光光谱的等高线图，可以观察到在中间圆形等高线的两侧有两条接近直线的等高线数据，中间圆形等高线为没被散射影响部分的荧光数据，两条接近直线的等高线为拉曼散射和瑞利散射数据。

第二步 保留数据中能量值高的点。

如图3所示，受散射影响的三维光谱图中，一般，无用散射的能量值高于目标散射的荧光能量值，在三维荧光光谱上呈现为近似直线区域。取散射部分最高能量点，拟合散射位置和能量值的直线表达式。瑞利散射直线表达式的斜率和截距在一定范围内的，瑞利散射的位置定义为激发波长等于发射波长，即斜率(k表示斜率)范围在1左右，截距(b表示截距)小于0。为了方便判段斜率的取值范围，且斜率在1左右，误差不超过0.5，所以斜率在0.5-1.5之间有概率存在瑞利散射。即瑞利散射直线表达式截距和斜率的范围为b＜0,0.5<k<1.5。

通过拟合直线表达式的斜率和截距进行判断。在范围内的表达式为散射位置和能量值的直线表达式(其中y对应激发波长em的值，x对应发射波长ex的值)。具体操作流程如下：

首先，保留部分高能量的点。高能量点的数据保留，高能量以外的数据点设置为0，如

将只保留高能量点的荧光矩阵，通过最小二法拟合直线y＝kx+b，

其中n为坐标数。判断拟合出直线的斜率k以及截距b是否b＜0,0.5<k<1.5范围内。

若在b＜0,0.5<k<1.5范围内，直接输出拟合出的直线斜率与截距；

若不在b＜0,0.5<k<1.5范围内，即b≥0或k≥1.5或者k≤0.5，增加高能量点的数量，增加高能量点的过程以每次迭代拟合增加光谱中点的数量nat％来实现。Nat定义为1～5，如

若仍然不在斜率与截距的b＜0,0.5<k<1.5范围内，继续递加0.001％高能量点，如，

直至斜率和截距保持b＜0,0.5<k<1.5范围内，输出高能量点与拟合出的直线斜率k与截距b。

第三步 拟合出散射的直线表达式，找到散射的最高能量位置。

根据拟合出的无用散射对应范围的直线的斜率k和截距b，组成散射的直线表达式y＝kx+b。其中，y和x为直线表达式的两个变量，是散射能量最高点所在的坐标的位置，分别为输入数据矩阵的横坐标ex和纵坐标em。

为了进一步去除散射，需提取每一组受散射影响的数据。将发射波长ex(对应图中x轴)设为固定值，激发波长em(对应图中y轴)和荧光强度数值(对应图中z轴)组成一组数据。再根据散射的直线表达式y＝kx+b，确定当发射波长ex为固定值时，激发波长em的值，得到的激发波长em为该组数据内散射最高能量的位置。如图2所示，提取x轴坐标固定为20时，y轴坐标和z轴坐标相对应的数；接下来，根据瑞利散射的直线表达式y＝kx+b，确定x＝20时最高能量点激发波长最高能量点的位置y，图2中当x＝20时，y大约为20左右。

(2)根据散射位置去除瑞利散射

被散射影响的荧光数据集一般能够出现两种类型。如图5，图8所示，通过二维曲线的形式将激发波长em数据进行可视化处理。其中横轴是激发波长em；纵轴是确定发射波长ex下激发波长的强度(例如图2中，x轴等于20时，y轴与z轴的变化曲线)；其中虚线为实际的荧光数据集，粗线为散射的影响数据集，细线为理想的没有散射影响的荧光数据集。本发明的目的就是根据虚线的实际荧光数据集，去除粗线的散射影响，最终生成近似于细线的去除散射后的数据集。

一种情况为图5所示，散射的影响数据集(粗线)和无散射的荧光数据集(细线)距离较近，两者之间的数据有重叠部分。另一种情况为图8所示，散射影响的数据集(粗线)和无散射的荧光数据集(细线)距离比较远，两者之间的数据没有重叠部分。将两种数据类型进行处理。

第四步 平滑处理曲线。

图4中，图4a表示通过二维曲线的形式将激发波长em数据进行可视化处理，其中横轴是激发波长em，纵轴是确定发射波长ex下，激发波长的强度的原始数据；图4b则为在图4a的基础上对二维曲线进行平滑处理后的数据曲线表示；图4c是在图4b的数据下对数据进行一阶求导，三角形标识的位置为导数为0的点；图4d是在图4b的数据下对数据进行二阶求导，圆形标识的位置为二阶导数为0的点。

由于实际的光谱中会存在一定的波动，如图4a所示，影响计算结果，在不影响原始数据的基础上对每一条被散射影响的曲线进行平滑处理，如图4b所示。可以采用移动均值滤波器对荧光数据列向量z_n进行平滑处理，窗宽为kk＝2m+1，m为节点数。返回去噪后数据列向量Z_n。当原数据点z_n(w)前后累加元素总和可以等于窗宽kk时，

当原数据点z_n(w)前后累加元素不等于窗宽kk时，取除z_n(w)外最近对称元素，将包括z_n(w)在内的所取元素累加，再除以元素总个数。例，移动均值滤波的窗宽为5，Z_n中的元素计算方法如下：

Z_n(1)＝z_n(1)

Z_n(2)＝(z_n(1)+z_n(2)+z_n(3))/3

Z_n(3)＝(z_n(1)+z_n(2)+z_n(3)+z_n(4)+z_n(5))/5

Z_n(4)＝(z_n(2)+z_n(3)+z_n(4)+z_n(5)+z_n(6))/5

Z_n(5)＝(z_n(3)+z_n(4)+z_n(5)+z_n(6)+z_n(7))/5

…

将曲线z_n中的每一数据点进行降噪处理保存为Z_n＝[Z_n(1),Z_n(2),Z_n(3)…Z(w)]。

第五步 找到散射两侧的波谷。

根据每一条曲线散射最高能量的位置。确定散射的区间位置，发现最高能量点两侧波谷，通过实际操作发现图5和图8的两种情况都可以快速的找到波谷的位置。具体操作如下：

a)计算出曲线中所有的波谷

Z_n(w)是平滑降噪处理后数据的集合，Z_n(w)表现为一条平滑的曲线，则可以对Z_n(w)进行求导得到Z′_n(w)(一阶导数)，并求出Z_n″(w)＝0(二阶导数等于零)时的根，判断Z′_n(w)左右两侧一阶导数的正负。如果Z′_n(w-1)＜0Z′_n(w+1)＞0那么Z_n(w)在这个根处取得极小值，即为波谷。

如图4所示，图4a为原数据集；图4b为平滑后的数据集，其点标识为所有波谷所在位置，两侧箭头的位置为波峰两侧的波谷位置；图4-3为该数据集一阶求导，其点标识为二阶导数为0且Z′_n(w-1)＜0Z′_n(w+1)＞0的点；图4d为二阶求导，其点标识为二阶导数为0的点。在二阶求导图(图4d)中确定Z″_n(w)＝0的根，如图4d中的点标识；然后判断一阶求导图(图4c)中Z′_n(w-1)＜0Z′_n(w+1)＞0时m的值，如图4c中的点标识；最后确定数据集中波谷的位置，图4b中的方框标识为波谷位置。

b)计算出前后两个波谷

在确定所有的波谷后，通过最高能量点的位置，确定对应距离该点最近的左右两个波谷。即根据

其中(x₁,y₁)为波谷的位置坐标，波峰的位置是(x₂,y₂)，确定波峰左右两侧波谷的R最小值。如图4b中的箭头标识为左右两波谷位置。

这两个波谷内的区域为需要去除的散射区域。例如，若数据集为[Z_n(1),Z_n(2),Z_n(3),Z_n(4),Z_n(5),Z_n(6),Z_n(7),Z_n(8),Z_n(9),Z_n(10)]，两波谷分别为Z_n(4)和Z_n(9)，则Z_n(4)和Z_n(9)以及之间的范围为需去除的散射区域。

第六步 去除散射的影响

确定前后波谷位置后，两侧波谷范围内的数据设置为空集值。例如，若数据集为[Z_n(1),Z_n(2),Z_n(3),Z_n(4),Z_n(5),Z_n(6),Z_n(7),Z_n(8),Z_n(9),Z_n(10)]，两波谷分别为Z_n(4)和Z_n(9)，则该条曲线去散射影响后的数据集为[Z_n(1),Z_n(2),Z_n(3),NAN,NAN,NAN,NAN,NAN,NAN,Z_n(10)]，NAN表示为非数据的类型(not a number)。

如图6，图9所示为在图5，图8数据基础上去除散射部分后的图像，其中横轴是激发波长em；纵轴是确定发射波长ex下，激发波长的强度；虚线为设置空集的数据集，粗线为散射的影响数据集，细线为理想的没有散射影响的荧光数据集，这时散射的数据已经被去除。

第七步 对去除后的曲线进行插值。

去除无用散射的数据后，根据剩余的数据判断散射数据的变化趋势。得到一组完整的数据。如图7，图10所示，为在图6，图9数据基础上对去除无用数据后插值的图像，其中横轴是激发波长em；纵轴是确定发射波长ex下，激发波长的强度；虚线为插值处理后的数据集，粗线为散射的影响数据集，细线为理想的没有散射影响的荧光数据集。

具体公式如下：

假定有n+1个数据节点(y₀,z₀),(y₁,z₁),(y₂,z₂)......(y_n,z_n)，y与z对应激发波长em与荧光强度z的值，n为发射波长的点的个数，计算步长h_i＝y_i+1-y_i(i＝0,1,…,n-1)，将数据节点和指定的首位端点条件带入矩阵方程。解矩阵方程，求得二次微分值m_i。该矩阵为三对角矩阵。计算样条曲线的系数：

a_i＝z_i

其中i＝0,1,…,n-1。在每个子区间x_i≤x≤x_i+1中，创建方程

g_i(y)＝a_i+b_i(y-y_i)+c_i(y-y_i)²+d_i(y-y_i)³，a_i b_i c_i d_i为插值后曲线对应的系数，g_i(y)为插值之后的曲线方程。

在具体计算实践中发现，无论散射峰值和荧光数据距离较近还是较远，处理散射的方法是一致的，所以两种情况均可

以采用本算法。

第八步 整合得到去除散射的三维光谱。

对每组数据进行单独处理，得到g_i(y)(i＝0,1,…,n-1)，再将处理完的数据重整合成一个三维光谱的矩阵

画出图像，此时瑞利散射已经去除完全了。

(3)去除拉曼散射

部分的荧光光谱只存在一种散射，而另一部分荧光光谱具备拉曼和瑞利两种散射，因此，需对光谱进行判断。判断其是否存在两种散射，若存在两种散射，需对两种散射分别进行处理；而若只存在一种散射，直接输出上一步的三维光谱的矩阵

即可。

第九步 判断是否有拉曼散射

判断是否存在拉曼散射，通过拟合最高能量点，判断散射的直线表达式的斜率(v表示斜率)以及截距(u表示截距)的范围，拉曼散射的位置位于激发波长是发射波长的2倍处，即斜率2处，为了方便判段斜率的取值范围，且斜率在2左右，误差不超过0.5，所以斜率在1.5-3之间有概率存在拉曼散射。即拉曼散射直线表达式斜率的范围为，1.5<v<3。若不在范围内，则说明该光谱只存在瑞利散射，输出三维荧光光谱的矩阵

若在范围内，则说明存在拉曼散射，运用上一步算法处理出三维荧光光谱的矩阵输出

(i＝0,1,…,n-1)。

(4)输出去散射之后的图像

处理散射之后，输出数据矩阵G，该矩阵的横坐标为矩阵的列数代表发射波长ex，纵坐标为矩阵的行数代表激发波长em，具体的行列对应的数值为该点的荧光强度数据

通过画图，分析去除散射后三维荧光光谱图，如图11所示，为图2去除无用散射后的三维荧光光谱图。其中横轴是发射波长ex；纵轴是激发波长em；等高线的值表示荧光强度的大小，从0.2开始，间隔为0.2，至1.4结束，等高线的颜色越深代表荧光强度越大，图中的拉曼散射和瑞利散射已全部去除。

图12-图18为通过等高线的形式将三维荧光光谱强度数据进行可视化处理。其中横轴是发射波长ex；纵轴是激发波长em；等高线的值表示荧光强度的大小，等高线的颜色越深代表荧光强度越大。左侧a图是原始数据图像，右侧b图为去除散射后的三维荧光光谱图像。此时观察去除无用散射后的三维荧光光谱图像，近似直线的拉曼散射和瑞利散射已去除。

本发明可以快速的自动识别出无用散射的位置，不需要人为设定散射能量的宽度，在可接受的误差范围内，最大限度的保留数据，这大大的提高了去除散射过程中的准确性和实现效率。这为实现大量含散射的光谱数据的批量处理提供可能。

Claims (5)

1.荧光光谱无用散射自动去除方法，其特征在于，包括以下步骤：

1)通过荧光光谱仪获取有机物中的荧光光谱数据，并对荧光光谱数据进行预处理；

2)取出预处理后的荧光光谱数据中荧光强度大于阈值的若干个数据，并对其进行拟合，得到拟合参数；

3)根据拟合参数，构建荧光光谱数据的表达式：y＝kx+b，其中，y为荧光光谱的激发波长，x为荧光光谱的发射波长，k、b为拟合参数；

4)取x为某一值时，y所对应的荧光光谱范围中，荧光强度的最大值；

5)对步骤4)中的荧光强度最大值进行处理，去除该发射波长下的散射；

6)重复步骤4)～步骤5)，直至对全部荧光光谱完成处理，对每一个发射波长下处理后的荧光光谱进行整合，得到完整的荧光光谱；

7)判断处理后的荧光光谱是否还存在散射，若不存在，则输出荧光光谱至上位机，使其可视化；若存在，则重复步骤2)～步骤6)继续对其处理，并将处理后的荧光光谱输出至上位机，使其可视化。

2.根据权利要求1所述的荧光光谱无用散射自动去除方法，其特征在于，所述散射包括：瑞利散射以及拉曼散射，其中，瑞利散射与拉曼散射的拟合参数不同，所述步骤7)中若还存在散射，则存在的散射为拉曼散射。

3.根据权利要求1所述的荧光光谱无用散射自动去除方法，其特征在于，采用最小二乘法对预处理后的荧光光谱数据进行拟合。

4.根据权利要求1所述的荧光光谱无用散射自动去除方法，其特征在于，所述对步骤5)中的荧光强度最大值进行处理，具体为：

对荧光强度曲线进行平滑处理，并找到该曲线中与波峰左右相邻的两个波谷；

去除两个波谷之间的数值；

对去除数据的曲线进行插值处理，得到去除散射的荧光强度曲线。

5.根据权利要求4所述的荧光光谱无用散射自动去除方法，其特征在于，使用均值滤波器对荧光强度曲线进行平滑处理。

Images