1、地理统计 (Geo-statistics)简介浙江大学农业与生物技术学院 唐启义（浙江大学华家池校区 杭州 310029）当一个变量呈空间分布时，就称之为区域化变量。这种变量反映了空间某种属性的分布特征。矿产、地质、海洋、土壤、气象、水文、生态、温度、浓度等领域都具有某种空间属性。区域化变量具有双重性，在观测前区域化变量 Z(x)是一个随机场，观测后是一个确定的空间点函数值。区域化变量具有两个重要的特征。一是区域化变量 Z(x)是一个随机函数,它具有局部的、随机的、异常的特征；其次是区域化变量具有一般的或平均的结构性质，即变量在点 X 与偏离空间距离为 h 的点 xh 处的随机量 Z(x)与
2、Z(x+h)具有某种程度的自相关，而且这种自相关性依赖于两点间距离 h 与变量特征。在某种意义上说这就是区域化变量的结构性特征。一、实验半变异函数变异函数又称变差函数、变异矩，是地统计分析所特有的基本工具。在一维条件下变异函数定义为，当空间点 x 在一维 x 轴上变化时，区域化变量 Z(x)在点 x 和 x+h 处的值 Z(x)与Z(x+h)差的方差的一半为区域化变量 Z(x)在 x 轴方向上的变异函数，记为 (h)，即在二阶平稳假设条件下，对任意的 h 有，EZ (x+h)=Ez(x)。因此上式可以改写为：从上式可知，变异函数依赖于两个自变量 x 和 h，当变异函数 (x,h) 仅仅依赖于距
3、离h 而与位置 x 无关时，可改写成 (h)，即设 Z(x)是系统某属性 Z 在空间位 x 处的值，Z(x) 为一区域化随机变量，并满足二阶平稳假设，h 为两样本点空间分隔距离，Z( xi)和 Z(xi+h)分别是区域化变量在空间位置 xi 和xi+h 处的实测值i=1,2,., N(h)，那么根据上式的定义，变异函数 (h) 的离散公式为：变异函数揭示了在整个尺度上的空间变异格局，而且变异函数只有在最大间隔距离 1/2处才有意义。DPS 数据处理系统提供的求二维实验变差函数值的功能。在 DPS 中，可以计算 0， 45，90 和 135角方向上的实验变差函数。实验变差函数的滞后距 h，在程序
4、运行时，首先由程序根据样本点分布情况，分别对 X 坐标和 Y 坐标进行排序,自动算出。由屏幕提示、用户可以修改。在实际工作中，因观测点并不完全呈规则的矩形网格分布，所以在求某个方向上的实验变差函数时， 由于观测点不一定完全位于这个方向的同一条直线上，本程序采用角度允许误差限和距离允许误差限的方式进行调整。角度允许误差限在数据处理程序中给定。若角度允许误差限为 10。要计算 90方向上的实验变差函数，则从某一点出发，对位于80100之间扇形区域内的任点都可以看成是在 90 方向上的点。距离允许误差限取为基本滞后距的一半，若基本滞后距为 400m，要求两点问的距离为 1000m 的实验变差函数值，
5、则接点相距 8001200m 范围内的任点都可以近似看成是这两点间的距离为1000m。DPS 系统中计算半变异函数的操作是在电子表格中，一行一个样本，每个样本包含经度（X 轴方向） 、纬度（Y 轴方向）和相应的观察值。然后选定为数据块（图 1） 。A B C1 经度 50-纬度 观察指标2 119.1 22 1.123 119 22.4 04 119.4 21.6 12.485 119.6 22.1 06 119.5 21.9 3.16 67 121.3 20.9 30.0168 121.8 20.6 54.6669 121.4 21.4 83.42图 1 实验半变异函数计算数据编辑、定义格式
6、计算时，系统将会给出如图 2 所示的用户操作界面。在这里，用户可以自行调整有关参数，如基本滞后距离，距离允许误差和角度允许误差。调整好之后，双击“画图”按钮，就可以得到当前设置下的半变异函数，并以图形方式显示出来。图 2 半变异函数计算的用户工作界面最后确认后，点击图 2 右上角的关闭按钮，即返回操作的电子表格界面，并输出结果如下：数据个数=69 平均值= 32.6932 方差= 779.1324基本滞后=0.45 距离偏差限=0.23 角度偏差限=22.50 半变异函数计算结果 Alpha= 0 No. 滞后距离 Variogram 数据对数 1 0.20000 275.01050 32 0
7、.48918 299.85552 583 0.92052 473.26980 994 1.36348 585.90474 1055 1.79573 601.68676 1086 2.25449 612.95724 767 2.68121 519.50553 568 3.10846 873.19956 22Alpha= 45 No. 滞后距离 Variogram 数据对数 1 0.19621 56.35225 122 0.49794 242.56537 483 0.93135 397.98511 874 1.35429 508.27885 1135 1.79470 351.03178 936 2
8、.22218 327.75980 707 2.61304 355.80146 24Alpha= 90 No. 滞后距离 Variogram 数据对数 1 0.15000 236.22365 22 0.48344 357.32835 633 0.92338 674.38745 1044 1.35888 789.73275 1075 1.79870 947.65614 1206 2.26608 1102.45099 927 2.67353 1678.68498 698 3.09238 1789.68990 409 3.52657 2738.54272 14Alpha= 135 No. 滞后距离 V
9、ariogram 数据对数 1 0.20013 108.98754 72 0.51359 251.43343 603 0.91273 528.13012 914 1.36563 748.56560 1355 1.82181 950.57698 1356 2.24840 1349.65597 1207 2.67823 1739.60956 1028 3.14662 2016.74519 579 3.61276 1111.98072 20二、协方差函数及相关系数协方差又称半方差，是用来描述区域化随机变量之间的差异的参数。在概率理论中，随机向量 X1与 X2 的协方差被定义为：区域化变量 Z(x)=
10、Z(xu,xv,xw) 在空间点 x 和 x+h 处的两个随机变量 Z(x)和 Z(x+h)的二阶混合中心矩定义为 Z(x)的自协方差函数，即 )()(, ExZEhCo 区域化变量 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)，根据协方差函数的定义公式，可得到协
11、方差函数的计算公式为： )()()()(1# xZxZxhNCiiiii 在上面的公式中，N(h)是分隔距离为 h 时的样本点对的总数， 和 分别)(i)(hxZi为 Z(xi)和 Z(xi+h)的样本平均数，即NiiixZx1)()( Niii hxZh1)()(在公式中 N 为样本单元数。一般情况下 （特殊情况下可以认为近似)()xii相等） 。若 常数） ，协方差函数可改写为如下：mhxZii)() )()(12#mhxZhNCii式中，m 为样本平均数，可由一般算术平均数公式求得，即 。由于协方差NiixZ1)(与相关系数有密切的关系，因此通过协方差函数可以导出相关函数的计算公式为，式
12、中 C#(0)为先验方差，可以根据公式)0()(#h )()(102mxZhNii计算，式中 m 和上面协方差计算中的相同，为样本平均数，N 为样本单元数。协方差函数的计算，和前面的半变异函数计算的数据格式相同。实际上，在计算协方差函数的同时，也计算出了半变异函数的计算结果。三、变异函数理论模型的最优拟合前面计算的变异函数是实验变异函数(Experimental variogram)。该变异函数可用变异曲线来表示。它是一定滞后距 a 的变异函数值 (h)与该 h 的对应图，一个理想化的变异曲线应如图 3 所示图 3 变异曲线示意图图 3 中的曲线包括如下几部分:C o 称为块金效应(nugge
13、t effect)，它表示 h 很小时两点间观察指标的变化；a 称为变程(range),当 ha 时，任意两点间的观测值有相关性，这个相关性随 h 的变大而减小，当 ha 时就不再具相关性，a 的大小反映了研究对象中某一区域化变量( 如品位) 的变化程度，从另一个意义看，a 反映了影响范围，估计 C 称为总基台值，它反映某区域化变量在研究范围内变异的强度，它是最大滞后距的可迁性变异函数的极限值。当 h 趋于无穷大时：()=C(0)=VarZ(x)=C即当 h 趋于无穷大时，变异函数值近于先验方差 C(0)，当无块金效应(常数) C。时CC ，当有块金效应时，C C +C。 ；而 C 称为基台值
14、(sill)，它是先验方差与块金效应(常数)之差： CC-C 0。尽管变异函数有助于解决一个矿床某区域化变量(例如品位) 的变化特征及结构性状，但它纯粹是一个数据的概括技术，当定量地描述全矿床的特征时，有关全矿床的变异结构还必须借助于推断，这个过程类似于用样品值构制直方图，再从直方团推断全矿床的理论分布。如果已绘制了试验半变异曲线，那么，为了要得到最后的结论，就必须给试验半变异曲线配以相应的理论模型，这些理论模型将直接参与克立格计算或其他地理统计学研究。如同经典统计学那样，理论变异函数也几个简单的模型。DPS 提供了如下 6 个常用的理论模型，其基本公式简介如下：1、 球状模型：一般公式为：该
15、模型在原点处(h0)，切线的斜率为比 3C/2a，切线到达 C 值的距离为 2a/3。2、高斯模型：其公式为：3、指数模型：4、球状套合模型：5、球状指数套合模型 111332102)(0)(0)( ahahCeehhah6、线性有基台模型 ahAh00)(在 DPS 中，上述半变异函数理论模型参数采用加权的非线性最小二乘法(麦夸特法)进行估计，并在屏幕上绘出理论变差函数拟台曲线图(图 4)。理论半变异模型的统计检验，可类似线性回归分析，将总平方和分解成残差平方和和回归平方和两部分，即 niinini 121212这样可以构造一个 F 统计量 niiiikn12)/(/进行统计检验，并可计算出
16、拟合结果的决定系数 R2。这些指标都在计算过程中及时反馈到用户的工作界面(图 4)，用户可以根据理论变差函数曲线的拟合情况及有关统计检验结果，选择适当的半变异函数模型、直到得出理想的理论变差函数。图 4 理论变异函数拟合数据编辑、定义及操作的用户界面四、交叉验证 （Cross Validation）由于实验变异函数并不是连续的图形，且无法满足条件半正定的性质，因此在实际应用用时必須套配成连续性的理论变异函数模式后方可应用。理论变异函数模式的套配即是利用数值方法来选取固定模式下的最佳参数，因此为了验证所选取的变异函数模型的优劣及利用克立格估计法推估时所做的假設是否合理，可以利用交叉验证法(Cro
17、ss Validation)来进行检定。半变异函数模型的验证往往因为观察资料的数量太少，故常常采用用交叉验证法來进行检定，其步骤如下：(1)利用所有 n 个观察资料可求得一代表性的半变异函数模型。(2)自观察资料点中任意取出一观察值 ，以其余 n-1 个观察资料来进行克立格推估，1Z则可求得该观察点 的估计值。1Z(3)计算该观察点的估计误差及克立格变异方差。(4)将该点 放回，并取出另一个观察点，记为 ；重复上述步骤，直到所有观察点1 2Z均完成克立格估计及其估计误差与克立格方差的计算。(5)利用上述计算结果分別检定其无偏性与一致性。无偏性与一致性的判断准则如下：检定研究中所配套的理论半变异
18、函数模型的优劣及利用克立格估计法推估时所做的假定是否合理，可利用推估结果的无偏性及一致性来检定。(1)无偏性Unbiasness当克立格平均误差，即估计值与观察值的期望值愈接近 0 时，表示其无偏估计的效果越佳。无偏性（KAE）的计算公式为： 1*niiiZKAE(2)一致性Coherence 若克立格估计误差的平方与克立格方差的比值的期望值越趋近于 1 时，则表示模型愈佳。其一致性（KRMSE ）的计算公式为： 12*nikiiZKRMSE该程序根据上述原理，用普通克立格法、对观测点上的数据进行估值、求出每个观测点上的估计慎、估计方差和偏差估计值与观测值之差 ，并绘出偏差直方图、及当前的KA
19、E 值和 KRMSE 值。并可调整有系统自动优化有关模型参数，使之更为接近。在搜索信息点时，采用等分象限的技术，即将每个象限分成若干份，在每一份中最多只取 n0 个点作为信息点参与估值。程序运行时、用尸可以根据屏幕提示，确定估值的搜索方位数(4 、8、12 或 16)，输入象限等份(角) 数(1、2、3 或 4)。即将每个象限等份、表示4 方位搜索；如将每个象限二等份，表示 8 方位搜索。程序设计成最多可以将每个象限 4 等份，即 16 方位搜索。DPS 系统中进行交叉验证分析，其操作是在电子表格中，一行一个样本，每个样本包含经度(X 轴方向 )、纬度(Y 轴方向)和相应的观察值，然后选定为数
20、据块( 数据块的左边一列可以放各个样本点的名称，但不要用鼠标选进来，系统会自动将名称放到交叉验证的结果里面去的) 。如果将半变异函数模型参数放在电子表格里的话，第一行放 X 轴的值，第二行放 Y 轴的值。每行分别放入参数 a, C 和 C0 的值，然后在按下 Ctrl 键的状态下用鼠标选定。计算时，系统将会给出如图所示的用户操作界面(图 5)。在这里，用户可以自行调整有关参数，调整好之后，双击“理论模型优化”按钮，就可以得到当前设置下的新的理论半变异函数参数、交叉验证残差直方图、及无偏性指标 KRE 和一致性指标 KRMSE 的值。图 5 交叉验证分析的数据编辑、定义及操作的用户界面完成调整、
21、模型优化之后，点击用户操作界面右上方的关闭按钮，即输出结果如下。从输出结果来看，模型优化后，能较好地改进模型参数估计的无偏性，如 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
22、 偏差均值 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 青田县 泰顺县
23、 龙泉市 云和县 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 是与
24、样点观察 Z(xi)有关的加权系数，用来表示各个样点值 Z(xi)对估计值 Zv*的贡献。对于任一给定的区域和数据信息 Z(xi)，i=1,2, n，存在一组加权系数 I。如果使得估计方差为最小，其区域内的真值就能在最小的可能置信区间内产生。这种给出最佳线性无偏估计权系数的即为克立格法。 在 DPS 系统中进行克立格插值，其操作是在电子表格中，一行一个样本，每个样本包含经度（X 轴方向） 、纬度（Y 轴方向）和相应的观察值。然后选定为数据块（图 1） 。计算时，系统将会给出如图 7 所示的用户操作界面。在这里，用户可以自行调整有关参数，如估计所用的信息点数，插值后生成的数据矩阵的行、列数，并以图形方式显示出来。图 7 克立格插值用户操作界面（a. 估计值）图 7 克立格插值用户操作界面（b.估计方差）
