一種結(jié)合直接正交信號(hào)校正與蒙特卡羅的波長(zhǎng)選擇方法

謝林江， 洪明堅(jiān)， 余志榮

重慶大學(xué)大數(shù)據(jù)與軟件學(xué)院， 重慶 401331

引 言

近紅外光譜分析技術(shù)因?yàn)槠浞治鏊俣瓤臁?無(wú)損、 無(wú)污染和成本低等優(yōu)點(diǎn)， 在農(nóng)業(yè)、 制藥、 石化、 環(huán)境和醫(yī)學(xué)等領(lǐng)域有著廣泛的應(yīng)用[1]。 近紅外光譜也有波長(zhǎng)點(diǎn)多、 譜峰重疊、 吸光強(qiáng)度低、 波長(zhǎng)點(diǎn)之間共線(xiàn)性嚴(yán)重等特點(diǎn)。 采用全光譜數(shù)據(jù)進(jìn)行分析和建模， 不但會(huì)增加模型的運(yùn)算時(shí)間， 也會(huì)增加模型的復(fù)雜度， 其中部分波長(zhǎng)點(diǎn)還會(huì)降低模型的預(yù)測(cè)質(zhì)量[2]。 選擇出有效的波長(zhǎng)點(diǎn)參與建模有著重要的實(shí)際意義。

為了有效選出合適的波長(zhǎng)點(diǎn)， 基于不同原理的波長(zhǎng)選擇算法先后被提出。 主要有蒙特卡羅無(wú)信息變量消除法(Monte Carlo uninformative variable elimination, MCUVE)[3]、 遺傳算法(genetic algorithms, GA)[4]、 連續(xù)投影算法(successive projections algorithms, SPA)[5]、 協(xié)同間隔偏最小二乘法(synergy interval PLS, SiPLS)[6]、 向后間隔偏最小二乘法(backward interval partial least square, B-iPLS)[7]、 變量迭代空間收縮算法(variable iterative space shrinkage approach, VISSA)[8]、 競(jìng)爭(zhēng)性自適應(yīng)權(quán)重取樣法(competitive adaptive reweighted sampling， CARS)[9]、 隨機(jī)蛙跳(random frog, RF)[10]和最小絕對(duì)收縮和選擇算法(least absolute shrinkage and selection operator, LASSO)等。 基于PLS回歸系數(shù)為波長(zhǎng)重要性指標(biāo)進(jìn)行波長(zhǎng)篩選是一類(lèi)常用方法[11-13]， 如MCUVE， VISSA， CARS等。 但是， 基于回歸系數(shù)“重要性”進(jìn)行選擇的方法存在兩個(gè)問(wèn)題： ①由于光譜容易受到噪聲和儀器測(cè)量誤差的影響， 回歸系數(shù)并非總是體現(xiàn)波長(zhǎng)“重要性”的真實(shí)信息。 ②利用回歸系數(shù)選擇波長(zhǎng)有著主因子個(gè)數(shù)選取的困難， 不同的主因子個(gè)數(shù)對(duì)應(yīng)的回歸系數(shù)往往有著較大的不穩(wěn)定性。

針對(duì)上述問(wèn)題， 提出一種通過(guò)濾除“不重要”波長(zhǎng)進(jìn)行選擇的新方法。 首先用直接正交信號(hào)校正(direct orthogonal signal correction， DOSC)提取光譜中與物質(zhì)濃度陣無(wú)關(guān)的信息， 將其權(quán)重向量絕對(duì)值歸一化后作為波長(zhǎng)被篩除的概率， 然后用Monte Carlo方法進(jìn)行迭代， 逐步剔除不重要的波長(zhǎng)。 通過(guò)玉米和汽油數(shù)據(jù)集對(duì)MC-DOSC進(jìn)行了驗(yàn)證， 并與MCUVE, GA和CARS三種算法進(jìn)行了對(duì)比， 驗(yàn)證MC-DOSC的有效性。

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

1.1 數(shù)據(jù)來(lái)源

1.1.1 玉米數(shù)據(jù)集

玉米數(shù)據(jù)集包含玉米近紅外光譜及其植物油(oil)的含量， 此數(shù)據(jù)集是80個(gè)玉米樣本在m5， mp5和mp6三臺(tái)不同光譜儀采集得到。 實(shí)驗(yàn)采用的光譜由mp5光譜儀采樣得到。 波長(zhǎng)采集范圍為1 100～2 498 nm， 采集間隔為2 nm, 共700個(gè)波長(zhǎng)點(diǎn)。 使用duplex方法將樣本分成60個(gè)訓(xùn)練集樣本和20個(gè)驗(yàn)證集樣本。 數(shù)據(jù)集可從http://www.eigenvector.com/data下載。

1.1.2 汽油數(shù)據(jù)集

汽油(gasoline)數(shù)據(jù)集[14]包含汽油的近紅外光譜及其辛烷值(octane)。 此數(shù)據(jù)集含有60個(gè)樣本， 波長(zhǎng)采集范圍為900～1 700 nm， 采集間隔為2 nm， 由于其前100個(gè)波長(zhǎng)(900～1 100 nm)幾乎不攜帶有效信息， 故提前將其刪除。 使用duplex方法將樣本分成45個(gè)訓(xùn)練集樣本和15個(gè)驗(yàn)證集樣本。

1.2 正交信號(hào)校正

正交信號(hào)校正方法(OSC)通常用于光譜預(yù)處理， 核心原理是將光譜陣X和濃度陣Y進(jìn)行正交， 扣除光譜中與濃度無(wú)關(guān)的信息。 當(dāng)光譜和濃度的相關(guān)性不大或者光譜中背景噪聲太大時(shí)， 使用OSC方法可以有效減少PLS模型的主因子個(gè)數(shù)， 降低模型的復(fù)雜度， 增強(qiáng)模型的穩(wěn)健性以及預(yù)測(cè)能力。

根據(jù)扣除方式的不同， 衍生了多種OSC算法， 包括Wold-OSC、 DO(direct orthogonal)、 Fearn-OSC、 O-PLS(orthogonal projects to latent structure)、 DOSC等。 大部分OSC算法， 需要選擇主因子個(gè)數(shù)， 然后通過(guò)迭代逐步扣除。 而DOSC算法可以只用一個(gè)主因子扣除光譜中與濃度無(wú)關(guān)的信息， 避免了主因子個(gè)數(shù)選擇的問(wèn)題[15]。 因此， 采用DOSC提取波長(zhǎng)不重要的度量信息。

DOSC首先提取濃度陣Y在光譜陣X所張開(kāi)的線(xiàn)性空間的投影M， 如圖1所示

M=XT((XT)-1)TY

(1)

計(jì)算X在M的正交補(bǔ)空間的投影Z， 見(jiàn)式(2)

Z=X-MM-1X

(2)

再對(duì)ZZT進(jìn)行主成分分析， 提取前k個(gè)需正交處理的主成分得分矩陣T， 再計(jì)算權(quán)重矩陣W， 見(jiàn)式(3)

W=X-1T

(3)

得到權(quán)重矩陣W。 由DOSC只需一個(gè)主因子可知， 僅需取W中第一個(gè)權(quán)重向量w來(lái)度量波長(zhǎng)的“不重要”性， 與PLS利用回歸系數(shù)β度量波長(zhǎng)的重要性有互補(bǔ)關(guān)系。 將DOSC權(quán)重向量w與PLS回歸系數(shù)β作對(duì)比， 如圖1所示。 可以看出， DOSC的權(quán)重向量和PLS的回歸系數(shù)具有明顯的互補(bǔ)關(guān)系。

圖1 DOSC權(quán)重向量和PLS回歸系數(shù)對(duì)比， 其中w和β來(lái)源于玉米數(shù)據(jù)集

1.3 MC-DOSC波長(zhǎng)選擇算法

將DOSC的權(quán)重向量w作為度量波長(zhǎng)“不重要”性的依據(jù)， 并以此選擇波長(zhǎng)。 由于近紅外光譜受到吸光度低、 譜峰重疊和噪聲等因素的影響[16]， 嚴(yán)格根據(jù)DOSC權(quán)重向量w的絕對(duì)值大小進(jìn)行波長(zhǎng)篩選一般不會(huì)得到最佳波長(zhǎng)點(diǎn)子集。 針對(duì)這個(gè)問(wèn)題， 將權(quán)重向量w絕對(duì)值歸一化后作為波長(zhǎng)點(diǎn)被濾除的概率， 絕對(duì)值越小， 對(duì)應(yīng)的波長(zhǎng)點(diǎn)被保留的概率越大。 選擇正態(tài)隨機(jī)分布函數(shù)進(jìn)行大量的蒙特卡羅模擬， 得到RMSECV最小時(shí)對(duì)應(yīng)的波長(zhǎng)點(diǎn)作為備選子集， 再對(duì)以上過(guò)程進(jìn)行迭代， 得到一系列備選子集， 然后從中選出RMSECV最小的作為最佳子集。 MC-DOSC波長(zhǎng)選擇算法的主要步驟如下：

Step1: 將光譜陣X與濃度陣Y進(jìn)行DOSC， 得到權(quán)重向量w；

Step2: 對(duì)權(quán)重向量w求絕對(duì)值并歸一化得到相應(yīng)的概率pj；

Step3: 初始化k=1， 生成正態(tài)分布隨機(jī)向量R～N(0, 1)， 將R求絕對(duì)值并歸一化得到rj；

Step4: 濾除波長(zhǎng)點(diǎn){j|pj>rj,j=1, 2…,n}， 其余波長(zhǎng)點(diǎn)用PLS建模， 使用十折交叉驗(yàn)證(10-fold cross-validation)計(jì)算RMSECV， 并確定主因子數(shù)；

Step5: 當(dāng)k

Step6: 在N個(gè)RMSECV中選擇最小的RMSECV所對(duì)應(yīng)的波長(zhǎng)點(diǎn)作為備選子集；

Step7: 將備選子集作為新的光譜陣Xnew代回Step1進(jìn)行迭代， 直到備選子集波長(zhǎng)點(diǎn)個(gè)數(shù)等于1， 執(zhí)行Step8;

Step8: 從所有的備選子集中， 選出RMSECV最小的備選子集作為最佳波長(zhǎng)子集。

在該算法中， Step5中的蒙特卡羅模擬次數(shù)N取500； 最大主因子數(shù)由全譜的蒙特卡羅交叉驗(yàn)證(MCCV)確定。 為了提高算法運(yùn)行效率， 同時(shí)盡量保留關(guān)鍵波長(zhǎng)， 在第一次迭代時(shí)， 進(jìn)行波長(zhǎng)粗篩， 直接剔除掉“不重要”性pj大于δ的波長(zhǎng)點(diǎn)， 為了確定δ， 以0.05為間隔， 取值范圍為0到1， 通過(guò)RMSECV確定最優(yōu)值； 在玉米數(shù)據(jù)集上，δ取0.5； 在汽油數(shù)據(jù)集上，δ取0.6； 之后進(jìn)行波長(zhǎng)精選， 用蒙特卡羅方法濾除波長(zhǎng)點(diǎn)。

1.4 模型建立及評(píng)價(jià)

將MC-DOSC得到的最佳波長(zhǎng)子集建立PLS模型， 并同MCUVE， GA和CARS三種方法所建PLS模型進(jìn)行對(duì)比， 以此驗(yàn)證MC-DOSC算法的有效性。 采用預(yù)測(cè)集相關(guān)系數(shù)R和預(yù)測(cè)均方根誤差(RMSEP)來(lái)評(píng)價(jià)模型的性能。 所有算法均在Matlab 2015b軟件中實(shí)現(xiàn)。

2 結(jié)果與討論

2.1 MC-DOSC算法的波長(zhǎng)選擇結(jié)果

2.1.1 玉米數(shù)據(jù)集

以玉米數(shù)據(jù)集的MC-DOSC波長(zhǎng)選擇過(guò)程為例， 如圖2所示。 圖2(a)展示了每一次迭代中濾除“不重要”性較大的波長(zhǎng)點(diǎn)這一過(guò)程， 第一次篩選直接剪掉“不重要”性大于δ=0.5的波長(zhǎng)點(diǎn)， 達(dá)到粗選的目的； 然后用蒙特卡羅方法進(jìn)行細(xì)選， 直到波長(zhǎng)點(diǎn)不能再減少為止。 從圖2(b)可以看出， 隨著迭代的不斷進(jìn)行， 備選子集所建模型的RMSECV逐步減小到最小值， 然后在34次迭代后迅速增大。 分析認(rèn)為RMSECV達(dá)到最小值前， 濾除的是干擾波長(zhǎng)點(diǎn)或者對(duì)模型貢獻(xiàn)很小的波長(zhǎng)點(diǎn)； 在達(dá)到最小值后， 濾除了對(duì)模型貢獻(xiàn)很大的關(guān)鍵波長(zhǎng)點(diǎn)。 將第34次迭代的備選子集作為最佳波長(zhǎng)子集。 圖2(c)為最佳波長(zhǎng)子集所建PLS模型與全光譜所建PLS模型的RMSEP對(duì)比圖， 由圖2(c)可知， 經(jīng)過(guò)MC-DOSC波長(zhǎng)篩選后， 模型的預(yù)測(cè)能力大幅度提高了。

2.1.2 汽油數(shù)據(jù)集

在汽油數(shù)據(jù)集中， 波長(zhǎng)選擇過(guò)程與玉米數(shù)據(jù)集類(lèi)似， 如圖3(a)、 圖3(b)所示， 在第16次迭代后， RMSECV達(dá)到最小值， 最終選擇第16次迭代時(shí)對(duì)應(yīng)的波長(zhǎng)點(diǎn)建立PLS模型。 從圖3(c)可以看出， 相比于全譜的PLS模型， 經(jīng)過(guò)MC-DOSC波長(zhǎng)篩選后， 模型的預(yù)測(cè)能力同樣提高了。

2.2 不同方法建模結(jié)果的對(duì)比及分析

為了驗(yàn)證MC-DOSC算法的有效性， 將MCUVE， GA， CARS和MC-DOSC 4種不同方法選擇的波長(zhǎng)點(diǎn)建立PLS模型， 然后對(duì)比不同方法所建模型預(yù)測(cè)能力的差異， 以驗(yàn)證MC-DOSC算法的有效性。

2.2.1 玉米數(shù)據(jù)集

4種方法在玉米數(shù)據(jù)集選擇的波段如圖4所示， 可以看出， 4種方法均選擇了1 700和2 300 nm附近的波長(zhǎng)點(diǎn)， 對(duì)應(yīng)著C-H鍵的伸縮振動(dòng)區(qū)域。 GA， CARS和MC-DOSC幾乎都僅選擇了1 650～2 498 nm波段的波長(zhǎng)點(diǎn)， 此波段對(duì)應(yīng)著C—H鍵的一級(jí)倍頻與合頻。 說(shuō)明對(duì)于玉米油分含量的預(yù)測(cè)， 1 650～2 498 nm波段有著關(guān)鍵作用。

圖4 4種方法選擇的波長(zhǎng)點(diǎn)

從表1可知， 相比于全光譜， 用4種不同方法篩選波長(zhǎng)后， 模型的預(yù)測(cè)能力都有一定程度的提升。 其中MCUVE選擇的波長(zhǎng)點(diǎn)最多， 預(yù)測(cè)結(jié)果也最差。 這主要是因?yàn)椋?MCUVE濾除了噪聲波長(zhǎng)點(diǎn)， 保留了不是噪聲但對(duì)模型沒(méi)有貢獻(xiàn)的波長(zhǎng)點(diǎn)。 CARS選擇了較少的波長(zhǎng)點(diǎn)， 預(yù)測(cè)結(jié)果也較好。 MC-DOSC選擇了15個(gè)波長(zhǎng)點(diǎn)， 略少于CARS。 MC-DOSC預(yù)測(cè)結(jié)果是4種方法中最好的， 這是因?yàn)椋?MC-DOSC濾除了噪聲波長(zhǎng)點(diǎn)， 也濾除了對(duì)模型沒(méi)有貢獻(xiàn)的波長(zhǎng)點(diǎn)。 相較于全光譜的PLS模型， 經(jīng)過(guò)MC-DOSC算法篩選波長(zhǎng)后， PLS模型的預(yù)測(cè)能力大幅度提高。 驗(yàn)證相關(guān)系數(shù)Rp從0.828 2提高到0.931 4， RMSEP從0.109 8減少到0.071 3。

表1 不同模型在玉米數(shù)據(jù)集預(yù)測(cè)能力的對(duì)比

2.2.2 汽油數(shù)據(jù)集

4種方法在汽油數(shù)據(jù)集選擇的波段如圖5所示， 可以看出， 4種方法選擇的波長(zhǎng)點(diǎn)數(shù)相差不大， 均選擇了1 410 nm附近的波段， 這對(duì)應(yīng)著芳烴和C—H的合頻吸收區(qū)域。 MCUVE選擇的波段較為集中， 幾乎沒(méi)有選擇1 450～1 700 nm波段的波長(zhǎng)點(diǎn)， 而這個(gè)區(qū)域是各種C—H基團(tuán)以及N—H基團(tuán)、 O—H基團(tuán)的伸縮振動(dòng)的一級(jí)倍頻吸收區(qū)域， 因此MCUVE對(duì)映的預(yù)測(cè)誤差也最大。 而CARS和MC-DOSC都選擇了其他方法沒(méi)有選擇的1 700 nm附近的波段， 這對(duì)映著甲基C—H的一級(jí)倍頻吸收區(qū)域。

圖5 4種方法選擇的波長(zhǎng)點(diǎn)

從表2可知， 盡管MCUVE選擇了最多的波長(zhǎng)， 但預(yù)測(cè)結(jié)果并不理想。 MC-DOSC選擇了最少的波長(zhǎng)， 預(yù)測(cè)結(jié)果也是最好的。 相較于全光譜PLS模型， 經(jīng)過(guò)MC-DOSC算法篩選波長(zhǎng)后， 驗(yàn)證相關(guān)系數(shù)Rp從0.987 5提高到0.993 9， RMSEP從0.255 5減小到0.178 8。

表2 不同模型在汽油數(shù)據(jù)集預(yù)測(cè)能力的對(duì)比

3 結(jié) 論

針對(duì)近紅外光譜的特性與光譜分析存在的問(wèn)題， 提出一種新的波長(zhǎng)選擇算法MC-DOSC。 該算法利用直接正交校正(DOSC)得到的權(quán)重向量w度量波長(zhǎng)的“不重要”性， 再結(jié)合蒙特卡羅方法進(jìn)行波長(zhǎng)篩選。 用MCUVE， GA和CARS 在兩個(gè)數(shù)據(jù)集上的實(shí)驗(yàn)進(jìn)行了對(duì)比， 驗(yàn)證了MC-DOSC算法的有效性。 實(shí)驗(yàn)結(jié)果表明， MC-DOSC是一種有效的波長(zhǎng)選擇算法， 具有廣泛的應(yīng)用價(jià)值。