基于泛協(xié)克里格法的空間插值方法研究

程勖　張明會

摘要：分析區(qū)域化變量的特征，采用合理的插值方法是觀察靶區(qū)樣品分布及估計(jì)靶區(qū)樣品儲量的重要手段。該文在分析趨勢模型和線性估計(jì)無偏條件的基礎(chǔ)上，討論了漂移向量對求解泛協(xié)克里格方程組的影響，描述了區(qū)域化向量的估計(jì)過程。以某礦區(qū)樣本為例提出該礦區(qū)的泛協(xié)克里格法可視化計(jì)算流程， 不但完成了對礦體的內(nèi)部屬性建模，并且較準(zhǔn)確地計(jì)算了礦體儲量，誤差為100KG左右。實(shí)驗(yàn)證明，泛協(xié)克里格法可為研究礦體的空間分布提供科學(xué)依據(jù)。

關(guān)鍵詞：泛協(xié)克里格；無偏估計(jì)；空間插值；克里格方程組

中圖分類號：TP202 文獻(xiàn)標(biāo)識碼：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

在分析地質(zhì)樣本數(shù)據(jù)過程，經(jīng)常遇到樣本數(shù)據(jù)不足，無法顯示地質(zhì)空間分布的特征，為解決這一問題，常采用常規(guī)的插值方法如：反距離加權(quán)插值法、Shepard法、多項(xiàng)式插值法等[1-2]。但是由于樣本在空間分布的雜亂性，使得插值結(jié)果易受到采樣點(diǎn)位置、樣本數(shù)量、樣本數(shù)值大小的影響?？死锔穹椒ú捎秒S機(jī)處理的思路，針對樣本數(shù)據(jù)空間分布的連續(xù)性，建立線性無偏估計(jì)值得有效方法，不僅解決局部靶區(qū)空間數(shù)據(jù)的分析特性問題，而且在區(qū)域尺度中可分析樣本的空間屬性。但是，線性估計(jì)由于在解決連續(xù)空間數(shù)據(jù)中會出現(xiàn)局部異常，導(dǎo)致估值結(jié)構(gòu)失真[3]。協(xié)克里格方法基于因子耦合分析，既考慮變化的趨勢區(qū)域化向量，又充分利用了多變量間的相關(guān)性，從而保證了針對靶區(qū)數(shù)據(jù)采集較少，局部異常的問題。但是，由于數(shù)據(jù)空間分布的雜亂性，在采用協(xié)變差函數(shù)擬合空間數(shù)據(jù)的過程中均處理水平位置，限制了插值的連續(xù)性。泛協(xié)克里格法在引入泛協(xié)變差函數(shù)后，解決了空間不同位置插值的難題，拓寬了克里格的應(yīng)用領(lǐng)域[4，5]。本文在研究泛協(xié)克里格方法的基礎(chǔ)上，研究了該方法的理論基礎(chǔ)及實(shí)際意義，并以某一靶區(qū)為例，深入探討了該方法的具體實(shí)施。

1 趨勢模型與線性估計(jì)的無偏條件

設(shè)有P+1個基本的趨勢函數(shù)[flx，l=0，1，…，p]和單變量泛Kriging的情形一樣，他們可以取為坐標(biāo)向量x的多項(xiàng)式中的各項(xiàng)，特別地，恒假定[f0x≡1]。設(shè)由這些趨勢函數(shù)組成的P+1維向量為[fx=f0x，f1x，…，fpx']，它使[ZX]的任意分量的數(shù)學(xué)期望都可表成基本趨勢函數(shù)的線性組合，即有系數(shù)矩陣[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拉直運(yùn)算的結(jié)果，即若[A=a0，…，ap]，則有[A=a0a1?ap]，這里用到關(guān)于Kronecker乘積和拉直運(yùn)算的一般為公式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）

當(dāng)k=1時，A為單行矩陣[a']，公式（1）成為：

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

這與文獻(xiàn)[6]中的趨勢模型一致。在有趨勢模型的情況下，我們還要求出[Z（x0）]的數(shù)學(xué)期望[m（x0）]的如下形式的估計(jì)公式6：

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

其中[Δα=δijk×k]為待求的估計(jì)系數(shù)矩陣，為使其給出估計(jì)具有無偏性，應(yīng)有[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方法的無偏條件非常相似，即：將系數(shù)[δα]換成k階方陣[Δα]，將F和f（x0）中的每個元素 都乘以k階單位矩陣I，就變成泛協(xié)克里格方法的無偏條件。

2 漂移向量[mx0]的估計(jì)

因[Em?（x0）=mx0]，公式6的估計(jì)方差為公式（7）：

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

為使估計(jì)方差最小，Lagrange不定乘數(shù)的目標(biāo)函數(shù)為：

[LΔ，θ=tr（Δ'CΔ）+2trθ'[（F'?I）Δ-f（x0）?I]]，其中的[θ]為由不定乘數(shù)組成的[（p+1）k×k]階矩陣，對目標(biāo)函數(shù)求導(dǎo)，并令導(dǎo)數(shù)為[O]，得估計(jì)方差極小化Krige方程組1：

[CF?IF'?IOΔθ=Ofx0?I]，我們假定方程組得（n+p+1）k階系數(shù)矩陣為正定矩陣，由于已假定其中的子塊C正定，只要f（x）中的基函數(shù)選的恰當(dāng)，這條件一般能夠滿足。通過解方程組，我們可以求出公式6中的諸系數(shù)矩陣[Δα（α=1，…，n）]，從而求出漂移估計(jì)。

將方程組1的解代入公式7，聯(lián)立上述無偏條件，即可得到最小估計(jì)方差，即泛協(xié)Krige方差公式8：

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

3 漂移系數(shù)矩陣A的估計(jì)

在討論區(qū)域化向量的泛協(xié)克里格中，我們可以推廣單變量的UK中漂移系數(shù)向量的估計(jì)方法，即求公式1中的漂移系數(shù)矩陣A的估計(jì)。為此，只需求[A]的如下形式估計(jì)：[A?=GZ]，其中[A?]和[A]都是（p+1）k維列向量，[Z]是由已知的RV組成的nk維列向量，G是待求的[（p+1）k×nk]階估計(jì)系數(shù)矩陣，它應(yīng)具有如下性質(zhì)：為了保證估計(jì)的無偏性，由公式1，應(yīng)有：[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方法的一般數(shù)學(xué)模型

各種Krige方法都可看成是線性回歸分析，它們所研究的是多個RV之間的線性依賴關(guān)系[7]。按RF的漂移（數(shù)學(xué)期望）的變化情況，所研究的問題和方法可分為平穩(wěn)（OK，OCK）和非平穩(wěn)（UK，UCK）兩種[8]?？傊鶕?jù)單變量或多變量，平穩(wěn)或非平穩(wěn)，可以組成四種線性Krige估計(jì)方法。它們的共同目的是求未知量的線性、無偏和方差最小估計(jì)，因此每種模型都要包含線性估計(jì)式、無偏條件和估計(jì)方差三個要素，作為問題的解答都要包含Krige方程組和Krige方差[9]。泛協(xié)克里格方法的數(shù)學(xué)模型最為復(fù)雜。在二階非平穩(wěn)情形下，它的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階方陣。系數(shù)矩陣是（n+p+1）k階方陣，未知元和右端項(xiàng)是[（n+p+1）k×k]階矩陣。當(dāng)k=1，p>0時，各自塊都退化成一個值，這方程組就成為UK的Krige方程組；當(dāng)p=0，k>1時，這時OCK的Krige方程組；當(dāng)p=0，且k=1時，它是OK的Krige方程組。

5 泛協(xié)克里格法應(yīng)用

研究從3個方面的原始資料建立三維模型， 主要是建立礦區(qū)的鉆孔三維模型 ： 一是鉆孔的空間總體位置信息， 即鉆孔的測量數(shù)據(jù)， 包括鉆孔在三維空間的起點(diǎn)坐標(biāo)（ X ， Y， Z） 以及鉆孔的長度， 見表1； 二是鉆孔在空間的位置變化信息， 即鉆孔在空間的傾斜方向和傾角， 這2個關(guān)于鉆孔空間位置信息的資料描述了鉆孔在空間的形態(tài)， 見表2； 三是對鉆孔的操作及有關(guān)的地質(zhì)描述， 即采樣信息，包括采樣位置、樣品編號、樣品長度、巖性代號。最后由二維的采樣信息表、鉆孔位置表，經(jīng)過投影變換和坐標(biāo)變換生成三維鉆孔立體圖，如圖1所示。

三維地質(zhì)體的泛協(xié)克里格計(jì)算步驟如下：

1） 數(shù)據(jù)分析。主要通過散點(diǎn)圖、頻率分布圖等對數(shù)據(jù)分布特征分析，挖掘特異值并對其處理。

2） 標(biāo)準(zhǔn)間距的樣品組合。進(jìn)行樣品的歸一化，針對樣品值的分布特點(diǎn)，按指定長度對其進(jìn)行加權(quán)平均，將其組合成等長的信息樣品。

3） 變量A ， B 的交叉協(xié)方差建模模擬。進(jìn)行實(shí)驗(yàn)變異函數(shù)的計(jì)算及理論變異函數(shù)的擬合，過程如下：

? 分別計(jì)算A ， B 的半變異函數(shù)。

? 如果兩個變量在數(shù)量級別上有差異， 首先要對兩個變量進(jìn)行歸一化， 使得這兩個變量的數(shù)量級一致。

? 三維空間中在有A ， B 兩個采樣值的位置上計(jì)算出第3個變量C=A+B。

? 計(jì)算C 的半變異函數(shù)。

? 用諸如球狀模型來對A， B， C 的半變異函數(shù)作模型化。

? 使用如下的公式將半變異函數(shù)模型值轉(zhuǎn)化為協(xié)變異函數(shù)值。

? 計(jì)算A ， B 參數(shù)之間的交叉協(xié)相關(guān)。

4） 泛協(xié)克里格計(jì)算。通過泛協(xié)克里格方法對采集樣品在整個靶區(qū)進(jìn)行品位估計(jì)，最后得到空間插值的地質(zhì)圖，如圖2所示。

5） 儲量計(jì)算。根據(jù)估計(jì)品位的結(jié)果、礦體的體積，計(jì)算礦體的礦石量、金屬量等數(shù)據(jù)。表3為泛協(xié)克里格計(jì)算結(jié)果與勘探數(shù)據(jù)報告結(jié)果的對比效果。

6 結(jié)論

泛協(xié)克里格方法是一種可以包含多種變量信息的插值方法，它可以同時結(jié)合較粗分辨率的空間信息和其他一些較細(xì)分辨率的空間信息進(jìn)行插值估計(jì)。與其他一些插值方法相比，泛協(xié)克里格提供了一種無偏最小訪查估計(jì)。本文在分析了趨勢模型和線性估計(jì)的無偏條件的基礎(chǔ)上，給出了漂移向量m（x0）的估計(jì)計(jì)算方法，在討論區(qū)域化向量[Zx0]過程中，得到泛協(xié)克里格方程組，并且指出泛協(xié)克里格方程組如何演變?yōu)槠渌死锔穹匠探M得特殊形式。以某靶區(qū)礦體為例，提出泛協(xié)克里格法的計(jì)算流程，不但完成了對礦體的內(nèi)部屬性建模，并且較準(zhǔn)確地計(jì)算了礦體儲量，誤差為100KG左右。實(shí)驗(yàn)證明，泛協(xié)克里格法可為研究礦體的空間分布提供科學(xué)依據(jù)。

參考文獻(xiàn)：

[1] 張志才，陳喜，王文，等.貴州降雨變化趨勢與極值特征分析[J].地球與環(huán)境，2007，35（4）：351-356.

[2] 何亞群，左蔚然，張書敏，等.基于地質(zhì)統(tǒng)計(jì)學(xué)的煤田煤質(zhì)插值方法比較[J].煤炭學(xué)報，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] 余先川，王世稱，王桂安.穩(wěn)健協(xié)同克立格因子分析及其在化探中的應(yīng)用[J].地球科學(xué)，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.