楊 曉 宇，包 妮 沙，曹 粵，劉 善 軍

(1.東北大學資源與土木工程學院，遼寧 沈陽 110819；2.煤炭科學研究總院礦山大數(shù)據(jù)研究院/煤炭資源高效開采與潔凈利用國家重點實驗室，北京 100013)

0 引言

土壤有機質(zhì)(Soil Organic Matter,SOM)及土壤全氮(Soil Total Nitrogen,STN)含量是土壤質(zhì)量和肥力的重要指標，也是精準農(nóng)業(yè)監(jiān)測的重要部分，其空間分布可為作物種植者提供重要的參考信息，從而提高田間管理水平[1]。傳統(tǒng)的土壤養(yǎng)分含量主要通過野外采樣和實驗室化學分析獲取，雖然精度較高，但費時費力且會對土壤造成破壞，同時監(jiān)測結(jié)果具有滯后性，無法滿足土壤質(zhì)量動態(tài)監(jiān)測的需求。高光譜遙感為野外土壤信息的快速獲取提供了新方式，也成為傳統(tǒng)土壤理化測試分析向土壤野外實時監(jiān)測方向發(fā)展的重要技術(shù)支撐[2]。

為探討如何選擇有效特征波段和反演模型提高土壤養(yǎng)分的反演和制圖精度問題，本文以松遼平原秋收后的農(nóng)田土壤為研究對象，針對無人機高光譜成像儀獲取的圖像，對比選擇最優(yōu)無人機圖像光譜去噪方法，引入競爭自適應(yīng)重加權(quán)采樣(Competitive Adaptive Reweighted Sampling,CARS)和連續(xù)投影方法(Successive Projections Algorithm,SPA)，篩選特征波段并消除波段間的共線性，明確土壤有機質(zhì)含量和全氮含量最佳光譜波段，進一步引入PSO-ELM并與PLSR和SVM對比，建立針對東北農(nóng)田土壤有機質(zhì)和全氮的高精度估測模型，分析無人機成像光譜技術(shù)農(nóng)田土壤養(yǎng)分估測能力，為田間管理提供技術(shù)指導。

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

研究區(qū)位于遼寧省沈陽市法庫縣農(nóng)田，遙感影像及采樣點分布如圖1所示。法庫縣地處松遼平原，地勢北高南低，海拔40～443 m，遼河及其支流秀水河、拉馬河貫穿全境，在東、南、西及中部形成較大面積的沖積平原，全縣地貌特征為“三山一水六分田”，耕地資源豐富。該區(qū)域?qū)儆诒睖貛Т箨懶约撅L氣候，年均氣溫6.7 ℃，年均降水量600 mm，主要土壤類型是棕壤和草甸土，富含豐富的有機質(zhì)和氮，農(nóng)作物包括玉米和花生。

圖1 研究區(qū)遙感影像及采樣點空間分布Fig.1 RS image of the study area and spatial distribution of sampling points

1.2 研究數(shù)據(jù)

(1)無人機高光譜圖像數(shù)據(jù)。圖像拍攝時間為2019年11月15日，天氣晴朗無云，前5天無降雪，地表無植被和積雪覆蓋，也無秸稈殘留，土壤含水率為2%，土壤顆粒大小均勻，對土壤光譜的可見光、近紅外波段影響極小[21,22]。采用大疆M600 Pro無人機搭載Resonon Pika L 高光譜成像儀，覆蓋可見光—近紅外區(qū)域中的400～1 000 nm共150個波段，光譜分辨率為2.1 nm，無人機飛行高度為100 m，圖像的空間分辨率為0.1 m。對原始高光譜圖像進行輻射校正、幾何校正預處理，并通過Sigma濾波消除田壟陰影等影響，最后剔除噪聲波段，保留403～759 nm和783～900 nm共112個波段作為原始反射率圖像。

(2)土壤有機質(zhì)和全氮含量數(shù)據(jù)。采用分區(qū)隨機劃分方法，在1 m×1 m的樣方內(nèi)采集0～10 cm深度處的土壤，每一樣點采用四分法取樣，共采集68個土樣，并用手持GPS記錄其坐標。對土樣進行風干、研磨及過篩(0.25 mm)處理，進而通過重鉻酸鉀—硫酸溶液氧化法測定其有機質(zhì)(SOM)含量，通過凱式法測定其全氮(STN)含量。分別對SOM和STN含量進行統(tǒng)計分析(表1)：SOM含量介于11.4～30.6 g/kg之間，平均值為19.72 g/kg，變異系數(shù)為24%；STN含量介于0.84～2.08 g/kg之間，平均值為1.38 g/kg，變異系數(shù)為22%；SOM和STN在全樣本集、建模集、驗證集上的均值和標準差均較接近，且均服從正態(tài)分布。依據(jù)SOM和STN含量，通過分層采樣法將68個土壤樣本按4∶1的比例[23]分為55個訓練樣本、13個驗證樣本，并考慮空間代表性對部分樣本進行微調(diào)(圖 1b)。

表1 采樣點有機質(zhì)含量及全氮含量統(tǒng)計Table 1 Statistics of SOM and STN content for sampling points

2 研究方法

2.1 光譜去噪方法

光譜預處理可有效去除光譜噪聲，應(yīng)用最廣泛的方法包括散射校正和光譜微分[24]。散射校正主要有標準正態(tài)變量變換(Standard Normal Variate transformation，SNV)和多元散射校正(Multiple Scattering Correction，MSC)：SNV可以削弱土壤顆粒大小不均勻以及粒子表面非特異性散射的影響，其前提是假定光譜每個波段的反射率均滿足特定的分布特征(如正態(tài)分布)，由原始光譜減去該條光譜的平均值，再除以該條光譜的標準差計算；MSC可有效消除土壤顆粒間的散射效應(yīng)，增強光譜數(shù)據(jù)中與成分含量相關(guān)的光譜信息，其首先將所有樣本的平均光譜作為標準光譜，對每個樣本光譜與標準光譜進行一元線性回歸，然后由每個樣本光譜減去線性平移量，再除以傾斜偏移量計算。光譜微分常用光譜一階微分(First Derivate,FD)和二階微分(Second Derivate,SD)，其不僅可消除光譜中的基線平移、大氣散射及其他背景的干擾，還可放大光譜曲線斜率的細微變化，提高分辨率和靈敏度。

2.2 特征篩選及建模方法

競爭自適應(yīng)重加權(quán)采樣(Competitive Adaptive Reweighted Sampling,CARS)是一種基于蒙特卡洛采樣和偏最小二乘回歸(PLSR)的特征波段選擇方法[25]，通過有監(jiān)督的特征選擇，可選出對預測變量敏感的波段子集[26]，但仍可能存在波段間相關(guān)性高的特征冗余問題[27]。連續(xù)投影算法(Successive Projections Algorithm,SPA)是一種使矢量空間共線性最小化的特征選擇算法[28]，本文在CARS的基礎(chǔ)上進一步用SPA消除波段間的共線性[29]，可實現(xiàn)特征波段集包含的信息量最多且內(nèi)部相似性最小[30-33]。

PLSR將原始自變量空間中的信息通過主成分分析投影到部分潛在變量中，進而基于少量的潛在變量實現(xiàn)線性回歸。SVM是基于核函數(shù)的統(tǒng)計機器學習方法，通過設(shè)計非線性核函數(shù)將原始的低維數(shù)據(jù)映射到高維特征空間，進而利用SVM對高維數(shù)據(jù)進行線性回歸。 ELM隨機選取輸入層權(quán)重和閾值，輸出層權(quán)重通過最小化損失函數(shù)依據(jù)Moore-Penrose(MP)廣義逆矩陣理論進行計算求解，其損失函數(shù)由訓練誤差項和輸出層權(quán)重范數(shù)的正則項構(gòu)成[34]。粒子群(PSO)算法是一種模擬鳥類覓食行為的群智能算法[35]，利用PSO算法優(yōu)化選擇極限學習機輸入層權(quán)重和閾值，可減少ELM所需的隱含層節(jié)點數(shù)量，提高模型泛化能力；該算法通過比較自身歷史最優(yōu)位置和種群最優(yōu)位置更新當前位置和速度，使粒子不斷向最優(yōu)解逼近[36]。粒子群優(yōu)化的極限學習機算法(PSO-ELM)實現(xiàn)過程(圖2)為：首先對PSO參數(shù)進行初始化，慣性權(quán)值設(shè)置為[0.4,0.8][37,38]，最大迭代次數(shù)設(shè)置為100，學習因子分別設(shè)置為2.4 和1.6[39,40]；然后將每個粒子對應(yīng)的輸入權(quán)值和閾值代入ELM模型，將預測值與實測值的均方誤差(MSE)作為PSO的適應(yīng)度，計算個體極值和全局極值；通過迭代比較更新粒子位置和速度，計算粒子適應(yīng)值，更新粒子極值和全局極值，直到取得最小誤差或達到最大迭代次數(shù)，輸出最優(yōu)輸入權(quán)值和閾值，代入ELM模型進行預測。

圖2 PSO-ELM算法流程Fig.2 Flowchart of PSO-ELM

將篩選的高光譜特征波段作為模型輸入數(shù)據(jù)，分別通過PLSR、SVM、ELM和PSO-ELM訓練并反演土壤有機質(zhì)和全氮含量，以驗證集樣本反演結(jié)果的決定系數(shù)(R2)、平均相對誤差(MAPE)、相對分析誤差(RPD)為指標評價模型精度?，F(xiàn)有研究一般認為，RPD>2的模型可實現(xiàn)準確預測，RPD介于1.4～2之間的模型預測能力中等，RPD<1.4的模型幾乎沒有預測能力[41]。

3 結(jié)果與分析

3.1 光譜去噪效果分析

對原始反射率圖像分別進行多元散射校正(MSC)、標準正態(tài)變量變換(SNV)、光譜一階微分(FD)、光譜二階微分(SD)。由圖3可知，多元散射校正和標準正態(tài)變量變換能保留地物特征信息(圖 3b、圖 3c)，而光譜微分處理后的圖像光譜信息基本丟失，無法有效區(qū)分地物的紋理和平滑區(qū)域等典型特征(圖 3d、圖 3e)。通過圖像信噪比定量評價不同方法的去噪效果(圖4)，可見多元散射校正方法和標準正態(tài)變量變換方法提高了原始反射率圖像的信噪比，前者在514 nm處將信噪比從31 dB 提高到160 dB，后者在715 nm處將信噪比從33 dB 提高到99 dB，而光譜微分方法降低了信噪比。由原始光譜、不同方法去噪后光譜反射率與SOM和STN的相關(guān)性分析結(jié)果(圖 5)可知，4種光譜去噪方法均提高了光譜反射率與土壤屬性的相關(guān)性，其中多元散射校正方法提高的幅度最大，將光譜反射率與SOM的最大相關(guān)系數(shù)絕對值從0.24提高到0.54，將光譜反射率與STN的最大相關(guān)系數(shù)絕對值從0.28提高到0.56。

圖3 不同光譜去噪方法的圖像去噪結(jié)果Fig.3 UAV hyperspectral images denoised by different spectral denoising methods

圖4 不同光譜去噪方法下的圖像信噪比Fig.4 Signal-to-noise ratio of images denoised by different spectral denoising methods

圖5 不同光譜去噪方法下樣本光譜反射率與土壤屬性之間的最大皮爾遜相關(guān)系數(shù)絕對值|r|Fig.5 Absolute value of the maximum Pearson′s correlation coefficient (|r|) between sample spectral reflectance and soil properties based on different spectral denoising methods

綜上，本文嘗試的4種光譜去噪方法中，MSC的去噪效果最好，能在增強光譜反射率與土壤屬性相關(guān)性基礎(chǔ)上，最大限度抑制影像噪聲，提高土壤制圖精度。Casa等[41]的研究也指出，MSC可降低原始光譜的噪聲，增強光譜反射率與土壤性質(zhì)的關(guān)系，并為土壤定量分析引入合適的預測因子。因此，本文選擇MSC作為最佳光譜去噪方法，并基于MSC處理后的光譜數(shù)據(jù)進行后續(xù)特征選擇和建模。

3.2 光譜特征選擇

圖6 不同方法篩選的特征波段Fig.6 Feature bands selected by different methods

3.3 模型精度評價

基于篩選的高光譜特征波段，分別采用PLSR、SVM、ELM和PSO-ELM對SOM和STN建立反演模型，驗證集精度如表2所示。分析結(jié)果表明，基于全波段的PLSR模型預測SOM的R2=0.55、MAPE=15.4%、RPD=1.31，預測STN的R2=0.54、MAPE=16.2%、RPD=1.44，而基于全波段的SVM、ELM及PSO-ELM模型預測SOM和STN的精度均不理想。由于PLSR模型通過主成分分析嵌入了特征降維和篩選過程，因此，本文對CARS和CARS-SPA特征波段集不進行PLSR建模。CARS去除了全光譜中的無關(guān)變量，提高了SVM、ELM和PSO-ELM模型的R2值，但RPD值仍小于或稍大于1.4，預測能力較差。采用CARS-SPA方法二次篩選特征波段，SVM、ELM及PSO-ELM模型預測精度均進一步提高。ELM模型的精度低于SVM模型，但引入粒子群優(yōu)化算法后有效提高了ELM模型的精度，PSO-ELM模型預測SOM的R2=0.73、MAPE=12.6%、RPD=1.91，預測STN的R2=0.63、MAPE=12.6%、RPD=1.53，預測精度總體上高于SVM和PLSR模型。

表2 不同特征選擇方法的模型反演SOM、STN的精度Table 2 Inversion accuracies of SOM and STN with different feature selection methods and different models

STN的預測精度整體低于SOM，主要原因在于：首先，土壤氮素含量變化對光譜的影響微弱，含氮官能基在吸收帶中多以二階以上倍頻或合頻存在，吸光系數(shù)偏小，宋雪等[48]在構(gòu)建太湖地區(qū)土壤全氮高光譜反演模型的研究中，人工神經(jīng)網(wǎng)絡(luò)模型反演結(jié)果的精度僅為R2=0.69、RPD=1.46；其次，本文采集的土壤樣本全氮含量較低，平均值為1.38 g/kg，因此對光譜的響應(yīng)也較低，已有研究中，一般STN含量大于5 g/kg的反演結(jié)果能達到RPD>2的預測精度[49,50]。因此，高光譜用于低含量的土壤全氮估算時，精度會受到一定限制。

3.4 土壤養(yǎng)分含量空間分布特征

利用基于CARS-SPA算法篩選的特征波段集構(gòu)建的PSO-ELM模型對研究區(qū)逐像元反演SOM和STN。如圖7所示，SOM和STN值分別介于12～30 g/kg和0.8～2 g/kg之間，SOM和STN的高值區(qū)主要位于東北部的玉米種植區(qū)。按照農(nóng)田土壤養(yǎng)分分級標準[51]，分別統(tǒng)計玉米和花生種植區(qū)各養(yǎng)分等級的面積占比(圖8)，可知玉米種植區(qū)的有機質(zhì)含量豐缺度為中等而花生種植區(qū)為稍缺，兩種作物種植區(qū)的全氮含量豐缺度均為中等。

圖7 基于最優(yōu)模型反演的土壤有機質(zhì)及全氮含量空間分布Fig.7 Spatial distribution of SOM and STN content based on optimal regression models

圖8 兩種作物種植區(qū)土壤有機質(zhì)及全氮含量等級分布Fig.8 Grade distribution of SOM and STN content in two crop planting areas

4 結(jié)論

本文基于無人機高光譜圖像數(shù)據(jù)，通過改進的極限學習機算法建立農(nóng)田土壤有機質(zhì)和全氮含量高光譜反演模型，結(jié)論如下：1)4種光譜去噪方法的對比表明，無人機高光譜圖像噪聲主要來自土壤間散射效應(yīng)，相比原始光譜反射率，多元散射校正(MSC)能降低原始光譜噪聲，增強光譜反射率與土壤屬性間的相關(guān)性，將圖像信噪比從31 dB 提高到160 dB，將光譜反射率與土壤有機質(zhì)(SOM)的最大相關(guān)系數(shù)絕對值從0.24提高到0.54，將光譜反射率與土壤全氮(STN)的最大相關(guān)系數(shù)絕對值從0.28提高到0.56。2)通過競爭自適應(yīng)重加權(quán)采樣(CARS)算法及其與連續(xù)投影(SPA)組合(CARS-SPA)算法對全波段進行光譜特征選擇，篩選后的特征波段較全波段均具有更強的預測能力，其中CARS-SPA算法能有效去除全波段中與SOM和STN的無關(guān)變量，且可進一步去除波段間共線性高的冗余變量，將SOM和STN的平均相對誤差分別降低了12.3%和6.6%。3)通過粒子群算法(PSO)對極限學習機(ELM)模型的輸入層權(quán)重和閾值進行優(yōu)化，有效彌補了ELM模型泛化性能低這一不足，將SOM和STN的平均相對誤差分別降低了2.9%和3.2%，同時模型精度高于偏最小二乘回歸(PLSR)和支持向量機(SVM)。PSO-ELM預測SOM的決定系數(shù)R2為0.73，相對分析誤差RPD為1.91，預測STN的決定系數(shù)R2為0.63，相對分析誤差RPD為1.53，利用該最優(yōu)模型繪制研究區(qū)域的SOM和STN空間分布圖，二者高值區(qū)主要位于東北部的玉米種植區(qū)，表明無人機高光譜圖像通過合適的光譜去噪方法和光譜特征選擇可應(yīng)用于田間尺度的土壤養(yǎng)分預測。

本研究使用的是400～1 000 nm成像高光譜數(shù)據(jù)，隨著成像高光譜傳感器發(fā)展和普及，1 000～2 500 nm波段數(shù)據(jù)的增加將進一步提高有機質(zhì)等土壤養(yǎng)分的預測精度，下一步將考慮地形及土壤粗糙度對光譜的影響，并在大尺度區(qū)域范圍內(nèi)和更多的田間樣本基礎(chǔ)上驗證算法的普適性。