張永彬，朱丹丹，陳 穎*，劉 喆，段瑋靚，李少華

1.燕山大學(xué)電氣工程學(xué)院河北省測(cè)試計(jì)量技術(shù)及儀器重點(diǎn)實(shí)驗(yàn)室，河北 秦皇島 066004 2.河北先河環(huán)?？萍脊煞萦邢薰?，河北 石家莊 050000

引 言

近年來我國(guó)工農(nóng)業(yè)發(fā)展迅速，海水富營(yíng)養(yǎng)化加劇，水體中浮游藻類大量繁殖，頻繁引發(fā)藻華現(xiàn)象，褐潮成因藻抑食金球藻、常見藻華中的小球藻和細(xì)長(zhǎng)聚球藻已經(jīng)成為海洋生態(tài)環(huán)境的主要污染，影響人類健康安全和環(huán)境發(fā)展，所以藻類種類和數(shù)量的測(cè)量對(duì)生態(tài)環(huán)境檢測(cè)具有重要意義[1-2]。浮游植物中含有的色素具有熒光效應(yīng)，因此熒光光譜分析是一種有效的測(cè)量方法，其中三維熒光光譜法在藻類測(cè)量中應(yīng)用廣泛，如何選擇光譜波長(zhǎng)變量從而提高三維熒光光譜的信噪比是目前一個(gè)重要研究方向。

三維熒光光譜是在已設(shè)置的激發(fā)波長(zhǎng)和發(fā)射波長(zhǎng)范圍下掃描所有波長(zhǎng)點(diǎn)而得到的數(shù)據(jù)矩陣，其數(shù)據(jù)較多，但并不是光譜區(qū)域下的所有波長(zhǎng)變量都含有有用信息，只有少數(shù)光譜區(qū)域含有物質(zhì)特征。若冗余區(qū)域波長(zhǎng)變量參與建模會(huì)影響預(yù)測(cè)精度，增加模型的計(jì)算，光譜預(yù)處理雖然可以提高信噪比，但并沒有去除不相關(guān)的波長(zhǎng)變量，不能最大程度提高模型的精度，沒有減小模型復(fù)雜度。

傳統(tǒng)的波長(zhǎng)選擇方法是峰值法，對(duì)固定波長(zhǎng)點(diǎn)下的熒光峰進(jìn)行分析，但該方法所選擇的波長(zhǎng)變量過少，原光譜的數(shù)據(jù)信息沒有得到充分利用，無法有效反映整體區(qū)域熒光強(qiáng)度變化情況[3]。目前的波長(zhǎng)選擇方法有無信息變量消除法(uninformative variables elimination，UVE)[4]、連續(xù)投影算法(successive projections algorithm，SPA)[5-6]、遺傳算法(genetic algorithm，GA)[7-8]等。SPA通過投影方式選取線性關(guān)系最小的波長(zhǎng)組合，減少光譜信息中的波長(zhǎng)變量，但SPA所選擇的波長(zhǎng)變量的個(gè)數(shù)不能超過所用樣品數(shù)量，其受樣品數(shù)量的限制較大。UVE將隨機(jī)噪聲矩陣加入光譜矩陣并建立偏最小二乘(partial least-square，PLS)回歸模型，通過回歸系數(shù)的穩(wěn)定性進(jìn)行波長(zhǎng)變量選擇。GA是對(duì)初始群體通過自然遺傳機(jī)制進(jìn)行選擇、交換、突變等算子操作產(chǎn)生新群體的過程，過程中適應(yīng)度值較優(yōu)的變量被保留。由于UVE的人工隨機(jī)噪聲矩陣和GA的隨機(jī)初始種群的隨機(jī)性，導(dǎo)致了每次波長(zhǎng)變量選擇的結(jié)果不同，且經(jīng)過UVE波長(zhǎng)選擇后剩余的波長(zhǎng)變量依然較多[9]。

熒光區(qū)域積分結(jié)合凸點(diǎn)提取的波長(zhǎng)選擇方法是根據(jù)整體區(qū)域熒光強(qiáng)度變化選擇出能夠表征藻類色素的熒光特性的波長(zhǎng)變量，同時(shí)避免選擇過程的隨機(jī)性。本文以抑食金球藻、小球藻、細(xì)長(zhǎng)聚球藻為研究對(duì)象，在對(duì)3種藻的原始光譜進(jìn)行預(yù)處理后，結(jié)合熒光區(qū)域積分和二元凸函數(shù)判定方法，進(jìn)行特征區(qū)域波長(zhǎng)變量篩選，并將篩選后的光譜數(shù)據(jù)作為輸入建立PLS回歸模型，驗(yàn)證其有效性。

1 基于凸點(diǎn)提取的波長(zhǎng)選擇原理

特征區(qū)域下的凸點(diǎn)提取是通過在主要熒光區(qū)域中提取曲面峰值點(diǎn)及峰值周圍的點(diǎn)來達(dá)到去除熒光光譜冗余信息的目的。其整體過程分為兩步，首先通過熒光區(qū)域積分[10]找到熒光特征區(qū)域，通過閾值改變熒光特征區(qū)域的大?。蝗缓笤谔卣鲄^(qū)域內(nèi)投票統(tǒng)計(jì)所有波長(zhǎng)點(diǎn)進(jìn)行凸點(diǎn)提取后的得分，設(shè)定閾值衡量得分可靠性。調(diào)節(jié)特征區(qū)域和凸點(diǎn)投票統(tǒng)計(jì)后得分的閾值，大于閾值的波長(zhǎng)變量被認(rèn)為利于模型建立，予以保留。其實(shí)現(xiàn)過程包括：

Step1：識(shí)別光譜數(shù)據(jù)中的凸點(diǎn)。三維熒光光譜可以看作是一個(gè)關(guān)于發(fā)射波長(zhǎng)(x軸)、激發(fā)波長(zhǎng)(y軸)、熒光強(qiáng)度(z軸)的二元函數(shù)z=f(x,y)，采用Savitzky-Golay多項(xiàng)式擬合微分法求偏導(dǎo)，計(jì)算x的偏導(dǎo)數(shù)時(shí)固定y軸，計(jì)算y的偏導(dǎo)數(shù)時(shí)固定x軸。根據(jù)二元凸函數(shù)[11]判定定理，三維光譜的任一點(diǎn)(x,y)為凸點(diǎn)需要滿足以下條件

(1)

(2)

對(duì)所有樣本三維熒光光譜的每個(gè)波長(zhǎng)點(diǎn)的凸點(diǎn)進(jìn)行統(tǒng)計(jì)，將得到的統(tǒng)計(jì)值與樣本總數(shù)的比值作為個(gè)波長(zhǎng)點(diǎn)凸點(diǎn)得分。

Step2：將整個(gè)熒光區(qū)域劃分為若干個(gè)單位區(qū)域，由于三維熒光光譜數(shù)據(jù)是離散數(shù)據(jù)，其單位區(qū)域的體積積分為

(3)

式(3)中，I(λEX,λEM)為在激發(fā)波長(zhǎng)λEX、發(fā)射波長(zhǎng)λEM處的熒光強(qiáng)度，ΔλEX為激發(fā)波長(zhǎng)間隔(取5 nm)，ΔλEM為發(fā)射波長(zhǎng)間隔(取5 nm)。

Step3：求光譜數(shù)據(jù)的凸點(diǎn)得分與特征區(qū)域的交集，得到特征區(qū)域下的光譜數(shù)據(jù)凸點(diǎn)得分結(jié)果。

Step4：設(shè)定特征區(qū)域下的光譜數(shù)據(jù)凸點(diǎn)得分閾值ηt(0.3<ηt<1), 小于設(shè)定閾值，則認(rèn)為該凸點(diǎn)是不利于建模的波長(zhǎng)變量或者是由于噪聲引起的波長(zhǎng)變量。

Step5：將選擇的波長(zhǎng)變量建立偏最小二乘回歸模型，采用留一交叉驗(yàn)證法計(jì)算內(nèi)部交叉驗(yàn)證均方根誤差(RMSECV)。

Step6：遍歷所有ηq值，重復(fù)Step2—Step5過程，當(dāng)RMSECV值最小時(shí)所選擇的波長(zhǎng)變量為最終的波長(zhǎng)變量選擇結(jié)果

2 實(shí)驗(yàn)部分

所用的抑食金球藻和小球藻由秦皇島海洋環(huán)境監(jiān)測(cè)中心站提供，細(xì)長(zhǎng)聚球藻購(gòu)自中國(guó)科學(xué)院淡水藻種庫(kù)。采用FS920三維熒光光譜儀對(duì)三藻種進(jìn)行光譜測(cè)量，測(cè)量條件：激發(fā)波長(zhǎng)范圍為400～650 nm，發(fā)射波長(zhǎng)范圍為630～730 nm，激發(fā)步長(zhǎng)和發(fā)射步長(zhǎng)均為5 nm，狹縫寬度為 5 nm，設(shè)置發(fā)射波長(zhǎng)滯后激發(fā)波長(zhǎng)10 nm，用1 cm石英比色池每天取樣一次進(jìn)行光譜測(cè)量。

實(shí)驗(yàn)中所選藻種的培養(yǎng)周期均為14 d，為采集到更多豐度下的光譜信息，每種藻培養(yǎng)3份，每份保持不同起始豐度。經(jīng)過多周期測(cè)量，最終獲得藻種光譜數(shù)據(jù)共306個(gè)。不同色素有不同的熒光特性，藻類中色素的種類決定熒光峰的位置，色素的含量決定熒光峰的強(qiáng)弱，圖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分別為細(xì)長(zhǎng)聚球藻中Chla的次熒光鋒和主熒光峰，620/655和620/680 nm處的熒光峰均由藻藍(lán)蛋白(PC)引起，但PC對(duì)光的吸收和釋放特性主要主要表現(xiàn)在620/655 nm處的強(qiáng)熒光峰。Chla存在于三種藻中，但在細(xì)長(zhǎng)聚球藻中Chla的熒光發(fā)射波長(zhǎng)為680 nm，可能是由于細(xì)長(zhǎng)聚球藻中的Chla與蛋白質(zhì)結(jié)合，改變了熒光的響應(yīng)特性，從而引起發(fā)射光波長(zhǎng)的偏移。

圖1 三藻種等高線光譜圖

3 結(jié)果與討論

3.1 光譜數(shù)據(jù)去噪

由于一些因素的影響，通過光譜儀獲得的熒光光譜既有有用信號(hào)，同時(shí)也有噪聲信號(hào)。本文通過Savitzky-Golay卷積平滑法去除噪聲信號(hào)，平滑窗口寬度為5，多項(xiàng)式次數(shù)為2。圖2和圖3為三種藻類樣本激發(fā)光譜平滑前后對(duì)比圖，可看出各藻種光譜的大體形狀沒有發(fā)生變化，消除了光譜上的毛刺，將誤差信息重新分配到整個(gè)光譜中，使光譜更加平滑，降低了噪聲干擾。

圖2 三藻種平滑去噪前的激發(fā)光譜圖

圖3 三藻種平滑去噪后的激發(fā)光譜圖

3.2 異常樣本數(shù)據(jù)剔除

在光譜采集過程中，由于測(cè)量人員、測(cè)量環(huán)境、測(cè)量?jī)x器等因素的影響，導(dǎo)致所測(cè)量的樣本中存在一些光譜異常樣本和化學(xué)值異常樣本，這些異常數(shù)據(jù)會(huì)影響模型的精度，需要在建模之前剔除[12]。

光譜異常樣本通過馬氏距離法剔除。計(jì)算出各個(gè)樣本到平均光譜的馬氏距離，通過平均值和標(biāo)準(zhǔn)差來設(shè)定閾值, 調(diào)節(jié)閾值權(quán)重系數(shù)來改變閾值范圍，剔除不利于建模的樣本。光譜異常樣本剔除結(jié)果如圖4所示，本文中抑食金球藻、小球藻、細(xì)長(zhǎng)聚球藻的閾值權(quán)重系數(shù)分別取0.5，0和0.5，剔除樣本數(shù)分別為17個(gè)、19個(gè)、13個(gè)。

圖4 三藻種馬氏距離分布圖

剩余的樣本中可能存在化學(xué)值異常樣本，本實(shí)驗(yàn)中化學(xué)值指藻種濃度，可用濃度殘差法剔除。通過留一交叉驗(yàn)證法依次獲得數(shù)據(jù)集各樣本的濃度殘差組成濃度殘差向量，設(shè)置濃度殘差的F統(tǒng)計(jì)性檢驗(yàn)的顯著性水平閾值，濃度殘差大于該閾值的樣本為化學(xué)值異常樣本。濃度殘差剔除異常樣本的結(jié)果如圖5，抑食金球藻、小球藻、細(xì)長(zhǎng)聚球藻的顯著性水平均為0.5，剔除的濃度異常樣本分別為32個(gè)、34個(gè)、34個(gè)，最終可用于建模的樣本數(shù)分別為53個(gè)、55個(gè)、49個(gè)。

圖5 三藻種濃度殘差分布圖

3.3 基于熒光區(qū)域積分和凸點(diǎn)提取的波長(zhǎng)選擇方法

3.3.1 全波長(zhǎng)凸點(diǎn)提取

對(duì)各藻類樣本數(shù)據(jù)集下的全譜波長(zhǎng)點(diǎn)依次進(jìn)行凸點(diǎn)統(tǒng)計(jì)，如果滿足凸點(diǎn)條件，則將該波長(zhǎng)點(diǎn)保留并將凸點(diǎn)數(shù)加1，以每個(gè)波長(zhǎng)點(diǎn)統(tǒng)計(jì)的凸點(diǎn)數(shù)與樣本總數(shù)的比值作為該點(diǎn)的凸點(diǎn)統(tǒng)計(jì)得分，圖6為各藻類凸點(diǎn)得分統(tǒng)計(jì)圖，抑食金球藻、小球藻、細(xì)長(zhǎng)聚球藻的全譜波長(zhǎng)點(diǎn)均為1 071個(gè)，被判定為凸點(diǎn)的波長(zhǎng)點(diǎn)個(gè)數(shù)分別為369個(gè)、422個(gè)和501個(gè)。

3.3.2 初始特征區(qū)域的確定

計(jì)算各藻類全譜范圍下所有單位熒光區(qū)域的體積積分，初始特征區(qū)域閾值取最大單位區(qū)域體積積分的1/10，大于該閾值的單位體積積分的并集作為初始特征區(qū)域，圖7為各藻類提取到的初始候選特征區(qū)域結(jié)果圖，抑食金球藻、小球藻、細(xì)長(zhǎng)聚球藻在初始特征區(qū)域下的波長(zhǎng)點(diǎn)個(gè)數(shù)分別為340個(gè)、409個(gè)、525個(gè)，由圖6和圖7可知雖然主要的色素峰值點(diǎn)既被判定為凸點(diǎn)，又包含在初始特征區(qū)域內(nèi)，但在初始候選特征區(qū)域外，有些波長(zhǎng)點(diǎn)既不是色素的特征峰又無其他明顯的光譜特征信息，卻多次被判定為凸點(diǎn)。

圖6 三藻種凸點(diǎn)統(tǒng)計(jì)圖

3.3.3 特征區(qū)域和凸點(diǎn)提取閾值的確定

圖8為初始特征區(qū)域下的凸點(diǎn)得分，抑食金球藻、小球藻、細(xì)長(zhǎng)聚球藻在該區(qū)域下被判定為凸點(diǎn)的波長(zhǎng)個(gè)數(shù)分別為131個(gè)、142個(gè)、259個(gè)，已有效去除不含特征信息的凸點(diǎn)，但此時(shí)由噪聲引起的光譜凸點(diǎn)波長(zhǎng)變量包含在內(nèi)，同時(shí)由于特征區(qū)域閾值選取問題，初始特征區(qū)域并不精確，區(qū)域內(nèi)可能存在與藻種信息相關(guān)性不大的波長(zhǎng)點(diǎn)，依舊有光譜冗余。圖9為閾值nq和nt取不同值時(shí)RMSECV的變化情況，RMSECV值越小，模型的預(yù)測(cè)能力越好，合適的nq和nt值可以提高模型的預(yù)測(cè)能力。圖9(a)中nq取1/8，1/7，1/6和1/5時(shí)RMSECV值相等，圖9(b)中nq取1/4和1/3時(shí)RMSECV值相等，是由于在這幾個(gè)特征區(qū)域閾值下，光譜的積分區(qū)域發(fā)生了變化，但波長(zhǎng)變量的數(shù)量并沒有改變，這時(shí)模型的預(yù)測(cè)能力由凸點(diǎn)的閾值nt決定。當(dāng)抑食金球藻、小球藻、細(xì)長(zhǎng)聚球藻的候選特征區(qū)域閾值分別為1/10，1/6和1/8，凸點(diǎn)得分閾值分別為0.4，0.5和0.9時(shí)RMSECV最佳。圖10為最佳候選特征區(qū)域閾值和最佳凸點(diǎn)得分閾值下各藻種篩選波長(zhǎng)的最終結(jié)果，抑食金球藻、小球藻、細(xì)長(zhǎng)聚球藻剩余波長(zhǎng)點(diǎn)數(shù)分別為77，75和67。

圖9 ηq和ηt取不同值時(shí)RMSECV的變化曲線

圖10 特征區(qū)域結(jié)合凸點(diǎn)提取波長(zhǎng)選擇結(jié)果

3.4 建模方法及模型評(píng)價(jià)

分別用全譜波長(zhǎng)數(shù)據(jù)、特征區(qū)域結(jié)合凸點(diǎn)提取后波長(zhǎng)數(shù)據(jù)、UVE提取后波長(zhǎng)數(shù)據(jù)建立PLS回歸模型，采用留一交叉驗(yàn)證法進(jìn)行內(nèi)部交叉驗(yàn)證，以內(nèi)部交叉驗(yàn)證決定系數(shù)(R2)、內(nèi)部交叉驗(yàn)證均方根誤差(RMSECV)作為評(píng)價(jià)指標(biāo)評(píng)估波長(zhǎng)選擇的效果，其中R2反映模型驗(yàn)證的穩(wěn)定性，RMSECV反映模型的預(yù)測(cè)能力，R2越接近1，RMSECV越接近0，模型的預(yù)測(cè)性能及穩(wěn)定性越好。圖11為經(jīng)過UVE提取的波長(zhǎng)變量結(jié)果，表1為各PLS回歸模型的R2、RMSECV值。抑食金球藻經(jīng)該方法波長(zhǎng)篩選后的R2相對(duì)全譜提高了0.016 4，RMSECV降低了1.8×105，波長(zhǎng)點(diǎn)由全譜的1 071個(gè)減少到77個(gè)，經(jīng)UVE波長(zhǎng)篩選后的R2相對(duì)全譜提高了0.001 9，RMSECV降低了3.0×104，波長(zhǎng)點(diǎn)由全譜的1 071個(gè)減少到676個(gè)；小球藻經(jīng)該方法波長(zhǎng)篩選后的R2相對(duì)全譜提高了0.002，RMSECV降低了2.0×105，波長(zhǎng)點(diǎn)由全譜的1 071個(gè)減少到75個(gè)，經(jīng)UVE波長(zhǎng)篩選后的R2相對(duì)全譜提高了0.001 6，RMSECV降低了1.3×105，波長(zhǎng)點(diǎn)由全譜的1 071個(gè)減少到432個(gè)；細(xì)長(zhǎng)聚球藻經(jīng)該方法波長(zhǎng)篩選后的R2相對(duì)全譜提高了0.032 4，RMSECV降低了2.6×105，波長(zhǎng)點(diǎn)由全譜的1 071個(gè)減少到67個(gè)，經(jīng)UVE波長(zhǎng)篩選后的R2相對(duì)全譜提高了0.010 1，RMSECV降低了1.0×105，波長(zhǎng)點(diǎn)由全譜的1 071個(gè)減少到384個(gè)。由以上分析可知，三藻種在經(jīng)過特征區(qū)域結(jié)合凸點(diǎn)提取后R2均得到提高，其中小球藻的R2提高效果并不顯著，但在沒有降低模型預(yù)測(cè)能力的基礎(chǔ)上波長(zhǎng)變量減小到原來的7%，該波長(zhǎng)提取方法是有效的。相比UVE，特征區(qū)域結(jié)合凸點(diǎn)提取所得到的模型精度更高，所選擇的波長(zhǎng)變量更少，且均分布在含有熒光特征的區(qū)域，保留了所含色素的特征峰。

表1 三種方法的R2和RMSECV

圖11 UVE波長(zhǎng)選擇結(jié)果

4 結(jié) 論

利用實(shí)驗(yàn)得到的三種藻類不同濃度下的三維熒光光譜數(shù)據(jù)，提出了特征區(qū)域結(jié)合凸點(diǎn)提取的波長(zhǎng)變量選擇方法。根據(jù)藻類光譜找到包含熒光峰的初始特征區(qū)域，利用二元凸函數(shù)判定定理統(tǒng)計(jì)初始特征區(qū)域下的凸點(diǎn)，設(shè)定特征區(qū)域和凸點(diǎn)統(tǒng)計(jì)結(jié)果的不同閾值，通過計(jì)算PLS回歸模型的RMSECV值衡量特征區(qū)域下波長(zhǎng)變量凸點(diǎn)提取的可靠性從而選擇波長(zhǎng)去留。結(jié)果表明，用該方法得到的波長(zhǎng)變量建立的PLS模型優(yōu)于用全譜波長(zhǎng)變量建立的PLS模型，且選擇的波長(zhǎng)變量在色素?zé)晒夥甯浇?，能夠很好地表征藻類色素的熒光特性，為三維熒光光譜分析提供了依據(jù)。