孫顯根 麻萬(wàn)諸 林陽(yáng)春 朱康瑩 卓志清# 施加春 陳千千

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

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

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

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

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

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

1.1 研究區(qū)概況

剡溪是錢塘江流域富春江段的主要溪流之一，發(fā)源于杭州市富陽(yáng)區(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)總體地勢(shì)自東南向西北傾斜，最高海拔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)田分布特征，采樣點(diǎn)集中布設(shè)在中、下游平原地區(qū)，并在下游人口密集區(qū)和工業(yè)集中分布區(qū)加密布點(diǎn)；根據(jù)流域內(nèi)的土壤類型及分布，結(jié)合網(wǎng)格布點(diǎn)法確定采樣點(diǎn)數(shù)量，確保每個(gè)土壤類型均有采樣點(diǎn)，綜合考慮上述原則共確定126個(gè)采樣點(diǎn)(見(jiàn)圖1)。土壤樣品采集時(shí)間為2021年9—10月，各采樣點(diǎn)按梅花法采集5個(gè)重復(fù)土樣，取表層(0～20 cm)土壤約1 kg，經(jīng)混勻后用四分法再取1 kg作為混合土樣；記錄采樣點(diǎn)坐標(biāo)，在樣品采集和運(yùn)轉(zhuǎn)過(guò)程中避免與金屬器皿直接接觸。采集的樣品經(jīng)室內(nèi)風(fēng)干、剔除植物殘?bào)w和石塊、研磨過(guò)100目尼龍網(wǎng)篩后分兩份保存?zhèn)溆?，一份用于土壤Cd含量測(cè)定，一份用于土壤光譜測(cè)試。土壤Cd總量用王水提取，電感耦合等離子體質(zhì)譜法(ICP-MS)測(cè)定。

1.3 土壤光譜測(cè)定

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

2 研究方法

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

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

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

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

2.3 模型建立與驗(yàn)證

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

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

3 結(jié)果與分析

3.1 土壤Cd含量與光譜特征

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

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

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

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

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

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

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

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

4 討 論

本研究基于小流域土壤重金屬元素的遷移和富集特征，采用FCM聚類方法對(duì)采樣點(diǎn)土壤反射光譜數(shù)據(jù)進(jìn)行分類，從不同類型反射光譜的視角分析土壤Cd含量的差異，從而實(shí)現(xiàn)對(duì)土壤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對(duì)多種預(yù)處理方法下的光譜數(shù)據(jù)特征變量進(jìn)行了篩選，但篩選后的主成分變量?jī)H作為模糊聚類中的輸入變量；而在PLSR構(gòu)建的反演模型中，以土壤的全波段光譜反射率作為輸入變量。因此，在后續(xù)研究中可以嘗試其他數(shù)據(jù)降維方法或敏感波段篩選方法。

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

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

5 結(jié) 論

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

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

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