北京市平谷區(qū)桃林土壤重金屬元素含量反演估測

劉泓君，牛 騰，于 強(qiáng)*, 蘇 凱，楊林哲，劉 維，王慧媛

1. 北京林業(yè)大學(xué)林學(xué)院，北京 100083 2. 廣西大學(xué)林學(xué)院，廣西 南寧 530005

引 言

近些年來，礦產(chǎn)資源的開采及工業(yè)污水的排放對礦區(qū)周邊自然環(huán)境的可持續(xù)發(fā)展和人體健康造成了嚴(yán)重威脅。特別在種植經(jīng)濟(jì)作物的地區(qū)，當(dāng)重金屬含量嚴(yán)重超標(biāo)后進(jìn)入種植物體內(nèi)，進(jìn)而造成植物光合作用減弱、 葉片呈現(xiàn)棕色斑塊等癥狀。人類直接或間接食用此類作物果實后，易引起重金屬中毒，重則危害生命[1-2]。北京東北部平谷區(qū)，栽培桃樹面積世界第一，但其礦石的冶煉及采選造成重金屬物質(zhì)大量殘留在土壤中，對周圍桃林造成了嚴(yán)重的污染[3]，因此了解平谷區(qū)重金屬污染情況極為重要。傳統(tǒng)以“點采樣+實驗室分析”、 電感耦合等離子體發(fā)射光譜法等方法分析重金屬含量，但此類方法依賴于大量的外業(yè)數(shù)據(jù)，費時費力，且不具有推廣性。多光譜遙感具有高分辨率、 高精度、 快速無損、 大面積監(jiān)測的優(yōu)勢，為宏觀獲取植被覆蓋區(qū)土壤重金屬元素污染情況提供了新的途徑[4]。

目前，基于遙感技術(shù)監(jiān)測土壤重金屬污染取得了豐碩的成果。Anne等基于高光譜影像，提取紅樹林植被反射光譜，構(gòu)建差值、 比值植被指數(shù)，基于偏最小二乘模型，成功反演出土壤中的有機(jī)質(zhì)、 黏土礦物等含量[5]。高偉等以玉米為研究對象，通過運用包絡(luò)線去除處理葉片光譜，構(gòu)建作物銅鉛污染判別規(guī)則[6]；楊靈玉等基于高光譜影像植被光譜反射率間接估算出土壤Zn和Cd元素含量[7]。綜上所述，前人研究多集中在運用幾個敏感的光譜波段或者常用植被指數(shù)建立反演模型，而對于基于敏感波段構(gòu)建植被指數(shù)來進(jìn)行反演的研究方法較少，且以桃葉為研究對象的更少。

以北京平谷區(qū)為例，將桃林作為監(jiān)測目標(biāo)，以野外采集的土壤樣本、 桃葉光譜數(shù)據(jù)、 Sential-2多光譜影像為基礎(chǔ)，探究受重金屬脅迫影響較為敏感的特征波段并構(gòu)建植被指數(shù)，同時將其與紅色歸一化植被指數(shù)、 類胡蘿卜素反射指數(shù)2、 結(jié)構(gòu)不敏感色素指數(shù)、 歸一化水指數(shù)四種常用植被指數(shù)做對比分析，驗證該指數(shù)在監(jiān)測重金屬方面的有效性與優(yōu)越性；通過植被指數(shù)與土壤重金屬Cd，AS和Pb的擬合建立反演模型，最終獲取果園重金屬污染空間分布情況，為快速有效地監(jiān)測經(jīng)濟(jì)作物區(qū)域土壤生態(tài)狀況提供技術(shù)支持。

1 實驗部分

1.1 研究區(qū)概況

北京市平谷區(qū)(116°55′—117°24′E，40°1′—40°22′N)，西距北京市區(qū)70 km，總面積1 075 km2，約有40萬人。屬于溫帶大陸性季風(fēng)氣候，年降水量為629.4 mm，降雨主要集中在6月—9月，僅28%的降雨量出現(xiàn)在其他月份, 平均最高氣溫為17.3 ℃。平谷是中國大桃第一區(qū)，大桃種植面積達(dá)16.8萬畝，產(chǎn)量1.5億公斤，其品種達(dá)到130多個。平谷區(qū)擁有大量豐富的金屬礦和非金屬礦等礦產(chǎn)資源，礦床、 礦點、 礦化點共79處。該地區(qū)礦區(qū)開采過程中，礦石礦渣中的重金屬經(jīng)雨水沖刷、 淋溶滲入周邊土壤中或河流中，對周邊生態(tài)環(huán)境造成污染。

圖1 研究區(qū)概況及采樣點分布位置Fig.1 Overview of the study area and distributionlocation of sampling points

1.2 樣本采集

依據(jù)當(dāng)?shù)氐牡乩憝h(huán)境特征，并遵循均勻性、 代表性、 規(guī)律性的原則來布置采樣點，共87個采樣點的空間分布如圖1所示。選擇無風(fēng)晴朗天氣，于2021年7月4日至25日10：00—14：00時間段，在每個布設(shè)點選取樹齡為9年的京艷桃樹(2株)采集相同數(shù)量的桃葉(共20片)、 使用GPS定位，詳細(xì)記錄采樣點編號、 地理坐標(biāo)和其所處空間類型等信息。同時，在礦區(qū)周圍每個布設(shè)點處選取5～10個分樣點采集土壤，采樣深度為0～40 cm，等量多點采集后將其均勻混合在一起，從中抽取1 kg混合樣品封裝，編號帶回實驗室。

1.3 遙感數(shù)據(jù)及處理

遙感數(shù)據(jù)來自歐空局哥白尼數(shù)據(jù)中心(https://scihub.copernicus.eu/dhus/#/home)，下載2021年七月份云量在10%以下的Sential-2影像。Sential-2具有13個光譜波段，在紅邊范圍內(nèi)含有三個波段的數(shù)據(jù)，有益于監(jiān)測指數(shù)健康信息。地面分辨率分別為10，20和60 m。利用ENVI對影像進(jìn)行輻射定標(biāo)、 大氣校正、 幾何校正、 拼接、 裁剪、 重采樣，預(yù)處理后得到研究區(qū)地表反射率數(shù)據(jù)。

1.4 葉片光譜反射率與重金屬含量測定

葉片光譜的采集使用ASD FieldSpec 4 便攜式地物光譜儀，波長范圍為350～2 500 nm，光譜分辨率為3 nm，光譜帶寬為2.2 nm，采樣時間為200 ms。測定前以白色參考板對光譜儀進(jìn)行標(biāo)準(zhǔn)化，測量過程中，選擇50 W的鹵素?zé)魹楣庠?，光源照射方向與樣品呈15°夾角，葉子置于水平臺距光源30 cm，距探頭5 cm處。采集每個樣點20片桃葉，剔除350～399 nm信噪比低、 噪聲大的邊緣波段，最后取平均反射率作為該樣點桃葉實際光譜反射數(shù)據(jù)。

將采集的土壤樣品碾碎、 風(fēng)干，過10目尼龍篩去除雜質(zhì)后，置于60 ℃烘箱烘干后，再次過16目尼龍篩，密封裝袋后用于實驗室化學(xué)分析。土壤平均pH值為5.78，平均有機(jī)質(zhì)含量為17.62 g·kg-1，黏粒含量為25.6%。其中，土壤As，Cr和Zn采用電感耦合等離子體質(zhì)譜儀進(jìn)行測試，使用xSPECTRAA-220型原子吸收光譜儀測量土壤Cd，Pb和Cu含量。采用陳同斌等[8]提出的北京市土壤重金屬的背景值作為本研究的背景值。如表1所示，平谷區(qū)土壤中重金屬平均含量與背景值比值大小依次是Cr>As>Cu>Cd>Pb>Zn。根據(jù)王威等[9]對平谷區(qū)重金屬危害程度的研究，綜合考慮比值與危害程度的大小，選取重金屬Cd，AS和Pb元素作為反演對象。

表1 土壤重金屬含量描述性統(tǒng)計Table 1 Descriptive statistics of heavy metal content in soil

1.5 方法

1.5.1 特征變量選取

當(dāng)植物中的重金屬濃度非常低時，反演難度大[10-11]。單一的影像光譜反射率很難突出其差異，為增強(qiáng)光譜響應(yīng)特征波段，采用一階微分(first deribative，F(xiàn)D)、 二階微分(second derivative，SD)、 倒數(shù)對數(shù)(reciprocal logarithm，RL)、 連續(xù)統(tǒng)去除法(continum removal，CR)等變換對光譜數(shù)據(jù)進(jìn)行預(yù)處理，增強(qiáng)有效波譜信息。為確定特征波段，采用相關(guān)性分析，選取在0.01水平下極顯著相關(guān)的波段作為重金屬特征波段。

1.5.2 基于光譜特征的植被指數(shù)計算

由于重金屬污染對植被光譜特征有一定的影響，植被指數(shù)能較好的反應(yīng)生長狀況。依據(jù)前人研究，重金屬對植被的脅迫主要體現(xiàn)在影響葉綠素的合成和阻礙植被中水分輸送。固選擇紅色歸一化植被指數(shù)(NDVI705)、 類胡蘿卜素反射指數(shù)2(CRL2)、 結(jié)構(gòu)不敏感色素指數(shù)(SIPI)、 歸一化水指數(shù)(NDWI)以及基于特征波段構(gòu)建的植被指數(shù)構(gòu)建模型反演葉片的重金屬含量。

2 結(jié)果與討論

2.1 桃葉光譜對Cd元素的響應(yīng)

將采集的葉片光譜曲線按照重金屬含量分為兩組，土壤中Cd，As和Pb含量小于背景值的組，代表正常葉片；其余樣本的光譜曲線，代表受重金屬脅迫的葉片。將兩組光譜曲線分別取其平均值，圖2顯示了正常桃葉與受重金屬脅迫桃葉光譜反射率曲線的差異。桃葉總體變化趨勢基本一致，在450～750 nm范圍內(nèi)有兩個明顯的光譜吸收帶和反射波峰與葉綠素吸收藍(lán)光、 紅光，反射綠光相關(guān)，由于植物吸收水分在1 450和1 900 nm形成兩個吸收帶，具有典型植被光譜曲線的特征。但受重金屬脅迫葉片的平均光譜反射率高于正常葉片，且740～1 300 nm波長范圍內(nèi)差值明顯，其原因是植被遭重金屬脅迫時，葉綠素大量減少，葉黃素、 葉紅素相對增加，導(dǎo)致綠光波段反射率變高。

圖2 受重金屬脅迫葉片與正常葉片反射光譜曲線(a) 及其一階微分(b)

植被的藍(lán)邊、 黃邊和紅邊位置可以反映植被的生長態(tài)勢和健康狀況，紅邊位置可以作為植物受污染后的重要指示波段。提取680～750 nm范圍內(nèi)反射率一階導(dǎo)數(shù)的最大值作為紅邊的斜率，其最大值處所對應(yīng)的波長代表紅邊位置，該波段范圍內(nèi)一階導(dǎo)數(shù)所包圍的面積作為紅邊面積。同理，用相同的方法確定藍(lán)邊(490～530 nm)及黃邊(550～580 nm)的位置。由表2顯示，受重金屬脅迫的葉片相較于正常葉片紅邊位置發(fā)生了“藍(lán)移”現(xiàn)象，藍(lán)邊和黃邊位置變化不大，表現(xiàn)出了較強(qiáng)的抗干擾能力。紅邊斜率及紅邊面積隨重金屬污染程度加劇而減小。

表2 藍(lán)邊、 黃邊、 紅邊參數(shù)Table 2 Parameter of blue edge, yellow edge and red edge

2.2 桃葉光譜特征波段的選取

為了更好增強(qiáng)信噪比、 提高分辨率、 突出葉片的光譜特性，提高葉片光譜數(shù)據(jù)與土壤各參數(shù)間的相關(guān)性，將原始光譜進(jìn)行一階微分、 二階微分、 倒數(shù)對數(shù)、 連續(xù)統(tǒng)去除法等變換預(yù)處理，將不同預(yù)處理后的光譜分別與三種重金屬含量進(jìn)行Pearson相關(guān)性分析。圖3依次表示為一階微分、 二階微分、 倒數(shù)對數(shù)、 連續(xù)統(tǒng)去除法與三種重金屬含量的相關(guān)性曲線。結(jié)果表明，預(yù)處理后的光譜與三種重金屬含量間的相關(guān)系數(shù)曲線均有明顯峰值，而二階微分處理后光譜與三種重金屬含量相關(guān)系數(shù)曲線峰值分布較多且分散。選取每種變換中六個相關(guān)系數(shù)較高的峰值波段，相關(guān)系數(shù)如表3所示。

圖3 四種變換形式的光譜數(shù)據(jù)Fig.3 Pearson correlations of spectral data using four transformations

表3 三種重金屬含量與光譜變換特征吸收波段相關(guān)系數(shù)Table 3 Correlation coefficients between the contents of three heavy metals andthe characteristic absorption bands of four spectral transformation

相關(guān)系數(shù)篩選出的光譜特征波段雖與重金屬含量有相關(guān)性，但還需要擬合回歸方程確定對重金屬含量預(yù)測起重要作用的光譜特征波段，通過剔除回歸建模中貢獻(xiàn)率不顯著的波普，得到貢獻(xiàn)率顯著的波段，如表4所示。

將逐步回歸篩選的重復(fù)特征波段進(jìn)行有效整合，最終確定三種重金屬的特征光譜變量為：780，945和1 375 nm。

2.3 植被指數(shù)的構(gòu)建與篩選

植被指數(shù)可以很好反應(yīng)植被物理參數(shù)，監(jiān)測植被生長狀況。為了充分利用780，945和1 375 nm三個特征波段的信息，便于直接分析特征波段與重金屬含量，構(gòu)建重金屬脅迫植被指數(shù)(HMSVI)

HMSVI=(R1 375-R945)/(R1 375-R780)

(1)

并選取應(yīng)用于植被脅迫性探測的指數(shù)：紅色歸一化植被指數(shù)(NDVI705)、 類胡蘿卜素反射指數(shù)2(CRL2)、 結(jié)構(gòu)不敏感色素指數(shù)(SIPI)、 歸一化水指數(shù)(NDWI)。將以上五種植被指數(shù)分別與三種不同重金屬含量做相關(guān)分析，相關(guān)系數(shù)如表5所示。NDVI705，SIPI和HMSVI與三種重金呈極顯著正相關(guān)性；NDWI與三種重金屬呈負(fù)相關(guān)性。四種植被指數(shù)均能較好反應(yīng)葉片重金屬污染程度，但NDVI705及SIPI與重金屬的相關(guān)性較弱，CRL2與三種重金屬都沒有顯著相關(guān)性；而HMSVI與三種重金屬相關(guān)性系數(shù)均較高，固HMSVI用于三種重金屬建模預(yù)測效果最佳。

2.4 重金屬Cd，As和Pb反演模型

以植被指數(shù)為自變量，以重金屬含量為因變量，分別建立線性、 二次多項式、 對數(shù)、 指數(shù)、 冪的回歸模型。通過擬合精度(R2)、 均方根誤差(RMSE)檢驗上述植被指數(shù)與重金屬回歸模型的精度，選取R2較高的公式作為各個重金屬擬合模型，結(jié)果如表6所示。

表4 逐步回歸確定的特征波段Table 4 Characteristic band determined bystepwise regression

表5 五種植被指數(shù)及其與三種重金屬含量相關(guān)系數(shù)Table 5 Vegetation indexes and its correlation coefficientswith the contents of three heavy metals

表6 三種重金屬含量的光譜模型Table 6 Spectral models of three heavy metals

2.5 重金屬空間分布的反演與驗證

以Sential-2影像為基礎(chǔ)，利用三種重金屬預(yù)測模型進(jìn)行反演，并將重金屬含量劃分為5個等級，三種金屬含量預(yù)測值的空間分布見圖4，研究發(fā)現(xiàn)：

(1)由反演結(jié)果統(tǒng)計可知，Cd含量變化范圍為0～0.85 mg·kg-1，As含量變化范圍為0～85 mg·kg-1，Pb含量變化范圍為0～98 mg·kg-1；Cd，As和Pb三種重金屬污染較為嚴(yán)重的區(qū)域面積占比分別為: 28.06%，16.82%和19.35%；表明桃林土壤受污染程度較為嚴(yán)重。

(2)圖4顯示，Cd金屬污染嚴(yán)重的區(qū)域主要分布在平谷西部劉家店鎮(zhèn)、 東北部羅營鎮(zhèn)、 東部黃松峪鄉(xiāng)及金海湖鎮(zhèn)、 東南部南獨樂河鎮(zhèn)等區(qū)域；As金屬污染嚴(yán)重的區(qū)域主要分布在西部劉家店鎮(zhèn)及峪口鎮(zhèn)、 東部金海湖鎮(zhèn)；Pb金屬污染嚴(yán)重的區(qū)域主要分布在西部劉家店鎮(zhèn)及峪口鎮(zhèn)、 東部金海湖鎮(zhèn)、 東南部南獨樂河鎮(zhèn)等區(qū)域。三種重金屬高值區(qū)均大面積的分布在西部萬莊尾礦庫及劉店尾礦庫、 東部黃松峪鄉(xiāng)及金海湖鎮(zhèn)附近，西部尾礦區(qū)比東部尾礦區(qū)重金屬污染更為嚴(yán)重。

隨機(jī)選取150組反演值和實測數(shù)據(jù)，利用均方根誤差(RMSE)和決定系數(shù)R2對三類金屬反演模型的結(jié)果進(jìn)行精度分析。RMSE越接近0，R2越大，表示模型反演效果越好，精度越高。結(jié)果如圖5所示，三種重金屬預(yù)測模型RMSE分別為0.019，0.673和5.412，R2分別為0.668，0.794和0.834。表明三種重金屬反演模型能準(zhǔn)確反映北京平谷區(qū)桃林三種重金屬空間分布。

3 結(jié) 論

以北京平谷桃林區(qū)為研究對象，通過光譜變換對桃葉片光譜數(shù)據(jù)進(jìn)行處理，選取特征波段，構(gòu)建植被指數(shù)，基于R2最大準(zhǔn)則確定了三種金屬元素的反演模型，并利用Sential-2影像對土壤重金屬含量進(jìn)行反演及空間分布成圖，得到了以下結(jié)論:

(1)研究區(qū)桃葉遭重金屬脅迫時，葉綠素大量減少，綠光波段反射率變高，其平均光譜反射率高于正常葉片，740～1 300 nm波長范圍內(nèi)差值明顯。隨重金屬污染程度加劇，紅邊位置發(fā)生了“藍(lán)移”現(xiàn)象，藍(lán)邊和黃邊位置變化不大，紅邊面積及紅邊斜率減小。

(2)將桃葉光譜曲線進(jìn)行光譜特征變換后，篩選確定780，945和1 375 nm為特征波段，構(gòu)建植被指數(shù)HMSVI。HMSVI相比于NDVI705，CRL2，SIPI，NDWI與重金屬Cd，As，Pb的相關(guān)性更高。

圖4 三種重金屬元素含量空間分布圖Fig.4 Spatial distribution of three heavy metal elements

圖5 Cd，As和Pb重金屬含量預(yù)測模型精度驗證Fig.5 Accuracy verification of heavy metal content prediction model for Cd, As and Pb

(3)以HMSVI為自變量，土壤Cd，As，Pb元素含量為因變量建立線性、 二次多項式、 對數(shù)和指數(shù)形式的預(yù)測模型，其最終預(yù)測模型分別為Y=0.44X+0.193，Y=7.436lnX+13.161，Y=-15.359X+13.583X2+23.541，R2均達(dá)到0.6以上，RMSE分別為0.22，1.394和2.403，模型穩(wěn)定、 適應(yīng)性強(qiáng)。

(4)利用實測數(shù)據(jù)驗證了反演結(jié)果的合理性, 結(jié)果表明：北京平谷區(qū)Cd，As和Pb含量變化范圍分別為0～0.85，0～85和0～98 mg·kg-1；三種重金屬高值區(qū)均大面積的分布在西部萬莊尾礦庫及劉店尾礦庫附近，東部黃松峪鄉(xiāng)及金海湖鎮(zhèn)；西部尾礦區(qū)比東部尾礦區(qū)重金屬污染更為嚴(yán)重，且在東北部羅營鎮(zhèn)存在較為嚴(yán)重的Cd污染。金屬反演結(jié)果精度結(jié)果顯示：RMSE分別為0.019，0.673和5.412，R2分別為0.668，0.794和0.834，三種模型均具有一定的預(yù)測能力。

結(jié)果表明構(gòu)建的HMSVI指數(shù)，能較好的反演出桃林區(qū)土壤重金屬Cd，As和Pb元素的分布，可用于監(jiān)測桃林土壤重金屬污染，并為植物覆蓋區(qū)間接反演土壤重金屬污染研究提供一定的科學(xué)支持。