Kriging方法的公式推导

地统计（Geostatistics）又称地质统计，它是以区域化变量为基础，借助变异函数，研究既具有随机性又具有结构性，或空间相关性和依赖性的自然现象的一门科学。凡是与空间数据的结构性和随机性，或空间相关性和依赖性，或空间格局与变异有关的研究，并对这些数据进行最优无偏内插估计，或模拟这些数据的离散性、波动性时，皆可应用地统计学的理论与方法。地统计方法的条件：(1)随机过程(2)正态分布（3)平稳性：①均值平稳②协方差函数有关的二阶平稳和与半变异函数有关的内蕴平稳Kriging插值的条件：变异函数和相关分析的结果表明区域化变量存在空间相关性。首先假设区域化变量满足二阶平稳假设和本征假设，其数学期望为m，协方差函数及变异函数存在。即假设在待估计点（x）的临域内共有n个实测点，即x1，x2，…，xn，其样本值为。那么，普通克里格法的插值公式为)(xZc(h))(hmxZE)]([2)]()([)(mhxZxZEhc2)]()([21)(hxZxZEhniiixZxZ1*)()()(ixZ(4.2.22)其中为权重系数，表示各空间样本点处的观测值对估计值的贡献程度。可见，克立格插值的关键就是计算权重系数。显然，权重系数的求取必须满足两个条件：一是使的估计是无偏的，即偏差的数学期望为零；二是最优的，即使估计值和实际值之差的平方和最小。为此，需要满足以下两个条件：iix)(ixZ)(*xZi)(*xZ)(*xZ)(ixZ(1)无偏性。要使成为的无偏估计量，即当时，也就是当时，则有这时，为的无偏估计量。（2）最优性。在满足无偏性条件下，估计方差为mxZExZEniiiniii11)]([)]([)]([)]([*xZExZE)(*xZ)(ixZmxZE)]([nii11)(*xZ)(ixZ212*2])()([)]()([niiiExZxZExZxZE使用协方差函数表达，它可以进一步写为(4.2.24)为使估计方差最小，根据拉格朗日乘数原理，令(4.2.25)求F对和的偏导数，并令其为0，得克立格方程组(4.2.26)ninjniiijijiExxcxxcxxc1112),(2),(),()1(212niiEFiniiijinjjiFxxcxxcF110)1(202),(2),(2niiiExxcxxcσ12),(),(niiijinjjxxcxxc111),(),((4.2.27)(4.2.28)整理后得解线性方程组(4.2.27)式,求出权重系数λi和拉格朗日系数μ,代入公式(4.2.24),可得克立格估计方差上述过程也可用矩阵形式表示，令则普通克立格方程组为其估计方差为1),(),(),(,,01111112121212222111211xxcxxcxxcDcccccccccKnnnnnnnnDKDK1DxxcTK),(2在变异函数存在的条件下，根据协方差与变异函数的关系：，也可以用变异函数表示普通克立格方程组和克立格估计方差，即(4.2.29)解线性方程组(4.2.27)式，求出权重系数和拉格朗日乘数μ，代入公式(4.2.24)，可得克立格估计方差，即)()0()(hcchniiijinjjxxxx111),(),(),(),(12xxxxniiiK(4.2.30)也可以将克立格方程组和估计方差用变异函数写成上述矩阵形式。令在以上的介绍中，区域化变量的数学期望可以是已知或未知的。如果m是已知常数，称为简单克立格法；如果m是未知常数，称为普通克立格法。不管是哪一种方法，均可根据方法计算权重系数和克立格估计量。1),(),(),(,,01111112121212222111211xxxxxxDKnnnnnnnnDKDK1),(2xxDTK)(xZmxZE)]([(4.2.34)(4.2.35)(4.2.36)克里格插值基础克里格方法概述克里格方法（Kriging）又称空间局部插值法，是以变异函数理论和结构分析为基础，在有限区域内对区域化变量进行无偏最优估计的一种方法，是地统计学的主要内容之一。其实质是利用区域化变量的原始数据和变异函数的结构特点，对未知样点进行线性无偏、最优估计。无偏是指偏差的数学期望为0，最优是指估计值与实际值之差的平方和最小。克里格插值基础2.克里格方法的具体步骤导入数据数据分析是否服从正态分布是否是否存在趋势否是数据变换泛克里格方法根据数据选择合适的方法计算样点间的距离矩阵计算样点间的属性方差按距离分组按组统计平均距离及对应的平均方差绘制方差变异云图绘制经验半变异函数图拟合理论半变异函数图计算克里格系数进行预测图10.26克里格方法流程图例如:某地区降水量是一个区域化变量，其变异函数的实测值及距离h的关系见下表，下面我们试用回归分析方法建立其球状变异函数模型。)(h实测值γ(h)距离h实测值γ(h)距离h2.10.69.24.94.31.110.35.15.72.210.56.26.52.510.97.57.83.111.29.58.83.812.49.8从上面的介绍和讨论，我们知道，球状变异函数的一般形式为当时，有ahccahahahcchh03300)223(00)(ah0330)2()23()(hachacch如果记，则可以得到线性模型根据表中的数据，对上式进行最小二乘拟合它是一种数学上的近似和优化，利用已知的数据得出一条直线或者曲线，使之在坐标系上与已知数据之间的距离的平方和最小。，得到(4.2.20)32132100,,21,23,),(hxhxacbacbcbhy22110xbxbby2192007.0731.1048.2xxy(4.2.19)比较(4.2.20)式与(4.2.19)式，并做简单计算可知：c0=2.048，c=1.154，a=8.353，所以，球状变异函数模型为535.8202.3535.80)535.821535.823(154.1048.200)(33*hhhhhh(4.2.21)4个观测点x1，x2，x3，x4的观测值分别为Z(x1)=37、Z(x2)=42、Z(x3)=36、Z(x4)=35，如果假设降水量的变异函数是向同性（即变异函数在各个方向的变化都相同）的二维球状模型，其具体形式为（4.2.21）式。现在，我们用普通克立格法估计观测点x0的降水量值Z(x0)。当时，根据克立格矩阵的对称性，当时，，由此计算可得ji202.3154.1048.2)0(044332211cccccccjijijiijxxxxcc202.3)(870.0)]535.8)2(21535.8223(154.1048.2[202.3)2(202.3)11(202.33322042112ccc542.0)13(202.3223113cc711.0)12(202.322024114ccc601.0)22(202.3223223cc383.0)14(202.3224334cc466.0)23(202.3224224cc将以上计算结果代入克立格方程组（4.2.31），得952.0)1(202.3201c571.0)3(202.3203c473.0301.0202.0210.0287.01870.0571.0711.0952.0011111202.3383.0466.0711.01383.0202.3601.0542.01466.0601.0202.3870.01711.0542.0870.0202.314321即克立格权重系数分别为：λ1=0.287，λ2=0.210,λ3=0.202,λ4=0.301，μ=-0.473，所以观测点的降水量的克立格估计值为：根据普通克立格法的基本原理，我们知道，Z(x0)估计的基本公式应该是37.25（mm）。)(473.0)(202.0)(210.0)(287.04321*0xZxZxZxZZ35473.036202.042210.037287.0