基于高光譜數據的鹽漬化土壤含鹽量預測模型研究

關鍵詞：高光譜；土壤含鹽量；鹽漬化土壤；預測模型；多元逐步線性回歸

0 引言

土壤鹽漬化嚴重制約農業生產和發展，已成為全球普遍存在的環境問題，是一種常見的土地退化過程[1]。土壤鹽漬化是一個復雜的過程，分析鹽漬化土壤在時間與空間上的動態變化，能夠更好地掌握當地鹽漬化的程度與狀態，為農業生產、生態環境保護提供理論依據[2-3]。近年來研究人員研究了許多精確和高效的土壤鹽分定量分析方法，這些方法對土壤鹽漬化的改善和治理、土地資源的保護、農業生產效率的提高作用顯著，但對于大范圍區域很難做到及時、準確地監測[4-5]。因此，廣泛開展鹽化土壤動態監測、土壤含鹽量精確估算等研究備受關注[6-7]。

在遙感技術日新月異的今天，地物的細微之處在光譜數據中都能得到準確地反映[8]。王飛等[9]對新疆維吾爾自治區天山典型綠洲進行深入研究，利用線性與非線性模型對土壤鹽分敏感的原始波段、鹽分指數、植被指數等指標進行對比分析，發現SI-T、TGDVI、EEVI、BAND2和BAND5等指標對鹽分特別敏感。WHITNEYK等[10]通過對土壤含鹽量時空變化的遙感光譜指數和冠層溫度進行動態監測，對土壤大尺度鹽分時空序列進行數據分析。王麗娜等[11]利用主要成分回歸分析方法，基于對鹽漬土壤光譜曲線變化規律的研究進行分析，確定了黃河三角洲鹽分敏感波段1490～1608和1900～1950nm兩個波長范圍，并通過估算得到了區域土壤含鹽量。在土壤修復過程中，朱赟等[12]通過對相關系數極值和不同相關系數范圍分析，從8種光譜數據集中篩選最優波段，運用偏最小二乘法，建立了光譜反演模型。利用野外敏感波段，張同瑞等[13]建立了光譜指數，進一步構建18個模型并驗證選定，在黃河三角洲地區鹽分條件下，以土壤調節植被指數的線性模型為基礎，確定最優的反演分析模型。蔡東全[14]對土壤鹽化程度不同的光譜特征規律進行了分析，利用相關分析的方法獲得特征波段，以多元線性回歸分析方法建立高光譜預測模型，通過波譜角分類法得到了土壤鹽化信息，建立了土壤波譜庫。黃帥等[15]通過高光譜指數與土壤含鹽量的相關性分析，以最小二乘回歸和多元漸線性回歸的方法，篩選出土壤含鹽量的光譜特征波段，建立土壤鹽分動態監測模型。蒲智等[16]確定了最優光譜指數和定量預測模型，分析了土壤含鹽量與5種高光譜指數的定量關系。亞森江·喀哈爾[17]利用室內高光譜數據建立了渭干河?庫車河三角綠洲的PLSR和SVM土壤含鹽量反演模型，R2分別為0.762和0.637。為提升模型的反演精度，張曉光等[18]在原始光譜反射率的基礎上進行整數階微分變換，結果表明，一階和二階微分變換比原始光譜反射率對土壤含鹽量更加敏感，建立的模型精度更高。本研究以內蒙古自治區（簡稱內蒙古）科爾沁右翼中旗巴彥淖爾蘇木為研究區，根據室內分析結果，對鹽漬化土壤進行分析，依據高光譜數據，將光譜數據進行多種數學變換，并與土壤含鹽量進行相關分析，從而對響應波段進行優選提取。采用地統計分析方法構建鹽漬化土壤遙感監測模型，模擬測定含鹽數量與光譜數據之間的數學關系，為大規模定量監測鹽漬化土壤提供科學依據。

1 材料與方法

1.1 研究區概況

巴彥淖爾蘇木隸屬于內蒙古興安盟科爾沁右翼中旗，地處東經121°46′27.5″～122°12′54.8″、北緯44°30′22.3″～44°49′2.7″，海拔高度184m，面積527.1km2。巴彥淖爾蘇木屬溫帶大陸性季風氣候，春季少雨多風、夏季多雨高溫、秋季風大干燥。巴彥淖爾蘇木地處科爾沁腹地，地形屬沙坨平原地區，風沙土和草甸為主。巴彥淖爾蘇木耕地面積1萬hm2，由于生態破壞、土地退化、原生與次生鹽堿化并存，根據全國第四次沙化監測，巴彥淖爾蘇木鹽堿地面積1.429萬hm2，其中林地面積0.016萬hm2、草原面積1.413萬hm2，鹽堿地面積大、程度重，嚴重制約當地經濟、生態和農牧業的發展[19]。

1.2 數據采集與處理

1.2.1 土壤樣品采集與處理

2020年6月中旬在研究區內開展野外調查，按照梅花樁取樣法（重復采集5次，然后對土壤進行混勻）對研究區表層（0～10cm）土壤樣本進行采集，共采集土壤樣本100個。樣品在實驗室內自然風干，充分研磨后過2mm篩。每個樣品分2份，用于土壤光譜檢測1份，用于水溶性鹽總量指標測定1份。土壤水溶性鹽用水浴蒸干法測定總量：吸取一定量的土壤浸出液放在瓷蒸發皿上，在水浴上蒸干，用過氧化氫氧化有機質，在105～110°C烘箱中烘干，稱質量。

1.2.2 土壤樣本光譜測定

光譜檢測采用FieldSpec4光譜儀（analyticalspectraldevices），頻段范圍350～2500nm，采樣間隔1nm，輸出波段數2151個。在暗室內進行光譜測定，先校正后進行光譜測定。在暗室內對土壤樣本進行光譜檢測，檢測之前進行校正。將處理好的土樣裝滿盛樣皿，表面刮平后開始檢測，每個樣品檢測10次，以算術平均值作為實際光譜反射率。

1.3 數據處理

1.3.1 光譜數據平滑去噪

運用Origin8.6軟件中的Savitizky-Golay平滑法進行9點加權平均，對測得的光譜反射率進行平滑處理。光譜曲線中噪聲較大的波段范圍經過平滑處理后被刪除，400～2500nm波段的光譜信息得以保留。

1.3.2 實測光譜數據變換

由于原始光譜對噪聲的敏感性不強，原始光譜經過數學變換處理后，再與土壤含鹽量構建模型，對土壤含鹽量進行估算[20]。為了更好地反映土壤含鹽量與光譜反射率的相關關系，對原始光譜反射率R分別進行一階導數R′、二階導數R″、平方根、平方根一階導數（）′、平方根二階導數（）″、對數logR′、對數一階導數（logR）′、對數二階導數（logR）″、對數倒數[log（1/R）]、對數倒數一階導數[log（1/R）]′、對數倒數二階導數[log（1/R）]″、倒數（1/R）、倒數一階導數（1/R）′和倒數二階導數（1/R）″14種常見光譜變換，以選擇最佳的光譜處理方法。光譜微分處理計算公式為

1.3.3 相關性分析

對樣本光譜數據與土壤含鹽量的相關系數進行逐波段測算，對特征波段進行提取。計算相關系數公式為

1.4 建模與檢驗

1.4.1 模型建立方法

通過原始光譜反射率及其數學變換形式與土壤含鹽量進行相關性分析，提取特征波段。供試土壤樣本數為100個，以供試土壤樣本的70%用于模型建立，30%用于模型的檢驗，采用多元逐步線性回歸方法建立土壤含鹽量高光譜預測模型。

1.4.2 模型檢驗

利用回歸決定系數（R2）和均方根誤差ERMSE對建模集和驗證集進行檢驗。R2越大、ERMSE越小，表明模型穩定性和預測效果越好。計算公式為

2 結果與分析

2.1 土壤鹽漬化等級劃分

測得的土壤樣本指標值統計特征如表1所示。樣本土壤含鹽量最大值55.34g/kg、最小值0.52g/kg、平均值13.10g/kg及變異系數97.25%，屬強變異，表明研究區土壤含鹽量分布很不均勻。根據研究區實測點的情況，參考《中國鹽漬土》中土壤分級方法，確定研究區土壤鹽漬化分級標準，如表2所示[21]。

2.2 土壤光譜特征分析

將100個土壤樣本的反射率依照表2進行分類，在不同鹽漬化土壤樣本中，取其平均值作為光譜反射率曲線，如圖1所示。5條曲線的變化趨勢基本一致，供試土壤含鹽量越高，光譜的反射率也越大。光譜反射率在350～1300nm范圍內呈上升趨勢，增加幅度較大，在1400～1850nm范圍內反射率不斷升高，但增加幅度較緩，1900～2150nm范圍內光譜反射率增加且增加速度很快，從2150～2500nm范圍內，光譜反射率先逐漸減小，后又快速增加再急劇下降，在接近1400、1900和2200nm處出現3個明顯的水分吸收谷。

2.3 最佳波段選取

由于供試土壤含鹽量變化較大，對最佳波段的選取會造成一定影響，根據土壤鹽漬化分級標準將參試土壤分為兩類：非鹽土（土壤含鹽量≤7g/kg）和鹽土（土壤含鹽量>7g/kg），計算土壤含鹽量與15種數學變換的光譜數據的相關系數，如圖2所示。非鹽土光譜反射率與土壤含鹽量的相關系數最大0.3646，鹽土的原始光譜反射率與土壤含鹽量的相關系數最大0.4545，經過數學變換后，與土壤含鹽量的正負相關性明顯增強。非鹽土的對數一階微分處理反射率與土壤含鹽量在1490nm處相關性最高，相關系數最?0.5898，鹽土的對數一階微分反射率與土壤含鹽量在727nm處相關性最高，相關系數?0.5591。數據中的細微信息經數學變換后得以放大，可以據此篩選出提取pRpR土壤含鹽量的特征波段。根據圖2中相關分析，R′、R″、（）′、（）″、（logR）′、（logR）″、[log（1/R）]′、[log（1/R）]″、（1/R）′和（1/R）″10組變換下顯著相關的波段數較多，而其他變換沒有或極少存在與土壤含鹽量顯著相關的特征波段。因此從以上10種變換中選擇特征波段，如表3和表4所示。

2.4 土壤含鹽量估測模型構建

根據相關分析結果得到的特征光譜波段，建立土壤含鹽量的多元逐步線性回歸模型。運用DPS數據處理軟件將非鹽土47個樣本、鹽土53個樣本中的70%（非鹽土32個樣本，鹽土36個樣本）用于建模，30%（非鹽土15個樣本，鹽土17個樣本）用于模型驗證。以實測土壤含鹽量作為因變量，以特征波段的光譜反射率作為自變量，在保證R2盡可能大的前提下，建立土壤含鹽量的多元逐步線性回歸模型，如表3和表4所示。

利用多元逐步線性回歸分析方法建立的鹽漬化土壤含鹽量預測模型精度較高，非鹽土用于建模的R2平均值0.7306、最大值0.7685、最小值0.6514，ERMSE平均值0.97、最大值1.10、最小值0.90；鹽土用于建模的R2平均值0.5557、最大值0.7682、最小值0.3466，ERMSE平均值7.76、最大值9.50、最小值5.66。根據R2越大、ERMSE越小，模型越穩定的原則，非鹽土的原始二階方程建模的R2=0.7685、ERMSE=0.90，模型穩定性最好；鹽土的對數二階方程建模的R2=0.7682、ERMSE=5.66，模型穩定性最好。

2.5 模型精度檢驗

模型的驗證主要考慮穩定性和預判能力兩個方面。測試模型的預測能力，以未參加建模的土壤樣本中土壤含鹽量為準真值。由表3和表4可知，非鹽土檢驗模型檢測的R2最大值0.7292、最小值0.4059，ERMSE最大值1.91、最小值0.95；非鹽土的原始二階方程檢測的R2=0.7292、ERMSE=0.95，模型穩定性最好。鹽土檢驗模型檢測的R2最大值0.8718、最小值0.3944，ERMSE最大值9.77、最小值5.80，鹽土的對數二階方程檢測的R2=0.7682、ERMSE=5.66，模型穩定性最好。

由圖3可知，建模樣本基本上聚集在對角線附近，非鹽土原始光譜反射率二階微分預測模型實測值與預測值間的樣本擬合系數R2=0.8246，鹽土原始光譜反射率對數二階微分預測模型實測值與預測值之間的樣本擬合系數R2=0.7666，預測值和實測值具有較好的一致性。本研究建立的土壤鹽漬化遙感監測模型是有效的，經檢驗合格，可在一定程度上用于預測土壤含鹽量。

3 討論

在構建模型時發現，當對全部供試土壤反射率及其數學變換形式進行整體構建模型時，回歸決定系數（R2）整體不高，模型穩定性較差，預測效果不好，當將供試土壤按照土壤含鹽量進行分類，再對其反射率及其數學變換形式進行建模時，R2值整體升高，模型的穩定性和預測效果增強。這是由于供試土壤實測土壤含鹽量最大值55.34g/kg、最小值0.52g/kg、平均值13.10g/kg及變異系數97.25%，研究區域內土壤含鹽量分布極不均勻，導致在特征波段選取上存在差別。

本研究是對采自內蒙古自治區科爾沁右翼中旗巴彥淖爾蘇木典型鹽漬化土壤進行光譜反射率測定，并對其進行均方根、對數、對數倒數、倒數及其一階、二階微分等數學變換，利用多元逐步線性回歸分析方法，建立鹽漬化土壤含鹽量預測模型，并對預測模型進行對比。研究表明，非鹽土的原始光譜對數一階微分預測模型最優，鹽土的對數一階微分預測模型最優，這也說明選取不同的數學變換形式在模型預測上是不同的。因此本研究篩選的光譜變換形態所建立的預測模型，對于其他地區鹽土的適用性尚待考證。

4 結束語

（1）從形態上看，不同鹽漬化程度的土壤光譜特征曲線基本一致；土壤含鹽量越高的土壤樣本，其光譜反射率就越高。

（2）光譜反射率經數學變換后與土壤含鹽量的正負相關性增強，特別是經過一階和二階微分變換后，相關性顯著增強：非鹽土的原始光譜對數一階微分反射率與土壤含鹽量相關性最高，在1490nm處相關系數最大為?0.5898，在727nm處鹽土原始光譜對數一階微分反射率最大為?0.5591。

（3）利用多元逐步線性回歸建立的預測模型，非鹽土以原始光譜的二階微分模型為最優，檢測R2=0.7292，鹽土以對數二階微分模型為最優，檢測R2=0.8718，可快速測定鹽堿土土壤含鹽量。