基于FOD和SVMDA-RF的土壤有機質含量高光譜預測

張智韜 勞聰聰 王海峰 ARNON Karnieli 陳俊英 李 宇

(1.西北農林科技大學水利與建筑工程學院， 陜西楊凌 712100；2.西北農林科技大學旱區農業水土工程教育部重點實驗室， 陜西楊凌 712100；3.本古里安大學Blaustein沙漠研究所， 思德博克 84990)

0 引言

土壤有機質含量(Soil organic matter，SOM)是評價土壤肥力的最主要指標，也是了解土壤退化狀態的主要指標，監測SOM是掌握土壤肥力和土壤退化情況的重要手段。高光譜遙感可以快速、無損監測土壤有機質含量，為荒漠土壤整治提供依據。然而土壤高光譜具有非特異性，波長相鄰的反射率互相重疊，且具有較強的相關性，當土壤肥力較低(SOM小于2%)時，土壤質地等特性會對光譜產生較大影響，在一定程度上掩蓋了SOM對光譜的作用，使得高光譜反演荒漠土壤的SOM具有一定困難[1]。

國內外學者對利用高光譜技術反演SOM已進行了大量研究[2-8]。但是上述研究均是以SOM很高或較高的地區為研究區域，其成果無法很好應用于SOM較低的干旱荒漠地區。侯艷軍等[1]利用高光譜反演新疆準噶爾盆地東部的荒漠土壤有機質，發現基于一階微分光譜的PLSR最優，其R2可達0.78。王海峰等[9]通過對以色列Sde Boker地區采集的荒漠化土壤進行研究得出，利用灰度關聯分析結合標準正態變換光譜建立的嶺回歸模型可以較好預測SOM，R2可達0.866。這說明高光譜遙感監測荒漠土壤SOM具有一定的可行性，但其所研究的模型均為線性模型，而線性模型無法較好解決非線性問題。因此，研究非線性組合模型反演荒漠土壤有機質含量十分必要。

此外，分數階微分[10](Fractional-order derivative，FOD)與一階、二階微分相比，能以較小微分間隔確保光譜曲線形狀緩慢變化，可以檢測到光譜信號的更多特征，而敏感光譜指數是與對應SOM經二維相關分析篩選確定，既考慮了波長間的關系，也考慮SOM與波長之間的關系，可減少無關波長的影響。根據HONG等[11]研究，將兩者耦合可進一步提升效果，但現有研究[11-12]一般取相關系數最高的指數建模，未完全利用經分數階微分優化的光譜指數。

本文以在以色列荒漠地區采集的砂質土與黏壤土為對象，在室內進行理化分析與光譜采集，利用平均光譜反射率建立支持向量機分類(Support vector machines discriminant analysis，SVMAD)模型，進行分數階微分后建立歸一化光譜指數(Normalized difference index，NDI)，再對NDI與SOM進行二維相關性分析，并以0階NDI的最高決定系數作為閾值篩選敏感指數，結合隨機森林模型(Random forest，RF)，構建SOM的SVMAD-RF模型，研究其對荒漠土壤的監測精度。

1 材料與方法

1.1 研究區概況

研究區為Sde Boker等地區(34°26′～34°47′E，30°17′～30°52′N)，該區域位于以色列南部，海拔640 m，屬于干旱半干旱氣候，年降水量介于100～200 mm之間，年均蒸發量達2 500 mm。尤其以色列南部全年降水稀少，占以色列一半以上的面積，被荒漠所覆蓋。故而此地的荒漠土壤占大部分且有機質含量低(0.5～18 g/kg)[13]，主要的土壤類型為黏壤土(CLS)與砂質土(SS)。圖1為取樣分布示意圖，其中圖1e所示樣本為砂質土，圖1b、1c、1d所示樣本均為黏壤土。

圖1 取樣點分布示意圖Fig.1 Sketch of study area and soil sampling locations

1.2 試驗材料及土壤樣本有機質含量的測定

通過十字法在研究區域土層(0～20 cm)采集134份土樣。土樣經過風干、磨細、經孔徑2 mm篩處理后，將每個土樣分成2份，分別用于光譜測試和理化性質分析。由激光粒度儀測得土樣粒度分布，根據國際土壤顆粒分類標準，將所取樣本劃分為兩類，其中砂質土樣本78個、黏壤土樣本56個。土樣SOM通過重鉻酸鉀外加熱法測定。

圖2 平滑后的土樣反射率曲線Fig.2 Reflectance curves of soil samples after smoothing

1.3 光譜測定

在照明控制的暗室，利用ASD FieldSpec 3型光譜儀測量土壤光譜反射率。將制備的134份土樣裝入直徑為10 cm、深度為2 cm的黑色容器中，土樣填滿容器后將表面刮平。光譜儀性能參數、光譜測試的幾何參數、光譜數據的采集方式等參照文獻[14]。

土樣光譜通過Savitzky-Golay算法(多項式階數為2，點數為10)進行平滑，并將噪聲較大的波段(350～399 nm、2 401～2 500 nm)剔除。為了降低冗余度，按10 nm為間隔重采樣，所得400～2 400 nm的光譜數據用于后續分析。去噪后的砂質土和黏壤土樣品光譜曲線如圖2(圖中灰色區域表示光譜的標準偏差，下同)所示。

1.3.1平均反射率

為提高對土壤質地的分類正確率，嘗試以平均反射率作為標準進行分類。全部波段(400～2 400 nm)的平均反射率計算公式為

(1)

式中Ri——波長i的土壤光譜反射率

圖3 不同有機質含量土壤平均反射率散點圖Fig.3 Scatter plots of average reflectance values of soils with different organic matter contents

1.3.2分數階微分

FOD是數學理論的分支，其將經典的整數階微分拓展至任意階微分。而Grünwald-Letnikov是一種離散形式的定義，便于進行數值計算，且運算效率較高[10]，故本文采用此表達式，具體為

(2)

式中 d——微分函數q——微分階數

h——步長t、a——微分的上限和下限

其中Gamma(Γ)函數[15]為

(3)

式中β——任意變量(本文指微分階數)

(4)

其中ν為1或2時，式(4)與一、二階微分的方程相同；而當ν=0時，表示未對光譜數據進行處理。本文利用式(4)分別計算兩種土壤的0～2階微分(間隔0.2階)。

1.3.3光譜指數

基于經分數階微分處理后的光譜數據，構建歸一化指數(Normalized difference index，NDI)，其計算式為

(5)

式中Aλ1、Aλ2—— 400～2 400 nm中任意兩波長，且Aλ1≠Aλ2

分析不同分數階微分下的NDI與SOM的相關性，設立閾值，確定估算SOM的敏感指數。

1.4 支持向量機分類和隨機森林模型

SVMDA具有強大的小樣本學習能力以及在高維空間中良好的擴展性能，因而廣泛用于分類[16]。RF是BREIMAN[17]提出的一個組合分類器算法，其主要思想是用自助法重采樣技術，基于原始樣本集，生成若干個自助樣本集，并將每個自助樣本集作為每棵分類樹的所有訓練數據。該模型是多個分類樹的集合，可以估計獨立變量之間復雜的非線性關系并且響應變量。有研究表明，RF預測性能較強，特別在高光譜技術估算土壤相關特征含量方面[18]。因而本文采用此法來預測SOM。

1.5 建模集和驗證集的劃分

SVMDA建模，以隨機抽樣法取89個樣本作為建模集，剩余45個樣本作為驗證集。RF建模采用Kennard-Stone(K-S)算法劃分建模集、驗證集,結果如圖4所示，圖中CV為變異系數，Mean為平均值，SD為標準差。

圖4 土壤有機質含量描述性統計特征箱線圖Fig.4 Box plots of density and descriptive statistics of SOM

1.6 模型建立與驗證

以土壤平均光譜反射率作為自變量，土壤質地類型作為因變量，采用非線性的SVMDA，建立土壤質地判別模型，以leave one out法進行交叉驗證。再對不同質地的土壤SOM分別建模，通過0階NDI的最高決定系數作為閾值獲取敏感指數，以所獲敏感指數作為自變量，SOM作為因變量，建立不同分數階的非線性RF模型。

1.7 數據分析軟件

Savitzky-Golay算法、平均反射率、分數階微分、光譜指數計算和篩選、K-S算法及SVMDA均在Matlab R2018a平臺進行。圖表的制作在OriginPro 2017C中實現。RF通過R3.5.2軟件中的Random Forest軟件包實現和優化，具體參數設置有分類樹數目(ntree)、訓練子集數目(mtry)。本文中ntree為500，mtry設置為1/3總變量數。

2 結果與分析

2.1 曲線特征分析

2.1.1土壤樣本高光譜曲線特征

圖5 不同土壤有機質含量光譜反射率變化曲線Fig.5 Reflectance spectra curves of soils with different SOM

圖5為2種土壤不同SOM的光譜反射率變化曲線，可以看出，在400～1 800 nm內，砂質土反射率曲線斜率大，變化快，黏壤土反射率曲線斜率較小，變化慢；土樣雖已風干，但殘留水分對光譜的作用依然顯著，在1 400、1 900、2 200 nm附近存在水分吸收谷[21]；在SOM基本相同時，黏壤土的反射率明顯小于砂質土；在全波段范圍內，就同一種土壤而言，隨著SOM的增加，光譜反射率逐漸降低。

2.1.2分數階微分光譜曲線特征

利用式(4)計算2種土壤0～2階分數階微分(間隔0.2階)，經預處理后的光譜如圖6、7所示。圖6a～6f、7a～7f顯示，當階數從0增加到1時，在580 nm處的負峰有所增加，此峰主要受SOM的影響，水分子振動作用逐漸凸顯，在1 400、1 900 nm處的吸收谷更加明顯[22]。由圖6d～6f、7d～7f可知，可見光區域的光譜吸收特征不斷增強，特別在430、550 nm處有2個正峰，而在480 nm處有1個負峰。由圖6e～6k、7e～7k可知，在0.8階以后，與水分有關的2 200 nm峰值較快增加，此處峰值也受SOM的影響，因為SOM越大，土壤吸附水分的能力就越強[23]。分數階微分光譜與原始光譜相比，可以改善光譜曲線分辨率，深度挖掘與SOM有關的潛在信息。

圖6 砂質土光譜反射率分數階微分曲線Fig.6 Spectral reflectance fractional differential curves of sandy soil

圖7 黏壤土光譜反射率分數階微分曲線Fig.7 Spectral reflectance fractional differential curves of clay loam soil

2.2 建模與分析

2.2.1土壤質地分類準確率分析

SVMDA對土壤質地的分類判別效果中，建模集和驗證集準確率都達到100%。這也表示土壤質地對土壤光譜反射率有顯著的影響。而平均反射率可以較好地表示砂質土和黏壤土對光譜響應的差異，用此模型可以精確地區分砂質土和黏壤土，為后續針對不同土壤質地建立SOM模型，并進行組合，提供了有效支撐。

2.2.2SOM與不同階數NDI的二維相關性分析

由圖8可知，砂質土決定系數高于黏壤土，這可能與土壤質地有較大關系，黏壤土的反射率較小，導致光譜中所含的SOM信息無法凸顯。

分數階微分對于反射率較小的黏壤土作用更加明顯，將決定系數從0.763提高到0.820。而對于反射率較大的砂質土作用較小，只將決定系數從0.901提高到0.917。在NDI指數中相關性較好的波段主要集中在可見光和短波近紅外區域(400～1 000 nm)(圖8)，這與HONG等[23]的研究結果一致。而砂質土、黏壤土的NDI指數與SOM的決定系數分別在0.8、0.6階時，達到最大值，但在1.4階以后，分數階微分對NDI的效果降低，這與WANG等[12]的研究結果基本一致。提取兩種土壤的每個分數階最佳組合波長，并計算NDI與SOM的決定系數，結果見表1。砂質土與黏壤土NDI的最佳組合波長都不相同，但砂質土的組合波長都稍大于黏壤土。在0.4～1.2階微分時，黏壤土最佳組合波長都集中在可見光區域。而砂質土的所有最佳組合波長中，均含有至少1個紅外波段。這表明黏壤土SOM對可見光區域更敏感，而砂質土SOM對紅外區域更敏感。

圖8 砂質土和黏壤土光譜反射率各階微分的歸一化光譜指數與有機質含量的決定系數Fig.8 Determination coefficients between SOM and normalized spectral index of sandy soil and clay loam soil reflectance treated by fractional order derivatives

2.2.3敏感光譜指數篩選

為獲得通過分數階微分所提升的敏感指數NDI，以0階最佳NDI所對應的決定系數作為閾值，篩選富含SOM信息的敏感指數。由表1可知，砂質土的篩選閾值為0.901，組合波長為R2210、R2170；黏壤土的篩選閾值為0.763，組合波長為R860、R850。

篩選結果如圖9所示，兩種土壤均在0.6階達到峰值，且黏壤土提升效果優于砂質土，差異的造成與土壤質地有較大關系。這說明分數階微分耦合光譜指數均能夠有效放大高光譜與SOM有關的潛在信息，尤其當階數為0.6時，可以獲得最多高相關性的自變量，以此建?？蛇M一步提高模型性能，且此法對提取被其他土壤屬性遮蔽的SOM信息十分有效，且對遮蔽作用越大(決定系數小)的土壤，其提升效果越明顯。

2.2.4敏感光譜指數建模

基于不同階數的敏感NDI所建RF模型參數見表2，表中N為模型自變量數目。為與敏感指數RF模型對比，建立不同土質的全波段RF模型(Full spectrum，FS)，其結果見表2。敏感NDI所建模型均具有較好的效果，均遠高于全波段模型。總體而言，利用敏感指數建模，可以一定程度上提高預測效果。RPD隨階數先增加后減小，但不同質地的土壤達到峰值的階數不同，砂質土1.2階模型(SS-1.2-NDI)的RPD最大，可達3.658；黏壤土0.6階模型(CLS-0.6-NDI)RPD最大，可達4.316。此差異可能與土壤質地有關。而基于全波段的不同土質FS模型中，砂質土模型較優，RPD為1.579，與敏感指數所建模型有較大差距。

將SS-1.2-NDI與CLS-0.6-NDI組合，建立SVMDA-RF模型，并計算其模型評價參數，結果見表3。為便于比較，建立將SS-FS、CLS-FS組合的SVMDA-RF模型，以及不區分土質(Unclassified samples，US)的全波段RF模型。從表3可以看出，SVMDA-RF(SS-1.2-NDI、CLS-0.6-NDI)效果最好， SVMDA-RF(SS-FS、CLS-FS)次之，RF(US)效果最差。且SVMDA-RF(SS-1.2-NDI、CLS-0.6-NDI)其RPD可達7.004，能精確預測SOM。而RF(US)的RPD只有2.228，只具有一定的預測能力。這表明區分土質分別建模并二次組合可以提高建模效果，結合敏感指數，能使模型具有較高精度，實現精確預測SOM。

SVMDA-RF(SS-1.2-NDI、CLS-0.6-NDI)的建立不僅可以快速區分土壤質地類型，還能夠針對不同質地的土壤精準預測土壤SOM。但如果土壤種類過多時，SVMDA區分土壤質地的正確率可能會下降，這將影響SVMDA-RF模型的穩定及其精度，因此需要尋找不同土質對光譜影響的特點，以此作為分類依據，為后續的RF模型精度提供保障。

高光譜數據具有高維度和數據冗余的特點，導致難以從中提取到與SOM相關性高的敏感數據。而且在SOM低于20 g/kg時，土壤中的其他組成物質(水分、鹽分、礦物等)和土壤質地的光譜反射將占主導，其中土壤含水率對光譜的作用較強。本文已將土樣經過風干處理、光譜重采樣和設置閾值剔除低相關性的指數，一定程度上消除了水分對光譜特征的影響。如果預測變量和響應變量之間存在非線性關系，則RF這類非線性模型通常會有較好的擬合效果，產生優異的估計精度。如WANG等[24]所建土壤鹽分的RF與偏最小二乘模型，RF優于線性的偏最小二乘模型，最佳模型的驗證效果較好，其RPD可達2.78。

圖10中，RF(US)模型的預測值與實測值在1∶1線上偏差較大，點的離散程度高。但利用SVMDAD對土質分類，再建模其精度得到明顯提升。由圖10c、10f可知，SVMDA-RF(SS-1.2-NDI、CLS-0.6-NDI) 的實測值與預測值的散點擬合線與1∶1直線基本重合，達到極高的精度，可以為高光譜遙感監測土壤荒漠程度提供較為可靠的依據。由圖10a、10b、10d、10e可知，黏壤土比砂質土的點更為分散，這與黏壤土的吸附特性優于砂質土有關。雖然經長時間的風干，但黏壤土可能還會吸附一些水分等物質。這些物質對于光譜的影響，更加掩蓋低含量SOM對于光譜的作用。故而有些研究者[1]會將1 400、1 900、2 200 nm附近存在的3個水分吸收谷作水汽吸收帶的剔除處理，以減少對有機質反演的影響。但是其附近可能含有SOM信息的波段，利用這些波長可能會提升有機質反演的精度。如本文砂質土的1.2階最佳NDI含有1 420 nm波長，然而利用其建模效果較好，RPD可達3.65。

3 討論

高光譜數據具有高維度和數據冗余的特點，導致難以從中提取到與SOM相關性高的敏感數據。而且在SOM低于20 g/kg時，土壤中的其他組成物質(水分、鹽分、礦物等)和土壤質地的光譜反射將占主導，其中土壤質地和黏粒含量是影響土壤光譜反射的重要因素之一[25]。故而為了減小土壤質地對定量反演SOM的影響，利用平均反射率建立SVMDA模型，其對所研究區域的土壤質地分類準確率可達100%。雖然以往研究表明利用光譜技術可以較好分辨不同的土壤質地，但是所建模型自變量數目多，較為復雜。如宋海燕等[26]通過正交信號校正(OSC)處理砂土、壤土、黏土3類土壤光譜，并結合偏最小二乘(PLS)對土壤質地分類，其正確率最高達93.3%，但無法很好區分壤土與砂土；WANG等[25]利用土壤吸光度對基于垂直干旱指數(PDI)分組的樣本，建立PLS分類模型，可以較好區分砂土與黏土，并預測黏粒含量。而本文利用平均反射率建立相對簡單的分類模型，且分類效果更好。這表明平均反射率用于砂質土和黏壤土分類是可行的。

表1 砂質土與黏壤土在不同階數下的最佳歸一化指數及其決定系數Tab.1 Determination coefficients between SOM and optimal normalized spectral index for sandy soil and clay loam under different fractional order derivatives

注：R2210表示波長為2 210 nm，以此類推。

圖9 敏感指數數量隨階數的變化曲線Fig.9 Changing curves of number of sensitive indices under different fractional order derivatives

表2 基于RF的SOM建模與預測結果Tab.2 Calibration and validation results of SOM based on RF

表3 分類組合模型與不分類模型的SOM建模與預測結果Tab.3 Calibration and validation results of SOM based on classification combination models and non-classification model

圖10 基于RF的SOM預測值與實測值的比較Fig.10 Comparison of measured and predicted SOM values based on RF

洪永勝等[27]分析基于4種不同光譜變換的歸一化光譜指數與SOM的二維相關性，發現R2最高達0.55。郭燕等[28]分析NDI、DI、RI與SOM的二維相關性，得到最大相關系數均在0.5～0.7之間。HONG等[23]分析基于6種不同光譜變換的光譜指數與SOM的二維相關性，得到最大相關系數為0.81。雖然二維相關分析能從二維光譜空間中更精準地捕捉到與SOM有關的光譜特征，有效減小土壤中其他組成物質(包括水分)對預測SOM的影響，但上述研究均未結合分數階微分，相關系數有待提高。而本文分析基于分數階微分的光譜指數與SOM的相關性，得到決定系數最高可達0.917。因此光譜數據先經過分數階微分處理，可以擴大二維光譜的空間，為獲取高相關性敏感指數提供了更大可能。且以0階決定系數篩選指數可以獲得更多特征信息，有利于提高建模的效果。

高光譜遙感反演的難點在于反演變量不完全是主導遙感光譜信息的因子，特別是在反演變量含量較低時，僅可以為遙感光譜提供微弱信號。而本文利用FOD聯合二維相關分析的方法可以得到高相關性的敏感指數。不過閾值的設立具有一定主觀性，需要更嚴謹客觀的方法進行篩選。今后的研究可以利用變量投影重要性分析法[30]、遺傳算法[31]等來篩選敏感指數。所建立的SVMDA-RF(SS-1.2-NDI、CLS-0.6-NDI)和SVMDA-RF(SS-FS、CLS-FS)都具有極好的預測性能，但是模型的普適性還需要大量其他地區荒漠土壤SOM數據來驗證。此外，本文只進行一種非線性模型(RF)的嘗試，以后的研究中可以建立不同非線性模型進行組合，比較獲得適應性更佳的模型。

4 結論

(1)基于平均反射率的SVMDA可以對土壤質地精確分類，將其與RF結合，可以提升荒漠土壤SOM的預測精度。

(2)0階NDI的最高決定系數作為閾值篩選敏感指數，可以獲得經分數階微分提升的NDI，其數量達到峰值的階數均為0.6階，但峰值與土壤質地有關。

(3)利用分數階微分結合NDI可獲得高相關性的敏感指數，以此建模可以極大提升模型的整體效果。