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

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

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

0 引言

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

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

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

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

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

1.1 研究區(qū)概況

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

圖1 研究區(qū)采樣點(diǎn)分布Fig.1 Sampling point distribution in the study area

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

1.2 土壤光譜采集

土壤反射光譜測(cè)量在南京師范大學(xué)地理科學(xué)學(xué)院實(shí)驗(yàn)室完成。土壤反射率采用美國(guó)ASD公司生產(chǎn)的FieldSpec Pro便攜式地物光譜儀測(cè)定，光譜采集范圍為350～2 500 nm。光譜采集工作在室內(nèi)進(jìn)行，把鹵素?zé)糇鳛槲ㄒ还庠础⑼寥罉悠分糜诜旁诤谏q布上的直徑10 cm、深2 cm玻璃培養(yǎng)皿內(nèi)，用工具將表面刮平。光源入射角為45°，距離土壤樣品30 cm，采集槍垂直于土壤樣品并保持15 cm的距離。每次測(cè)定之前都需要進(jìn)行白板校正，每個(gè)樣品測(cè)10次光譜曲線，取其平均值[16]。受檢測(cè)周期長(zhǎng)、土樣殘留水分以及光譜測(cè)量背景環(huán)境等影響，全光譜范圍兩端噪聲較大，剔除光譜兩端和水汽吸收影響較大的波段，最終保留1 640個(gè)有效波段進(jìn)行后續(xù)建模研究。

1.3 光譜數(shù)據(jù)預(yù)處理

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

2 研究方法

2.1 建模方法選擇

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

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

2.2 模型精度評(píng)價(jià)

Rank-KS方法可有效提升建模樣本與驗(yàn)證樣本選擇的合理性[19]。采用Rank-KS方法將69個(gè)樣本分為2組，其中46個(gè)樣本用來(lái)建模，23個(gè)樣本用來(lái)驗(yàn)證模型精度。在模型精度分析過(guò)程中，主要參考擬合優(yōu)度調(diào)節(jié)參數(shù)R2，并結(jié)合均方根誤差(root mean square error，RMSE)和相對(duì)分析誤差(relative percent deviation，RPD)等參數(shù)進(jìn)行模型對(duì)比。R2和RPD越高，RMSE越小，說(shuō)明模型的擬合程度越好。RPD為樣本標(biāo)準(zhǔn)差和RMSE的比值，可以用來(lái)判斷模型的預(yù)測(cè)能力。當(dāng)RPD<1.4時(shí)，認(rèn)為該模型不具有預(yù)測(cè)能力；當(dāng)1.4≤RPD<2時(shí)，可以對(duì)樣本進(jìn)行粗略預(yù)測(cè)；當(dāng)RPD≥2時(shí)，模型具有極好的預(yù)測(cè)能力[20]。

3 結(jié)果分析

3.1 Cd特征量統(tǒng)計(jì)與污染現(xiàn)狀分析

表1為Cd的特征量統(tǒng)計(jì)值。所采集土壤樣本中，Cd的含量分布范圍為0.138～0.359 mg/kg，平均值達(dá)0.220 mg/kg。總體樣本、建模樣本和驗(yàn)證樣本的各統(tǒng)計(jì)量值都較為接近，標(biāo)準(zhǔn)差較小，說(shuō)明Rank-KS選取樣本較為合理，分類(lèi)數(shù)據(jù)具有一定的代表性。

表1 Cd特征量統(tǒng)計(jì)Tab.1 Statistical characteristics of heavy metal Cd in soil

采用單因子污染指數(shù)評(píng)價(jià)法[21]對(duì)研究區(qū)Cd的污染現(xiàn)狀進(jìn)行評(píng)價(jià)分析，其計(jì)算公式為

Pi=Ci/Si，

(1)

式中：Pi為土壤中污染物的環(huán)境質(zhì)量指數(shù)；Ci為實(shí)測(cè)值；Si為背景值。

石家莊地區(qū)Cd的背景值為0.08 mg/kg，當(dāng)污染指數(shù)P≥1時(shí)，可以認(rèn)定該地區(qū)存在Cd污染，指數(shù)越大污染狀況越嚴(yán)重[22]。根據(jù)統(tǒng)計(jì)可以得出所有樣本的污染指數(shù)P>1，92.8%的樣本的污染指數(shù)P>2。因此，土壤樣本采集區(qū)域存在較為嚴(yán)重的Cd污染，需要加強(qiáng)對(duì)該地區(qū)Cd含量的全面監(jiān)測(cè)。

表2 土壤重金屬Cd的污染指數(shù)統(tǒng)計(jì)Tab.2 Statistical analysis of pollution index of heavy metal Cd in soil

3.2 MLSR模型分析

經(jīng)計(jì)算，總體樣本的重金屬Cd和有機(jī)質(zhì)含量的相關(guān)系數(shù)達(dá)到0.7，并通過(guò)了0.01水平的顯著性檢驗(yàn)，說(shuō)明兩者之間存在一定的吸附賦存關(guān)系[12,14]。根據(jù)吸附機(jī)理，將有機(jī)質(zhì)含量與各光譜變換指標(biāo)進(jìn)行相關(guān)性分析，分析結(jié)果如表3所示。

表3 土壤有機(jī)質(zhì)含量與光譜變量的最大相關(guān)系數(shù)Tab.3 Maximum correlation coefficients of soil organic matter content and spectral variables

①R為光譜反射率；②**表示通過(guò)0.01水平的顯著性檢驗(yàn)。

從表3可以看出，原始光譜反射率在797 nm波段處與有機(jī)質(zhì)含量存在最大相關(guān)性，與前人研究結(jié)果一致[23]，ATFD與有機(jī)質(zhì)含量的相關(guān)性最大，且呈負(fù)相關(guān)，F(xiàn)D與有機(jī)質(zhì)含量存在最大的正相關(guān)關(guān)系。根據(jù)有機(jī)質(zhì)含量與光譜變量的相關(guān)分析，將各敏感波段對(duì)應(yīng)的光譜變換指標(biāo)作為Cd估算模型的自變量因子，Cd含量實(shí)測(cè)值作為因變量，構(gòu)建MLSR模型。經(jīng)MLSR分析，保留對(duì)模型精度影響較大的變換指標(biāo)(ATSD1409和FD1396)，模型散點(diǎn)圖及驗(yàn)證結(jié)果如圖2和表4所示，具體模型為

Y=0.432+1 389.565XATSD-66.253XFD。

(2)

(a)建模樣本 (b)驗(yàn)證樣本

圖2 MLSR模型散點(diǎn)圖

Fig.2ScatterplotsofthetwosetsofsamplesoffittingandtestingforMLSR

表4 MLSR模型結(jié)果Tab.4 Results of MLSR model

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

3.3 U-PLSR模型分析

經(jīng)過(guò)3.2節(jié)MLSR分析，篩選出對(duì)模型精度影響較大的變換指標(biāo)(ATSD和FD)。分別采用ATSD和FD光譜指標(biāo)作為自變量，構(gòu)建U-PLSR模型。ATSD-U-PLSR模型散點(diǎn)圖及驗(yàn)證結(jié)果如圖3和表5所示，具體模型為

Y=0.445 2+1 271.535XATSD。

(3)

(a)建模樣本 (b)驗(yàn)證樣本

圖3 ATSD-U-PLSR模型散點(diǎn)圖

Fig.3ScatterplotsofthetwosetsofsamplesoffittingandtestingforATSD-U-PLSR

表5 ATSD-U-PLSR模型結(jié)果Tab.5 Results of ATSD-U-PLSR model

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

FD-U-PLSR模型散點(diǎn)圖及驗(yàn)證結(jié)果如圖4和表6所示，具體模型為

Y=0.287 9+151.316 6XFD。

(4)

(a)建模樣本 (b)驗(yàn)證樣本

圖4 FD-U-PLSR模型散點(diǎn)圖

Fig.4ScatterplotsofthetwosetsofsamplesoffittingandtestingforFD-U-PLSR

表6 FD-U-PLSR模型結(jié)果Tab.6 Results of FD-U-PLSR model

從圖4和表6可以看出，F(xiàn)D-U-PLSR模型建模樣本與驗(yàn)證樣本的擬合度均過(guò)低，R2僅分別為0.24和0.42。雖然總體樣本FD光譜變量與有機(jī)質(zhì)含量存在最大的正相關(guān)關(guān)系，但是相對(duì)于集成多光譜變換指標(biāo)構(gòu)建的MLSR模型，F(xiàn)D-U-PLSR模型估算誤差較大，反映出單光譜變換指標(biāo)估算模型的不穩(wěn)定性。

3.4 M-PLSR模型分析

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

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

(5)

(a)建模樣本 (b)驗(yàn)證樣本

圖5 M-PLSR模型散點(diǎn)圖

Fig.5ScatterplotsofthetwosetsofsamplesoffittingandtestingforM-PLSR

表7 M-PLSR模型結(jié)果Tab.7 Results of M-PLSR model

M-PLSR模型建模樣本與驗(yàn)證樣本R2都超過(guò)了0.8，RPD都超過(guò)了2，RMSE也較低，并且驗(yàn)證樣本與建模樣本的預(yù)測(cè)精度相近。與MLSR模型和U-PLSR模型相比，M-PLSR模型擬合優(yōu)度有所提高，模型穩(wěn)定度得到明顯增強(qiáng)。

4 結(jié)論與討論

以石家莊市水源保護(hù)區(qū)褐土為研究對(duì)象，采用偏最小二乘法構(gòu)建土壤重金屬Cd的高光譜間接反演模型，具體結(jié)論如下：

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

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

3)有機(jī)質(zhì)的原始光譜反射率對(duì)應(yīng)的敏感波段為797 nm，各種光譜變換中倒數(shù)對(duì)數(shù)的一階微分與有機(jī)質(zhì)含量的相關(guān)性最大，可達(dá)到-0.766，一階微分與有機(jī)質(zhì)含量存在最大的正相關(guān)關(guān)系，達(dá)到0.721。

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

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

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