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

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

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

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

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

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

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

1 數據采集與測試

1.1 研究區(qū)概況

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

1.2 樣品采集與分析

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

1.3 土壤光譜測定

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

2 研究方法

2.1 光譜數據預處理

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

2.2 光譜數據分類方法

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

2.3 模型建立與驗證

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

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

3 結果與分析

3.1 土壤Cd含量與光譜特征

研究區(qū)126個土壤樣品的Cd質量濃度為0.13～2.92 mg/kg，中位值為1.37 mg/kg，偏度和峰度分別為0.29、2.58，可知Cd含量存在偏高值，但樣本數據總體分布均衡，具有代表性；Cd含量變異系數為38.50%，屬于中等程度的空間變異?；谒袠悠返墓庾V反射率繪制光譜曲線，為量化土壤光譜反射率對土壤Cd含量變化的響應程度，依據土壤Cd含量將其分為三級，分別繪制一級(25%分位數)、二級(平均值)和三級(75%分位數)含量的光譜反射率曲線(見圖2)。研究表明，隨著土壤Cd含量的增加，光譜反射率逐漸減小，說明光譜反射率與土壤Cd含量具有一定的相關性。不同土壤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 土壤光譜特征的模糊聚類分析

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

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

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

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

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

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

4 討 論

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

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

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

5 結 論

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

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

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