基于連續(xù)投影算法的博斯騰湖西岸湖濱綠洲土壤有機碳含量的高光譜估算

牛芳鵬，李新國*，麥麥提吐爾遜·艾則孜，趙慧

（1.新疆師范大學地理科學與旅游學院，烏魯木齊 830054；2.新疆干旱區(qū)湖泊環(huán)境與資源實驗室，烏魯木齊 830054）

土壤是碳在陸地生態(tài)系統(tǒng)中最大的儲存庫，并影響其生產力的發(fā)展，其中土壤有機碳（soil organic carbon,SOC）是構成土壤碳庫的重要部分[1]。干旱半干旱區(qū)是穩(wěn)定全球生態(tài)環(huán)境變化的重要環(huán)節(jié)之一，綠洲及其生態(tài)系統(tǒng)在對抗與緩解干旱區(qū)氣候方面起著重要作用，對土壤碳的研究特別是綠洲土壤有機碳的研究對全球碳循環(huán)及生態(tài)環(huán)境保護具有重要意義[2]。土壤光譜中包含豐富的土壤信息，便攜快速、無損與高精度的波譜特性使高光譜遙感技術在估算SOC含量的應用中，從定性到定量越來越成熟[3]。篩選具有代表性的光譜響應波段是構建高光譜估算模型的基礎，學者們在進行光譜數(shù)據特征波段的優(yōu)選方面已做了大量研究[4-5]。WANG 等研究表明，利用連續(xù)投影算法（successive projection algorithm,SPA）選擇的光譜特征波段可以有效地提高偏最小二乘回歸（partial least square regression,PLSR）模型的決定系數(shù)（coefficient of determination,R2）和殘余預測誤差（residual prediction deviation,RPD）性能[6]；王濤等采用相關性分析法和SPA 篩選特征波長，構建的最優(yōu)模型精度R2=0.98，均方根誤差（root mean square error,RMSE）=0.30 g/kg，RPD=9.32[7]；章海亮等采用遺傳算法結合SPA 挑選出18個特征波段并構建PLSR模型，其預測集的R2=0.83，RMSE=0.20 g/kg，RPD=2.45[8]。然而，利用土壤光譜反射率對SOC含量的估算，大多是通過建立線性方程來進行，這些傳統(tǒng)的回歸方法受數(shù)據本身影響較大；而支持向量機（support vector machine,SVM）的發(fā)展，有效解決了樣本維數(shù)問題[9]。張子鵬等通過構建SVM 模型，比較了不同維度光譜參數(shù)對SOC含量估測的準確性[10]；SUN 等對光譜數(shù)據進行SG（Savitzky-Golay）平滑與多種數(shù)據變換后，利用PLSR 方法構建了復墾土地SOC 含量估算模型，其驗證集R2=0.78，RMSE=1.81 g/kg，RPD=2.69[11]；HONG 等利用SVM 算法對華中地區(qū)土壤有機質（soil organic matter, SOM）含量進行估算，結果表明，用1.5 階微分構建的SVM 模型精度較高，R2達0.88，RPD=2.89[12]。

本研究以博斯騰湖西岸湖濱綠洲為研究區(qū)，通過測量采集的255 個樣本的SOC 含量與高光譜數(shù)據，并對光譜反射率進行SG平滑、標準化正態(tài)變換（standard normal variate, SNV）與一階微分（first derivative,1st Der）預處理，運用SPA 與相關系數(shù)法（correlation coefficient,CC）從原始光譜和預處理后的光譜中篩選特征波段，對比用PLSR 與SVM 2 種方法構建的SOC含量的高光譜估算模型的差異，為利用可見- 近紅外光譜（visible-near infrared spectroscopy,Vis-NIR）技術快速無損地估算湖濱綠洲SOC含量提供方法支持。

1 研究區(qū)概況與研究方法

1.1 研究區(qū)概況

博斯騰湖西岸湖濱綠洲位于新疆焉耆盆地，行政區(qū)隸屬于新疆博湖縣，地理位置為41°45′—42°10′N，86°15′—86°55′E，海拔1 047～1 048 m，地勢由西北向東南傾斜，為典型的山前湖泊綠洲。夏季月平均氣溫22.8 ℃，干旱少雨，蒸發(fā)強烈；冬季月平均氣溫9.0 ℃，寒冷干燥，無霜期176～200 d，多年平均降水量83.55 mm，年均溫8.0～8.6 ℃，季節(jié)過渡快，屬于大陸性荒漠氣候區(qū)；主要的自然植被類型有胡楊、檉柳和梭梭等；主要的土壤類型有綠洲潮土、草甸土、荒漠林土和沼澤土等[13]。

1.2 土樣采集與制備

綜合考慮到研究區(qū)范圍內的主要土地利用類型、植被覆蓋類型和微地形等因素，土樣采集按照“S”形線路隨機均勻布點，采樣時間為2019 年10月，共布設51 個典型土壤剖面，并進一步分割出5個不同深度（0～10、＞10～20、＞20～30、＞30～40和＞40～50 cm）土層進行取樣，去除雜物后采用四分法選取200 g土裝袋，共計255個樣品。帶回實驗室自然風干后進行研磨和過60目孔篩，一部分用作土壤高光譜數(shù)據的測試；另一部分用重鉻酸鉀-外加熱法測定SOC含量[14]。

1.3 土壤光譜的測定及預處理

光譜數(shù)據通過ASD FieldSpec3 地物光譜儀于室外采集，光譜波長為350～2 500 nm。選擇云量小于5%、無風或風力低于3 級的晴朗天氣，采集時間為12：00—14：00；光譜采集前先預熱儀器30 min并進行白板校正；光纖探頭視場角不超過25°，垂直于土樣表面15 cm處，使用五點梅花采樣法采集，在每個位置測量3條光譜曲線，共15條光譜曲線記錄存檔，每測完一組剖面土樣采集一次暗電流，同時進行白板優(yōu)化校正，減小誤差[15]。由于環(huán)境不可控因素，統(tǒng)一去除噪聲較大的尾部波段（2 451～2 500 nm）以及受環(huán)境水汽影響的1 300～1 450和1 800～1 950 nm 波段；此外，為減少高頻噪聲的影響，提升光譜數(shù)據信噪比，對原始光譜曲線進行窗口數(shù)為5的二次多項式SG平滑處理，并聯(lián)合使用SNV與1st Der 進行光譜數(shù)據預處理[16]。圖1 為預處理后的土壤光譜反射率曲線圖。

圖1 SG 平滑及SG-SNV-1st Der 預處理后的土壤光譜反射率曲線Fig.1 Soil spectral reflectance curve after SG smoothing and SG-SNV-1st Der pretreatments

1.4 連續(xù)投影算法的特征變量選擇

連續(xù)投影算法（SPA）是一種前向變量選擇算法，運用向量投影分析選擇最大向量，最終通過校正模型提取幾個特征波長[17]。其優(yōu)點是從光譜矩陣中選擇最小共線性的變量組合，從而降低模型的冗余度，提高模型的穩(wěn)定性和準確性。SPA 的具體步驟如下：

記xk(0)和N分別為初始迭代向量與需要提取的變量個數(shù)，光譜矩陣為J列。

1）任選光譜矩陣的1列（第j列），把建模集的第j列賦值給xj，記為xk(0)。2）將未選入的列向量位置的集合記為s，

3）用xj分別對剩余的列向量投影進行計算，

4）選取向量投影最大的光譜波長，

5）令xj=Px，j∈s。

6）n=n+1，如果n＜N，則按公式（1）循環(huán)計算。

最后，提取出的變量為{xk(n)=0，…，N-1}，分別構建關于每一次循環(huán)中k(0)和N的多元線性回歸（multiple linear regression, MLR）模型，得到建模集交互驗證的RMSE，以及對應的不同候選子集，其中最小的RMSE值對應的k(0)和N就是最優(yōu)值。

1.5 集合劃分與模型構建

采用基于聯(lián)合X-Y距離（sample set partitioning based on jointX-Ydistances, SPXY）方法進行建模樣本的劃分，該方法充分考慮了X和Y的空間可變性，計算樣本間的距離時能同時考慮SOC含量的影響與土壤光譜信息特征[18]。SPA、PLSR模型與SVM模型構建運用MatLab R2020a 軟件完成。PLSR 模型中采用“留一法”交叉驗證來決定最佳主成分數(shù)；SVM 模型選擇Poly 核函數(shù)進行分析。模型精度選用R2、RMSE、RPD 和四分位數(shù)間距性能比（ratio of performance to interquartile distance, RPIQ）進行評價，其中：R2的取值范圍為0～1.00，R2越大，RMSE越接近0，模型擬合度越高。當RPD＜1.40時，模型預測能力較差；當1.40≤RPD＜2.00 時，模型預測效果一般；當RPD≥2.00 時，模型預測能力較好[19]。RPIQ為樣本觀測值第三、四分位數(shù)Q3和第一、四分位數(shù)Q1的差與RMSE的比值，RPIQ越大，表示模擬結果分布的偏態(tài)程度越小[20]。

2 結果與分析

2.1 土壤有機碳含量描述

對255個樣本光譜反射率數(shù)據進行主成分分析后，剔除11 個異常樣本，運用SPXY 方法對剩余樣本進行訓練集與驗證集的劃分。由表1 可知：SOC質量分數(shù)變化范圍為0.75～48.13 g/kg，均值為13.31 g/kg，呈中等變異性，變異系數(shù)為63.19%；訓練集與驗證集的SOC質量分數(shù)平均值分別為12.67和12.16 g/kg，分別為右偏平頂峰與右偏尖頂峰，呈中等變異性，變異系數(shù)分別為52.33%和58.56%。

表1 土壤有機碳含量描述性統(tǒng)計Table 1 Descriptive statistics of soil organic carbon(SOC)content

2.2 土壤有機碳含量光譜曲線特征

根據研究區(qū)SOC含量狀況，并結合前人研究成果[21]，采用K-均值（K-means）聚類分析方法將SOC質量分數(shù)劃分為≤0.75、＞0.75～17.40、＞17.40～32.72、＞32.72 g/kg 4 類，圖2 為4 種不同SOC 含量的平均光譜曲線進行基線校正后的效果圖。從中可見，4 類光譜曲線形狀變化基本一致。在350～1 800 nm之間同一波段的不同SOC含量光譜反射率差異較明顯，SOC含量越高，土壤光譜反射率越低；在1 950～2 450 nm之間，光譜反射率差異較小。當SOC質量分數(shù)≤0.75 g/kg時，其反射率均值為0.38；當SOC質量分數(shù)＞0.75～17.40 g/kg時，其反射率均值為0.37；當SOC 質量分數(shù)＞17.40～32.72 g/kg 時，其反射率均值為0.35；當SOC質量分數(shù)＞32.72 g/kg時，其反射率均值為0.29。在350～569 nm 之間光譜反射率增幅較大，在570～1 299 nm 之間光譜反射率變化趨勢趨于平緩，在1 451～2 450 nm之間光譜反射率的波動較大；在923～955、1 109～1 141、2 148～2 240與2 333～2 358 nm波段內存在4個較明顯的吸收谷。

圖2 土壤有機碳含量的光譜反射率曲線特征Fig.2 Characteristics of spectral reflectance curve of soil organic carbon content

2.3 土壤有機碳含量與光譜曲線的相關性分析

由圖3 可知：SOC 含量與原始光譜反射率表現(xiàn)為負相關性，－0.62＜相關系數(shù)（r）＜－0.07，有1 160 個波段通過極顯著性檢驗（P＜0.01），主要集中在524～1 299、1 469～1 790 與1 973～2 056 nm之間，在661 nm 波長處相關性最高，相關系數(shù)∣r∣為0.62。經SG-SNV-1st Der預處理后，光譜反射率曲線呈正負波動，放大了原始光譜曲線的細微變化；通過極顯著性檢驗（P＜0.01）的波段數(shù)量縮減到414個，相關性較高波段主要集中在487～575、725～998 和1 464～1 514 nm 范圍內，在788、800 與1 768 nm波長處相關性最高，r均大于0.80。

圖3 土壤有機碳含量與光譜反射率的相關系數(shù)曲線Fig.3 Correlation coefficient curve between soil organic carbon content and spectral reflectance

2.4 基于SPA 對特征變量的篩選結果

由圖4 可知：利用SPA 對預處理后的光譜數(shù)據進行特征波長的篩選時，隨著篩選變量數(shù)量的增加，RMSE 先是迅速下降，當變量數(shù)為14 時，RMSE趨于穩(wěn)定狀態(tài)，其值為5.010 1 g/kg。通過SPA運算后得到14 個特征波長，僅占全光譜數(shù)據的0.78%，分別為399、1 011、1 046、1 061、1 073、1 596、1 632、1 667、1 749、2 012、2 103、2 268、2 305、2 341 nm，極大地縮減了光譜信息中的冗余變量。同理，運用SPA對原始光譜數(shù)據進行篩選，最終得到19個最優(yōu)特征變量，占全波段光譜數(shù)據的1.06%。

圖4 利用SPA篩選預處理光譜特征波長Fig.4 Screening of preprocessing spectral characteristic wavelength by SPA

2.5 基于PLSR 的土壤有機碳含量估算模型

由表2 可知，光譜預處理后，利用PLSR 模型獲得的RPD 均大于1.40，可以較好地估算SOC 含量。原始光譜通過SPA 構建的PLSR 模型驗證集R2為0.75，RMSE 為3.98 g/kg，獲得的RPD 為1.79，RPIQ為2.01，模型精度高于全波段（full-band）建模（R2=0.65）；通過相關系數(shù)法篩選出相關系數(shù)最高的5個波段進行建模，驗證集R2為0.70，RMSE 為6.14 g/kg，RPD 為1.16，RPIQ 為1.30，建模效果略好于全波段。光譜預處理后SPA模型精度最高，訓練集R2為0.79，RMSE為5.73 g/kg，驗證集R2為0.79，RMSE為3.58 g/kg，RPD為1.99，RPIQ為2.23。

表2 2種光譜模式的PLSR建模結果Table 2 PLSR modeling results of the two spectral modes

綜合比較2 種光譜模式下3 種變量的PLSR 模型精度，其模型估算能力表現(xiàn)為SPA＞相關系數(shù)法（CC）＞全波段法。由圖5可知：光譜預處理后構建的PLSR 模型較原始光譜PLSR 模型數(shù)據點更靠近1∶1 線；通過SPA 構建的PLSR 模型樣本的SOC 含量預測值小于實測值，大多較均勻地分布在1∶1 線下方。

圖5 2種光譜模式下PLSR估算模型結果圖Fig.5 Scatter plot of PLSR model of full-band and characteristic bands under the two spectral modes

2.6 基于SVM 的土壤有機碳含量估算模型

由表3 可知，基于全波段光譜構建的2 種SVM模型，其驗證集R2分別為0.68和0.76，RMSE分別為4.21 和4.26 g/kg，RPD 分別為1.69 和1.67，RPIQ 分別為1.90 和1.88，可粗略對樣本進行預測。通過相關系數(shù)法（CC）與SPA 進行特征變量篩選后構建的SVM 模型的RPD 均大于2.00，擬合程度較好，模型效果有明顯提高。原始光譜基于相關系數(shù)法構建的SVM模型訓練集和驗證集R2分別為0.69和0.70，RPD 為2.00，RPIQ 為2.25；光譜預處理后基于相關系數(shù)法構建的SVM 模型訓練集和驗證集R2分別為0.80 和0.77，RPD 為2.13，RPIQ 為2.39。原始光譜與預處理光譜基于SPA 構建的模型驗證集R2分別為0.73 和0.81，RMSE 分別為3.35 和3.16 g/kg，RPD分別為2.13和2.25，RPIQ分別為2.39和2.53。綜合比較3種不同變量構建的SVM模型效果，對SOC含量估算的效果依次為SPA＞相關系數(shù)法＞全波段法。

表3 2種光譜模式的SVM建模結果Table 3 SVM modeling results of the two spectral modes

由圖6 可知：光譜預處理后基于SPA 與相關系數(shù)法構建的SOC 含量的SVM 模型估算結果較好，樣本均勻地接近1∶1線；全波段SVM模型及原始光譜SPA 模型驗證樣本逐漸偏離1∶1 線，驗證樣本數(shù)據點分布比較分散，模型預測效果較差，與前文分析一致。

圖6 2種光譜模式下SVM估算模型結果圖Fig.6 Scatter plot of SVM model of full-band and characteristic bands under the two spectral modes

3 討論

關于室內與室外土壤光譜反射率的采集方法始終是土壤研究的主要課題之一。室內測試是為了研究土壤中的某個因子對光譜反射率特性的影響，室外測試能較好地反映自然景觀的真實性，描述表面反射特征，以便為航空和航天傳感器定標[22-23]。徐彬彬等認為，在土壤光譜反射特性研究中，應當注重野外實測，盡管在野外測試中受到當前所用儀器的限制和環(huán)境條件的影響，但它還是能較好地反映自然界的部分真實情況[24]；馬利芳等于野外采集高光譜數(shù)據后，研究了新疆阜康市土壤鹽分離子的高光譜特征，為區(qū)域尺度的土壤鹽分主要離子含量估算提供了良好的支撐[25]。本研究采用室外光譜采集的方法，使獲得的數(shù)據能更好地接近研究區(qū)自然環(huán)境條件下土壤有機碳高光譜信息。

原始光譜經SG-SNV-1st Der 預處理后，通過全波段、相關系數(shù)法與SPA 構建的PLSR 模型驗證集R2分別提高了6.15%、5.71%和5.33%；SVM 模型驗證集R2分別提高了11.76%、10.00%和10.96%。這與張子鵬等[16]和李冠穩(wěn)等[26]的研究結果基本一致。在光譜建模之前對光譜數(shù)據進行預處理，既能突出光譜的特征波段，還可以提高模型的擬合效果。本研究運用相關系數(shù)法與SPA分別篩選出5和19個特征波段，僅占全波段數(shù)據的0.28%和1.06%，驗證集R2高于全波段建模1.31%和6.58%。這與韓建等[27]和VISCARRA ROSSEL等[28]的研究結果相一致。SPA 可以有效地消除波段之間的共線性影響，剔除不相關變量，降低數(shù)據冗余度，提高模型的建模精度。光譜預處理后，基于3 種變量方法構建的SVM 模型較PLSR 模型的驗證集R2分別提高了10.14%、4.05%和2.53%，說明SVM模型能在一定程度上彌補PLSR模型在解決非線性關系問題上的缺陷。這與楊愛霞等[29]、曾胤等[30]的研究結果相一致。本文構建的SOC 含量高光譜估算模型尚未考慮土地利用類型、土壤結構和土壤水分等因素，其對模型的影響還需進一步驗證。

4 結論

1）研究區(qū)SOC 質量分數(shù)變化范圍為0.75～48.13 g/kg，平均值為13.31 g/kg，呈中等變異性，變異系數(shù)為63.19%；同一波段內，隨著土層深度的增加，光譜反射率越高，且SOC含量越高，土壤光譜反射率越低。當SOC 質量分數(shù)≤0.75 g/kg 時，其反射率均值為0.38，當SOC 質量分數(shù)＞32.72 g/kg 時，其反射率均值為0.29。

2）通過光譜變換可以明顯提高相關系數(shù)，SOC含量與原始光譜反射率呈負相關性，－0.62＜相關系數(shù)（r）＜－0.07；經SG-SNV-1st Der 預處理后，通過極顯著性檢驗（P＜0.01）的波段數(shù)達到414個，主要集中在487～575、725～998和1 464～1 514 nm范圍內，在788、800 與1 768 nm 波長處的相關系數(shù)均大于0.80。

3）SPA 的降維效果優(yōu)于相關系數(shù)法，光譜建模時SNV 非線性方法的效果優(yōu)于PLSR 線性回歸模型。光譜經SG-SNV-1st Der預處理后，運用SPA結合SVM 模型能很好地估算研究區(qū)SOC 含量，其訓練集與驗證集R2分別為0.79 和0.81，RMSE 分別為5.61和3.16 g/kg，RPD為2.25，RPIQ為2.53。