基于偏最小二乘的土壤重金屬鎘間接反演模型

賀軍亮，韓超山，韋 銳，周智勇，東啟亮

(1.石家莊學院資源與環境科學學院，石家莊 050035；2.河北省水文工程地質勘查院，石家莊 050021)

0 引言

土壤是不可再生的重要資源，是人類賴以生存的物質基礎，保護土壤資源，防治土壤環境污染，是推進生態文明建設和維護國家生態安全的重要內容[1-2]。隨著工業的發展，過量重金屬通過污水排放、大氣降塵等途徑進入土壤，造成土壤重金屬污染[3]。土壤重金屬不能被土壤微生物分解，具有殘留時間長、難遷移、易累積的特點，而且可通過食物鏈以有害的濃度在人體內蓄積，嚴重危害人體健康和環境安全[4-5]。

土壤重金屬含量的調查監測是進行土壤污染有效防治的前提。高光譜遙感具有快速、無損、光譜分辨率高等特點，可以在土壤定量遙感監測中發揮重要作用[6-7]。Kemper等[8]利用多元線性逐步回歸等方法基于土壤反射光譜與重金屬的相關分析構建模型，反演了Aznalcollar礦區土壤Pb，As，Fe和Hg的含量；郭云開等[9]基于水稻冠層光譜變化特征，利用最小二乘方法擬合建立了Zn，Pb和Cd等土壤重金屬含量的反演模型；程先鋒等[10]在蘭坪鉛鋅礦區利用逐步回歸方法構建了土壤重金屬含量的高光譜估算模型。

一般情況下，土壤重金屬含量較低，屬于痕量級，即使在重污染狀態下，重金屬元素的直接光譜響應也非常微弱[11]。但是，當土壤有機質含量高時，對重金屬的富集具有一定的吸附作用，使得利用有機質光譜特征間接反演重金屬含量成為可能[12]。賀軍亮等[13]利用水稻土有機質的敏感波段構建了有機質光譜診斷指數，用于重金屬Cu和Pb含量間接反演模型的構建；蘭澤英等[14]在樂安河流域采用土壤高光譜數據間接反演了Cu，Zn和Pb的含量，為該區域土壤生態環境監測提供了相關參考。除原始光譜反射率外，微分、倒數和對數等多種光譜變換指標能夠從不同角度和程度上突出反演對象的光譜特征，降低背景噪聲的影響[13-15]。前人所構建的重金屬最優反演模型大多采用其中一種指標作為自變量建模，多種光譜變換指標的集成建模有可能會提高模型的估算精度和穩定度。

本文以石家莊市地表水源保護區土壤重金屬Cd為研究對象，基于間接反演的研究機理，嘗試通過多種光譜變換和相關性分析方法提取有機質光譜診斷特征，并建立多指標集成估算模型來間接反演重金屬Cd的含量，進一步豐富土壤重金屬高光譜反演研究的案例，以期為土壤定量遙感研究提供科學參考。

1 研究區及數據源

1.1 研究區概況

石家莊市地表水源保護區位于石家莊市西部的平山縣和井陘縣境內，保護區范圍為N38°01′～38°45′，E113°34′～114°18′。一級保護區包括滹沱河干流、崗南水庫和黃壁莊水庫。一級保護區之外按緩沖距離和行政區劃設二級和三級保護區。整個保護區地勢西高東低，高差懸殊。雖然該保護區是石家莊市生活用水的主要供應地，也是北京市備用水源地，但是該區礦產開發歷史悠久，人地矛盾突出，鐵礦和石灰礦等礦產資源開采存在尾礦威脅、植被破壞、水土流失和土壤污染等狀況，對保護區生態安全造成了一定影響。考慮到土壤空間分布的不均勻性及研究區地形的影響，參考研究區內主要采礦點和主要冶煉企業的區位，在研究區共采集69個土壤樣品，土壤類型以褐土為主，采集深度為0～20 cm，同時利用GPS對采樣點進行了定位(圖1)。

圖1 研究區采樣點分布Fig.1 Sampling point distribution in the study area

所有樣點采集的土壤樣品帶回室內后經風干、研磨、過篩(100目)后備用。土壤有機質含量采用重鉻酸鉀法測定，重金屬Cd全量利用電感耦合等離子體質譜儀ICP-MS檢測方法測定。

1.2 土壤光譜采集

土壤反射光譜測量在南京師范大學地理科學學院實驗室完成。土壤反射率采用美國ASD公司生產的FieldSpec Pro便攜式地物光譜儀測定，光譜采集范圍為350～2 500 nm。光譜采集工作在室內進行，把鹵素燈作為唯一光源。將土壤樣品置于放在黑色絨布上的直徑10 cm、深2 cm玻璃培養皿內，用工具將表面刮平。光源入射角為45°，距離土壤樣品30 cm，采集槍垂直于土壤樣品并保持15 cm的距離。每次測定之前都需要進行白板校正，每個樣品測10次光譜曲線，取其平均值[16]。受檢測周期長、土樣殘留水分以及光譜測量背景環境等影響，全光譜范圍兩端噪聲較大，剔除光譜兩端和水汽吸收影響較大的波段，最終保留1 640個有效波段進行后續建模研究。

1.3 光譜數據預處理

光譜變換可以很大程度上消除土壤背景的影響，進一步提高光譜信噪比，突出光譜的吸收和反射特征[17]。本文光譜實測值為土壤光譜反射率，采用以下方法對原始光譜反射率進行變換：一階微分(first derivative，FD)、二階微分(second derivative，SD)、倒數變換(reciprocal transformation，RT)、倒數的一階微分(reciprocal transformation and first derivative，RTFD)、倒數的二階微分(reciprocal transformation and second derivative，RTSD)、倒數的對數變換(absorbance transformation，AT)、倒數對數的一階微分(absorbance transformation and first derivative，ATFD)、倒數對數的二階微分(absorbance transformation and second derivative，ATSD)、連續統去除(continuum removal，CR)。將有機質含量實測值與以上各光譜數據進行相關性分析，在每項變換通過0.01水平顯著性檢驗的波段中，找出相關系數絕對值最大處所對應的波段，作為該變換下有機質的敏感波段。光譜變換工作在Origin2017和ENVI5.3軟件中進行，相關性分析在SPSS24軟件中進行。

2 研究方法

2.1 建模方法選擇

FD，SD，RT和AT等多種光譜變換指標能夠從不同角度和程度上突出反演對象的光譜特征，考慮間接反演誤差累積的可能影響，有必要對以上各光譜變換指標進行篩選。篩選方法采用多元線性逐步回歸方法(multiple linear stepwise regression，MLSR)，保留對模型精度影響較大的變換指標。

反演模型的構建采用偏最小二乘回歸方法。偏最小二乘法是一種數學優化技術，通過最小化誤差的平方和找到一組數據的最佳函數匹配[18]。前人所構建的重金屬最優反演模型大多只采用一種光譜指標作為自變量建模，即單光譜變換指標偏最小二乘(univariate partial least squares regression，U-PLSR)模型。多種光譜變換數據集成建模能夠綜合反映各種光譜變換形式的特點，通過嘗試建立多光譜變換指標偏最小二乘(multivariate partial least squares regression，M-PLSR)模型，并與U-PLSR模型進行對比分析。

2.2 模型精度評價

Rank-KS方法可有效提升建模樣本與驗證樣本選擇的合理性[19]。采用Rank-KS方法將69個樣本分為2組，其中46個樣本用來建模，23個樣本用來驗證模型精度。在模型精度分析過程中，主要參考擬合優度調節參數R2，并結合均方根誤差(root mean square error，RMSE)和相對分析誤差(relative percent deviation，RPD)等參數進行模型對比。R2和RPD越高，RMSE越小，說明模型的擬合程度越好。RPD為樣本標準差和RMSE的比值，可以用來判斷模型的預測能力。當RPD<1.4時，認為該模型不具有預測能力；當1.4≤RPD<2時，可以對樣本進行粗略預測；當RPD≥2時，模型具有極好的預測能力[20]。

3 結果分析

3.1 Cd特征量統計與污染現狀分析

表1為Cd的特征量統計值。所采集土壤樣本中，Cd的含量分布范圍為0.138～0.359 mg/kg，平均值達0.220 mg/kg。總體樣本、建模樣本和驗證樣本的各統計量值都較為接近，標準差較小，說明Rank-KS選取樣本較為合理，分類數據具有一定的代表性。

表1 Cd特征量統計Tab.1 Statistical characteristics of heavy metal Cd in soil

采用單因子污染指數評價法[21]對研究區Cd的污染現狀進行評價分析，其計算公式為

Pi=Ci/Si，

(1)

式中：Pi為土壤中污染物的環境質量指數；Ci為實測值；Si為背景值。

石家莊地區Cd的背景值為0.08 mg/kg，當污染指數P≥1時，可以認定該地區存在Cd污染，指數越大污染狀況越嚴重[22]。根據統計可以得出所有樣本的污染指數P>1，92.8%的樣本的污染指數P>2。因此，土壤樣本采集區域存在較為嚴重的Cd污染，需要加強對該地區Cd含量的全面監測。

表2 土壤重金屬Cd的污染指數統計Tab.2 Statistical analysis of pollution index of heavy metal Cd in soil

3.2 MLSR模型分析

經計算，總體樣本的重金屬Cd和有機質含量的相關系數達到0.7，并通過了0.01水平的顯著性檢驗，說明兩者之間存在一定的吸附賦存關系[12,14]。根據吸附機理，將有機質含量與各光譜變換指標進行相關性分析，分析結果如表3所示。

表3 土壤有機質含量與光譜變量的最大相關系數Tab.3 Maximum correlation coefficients of soil organic matter content and spectral variables

①R為光譜反射率；②**表示通過0.01水平的顯著性檢驗。

從表3可以看出，原始光譜反射率在797 nm波段處與有機質含量存在最大相關性，與前人研究結果一致[23]，ATFD與有機質含量的相關性最大，且呈負相關，FD與有機質含量存在最大的正相關關系。根據有機質含量與光譜變量的相關分析，將各敏感波段對應的光譜變換指標作為Cd估算模型的自變量因子，Cd含量實測值作為因變量，構建MLSR模型。經MLSR分析，保留對模型精度影響較大的變換指標(ATSD1409和FD1396)，模型散點圖及驗證結果如圖2和表4所示，具體模型為

Y=0.432+1 389.565XATSD-66.253XFD。

(2)

(a)建模樣本 (b)驗證樣本

圖2 MLSR模型散點圖

Fig.2ScatterplotsofthetwosetsofsamplesoffittingandtestingforMLSR

表4 MLSR模型結果Tab.4 Results of MLSR model

建模樣本MLSR模型R2為0.81，說明該模型對數據具有良好的解釋能力。其RMSE僅為0.02，RPD為2.31，說明模型誤差較小并具有良好的預測能力。而驗證樣本在R2和RMSE與建模樣本相近的情況下，RPD卻只有1.23，模型預測能力較差。通過圖2可以看出，建模樣本MLSR模型的趨勢線與1∶1線夾角較小，模型預測值與實測值較為接近；驗證樣本的趨勢線與1∶1線夾角較大，說明模型存在一定誤差，模型穩定度有待進一步提高。

3.3 U-PLSR模型分析

經過3.2節MLSR分析，篩選出對模型精度影響較大的變換指標(ATSD和FD)。分別采用ATSD和FD光譜指標作為自變量，構建U-PLSR模型。ATSD-U-PLSR模型散點圖及驗證結果如圖3和表5所示，具體模型為

Y=0.445 2+1 271.535XATSD。

(3)

(a)建模樣本 (b)驗證樣本

圖3 ATSD-U-PLSR模型散點圖

Fig.3ScatterplotsofthetwosetsofsamplesoffittingandtestingforATSD-U-PLSR

表5 ATSD-U-PLSR模型結果Tab.5 Results of ATSD-U-PLSR model

與MLSR模型結果對比，ATSD-U-PLSR建模樣本的R2降低了0.03，RMSE提高了0.017，重金屬Cd實測值和預測值的擬合度下降，但RPD提高了12.09；驗證樣本R2，RMSE和RPD分別提高了0.05，0.024和7.87，相比MLSR模型數據預測能力有大幅提升。

FD-U-PLSR模型散點圖及驗證結果如圖4和表6所示，具體模型為

Y=0.287 9+151.316 6XFD。

(4)

(a)建模樣本 (b)驗證樣本

圖4 FD-U-PLSR模型散點圖

Fig.4ScatterplotsofthetwosetsofsamplesoffittingandtestingforFD-U-PLSR

表6 FD-U-PLSR模型結果Tab.6 Results of FD-U-PLSR model

從圖4和表6可以看出，FD-U-PLSR模型建模樣本與驗證樣本的擬合度均過低，R2僅分別為0.24和0.42。雖然總體樣本FD光譜變量與有機質含量存在最大的正相關關系，但是相對于集成多光譜變換指標構建的MLSR模型，FD-U-PLSR模型估算誤差較大，反映出單光譜變換指標估算模型的不穩定性。

3.4 M-PLSR模型分析

在MLSR和U-PLSR模型對比分析的基礎上，考慮不同光譜變換形式提高光譜信噪比的能力差異，分別基于建模樣本和驗證樣本數據集，進一步建立M-PLSR模型。所建模型解釋變量與MLSR模型一致，仍為ATSD與FD。M-PLSR模型散點圖及驗證結果如圖5和表7所示，具體模型為

Y=0.452 1+1 379.77XATSD-34.186 1XFD。

(5)

(a)建模樣本 (b)驗證樣本

圖5 M-PLSR模型散點圖

Fig.5ScatterplotsofthetwosetsofsamplesoffittingandtestingforM-PLSR

表7 M-PLSR模型結果Tab.7 Results of M-PLSR model

M-PLSR模型建模樣本與驗證樣本R2都超過了0.8，RPD都超過了2，RMSE也較低，并且驗證樣本與建模樣本的預測精度相近。與MLSR模型和U-PLSR模型相比，M-PLSR模型擬合優度有所提高，模型穩定度得到明顯增強。

4 結論與討論

以石家莊市水源保護區褐土為研究對象，采用偏最小二乘法構建土壤重金屬Cd的高光譜間接反演模型，具體結論如下：

1)石家莊市水源保護區內土壤樣本采集區重金屬Cd含量的平均值為0.220 mg/kg，高于背景值，并且總體樣本的污染指數均大于1，存在較為嚴重的Cd污染。

2)研究區土壤重金屬Cd和有機質含量在0.01水平上顯著相關，相關系數達到0.7，兩者之間存在一定的吸附賦存關系，基于土壤有機質的光譜特征間接反演重金屬Cd含量具有一定的合理性。

3)有機質的原始光譜反射率對應的敏感波段為797 nm，各種光譜變換中倒數對數的一階微分與有機質含量的相關性最大，可達到-0.766，一階微分與有機質含量存在最大的正相關關系，達到0.721。

4)采用多元線性逐步回歸方法篩選確定的模型解釋變量為ATSD1409和FD1396。對于所建多光譜變換指標偏最小二乘模型，其建模樣本和驗證樣本的R2分別為0.83和0.80，模型預測值與重金屬Cd含量實測值較為接近，整體擬合效果優于單光譜變換指標模型。

與以往基于單光譜變換指標的建模方法研究不同，采用多種光譜變換指標進行集成建模，具有更好的建模效果和預測能力，可以對重金屬Cd含量起到較好的預測作用。然而，由于影響土壤重金屬反射光譜特征的因素不是單一的，僅僅考慮重金屬與有機質的光譜響應特征還不全面。因此，在進一步研究工作中，可以嘗試考慮其他因素如鐵錳氧化物的綜合影響，提取多因素特征波段構建多光譜變換指標估算模型。

此外，由于氣候、地質地貌等條件的影響，不同地區土壤具有其獨特的區域性。選取石家莊市水源保護區褐土為研究樣本，土樣經過研磨、風干處理后，雖然基本消除了土壤質地、土壤濕度等對土壤光譜的影響，但是所建立的土壤重金屬Cd含量的估算模型是否適用于其他地區，有待進一步研究。