普通克立格方法的奇异性问题研究

张培德，罗晓春

(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］成秋明.多重分形与地质统计学方法用于勘查地球化学异常空间结构和奇异性分析［J］.地球科学－中国地质大学学报，2001(2):161.

［8］孙玉建.非参数地质统计学进展［J］.地质与勘探，2007(4):79.

［9］高彦伟，戴经隆，马瑞杰.时空域克里格方法在地下水水质评价中的应用［J］.工程勘察，2007(10):381.

［10］谭继强，丁明柱.空间数据插值方法的评价［J］.测绘与空间地理信息，2004(4):11.

［11］董辉，高光明，刘碧虹，杨自安.基于空间散乱点插值的曲面重构［J］.地质与勘探，2004(2):80.