克里金插值法的步骤是怎样的？

克里金插值法的步骤是怎样的？

正确答案：（1）输入原始数据。（2）数据检验与分析，删去明显偏离实际的采样数据点。（3）确定适当的函数，绘制半方差图，了解空间变量的集聚范围与方向。（4）进行克里金插值。

各种Krige估计方法都是多元线性回归分析的特例，所要解决的问题是：根据RFZ(u）在n个取样点uα（α＝1，…，n）处的已知观测值求某一点的估计值Z*(u）。估计值的一般表达式是

地质勘探三维可视化技术及系统开发

要求估计值满足无偏条件

地质勘探三维可视化技术及系统开发

并使估计方差 最小。

在不同的具体条件下，求满足上述要求的估计系数λα（α＝1，…，n）及估计值Z*(u）就是各种Krige方法所要解决的问题。

这组程序有一般克立格、协同克立格和指示克立格，设置不同的参数可以完成以下功能。

1.简单克立格（SK）

这里要求已知RFZ(u）的数学期EZ(u)=m(u），并不要求Z(u）是二阶平稳的。这时无偏条件（7.2.2）将（7.2.1）变成

地质勘探三维可视化技术及系统开发

极小化估计方差 得到估计系数所应满足的Krige方程组

地质勘探三维可视化技术及系统开发

其中

地质勘探三维可视化技术及系统开发

在二阶平稳的情况下，上述协方差可表为

地质勘探三维可视化技术及系统开发

最小估计方差，即Krige方差为

地质勘探三维可视化技术及系统开发

2.普通克立格（OK）

这里不要求Z(u）的数学期望已知，但却要求Z(u）是二阶平稳的，即有EZ(u)=m（常数）。这时的无偏条件（7.2.2）成为

地质勘探三维可视化技术及系统开发

从而（7.2.1）成为

地质勘探三维可视化技术及系统开发

使估计方程极小化的Krige方程组为

地质勘探三维可视化技术及系统开发

Krige方差为

地质勘探三维可视化技术及系统开发

3.泛克立格（含趋势模型的克立格.KT）

这里假定RFZ(u）的均值可以表成如下形式

地质勘探三维可视化技术及系统开发

其中ak是未知常数，fk(u）是已知的基函数，通常假设它们为u=(u1,u2,u3）的多项函数中的各项，f0(u)=1，例如

地质勘探三维可视化技术及系统开发

当K=0时，（7.2.9）式成为Z(u）的平稳条件m(u)=a0，这种方法就成为普通克立格，因此KT是OK的推广。在一般情况下，无偏条件（7.2.5）推广为

地质勘探三维可视化技术及系统开发

估计系数λβ除满足该式以外，还要满足如下方程

地质勘探三维可视化技术及系统开发

其中μk是Lagrange不定乘数，λk是Z(u）的估计值的系数。

此外，我们还可以求出均值m(u）的估计值

地质勘探三维可视化技术及系统开发

这里的估计系数 ）要满足与（7.2.12）稍有不同的方程组

地质勘探三维可视化技术及系统开发

程序可用于点和块段的3D克里金，它包括SK,OK,KT三种算法。程序中所涉及的变差函数和协方差函数的计算另有专门程序可用，关于数据搜索和从点克里金到块段克里金的过度，我们已作了专门的叙述。这里仅就克立格方法的核心内容画出流程图（图7-7、图7-8），自然就显得简单了。这里的说明也适于后面的流程图。

4.协克立格（Cokriging,CK）

为了得到Z(u）的更精确的估计，特别是当Z(u）的已知样品数不多时，可用与之相关的另外一些RF的观测值对主变量Z(u）进行估计。例如，含一个次变量Y(u）的观测值的估计式为

地质勘探三维可视化技术及系统开发

这是一种协Krige估计式，为求估计式系数λα1，λα2，所用的协方差矩阵应包括Z,Y各自的协方差函数cZ(h）,cY(h），以及Z-Y交义协方差函数CZY(h)＝Cov｛Z(u）,Y(u+h)｝和Y-Z交叉协方差函数CYZ(h）。一般地，如果协Krige估计中包含K个变量，则协方差矩阵应包含K2个协方差函数。

图7-7 计算流程（SK）

图7-8 计算流程（OK,KT)

我们的程序要求包含如下三种常用的特殊类型：

a.传统的普通协克立格。由无偏条件可知，主变量的系数和为1，其余的次变量的系数和为0，在两个变量的情形，即为

地质勘探三维可视化技术及系统开发

这就可以严格地限制次变量对估计的影响。相应的克立格方程组除（7.2.15）外，还用下式：

地质勘探三维可视化技术及系统开发

b.标准化普通协克立格。如果所用的次变量与主变量的平稳均值相等，即

地质勘探三维可视化技术及系统开发

则（7.2.14）的无偏条件成为

地质勘探三维可视化技术及系统开发

这时克立格方程组由（7.2.16）和（7.2.17）组成。

c.简单协克立格。这里假定所有变量都标准化了，使之均值为0。这时像SK一样，对权系数没有限制，求估计系数仅用（7.2.16）即可。当所涉及的变量是正态得分变换时，即可用这种方法。

除了传统的普通协克立格，所有方法只能用协方差函数，而不能用变差函数或交叉变差函数。

5.指示克立格（IK）

借用K个阈值z1…，zk，我们可将连续型的RFZ(u）转换成K个指示变量

地质勘探三维可视化技术及系统开发

相应的观测值z(u）也可以转换为指示值

地质勘探三维可视化技术及系统开发

这时，I(u,zk）的数学期望

地质勘探三维可视化技术及系统开发

地质勘探三维可视化技术及系统开发

它是RV Z(u）的分布函数（cdf）在zk处的值。

对指示RFI(u,zk）作Krige估计（zk固定）

地质勘探三维可视化技术及系统开发

图7-9 计算流程（CK)

所得结果既是观测值i(u,zk）的估计值，也是在已给观测值i(uα，zk）（α＝1，…，n，）的条件下，I(u,zk）的条件数学期望的估计值，将上述条件记为（n），这时有

地质勘探三维可视化技术及系统开发

这个估计实质上是RF Z(u）的条件分布函数（ccdf）的估计值。

对于离散的RF Z(u），也可以很容易地将它转换为指示变量。设对于 能取到，且仅能取到K个状态sk(k=1,2，…，K）中的一个。这时对可定义指示RF

地质勘探三维可视化技术及系统开发

这时就有

地质勘探三维可视化技术及系统开发

指示克立格的主要目的是对Z(u）在非取样点的不确定性做出估计。具体的估计方法可分为如下几种情形：

a.带有等均值的简单IK。对任意确定的阈值：，设指示RFI(u,z）是平稳的，即它的用（7.2.18）式给出的均值F(u,z）与u无关，即可表为 。再设 是已知的，用 代替（7.2.3）中的m(u）及m(uα），则有如下的简单IK的估计式

地质勘探三维可视化技术及系统开发

其中的权系数λα（z）满足如下的与（7.2.4）类似的方程组

地质勘探三维可视化技术及系统开发

其中c1(h,z)=Cov(I(u,z）,I(u+h,z）是相应于阈值z的指示协方差函数。如果有K个阈值，这就要求有K个指示协方差函数,K个分布函数（cdf）值F(zk），并解K个方程组。所得的估计值是Z(u）在u处的条件分布函数的估计值。

图7-10 计算流程（IK)

b.带有不等均值的简单IK。

如果均值EI(u,z)=F(u,z）对于不同的u是不等的，但却能用某种方法事先估算出来，即可认为是已知的。这时，可用不平稳的SK估计式（7.2.3）来对指示值作出估计。例如，设I(u,sk）是由离散的RF Z(u）转化而来的，它表示Z(u）是否处于状态skk,设p(u,sk）表示处于这状态的概率，则有

地质勘探三维可视化技术及系统开发

则（7.2.3）式成为

地质勘探三维可视化技术及系统开发

这里的估计系数同样满足（7.2.18）。

c.普通IK。当数据量足够充分时，可不假定均值EI(u,z)=F(u,z）是已知的。只须在可变动的邻域内用OK方法进行估计。这时的估计式与（7.2.6）类似，估计系数满足的方程组与（7.2.7）类似。

地质勘探三维可视化技术及系统开发

d.中位数（median)IK。

在一般的IK中，K个阈值zk的选取，通常要求使指示协方差函数c 1(h,zk）有明显的差异。但有那种情形：当样本指示协方差函数（或变差函数）互相成比例时，样本互相关系数就都是相似的。这时相应的连续RF Z(u）就叫马塞克模型：

地质勘探三维可视化技术及系统开发

其中ρZ和 ´分别是Z(u）的互相关系数和指示交叉互相关系数。

于是单一的互相关系数可以从Z的样本互相关系数或样本指示互相关系数较好地估算出来，这里要取一个中位数阈值zk=M，它使F(M)=0.5。这里是用一个中位数阈值代替所有的K个阈值，而且如果对所有阈值的指示数据的构形是相同的，只须解一个IK方程组，所得的权系数就可用于所有的阈值。

6.块段克立格（BK）

由于Krige算法都是线性的，我们可以对Z(u）的线性平均值进行估计。例如待估值为

地质勘探三维可视化技术及系统开发

这里是用Z(u）在V(u）内的离散平均值代替积分平均值。我们可以用离散点处的估计值Z＊ 的平均值作为平均值ZV(u）的估计值

地质勘探三维可视化技术及系统开发

如果N个点的Krige估计所用的是相同的数据，则有关的Krige方程组可以平均成一个块段Krige方程组，只须将右端项的点对点的协方差换成点对块段的协方差

地质勘探三维可视化技术及系统开发

7.泛克立格（UK）

泛克立格计算需要有网格参数、搜索参数、漂移参数和变差函数参数等参数。

a.网格参数。首先要对待估值的点位置进行确定，这些点排列成规则的网格。首点坐标是网格点中坐标值最小点的坐标。格间距是指定方向上的行距。网格数为指定方向上的行数。离散点数是为块段克立格设置的参数，它将一个单元格继续细分成更小的网格，估计小格点上的值，然后计算平均值，作为大格点所代表块段的值。

b.搜索参数。要估计一点的值，不需要把所有的样本点的信息都用上，而只用到待估点附近的样本点就可以了。

最多点数限制在指定范围内搜索到的样本点数，当搜索到的样本点多于最多点数时，只取更靠近待估点的那些样本点进行计算。

最少点数是做估计的最少样本点数，当在指定区域内的样本点数少于此值时，放弃对待估点值的估计。

图7-11 泛克立格主对话框

区最多点是为所用样本点在待估点周围分布均匀而设置的。以待估点为原点建立坐标系，每个挂限为一个区。当此参数为0时，则此参数失效。

图7-12 搜索参数设置对话框

角1、角2、角3和主半径、副半径1、副半径2一起定义了一个空间的椭球体。主轴方向为坐标系绕z轴逆转角1度，再绕旋转后的x轴旋转角2度，新y轴所指的方向。再绕新y轴逆转角3度后，新x轴所指方向为副轴1，新z轴所指方向为副轴2。椭球在主轴方向的半轴长为主半径，在副轴1方向的半轴长为副半径1，在副轴2方向的半轴长为副半径2。

若待估点在上面定义椭球的中心，则程序只用椭球内部的样本点来估计待估点的值。

c.漂移参数。程序假设漂移可以用二次多项式来表示，用户可以设定某些项的系数为0。图7-13对话框上的x、xx、xy等分别表示二次项式中x的一次项、x的平方项和x与y的乘积项。对话中未选中的项，在二次多项式中的系数为0。

图7-13 漂移参数对话框

d.变差函数参数。变差函数中子变量为0时的值叫块金效应。一个变差函数可以分解成块金效应和多个理论变差模型，每个理论模型称为一个结构。

前面介绍的变差函数模型是单方向上的变差函数，是一元的。在空间上的变差函数应该是一个三元函数，但在此对话框（图7-14）中使用另外的表达方式。

圈7-14 变差函数设置对话框

对话框中的类型号、准拱高、准变程、副变程1、副变程2和角1、角2、角3表示一个结构。准变程、副变程1、副变程2和角1、角2、角3对应于搜索参数中的主半径、副半径1、副半径2和角1、角2、角3（参看搜索参数）定义了一个椭球，准变程、副变程1、副变程2是各自对应方向的准变程，准拱高为所有方向的准拱高，类型号为理论模型序号。

---