儀器噪聲在近紅外光譜偏最小二乘模型內的傳播效應

張若秋，杜一平

(華東理工大學 上海市功能性材料化學重點實驗室 化學與分子工程學院，上海 200237)

近紅外(Near-infrared,NIR)光譜分析是近年來最引人關注的分析技術之一，在農業、制藥、化工等眾多領域得到廣泛應用[1-2]。目前NIR技術的發展重點是應用，簡便、快速、無損、準確等優點正是其受重視并被廣泛應用于國民經濟發展各個領域的本質原因。但在使用NIR技術的基層單位，技術隊伍的研究能力往往不如研發單位，更不及高校和專門的科研機構。而NIR技術使用和推廣的難點是近紅外光譜分析模型的建立需要一些專門的技術，尤其是依賴于化學計量學。因此，模型的建立和維護是近紅外光譜技術發展的一個突出難點和瓶頸。

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

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

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

X=X0+E

(1)

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

1.1 偏最小二乘算法實現

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

Fora=1,…,A;

a=1:對S0進行奇異值分解;

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

計算得到權重向量r=第一個左奇異向量;

計算得到得分向量t=X0r;

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

End

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

1.2 儀器噪聲傳遞

當PLS模型的隱變量數量(nLVs)為k個時，在不同隱變量下所形成的空間和由Krylov序列形成的空間滿足如下等式[5]：

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

(2)

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

(3)

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

(4)

(5)

因此，通過上述公式發現偏最小二乘模型將從第一個隱變量開始發生偏離真模型的現象，噪聲同時被引入到模型中，并在隨后隱變量的計算過程中被不斷傳遞和累積：

(6)

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

(7)

βpls=RkRkTS0

(8)

上述噪聲傳遞過程同時也會對偏最小二乘模型的預測能力產生影響。預測殘差可表示為：

(9)

(10)

(11)

(12)

圖1 噪聲傳遞和偏最小二乘預測誤差關系示意圖Fig.1 The scheme of the relationship between noise propagation and prediction error

噪聲傳遞和偏最小二乘的預測誤差關系如圖1所示，不同噪聲水平下的最優隱變量數也將發生改變，噪聲水平越高，最優隱變量的數值越小。這說明噪聲傳遞不僅影響偏最小二乘模型的預測能力，還會給模型的選擇帶來影響。

2 數據集與計算方法

2.1 數據集模擬

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

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

2.2 指標與計算軟件評價

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

(13)

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

為準確量化偏最小二乘模型的噪聲水平，采用Durbin-Watson(DW)[12-16]統計量評價噪聲對模型的影響程度。DW統計量由下式計算得到：

(14)

xi和xi-1是p維向量中的兩個連續元素。如果整個向量中連續元素間的相關性很弱，例如該向量中每個元素均為隨機變量，則DW統計量的大小將收斂于2。一個擁有至少100個元素的隨機向量的DW統計量的95%置信區間在1.7和2.3之間。

3 結果與討論

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

首先使用DW統計量對3種不同噪聲水平的光譜所建立的偏最小二乘模型的權重向量r、載荷向量p和回歸系數向量βpls進行噪聲程度的估計，其結果見圖2。

由圖2可知，偏最小二乘模型中這3種向量的DW統計量均隨著隱變量數的增加而增加，這反映了噪聲在偏最小二乘模型中的傳遞和累積現象，與“1.2”理論推導的結果相符合。對于無噪光譜所建立的PLS模型，在1～41個隱變量范圍內，r、p和βpls的DW統計量遠小于有噪光譜建模得到的各中間向量的DW值。r和p的DW統計量在第42個隱變量處突然增至2左右，說明此時這兩個向量已無法包含任何有用的光譜信息，因為模擬數據集的真實組分數是41。然而βpls的DW在第42個隱變量處卻增加不明顯，這是因為βpls是由全部權重向量r的線性組合，其DW的大小應該由這些權重向量的噪聲程度共同決定。對于包含了儀器噪聲的光譜所建立的PLS模型來說，在前25個隱變量下r、p和βpls的DW與無噪光譜PLS模型相比其對應的向量DW相差不大，但在第25個隱變量之后它們的差異顯著增加，這是因為噪聲的傳遞和累積的顯著程度會隨著隱變量數的增加而增加。

圖2 模擬數據集在隱變量1～42下建立的偏最小二乘模型的r(A)、p(B)和βpls(C)的DW統計量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 噪聲傳遞對模型預測能力的影響

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

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

表1 模擬數據集的交互檢驗和RMSEP結果Table 1 Results of cross-validation and RMSEP to simulated dataset

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

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

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

3.3 去噪算法對噪聲傳遞的驗證

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

4 結 論

本文闡述了儀器噪聲如何從第一個隱變量的計算開始被引入偏最小二乘模型中并通過隨后的隱變量計算被不斷傳遞和累積，通過對偏最小二乘的計算過程相關公式的推導和對一套模擬光譜數據集的詳細研究，總結并論證了與無噪的光譜相比，噪聲的引入將使得偏最小二乘模型的預測能力變差且模型最優隱變量數減少，且光譜中包含的噪聲水平越高，則其對模型預測能力和最優隱變量的改變越大的結論。