向后迭代區(qū)間選擇算法及其在近紅外光譜模型轉移中的應用

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

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

引 言

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

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

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

1 算法原理

1.1 直接校正算法

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

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

①直接用At對Bt進行多元線性回歸， 進而獲得轉移矩陣T

(1)

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

Bnew=Bp×T

(2)

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

IIBS算法的主要步驟如下：

①構造主光譜和從光譜的變量重要性信息向量：

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

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

(3)

(4)

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

(5)

(6)

(7)

(8)

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

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

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

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

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

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

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

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

IIBS算法的詳細流程如圖1所示。

圖1 IIBS算法的流程圖

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

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

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

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

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

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

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

2 結果與討論

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

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

在通過驗證集確定參數(shù)后， 需要用獨立測試集檢測相應參數(shù)下， 模型的計算結果， 結果如表1所示。

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

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

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

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

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

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

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

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

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

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

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

3 結 論

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