基于特征区域下凸点提取的藻类荧光光谱波长选择方法

张永彬，朱丹丹，陈 颖*，刘 喆，段玮靓，李少华

1.燕山大学电气工程学院河北省测试计量技术及仪器重点实验室，河北 秦皇岛 066004 2.河北先河环保科技股份有限公司，河北 石家庄 050000

引 言

近年来我国工农业发展迅速，海水富营养化加剧，水体中浮游藻类大量繁殖，频繁引发藻华现象，褐潮成因藻抑食金球藻、常见藻华中的小球藻和细长聚球藻已经成为海洋生态环境的主要污染，影响人类健康安全和环境发展，所以藻类种类和数量的测量对生态环境检测具有重要意义[1-2]。浮游植物中含有的色素具有荧光效应，因此荧光光谱分析是一种有效的测量方法，其中三维荧光光谱法在藻类测量中应用广泛，如何选择光谱波长变量从而提高三维荧光光谱的信噪比是目前一个重要研究方向。

三维荧光光谱是在已设置的激发波长和发射波长范围下扫描所有波长点而得到的数据矩阵，其数据较多，但并不是光谱区域下的所有波长变量都含有有用信息，只有少数光谱区域含有物质特征。若冗余区域波长变量参与建模会影响预测精度，增加模型的计算，光谱预处理虽然可以提高信噪比，但并没有去除不相关的波长变量，不能最大程度提高模型的精度，没有减小模型复杂度。

传统的波长选择方法是峰值法，对固定波长点下的荧光峰进行分析，但该方法所选择的波长变量过少，原光谱的数据信息没有得到充分利用，无法有效反映整体区域荧光强度变化情况[3]。目前的波长选择方法有无信息变量消除法(uninformative variables elimination，UVE)[4]、连续投影算法(successive projections algorithm，SPA)[5-6]、遗传算法(genetic algorithm，GA)[7-8]等。SPA通过投影方式选取线性关系最小的波长组合，减少光谱信息中的波长变量，但SPA所选择的波长变量的个数不能超过所用样品数量，其受样品数量的限制较大。UVE将随机噪声矩阵加入光谱矩阵并建立偏最小二乘(partial least-square，PLS)回归模型，通过回归系数的稳定性进行波长变量选择。GA是对初始群体通过自然遗传机制进行选择、交换、突变等算子操作产生新群体的过程，过程中适应度值较优的变量被保留。由于UVE的人工随机噪声矩阵和GA的随机初始种群的随机性，导致了每次波长变量选择的结果不同，且经过UVE波长选择后剩余的波长变量依然较多[9]。

荧光区域积分结合凸点提取的波长选择方法是根据整体区域荧光强度变化选择出能够表征藻类色素的荧光特性的波长变量，同时避免选择过程的随机性。本文以抑食金球藻、小球藻、细长聚球藻为研究对象，在对3种藻的原始光谱进行预处理后，结合荧光区域积分和二元凸函数判定方法，进行特征区域波长变量筛选，并将筛选后的光谱数据作为输入建立PLS回归模型，验证其有效性。

1 基于凸点提取的波长选择原理

特征区域下的凸点提取是通过在主要荧光区域中提取曲面峰值点及峰值周围的点来达到去除荧光光谱冗余信息的目的。其整体过程分为两步，首先通过荧光区域积分[10]找到荧光特征区域，通过阈值改变荧光特征区域的大小；然后在特征区域内投票统计所有波长点进行凸点提取后的得分，设定阈值衡量得分可靠性。调节特征区域和凸点投票统计后得分的阈值，大于阈值的波长变量被认为利于模型建立，予以保留。其实现过程包括：

Step1：识别光谱数据中的凸点。三维荧光光谱可以看作是一个关于发射波长(x轴)、激发波长(y轴)、荧光强度(z轴)的二元函数z=f(x,y)，采用Savitzky-Golay多项式拟合微分法求偏导，计算x的偏导数时固定y轴，计算y的偏导数时固定x轴。根据二元凸函数[11]判定定理，三维光谱的任一点(x,y)为凸点需要满足以下条件

(1)

(2)

对所有样本三维荧光光谱的每个波长点的凸点进行统计，将得到的统计值与样本总数的比值作为个波长点凸点得分。

Step2：将整个荧光区域划分为若干个单位区域，由于三维荧光光谱数据是离散数据，其单位区域的体积积分为

(3)

式(3)中，I(λEX,λEM)为在激发波长λEX、发射波长λEM处的荧光强度，ΔλEX为激发波长间隔(取5 nm)，ΔλEM为发射波长间隔(取5 nm)。

Step3：求光谱数据的凸点得分与特征区域的交集，得到特征区域下的光谱数据凸点得分结果。

Step4：设定特征区域下的光谱数据凸点得分阈值ηt(0.3<ηt<1), 小于设定阈值，则认为该凸点是不利于建模的波长变量或者是由于噪声引起的波长变量。

Step5：将选择的波长变量建立偏最小二乘回归模型，采用留一交叉验证法计算内部交叉验证均方根误差(RMSECV)。

Step6：遍历所有ηq值，重复Step2—Step5过程，当RMSECV值最小时所选择的波长变量为最终的波长变量选择结果

2 实验部分

所用的抑食金球藻和小球藻由秦皇岛海洋环境监测中心站提供，细长聚球藻购自中国科学院淡水藻种库。采用FS920三维荧光光谱仪对三藻种进行光谱测量，测量条件：激发波长范围为400～650 nm，发射波长范围为630～730 nm，激发步长和发射步长均为5 nm，狭缝宽度为 5 nm，设置发射波长滞后激发波长10 nm，用1 cm石英比色池每天取样一次进行光谱测量。

实验中所选藻种的培养周期均为14 d，为采集到更多丰度下的光谱信息，每种藻培养3份，每份保持不同起始丰度。经过多周期测量，最终获得藻种光谱数据共306个。不同色素有不同的荧光特性，藻类中色素的种类决定荧光峰的位置，色素的含量决定荧光峰的强弱，图1给出了三种藻类的等高线谱。图1(a)中440/685，450/685和470/685 nm三处的荧光峰分别由抑食金球藻的叶绿素a(Chla)、叶绿素c(Chlc)和叶绿素b(Chlb)造成；图1(b)中440/685，470/685，480/685 nm三处的荧光峰分别由小球藻的Chla、Chlb和类胡萝卜素造成，420/685 nm处的荧光峰为Chla的次吸收带引起的次荧光峰；图1(c)中420/680和440/680 nm分别为细长聚球藻中Chla的次荧光锋和主荧光峰，620/655和620/680 nm处的荧光峰均由藻蓝蛋白(PC)引起，但PC对光的吸收和释放特性主要主要表现在620/655 nm处的强荧光峰。Chla存在于三种藻中，但在细长聚球藻中Chla的荧光发射波长为680 nm，可能是由于细长聚球藻中的Chla与蛋白质结合，改变了荧光的响应特性，从而引起发射光波长的偏移。

图1 三藻种等高线光谱图

3 结果与讨论

3.1 光谱数据去噪

由于一些因素的影响，通过光谱仪获得的荧光光谱既有有用信号，同时也有噪声信号。本文通过Savitzky-Golay卷积平滑法去除噪声信号，平滑窗口宽度为5，多项式次数为2。图2和图3为三种藻类样本激发光谱平滑前后对比图，可看出各藻种光谱的大体形状没有发生变化，消除了光谱上的毛刺，将误差信息重新分配到整个光谱中，使光谱更加平滑，降低了噪声干扰。

图2 三藻种平滑去噪前的激发光谱图

图3 三藻种平滑去噪后的激发光谱图

3.2 异常样本数据剔除

在光谱采集过程中，由于测量人员、测量环境、测量仪器等因素的影响，导致所测量的样本中存在一些光谱异常样本和化学值异常样本，这些异常数据会影响模型的精度，需要在建模之前剔除[12]。

光谱异常样本通过马氏距离法剔除。计算出各个样本到平均光谱的马氏距离，通过平均值和标准差来设定阈值, 调节阈值权重系数来改变阈值范围，剔除不利于建模的样本。光谱异常样本剔除结果如图4所示，本文中抑食金球藻、小球藻、细长聚球藻的阈值权重系数分别取0.5，0和0.5，剔除样本数分别为17个、19个、13个。

图4 三藻种马氏距离分布图

剩余的样本中可能存在化学值异常样本，本实验中化学值指藻种浓度，可用浓度残差法剔除。通过留一交叉验证法依次获得数据集各样本的浓度残差组成浓度残差向量，设置浓度残差的F统计性检验的显著性水平阈值，浓度残差大于该阈值的样本为化学值异常样本。浓度残差剔除异常样本的结果如图5，抑食金球藻、小球藻、细长聚球藻的显著性水平均为0.5，剔除的浓度异常样本分别为32个、34个、34个，最终可用于建模的样本数分别为53个、55个、49个。

图5 三藻种浓度残差分布图

3.3 基于荧光区域积分和凸点提取的波长选择方法

3.3.1 全波长凸点提取

对各藻类样本数据集下的全谱波长点依次进行凸点统计，如果满足凸点条件，则将该波长点保留并将凸点数加1，以每个波长点统计的凸点数与样本总数的比值作为该点的凸点统计得分，图6为各藻类凸点得分统计图，抑食金球藻、小球藻、细长聚球藻的全谱波长点均为1 071个，被判定为凸点的波长点个数分别为369个、422个和501个。

3.3.2 初始特征区域的确定

计算各藻类全谱范围下所有单位荧光区域的体积积分，初始特征区域阈值取最大单位区域体积积分的1/10，大于该阈值的单位体积积分的并集作为初始特征区域，图7为各藻类提取到的初始候选特征区域结果图，抑食金球藻、小球藻、细长聚球藻在初始特征区域下的波长点个数分别为340个、409个、525个，由图6和图7可知虽然主要的色素峰值点既被判定为凸点，又包含在初始特征区域内，但在初始候选特征区域外，有些波长点既不是色素的特征峰又无其他明显的光谱特征信息，却多次被判定为凸点。

图6 三藻种凸点统计图

3.3.3 特征区域和凸点提取阈值的确定

图8为初始特征区域下的凸点得分，抑食金球藻、小球藻、细长聚球藻在该区域下被判定为凸点的波长个数分别为131个、142个、259个，已有效去除不含特征信息的凸点，但此时由噪声引起的光谱凸点波长变量包含在内，同时由于特征区域阈值选取问题，初始特征区域并不精确，区域内可能存在与藻种信息相关性不大的波长点，依旧有光谱冗余。图9为阈值nq和nt取不同值时RMSECV的变化情况，RMSECV值越小，模型的预测能力越好，合适的nq和nt值可以提高模型的预测能力。图9(a)中nq取1/8，1/7，1/6和1/5时RMSECV值相等，图9(b)中nq取1/4和1/3时RMSECV值相等，是由于在这几个特征区域阈值下，光谱的积分区域发生了变化，但波长变量的数量并没有改变，这时模型的预测能力由凸点的阈值nt决定。当抑食金球藻、小球藻、细长聚球藻的候选特征区域阈值分别为1/10，1/6和1/8，凸点得分阈值分别为0.4，0.5和0.9时RMSECV最佳。图10为最佳候选特征区域阈值和最佳凸点得分阈值下各藻种筛选波长的最终结果，抑食金球藻、小球藻、细长聚球藻剩余波长点数分别为77，75和67。

图9 ηq和ηt取不同值时RMSECV的变化曲线

图10 特征区域结合凸点提取波长选择结果

3.4 建模方法及模型评价

分别用全谱波长数据、特征区域结合凸点提取后波长数据、UVE提取后波长数据建立PLS回归模型，采用留一交叉验证法进行内部交叉验证，以内部交叉验证决定系数(R2)、内部交叉验证均方根误差(RMSECV)作为评价指标评估波长选择的效果，其中R2反映模型验证的稳定性，RMSECV反映模型的预测能力，R2越接近1，RMSECV越接近0，模型的预测性能及稳定性越好。图11为经过UVE提取的波长变量结果，表1为各PLS回归模型的R2、RMSECV值。抑食金球藻经该方法波长筛选后的R2相对全谱提高了0.016 4，RMSECV降低了1.8×105，波长点由全谱的1 071个减少到77个，经UVE波长筛选后的R2相对全谱提高了0.001 9，RMSECV降低了3.0×104，波长点由全谱的1 071个减少到676个；小球藻经该方法波长筛选后的R2相对全谱提高了0.002，RMSECV降低了2.0×105，波长点由全谱的1 071个减少到75个，经UVE波长筛选后的R2相对全谱提高了0.001 6，RMSECV降低了1.3×105，波长点由全谱的1 071个减少到432个；细长聚球藻经该方法波长筛选后的R2相对全谱提高了0.032 4，RMSECV降低了2.6×105，波长点由全谱的1 071个减少到67个，经UVE波长筛选后的R2相对全谱提高了0.010 1，RMSECV降低了1.0×105，波长点由全谱的1 071个减少到384个。由以上分析可知，三藻种在经过特征区域结合凸点提取后R2均得到提高，其中小球藻的R2提高效果并不显著，但在没有降低模型预测能力的基础上波长变量减小到原来的7%，该波长提取方法是有效的。相比UVE，特征区域结合凸点提取所得到的模型精度更高，所选择的波长变量更少，且均分布在含有荧光特征的区域，保留了所含色素的特征峰。

表1 三种方法的R2和RMSECV

图11 UVE波长选择结果

4 结 论

利用实验得到的三种藻类不同浓度下的三维荧光光谱数据，提出了特征区域结合凸点提取的波长变量选择方法。根据藻类光谱找到包含荧光峰的初始特征区域，利用二元凸函数判定定理统计初始特征区域下的凸点，设定特征区域和凸点统计结果的不同阈值，通过计算PLS回归模型的RMSECV值衡量特征区域下波长变量凸点提取的可靠性从而选择波长去留。结果表明，用该方法得到的波长变量建立的PLS模型优于用全谱波长变量建立的PLS模型，且选择的波长变量在色素荧光峰附近，能够很好地表征藻类色素的荧光特性，为三维荧光光谱分析提供了依据。