譚 洋，姜琦剛，劉驊欣，劉 斌，高 鑫，張 博

吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院，吉林 長(zhǎng)春 130026

引 言

土壤是由礦物質(zhì)、有機(jī)質(zhì)(soil organic matter，SOM)、水分(soil moisture content，SMC)和空氣等多種成分構(gòu)成的復(fù)雜系統(tǒng)。土壤成分檢測(cè)和分析可用于提高土壤條件評(píng)估，指導(dǎo)精細(xì)化生產(chǎn)的農(nóng)業(yè)活動(dòng)。其中，SOM和SMC是土壤及農(nóng)作物生長(zhǎng)的重要因素，關(guān)系土壤肥力，影響農(nóng)作物長(zhǎng)勢(shì)；總鐵(Fe)是土壤礦物中的主要元素之一，包括Fe2+，F(xiàn)e3+，復(fù)合其他化學(xué)成分的鐵，如EDTA，以及各種鐵的氧化物，是植物生命活動(dòng)必須的微量元素；土壤pH值描述了堿基飽和水平，與微量元素的有效性、以及與潛在有毒元素，如鋁，有密切關(guān)系。傳統(tǒng)基于樣本的化學(xué)分析方法成本高、檢測(cè)周期長(zhǎng)，反射率光譜法已成為快速、可靠地獲取土壤組成信息的一種替代方法?；诓煌ㄩL(zhǎng)反射率與土壤成分的響應(yīng)關(guān)系，可見光-近紅外反射光譜已被廣泛應(yīng)用于預(yù)測(cè)不同土地利用類型土壤的SOM、SMC[1-2]等。

在光譜數(shù)據(jù)處理中，選擇一種適合的信號(hào)處理方法至關(guān)重要，有效的光譜處理方法可以極大地提高估測(cè)模型的精度。對(duì)于不同土壤成分的提取，通?？刹捎貌煌臄?shù)據(jù)處理方法獲得理想的精度[3-4]，如利用多元散射校正(multiple scatter correction，MSC)方法提高總鐵、游離鐵和無定形鐵的估算精度，平滑濾波(savitzky golay，SG)結(jié)合MSC提高土壤pH和其他土壤成分的估算精度[5]，分?jǐn)?shù)階導(dǎo)數(shù)處理提高鹽漬化土壤SMC估測(cè)精度[6]，利用離散小波變換提高土壤重金屬元素的估測(cè)精度[7]。大部分研究為提高土壤成分的預(yù)測(cè)精度選擇烘干樣本測(cè)量反射率，或利用外部參數(shù)正交化法[8]等去除水分影響。且常規(guī)土壤光譜數(shù)據(jù)處理方法大多著重于突出單一估測(cè)目標(biāo)的變化特征，削弱其他目標(biāo)的影響[9]。為同時(shí)提高多種土壤成分的預(yù)測(cè)精度，需要減弱顆粒大小不均和粒子表面帶來的散射影響對(duì)光譜反射率測(cè)量帶來的誤差，并增強(qiáng)不同成分響應(yīng)特征的細(xì)節(jié)信息。

標(biāo)準(zhǔn)正態(tài)變量變換(standard normal variable，SNV)針對(duì)每一樣本的光譜數(shù)據(jù)，降低樣本表面的非特異性散射噪聲。多尺度連續(xù)小波變換(continuous wavelet transform，CWT)在不同程度上放大光譜數(shù)據(jù)的細(xì)節(jié)特征[10]，可有效探測(cè)化學(xué)元素，樹冠水分[11]，葉片的氯離子濃度[12]等。本文提出使用SNV-CWT多尺度分解法處理濕土樣本高光譜反射率數(shù)據(jù)，以隨機(jī)森林(RF)預(yù)測(cè)模型的精度為標(biāo)準(zhǔn)，對(duì)比7種常用處理方法，確定不同成分的最佳分解尺度。利用三種方法計(jì)算和分析光譜分解前后與各成分相關(guān)關(guān)系，并實(shí)現(xiàn)各成分特征波段的提取與建模分析。

1 實(shí)驗(yàn)部分

1.1 研究區(qū)概況

東北黑土區(qū)主要分布在吉林省和黑龍江省，少量分布在遼寧省北部和內(nèi)蒙古東部，如圖1所示。根據(jù)以往研究，黑土區(qū)土壤肥沃，SOM含量豐富，且自南向北遞增，極少部分SOM含量少于2%。采樣點(diǎn)位于松遼平原北部的黑土帶內(nèi)，主要分布在黑龍江省中部，海倫市(44°04′—44°45′N，125°26′—126°31′E)和吉林省中部，德惠市(47°00′—47°41′N，126°25′—127°06′E)。

圖1 研究區(qū)及采樣點(diǎn)分布

1.2 土壤樣本采集、成分測(cè)定及光譜測(cè)量

土壤樣本采集于2017年5月4日至9日，正值春耕時(shí)期，土壤表面冰雪完全消融，且無植被覆蓋。共有74個(gè)采樣點(diǎn)，每個(gè)采樣點(diǎn)在5 m×5 m范圍采集5份土壤，采集深度為0～15 cm，所有樣本均被獨(dú)立裝入自封袋并標(biāo)記樣本號(hào)，采樣點(diǎn)中心地理坐標(biāo)由便攜式GPS記錄。

在實(shí)驗(yàn)室條件下將每個(gè)采樣點(diǎn)的5份土壤混合均勻，過2 mm篩剔除碎石子殘留秸稈等雜質(zhì)。每份樣本共取出100 g土壤分別用于測(cè)定各土壤成分含量，其中，利用重力法測(cè)定SMC含量，重鉻酸鉀加熱法測(cè)定SOM含量，點(diǎn)位法測(cè)定土壤pH值，等離子光譜儀測(cè)定Fe含量。土壤樣本成分含量統(tǒng)計(jì)如表1所示。其中SOM含量介于2.0%～6.3%，平均值約為3.5%，整體含量較高；土壤樣本SMC近似于野外土壤水分含量，值域分布區(qū)間較大，介于3.8%～24.4%；Fe含量介于2.3%～3.4%之間，且標(biāo)準(zhǔn)偏差(SD)與變異系數(shù)(CV)均較小(SD=2.12，CV=7.81)，總體含量差異不明顯；pH值介于4.7～8.1之間，均值為5.7，標(biāo)準(zhǔn)偏差較小(SD=0.77)，總體偏離平均值較小。

表1 研究區(qū)土壤樣本成分值統(tǒng)計(jì)

室內(nèi)光譜測(cè)量在黑暗的實(shí)驗(yàn)室進(jìn)行，將過篩土壤樣本置于深度1 cm的鋁制測(cè)試盒，刮平表面，用90 W的鹵素?zé)粽障虮砻?。采用ASD FieldSpec3光譜儀測(cè)定光譜反射率，每個(gè)樣本在不同位置共計(jì)測(cè)量20條光譜曲線。利用ViewSpecPro軟件計(jì)算平均光譜作為該樣本的原始反射光譜曲線。

1.3 光譜數(shù)據(jù)處理

光譜數(shù)據(jù)首先以原始光譜反射率為輸入矩陣，進(jìn)行標(biāo)準(zhǔn)正態(tài)變量變換。逐一將光譜反射率數(shù)據(jù)X=(X1，X2，…，Xp)在p波長(zhǎng)處轉(zhuǎn)換為Z=(Z1，Z2，…，Zp)，其中

Zi=(Xi-m)/s

(1)

式(1)中，m和s分別表示X1，X2，…，Xp的均值和標(biāo)準(zhǔn)差[13]。以變換后的數(shù)據(jù)Zn×p(n為樣本數(shù))作為輸入矩陣進(jìn)行多尺度連續(xù)小波變換。光譜反射率的連續(xù)小波變換可以表示為小波基函數(shù)Ψ(λ)經(jīng)過比例縮放和平移所產(chǎn)生的一系列系數(shù)，用式(2)表示為[12]

(2)

式(2)中，Wf(a,b)為小波系數(shù)，a和b均為實(shí)數(shù)，a為比例因子(小波基的寬度)，b為平移因子(小波基的位置)。小波系數(shù)Wf(a，b)為小波基與光譜反射率在不同尺度下的關(guān)系。

參與對(duì)比的數(shù)據(jù)處理方法包括一階導(dǎo)數(shù)(first derivative, FD)，二階導(dǎo)數(shù)(second derivative，SD)，連續(xù)同去除(continuum removal，CR)，標(biāo)準(zhǔn)正態(tài)變量變換(SNV)，多元散射矯正(MSC)，高斯濾波法(Gaussian filter，GS)，以及對(duì)數(shù)變換[logarithm transform，log(1/R)]。其中FD和SD變換，首先對(duì)離散光譜反射率進(jìn)行具有11個(gè)平滑點(diǎn)的二階多項(xiàng)式擬合，CR通過調(diào)用Tripack R包中的Delaunay三角測(cè)量的凸包實(shí)現(xiàn)光譜擬合；其余處理方法均利用Python實(shí)現(xiàn)，且均以濕土壤樣本的原始光譜反射率為輸入數(shù)據(jù)。為減少反射率邊緣的噪聲，原始光譜反射率僅保留400～2 500 nm范圍的值。為了可視化不同處理方法的差異，所有的信號(hào)變換都繪制在圖2中，其中SNV-CWT處理方法繪制的是10尺度(210)的分解系數(shù)。

圖2 經(jīng)過不同處理方法的光譜曲線

1.4 特征分析、建模與精度驗(yàn)證

皮爾森相關(guān)系數(shù)(Pearson correlation coefficients，PCC)常用來衡量特征與目標(biāo)變量間的線性相關(guān)程度，定義為兩個(gè)變量之間的協(xié)方差與標(biāo)準(zhǔn)差的商?；谀Ｐ偷南禂?shù)(model based coefficients，MBC)是一種基于模型的相關(guān)關(guān)系指標(biāo)，與傳統(tǒng)的線性方法相比，可以在建模之前判斷特征變量與目標(biāo)變量間基于模型的非線性關(guān)系?；疑P(guān)聯(lián)分析(grey relation analysis，GRA)是一種多因素統(tǒng)計(jì)分析方法，通過計(jì)算參考數(shù)列和比較數(shù)列變化曲線幾何形狀相似程度判斷因素間的關(guān)聯(lián)程度，由此方法計(jì)算得出的灰色關(guān)聯(lián)度(GRD)可作為判定特征數(shù)據(jù)與相關(guān)因素序列的相關(guān)性指標(biāo)。研究中采用三種方法對(duì)SNV-CWT處理前后的光譜數(shù)據(jù)進(jìn)行特征分析，判斷不同土壤屬性的響應(yīng)波段范圍。

隨機(jī)森林(random forest，RF)回歸是一種基于決策樹的機(jī)器學(xué)習(xí)算法，采用bootstrap抽樣技術(shù)抽取訓(xùn)練集，并建立每個(gè)訓(xùn)練集的CART樹，構(gòu)成森林，具有強(qiáng)大的預(yù)測(cè)能力，已被廣泛應(yīng)用于不同科學(xué)領(lǐng)域。該方法學(xué)習(xí)速度快，可以處理大量的輸入變量，即使在特征變量數(shù)遠(yuǎn)大于樣本數(shù)據(jù)時(shí)也可以很好地避免過擬合，且模型抗粗差能力較強(qiáng)。

采用決定系數(shù)(coefficient of determination，R2)、均方差(mean squared error，Mse)和相對(duì)分析誤差，即標(biāo)準(zhǔn)差與預(yù)測(cè)均方根誤差比(residual predictive deviation，RPD)這3個(gè)指標(biāo)檢測(cè)與評(píng)價(jià)模型的穩(wěn)定性和預(yù)測(cè)能力。其中R2越接近1說明模型擬合程度越好；當(dāng)RPD<1.4時(shí)，模型質(zhì)量較差，無法對(duì)樣本進(jìn)行預(yù)測(cè)；當(dāng)1.42.5時(shí)，模型具有極好的預(yù)測(cè)能力[14]。檢驗(yàn)與校正模型的統(tǒng)計(jì)指標(biāo)(statistical indicators，SI)計(jì)算公式如式(3)—式(5)

(3)

(4)

(5)

2 結(jié)果與討論

2.1 原始數(shù)據(jù)分析

原始光譜反射率曲線在1 400，1 900及2 200 nm附近均具有明顯的拐點(diǎn)。濕土壤條件下的光譜反射率較干土相比，在1 400和1 900 nm附近具有較大程度的波動(dòng)。如圖3所示，兩條曲線的有機(jī)質(zhì)含量相近，分別為6.0%和6.1%，但含水量差異較大，分別為15.8%和5.4%。干土情況下兩樣本的光譜反射率曲線形狀相似，含水量不同時(shí)，濕土樣本兩條光譜反射率形狀具有較大差異，說明水分對(duì)反射率具有顯著影響。而Fe、pH與反射率的最大PCC分別為0.18與0.25，與SOM的最大值(0.63)相比對(duì)反射率的影響較低。

圖3 干濕土壤樣本原始光譜反射率示意圖

2.2 數(shù)據(jù)處理方法的有效性分析

表2 不同數(shù)據(jù)處理方法下的SOM，SMC，F(xiàn)e和pH驗(yàn)證模型精度

2.3 多尺度SNV-CWT的有效性及各成分最佳分解尺度

研究中以Mexh為小波基函數(shù)，進(jìn)行10個(gè)尺度(21，22，23，…，210)的SNV-CWT分解，將一維光譜反射率分解為10個(gè)尺度的小波系數(shù)。以原始光譜數(shù)據(jù)集74×2 101維矩陣為輸入，經(jīng)SNV-CWT處理后轉(zhuǎn)換為74×10×2 101維矩陣。針對(duì)四種土壤成分，選擇每個(gè)尺度置信度小于0.05的小波系數(shù)作為輸入變量，建立RF回歸模型，以驗(yàn)證模型獲得的決策系數(shù)R2為標(biāo)準(zhǔn)，確定各土壤成分最佳分解尺度。如圖4所示，SMC的模型精度在不同尺度上整體優(yōu)于SOM，且SOM和SMC的建模精度較Fe和pH更好。這是由于SMC和SOM對(duì)光譜反射率影響起主導(dǎo)作用，對(duì)Fe與pH驗(yàn)證模型信息具有掩蓋作用，且Fe的含量變化區(qū)間小。SOM與SMC的最佳分解尺度分別為7和8；Fe與pH驗(yàn)證模型的最佳分解尺度為1和10。每種土壤成分對(duì)應(yīng)于不同的分解尺度有明顯的增強(qiáng)信息，相比于單一尺度分解，多尺度SNV-CWT分解對(duì)多屬性探測(cè)具有更大優(yōu)勢(shì)。

圖4 SNV-CWT多尺度情況下的不同目標(biāo)變量的模型決策系數(shù)

2.4 特征分析與建模

分別計(jì)算SOM，SMC，F(xiàn)e和pH與原始反射率以及在最佳分解尺度系數(shù)的PCC，MBC和GRD，其中PCC為線性相關(guān)系數(shù)的平方，MBC為單變量RF模型決策系數(shù)，GRD計(jì)算時(shí)取分辨系數(shù)ρ=0.5。從圖5中可以看出，經(jīng)過SNV-CWT變換后，三種方法計(jì)算得出的相關(guān)系數(shù)均較原始反射率有所提高，且SOM與SMC相關(guān)系數(shù)均高于Fe與pH模型。這也充分解釋了多維特征建模時(shí)Fe與pH的模型精度在不同尺度上均不及SOM與SMC。與土壤成分對(duì)應(yīng)的相關(guān)波段通過分析可分為兩類：一類是三種方法計(jì)算的相關(guān)系數(shù)曲線在相同的波長(zhǎng)范圍均有較明顯的峰值，這種波段可以認(rèn)為是該成分的主要特征波段；另一類是一種或兩種分析方法計(jì)算的相關(guān)系數(shù)較高而對(duì)應(yīng)的其他方法相關(guān)系數(shù)較低，這種波段屬于該成分的特征波段。SOM在450，580和820 nm附近的波段屬于主要特征波段范圍；超過1 400 nm，例如1 495，2 210以及2 420 nm左右，屬于SOM的特征波段。而SMC的特征波段主要分布在大于1 400 nm，從PCC和MBC曲線可以看出波段均表現(xiàn)出較高的相關(guān)性。1 460和1 950 nm附近是明顯的主要特征波段，與前人總結(jié)的SMC特征波段一致。Fe在整個(gè)可見光近紅外光譜范圍區(qū)間，并沒有很明顯的主要特征峰，但較原始反射率相比，三種分析方法計(jì)算得出的所有波段相關(guān)度有所提高，PCC的最大值從0.18提升至0.47。從GRD曲線看，特征間的相關(guān)系數(shù)大小差異不明顯；從PCC與MBC曲線判斷，730，1 360，1 750以及2 110 nm附近屬于Fe的特征波段。pH的特征波段分散在整個(gè)范圍內(nèi)，從PCC與GRD曲線來看，特征波段集中在460，920，2 240和2 430 nm附近；從MBC來看，還集中在1 100和1 740 nm附近。

圖5 SNV-CWT處理前后的反射率與不同成分的關(guān)系

總體分析，對(duì)反射率具有主導(dǎo)影響的屬性，如SOM與SMC，不同分析方法均可篩選出其對(duì)應(yīng)的主要特征波段。而影響較弱或含量較低的屬性，如pH與Fe，主要特征波段不明顯，分析方法不同時(shí)所對(duì)應(yīng)的特征波段會(huì)有所差別。從Fe和pH的相關(guān)關(guān)系圖中可以看出，對(duì)于主要特征波段不明顯的土壤成分，不同方法計(jì)算的波段相關(guān)系數(shù)排序不同，由此所確定的特征波段不同。因此使用三種分析方法計(jì)算波段與土壤成分的相關(guān)系數(shù)為指標(biāo)，采用過濾式方法篩選波段作為模型輸入變量并進(jìn)行模型精度分析。

表3和表4分別記錄了以PCC和GRD為過濾式篩選標(biāo)準(zhǔn)的具體建模指標(biāo)。從表中可以看出，PCC和GRD對(duì)四種成分的計(jì)算精度相似；除SOM外，其他三種成分所需的波段數(shù)目均較多，且與多維特征建模(表2)相比精度近似。表中N(Numbers)表示樣本數(shù)量，T(Threshold)表示相關(guān)系數(shù)閾值。

表3 線性相關(guān)法篩選特征及建模精度

表4 灰色關(guān)聯(lián)分析法篩選特征及建模精度

圖6為根據(jù)MBC值逐步添加特征得到的RF驗(yàn)證模型R2曲線。與基于PCC和GRD指標(biāo)分析的特征篩選方法相比，此方法計(jì)算得出的相關(guān)關(guān)系更適于提高預(yù)測(cè)模型的精度。從子圖中曲線的變化趨勢(shì)可以看出，對(duì)于SOM和SMC，波段數(shù)量的增加對(duì)模型精度的影響并不大，基于第6尺度分解系數(shù)的SOM與第7尺度分解系數(shù)的SMC驗(yàn)證模型R2基本位于0.9～0.95之間，且可以在較少波段的情況下獲得較優(yōu)的預(yù)測(cè)模型。

圖6 基于模型的相關(guān)關(guān)系篩選特征及建模精度

Fe的模型精度隨著波段數(shù)目的增加明顯降低。在可見光近紅外波長(zhǎng)范圍內(nèi)，由于SOM含量較高，在1 400 nm前對(duì)反射率具有主導(dǎo)影響，且SMC對(duì)光譜的吸收作用，干擾反射率曲線，使1 400 nm后的反射率曲線扭曲，致使Fe的特征波段信息被掩蓋。同時(shí)由于Fe包含的成分復(fù)雜，易受到外界環(huán)境及有機(jī)質(zhì)的影響[5]，較難估算。相關(guān)研究經(jīng)過100目研磨及MSC處理后的驗(yàn)證模型精度(R2=0.6，RPD=0.86)相比，圖6中Fe的建模當(dāng)篩選波段數(shù)目為8時(shí)，最佳模型R2為0.67，此時(shí)RPD為1.76(1.4<1.76<2)，模型較好可用于粗略估算。pH的預(yù)測(cè)模型在入選波段數(shù)目為41時(shí)達(dá)到最佳，此時(shí)驗(yàn)證模型R2=0.80，RPD=2.24。

3 結(jié) 論

多尺度SNV-CWT分解是一種能夠有效提高多種土壤理化指標(biāo)預(yù)測(cè)精度的方法。常規(guī)數(shù)據(jù)處理方法僅能在一定程度上提高某種成分的預(yù)測(cè)精度，而以Mexh為小波基函數(shù)的多尺度SNV-CWT分解方法可全面提高各種成分的信息提取質(zhì)量。與GS，F(xiàn)D和CR等7種常用的處理方法相比，四種成分預(yù)測(cè)模型精度在SNV-CWT分解下均有提高且優(yōu)于其他方法。通過對(duì)比處理前后反射率與土壤成分的PCC，MBC以及GRD，可以看出多尺度分解系數(shù)能夠有效地在不同尺度上放大成分的細(xì)節(jié)信息，提高土壤成分與反射率的相關(guān)系數(shù)。以MBC為指標(biāo)進(jìn)行的過濾式特征篩選為輸入變量的預(yù)測(cè)模型具有最佳的精度，其中SOM和SMC的R2=0.94，Mse分別為0.09%和1.1%，可以獲得極好的預(yù)測(cè)精度。Fe的R2=0.67，Mse=0.01%，RPD=1.76；pH的R2=0.80，Mse=0.1，RPD=2.24，可以達(dá)到用于定量預(yù)測(cè)的程度。土壤屬性復(fù)雜，在濕土壤情況下同時(shí)獲得多種土壤成分的含量，可以大大提高土壤質(zhì)量評(píng)估效率，快速進(jìn)行精準(zhǔn)農(nóng)業(yè)指導(dǎo)。