云南山地膠園土壤有機質高光譜估算

陳桂良，劉忠妹，許木果，黎小清，丁華平，楊春霞

（云南省熱帶作物科學研究所，云南景洪 666100）

0 引言

云南是中國種植面積最大、產膠最多、單產最高的優(yōu)質天然橡膠生產基地[1]。云南橡膠園為典型的山地膠園，地形地貌復雜，小氣候復雜多變，成土母質豐富，土壤養(yǎng)分空間差異明顯[2-3]，植膠初期土壤肥力較高，但在長期植膠后，土壤退化嚴重[4-5]。土壤有機質是反映橡膠園土壤養(yǎng)分狀況的重要指標，其含量變化受到廣泛關注[6-8]。傳統(tǒng)化學分析方法存在周期長、成本高、污染環(huán)境等不足，而高光譜技術憑借快速、簡便、無污染等優(yōu)勢，為土壤有機質含量的快速測定開辟了新的途徑[9-10]。土壤有機質含量快速、高精度估算模型的建立，可以更好地指導云南山地膠園精細化生產管理，對云南天然橡膠產業(yè)的高質量發(fā)展意義重大。

國內外學者普遍認為土壤光譜信息中存在有機質敏感波段，利用土壤有機質敏感波長建立土壤有機質含量估算模型，可以簡化模型，保持甚至提升模型精度[11-13]。為突出土壤有機質光譜信號，通常會對原始光譜反射率進行適當變換，然后進行特征波長優(yōu)選，再建立相應的預測模型，以提高模型的穩(wěn)定性和精準度[14-17]。有關光譜反演土壤有機質含量的研究中，模型的建立大多數(shù)以經驗統(tǒng)計方法為主，由于不同區(qū)域土壤類型多樣且成因復雜，使得高光譜特征和估算模型的差異很大，難以建立統(tǒng)一的有機質估算模型[18]。為了提高模型的有效性，針對特定區(qū)域建立區(qū)域化的有機質光譜預測模型是常用手段[19]。很多學者對國內典型區(qū)域進行了土壤有機質高光譜估算研究，如東北黑土[20]、干旱區(qū)典型綠洲[21]、紅壤區(qū)[22]、華南地區(qū)[23]、黃河三角洲[24]等。目前還未見利用高光譜技術估算云南山地膠園土壤有機質的研究報道。

筆者以云南典型山地膠園東風農場為研究區(qū)，基于獲取的土壤有機質和光譜反射率數(shù)據，在優(yōu)選光譜變換模式和特征波長的基礎上，采用多元線性回歸、偏最小二乘回歸和支持向量回歸，構建云南山地膠園土壤有機質含量的高光譜最優(yōu)估算模型，以期實現(xiàn)云南山地膠園土壤有機質含量快速檢測。

1 材料與方法

1.1 研究區(qū)與土壤樣品采集

云南省景洪市東風農場位于景洪市西南部，在21°30′00″—21°46′12″N、100°34′48″—100°49′48″E 之間，海拔600～900 m 之間，成土母質主要是花崗巖、千枚巖、砂頁巖和老沖積物，土壤主要為磚紅壤和磚紅壤性紅壤，下轄6 個分場，現(xiàn)有天然橡膠1.02 萬hm2。2019 年9 月，在東風農場選取了103 個GPS 定位采樣點，根據地形的差異，在每個GPS定位采樣點采集1～3個橡膠樹保護帶土壤樣品(0～20 cm)，共采集到土壤樣品225個。圖1為研究區(qū)位置及采樣點分布圖。

圖1 研究區(qū)位置及采樣點分布圖

1.2 土壤有機質含量與光譜反射率測定

土壤樣品帶回實驗室后，將每個土壤樣品混合均勻，剔除樹根、石塊等雜質，經室內自然風干后研磨并過0.15 mm 孔徑篩，用于土壤有機質和光譜反射率測定。采用重鉻酸鉀容量法測定土壤有機質[25]。采用FieldSpec4光譜儀（美國ASD公司產）進行土壤光譜反射率測定，整個測定工作在一個同等于暗室的實驗室進行。用培養(yǎng)皿（半徑10 cm、深1.5 cm）裝滿土壤樣品，并將表面刮平，放在反射率近似為0的黑色橡膠墊上。光源為50 W鹵素燈，光源入射角度為45°，距土壤樣品表面70 cm。傳感器探頭采用25°視場角，位于土壤樣品表面垂直上方15 cm處。光譜測定之前先進行白板校正，測定期間每隔15 min 進行白板校正，每個土壤樣品采集10條光譜曲線，算術平均后得到該土壤樣品的光譜反射率數(shù)據。

1.3 研究方法

1.3.1 光譜預處理及建模集劃分 剔除噪聲較大的350～399 nm波段，剩余的400～2500 nm光譜反射率數(shù)據用于本研究的建模研究。將光譜反射率重采樣為5 nm間隔(400、405、410、…、2495、2500 nm)。 采用Kennard-Stone(KS)算法對225 個土壤樣本進行優(yōu)選，選出150個作為校正集，剩余75個作為驗證集，校正集和驗證集的土壤樣品有機質含量統(tǒng)計見表1。

表1 用于模型校正和驗證的土壤樣品有機質含量統(tǒng)計

1.3.2 光譜變換為了消除光譜數(shù)據無關信息和噪聲，提高建模精度，首先采用3 種方法對重采樣后的光譜反射率R進行光譜變換處理，即倒數(shù)對數(shù)變換[log(1/R)]、多元散射校正(MSC)、標準正態(tài)變換(SNV)，然后對光譜反射率R以及3 種變換形式光譜數(shù)據進行Savitzky-Golay(SG)平滑或導數(shù)變換，以得到最佳的光譜變換模式。

1.3.3 特征波長篩選方法

（1）競爭適應重加權采樣。競爭適應重加權采樣(competitive adaptive reweighted sampling，CARS)，通過優(yōu)選偏最小二乘回歸模型回歸系數(shù)絕對值相對較大的波長變量，剔除模型回歸系數(shù)絕對值相對較小的波長變量，獲得模型交叉驗證均方根誤差最小的優(yōu)選變量子集[26-27]。

（2）連續(xù)投影算法。連續(xù)投影算法(successive projections algorithm，SPA)，是在光譜矩陣中應用變量投影操作尋找含有冗余信息最低、共線性最少的光譜特征變量組，最大程度地避免光譜信息重疊、簡化模型結構、提高建模的速度和效率[28]。

（3）CARS-SPA。將CARS與SPA 2種方法結合用于特征波長篩選，首先用CARS 方法獲得優(yōu)選的變量子集，然后再應用SPA 方法對CARS 方法獲得的變量子集進一步優(yōu)選。

1.3.4 建模方法

（1）多元線性回歸。多元線性回歸(multiple linear regression，MLR)是由多個自變量的最優(yōu)組合通過回歸預測因變量的一種回歸方法，在已知多組自變量和所對應的單組因變量時，可以對它們之間的關系進行很好的擬合和預測[29]。

（2）偏最小二乘回歸。偏最小二乘回歸(partial least squares regression，PLSR)廣泛用于高光譜數(shù)據回歸建模，是建立穩(wěn)健線性光譜定量校正模型的通用方法。PLSR 集主成分分析、線性回歸分析和典型相關分析的優(yōu)點于一身，適于處理自變量存在多重共線性的回歸問題[30]。

（3）支持向量回歸。支持向量回歸(support vector regression，SVR) 是支持向量機(support vector machine，SVM)的重要應用分支，用于回歸模型的構建。SVM 通過核函數(shù)將低維非線性問題轉換成高維的線性問題，通過結構風險最小化原則提高泛化能力，在保證最小化樣本的同時縮小了模型泛化誤差的上界，廣泛應用于分類與回歸[31-33]。

1.3.5 模型驗證模型驗證采用決定系數(shù)(R2)、均方根誤差(RMSE)和相對分析誤差(RPD)，R2越接近1，RMSE越小，RPD越大，說明模型預測效果越好。1.5＜RPD＜2 表明模型只能對樣品進行粗略估算，2.0＜RPD＜2.5表明模型具有較好的估算能力，2.5＜RPD＜3.0 表明模型具有很好的估算能力，RPD＞3.0 表明模型具有極好的估算能力[18]。

2 結果與分析

2.1 土壤有機質含量與光譜反射率

參照GB/T 29570—2013 附錄C 膠園土壤養(yǎng)分含量正常指標[34]，將所有樣本按有機質含量高低劃分為高(＞25 g/kg)、中(20～25 g/kg)、低(＜20 g/kg)3 種有機質含量水平，并計算每種有機質含量水平土壤樣品的平均反射率光譜，發(fā)現(xiàn)400～2500 nm 波段范圍內反射率與有機質含量負相關，且不同有機質含量水平土壤的反射率光譜曲線有相似的反射特征，總體呈現(xiàn)先增加后降低的拋物線型（圖2）。在可見光波段，受土壤發(fā)色團和有機質本身黑色的影響，存在較寬的吸收波段，反射率總體較低，隨著波長的增加而急劇上升；在900 nm 附近出現(xiàn)典型的Fe3+的吸收谷；在1400、1900、2200 nm 附近有3 個明顯的吸收峰，深度略有差別，這可能是黏土礦物中含有的水分子和羥基的吸收帶[35]。

圖2 3種有機質含量水平土壤的平均原始反射率光譜

2.2 光譜變換模式篩選

對重采樣后的光譜反射率R進行光譜變換處理后，對光譜反射率R以及3 種變換形式光譜數(shù)據進行SG 平滑或導數(shù)變換模式優(yōu)選。優(yōu)選方法如下：基于SG 平滑或導數(shù)變換后的全波段光譜數(shù)據與有機質含量數(shù)據，采用PLSR 及留一交叉驗證方法建模，其中，SG 平滑或導數(shù)變換的濾波窗口在3～101 的奇數(shù)中篩選，多項式次數(shù)1～9 中篩選，導數(shù)階數(shù)在0、1、2 中優(yōu)選，PLSR主成分數(shù)在1～20中篩選，按照交叉驗證均方根誤差(root mean square error of cross-validation，RMSECV)選擇最佳的光譜變換模式。

log(1/R)、MSC 和SNV 變換在Unscrambler 9.7 中實現(xiàn)，SG 平滑或導數(shù)變換和PLSR 模型的建立通過MATLAB 2015 編程實現(xiàn)。研究發(fā)現(xiàn)MSC 和SNV 變換對模型影響不大，log(1/R)變換能夠顯著提升模型性能，log(1/R)變換后再進行SG平滑或導數(shù)變換，可以進一步提升模型預測能力，但不同的導數(shù)階數(shù)、濾波窗口大小對模型的預測能力影響較大（圖3）。log(1/R)結合SG平滑變換是最佳的光譜變換模式，其中，SG平滑變換模式為導數(shù)階數(shù)0，SG 濾波窗口5，多項式次數(shù)2 或3（表2）。

表2 不同光譜變換和不同導數(shù)階數(shù)下最優(yōu)模型預測效果

圖3 不同光譜變換下各階導數(shù)不同濾波窗口大小對應的最小交叉驗證均方根誤差

采用PLSR 模型和留一交叉驗證法，分別基于全波段的光譜反射率R和最佳光譜變換光譜數(shù)據，建立土壤有機質高光譜估算模型，利用驗證集對模型進行驗證，相比光譜反射率R，基于最佳光譜變換光譜數(shù)據建立的模型精度顯著提高，RPD從2.433 增加到2.795（圖4）。表明應用本研究的最佳光譜變換模式對光譜反射率R進行光譜變換處理，能夠有效突出土壤有機質光譜信號，有利于提高模型的穩(wěn)定性和精準度。

圖4 基于全波段的光譜反射率R(a)和最佳光譜變換光譜數(shù)據(b)建立的PLSR模型預測效果

2.3 特征波長篩選

為進一步提高模型的穩(wěn)定性和精準度，基于篩選得到的最佳光譜變換光譜數(shù)據與土壤有機質含量數(shù)據，選擇CARS[36]、SPA[37]、CARS-SPA 3 種特征波長篩選方法，提取得到3個特征波長組合（表3）。

表3 不同特征波長篩選方法提取的特征波長

2.4 模型構建與驗證

基于校正集的土壤有機質含量與最佳光譜變換模式下的光譜數(shù)據，分別選用表3所示的3種特征波長組合對應的光譜特征變量作為自變量，土壤有機質含量作為因變量，并采用MLR、PLSR和SVR 3種方法構建土壤有機質高光譜估算模型。其中，MLR 和PLSR 模型在Unscrambler 9.7 中實現(xiàn)，結合留一交叉驗證方法進行建模。SVR 模型利用LIBSVM 軟件[38]實現(xiàn)，選取SVR 類型為epsilon-SVR，核函數(shù)類型為radial basis function，通過網格搜索和留一交叉驗證進行參數(shù)優(yōu)選。利用驗證集對所建模型進行精度驗證，模型的預測效果如圖5 所示。從模型預測效果來看，CARSSVR最優(yōu)。

圖5 不同模型的預測效果

3 結論

本研究嘗試對土壤光譜反射率R進行多種光譜變換后，再結合SG平滑或導數(shù)變換模式篩選，得到最佳的光譜變換模式為log(1/R)結合SG 平滑變換，其中，SG 平滑變換模式為導數(shù)階數(shù)0、SG 濾波窗口5、多項式次數(shù)2 或3。基于最佳光譜變換光譜數(shù)據與土壤有機質含量數(shù)據，選擇CARS、SPA、CARS-SPA 3 種方法進行特征波長篩選，并采用MLR、PLSR和SVR 3種方法構建土壤有機質高光譜估算模型。結果表明，MLR、PLSR、SVR 3 種模型類型中，最佳模型分別為CARS-MLR、SPA-PLSR、CARS-SVR，RPD分別為2.745、2.617、2.947。CARS-SVR 模型估算精度最高，R2、RMSE、RPD分別為0.897、3.990 g/kg、2.947，該模型RPD位于2.5～3.0之間，具有很好的估算能力。本研究應用高光譜技術，對云南省景洪市東風農場典型山地膠園土壤有機質含量進行估算，為云南山地膠園土壤有機質含量的快速估算提供了參考。

4 討論

本研究采用3種方法[log(1/R)、MSC、SNV]對5 nm間隔光譜反射率R(400、405、…、2495、2500 nm)進行光譜變換處理，并對光譜反射率R以及3 種變換形式光譜數(shù)據進行SG 平滑或導數(shù)變換模式優(yōu)選，得到最佳的光譜變換模式為log(1/R)結合SG 平滑變換，其中，SG 平滑變換模式為導數(shù)階數(shù)0、SG 濾波窗口5、多項式次數(shù)2 或3。研究發(fā)現(xiàn)，log(1/R)變換能夠明顯提升土壤有機質估算模型性能，這與方少文等[22]的研究一致。log(1/R)變換后再進行SG 平滑或導數(shù)變換，可以進一步提升模型預測能力，但不同的導數(shù)階數(shù)、濾波窗口大小對模型的預測能力影響較大。

基于全波段的最佳光譜變換光譜數(shù)據建立土壤有機質PLSR模型，并利用驗證集對模型進行驗證，RPD達2.795?；诤Y選得到的最佳光譜變換光譜數(shù)據與土壤有機質含量數(shù)據，選擇CARS、SPA、CARS-SPA 3種方法進行特征波長篩選，并采用MLR、PLSR 和SVR 3 種方法構建土壤有機質高光譜估算模型，利用驗證集對各模型進行驗證，結果表明，MLR、PLSR、SVR 3 種模型類型中，最佳模型分別為CARS-MLR、SPA-PLSR、CARS-SVR，RPD分別為2.745、2.617、2.947。研究發(fā)現(xiàn)，PLSR 作為一種建立穩(wěn)健線性光譜定量校正模型的通用方法，在全波段模型中具有一定優(yōu)勢，但在特征波段模型中，不如SVR，甚至不如MLR，與王濤等[39]的研究結果一致，這可能也與特征波長提取方法有關。SVR作為一種非線性建模方法，特別是基于特征波長建模，在土壤有機質高光譜估算研究中具有較大潛力。

為了盡可能地提高光譜數(shù)據與土壤有機質含量之間的相關性，有效突出土壤有機質光譜信號，本研究嘗試對土壤光譜反射率R進行多種光譜變換后，再結合SG平滑或導數(shù)變換模式篩選，這為土壤有機質含量高光譜估算模型的建立提供了新思路。本研究采集的土壤樣品包括花崗巖、千枚巖、砂頁巖和老沖積物等成土母質，也是云南山地膠園主要母質類型[5]。由于建立的有機質估算模型主要針對云南山地膠園，為了模型的區(qū)域普適性，并未按母質分類建模，而是將全部土壤樣品混合建模，但仍然取得了較好的效果。本研究僅用同一農場的土壤樣本數(shù)據對模型進行了驗證，在后續(xù)研究中，還需在云南山地膠園開展更為廣泛的模型驗證及優(yōu)化。