【摘要】 針對高維小樣本光譜數據所顯現的函數型數據（Functional data）特性、與性質參數的非線性關系及變量間存有的嚴重共線性，采用了樣條變換集成罰函數偏最小二乘回歸新技術。它首先以三次B基樣條變換實現非線性光譜數據的線性化重構，隨后將重構的新光譜矩陣交由罰函數偏最小二乘法（Penalized PLS）構建其與性質參變量間的校正模型，其中罰函數中的光滑因子由交叉驗證優化確定以調控模型的擬合精度。最后，通過小麥樣品水分含量的近紅外光譜定量分析，結果顯示該技術光譜數據重構穩健，去噪明顯，并有效解決高維小樣本的過擬合和變量間的共線性，而預測集的均方根誤差（RMSEP）為0.1808%，方法的非線性校正模型預測能力得到了明顯提高。

【關鍵詞】 樣條函數， 偏最小二乘， 粗糙懲罰， 近紅外光譜， 定量分析， 小麥

1 引言 現代光譜以其分析速度快、重現性好、成本低、不消耗樣品、易于實現在線分析等特點而得到廣泛應用。而光譜化學計量學是近代紅外光譜分析技術的重要組成部分，它通過多變量校正技術來進行數據（樣本光譜和其性質參數）處理，以獲得準確的分析結果[1，2]。考慮到近紅外光譜數據通常呈多變量、強相關性，并與樣品性質參變量間的非線性關系，適宜選用非線性偏最小二乘法（Nonlinear PLS，NLPLS）。目前，NLPLS實現方式有3種：一是基于樣本矩陣的非線性變換，即在建模自變量中引入某些原始變量的非線性項，如二次項、交叉項等[3]；二是將建模變量投影到低維的曲線或曲面上得到非線性特征向量，再建立輸入輸出特征向量間的非線性關系[4]，但該方法計算復雜，建模受初值影響大；三是保留PLS的線性外部模型，而內部模型采用多項式、樣條函數、模糊規則、神經網絡、支持向量機等非線性形式[5～9]，該方法缺乏對建模物理變量的直觀解釋能力。 鑒于光譜變量與性質參變量間的具體非線性依存關系不明確，及樣本個體光譜數據顯現為波長變量的函數型數據（Functional data）特性[10]，本研究采用樣條（Spline）變換集成罰函數偏最小二乘（Penalized PLS）回歸新技術，記為SplinePPLS方法。首先利用樣條基函數將光譜自變量與性質因變量之間的未知非線性關系按照各維自變量與因變量的擬線性關系相加展開[11]。由于樣條函數分段擬合、可按需要裁剪以適應任意曲線連續變化的特點，使光譜的重構函數曲線適應光譜數據局部敏感特性的同時保持了函數的光滑性和連續性，從而可削減原始數據中的噪聲。隨后，考慮到光譜矩陣經樣條變換后變量維數顯著增加，將重構的新光譜矩陣交由罰函數偏最小二乘法構建其與性質參變量間的定量線性校正模型，其中基于轉換權向量二階導數的罰函數用以調控模型的擬合精度。為考察SplinePPLS方法的有效性及性能，對小麥近紅外光譜數據進行了研究。

2 SplinePPLS方法的構建

2.1 B基樣條曲線 設變量λ與x滿足如下隨機模型：x=s(λ)＋ε，Ｅ（ε）＝０，Ｖar(ε)＝σ２（１）若λ在區間［a，b］上的一個Ｍ段劃分π∶a=ξ０＜ξ１＜…＜ξＭ＝b， 則式(1)中s(λ)的三次B基樣條逼近曲線方程[11]為：s(λ)=∑M+2l=0clΩ3λ－ξl－1h， a≤λ≤b（2）式中Ω3λ－ξl－1h=13！h3∑4k=0(－１)k4

k(λ－ξl+k－３)３＋，是以ξl+k－３（k=0，1，2，3，4）為內控節點、步長為h的三次B基樣條函數，它與x呈線性關系。對于分點ξl－１及其內控節點ξl+k－３位于劃分的兩側，本研究取 ξ－３＝ξ－２＝ξ－１＝ξ０和ξＭ＝ξＭ＋１＝ξＭ＋２＝ξＭ＋３［１２］。 取觀測位置λ１， λ2， λ…， λp相應數據點x1， x2， …， xp與它們在樣條曲線上插值映射點的距離平方和最小為目標函數，優化求取各基函數的線性加權系數cl(l=0， 1， …， M＋２)，即共有Ｋ＝Ｍ＋３個基函數，從而完成該序列數據形如式（2）的B基樣條逼近曲線構造。

2.2 罰函數偏最小二乘方法（Penalized PLS） 數據點xj(j=1，2，…，p)在B基樣條映射的線性空間中的插值映射點zj的分坐標定義為：zj，0=Ω３（λj－ξ－１h）， zj，l=Ω3(λj－ξ0h， …，，zj，M+2=Ω3(λj－ξＭ＋１h)（3）現將樣本個體自變量x在p維變量空間中的取值，即x=(x1，x2，…，xp)T，計算其所有插值映射點zj各分坐標并加以組合，即可得到x的B基插值映射點矢量 z，即為z=(z1，0， z1，1， …， z1，Ｍ＋２， z２，0， z２，1，…， z２，Ｍ＋２， …， zp，0， zp，1， …，zp，M+2)T（4）

由于每一維變量xj有K個映射分坐標，故矢量z的空間維數將為p×Ｋ。 對于樣本容量n的自變量矩陣Xn×p，欲構建其與性質矢量y間的非線性校正模型，則先實施X的每一樣本個體xi(i=1，2，…，n)的B基樣條變換（各樣本個體選取基函數個數應相同，即K1=K2=…＝Ｋn=K），得到映射樣本矩陣Zn×（p×K），而其與y已演變為擬線性關系。再實施Z與y的線性PLS算法。 鑒于Z較X的變量維數顯著增加，變量間的相關性更為嚴重，本研究采用罰函數偏最小二乘（Penalized PLS，PPLS）方法[10]構建 Z與y間的校正模型。PPLS方法的目標函數為

arg maxwwTZTyyTZwwTw+P(w)(5)

上式分母中基于轉換權向量w的罰函數P(w)=wT(ΔφK2)w， 其中K2=（DK－１ＤＫ）ＴＤＫ－１ＤＫ，而DK為(K－１)×Ｋ維的w一階差分算子矩陣。另外， Δφ＝diag(φ１，φ２，…，φp)為各初始自變量的光滑因子對角陣。Pw實為二階導數罰函數，意在增強w平滑其特征向量t=Zw的能力，以提高模型的穩健性。 由上述PPLS算法思想可見，它是以放棄擬合精度為代價尋求預測性能更優的一種改進PLS方法。同時，PPLS對病態數據的耐受性遠強于普通PLS方法。現將基于初始樣本陣X，y及各自變量光滑因子對角陣Δφ的B基樣條變換的PPLS算法（SplinePPLS）步驟[10]歸結如下：（1）實施Ｘ的n個B基樣條逼近曲線的優化構造，并得其映射矩陣Z；（2）令h=1， Zh=Z，并計算P=ΔφK2及M=(Ip+P)－１，其中Ｉp為p維單位陣；（3）計算轉換權向量wh=MZThy， 并規一化wh=wh/‖wh‖；（4）計算特征向量th=Zhwh，并規一化 th=th/‖th‖；（5）記Th=［t1，t2，…，th］，計算其正交投影矩陣Qh=Th(TThTh)+TTh， 式中“+”為矩陣廣義逆；（6）計算剩余矩陣Zh+1=Zh－ＱhZh；（7）令h=h+1， 重復步驟（3）～（7），直至由交叉驗證(Cross validation )法[13] 確定所需提取的最優成分數h后，將得到轉換權矩陣W=[w1，w2，…，wh]， 進而可計算Ｚ與y間的線性回歸系數β＝（Ｗ（ＷＴＺＴＺＷ）－１ＷＴＺＴy）zy。

3 SplinePPLS方法為小麥近紅外光譜重構及定量分析

3.1 樣本數據說明 小麥樣品數據取自文獻[14]，自變量取其近紅外光譜在波長1100～2500 nm、掃描分辨率為2 nm的若干波長處的吸光度值log(1/R)， 即維數p=701，其中R為樣本的反射率，樣本容量n=100，譜圖如圖1所示。性質參變量為小麥水分質量百分含量，數值范圍在12.45%～17.36%。從原始數據集中隨機劃出80個構成訓練集用于光譜B基樣條變換的優化確定及后繼校正模型建立，其余20個組成獨立測試集，用于檢驗B基樣條的光譜插值重構能力及模型的預測性能。

Fig.1 NIR diffuse reflectance spectra of wheat samples

3.2 實驗方式與性能評價指標 為檢驗SplinePPLS方法的性能，先將訓練樣本Xntrain×p以“變量‘留一’交叉驗證”[13]選定K個樣條變換的B基函數，其優化評定指標為式(6)中的RMSECVspline；再以“樣本個體‘留一’交叉驗證”選定各變量的光滑因子φ1，φ2，…，φp及校正模型所需PLS最優成分數h，它們的優化評定指標則為式(6)中的RMSECVppls。而光譜B基樣條變換的插值重構性能及校正模型的預測性能，則交由測試樣本Xntest×p計算，它們的評價指標分別為式(7)中的RMSEPspline和RMSEPppls。RMSECVspline=∑pj=1∑ntraini=1(xij－ij)2/(ntrain×p)， RMSECVppls=∑ntraini=1(yi－i)2/ntrain（6）

RMSEＰspline=∑pj=1∑ntesti=1(xij－ij)2/(ntest×p)， RMSEＰppls=∑ntesti=1(yi－i)2/ntest（7）式中xij和ij分別為第i樣本個體、第j波長變量下吸光度的實驗值和B基樣條曲線的插值；yi和i則分別為第i樣本個體性質參變量的實驗測試值和模型預報值。

3.3 結果與分析 3.3.1 光譜數據的B基樣條變換 現基于光譜陣Xntest×p，構造它們的３次B基樣條逼近曲線，據此完成測試光譜陣Xntest×p的三次B基樣條變換。即，先將光譜波長1100～2500 nm作預設的M段劃分，再實施Xntrain×p的“變量‘留一’交叉驗證”實驗，即依次留用一個測量位置(變量)λj的觀測數據矢量x·j=[x1j，x2j，…，xntrainj]T作內部驗證，而剩余p－１列觀測數據用于樣條逼近曲線的最小二乘擬合，可得ntrain個樣條曲線各K=M+3個基函數的線性加權系數ci，j(l=0，1，…，M+2，i=1，2，…，ntrain)， 并據此實現x·j的插值估計·j。最后，將p輪循環得到的所有·j(j=0，1，…，p)代入式(6)，即可計算該M取值下的RMSECVspline。改變M取值，并依據對插值精度RMSECVspline的要求，即可選定所需的基函數個數K。 圖2 B基樣條曲線的基函數個數優化確定及第1訓練樣本個體的光譜重構結果（略）

Fig.2 Selection of basic functions number and the rebuilding spectrum curve for the first sample data

a. Rootmean squared error at different numbers of basic functions; b. 28 Basic functions and its rebuilding spectrum curve.

圖2a顯示了RMSECVspline與K的相關關系，隨著K的增多，RMSECVspline總體呈下降趨勢，在K=28時，RMSECVspline已很小，而其后的RMSECVspline值下降有限。考慮到參數K取值越大，經B基樣條變換后的擬線性變量個數將越多，且變量間將出現更多復共線性， B基樣條逼近曲線易出現對訓練數據的過擬合，而對包含于數據中的噪音削減不夠；若K取值太小，B基樣條逼近曲線對數據的插值能力將下降。因此，本研究選定K=28。圖2b顯示了這28個基函數及由它們所重構的第1訓練樣本個體光譜數據的三次B基樣條逼近曲線。由圖2b可見，該樣條曲線實現了對光譜數據較高精度的插值擬合。另外，將這28個基函數用于Xntest×p的插值重構，其精度指標RMSEＰspline，表明B基樣條具有極強的插值重構能力。

3.3.2 SplinePPLS方法中參數的選擇 對于模型而言，預報性能最為重要。影響SplinePPLS模型預報精度的主要因素有：B基樣條變換基函數個數K、PPLS模型中各自變量光滑因子φ1，φ2，…，φp及PLS最優成分數h等。其中，K的優選過程見3.3.1節，并將其取值為K=28。而表1則為PPLS方法施于樣本陣｛Ｘ，y｝優化選擇φ1， φ2， …， φp及h的過程結果。其中，為減少搜索空間維數，將各光滑因子簡單取值相同φ=φ1=φ2=…，＝φp。

表1 PPLS方法中參數的優化選擇（略）

Table 1 Parameter selection of the penalized PLS method

從表1可見，模型性能指標ＲＭＳＥＣＶppls在參數φ和h的二維搜索格子點φ＝１０００００， h=5位置達最小，由此選定它們為PPLS模型相應參數的最優取值。同時還發現，該位置ＲＭＳＥＣＶppls指標小于φ＝０， h=5的PLS方法的結果。由2.2節PPLS方法的目標函數式(5)可知，φ取值的大小將影響w平滑其特征向量t=Zw的程度，并進一步作用于PPLS方法的模型系數β＝（Ｗ（ＷＴＸＴＸＷ）－１ＷＴＴＴy）x，y。圖3即為PLS和PPLS方法施于樣本陣｛Ｘ，y｝的β結果比較。其中PPLS模型系數β曲線受到變量光滑因子φ的粗造懲罰而得到了平滑。這樣，它既可最大限度保證分析信號不失真，又能進一步削除噪音。

圖3 PLS和PPLS模型回歸系數比較（略）

Fig.3 Comparison of model coefficients for the PLS and Penalized PLS method

3.3.3 模型精度比較和分析 表2為多元線性回歸MLR（Multiple linear regression）， PLS， PPLS， SplinePLS及SplinePPLS的建模結果，各方法的參數最優值由交叉驗證法[13]確定。另外，SplinePLS和SplinePPLS方法中Ｋ＝２８。先分析校正方法對模型精度的影響，SplinePLS和SplinePPLS屬NLPLS校正方法，它們的RMSECV和RMSEP分別小于同屬線性校正的PLS和PPLS方法，但它們提取的PLS最優成分數h均多于后兩種線性方法。由此說明，B基樣條變換在一定程度上實現了光譜數據與樣品性質變量間的內在非線性關系，但需通過增加PLS成分將包含在擬線性化變量項中的非線性信息帶入校正模型。而MLR方法未進行光譜數據的噪音削減，以及消除變量間的復共線性，RMSECV和MSSEP值均顯著高于其它4種方法。PPLS和SplinePPLS方法的RMSECV和RMSEP值均分別小于未進行粗糙懲罰的PLS和SplinePLS方法的。由此說明，在用于PLS特征向量提取的目標函數中集成轉換權向量二階導數的罰函數，可平滑特征向量和進一步削減噪音，從而使模型的預測能力和穩定性得以提高。

表2 5種不同校正模型的性能比較（略）

Table 2 Comparison of model performance for five methods