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

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

0 引言

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

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

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

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

1 材料與方法

1.1 土壤樣品采集

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

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

1.2 土壤有機(jī)質(zhì)含量測定

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

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

1.3 土壤光譜測試

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

1.4 光譜預(yù)處理

剔除波段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種形式對平滑后的光譜曲線進(jìn)行處理。

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

1.6 模型建立與驗證

SOM反演模型構(gòu)建采用偏最小二乘回歸(Partial least square regression，PLSR)，其有利于解決自變量之間多重相關(guān)的現(xiàn)象，能夠避免光譜建模中的過擬合問題。選取決定系數(shù)(Determination coefficients，R2)和均方根誤差(Root mean squared error，RMSE)從模型的穩(wěn)定性和預(yù)測能力兩方面對模型的精度進(jì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 結(jié)果與討論

2.1 土壤光譜特征分析

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

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

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

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

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

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

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

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

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

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

CWT處理后的潮土光譜與SOM含量的相關(guān)性最強(qiáng)，決定系數(shù)最高可達(dá)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，對應(yīng)的尺度分別為第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，對應(yīng)的尺度分別為第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，對應(yīng)的尺度分別為第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，對應(yīng)的尺度分別為第7、7、7、7、8、8尺度。整體而言，不同類型土壤的SOM敏感波段基本上仍保持在相同的位置，并未發(fā)生大幅變動，且主要分布在近紅外范圍內(nèi)，可能是因為東北旱作農(nóng)田SOM主要源于農(nóng)作物殘體，其由糖類化合物、含氮化合物、纖維素等組成，這些成分中C—H鍵、C—O鍵、N—O鍵、N—H鍵等的光譜響應(yīng)波段位于近紅外區(qū)域[24]。

2.4 土壤有機(jī)質(zhì)敏感波段篩選

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

圖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 土壤有機(jī)質(zhì)含量PLSR反演模型Tab.2 PLSR prediction model of soil organic matter content

2.5 土壤有機(jī)質(zhì)高光譜反演模型

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

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

3 結(jié)論

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

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

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

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