本发明属于金属污染源定位技术领域，涉及一种确定城市表层土重金属污染源位置的方法。
背景技术：
2005年4月至2013年12月，环保部与国土资源部首次开展了全国范围内的土壤调查。结果显示，全国土壤总点位超标率为16.1％，其中轻微污染占11.2％，轻度污染占2.3％，中度污染1.5％，重度污染1.1％。镉、汞、砷、铅4种无机污染物含量分布呈现从西北到东南、从东北到西南方向逐渐升高的态势。镉、汞、砷、铜、铅、铬、锌、镍8种无机污染物点位超标率分别为7.0％、1.6％、2.7％、2.1％、1.5％、1.1％、0.9％、4.8％。
城市土壤是城市生态环境的重要组成部分，对城市生态系统稳定性起着关键作用。城市土壤由于人类活动导致环境变化，使得自然生态系统下的成土过程向人工生态系统下成土过程转变，其中城市重金属污染是城市土壤研究的热点。城市土壤重金属分布在时间和空间上存在显著差异，在城市布局和产业结构上也有很大差异，导致重金属污染源位置难以确定。关于城市重金属污染源位置确定模型是国内外研究的热点问题，但至今没有得到特别好的结果。
重金属的传播特征可理解为它的空间分布和动态传播规律。元素分布特点研究可通过地理信息系统(gis)完成，但gis中显示的重金属浓度最高点未必是污染源，而重金属动态变化的难度在于影响因子的多样性、土壤的不均匀性及区域土壤不同等。
技术实现要素：
为实现上述目的，本发明提供一种确定城市表层土重金属污染源位置的方法，解决了现有技术中重金属污染源位置难以确定，地理信息系统显示不准确的问题。
本发明所采用的技术方案是，确定城市表层土重金属污染源位置的方法，具体按照以下步骤进行：
步骤一、模型假设：假设城区的重金属污染浓度是一种连续的稳定场；污染源以点/面污染源为主；污染源扩散方式为同心圆或扇形，且有一定的重金属扩散半径；污染物的浓度变化方式不受高度影响，并且以相同的梯度减小；
步骤二、模型建立：首先制作污染元素浓度分布图，其次定义相关参数，制作一维矢量方程，并判定污染源范围，最后对边界范围污染程度进行估计，从而确定城市表层土重金属污染源的位置。
进一步的，所述步骤二中污染元素浓度分布图的制作步骤是：根据元素分布特点，建立污染元素浓度分布图，污染元素浓度分布图可通过地理信息系统、matlab或surfer软件完成。
进一步的，所述步骤二定义相关参数的步骤是：污染元素浓度分布图制作完成后，直观的看出重金属的浓度在不同位置的分布，在重金属迁移过程中，li为重金属在表层土迁移的垂直最大距离；χ为重金属行走极限点，即某一元素达到土壤环境质量一级标准的背景值点，背景值点的浓度为ρχ；χ到最大浓度点的直线为la，la的长度为l；ω为扩散能量极大点，即扩散能力能够穿透10cm地面的点；ρ为重金属浓度，则任一点重金属浓度为ρi其中，i＝1，2，3，...，n，ρmax为最大浓度点的浓度。
进一步的，所述步骤二制作一维矢量方程的步骤是：选取背景值点坐标(x1,y1)，再选取周围浓度大于背景值的点(x2,y2)，(x3,y3)……(xn-1,yn-1)，(xn,yn)对其做一元回归方程y＝β0+β1x，其中，β0为节距，β1为斜率，一元回归方程为χ→ω的一维矢量方程。
进一步的，所述步骤二确定污染源范围的步骤是：在污染元素浓度分布图中取任意元素的局部等浓梯度变化图，找出该元素的背景值浓度ρχ，若污染源重金属最大浓度点落在一维矢量方程直线上，则该点就为所要找的污染源点；若污染源重金属最大浓度点不在一维矢量方程上，以最大浓度点为圆心，圆心点坐标(x0，y0)，向一维矢量方程引垂线，垂线为半径r作圆，圆与直线la相交的点称为边界点，垂线与一元回归方程的交点是扩散能量极大点ω，r为最可疑半径，最终确定污染源就在最可疑半径内。
进一步的，所述半径
进一步的，所述步骤二边界范围污染程度估计的步骤是：污染源重金属最大浓度点不在一维矢量方程上，确定污染源最小范围后，根据模型评估污染源的危害程度，污染物浓度扩散的分布是均等分布的，污染物的浓度和扩散梯度土壤剖面的理想面积成正比，一个微理想面积等于一个lidx，扩散梯度土壤剖面的理想面积由若干个微理想面积lidx，dx为任意相邻两点之间的间隔，组成重金属行走极限点χ形成的浓度梯度角为
当取x1时，l1＝x1tanθ，面积
当取x2时，l2＝x2tanθ，面积
……
当取xi时，li＝xitanθ，面积
扩散梯度土壤剖面总面积和浓度的关系：
由于ρmax为已知浓度，边界点浓度为
进一步的，所述边界点浓度和最大浓度的平均值
本发明的有益效果是，与现有技术相比，本发明基于几何和区域散点线性关系，建立的重金属动态传播模型，能够将重金属污染源区域缩小到更小的范围，并能通过随机取点来找污染源。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是污染源重金属元素浓度理想分布图；
图2是as污染元素浓度分布图；
图3是zn污染元素浓度分布图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
同一污染源，不同元素的迁移变化趋势是不同的。确定城市表层土重金属污染源位置的方法，具体按照以下步骤进行：
步骤一、模型假设：假设城区的重金属污染浓度是一种连续的稳定场；污染源以点/面污染源为主；污染源扩散方式为同心圆或扇形，且有一定的重金属扩散半径；污染物的浓度变化方式不受高度影响，并且以相同的梯度减小；
步骤二、模型建立：
s1、污染元素分布制作：根据元素分布特点，建立污染元素浓度分布图，污染元素浓度分布图可通过地理信息系统(gis)、matlab或surfer软件完成；
s2、定义相关参数：污染元素浓度分布图制作完成后，直观的看出重金属的浓度在不同位置的分布，由于污染源的不确定性，gis中显示的重金属最大浓度点未必是污染源，假设该地区的污染源在一个圆形区域内，在重金属迁移过程中，li为重金属在表层土迁移的垂直最大距离(规定表层土厚度是10～20cm的土层，本发明实施例选择10cm作为计算依据)；χ为重金属行走极限点，即某一元素达到土壤环境质量一级标准的背景值点，背景值点的浓度为ρχ，当重金属元素到达重金属行走极限点χ的时候，不在向外扩散；χ到最大浓度点的直线为la，la的长度为l；ω为扩散能量极大点，即扩散能力能够穿透10cm地面的点；ρ为重金属浓度，则任一点重金属浓度为ρi，其中i＝1，2，3，...，n，ρmax为最大浓度点的浓度；
s3、制作一维矢量方程：选取背景值点坐标(x1,y1)，再选取周围浓度大于背景值的点(x2,y2)，(x3,y3)……(xn-1,yn-1)，(xn,yn)对其做一元回归方程y＝β0+β1x，其中，β0为节距，β1为斜率，如图1所示；一元回归方程为χ→ω的一维矢量方程；
s4、确定污染源范围：在污染元素浓度分布图中取任意元素的局部等浓梯度变化图，找出该元素的背景值浓度ρχ，χ点向比它浓度更大的地方做χ→ω行走y＝β0+β1x，其中，x和y是有界域的值，若污染源重金属最大浓度点落在一维矢量方程直线上，则该点就为所要找的污染源点；若污染源重金属最大浓度点不在一维矢量方程上，如图1所示，以最大浓度点为圆心，圆心点坐标(x0，y0)，向一维矢量方程引垂线，垂线为半径r作圆，其中，
s5、边界范围污染程度估计：污染源重金属最大浓度点不在一维矢量方程上，确定污染源最小范围后，根据模型评估污染源的危害程度，污染物浓度扩散的分布是均等分布的，污染物的浓度和扩散梯度土壤剖面的理想面积成正比，一个微理想面积等于一个lidx，就是一个微积分，理想面积由若干个微理想面积lidx组成，dx为任意相邻两点之间的间隔，dx的大小可根据《关于调整城市规模划分标准的通知》(国发2014第51号文件)，通过人口对城市规模进行分类，超大城市：城区常住人口1000万以上，规定dx＝100m；特大城市：城区常住人口500至1000万，规定dx＝80m；大城市：城区常住人口100至500万，规定dx＝40m；中等城市：城区常住人口50至100万，规定dx＝20m；小城市：城区常住人口50万以下，规定dx＝1m；重金属行走极限点χ形成的浓度梯度角为
当取x1时，l1＝x1tanθ，面积
当取x2时，l2＝x2tanθ，面积
……
当取xi时，li＝xitanθ，面积
扩散梯度土壤剖面总面积和浓度的关系：
由于ρmax为已知浓度，边界点浓度为
表1土壤环境质量评价标准限值单位mg/kg
其中，a级标准为土壤环境质量目标值，代表了土壤未受污染的环境水平，符合a级标准的土壤可以适用于各土壤利用类型；b级标准为土壤修复行动值，当某场地土壤污染物检测值超过了b级标准限值，该场地必须实施土壤修复工程，使之符合a级标准。
最后，取边界点浓度和最大浓度的平均值
实施例1
s1、污染元素分布制作：
某县人口为16.5万，属于中、小城市，元素砷(as)分布图通过matlab进行绘制，如图2所示，其中三角形点表示浓度小于背景值的点，圆形点表示浓度大于背景值的点；
s2、制作一维矢量方程，确定污染源范围：
如图2所示，浓度对比图可知as的背景值浓度为5.4ug/g，选取背景值点χ，χ在图中为4，坐标(5640，7998)，再选取周围浓度大于背景值浓度的点1的坐标(4750，7400)、点2的坐标(4780，7631)和点3的坐标(4930，7250)，对其做出一元回归方程y＝0.622x+4442，相关系数r2＝0.636，则该方程为χ→ω的一维矢量方程，看出最大浓度点2并不在一维矢量方程上，因此以点2为圆形，向一维矢量方程引垂线，如图2所示，则r的距离为：
s3、边界范围污染程度估计：
污染源重金属最大浓度点不在一维矢量方程上，确定污染源最小范围后，根据模型假设评估污染源的危害程度，根据距离公式计算背景值点和最大浓度点之间的距离，背景值点坐标为(5640，7998)，最大浓度点坐标为(4780，7631)，最大浓度点和背景值点的距离
根据扩散梯度土壤剖面总面积和浓度公式：
背景点到边界点的距离为935-183＝752m，其中xmax、xmax+1、xi、xi+1均为l中的一个微理想长度，而xmax、xmax+1间隔为1m，xi、xi+1间隔为1m，xmax+1＝935，那么xmax＝934，xi+1＝752，那么xi＝751，则
由浓度对比图可知ρmax的浓度为21.87mg/kg，则边界点浓度为：
可通过《展览会用地土壤环境质量评价标准(暂行)》(hj350—20072007-08-01实施)中的(表1)进行评价，边界浓度和最高浓度平均值为(21.87mg/kg+17.64mg/kg)/2＝19.805mg/kg<20mg/kg，因此该区域as浓度符合a级。
搜查结果表明，重金属as主要分布于工业园区和道路交通较为发达的地区，最终发现该区块污染源为一排水管道，该管道是附近一家化肥厂所有。
可看出，污染源的范围从半径为935m的圆的范围缩小到了半径为183m的圆，在污染源的搜索上大大节省了人力物力。
实施例2
s1、污染元素分布制作：
某市人口为700万，属于特大城市，元素锌(zn)分布图通过matlab进行绘制，如图3所示：
s2、制作一维矢量方程，确定污染源范围：
如图3所示，浓度对比图可知zn的背景值为32.86mg/kg，选取三个重金属独立分布的区域，选取图3中1区、2区、3区的背景值点，选取周围浓度大于背景值的点分别做出三个区域的一元回归方程：
第1区块：y1＝1.24x1+1800；
第2区块：y2＝3.41x2-20460；
第3区块：y3＝-0.75x3+16875
则以上方程为χ→ω的一维矢量方程，由于最大浓度点并不在一维矢量方程上，分别以3个区块的最大浓度点为圆心，向各个区域的一维矢量方程引垂线，根据点到直线的距离公式，
s3、边界范围污染程度估计：
污染源重金属最大浓度点不在一维矢量方程上，确定污染源最小范围后，根据模型假设评估污染源的危害程度，根据距离公式计算背景值点和最大浓度点之间的距离，计算最大浓度点和背景值点的距离l1＝7660m，l2＝4320m，l3＝5000m。第1区块的背景值点到边界点的距离为7600-754＝6846m；第2区块的背景值点到边界点的距离为4320-258＝4062m；第3区块的背景值点到边界点的距离为5000-433＝4567m。
根据公式：
其中xmax、xmax+1、xi、xi+1均为l中的一个微理想长度，xmax、xmax+1间隔为80m，xi、xi+1间隔为80m，则
第一区块：
其中，xmax+1＝7600，那么xmax＝7520，xi+1＝6846，那么xi＝6766；
第二区块：
其中，xmax+1＝4320，那么xmax＝4240，xi+1＝4063，那么xi＝3983；
第三区块：
其中，xmax+1＝5000，那么xmax＝4920，xi+1＝4567，那么xi＝4487；
由浓度对比图可知：
第1区块：ρmax＝726.56mg/kg；
第2区块：ρmax＝2893.78mg/kg；
第3区块：ρmax＝3760.15mg/kg；
则边界点浓度为：
第1区块：
第2区块：
第3区块：
可通过《展览会用地土壤环境质量评价标准(暂行)》(hj350—20072007-08-01实施)中的表1进行评价：
第1区块：边界浓度和最高浓度平均值为(621.4mg/kg+726.56mg/kg)/2＝673.98>200mg/kg，因此该区域zn浓度不符合a级。
第2区块：边界浓度和最高浓度平均值为(2588.1mg/kg+2893.78mg/kg)/2＝2740.94>200mg/kg，因此该区域zn浓度不符合a级。
第3区块：边界浓度和最高浓度平均值为(3161.3mg/kg+3760.15mg/kg)/2＝3460.73>200mg/kg，因此该区域zn浓度不符合a级。
搜索结果表明，重金属zn主要分布于工业园区和道路交通较为发达的地区，最终发现第1区块重金属污染来自于道路交通网络地区，第2区块重金属污染来自于一家仪器仪表场所排放的污水，第3区块重金属污染来自于工业园区内的一家农药厂。
以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均包含在本发明的保护范围内。