吳 兵，楊可明，高 偉，李艷茹，韓倩倩，張建紅

中國礦業(yè)大學(xué)(北京)地球科學(xué)與測繪工程學(xué)院，北京 100083

引 言

工業(yè)化和城市化發(fā)展以及一些自然或人為干擾，使人類賴以生存的土壤中重金屬濃度逐漸增大，當(dāng)土壤中重金屬含量超過了當(dāng)?shù)丨h(huán)境背景值時，其累積效應(yīng)會影響土地承載力，最終導(dǎo)致土壤重金屬污染[1-2]，引發(fā)生態(tài)環(huán)境惡化和糧食安全危機。重金屬污染物帶有毒性并且不可降解，一旦通過生態(tài)系統(tǒng)或食物鏈進入人體，就會對人類身體健康構(gòu)成嚴重威脅[3-4]。傳統(tǒng)的重金屬污染檢測方法通常是對野外采集的樣品開展實驗室化學(xué)分析，該方法測量精度高、準確性強，但由于檢測環(huán)節(jié)多、耗時長、成本高，很難快速獲取大區(qū)域的污染物含量及分布信息[5]?？焖倥袆e重金屬污染的元素種類以及污染程度，為污染區(qū)的針對性治理計劃及時提供依據(jù)，成為當(dāng)前亟待解決的關(guān)鍵步驟[6-7]。

由于高光譜遙感具有光譜信息量豐富、技術(shù)成本低、理化反演能力強，分析操作實時便捷、可非接觸式無損檢測等優(yōu)點[8]，從而如何將光譜分析技術(shù)應(yīng)用于重金屬污染檢測已成為了現(xiàn)今遙感領(lǐng)域研究的一個熱點，并取得了一定進展。Hou等利用可見-近紅外光譜數(shù)據(jù)分析山東省鄒城市某煤礦周邊的土壤重金屬濃度，發(fā)現(xiàn)利用光譜數(shù)據(jù)建立的偏最小二乘回歸模型預(yù)測土壤中Cu和Pb重金屬濃度精度較高[9]。也有一些學(xué)者研究表明，土壤中的重金屬能對植物的生理結(jié)構(gòu)特征產(chǎn)生影響，當(dāng)土壤中的重金屬含量過高時，會導(dǎo)致植物葉片變黃、分枝、落葉甚至死亡[10]，從而改變植物的光譜特性。在農(nóng)作物重金屬污染方面，楊可明等利用玉米葉片光譜的分數(shù)階微分(fractional order differential，F(xiàn)OD)信息構(gòu)建玉米葉片光譜紅邊位置的銅鉛敏感指數(shù)集群，實現(xiàn)銅鉛脅迫識別[11]；黃鐘霆等以湖南省某典型關(guān)停錳礦區(qū)為研究對象，采集礦區(qū)周邊的農(nóng)作物及其對應(yīng)的土壤樣品，探討土壤及對應(yīng)農(nóng)作物間重金屬遷移規(guī)律，并進行生態(tài)風(fēng)險評價[12]。上述研究在重金屬元素含量預(yù)測以及識別方面均得到了較好的結(jié)果，但光譜處理繁瑣，算法較為復(fù)雜，計算量大，不易復(fù)現(xiàn)。

以不同銅鉛離子濃度梯度脅迫生長下的玉米葉片光譜為研究對象，使用ASD FieldSpec 4便攜式地物光譜儀獲取不同脅迫梯度下穗期玉米葉片光譜數(shù)據(jù)，基于常規(guī)的包絡(luò)線去除(continuum removed，CR)，光譜比值(spectral ratio，SR)以及分數(shù)階微分(FOD)光譜預(yù)處理方法構(gòu)建銅鉛元素識別指數(shù)(copper and lead identification index，CLI)，聯(lián)合歐式聚類(euclidean cluster，EC)與垂直平分線(perpendicular bisector，PB)法建立二維(2D)與三維(3D)銅、鉛污染元素判別規(guī)則，實現(xiàn)依據(jù)被污染玉米葉片光譜的銅鉛元素準確識別，為重金屬元素快速識別，污染區(qū)域治理提供新思路。

1 實驗部分

1.1 試驗設(shè)計

玉米是我國重要的農(nóng)作物之一，也是我國種植面積最大的農(nóng)作物。為了探究重金屬對農(nóng)作物的影響，試驗以盆栽種植的玉米植株為研究對象，設(shè)置CuSO4·5H2O和Pb(NO3)2含量分別為0，50，100，150，200，300，400，600，800和1 000 μg·g-1的污染脅迫梯度土壤用于培育植株生長，培育場所室溫可控，水氣環(huán)境適宜，營養(yǎng)元素充足，不同重金屬脅迫梯度設(shè)置3個平行試驗樣本組。光譜數(shù)據(jù)采集在暗室進行，使用波段范圍為350～2 500 nm的ASD FieldSpec 4便攜式地物光譜儀獲取不同脅迫梯度下穗期玉米葉片光譜數(shù)據(jù)，每個樣本測量3次光譜后剔除異常信息取平均光譜，使用電感耦合等離子發(fā)射光譜儀測定玉米葉片Cu2+和Pb2+含量(如表1所示)并用于研究結(jié)果驗證。

表1 不同濃度梯度下玉米葉片中Cu2+，Pb2+含量

本試驗共設(shè)置9個不同銅鉛濃度梯度脅迫下生長的玉米植株，每個濃度梯度設(shè)置3個平行組，共54組樣本數(shù)據(jù)，由Python中的sklearn庫隨機劃分樣本數(shù)據(jù)的70%(38組)為訓(xùn)練集數(shù)據(jù)，30%(16組)為驗證集數(shù)據(jù)。

1.2 光譜變換

(1)包絡(luò)線去除(CR)。CR也稱連續(xù)統(tǒng)去除，是一種突出光譜曲線的吸收和反射特征，并將反射率歸一化的光譜預(yù)處理方法。根據(jù)Clark提出的外殼系數(shù)法[13]，通過計算光譜曲線上的最大極大值點作為包絡(luò)線的一個端點，計算該點與波長增長方向各個極大值點連線的斜率，以斜率最大點作為包絡(luò)線的下一個端點，以此循環(huán)，將所有端點連接，形成光譜曲線的包絡(luò)線，用實際光譜曲線的反射率值除以包絡(luò)線上相應(yīng)波段的反射率值即可得到原始光譜曲線去除包絡(luò)線后的結(jié)果，如式(1)和式(2)所示。

(1)

(2)

式中，Si為波段i的包絡(luò)線去除值，K為極大值始點和端點之間的斜率，ρi為波段i的原始光譜反射率，ρs和ρe為極大值始點和端點的原始光譜反射率，λs和λe為ρs和ρe對應(yīng)的波段值。本研究中所有玉米葉片光譜數(shù)據(jù)均進行了CR處理。

(2)光譜比值(SR)。SR最早應(yīng)用于消除遙感圖像中的陰影。隨著高光譜技術(shù)的成熟，比值技術(shù)在混合高光譜信號解混中表現(xiàn)出巨大潛力，其核心思想是通過相同波段不同光譜數(shù)據(jù)反射率的比值進行SR變換，SR技術(shù)可以抑制作為分母的光譜特征，而突出分子光譜的影響。本研究中為了突出重金屬元素銅鉛的光譜特征，以相同環(huán)境下受不同程度銅鉛污染的玉米葉片光譜作為分子，無污染的正常的玉米葉片光譜作為分母進行SR處理。

(3)分數(shù)階微分(FOD)。光譜微分可以消除光譜數(shù)據(jù)之間的系統(tǒng)誤差，減弱大氣輻射、散射和大氣吸收對待測物體光譜的影響，突出待測光譜曲線特征吸收峰和斜率的微小變化，光譜的FOD是整數(shù)階微分的推廣，根據(jù)Grünwald-Letnikov定義形式[14]，F(xiàn)OD計算方式如式(3)所示。

(3)

式(3)中，λ表示光譜曲線某一波段值，[m,n]為波段區(qū)間，λ∈[m,n]，h為光譜采樣間隔，a為常數(shù)參數(shù)，q為任意階數(shù)，Г為Gamma函數(shù)，其積分定義如式(4)所示。

(4)

本研究中對經(jīng)過CR、SR處理后的玉米葉片光譜按照0.1的間隔進行0～2的各階次FOD處理。

1.3 銅鉛污染元素判別特征點(CLDFP)

作物中所含的重金屬元素屬于微量元素，在光譜曲線中即使經(jīng)過了增強處理也很難僅通過光譜曲線進行重金屬元素識別，因此需要建立一種判別規(guī)則對作物樣本中的重金屬元素進行準確的識別。本研究基于改進紅邊比值指數(shù)(modified red edge simple ratio index，MSR)[15]建立一種銅鉛元素識別指數(shù)(CLI)如式(5)所示。

(5)

式(5)中，ρ(λa)和ρ(λb)分別表示波長為λa和λb處的光譜值。本研究中將經(jīng)過CR-SR-FOD處理后的各組玉米葉片光譜代入CLI中，挑選出與銅鉛元素種類相關(guān)性最強的三個階數(shù)的CLI，分別為CLI1，CLI2和CLI3。以CLI1和CLI2為坐標分量構(gòu)建2D坐標系，(CLI1為x軸，CLI2為y軸)；以CLI1，CLI2和CLI3為坐標分量構(gòu)建3D坐標系，(CLI1為x軸，CLI2為y軸，CLI3為z軸)構(gòu)建銅鉛元素判別特征點(copper and lead discriminantfeature points，CLDFP)，如式(6)和式(7)所示。

CLDFP2D=(CLI1，CLI2)

(6)

CLDFP3D=(CLI1，CLI2，CLI3)

(7)

1.4 銅鉛污染元素判別規(guī)則

歐式聚類(EC)[16]即基于歐幾里得距離的一種聚類方法，針對CLDFP需要確定不同元素各樣本點內(nèi)最大的兩個歐式距離d1和d2，2D與3D坐標系下的歐式距離公式如式(8)和式(9)所示。

(8)

式(8)中，x1，y1，x2和y2分別為2D坐標系下不同污染元素各樣本點內(nèi)最大歐式距離的CLI值。

(9)

式(9)中，x1，y1，z1，x2，y2和z2分別為3D坐標系下不同污染元素各樣本點內(nèi)最大歐式距離的CLI值。

在2D坐標系中，以最大歐式距離d2D1和d2D2為直徑作圓(圓心為兩點的中點)，將訓(xùn)練集樣本分為銅污染和鉛污染兩類，連接兩圓圓心，以圓心連線的垂直平分線(PB)作為銅鉛元素的判別規(guī)則線(copper and lead discriminant rule lines，CLDRL2D)；同理，在3D坐標系中，以最大歐式距離d3D1和d3D2為直徑作球(球心為兩點的中點)，將訓(xùn)練集樣本分為銅污染和鉛污染兩類，連接兩球球心，在3D坐標系中將垂直平分線進一步推廣為垂直平分面，以球心連線的垂直平分面作為銅鉛元素的判別規(guī)則面(copper and lead discriminant rule planes，CLDRP3D)實現(xiàn)玉米葉片光譜重金屬銅鉛元素種類的準確識別。

2 結(jié)果與討論

2.1 光譜變換處理與相關(guān)性分析

不同脅迫梯度下的玉米葉片原始光譜數(shù)據(jù)如圖1所示。通過圖1可以看出不同脅迫梯度的玉米葉片原始光譜曲線均表現(xiàn)出典型植被光譜特征，不同Cu2+和Pb2+濃度脅迫下僅在反射率中表現(xiàn)出部分差異，規(guī)律性不強，很難通過原始數(shù)據(jù)實現(xiàn)銅鉛元素類別的區(qū)分，因此需要進行不同的光譜變換處理增強光譜特征以及不同污染元素光譜曲線之間的差異性。

圖1 不同脅迫梯度下的玉米葉片原始光譜數(shù)據(jù)

以土壤中CuSO4·5H2O含量為50 μg·g-1下培育的玉米植株葉片光譜數(shù)據(jù)為例進行不同光譜變換處理，結(jié)果如圖2所示。通過圖2可以看出經(jīng)過包絡(luò)線去除(CR)處理后的光譜曲線峰谷特征表現(xiàn)的更加明顯；光譜比值(SR)處理后的光譜曲線已無植被的光譜特征，進一步突出了銅元素的光譜特征；1.2階分數(shù)階微分(FOD)處理后的光譜曲線消除了光譜數(shù)據(jù)中的無用信息，使得光譜數(shù)據(jù)中的有用信息能夠更有效的被利用。

圖2 不同光譜變換處理的光譜

相關(guān)性可以反映數(shù)據(jù)間的關(guān)聯(lián)程度。通過計算不同光譜變換處理后各波段光譜與葉片中重金屬銅鉛元素種類之間的相關(guān)系數(shù)，反映光譜信息與銅鉛元素種類的關(guān)聯(lián)性。分析得出：不作任何變形處理的原始光譜數(shù)據(jù)與葉片中銅鉛元素種類之間相關(guān)系數(shù)絕對值的最大值為0.21，平均值為0.09；而CR處理后的光譜信息與葉片中銅鉛元素種類之間相關(guān)系數(shù)絕對值的最大值為0.51，平均值為0.17；SR是在CR基礎(chǔ)上的變換，其相關(guān)性與CR相同；0.1～2.0各階次FOD處理后的光譜信息與葉片中銅鉛元素種類之間的相關(guān)系數(shù)如表2所示。

表2 各階次FOD與重金屬銅鉛元素種類相關(guān)系數(shù)

通過觀察不同光譜變換處理后各波段光譜數(shù)據(jù)與葉片中銅鉛元素種類之間的相關(guān)系數(shù)，可以發(fā)現(xiàn)經(jīng)過光譜變換處理后相關(guān)系數(shù)最大值、平均值與原始數(shù)據(jù)相比均有較為明顯的增加，說明光譜變換處理增加了光譜數(shù)據(jù)與銅鉛元素種類之間的關(guān)聯(lián)程度；不同階次FOD處理后相關(guān)系數(shù)各不相同，為了更好地區(qū)分銅鉛元素種類，需挑選出相關(guān)性較高的階次用于后續(xù)研究。

2.2 建立與挑選CLI

訓(xùn)練集樣本數(shù)據(jù)經(jīng)過CR、SR以及各階次FOD處理后代入式(5)計算銅鉛元素識別指數(shù)CLI并計算0.1～2.0各階次FOD對應(yīng)的CLI與重金屬銅鉛元素種類之間的相關(guān)性，相關(guān)系數(shù)隨波長變化情況如圖3所示。

圖3 各階次CLI與銅鉛元素種類相關(guān)性

其中，0.1—0.4階FOD對應(yīng)的CLI與銅鉛元素種類相關(guān)系數(shù)值在-0.60～0.63之間；0.5—0.7階FOD對應(yīng)的CLI與銅鉛元素種類相關(guān)系數(shù)值在-0.70～0.69之間；0.8—0.9階FOD對應(yīng)的CLI與銅鉛元素種類相關(guān)系數(shù)值在-0.67～0.66之間；1.0—1.2階FOD對應(yīng)的CLI與銅鉛元素種類相關(guān)系數(shù)值在-0.72～0.70之間；1.3—1.5階FOD對應(yīng)的CLI與銅鉛元素種類相關(guān)系數(shù)值在-0.67～0.69之間；1.6—2.0階FOD對應(yīng)的CLI與銅鉛元素種類相關(guān)系數(shù)值在-0.66～0.67之間。通過各階次FOD對應(yīng)的CLI與銅鉛元素種類相關(guān)系數(shù)可以發(fā)現(xiàn)，隨著階次的增加，相關(guān)性呈現(xiàn)先遞增，后遞減的趨勢，其中相關(guān)系數(shù)最高的三個階次分別為1.2階、0.7階、1.0階。

以1.2階CLI作為CLI1，其構(gòu)成波段為639 nm(λa)和1 150 nm(λb)，與銅鉛元素種類相關(guān)系數(shù)為-0.719；以0.7階CLI作為CLI2，其構(gòu)成波段為567 nm(λa)和393 nm(λb)，與銅鉛元素種類相關(guān)系數(shù)為-0.701；以1.0階CLI作為CLI3，其構(gòu)成波段為2 009 nm(λa)和636 nm(λb)，與銅鉛元素種類相關(guān)系數(shù)為0.700。

2.3 識別效果與驗證

(1)銅鉛污染元素種類的判別規(guī)則線(CLDRL)

分別將CLI1和CLI2對應(yīng)的訓(xùn)練集樣本值帶入2D坐標系，形成二維銅鉛元素判別特征點CLDFP2D，逐點尋找不同元素各樣本點內(nèi)最大兩點的歐式距離d2D1和d2D2進行銅污染與鉛污染的分類，其中銅污染判別特征點間最大距離點為(1.267，0.701)與(-0.527，0.329)，兩點間歐式距離d2D1為1.833，即以圓心為(0.370，0.515)(最大距離兩點中點)，半徑為0.916(最大距離的二分之一)作銅污染判別特征類；鉛污染判別特征點間最大距離點為(-0.424，1.665)與(2.970，1.457)，兩點間歐式距離d2D2為3.400，即以圓心為(1.273，1.561)(最大距離兩點中點)，半徑為1.700(最大距離的二分之一)作鉛污染判別特征類，連接兩圓圓心，圓心連線的垂直平分線即為CLDRL2D，其表達式如式(10)所示，2D坐標系下銅鉛判別特征點與判別規(guī)則線分布情況如圖4所示。

CLDRL2D：y=-0.863 0x+1.746 6

(10)

通過圖4可以明顯看出CLDFP2D根據(jù)兩個判別特征類分成了不同的區(qū)域，大部分CLDFP2D可以由CLDRL2D正確區(qū)分，總訓(xùn)練集樣本為38個，判別正確樣本為30個，判別正確率達到了78.95%。

圖4 訓(xùn)練集中2D銅鉛判別特征點與判別規(guī)則線分布

將CLI1和CLI2對應(yīng)的驗證集樣本值帶入2D坐標系，以相同的CLDRL2D進行判別，結(jié)果如圖5所示，總驗證集樣本為16個，判別正確樣本為12個，判別正確率為75.0%。

圖5 驗證集中2D銅鉛判別特征點與判別規(guī)則線分布

(2)銅鉛污染元素種類的判別規(guī)則面(CLDRP)

為了增加樣本的直觀性，對2D判別規(guī)則進一步推廣，將CLI1，CLI2和CLI3對應(yīng)的訓(xùn)練集樣本值帶入3D坐標系，構(gòu)建三維銅鉛元素判別特征點CLDFP3D，逐點尋找不同元素各樣本點內(nèi)最大兩點的歐式距離d3D1和d3D2進行銅污染與鉛污染的分類，其中銅污染判別特征點間最大距離點為(0.169，1.689，-0.831)與(0.447，0.517，2.266)，兩點間歐式距離d3D1為3.323，即以球心為(0.308，1.103，0.717)(最大距離兩點中點)，半徑為1.662(最大距離的二分之一)作銅污染判別特征類；鉛污染判別特征點間最大距離點為(0.659，1.578，3.290)與(2.114，0.670，-0.904)，兩點間歐式距離d3D2為4.531，即以球心為(1.387，1.124，1.193)(最大距離兩點中點)，半徑為2.266(最大距離的二分之一)作鉛污染判別特征類，連接兩球球心，球心連線的垂直平分面即為CLDRP3D，其表達式如式(11)所示，3D坐標系下銅鉛判別特征點與判別規(guī)則面分布情況如圖6所示。

圖6 訓(xùn)練集中3D銅鉛判別特征點與判別規(guī)則面分布

CLDRP3D：z=-2.268 0x-0.043 0y+2.924 6

(11)

通過圖6可以發(fā)現(xiàn)CLDFP3D同樣根據(jù)兩個判別特征類形成了不同的區(qū)域，判別效果更為直觀，由于新增加了一個維度，誤差也在相應(yīng)累加，總訓(xùn)練集樣本為38個，判別正確樣本為29個，判別正確率為76.32%。

將CLI1，CLI2和CLI3對應(yīng)的驗證集樣本值帶入3D坐標系，以相同的CLDRP3D進行判別，結(jié)果如圖7所示，總驗證集樣本為16個，判別正確樣本為12個，判別正確率為75.0%。

圖7 驗證集中3D銅鉛判別特征點與判別規(guī)則面分布

經(jīng)2D、3D坐標系下訓(xùn)練集與驗證集的判別結(jié)果可知，通過光譜變換與CLI形成的CLDFP包含了其所受銅鉛元素種類相關(guān)的信息，使得EC-PB法建立的CLDRL2D與CLDRP3D能夠較好的區(qū)分銅鉛污染下的玉米葉片光譜數(shù)據(jù)，其中2D數(shù)據(jù)判別精度略高于3D數(shù)據(jù)；3D數(shù)據(jù)可從多個角度查看，判別效果更為直觀、易懂。

3 結(jié) 論

(1)玉米葉片光譜數(shù)據(jù)中的銅鉛元素污染信息難以通過原始光譜曲線進行準確判別，原始光譜數(shù)據(jù)與銅鉛元素種類之間的相關(guān)性較低；CR，SR以及FOD光譜變換處理均增加了玉米葉片光譜數(shù)據(jù)與銅鉛元素種類之間的相關(guān)性；不同階次FOD處理后的光譜數(shù)據(jù)與銅鉛元素種類之間的相關(guān)性各不相同。

(2)各階次FOD對應(yīng)的CLI與重金屬銅鉛元素種類相關(guān)系數(shù)各不相同，隨著階次的增加，相關(guān)性呈現(xiàn)先遞增，后遞減的趨勢，其中相關(guān)系數(shù)最高的三個階次分別為1.2階，0.7階，1.0階，其對應(yīng)的相關(guān)系數(shù)為-0.719，-0.701，0.700。

(3)基于EC-PB建立的2D、3D坐標系下的銅鉛元素判別規(guī)則線CLDRL2D與銅鉛元素判別規(guī)則面CLDRP3D能夠較好地區(qū)分玉米葉片光譜數(shù)據(jù)中的銅鉛污染元素種類，在2D坐標系下訓(xùn)練集判別正確率為78.95%，驗證集判別正確率為75.0%；在3D坐標系下訓(xùn)練集判別正確率為76.32%，驗證集判別正確率為75.0%。其中2D數(shù)據(jù)判別精度略高于3D數(shù)據(jù)；3D數(shù)據(jù)可從多個角度查看，判別效果更為直觀、易懂。

利用EC-PB能較為理想地取得玉米葉片污染的銅鉛元素判別與識別效果，同時也存在一定的局限性，如模型的普適性、除銅鉛外其他重金屬元素的污染識別以及野外大面積農(nóng)作物重金屬污染元素種類識別還需要做進一步的研究。