張亞龍，黃照強(qiáng)，倪 斌，江 淼，朱富曉

（中國(guó)冶金地質(zhì)總局礦產(chǎn)資源研究院，北京 101300）

0 引言

鉛、鋅是普遍存在于土壤中的重金屬元素，由于金屬冶煉活動(dòng)會(huì)對(duì)周邊地區(qū)造成重金屬污染，對(duì)于人體健康存在重大風(fēng)險(xiǎn)。中國(guó)地質(zhì)調(diào)查局對(duì)全國(guó)土地質(zhì)量地球化學(xué)調(diào)查中發(fā)現(xiàn)，我國(guó)約230 萬(wàn)hm2的耕地屬于重金屬中-重度污染[1]。河北省地質(zhì)調(diào)查院在對(duì)河北省平原區(qū)開(kāi)展地球化學(xué)調(diào)查中發(fā)現(xiàn)，保定市東南部土壤存在重金屬污染[2]。開(kāi)展該地區(qū)鉛鋅重金屬污染情況的分析研究對(duì)于了解雄安新區(qū)土壤污染情況和土壤污染評(píng)價(jià)具有重要作用。

傳統(tǒng)的土壤重金屬研究需要進(jìn)行大量樣品的采樣和分析，需要投入巨大的人力和時(shí)間，開(kāi)展大范圍土壤重金屬研究需要的周期較長(zhǎng)[3]。近幾年，隨著高光譜遙感的快速發(fā)展，大尺度重金屬污染越來(lái)越受到廣大國(guó)內(nèi)外學(xué)者的重視，其基于高光譜遙感數(shù)據(jù)和土壤重金屬元素化驗(yàn)數(shù)據(jù)，依據(jù)鉛元素、鋅元素與光譜信息之間的響應(yīng)關(guān)系對(duì)土壤中鉛、鋅兩種重金屬元素污染進(jìn)行了大量研究。解憲麗等[4]研究土壤重金屬含量與光譜反射率之間的關(guān)系，得出Pb、Zn、Co、Ni、Fe 的最高相關(guān)系數(shù)達(dá)到高度相關(guān)（|r|＞0.80）；王璐等[5]基于土壤反射光譜利用偏最小二乘法構(gòu)建了Cd、Pb、Hg 含量預(yù)測(cè)模型；陳銀瑩等[6]通過(guò)相關(guān)分析研究土壤銅、鉻、錳、鉛、鋅的特征波段，利用元素含量、光譜反射數(shù)據(jù)變換和特征波段的關(guān)系構(gòu)建土壤重金屬高光譜反演模型；譚琨等[7]基于ASD 獲取的土壤高光譜數(shù)據(jù)，建立了As、Pb、Cr、Cu、Zn 五種重金屬元素的定量反演模型；AND 等[8]通過(guò)重金屬污染區(qū)的高光譜數(shù)據(jù)，利用多元線性回歸方法和人工神經(jīng)網(wǎng)絡(luò)方法對(duì)重金屬元素光譜數(shù)據(jù)進(jìn)行建模分析，得到的模型預(yù)測(cè)精度較高；單海斌等[9]通過(guò)六種光譜變換，研究土壤有機(jī)質(zhì)含量和土壤光譜特征的關(guān)系，建立了研究區(qū)灰漠土土壤有機(jī)質(zhì)含量高光譜反應(yīng)模型；馬磊等[10]建立了地物光譜和Landsat8 數(shù)據(jù)的偏最小二乘法土壤鉛含量反演模型，通過(guò)該模型能夠粗略預(yù)測(cè)土壤中鉛含量；阿依努爾·麥提努日等[11]對(duì)吐魯番盆地葡萄園土壤光譜進(jìn)行15種變換，并研究其與Pb 元素含量的關(guān)系，建立偏最小二乘回歸模型和地理加權(quán)重回歸模型，估算土壤中Pb 的含量；黃啟廳等[12]利用偏最小二乘法定量分析了土壤中Pb 的含量。

綜上所述，利用高光譜數(shù)據(jù)與土壤重金屬實(shí)測(cè)含量建立的反演模型可以有效預(yù)測(cè)重金屬含量。本研究選取保定市西南部三個(gè)縣城交界處土壤重金屬中Pb 元素、Zn 元素為研究對(duì)象，在室內(nèi)使用SVC（HR-1 024）地物光譜儀獲取土壤高光譜數(shù)據(jù)，與實(shí)驗(yàn)室測(cè)量的土壤Pb、Zn 含量數(shù)據(jù)相結(jié)合，分析研究區(qū)土壤的光譜特征，采用偏最小二乘法（PLS）和反向傳播神經(jīng)網(wǎng)絡(luò)法（BPNN）建立Pb 元素、Zn 元素預(yù)測(cè)估算模型，最終確認(rèn)研究區(qū)Pb、Zn 兩種金屬元素的最優(yōu)反演方法和估算模型。

1 研究區(qū)概況

研究區(qū)位于保定市清苑縣、安新縣和高陽(yáng)縣三縣交接處，主要行政區(qū)劃包含老河頭鎮(zhèn)、蘆莊鄉(xiāng)、何橋鄉(xiāng)、同口鎮(zhèn)等鄉(xiāng)鎮(zhèn)，面積約447 km2。研究區(qū)內(nèi)有色金屬、冶煉、紡織業(yè)、紙業(yè)、蓄電池等產(chǎn)業(yè)發(fā)達(dá)。該地區(qū)在歷史上存在大量的多金屬冶煉廠，承擔(dān)著鋼鐵冶煉功能，工業(yè)的粗放式發(fā)展和廢棄冶煉廠尾礦渣的堆放，隨著降雨、大風(fēng)等原因造成部分耕地土壤環(huán)境鉛鋅重金屬元素超標(biāo)[13-16]。

2 研究方法

2.1 樣品采集與化驗(yàn)分析

2.1.1 土壤樣品采集

采樣時(shí)間為2020 年9 月—11 月，樣品采集深度為0～20 cm 的耕作層土壤，共采集438 件，采樣點(diǎn)位分布如圖1 所示。樣品采集采用S 形布設(shè)，由4 處子采樣點(diǎn)樣品組合成一個(gè)混合樣。采樣時(shí)使用手持GPS 對(duì)采樣點(diǎn)經(jīng)緯度進(jìn)行采集，并且將采樣點(diǎn)詳細(xì)信息錄入土壤信息采集系統(tǒng)。

圖1 研究區(qū)及采樣點(diǎn)位置圖Fig.1 Location map of study area and sampling points

2.1.2 樣品化驗(yàn)

本次采集的土壤樣品由中國(guó)冶金地質(zhì)總局物勘院測(cè)試中心進(jìn)行重金屬元素的化驗(yàn)，樣品經(jīng)風(fēng)干（確保無(wú)污染）、過(guò)篩（孔徑0.074 mm 標(biāo)準(zhǔn)篩）、送樣（采用四分法混合均勻），每個(gè)采樣點(diǎn)樣品分成兩份，一份用于化學(xué)分析，采用X 射線熒光光譜儀，利用波長(zhǎng)色散型X 射線熒光光譜方法進(jìn)行檢測(cè)；另一份用于室內(nèi)光譜測(cè)量。

2.2 樣品光譜測(cè)定及處理

樣品光譜測(cè)量?jī)x器為SVC HR-1024 型地物光譜儀，光譜波長(zhǎng)范圍350～2 500 nm。樣品光譜測(cè)定在暗室中進(jìn)行，盛放于黑色器皿中，將表面刮平后，進(jìn)行光譜測(cè)定。以接觸式高密度反射探頭測(cè)定，該種光譜測(cè)定受外界雜散光干擾較少，獲得的土壤光譜數(shù)據(jù)更加精確。測(cè)量中始終確保土樣和探頭之間完全接觸，在測(cè)定下一個(gè)樣本時(shí)，將接觸探頭上的殘留土樣清除干凈，盡量避免各土樣間相互混合、相互干擾。每個(gè)樣品進(jìn)行5 次光譜測(cè)量，取平均值，以減少誤差，提高精度，每測(cè)定5 個(gè)樣本，需要做一次白板校準(zhǔn)，確保每個(gè)土樣同樣的測(cè)定標(biāo)準(zhǔn)。利用ENVI 軟件剔除異常曲線后，取光譜反射率平均值作為該樣品原始實(shí)測(cè)光譜數(shù)據(jù)。

對(duì)土壤樣品的原始實(shí)測(cè)光譜進(jìn)行平滑處理，再對(duì)其分別進(jìn)行一階導(dǎo)數(shù)（FD）、二階導(dǎo)數(shù)（SD）、連續(xù)統(tǒng)去除（CR）、多元線性回歸（MSC）、歸一化（Normalize）和標(biāo)準(zhǔn)正態(tài)變換（SNV）六種光譜變換處理，以減少背景噪聲。

2.3 單因子污染指數(shù)法

單因子污染指數(shù)法[17]是計(jì)算土壤中污染物指標(biāo)i的單項(xiàng)污染指數(shù)Pi，計(jì)算見(jiàn)式（1）。

式中：Ci為土壤中污染物指標(biāo)i的實(shí)測(cè)質(zhì)量分?jǐn)?shù)，mg/kg；Si為土壤中污染物指標(biāo)i在《土壤環(huán)境質(zhì)量農(nóng)用地土壤污染風(fēng)險(xiǎn)管控標(biāo)準(zhǔn)（試行）》（GB 15618—2018）中給出的二級(jí)標(biāo)準(zhǔn)值，mg/kg。Pi≤1 時(shí)土壤重金屬含量在土壤背景值范圍內(nèi)，土壤未受到人類活動(dòng)污染；Pi＞1 時(shí)土壤受到人為污染，指數(shù)值越大土壤重金屬累積污染程度越大，土壤環(huán)境地球化學(xué)等級(jí)劃分見(jiàn)表1[17-18]。

表1 土壤環(huán)境地球化學(xué)等級(jí)劃分Table 1 Geochemistry classification of soil environment

2.4 反演模型建立與分析

2.4.1 相關(guān)性分析

皮爾遜（Pearson）相關(guān)系數(shù)是評(píng)價(jià)相關(guān)性的常用指標(biāo)之一，表示兩變量間的協(xié)方差與標(biāo)準(zhǔn)差的商，取值范圍為[-1,1]。rXY越大，變量X與變量Y相關(guān)性越高，1 表示正相關(guān)，-1 表示負(fù)相關(guān)，0 表示不相關(guān)。

2.4.2 偏最小二乘法

偏最小二乘法[19-20]（Partial Least Squares，PLS）是一種多自變量對(duì)多因變量的線性回歸建模方法，可同時(shí)實(shí)現(xiàn)提取變量特征、分析變量間相關(guān)性和回歸建模，適合多元數(shù)據(jù)的統(tǒng)計(jì)分析。本文中PLS 模型運(yùn)算在The UnscramblerX10.4 軟件中完成。

2.4.3 反向傳播神經(jīng)網(wǎng)絡(luò)法

反向傳播神經(jīng)網(wǎng)絡(luò)[21]（Back Propagation Neural Network，BPNN）是一種有效的多層神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)方法，通過(guò)反向傳播不斷調(diào)整網(wǎng)絡(luò)的權(quán)值和閾值，使網(wǎng)絡(luò)的最終輸出盡可能接近預(yù)期輸出，從而達(dá)到訓(xùn)練目的。當(dāng)網(wǎng)絡(luò)處于學(xué)習(xí)過(guò)程時(shí)，信號(hào)從輸入層傳輸?shù)捷敵鰧印Ｈ绻敵鼋Y(jié)果不符合目標(biāo)，則將梯度反饋到網(wǎng)絡(luò)中，以調(diào)整每個(gè)神經(jīng)元的權(quán)重和偏差，并最小化預(yù)測(cè)數(shù)據(jù)和實(shí)際數(shù)據(jù)之間的誤差。本文中BPNN 算法是基于軟件Matlab 2020 版本編寫(xiě)和實(shí)現(xiàn)。

2.5 模型驗(yàn)證

模型的穩(wěn)定性用決定系數(shù)R2的大小來(lái)檢驗(yàn)，模型的精度用均方根誤差（Root Mean Square Error，RMSE）來(lái)檢驗(yàn)，R2越大、RMSE值越小，模型性能越好。

3 結(jié)果與討論

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

本次對(duì)438 件土壤樣品進(jìn)行統(tǒng)計(jì)分析，樣品中Pb 含量、Zn 含量統(tǒng)計(jì)分析見(jiàn)表2。研究區(qū)Pb、Zn 兩種元素含量平均值均高于冀中平原背景值和中國(guó)土壤背景值，分別為3.6 倍、2.3 倍和3.1 倍、2.2 倍，兩個(gè)元素在研究區(qū)存在一定的富集，Pb 元素的富集程度較高于Zn 元素。與土壤污染風(fēng)險(xiǎn)篩選值相比，Pb元素、Zn 元素平均值均低于篩選值，但最大值遠(yuǎn)高于篩選值，分別為4.3 倍和3.5 倍，兩種金屬在研究區(qū)存在局部超標(biāo)的情況。從變異系數(shù)來(lái)看，兩種金屬的變異系數(shù)在0.9 以上，Pb 元素變異程度稍高于Zn元素。

表2 重金屬Pb 元素、Zn 元素實(shí)測(cè)含量統(tǒng)計(jì)Table 2 Statistics of measured content of heavy metal Pb and Zn elements

根據(jù)土壤重金屬污染的單因子指數(shù)法測(cè)算，樣品中Pb 元素有413 個(gè)屬于清潔，占比94.29%；18 個(gè)屬于輕微污染，占比4.11%；2 個(gè)屬于輕度污染，占比0.46%；5 個(gè)屬于中度污染，占比1.14%。樣品中Zn 元素有398 個(gè)屬于清潔，占比90.87%；25 個(gè)屬于輕微污染，占比5.71%；10 個(gè)屬于輕度污染，占比2.28%；5 個(gè)屬于中度污染，占比1.14%（表1 和表3）?？傮w來(lái)說(shuō)，研究區(qū)土壤中重金屬Pb 元素、Zn 元素存在一定污染。

表3 重金屬Pb 元素、Zn 元素污染指數(shù)Pi 統(tǒng)計(jì)分析Table 3 Statistical analysis of pollution index Pi for heavy metal Pb and Zn elements單位：個(gè)

3.2 空間分布特征分析

為反映研究區(qū)土壤重金屬Pb 元素、Zn 元素空間分布特征，利用ArcFIS 軟件地統(tǒng)計(jì)學(xué)模塊對(duì)Pb 元素、Zn 元素污染指數(shù)Pi進(jìn)行反距離權(quán)重插值（IDW），得到Pb 元素、Zn 元素的污染指數(shù)分布圖，并按照土壤環(huán)境地球化學(xué)等級(jí)劃分標(biāo)準(zhǔn)，以不同顏色對(duì)應(yīng)不同土壤環(huán)境加以界定劃分，結(jié)果如圖2 所示。

圖2 重金屬元素污染指數(shù)分布圖Fig.2 Pollution index distribution of heavy metal elements

由圖2 可知，研究區(qū)土壤中Pb、Zn 污染主要集中在安新縣和清苑縣交界處，其中，Pb 污染區(qū)分布于工作區(qū)中部，主要在安新縣老河頭鎮(zhèn)和清苑縣望亭鎮(zhèn)、何橋鄉(xiāng)周邊；Zn 污染區(qū)主要分布于工作區(qū)西南部，主要在清苑縣何橋鄉(xiāng)周邊。

3.3 Pb、Zn 含量與光譜變換后反射率的相關(guān)性分析

將原始光譜進(jìn)行FD、SD、CR、MSC、Normalize和SNV 六種光譜變換，再分別對(duì)Pb 元素、Zn 元素含量和六種光譜變換進(jìn)行皮爾遜（Pearson）相關(guān)性分析，比較不同光譜變換形式對(duì)相關(guān)程度的影響，得出兩種元素與六種光譜變換后最大相關(guān)系數(shù)及對(duì)應(yīng)的特征波段，見(jiàn)表4。

表4 Pb、Zn 含量與光譜反射率間最大相關(guān)系數(shù)Table 4 Maximum correlation coefficient between Pb,Zn contents and spectral reflectance

由表4 可知，Pb 元素與六種光譜變換的相關(guān)系數(shù)中，Normalize 變換相關(guān)系數(shù)絕對(duì)值不足0.400，其余相關(guān)系數(shù)絕對(duì)值均高于0.400，其中，MSC 變換在波段553.2 nm 處相關(guān)性最高，為-0.451。Zn 元素與六種光譜變換的相關(guān)系數(shù)中，Normalize 變換相關(guān)系數(shù)絕對(duì)值不足0.400，其余相關(guān)系數(shù)絕對(duì)值均高于0.450，其中，SD 變換在波段746.5 nm 處相關(guān)性最高，為-0.564。

3.4 PLS 分析

將438 個(gè)樣本按實(shí)測(cè)重金屬元素含量從高到低排序，每3 個(gè)樣本抽取1 個(gè)驗(yàn)證樣本，最終選取292個(gè)樣本作為建模樣本，剩余146 個(gè)樣本用作模型的驗(yàn)證樣本。選取六種光譜變換后的反射率作為模型的自變量，土壤重金屬含量實(shí)測(cè)值為因變量，建立PLS模型。

3.4.1 Pb 元素PLS 分析

基于六種光譜變換，通過(guò)PLS 建模和驗(yàn)證結(jié)果發(fā)現(xiàn)，以SD 為自變量的PLS 模型的建模精度R2最高，為0.441，RMSE值為60.848；以FD 變換為自變量的PLS 模型的驗(yàn)證精度R2最 高，為0.529，RMSE值為48.276（表5）。以驗(yàn)證樣本Pb 含量的實(shí)測(cè)值作為橫坐標(biāo)，PLS 模型中Pb 含量預(yù)測(cè)值作為縱坐標(biāo)，得出實(shí)測(cè)值與預(yù)測(cè)值的散點(diǎn)圖（圖3）。對(duì)比六種光譜變換的模型精度，綜合建模精度、驗(yàn)證精度和均方根誤差，其中，基于FD 和SD 的PLS 模型預(yù)測(cè)效果比較相近，效果較好。

表5 Pb 元素PLS 模型結(jié)果Table 5 Results of PLS model for Pb element

圖3 基于PLS 建立的Pb 元素模型精度評(píng)價(jià)Fig.3 Accuracy evaluation of Pb element model based on PLS

3.4.2 Zn 元素PLS 分析

基于六種光譜變換通過(guò)PLS 建模和驗(yàn)證結(jié)果發(fā)現(xiàn)，以SD 為自變量的PLS 模型建模集精度和驗(yàn)證集精度最高，分別為0.632 和0.445，RMSE值分別為94.111 和109.289（表6）。以驗(yàn)證樣本Zn 含量的實(shí)測(cè)值作為橫坐標(biāo)，PLS 模型中Zn 含量預(yù)測(cè)值作為縱坐標(biāo)，得出實(shí)測(cè)值與預(yù)測(cè)值的散點(diǎn)圖（圖4）。對(duì)比六種光譜變換模型精度，綜合建模精度、驗(yàn)證精度和均方根誤差，其中，基于SD 的PLS 散點(diǎn)圖基本分布在1∶1 線附近，模型預(yù)測(cè)效果相對(duì)比較理想。

表6 Zn 元 素PLS 模型結(jié)果Table 6 Results of PLS model for Zn element

圖4 基于PLS 建立的Zn 元素模型精度評(píng)價(jià)Fig.4 Accuracy evaluation of Zn element model based on PLS

3.5 BPNN 分析

3.5.1 Pb 元素BPNN 分析

基于六種光譜變換，通過(guò)BPNN 建模和驗(yàn)證發(fā)現(xiàn)，其中，以SNV 變換為自變量的BPNN 模型建模精度R2較高，為0.949，RMSE值為18.531；以FD 變換為自變量的BPNN 模型驗(yàn)證精度R2較高，為0.991，RMSE值為8.461（表7）。以驗(yàn)證樣本Pb 含量的實(shí)測(cè)值作為橫坐標(biāo)，BPNN 模型中Pb 含量預(yù)測(cè)值作為縱坐標(biāo)，得出實(shí)測(cè)值與預(yù)測(cè)值的散點(diǎn)圖（圖5）。對(duì)比六種光譜變換模型精度，綜合建模精度、驗(yàn)證精度和均方根誤差，基于SNV 變換的BPNN 散點(diǎn)圖都基本分布在1∶1 線附近，模型預(yù)測(cè)效果最理想。

表7 Pb 元素BPNN 模型結(jié)果Table 7 Results of BPNN model for Pb element

圖5 基于BPNN 建立Pb 元素模型精度評(píng)價(jià)Fig.5 Accuracy evaluation of Pb element model based on BPNN

3.5.2 Zn 元素BPNN 分析

基于六種光譜變換通過(guò)BPNN 建模和驗(yàn)證發(fā)現(xiàn)，其中，以CR 變換為自變量的BPNN 模型建模精度和驗(yàn)證精度較高，建模精度R2為0.874，RMSE值為53.183，驗(yàn)證精度R2為0.957，RMSE值為35.721（表8）。以驗(yàn)證樣本Zn 含量的實(shí)測(cè)值作為橫坐標(biāo)，BPNN 模型中Zn 含量預(yù)測(cè)值作為縱坐標(biāo)，得出實(shí)測(cè)值與預(yù)測(cè)值的散點(diǎn)圖（圖6）。對(duì)比六種光譜變換模型精度，綜合建模精度、驗(yàn)證精度和均方根誤差，基于CR 變換的BPNN 散點(diǎn)圖都基本分布在1∶1 線附近，模型預(yù)測(cè)效果最理想。

表8 Zn 元素BPNN 模型結(jié)果Table 8 Results of BPNN model for Zn element

圖6 基于BPNN 建立Zn 元素模型精度評(píng)價(jià)Fig.6 Accuracy evaluation of Zn element model based on BPNN

4 結(jié)論

1）研究區(qū)土壤中Pb 元素、Zn 元素的平均含量均遠(yuǎn)高于冀中平原土壤和中國(guó)土壤的背景值，Pb 元素屬于強(qiáng)變異，Zn 元素屬于中等變異，均存在局部重金屬超標(biāo)的情況，Pb、Zn 兩種重金屬元素已在研究區(qū)內(nèi)造成了一定的重金屬污染。

2）通過(guò)對(duì)研究區(qū)土壤Pb 元素、Zn 元素含量統(tǒng)計(jì)分析發(fā)現(xiàn)，土壤中Pb、Zn 污染占比為分別為5.71%和9.13%，其中大部分為輕微污染和輕度污染，均不存在重度污染。

3）通過(guò)PLS 模型和BPNN 模型對(duì)研究樣品進(jìn)行建模分析，在PLS 方法中，Pb 元素FD 模型和SD 模型精度較好，Zn 元素SD 模型精度較好；在BPNN 方法中，Pb 元素SNV 模型精度較好，Zn 元素CR 模型精度較好。

4）通過(guò)兩種建模方法對(duì)于Pb、Zn 兩種元素的建模精度對(duì)比，BPNN 的模型精度要高于PLS。其中，Pb 元素應(yīng)用BPNN 基于SNV 建立的估算模型中建模精度R2為0.949，驗(yàn)證精度R2為0.980，為最佳估算模型；Zn 元素應(yīng)用BPNN 基于CR 建立的估算模型中建模精度R2為0.874，驗(yàn)證精度R2為0.957，為最佳估算模型。