基于CWT-sCARS的東北旱作農田土壤有機質高光譜反演

勾宇軒 趙云澤 李 勇 卓志清 曹 夢 黃元仿

(1.中國農業(yè)大學土地科學與技術學院， 北京 100193； 2.自然資源部農用地質量與監(jiān)控重點實驗室， 北京 100193；3.浙江大學農業(yè)遙感與信息技術應用研究所， 杭州 310058；4.中國建筑一局(集團)有限公司生態(tài)園林分公司， 北京 100071)

0 引言

土壤有機質(Soil organic matter, SOM)作為反映土壤肥力狀況和退化程度的重要指標，其含量的動態(tài)變化受到研究者的廣泛關注[1-2]。傳統(tǒng)SOM含量的測定主要通過濕燒法、干燒法、重鉻酸鉀容量法等化學分析法來實現(xiàn)，操作復雜、周期長、成本高，且為點狀信息數(shù)據(jù)，不適合快速、大面積測定，難以滿足土壤退化動態(tài)監(jiān)測和精準防治的要求。高光譜技術憑借快速、簡便、無污染、不破壞等優(yōu)勢，在土壤屬性預測研究中得到了廣泛應用。有機質含量不同的土壤在可見光和近紅外光譜區(qū)域存在獨特的光譜特征，為SOM含量的快速測定開辟了新的途徑[3-4]。近年來，國內外學者利用光譜技術進行SOM的預測取得了較好的效果[5-8]。

然而，由于受到土壤樣本表面、光譜測試環(huán)境、光譜高頻隨機噪聲、雜散光等干擾及光譜數(shù)據(jù)的信息冗余、多重共線性、吸收峰重疊等現(xiàn)象影響，模型的性能和精度易受到嚴重影響。前人研究多采用倒數(shù)、對數(shù)、微分、吸收峰深度等傳統(tǒng)變換來消除土壤光譜噪聲，通過無信息變量消除(Uninformative variables elimination，UVE)、連續(xù)投影算法(Successive projections algorithm，SPA)、競爭性自適應重加權采樣(Competitive adaptive reweighted sampling，CARS)等方法濾除無效變量或冗余變量[9-10]，優(yōu)選出SOM的敏感波段，以降低模型復雜度，提升模型的預測精度和穩(wěn)定性。如何對高光譜數(shù)據(jù)進行處理，消除干擾，濾除冗余、共線信息，篩選敏感波段，已成為當前研究的熱點之一。

連續(xù)小波變換(Continuous wavelet transform，CWT)通過多尺度分解，得到光譜信號中的近似特征和細節(jié)信息，在高光譜數(shù)據(jù)的噪聲去除和數(shù)據(jù)壓縮上具有獨特優(yōu)勢[11-13]。CARS算法在敏感波段的篩選上具有很好的性能，但只考慮了變量回歸系數(shù)，并不能總是反映變量重要性的真實信息。ZHENG等[14]構建了穩(wěn)定性競爭自適應重加權采樣(Stability competitive adaptive reweighted sampling, sCARS)算法克服了這一缺點。劉國海等[15]提出了基于sCARS策略的特征波長篩選方法，與UVE、CARS等現(xiàn)有方法相比，能夠有效地增強光譜特征波段選擇的準確性和穩(wěn)定性。何東健等[16]比較了sCARS、SPA、CARS和RF(隨機蛙跳算法)4種波長選擇方法對建模效果的影響，發(fā)現(xiàn)sCARS的結果最佳。李冠穩(wěn)等[17]使用sCARS算法結合隨機森林模型，實現(xiàn)了SOM含量的快速、無損、精準估測。

東北地區(qū)土層深厚，土壤性狀良好，是我國重要的商品糧基地，但是由于長期的高強度利用，加之侵蝕作用的影響，其SOM含量顯著下降，建立可靠的SOM含量動態(tài)監(jiān)測方法，對東北土壤退化防治和耕地質量提升有重要的現(xiàn)實意義。本文探究東北旱作農田典型土壤類型的光譜反射特性，采用CWT對高光譜數(shù)據(jù)進行分析，利用sCARS算法篩選敏感波段，提取有益信息，提高SOM預測模型的精度和穩(wěn)定性，為建立可靠的土壤退化監(jiān)測網(wǎng)絡提供有力的支撐，進一步推動高光譜技術在SOM含量預測中的應用。

1 材料與方法

1.1 土壤樣品采集

在黑龍江、吉林、遼寧3省地形坡度小于5°、1 km2網(wǎng)格旱地占耕地面積40%以上的區(qū)域內，以15 km×15 km網(wǎng)格布點，結合土壤類型進行分層抽樣，綜合考慮種植體系、分布面積和集中程度等因素，共確定118個采樣點(圖1)，其中黑土采樣點37個，黑鈣土采樣點33個，潮土采樣點28個，棕壤采樣點20個。土壤樣品采集時間為2017年5—6月，采樣深度為0～20 cm。采集的土壤樣品自然風干、研磨、過2 mm篩后分為兩份，一份用于土壤有機質含量測定，一份用于土壤光譜測試。

圖1 研究區(qū)及采樣點地理位置Fig.1 Study area and sampling point

1.2 土壤有機質含量測定

采用重鉻酸鉀-外加熱法測定樣品的SOM含量，統(tǒng)計特征見表1。東北旱作農田SOM含量均值為26.70 g/kg，整體處于較高水平，數(shù)據(jù)離散程度一般，屬于中等變異強度(15%～100%)。黑土SOM含量最高，且變異系數(shù)最小，離散程度最低，潮土和棕壤SOM含量最低。

表1 土壤有機質含量的統(tǒng)計特征Tab.1 Statistical characteristics of SOM content

1.3 土壤光譜測試

在暗室中使用ASD FieldSpec 4型便攜式地物光譜儀進行土壤光譜數(shù)據(jù)測試。將土壤樣本放置于深1.5 cm、半徑5 cm的玻璃盛樣皿內，選擇50 W鹵素燈為光源，保持15°頂角，距離土壤樣品表面30 cm，采用8°視場角，使探頭垂直位于土壤樣本正上方15 cm處，使用前先進行白板校正，儀器穩(wěn)定后再進行實驗。每個土樣采集10條光譜曲線，算術平均后得到該土樣的光譜反射率。儀器的光譜波長為350～2 500 nm，重采樣間隔為1 nm，共輸出2 151個波段。

1.4 光譜預處理

剔除波段350～399 nm和2 401～2 500 nm，以消除邊緣噪聲的干擾。通過Savitzky-Golay算法(多項式階數(shù)為2，點數(shù)為11)平滑400～2 400 nm的土壤光譜曲線，并采用倒數(shù)對數(shù)、一階微分、連續(xù)統(tǒng)去除(Continuum removal，CR)和CWT 4種形式對平滑后的光譜曲線進行處理。

CWT源于傅里葉算法，可同時在頻率和時間開展分析，從信號中提取有效信息，其通過小波基函數(shù)將一維的土壤高光譜數(shù)據(jù)分解成包含波段和不同分解尺度的二維小波系數(shù)，計算式為

(1)

(2)

式中Wf(a,b)——小波系數(shù)，包含i、j兩維，分別是分解尺度(i=1,2,…,m)和波段(j=1,2,…,n)組成的m×n矩陣

f(t)——土壤高光譜反射率

t——光譜波段(400～2 400 nm)

ψa,b——小波基函數(shù)

a——尺度因子

b——平移因子

1.5 敏感波段篩選

sCARS算法首先計算各波長變量的穩(wěn)定性值，然后利用自適應重加權采樣技術(Adaptive reweighted sampling，APS)和指數(shù)衰減函數(shù)(Exponentially decreasing function，EDF)篩選出回歸系數(shù)絕對值大且穩(wěn)定性高的變量，多次循環(huán)迭代此過程。最終通過十折交互檢驗對每次循環(huán)后所得的變量子集進行檢驗，優(yōu)選出交互驗證均方根誤差(RMSECV)最小的變量子集[18]。

1.6 模型建立與驗證

SOM反演模型構建采用偏最小二乘回歸(Partial least square regression，PLSR)，其有利于解決自變量之間多重相關的現(xiàn)象，能夠避免光譜建模中的過擬合問題。選取決定系數(shù)(Determination coefficients，R2)和均方根誤差(Root mean squared error，RMSE)從模型的穩(wěn)定性和預測能力兩方面對模型的精度進行驗證。

光譜數(shù)據(jù)的平滑、PLSR建模分析在Unscrambler X(CAMO Inc., 挪威)軟件中實現(xiàn)，CR處理通過ENVI 5.3(ESRI Inc., 美國)軟件實現(xiàn)，CWT處理和sCARS波段篩選在Matlab R2018a(The Mathworks Inc., 美國)環(huán)境下完成。

2 結果與討論

2.1 土壤光譜特征分析

按照SOM含量將所有土壤樣品分為6個等級，計算各等級所有土壤光譜反射率平均值，得到不同有機質含量下土壤的光譜反射率曲線(圖2)。研究區(qū)土壤光譜反射率隨著SOM含量的增加而減小，不同有機質含量等級的土壤光譜反射率曲線整體變化趨勢相近。整體來看，光譜曲線表現(xiàn)平緩，光譜反射率隨波長而增大，在400～900 nm波段區(qū)域反射率增加較快，斜率較大，而在近紅外的900～2 400 nm波段反射率增加較為平緩。

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

由不同土壤類型的光譜反射率曲線變化特征可知(圖3)，黑土和黑鈣土的光譜反射率較低，這是由于黑土、黑鈣土SOM含量均值較高，顏色較深，土壤質地粘細，其光譜曲線屬于富含有機質類型[19]。一般情況下，隨著SOM含量的增加，光譜反射率減小。潮土和棕壤雖然SOM含量均值相近，但與棕壤相比，潮土的光譜反射率相對較低，可能受到氧化鐵含量等其他因素的影響[20]。整體來看，不同土壤類型的光譜特征吸收帶出現(xiàn)的位置總體上一致，只是深度有所不同。在可見光波段，主要受土壤發(fā)色團和有機質本身顏色的影響，存在較寬的吸收波段；在近紅外波段，則主要受到N—H、C—H和C—O等基團的影響。在波長1 400、1 900、2 200 nm處土壤光譜曲線有3個明顯的水分吸收谷，這可能是黏土礦物中所含的水分子和羥基的吸收帶[21]。

圖3 不同類型土壤的光譜反射率曲線Fig.3 Reflectance spectra of different soils

2.2 傳統(tǒng)光譜變換的相關性分析

3種傳統(tǒng)光譜變換對相關性的提升效果并不理想(圖4)。原始光譜與SOM含量的決定系數(shù)最大為0.41，倒數(shù)對數(shù)變換后，土壤光譜與SOM含量的相關性無明顯變化；經(jīng)過連續(xù)統(tǒng)去除處理，相關性有所降低，決定系數(shù)最大僅為0.36；一階微分處理后相關性有一定的提升，決定系數(shù)最大可達0.46，提高了0.05，這是由于微分技術對分離重疊樣品、消除平行噪聲的干擾有一定的效果。原始光譜和倒數(shù)對數(shù)R2>0.4的波段主要分布在670～890 nm，并在794 nm處達到最大值；一階微分R2>0.4的波段主要分布在530～630 nm和1 920～1 950 nm，并在594 nm和1 930 nm處達到最大值。

圖4 不同光譜變換形式的相關性分析Fig.4 Correlation analysis of different spectral transformation forms

2.3 連續(xù)小波變換的相關性分析

由于土壤光譜曲線特征與Gaussian函數(shù)近似，選取Gaussian4函數(shù)作為小波基函數(shù)[22]，CWT分解尺度設定為21、22、…、210，即分解尺度1～10。將變換后的小波系數(shù)與SOM含量進行相關性分析，得到不同類型土壤小波系數(shù)與SOM含量的決定系數(shù)矩陣(圖5)，圖中紅色區(qū)域代表相關性強的區(qū)城，藍色區(qū)域代表相關性弱的區(qū)城。

由圖5可以看出，有效光譜信息主要集中于中高尺度分解。在分解尺度為1～3時，信息相對均一，部分光譜細節(jié)特征消失，有效信息比較匱乏；隨著分解尺度增大，光譜對SOM的敏感性也逐步增大，CWT處理后土壤光譜與SOM含量的相關性明顯高于原始光譜，這表明土壤的部分有效光譜信息被隱藏，利用連續(xù)小波分析，適度增大分解尺度可有效挖掘其內的隱含信息[23]。

CWT處理后的潮土光譜與SOM含量的相關性最強，決定系數(shù)最高可達0.69，而黑土較差。將決定系數(shù)最高的5%區(qū)域作為SOM的敏感區(qū)域，黑土的敏感區(qū)域主要分布在563～595 nm、563～612 nm、788～987 nm、2 084～2 126 nm、425～468 nm、755～902 nm、857～947 nm、1 282～1 630 nm，對應的尺度分別為第5、6、6、7、8、8、9、9尺度；黑鈣土的敏感區(qū)域主要分布在1 384～1 425 nm、2 319～2 359 nm、479～532 nm、872～912 nm、1 416～1 469 nm、2 166～2 207 nm、2 257～2 301 nm、1 498～1 703 nm、1 125～1 265 nm，對應的尺度分別為第4、5、6、6、7、7、7、9、10尺度；潮土的敏感區(qū)域主要分布在2 117～2 173 nm、2 252～2 282 nm、1 548～1 609 nm、431～542 nm、1 035～1 167 nm、1 823～2 090 nm、400～537 nm，對應的尺度分別為第7、7、8、9、9、9、10尺度；棕壤的敏感區(qū)域主要分布在419～525 nm、607～737 nm、1 621～1 694 nm、1 751～2 125 nm、417～507 nm、717～829 nm，對應的尺度分別為第7、7、7、7、8、8尺度。整體而言，不同類型土壤的SOM敏感波段基本上仍保持在相同的位置，并未發(fā)生大幅變動，且主要分布在近紅外范圍內，可能是因為東北旱作農田SOM主要源于農作物殘體，其由糖類化合物、含氮化合物、纖維素等組成，這些成分中C—H鍵、C—O鍵、N—O鍵、N—H鍵等的光譜響應波段位于近紅外區(qū)域[24]。

2.4 土壤有機質敏感波段篩選

以黑鈣土原始光譜為例，采用sCARS算法篩選SOM敏感波段的過程見圖6，可以看出，隨著迭代次數(shù)的增加，所保留的敏感波段數(shù)量逐漸減少，RMSECV呈現(xiàn)先減小再增大的變化趨勢，當運行次數(shù)為73次時，RMSECV最小(4.86 g/kg)。這是由于在1～73次變量篩選運行過程中，不斷剔除與SOM含量相關性較弱，對建模結果影響不大的波段，使RMSECV減小；而在73次之后，開始剔除與SOM含量相關性強的波段，導致RMSECV增大。當運行次數(shù)為73次時，得到的敏感波段子集最佳，共篩選了13個敏感波段。圖7為篩選后13個黑鈣土敏感波段的分布情況，其主要分布在小于1 200 nm和大于2 200 nm的波段。不同處理篩選后敏感波段的數(shù)量(表2)僅占初始波段的0.07%～2.85%，有效去除了冗余的信息變量，提取了與SOM相關的重要特征信息變量，減少了參與建模的變量數(shù)量，降低了模型的復雜度。

圖6 sCARS敏感波段篩選過程Fig.6 Screening process of sensitive bands by sCARS method

圖7 sCARS篩選后黑鈣土敏感波段分布Fig.7 Distribution of sensitive band of chernozem by sCARS method

表2 土壤有機質含量PLSR反演模型Tab.2 PLSR prediction model of soil organic matter content

2.5 土壤有機質高光譜反演模型

將篩選后的敏感波段作為反演模型的輸入變量，采用PLSR方法構建SOM含量反演模型，進行內部交叉驗證(表2)?？傮w來看，各類型土壤的有機質反演效果均較理想，僅黑土原始光譜(R)、倒數(shù)對數(shù)(lg(1/R))和連續(xù)統(tǒng)去除(CR)模型驗證R2低于0.50。黑土、黑鈣土、潮土和棕壤最大R2分別達到0.83、0.88、0.93和0.93，潮土和棕壤的模型精度及穩(wěn)定性略優(yōu)于黑土和黑鈣土。

比較分析4種處理的模型效果，一階微分(R′)和連續(xù)小波變換(CWT)模型的效果明顯優(yōu)于倒數(shù)對數(shù)和連續(xù)統(tǒng)去除，與相關性分析結果一致。不同類型土壤建模集最佳模型均為CWT模型，驗證集黑土、潮土的最佳模型為CWT模型，黑鈣土、棕壤一階微分模型驗證效果最好。與原始光譜相比，倒數(shù)對數(shù)處理后，僅潮土的R2提升了0.04，RMSE下降了0.46 g/kg，其他類型土壤的模型精度和穩(wěn)定性都略有降低；連續(xù)統(tǒng)去除處理后，不同類型土壤的R2提升了0.01～0.06，模型效果得到了一定程度提高，但效果不明顯。而一階微分和連續(xù)小波變換處理后，不同土壤類型的模型預測能力和穩(wěn)定性均有較大提高，一階微分處理后，建模集、驗證集的R2最大提升了0.11和0.16，RMSE最大降低了1.33、1.44 g/kg；連續(xù)小波變換處理后，建模集、驗證集的R2最大提升了0.13、0.28，RMSE最大降低了2.48、2.40 g/kg。這是因為連續(xù)小波技術不僅能有效地挖掘土壤光譜內隱含的有效信息，還能夠壓制噪聲干擾的負面影響，提升光譜信息的信噪比，從而提升對SOM的估測精度與穩(wěn)定性，而一階微分處理雖然有效消除了基線和低頻噪聲的干擾，但同時也引入了高頻噪聲。

3 結論

(1)東北旱作農田黑土和黑鈣土有機質含量豐富，光譜反射率較低；潮土和棕壤雖然有機質含量均值相近，但受氧化鐵含量等其他因素的影響，潮土的光譜反射率低于棕壤。

(2)一階微分對于分離重疊樣品、消除平行噪聲干擾有一定效果；CWT不僅可以抑制背景和噪聲的干擾，還能夠挖掘土壤光譜內隱含的有效信息，提升光譜信息的信噪比，可以有效提升土壤光譜與有機質含量的相關性。

(3)sCARS算法能夠提取與SOM相關的重要特征信息變量，去除冗余、重疊的光譜信息，提高建模效率。

(4)CWT構建的SOM含量反演模型預測精度和穩(wěn)定性最佳，黑土、黑鈣土、潮土和棕壤R2分別達到0.83、0.88、0.93和0.93，CWT技術結合sCARS算法，能夠較好地反演東北旱作農田典型土壤類型的SOM含量，為SOM含量的高光譜預測提供了新途徑。