基于變量優(yōu)選和近紅外光譜技術(shù)的紅富士蘋果產(chǎn)地溯源

張立欣，楊翠芳，陳杰，張曉果，張楠楠，張曉*

1(塔里木大學 信息工程學院，新疆維吾爾自治區(qū) 阿拉爾，843300)2(南京理工大學 理學院，江蘇 南京，210094) 3(河南城建學院 數(shù)理學院，河南 平頂山，467036)

中醫(yī)認為蘋果具有生津止渴、潤肺除煩、健脾益胃、養(yǎng)心益氣等功效，并且它的味道酸甜適口，營養(yǎng)豐富，因此，成為老幼皆宜的水果之一。由于環(huán)境因素、土壤特征等的不同，各產(chǎn)地蘋果的品質(zhì)存在差異，其口感也不盡相同，而這些差異無法通過肉眼直接準確地辨別，化學鑒定法費時費力，不適合進行大批量操作。

近紅外光譜檢測以其方便、高效、快速、無污染、無需對樣品預處理等優(yōu)點，被廣泛應用到現(xiàn)代農(nóng)業(yè)檢測分析中。主要應用領(lǐng)域包括：農(nóng)產(chǎn)品成熟度的鑒別[1]、復溶果品損傷的鑒別[2]、復溶同一農(nóng)產(chǎn)品不同基因型的鑒別[3]、復溶農(nóng)產(chǎn)品產(chǎn)地的鑒別[4]、農(nóng)產(chǎn)品新鮮度[5]、貨架期[6]、霉變程度[7]、摻假的檢測[8]等。在光譜分析中，經(jīng)常會受到背景等隨機因素的干擾，因此需對光譜數(shù)據(jù)進行預處理[9]，常用的光譜預處理方法有：一階導數(shù)(first derivative，1-DER)[10]、二階導數(shù)(second derivative，2-DER)[11]、標準正態(tài)變換(standard normal transformation，SNV)[12]，多元散射校正(multivariate scatter correction，MSC)[13]、平滑變換(smooth transformation，SG)[14]、標準化[15]、歸一化(normalization，NOR)[16]、中心化(centralization,CEN)[17]等。為降低模型的復雜度，減少共線性的干擾，需要提取特征波段[18]，常用的方法有：連續(xù)投影算法(successive projection algorithm，SPA)[19]、競爭性自適應重加權(quán)算法(competitive adaptive reweighted sampling，CARS)[2]、主成分分析(principal component analysis，PCA)[20]、隨機蛙跳算法(random frog，RF)[21]等，也有將幾種方法聯(lián)合起來選擇特征變量的。程介虹等[22]提出一種改進聯(lián)合區(qū)間的RF選擇特征波長，通過聯(lián)合區(qū)間偏最小二乘法對全譜進行變量初選，此時得到的波長對目標變量變化最為敏感，將其作為RF的初始變量子集，以解決其運行時間較長、效率較低的問題。袁凱等[23]采用3步混合策略，提出了間隔偏最小二乘、區(qū)間變量迭代空間收縮法和迭代保留信息變量聯(lián)用的特征變量選擇方法，對生鮮雞胸肉的近紅外光譜進行特征波長選擇，建立了雞肉水分R檢測模型。結(jié)果表明，建模波長數(shù)量經(jīng)3步選擇后減少為全光譜建模的0.76%，但模型精確度和穩(wěn)定性逐步提高。FANG等[2]將SPA、CARS、過濾式特征選擇。3種方法選取的特征變量組合起來建模，取得了很好的預測效果。

研究表明，光譜技術(shù)在農(nóng)產(chǎn)品的檢測分析中具有廣泛的應用，但是，對蘋果產(chǎn)地溯源的研究相對較少。目前見刊的有馬永杰等[24]采用多種數(shù)據(jù)降維方法，建立k近鄰(k-nearest neighbor, KNN)模型，對紅富士蘋果進行產(chǎn)地溯源，建模集和預測集的正確率分別達到97.3%和92.3%，模型的正確率有待進一步提高。本研究在前人研究的基礎(chǔ)上，基于近紅外光譜技術(shù)，以新疆阿克蘇、甘肅靜寧、河南靈寶和山東煙臺的紅富士蘋果為研究對象，利用光譜分析的理論和方法，借助于光譜預處理算法，特征波長選擇方法，建立概率神經(jīng)網(wǎng)絡(luò)模型對蘋果的產(chǎn)地進行判別分析，重點研究不同的光譜預處理方法、不同的變量篩選方法對預測模型的影響，為實現(xiàn)紅富士蘋果的產(chǎn)地溯源提供理論參考。

1 材料與方法

1.1 材料與儀器

以阿克蘇、靜寧、靈寶和煙臺的紅富士蘋果為實驗對象，在實驗中所使用的蘋果均為2020年10月份在各蘋果產(chǎn)地郵寄所得。挑選表面沒有缺陷、直徑范圍為65～85 cm、大小均勻的蘋果樣品256個，去除表面的污垢，放置在冰柜內(nèi)保存，溫度控制在4 ℃，實驗前分批拿出，待其恢復到室溫(20～25 ℃)環(huán)境后開始實驗。

實驗中所用的推掃式高光譜分選系統(tǒng)(Hyperspectral Sorting System)為北京卓立漢光公司生產(chǎn)。光譜測定的范圍為900～1 700 nm(實際可測量到1 750 nm)，光譜分辨率5 nm，光譜采樣點4 nm，選取果身中心前后左右4個方位，提取大小為20像素×15像素，4個面均進行提取，共1 200像素點，選取平均值為該樣本反射率。通過自帶的ENVI5.3軟件提取ROI的光譜值，最后導出為Excel文件。

1.2 光譜數(shù)據(jù)的采集和校正

為了得到清晰的圖像，在采集光譜數(shù)據(jù)前需要多次的相機聚焦和移動平臺的速度測試，反復嘗試后，確定平臺的移動速度為0.35 cm/s，相機曝光時間為0.09 s。為了減少光照不均勻和暗電流對實驗的影響，需要對采集到的光譜數(shù)據(jù)進行黑白校正，校正公式如公式(1)所示：

R=(I-B)/(W-B)

(1)

式中：R是校正后的光譜數(shù)據(jù)，I是原始光譜數(shù)據(jù)，W為對準白板采集到的數(shù)據(jù)，B是蓋上相機鏡頭采集到的數(shù)據(jù)。

1.3 數(shù)據(jù)處理

1.3.1 光譜數(shù)據(jù)的預處理

在光譜檢測的過程中，會受到樣品背景等隨機因素的影響，導致光譜數(shù)據(jù)中含有噪聲，為提高模型的準確性和穩(wěn)健性，需要對數(shù)據(jù)進行預處理，采用的方法有NOR、CEN、1-DER、2-DER、SNV、MSC、小波變換(wave transformation，WT)、SG、傅里葉變換(Fourier transformation，F(xiàn)T)。

1.3.2 特征波長提取

光譜能夠體現(xiàn)物質(zhì)所含成分及含量，但同時包含大量的冗余信息，為降低模型的復雜性，減少共線性的影響，因此，需要提取特征波長。采用的方法有PCA、SPA、CARS和RF。

1.3.3 判別分析模型

概率神經(jīng)網(wǎng)絡(luò)(probabilistic neural network，PNN)是徑向基網(wǎng)絡(luò)的一個分支，屬于前饋網(wǎng)絡(luò)的一種。它具有學習過程簡單、訓練速度快、分類更準確、容錯性好等優(yōu)點。

1.3.4 模型驗證

將樣本以2∶1的比例間隔分為訓練集和測試集，依靠訓練集建立上述判別分析模型，測試集將通過已經(jīng)建立好的PNN模型進行驗證，以正確率為標準來評判各類方法的優(yōu)劣。

2 結(jié)果與分析

采集的紅富士蘋果光譜數(shù)據(jù)中，剔除異常值后，共得到阿克蘇、靜寧、靈寶、煙臺4個產(chǎn)地的蘋果樣本各60個，其原始光譜曲線如圖1所示。

圖1 原始光譜圖Fig.1 Original spectrogram

不同產(chǎn)地的紅富士蘋果樣本的光譜曲線變化趨勢大致相同，反映不同產(chǎn)地的蘋果之間也有著極大的相似性。但是在一些波峰和波谷處存在明顯的偏離，這是由于不同產(chǎn)地的蘋果內(nèi)部成分含量的多少存在差異。在1 060 nm處也有1個明顯的波峰，這與C—O—C基團有關(guān)；在1 440 nm處有1個明顯的波谷，這與H2O的二倍頻吸收帶有關(guān)。成分含量上的差異會導致光譜曲線的差異，這為基于光譜技術(shù)對蘋果產(chǎn)地的識別研究提供了信息。

阿克蘇、靜寧、靈寶、煙臺4個產(chǎn)地的蘋果樣本均以2∶1的比例間隔產(chǎn)生訓練集和測試集，得訓練集中4個產(chǎn)地的蘋果樣本各為40個，即160個蘋果樣本組成訓練集，測試集中4個產(chǎn)地的蘋果樣本各為20個，即80個蘋果樣本組成測試集。

2.1 光譜數(shù)據(jù)的預處理

在光譜檢測的過程中，會受到樣品背景等隨機因素的影響，導致光譜數(shù)據(jù)中含有噪聲，為提高模型的準確性和穩(wěn)健性，需要對數(shù)據(jù)進行預處理。分別采用NOR、CEN、1-DER、2-DER、SNV、MSC、WT、SG、FT等9種預處理方法對原始光譜進行預處理，以消除各類隨機因素對模型性能的干擾。分別以沒有經(jīng)過預處理的原始光譜(no pretreat, NO)和經(jīng)過預處理之后的光譜數(shù)據(jù)作為輸入自變量，分別建立線性的偏最小二乘判別分析(partial least squares discrimination analysis，PLSDA)模型和非線性的PNN模型，對蘋果的產(chǎn)地進行判別，其結(jié)果如表1和表2所示。

表1 不種預處理方法下PLSDA模型的準確率 單位：%

由表1可知，不同的預處理方式下，PLSDA模型的總準確率為60%～85%，以總準確率最高為標準，最優(yōu)預處理方法為MSC，此時阿克蘇、靜寧、靈寶、煙臺4個產(chǎn)地的蘋果識別率分別為100%、95%、85%、60%。由表2可知，PNN模型總準確率為86.25%～97.50%，最優(yōu)預處理方法亦為MSC，此時阿克蘇、靜寧、靈寶、煙臺4個產(chǎn)地的蘋果識別率分別為100%、100%、90%、100%。

表2 不種預處理方法下PNN模型的準確率 單位：%

比較表1和表2可以發(fā)現(xiàn)，在相同的預處理方式下，PNN模型的總準確率高于PLSDA模型的總準確率，這是因為蘋果內(nèi)部的結(jié)構(gòu)復雜，光在蘋果內(nèi)部的傳輸是一種復雜的結(jié)構(gòu)，擁有非線性判別能力的PNN優(yōu)于線性的PLSDA。為提高模型的識別率，后續(xù)的判別分析模型中，均采用MSC預處理之后的光譜數(shù)據(jù)，建立PNN模型。

2.2 特征波長的選取

對MSC預處理之后的光譜數(shù)據(jù)，分別采用PCA、SPA、CARS提取光譜數(shù)據(jù)中的重要變量作為建模輸入自變量。

2.2.1 PCA選取特征變量

采用PCA算法提取主成分，各主成分解釋的總方差如表3所示。

表3 解釋的總方差Table 3 Eexplained total variance

前3個主成分累計貢獻率達到95%以上，因此，選取前3個主成分，作為下一步判別分析模型的輸入自變量。不同產(chǎn)地的蘋果光譜數(shù)據(jù)前3個主成分的得分分別如圖2所示。

圖2 主成分圖Fig.2 Principal component diagram

由圖2可知，阿克蘇紅富士蘋果的前3個主成分可以和其他產(chǎn)地的分開，而其他3個產(chǎn)地的蘋果光譜數(shù)據(jù)之間重疊比較多，這將影響模型的識別效果。

2.2.2 SPA選取特征變量

SPA進行波長變量選擇，指定波長數(shù)為1～20，采用均方根誤差最小來確定最終參與建模的波長變量個數(shù)，選取過程如圖3所示。

圖3 SPA選取變量過程Fig.3 Variable selection process

由圖3可知，隨著所選的特征波長數(shù)的增加，參與建模的波長變量數(shù)增加，而均方誤差整體呈現(xiàn)遞減的趨勢。當選取特征波長變量數(shù)為14時，均方根誤差為0.318 4，之后均方誤差逐漸趨于平緩，此時，若再多選取變量，過多的波長變量參與建模會增加模型的復雜性，而均方根誤差并不會有太大的變化。因此，最終選取14個波長變量參與建模，選取的特征波長如圖4所示。

圖4 SPA選取的特征變量Fig.4 Selected variable

最終選取的14個特征波長為圖4中的小方塊對應的橫坐標，即對應波長分別為911.06、923.53、929.78、980.02、1 092.04、1 175.88、1 368.46、1 402.29、1 453.03、1 529.06、1 693.72、1 704.37、1 711.48、1 715.03 nm。

2.2.3 CARS選取特征變量

采用CARS方法進行特征波長的選取，蒙特卡羅方法抽樣迭代200次，抽取過程如圖5所示。

a-變量優(yōu)化過程；b-均方根誤差變化趨勢；c-回歸系數(shù)變化圖5 CARS 選取變量結(jié)果Fig.5 Variable selection results

從抽取結(jié)果來看，到第82次迭代時，均方根誤差達到最小為0.290 1，此時選出35個波長變量，對應的波長為：914.17、920.41、929.78、932.90、967.41、1 127.08、1 133.56、1 136.81、1 146.55、1 375.21、1 378.59、1 405.68、1 409.08、1 412.47、1 415.87、1 419.27、1 456.80、1 484.23、1 487.67、1 491.11、1 494.55、1 497.99、1 501.43、1 504.88、1 508.33、1 515.23、1 522.14、1 529.06、1 581.17、1 584.66、1 588.15、1 591.64、1 697.27、1 700.82、1 736.41 nm。

2.2.4 RF選取特征變量

采用RF算法提取特征波長變量，迭代1 000次，每個波長變量被選中的頻率如圖6所示。

圖6 每個波長變量被選中的概率圖Fig.6 Probability diagram of each wavelength variable being selected

由圖6可知，波長變量被選中的概率范圍為0～0.941 0，優(yōu)先選擇率較大的波長變量，采用交叉驗證的均方根誤差最小來確定最終選取的變量個數(shù)，均方根誤差和所選變量個數(shù)的關(guān)系如圖7所示。

圖7 RF選取變量過程Fig.7 Process of selecting variables by RF method

由圖7可知，隨著所選變量個數(shù)的增加，均方根誤差迅速下降，當所選變量個數(shù)為82時，交叉驗證的均方根誤差最小，為0.272 2。所選取的變量為圖6中概率大于0.088 0的波長，即圖6中概率在水平線上方的波長變量，所對應的波長分別為：911.06、914.17、917.23、920.41、923.53、929.78、932.90、945.43、957.98、967.41、970.56、976.86、980.02、989.49、998.98、1 002.15、1 008.48、1 014.83、1 049.87、1 065.87、1 078.70、1 101.23、1 114.14、1 130.32、1 136.81、1 143.30、1 146.55、1 149.80、1 153.06、1 156.31、1 159.57、1 166.09、1 169.35、1 225.09、1 234.98、1 238.28、1 241.59、1 244.89、1 278.04、1 308.03、1 344.89、1 348.25、1 351.61、1 365.09、1 371.83、1 375.21、1 385.35、1 398.90、1 405.68、1 412.47、1 415.87、1 419.27、1 484.23、1 491.11、1 494.55、1 497.99、1 501.43、1 504.88、1 508.33、1 511.78、1 515.23、1 522.14、1 529.06、1 532.52、1 549.85、1 570.71、1 591.64、1 605.63、1 658.34、1 665.41、1 672.47、1 676.01、1 679.55、1 683.09、1 686.63、1 690.17、1 693.72、1 697.27、1 700.82、1 711.48、1 715.03、1 736.41 nm。

由于CARS和RF算法選出的 特征變量的個數(shù)仍然比較多，變量之間可能還存在共線性，為簡化模型，將CARS和RF選出的特征變量再進一步優(yōu)選。

2.2.5 CARS-SPA選取特征變量

對CARS方法選取的特征波長變量，采用SPA算法進一步優(yōu)選，變量的選取過程如圖8所示。

圖8 CARS-SPA優(yōu)選變量過程Fig.8 Process of optimizing variables by CARS-SPA method

隨著所選變量個數(shù)的增加，均方根誤差呈現(xiàn)遞減的趨勢，當選取變量個數(shù)為16時，均方根誤差達到最小0.261 9。此時選取的特征波長變量如圖9中小方塊所對應的橫坐標。

圖9 CARS-SPA優(yōu)選的變量Fig.9 Preferred variables by CARS-SPA method

最后選出16個特征波長變量，對應的波長為：920.41、929.78、967.41、1 127.08、1 133.56、1 146.55、1 375.21、1 378.59、1 405.68、1 412.47、1 415.87、1 504.88、1 581.17、1 591.64、1 700.82、1 736.41 nm。

2.2.6 RF-SPA選取特征變量

對RF方法選取的特征波長變量，采用SPA算法進一步優(yōu)選，變量的選取過程如圖10所示。

圖10 RF-SPA優(yōu)選變量過程Fig.10 Process of optimizing variables by RF-SPA method

隨著所選變量個數(shù)的增加，均方根誤差整體下滑，綜合考慮均方根誤差和模型的復雜性，當所選變量為17時，均方根誤差為0.260 9，之后，隨著變量個數(shù)的增加，均方根誤差并無明顯減少的趨勢，因此，選取17個特征波長變量。所選取的變量波長如圖11中小方塊對應的橫坐標。

圖11 RF-SPA 優(yōu)選的變量Fig.11 Preferred variables by RF-SPA method

最后選出的17個特征波長變量對應的波長為：911.06、914.17、920.41、929.78、932.90、980.02、1 114.14、1 169.35、1 371.83、1 405.68、1 453.03、1 515.23、1 591.64、1 658.34、1 693.72、1 700.82、1 736.41 nm。

2.3 建模結(jié)果分析

以選出的特征波長變量作為輸入變量，建立PNN判模型，對蘋果的產(chǎn)地進行判別分析。模型識別的準確率如表4所示。

表4 不同的變量選取方式下模型的準確率 單位：%

由表4可知，總準確率范圍為71.25%～100%，PCA方法選取的特征變量建模識別率最低，這是因為PCA在選取過程中，只是對樣品的自變量進行重新組合，沒有考慮因變量的影響。SPA算法可以最大程度地降低被選中變量間的共線性問題，但是對于光譜而言，有效變量之間的投影距離不一定最大，因此，篩選出的變量子集中可能包含一些無用信息，甚至是干擾信息，降低模型的泛化能力。CARS和RF算法選出的特征變量個數(shù)比較多，無法克服變量之間的共線性的影響，這都將會影響模型的泛化能力。CARS-SPA和RF-SPA在初選特征變量的基礎(chǔ)上再進一步優(yōu)選，選出的變量個數(shù)適中，消除了變量之間的共線性。其中CARS-SPA選出的特征變量建模，阿克蘇、靜寧、靈寶、煙臺4個產(chǎn)地的紅富士蘋果識別率分別達到了100%、100%、95%、100%，總識別率達到了98.75%。綜合考慮識別率和模型的復雜性，最優(yōu)模型為MSC-CARS-SPA-PNN模型。

2.4 與KNN模型的比較

采用馬永杰等[24]的K最近鄰算法，分別以PCA、SPA、CARS、RF、CARS-SPA、RF-SPA選出的特征波長變量作為輸入變量，建立KNN模型。k的取值范圍為3～9，通過交叉驗證的方式確定最佳的k值，識別結(jié)果如表5所示。

表5 KNN的識別結(jié)果 單位:%

比較表4和表5可知，在PCA、SPA、CARS、CARS-SPA變量選取方式下，PNN模型的總準確率均高于KNN模型的總準確率；在RF變量選取方式下，PNN模型和KNN模型有相同的總準確率；在RF-SPA變量選取方式下，PNN模型的總準確率略低于KNN模型的總準確率。但是，在CARS-SPA變量選取方式下，PNN模型的總準確率達到最高98.75%。

3 討論與結(jié)論

采集阿克蘇、靜寧、靈寶和煙臺的紅富士蘋果近紅外光譜數(shù)據(jù)，分別進行歸一化、中心化、一階導數(shù)、二階導數(shù)、標準正態(tài)變換、多元散射校正、小波變換、平滑變換、傅里葉變換共9種預處理方法，建立概率神經(jīng)網(wǎng)絡(luò)模型對蘋果產(chǎn)地進行識別，結(jié)果表明，多元散射校正預處理的效果最優(yōu)，4個產(chǎn)地的蘋果樣品識別率分別為100%、100%、90%、100%，總識別率為97.5%。

為了簡化模型提高識別率，將多元散射校正預處理后的光譜數(shù)據(jù)分別采用主成分法、連續(xù)投影算法、競爭性自適應重加權(quán)算法、隨機蛙跳算法、競爭性自適應重加權(quán)-連續(xù)投影算法、隨機蛙跳-連續(xù)投影算法選取特征變量建模。研究結(jié)果表明，競爭性自適應重加權(quán)-連續(xù)投影算法效果最佳，MSC-CARS-SPA-PNN模型的總識別率達到了98.75%，4個產(chǎn)地的紅富士蘋果識別率分別達到了100%、100%、95%、100%。因此，近紅外光譜技術(shù)能夠快速、無損地識別蘋果的產(chǎn)地，為紅富士蘋果的產(chǎn)地溯源提供理論參考。