普通克立格方法的奇異性問題研究

張培德，羅曉春

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

0 前言

地質統計學中的克立格技術，已經廣泛用于各種地質勘測和開發中。時空域的克立格技術，是一種用于對時空域未知點信息進行估計的線性插值技術，它具有以下特征:

(1)線性:

式中λi是時空域的估計點x(si，ti)的系數。

(2)無偏:

(3)估計方差最小:

普通克立格方法是一種十分常用的克立格方法。時空域的普通克立格方法［3］(簡稱時空普通克立格)可以簡單用以下線性估計量，來估計具有已知均值的時空域平穩隨機過程RF X(s，t):

其估計方差可表示為:

運用拉格朗日乘數法，則可得到以下時空普通克立格方程:

其相應的矩陣方程為

其中C表示協方差陣;λ表示權系數向量;θ表示方程式(6)右邊的協方差向量。

權系數向量可由以下方程得到:

相應的普通克立格方差可表述如下:

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

1.1 協方差函數的嚴格正定要求

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

一個時空域的協方差函數C(h，τ)被認為是嚴格正定的，當對所有實數值a1、a2、…、aN不全為零時，有

協方差函數的嚴格正定特征，保證了普通克立格方程的非奇異性。

為了便于討論，在此假定時空隨機過程RF在時空域是平穩的，其時空域協方差函數可用交叉距離函數表示。我們容易推出，在Rn×T域上的協方差函數的正定性，等同于在Rn+1域上的協方差函數的正定性［6］。

對于一個在Rn+1域上各向同性的隨機過程RF，其協方差函數C(h)可由其譜函數S(λ)表示如下:

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

1.2 實例

(1)實例1。協方差函數的指數模型

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

(2)實例2。協方差函數的高斯模型

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

(3)實例3。協方差函數的球狀函數模型

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

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

其中在R5空間的譜函數為［2］

由于J5/2是5/2階的貝葉斯函數，所以該函數在時空域上是嚴格正定的。

1.3 協方差函數嚴格正定的檢驗準則

以下準則，可以用來檢驗一個協方差函數是否是嚴格正定的。

(1)準則1。一個協方差函數不是嚴格正定的，如果其譜函數是一系列德爾塔函數的線性組合:

式中bk為正常數;ck為常數向量，而德爾塔函數為

因為嚴格正定函數的譜函數，必須是零值有限，而德爾塔函數的零值無限，而且德爾塔函數的線性組合也是零值無限的。

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

“準則1”對時空域的各向異性結構模型的研究也很有用。在時空域R2×T的協方差函數，可以用以下基本模型表述:

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

其中S0為非負常量;S1、…、S7分別是相應于C1、…、C7的譜函數。

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

2 地質取樣點的空間分布對普通克立格方法的影響

根據前面的討論可知，如果所選的協方差函數不是嚴格正定的，則在應用普通克立格方法時，就可能遇到克立格方程奇異性無解問題，這是我們在實際應用中應當盡量避免的。

不過在一些特殊的應用實例中，某些協方差函數雖然是非嚴格正定，但它卻對實際數據的二階統計特征擬合得很好。在這種情況下，人們也常常使用這些協方差函數來刻畫和解決實際問題。例如，cos函數就時常用來刻畫具有一定周期性波動的二階統計特征的實際問題［5］。

下面，我們就以cos函數為協方差函數為例，來分析地質取樣點的空間分布對普通克立格方程奇異性的影響。

假設協方差函數為

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

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

則有:

于是其協方差陣為:

容易看出，此矩陣是奇異的，因為該矩陣的第一行和第二行是相同的。而且這種現象可以推而廣之，即任意再加一個取樣點xi進來，由于它與x1和x2間的協方差是相同的:其中i=3，4，…。

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

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

我們還從各向異性公式(16)中得到同樣有趣的結果。假設我們把各向異性公式(16)簡化為以下純各向異性結構:

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

假設有四個地質取樣點，其位置形成一個如圖2所示的四邊形，則其相應的協方差為:

圖2 在R2×T時空域上的一個“四方形”分布結構Fig.2 A"square"distribution and structure under R2×T time-space domain

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

由于其第一行與第四行之和等于第二行與第三行之和，所以此協方差陣是奇異的。進而推之，如果還有另一個取樣點，其時空域位置為(x5，y5，t5)，于是相應的協方差為

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

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

在實際應用中，以上提及的取樣點時空域“四邊形”分布結構是經常遇到的。比如，在兩個污染觀測站上，在兩個相同時間上進行污染取樣;或是在兩口鉆井中，在兩個相同時間上進行巖芯取樣等等。

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

我們還可以在時空域R3×T上得到與上類似的結論。實際上，奇異性風險是隨著時空域維數的增加而增加。因為R2×T只是R3×T的一個子空間，因而任何在R2×T上出現的奇異性問題都會在R3×T上出現。

3 結論

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

［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］成秋明.多重分形與地質統計學方法用于勘查地球化學異?？臻g結構和奇異性分析［J］.地球科學－中國地質大學學報，2001(2):161.

［8］孫玉建.非參數地質統計學進展［J］.地質與勘探，2007(4):79.

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

［10］譚繼強，丁明柱.空間數據插值方法的評價［J］.測繪與空間地理信息，2004(4):11.

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