普通克立格方法的奇異性問(wèn)題研究

張培德，羅曉春

(1.成都理工大學(xué)信息管理學(xué)院，四川成都610059;2.221 Trinidad，Dollard－Des－Ormeaux，Quebec Canada H9G 2X3)

0 前言

地質(zhì)統(tǒng)計(jì)學(xué)中的克立格技術(shù)，已經(jīng)廣泛用于各種地質(zhì)勘測(cè)和開(kāi)發(fā)中。時(shí)空域的克立格技術(shù)，是一種用于對(duì)時(shí)空域未知點(diǎn)信息進(jìn)行估計(jì)的線性插值技術(shù)，它具有以下特征:

(1)線性:

式中λi是時(shí)空域的估計(jì)點(diǎn)x(si，ti)的系數(shù)。

(2)無(wú)偏:

(3)估計(jì)方差最小:

普通克立格方法是一種十分常用的克立格方法。時(shí)空域的普通克立格方法［3］(簡(jiǎn)稱時(shí)空普通克立格)可以簡(jiǎn)單用以下線性估計(jì)量，來(lái)估計(jì)具有已知均值的時(shí)空域平穩(wěn)隨機(jī)過(guò)程RF X(s，t):

其估計(jì)方差可表示為:

運(yùn)用拉格朗日乘數(shù)法，則可得到以下時(shí)空普通克立格方程:

其相應(yīng)的矩陣方程為

其中C表示協(xié)方差陣;λ表示權(quán)系數(shù)向量;θ表示方程式(6)右邊的協(xié)方差向量。

權(quán)系數(shù)向量可由以下方程得到:

相應(yīng)的普通克立格方差可表述如下:

1 普通克立格方法的奇異性分析

1.1 協(xié)方差函數(shù)的嚴(yán)格正定要求

我們可以很容易從方程(6)得出如下結(jié)論:普通克立格方程有唯一解的充要條件，是協(xié)方差矩陣為正定陣。也就是說(shuō)，其協(xié)方差函數(shù)必須是嚴(yán)格正定的［1］。如果協(xié)方差矩陣不是正定陣，則協(xié)方差矩陣的行列式就可能為零或趨近為零，這就是所謂的普通克立格方程的奇異性問(wèn)題，其結(jié)果將導(dǎo)致普通克立格方程無(wú)解。這是一個(gè)在運(yùn)用普通克立格方法，解決實(shí)際地質(zhì)問(wèn)題的過(guò)程中需要避免的問(wèn)題。

一個(gè)時(shí)空域的協(xié)方差函數(shù)C(h，τ)被認(rèn)為是嚴(yán)格正定的，當(dāng)對(duì)所有實(shí)數(shù)值a1、a2、…、aN不全為零時(shí)，有

協(xié)方差函數(shù)的嚴(yán)格正定特征，保證了普通克立格方程的非奇異性。

為了便于討論，在此假定時(shí)空隨機(jī)過(guò)程RF在時(shí)空域是平穩(wěn)的，其時(shí)空域協(xié)方差函數(shù)可用交叉距離函數(shù)表示。我們?nèi)菀淄瞥?，在Rn×T域上的協(xié)方差函數(shù)的正定性，等同于在Rn+1域上的協(xié)方差函數(shù)的正定性［6］。

對(duì)于一個(gè)在Rn+1域上各向同性的隨機(jī)過(guò)程RF，其協(xié)方差函數(shù)C(h)可由其譜函數(shù)S(λ)表示如下:

將式(9)代入公式(8)，得到:

1.2 實(shí)例

(1)實(shí)例1。協(xié)方差函數(shù)的指數(shù)模型

在R5空間的譜函數(shù)為S(λ)，所以它在時(shí)空域上是嚴(yán)格正定的。

(2)實(shí)例2。協(xié)方差函數(shù)的高斯模型

在R5空間的譜函數(shù)為S(λ)，所以它在時(shí)空域上也是嚴(yán)格正定的。

(3)實(shí)例3。協(xié)方差函數(shù)的球狀函數(shù)模型

在R3空間的譜函數(shù)為S(λ)，所以它在R2×T時(shí)空域上是嚴(yán)格正定的。

但是，當(dāng)n＞3時(shí)球狀函數(shù)被禁用于Rn，也即是說(shuō)，球狀函數(shù)不能用于R3×T［2］。但以下模型在R5空間上是允許的［4］。

其中在R5空間的譜函數(shù)為［2］

由于J5/2是5/2階的貝葉斯函數(shù)，所以該函數(shù)在時(shí)空域上是嚴(yán)格正定的。

1.3 協(xié)方差函數(shù)嚴(yán)格正定的檢驗(yàn)準(zhǔn)則

以下準(zhǔn)則，可以用來(lái)檢驗(yàn)一個(gè)協(xié)方差函數(shù)是否是嚴(yán)格正定的。

(1)準(zhǔn)則1。一個(gè)協(xié)方差函數(shù)不是嚴(yán)格正定的，如果其譜函數(shù)是一系列德?tīng)査瘮?shù)的線性組合:

式中bk為正常數(shù);ck為常數(shù)向量，而德?tīng)査瘮?shù)為

因?yàn)閲?yán)格正定函數(shù)的譜函數(shù)，必須是零值有限，而德?tīng)査瘮?shù)的零值無(wú)限，而且德?tīng)査瘮?shù)的線性組合也是零值無(wú)限的。

(2)實(shí)例4。一個(gè)反例是cos函數(shù)C(h)=a+cos(bh)(a和b為常數(shù)，b≠0)，其在R1空間的譜函數(shù)為S(λ)，所以cos函數(shù)不是嚴(yán)格正定的。

“準(zhǔn)則1”對(duì)時(shí)空域的各向異性結(jié)構(gòu)模型的研究也很有用。在時(shí)空域R2×T的協(xié)方差函數(shù)，可以用以下基本模型表述:

式中C0為非負(fù)常量;C1、C2、…、C7分別為在R2×T、R2、R1×T、R1和T空間允許的協(xié)方差函數(shù)，其譜函數(shù)可表為

其中S0為非負(fù)常量;S1、…、S7分別是相應(yīng)于C1、…、C7的譜函數(shù)。

由于這里沒(méi)有對(duì)C1、C2、…、C7作任何限定，所以由公式(16)所表述的協(xié)方差函數(shù)在時(shí)空域R2×T上是嚴(yán)格正定的，且其子函數(shù)C1(hx，hy，τ)是嚴(yán)格正定的。

2 地質(zhì)取樣點(diǎn)的空間分布對(duì)普通克立格方法的影響

根據(jù)前面的討論可知，如果所選的協(xié)方差函數(shù)不是嚴(yán)格正定的，則在應(yīng)用普通克立格方法時(shí)，就可能遇到克立格方程奇異性無(wú)解問(wèn)題，這是我們?cè)趯?shí)際應(yīng)用中應(yīng)當(dāng)盡量避免的。

不過(guò)在一些特殊的應(yīng)用實(shí)例中，某些協(xié)方差函數(shù)雖然是非嚴(yán)格正定，但它卻對(duì)實(shí)際數(shù)據(jù)的二階統(tǒng)計(jì)特征擬合得很好。在這種情況下，人們也常常使用這些協(xié)方差函數(shù)來(lái)刻畫和解決實(shí)際問(wèn)題。例如，cos函數(shù)就時(shí)常用來(lái)刻畫具有一定周期性波動(dòng)的二階統(tǒng)計(jì)特征的實(shí)際問(wèn)題［5］。

下面，我們就以cos函數(shù)為協(xié)方差函數(shù)為例，來(lái)分析地質(zhì)取樣點(diǎn)的空間分布對(duì)普通克立格方程奇異性的影響。

假設(shè)協(xié)方差函數(shù)為

我們來(lái)分析一維空間R1的情況。如圖1所示，假設(shè)有三個(gè)地質(zhì)取樣點(diǎn)x1、x2、x3，其中x1和x2間的距離為2kπ(此處k為任意正整數(shù)):

圖1 一種取樣點(diǎn)在一維空間分布的情況Fig.1 The distribution of one kind of sample points locates in one-dimensional space

則有:

于是其協(xié)方差陣為:

容易看出，此矩陣是奇異的，因?yàn)樵摼仃嚨牡谝恍泻偷诙惺窍嗤?。而且這種現(xiàn)象可以推而廣之，即任意再加一個(gè)取樣點(diǎn)xi進(jìn)來(lái)，由于它與x1和x2間的協(xié)方差是相同的:其中i=3，4，…。

因而所得協(xié)方差矩陣的前兩行，仍然是相同的，也即該矩陣仍是奇異的。于是我們可以得到以下準(zhǔn)則:

準(zhǔn)則2。對(duì)于在一維空間R1上的協(xié)方差函數(shù)(cos(h)+1)/2，只要在地質(zhì)取樣點(diǎn)中有兩相鄰點(diǎn)之間距離是2kπ(k為任意正整數(shù))，則所得到的普通克立格方程將是奇異的。

我們還從各向異性公式(16)中得到同樣有趣的結(jié)果。假設(shè)我們把各向異性公式(16)簡(jiǎn)化為以下純各向異性結(jié)構(gòu):

其中C2(hx，hy)表示二維空間R2上的任意協(xié)方差函數(shù);C7(τ)表示時(shí)間域上的任意協(xié)方差函數(shù)。

假設(shè)有四個(gè)地質(zhì)取樣點(diǎn)，其位置形成一個(gè)如圖2所示的四邊形，則其相應(yīng)的協(xié)方差為:

圖2 在R2×T時(shí)空域上的一個(gè)“四方形”分布結(jié)構(gòu)Fig.2 A"square"distribution and structure under R2×T time-space domain

設(shè)a=C2(|x2－x1|，|y2－y1|)，b=C7(|t2－t1|)，c=C2(0，0)，d=C7(0)，則相應(yīng)的協(xié)方差矩陣為:

由于其第一行與第四行之和等于第二行與第三行之和，所以此協(xié)方差陣是奇異的。進(jìn)而推之，如果還有另一個(gè)取樣點(diǎn)，其時(shí)空域位置為(x5，y5，t5)，于是相應(yīng)的協(xié)方差為

這樣，協(xié)方差陣C就多增加了一行/列(C(h15，τ15)，C(h25，τ25)，C(h35，τ35)，C(h45，τ45)，c+d)T)。由于

所以協(xié)方差陣C仍然是奇異的。如果我們選取純各向異性公式(17)作為協(xié)方差函數(shù)，那么若是在地質(zhì)取樣點(diǎn)中有四點(diǎn)在時(shí)空域R2×T上的位置，構(gòu)成了如圖2所示的四邊形分布，則得到的普通克立格方程將是奇異的。

在實(shí)際應(yīng)用中，以上提及的取樣點(diǎn)時(shí)空域“四邊形”分布結(jié)構(gòu)是經(jīng)常遇到的。比如，在兩個(gè)污染觀測(cè)站上，在兩個(gè)相同時(shí)間上進(jìn)行污染取樣;或是在兩口鉆井中，在兩個(gè)相同時(shí)間上進(jìn)行巖芯取樣等等。

值得注意的是，在以上討論中對(duì)兩個(gè)協(xié)方差函數(shù)C2(hx，hy)和C7(τ)并未作任何限制，它們可以是在二維空間域R2上和時(shí)間域T上的任何協(xié)方差函數(shù)。例如，C2(hx，hy)可以是在R2上的指數(shù)模型，而C7(τ)可以是在T上的球狀模型，所以它們分別在R2和T上是嚴(yán)格正定的，而其線性組合在R2×T上卻是非嚴(yán)格正定的。這意味著，普通克立格方法中的純各向異性模型，具有很高的奇異性風(fēng)險(xiǎn)。

我們還可以在時(shí)空域R3×T上得到與上類似的結(jié)論。實(shí)際上，奇異性風(fēng)險(xiǎn)是隨著時(shí)空域維數(shù)的增加而增加。因?yàn)镽2×T只是R3×T的一個(gè)子空間，因而任何在R2×T上出現(xiàn)的奇異性問(wèn)題都會(huì)在R3×T上出現(xiàn)。

3 結(jié)論

為了避免在實(shí)際應(yīng)用普通克立格方法中，遇到克立格方程奇異無(wú)解的問(wèn)題，我們應(yīng)當(dāng)在選取協(xié)方差函數(shù)時(shí)盡量避免選取非嚴(yán)格正定函數(shù)。值得注意的是，一個(gè)協(xié)方差函數(shù)是否嚴(yán)格正定，與實(shí)際應(yīng)用中的空間維數(shù)密切相關(guān)。一個(gè)在高維空間上嚴(yán)格正定的函數(shù)必定在低維空間上嚴(yán)格正定，而一個(gè)在低維空間上嚴(yán)格正定的函數(shù)在高維空間上卻很可能非嚴(yán)格正定，球狀函數(shù)就是一個(gè)很好的例子。如果在實(shí)際應(yīng)用中需要選取某些非嚴(yán)格正定函數(shù)作為協(xié)方差函數(shù)，那么就需要注意取樣點(diǎn)的空間分布，盡量避免選取造成奇異性的取樣點(diǎn)空間分布結(jié)構(gòu)，如純各向異性模型中的“四邊形”結(jié)構(gòu)。綜上所述，普通克立格方法中的奇異性無(wú)解的問(wèn)題，是有因可循的，也是可以避免的。

［1］BERG C.Theory of Positive Definite and Related Functions［M］.Springer－Verlog，New York，1984.

［2］CHRISTAKOSG.On the problem of permissible covariance and variogram models［J］.Water Resources Research，1984，20(2):251.

［3］CHRISTAKOSG.Random Field Models in Earth Sciences［M］.Academic Press，New York，1992.

［4］MATERN B.Spatial Variation.Medd.Stens Skogsforskniginst［M］.Swed，1960，

［5］MELKUMYAN A，RAMOSF.A sparse covariance function for exactgaussian process inference in large datasets［M］.In The21st International Joint Conference on Artificial Intelligence，2009.

［6］LUO X.Developing Spatiotemporal Random Field Models for Earth Science and Engineering Applications［M］.Ph.D.thesis of McGill University，1998.

［7］成秋明.多重分形與地質(zhì)統(tǒng)計(jì)學(xué)方法用于勘查地球化學(xué)異?？臻g結(jié)構(gòu)和奇異性分析［J］.地球科學(xué)－中國(guó)地質(zhì)大學(xué)學(xué)報(bào)，2001(2):161.

［8］孫玉建.非參數(shù)地質(zhì)統(tǒng)計(jì)學(xué)進(jìn)展［J］.地質(zhì)與勘探，2007(4):79.

［9］高彥偉，戴經(jīng)隆，馬瑞杰.時(shí)空域克里格方法在地下水水質(zhì)評(píng)價(jià)中的應(yīng)用［J］.工程勘察，2007(10):381.

［10］譚繼強(qiáng)，丁明柱.空間數(shù)據(jù)插值方法的評(píng)價(jià)［J］.測(cè)繪與空間地理信息，2004(4):11.

［11］董輝，高光明，劉碧虹，楊自安.基于空間散亂點(diǎn)插值的曲面重構(gòu)［J］.地質(zhì)與勘探，2004(2):80.