As Cd Cr Cu Ni Pb Zn 1 一级 三级 二级 二级 二级 二级 三级 2 二级 三级 二级 四级 五级 三级 四级 3 一级 一级 一级 一级 三级 二级 一级 4 一级 三级 二级 二级 五级 二级 三级 5 一级 三级 一级 二级 二级 二级 二级 与此同时，我们在附录中还提供了一种单项污染指数模型的修正模型，希望对该模型的建立有所帮助。 问题二
对于问题二，意在说明重金属污染的主要原因。 4.2.1基于因子分子法判断重金属污染的主要原因
多变量研究中，由于变量的个数很多，并且彼此往往存在一定的相关性，因此使观察的数据反映的信息在一定程度上重叠。因子分析则是通过一种降维方法进行简化得到综合指标。综合指标之间既互不相关，又能反映原来的观察指标的信息。
因子分析从变量的相关矩阵出发将一个m 维的随机向量X分解成低于m个且有代表性的公因子和一个特殊的m维向量，使其公因子数取得最佳的个数，从而使对m维随机向量的研究转化成对较少个数的公因子的研究。[3]
题中土壤单点样重金属元素含量的数据特征完全符合因子分析的要求，在这里以Hg、Cd、Pb、As、Cu、Cr、Ni、Zn八种重金属元素指标作因子分析，这样在解释各指标变化异常时可以着重讨论综合指标因子，同时为题中重金属污染成因的解释提供一定的理论依据。下面对各功能区土壤采点样重金属元素含量的数据标准化处理后，经 for windows统计软件进行因子分析，可得出以下结果。 1）给出题中表层土壤Hg、Cd、Pb、As、Cu、Cr、Ni、Zn八种重金属原始含量数据的相关系数矩阵。
表8 相关系数矩阵
表中可见，Ni和Cr的相关系数最大，为，相关性最好，从成因上来分析，相关性较好的元素可能在成因和来源上有一定的关联。其原因肯能是： 1.工业生产的相互关联性，镍和铬合金是工业上常用的合金。
2.考虑到表7中山区也含有很高浓度的镍，可以推断该城区富含镍矿，而对镍矿的开采和镍铬合金的冶炼是其污染的主要原因。
其次为Pb和Cd，相关系数为0．660，这确实出乎我们的意料，因为一般镍镉合金出现在一起比较多，但了解了Pb和Cd的物理性质后，我们不难发现这两种金属的一些共同点： 1.它们都广泛用于蓄电池的制造且易于以气态的方式进行传播
2.它们都是工业废水的主要组成部分[4] 。所以推断其污染原因是工业废水的排放。另外，我们从表7中可以看到Pb的污染等级很平均，那是因为像铅在石油中含量很大，故在工业区中其污染原因是石油能源的燃烧与化工产品的生产；而汽车也以石油为动力源，故在交通区，其也有广泛的分布；而化工产品的使用，生活区的铅含量有时会很大，这在附件数据中也有反映，铅的最大浓度出现在生活区。
以下依次是Cr和Cu，Pb和Cu，相关系数分别为，，其它元素之间的相关性并不是很好。
2）利用相关系数矩阵求出相应的因子的特征值和累计贡献率，见表9
表9 特征值和累计贡献率
从表中可见，在累积方差为％ (>85％ )的前提下，分析得5个主因子，可以看到5个主因子提供了源资料％的信息，满足因子分析的原则。从上表可以看出旋转前后总的累计
贡献率有发生变化，即总的信息量没有损失。而且，表中显示旋转之后，主因子1和主因子2的方差贡献率均为22％左右，主因子2到主因子5的方差贡献率的范围为％到％之间。这可以解释为因子1和因子2可能为本市土壤重金属污染的最重要的污染源，从而推测出金属污染的原因主要是本市重金属污染的贡献最大，因子3、因子4、因子5对本市重金属污染有重要作用。而从上面分析和下面的因子组成可知因子1和因子2主要是工业污染和交通废气污染，故更加确定了污染原因。 3）因子的相关分析
因子分析的主要目的是将具有相近的因子荷载的各个变量置于一个公因子之下，正交方差最大旋转使每一个主因子只与最少个数的变量有相关关系，而使足够多的因子负荷均很小，以便对因子的意义作出更合理的解释。输出结果见表。
表10正交方差最大旋转后的矩阵
基于旋转理论：变量与某一个因子的联系系数绝对值(载荷)越大，则该因子与变量关系越近。正交因子解说明：因子1为Cr、Ni和Cu的组合，这说明这几种土壤重金属污染物可能是同一来源或相似来源；因子2为Cd和Pb的组合表明两者可能有相似的来源；因子3为Hg；因子4为As；因子5为Zn。
这里，我们对Hg的污染原因作出分析。Hg在土壤的污染严重，污染区域大多处于某一工业区或某大型污染企业，其来源比较单一，废气和废水是其污染的主要来源。[5]
表11中国土壤重金属污染物来源[6] 来源
矿产开采、冶炼、加工排放的废气、废水和废渣
煤和石油燃烧过程中排放的飘尘
电镀工业废水
塑料、电池、电子工业排放的废水
Hg工业排放的废水
染料、化工制革工业排放的废水
汽车尾气 农药、化肥
重金属
Cr Hg As Pb Ni Mo
Cr Hg As Pb Cr Cd Ni Pb Cu Zn Hg Cd Pb Ni Zn
Hg Cr Cd Pb As Cu Cd
问题三 4.3.1 模型综述
对于问题的第三问我们主要采取了两个模型：
一是以全体坐标为研究对象，通过假设污染源来进行比较和检索，从而确定污染源的范围。
二是以第一问中的图像和污染源范围为基础，利用已经筛选了的数据，基于最小二乘法的原理，作出其污染源的目标函数，直接求出污染源的空间坐标。 4.3.2 模型一
根据第一问所得出的8种主要重金属元素在该城区的空间分布图，以及统计得出的八种元素的浓度排序，首先在元素浓度集中分布区选择一个浓度较高点L(xi,yi,zi,qi)，其中
x,y,z为其空间坐标，q为其所在功能区，由于题目中所提到的采样特点，我们可以找到其相邻的两个采样点L(xi?1,yi?1,zi?1,qi?1)和 L(xi?1,yi?1,zi?1,qi?1)，我们根据污染物扩散的基本原理，认为一般情况下：污染物浓度与距污染源的空间距离成反比关系，反映到所选的浓度最高点上，即：
我们在这里设这个最初距离为Ri?1。
所以我们要做的第一步是：若现实情况下Ci?1?Ci，我们则保留L(xi,yi,zi,qi)，若不然，我们就用Ci?1中较大的点代替L(xi,yi,zi,qi)，并重复上面的比较。
接着我们进行第二步，将Ri?1中的较大值赋值给R作为基准值，其编号赋值给m随后顺藤摸瓜，找到另一个与L(xi?1,yi?1,zi?1,qi?1)相邻的采样点L(xm?1,ym?1,zm?1,pm?1)，计算它与
L(xi,yi,zi,qi)的距离r。
在第三步中，我们对R与r进行比较，若R?r：就对Cm与Cm?1做比较，若Cm>Cm?1，则将m?1赋值给m，r的值赋值给R，继续顺藤摸瓜。若Cm
m?1赋值给m，r的值赋值给R，继续顺藤摸瓜，若不然则重复第一步中有关置换浓度最高点的方法。
第四步中，我们要对这个循环算法选一个终止条件，在参照了第一问图中污染物聚集区的大体范围后，我们以i?m?3，即以一个在其相邻上下三个采样点中浓度保持递减的采样点为污染源范围基准点。同时，我们以rr为半径来划定污染源的范围。 下面，我们用一个流程图来形象的描述我们的求解过程： 图3 流程图（见下页）
由于要输投入的数据量较大，在这里，我们仅以As为例作出其污染源范围，相关程序附于附件上。所得到的污染源范围是以（74，781，5）为原点，以1042m为半径的区域范围。
4.3.3 模型二
我们建立这个模型的基本思路是：依照各元素的浓度与坐标点的关系，利用最小二乘法拟合曲面方程。同时，通过对所求的曲面方程求极值，得到的极值点一方面可以在空间上反映重金属元素的分布，另一方面可结合极值点附近的元素浓度分布，来判断污染源的位置。
1.曲面方程的规范化形式及其系数
本题中，我们采用二次曲面方程表示各浓度与坐标点的关系。二次曲面方程可表示为[7]
Z?C20x2?C11xy?C02y2?C10x?C01y?C00??(?Ci,j?ixiyj?i)，
j?2i?j002.曲面方程的拟合
（xk,yk,zk）对于给定的一组数据点，k=1,2,…,N，求作m次(m=N)多项式
使总误差
为最小。
这里Q可看作是关于Cij?(i,j?0,1,?,m)的函数，共有(m+1)(m+2)/2个未知量Cij,所以上述拟合多项式的构造问题可以归结为求多元函数的极值问题，即Ci,j-i应满足而得
N?0?N??IJ?II?iyJ?j?I?i???Ci,j-ixzx??????kkykk化简并交换求和符号后得，k???k?1?j?m?i?j?k?10[8]
?Q?0，因?Ci,j?i
本题中以Cu元素为例，在其空间分布图的基础上，缩小范围拟合其浓度与坐标的关系。
表 12 Cu元素集中分布数组 编号 x(m) y(m) 海拔(m) Cu (μg/g) 13 2427 3971 2 21 2486 5999 2 7 2883 3617 15 150 4020 2990 27 23 3573 6213 5 14 3526 4357 7 6 1647 2728 6 22 3299 6018 4 8 2383 3692 7 利用最小二乘法求二次曲面方程的过程如下（X 、Y单位取km）： （1）求出坐标x、y立方、四次方对应的值（见附录）； （2）计算z 与坐标x、y平方的乘积（见附录）；
（3）列出包含6个未知系数的方程组，利用matlab求解出的系数为： C=+003 *
即所拟合的二次曲面方程为：
Z?1.0e?003*(?0.7776x2?0.2040xy?0.1115y2?3.5721x?1.8003y?0.0071)
对于污染源，我们可以考虑对曲面方程求极大值，相当于对
Z'?1.0e?003*(0.7776x2?0.2040xy?0.1115y2?3.5721x?1.8003y?0.0071)
求极小值。
也即minZ''?0.7776x2?0.2040xy?0.1115y2?3.5721x?1.8003y （忽略了共同比例和常数项） 问题四
4.4.1模型的优缺点
在此我们针对第三问中的模型进行评判。 对于模型一，其优点是： 1)层次分明，易于理解； 2）覆盖的数据全面；
3）考虑因素全面，比如考虑进了海拔的影响。
其缺点是：
1）数据量大，输入繁琐枯燥；
2）污染源表示为以采样点为中心的圆形区域，无法精确到点。 对于模型二，其优点是： 1）数学模型清晰；
2）可以将污染源确定到点，进一步缩小了污染源范围； 3）数学原理明确；
4）需要处理的数据量小。 其缺点是：
1）所取的拟合点有限，且模型不可能通过所有采样点； 2）忽略了海拔对采样点的影响。
对于这两个模型，由于受资料的限制，无法讨论污染源扩散与时间，天气，温度的关系。 4.4.2 地质环境演变模型
研究城市地质环境的演变环境首先要对该城市在不同时期内的水陆分布有初步的坐标定位以便作出比较。从土壤分类上可分为饱和土壤和非饱和土壤。饱和土壤的意思是土壤中所有孔隙全部充满水。一般情况下，土壤都是出于水分不饱和状态的。只有在淹渍条件下，才会出于水分饱和状态。所以假设城市土壤绝大部分为非饱和土壤。将土壤视为非饱和多孔介质, 描述非饱和多孔介质中多种迁移场量耦合的热质迁移机制, 与污染物在土壤中对流- 弥散迁移规律相结合, 并且考虑土壤固体骨架对重金属污染物吸附的影响因素, 建立相关的数学模型.得出有益的结论.总体上，可以将土壤中重金属迁移方式分为以下几种：1）土壤中迁移 2）水流中迁移 3）大气中沉降及降雨。以上三种形式中重金属分别是夹杂在气液固三种形式的水载体上被土壤吸附。 第一步，根据性质建立方程： （1）连续性方程
非饱和土壤中, 存在液相的蒸发或水蒸汽的冷凝, 因而在连续性方程中出现蒸发率或冷凝率m；混合气体中包含空气和水蒸汽; 水蒸汽的迁移是在整体气体迁移了蒸汽的扩散运动
?V?Vg?上，叠加
V:，从而呈现
VVg???V
????液相：
??ll??Ng???V?ll?l??m ； （1）
?????混合气体：
??g??N:???V?gg?g:??m ； （2）
?????水蒸气：
??Vg??N??????????V?g????m ； （3） VVV???g?（2）动量方程
在非饱和土壤中, 液相流动受到毛细抽吸力、Darcy 阻力和重力的作用; 将孔隙中的空气和蒸汽看作为理想混合气, 不仅考虑混合气的整体迁移, 而且也考虑蒸汽相对于空气
全国数学建模大赛A题获奖论文
