程介虹, 陳爭光*, 張慶華

(1.黑龍江八一農墾大學電氣與信息學院, 黑龍江 大慶 163319 2.大慶技師學院計算機工程系, 黑龍江 大慶 163524)

土壤作為地球表面五大典型地物(水域、農田、山脈、城市、土壤)之一，其可見光-近紅外反射光譜特性研究對當前精準農業(yè)、土壤數(shù)字制圖、土壤資源遙感調查等工作有著重要作用, 備受眾多科研人員的關注[1]。土壤有機質(soil organic matter,SOM)泛指土壤中來源于生命的物質，是土壤固態(tài)部分的重要組成成分，能促進植物的生長發(fā)育，改善土壤的物理性質。研究表明，近紅外光譜對土壤生物量碳和生物量氮、全碳和全氮以及有機質含量均有較好的評估效果[2]。由于近紅外光譜檢測技術具有快速、成本低、無污染等優(yōu)點，所以被用作一種“快速無損”的技術來檢測SOM含量，可以避免大量的土壤取樣。

近年來，王昶等[2]結合偏最小二乘法對長江中下游糧食生產區(qū)水稻土的土壤質量(全碳、全氮、碳氮比和土壤pH)進行評估，決定系數(shù)均大于0.9，效果較好，可進行快速測定；李偉等[3]用偏最小二乘和人工神經網絡方法對近紅外光譜建模，可測定土壤堿解氮、速效磷和速效鉀含量；吳金卓等[4]結合偏最小二乘和主成分分析法，建立土壤全氮和堿解氮含量測定的定標模型，在野外快速測定土壤全氮和堿解氮含量。徐夕博等[5]通過主成分波段信息判別提取6個主成分，建立MLR和BPN預測模型，發(fā)現(xiàn)BPN模型預測精度優(yōu)于MLR模型。單海斌等[6]通過對土壤光譜進行倒數(shù)、對數(shù)、一階微分、倒數(shù)的一階微分、對數(shù)的一階微分變換，運用單相關分析法提取土壤光譜特征波段，建立多元回歸方程。Ge等[7]通過外部參數(shù)正交化建立偏最小二乘回歸模型，對不同含水率土壤樣品的粘土和有機碳含量進行了預測。Kirshnan等[8]采用逐步多元線性回歸對土壤反射率數(shù)據(jù)進行分析，確定土壤有機質含量預測的最佳波長為564和624 nm。 Brown等[9]以PLSR和BRT回歸證實大量特征波段都可用來對SOM進行估算。Wang等[10]利用不同預處理方式改善光譜與土壤性質間相關性。Bilgili等[11]使用PLSR和MARS對土壤有機質含量進行建模研究，取得良好的預測精度。

上述研究多基于全譜或特征波段建模來估測土壤的養(yǎng)分信息，但近紅外光譜所采集的樣本變量數(shù)目龐大，影響建模速度，并且近紅外波段是合頻和倍頻的吸收區(qū)域，光譜信息重疊嚴重。若選用波段信息，變量中會出現(xiàn)干擾、共線性變量，所以僅僅建立一部分波段與土壤某個屬性之間的線性關系無法達到最優(yōu)估算結果[5]。

為了解決這一問題，波長選擇成為近紅外光譜分析中非常重要的一步。特征波長選擇是從全譜數(shù)據(jù)中提取部分涵蓋有用信息的光譜，去除噪聲光譜及無用信息，建立一個更為簡約、穩(wěn)定的近紅外光譜模型。波長選擇可大大減少變量數(shù)目，加快模型的計算效率，提高模型的穩(wěn)健性[12]。本文以土壤為研究對象，Rinnan等[13]已證明近紅外光譜與SOM含量密切相關，相關系數(shù)超過0.9。為提高其預測精度，本文選取連續(xù)投影算法(successive projections algorithm,SPA)消除共線性信息變量，偏最小二乘法(interval partial least squares,IPLS)選擇特征波長區(qū)間間隔，競爭自適應重加權采樣法濾除冗余信息變量(competitive adaptive reweighted sampling, CARS),從三個不同方面有代表性算法進行特征波長選擇，建立模型，比較其預測性能。

1 材料與方法

1.1 樣本數(shù)據(jù)

土壤樣本數(shù)據(jù)來自于Quality & Technology[13]，是Havstrom等在瑞典北部阿比斯科(68°21′N，18°49′E)進行長期野外試驗時，模擬了亞北極氣候的預測變化對SOM的影響而得到的樣本數(shù)據(jù)。該數(shù)據(jù)為一組公開的數(shù)據(jù)集，基于此數(shù)據(jù)的方法改進和所建模型具有一定的可比性。樣本數(shù)據(jù)包含兩部分，樣本的NIR光譜及其化學性質。共計108個樣本。

樣本光譜的波長范圍為400～2 500 nm，采樣間隔為2 nm。樣本的化學特性為土壤有機質含量。本文采用SOM含量作為因變量進行近紅外光譜數(shù)據(jù)建模預測分析。

1.2 變量選擇算法

1.2.1SPA算法步驟 SPA是由Araujo等[14]提出的一種消除變量間共線性的波長選擇算法，該方法通過提取全譜的幾個特征波長消除原始光譜中冗余的信息[15]。其主要原理是利用向量的投影分析，假設已給出初始波長k(0)和所需提取波長數(shù)目N，算法步驟[14,16]如下。

Step 0:在第一次迭代(n=1)之前，將校正集Xcal的第j列光譜數(shù)據(jù)賦值給xj,j=1,…,J，J為總波長數(shù)。

Step 1:沒有被選擇的列向量的集合記為S,S={j,1≤j≤J,j?{k(0),…,k(n-1)}}

Step 3:令k(n)=arg(max‖Pxj‖),j∈S

Step 4:令xj=Pxj,j∈S

Step 5:令n=n+1，如果n

End:最后得到的波長為{k(n);n=0,…,N-1}

其中，k(0)和N的選擇是很關鍵的一個步驟。k(0)取值在1～J之間變化，而N的變化范圍為1≤N≤Mcal，Mcal為校正集樣本數(shù)，這是可以由SPA選擇的最大波長數(shù)。

1.2.2CARS算法步驟 CARS是由Li等[17]在2009年提出來的一種結合蒙特卡洛(Monte Carlo,MC)采樣與偏最小二乘(partial least squares,PLS)模型回歸系數(shù)的一種波長選擇方法，通過消除變量來建立一個高性能的校準模型。該方法可以對無信息變量進行有效去除，最終選擇出對預測目標起關鍵作用的變量。

具體步驟為：首先，隨機抽取80%的樣本做為校正集建立PLS模型，計算回歸系數(shù)的絕對值和每個波長點對應的權重，然后，利用指數(shù)衰減函數(shù)(exponentially decreasing function，EDF)和自適應重加權采樣法(adaptive reweighted sampling，ARS)對變量進行選擇，通過交叉驗證的方法計算交叉驗證均方根誤差(RMSECV)。N次MC采樣后選擇N個子集，得到N個RMSECV，選擇RMSECV最小的波長子集，該子集所包含的變量即為最優(yōu)變量組合[18]。

1.2.3IPLS算法步驟 IPLS是由Norgaard等[19]于2000年提出的一種特征波長區(qū)間選擇方法，它的主要作用是在不同的光譜分支中提供有關信息的整體圖像，從而聚焦于重要的光譜區(qū)域，并消除來自其他區(qū)域的干擾[19]。具體步驟為：將全譜分成若干相同寬度的光譜子區(qū)間，然后在每個子區(qū)間建立偏最小二乘回歸模型，通過比較區(qū)間模型和全譜模型的RMSECV來比較預測性能，RMSECV最低時的區(qū)間模型即為最優(yōu)波段。

1.3 模型建立與評價

偏最小二乘回歸(partial least squares regression,PLSR)建模方法在近紅外光譜分析中使用較為廣泛，它考慮光譜數(shù)據(jù)與性質之間的內在聯(lián)系，模型更加穩(wěn)健，但模型建立過程復雜、抽象，計算速度較慢且繁瑣。當各變量高度線性相關時，用PLSR法非常有效[20]。而多元線性回歸(multiple linear regression,MLR)的特點適用于線性較好的數(shù)據(jù)，不考慮參數(shù)之間的相互干擾、計算簡單、公式含義清楚，產生的模型比PLSR更簡單，更容易解釋。故本文基于所選取的特征波長，通過建立MLR模型分析其預測能力。

MLR是一種常規(guī)的校正方法，直觀簡單，且具有良好的統(tǒng)計特性，應用非常普遍。MLR采用最小二乘法進行回歸計算，模型一般為

近紅外光譜模型的預測能力主要通過變量數(shù)(n)、校正相關系數(shù)(Rc)、校正均方根誤差(RMSEC)、預測相關系數(shù)(Rp)、預測均方根誤差(RMSEP)指標來評價。其中，R取值越接近1，RMSEC和RMSEP越接近0，模型的擬合性越好，預測精度越高[4]。

1.4 數(shù)據(jù)分析工具

所用軟件包括MATLAB R2015b和The Unscrambler X 10.3 (64-bit)。The Unscrambler軟件具有廣泛的數(shù)據(jù)預處理方案,確保數(shù)據(jù)都適合進行多因素分析，并且具有先進的回歸、分類和預測建模工具。本文光譜數(shù)據(jù)的預處理、建模分析在Unscrambler軟件中實現(xiàn)，各變量選擇方法、圖形的繪制在MATLAB中實現(xiàn)。

2 結果與分析

2.1 光譜數(shù)據(jù)特征

各樣本400～2 500 nm區(qū)間原始近紅外光譜變化趨勢相同，無明顯異常樣本。為了消除基線漂移、消除背景影響、提高光譜分辨率，使用窗口為5的二階Savitzky-Golay(S-G)求導法對光譜進行預處理，結果如圖1所示。

圖1 原始光譜及預處理后光譜Fig.1 Original spectrum and pre-processed spectrum

將108個樣本通過SPXY(sample set portioning based on joint x-y distance)算法分為75%訓練集和25%預測集，建模集包含81個樣本，預測集包含27個樣本， 土壤有機質含量統(tǒng)計數(shù)據(jù)結果如表1所示。劃分后建模集的SOM含量范圍涵蓋預測集的SOM含量，證明建模集具有代表性。

圖1 原始光譜及預處理后光譜Fig.1 Original spectrum and pre-processed spectrum

分別對原始光譜和預處理后的光譜建立PLSR模型，結果如表2所示。經S-G求導預處理后模型的預測性能有所提高，故將預處理后的數(shù)據(jù)作為后續(xù)波長選擇的數(shù)據(jù)。

表2 原始光譜和預處理后的光譜PLSR模型參數(shù)對比Table 2 Comparison of PLSR models of the raw spectrum and pretreated spectrum

2.2 特征波長選取

2.2.1SPA變量選擇結果 由SPA方法選擇特征波長建立MLR模型，選取RMSEP最小值對應的波長個數(shù)為最終的特征波長個數(shù)。圖2A中的正方形標記所示為SPA多元線性回歸模型選擇的變量數(shù)，RMSEP的最小值為1.214 4，之后RMSEP值基本達到穩(wěn)定，此時選擇特征波長數(shù)為6個。圖2B中方塊所對應的6個點即為SPA選擇的最佳特征波長,分別為636、688、1 934、2 224、2 330、2 426 nm。

2.2.2CARS變量選擇結果 利用CARS 算法對光譜數(shù)據(jù)進行特征波長選擇，MC抽樣的次數(shù)設置為100次。圖3A所示為變量變化的趨勢，由于指數(shù)衰減函數(shù)的作用，在剛開始進行MC采樣時進行快速選擇，變量減少的速度非?？?。然后在精細選擇中，使用自適應重加權采樣法來選擇基于回歸系數(shù)的關鍵變量，變量減少的速度較慢，直到得到最優(yōu)子集。圖3B所示為十折交叉驗證RMSECV的變化趨勢，隨著MC采樣次數(shù)的遞增，PLS模型的十折交叉驗證RMSECV的值逐漸減少，可知剔除了大量無信息變量。但當達到最小值后(虛線處)又逐漸增大，表明光譜數(shù)據(jù)中一些與預測相關的信息被剔除，導致模型的性能變差。圖3C所示為回歸系數(shù)路徑，星號垂線處是為第55次抽樣，此時的RMSECV值最小為2.032 6，有34個回歸系數(shù)波長為非零值，即特征波長共計34個變量。

A:模型包含的變量數(shù)和模型RMSEP的關系;B：特征波長點。A: Relationship between the number of variables and RMSEP of the model; B: Characteristic wavelength points.圖2 SPA的特征波長選擇Fig.2 Characteristic wavelength selection of SPA

A:變量變化趨勢；B:十折交叉驗證RMSECV的變化趨勢；C:回歸系數(shù)路徑。A: Variation trend of variables; B: Trend of of RMSECV under 10-fold cross validation; C:Regression coefficients path.圖3 CARS的特征波長選擇Fig.3 Characteristic wavelengths selection of CARS

2.2.3IPLS變量選擇結果 利用IPLS算法對光譜數(shù)據(jù)進行特征波長選擇，將預處理后的光譜數(shù)據(jù)依次等分成20、25、30、35、40個區(qū)間，然后在每個子區(qū)間建立PLS模型，比較這些區(qū)間模型和全譜模型的RMSECV，選擇效果最佳的子區(qū)間。IPLS選擇結果如表3所示。

從表3可以看出，當光譜區(qū)間劃分為30時，對應的RMSECV最小。因此,將光譜數(shù)據(jù)等分成30個區(qū)間進行特征波長選擇。由圖4A可知，當全譜的主成分數(shù)為3時，RMSECV基本平穩(wěn)。由圖4B可知,建立在第15區(qū)間的PLS模型的RMSECV最低。因此選擇1 392～1 460 nm區(qū)間，共計35個波長點作為波長選擇結果。

表3 IPLS子區(qū)間優(yōu)選結果Table 3 Results of IPLS subinterval optimization

A:RMSECV與全譜模型的主成分數(shù)；B:各子區(qū)間模型和全譜模型RMSECV對比，各柱狀圖上斜體數(shù)字代表子區(qū)間的主成分數(shù)。A:RMSECV versus PLS components for global model；B:Comparison of RMSECV between each interval model and full spectrum model, and the italic number on each rectangle represents the number of principal components of the subinterval.圖4 IPLS的特征波長選擇 Fig.4 Characteristic wavelengths selection of IPLS

2.3 模型建立與比較

本研究擬利用MLR對全譜以及經三種特征波長選擇算法后得到的數(shù)據(jù)進行建模，比較預測精度，但由于全譜含有1 050個波長點，數(shù)據(jù)過于龐大，而MLR只適用于對波長點少于樣本數(shù)的數(shù)據(jù)進行建模，所以對全譜進行PLS回歸；而使用IPLS方法進行特征波長選擇，所選擇的波長點是連續(xù)的，存在共線性，所以對IPLS選擇的35個特征波長點進行PLS回歸；對SPA選擇的6個特征波長點和CARS選擇的34個特征波長點依然建立MLR模型，得到的模型的校正、預測相關系數(shù)和校正、預測均方根誤差的值如表4所示。

表4 不同波長選擇方法所得模型對比Table 4 Comparison of models obtained by different wavelength selection methods

由表4可以發(fā)現(xiàn)，經過特征波長提取后的預測精度均高于原始光譜，進一步證明了對全譜進行特征波長選擇的重要性。通過對比三種特征波長選擇算法的預測效果，發(fā)現(xiàn)SPA僅挑選全譜的6個波長點，RMSEP為1.214 4，預測效果最好。SPA不僅降低了建模的復雜度，還大大提高了模型的預測精度，這種方法可以有效地應用于特征波長的提取中，簡化模型復雜度，提高模型計算效率。

圖5分別為對預測集樣本不同模型的土壤有機質實測值和預測值相關圖。其中，SPA-MLR模型得到的方程如下。

圖5 不同模型下土壤有機質的實測值和預測值相關性Fig.5 Correlation between measured and predicted values of soil organic matter under different models

y=-12.00-5 410.73x(636 nm)-375.76x(688 nm)-15 957.55x(1 934 nm)+23 717.22x(2 224 nm)-50 740.52x(2 330 nm)+45 713.68x(2 426 nm)

3 討論

王巖[21]利用此數(shù)據(jù)，通過IPLSR、PLSR對土壤有機質含量建立模型，預測模型的RMSEP分別為2.787 3和2.709 8，而本文使用三種變量選擇方法建立模型的預測能力均得到不同程度的改善，進一步證明了對全譜進行特征波長選擇的重要性。文中SPA預測能力最優(yōu)，CARS次之。基于SPA選擇的特征波長MLR模型預測集的相關系數(shù)為0.970 2，RMSEP為1.214 4。并且僅選取全譜的6個波長點，大大降低模型復雜度，提高模型效率。這與吳龍國等[22]應用不同波長選擇方法對土壤含水率進行建模研究時得到的結論一致，相比于CARS，SPA-MLR模型的預測精度更優(yōu)。Chen等[23]利用IPLS和SPA結合MLR建模，對谷物的蛋白質、總碳水化合物和粗脂肪的含量進行預測，發(fā)現(xiàn)SPA-MLR的預測精度優(yōu)于IPLS，這也與本文結果一致，說明SPA-MLR在某些光譜模型中具有一定的有效性。

在SPA方法中，通過不斷提高MLR模型的預測能力來選擇變量。與PLS方法相比，SPA投影只用于選擇變量，不會修改原始數(shù)據(jù)，因此保留了光譜數(shù)據(jù)和樣本化學性質間的關系。SPA在進行波長選擇前，對各變量進行標準化處理，各波長點向量模相同，因此進行投影時選擇最大投影向量的過程其實就是選擇相互正交的向量的過程，所以SPA可以將變量之間的共線性降到最低。而MLR對線性較好的數(shù)據(jù)回歸效果好，所以可以解釋SPA-MLR的預測性能最好。

土壤的NIR與SOM之間存在較好的相關性，其預測集的相關系數(shù)大于0.94，更加證明了土壤的光譜特性可以對SOM含量進行預測。本研究通過三種不同特征波長提取算法對土壤有機質含量進行建模，取得了較滿意的成果?，F(xiàn)有研究大多對土壤有機質、重金屬及氮磷鉀的含量進行測定，土壤的微生物含量對農業(yè)生產也有著影響，今后應對土壤的微生物含量的預測進行進一步研究。