富春江典型農(nóng)業(yè)小流域土壤Cd的光譜特征分類及含量估算*

孫顯根 麻萬諸 林陽春 朱康瑩 卓志清# 施加春 陳千千

(1.杭州市生態(tài)環(huán)境局富陽分局，浙江 杭州 311400；2.浙江省農(nóng)業(yè)科學(xué)院數(shù)字農(nóng)業(yè)研究所，浙江 杭州 310021；3.浙江省環(huán)境監(jiān)測工程有限公司，浙江 杭州 310012；4.浙江大學(xué)環(huán)境與資源學(xué)院,土水資源與環(huán)境研究所，浙江 杭州 310058；5.浙江大學(xué)環(huán)境與資源學(xué)院,農(nóng)業(yè)遙感與信息技術(shù)應(yīng)用研究所，浙江 杭州 310058)

土壤重金屬污染因具有持久性、毒性和對區(qū)域環(huán)境影響的難可逆性，導(dǎo)致土壤質(zhì)量退化和生態(tài)環(huán)境惡化，進(jìn)而直接或間接威脅糧食安全和人體健康。重金屬Cd是一種對人體有害且非必需元素，在土壤中具有較強遷移性和較高生物累積性，易通過食物鏈從環(huán)境中富集傳遞到人體，從而引起慢性中毒。據(jù)研究表明，重金屬是影響我國農(nóng)田土壤環(huán)境質(zhì)量的主要無機污染物，其中Cd是首要污染物，在糧食主產(chǎn)區(qū)土壤Cd的點位超標(biāo)率最高達(dá)到17.39%[1]。長期以來，提升土壤重金屬含量監(jiān)測的時效性和準(zhǔn)確性是農(nóng)田土壤重金屬污染防治的研究重點。傳統(tǒng)土壤重金屬含量獲取需進(jìn)行大批量野外取樣和長周期化驗分析，難以滿足農(nóng)田土壤重金屬污染動態(tài)監(jiān)測和精準(zhǔn)防控的需求[2-3]，而現(xiàn)階段光譜技術(shù)反演以其快速、精準(zhǔn)、無損等優(yōu)勢，在土壤重金屬污染狀況監(jiān)測中得到廣泛應(yīng)用[4]。

目前，在土壤Cd含量反演方面，光譜技術(shù)使用的主要波段為可見光-近紅外(VNIR)和中紅外(MIR)兩個波段，VNIR尤為常用。KEMPER等[5]在2002年利用土壤反射光譜估測了西班牙Aznalcollar礦區(qū)重金屬As、Cd、Cu及Zn的含量，認(rèn)為500、1 400、2 200 nm是重金屬元素含量預(yù)測的重要波段，但其對土壤Cd含量的反演精度并不高。為了提升模型的精度和穩(wěn)定性，眾多研究針對光譜數(shù)據(jù)預(yù)處理、特征波段篩選及建模方法等進(jìn)行了大量探索。夏芳等[6]針對農(nóng)田土壤Cd的反演研究發(fā)現(xiàn)，采用一階微分處理后土壤Cd含量的敏感波段與有機質(zhì)的敏感波段有較高重疊，主要集中在400～880、1 350～1 500、2 130～2 350 nm幾個波段；TU等[7]使用競爭自適應(yīng)重加權(quán)采樣算法篩選敏感波段，并采用偏最小二乘回歸(PLSR)預(yù)測了礦區(qū)土壤Cd含量，預(yù)測值與實測值間的相關(guān)系數(shù)達(dá)0.71。此外，也有學(xué)者基于VNIR、MIR及X射線熒光光譜(XRF)等光譜融合技術(shù)對土壤Cd含量進(jìn)行反演，發(fā)現(xiàn)通過光譜數(shù)據(jù)融合能夠提高土壤Cd含量的反演精度，但模型的遷移性仍有不足[8]。

目前針對農(nóng)田土壤Cd含量光譜反演研究的對象多為礦區(qū)及周邊農(nóng)田、金屬回收冶煉區(qū)及周邊農(nóng)田和污灌區(qū)及周邊農(nóng)田等受重金屬污染風(fēng)險較高的區(qū)域，主要存在兩個特點：一是區(qū)域重金屬環(huán)境背景值較高或受污染程度較嚴(yán)重；二是研究中用于構(gòu)建模型的樣本集重金屬含量整體梯度變化較小[9]。基于上述特點所構(gòu)建的土壤Cd含量預(yù)測模型雖然建模效果較好，但多適用于重金屬污染嚴(yán)重的特定區(qū)域，而對于一般污染區(qū)或無污染區(qū)，模型的普適性不足。農(nóng)業(yè)小流域是農(nóng)業(yè)面源污染的基本輸出單元，流域內(nèi)重金屬元素自然遷移及人為活動造成的局地富集以及自然與人為因素的疊加作用，導(dǎo)致土壤Cd含量不僅具有時空分異性，還具有從上游到下游的梯度變化性[10]。因此，針對農(nóng)業(yè)小流域土壤Cd含量開展光譜反演研究，在區(qū)域農(nóng)田土壤重金屬動態(tài)監(jiān)測和污染精準(zhǔn)防控中具有典型代表意義。

剡溪小流域作為富春江的主要集水區(qū)之一，同時也是當(dāng)?shù)刂饕募Z食、蔬菜生產(chǎn)區(qū)，由于長期受自然和人為因素的共同影響，流域內(nèi)的景觀格局發(fā)生了較大變化，從而直接或間接地改變了農(nóng)田土壤重金屬元素的含量及空間分布。本研究以剡溪小流域為研究對象，分析流域內(nèi)農(nóng)田土壤光譜反射特征，探索利用土壤反射光譜強度對土壤Cd含量進(jìn)行聚類的方法，優(yōu)化土壤Cd含量的光譜預(yù)測模型，推動光譜技術(shù)在農(nóng)田土壤重金屬含量預(yù)測中的應(yīng)用。

1 數(shù)據(jù)采集與測試

1.1 研究區(qū)概況

剡溪是錢塘江流域富春江段的主要溪流之一，發(fā)源于杭州市富陽區(qū)南部龍門山，經(jīng)上官、龍門在環(huán)山鄉(xiāng)注入富春江。流域匯水面積87.10 km2。該區(qū)域?qū)俦眮啛釒Ъ撅L(fēng)氣候區(qū)，年平均氣溫16.27 ℃，年平均降水量1 235.30 mm，主要集中在每年5—6月。流域內(nèi)總體地勢自東南向西北傾斜，最高海拔1 057 m；東南部是以流紋角質(zhì)凝灰?guī)r和砂質(zhì)泥巖為主要構(gòu)成的山地和丘陵，西北部主要為沖洪積物堆積的平原和河漫灘。山地和丘陵區(qū)多為砂質(zhì)紅壤或礫石紅壤，平原區(qū)土壤類型以淹育水稻土和滲育水稻土為主。流域內(nèi)主要農(nóng)用地類型為水田、旱地等，農(nóng)田集中分布在下游平原地區(qū)，主要種植水稻、蔬菜等，是當(dāng)?shù)刂匾募Z食和蔬菜生產(chǎn)功能區(qū)[11]。

1.2 樣品采集與分析

依據(jù)流域內(nèi)水系和農(nóng)田分布特征，采樣點集中布設(shè)在中、下游平原地區(qū)，并在下游人口密集區(qū)和工業(yè)集中分布區(qū)加密布點；根據(jù)流域內(nèi)的土壤類型及分布，結(jié)合網(wǎng)格布點法確定采樣點數(shù)量，確保每個土壤類型均有采樣點，綜合考慮上述原則共確定126個采樣點(見圖1)。土壤樣品采集時間為2021年9—10月，各采樣點按梅花法采集5個重復(fù)土樣，取表層(0～20 cm)土壤約1 kg，經(jīng)混勻后用四分法再取1 kg作為混合土樣；記錄采樣點坐標(biāo)，在樣品采集和運轉(zhuǎn)過程中避免與金屬器皿直接接觸。采集的樣品經(jīng)室內(nèi)風(fēng)干、剔除植物殘體和石塊、研磨過100目尼龍網(wǎng)篩后分兩份保存?zhèn)溆茫环萦糜谕寥繡d含量測定，一份用于土壤光譜測試。土壤Cd總量用王水提取，電感耦合等離子體質(zhì)譜法(ICP-MS)測定。

1.3 土壤光譜測定

使用FieldSpec 4便攜式地物光譜儀(美國ASD公司)測定研究區(qū)土壤樣品的VNIR光譜(350～2 500 nm)，采樣間隔為1 nm，共輸出2 150個波段。將土壤樣品平鋪于高1 cm的玻璃盛樣皿內(nèi)，用地物光譜儀中內(nèi)置光源的高強度接觸式探頭測量，每次測量前用白板校正。為提高土壤光譜數(shù)據(jù)的測量精度，每個樣品測10次重復(fù)，以平均值作為該土樣的光譜反射率。

2 研究方法

2.1 光譜數(shù)據(jù)預(yù)處理

在VNIR波段中剔除2個噪聲較大的波段(350～399、2 401～2 500 nm)，為減少實驗室光學(xué)環(huán)境場差異和磨樣過篩的影響，采用Savitzky-Golay(SG)算法對400～2 400 nm的原始土壤光譜(R)曲線進(jìn)行平滑處理，并采用4種變換形式對平滑后的光譜曲線進(jìn)行預(yù)處理，包括一階微分(R′)、二階微分(R′′)、倒數(shù)對數(shù)(lg(1/R))和連續(xù)統(tǒng)去除(CR)。為壓縮土壤光譜數(shù)據(jù)，采用主成分分析(PCA)方法對原始光譜數(shù)據(jù)降維，并對轉(zhuǎn)換后的PCA結(jié)果做進(jìn)一步分析。

2.2 光譜數(shù)據(jù)分類方法

本研究采用模糊C均值(FCM)聚類法進(jìn)行光譜數(shù)據(jù)分類，其主要方法是將一個數(shù)據(jù)集劃分為多個類別，尋找目標(biāo)函數(shù)的迭代最小化，給出最佳分類數(shù)目的量化方法[12]。研究引入分類系數(shù)(PC)和分類熵(PE)來確定最佳光譜分類數(shù)。PC和PE取值均為0～1，PC越接近1，表示聚類時共用數(shù)據(jù)越少，分類的劃分越明顯。PE越接近0，分區(qū)內(nèi)的數(shù)據(jù)相似程度越高，聚類效果越好[13]。本研究通過對不同分類數(shù)目進(jìn)行聚類分析，以PC和PE同時達(dá)到最大值和最小值時的聚類數(shù)為最佳分類數(shù)。

2.3 模型建立與驗證

采用PLSR構(gòu)建土壤Cd含量反演模型，該方法可簡化光譜數(shù)據(jù)結(jié)構(gòu)，有利于解決自變量之間多重相關(guān)的現(xiàn)象，能避免建模中的過擬合問題。預(yù)測模型采用Leave-one-out方法進(jìn)行檢驗，并選取決定系數(shù)(r2)和均方根誤差(RMSE)估算模型的反演精度。r2越大，RMSE越小，說明預(yù)測效果越好。

研究中針對光譜數(shù)據(jù)的微分轉(zhuǎn)換、主成分分析、FCM分類、PLSR建模和驗證等處理分別在Unscrambler X10.1 軟件和R軟件(v4.2.0)中完成。

3 結(jié)果與分析

3.1 土壤Cd含量與光譜特征

研究區(qū)126個土壤樣品的Cd質(zhì)量濃度為0.13～2.92 mg/kg，中位值為1.37 mg/kg，偏度和峰度分別為0.29、2.58，可知Cd含量存在偏高值，但樣本數(shù)據(jù)總體分布均衡，具有代表性；Cd含量變異系數(shù)為38.50%，屬于中等程度的空間變異。基于所有樣品的光譜反射率繪制光譜曲線，為量化土壤光譜反射率對土壤Cd含量變化的響應(yīng)程度，依據(jù)土壤Cd含量將其分為三級，分別繪制一級(25%分位數(shù))、二級(平均值)和三級(75%分位數(shù))含量的光譜反射率曲線(見圖2)。研究表明，隨著土壤Cd含量的增加，光譜反射率逐漸減小，說明光譜反射率與土壤Cd含量具有一定的相關(guān)性。不同土壤Cd含量下土壤光譜反射率曲線整體變化趨勢相近，在400～900 nm曲線斜率較大，隨著波長的增加，反射率呈現(xiàn)快速上升趨勢；在1 000 nm以后，曲線整體呈現(xiàn)平緩上升的趨勢。在波長1 410～1 415、1 910～1 920、2 205～2 210 nm附近土壤光譜反射率曲線存在3個明顯的水分吸收谷；在2 100～2 140 nm 波段出現(xiàn)反射峰。

3.2 土壤光譜特征的模糊聚類分析

由于土壤樣本光譜測試數(shù)據(jù)量較大，且各波段之間存在多重相關(guān)性，因此先對全部樣品的光譜反射率數(shù)據(jù)進(jìn)行數(shù)據(jù)轉(zhuǎn)換和平滑處理，通過對后期預(yù)測效果的比較，選用lg(1/R)進(jìn)行數(shù)據(jù)處理，并采用PCA方法對預(yù)處理后的光譜數(shù)據(jù)進(jìn)行降維，得到2個主成分，第一主成分(PCA1)貢獻(xiàn)率86.76%，第二主成分(PCA2)貢獻(xiàn)率7.09%，兩者合計貢獻(xiàn)率達(dá)93.85%。運用FCM聚類法對土壤光譜數(shù)據(jù)進(jìn)行分類，并確定最佳分類數(shù)目。對前兩個主成分?jǐn)?shù)據(jù)進(jìn)行歸一化處理，之后作為輸入變量用于FCM聚類分析。為了找出最佳分類數(shù)目，分別產(chǎn)生2、3、4、5、6、7、8、9、10個類別，計算并比較不同分類數(shù)目下的PC和PE。結(jié)果表明，當(dāng)分類數(shù)目為2時，PC和PE同時達(dá)到最大值和最小值(見圖3)，因此研究區(qū)土壤反射光譜數(shù)據(jù)的最佳分類數(shù)目為2。

根據(jù)確定的土壤光譜數(shù)據(jù)最佳分類數(shù)目，將研究區(qū)內(nèi)的土壤樣品對應(yīng)分成兩類，各類別的聚類中心和主成分值如圖4所示。結(jié)果表明，PCA1是土壤光譜分類主要依據(jù)，即可以表征研究區(qū)土壤反射光譜強度；PCA2是表征反射光譜曲線形狀特征的主要參數(shù)。類型1聚類中心的PCA2為-0.015，類型2聚類中心的PCA2為0.009，均在0附近，可見二者譜線的形狀總體差異不大。

不同類型土壤的光譜反射率曲線見圖5。可以看出類型1的光譜反射率整體大于類型2，類型1所屬樣本的土壤Cd均值(0.33 mg/kg)低于類型2(1.96 mg/kg)。從研究區(qū)土壤樣本類型的空間分布來看，類型1的樣本主要為小流域上游采樣點，所處地形以丘陵、山地為主，農(nóng)田分布面積較小且零散，區(qū)域內(nèi)土壤環(huán)境受人為活動干擾較少，土壤Cd含量較低。類型2的樣本主要分布在小流域下游沖積平原，地形平坦，耕地集中分布；另外，該區(qū)域歷史上分布較多金屬加工冶煉企業(yè)，周邊農(nóng)田土壤受含重金屬廢水、大氣沉降的影響較大，土壤Cd含量相對較高(見圖6)。

3.3 土壤Cd含量的光譜反演建模

以采樣點土壤光譜全波段反射率作為自變量，土壤Cd含量作為因變量，基于光譜反射率原始數(shù)據(jù)和其他4種預(yù)處理方法數(shù)據(jù)，采用PLSR構(gòu)建土壤Cd含量反演模型，并進(jìn)行內(nèi)部交叉驗證，結(jié)果見表1。結(jié)果表明，基于lg(1/R)處理后建立的PLSR模型預(yù)測精度最高，訓(xùn)練集和驗證集的r2均達(dá)到0.60以上，RMSE在1.25 mg/kg以下；其次為一階微分處理，訓(xùn)練集和驗證集的r2分別為0.60和0.57，RMSE分別為1.18、1.31 mg/kg。與原始數(shù)據(jù)相比，經(jīng)lg(1/R)處理后建模集、驗證集的r2分別提升了0.16、0.17，RMSE分別降低了0.44、0.45 mg/kg，說明該模型具有較好的預(yù)測能力。

為了進(jìn)一步分析模型的預(yù)測效果，對比基于R、lg(1/R)構(gòu)建的土壤Cd含量PLSR反演模型，分別以土壤Cd含量的預(yù)測值與實測值繪制散點圖。由圖7可知，光譜原始數(shù)據(jù)的預(yù)測值和實測值基本分布在1∶1線附近，但基于lg(1/R)構(gòu)建的反演模型，與1∶1線更為貼近，該結(jié)果也表明，經(jīng)lg(1/R)處理后的PLSR反演模型具有更高的估測精度和穩(wěn)定性，可作為研究區(qū)土壤Cd含量的優(yōu)選估測方法。

4 討 論

本研究基于小流域土壤重金屬元素的遷移和富集特征，采用FCM聚類方法對采樣點土壤反射光譜數(shù)據(jù)進(jìn)行分類，從不同類型反射光譜的視角分析土壤Cd含量的差異，從而實現(xiàn)對土壤Cd含量直接反演。研究結(jié)果表明，土壤光譜反射率隨Cd含量的增加而減小，光譜反射率與Cd含量具有一定的相關(guān)性。已有研究發(fā)現(xiàn)土壤Cd的光譜特征微弱，敏感波段不易確定，進(jìn)而導(dǎo)致土壤Cd含量與原始光譜反射率相關(guān)性不明顯[14]。諸多研究利用多種光譜預(yù)處理手段結(jié)合特征波段提取從而獲得土壤Cd含量的敏感波段，但受區(qū)域土壤類型、成土因素、重金屬含量及光譜數(shù)據(jù)處理方法等影響，所提取的敏感波段存在較大差異[15-16]。本研究中雖然采用PCA對多種預(yù)處理方法下的光譜數(shù)據(jù)特征變量進(jìn)行了篩選，但篩選后的主成分變量僅作為模糊聚類中的輸入變量；而在PLSR構(gòu)建的反演模型中，以土壤的全波段光譜反射率作為輸入變量。因此，在后續(xù)研究中可以嘗試其他數(shù)據(jù)降維方法或敏感波段篩選方法。

表1 土壤Cd的PLSR反演模型Table 1 PLSR prediction model of soil Cd

本研究通過土壤光譜數(shù)據(jù)直接估測了土壤Cd含量，基于lg(1/R)處理后建立的PLSR模型預(yù)測精度最高，模型訓(xùn)練集和驗證集的r2均達(dá)到0.60以上，模型具有一定適用性，但與前人研究相比r2相對較低。這主要因為前人研究多關(guān)注土壤重金屬污染嚴(yán)重區(qū)域，樣本的Cd含量差異較小，直接反演易于實現(xiàn)[17]。由于受樣本數(shù)量和研究區(qū)面積的限制，本研究只采用了普適性較強的PLSR構(gòu)建Cd含量的光譜反演模型，未來需針對大樣本量的土壤Cd含量光譜進(jìn)行反演建模，可探索隨機森林、神經(jīng)網(wǎng)絡(luò)等機器學(xué)習(xí)方法以提高模型的預(yù)測精度。另外，土壤組分中有機質(zhì)、黏土礦物及氧化鐵等對土壤反射光譜形態(tài)具有明顯影響[18]，這些土壤組分因參與土壤中重金屬的吸附、絡(luò)合及氧化還原等一系列反應(yīng)，均對土壤Cd含量的反演精度有一定影響[19-20]。有研究表明，基于土壤重金屬含量與土壤理化屬性之間的相關(guān)性，可間接反演土壤Cd含量，從而提升反演精度[21]，因此在后續(xù)研究中可針對流域土壤Cd含量間接反演展開進(jìn)一步研究。

5 結(jié) 論

(1) 富春江剡溪小流域土壤光譜反射率與土壤Cd含量具有一定的相關(guān)性，隨土壤Cd含量的增加而減小。不同土壤Cd含量條件下光譜反射率曲線的整體變化趨勢相近。通過對原始光譜數(shù)據(jù)進(jìn)行多種光譜變換，可以有效地去除噪聲，提高相關(guān)性。

(2) 采用FCM方法可有效量化研究區(qū)土壤反射光譜的最佳分類數(shù)目，實現(xiàn)對小流域內(nèi)土壤Cd含量的分類。其中類型1的光譜反射率整體大于類型2，類型1所屬樣本的Cd含量低于類型2，兩種光譜類型所屬采樣點呈現(xiàn)明顯空間集聚特征。

(3) 與其他幾種預(yù)處理方法相比，基于lg(1/R)處理后建立的PLSR預(yù)測模型精度最高，訓(xùn)練集和驗證集的r2均達(dá)到0.60以上；RMSE在1.25 mg/kg以下，模型具有較好的預(yù)測能力和穩(wěn)定性，可優(yōu)先作為典型農(nóng)業(yè)小流域土壤Cd含量的估測方法。