克里金插值法的详细介绍。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.估计误差的方差达到最小。
克里金插值法的详细介绍。kriging。
kriging 插值作为地统计学中的一种插值方法由南非采矿工程师D.G.Krige于1951年首次提出，是一种求最优、线形、无偏的空间内插方法。
在充分考虑观测资料之间的相互关系后，对每一个观测资料赋予一定的权重系数，加权平均得到估计值。
这里介绍普通Kriging插值方法的基本步骤：1.该方法中衡量各点之间空间相关程度的测度是半方差，其计算公式为：h为各点之间距离，n 是由h 分开的成对样本点的数量，z 是点的属性值。
2.在不同距离的半方差值都计算出来后，绘制半方差图，横轴代表距离，纵轴代表半方差。
半方差图中有三个参数nugget（表示距离为零时的半方差）,sill（表示基本达到恒定的半方差值）,range(表示一个值域范围，在该范围内半方差随距离增加，超过该范围，半方差值趋于恒定)。
利用做出的半方差图找出与之拟合的最好的理论变异函数模型（这是关键所在），可用于拟合的模型包括高斯模型、线性模型、球状模型、指数模型、圆形模型。
----球状模型，球面模型空间相关随距离的增长逐渐衰减，当距离大于球面半径后，空间相关消失。
3.用拟合的模型计算出三个参数。
例如球状模型中nugget为c0，range为a，sill为c。
4．利用拟合的模型估算未知点的属性值，方程为：，z0为估计值，zx是已知点的值，wx为权重，s是用来估算未知点的已知点的数目。
假如用三个点来估算，则有这样权重就可以求出，然后估算未知点。
（上述内容根据《地理信息系统导论》（Kang-tsung Chang著；陈健飞等译，科学出版社，2003）第十三章内容进行总结，除球状模型公式外其余公式皆来自此书）下面是本人自己编写的利用海洋中断面上观测站点的实测温度值来估算未观测处的温度的Fortran程序，利用距离未知点最近的五个观测点来估算未知点的温度，选用模型为球状模型。
do ii=1,nxif(tgrid(ii,1)==0.)thendo i=1,dsite(ii)!首先寻找距离最近的五个已知点位置do j=1,nhif(d(mm(ii),j).ne.0.or.j==1)thenhmie(j)=d(mm(ii),j)-dgrid(i)elsehmie(j)=9999end ifhmid(j)=abs(hmie(j))end dodo j=1,nhdo k=j,nhif(hmid(j)<hmid(k))thenelsem1=hmid(j)hmid(j)=hmid(k)hmid(k)=m1end ifend doend dodo j=1,5do k=1,nhif(abs(hmie(k))==hmid(j))thenlocat(j)=kend ifend doend dodo j=1,4do k=j+1,5if(locat(j)==locat(k))thendo i3=1,nhif(abs(hmie(i3))==abs(hmie(locat(j))).and.i3.ne.locat(j))thenlocat(j)=i3exitend ifenddoendifenddoenddo！然后求各点间距离，并求半方差do j=1,5do k=1,5hij(j,k)=abs(d(mm(ii),locat(j))-d(mm(ii),locat(k)))/1000.end doend dodo j=1,5hio(j)=sqrt(hmid(j)**2+(abs(latgrid(ii)-lonlat(mm(ii),2))*llat)**2 $ +(abs(longrid(ii)-lonlat(mm(ii),1))*(1.112e5* $ cos(0.017*(latgrid(ii)+lonlat(mm(ii),2))/2)))**2)/1000.end dodo j=1,5do k=1,5if(hij(j,k).eq.0.)thenrleft(j,k)=0.elserleft(j,k)=sill*(1.5*hij(j,k)/range-0.5*hij(j,k)**3/range**3)end ifif(hio(j).eq.0.)thenrrig(1,j)=0.elserrig(1,j)=sill*(1.5*hio(j)/range-0.5*hio(j)**3/range**3)end ifend doend dorrig(1,6)=1.rleft(6,6)=0.do j=1,5rleft(6,j)=1.rleft(j,6)=1.end dotry=rleftcall brinv(rleft,nnn,lll,is,js)ty1=matmul(try,rleft)！求权重wq=matmul(rrig,rleft)!插值所有格点上t,sdo j=1,5tgrid(ii,i)=tgrid(ii,i)+wq(1,j)*t(mm(ii),locat(j)) sgrid(ii,i)=sgrid(ii,i)+wq(1,j)*s(mm(ii),locat(j))end doenddoendifenddo。
克里金插值(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法（克里金法，克立格 法）:“根据样品空间位置不同、样 品间相关程度的不同，对每个样品 品位赋予不同的权，进行滑动加权 平均，以估计中心块段平均品位”
克里金插值(kriging)
随机场：
P
当随机函数依赖于多个
自变量时，称为随机场。
如具有三个自变量(空间
点的三个直角坐标)的随
机场
随机函数的特征值
协方差(Variance)： 二个随机变量ξ，η的协方差为二维随机变量(ξ，
η)的二阶混合中心矩μ11，记为Cov(ξ，η)，或σξ,η。
Cov(ξ,η) = σξ,η = E[ξ-E(ξ)][η-E(η)]
块金效应的尺度效应
如果品位完全是典型的随机变量，则不论 观测尺度大小，所得到的实验变差函数曲线总 是接近于纯块金效应模型。
当采样网格过大时，将掩盖小尺度的结构， 而将采样尺度内的变化均视为块金常数。这种 现象即为块金效应的尺度效应。
1
3
3
3
1
2
3
1
1
(h) = C(0) – C(h)
基台值(Sill)：代表变量在空间上的总变异性大小。即为变 差函数在h大于变程时的值，为块金值c0和拱高cc之和。 拱高为在取得有效数据的尺度上，可观测得到的变异性幅 度大小。当块金值等于0时，基台值即为拱高。
对于单变量而言：
P
F(u;z)F(uh;z)
可从研究区内所有数据的累积直方图推断而得 （将邻近点当成重复取样点）
太强的假设，不符合实际
二阶平稳
当区域化变量Z(u)满足下列二个条件时，则称其 为二阶平稳或弱平稳：
① 在整个研究区内有Z(u)的数学期望存在， 且等于常数，即： E[Z(u)] = E[Z(u+h)] = m(常数) x h
相当于要求：Z(u)的变差函数存在且平稳。
可出现协方差函数不存在，但变差函数存在的情况。
例：物理学上的著名的布朗运动是一种呈现出无限 离散性的物理现象，其随机函数的理论模型就是维 纳-勒维(Wiener-Levy)过程(或随机游走过程)。
克里金插值
克里金(Kriging)插值克里金(Kriging)插值法又称空间自协方差最佳插值法，它是以南非矿业工程师D．G．Krige的名字命名的一种最优内插法。
克里金法广泛地应用于地下水模拟、土壤制图等领域，是一种很有用的地质统计格网化方法它首先考虑的是空间属性在空间位置上的变异分布．确定对一个待插点值有影响的距离范围，然后用此范围内的采样点来估计待插点的属性值。
该方法在数学上可对所研究的对象提供一种最佳线性无偏估计(某点处的确定值)的方法。
它是考虑了信息样品的形状、大小及与待估计块段相互间的空间位置等几何特征以及品位的空间结构之后，为达到线性、无偏和最小估计方差的估计，而对每一个样品赋与一定的系数，最后进行加权平均来估计块段品位的方法。
但它仍是一种光滑的内插方法在数据点多时，其内插的结果可信度较高。
克里金法类型分常规克里金插值(常规克里金模型／克里金点模型)和块克里金插值。
常规克里金插值其内插值与原始样本的容量有关，当样本数量较少的情况下，采用简单的常规克里金模型内插的结果图会出现明显的凹凸现象；块克里金插值是通过修改克里金方程以估计子块B内的平均值来克服克里金点模型的缺点，对估算给定面积实验小区的平均值或对给定格网大小的规则格网进行插值比较适用。
块克里金插值估算的方差结果常小于常规克里金插值，所以，生成的平滑插值表面不会发生常规克里金模型的凹凸现象。
按照空间场是否存在漂移(drift)可将克里金插值分为普通克里金和泛克里金，其中普通克里金(Ordinary Kriging简称OK法)常称作局部最优线性无偏估计．所谓线性是指估计值是样本值的线性组合，即加权线性平均，无偏是指理论上估计值的平均值等于实际样本值的平均值，即估计的平均误差为0，最优是指估计的误差方差最小。
在科学计算领域中，空间插值是一类常用的重要算法，很多相关软件都内置该算法，其中GodenSoftware 公司的Surfer软件具有很强的代表性，内置有比较全面的空间插值算法，主要包括：Inverse Distance to a Power（反距离加权插值法）Kriging（克里金插值法）Minimum Curvature（最小曲率）Modified Shepard's Method（改进谢别德法）Natural Neighbor（自然邻点插值法）Nearest Neighbor（最近邻点插值法）Polynomial Regression（多元回归法）Radial Basis Function（径向基函数法）Triangulation with Linear Interpolation（线性插值三角网法）Moving Average（移动平均法）Local Polynomial（局部多项式法）下面简单说明不同算法的特点。
克里金插值方法介绍
特殊地，当h=0时，上式变为 Var[Z(u)]=C(0), 即方差存在且为常数。
u+h u
本征假设 intrinsic hypothese
（比二阶平稳更弱的平稳假设）
当区域化变量Z(u)的增量[Z(u)-Z(u+h)]满足下列二 条件时，称其为满足本征假设或内蕴假设。
①在整个研究区内有 E[Z(u)-Z(u+h)] = 0
半变差函数(或半变异函数)
在二阶平稳假设，或作本征假设，此时：
E[Z(x)-Z(x+h)] = 0 h
则：
(x,h) =
1
2 Var[Z(x)-Z(x+h)]
=
1 2
E[Z(x)-Z(x+h)]2-{E[Z(x)-Z(x+h)]}2
(x,h)
=
1 2
E[Z(x)-Z(x+h)]2
地质统计学中最常用 的基本公式之一。
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
当随机函数不满足二阶平稳,而满足内蕴（本征）假设时, 可用变差函数来表示克里金方程组如下:
•在实际变程处，变差函 数为0.95c。
•模型在原点处为抛物线。
幂函数模型：
h c.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值上下波动。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
kriging 插值作为地统计学中的一种插值方法由南非采矿工程师D.G.Krige于1951年首次提出，是一种求最优、线形、无偏的空间内插方法。
在充分考虑观测资料之间的相互关系后，对每一个观测资料赋
予一定的权重系数，加权平均得到估计值。
这里介绍普通Kriging插值方法的基本步骤：1.该方法中衡量各点之间空间相关程度的测度是半方差，
其计算公式为：
h为各点之间距离，n 是由h 分开的成对样本点的数量，z 是点的属性值。
2.在不同距离的半方差值都计算出来后，绘制半方差图，横轴代表距离，纵轴代表半方差。
半方差图中有三个参数nugget（表示距离为零时的半方差）,sill（表示基本达到恒定的半方差值）,range(表示一个值域范围，在该范围内半方差随距离增加，超过该范围，半方差值趋于恒定)。
利用做出的半方差图找出与之拟合的最好的理论变异函数模型（这是关键所在），可用于拟合的模型包括高斯模型、线性模型、球状
模型、指数模型、圆形模型。
----球状模型，球面模型空间相关随距离的增长逐渐衰减，当距离大于球面半径后，空间相关消失。
3.用拟合的模型计算出三个参数。
例如球状模型中nugget为c0，range为a，sill为c。
4．利用拟合的模型估算未知点的属性值，方程为：
，z0为估计值，zx是已知点的值，wx为权重，s是用来估算未知点的
已知点的数目。
假如用三个点来估算，则有
这样权重就可以求出，然后估算未知点。
（上述内容根据《地理信息系统导论》（Kang-tsung Chang著；陈健飞等译，科学出版社，2003）第十三章内容进行总结，除球状模型公式外其余公式皆来自此书）
下面是本人自己编写的利用海洋中断面上观测站点的实测温度值来估算未观测处的温度的Fortran程序，利用距离未知点最近的五个观测点来估算未知点的温度，选用模型为球状模型。
do ii=1,nx
if(tgrid(ii,1)==0.)then
do i=1,dsite(ii)
!首先寻找距离最近的五个已知点位置
do j=1,nh
if(d(mm(ii),j).ne.0.or.j==1)then
hmie(j)=d(mm(ii),j)-dgrid(i)
else
hmie(j)=9999
end if
hmid(j)=abs(hmie(j))
end do
do j=1,nh
do k=j,nh
if(hmid(j)<hmid(k))then
else
m1=hmid(j)
hmid(j)=hmid(k)
hmid(k)=m1
end if
end do
end do
do j=1,5
do k=1,nh
if(abs(hmie(k))==hmid(j))then
locat(j)=k
end if
end do
end do
do j=1,4
do k=j+1,5
if(locat(j)==locat(k))then
do i3=1,nh
if(abs(hmie(i3))==abs(hmie(locat(j))).and.i3.ne.locat(j))then
locat(j)=i3
exit
end if
enddo
endif
enddo
enddo
！然后求各点间距离，并求半方差
do j=1,5
do k=1,5
hij(j,k)=abs(d(mm(ii),locat(j))-d(mm(ii),locat(k)))/1000.
end do
end do
do j=1,5
hio(j)=sqrt(hmid(j)**2+(abs(latgrid(ii)-lonlat(mm(ii),2))*llat)**2 $ +(abs(longrid(ii)-lonlat(mm(ii),1))*(1.112e5* $ cos(0.017*(latgrid(ii)+lonlat(mm(ii),2))/2)))**2)/1000.
end do
do j=1,5
do k=1,5
if(hij(j,k).eq.0.)then
rleft(j,k)=0.
else
rleft(j,k)=sill*(1.5*hij(j,k)/range-0.5*hij(j,k)**3/range**3)
end if
if(hio(j).eq.0.)then
rrig(1,j)=0.
else
rrig(1,j)=sill*(1.5*hio(j)/range-0.5*hio(j)**3/range**3)
end if
end do
end do
rrig(1,6)=1.
rleft(6,6)=0.
do j=1,5
rleft(6,j)=1.
rleft(j,6)=1.
end do
try=rleft
call brinv(rleft,nnn,lll,is,js)
ty1=matmul(try,rleft)
！求权重
wq=matmul(rrig,rleft)
!插值所有格点上t,s
do j=1,5
tgrid(ii,i)=tgrid(ii,i)+wq(1,j)*t(mm(ii),locat(j)) sgrid(ii,i)=sgrid(ii,i)+wq(1,j)*s(mm(ii),locat(j))
end do
enddo
endif
enddo。