基于高光譜的黑土區(qū)土壤重金屬含量估測

林 楠 劉翰霖 孟祥發(fā) 劉海琪 楊佳佳

(1.吉林建筑大學測繪與勘查工程學院，長春 130018； 2.東北大學資源與土木工程學院，沈陽 110819；3.中國地質調查局沈陽地質調查中心，沈陽 110034)

0 引言

東北黑土區(qū)土壤肥沃、養(yǎng)分供應能力強，是我國寶貴的耕地資源。黑土中銅(Cu)、鋅(Zn)、錳(Mn)等重金屬元素在作物生長過程中起著至關重要的作用，是土壤質量監(jiān)測的重要指標。重金屬元素過量積累會影響土壤微生物多樣性，從而破壞土壤肥力和緩沖能力，造成土壤質量下降。因此，快速、準確地估測黑土中的重金屬含量對黑土區(qū)科學管理和持續(xù)保護具有重要意義[1-2]。近年來，高光譜遙感以其快速、準確的優(yōu)點成功應用于土壤元素含量的估測計算中，并取得了一定成果[3-9]。

高光譜遙感通過連續(xù)緊密的波段記錄土壤光譜信息，其數(shù)據(jù)量大、信息豐富。由于波段數(shù)量較多，使波段之間存在信息冗余現(xiàn)象，從而影響了建模預測結果[10-11]。目前，大多數(shù)研究都是首先計算光譜波段與元素含量的相關系數(shù)，然后通過設定相關系數(shù)的閾值來選擇敏感波段。陳元鵬等[12]以工礦復墾區(qū)為試驗區(qū)域，通過計算相關系數(shù)進行實測高光譜數(shù)據(jù)的優(yōu)選，結合實測的土壤重金屬含量，采用偏最小二乘、隨機森林、支持向量機等回歸分析模型進行土壤重金屬含量反演試驗研究。王金鳳等[13]以喀斯特流域為研究區(qū)，基于測定土壤樣品的鋅元素含量和采集的土壤光譜數(shù)據(jù)，利用相關分析進一步篩選特征變量，并采用非線性算法揭示光譜敏感波段反射率與重金鋅元素含量之間的映射關系，從而進行土壤重金屬含量估測。由于高光譜測量數(shù)據(jù)的高維非線性特征，通過設置相關系數(shù)閾值提取后的敏感波段仍然可能存在信息冗余，進而影響模型估測精度。核主成分分析(Kernel principal component analysis, KPCA)是一種有效的特征識別和提取方法，它結合了核函數(shù)與主成分分析(Principal component analysis，PCA)的特點，具有極強的非線性特征提取能力，常用于高維數(shù)據(jù)分類和異常識別[14]。本研究對黑土中銅(Cu)、鋅(Zn)、錳(Mn)重金屬元素的光譜反射率及其特征進行分析，分別與各元素含量進行相關性計算，引入KPCA方法，進行敏感特征波段選取，與極限學習機(Extreme learning machine，ELM)模型相結合構建元素含量預測模型，并對KPCA降維前后模型的估測精度進行對比分析，以期為黑土區(qū)重金屬含量高光譜估測研究提供一種高效的方法。

1 材料與方法

1.1 研究區(qū)概況

土壤樣品采集試驗區(qū)位于黑龍江省訥河市，該區(qū)地處松遼平原北端，大小興安嶺南緣，是黑土分布的典型區(qū)域。區(qū)內耕地面積廣闊，農(nóng)產(chǎn)品種類豐富，是甜菜和馬鈴薯的盛產(chǎn)地。2019年4月在試驗區(qū)內挑選典型黑土地塊，共采集土壤樣本80個(圖1)，采集樣品過程中確保采樣點距離公路至少150 m，以每個采樣點位置為中心，在其周圍5 m×5 m范圍內進行樣品采集，每個采樣點位置共采集5份黑土樣本，取樣深度為15 cm，將樣本充分混合后裝入采樣袋。

1.2 元素含量測量及樣本劃分

將樣品中的秸稈、砂礫等雜物剔除，風干后研磨過篩，使得土壤粒徑小于0.25 mm，將合格的樣品分成2份，1份用于測定土壤元素含量，另1份用于室內高光譜測量[15]。結合多目標樣品分析質量要求等技術規(guī)范，在滿足規(guī)范要求樣品分析的檢出限、準確度、精密度等技術條件的前提下，參照《土壤農(nóng)業(yè)化學分析方法》，選用X射線熒光光譜法進行土壤銅(Cu)、鋅(Zn)、錳(Mn)3種重金屬元素含量的測定。對元素含量的測量分析結果進行統(tǒng)計(表1)，將80份土壤樣品按元素含量由低到高分為20組，每組隨機抽取1個樣品放入驗證集，共20個樣品作為驗證集，其余60個樣品為訓練集。

表1 土壤樣本重金屬含量信息Tab.1 Statistics of soil samples heavy metal contents

1.3 光譜采集及特征變化

在暗室中進行土壤光譜測量，本次試驗測量儀器選用ASD FiledSpec4型地物波譜儀，為提高光譜測量數(shù)據(jù)的精度，測量中取5次光譜測量的均值作為土壤樣品的反射光譜數(shù)據(jù)。由于光譜分辨率較高，波段數(shù)較多，相鄰波段之間可能存在信息重疊，更容易受到噪聲的影響，因此，將光譜數(shù)據(jù)進行重采樣和去噪處理，設置采樣間隔為10 nm。在去噪和重采樣處理的基礎上，對原始光譜反射率(R)進行光譜一階微分交換(FD)、光譜連續(xù)統(tǒng)去除變換(CR)、光譜連續(xù)統(tǒng)去除一階微分交換(CR-FD)，能夠有效突出光譜曲線峰谷，可以快速準確地確定敏感波段(圖2)。

從圖2可以看出，采集的黑土樣品原始光譜反射率在0～0.9之間，各樣品光譜曲線波動形狀相似，在可見光波段，隨著波長的增加，反射率逐漸增加，直到1 200 nm處反射率趨于穩(wěn)定。近紅外區(qū)土壤光譜反射率總體高于可見光區(qū)，兩個明顯的波谷分布在1 400 nm和1 900 nm附近，主要是土壤中殘留的水分和空氣中水蒸氣吸收造成的，而在2 200 nm處有輕微的凹陷，主要是受土壤中存在粘土礦物的影響。圖2b～2d為原始光譜經(jīng)過3種不同變換后的光譜曲線，可以看出3種變換均可以對原始光譜特征起到放大作用，經(jīng)變換后反射率在1 400、1 900、2 200 nm附近波動更大。

2 模型原理

2.1 極限學習機

極限學習機(Extreme learning machine, ELM)是單隱層前饋神經(jīng)網(wǎng)絡(Single-hidden layer feedforward neural network, SLFN)學習訓練的新模型，相比于傳統(tǒng)的人工神經(jīng)網(wǎng)絡，ELM算法的初始化參數(shù)少、學習訓練速度快，建模過程不需要調整隱含層參數(shù)，相當于一個線性參數(shù)模式，可以簡單地求出系統(tǒng)的最小二乘解[16-17]。極限學習機的輸出函數(shù)定義為[18-20]

(1)

其中

β=(β1,β2,…,βL)

(2)

h(x)=(h1(x),h2(x),…,hL(x))T

(3)

h(wixj+bi)=G(w,b,x)

(4)

式中wi——隱含層權重

bi——隱含層偏置參數(shù)

β——L個節(jié)點的隱含層到m個輸出節(jié)點之間的輸出權向量

h(x)——極限學習機非線性特征映射關系

h(wixj+bi)——第i個隱藏節(jié)點的輸出函數(shù)

G(w,b,x)——連續(xù)函數(shù)的非線性映射，即激活函數(shù)

β可以通過最小二乘解獲得，即

(5)

式中H——隱含層輸出矩陣(隨機矩陣)

T——訓練數(shù)據(jù)目標矩陣

最小二乘解為

=H+T

(6)

式中H+——H的Moore-Penrose廣義逆

通常采用奇異值分解法(Singular value decomposition, SVD)計算矩陣的Moore-Penrose廣義逆。由于ELM算法隨機定義訓練數(shù)據(jù)的數(shù)量，并且每次計算只更新權重而不重新訓練樣本，使得ELM算法在訓練時間和精度上表現(xiàn)十分優(yōu)異。

2.2 核主成分分析

核函數(shù)是處理非線性多維數(shù)據(jù)統(tǒng)計的有效方法，核主成分分析是主成分分析改進，基本思想是將核函數(shù)與PCA相結合，有效減少非線性數(shù)據(jù)的維數(shù)。KPCA的核心思想是將輸入的非線性數(shù)據(jù)轉換到高維特征空間F，進行主成分分析，特征空間F通常為高維歐氏空間或無窮維Hilbert空間[21-23]。

令x1,x2，…，xn∈Rd為KPCA學習的n個訓練樣本，得到第i個經(jīng)KPCA變換后的特征向量為

(7)

其中

(8)

式中k(xi,x)——特征空間F中任意兩個變量經(jīng)核函數(shù)映射后的內積

γi——核矩陣K中特征值對應的標準正交化特征向量

3 基于KPCA-ELM模型的土壤重金屬含量估測

3.1 相關性分析

分別計算了土壤銅(Cu)、鋅(Zn)、錳(Mn)含量與土壤原始光譜反射率及其變換的相關系數(shù)，并繪制了相關系數(shù)曲線(圖3)。從圖3a可以看出，土壤鋅(Zn)和錳(Mn)元素含量與原始光譜反射率的相關系數(shù)均為正值，而銅(Cu)元素含量與原始光譜反射率的相關系數(shù)則正負均有。與原始光譜反射率相比，經(jīng)特征變化后的光譜數(shù)據(jù)與土壤重金屬元素含量的相關性更高，且特征變化后光譜數(shù)據(jù)與重金屬元素含量的相關系數(shù)表現(xiàn)為正、負交叉，同時波峰和波谷的個數(shù)增多，峰值較多，各元素的最高相關系數(shù)顯著提高。

試驗中對計算的相關系數(shù)進行統(tǒng)計，并設定相關系數(shù)閾值進行敏感波段提取，將提取的敏感波段作為預測模型的樣本輸入數(shù)據(jù)(表2)。由表2可以看出，經(jīng)過不同的光譜特征變換，元素含量與光譜反射率的相關系數(shù)計算結果差別較大，其中，土壤銅(Cu)元素的最佳變換形式為連續(xù)統(tǒng)去除變換，與光譜反射率的相關系數(shù)最高為0.527，波長在2 310 nm附近；土壤鋅(Zn)元素的最佳變換形式為連續(xù)統(tǒng)去除一階微分變換，相關系數(shù)最高為0.657，波長在440 nm附近；土壤錳(Mn)元素對應的最佳變換形式為連續(xù)統(tǒng)去除變換，相關系數(shù)最高為-0.513，波長在2 230 nm附近。

表2 土壤元素最大相關系數(shù)和敏感波段Tab.2 Maximum correlation coefficients and sensitive bands

3.2 核主成分特征波段選擇

從相關系數(shù)計算結果可以看出(表2)，設定提取閾值后，3種重金屬元素敏感波段的數(shù)量依然較多，光譜波段之間仍然存在數(shù)據(jù)冗余，試驗采用KPCA方法對光譜敏感波段進行特征提取，得到相互獨立的特征變量用于建模預測。利用KPCA方法提取特征敏感信息具體步驟如下：①將選取60個訓練樣品的最佳光譜變換形式對應的波段作為樣本輸入數(shù)據(jù)，試驗中土壤銅(Cu)元素輸入波段數(shù)為44個，土壤鋅(Zn)元素輸入波段數(shù)為61個，土壤錳(Mn)元素輸入波段數(shù)為51個，對各元素輸入波段數(shù)據(jù)進行歸一化處理。②選擇KPCA的核函數(shù)，計算核矩陣K，本次試驗核函數(shù)選擇徑向基核函數(shù)。③求核矩陣K的特征值和特征向量。④計算特征值的累計貢獻率Ci，設定閾值K，當累計貢獻率大于閾值時提取對應的前i個主成分，累計貢獻率閾值K設定為90%。

根據(jù)以上步驟進行土壤重金屬元素光譜特征主成分提取，計算結果如表3所示，從表中可以看出，KPCA方法特征提取效果顯著，銅(Cu)元素用于建模的變量從44個降為17個，鋅(Zn)元素從61個降為24個，錳(Mn)元素從51個降為19個(圖4)。

表3 KPCA特征提取結果Tab.3 Results of KPCA feature extraction

3.3 基于極限學習機的重金屬含量估測模型

ELM模型基于Python軟件平臺實現(xiàn)，分別將KPCA方法提取的特征主成分分量和土壤元素含量作為模型建模的輸入和輸出端，構建訓練樣本數(shù)據(jù)集[24-25]。極限學習機模型初始化參數(shù)包括隱含層的權重w、偏置參數(shù)b、神經(jīng)元個數(shù)N以及輸出類型，試驗選擇Guass函數(shù)作為輸出函數(shù)，權重和偏置參數(shù)初始化參數(shù)敏感性不強，選擇默認缺省參數(shù)進行初始化。隱層神經(jīng)元個數(shù)N是ELM預測模型的重要參數(shù)，一般來說隱層神經(jīng)元個數(shù)N太小會導致模型擬合不充分，N太大模型又可能過擬合，導致隱層神經(jīng)元輸出矩陣為稀疏矩陣，所以需要對其進行優(yōu)化選擇，選擇決定系數(shù)(R2)作為適應度評價標準。通過迭代計算設置不同神經(jīng)元個數(shù)時預測結果的決定系數(shù)(圖5)，選取決定系數(shù)最大時對應的神經(jīng)元個數(shù)為模型建模的輸入值。通過計算得到，銅(Cu)元素預測模型的隱含層神經(jīng)元個數(shù)N為35時，訓練模型決定系數(shù)達到最大值0.744；鋅(Zn)元素預測模型的隱含層神經(jīng)元個數(shù)N為49時，訓練模型決定系數(shù)達到最大值0.938；錳(Mn)元素預測模型的隱含層神經(jīng)元個數(shù)N為40時，訓練模型決定系數(shù)達到最大值0.767。

3.4 估測結果精度分析

表4 元素含量預測模型結果對比Tab.4 Comparison of element content estimation model results

為進一步對比分析3種元素模型擬合效果，利用模型預測值和實測值繪制預測結果(圖6)。從圖6可以看出，利用KPCA進行特征提取后的ELM模型的預測值更緊密地分布在1∶1左右，數(shù)據(jù)擬合能力和穩(wěn)定性均有一定程度的提高。同時可以看出，訓練集的擬合比驗證集效果好，驗證集樣本點分布相對分散，說明極限學習機模型對訓練樣本的可靠性和一致分布性要求較高。

4 結論

(1)對光譜曲線進行不同形式的特征變換可以有效地增強光譜特征、提高光譜反射率與土壤元素含量的相關系數(shù)，光譜特征變換后土壤銅(Cu)、鋅(Zn)、錳(Mn)元素分別在波長2 310、440、2 230 nm處相關系數(shù)最大。

(2)與傳統(tǒng)的閾值法提取敏感波段方法相比，KPCA方法可以有效降低光譜數(shù)據(jù)的冗余性，實現(xiàn)光譜數(shù)據(jù)特征信息的有效提取。將KPCA方法與ELM模型相結合，構建了KPCA-ELM土壤含量估測模型，該組合模型可以有效選擇最佳變量子集合，提高了模型預測的穩(wěn)定性，特征提取后的ELM模型決定系數(shù)增大、均方根誤差變小，估測精度明顯提高，說明KPCA-ELM模型在土壤重金屬含量高光譜估測中具有一定的適用性，是一種高效的土壤元素含量優(yōu)化估測模型。

(3)基于野外采集黑土樣品和實測高光譜數(shù)據(jù)進行土壤重金屬含量估測，適用于小區(qū)域內土壤重金屬含量的定量反演，由于衛(wèi)星遙感影像的光譜信息受到大氣、地形、地面植被、土壤粗糙度、水分等多因素影響，導致利用遙感影像進行大范圍土壤重金屬含量估算的精度低于室內實測高光譜估測精度，所以針對大范圍的土壤重金屬含量遙感估測還需進一步開展相關分析和研究。