基于CARS算法的不同類型土壤有機質高光譜預測

唐海濤，孟祥添，蘇循新，馬 濤，劉煥軍,4，鮑依臨，張美薇，張新樂※，霍海志

（1. 東北農業大學公共管理與法學院，哈爾濱 150030；2. 黑龍江省地質資料檔案館，哈爾濱 150030；3. 黑龍江省第五地質勘察院，哈爾濱 150030；4. 中國科學院東北地理與農業生態研究所，長春 130012）

0 引言

土壤有機質（Soil Organic Matter，SOM）可以通過生物合成和分解，改善土壤的物理、化學和生物特性[1]，在控制土壤功能和質量、抵消溫室氣體排放、完善全球碳循環系統信息等方面發揮著重要作用[2]。高光譜預測模型為實現SOM等土壤屬性速測與遙感反演以及表層碳庫估算等提供數據信息[3]，并為SOM速測儀器研制、土壤制圖與退化監測、精準農業實施等提供數據與技術支持[4]。高光譜技術具有精細的光譜分辨率，可獲取地物納米級的連續光譜信息，SOM具有多種官能團（如羥基、羧基等），分別在紅外光譜區域有特征性吸收，且不同波段的吸收強度與該物質的分子結構及濃度存在對應關系，因此，紅外光譜可以反映SOM含量，為其定量估算提供了一種有效的手段，為預測SOM提供了可能[5]。黑龍江省海倫市位于世界三大黑土地分布區之一的松嫩平原東北端，土壤類型多樣，其中黑土面積達到全市面積1/2以上，且是中國重要的商品糧基地，了解其SOM的分布情況、空間變化規律，有利于科學評價土壤的質量情況并對農場合理施肥提供指導，對耕地資源的可持續利用具有十分重要的實際意義，可為海倫市耕地的可持續利用和土壤質量保護監測提供技術支持，為將來海倫市土地管理建立完整的空間土壤信息系統提供框架。

以往室內高光譜對于SOM的輸入變量研究多停留在以全波段反射率或對應的數學變換上，選取相關系數較大的波段進行建模，該方法僅考慮了SOM與光譜間的關系，并沒有考慮光譜間的重疊吸收或相互影響[6]。利用光譜指數技術預測SOM的研究成為當前熱點，光譜指數是由幾個窄波段或寬波段組合而成，可通過分析特定波段間的相互作用，提高對待測屬性的敏感程度[7]，有助于挖掘波段間的隱晦信號[8]。SOM空間分布特征受到高程、坡度、坡向等地形因子不同程度的影響，地形條件影響其物質循環過程和強度[9]，通過數字高程模型（Digital Elevation Model，DEM）提取高程作為模型輔助變量參與建模。同時特征波段選擇是進行SOM含量預測的一個重要方面，已經引起了越來越多學者的關注。土壤光譜反射數據通過競爭自適應重加權采樣（Competitive Adaptive Reweighted Sampling，CARS）篩選出的特征波段不僅將輸入波段壓縮至全波段數目的一半以下，同時提升了模型估測精度，降低了變量維度和模型復雜度[10]，Vohland等[11]發現，在60個農業樣品的土壤屬性預測中，CARS算法減少了建模時間，且能夠合理、精確、有效的確定特征波段在全波段中的位置。以往的學者多以一種類型的土壤為對象，進行SOM高光譜響應特性研究，但是由于土壤的光譜反射率是土壤內在理化特性光譜行為的綜合反應，不同類型土壤的光譜特征不同[12]，因此模型普適性較弱。盧艷麗等[13]利用不同土壤類型分組試驗分析了東北平原土壤光譜反射率曲線形狀變化，確定了8種不同類型土壤與原始光譜反射率的相關敏感波段并建立了同質性SOM預測線性模型，從而達到簡化SOM預測模型的目的。Bao等[14]對比了多種土壤分組策略下SOM的預測精度，同時引入競爭自適應重加權采樣方法進行模型輸入量的篩選，證實了土壤分類的優勢與多輸入量降維的必要性。因此，不同類型土壤分別提取輸入變量進行高光譜SOM預測將有利于分析各類土壤的理化性質，從而提高預測精度。

已有SOM高光譜預測研究常基于一種土壤類型建立模型或者多種土壤類型進行全局回歸建模，且輸入變量的類型較為單一，存在SOM預測精度不高的情況[15]。為了充分考慮土壤光譜信息及影響因素，本研究以海倫市為研究區域，根據全國第二次土壤普查結果及對采樣點的地理位置對土樣進行分類。在土壤分類的前提下，以土壤光譜反射率數據、DEM數據以及光譜指數作為輸入變量，建立基于隨機森林算法（Random Forest，RF）的分類高光譜SOM預測模型。為了降低輸入量之間的共線性，引入CARS算法篩選特征波段，提高不同類型SOM預測的精度，以期實現動態快速預測SOM含量。

1 材料與方法

1.1 研究區概況

海倫市位于松嫩平原的中心地帶，地理位置在46°58"N～47°52"N，126°14"E～127°45"E之間，屬溫帶大陸性季風氣候，地勢平坦，土質肥沃，耕地面積廣闊，是國家重要的商品糧基地。其土壤類型主要為黑土、草甸土和沼澤土，在該研究區內還有少量的水稻土、暗棕壤及白漿土。黑土土層深厚，結構良好，富含SOM和腐殖質，自然肥力高。沼澤土所處的地勢大都比較低洼，SOM累積明顯。由于該區地形高程差較大，加上耕地的長期粗放利用導致土壤侵蝕嚴重，降水將地勢較高的土壤沖積到地勢較低的草甸土表面，導致表層草甸土性質較為復雜多樣[16]。海倫市主要土壤類型（全國第二次土壤普查結果）及采樣點分布圖和海倫市30 m空間分辨率的DEM數據見圖1。

1.2 土樣樣品

2019年5月15—20日，于作物出苗前，沿主要鄉級以上道路，在黑龍江省海倫市全市進行樣本采集。選擇土壤裸露的地區作為樣區，考慮土地利用類型和土壤類型采集0～20 cm耕層土壤。為保證采樣點的有機質含量能夠代表采樣點附近一定空間內的SOM水平，采用四分法收集樣品，同時利用GPS記錄采樣點經緯度，總共采集土壤樣本548個。采集的樣品經過風干，研磨，過2 mm篩。每個樣品分2份，一份用于光譜測量；一份用于SOM含量分析。SOM含量用高溫外熱重鉻酸鉀氧化容量法測定[17]。

1.3 光譜測量及數據預處理

采用ASD FieldSpec○R3便攜式光譜儀在暗室內對風干土進行光譜測試。光譜測試流程詳見文獻[18]。由于反射率波譜在400～430和2400～2500 nm范圍內噪聲較為強烈，為減少高頻噪聲的干擾，本文選取光譜反射率波譜范圍為430～2400 nm，并對其進行9點平滑、10 nm重采樣處理，此過程分別在EXCEL和ENVI 5.3中實現。

不考慮土壤空間差異性，將整個土壤樣本作為全局回歸預測數據集。同時，土壤樣本根據全國第二次土壤普查圖，利用ArcGIS 10.1中的工具箱提取每個土壤樣本的土壤類型，將土壤樣本劃分為不同土壤類型，同一種土壤具有相同光譜表現特征的土壤樣本集。根據中國土壤分類，土壤類型可分為黑土、草甸土、沼澤土，然后針對不同分類樣本進行局部回歸預測建模。

1.4 輸入量提取

國內外學者進行SOM高光譜估測時，輸入量多選擇為高光譜反射率或光譜吸收特征建立模型，輸入變量類型結構單一，容易忽略土壤高光譜反射率之間的高度共線性[19]。本研究通過CARS算法挑選的特征變量、光譜指數結合DEM數據作為模型輸入變量。

1.4.1 光譜指數

在高光譜數據預測SOM的研究中，為了確定敏感的波段，必須從SOM含量信息中獲取深度信號，因此光譜指數常作為一個重要指標[20]。本文探討歸一化指數（Normalized Difference Index，NDI）、再歸一化指數（Renormalized Difference Vegetation Index，RDVI）、比值指數（Ratio Index，RI）與SOM含量之間的關系。

表1 光譜指數及公式 Table 1 Spectral indices and formula

1.4.2 地形因素

地表微氣候、土壤中的水分運動以及物質的重新分配進程，都受到地形的影響[25]。在美國地質勘探局網站（http://www.usgs.gov/）下載DEM數據，其空間分辨率為30 m。在ArcGIS 10.1中，利用Spatial Analyst Tools中的Extract Multi Values to Points工具，提取出每個采樣點的DEM，將DEM作為模型的輸入變量。

1.4.3 競爭性自適應加權算法

土壤高光譜數據量大、存在光譜信息冗余和重疊現象，通過CARS算法挑選特征變量可以降低光譜波段之間的高度共線性問題，從而提高預測模型的精度及速度。CARS算法將各波段變量作為單一個體，在進行個體選擇的過程中，保留具有較強適應能力的個體。其具體步驟為：首先，隨機抽取固定比率的樣本作為校正集建立PLS模型，計算回歸系數的絕對值和每個波段點對應的權重，然后利用指數衰減函數（Exponentially Decreasing Function，EDP）和自適應重加權采樣法（Adaptive Reweighted Sampling，ARS）對變量進行選擇，通過交叉驗證的方法計算交叉驗證均方根誤差（Root Mean Square Error of Cross-Validation，RMSECV），N次蒙特卡羅采樣后選擇N個子集，得到N個RMSECV，選擇RMSECV最小的波段子集，該子集所包含的變量即為最優變量組合[14,26]。本次試驗在MATLAB 2014a軟件環境中運行CARS算法。由蒙特卡羅交叉驗證法選擇最優潛在波段變量，其中將蒙特卡羅采樣次數設定為100，對采樣次數進行反復迭代，通過對比各次采樣的RMSECV值，當其值最小時，相應采樣次數的變量被篩選為最優變量子集。

1.5 模型構建與驗證

RF是基于決策樹分類集成算法，其中每一棵樹都依賴于一個隨機向量，通過對數據集的列變量和行變量觀測進行隨機化，生成多個分類樹，最終將分類樹結果進行匯總。RF對于非線性問題有很好的解釋能力，降低了運算量的同時也提高了預測精度[27]。本試驗在R語言中，利用‘Random Forest’工具包進行預測，在進行擬合前，分別對需要生成樹的數量（ntree）參數設定為500，每個節點用于分割節點的預測變量樹（mtry）參數設定為1/3總變量數[28]。

模型構建按照建模集與驗證集2∶1的比例選取樣本。以CARS篩選后土壤高光譜反射率數據、DEM以及光譜指數為自變量，SOM含量作為因變量，運用RF，構建SOM預測模型。使用調整后決定系數（R2adj）、均方根誤差（RMSE）以及性能與四分位間隔距離的比率（Ratio of Performance to Interquartile distance，RPIQ）為精度評價指標。R2adj越大、表明模型越穩定；RMSE越小、表明模型精度越高；RPIQ同時考慮了預測誤差和觀測值的變化，提供了一個更客觀、更容易在模型驗證研究中進行比較的模型有效性度量。RPIQ越大，模型的預測能力越強。與殘差預測偏差不同，RPIQ對觀測值的分布沒有任何假設[29]，其公式如下：

式中IQ是第三和第一個四分位數之間的差值。

2 結果與分析

2.1 SOM描述統計

土壤樣本SOM含量統計特征見表2，質量分數最大值為11.38%，最小值為0.98%，土壤樣品SOM差異較大，這為全面解析SOM反射光譜特性研究提供了較完整的樣本數據。根據土壤樣本SOM描述統計表的偏度和峰度值可以判斷SOM含量數據呈現非正態分布。在SOM相關的研究中可知SOM質量分數達到2%以上，對土壤光譜特征起主導作用[30]，SOM質量分數小于2%的土壤，其光譜曲線特征易受其他母質等成分的影響，而本次研究中SOM平均含量（質量分數）4.5%以上，能夠充分說明SOM的含量決定了土壤光譜的特征。

表2 土壤樣本有機質含量統計結果 Table 2 Statistical results of organic matter content in soil samples

2.2 CARS算法篩選特征波段

3種土壤類型以及未分類整體在指數衰減函數的作用下，優選變量的數量均隨迭代次數的增加呈指數減少，其RMSECV值整體均呈現先減后升的趨勢。以黑土為例（圖 2），從圖2a可以看出，隨著運行次數增加，被優選出的波段變量數逐漸減少，前5次采樣過程有明顯遞減，此后逐漸平穩。圖2b 整體上在1～47次采樣中，RMSECV值不斷降低，表明篩選過程中剔除的變量與SOM去除量無關，而47次采樣迭代以后，RMSECV值呈回升趨勢，表明反射率光譜中與SOM無關的大量信息或噪聲被添加，從而導致RMSECV值上升。圖2c為所有變量在每次采樣過程中的回歸系數路徑變化圖，圖中各線表示隨著運行次數的增加各波段變量回歸系數的變化趨勢。結合圖2b分析發現當采樣次數為第47次時，RMSECV值最小即所選擇的光譜變量子集最優。草甸土、沼澤土以及未分類整體的RMSECV最小值、相應運行次數及特征波段見表3。

表3 CARS下基于不同土壤類型的特征波段，運行次數和最小交叉驗證均方根誤差 Table 3 Characteristic wavebands, number of sampling runs and minimal RMSECV of different soil types under CARS

從表3可知，通過CARS算法，黑土、草甸土、沼澤土以及整體未分類分別篩選出23、30、21和9個特征波段，輸入波段壓縮至全波段數目的16%以下。黑土特征波段的分布主要在1280～2230 nm近紅外光譜區域，這是由于受到NH，CH和CO等基團的分子振動的倍頻與合頻吸收影響[31]，草甸土在可見光-近紅外光譜區域均有波段選中，其中1700～1790 nm處SOM響應可能是由氧化鋁影響的光譜變化引起的。沼澤土篩選的特征波段在1300～2000 nm比較均勻分布，這主要是由于沼澤土中的大量三氧化物被還原。值得注意的是，波段1450、1470、2150 nm在3種土壤類型中均被選擇，這是由于SOM在1400 nm附近受到土壤黏土礦物質中所含羥基的影響，2220 nm附近存在一個與SOM相關的烷烴特征峰和存在的氫氧化鋁黏土礦物吸收帶影響[32]。沼澤土、草甸土篩選的430、440、530、550、670 nm少量特征波段位于可見光波段，這是由于受到了土壤發色團和SOM本身黑色的影響，可見光波段存在較寬的吸收波段。

2.3 光譜指數的選取依據

光譜指數通過迭代運算，充分考慮波段之間的協同作用，同時最小化無關波段的影響[33]。研究選取的光譜指數是通過文獻查閱，選擇可用來估測SOM的一系列物理和化學參數的相關光譜指數，并結合本次實際采樣點數據進行相關性計算得出。3種土壤類型原始反射率數據與SOM之間的NDI、RDVI、RI指數的相關性均較高，且均通過了P=0.01水平上的極顯著性檢驗（表4）。黑土RI指數與SOM的相關性最高，相關系數為0.757，草甸土RDVI指數與SOM的相關性最高，相關系數為-0.784，沼澤土RDVI指數與SOM的相關性最高，相關系數為0.922。圖3是不同土壤類型的3種光譜指數與SOM含量的二維相關系數矩陣圖。3種土壤類型的SOM敏感波段區域主要集中于短波紅外部分，主要集中在1000、1900和2200 nm附近。

表4 土壤有機質含量與最佳光譜指數的關系 Table 4 Relationship between soil organic matter content and optimal spectral index

2.4 SOM光譜預測模型

由表5可知，黑土、草甸土、沼澤土的驗證集調整后決定系數依次為0.678、0.674、0.768，其中沼澤土精度最高，草甸土精度最低，這是由于沼澤土在積水條件下，空氣隔絕，微生物活動受到強烈抑制，植物殘體不能充分分解，而以粗SOM和半腐爛SOM的形式積累于地表。全局回歸模型R2adj達到0.742，局部回歸模型R2adj達到0.777。通過局部回歸，在一定程度上提高了SOM的預測精度。無論是單一土壤類型，還是整體SOM預測，其R2adj均達到0.67以上，RPIQ均大于1.8，表明該模型能較好實現SOM預測。

表5 不同土壤類型隨機森林預測模型精度 Table 5 Prediction model accuracy of random forest for different soil types

在高光譜SOM預測研究中，波段篩選是一個關鍵方面。本研究通過CARS算法篩選波段與已往學者利用相關分析取相關系數大于0.65篩選出的波段[34]進行建模比較，研究發現CARS算法不僅極大地降低土壤高光譜變量維度和計算復雜程度，驗證集R2adj提高了0.167，精度有一定程度的提升。

表6 不同波段篩選方式隨機森林預測模型精度 Table 6 Accuracy of random forest prediction model with different band screening methods

3 討論

本研究將不同土壤類型（黑土、草甸土、沼澤土）分別進行SOM的預測，取得了較高精度。通過土壤分類進行SOM預測，消除了不同土壤類型由于“向鄰性”導致的反射光譜曲線相似的影響，從而有利于提高預測精度。由于不同類型土壤中礦物成分與SOM含量的差異，造成反射光譜間存在顯著的區別，通過土壤分類，將有利于提取不同類型土壤光譜參數進行SOM預測。陸龍妹等[35]通過全局回歸與局部回歸進行SOM預測比較，依照傳統土壤類型建立各自的有機質光譜預測模型精度并不好，這是由于砂姜黑土和黃褐土2種土壤類型的黏土礦物都存在蒙脫石且含量較高，SOM含量接近，所以2種土壤類型之間光譜曲線特征相似，造成SOM全局回歸精度低。而黑土、草甸土、沼澤土之間黏土礦物存在著較大的差異，因此通過全局回歸與局部回歸比較，全局回歸能夠提高有效信息的獲取程度提高模型精度。其沼澤土的預測精度高于草甸土，這是由于沼澤土土壤濕、土層緊且富有彈性，有機質含量豐富、土體酸堿度從微酸到堿性、土壤顏色深，而草甸土土壤表層砂礫化、有浮沙覆蓋、有機質含量較低、土體呈堿性、質地較粗、細粒物質少、土壤色澤淺有一定的關系。

以往許多學者們采用相關分析法研究SOM與土壤光譜反射率（或其不同數學變換形式）的關系，將相關系數高的波段作為SOM敏感波段。而后，越來越多的學者采用CARS變量優選方法，從全波段中濾除無效變量或冗余變量，優選出敏感波段。本研究基于CARS算法，黑土、草甸土、沼澤土分別選擇23、30、21個特征變量，占全波段數目的11.6%、15.2%、10.6%，極大地縮減了波段信息，解決了SOM預測研究中波段數目多，計算任務繁重的問題。結果表明，CARS篩選的最優子集存在一定的規律性，波段主要集中在1100～2400 nm之間，這主要由于受到羰基、酰胺和羥基等基團的分子振動的倍頻與合頻吸收影響。其中，黑土篩選的特征波段少位于1000 nm以下，這是由于CARS是通過利用線性模型偏最小二乘法作為適應度函數，及交叉驗證不斷優化計算，最終選擇出最優子集而不是常用的相關性分析確定特征波段。已有的相關研究表明：SOM在整個NIR-SWIR范圍比較敏感，李穩冠等[26]將栗鈣土、黑鈣土、灰鈣土、山地草甸土等土壤光譜曲線通過CARS挑選的特征波段，變量主要分布在1900～2400 nm的近紅外光譜區域，在可見-近紅外光譜區域均有分布。CARS對原始光譜進行特征變量篩選，在保證模型精度的同時顯著減少構建模型的變量數。Bao等[14]對黑土、黑鈣土、風沙土、草甸土4種土壤類型通過CARS算法篩選最優變量子集，其波段大多位于1350～2400 nm范圍內，少量位于400～1200 nm。因此，通過CARS算法篩選的特征波段，與已有研究SOM的反射光譜響應波段相吻合。不同土壤類型通過CARS篩選的最優子集也存在差異，其選擇的特征變量具有不穩定性。

通過耦合敏感波段的反射率數值進行數學變換所計算得到的光譜指數，避免了由于原始反射率作為輸入量所造成的數據冗余，以及產生的共線性問題。黑土篩選出的波段主要為1030、1910、1940、1950 nm，草甸土在1420、1340、2150、2230 nm，沼澤土集中在1920和1930 nm。3種土壤類型的篩選的波段都位于NIR-SWIR范圍，這是由于羰基基團的基頻振動和其在NIR-SWIR范圍所對應的酰胺、羥基等基團倍頻和合頻吸收影響，也與以往的研究一致[36]。因此通過將不同類型土壤分別，以CARS篩選的特征波段、DEM數據和光譜指數作為數據源，建立的RF模型能夠有效實現SOM預測，使精度有著顯著的提升。然而，本次研究仍存在不足之處：土壤的光譜反射率還會受到土壤的成土母質、礦物成分、土壤表面粗糙度、粒徑、水分等因素的影響，因此，后續研究在原土室外光譜的基礎上，將考慮更多的影響因素，加強原土室外光譜SOM的估測模型研究，以提升SOM的預測精度。

4 結論

為了解決不同類型土壤預測有機質（Soil Organic Matter，SOM）輸入量類型單一造成精度偏低的問題，本文以海倫市3種土壤類型（黑土、草甸土、沼澤土）的室內光譜反射率為研究對象，結合數字高程模型（Digital Elevation Model，DEM）以及光譜指數作為輸入量，運用隨機森林算法（Random Forest，RF）進行SOM預測，得出以下結論：

1）通過競爭自適應重加權采樣（Competitive Adaptive Reweighted Sampling，CARS）算法，篩選出的特征波段不僅將輸入波段壓縮至全波段數目的16%以下，而且能夠在很大程度上降低土壤高光譜變量維度和計算復雜程度，從而提高了模型的預測能力。光譜變量經CARS算法篩選后模型調整后決定系數提高0.167，估測效果更好。說明CARS算法在提取特征關鍵波段變量、優化模型結構方面起到關鍵作用。

2）通過土壤分類進行SOM預測，不同土壤類型的SOM調整后決定系數存在差異，沼澤土的調整后決定系數最高為0.768，黑土次之，草甸土的預測精度最低，只有0.674，運用RF對3類土壤的SOM預測性能與四分位間隔距離的比率均大于1.8，說明無論是黑土、草甸土還是沼澤土，該模型都有一定的可信度，具有較好的預測能力。

3）通過將CARS篩選的特征波段、DEM以及光譜指數作為輸入量，運用RF模型，SOM的局部回歸模型驗證集精度最優，調整后決定系數為0.777，且RPIQ達到2.689，與全局回歸模型相比，模型的驗證精度提高了0.035。研究表明，3種類型的輸入量，進行單一土壤類型分別建模和全局回歸建模，其均具有較好的預測能力，在一定程度上可為以后不同土壤類型SOM預測時輸入量的選擇提供幫助，從而促進區域不同類型土壤進行SOM預測研究的進展，為農業和環境領域SOM的動態監測和建模提供理論支撐。