張家港市農田土壤重金屬含量高光譜遙感監(jiān)測模型構建

錢家煒，劉曉青，張靜靜，周衛(wèi)紅，2，李建龍，*

(1.南京大學 生命科學學院 應用生態(tài)研究所，江蘇 南京 210093； 2.江蘇科技大學 蘇州理工學院，江蘇 張家港 215600)

近年來，隨著經(jīng)濟和城市化的發(fā)展，重金屬污染問題日益嚴重。重金屬污染會嚴重影響土壤的物理、化學性質和生物學特性，還會抑制生物酶的活性[1]，并通過土壤-作物系統(tǒng)危害農產(chǎn)品質量安全與人體健康[2]。重金屬含量預測是預防、評估和治理重金屬污染的關鍵環(huán)節(jié)。目前，常規(guī)的重金屬含量預測技術多是通過采集樣本和實驗室化學分析，借助于ArcGIS軟件進行空間插值完成的。當采集的樣本數(shù)量有限時，其結果并不能代表整個研究區(qū)；但大規(guī)模采集樣本又不切實際，耗時，且費用昂貴，而且實驗室使用強酸等物質進行分析時也存在安全隱患[3]。遙感作為一種新型的探測手段，可以方便、快捷、動態(tài)地獲得空間上連續(xù)的地物光譜信息[4]，其中，可見光和近紅外反射光譜已被廣泛用于測量土壤的物理和化學性質[5]。

國內外許多學者利用高光譜遙感監(jiān)測對土壤重金屬含量進行研究，并取得了很好的成果。Ma等[6]利用超限學習機、支持向量機等方法建模，估測土壤中5種重金屬的含量。劉華等[7]根據(jù)EO-1(Earth Observing-1，地球觀察者1號)衛(wèi)星Hyperion高光譜儀的波段設置，結合實測數(shù)據(jù)，基于光譜反射率運用偏最小二乘法很好地估測了土壤重金屬含量。袁中強等[8]分析了環(huán)境一號衛(wèi)星上搭載的超光譜成像儀的地表反射率數(shù)據(jù)與實測的土壤重金屬元素含量的關系，建立了預測土壤中5種重金屬含量的原始反射率、一階導數(shù)、倒數(shù)對數(shù)3個回歸模型。St Luce等[9]利用偏最小二乘法建立可見近紅外光譜與土壤中總鋅、總鎘的模型。Choe等[10]研究發(fā)現(xiàn)，可見和近紅外范圍內610 nm和500 nm光譜反射率的比值、2 200 nm處的吸收面積，以及2 200 nm處吸收特征的不對稱性分別與Pb、Zn、As濃度存在顯著相關性。綜上所述，利用高光譜與土壤重金屬含量建立的反演模型可以預測重金屬含量。但前人的研究多是用不同變換方式的光譜反射率建立多個模型，然后通過比較精度篩選出最優(yōu)模型。與此不同，本文擬在建立模型的過程中，分析5種變換方式與8種重金屬含量的相關性，先篩選出每種重金屬對應的最優(yōu)變換方式，再建立模型。為此，特以長三角地區(qū)張家港市為研究對象，選取對人體有傷害、且存在范圍比較廣的8種重金屬做總結性研究，在室內用ASD地物光譜儀(美國Analytica Spectra Devices公司)獲取土壤高光譜數(shù)據(jù)，與經(jīng)電感耦合等離子體原子發(fā)射光譜(ICP-AES)法測量的土壤As、Cd、Cr、Cu、Zn、Ni、Pb、Hg含量相結合，光譜經(jīng)過各種形式的變換，根據(jù)相關性分析選出每種重金屬對應的敏感特征波段，利用逐步回歸法建立土壤重金屬的高光譜反演模型，預測土壤重金屬含量。

1 材料與方法

1.1 研究區(qū)概況

張家港市位于江蘇省蘇州市北部，長江下游南岸，是沿海開放地區(qū)和長三角經(jīng)濟帶交匯區(qū)較早開放的工業(yè)城市，地處120°21′～120°52′E、31°43′～32°02′N，年平均氣溫17 ℃，年降水量1 200 mm以上，年日照時數(shù)在1 600 h以上，屬北亞熱帶南部濕潤性氣候。張家港市現(xiàn)有耕地面積30 646.7 hm2，常年種植水稻18 133.3 hm2、小麥19 000.0 hm2、油菜1 600.0 hm2、蔬菜3 666.7 hm2、果品1 400 hm2。張家港市下轄8鎮(zhèn)2區(qū)，分別為楊舍鎮(zhèn)、塘橋鎮(zhèn)、金港鎮(zhèn)、錦豐鎮(zhèn)、樂余鎮(zhèn)、鳳凰鎮(zhèn)、南豐鎮(zhèn)、大新鎮(zhèn)，及常陰沙現(xiàn)代農業(yè)示范園區(qū)、雙山島旅游度假區(qū)。

圖1 采樣點分布圖Fig.1 Distribution of soil samples

1.2 樣品采集與制備

在地理空間數(shù)據(jù)云上下載張家港市的TM數(shù)據(jù)，并將數(shù)據(jù)在ENVI 5.3軟件中進行監(jiān)督分類，提取農田土壤，以混淆矩陣法算出分類結果精度為0.87。在張家港市8個鎮(zhèn)和1個現(xiàn)代農業(yè)示范園區(qū)的農田土壤中均勻地設置23個采樣點。依據(jù)樣點分布(圖1)，用GPS精準定位后去野外采集0～20 cm的農田土壤樣品。將收集的土壤樣本帶回實驗室自然風干，粗略去除砂礫和植物殘體，研磨、過100目塑料尼龍篩后混合均勻。采用四分法取樣，一式2份，一份用于測定重金屬含量，另一份用于光譜測定。土壤As、Cd、Cr、Cu、Zn、Ni、Pb用硝酸-鹽酸-高氯酸消解，Hg經(jīng)硝酸-鹽酸混合試劑在沸水浴中加熱消解，消解后均采用ICP-AES法[11]進行測定。

1.3 光譜測定

在室內條件下，采用ASD光譜儀測定經(jīng)過處理的土壤樣品的光譜反射率，光譜波長范圍在350～2 500 nm。將土壤盛裝在直徑10 cm、深3 cm的黑色器皿中。在進行光譜測定之前，先對土壤表面做刮平處理[12]。光譜測定在暗室中進行，以功率50 W的鹵素燈作為唯一光源，測量方法嚴格按照有關規(guī)范操作：光源入射角為45°，光源距離土樣表面中心30 cm，探頭距離土樣15 cm。設置40 cm×40 cm的參考白板，用于獲取絕對反射率。每個土樣測量10次光譜曲線，取平均值，以減少誤差，提高精度。利用View Spec Pro軟件剔除異常曲線后，取光譜反射率平均值作為初始反射率光譜值[13]。

1.4 光譜變換

在對目標進行光譜測定的過程中，背景、大氣、儀器、光照條件等因素影響都可能覆蓋目標自身的光譜信息特征。因此，嘗試進行多種形式的光譜變換，以減弱背景噪聲的影響[4]，增強目標光譜信息，提高信噪比。為描述方便，本文將經(jīng)過平滑處理后的反射光譜稱為原始反射光譜。對原始反射光譜分別進行一階導數(shù)、倒數(shù)一階導數(shù)、倒數(shù)的對數(shù)一階導數(shù)、平方根一階導數(shù)和連續(xù)統(tǒng)去除[14]共5種形式的變換處理。

1.5 模型建立與精度分析

1.5.1 相關性分析

在高光譜數(shù)據(jù)與農田土壤重金屬含量建模的過程中，篩選與重金屬含量相關性高的波段作為特征敏感波段，相關性越高，波段響應越敏感[15]。將經(jīng)過變換后的光譜反射率與土壤重金屬含量進行皮爾遜(Pearson)相關分析，考慮不同重金屬光譜反演的需要，根據(jù)每種重金屬的光譜特征，篩選出相關系數(shù)通過P<0.05水平顯著性檢驗的波段，作為建立高光譜定量估算模型的自變量。

1.5.2 建模方法

逐步回歸的基本思想簡述如下：(1)將變量逐個引入模型；(2)每引入一個解釋變量后都要進行F檢驗，并對已選入的解釋變量逐個進行t檢驗；(3)當原變量因為新解釋變量的引入變得不再顯著時，將其剔除，以確保每次引入新的變量之前回歸方程中只包含主動變量；(4)重復進行這一過程，直到既沒有顯著的解釋變量選入回歸方程，也沒有不顯著的解釋變量從回歸方程中剔除為止，以保證最后所得到的解釋變量集最優(yōu)[16]。

依據(jù)上述思想，利用逐步回歸篩選并剔除引起多重共線性的變量[17]，步驟如下：先用被解釋變量對每一個所考慮的解釋變量做簡單回歸，然后以對被解釋變量貢獻最大的解釋變量所對應的回歸方程為基礎，再逐步引入其余解釋變量[18]。經(jīng)過逐步回歸，使得最后保留在模型中的解釋變量既是重要的，同時又可減小多重共線性[19]。

2 結果與分析

2.1 土壤樣品統(tǒng)計分析

土壤重金屬含量測定結果見表1。從標準差來看，Hg、Cd、Cr的標準差都小于1 mg·kg-1，其余5種重金屬元素的標準差為1.173～5.568 mg·kg-1。從變異系數(shù)來看，Cd、Hg的變異系數(shù)分別為0.342、0.315，其余的變異系數(shù)均在0.1左右，變異程度為中等變異性[20]。從平均值來看，Cd、Hg均超過背景值，其余重金屬元素的平均含量則均低于背景值，說明該區(qū)域土壤環(huán)境質量總體良好；但從最大值看，有部分樣點的土壤重金屬含量與背景值相近，部分樣點的Cd、Hg、Cu、Zn含量已超出背景值。根據(jù)土壤重金屬污染的單因子指數(shù)法測算，As、Cr、Pb、Ni的污染指數(shù)均小于0.7，含量在安全范圍內；Zn、Cu的污染指數(shù)接近1，存在潛在危害；Cd的污染指數(shù)為1.324，屬于輕度污染；Hg的污染指數(shù)為2.724，屬于中度污染。因此，張家港市應該加強土壤質量調查與動態(tài)監(jiān)測，以便及時發(fā)現(xiàn)并控制耕地土壤的重金屬污染問題。

表1 土壤重金屬含量統(tǒng)計特征

2.2 光譜變換分析

如圖2所示，將原始光譜進行一階導數(shù)、連續(xù)統(tǒng)去除、平方根一階導數(shù)、倒數(shù)一階導數(shù)、倒數(shù)的對數(shù)一階導數(shù)變換后發(fā)現(xiàn)，除了倒數(shù)一階導數(shù)減少了原來的光譜特征信息外，一階導數(shù)、連續(xù)統(tǒng)去除、平方根一階導數(shù)、倒數(shù)的對數(shù)一階導數(shù)變換均能夠將原始光譜波峰、波谷等信號特征放大，產(chǎn)生更大幅度的波動，增大與背景光譜的差異。其中，一階導數(shù)和連續(xù)統(tǒng)去除變換方法都在相同的幾個波段處變化劇烈，放大了原始光譜信息，在消除或減弱土壤背景噪聲、提高信噪比、增強目標光譜信息方面作用明顯，更容易從光譜數(shù)據(jù)中提取出能夠反映土壤特性的微弱信號。

a，原始光譜；b，連續(xù)統(tǒng)去除；c，一階導數(shù)；d，平方根一階導數(shù)；e，倒數(shù)一階導數(shù)；f，導數(shù)的對數(shù)一階導數(shù)。a， Raw spectrum; b， Continuum removal; c， First derivative; d， Square root of first derivative; e， Reciprocal of first derivative; f， Logarithm of reciprocal of first derivative.圖2 光譜變換結果Fig.2 Spectral transformation result

2.3 土壤重金屬含量與光譜反射率的相關性分析

分別計算每種重金屬元素含量和對應原始光譜的一階導數(shù)、連續(xù)統(tǒng)變換的相關系數(shù)，比較不同光譜變換形式對該相關程度的影響，找出特征波段。

將土壤重金屬含量與原始光譜的一階導數(shù)進行相關性分析，得到每種重金屬與對應光譜反射率的相關系數(shù)曲線，并對其進行顯著性檢驗(顯著性水平設定為P<0.05)。從圖3可以看出，As、Cd、Cr、Pb、Zn、Ni、Hg含量與光譜反射率相關系數(shù)曲線中的許多波段都通過了P<0.05水平的顯著性檢驗。其中，As、Cr、Pb、Zn、Ni、Hg的相關系數(shù)在900～2 500 nm通過顯著性檢驗的波段比較多；Cd相關系數(shù)通過顯著性檢驗的波段數(shù)少于除Cu以外的其他6種重金屬；Cu與光譜反射率的相關系數(shù)較低，未能通過P<0.05水平的顯著性檢驗。

從圖4可以看出，Cr、Zn、Hg含量與反射率連續(xù)統(tǒng)去除的相關系數(shù)通過顯著性檢驗的波段很少；Ni相關系數(shù)通過顯著性檢驗的波段集中在360～600 nm，且與反射率呈現(xiàn)負相關；As相關系數(shù)通過顯著性檢驗的波段集中在355～585、620～780 nm；Pb相關系數(shù)通過顯著性檢驗的波段集中在995～1 305、1 490～1 790 nm；Cd、Cu通過顯著性檢驗的波段較多，Cd相關系數(shù)通過顯著性檢驗的波段集中在360～800、850～935、1 800～1 840 nm，Cu相關系數(shù)通過顯著性檢驗的波段集中在860～1 340、1 450～1 660、2 180～2 485 nm。綜上所述，As、Cr、Zn、Pb、Ni、Hg含量與原始反射率一階導數(shù)的相關性高，Cd、Cu含量與原始反射率連續(xù)統(tǒng)去除的相關性較高。

2.4 高光譜模型建立與驗證

將23個樣本分成2組，從第1個樣本開始，每5個抽取1個樣本用于模型檢驗；因此，模型構建樣本數(shù)為18，模型檢驗樣本數(shù)為5。建立模型時，選擇與重金屬含量相關性高的光譜變換形式的反射率作為自變量(x)，以重金屬含量作為因變量(y)。在SPSS 20.0中進行逐步回歸分析，默認設置F值概率(P)小于0.05進入，大于0.1刪除，得到高光譜重金屬含量的定量估算模型。

圖3 反射率的一階導數(shù)與土壤重金屬含量的相關系數(shù)Fig.3 Correlation coefficient between first derivative of reflectance and soil heavy metal content

圖4 反射率的連續(xù)統(tǒng)去除與土壤重金屬含量的相關系數(shù)Fig.4 Correlation coefficient between continuum removal of reflectance and soil heavy metal content

基于最優(yōu)光譜形式組合和重金屬含量，用逐步回歸法建立模型。由表2可以看出，各種重金屬和光譜建立的線性模型中P值均小于0.01，說明所建模型可以反映光譜與重金屬含量間的關系。As、Cd、Cr和Hg定量估算模型的擬合度(R2)均達到0.7以上，其中Hg的估算模型的擬合度最高(R2=0.777)，Zn、Ni的擬合度達到0.6以上，Cu、Pb的擬合度在0.5～0.6。8種重金屬估算模型的F統(tǒng)計量均較大，除了Pb的估算模型F值為9.654以外，其他模型的F值都大于10，其中Cd的估算模型的F值最大，為27.710。綜上，建立的8種重金屬估算模型可以在一定程度上反映研究區(qū)的土壤重金屬含量。

如圖5所示，用之前的5個驗證樣本對所構建的模型進行驗證，得到模型實際值與驗證值的擬合曲線。從擬合曲線可以看出，Cd、Hg估算模型的實際值與驗證值的擬合度大于0.8，分別為0.874、0.879，Cr估算模型的擬合度為0.800，As、Cu、Zn、Ni、Pb估算模型的擬合度也都大于0.5，說明8種重金屬的定量估算模型驗證精度較高，具有預測張家港市農田土壤重金屬的能力，能進一步用于大面積遙感估算。

表2 所建立的光譜與重金屬含量模型

3 討論

基于本研究試驗數(shù)據(jù)，張家港市農田土壤重金屬As、Cd、Cr、Cu、Pb、Zn、Ni、Hg的含量在0.053～99.258 mg·kg-1，部分樣點的Cd、Hg、Cu、Zn含量超出背景值。依土壤重金屬污染的單因子指數(shù)法測算：Hg的污染指數(shù)為2.724，屬于中度污染；Cd的污染指數(shù)為1.324，屬于輕度污染；Zn、Cu的污染指數(shù)都接近1，存在潛在危害；As、Cr、Pb、Ni的污染指數(shù)小于0.7，重金屬含量在安全范圍內。

通過土壤重金屬與光譜反射率的相關性分析，選取As的敏感波段為2 069、2 201、2 267 nm，Cd的敏感波段為523、1 809 nm，Cr的敏感波段為2 065、2 413、1 634 nm，Cu的敏感波段為2 225、2 431 nm，Zn的敏感波段為2 332、2 060 nm，Ni的敏感波段為2 332、1 922 nm，Pb的敏感波段為2 332、2 038 nm，Hg的敏感波段為2 312、1 573、997 nm。選其用作張家港市土壤重金屬含量反演模型的自變量。

8種重金屬和光譜建立的線性模型，P值均小于0.05，說明構建的逐步回歸模型整體顯著。在土壤重金屬含量和高光譜波段建立的模型中，Cu、Pb的擬合優(yōu)度達到0.5以上，Zn、Ni的擬合優(yōu)度達到0.6以上，As、Cd、Cr、Hg的擬合優(yōu)度達到0.7以上，擬合優(yōu)度較高。從模型實際值與驗證值的擬合曲線得出，Cd、Hg實際值與驗證值的擬合度分別為0.874、0.879，As、Cr、Cu、Zn、Ni、Pb的擬合度也都大于0.5。由此可知，8種重金屬估算模型驗證精度良好。以上結果說明，對As、Cr、Zn、Pb、Ni和Hg光譜反射率進行一階求導處理，對Cd和Cu光譜反射率進行連續(xù)統(tǒng)去除處理后，所建立的8種重金屬光譜模型經(jīng)驗證效果理想，可有效預測土壤重金屬含量，是替代傳統(tǒng)重金屬預測方法的優(yōu)化選擇，可在實際生產(chǎn)中加以應用。

本研究采用多元線性逐步回歸方法對8種重金屬含量進行了建模分析，但由于只設置了23個采樣點，建模集和驗證集數(shù)量均有一定不足，因此，今后應進一步擴大采樣范圍以優(yōu)化模型性能。與已有的基于高光譜方法預測土壤重金屬含量的研究相比，李瓊瓊等[21]對土壤Cu、Pb、Zn含量進行反演研究的結果表明，最小二乘回歸方法的建模精度優(yōu)于多元線性逐步回歸；王金鳳等[22]運用隨機森林、支持向量機、偏最小二乘3種算法對土壤Zn含量進行建模，發(fā)現(xiàn)基于二階微分變換光譜的隨機森林算法準確度最高；許吉仁等[23]用支持向量機方法對土壤Cd含量建立高光譜監(jiān)測模型，R2為0.947。為了探究高光譜方法預測土壤重金屬含量的更優(yōu)方法，未來可考慮引入多元散射校正、標準正態(tài)化等光譜預處理方法和偏最小二乘、隨機森林的建模方法來進行比選優(yōu)化。