基于航空高光譜反演的黑土地質量綜合評價研究

趙寧博， 楊佳佳， 趙英俊， 秦凱*， 楊越超， 李明

(1.核工業北京地質研究院，遙感信息與圖像分析技術國家級重點實驗室，北京 100029；2.中國地質調查局沈陽地質調查中心，沈陽 110034)

土壤質量通常指土壤在生態系統中保持生物生產力、維持環境質量及促進動植物健康的能力[1]。土壤質量評價主要從土壤肥力質量、環境質量和健康質量三個維度[2-3]開展，評價指標涵蓋土壤物理、化學和生物學指標。土壤質量評價結果可以直接指導土地利用規劃、農業生產、土壤環境保護工作，同時也是污染土壤治理修復、土壤污染生態效應評價及地球化學災害預測研究的基礎[4]。

我國東北地區的黑土地性狀好、肥力高，非常適合農作物生長，是我國重要的糧食生產基地。但近年來黑土地腐殖層變薄、肥力下降、水土流失等現象已經引起了廣泛關注[5-7]，黑土地的退化將直接影響農業生態系統和農業可持續發展，因此必須科學及時地開展土壤質量評價工作，為黑土地的保護與可持續利用提供有效決策依據。

目前，我國黑土地的土壤質量評價主要采用地球化學調查手段，局限性表現為地面采樣耗費較多人力物力，樣品分析工作周期較長，在大比例尺調查時局限性表現的更為明顯。遙感技術憑借時效性強、數據覆蓋面廣等優勢，能夠對地球化學調查進行有效補充，逐漸在評價工作中發揮作用。根據遙感數據的特點，在土壤質量評價中主要有兩方面應用，一是利用多光譜影像提取植被指數(normalized differential vegetation index，NDVI)[8]、土壤退化指數(ratio vegetation index，RVI)[9]、土壤水分指數(difference vegetation index，DVI)[9]、地形因子[10]、土地利用[11]等信息；二是利用高光譜數據提取土壤有機質[12-15]、養分元素[16]、重金屬元素[17-18]、含鹽量[19]、質地[20]等指標。

相比于多光譜遙感，高光譜遙感提取土壤理化性質的能力更強，能夠獲得更為全面的土壤質量指標數據。目前，土壤質量調查中主要采用的是地面高光譜技術，數據為散點形式且工作周期較長，航空高光譜技術改善了地面調查方式下數據獲取周期長及覆蓋面不足的劣勢，且具有“圖譜合一”的特征，提高了在土壤質量評價中的實際應用能力。本文以黑龍江省海倫典型黑土區為例，利用航空高光譜數據從土壤肥力、土壤環境和農作物長勢幾個維度提取相關信息，構建土壤質量綜合評價模型，初步探討航空高光譜技術在黑土地質量評價中的應用能力和效果，為田塊尺度的黑土地質量快速評價提供技術支撐。

1 材料與方法

1.1 研究區概況

研究區位于我國黑龍江省海倫地區，坐標126°29′10″—127°4′40″E、47°31′10″—47°35′40″N。海倫地處松嫩平原東北端，小興安嶺西麓，地形為丘陵、漫崗，屬中溫帶大陸性氣候。

海倫地區位于我國東北黑土區的中心區域，主要土壤類型為黑土、暗棕壤和草甸土等，是黑土地研究領域中的代表性區域。該地區已經開展了1∶5萬多目標地球化學研究，航空高光譜與地球化學的有機結合能夠為黑土地關鍵帶的綜合建模提供支撐，便于開展黑土質量評價工作。

1.2 數據獲取與預處理

1.2.1航空高光譜數據 設備采用CASI-1500和SASI-600高光譜成像系統，由加拿大Itres公司生產。CASI-1500傳感器譜段范圍為380～1 050 nm，空間分辨率1.5 m，波段數72個；SASI-600傳感器譜段范圍為950～2 450 nm，空間分辨率為3.75 m，波段數100個。為了保證航空高光譜數據質量，挑選晴朗無云的天氣飛行，飛行時間設定為當天10：00—14：00。數據分兩期獲?。旱谝黄诖螌嶋H獲取日期為2018年5月29日，地面處于裸土期，用于反演土壤質量參數；第二期次實際獲取時間為2018年7月31日，為夏季農作物成熟期，用于反演農作物質量參數。

數據預處理主要包括輻射定標、幾何校正、大氣糾正和光譜重建。航空成像光譜測量系統配備了輻射定標軟件RCX(radiometric calibration xpress) 9.3.5.1和幾何校正軟件Geocor 3.0。首先利用RCX軟件進行輻射定標，然后利用Geocor軟件對原始GPS數據和基站數據進行差分處理，提取GPS時間并處理傳感器姿態數據，對傳感器姿態數據與GPS定位數據進行時間同步與集成，最后利用幾何校正信息對輻射定標后的數據進行幾何校正。然后在ENVI5.3軟件中繼續進行幾何精校正處理，利用高分2號影像作為參考數據，利用控制點校正的方法完成幾何精校正，進一步提升數據的幾何精度。

在此基礎上，利用ENVI5.3軟件中大氣輻射傳輸模型和地空回歸方法對航空高光譜遙感數據重建光譜。首先在ENVI的FLAASH模塊中利用大氣輻射傳輸模型進行校正，然后利用地面黑、白兩種定標布的光譜測量數據在ENVI的empirical line模塊中進行地空回歸校正，進一步消除因輻射定標、波段間相對定標、波段配準、大氣參數選取等誤差因素而造成的光譜誤差，最終獲得地物光譜反射率。

為了更全面地提取研究指標的光譜特征，利用ENVI5.3軟件擴展工具中的Image Derivate模塊反射率數據進行一階導數變換，Continuum Removal模塊進行去連續統變換。

本次數據反演及質量評價僅限于研究區內的旱田范圍，因此在高光譜影像中提取旱田范圍用于后續數據處理，地物分類提取方法采用ENVI軟件擴展工具中的隨機森林法。首先在影像內建立感興趣區(region of interest，ROI)，ROI均勻分布于全區范圍，保證覆蓋到各種地物類型，然后利用隨機森林模型對選取的ROI進行訓練，獲得研究區的地物分類數據，最后利用ENVI軟件中的Majority Analysis工具分類數據進行后處理，去除細小斑塊。

1.2.2地面數據 在航空數據獲取時同步開展地面調查，調查點包括38個土壤調查點及22個農作物調查點。地面獲取的數據指標類型即為參與土壤質量評價的指標類型，地面數據用于后續的航空高光譜反演建模及驗證。指標選取的依據包括：①以土地質量地球化學評價規范[21]為基礎，選取通用的土壤養分和環境指標，保證評價結果的科學性和通用性。②突出航空高光譜的優勢，結合評價目標拓展土壤理化性質和農作物長勢方面的指標。

根據上述指標選取依據，地面共獲取了包括土壤肥力、土壤環境和農作物長勢三個類別共15項指標： ①土壤肥力指標。包括有機質(SOM)、全氮(N)、全磷(P)、全鉀(K)、硒(Se)、鍺(Ge)、陽離子交換量(cation exchange capacit，CEC)；②土壤環境指標。土壤環境指標包括重金屬污染和土壤退化兩類。重金屬元素參照規范[22]選取了砷(As)、鉻(Cr)、鎘(Cd)、汞(Hg)、鉛(Pb)五種元素。土壤退化指標選取了全鹽量，代表土壤鹽堿化程度。③農作物長勢指標。包括葉綠素含量和葉面積指數，在第二期次農作物航空高光譜調查時同步測量。

1.2.3土壤指標化學分析 SOM采用硫酸亞鐵銨容量法[23]滴定，儀器為酸式滴定管；采用凱氏定氮法[24]測定N含量，儀器為上海市歐公司生產的全自動凱氏定氮儀；采用X射線熒光光譜法[25]測量P和K含量，儀器為荷蘭帕納科公司生產的Axios-mAX 波長色散型X射線熒光光譜儀；采用電感耦合等離子體發射光譜法[26]測量Se、Ge、Hg、As、Pb、Cd、Cr含量，儀器為美國賽默飛世爾公司生產的ELEMENT XR 等離子體質譜儀；采用乙酸銨法[27]測定CEC，儀器為上海赫田公司生產的電動離心機；采用浙江托普儀器公司生產的SPAD502PLUS葉綠素儀測定葉綠素含量；采用英國Delta-T公司生產的SUNSCAN冠層分析儀測定葉面積指數。

1.3 土壤質量參數反演

本次評價所采用的各類土壤質量參數數據均為航空高光譜數據反演所得。航空數據的光譜分辨率較高，通過不同土壤質量參數的光譜特征對其含量進行反演。將航空數據與地面采樣點進行空間疊加，在航空數據中以采樣點位置為中心，周圍3×3個像元范圍的光譜平均后作為該點的光譜數據。

采用偏最小二乘法建模。在建模時，特征波段選擇直接關系著模型的精度和穩定性，由于土壤成分和性質較為復雜，針對同一指標，土壤類型或區域不同時相應的特征波段也不盡相同。在選擇特征波段時，首先利用全部波段參與建模，然后根據模型系數曲線選擇峰值區域的波段再次建模，不斷優化調整特征波段的，最終根據建模和驗證精度選擇最優模型。

將總樣本數量的2/3用于建模，1/3用于驗證，根據樣本的空間分布均勻挑選建模和驗證樣本，以減少模型誤差。

1.4 土壤質量綜合評價

以土壤及農作物多項質量參數的航空高光譜反演數據為基礎，通過建立評價模型最終獲得研究區的土壤質量綜合評價指數。選取層次分析法(analytical hierarchy process)作為評價方法將土壤質量評價系統結構化，形成遞階層次結構模型，然后對每一層次中的各因素進行較客觀的逐對比較和判斷，最后通過排序結果來分析和解決問題。層次分析法所采用的軟件為YAAHP10.1。

1.4.1評價層次結構 根據選取的評價指標建立層次結構(圖1)，分為目標層、中間層和指標層三個層次。目標層為土壤質量，即評價工作的最終目標。中間層包括土壤肥力、土壤環境和農作物長勢三個大類，土壤肥力又包括大量元素和微量有益元素類別；土壤環境包括重金屬污染和土壤退化兩類。指標層即上文中選取的SOM、N、P、K等各項土壤養分、環境、農作物長勢類別中的具體指標。

圖1 評價模型層次結構

1.4.2權重計算 權重計算是建立評價模型的重要步驟之一。根據咨詢專家、參考現有規范和統計資料對判斷矩陣中各指標的重要性進行兩兩比較，按照9標度方法[28]進行打分，構建出比較判斷矩陣。根據建立的判斷矩陣，采用和積法分別計算出最大特征值所對應的特征向量，所求特征向量即為相對權重值。權重計算完成后需要計算一致性比率(consistent ratio，CR)，當CR=0時，表明權重分配具有完全一致性；當CR<0.1時,具有滿意一致性,說明權重分配合理；否則就要調整判斷矩陣，直到具有滿意的一致性為止。本次判斷矩陣的CR計算結果為0.017 6，通過一致性檢驗。最終權重計算結果見表1。

表1 模型權重計算結果

1.4.3評價單元的確定 通常根據評價數據比例尺的不同，評價單元可以按照規則網格、土地利用現狀調查圖斑、實際地塊等不同方式劃分，本研究憑借航空數據的高空間分辨率優勢，評價單元采用最高精度劃分方式，即研究區的實際地塊。地塊數據以研究區第二次全國土地調查數據為基礎，并與航空高光譜影像進行疊加后檢查地塊空間劃分方式的變化，經檢查后整體上地塊劃分無明顯變化，并對個別地塊發生改變的區域進行了修改完善。

1.4.4指標評分 通過制定打分標準對每種指標進行評分，分值區間為0～1分。評分標準制定以相關規范中的標準為基礎，并考慮不同指標特征進行科學劃分，具體評分標準如下。

①土壤肥力指標評價標準。肥力指標中SOM、N、P、K的評價標準參照規范[21]中的規定，評級由高至低包括五類等級：一等(1.0分)、二等(0.8分)，三等(0.6分)、四等(0.4分)、五等(0.2分)。需要說明的是，P和K采用的是全量，而不是有效態含量。速效磷和速效鉀是能夠被植物吸收的部分，與土壤肥力的關系比全量更為密切，但是通常情況下有效態含量占全量的比重很小，過低的含量會大大提高其光譜特征的獲取難度，從而影響反演精度。本研究的重點是在保證反演精度的前提下構建綜合評價體系，后續將進一步開展對速效磷和速效鉀等指標的高光譜反演研究。

有益元素Se和Ge的評價標準采用隸屬度函數模型。隸屬度函數模型針對不同類型評價目標包分為戒上型、戒下型或峰值型三類。Se和Ge均具有生物學的兩面性，在一定含量范圍內對人體和動植物生長有益，過量則會造成損害，基于該特性，如果研究區指標含量中存在明顯過量級別時宜采用峰值型函數，如果沒有過量級別存在則采用戒上型函數。Se和Ge元素的過量臨界值標準分別采用3.0和5.8 mg·kg-1[29]，根據研究區統計結果，Se最大值為0.43 mg·kg-1，Ge最大值為1.69 mg·kg-1，均低于過量的臨界值，因此隸屬度函數模型采用戒上型，通過隸屬度函數計算得到Se和Ge的評分。

目前CEC較為常用的評價標準分為三個等級，CEC≥20 cmol·kg-1代表保肥能力強；1017～20)，0.8分；三等(>13～17)，0.6分；四等(>10～13)，0.4分；五等(≤10)，0.2分。

②土壤環境指標評價標準。土壤重金屬元素As、Cr、Cd、Hg、Pb的分級標準同樣參照規范[21]。

土壤全鹽量數值通過電導率法測定，因此參照相關電導率數值的分級標準[30]制定評分規則。

③農作物長勢指標評價標準

農作物長勢指標包括葉綠素含量和葉面積指數，這兩種指標含量均與農作物長勢呈正相關，因此采用戒上型隸屬度函數模型進行計算評價。

1.4.5土壤質量分級 評價單元賦值包括兩個步驟：①在ArcGIS軟件中將評價所采用的矢量點文件與指標反演的柵格文件進行空間分析，將矢量點所在空間位置對應的柵格數值賦予該點；②將矢量點文件與地塊圖斑進行空間分析，同一地塊中所有矢量點數據取平均值作為該地塊的數值。

將所有評價指標的數據賦予地塊圖斑后，依照上文所述的相應評價標準進行單指標得分計算，然后依照如下計算方式獲得綜合評價得分。

式中，F綜為綜合評價得分，Fi為第i個指標評價得分，Wi為第i個指標權重。

綜合評價得分的分值區間為0～1.0，共分為五個等級：一等(優質)土壤，>0.9～1.0；二等(良好)土壤，>0.7～0.9；三等(中等)土壤，>0.5～0.7；四等(差等)土壤，>0.3～0.5；五等(劣等)土壤，≤0.3。

2 結果與分析

2.1 模型關鍵參數分析

2.1.1光譜數學變換形式對建模的影響 分別采用光譜反射率、反射率一階導數和反射率去連續統變換數據作為模型自變量參與計算，各指標對應的模型決定系數(R2)見表2。針對每一種指標，反射率一階導數變換對應的建模和驗證精度均遠低于反射率相對應的精度；基于反射率去連續統建立的模型中，N、Ge、Cr和Hg的建模精度優于反射率對應的模型精度，但是所有指標的驗證精度均低于反射率對應的精度。

表2 基于不同光譜形式建立的模型R2

一階導數和去連續統變換均有提升微弱光譜特征的作用，但此次建模效果并不如原始光譜反射率，尤其是基于一階導數變換的模型精度最低。原因是航空光譜受到大氣傳輸等因素的影響，數據信噪比要低于地面光譜，一階導數和去連續統變換在增強光譜特征的同時也提高了噪聲的影響，因此最終選用原始反射率數據進行建模。

2.1.2特征波段分析 首先利用全部波段參與運算，然后基于回歸系數曲線挑選具有代表性的峰值波段作為特征波段再次運算，根據模型精度進行波段的優化調整，直至模型達到最優。以SOM為例，根據圖2挑選了570、943、1 085、1 175、1 505、1 670、1 775、2 105、2 225和2 390 nm作為特征波段，利用上述特征波段再次建模并觀察波段回歸系數，最終將570 nm剔除波段組合后模型精度達到最優。根據上述方法對各指標進行特征波段選擇與模型優化，各指標的特征波段及最終模型精度見表3，R2均達到較高水平，說明優化后的模型精度得到提升。

圖2 SOM模型回歸系數

2.2 反演模型精度分析

評價指標反演的可靠性是后續土壤質量評價的基礎。從表3可以看出，SOM的建模R2為0.813，驗證R2為0.726，兩項均為所有指標中最高；N、P、K、Se、CEC的建模R2均超過了0.7，驗證R2均超過了0.6；養分指標中Ge元素的建模R2稍低，為0.667；As、Cr、Cd、Hg、Pb幾種重金屬元素及全鹽量的建模R2處于0.6～0.7區間，驗證R2處于0.5～0.6區間。葉綠素和葉面積指數兩種農作物長勢指標的建模R2均超過了0.7，驗證R2均超過了0.6。

表3 各指標特征波段與模型最終精度

從全部指標反演精度的橫向對比來看，整體上土壤養分和農作物長勢指標的模型精度要稍高于重金屬元素，一方面因為重金屬元素不像有機質等成分具有明顯的光譜特征，另一方面研究區內重金屬元素含量較低，一定程度上影響了通過土壤礦物等組分進行間接反演的效果。對比各指標的建模R2和驗證R2，沒有出現驗證R2顯著降低的情況，說明反演沒有明顯的過擬合現象。

2.3 反演結果分析

通過反演結果(圖3)觀察研究區各指標的空間分布狀態?？梢钥闯?，SOM高值區主要集中于愛民鄉，低值區主要位于海北鎮，其他區域處于中間狀態；N和SOM的空間分布規律接近，這是由于土壤中的氮素絕大多數是貯藏在有機質中的有機態含氮化合物，其次是被黏土礦物吸附的銨態氮、硝態氮和亞硝態氮，因此兩者關系密切；P高值區主要分布于研究區南部，向榮鄉-長發鎮-東林鄉一帶；K的分布較為均勻，空間變異程度較低。Se元素高值區主要分布于長發鎮、向榮鄉和愛民鄉；Ge元素分布狀態較為均勻，在海北鎮稍高；陽離子交換量的分布形態與有機質類似。

圖3 部分指標含量反演分布

重金屬元素As、Cr、Cd、Hg、Pb的含量均低于規范[22]規定的風險管控值，而在空間分布的相對趨勢上，Hg元素的分布較為均勻，Pb及其他三種元素的相對高值區均位于研究區南部。葉綠素和葉面積指數的分布規律較為接近，整體上分布狀態較為均勻，高值區相對集中于長發鎮和東林鄉。

2.4 土壤質量綜合評價結果分析

評價范圍為研究區的旱田范圍，將所研究的地塊按照等級進行分級顯示，其余區域(水田、居民地等)以遙感影像作為底圖，最終評價結果如圖4所示。結果統計顯示，研究區地塊等級均在二等及以上，可見該區土壤質量整體較高。其中一等(優質)地塊面積占全區的58.63%，二等(良好)地塊面積占全區的41.37%。在空間分布上，兩個等級地塊分布較為均勻，呈交錯分布狀態。

圖4 研究區土地質量綜合評價結果

通過單指標評價結果進一步分析土壤綜合等級的影響因素。根據單指標評分，SOM、N、K和CEC四種養分指標評級均以一等為主，其中SOM一等地塊占比98.23%，N一等地塊占比96.34%，P一等地塊占比100%，CEC一等地塊占比98.65%；P和Se以二等地塊為主，分別占比95.31%和97.45%。

五種重金屬元素及全鹽量的評級中一等地塊均占比100%，表明研究區土壤重金屬污染及鹽堿化的威脅較小，土壤環境質量較高。農作物數據由于是單期影像，因此主要評價空間上不同地塊之間的長勢差異性，評價結果顯示葉綠素一等地塊面積占比53.26%，二等地塊占比46.74%，葉面積指數一等地塊面積占比46.35%，二等地塊占比53.65%。

根據地面數據驗證航空高光譜評價結果的準確性。由于評價單元達到了地塊級，所以挑選與評價單元相匹配的最大比例尺的地面數據，地面驗證數據比例尺達到1∶10 000，分布區域主要位于長發鎮，面積占研究區總面積的21.3%，該部分地面數據并未參與前期的建模及土地質量評價，因此可以保證驗證的客觀性。地面數據的綜合評價方法與航空數據評價方法完全一致，以地塊為單元對比兩類數據評價等級的一致性，二者一致率達到97.6%，表明航空高光譜評價結果可靠。

3 討論

本研究基于航空高光譜技術構建土壤質量綜合評價模型，將航空高光譜與地球化學調查進行有機結合，發揮了航空高光譜技術空間及光譜分辨率高、時效性強等優勢。土地質量地球化學評價工作在大尺度(地塊級)調查中的比例尺通常需達到1∶10 000及以上，而本次航空高光譜數據空間分辨率為3.75 m，具備數據“全覆蓋”的特征，數據精度有較大優勢。土壤質量評價中單純依靠地面調查工作將會耗費大量的人力、物力，工期較長，而航空高光譜調查能夠有效縮減工作周期，一次飛行即可獲得全區數據，能夠有效減弱由于不同數據獲取周期引起的系統誤差。同時，所處研究區此前開展的地面調查數據能夠用于綜合研究，提升數據的利用率。

數據反演質量是評價的基礎，目前常用的反演建模方法[31-37]包括多元逐步線性回歸、支持向量機、偏最小二乘法、人工神經網絡、隨機森林等多種方法。本研究選擇建模方法的原則是方法應用效果較為成熟并且樣本需求量較少，以便于后續土壤航空高光譜調查工作的工程性推廣，最終選用了偏最小二乘法。偏最小二乘法結合了建模類型的預測分析方法和非模型式的數據內側分析方法，在一個算法下同時進行數據結構簡化(主成分分析)、回歸建模(多元線性回歸)以及兩組變量之間的相關性分析(單相關分析)，簡化了多維復雜數據的分析難度。此次評價指標中有機質、養分和農作物指標的建模總體上獲得了較好的反演效果，各項指標的驗證精度與建模精度對比沒有明顯的下降，表明建模沒有出現明顯的過擬合現象，為評價工作提供了較可靠的數據基礎。

利用航空高光譜反演獲取了土壤肥力、環境與農作物長勢等多項土壤質量指標，基于層次分析法構建了土壤質量綜合評價模型。綜合評價模型考慮了影響土壤質量的多方面因素，在為地方農業生產及規劃提供決策依據時兼顧了全面性與簡潔性，使整個調查系統更為完善。

目前基于航空高光譜技術開展土壤質量的綜合評價尚處于起步階段，后續工作可對以下方面進行改進和完善：①土地質量參數反演方法。不同質量參數的光譜特征不盡相同，機理較為復雜，可借助深度學習方法充分研究光譜機理，進一步提升反演精度；②評價指標補充完善。本研究篩選了一系列指標進行綜合建模，但還需要充分挖掘航空高光譜對土壤質量參數的反演能力，擴充及優化評價指標和模型結構，提升評價的綜合性和客觀性。