城市表层土壤重金属污染分析
摘要
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。本文针对某城市八种重金属进行污染浓度分布建模分析，研究其不同地区受污染程度和污染原因，并确定污染源。
针对问题一，首先对附件一中数据进行处理，然后利用MATLAB 插值拟合曲线，得到八种重金属元素在该城市的空间分布图，运用单因子污染指数法和内梅罗综合污染指数法分别确定单种重金属元素的环境污染指数和元素的综合污染指数后进行污染等级评价，最总得出结论，在五大功能区中，工业区污染较为严重，山区污染污染程度不高，较清洁，且该城市受Z n 、H g 、Cu 元素污染较严重。 针对问题二，重金属元素污染之间存在必然的联系，同一区域的污染元素之间有较强的相关性，因此我们采用因子分析法来研究多个变量的相关性。通过对原始数据的处理，将成因进行归类，总结出几个比较客观的成因线索，最终导出成因结果。本文使用SPSS 软件对八种重金属元素进行相关性分析，最终得出结论：该城市的重金属污染主要原因为以下方面：工业污染源、交通污染源、燃煤污染源、农药污染源。
针对问题三，讨论重金属元素在土壤中的传播问题，我们建立抛物型的偏微分方程模型，解出这个偏微分方程的Cauchy 问题的解。由于污染源污染范围有限，所以我们挑选出污染较为严重的一处，通过反演的方法拟合出污染源的位置，本文仅针对Z n 、H g 、Cu 这类污染较严重的元素，得到一处或多处（x 0,y 0,z 0）
污染源。
针对问题四，我们首先提出模型的优缺点，然后根据优缺点提出模型的改进建议。考虑到各区域某一重金属元素污染下的地质累积污染程度强弱状况、地质演变因素，我们收集材料，在条件允许下，建立优化模型。例如，考虑地质累积指数法，建立更为精确的污染程度分析，使得问题一的结果更加可信；考虑大气干湿度对重金属土壤污染扩散的影响，建立相关数学模型，能够估计n 年后该城市的土壤污染程度。
关键词：单因子污染指数法、内梅罗综合污染指数法、因子相关性分析法、抛物型的偏微分方程模型
一、问题重述
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。
按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。
现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并用GPS 记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。
附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。
现要求你们通过数学建模来完成以下任务：
(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
(2) 通过数据分析，说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。 (4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？
二、模型假设
1、假设该城市各个功能区相互独立，不交叉。
2、假设该城市仅受这八种重金属污染，无其他污染物干扰。 3、假设八种重金属污染源之间无相互干扰。
4、假设八种金属污染物在分析期间分布浓度无明显变化。
三、符号说明
2,3,4,5,6,7,8 —— 环境质量指数 P i i=1，
2,3,4,5,6,7,8 —— 评价对象污染物i 的实测质量分数 C i i=1，
2,3,4,5,6,7,8 —— 污染物i 的评价标准，一般取二类标准 S i i=1，
—— 内梅罗综合污染指数 P 综
P imax —— 元素环境质量指数最大值 P iave —— 元素环境质量指数平均值
u (x , y , z , t ) —— t 时刻(x , y , z ) 处的一种重金属元素浓度
a 2, b 2, c 2 —— 分别为沿x , y , z 方向上的扩散系数
k 2 —— 衰减系数
(x i , y i , z i , u i ), i =1,2,..., n —— 取样点的横坐标、纵坐标、高度及一种元素的浓度
四、问题一 重金属污染分布及污染程度
4.1 问题分析
利用附件一给出的取样点位置以及所属功能区，附件二给出的八种重金属在对应位置的浓度。运用指数法计算计算出八种重金属元素在指数法条件下的各个指标值，从而确定区域污染程度。
初步分析计算得到如下指标：
表一:该城市八种重金属含量参数统计表
元素 背景值 最大值 最小值 平均值 标准差
As （μg /g ） Cd （μg /g ） C r （μg /g ） Cu （μg /g ） Hg （μg /g ） N i （μg /g ） （μg /g ）Pb Z n （μg /g ）
3.6 130 31 13.2 35 12.3 31 69
30.13 1.62 920.84 2528.48 16 142.5 472.48 3760.82
1.61 0.04 15.32 2.29 0.01 4.27 19.68 32.86
5.68 0.3 53.51 55.02 0.3 17.26 61.74 201.20
3.02 0.22 69.89 162.66 1.63 9.93 49.98 338.70
4.2 数据处理
1）首先将中金属浓度单位统一为μg /g 。
2）根据单因子污染指数法求出各重金属元素在不同浓度的环境质量指数
P i ，根据标准得出各个重金属元素质量评价等级。
3）利用P i 求得元素环境质量指数最大值P imax 和平均值P iave ，根据综合污染指数法的计算公式，求得内梅罗综合污染指数，确定综合元素质量评价等级。
4.3 模型建立
4.3.1 单因子污染指数法
单因子污染指数法是指在某一污染物影响下的环境污染指数。其能够反映污染物在该环境中的污染程度。
根据下式计算单因子污染程度等级，并进行分级。
P i =
C i
, i =1,2,3,4,5,6,7,8 （1）
S i
其中，国家土壤环境质量标准的二级标准经查阅为
表二:国家二级标准单因子污染评价各污染标准(mg /kg )
表三:城市表层土壤单因子污染分级标准
污染分指数 质量等级
P i
1≤P i
潜在汚染
2≤P i
P i ≥3 重污染
4.3.2 内梅罗综合污染指数法
前文对该城市单个元素的污染做了评价，但不能反映各元素共同作用对城市表层土壤的复合污染。鉴于此为突出环境要素中含量最大的污染物随环境质量的影响，因此采用内梅罗综合污染指数法进行综合评价，其公式为：
P 综（2） 内梅罗综合污染指数分级标准见下表：
表四:内梅罗综合污染指数分级标准
等级 Ⅰ Ⅱ Ⅲ Ⅳ
Ⅴ
综合指数
P 综≤0.7 1
3
污染等级 清洁 警戒限 轻度污染 重度污染 重污染
4.4 模型求解
根据附件一中给出的三维空间坐标进行插值处理，得到该城市功能区分布图。
图一 该城市功能区分布图
4.4.1 重金属元素的空间分布
依据附件所给数据，利用MATLAB 编程，得出八种重金属元素在该城市的空间分布等高线图，如下，
图二 As 的分布图 图三 Cd 的分布图
图四 Cr 的分布图 图五 Cu 的分布图
图六 Hg 的分布图 图七 Ni 的分布图
图八 Pb 的分布图 图九 Zn 的分布图
4.4.2 不同区域重金属污染程度分析 依据单因子污染指数公式，得出各重金属元素污染指数，应用单因子污染指数法对这八种重金属元素进行污染程度评价，得到下表：
从表七看出，该城市表层土壤污染并不严重，但有部分地区受Cd ，Cu ，Hg ，Zn 元素重污染，且Zn 重污染区最多达到18个。
下面给出不同功能区的受污染程度分析如下，
依据内梅罗综合污染指数法的计算公式得到五个功能区的内梅罗综合污染指数P （取平均值），建立表格： 综
依据上面几个表格得出结论，工业区污染状况最为严重，其次是生活区和交通区，然后是公园绿地区，而山区污染状况最好，污染程度都为清洁。
五、问题二 重金属污染主要原因
5.1 问题分析
问题二要求我们通过数据分析得出重金属污染的主要原因，由问题一得出八种元素在该城市空间分布以及不同功能区的污染程度。在不同功能区中，重金属的来源可能不同，但也可能来源于同一污染源，因此每处土壤中的污染元素存在一定的相关性，不互相独立，且相互影响。因此使用因子分析算法来解决这个问题，通过因子分析法并利用SPSS 软件各个因子空间分布，确定各因子集区域，从而得出污染主要原因。 5.2 模型建立及求解
通过对因子分析算法的了解, 我们发现该城市重金属元素含量的数据特征符合因子分析的要求。在该城市环境质量评价基础上，对存在污染的土壤的污染来源及分布进行分析。运用SPSS 22软件进行因子分析，得到如下结论：首先给出八种重金属元素的相关系数矩阵：
由表可见，Ni 和Cr 的相关性最好，相关系数最大为0.716。其次为Pb 和Cd ，相关系数为0.636。以下依次是Cu 和Cr 相关系数为0.509，而其它元素之间的相关性并不好。从成分上分析，相关信号的元素在成因和来源上有一定的联系。其原因分析：
1）工业生产中，镍铬合金是常用合金。
2）铅和镉元素在蓄电池的制造中会以气态形式散发。且他们都是工业排放的废物的主要元素。从问题一中各元素在不同功能区中的污染浓度看出，铅元素不仅在工业区，交通区和生活区都占有较大比重，说明汽车尾气污染也很严重。 3）绿地公园区的污染不大，其重金属遗留原因应是农药残留。
因子分析的关键在于利用相关系数矩阵求出相应的因子特征值和累计贡献率，依然可用用SPSS 22软件计算得出，结果如下：
从表中看出在累计方差达87.756%的前提下，分析得五个主因子，可以看出五个主因子为源资料提供了87.756%的信息，满足因子分析原则。
因子分析的主要目的是将具有相近的因子荷载的各个变量置于一个公因子之下，正交方差最大旋转使每一个主因子只与最少个数的变量有相关关系，而使足够多的因子负荷均很小，以便对因子的意义作出更合理的解释，输出结果见表。
旋轉元件矩陣
元件
a
As (μg/g) Cd (ng/g) Cr (μg/g) Cu (μg/g) Hg (ng/g) Ni (μg/g) Pb (μg/g) Zn (μg/g)
1 .452 .350 .846 .488 -.135 .883 .369 .551
2 .130 .679 .149 .591 .776 .088 .739 .431
根据软件得出结果，基于旋转理论:变量与某一个因子的联系系数绝对值(载荷) 越大, 则该因子与变量关系越近。正交因子解说明:因子1为Cr 、Ni 和Cu 的组合, 这说明这几种土壤重金属污染物可能是同一来源或相似来源; 因子2为Cd 和Pb 的组合表明两者可能有相似的来源。
Hg在土壤的污染严重, 污染区域大多处于某一工业区或某大型污染企业, 其来源比较单一。
表十四:中国土壤重金属污染来源表 来源 重金属
矿产开采、冶炼、加工排放的废气、废
Cr Hg As Pb Ni
水和废渣
煤和石油燃烧过程中排放的飘尘 Cr Hg As Pb
电镀工业废水 Cr Cd Ni Pb Cu Zn
塑料、电池、电子工业排放的废水 Hg Cd Pb Ni Zn
工业排放的废水 Hg
染料、化工制革工业排放的废水 Cr Cd
汽车尾气 Pb 农药、肥料 As Cu Cd
最后给出结论：该城市重金属污染来源主要有四个方面:交通污染源, 燃煤污染源, 工业污染源, 化肥农药污染源, 其中Zn 、Pb 、Cu 主要来自汽车尾气的排放, 即交通污染源;Hg 主要来自工业燃煤及居民采暖用煤;As 污染主要来自化工、冶金、炼焦、火力发电、造纸、玻璃、毛革、电子工业等工业污染；Cd 、Cr 、Ni 污染主要来自金属冶炼。
六、问题三 重金属污染源确定
6.1 问题分析
问题三要求我们确定重金属元素的污染源，因此我们要先了解污染物的传播特征，考虑到金属元素在土壤中并不固定，存在迁移和变化规律，建立迁移扩散模型。再此，我们建立非线性抛物型偏微方程作为模型来定性解决问题，考虑到扩散，采用衰减的扩散过程数学模型。随后利用反演的方法逐步找到污染源。
6.2 模型建立
6.2.1 有衰减的扩散过程数学模型
考虑到元素扩散和衰减的合作用，有物质不灭定理得：
222
∂u 2∂u 2∂u 2∂u 2
（1） =a +b +c -k u 222
∂t ∂x ∂y ∂z
方程（1）是常系数线性抛物方程，它就是有衰减的扩散过程的数学模型（对
于具体问题，尚需匹配求解）。在本题中就是要通过初始条件t =0时污染源的条件来进行求解。
6.2.2 线性抛物方程的Cauchy 问题求解
设污染源在点(x 0, y 0, z 0) 处，则此扩散模型满足Cauchy 问题：
222
⎧∂u 2∂u 2∂u 2∂u 2
=a +b +c -k u ⎪222
⎨ （2）∂t ∂x ∂y ∂z
⎪u (x , y , z ,0) =M δ(x -x ) δ(y -y ) δ(z -z )
000⎩
⎧0, x ≠0+∞
其中δ(x ) 为Dirac 函数，即δ(x ) =⎨，⎰δ(x ) dx =1，
⎩∞, x =0-∞解上述方程得到，
⎧(x -x 0) 2(y -y 0) 2(z -z 0) 2M 2⎫（3） u (x , y , z , t ) =exp ⎨----k t ⎬ 222
8πtabc *πt 4a t 4b t 4c t ⎩⎭
6.2.3 参数估计
在解（3）中存在四个未知参数，a 、b 、c 、k ，他们分别是扩散与衰减过程中的扩散系数和衰减系数的算术平方根。令取样时刻的t =1，对（3）式两边取对数得
⎧(x -x 0) 2(y -y 0) 2(z -z 0) 2M 2⎫ ln u (x , y , z , t ) =ln -⎨+++k t ⎬ （4）222
8πtabc *πt ⎩4a t 4b t 4c t ⎭
(x -x 0) 2(y -y 0) 2(z -z 0) 2111, Y =, Z =, α=-2, β=-2, γ=-2 令X =
444a b c
ε=ln
M 2-ln abc -k 3
(2π)
则（4）变为
W =ln u (x , y , z ,1) =αX +βY +γZ +ε （5）
利用已知数据(x i , y i , z i , u i ), i =1,2,..., n ，可求得(X , Y , Z , W ), i =1,2,..., n 的数据，用三元回归分析方法求出α、β、γ、ε的估值如下，
ε=W -（αX +βY +γZ ） (6)
1n 1n 1n 1n
其中W =∑W k , X =∑X k , Y =∑Y k , Z =∑Z k
n k =1n k =1n k =1n k =1
∧∧∧∧
α、β、γ满足方程
∧∧
⎧∧
⎪l 11α+l 12β+l 13γ=l 10
∧∧⎪∧
⎨l 21α+l 22β+l 23γ=l 20 （7）
⎪∧∧∧
⎪l 31α+l 32β+l 33γ=l 30⎩
n
n
n
∧∧∧
l 10=∑(X k -X )(W k -W ), l 20=∑(Y k -Y )(W k -W ), l 30=∑(Z k -Z )(W k -W )
k =1n
k =1
k =1
l 11=∑(X k -X ) , l 22=∑(Y k -Y ) , l 33=∑(Z k -Z ) 2 （8）
2
2
k =1
k =1
k =1
n n
l 12=∑(X k -X )(Y k -Y ), l 13=∑(Y k -Y )(Z k -Z ),
k =1n
k =1
n n
l 23=∑(Y k -Y )(Z k -Z ), l 12=l 21, l 31=l 13, l 32=l 23
k =1∧
1
β、γ可求得a 、b 、c 的估值，即利用α=-∧2, β=-∧2, γ=-∧2， 由α、
a b c
∧
∧
2
2
2
∧
∧
∧
11
又由于
k 2=ε+ln abc -ln
∧
∧
∧
∧
M
（9） (2π) 3
b c 带入（9）式得， 由（5）式可得ε，再把a 、、
k =ε+ln abc -ln
∧2
∧2
∧2
∧∧∧∧
M
（10） (2π) 3
∧2
∧2
b 、c 、k ，把他们带入（2）式分至此得到参数a 、b 、c 、k 的估计值a 、
2222
别代替a 2、b 2、c 2、k 2，则得不含未知数的解u (x , y , z , t ) 得近似表达式。 6.3 模型求解
通过第一问求解，我们知道该城市主要受Zn 、Cd 、Hg ，这里主要讨论Zn 元素的某一污染区域，讨论用反演法求污染源。从第一问中图九发现Zn 主要污染区域在{(x , y ) |x ∈[12000,15000],y ∈[8500,11000]}。
(12855,8945,18)、(13797,9621,18)、(14325,8666,23)、
1）筛选后区域内的点有(12641,9560,11)、 (14000,8970,14)、(14207,9980,14)、
(14065,10987,25)、(12734,10344,32)
2）初步选择(13797,9621,18)为污染源，将这个数据及t=1带入（3）式，这时u 是关于x ，y ，z 的方程且有参数a ，b ，c ，k ；
3）计算得a=425.73，b=472.8，c=16.655，k=23.56； 4）换一个假设污染源，重复上述步骤得误差最小。
最终得到污染源位置(13522.3,9448.0,23.9) 得出结论见下表
表十五:污染较严重重金属污染源表 元素 污染源坐标
（13522.3,9448.0,23.9）
Zn
（945.2,4562.9,69.3） （21456.4，11433.2，164.3）
Cd （2017.2,3112.3,12.3）
（15387.1,9190.0,25.0） （18522.8,3563.9,25.1）
Hg
（2501.1,3095.5,48.2）
七、问题四 模型优缺点
7.1 模型优缺点
7.1.1 优点：
1）问题一用单因子污染指数法对八种重金属的污染指数进行分析，得出单个元素的污染指数，并用内梅罗综合污染指数法求取了综合污染指数，对不同功能区的污染程度进行了总结，使得评价更客观，然后用MATLAB 曲线拟合方法画出了空间分布图，能更加直观看见金属分布状况。
2）问题2中运用了因子分析算法, 考虑到每个重金属元素污染成因并不是相互独立的, 运用SPSS 22软件建立了每个重金属元素的相关系数矩阵，通过主因子的属性确定重金属元素的污染原因。优点在于将每个重金属元素联系起来而不是单独分析每个元素的产生原因, 这样更贴近实际。
3）问题3中的抛物型偏微分方程形式很好, 每一个系数都有其代表的实际意义, 符合重金属污染物在土壤中的迁移规律。
7.1.2 缺点
1）城市土壤污染分析复杂，本文仅对表层土壤样品进行分析，其评价结果不能十分客观和完全真实地反映城市的重金属污染情况。
2）土壤环境质量评价应注意其历史变化过程, 不仅要对土壤环境质量现状作出评价, 也要对土壤环境质量回顾和土壤环境质量预测影响作出评价。本文仅对现状进行了分析。
3）问题3中, 利用反演的方法寻找污染源点肯定有所误差。
4）没有考虑自然地质变化过程及人类活动影响。 7.2 收集资料
地质环境演变
1）重金属元素不仅来自人类活动的影响, 同时大气干湿沉降是重金属的主要来源, 同时我们假设重金属大气干湿沉降速率此后保持不变, 在上述假设条件下, 我们需要知道各元素的年平均大气干湿沉降速率, 分别对重金属大气干湿沉降对表层土壤中重金属含量的累积影响进行分析，利用下列公式计算：
v*t
C t =C 0+
230
其中C t 表示t 年后表层重金属含量，C 0表示表层土壤重金属含量现状，v 表示重金属大气干湿沉降速率，t 表示预测年限。利用该计算公式可粗略估计t 年后表层土壤中重金属含量。
2）若给出了城市主要工厂的具体位置或该工厂的污染排放是否达标等信息, 则可以更为准确地分析该城市的污染原因及污染源的具体位置。
3）我们需要知道上述八种重金属元素在深层土壤中的浓度数据, 再利用问题一所述模型对各功能区深层土壤的重金属污染状况作物污染评价。
八、模型改进
1) 问题一中除了采用以上三种污染指数分析法外, 还可采用潜在生态危害指数法, 研究沉积物中重金属对环境的影响, 该种方法不仅反映了某一特定环境中的各种污染物的影响, 而且也反映了多种污染物的综合影响, 并且用定量的方法划分出潜在生态危害的程度。
2) 主因子分析的步骤:不仅需要对题目给出的所有数据进行主因子分析, 最好还应有针对各个功能区的主因子分析。
九、参考文献
[1] 姜启源, 谢金星, 叶俊, 《数学模型(第三版) 》, 北京:高等教育出版社,2003
附录
（1）城市取样点分布及等高线示意图 clear
jinshu=xlsread('wurannongdu.xls' ); %金属浓度 pos=xlsread('zuobiao.xls' ); %坐标 x=pos(:,1); y=pos(:,2); k=pos(:,4); z=jinshu(:,2);
[xx,yy]=meshgrid(0:500:28654,0:500:188449); zz=griddata(x,y,z,xx,yy,'v4' ); figure; subplot(2,1,1); hold on ; mesh(xx,yy,zz); hold on ;
plot3(x,y,z,'*');
axis([0 28654 0 18449 0 408]); subplot(2,1,2); contour(xx,yy,zz,20); hold on ;
zx=zz(1:end-1,2:end)-zz(1:end-1,1:end-1); zy=zz(2:end,1:end-1)-zz(1:end-1,1:end-1);
quiver(xx(1:end-1,1:end-1),yy(1:end-1,1:end-1),zx,zy);
（2）功能区分布图
[v,c]=voronoin([x,y]); figure; hold on ;
axis([0 28654 0 18449]); for i=1:length(c)
patch(v(c{i},1),v(c{i},2),k(i)); plot(x(i),y(i)); end
（3）重金属污染分布示意图（以As 为例） clear
jinshu=xlsread('wurannongdu.xls' ); %金属浓度 pos=xlsread('zuobiao.xls' ); %坐标 x=pos(:,1);
y=pos(:,2); z=jinshu(:,1);
[xx,yy]=meshgrid(0:500:28654,0:500:18449); zz=griddata(x,y,z,xx,yy,'v4' ); figure;
contour(xx,yy,zz,15); %hold on;
%plot(9328,4311,'*'); %hold on;
%plot(13797,9621,'+');
（4）问题三反演法求参数a b c k clear;
data=[12855 8945 18 102.86 13797 9621 18 3760.82 14325 8666 23 85.61 12641 9560 11 90.07 14000 8970 14 94.34 14207 9980 14 117.87]; x0=13797; y0=9621; z0=18;
M=1.60*1.0e+08;
X=((data(:,1)-x0).^2)/4; x=mean(X);
Y=((data(:,2)-y0).^2)/4; y=mean(Y);
Z=((data(:,3)-z0).^2)/4; z=mean(Z);
W=log(data(:,4)); w=mean(W);
l(1,1)=sum((X-x).^2); l(2,2)=sum((Y-y).^2); l(3,3)=sum((Z-z).^2); l(1,2)=sum((X-x).*(Y-y)); l(1,3)=sum((X-x).*(Z-z)); l(2,3)=sum((Y-y).*(Z-z));
l(2,1)=l(1,2); l(3,1)=l(1,3); l(3,2)=l(2,3); m(1)=sum((X-x).*(W-w)); m(2)=sum((Y-y).*(W-w)); m(3)=sum((Z-z).*(W-w)); zm=inv(l)*m'; a=(-1/zm(1)).^0.5; b=(-1/zm(2)).^0.5; c=(-1/zm(3)).^0.5;
ib=w-zm'*[x y z]';
k=(ib+log(a*b*c)-log(M/(2*pi^0.5)^3))^0.5; a b c k
（5）检验拟合参数a b c k准确性 function f=nongdu(x,y,z) a=534.59; b=324.72; c=2.89; k=2.21;
M=1.60*1.0e+08; x0=13797; y0=9621; z0=18;
f=M/(8*pi*a*b*c*(pi)^0.5)*exp(-(x-x0)^2/(4*a*a)-(y-y0)^2/(4*b*b)?-(z-z0)^2/(4*c*c)-k*k);
城市表层土壤重金属污染分析
摘要
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。本文针对某城市八种重金属进行污染浓度分布建模分析，研究其不同地区受污染程度和污染原因，并确定污染源。
针对问题一，首先对附件一中数据进行处理，然后利用MATLAB 插值拟合曲线，得到八种重金属元素在该城市的空间分布图，运用单因子污染指数法和内梅罗综合污染指数法分别确定单种重金属元素的环境污染指数和元素的综合污染指数后进行污染等级评价，最总得出结论，在五大功能区中，工业区污染较为严重，山区污染污染程度不高，较清洁，且该城市受Z n 、H g 、Cu 元素污染较严重。 针对问题二，重金属元素污染之间存在必然的联系，同一区域的污染元素之间有较强的相关性，因此我们采用因子分析法来研究多个变量的相关性。通过对原始数据的处理，将成因进行归类，总结出几个比较客观的成因线索，最终导出成因结果。本文使用SPSS 软件对八种重金属元素进行相关性分析，最终得出结论：该城市的重金属污染主要原因为以下方面：工业污染源、交通污染源、燃煤污染源、农药污染源。
针对问题三，讨论重金属元素在土壤中的传播问题，我们建立抛物型的偏微分方程模型，解出这个偏微分方程的Cauchy 问题的解。由于污染源污染范围有限，所以我们挑选出污染较为严重的一处，通过反演的方法拟合出污染源的位置，本文仅针对Z n 、H g 、Cu 这类污染较严重的元素，得到一处或多处（x 0,y 0,z 0）
污染源。
针对问题四，我们首先提出模型的优缺点，然后根据优缺点提出模型的改进建议。考虑到各区域某一重金属元素污染下的地质累积污染程度强弱状况、地质演变因素，我们收集材料，在条件允许下，建立优化模型。例如，考虑地质累积指数法，建立更为精确的污染程度分析，使得问题一的结果更加可信；考虑大气干湿度对重金属土壤污染扩散的影响，建立相关数学模型，能够估计n 年后该城市的土壤污染程度。
关键词：单因子污染指数法、内梅罗综合污染指数法、因子相关性分析法、抛物型的偏微分方程模型
一、问题重述
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。
按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同。
现对某城市城区土壤地质环境进行调查。为此，将所考察的城区划分为间距1公里左右的网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并用GPS 记录采样点的位置。应用专门仪器测试分析，获得了每个样本所含的多种化学元素的浓度数据。另一方面，按照2公里的间距在那些远离人群及工业活动的自然区取样，将其作为该城区表层土壤中元素的背景值。
附件1列出了采样点的位置、海拔高度及其所属功能区等信息，附件2列出了8种主要重金属元素在采样点处的浓度，附件3列出了8种主要重金属元素的背景值。
现要求你们通过数学建模来完成以下任务：
(1) 给出8种主要重金属元素在该城区的空间分布，并分析该城区内不同区域重金属的污染程度。
(2) 通过数据分析，说明重金属污染的主要原因。
(3) 分析重金属污染物的传播特征，由此建立模型，确定污染源的位置。 (4) 分析你所建立模型的优缺点，为更好地研究城市地质环境的演变模式，还应收集什么信息？有了这些信息，如何建立模型解决问题？
二、模型假设
1、假设该城市各个功能区相互独立，不交叉。
2、假设该城市仅受这八种重金属污染，无其他污染物干扰。 3、假设八种重金属污染源之间无相互干扰。
4、假设八种金属污染物在分析期间分布浓度无明显变化。
三、符号说明
2,3,4,5,6,7,8 —— 环境质量指数 P i i=1，
2,3,4,5,6,7,8 —— 评价对象污染物i 的实测质量分数 C i i=1，
2,3,4,5,6,7,8 —— 污染物i 的评价标准，一般取二类标准 S i i=1，
—— 内梅罗综合污染指数 P 综
P imax —— 元素环境质量指数最大值 P iave —— 元素环境质量指数平均值
u (x , y , z , t ) —— t 时刻(x , y , z ) 处的一种重金属元素浓度
a 2, b 2, c 2 —— 分别为沿x , y , z 方向上的扩散系数
k 2 —— 衰减系数
(x i , y i , z i , u i ), i =1,2,..., n —— 取样点的横坐标、纵坐标、高度及一种元素的浓度
四、问题一 重金属污染分布及污染程度
4.1 问题分析
利用附件一给出的取样点位置以及所属功能区，附件二给出的八种重金属在对应位置的浓度。运用指数法计算计算出八种重金属元素在指数法条件下的各个指标值，从而确定区域污染程度。
初步分析计算得到如下指标：
表一:该城市八种重金属含量参数统计表
元素 背景值 最大值 最小值 平均值 标准差
As （μg /g ） Cd （μg /g ） C r （μg /g ） Cu （μg /g ） Hg （μg /g ） N i （μg /g ） （μg /g ）Pb Z n （μg /g ）
3.6 130 31 13.2 35 12.3 31 69
30.13 1.62 920.84 2528.48 16 142.5 472.48 3760.82
1.61 0.04 15.32 2.29 0.01 4.27 19.68 32.86
5.68 0.3 53.51 55.02 0.3 17.26 61.74 201.20
3.02 0.22 69.89 162.66 1.63 9.93 49.98 338.70
4.2 数据处理
1）首先将中金属浓度单位统一为μg /g 。
2）根据单因子污染指数法求出各重金属元素在不同浓度的环境质量指数
P i ，根据标准得出各个重金属元素质量评价等级。
3）利用P i 求得元素环境质量指数最大值P imax 和平均值P iave ，根据综合污染指数法的计算公式，求得内梅罗综合污染指数，确定综合元素质量评价等级。
4.3 模型建立
4.3.1 单因子污染指数法
单因子污染指数法是指在某一污染物影响下的环境污染指数。其能够反映污染物在该环境中的污染程度。
根据下式计算单因子污染程度等级，并进行分级。
P i =
C i
, i =1,2,3,4,5,6,7,8 （1）
S i
其中，国家土壤环境质量标准的二级标准经查阅为
表二:国家二级标准单因子污染评价各污染标准(mg /kg )
表三:城市表层土壤单因子污染分级标准
污染分指数 质量等级
P i
1≤P i
潜在汚染
2≤P i
P i ≥3 重污染
4.3.2 内梅罗综合污染指数法
前文对该城市单个元素的污染做了评价，但不能反映各元素共同作用对城市表层土壤的复合污染。鉴于此为突出环境要素中含量最大的污染物随环境质量的影响，因此采用内梅罗综合污染指数法进行综合评价，其公式为：
P 综（2） 内梅罗综合污染指数分级标准见下表：
表四:内梅罗综合污染指数分级标准
等级 Ⅰ Ⅱ Ⅲ Ⅳ
Ⅴ
综合指数
P 综≤0.7 1
3
污染等级 清洁 警戒限 轻度污染 重度污染 重污染
4.4 模型求解
根据附件一中给出的三维空间坐标进行插值处理，得到该城市功能区分布图。
图一 该城市功能区分布图
4.4.1 重金属元素的空间分布
依据附件所给数据，利用MATLAB 编程，得出八种重金属元素在该城市的空间分布等高线图，如下，
图二 As 的分布图 图三 Cd 的分布图
图四 Cr 的分布图 图五 Cu 的分布图
图六 Hg 的分布图 图七 Ni 的分布图
图八 Pb 的分布图 图九 Zn 的分布图
4.4.2 不同区域重金属污染程度分析 依据单因子污染指数公式，得出各重金属元素污染指数，应用单因子污染指数法对这八种重金属元素进行污染程度评价，得到下表：
从表七看出，该城市表层土壤污染并不严重，但有部分地区受Cd ，Cu ，Hg ，Zn 元素重污染，且Zn 重污染区最多达到18个。
下面给出不同功能区的受污染程度分析如下，
依据内梅罗综合污染指数法的计算公式得到五个功能区的内梅罗综合污染指数P （取平均值），建立表格： 综
依据上面几个表格得出结论，工业区污染状况最为严重，其次是生活区和交通区，然后是公园绿地区，而山区污染状况最好，污染程度都为清洁。
五、问题二 重金属污染主要原因
5.1 问题分析
问题二要求我们通过数据分析得出重金属污染的主要原因，由问题一得出八种元素在该城市空间分布以及不同功能区的污染程度。在不同功能区中，重金属的来源可能不同，但也可能来源于同一污染源，因此每处土壤中的污染元素存在一定的相关性，不互相独立，且相互影响。因此使用因子分析算法来解决这个问题，通过因子分析法并利用SPSS 软件各个因子空间分布，确定各因子集区域，从而得出污染主要原因。 5.2 模型建立及求解
通过对因子分析算法的了解, 我们发现该城市重金属元素含量的数据特征符合因子分析的要求。在该城市环境质量评价基础上，对存在污染的土壤的污染来源及分布进行分析。运用SPSS 22软件进行因子分析，得到如下结论：首先给出八种重金属元素的相关系数矩阵：
由表可见，Ni 和Cr 的相关性最好，相关系数最大为0.716。其次为Pb 和Cd ，相关系数为0.636。以下依次是Cu 和Cr 相关系数为0.509，而其它元素之间的相关性并不好。从成分上分析，相关信号的元素在成因和来源上有一定的联系。其原因分析：
1）工业生产中，镍铬合金是常用合金。
2）铅和镉元素在蓄电池的制造中会以气态形式散发。且他们都是工业排放的废物的主要元素。从问题一中各元素在不同功能区中的污染浓度看出，铅元素不仅在工业区，交通区和生活区都占有较大比重，说明汽车尾气污染也很严重。 3）绿地公园区的污染不大，其重金属遗留原因应是农药残留。
因子分析的关键在于利用相关系数矩阵求出相应的因子特征值和累计贡献率，依然可用用SPSS 22软件计算得出，结果如下：
从表中看出在累计方差达87.756%的前提下，分析得五个主因子，可以看出五个主因子为源资料提供了87.756%的信息，满足因子分析原则。
因子分析的主要目的是将具有相近的因子荷载的各个变量置于一个公因子之下，正交方差最大旋转使每一个主因子只与最少个数的变量有相关关系，而使足够多的因子负荷均很小，以便对因子的意义作出更合理的解释，输出结果见表。
旋轉元件矩陣
元件
a
As (μg/g) Cd (ng/g) Cr (μg/g) Cu (μg/g) Hg (ng/g) Ni (μg/g) Pb (μg/g) Zn (μg/g)
1 .452 .350 .846 .488 -.135 .883 .369 .551
2 .130 .679 .149 .591 .776 .088 .739 .431
根据软件得出结果，基于旋转理论:变量与某一个因子的联系系数绝对值(载荷) 越大, 则该因子与变量关系越近。正交因子解说明:因子1为Cr 、Ni 和Cu 的组合, 这说明这几种土壤重金属污染物可能是同一来源或相似来源; 因子2为Cd 和Pb 的组合表明两者可能有相似的来源。
Hg在土壤的污染严重, 污染区域大多处于某一工业区或某大型污染企业, 其来源比较单一。
表十四:中国土壤重金属污染来源表 来源 重金属
矿产开采、冶炼、加工排放的废气、废
Cr Hg As Pb Ni
水和废渣
煤和石油燃烧过程中排放的飘尘 Cr Hg As Pb
电镀工业废水 Cr Cd Ni Pb Cu Zn
塑料、电池、电子工业排放的废水 Hg Cd Pb Ni Zn
工业排放的废水 Hg
染料、化工制革工业排放的废水 Cr Cd
汽车尾气 Pb 农药、肥料 As Cu Cd
最后给出结论：该城市重金属污染来源主要有四个方面:交通污染源, 燃煤污染源, 工业污染源, 化肥农药污染源, 其中Zn 、Pb 、Cu 主要来自汽车尾气的排放, 即交通污染源;Hg 主要来自工业燃煤及居民采暖用煤;As 污染主要来自化工、冶金、炼焦、火力发电、造纸、玻璃、毛革、电子工业等工业污染；Cd 、Cr 、Ni 污染主要来自金属冶炼。
六、问题三 重金属污染源确定
6.1 问题分析
问题三要求我们确定重金属元素的污染源，因此我们要先了解污染物的传播特征，考虑到金属元素在土壤中并不固定，存在迁移和变化规律，建立迁移扩散模型。再此，我们建立非线性抛物型偏微方程作为模型来定性解决问题，考虑到扩散，采用衰减的扩散过程数学模型。随后利用反演的方法逐步找到污染源。
6.2 模型建立
6.2.1 有衰减的扩散过程数学模型
考虑到元素扩散和衰减的合作用，有物质不灭定理得：
222
∂u 2∂u 2∂u 2∂u 2
（1） =a +b +c -k u 222
∂t ∂x ∂y ∂z
方程（1）是常系数线性抛物方程，它就是有衰减的扩散过程的数学模型（对
于具体问题，尚需匹配求解）。在本题中就是要通过初始条件t =0时污染源的条件来进行求解。
6.2.2 线性抛物方程的Cauchy 问题求解
设污染源在点(x 0, y 0, z 0) 处，则此扩散模型满足Cauchy 问题：
222
⎧∂u 2∂u 2∂u 2∂u 2
=a +b +c -k u ⎪222
⎨ （2）∂t ∂x ∂y ∂z
⎪u (x , y , z ,0) =M δ(x -x ) δ(y -y ) δ(z -z )
000⎩
⎧0, x ≠0+∞
其中δ(x ) 为Dirac 函数，即δ(x ) =⎨，⎰δ(x ) dx =1，
⎩∞, x =0-∞解上述方程得到，
⎧(x -x 0) 2(y -y 0) 2(z -z 0) 2M 2⎫（3） u (x , y , z , t ) =exp ⎨----k t ⎬ 222
8πtabc *πt 4a t 4b t 4c t ⎩⎭
6.2.3 参数估计
在解（3）中存在四个未知参数，a 、b 、c 、k ，他们分别是扩散与衰减过程中的扩散系数和衰减系数的算术平方根。令取样时刻的t =1，对（3）式两边取对数得
⎧(x -x 0) 2(y -y 0) 2(z -z 0) 2M 2⎫ ln u (x , y , z , t ) =ln -⎨+++k t ⎬ （4）222
8πtabc *πt ⎩4a t 4b t 4c t ⎭
(x -x 0) 2(y -y 0) 2(z -z 0) 2111, Y =, Z =, α=-2, β=-2, γ=-2 令X =
444a b c
ε=ln
M 2-ln abc -k 3
(2π)
则（4）变为
W =ln u (x , y , z ,1) =αX +βY +γZ +ε （5）
利用已知数据(x i , y i , z i , u i ), i =1,2,..., n ，可求得(X , Y , Z , W ), i =1,2,..., n 的数据，用三元回归分析方法求出α、β、γ、ε的估值如下，
ε=W -（αX +βY +γZ ） (6)
1n 1n 1n 1n
其中W =∑W k , X =∑X k , Y =∑Y k , Z =∑Z k
n k =1n k =1n k =1n k =1
∧∧∧∧
α、β、γ满足方程
∧∧
⎧∧
⎪l 11α+l 12β+l 13γ=l 10
∧∧⎪∧
⎨l 21α+l 22β+l 23γ=l 20 （7）
⎪∧∧∧
⎪l 31α+l 32β+l 33γ=l 30⎩
n
n
n
∧∧∧
l 10=∑(X k -X )(W k -W ), l 20=∑(Y k -Y )(W k -W ), l 30=∑(Z k -Z )(W k -W )
k =1n
k =1
k =1
l 11=∑(X k -X ) , l 22=∑(Y k -Y ) , l 33=∑(Z k -Z ) 2 （8）
2
2
k =1
k =1
k =1
n n
l 12=∑(X k -X )(Y k -Y ), l 13=∑(Y k -Y )(Z k -Z ),
k =1n
k =1
n n
l 23=∑(Y k -Y )(Z k -Z ), l 12=l 21, l 31=l 13, l 32=l 23
k =1∧
1
β、γ可求得a 、b 、c 的估值，即利用α=-∧2, β=-∧2, γ=-∧2， 由α、
a b c
∧
∧
2
2
2
∧
∧
∧
11
又由于
k 2=ε+ln abc -ln
∧
∧
∧
∧
M
（9） (2π) 3
b c 带入（9）式得， 由（5）式可得ε，再把a 、、
k =ε+ln abc -ln
∧2
∧2
∧2
∧∧∧∧
M
（10） (2π) 3
∧2
∧2
b 、c 、k ，把他们带入（2）式分至此得到参数a 、b 、c 、k 的估计值a 、
2222
别代替a 2、b 2、c 2、k 2，则得不含未知数的解u (x , y , z , t ) 得近似表达式。 6.3 模型求解
通过第一问求解，我们知道该城市主要受Zn 、Cd 、Hg ，这里主要讨论Zn 元素的某一污染区域，讨论用反演法求污染源。从第一问中图九发现Zn 主要污染区域在{(x , y ) |x ∈[12000,15000],y ∈[8500,11000]}。
(12855,8945,18)、(13797,9621,18)、(14325,8666,23)、
1）筛选后区域内的点有(12641,9560,11)、 (14000,8970,14)、(14207,9980,14)、
(14065,10987,25)、(12734,10344,32)
2）初步选择(13797,9621,18)为污染源，将这个数据及t=1带入（3）式，这时u 是关于x ，y ，z 的方程且有参数a ，b ，c ，k ；
3）计算得a=425.73，b=472.8，c=16.655，k=23.56； 4）换一个假设污染源，重复上述步骤得误差最小。
最终得到污染源位置(13522.3,9448.0,23.9) 得出结论见下表
表十五:污染较严重重金属污染源表 元素 污染源坐标
（13522.3,9448.0,23.9）
Zn
（945.2,4562.9,69.3） （21456.4，11433.2，164.3）
Cd （2017.2,3112.3,12.3）
（15387.1,9190.0,25.0） （18522.8,3563.9,25.1）
Hg
（2501.1,3095.5,48.2）
七、问题四 模型优缺点
7.1 模型优缺点
7.1.1 优点：
1）问题一用单因子污染指数法对八种重金属的污染指数进行分析，得出单个元素的污染指数，并用内梅罗综合污染指数法求取了综合污染指数，对不同功能区的污染程度进行了总结，使得评价更客观，然后用MATLAB 曲线拟合方法画出了空间分布图，能更加直观看见金属分布状况。
2）问题2中运用了因子分析算法, 考虑到每个重金属元素污染成因并不是相互独立的, 运用SPSS 22软件建立了每个重金属元素的相关系数矩阵，通过主因子的属性确定重金属元素的污染原因。优点在于将每个重金属元素联系起来而不是单独分析每个元素的产生原因, 这样更贴近实际。
3）问题3中的抛物型偏微分方程形式很好, 每一个系数都有其代表的实际意义, 符合重金属污染物在土壤中的迁移规律。
7.1.2 缺点
1）城市土壤污染分析复杂，本文仅对表层土壤样品进行分析，其评价结果不能十分客观和完全真实地反映城市的重金属污染情况。
2）土壤环境质量评价应注意其历史变化过程, 不仅要对土壤环境质量现状作出评价, 也要对土壤环境质量回顾和土壤环境质量预测影响作出评价。本文仅对现状进行了分析。
3）问题3中, 利用反演的方法寻找污染源点肯定有所误差。
4）没有考虑自然地质变化过程及人类活动影响。 7.2 收集资料
地质环境演变
1）重金属元素不仅来自人类活动的影响, 同时大气干湿沉降是重金属的主要来源, 同时我们假设重金属大气干湿沉降速率此后保持不变, 在上述假设条件下, 我们需要知道各元素的年平均大气干湿沉降速率, 分别对重金属大气干湿沉降对表层土壤中重金属含量的累积影响进行分析，利用下列公式计算：
v*t
C t =C 0+
230
其中C t 表示t 年后表层重金属含量，C 0表示表层土壤重金属含量现状，v 表示重金属大气干湿沉降速率，t 表示预测年限。利用该计算公式可粗略估计t 年后表层土壤中重金属含量。
2）若给出了城市主要工厂的具体位置或该工厂的污染排放是否达标等信息, 则可以更为准确地分析该城市的污染原因及污染源的具体位置。
3）我们需要知道上述八种重金属元素在深层土壤中的浓度数据, 再利用问题一所述模型对各功能区深层土壤的重金属污染状况作物污染评价。
八、模型改进
1) 问题一中除了采用以上三种污染指数分析法外, 还可采用潜在生态危害指数法, 研究沉积物中重金属对环境的影响, 该种方法不仅反映了某一特定环境中的各种污染物的影响, 而且也反映了多种污染物的综合影响, 并且用定量的方法划分出潜在生态危害的程度。
2) 主因子分析的步骤:不仅需要对题目给出的所有数据进行主因子分析, 最好还应有针对各个功能区的主因子分析。
九、参考文献
[1] 姜启源, 谢金星, 叶俊, 《数学模型(第三版) 》, 北京:高等教育出版社,2003
附录
（1）城市取样点分布及等高线示意图 clear
jinshu=xlsread('wurannongdu.xls' ); %金属浓度 pos=xlsread('zuobiao.xls' ); %坐标 x=pos(:,1); y=pos(:,2); k=pos(:,4); z=jinshu(:,2);
[xx,yy]=meshgrid(0:500:28654,0:500:188449); zz=griddata(x,y,z,xx,yy,'v4' ); figure; subplot(2,1,1); hold on ; mesh(xx,yy,zz); hold on ;
plot3(x,y,z,'*');
axis([0 28654 0 18449 0 408]); subplot(2,1,2); contour(xx,yy,zz,20); hold on ;
zx=zz(1:end-1,2:end)-zz(1:end-1,1:end-1); zy=zz(2:end,1:end-1)-zz(1:end-1,1:end-1);
quiver(xx(1:end-1,1:end-1),yy(1:end-1,1:end-1),zx,zy);
（2）功能区分布图
[v,c]=voronoin([x,y]); figure; hold on ;
axis([0 28654 0 18449]); for i=1:length(c)
patch(v(c{i},1),v(c{i},2),k(i)); plot(x(i),y(i)); end
（3）重金属污染分布示意图（以As 为例） clear
jinshu=xlsread('wurannongdu.xls' ); %金属浓度 pos=xlsread('zuobiao.xls' ); %坐标 x=pos(:,1);
y=pos(:,2); z=jinshu(:,1);
[xx,yy]=meshgrid(0:500:28654,0:500:18449); zz=griddata(x,y,z,xx,yy,'v4' ); figure;
contour(xx,yy,zz,15); %hold on;
%plot(9328,4311,'*'); %hold on;
%plot(13797,9621,'+');
（4）问题三反演法求参数a b c k clear;
data=[12855 8945 18 102.86 13797 9621 18 3760.82 14325 8666 23 85.61 12641 9560 11 90.07 14000 8970 14 94.34 14207 9980 14 117.87]; x0=13797; y0=9621; z0=18;
M=1.60*1.0e+08;
X=((data(:,1)-x0).^2)/4; x=mean(X);
Y=((data(:,2)-y0).^2)/4; y=mean(Y);
Z=((data(:,3)-z0).^2)/4; z=mean(Z);
W=log(data(:,4)); w=mean(W);
l(1,1)=sum((X-x).^2); l(2,2)=sum((Y-y).^2); l(3,3)=sum((Z-z).^2); l(1,2)=sum((X-x).*(Y-y)); l(1,3)=sum((X-x).*(Z-z)); l(2,3)=sum((Y-y).*(Z-z));
l(2,1)=l(1,2); l(3,1)=l(1,3); l(3,2)=l(2,3); m(1)=sum((X-x).*(W-w)); m(2)=sum((Y-y).*(W-w)); m(3)=sum((Z-z).*(W-w)); zm=inv(l)*m'; a=(-1/zm(1)).^0.5; b=(-1/zm(2)).^0.5; c=(-1/zm(3)).^0.5;
ib=w-zm'*[x y z]';
k=(ib+log(a*b*c)-log(M/(2*pi^0.5)^3))^0.5; a b c k
（5）检验拟合参数a b c k准确性 function f=nongdu(x,y,z) a=534.59; b=324.72; c=2.89; k=2.21;
M=1.60*1.0e+08; x0=13797; y0=9621; z0=18;
f=M/(8*pi*a*b*c*(pi)^0.5)*exp(-(x-x0)^2/(4*a*a)-(y-y0)^2/(4*b*b)?-(z-z0)^2/(4*c*c)-k*k);