沱江上源支流土壤重金属污染空间相关性及变异解析

任加国，王 彬，2，3，师华定，武倩倩

（1.山东科技大学地球科学与工程学院，山东 青岛 266590；2.生态环境部土壤与农业农村生态环境监管技术中心，北京100012；3.中国环境科学研究院土壤与固体废物环境研究所，北京 100012）

土壤是生态系统的基本要素之一，它不仅直接影响到国民经济发展和国土资源环境安全，而且直接关系到农产品安全和人体健康[1-2]。土壤重金属污染具有隐蔽性、滞后性、累积性、不可逆转性和难治理性，其污染成因错综复杂。因此，辨识土壤重金属污染空间变异特征，解析其变异成因就显得尤为重要[3-4]。

近些年来，国内外学者在土壤重金属空间相关性及变异性方面进行了大量研究[5-11]，主要是利用Moran's I指数、变异系数、半变异函数理论、尺度方差分析法、地统计学方法等对研究区土壤重金属的空间相关性、空间结构、变异规律及成因进行分析。其中，Moran's I指数和变异系数均可与半变异函数理论结合在一起使用，例如秦治恒等[5]利用Moran's I指数和半变异函数理论揭示了研究区土壤Cd含量的空间相关性及其在空间结构上的变异规律，而李传章等[6]和Chandra等[7]利用变异系数和半变异函数理论探究了研究区土壤重金属变异强弱及空间自相关性范围。尺度方差分析法是以半变异函数理论为基础，对不同空间等级尺度上的土壤重金属变异进行探究，例如刘伟等[8]为了研究土壤重金属的空间分异性及其空间结构特征，采用半方差分析法和尺度方差分析法定量刻画了土壤Cd空间变异的等级结构和特征尺度。地统计学方法也是分析空间变异的有力工具之一，例如Wu等[9]、Lin等[10]和叶露萍等[11]利用该方法对土壤中元素含量的空间分布进行插值预测，并通过分析土壤本身性质及外界影响因素，对其空间分布特征及变异规律进行定量分析。通过以上文献分析可知，Moran's I指数用于度量不同位置上样点数据之间的相互依赖程度，而半变异函数理论可以进一步解释样点数据在不同距离和方向上的变异程度，且利用地统计学方法对样点数据进行插值分析时也需根据半变异函数理论对其进行拟合，因此三者可以结合起来对样点数据进行不同尺度的分析。另外，在以上研究中均没有考虑到局部范围的聚类特征，且对于利用Anselin Local Moran's I指数来分析土壤重金属含量空间聚类特征方面的研究仍未见报道，仅有部分学者利用该指数来分析动植物及其他方面的空间聚类特征[12-14]。因此本文以沱江上源流域土壤为研究对象，将Global Moran's I指数、Anselin Local Moran's I指数、半变异函数理论以及地统计学方法等结合在一起，充分发挥不同方法的功能性，全面分析研究区土壤重金属空间相关性及其变异规律，旨在为区域土壤后续的污染控制及治理提供理论基础。

1 材料与方法

1.1 研究区概况

沱江是长江上游支流，位于中国四川省中部，发源于九顶山，上源分三条支流即东源绵远河、中源石亭江以及西源湔江。研究区主要土地利用类型为农田，部分区域为成林和灌木林，属于亚热带季风气候和半湿润气候，雨量充沛，特别适合农作物的生长，且沿河两岸十几万亩（1亩=667 m2）农田均由沱江灌溉。沱江是四川省工业城市最集中的河流，为众多企业提供生产、生活用水。

1.2 样品的采集、制备与测定

利用ArcGIS 10.5软件进行研究区点位布设，原则上按照每5000 m×5000 m网格布设一个点位，网格的中心点作为土壤采样点，根据地形及周边企业情况作适当调整。共布设137个表层土壤点位，绝大部分点位布设集中在农田，仅少数点位布设于成林，如图1。根据网格中心采样点经纬度，利用GPS对每个计划采样点进行精准定位，一般误差不超过100 m，采样人员到达目标点位后，必须观察其是否符合土壤采样的代表性要求，在允许范围内优选采样点，位移距离一般不超过50 m，且陡坡地、低洼积水地、住宅、道路、沟渠、粪坑附近等不设采样点。

样品采集以确定点位为中心划定采样区域，一般为100 m×100 m，采用梅花形混合采样法，采样时先用铁铲切割一个大于取土量的20 cm深的土方，再用木铲去掉铁铲接触面后装入样品袋，不斜向挖土，尽可能做到采样量上下一致，之后将采集完的土壤样品置于阴凉处保存。土壤样品制备时，首先取适量新鲜土壤样品平铺在干净的搪瓷盘或玻璃板上，避免阳光直射，且环境温度不超过40℃，自然风干，去除石块、树枝等杂质，将土块粉碎后过100目尼龙筛，并混合均匀待测。采用电感耦合等离子体质谱法（ICP-MS）测定Cd、Pb和Cr的含量，原子荧光法测定Hg和As的含量，分析方法精密度和准确度采用国家土壤标准物质GSS-15和室内平行样品进行质量控制，回收率在99%～105%。

图1 研究区采样点位布设及土地利用类型Figure 1 Sample site laying and land use types in the study area

1.3 研究方法

1.3.1 空间自相关分析

在地统计学上完全随机的样点是没有意义的，而空间上某一点土壤重金属含量也并非是随机变化的，它不仅与样点的距离、位置有关系，而且与其传播途径和方式具有重大的联系[15-17]。在判断空间自相关性时，除考虑整个研究区大尺度趋势外，也要兼顾局部效应所反映的问题，因此需要用Global Moran's I指数来评估全域范围内的空间自相关性，并利用Anselin Local Moran's I指数找出局部范围内的热点、冷点以及异常值[18]，其计算公式分别为：

式中：n为空间要素总数，wij为空间权重矩阵，xi和xj分别为空间变量x在不同位置i和j上的实际观测值，x为空间变量x的平均值。当I＞0时，表示空间呈正相关，其值越大，空间相关性越明显；当I=0时，表示空间不具有相关性；当I＜0时，表示空间呈负相关，其值越小，空间差异越明显。而根据Ii可区分具有统计学上的显著性（P=0.05水平）的高值（HH）聚类、低值（LL）聚类、由高值围绕的低值异常值（LH）聚类以及由低值围绕的高值异常值（HL）聚类等。

1.3.2 半变异函数分析

半变异函数又称半变差函数，是地统计分析的特有函数，可以对区域化变量的空间变异特征进行较为准确的描述[19]。半变异函数可以揭示区域化变量的内在联系，直接测定和分析变量的空间相关性和依赖性，研究变量的空间分布规律[20-21]，且能定量地呈现随着距离的变化而不断改变的空间性质[22]。其计算公式为：

式中：N（h）表示空间上具有相同步长向量h时的离散点对数量，Z（xi）和Z（xi+h）分别表示位置在xi和xi+h上的区域化变量值。半变异函数值随着距离的加大而增加，而在离散点对间隔为0和半变异函数趋于平稳时会产生四个参数，即块金值、变程、偏基台值、基台值。当步长h为0时，其半变异函数也应为0，但由于随机因素的影响，使得半变异函数r（h）不为0，此时r（h）为块金值；随着步长h的增大，离散点之间的相关性越来越小，半变异函数r（h）趋于稳定状态，此时r（h）为基台值，h为变程；偏基台值为基台值和块金值之差。

2 结果与讨论

2.1 土壤重金属含量统计

研究区土壤重金属含量统计情况见表1。对137个土壤样品的分析结果表明，土壤pH的平均值和中位值分别为6.75和6.88，说明研究区土壤整体上呈弱酸性。Cd、Hg、As、Pb、Cr 5种重金属含量的平均值均未超过国家土壤环境质量标准（GB 15618—2018）风险筛选值，且高于中位值，说明研究区土壤重金属含量大部分处于较低水平。但Cd、Hg、As三种重金属含量的最大值均超过了风险筛选值但低于风险管制值，说明这三种重金属对研究区的污染相对较重；而Pb、Cr的含量最大值均未超过风险筛选值，说明这两种重金属污染较轻。从峰度来看，Cd、Hg、As、Pb、Cr 5种重金属数据的分布均比正态分布高耸狭窄，说明其数据比正态分布更集中于平均值附近；从偏度来看，Cd、Hg、As、Pb、Cr 5种重金属数据分布为正偏度（右偏度），但经过Log变换之后均符合正态分布。

表1 土壤重金属含量统计分析Table 1 Statistical analysisof heavy metals content in soil

2.2 相关性分析

2.2.1 全域范围内相关性分析

利用ArcGIS 10.5软件计算得到空间关联指数Global Moran's I及相关参数，见表2。从Global Moran's I指数来看，研究区内土壤 Cd、Hg、As、Pb、Cr 5种重金属均呈现空间正相关的特点，即空间上越相近，其相关性越强，按照相关性大小依次排序为：Pb＞Cd＞Cr＞As＞Hg；从Z得分（表示标准差的倍数）和P值（产生随机模式的概率）来看，只有土壤Hg元素含量在空间上可能会产生随机分布，其余元素含量随机分布的可能性均为0，尤其是土壤Pb、Cd的Z得分较高，说明其空间分布存在明显的聚类特征。

表2 Global Moran's I指数及相关参数Table 2 Global Moran's I index and related parameters

2.2.2 局部冷热点及异常值分析

与Global Moran's I指数相比，Local Moran's I指数更侧重于分析局部范围内观测值的聚类特征以及推断异常值的位置，二者互为补充和说明，彼此之间具有一定的相似性及其特殊性。Bivand等[23]对比了全局关联法和局部关联法在度量空间自相关时所存在的差异，指出使用全局关联法度量时应考虑到整个数据集，而局部关联法是针对每个独立体且可以对统计结果进行推断；Boots等[24]研究了Global Moran's I指数和Local Moran's I指数在不同空间模式下的特征，说明了Global Moran's I指数适合研究区域的总体度量，而Local Moran's I指数用于特定的区域及其邻近区域，适合局部范围观测值特征的确定。本文在分析研究区全域范围内土壤重金属空间相关性的前提下，为进一步探索其局部空间分布的特征，利用ArcGIS 10.5软件ArcToolbox工具箱中的聚类分布制图，运行Anselin Local Moran's I分析模块，并设置相应参数，得到研究区土壤重金属含量聚类分布，见图2。从高值聚类点（热点）来看，Cd、Hg、As、Pb、Cr 5种重金属均存在，主要分布在绵远河西北方向及湔江以西方向，而绵远河和石亭江之间也有少量高值聚类点；从低值聚类点（冷点）来看，Cd、As、Cr三种重金属分布数量较多且主要分布在石亭江两岸及湔江西南方向，而Hg主要分布在绵远河周围，Pb分布比较零星；从被高值围绕的低值异常点来看，Cd、Hg、As、Pb、Cr五种重金属分布数量较少，湔江周围五种重金属均有分布，石亭江附近仅有Cd、Hg分布，而绵远河附近有Cd、Hg、As分布；从被低值围绕的高值异常点来看，Cd、Hg、As、Pb、Cr 5种重金属分布数量均最少，绵远河附近有少量Hg、As、Cr的此类异常点分布，石亭江有Cd、As、Pb、Cr 4种重金属的此类异常点各1个，湔江有Hg、As、Pb 3种重金属的此类异常点各1个。

2.3 空间结构分析

2.3.1 半变异函数云分析

大部分地理现象均具有空间相关特性，即距离越近的事物其特性越相近，而半变异函数云就是这种相似性的定量化表示[25]。利用ArcGIS10.5软件对研究区土壤Cd、Hg、As、Pb、Cr 5种重金属含量进行半变异函数云分析，由于空间各向异性较为普遍，需要在搜索方向为0°（正南北方向）、45°（东北-西南方向）、90°（正东西方向）以及135°（西南-东北）方向上（角度容差为±22.5°）观察半变异函数云的变化情况，但通过比较，5种重金属含量的半变异函数云仅在0°（正南北方向）和90°（正东西方向）方向上差异略大，如图3。研究发现，半变异函数云中距离较近（靠近X轴左侧）的样点对具有更大的相似性（Y轴下方）；反之，半变异函数云中距离较远（朝X轴的右侧移动）的样点对的方差更大（朝Y轴的上方移动），相似性较低；而当样点对形成一条水平直线时，则数据将不存在空间相关性。从5种重金属含量的半变异函数云的样点对分布状况可知，它们的大部分样点对处于聚集的状态，仅存在部分的离散样点对，说明它们在一定的距离内（步长）具有一定的空间相似性，而随着距离（步长）的增加，其空间相似性逐渐降低；从5种重金属含量的半变异函数云的样点对的方向上来看，它们在90°方向上（正东西）比0°方向上（正南北）均具有更远距离（步长）的空间相关性。

2.3.2 半变异函数拟合结果分析

图2 研究区土壤重金属含量聚类分布Figure 2 Cluster distributions of heavy metals content in soil of study area

对样点数据进行半变异函数拟合，并选择最优的拟合模型，对于判断样点的空间相关性强弱及范围尤为重要。于冬雪等[26]利用半变异函数对黄土区不同土壤深度的容重进行拟合，得出不同深度土壤容重的变程即空间相关性范围为22～780 km；席小康等[27]基于半方差函数理论，对锡林河流域不同深度土壤有机碳含量进行了半变异函数拟合，研究表明随着土层深度的增加，土壤有机碳空间分布相关性增强的结论。本文利用GS+for windows9.0软件对样点数据进行半变异函数拟合，而半变异函数只有在最大间隔1/2以内才有意义，研究区中最大的采样点距离为72 133.89 m，则设定有效滞后距离为36 066.94 m，数据设定为10组，组距为4000 m，选择决定系数较高且残差最小的拟合模型作为最优模型，模拟结果见表3。从拟合结果来看，除了Hg元素外其余元素的决定系数均高于0.8，说明各模型拟合精度较好满足要求。Cd、Hg、As、Pb四种重金属的块金值较低，说明随机部分（样品采集、测量误差等）引起的空间变异可以忽略不计。Cd、Pb元素的基底效应低于0.25，说明其空间连续性较好，空间相关性程度强烈；As、Cr元素的基底效应高于0.25但低于0.75，说明其空间连续性一般，空间相关性程度为中等；Hg元素的基底效应最大，说明其空间连续性较差，空间相关性程度弱，其空间分布很可能受到了一定的外界干扰，随机分布的可能性较大[28]。变程表示空间自相关范围的大小，5种重金属变程从大到小为Pb＞Cd＞As＞Cr＞Hg，其中Pb、Cd元素的变程超过50 000 m，而As、Cr元素的变程也超过30 000 m，说明这4种元素在较大范围内都具有一定的相关性；Hg元素的变程相对较小且低于10 000 m，说明该元素只在较小的范围内存在空间相关性。

图3 研究区土壤重金属半变异函数云Figure 3 Semi-variation function cloud of heavy metals in soil of study area

续图3研究区土壤重金属半变异函数云Continued figure 3 Semi-variation function cloud of heavy metals in soil of study area

表3 半变异函数理论模型及参数Table 3 Theoretical models and parameters of semi-variational function

2.4 空间分布预测及变异成因分析

2.4.1 全局趋势分析

全局趋势是指空间数据在某一特定方向上的变化趋势，反映了空间变量在区域上变化的主体特征，用来揭示空间变量的总体规律。地统计学应用的条件是必须满足二阶平稳假设，而全局趋势会破坏该假设成立的条件，因此在利用地统计学进行插值分析时需要剔除全局趋势[29]。利用ArcGIS10.5软件的全局趋势分析工具将研究区二维平面上的点以属性值为高转化成三维视图，并将点投影在两个互相垂直的平面上，并运用多项式对投影进行拟合，如图4。从图中可以看出，Cd元素的拟合曲线较为平直，不存在全局趋势；Hg元素在东西方向上存在微弱的全局趋势；Pb、Cr两种元素在南北方向上呈U型二次拟合曲线，存在明显的全局趋势；As元素在南北方向上呈U型二次拟合曲线，也存在明显的全局趋势。

2.4.2 污染空间分布预测

克里金插值是以变异函数理论和结构分析为基础，在有限的区域内对区域化变量进行无偏最优估计的一种方法[30-31]。由上文可知，研究区土壤重金属污染空间在距离和方向上均存在自相关性，因此可以采用普通克里金方法进行内插或外推，从而对研究区污染空间分布进行预测。打开ArcGIS10.5软件的地统计向导工具，选择对应的数据集和普通克里金插值方法之后，对数据进行Log转换并剔除全局趋势，最后根据前文对半变异函数云及拟合模型的分析，选用相应的参数，得到污染空间分布预测结果，见图5。从Cd元素来说，其高值区域主要分布于绵远河、湔江、石亭江流域的西北部，而低值区域与之相反；从Hg元素来说，其高值区域主要分布于绵远河、石亭江两流域之间，而低值区域主要分布于湔江流域周围；从As元素来说，其高值区域沿着研究区的西南-东北方向穿越湔江、石亭江和绵远河流域，低值区域则被高值区域所包围；从Cr元素来说，其高值区域主要分布于石亭江、绵远河两流域之间及三条流域的交汇处（东南部），而大部分低值区域被高值区域所包围，部分低值区域分布于湔江流域；从Pb元素来说，其高值区域仅存在于石亭江和绵远河流域之间，而低值区域分布较散，并在多处区域形成聚集区。整体上来说，五种重金属元素含量均存在聚集现象，尤其是低值区域的聚集较为频繁，高值区域的聚集相对较少，说明研究区整体污染水平较低，土壤环境质量良好，但对于高值聚集区域也应引起一定的重视。

图4 土壤重金属元素全局趋势分析Figure 4 Analysis of the global trend of heavy metals in soil

2.4.3 变异成因分析

图5 土壤重金属污染空间分布预测Figure 5 Prediction of spatial distribution of heavy metals pollution in soil

由前文可知，该研究区土壤重金属有较为明显的聚类现象和空间异质性，且存在高值和低值异常点，因此需要进一步对其变异成因进行分析。本研究搜集了研究区相关企业信息及实际调研情况，见图6。有研究表明[32-36]，土壤重金属污染相对较重的区域其主要特征是企业分布数量较多、交通发达、污水灌溉严重、施肥过量、大气沉降集中或者属于地质高背景区域。该研究区土壤重金属含量的高值主要聚集在绵远河和石亭江流域周围以及三条流域的交汇处（东南部），而该区域的企业分布也最为密集，主要有化学原料和化学制品制造业、医药制造业、金属制品业、黑色金属冶炼和压延加工业以及皮革、毛皮、羽毛及制鞋业等，说明土壤Cd、Hg、As、Pb、Cr 5种重金属的空间分布与此类型企业分布具有一定的正相关性，尤其是化学原料和化学制品制造业与污染空间分布预测的高值区域位置吻合度较高，对于该企业类型应当引起一定的重视，且有相关研究表明[37-38]，该流域的水质受到了一定程度的重金属污染，由于研究区土地利用类型主要为农田，考虑到居民对农田进行施肥并利用该流域的水资源进行农田灌溉，污水灌溉和大量施肥也可能是引起高值区域聚集的重要因素；低值区域主要分布于湔江流域周围，发现其周围企业分布数量较少，仅有少量的化学原料和化学制品制造业和金属制品业，该区域土壤重金属含量较低的原因可能是其地质背景值较低且受外界因素的影响较小；从研究区异常值来说，其分布比较分散，发现Cd、As两种重金属的低异常值点出现在企业密度较大的石亭江流域的东南部和绵远河流域的东北部，Cr的低异常值出现在仅有一家企业的湔江流域西部，而Hg、Pb两种元素的高异常点却出现在企业密度较小的湔江流域附近，这些异常值的出现与企业并不呈现一定的相关性，因此考虑到自然因素对异常值的影响，其分布很可能与成土母质、地形地貌、大气沉降、交通要素等有关，而与企业的分布情况并没有直接的关系。笔者认为，空间变异的成因相对比较复杂，由于土壤污染源扩散的不确定性以及污染源汇入能力的不一致性，它不仅与外界因素息息相关，而且与土壤内物质本身的迁移规律存在一定的关系，因此在之后的土壤重金属空间变异的研究当中，应基于土壤污染的源-汇关系来进一步分析空间变异的成因。

图6 研究区周围企业分布状况Figure 6 Distribution of enterprises around the study area

3 结论

（1）研究区土壤Cd、Hg、As、Pb、Cr 5种重金属含量大部分处于较低水平，其平均值未超过GB 15618—2018国家土壤环境质量标准风险筛选值。

（2）相关性分析表明，研究区土壤Cd、Hg、As、Pb、Cr 5种重金属在空间上均呈现正相关的特点，即空间上越相近，其相关性越强；除Hg外其他重金属产生聚类特征的可能性均很大，但从局部冷热点及异常值分析可知，Hg也存在聚类特征。

（3）空间结构分析表明，研究区土壤Cd、Hg、As、Pb、Cr 5五种重金属90°（正东西）方向比0°方向上（正南北）均具有更远距离（步长）的空间相关性，其中Cd、Pb的空间连续性较好、相关性强烈且变程最长，As、Cr次之，Hg最差。

（4）研究区空间污染的高值区域与周围企业分布的数量具有一定的正相关性，且污水灌溉和大量施肥也是引起高值区域产生的重要因素；低值区域与企业分布数量呈现负相关性；而异常值区域的分布较为分散，与周围企业的分布没有直接的关系，其成因很可能与自然因素（如成土母质、地形地貌、大气沉降、交通要素等）有关。