蔡太義，王志剛，楊劉帥，2，王群，黃會娟，于海洋，張傳忠，張燦，劉鵬，馮玉慶，賀成龍，張合兵*

（1.河南理工大學(xué)測繪與國土信息工程學(xué)院，河南 焦作 454000；2.自然資源部第四地形測量隊(duì)，哈爾濱 150000；3.河南農(nóng)業(yè)大學(xué)農(nóng)學(xué)院/河南省農(nóng)田生態(tài)系統(tǒng)野外科學(xué)觀測研究站，鄭州 450046；4.河南省土壤調(diào)理與修復(fù)工程技術(shù)研究中心，河南 商丘 476000；5.中向旭曜科技有限公司，江蘇 無錫 214400；6.河南省地質(zhì)局礦產(chǎn)資源勘查中心,鄭州 450053）

土壤重金屬是土壤復(fù)雜生態(tài)系統(tǒng)中重要的組成成分，影響著土壤環(huán)境質(zhì)量[1]。鋅（Zn）是一種重要的土壤重金屬元素[2]，是在大氣沉降、工業(yè)污染和施肥等過程中緩慢積累，且大部分處于結(jié)合狀態(tài)的金屬離子，并且其不易降解[3]。我國土壤中Zn的環(huán)境背景含量在2.60～593.00 mg·kg-1之間[4]。Zn含量過低會減少作物激素的生成，過高則會對作物產(chǎn)生毒害作用[5]。因此，Zn含量的研究是高標(biāo)準(zhǔn)農(nóng)田建設(shè)過程中的一個(gè)重要內(nèi)容。傳統(tǒng)重金屬測量與監(jiān)測方法具有耗時(shí)長、范圍窄和實(shí)時(shí)性差的缺點(diǎn)，不能滿足反映生態(tài)系統(tǒng)結(jié)構(gòu)變化的要求[6]。高光譜遙感數(shù)據(jù)可以快速檢測土壤重金屬屬性信息[7]，能節(jié)省人力與物力，同時(shí)具有成本低、無污染和無損的優(yōu)勢[8]，是土壤重金屬含量研究的新趨勢[9]。

前人對Zn元素含量研究已取得重大進(jìn)展。陶超等[10]利用定量回歸模型和定性反演模型比較采樣點(diǎn)重金屬鉛（Pb）與Zn的可遷移能力，證明在快速檢測土壤重金屬污染狀況的問題上，定性分類是一種更加切實(shí)可行的方式。王金鳳等[11]運(yùn)用電感耦合等離子體質(zhì)譜和分光光度計(jì)測定土壤樣品，分析土壤光譜和Zn含量的映射關(guān)系，實(shí)現(xiàn)了Zn含量高效快速反演，一定程度上揭示了喀斯特地區(qū)土壤重金屬Zn的光譜吸收特性。程先鋒等[12]通過構(gòu)建土壤中Zn、Pb、砷（As）、鎘（Cd）4種重金屬含量及土壤光譜之間的多元線性回歸模型，證明土壤重金屬含量與光譜特征具有較強(qiáng)的相關(guān)性。宋婷婷等[13]建立了基于乘積變換的Zn含量最優(yōu)預(yù)測模型并制作了污染分布圖，表明Zn污染和礦業(yè)活動(dòng)間具有顯著相關(guān)性。TAN等[14]在伊通縣礦區(qū)對選定的光譜特征和土壤樣品的重金屬數(shù)據(jù)進(jìn)行反演并建立反演模型，結(jié)果表明CARS疊加法能較好地預(yù)測研究區(qū)As、鉻（Cr）、Pb和Zn的含量。LIN等[15]利用黑龍江省黑土樣本的室內(nèi)光譜數(shù)據(jù)，分析了土壤有機(jī)質(zhì)、鐵（Fe）、Zn含量與6種光譜反射率的相關(guān)性；然后在α=0.01水平進(jìn)行相關(guān)系數(shù)顯著性檢驗(yàn)，采用主成分分析（KPCA）算法，結(jié)合隨機(jī)森林（RF）和支持向量機(jī)（SVM）組合模型應(yīng)用于Zn、Fe含量的定量反演。結(jié)果表明在土壤重金屬Zn、Fe的預(yù)測中，KPCA-RF模型優(yōu)于RF模型，其在土壤成分含量的定量反演中具有廣闊的應(yīng)用前景。

目前，對于土壤Zn的研究，國內(nèi)外多側(cè)重于礦區(qū)Zn的污染及土壤Zn的回歸模型的精度與遷移能力，而較少關(guān)注耕地質(zhì)量提升與治理中Zn的反演研究，從而制約了國家土壤健康及高標(biāo)準(zhǔn)農(nóng)田建設(shè)的穩(wěn)步推進(jìn)。因此，本研究以黃淮海典型區(qū)域（西平縣）高標(biāo)準(zhǔn)農(nóng)田土壤重金屬Zn作為研究對象，于實(shí)驗(yàn)室測得野外采集土樣的Zn含量，并利用ASD Field Spec便攜式光譜儀采集土壤光譜數(shù)據(jù)；運(yùn)用5種方式對反射率光譜進(jìn)行變換以增強(qiáng)光譜特征信息，并與土壤Zn含量實(shí)測值進(jìn)行皮爾遜（Pearson）相關(guān)分析，初次篩選出顯著特征波段；再利用連續(xù)投影算法（SPA）二次篩選出最佳特征波段，以構(gòu)建Zn含量偏最小二乘法回歸模型（PLS）；檢驗(yàn)5個(gè)模型的均方根誤差（RMSEC、RMSECV、RMSEP）、校正決定系數(shù)（R2）和相對分析誤差（RPD）并進(jìn)行比較，選取Zn含量最佳反演模型。本研究旨在為預(yù)測高標(biāo)準(zhǔn)農(nóng)田建設(shè)區(qū)域土壤重金屬Zn含量及其治理提供理論參考。

1 材料與方法

1.1 研究區(qū)概況

駐馬店市西平縣位于北緯33°10'～33°32'，東經(jīng)113°36'～114°13'（圖1），河南中南部區(qū)域，其土壤類型分布及面積占比見表1?？h境總面積為1 089.77 km2，劃分為3個(gè)街道辦事處、8鎮(zhèn)、8鄉(xiāng)和老王坡管理委員會等行政區(qū)。該地氣候?qū)賮啙駶櫞箨懶约撅L(fēng)型，全年氣候溫和。西平縣資源豐富，農(nóng)作物物種豐盛，是我國糧食生產(chǎn)先進(jìn)縣與糧食主產(chǎn)區(qū)之一，同時(shí)其西部山區(qū)擁有森林覆蓋率95%以上的國家級森林公園。

圖1 研究區(qū)位置示意圖Figure 1 Location diagram of the study area

表1 西平縣土壤類型分布及面積占比Table 1 Distribution of soil types and area share in Xiping County

1.2 土樣采集

土樣采集前對西平縣研究區(qū)的土壤類型、空間特征、行政單元和地形因素進(jìn)行分析，以設(shè)定精確定位的采樣點(diǎn)。利用對角線采樣法采樣，野外采樣深度為0～30 cm土壤層，同時(shí)記錄詳細(xì)地形特征。從田里采集的土樣及時(shí)攤開風(fēng)干，倒在干凈的容器內(nèi)，去除石塊、植物根莖殘?bào)w及其余雜物。取研磨均勻后通過1 mm孔徑孔篩的土壤樣本，利用四分法取樣，用于土壤光譜數(shù)據(jù)的采集和Zn含量的測定[16]。

1.3 土樣光譜室內(nèi)測定

1.4 光譜變換方法

倒數(shù)對數(shù)光譜變換可以增強(qiáng)光譜反射率數(shù)據(jù)在可見光區(qū)域的特征差異，并能減少因光照強(qiáng)度和地形起伏變化等因素對光譜特征曲線引起的負(fù)面影響[18]，其變換公式見表2。

光譜微分技術(shù)是一種簡單而有效的確定光譜曲線極值點(diǎn)波長位置的方法，是通過對反射率光譜數(shù)據(jù)進(jìn)行階數(shù)不同的數(shù)學(xué)微分模擬來提取特征光譜參數(shù)，該技術(shù)能夠消除平緩背景部分環(huán)境或基線漂移的影響，擴(kuò)大樣本之間的光譜特征差異，有助于吸收特征的提取[19]（表2）。

包絡(luò)線是利用直線連接每個(gè)隨波長不斷變化的吸收谷或反射峰凸出的特征點(diǎn)，并使連接線的外角大于180°的線[20]。去包絡(luò)線是以相應(yīng)波段上實(shí)際光譜與包絡(luò)線兩者反射率值進(jìn)行比值計(jì)算，從而使之歸一化[21]（表2）。

表2 光譜變換公式Table 2 Spectral transformation equations

1.5 相關(guān)分析與顯著性檢驗(yàn)

相關(guān)分析是用于判斷多個(gè)隨機(jī)變量間是否存在統(tǒng)計(jì)學(xué)關(guān)聯(lián)的有效統(tǒng)計(jì)分析方法[22]，且能深入分析存在關(guān)聯(lián)的變量間的關(guān)聯(lián)強(qiáng)度與方向，反映事物間的統(tǒng)計(jì)學(xué)相關(guān)程度。兩個(gè)變量之間相關(guān)關(guān)系的強(qiáng)弱程度，可以用數(shù)值的形式揭示，其數(shù)值又稱為相關(guān)系數(shù)。皮爾遜相關(guān)系數(shù)（積差相關(guān)系數(shù)）是應(yīng)用最廣泛的相關(guān)系數(shù)，其公式如下：

式中：xi、yi（i=1，2…，n）為兩個(gè)隨機(jī)變量x和y的樣本值；xˉ、yˉ為兩個(gè)隨機(jī)變量樣本值的平均值；n為兩個(gè)變量樣本值的個(gè)數(shù)。R為相關(guān)系數(shù)，且-1≤R≤1；||R越接近1，說明兩變量間線性相關(guān)性越強(qiáng)，||R=1，為完全線性相關(guān)越接近0，說明兩變量間線性相關(guān)性越弱，為完全不具有相關(guān)性；R＞0，說明兩變量為正相關(guān)；R＜0，說明兩變量為負(fù)相關(guān)。在科研應(yīng)用中，可以劃分幾個(gè)區(qū)間取值來表示不同的相關(guān)程度。

樣本間是否存在顯著線性相關(guān)關(guān)系不能直接通過相關(guān)性分析得到的系數(shù)值確定，其原因在于隨機(jī)抽樣的性質(zhì)或樣本數(shù)量的不足等，因此必須使用假設(shè)檢驗(yàn)的方式對總體的顯著性進(jìn)行相關(guān)分析顯著性檢驗(yàn)[23]。

1.6 連續(xù)投影算法

連續(xù)投影算法（SPA）是將向矢量空間共線性最小化的前向變量選擇算法[24]。SPA通過投影將光譜數(shù)據(jù)的波長映射到其他波長上進(jìn)行向量的投影分析，從而比較出所投影向量的大小，選擇投影向量值最大的波長作為優(yōu)選波長，然后經(jīng)過矯正模型選取出所需要的特征波長[25]。

1.7 偏最小二乘法回歸模型反演與精度檢驗(yàn)方法

在構(gòu)建偏最小二乘法回歸模型（PLS）[26]中，從自變量中提取的主成分要求盡可能體現(xiàn)出光譜曲線的原特征信息[27]，另外從因變量中提取的獨(dú)立成分要求與主成分存在最大的相關(guān)性。將提取出的成分與因變量通過多元回歸方法構(gòu)建回歸方程[28]。RMSE、R2和PRD[29]的公式及定義見表3。

陳樹華所說的黃龍病，由一種限于韌皮部內(nèi)寄生的革蘭氏陰性細(xì)菌，是世界柑橙生產(chǎn)上的毀滅性病害?！包S龍病雖不可治，但是可防可控?！标悩淙A介紹。

表3 精度校驗(yàn)方法公式Table 3 Precision calibration method formula

2 結(jié)果與分析

2.1 土壤Zn屬性特征

土壤重金屬Zn含量范圍在57.13～306.40 mg·kg-1之間（表4），變異系數(shù)[30]為0.401，為中等變異性，適宜構(gòu)建模型。Zn元素的質(zhì)量分?jǐn)?shù)偏態(tài)分布密度曲線表現(xiàn)出土壤Zn含量多集中于125～130 mg·kg-1之間（圖2），極端值分布較少?！锻寥拉h(huán)境質(zhì)量 農(nóng)用地土壤污染風(fēng)險(xiǎn)管控標(biāo)準(zhǔn)（試行）》（GB 15618—2018）規(guī)定[31]，當(dāng)土壤pH≤6.5時(shí)，土壤重金屬Zn含量的風(fēng)險(xiǎn)篩選值為200 mg·kg-1。調(diào)查研究區(qū)土壤平均pH值為5.9，樣本數(shù)據(jù)表明該區(qū)整體Zn污染風(fēng)險(xiǎn)較低，少數(shù)樣本極端值表明該區(qū)局部Zn污染風(fēng)險(xiǎn)較高，應(yīng)當(dāng)加強(qiáng)土壤環(huán)境監(jiān)測和農(nóng)產(chǎn)品協(xié)同檢測。

圖2 Zn質(zhì)量分?jǐn)?shù)描述統(tǒng)計(jì)圖Figure 2 Descriptive statistics of Zn mass fraction

表4 土壤重金屬Zn含量分析Table 4 Analysis of soil properties for heavy metal Zn content

2.2 土壤屬性光譜反射特征

土壤光譜曲線一般為突出弧線，土壤樣本間的光譜差異性主要體現(xiàn)為相同波段上的光譜反射率數(shù)值不同（圖3）。土壤樣本在可見光帶400～780 nm范圍處的反射率較低，但呈大幅度上升的趨勢，光譜曲線陡峭，土壤有機(jī)質(zhì)在600、800 nm附近有微弱的反射峰，在780～2 300 nm處近紅外波段區(qū)域曲線較平緩，在2 300～2 400 nm處逐漸下降。1 850～2 140 nm光譜曲線呈現(xiàn)明顯的“V”字形。圖中3個(gè)水分吸收帶表現(xiàn)出的波谷與黏土中的多種礦物質(zhì)組成有關(guān)，1 420 nm為羥基（OH）譜帶、1 920 nm附近為H2O譜帶、2 200 nm波段范圍為羥基伸縮振動(dòng)與和OH彎曲振動(dòng)的合譜帶[32]。4條平均光譜曲線整體趨勢相似，丘崗地區(qū)的平均光譜曲線較低，山地丘陵區(qū)的平均光譜曲線最高，而縣區(qū)與平原區(qū)的平均光譜曲線居中，幾乎吻合。

圖3 西平縣不同地類土壤平均光譜反射率曲線Figure 3 Average spectral reflectance curves of different land types of soil in Xiping County

2.3 相關(guān)分析與顯著性檢驗(yàn)

從5種光譜變換后的光譜數(shù)據(jù)與土壤Zn含量的相關(guān)系數(shù)折線圖（圖4）可以看出，經(jīng)過光譜變換，各波段與土壤Zn含量的相關(guān)系數(shù)也有所不同，呈現(xiàn)多元化的變化趨勢。倒數(shù)對數(shù)光譜變換（LOG）與平滑光譜變換（SG）的相關(guān)系數(shù)曲線較為平滑，且呈現(xiàn)以數(shù)值0為對稱軸的對稱顯示，顯著相關(guān)系數(shù)的波段數(shù)、最大相關(guān)波段基本一致；二者中一者遞增，另一者則遞減，最顯著波段相關(guān)系數(shù)呈現(xiàn)相反關(guān)系，其絕對值十分相近，并均接近0.3。

圖4 5種變換光譜反射率相關(guān)系數(shù)顯著性曲線Figure 4 Significance curves for the correlation coefficients of the five transformed spectral reflectances

去包絡(luò)線光譜變換（CR）的相關(guān)系數(shù)曲線出現(xiàn)上下頻繁波動(dòng)，在600～1 400、1 900～2 400 nm范圍內(nèi)有多處波峰和波谷，數(shù)值變化范圍介于-0.3～0.3，曲線變化形似于一階微分變換（FD）后的光譜相關(guān)系數(shù)曲線。SG光譜反射率相關(guān)系數(shù)曲線與CR光譜相關(guān)系數(shù)曲線擁有相似的波峰位置。FD與二階微分變換（SD）光譜的相關(guān)系數(shù)數(shù)值變化折線圖不再呈現(xiàn)近似單一數(shù)值的平滑線，而是表現(xiàn)為在正負(fù)值上下頻繁的波動(dòng)，F(xiàn)D光譜的相關(guān)系數(shù)曲線在580 nm附近出現(xiàn)一個(gè)波峰，在1 900 nm附近出現(xiàn)一個(gè)波谷，SD的相關(guān)系數(shù)曲線變化較為劇烈，連續(xù)波段相鄰的相關(guān)系數(shù)在0上下跳躍。

FD、SD、CR變換光譜與Zn含量顯著相關(guān)的波段數(shù)量較SG交換光譜的顯著相關(guān)波段數(shù)明顯減少，其中SD變換光譜相關(guān)系數(shù)曲線的顯著相關(guān)波段數(shù)最少。5種變換光譜對應(yīng)的最顯著波段中，SD變換光譜的相關(guān)系數(shù)最大，F(xiàn)D變換和CR變換光譜的相關(guān)系數(shù)次之，LOG變換和SG變換光譜對應(yīng)的相關(guān)系數(shù)最小（表5）。

表5 光譜變換后反射率與Zn含量顯著相關(guān)的波段Table 5 The wavelengths with significant correlation between reflectivity and Zn content after spectral transformation

綜上所述，經(jīng)過LOG、FD、SG等5種光譜變換后的光譜反射率能有效突出光譜反射率的變化特征，經(jīng)檢驗(yàn)光譜變換下通過α=0.01水平顯著性檢驗(yàn)的相關(guān)波段可作為第一次特征波段數(shù)據(jù)，用于進(jìn)行第二次特征波段的篩選，從而構(gòu)建土壤Zn含量的反演模型。

2.4 最佳特征波段的選取

利用Matlab軟件結(jié)合連續(xù)投影算[33]法分別對SG、LOG、FD、SD和CR變換后的光譜數(shù)據(jù)進(jìn)行波段的第二次篩選，得出最佳特征波段見表6。

表6 土壤Zn屬性最佳光譜特征波段Table 6 Best spectral feature wavelengths for soil Zn properties

實(shí)驗(yàn)沒有考慮其他重金屬元素會干擾光譜反演分析過程，本研究認(rèn)為每種重金屬都有很多敏感波段且不盡相同，但在通過優(yōu)選特征波段的組合后可逐步降低與其他重金屬敏感波段的重疊程度。通過第一次篩選和第二次篩選后，將從300～2 500 nm整譜波段中篩選出的5～10個(gè)最佳特征波段進(jìn)行組合，使光譜數(shù)據(jù)與Zn含量的相關(guān)性進(jìn)一步上升，很大程度上降低了與其他重金屬的關(guān)聯(lián)性。

2.5 土壤重金屬Zn元素反演模型

Zn分別以FD、LOG、SG、SD、CR光譜變換為自變量的模型[34]整體上具有較高的Rˉ2，其值均大于0.60，最高為0.70，RMSEC都在20范圍上下（表7）。Zn以SG光譜變換為自變量的模型交叉驗(yàn)證系數(shù)Rˉ2最小，為0.61，以CR光譜變換的交叉驗(yàn)證系數(shù)Rˉ2最高，為0.68。

表7 土壤Zn含量的偏最小二乘法模型的建模與驗(yàn)證Table 7 Modelling and validation of partial least squares models for soil Zn content

對比建模集結(jié)果和交叉驗(yàn)證結(jié)果，Zn含量以CR光譜變換為自變量的模型建模效果較好。檢驗(yàn)的均在驗(yàn)證集檢驗(yàn)的結(jié)果中有不同的變化，但建模效果比較好的模型在驗(yàn)證集的驗(yàn)證中精度仍然較高。進(jìn)行RPD數(shù)據(jù)分析發(fā)現(xiàn)，5種光譜變換模型中CR光譜變換模型的RPD值最高為2.29，其余4種回歸模型的RPD值介于1.71～2.02之間，表明以SD、LOG和SG光譜變換構(gòu)建的偏最小二乘法回歸模型具有定量預(yù)測能力，而CR光譜變換模型的RPD值最高最大，RMSEP最小，表明CR具有很好的模型預(yù)測能力。在Zn含量實(shí)測值與預(yù)測值的擬合散點(diǎn)圖（圖5）中也可以看出，擬合點(diǎn)更多位于y=x直線下方，說明模型的預(yù)測值更多的是小于實(shí)際值；而基于CR變換構(gòu)建的模型中預(yù)測值與實(shí)測值的擬合直線（圖5a）與y=x直線夾角最小，說明其是5種模型中的最優(yōu)模型。

圖5 土壤Zn含量實(shí)測值與預(yù)測值的擬合散點(diǎn)圖Figure 5 Fitting scatter plots of measured and predicted soil Zn content

3 結(jié)論

（1）研究區(qū)土壤Zn含量呈一定程度的偏態(tài)分布，大部分土壤樣品Zn含量介于60～160 mg·kg-1，均值為123.86 mg·kg-1，表明研究區(qū)整體Zn污染風(fēng)險(xiǎn)較低；但存在極端值大于300 mg·kg-1的樣品，說明研究區(qū)內(nèi)局部地點(diǎn)Zn污染風(fēng)險(xiǎn)高。

（2）通過光譜反射率平滑處理消除儀器不同探測元件對光譜數(shù)據(jù)的影響，使反射率曲線變平滑，提高信噪比。通過光譜變換即倒數(shù)對數(shù)（LOG）、一階微分（FD）、二階微分（SD）和去包絡(luò)線（CR）等方法提高提取的特征波段的準(zhǔn)確性和模型的預(yù)測能力，有效突出光譜曲線的吸收和反射特征，從而提高了光譜的靈敏性并增強(qiáng)了有用信息。

（3）以5種反射率變換光譜數(shù)據(jù)為基礎(chǔ)，提取出最佳特征波段作為模型自變量，對應(yīng)土壤Zn含量作為模型因變量，構(gòu)建了基于偏最小二乘法的西平縣高標(biāo)準(zhǔn)基本農(nóng)田建設(shè)區(qū)域土壤Zn屬性反演模型。土壤Zn的最佳模型是以CR光譜變換為最佳的偏最小二乘模型。

本研究僅比較了通過不同光譜變換而挑選出的特征波段經(jīng)同一建模方法建立反演模型的精度，沒有進(jìn)行多種光譜變換方式、多種不同建模方式間的組合比較，土壤光譜數(shù)據(jù)還有待深入挖掘。此外本模型精度有待進(jìn)一步提高，且當(dāng)前模型僅適用于研究區(qū)，而對其他地區(qū)的適應(yīng)性有待研究驗(yàn)證。