西平縣高標準農田Zn含量高光譜反演研究

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

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

駐馬店市西平縣位于北緯33°10'～33°32'，東經(jīng)113°36'～114°13'（圖1），河南中南部區(qū)域，其土壤類型分布及面積占比見表1。縣境總面積為1 089.77 km2，劃分為3個街道辦事處、8鎮(zhèn)、8鄉(xiāng)和老王坡管理委員會等行政區(qū)。該地氣候屬亞濕潤大陸性季風型，全年氣候溫和。西平縣資源豐富，農作物物種豐盛，是我國糧食生產先進縣與糧食主產區(qū)之一，同時其西部山區(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ū)的土壤類型、空間特征、行政單元和地形因素進行分析，以設定精確定位的采樣點。利用對角線采樣法采樣，野外采樣深度為0～30 cm土壤層，同時記錄詳細地形特征。從田里采集的土樣及時攤開風干，倒在干凈的容器內，去除石塊、植物根莖殘體及其余雜物。取研磨均勻后通過1 mm孔徑孔篩的土壤樣本，利用四分法取樣，用于土壤光譜數(shù)據(jù)的采集和Zn含量的測定[16]。

1.3 土樣光譜室內測定

1.4 光譜變換方法

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

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

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

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

1.5 相關分析與顯著性檢驗

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

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

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

1.6 連續(xù)投影算法

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

1.7 偏最小二乘法回歸模型反演與精度檢驗方法

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

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

表3 精度校驗方法公式Table 3 Precision calibration method formula

2 結果與分析

2.1 土壤Zn屬性特征

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

圖2 Zn質量分數(shù)描述統(tǒng)計圖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范圍處的反射率較低，但呈大幅度上升的趨勢，光譜曲線陡峭，土壤有機質在600、800 nm附近有微弱的反射峰，在780～2 300 nm處近紅外波段區(qū)域曲線較平緩，在2 300～2 400 nm處逐漸下降。1 850～2 140 nm光譜曲線呈現(xiàn)明顯的“V”字形。圖中3個水分吸收帶表現(xiàn)出的波谷與黏土中的多種礦物質組成有關，1 420 nm為羥基（OH）譜帶、1 920 nm附近為H2O譜帶、2 200 nm波段范圍為羥基伸縮振動與和OH彎曲振動的合譜帶[32]。4條平均光譜曲線整體趨勢相似，丘崗地區(qū)的平均光譜曲線較低，山地丘陵區(qū)的平均光譜曲線最高，而縣區(qū)與平原區(qū)的平均光譜曲線居中，幾乎吻合。

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

2.3 相關分析與顯著性檢驗

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

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

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

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

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

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

2.4 最佳特征波段的選取

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

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

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

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

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

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

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

圖5 土壤Zn含量實測值與預測值的擬合散點圖Figure 5 Fitting scatter plots of measured and predicted soil Zn content

3 結論

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

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

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

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