馬 馳

(遼寧省交通高等?？茖W(xué)校 測繪工程系，遼寧 沈陽 110122)

0 引言

氮元素是土壤養(yǎng)分的重要組成部分，是植物生長的必要元素，土壤中氮元素的含量直接影響植物的長勢及農(nóng)作物產(chǎn)量。因此，精確、快速量測土壤氮元素含量對于改善土壤質(zhì)量、提高農(nóng)作物產(chǎn)量乃至精細(xì)農(nóng)業(yè)的實(shí)施均具有重要意義[1-2]。傳統(tǒng)的利用化學(xué)試劑測試土壤氮元素含量的方法耗費(fèi)大量的人力與物力，難以實(shí)現(xiàn)土壤中氮元素含量的實(shí)時(shí)監(jiān)測。遙感技術(shù)具有數(shù)據(jù)獲取時(shí)間短、遙感影像覆蓋范圍廣、所含信息量豐富等優(yōu)勢。近年來逐漸成熟的高光譜技術(shù)由于具有更多的光譜通道和更高的光譜分辨率，能更加精細(xì)地反映地物的光譜特征和光譜差異，為土壤成分快速監(jiān)測提供了更有效的數(shù)據(jù)源[3-4]。

近年來，國內(nèi)外諸多學(xué)者相繼展開了利用遙感技術(shù)探測土壤全氮含量等方面的研究。Brunet D等[5]研究表明，將土壤光譜數(shù)據(jù)進(jìn)行一階微分處理后可以有效地提升土壤氮元素含量的反演精度；盧艷麗等[6]利用ASD2500 高光譜儀在實(shí)驗(yàn)室測繪土壤樣品的光譜曲線，并將光譜反射率進(jìn)行對數(shù)的一階微分變換后與土壤全氮相關(guān)系數(shù)在556，1642， 2 491 nm處出現(xiàn)峰值；王世東等[7]利用FieldSpec 3型光譜儀測定永城市礦區(qū)土壤光譜曲線，對土壤光譜進(jìn)行一階微分、二階微分、連續(xù)統(tǒng)去除數(shù)學(xué)變換，并利用偏最小二乘分析的方法建立研究區(qū)土壤全氮含量的反演模型，模型的決定系數(shù)R2達(dá)到了0.92；方向等[8]使用OFS-1700型便捷式地物非成像光譜儀對黃山市黃山區(qū)和池州市石臺縣土壤樣品進(jìn)行光譜分析，對土壤光譜進(jìn)行一階導(dǎo)數(shù)、二階導(dǎo)數(shù)、對數(shù)、去趨勢校正等數(shù)學(xué)變換，并利用偏最小二乘分析方法建立研究區(qū)土壤速效氮含量的估測模型，模型的決定系數(shù)R2達(dá)到了0.94。王文俊等[9]以山西典型褐土為研究對象，利用Starter Kit光譜儀在實(shí)驗(yàn)室繪制土壤樣本的近紅外高光譜曲線，采用平均光譜曲線和平均光譜曲線的一階微分等方法對光譜數(shù)據(jù)進(jìn)行預(yù)處理，并利用偏最小二乘的方法建立研究區(qū)建立土壤全氮含量預(yù)測模型，研究結(jié)果表明，使用一階微分進(jìn)行建模能夠獲得更好的預(yù)測效果。

綜上，當(dāng)前利用遙感技術(shù)研究土壤中全氮含量主要集中于氮元素的光譜特征方面，所使用的光譜數(shù)據(jù)多為利用便攜式光譜儀在實(shí)驗(yàn)室采集的土壤光譜數(shù)據(jù)，將其應(yīng)用于土壤元素的遙感制圖仍需輔以其他遙感影像數(shù)據(jù)，限制了其實(shí)用性。2008年9月發(fā)射升空的我國自行研制的環(huán)境災(zāi)害小衛(wèi)星(HJ-1A)所搭載的超光譜成像儀(Hyperspectral Imaging Radiometer,HSI)，與多光譜遙感影像相比，擁有多達(dá)115個(gè)探測波段和精細(xì)的光譜分辨率(僅為4.3 nm)，在土壤成分探測等方面已獲得了廣泛應(yīng)用[10-13]。然而，由于HJ-1A高光譜影像波段過多而造成各波段間的數(shù)據(jù)冗余度較高、對其進(jìn)行數(shù)據(jù)處理的過程也較復(fù)雜；影像的空間分辨率較低(100 m)，難以反映地表細(xì)部信息。2016年6月升空的Sentinel-2A遙感衛(wèi)星搭載的多光譜成像儀，與Landsat TM、ETM等多光譜成像儀相比，具有更高的空間分辨率、光譜分辨率和更短的重訪周期，必將在土壤成分探測等方面呈現(xiàn)出巨大優(yōu)勢。

為此，本文試驗(yàn)以中國的環(huán)境災(zāi)害小衛(wèi)星高光譜影像(HJ-1A HSI)和Sentinel-2A遙感影像為數(shù)據(jù)源，輔以研究區(qū)土壤采樣的全氮含量實(shí)驗(yàn)室化驗(yàn)數(shù)據(jù)，利用SPSS軟件統(tǒng)計(jì)光譜反射率與全氮含量的相關(guān)系數(shù)，篩選HJ-1A HSI影像全氮含量的敏感波段，建立兩種遙感影像土壤全氮含量的反演模型；通過比較兩種遙感影像反演研究區(qū)土壤全氮含量的模型精度，分析兩種遙感影像在土壤全氮含量方面的探測能力，為Sentinel-2A遙感影像在土壤元素探測方面的應(yīng)用研究提供參考。

1 材料與方法

1.1 研究區(qū)概況

農(nóng)安縣隸屬吉林省長春市，東經(jīng)124°31′～125°45′，北緯43°55′～44°55′，屬典型大陸性季風(fēng)氣候，年平均溫度為4.9 ℃、年平均降雨量為500 mm。農(nóng)安縣地勢平坦，平均海拔約250 m，地貌為沖積、湖積平原區(qū)，土壤類型主要包括黑土、黑鈣土、草甸土、風(fēng)沙土、鹽堿土等。

1.2 土壤采樣及化驗(yàn)

2017年4月29日—4月30日在農(nóng)安縣進(jìn)行土壤采樣，共采集土壤樣品82個(gè)，其采樣點(diǎn)分布如圖1所示。土壤采樣過程中選擇裸土區(qū)域，利用GNSS接收機(jī)測量采樣位置的經(jīng)緯度坐標(biāo)，在采樣位置20 m×20 m范圍內(nèi)采集表層土壤的多個(gè)土樣并混合，收集約0.5 kg土壤樣品裝入土壤采集袋；將土壤樣品置于干燥陰涼處風(fēng)干，剔除土壤樣品中的植物莖葉與根須、昆蟲軀殼、小石子等雜物并過篩。土壤全氮含量測定采用凱氏蒸餾法。

圖1 土壤采樣點(diǎn)分布Fig.1 Distribution of soil sampling points

1.3 遙感數(shù)據(jù)的選取與預(yù)處理

選取HJ-1A HSI遙感影像6景(成像時(shí)間為2017年4月30日)、Sentinel-2A遙感影像2景(成像時(shí)間為2017年4月29日)，遙感影像的云覆蓋率均小于2%。將HJ-1A HSI和Sentinel-2A兩種遙感影像進(jìn)行FLAASH輻射校正和大氣校正，用以消除兩種遙感衛(wèi)星的傳感器在成像過程中大氣的輻射誤差；分別對HJ-1A HSI和Sentinel-2A遙感影像進(jìn)行幾何精校正、遙感影像的拼接和裁剪等工作，獲得覆蓋農(nóng)安縣的遙感影像；根據(jù)土壤采樣時(shí)利用GNSS接收機(jī)測量獲得的采樣點(diǎn)坐標(biāo)，在HJ-1A HSI和Sentinel-2A遙感影像中提取采樣點(diǎn)反射率數(shù)據(jù)。

1.4 相關(guān)性分析

將HJ-1A HSI和Sentinel-2A遙感影像的反射率與研究區(qū)土壤全氮含量逐波段進(jìn)行相關(guān)性分析如式(1)所示，用以提取研究區(qū)土壤全氮含量敏感波段，并獲得建立全氮含量反演模型的最佳組合波段。

(1)

為了降低遙感成像中的噪聲對土壤光譜的影響，提高HJ-1A HSI和Sentinel-2A遙感影像反射率與土壤元素含量的相關(guān)性[5-9，12-14]，本文將兩種遙感影像的反射率進(jìn)行適當(dāng)?shù)臄?shù)學(xué)變換，包括對數(shù)lnR、倒數(shù)1/R、一階微分R′等，并與研究區(qū)土壤全氮含量在SPSS軟件中進(jìn)行相關(guān)性分析，再篩選出研究區(qū)土壤全氮含量的敏感波段。本試驗(yàn)中遙感影像反射率的一階微分采用光譜差分的方法估算如下：

(2)

式中，R′i為第i波段反射率的一階微分；Ri為第i波段的反射率；Δλ為第i波段至第i+1波段的光譜間隔。

1.5 土壤全氮含量反演模型的建立與模型檢驗(yàn)

將82個(gè)土壤樣本中的72個(gè)作為建模樣本、其余的10個(gè)作為模型精度的檢驗(yàn)樣本。將相關(guān)性分析所提取的敏感波段反射率與農(nóng)安縣土壤全氮含量進(jìn)行逐步回歸分析，建立農(nóng)安縣裸土全氮含量的反演模型。為了對反演模型進(jìn)行檢驗(yàn)，試驗(yàn)引入模型的決定系數(shù)R2以衡量模型的穩(wěn)定性，引入模型的均方根誤差(Root Mean Square Erro，RMSE)以衡量模型的精度。選擇具有較大模型決定系數(shù)和較小的模型均方根誤差的土壤全氮含量反演模型作為研究區(qū)最優(yōu)反演模型，用以反演研究區(qū)裸土的全氮含量。

2 結(jié)果與分析

2.1 相關(guān)性分析與敏感波段的選取

將農(nóng)安縣土壤全氮含量與HJ-1A HSI影像反射率(反射率變換)逐波段進(jìn)行相關(guān)性分析如圖2和圖3所示，獲得農(nóng)安縣裸土全氮含量的敏感波段如表1所示。相關(guān)性分析結(jié)果顯示,研究區(qū)土壤的全氮含量與HJ-1A HSI影像反射率呈負(fù)相關(guān)，并在第28波段(中心波長為524 nm)達(dá)到峰值，為r=-0.499。將反射率進(jìn)行適當(dāng)?shù)臄?shù)學(xué)變換后，與研究區(qū)裸土全氮含量的相關(guān)性有所提升。其中，反射率的對數(shù)變換與研究區(qū)土壤全氮含量的相關(guān)系數(shù)在第92波段(中心波長為782 nm)達(dá)到峰值，為r=0.605。反射率(反射率變換)經(jīng)過一階微分后與農(nóng)安縣裸土全氮含量的相關(guān)性進(jìn)一步得到改善，HJ-1A HSI影像反射率的一階微分變換與土壤全氮含量的相關(guān)系數(shù)在第92波段(中心波長為782 nm)達(dá)到峰值，相關(guān)系數(shù)為r=0.628，遙感影像的反射率指數(shù)一階微分變換與研究區(qū)土壤全氮含量的相關(guān)系數(shù)在第92波段達(dá)到峰值，相關(guān)系數(shù)達(dá)到0.695。

(a) 反射率

(a) 反射率的一階微分

表1 建模波段反射率(反射率變換)與裸土全氮含量的相關(guān)系數(shù)

2.2 反演模型的建立

根據(jù)相關(guān)性分析結(jié)果，選擇建模波段建立研究區(qū)土壤全氮含量的反演模型。建模波段的選取原則為：波段反射率(反射率變換)與研究區(qū)土壤全氮含量相關(guān)性較好、波段所含信息量豐富；波段間相距較遠(yuǎn)、波段間的相關(guān)性較小、數(shù)據(jù)的冗余度較低。本試驗(yàn)依據(jù)建模波段的選取原則，參照HJ-1A HSI影像反射率(反射率變換)與農(nóng)安縣裸土全氮含量相關(guān)性分析結(jié)果，選取全氮含量反演模型的建模波段(見表1)。將隨機(jī)選取的72個(gè)建模樣本所對應(yīng)遙感影像的建模波段反射率(反射率變換)作為自變量，以建模樣本全氮含量的化驗(yàn)值為因變量，利用SPSS軟件進(jìn)行逐步回歸分析，建立農(nóng)安縣裸土區(qū)全氮含量的反演模型如表2所示。

表2 HJ-1A HSI土壤全氮含量的反演模型

對比研究區(qū)土壤全氮含量反演模型發(fā)現(xiàn)，利用HJ-1A HSI影像反射率建立的反演模型的決定系數(shù)R2為0.568，RMSE為0.296 g/kg ；利用反射率的數(shù)學(xué)變換建立的全氮含量反演模型較反射率模型，模型的決定系數(shù)普遍有所提升。其中，利用反射率的冪變換建立的土壤全氮含量反演模型，模型決定系數(shù)R2達(dá)到0.754，模型的RMSE為0.218 g/kg。反射率(反射率變換)經(jīng)過一階微分后建立的反演模型的決定系數(shù)普遍高于非微分形式反演模型的決定系數(shù)。其中，反射率冪變換的一階微分反演模型的決定系數(shù)R2達(dá)到了0.806，模型的RMSE減小為0.185 g/kg。

2.3 模型精度檢驗(yàn)

利用反射率冪的一階微分反演模型Y=2.45+0.517X15+4.062X85+7.281X92+2.575X105計(jì)算檢驗(yàn)樣本全氮含量的反演值，與檢驗(yàn)樣本全氮含量的化驗(yàn)值建立散點(diǎn)圖。散點(diǎn)圖4顯示，研究區(qū)土壤檢驗(yàn)樣本全氮含量的反演值和實(shí)測值較均勻的分布于1∶1直線兩側(cè)，其擬合模型為Y=1.062X-0.097，擬合模型的決定系數(shù)R2=0.912，模型的RMSE=0.162 g/kg。

圖4 反演值與實(shí)測值散點(diǎn)圖Fig.4 Scatter diagram of inversion value and measured value

相對誤差散點(diǎn)圖5顯示，10個(gè)檢驗(yàn)樣本中有6個(gè)檢驗(yàn)樣本的相對誤差δ在(-0.1,0.1)之間，占檢驗(yàn)樣本總數(shù)的60%，有2個(gè)檢驗(yàn)樣本的相對誤差δ在(-0.2,-0.1)或(0.1,0.2)之間，表明HJ-1A HSI反射率冪的一階微分反演模型Y=2.45+0.517X15+4.062X85+7.281X92+2.575X105具有較好的預(yù)測效果，可以應(yīng)用于研究區(qū)土壤全氮含量的反演。

圖5 相對誤差散點(diǎn)圖Fig.5 Scatter diagram of relative error

2.4 多光譜建模

根據(jù)表1中HJ-1A HSI影像的土壤全氮含量敏感波段的分析結(jié)果，參考Sentinel-2A影像各波段的中心波長及光譜范圍，選擇建模波段并采用回歸分析的方法，建立農(nóng)安縣裸土全氮含量的Sentinel-2A多波段反演模型，如表3所示。建模結(jié)果顯示，非微分模型中，Sentinel-2A遙感影像的對數(shù)模型精度較高，反演模型的決定系數(shù)R2=0.634，RMSE=0.235 g/kg；而一階微分反演模型按反射率倒數(shù)模型、反射率模型、反射率指數(shù)模型、反射率對數(shù)模型、反射率冪模型的精度順序遞增，且模型間的反演精度差異較大，即：反射率冪的一階微分反演模型的預(yù)測能力最強(qiáng)，模型決定系數(shù)R2最大，達(dá)到0.768，模型的RMSE最小，為0.195 g/kg；反射率倒數(shù)的一階微分反演模型的預(yù)測能力最差，模型的決定系數(shù)R2僅為0.541。

表3 Sentinel-2A土壤全氮含量的反演模型

利用Sentinel-2A遙感影像反射率冪的一階微分反演模型Y=-0.254+2.6687X2+0.508X6-5.164X7-0.150X8A反演農(nóng)安縣土壤檢驗(yàn)樣本的全氮含量，與實(shí)驗(yàn)室土壤全氮含量的實(shí)測值建立散點(diǎn)圖,如圖6所示。由圖6可以看出，檢驗(yàn)樣本全氮含量的反演值與實(shí)測值較均勻的分布于1∶1直線兩側(cè)，擬合模型為Y=0.918X+0.264，擬合模型的決定系數(shù)R2=0.892，RMSE=0.196 g/kg。相對誤差散點(diǎn)圖(圖6)顯示，10個(gè)檢驗(yàn)樣本中有5個(gè)樣本的相對誤差δ在(-0.1,0.1)之間，占檢驗(yàn)樣本總數(shù)的50%，有4個(gè)檢驗(yàn)樣本的相對誤差δ在(-0.2,-0.1)或(0.1,0.2)之間，表明利用Sentinel-2A多光譜遙感影像反演研究區(qū)土壤全氮含量是可行的。

圖6 反演值與實(shí)測值散點(diǎn)圖Fig.6 Scatter diagram of inversion value and measured value

圖7 相對誤差散點(diǎn)圖Fig.7 Scatter diagram of relative error

2.5 土壤全氮含量制圖

根據(jù)表2和表3建模結(jié)果，選擇HJ-1A HSI及Sentinel-2A土壤全氮含量的最優(yōu)反演模型，反演研究區(qū)土壤全氮含量，如圖7和圖8所示。研究區(qū)全氮含量的反演結(jié)果圖顯示，利用HJ-1A HSI和Sentinel-2A遙感影像反演的研究區(qū)土壤全氮含量具有相似的空間分布，即：土壤全氮含量大于3 g/kg的區(qū)域主要集中于農(nóng)安縣東部和南部，土壤全氮含量小于2 g/kg的區(qū)域主要集中于農(nóng)安縣西部，農(nóng)安縣的中部和北部地區(qū)土壤全氮含量多介于2～3 g/kg。對研究區(qū)進(jìn)行調(diào)研，土壤采樣時(shí)發(fā)現(xiàn)，農(nóng)安縣東部、南部與我國東北黑土帶毗鄰，土質(zhì)肥沃、土壤全氮含量較高，農(nóng)安縣西臨松原市，為黑土、黑鈣土與鹽堿土過渡地帶，土壤含鹽量升高而全氮含量降低，在遙感影像中表現(xiàn)出較高的反射率。

圖8 HJ-1A HSI 全氮含量反演結(jié)果Fig.8 Total nitrogen inversion results chart of HJ-1A HSI

圖9 Sentinel-2A 全氮含量反演結(jié)果Fig.9 Total nitrogen inversion results chart of Sentinel-2A

3 討論

研究結(jié)果顯示，HJ-1A HSI遙感影像在可見光波譜范圍與研究區(qū)土壤全氮含量具有良好的相關(guān)性，且在紅光、綠光、藍(lán)光區(qū)域均存在峰值，與文獻(xiàn)[4，6，15]的研究結(jié)論相同或相近。將HJ-1A HSI遙感影像反射率進(jìn)行倒數(shù)、對數(shù)、指數(shù)、冪、一階微分等數(shù)學(xué)變換后可以有效改善與土壤全氮含量的相關(guān)性，與文獻(xiàn)[5-9，16-18]的研究結(jié)論相同。其中，HJ-1A HSI遙感影像反射率的對數(shù)變換與研究區(qū)土壤全氮含量的相關(guān)性在第92波段(中心波長為782 nm)達(dá)到最大值，為r=0.605；經(jīng)過一階微分變換后的遙感影像反射率后可以有效提升與農(nóng)安縣裸土全氮含量的相關(guān)性。其中，反射率指數(shù)一階微分與土壤全氮含量在第92波段達(dá)相關(guān)性最好，為r=0.695。利用回歸分析方法建立的HJ-1A HSI影像反射率冪的一階微分反演模型(模型的判定系數(shù)為R2=0.806，RMSE=0.185 g/kg)，其模型精度略優(yōu)于Sentinel-2A影像反射率冪的一階微分反演模型(模型判定系數(shù)R2=0.768，RMSE=0.195 g/kg)，究其原因，可能與HJ-1A HSI遙感影像的光譜分辨率遠(yuǎn)高于Sentinel-2A遙感影像有關(guān)，例如，HJ-1A HSI影像的全氮含量敏感波段第28波段，波譜范圍為5 nm，而與之相對應(yīng)的Sentinel-2A遙感影像綠光波段的波譜范圍為35 nm，HJ-1A HSI影像的精細(xì)光譜大大提高了其對土壤全氮的識別能力。

利用HJ-1A HSI遙感影像反射率冪的一階微分建立的研究區(qū)土壤全氮含量的最優(yōu)反演模型Y=2.45+0.517X15+4.062X85+7.281X92+2.575X105與參考HJ-1A HSI全氮含量敏感波段而建立的Sentinel-2A全氮含量的最優(yōu)反演模型Y=-0.254+2.6687X2+0.508X6-5.164X7-0.150X8A均有較高的模型精度，繪制的全氮含量反演結(jié)果圖具有相似的全氮含量空間分布。究其原因，可能與以下因素有關(guān)：① HJ-1A HSI遙感影像具有較高的光譜分辨率，大大提高了對土壤成分的識別能力；② 土壤采樣時(shí)間為4月29日—30日，與遙感影像的獲取時(shí)間幾乎同步，此時(shí)刻研究區(qū)地表已無冰雪及綠色植被，遙感影像能夠真實(shí)的反映土壤采樣時(shí)刻的裸土信息；③ 通過對HJ-1A HSI、Sentinel-2A遙感影像進(jìn)行了大氣校正，消除了傳感器成像時(shí)噪聲對反射率的影響；④ Sentinel-2A遙感影像光譜分辨率雖然低于HJ-1A HSI遙感影像，但Sentinel-2A遙感影像具有更高的空間分辨率，更容易顯示地表的細(xì)部特征。

然而，本試驗(yàn)中的檢驗(yàn)樣本全氮含量反演值與實(shí)測值之間仍然存在一定誤差，可能與以下因素有關(guān)：① 在利用FLAASH模型對兩種遙感影像進(jìn)行大氣校正時(shí)，由于所輸入的大氣能見度、大氣氣溶膠濃度、研究區(qū)地表高程等參數(shù)與實(shí)際存在一定誤差，使校正后影像的反射率與地表真實(shí)反射率存在誤差；② 本試驗(yàn)未考慮土壤中有機(jī)質(zhì)、鹽分等對土壤光譜的影響，使土壤全氮含量反演值存在一定誤差，將在以后研究中進(jìn)一步探究。

4 結(jié)束語

本文試驗(yàn)利用Sentinel-2A遙感影像代替HJ-1A HSI遙感影像反演農(nóng)安縣土壤全氮含量，獲得了以下結(jié)論：

(1) HJ-1A HSI遙感影像在可見光波段的反射率與農(nóng)安縣土壤全氮含量具有較強(qiáng)的負(fù)相關(guān)性，并在第28波段(中心波長為524 nm)相關(guān)性最好。

(2) 經(jīng)過適當(dāng)?shù)臄?shù)學(xué)變換后，HJ-1A HSI遙感影像反射率與農(nóng)安縣裸土全氮含量的相關(guān)性，較變換前有所改善；反射率(反射率變換)經(jīng)過一階微分變換后可以進(jìn)一步提升與裸土全氮含量的相關(guān)性。

(3) HJ-1A HSI遙感影像反射率經(jīng)過冪的一階微分變換后所建立的農(nóng)安縣土壤全氮含量反演模型，其模型精度最高，模型判定系數(shù)R2達(dá)到0.806；參考HJ-1A HSI遙感影像的土壤全氮含量敏感波段，建立Sentinel-2A遙感影像的土壤全氮含量反演模型，最優(yōu)模型的判定系數(shù)R2達(dá)到0.768，2個(gè)模型具有相近的模型精度，表明用于反演研究區(qū)土壤全氮含量的Sentinel-2A遙感影像可以代替HJ-1A HSI遙感影像。

(4) 利用HJ-1A HSI與Sentinel-2A遙感影像的土壤全氮含量最優(yōu)反演模型，反演研究區(qū)土壤全氮含量并制圖，均獲得了較好的反演效果。兩幅反演結(jié)果圖顯示，研究區(qū)土壤全氮含量具有相同的空間分布，即：農(nóng)安縣東部、南部土壤全氮含量較高，中部、西部土壤全氮含量較低。