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

承诺书
我们仔细阅读了中国大学生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式〔包括、电子邮件、网上咨询等〕与队外的任何人〔包括指导教师〕研究、讨论与赛题有关的问题。

我们知道，抄袭别人的成果是违反竞赛规则的, 如果引用别人的成果或其他公开的资料〔包括网上查到的资料〕，必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。

我们郑重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。

如有违反竞赛规则的行为，我们将受到严肃处理。

我们参赛选择的题号是〔从A/B/C/D中选择一项填写〕： A
我们的参赛报名号为〔如果赛区设置报名号的话〕：
所属学校〔请填写完整的全名〕：同济大学
参赛队员(打印并签名) ：1. 冯建设
2. 赵云波
3. 刘雄飞
指导教师或指导教师组负责人(打印并签名)：
日期： 2011 年 9 月 11日赛区评阅编号〔由赛区组委会评阅前进行编号〕：
编号专用页
赛区评阅编号〔由赛区组委会评阅前进行编号〕：
全国统一编号〔由赛区组委会送交全国前编号〕：全国评阅编号〔由全国组委会评阅前进行编号〕
城市表层土壤重金属污染分析
摘要：
本文中，城市表层土壤金属污染分析需要综合不同区域，不同金属的综合影响，根据随机的数据采样点，通过统计与插值的分析方法进行处理。

首先可以考虑对采样点进行网格化的数据处理，然后通过Kriging 方法进行空间散乱点的插值处理。

通过对函数的插值结果的观察，与原始采样数据的空间分布相比较，可以发现二者具有较好的吻合度，结果令人满意。

对城区内不同区域重金属污染程度的综合评价，要对各重金属评价指标分别加权。

利用熵权法来确定各重金属评价指标的权重系数，熵权法的实际意义在这里表达得尤为明显，根据熵权法得到的相关系数均为正值，这一点也验证了熵权法在寻找个金属污染物权重时的正确性，然后由综合权重进行线性加和，得到各个区域的综合评定指标，同时根据金属含量背景值进行等级标准的划分，从而确定不同区域的污染程度，结果与实际完全符合，说明熵权法的运用是正确的，从
在建立重金属污染的传播特征模型，先假设了污染源的位置，然后考虑根据扩散定律建立模型，根据一维扩散方程建立模型，但是这样的话在数据处理上必须首先根据采样点浓度特征大致确定污染源的位置，然后建立方程，根据采样数据提取信息求解，由于采样点较多，本问题的处理可能会遗漏部分有用信息。

如果与实际偏差较大，没有得到有效结果，可以进一步寻求其他解决方法。

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

假设引入时间变量，将扩散模型修改为：
2222(,)(,)()(,,)u u u F x y k x y f x y u x y t ∂∂∂++-=∂∂∂
如此一来，还需要知道浓度的时间变化率，即不同时刻采样点的浓度数值。

然后要考虑二维(或三维) 变差函数。

但它们都是以一维函数为基础的，只不过要考虑各向同性或各向异性，以及结构的套合。

或者可以进一步的加强考虑海拔Z 的影响，得到四元的方程模型，然后建立求解。

目录
一问题重述 (3)
二、问题分析 (4)
三模型假设 (5)
四变量与符号说明 (5)
五模型建立与求解 (6)
1 随机场和区域化变量 (6)
2 变差函数的构造 (6)
3 平稳性假设和本征假设 (6)
4 实验变差函数 (7)
5 变差函数的理论模型 (7)
6 Kriging 方法插值 (8)
7 确定各重金属评价指标的权重系数 (12)
8 传播特征 (16)
六、模型评价与改良 (17)
1 模型评价 (17)
2 改良方向 (18)
参考文献 (18)
附件 (19)
附件一重金属空间分布模型的建立与求解的matlab实现 (19)
附件二利用熵权法进行各区域污染程度的综合评定的matlab实现 (22)
一问题重述
随着城市经济的快速发展和城市人口的不断增加，人类活动对城市环境质量的影响日显突出。

对城市土壤地质环境异常的查证，以及如何应用查证获得的海量数据资料开展城市环境质量评价，研究人类活动影响下城市地质环境的演变模式，日益成为人们关注的焦点。

按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为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个评价值的比重
e第i个因素的熵值
i
W 各金属权重矩阵
d j第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 (x i ，x j) = y (x i- x j)。

这里，y i
，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 项评价指标，形成原始指标数据矩阵
()i j m n X x ⨯=〔这里m=319，n=8〕，对于某项指标j x ，指标值ij X 的值越大，则该指标在综合
评价中所起的作用。

在信息论中信息熵
1ln m i ij ij
j e k P p ==-∑ 表示系统的有序程度，一个系统的有序程度越高，则信息熵越大，反之，一个系
统的无序程度越高，则信息熵越小。

所以，可以根据各项指标的指标值的差异程度，利用信息熵的这个工具计算出各指标的权重。

1. 数据处理
111212122
212
,,,,,,,,,n n i j m m mn X X X X X X X X X X ⋯⎡⎤⎢⎥
⋯⎢⎥=⎢⎥⎢⎥⋯⎣⎦ (i=1,2,…,m;j= 1,2,…,n,)，这里取n=8，m=319。

由于参与评价的各项指标有越大越优型、越小越优型，故需对矩阵中的特征值进行归一化处理，方法如下：
据此得到归一化矩阵X ’:
111212122
212
',',,'',',,''',',,'n n i j m m mn X X X X X X X X X X ⋯⎡⎤⎢⎥
⋯⎢⎥=⎢⎥⎢⎥⋯⎣⎦ (i=1,2,…,m;j= 1,2,…,n,)
这里以越大越优型进行计算求解。

i 个因素下第j 个评价值的比重ij P
1ij
ij m ij
j x P x
==∑ 3. 计算第i 个因素的熵值i e
1ln ,
1,01ln m
i ij ij j i e k P p e m
==-≤≤∑若取k=则
i 个因素的差异系数i g
对于给定的i e 越大，因素评价值的差异性越小，则因素在综合评价中所起的作用越小。

定义差异系数1i i g e =-，则当因素i g 越大时，因素越重要。

1i ij m j
j g w g
==∑，则j w 就是熵权法确定的权重。

记As 、Cd 、Cr 、Cu 、Hg 、Ni 、Pb 、Zn 的权重矩阵为
W=[w 1, w 2, w 3, w 4, w 5, w 6, w 7, w 8]’
计算得到：
0.49028 0.022241 0.044582 0.123778]’
得到各取样点加权综合浓度指标
D=X W=[d 1,d 2,…,d 319]’；
进一步得到i 类区加权综合浓度指标
11i n i j j i d d n ==∑，
其中，d j 为所有属于i 类区的d ，n i 为其采样点个数。

6.得到评价结果
得到五个区的加权综合浓度指标如下：
区
1 2 3 4 5 加权综合浓度指标值
对各重金属浓度背景值及标准差做相同归一化处理，得到归一化之后的背景平均值及标准差，然后对8种金属归一化之后的背景平均值和标准差进行加权求和，求得和各点加权综合浓度指标相对应的综合背景平均值和标准差如下：
平均值对应指标 标准差对应指标
设
轻度污染标准=平均值对应指标+3*标准差对应指标；
中度污染标准=平均值对应指标+8*标准差对应指标；
重度污染标准=平均值对应指标+15*标准差对应指标；
污染程度
无 轻度污染 中度污染 重度污染 得分
故可得到各区的污染等级如下表：
区 1
生活区2
工业区
3
山区
4
主干道
5
公园绿地
污染等级轻度污染重度污染无污染重度污染轻度污染
由题设条件可知，按照功能划分，城区一般可分为生活区、工业区、山区、主干道路区及公园绿地区等，分别记为1类区、2类区、……、5类区，不同的区域环境受人类活动影响的程度不同，通过以上的计算，得出的结论为：工业区与主干道严重污染，生活区与公园绿地轻度污染，山区几乎不受污染，这与实际情况符合的较好。

根据该金属元素空间分布模型，可以确定工业区污染较为严重，故可以初步确定重金属污染的主演原因在于工业区活动对环境的影响。

得到污染源分布如下：
8 传播特征
第一问中已经确定了污染物的浓度场，记为u(x,y)，本例中，按稳定场处理。

由扩散定律，得到扩散方程为：
2222(,)(,)()(,,)0u u F x y k x y f x y u x y ∂∂++-=∂∂
其中，F(x,y)为污染源强度，
K(x,y)为〔x ，y 〕处土壤传播污染物的特性，
f(x,y,u)为耗散强度。

建立离散模型：
设污染源为点源F(x i ,y j )只有在源(x i ,y j )处等于污染物的排放强度，在源处，恒为零。

假设各分区内的土壤特性相同。

则在非源处，
2222(,)1()i j i j u u k x y x y ∂∂=+∂∂
则各区的
1j i k k n =∑
其中，k j 为j 类区的平均土壤特性，
K i 为i 采样点处的土壤特性。

因为城区内土壤特性相同，故f(x,y,u)可以简写为f(u)。

取浓度场上的节点，建立离散方程组，通过求解方程组，可以分别得到F(x,y)，f(u)的一组数值，进行拟合，得出解析式，将解析式带回扩散方程，即可得出扩散规律。

六、模型评价与改良
1 模型评价
1．为了给出8种主要重金属的空间分布，并分析城区不同区域的污染程度。

首先对空间位置比较散乱的采样点此进行网格化的数据处理，然后通过Kriging方法进行空间散乱点的插值处理，在空间数据场满足给定的统计分布特征的前提假设下进行插值。

在求得相关参数后建立显示的高斯模型，通过对函数的插值结果的观察，与原始采样数据的空间分布相比较，可以发现二者具有较好的吻合度，结果令人满意。

2.对城区内不同区域重金属污染程度的综合评价，要对各重金属评价指标分别加权。

利用熵权法来确定各重金属评价指标的权重系数，熵权法的实际意义在这里表达得尤为明显，根据熵权法得到的相关系数均为正值，这一点也验证了熵权法在寻找个金属污染物权重时的正确性，然后由综合权重进行线性加和，得到各个区域的综合评定指标，同时根据金属含量背景值进行等级标准的划分，从而确定不同区域的污染程度，结果与实际完全符合，说明熵权法的运用是正确的，从而便于找到重金属污染的主要原因。

2 改良方向
1. 在建立重金属污染的传播特征模型，先假设了污染源的位置，然后考虑根据扩散定律建立模型，根据一维扩散方程建立模型，但是这样的话在数据处理上必须首先根据采样点浓度特征大致确定污染源的位置，然后建立方程，根据采样数据提取信息求解，由于采样点较多，本问题的处理可能会遗漏部分有用信息。

如果与实际偏差较大，没有得到有效结果，可以进一步寻求其他解决方法。

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

假设引入时间变量，将扩散模型修改为：
2222(,)(,)()(,,)u u u F x y k x y f x y u x y t ∂∂∂++-=∂∂∂
如此一来，还需要知道浓度的时间变化率，即不同时刻采样点的浓度数值。

然后要考虑二维(或三维) 变差函数。

但它们都是以一维函数为基础的，只不过要考虑各向同性或各向异性， 以及结构的套合。

或者可以进一步的加强考虑海拔Z 的影响，得到四元的方程模型，然后建立求解。

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