基于遥感数据和随机森林算法的土壤重金属空间分布模拟

摘要 为评估土壤重金属的富集状态及空间分异态势，选取山东省章丘市为研究区，系统采集425处土壤样品，测定土壤中铬（Cr）元素含量，采用描述性统计特征评估重金属在土壤中的富集状态；获取与土壤采样同期的Landsat-8 OLI遥感数据，将土壤重金属的环境要素作为自变量，测定的土壤Cr元素含量为因变量，构建基于随机森林算法的土壤重金属空间模拟模型，完成土壤中的重金属含量预测和空间分布模拟。结果表明，土壤重金属Cr含量均值高出土壤元素背景值37.22%，但低于农用地土壤污染风险筛选值，表明土壤中Cr的富集在可管控范围内；随机森林算法支持的空间模拟模型具有较好的精度和稳定度，精度系数R2和RMSE值分别为0.87和7.19，优于普通克里格法（R2=0.66，RMSE=13.15）对土壤重金属的空间分布模拟。

关键词 土壤重金属；随机森林算法；遥感；空间分布模拟；铬

中图分类号 X53  文献标识码 A  文章编号 0517-6611（2023）14-0051-04

doi：10.3969/j.issn.0517-6611.2023.14.013

基金项目 环境演变与自然灾害教育部重点实验室开放课题（2022-KF-14）。

作者简介 周忠科（1997—），男，山东枣庄人，从事生态遥感研究。*通信作者，副教授，硕士，从事生态旅游研究。

收稿日期 2022-07-30

土壤是地表生物正常进行各类生命活动的基础和媒介。重金属是土壤的重要组成部分，因其具备毒性大、富集性强、残留周期长的特点，极易对人类和动植物造成危害。随着我国改革开放以来多年的工业化发展，土壤重金属污染愈发严重［ 1-3］。土壤中的重金属元素会污染正常耕作的耕地土壤，威胁动植物正常生长，并在食物链的累积作用下危及人类的健康［ 4-9］。

实现对土壤重金属的含量特征和空间分布的监测是开展土壤污染治理和修复的关键。由于土壤中重金属元素的种类、含量不是一成不变的，是不断在累积的，因此传统方法便是在实验室内不断对土壤重金属含量、种类进行调查和检测。例如，贾佳瑜等［ 10］采用室内实测和地统计法相结合的方式对土壤重金属含量及其空间分布特征进行分析；黄志伟等［ 1］对表层土壤重金属含量开展实地调查的方法，进行土壤重金属污染特征及潜在风险评价；SONG等［ 11］利用地理信息系统的方法，来绘制重金属空间分布热点图，以此来建立空间分布与工业布局之间的相互关系，以判别土壤重金属的空间分布受到人类活动干扰的情况。上述研究，较好地完成了对研究地区土壤重金属含量特征和空间分布评估，但是此种方式费时费力，测试过程中也产生污染废物，会对测试人员和周边环境造成二次污染。所以寻找一种高效、准确的调查与监测方法是非常有必要的。近些年，随着遥感技术的快速发展，卫星影像的空间分辨率和对地物的光谱敏感性有了显著的提升，为各种环境相关信息进行光谱定量分析提供了数据支撑和技术支持［ 12］。尽管重金属元素在土壤中含量较低，直接建立重金属与卫星光谱数据之间的联系较为困难；但是可以利用基于卫星数据获得的多种环境要素信息（例如高程、坡度和植被指数等），间接地表征土壤重金属的含量和富集状态，并以此为输入变量，进行区域建模，完成土壤重金属空间分布模拟。相对比传统方式，土壤重金属污染的监测调查更为方便快捷，同时节省大量的费用和成本。

笔者选取山东省章丘市为研究区，对 425 处土壤样品进行系统采集，并对土壤中Cr元素含量进行化学分析，通过分析重金属含量的描述性特征，评估重金属的富集状态，进一步建立基于随机森林算法的空间模拟模型，实现对土壤重金属的含量预测和空间分布模拟，为保护周边居民身体健康、动植物健康生长和土壤修复提供技术和理论支持。

1 材料与方法

1.1 研究区概况

章丘区位于山东省中部，地理坐标在 36°25′～37°09′N 、 117°10′～117°35′E，是省会济南市的一个辖区，总面积约为1 719 km2（图1），人口约107.6万。章丘地形南高北低，地形分别为山区、丘陵、平原、洼地，其中山区和丘陵占比达到 56.7%左右，黄河流经北境。章丘区属于温带半湿润大陆性季风气候，四季分明，雨热同季，春旱多风，夏季多雨，秋季温爽，冬季干冷。年平均气温12.8 ℃，年降水量600.8 mm。章丘区内铁路和公路横纵贯通，交通便捷［ 13］。区域内涵盖多种产业的发展，尤其是第二产业交通装备、机械制造、生物食品和化工材料等，为区域经济发展奠定了坚实的基础［ 14］。

1.2 样品采集与实验室分析

表层（0～20 cm）土壤样本在章丘区实地采集，在选取和设计时，综合考虑研究区的土地利用和水文地质条件，基于多点取样、等量混合、与抽样统计原理相结合的原则，并结合已知的道路信息，确保顺利完成采样。在每个采样点，按照100 m范围为半径，将表层土壤混合至1 kg左右的综合土壤样本，密封到聚乙烯塑料袋中，送往实验室待测［ 15-16］。同时利用手持GPS，记录采样点的真实坐标，将完成的425处采样点做好记录。送回实验室的土壤样本，简单去除杂草、根、枯枝残留物和小石子等明显杂质，将样品放置在25 ℃的条件下自然风干，其后采用研钵进行研磨，过0.2 mm筛；最后在实验室内采用石墨炉原子吸收光谱法分析检测425份样品的Cr含量，技术规范和回收率符合预期监控要求。

1.3 多源环境遥感数据及预处理

该研究采用的环境遥感数据主要包括高程、坡度、坡向、黏土矿物比值（CMR）、改进的归一化差异水体指数（MNDWI）、 归一化差异植被指数（NDVI）、改进的土壤调节植被指数（MSAVI）、增强植被指数（EVI）和红绿比率指数（RGRI）。高程、坡度和坡向的原始数据获取自ASTER GDEM数字高程的数据产品（下载自地理空间数据云https：//www.gscloud.cn/），其后采用ArcGIS 10.2的分析工具，完成区域的高程数据图、坡度数据图和高差斜坡数据图［ 17］。其余环境指数的计算是基于Landsat-8多光谱遥感数据，原始数据免费从美国地质勘探局网站（https：//glovis.usgs.gov）获取得到，其后经过图像裁剪、辐射校正、辐射亮度值的提取以及光谱信息加强等操作后，运用数学公式（1）～（6）计算得到相应的环境指数［ 18-19］；最后取样点的425个经纬度坐标分别导入相应的高程数据图、坡度数据图、高差数据图和各类环境指数图，在上述数据图中分别完成提取，得到各个采样点对应的值。上述数据的预处理和求算均在ENVI 5.3.1软件中完成。

CMR=SWIR1SWIR2（1）

MNDWI=Green-SWIR1Green-SWIR2（2）

NDVI=NIR-RedNIR+Red（3）

MSAVI=2×NIR+1-（2×NIR+1）2-8×（NIR-Red）2（4）

EVI=2.5×（NIR-Red）NIR+6×Red-7.5×Blue+1（5）

RGRI=RedGreen（6）

式中，Blue波段、Green波段、Red波段、NIR波段、SWIR1波段、SWIR2波段分别指Landsat多光谱影像中的第2、3、4、5、6和7波段；CMR是指通过2个短波红外波段的比值，环境意义表示的是岩石或土壤中黏土矿物含量情况；

MNDWI是可以很容易地区分建筑用地背景下的阴影和水体，更好地提取水体特征。利用MNDWI解决了在建筑背景下提取水体信息效果不佳的问题，从而使建筑背景下的水体信息提取更精确。NDVI作为一种常见的植被指数，其取值区间为［-1，1］，-1、0、1分别表示地面覆盖为水、雪或有云的天空，地表有岩石或裸土等情况，对可见光的反射率较高；当地表有植被覆盖时，且覆盖密度越大，NDVI越大；但NDVI也具有缺点，就是对高植被覆盖地区的灵敏度不高。MSAVI克服了传统的土壤调节植被指数难以适应不断变化的土壤条件，该指数能够很好地解释背景的光学特征变化，并且土壤背景的灵敏度也在一定程度上得到提高。EVI可以降低大气和土壤噪声以及水汽的影响，对于检测区的植被情况具有一个稳定的反映；用于植被茂密区，但对植被稀疏区也具备一定的探测能力，工作原理是通过将红光和近红外波段的范围设定变得更窄来实现，与此同时引入蓝光波段来纠正气溶胶的散射与土壤背景。红绿比率指数（RGRI）是指将红波段与绿波段相比较，得到一个比值，可以显著区别茂密与稀疏林地。为便于接下来土壤重金属空间模拟模型的建立，所有环境指数图在建立之后，均进行30 m分辨率的空间重采样。

1.4 空间模拟模型的建立与精度评估

土壤重金属与遥感数据之间有着极其复杂的关系，选择与土壤重金属最为密切的遥感数据作为模型输入自变量，有助提高模型的效率和计算精度。皮尔逊相关系数（Pearson correlation coefficient）是2个变量的协方差与标准差之间的一个比值，可以表示2个变量之间存在的相关关系，应用广泛［ 20］，其公式如下：

r（x，y）在［-1，1］之间取值，当其值为正时，表示正相关；当其值为负时，表示负相关；当其值为0时表示不相关。因此，r的绝对值越大，说明变量之间存在的相关性越强，否则相关性就越弱。该研究采用皮尔逊相关系数分析遥感数据与土壤重金属之间的相关性。

随机森林模型（random forest，RF）是由一定数量决策树整合而成的集成模型的一种，通过对策树投票或平均数据集的划分，得到了分类或回归的随机森林输出结果［ 21］。随机森林算法具有训练速度快、预测精度佳和适合不平衡数据集的特点［ 22-24］。决策树数量和分裂节点数目是构建基于随机森林算法的土壤重金属空间模型的重要参数，经过多次测试，分别设置为1 000和3。

采集的425个土壤样品采用随机抽样法选取383份数据用作训练集来构建随机森林模型，剩下42个数据作为独立数据集用来评估模型的精度。决定系数（R2）和均方根误差（RMSE）是评价空间模拟模型的常用指标，R2表示因变量被自变量完全解释信息的程度，RMSE表示模型估算值与实测值差值的大小，较大的R2和较小的RMSE表示模型具有更好的预测精度和稳定性［ 25］。计算公式如下：

2 结果与分析

2.1 土壤重金属的描述性统计特征

从表1可以看出，土壤中Cr含量超出元素背景值［ 26］37.22%，但是没有超出《土壤环境质量农用地土壤污染风险管控标准》（GB 15618—2018）的筛选值，检测结果均在风险管控值范围之内。值得注意的是，Cr元素的中位数值超过背景值，Cr元素可能为该区域内具有较高等级的潜在风险元素。变异系数主要用来评价数据离散程度的相对大小［ 15］，该研究测得Cr变异系数为0.03～0.06，离散程度较小，属于低度变异。峰度和偏度主要用来衡量数据的分布状态，该研究的总样本集、建模集和验证集的偏度值均大于1，处在正偏状态，表明土壤中Cr元素含量的分布受到不同程度的外部因素干扰。

2.2 土壤重金属与多源环境遥感变量相关分析

研究区域的地质地貌、植被覆盖、气候条件和人类活动对土壤中重金属元素含量的影响很大，选取建模因子对反演模型质量的好坏影响很大。在Cr元素建模时需要综合考虑光谱波段、衍生光谱指数和地形等因素，这样建立的反演模型得出的反演结果最优。综合考虑光谱波段、光谱衍生指数和地形因素，在数据齐全的情况下，筛选出相关系数最高的因子，结果表明，Cr和高程、坡度、坡向、近红外、短波红外2、增强植被指数和短波红外1的相关系数分别为0.231、-0.194、-0.153、-0.112、0.090、0.077和0.060，均具有明显相关性，因而选为土壤重金属的敏感环境变量，用于模型的输入自变量。

2.3 空间分布模拟模型建立与对比

以选择得到的重金属敏感的环境要素作为输入变量，测得的Cr元素含量为因变量，分别训练得到随机森林和普通克里格估算模型，用于土壤重金属的含量预测和空间分布模拟。模型在验证集中预测的效果如图2所示。从图2可以看出，普通克里格模型（R2=0.66 和RMSE=13.15）的估算值和实测值基本保持在1∶1 线附近，具备一定的估算能力，但普通克里格在估算过程中出现高估和低估误差偏离点较多。随机森林模型的重金属预测精度指标R2和RMSE分别为0.87和7.19，相对普通克里格模型R2上升了31.82%，RMSE 下降了45.32%；随机森林模型的估算值与实测值在标准1∶1 线的两侧，预测值与实测值偏差较小，估算效果最佳。