基于泛协克里格法的空间插值方法研究

程勖　张明会

摘要：分析区域化变量的特征，采用合理的插值方法是观察靶区样品分布及估计靶区样品储量的重要手段。该文在分析趋势模型和线性估计无偏条件的基础上，讨论了漂移向量对求解泛协克里格方程组的影响，描述了区域化向量的估计过程。以某矿区样本为例提出该矿区的泛协克里格法可视化计算流程， 不但完成了对矿体的内部属性建模，并且较准确地计算了矿体储量，误差为100KG左右。实验证明，泛协克里格法可为研究矿体的空间分布提供科学依据。

关键词：泛协克里格；无偏估计；空间插值；克里格方程组

中图分类号：TP202 文献标识码：A 文章编号：1009-3044（2015）12-0220-03

Research on Spatial Interpolation Methods of Universal Cokriging

CHENG Xu1， ZHANG Ming-hui2

（1. School of Management， Dalian Polytechnic University， Dalian 116000， China； 2. Department of Computer Science and Technology，Dalian Neusoft University of Information， Dalian 116023， China）

Abstract： Based on the analysis of regionalized variable characteristics ，adopting rational interpolation method has become one important means for observe the target sample distribution and estimation of the target zone.During analysis of trend model and linear estimate unbiased condition， discusses some drift factors that influence the effects for solving Universal Cokriging equations， describes the process of estimating regional vector. This paper puts forward the visualization computing process of the Universal Cokriging on some mining area. Not only finished on the internal attribute modeling， and accurately calculate the ore reserves， the error is about 100KG. The experiment proved， the Universal Cokriging provide a scientific basis for the study of the spatial distribution of ore bodies.

Key words： Universal Cokriging； Unbiased Estimate； Spatial Interpolation； Kriging Equations

在分析地质样本数据过程，经常遇到样本数据不足，无法显示地质空间分布的特征，为解决这一问题，常采用常规的插值方法如：反距离加权插值法、Shepard法、多项式插值法等[1-2]。但是由于样本在空间分布的杂乱性，使得插值结果易受到采样点位置、样本数量、样本数值大小的影响。克里格方法采用随机处理的思路，针对样本数据空间分布的连续性，建立线性无偏估计值得有效方法，不仅解决局部靶区空间数据的分析特性问题，而且在区域尺度中可分析样本的空间属性。但是，线性估计由于在解决连续空间数据中会出现局部异常，导致估值结构失真[3]。协克里格方法基于因子耦合分析，既考虑变化的趋势区域化向量，又充分利用了多变量间的相关性，从而保证了针对靶区数据采集较少，局部异常的问题。但是，由于数据空间分布的杂乱性，在采用协变差函数拟合空间数据的过程中均处理水平位置，限制了插值的连续性。泛协克里格法在引入泛协变差函数后，解决了空间不同位置插值的难题，拓宽了克里格的应用领域[4，5]。本文在研究泛协克里格方法的基础上，研究了该方法的理论基础及实际意义，并以某一靶区为例，深入探讨了该方法的具体实施。

1 趋势模型与线性估计的无偏条件

设有P+1个基本的趋势函数[flx，l=0，1，…，p]和单变量泛Kriging的情形一样，他们可以取为坐标向量x的多项式中的各项，特别地，恒假定[f0x≡1]。设由这些趋势函数组成的P+1维向量为[fx=f0x，f1x，…，fpx']，它使[ZX]的任意分量的数学期望都可表成基本趋势函数的线性组合，即有系数矩阵[A=ailk×p+1]，得到公式1：

[m（x）=EZ（x）=Af（x）=[f（x）'?I]A，x∈G] （1）

其中的I为k阶单位矩阵，[?]为矩阵的Kronecker乘积符号，[k（p+1）]维向量[A]是矩阵A拉直运算的结果，即若[A=a0，…，ap]，则有[A=a0a1?ap]，这里用到关于Kronecker乘积和拉直运算的一般为公式2：[A1A2A3=（A'3?A1）A2] （2）

于是得到公式3：

[m（x）=Af（x）=IAf（x）=[f（x）'?I]A] （3）

将公式（3）代入公式（1），即：

[Af（x）=l=0pfl（x）al=[f0（x）I，…，fp（x）I]A=[fx'?I]A] （4）

当k=1时，A为单行矩阵[a']，公式（1）成为：

[m（x）=EZ（x）=a'f（x）=f（x）'a] （5）

这与文献[6]中的趋势模型一致。在有趋势模型的情况下，我们还要求出[Z（x0）]的数学期望[m（x0）]的如下形式的估计公式6：

[m*（x0）=α=1nΔ'αZ（xα）=Δ'Z] （6）

其中[Δα=δijk×k]为待求的估计系数矩阵，为使其给出估计具有无偏性，应有[m（x0）=Em*（x0）=α=1nΔ'αEZ（xα）=α=1nΔ'αm（xα）]，将其代入公式（1），即：[Af（x0）=α=1nΔ'αAf（xα）]，或[[fx0'?I]A=α=1nΔ'α[fxα'?I]A]，其中I是k阶单位矩阵。由此得到无偏条件：[fx0'?I=α=1nΔ'αfxα'?I]。可以看出，此无偏条件与单变量情形UK方法的无偏条件非常相似，即：将系数[δα]换成k阶方阵[Δα]，将F和f（x0）中的每个元素 都乘以k阶单位矩阵I，就变成泛协克里格方法的无偏条件。

2 漂移向量[mx0]的估计

因[Em?（x0）=mx0]，公式6的估计方差为公式（7）：

[δ2E=trCov[m?x0-mx0]，[m?x0-mx0]'=trΔ'CovZ，Z'Δ=trΔ'CΔ] （7）

为使估计方差最小，Lagrange不定乘数的目标函数为：

[LΔ，θ=tr（Δ'CΔ）+2trθ'[（F'?I）Δ-f（x0）?I]]，其中的[θ]为由不定乘数组成的[（p+1）k×k]阶矩阵，对目标函数求导，并令导数为[O]，得估计方差极小化Krige方程组1：

[CF?IF'?IOΔθ=Ofx0?I]，我们假定方程组得（n+p+1）k阶系数矩阵为正定矩阵，由于已假定其中的子块C正定，只要f（x）中的基函数选的恰当，这条件一般能够满足。通过解方程组，我们可以求出公式6中的诸系数矩阵[Δα（α=1，…，n）]，从而求出漂移估计。

将方程组1的解代入公式7，联立上述无偏条件，即可得到最小估计方差，即泛协Krige方差公式8：

[δ2UCK=-trΔ'（F?I）θ=-tr[f（x0）'?I]θ] （8）

3 漂移系数矩阵A的估计

在讨论区域化向量的泛协克里格中，我们可以推广单变量的UK中漂移系数向量的估计方法，即求公式1中的漂移系数矩阵A的估计。为此，只需求[A]的如下形式估计：[A?=GZ]，其中[A?]和[A]都是（p+1）k维列向量，[Z]是由已知的RV组成的nk维列向量，G是待求的[（p+1）k×nk]阶估计系数矩阵，它应具有如下性质：为了保证估计的无偏性，由公式1，应有：[A=EA?=GEZ=Gm（x1）?m（xn）=Gf（x1）'?Ik?f（xn）'?IkA=G（F?Ik）A]。于是得到无偏性的充分条件：[G（F?Ik）=I（P+1）k]。

4 Krige方法的一般数学模型

各种Krige方法都可看成是线性回归分析，它们所研究的是多个RV之间的线性依赖关系[7]。按RF的漂移（数学期望）的变化情况，所研究的问题和方法可分为平稳（OK，OCK）和非平稳（UK，UCK）两种[8]。总之，根据单变量或多变量，平稳或非平稳，可以组成四种线性Krige估计方法。它们的共同目的是求未知量的线性、无偏和方差最小估计，因此每种模型都要包含线性估计式、无偏条件和估计方差三个要素，作为问题的解答都要包含Krige方程组和Krige方差[9]。泛协克里格方法的数学模型最为复杂。在二阶非平稳情形下，它的Krige方程组可以具体写成如下形式：

[C（x1-x1）…C（x1-xn）f0（x1）I…fp（x1）I…C（xn-x1）…C（xn-xn）f0（xn）I…fp（xn）If0（x1）I…f0（x1）IO…O…fp（x1）I…fp（xn）IO…OΛ1?Λnθ0?θp=C（x1-x0）?C（xn-x1）f0（x0）I?fl（x0）I]

其中的各个子块都是k阶方阵。系数矩阵是（n+p+1）k阶方阵，未知元和右端项是[（n+p+1）k×k]阶矩阵。当k=1，p>0时，各自块都退化成一个值，这方程组就成为UK的Krige方程组；当p=0，k>1时，这时OCK的Krige方程组；当p=0，且k=1时，它是OK的Krige方程组。

5 泛协克里格法应用

研究从3个方面的原始资料建立三维模型， 主要是建立矿区的钻孔三维模型 ： 一是钻孔的空间总体位置信息， 即钻孔的测量数据， 包括钻孔在三维空间的起点坐标（ X ， Y， Z） 以及钻孔的长度， 见表1； 二是钻孔在空间的位置变化信息， 即钻孔在空间的倾斜方向和倾角， 这2个关于钻孔空间位置信息的资料描述了钻孔在空间的形态， 见表2； 三是对钻孔的操作及有关的地质描述， 即采样信息，包括采样位置、样品编号、样品长度、岩性代号。最后由二维的采样信息表、钻孔位置表，经过投影变换和坐标变换生成三维钻孔立体图，如图1所示。

三维地质体的泛协克里格计算步骤如下：

1） 数据分析。主要通过散点图、频率分布图等对数据分布特征分析，挖掘特异值并对其处理。

2） 标准间距的样品组合。进行样品的归一化，针对样品值的分布特点，按指定长度对其进行加权平均，将其组合成等长的信息样品。

3） 变量A ， B 的交叉协方差建模模拟。进行实验变异函数的计算及理论变异函数的拟合，过程如下：

? 分别计算A ， B 的半变异函数。

? 如果两个变量在数量级别上有差异， 首先要对两个变量进行归一化， 使得这两个变量的数量级一致。

? 三维空间中在有A ， B 两个采样值的位置上计算出第3个变量C=A+B。

? 计算C 的半变异函数。

? 用诸如球状模型来对A， B， C 的半变异函数作模型化。

? 使用如下的公式将半变异函数模型值转化为协变异函数值。

? 计算A ， B 参数之间的交叉协相关。

4） 泛协克里格计算。通过泛协克里格方法对采集样品在整个靶区进行品位估计，最后得到空间插值的地质图，如图2所示。

5） 储量计算。根据估计品位的结果、矿体的体积，计算矿体的矿石量、金属量等数据。表3为泛协克里格计算结果与勘探数据报告结果的对比效果。

6 结论

泛协克里格方法是一种可以包含多种变量信息的插值方法，它可以同时结合较粗分辨率的空间信息和其他一些较细分辨率的空间信息进行插值估计。与其他一些插值方法相比，泛协克里格提供了一种无偏最小访查估计。本文在分析了趋势模型和线性估计的无偏条件的基础上，给出了漂移向量m（x0）的估计计算方法，在讨论区域化向量[Zx0]过程中，得到泛协克里格方程组，并且指出泛协克里格方程组如何演变为其它克里格方程组得特殊形式。以某靶区矿体为例，提出泛协克里格法的计算流程，不但完成了对矿体的内部属性建模，并且较准确地计算了矿体储量，误差为100KG左右。实验证明，泛协克里格法可为研究矿体的空间分布提供科学依据。

参考文献：

[1] 张志才，陈喜，王文，等.贵州降雨变化趋势与极值特征分析[J].地球与环境，2007，35（4）：351-356.

[2] 何亚群，左蔚然，张书敏，等.基于地质统计学的煤田煤质插值方法比较[J].煤炭学报，2008，33（5）：514-517.

[3] Lark R M， Bellamy H， Rawlins B G. Spatial-timporal variability of some metal concentrations in the soil of eastern England， and implications for soil monitoring[J]， Geodma， 2005，133（3）：363-379.

[4] 余先川，王世称，王桂安.稳健协同克立格因子分析及其在化探中的应用[J].地球科学，1997，23（2）：171-174.

[5] Lu D T， Zhang T， Yang J Q， et al. A reconstruction method of porous media integrating of data with hard data[ J] . Chinese Sci Bull.2009， 54 （11）： 1876-1885.

[6] Zhang R， Shouse P， Yates S. User of pseudo-cross variograms and cokriging to improve estimates of soil solute concentrations[J].Soil Science Society of America Jouranl，1997，61（5）：1342-1347.

[7] Zhang T， Lu D T， L i D L. Porous media reconstruct ion using across-sect ion image and multiple-point geostatistics[C]// Proceedings of 2009 International Conference on Advanced Computer Control， Singapore： IEEE Press， 22-24 Jan. ， 2009， 24-29.

[8] Zhang T， Lu D T， Li D L. A statistical information reconstruction method o f images based on multiple-point geostatistics integrating soft data with hard data[C]// Proceedings of International Symposium on Computer Science and Computational Technology.Shanghai IEEE Press， 2008，21-22： 573-578.

[9] Eulogio P I， Peter M A. M odelling the semivariograms and crosses mivariograms required in down scaling cokriging by numerical convolution-deconvolution[J] . Computers & Geo sciences， 2007，33（10）： 1273-1284.