随着区域工农业快速发展，金属矿山开采、化肥农药施用和污水灌溉成为威胁土壤环境质量的重要因素. 金属矿山在开采过程中，矿体的暴露往往使得伴生元素被释放到环境中；矿渣、尾矿无序堆放与不当处理也会造成重金属在土壤中的积累. 积累初期重金属的毒理效应并不明显，当累积量超过土壤自净能力时，便会造成土壤重金属污染，并通过食物链被人体吸收，进而威胁人体健康^[1]. 因此对矿区土壤重金属进行空间分布特征和环境风险研究对重金属污染防控和土壤修复有重要意义.

目前，空间统计分析，包括空间自相关与地统计，与GIS被广泛应用于土壤重金属空间结构特征分析与制图中，许多学者对此开展了卓有成效的研究^[1~6]. 空间自相关分析不仅可以研究整个区域的空间自相关程度与尺度，还可以通过区域自相关探索重金属含量的空间集聚区，预估潜在的污染风险. 常见的内梅罗指数、地累积指数、污染负荷指数、潜在生态危害指数等土壤重金属污染评估方法均以土壤元素的自然背景值为参考，仅能对采样点的污染状况进行评价，而无法推广到未采样地区^[7~9]. 为了获得未采样地区的土壤环境质量状况，通常会对评价结果进行空间插值^{[7, 10~12]}，而指示克里金法为区域土壤重金属潜在污染风险的评估提供了新思路.

焦家式金矿为山东省典型金矿区，由于金矿体内与冶炼废弃物中多含有Pb、Zn、Cu和Cr等元素，易于被人为富集和扩散，金矿的集中与城镇人口的密集分布也会在一定程度上加速重金属在土壤中的积累^[13]. 但目前尚未有针对研究区土壤重金属空间特征及环境风险评估的研究，因此本文依托空间统计分析与GIS技术，对金矿区内Cu、Pb、Zn和Cr的空间结构特征与分布规律进行分析；利用热点分析寻找高值集聚点与低值集聚点，探索潜在污染风险区位置；以研究区内土壤元素自然背景值为基准，通过指示克里金技术获取整个区域重金属超标的概率，并验证与热点分析得到的风险区的一致性，以期为研究区土壤环境质量监测和预警提供科学参考.

研究区位于山东省烟台市西北部，地跨莱州和招远两个县级市，主要涉及金城、朱桥、辛庄和蚕庄这4个镇，地理坐标为120°05′~120°15′E，37°20′~37°32′N，属于温带季风性气候，年平均降水约为600 mm；土壤类型以棕壤土和褐土为主；地势东南高西北低，东南多丘陵，西北多为平原，海拔介于-2~298 m之间. 研究区位于我国焦家式金矿带^[14]，受北东向焦家断裂严格控制，多期多因的强烈岩浆活动使得区域内金矿资源十分丰富，开采和保有储量均为全国第一；其中焦家金矿早在1980年就建成并投入生产，是我国黄金行业生产规模最大的现代化矿山之一，生产规模日处理矿量可达1 500 t.

根据研究区土地利用类型及矿区分布情况，同时考虑道路的通达性和样点分布的均匀性，利用研究区遥感影像图和手持GPS进行野外调查和选点. 本研究共采集了77个土壤样品，采样工作分两次完成，2014年11月采集35个土样，2015年3月又增加了42个土样. 研究区及采样点位置见图 1. 样方内均采用多点混合法采样，即沿样方对角线采点，在对角线4等分点处各选一个采样点，共5个采样点，每点用木铲取0~20 cm表层土，剔除杂物，均匀混合后装入聚乙烯塑料袋中，采样同时记录坐标.

图 1

Fig. 1

图 1 研究区与采样点位置示意 Fig. 1 Location of the study area and the sampling sites

将采集的土样带回实验室后，置于木盘中通风阴干，用玛瑙研钵压碎并过100目尼龙筛. 土壤样品经HNO₃-HCL-HF-HCLO₄消解，采用电感耦合等离子体原子发射光谱法(ICP-AES)对Cu、Pb、Zn和Cr的含量进行测定. 分析所用试剂均为优级纯，测试精度采用国家标准土壤样品(GBW07401)检验，回收率介于100%±10%之间，测试结果符合研究要求.

空间自相关分析是检验具有空间位置的变量，其观测值是否与相邻空间点观测值显著相关的方法^[15]，通常分为全局性与区域性. 全局性分析是对空间自相关的全局评估，但往往会掩盖反常的局部状况与小范围的不稳定性，而区域性分析能更好地发掘随机因素影响下的空间分布模式^[16]. 热点分析(getis-ord G_i^*)是局部空间自相关分析的方法之一，以全局空间集聚为前提，通过计算G_i^*统计量揭示集聚模式并反映高值或低值要素在空间上发生聚类的位置. 显著的正G_i^*表示在该区域单元周围，高观测值的区域单元趋于空间集聚，即“热点”(hotspot)；而显著的负G_i^*值表示低观测值的区域单元趋于空间集聚，即“冷点”(coldspot).

地统计以区域化变量理论为基础，以半方差函数为工具，是研究土壤空间变异的范式性方法. 普通克里金是最常用的地统计方法，可以分析土壤重金属含量的空间连续性变异；而指示克里金作为非参数地统计方法，可以估计土壤重金属在给定位置超出或小于规定阈值的概率(P)^[17]，进而表征其环境风险. 普通克里金要求数据符合正态分布，因此利用Kolmogorov-Smirnov(K-S)检验法检验数据的正态分布情况.

本研究利用SPSS 20.0软件(IBM Inc.)对研究区表层土壤内重金属含量进行描述性统计分析、正态分布检验以及Pearson相关性分析；空间结构分析在GS+9.0(Gamma Inc.)中完成；重金属含量空间分布、热点分析以及风险概率评估均借助ArcGIS 10.2(Esri Inc.)中的地统计模块(Geostatistial Analyst)实现.

土壤重金属含量描述性统计结果(表 1)表明，Zn、Pb和Cr符合正态分布，但Pb和Cr的K-S检验显著性较差，含量值正偏态明显，特别是Cr，偏度达3.3；Cu为非正态分布，经对数转换(logarithm transformation)后Asymp.Sig.为0.253，达到正态分布标准. 当数据本身为正态分布时，以算术平均值作为数据均值；当数据变换后符合正态分布的则使用几何平均值，因此Cu、Pb、Zn和Cr含量平均值分别为19.41、27.32、49.81和39.27 mg·kg^-1. 除Pb外，其余重金属含量平均值均低于研究区自然背景值^[13]，表现贫化，Pb相对富集. 重金属Cr和Zn含量值的分布范围较大，对应极差达到134.23 mg·kg^-1和80.5 mg·kg^-1，相应标准差也表现出两者含量的高离散程度；Cu和Pb的含量数据则分布相对均匀.

表 1 (Table 1)

表 1 土壤重金属含量描述性统计分析1)/mg·kg-1 Table 1 Descriptive statistical analysis of heavy metal contents in soils /mg·kg-1 重金属 最大值 最小值 极差 均值 标准差 背景值 峰度 偏度 K-S检验(Asymp.Sig.) Cu 55.5 7.1 48.4 19.41a 8.41 24.5 5.20 1.93 0.01 Pb 50.3 17.1 33.2 27.32b 5.18 27.2 4.75 1.39 0.13 Zn 100.3 19.8 80.5 49.81b 13.13 60.4 2.71 0.77 0.49 Cr 153.4 18.8 134.6 39.27b 18.32 55.8 18.88 3.30 0.14 1)a代表几何平均值，b代表算术平均值

表 1 土壤重金属含量描述性统计分析1)/mg·kg-1 Table 1 Descriptive statistical analysis of heavy metal contents in soils /mg·kg-1

Pearson相关性分析(表 2)表明，Cu、Pb、Zn和Cr之间均存在0.01显著性水平上的相关关系；Pb-Zn、Pb-Cu相关系数大于0.5，表现为较强的相关性，Chen等^[18]、解宪丽等^[19]和戴彬等^[20]等学者研究同样发现三者相关系数较高，表明可能拥有共同的来源，通常源于交通及工业活动等. Cr与Pb、Cu虽显著相关，但与Zn的相关系数更高(r=0.61)，两者更有可能受相同因素控制；结合描述性统计结果(表 1)可知，Cr和Zn的平均值低于自然背景值，受人类活动干扰小；Borůvka等^[21]、赵彦锋等^[22]的研究结果表明Cr和Zn受地质背景控制，属于自然来源.

表 2 (Table 2)

表 2 土壤重金属相关性分析1) Table 2 Correlation analyisis between heavy metals Cu Pb Zn Cr Cu 1 Pb 0.552** 1 Zn 0.450** 0.528** 1 Cr 0.338** 0.348** 0.611** 1 1)**代表在0.01水平下显著相关

表 2 土壤重金属相关性分析1) Table 2 Correlation analyisis between heavy metals

土壤重金属的空间变异结构利用半方差函数进行拟合. 变异模型参数中，块金值(Nugget)由测定误差和小于最小采样尺度的非连续性变异引起，属于随机变异；基台值(Sill)通常表示系统内的总变异，它是结构性变异和随机性变异之和. 块金值与基台值之比是反映区域化变量空间异质性程度的重要指标，又称为块金效应(Nugget effect)，该比值用以反映空间变异影响因素中区域因素(自然因素)和非区域因素(人为因素)的作用. 当C₀/(C₀+C)<0.25时，变量以结构性变异为主，变量具有强烈的空间自相关；当比值大于0.25而小于0.75时，变量为中等程度空间自相关；当系数>0.75时，则以随机变异为主，变量的空间自相关性微弱^[23]. 而变程(Range)表示了在某种观测尺度下，空间自相关的作用范围，即在变程范围内，变量有空间自相关性，反之则不存在.

本研究中首先对所有重金属含量进行了对数转换，并根据决定系数(R²)最大和残差平方和(RSS)最小的原则，选取最佳空间变异模型. 由半方差函数理论模型拟合结果(表 3)可知，所有模型的均RSS较小且R²达到0.8以上，效果理想；Cu、Zn和Cr符合球状模型，Pb则符合指数模型；Cu和Pb的块金系数介于0.25和0.75之间，表现为中等程度的空间自相关，其含量可能受到自然和人为因素的双重作用；Zn和Cr的块金系数均小于0.25，以结构性空间变异为主，说明两种元素在土壤中的含量受内在因素(如土壤类型、母岩、地形等)控制，这与Pearson相关性分析结果一致. 各变量的变程介于8 050~19 260 m之间，大于实际最大采样的间隔2 284 m，说明各样点均在空间变异范围内，能够反映重金属的真实变异特征^[24]. Pb的变程最大，表明在更大尺度上存在空间自相关性，其次为Cu、Cr、Zn，属于中等尺度空间自相关.

表 3 (Table 3)

表 3 土壤重金属半方差函数理论模型拟合结果 Table 3 Fitting results of semivariance models of heavy metals in soils 项目 理论模型 块金值(C0) 基台值(C0+C) 块金系数[(C0+C)/C0] 决定系数 残差平方和 变程/m lg Cu Spherical 0.076 0.153 0.497 0.849 0.001 20 13 450 lg Pb Exponential 0.017 0.042 0.394 0.946 0.000 02 19 260 lg Zn Spherical 0.017 0.085 0.204 0.831 0.001 00 8 050 lg Cr Spherical 0.010 0.173 0.059 0.953 0.001 45 9 050

表 3 土壤重金属半方差函数理论模型拟合结果 Table 3 Fitting results of semivariance models of heavy metals in soils

重金属含量空间分布情况如图 2所示，Cu和Pb呈岛状的空间分布特征，Zn和Cr含量表现为条带状和岛状结合的特点. Pb含量的高值区主要位于金城镇和辛庄镇东部；低值区则位于研究区南部、蚕庄镇中部及辛庄镇西部等部分地区. Cu空间分布特点与Pb基本一致，但Cu元素在朱桥镇北部存在明显高值区，且岛状区域多而零散，这主要与Cu元素空间变异性明显有关. Zn和Cr含量全局性规律十分明显，表现为由研究区中部向两侧递减的梯度变化，高值区位于金城镇和辛庄镇与蚕庄镇的交界带，低值区位于辛庄镇北部与研究区东南部.

图 2

Fig. 2

图 2 土壤重金属含量空间分布 Fig. 2 Spatial distributions of heavy metals in soils

Cu、Pb、Zn和Cr含量超过自然背景值的区域面积分别占研究区总面积的82%、50%、96%和91%，含量总体较低. 一般高值区均位于研究区中部的金矿条带区域，而在条带两侧，特别是东南部丘陵区含量则较低；金矿越密集的区域重金属含量也越高. 这与Rashed^[25]、徐友宁等^[26]、高彦鑫等^[27]的研究结果一致，可以认为金矿开采和冶炼废弃物排放对土壤元素的含量和分布存在一定影响.

Cu、Pb、Zn和Cr的全局Moran's Ⅰ指数分别为0.06、0.17、0.24和0.18，均表现为集聚的空间模式，即相似值(高值或低值)在空间上趋于集中分布. 为进一步区分集聚方式，分别对77个采样点上4种重金属含量的G_i^*统计量进行计算，得到了不同显著性水平下的“热点”与“冷点”(图 3). 重金属含量的热点区域为土壤环境研究的重点区域，存在一定的潜在污染风险. 根据热点分析结果，Cu的热点数量为12、Pb为7、Zn为15，Cr为28，Cr的污染风险点最多且0.01显著性水平的点所占比例最高；Cu和Pb的高-高集聚点全部位于金城镇内，低-低集聚点多分布于蚕庄镇内：Cu有6个冷点位于其中部和南部，Pb有3个冷点位于南部边界地区；Zn和Cr的热点分布范围较Cu和Pb更为广泛，呈条带状分布于研究区中部：包括金城镇以及辛庄镇和蚕庄镇的交界地带. 总体来说，热点均位于研究区中部金矿条带内，同样与金矿分布位置与密度表现出较高的相关性.

图 3

Fig. 3

图 3 土壤重金属热点分布与环境风险 Fig. 3 Sptial distributions of hotspots and risk probability of heavy metal contamination in soil

土壤重金属环境风险概率如图 3所示，本研究依据等间距法^{[1, 28]}将环境风险分为无风险(P≤0.2)、低风险(0.2<P≤0.4)、中等风险(0.4<P≤0.6)、较高风险(0.6<P≤0.8)与高风险(0.8<P≤1)5类. Cu元素低风险区位于金城镇、辛庄镇南与蚕庄镇北部，其他区域为无风险区. Zn与Cu的低风险区分布相似，集中于研究区中部，但Cu的低风险区面积较Zn更大；Zn元素存在一定的中等风险，表现在金城镇东，面积仅为532 hm²，占总面积的2%. Cr在研究区内基本不存在环境风险，仅有一处低风险区位于金城镇东部. Pb元素在金城镇西北部存在高风险等级，更有中等风险以上的区域分布于金城镇与研究区东部，风险概率小于0.4的区域则集中分布于研究区中部条带. Pb的风险概率高且中高等级环境风险区域在研究区内广泛分布，应引起警惕.

通过指示克里金法得到的风险概率结果与热点分析结果基本一致. 热点及周边地区一般存在低等级及以上的环境风险，而冷点区域均为无风险区. Pb元素的高风险区位置与热点所在位置有所差异，表现在辛庄镇东北部，这可能与热点分析所采用的空间距离有关. 因此，热点分析与指示克里金均能很好地反映研究区重金属的土壤环境风险，不同的是热点分析仅能评估点且必须为已知点数据，而指示克里金能够对未采样区域进行潜在污染风险评估，弥补了热点分析的不足.

环境风险分析结果与土壤重金属的空间分布状况也较为吻合，说明了潜在污染风险评估结果的可信性. 但Cr元素两者分析结果存在偏差，特别是辛庄镇东南部的岛状区域含量值介于108.5~153.4 mg·kg^-1之间，但表现为无风险等级. 这可能因在辛庄镇东南角仅有一个样点表现出异常高值，而该样点周围的样本含量均为低值且在阈值之下，使得该岛状区域及其邻接单元表现为高-低的空间集聚模式，在插值时被视为异常点而予以剔除. 同时，重金属的潜在污染风险区也与金矿的密集区呈现出密切的空间相关性，说明在研究区内，采矿等工业活动对土壤中Cu、Pb、Zn和Cr含量的现状与发展产生着不可忽视的影响.

(1) Pb元素含量介于17.1~50.3 mg·kg^-1之间，平均值与自然背景值相比呈现富集状态，说明Pb在研究区土壤中存在一定积累；Cu、Zn和Cr的含量表现相对贫化.

(2) Cu、Zn、Cr的最优变异模型为球状模型，Pb为指数模型. Cu和Pb属于中等空间自相关，两者含量可能受到自然和人为因素的共同影响；而Zn和Cr的块金系数小于0.25，表现为强烈的空间自相关，为自然来源.

(3) 4种重金属含量的空间分布基本一致，高值区位于金城镇、辛庄镇与蚕庄镇的交界地带，与金矿带位置吻合；而东南部丘陵区4种重金属含量总体较低，这主要与东南部地势相对高且没有金矿存在，受工业活动等因素影响小有关.

(4) 潜在污染风险评估结果显示Cu、Zn和Cr在研究区表层土壤中环境状况良好，不存在中高等级的环境风险；重金属Pb的超标概率高，风险概率在40%以上的区域分布广泛. 总体来说，金城镇、辛庄镇与蚕庄镇交界带的土壤环境质量较差，应引起足够的重视，而研究区南部可以认为是污染风险较小的环境安全区.