本程序主要是求二维实验变差函数值可以分别求doc.docx
《本程序主要是求二维实验变差函数值可以分别求doc.docx》由会员分享，可在线阅读，更多相关《本程序主要是求二维实验变差函数值可以分别求doc.docx（21页珍藏版）》请在冰豆网上搜索。
本程序主要是求二维实验变差函数值可以分别求doc
地理统计（Geo-statistics）简介
浙江大学农业与生物技术学院唐启义
（浙江大学华家池校区杭州310029）
当一个变量呈空间分布时，就称之为区域化变量。
这种变量反映了空间某种属性的分布
特征。
矿产、地质、海洋、土壤、气象、水文、生态、温度、浓度等领域都具有某种空间属
性。
区域化变量具有双重性，在观测前区域化变量Z（x）是一个随机场，观测后是一个确定的
空间点函数值。
区域化变量具有两个重要的特征。
一是区域化变量Z（x）是一个随机函数,它具有局部的、
随机的、异常的特征；其次是区域化变量具有一般的或平均的结构性质，即变量在点X与
偏离空间距离为h的点x＋h处的随机量Z（x）与Z（x+h）具有某种程度的自相关，而且这种自相关性依赖于两点间距离h与变量特征。
在某种意义上说这就是区域化变量的结构性特征。
一、实验半变异函数
变异函数又称变差函数、变异矩，是地统计分析所特有的基本工具。
函数定义为，当空间点x在一维x轴上变化时，区域化变量Z（x）在点x和
Z（x+h）差的方差的一半为区域化变量Z（x）在x轴方向上的变异函数，记为
在一维条件下变异
x+h处的值Z（x）与
（h），即
在二阶平稳假设条件下，对任意的h有，E[Z（x+h）]=E[z（x）]。
因此上式可以改写为：
从上式可知，变异函数依赖于两个自变量
x和
h，当变异函数
（x,h）
仅仅依赖于距离
h
而与位置x无关时，可改写成（h），即
设Z（x）是系统某属性Z在空间位
x处的值，Z（x）为一区域化随机变量，并满足二阶平稳
假设，h为两样本点空间分隔距离，
Z（xi）和Z（xi+h）分别是区域化变量在空间位置
xi和xi+h
处的实测值[i=1,2,...,N（h）]，那么根据上式的定义，变异函数
（h）的离散公式为：
变异函数揭示了在整个尺度上的空间变异格局，而且变异函数只有在最大间隔距离
处才有意义。
DPS数据处理系统提供的求二维实验变差函数值的功能。
在DPS
1/2
中，可以计
算0，45，90和135角方向上的实验变差函数。
实验变差函数的滞后距h，在程序运行时，
首先由程序根据样本点分布情况，分别对X坐标和Y坐标进行排序,自动算出。
由屏幕提示、
用户可以修改。
在实际工作中，因观测点并不完全呈规则的矩形网格分布，所以在求某个方向上的实
验变差函数时，由于观测点不一定完全位于这个方向的同一条直线上，本程序采用角度允
许误差限和距离允许误差限的方式进行调整。
角度允许误差限在数据处理程序中给定。
若角
度允许误差限为10。
要计算90方向上的实验变差函数，则从某一点出发，对位于80—100
之间扇形区域内的任—点都可以看成是在90方向上的点。
距离允许误差限取为基本滞后
距的一半，若基本滞后距为400m，要求两点问的距离为1000m的实验变差函数值，则接点
相距800—1200m范围内的任—点都可以近似看成是这两点间的距离为1000m。
DPS系统中计算半变异函数的操作是在电子表格中，一行一个样本，每个样本包含经
度（X轴方向）、纬度（Y轴方向）和相应的观察值。
然后选定为数据块（图1）。
A
B
C
1
2
3
4
5
6
经度
119.1
119
119.4
119.6
119.5
50-纬度
22
22.4
21.6
22.1
21.9
观察指标
1.12
0
12.48
0
3.16
67
121.3
20.9
30.01
68
121.8
20.6
54.66
69
121.4
21.4
83.42
图1实验半变异函数计算数据编辑、定义格式
计算时，系统将会给出如图2所示的用户操作界面。
在这里，用户可以自行调整有关参数，如基本滞后距离，距离允许误差和角度允许误差。
调整好之后，双击“画图”按钮，就
可以得到当前设置下的半变异函数，并以图形方式显示出来。
图2半变异函数计算的用户工作界面
最后确认后，点击图2右上角的关闭按钮，即返回操作的电子表格界面，并输出结果如下：
数据个数=69平均值=32.6932方差=779.1324
基本滞后=0.45
距离偏差限=0.23
角度偏差限=22.50
半变异函数计算结果
Alpha=
0
No.
滞后距离
Variogram
数据对数
1
0.20000
275.01050
3
2
0.48918
299.85552
58
3
0.92052
473.26980
99
4
1.36348
585.90474
105
5
1.79573
601.68676
108
6
2.25449
612.95724
76
7
2.68121
519.50553
56
8
3.10846
873.19956
22
Alpha=
45
No.
滞后距离
Variogram
数据对数
1
0.19621
56.35225
12
2
0.49794
242.56537
48
3
0.93135
397.98511
87
4
1.35429
508.27885
113
5
1.79470
351.03178
93
6
2.22218
327.75980
70
7
2.61304
355.80146
24
Alpha=
90
No.
滞后距离
Variogram
数据对数
1
0.15000
236.22365
2
2
0.48344
357.32835
63
3
0.92338
674.38745
104
4
1.35888
789.73275
107
5
1.79870
947.65614
120
6
2.26608
1102.45099
92
7
2.67353
1678.68498
69
8
3.09238
1789.68990
40
9
3.52657
2738.54272
14
Alpha=
135
No.
滞后距离
Variogram
数据对数
1
0.20013
108.98754
7
2
0.51359
251.43343
60
3
0.91273
528.13012
91
4
1.36563
748.56560
135
5
1.82181
950.57698
135
6
2.24840
1349.65597
120
7
2.67823
1739.60956
102
8
3.14662
2016.74519
57
9
3.61276
1111.98072
20
二、协方差函数及相关系数
协方差又称半方差，是用来描述区域化随机变量之间的差异的参数。
在概率理论中，随机向量X1与X2的协方差被定义为：
区域化变量
Z（x）=Z（xu,xv,xw）在空间点x和x+h处的两个随机变量
Z（x）和Z（x+h）的二阶
混合中心矩定义为
Z（x）的自协方差函数，即
CovZ（x）,Z（xh）E[Z（x）Z（xh）]
E[Z（x）]E[Z（x
h）]
区域化变量Z（x）的自协方差函数也简称为协方差函数。
一般来说，它是一个依赖于空间
点x和向量h的函数。
设Z（x）为区域化随机变量，并满足二阶平稳假设，即随机函数Z（x）的空间分布规律不因
位移而改变，h为两样本点空间分隔距离或距离滞后，Z（xi）为Z（x）在空间位置xi处的实测值，Z（xi+h）是Z（x）在xi处距离偏离h的实测值（i=1,2,3,,N（h）），根据协方差函数的定义公式，可
得到协方差函数的计算公式为：
C#（h）
1
N（h）
[Z（x
）
Z（x
）][Z（xh）Z（x
i
h）]
N（h）
i1
i
i
i
在上面的公式中，
N（h）是分隔距离为
h时的样本点对的总数，
Z（xi）和Z（xi
h）分别
为Z（xi）和Z（xi+h）的样本平均数，即
Z（xi
1
N
Z（xi
h）
1N
）
Z（xi）
Z（xih）
Ni1
Ni1
在公式中N为样本单元数。
一般情况下
Z（xi）
Z（xi
h）（特殊情况下可以认为近似相
等）。
若Z（xi）Z（xih）m常数），协方差函数可改写为如下：
C#（h）
1
N（h）
m2]
[Z（x）Z（x
i
h）
N（h）
i
i1
式中，m为样本平均数，可由一般算术平均数公式求得，即
m
1
N
Z（xi）。
由于协方差与
Ni1
相关系数有密切的关系，因此通过协方差函数可以导出相关函数的计算公式为
#（h）C#（h），式中C#（0）为先验方差，可以根据公式
C#（0）
1
N
[Z（xi）]2
m2]
C#（0）
N（h）i
1
计算，式中m和上面协方差计算中的相同，为样本平均数，
N为样本单元数。
协方差函数的计算，
和前面的半变异函数计算的数据格式相同。
实际上，在计算协方差
函数的同时，也计算出了半变异函数的计算结果。
三、变异函数理论模型的最优拟合
前面计算的变异函数是实验变异函数（Experimentalvariogram）。
该变异函数可用变异曲
线来表示。
它是一定滞后距a的变异函数值（h）与该h的对应图，一个理想化的变异曲线应如图3所示
图3变异曲线示意图
图3中的曲线包括如下几部分:
Co称为块金效应（nuggeteffect），它表示h很小时两点间观察指标的变化；a称为变程（range）,当ha时，任意两点间的观测值有相关性，这个相关性随h的变大而减小，当h>a时就不再具相关性，a的大小反映了研究对象中某一区域化变量
（如品位）的变化程度，从另一个意义看，
a反映了影响范围，估计
C称为总基台值，它反映
某区域化变量在研究范围内变异的强度，它是最大滞后距的可迁性变异函数的极限值。
当趋于无穷大时：
h
（
）=C（0）=Var[Z（x）]=C
即当h趋于无穷大时，变异函数值近于先验方差C，当有块金效应时，C＝C+C。
；而C称为基台值差：
C＝C-C0。
C（0），当无块金效应（常数）C。
时C＝（sill），它是先验方差与块金效应（常数）之
尽管变异函数有助于解决一个矿床某区域化变量（例如品位）的变化特征及结构性状，但
它纯粹是一个数据的概括技术，当定量地描述全矿床的特征时，有关全矿床的变异结构还必
须借助于推断，这个过程类似于用样品值构制直方图，再从直方团推断全矿床的理论分布。
如果已绘制了试验半变异曲线，那么，为了要得到最后的结论，就必须给试验半变异曲线配
以相应的理论模型，这些理论模型将直接参与克立格计算或其他地理统计学研究。
如同经典统计学那样，理论变异函数也几个简单的模型。
DPS提供了如下6个常用的
理论模型，其基本公式简介如下：
1、球状模型：
一般公式为：
该模型在原点处（h＝0），切线的斜率为比3C/2a，切线到达C值的距离为2a/3。
2、高斯模型：
其公式为：
3、指数模型：
4、球状套合模型：
5、球状＋指数套合模型
0
h
h
0
C（3h
1h3
C
）C
（1ea1
）
0ha
（h）
0
1
2a1
2a13
h
2
1
C
CC（1ea1）
ah3a
0
1
2
1
1
C0
C1
C2
h
3a1
6、线性有基台模型
C0
h
0
（h）C0
Ah
0h
a
C0
C
h
a
在DPS中，上述半变异函数理论模型参数采用加权的非线性最小二乘法
（麦夸特法）进
行估计，并在屏幕上绘出理论变差函数拟台曲线图
（图4）。
理论半变异模型的统计检验，可
类似线性回归分析，将总平方和分解成残差平方和和回归平方和两部分，即
n
n
?
2
n
i
2
i
?
2
i
i
i
i1
i
1
i
1
这样可以构造一个F统计量
n
?
i
2/（k
1）
i
F
i
1
n
?
i
2/（n
k）
i
i1
进行统计检验，并可计算出拟合结果的决定系数
R2。
这些指标都在计算过程中及时反
馈到用户的工作界面（图4），用户可以根据理论变差函数曲线的拟合情况及有关统计检验结果，选择适当的半变异函数模型、直到得出理想的理论变差函数。
图4理论变异函数拟合数据编辑、定义及操作的用户界面
四、交叉验证（CrossValidation）
由于实验变异函数并不是连续的图形，且无法满足条件半正定的性质，因此在实际应用
用时必須套配成连续性的理论变异函数模式后方可应用。
理论变异函数模式的套配即是利用
数值方法来选取固定模式下的最佳参数，因此为了验证所选取的变异函数模型的优劣及利用
克立格估计法推估时所做的假設是否合理，可以利用交叉验证法（CrossValidation）来进行检
定。
半变异函数模型的验证往往因为观察资料的数量太少，
故常常采用用交叉验证法來进行
检定，其步骤如下：
（1）利用所有n个观察资料可求得一代表性的半变异函数模型。
（2）自观察资料点中任意取出一观察值Z1，以其余n-1个观察资料来进行克立格推估，
则可求得该观察点Z1的估计值。
（3）计算该观察点的估计误差及克立格变异方差。
（4）将该点Z1放回，并取出另一个观察点，记为Z2；重复上述步骤，直到所有观察点均
完成克立格估计及其估计误差与克立格方差的计算。
（5）利用上述计算结果分別检定其无偏性与一致性。
无偏性与一致性的判断准则如下：
检定研究中所配套的理论半变异函数模型的优劣及利用克立格估计法推估时所做的假定是否合理，可利用推估结果的无偏性及一致性来检定。
（1）无偏性﹝Unbiasness﹞
当克立格平均误差，即估计值与观察值的期望值愈接近0时，表示其无偏估计的效果
越佳。
无偏性（KAE）的计算公式为：
KAE
1
n
n
Zi*Zi0
i1
（2）一致性﹝Coherence﹞
若克立格估计误差的平方与克立格方差的比值的期望值越趋近于1时，则表示模型愈
佳。
其一致性（KRMSE）的计算公式为：
KRMSE
1
n
n
*
Zi
2
Zi
1
i1
2
k
该程序根据上述原理，用普通克立格法、对观测点上的数据进行估值、求出每个观测
点上的估计慎、估计方差和偏差〔估计值与观测值之差〕，并绘出偏差直方图、及当前的KAE值和KRMSE值。
并可调整有系统自动优化有关模型参数，使之更为接近。
在搜索信息点时，采用等分象限的技术，即将每个象限分成若干份，在每一份中最多
只取n0个点作为信息点参与估值。
程序运行时、用尸可以根据屏幕提示，确定估值的搜索
方位数（4、8、12或16），输入象限等份（角）数（1、2、3或4）。
即将每个象限—等份、表示
4
方位搜索；如将每个象限二等份，表示
8方位搜索。
程序设计成最多可以将每个象限
4
等份，即16方位搜索。
DPS系统中进行交叉验证分析，其操作是在电子表格中，一行一个样本，每个样本包
含经度（X轴方向）、纬度（Y轴方向）和相应的观察值，然后选定为数据块
（数据块的左边一列
可以放各个样本点的名称，
但不要用鼠标选进来，
系统会自动将名称放到交叉验证的结果里
面去的）。
如果将半变异函数模型参数放在电子表格里的话，第一行放
X轴的值，第二行放
Y轴的值。
每行分别放入参数
a,C和C0的值，然后在按下Ctrl键的状态下用鼠标选定。
计算时，系统将会给出如图所示的用户操作界面
（图5）。
在这里，用户可以自行调整有
关参数，调整好之后，双击“理论模型优化”按钮，就可以得到当前设置下的新的理论半变
异函数参数、交叉验证残差直方图、及无偏性指标
KRE和一致性指标
KRMSE的值。
图5
交叉验证分析的数据编辑、定义及操作的用户界面
完成调整、模型优化之后，点击用户操作界面右上方的关闭按钮，输出结果来看，模型优化后，能较好地改进模型参数估计的无偏性，如
即输出结果如下。
从
KRE指标从-0.5015
下降到0，残差方差从244.17下降到222.44（图6）。
图6交叉验证分析结果
最后输出的结果如下：
变差函数参数
球形模型
X方向:
a=2.17
基台c+c0=1163.4246
块金=1051.7414
Y方向:
a=2.07
基台c+c0=1163.4939
块金=1051.8239
套合:
a=2.17
基台c+c0=1163.4592
块金=1051.7827
每个象限分
1等分角
样本个数=69
偏差均值KRE=0.0000
偏差标准差=222.4425
一致性系数KRMSE=0.9855
克立格估计结果
序号
样本号
Zi值
估计值
估计方差
偏差
参与估计值的样本
1
龙泉市
1.1200
5.13473
1438.74238
-4.01473
云和县
庆元县
松阳县
2
庆元县0.0000
0.56913
1641.75198
-0.56913
景宁县
龙泉市
3
松阳县
12.4800
10.43051
1344.99388
2.04949
丽水市
云和县
遂昌县
武义县
4
景宁县
0.0000
2.20430
1343.25375
-2.20430
青田县
泰顺县
龙泉市
云和县
5
云和县
3.1600
9.82873
1337.65271
-6.66873
丽水市
景宁县
龙泉市
松阳县
67
三门县
30.0100
32.87049
1338.15438
-2.86049
象山县
临海市
天台县
宁海县
68
象山县
54.6600
56.52749
1453.87102
-1.86749
椒江市
宁海县
鄞县
69
椒江市
83.4200
44.37612
1428.12143
39.04388
温岭市
黄岩市
三门县
五、克立格插值
克立格法是将任一个点的估计值通过该点影响范围内的n个有效样本值Z（xi）的线性组
合得到，即
式中，i是与样点观察Z（xi）有关的加权系数，用来表示各个样点值
Z（xi）对估计值Zv*的
贡献。
对于任一给定的区域和数据信息
Z（xi），i=1,2,,n，存在一组加权系数I。
如果使得估
计方差为最小，其区域内的真值就能在最小的可能置信区间内产生。
这种给出最佳线性无偏
估计权系数的即为克立格法。
在DPS系统中进行克立格插值，其操作是在电子表格中，一行一个样本，每个样本包
含经度（X轴方向）、纬度（Y轴方向）和相应的观察值。
然后选定为数据块（图1）。
计算
时，系统将会给出如图7所示的用户操作界面。
在这里，用户可以自行调整有关参数，如估计所用的信息点数，插值后生成的数据矩阵的行、列数，并以图形方式显示出来。
图7克立格插值用户操作界面（a.估计值）
图7克立格插值用户操作界面（b.估计方差）
升级会员 