鄭開(kāi)逸， 馮雨航， 張 文， 黃曉瑋， 李志華， 張 迪， 石吉勇， 鄒小波

江蘇大學(xué)食品與生物工程學(xué)院， 江蘇 鎮(zhèn)江 212013

引 言

近紅外光譜(near infrared spectra, NIR)是一種位于760～2 500 nm的電磁波。 NIR分析方法具有檢測(cè)快速、 不破壞樣品的特點(diǎn)， 故廣泛用于食品[1-2]、 石油[3]、 制藥[4]等領(lǐng)域。 但是， NIR光譜產(chǎn)生機(jī)理復(fù)雜， 干擾較多， 因此它不可能像紫外-可見(jiàn)光譜(ultraviolet-visible spectra, UV-Vis)那樣遵循嚴(yán)格的Lambert-Beer定律， 也就難以建立理論模型來(lái)描述光譜和濃度的關(guān)系。 對(duì)此， 可以用統(tǒng)計(jì)模型建立NIR光譜(X)和理化指標(biāo)(y)之間的關(guān)系。 故在NIR分析中， 模型的建立和維護(hù)對(duì)分析結(jié)果的正確性至關(guān)重要[5-6]。

然而， 當(dāng)某一測(cè)量條件發(fā)生了變化， 就算是測(cè)量具有相同理化指標(biāo)的樣品， 得到的光譜也是有很大區(qū)別的， 這就導(dǎo)致建立好的模型預(yù)測(cè)不了發(fā)生變化的測(cè)量條件下的樣品。 常見(jiàn)導(dǎo)致相同理化指標(biāo)的光譜發(fā)生變化的原因有如下幾點(diǎn)： (1)樣品性狀的改變， 即樣品中與理化指標(biāo)無(wú)關(guān)的成分發(fā)生變化。 (2)儀器對(duì)理化指標(biāo)函數(shù)關(guān)系的變化。 (3)諸如濕度和溫度等環(huán)境因素的改變。 為了解決這一矛盾， 人們提出了模型轉(zhuǎn)移。 模型轉(zhuǎn)移， 是指在不重新建模的情況下， 通過(guò)一定算法校正新光譜的偏移， 進(jìn)而使得校正后的光譜能被原有的模型準(zhǔn)確預(yù)測(cè)。

在模型轉(zhuǎn)移中， 主光譜(A)為用于建模的那組光譜， 其在模型轉(zhuǎn)移中起主導(dǎo)作用。 而發(fā)生偏移， 需要通過(guò)模型轉(zhuǎn)移算法將其校正成類(lèi)似于主光譜的光譜被稱(chēng)為從光譜(B)[7-9]。 在模型轉(zhuǎn)移過(guò)程中， 光譜的變量數(shù)往往遠(yuǎn)大于樣本數(shù)。 這些過(guò)多的變量會(huì)增加計(jì)算的負(fù)擔(dān)， 降低預(yù)測(cè)精度。 故必須要對(duì)模型轉(zhuǎn)移中的光譜做變量選擇。 以往模型轉(zhuǎn)移使用的變量選擇算法， 大多是對(duì)主光譜進(jìn)行變量選擇， 然后從主光譜和從光譜中選擇相同的波段實(shí)現(xiàn)模型轉(zhuǎn)移。 這種方法只考慮了主光譜的有信息區(qū)段而未考慮從光譜的區(qū)段。 在實(shí)際應(yīng)用中， 由于主光譜和從光譜的差異性， 主光譜的有信息區(qū)段并非從光譜的有信息區(qū)段。 此外， 有時(shí)候主光譜和從光譜并非具有相同的波段(例如主光譜為1 100～2 500 nm， 從光譜為800～1 100 nm)， 甚至主從光譜并非同一種類(lèi)型的的光譜(例如主光譜為NIR區(qū)段， 從光譜為可見(jiàn)光譜區(qū)段)。 此時(shí)， 我們無(wú)法從主光譜和從光譜中選擇相同的波段。 為此提出采用向后迭代區(qū)間選擇法(iterative interval backward selection, IIBS)， 基于主光譜和從光譜的重要性信息， 對(duì)主光譜和從主光譜同時(shí)進(jìn)行變量選擇， 進(jìn)而獲得建模能力較強(qiáng)的波段。

1 算法原理

1.1 直接校正算法

基于光譜校正的模型轉(zhuǎn)移主要是通過(guò)建立主光譜與從光譜之間的一個(gè)轉(zhuǎn)移矩陣T來(lái)實(shí)現(xiàn)模型轉(zhuǎn)移。 主要的操作步驟是： 在主光譜和從光譜中分別找到一組濃度相同的樣本(轉(zhuǎn)移集)， 設(shè)為At (m×n1)和Bt(m×n2)， 然后通過(guò)矩陣運(yùn)算， 獲得T。 獲得T之后， 將要預(yù)測(cè)的從光譜數(shù)據(jù)乘以T， 這樣就可以得到一個(gè)類(lèi)似主光譜的光譜。 這樣， 校正后的從光譜就可以通過(guò)由主光譜建立的模型來(lái)預(yù)測(cè)。

通常， 矩陣T通過(guò)直接校正法(direct satandardization， DS)[10-13]實(shí)現(xiàn)， 故本工作就用DS算法實(shí)現(xiàn)模型轉(zhuǎn)移， 具體的算法是：

①直接用At對(duì)Bt進(jìn)行多元線(xiàn)性回歸， 進(jìn)而獲得轉(zhuǎn)移矩陣T

(1)

②對(duì)于Bp， 可以按照式(2)直接地校正成Bnew

Bnew=Bp×T

(2)

此時(shí)，Bnew就可以直接地用主光譜的模型預(yù)測(cè)。 DS算法的優(yōu)勢(shì)是， 其用光譜矩陣的整體信息進(jìn)行模型轉(zhuǎn)移， 計(jì)算較為方便。 同時(shí)， 其矩陣乘法可以用于校正變量數(shù)不同的兩個(gè)光譜矩陣。

IIBS算法的主要步驟如下：

①構(gòu)造主光譜和從光譜的變量重要性信息向量：

主光譜和從光譜均可構(gòu)造重要性信息向量， 在此用β， Res以及VIP數(shù)值分別構(gòu)造有信息向量[14-15]。 從光譜可以通過(guò)對(duì)其轉(zhuǎn)移集的PLS建模獲得相應(yīng)的變量重要性指標(biāo)， 諸如β， Res， VIP數(shù)值等。 其有信息向量簡(jiǎn)介如下：

β為回歸系數(shù)向量， 主光譜和從光譜的回歸系數(shù)可以通過(guò)PLS擬合獲得

(3)

(4)

如式(3)和式(4)所示。 在PLS模型中， 其β的絕對(duì)值大小可以作為變量選擇的指標(biāo)。 如果β的絕對(duì)值較大， 其建模能力較強(qiáng)， 故這些變量需要被選取。 因此， 通過(guò)比較主光譜和從光譜β絕對(duì)值的大小， 選擇絕對(duì)值較大的變量， 建立模型， 即可實(shí)現(xiàn)變量選擇， 降低預(yù)測(cè)誤差。

(5)

(6)

(7)

(8)

其中ej表示E矩陣中第j個(gè)列向量。 在式(7)和式(8)中可以看出， 如果第j列的殘差平方和越小，qj值就越大， 該變量有信息成分占的比例就越大， 因此該變量也就越重要。 可以選擇一些qj值較大的變量， 然后將這些變量組成一個(gè)集合， 這樣就可以提高模型的準(zhǔn)確度。 這樣， 各個(gè)變量的的q值便構(gòu)成了一條殘差向量(Res)。

VIP為變量重要性投影， 它也是通過(guò)PLS成分計(jì)算得到的向量， 其長(zhǎng)度表示的就是變量數(shù)。 一般通過(guò)設(shè)定一個(gè)閾值， 然后VIP大于這個(gè)閾值的變量就可以視為重要的變量， 進(jìn)而被選擇并建立模型。 也可以將變量按照VIP值大小排序， 選擇具有較大的VIP值的變量并組成集合， 進(jìn)而建模， 以便提高模型的精度。

②構(gòu)造光譜區(qū)間的重要性向量：

考慮到主光譜和從光譜變量數(shù)均較多， 如果直接模仿主光譜校正集的基于單個(gè)變量的變量選擇算法， 其主光譜和從光譜的變量子集合的組合將會(huì)非常多。 所以采用變量區(qū)間代替單個(gè)變量選擇的方法來(lái)提高運(yùn)算速度。 此外， 相對(duì)于離散的光譜數(shù)據(jù)點(diǎn)， 光譜波段更能反映光譜的化學(xué)信息[16]。 因此， 我們將整個(gè)光譜集就分成多個(gè)區(qū)間， 每個(gè)區(qū)間的重要性以該區(qū)間每個(gè)變量的重要性的平均數(shù)來(lái)表征。

考慮到每個(gè)變量的重要性指標(biāo)均大于零， 而且有時(shí)候區(qū)間中某個(gè)重要性較大變量可能會(huì)掩蓋重要性較小的變量。 故我們選擇幾何平均數(shù)而非算術(shù)平均數(shù)， 因?yàn)閹缀纹骄鶖?shù)既可以總體反映變量區(qū)間中各個(gè)變量的重要性信息， 也可以保證變量重要性信息受到異常的大值的影響較小。

(3)選擇重要性較大的光譜區(qū)間：

按照重要性順序排列， 將主光譜和從光譜的區(qū)間， 按照其區(qū)間的重要性排序， 選擇重要性較大的區(qū)間。 考慮到變量區(qū)間的重要性指標(biāo)會(huì)隨著變量區(qū)間數(shù)的縮小發(fā)生細(xì)微的變化， 故我們計(jì)劃用逐步刪除的辦法， 每一次迭代， 刪除一個(gè)重要性最差的區(qū)間， 最后將重要性數(shù)值較大的區(qū)間保留下來(lái)， 同時(shí)重新計(jì)算每個(gè)區(qū)段的重要性并進(jìn)行新的一輪迭代， 直到剩下最后一個(gè)區(qū)間。

考慮到光譜信息的復(fù)雜性， 有的化學(xué)信息往往在多個(gè)波段中均有體現(xiàn)。 故在區(qū)間優(yōu)化的過(guò)程中， 如果單純優(yōu)選出變量重要性最大的區(qū)間， 其建模能力也可能不是最優(yōu)， 因?yàn)橛袝r(shí)若干個(gè)建模能力較弱的區(qū)間， 其信息具有互補(bǔ)性， 其組合建模的預(yù)測(cè)效果好。 故為了提高區(qū)間選擇的效果， 用驗(yàn)證均方根誤差(root mean squared error of validation, RMSEV)來(lái)評(píng)價(jià)區(qū)間組合的建模能力。

IIBS算法的詳細(xì)流程如圖1所示。

圖1 IIBS算法的流程圖

在圖1中， IIBS先以區(qū)間重要性指標(biāo)， 通過(guò)向后篩選法構(gòu)造主光譜的一系列區(qū)間子集合。 然后對(duì)主光譜的每個(gè)集合， 計(jì)算從光譜的區(qū)間重要性指標(biāo)， 以向后篩選的方法建立從光譜的一系列區(qū)間子集合。 最后比較這些主從光譜子集合組合后模型轉(zhuǎn)移的RMSEV值， 選擇RMSEV最小的主光譜和從光譜的子集合組合。

1.2 數(shù)據(jù)集

數(shù)據(jù)集被分為四部分， 轉(zhuǎn)移集， 校正集， 驗(yàn)證集， 獨(dú)立測(cè)試集。 轉(zhuǎn)移集用于模型轉(zhuǎn)移； 校正集用于建立模型； 驗(yàn)證集用于計(jì)算驗(yàn)證誤差， 優(yōu)化參數(shù)， 進(jìn)而獲得最佳的變量集； 獨(dú)立測(cè)試集不參與模型優(yōu)化， 只用于檢驗(yàn)變量?jī)?yōu)選后模型的預(yù)測(cè)最終結(jié)果。 主光譜的轉(zhuǎn)移集用Kennard-Stone方法從主光譜的校正集中選出， 然后從光譜中和主光譜相同濃度的樣品作為轉(zhuǎn)移集的從光譜。

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

玉米數(shù)據(jù)下載于： http://www.eigenvector.com/data/Corn/index.html。 這套光譜里有三組數(shù)據(jù)集m5， mp5， mp6； 波長(zhǎng)范圍均是1 100～2 498 nm (700個(gè)波長(zhǎng)點(diǎn))。 選擇mp6作為主光譜， m5作為從光譜， 取水分?jǐn)?shù)據(jù)作為y值。 將y濃度從小到大排序， 每4個(gè)連續(xù)的樣本中取出第一個(gè)樣本， 這樣20個(gè)樣本就被取出， 剩下60個(gè)樣本為校正集。 取出的20個(gè)樣本中， 按照濃度排序， 每?jī)蓚€(gè)樣本中第一個(gè)為驗(yàn)證集， 第二個(gè)為獨(dú)立測(cè)試集。 因此驗(yàn)證集與獨(dú)立測(cè)試集的樣本數(shù)均為10。

1.2.2 小麥數(shù)據(jù)集

小麥的數(shù)據(jù)取自： http://www.wiley.com/legacy/wileychi/chemometrics/datasets.html， 這套數(shù)據(jù)有775個(gè)樣本， 1 050個(gè)波數(shù)點(diǎn)， 波長(zhǎng)范圍是400～2 498 nm。 其中蛋白質(zhì)含量作為待測(cè)指標(biāo)。 為了研究IIBS算法處理不同波數(shù)點(diǎn)數(shù)據(jù)集的能力， 我們將該數(shù)據(jù)集分為兩個(gè)部分： 可見(jiàn)-短波NIR和長(zhǎng)波NIR。 其中可見(jiàn)-短波NIR的數(shù)據(jù)點(diǎn)包括350個(gè)波長(zhǎng)點(diǎn)(400～1 098 nm)， 長(zhǎng)波NIR包括700個(gè)波長(zhǎng)點(diǎn)(1 100～2 498 nm)。 主光譜選擇長(zhǎng)波NIR， 從光譜選擇可見(jiàn)-短波NIR。 其中400條光譜作為校正集， 50條光譜作為驗(yàn)證集， 325條光譜作為獨(dú)立測(cè)試集合。

2 結(jié)果與討論

2.1 玉米數(shù)據(jù)的結(jié)果

IIBS算法中， 兩個(gè)參數(shù)對(duì)建模非常重要： 轉(zhuǎn)移集的樣本數(shù)(m)以及區(qū)間的長(zhǎng)度(n)。 選擇不同的m和n組合， 可以獲得不同組合下RMSEV的值， 同時(shí)計(jì)算不同m值下， 全光譜的RMSEV值。 以β做為變量重要性指標(biāo)， 經(jīng)過(guò)搜索， 發(fā)現(xiàn)在m=30,n=14時(shí)， 所選擇的變量可以取得較小的RMSEV值。 故選擇m=30,n=14。

在通過(guò)驗(yàn)證集確定參數(shù)后， 需要用獨(dú)立測(cè)試集檢測(cè)相應(yīng)參數(shù)下， 模型的計(jì)算結(jié)果， 結(jié)果如表1所示。

表1 玉米數(shù)據(jù)不同重要性指標(biāo)變量選擇結(jié)果

在表1中， 與全光譜變量選擇的結(jié)果相比， 基于β的IIBS算法不僅可以使得驗(yàn)證集獲得較小的RMSEV值， 而且可以使獨(dú)立測(cè)試集獲得較小的RMSEP值。 和β一樣， VIP值也可以選擇合適的變量獲得較低的RMSEV和RMSEP值。 Res雖然也可以降低RMSEV數(shù)值， 但是Res的RMSEP數(shù)值反而大于全光譜的RMSEP值。 其原因可能是Res雖然選擇了較少的變量， 而這些變量只利于校正集， 導(dǎo)致了過(guò)擬合， 反而增大了獨(dú)立測(cè)試集的RMSEP值。 故β， Res以及VIP可以被看作模型轉(zhuǎn)移中的變量重要性向量， 用于變量選擇。

為了研究變量選擇結(jié)果的化學(xué)意義， 基于β， Res， VIP的IIBS算法選擇的變量如圖2所示。

圖2 主光譜(a， c)與從光譜(b， d)光譜圖以及IIBS不同重要性向量(β， Res和VIP)的變量選擇的結(jié)果

從圖2中可以看出， 雖然β， Res， VIP三種指標(biāo)選擇的變量互不相同， 但是β和VIP選擇的變量位置相似性較高， 諸如二者的主光譜均選擇了1 450和2 300 nm附近的吸收峰， 且二者的從光譜也選擇了1 450， 1 950以及2 300 nm附近的吸收峰。 這些吸收峰都與水的吸收密切相關(guān)。 1 950 nm附近的吸收可以稱(chēng)為水的吸收Ⅰ區(qū)1 450 nm附近的吸收可以稱(chēng)為水的吸收Ⅱ區(qū)， 均與O—H的伸縮振動(dòng)有關(guān)[17-18]。 此外， 2 300 nm附近的吸收也與水的吸收有關(guān)[18-19]。 而基于Res的變量選擇算法則與β以及VIP有較大的區(qū)別， 首先， Res從主光譜中選擇了280個(gè)變量， 從從光譜選擇了462個(gè)變量。 其次， 它選擇了一些與水相關(guān)性較小的區(qū)段， 諸如1 150～1 350 nm。 此外， 它沒(méi)有選中一些與水相關(guān)性較強(qiáng)的區(qū)段， 例如： 主光譜沒(méi)有選擇2 300 nm附近的吸收峰， 從光譜沒(méi)有選擇1 950 nm附近的吸收峰。 這可能是導(dǎo)致Res選擇的變量具有較高的誤差的原因。

為了更加深入研究變量選擇的結(jié)果， 我們將數(shù)據(jù)進(jìn)行隨機(jī)分類(lèi)， 隨機(jī)生成校正集(60個(gè)樣本)， 驗(yàn)證集(20個(gè)樣本)， 獨(dú)立測(cè)試集(20個(gè)樣本)。 然后用β(n=14)， Res(n=14)以及VIP(n=20)進(jìn)行變量選擇， 利用驗(yàn)證集篩選出好的變量， 然后將其代入獨(dú)立測(cè)試集中或的預(yù)測(cè)誤差。 重復(fù)上述步驟100次， 獲得的誤差均值如圖3所示。

從圖3中可以看出， 基于β以及VIP的IIBS選擇的變量， 其RMSEP均值明顯地小于全波段的RMSEP均值， 這證明了上述算法的有效性。 選擇Res的IIBS， 其計(jì)算結(jié)果的RMSEP均值和全光譜的均值相近， 甚至在一些m值下， 其誤差反而大于全波段的RMSEP均值。 故對(duì)于玉米數(shù)據(jù)， IIBS結(jié)合β以及VIP可以選擇出較好的變量， 并獲得較低的RMSEP值。

圖3 β(a)， Res(b)以及VIP(c)不同m值條件下的Monte Carlo抽樣下的玉米數(shù)據(jù)計(jì)算結(jié)果

2.2 小麥數(shù)據(jù)的計(jì)算結(jié)果

與前者類(lèi)似， 小麥數(shù)據(jù)也被隨機(jī)分成三部分： 校正集400條光譜， 驗(yàn)證集50條光譜， 獨(dú)立測(cè)試集325條光譜。 在IIBS算法中，β， Res以及VIP算法的n值均為20。 將上述方法重復(fù)運(yùn)行100次， 其計(jì)算結(jié)果如圖4所示。

從圖4中可以看出， 與全波段建模比較， 基于β， Res， VIP的IIBS算法均可以降低獨(dú)立測(cè)試集的RMSEP值。 這證明了變量選擇的有效性。 同時(shí)， 在圖4中可以得出基于β以及Res的計(jì)算誤差要顯著小于基于VIP的計(jì)算誤差。 故對(duì)于小麥數(shù)據(jù)，β， Res可以獲得較好的變量集合。

3 結(jié) 論

向后迭代區(qū)間選擇法(iterative interval backward selection, IIBS)通過(guò)多次迭代， 每次迭代刪去重要性最小的區(qū)間， 最終獲得主光譜和從光譜模型轉(zhuǎn)移誤差最小的區(qū)間。 玉米、 小麥NIR數(shù)據(jù)測(cè)試了IIBS算法。 結(jié)果顯示， 相對(duì)于全波段， IIBS算法可以有效地從主光譜以及從光譜中同時(shí)篩選出各自有意義的波段， 實(shí)現(xiàn)降低誤差， 提高預(yù)測(cè)精度。 同時(shí)， 在選擇不同的變量重要性向量方面， 基于回歸系數(shù)的IIBS算法可以獲得較小的預(yù)測(cè)誤差。 因此， IIBS可以用于模型轉(zhuǎn)移中的變量選擇， 進(jìn)而獲得較小的誤差。