基于多尺度SNV-CWT特征的黑土有機質(zhì)、水分、總鐵及pH值估測

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

吉林大學地球探測科學與技術(shù)學院，吉林 長春 130026

引 言

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

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

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

1 實驗部分

1.1 研究區(qū)概況

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

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

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

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

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

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

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

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

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

Zi=(Xi-m)/s

(1)

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

(2)

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

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

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

1.4 特征分析、建模與精度驗證

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

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

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

(3)

(4)

(5)

2 結(jié)果與討論

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

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

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

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

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

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

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

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

2.4 特征分析與建模

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

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

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

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

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

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

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

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

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

3 結(jié) 論

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