2023年高教社杯全国大学生数学建模竞赛.doc

高教社杯全国大学生数学建模竞赛

承 诺 书

咱们仔细阅读了中华人民共和国大学生数学建模竞赛竞赛规则.
咱们完全明白，在竞赛开始后参赛队员不能以任何方式（波及 、电子邮件、网上征询等）与队外任何人（波及指导教师）研究、讨论与赛题有关问题。
咱们懂得，抄袭他人成果是违反竞赛规则，假如引用他人成果或其她公开资料（波及网上查到资料），必要按照规定参照文献表述方式在正文引用处和参照文献中明确列出。
咱们郑重承诺，严格遵守竞赛规则，以保证竞赛公正、公平性。如有违反竞赛规则行为，咱们将受到严厉处理。

咱们参赛选用题号是（从A/B/C/D中选用一项填写）： A
咱们参赛报名号为（假如赛区设置报名号话）：
所属学校（请填写完整全名）： 同济大学
参赛队员 (打印并签名) ：1. 冯建设
2. 赵云波
3. 刘雄飞
指导教师或指导教师组负责人 (打印并签名)：

日期： 年 9 月 11日

赛区评阅编号（由赛区组委会评阅前进行编号）：

高教社杯全国大学生数学建模竞赛

编 号 专 用 页

赛区评阅编号（由赛区组委会评阅前进行编号）：

赛区评阅记录（可供赛区评阅时使用）：

评
阅
人

评
分

备
注

全国统一编号（由赛区组委会送交全国前编号）：

全国评阅编号（由全国组委会评阅前进行编号）

都市表层土壤重金属污染分析

摘要：
本文中，都市表层土壤金属污染分析需要综合不一样区域，不一样金属综合影响，根据随机数据采样点，通过记录与插值分析措施进行处理。
首先可以考虑对采样点进行网格化数据处理，然后通过Kriging措施进行空间散乱点插值处理。通过对函数插值成果观测，与原始采样数据空间分布相比较，可以发现两者具有很好吻合度，成果令人满意。
对城区内不一样区域重金属污染程度综合评价，要对各重金属评价指标分别加权。运用熵权法来确定各重金属评价指标权重系数，熵权法实际意义在这里体现得尤为明显，根据熵权法得到有关系数均为正值，这一点也验证了熵权法在寻找个金属污染物权重时对旳性，然后由综合权重进行线性加和，得到各个区域综合评估指标，同步根据金属含量背景值进行级别原则划分，从而确定不一样区域污染程度，成果与实际完全符合，阐明熵权法运用是对旳，从而找到重金属污染重要原因。
区号
区类型
1
生活区
2
工业区
3
山区
4
主干道
5
公园绿地
污染级别
轻度污染
重度污染
无污染
重度污染
轻度污染

在建立重金属污染传播特性模型，先假设了污染源位置，然后考虑根据扩散定律建立模型，根据一维扩散方程建立模型，不过这样话在数据处理上必要首先根据采样点浓度特性大体确定污染源位置，然后建立方程，根据采样数据提取信息求解，由于采样点较多，本问题处理也许会遗漏某些有用信息。假如与实际偏差较大，没有得到有效成果，可以深入寻求其她处理措施。
对于问题四，可以基于问题一二三得到模型进行综合考虑，从数据处理方式，以及地质演变过程中时间变量影响等原因入手，寻求更好建立模型方式。若引入时间变量，将扩散模型修改为：

如此一来，还需要懂得浓度时间变化率，即不一样步刻采样点浓度数值。然后要考虑二维(或三维) 变差函数。 但它们都是以一维函数为基本，只不过要考虑各向同性或各向异性，以及构造套合。或者可以深入加强考虑海拔Z影响，得到四元方程模型，然后建立求解。

目 录

一 问题重述 3
二、问题分析 4
三 模型假设 5
四 变量与符号阐明 5
五 模型建立与求解 6
1 随机场和区域化变量 6
2 变差函数构造 6
3 平稳性假设和本征假设 6
4 试验变差函数 7
5 变差函数理论模型 7
6 Kriging 措施插值 8
7 确定各重金属评价指标权重系数 13
8 传播特性 16
六、模型评价与改善 17
1 模型评价 17
2 改善方向 17
参照文献 18
附件 19
附件一 重金属空间分布模型建立与求解matlab实现 19
附件二 运用熵权法进行各区域污染程度综合评估matlab实现 21

一 问题重述

伴随都市经济迅速发展和都市人口不停增长，人类活动对都市环境质量影响日显突出。对都市土壤地质环境异常查证，以及怎样应用查证获得海量数据资料开展都市环境质量评价，研究人类活动影响下都市地质环境演变模式，日益成为人们关注焦点。
按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不一样区域环境受人类活动影响程度不一样。
现对某都市城区土壤地质环境进行调查。为此，将所考察城区划分为间距1公里左右网格子区域，按照每平方公里1个采样点对表层土（0~10 厘米深度）进行取样、编号，并用GPS记录采样点位置。应用专门仪器测试分析，获得了每个样本所含多种化学元素浓度数据。另首先，按照2公里间距在那些远离人群及工业活动自然区取样，将其作为该城区表层土壤中元素背景值。
附件1列出了采样点位置、海拔高度及其所属功能区等信息，附件2列出了8种重要重金属元素在采样点处浓度，附件3列出了8种重要重金属元素背景值。
现规定通过数学建模来完毕如下任务：
(1) 给出8种重要重金属元素在该城区空间分布，并分析该城区内不一样区域重金属污染程度。
(2) 通过数据分析，阐明重金属污染重要原因。
(3) 分析重金属污染物传播特性，由此建立模型，确定污染源位置。
(4) 分析你所建立模型优缺陷，为更好地研究都市地质环境演变模式，还应搜集什么信息？有了这些信息，怎样建立模型处理问题？

二、问题分析

本题中，都市表层土壤金属污染分析需要综合不一样区域，不一样金属综合影响，根据随机数据采样点，考虑通过记录与插值分析措施进行处理。

对于问题一，需要给出8种重要重金属空间分布，并分析城区不一样区域污染程度。由于采样数据给出采样点相对来说空间位置比较散乱，首先可以考虑对此进行网格化数据处理，然后通过Kriging措施进行空间散乱点插值处理。
Kriging是一种距离加权插值措施，从地质记录学中借鉴而来，在空间数据场满足给定记录分布特性前提假设下进行插值。 通过对函数插值成果观测，与原始采样数据空间分布相比较， 深入确定该措施可以获得插值效果与否有效。
对城区内不一样区域重金属污染程度综合评价，要对各重金属评价指标分别加权。权重是衡量因子集中某一因子对土壤污染程度影响相对大小量,权重系数越大,则该因子对土壤影响程度越大，考虑运用熵权法来确定各重金属评价指标权重系数，在信息论中信息熵体现系统有序程度，一种系统有序程度越高，则信息熵越大，反之，一种系统无序程度越高，则信息熵越小。因此，可以根据各项指标指标值差异程度，运用信息熵这个工具计算出各指标权重。然后由综合权重进行线性加和，可以得到各个区域综合评估指标，同步根据金属含量背景值进行级别原则划分，从而确定不一样区域污染程度。

对于问题二，可以结合问题一中得到模型，同步对采样数据进行简朴分析，根据各区域污染程度不一样，可以感性得到重金属污染重要原因必然和污染最严重区域有直接关系。

对于问题三，需要建立重金属污染传播特性模型，并确定污染源位置。可以考虑根据扩散定律建立模型，首先借鉴一维扩散方程建立模型，不过这样话在数据处理上必要首先根据采样点浓度特性大体确定污染源位置，然后建立方程，根据采样数据提取信息求解，由于采样点较多，本问题处理也许会遗漏某些有用信息。假如与实际偏差较大，可以深入寻求其她处理措施。

对于问题四，可以基于问题一二三得到模型进行综合考虑，从数据处理方式，以及地质演变过程中时间变量影响等原因入手，寻求更好建立模型方式。

三 模型假设

1． 金属污染浓度场按稳定场处理，即各坐标点浓度不随时间变化。
2． 假设各区域土壤特性相似。
3． 重力对金属污染影响忽视。

四 变量与符号阐明

变量符号
符号阐明
X i
取样点横坐标
y i
取样点纵坐标

点x处半差函数

各点对应权重系数

估计值点与i点之间变差函数值

i点与j点之间变差函数值

某金属浓度值矩阵

某金属浓度值归一化矩阵

第 i个原因下第 j个评价值比重

第个原因熵值
W
各金属权重矩阵
dj
第j区加权综合浓度指标

五 模型建立与求解

模型需要给出8种重要重金属空间分布，并分析城区不一样区域污染程度。由于采样数据给出采样点相对来说空间位置比较散乱，首先可以考虑对此进行网格化数据处理，然后通过Kriging措施进行空间散乱点插值处理。
Kriging是一种距离加权插值措施，从地质记录学中借鉴而来，在空间数据场满足给定记录分布特性前提假设下进行插值。
运用Kriging措施对各浓度散乱点进行插值处理，建立起浓度分布空间曲面（曲面高度值即用来表征重金属浓度值），详细简介如下：Kriging 措施就是对空间数据进行加权插值权值设计措施。Kriging措施通过引进以距离为自变量变差函数来计算权值。由于变差函数既可以反应变量空间构造特性，又可以反应变量随机分布特性，因此运用Kriging措施进行空间数据插值往往可以获得理想效果。此外，通过设计变差函数，Kriging措施很轻易实现局部加权插值，这样就克服了一般距离加权插值措施插值成果不稳定性。

1 随机场和区域化变量
首先，重金属浓度数据场可体现为分布于空间单值函数S = f (x ，y ，z )，由题意知S为标量，则数据场为标量场。运用记录学措施来研究该数据场，首先将f 当作随机函数，记为Z，依赖于多种自变量随机函数，称为随机场。以空间点x 直角坐标为自变量随机场称为一种区域化变量。区域化变量在观测前，可以看作是随机场；观测后就得到随机场一种实现(一般都记作Z (x ) ，写法上不加区别)。浓度区域化变量同步反应地质变量构造性和随机性特性。从地质学观点来看，区域化变量可反应地质变量如下特性：局部性、持续性、异向性、可迁性。

2 变差函数构造
假设空间点x 只在一维x 轴上变化， 咱们把区域化变量Z (x ) 在x 和x + h (h 为与x 具有相似维数距离向量) 2个点处值之差方差之半定义为Z (x ) 在x 方向上变差函数，记为y (x，h)，即

深入，由于点x 和h 是在二维(或三维) 空间中变化， 因此要考虑二维(或三维) 变差函数。但它们都是以一维变差函数为基本，只不过要考虑各向同性或各向异性，还要考虑构造套合。这里临时先以各向同性进行数据分析处理。

3 平稳性假设和本征假设
当区域化变量Z (x ) 满足下列条件时，则称Z (x ) 满足二阶平稳(或弱平稳)。
(1) 在整个研究区域内， 区域化变量Z (x ) 数学期望存在， 且等于常数， 即E [Z (x ) ]= m (常数)，P x；
(2) 在整个研究区域内，区域化变量Z (x ) 协方差函数存在且相似(即只依赖于滞后h，而与x 无关)，即

在实际工作中常常连二阶平稳假设也不能满足，故提出本征假设。
当区域化变量Z (x ) 增量[Z (x ) - Z (x + h) ]满足下列条件时，则称Z (x ) 满足本征假设，或说Z (x ) 是本征。
(1) 在整个研究区域内有，区域化变量Z (x ) 数学期望存在，且等于常数， 即E [Z (x ) - Z (x + h) ]= 0，P x，P h；
(2) 在整个研究区域内，增量[Z (x ) - Z (x + h) ]方差函数存在且平稳(不依赖于x ) ，即

4 试验变差函数
试验变差函数就是根据观测数据构造变差函数y (h) 估计值。有了二阶平稳假设或本征假设，重金属浓度区域化变量Z (x ) 增量[Z (x ) - Z (x + h) ]只依赖于分隔它们h，而不依赖于详细位置x。这样，被向量h 分割每一种数据对{Z (x i),Z (x i+ h) } （ i= 1， 2， ⋯， N (h ) ） 可以当作是{Z (x ),Z (x + h) }一次不一样实现（N (h ) 是被向量h 相隔数据对个数）。 这样便可以用求所有这些观测值算术平均措施来计算。于是得到

这就是试验变差函数基本计算公式。

5 变差函数理论模型
为了对区域化变量未知值作出估计，需要将试验变差函数拟合成对应理论变差函模型，这些模型将直接参与Kriging措施计算。变差函数模型可以分为有基台值和无基台值两大类，这里运用常用有基台值模型进行计算。3种常用有基台值模型如下：
(1) 球状模型(亦称马特隆模型，在原点处为线性型)。
球状模型几何解释，是由于它来源于2个半径为a 且球心距
为2h 球体重叠某些体积计算公式。它一般体现为

其中，为块金常数，+ C 为基台值，C 为拱高，a 为变程。
当= 0，C= 1时，称为原则球状模型。其中，为块金常数，+ C 为基台值， C 为拱高，a 为变程。
(2) 指数函数模型(在原点处为线性型)。它一般公式为
此处a 不是变程。由于当h= 3a 时，有1- e- 3≈ 0.95≈ 1。因此y (h)≈ + C， 故其变程为3a。当= 0，C = 1时称为原则指数函数模型。
(3) 高斯模型(在原点处为抛物线型)。它一般公式为

此处a 亦不是变程。 由于当h= a 时，有1- e- 3≈ 0.95≈ 1，因此y (h)≈ + C，故其变程为3a。当= 0，C= 1时称为原则指数函数模型。

6 Kriging 措施插值
设Z (x ) 是点x承载区域化变量，且是二阶平稳(或本征) 。Z i (i= 1， 2， ⋯， n)是一组离散信息样品数据，即重金属浓度，它们是定义在点承载x i ( i= 1， 2， ⋯， n) 上。现要对点x 0承载处区域化变量进行估计，所用估计量为 ，它是n 个数值加权线性组合。 Kriging 措施原则， 就是在保证这个估计量是无偏，且估计方差最小前提下，求出n 个权值系数。
(1) 无偏性条件
若要使 为无偏估计量，即规定E [ - ]= 0，由此得到无偏条件：
(2) Kriging 方程组
区域化变量在满足二阶平稳条件下推导，可以得到估计方差计算公式

估计方差对 偏导数为
。
在无偏性条件下，为了使估计方差最小，这是个求条件极值问题， 要用到拉格朗日乘子法。
令。
这里F 是n 个权系数 和(n+ 1) 元函数，求出F对 ( i= 1， 2， ⋯， n) 和偏导数， 并令其为零， 便得到下列Kriging 方程组

整顿得

假如区域化变量只满足本征假设，而不满足二阶平稳假设，则运用协方差函数和变差函数关系C (h) = C (0) - y (h) ， 可得用变差函数体现Kriging方程组

这里， y i， j= y (x i ， x j ) = y (x i- x j )。
详细算法实现过程通过数学软件matlab编程实现，（matlab代码见附件）。

(3) kriging模型详细实现
随机取100个观测点作为原始数据，其分布如下图所示：

值得指出是，在初步数据处理中，可以发现该城区地形并不是规则长方形，不过为了以便数据处理，这里在一种30*20区域内随机取点，刚好可以完全涵盖该城区，不过这并不影响最终成果，由于可以再得到模型之后再将不属于城区区域挖去，这在matlab软件中是轻易实现。
根据试验变差函数成果，作数据拟合，选用恰当变差函数模型进行插值， 确定变差函数计算公式。这里有3种常用变差函数模型：球状模型、指数模型和高斯模型。不过这3种模型有3个关键参数需要确定，变程a、拱高C 和块金常数C0。本题中块金常数C0=0，拱高C 和变程a 由试验变差函数图人工确定。
对于第一种重金属As，由试验变差函数100个观测点观测值得到变差函数，可以计算得到如图：

可以求得近似值：C= 20和a= 3。然后确定变差函数模型，选定变差函数模型后，变差函数计算公式可以显式写出。这样Kriging 方程组系数矩阵与增广矩阵都已经确定。然后进行Kriging 插值，也就是求解Kriging 方程组，确定加权系数，然后进行线性加权，即可得到该模型下Kriging 插值成果。（详细实现见附件matlab代码）
经画图比较，在本题中高斯变差函数模型符合程度很好，建立高斯变差函数模型之后得到如下图所示函数图像：

1号金属As分布原始数据成果 1号金属As分布kriging插值成果

同样方式，对第二种重金属Cd浓度做变差处理，值得指出是，为以便比较，这里不再进行随机取100个观测点工作，而是直接运用在计算第一种金属As变差函数时候随机点数据，在此基本上得到变差函数观测数据，如下图：

由试验变差函数观测图可以近似求得 C= 126000和a= 4。然后确定变差函数模型，选定变差函数模型后， 变差函数计算公式可以显式写出。这样Kriging 方程组系数矩阵与增广矩阵都已经确定。然后进行Kriging 插值，也就是求解Kriging 方程组，确定加权系数，然后进行线性加权，即可得到该模型下Kriging 插值成果。经画图比较，在本题中高斯变差函数模型符合程度很好，建立高斯变差函数模型之后得到如下图所示图像：

2号金属分布原始数据 2号金属kriging插值数据

对应，运用kriging措施也可以深入得到其六种元素空间分布模型，作图如下：

7 确定各重金属评价指标权重系数

在综合评估该城区内不一样区域重金属污染程度时，还要对各重金属评价指标分别加权。权重是衡量因子集中某一因子对土壤污染程度影响相对大小量,权重系数越大,则该因子对土壤影响程度越大，在这里运用熵权法来确定各重金属评价指标权重系数，详细原理及操作如下：
措施长处：客观赋权法。
背景：设有m个待评方案，n项评价指标，形成原始指标数据矩阵（这里m=319，n=8），对于某项指标，指标值值越大，则该指标在综合评价中所起作用。
在信息论中信息熵

体现系统有序程度，一种系统有序程度越高，则信息熵越大，反之，一种系统无序程度越高，则信息熵越小。因此，可以根据各项指标指标值差异程度，运用信息熵这个工具计算出各指标权重。
1. 数据处理

(i=1,2,…,m;j= 1,2,…,n,)，这里取n=8，m=319。
由于参与评价各项指标有越大越优型、越小越优型，故需对矩阵中特性值进行归一化处理，措施如下：

据此得到归一化矩阵X’:

(i=1,2,…,m;j= 1,2,…,n,)
这里以越大越优型进行计算求解。
2.计算第个原因下第个评价值比重

3. 计算第个原因熵值

4.计算第个原因差异系数
对于给定越大，原因评价值差异性越小，则原因在综合评价中所起作用越小。定义差异系数，则当原因越大时，原因越重要。
5.定义权数，则就是熵权法确定权重。
记As、Cd、Cr、Cu、Hg、Ni、Pb、Zn权重矩阵为
W=[w1，w2，w3，w4，w5，w6，w7，w8]’
计算得到：
W=[0.023335 0.047536 0.065113 0.183134
0.49028 0.022241 0.044582 0.123778]’

得到各取样点加权综合浓度指标
D=X W=[d1,d2,…,d319]’；
深入得到i类区加权综合浓度指标
，
其中，dj为所有属于i类区d，ni为其采样点个数。

6.得到评价成果

得到五个区加权综合浓度指标如下：
区
1
2
3
4
5
加权综合浓度指标值
0.002167
0.005942
0.000953
0.004096
0.001787

对各重金属浓度背景值及原则差做相似归一化处理，得到归一化之后背景平均值及原则差，然后对8种金属归一化之后背景平均值和原则差进行加权求和，求得和各点加权综合浓度指标相对应综合背景平均值和原则差如下：

平均值对应指标
原则差对应指标
0.000799
0.000195
在此基本上得到对污染程度级别进行分类值如下：
设
轻度污染原则=平均值对应指标+3*原则差对应指标；
中度污染原则=平均值对应指标+8*原则差对应指标；
重度污染原则=平均值对应指标+15*原则差对应指标；
即：
污染程度
无
轻度污染
中度污染
重度污染
得分
<0.001385
0.001385~0.002361
0.002361~0.003727
>0.003727

故可得到各区污染级别如下表：
区
1
生活区
2
工业区
3
山区
4
主干道
5
公园绿地
污染级别
轻度污染
重度污染
无污染
重度污染
轻度污染

由题设条件可知，按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不一样区域环境受人类活动影响程度不一样，通过以上计算，得出结论为：工业区与主干道严重污染，生活区与公园绿地轻度污染，山区几乎不受污染，这与实际状况符合很好。
根据该金属元素空间分布模型，可以确定工业区污染较为严重，故可以初步确定重金属污染主演原因在于工业区活动对环境影响。

得到污染源分布如下：

8 传播特性

第一问中已经确定了污染物浓度场，记为u(x,y)，本例中，按稳定场处理。
由扩散定律，得到扩散方程为：

其中，F(x,y)为污染源强度，
K(x,y)为（x，y）处土壤传播污染物特性，
f(x,y,u)为耗散强度。
建立离散模型：
设污染源为点源F(xi,yj)只有在源(xi,yj)处等于污染物排放强度，在源处，恒为零。
假设各分区内土壤特性相似。 则在非源处，

则各区

其中，kj为j类区平均土壤特性，
Ki为i采样点处土壤特性。
由于城区内土壤特性相似，故f(x,y,u)可以简写为f(u)。
取浓度场上节点，建立离散方程组，通过求解方程组，可以分别得到F(x,y)，f(u)一组数值，进行拟合，得出解析式，将解析式带回扩散方程，即可得出扩散规律。

六、模型评价与改善

1 模型评价

1．为了给出8种重要重金属空间分布，并分析城区不一样区域污染程度。首先对空间位置比较散乱采样点此进行网格化数据处理，然后通过Kriging措施进行空间散乱点插值处理，在空间数据场满足给定记录分布特性前提假设下进行插值。 在求得有关参数后建立显示高斯模型，通过对函数插值成果观测，与原始采样数据空间分布相比较， 可以发现两者具有很好吻合度，成果令人满意。
2.对城区内不一样区域重金属污染程度综合评价，要对各重金属评价指标分别加权。运用熵权法来确定各重金属评价指标权重系数，熵权法实际意义在这里体现得尤为明显，根据熵权法得到有关系数均为正值，这一点也验证了熵权法在寻找个金属污染物权重时对旳性，然后由综合权重进行线性加和，得到各个区域综合评估指标，同步根据金属含量背景值进行级别原则划分，从而确定不一样区域污染程度，成果与实际完全符合，阐明熵权法运用是对旳，从而便于找到重金属污染重要原因。

2 改善方向
1. 在建立重金属污染传播特性模型，先假设了污染源位置，然后考虑根据扩散定律建立模型，根据一维扩散方程建立模型，不过这样话在数据处理上必要首先根据采样点浓度特性大体确定污染源位置，然后建立方程，根据采样数据提取信息求解，由于采样点较多，本问题处理也许会遗漏某些有用信息。假如与实际偏差较大，没有得到有效成果，可以深入寻求其她处理措施。

2.对于问题四，可以基于问题一二三得到模型进行综合考虑，从数据处理方式，以及地质演变过程中时间变量影响等原因入手，寻求更好建立模型方式。若引入时间变量，将扩散模型修改为：

如此一来，还需要懂得浓度时间变化率，即不一样步刻采样点浓度数值。然后要考虑二维(或三维) 变差函数。 但它们都是以一维函数为基本，只不过要考虑各向同性或各向异性， 以及构造套合。或者可以深入加强考虑海拔Z影响，得到四元方程模型，然后建立求解。

参照文献

[1]
周晓云,朱心雄. 散乱数据点三角剖分措施综述[J]. 工程图学学报，1993,(01) .
[2]
王靖波,潘懋,张绪定. 基于Kriging措施空间散乱点插值[J]. 计算机辅助设计与图形学学报，1999,(06) .
[3] 乔家君. 改善熵值法在河南省可持续发展能力评估中应用[J]. 资源科学，,(01) .

附件

附件一 重金属空间分布模型建立与求解matlab实现
注：为减少篇幅，此处代码仅仅实现了对第一种金属元素含量空间分布求解，别旳个元素求解方式类似可以得到。

%坐标网格化
%zuobiao为采样点坐标矩阵
%xx将坐标单位化为公里
x=zuobiao(:,1)';
y=zuobiao(:,2)';
z=zuobiao(:,3)';
xx=x/1000;
yy=y/1000;
[X,Y]=meshgrid(0:0.5:29,0:0.5:19);
Z=griddata(xx,yy,z,X,Y,'v4');
surf(X,Y,Z);
shading interp

%金属浓度数据网格化
jinshu_1=jinshu(:,1)';
J_1=griddata(xx,yy,jinshu_1,X,Y,'v4');
jinshu_2=jinshu(:,2)';
J_2=griddata(xx,yy,jinshu_2,X,Y,'v4');
jinshu_3=jinshu(:,3)';
J_3=griddata(xx,yy,jinshu_3,X,Y,'v4');
jinshu_4=jinshu(:,4)';
J_4=griddata(xx,yy,jinshu_4,X,Y,'v4');
jinshu_5=jinshu(:,5)';
J_5=griddata(xx,yy,jinshu_5,X,Y,'v4');
jinshu_6=jinshu(:,6)';
J_6=griddata(xx,yy,jinshu_6,X,Y,'v4');
jinshu_7=jinshu(:,7)';
J_7=griddata(xx,yy,jinshu_7,X,Y,'v4');
jinshu_8=jinshu(:,8)';
J_8=griddata(xx,yy,jinshu_8,X,Y,'v4');

%注：产生100个随机坐标代码只在对第一号金属进行变差函数观测时候运行，在对其
%她金属进行变差函数观测求解时候只需运用目前100个观测点即可，该代码不再运行
%产生100个随机坐标
%s_suiji即为随机观测点得坐标
%for i=1:100
% x_suiji(i,1:2)=[ceil(rand*30),ceil(rand*20)];
%end

%求变差函数原始数据
%bianchahanshu为变差函数矩阵
bianchahanshu=[0,0];
for j=1:100
for k=1:100
bianchahanshu=[bianchahanshu;...
[sqrt((x_suiji(j,1)-x_suiji(k,1)).^2+(x_suiji(j,2)-x_suiji(k,2)).^2),...
(J_1(x_suiji(j,2),x_suiji(j,1))-J_1(x_suiji(k,2),x_suiji(k,1))).^2]];
end
end

%对变差函数矩阵进行处理，得到最终变差函数矩阵
for m=1:600
flag=find(bianchahanshu(:,1)==bianchahanshu(m,1));
bianchahanshu(m,:)=sum(bianchahanshu(flag(:,1),:))/sum(flag==flag);
bianchahanshu(flag(2:end,1),:)=[];m
end
bianchahanshu(:,2)=bianchahanshu(:,2)/2;

%做出变差函数观测点图像
scatter(bianchahanshu(:,1),bianchahanshu(:,2))

%根据变差函数，求每一种估计点权系数系数矩阵以及增广矩阵
%C为系数矩阵
%C_plus为增广矩阵
C=zeros(100,100);
for p=1:100
for q=1:100
h=sqrt((x_suiji(p,1)-x_suiji(q,1)).^2+(x_suiji(p,2)-x_suiji(q,2)).^2);
C(p,q)=20*(1-exp(-h/3));
end
end

C_plus=zeros(101,101);
C_plus(1:100,1:100)=C;
C_plus(101,1:100)=1;
C_plus(1:100,101)=1;
C_plus(101,101)=0;

%求解每一种估计点权重系数矩阵，得到空间模型
%suijinongduJ_1为100个随机观测点第一重金属浓度观测值
%lamda为权重系数矩阵
% jisuannongduJ_1即为根据模型得到每一种需要估计点浓度值
for n=1:100
suijinongduJ_1(n,1)=J_1(x_suiji(n,2),x_suiji(n,1));
end
lamda=zeros(101,1);
for ii=1:30
for jj=1:20
for i=1:100
h0=sqrt((x_suiji(i,1)-ii).^2+(x_suiji(i,2)-jj).^2);
y0(i,1)=20*(1-exp(-h0/3));
end
lamda=C_plus\[y0;1];
jisuannongduJ_1(jj,ii)=suijinongduJ_1'*lamda(1:end-1,1);
end
end
surf(X,Y,jisuannongduJ_1);shading interp

附件二 运用熵权法进行各区域污染程度综合评估matlab实现

%运用熵权法求8种重金属评价指标权重系数
% jinshuguiyi_1为归一化之后数据
%g为得到权重系数矩阵
jinshuguiyi_1(:,:)=jinshu(:,:)/sum(jinshu);
H=-1/log(319).*sum(jinshuguiyi_1.*log(jinshuguiyi_1));
g=(1-H)/(8-sum(H));

%对各不一样区域数据分开处理
% jinshu_1_1_1为1区平均化金属浓度矩阵
% result_1为1区加权综合浓度指标，下面各区类似名具有类似意义
a1=find(jinshu(:,9)==1);
jinshu_1_1=jinshuguiyi_1(a1(:,1),:);
jinshu_1_1_1=sum(jinshu_1_1)/sum(a1==a1);
result_1=g*jinshu_1_1_1';

a2=find(jinshu(:,9)==2);
jinshu_1_2=jinshuguiyi_1(a2(:,1),:);
jinshu_1_2_1=sum(jinshu_1_2)/sum(a2==a2);
result_2=g*jinshu_1_2_1';

a3=find(jinshu(:,9)==3);
jinshu_1_3=jinshuguiyi_1(a3(:,1),:);
jinshu_1_3_1=sum(jinshu_1_3)/sum(a3==a3);
result_3=g*jinshu_1_3_1';

a4=find(jinshu(:,9)==4);
jinshu_1_4=jinshuguiyi_1(a4(:,1),:);
jinshu_1_4_1=sum(jinshu_1_4)/sum(a4==a4);
result_4=g*jinshu_1_4_1';

a5=find(jinshu(:,9)==5);
jinshu_1_5=jinshuguiyi_1(a5(:,1),:);
jinshu_1_5_1=sum(jinshu_1_5)/sum(a5==a5);
result_5=g*jinshu_1_5_1';