环境科学  2025, Vol. 46 Issue (6): 3902-3912   PDF    
引用本文
王玮, 王政, 孔祥意, 吕建树. 基于INLA-SPDE模型的区域土壤硒元素空间预测及富硒区优选[J]. 环境科学, 2025, 46(6): 3902-3912.
WANG Wei, WANG Zheng, KONG Xiang-yi, LÜ Jian-shu. Spatial Prediction of Selenium in Soils by Using INLA-SPDE Approach and the Delimitation of Selenium-enriched Land with Low Heavy Metal Pollution Risk[J]. Environmental Science, 2025, 46(6): 3902-3912.
基于INLA-SPDE模型的区域土壤硒元素空间预测及富硒区优选
王玮1, 王政2, 孔祥意1, 吕建树1     
1. 山东师范大学地理与环境学院，济南 250300;
2. 浙江大学资源与环境学院，杭州 310058
收稿日期: 2024-05-11; 修订日期: 2024-08-01
基金项目: 山东省自然科学基金优秀青年基金项目（ZR2020YQ31）
作者简介: 王玮（1999~），女，硕士研究生，主要研究方向为土壤污染物空间建模，E-mail：wangwei172211@163.com
通信作者: 吕建树, E-mail：lvjianshu@126.com
摘要: 土壤是动植物和人类获取生长发育所需硒元素（Se）的主要介质，但经济的快速发展造成了严重的土壤重金属污染，土壤环境和人类身体健康受到威胁. 因此，预测土壤Se空间分布，划定低重金属污染风险的“清洁富硒土地”，可为工农业健康发展提供参考. 以淄博市部分地区为研究区，系统采集212个土壤样品进行Se及重金属元素测试；采用集成嵌套拉普拉斯逼近-随机偏微分方程（INLA-SPDE），整合距道路距离和距工厂距离、土壤类型、降雨量和土壤粒度这5种环境协变量，对土壤Se元素空间分布进行预测，并结合有限混合分布模型（FMDM）确定的污染阈值，进行富硒区优选. 结果表明：①研究区土壤ω（Se）平均值为0.34 mg∙kg-1，高于山东省土壤Se背景值和全国土壤Se含量平均值. ② INLA-SPDE方法预测结果显示，土壤Se含量空间变异与自然因素和人类活动相关；富硒土地集中分布研究区中部，面积为423.49 km2. 空间预测标准差分布表明，采样点密度大，预测不确定性小. ③依据FMDM模型提供的污染临界值，划定低于中污染阈值和高污染阈值的富硒土地面积分别是2.22 km2和291.81 km2. 富硒土地的优选可为持续开发、土壤环境管控和人类健康保护提供参考.
关键词: 富硒土地      集成嵌套拉普拉斯逼近-随机偏微分方程      土壤重金属      低重金属污染富硒区      淄博市     
Spatial Prediction of Selenium in Soils by Using INLA-SPDE Approach and the Delimitation of Selenium-enriched Land with Low Heavy Metal Pollution Risk
WANG Wei1 , WANG Zheng2 , KONG Xiang-yi1 , LÜ Jian-shu1     
1. College of Geography and Environment, Shandong Normal University, Jinan 250300, China;
2. College of Environment and Resource Science, Zhejiang University, Hangzhou 310058, China
Abstract: Soil is the main medium for plants, animals, and humans to obtain the selenium element required for growth and development. Rapid economic development has caused serious soil heavy metal pollution, threatening soil environment and human health. Therefore, predicting the spatial distribution of soil selenium elements and delineating "clean and selenium-rich land" with low heavy metal pollution risk can provide a reference for the healthy development of agriculture and industry. A total of 212 soil samples were collected in part of Zibo City for selenium and heavy metals testing. The integrated nested Laplace approximation-stochastic partial differential equation (INLA-SPDE) method was used to integrate five environmental covariates, including distance to roads and factories, soil type, rainfall, and soil particle size, to predict the spatial distribution of soil selenium elements. Combined with the pollution thresholds determined by the finite mixture distribution model (FMDM), the selenium-rich areas were optimized. The results show that: ① The average selenium content in the study area was 0.34 mg∙kg-1, higher than the background value of selenium in Shandong soil and the national average. ② The INLA-SPDE method prediction results showed that the spatial variation of soil selenium content was related to natural factors and human activities; the selenium-rich land was concentrated in the middle of the study area, with an area of 423.49 km2. The distribution of spatial prediction standard deviation showed that the sampling density was high, and the prediction uncertainty was small. ③ According to the pollution threshold provided by the FMDM model, the areas with selenium-rich land lower than the medium pollution threshold and high pollution threshold were 2.22 km2 and 291.81 km2, respectively. The optimization of selenium-rich land can provide a reference for sustainable development, soil environment control, and human health protection.
Key words: selenium-enriched land      integrated nested Laplace approximation-stochastic partial differential equation      soil heavy metals      selenium-enriched land with low heavy metals pollution      Zibo City     

硒元素（Se）是人类以及动植物生长发育所必须的元素[1~5]. 土壤-植物系统是人类获取Se的主要途径，而土壤中Se含量的空间分布存在异质性，合理开发富硒土地成为解决Se含量空间差异的重要方式之一[6~10]. 但随着经济的快速发展，土壤环境不仅受自然因素的影响，同时受到人类活动的干扰. 据《全国土壤污染状况调查公报》，土壤污染类型以重金属无机污染为主，重金属污染日趋严重[11~15]. 土壤中的重金属元素会通过食物链进入人体，威胁人类身体健康. 同时，Se和重金属元素的协同作用可能会使植物吸收更多的重金属. 因此，依据土壤重金属污染程度进行富硒土地优选，对富硒特色农业的可持续发展具有重要的意义.

地统计学是用于刻画空间变量异质性的范式性工具，包括克里金和随机模拟两类[16，17]. 克里金是最早且应用最广泛的预测变量空间分布的方法[18~20]. 然而，克里金在预测过程中产生的平滑效应限制了空间预测精度的提升[21]. 随机模拟方法通过实现多个等概率实现克服了克里金的不足，并可进行不确定性分析[22~24]. 在土壤中，Se的空间分布受到多重环境因素的影响[25~27]，因而需要将环境因素作为协变量引入到预测算法中，以期能够获得更好的预测效果并完成不确定性分析. 贝叶斯地统计学方法中的集成嵌套拉普拉斯逼近-随机偏微分方程（INLA-SPDE）则可整合环境因素作为辅助变量对化学元素进行空间预测和不确定性分析，并能够提供固定效应参数和随机效应参数，为化学元素来源提供参考.

山东省淄博市的淄川区、博山区和周村区富硒土地资源丰富，富硒特色农业快速发展，同时也是具有百年历史的老工业区. 农业和工业等人类活动聚集必然对土壤环境产生强烈的影响. 因此，本研究在此范围内系统采集212个表层土壤样品（0~20 cm），测定了土壤Se和As、Cd、Co、Cr、Cu、Hg、Mn、Ni、Pb和Zn的含量，采用集成嵌套拉普拉斯逼近-随机偏微分方程（INLA-SPDE）方法，集成环境辅助变量，对研究区土壤Se空间分布进行预测，依据硒含量划分低硒、中硒和富硒土地. 结合有限混合分布模型（FMDM）确定的土壤重金属中污染和高污染阈值，对研究区富硒土地进行优选，划分出低重金属污染富硒土地.

1 材料与方法 1.1 研究区概况

研究区主要包括位于山东省淄博市中部的淄川区、博山区和周村区[图 1（a）]，地理位置介于117°40′~118°10′E，36°10′~36°50′N之间，总面积约为1 800 km2. 研究区整体地势呈现出南高北低的特征，北部及中部主要为平原，南、东和西部以山地丘陵地形为主. 土壤以粗骨土和褐土为主，还包括棕壤和河潮土[图 1（a）]. 成土母质类型包括泥岩、灰岩、白云岩、砂岩、花岗岩、砂砾岩和第四纪洪冲积物[图 1（b）]. 土地利用类型主要有耕地、林地、草地、水体和人造地表[图 1（c）]. 研究区为温带大陆性季风气候，年平均气温为13.7℃，年平均降水量为660 mm. 研究区工业发展历史悠久，部门齐全，主导产业为建材冶金、机械制造、医药化工和纺织服装，区域内分布有2座燃煤发电厂和多个工业园区[图 1（a）]. 此外，研究区内依托天然富硒土地资源，大力发展富硒特色农业和采摘园经济，并建立起了大片集中的果园. 淄川区和博山区因富硒农业的发展分别获得了“全国富硒农业示范区”和“中国富硒农业示范基地”的荣誉称号.

图 1 研究区位置、成土母质和土地利用概况 Fig. 1 Location of the study area with sampling sites, soil types, parent materials, and land use types

1.2 土壤样品采集与分析测试

本研究采用k-means聚类方法进行室内采样点设计[28~30]，选取成土母质、土壤类型、降雨量、土地利用、植被归一化指数、海拔、复合地形指数、坡度、距工业点距离和距道路距离等10种环境变量进行空间聚类. 将上述变量重采样为100 m×100 m的栅格. 利用k-means聚类方法建立220个聚类单元，聚类的数量和样本的数据相等. 在聚类过程中，环境变量被缩放，方差逐渐变为1，采样点最终确定为环境变量和聚类单元中心距离最短的栅格单元. 为验证获取的采样点的代表性，采用卡方检验的方法，将采样点的位置的协变量和原协变量进行验证. 显著性均大于0.05，说明采样点具有较好的代表性. 但在实际采样的过程中，由于受到交通可达性和经济成本等因素的限制，在尽可能靠近预设采样点的原则下，最终获取212个土壤样品[图 1（a）]. 在每个点的采样过程中，采用多点采样法，在预设采样点周围的100 m的半径范围内，在样点周围获取3~5个土壤子样品，每个子采样品垂直采集0~20 cm深度的表层土壤，等量混合后放入聚乙烯自封袋中，样品原始重量大于10kg. 同时，记录采样序号、样品编号、距工厂、道路的距离以及土地利用类型采样点信息. 利用10 m精度手持式GPS确定采样点的实际坐标位置.

本研究将测试土壤中的As、Cd、Co、Cr、Cu、Hg、Mn、Ni、Pb、Zn和Se共11种化学元素. 具体测试方法如下：首先将样品在实验室内烘干，取样品50 g并用木棒敲碎，过2 mm筛以去除较大的砂砾，并混合均匀. 使用玛瑙研钵将过2 mm筛后的样品研磨至能够通过0.149 mm筛并取样品约0.2 g，使用HClO4-HNO3-HF进行消解. 利用电感耦合等离子体发射光谱仪（Agilent 5100）测试Cd、Co、Cr、Cu、Mn、Ni、Pb和Zn等元素；称取0.5 g样品，采用H2SO4-HNO3-KMnO4消解后，As、Hg和Se元素采用原子荧光光度计（AFS-2202）测定. 实验流程符合国家相关规定，所有元素的回收率均在（100±10）%之间，说明测试精度高，符合本研究的要求.

1.3 INLA-SPDE

贝叶斯方法最重要的是数据的后验边缘分布的计算. 虽然马尔科夫链-蒙特卡洛模拟（MCMC）可以作为后验分布计算的工具，但在计算收敛方面存在问题，并且难以快速处理大数据. INLA则可作为MCMC的替代工具，目的是在隐高斯模型中进行贝叶斯推理，并获得后验边缘分布. INLA采用解析近似和稀疏矩阵相结合的算法，以解析式的方式表达后验分布，这种方式计算速度快、逼近精度高，避免了样本收敛和混合的问题，适合大数据集的研究. 以本研究中土壤Se为例，同时假设研究区D为二维空间R2中的一个子集，点s在D中连续变化，随机变量X（s）表示为在s处的土壤Se含量，当在D中存在有n个点时，那么就存在有X（s1），X（s2），X（s3），…，X（sn ）组成含有n个随机变量的随机向量，空间自相关将随机变量联系在一起，构建起空间随机场. 当这个随机场表现为高斯分布时，可称之为高斯随机场. 假设土壤Se满足高斯分布条件，通过定义空间协方差函数和平均值即可建立一个高斯随机场. SPDE方法采用Matérn协方差函数模型，公式表达为：

(1)

式中，si - sj表示si和sj两个位的欧氏距离；X表示目标变量；v表示协方差的非负参数；σx2表示为方差；κ表示与变程相关的尺度参数；Κv表示第二类贝塞尔模型；Γ表示伽马函数；σe2表示噪声方差（noise variance），也可称为块金值；δij表示即克罗内克（Kronecker delta）.

当一个高斯场满足Matérn协方差函数模型，这个高斯场可称为随机偏微分方程的一个平稳解，随机偏微分方程可表示为：

(2)

式中，W（s）表示为白噪声；Δ表示拉普拉斯算子且，s1和s2表示二维空间中的坐标，α表示与平滑参数v有关的正整数且α=v+d/2. 此外，SPDE可将连续的空间随机变量X（s）近似为分段的线性函数，这可使计算量大大减少.

本研究所使用的INLA-SPDE方法借助R语言3.6.3版本（R语言开发团队，2017）提供的INLA包搭建（http://www.r-inla.org/）. 有关INLA-SPDE方法的更多理论细节可参考相关文献及INLA包的使用说明书[30~34].

1.4 有限混合分布模型

土壤重金属元素的含量可被视作由k个组分组成的混合分布：

(3)

式中，fi(z)表示随便变量的概率密度函数，可被视作由自然背景和污染分布组成的混合重金属数据. πi表示fi(z)的混合分布权重，需要满足以下约束条件：

(4)

式中，μi、σi和πi采用最大似然估计方法，并结合牛顿型（Newton-type）算法和EM算法进行估计.

代表污染阈值的截断值（cuttoff value）采用如下方法确定：

(6)

FMDM模型采用R语言3.6.3版本提供的Mixdist包进行计算.

1.5 模型预测性能评估

建模精度是对空间预测评估的重要方法. 在本研究中对研究区内212个采样点进行7∶3的随机划分，70%的数据用于建模，30%的数据进行精度验证. 主要使用平均绝对误差（mean absolute error，MAE）、均方误差（mean Square Error，MSE）和均方根误差（root mean square error，RMSE）来评价INLA-SPDE方法对土壤Se空间预测的精度，其公式分别为：

(7)
(8)
(9)

式中，n表示采样点数量；yi 表示真实值；表示预测值.

1.6 土壤硒等级划分标准

本研究参考2019年11月自然资源部中国地质调查局发布的《天然富硒土地划定与标识（试行）》（DD-2019-10）中富硒土地的阈值[35]和李嘉熙等[36]对低硒和中硒土地的划分阈值，对本研究富硒土地资源进行划分：大于0.3 mg∙kg-1划分为富硒土地；0.3~0.2 mg∙kg-1划分为中硒土地；小于0.2 mg∙kg-1划分为低硒土地.

2 结果与讨论 2.1 不同母质和土壤下土壤Se含量

成土母质是土壤元素的主要来源，是导致土壤化学元素呈现空间异质性的主要原因[37，38]. 不同土壤类型下，土壤Se含量同样会呈现出差异性. 表 1为研究区不同成土母质和土壤类型背景下土壤Se含量的描述性统计. 土壤ω（Se）的平均值为0.34 mg∙kg-1，是山东省A层土壤Se背景值的2.6倍，是全国A层土壤中Se含量的1.17倍[39].

表 1 研究区及不同成土母质和土壤类型土壤Se含量描述性统计1） Table 1 Descriptive data of selenium in soils in different parent materials and soil types

在研究区白云岩、花岗岩、灰岩、泥岩、砂岩、砂砾岩和第四纪冲洪积物这7种母质背景下，土壤ω（Se）的平均值分别为0.28、0.17、0.23、0.63、0.63、0.29和0.39 mg∙kg-1. 化学元素在土壤中的含量受到粒径效应的影响，土壤粒径越小，吸附能力越强[40]. 泥岩母质下的土壤粒径小，花岗岩因其抗风化能力强，土壤粒径大. 所以Se在泥岩背景下的土壤中含量最高，在花岗岩母质背景下含量最低. 在棕壤、粗骨土、褐土和河潮土这4种土壤类型中土壤ω（Se）的平均值分别为0.15、0.30、0.39和0.23 mg∙kg-1，在褐土中土壤Se含量最高，在棕壤中含量最少. 研究区内褐土覆盖的成土母质包括泥岩、第四纪冲洪积物和白云岩等岩石组成粒径较小的母质类型，而棕壤在研究区内的母质多为花岗岩和砂岩等抗风化能力较强的母质，土壤粒径较大，因而导致研究区内褐土中Se含量平均值高，棕壤中含量平均值最低.

2.2 INLA-SPDE方法土壤Se空间分布预测 2.2.1 土壤Se空间结构

土壤Se含量和经过正态变化的环境变量进行变异函数拟合，拟合参数如表 2所示. Se在空间中的结构可用高斯理论模型进行拟合，块金值为0.53，基台值为1.10，块金值和基台值的比值为0.48，表明Se属于中度的空间相关性，变程为24 491 m，残差平方和决定系数分别为0.03和0.96，表明拟合程度较高.

表 2 土壤Se变异函数拟合参数 Table 2 Variogram fitting for selenium in soils

2.2.2 INLA-SPDE空间预测模型

以往的地统计方法在土壤重金属空间预测方面已经取得了较大的发展，但是大多数方法无法同时实现整合环境变量的预测和不确定性分析. INLA-SPDE方法具有整合多个环境变量到空间预测算法中并可根据固定效应参数为元素来源提供参考的优势. 因此，本研究首先选取了距道路距离、距工业点距离、土壤类型、降雨量、植被归一化指数（NDVI）、海拔、坡度、复合地形指数（CTI）、母质、土地利用类型和土壤粒径这11种环境变量，经过多重共线性检验去除冗余变量后，最终使用了距道路距离、距工业点距离、土壤类型、降雨量和土壤粒径5种环境变量，引入到INLA-SPDE空间模型中. 图 2为土壤Se空间预测所建立的三角网格. 表 3为土壤Se空间建模中所得固定效应参数和超参数描述性统计. 结果显示距工业点距离固定效应参数的绝对值为2.27，在5种环境变量参数中最大，表明土壤中的Se含量与距工业点的远近有着密切的关系，距离越小土壤Se含量越高，2.5%分位数和97.5%分位数之差为1.60，这表明固定效应的不确定性较大. 边缘方差和块金方差的平均值分别为0.47和0.53，这与空间结构中块金值/基台值的比值所示结果相近，表现为中度空间相关性. 变程的平均值为3 104 m. 表 4为利用INLA-SPDE方法制图后的精度验证指标. 结果显示将环境变量纳入算法中，固定参数可以对目标变量的空间预测具有较高的预测精度.

图 2 INLA-SPDE方法所构建和使用三角网 Fig. 2 Constructed mesh for the INLA-SPDE approach

表 3 INLA-SPDE方法的固定效应的参数及超参数后验边缘分布描述统计 Table 3 Descriptive statistics of fixed effect and marginal posterior distribution of hyperparameters for INLA-SPDE approach

表 4 INLA-SPDE方法的空间预测方法精度检验 Table 4 Accuracy test of spatial prediction method for INLA-SPDE approach

2.2.3 土壤Se空间分布预测及富硒土地划定

图 3（a）为土壤Se含量空间预测分布. 总体上，土壤Se空间变异程度较大，ω（Se）范围介于0.13~0.80 mg∙kg-1之间. 土壤Se高值区呈东西方向集中分布于研究区中部，在研究区西北部王村镇和东南部的太河镇分布也有小范围的高值区. 泥岩在研究区中部横贯东西，泥岩背景下所形成的土壤颗粒较细，对化学元素具有较强的吸附作用. 研究区北部分布有东西方向的煤田，煤系地层是土壤Se含量高的重要因素. 此外，在工业及交通线集中的城区，同样分布有高值区，表明Se在土壤中的分布受到自然因素和人类活动的双重影响.

图 3 研究区土壤硒含量INLA-SPDE预测空间分布、土壤硒等级划分和空间预测标准差分布 Fig. 3 Distribution of selenium in soils predicted by INLA-SPDE approach, the spatial distribution of selenium contents with different levels, and the standard deviation distribution of spatial prediction

按照本研究使用的土壤Se含量等级划分标准，研究区大部分属于中硒土地[图 3（b）]，ω（Se）介于0.2~0.3 mg∙kg-1之间，面积约为1 323.5 km2，占研究区总面积的73.28%；其次为富硒土地，ω（Se）大于0.3 mg∙kg-1，面积约为423.49 km2，占研究区总面积的23.45%；低硒土地面积最小，仅59.09 km2，为研究区总面积的3.27%. 总体上，研究区富硒土地资源面积较大且成片分布，为富硒特色产业的发展提供了坚实基础.

图 3（c）为INLA-SPDE空间预测标准差分布，通常用来表征空间预测方法的不确定性. 结果显示，标准差的变化范围在0.372~0.745之间，采样点密度越大，标准差越小，空间预测的不确定性越小，这表明空间预测的不确定性与采样点密度密切相关.

2.3 结合FMDM污染临界值的低重金属污染富硒区划定 2.3.1 土壤Se与重金属的相关分析

Se在自然和人类活动的长期影响下，在土壤中会发生迁移、聚集和分散的效应，并和一些元素呈现出有规律的组合，而这种组合可作为识别来源和分析变化的重要参考. 表 5为Se与As、Cd、Co、Cr、Cu、Hg、Mn、Ni、Pb和Zn的相关系数. 结果显示，Se与Cd、Hg、Mn和Zn呈显著正相关. Cd和Zn可指示农药和化肥等农业活动[41~45]；土壤中Hg和Mn主要来源于煤炭燃烧和工业活动[46，47]，说明研究区在形成富硒土地的过程中也可能伴随着人类活动导致的重金属污染. 此外，Se与Cr同样呈显著相关，Cr通常与母质及成土过程相关，受人类活动影响小[48，49]. 同时，有研究表明煤系地层会导致Se在土壤中富集[26，50]，研究区内潟湖相石炭二叠纪煤系地层的分布使Se含量富集[51，52]. 相关分析结果表明，土壤中Se受到自然背景和人类活动的双重影响.

表 5 土壤Se与10种重金属皮尔逊相关性分析 Table 5 Pearson correlation analysis of soil Se and ten heavy metals

2.3.2 土壤重金属的FMDM拟合及污染阈值确定

土壤重金属含量数据的混合分布通常可由2~3个组分的对数正态分布拟合. 2个组分的对数正态混合分布由自然背景值分布和中污染分布组成；3个组分的对数正态混合分布由自然背景值分布、中污染分布和高污染分布组成. 而分布的组分个数则通过观察重金属含量的分布特征确定[53，54].

研究区10种重金属FMDM模型拟合结果（表 6和图 4），As、Co、Cr、Cu、Mn、Ni和Zn呈双组分对数正态分布，自然背景值组分占比在98.8%~83.2%之间，表明以上7种重金属的含量主要受到自然环境因素的影响. FMDM自然背景组平均值与山东省背景值和淄博市背景值相比：As、Cr、Cu和Ni这4种重金属背景值要高于山东省背景值但低于淄博市背景值；Co元素背景值要低于山东省背景值；Mn和Zn平均值高于淄博市背景值. 以上7种重金属的中污染临界值，均高于山东省背景值和淄博市背景值. 综合以上信息，As、Co、Cr、Cu、Mn、Ni和Zn在土壤中的含量总体上与自然环境相关，少量来自于人类活动等外部因子输入，较好地保持了自然本底含量.

表 6 FMDM拟合参数和截断值及山东省和淄博市土壤背景值1） Table 6 FMDM fitting parameters and cut-off values and soil background values of Shandong Province and Zibo City

图 4 As、Cd、Co、Cr、Cu、Hg、Mn、Ni、Pb和Zn的FMDM模型拟合结果 Fig. 4 FMDM fitting of As, Cd, Co, Cr, Cu, Hg, Mn, Ni, Pb, and Zn in soils

Cd、Hg和Pb呈3组分对数正态分布，自然背景组分占比分别为33.7%、31.3%和63.6%，中污染组分占比则达到58.9%、58.1%和29.4%，高污染组分占比在7.0%~10.6%之间，这表明以上3种重金属在土壤中的含量受到人类活动的深刻影响. 但是，通过观察Cd、Hg和Pb的中污染和高污染临界值发现，相同元素的2个临界值之间相差较小. 同时，中污染临界值虽然要高于山东省整体背景值，但是却与淄博市的背景值十分相近，而且该区域分布有较大面积的煤系地层，这会造成Cd、Hg和Pb等元素在成土母质中富集从而产生较高的背景值[56，57]. 因此推断研究区内Cd、Hg和Pb在土壤中的自然本底赋存较高.

2.3.3 低于重金属污染阈值的富硒土地

Se与重金属元素常有一定的协同作用[58]，为了合理利用清洁安全的富硒土地，结合FMDM模型提供的研究区土壤重金属的中污染和高污染临界值和重金属的空间分布，获得了低于土壤重金属中污染阈值和高污染阈值的富硒土地范围和面积（图 5和表 7）.

图 5 低于FMDM模型计算所得中污染阈值和高污染阈值的富硒土地优选分布 Fig. 5 Distribution of selenium-enriched land optimization according to the results of moderate and high pollution threshold calculated by FMDM model

表 7 低于不同污染阈值的富硒土地面积 Table 7 Selenium-rich land area below different pollution thresholds

研究区内低于10种重金属中污染阈值的富硒土地分布于城东街道和八陡镇的交界处，面积仅为2.22 km2，占研究区总面积的0.13%. 但是，基于本研究区内Cd、Hg和Pb的高背景值特征，以上3种重金属取高污染阈值，其他7种元素取中污染阈值的前提下，获得低于重金属高污染阈值富硒土地结果，主要分布在研究区西北的王村镇，以及研究区东部，呈东西方向延伸. 同时，在研究区其他部分也有零散分布，总面积达到291.81 km2，占研究区总面积的17.60%. 与不考虑重金属污染的富硒土地面积减少了约131.68 km2. 以上富硒土地可称之为“清洁富硒土地”，这对防止重金属进入人体，保护人类身体健康具有重要的作用与意义.

3 结论

（1）研究区土壤ω（Se）平均值为0.34 mg·kg-1，明显高于全国和山东省A层土壤背景值. 不同母质背景下，在砂岩中含量最高，其次为泥岩，在花岗岩中的含量最少；不同土壤类型背景下，在褐土中土壤Se含量最多，棕壤中含量最少. Se与10种重金属元素的相关系数表明，土壤Se主要来源于自然背景和人类活动.

（2）本文采用INLA-SPDE方法，整合土壤Se和环境协变量进行空间分布预测. 结果显示土壤Se高值区集中分布于研究区中部，含量受到自然因素和人类活动的影响. 空间预测的不确定性与采样点的布设密度有关. 富硒土地面积为423.49 km2，主要分布在研究区中部.

（3）FMDM模型将As、Co、Cr、Cu、Mn、Ni和Zn划分为背景值和中污染值双组分混合分布，且以上7种元素主要受到自然背景的影响；将Cd、Hg和Pb划分为背景值、中污染值和高污染值这3组分混合分布，来源不仅受到自然背景影响，而且受到人类活动的深刻干扰. 整合富硒土地分布和重金属空间分布，以FMDM所得污染临界值为标准，划分出低重金属污染的“清洁富硒土地”，低于中污染阈值的富硒土地面积和低于高污染阈值的富硒土地面积分别是2.22 km2和291.81 km2.

参考文献
[1] 何亚琳. 贵州省土壤含硒量及其分布[J]. 土壤学报, 1996, 33(4): 391-397.
He Y L. Se contents and distribution in soils of Guizhou Province[J]. Acta Pedologica Sinica, 1996, 33(4): 391-397.
[2] 赵成义, 任景华, 薛澄泽. 紫阳富硒区土壤中的硒[J]. 土壤学报, 1993, 30(3): 253-259.
Zhao C Y, Ren J H, Xue C Z. Selenium in soils of selenium-rich areas in Ziyang County[J]. Acta Pedologica Sinica, 1993, 30(3): 253-259.
[3] Lv Y Y, Yu T, Yang Z F, et al. Constraint on selenium bioavailability caused by its geochemical behavior in typical Kaschin-Beck disease areas in Aba, Sichuan Province of China[J]. Science of the Total Environment, 2014, 493: 737-749.
[4] Shi Z M, Pan P J, Feng Y W, et al. Environmental water chemistry and possible correlation with Kaschin-Beck Disease (KBD) in northwestern Sichuan, China[J]. Environment International, 2017, 99: 282-292.
[5] 范健, 任静华, 廖启林, 等. 苏南典型区农田土壤硒-镉拮抗作用研究[J]. 土壤, 2021, 53(5): 1023-1032.
Fan J, Ren J H, Liao Q L, et al. Antagonism between Se and Cd in typical farmland soil in southern Jiangsu province[J]. Soils, 2021, 53(5): 1023-1032.
[6] 张建东, 王丽, 赫栗涛, 等. 岚皋县岩石、土壤和农产品中硒分布规律研究[J]. 土壤通报, 2022, 53(1): 195-203.
Zhang J D, Wang L, He L T, et al. Distribution characteristics of selenium in rocks, soils and agricultural products of Langao[J]. Chinese Journal of Soil Science, 2022, 53(1): 195-203.
[7] 周菲, 彭琴, 王敏, 等. 土壤-植物体系中硒生物有效性评价研究进展[J]. 科学通报, 2022, 67(6): 461-472.
Zhou F, Peng Q, Wang M, et al. Advances in the evaluation of selenium bioavailability in soil-plant system[J]. Chinese Science Bulletin, 2022, 67(6): 461-472.
[8] 覃建勋, 付伟, 郑国东, 等. 广西岩溶区表层土壤硒元素分布特征与影响因素探究——以武鸣县为例[J]. 土壤学报, 2020, 57(5): 1299-1310.
Qin J X, Fu W, Zheng G D, et al. Selenium distribution in surface soil layer of karst area of Guangxi and its affecting factors: a case study of Wuming County[J]. Acta Pedologica Sinica, 2020, 57(5): 1299-1310.
[9] 周国华. 富硒土地资源研究进展与评价方法[J]. 岩矿测试, 2020, 39(3): 319-336.
Zhou G H. Research progress of selenium-enriched land resources and evaluation methods[J]. Rock and Mineral Analysis, 2020, 39(3): 319-336.
[10] 邵亚, 王毅伟, 蔡崇法, 等. 西南典型岩溶区土壤硒空间分布预测[J]. 农业工程学报, 2016, 32(22): 178-183.
Shao Y, Wang Y W, Cai C F, et al. Prediction on spatial distribution of soil selenium in typical karst area of Southwest China[J]. Transactions of the Chinese Society of Agricultural Engineering, 2016, 32(22): 178-183.
[11] 环境保护部, 国土资源部. 全国土壤污染状况调查公报[R]. 北京: 环境保护部, 2014. 1-3.
[12] Cai L M, Huang L C, Zhou Y Z, et al. Heavy metal concentrations of agricultural soils and vegetables from Dongguan, Guangdong[J]. Journal of Geographical Sciences, 2010, 20(1): 121-134.
[13] Hu W Y, Wang H F, Dong L R, et al. Source identification of heavy metals in peri-urban agricultural soils of Southeast China: an integrated approach[J]. Environmental Pollution, 2018, 237: 650-661.
[14] Jiang H H, Cai L M, Wen H H, et al. Characterizing pollution and source identification of heavy metals in soils using geochemical baseline and PMF approach[J]. Scientific Reports, 2020, 10(1). DOI:10.1038/s41598-020-63604-5
[15] Sun X F, Zhang L X, Lv J S. Spatial assessment models to evaluate human health risk associated to soil potentially toxic elements[J]. Environmental Pollution, 2021, 268. DOI:10.1016/j.envpol.2020.115699
[16] Wackernagel H. Multivariate geostatistics: an introduction with applications[M]. Berlin: Springer, 2003.
[17] Webster R, Oliver M A. Geostatistics for environmental scientists[M]. Chichester: John Wiley & Sons, 2007.
[18] 袁宏, 赵利, 王茂丽, 等. 拉萨河流域曲水一带土壤锗含量地统计学分析与评价[J]. 土壤通报, 2019, 50(5): 1079-1084.
Yuan H, Zhao L, Wang M L, et al. Geostatistical analysis and evaluation of soil germanium content in the Qushui area of Lhasa River Basin[J]. Chinese Journal of Soil Science, 2019, 50(5): 1079-1084.
[19] 高文武, 姜燕, 赵晋陵. 基于协同克里金插值法的土壤锰元素含量预测[J]. 地理与地理信息科学, 2018, 34(3): 119-124.
Gao W W, Jiang Y, Zhao J L. Predicting and mapping the Mn content in soil based on Cokriging[J]. Geography and Geo-Information Science, 2018, 34(3): 119-124.
[20] 陈琳, 任春颖, 王宗明, 等. 基于克里金插值的耕地表层土壤有机质空间预测[J]. 干旱区研究, 2017, 34(4): 798-805.
Chen L, Ren C Y, Wang Z M, et al. Prediction of spatial distribution of topsoil organic matter content in cultivated land using Kriging methods[J]. Arid Zone Research, 2017, 34(4): 798-805.
[21] Lv J S, Wang Y M. PMF receptor models and sequential Gaussian simulation to determine the quantitative sources and hazardous areas of potentially toxic elements in soils[J]. Geoderma, 2019, 353: 347-358.
[22] 吕建树. 烟台海岸带土壤重金属定量源解析及空间预测[J]. 地理学报, 2021, 76(3): 713-725.
Lyu J S. Source apportionment and spatial prediction of heavy metals in soils of Yantai coastal zone[J]. Acta Geographica Sinica, 2021, 76(3): 713-725.
[23] 徐辰星, 濮励杰, 朱明, 等. 基于序贯高斯条件模拟的土壤重金属含量预测与不确定性评价——以宜兴市土壤Hg为例[J]. 土壤学报, 2018, 55(4): 999-1006.
Xu C X, Pu L J, Zhu M, et al. Prediction of soil heavy metals content based on sequential gaussian simulation and evaluation of its uncertainties: a case study of soil Hg content in Yixing[J]. Acta Pedologica Sinica, 2018, 55(4): 999-1006.
[24] Lin Y P, Chang T K, Teng T P. Characterization of soil lead by comparing sequential Gaussian simulation, simulated annealing simulation and kriging methods[J]. Environmental Geology, 2001, 41(1-2): 189-199.
[25] 孙朝, 侯青叶, 杨忠芳, 等. 典型土壤环境中硒的迁移转化影响因素研究——以四川省成都经济区为例[J]. 中国地质, 2010, 37(6): 1760-1768.
Sun Z, Hou Q Y, Yang Z F, et al. Factors controlling the transport and transformation of selenium in typical soil environments: a case study of the Chengdu economic zone in Sichuan Province[J]. Geology in China, 2010, 37(6): 1760-1768.
[26] 张春来, 杨慧, 黄芬, 等. 广西马山县岩溶区土壤硒含量分布及影响因素研究[J]. 物探与化探, 2021, 45(6): 1497-1503.
Zhang C L, Yang H, Huang F, et al. Distribution and influencing factors of selenium content in soil in karst areas in Mashan County, Guangxi, China[J]. Geophysical and Geochemical Exploration, 2021, 45(6): 1497-1503.
[27] 陈东平, 张金鹏, 聂合飞, 等. 粤北山区连州市土壤硒含量分布特征及影响因素研究[J]. 环境科学学报, 2021, 41(7): 2838-2848.
Chen D P, Zhang J P, Nie H F, et al. Selenium distribution in soils of Lianzhou City, mountain area of northern Guangdong Province and its influencing factors[J]. Acta Scientiae Circumstantiae, 2021, 41(7): 2838-2848.
[28] Brus D J. Sampling for digital soil mapping: a tutorial supported by R scripts[J]. Geoderma, 2019, 338: 464-480.
[29] Walvoort D J J, Brus D J, De Gruijter J J. An R package for spatial coverage sampling and random sampling from compact geographical strata by k-means[J]. Computers & Geosciences, 2010, 36(10): 1261-1267.
[30] 朱钰, 郑屹然, 尹默. 统计学意义下的多重共线性检验方法[J]. 统计与决策, 2020, 36(7): 34-36.
Zhu Y, Zheng Y R, Yin M. Multicollinearity test under statistical significance[J]. Statistics and Decision, 2020, 36(7): 34-36.
[31] Molla A, Zuo S D, Zhang W W, et al. Optimal spatial sampling design for monitoring potentially toxic elements pollution on urban green space soil: a spatial simulated annealing and k-means integrated approach[J]. Science of the Total Environment, 2022, 802. DOI:10.1016/j.scitotenv.2021.149728
[32] Lindgren F, Rue H, Lindström J. An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach[J]. Journal of the Royal Statistical Society Series B: Statistical Methodology, 2011, 73(4): 423-498.
[33] Lindgren F, Rue H. Bayesian spatial modelling with R-INLA[J]. Journal of Statistical Software, 2015, 63(19): 1-25.
[34] Poggio L, Gimona A, Spezia L, et al. Bayesian spatial modelling of soil properties and their uncertainty: the example of soil organic matter in Scotland using R-INLA[J]. Geoderma, 2016, 277: 69-82.
[35] DD 2019-10, 天然富硒土地划定与标识(试行)[S].
[36] 李家熙, 黄怀曾, 刘晓端. 环境地球化学在农业和生命科学上的应用研究[J]. 第四纪研究, 1999, 19(3): 224-230.
Li J X, Huang H Z, Liu X D. The application of environmental geochemistry to agriculture and life science[J]. Quaternary Sciences, 1999, 19(3): 224-230.
[37] 严明书, 龚媛媛, 杨乐超, 等. 重庆土壤硒的地球化学特征及经济意义[J]. 物探与化探, 2014, 38(2): 325-330.
Yan M S, Gong Y Y, Yang L C, et al. Geochemical characteristics and economic significance of the Se-Rich soil in ChongQing[J]. Geophysical and Geochemical Exploration, 2014, 38(2): 325-330.
[38] 郦逸根, 董岩翔, 郑洁, 等. 浙江富硒土壤资源调查与评价[J]. 第四纪研究, 2005, 25(3): 323-330.
Li Y G, Dong Y X, Zheng J, et al. Selenium: abundant soil survey and assessment in Zhejiang[J]. Quaternary Sciences, 2005, 25(3): 323-330.
[39] 中国土壤元素背景值[M]. 北京: 中国环境科学出版社, 1990.
[40] Wang Z, Chen X M, Yu D Q, et al. Source apportionment and spatial distribution of potentially toxic elements in soils: a new exploration on receptor and geostatistical models[J]. Science of the Total Environment, 2021, 759. DOI:10.1016/j.scitotenv.2020.143428
[41] Wang S, Cai L M, Wen H H, et al. Spatial distribution and source apportionment of heavy metals in soil from a typical county-level city of Guangdong Province, China[J]. Science of the Total Environment, 2019, 655: 92-101.
[42] Nacke H, Gonçalves A C, Schwantes D, et al. Availability of heavy metals (Cd, Pb, And Cr) in agriculture from commercial fertilizers[J]. Archives of Environmental Contamination and Toxicology, 2013, 64(4): 537-544.
[43] Franco-Uría A, López-Mateo C, Roca E, et al. Source identification of heavy metals in pastureland by multivariate analysis in NW Spain[J]. Journal of Hazardous Materials, 2009, 165(1-3): 1008-1015.
[44] Song H Y, Hu K L, An Y, et al. Spatial distribution and source apportionment of the heavy metals in the agricultural soil in a regional scale[J]. Journal of Soils and Sediments, 2018, 18(3): 852-862.
[45] Lee D H, Kim J H, Mendoza J A, et al. Characterization and source identification of pollutants in runoff from a mixed land use watershed using ordination analyses[J]. Environmental Science and Pollution Research, 2016, 23(10): 9774-9790.
[46] Jin Z, Lv J S. Integrated receptor models and multivariate geostatistical simulation for source apportionment of potentially toxic elements in soils[J]. CATENA, 2020, 194. DOI:10.1016/j.catena.2020.104638
[47] Herndon E M, Jin L X, Brantley S L. Soils reveal widespread manganese enrichment from industrial inputs[J]. Environmental Science & Technology, 2011, 45(1): 241-247.
[48] 王若锦, 邵天杰, 卫佩茹. 环青海湖地区表层土壤重金属富集含量及其生态风险评价[J]. 干旱区研究, 2021, 38(2): 411-420.
Wang R J, Shao T J, Wei P R. Enrichment content and ecological risk assessment of heavy metal in surface soil around Qinghai Lake[J]. Arid Zone Research, 2021, 38(2): 411-420.
[49] 孙雪菲, 张丽霞, 董玉龙, 等. 典型石化工业城市土壤重金属源解析及空间分布模拟[J]. 环境科学, 2021, 42(3): 1093-1104.
Sun X F, Zhang L X, Dong Y L, et al. Source apportionment and spatial distribution simulation of heavy metals in a typical petrochemical industrial city[J]. Environmental Science, 2021, 42(3): 1093-1104.
[50] 杨泽, 刘国栋, 戴慧敏, 等. 黑龙江兴凯湖平原土壤硒地球化学特征及富硒土地开发潜力[J]. 地质通报, 2021, 40(10): 1773-1782.
Yang Z, Liu G D, Dai H M, et al. Selenium geochemistry of soil and development potential of Se-rich soil in Xingkai Lake Plain[J]. Geological Bulletin of China, 2021, 40(10): 1773-1782.
[51] 邵龙义, 徐小涛, 王帅, 等. 中国含煤岩系古地理及古环境演化研究进展[J]. 古地理学报, 2021, 23(1): 19-38.
Shao L Y, Xu X T, Wang S, et al. Research progress of palaeogeography and palaeoenvironmental evolution of coal-bearing series in China[J]. Journal of Palaeogeography, 2021, 23(1): 19-38.
[52] 常嘉, 陈世悦, 鄢继华. 淄博博山地区晚古生代煤系层序地层与聚煤作用[J]. 沉积学报, 2019, 37(5): 968-980.
Chang J, Chen S Y, Yan J H. Sequence stratigraphy and coal accumulation in Late Paleozoic coal-bearing strata in Zibo Boshan area[J]. Acta Sedimentologica Sinica, 2019, 37(5): 968-980.
[53] Lv J S, Liu Y. An integrated approach to identify quantitative sources and hazardous areas of heavy metals in soils[J]. Science of the Total Environment, 2019, 646: 19-28.
[54] Lin Y P, Cheng B Y, Shyu G S, et al. Combining a finite mixture distribution model with indicator kriging to delineate and map the spatial patterns of soil heavy metal pollution in Chunghua County, central Taiwan[J]. Environmental Pollution, 2010, 158(1): 235-244.
[55] 庞绪贵, 代杰瑞, 陈磊, 等. 山东省17市土壤地球化学背景值[J]. 山东国土资源, 2019, 35(1): 46-56.
Pang X G, Dai J R, Chen L, et al. Soil geochemical background value of 17 cities in Shandong province[J]. Shandong Land and Resources, 2019, 35(1): 46-56.
[56] 王磊, 唐文春, 秦兵, 等. 四川龙门山地区磷矿、煤矿开采对水系沉积物Cd等元素影响调查[J]. 地质科技情报, 2007, 26(6): 36-41.
Wang L, Tang W C, Qin B, et al. The survey of Cd and other elements in river sediments affected by phosphorite deposit and coal mine in the area of longmenshan mountain, Sichuan Province[J]. Geological Science and Technology Information, 2007, 26(6): 36-41.
[57] 张蒜, 陈健, 张博斐, 等. 安徽淮南煤田煤系地层铝质泥岩的地球化学特征[J]. 中国煤炭地质, 2020, 32(8): 20-25, 47.
Zhang S, Chen J, Zhang B F, et al. Geochemical features of coal measures strata aluminous mudstone in Huainan coalfield, Anhui[J]. Coal Geology of China, 2020, 32(8): 20-25, 47.
[58] 冯辉, 张学君, 张群, 等. 北京大清河流域生态涵养区富硒土壤资源分布特征和来源解析[J]. 岩矿测试, 2019, 38(6): 693-704.
Feng H, Zhang X J, Zhang Q, et al. Distribution characteristics and sources identification of selenium-rich soils in the ecological conservation area of the daqinghe river watershed, Beijing[J]. Rock and Mineral Analysis, 2019, 38(6): 693-704.