地质数据处理_插值方法

二维数据场的插值方法
1．二维数据场描述及处理目的
数据场数据
{(xi,yi,zi), i=1,…,n}, 即某特征在二维空间中的n个预测值列表：
处理目的
了解该数据场的空间分布情况
处理思路
网格化 绘制等值线图
网格化方法：
二维数据插值
2．空间内插方法
Surfer8.0中常用的插值方法
Gridding Methods
Inverse Distance to a Power(距离倒数加权)
Kriging（克立格法）
Minimum Curvature（最小曲率法）
Modified Shepard's Method（改进Shepard方法）
Natural Neighbor（近邻法）
Nearest Neighbor（最近邻法）
Polynomial Regression（多项式回归法）
Radial Basis Function（径向基函数法）
Triangulation with Linear Interpolation（线性插值三角形法） Moving Average（移动平均法）
Data Metrics（数据度量方法）
Local Polynomial（局部多项式法）
Geostatistics Analyst Model in ArcGIS 9
2.1反距离加权插值
反距离加权插值（Inverse Distance Weighting ，简称IDW ），反距离加权法是最常用的空间内插方法之一。

它的基本原理是：空间上离得越近的物体其性质越相似，反之亦然。

这种方法并没有考虑到区域化变量的空间变异性，所以仅仅是一种纯几何加权法。

反距离加权插值的一般公式为：
∑==n
i i i i y x Z y x Z 1),(),(λ
其中，0Z(x )为未知点0x 处的预测值，i Z(x )为已知点i x 处的值，n 为样点的数量，λ为样点的权重值，其计算公式为：
n
p p i i0
i0i 1d /d λ--==∑ 式中i0d 为未知点与各已知点之间的距离，p 是距离的幂。

样点在预测过程中受参数p 的影响，幂越高, 内插的平滑效果越佳。

尽管反距离权重插值法很简单，易于实现，但它不能对内插的结果作精度评价，所得结果可能会出现很大的偏差，人为难以控制。

2.2全局多项式插值（趋势分析法）
根据有限的样本数据拟合一个表面来进行内插，称之为全局多项式内插方法。

一般多采用多项式来进行拟合，求各样本点到该多项式的垂直距离的和，通过最小二乘法来获得多项式的系数，这样所得的表面可使各样本点到表面之间距
离的平方和最小。

),(),(y x f y x Z =
如果表面平滑、无弯曲，使用一次多项式拟合；有一处弯曲的表面则用二次多项式进行拟合；若有两处弯曲则需使用三次多项式，依次类推。

全局多项式内插一般适用于表面变化平缓的研究区域，或者仅研究区域内全局性趋势的情况
[3]。

2.3局部多项式内插
局部多项式内插与全局多项式内插相对应，是用多个多项式拟合表面的一种方法，它更多地用来表现研究区域西部的变异情况。

其基本原理与全局多项式内插相同。

The Local Polynomial gridding method assigns values to grid nodes by using a weighted least squares fit with data within the grid node's search ellipse.
2.4径向基函数方法
径向基函数法属于人工神经网络方法，该方法所拟合的表面都必须经过所有样本数据。

径向基函数以某个已知点为中心按一定距离变化的函数，因此在每个数据点都会形成径向基函数，即每个基函数的中心落在某一个数据点上。

径向基函数适合于非常平滑的表面，要求样本数据量大，如果数据点少，则内插效果不佳[3]。

同时，径向基函数难以对误差进行估计，也是其缺点之一。

常用的径向基函数法，它们分别是：
薄盘样条函数（thin-plate spline ）：2222B(h)(h R )ln(h R )=⨯⨯
张力样条函数（spline with tension ）： 20E R h B(h)ln(
)K (R h)C 2
⋅=+⋅+ 规则样条函数（completely regularized spline ）： n 2n 221E n 1(-1)r R h R h B(h)ln()E ()C n!n 22∞
=⋅⋅⋅=-=++∑ 高次曲面样条函数（multiquadric spline ）
：B(h)=反高次曲面样条函数（inverse multiquadric spline ）
：B(h)=各式子中h 为表示由点(x ，y)到第i 个数据点的距离，R 参数是用户指定的平滑因子，0K 为修正贝塞尔函数，1E 为指数积分函数，E C 为 Euler 常数，其值约为0.577215。

Radial Basis Function interpolation is a diverse group of data interpolation methods. In terms of the ability to fit your data and to produce a smooth surface, the Multiquadric method is considered by many to be the best. All of the Radial Basis Function methods are exact interpolators, so they attempt to honor your data. You can introduce a smoothing factor to all the methods in an attempt to produce a smoother surface.
Function Types
The basis kernel functions are analogous to variograms in Kriging. The basis kernel functions define the optimal set of weights to apply to the data points when interpolating a grid node. The available basis kernel functions are listed in the Type drop-down list in the Radial Basis Function Options dialog.
Inverse Multiquadric
Multilog
Multiquadratic
Natural Cubic Spline
Thin Plate Spline
where:
h is the anisotropically rescaled, relative distance from the point to the node
R2is the smoothing factor specified by the user
Default R2 Value
The default value for R2 in the Radial Basis Function gridding algorithm is calculated as follows:
(length of diagonal of the data extent)2 / (25 * number of data points)
Specifying Radial Basis Function Advanced Options
1. Click on Grid | Data.
2. In the Open dialog, select a data file and then click the Open button.
3. In the Grid Data dialog, choose Radial Basis Function in the Gridding Method group.
4. Click the Advanced Options button to display the Radial Basis Advanced
Options dialog.
5. In the General page, you can specify the function parameters for the
gridding operation.
▪The Basis Function list specifies the basis kernel function to use during gridding. This defines the optimal weights applied to the data points during the interpolation. The Basis
Function is analogous to the variogram in Kriging. Experience indicates that the
Multiquadric basis function works quite well in most cases. Successful use of the Thin Plate
Spline basis function is also reported regularly in the technical literature.
▪The R2Parameter is a shaping or smoothing factor. The larger the R2Parameter shaping factor, the rounder the mountain tops and the smoother the contour lines. There is no universally
accepted method for computing an optimal value for this factor. A reasonable trial value
for R2Parameter is between the average sample spacing and one-half the average sample
spacing.
Triangulation with Linear Interpolation
The Triangulation with Linear Interpolation method in Surfer uses the optimal Delaunay triangulation. The algorithm creates triangles by drawing lines between data points. The original points are connected in such a way that no triangle edges are intersected by other triangles. The result is a patchwork of triangular faces over the extent of the grid. This method is an exact interpolator.
Each triangle defines a plane over the grid nodes lying within the triangle, with the tilt and elevation of the triangle determined by the three original data points defining the triangle. All grid nodes within a given triangle are defined by the triangular surface. Because the original data are used to define the triangles, the data are honored very closely.
Triangulation with Linear Interpolation works best when your data are evenly distributed over the grid area. Data sets that contain sparse areas result in distinct triangular facets on the map.
2.5最小曲率法
Minimum Curvature is widely used in the earth sciences. The interpolated surface generated by Minimum Curvature is analogous to a thin, linearly elastic plate passing through each of the data values with a minimum amount of bending.
The Minimum Curvature gridding algorithm is solves the specified partial
differential equation using a successive over-relaxation algorithm. The
interior is updated using a "chessboard" strategy, as discussed in Press, et al.
(1988, p. 868). The only difference is that the biharmonic equation must have nine different "colors," rather than just black and white.
Minimum Curvature generates the smoothest possible surface while attempting to honor your data as closely as possible. Minimum Curvature is not an exact interpolator, however. This means that your data are not always honored exactly.
Minimum Curvature produces a grid by repeatedly applying an equation over the grid in an attempt to smooth the grid. Each pass over the grid is counted as one iteration. The grid node values are recalculated until successive changes in the values are less than the Maximum Residuals value, or the maximum number of iterations is reached (Maximum Iteration field).
The Maximum Residual parameter has the same units as the data, and an appropriate value is approximately 10% of the data precision. If data values are measured to the nearest 1.0 units, the Maximum Residual value should be set at 0.1. The iterations continue until the maximum grid node correction for the entire iteration is less than the Maximum Residual value. The default Maximum Residual value is given by:
Default Max Residual = 0.001 (Z max - Z min) The Maximum Iteration parameter should be set at one to two times the
number of grid nodes generated in the grid file. For example, when
generating a 50 by 50 grid using Minimum Curvature, the Maximum Iteration value should be set between 2,500 and 5,000.
The Internal Tension and Boundary Tension ,Qualitatively, the Minimum
Curvature gridding algorithm is attempting to fit a piece of sheet metal
through all of the observations without putting any creases or kinks in the surface. Between the fixed observation points, the sheet bows a bit. The
Internal Tension is used to control the amount of this bowing on the interior: the higher the tension, the less the bowing. For example, a high tension makes areas between observations look like facets of a gemstone. The Boundary
Tension controls the amount of bowing on the edges. The range of values for Internal Tension and Boundary Tension are 0 to 1. By default, the Internal Tension and the Boundary Tension are set to 0.
the Relaxation Factor,
The Relaxation Factor is as described in Press et al. (1988). In general, the Relaxation Factor should not be altered. The default value (1.0) is a good
generic value. Roughly, the higher the Relaxation Factor (closer to two) the faster the Minimum Curvature algorithm converges, but the more likely it will not converge at all. The lower the Relaxation Factor (closer to zero) the more likely the Minimum Curvature algorithm will converge, but the
algorithm is slower. The optimal Relaxation Factor is derived through trial and error.
2.6近邻法
The Natural Neighbor gridding method is quite popular in some fields. What is Natural Neighbor interpolation?
Consider a set of Thiessen polygons (the dual of a Delaunay triangulation). If a new point (target) were added to the data set, these Thiessen polygons would be modified. In fact, some of the polygons would shrink in size, while none would increase in size. The area associated with the target's Thiessen polygon that was taken from an existing polygon is called the "borrowed area."
The Natural Neighbor interpolation algorithm uses a weighted average of the neighboring observations, where the weights are proportional to the "borrowed area."
2.7最近邻法
The Nearest Neighbor gridding method assigns the value of the nearest point to each grid node. This method is useful when data are already evenly spaced, but need to be converted to a Surfer grid file. Alternatively, in cases where the data are nearly on a grid with only a few missing values, this method is effective for filling in the holes in the data.
2.8多项式回归方法
You can select the type of polynomial regression to apply to your data from the Surface Definition group. As you select the different types of polynomials, a generic polynomial form of the equation is presented in the dialog, and the values in the Parameters group change to reflect the selection. The available choices are:
Simple planar surface\Bi-linear saddle\Quadratic surface\Cubic surface\User
defined polynomial
The Parameters group allows you to specify the maximum powers for the X and Y component in the polynomial equation. As you change the Parameters values, the options are changed in the Surface Definition group to reflect the defined parameters.
The Max X Order specifies the maximum power for the X component in the polynomial equation.
The Max Y Order specifies the maximum power for the Y component in the polynomial equation.
The Max Total Order specifies the maximum sum of the Max X Order and Max Y Order powers. All of the combinations of the X and Y components are included in the polynomial equation as long as the sum of the two powers does not exceed the Max Total Order value.
Data Metrics
The collection of data metrics gridding methods creates grids of information about the data on a node-by-node basis. The data metrics gridding methods are not, in general, weighted average interpolators of the Z-values. For example, you can obtain information such as:
the number of data points used to interpolate each grid node. If the number of data points used are fairly equal at each grid node, then the quality of the grid at each grid node can be interpreted.
the standard deviation, variance, coefficient of variation, and median
absolute deviation of the data at each grid node. These are measures of the variability in space of the grid, and are important information for statistical analysis.
the distance to the nearest data point. For example, if the XY values of a data set are sampling locations, use the Distance to Nearest data metric to determine
locations for new sampling locations. A contour map of the distance to the nearest data point, quantifies where higher sampling density may be desired.
Data metrics use the local data set including breaklines, for a specific grid node for the selected data metric. The local data set is defined by the search parameters. These search parameters are applied to each grid node to determine the local data set. In the following descriptions, when computing the value of a grid node (r,c), the local data set S(r,c) consists of data within the specified search parameters centered at the
specific grid node only. The set of selected data at the current grid node (r,c), can be represented by S(r,c), where}
yi
zi
xi
r
=
i
c
S=
,(n
1
,...,
,
),
,
)
{(
where n = number of data points in the local data set
The Z(r,c) location refers to a specific node within the grid.
There are five groups of data metrics, Z Order Statistics, Z Moment Statistics, Other Z Statistics, Data Location Statistics, and Terrain Statistics.
Moving Average
The Moving Average gridding method assigns values to grid nodes by averaging the data within the grid node's search ellipse.
Search Ellipse
To use Moving Average, define a search ellipse and specify the minimum number of data to use. For each grid node, the neighboring data are identified by centering the search ellipse on the node. The output grid node value is set equal to the arithmetic average of the identified neighboring data. If there are fewer than the specified minimum number of data within the neighborhood, the grid node is blanked.
2.9多项式回归方法
2.6克立格内插方法
2.5.1基本理论
克立格（Kriging ）内插方法最早始于采矿领域，它将传统矿产储量计算方法和统计学结合起来，对空间域或时空域上的区域化变量（Regionalized Variable ）进行研究，其研究的主体内容是对数据的结构分析和插值表面的建立及质量评价。

2.5.1.1 基本假设条件
介绍地质统计学中结构分析研究工具——变异函数之前，先来看看地质统计学中提出的两个基本假设条件[1]：（1）二阶平稳假设（Second-Order Stationary ），它必须满足两个条件，一是研究区域内区域化变量的数学期望存在，且等于常数，二是区域化变量的协方差存在且相同，即与变量位置无关，仅依赖于变量间的距离，即：
E[Z(x)]m E[Z(x)-m][Z(x h)-m]C(h)=⎧⎨+=⎩
其中，Z(x)为区域化变量，h 为变量间的距离，m 为常数；（2）内蕴假设（Intrinsic Hypothesis ），当区域化变量Z(x)的增量)]h x (Z )x (Z [+-满足下列两个条件：一是在整个研究区域内对任意x 和h 都有0)]h x (Z )x (Z [E =+-，二是增量)]h x (Z )x (Z [+-的方差函数存在且平稳（不依赖于x ），即：
222)]h x (Z )x (Z [E )]}h x (Z )x (Z [E {)]h x (Z )x (Z [E )]
h x (Z )x (Z [Var +-=+--+-=+-
则称Z(x)满足内蕴假设。

在地质统计学中，普通克立格法和简单克立格法一般都要求其满足上述假设条件中的一个，作为区域化变量结构分析的基础。

2.5.1.2变异函数
地质统计学的主要技术为克立格方法，它以区域化变量理论为基础，因研究对象具两重性（随机性和结构性），所以以变异函数（V ariogram ）为主要工具来分析处理数据。

所谓变异函数是指：假设空间点x 只在一维x 轴上变化，我们把区域化变量Z(x)在x ，x+h 两点处的值之差的方差之半定义为Z(x)在x 方向上的变异函数，记为)h ,x (γ，即： 22)]}h x (Z )x (Z [E {2
1)]h x (Z )x (Z [E 21)]h x (Z )x (Z [Var 21)h ,x (+--+-=+-=γ 根据2.5.1.1中两个假设条件可知，凡满足这两个条件之一的区域化变量，其变异函数的计算公式可变为：
2)]h x (Z )x (Z [E 2
1)h (+-=γ 因此，变异函数值可由区域化变量间其增量平方的数学期望来表示。

据此，又可得出实验变异函数的计算公式：
∑=*
+-=)
h (N 1i 2i i )]h x (Z )x (Z [)h (2N 1)h (γ 其中，)h (N 为样本数据的个数[1]。

如果以距离h 为横坐标，半变异函数值为纵坐标作图，即可得半变异函数图，如下所示：
h C 0
a 图5 半变异函数图
该半变异函数图中有三个重要的参数：块金效应（Nugget Effect ）、基台值（Sill ）、Range （变程或自相关阈值）。

其中，块金效应主要来源于远小于抽样间距的空间尺度上所存在的差异，它将直接影响空间内插的精度；基台值是理论上变异函数的最大值，即空间上随机变
量的先验方差；自相关阈值是数据空间相关性的界限，超过该值的两变量间无空间相关性[6]。

由于在变异值的计算中有可能出现负值，所以需要对实验变异函数进行拟合，用于此拟合的数学模型很多，就目前的研究和应用领域来看，这些模型已基本能满足我们的需求，
主要包括[1，4]：
球状模型：20300 r 03r 1r (r)C C() 0r a 2a 2a C C r a
γ=⎧⎪⎪=+⋅-⋅<≤⎨⎪+>⎪⎩，其中，C 0为块金常数，C C 0+为基台值，C 为拱高，a 为变程。

它所表述的空间相关随距离的增长逐渐衰减，当距离大于a 时，空间相关消失。

指数模型：)e 1(C C )r (a 3r
0--+=γ，各参数含义与球状模型相同。

其空间相关性随
距离的增长以指数形式衰减，相关性消失于无穷远。

高斯模型：)e C(1C )r (22
a 3r 0--+=γ，各参数含义与球状模型相同。

其空间相关性随
距离的增长而衰减，相关性消失于无穷远。

幂函数模型：(r)r 02θγλθ=<<，其中，λ为常数。

最常用的是1=θ的情况（即
线性模型）。

对于线性模型来说，其相关性随距离的增长而线性递增。

对数函数模型：)r (log )r (=γ，一般不用于点数据变异函数的拟合。

纯块金效应模型：0
0 r=0(r)C r>0γ⎧=⎨
⎩，这种模型仅适用于纯随机变量，即无空间相关性的数据。

孔穴效应模型：⎥⎦⎤⎢⎣⎡⋅-+=-)b r 2(cos e 1C C )r (a r 0πγ（一维情况下），一般当变异函数)r (γ不是单调递增，而呈周期波动时使用。

图6 含有异常值的变异函数云图，其中黄线为拟合的模型
2.5.1.3质量评定之方差公式
现实世界是复杂多变的，要想从有限的已知数据来模拟真实世界是很困难的。

使用各种建模方法难免会产生误差，研究人员在各方面的辛勤工作无非是想找到一种使误差达到最小的方法，即一种好的内插方法应是方差最小的无偏估计。

由于方差的推导很复杂，因篇幅有限，这里就不赘述了，只给出其最后的计算公式：
)v ,v ()V ,V ()v ,V (22E γγγσ--=
其中：V 为整个待估区域，
v 为某一具体的计算区域，
∑⎰='-=n 1i V
i dy )x y (V 1n 1)v ,V (γγ， ⎰⎰-=V V
2dydt )t y (V 1
)V ,V (γγ， ∑∑=='-'=n 1i n
1j j
i 2)x x (n 1)v ,v (γγ。

估计方差越小，其估计的效果越好。

从上面的公式可以看出，只要求出变异函数及其平均值，就可以求得方差。

影响方差的因素很多，一般包括：区域化变量的结构特征和空间连续性、研究区域的几何特征（大小、形状等）、样本数据的特性，即它的数量、空间分布等[1]。

克立格法是一种最优无偏的估计方法，随着其自身的不断发展和完善，在不同的研究
区域和研究尺度下，可使用不同的克立格方法来进行数据的处理和分析[6]。

下面将作具体介绍：
2.5.2普通克立格法（Ordinary Kriging ）
普通克立格法是克立格家族中最常用的一种最优无偏估计方法，其计算公式与反距离加权插值类似，具体为：
n
0i i i 1z(x )z(x )λ==∑
其中i z(x )(i=0,1,2,,n)为数据点i x 的值，0x 为未知点，其他为已知点，i λ为权。

此处的权系数不是简单的由距离来决定，它是在无偏性和最小方差性的条件下，依赖于变异函数的计算结果而确定的，此特征正是与反距离加权插值不同之处，也是其优势的体现。

普通克立格法适用的条件是：所研究的区域化变量，其数学期望未知，但为一常数。

为了满足无偏条件，需要使权系数的和为1，即1n 1i i =∑=λ。

在这个条件下，为了使估计方差
最小，就要使用求条件极值中的拉格朗日乘数法，最终得到普通克立格法方程组：
n
i i j i j 1n i i 1
(x ,x )(x ,V) (i 1,2,,n)1 λγμγλ==⎧+==⎪⎪⎨⎪=⎪⎩∑∑ 其方差公式为：
∑=+=n 1
i i i 2
K
V)(V,-V),x (μγγλσ 其中μ为拉格朗日乘数。

图7 普通克立格法内插的表面
2.5.3简单克立格法（Simple Kriging ）
当区域化变量的数学期望是已知的常数时，就要使用简单克立格法了。

平均值为已知的常数，将为克立格的研究带来很大的便利，计算起来更简单快捷。

在其计算公式n
0i i i 1z(x )z(x )λ==∑中，无需权系数的和为1，就能满足无偏条件。

其克立格方程组与普通克立格法略有不同，因简单克立格法中所研究的区域化变量满足二阶平稳假设条件，所以协方差存在；又因为克立格计算公式自动满足无偏条件，只需在满足最小方差条件下求得各权重系数值，故其方差公式为：
∑=-=n
1i i i 2
K
)V ,x (C )V ,V (C λσ 2.5.4泛克立格法（Universal Kriging ）
上面介绍的两种方法都要求区域化变量满足一个假设条件：二阶平稳或内蕴假设，但在实际情况中，有许多区域化变量的数学期望都不是固定值，而是随空间位置的变化而变化的，即)x (m )]x (Z [E =。

在地质统计学中，m(x)被叫做漂移（Drift ）。

漂移的形式可以是一维或二维情况，也可以是三维或三维以上的，一般一维和二维漂移就能满足要求。

如果用多项式来表示漂移，则：
一维： +++=2210x a x a a )x (m
二维： ++++++=25423210y a xy a x a y a x a a )y ,x (m
泛克立格方法就是用来处理存在漂移的区域化变量对象的，它有两种类型：一为已知协方差函数的泛克立格法，二为存在变异函数的泛克立格法。

由于泛克立格法计算复杂、量大，且有一定的局限性，Georges Matheron 提出了IRF-K （即
K 阶内蕴随机函数）法[7]。

该方法可以避过漂移问题，用区域化变量的增量来代
替变量本身，并允许权系数的和为0，以便将非平稳的变量过渡到平稳的变量中来。

普通克立格法、简单克立格法、泛克立格法均为线性克立格插值方法。

一般，当区域化变量满足二阶平稳条件时使用简单克立格法；满足内蕴假设时使用普通克立格法；当数据存在漂移时则用泛克立格法。

下面表中是对这三种线性克立格法的比较：
2.5.5指示克立格法（Indicator Kriging ）
指示克立格法是一种非参数方法，无需了解数据的分布类型，它是由Journel 教授于1982年提出的。

该方法因其自身的特点可以将异常值对估值的影响降到最低，因此也是人们常使用的方法之一。

其指示函数为：
⎩
⎨⎧>≤=Z )X (Z 0Z )X (Z 1)Z ;X (I 其中Z 为临界值[7]。

研究人员正是通过这个指示函数来估计未知点的指示值，即：
)L ,,2,1l ()Z ;X (I )Z ;X (I n
1
l l ==∑=*
αααλ 式中αλ由克立格方程组给出。

通过对具体位置X 处指示平均值的估计，最终可以得到)X (Z 的估计值。

2.5.6概率克立格法（Probability Kriging ）
概率克立格法是指示克立格法的一种改进，它不仅具有指示克立格法的优点，即非参数和无分布特性，同时也减小了估计方差，提高了插值精度，降低了指示克立格法的平滑作用。

因其对指示值做出了概率解释，即：
}Z )X (Z Z )X (Z {Prob )Z ;X (I αααα=≤=
所以称其为概率克立格法。

2.5.7析取克立格法（Disjunctive Kriging ）
线性地质统计学是对数据线性组合的一种研究，一般只能用来估计)X (Z 的值，具有一定局限性。

析取克立格法于1972年由Matheron 提出，是最早的一种非线性地质统计学方法，也是介于线性克立格估计量与条件数学期望之间切实可行的一种中间估计量。

其估计量的一般形式为：
∑==*n
1DK )Z (f Z ααα
其中)Z (f αα是每一个有效数据变量αZ 的函数，n ,,2,1 =α。

使用析取克立格法，不仅可以估计某一点的待估值，还可以对某点的指示值做出估计。

但它对数据的要求更严格，需数据服从二元正态分布。

同时，在计算析取克立格方程组和其方差时，须利用埃尔米特多项式（hermite polynomials ）的有限展开式，该多项式是特殊函数中的一种正交多项式，用途很广[1，3]。

尽管
一般情况下析取克立格法比普通克立格法的应用效果更佳，但有利比有弊，正由于上面所陈述的原因，析取克立格法在计算机中进行处理时困难很大，计算方法非常复杂，所以应谨慎使用。

2.5.8协同克立格内插
协同克立格法是克立格法的直接扩展，它通过研究两个或两个以上变量的自相关性和变量间的交叉相关性，获得变异函数和交叉变异函数模型，以此来优化插值估计。

对样品数据不足的情况尤其使用。

相对于克立格法来说，协同克立格法多了一个新的假定条件，即使变量间的差值所得方差最小。

该方法的核心还是对交叉变异函数的估计，其实验交叉变异函数的计算公式为：
n k
k k i i i i i 11(h)[z(x )z(x h)][z (x )z (x h)]2n γ-=-+-+∑ 式中：k z (x)是与估计值z(x)相关的第k 个变量[4]。

在ArcGIS 9软件的Geostatistics Analyst 模块中，协同克立格法有如下几个类型：普通、简单、泛、指示、概率和析取协同克立格法，您可以根据不同的情况选择使用上面的这几种协同克立格法。

4.结论
ArcGIS 9软件的Geostatistics Analyst 模块开创性地将地质统计学与地理信息系统融合在一起，实现了多种功能。

地统计学的结构分析、插值方法与GIS 中空间分析功能（如数据的显示、叠加等）相结合，结合各自的优势可以完成对GIS 空间数据的获取、存储、显示、分析、查询等一系列空间数据处理功能。

该分析模块中的内插方法基本上包含了目前一些主流的地质统计学插值方法，应用时可根据研究区域和数据的分布情况，选择最适宜的插值方法。

尽管如此，地质统计分析模块并没有涵盖所有的地质统计学内容，其中也有一些比较重要的地质统计方法没有涉及到，如
条件模拟、时空域上的地质统计学模型[9，10]、对数正态克立格法、因子克立格法等。

且该分
析模块中应用界面的可利用性不高，难以进行更深入地使用，所以还有待于人们进一步的研究与开发。