CN108918815A - 一种土壤重金属风险预测方法 - Google Patents

Abstract

本发明提供一种土壤重金属风险预测方法，基于序贯条件模拟方法对目标区域内土壤的各种类重金属含量进行估计，并基于Hakanson潜在生态风险指数法综合各种重金属含量，生成可视化土壤重金属风险图。与传统的克里格插值法而言，序贯条件模拟方法更加能体现重金属分布的不确定性，克服传统克里格插值方法的平滑作用，具有更加良好的预测性；Hakanson潜在生态风险指数法可科学的综合各种土壤重金属含量，得出较为准确的评估风险；可较为直观的从可视化土壤重金属风险图中观测出高风险区域，并及时作出预防措施，具有良好的实用性。

Description

一种土壤重金属风险预测方法

技术领域

本发明属于地统计学领域，具体涉及到一种土壤重金属风险预测方法。

背景技术

随着城镇化、工矿业的快速发展以及农药有机肥的大量使用，越来越多的重金属通过污水灌溉、大气沉降、降水径流等途径进入农田土壤环境。农田中重金属的含量日益累积，而重金属难以被土壤中的微生物降解，重金属含量一旦超过特定阈值，将对农作物安全、生态环境和人类健康造成严重危害。针对农田土壤重金属的危害胁迫，预先识别农田重金属的风险区域，能对区域的农田环境保护、重金属污染预警和风险管控等提供科学的指导和依据。

现有方法在判断土壤重金属区域风险时一般分为两步，首先是对土壤重金属含量空间分布的预测，然后再根据评价方法确定风险的阈值，进而判断该地区是否有风险。

农田重金属风险识别需要区域内土壤重金属精确的空间分布制图作为支撑，通过空间分布及趋势来分析其聚集特点，对比预测含量与风险阈值来判断土壤是否存在风险性。当前在土壤属性预测研究中最常用的是克里格插值法，识别土壤重金属风险一般有两种途径：一是直接利用插值方法对土壤重金属含量进行空间预测，获得重金属含量的空间分布后和阈值进行对比，多用线性克里格模型，如普通克里格(Ordinary Kriging,OK)；二是利用不确定性评估来确定区域重金属含量超过阈值概率的空间分布，多用非线性克里格模型，如指示克里格(Indicator Kriging,IK)。

污染评价可以界定土壤受重金属污染的程度，通过分级确定重金属污染的阈值，从而划分土壤中重金属的危害等级。对重金属污染的评级方法和标准有很多，主要可分为指数法、模型指数法和其他评价方法。

现有常用的预测方法存在一定的限制和缺陷，一方面克里格插值是对未采样点的最优无偏估计，其对数据的分布有较为严格的要求。但土壤中重金属的含量不仅与土壤本身质地有关，还与人为因素的影响有关，这导致土壤重金属的分布呈现极高的空间异质性。在实际采样中也存在着异常值，使数据成偏态分布，具有较大的变异系数和偏度值。这使得数据经常不完全符合克里格的要求。一般通过非线性变换来使数据服从正态分布，但难以转换为原尺度；或者利用变异云图、三倍标准差准则等来剔除异常值，但异常值客观存在，对其直接忽视可能会对预测的结果产生较大的影响。同时克里格插值还存在平滑效应，平滑程度受到采样点分布影响，采样点相距越远其插值的平滑效应越强烈。而在大的尺度下，由于农田分布不均、空间、时间和人力的成本等的限制性，很难做到均匀密集的采样，产生的平滑效应可能导致土壤重金属含量异常区的重要信息丢失。

另一方面，指数法形式简单易懂，容易掌握与操作，但还是存在一定的问题:例如单因子指数法难以全面、综合地展现土壤的污染程度。内梅罗综合污染指数法，人为估算因子带有主观性，容易夸大高浓度重金属污染的影响。地累计指数法忽略了元素间的相关性，各金属的不同污染能力和生物有效性，难以比较元素间或区域间环境质量，在选取k值时会带有主观性。

发明内容

本发明提供了一种土壤重金属风险预测方法，基于序贯条件模拟方法对目标区域内土壤的各种类重金属含量进行估计，并基于Hakanson潜在生态风险指数法综合各种重金属含量，生成可视化土壤重金属风险图。与传统的克里格插值法而言，序贯条件模拟方法更加能体现重金属分布的不确定性，克服传统克里格插值方法的平滑作用，具有更加良好的预测性；Hakanson潜在生态风险指数法可科学的综合各种土壤重金属含量，得出较为准确的评估风险；可较为直观的从可视化土壤重金属风险图中观测出高风险区域，并及时作出预防措施，具有良好的实用性。

本发明提供了一种土壤重金属风险预测方法,所述方法包括以下步骤：

确定土壤重金属风险预测目标区域；

在所述目标区域内选择采样点并进行土壤采样；

测定所述采样点的土壤重金属含量；

基于所述采样点的土壤重金属含量数据，推导未采样点的土壤重金属含量数据并生成所述目标区域重金属含量的空间分布图；

基于所述空间分布图，计算所述目标区域的生态风险指数并生成土壤重金属风险评价图。

优选的实施方式，所述目标区域为同时存在农业生产和工业生产的地区。

优选的实施方式，基于梅花形布点法对所述采样点进行土壤采样。

优选的实施方式，所述土壤重金属包括铜、锌、铅、镉、铬、砷、汞。

优选的实施方式，基于序贯指示模拟方法，根据所述采样点的土壤重金属含量数据推导未采样点的土壤重金属含量数据。

优选的实施方式，所述序贯指示模拟方法包括以下步骤：

对所述采样点的重金属含量数据进行指示变换；

分别计算重金属含量数据在各阈值条件下的指示半变异函数；

基于所述指示半变异函数建立先验条件累积分布函数；

将所述目标区域划分为同一分辨率的格网，定义一条经过所有格网的随机路径，并在第一个位置格网处从条件累积分布函数中随机抽取一个值作为模拟值；

将所述模拟值用于下一位置格网的先验条件累积分布函数，从下一位置格网的先验条件累积分布函数中随机抽取一个值作为模拟值，重复执行该步骤直至所有格网模拟完毕。

优选的实施方式，所述计算所述目标区域的生态风险指数包括以下步骤：

计算所述目标区域土壤中单种重金属的潜在生态风险指数；

计算所述目标区域土壤中多种重金属的综合潜在生态风险指数。

优选的实施方式，所述单种金属的潜在生态风险指数计算公式为

为单项重金属风险系数，C_i为表层土壤重金属浓度实测值，为参比值。

优选的实施方式，所述多种重金属的综合潜在生态风险指数计算公式为

其中，为第i种重金属潜在生态危害系数，为第i种重金属的毒性响应系数；RI为综合潜在生态危害指数。

本发明提供了一种土壤重金属风险预测方法，基于序贯条件模拟方法对目标区域内土壤的各种类重金属含量进行估计，并基于Hakanson潜在生态风险指数法综合各种重金属含量，生成可视化土壤重金属风险图。与传统的克里格插值法而言，序贯条件模拟方法更加能体现重金属分布的不确定性，克服传统克里格插值方法的平滑作用，具有更加良好的预测性；Hakanson潜在生态风险指数法可科学的综合各种土壤重金属含量，得出较为准确的评估风险；可较为直观的从可视化土壤重金属风险图中观测出高风险区域，并及时作出预防措施，具有良好的实用性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。

图1示出了本发明实施例的土壤重金属风险预测方法流程图；

图2示出了本发明实施例梅花形布点法示意图；

图3示出了本发明实施例的采样点分布图；

图4示出了本发明实施的采样点土壤重金属描述性统计图；

图5示出了半变异函数图像示意图；

图6示出了各阈值条件下Cu的指示半变异函数图象；

图7示出了各阈值条件下指示半变异函数的最佳参数表；

图8示出了十分位下Cu的累积分布函数；

图9示出了本发明实施例序贯指示模拟的工作界面及结果示意图；

图10示出了本发明实施例的土壤重金属空间分布图；

图11示出了本发明实施例的土壤重金属潜在生态风险分级标准；

图12示出了本发明实施例的土壤重金属潜在风险指数；

图13示出了本发明实施例的SISIM与IK的误判率对比数据表。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

本发明实施例提供了一种土壤重金属风险预测方法，解决传统方法在重金属含量空间分布预测中存在的采样异常值筛除和平滑效应严重的缺陷，通过前期对数据的指示化处理，将异常值指示化进行指示半变异函数的计算，同时利用序贯指示模拟对研究区的格网化，降低平滑效应的影响；对预测结果采用考虑因素广泛的Hakanson潜在生态风险指数法进行风险区域识别并以可视化手段进行呈现，具有良好的准确性、直观性和实用性。

在统计学中，当变量在空间一些位置的值依赖于该变量在其他位置的值时，就存在着空间自相关。在地统中，有些变量，在考虑空间或时间时，变量往往不是随机的，因而，计算这些变量的特征时，除了变量的均值、方差等统计量，还需要计算变量的空间结构。依赖于空间分布位置的变量称为区域化变量。对于一个变量Z，如本发明实施例的某种重金属含量Z(i)，随不同空间位置i而变化，该变量的变化决定于三个部分：大尺度范围内的总体值f(i)，用于表示变化趋势、局部小范围内的空间依赖性s(i)和误差ε，可用公式Z(i)＝f(i)+s(i)+ε表示。

区域化变量的两大特点是随机性和结构性，随机性是指在局部的某一点而言，区域化变量的取值时随机的；结构性是指，对整个区域而言，存在一个总体或平均的结构，相邻的区域化变量的取值具有该结构所表达的相关关系。

例如，在小范围中，重金属i的分布与大区域是不同的，小范围内的重金属含量可能高于、低于大区域的均值，而这是随机的。具体的，该小范围可以用Z(i₁)表示某种状态，每个观测值可以认为是Z(i₁)的一个实现；同样，另外一个小范围可以用Z(i₂)表示，每个观测值可认为是Z(i₂)的一个实现；同时，Z(i₁)和Z(i₂)存在着某种数量关系。

具体的，对于随机过程，要完整地描述其特征，必须给出它的分布函数或分布密度，但是这种分布一般很难求得。所以研究的重点是其数学特征，这些数学特征主要包括：数学期望、相关函数、方差、协方差均方值；其中，数学期望是一阶矩，后面四个数字特征都是二阶矩。

当一个随机过程是二阶矩过程时，我们可以研究它的二阶矩数字特征，从而对随机过程进行相应的分析，比如判断是否平稳等。

区域化变量完全的分布函数或分布密度函数的求取也没有必要，区域化变量分布函数或分布密度函数的头两个矩，一阶矩和二阶矩足以提供大多数情况下所研究问题的近似解，如本发明实施例所涉及的重金属含量分布数据。另一方面，实地所测的数据也不足以求得区域化变量完全的分布函数或分布密度函数，所以区域化变量完全的分布函数或分布密度函数的求取是没有必要的。

此外，如果只使用随机函数的头两个矩，如果两个随机函数的一阶矩和二阶矩是相同的，那么可以认为这两个随机函数是相同的。

一阶矩为一个区域化变量的平均值函数，定义为E[Z(x)]；

一个区域化变量中的二阶矩共有三个，分别为区域化变量的方差函数、区域化变量的协方差函数、变差函数；其中，变差函数又称为半方差，为一个区域化变量两点差值方差的一半。

在实际应用中，以上函数往往通过若干测定值作出估计。例如，半方差函数需要通过数学期望E[Z(x)-Z(x+h)]的值进行计算，因此，必要有Z(x)和Z(x+h)的若干实现数据。

然而，在土壤空间研究中，在点x和点x+h往往只能得到一对数据，不可能在空间同一点再去得到第二对数据，因此，为了克服这个困难，需要对区域化变量提出一些假设及限定。

假定在局部分为内，变量的均值为一常数，不随位置而变化；假定样本点x和y的协方差C[Z(x),Z(y)]存在，且只取决于样本点x和y的距离，则Z是二阶平稳的，结果是Z的均值和协方差在空间上不发生变化。

这样，在空间某一局部范围内，对空间某一点x₀，相距为h的多个点，可以看做是Z(x₀)的多个实现，从而可以进行统计推断和估值预测。

进一步的，如果随机函数Z(x)具有以下内蕴假设：

均值存在且不取决于x，即E[Z(x)]＝m，对于任意E[Z(x+h)-Z(x)]＝0；

对于任意距离h，变量[Z(x+h)-Z(x)]具有一个有限的方差，该方差不取决于x；

则对于任意的x和h，下式成立：

Var[Z(x+h)-Z(x)]＝E{[Z(x+h)-Z(x)]²}＝2γ(h)，其中，γ(h)称为半变异函数，半变异函数用来表征随机变量的空间变异结构，或空间连续性。

因此，平稳性主要包括两类，其中一类指均值平稳，假设均值是不变的，并且与位置围观；另一类适于协方差函数有关的二阶平稳和与半变异函数有关的内蕴平稳；协方差平稳是假设具有相同的距离和方向的任意两点的协方差是相同的，协方差只与两点的值相关而与它们的位置无关；内蕴平稳是假设具有相同距离和方向的任意两点的方差是相同的。

本发明实施例提供的土壤重金属危害预测方法在模型预测阶段，主要利用采样点的数据和半变异函数的结构性，采用序贯指示克里格的方法，对未采样点的数据进行估计，最终通过Hakanson潜在生态风险指数法生成土壤重金属危害风险图，从而判断重金属风险区域。

图1示出了本发明实施例的土壤重金属风险预测方法流程图。本发明实施例提供了一种土壤重金属危害预测方法主要包括数据采集、模型预测和风险评估三个方面的步骤，具体实施步骤如下：

S101：确定土壤重金属危害预测的目标区域；

重金属污染指由重金属或其化合物造成的环境污染，主要由采矿、废气排放、污水灌溉和使用重金属超标制品等人为因素所致，因此，本发明实施例提供的土壤重金属危害风险预测方法应用的目标区域主要为同时存在农业生产和工业生产的地区。

考虑到重金属危害风险预测是用于为重金属污染的风险管控等提供依据，目标区域以行政区划来选择有利于成果的应用。本发明实施例的目标区域选择为广州市增城区，它位于广东省中东部，总面积1616.47km²，其地貌多以南部的三角洲与河谷平原为主，耕地土壤以赤红壤和渗育型水稻土为主。

增城区是珠江三角洲粮食等农产品主要生产基地，同时拥有零部件制造业、服装制造业等支柱产业。研究该区域的农田土壤重金属含量及其风险的空间分布情况，对于保护该地区的农产品安全有重要的参考意义。

S102：从所述目标区域内选择采样点并进行土壤采样；

确定土壤重金属危害预测目标区域后，需要收集目标区域内的土壤重金属含量数据。考虑到目标区域的范围大小，对目标区域内的土壤完全测定重金属含量是不现实的，因此，需要对目标区域内进行采样。

具体实施中，分为采样点的确定和对采样点进行土壤采样两个步骤，通过对土壤采样样品的测定获得采样点的重金属含量数据。

在本发明实施例中，首先对增城区的土地利用数据进行收集，根据研究区内农田的分布情况以及工矿企业的分布情况来设计采样点。采样点尽可能覆盖整个研究区，尤其是人类工业活动较为密集的场所，主要包括原材料开采区域、工厂集中区域、物流集中区域、废物处理区域等场所，结合现代生产的特点，大部分人类活动密集区域多以路网为基础的，或可以认为大部分人类活动密集区域的路网条件较为良好，因此，路网及路网周围的重金属含量较其余区域较高，具有一定的代表性，可以以目标区域内的路网作为采样点的干线，在路网附近进行采样。

在土壤采样实施过程中，样品采自表层土壤，表层土壤是指距地表高度0～20cm的土壤，采样方法选择梅花形布点法。

图2示出了本发明实施例梅花形布点法示意图。梅花形布点法是指在一点取完土壤样品后，基于该点向四周辐射10m选四个点进行采样，将五个点位的土壤样品制成混合样品，作为最初采样点的采样结果，并对其进行GPS坐标定位，确认其经纬度坐标。通过梅花形布点法采集的土壤样品可避免了局部极小区域的特异性，避免数据产生极端值。

图3示出了本发明实施例的采样点分布图。具体的，本发明实施例收集了2010年增城区的土地利用数据，分析其农田及工矿用地分布情况，同时考虑采样的可操作性、时间及人力成本等，最终决定以增城区路网为基础，在增城区共布置204个农用地采样点，采样所得的土壤样本数共204份。

S103：测定所述采样点的样品重金属含量；

土壤采样后需要对样品的重金属含量进行测定，具体实施中，对风干后的样品，过2mm的筛制成土壤样品，分别测定总铜、总锌、总铅、总镉、总铬、总砷、总汞含量，需要说明的是，由于本发明实施例的研究区内并无特殊矿产地、特殊矿产中转地或特殊矿产加工地等场所，因此，只需测定常见的重金属元素即可。针对不同矿产地或具有矿产中转、加工产业的区域，可针对特定的重金属元素进行专门的检测，本发明实施例不一一介绍其测定方法。

进一步的，砷As和汞Hg采用1:1的HC_l-HNO₃制备待测液，由还原气化—原子荧光光度计法测定；铜Cu、锌Zn、铅Pb、镉Cd、铬Cr采用HF-HNO₃-HClO₄消化法制备待测液，除Cd、Pb采用石墨炉原子吸收法测定外，其余采用火焰原子吸收分光光度法测测定；具体的，原子吸收光谱仪(石墨炉)可采用日本HITACHI公司的Z-2000，火焰原子吸收光谱仪可采用日本HITACHI公司的Z-5300；检测时，加入国家土壤标准物质GSS-1用于质量控制。

经测定，各重金属回收率在92％～104％范围内，相对标准偏差在允许范围±10％之内，成分分析结果较为可靠，可供后续步骤采用。

需要进行说明的是，最终得出每组采样点的数据至少包括采样点坐标和总铜、总锌、总铅、总镉、总铬、总砷、总汞的重金属含量，重金属含量的单位一般为mg·kg-¹。

图4示出了本发明实施的采样点土壤重金属描述性统计图。对采样点的数据进行初步统计分析，本研究对7种重金属Cu、Zn、Pb、Cd、Cr、As、Hg分别利用SPSS统计软件统计计算了最大值、最小值、均值、标准差、变异系数、偏度、峰度和K-S检验值。对增城区农田土壤7种重金属元素进行描述性统计分析。Cu、Zn、Pb、Cd、Cr、As、Hg的平均含量分别为12.71mg·kg^-1、39.63mg·kg^-1、42.56mg·kg^-1、0.08mg·kg^-1、35.44mg·kg^-1、8.84mg·kg^-1和0.11mg·kg^-1。

其中，变异系数的大小可粗略估计变量的变异程度，能在一定程度上反映样品受人为影响的程度，变异系数小于10％为弱变异性，变异系数在10％—100％之间为中等变异性，变异系数大于100％为强变异性。

偏度反映了地理数据分布的对称性，其不为0时，表示数据存在极端值。

峰度反映了地理数据在均值附件的集中程度。

K-S检验是一种分布拟合优度的检验，用来检验地理数据是否符合正态分布的方法。

由SPSS统计软件进行描述性统计可得，除Zn、Pb外其他重金属的分布都存在一定的右偏及异常值。研究区内土壤的7种重金属元素的变异系数在29％-87％之间，属于中等变异。重金属Cu、Cd、As和Hg的变异系数相对较大，说明这4种重金属的分布受人为因素影响较大。研究区土壤环境的pH平均为5.4，通过比较《土壤重金属风险评价筛选值珠江三角洲》(DB 44/T1415-2014)与《土壤环境质量标准》(GB15618-1995)发现，7种重金属元素都未超过土壤环境质量二级标准，但均略高于土壤背景值，说明研究区农田土壤重金属对作物和生态环境存在潜在的风险性。

以上步骤为数据采集阶段，通过数据采集后，获得目标区域内采样点的重金属含量数据。接下来执行的步骤为有关目标区域模型预测的步骤，以已获取的采样点重金属含量数据，对目标区域内未采样点的重金属含量进行估计从而做出风险预测。

本发明实施例采用空间模拟中的序贯指示模拟方法来进行重金属的空间分布预测，它综合了序贯随机模拟算法与指示克里格方法的各自的优点，在随机模拟方法中适应性较强。

序贯指示模拟方法的基本思路为首先将目标区划分为统一分辨率的格网，对未采样的格网利用指示克里格方法获得先验信息，而指示克里格不参与模型随机模拟计算，根据样本数据的条件累积分布函数(CCDF)的条件概率分布，构建后验概率模型进行空间模拟，最后从该分布中随机取值作为模拟现实得出模拟值，利用序贯算法一直模拟到最后一个格网，通过进行给定的n次模拟实现，得到不确定性评估的结果。

具体的实施步骤如下：

S104：对所述采样点的重金属含量数据进行指示变换；

具体的，根据指示克里格的思想得到条件累积分布函数(CCDF)这个先验信息。

指示变换步骤如下：设{Z(x_i)，i＝1，2，...，n}是一组采样数据，给定K个阈值，通常采用分位数来确定，分位数越多其CCDF的重建越精确，编码成0和1的指示变量，其公式为：

具体的，本发明实施例选取十分位数，即设定9个阈值，对采样点的重金属含量数据进行指示变换，利用SPSS统计软件取各重金属含量0.1～0.9的分位数作为阈值。

具体的，本发明实施例统计采样点的重金属铜的采样数据，其最大值和最小值分别为37.16mg·kg^-1和2.61mg·kg^-1,本发明实施例取k＝9个阈值，用分位数为0.1、0.2、0.3、0.4、0.5、0.6、0.7、0.8、0.9，则将2.61至37.16之间的范围区间按分位数取阈值，分别为Z_k＝1＝6.30、Z_k＝2＝9.73、Z_k＝3＝13.21、Z_k＝4＝16.60、Z_k＝5＝20.07、Z_k＝6＝23.40、Z_k＝7＝26.94、Z_k＝8＝30.32、Z_k＝9＝33.71，单位为mg·kg^-1；需要注意的是，目标区域中同一次执行指示变换公式时，阈值应保持不变；针对不同的重金属种类，分位数不变，阈值取值空间改变，同一分位数对应的阈值也会发生改变。

确认阈值后，即可进行指示变换，公式为：

具体实施中，还可用条件概率来描述指示函数：

当X_a，a＝1,2,3…,n表示采样点，本发明实施例n＝204，对于采样点，

I(X_a；Z)＝P{Z(X_a)≤Z|Z(X_a)＝Z_a}

此时，某待估测点X的指示函数估计值可以表示为：

I^*(X；Z)＝P{Z(X)≤Z|Z(X_a)＝Z_a,a＝1,2,L,n}

对于采样点来说,指示值可解释为已知该点的实测值为Za时,该点的真实值小于等于阈值的概率,而对于待估测点,其指示函数估计值可解释为已知待估点周围信息，即采样样本的实测值时,该点的真实值小于等于阈值的概率。

具体的，可通过采用交叉验证来检验模拟精度，采用2:8的比例随机抽样，选取40个样点作为独立验证点不参与建模，其余164个点作为训练样本。

S105：分别计算重金属含量数据在各阈值条件下的指示半变异函数；

图5示出了半变异函数示意图。半变异函数是区域化变量空间变异性的一种度量，反映了空间变异程度随距离而变化的特征，从而可定量的描述区域化变量的空间相关性。

半变异函数γ(h)的数学表达式为：

其中，x为目标区域中其中一点，x+h为目标区域中距离x距离为h的一点；Z(x)为x除的取值，即本发明实施例的某种重金属含量数据，Z(x+h)为x+h处的取值，E为一常量。一般的，h命名为步长，在同一方向上，对不同步长hi(i＝1，2…，n)的半变异函数进行统计，可得到一组不同的实验半变异函数值γ(h_i)。以h为横坐标，γ(h_i)为纵坐标所得到的一组{h,γ(h_i)}点称为半变异函数图。

半变异函数图中的几个主要参数分别为a、co、c。其中，a表示变程(range)，反映区域化变量在空间上具有相关性的范围，在变程范围之内数据具有相关性，在变程范围之外数据互不相关，针对不同的模型而有所不同，如对于球状和线性模型，A表明土壤性质存在空间变异结构的最大相关距离；对于高斯模型和指数模型最大相关距离则分别为1.73A和3A。变程大小受观测尺度限定，在变程范围内，样点间的距离越小，其相似性，即空间相关性越大。当h>A时，区域化变量z(x)的空间相关性不存在，即当某点与已知点的距离大于变程时，该点数据不能用于内插或外推。

co表示块金效应(nugget effect)，用以描述区域化变量在很小的距离内发生的突变程度。理论上，当采样点间的距离为0时，半变异函数值应为0.由于存在测量误差和空间变异，使得采样点非常接近时，他们的半变异函数值不为0，即存在块金值，块金效应是非空间性质观测值的变异。在数学上，块金效应相当于变量的纯随机部分。

c为基台值(sill)，当h趋于核实的间距(A)时，半方差的极限值，等于变量的方差，超过此变程可以认为属性变量空间独立。反映变量在空间上的总变异性大小，基台值越大说明数据的波动程度越大，参数变化的幅度越大。

块金与基台的比值叫基底效应，此值越大，说明空间变异更多的是随机成分引起的，否则则是由特定的地理过程或多个过程综合引起的

地统计学将变异函数的理论模型分为三大类：一类是有基台值模型，包括球状模型、指数模型、高斯模型、线性有基台值模型和纯块金效应模型；一类是无基台模型，包括幂函数模型、线性无基台值模型、抛物线模型；还有一类为孔穴效应模型；具体的函数图可参照现有技术进行参考。

图6示出了各阈值条件下Cu的指示半变异函数图象，图7示出了各阈值条件下指示半变异函数的最佳参数表，图8示出了十分位下Cu的累积分布函数。半变异函数的计算需要分别设定步长、步长数、步长容差以及搜索方向和方向容差。步长为采样点的空间间隔距离，步长数为最大变程a下包含的步长数量，步长容差为允许的搜索半径误差；步长和步长数确定如何对半半变异函数值进行分组。由于采样点分布是非均匀的，采用平均最近邻方法确定步长总数，即测量每个点位与其最近邻点位置之间的距离，然后计算所有这些最近邻距离的平均值。根据点对间最大距离的一半除去步长的经验法则得到步长数，一般取步长的一半作为步长容差。搜索方向是指针对空间数据插值存在的各向异性而设定的搜索方向，一般用E-W方向、S-N方向、SW-NE方向以及SE-NW方向来进行搜索，方向容差是指允许的搜索方向误差。

可得步长约为2500m，步长数为10，步长容差为1250m。分别计算0°、45°、90°和135°方向以及各向同性下的指示半变异函数，方向容差为22.5°，并调整所得到的指示半变异函数参数采用合适模型使拟合结果最佳，利用GS+9.0计算小于各阈值条件的指示克里格估计，重建CCDF的F_z(x)。

S106：将目标区域划分为同一分辨率的格网，定义一条经过所有格网节点的随机路径并在第一个位置格网处从F_z(x)中随机抽取一个值作为模拟值；

条件累积分布函数在目标格网节点的估计值依赖指示半变异函数信息的简单克里格方法得到：

λ_i为权重，表示z_k阈值下样品的期望值，权重通过简单克里格方法确定：

γ^I(x₀，x_i；z_k)和γ^I(x₀，x_i；z_k)分别是z_k阈值处点x_i和点x_j的之间指示半变异值和点x₀和点x_j的之间指示半变异值。

S107：将所述模拟值加入到条件数据集中，在新的数据集的条件下，加入到下一个位置的格网先验条件累积分布函数Fz(x)的建模中，在格网处继续从条件分布函数中抽取一个值，重复该过程到所有格网节点被模拟完，即完成一次模拟实现。

S108：重复执行步骤S106和步骤S107，直至模拟次数达到预设值；

根据每个格网统计n次中大于给定阈值的次数n1，n1与n的比值就是该格网在n次模拟中大于阈值的概率。

图9示出了本发明实施例序贯指示模拟的工作界面及结果示意图。

本发明实施例考虑研究区的大小，面积大约为60km×50km，以及栅格划分计算和可视化运算时间，将研究区划分为600×500的格网，每个格网栅格的大小为100m×100m。通过SGeMSv2.5b软件中SISIM模块，设定一条经过所有栅格的随机路径，利用指示半变异函数参数输入，包括指示阈值、拟合模型、变程a等，对栅格进行1200次模拟实现。

S109：生成目标区域的重金属含量的空间分布图；

图10示出了本发明实施例的土壤重金属空间分布图。通过后处理中的E-type估计得到栅格预测的重金属含量值，利用ASCII文件导入Arcgis10.2转为栅格数据，得到重金属含量的空间分布情况。

通过步骤S105～S1109，可得到目标区域的各种重金属空间分布估测数据，以下需要进行的为风险评估相关步骤。风险评价可以根据重金属对生态环境、人类健康等评价对象的毒性机理界定土壤受重金属危害的程度，通过分级确定重金属风险阈值，从而划分土壤中重金属的危害等级。

Hakanson潜在生态风险指数法是一种有效划分出重金属潜在的危害程度的方法，对于非典型研究区而言，它能考虑重金属的环境、生态等效应、迁移规律以及毒理学原理；它既关注了单种重金属的危害性，又考虑了多种重金属的综合危害；结合研究区背景值，弱化了区域差异的影响，适合于区域尺度的农田土壤重金属风险程度研究。

S110：目标区域生态风险指数计算；

根据Hakanson指数法，研究区域内土壤中单种重金属的潜在生态风险指数以及土壤中多种重金属的综合潜在生态风险指数可分别用和RI表示，计算公式为：

式中：为单项重金属风险系数，C_i为表层土壤重金属浓度实测值，为参比值(一般为土壤背景值或二级标准)；为第i种重金属潜在生态危害系数，为第i种重金属的毒性响应系数；RI为综合潜在生态危害指数。

图11示出了本发明实施例的土壤重金属潜在生态风险分级标准，图12示出了本发明实施例的土壤重金属潜在风险指数。对研究区204个样点分别计算单种重金属生态危害指数，土壤参比值采用珠江三角洲土壤重金属元素背景值上限，毒性系数采用已有研究计算的重金属毒性系数(Hg＝40，Cd＝30，As＝10，Cu＝5，Pb＝5，Cr＝2，Zn＝1)为评价依据。观察结果可发现，本发明实施例的风险指数Er：Hg>Cd>As>Pb>Cu>Cr>Zn，Hg和Cd两种重金属元素存在中等以上潜在生态风险的可能。

S111：生成土壤重金属风险评价图。

图12示出了本发明实施例的土壤重金属风险综合评价图。利用序贯指示模拟预测的农田土壤重金属结果的栅格格网值，计算每个格网的Hg、Cd以及7种重金属的综合风险指数值，得到研究区土壤重金属风险评价图，以此作为农田土壤重金属风险区识别的依据。

具体的，由于通过模拟生成的格网数据和采样点数据可能存在偏差，导致最终的模拟预测结果可能出现误判，可以使用误判率(Error Rate,ER)来评价识别风险的准确度。其误判情况一般分为两种，一是将非风险区判定为风险区，二是将某一级别风险区判定为非风险区或其他级别风险区。

式中：N₁为将非风险区判定为风险区的样点数，N₂为将某一级别风险区判定为非风险区或其他级别风险区的样点数，N为总样点数。

图13示出了本发明实施例的SISIM与IK的误判率对比数据表。通过将41个验证样点数据分布，与序贯指示模拟(SISIM)的模拟结果对比，找出分类错误的点，分别计算它们的误判率(ER)。

由序贯指示模拟(SISIM)计算出的Hg潜在生态风险指数、Cd潜在生态风险指数、综合潜在生态风险指数介于4.88-17.07％之间，较优于指示克里格(IK)的9.76％-19.51％，其精度属于可接受的范围。

本发明实施例提供了一种土壤重金属风险预测方法，采用序贯指示模拟方法与GIS相结合，预测了农田土壤重金属含量的空间分布，和常规的方法相比，普通克里格插值中预测点的结果受到邻域的采样点和到采样点的距离影响，同时产生平滑效应。而序贯指示模拟的每次实现只是利用了采样点的条件概率分布的信息，平滑效应低，没有忽略异常值对含量高值区域的贡献，这使得它的模拟结果更加接近真实情况；根据Hakanson潜在生态风险指数法划定了农田土壤重金属风险区域并分级，对农田土壤重金属含量的空间分布进行了可视化呈现，清晰直观的观测到土壤重金属风险区域，具有良好的实用性。

以上对本发明实施例所提供的土壤重金属风险预测方法进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (9)

1.一种土壤重金属风险预测方法,其特征在于，所述方法包括以下步骤：

确定土壤重金属风险预测目标区域；

在所述目标区域内选择采样点并进行土壤采样；

测定所述采样点的土壤重金属含量；

基于所述采样点的土壤重金属含量数据，推导未采样点的土壤重金属含量数据并生成所述目标区域重金属含量的空间分布图；

基于所述空间分布图，计算所述目标区域的生态风险指数并生成土壤重金属风险评价图。

2.如权利要求1所述的土壤重金属风险预测方法,其特征在于，所述目标区域为同时存在农业生产和工业生产的地区。

3.如权利要求1所述的土壤重金属风险预测方法,其特征在于，基于梅花形布点法对所述采样点进行土壤采样。

4.如权利要求1所述的土壤重金属风险预测方法,其特征在于，所述土壤重金属包括铜、锌、铅、镉、铬、砷、汞。

5.如权利要求1所述的土壤重金属风险预测方法,其特征在于，基于序贯指示模拟方法，根据所述采样点的土壤重金属含量数据推导未采样点的土壤重金属含量数据。

6.如权利要求5所述的土壤重金属风险预测方法，其特征在于，所述序贯指示模拟方法包括以下步骤：

对所述采样点的重金属含量数据进行指示变换；

分别计算重金属含量数据在各阈值条件下的指示半变异函数；

基于所述指示半变异函数建立先验条件累积分布函数；

将所述目标区域划分为同一分辨率的格网，定义一条经过所有格网的随机路径，并在第一个位置格网处从条件累积分布函数中随机抽取一个值作为模拟值；

将所述模拟值用于下一位置格网的先验条件累积分布函数，从下一位置格网的先验条件累积分布函数中随机抽取一个值作为模拟值，重复执行该步骤直至所有格网模拟完毕。

7.如权利要求5所述的土壤重金属风险预测方法，其特征在于，所述计算所述目标区域的生态风险指数包括以下步骤：

计算所述目标区域土壤中单种重金属的潜在生态风险指数；

计算所述目标区域土壤中多种重金属的综合潜在生态风险指数。

8.如权利要求7所述的土壤重金属风险预测方法，其特征在于，所述单种金属的潜在生态风险指数计算公式为

为单项重金属风险系数，C_i为表层土壤重金属浓度实测值，为参比值。

9.如权利要求8所述的土壤重金属风险预测方法，其特征在于，所述多种重金属的综合潜在生态风险指数计算公式为

其中，为第i种重金属潜在生态危害系数，为第i种重金属的毒性响应系数；RI为综合潜在生态危害指数。