基于蒙特卡羅頻率法的葡萄籽總酚含量高光譜測量變量選擇

成云玲，楊蜀秦，*

（1.西北農(nóng)林科技大學(xué) 機(jī)械與電子工程學(xué)院，咸陽712100； 2.農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)物聯(lián)網(wǎng)重點(diǎn)實(shí)驗(yàn)室，咸陽712100）

近紅外高光譜技術(shù)是利用物質(zhì)在近紅外光譜區(qū)特定的吸收特性，對樣品中一種或多種化學(xué)成分進(jìn)行快速檢測的方法［1-2］。由于其成本低、速度快、無損檢測等優(yōu)點(diǎn)，已被廣泛應(yīng)用于食品領(lǐng)域［3-5］。近紅外高光譜數(shù)據(jù)具有譜帶寬、信號弱和重疊嚴(yán)重的特點(diǎn)，一般由幾百到幾千個波段組成，相鄰波段之間共線性嚴(yán)重，并且包含有大量的冗余信息［6-7］。因此，特征波長選擇對于簡化模型，提高模型的預(yù)測精度和魯棒性具有重要意義［8-9］，這也使得波長選擇成為近紅外分析領(lǐng)域的一個熱點(diǎn)研究課題。

目前常用的變量選擇方法有連續(xù)投影算法（Successive Projections Algorithm，SPA）［10］、無信息變量消除（Uninformative Variable Elimination，UVE）法［11-12］等。SPA是一種基于變量投影比較的特征波長選擇方法。其通過比較某波長在其他波長上的投影，選擇投影向量最大的波長作為待選波長，然后基于校正模型的均方根誤差（Root Mean Square Error，RMSE）從待選波長集合中選擇最終的特征波長。UVE在保留特征波長的同時(shí)消除無信息變量，是一種基于偏最小二乘（Partial Least Squares，PLS）模型回歸系數(shù)的波長選擇方法。該方法引入穩(wěn)定性來評價(jià)模型中各變量的可靠性，從而確定最終選擇的變量，在光譜變量的選擇中得到了廣泛的應(yīng)用。這些方法的一個共同特點(diǎn)是，它們試圖為給定的數(shù)據(jù)集選擇一個固定的變量子集，而不考慮樣本變化對變量選擇的影響。

結(jié)合蒙特卡羅采樣（Monte Carlo Sampling，MCS）技術(shù)建立變量選擇方法可以有效地解決這一問題。例如，競爭自適應(yīng)重采樣（Competitive Adaptive Reweighted Sampling，CARS）［13］基于MCS和PLS回歸系數(shù)選擇特征變量。其首先通過MCS建立PLS模型，然后通過自適應(yīng)加權(quán)采樣保留模型中回歸系數(shù)絕對值權(quán)重較大的波長作為新的子集，基于新的波長子集重新建立PLS模型，經(jīng)過多次計(jì)算，選擇交叉驗(yàn)證均方根誤差（Root Mean Square Error of Cross Validation，RMSECV）最小的波長子集作為特征波長，降維性能較好。此外，其他變量選擇方法包括蒙特卡羅無信息變量消除（Monte Carlo-Uninformative Variable Elimination，MC-UVE）［14］、模型種群分析（Model Population Analysis，MPA）［15］和變量互補(bǔ)網(wǎng)絡(luò)（Variable Complementary Network，VCN）［16］等方法表明，結(jié)合MCS進(jìn)行變量選擇可以得到更好的預(yù)測結(jié)果。

本文搭建了葡萄籽總酚含量近紅外高光譜預(yù)測系統(tǒng)，根據(jù)模型集群分析的思想，提出了將MCS和波長出現(xiàn)頻次結(jié)合選擇特征波長的方法，簡稱蒙特卡羅頻率法（Monte Carlo Frequency，MCF）。該方法能夠減少建模過程中的無信息變量及干擾變量，為開發(fā)葡萄籽總酚含量檢測設(shè)備提供理論依據(jù)。

1 材料與方法

1.1 樣本采集

試驗(yàn)樣品來自陜西省楊凌盛唐酒莊，包括霞多麗、貴人香、8802、8803和雷司令5個白葡萄品種。采摘工作從葡萄轉(zhuǎn)色期至成熟期進(jìn)行，由于不同品種釀酒葡萄的成熟時(shí)間存在差異，采摘時(shí)間為2015年7月中旬至9月中旬。每個品種從轉(zhuǎn)色期一周后開始，15天作為一個采摘周期，共采摘3次。每次采摘4組葡萄，每組20個葡萄（籽）作為一個樣本。因此，每個品種包括12個樣本，總計(jì)60個樣本。將采摘的樣本快速運(yùn)送到實(shí)驗(yàn)室，手工將葡萄籽分離出來，拍攝其近紅外高光譜圖像，隨后用蛋白質(zhì)沉淀法［17］測量其總酚含量。

原始圖像具有256個波長，去除兩端包含噪音的波段，最終選擇950～1600 nm之間196個波長對應(yīng)的光譜數(shù)據(jù)。為了去除高頻隨機(jī)噪聲等干擾因素［18］，采用S-G濾波［19］對數(shù)據(jù)進(jìn)行預(yù)處理，預(yù)處理后的平均光譜如圖1所示。隨機(jī)將樣本按4∶1的比例分為訓(xùn)練集和測試集，葡萄籽總酚含量分布如表1所示。

圖1 預(yù)處理后的5個品種葡萄籽平均光譜Fig.1 Average spectra of five types of pretreated grape seeds

表1 葡萄籽總酚含量分布統(tǒng)計(jì)Table 1 Distribution statistics of total phenol content in grape seeds

1.2 儀器設(shè)備

本文采用HyperSIS高光譜成像系統(tǒng)，包括IMSpector N17E型近紅外成像光譜儀（Spectral Imaging Ltd.，F(xiàn)inland）、320像素×256像素的XEVA3616型面陣CCD相機(jī)（XenICs Ltd.，Belgium）、白光漫反射型鹵鎢白熾燈（4個）、高精度電動平移載物臺和一臺計(jì)算機(jī)。該系統(tǒng)的光譜分辨率為2.8 nm，采集的波長范圍為900～1700 nm，平臺移動速度為20 mm/s，相機(jī)曝光時(shí)間為10 ms。本文方法均在MATLAB R2017a中實(shí)現(xiàn)。

1.3 預(yù)測模型

本文采用支持向量回歸（Support Vector Regression，SVR）建立葡萄籽總酚的預(yù)測模型，用于比較MCF和其他常用變量選擇方法的性能。SVR是一種適用于解決小樣本、非線性及高維數(shù)據(jù)問題的方法［20-21］。SVR通過核函數(shù)，將低維空間向量映射到高維空間，在高維空間中構(gòu)造線性決策函數(shù)來實(shí)現(xiàn)原空間中的非線性決策。考慮到高光譜數(shù)據(jù)和預(yù)測變量總酚含量之間映射關(guān)系的復(fù)雜性和非線性，利用SVR建模在一定程度上規(guī)避了過擬合風(fēng)險(xiǎn)。本文使用Libsvm［22］實(shí)現(xiàn)SVR方法。

1.4 評價(jià)指標(biāo)

本文采用相關(guān)系數(shù)R2和RMSE作為回歸模型的評價(jià)指標(biāo)。R2越高，RMSE越低，表明模型的效果越好。R2和RMSE的計(jì)算公式分別為

式中：y′i和yi分別為第i個樣本的預(yù)測值和真實(shí)值為樣本的平均值；n為樣本個數(shù)。

2 基于蒙特卡羅頻率法選擇波長變量

蒙特卡羅方法是一種基于隨機(jī)數(shù)和概率統(tǒng)計(jì)來研究問題的技術(shù)。本文提出的MCF是一種基于MCS和波長出現(xiàn)頻次的變量選擇方法，該方法可以和多種回歸方法結(jié)合，能夠有效選擇特征變量。

2.1 模型集群分析基本思想

模型集群分析［15］方法是首先通過MCS獲取數(shù)據(jù)子集；然后針對每個子數(shù)據(jù)集，建立一個子模型；最后從樣本空間、變量空間、參數(shù)空間或者模型空間中，對所有建立的集群子模型的參數(shù)進(jìn)行統(tǒng)計(jì)分析，以獲得有用信息。

2.2 基本原理

MCF選擇特征波段主要采用MCS選擇波長子集，然后利用波長子集建立大量回歸子模型；選擇RMSE較小的子模型，統(tǒng)計(jì)每個波長出現(xiàn)的頻次；根據(jù)指數(shù)遞減函數(shù)選擇波長個數(shù)，選取頻次最高的波長作為特征波長，具體步驟如下：

1）MCS選擇波長子集。設(shè)樣本的光譜矩陣X為n×q，表示矩陣由n個樣本和q個波長組成；化學(xué)值Y由向量n×1表示。根據(jù)模型集群分析的建模思想，首先采用MCS對所有波長進(jìn)行采樣，每次隨機(jī)選擇p（p＜q）個波長，可得到n×p的子光譜矩陣。將該過程重復(fù)N次（N＞1 000），得到N個子數(shù)據(jù)集（Xsub，Y）i，i=1，2，…，N。此過程不僅能得到N組不同變量的組合，還能確保每個變量具有相同的采樣頻率。

2）將子數(shù)據(jù)集按4∶1隨機(jī)分為訓(xùn)練集和預(yù)測集，建立N個回歸子模型。SVR是一種適用于解決小樣本及高維數(shù)據(jù)問題的最常用的建模方法，因此本文采用SVR建立預(yù)測模型，然后計(jì)算N個子模型預(yù)測集的RMSE。

3）計(jì)算波長出現(xiàn)的頻率。將上述所有子模型按照預(yù)測RMSE從小到大進(jìn)行排序，只保留預(yù)測結(jié)果較好的前K個子模型，計(jì)算這些模型中各個波段出現(xiàn)的頻次f。一般波段出現(xiàn)頻次越高，則認(rèn)為該波段和化學(xué)值相關(guān)性越高，根據(jù)頻次對波長進(jìn)行重要性排序。f的計(jì)算公式為

式中：i為波段序號；j表示保留的子模型；Fi表示波段i是否出現(xiàn)在模型j中，若出現(xiàn)則為1，否則為0；K為保留子模型的個數(shù)。

4）根據(jù)指數(shù)遞減函數(shù)選擇波長個數(shù)。建立m個SVR回歸模型，根據(jù)模型的預(yù)測RMSE選擇最佳的特征波長個數(shù)。指數(shù)遞減函數(shù)［13］定義為

式中：ri為第i次選擇的波長個數(shù)；?和k是由以下2個條件決定常數(shù)：

1）在第一次運(yùn)行中，所有q個波長都被用來建模。由于本文共采用了波段裁剪后的196個原始波長，因此r1=196。

2）在第m（本文取m=40）次運(yùn)行中，只保留2個波長，即r40=2。

在這2個條件下，?和k的計(jì)算公式分別為

圖2所示為指數(shù)遞減函數(shù)選擇波長個數(shù)的過程。可以看出，波長個數(shù)呈遞減趨勢，并分為2個階段，第1階段波長數(shù)量下降較快，可快速去除出現(xiàn)頻次少的波長；第2階段波長數(shù)量下降緩慢，可有效保留出現(xiàn)頻次較高的波段。

圖2 MCF特征波長個數(shù)選擇Fig.2 Selection of number of characteristic wavelengths by MCF

3 結(jié)果與分析

3.1 不同波長選擇方法的預(yù)測結(jié)果

為了驗(yàn)證提出方法的有效性，將MCF與SPA、CARS等2種方法分別結(jié)合SVR方法，預(yù)測葡萄籽中的總酚含量。

MCF通過MCS，每次從訓(xùn)練集中隨機(jī)選擇100個波段構(gòu)建回歸子模型。重復(fù)2 000次，并計(jì)算各子模型的RMSE。然后，對所有子模型進(jìn)行排序，根據(jù)RMSE值由小到大，分別選擇前10%、20%、30%、40%和50%的子模型，計(jì)算模型中每個波長出現(xiàn)的頻率。根據(jù)式（4）的指數(shù)遞減函數(shù)選擇RMSE較小的波長數(shù)作為特征波長個數(shù)（見圖2），并找出出現(xiàn)頻率最高的波長子集作為特征波長。實(shí)驗(yàn)表明，前30%的子模型預(yù)測性能最佳，因此，本文保留前30%子模型用于計(jì)算波長頻率。各波段頻次分布如圖3所示，最終MCF選擇9個特征波長。

在SPA降維算法中，設(shè)置波長個數(shù)范圍為2～49，根據(jù)訓(xùn)練集RMSE值確定最佳的光譜變量總數(shù)。當(dāng)波長數(shù)量較少時(shí)，RMSE值較大，隨著波長個數(shù)的增加，RMSE開始呈下降趨勢，當(dāng)選取18個波長時(shí)達(dá)到最小值。因此，SPA最終選擇的波長個數(shù)為18。

圖3 波段頻次分布Fig.3 Frequency distribution of spectral bands

圖4 CARS方法波長的系數(shù)變化Fig.4 Coefficient variation of wavelength by CARS

圖4為采用CARS進(jìn)行特征波長選擇后，各波長的回歸系數(shù)路徑，設(shè)置采樣次數(shù)為50次。每條線反映一個波長系數(shù)的變化。星號線處的臨界點(diǎn)表示RMSECV的最優(yōu)子集，星號之后由于有效波長的去除，RMSECV值開始增大。根據(jù)RMSECV值最小的原則，CARS共選擇了7個特征波長。

圖5為3種方法選擇的變量分布，直觀地給出了方法所選變量的波長分布。可以看出，3種方法波長選擇區(qū)間大致相同，主要集中在950～1 400 nm。由于高光譜圖像光譜分辨率高，光譜曲線幾乎連續(xù)分布，相鄰波長之間數(shù)據(jù)相關(guān)性強(qiáng)，而MCF的特征波長分布較為均勻，說明該方法在去除冗余信息方面具有優(yōu)勢。此外，SPA 選取18個特征波長，波長個數(shù)最多。CARS選擇波長個數(shù)最少，其光譜包含的信息量少，因此可能導(dǎo)致模型效果不理想。可進(jìn)一步根據(jù)回歸模型的R2和RMSE比較3個方法的優(yōu)劣。

以總酚含量為因變量（Y），基于波長選擇的光譜為自變量（矩陣X）構(gòu)建SVR模型。采用高斯徑向基函數(shù)（Radial Basis Function，RBF）作為SVR的核函數(shù)，通過網(wǎng)格尋優(yōu)算法找到使分類模型最佳的懲罰函數(shù)c和核函數(shù)參數(shù)g，c和g尋優(yōu)范圍取［2-8，216］。通過訓(xùn)練集和預(yù)測集的R2和RMSE對模型性能進(jìn)行評價(jià)。

圖5 3種方法選擇的變量分布Fig.5 Distribution of variables selected by three methods

表2 不同降維方法的總酚預(yù)測結(jié)果比較Table 2 Comparison of total phenol pr ediction results with different dimensionality reduction methods

不同降維方法對葡萄籽總酚含量預(yù)測結(jié)果如表2所示。可以看出，全光譜模型的預(yù)測R2和RMSE分別約為0.90和0.42，表明總酚含量與高光譜數(shù)據(jù)高度相關(guān)。對比3種方法，MCF降維后模型具有最大的預(yù)測R2（0.91）和最小的RMSE（0.37），預(yù)測結(jié)果最好。該方法分別選擇了958、1 044、1 091、1 127、1 230、1 264、1 280、1 317和1323 nm處的特征波長。SPA選取波段個數(shù)最多，其模型結(jié)果略低于MCF，預(yù)測R2達(dá)到0.89。CARS選取波長個數(shù)最少，同時(shí)預(yù)測效果最差，預(yù)測相關(guān)系數(shù)均小于0.80。此外，MCF降維后的波長預(yù)測結(jié)果優(yōu)于全波段，說明該波長選擇方法可以提高模型的預(yù)測準(zhǔn)確度。

3.2 MCF性能影響因素

3.2.1 采樣次數(shù)

為了研究MCS次數(shù)對MCF性能的影響，將波長采樣次數(shù)分別設(shè)置為1 000、2 000、3 000、4 000和5 000次，建立SVR子模型并統(tǒng)計(jì)各模型的預(yù)測RMSE值，箱型圖如圖6所示。由圖可得，不同采樣次數(shù)下N個模型的RMSE最大值、最小值和中值接近，RMSE分布無明顯差別。結(jié)果表明，MCS次數(shù)對MCF的性能沒有顯著影響。因此，本文采用2 000次作為默認(rèn)波長采樣次數(shù)。

圖6 MCF不同采樣次數(shù)的箱型圖Fig.6 Box graph of MCF with different sampling times

3.2.2 MCF結(jié)合不同回歸方法

對MCF結(jié)合不同回歸方法的性能進(jìn)行比較。除SVR之外，還采用最小二乘回歸（Partial Least Squares Regression，PLSR）法、RBF神經(jīng)網(wǎng)絡(luò)建立子模型選擇特征波段。為了比較不同回歸方法的波段選擇效果，采用蒙特卡羅對波長采樣2 000次，分別用SVR、PLSR和RBF這3種方法建立回歸子模型，選擇出現(xiàn)頻次最高的前9個波長作為特征波長。用特征波長建立葡萄籽總酚的SVR預(yù)測模型，表3為不同模型的預(yù)測結(jié)果。MCF結(jié)合3種回歸方法進(jìn)行波段選擇，預(yù)測R2達(dá)到0.85～0.91，RMSE約為0.37～0.55。而其中采用SVR建立子模型進(jìn)行波段選擇時(shí)，預(yù)測效果最佳。

表3 MCF結(jié)合不同回歸方法的總酚預(yù)測結(jié)果比較Table 3 Comparison of total phenol prediction results of MCF combined with different regression methods

4 結(jié) 論

本文提出了一種基于MCS和波長出現(xiàn)頻次結(jié)合的變量選擇方法，簡稱蒙特卡羅頻率法（MCF）。

1）針對葡萄籽總酚近紅外高光譜，利用MCF進(jìn)行特征波長選擇，波長數(shù)目由196個減少到9個。

2）采用SVR建立總酚的回歸模型，預(yù)測R2和RMSE分別約為0.91和0.37。與其他變量選擇方法相比，MCF在減少無信息變量和干擾變量的同時(shí)提高了模型的預(yù)測結(jié)果。

3）討論了MCS次數(shù)和不同回歸方法對MCF性能的影響。結(jié)果表明，采樣次數(shù)對MCF波長選擇無顯著影響，采用SVR建立子模型進(jìn)行波段選擇時(shí)，模型效果最佳。

因此，MCF可以作為一種有效的波長選擇工具應(yīng)用于高光譜數(shù)據(jù)分析，具有良好的預(yù)測性能。

致謝感謝西北農(nóng)林科技大學(xué)葡萄酒學(xué)院劉旭副教授及其團(tuán)隊(duì)在樣本采集和葡萄籽總酚含量測量中的貢獻(xiàn)。