西平县高标准农田Zn含量高光谱反演研究

蔡太义，王志刚，杨刘帅，2，王群，黄会娟，于海洋，张传忠，张灿，刘鹏，冯玉庆，贺成龙，张合兵*

（1.河南理工大学测绘与国土信息工程学院，河南 焦作 454000；2.自然资源部第四地形测量队，哈尔滨 150000；3.河南农业大学农学院/河南省农田生态系统野外科学观测研究站，郑州 450046；4.河南省土壤调理与修复工程技术研究中心，河南 商丘 476000；5.中向旭曜科技有限公司，江苏 无锡 214400；6.河南省地质局矿产资源勘查中心,郑州 450053）

土壤重金属是土壤复杂生态系统中重要的组成成分，影响着土壤环境质量[1]。锌（Zn）是一种重要的土壤重金属元素[2]，是在大气沉降、工业污染和施肥等过程中缓慢积累，且大部分处于结合状态的金属离子，并且其不易降解[3]。我国土壤中Zn的环境背景含量在2.60～593.00 mg·kg-1之间[4]。Zn含量过低会减少作物激素的生成，过高则会对作物产生毒害作用[5]。因此，Zn含量的研究是高标准农田建设过程中的一个重要内容。传统重金属测量与监测方法具有耗时长、范围窄和实时性差的缺点，不能满足反映生态系统结构变化的要求[6]。高光谱遥感数据可以快速检测土壤重金属属性信息[7]，能节省人力与物力，同时具有成本低、无污染和无损的优势[8]，是土壤重金属含量研究的新趋势[9]。

前人对Zn元素含量研究已取得重大进展。陶超等[10]利用定量回归模型和定性反演模型比较采样点重金属铅（Pb）与Zn的可迁移能力，证明在快速检测土壤重金属污染状况的问题上，定性分类是一种更加切实可行的方式。王金凤等[11]运用电感耦合等离子体质谱和分光光度计测定土壤样品，分析土壤光谱和Zn含量的映射关系，实现了Zn含量高效快速反演，一定程度上揭示了喀斯特地区土壤重金属Zn的光谱吸收特性。程先锋等[12]通过构建土壤中Zn、Pb、砷（As）、镉（Cd）4种重金属含量及土壤光谱之间的多元线性回归模型，证明土壤重金属含量与光谱特征具有较强的相关性。宋婷婷等[13]建立了基于乘积变换的Zn含量最优预测模型并制作了污染分布图，表明Zn污染和矿业活动间具有显著相关性。TAN等[14]在伊通县矿区对选定的光谱特征和土壤样品的重金属数据进行反演并建立反演模型，结果表明CARS叠加法能较好地预测研究区As、铬（Cr）、Pb和Zn的含量。LIN等[15]利用黑龙江省黑土样本的室内光谱数据，分析了土壤有机质、铁（Fe）、Zn含量与6种光谱反射率的相关性；然后在α=0.01水平进行相关系数显著性检验，采用主成分分析（KPCA）算法，结合随机森林（RF）和支持向量机（SVM）组合模型应用于Zn、Fe含量的定量反演。结果表明在土壤重金属Zn、Fe的预测中，KPCA-RF模型优于RF模型，其在土壤成分含量的定量反演中具有广阔的应用前景。

目前，对于土壤Zn的研究，国内外多侧重于矿区Zn的污染及土壤Zn的回归模型的精度与迁移能力，而较少关注耕地质量提升与治理中Zn的反演研究，从而制约了国家土壤健康及高标准农田建设的稳步推进。因此，本研究以黄淮海典型区域（西平县）高标准农田土壤重金属Zn作为研究对象，于实验室测得野外采集土样的Zn含量，并利用ASD Field Spec便携式光谱仪采集土壤光谱数据；运用5种方式对反射率光谱进行变换以增强光谱特征信息，并与土壤Zn含量实测值进行皮尔逊（Pearson）相关分析，初次筛选出显著特征波段；再利用连续投影算法（SPA）二次筛选出最佳特征波段，以构建Zn含量偏最小二乘法回归模型（PLS）；检验5个模型的均方根误差（RMSEC、RMSECV、RMSEP）、校正决定系数（R2）和相对分析误差（RPD）并进行比较，选取Zn含量最佳反演模型。本研究旨在为预测高标准农田建设区域土壤重金属Zn含量及其治理提供理论参考。

1 材料与方法

1.1 研究区概况

驻马店市西平县位于北纬33°10'～33°32'，东经113°36'～114°13'（图1），河南中南部区域，其土壤类型分布及面积占比见表1。县境总面积为1 089.77 km2，划分为3个街道办事处、8镇、8乡和老王坡管理委员会等行政区。该地气候属亚湿润大陆性季风型，全年气候温和。西平县资源丰富，农作物物种丰盛，是我国粮食生产先进县与粮食主产区之一，同时其西部山区拥有森林覆盖率95%以上的国家级森林公园。

图1 研究区位置示意图Figure 1 Location diagram of the study area

表1 西平县土壤类型分布及面积占比Table 1 Distribution of soil types and area share in Xiping County

1.2 土样采集

土样采集前对西平县研究区的土壤类型、空间特征、行政单元和地形因素进行分析，以设定精确定位的采样点。利用对角线采样法采样，野外采样深度为0～30 cm土壤层，同时记录详细地形特征。从田里采集的土样及时摊开风干，倒在干净的容器内，去除石块、植物根茎残体及其余杂物。取研磨均匀后通过1 mm孔径孔筛的土壤样本，利用四分法取样，用于土壤光谱数据的采集和Zn含量的测定[16]。

1.3 土样光谱室内测定

1.4 光谱变换方法

倒数对数光谱变换可以增强光谱反射率数据在可见光区域的特征差异，并能减少因光照强度和地形起伏变化等因素对光谱特征曲线引起的负面影响[18]，其变换公式见表2。

光谱微分技术是一种简单而有效的确定光谱曲线极值点波长位置的方法，是通过对反射率光谱数据进行阶数不同的数学微分模拟来提取特征光谱参数，该技术能够消除平缓背景部分环境或基线漂移的影响，扩大样本之间的光谱特征差异，有助于吸收特征的提取[19]（表2）。

包络线是利用直线连接每个随波长不断变化的吸收谷或反射峰凸出的特征点，并使连接线的外角大于180°的线[20]。去包络线是以相应波段上实际光谱与包络线两者反射率值进行比值计算，从而使之归一化[21]（表2）。

表2 光谱变换公式Table 2 Spectral transformation equations

1.5 相关分析与显著性检验

相关分析是用于判断多个随机变量间是否存在统计学关联的有效统计分析方法[22]，且能深入分析存在关联的变量间的关联强度与方向，反映事物间的统计学相关程度。两个变量之间相关关系的强弱程度，可以用数值的形式揭示，其数值又称为相关系数。皮尔逊相关系数（积差相关系数）是应用最广泛的相关系数，其公式如下：

式中：xi、yi（i=1，2…，n）为两个随机变量x和y的样本值；xˉ、yˉ为两个随机变量样本值的平均值；n为两个变量样本值的个数。R为相关系数，且-1≤R≤1；||R越接近1，说明两变量间线性相关性越强，||R=1，为完全线性相关越接近0，说明两变量间线性相关性越弱，为完全不具有相关性；R＞0，说明两变量为正相关；R＜0，说明两变量为负相关。在科研应用中，可以划分几个区间取值来表示不同的相关程度。

样本间是否存在显著线性相关关系不能直接通过相关性分析得到的系数值确定，其原因在于随机抽样的性质或样本数量的不足等，因此必须使用假设检验的方式对总体的显著性进行相关分析显著性检验[23]。

1.6 连续投影算法

连续投影算法（SPA）是将向矢量空间共线性最小化的前向变量选择算法[24]。SPA通过投影将光谱数据的波长映射到其他波长上进行向量的投影分析，从而比较出所投影向量的大小，选择投影向量值最大的波长作为优选波长，然后经过矫正模型选取出所需要的特征波长[25]。

1.7 偏最小二乘法回归模型反演与精度检验方法

在构建偏最小二乘法回归模型（PLS）[26]中，从自变量中提取的主成分要求尽可能体现出光谱曲线的原特征信息[27]，另外从因变量中提取的独立成分要求与主成分存在最大的相关性。将提取出的成分与因变量通过多元回归方法构建回归方程[28]。RMSE、R2和PRD[29]的公式及定义见表3。

陈树华所说的黄龙病，由一种限于韧皮部内寄生的革兰氏阴性细菌，是世界柑橙生产上的毁灭性病害。“黄龙病虽不可治，但是可防可控。”陈树华介绍。

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

2 结果与分析

2.1 土壤Zn属性特征

土壤重金属Zn含量范围在57.13～306.40 mg·kg-1之间（表4），变异系数[30]为0.401，为中等变异性，适宜构建模型。Zn元素的质量分数偏态分布密度曲线表现出土壤Zn含量多集中于125～130 mg·kg-1之间（图2），极端值分布较少。《土壤环境质量 农用地土壤污染风险管控标准（试行）》（GB 15618—2018）规定[31]，当土壤pH≤6.5时，土壤重金属Zn含量的风险筛选值为200 mg·kg-1。调查研究区土壤平均pH值为5.9，样本数据表明该区整体Zn污染风险较低，少数样本极端值表明该区局部Zn污染风险较高，应当加强土壤环境监测和农产品协同检测。

图2 Zn质量分数描述统计图Figure 2 Descriptive statistics of Zn mass fraction

表4 土壤重金属Zn含量分析Table 4 Analysis of soil properties for heavy metal Zn content

2.2 土壤属性光谱反射特征

土壤光谱曲线一般为突出弧线，土壤样本间的光谱差异性主要体现为相同波段上的光谱反射率数值不同（图3）。土壤样本在可见光带400～780 nm范围处的反射率较低，但呈大幅度上升的趋势，光谱曲线陡峭，土壤有机质在600、800 nm附近有微弱的反射峰，在780～2 300 nm处近红外波段区域曲线较平缓，在2 300～2 400 nm处逐渐下降。1 850～2 140 nm光谱曲线呈现明显的“V”字形。图中3个水分吸收带表现出的波谷与黏土中的多种矿物质组成有关，1 420 nm为羟基（OH）谱带、1 920 nm附近为H2O谱带、2 200 nm波段范围为羟基伸缩振动与和OH弯曲振动的合谱带[32]。4条平均光谱曲线整体趋势相似，丘岗地区的平均光谱曲线较低，山地丘陵区的平均光谱曲线最高，而县区与平原区的平均光谱曲线居中，几乎吻合。

图3 西平县不同地类土壤平均光谱反射率曲线Figure 3 Average spectral reflectance curves of different land types of soil in Xiping County

2.3 相关分析与显著性检验

从5种光谱变换后的光谱数据与土壤Zn含量的相关系数折线图（图4）可以看出，经过光谱变换，各波段与土壤Zn含量的相关系数也有所不同，呈现多元化的变化趋势。倒数对数光谱变换（LOG）与平滑光谱变换（SG）的相关系数曲线较为平滑，且呈现以数值0为对称轴的对称显示，显著相关系数的波段数、最大相关波段基本一致；二者中一者递增，另一者则递减，最显著波段相关系数呈现相反关系，其绝对值十分相近，并均接近0.3。

图4 5种变换光谱反射率相关系数显著性曲线Figure 4 Significance curves for the correlation coefficients of the five transformed spectral reflectances

去包络线光谱变换（CR）的相关系数曲线出现上下频繁波动，在600～1 400、1 900～2 400 nm范围内有多处波峰和波谷，数值变化范围介于-0.3～0.3，曲线变化形似于一阶微分变换（FD）后的光谱相关系数曲线。SG光谱反射率相关系数曲线与CR光谱相关系数曲线拥有相似的波峰位置。FD与二阶微分变换（SD）光谱的相关系数数值变化折线图不再呈现近似单一数值的平滑线，而是表现为在正负值上下频繁的波动，FD光谱的相关系数曲线在580 nm附近出现一个波峰，在1 900 nm附近出现一个波谷，SD的相关系数曲线变化较为剧烈，连续波段相邻的相关系数在0上下跳跃。

FD、SD、CR变换光谱与Zn含量显著相关的波段数量较SG交换光谱的显著相关波段数明显减少，其中SD变换光谱相关系数曲线的显著相关波段数最少。5种变换光谱对应的最显著波段中，SD变换光谱的相关系数最大，FD变换和CR变换光谱的相关系数次之，LOG变换和SG变换光谱对应的相关系数最小（表5）。

表5 光谱变换后反射率与Zn含量显著相关的波段Table 5 The wavelengths with significant correlation between reflectivity and Zn content after spectral transformation

综上所述，经过LOG、FD、SG等5种光谱变换后的光谱反射率能有效突出光谱反射率的变化特征，经检验光谱变换下通过α=0.01水平显著性检验的相关波段可作为第一次特征波段数据，用于进行第二次特征波段的筛选，从而构建土壤Zn含量的反演模型。

2.4 最佳特征波段的选取

利用Matlab软件结合连续投影算[33]法分别对SG、LOG、FD、SD和CR变换后的光谱数据进行波段的第二次筛选，得出最佳特征波段见表6。

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

实验没有考虑其他重金属元素会干扰光谱反演分析过程，本研究认为每种重金属都有很多敏感波段且不尽相同，但在通过优选特征波段的组合后可逐步降低与其他重金属敏感波段的重叠程度。通过第一次筛选和第二次筛选后，将从300～2 500 nm整谱波段中筛选出的5～10个最佳特征波段进行组合，使光谱数据与Zn含量的相关性进一步上升，很大程度上降低了与其他重金属的关联性。

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

Zn分别以FD、LOG、SG、SD、CR光谱变换为自变量的模型[34]整体上具有较高的Rˉ2，其值均大于0.60，最高为0.70，RMSEC都在20范围上下（表7）。Zn以SG光谱变换为自变量的模型交叉验证系数Rˉ2最小，为0.61，以CR光谱变换的交叉验证系数Rˉ2最高，为0.68。

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

对比建模集结果和交叉验证结果，Zn含量以CR光谱变换为自变量的模型建模效果较好。检验的均在验证集检验的结果中有不同的变化，但建模效果比较好的模型在验证集的验证中精度仍然较高。进行RPD数据分析发现，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种模型中的最优模型。

图5 土壤Zn含量实测值与预测值的拟合散点图Figure 5 Fitting scatter plots of measured and predicted soil Zn content

3 结论

（1）研究区土壤Zn含量呈一定程度的偏态分布，大部分土壤样品Zn含量介于60～160 mg·kg-1，均值为123.86 mg·kg-1，表明研究区整体Zn污染风险较低；但存在极端值大于300 mg·kg-1的样品，说明研究区内局部地点Zn污染风险高。

（2）通过光谱反射率平滑处理消除仪器不同探测元件对光谱数据的影响，使反射率曲线变平滑，提高信噪比。通过光谱变换即倒数对数（LOG）、一阶微分（FD）、二阶微分（SD）和去包络线（CR）等方法提高提取的特征波段的准确性和模型的预测能力，有效突出光谱曲线的吸收和反射特征，从而提高了光谱的灵敏性并增强了有用信息。

（3）以5种反射率变换光谱数据为基础，提取出最佳特征波段作为模型自变量，对应土壤Zn含量作为模型因变量，构建了基于偏最小二乘法的西平县高标准基本农田建设区域土壤Zn属性反演模型。土壤Zn的最佳模型是以CR光谱变换为最佳的偏最小二乘模型。

本研究仅比较了通过不同光谱变换而挑选出的特征波段经同一建模方法建立反演模型的精度，没有进行多种光谱变换方式、多种不同建模方式间的组合比较，土壤光谱数据还有待深入挖掘。此外本模型精度有待进一步提高，且当前模型仅适用于研究区，而对其他地区的适应性有待研究验证。