Kriging方法的公式推导
克里金(kriging)插值的原理与公式推导
克里金(kriging)插值的原理与公式推导
克里金插值是一种空间插值方法，用于估计未知区域的数值，其
原理是基于空间数据的空间相关性来进行插值。
具体来说，克里金插
值假设空间数据在不同位置之间具有一定的相关性，即在空间上相邻
的点具有相似的数值。
克里金插值利用这种相关性来进行插值，从而
可以更准确地估计未知位置的数值。
克里金插值的公式推导涉及到半变异函数的定义，通常使用高斯
模型、指数模型或球形模型来描述数据的空间相关性。
在推导过程中，会利用已知数据点的数值和位置信息，以及半变异函数的参数来构建
插值模型，进而估计未知位置的数值。
克里金插值的公式可以表示为：
\[Z(u) = \sum_{i=1}^{n} \lambda_i \cdot Z(u_i)\]
其中，\(Z(u)\)为未知位置的数值，\(Z(u_i)\)为已知数据点的
数值，\(\lambda_i\)为插值权重，通过半变异函数及数据点之间的空
间距离计算得出。
除了基本的克里金插值方法外，还有一些相关的扩展方法，如普通克里金、泛克里金等，这些方法在建模和插值的过程中考虑了更多的因素，如均值趋势、空间方向等，使得插值结果更加准确和可靠。
总的来说，克里金插值是一种常用的空间插值方法，适用于各种地学环境下的数据分析与建模。
在实际应用中，需要根据具体数据的特点选择合适的插值方法和模型参数，以获得准确的插值结果。
克里金插值-Kriging插值-空间统计-空间分析
克里金插值方法-Kriging 插值-空间统计-空间分析1.1 Kriging 插值克里金插值（Kriging 插值）又称为地统计学，是以空间自相关为前提，以区域化变量理论为基础，以变异函数为主要工具的一种空间插值方法。
克里金插值的实质是利用区域化变量的原始数据和变异函数的结构特点，对未采样点的区域化变量的取值进行线性无偏、最优估计。
克里金插值包括普通克里金插值、泛克里金插值、指示克里金插值、简单克里金插值、协同克里金插值等，其中普通克里金插值是最为常用的克里金插值方法。
以下介绍普通克里金插值的原理。
包括普通克里金方法在内的各种克里金插值方法的使用前提是空间数据存在着显著的空间相关性。
判断数据空间相关性是否显著的工具是半变异函数（semi-variogram ），该函数以任意两个样本点之间的距离h 为自变量，在h 给定的条件下，其函数值估计方法如下：2||||1()[()()]2()i j i j s s h h z s z s N h γ-==-∑其中()N h 是距离为h 的样本点对的个数。
()h γ最大值与最小值的差m a x m i n γγ-可以度量空间相关性的强度。
max min γγ-越大，空间相关性越强。
如果()h γ是常数，即max min 0γγ-=，则说明无论样本点之间的距离是多少，样本点之间的差异不变，也就是说样本点上的值与其周围样本点的值无关。
在实际操作中，会取一些离散的h 值，当||s s ||i j -接近某个h 时，即视为||||i j s s h -=。
然后会通过这些离散点拟合成连续的半变异函数。
拟合函数的形式有球状、指数、高斯等。
在数据存在显著的空间相关性的前提下，可以采用普通克里金方法估计未知点上的值。
普通克里金方法的基本公式如下：01ˆ()()()n i ii Z s w s Z s ==∑普通克里金方法的基本思想是：通过调整i s 的权重()i w s ，使未知点的估计值0ˆ()Z s 满足两个要求：1.0ˆ()Z s 是无偏估计，即估计误差的期望值为0，2.估计误差的方差达到最小。
克里金法
Z ( x ) i Z ( x i )
* i 1
n
i 为权重系数，表示各空间样本点处的观测值对估值的影响度或者贡
献程度。 显然，克里格估值的关键问题就是在于求解 i 的值，同时根据估值 的基本原则，即无偏性和估计方差最小（最优性）的要求，具体就是要满 足以下条件：
整理后得：
n j c( xi , x j ) c( xi , x) j 1 n 1 i i=1，2，3…… i 1
解上式线性方程组,求出权重系数λi和拉格朗日系数μ,代入公式
2 E c( x, x) i j c( xi , x j ) 2 i c( xi , x) i 1 j 1 i 1 n n n
可得克里格估计方差
2 σ E c( x, x) i c( xi , x) i 1 n
上述过程也可用矩阵形式表示，令
c11 c12 c 21 c22 K cn1 cn 2 1 1
c1n c2 n cnn 1
1 1 , 1 0
1 2 , n
c( x1 , x) c ( x , x ) 2 D c ( xn , x ) 1
首先，假设区域变化变量为Z(x)，其满足内蕴假设条件和 二阶平稳条件，数学期望为m，协方差函数c(h)及变异函数 (h) 存在，即：
E[ Z ( x)] m c(h) E[ Z ( x) Z ( x h)] m 2 1 (h) E[ Z ( x) Z ( x h)]2 2
克里金插值法
克里金插值法克里金插值法又称空间局部插值法，是以变异函数理论和结构分析为基础，在有限区域内对区域化变量进行无偏最优估计的一种方法，是地统计学的主要内容之一，由南非矿产工程师D. Matheron 于1951年在寻找金矿时首次提出，法国著名统计学家G. Matheron 随后将该方法理论化、系统化，并命名为Kriging ，即克里金插值法。
1 克里金插值法原理克里金插值法的适用范围为区域化变量存在空间相关性，即如果变异函数和结构分析的结果表明区域化变量存在空间相关性，则可以利用克里金插值法进行内插或外推。
其实质是利用区域化变量的原始数据和变异函数的结构特点，对未知样点进行线性无偏、最优估计，无偏是指偏差的数学期望为0，最优是指估计值与实际值之差的平方和最小[1]。
因此，克里金插值法是根据未知样点有限领域内的若干已知样本点数据，在考虑了样本点的形状、大小和空间方位，与未知样点的相互空间关系，以及变异函数提供的结构信息之后，对未知样点进行的一种线性无偏最优估计。
假设研究区域a 上研究变量Z （x ），在点x i ∈A （i=1，2，……，n ）处属性值为Z （x i ），则待插点x 0∈A 处的属性值Z （x 0）的克里金插值结果Z*（x 0）是已知采样点属性值Z （x i ）（i=1，2，……，n ）的加权和，即：)()(10*i ni i x Z x Z ∑==λ （1） 式中i λ是待定权重系数。
其中Z(x i )之间存在一定的相关关系，这种相关性除与距离有关外，还与其相对方向变化有关，克里金插值方法将研究的对象称“区域化变量”针对克里金方法无偏、最小方差条件可得到无偏条件可得待定权系数i λ (i=1，2，……，n)满足关系式：11=∑=n i i λ（2）以无偏为前提，kriging 方差为最小可得到求解待定权系数i λ的方程组：⎪⎪⎩⎪⎪⎨⎧=⋯⋯==+∑∑==1)n ,2,1)(,(),(101n i i j j i n i i j x x C x x C λμλ， （3） 式中，C （x i ，x j ）是Z(x i )和Z(x j )的协方差函数。
python普通克里金（Kriging）法的实现
python普通克⾥⾦（Kriging）法的实现克⾥⾦法时⼀种⽤于空间插值的地学统计⽅法。
克⾥⾦法⽤半变异测定空间要素，要素即⾃相关要素。
半变异公式为：其中γ(h) 是已知点 xi 和 xj 的半变异，***h***表⽰这两个点之间的距离，z是属性值。
假设不存在漂移，普通克⾥⾦法重点考虑空间相关因素，并⽤拟合的半变异直接进⾏插值。
估算某测量点z值的通⽤⽅程为：式中，z0是待估计值，zx是已知点x的值，Wx是每个已知点关联的权重，s是⽤于估计的已知点数⽬。
权重可以由⼀组矩阵⽅程得到。
此程序对半变异进⾏拟合时采⽤的时最简单的正⽐例函数拟合数据为csv格式保存格式如下：第⼀⾏为第⼀个点以此类推最后⼀⾏是待求点坐标，其中z为未知值，暂且假设为0代码如下：import numpy as npfrom math import*from numpy.linalg import *h_data=np.loadtxt(open('⾼程点数据.csv'),delimiter=",",skiprows=0) print('原始数据如下(x,y,z):\n未知点⾼程初值设为0\n',h_data)def dis(p1,p2):a=pow((pow((p1[0]-p2[0]),2)+pow((p1[1]-p2[1]),2)),0.5)return adef rh(z1,z2):r=1/2*pow((z1[2]-z2[2]),2)return rdef proportional(x,y):xx,xy=0,0for i in range(len(x)):xx+=pow(x[i],2)xy+=x[i]*y[i]k=xy/xxreturn kr=[];pp=[];p=[];for i in range(len(h_data)):pp.append(h_data[i])for i in range(len(pp)):for j in range(len(pp)):p.append(dis(pp[i],pp[j]))r.append(rh(pp[i],pp[j]))r=np.array(r).reshape(len(h_data),len(h_data))r=np.delete(r,len(h_data)-1,axis =0)r=np.delete(r,len(h_data)-1,axis =1)h=np.array(p).reshape(len(h_data),len(h_data))h=np.delete(h,len(h_data)-1,axis =0)oh=h[:,len(h_data)-1]h=np.delete(h,len(h_data)-1,axis =1)hh=np.triu(h,0)rr=np.triu(r,0)r0=[];h0=[];for i in range(len(h_data)-1):for j in range(len(h_data)-1):if hh[i][j] !=0:a=h[i][j]h0.append(a)if rr[i][j] !=0:a=rr[i][j]r0.append(a)k=proportional(h0,r0)hnew=h*ka2=np.ones((1,len(h_data)-1))a1=np.ones((len(h_data)-1,1))a1=np.r_[a1,[[0]]]hnew=np.r_[hnew,a2]hnew=np.c_[hnew,a1]print('半⽅差联⽴矩阵：\n',hnew)oh=np.array(k*oh)oh=np.r_[oh,[1]]w=np.dot(inv(hnew),oh)print('权阵运算结果：\n',w)z0,s2=0,0for i in range(len(h_data)-1):z0=w[i]*h_data[i][2]+z0s2=w[i]*oh[i]+s2s2=s2+w[len(h_data)-1]print('未知点⾼程值为：\n',z0)print('半变异值为：\n',pow(s2,0.5))input()运算结果python初学，为了完成作业写了个⼩程序来帮助计算，因为初学知识有限，有很多地⽅写的很复杂，可以优化的地⽅很多。
kriging(克里金方法,克里金插值)
(h) C(0) C(h)
（二阶平稳假设条件下变差函数与协方差的关系）
变程(Range) ：指区域化变量在空间上具有相关性的 范围。在变程范围之内，数据具有相关性；而在变 程之外，数据之间互不相关，即在变程以外的观测 值不对估计结果产生影响。
具不同变程 的克里金插 值图象
块金值(Nugget) ：变差函数如果在原点间断，在地质统计学中称 为“块金效应”，表现为在很短的距离内有较大的空间变异性， 无论h多小，两个随机变量都不相关 。它可以由测量误差引起， 也可以来自矿化现象的微观变异性。在数学上，块金值c0相当于 变量纯随机性的部分。
min
应用拉格朗日乘数法求条件极值
j
E
Z *x0 Zx0 2
2
n
j
0,
i1
j 1,, n
Z*(x0)
进一步推导，可得到n+1阶的线性方程组, 即克里金方程组
n
i 1
C
xi
xj
i
C
x0
n
xj
i 1
i 1
j 1,, n
当随机函数不满足二阶平稳,而满足内蕴（本征）假设时, 可用变差函数来表示克里金方程组如下:
①在整个研究区内有 E[Z(u)-Z(u+h)] = 0
可出现E[Z(u)]不存在， 但E[Z(u)-Z(u+h)]存在并为零的情况
E[Z(u)]可以变化,但E[Z(u)-Z(u+h)]=0
② 增量[Z(u)-Z(u+h)]的方差函数 (变差函数,Variogram)
存在且平稳 (即不依赖于u)，即：
Var[Z(u)-Z(u+h)] = E[Z(u)-Z(u+h)]2-{E[Z(u)-Z(u+h)]}2 = E[Z(u)-Z(u+h)]2 = 2γ(u,h) = 2γ(h),
克里金插值(kriging)
二、统计推断与平稳要求
任何统计推断（cdf，数学期望等）均要求重复取样。 但在储层预测中，一个位置只能有一个样品。 同一位置重复取样，得到cdf，不现实
P
考虑邻近点，推断待估点
区域化变量： 能用其空间分布来表征一个自然现象的变量。
（将空间位置作为随机函数的自变量）
空间一点处的观测值可解释为一个随机变量在该点
P
F(u; z) F(u h; z)
可从研究区内所有数据的累积直方图推断而得 （将邻近点当成重复取样点）
太强的假设，不符合实际
二阶平稳
当区域化变量Z(u)满足下列二个条件时，则称其 为二阶平稳或弱平稳：
① 在整个研究区内有Z(u)的数学期望存在， 且等于常数，即： E[Z(u)] = E[Z(u+h)] = m(常数) x h
为相应的观测值。区域化变量在 x0处的值 z* x0 可
采用一个线性组合来估计：
n
z*x0 i zxi i 1
无偏性和估计方差最小被作为 i 选取的标准
无偏 E Zx0 Z * x0 0 最优 Var Zx0 Z * x0 min
绝对收敛，则称它为ξ的数学期望，记为E(ξ)。
E(ξ) =
xp( x)dx
数学期望是随机变量的最基本的数字特征，
相当于随机变量以其取值概率为权的加权平均数。
从矩的角度说，数学期望是ξ的一阶原点矩。
对于一组样本：
N
( zi )
m i1 N
(2)方差 为随机变量ξ的离散性特征数。若数学期望
随机函数在空间上的变化没有明显趋势， 围绕m值上下波动。
克里金插值法
克里金插值法及其合用范围克里金插值法又称空间局部插值法，是以变异函数理论和构造剖析为基础，在有限地区内对地区化变量进行无偏最优预计的一种方法，是地统计学的主要内容之一，由南非矿产工程师 D.Matheron于1951年在找寻金矿时首次提出，法国着名统计学家G.Matheron随后将该方法理论化、系统化，并命名为Kriging，即克里金插值法。
克里金插值法原理克里金插值法的合用范围为地区化变量存在空间有关性，即假如变异函数和构造剖析的结果表示地区化变量存在空间有关性，则能够利用克里金插值法进行内插或外推。
其本质是利用地区化变量的原始数据和变异函数的构造特色，对未知样点进行线性无偏、最优预计，无偏是指偏差的数学希望为0，最优是指预计值与本质值之差的平方和最小[1]。
所以，克里金插值法是依据未知样点有限领域内的若干已知样本点数据，在考虑了样本点的形状、大小和空间方向，与未知样点的互相空间关系，以及变异函数供给的构造信息以后，对未知样点进行的一种线性无偏最优预计。
假定研究地区a上研究变量Z（x），在点xi A（i=1，2，，n）处属性值为Z（xi），则待插点x0 A处的属性值Z（x0）的克里金插值结果 Z*（x0）是已知采样点属性值Z（xi）（i=1，2，，n）的加权和，即：Z*(x0)ni Z(x i)（1）i1式中i是待定权重系数。
此中Z(xi)之间存在必定的有关关系，这类有关性除与距离有关外，还与其相对方向变化有关，克里金插值方法将研究的对象称“地区化变量”针对克里金方法无偏、最小方差条件可获取无偏条件可得待定权系数i (i=1，2，，n)知足关系式：ni1（2）i1以无偏为前提，kriging方差为最小可获取求解待定权系数i的方程组：ni C(x i,x j) C(x0,x j)(j 1,2, ，n)i 1ni 1i1（3）式中，C（xi，xj）是Z(xi)和Z(xj)的协方差函数。
2国内外研究进展从克里金方法被提出到此刻已有完美的理论，并在好多领域获取了本质的应用，在某些领域的应用又推进了克里金理论的发展[3]。
克里金插值方法介绍 武汉大学 高等水文学
(h) C(0) C(h)
（二阶平稳假设条件下边查函数与写防查的关系）
变程(Range) ：指区域化变量在空间上具有相关性的 范围。在变程范围之内，数据具有相关性；而在变 程之外，数据之间互不相关，即在变程以外的观测 值不对估计结果产生影响。
具不同变程 的克里金插 值图象
块金值(Nugget) ：变差函数如果在原点间断，在地质统计学中称 为“块金效应”，表现为在很短的距离内有较大的空间变异性， 无论h多小，两个随机变量都不相关 。它可以由测量误差引起， 也可以来自矿化现象的微观变异性。在数学上，块金值c0相当于 变量纯随机性的部分。
E[ξ-E(ξ)]2存在，则称它为ξ的方差，记为D(ξ)， 或Var(ξ)，或σξ2。
D(ξ)= E[ξ-E(ξ)]2 其简算公式为
D(ξ)=E(ξ2) –[E(ξ)]2
方差的平方根为标准差，记为σξ
σξ=
D( ) E[ - E( )]2 E( 2) -[E( )]2
从矩的角度说，方差是ξ的二阶中心矩。
相当于要求：Z(u)的变差函数存在且平稳。
可出现协方差函数不存在，但变差函数存在的情况。
例：物理学上的著名的布朗运动是一种呈现出无限 离散性的物理现象，其随机函数的理论模型就是维 纳-勒维(Wiener-Levy)过程(或随机游走过程)。
布朗运动:
既不能确定验前方差，也不能确定协方差函数。
但是其增量却具有有限的方差： Var[Z(x)-Z(x+h)] = 2 (h)= A·|h| (其中，A是个常数)，
如具有三个自变量(空间
点的三个直角坐标)的随
机场
随机函数的特征值
协方差(Covariance)： 二个随机变量ξ，η的协方差为二维随机变量(ξ，
克里格法Kriging——有公式版
克里格法（Kriging）——有公式版二、克里格法（Kriging）克里格法（Kriging）是地统计学的主要内容之一，从统计意义上说，是从变量相关性和变异性出发，在有限区域内对区域化变量的取值进行无偏、最优估计的一种方法；从插值角度讲是对空间分布的数据求线性最优、无偏内插估计一种方法。
克里格法的适用条件是区域化变量存在空间相关性。
克里格法，基本包括普通克里格方法（对点估计的点克里格法和对块估计的块段克里格法）、泛克里格法、协同克里格法、对数正态克里格法、指示克里格法、折取克里格法等等。
随着克里格法与其它学科的渗透，形成了一些边缘学科，发展了一些新的克里金方法。
如与分形的结合，发展了分形克里金法；与三角函数的结合，发展了三角克里金法；与模糊理论的结合，发展了模糊克里金法等等。
应用克里格法首先要明确三个重要的概念。
一是区域化变量；二是协方差函数，三是变异函数一、区域化变量当一个变量呈空间分布时，就称之为区域化变量。
这种变量反映了空间某种属性的分布特征。
矿产、地质、海洋、土壤、气象、水文、生态、温度、浓度等领域都具有某种空间属性。
区域化变量具有双重性，在观测前区域化变量Z(X)是一个随机场，观测后是一个确定的空间点函数值。
区域化变量具有两个重要的特征。
一是区域化变量Z(X)是一个随机函数，它具有局部的、随机的、异常的特征；其次是区域化变量具有一般的或平均的结构性质，即变量在点X 与偏离空间距离为h的点X＋h处的随机量Z(X)与Z(X+h)具有某种程度的自相关，而且这种自相关性依赖于两点间的距离h与变量特征。
在某种意义上说这就是区域化变量的结构性特征。
二、协方差函数协方差又称半方差，是用来描述区域化随机变量之间的差异的参数。
在概率理论中，随机向量X与Y的协方差被定义为：区域化变量在空间点x 和x+h处的两个随机变量Z(x) 和Z(x+h) 的二阶混合中心矩定义为Z(x) 的自协方差函数，即区域化变量Z(x) 的自协方差函数也简称为协方差函数。
克里格法Kriging——有公式版
克里格法（Kriging）——有公式版二、克里格法（Kriging）克里格法（Kriging）是地统计学的主要内容之一，从统计意义上说，是从变量相关性和变异性出发，在有限区域内对区域化变量的取值进行无偏、最优估计的一种方法；从插值角度讲是对空间分布的数据求线性最优、无偏内插估计一种方法。
克里格法的适用条件是区域化变量存在空间相关性。
克里格法，基本包括普通克里格方法（对点估计的点克里格法和对块估计的块段克里格法）、泛克里格法、协同克里格法、对数正态克里格法、指示克里格法、折取克里格法等等。
随着克里格法与其它学科的渗透，形成了一些边缘学科，发展了一些新的克里金方法。
如与分形的结合，发展了分形克里金法；与三角函数的结合，发展了三角克里金法；与模糊理论的结合，发展了模糊克里金法等等。
应用克里格法首先要明确三个重要的概念。
一是区域化变量；二是协方差函数，三是变异函数一、区域化变量当一个变量呈空间分布时，就称之为区域化变量。
这种变量反映了空间某种属性的分布特征。
矿产、地质、海洋、土壤、气象、水文、生态、温度、浓度等领域都具有某种空间属性。
区域化变量具有双重性，在观测前区域化变量Z(X)是一个随机场，观测后是一个确定的空间点函数值。
区域化变量具有两个重要的特征。
一是区域化变量Z(X)是一个随机函数，它具有局部的、随机的、异常的特征；其次是区域化变量具有一般的或平均的结构性质，即变量在点X 与偏离空间距离为h的点X＋h处的随机量Z(X)与Z(X+h)具有某种程度的自相关，而且这种自相关性依赖于两点间的距离h与变量特征。
在某种意义上说这就是区域化变量的结构性特征。
二、协方差函数协方差又称半方差，是用来描述区域化随机变量之间的差异的参数。
在概率理论中，随机向量X与Y的协方差被定义为：区域化变量在空间点x 和x+h处的两个随机变量Z(x) 和Z(x+h) 的二阶混合中心矩定义为Z(x) 的自协方差函数，即区域化变量Z(x) 的自协方差函数也简称为协方差函数。
克里金插值法
克里金插值法克里金插值法又称空间局部插值法，是以变异函数理论和结构分析为基础，在有限区域内对区域化变量进行无偏最优估计的一种方法，是地统计学的主要内容之一，由南非矿产工程师D. Matheron 于1951年在寻找金矿时首次提出，法国著名统计学家G. Matheron 随后将该方法理论化、系统化，并命名为Kriging ，即克里金插值法。
1 克里金插值法原理克里金插值法的适用范围为区域化变量存在空间相关性，即如果变异函数和结构分析的结果表明区域化变量存在空间相关性，则可以利用克里金插值法进行内插或外推。
其实质是利用区域化变量的原始数据和变异函数的结构特点，对未知样点进行线性无偏、最优估计，无偏是指偏差的数学期望为0，最优是指估计值与实际值之差的平方和最小[1]。
因此，克里金插值法是根据未知样点有限领域内的若干已知样本点数据，在考虑了样本点的形状、大小和空间方位，与未知样点的相互空间关系，以及变异函数提供的结构信息之后，对未知样点进行的一种线性无偏最优估计。
假设研究区域a 上研究变量Z （x ），在点x i ∈A （i=1，2，……，n ）处属性值为Z （x i ），则待插点x 0∈A 处的属性值Z （x 0）的克里金插值结果Z*（x 0）是已知采样点属性值Z （x i ）（i=1，2，……，n ）的加权和，即：)()(10*i ni i x Z x Z ∑==λ （1） 式中i λ是待定权重系数。
其中Z(x i )之间存在一定的相关关系，这种相关性除与距离有关外，还与其相对方向变化有关，克里金插值方法将研究的对象称“区域化变量”针对克里金方法无偏、最小方差条件可得到无偏条件可得待定权系数i λ (i=1，2，……，n)满足关系式： 11=∑=n i i λ（2）以无偏为前提，kriging 方差为最小可得到求解待定权系数i λ的方程组：⎪⎪⎩⎪⎪⎨⎧=⋯⋯==+∑∑==1)n ,2,1)(,(),(101n i i j j i n i i j x x C x x C λμλ， （3） 式中，C （x i ，x j ）是Z(x i )和Z(x j )的协方差函数。
克里金模型详细推导
克里金模型详细推导## Kriging Model Derivation.Covariance Function:The kriging model is based on a statistical framework that incorporates spatial autocorrelation. This autocorrelation is captured through a covariance function, denoted by C(s, t), which specifies the correlation between the values at locations s and t. Common covariance functions include:Spherical:Exponential:Gaussian:\( C(s,t) = \sigma^2 \exp\left(-\frac{\|s-t\|^2}{2a^2}\right) \)。
where σ^2 is the process variance, and a is the range parameter that controls the distance at which values become uncorrelated.Kriging Equation:The kriging equation is used to predict the value of the underlying field at an unobserved location s0:\( Z_{s_0} = \sum_{i=1}^n \lambda_i Z(s_i) \)。
where:Z(s_i) are the observed values at locations s_i.λ_i are weights that sum to 1。
高斯过程回归模型 (kriging)
高斯过程回归模型 (kriging)高斯过程回归模型，也被称为kriging方法，是一种基于高斯过程的非参数回归技术。
它通过利用高斯过程对未知函数进行建模，并根据已观测到的数据点来估计未观测到的数据点的值。
在很多实际应用中，高斯过程回归模型被广泛应用于空间插值、地质建模、地理信息系统、环境工程、农业科学等领域。
高斯过程回归模型的基本假设是：给定任意输入x，对应的输出y满足一个联合高斯分布，即y ~ N(m(x), k(x, x'))，其中m(x)是均值函数，k(x, x')是协方差函数。
均值函数描述了数据的全局趋势，协方差函数描述了不同点之间的相关性。
在高斯过程回归模型中，对未观测到的数据点进行预测时，首先需要估计均值函数和协方差函数的参数。
常用的估计方法包括最大似然估计和贝叶斯推断。
通过优化似然函数，可以得到均值函数和协方差函数的最优参数。
然后，根据已观测到的数据点和估计得到的参数，可以通过贝叶斯推断方法，计算未观测数据点的后验分布，并进行预测。
在具体的算法实现中，高斯过程回归模型通常分为两个步骤：训练和预测。
在训练阶段，首先根据已知的输入和输出数据点，利用最大似然估计或贝叶斯推断方法，估计均值函数和协方差函数的参数。
然后，根据估计得到的参数，计算数据点之间的协方差矩阵，并将其分解为一个低秩矩阵和一个对角矩阵，以减少计算复杂度。
在预测阶段，根据已知的输入和输出数据点，利用训练阶段得到的参数，计算未观测数据点的条件分布，并进行预测。
高斯过程回归模型的优点之一是它能够提供预测结果的不确定性估计。
由于高斯过程的后验分布是一个高斯分布，可以通过计算均值和方差来描述预测结果的中心和离散程度。
这对于决策制定者来说非常重要，因为他们可以据此评估预测结果的可信度。
另一个优点是高斯过程回归模型的灵活性。
通过选择不同的均值函数和协方差函数，可以适应不同的数据特征和模型假设。
常用的协方差函数包括常值函数、线性函数、指数函数、高斯函数等。
kriging(克里金方法_克里金插值)[1]
![kriging(克里金方法_克里金插值)[1]](https://img.taocdn.com/s3/m/f5fb9b4e79563c1ec5da7181.png)
(h) C(0) C(h)
（二阶平稳假设条件下边查函数与写防查的关系）
变程(Range) ：指区域化变量在空间上具有相关性的 范围。在变程范围之内，数据具有相关性；而在变 程之外，数据之间互不相关，即在变程以外的观测 值不对估计结果产生影响。
具不同变程 的克里金插 值图象
块金值(Nugget) ：变差函数如果在原点间断，在地质统计学中称 为“块金效应”，表现为在很短的距离内有较大的空间变异性， 无论h多小，两个随机变量都不相关 。它可以由测量误差引起， 也可以来自矿化现象的微观变异性。在数学上，块金值c0相当于 变量纯随机性的部分。
Z*(x0)
（1）无偏条件
从本征假设出发, 可知 EZx为常数,有
EZ * x0 Zx0
E n i Z xi Z x0
i1
n i m m 0 i1
（在搜寻邻域内为 常数，不同邻域可 以有差别）
可出现E[Z(u)]不存在， 但E[Z(u)-Z(u+h)]存在并为零的情况
E[Z(u)]可以变化,但E[Z(u)-Z(u+h)]=0
② 增量[Z(u)-Z(u+h)]的方差函数 (变差函数,Variogram)
存在且平稳 (即不依赖于u)，即：
Var[Z(u)-Z(u+h)] = E[Z(u)-Z(u+h)]2-{E[Z(u)-Z(u+h)]}2 = E[Z(u)-Z(u+h)]2 = 2γ(u,h) = 2γ(h),
发表了专著《应用地质统计学论》。
阐明了一整套区域化变量的理论，
为地质统计学奠定了理论基础。
1977年我国开始引入
区域化变量理论 克里金估计 随机模拟
kriging三种模式
Homework#5 內插法黃琦聆M9505163（一）使用IDW，Spline，kriging三種模式，內插出三張GDP資料。
利用ARCGIS 可求出其內插值，其內插值與真實值之相關係數(Correation)與相關系數圖分別如下：1. IDW=0.9999852.Spline=0.999264（二） 內插法的model(計算方式)及寫出各個方法的設定係數。
1. IDW((Inverse Distance Weight ，距離反比權重法)針對每一個未知的數值推估，距離反比權重法是利用它鄰近的以經典之數值來進行加權運算，所給的權重依照距離遠近來計算，公式如下(數值等高線內插之比較研究，1996)：()⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡=∑∑==N i i N i i i d w z d w y x f 11)()(, 其中)(i d w 是權重方程，ｚｉ是第ｉ個已知點的數值，ｄｉ是ｉ點到未知點之間的距離。
在此分別設定power=2、Search radius 為variable 、Number of point=12、Output cell size=500。
2. Spline 法Spline 法是一個應用在選擇樣點後，利用多項式方進行內插法，產生平滑曲面的內插法。
ArcGIS 系統中，在進行Spline 模擬時需進行三個參數設定，包括了權重(Weight)、點數(Number of point)及曲面之型態(Type)；在此分別設定為Weight=0.1、Number of point=12、Type 為Regularized 、Output cell size=500，其中Regularized 的型式所模擬出的結果較使用另一type 為Tension 所模擬出的結果可產生較平滑的曲面(應用地理統計於土壤重金屬污染物之空間分佈探討，2004)。
Spline 內插法公式如下：(依據ARCGIS 公式說明))(),(),(1j Nj j r R y x T y x S ∑=+=λ其中j = 1, 2, ..., N (N is the number of points)λj are coefficients found by the solution of a system of linear equations--係線性方程式解的係數r j is the distance from the point (x,y) to the j th point--係點(x,y)到點 j th 的距離while T(x,y) and R(r) are defined differently, depending upon the selected option. For the REGULARIZED option:T(x,y) = a 1 + a 2x + a 3y⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡⎪⎭⎫ ⎝⎛++⎪⎭⎫ ⎝⎛+⎥⎦⎤⎢⎣⎡-+⎪⎭⎫ ⎝⎛=πτττπ2ln 12ln 421)(022r c r K c r r r R and for the TENSION option:T(x,y) = a 1()⎥⎦⎤⎢⎣⎡++⎪⎭⎫ ⎝⎛-=ϕϕπϕr K c r r R 022ln 21)( where,τ2and φ2are the parameters entered at the command line.r is the distance between the point ant the sample.K 0is the modified Bessel function.c is a constant equal to 0.577215.a i are coefficients found by the solution of a system of linear equations.3. Kriging 法Kriging 內插法是空間統計上的一項重要方法。
克里金插值公式推导
克里金插值公式推导克里金插值（Kriging Interpolation）是一种空间插值方法，它是由法国数学家达卡斯特罗（Georges Matheron）在1951年提出的。
克里金插值在地质学、环境科学、地理信息系统等领域有广泛的应用。
克里金插值的基本思想是通过已知离散点的观测值，推断和估计未知位置处的值。
它的特点是能够提供具有空间连续性的插值结果，并且能够提供对预测值的误差估计。
克里金插值的推导基于统计学中的协方差函数和高斯过程。
假设我们有 n 个观测点（x1, y1, z1）、（x2, y2, z2）、..、（xn, yn, zn），其中 (xi, yi) 为观测点的坐标，zi 为观测值。
我们需要推断和估计未知位置 (x0, y0) 处的值 z0。
首先，我们需要定义一个协方差函数 C(h)，其中 h 为两个点之间的距离。
协方差函数用来描述两个点之间的相关性，通常采用指数型（Exponential）、高斯型（Gaussian）或球型（Spherical）等函数形式。
接下来，我们可以利用协方差函数构建协方差矩阵 K。
协方差矩阵是一个对称正定矩阵，其元素 kij 表示点 i 和点 j 之间的协方差。
然后，我们需要定义一个权重函数W(x0,y0)，其中(x0,y0)是未知位置的坐标。
权重函数的作用是为未知位置处的值z0分配权重，权重与样本点之间的距离以及协方差函数的取值相关。
权重函数的形式可以根据具体问题的需求进行选择，常见的有逆距离权重法（Inverse Distance Weighting）和克里金权重法（Kriging Weighting）。
逆距离权重法主要考虑了样本点与未知位置之间的距离，而克里金权重法则同时考虑了距离和协方差。
最后，我们可以利用权重函数和观测值计算未知位置处的值z0。
根据克里金插值的思想，插值结果是观测值的加权平均，权重由权重函数给出。
具体的计算公式如下：z0=∑(Wi*Zi)其中，Wi表示未知位置与观测点i之间的权重，Zi表示观测点i的观测值。
克里金插值(kriging)
2021/6/16
12
二、统计推断与平稳要求
•任何统计推断（cdf，数学期望等）均要求重复取样。 •但在储层预测中，一个位置只能有一个样品。 •同一位置重复取样，得到cdf，不现实
P
2021/6/16
13
考虑邻近点，推断待估点
区域化变量： 能用其空间分布来表征一个自然现象的变量。
（将空间位置作为随机函数的自变量）
•空间一点处的观测值可解释为一个随机变量在该点
处的一个随机实现。
• 空间各点处随机变量的集合构成一个随机函数。
（可以应用随机函数理论解决插值和模拟问题）
2021/6/16
14
考虑邻近点，推断待估点 ----空间统计推断要求平稳假设
方差的平方根为标准差，记为σξ
σξ=
D ()E [-E ()]2 E (2 )-[E ()]2
•从矩的角度说，方差是ξ的二阶中心矩。
2021/6/16
10
2. 随机函数
研究范围内的一组随机变量。
{Z(u),u研究范}围 简记为 Z (u )
条件累积分布函数（ccdf）
F ( u 1 , , u K ; z 1 , , z K |( n ) P ) o { Z ( u r 1 ) b z 1 , , Z ( u K ) z K |( n )}
2021/6/16
2
H. S. Sichel (1947) D.G. Krige (1951)
应用统计学方法研究金矿品位
Kriging法（克里金法，克立格 法）:“根据样品空间位置不同、样 品间相关程度的不同，对每个样品 品位赋予不同的权，进行滑动加权 平均，以估计中心块段平均品位”
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
首先假设区域化变量 Z ( x)满足二阶平稳假 设和本征假设，其数学期望为 m ，协方差函数 及变异函数 c(h) 存在。 (h) 即
E[Z ( x)] m
c(h) E[Z ( x)Z ( x h)] m2
1 (h) E[ Z ( x) Z ( x h)] 2 2
10.3 10.5 10.9 11.2
4.9
5.1 6.2 7.5 9.5
8.8
3.8
12.4
9.8
从上面的介绍和讨论，我们知道，球状 变异函数的一般形式为
0 3h h 3 (h) c 0 c( 3 ) 2a 2a c 0 c h0 0ha ha
当 0 h a 时，有
11 12 21 22 K n1 n 2 1 1 1n 2n nn 1 1 1 , 1 0 1 2 , n ( x1 , x) ( x , x) 2 D ( x , x ) n 1
2 E i 1 j 1 i 1 n n n
(4.2.24) 为使估计方差最小，根据拉格朗日乘数原 理，令 n 2 F E 2 ( i 1) (4.2.25) i 1 求F对 i 和 的偏导数，并令其为0，得克 立格方程组
n F 2 j c( xi , x j ) 2c( xi , x) 2 0 j 1 i n F 2( 1) 0 i i 1
K D
(4.2.34)
(4.2.35) K 1 D 2 (4.2.36) K T D ( x, x) 在以上的介绍中，区域化变量 Z ( x) 的数学 期望E[Z ( x)] m可以是已知或未知的。如果m是 已知常数，称为简单克立格法；如果m是未 知常数，称为普通克立格法。不管是哪一种 方法，均可根据方法计算权重系数和克立格 估计量。
3c c ( h) c 0 ( ) h ( 3 ) h 3 2a 2a
3c 1 c ,b , x h, x h ，则可以 如果记 y (h), b c , b 2 a 2a 得到线性模型 y b0 b1 x1 b2 x2 (4.2.19)
(4.2.29)
解线性方程组 (4.2.27) 式，求出权重系数 和拉格朗日乘数μ ，代入公式(4.2.24)，可得克立 格估计方差 ，即
2 K i ( xi , x) ( x, x) i 1 n
(4.2.30矩阵形式。令
* Z0 0.287Z ( x1 ) 0.210Z ( x2 ) 0.202Z ( x3 ) 0.473Z ( x4 )
0.287 37 0.210 42 0.202 36 0.473 35 37.25（mm）。
3 0 0 1 2 3 1 2
根据表中的数据，对上式进行最小二乘拟合
它是一种数学上的近似和优化，利用已知的数据得 出一条直线或者曲线，使之在坐标系上与已知数据之 间的距离的平方和最小。，得到
y 2.048 1.731x1 0.007 92x2
(4.2.20)
比较 (4.2.20) 式与 (4.2.19) 式，并做简 单计算可知： c0=2.048 ， c=1.154 ， a=8.353 ， 所以，球状变异函数模型为
地统计（Geostatistics）又称地质统计，它是以区域化变
量为基础，借助变异函数，研究既具有随机性又具有结构 性，或空间相关性和依赖性的自然现象的一门科学。凡是 与空间数据的结构性和随机性，或空间相关性和依赖性， 或空间格局与变异有关的研究，并对这些数据进行最优无 偏内插估计，或模拟这些数据的离散性、波动性时，皆可 应用地统计学的理论与方法。 地统计方法的条件：(1)随机过程 (2)正态分布 （3)平稳性： ①均值平稳 ②协方差函数有关的二阶平稳和与半变异函数有关的内蕴 平稳 Kriging 插值的条件：变异函数和相关分析的结果表明区域 化变量存在空间相关性。
1
0.952 0.287 0.711 0.210 0.571 0.202 0 . 870 0 . 301 1 0.473
即 克 立 格 权 重 系 数 分 别 为 ： λ1=0.287 ， λ2=0.210,λ3=0.202,λ4=0.301 ， μ= -0.473 ，所以观 测点的降水量的克立格估计值为：根据普通克立格 法的基本原理，我们知道， Z(x0) 估计的基本公式 应该是
2 σE c( x, x) i c( xi , x) i 1
n
(4.2.28)
上述过程也可用矩阵形式表示，令
c11 c 21 K cn1 1 c12 c1n c22 c2 n cn 2 cnn 1 1 1 1 , 1 0 1 2 , n c( x1 , x) c ( x , x ) 2 D c ( x , x ) n 1
克里格插值基础
克里格方法概述
克里格方法（Kriging）又称空间局部插值法， 是以变异函数理论和结构分析为基础，在有限 区域内对区域化变量进行无偏最优估计的一种 方法，是地统计学的主要内容之一。其实质是 利用区域化变量的原始数据和变异函数的结构 特点，对未知样点进行线性无偏、最优估计。 无偏是指偏差的数学期望为0，最优是指估计值 与实际值之差的平方和最小。
c24 c42 3.202 ( 32 2 2 ) 0.466
c 01 3.202 ( 12 ) 0.952
c 03 3.202 ( 3 2 ) 0.571
将以上计算结果代入克立格方程组 （4.2.31 ），得
1 3.202 0.870 2 3 0.542 4 0.711 1 0.870 3.202 0.601 0.466 1 0.542 0.601 3.202 0.383 1 0.711 0.466 0.383 3.202 1 1 1 1 1 0
c13 c31 3.202 ( 3 2 12 ) 0.542
c14 c 41 c 02 3.202 ( 2 2 12 ) 0.711
c 23 c32 3.202 ( 2 2 2 2 ) 0.601
c34 c 43 3.202 ( 4 2 12 ) 0.383
(1)无偏性。要使 Z * ( x) 成为 Z ( xi ) 的无偏估计量， 即
E[Z * ( x)] E[Z ( x)]
当 E[Z ( x)] m 时，也就是当 E [ i Z ( xi )] i E[Z ( xi )] m i 1 i 1 时，则有
n
n
i 1
0 3 3 h 1 h * (h) 2.048 1.154( ) 3 2 8.535 2 8.535 3.202
(4.2.21)
h0 0 h 8.535 h 8.535
4个观测点x1，x2，x3，x4的观测值分别为 Z(x1)=37、Z(x2)=42、Z(x3)=36、Z(x4)=35，如果 假设降水量的变异函数是向同性（即变异函数 在各个方向的变化都相同）的二维球状模型， 其具体形式为（4.2.21）式。现在，我们用普 通克立格法估计观测点x0的降水量值Z(x0)。
(4.2.26)
整理后得
n j c( xi , x j ) c( xi , x) j 1 n i 1 i 1
(4.2.27)
解线性方程组(4.2.27)式,求出权重系数λ i和拉格朗日系数μ ,代入公 式(4.2.24),可得克立格估计方差
假设在待估计点（x）的临域内共有n个 实测点，即x1，x2 ，… ，xn ，其样本值为 Z ( xi )。 那么，普通克里格法的插值公式为
Z * ( x) i Z ( xi )
i 1 n
(4.2.22)
其中 i 为权重系数 ，表示各空间样本点 x i 处的观测值Z ( xi ) 对估计值Z * ( x) 的贡献程度。 可见，克立格插值的关键就是计算权重 系数 i 。显然，权重系数的求取必须满足两个 条件： 一是使Z * ( x)的估计是无偏的，即偏差的 数学期望为零； 二是最优的，即使估计值 Z * ( x) 和实际值Z ( xi ) 之差的平方和最小。 为此，需要满足以下两个条件：
图10.26 克里格方法流程 图
例如:某地区降水量是一个区域化变量，其变 异函数 (h)的实测值及距离h的关系见下表， 下面我们试用回归分析方法建立其球状变异 函数模型。
实测值γ(h) 距离h 实测值γ(h) 距离h
2.1
4.3 5.7 6.5 7.8
0.6
1.1 2.2 2.5 3.1
9.2
克里格插值基础
2.克里格方法的具体步骤
导入数据 进行预测
数据分析
否
计算克里格系数 数据变换
是否服从 正态分布 是 是否存 在趋势 否 根据数据选择 合适的方法
拟合理论半变 异函数图 绘制经验半变 异函数图
是 泛克里格方法
绘制方差变 异云图 按组统计平均距离 及对应的平均方差 按距离分组
计算样点间的 距离矩阵 计算样点间的 属性方差
n
i
1
这时， Z * ( x) 为Z ( xi ) 的无偏估计量。
（2）最优性。在满足无偏性条件下，估计方差
为
E[Z ( x) Z ( x)] E[Z ( x) i Z ( xi )] 2
2 E * 2 i 1 n
使用协方差函数表达，它可以进一步写 为