陳嘉偉,周德強(qiáng)*,崔晨昊,任志俊,左文娟
1. 江南大學(xué)機(jī)械工程學(xué)院,江蘇 無錫 214122 2. 江蘇省食品先進(jìn)制造裝備技術(shù)重點實驗室,江蘇 無錫 214122 3. 布勒中國創(chuàng)新中心,江蘇 無錫 214111
用小麥粉加工出的食品是人們?nèi)粘o嬍车闹匾M成部分,而不同的小麥?zhǔn)称穼ψ鳛樵系男←湻塾胁煌钠焚|(zhì)要求,其中小麥粉的粉質(zhì)特性是主要參考指標(biāo)之一。粉質(zhì)特性一般通過粉質(zhì)儀測量,小麥粉在儀器中依據(jù)預(yù)定的要求形成面團(tuán),粉質(zhì)儀依據(jù)形成面團(tuán)的力學(xué)特性繪制粉質(zhì)曲線[1]。粉質(zhì)曲線可以直觀的反映小麥粉的粉質(zhì)特性,其中重要的參數(shù)包括吸水率、形成時間、穩(wěn)定時間、弱化度[2]。由于粉質(zhì)特性會直接影響面團(tuán)加工品質(zhì),所以大量研究致力于通過改變面粉粉質(zhì)特性來改善面粉產(chǎn)品的最終品質(zhì)[3-5]。
近紅外光譜被廣泛應(yīng)用于預(yù)測小麥成分參數(shù),如水分、蛋白質(zhì)、灰分和濕面筋含量,而這些成分也影響著小麥粉的粉質(zhì)特性[6-7],有研究單獨對紅色硬質(zhì)春小麥的粉質(zhì)和烘焙品質(zhì)進(jìn)行了建模和評估,使用了包括近紅外光譜在內(nèi)的多種預(yù)測變量預(yù)測了吸水率、穩(wěn)定時間等多個數(shù)據(jù)[8]。很多研究雖然只是針對單一的小麥品種,但其結(jié)果也說明了近紅外技術(shù)有預(yù)測小麥粉粉質(zhì)特性的能力[9-10],由于不同年份和產(chǎn)地的小麥之間存在差異[11],這些不同的小麥制作的小麥?zhǔn)称返钠焚|(zhì)、理化性質(zhì)、保值時間等有顯著不同,這些性質(zhì)也與粉質(zhì)特性相關(guān)[12-13]。已有研究說明近紅外技術(shù)能對粉質(zhì)特性進(jìn)行有效檢測,但很少有研究對包含了多品種、不同產(chǎn)地的小麥粉樣本進(jìn)行建模分析,本研究對一個多品種,多產(chǎn)地的大規(guī)模數(shù)據(jù)集進(jìn)行建模分析,并通過將分類與回歸方法進(jìn)行融合以提高預(yù)測的精確度。
使用包括標(biāo)準(zhǔn)正態(tài)變換、多元散射矯正在內(nèi)的多種光譜預(yù)處理方法處理原始光譜[14],使用偏最小二乘法[15]、主成分分析[16]、高斯過程回歸[17]等方法對小麥粉的粉質(zhì)特性進(jìn)行回歸分析,由于研究發(fā)現(xiàn)直接應(yīng)用高斯過程回歸(GPR)時存在明顯的過擬合現(xiàn)象,所以在建模過程中只把它的結(jié)果作為分類的參考,使用高斯過程回歸加上Sigmoid函數(shù)對樣本進(jìn)行模糊分類,并在單個分簇中建立PLSR模型進(jìn)行預(yù)測。研究發(fā)現(xiàn)兩階段的回歸方法對粉質(zhì)特性的預(yù)測精度有明顯提升。
小麥粉的近紅外光譜數(shù)據(jù)采集自布勒(BUHLER)的MYRG在線近紅外分析儀,儀器的光譜范圍為850~1 650 nm,采樣間隔為5 nm,每個樣本由四個近紅外探頭平行采樣,并取四個探頭下數(shù)據(jù)的平均值作為一次采樣結(jié)果,小麥粉樣本分別在10、25和40 ℃下各進(jìn)行一次采集,最終光譜采樣數(shù)據(jù)包含269個10 ℃樣本、1106個25 ℃樣本、270個40 ℃樣本。為了避免溫度對實驗結(jié)果造成影響,同時剔除缺失信息的數(shù)據(jù),以25 ℃的968個樣本為研究對象。
968份小麥粉實驗樣本來自不同的國家和地區(qū),包括208份瑞士、51份德國、50份奧地利、46份美國、35份巴基斯坦和578份來自83個其他國家和地區(qū)。這些樣本均由布勒(BUHLER)提供。
建模所使用的開發(fā)環(huán)境為python(版本號為3.9.7),建模計算所使用的工具為scikit-learn(版本號為0.24.2)和scipy(版本號為1.7.1)。
小麥粉粉質(zhì)特性包括吸水率(water absorption)、形成時間(development time)、穩(wěn)定時間(stability time)和弱化度(degree of softening),這些數(shù)據(jù)由Brabender粉質(zhì)儀(德國)獲得。樣本集四個粉質(zhì)參數(shù)測量值的分布情況如圖1所示。
圖1 粉質(zhì)特性的頻率分布直方圖(a):吸水率;(b):形成時間;(c):穩(wěn)定時間;(d):弱化度Fig.1 Frequency distribution histogram of farinograph characteristics of wheat flour(a):Water absorption;(b):Development time;(c):Stability time;(d):Degree of softening
直接采集到的近紅外光譜數(shù)據(jù)有噪聲、基線漂移等因素的干擾,需要對其進(jìn)行預(yù)處理。預(yù)處理方法包括標(biāo)準(zhǔn)正態(tài)變換(SNV)、線性去趨勢(Detrend)、多元散射矯正(MSC)、Savitzky-Golay一階卷積求導(dǎo)(Sav-Gol D1)。嘗試使用不同的預(yù)處理組合,為不同的模型找到最優(yōu)的組合方式。近紅外光譜經(jīng)過不同預(yù)處理后的情況如圖2所示。
圖2 預(yù)處理后的近紅外光譜(a):原始光譜;(b):SNV+deternd;(c):SNV+MSC(d):SNV+Sav-Gol D1Fig.2 Preprocesses diagram of near infrared spectra(a):Raw spectra;(b):SNV+detrend;(c):SNV+MSC;(d):SNV+Sav-Gol D1
研究分別對小麥粉的4個粉質(zhì)參數(shù)建立模型,使用的回歸算法包括多元線性回歸(MLR)、偏最小二乘回歸(PLSR)、高斯過程回歸(GPR),降維算法包括主成分分析(PCA)。
1.4.1 多元線性回歸
多元線性回歸(MLR)是一種最直接的線性回歸形式,通過輸入的多個自變量直接得出因變量,其數(shù)學(xué)表達(dá)式為
Y=β0+β1X1+β2X2+…+βpXp
(1)
式(1)中,Y是因變量,X1,X2,…,Xp是自變量,β0是截距,β1,β2,…,βp為每個指標(biāo)相應(yīng)的系數(shù)。
1.4.2 偏最小二乘回歸
偏最小二乘回歸(PLSR)是一種對自變量X線性降維后再進(jìn)行回歸的建模方法。
首先對自變量矩陣和因變量矩陣進(jìn)行分解
X=DPT+E
(2)
Y=UQT+F
(3)
式中,X∈RN×M為近紅外光譜矩陣,N和M為樣本和波長的變量數(shù),D∈RN×K為主因子得分矩陣,P∈RM×K為載荷矩陣,K表示選取的主因子個數(shù),E∈RN×M表示擬合殘差矩陣,Y∈RN×K表示標(biāo)簽矩陣,U∈RN×M對應(yīng)D∈RN×K,Q∈RK×M對應(yīng)P∈RM×K,F∈RN×K對應(yīng)P∈RM×K。
矩陣D和U的線性關(guān)系為
U=DB
(4)
Ypre=DpreBQ
(5)
式(4)和式(5)中,B=(DTD)-1DTU為回歸系數(shù)矩陣。之后輸入待測樣本的光譜矩陣,根據(jù)式(2)得出其主因子得分矩陣Dpre,再根據(jù)式(5)得到待測樣本標(biāo)簽矩陣Ypre。通過對回歸系數(shù)矩陣的投影,可以將PLSR的模型表示為式(1)中所描述的形式。
1.4.3 主成分分析
主成分分析(PCA)可以將n維特征映射到k維上,這k維是全新的正交特征也被稱為主成分,是在原有n維特征的基礎(chǔ)上重新構(gòu)造出來的k維特征,包括以下步驟。近紅外數(shù)據(jù)矩陣X∈Rn×m,代表m個有n維特征變量的樣本,式(6)為近紅外數(shù)據(jù)矩陣。
(6)
式(6)中,xij為第m個樣本的第n維變量,對xij標(biāo)準(zhǔn)化處理,得
(7)
依據(jù)標(biāo)準(zhǔn)化處理后的矩陣計算相關(guān)系數(shù)矩陣A,如式(8)所示
A=(rij)m×m
(8)
最后計算累計解釋率Wi,式(9)所示
(9)
式(9)中,λi為特征矩陣對應(yīng)的特征值。
1.4.4 高斯過程回歸
高斯過程回歸(GPR)是使用高斯過程先驗對數(shù)據(jù)進(jìn)行回歸分析的非參數(shù)模型,它的定義如下。
對于所有光譜數(shù)據(jù)x=[x1,x2,…,xn],f(x)=[f(x1),f(x2),…,f(xn)]都服從多元高斯分布[18],則f為一個高斯過程,表示為
(f(x)~N[μ(x),κ(x,x)]
(10)
(11)
式(10)和式(11)中,μ(x)表示光譜數(shù)據(jù)各個維度的均值,κ(x,x)為協(xié)方差函數(shù),即核函數(shù),此處用到的核函數(shù)是高斯核函數(shù),基本形式為式(11),其中σ和l是它的超參數(shù)。
1.4.5 Sigmoid函數(shù)
Sigmoid函數(shù)是機(jī)器學(xué)習(xí)中的一個常用函數(shù),函數(shù)輸出范圍在0到1之間,在此把它作為二分類的概率預(yù)測,其數(shù)學(xué)公式為
(10)
式(10)中,x=σ為函數(shù)的閾值點,對應(yīng)輸出值為0.5,ω為函數(shù)的權(quán)值,權(quán)值的增大會使函數(shù)斜率增大,改變σ和ω可以改變函數(shù)在坐標(biāo)軸上的位置和形狀。
1.4.6 建模方法
為了更好的觀察模型效果的提升,用了三種不同的建模方法。
(1)采用偏最小二乘回歸(PLSR)直接建模。
(2)先用主成分分析(PCA)對近紅外數(shù)據(jù)降維,再使用PCA主成分做多元線性回歸(MLR)的模型。后文記為PCR。
(3)用主成分分析(PCA)對近紅外數(shù)據(jù)降維,根據(jù)PCA主成分對原數(shù)據(jù)的解釋度,選取若干個主成分作為高斯過程回歸(GPR)的輸入數(shù)據(jù)。在訓(xùn)練階段,把GPR模型的回歸結(jié)果作為參考值,并設(shè)定一個閾值將樣本分成高于閾值或者低于閾值的兩類,在兩類樣本上分別建立PLSR模型。在預(yù)測階段,將GPR模型的回歸結(jié)果代入以訓(xùn)練階段設(shè)定的分類閾值為中心的Sigmoid函數(shù)中,判斷預(yù)測樣本在兩類PLSR模型中的概率進(jìn)行模糊分類,并分別計算兩類PLSR模型的回歸結(jié)果,以分類概率進(jìn)行加權(quán)平均。此方法后文記為GPR-PLSR。
GPR-PLSR模型的樣本回歸結(jié)果公式為
(11)
在尋找最佳模型時,綜合考慮了PCA主成分的累計解釋度、PLSR模型系數(shù)和建模結(jié)果,在比較不同模型的性能時,使用校正均方根誤差(RMSEC)作為評估指標(biāo)用來評估模型預(yù)測能力,并使用相對均方根誤差(rRMSE)來協(xié)助判斷模型效果。
均方根誤差的公式為
(12)
相對均方根誤差的公式為
(13)
表1為PLSR方法和PCR方法建立的模型,針對四個粉質(zhì)參數(shù)分別建立模型,表中給出了最優(yōu)情況下的模型參數(shù)。從中得出,SNV+Sav-GolD1的預(yù)處理方法為除了形成時間外的三個參數(shù)的預(yù)測提供了最好的效果,對比兩種建模方法,PLSR模型使用了更少的因子數(shù)而建模效果卻略好于PCR模型,可以看出PLSR的主因子比PCA的主成分擁有更強(qiáng)的解釋能力[19]。
表1 PLSR和PCR的建模結(jié)果Table 1 Modeling results of PLSR and PCR
PLSR模型的回歸系數(shù)如圖3所示,其中紅色曲線代表表1中最優(yōu)因子數(shù)的PLSR模型的回歸系數(shù)曲線,綠色曲線比紅色曲線少一個因子數(shù),藍(lán)色曲線比紅色曲線多一個因子數(shù)。在吸水率的PLSR模型上,三個不同因子數(shù)的曲線較為相似。選擇RMSEC最小的紅色曲線模型作為最優(yōu)模型,而在其他三個粉質(zhì)參數(shù)模型上,雖然藍(lán)色曲線代表的模型得到了更小的RMSEC,但是它與其他兩條曲線相比噪聲顯著增加,模型魯棒性更低,而紅色和綠色曲線較為相似,但紅色曲線的RMSEC更小,最終選擇了紅色曲線代表的模型作為最優(yōu)模型。
圖3 PLSR模型的回歸系數(shù)(a):吸水率;(b):形成時間;(c):穩(wěn)定時間;(d):弱化度Fig.3 Regression coefficients of PLSR model(a):Water absorption;(b):Development time;(c):Stability time;(d):Degree of softening
表2為GPR-PLSR方法建立的模型,可以得出,SNV+MSC的預(yù)處理方法為除了吸水率外的三個參數(shù)的預(yù)測提供了最好的效果,這與前文PLSR方法和PCR方法建立的模型不同。在GPR中,使用PCA的主成分作為變量輸入。表中給出了這些變量的累計解釋度,由于GPR方法的擬合能力很強(qiáng)大,它的預(yù)測結(jié)果表現(xiàn)出明顯的過擬合現(xiàn)象,所以只用它的預(yù)測結(jié)果作分類的參考。
表2 GPR-PLSR的建模結(jié)果Table 2 Modeling results of GPR-PLSR
通過調(diào)整Sigmoid函數(shù)的閾值σ和權(quán)值ω使其與四個粉質(zhì)參數(shù)在坐標(biāo)軸上的分布相適配,并且閾值σ將GPR的預(yù)測結(jié)果分成兩簇,第一簇數(shù)據(jù)的GPR預(yù)測結(jié)果低于閾值σ,第二簇數(shù)據(jù)的GPR預(yù)測結(jié)果高于閾值σ,圖4所示為Sigmoid函數(shù)結(jié)果,橫坐標(biāo)是樣本的真實標(biāo)簽,縱坐標(biāo)是模型預(yù)測在兩個簇上的概率,可以看出樣本在四個參數(shù)上的分類結(jié)果都能較好的擬合出Sigmoid函數(shù)的形狀[20]。
圖4 Sigmoid函數(shù)概率圖(a):吸水率;(b):形成時間;(c):穩(wěn)定時間;(d):弱化度Fig.4 Probabilities of Sigmoid function(a):Water absorption;(b):Development time;(c):Stability time;(d):Degree of softening
用分成兩簇的數(shù)據(jù)分別建立兩個PLSR模型,這兩個PLSR模型的主因子數(shù)與上文PLSR方法的建模結(jié)果一致,圖5所示為這兩個PLSR模型的擬合效果,紅綠兩色的樣本代表GPR分類器預(yù)測出的類別。樣本分布基本遵循預(yù)定的閾值分類規(guī)則。最后,將整個數(shù)據(jù)集放入兩個模型中得到不同的兩個預(yù)測結(jié)果,再以Sigmoid函數(shù)得出的樣本落在兩個模型中的概率為權(quán)重,依據(jù)式(11)得到GPR-PLSR模型的預(yù)測結(jié)果。四個粉質(zhì)參數(shù)模型的校正均方根誤差(RMSEC)如表2所示,其預(yù)測值與真實值的分布如圖6所示,其中紅色散點代表GPR-PLSR模型預(yù)測值和真實值的分布,作為對照,藍(lán)色散點代表2.1節(jié)中描述的PLSR模型預(yù)測值與真實值的分布。
圖5 兩個分簇后PLSR模型的預(yù)測回歸圖(a):吸水率;(b):形成時間;(c):穩(wěn)定時間;(d):弱化度Fig.5 Predictions of PLSR model based on clustered regression(a):Water absorption;(b):Development time;(c):Stability time;(d):Degree of softening
圖6 PLSR模型和GPR-PLSR模型的預(yù)測回歸圖(a):吸水率;(b):形成時間;(c):穩(wěn)定時間;(d):弱化度Fig.6 Predictions of PLSR and GPR-PLSR regression models(a):Water absorption;(b):Development time;(c):Stability time;(d):Degree of softening
由于PLSR方法建立的模型與PCA方法建立的模型預(yù)測能力相當(dāng),并且前者RMSEC略好于后者,所以將PLSR模型與GPR-PLSR模型進(jìn)行比較。圖6中吸水率模型的擬合程度最高,兩種模型的散點基本沿著y=x直線線性分布,這兩種模型對吸水率的預(yù)測準(zhǔn)確度相當(dāng),從表1和表2也能看出,PLSR模型的RMSEC為2.039,GPR-PLSR模型為1.876。其他三個粉質(zhì)參數(shù)的模型,GPR-PLSR模型的擬合能力明顯更強(qiáng),相對于PLSR模型,GPR-PLSR模型的散點明顯更貼近的分布在y=x直線附近,模型的RMSEC也從1.838、4.037和21.693降低到了1.160、2.459和14.449,同時,rRMSE也反映了模型擬合能力的增強(qiáng)。
通過引入基于高斯過程的二階段回歸模型,來提高近紅外光譜對小麥粉粉質(zhì)特性的預(yù)測精度,并在一個大規(guī)模近紅外光譜數(shù)據(jù)集上進(jìn)行了驗證和對比。包括PLSR模型在內(nèi)的線性模型在近紅外光譜分析上有廣泛的應(yīng)用,但是在樣本種類豐富,組成復(fù)雜的數(shù)據(jù)集上,并不能得到較為精確的預(yù)測結(jié)果。提出的GPR-PLSR模型,先基于高斯過程回歸結(jié)果對樣本進(jìn)行模糊分類,再對不同區(qū)間內(nèi)的PLSR子模型的預(yù)測結(jié)果進(jìn)行加權(quán)平均作為模型最后的預(yù)測結(jié)果,相比于經(jīng)典PLSR模型有較大的提升。吸水率、形成時間、穩(wěn)定時間和弱化度四個小麥粉粉質(zhì)特性的預(yù)測,GPR-PLSR模型的RMSEC從2.039、1.838、4.037和21.693降低到了1.876、1.160、2.459和14.449??梢詾榻t外預(yù)測小麥粉粉質(zhì)特性提供一些技術(shù)參考。