雷浩川 ， 劉尊方， 于曉晶， 曹金蓮， 吳笑天

(1.青海省基礎(chǔ)地理信息中心，青海 西寧 810016； 2.青海大學(xué)地質(zhì)工程系，青海 西寧 810016)

土壤是人類生產(chǎn)生活中重要的自然資源，是人類生存環(huán)境的重要載體。土壤作為土地資源的主要物質(zhì)基礎(chǔ)，其固有的不可再生性決定了承載容量的有限性[1,2]。在數(shù)字化農(nóng)業(yè)快速發(fā)展的時(shí)代，準(zhǔn)確、快速、動(dòng)態(tài)按需獲取土壤信息是現(xiàn)代精準(zhǔn)農(nóng)業(yè)的保證。土壤全氮含量是土壤氮素供應(yīng)狀況的重要指標(biāo)，也是衡量土壤肥力的重要因素之一。氮素是植被生長(zhǎng)過程中所需較多的營(yíng)養(yǎng)元素[3]。隨著現(xiàn)代農(nóng)業(yè)的發(fā)展，如何在有限的時(shí)間里獲取所需的土壤信息，進(jìn)而及時(shí)評(píng)價(jià)土地肥力，指導(dǎo)農(nóng)業(yè)生產(chǎn)的測(cè)土配方科學(xué)施肥，已是研究人員面臨的重大課題。如今，實(shí)時(shí)獲取土壤全氮含量，利用遙感影像進(jìn)行大范圍土壤全氮反演是現(xiàn)代精準(zhǔn)農(nóng)業(yè)的迫切要求。從 20 世紀(jì) 70年代開始，土壤光譜與氮素的關(guān)系已經(jīng)是國(guó)內(nèi)外許多學(xué)者的研究熱點(diǎn)[4]。徐麗華等[5]構(gòu)建PLSR模型反演土壤全氮含量，結(jié)果顯示土壤全氮含量實(shí)測(cè)值與預(yù)測(cè)值的相關(guān)系數(shù)r為0.672，預(yù)測(cè)效果較好。李燕麗等[6]研究遙感因子與土壤全氮含量變化的關(guān)系，研究結(jié)果顯示，土壤全氮含量預(yù)測(cè)的RMSE為0.35，表明土壤全氮含量反演模型的預(yù)測(cè)精度較高；王世東等[7]建立了PLSR、BPNN、RF等模型，并將PLSR與其他兩個(gè)模型相結(jié)合共同反演土壤全氮含量，結(jié)果表明，PLSR-BPNN模型的回歸決定系數(shù)R2可以達(dá)到0.92，說明模型具有較好的反演能力；祁亞琴等[8]構(gòu)建了土壤全氮含量的指數(shù)函數(shù)估算模型，該模型的預(yù)測(cè)效果最好，其回歸決定系數(shù)R2最高，為0.7982；楊越超等[9]構(gòu)建PLSR模型和BPNN模型反演了黑土土壤全氮含量，BPNN模型的土壤全氮含量預(yù)測(cè)精度比PLSR模型的預(yù)測(cè)精度高6.5%；張強(qiáng)等[10]基于NDI估算土壤全氮含量，其中預(yù)測(cè)效果最好的模型是指數(shù)函數(shù)(回歸決定系數(shù)R2為0.7982)；袁石林等[11]采用的最小二乘支持向量機(jī)(LS-SVM)的耦合模型對(duì)土壤氮素含量的反演精度較高；張娟娟等[12]研究了我國(guó)中、東部地區(qū)5種土壤，結(jié)果表明，經(jīng)過Norris濾波平滑之后，根據(jù)反射率的一階導(dǎo)數(shù)光譜構(gòu)建的模型預(yù)測(cè)效果最好；邱壑[13]研究反演了SOM、速效氮等養(yǎng)分含量，通過相關(guān)性分析提取相關(guān)波段，并根據(jù)實(shí)測(cè)值與預(yù)測(cè)值的偏差選出最優(yōu)的反演模型；肖云飛[14]的研究表明，隨機(jī)森林模型預(yù)測(cè)土壤全氮含量的精度高于PLSR模型，PLSR模型建模具有一定的局限性。

綜上，基于衛(wèi)星遙感影像的土壤全氮含量預(yù)測(cè)反演仍有很大的研究空間。本文以Landsat 5 TM多光譜遙感影像為數(shù)據(jù)源，分析了青海省西寧市大通回族土族自治縣的土壤表層全氮空間分布格局。利用PLSR模型和BP神經(jīng)網(wǎng)絡(luò)模型研究該地區(qū)的土壤全氮含量，為土壤碳氮循環(huán)研究、土壤質(zhì)量評(píng)價(jià)、糧食估產(chǎn)等工作提供數(shù)據(jù)支撐和理論參考。

1 材料與方法

1.1研究區(qū)概況青海省西寧市大通回族土族自治縣位于青海省東部河湟谷地(100°51′～101°56′E，36°43′～37°23′N)，祁連山南麓，湟水河上游北川河流域，是青藏高原和黃土高原的過渡地帶。海拔2 280～4 622 m，地勢(shì)西北高、東南低。研究區(qū)屬于高原大陸性氣候，年平均氣溫4.9 ℃，最值氣溫分別達(dá)到35.6 ℃和-26.1 ℃，8月份降水量最多。研究區(qū)行政區(qū)劃圖及采樣點(diǎn)分布圖如圖1所示。

圖1 研究區(qū)行政區(qū)劃圖及采樣點(diǎn)分布圖Fig.1 Administrative division and sampling point distribution of the study area

1.2樣本采集在研究區(qū)的農(nóng)用地塊以0～20 cm 深度采集土壤樣品73個(gè)，各約2 kg。在實(shí)驗(yàn)室將樣品自然風(fēng)干、磨細(xì)。每個(gè)樣品分為2份，分別用于光譜測(cè)試和土壤成分化學(xué)分析。土壤全氮質(zhì)量分?jǐn)?shù)的測(cè)定采用半微量凱氏定氮法[15]。53個(gè)樣品用于建立模型，剩余的20個(gè)樣品用于驗(yàn)證模型。53個(gè)建模樣品的全氮質(zhì)量分?jǐn)?shù)最大值為3.375 g/kg，最小值為0.639 g/kg，平均值為 1.857 g/kg，標(biāo)準(zhǔn)差為0.578 g/kg，變異系數(shù)為18%。

1.3影像獲取及預(yù)處理本研究所用的技術(shù)流程如圖2所示，Landsat 5 TM遙感影像數(shù)據(jù)從谷歌地球引擎(GEE)遙感云平臺(tái)快速獲取。觀測(cè)參數(shù)及波段特征如表1所示。影像成像時(shí)間分別是2008年9月11日和2008年9月20日，與土壤樣本的采集時(shí)間基本一致。9月左右農(nóng)作物已經(jīng)收割完畢，土壤表面沒有大量植被和作物遮擋，獲取土壤光譜信息較為方便。通過實(shí)地采樣養(yǎng)分化驗(yàn)數(shù)據(jù)和影像各個(gè)波段的相關(guān)性分析，得到土壤全氮含量的最佳反演波段。

圖2 技術(shù)路線圖Fig.2 Technology roadmap

表1 Landsat 5 TM觀測(cè)參數(shù)及波段特征Tab.1 Observation parameters of Landsat 5 TM and band characteristics

幾何校正、輻射校正(輻射定標(biāo)、大氣校正等)、影像去云、影像拼接和影像裁剪等是遙感影像預(yù)處理的主要過程。經(jīng)過影像預(yù)處理，原始影像的數(shù)字量化值(Digital Number，DN)轉(zhuǎn)化為地表真實(shí)的反射率(Surface Reflectance，SR)。谷歌地球引擎(GEE)遙感云平臺(tái)可以提供大氣表觀反射率TOA影像集和地表反射率SR影像集，為處理數(shù)據(jù)節(jié)省大量時(shí)間。影像預(yù)處理階段從GEE平臺(tái)使用函數(shù)ee.ImageCollection調(diào)用Landsat 5地表反射率影像(LANDSAT/L5/C01/T1_SR)，對(duì)于調(diào)用的影像只進(jìn)行影像去云操作、影像拼接和影像裁剪的預(yù)處理。依據(jù)前人的研究結(jié)論可知，僅僅通過原始波段的反射率提取特征波段有一定的困難，但是通過處理原始反射率可以有效提取全氮的特征波段(光譜預(yù)處理包括反射率的對(duì)數(shù)(logR)、倒數(shù)(1/R)、倒數(shù)的對(duì)數(shù)log(1/R)和對(duì)數(shù)的倒數(shù)1/logR等處理。)

土壤全氮含量與波段反射率及變換形式間的關(guān)系如表2所示，可以看出土壤全氮含量與波段原始反射率的相關(guān)性并不顯著，最大的負(fù)相關(guān)系數(shù)-0.534出現(xiàn)在B2波段，并且可見光波段的反射率與土壤全氮含量的相關(guān)性高于其他波段。不過經(jīng)過數(shù)學(xué)變換之后，原始反射率與土壤全氮含量的相關(guān)性明顯有所提升，尤其是經(jīng)過倒數(shù)處理的反射率與土壤全氮含量相關(guān)性最高，最大相關(guān)系數(shù)r達(dá)到了0.584。通過以上分析可知，與土壤全氮含量最敏感的是可見光波段(B1、B2和B3)。

表2 土壤全氮含量與波段反射率及變換形式間的關(guān)系Tab.2 Relationship between soil total nitrogen and band reflectance and transformation forms

1.4 模型構(gòu)建與驗(yàn)證

1.4.1 偏最小二乘回歸模型 Wold等[16]首次提出了偏最小二乘回歸(Partial Least Squares Regression，PLSR)算法。此算法涵蓋了普通的多元回歸分析、主成分分析和相關(guān)性分析，同時(shí)保留了這三種回歸分析的優(yōu)點(diǎn)，是以往線性回歸的優(yōu)化算法。先以TM影像多光譜波段的原始反射率或經(jīng)過數(shù)學(xué)變換后的反射率為基礎(chǔ)，通過KL主成分變換得到各個(gè)主成分，再以研究區(qū)土壤全氮含量為因變量，53個(gè)采樣點(diǎn)的地表真實(shí)反射率值或經(jīng)過數(shù)學(xué)變換之后的反射率為自變量，借助SPSS軟件建立PLSR模型，偏最小二乘回歸方程如式(1)所示：

Y=a1X1+a2X2+…+a2Xm

(1)

式中:Y代表土壤全氮含量，X代表回歸方程的自變量，系數(shù)ai對(duì)應(yīng)本文所研究的遙感影像多光譜波段的原始反射率和經(jīng)過數(shù)學(xué)變換之后的反射率，取自每個(gè)主成分對(duì)應(yīng)的系數(shù)值。對(duì)每個(gè)主成分進(jìn)行回歸分析，得到回歸決定系數(shù)R2和均方根誤差RMSE，取R2最大及RMSE最小時(shí)對(duì)應(yīng)的主成分建立PLSR回歸模型。

1.4.2 神經(jīng)網(wǎng)絡(luò)模型 神經(jīng)網(wǎng)絡(luò)(Neural Network)是以工程技術(shù)手段模擬人腦神經(jīng)元網(wǎng)絡(luò)結(jié)構(gòu)與功能的系統(tǒng),由大量簡(jiǎn)單的非線性處理單元組成。人工神經(jīng)網(wǎng)絡(luò)通過網(wǎng)絡(luò)中神經(jīng)元群體的相互作用體現(xiàn)自身的處理功能,可以處理模糊的、非線性的、含有噪聲的資料，特別適用于處理非線性問題，因而在模式識(shí)別、圖像處理和自動(dòng)控制等方面獲得了廣泛應(yīng)用。神經(jīng)網(wǎng)絡(luò)模型由輸入層、隱藏層和輸出層組成，其中輸入層和輸出層各1個(gè)[17]，隱藏層可以為多個(gè)。本文采用的BP神經(jīng)網(wǎng)絡(luò)模型是比較常用的神經(jīng)網(wǎng)絡(luò)模型之一，其結(jié)構(gòu)示意圖如圖3所示。

圖3 BP神經(jīng)網(wǎng)絡(luò)模型Fig.3 BP neural network model

本文構(gòu)建的BP神經(jīng)網(wǎng)絡(luò)模型包含3層，其中只包含1個(gè)隱藏層。輸入層的神經(jīng)元個(gè)數(shù)為通過相關(guān)性分析得到特征波段數(shù)目，以土壤全氮含量分析得到的特征波段確定神經(jīng)網(wǎng)絡(luò)的輸入層的神經(jīng)元的節(jié)點(diǎn)數(shù)，本研究中輸出層的神經(jīng)元數(shù)為1個(gè)，即土壤全氮含量，輸入層節(jié)點(diǎn)數(shù)為光譜波段個(gè)數(shù)，即6個(gè)神經(jīng)元數(shù)，隱藏層的節(jié)點(diǎn)數(shù)為10個(gè)。設(shè)置模型學(xué)習(xí)率為0.01，最大8 000次的迭代次數(shù)，BP神經(jīng)網(wǎng)絡(luò)模型的傳遞函數(shù)采用sigmoid函數(shù)計(jì)算隱藏層的輸出和最終輸出層的結(jié)果。BP神經(jīng)網(wǎng)絡(luò)模型的各個(gè)參數(shù)需要經(jīng)過反復(fù)試驗(yàn)后比較運(yùn)行結(jié)果而確定。

1.4.3 精度評(píng)價(jià)指標(biāo) 本文選用標(biāo)準(zhǔn)回歸評(píng)價(jià)中的回歸決定系數(shù)R2和誤差指數(shù)評(píng)價(jià)中的均方根誤差RMSE評(píng)價(jià)模型的穩(wěn)定性和預(yù)測(cè)精度。R2取值為[0，1]，模型的擬合程度越高，即實(shí)測(cè)值與預(yù)測(cè)值越接近，R2越接近于1；RMSE的值越小，說明模型預(yù)測(cè)值與實(shí)測(cè)值偏差較小，證明預(yù)測(cè)精度越高，土壤全氮含量反演誤差越小[18,19]。具體公式如(式2和式3)。

(2)

(3)

2 結(jié)果與分析

2.1 偏最小二乘回歸

最佳的主成分?jǐn)?shù)、模型建模樣本及驗(yàn)證樣本的回歸決定系數(shù)R2和均方根誤差RMSE如表3所示。本文構(gòu)建的兩個(gè)模型R2都在0.5～0.6，根據(jù)反射率倒數(shù)構(gòu)建的模型R2最大為0.604，RMSE為0.285，通過分析得出的全氮反演模型可以對(duì)土壤全氮含量進(jìn)行預(yù)測(cè)。

表3 全氮PLSR模型評(píng)價(jià)結(jié)果比較Tab.3 Comparison of evaluation results of total nitrogen PLSR model

通過以上分析繪制出土壤全氮含量實(shí)測(cè)值與根據(jù)反射率倒數(shù)構(gòu)建的模型反演預(yù)測(cè)值的散點(diǎn)圖，如圖4所示。根據(jù)各個(gè)樣本點(diǎn)的離散程度可知，以1/R為自變量構(gòu)建的模型預(yù)測(cè)效果最佳，對(duì)于土壤全氮含量的預(yù)測(cè)有較高的精度。

圖4 PLSR 1/R模型土壤全氮含量預(yù)測(cè)值與實(shí)測(cè)值比較Fig.4 Comparison between predicted and measured total nitrogen content of PLSR 1/R model

2.2 BP神經(jīng)網(wǎng)絡(luò)

全氮BP神經(jīng)網(wǎng)絡(luò)模型反演評(píng)價(jià)結(jié)果如表4所示。BP神經(jīng)網(wǎng)絡(luò)模型的預(yù)測(cè)能力較好，其均方根誤差RMSE最小值為0.246。根據(jù)4種反射率構(gòu)建的模型回歸決定系數(shù)R2均大于0.6，可以對(duì)土壤全氮含量進(jìn)行定量反演。驗(yàn)證樣本的回歸決定系數(shù)R2均略低于建模樣本的回歸決定系數(shù)R2，原因可能是采樣點(diǎn)的土壤全氮含量分布不均勻。總體而言1/R模型的預(yù)測(cè)能力高，運(yùn)用該模型繪制了土壤全氮含量預(yù)測(cè)值與實(shí)測(cè)值的散點(diǎn)圖(圖5)。相較而言，BP神經(jīng)網(wǎng)絡(luò)模型具有很大的優(yōu)勢(shì)，其預(yù)測(cè)能力和穩(wěn)定性都高于PLSR，這與肖云飛[14]的研究結(jié)果一致。

圖5 BP神經(jīng)網(wǎng)絡(luò)1/R模型土壤全氮含量預(yù)測(cè)值與實(shí)測(cè)值比較Fig.5 Comparison between predicted and measured total nitrogen content of BP neural network 1/R model

表4 全氮BP神經(jīng)網(wǎng)絡(luò)模型反演評(píng)價(jià)結(jié)果Tab.4 Inversion evaluation results of total nitrogen BP neural network model

2.3 研究區(qū)土壤全氮含量反演

利用所建立的偏最小二乘回歸模型和BP神經(jīng)網(wǎng)絡(luò)模型比較預(yù)測(cè)精度，結(jié)合遙感影像并選擇最優(yōu)模型反演研究區(qū)土壤全氮含量，以精度最高的BP神經(jīng)網(wǎng)絡(luò)模型對(duì)土壤全氮含量反演，如圖6所示。從圖中可以看出研究區(qū)西北部的土壤全氮含量較高，南部城區(qū)的土壤全氮含量相對(duì)較少。

圖6 研究區(qū)土壤全氮含量分布圖Fig.6 Distribution of soil total nitrogen content

3 討論與結(jié)論

土壤全氮含量是衡量土壤肥力的重要因素之一。實(shí)時(shí)獲取土壤全氮含量，利用遙感影像進(jìn)行大范圍土壤全氮反演，對(duì)于滿足現(xiàn)代精準(zhǔn)農(nóng)業(yè)的要求、為土壤質(zhì)量評(píng)價(jià)和糧食估產(chǎn)等工作提供基礎(chǔ)數(shù)據(jù)具有重要意義。近年來，國(guó)內(nèi)外許多學(xué)者關(guān)注土壤全氮研究以及土壤光譜與氮素的關(guān)系。高燈州等[20]研究了閩江口濕地的土壤全氮含量，發(fā)現(xiàn)在波長(zhǎng)0.5 μm附近土壤全氮含量的負(fù)相關(guān)系數(shù)最高；吳明珠等[21]研究亞熱帶地區(qū)土壤全氮含量與光譜反射率之間的相關(guān)性，結(jié)果表明在波長(zhǎng)0.35～2.5 μm光譜反射率與土壤全氮含量均呈現(xiàn)負(fù)相關(guān)關(guān)系。彭杰等[22]分析研究了4種土壤，所構(gòu)建的多元逐步回歸模型預(yù)測(cè)精度高于一元線性回歸模型，其中根據(jù)對(duì)數(shù)倒數(shù)的導(dǎo)數(shù)所構(gòu)建的模型土壤全氮含量預(yù)測(cè)值與實(shí)測(cè)值的回歸決定系數(shù)R2達(dá)到0.837 2。本文以大通縣Landsat 5遙感影像和采樣地?cái)?shù)據(jù)為基礎(chǔ)，基于多元線性回歸分析方法和主成分分析法，分別構(gòu)建了偏最小二乘回歸模型和BP神經(jīng)網(wǎng)絡(luò)模型，反演出整個(gè)研究區(qū)的土壤全氮含量空間分布狀況；利用統(tǒng)計(jì)軟件對(duì)采樣點(diǎn)表層土壤全氮含量與對(duì)應(yīng)點(diǎn)的TM影像各波段SR值進(jìn)行了相關(guān)分析，并得出了以下主要結(jié)論：

共有7個(gè)參與相關(guān)性分析的Landsat 5光譜響應(yīng)波段，其中B2波段最佳，對(duì)原始光譜反射率進(jìn)行數(shù)學(xué)變換預(yù)處理之后，B2波段的反射率1/R的相關(guān)系數(shù)r達(dá)到了0.584，P<0.01；對(duì)比偏最小二乘回歸模型和BP神經(jīng)網(wǎng)絡(luò)模型可知，非線性模型的預(yù)測(cè)精度優(yōu)于線性模型，根據(jù)反射率倒數(shù)構(gòu)建的BP神經(jīng)網(wǎng)絡(luò)模型回歸決定系數(shù)R2為0.792，均方根誤差RMSE為0.246，明顯優(yōu)于偏最小二乘回歸模型，預(yù)測(cè)能力較好。