富春江典型农业小流域土壤Cd的光谱特征分类及含量估算*

孙显根 麻万诸 林阳春 朱康莹 卓志清# 施加春 陈千千

(1.杭州市生态环境局富阳分局，浙江 杭州 311400；2.浙江省农业科学院数字农业研究所，浙江 杭州 310021；3.浙江省环境监测工程有限公司，浙江 杭州 310012；4.浙江大学环境与资源学院,土水资源与环境研究所，浙江 杭州 310058；5.浙江大学环境与资源学院,农业遥感与信息技术应用研究所，浙江 杭州 310058)

土壤重金属污染因具有持久性、毒性和对区域环境影响的难可逆性，导致土壤质量退化和生态环境恶化，进而直接或间接威胁粮食安全和人体健康。重金属Cd是一种对人体有害且非必需元素，在土壤中具有较强迁移性和较高生物累积性，易通过食物链从环境中富集传递到人体，从而引起慢性中毒。据研究表明，重金属是影响我国农田土壤环境质量的主要无机污染物，其中Cd是首要污染物，在粮食主产区土壤Cd的点位超标率最高达到17.39%[1]。长期以来，提升土壤重金属含量监测的时效性和准确性是农田土壤重金属污染防治的研究重点。传统土壤重金属含量获取需进行大批量野外取样和长周期化验分析，难以满足农田土壤重金属污染动态监测和精准防控的需求[2-3]，而现阶段光谱技术反演以其快速、精准、无损等优势，在土壤重金属污染状况监测中得到广泛应用[4]。

目前，在土壤Cd含量反演方面，光谱技术使用的主要波段为可见光-近红外(VNIR)和中红外(MIR)两个波段，VNIR尤为常用。KEMPER等[5]在2002年利用土壤反射光谱估测了西班牙Aznalcollar矿区重金属As、Cd、Cu及Zn的含量，认为500、1 400、2 200 nm是重金属元素含量预测的重要波段，但其对土壤Cd含量的反演精度并不高。为了提升模型的精度和稳定性，众多研究针对光谱数据预处理、特征波段筛选及建模方法等进行了大量探索。夏芳等[6]针对农田土壤Cd的反演研究发现，采用一阶微分处理后土壤Cd含量的敏感波段与有机质的敏感波段有较高重叠，主要集中在400～880、1 350～1 500、2 130～2 350 nm几个波段；TU等[7]使用竞争自适应重加权采样算法筛选敏感波段，并采用偏最小二乘回归(PLSR)预测了矿区土壤Cd含量，预测值与实测值间的相关系数达0.71。此外，也有学者基于VNIR、MIR及X射线荧光光谱(XRF)等光谱融合技术对土壤Cd含量进行反演，发现通过光谱数据融合能够提高土壤Cd含量的反演精度，但模型的迁移性仍有不足[8]。

目前针对农田土壤Cd含量光谱反演研究的对象多为矿区及周边农田、金属回收冶炼区及周边农田和污灌区及周边农田等受重金属污染风险较高的区域，主要存在两个特点：一是区域重金属环境背景值较高或受污染程度较严重；二是研究中用于构建模型的样本集重金属含量整体梯度变化较小[9]。基于上述特点所构建的土壤Cd含量预测模型虽然建模效果较好，但多适用于重金属污染严重的特定区域，而对于一般污染区或无污染区，模型的普适性不足。农业小流域是农业面源污染的基本输出单元，流域内重金属元素自然迁移及人为活动造成的局地富集以及自然与人为因素的叠加作用，导致土壤Cd含量不仅具有时空分异性，还具有从上游到下游的梯度变化性[10]。因此，针对农业小流域土壤Cd含量开展光谱反演研究，在区域农田土壤重金属动态监测和污染精准防控中具有典型代表意义。

剡溪小流域作为富春江的主要集水区之一，同时也是当地主要的粮食、蔬菜生产区，由于长期受自然和人为因素的共同影响，流域内的景观格局发生了较大变化，从而直接或间接地改变了农田土壤重金属元素的含量及空间分布。本研究以剡溪小流域为研究对象，分析流域内农田土壤光谱反射特征，探索利用土壤反射光谱强度对土壤Cd含量进行聚类的方法，优化土壤Cd含量的光谱预测模型，推动光谱技术在农田土壤重金属含量预测中的应用。

1 数据采集与测试

1.1 研究区概况

剡溪是钱塘江流域富春江段的主要溪流之一，发源于杭州市富阳区南部龙门山，经上官、龙门在环山乡注入富春江。流域汇水面积87.10 km2。该区域属北亚热带季风气候区，年平均气温16.27 ℃，年平均降水量1 235.30 mm，主要集中在每年5—6月。流域内总体地势自东南向西北倾斜，最高海拔1 057 m；东南部是以流纹角质凝灰岩和砂质泥岩为主要构成的山地和丘陵，西北部主要为冲洪积物堆积的平原和河漫滩。山地和丘陵区多为砂质红壤或砾石红壤，平原区土壤类型以淹育水稻土和渗育水稻土为主。流域内主要农用地类型为水田、旱地等，农田集中分布在下游平原地区，主要种植水稻、蔬菜等，是当地重要的粮食和蔬菜生产功能区[11]。

1.2 样品采集与分析

依据流域内水系和农田分布特征，采样点集中布设在中、下游平原地区，并在下游人口密集区和工业集中分布区加密布点；根据流域内的土壤类型及分布，结合网格布点法确定采样点数量，确保每个土壤类型均有采样点，综合考虑上述原则共确定126个采样点(见图1)。土壤样品采集时间为2021年9—10月，各采样点按梅花法采集5个重复土样，取表层(0～20 cm)土壤约1 kg，经混匀后用四分法再取1 kg作为混合土样；记录采样点坐标，在样品采集和运转过程中避免与金属器皿直接接触。采集的样品经室内风干、剔除植物残体和石块、研磨过100目尼龙网筛后分两份保存备用，一份用于土壤Cd含量测定，一份用于土壤光谱测试。土壤Cd总量用王水提取，电感耦合等离子体质谱法(ICP-MS)测定。

1.3 土壤光谱测定

使用FieldSpec 4便携式地物光谱仪(美国ASD公司)测定研究区土壤样品的VNIR光谱(350～2 500 nm)，采样间隔为1 nm，共输出2 150个波段。将土壤样品平铺于高1 cm的玻璃盛样皿内，用地物光谱仪中内置光源的高强度接触式探头测量，每次测量前用白板校正。为提高土壤光谱数据的测量精度，每个样品测10次重复，以平均值作为该土样的光谱反射率。

2 研究方法

2.1 光谱数据预处理

在VNIR波段中剔除2个噪声较大的波段(350～399、2 401～2 500 nm)，为减少实验室光学环境场差异和磨样过筛的影响，采用Savitzky-Golay(SG)算法对400～2 400 nm的原始土壤光谱(R)曲线进行平滑处理，并采用4种变换形式对平滑后的光谱曲线进行预处理，包括一阶微分(R′)、二阶微分(R′′)、倒数对数(lg(1/R))和连续统去除(CR)。为压缩土壤光谱数据，采用主成分分析(PCA)方法对原始光谱数据降维，并对转换后的PCA结果做进一步分析。

2.2 光谱数据分类方法

本研究采用模糊C均值(FCM)聚类法进行光谱数据分类，其主要方法是将一个数据集划分为多个类别，寻找目标函数的迭代最小化，给出最佳分类数目的量化方法[12]。研究引入分类系数(PC)和分类熵(PE)来确定最佳光谱分类数。PC和PE取值均为0～1，PC越接近1，表示聚类时共用数据越少，分类的划分越明显。PE越接近0，分区内的数据相似程度越高，聚类效果越好[13]。本研究通过对不同分类数目进行聚类分析，以PC和PE同时达到最大值和最小值时的聚类数为最佳分类数。

2.3 模型建立与验证

采用PLSR构建土壤Cd含量反演模型，该方法可简化光谱数据结构，有利于解决自变量之间多重相关的现象，能避免建模中的过拟合问题。预测模型采用Leave-one-out方法进行检验，并选取决定系数(r2)和均方根误差(RMSE)估算模型的反演精度。r2越大，RMSE越小，说明预测效果越好。

研究中针对光谱数据的微分转换、主成分分析、FCM分类、PLSR建模和验证等处理分别在Unscrambler X10.1 软件和R软件(v4.2.0)中完成。

3 结果与分析

3.1 土壤Cd含量与光谱特征

研究区126个土壤样品的Cd质量浓度为0.13～2.92 mg/kg，中位值为1.37 mg/kg，偏度和峰度分别为0.29、2.58，可知Cd含量存在偏高值，但样本数据总体分布均衡，具有代表性；Cd含量变异系数为38.50%，属于中等程度的空间变异。基于所有样品的光谱反射率绘制光谱曲线，为量化土壤光谱反射率对土壤Cd含量变化的响应程度，依据土壤Cd含量将其分为三级，分别绘制一级(25%分位数)、二级(平均值)和三级(75%分位数)含量的光谱反射率曲线(见图2)。研究表明，随着土壤Cd含量的增加，光谱反射率逐渐减小，说明光谱反射率与土壤Cd含量具有一定的相关性。不同土壤Cd含量下土壤光谱反射率曲线整体变化趋势相近，在400～900 nm曲线斜率较大，随着波长的增加，反射率呈现快速上升趋势；在1 000 nm以后，曲线整体呈现平缓上升的趋势。在波长1 410～1 415、1 910～1 920、2 205～2 210 nm附近土壤光谱反射率曲线存在3个明显的水分吸收谷；在2 100～2 140 nm 波段出现反射峰。

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)，因此研究区土壤反射光谱数据的最佳分类数目为2。

根据确定的土壤光谱数据最佳分类数目，将研究区内的土壤样品对应分成两类，各类别的聚类中心和主成分值如图4所示。结果表明，PCA1是土壤光谱分类主要依据，即可以表征研究区土壤反射光谱强度；PCA2是表征反射光谱曲线形状特征的主要参数。类型1聚类中心的PCA2为-0.015，类型2聚类中心的PCA2为0.009，均在0附近，可见二者谱线的形状总体差异不大。

不同类型土壤的光谱反射率曲线见图5。可以看出类型1的光谱反射率整体大于类型2，类型1所属样本的土壤Cd均值(0.33 mg/kg)低于类型2(1.96 mg/kg)。从研究区土壤样本类型的空间分布来看，类型1的样本主要为小流域上游采样点，所处地形以丘陵、山地为主，农田分布面积较小且零散，区域内土壤环境受人为活动干扰较少，土壤Cd含量较低。类型2的样本主要分布在小流域下游冲积平原，地形平坦，耕地集中分布；另外，该区域历史上分布较多金属加工冶炼企业，周边农田土壤受含重金属废水、大气沉降的影响较大，土壤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反演模型具有更高的估测精度和稳定性，可作为研究区土壤Cd含量的优选估测方法。

4 讨 论

本研究基于小流域土壤重金属元素的迁移和富集特征，采用FCM聚类方法对采样点土壤反射光谱数据进行分类，从不同类型反射光谱的视角分析土壤Cd含量的差异，从而实现对土壤Cd含量直接反演。研究结果表明，土壤光谱反射率随Cd含量的增加而减小，光谱反射率与Cd含量具有一定的相关性。已有研究发现土壤Cd的光谱特征微弱，敏感波段不易确定，进而导致土壤Cd含量与原始光谱反射率相关性不明显[14]。诸多研究利用多种光谱预处理手段结合特征波段提取从而获得土壤Cd含量的敏感波段，但受区域土壤类型、成土因素、重金属含量及光谱数据处理方法等影响，所提取的敏感波段存在较大差异[15-16]。本研究中虽然采用PCA对多种预处理方法下的光谱数据特征变量进行了筛选，但筛选后的主成分变量仅作为模糊聚类中的输入变量；而在PLSR构建的反演模型中，以土壤的全波段光谱反射率作为输入变量。因此，在后续研究中可以尝试其他数据降维方法或敏感波段筛选方法。

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

本研究通过土壤光谱数据直接估测了土壤Cd含量，基于lg(1/R)处理后建立的PLSR模型预测精度最高，模型训练集和验证集的r2均达到0.60以上，模型具有一定适用性，但与前人研究相比r2相对较低。这主要因为前人研究多关注土壤重金属污染严重区域，样本的Cd含量差异较小，直接反演易于实现[17]。由于受样本数量和研究区面积的限制，本研究只采用了普适性较强的PLSR构建Cd含量的光谱反演模型，未来需针对大样本量的土壤Cd含量光谱进行反演建模，可探索随机森林、神经网络等机器学习方法以提高模型的预测精度。另外，土壤组分中有机质、黏土矿物及氧化铁等对土壤反射光谱形态具有明显影响[18]，这些土壤组分因参与土壤中重金属的吸附、络合及氧化还原等一系列反应，均对土壤Cd含量的反演精度有一定影响[19-20]。有研究表明，基于土壤重金属含量与土壤理化属性之间的相关性，可间接反演土壤Cd含量，从而提升反演精度[21]，因此在后续研究中可针对流域土壤Cd含量间接反演展开进一步研究。

5 结 论

(1) 富春江剡溪小流域土壤光谱反射率与土壤Cd含量具有一定的相关性，随土壤Cd含量的增加而减小。不同土壤Cd含量条件下光谱反射率曲线的整体变化趋势相近。通过对原始光谱数据进行多种光谱变换，可以有效地去除噪声，提高相关性。

(2) 采用FCM方法可有效量化研究区土壤反射光谱的最佳分类数目，实现对小流域内土壤Cd含量的分类。其中类型1的光谱反射率整体大于类型2，类型1所属样本的Cd含量低于类型2，两种光谱类型所属采样点呈现明显空间集聚特征。

(3) 与其他几种预处理方法相比，基于lg(1/R)处理后建立的PLSR预测模型精度最高，训练集和验证集的r2均达到0.60以上；RMSE在1.25 mg/kg以下，模型具有较好的预测能力和稳定性，可优先作为典型农业小流域土壤Cd含量的估测方法。