基于高光譜反演的復墾區土壤重金屬含量經驗模型優選

陳元鵬 張世文 羅 明 鄖文聚 鞠正山 李少帥

(1.自然資源部國土整治中心， 北京 100035； 2.安徽理工大學地球與環境學院， 淮南 232001)

0 引言

作為地球生命的重要載體、國土資源的基本元素，土壤不僅是農業生產的根基，也是國家糧食安全的重要保證，同時是水資源、生態環境優化的重要基礎[1]。加強國土資源數量、質量、生態“三位一體”監測與管護，對土壤進行快速準確的調查識別至關重要。然而，傳統土壤調查、評價、測量方法成本高、效率低，無法實現大范圍的快速檢測[2]，與之相比，光譜技術快速、便捷，對土地無污染、無破壞，優勢明顯，僅利用單一光譜便可獲取多種土壤理化屬性信息，同時還可以在野外直接進行原位光譜采集[3-4]，節省大量的人財物力，增強了大范圍土壤屬性數據測度的技術可行性，是土壤屬性信息獲取的重要技術手段。而光譜檢測中的高光譜遙感技術更具優勢，其光譜分辨率高、波段連續性強，能獲取連續的地物光譜曲線，使地物信息更加精細地體現在光譜曲線上，實現更準確的反演和監測。隨著技術的不斷進步發展，高光譜遙感目前已被廣泛應用于土壤資源調查、分類、評價、制圖等工作中[5]。

《土地復墾條例》(國務院令第592號)明確規定“復墾土地需要開展五年監測”。但基于傳統方法的長期跟蹤監測需要消耗大量的人財物力，為改變這一現狀迫切需要引入高光譜遙感技術，以更加高效、便捷、成本節約的優勢，對土地復墾項目開展土壤屬性數據的檢測和監測，以新理論、新技術來突破傳統方法的弊端，為土壤調查監測提供一個新的研究和工作方向。現階段，高光譜數據已能夠準確反映地物光譜的細微特征，變換后的光譜在消除背景和噪聲影響，放大提取光譜的吸收和反射特征等方面具有較優效果[6]，對挖掘光譜信息，構建精度高、穩健性好、泛化能力強的模型具有重要作用[7-10]。不少學者利用高光譜遙感數據針對不同區域的As、Cr、Cd、Cu、Fe、Zn、Ni等重金屬含量進行建模反演，獲得了較好的反演精度。但目前利用高光譜數據反演工礦復墾區土壤重金屬含量的研究相對較少，且現有研究中對于反演模型、特征選擇方法的研究也有待進一步深入。

本文基于ASD FieldSpec 4高光譜遙感數據，結合土壤重金屬屬性數據，利用機器學習方法包括線性(偏最小二乘回歸)與非線性(隨機森林與支持向量機)回歸方法，開展基于高光譜數據的土壤重金屬含量反演研究與實驗，采用3種特征選擇方法包括粒子群算法、遺傳算法、Relief F 算法優化經驗模型，以減少模型自變量數量并提高建模精度。

1 研究區概況

研究區為歷史遺留工礦廢棄地復墾項目區，位于四川省瀘州市古藺縣，地理坐標為27°59′43″～28°7′24″N，105°55′41″～106°4′47″E，區域內海拔340～1 207 m，中亞熱帶氣候，年平均氣溫17.1～18.5℃、平均降雨量748.4～1 184.2 mm[11-13]。該項目分為5個復墾區，項目區域東西寬3.4 km，南北長2.8 km，總面積297 hm2。研究區地理區位及遙感影像數據如圖1所示。

2 數據獲取與分析

2.1 土壤樣本數據采集

在研究區域內，綜合考慮復墾方向和復墾措施，采用網格布點、分層抽樣，共采集42個土壤樣品，實際有效采集了表層土壤樣品41個，剖面土壤樣品10個，樣點分布如圖2所示。

土壤樣品采集時間為2017年8月。表層土壤采樣深度為0～25 cm。剖面土壤采集0～25 cm、25～50 cm、50～75 cm 3個層次，采集時間在上茬作物玉米成熟或收獲以后，下茬作物尚未施用底肥和種植以前，該時間段能夠反映采樣地塊的真實養分狀況和供肥能力。野外采樣采用鐵鍬挖采樣坑，挖完坑后，用竹片去除與金屬采樣器接觸的土壤，再采集樣品。

圖1 研究項目區地理位置和遙感影像Fig.1 Location of study area and remote sensing image

圖2 樣點分布圖Fig.2 Map of samples distribution

在采樣點的布設上，以GPS定位點為中心，向四周輻射確定3～5個分樣點，組合成一個混合樣。每個分樣點的采土部位、深度均一致。采樣地塊為長方形時，采用“S”形布設分樣點；采樣地塊近似正方形時，采用“X”形或棋盤形布設分樣點。采樣時均避開了溝渠、林帶、田埂、路邊、舊房基、糞堆及微地形高低不平無代表性地段。采集各分樣點時將土壤掰碎，挑出根系、秸稈、石塊、蟲體等雜物，充分混合后，四分法留取1.0～1.5 kg裝入樣品袋。

土壤樣品在日光下干燥。在樣品干燥過程中采用揉搓法，以免膠結。干燥后的樣品在加工(過篩)前均用木槌輕輕敲打，以便使土壤樣品恢復至自然粒級狀態。

2.2 樣本土壤理化指標分析

所有土壤樣本經風干、研磨并通過2 mm孔徑篩，采用四分法分成2份，分別供化學分析和光譜測試用。采用內標法定量處理樣土數據，結果經儀器校對、回收率校正，土壤理化性質的測試分析方法見表1，表層土壤樣品測定結果見表2、3。

表1 土壤樣品分析方法檢出限Tab.1 Detection methods and limit of soil samples

表2 表層土壤元素含量特征值Tab.2 Eigenvalue of surface soil element contents

表3 項目區土壤本底元素含量特征值Tab.3 Eigenvalue statistics of soil background element content in project area

由表2、3可知： ①復墾區土壤環境指標中重金屬元素鎘(Cd)、鉻(Cr)、砷(As)、鎳(Ni)有明顯的富集，鎘(Cd)、鉻(Cr)、砷(As)、鎳(Ni)含量分別是區域背景值的5.52、2.01、1.82、2.84倍，其均值含量均是區域背景值的1.5倍以上，特別是鎘(Cd)，高達5倍。與區域土壤背景值相比，項目所在地區土壤本底元素鎘(Cd)、鉻(Cr)、砷(As)、鎳(Ni)含量分別是區域背景值的1.89、1.06、0.98、1.48倍，其中砷(As)含量沒有超過區域背景值。②復墾區內監測指標的變異系數均較高，鎘(Cd)、鉻(Cr)、砷(As)、鎳(Ni)變異系數分別是92.39%、33.16%、43.35%、25.23%，特別是鎘(Cd)變異系數最大。這主要是先前采礦活動對環境擾動較大，不確定性因素多，致使鎘(Cd)分布發生突變。而項目區土壤本底元素中也是鎘(Cd)的變異系數最大，為83.33%，而鉻(Cr)、砷(As)、鎳(Ni)變異系數分別是26.38%、45.88%、20.51%。

2.3 樣本土壤光譜數據獲取

土壤光譜數據測量采用美國ASD公司生產的ASD FieldSpec 4地物光譜儀。在室內條件下用標準白板和鹵素光源獲取土壤樣本的反射光譜曲線，光譜波段范圍為350 ～2 500 nm。該光譜儀在不同的波段間有不同的間隔和分辨率，其中350～1 000 nm之間的采樣間隔為1.4 nm，光譜分辨率為3 nm，1 000～2 500 nm范圍內的采樣間隔為2 nm，采樣光譜分辨率為10 nm。經過間隔為1 nm光譜重采樣，兩者的輸出波段數均為2 151。由于光譜很容易受到外界光線的干擾，因此測試選擇在無光的暗室中進行，選擇50 W的鹵素燈為光源、25°裸光纖鏡頭接收反射波段。測量時將樣本放入直徑為100 mm、高2 mm的透明玻璃培養皿中，光源距離樣品50 cm，光線與樣品呈15°，探頭距樣本10 cm位于光源對面，探頭光纖末端位于土壤樣本正上方。為保證測量的準確性，以及最大限度的提高儀器的優化性能，在測量前，先開機預熱30 min，再按照暗電流(DC)采集、儀器優化(OPT)、白板校正(WR)順序對儀器進行調整和校準。為防止測量過程中有其他因素影響光譜信息收集的準確性，每隔10個樣本進行一次白板校正。在光譜采樣過程中由于光譜波段兩側產生不穩定的噪聲區，為此，實驗去除了350～399 nm和2 351～2 500 nm范圍內的波段數據，獲得1 950個波段數據。

2.4 光譜數據預處理方法

光譜特性由土壤性質決定，但土壤是一種復雜的混合物，具有極其復雜的多孔體系，由不同的礦物質、水分、氣體和土壤有機質組成，影響土壤光譜反射特性的主要因素包括土壤有機質、重金屬元素含量、含水率以及土壤質地和母質等。對本次研究采集的41個土壤樣本光譜反射率測定數據進行制圖，利用Savitzky-Golay(S-G)卷積平滑法[14]以10為窗口平滑處理，結果如圖3所示。

圖3 土壤樣點光譜曲線Fig.3 Spectracurves of soil samples

由圖3可知，雖然土壤樣本光譜反射率有所不同，但整體趨勢呈現共同特征，反射率介于0.2～0.8之間，各樣本土壤光譜曲線在波長范圍內呈相似波動性，光譜反射率隨波長增加而增大，在400～780 nm可見光波段范圍內反射率增加較快，在近紅外的780～2 350 nm波段范圍內反射率增加減緩，800～1 400 nm范圍內的近紅外光譜區的光譜反射率高于可見光波段，在1 400～1 500 nm、1 900～2 000 nm和2 200～2 300 nm 3個范圍內存在明顯的水分吸收區間。綜上，光譜數據的波形、波動性、吸收峰等特征與已有的研究結論相符，表明光譜數據的采集情況良好。

為更好地反映土壤重金屬含量與光譜反射率之間的對應關系，消除噪聲干擾，分離重疊樣本，更加準確地尋找特征波段，對光譜曲線進行一階微分變換、對數一階微分變換以及對數倒數的一階微分變換，結果如圖4～6所示。由圖可知，經過一階微分變換后的光譜，有效地消除了基線和背景的干擾，提高了部分波段的分辨率和靈敏度，使重疊樣本得到分離，進一步提取了原始數據中差異性不顯著的光譜信息，使光譜特征波段峰值更加明顯。

圖4 土壤樣點光譜一階微分變換Fig.4 First order differential transformation of soil sample spectra

圖6 土壤樣點光譜對數倒數一階微分變換Fig.6 Logarithmic inverse first order differential transformation of soil sample spectra

圖5 土壤樣點光譜對數一階微分變換Fig.5 Logarithmic first order differential transformation of soil sample spectra

2.5 土壤光譜相關性分析

為進一步分析土壤光譜反射率與重金屬元素之間的關聯，探究土壤重金屬含量與原始、一階微分變換、對數一階微分變換以及對數倒數一階微分變換光譜數據的相關性，運用相關系數法對其進行分析，相關系數計算公式為

(1)

式中ri——土壤重金屬元素含量與高光譜反射率的相關系數

Rni——第n個土壤樣本第i波段的光譜反射率

y——第n個土壤樣本重金屬含量

表4為土壤重金屬元素與光譜最大相關系數絕對值和對應波段。由表4可知，只經過S-G平滑處理的光譜曲線相關系數絕對值在0.16～0.28之間；一階微分變換后的光譜曲線相關系數絕對值在0.36～0.47之間；對數一階微分變換后的光譜曲線相關系數絕對值在0.40～0.55之間；對數倒數一階微分變換后的光譜曲線相關系數絕對值在0.46～0.54之間。土壤重金屬元素含量與光譜曲線的相關性逐步增強，說明土壤光譜經過預處理變換后，均能在一定程度上消除背景因素或系統噪聲等影響，對于提高建模精度能夠起到一定作用。

相關系數絕對值中，不低于0.5的情況分別有：光譜對數一階微分變換下的鎘(Cd)元素相關系數、光譜對數倒數一階微分變換下的鉻(Cr)元素相關系數、光譜對數倒數一階微分變換下的鎳(Ni)元素相關系數。所以本文選擇對數一階微分、對數倒數一階微分變換后的光譜進行建模。

表4 土壤重金屬與光譜最大相關系數絕對值和對應波段Tab.4 Absolute value of maximum correlation coefficient between heavy metal elements and spectra and corresponding band

3 模型選擇與評價

利用41個土壤樣點光譜的對數一階微分、對數倒數一階微分變換后光譜數據進行重金屬的反演建模，分別選擇線性和非線性回歸模型開展建模實驗。

3.1 線性回歸模型

偏最小二乘回歸(Partial least squares regression, PLSR)集成了主成分、典型相關和線性回歸分析方法的優點，其提供一種多對多的線性回歸建模的方法。PLSR是將因子分析和回歸分析結合的方法。

(1)因子分析。將X和Y作如下分解

X=TPT+E

(2)

Y=UQT+F

(3)

式中X——樣本光譜矩陣

Y——樣本待檢測屬性矩陣

T——X的得分矩陣

U——Y的得分矩陣

P——X的載荷(即主成分矩陣)

Q——Y的載荷(即主成分矩陣)

E——用PLSR模型擬合X時所引進的殘差矩陣

F——用PLSR模型擬合Y時所引進的殘差矩陣

(2)回歸分析。將T和U作線性回歸，得

U=TB

(4)

Y1=T1BQ=X1PTBQ

(5)

式中B——關聯系數矩陣

預測中，由未知樣品矩陣X1和校正得到的P，求出未知樣品X1的矩陣T1。

3.2 非線性回歸模型

3.2.1隨機森林回歸

隨機森林(Random forest regression, RFR)是一種基于決策樹的高效機器學習算法，可用于對樣本進行分類，也可用于回歸分析。它屬于非線性分類器，因此可挖掘變量之間復雜的非線性相互依賴關系[15-19]。

3.2.2支持向量機回歸

支持向量機(Support vector machine regression, SVMR)是一種在分類和回歸問題中應用較為廣泛的模型，在分類和回歸問題中均能得到較優效果，支持向量機在回歸分析中，目標是在有限的噪聲樣本基礎上預測出未知的連續函數。支持向量機回歸包括線性回歸和非線性回歸。對于非線性回歸，只需要引入核函數，通過非線性映射將輸入空間映射到高維的特征空間，在高維空間上進行線性回歸。常用的核函數有：徑向基核函數(RBF)、Sigmoid核函數與多項式核函數等。本文選用RBF作為SVM的核函數[20-27]。

3.3 光譜特征選擇方法

基于高光譜數據進行回歸分析，通常自變量較多，過多的樣本數量或存在冗余信息導致建模精度降低，為此需要利用特征選擇算法進行數據降維，減少自變量個數的同時保證建模精度不降低。為進一步降低數據維度、減少冗余信息、提高運算效率、提升模型穩定性與適應性，本文選取3種光譜特征選擇方法來優化模型預測能力，實現模型預測準確性和穩定性的進一步提高[28-30]。

3.3.1粒子群算法

粒子群優化算法(Particle swarm optimization, PSO)模擬的是鳥群覓食行為，是一種基于集群智能的隨機尋優算法，鳥群通過自身經驗和種群之間的交流調整自己的搜尋路徑，從而找到食物最多的地點。PSO算法步驟如下：初始化一個規模為m的粒子群，設定初始位置和速度；計算每個粒子的適應值；對每個粒子，將其適應值和其經歷過的最好位置的適應值進行比較，若較好，則將其作為當前最好位置；分別對粒子的速度和位置進行更新，以達到滿足終止條件[31-33]。PSO算法涉及的參數包括：種群數量、迭代次數、慣性權重、學習因子、空間維數、位置限制和速度限制等。

3.3.2遺傳算法

遺傳算法(Genetic algorithm, GA)是模擬達爾文生物進化論的自然選擇和遺傳學機理的生物進化過程的計算模型，是一種通過模擬自然進化過程搜索最優解的方法[34-36]。GA從一組隨機產生的“種群”(初始可行解)開始搜索過程。其中可行解即被稱為 “染色體”，一個可行解一般由多個元素構成，那么每一個元素就被稱為染色體上的一個“基因”。在算法執行過程中，染色體通常在編碼后進行迭代計算，每次迭代都會生成若干條染色體。采用適應度函數分別計算每一條染色體的適應程度，并根據適應程度計算每一條染色體在下一次進化中被選中的概率，對染色體進行優勝劣汰。在未達到算法終止條件時，通過選擇、交叉和變異3個算子產生下一代染色體。重復此過程直到算法收斂從而得到最優解(最優染色體)。

3.3.3Relief F 算法

Relief F算法是一種考慮了多變量和特征互相依賴作用的過濾式(Filter)特征選擇算法[37]，算法通過“假設間隔”對特征的分類能力進行評價，綜合考慮類間距與類內距，若類間距大于類內距，則增加其權值，若類間距小于類內距，則降低其權值，通過類間距與類內距不斷更新其權值，并根據計算的最終權值進行特征選擇，權值更新公式為

(6)

假設間隔θ是指當保持樣本類別不變時分類決策面可移動的最大距離，即

(7)

式中diff()——不同樣本間的距離

n——樣本數量f——評價特征

mclass(xi)——樣本抽樣次數

i——隨機抽中的樣本

H(x)、M(x)——樣本x的同類和異類最近鄰樣本點

4 結果分析

4.1 線性擬合結果

線性模型選擇偏最小二乘回歸(PLSR)模型，將樣本按照3∶1的比例劃分進行建模驗證(建模樣本30個，驗證樣本11個)，回歸分析結果如表5所示，表中決定系數R2和均方根誤差(RMSE)均為驗證集的評價系數。其中鎘(Cd)元素與對數一階微分變換后的光譜曲線建模擬合效果最好，決定系數R2最大，達0.76，均方根誤差RMSE最小，為0.62 mg/kg，說明基于室內處理的光譜在預測土壤中的重金屬元素鎘(Cd)時能夠達到良好效果；而重金屬元素鎳(Ni)在與經對數一階微分變換的光譜曲線建模擬合中，最優決定系數R2僅為0.29，均方根誤差RMSE為10.6 mg/kg，在與對數倒數一階微分變換的光譜曲線建模擬合中R2則更低，說明在研究區域中，高光譜數據無法滿足精確建模預測土壤重金屬元素鎳(Ni)含量的基本要求。元素砷(As)、鉻(Cr)與對數倒數一階微分變換后的光譜建模，決定系數R2分別為0.54、0.65，但均方根誤差RMSE差別較大，分別為1.94、16.1 mg/kg，雖然決定系數R2基本能夠滿足預測精度，但鉻(Cr)的RMSE為16.1 mg/kg，說明模型穩定性和泛化能力有待進一步增強。

表5 PLSR模型反演結果評價Tab.5 Assessment of PLSR inversion results

4.2 非線性擬合結果

非線性模型選擇隨機森林(RFR)和支持向量機(SVMR)模型，同樣將樣本按照3∶1的比例劃分進行建模驗證(建模樣本30個，驗證樣本11個)，回歸分析結果如表6、7所示。由表6、7可知，利用隨機森林和支持向量機方法的反演結果中，依然是鎘(Cd)元素的決定系數R2最高，鎳(Ni)元素的R2最低，均方根誤差RMSE中鎘(Cd)元素分別為1.69、1.75 mg/kg，鉻(Cr)元素分別為18.20、18.70 mg/kg，結果呈現的整體趨勢與PLSR反演結果較為一致。與PLSR相比較，非線性模型對各重金屬元素擬合精度偏低，究其原因或是因為建模樣本數量過少，導致非線性建模效果不夠理想，在大樣本容量下，預測精度和模型的穩定性將會有所提高。

表6 RFR模型反演結果評價Tab.6 Assessment of RFR inversion results

圖7 基于不同特征選擇的Cd含量反演結果Fig.7 Cd content inversion results based on feature selection

金屬元素最優光譜變換方法R2RMSE/(mg·kg-1)砷(As)對數倒數一階微分0.407.87鎘(Cd)對數一階微分0.481.75鉻(Cr)對數倒數一階微分0.3318.70鎳(Ni)對數一階微分0.2513.65

4.3 基于特征選擇的模型優化

由線性(PLSR)和非線性(RFR、SVMR)模型分析結果可知，在PLSR模型下鎘(Cd)元素與對數一階微分變換后的光譜曲線建模擬合效果最優，精度最高，決定系數R2最大，達到0.76，所以選擇PLSR模型下鎘(Cd)元素的反演作為優化目標，以 PSO、GA、Relief F 3種特征選擇方法對高光譜數據進行特征波段選擇后再進行PLSR建模，以驗證特征選擇、模型優化的效果。

PSO方法參數設置：初始種群數量20、迭代次數20、慣性權重0.33；GA方法參數設置：初始種群數量20，最大繁殖代數20、交叉概率0.6、變異概率0.03。基于以上兩種方法選擇的特征波段數量趨于一致，分別為535、537個。同時基于 Relief F 方法，選擇了535個特征變量。

以PSO、GA、Relief F 3種方法選擇的特征波段為自變量，基于 PLSR 模型對鎘(Cd)元素進行回歸分析，反演結果如表8所示。由表8可知，基于PSO特征選擇的反演結果最優，R2為0.84；其次為基于GA的特征選擇的反演結果，R2為0.64；基于Relief F特征選擇的反演結果R2只有0.58。與未經過特征選擇的反演結果比較，只有基于PSO特征選擇方法的反演精度有所提高，R2由0.76提高到0.84，基于GA、Relief F特征選擇方法的反演精度均出現了一定程度的降低，R2由0.76分別降至0.64和0.58。所以，在3種方法中，基于PSO優化的PLSR模型能夠進一步提高對實驗區內鎘(Cd)元素反演的精度，并有效減少了特征波段的數量。

表8 基于不同特征選擇方法的PLSR模型反演結果評價Tab.8 Assessment of PLSR inversion results based on feature selection

圖7為鎘(Cd)元素反演模型的預測值散點圖與擬合結果。

5 結論

(1)以工礦復墾區為實驗區域，進行了基于高光譜數據的土壤重金屬含量反演研究和經驗模型優選。實驗結果表明，通過對光譜曲線進行一階微分、對數一階微分以及對數倒數的一階微分等數學變換能夠有效提高光譜數據與土壤重金屬含量的相關性；針對實驗區數據，PLSR對復墾區土壤重金屬含量的反演精度最優，尤其對區域內主要障礙因子鎘(Cd)元素含量的反演效果最佳，驗證集的決定系數R2達到0.76；采用的 PSO、GA、Relief F 3種特征選擇方法中，PSO 在降低特征波段變量維度的基礎上，能夠有效提高 PLSR 模型對復墾區土壤重金屬含量反演的精度，使決定系數R2由0.76提至0.84。

(2)運用線性(PLSR)和非線性(RFR、SVMR)機器學習方法開展建模實驗，采用基于隨機搜索策略的3種特征選擇方法(PSO、GA、Relief F)進行模型優化，通過對3種機器學習方法和3種特征選擇方法分析優選，得到基于 PSO-PLSR 的土壤重金屬含量高光譜反演方法，可以獲得較佳建模精度。