陳昊宇, 楊光, 韓雪瑩, 劉昕, 劉峰, 王寧

(內(nèi)蒙古農(nóng)業(yè)大學(xué)沙漠治理學(xué)院， 內(nèi)蒙古自治區(qū)風(fēng)沙物理與防沙治沙工程重點實驗室， 呼和浩特 010010)

精準(zhǔn)農(nóng)業(yè)作為目前農(nóng)業(yè)發(fā)展的主要方向，是一種基于信息和知識管理的現(xiàn)代化生產(chǎn)系統(tǒng)，主要是通過3S(GPS、GIS和RS)技術(shù)與現(xiàn)代農(nóng)業(yè)相結(jié)合，最大限度地提高農(nóng)業(yè)生產(chǎn)力。所以快速、無損、精確地獲取土壤中水分、養(yǎng)分的空間分布成為了實現(xiàn)精準(zhǔn)農(nóng)業(yè)的關(guān)鍵環(huán)節(jié)，近年來，光譜分析在土壤化學(xué)分析領(lǐng)域得到了迅猛發(fā)展，為實現(xiàn)土壤養(yǎng)分的快速診斷提供了新思路[1]。有機質(zhì)是土壤養(yǎng)分供應(yīng)能力和肥力的重要指標(biāo)之一，在全球碳循環(huán)中發(fā)揮著重要作用。因此，快速準(zhǔn)確地估測土壤有機質(zhì)含量對于發(fā)展精準(zhǔn)農(nóng)業(yè)具有重要意義[2]。

傳統(tǒng)的土壤有機質(zhì)測定方法雖然精度比較高，但周期較長、成本較高，只能達(dá)到瞬測量，很難進(jìn)行長時間大面積測量。高光譜遙感具有波段多、波段窄、信息豐富和實時高效等特點，為快速測量土壤有機質(zhì)含量提供了一種新的方法和手段[3]。

國內(nèi)外已經(jīng)有大量研究表明，通過對光譜數(shù)據(jù)進(jìn)行不同的數(shù)學(xué)變換(主要通過對光譜進(jìn)行倒數(shù)、對數(shù)、微分、平方根、吸收峰深度、包絡(luò)線去除等方法)可以有效提高光譜數(shù)據(jù)與土壤有機質(zhì)含量之間的相關(guān)系數(shù)，有效篩選出光譜信息中的敏感波段[4]?，F(xiàn)在各學(xué)者主要將研究重心放到了模型建立上[5]，普遍運用的線性模型有多元逐步回歸與偏最小二乘回歸[6]；常見的非線性模型包括BP神經(jīng)網(wǎng)絡(luò)[7]、支持向量機[8]、決策樹[9]等，而且隨著非線性模型算法的逐步改良與完善，在土壤有機質(zhì)含量估算中已經(jīng)成為不可取代的一部分。隨著小波算法的改進(jìn)與發(fā)展，最初僅運用于植物葉綠素、冠層成分含量預(yù)測中[10-11]，目前已成為土壤養(yǎng)分預(yù)測的熱點問題[12-13]，連續(xù)小波變換是目前被廣泛應(yīng)用的一種方法。王祥浩[14]選擇土地裸露地區(qū)為樣區(qū)，利用神經(jīng)網(wǎng)絡(luò)算法對光譜連續(xù)小波變換、一階導(dǎo)數(shù)、對光譜的平均值處理、光譜背景及深度4種方法建模，模型結(jié)果表明，小波變換方法得到的神經(jīng)網(wǎng)絡(luò)模型精度最高；包青嶺等[15]選擇渭干河-庫車河三角洲具有代表性的干旱區(qū)綠洲為研究區(qū)，對光譜進(jìn)行8層分解，結(jié)果表明小波變換不同分解層，從低頻到高頻范圍內(nèi)與土壤有機質(zhì)含量的相關(guān)性呈現(xiàn)先減后增的趨勢，結(jié)合隨機森嶺模型可以對干旱區(qū)土壤有機質(zhì)含量進(jìn)行有效的估算；王延倉等[16]以北京東部區(qū)潮土為例，對不同梯度重采樣的光譜進(jìn)行連續(xù)小波變換后，利用偏最小二乘法建立模型，結(jié)果表明連續(xù)小波分析算法可深入挖掘土壤光譜內(nèi)的有益信息，提升對有機質(zhì)含量的估測能力，與土壤高光譜反射率相比，經(jīng)連續(xù)小波技術(shù)處理后，模型精度得到了有效的提升；葉紅云等[17]同樣針對干旱區(qū)土壤，通過對兩種常用光譜變換R′、Ln(1/R)進(jìn)行連續(xù)小波變換建立偏最小二乘模型，結(jié)果表明連續(xù)小波變換不會因人類干擾程度的提高而使模型精度大幅度降低，更加適用于干旱區(qū)有機質(zhì)含量的預(yù)測；林鵬達(dá)等[18]通過解決黑土有機質(zhì)高光譜野外反演的困難，同樣證明了連續(xù)小波變換可有效提升模型精度。小波技術(shù)在土壤有機質(zhì)高光譜反演研究中逐漸趨于成熟，但目前學(xué)者的研究多數(shù)都在同一土壤類型下或同一區(qū)域內(nèi)，對于不同土壤類型及土地利用下土壤有機質(zhì)高光譜反演是否存在影響的研究目前并不多。本文研究區(qū)內(nèi)土壤類型主要包括3類：沙壤土、栗鈣土、鹽堿土，且部分區(qū)域土壤鹽漬化程度嚴(yán)重，導(dǎo)致土壤養(yǎng)分空間分布上存在較大差異，取樣表層土地利用類型主要包括：耕地、林地、草地、鹽漬地、荒地。

通過對原始光譜(R)、原始光譜倒數(shù)(1/R)、原始光譜對數(shù)(LnR)以及原始光譜一階微分(R′)4種不同情況進(jìn)行連續(xù)小波變換，利用BP神經(jīng)網(wǎng)絡(luò)以及支持向量機2種模型，探究了不同土壤類型與不同土地利用類型下是否會對土壤有機質(zhì)高光譜反演模型產(chǎn)生影響，小波變換前后土壤有機質(zhì)反演模型的精度，旨為區(qū)域土壤有機質(zhì)含量監(jiān)測及實現(xiàn)精準(zhǔn)農(nóng)業(yè)提供理論與技術(shù)支持。

1 材料與方法

1.1 研究區(qū)概況

托克托縣隸屬于內(nèi)蒙古自治區(qū)呼和浩特市，位于自治區(qū)中部、大青山南麓、黃河上中游分界處北岸的土默川平原上(圖1)。地理坐標(biāo)東經(jīng)111°2′30″—111°32′21″、北緯40°5′55″—40°35′15″，總面積1 409.67 km2，平均海拔1 117 m，屬于溫帶大陸性干旱氣候，年均氣溫7.3 ℃，年均降雨362 mm。托克托縣耕地總面積達(dá)400 km2，其中古城鎮(zhèn)、新營子鎮(zhèn)和五申鎮(zhèn)的耕地較多，占全縣耕地面積的60%以上[19]，主要作物包括小麥、玉米、莜麥。工農(nóng)業(yè)及生產(chǎn)生活用水主要來源于大黑河和黃河水資源，整個地形以大黑河為軸，呈現(xiàn)由丘陵向平原過渡的趨勢，地勢為東南高、西北和西南低。東南向西北土壤類型依次為栗鈣土、砂壤石灰性沖積土、鹽漬化石灰性沖積土[20]，土壤類型的不同導(dǎo)致土壤養(yǎng)分存在差異性分布。植被類型從西向東依次為草甸草原、干草原和退化灌叢草原分布。以Landsat8OLI影像為基礎(chǔ)數(shù)據(jù)源，運用人工目視解譯與BP神經(jīng)網(wǎng)絡(luò)分類法得到托克托縣2019年7月份土地利用數(shù)據(jù)，其中耕地面積最大為730.12 km2，占51.79%；林草地338.7 km2，占24.02%；鹽堿地141.1 km2，占10.00%。詳細(xì)土地利用空間分布見圖1。

圖1 土樣采集點及土地利用空間分布

1.2 研究方法

1.2.1土樣采集與處理土壤樣本點均勻地分布在托克托縣境內(nèi)，采集方法為五點采樣法，采集深度為0—20 cm，共采集120個點。采集的土樣置于通風(fēng)干燥室內(nèi)進(jìn)行自然風(fēng)干、研磨，過10目篩，進(jìn)行土壤光譜測定；過100目篩，采用重鉻酸鉀外加熱法進(jìn)行土壤有機質(zhì)含量測定。

1.2.2光譜測量及光譜處理土壤光譜于暗室內(nèi)測量，采用SVC HR-1024(北京東方佳氣科技有限公司)便攜式光譜儀，光譜范圍在350～2 500 nm。在350～1 000 nm波段之間光譜分辨率≤3.5 nm；在1 000～1 850 nm波段之間，光譜分辨率≤9.5 nm；在1 850～2 500 nm波段之間，光譜分辨率≤6.5 nm。光源采用與太陽光接近的50 W鹵素?zé)?，將土壤樣品放入? cm、寬10 cm的黑色器皿內(nèi)，用直尺將土壤表面刮平，探頭距離土樣10 cm，光源距離土壤表面30 cm，天頂角為15°。測量前用白板進(jìn)行標(biāo)定，每個土樣采集5條光譜作為該土樣的光譜數(shù)據(jù)。

由于受噪音與儀器暗電流的的影響，導(dǎo)致光譜數(shù)據(jù)混入噪音等信息，因此刪除350～399 nm和2 400～2 500 nm的波段，采用五點平滑法對光譜進(jìn)行平滑處理，并將光譜重采樣至5 nm，同時對原始光譜(R)進(jìn)行一階微分(R′)、倒數(shù)(1/R)、對數(shù)(LnR)等傳統(tǒng)數(shù)學(xué)變換。

1.2.3連續(xù)小波變換采用連續(xù)小波變換，并用Mexh小波母函數(shù)對原始光譜、原始光譜的倒數(shù)、對數(shù)、一階微分進(jìn)行10層小波變換，生成一系列小波系數(shù)。

(1)

式中，a為伸縮因子，b為平移因子，λ為土壤高光譜數(shù)據(jù)的波段數(shù)。

(2)

式中，f(λ)為土壤光譜反射率，小波系數(shù)Wf(a，b)包含二維，分別為波長(350～2 500)與分解尺度(1，2，3…10)， 故小波系數(shù)行為尺度數(shù)，列為波長數(shù)的矩陣[16]。

1.2.4模型及精度驗證采用BP神經(jīng)網(wǎng)絡(luò)與支持向量機模型(support vector machine，SVM)建立土壤有機質(zhì)預(yù)測模型，支持向量機采用線性核函數(shù)，相對于徑向基函數(shù)(radial basis function, RBF)來說計算高效,不易過擬合。BP神經(jīng)網(wǎng)絡(luò)的迭代次數(shù)設(shè)置為1 000，學(xué)習(xí)率0.01，訓(xùn)練的均方根誤差(root mean square error，RMSE)小于0.001。

依據(jù)相關(guān)系數(shù)篩選的特征波段以及小波系數(shù)作為自變量，土壤有機質(zhì)含量為因變量，分別建立模型，模型精度采用決定系數(shù)(R2)、均方根誤差(RMSE)、相對分析誤差(relative percent deviation，RPD)以及1∶1線共同評價。R2表征模型的穩(wěn)定性，越接近于1模型越穩(wěn)定，擬合程度越好。均方根誤差(RMSE)用來檢驗?zāi)Ｐ偷念A(yù)報能力，RMSE越小則表明模型的估測能力越好。RPD是樣本的標(biāo)準(zhǔn)差與RMSE的比值，RPD<1.4時，模型無法對樣品進(jìn)行預(yù)測；1.4≤RPD<2時，模型效果一般，可以用來對樣品進(jìn)行粗略評估；RPD≥2時，模型具有極好的預(yù)測能力。1∶1線表示實測值與預(yù)測值構(gòu)成的點偏離y=x線的程度[21]。

2 結(jié)果與分析

2.1 土壤有機質(zhì)含量統(tǒng)計分析

建模樣品集、不同土地利用方式、不同土壤類型下土壤有機質(zhì)含量描述性統(tǒng)計見表1。本研采樣點內(nèi)土地利用方式主要包括林地、草地、耕地、鹽漬地，土壤有機質(zhì)在草地內(nèi)均值含量最大(0.80%)，其次為林地(0.72%)、耕地(0.67%)、鹽漬地有機質(zhì)含量最低(0.63%)；土壤有機質(zhì)含量最大值位于耕地(1.28%)，最小值位于林地(0.19%)。采樣點內(nèi)主要土壤類型為栗鈣土、沙壤土、鹽堿土，沙壤土有機質(zhì)含量最高(0.77%)，其次為鹽堿土(0.68%)和栗鈣土(0.67%)，土壤有機質(zhì)含量最大值位于沙壤土內(nèi)(1.28%)，最小值位于鹽堿土內(nèi)(0.19%)。

表1 土壤有機質(zhì)含量描述性統(tǒng)計結(jié)果

2.2 土壤反射光譜特征

對R、1/R、LnR、R′進(jìn)行小波變換，變換結(jié)果如圖2所示，R、1/R、LnR光譜曲線較為平滑，分解曲線隨波峰波谷變化.R′其光譜曲線并不規(guī)則存在較多波峰波谷，分解小波系數(shù)與前三者不同。R、1/R、LnR、R′分解后，小波系數(shù)均隨分解尺度的增加而增加，同時可以看出，由Mexh小波母函數(shù)進(jìn)行的連續(xù)小波變換，對于光譜波峰與波谷有較高的敏感性，對于放大、挖掘光譜信息有著顯著的作用。

圖2 連續(xù)小波變換光譜特性

2.3 相關(guān)性分析

2.3.1不同導(dǎo)數(shù)變換光譜與土壤有機質(zhì)含量相關(guān)性土壤有機質(zhì)含量與光譜相關(guān)性曲線及敏感波段見圖3。R與土壤有機質(zhì)含量呈負(fù)相關(guān)關(guān)系(相關(guān)系數(shù)r=-0.463)，主要集中于735～780 nm處波段；1/R與土壤有機質(zhì)的相關(guān)性則與R相反，呈正相關(guān)關(guān)系(r=0.462)，集中于600～800 nm與1 800～2 200 nm處波段；LnR的相關(guān)性曲線圖與R相關(guān)性曲線類似，總體呈現(xiàn)負(fù)相關(guān)關(guān)系，相關(guān)系數(shù)(r=-0.465)，主要集中于745～795 nm處的波段；R′相關(guān)性在500 nm(r=-0.589)與1 400 nm(r=-0.411)處為負(fù)相關(guān)，在8 00 nm(r=0.408)與1 380 nm(r=0.412)處為正相關(guān)，相關(guān)系數(shù)曲線變換趨勢與前三者不同，呈無規(guī)律變化。

圖3 土壤光譜相關(guān)性曲線及敏感波段

2.3.2不同分解尺度小波系數(shù)與土壤有機質(zhì)含量的相關(guān)性圖4為不同光譜變換方式經(jīng)過連續(xù)小波變換后與土壤有機質(zhì)含量的相關(guān)系數(shù)矩陣圖，其中紅色代表相關(guān)性高的區(qū)域，藍(lán)色代表相關(guān)性低的區(qū)域。R在800～1 000、1 400～1 600 nm處相關(guān)性明顯增加，在500、800、2 200 nm波段處相關(guān)系數(shù)達(dá)到最大值(r=0.667)；1/R在800～1 200 nm處相關(guān)系數(shù)達(dá)到最大值(r=0.552)，在2 400～2 500 nm處相系數(shù)達(dá)到0.4，受噪音和儀器本身的影響，此波段的相關(guān)系數(shù)不進(jìn)行相關(guān)性參考；LnR在分解尺度1下相關(guān)性較低，在2～10尺度下，相關(guān)性出現(xiàn)最大值(r=0.664)；R′相關(guān)性主要集中在500～900、1 200～1 600、2 100～2 300 nm處。篩選的敏感波段與尺度如表2所示。有效的光譜信息主要存在于低分解尺度，隨分解尺度的增加呈遞減趨勢，相關(guān)性最大值較未處理前分別增加了0.204、0.09、0.199、0.252，對于挖掘潛在光譜信息有著重要意義。

圖4 土壤有機質(zhì)與小波系數(shù)相關(guān)性

表2 篩選的敏感波段

2.4 土壤有機質(zhì)高光譜模型建立

2.4.1BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型采用BP神經(jīng)網(wǎng)絡(luò)構(gòu)建反演模型，結(jié)果如表3所示。未進(jìn)行連續(xù)小波變換處理的模型中，BP-R與BP-R′效果較好，R2分別為0.69和0.73，RPD為1.45與1.53，模型能粗略估算土壤有機含量，BP-LnR與BP-1/R樣本外預(yù)測能力較差，同時RPD未達(dá)到1.4以上，不能對土壤有機質(zhì)未能進(jìn)行有效預(yù)測；連續(xù)小波變換處理之后的模型，僅BP-CWT-1/R模型RPD未達(dá)到預(yù)測水平，其余3種模型R2與RPD較未處理前均有所增加，RMSE均減少，其中BP-CWT-LnR模型預(yù)測效果較好，RPD達(dá)到2.12可以有效地對土壤有機質(zhì)進(jìn)行預(yù)測。將BP-CWT處理的4個模型的實測值與預(yù)測值進(jìn)行1∶1線分析。由圖5可知，除BP-CWT-1/R模型外，其余模型的實測值與預(yù)測值樣點基本分布在1∶1線附近，BP-CWT-LnR效果較為明顯，且估算精度高，可較好地進(jìn)行土壤有機質(zhì)含量的估算。

表3 土壤有機質(zhì)BP神經(jīng)網(wǎng)絡(luò)估測模型結(jié)果

圖5 BP-CWT模型土壤實測值與預(yù)測值對比

2.4.2支持向量機預(yù)測模型SVM構(gòu)建反演模型，結(jié)果如表4所示。未經(jīng)過連續(xù)小波處理的光譜特征波段未能較好地對土壤有機質(zhì)進(jìn)行預(yù)測反演，經(jīng)過CWT后模型SVM-CWT-R與SVM-CWT-R′預(yù)測結(jié)果較之前有較大的提升，R2分別達(dá)到了0.50與0.56，二者RPD均達(dá)到1.4以上，可以粗略地對土壤有機質(zhì)進(jìn)行預(yù)測。同時根據(jù)圖6，SVM-CWT模型進(jìn)行1：1線分析，二者實測值與預(yù)測值分布情況在4種模型下較好，雖然模型SVM-CWTLnR分布同樣較為集中，但其樣本外預(yù)測情況較差(RPD=1.38)，綜合考慮不對其進(jìn)行土壤有機質(zhì)預(yù)測。結(jié)合表3和表4的結(jié)果分析，連續(xù)小波變換能夠有效地提升模型精度與模型泛化能力，對于光譜信息挖掘有著重要意義，BP神經(jīng)網(wǎng)絡(luò)與支持向量機對CWT-R與CWT-R′都能夠提升R2減少RMSE，可對土壤有機質(zhì)做出較好的預(yù)測。雖然BP神經(jīng)網(wǎng)絡(luò)與支持向量機在處理非線性回歸問題中有較強的能力，但本身模型中存在不穩(wěn)定性，對模型的環(huán)境設(shè)置同樣要求較高，所以未能對所有數(shù)據(jù)集進(jìn)行良好的預(yù)測。

表4 土壤有機質(zhì)支持向量機估測模型結(jié)果

圖6 利用SVM-CWT模型土壤實測值與預(yù)測值的對比

3 討論

本研究采用連續(xù)小波變換對光譜進(jìn)行處理，用BP神經(jīng)網(wǎng)絡(luò)與支持向量機(SVM)兩種模型對土壤有機質(zhì)含量進(jìn)行反演預(yù)測。未經(jīng)過連續(xù)小波變換前，R、1/R、LnR、R′與土壤有機質(zhì)的相關(guān)系系數(shù)最大值分別為-0.463、0.462、-0.465、0.589，可以看出，R′與土壤有機質(zhì)的相關(guān)系數(shù)最高，與吳倩等[22]、張新樂等[23]的研究結(jié)果相同；經(jīng)過連續(xù)小波變換后，CWT-R、CWT-1/R、CWT-LnR、CWT-R′相關(guān)系數(shù)最大值分別為0.667、0.552、0.664、0.662，較之前分別增加了0.20、0.09、0.19、0.07。王延倉等[1]、于雷等[4]、葉紅云等[17]等同樣證明連續(xù)小波變換可有效提高與土壤有機質(zhì)含量的相關(guān)系數(shù)。不同分解尺度對于光譜數(shù)據(jù)的深度挖掘有著重要意義，本研究只利用Mexh小波母函數(shù)進(jìn)行處理，未對其他函數(shù)進(jìn)行考慮，分解層數(shù)同樣是根據(jù)前人經(jīng)驗所得[4，10]，小波技術(shù)的研究與發(fā)展仍然有很大的探索空間。

相對于兩種模型來看，未進(jìn)行連續(xù)小波處理的支持向量機模型中，只有SVM-R′模型R2最高達(dá)到0.43，其余三者均未到達(dá)0.4。綜合多種模型評價方法，由于其RPD未達(dá)到1.4以上，無法對土壤有機質(zhì)含量進(jìn)行預(yù)測。經(jīng)過連續(xù)小波處理后，各模型的R2有明顯提高，其中SVM-CWT-R與SVM-CWT-R′模型效果較好，R2分別提高了0.29、0.13，RPD達(dá)到1.62與1.53實現(xiàn)了對土壤有機質(zhì)有效的預(yù)測，但預(yù)測結(jié)果較BP神經(jīng)網(wǎng)絡(luò)較低。在BP神經(jīng)網(wǎng)絡(luò)預(yù)測模型中，未進(jìn)行連續(xù)小波變換前，BP-R與BP-R′預(yù)測效果較好，R2達(dá)到0.69與0.73，RPD為1.45與1.53；進(jìn)行連續(xù)小波處理后，除SVM-CWT-1/R模型未到達(dá)預(yù)測效果，其余3種模型預(yù)測結(jié)果較之前均有明顯改善，可實現(xiàn)對土壤有機質(zhì)較好的預(yù)測，其中BP-CWT-LnR預(yù)測模型效果最佳R2達(dá)到0.76，RPD達(dá)到2.12。根據(jù)1：1線分析圖也可看出，其實測值與預(yù)測值分布較為集中，于雷等[4]、葉紅云等[17]、林鵬達(dá)等[18]同樣通過連續(xù)小波變換有效提升了模型的精度與泛化能力。

針對土壤有機質(zhì)高光譜反演研究中，姚聰[24]對耕層土壤通過BP神經(jīng)網(wǎng)絡(luò)與支持向量機模型，反演精度R2分別為0.42與0.67；葉紅云等[17]采用連續(xù)小波變換對干旱區(qū)土壤有機質(zhì)反演，模型精度R2=0.75、EMSE=0.71；謝文[25]在森林土壤有機質(zhì)反演研究中，BP神經(jīng)網(wǎng)絡(luò)模型R2=0.78、EMSE=0.77，支持向量機模型R2=0.87、EMSE=0.76。本研究對耕地、林草地、鹽堿地、栗鈣土、沙壤土、鹽漬土等不同土地利用類型與土壤類型進(jìn)行綜合反演，最佳反演模型為BP-CWTLnR，R2=0.76、EMSE=0.15、RPD=2.12，與前人研究的結(jié)果基本相符，證明通過連續(xù)小波變換處理，不同土壤類型與土地利用類型未對土壤反演模型精度產(chǎn)生影響。所以采用連續(xù)小波變換進(jìn)行光譜數(shù)據(jù)挖掘，采用BP-CWT-LnR神經(jīng)網(wǎng)絡(luò)建立反演模型，可對不同土地利用于土壤類型條件下土壤有機質(zhì)高光譜反演提供一定的理論支持與應(yīng)用價值。