趙靜遠(yuǎn)， 熊智新*， 寧井銘， 謝德紅

(1.南京林業(yè)大學(xué)輕工與食品學(xué)院，江蘇南京 210037；2.安徽農(nóng)業(yè)大學(xué)茶樹生物學(xué)與資源利用國(guó)家重點(diǎn)實(shí)驗(yàn)室，安徽合肥 230036)

茶葉中的咖啡堿是茶葉生物堿的主體成分，約占干物質(zhì)質(zhì)量的2%～4%，其有興奮、利尿、解毒、強(qiáng)心，促進(jìn)血液循環(huán)的作用[1 - 3]。咖啡堿可以作為茶葉的特征性化學(xué)物質(zhì)之一，可用于判斷茶葉的真假[4]。近紅外光譜(Near Infrared Spectroscopy,NIRS)作為一種高效、快速、低成本、無(wú)污染的在線檢測(cè)的分析技術(shù)，在農(nóng)業(yè)[5]、醫(yī)藥[6]、煙草[7]和食品[8]等行業(yè)得到廣泛地重視和應(yīng)用。目前，國(guó)內(nèi)外學(xué)者已利用近紅外光譜對(duì)茶葉進(jìn)行了定性分類和定量研究[9 - 13]。利用近紅外光譜定量分析咖啡堿依賴于模型的可靠性，但是近紅外光譜信息復(fù)雜和全譜建模的預(yù)測(cè)能力不足會(huì)影響模型的穩(wěn)定，通常采用波長(zhǎng)變量選擇方法，篩選出信息強(qiáng)且對(duì)外界影響因素不敏感的波長(zhǎng)，從而建立更穩(wěn)健的模型[14]。已有文獻(xiàn)報(bào)道[15 - 18]連續(xù)投影算法(Successive Projections Algorithm,SPA)能夠有效消除眾多波長(zhǎng)變量間的共線性影響，降低模型復(fù)雜度，以其簡(jiǎn)便、快速的特點(diǎn)得到廣泛應(yīng)用，在多種樣品組分的特征波長(zhǎng)選取中取得了很好效果。

近紅外漫反射光譜由于受到儀器條件、測(cè)量環(huán)境、樣品狀態(tài)等諸多因素影響，其中含有各種高頻的隨機(jī)噪聲、基線漂移、信號(hào)本體、樣品不均勻和光散射等干擾[19]。已有研究多采用多元散射校正等[20，21]方法盡可能消除這些因素對(duì)光譜的干擾，但究竟應(yīng)該采用何種預(yù)處理方法更有利于SPA算法挑選波長(zhǎng)，目前所見相關(guān)報(bào)道較少。本文研究重點(diǎn)即采用小波變換(Wavelet Transform,WT)、多元散射校正(Multiplicative Scatter Correction,MSC)、標(biāo)準(zhǔn)正態(tài)變量變換(Standard Normal Variate Transformation,SNV)以及一階微分(First Derivative,1stD)等預(yù)處理方法單獨(dú)或組合對(duì)光譜進(jìn)行預(yù)處理，再結(jié)合SPA算法建立茶葉中咖啡堿的偏最小二乘回歸(Partial Least Square Regression,PLSR)預(yù)測(cè)模型。通過(guò)比較各模型的精度，探討SPA算法和各種預(yù)處理方法結(jié)合的最佳途徑，以期在消除各種低頻和高頻噪聲干擾的基礎(chǔ)上，充分發(fā)揮SPA波長(zhǎng)選擇能力，獲得最佳共線性的波長(zhǎng)變量子集，提高所建模型的質(zhì)量。

1 實(shí)驗(yàn)部分

1.1 儀器及試劑

德國(guó)布魯克公司生產(chǎn)的MPA型多功能傅里葉變換近紅外光譜儀，帶有RT-PbS檢測(cè)器、OPUS/OVP(OPUS validation program)自檢功能，附有內(nèi)置鍍金漫反射積分球、樣品旋轉(zhuǎn)器和石英樣品杯等配件。

茶葉樣本包括2007年度不同地區(qū)的高、中、低檔茶葉樣本158個(gè)，由寧波市進(jìn)出口檢驗(yàn)檢疫局、中國(guó)農(nóng)業(yè)科學(xué)院茶葉研究所和鎮(zhèn)江市林業(yè)局提供。其中包括龍井、碧螺春等優(yōu)質(zhì)名茶、大宗茶，市售一般綠茶等。由于樣本大多數(shù)取自于各茶廠的送檢樣品，故每份樣品只編號(hào)未標(biāo)注具體的茶葉品種和送檢單位。同時(shí)為擴(kuò)大模型的適用性，實(shí)驗(yàn)中未細(xì)分綠茶品種。所有送檢樣品均已密封保存4個(gè)月。

1.2 實(shí)驗(yàn)方法

1.2.1 咖啡堿含量測(cè)定稱取茶葉樣品30±1 g，置于HY-02小型磨碎機(jī)(北京環(huán)亞)中研磨10 s，按照我國(guó)國(guó)家標(biāo)準(zhǔn)(GB/T 8303-2012)[22]，將磨碎樣品過(guò)100目篩(150 μm)，按照國(guó)家標(biāo)準(zhǔn)(GB/T 8312-2013)[23]規(guī)定方法測(cè)定茶葉樣品中咖啡堿的含量，每個(gè)樣品進(jìn)行3次平行檢測(cè)。

1.2.2 光譜采集茶葉樣品光譜采集范圍：12 800～4 000 cm-1，分辨率：4 cm-1。實(shí)驗(yàn)溫度：20 ℃，相對(duì)濕度：40%。均以儀器內(nèi)置鍍金為背景。因在8 000～4 000 cm-1波段茶葉近紅外光譜圖中的波峰較多且比較明顯，故截取該波長(zhǎng)范圍內(nèi)共1 038個(gè)波長(zhǎng)點(diǎn)用于分析研究。圖1為部分樣品的近紅外光譜圖。

圖1 部分茶葉樣品的近紅外光譜圖Fig.1 Near-infrared spectra of tea samples

1.2.3 樣品集的劃分本文采用主成分分析(Principal Component Analysis,PCA)，利用提取的光譜主成分得分向量來(lái)代替光譜向量計(jì)算樣本的馬氏距離(Mahalanobis Distance)[24]，將超出設(shè)定的馬氏距離閾值的異常樣本進(jìn)行剔除，最終從158個(gè)樣本中剔除兩個(gè)異常樣品，然后采用SPXY[14]將剩余的156個(gè)樣本劃分為校正集和預(yù)測(cè)集。SPXY在計(jì)算樣品距離時(shí)同時(shí)將變量x和變量y考慮在內(nèi)，分別計(jì)算各變量和變量之間的距離，分別見式(1)和式(2)：

(1)

(2)

在式(1)和式(2)中，n和m表示任意兩個(gè)樣品，Z是總的樣本數(shù)，l為光譜的波數(shù)點(diǎn)，dx(n,m)表示兩條光譜的空間距離，dy(n,m)表示對(duì)應(yīng)咖啡堿之間的距離。

兩個(gè)樣品間的SPXY標(biāo)準(zhǔn)距離計(jì)算公式如式(3)所示。

(3)

最終156個(gè)樣本用SPXY法以3∶2的比例，劃分為93個(gè)校正集樣本和63個(gè)預(yù)測(cè)集樣本。校正集用于建立茶葉咖啡堿近紅外模型，預(yù)測(cè)集用于驗(yàn)證所建模型的準(zhǔn)確性和可靠性。樣品集劃分結(jié)果見表1。

表1 茶葉樣本中咖啡堿的含量分布Table 1 Content of caffeine in tea samples

1.2.4 光譜預(yù)處理方法分別采用WT、MSC、SNV以及1stD對(duì)光譜進(jìn)行預(yù)處理，探討各種方法及它們之間的組合對(duì)模型精度的影響。小波變換的實(shí)質(zhì)[25]是將信號(hào)x(t)投影到小波基函數(shù)ψa,b(t)上，即x(t)與ψa,b(t)的內(nèi)積，得到小波變換系數(shù)。按照分析的需要對(duì)小波變換系數(shù)進(jìn)行處理，然后利用處理后的小波系數(shù)得到反變換信號(hào)。小波變換用于處理光譜原始數(shù)據(jù)的一般方法為：首先對(duì)原始光譜進(jìn)行小波分解變換分解得到高頻和低頻小波系數(shù)；然后通過(guò)閾值法去除小波系數(shù)中被認(rèn)為是表示噪聲的元素(濾噪)，或去除小波系數(shù)中的高頻(低尺度)元素；最終用經(jīng)過(guò)處理系數(shù)進(jìn)行反變換即可得到濾噪后的光譜信號(hào)。

由于小波變換中用到的小波函數(shù)不具有唯一性，同一問(wèn)題用不同的小波函數(shù)進(jìn)行分析有時(shí)結(jié)果相差甚遠(yuǎn)，因此小波函數(shù)的選用是小波變換實(shí)際應(yīng)用中的一個(gè)難點(diǎn)，目前都是基于經(jīng)驗(yàn)或不斷實(shí)驗(yàn)后通過(guò)結(jié)果比較來(lái)選擇最佳的小波函數(shù)[14]。本文即采取實(shí)驗(yàn)方法選擇不同小波函數(shù)處理茶葉近紅外光譜并建模，結(jié)果表明采用db3小波函數(shù)并分解為7層能較好地分辨光譜低頻和高頻信號(hào)并予以處理，建模精度較高。本文后續(xù)即選擇db3小波7層分解開展光譜數(shù)據(jù)處理研究。

1.2.5 變量篩選方法SPA是一種前向循環(huán)變量選擇方法[26]。從選定一個(gè)波長(zhǎng)開始，每次循環(huán)都計(jì)算在未入選波長(zhǎng)上的投影，將投影向量最大的波長(zhǎng)引入到波長(zhǎng)組合。每一個(gè)新入選的波長(zhǎng)，都與前一個(gè)線性關(guān)系最小。在波長(zhǎng)選擇中，假設(shè)n為初始波長(zhǎng)數(shù)，為m波長(zhǎng)點(diǎn)個(gè)數(shù)，k(0)為初始波長(zhǎng)點(diǎn)。SPA的計(jì)算步驟如下：

(1)第一次迭代(Z=1)開始前，在光譜矩陣中任選一列波長(zhǎng)向量Xj。記為Xk(0)，即k(0)=j，j∈1,…,n；

(2)將未入選的列向量位置集合記為S，S={j,1≤j≤n,j?{k(0),…k(Z-1)}}；

(4)提取最大的投影值的波長(zhǎng)變量序號(hào)：k(Z)=arg[max(‖Pxj‖)]，j∈S；

(5)令xj=Pxj，j∈S；

(6)Z=Z+1，如果Z

{k(Z),Z=0,…,m-1}為最終選取的波長(zhǎng)變量組合。算法對(duì)每次選擇的結(jié)果進(jìn)行建模預(yù)測(cè)分析，以最小均方根誤差(RMSEP)來(lái)決定所建模型的優(yōu)劣，它所對(duì)應(yīng)的即為波長(zhǎng)集合的最佳初始波長(zhǎng)點(diǎn)n和最佳波長(zhǎng)點(diǎn)個(gè)數(shù)m。

2 結(jié)果與討論

2.1 不同預(yù)處理方法結(jié)果分析

研究中采用粉末狀茶葉樣品測(cè)量近紅外光譜，因此由茶粉顆粒分布不均勻及顆粒大小不同產(chǎn)生的散射誤差將是信號(hào)中的主要干擾，其通常表現(xiàn)為基線的平移、旋轉(zhuǎn)(低頻)和二次和高次曲線(高頻)等[25]。在各種預(yù)處理方法中，SNV、MSC通常用于校正低頻散射信號(hào)、1stD可以消除信號(hào)的旋轉(zhuǎn)，而WT基于多分辨分析原理，則可以同時(shí)消除光譜信號(hào)中的高頻和低頻噪聲。分別以WT、MSC、SNV和1stD預(yù)處理光譜后采用PLSR建模統(tǒng)計(jì)結(jié)果列于表2。使用驗(yàn)證集相對(duì)分析誤差(RPD)對(duì)模型進(jìn)一步評(píng)價(jià)，如果RPD≥3.0，說(shuō)明定標(biāo)效果良好，建立的定標(biāo)模型可以用于實(shí)際檢測(cè);如果2.5

表2 采用不同預(yù)處理方法的全波段PLSR模型性能Table 2 Performance of full-band PLSR models using different preprocessing methods

2.2 不同預(yù)處理方法結(jié)合波長(zhǎng)選擇結(jié)果分析

分別用WT、MSC、SNV和1stD處理近紅外光譜，再利用SPA算法挑選波長(zhǎng)并建立相應(yīng)的PLSR預(yù)測(cè)模型，模型的預(yù)測(cè)結(jié)果見表3。從表3可看出，通過(guò)三種預(yù)處理方法結(jié)合SPA波長(zhǎng)選擇方法，能有效地減少入選波長(zhǎng)數(shù)量(從原始1 038個(gè)波長(zhǎng)減少到31個(gè)以下)，且模型的精度總體上也有了一定地提升，其中1stD入選波長(zhǎng)最少(11個(gè))。

表3 不同預(yù)處理方法結(jié)合SPA的PLSR模型預(yù)測(cè)結(jié)果Table 3 PLSR model prediction results of different pretreatment methods combined with SPA

圖2分別為WT、MSC和SNV分別結(jié)合1stD預(yù)處理方后的光譜圖，為便于比對(duì)，圖中還列出先1stD再WT處理結(jié)果(圖2(b))。表4則是在1stD+SPA基礎(chǔ)上分別與MSC、SNV、WT等預(yù)處理方法結(jié)合后的建模效果。

表4 小波變換結(jié)合1stD+SPA的PLSR模型預(yù)測(cè)結(jié)果及比較Table 4 Prediction results and comparison of PLSR models combined with wavelet transform and 1stD+SPA

從圖2中可以看出，茶粉近紅外光譜經(jīng)WT-1stD或1stD -WT預(yù)處理后在6 000～4 000 cm-1波段有較明顯的吸收峰，幾乎不受噪聲干擾，而MSC-1stD和SNV-1stD處理后峰信息則完全淹沒(méi)在噪聲信號(hào)中。從表4中可以看出，WT結(jié)合1stD(WT-1stD或1stD -WT)最終建模效果總體優(yōu)于其它幾種預(yù)處理方法，而又以WT-1stD -SPA-PLSR組合所建模型最好，且只用10個(gè)入選波長(zhǎng)建模。這說(shuō)明利用WT預(yù)先消除光譜信號(hào)中的大部分高頻噪聲，有利于1stD方法提高分辨率和提供更清晰的光譜輪廓變化，減少近紅外原始重疊譜帶波長(zhǎng)之間高共線性，挑選出更有代表性的特征波長(zhǎng)。其挑選出的波長(zhǎng)數(shù)10和通過(guò)交叉驗(yàn)證確定的最佳PLS成分?jǐn)?shù)7也比較接近，表明WT-1stD預(yù)處理后采用SPA方法所選波長(zhǎng)具有較好的代表性和相對(duì)獨(dú)立性，由此建立的模型比較穩(wěn)健，故其校正精度和預(yù)測(cè)精度均較高。

從表4還可以看出，1stD和MSC、SNV、WT結(jié)合處理先后順序?qū)δＰ徒Y(jié)果影響明顯，而且先1stD后進(jìn)行其它預(yù)處理總體建模效果略差。其原因可能在于原始信號(hào)經(jīng)1stD處理后，有效信號(hào)和被放大的高頻噪聲混疊嚴(yán)重，后續(xù)的WT、SNV和MSC分別采用頻域方法和統(tǒng)計(jì)方法抑制高頻信號(hào)的同時(shí)也削弱了有效信號(hào)。比較圖2(a)和圖2(b)可以看出，對(duì)于1stD -WT處理，由于WT對(duì)1stD放大的噪聲以及微分信號(hào)都進(jìn)行了去高頻處理，故其信號(hào)較WT-1stD平滑，可能損失部分譜峰信息而降低建模效果。

圖2 茶葉近紅外光譜分別經(jīng)WT、MSC和SNV并結(jié)合1stD預(yù)處理圖Fig.2 The near-infrared spectra of tea pretreated by WT,MSC and SNV combined with 1stD

圖3是不同方法與1stD及SPA組合預(yù)處理光譜并建模計(jì)算RMSEP值的結(jié)果。由圖3(a)可以看出，經(jīng)過(guò)對(duì)校正集數(shù)據(jù)進(jìn)行WT-1stD預(yù)處理后，在SPA選擇過(guò)程中均方根誤差的變化情況。在選擇4個(gè)光譜波長(zhǎng)變量之前曲線下降趨勢(shì)很大，表明茶葉中咖啡堿含量的光譜波長(zhǎng)變量應(yīng)該至少選擇4個(gè)以免產(chǎn)生過(guò)擬合的問(wèn)題；選擇4個(gè)到10個(gè)光譜波長(zhǎng)變量時(shí)均方根誤差曲線下降至最低值；之后選擇10個(gè)以上的光譜波長(zhǎng)變量時(shí)，曲線雖有升降但總體趨于平緩。最終選擇10個(gè)光譜波長(zhǎng)變量，其均方根誤差(RMSEP)達(dá)到最小為0.2794。其它幾種預(yù)處理方法組合也有類似趨勢(shì)，但最終入選波長(zhǎng)均大于10個(gè)，但大都少于未用1stD處理入選波長(zhǎng)數(shù)(表3)。這都說(shuō)明正是1stD提高了光譜分辨率，利于后續(xù)SPA算法挑選出少數(shù)特征波長(zhǎng)。

圖3 WT、MSC和SNV結(jié)合1stD及SPA選取特征波長(zhǎng)數(shù)及RMSEP值Fig.3 The number of selected characteristic wavelength and RMSEP by WT,MSC and SNV combined with 1stD and SPA

圖4譜圖上標(biāo)出了經(jīng)WT-1stD -SPA預(yù)處理并建立咖啡堿近紅外分析PLSR模型所選擇的10個(gè)波長(zhǎng)點(diǎn)(4 003、4 350、4 473、4 627、4 812、5 175、5 244、5 275、5 583和5 984 cm-1)位置?？Х葔A又名1,3,7-三甲基黃嘌呤，其碳氮嘌呤環(huán)結(jié)構(gòu)由一個(gè)嘧啶環(huán)和一個(gè)咪唑環(huán)稠合而成，雜環(huán)上含有三個(gè)-CH3和兩個(gè)C=O[4]。由于茶葉作為組分復(fù)雜的自然產(chǎn)物，各有機(jī)組分官能團(tuán)在近紅外光譜帶重疊嚴(yán)重，一般很難清晰地指出各譜峰對(duì)應(yīng)的具體基團(tuán)信息，但根據(jù)近紅外主要譜帶歸屬仍可對(duì)所選波長(zhǎng)合理性做出必要分析。由于甲基吸收譜帶位于4 200～4 440 cm-1，為C-H鍵伸縮振動(dòng)與彎曲振動(dòng)的合頻區(qū)；羰基吸收譜帶位于5 050～5 210 cm-1，為C=O鍵伸縮振動(dòng)的二級(jí)倍頻區(qū)[14]。而且由于核酸中的嘌呤核苷酸，蛋白質(zhì)中的蛋氨酸、甘氨酸、谷酰胺、天冬氨酸等都和咖啡堿的合成有關(guān)[4]，因此可以推斷在蛋白質(zhì)的特征吸收區(qū)4 545～4 655 cm-1、4 854～4 878 cm-1，以及蛋白質(zhì)仲酰胺(CONH)中的C=O伸縮振動(dòng)二級(jí)倍頻吸收(5 208 cm-1)等[28]波段附近應(yīng)該存在茶葉咖啡堿的特征吸收。對(duì)照已選出的10個(gè)波長(zhǎng)點(diǎn)可以看出，它們大部分都落在以上譜帶區(qū)域內(nèi)或附近，表明挑選出的各峰位較好地反映了咖啡堿的吸收特征[29,30]。

圖4 最佳模型選取的最優(yōu)特征波長(zhǎng)點(diǎn)Fig.4 Optimal characteristic wavelength points selected by the best model

3 結(jié)論

(1)利用全波段建立茶葉咖啡堿近紅外分析模型，處理數(shù)據(jù)量大，建模時(shí)間長(zhǎng)，同時(shí)可能引入更多的背景干擾。SPA方法則可以通過(guò)挑選特征波長(zhǎng)，簡(jiǎn)化模型，并消除多重共線性，有利于提高模型的簡(jiǎn)潔性和精度，但該方法易受到重疊峰、高低頻噪聲和基線漂移的影響。

(2)采用WT-1stD相結(jié)合預(yù)處理光譜，可實(shí)現(xiàn)兩者優(yōu)勢(shì)互補(bǔ)，即WT有去高、低頻噪聲及基線漂移的作用，1stD則進(jìn)一步提高重疊峰的分辨率，方便后續(xù)SPA算法挑選更具代表性的最小共線性的波長(zhǎng)變量子集。建模實(shí)驗(yàn)表明，WT-1stD預(yù)處理方法顯著地提高了茶葉咖啡堿近紅外SPA+PLSR預(yù)測(cè)模型的精度。

(3)通過(guò)WT-1stD預(yù)處理后，SPA算法提取特征波長(zhǎng)較好地分布于茶葉咖啡堿主要官能團(tuán)吸收區(qū)域，提高了咖啡堿近紅外分析模型的可解釋性。