張若秋，杜一平

(華東理工大學(xué) 上海市功能性材料化學(xué)重點(diǎn)實(shí)驗(yàn)室 化學(xué)與分子工程學(xué)院，上海 200237)

近紅外(Near-infrared,NIR)光譜分析是近年來(lái)最引人關(guān)注的分析技術(shù)之一，在農(nóng)業(yè)、制藥、化工等眾多領(lǐng)域得到廣泛應(yīng)用[1-2]。目前NIR技術(shù)的發(fā)展重點(diǎn)是應(yīng)用，簡(jiǎn)便、快速、無(wú)損、準(zhǔn)確等優(yōu)點(diǎn)正是其受重視并被廣泛應(yīng)用于國(guó)民經(jīng)濟(jì)發(fā)展各個(gè)領(lǐng)域的本質(zhì)原因。但在使用NIR技術(shù)的基層單位，技術(shù)隊(duì)伍的研究能力往往不如研發(fā)單位，更不及高校和專(zhuān)門(mén)的科研機(jī)構(gòu)。而NIR技術(shù)使用和推廣的難點(diǎn)是近紅外光譜分析模型的建立需要一些專(zhuān)門(mén)的技術(shù)，尤其是依賴(lài)于化學(xué)計(jì)量學(xué)。因此，模型的建立和維護(hù)是近紅外光譜技術(shù)發(fā)展的一個(gè)突出難點(diǎn)和瓶頸。

偏最小二乘(Partial least squares,PLS)是近紅外光譜多元校正領(lǐng)域使用最為廣泛的算法之一[3-4]。在光譜多元校正的實(shí)際應(yīng)用中，有很多因素會(huì)影響該算法的建模過(guò)程和模型的預(yù)測(cè)能力，其中包括光譜儀器噪聲。然而，目前絕大多數(shù)研究都針對(duì)如何使用濾波器或平滑算法減少光譜噪聲[5-6]，鮮有針對(duì)儀器噪聲如何影響偏最小二乘建模過(guò)程及其預(yù)測(cè)能力的研究。本文重點(diǎn)闡述和討論儀器噪聲怎樣通過(guò)第一個(gè)隱變量的計(jì)算被引入模型中，并隨著后續(xù)隱變量不斷在模型中傳遞和累積從而對(duì)整個(gè)建模過(guò)程產(chǎn)生影響，可為今后進(jìn)一步研究在建模中減少或抑制噪聲在模型內(nèi)的傳播和放大提供理論依據(jù)。

1 儀器噪聲對(duì)偏最小二乘建模的影響

樣本的測(cè)量光譜用X表示，組分含量用y表示，多元校正關(guān)系模型為y=Xβ+e。考慮測(cè)量光譜存在誤差，將X表示為：

X=X0+E

(1)

其中，X0為不含測(cè)量誤差的真光譜，E表示儀器噪聲，假設(shè)矩陣中的每個(gè)元素均滿(mǎn)足均值為0且方差相同的高斯分布。類(lèi)似地，含量真值用y0表示(本文只考慮光譜存在誤差，含量無(wú)誤差，即y0=y)。

1.1 偏最小二乘算法實(shí)現(xiàn)

PLS的主流實(shí)現(xiàn)方式有非線(xiàn)性迭代偏最小二乘(NIPALS)[7]和SIMPLS兩種，為方便理論推導(dǎo)本文使用SIMPLS，下文給出SIMPLS的偽代碼[8]。SIMPLS算法流程如下：

Fora=1,…,A;

a=1:對(duì)S0進(jìn)行奇異值分解;

a>1:對(duì)(I-P(PTP)-1PT)S0進(jìn)行奇異值分解;

計(jì)算得到權(quán)重向量r=第一個(gè)左奇異向量;

計(jì)算得到得分向量t=X0r;

將r、t和p分別存入R、T和P;

End

無(wú)噪光譜矩陣X0含有p個(gè)變量和n個(gè)樣本，濃度向量y0含有n個(gè)樣本。R、T和P分別代表權(quán)重矩陣、得分矩陣和載荷矩陣，而βpls則為偏最小二乘模型的回歸系數(shù)向量。為方便，本文使用LSV1(S0)表示S0的第一個(gè)左奇異向量。

1.2 儀器噪聲傳遞

當(dāng)PLS模型的隱變量數(shù)量(nLVs)為k個(gè)時(shí)，在不同隱變量下所形成的空間和由Krylov序列形成的空間滿(mǎn)足如下等式[5]：

{Rk}={S0,XTXS0,(XTX)2S0,…,(XTX)k-1S0}

(2)

{Pk}={XTXS0,(XTX)2S0,…,(XTX)kS0}

(3)

以上兩個(gè)等式表明協(xié)方差矩陣S0(當(dāng)預(yù)測(cè)指標(biāo)只有1個(gè)時(shí)是協(xié)方差向量)和XTX完全確定了βpls。前者確定模型的初始方向，而后者則影響每個(gè)隱變量的實(shí)際計(jì)算。當(dāng)光譜數(shù)據(jù)中包含了儀器噪聲時(shí)，S0和pk的計(jì)算可改寫(xiě)為如下形式：

(4)

(5)

因此，通過(guò)上述公式發(fā)現(xiàn)偏最小二乘模型將從第一個(gè)隱變量開(kāi)始發(fā)生偏離真模型的現(xiàn)象，噪聲同時(shí)被引入到模型中，并在隨后隱變量的計(jì)算過(guò)程中被不斷傳遞和累積：

(6)

Rk=[LSV1(S1),LSV1(S2),…,LSV1(Sk)]

(7)

βpls=RkRkTS0

(8)

上述噪聲傳遞過(guò)程同時(shí)也會(huì)對(duì)偏最小二乘模型的預(yù)測(cè)能力產(chǎn)生影響。預(yù)測(cè)殘差可表示為：

(9)

(10)

(11)

(12)

圖1 噪聲傳遞和偏最小二乘預(yù)測(cè)誤差關(guān)系示意圖Fig.1 The scheme of the relationship between noise propagation and prediction error

噪聲傳遞和偏最小二乘的預(yù)測(cè)誤差關(guān)系如圖1所示，不同噪聲水平下的最優(yōu)隱變量數(shù)也將發(fā)生改變，噪聲水平越高，最優(yōu)隱變量的數(shù)值越小。這說(shuō)明噪聲傳遞不僅影響偏最小二乘模型的預(yù)測(cè)能力，還會(huì)給模型的選擇帶來(lái)影響。

2 數(shù)據(jù)集與計(jì)算方法

2.1 數(shù)據(jù)集模擬

本文用模擬數(shù)據(jù)對(duì)上述討論進(jìn)行驗(yàn)證。模擬數(shù)據(jù)集包含了2 000個(gè)建模樣本和2 000個(gè)獨(dú)立預(yù)測(cè)樣本，其變量數(shù)為1 000。數(shù)據(jù)集中包含1個(gè)待測(cè)組分和40個(gè)干擾組分，共41個(gè)組分。每個(gè)模擬組分的純光譜信號(hào)均滿(mǎn)足一個(gè)隨機(jī)均值和方差的高斯分布。建模集中待測(cè)組分和干擾組分的濃度值服從0到1的均勻分布，而在預(yù)測(cè)集中服從0.05到0.95的均勻分布。模擬光譜由濃度矩陣和純光譜信號(hào)矩陣相乘得到。

對(duì)于每個(gè)固定的光譜信號(hào)xi(i=1,2,…,n)均產(chǎn)生一個(gè)服從標(biāo)準(zhǔn)高斯分布的噪聲向量ei，噪聲水平的大小定義為ei的最大值和xi的比值。更高的噪聲水平意味著更低的信噪比(S/N)，反之亦然。3個(gè)不同水平的噪聲(0.3%、0.75%和1%)用于給無(wú)噪光譜數(shù)據(jù)添加人為噪聲。

2.2 指標(biāo)與計(jì)算軟件評(píng)價(jià)

評(píng)價(jià)所建立的NIR模型通常考慮模型誤差和相關(guān)系數(shù)兩個(gè)指標(biāo)，本文只考慮模型誤差指標(biāo)。建模時(shí)一般將數(shù)據(jù)集劃分為校正集(C)、交互檢驗(yàn)集(CV)和預(yù)測(cè)集(P)，因此模型誤差也包括校正誤差、交互檢驗(yàn)誤差和預(yù)測(cè)誤差。模型誤差用均方根誤差(Root mean squared error,RMSE)表示：

(13)

已有不少研究表明蒙特-卡洛采樣(Monte-carlo sampling,MCS)能夠降低模型過(guò)擬合的風(fēng)險(xiǎn)[9-10]，因此本文采用基于MCS的方法，即采樣誤差分布分析的交互檢驗(yàn)(Cross-validation based on sampling error profile analysis,SEPA-CV)[11]確定偏最小二乘的最優(yōu)隱變量大小。MCS的采樣數(shù)和建模樣本比例分別設(shè)定為2 000和0.8。所有計(jì)算均使用MATLAB完成(Version 2010a,The MathWorks,USA)。

為準(zhǔn)確量化偏最小二乘模型的噪聲水平，采用Durbin-Watson(DW)[12-16]統(tǒng)計(jì)量評(píng)價(jià)噪聲對(duì)模型的影響程度。DW統(tǒng)計(jì)量由下式計(jì)算得到：

(14)

xi和xi-1是p維向量中的兩個(gè)連續(xù)元素。如果整個(gè)向量中連續(xù)元素間的相關(guān)性很弱，例如該向量中每個(gè)元素均為隨機(jī)變量，則DW統(tǒng)計(jì)量的大小將收斂于2。一個(gè)擁有至少100個(gè)元素的隨機(jī)向量的DW統(tǒng)計(jì)量的95%置信區(qū)間在1.7和2.3之間。

3 結(jié)果與討論

3.1 噪聲傳遞對(duì)偏最小二乘模型回歸向量的影響

首先使用DW統(tǒng)計(jì)量對(duì)3種不同噪聲水平的光譜所建立的偏最小二乘模型的權(quán)重向量r、載荷向量p和回歸系數(shù)向量βpls進(jìn)行噪聲程度的估計(jì)，其結(jié)果見(jiàn)圖2。

由圖2可知，偏最小二乘模型中這3種向量的DW統(tǒng)計(jì)量均隨著隱變量數(shù)的增加而增加，這反映了噪聲在偏最小二乘模型中的傳遞和累積現(xiàn)象，與“1.2”理論推導(dǎo)的結(jié)果相符合。對(duì)于無(wú)噪光譜所建立的PLS模型，在1～41個(gè)隱變量范圍內(nèi)，r、p和βpls的DW統(tǒng)計(jì)量遠(yuǎn)小于有噪光譜建模得到的各中間向量的DW值。r和p的DW統(tǒng)計(jì)量在第42個(gè)隱變量處突然增至2左右，說(shuō)明此時(shí)這兩個(gè)向量已無(wú)法包含任何有用的光譜信息，因?yàn)槟M數(shù)據(jù)集的真實(shí)組分?jǐn)?shù)是41。然而βpls的DW在第42個(gè)隱變量處卻增加不明顯，這是因?yàn)棣聀ls是由全部權(quán)重向量r的線(xiàn)性組合，其DW的大小應(yīng)該由這些權(quán)重向量的噪聲程度共同決定。對(duì)于包含了儀器噪聲的光譜所建立的PLS模型來(lái)說(shuō)，在前25個(gè)隱變量下r、p和βpls的DW與無(wú)噪光譜PLS模型相比其對(duì)應(yīng)的向量DW相差不大，但在第25個(gè)隱變量之后它們的差異顯著增加，這是因?yàn)樵肼暤膫鬟f和累積的顯著程度會(huì)隨著隱變量數(shù)的增加而增加。

圖2 模擬數(shù)據(jù)集在隱變量1～42下建立的偏最小二乘模型的r(A)、p(B)和βpls(C)的DW統(tǒng)計(jì)量Fig.2 DW values of r(A),p(B) and βpls(C) from simulated dataset with different noise levels at nLVs from 1 to 42

3.2 噪聲傳遞對(duì)模型預(yù)測(cè)能力的影響

圖3 模擬數(shù)據(jù)集中不同噪聲水平光譜所建立 模型的預(yù)測(cè)誤差示意圖Fig.3 The prediction error of PLS models built by simulated dataset with different noise levels

采用上述模擬數(shù)據(jù)驗(yàn)證噪聲傳遞對(duì)PLS模型預(yù)測(cè)能力的影響，結(jié)果見(jiàn)圖3和表1。Mp為無(wú)噪光譜所建立的偏最小二乘模型的RSS，隨著隱變量的增加，模型復(fù)雜度不斷增加，Mp則不斷降低，直至隱變量數(shù)量和模型的真實(shí)組分?jǐn)?shù)(41)達(dá)到一致，此時(shí)模型的復(fù)雜度最優(yōu)。使用無(wú)噪光譜建立的PLS模型的預(yù)測(cè)誤差并不包含Np項(xiàng)，因?yàn)镹p是由儀器噪聲產(chǎn)生。由圖3可知，Np的大小受兩方面影響，一是其隨著PLS隱變量數(shù)量的增加而增加，二是其隨著噪聲水平的增加而增加，這與“1.2”的理論分析一致。由圖4和表1可知，Np的存在和增大會(huì)使模型的預(yù)測(cè)能力下降。

表1 模擬數(shù)據(jù)集的交互檢驗(yàn)和RMSEP結(jié)果Table 1 Results of cross-validation and RMSEP to simulated dataset

圖4和表1描述了模擬數(shù)據(jù)集交互檢驗(yàn)和RMSEP的結(jié)果。對(duì)于無(wú)噪光譜建立的模型，最優(yōu)隱變量數(shù)量和期望的情況一致(均為41)，且在該隱變量下的RMSECV和RMSEP均很接近于0(1.190 2×10-6和2.774 5×10-6)，此時(shí)該模型是這套模擬數(shù)據(jù)集的真實(shí)模型。隨著光譜的噪聲水平增加，模型的RMSECV和RMSEP相應(yīng)增加，但更重要的是模型的最優(yōu)隱變量也不斷發(fā)生變化(例如：對(duì)于噪聲程度為0.3%和1%的光譜來(lái)說(shuō)，最優(yōu)隱變量分別為33和31)。此外，不同噪聲水平下所建立模型的RMSECV和RMSEP的差異隨著隱變量的增加而增加，這證明了“1.2”的理論推導(dǎo)，即噪聲會(huì)在PLS模型中不斷傳遞和累積。由圖4還能發(fā)現(xiàn)，無(wú)噪光譜建立的偏最小二乘模型的預(yù)測(cè)誤差在隱變量超過(guò)41后會(huì)突然顯著增加，這意味著模型的過(guò)擬合。

圖4 模擬數(shù)據(jù)集包含不同水平噪聲的光譜所建立模型的RMSECV(A)和RMSEP(B)Fig.4 RMSECV(A) and RMSEP(B) values at nLVs from 1 to 42 for simulated dataset

圖5 第一個(gè)隱變量下PLS和PoLiSh模型的回歸系數(shù) 及二者差值Fig.5 Regression coefficients of PLS and PoLiSh and their differences

3.3 去噪算法對(duì)噪聲傳遞的驗(yàn)證

采用PoLiSh算法對(duì)人為添加1%噪聲的模擬數(shù)據(jù)進(jìn)行噪聲傳遞的驗(yàn)證。PoLiSh算法的原理為，在使用NIPALS計(jì)算每一個(gè)隱變量下權(quán)重向量時(shí)，對(duì)其使用Savizky-Golay平滑[17]以消除噪聲的積累。本文中Savizky-Golay平滑窗口大小為21，擬合階數(shù)為2，使用PoLiSh所建立模型的最優(yōu)隱變量數(shù)為35，相應(yīng)RMSEP值為0.076 0。與表1中引入了1%噪聲的光譜所建立的PLS模型相比，RMSEP降低且最優(yōu)隱變量的數(shù)目增加，說(shuō)明PoLiSh一定程度上減少了噪聲的影響。圖5反映了第一個(gè)隱變量的部分PLS模型和PoLiSh模型，可以明顯看出使用包含儀器噪聲的光譜建立的PLS模型從第一個(gè)隱變量開(kāi)始包含噪聲，而PoLiSh算法建立的模型可在一定程度上減弱噪聲的影響。

4 結(jié) 論

本文闡述了儀器噪聲如何從第一個(gè)隱變量的計(jì)算開(kāi)始被引入偏最小二乘模型中并通過(guò)隨后的隱變量計(jì)算被不斷傳遞和累積，通過(guò)對(duì)偏最小二乘的計(jì)算過(guò)程相關(guān)公式的推導(dǎo)和對(duì)一套模擬光譜數(shù)據(jù)集的詳細(xì)研究，總結(jié)并論證了與無(wú)噪的光譜相比，噪聲的引入將使得偏最小二乘模型的預(yù)測(cè)能力變差且模型最優(yōu)隱變量數(shù)減少，且光譜中包含的噪聲水平越高，則其對(duì)模型預(yù)測(cè)能力和最優(yōu)隱變量的改變?cè)酱蟮慕Y(jié)論。