土壤重金屬砷的高光譜估算模型

李志遠(yuǎn)， 鄧 帆， 賀軍亮， 魏 薇

1. 長(zhǎng)江大學(xué)地球科學(xué)學(xué)院,湖北 武漢 430100 2. 石家莊學(xué)院資源與環(huán)境科學(xué)學(xué)院,河北 石家莊 050035

引 言

土壤重金屬監(jiān)測(cè)在農(nóng)業(yè)生產(chǎn)中至關(guān)重要， 尤其工業(yè)和農(nóng)業(yè)密集地區(qū)存在重金屬擴(kuò)散到食物鏈的風(fēng)險(xiǎn)[1]。 砷(As)是人體所需微量元素之一， 但若吸入出現(xiàn)超標(biāo)， 會(huì)對(duì)人體健康造成極大危害。 含砷三廢的排放和農(nóng)藥化肥的過(guò)度使用人為加劇了重金屬砷污染[2]。

傳統(tǒng)的土壤特性檢測(cè)方法能夠準(zhǔn)確測(cè)得重金屬含量， 但是花費(fèi)高、 步驟繁瑣、 耗時(shí)長(zhǎng)， 無(wú)法做到實(shí)時(shí)高效[3]。 高光譜技術(shù)具有成本低、 效率高、 用時(shí)短等特點(diǎn)， 被廣泛用于重金屬含量估算研究中[4]。 目前， 用于土壤重金屬As高光譜反演的線(xiàn)性建模方法主要有單變量回歸、 多元線(xiàn)性逐步回歸、 主成分回歸、 偏最小二乘回歸、 地理加權(quán)回歸等， 這些方法在對(duì)應(yīng)的研究區(qū)域取得了較好的效果。

多元線(xiàn)性逐步回歸(multiple linear stepwise regression,MLSR)和偏最小二乘回歸(partial least squares regression,PLSR)在估算重金屬As含量中運(yùn)用非常普遍， 它盡可能選出對(duì)重金屬As含量影響較大的自變量， 但是容易出現(xiàn)過(guò)度擬合情況， 導(dǎo)致模型穩(wěn)定性降低[5]； 偏最小二乘回歸綜合了主成分回歸和多元線(xiàn)性逐步回歸方法， 能夠在波段個(gè)數(shù)較多且自相關(guān)嚴(yán)重時(shí)建立回歸模型， 避免共線(xiàn)性問(wèn)題[5]。

以采礦地和冶煉企業(yè)較多的石家莊市地表水源地保護(hù)區(qū)褐土為研究對(duì)象， 主要根據(jù)重金屬As含量與光譜指標(biāo)進(jìn)行相關(guān)分析選出特征波段， 運(yùn)用多元線(xiàn)性逐步回歸和偏最小二乘回歸方法來(lái)估算土壤重金屬As含量， 找出反演重金屬As的優(yōu)選方法， 對(duì)研究該區(qū)域土壤As污染的快速監(jiān)測(cè)提供借鑒。

1 實(shí)驗(yàn)部分

1.1 研究區(qū)概況

研究區(qū)位于石家莊市滹沱河中游的平山縣和井陘縣境內(nèi)， 主體(一級(jí)保護(hù)區(qū))包括滹沱河干流、 黃壁莊水庫(kù)和崗南水庫(kù)及其支流， 主體外設(shè)有二級(jí)和三級(jí)保護(hù)區(qū)(如圖1)。 總體呈西北向東南逐漸降低， 地勢(shì)差距較大。 由于樣本采集區(qū)域開(kāi)發(fā)較早， 人地矛盾突出， 礦產(chǎn)的開(kāi)采造成了土壤重金屬污染問(wèn)題， 對(duì)研究區(qū)現(xiàn)有生態(tài)安全造成一定的負(fù)面影響[6]。

圖1 研究區(qū)和土壤采樣點(diǎn)分布圖Fig.1 Study area and the distributionof the soil sampling points

1.2 樣本采集

基于研究區(qū)地形和土壤空間分布不均勻的特性， 依據(jù)研究區(qū)內(nèi)主要采礦點(diǎn)和冶煉企業(yè)的分布， 按照地表徑流路徑設(shè)置三條樣帶， 采集土壤樣本48個(gè)， 采集深度為0～20 cm， 同時(shí)使用GPS定位， 采樣點(diǎn)位置如圖1所示。 采集的所有土壤樣本帶回室內(nèi)， 經(jīng)風(fēng)干、 研磨、 過(guò)100目篩， 變成粉狀土， 然后將每份樣本一分為二， 分別進(jìn)行土壤基本理化性質(zhì)和高光譜數(shù)據(jù)的測(cè)定。 采用微波消解儀(CEM-Mars)消解， 利用電感耦合等離子體質(zhì)譜儀(PE-Nexion 300 ICP-MS)測(cè)定土壤中的砷含量全量。 表1為48個(gè)土壤樣本中砷含量的描述性統(tǒng)計(jì)值。

表1 土壤中砷含量的描述性統(tǒng)計(jì)值(mg·kg-1)Table 1 Descriptive statistics of the arseniccontent (mg·kg-1) in soil

研究區(qū)土壤中重金屬As含量分布在10.844～15.125 mg·kg-1之間， 平均值為13.230 mg·kg-1， 接近河北省背景值[7]， 說(shuō)明采集的部分土壤樣本中重金屬As含量已經(jīng)富集到了一定程度。 土壤重金屬污染程度采用單因子評(píng)價(jià)法[8]來(lái)分析， 計(jì)算公式為

Pi=Ci/Si

(1)

其中Pi表示土壤中重金屬As的污染指數(shù)， Ci表示土壤中重金屬As的實(shí)測(cè)值， Si表示土壤中重金屬As的河北省背景值。 Pi≤0.7為安全范圍， 0.7

由表2統(tǒng)計(jì)得到， 污染指數(shù)的最小值為0.797， 平均值為0.973， 最大值為1.112。 有29個(gè)樣點(diǎn)處于警戒范圍； 有19個(gè)樣點(diǎn)受到重金屬As的輕度污染[9]。 綜合看來(lái)， 約40%土壤樣本受到重金屬As的輕度污染， 受到了采礦點(diǎn)和冶煉企業(yè)的影響， 有必要對(duì)重金屬As含量進(jìn)行長(zhǎng)期動(dòng)態(tài)監(jiān)測(cè)。

表2 土壤中砷的單因子污染指數(shù)統(tǒng)計(jì)結(jié)果Table 2 Statistical results of single factorpollution index of the soil arsenic

1.3 光譜測(cè)定

采用FieldSpec Pro便攜式光譜儀測(cè)定土壤光譜反射率， 波段范圍為350～2 500 nm， 采樣間隔為1.4 nm(350～1 000 nm)和2 nm(1 000～2 500 nm)， 重采樣間隔為1 nm。 光譜測(cè)量在可調(diào)節(jié)光照的暗室內(nèi)進(jìn)行， 為盡可能減少外來(lái)光線(xiàn)的影響， 選用鹵素?zé)糇鳛槲ㄒ还庠础?將處理好的土壤樣本放入玻璃培養(yǎng)皿(直徑10 cm、 深2 cm)內(nèi)， 置于黑色絨布上， 用工具將土壤表面刮平。 光源天頂角為45°， 距土樣表面30 cm， 探頭置于距土壤樣本表面垂直上方的15 cm處， 光纖探頭視場(chǎng)角為5°。 每次測(cè)試前進(jìn)行白板標(biāo)定， 每個(gè)土壤樣本測(cè)試10次光譜曲線(xiàn)， 取算數(shù)平均值作為該樣本實(shí)際光譜數(shù)據(jù)。

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

為消除儀器信噪比、 樣本風(fēng)干和雜散光等因素的影響， 采用Savitzky-Golay 7點(diǎn)平滑， 去除噪聲較大的異常波段， 保留1 640個(gè)波段用于后續(xù)建模研究。 為了提高光譜信噪比， 去除背景噪聲， 突出光譜特征， 通過(guò)原始光譜反射率(SR)分別計(jì)算一階微分(FD)、 二階微分(SD)、 倒數(shù)(RT)、 倒數(shù)一階微分(RTFD)、 倒數(shù)二階微分(RTSD)、 倒數(shù)對(duì)數(shù)(AT)、 倒數(shù)對(duì)數(shù)一階微分(ATFD)、 倒數(shù)對(duì)數(shù)二階微分(ATSD)、 連續(xù)統(tǒng)去除(CR)。 將重金屬實(shí)測(cè)值與以上經(jīng)光譜變換后的數(shù)據(jù)進(jìn)行相關(guān)性檢驗(yàn)， 得到相關(guān)性最高的波段。 光譜變換由Origin2020b和ENVI5.1軟件進(jìn)行， 結(jié)果見(jiàn)圖2。 相關(guān)性分析及后續(xù)模型建立在Minitab18軟件中完成。

圖2 光譜變換后的土壤樣本光譜反射率曲線(xiàn)Fig.2 Transformed soil reflectance spectra after

1.5 多變量反演方法

(1) MLSR

多元線(xiàn)性逐步回歸是通過(guò)一個(gè)因變量和多個(gè)自變量組合表達(dá)它們之間的線(xiàn)性關(guān)系。 回歸方程為

y=β0+β1x1+β2x2+…+βkxk+ε

(2)

其中， y表示因變量， x表示自變量， βi(i=0,1,2,…,k)表示回歸系數(shù)， ε表示隨機(jī)誤差， k表示自變量的數(shù)量[10]。 將土壤重金屬As實(shí)測(cè)含量作為建模樣本的因變量， 將光譜指標(biāo)作為自變量， 根據(jù)貢獻(xiàn)率高低來(lái)建立MLSR模型。

(2) PLSR

偏最小二乘回歸是一種結(jié)合主成分回歸和MLSR的分析方法， 由Herman Wold于1966年提出， 目前已被廣泛應(yīng)用于光譜數(shù)據(jù)處理。

PLSR的主要研究?jī)?nèi)容是在大量的兩組高線(xiàn)性相關(guān)變量的情況下， 建立自變量的線(xiàn)性模型， 以解決樣本數(shù)小于變量數(shù)的問(wèn)題， 避免出現(xiàn)過(guò)擬合， 因此PLSR是對(duì)MLSR的改進(jìn)。 PLSR的原理是： 首先從自變量(x1,x2,…,xm)中提取相互獨(dú)立的分量Th(h=1,2,…)， 并且提取的主成分?jǐn)y帶盡可能多的原始成分； 然后從因變量(y1,y2,…,ym)中提取獨(dú)立分量Uh(h=1,2,…)， 要求Th和Uh之間的協(xié)方差最大化， 并利用MLSR得出的因變量建立回歸模型。 偏最小二乘回歸的基本模型是

X=ThPT+E

(3)

Y=UhQT+F

(4)

其中,P和Q分別是m×h的正交載荷矩陣，E和F是誤差項(xiàng)(服從正態(tài)分布的隨機(jī)變量)[11]。 以土壤重金屬As實(shí)測(cè)含量為因變量， 選用單一光譜指標(biāo)的特征波段作為自變量建立單光譜變換指標(biāo)偏最小二乘回歸(univariate partial least squares regression,U-PLSR)模型， 用MLSR篩選出的多種光譜指標(biāo)作為自變量建立多光譜變換指標(biāo)偏最小二乘回歸(multivariate partial least squares regression,M-PLSR)模型。

1.6 模型建立與檢驗(yàn)

隨機(jī)將48個(gè)樣本分為兩部分， 32個(gè)用作建立模型， 16個(gè)用來(lái)驗(yàn)證模型。 采用相關(guān)系數(shù)r、 均方根誤差(root mean square error,RMSE)和統(tǒng)計(jì)值F來(lái)衡量建模效果。F值和r越大， RMSE越小， 建模樣本和驗(yàn)證樣本的散點(diǎn)分布越靠近1∶1線(xiàn)， 建模效果越好， 由此確定最佳回歸模型。

2 結(jié)果與討論

2.1 研究區(qū)土壤光譜特性

原始光譜曲線(xiàn)經(jīng)過(guò)平滑之后得到的曲線(xiàn)如圖3所示， 48個(gè)土壤樣本的光譜曲線(xiàn)形態(tài)基本相同。 在可見(jiàn)光波段， 光譜反射率隨波長(zhǎng)增長(zhǎng)而增大， 在近紅外波段， 整體趨于平緩。 光譜曲線(xiàn)在600和800 nm附近有機(jī)質(zhì)的反射峰； 1 400， 1 900和2 200 nm波段附近存在三處明顯的水分吸收帶， 這是受到土壤中硅酸鹽礦物和粘土礦物有關(guān)[12]。 FD， SD和ATFD變換光譜曲線(xiàn)起伏較大； RT和AT變換呈現(xiàn)從起始波段迅速下降， 到600 nm附近降幅變緩趨勢(shì)； RTFD， RTSD和ATSD變換后反射率變化主要出現(xiàn)在500 nm之前； 經(jīng)CR變換， 光譜曲線(xiàn)歸一化在0～1之間， 吸收谷出現(xiàn)在500， 1 010， 1 400， 1 900和2 200 nm附近(圖3)。

圖3 平滑后的土壤樣本光譜反射率曲線(xiàn)Fig.3 Smoothed reflectance spectra of soil samples

2.2 特征波段提取

為了探究不同光譜指標(biāo)與重金屬砷含量的關(guān)系， 將As含量與各種光譜變換進(jìn)行Pearson相關(guān)性分析， 表3表示最大相關(guān)系數(shù)及其對(duì)應(yīng)的敏感波段。

表3 土壤重金屬砷含量與光譜指標(biāo)的最大相關(guān)系數(shù)Table 3 Maximum correlation coefficients of heavymetal As content and spectral indexes

可以發(fā)現(xiàn)， 經(jīng)CR光譜變換后682 nm與重金屬As含量的相關(guān)系數(shù)出現(xiàn)明顯峰值， 達(dá)到0.603(如圖4)， 這可能與重金屬As受到土壤有機(jī)質(zhì)的吸附有關(guān)。 FD的1 211 nm與土壤重金屬As含量存在最大負(fù)相關(guān)性。 其中通過(guò)0.01(雙尾)顯著性檢驗(yàn)的波段還有SD和ATSD的2 208 nm、 RTSD的1 065 nm、 ATFD的464 nm， SR的2 342 nm， 相關(guān)系數(shù)均超過(guò)0.370， 表示相關(guān)性顯著， 而RT和AT的1 740 nm相關(guān)系數(shù)大于0.28， 通過(guò)了0.05(雙尾)顯著性檢驗(yàn)， 表示相關(guān)性比較顯著。 RTFD未能通過(guò)任何相關(guān)性檢驗(yàn)， 因此在最后建立模型時(shí)不予考慮。

圖4 連續(xù)統(tǒng)去除相關(guān)系數(shù)曲線(xiàn)圖Fig.4 Continuum removal correlation coefficient curve

2.3 MLSR

MLSR是篩選對(duì)因變量影響較大的自變量建立回歸模型的方法。 根據(jù)自變量對(duì)因變量的作用大小， 剔除不顯著的自變量， 由大到小逐個(gè)引入回歸方程， 由于考慮了多個(gè)變量因素， 建模結(jié)果具有一定的可靠性。

根據(jù)土壤重金屬砷含量與光譜指標(biāo)的相關(guān)分析， 將表3中各光譜指標(biāo)對(duì)應(yīng)的敏感波段作為As含量估算模型的自變量， 以As實(shí)測(cè)值為因變量， 建立MLSR模型。

通過(guò)建模分析， 將FD1 211， RTSD1 065， ATFD464和CR6824個(gè)對(duì)模型擬合度影響較大的光譜指標(biāo)保留。 具體模型為

Y=-638.304-6 629.954XFD+1 110.784XRTSD-

391.284XATFD+651.907XCR

(5)

MLSR模型的散點(diǎn)圖和驗(yàn)證結(jié)果如表4和圖5所示， 建模樣本r超過(guò)0.8， 基本分布在1∶1線(xiàn)兩側(cè)， 說(shuō)明該模型對(duì)建模數(shù)據(jù)具有良好的解釋能力。 而驗(yàn)證樣本數(shù)據(jù)分布相較于建模樣本距離1∶1線(xiàn)存在較大偏離， 說(shuō)明模型預(yù)測(cè)效果不如建模效果， 模型穩(wěn)定度需要進(jìn)一步探究。

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

圖5 MLSR模型建模樣本和驗(yàn)證樣本散點(diǎn)圖(a)： 建模樣本； (b)： 驗(yàn)證樣本Fig.5 Scatter plots of fitting and testing sets for MLSR(a)： Fitting set； (b)： Testing set

2.4 U-PLSR

依照2.3分析， 將篩選出來(lái)模型精度影響較大的光譜指標(biāo)(FD1 211， RTSD1 065， ATFD464和CR682)分別作為自變量， 構(gòu)建U-PLSR模型。 為避免發(fā)生過(guò)度擬合現(xiàn)象， 采用逐一剔除法進(jìn)行交叉驗(yàn)證(下同)。 U-PLSR模型結(jié)果見(jiàn)表5。

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

依照模型的r和F值最大， RMSE最小原則， 經(jīng)CR光譜變換的建模和驗(yàn)證效果最好， 其次是FD， ATFD效果最差。 相比MLSR建模結(jié)果， 4種U-PLSR模型的r均較低， 最高僅達(dá)到0.381， 這表明U-PLSR模型穩(wěn)定性較差， 不能做到準(zhǔn)確有效的估算重金屬As含量。

2.5 M-PLSR

基于MLSR和U-PLSR模型， 為了尋找更為穩(wěn)定的模型， 進(jìn)一步建立M-PLSR模型， 所用模型光譜指標(biāo)和MLSR模型相同。

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

圖6 M-PLSR模型建模樣本和驗(yàn)證樣本散點(diǎn)圖(a)： 建模樣本； (b)： 驗(yàn)證樣本Fig.6 Scatter plots of the two set of samplesof fitting and testing for M-PLSR(a)：Fitting samples； (b)： Testing samples

表6 M-PLSR模型結(jié)果Table 6 Results of M-PLSR model

Y=-690.390-6 011.343XFD+1 174.766XRTSD-391.829XATFD+643.903XCR

(6)

建模樣本和預(yù)測(cè)樣本均超過(guò)0.8， 分布更為均勻， 相較于MLSR模型更靠近1∶1線(xiàn)。

M-PLSR模型的建模樣本r，F(xiàn)值和RMSE分別達(dá)到0.852， 32.384和0.147， 驗(yàn)證樣本r為0.812，F(xiàn)值為8.598， RMSE為0.312， 比MLSR模型相比， 建模r增加0.038，F(xiàn)值增加2.947， RMSE降低0.012， M-PLSR建模效果好于MLSR， 預(yù)測(cè)模型也有所提升。 綜合來(lái)看， M-PLSR模型擬合度更高， 建模效果更佳， 模型穩(wěn)定度也更強(qiáng)， 更適合用于估算研究區(qū)重金屬砷含量。

3 結(jié) 論

通過(guò)實(shí)地采樣分析， 研究區(qū)重金屬As含量平均值為13.230 mg·kg-1， 接近河北省背景值。 采礦和冶煉企業(yè)對(duì)土壤重金屬As含量有一定影響， 使得部分土壤樣本As含量已經(jīng)存在輕度污染情況， 大部分樣本略低于背景值， 處于污染的警戒水平。

研究區(qū)域和土壤背景的不同使得重金屬As的特征波段存在差別。 某礦區(qū)農(nóng)田潴育性水稻土重金屬As的特征波段有450， 700， 1 200和2 200 nm等[13]； 三江源地區(qū)以高山草甸土和高山草原土為主的特征波段為400， 429， 470， 550， 600， 815， 900， 1 415， 1 915， 2 205， 2 355和2 455 nm等[12]； 商洛黃褐土的特征波段包括450， 470， 640～660， 1 030～1 060， 1 080～1 110， 1 230和2 190～2 220 nm等[11]。 這些特征波段與本研究所得結(jié)果基本一致(464， 682， 1 065， 1 211和2 208 nm)， 但并不完全相同， 其中本研究與商洛黃褐土重合率最高， 可能因?yàn)閮煞N土壤類(lèi)型相似性較高。

土壤重金屬As含量的研究多采用單光譜變換指標(biāo)建模， 突出單一光譜指標(biāo)的特征， 模型可能存在誤差。 本工作選用了10種光譜指標(biāo)， CR與土壤重金屬As含量相關(guān)性最高， 且在U-PLSR模型中CR建模效果最好， 這主要是因?yàn)镃R能夠有效的消除研究區(qū)潛在因素對(duì)光譜的影響， 可以增強(qiáng)光譜對(duì)As的敏感性。 多光譜變換指標(biāo)集成建模結(jié)合多種最佳特征波段， 一定程度上避免了單一光譜指標(biāo)的不足， 達(dá)到提升建模效果和模型穩(wěn)定性的目的。

重金屬估算模型受到區(qū)域位置和土壤類(lèi)型甚至污染水平的影響， 構(gòu)建的最佳模型是否適用其他土壤類(lèi)型和其他地區(qū)土壤還有待深入研究。

針對(duì)石家莊市水源地保護(hù)區(qū)褐土重金屬As含量的高光譜估測(cè)， 通過(guò)相關(guān)分析， 建立了土壤重金屬砷含量MLSR， U-PLSR和M-PLSR估測(cè)模型， 具體結(jié)論如下：

(1)重金屬As與經(jīng)連續(xù)統(tǒng)去除變換后的光譜特征的相關(guān)性最大， 達(dá)到0.603， 與一階微分含量存在最大負(fù)相關(guān)性， 達(dá)到-0.527。

(2)相較于MLSR和U-PLSR， M-PLSR方法估算土壤重金屬As含量建模和驗(yàn)證樣本的r分別達(dá)到0.852和0.812， 模型預(yù)測(cè)值和重金屬As實(shí)測(cè)值整體擬合度最高。

(3)多光譜變換指標(biāo)建模效果顯著， 相比單光譜變換指標(biāo)模型， 它的建模結(jié)果和預(yù)測(cè)能力表現(xiàn)更好， 可以對(duì)估算重金屬As含量起到較好的作用。