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

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

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

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

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

1 材料和方法

1.1 研究區(qū)概況

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

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

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

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

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

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

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

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

1.4 光譜指標提取

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

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

(1)

(2)

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

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

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

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

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

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

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

(3)

(4)

(5)

2 結果與分析

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

本研究區(qū)內41個采樣點的冬小麥重金屬含量描述性統(tǒng)計結果見表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)計結果

2.2 冬小麥冠層各光譜指標與6種重金屬含量的相關性分析

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

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

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

2.3 冬小麥葉片重金屬高光譜反演模型的構建

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

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

表3 冬小麥冠層各光譜指標與6種重金屬含量的建模結果比較

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

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

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

表4 冬小麥冠層各光譜指標與6種重金屬含量反演模型的驗證結果比較

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

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

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

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

在每種重金屬的最優(yōu)反演模型選出來后，對20個反演樣本的重金屬含量進行反演并根據(jù)反演出的重金屬含量以及實測的重金屬含量進行普通克里金插值(圖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 結論與討論

經污水灌溉的農作物會受到重金屬的污染，重金屬在植物葉片內的積累可通過微弱的光譜波動體現(xiàn)出來[32-33]。本研究利用冠層高光譜信息通過建立回歸模型對污灌區(qū)冬小麥的重金屬含量進行估算。結果表明， 對于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)；對于Cr，基于反射率二階微分的偏最小二乘回歸模型(SDR-PLSR)為最優(yōu)模型(R2=0.846，RPD=2.013)；對于Ni、Hg，基于光譜參數(shù)的偏最小二乘回歸模型(SP-PLSR)為最優(yōu)模型(Ni：R2=0.887，RPD=1.872；Hg：R2=0.819，RPD=1.684)。

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

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