空间分析-第6讲：空间插值与地统计

1、第六讲 空间插值与地统计引言一阶效应方法距离反比加权(IDW)插值趋势面分析二阶效应方法区域化变量变异函数地统计：克里金(Kriging)插值普通克里金(Ordinary Kriging, OK)简单克里金(Simple Kriging, SK) 一、引言复习：1.2节的空间数据类型：对象观点v.s.场观点2.3节，尤其是图2.3、2.4和2.5，以及与近邻多边形图2.3 一组点事件的近邻多边形图2.4 近邻多边形的构建. 多边形的边都是两点之间连线的垂直平分线.图2.5 从近邻多边形导出Delaunay三角网场及其用途的一个例子地理学中场数据的最好例子是地表高度。这种数据可用于：制作地图和其

2、他的地形可视化，用于导航和一般应用在水文学中，如河流的流域与泄流网络在生态学中，如计算坡度和坡向移动电话公司或军队的无线电和雷达传播研究在仿真模拟中，如飞行模拟在地面制导的导航系统中操纵飞行器，如巡航导弹在园林建筑学中，图示可视域在地理科学中，作为一种系统输入，用于预测能量到达（例如，太阳光）和穿越（例如，水）地表的连续变化地面高度场切实存在你几乎一直站在它上面。 空间连续表面（场）：抽样与存储方法 教材p.166-169随机抽样非随机抽样简单随机抽样系统抽样分层抽样分群抽样多步抽样配额抽样雪球抽样便利抽样自选抽样目的抽样极端抽样同质抽样不均匀抽样典型抽样关键抽样空间插值 教材p.169Est

3、imating a point here: interpolation样本点A：内插Estimating a point here: extrapolationB：外推空间插值常用于将测量的离散点数据转换为连续的数据面，包括空间内插和外推两种算法。空间内插：通过已知点的数据推求同一区域未知点数据。空间外推：通过已知区域的数据，推求其它区域数据。Elevation profileSample dataIf A = 8 feet and B = 4 feet thenC = (8 + 4) / 2 = 6 feet线性内插ABC？Elevation profileABC？非线性内插Often re

4、sults in a more realistic interpolation but estimating missing data values is more complex全局内插局部内插Uses a neighborhood of sample points to estimate a value at an un-sampled locationUses all known sample points to estimate a value at an unsampled location内插：全局 & 局部插值插值 空间内插的理论基础：地理学第一定律(FLG)Everything

5、 is related to everything else, but near things are more related than distant things. Generally true with discrete data; Definitely true with continuous data鉴于徒手插值的困难，我们很自然地要问，是否有可能设计计算机算法，以一种一致和可重复的方式进行插值？下述考虑未知位置上的场值预测问题是有用的：如果没有关于控制点位置的信息，关于新样本的可能值，你将做出怎样的估计？基础统计学告诉我们，最好的估计是样本点数据的简单平均值。用空间术语，这相当于

6、图9.5所示的情况。我们做出一种假设，所有未知的场高度值具有等于平均值的相同值，因而它们在研究区域上形成一个单一的水平面。插值(1): 近邻多边形对采用均值进行预测的一种简单改进是应用近邻多边形，赋予每个未抽样点最近控制点的值。用地理学第一定律来说，这意味着我们把邻近概念表示为最近的意思。该方法通过构建控制点位置的近邻多边形而实现，然后假设每个多边形都具有等于控制点值的相同高度值。该技术在几乎一百年以前由Thiessen (1911)介绍，他想用雨量站的记录值估计区域的总降雨量。插值(2):局部空间平均计算样本点数据的局部空间平均。事实上这种插值方法是一种局部统计量。这里的思想是假设地理学第一

7、定律是合理的，并应用邻近位置的均值预测未抽样位置上的值。关键问题是：哪些位置是邻近的？ 近邻多边形方法已经提供了一种答案应用最近的位置。我们可以采用到待确定位置某一固定距离内的点，而不是仅应用最近的控制点。限制半径和最近邻插值有两个共同特征：第一，无论是距离还是近邻数量，限度的选择是主观的。第二，当估计所依赖的值的集合增大时，插值表面中的不连续阶梯变得比较小，产生一种比较平滑的外观。但平滑外观是不真实的，因为控制点仍然在局部均值的计算中突然进入和退出。其中不好的一面是难以在得到的插值结果表面画等高线。而且，用于估计的控制点集越大，插值表面变得越像水平面。稍作思考，那些熟悉这种概念的读者应该能看

8、到，这是经典统计学中心极限定理的空间形式。一、一阶效应方法距离反比加权（IDW）插值趋势面分析距离反比加权(IDW) 插值方法假定每个输入点都有着局部影响，这种影响随着距离增加而减弱。步骤:计算未知点到所有点的距离；计算每个点的权重： 权重是距离倒数的函数。计算结果：1. 距离反比加权(IDW)插值Sum of =14 is the highest6 is the lowest示例：IDW插值（求图中0点的值）示例：IDW插值IDW插值的一般模型 The values are computed for all or selected iThe value for alpha is typica

9、lly 1 or 2The division by the sum of weighted distances ensures that the weights add up to 12022/7/29河南大学环境与规划学院 zhaoy 432100160200例：倒距离加权(IDW)插值结果IDW插值的缺点IDW cannot make estimates above the maximum or below the minimum sample values. For an elevation surface, this has the effect of flattening peaks

10、 and valleys (unless their high and low points are part of the sample). Because the estimated values are averages, the resulting surface will not pass through the sample points. Source dataIDW : alpha=2, no smoothing示例： IDW插值结果IDW插值效果Bulls eye effecta = 2grey: a = 0.1 red: a = 2green: a = 10 red: a

11、= 2美国爱达何洲（Idaho）105个气象台30年平均降水量。中部和西南地区为山区，降水量小。示例IDW插值在ArcMap中的实现IDW插值在ArcMap中的实现Idaho州降雨量等直线图IDW插值在ArcMap中的实现2 趋势面分析趋势面分析(Trend Surface Analysis，TSA)：用多项式方程近似地拟合已知数据的点（Bailey等1995），概括出场的主要特征或趋势描述地理数据的空间分布及其区域性变化的方法。趋势面分析基本思想趋势面分析本质上是多元线性回归分析的简单扩展。基本思想是用多项式表示线或面，按最小二乘法原理对数据点进行拟合，拟合时假定数据点的空间坐标X、Y为独立

12、变量，而表示特征值的Z坐标为因变量。deterministic interpolationfunction fitting: polynomialse.g y = ax2+bx+clocal (piecewise)globallocal 1st orderglobal 1st order趋势面分析：示例1把坐标作为两个自变量，完全按照上述的多元线性回归进行。因此，增广矩阵X 为即，每个点上的残差是该点上的表面观测值与拟合值之差。残差图是探索数据的一种有效方式，以揭示趋势面中没有包含的局部因素。序号降水量z(mm)横坐标x(104m)纵坐标y（104m）132.81.042.13241.20.2

13、3.50326.73.431.72424.62.850.48551.20.650.27638.82.683.82742.73.522.65835.41.642.95931027.93.270.851138.00.381.121240.11.923.191329.60.742.141446.62.650.981537.50.273.651671725.72.412.861801948.93.850.462042.80.641.63趋势面分析：示例 2某地区降水量及地理位置坐标某地区降水量二、三次次多项式趋势面将自变量x和y的值代入

14、以上两个趋势面模型，得到一系列趋势值z，输出结果见下图。根据表1中的数据，运用最小二乘法原理构造趋势面方程，得到二次、三次趋势面多项式分别如下：某地区降水量二、三次次多项式趋势面多项式方程描述趋势面因为任何函数在一定范围内总可以用多项式来逼近，并可调整多项式的次数来满足趋势面分析的需要，一般来说，多项式的次数越高则趋势值越接近于观测值，而剩余值越小。通常使用的是一次、二次、三次趋势面，过高次的趋势面不利于反映空间趋势，并可能存在趋势面的畸变。多项式趋势面的数学模型趋势面分析在ArcMap中的实现趋势面分析在ArcMap中的实现趋势面分析在ArcMap中的实现Idaho州降雨量等直线图无论趋势面

15、分析的优点是什么，它显然是一种相对不明智的方法：通常没有令人信服的理由去假设感兴趣的现象以这种简单的形式进行变化，不管是用空间坐标，或者甚至用坐标平方、立方的组合方式。在实践中，残差几乎总是存在空间自相关。这表明我们的模型是误设的，暗示不能使用该拟合可靠地对结果的统计显著性进行解释。虽然通过最小二乘多元回归，控制点数据用于所选定趋势模型的拟合，但除了它们所呈现的模式的简单可视化外，数据对模型选择没有提供什么帮助。三、二阶效应方法区域化变量平方根差云图和半变异函数地统计：克里金（Kriging）插值随机变量随机函数随机过程随机场区域化变量与时间有关的随机函数带有多个(2个以上)自变量的随机函数以

16、空间点的三个直角坐标为自变量区域化变量当一个变量的取值与其空间位置有关时，就称为区域化变量（regionalized variable）。区域化变量常常反映某种空间现象的特征，用它来描述的现象称之为区域化现象。区域化变量，Matheron(1963)将它定义为以空间点x的三个直角坐标为自变量的随机场区域化变量具有两个最显著，也是最重要的特征：随机性和结构性。http:/cg.ensmp.fr/Presentation/Matheron/Matheron_en.shtml Professor Georges Matheron (1930-2000.8.7)法国数学家和地质学家 区域化变量功能由于

17、区域化变量是一种随机函数，因而能同时反映空间变量的结构性和随机性。一方面，当空间点x固定后，Z(x)就是一个随机变量，这体现了其随机性。另一方面，在空间两个不同点x与x+h处的区域化变量值具有某种程度的相关性，这体现了其结构性。区域化变量的组成部分 数据点结构性 可以用均值和常数趋势表示空间相关 数据通常呈现正空间相关性 随机性 测量误差，其他误差 对于具有n 个抽样控制点的表面，描绘其空间自相关的一种自然方法是根据控制点对之间的距离差异和对应的高度值差异进行做图。精确地做到这一点的是方差云图（variogramcloud）。这些数据是在一个310英尺310英尺的调查区域中收集到的控制点(高程

18、点)。其中加入了手绘等高线，以给出对整体空间结构的一些认识。坡度从北向南有整体上升的趋势，而地图区域的南部在某种程度上更混乱一些。2、平方根差云图和半变异函数p.195对这些数据中每一个可能的点对(有52个高程点，因此有1326个点对)，对他们之间的距离与高度差异的平方根进行做图，得到平方根差云图(square root differences cloud)，如图10.4所示。半方差/半变异函数(semi-variogram) 区域化变量的基本研究工具 半变异函数就是区域化变量增量平方的 数学期望之半区域化变量在i、i+h点的值步长为h的样品对数步长(h)：在一定方向上，距离为h的矢量Vari

19、ance:Variogram:h半变异云图（semivariogram cloud）分组(Binning)如果数据很大，样点对的数目将迅速增加并且变得难以操作。 因此，可以将样点对分组，也就是步长(h)分组。如，可将距离在01之间的作为第一组；h在12之间的作为第二组；依此类推。 Binning is a process that averages semivariance data by distance and direction (hey, its weighted!)Binned semivariogram例：假设某地区降水量Z(x)（单位：mm）是二维区域化随机变量，满足二阶平稳假设

20、，其观测值的空间正方形网格数据如图1所示（点与点之间的距离为h=1 km）。试计算其南北方向及西北和东南方向的变异函数。图1 空间正方形网格数据（点间距h=100m） 从图1可以看出，空间上有些点，由于某种原因没有采集到。如果没有缺失值，可直接对正方形网格数据结构计算变异函数；在有缺失值的情况下，也可以计算变异函数。只要跳过缺失点位置即可（图2）。首先计算南北方向上的变异函数值，由变异函数的计算公式可得 =385/72=5.35 图2 缺失值情况下样本数对的组成和计算过程为缺失值 同样计算出最后，得到南北方向和西北东南方向上的变异函数（下表）。同样可以计算东西方向上的变异函数。 方向：南北 方

21、向：西北东南 h1h2h3h4h5hh1.41h2.82h4.24h5.65h7.07hN(h) 36 27 21 13 5 N(h) 322113825.359.2617.5525.6922.907.0612.9530.8558.1350.00semivariogram南北向semivariogram西北东南向交叉验证（Cross Validation）Re-estimate the samples to find errors in the model Error Statistics Error Histogram Erro Spatial diagram observed x esti

22、mated valueOK?Variogram ModelYesNO?12345IDW插值和趋势面方法的缺陷IDW： 距离权重函数的选择和邻居的定义是其致命缺陷（Achilles heel）。趋势面分析：In a sense, trend surface analysis lets the data speak for themselves, whereas IDW interpolation forces a set structure onto them.以某种方式结合这两种方法是明智之举。通过应用一种距离加权方法，但同时采用可能有助于函数、权重和邻域选择的最好方式,让样本数据自己说话。3

23、、克里金（Kriging）插值Kriging is a group of geostatistical techniques to interpolate the value of a random field (e.g., the elevation, z, of the landscape as a function of the geographic location) at an unobserved location from observations of its value at nearby locations.The theory behind interpolation a

24、nd extrapolation by Kriging was developed by the French mathematician Georges Matheron based on the Masters thesis of Daniel Gerhardus Krige, the pioneering plotter of distance-weighted average gold grades at the Witwatersrand reef complex in South Africa. The English verb is to krige and the most c

25、ommon noun is kriging.克里金（Kriging）插值法，又称空间局部估计或空间局部插值法，是地统计学的主要内容之一。克里金方法适用的条件是，如果变异函数和相关分析的结果表明区域化变量存在空间相关性。克里金方法建立在变异函数理论及结构分析基础之上，利用区域化变量的原始数据和变异函数的结构特点，对未采样点的区域化变量值进行线性无偏最优估计。 本质上，克里金方法所有要做的就是把控制点数据作为一个样本，寻求每个未知位置插值中所包含数据值的最优权重。克里金(Kriging)插值可以分为三步：第1步：描述样本控制点数据中的空间变异经验半变异函数分析。第2步，用适当的数学函数概括这种空间

26、变异理论半变异函数。第3步：应用理论半变异函数模型确定插值权重，对未知点值进行预测。第1步：经验半变异函数：已经完成！第2步：理论半变异函数 半变异函数必须拟合一个数学函数或模型，以用于估计任何给定距离的半方差。理论变异函数模型实践中，常用的是变异函数图：偏基台值:C(partial sill)块金值:C0(nugget)变程: a(range)h基台值(sill)not related anymore变程范围内才有结构性变化（有规律的变化）反映随机性大小：主要来源于区域化变量Z(x)在小于抽样尺度h时所具有的内部变异；另外还有抽样分析误差。变异函数是一个单调不减函数。当h超过某一个范围，例如

27、变程，变异函数不再增大，而是趋于一个极限值，即为基台值。实际上等于区域化变量的先验方差。即，即基台值与块金值之差，表示数据中存在空间相关性引起的方差变化范围。 Spherical ModelC = Co+ C1CohgC1aSill理论半方差模型的类型 CohagC1 Gaussian ModelArcGIS：Geostatistical Analyst provides the following functions to choose from to model the empirical semivariogram: Circular Spherical Tetraspherical P

28、entaspherical Exponential Gaussian Rational Quadratic Hole Effect K-Bessel J-Bessel Stable 理论半方差模型的类型 左侧图是一系列表面的截面，其中属性值的局部变化逐渐增大。右侧图是可能的半变异模型。随着表面局部变异的增加，变程减小且块金值增加。因为在所有情况中数据值的整体变异相似，所以在所有三个半变异中基台值也相似。在这些图中，最明显的区别是变程捕捉数据在空间上局部变异程度的不同。这个简单的模型类似于倒距离权重法(IDW)，不同的是IDW的权重只与距预测点的距离有关，而普通克里金(OK)方法的权重则取决于半

29、变异图、距预测点的距离和预测点周围观测值的空间关系。x1x2x3x4x0Inverse Distance WeightsKrigingLocal Means第3步：应用理论半变异函数模型确定插值权重最优(Best)：要求 尽可能小，从而预测值将尽可能的接近未知值线性(Linear) : weighted linear combinations of the data无偏性(Unbiased)：一些点的预测值比真值高，一些比真值低，平均起来预测值和真值的差为0。要求权重系数值和为1。Estimation 普通克里金（OK）的目标：BLUE包含所有样点对的半变异值包含每一个观测点与预测点之间的半变

30、异值待定基于点i和j之间距离的半变异值基于第i个观测点与预测点之间距离的半变异值普通克里金（OK）的计算过程计算样点对之间的距离&方差，理论半变异值&平均理论半变异值拟合模型计算权系数与预测克里金方差计算样点对之间的距离和方差半变异值计算结果位置距离方差(Zi-Zj)2半变异值经验半变异值 = 1/2方差示例：普通克里金方法距离计算如果数据很大，样点对的数目将迅速增加并且变得难以操作。因此，可以将样点对分组，也就是将步长(h)分组。在该例中，首先将距离在12之间的作为第一组；h在23之间的作为第二组；依此类推（下表）。步长间距样点对间距平均距离半变异值平均半变异值半变异值的步长分组表拟合理论半

31、变异函数模型为了预测，需要用理论半方差模型拟合经验半方差，以描述其变程、基台值和块金(range, sill, & nugget )。 四个常用模型：Gaussian:Linear:Spherical:Exponential:Where: c0 = nuggetb = regression slopea = rangec0+ c = sillAssumes no sill or range用平均半变异值为纵坐标，步长（抽样间距）为横坐标，做理论半变异函数模型图。观测值拟合值距离半变异值理论半变异值 = 斜率距离 = 13.5h生成 A 矩阵：如样点(1,5)与样点(3,4)的半变异值为：13.

32、52.236=30.19A矩阵最后一行的1和0根据无偏估计的限制条件求得。A矩阵的逆矩阵A矩阵A矩阵计算权重系数与预测普通克里金的权重系数矩阵为：向量b由未知点生成。计算未知点(1,4)与每个观测点的距离，如(1,4)与(1,5)、 (3,4)、 (1,3)、 (4,5)、 (5,1)的距离。预测点(1,4)的b向量计算结果如下表。半变异值 = 13.5hb计算权重：权重观测值乘积克里金预测结果正如预期，权系数随距离的增加而减小，且权重之和等于1。但由于将各点的空间分布也考虑进去，其结果比直接的距离权重要好。权系数为负值的处理Deutsch Clayton V. Correcting for

33、negative weights in ordinary krigingJ. Computers & Geosciences, 1996, 22(7): 765-773.克里金方差为了估计预测结果的不确定性，将权系数向量的每一行和g向量的每一行相乘，然后将结果相加，得到预测克里金方差，其平方根就是克里金标准差。b向量权系数乘积克里金方差克里金标准差如果假设误差是正态分布的，则95%的置信区间为：克里金预测值1.96克里金标准差= 95.49, 109.75(w)这仅是一个小例子，但从中可以看出几个重要特点：Kriging是计算密集型的方法。需要一个合适的软件。虽然一些GIS软件提供半方差估计、

34、建模和Kriging，该邻域大多严肃的工作者使用特殊软件如GSLIB(GeostatisticalSoftwareLIBrary)、Variowin，或GS+。所有结果依赖于为从样点数据估计的半方差所拟合的模型，以及相关假设。其中包含一些主观决定（多少个距离组合？拟合的模型是什么？基台值和块金值取多少？）。 不幸的是，甚至一个相当粗略的权重集，当应用于数据时会给出极好的结果 (Chils and Delfiner,1999,p.175; Cressie,1993,p.289-299)。除了所讨论的普通克里金，还有其他不同类型的克里金方法，如当存在均值漂移（drift）时采用Universal

35、Kriging；CoKriging扩展到同时考虑两个或更多变量。Variogram Model Parameters变异函数模型参数对权重结果（如普通克里金，OK）的影响。We now look at how parameters of a variogram (covariance) model affect the (Ordinary Kriging, OK) weights基台值（Sill）, 形状（shape）, 块金值（nugget）, 变程（range）, 和各向异性（anisotropy）基台值（Sill）Sill of 10 vs. 20Sill = 10Sill = 20克里

36、金权重和估计不变，但克里金方差变化相同的比例（如(a)方差是(b)的20/10=2倍）。Exponential vs. GaussianExponential modelGaussian model理论模型形状（Shape）Shape高斯模型赋予较近样本更大权重The Gaussian variogram model assigns more weight to the closer samplesScreen effect - a sample falls behind another sample that is closer to the unknown. It receives les

37、s (or negative) weights, sample 5 vs. 6. 高斯模型比指数模型具有更强的 screen effect Weights that are less than 0 or greater than 1 can produce estimates larger than the largest sample value or smaller than the smallest. Weights within 0,1 produce estimates only within the min and max of sample values Negative wei

38、ghts may produce negative estimates, although in most science applications values should be positive 块金值（Nugget ）Nugget = 0 vs. =1/2 sillThe Nugget EffectThe nugget effect makes weights become more similar to each other and results in higher kriging variance纯块金效应模型表明完全缺乏空间相关性，或在小于最小抽样间距的尺度上才具有空间依赖性。

39、变程（Range ）Range of h vs. 0.5hRange = 10Range = 20The Effect of Range变程的减小导致克里金方差的增大。A decrease in the range raises the kriging variance 如果变程太小，所有样本与待估计点的距离差不多相等。那么估计结果变得与简单平均相似。If the range becomes too small, then all samples appear to be equally far away from the point being estimated. Then the est

40、imation becomes similar to the simple average of the samples with same weight, 1/nAnisotropy Directional variograms and covariance functionsEffect of Anisotropy More weights are given to samples lie in the direction of maximum continuityWeights given to the samples in the maximum spatial continuity

41、would increase as the anisotropy ratio becomes largerKriging在ArcMap中的实现：探索性数据分析Kriging在ArcMap中的实现：探索性数据分析Kriging在ArcMap中的实现：探索性数据分析Kriging在ArcMap中的实现：探索性数据分析Kriging在ArcMap中的实现Kriging在ArcMap中的实现Kriging在ArcMap中的实现Idaho州降雨量等直线图Kriging在ArcMap中的实现KrigingIDW插值趋势面：5次多项式地形插值比较Voronoi图TINIDWKriging样条规则样条克里金方

42、法的基本形式确定趋势项空间相关随机误差项克里金家族的其他成员对误差项的假设：期望值为0，并且 和 之间的自相关不取决于s点的位置，而取决于位移量h。为确保自相关方差有解，必须允许某两点间的自相关可以相等。如下面有箭头相连的两对位置点假设具有相同的自相关性。趋势值 可以被简单地赋予一个常量，即，在任何位置处 。理解不同的克里金模型如果 未知，就是普通克里金模型。如果在任何时候趋势 已知，无论趋势是否是常量，都形成简单克里金模型。趋势也可以表示为： 若趋势中的系数未知，就是泛克里金模型。理解不同的克里金模型可以对Z(s)进行一下变换。例如，可以把它换为一个指示变量，若Z(s)低于一定的阈值（如空气

43、中的溴氧浓度值0.12ppm），则变为0；若Z(s)高于一定的阈值，则变为1。然后对高于阈值的情况进行预测。便形成了指示克里金模型。现在看一下方程 的左边部分。普通克里金（Ordinary Kriging，OK）简单克里金（Simple Kriging，SK）泛克里金（Universal Kriging，UK）指示克里金（Indicator Kriging，IK）协克里金（cokriging）Assumes that and are unknown constantsThe variable of interest is Z1 , and both autocorrelation for Z1

44、 and cross-correlations between Z2 are used.Types of KrigingSimple kriging is optimal estimation of a random field, e.g. T(x), with a known mean, m(x), and a known covariance PTT(x,x).Ordinary kriging is optimal estimation of a random field, e.g. T(x), with an unknown constant or linearly trending m

45、ean, but a known semivariogram gTT(x,x).Universal kriging is optimal estimation of a random field, e.g. T(x), with an unknown polynomial trending mean, but a known semivariogram gTT(x,x).Simple KrigingOdinary KrigingUniversal KrigingOdinary Kriging误差图Universal KrigingSimple KrigingSimple KrigingOdin

46、ary KrigingUniversal KrigingArcMap：普通克里金的四张图Demo： package gstat/ Meuse data set155 samples taken on a support of 10 x10 m from the top 0-20 cm of alluvial soils in a 5x2 km part the floodplain of the Maas (Meuse) near Stein (NL).id point numberx, y coordinates E and N, in mcadmium concentration in t

47、he soil, in mg kg-1copper concentration in the soil, in mg kg-1lead concentration in the soil, in mg kg-1zinc concentration in the soil, in mg kg-1elev elevation above local reference level, in mom organic matter loss on ignition, in percentffreq flood frequency class, 1: annual, 2: 2-5 years, 3: ev

48、ery 5 yearssoil soil class, codedlime has the land here been limed? 0 or 1 = F or Tlanduse land use, codeddist.m distance from main River Maas channel, in mlibrary(gstat)library(sp) #meuse数据data(meuse)meusesummary(meuse)data(meuse.grid)g - gstat(formula=log(zinc)1,locations=x+y,data=meuse,model=vgm(1,Exp,300)x - predict(g, meuse.grid)image(x, 4, main=kriging variance and data points)points(meuse$x, meuse$y, pch = +) # 克里金方差和样本点coordinates(meuse) = x+yplot(variogram(log(zinc)1, meuse, cloud=TRUE) # 云图plot(variogram(log(zinc)1, meuse),col=2) # 经验图x - vario