張俊華 尚天浩 陳睿華 王怡婧 丁啟東 李小林
(1.寧夏大學生態環境學院, 銀川 750021;2.西北土地退化與生態恢復國家重點實驗室培育基地, 銀川 750021;3.西安煤航遙感信息有限公司,西安 710199;4.寧夏大學地理科學與規劃學院,銀川 750021)
土壤有機質(Soil organic matter,SOM)是植物營養物質的主要來源之一,在土壤肥力質量評價中具有重要地位[1]。SOM成分復雜,光譜特性差異大,但總體呈現出SOM含量升高,整個譜線反射率逐漸降低的趨勢[2]。但當SOM含量較低,光譜波段被吸收部分的能量過低時,光譜特性分析效果不佳[3]。ABERGAZ等[4]指出,SOM含量(質量比)小于20 g/kg時,光譜曲線易受其他成分的影響,不能準確反映由官能團引起的吸收特征。因此,如何提高低含量SOM的反演精度是相關領域的研究熱點之一。
利用不同的反射率變換方式篩選SOM敏感波段和采用不同建模方法都是提高低含量SOM反演模型精度的有效手段。XIE等[5]研究指出,與反射率的一階微分、倒數的一階微分、平方根的一階微分和對數的一階微分變換方式相比,特征波段選擇一階微分可以顯著提高低含量SOM的反演精度。LIU等[6]在一階微分、二階微分和吸光度對數、吸光度對數的一階和二階微分及連續統去除等變換方式的基礎上,采用競爭自適應加權算法篩選出15~40個不等的特征波段,發現采用連續統去除-隨機森林建立反演陜西省靖邊縣SOM模型精度達0.96,相對分析誤差(Residual predictive derivation,RPD)為3.02。在建模方法中,最小二乘法(Least square,LS)[7-8]、支持向量機(Support vector machine,SVM)[7-8]、支持向量機分類(Support vector machine classification,SVMC)[9]、遺傳算法(Genetic algorithm,GA)[10]、最小二乘回歸(Least square regression,LSR)[11]、反向傳播神經網絡(Back propagation neural network,BPNN)[12-13]、隨機森林(Random forest,RF)[14-15]等一種或多種方法被用于巴西、以色列荒漠區以及我國河北、新疆、西藏、寧夏等地區低含量SOM的反演,都取得了較好的效果。此外,引入其他變量也可以提高SOM反演精度[13,16-17]。雖然變量的引入可以顯著提高SOM含量的反演精度,但這些變量的獲取技術難度較大,所以其應用和推廣有一定困難。
為提高模型運算效率,減少建模變量輸入個數,有學者在敏感波段篩選的基礎上,用兩個或兩個以上特征波段分別進行波段耦合運算,構建比值指數(Ratio index,RI)、差值指數(Difference index,DI)、土壤沙化指數(Soil desertification index,SDI1)、土壤退化指數(Soil degradation index,SDI2)、修正歸一化差分指數(Corrected normalized difference index,CNDI)等[18-20],也有學者嘗試在二維空間建立光譜數據與SOM間的全波段相關性運算,以解決單一敏感波段在一維空間構建光譜指數所造成的光譜信息缺失問題[21-22]。二維空間所建模型精度較一維空間所建模型精度得到顯著提升,但因研究區域不同,導致所選最佳建模光譜指數存在較大差異。綜合來看,已有研究或基于不同二維光譜指數分別建模,或以分數階微分(Fractional-order derivarives,FOD)、整數階微分聯合單一光譜指數建模,將FOD聯合優化指數進行較低SOM含量反演的報道相對較少。

銀川平原(37°50′~39°23′N,104°17′~107°39′E)地處寧夏北部,是我國糧食主產區之一,屬溫帶干旱氣候,年均降水量約190 mm,全年蒸發量約 1 750 mm。土壤類型多以灌淤土、灰鈣土和鹽堿土為主[24],SOM含量普遍偏低。
根據銀川平原氣候、水文、地質地貌和土壤等基礎資料數據,主要以農田土壤為研究對象,利用ArcGIS 10.4預設5 km×5 km 網格。于2019年5月以梅花五點法采集裸露地表0~20 cm的土樣(以光譜測定點為中心,半徑為1 m距離4個方向處各選1個點,共5個點),混合均勻后保留0.5 kg作為該點土樣。采樣過程中受建筑、道路和水體的影響,適當調整采樣點(圖1)。土樣風干后采用外加熱-重鉻酸鉀法測定SOM含量,共獲得有效土壤樣本187個。

圖1 研究區及樣點分布
采用ASD FieldSpc4型地物光譜儀進行野外原位土壤的高光譜測定,波段為350~2 500 nm,每次測定前均進行白板校正,光譜探頭始終保持在距地面約80 cm,每一樣點重復測定5次,取平均值作為該樣點的光譜反射率。測定時間為10:00—14:00,天氣狀況良好,晴朗無風。
為消除儀器噪聲和環境背景干擾,去除噪聲過大的邊緣波段(350~399 nm和2 401~2 500 nm)。在Matlab 2018b軟件中對400~2 400 nm范圍內的高光譜數據作S-G平滑去噪處理。FOD作為一種光譜信號處理技術[22],通過細化階數改變光譜曲線中的斜率和曲率,使不同分數階下光譜曲線間的特征波段處微弱差異得到明顯放大[25]。本研究以整數階微分定義推廣出的G-L算法進行分數階微分計算,計算式為


(1)
式中n——上下限微分階數之差
f(λ)——以波長λ處反射率為自變量的函數
v——階數
本研究分數階設置為0~2階,以0.20階為間隔,共11階。
基于高光譜遙感信息進行表層土壤屬性快速反演,其中光譜信息多以連續窄波段為主。本文在全波段400~2 400 nm范圍,通過對差值指數、比值指數、歸一化指數(Normalized difference index,NDI)、再歸一化差值指數(Renormalized difference index,RDI)和廣義差值指數(Generalized difference index,GDI)進行兩兩算法優化(表1),建立全波段(400~2 400 nm)與SOM含量間的二維相關性圖。

表1 優化光譜指數及公式
SVM模型基于結構風險最小化處理,能夠很好地解決模型中存在的局部最優問題[26]。本研究選用決定系數R2、均方根誤差(Root mean square error,RMSE)和RPD作為所建模型的精度和預測效果評價。一般將RPD分為3類來評價模型的可靠性:RPD小于等于1.4,表明模型預測能力極差;RPD大于1.4、小于等于2.0,表明模型能夠具有較為可靠預測能力;RPD大于2.0,表明模型預測能力極好[23]。
根據全國第二次土壤普查推薦的土壤肥力分級標準,將銀川平原187個土樣按SOM含量共分為6個等級(表2):SOM含量為2.81~56.14 g/kg,均值12.48 g/kg、變異系數57.26%,其中處于四級~六級(SOM含量小于20 g/kg)的土樣有174個,占總數的93.05%,該范圍SOM含量平均值為11.06 g/kg,變異系數為39.30%,本文選取這174個土樣為研究對象。

表2 銀川平原土壤SOM統計特征
利用K-S算法[27]中的歐氏距離劃分總樣本集,其中建模和驗證樣本分別為117和57(表3)[23],所劃分的建模數據和驗證數據較總樣本SOM含量區間值和標準差較為一致,變異系數說明SOM均屬中度變異。

表3 建模與驗證土壤樣本有機質含量統計特征
由圖2可知,原始光譜反射率吸收帶整體較寬,吸收特征差異明顯,其中1 400、1 900 nm處存在明顯吸收峰。在400~1 800 nm范圍內,光譜反射率隨波長的增加而增大;在1 900~2 200 nm范圍內,光譜重疊嚴重。隨著分數階的不斷增加,重疊峰和基線漂移逐漸消除,光譜強度的幅度下降,光譜反射率不斷趨近于0,其中1 900 nm處的吸收谷逐漸由正峰變為負峰;近紅外波段玫紅色區域的光譜標準差差異明顯大于可見光區域。當FOD從0階增加到1.6階時,土壤光譜反射率逐漸降低并最終穩定在-0.07~0.07;當從1.6階增加到2.0階時,土壤光譜反射率在-0.06~0.06穩定不變,且高分數階微分的吸收谷明顯小于低分數階微分。

圖2 土壤反射率分數階微分曲線
比較土壤NDI/RDI、DI/NDI、DI/RDI、RDI/NDI、DI/GDI和RI/GDI在不同分數階微分下的最大相關系數絕對值(Maximum absolute correlation coefficient,MACC)可以看出(圖3),優化指數DI/RDI與SOM的MACC顯著高于其他指數,平均值高達0.997 6,其次為NDI/RDI和RDI/NDI,均值分別為0.740 6和0.730 2,DI/NDI最低,MACC均值為0.309 3。同時,DI/RDI、NDI/RDI和RDI/NDI與SOM含量間的MACC在0.8階時達到峰值后逐漸下降,此時MACC分別為0.998 6、0.762 1和0.799 3。

圖3 不同分數階微分下光譜優化指數最大相關系數絕對值
為篩選最佳建模變量,本研究以MACC大于等于0.800 0為閾值進行優選。土壤DI/NDI、DI/GDI、RI/GDI、NDI/RDI和RDI/NDI在0~2階范圍內并無MACC大于等于0.800 0的敏感波段組合存在,故在后續建模中不作考慮。DI/RDI在0.2~2.0階范圍內的MACC為0.996 5~0.998 6,全部高于0.800 0,故在后續建模中采用DI/RDI指數。
圖4為SOM含量與優化光譜指數DI/RDI的二維相關系數分布圖,該圖能夠清晰展示相關系數較高特征指數的波段分布區間。與優化指數NDI/RDI和RDI/NDI的最大MACC相比(分布圖未展示),DI/RDI最大相關系數絕對值整體提高0.25,其敏感波段主要集中在1 450~1 750 nm和 2 100~2 400 nm之間(圖4),這也進一步說明分數階微分聯合優化指數在較低含量SOM的敏感波段篩選中,能夠有效消除土壤水分敏感波段在1 400 nm 和1 900 nm處的干擾[9]。

圖4 最優DI/RDI指數與SOM含量的二維相關系數


表4 基于DI/RDI-SVM的SOM建模與預測結果
在ArcGIS 10.2中,利用模型殘差結合普通克里格對實測和0.2階下DI/RDI-SVM模型的結果分別作反演插值(圖5),可以看出,銀川平原南部的吳忠市利通區、青銅峽市和巴音陶亥鎮黃河以西地區SOM含量明顯高于其他地區。賀蘭縣整體及惠農區東北部SOM含量較高,而銀川市西夏區西南部、永寧縣東西兩側及石嘴山市大武口區SOM含量整體較低,這是由于該區域地勢低洼、蒸發強烈,加之多年引黃灌溉,導致該區域土壤鹽漬化較重,土壤結構差,保水、保肥性能下降,故SOM含量低下。

圖5 銀川平原土壤有機質含量實測值和0.2階下DI/RDI-SVM反演值插值圖
與實測值相比,四級水平(SOM含量 10~20 g/kg)SOM的反演值所占面積比實測值高2.80%(表5),而五級(SOM含量6~10 g/kg)和六級(SOM含量小于6 g/kg)水平面積低了2.07%和0.73%,說明該方法反演的銀川平原SOM在相對較高水平反演值略高于實測值,而低水平反演值則略低于實測值。

表5 銀川平原不同等級有機質實測值與反演值所占面積

光譜指數能夠從二維光譜空間凸顯SOM的響應特征,充分利用高光譜信息,減少其他土壤信息對SOM的影響,降低模型的復雜性、去除冗余信息變量[29]。郭燕等[21]分析了NDI、DI、RI與SOM的二維相關性,最大相關系數在0.5~0.7之間;HONG等[22]基于不同光譜變換的光譜指數與SOM的二維相關性,得到最大相關系數為0.810 0,比一維光譜與SOM的相關性有了顯著提升。尚天浩等[23]采用單一光譜指數對銀川平原SOM估算,發現RDI的MACC最高,MACC可達0.801 0,RDI-SVM模型估測精度最高,RPD為2.32。目前大多應用RI、DI和NDI等單一光譜指數來反演土壤屬性,而兩兩組合很少使用[15]。本研究基于5種光譜指數的兩兩算法優化,發現NDI/RDI和RDI/NDI的MACC最大值分別為0.762 1和0.769 3,說明光譜指數NDI和RDI無論作何算法,其比值優化波段指數間的MACC范圍都差異不大;與此相比,DI/RDI的MACC顯著提升20%以上,最高達0.998 6,表明在土壤低有機質含量條件下,DI/RDI指數作算法優化時的作用優于單一NDI或RDI指數。這也說明SOC對優化光譜指數的靈敏度明顯高于使用單一光譜指數時的靈敏度,從而突出了使用光譜指數組合的優勢[7]。在SOM估測建模方法中,SVM作為一種有效監督機器學習的方法,通過對數據的智能化處理和數據價值的充分挖掘,在一定程度上解決了數據處理過程中的“過學習”和“離散值多”問題[23]。本研究利用分數階微分聯合不同優化指數所建SVM模型,預測RPD普遍高于2.00,其中DI/RDI-SVM模型RPD高達4.31,這是因為SVM為機器算法,其模型在數據驅動中具有極強的自學習能力,能夠準確獲取土壤光譜中極為復雜的非線性特征[33],SVM模型依靠核函數,經過多次訓練后,將輸入數據繪制到新的超空間中,在該超空間中執行分離,最終使用ε不敏感損失函數獲得用于數據擬合和預測的最優超空間[17],該超空間可以容忍小于常數ε集作為閾值的誤差[34],使得模型擬合精度不斷接近100%[35]。ZHANG等[20]還基于波段優化算法和兩波段指數形式構建了一個新的三波段指數——修正歸一化差分指數,并指出該三維指數在利用可見-近紅外光譜估計土壤其他生化參數方面具有很強的應用潛力。本研究后期也可以嘗試采用三維光譜指數建立SOM估算模型。此外,還需考慮銀川平原土壤屬性中的水分、鹽分、pH值和養分等因素對SOM光譜反演的影響。
(1)銀川平原土壤樣本SOM含量平均值為12.48 g/kg,屬于高度變異,其中93.05%的土樣SOM處于四級~六級水平。
(2)土壤野外原始光譜反射率吸收帶整體較寬,吸收特征差異明顯,在1 400 nm和1 900 nm處有明顯吸收峰。隨著分數階的不斷增加,光譜反射率不斷趨近于0;高分數階微分的吸收谷明顯小于低分數階微分。
(3)優化光譜指標MACC從大到小依次為DI/RDI、NDI/RDI、RDI/NDI、DI/GDI、DI/NDI和RI/GDI,敏感波段主要集中在1 450~1 750 nm和2 100~2 400 nm。
