熊 芩，張若秋，李 輝，陳萬超，杜一平*

(1.華東理工大學(xué) 化學(xué)與分子工程學(xué)院 上海市功能性材料化學(xué)重點(diǎn)實(shí)驗(yàn)室，上海 200237；2.德宏師范高等專科學(xué)校 理工系，云南 德宏 678400)

近年來，近紅外光譜(Near infrared spectrum，NIR)技術(shù)得到廣泛應(yīng)用，但由于NIR數(shù)據(jù)通常存在嚴(yán)重的共線性問題，對其進(jìn)行降維和變量選擇的研究一直是該領(lǐng)域的研究熱點(diǎn)。經(jīng)典的降維方法，如主成分回歸(PCR)和偏最小二乘回歸(PLS)，雖能實(shí)現(xiàn)維度的大幅度降低，但模型中仍包括所有變量(波長)的信息，模型缺乏可解釋性。有效解決變量共線性的另一方式為變量選擇，近年提出了很多變量選擇算法，如：無信息變量消除法(UVE)、蒙特卡洛無信息變量消除法(MCUVE)[1]、移動窗口偏最小二乘法(MWPLS)、競爭自適應(yīng)重加權(quán)采樣法(CARS)[2]、穩(wěn)定性競爭自適應(yīng)重加權(quán)采樣法(SCARS)[3]、遺傳偏最小二乘法(GA-PLS)等，這些方法已被廣泛應(yīng)用于光譜變量的選擇中。

最小角回歸(Least angle regression，LAR)算法由Efron等[4]提出，常用于回歸分析和高維數(shù)據(jù)的特征選擇[5-7]，在光譜建模中應(yīng)用較少。該算法在前向梯度算法的基礎(chǔ)上，逐步前進(jìn)、步長適中，無需進(jìn)行大量迭代運(yùn)算，大大降低了計(jì)算復(fù)雜度，提高了運(yùn)行速度。Liu等[8]利用LAR對臍橙的近紅外數(shù)據(jù)集進(jìn)行了回歸預(yù)測，相較于PLS模型具有更好的預(yù)測性能和解釋性，同時運(yùn)行速度快于最小二乘支持向量機(jī)(LS-SVM)。顏勝科等[9]通過LAR將系數(shù)不為0的變量作為初步篩選的變量，進(jìn)一步結(jié)合GA-PLS再次篩選波長，從光譜數(shù)據(jù)中提取出信息變量，預(yù)測精度優(yōu)于全光譜模型。

建立NIR模型通常需將樣本集劃分為校正集、交互檢驗(yàn)集、驗(yàn)證集等，而不同的劃分會得到完全不同的模型。針對這一問題，本課題組[10]基于模型集群分析(MPA)[11]提出采樣誤差分布分析(Sampling error profile analysis，SEPA)算法。SEPA通過蒙特卡洛采樣(MCS)獲得多個(比如800或2 000)數(shù)據(jù)劃分，建立多個子模型，對子模型的誤差分布進(jìn)行統(tǒng)計(jì)分析，從而解決了單次劃分導(dǎo)致的偶然性問題，并成功應(yīng)用于PLS模型優(yōu)化和評價[10]以及模型傳遞[12]。Zhang等[13]還將SEPA與最小絕對收縮和選擇算子(LASSO)[14]相結(jié)合，提出了SEPA-LASSO的變量選擇方法，相較于SCARS、MCUVE、MWPLS等方法，SEPA-LASSO選擇更少的變量，獲得了更優(yōu)異的預(yù)測性能。本文采用SEPA的思路，基于LAR回歸系數(shù)提出了一種SEPA-LAR的變量逐步篩選方法，用于建立穩(wěn)健的近紅外分析模型。通過3套近紅外數(shù)據(jù)進(jìn)行建模分析，并與PLS、MCUVE、MWPLS、CARS等方法進(jìn)行了比較。

1 理論與算法

1.1 最小角回歸

矩陣Xn×p表示NIR光譜矩陣，n為樣本數(shù)，p為變量數(shù)，yn×1表示指標(biāo)向量。

最小角回歸算法與經(jīng)典的前向梯度回歸算法(Forward stepwise regression)類似，首先將所有變量回歸系數(shù)的初值賦為0，然后找出與因變量y相關(guān)性最大的變量，假設(shè)為x1。不同于前向梯度算法，LAR將直接沿x1的方向?qū)逼近，直至下一個變量x2出現(xiàn)，使得x1、x2與當(dāng)前殘差具有相同的相關(guān)性，再沿x1與x2的角平分線方向找到第3個變量，以此類推。LAR算法的每一步路徑均保持當(dāng)前殘差與所有入選變量的相關(guān)性相同。

圖1 變量數(shù)為2時LAR算法路徑示意圖Fig.1 The path of LAR algorithm in case of 2 variables

LAR的具體運(yùn)算步驟可參考文獻(xiàn)[4]。

1.2 采樣誤差分布分析(SEPA)

以交互檢驗(yàn)的計(jì)算過程為例，按3個步驟對SEPA(SEPA-CV)進(jìn)行說明：設(shè)采樣次數(shù)為N，首先通過MCS將樣本集合劃分為校正集和交互驗(yàn)證集，獲得N個不同的樣本集合劃分；然后分別對各樣本集合建立校正模型，并計(jì)算其交互驗(yàn)證預(yù)測均方根誤差(RMSECV)，對所有得到的N個RMSECV誤差數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析；最終綜合考慮RMSECV中位數(shù)及其標(biāo)準(zhǔn)偏差(STD)，以及誤差分布輪廓(傾向于正態(tài)分布的優(yōu)先選擇)來選擇隱變量。由于采用隨機(jī)選擇樣本集劃分，并多指標(biāo)綜合考慮誤差分布，因此所選擇的參數(shù)和建立的模型更加穩(wěn)健。

1.3 最小角回歸結(jié)合采樣誤差分布分析(SEPA-LAR)建立穩(wěn)健模型

SEPA-LAR的基本思路是通過N次Monte Carlo采樣得到N個劃分，對每個劃分進(jìn)行LAR計(jì)算獲得不同的變量選擇集合以及各變量集合下的LAR回歸系數(shù)，對同一變量數(shù)集合，N個劃分將得到N個變量選擇結(jié)果。對相同變量數(shù)集合下的回歸系數(shù)進(jìn)行統(tǒng)計(jì)分析，統(tǒng)計(jì)各變量回歸系數(shù)的絕對值之和，其值越大說明該變量越重要，因此選擇回歸系數(shù)絕對值之和最大的若干個變量作為候選的變量選擇集合。利用這些變量集合建立PLS模型進(jìn)行交互檢驗(yàn)，以確定最佳變量集合進(jìn)行建模。

SEPA-LAR步驟如下：

①通過N次MCS，得到N個不同的樣本子集(校正集和驗(yàn)證集)，其中校正集的樣本數(shù)均為n。對每個子集進(jìn)行LAR計(jì)算，得到n-1個變量選擇集合(其變量數(shù)i為1～n-1)。N個子集的所有變量選擇集合共計(jì)N×(n-1)；

②對每一變量數(shù)i的所有變量選擇集合(N個)，統(tǒng)計(jì)每個變量的回歸系數(shù)絕對值之和，將該值進(jìn)行降序排列，選擇前i個變量?？色@得i從1到n-1的n-1個變量子集；

③采用SEPA和PLS進(jìn)行交互驗(yàn)證，計(jì)算各變量子集下的RMSECV，以RMSECV中位數(shù)和STD為依據(jù)選擇最佳變量子集；

④在最佳變量子集下，利用SEPA-PLS選擇最佳隱變量數(shù)；

⑤用最終選定的最佳變量子集以及隱變量數(shù)建立PLS模型，并用獨(dú)立驗(yàn)證集對模型進(jìn)行評價。

1.4 計(jì)算方法

所有計(jì)算在MATLAB(Version 2010a，The MathWork，USA)環(huán)境中自編程序完成，所有數(shù)據(jù)均進(jìn)行均值中心化處理，以消除光譜的絕對吸收值，突出樣品間差異。對比方法PLS、MWPLS、MCUVE以及CARS的隱變量數(shù)目由十折交互檢驗(yàn)確定，其中，MWPLS窗口大小設(shè)為15，MCUVE迭代次數(shù)設(shè)為100，CARS迭代次數(shù)設(shè)為100。

對整個數(shù)據(jù)用Kennard-Stone(K-S)方法[15]劃分出一定數(shù)目的樣本組成獨(dú)立驗(yàn)證集，用以評價模型。其余樣本進(jìn)行SEPA-LAR計(jì)算和PLS建模。在SEPA-LAR算法中，MCS采樣次數(shù)N均設(shè)為2 000，隨機(jī)抽取的校正集和驗(yàn)證集比例設(shè)為60∶40應(yīng)用在“1.3”模型的步驟①和③。

2 數(shù)據(jù)集

2.1 玉米數(shù)據(jù)集

數(shù)據(jù)集下載自http://www.eigenvector.com/data/Corn/index.html，包含80個玉米樣本的近紅外光譜數(shù)據(jù)，波長范圍1 100～2 498 nm，分辨率為2 nm，共700個波長點(diǎn)。本文采用m5儀器測得的光譜數(shù)據(jù)和玉米樣品的濕度數(shù)據(jù)，濕度在9.377%～10.993%之間，獨(dú)立驗(yàn)證集有20個樣本。

2.2 柴油數(shù)據(jù)集

數(shù)據(jù)下載自http://www.eigenvector.com/data/SWRI/index.html，為柴油樣本的近紅外光譜數(shù)據(jù)，波長范圍750～1 550 nm，分辨率2 nm，共401個波長點(diǎn)。本文采用密度數(shù)據(jù)，密度為0.781 8～0.872 8 g/mL。共395個樣本，獨(dú)立驗(yàn)證集195個樣本。

2.3 奶酪數(shù)據(jù)集

奶酪樣本購于上海某大型連鎖超市，產(chǎn)自8個不同品牌廠家，共107個樣本。用微型近紅外光譜儀(INSION近紅外檢測器，德國)采集光譜，波長范圍909～1 860 nm，分辨率為8 nm，共116個波長點(diǎn)。奶酪中脂肪含量采用GB 5009.6-2016堿水解法[16]進(jìn)行測定，含量范圍為8.21～26.65 g/100 g，平均含量為19.51 g/100 g，獨(dú)立驗(yàn)證集37個樣本。

3 結(jié)果與討論

3.1 玉米數(shù)據(jù)集預(yù)測

利用K-S方法選擇的60個樣本進(jìn)行SEPA-LAR以及PLS建模，其余20個樣本為獨(dú)立驗(yàn)證集。

圖2 基于SEPA的不同變量數(shù)的交互檢驗(yàn)結(jié)果Fig.2 Variation of RMSECVs with the number of variables based on SEPA

圖3 玉米數(shù)據(jù)中SEPA-LAR、MCUVE、MWPLS以及CARS選擇的變量Fig.3 The variables selected by SEPA-LAR,MCUVE,MWPLS and CARS on corn dataset

3.1.1變量選擇采用MCS方法從60個樣本中隨機(jī)抽取24個樣本作為驗(yàn)證集，剩余36個樣本組成校正集。經(jīng)LAR計(jì)算得到35個變量子集(變量數(shù)i=1,2,…,35)。不同變量選擇下的交互檢驗(yàn)結(jié)果見圖2，其中橫坐標(biāo)表示變量數(shù)。由圖2可見，當(dāng)選擇的變量數(shù)低于20時，RMSECV隨著變量數(shù)的增加明顯減小，當(dāng)變量數(shù)達(dá)到20之后，RMSECV下降趨于平緩，同時STD(圖2中的誤差棒)很小。由此可見，所選的20個變量對于預(yù)測玉米濕度具有重要作用，因此最終選擇20個變量用于PLS建模。當(dāng)變量數(shù)為20時，各變量回歸系數(shù)絕對值之和的統(tǒng)計(jì)結(jié)果表明，大部分的波長變量未被LAR選中，而某些特征變量可被LAR選中。LAR能夠篩選出相關(guān)性大的變量并賦予其回歸系數(shù)，而將其余變量的回歸系數(shù)直接賦值為0。由于集合劃分的不同，不同樣本子集下，LAR并非選中同樣的波長變量，SEPA在建模中充分考慮了集合拆分的影響，最終選擇的特征變量是基于大量重復(fù)采樣后的統(tǒng)計(jì)結(jié)果，使所建模型更合理與可靠。

圖3給出了SEPA-LAR、MCUVE、MWPLS以及CARS選擇的波長變量。1 900 nm與2 100 nm左右屬于水分子O—H的組合頻吸收譜帶，1 410 nm左右屬于水分子O—H的一級倍頻吸收譜帶[2]，該波長區(qū)域?qū)τ谟衩诐穸鹊念A(yù)測具有重要作用。CARS選擇了1 908、2 108 nm兩個特征波長，SEPA-LAR和MCUVE均選擇到該重要吸收譜帶內(nèi)的波長，且兩種方法所選波長基本一致，但SEPA-LAR選中的波長變量更少，其預(yù)測精度更高。MWPLS則選擇了更多的波長變量，未選中包含2 100 nm左右的特征波長。

3.1.2模型的評價利用SEPA-LAR選擇的變量以及SEPA-CV確定隱變量數(shù)，建立PLS模型，同時分別采用全譜以及MWPLS、MCUVE、CARS波長選擇方法篩選出特征變量建立PLS模型，利用獨(dú)立驗(yàn)證集評價上述模型，表1給出了SEPA-LAR及其他4種方法建模結(jié)果的比較。由表1可知，對于獨(dú)立驗(yàn)證集，全譜PLS模型的預(yù)測均方根誤差(RMSEP)為0.015 2%，MWPLS-PLS、MCUVE-PLS、CARS-PLS以及SEPA-LAR的RMSEP分別為0.029 0%、0.002 76%、0.000 314%以及0.001 44%。相較于全譜PLS建模，采用MCUVE、CARS和SEPA-LAR方法對變量進(jìn)行篩選均能提高對玉米濕度的預(yù)測性能。相比MCUVE方法，SEPA-LAR通過篩選出更少的變量數(shù)，在提高建模效率的同時獲得了更高的預(yù)測精度，表明SEPA-LAR對玉米濕度數(shù)據(jù)具有優(yōu)異的預(yù)測性能。在玉米濕度這類變量特征性強(qiáng)的簡單體系中，CARS方法表現(xiàn)優(yōu)異，僅篩選出兩個特征變量，預(yù)測精度優(yōu)于本方法及其他對比方法。但當(dāng)樣品體系復(fù)雜時，CARS選擇的波長子集波動較大，建模穩(wěn)定性變差。而SEPA-LAR通過多指標(biāo)綜合考慮誤差分布，模型更穩(wěn)健。

表1 5種方法對玉米濕度的預(yù)測結(jié)果Table 1 Results of five methods on the corn moisture dataset

nVAR represents the number of selected variables；nLV represents the number of latent variables；RMSEC represents the root-mean-square error of calibration；RMSEP represents the root-mean-square error of prediction

3.2 柴油數(shù)據(jù)集與奶酪數(shù)據(jù)集的預(yù)測

利用SEPA-LAR對柴油數(shù)據(jù)和奶酪數(shù)據(jù)進(jìn)行了變量選擇和模型建立，并與全譜、MWPLS、MCUVE、CARS方法進(jìn)行比較。表2為柴油密度和奶酪脂肪數(shù)據(jù)的建模結(jié)果。相較于其他4種方法，SEPA-LAR表現(xiàn)出更優(yōu)異的建模能力，且選擇的變量最少，模型的預(yù)測誤差RMSEP最小。對于柴油密度數(shù)據(jù)，SEPA-LAR僅選擇了28個波長，RMSEP低至0.001 58 g/mL。將SEPA-LAR應(yīng)用于奶酪脂肪數(shù)據(jù)進(jìn)行建模，所建模型僅含有25個波長，RMSEP為1.13 g/100 g。

表2 5種方法對柴油密度和奶酪脂肪的預(yù)測結(jié)果Table 2 Results of five methods on the diesel density dataset and cheese fat dataset

圖4給出了柴油數(shù)據(jù)和奶酪數(shù)據(jù)的近紅外光譜，以及4種方法所選擇的波長。由圖4A可見，對于柴油數(shù)據(jù)，4種方法選擇的波段近似，所選波長范圍集中在750～850、1 000～1 100、1 200～1 300、1 400～1 550 nm，這些區(qū)域?qū)?yīng)柴油的近紅外主要吸收峰位置。與MCUVE和CARS相比，SEPA-LAR選擇的變量大都被MCUVE和CARS選中，但SEPA-LAR最終選擇的變量更少，模型的預(yù)測性能更優(yōu)異。

由圖4B可見，4種變量選擇方法選擇了大體相同的波長區(qū)域，即900～1 000、1 100～1 350、1 580～1 750 nm。其中SEPA-LAR、MCUVE和CARS選擇的具體波長差別較大，而SEPA-LAR的所選波長幾乎均被MWPLS選中。CARS選擇的波長更少，但其預(yù)測性能不如SEPA-LAR。相較于MCUVE，SEPA-LAR則選擇更少的變量，模型的誤差更低。研究結(jié)果表明，SEPA-LAR可在選擇較少變量參與建模的同時，提高模型的預(yù)測性能。

4 結(jié) 論

本文提出了一種基于最小角回歸算法和SEPA算法的變量逐步篩選方法，用于構(gòu)建穩(wěn)定的近紅外光譜分析模型。以玉米濕度、柴油密度以及奶酪脂肪數(shù)據(jù)集為例，首先基于大量的蒙特卡洛采樣計(jì)算得到的LAR回歸系數(shù)統(tǒng)計(jì)結(jié)果得出各個變量子集，再采用SEPA框架的思路進(jìn)行交互檢驗(yàn)確定最佳變量子集，最后利用得到的最佳變量子集通過SEPA-CV確定PLS模型的隱變量數(shù)，從而建立穩(wěn)健合理的近紅外光譜模型。實(shí)驗(yàn)結(jié)果表明，相較于全譜以及MCUVE、MWPLS方法，SEPA-LAR能有效篩選出更少的特征變量，消除原始變量的共線性問題，所建模型具有優(yōu)異的預(yù)測性能，穩(wěn)定性優(yōu)于CARS方法，相較于全譜PLS模型擁有更好的可解釋性，有效提高了模型的準(zhǔn)確性和穩(wěn)定性。