环境科学  2023, Vol. 44 Issue (4): 2177-2191   PDF    
基于地理探测器和多源数据的耕地土壤重金属来源驱动因子及其交互作用识别
张宏泽1,2, 崔文刚1,2, 刘绥华1,2, 崔瀚文1,2, 黄月美1,2     
1. 贵州师范大学地理与环境科学学院, 贵阳 550025;
2. 贵州省山地资源与环境遥感应用重点实验室, 贵阳 550001
摘要: 明确耕地土壤重金属来源对土壤健康合理管护及其可持续发展至关重要.利用正定矩阵因子分解(PMF)模型源解析结果(源成分谱和源贡献)、历史调查数据和时序遥感数据,整合地理探测器(GD)、最佳参数地理探测器(OPGD)、空间关联探测器(SPADE)和空间关联交互探测器(IDSA)模型,探讨了土壤重金属来源空间分层异质性研究中可塑性面积单元问题(MAUP),分别识别了分类变量和连续变量中控制土壤重金属来源空间分异的驱动因子及其交互作用.结果表明,中小尺度下土壤重金属来源的空间分异性研究受空间尺度影响,0.08 km2能够作为区域土壤重金属来源空间分异性探测的可选空间单元.连续变量受空间相关性与离散化水平影响,分位数法和中断数为10的离散化参数组合能够减少连续变量在土壤重金属来源空间分异性探测中所受分区效应的影响.分类变量中地层(PD 0.12~0.48)控制土壤重金属来源空间分异,地层和流域的交互作用对各来源解释力为27.28%~60.61%,各来源的高风险区域分别分布在地层中的下震旦系、上白垩系、土地利用中采矿用地和土壤种类中典型强淋溶土区域.连续变量中人口(PSD 0.40~0.82)控制土壤重金属来源空间分异,连续变量的空间组合对各来源解释力为61.77%~78.46%,各来源的高风险区域分别分布在蒸腾量(41.2~43 kg ·m-2)、距河流距离(315~398 m)、增强型植被指数(0.796~0.995)和距河流距离(499~605 m)区域.研究结果为耕地土壤重金属来源驱动因子及其交互作用的相关研究提供参考,为喀斯特地区耕地土壤健康合理管护及其可持续发展提供重要的科学依据.
关键词: 耕地土壤      重金属      多源数据      地理探测器      驱动因子      可塑性面积单元问题     
Identifying Driving Factors and Their Interacting Effects on Sources of Heavy Metal in Farmland Soils with Geodetector and Multi-source Data
ZHANG Hong-ze1,2 , CUI Wen-gang1,2 , LIU Sui-hua1,2 , CUI Han-wen1,2 , HUANG Yue-mei1,2     
1. School of Geography and Environmental Science, Guizhou Normal University, Guiyang 550025, China;
2. Guizhou Mountain Resources and Environmental Remote Sensing Application Laboratory, Guiyang 550001, China
Abstract: The identification of heavy metal sources in farmland soils is essential for the rational health condition management and sustainable development of soil. Using source resolution results(source component spectrum and source contribution)of a positive matrix factorization(PMF)model, historical survey data, and time-series remote sensing data, integrating a geodetector(GD), an optimal parameters-based geographical detector(OPGD), a spatial association detector(SPADE), and an interactive detector for spatial associations(IDSA)model, this study explored the modifiable areal unit problem(MAUP) of spatial heterogeneity of soil heavy metal sources and identified the driving factors and their interacting effects on the spatial heterogeneity of soil heavy metal sources in categorical and continuous variables, respectively. The results showed that the spatial heterogeneity of soil heavy metal sources at small and medium scales was affected by the spatial scale, and the optional spatial unit was 0.08 km2 for detecting spatial heterogeneity of soil heavy metal sources in the study region. Considering spatial correlation and discretization level, the combination of the quantile method and discretization parameters with an interruption number of 10 could be implied to reduce the partitioning effects on continuous variables in the detection of spatial heterogeneity of soil heavy metal sources. Within categorical variables, strata(PD 0.12-0.48) controlled the spatial heterogeneity of soil heavy metal sources, the interaction between strata and watersheds explained 27.28%-60.61% of each source, and the high-risk areas of each source were distributed in the lower sinian system, upper cretaceous in strata, mining land in land use, and haplic acrisols in soil types. Within continuous variables, population (PSD 0.40-0.82) controlled the spatial variation in soil heavy metal sources, and the explanatory power of spatial combinations of continuous variables for each source ranged from 61.77% to 78.46%. The high-risk areas of each source were distributed in evapotranspiration (41.2-43 kg·m-2), distance from the river (315-398 m), enhanced vegetation index (0.796-0.995), and distance from the river (499-605 m). The results of this study provide a reference for the research of the drivers of heavy metal sources and their interactions in arable soils and provide an important scientific basis for the management of arable soil and its sustainable development in karst areas.
Key words: farmland soil      heavy metals      multi-source data      geodetector      driving factor      modifiable areal unit problem     

土壤是人类生存和发展的重要自然资源, 其在联合国可持续发展目标中得到前所未有的重视[1].然而, 存在于地球关键带上的表层土壤易受到不同来源和程度的重金属污染[2~5].土壤重金属污染是导致土壤环境遭到破坏的主要原因, 现已成为全球范围内的一个严重问题, 且粮食安全与土壤环境质量息息相关[4~7].现阶段土壤环境质量不容乐观, 高强度和超负荷利用致使耕地土壤健康状况堪忧[5, 8].喀斯特地区独特的水力和水文地质特征使表层土壤极易受到人类活动的影响[9, 10].因此, 对控制喀斯特地区耕地土壤重金属来源的驱动因子及其交互作用的识别研究具有明确的科学与现实意义.

现阶段, 对土壤重金属源解析的研究主要分为定性识别和定量解析两类.定性识别土壤中重金属来源的方法主要以因子分析、主成分分析、聚类分析、变量逻辑回归和空间分析等为主[11, 12]; 定量解析各来源对各土壤重金属采样点位含量贡献的方法多采用源排放清单、化学质量平衡、UN-MIX模型、绝对主成分得分-多元线性回归、正定矩阵因子分解(positive matrix factorization, PMF)、有限混合分布、条件树推断和随机森林等[11, 12]; 此外, 越来越多的学者结合地理信息系统和多元统计分析技术揭示区域土壤重金属空间特征, 识别其来源[6].空间相关性[13, 14]作为空间数据最重要的特征, 在空间数据层内通过Moran's I、Geary's C、Getis-Ord General G、Ripley's K、半变异函数和空间计量模型等描述, 在层间以地理加权回归描述[14, 15].有研究者基于上述方法或模型进行土壤重金属源解析的适应性研究, 合理分析出人为或自然因素的单方面影响[11, 12].然而, 土壤重金属受到自然和人类活动的双重影响, 各种环境因子的多样性、复杂性、非线性和异质性共同影响土壤重金属累积的整个过程[16].传统统计模型在计算交互作用时需假设线性和受多变量共线影响, 无法准确量化土壤重金属各来源之间的交互作用[10, 17].

受上述问题困扰, 许多学者采用地理探测器(geodetector, GD)模型进行土壤重金属空间分异性研究[18~23].GD模型基于空间分层异质性理论探索地理因素[24, 25], 其广泛应用于气候与环境、土地利用、污染与人类健康、水资源与动力学等方面决定因素的探测中[26].然而, 目前基于GD模型进行空间分异性研究时易忽略地理变量的空间相关性和可塑性面积单元问题(modifiable areal unit problem, MAUP), 其是地理学中长期存在和影响深远的问题之一, 分为尺度效应和分区效应[27~31].空间数据受空间相关性和离散化水平影响, 在应用时需进行适应性分析, 已有研究表明统计结果因规模和聚合而不同[15, 27, 31].伴随GD模型的应用, 基于最佳参数的空间数据离散化、空间尺度效应、空间相关性和基于空间相关性的交互作用探测被整合到GD模型中[17].现阶段, 缺乏集合空间尺度、空间相关性和离散化水平进行土壤重金属来源空间分异性的适应性研究.采用多源数据, 尤其是引入时序遥感数据进行土壤重金属空间分异性的研究相对较少.

PMF模型在求解过程中对因子载荷和因子得分均做非负约束, 所提取因子的空间分布与其来源的空间分布相似[23], 源贡献具有可解释性和明确的物理意义[32].本研究利用地理探测器(整合空间尺度、空间相关性、空间离散化)和多源数据识别土壤重金属来源与环境因子之间的空间关系, 揭示多种环境因子如何共同影响土壤重金属来源的累积, 通过指出各重金属来源进入土壤的途径, 以期为耕地土壤健康合理管护及其可持续发展提供科学依据.

1 材料与方法 1.1 数据来源和处理

本研究以贵阳市乌当区北部区域为研究区, 以378个耕地土壤样品中Cd、Hg、As、Pb和Cr的含量数据为基础数据.研究区概况、样品采集及预处理、污染评价方法和PMF模型解析结果见文献[33].本研究在前期研究基础上, 分别统计整理: ①378个采样点处酸碱度(pH); ②各采样点综合污染负荷指数(PLI); ③PMF解析4个土壤重金属来源(F1~F4)对各采样点处含量贡献.对以上数据采用普通克里金法进行空间插值, 各土壤重金属来源含量贡献和综合污染负荷指数等值区域图见图 1.

图 1 各土壤重金属来源含量贡献和综合污染负荷指数等值区域 Fig. 1 Equivalent regional map of the content contribution of each soil heavy metal source and the integrated pollution load index

1.1.1 环境因子选取

环境因子的选取对探测土壤中重金属来源的空间分异性至关重要[6].以土壤成土因素和前人关于土壤重金属来源的研究成果为依据[10, 34], 选择土壤性质、地形、生态、社会经济和距离这5种类型作为解释变量.

距离因子空间分布见图 2, 土壤性质、地形、生态和社会经济因子空间分布见图 3.

图 2 距离因子空间分布示意 Fig. 2 Distance factors spatial distribution map

图 3 解释变量 Fig. 3 Explanatory variables

1.1.2 环境因子来源

(1) 土壤性质

土壤中的重金属受成土母质控制, 不同发育类型土壤中重金属含量差异大[6, 10].土壤质地和土壤种类数据来源于国家青藏高原科学数据中心的基于世界土壤数据库(harmonized world soil database, HWSD)的中国土壤数据集[35].地层由中国地质调查局(https://www.cgs.gov.cn)公布的1∶200 000地质图矢量化处理所得.

(2) 地形因子

地形是区域和局部空间上地球动力学过程的重要驱动因素[36].数字高程模型(DEM)来源于航天飞机雷达地形任务数据(https://asf.alaska.edu), 分辨率为30 m.坡度、坡向、平面曲率、剖面曲率、标准曲率和流域由DEM经自动化地学分析系统(system for automated geoscientific analyses, SAGA)处理所得[36].

(3) 生态因子

遥感数据常作为评估区域生态环境和社会经济的有效数据源, 谷歌云计算平台(google earth engine, GEE)常用于快速获取并处理遥感数据[34, 37~39].归一化植被指数(NDVI)和增强型植被指数(EVI)来源于MOD13Q1数据集, 分别通过GEE收集23幅影像; 叶面积指数(LAI)基于MOD15A2H数据集, 共收集46幅影像; 蒸腾量(ET)基于MOD16A2数据集, 共收集45幅影像; 年度净初级生产量(NPP)基于MOD17A3HGF数据集, 共收集1幅影像; 遥感生态指数(RSEI)[39]基于USGS Landsat8 Level2 Collection2 Tier1数据集, 由云和水掩膜后的影像利用主成分分析技术并集成绿度、湿度、热度和干度所计算的区域生态环境质量.其中, 基于空间相关性的时间不变效应来降低数据源获取时间不同所带来的影响, 所收集遥感数据的时间分辨率选择在2019年; 为突出研究区生态状况和抑制云层、大气和太阳高度等方面所导致的噪声, NDVI、EVI、LAI和ET经最大合成法进行年度数据的制作.

(4) 社会经济因子

人类的生产和生活等经营性活动影响土壤中微量元素的积累[6, 10].土地利用数据来源于贵阳市乌当区2018年土地利用现状数据; 人口数据来源于WorldPop(https://worldpop.org)的开放空间人口统计数据集, 其以100 m×100 m网格评估2019年人口规模; 平均夜间灯光指数(VIIRS)基于VIIRS Nighttime Day/Night Band Composites Version1数据集, 是在GEE平台经时间序列合成的地面分辨率500 m的2019年度影像1幅.

(5) 距离因子

交通和工农业活动过程中产生的废气、废水和废渣会对周边一定范围内的土壤重金属含量产生影响[40].道路、铁路、河流和开阳磷矿位置以ArcGIS世界影像图为底图进行绘制, 基于ArcGIS10.2计算各采样点距道路、铁路、河流和开阳磷矿的欧式距离.

1.1.3 空间离散化

(1) 分类变量

土壤质地和土壤种类基于HWSD数据集进行分类; 土地利用依据土地利用现状分类标准进行划分; 地层按成层岩石所形成的年代进行划分; 流域按可能存在的集水区域进行划分.

(2) 基于专业经验的连续变量空间离散化

基于历史文献和相关标准将解释变量中的连续变量离散化为分类变量.pH和坡向基于专业经验进行划分; 坡度基于《第三次全国国土调查技术规程》(TD/T1055-2019)划分; 遥感生态指数参考2015年修订的《生态环境状况评价技术规范》(HJ/T192-2006)划分; 平面曲率、剖面曲率和标准曲率按加速(分散)为正, 减速(汇聚)为负, 线性为零划分[41].

(3) 基于最佳参数的连续变量空间离散化

数字高程模型、归一化植被指数、增强型植被指数、叶面积指数、蒸腾量、年度净初级生产量、人口、夜间灯光指数、距道路距离、距铁路距离、距河流距离和距开阳磷矿距离, 分别采用OPGD和SPADE模型进行最佳参数空间离散化.

1.2 研究方法 1.2.1 地理探测器(GD)模型

GD模型基于地理现象的空间分异定律探测空间分异性, 并揭示其背后驱动力的一组统计学方法, 分为因子探测、交互作用探测、风险探测和生态探测, 具体见文献[24, 25].

(1) 因子探测  GD模型的核心, 探测响应变量的空间分异性, 通过PD揭示变量的相对重要性, 表达式为:

(1)

式中, h=1, …, L为解释变量的分层, 即分类或分区; NhN分别为分层h和全区单元数; σh2σ2分别是分层h和全区的响应变量方差.PD范围为[0, 1], 以非中心F分布检验显著性.

(2) 交互作用探测  交互区域为两个解释变量空间上交集区域, 基于因子探测计算交互区域PD, 将其定义为PID.通过比较两个单一解释变量PD和交互作用PID探索交互作用.交互作用类型见文献[25].

(3) 风险探测  探测解释变量各分层的属性均值之间是否存在显著差异, 以t统计量检验显著性.

(4) 生态探测  探测解释变量对响应变量空间分布的影响是否存在显著差异, 以F统计量检验显著性.

1.2.2 基于最佳参数的地理检测器(OPGD)模型

OPGD模型是对GD模型的改进, 原理及计算公式一致, 改进在于连续数据的空间离散化.OPGD模型通过不同离散化方法(等间段、自然间段、分位数间段、几何间段和标准差间段)和中断数(3~6)的参数组合寻求连续变量的最大PD.对于MAUP效应, OPGD模型通过比较各空间单元PD第90%分位数大小来评估尺度效应[26].

1.2.3 空间自相关

Tobler's地理学第一定律指出: 所有的地理事物都是相互联系的, 但离得越近的事物彼此之间的联系越强[13].作为地理现象表现形式的连续数据具有空间相关性, 通过全局Moran's I统计量测量其要素位置和属性值之间的空间自相关性, 以Z得分和P值进行检验.原理和计算公式见文献[42].

1.2.4 空间关联检测器(SPADE)模型

SPADE模型考虑变量间空间联系和最小化离散化水平的影响, 其基于解释变量空间离散化水平对响应变量进行分区, 以空间加权叉积[43]计算的空间方差代替传统的方差来衡量彼此之间空间相关性.此外, 通过多级离散化和考虑离散化导致的信息损失来减少离散化水平对PSD影响[15].表达式为:

(2)

式中, Nh(h=1, …, L)为分层h单元数, N为全区单元数; ΓhΓ分别为分层h和全区的空间方差, Γ定义为空间加权叉积的平均值, 表达式为:

(3)

式中, Wij为位置i和位置j之间空间加权矩阵的权重值; Cij为位置ij之间响应变量的相似性, 以半平方差函数(xi-xj)2/2表示.

对于离散化会带来信息损失, 以响应变量和解释变量PSD之比来调整偏差, 调整后表达式为:

(4)
1.2.5 基于空间关联的交互检测器(IDSA)模型

IDSA模型用于检验响应变量和解释变量之间空间关联[17].IDSA模型通过最少数量的解释变量、空间重叠区和模糊划分区来解释响应变量, 步骤如下:

(1) 空间模糊叠加  空间模糊叠加是解释变量在空间上模糊关系的重叠过程.模糊关系假设解释变量之间有一定联系, 以模糊隶属度函数(以模糊数表示)量化.模糊数等于GD模型风险检测器中归一化后的平均风险值, 表达式为:

(5)

式中, Xi(i=1, …, m)为解释变量; fn(Xi)为模糊数; η(Xi)为平均风险值; g(X)为归一化函数.

(2) 具有交互作用的组合模糊数  基于模糊AND算子对模糊数进行积分, 估计解释变量间交互作用的组合模糊数, 表达式为:

(6)

式中, X1X2∩…∩Xm为解释变量之间的交互作用.

(3) 最佳交互作用的计算和偏差调整  交互作用的计算为具有交互作用的单个解释变量空间叠加区域的方差之和与单个解释变量的全局方差之比, 表达式为:

(7)

式中, Ni, h(i=1, …, m; h=1, …, Li)为单个解释变量Xi在分区h上与其他解释变量具有交互作用的单元数; Γi, h为解释变量Xi在分区h上空间方差; NiΓi分别为解释变量Xi单元数和全区空间方差.

对于离散化过程中信息缺失而导致的偏差, 以响应变量和解释变量PSD之比进行调整, 以t统计量检验显著性, 调整后表达式为:

(8)
1.3 采样

以研究区外接矩形为范围进行网格与样点制作.正六边形具有邻近空间概念清晰、视野偏差少和信息空间结构明确等特点[44].因此, 本研究选用六边形格网来减少边界效应触发的样本偏差, 共设计14组空间单元探究MAUP效应, 各空间单元面积和样点数见表 1.

表 1 各空间单元面积和样点数描述性统计 Table 1 Descriptive statistics of space unit area and sample number

1.4 数据处理

使用Excel 2019进行数据整理和计算; 利用GEE进行遥感数据收集和处理; 地形因子计算由SAGA完成; 采用ArcGIS10.2进行栅格数据处理、空间相关性检测、空间单元布设和采样; 基于R语言GD包、SPADE包和IDSA包进行数据分析与绘图.

2 结果与分析 2.1 最佳空间单元的确定

在空间分异性探讨方面, 确定合理空间单元的常用方法是基于OPGD模型比较各空间单元PD第90%分位数(PD90%)的大小[26].以PLI作为响应变量探究MAUP效应, 各空间单元尺度效应比较见图 4.由单变量PD变化趋势和PD90%局部加权回归(LOESS)模型分析结果可知, 随着空间单元面积增大, PD排序稳定, PD90%逐渐增加; 具有显著性(P < 0.05)的解释变量个数逐渐减少.单变量PD变化表明, 随着空间单元增大, 分类变量PD高而稳定, 在空间单元增大到0.08 km2后, 连续变量PD变化幅度大且一些连续变量未通过显著性检验.PD90%的LOESS模型表明, PD90%增长速率趋势为增长-下降-增长, 当空间单元从0.01 km2增大到0.16 km2时, 增长速率先增后降, 增长速率在0.08 km2时达到局部最大; 当空间单元从0.16 km2增大到1.28 km2时, 增长速率提升较快, PD90%变化波动强烈.从整体空间格局来看, 0.16 km2的空间单元具有相对较高的PD90%.空间单元从0.01 km2增大至0.08 km2时, PD90%增长速率稳步增大, PD90%在0.08 km2时达到局部最优解且能够去除无效因子.因此, 从局部格局的稳定性来看, 0.08 km2是本研究用于评估环境因子对土壤重金属源空间分布影响的最佳空间单元.

显著性P < 0.05 图 4 基于OPGD模型的各空间单元尺度效应比较 Fig. 4 OPGD-based comparison of size effects of spatial units

2.2 最佳空间单元下土壤重金属来源的驱动因子及其交互作用识别

在最佳空间单元下, 以F1~F4作为响应变量, 各环境因子作为解释变量, 从两方面探究环境因子中控制土壤重金属来源的驱动因子及其之间的交互作用: ①基于GD模型对分类变量进行识别; ②基于OPGD模型、SPADE模型和IDSA模型识别连续变量.

2.2.1 基于GD模型的分类变量驱动因子及其交互作用识别

由GD模型分析结果可知(图 5), 地层是分类变量中导致各源空间分异的主导因子.F1主要驱动因子分别为地层(45.78%)、pH(10.09%)和流域(8.06%), 地层为下震旦系的区域具有最高的源贡献率风险, 与其他类型的源贡献率均值之间差异最大, 下三叠系则为最低风险区域; 交互探测表明地层和其他解释变量的交互作用表现为增强, 其中地层∩流域对F1贡献率为58.57%, 两者呈非线性增强.F2主要驱动因子分别为地层(30.41%)、流域(14.54%)和坡度(5.16%), 地层为上白垩系的区域具有高风险, 地层∩流域呈非线性增强, 其对F2贡献率为49.28%.F3主要驱动因子分别为地层(11.71%)、pH(5.14%)和流域(3.39%), 土地利用为采矿用地的区域具有高风险, 地层∩流域呈非线性增强, 其对F3贡献率为27.28%.F4主要驱动因子分别为地层(48.47%)、流域(16%)和坡度(7.69%), 土壤种类为典型强淋溶土的区域具有高风险, 地层∩流域呈线性增强, 其对F4贡献率为60.61%.生态矩阵探测结果表明地层对各源贡献率的影响分别与土壤种类和土壤质地显著不同.

图 5 基于GD模型的分类变量对土壤重金属来源解释力探测 Fig. 5 GD-based categorical variables for soil heavy metal source explanatory power detection

2.2.2 基于OPGD模型、SPADE模型和IDSA模型的连续变量驱动因子及其交互作用识别

(1) 连续变量的空间相关性检测  SPADE模型要求连续变量具有空间自相关假设[15].基于反距离平方策略, 通过ArcGIS10.2进行连续变量全局Moran's I统计量计算.由表 2可知, 各连续变量全局Moran's I指数大于0.2, Z得分大于1.96, P值小于0.01, 具有统计学意义, 结果表明各连续变量空间分布呈显著聚集特征.

表 2 各连续变量全局Moran's I分析结果1) Table 2 Results of Global Moran's I analysis for continuous variables

(2) 基于SPADE模型的连续变量最佳参数离散化  由图 6可知, 大部分连续变量的PSD随中断数增加而增大, 连续变量对各源解释力的变化趋势一致.随中断数增加, PSD75%逐步增大, 增长速率总体上呈由高递减趋势.当中段数为7~18时, 各连续变量PSD排序稳定, 其中, 当中断数为10时, F1和F3的PSD75%增长率降到5%; F2和F4的PSD75%增长率先增(微变)后降, 最大增长率未达5%, 但在中断数9~10时增长率达到最大值.结合PSD变化趋势和PSD75%增长率确定各连续变量的最佳中断数为10.

图 6 基于SPADE模型的连续变量最佳空间离散化 Fig. 6 SPADE-based optimal spatial discretization of continuous variables

(3) 基于OPGD模型和SPADE模型的连续变量单因子作用和风险区探测  由图 7(a)可知, 人口、距开阳磷矿距离和距铁路距离的PSD较高, 表明其与各来源的空间关联度高, 而EVI、距道路距离和距河流距离与各源空间关联低.除F2中距道路距离的PSD较PD下降0.074 9外, 其余连续变量PSD高于PD, 大部分连续变量空间关联性增强.人口是连续变量中控制研究区土壤重金属源空间分布的主导因子.由图 7(b)可知, 连续变量划分后的各分层内各源贡献率不同.其中, F1、F2、F3和F4贡献率较高的区域分别分布在ET为41.2~43 kg·m-2、距河流距离为315~398 m、EVI为0.796~0.995和距河流距离为499~605 m的区域.

(a)基于OPGD和SPADE模型的单因子作用, (b)基于SPADE模型的风险区探测, (c)基于OPGD模型的交互作用探测, (d)基于IDSA模型的交互作用探测, (e)IDSA模型中叠加区内观测数密度分布, (f)IDSA模型中具有显著差异的叠加区对概率分布 图 7 基于OPGD、SPADE和IDSA模型的连续变量对土壤重金属来源解释力探测 Fig. 7 Continuous variables based on OPGD, SPADE and IDSA models for soil heavy metal source explanatory power detection

(4) 基于OPGD模型和IDSA模型的交互作用探测和模型性能评估  由表 3图 7(c)可知, 对于F1、F2和F4, 连续变量(DEM、LAI、NDVI、VIIRS、距道路距离和距开阳磷矿距离彼此之间存在交互作用)的PID分别为0.148 3、0.222 2和0.201 7; 对F3, 连续变量(DEM、ET、NDVI、NPP、距道路距离和距铁路距离彼此之间存在交互作用)的PID为0.063 6.由表 3可知, IDSA模型交互探测中12个连续变量交互作用的PIDIDSA分别为0.646 5、0.784 6、0.617 7和0.776 1, 图 7(d)表明, VIIRS∩ET∩NDVI∩NPP的空间交互对F1的解释力为57.89%; 距开阳磷矿距离∩VIIRS∩距河流距离的空间交互对F2的解释力为66.01%; 人口∩EVI∩ET∩NPP的空间交互对F3的解释力为60.99%; VIIRS∩距开阳磷矿距离∩DEM∩LAI的空间交互对F4的解释力为76.41%.

表 3 连续变量最佳交互作用探测和模型性能评估 Table 3 Optimal interaction detection and model performance evaluation for continuous variables

在模型性能评估中, OPGD模型划分的交互作用区域为262~423个, 叠加区对数量为34 191~89 253对, 单观测值分区出现40~94个; IDSA模型划分的交互作用区域为52~68个, 叠加区对数量为1 326~2 278对, 单观测值分区出现5~9个.此外, IDSA模型中连续变量的密度曲线[图 7(e)]近似钟形, 大部分值聚集在中间区域, 叠加区对之间存在显著差异的概率[图 7(f)]在18.72%~41.70%之间.

3 讨论

本研究基于空间分异性原理, 通过环境因子描述和解释PMF模型所解析出各来源的空间分异性, 揭示对各源起主导作用的环境因子及其之间的交互作用.先前土壤重金属来源空间分异性的研究大多是基于GD或OPGD模型来解释变量之间关联[18~23, 45, 46].然而, GD模型基于数据特征和OPGD模型基于最大PD原则所进行的空间数据离散化仅考虑连续变量分布, 忽略解释变量与响应变量之间相关性[15, 17].此外, MAUP存在于大规模和空间数据明确的任何方面, 研究者在评估地理现象时所选空间单元的面积和方法会对统计分析产生影响[27, 30, 31, 47].MAUP的分区效应体现在连续变量的离散化水平上, 太少和太多的离散化水平分别描述原始数据特征不足和产生过多冗余信息的情况, 阻碍分析和解释[48]. 有研究表明, 环境因子空间尺度的选取在探讨土壤重金属累积方面非常重要[10], 遥感数据分辨率与MAUP的尺度效应相对应, 不合理的空间单元难以准确响应自然系统的复杂性[27].也有研究表明, 岩溶土壤中微量元素累积的机制复杂, 研究区域规模和边界的具体选择会对导致土壤中微量元素累积主导因子的确定产生影响, 甚至可以极大地改变研究结果[10].由于喀斯特地区土壤的特异性, 确定最佳空间单元仍然具有难度.另一方面, 总是很难确定最佳离散化组合参数, 因为最佳离散化水平与数据特征和具体研究问题有关[48].研究结果表明, 相对于GD和OPGD模型, SPADE和IDSA模型明确考虑空间相关性和补偿离散化导致的信息损失, 两者联合使得土壤重金属来源空间分异性的探测结果更稳定和准确.因此, 本研究设计14组数据探究土壤重金属来源分异性研究中的MAUP效应具有明确的科学意义.

已有研究多从历史调查数据中提取环境因子[10, 16, 49, 50], 结合时序遥感数据进行分析的研究较少.本研究基于土壤成土因素和前人关于土壤重金属源解析的研究成果选取环境因子, 采用多源数据探测土壤重金属来源的空间分异性, 以期补充遥感数据在土壤重金属来源空间分异性研究中的应用.但本研究仍有不足.首先, 能够描述土壤性质的因子较少; 其次, 缺乏农业经营性活动和源清单数据.

对影响土壤重金属来源的主导因子及其交互作用的识别有助于区分重金属的来源和评估其风险[6, 10, 22].研究结果表明, 地层、流域、坡度、pH、人口、距开阳磷矿距离和距铁路距离是研究区内导致各来源贡献空间分异的主要因子.有研究表明, 喀斯特地区复杂的地质背景导致土壤重金属含量变异大, 碳酸盐岩相关地层是最易导致土壤重金属含量高的因素之一, 碎屑岩类岩石在氧化环境下会释放大量元素[11, 33, 10].也有研究表明, 喀斯特地区土地利用和流域的交互作用会显著增加土壤侵蚀, 从而导致土壤中微量元素随径流发生累积和转移[10, 49].由连续变量的空间分布可知, PSD较高的连续变量和土壤重金属来源的空间分布相似.有研究表明, 人类活动是影响土壤重金属含量的主要因素, 靠近河流和矿区的土壤更易受到人类经营性活动的影响[6, 11].此外, 对风险探测出的高风险区域仍需进一步基于实地情况和结合空间分析进行解释.

图 7表 3可知, IDSA模型对各来源贡献迭代计算中所需的连续变量个数分别为4、3、4和5, 其在多变量交互作用中产生单观测值分区的比例为8.06%~13.43%.研究结果表明其较OPGD模型在多变量交互作用的识别过程中更具优势, 但在交互作用的估计中仍然存在不确定性.因此, 接下来需结合其他更具创新的模糊算子用于空间模糊叠加以及结合机器学习算法、多变量回归和最大熵等模型进一步分析.

4 结论

(1) 基于各空间单元的尺度效应比较, 确定本研究最佳空间单元为0.08 km2.通过分位数法进行的多级离散化结果表明, 本研究连续变量的最佳中断数为10.

(2) 因子探测结果表明, 地层和人口分别是分类变量和连续变量中控制土壤重金属来源空间分异的主要驱动因子.

(3) 交互作用探测结果表明, 分类变量中地层和流域的空间交互对土壤重金属来源的解释力分别为F1(58.57%)、F2(49.28%)、F3(27.28%)和F4(60.61%); 连续变量的空间组合对土壤重金属来源空间格局的解释力分别为F1(64.65%)、F2(78.46%)、F3(61.77%)和F4(77.61%).

(4) 风险探测结果表明, 分类变量中各来源的高风险区域分别为地层中下震旦系、上白垩系区域、土地利用中采矿用地和土壤种类中典型强淋溶土区域; 连续变量中各来源的高风险区域分别为ET(41.2~43 kg·m-2)、距河流距离(315~398 m)、EVI(0.796~0.995)和距河流距离(499~605 m).

参考文献
[1] 沈仁芳, 颜晓元, 张甘霖, 等. 新时期中国土壤科学发展现状与战略思考[J]. 土壤学报, 2020, 57(5): 1051-1059.
Shen R F, Yan X Y, Zhang G L, et al. Status quo of and strategic thinking for the development of soil science in China in the new era[J]. Acta Pedologica Sinica, 2020, 57(5): 1051-1059.
[2] 陈卫平, 杨阳, 谢天, 等. 中国农田土壤重金属污染防治挑战与对策[J]. 土壤学报, 2018, 55(2): 261-272.
Chen W P, Yang Y, Xie T, et al. Challenges and countermeasures for heavy metal pollution control in farmlands of China[J]. Acta Pedologica Sinica, 2018, 55(2): 261-272.
[3] 陈能场, 郑煜基, 何晓峰, 等. 《全国土壤污染状况调查公报》探析[J]. 农业环境科学学报, 2017, 36(9): 1689-1692.
Chen N C, Zheng Y J, He X F, et al. Analysis of the Report on the national general survey of soil contamination[J]. Journal of Agro-Environment Science, 2017, 36(9): 1689-1692.
[4] 宋伟, 陈百明, 刘琳. 中国耕地土壤重金属污染概况[J]. 水土保持研究, 2013, 20(2): 293-298.
Song W, Chen B M, Liu L. Soil heavy metal pollution of cultivated land in China[J]. Research of Soil and Water Conservation, 2013, 20(2): 293-298.
[5] 曾思燕, 于昊辰, 马静, 等. 中国耕地表层土壤重金属污染状况评判及休耕空间权衡[J]. 土壤学报, 2022, 59(4): 1036-1047.
Zeng S Y, Yu H C, Ma J, et al. Identifying the status of heavy metal pollution of cultivated land for tradeoff spatial fallow in China[J]. Acta Pedologica Sinica, 2022, 59(4): 1036-1047.
[6] Hou D Y, O'Connor D, Nathanail P, et al. Integrated GIS and multivariate statistical analysis for regional scale assessment of heavy metal soil contamination: A critical review[J]. Environmental Pollution, 2017, 231: 1188-1200. DOI:10.1016/j.envpol.2017.07.021
[7] Zhao F J, Ma Y B, Zhu Y G, et al. Soil contamination in China: current status and mitigation strategies[J]. Environmental Science & Technology, 2015, 49(2): 750-759.
[8] 农业农村部. 2019年全国耕地质量等级情况公报[J]. 中华人民共和国农业农村部公报, 2020(4): 113-121.
[9] 张双银, 吴琳娜, 张广映, 等. 都柳江上游沿岸喀斯特地区土壤重金属污染及健康风险分析[J]. 环境科学学报, 2022, 42(5): 1-14.
Zhang S Y, Wu L N, Zhang G Y, et al. Soil heavy metal pollution analysis and health risk assessment in karst areas along the upper reaches of the Du Liujiang River[J]. Acta Scientiae Circumstantiae, 2022, 42(5): 1-14. DOI:10.13671/j.hjkxxb.2021.0474
[10] Tao H, Liao X Y, Li Y, et al. Quantifying influences of interacting anthropogenic-natural factors on trace element accumulation and pollution risk in karst soil[J]. Science of the Total Environment, 2020, 721. DOI:10.1016/j.scitotenv.2020.137770
[11] 于旦洋, 王颜红, 丁茯, 等. 近十年来我国土壤重金属污染源解析方法比较[J]. 土壤通报, 2021, 52(4): 1000-1008.
Yu D Y, Wang Y H, Ding F, et al. Comparison of analysis methods of soil heavy metal pollution sources in China in last ten years[J]. Chinese Journal of Soil Science, 2021, 52(4): 1000-1008. DOI:10.19336/j.cnki.trtb.2020101202
[12] 陈雅丽, 翁莉萍, 马杰, 等. 近十年中国土壤重金属污染源解析研究进展[J]. 农业环境科学学报, 2019, 38(10): 2219-2238.
CHEN Y L, WENG L P, MA J, et al. Review on the last ten years of research on source identification of heavy metal pollution in soils[J]. Journal of Agro-Environment Science, 2019, 38(10): 2219-2238. DOI:10.11654/jaes.2018-1449
[13] Tobler W R. A computer movie simulating urban growth in the detroit region[J]. Economic Geography, 1970, 46(S1): 234-240.
[14] Odland J. Spatial autocorrelation[EB/OL]. https://researchrepository.wvu.edu/rri-web-book, 2020-09-15.
[15] Cang X Z, Luo W. Spatial association detector (SPADE)[J]. International Journal of Geographical Information Science, 2018, 32(10): 2055-2075. DOI:10.1080/13658816.2018.1476693
[16] Qiao P W, Yang S C, Lei M, et al. Quantitative analysis of the factors influencing spatial distribution of soil heavy metals based on geographical detector[J]. Science of the Total Environment, 2019, 664: 392-413. DOI:10.1016/j.scitotenv.2019.01.310
[17] Song Y Z, Wu P. An interactive detector for spatial associations[J]. International Journal of Geographical Information Science, 2021, 35(8): 1676-1701. DOI:10.1080/13658816.2021.1882680
[18] 周伟, 李丽丽, 周旭, 等. 基于地理探测器的土壤重金属影响因子分析及其污染风险评价[J]. 生态环境学报, 2021, 30(1): 173-180.
Zhou W, Li L L, Zhou X, et al. Influence factor analysis of soil heavy metal based on geographic detector and its pollution risk assessment[J]. Ecology and Environmental Sciences, 2021, 30(1): 173-180.
[19] 龚仓, 王亮, 王顺祥, 等. 基于地理探测器的镇域尺度土壤重金属含量空间分异及其影响因素分析[J]. 环境科学, 2022, 43(10): 4566-4577.
Gong C, Wang L, Wang S X, et al. Spatial differentiation and influencing factors analysis of soil heavy metal content at scale on town level based on geographic detector[J]. Environmental Science, 2022, 43(10): 4566-4577.
[20] 顾高铨, 万小铭, 曾伟斌, 等. 焦化场地内外土壤重金属空间分布及驱动因子差异分析[J]. 环境科学, 2021, 42(3): 1081-1092.
Gu G Q, Wan X M, Zeng W B, et al. Analysis of the spatial distribution of heavy metals in soil from a coking plant and its driving factors[J]. Environmental Science, 2021, 42(3): 1081-1092.
[21] 高浩然, 周勇, 刘甲康, 等. 基于EBK插值预测和GDM模型的襄州区耕地土壤重金属时空分布及来源变化分析[J]. 环境科学, 2022, 43(11): 5180-5191.
Gao H R, Zhou Y, Liu J K, et al. Spatial and temporal distribution and source variation of heavy metals in cultivated land soil of Xiangzhou district based on EBK interpolation prediction and GDM model[J]. Environmental Science, 2022, 43(11): 5180-5191.
[22] Ren Z Q, Xiao R, Zhang Z H, et al. Risk assessment and source identification of heavy metals in agricultural soil: a case study in the coastal city of Zhejiang Province, China[J]. Stochastic Environmental Research and Risk Assessment, 2019, 33(11-12): 2109-2118.
[23] Guo G H, Li K, Zhang D G, et al. Quantitative source apportionment and associated driving factor identification for soil potential toxicity elements via combining receptor models, SOM, and geo-detector method[J]. Science of the Total Environment, 2022, 830. DOI:10.1016/j.scitotenv.2022.154721
[24] Wang J F, Zhang T L, Fu B J. A measure of spatial stratified heterogeneity[J]. Ecological Indicators, 2016, 67: 250-256.
[25] 王劲峰, 徐成东. 地理探测器: 原理与展望[J]. 地理学报, 2017, 72(1): 116-134.
Wang J F, Xu C D. Geodetector: Principle and prospective[J]. Acta Geographica Sinica, 2017, 72(1): 116-134.
[26] Song Y Z, Wang J F, Ge Y, et al. An optimal parameters-based geographical detector model enhances geographic characteristics of explanatory variables for spatial heterogeneity analysis: cases with different types of spatial data[J]. GIScience & Remote Sensing, 2020, 57(5): 593-610.
[27] Dark S J, Bram D. The modifiable areal unit problem(MAUP) in physical geography[J]. Progress in Physical Geography: Earth and Environment, 2007, 31(5): 471-479.
[28] 陈江平, 张瑶, 余远剑. 空间自相关的可塑性面积单元问题效应[J]. 地理学报, 2011, 66(12): 1597-1606.
Chen J P, Zhang Y, Yu Y J. Effect of MAUP in spatial autocorrelation[J]. Acta Geographica Sinica, 2011, 66(12): 1597-1606.
[29] Ye X, Rogerson P. The impacts of the modifiable areal unit problem (MAUP) on omission error[J]. Geographical Analysis, 2022, 54(1): 32-57.
[30] Gao F, Li S Y, Tan Z Z, et al. Understanding the modifiable areal unit problem in dockless bike sharing usage and exploring the interactive effects of built environment factors[J]. International Journal of Geographical Information Science, 2021, 35(9): 1905-1925.
[31] Duque J C, Laniado H, Polo A. S-maup: statistical test to measure the sensitivity to the modifiable areal unit problem[J]. PloS one, 2018, 13(11). DOI:10.1371/journal.pone.0207377
[32] Paatero P, Eberly S, Brown S G, et al. Methods for estimating uncertainty in factor analytic solutions[J]. Atmospheric Measurement Techniques, 2014, 7(3): 781-797.
[33] 张宏泽, 崔文刚, 黄月美, 等. 黔中喀斯特地区临近矿区耕地土壤重金属污染评价及其源解析[J]. 环境科学学报, 2022, 42(4): 412-421.
Zhang H Z, Cui W G, Huang Y M, et al. Evaluation and source analysis of heavy metal pollution of farmland soil around the miningarea of karst region of central Guizhou Province[J]. Acta Scientiae Circumstantiae, 2022, 42(4): 412-421.
[34] Shi T Z, Hu Z W, Shi Z, et al. Geo-detection of factors controlling spatial patterns of heavy metals in urban topsoil using multi-source data[J]. Science of the Total Environment, 2018, 643: 451-459.
[35] Fischer G, Nachtergaele F, Prieler S, et al. Global agro-ecological zones assessment for agriculture(GAEZ 2008)[R]. Rome, Italy: ⅡASA, Laxenburg, Austria and FAO, 2008.
[36] Conrad O, Bechtel B, Bock M, et al. System for automated geoscientific analyses (SAGA) v. 2.1.4[J]. Geoscientific Model Development, 2015, 8(7): 1991-2007.
[37] Kong D D, Zhang Y Q, Gu X H, et al. A robust method for reconstructing global MODIS EVI time series on the Google Earth Engine[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2019, 155: 13-24.
[38] Luo P, Song Y Z, Huang X, et al. Identifying determinants of spatio-temporal disparities in soil moisture of the Northern Hemisphere using a geographically optimal zones-based heterogeneity model[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2022, 185: 111-128.
[39] 徐涵秋. 城市遥感生态指数的创建及其应用[J]. 生态学报, 2013, 33(24): 7853-7862.
Xu H Q. A remote sensing urban ecological index and its application[J]. Acta Ecologica Sinica, 2013, 33(24): 7853-7862.
[40] Pagotto C, Rémy N, Legret M, et al. Heavy metal pollution of road dust and roadside soil near a major rural highway[J]. Environmental Technology, 2001, 22(3): 307-319.
[41] Olaya V, Conrad O. Chapter 12 Geomorphometry in SAGA[J]. Developments in Soil Science, 2009, 33: 293-308.
[42] 陈彦光. 基于Moran统计量的空间自相关理论发展和方法改进[J]. 地理研究, 2009, 28(6): 1449-1463.
Chen Y G. Reconstructing the mathematical process of spatial autocorrelation based on Moran's statistics[J]. Geographical Research, 2009, 28(6): 1449-1463.
[43] Getis A. Spatial interaction and spatial autocorrelation: a cross-product approach[J]. Environment and Planning A: Economy and Space, 1991, 23(9): 1269-1277.
[44] Carr D B, Olsen A R, White D. Hexagon mosaic maps for display of univariate and bivariate geographical data[J]. Cartography and Geographic Information Systems, 1992, 19(4): 228-236.
[45] 李雨, 韩平, 任东, 等. 基于地理探测器的农田土壤重金属影响因子分析[J]. 中国农业科学, 2017, 50(21): 4138-4148.
Li Y, Han P, Ren D, et al. Influence factor analysis of farmland soil heavy metal based on the geographical detector[J]. Scientia Agricultura Sinica, 2017, 50(21): 4138-4148.
[46] 陈艺, 蔡海生, 曾君乔, 等. 袁州区表层土壤重金属污染特征及潜在生态风险来源的地理探测[J]. 环境化学, 2021, 40(4): 1112-1126.
Chen Y, Cai H S, Zeng J Q, et al. Characteristics of heavy metal pollution in surface soil and geographical detection of potential ecological risk sources in Yuanzhou District[J]. Environmental Chemistry, 2021, 40(4): 1112-1126.
[47] Openshaw S. Ecological fallacies and the analysis of areal census data[J]. Environment and Planning A: Economy and Space, 1984, 16(1): 17-31.
[48] Meng X Y, Gao X, Lei J Q, et al. Development of a multiscale discretization method for the geographical detector model[J]. International Journal of Geographical Information Science, 2021, 35(8): 1650-1675.
[49] Liang S Z X, Fang H Y. Quantitative analysis of driving factors in soil erosion using geographic detectors in Qiantang River catchment, Southeast China[J]. Journal of Soils and Sediments, 2021, 21(1): 134-147.
[50] Ye H, Sun C G, Wang K, et al. The role of urban function on road soil respiration responses[J]. Ecological Indicators, 2018, 85: 271-275.