HJ-1A HSI與Sentinel-2A遙感數據土壤全氮含量反演精度的對比研究

馬 馳

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

0 引言

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

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

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

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

1 材料與方法

1.1 研究區概況

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

1.2 土壤采樣及化驗

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

圖1 土壤采樣點分布Fig.1 Distribution of soil sampling points

1.3 遙感數據的選取與預處理

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

1.4 相關性分析

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

(1)

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

(2)

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

1.5 土壤全氮含量反演模型的建立與模型檢驗

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

2 結果與分析

2.1 相關性分析與敏感波段的選取

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

(a) 反射率

(a) 反射率的一階微分

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

2.2 反演模型的建立

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

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

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

2.3 模型精度檢驗

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

圖4 反演值與實測值散點圖Fig.4 Scatter diagram of inversion value and measured value

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

圖5 相對誤差散點圖Fig.5 Scatter diagram of relative error

2.4 多光譜建模

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

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

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

圖6 反演值與實測值散點圖Fig.6 Scatter diagram of inversion value and measured value

圖7 相對誤差散點圖Fig.7 Scatter diagram of relative error

2.5 土壤全氮含量制圖

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

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

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

3 討論

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

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

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

4 結束語

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

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

(2) 經過適當的數學變換后，HJ-1A HSI遙感影像反射率與農安縣裸土全氮含量的相關性，較變換前有所改善；反射率(反射率變換)經過一階微分變換后可以進一步提升與裸土全氮含量的相關性。

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

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