城市表层土壤重金属污染分析
摘要
随着城市经济的快速发展和人口的不断增加，人类活动对城市环境质量的影响增大。日常生活废弃物的排放和工厂的制作加工活动向土壤排放了大量的重金属。当重金属的排放量超出正常范围，环境质量也随之恶化。因此，开展对土壤重金属污染的评价，分析重金属污染物的传播特征，确定污染源的位置等对于我们有效地控制土壤重金属含量，改善环境来说是十分重要的。
首先，利用surfer软件对附件中8种重金属元素的数据进行加工处理，分别画出它们在该城区的空间分布。发现污染浓度高的地区大多集中在工业区和主干道路区，而山区的重金属浓度大多小于其背景值。本文通过对单因子指数法和内梅罗综合污染指数法进行改进，建立了新的单因子污染指数及综合污染指数公式。得出结论：8种重金属中Hg和Cu的综合污染指数最高；工业区（污染指数为7.052）和主干道路区（4.566）为重度污染区，生活区（1.670）和公园绿地区（1.246）为中度污染区，山区为轻度污染区（0.325）。
根据数据的分析和第一问的结论（8种重金属元素中以Gu和Hg的污染指数为最高，分别为3.193和7.673，均属于重度污染物），可以得到城区重金属的污染主要是由于工业生产排放的废弃物，农药和肥料的施用，汽油的燃烧，汽车尾气的排放等。
然后，根据不同区域以及海拔高度在该城区的分布图，分析该城区的地势分布特点以及区域分布特点。在地表的重金属一方面随着地表水流动向四周传播，由地势高的地区传播到地势低的地区；另一方面，由于分子运动，重金属本身也会向四周扩散。以该城区污染最严重的重金属Hg和Cu为分析对象，根据重金属污染物的传播特性，建立一个二维对流-扩散方程。对该偏微分方程采用中心差分格式离散，通过matlab编程计算得到污染源的一个数值模拟结果。由此得出Cu的污染源在（2578.86±286.54，3320.82±184.49）以及（3151.94±286.54，5903.68±184.49），Hg的污染源在（16046.24±286.54，5±184.49）以及（2292.32±286.54，3874.29±184.49）。
最后，为了更好地研究城市地质环境的演变模型，还应收集不同金属浓度随时间的变化值，各点地表水的流速，重金属的降解系数以及扩散系数。对二维对流-扩散模型的缺点进行模型改进，得到一个更有说服力的污染源确定模型。
关键词：污染程度；单因子指数法；内梅罗综合污染指数法；二维对流-扩散方程；最临近插值法
一、问题重述
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日益突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。
在城市众多污染物中，以重金属污染最为严重。重金属污染指由重金属或其化合物造成的环境污染。主要由采矿、废气排放、污水灌溉和使用重金属制品等人为因素所致。因此人类活动导致环境中的重金属含量增加并超出正常范围，必定导致环境质量恶化。
通过对城市的分布范围以间距为1公里左右的网格进行划分，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并用GPS记录采样点的位置，得到8种主要重金属的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。
城区按照功能划分，一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。
根据附件1采样点的位置、海拔高度及其所属功能区等信息，附件2中8
种主要重金属元素在采样点处的浓度，以及附件3中8种主要重金属元素的背景值。通过建立数学模型来完成以下任务：
(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
(2) 通过数据分析，说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。
(4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？
二、问题分析
问题一，已知8种主要重金属在该城区的采样点的位置坐标和对应的浓度，可以用surfer软件画出主要重金属在该城区的空间分布。评价重金属的污染程度与不同重金属的浓度之间的关系，目前学术界主要使用的方法大致包括生态危害指数法、地积累指数法、单因子指数法和内梅罗综合污染指数法[1]、模糊贴近度方法、潜在生态危害指数法等等。这些评价方法在评价土壤重金属污染程度上都曾被运用过。其中潜在生态危害指数的计算方法，需要计算相关的含量条件，数量条件，毒性条件和敏感性条件。在金属毒性系数的确定中，应对“丰富原则”和“释放效应”加以讨论。然而附件所给条件不甚具体，如土地的物质构成对“丰富原则”有很大影响，所以限制性很大，无法采用这种方法。而对于地积累指数
法和生态危害指数法，其的确能够测量出各个重金属元素的毒性含量和对土壤的污染程度，然而无法将多个重金属元素对土壤的综合影响表示出来，也便无法阐述对土壤重金属污染的整体评价。因此，本文对单因子指数法和内梅罗综合污染指数法进行改进，得到新的单因子污染指数以及综合污染指数公式，得出不同区域中各种重金属的污染指数，以此数据得出该城区内不同区域重金属的污染程度。
问题二，根据问题一所得出的单个重金属对不同区域的污染指数大小，分析该区域污染指数相对较大的重金属元素污染的主要原因。计算各种金属对整个区域的污染指数大小，同时结合该城区内不同区域重金属的污染程度，得出在整个城区中哪几个元素的污染指数是最高的，分析其原因。
问题三，标出不同区域以及海拔高度在该城区的分布图，分析该城区的地势分布特点以及区域分布特点。一方面表层的重金属主要会随着地表水的流动向四周传播，由海拔高的地区向海拔低的地区传播；另一方面，由于分子运动，重金属也会向四周扩散。考虑上述几个问题所得出的对城区污染最严重的几个重金属为研究对象，根据以上关于重金属污染物的传播特性分析，可建立二维对流-扩散方程来描述这一过程。得到污染源的数值模拟结果。
问题四，由于附件中的数据没有随时间的一个变化值，因此二维对流-扩散方程模型是在不考虑浓度随之间变化这个变量值。建立重金属的降解系数以及扩散系数模型，根据上述模型的缺点进行模型改进，得到一个更有说服力的污染源确定模型。
三、模型假设
1、假设中国绿色食品发展中心《绿色食品产地环境质量现状评价纲要(试行)》中的污染程度分级标准是有效的；
2、假设附件中的数据是在重金属污染达到一个稳定状态时所采集的；
3、假设通过资料所查得的地表水的流速范围是有效的；
4、假设重金属降解量很少，可以忽略其降解的影响。
四、符号说明以及定义
五、模型建立与求解
5.1重金属的空间分布以及污染程度
5.1.1重金属的空间分布
通过附件中采样点的坐标以及对应浓度，用surfer 软件画出8种主要重金属元素在该城区的空间分布。同时，也将每一种重金属的浓度分成不同范围，得出不同范围内数据的个数以及位置，通过matlab 画出某一种重金属在不同浓度范围内的散点图，图见附录一。以下是8种主要重金属元素的空间分布图，左侧为三维立体图，右侧对应等值线图。
图1 As 空间分布图
由图可知浓度较低的As 分布在中部和东北部（上北下南，左西右东，下同）， 浓度<4（mg/kg ）,而浓度很高的大多分布在西部，北部以及南部个别地区如工业
N
区和主干交通区。
图
2 Cd 空间分布图
可以看到，高浓度的
Cd （浓度大于0.6（
mg/kg ）
）主要分布在西部主干交通区和生活区。相对地，中东部山区的Cd 浓度较低，范围大致在0（mg/kg ）~0.1（mg/kg ）间。东北部偶有高浓度Cd 地点出现。
图3 Cr 空间分布图
从图中可以看出Cr 的浓度分布较均匀，浓度在20（mg/kg ）至40（mg/kg ）的数量较多。浓度大于100（mg/kg ）的Cr 集中在西南部公园绿地区，极少在中部。
图4 Cu 空间分布图
图中Cu 的浓度分布很明显—高浓度Cu 集中分布在西南部公园绿地区，其浓度值大于70（mg/kg ）。而其他地区的Cu 浓度分布均匀,致小于20（mg/kg ）。
图5 Hg 空间分布图
Hg 的浓度在图中出现三个高峰，一个位于西南部，另外两个位于中东部，其他地区的浓度分布平均，大多处于0.07（mg/kg ）~1（mg/kg ）和0.02（mg/kg ）~0.04（mg/kg ）段，且大都分布在主干道路区，生活区和公园绿地区。
图6 Ni 空间分布图
结合原始数据，图中大部分Ni 浓度都符合标准，而且地理位置的分布较均匀，个别浓度高峰出现在西南部公园绿地区，浓度范围大致为16（mg/kg ）~25（mg/kg ）。
图7 Pb 空间分布图
由图可以明显地看出，高浓度的Pb 集中分布在西部，中部次之，东部浓度值最低，总体呈由西到东的阶梯状分布。东北部山区的浓度值很低。
图8 Zn 空间分布图
由上图知，浓度符合标准的Zn 点分布仍然集中分布在中部和西南部主干道路区和公园绿地区，最高值出现在中部地带，而大部分浓度很高的点也仍旧分布在西南部。
5.1.2重金属污染程度评价模型
评价重金属的污染程度与不同重金属的浓度之间的关系，目前学术界主要使用的方法大致包括生态危害指数法、地积累指数法、单因子指数法和内梅罗综合污染指数法、模糊贴近度方法、潜在生态危害指数法等等。这些评价方法在评价土壤重金属污染程度上都曾被运用过。
其中潜在生态危害指数的计算方法，需要计算相关的含量条件，数量条件，毒性条件和敏感性条件。在金属毒性系数的确定中，应对“丰富原则”和“释放效应”加以讨论。然而附件所给条件不甚具体，如土地的物质构成对“丰富原则”有很大影响，所以限制性很大，无法采用这种方法。
而对于地积累指数法和生态危害指数法，其的确能够测量出各个重金属元素的毒性含量和对土壤的污染程度，然而无法将多个重金属元素对土壤的综合影响表示出来，也便无法阐述对土壤重金属污染的整体评价。
因此，本文采用单因子指数法和内梅罗综合污染指数法对城市重金属污染程度进行评价，通过其计算式可以评价不同土地类型的土壤受到重金属污染的强度。
（一）单因子指数法
以第i 种重金属的实际测量的浓度与其背景值的比值来表示该重金属的污染指数，公式为
i
i i
C P S
。 式中，i P 为土壤中第i 种重金属元素的污染指数，i C 为土壤中第i 种重金属元素的实测浓度，i S 为区域中第i 种重金属的环境背景值。 （二）内梅罗综合污染指数法
P =
综 式中，n 为参与评价的重金属种类总数，P 综为土壤中重金属元素的综合污染指数。
通过对上述两种方法的分析，可以发现如果用实际测量值与背景值的差值占背景值的大小作为重金属的污染指数，比单单以实际测量值占背景值的大小作为污染指数更有说服力。因此，本文对原有的单因子指数法进行改进，得到如下计算公式：
1,0,ij j ij
ij j j j ij ij j C S C C S S S P C S -?=->?=??≤?
， 重金属的背景值是指重金属的正常含量。当重金属的实测值小于其背景值时，其污染指数为0。
在评价重金属的综合污染指数时，实际数据中可能会有异常值存在，若以
(
)i
i
C Max S 来计算，会发现数据中存在特别大的值，这样会使求得的综合评价指数偏大，误差较大。因此，本文对原有的内梅罗综合污染指数法进行改进，得到如下计算公式：
P =
i 综，
式中ij P 为第i 区域第j 种重金属的污染指数，ij C 为第i 区域第j 种重金属的实测浓度，P i 综
为第i 区域的综合评价污染指数。 等级划分 1
2 3 4
单因子污染指数 ij P =0
0
污染程度
无污染
轻度污染
中度污染
重度污染
表3 土壤综合污染程度分级标准
等级划分 1 2 3 4
综合污染指数P i综=0 0
2 污染程度无污染轻度污染中度污染重度污染
依据：中国绿色食品发展中心《绿色食品产地环境质量现状评价纲要》(1994年)
5.1.3模型处理
由于附件里同区域中的各重金属数据量比较大，因此，在对不同重金属污染指数的处理过程中，我们采用第i种重金属污染指数的平均值作为其污染指数。下表中的数据都是以实际值大于背景值进行处理的。
表4 生活区土壤中各重金属的污染指数（每种样本个数44个）
As Cd Cr Cu Hg Ni Pb Zn 背景值 3.6 0.13 31 13.2 0.035 12.3 31 69 和120.68 7.238 1693.6 1597.8 2.728 281.69 1695.2 7474 平均值 2.743 0.165 38.490 36.313 0.062 6.402 38.526 169.9 污染指数0.726 1.265 1.242 2.751 1.771 0.520 1.243 2.462 污染程度轻度中度中度重度中度轻度中度重度表5 工业区土壤中各重金属的污染指数（每种样本个数36个）
As Cd Cr Cu Hg Ni Pb Zn 背景值 3.6 0.13 31 13.2 0.035 12.3 31 69 和135.59 9.500 845.88 4116.6 21.945 287.05 2233.5 7534 平均值 3.766 0.264 23.497 114.35 0.610 7.974 62.041 209.3 污染指数 1.046 2.030 0.758 8.663 17.417 0.648 2.001 3.033 污染程度中度重度轻度重度重度重度重度重度表6 山区土壤中各重金属的污染指数（每种样本个数66个）
As Cd Cr Cu Hg Ni Pb Zn 背景值 3.6 0.13 31 13.2 0.035 12.3 31 69 和56.200 2.589 765.63 363.71 0.747 306.52 515.94 796.19 平均值0.852 0.039 11.600 5.511 0.011 4.644 7.817 12.063 污染指数0.237 0.302 0.374 0.417 0.323 0.378 0.252 0.175
续表6 山区土壤中各重金属的污染指数（每种样本个数66个）
As Cd Cr Cu Hg Ni Pb Zn
污染程度轻度轻度轻度轻度轻度轻度轻度轻度
表7 主干道路区土壤中各重金属的污染指数（每种样本个数138个）
As Cd Cr Cu Hg Ni Pb Zn 背景值306.83 32.239 3809.3 6765.3 57.308 788.0 4558.6 24208 和135.59 9.500 845.88 4116.6 21.945 287.1 2233.5 7534.
平均值 2.223 0.234 27.604 49.024 0.415 5.710 33.034 175.4 污染指数0.618 1.797 0.890 3.714 11.865 0.464 1.066 2.542 污染程度轻度中度轻度重度重度轻度中度重度
表8 公园绿地区土壤中各重金属的污染指数（每种样本个数35个）
As Cd Cr Cu Hg Ni Pb Zn 背景值 3.6 0.13 31 13.2 0.035 12.3 31 69 和95.020 5.402 469.84 601.39 2.941 124.63 1046.1 3111 平均值 2.715 0.154 13.424 17.183 0.084 3.561 29.890 88.91 污染指数0.754 1.187 0.433 1.302 2.401 0.290 0.964 1.288 污染程度轻度中度轻度中度重度轻度轻度中度
根据不同区域的不同重金属污染指数以及综合评价污染指数的公式，可以得出五个区域的综合评价污染指数，值如下表所示。
表9 五个区域的综合评价污染指数
区域山区公园绿地区生活区主干道路区工业区
综合评价污
0.325 1.240 1.670 4.566 7.052
染指数
污染程度轻度污染中度污染中度污染重度污染重度污染
从表中可以看出不同地区的污染程度从大到小依次为：工业区（7.052）、主干道路区（4.566）、生活区（1.670）、公园绿地区（1.240）、山区（0.325）。其中，工业区和主干道路区为重度污染区域，这与工业区和主干道路区的环境特点相符：在工业区集中了大部分该城区的工厂，必然会有大量的重金属排放。在主干道路区，来往的车辆所排放的有毒物质及汽车燃料的燃烧等也产生很多重金属。生活区和公园绿地区为中度污染区域，其与公园绿地区都是人类活动密集的地方，必然会产生一定量的重金属；山区为轻度污染区域：山区的人类活动相较较
少，重金属的浓度自然也就相对其他地区少。通过上述分析，本文所得到的不同区域的重金属污染程度与实际较符合。
5.2重金属污染的主要原因
土壤中重金属的来源是多途径的，除成土嚷本身含有重金属外，人类工农业生产活动也造成重金属对大气、水体和土壤的污染。
表10 城区各重金属的污染指数
As Cd Cr Cu Hg Ni Pb Zn
背景值 3.6 0.13 31 13.2 0.035 12.3 31 69 和714.3256.977584.2513444.785.671787.9210049.443124平均值 2.2390.17923.77542.1470.269 5.60531.503135.18污染指数0.6220 1.3737 0.7669 3.1929 7.6730 0.4557 1.0162 1.9592 污染程度轻度中度轻度中度重度轻度轻度中度
由表10所示，铜，汞的综合污染指数分别为3.1929和7.6730，大于2，属于重度污染。锌的污染指数为1.9592，接近于2。结合城区各区域分布图
图9 城区各区域分布图
(1)该城市的主干道路分布范围广而均匀，交通污染带来的影响大；工业区成条状带比较集中，污染量大；地势平坦，污染物容易扩散。
(2)锌污染：主要由于工业区的锌矿开采，冶炼加工，机械制造，仪器仪表；交通区的汽车轮胎磨损；生活区中煤燃烧产生的粉尘等。
(3)铜污染：主要由于工业区冶炼排放的烟尘，矿区的开采；交通区的废弃排放，含铅汽油的燃烧等。
(4)汞污染：主要来源于工业区的集中狭长排放了大量汞污染，如仪表厂，食盐电解，贵金属冶炼；生活区的化妆品、照明用灯、齿科材料、燃煤、水生生物；交通区的尾气排放等以及公园绿地区的污水灌溉等。
下面分区介绍主要污染元素和原因：
生活区的主要污染元素为汞，铜，锌。在农作物种植过程中，施用含有铅、汞、镉、砷等的农药和不合理地施用化肥，都可以导致土壤中重金属的污染。近
年来污水灌溉已成为农业灌溉用水的重要组成部分，污水灌溉导致土壤重金属Hg、Cd、Cr、As、Cu、Zn、Pb等含量的增加；生活重的污中的污泥施肥含有大量的有机质和氮、磷、钾等营养元素，但同时污泥中也含有大量的重金属,随着大量的市政污泥进入农田导致土壤中Cd、Hg、Cr、Cu、Zn、Ni、Pb含量的增加。生活区中也有含重金属废弃物如化妆品、电池、燃煤等的堆积也使汞，铜，锌的含量超过正常值。
工业区的主要污染元素是铜，锌，汞，铅等。一方面，工业排放了大量的含有铜、锌、汞的污水，并且重金属废弃物种类繁多，主要包含了弃旧的铜，锌，汞，铅等元素。另一方面，工业区金属矿山的开采、冶炼、重金属尾矿、冶炼废渣和矿渣堆放等，被酸溶出含重金属离子的矿山酸性废水，随着矿山排水和降雨使之直接进入土壤，致使土壤中铜、锌、汞、铅等含量增高。
山区属于轻度污染。因为山区的地势相对较高，离工业区较远，我国季风均为偏北风或偏南风，很少受到工业区影响。然而主干道路已经深入到山区，主干道路污染的扩散使得山区呈轻度污染。
交通区属于重度污染，主要污染原因是铜、锌、汞。含铅汽油的燃烧，汽车轮胎磨损产生的含锌粉尘含有铜，锌，汞，铅等物质。
公园绿地区重金属污染的主要因素是汞、锌、镉。图中所示的公园绿地区周围有交通枢纽也有生活区，镉和汞的产生主要来自于交通运输的过程。有机肥，化肥和农药的大量使用是土壤中锌污染的主要途径。
5.3确定污染源位置
5.3.1重金属污染物的传播特征
从5.2中已经知道不同的重金属中，尤其以Cu和Hg的污染指数最高，分别为3.19和7.67，远大于其它6种重金属元素，这也意味着这两种重金属元素对该城区的污染程度的大小起着确定性的作用。Cu和Hg的浓度不超过其背景值的大部分集中在山区，只有极小部分是分布在其它 4 个区域。浓度高的样本点大多集中在工业区和交通区，其中以工业区最为密集。
工业区的范围大致是东北至西南的部分地区，山区集中在中东部地区，而主干道路区被山区和工业区所包围，分布在该城市的各个区域。污染物的浓度大致以工业区和主干道路为中心向四周递减。
大气中的大多数重金属经过自然沉降和雨淋沉降进入土壤。同时，由于城市工业化的迅速发展，大量的工业废水涌入河道，使城市污水中含有的大量重金属离子随着污水灌溉而进入土壤。在分布上，靠近污染源头和城市工业区的土壤污染往往严重，远离污染源头和城市工业区的土壤污染程度较低，矿山排水和降雨将重金属带入水环境（如河流等）或直接进入土壤，这些都可以间接或直接地造成土壤重金属污染，又由于金属元素迁移能力减弱和水体自净化能力的适度恢复，金属化学污染强度逐渐降低。
在土壤中的各类重金属主要随着地下水和地表水向特定方向传播，同时，由于分子运动，重金属本身也会向四周扩散。重金属也会随着风向和风速的不同向不同的方位进行传播。但是，风向和风速时时在变，不容易进行定量和定性的分析，因此，忽略风向和风速的影响。
5.3.2二维对流-扩散方程
根据5.3.1对重金属的传播特征的分析，可以发现重金属在土壤中的传播与
地下水和地表水的水流速度以及其本身的扩散系数有关。由此可以得出，重金属污染物在土壤中的传播，其实相当于其在水中的传播。根据污染物质在河流中的一维对流-扩散方程[2][3]：
2
2(,)C C C v E f x t k C t x x
???+=+-???， (1) 式中，C 为河流中污染物质的浓度，v 为该段河流的水流速度，f 为外界输入的
污染物浓度，E 为污染物物质的扩散系数，k 为污染物质的降解系数。
根据这个方程，我们可以发现若以河流的长度为x 轴，河流在某一点0x 的流速与该点河宽的关系不大。同时，污染物质自身的扩散系数会对河流的周围造成影响，但是，河流的流速以及其走势大致决定了污染的扩散方向，污染物浓度会随着横向的发生变化，一般情况下是不会对纵向造成影响，因此，可以忽略纵向因素，只建立关于x 轴的一维对流扩散方程。
本文是对城市表层土壤重金属污染进行分析，因此只要研究地表水对传播的影响。由于地表水分布在该城区的各个角落，其流速在没有特殊情况下，一般会比河流的流速要小得多，而且其流向与地势有关。
图10 地势空间分布图
从图10中可以看出，山区的地势较高，其它区域的地势的高低情况也都不相同，区域中地势的走向一直在变化。地表水的流向是从高地势流向低地势，这也就意味着流速会随着地势的变化而发生改变。
通过上述的分析，本文对污染物质在河流中的一维对流扩散方程进行改进，考虑重金属浓度随着纵向和横向变化的影响，建立重金属在土壤中的二维对流-扩散方程。
2222
(,,)k k k k k k
k k k C h C h C C C v v E E f x y t C t x x y y x y λ???????--=++-???????， (2) 式中，k C 为第k 种重金属的浓度，v 为地表水的流速，k E 为第k 种重金属的扩散系数，(,,)f x y t 为外界输入的第k 种重金属污染物质在某一区域的浓度，
λ为降解系数，h 为地势的高度。地表水流速随着地势增高而减小。由于实际数据中重金属的浓度是在某一时刻所采取得，并没有随时间变化的浓度值，因此，
本文假设此次所采取的数据是在重金属的污染浓度达到一个稳定状态时的一组
数据。这也就意味着0k
C t
?=?。而重金属物的一个很重要的特性是其难降解，根据这一特性，本文假设λ=0。方程（2）简化为
2222
(,)k k k k
k k k h C h C C C v v E E f x y x x y y x y ??????--=++??????。 (3)
5.3.3差分方程反演模型的建立与求解
对对流-扩散方程（3）采用欧拉差分格式离散。进行处理，两边积分，积分区间为（1j x -，j x ），其中2,3,,j n = 。含积分项利用梯形公式，含偏导数项利用一阶向前差商，j=n 时对距离偏导的差分用向后差商，得到如下公式:
,1,
,,1,,2,3,k k k i j i j k k
k
i j i j C C C x
x i j n C C C y
y --?-?=?
???=?-??=???? ， (4) ,1,,,1
,,2,3,i j i j i j i j h h h x
x i j n h h h y
y ---??=
????=?-??=???? (5) 21,,1,222,1
,,1222,,2,3,,2k k k k i j i j i j k k k k
i j i j i j C C C C x
x i j n C C C C y
y +-+-?-+?=
????=?-+??=???? (6) 将方程（4），(5)，（6）代入（3）中，得到1j f -的反演式，
,1,,1,,,1,,1
(,)(
)k k k k i j i j
i j i j
i j i j i j i j k C C h h C C h h f x y v x
x
y y --------=-?
+?????
1,,1,
,
1,,
12
2
22(
)i j i j i j i
j i j i
j k C C C C C C E x y +-
+-
-+-+-+
??。
(7) 式中，,2,3,,i j n = ，将城区分布图分成网格状其中i 代表横向节点，j 代表纵向节点。
5.3.4模型求解
从5.1的图形和采样点的数据分析中，我们可以发现，该城区的区域范围是不规则图形，而二维对流-扩散模型是在规则模型的基础上求解的。因此，我们
根据最临近插值法，可以快速的根据原有的点的数据补充完整该不规则图形。
从5.2中已经知道在不同的重金属中，尤其以Cu 和Hg 的污染指数最高分别为3.19和7.67，远大于其它6种的重金属元素，这也意味着这两种重金属元素对该城区的污染程度的大小起着确定性的作用，因此，本文通过求解Cu 和Hg 的外界输入的污染物浓度，得出重金属的大致污染源的位置。大致求解过程如下。
（1）根据已有Cu 数据运用最临近插值法，将城区Cu 的(x ，y ，C 1)以及(x ，y ，h ）从原先的散乱分布补充成一个长方形的网格图，得出各个小网格节点对应的(x ，y ，C 1)和(x ，y ，h)值。
（2）将（1）所得的网格各节点数据，计算出任意两点之间的f 值，通过MATLAB 进行计算得出重金属Cu 污染源的位置图，得出其大概的一个位置。
（3）同理可以得出重金属Hg 污染源的位置图，得出其大概的一个位置。 通过资料的查找，由[4]得到地表水的流速大致为0.05m/s~0.2m/s ，我们发现现今的重金属的扩散系数并没有一个大概的值，而其值的大小又与浓度随时间的变化有关，但是没有实际的测量数据，无法通过模型得出。本文通过从0~1中取不同大小的重金属的传播系数，画出Cu 和Hg 的污染源的平面图(见附录三)，发现其区域几乎不随扩散系数的改变而发生改变，因此，本文取Cu 和Hg 的扩散系数为0.52/m s 。
第一种情况:当v=0.05m/s ，Cu 的污染源位置确定。
,1,,1,,,1,,1(,)0.05(
)k k k k i j i j
i j i j
i j i j i j i j k C C h h C C h h f x y x
x y y
--------=-?
+?????
1,,1,,1,,1
2
2
220.5(
)i j i j i j
i j i j i j C C C C C C x
y
+-+--+-+-+
??
图11 Cu 污染源的立体图和平面图
第二种情况：当v=0.2m/s ，Cu 的污染源位置确定
,1,,1,,,1,,1
(,)0.2(
)k k k k i j i j
i j i j
i j i j i j i j k C C h h C C h h f x y x
x y y
--------=-?
+?????
1,,1,,1,,1
2
2
220.5(
)i j i j i j
i j i j i j C C C C C C x
y
+-+--+-+-+
??
图12 Cu 污染源的立体图和平面图
从图11和图12重金属Cu 污染源的平面图中，可以发现虽然地表水的流速有发生变化，但是对于污染源的位置确定的影响几乎看不出。立体图中可以看出Cu 的外界输入的污染物的浓度有两个值是特别大的，其它的值都小于它，因此，可以去认为Cu 的污染源有两个。平面图中有两白色区域为其污染源。根据每个网格数所代表的长度，得出其污染源的大致位置为（2578.86±286.54，3320.82±184.49）以及（3151.94±286.54，5903.68±184.49）。
第三种情况：当v=0.05m/s ，Hg 的污染源位置确定
,1,,1,,,1,,1(,)0.05(
)k k k k i j i j
i j i j
i j i j i j i j k C C h h C C h h f x y x
x y y
--------=-?
+?????
1,,1,,1,,1
2
2
220.5(
)
i j i j i j
i j i j i j C C C C C C x y +-+--+-+-+
??。
图13 Hg 污染源的立体图和平面图
第四种情况:当v=0.2m/s 时，Hg 的污染源位置确定
,1,,1,,,1,,1(,)0.2(
)k k k k i j i j
i j i j
i j i j i j i j k C C h h C C h h f x y x
x y y
--------=-?
+?????
1,,1,,1,,1
2
2
220.5(
)i j i j i j
i j i j i j C C C C C C x
y
+-+--+-+-+
??。
图14 Hg污染源的立体图和平面图
从图13和图14重金属Hg污染源的平面图中，可以发现虽然地表水的流速有发生变化，但是对于污染源的位置确定的影响几乎看不出。立体图中可以看出Hg的外界输入的污染物的浓度有两个值是特别大的，其它的值都小于它，因此，可以去认为Hg的污染源有两个。平面图中有两白色区域为其污染源。根据每个网格数所代表的长度，得出其污染源的大致位置为（16046.24±286.54，5±184.49）以及（2292.32±286.54，3874.29±184.49）
5.4二维对流-扩散模型评价与优化
5.4.1二维对流-扩散模型评价
本文在建立模型确定污染源的位置时，通过分析重金属的传播途径，考虑到由于分子运动，重金属本身也会向四周扩散，在地表的重金属主要随着地表水流动由高低是向低地势传播。以该城区污染最严重的重金属Hg和Cu为分析对象，根据重金属污染物的传播特性，建立一个二维对流-扩散方程。这样就避免了若单单考虑横向或纵向时，出现片面性，不能全方位的综合考虑重金属在整个城区中的一个随位置不同所产生的浓度变化。而二维对流-扩散方程就避免这种情况的发生。但是，由于附件中所提供的数据是某一时刻的数据。因此该模型是在假设该时刻重金属传播已达到稳定状态条件下建立的，这样模型得出的最后结果存在一定偏差。
5.4.2二维对流-扩散模型优化
从模型的缺点出发考虑，为了更好地研究城市地质坏境的演变模式，还应收集：（1）在相同的采样点位置处采集不同时间段的各重金属的浓度；
（2）测量该城区采样点的位置处的地表水的流速；
（3）收集重金属元素Gu和Hg的扩散系数以及降解系数。
根据所收集的数据将二维对流-扩散模型进行优化。
22
22
(,,)
k k k k k
k
k k k
C h C h C C C
v v E E f x y t C
t y x y y x y
λ
???????
--=++-
???????
依旧采用中心差分法求解偏微分,得出如下的公式。
22
22
(,,)
k k k k k
k
k k k
C h C h C C C
v v E E f x y t C
t y x y y x y
λ
???????
--=++-
???????
（8）
依旧采用中心差分法求解偏微分,得出如下的公式。
,,1,,1
1,,1,,,1,,1,2,2,3,,1,2,22k k k i j m i j m k k
k i j m i j m k k k i j m i j m
C C C t t C C i j n
C m s x
x C C C y
y +-+-+-?-?=????
?-=??=?
=????-?=?
???? ，（9） 1,,1,,,1,,1,,2,3,2,1,2,2i j m i j m
i j m i j m
h h h i j n x
x h h m s h y
y +-+--??=?=????-=??=???? ，（10） 21,,,,1,,222,1,,,,1,222,2,3,,1,2,2k k k k i j m i j m i j m
k k k
k i j m i j m i j m
C C C C i j n x
x m s C C C C y
y +-+-?-+?=
?=????=-+??=???? 。（11） 方程（9），(10)，（11）代入（8）中，得到1j f -的反演式，
,,1,,1
1,,1,,1,,1,,,,,,(,,)(
222k k k k i j m i j m i j m i j m
i j m i j m
k
k i j m
i j m
C C h h C C f
x y t v
t
x
x
+-+-+----=
-?
???
,1,,1,,1,,1,1,,,,1,,2
2)(22k k k k k i j m i j m i j m i j m i j m i j m i j m
k h h C C C C C E y y x +-+-+----++?-??? ,1,,,,1,,,2
2)k k k i j m i j m i j m
k i j m C C C C y
λ+--++
+?
式中，m 为时间的变化值，其余变量的含义不变。这个模型考虑到了重金属
浓度随时间的变化而变化以及地表水随时间和地点的变化而变化的情况。与原模型的求解过程相似，反推出污染源的位置。该模型考虑的更全面，得出的结果会更符合实际。
六．模型评价与改进
6.1模型评价
在分析该城区内不同区域重金属的污染程度问题时，本文对单因子指数法和内梅罗综合污染指数法进行合理的改进。考虑了当实测浓度值小背景值似的情况，得到新的单因子污染指数以及综合污染指数公式，得出不同区域中各种重金属的污染指数，以此数据得出该城区内不同区域重金属的污染程度。但是在评价某一区域重金属的综合污染程度时，没有考虑到不同重金属的毒性不同，则所产生的影响力也不同。
本文在建立模型确定污染源的位置时，通过分析重金属的传播途径，考虑到
由于分子运动，重金属本身也会向四周扩散，在地表的重金属主要随着地表水流动由高低是向低地势传播。以该城区污染最严重的重金属Hg和Cu为分析对象，根据重金属污染物的传播特性，建立一个二维对流-扩散方程。这样就避免了若单单考虑横向或纵向时，出现片面性，不能全方位的综合考虑重金属在整个城区中的一个随位置不同所产生的浓度变化。而二维对流-扩散方程就避免这种情况的发生。但是，由于附件中所提供的数据是某一时刻的数据。因此该模型是在假设该时刻重金属传播已达到稳定状态条件下建立的，这样模型得出的最后结果存在一定偏差。
6.2模型改进
由于资源和时间有限，本文无法查到不同重金属的毒性系数，若能将重金属的毒性系数考虑到综合评价污染指数的公式，这样会使得到的综合污染指数更具有说服力。
如果能够找到重金属扩散系数与其降解系数，会使二维对流-扩散方程得出的污染源位置更准确。
七．参考文献
[1]郝欣欣，慕金波，山东省辖沭河水系地表水环境容量研究，
https://www.doczj.com/doc/577465677.html,/p-167008196.html, 2011-9-11。
[2] 徐祖信，廖振良，张锦平，基于数学模型的苏州河上游和支流水质对干流水
质的影响分析，水动力学研究与进展，6(19)：733-743, 2004。
[3] 李云峰，王勇，杨林，沈锦仁，长江水质的评价和预测，工程数学学报，7(22)：
26-34，205。
[4] 董元彦，磷钾在红壤中的迁移及扩散系数，华中农业大学学报，3(14)：247-252，
1995。
[5] 尹骏柳云龙，上海市城郊土壤重金属空间分布及其污染评价，现代农业科
技，10(3)：251-255，2010。
[6] 冯仲科，空间数据的最佳内插法(Kriging法)及其在GIS中应用的构想，测
绘科技动态，3(4)：22-26, 1995。
附录一
图1 As含量分布图
图2 Cd含量分布图
图3 Cr含量分布图