东北黑土区是我国重要的粮食生产基地, 随着工业化、城镇化和农业集约化的迅速发展, 黑土地土壤环境质量问题引起了科技工作者广泛关注^{[1, 2]}. 土壤重金属具有高生理毒性、可还原性、抗降解性和积累性等特点, 通过生物地球化学循环过程影响表层地球系统安全, 并经由食物链网威胁人类健康^[3~6]. 土壤重金属来源机制复杂多元, 既受到成土母质限制, 又有工农业生产等人类活动影响. 因此, 揭示区域土壤重金属积累现状与来源特征对于维护区域生态安全和保障人类健康至关重要^{[7, 8]}.

近年来, 学者们对土壤重金属污染开展了广泛研究. 在污染评价方法上, 常用的有单因子指数、内梅罗综合指数、地累积指数和潜在生态风险指数等方法, 可用于评价单一或综合土壤重金属的污染水平^[9]. 在污染来源解析层面, 常用的方法有多元统计模型、化学质量平衡模型、PMF模型和UNMIX模型等定性或定量方法^[10~13]. 徐祖奔等^[11]结合地累积指数和PMF模型研究发现黄河下游原阳滩区土壤重金属污染以农业源占主导. Fu等^[9]通过计算多个污染指数评价了海南岛东部海岸带土壤重金属的积累特征, 并结合地统计学和PMF模型解析出重金属污染为人类活动占主导的复合污染源. 总结诸多学者开展的相关研究, 发现土壤重金属污染主要为复合污染. 自然源主要包括成土母质来源, 即土壤形成过程中受风化和淋滤等作用影响, 会“继承”成土母质的地球化学特征^[2]. 人类活动来源如工业生产、化肥农药、交通工具尾气和电子垃圾等, 都会导致不同程度的土壤重金属污染^{[11, 14, 15]}. 此外, 土壤重金属污染及其来源存在显著的空间异质性. 本研究将结合多元统计分析模型与PMF模型对区域土壤重金属进行定性源识别与定量源解析, 筛选主要污染源, 进而针对性地防控重金属污染风险.

以往研究表明, 东北黑土区土壤环境质量总体较好, 但是由于近几十年工农业的快速发展, 导致黑土地出现了不同程度的重金属元素积累. 张淑珂等^[16]研究发现吉林白城市黑土区农田土壤重金属Cr、Cu和Cd存在部分超标区域. 李金凤等^[17]研究发现吉林省黑土区土壤主要受到Ni和As污染. 刘媚媚等^[18]揭示哈尔滨市阿城区海沟河小流域土壤重金属Cu、Zn、Hg、Ni和Cr均存在超标现象. 东北黑土区在保障国家粮食安全中具有举足轻重的地位, 其土壤环境风险防控非常重要^[19]. 松嫩平原西北缘莫旗地区是低山丘陵区黑土资源的主要分布区域, 大豆产量居全国第一, 是全国重要的商品粮生产基地. 本文依托1∶25万黑土地地表基质调查所获取的莫旗地区表层土壤重金属元素数据, 综合运用ArcGIS空间插值、多元统计分析方法和PMF模型, 开展黑土地土壤重金属污染特征研究, 定量解析低山丘陵区黑土地重金属来源, 揭示研究区土壤重金属空间分布格局, 评估土壤重金属污染风险水平, 以期为保障黑土地安全健康和黑土资源的可持续利用提供数据支撑.

莫旗地处大兴安岭与松嫩平原过渡地带, 介于北纬48°05′10″~49°50′50″、东经123°32′55″~125°16′14″之间, 面积10 358.94 km², 辖15个乡（镇）, 总人口约32.5万人. 地势由西北向东南倾斜, 海拔51~648 m, 以低山丘陵地貌为主. 莫旗属寒温带大陆性季风气候, 年均气温1.3 ℃, 年降水量在400~500 mm之间. 土壤类型丰富, 以暗棕壤、黑土和草甸土为主. 莫旗农牧业发达, 现有耕地面积56.94万hm², 占比91.76%, 大豆产量居全国首位；现有林地22.05万hm², 草地6.27万hm²（图 1）. 境内河流较多, 以嫩江水系为主, 水资源丰富. 莫旗地区地质背景以中生代花岗岩、碱性花岗岩、玄武岩、流纹岩和火山碎屑岩为主, 地质构造活动较为活跃. 中生代岩浆岩风化壳是最主要的成土母质, 另外现代河流阶地还发育有第四系堆积物.

图 1

Fig. 1

图 1 研究区土地利用类型与采样点分布示意 Fig. 1 Land use type and distribution of sampling sites in the research area

本研究依托1∶25万地表基质调查项目, 采用地质背景分区、地貌单元分类和采样深度分层次的方法部署调查点. 使用洛阳铲于研究区共采集0~20 cm表层土壤样品161件, 取样时去除地表腐殖质层, 剔除动植物残体、石块等杂物. 样品处理严格按照《多目标区域地球化学调查规范（1∶250 000）》中相关要求执行, 将样品在通风避光的室内自然风干并压碎, 干燥后用10目尼龙筛过筛, 然后将样品放在干净的塑料薄膜上, 混合均匀并铺成四方形, 划分对角线分成4份, 选择对角的两份样品送实验室分析. 土壤样品测试分析由华北有色（三河）燕郊中心实验室有限公司完成, Cu、Pb、Ni和Cd这4种元素测试方法为等离子体质谱法, 检出限分别为1 μg·g^-1、2 μg·g^-1、2 μg·g^-1和30 ng·g^-1；Cr元素测试方法为X射线荧光光谱法, 检出限为5 μg·g^-1；Zn元素测试方法为等离子体光学发射光谱法, 检出限为2 μg·g^-1；As元素测试方法为氢化物-原子荧光光谱法, 检出限为1 μg·g^-1；Hg元素测试方法为冷蒸气-原子荧光光谱法, 检出限为0.5 ng·g^-1. 抽取5%样品进行重复分析, 其相对偏差合格率均大于90%, 说明分析测试质量达到规范要求^[2].

单因子污染指数法（PI）是评价土壤污染程度或环境质量的无量纲指标, 是评估土壤重金属污染水平和估算人类活动影响的一种简单有效的方法^[20]. 采用1980年Hakanson^[21]提出的潜在生态风险指数法评价土壤重金属的潜在生态风险, 该方法综合考虑重金属含量、协同作用、生理毒性和环境对重金属的敏感性等因素, 被广泛应用于评价土壤生态危害^{[22, 23]}, 包括单元素生态风险评价（EI）和总生态风险评价（RI）. 两种污染指数计算公式与分级标准如下（表 1）：

表 1 (Table 1)

表 1 土壤重金属污染指数分级评价 Table 1 Grading evaluation of soil heavy metal pollution index PI 评价等级 EI 评价等级 RI 评价等级 PI < 0.7 清洁 EI < 40 轻微生态风险 RI < 110 轻微生态风险 0.7≤PI < 1 较清洁 40≤EI < 80 中等生态风险 110≤RI < 220 中等生态风险 1≤PI < 2 轻度污染 80≤EI < 160 较高生态风险 220≤RI < 440 较高生态风险 2≤PI < 3 中度污染 160≤EI < 320 高生态风险 RI≥440 高生态风险 PI≥3 重度污染 EI≥320 极高生态风险

表 1 土壤重金属污染指数分级评价 Table 1 Grading evaluation of soil heavy metal pollution index

式中, PI_i为土壤重金属元素i的单因子污染指数；C_i为元素i的实测含量数值（mg·kg^-1）；S_i为元素i的背景值, 本文选用松嫩平原土壤元素背景值^[24]. RI为总潜在生态风险指数；EI_i为单因素潜在生态风险指数；T_i为元素i的毒性响应系数, 表征该元素的环境危害程度, n为重金属种类数. 参考相关研究^[9], Cu、Pb、Ni、Cd、Cr、Zn、As和Hg的毒性响应系数分别为5、5、5、30、2、1、10和40.

主成分分析法（principal component analysis, PCA）可以用少数指标反映原始数据的信息, 被广泛运用于土壤重金属来源研究. 主成分分析通过将初始的因子载荷矩阵通过正交旋转变换进行“降维”处理, 将原始数据信息浓缩为较少的新变量, 再根据各变量的因子载荷情况来识别主因子^[25].

正定矩阵因子分解法（PMF）是美国国家环保署推荐使用的定量源解析方法. 该方法不依赖于污染源的化学成分谱分析, 而是将样品数据集当做一个矩阵, 将矩阵分解为源贡献率矩阵和源成分谱矩阵, 使用最小二乘法迭代运算, 不断分解原始矩阵以获得最优解^[26]. 基本计算公式如下：

式中, x_ij为样本i中j元素的含量；g_ik为第k个污染源对样本i的贡献率；f_kj为污染源k对j元素的贡献率；e_ij为残差矩阵.

莫旗地区161个表层土壤样品ω（Cu）、ω（Pb）、ω（Ni）、ω（Cd）、ω（Cr）、ω（Zn）、ω（As）和ω（Hg）的平均值分别为21.27、27.90、26.89、0.07、56.23、68.74、9.37和0.03 mg·kg^-1, 分别为松嫩平原土壤背景值的1.20、1.38、1.14、1.00、1.32、1.32、1.03和1.00倍（表 2）. 从全区看, 除Cd和Hg外, 其余重金属在研究区表层土壤中有着不同程度的富集. Cd和Hg的含量平均值虽未超过松嫩平原土壤背景值, 但其最大值分别是松嫩平原土壤背景值的2.43倍和3倍, 且两者的变异程度为强变异（CV≥0.36）, 反映Cd和Hg元素可能存在局部点位或小范围相对较严重污染. 其他重金属元素均为中等变异（0.15 < CV < 0.36）, 推测受到人类活动的影响有限. 与《土壤环境质量农用地土壤污染风险管控标准（试行）》（GB 15618-2018）规定的风险筛选值相比, 研究区8种重金属元素均不存在超标点, 土壤污染风险较低.

表 2 (Table 2)

表 2 研究区土壤重金属元素含量统计特征 Table 2 Statistical characteristics of heavy metal element content in the research area 项目 Cu Pb Ni Cd Cr Zn As Hg 最大值/mg·kg-1 34.90 55.30 50.30 0.17 103.00 112.00 20.30 0.09 最小值/mg·kg-1 6.77 14.00 5.37 0.02 19.10 35.50 3.28 0.01 平均值/mg·kg-1 21.27 27.90 26.89 0.07 56.23 68.74 9.37 0.03 中位数/mg·kg-1 22.10 26.80 26.60 0.07 55.60 66.60 9.37 0.03 标准偏差/mg·kg-1 5.11 5.81 7.35 0.03 14.74 13.04 3.23 0.01 变异系数 0.24 0.21 0.27 0.42 0.26 0.19 0.34 0.38 土壤污染风险筛选值/mg·kg-1 50 90 70 0.30 150 200 40 1.80 松嫩平原土壤背景值/mg·kg-1 17.78 20.23 23.65 0.07 42.46 52.05 9.14 0.03

表 2 研究区土壤重金属元素含量统计特征 Table 2 Statistical characteristics of heavy metal element content in the research area

采用克里格插值法对研究区8种土壤重金属的空间分布进行可视化分析, 结果显示Cu、Ni和Cr的空间分布相似（图 2）, 由东向西呈现较明显的梯度变化, 高值区主要集中在研究区中西部. 与研究区土地利用图对比发现, Cu、Ni和Cr的空间分布与土地利用分布特征吻合相对较好, 初步推测其含量与土地利用情况紧密相关. Pb和As的空间分布特征整体较为相似, 高值区主要分布在研究区东南部和东北部交通干道交汇处, 且该地区也是农业活动活跃区域, 故推测交通路网机动车辆尾气沉降和农业生产活动可能对Pb和As的积累有影响. Cd和Zn的空间分布特征整体较为相似, 高值区分布于研究区东北部和中西部地区, 为农业活动活跃地区, 初步推测耕地区农业种植活动对两者的积累有影响. Hg在研究区中部及北部区域含量较高, 极大值区零星分布, 与Hg元素变异程度较高有关. 总体上, 8种土壤重金属的积累量在研究区西部相对较高, 东南部相对较低.

图 2

Fig. 2

图 2 重金属元素空间分布 Fig. 2 Spatial distribution of heavy metal elements

为进一步确定土壤中8种重金属的污染程度, 以松嫩平原土壤重金属元素背景值为依据, 计算了单因子污染指数（PI）, PI平均值依次为：Pb（1.38） > Zn（1.32）=Cr（1.32） > Cu（1.20） > Ni（1.14） > As（1.03） > Cd（1.02） > Hg（1.01）. Pb、Zn和Cr三者PI值相对较高[图 3（a）], 为研究区主要土壤重金属污染元素. 按照污染等级分级标准, 研究区土壤重金属整体污染程度为轻度污染等级, Cd、Pb、Cr和Hg分别有4.35%、3.11%、3.11%和3.11%的调查点位为中度污染等级, 无重度污染等级调查点.

图 3

Fig. 3

图 3 重金属污染指数箱线图 Fig. 3 Box plot of heavy metal pollution index

采用Hakanson^[21]提出的重金属毒性反应系数作为标准, 评价了研究区土壤重金属的潜在生态风险. 研究区所有采样点显示Cu、Pb、Ni、Cr、Zn和As的生态风险较低[图 3（b）], 未超出生态风险筛选值, 仅Cd和Hg存在超出生态风险筛选值的采样点, 其原因在于Cd和Hg对生物具有相对较高的生理毒性. Cd有17.39%的采样点显示中等生态风险, Hg有36.65%和3.11%的采样点表现为中等和较高生态风险, 说明Cd和Hg在大部分土壤中的积累是安全的, 生态风险较低. 经计算, 研究区综合生态风险指数（RI）平均值为103.85, 68.32%的采样点未超出生态风险筛选值, 31.68%的采样点表现为中等生态风险. 因此, 研究区土壤重金属整体生态风险轻微, 但Cd和Hg因其高生理毒性导致局部出现的生态危害不容忽视^[27].

主成分分析的KMO值为0.608, 高于临界值0.5, 且Bartlett检验显示P < 0.001, 说明主成分分析结果可靠, 适合用于重金属的来源识别与判断. 通过最大方差法旋转矩阵使得各成分的系数载荷更为清楚, 通常将特征值超过1的成分作为主要成分. 主成分分析提取出4个特征值分别为2.96、1.56、1.51和1.04的主成分, 累计方差贡献率达88.30%, 在很大程度上可以解释研究区土壤重金属的污染来源.

第一主成分（PC1）方差贡献率为36.95%, 其中Cu、Ni和Cr具有较高的正载荷, 其权重系数分别为0.92、0.95和0.94, 且三者间呈极显著正相关关系（表 3和表 4）, 表现出高度的同源性. 三者超过背景值的采样点数分别为78.26%、73.29%和58.39%, 变异系数在0.24~0.27之间, 表明三者受人类活动的影响较大. 有研究表明, 在自然界中, Cr和Ni受成土母质的影响相对较大, 主要集中于超基性-基性岩中^{[28, 29]}. 本研究结果与之不同, 对比图 1和图 2发现, Cu、Ni和Cr在研究区耕地积累值较高, 表明表层土壤Cu、Ni和Cr的积累可能更多地受到了农业活动影响. Guan等^[30]和Cai等^[31]研究发现, Cu经常作为杀虫剂的有效成分添加到农药中, 如波尔多混合剂, 通过喷洒农药进入土壤. 此外, 有研究发现, 有机肥畜禽粪原料中也有部分重金属积累, 如Cu和Cr等^[32~34]. 因此, PC1可识别为以农药和肥料为主的农业活动源.

表 3 (Table 3)

表 3 重金属元素相关性分析1） Table 3 Correlation analysis of heavy metal elements Cu Pb Ni Cd Cr Zn As Pb 0.08 Ni 0.87** 0.08 Cd 0.11 0.11 0.05 Cr 0.80** 0.04 0.87** -0.05 Zn 0.50** 0.03 0.43** 0.51** 0.46** As 0.33** 0.52** 0.34** -0.21** 0.25** -0.11 Hg 0.26** 0.19* 0.05 0.40** 0.09 0.26** 0.07 1）*表示P < 0.05, **表示P < 0.01

表 3 重金属元素相关性分析1） Table 3 Correlation analysis of heavy metal elements

表 4 (Table 4)

表 4 土壤重金属主成分分析1） Table 4 Principal component analysis of heavy metals in soil 重金属 旋转后因子载荷 变量共同度 PC1 PC2 PC3 PC4 Cu 0.92 0.11 0.10 0.19 0.90 Pb -0.05 0.17 0.92 0.04 0.89 Ni 0.95 0.08 0.12 -0.04 0.92 Cd -0.07 0.88 0.02 0.26 0.84 Cr 0.94 0.03 0.03 -0.01 0.88 Zn 0.49 0.76 -0.08 0.03 0.82 As 0.31 -0.33 0.79 0.09 0.83 Hg 0.07 0.22 0.09 0.96 0.98 特征值 2.96 1.56 1.51 1.04 解释方差/% 36.95 19.48 18.90 12.98 累计解释方差/% 36.95 56.43 75.32 88.30 1）黑体字表示该主成分轴上的主要载荷元素

表 4 土壤重金属主成分分析1） Table 4 Principal component analysis of heavy metals in soil

第二主成分（PC2）中, Cd和Zn具有较高的载荷值, 权重系数分别为0.88和0.76, 累计方差贡献率为56.43%. 有研究表明, 化肥和农药的过度使用是土壤中Cd和Zn元素的重要来源之一^[35~37]. 特别的, 在无工业活动影响下, Cd被认为是农业活动的标志元素^[25]. 莫旗地区农林业发达, 工业化水平较低, 因而一定程度上可以排除工业生产来源影响. 同时, 近几十年来莫旗地区农业集约化程度不断增加, 土地利用强度增大, 农用化学品使用强度增加可能是土壤Cd和Zn相对富集的主要因素. 因此, 认为PC2可归纳为以农药和化肥为主的农业活动源.

第三主成分（PC3）方差贡献率为18.90%, 累计方差贡献率为75.32%. 其中Pb和As具有较高的载荷值, 其权重系数分别为0.92和0.79, 且两者间呈极显著正相关关系, 说明Pb和As可能具有相似的来源. 结合其空间分布状况发现, 交通干道与农业用地交汇区域是Pb和As的高值区, 故推测其可能受到交通与农业活动的共同控制. 研究表明, 汽车尾气中含有Pb和As等有害元素, 这些有害物质会通过大气沉降在周围土壤富集^[38~40]. 其中, Pb常被作为交通工具尾气污染的标志性元素. 莫旗地区大规模使用以柴油为主要燃料的农用机械, 其排放尾气与轮胎磨损均对土壤Pb富集有重要影响^{[41, 42]}. 此外, 农药与化肥中含有一定量的Pb和As元素, 如砷酸钙和砷酸钠常被用于除草剂和杀虫剂^[43]. 因此, 长期大量施用化肥、粪肥、除草剂和杀虫剂可能会造成农田土壤Pb和As污染. 综上, PC3可识别为交通与农业混合源.

第四主成分（PC4）方差贡献率为12.98%, 累计方差贡献率为88.30%. 其中仅Hg具有较高的载荷值, 权重系数为0.96. 研究区Hg的平均含量是背景值的1.01倍, 表明Hg元素积累的富集程度较低, 仅有少数调查点呈现出极高值. 结合本区功能定位来看, 研究区为农业生产功能区, 没有成规模的制造业与工业. 因此, 推测PC4为自然母质源.

本文在采用主成分分析判定污染源类别和主次关系的基础上, 通过PMF模型对研究区8种重金属来源进行定量解析, 并计算各源贡献率. 本文数据信噪比均大于3, 适合使用该模型. 本文PMF模型的最佳因子数为4, 各重金属的污染源成分谱情况见图 4.

图 4

Fig. 4

图 4 土壤重金属污染源成分谱 Fig. 4 Composition spectrum of soil heavy metal pollution sources

主成分分析表明, 研究区8种重金属存在4种污染源. 其中, Cu、Cr和Ni污染源主要为农业活动, 而污染源成分谱表明, 因子1对Cu、Cr和Ni的贡献率较其他因子大, 分别为50.06%、58.01%和58.77%, 因此, 因子1判定为农业活动源. 因子2对Pb和As的贡献率较高, 分别为55.09%和78.64%, 结合源识别结果判定因子2为交通与农业活动混合源. 因子3对Cd和Zn的贡献率较大, 分别为47.24%和40.05%, 故判定因子3为农业活动源. Hg污染来源主要为成土母质, 成分谱显示Hg污染来源中因子4的贡献率最大, 为73.36%, 故判断因子4为自然母质源.

从4种污染源对整个研究区土壤重金属污染的贡献率看（图 5）, 自然母质源仅占20.3%, 人类活动影响占比为79.7%, 整体上, 研究区人类活动源贡献率远大于自然源. 本研究认为, 研究区重金属元素来源复杂多样, 属于混合污染. 从图 2可以看出, 多数重金属元素在林地类型下含量更接近其背景值, 推测认为成土母质决定表层土壤重金属的初始分布与含量, 而长期农业种植活动会显著改变表层土壤重金属的积累, 即农化产品以及农业机械的使用已成为研究区表层土壤重金属污染的主要来源^{[44, 45]}.

图 5

Fig. 5

图 5 不同污染来源的总贡献率 Fig. 5 Contribution rates of different pollution sources

本研究结果基本符合研究区农业生产背景, 查明了研究区表层土壤8种重金属元素的积累情况, 但是尚未揭示重金属在土壤中的存在状态、迁移机制及其生物有效性等. 下一步应重点研究土壤重金属的赋存形态及其表生地质迁移过程, 结合成土母质、交通工具尾气、农化产品、煤炭等的重金属赋存特点进行更详细地分析, 更为精准地解析重金属来源及其迁移机制, 以便为防控黑土地土壤重金属污染和保障黑土地健康提供更切实的依据.

（1）莫旗地区表层土壤Cu、Pb、Ni、Cd、Cr、Zn、As和Hg重金属的含量平均值分别为21.27、27.90、26.89、0.07、56.23、68.74、9.37和0.03 mg·kg^-1. 与松嫩平原土壤背景值相比, 莫旗地区土壤重金属整体处于轻度污染水平. 相对而言, Pb（PI=1.38）、Zn（PI=1.32）和Cr（PI=1.32）为研究区主要土壤重金属污染元素. 与国家土壤环境质量标准相比, 8种重金属均未超出污染风险筛选值.

（2）潜在生态风险指数表明, 莫旗地区土壤重金属整体潜在生态风险（31.68%）轻微, 仅Cd和Hg分别有17.39%和39.76%的概率发生中等以上生态风险, 主要由其高生理毒性所导致. Cd和Hg变异系数较高, 在研究区东北及中部地区相对较高, 因此应注意防控其局部生态风险.

（3）定性与定量源解析识别归纳出3种污染源, 分别为农业活动源、交通与农业活动混合源和自然母质源. 自然源相对贡献率仅占20.3%, 远低于人为活动源79.7%. 其中, 研究区土壤Cu、Cr和Ni同源性显著, 主要来源于农业活动, 在耕地下积累较高；Pb和As主要受交通工具与农业活动的混合影响；Cd和Zn受农业活动的影响较高；Hg主要来源于自然母质.