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

成云玲，楊蜀秦，*

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

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

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

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

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

1 材料與方法

1.1 樣本采集

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

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

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

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

1.2 儀器設備

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

1.3 預測模型

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

1.4 評價指標

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

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

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

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

2.1 模型集群分析基本思想

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

2.2 基本原理

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

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

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

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

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

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

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

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

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

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

圖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 結果與分析

3.1 不同波長選擇方法的預測結果

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

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

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

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

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

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

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

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

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

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

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

3.2 MCF性能影響因素

3.2.1 采樣次數(shù)

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

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

3.2.2 MCF結合不同回歸方法

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

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

4 結 論

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

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

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

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

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

致謝感謝西北農林科技大學葡萄酒學院劉旭副教授及其團隊在樣本采集和葡萄籽總酚含量測量中的貢獻。