山東典型污灌區(qū)冬小麥葉片重金屬高光譜反演及空間分布特征

于 慶，吳泉源，姚 磊，徐夕博，周 旭

(山東師范大學(xué) 地理與環(huán)境學(xué)院， 山東 濟(jì)南 250358)

在我國(guó)水資源相對(duì)匱乏的北方地區(qū)，有利用處理后的工業(yè)污水灌溉的傳統(tǒng)。20世紀(jì)60年代以來(lái)，北方地區(qū)的污水灌溉面積迅速擴(kuò)大，約占全國(guó)污水灌溉總面積的90%。由于污水中含有豐富的有機(jī)質(zhì)懸浮物，采用污水灌溉可以緩解水資源短缺、節(jié)省肥料、提高土壤肥力，然而長(zhǎng)期的污水灌溉必然會(huì)導(dǎo)致Pb、Cr、Cd、Hg等重金屬在土壤中累積[1]，最終導(dǎo)致污水灌溉的農(nóng)作物受到重金屬的污染[2-3]。冬小麥?zhǔn)潜狈降貐^(qū)的主要糧食作物，產(chǎn)量約占全國(guó)小麥總產(chǎn)量的56%，冬小麥中重金屬的富集關(guān)系食品安全和人體健康[4-7]。因此，在污水灌溉區(qū)進(jìn)行冬小麥重金屬監(jiān)測(cè)具有重要意義。

傳統(tǒng)的重金屬污染檢測(cè)是實(shí)地采樣后進(jìn)行實(shí)驗(yàn)室分析，具有精度高的特點(diǎn)，但工作量大、成本高，檢驗(yàn)點(diǎn)不連續(xù)，只能做到以點(diǎn)代面，不能推廣到大范圍，極大地影響了農(nóng)業(yè)決策的全面性、時(shí)效性和客觀性[8]。高光譜遙感技術(shù)使快速進(jìn)行農(nóng)作物重金屬無(wú)損監(jiān)測(cè)成為可能，具有監(jiān)測(cè)面積大、易推廣、工作量小的特點(diǎn)，尤其是冠層高光譜[9-12]，其不但可以獲取作物體內(nèi)生化組分信息，而且還能反映作物下伏土壤的理化性質(zhì)。目前，利用高光譜對(duì)農(nóng)田重金屬污染進(jìn)行監(jiān)測(cè)的研究很多[13-23]，但多是對(duì)作物下伏土壤重金屬的高光譜監(jiān)測(cè)，并未直接對(duì)作物體內(nèi)的重金屬含量進(jìn)行反演估算，而重金屬在作物中的累積會(huì)因各種因素有所差異，作物下伏土壤的重金屬含量并不能間接反映作物自身重金屬的含量。

本研究選擇山東省的典型污水灌溉區(qū)龍口市，在龍口市北部平原污水灌溉區(qū)設(shè)61個(gè)采樣點(diǎn)，野外實(shí)測(cè)污灌區(qū)冬小麥冠層反射光譜，結(jié)合小麥重金屬含量數(shù)據(jù)，運(yùn)用逐步多元線性回歸(Stepwise multiple linear regression,SMLR)以及偏最小二乘回歸(Partial least squares regression,PLSR)的方法建立反演模型，選取各種重金屬含量的最優(yōu)反演模型，并通過(guò)反演及插值，分析該污水灌溉區(qū)冬小麥重金屬含量的空間分布特征，以期為該污水灌溉區(qū)農(nóng)作物重金屬的實(shí)時(shí)監(jiān)測(cè)提供模型理論基礎(chǔ)，并為該污水灌溉區(qū)冬小麥重金屬脅迫診斷提供理論依據(jù)。

1 材料和方法

1.1 研究區(qū)概況

本研究的污灌區(qū)位于山東省龍口市(120°13′14″～120°44′46″E、37°27′ 30″～37°47′24″N)北部平原區(qū)，屬溫帶季風(fēng)性氣候，冬無(wú)嚴(yán)寒，夏無(wú)酷暑，氣候宜人，總面積409 km2，是膠東半島的主要糧食產(chǎn)區(qū)，農(nóng)業(yè)生產(chǎn)條件發(fā)達(dá)，主要種植冬小麥和夏玉米等作物。該區(qū)污水灌溉歷史已有30 a，鋁金屬冶煉、塑膠、鈦業(yè)化工等工廠排放的污水以及生活廢水經(jīng)龍口市污水處理廠處理后直接用于農(nóng)田的灌溉，使得土壤中重金屬含量不斷積累，重金屬主要為Cr、Ni、Pb、Zn、Hg和Cd 6種。

1.2 冬小麥冠層光譜及葉片采集

冬小麥冠層光譜的采集采用美國(guó)ASD (Analytical spectral device)公司的Field Spec HH便攜式手持地物光譜儀，采樣間隔為1 nm，光譜分辨率為3 nm，其波長(zhǎng)為325～1 075 nm，采樣時(shí)間為冬小麥的拔節(jié)期和灌溉最佳時(shí)期的4月中旬(2016年4月13—15日)，采樣時(shí)段是晴朗無(wú)風(fēng)的中午(11:00—14:00)。根據(jù)本研究區(qū)的土地利用現(xiàn)狀圖以及實(shí)地勘察，在本研究區(qū)的冬小麥種植區(qū)隨機(jī)設(shè)立61個(gè)樣點(diǎn)(圖1)，采集61個(gè)樣點(diǎn)處的冬小麥冠層光譜數(shù)據(jù)并隨機(jī)采集41處冬小麥葉片，將其帶回實(shí)驗(yàn)室化驗(yàn)重金屬含量。樣點(diǎn)坐標(biāo)采用GPS載波相位差分技術(shù)RTK方法進(jìn)行定位。

1.3 光譜數(shù)據(jù)預(yù)處理和重金屬含量測(cè)定

在進(jìn)行冠層光譜測(cè)量前，先利用白板進(jìn)行校正，將光譜儀器探頭垂直向下，距離冬小麥冠層頂部0.5 m，視角為8°，光譜測(cè)量時(shí)，每個(gè)樣點(diǎn)重復(fù)測(cè)量10次。在ViewSpecPro軟件中將每個(gè)樣點(diǎn)重復(fù)測(cè)量的10條光譜中與其他光譜曲線有明顯區(qū)別的曲線去除掉，再將剩下的光譜曲線取平均值得到該采樣點(diǎn)實(shí)際光譜的反射率。為了從野外高光譜數(shù)據(jù)中提取有效光譜信息，需對(duì)野外獲得光譜信息進(jìn)行平滑處理，并用9點(diǎn)加權(quán)移動(dòng)平均法去除噪聲，原始高光譜經(jīng)處理后可以消除背景噪聲，增強(qiáng)差別，突出光譜特征[24]，見(jiàn)圖2。

圖1 污灌區(qū)冬小麥冠層光譜數(shù)據(jù)采樣點(diǎn)分布情況

圖2 去除噪聲以及平滑處理后的光譜

光譜采集后，將采樣點(diǎn)處的冬小麥連根拔起，裝入樣品袋中帶回實(shí)驗(yàn)室分析。將帶回實(shí)驗(yàn)室的小麥樣品使用去離子水洗凈后放入80 ℃烘箱中烘干，將葉子剪下，放入粉碎機(jī)進(jìn)行粉碎，冷卻后用1∶1的HNO3溶解灰樣，蒸干，用0.1 mol/L的HNO3定容，將定容完成的植物樣品溶液用氫化物發(fā)生-原子熒光光譜法(HG-AFS)檢測(cè)Hg含量，用電感耦合等離子體質(zhì)譜法(ICP-MS)檢測(cè)Cd含量，用電感耦合等離子體原子發(fā)射光譜法(ICP-OES)檢測(cè)pb、Zn、Ni、Cr含量。

1.4 光譜指標(biāo)提取

原始光譜反射率(Raw spectra reflectance,R)與各種重金屬含量的相關(guān)性較差，對(duì)重金屬含量信息的反映很微弱，而低階微分處理能平緩背景干擾產(chǎn)生的影響，使光譜數(shù)據(jù)的輪廓變得更加清晰，提高原始光譜分辨率并且具有放大微小峰值的效果[25]。光譜參數(shù)(Spectral parameters,SP)也廣泛應(yīng)用于植被光譜建模反演研究中，它可以綜合多個(gè)敏感波段的特征，加強(qiáng)對(duì)信息的提取能力，避免單一波段的偶然性和不精確性[26]。因此，本研究將原始光譜反射率、反射率一階微分 (First derivate reflectance,FDR)、反射率二階微分(Second derivate reflectance,SDR)和光譜參數(shù)分別作為反演指標(biāo)與原始反射光譜進(jìn)行建模對(duì)比，以期得到高精度的反演模型。

1.4.1 低階微分指標(biāo) 將原始光譜反射率進(jìn)行低階微分處理，反射率一階微分、反射率二階微分見(jiàn)圖3。反射率一、二階微分計(jì)算公式如公式(1)和公式(2)所示：

(1)

(2)

其中，R(λi+1)為第i+1條波段的反射率，R(λi-1)為第i-1條波段的反射率，Δλ為光譜的波段間隔。

1.4.2 光譜參數(shù) 植物在可見(jiàn)光波段具有很強(qiáng)的光吸收能力，根據(jù)不同波長(zhǎng)的光吸收強(qiáng)弱變化形成波峰和波谷，因此，可見(jiàn)光波段是研究的重點(diǎn)光譜區(qū)域[27-29]。本研究在可見(jiàn)光區(qū)域內(nèi)提取光譜參數(shù)，提取的光譜參數(shù)有各種光譜位置參數(shù)、光譜面積參數(shù)。各種參數(shù)的計(jì)算方法見(jiàn)表1。

1.5 冬小麥葉片重金屬高光譜反演的方法

偏最小二乘回歸法是一種多對(duì)多的回歸建模方法，比較適合處理變量多而樣本數(shù)少的問(wèn)題[30]，而逐步多元線性回歸法是常用的統(tǒng)計(jì)建模方法，可以解決數(shù)據(jù)維數(shù)帶來(lái)的波段冗余和高光譜數(shù)據(jù)繁多的問(wèn)題[31]，本研究樣本較小且光譜數(shù)據(jù)冗余，因此，選擇偏最小二乘回歸法及逐步多元線性回歸法進(jìn)行建模，并對(duì)比分析，旨在選出各種重金屬的最優(yōu)反演模型。本研究在對(duì)4種光譜指標(biāo)與6種重金屬含量進(jìn)行相關(guān)性分析的基礎(chǔ)上，分別以各種光譜指標(biāo)為自變量，重金屬含量為因變量進(jìn)行建模分析，建模過(guò)程如下。

圖3 原始光譜反射率微分變換

類型參數(shù)符號(hào)計(jì)算方法光譜位置參數(shù)紅邊位置λr紅邊內(nèi)(680～750 nm)最大一階微分對(duì)應(yīng)的波長(zhǎng)位置紅邊峰值Dr紅邊內(nèi)(680～750 nm)光譜一階微分的最大值紅谷深度Ro紅邊內(nèi)(680～750 nm)包絡(luò)線去除后光譜的最小值黃邊位置λy黃邊內(nèi)(550～650 nm)最小一階微分值對(duì)應(yīng)的波長(zhǎng)位置黃邊峰值Dy黃邊內(nèi)(550～650 nm)光譜一階微分的最小值藍(lán)谷位置λb藍(lán)邊內(nèi)(490～530 nm)最大一階微分值對(duì)應(yīng)的波長(zhǎng)位置藍(lán)谷深度Db藍(lán)邊內(nèi)(490～530 nm)光譜一階微分的最大值綠峰位置λg510～580 nm波段內(nèi)最大的波段發(fā)射率所在位置綠峰峰值Rg510～580 nm波段內(nèi)最大的波段發(fā)射率光譜面積參數(shù)紅邊面積SDr紅邊內(nèi)(680～750 nm)一階微分值的總和藍(lán)邊面積SDb藍(lán)邊內(nèi)(490～530 nm)一階微分值的總和黃邊面積SDy黃邊內(nèi)(550～650 nm)一階微分值的總和

將獲取的61份實(shí)測(cè)數(shù)據(jù)分為兩部分，一部分是測(cè)得冠層光譜并測(cè)得對(duì)應(yīng)冬小麥各種重金屬含量的41份樣本數(shù)據(jù)，這部分?jǐn)?shù)據(jù)進(jìn)行建模與驗(yàn)證，其中隨機(jī)選擇31個(gè)樣本數(shù)據(jù)用于建立反演模型，10個(gè)樣本為檢測(cè)樣本用于檢驗(yàn)?zāi)Ｐ?；另一部分剩余?0份只采集光譜數(shù)據(jù)而未測(cè)重金屬含量的樣本用來(lái)反演。采用模型決定系數(shù)(R2)、校正集的均方根誤差(RMSEC)、檢驗(yàn)集均方根誤差(RMSEV)以及相對(duì)分析誤差(RPD)、相對(duì)誤差(RE)來(lái)檢驗(yàn)?zāi)Ｐ途?，?jì)算公式如式(3)—(5)所示。其中，R2越大，RMSEC越小，建模效果越好；RMSEV、RE越小，RPD越大，模型的預(yù)測(cè)精度越高。

(3)

(4)

(5)

2 結(jié)果與分析

2.1 山東典型污灌區(qū)冬小麥葉片重金屬含量分析

本研究區(qū)內(nèi)41個(gè)采樣點(diǎn)的冬小麥重金屬含量描述性統(tǒng)計(jì)結(jié)果見(jiàn)表2，Cr、Ni、Pb、Zn、Hg、Cd含量的平均值分別為8.768、3.663、4.279、35.493、0.011、0.149 mg/kg；變異系數(shù)反映數(shù)據(jù)離散程度，重金屬含量變異系數(shù)表現(xiàn)為Ni>Cr>Pb>Cd>Hg>Zn，表明本研究區(qū)Cr、Ni含量相較Hg、Pb、Zn、Cd含量離散程度大，地區(qū)間分布不均。

表2 山東典型污灌區(qū)冬小麥葉片重金屬含量的描述性統(tǒng)計(jì)結(jié)果

2.2 冬小麥冠層各光譜指標(biāo)與6種重金屬含量的相關(guān)性分析

由圖4可知，原始光譜反射率對(duì)不同重金屬含量的波譜響應(yīng)不明顯，相關(guān)性曲線相對(duì)平緩，相關(guān)系數(shù)較低，基本在-0.25～0.25，重金屬Cr、Ni、Pb、Zn含量與700～900 nm波段相關(guān)系數(shù)較大，相關(guān)性較強(qiáng)。原始光譜反射率與重金屬Hg含量呈負(fù)相關(guān)，在380～700 nm的可見(jiàn)光范圍內(nèi)的相關(guān)性很強(qiáng)，相關(guān)系數(shù)為-0.45左右，而原始光譜反射率與重金屬Cd含量相關(guān)系數(shù)的高值區(qū)在750～1 100 nm，但總體上原始光譜反射率與重金屬含量的相關(guān)性較低，對(duì)重金屬含量信息的反映很微弱。

在經(jīng)過(guò)一階、二階微分變換后，各波段與重金屬含量的相關(guān)性明顯增強(qiáng)，光譜響應(yīng)劇烈，相關(guān)系數(shù)的變化幅度很大，在正負(fù)之間來(lái)回波動(dòng)，相關(guān)系數(shù)的絕對(duì)值基本都在0.35以上，重金屬Hg含量與變換后各波段的相關(guān)性較高，相關(guān)系數(shù)的絕對(duì)值在0.40～0.75。光譜參數(shù)與各種重金屬含量的相關(guān)系數(shù)絕對(duì)值基本在0.25～0.40(圖5)，相關(guān)性大部分都達(dá)到0.05顯著性水平，其中重金屬Hg含量與黃邊面積、紅邊面積的相關(guān)系數(shù)分別為0.78、-0.75，相關(guān)性較強(qiáng)，達(dá)到了極顯著水平(P<0.01)。

圖4 冬小麥冠層原始光譜反射率、反射率一階微分以及反射率二階微分與重金屬含量的相關(guān)系數(shù)

2.3 冬小麥葉片重金屬高光譜反演模型的構(gòu)建

由于光譜波段數(shù)目繁多，本研究對(duì)原始光譜反射率、反射率一階微分、反射率二階微分進(jìn)行重采樣，選出與各種重金屬含量相關(guān)系數(shù)大于0.35且呈顯著相關(guān)的波段進(jìn)行建模。建模過(guò)程中，以光譜參數(shù)及重采樣的原始光譜反射率、反射率一階微分、反射率二階微分作為模型的自變量，以各種冬小麥葉片重金屬含量作為因變量，建模比較結(jié)果見(jiàn)表3。從表3可以看出，利用偏最小二乘回歸法建立的回歸模型的R2都達(dá)到0.70以上，建模穩(wěn)定性較高，建模效果明顯優(yōu)于逐步多元線性回歸法，尤其對(duì)于光譜參數(shù)來(lái)說(shuō)，建立的偏最小二乘回歸模型的穩(wěn)定性遠(yuǎn)遠(yuǎn)大于建立的逐步多元線性回歸模型，誤差也小得多。相比原始光譜反射率建模，利用反射率低階微分建立的模型的R2總體明顯增大，RMSEC較小。對(duì)于重金屬Cr，SDR-PLSR模型效果最優(yōu)，R2可達(dá)到0.846，RMSEC較?。粚?duì)于重金屬Ni，SP-PLSR模型的R2最高，達(dá)0.887，RMSEC最小，為1.313，建模效果很好；對(duì)于重金屬Pb，F(xiàn)DR-PLSR建模效果最優(yōu)，R2為0.848，RMSEC較小，為1.964；對(duì)于重金屬Zn，F(xiàn)DR-PLSR模型的R2達(dá)0.790，RMSEC較小，為5.139，建模效果最優(yōu)；對(duì)于重金屬Hg，SP-PLSR模型的R2最高，為0.819，RMSEC較小，為0.002，建模效果很好；對(duì)于重金屬Cd，F(xiàn)DR-PLSR建模效果最優(yōu)，R2可達(dá)到0.868，模型穩(wěn)定性遠(yuǎn)遠(yuǎn)高于其他模型的穩(wěn)定性，RMSEC為0.085，與其他指標(biāo)差距較小。

圖5 冬小麥冠層光譜參數(shù)與重金屬含量的相關(guān)系數(shù)

表3 冬小麥冠層各光譜指標(biāo)與6種重金屬含量的建模結(jié)果比較

2.4 冬小麥葉片重金屬高光譜反演模型的檢驗(yàn)與最優(yōu)模型的選取

對(duì)回歸模型進(jìn)行檢驗(yàn)(表4)發(fā)現(xiàn)，對(duì)于各種重金屬，利用逐步多元線性回歸法建立的反演模型的RMSEV較大，RPD較小,在0.509～1.057，RE較大，說(shuō)明對(duì)重金屬含量的預(yù)測(cè)能力一般；而偏最小二乘回歸模型的RPD為0.717～2.406，RMSEV、RE較小，說(shuō)明偏最小二乘回歸模型對(duì)重金屬含量的預(yù)測(cè)誤差較小，預(yù)測(cè)能力較好。對(duì)于重金屬Cr，SDR-PLSR模型的RMSEV為1.136，RPD為2.013，RE為26.711%，與其他指標(biāo)相比，模型預(yù)測(cè)效果最優(yōu)；對(duì)于重金屬Ni，SP-PLSR模型的RMSEV最小，RPD和RE分別為1.872和22.277%，相比其他指標(biāo)，RE最小，RPD最大，模型預(yù)測(cè)效果最好；對(duì)于重金屬Hg，SP-PLSR模型的RPD為1.684，RE為8.182%，模型預(yù)測(cè)效果最好；而對(duì)于重金屬Pb、Zn、Cd，F(xiàn)DR-PLSR模型的RMSEV最小，RE也最小，RPD最大，模型預(yù)測(cè)效果最優(yōu)。

綜合建模與驗(yàn)證結(jié)果可以看出，重金屬Cr含量的最優(yōu)反演模型為SDR-PLSR模型，SP-PLSR模型是重金屬Ni、Hg含量的最優(yōu)反演模型，F(xiàn)DR-PLSR模型是重金屬Pb、Zn和Cd含量的最優(yōu)反演模型。

表4 冬小麥冠層各光譜指標(biāo)與6種重金屬含量反演模型的驗(yàn)證結(jié)果比較

續(xù)表4 冬小麥冠層各光譜指標(biāo)與6種重金屬含量反演模型的驗(yàn)證結(jié)果比較

為了更加直觀地表現(xiàn)各種重金屬的最優(yōu)模型擬合結(jié)果，繪制6種重金屬含量實(shí)測(cè)值和反演模型預(yù)測(cè)值的散點(diǎn)圖，如圖6所示。每種重金屬含量的實(shí)測(cè)值與預(yù)測(cè)值的擬合曲線都在y=x附近，所選取的最優(yōu)反演模型的擬合效果都較好。

圖6 山東污灌區(qū)冬小麥葉片重金屬含量實(shí)測(cè)值和最優(yōu)反演模型預(yù)測(cè)值的擬合情況

2.5 山東典型污水灌溉區(qū)重金屬含量空間分布特征

在每種重金屬的最優(yōu)反演模型選出來(lái)后，對(duì)20個(gè)反演樣本的重金屬含量進(jìn)行反演并根據(jù)反演出的重金屬含量以及實(shí)測(cè)的重金屬含量進(jìn)行普通克里金插值(圖7)。由圖7可知，冬小麥葉片中Cr含量為2.579～30.237 mg/kg，Ni含量為0.730～15.374 mg/kg,Cr、Ni含量的高值區(qū)都位于研究區(qū)的東南部，北部及西北部含量較低；Pb、Zn含量分布的高值區(qū)主要位于中部及南部地區(qū)，北部含量較低；Hg含量為0.002 5～0.021 8 mg/kg,在西北部較低；Cd含量為0.062～0.478 mg/kg，在中部、北部以及西北部較低，其他位置較高。

圖7 山東典型污水灌溉區(qū)冬小麥葉片重金屬的空間分布

3 結(jié)論與討論

經(jīng)污水灌溉的農(nóng)作物會(huì)受到重金屬的污染，重金屬在植物葉片內(nèi)的積累可通過(guò)微弱的光譜波動(dòng)體現(xiàn)出來(lái)[32-33]。本研究利用冠層高光譜信息通過(guò)建立回歸模型對(duì)污灌區(qū)冬小麥的重金屬含量進(jìn)行估算。結(jié)果表明， 對(duì)于Pb、Zn、Cd，基于反射率一階微分的偏最小二乘回歸模型(FDR-PLSR)為最優(yōu)模型(Pb：R2=0.848，RPD=1.598；Zn：R2=0.790，RPD=2.295；Cd：R2=0.868，RPD=2.406)；對(duì)于Cr，基于反射率二階微分的偏最小二乘回歸模型(SDR-PLSR)為最優(yōu)模型(R2=0.846，RPD=2.013)；對(duì)于Ni、Hg，基于光譜參數(shù)的偏最小二乘回歸模型(SP-PLSR)為最優(yōu)模型(Ni：R2=0.887，RPD=1.872；Hg：R2=0.819，RPD=1.684)。

從空間插值結(jié)果可以看出，冬小麥葉片中Cr、Ni含量在研究區(qū)東南部較高，北部及西北部較低；Pb、Zn含量在中部以及南部較高；Hg含量在西北部較低；Cd含量在中部、北部、西北部較低，其他區(qū)域較高。研究表明，選取的最優(yōu)模型能較精確地對(duì)污水灌溉區(qū)冬小麥葉片重金屬進(jìn)行估算，可實(shí)時(shí)且大面積地進(jìn)行冬小麥葉片重金屬污染監(jiān)測(cè)。

本研究所用的是野外實(shí)測(cè)冠層光譜數(shù)據(jù)，在采樣過(guò)程中，土壤背景對(duì)冠層光譜的影響不可避免，并且由于是區(qū)域?qū)嵉乇O(jiān)測(cè)，不同區(qū)域土壤質(zhì)地、水分狀況等因素差異較大，可能對(duì)光譜產(chǎn)生一定影響。在下一步的研究中，將研究土壤質(zhì)地以及含水量對(duì)冠層光譜的影響，并結(jié)合高空高光譜影像在更大范圍更加精確地進(jìn)行重金屬污染監(jiān)測(cè)。