李志遠， 鄧 帆， 賀軍亮， 魏 薇

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

引 言

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

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

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

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

1 實驗部分

1.1 研究區(qū)概況

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

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

1.2 樣本采集

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

表1 土壤中砷含量的描述性統計值(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]， 說明采集的部分土壤樣本中重金屬As含量已經富集到了一定程度。 土壤重金屬污染程度采用單因子評價法[8]來分析， 計算公式為

Pi=Ci/Si

(1)

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

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

表2 土壤中砷的單因子污染指數統計結果Table 2 Statistical results of single factorpollution index of the soil arsenic

1.3 光譜測定

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

1.4 光譜數據預處理

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

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

1.5 多變量反演方法

(1) MLSR

多元線性逐步回歸是通過一個因變量和多個自變量組合表達它們之間的線性關系。 回歸方程為

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

(2)

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

(2) PLSR

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

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

X=ThPT+E

(3)

Y=UhQT+F

(4)

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

1.6 模型建立與檢驗

隨機將48個樣本分為兩部分， 32個用作建立模型， 16個用來驗證模型。 采用相關系數r、 均方根誤差(root mean square error,RMSE)和統計值F來衡量建模效果。F值和r越大， RMSE越小， 建模樣本和驗證樣本的散點分布越靠近1∶1線， 建模效果越好， 由此確定最佳回歸模型。

2 結果與討論

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

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

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

2.2 特征波段提取

為了探究不同光譜指標與重金屬砷含量的關系， 將As含量與各種光譜變換進行Pearson相關性分析， 表3表示最大相關系數及其對應的敏感波段。

表3 土壤重金屬砷含量與光譜指標的最大相關系數Table 3 Maximum correlation coefficients of heavymetal As content and spectral indexes

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

圖4 連續(xù)統去除相關系數曲線圖Fig.4 Continuum removal correlation coefficient curve

2.3 MLSR

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

根據土壤重金屬砷含量與光譜指標的相關分析， 將表3中各光譜指標對應的敏感波段作為As含量估算模型的自變量， 以As實測值為因變量， 建立MLSR模型。

通過建模分析， 將FD1 211， RTSD1 065， ATFD464和CR6824個對模型擬合度影響較大的光譜指標保留。 具體模型為

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

391.284XATFD+651.907XCR

(5)

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

表4 MLSR模型結果Table 4 Results of MLSR model

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

2.4 U-PLSR

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

表5 U-PLSR模型結果Table 5 Results of U-PLSR model

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

2.5 M-PLSR

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

M-PLSR模型的散點圖和驗證結果如表6和圖6所示， 具體模型為

圖6 M-PLSR模型建模樣本和驗證樣本散點圖(a)： 建模樣本； (b)： 驗證樣本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模型結果Table 6 Results of M-PLSR model

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

(6)

建模樣本和預測樣本均超過0.8， 分布更為均勻， 相較于MLSR模型更靠近1∶1線。

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

3 結 論

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

研究區(qū)域和土壤背景的不同使得重金屬As的特征波段存在差別。 某礦區(qū)農田潴育性水稻土重金屬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]。 這些特征波段與本研究所得結果基本一致(464， 682， 1 065， 1 211和2 208 nm)， 但并不完全相同， 其中本研究與商洛黃褐土重合率最高， 可能因為兩種土壤類型相似性較高。

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

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

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

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

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

(3)多光譜變換指標建模效果顯著， 相比單光譜變換指標模型， 它的建模結果和預測能力表現更好， 可以對估算重金屬As含量起到較好的作用。