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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

張家港市位于江蘇省蘇州市北部，長江下游南岸，是沿海開放地區(qū)和長三角經(jīng)濟(jì)帶交匯區(qū)較早開放的工業(yè)城市，地處120°21′～120°52′E、31°43′～32°02′N，年平均氣溫17 ℃，年降水量1 200 mm以上，年日照時(shí)數(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)代農(nóng)業(yè)示范園區(qū)、雙山島旅游度假區(qū)。

圖1 采樣點(diǎn)分布圖Fig.1 Distribution of soil samples

1.2 樣品采集與制備

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

1.3 光譜測(cè)定

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

1.4 光譜變換

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

1.5 模型建立與精度分析

1.5.1 相關(guān)性分析

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

1.5.2 建模方法

逐步回歸的基本思想簡(jiǎn)述如下：(1)將變量逐個(gè)引入模型；(2)每引入一個(gè)解釋變量后都要進(jìn)行F檢驗(yàn)，并對(duì)已選入的解釋變量逐個(gè)進(jìn)行t檢驗(yàn)；(3)當(dāng)原變量因?yàn)樾陆忉屪兞康囊胱兊貌辉亠@著時(shí)，將其剔除，以確保每次引入新的變量之前回歸方程中只包含主動(dòng)變量；(4)重復(fù)進(jìn)行這一過程，直到既沒有顯著的解釋變量選入回歸方程，也沒有不顯著的解釋變量從回歸方程中剔除為止，以保證最后所得到的解釋變量集最優(yōu)[16]。

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

2 結(jié)果與分析

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

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

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

2.2 光譜變換分析

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

a，原始光譜；b，連續(xù)統(tǒng)去除；c，一階導(dǎo)數(shù)；d，平方根一階導(dǎo)數(shù)；e，倒數(shù)一階導(dǎo)數(shù)；f，導(dǎo)數(shù)的對(duì)數(shù)一階導(dǎo)數(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 光譜變換結(jié)果Fig.2 Spectral transformation result

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

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

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

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

2.4 高光譜模型建立與驗(yàn)證

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

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

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

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

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

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

3 討論

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

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

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