程介虹，陳爭光,*，衣淑娟

1.黑龍江八一農(nóng)墾大學(xué)信息與電氣工程學(xué)院，黑龍江 大慶 163319 2.黑龍江省水稻生態(tài)育秧裝置及全程機(jī)械化工程技術(shù)中心，黑龍江 大慶 163319

引 言

近紅外光譜區(qū)(800～2 500 nm)的含氫基團(tuán)的倍頻和合頻吸收帶較寬且嚴(yán)重重疊，全譜建模定量分析會(huì)存在多重共線性或無信息變量過多導(dǎo)致模型性能不佳。因此，對全譜進(jìn)行特征波長選擇在一定程度上可以減少數(shù)據(jù)冗余和多重共線性，提高模型的預(yù)測精度和預(yù)測效率。特征波長選擇一直是近紅外光譜分析中的熱點(diǎn)。

常用的近紅外光譜特征波長選擇方法包括：無信息變量消除法[1]、競爭性自適應(yīng)重加權(quán)采樣法[2]、間隔偏最小二乘法[3]、遺傳算法[4]、連續(xù)投影算法[5]、隨機(jī)蛙跳[6]、迭代保留信息變量法[7]等等。所有這些研究表明，使用特征波長代替全譜建?？梢垣@得更好的預(yù)測精度，這說明波長選擇的重要性。

眾所周知，在近紅外光譜的定量分析中，變量間的多重共線性問題是導(dǎo)致模型效果較差的一個(gè)關(guān)鍵問題。這是因?yàn)槿绻麅蓚€(gè)向量之間存在共線性，意味著兩個(gè)向量具有相似的趨勢并且有可能攜帶相似信息，這類變量的存在會(huì)降低模型的性能。變量間共線性可以通過相關(guān)性分析判斷，相關(guān)系數(shù)高于0.8，表明存在多重共線性。所以，以變量間相關(guān)系數(shù)為選擇標(biāo)準(zhǔn)，提出一種以所選變量之間共線性最小化，并且輸入變量對響應(yīng)變量影響最大化為主要目的的波長選擇方法，稱之為最小相關(guān)系數(shù)法(minimal correlation coefficient,MCC)。

我們將該算法應(yīng)用于兩組公開的近紅外光譜數(shù)據(jù)集，并與其他常用波長選擇方法，如：連續(xù)投影算法(successive projections algorithm,SPA)、競爭性自適應(yīng)重加權(quán)采樣法(competitive adaptive reweighted sampling,CARS)、隨機(jī)蛙跳算法(random frog,RF)、迭代保留信息變量法(iteratively retains informative variables,IRIV)進(jìn)行對比，以說明本方法的有效性。

1 算法原理

1.1 符號說明

矩陣用粗斜體大寫字母表示，向量由粗體小寫字母表示，變量(標(biāo)量)用斜體字母表示。光譜數(shù)據(jù)矩陣表示為XN×K，響應(yīng)(濃度)向量表示為yN×1，其中，N為樣本數(shù)，K為波長數(shù)。下標(biāo)變量i,j表示光譜矩陣的第i列、第j列，“^”符號代表預(yù)測值。

1.2 相關(guān)概念

1.2.1 向量間最小相關(guān)系數(shù)

在近紅外光譜的定量分析中，多重共線性問題是導(dǎo)致模型效果較差的一個(gè)關(guān)鍵問題。對于多元線性回歸模型來說，當(dāng)變量之間不存在線性相關(guān)性，即參與建模的變量中的任一變量與其他變量間均不存在線性相關(guān)性時(shí)，模型的預(yù)測性能較優(yōu)。反之，如果參與建模的變量中的任一變量與其他變量間存在線性相關(guān)性，則模型存在多重共線性。

皮爾森相關(guān)系數(shù)(Pearson correlation)是一種衡量向量之間線性相關(guān)性的指標(biāo)，計(jì)算方法如式(1)所示。

(1)

式中，cov(x,y)為向量x和y的協(xié)方差，σx，σy為向量x和y的標(biāo)準(zhǔn)差，E(x),E(y)為向量x和y的期望。

皮爾森相關(guān)系數(shù)ρ取值區(qū)間為[-1,1]，其絕對值是表征兩個(gè)向量相關(guān)性的一個(gè)度量(表1)。如果兩個(gè)向量之間的皮爾森相關(guān)系數(shù)ρ絕對值高于0.8，即可判定存在一定程度的多重共線性；而兩個(gè)向量之間的皮爾森相關(guān)系數(shù)ρ絕對值小于0.4，那么兩個(gè)向量弱相關(guān)。因此，向量之間皮爾森相關(guān)系數(shù)ρ的大小可以作為判斷向量間共線性程度的一個(gè)度量。

表1 變量間相關(guān)系數(shù)所對應(yīng)的相關(guān)強(qiáng)度Table 1 The strength of correlation corresponding to the correlation coefficient between variables

光譜數(shù)據(jù)矩陣XN×K可以看作是一個(gè)包含K個(gè)N維向量的向量組，每個(gè)向量代表了在特定波長下近紅外光譜的吸光度特性。為了后面算法描述方便，在此用運(yùn)算符r(·)表示兩個(gè)向量之間的相關(guān)系數(shù)絕對值，r(xi,xj)與皮爾森相關(guān)系數(shù)ρ的關(guān)系如式(2)所示。r(·)具有特性：r(xi,xj)=r(xj,xi)，r(xi,xi)=1，如果光譜矩陣進(jìn)行了標(biāo)準(zhǔn)化處理，那么，r(xi,xj)=|cov(xi,xj)|。

r(xi,xj)=|ρ(xi,xj)| (i,j=1,2,3,…,K)

(2)

由上述理論可知，如果某一個(gè)波長與其余波長之間的相關(guān)系數(shù)r(·)均較小，則意味著該波長不能由其他波長線性表示，需要作為建模變量保留下來，這個(gè)波長即為關(guān)鍵波長。反之，如果一個(gè)波長與其余波長間的相關(guān)系數(shù)均較大，則意味著該波長與其余波長間存在多重共線性，不能作為建模變量，需要剔除。

那么，如何判斷某一波長與其余所有波長之間的相關(guān)性均較小呢？首先計(jì)算矩陣XN×K中各列向量之間相關(guān)系數(shù)的絕對值，得到相關(guān)系數(shù)矩陣RK×K，其中第i列向量為XN×K中第i個(gè)波長xi與其他波長之間的相關(guān)系數(shù)(圖1)。然后計(jì)算矩陣RK×K中各列向量的平均值r和標(biāo)準(zhǔn)差σ，選取平均值ri和標(biāo)準(zhǔn)差σi均較小的列對應(yīng)的波長作為候選波長。這是因?yàn)?，對于RK×K的第i列向量，其對應(yīng)的ri較小不一定說明該列向量中的所有元素均較小，只有同時(shí)滿足標(biāo)準(zhǔn)差σi較小時(shí)，才能說明RK×K中第i列相關(guān)系數(shù)向量各元素值均較小，且分布較為集中。此時(shí)XN×K中第i個(gè)波長應(yīng)該保留下來作為建模候選波長。

圖1 相關(guān)系數(shù)平均值與標(biāo)準(zhǔn)差Fig.1 Average value and standard deviation of correlation coefficient

1.2.2 最大標(biāo)準(zhǔn)回歸系數(shù)

由于通過MCC方法消除了變量之間的相關(guān)性，因此，可以建立基于MCC方法的波長選擇結(jié)果的線性回歸方程。假設(shè)線性回歸方程的形式如式(3)

(3)

為了使不同變量的效應(yīng)大小具有可比性，通常對回歸系數(shù)βi進(jìn)行無量綱化處理，將其按式(4)轉(zhuǎn)化為一個(gè)無量綱的標(biāo)準(zhǔn)回歸系數(shù)。

(4)

向前選擇法是一種回歸模型的自變量選擇方法，其特點(diǎn)是把候選的自變量逐個(gè)引入回歸方程，故稱向前法。首先對符合最小相關(guān)系數(shù)準(zhǔn)則的波長按如前所述計(jì)算得到標(biāo)準(zhǔn)回歸系數(shù)bi，并按標(biāo)準(zhǔn)回歸系數(shù)絕對值|bi|大小進(jìn)行降序排列，然后按排列順序每次將一個(gè)標(biāo)準(zhǔn)回歸系數(shù)所對應(yīng)的波長引入模型，建立線性回歸模型，計(jì)算每個(gè)模型的均方根誤差，以獲得最小均方根誤差所引入的變量為準(zhǔn)得到一個(gè)最優(yōu)模型。

1.3 算法步驟

根據(jù)上述理論基礎(chǔ)，MCC具體算法步驟如下：

步驟1：利用通過式(2)計(jì)算光譜矩陣XN×K中各列向量之間相關(guān)系數(shù)的絕對值，得到相關(guān)系數(shù)矩陣RK×K；

步驟2：計(jì)算矩陣RK×K中各列數(shù)據(jù)除對角線元素之外的其他元素的平均值ri和標(biāo)準(zhǔn)差σi(圖1)；

步驟3：通過網(wǎng)格尋優(yōu)法(后面詳細(xì)說明)，選擇滿足某一閾值條件的波長組成候選波長集S；

步驟4：以集合S中波長建立的MLR方程[式(3)]，然后通過式(4)計(jì)算標(biāo)準(zhǔn)回歸系數(shù)bi。對集合S中的波長按bi的絕對值降序排序，得到集合S′;

步驟5：對集合S′中的波長通過向前選擇法建模獲得該閾值條件下最小均方根誤差模型；

步驟6：選擇網(wǎng)格中的下一對閾值，重復(fù)步驟3—步驟5，直至完成網(wǎng)格尋優(yōu)。

END：所有閾值條件下最小均方根誤差對應(yīng)的變量集合即為最終所選波長。

1.4 網(wǎng)格尋優(yōu)

算法步驟3中的網(wǎng)格尋優(yōu)法，具體操作表述如下：

設(shè)閾值范圍(rmin,rmax)=(min(ri),max(ri))，(σmin,σmax)=(min(σi),max(σi)),i=1,2,…,K，將閾值范圍(rmin,rmax)和(σmin,σmax)劃分為q×q均勻網(wǎng)格T，如圖2所示。以網(wǎng)格T中每一個(gè)點(diǎn)(trm,tσn)(m,n=1,2,…，q)作為閾值對，從矩陣R中選擇滿足條件(r

圖2 網(wǎng)格T閾值對圖解Fig.2 Graphic analysis of threshold in grid T

注意，當(dāng)m=1或n=1時(shí)，tr1=rmin，tσ1=σmin，因?yàn)閞min和σmin分別是相關(guān)系數(shù)的平均值的最小值和標(biāo)準(zhǔn)差的最小值，此時(shí)滿足條件(r

同樣，當(dāng)m=q或n=q時(shí)，trq=rmax，tσq=σmax時(shí)，則所有的波長都滿足條件(r

因此，取兩個(gè)閾值系數(shù)thr1>1和thr2<1，令(rmin,rmax)=(thr1*min(ri),thr2*max(ri))，(σmin,σmax)=(thr1*min(σi),thr2*max(σi))，在此基礎(chǔ)上生成網(wǎng)格矩陣T。這樣，當(dāng)滿足條件(r

根據(jù)經(jīng)驗(yàn)，閾值系數(shù)thr1和thr2的取值范圍分別為1.3～1.6和0.7～0.9，這主要是根據(jù)數(shù)據(jù)的相關(guān)性確定的，并通過實(shí)驗(yàn)得到，閾值系數(shù)對波長選擇的影響將在后面討論部分進(jìn)行詳細(xì)闡述。

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

2.1 樣本數(shù)據(jù)集

選取兩個(gè)公共數(shù)據(jù)集進(jìn)行試驗(yàn)：數(shù)據(jù)集1為一組柴油樣本的近紅外光譜數(shù)據(jù)，來自于Eigenvector網(wǎng)站，包含246個(gè)柴油樣本，每條光譜有401個(gè)波長變量，光譜范圍為750～1 550 nm。響應(yīng)變量為柴油50%回收率下的沸點(diǎn)(BP50)值。

數(shù)據(jù)集2為一組土壤樣本近紅外光譜數(shù)據(jù)，來自于Quality &Technology網(wǎng)站，包含108個(gè)土壤樣本，每條光譜有1 050個(gè)波長變量，光譜范圍為400～2 500 nm。響應(yīng)變量為土壤有機(jī)質(zhì)(soil organic matter，SOM)的含量。

2.2 樣本集劃分

近紅外光譜數(shù)據(jù)集被分為一個(gè)校正集和一個(gè)獨(dú)立預(yù)測集。校正集用于建立校正模型，將校正集進(jìn)一步劃分為訓(xùn)練集(calibration set)和驗(yàn)證集(validation set)，利用驗(yàn)證集對訓(xùn)練集的標(biāo)定模型的誤差進(jìn)行評估，即驗(yàn)證集是用來指導(dǎo)選擇候選波長子集變量個(gè)數(shù)的集合。獨(dú)立的預(yù)測集(prediction set)用來評估所生成模型的性能，它不用于校正和驗(yàn)證程序的任何步驟。通過SPXY(sample set portioning based on joint x-y distance)算法將數(shù)據(jù)集劃分為一個(gè)校正集(60%的樣本)、一個(gè)驗(yàn)證集(20%的樣本)和一個(gè)獨(dú)立的預(yù)測集(20%的樣本)。各數(shù)據(jù)集的樣本集劃分?jǐn)?shù)及樣本化學(xué)性質(zhì)統(tǒng)計(jì)結(jié)果如表2所示。

表2 樣本化學(xué)性質(zhì)含量統(tǒng)計(jì)結(jié)果Table 2 Results of sample chemical property

2.3 模型建立與評價(jià)

偏最小二乘回歸(partial least squares regression,PLSR)是光譜分析中常用的多元統(tǒng)計(jì)數(shù)據(jù)分析方法。PLSR適用于預(yù)測變量高度共線性，特別是當(dāng)預(yù)測變量大于樣本數(shù)時(shí)，該方法特別有效。然而PLSR通常存在潛在變量與原始變量相比難以解釋的問題。

多元線性回歸(multiple linear regression,MLR)是一種簡單、常用的校正方法。MLR采用最小二乘法進(jìn)行回歸計(jì)算，其優(yōu)點(diǎn)是方便解釋波長變量對因變量的影響，利用該方法所得到的模型比PLSR模型更易于解釋。但其缺點(diǎn)是存在多重共線性，當(dāng)變量數(shù)大于樣本數(shù)時(shí)無法實(shí)現(xiàn)，因此MLR一般在提取特征波長之后再建模。因?yàn)楸疚奶岢龅腗CC算法消除了變量之間的多重共線性，MCC波長選擇結(jié)果適合通過MLR方法建模。

波長選擇的一個(gè)最主要的目的就是提高模型的預(yù)測能力，而模型的預(yù)測能力主要通過模型的決定系數(shù)(R2)和均方根誤差(root mean squared error,RMSE)指標(biāo)來評價(jià)。其中，決定系數(shù)又稱為擬合優(yōu)度，取值范圍為[0-1]，R2越接近1，自變量對因變量的解釋程度越高。RMSE用來衡量實(shí)測值與模型預(yù)測值之間的偏差，RMSE的值越小，說明模型預(yù)測精確度越高。

3 結(jié)果與討論

3.1 柴油數(shù)據(jù)結(jié)果

柴油數(shù)據(jù)的原始近紅外光譜圖如圖3(a)所示，通過窗口寬度為11的Savitzky-Golay(S-G)一階求導(dǎo)法進(jìn)行預(yù)處理，預(yù)處理后的近紅外光譜如圖3(b)所示，后續(xù)的波長選擇和建模均在圖3(b)數(shù)據(jù)基礎(chǔ)上進(jìn)行。對原始光譜進(jìn)行S-G導(dǎo)數(shù)預(yù)處理能夠提高光譜分辨率，減少變量間的線性相關(guān)性，為后續(xù)的波長選擇奠定基礎(chǔ)。

圖3 柴油樣本原始近紅外光譜圖(a)和預(yù)處理后的近紅外光譜圖(b)Fig.3 Original (a)and Preprocessed (b)NIR spectra of diesel fuels sample

對于柴油數(shù)據(jù)，設(shè)置thr1=1.45，thr2=0.85，通過MCC方法選擇20個(gè)變量，選擇出的變量分布如圖4所示。分別為1 272，1 264，1 266，1 232，1 288，1 294，1 078，1 040，1 274，1 196，1 076，1 230，1 234，1 330，1 316，1 036，1 286，1 338，1 262和1 080 nm。20個(gè)波長中，任一波長與其他19個(gè)侯選波長間的相關(guān)系數(shù)的平均值和標(biāo)準(zhǔn)差如圖5所示，各波長與其余所選波長間的相關(guān)系數(shù)平均值在0.3～0.4之間，處于弱相關(guān)水平。標(biāo)準(zhǔn)差在0.2～0.3之間，相關(guān)系數(shù)相對集中。由此可見，所選變量之間弱相關(guān)，變量間的共線性較小。并且MCC選擇的波長集中在1 100和1 200～1 300 nm附近，柴油是輕質(zhì)石油產(chǎn)品，復(fù)雜烴類混合物，而烴類化合物是由碳與氫原子所構(gòu)成的化合物。MCC方法選取的特征波長都處于C—H鍵的三倍頻和二倍頻吸收區(qū)域。

圖4 MCC波長選擇結(jié)果Fig.4 Wavelength selection results by MCC

圖5 所選波長間的相關(guān)系數(shù)平均值與標(biāo)準(zhǔn)差Fig.5 Mean and standard deviation of correlation coefficient between selected wavelengths

為證明MCC波長選擇算法的有效性，我們對數(shù)據(jù)集通過SPA，CARS，IRIV及RF算法進(jìn)行波長選擇，建立相應(yīng)回歸模型，比較預(yù)測精度。通過文獻(xiàn)[8]可知，各方法的變量初始化、建模方法各不相同，如MCC是通過計(jì)算所有變量間的最小相關(guān)系數(shù)，建立MLR模型，選擇最優(yōu)變量；SPA是計(jì)算全部變量正交子空間上的最大投影值，建立MLR模型，選擇最優(yōu)變量；CARS是通過蒙特卡洛采樣隨機(jī)抽取80%的樣本建立PLS模型；IRIV是通過二進(jìn)制矩陣取樣，建立PLS模型；RF是通過蒙特卡洛采樣，建立PLS模型。根據(jù)各波長選擇算法特點(diǎn)，對SPA和MCC建立MLR模型，對CARS，IRIV和RF建立PLSR模型，會(huì)獲得較優(yōu)的預(yù)測能力。

表3為MCC、全譜及四種常用的波長選擇算法(SPA，CARS，IRIV和RF)所建模型驗(yàn)證集和測試集的模型結(jié)果。由表3可知，基于波長選擇的回歸模型的參數(shù)均不同程度地優(yōu)于基于全譜的PLSR模型(FULL-PLSR)，說明對全譜進(jìn)行特征波長選擇的重要性。對比可以發(fā)現(xiàn)，MCC相比于SPA，CARS和RF三種算法具有明顯的優(yōu)勢，IRIV較MCC的預(yù)測集均方根誤差略優(yōu)，但在選擇的變量個(gè)數(shù)上，MCC算法略有優(yōu)勢。這個(gè)結(jié)果進(jìn)一步表明，MCC變量選擇方法具有一定的有效性。

表3 基于柴油數(shù)據(jù)集不同波長選擇方法下模型的參數(shù)Table 3 Model results based on different wavelength selection methods for diesel fuels datasets

3.2 土壤數(shù)據(jù)

將土壤數(shù)據(jù)的原始近紅外光譜數(shù)據(jù)通過窗口寬度為11的S-G一階求導(dǎo)法預(yù)處理后，分別通過MCC，SPA，CARS，IRIV及RF進(jìn)行特征波長提取并建立回歸模型，驗(yàn)證集和測試集模型結(jié)果如表4所示。由表4可知，MCC的預(yù)測集均方根誤差相比于SPA，CARS和RF三種算法具有較優(yōu)的預(yù)測精度，IRIV與MCC的預(yù)測效果不相上下。與數(shù)據(jù)集1所得結(jié)果一致，可以證明MCC方法的有效性。

表4 基于土壤數(shù)據(jù)不同波長選擇方法下模型參數(shù)Table 4 Model results based on different wavelength selection methods for soil datasets

3.3 MCC方法優(yōu)勢分析

Zhang等[9]在土壤數(shù)據(jù)集的基礎(chǔ)上，利用SiPLS，SiPLS-GA和SiPLS-GA-SPA三種方法進(jìn)行波長選擇后的數(shù)據(jù)建立PLSR模型，其中基于SiPLS-GA-SPA波長選擇結(jié)果建立的PLSR模型具有最佳預(yù)測精度，RMSEP=1.42。用MCC波長選擇法建立的模型的RMSEP為1.032 3。相較而言，基于MCC的波長選擇算法更具有優(yōu)勢。根據(jù)以往的研究結(jié)果可知，基于SPA方法的波長選擇模型優(yōu)于GA[10]，CARS[11]和IPLS[12]等其他波長選擇方法。由此可見，MCC算法在波長選擇方面具有一定的優(yōu)勢。

從MCC算法可以看出，MCC和SPA的基本原理都是選擇最少冗余信息和最小共線性的變量組合，不同的是，SPA是利用向量的投影分析，通過將多個(gè)向量投影到某個(gè)超平面(前一個(gè)投影向量的法平面)上，比較投影向量的模，選擇模最小的向量作為待選向量，以此達(dá)到該向量與前一個(gè)向量的共線性最小。但是，SPA波長選擇只考慮相鄰兩次選擇波長之間的相關(guān)性，而不考慮所選擇的波長和其他所選波長之間的相關(guān)性，即，SPA波長選擇算法并沒有考慮所有特征波長之間的共線性。

MCC算法以某一波長和其他所有波長間相關(guān)系數(shù)最小為準(zhǔn)則選擇特征波長，使得所選特征波長和其他所有波長之間相關(guān)性最小，從而達(dá)到降低共線性的目的。由于減少了所有特征波長之間的共線性，因此MCC波長選擇的結(jié)果適合使用線性回歸方法建模，所以MCC-MLR模型的預(yù)測精度較優(yōu)。從前面兩個(gè)數(shù)據(jù)集建模所得的結(jié)果來看，相較于其他波長選擇算法所建模型的預(yù)測精度，MCC-MLR模型的預(yù)測精度較優(yōu)，是一種有效的模型。MCC算法在應(yīng)用于近紅外光譜數(shù)據(jù)的波長選擇具有一定的優(yōu)勢，本方法的設(shè)計(jì)思想可以為其他高維數(shù)據(jù)的特征選擇和降維提供一定的參考。

3.4 閾值系數(shù)的選取對MCC方法的影響

根據(jù)算法步驟3，為了保證候選波長集S在某一閾值下不為空集，閾值系數(shù)thr1應(yīng)大于1。另外，為了保證所選擇的波長更有效，閾值系數(shù)thr2應(yīng)小于1。為了獲得更好的thr1和thr2，通過選擇thr1和thr2的不同組合，得到MCC-MLR模型的RMSEV。圖6為以土壤數(shù)據(jù)集為例得到的不同的thr1和thr2閾值系數(shù)取值情況下的模型RMSE分布情況。

圖6 不同閾值系數(shù)下MCC-MLR模型的RMSEV值變化Fig.6 Changes in RMSEV value of MCC-MLR model under different threshold coefficients

從圖6可以看出，隨著thr2變大，RMSEV逐漸變小，當(dāng)thr2大于0.8時(shí)，RMSEV趨于平緩，MCC算法的thr2取值為0.85。隨著thr1的降低，RMSEV有下降趨勢，在thr1=1.45附近取得最小值，以此，MCC算法的thr1的取值為1.45。對于柴油數(shù)據(jù)集來說，當(dāng)thr1取值1.45，thr2取值為0.85時(shí)，MCC-MLR方法也能得到最小RMSE模型。因此，MCC設(shè)置thr1和thr2的默認(rèn)值分別為1.45和0.85。對于其他數(shù)據(jù)集，兩個(gè)閾值系數(shù)的取值可能有細(xì)微波動(dòng)，可以通過實(shí)驗(yàn)確定。

4 結(jié) 論

針對變量間的多重共線性問題，提出一種基于最小相關(guān)系數(shù)的近紅外光譜波長選擇方法，該方法以所選變量之間共線性最小化，并且輸入變量對響應(yīng)變量的影響最大化為主要目的，消除模型間的共線性，提高模型預(yù)測精度。為了驗(yàn)證所提出算法的有效性，利用兩組公開的近紅外光譜數(shù)據(jù)集對該方法進(jìn)行了測試。結(jié)果表明，MCC算法獲得了良好的預(yù)測性能，是一種有效的波長選擇算法。MCC方法可以為其他類型高維數(shù)據(jù)降維提供參考。