張培德,羅曉春
(1.成都理工大學(xué)信息管理學(xué)院,四川成都610059;2.221 Trinidad,Dollard-Des-Ormeaux,Quebec Canada H9G 2X3)
地質(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)的普通克立格方差可表述如下:
我們可以很容易從方程(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)實(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)格正定的。
以下準(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)格正定的。
根據(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)。
為了避免在實(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.