龙源期刊网 http://www.qikan.com.cn
基于泛协克里格法的空间插值方法研究
作者：程勖 张明会
来源：《电脑知识与技术》2015年第12期
摘要：分析区域化变量的特征，采用合理的插值方法是观察靶区样品分布及估计靶区样品储量的重要手段。该文在分析趋势模型和线性估计无偏条件的基础上，讨论了漂移向量对求解泛协克里格方程组的影响，描述了区域化向量的估计过程。以某矿区样本为例提出该矿区的泛协克里格法可视化计算流程， 不但完成了对矿体的内部属性建模，并且较准确地计算了矿体储量，误差为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]。但是由于样本在空间分布的杂乱性，使得插值结果易受到采样点位置、样本数量、样本数值大小的影响。克里格方法采用随机处理的思路，针对样本数据空间分布的连续性，建立线性
龙源期刊网 http://www.qikan.com.cn
无偏估计值得有效方法，不仅解决局部靶区空间数据的分析特性问题，而且在区域尺度中可分析样本的空间属性。但是，线性估计由于在解决连续空间数据中会出现局部异常，导致估值结构失真[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）
龙源期刊网 http://www.qikan.com.cn
=α=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方程