基于特征區(qū)域下凸點提取的藻類熒光光譜波長選擇方法

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

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

引 言

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

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

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

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

1 基于凸點提取的波長選擇原理

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

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

(1)

(2)

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

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

(3)

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

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

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

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

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

2 實驗部分

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

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

圖1 三藻種等高線光譜圖

3 結(jié)果與討論

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

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

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

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

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

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

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

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

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

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

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

3.3.1 全波長凸點提取

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

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

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

圖6 三藻種凸點統(tǒng)計圖

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

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

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

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

3.4 建模方法及模型評價

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

表1 三種方法的R2和RMSECV

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

4 結(jié) 論

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