基于PCA的土壤Cd含量高光譜反演模型對比研究

郭 飛，許 鎮(zhèn)，馬宏宏，劉秀金，楊 崢，唐世琪

1.中國地質(zhì)科學(xué)院地球物理地球化學(xué)勘查研究所，河北 廊坊 065000 2.中國地質(zhì)調(diào)查局土地質(zhì)量地球化學(xué)調(diào)查評價研究中心，河北 廊坊 065000 3.中國科學(xué)院空天信息創(chuàng)新研究院，北京 100101

引 言

鎘(Cadmium，Cd)是一種有毒重金屬[1]，它不僅會降低土壤微生物的生物活性，還易通過在植物可食部位的累積，進入食物鏈危害人體[2]。礦產(chǎn)資源的開采是造成其周邊農(nóng)用地土壤污染的重要原因之一[3]。如何快速有效測定土壤中Cd的含量及空間分布已成為目前亟待解決的問題。高光譜遙感由于光譜分辨率高、波段連續(xù)，能快速高效獲取精細的光譜信息等優(yōu)勢，成為快速查明土壤重金屬污染狀況的新技術(shù)方法之一[4]。

利用可見-近紅外光譜對土壤重金屬含量進行定量反演已成為國內(nèi)外熱點研究問題。Kemper等[5]利用可見-近紅外光譜，基于線性模型MLR建模預(yù)測As和Cd等重金屬含量，認為土壤重金屬含量與鐵、鐵氧化物相關(guān)；有研究采用SMLR，PLSR等線性建模方法建立了土壤重金屬含量反演模型；有報道基于逐步回歸和相關(guān)系數(shù)方法，篩選出對重金屬敏感的特征波段，將它們組合成綜合特征變量對研究區(qū)Cu元素進行了反演。盡管國內(nèi)外關(guān)于土壤重金屬含量估算相關(guān)研究逐漸增多，但是仍存在一些問題。例如，針對高光譜數(shù)據(jù)波段信息冗余的問題，多數(shù)研究選擇丟掉大量的波段，僅利用相關(guān)系數(shù)以及逐步回歸法篩選出了部分特征波段，損失了大量有用的信息。事實上，土壤中重金屬含量與光譜曲線之間的關(guān)系很難用幾個波段解釋。因此，選擇一種既可以保證波段主要信息量，又能減少輸入變量的特征參數(shù)尤為重要。此外，關(guān)于土壤重金屬含量估算模型的問題，絕大部分的研究主要采用線性回歸模型，如SMLR和PLSR等；而非線性回歸模型考慮較少。客觀上講，土壤中重金屬含量在光譜曲線上的響應(yīng)會受多種因素影響，二者之間關(guān)系非常復(fù)雜；而簡單線性回歸模型很難處理非線性、隨機性等復(fù)雜的問題。因此，在高光譜模型的選擇上應(yīng)對非線性模型加以考慮。

選擇湖北省黃石市礦山周邊農(nóng)用地土壤為研究對象，針對于高光譜反演中波段信息冗余等問題，提出了基于PCA的降維方法，結(jié)合多種高光譜反演模型，驗證PCA篩選主成分量可實現(xiàn)土壤重金屬含量的高精度反演，并通過不同高光譜模型的對比，確定了適合該研究區(qū)域Cd含量的最佳預(yù)測模型，從而實現(xiàn)了土壤Cd含量的快速、精確光譜檢測，為土壤重金屬反演提供新的思路。

1 實驗部分

1.1 研究區(qū)與土壤采樣

研究區(qū)位于湖北省東南部的黃石市(114°30′—115°30′E，29°30′—30°20′N)，地處長江中下游，具有典型的大陸性季風(fēng)氣候。地勢南高北低，東西平，海拔高度為120～200 m。研究區(qū)內(nèi)礦產(chǎn)資源豐富，有多個大中型礦床，礦山開采、冶煉生產(chǎn)對周邊土壤造成一定的重金屬污染。在研究區(qū)共采集0～20 cm表層土壤56件，采樣點(圖1)位于礦山周邊的農(nóng)用地，采集表層土樣初始質(zhì)量大于1 kg，樣品經(jīng)室內(nèi)自然風(fēng)干、研磨后過10目(孔徑2 mm)的尼龍篩，利用四分法分成兩份，分別用于室內(nèi)光譜測試和實驗室化學(xué)分析。

圖1 研究區(qū)采樣點分布示意圖Fig.1 The distribution of sampling point in the study area

1.2 光譜測定

土壤光譜數(shù)據(jù)獲取采用美國ASD公司的FieldSpec4地物光譜儀(光譜波段范圍350～2 500 nm)，利用鹵素光源和標(biāo)準(zhǔn)白板完成測量。該光譜儀采樣間隔為1.4 nm(350～1 000 nm)和2 nm(1 000～2 500 nm)，經(jīng)光譜重采樣后(間隔1 nm)，共輸出2 151個波段。測試在暗室進行，選擇一穩(wěn)固平臺，將土樣放入直徑90 mm，高19 mm的透明玻璃器皿，使其表面盡量平整，以50 W的鹵素?zé)魹楣庠?，光源與樣品保持50 cm距離，光源探頭位于樣本正上方7 cm高，光線與樣品保持15°的照射角度，保證測量時無陰影遮擋。開機預(yù)熱30 min后對儀器進行調(diào)整和校準(zhǔn)并開始測量。每個土壤樣本采集10條光譜曲線，取光譜反射率的平均值作為樣本的反射率光譜值，剔除350～399和2 450～2 500 nm信噪比低、噪聲大的邊緣波段，共獲得2 050個波段數(shù)據(jù)。

1.3 數(shù)據(jù)處理

1.3.1 光譜預(yù)處理

土壤樣品光譜數(shù)據(jù)測定過程中，由于光線亮度變化和土壤表面凹凸不平會對實驗結(jié)果產(chǎn)生影響，采用取光譜反射率倒數(shù)對數(shù)的方法來避免此影響。倒數(shù)對數(shù)[6]計算公式為

(1)

其中λi為光譜波長值，R(λi)為對應(yīng)光譜波段的反射率。

主成分分析(principal components analysis，PCA)是由Pearson于1901年提出的一種分析、簡化數(shù)據(jù)集的方法[7]。該方法的優(yōu)勢在于降低數(shù)據(jù)集維數(shù)，同時保證信息量最大，對于擁有大量波段信息的高光譜數(shù)據(jù)，通過一系列的矩陣變化，在測量空間尋找?guī)捉M正交向量，保留數(shù)據(jù)方差最大、信息量最多的組分，從而達到高光譜數(shù)據(jù)降維的目的。主要步驟如下：

(1)將波段數(shù)據(jù)組合成為矩陣，設(shè)隨機變量X1，X2，…，XP；其樣本均數(shù)為X1，X2，…，XP；樣本標(biāo)準(zhǔn)差記為S1，S2，…，SP。首先進行標(biāo)準(zhǔn)化變換

(2)

以此類推求得第三，第四，…，第p個主成分。保留主成分個數(shù)取決于累積方差在總方差所占百分比(貢獻率)。

1.3.2 反演模型

利用PCA對原始光譜以及倒數(shù)對數(shù)光譜進行波段降維，將累積貢獻率達到99.99%的主成分作為特征變量，選擇線性模型PLSR，以及非線性模型SVM，ANN和RF分別建立土壤Cd含量估算模型。PLSR是一種常用于高光譜反演土壤元素含量的新型多元統(tǒng)計方法[8]，它能夠很好地解決自變量間多重共線問題。SVM是以內(nèi)核統(tǒng)計學(xué)理論為基礎(chǔ)理論，它的優(yōu)勢主要體現(xiàn)在解決小樣本、非線性以及高維模式的識別[9]。ANN由一組相互連接的人工神經(jīng)元組成，利用大量神經(jīng)元之間的鏈接結(jié)構(gòu)進行分布式并行信息處理的數(shù)學(xué)模型，該模型基本架構(gòu)由輸入層、輸出層和隱藏層三部分構(gòu)成[10]。RF是一個組合分類器算法[11]，由一系列決策樹組成，利用自助法重采樣技術(shù)，在初始樣本數(shù)據(jù)集上生成多個自助樣本集，每個自助樣本集是每棵分類樹的全部訓(xùn)練數(shù)據(jù)，然后根據(jù)自主樣本集生成多個分類樹組成隨機森林。

1.3.3 精度評估

采用R2、RMSE和RPD評價指標(biāo)對估算模型的反演精度進行評估。R2和RPD越大，RMSE越小，說明預(yù)測效果越好，通常認為R2越趨近于1，模型的預(yù)測效果越佳。當(dāng)RPD>2時，模型具極好的預(yù)測能力；當(dāng)1.4

2 結(jié)果與討論

2.1 土壤重金屬含量統(tǒng)計分析

利用等離子體質(zhì)譜法測定Cd含量，Cd元素的描述性統(tǒng)計結(jié)果如表1所示，Cd均值為0.64 mg·kg-1。根據(jù)土壤環(huán)境質(zhì)量標(biāo)準(zhǔn)(GB15618—2018)，該區(qū)域的Cd含量高于農(nóng)用地污染風(fēng)險篩選值，而低于管制值，存在一定的土壤污染風(fēng)險。從空間分布看，其變異系數(shù)介于0.5～0.75之間，屬于中等變異，說明Cd在土壤中分布不均，空間變異較為顯著。將56個樣本數(shù)據(jù)按照7∶3比率隨機分割，訓(xùn)練樣本39個，用于篩選模型輸入變量。驗證樣本17個，用于對高光譜模型的評估。

表1 土壤Cd含量描述統(tǒng)計分析(mg·kg-1)Table 1 Descriptive statistics analysis of soil Cd content

2.2 光譜變換與PCA降維

所有土壤樣本原始光譜反射率曲線[圖2(a)]趨勢大致相同，在可見光區(qū)域反射率呈明顯上升，超過800 nm后光譜曲線趨于平緩。在1 400，1 900和2 200 nm附近有三個明顯凹陷的吸收峰，為土壤黏土礦物的吸收特征。經(jīng)倒數(shù)對數(shù)變化后[圖2(b)]的光譜曲線與原始曲線的變化趨勢基本相反。

圖2 (a)土壤樣本原始光譜反射率曲線；(b)土壤樣本倒數(shù)對數(shù)光譜反射率曲線Fig.2 (a) The original spectral reflectance curve of soil samples； (b) The reciprocal logarithmic spectral reflectance curve of soil samples

利用PCA算法，對原始光譜曲線以及變換后的倒數(shù)對數(shù)光譜曲線的2 050個波段進行降維，原始光譜曲線和倒數(shù)對數(shù)光譜曲線各主成分的貢獻率和累計貢獻率值如表2所示。選取經(jīng)PCA之后，累計貢獻率達到99.99%的主成分個數(shù)作為模型的輸入變量，其中，原始光譜累積貢獻率達到99.99%的主成分個數(shù)為12個，光譜變換之后累積貢獻率達到99.99%的主成分個數(shù)也為12個。將PCA降維選取的組分作為四種模型的輸入變量。

表2 主成分貢獻率Table 2 Principal component contribution rate

2.3 土壤重金屬含量反演模型建立與檢驗

將PCA降維選取的主成分作為模型的自變量(X)，土壤Cd含量為因變量(Y)，采用線性回歸PLSR模型，以及非線性回歸SVM，ANN和RF模型分別建模比對，驗證基于PCA篩選的特征變量對不同模型預(yù)測能力的影響，以及優(yōu)選出研究區(qū)Cd含量的最佳預(yù)測模型。

2.3.1 基于PCA原始光譜建模

利用PCA對原始光譜數(shù)據(jù)降維，選擇累計貢獻率達99.99%的12個主成分量作為模型輸入變量，運用PLSR，SVM，ANN和RF方法建模，土壤Cd含量的反演模型[圖3(a)]的精度評價如表3所示，根據(jù)圖3(a)與表3可知：PCA-RF的決定系數(shù)(R2=0.856)最高，RPD高達3.39，表明PCA-RF模型具有極好的預(yù)測能力，是預(yù)測土壤Cd含量的優(yōu)勢模型；PCA-ANN和PCA-SVM的RPD都高于2，其決定系數(shù)(R2)分別為0.621和0.581，兩種模型同樣具有好的預(yù)測能力；而PCA-PLSR的R2和RPD分別僅為0.484和1.8，該模型的預(yù)測能力一般。經(jīng)PCA降維選取的特征波段，使得模型均具有一定的預(yù)測能力。

2.3.2 基于PCA倒數(shù)對數(shù)光譜建模

運用四種方法對PCA降維后的倒數(shù)對數(shù)光譜進行建模[圖3(b)]，其反演精度如表3。由圖3(b)和表3可知：PCA-RF模型的預(yù)測能力在光譜變換后仍為最佳，其R2為0.855，RPD為3.39，表明模型仍具有極好的預(yù)測能力；PCA-ANN次之，其R2為0.623，RPD為2.12，模型同樣具有好的預(yù)測能力；PCA-SVM的R2為0.607，RPD為2.00，模型也具有好的預(yù)測能力，而PCA-PLSR的R2為0.535，RPD僅為1.89，模型預(yù)測能力一般。

表3 基于原始光譜-倒數(shù)對數(shù)光譜不同模型精度評價Table 3 Accuracy evaluation of different models based on original spectral-reciprocal logarithmic spectral

圖3 (a) 原始光譜不同預(yù)測模型散點圖；(b) 倒數(shù)對數(shù)光譜不同預(yù)測模型散點圖Fig.3 (a) Scatterplots of different prediction models based on original spectral data； (b) Scatter plots between different prediction models based on reciprocal logarithmic spectral

2.3.3 基于PCA原始光譜-倒數(shù)對數(shù)模型對比分析

四種模型的預(yù)測能力順序在光譜變換前后未發(fā)生改變(圖4)，光譜變換對于各模型的預(yù)測能力有所提升，其中提升效果最為顯著的是PCA-PLSR模型，該模型的R2提升了10.5%，RPD提升了5.0%，其次為PCA-SVM模型，該模型的R2提升了4.5%，RPD提升了2.5%，PCA-ANN模型，R2和RPD分別提升了1.8%和1.4%，而PCA-RF模型無明顯改變。

圖4 原始光譜-倒數(shù)對數(shù)對比分析圖Fig.4 The contrast analysis diagram between original spectral and reciprocal logarithmic

通過對比光譜變化前后各模型的預(yù)測精度可得，非線性模型的預(yù)測能力優(yōu)于線性模型，倒數(shù)對數(shù)光譜變換對于模型的預(yù)測能力有所提升，可弱化光譜數(shù)據(jù)測定時光線亮度和土壤表面凹凸產(chǎn)生影響。

3 結(jié) 論

以湖北省黃石市礦區(qū)周邊農(nóng)用地土壤為研究對象，利用PCA方法對光譜變化前后數(shù)據(jù)進行降維，選取特征變量，在此基礎(chǔ)上對比分析了不同反演模型對土壤Cd含量測定的反演精度，得出如下結(jié)論：

(1)經(jīng)倒數(shù)對數(shù)變換后的光譜，預(yù)測能力有所提升，PCA-PLSR模型的提升效果最為明顯，PCA-SVM和PCA-ANN稍有提高，倒數(shù)對數(shù)變換可弱化光譜測定中光強度變化和土壤表面凹凸的影響。

(2)利用PCA方法進行降維處理可以有效降低高光譜數(shù)據(jù)量，選取的12個主成分量對變化前后的光譜累計貢獻率可達99.99%，四種模型均具有一定的預(yù)測能力，保證模型具有極好的輸入變量。

(3)不同模型的反演精度順序為：PCA-RF>PCA-ANN>PCA-SVM>PCA-PLSR，非線性模型PCA-RF，PCA-ANN和PCA-SVM的RPD均大于2，具有極好的預(yù)測能力，其中PCA-RF模型的RPD超過3，說明模型具有較高穩(wěn)定性和預(yù)測精度。

本研究主要采用PCA對光譜數(shù)據(jù)進行降維，對比分析不同模型的反演能力，PCA-RF模型可為土壤重金屬含量反演提供很好的參考依據(jù)。PCA對高光譜數(shù)據(jù)特征變量選取具有顯著效果，但仍存在其他的降維方法，需要進一步深入研究。