克里格插值法

克里格插值什么是克里格插值？距离权重倒数插值和样条法插值被归类为确定性的插值方法，因为它们是直接基于周围已知点的值进行计算或是用指定的数学公式来决定输出表面的平滑度的插值方法。

而第二个插值方法家族包括的是一些地统计学的插值方法（如克里格插值），这些方法基于一定的包括诸如自相关（已知点间的统计关系）之类的统计模型。

因此，这些方法不仅有能力生成一个预测表面，而且还可以给出预测结果的精度或确定性的度量。

克里格插值与距离权重倒数插值相似之处在于给已知的样本点赋权重来派生出未知点的预测值。

这两种内插方法的通用公式如下，表达为数据的权重总和。

其中， Z(Si)是已测得的第i个位置的值；λi是在第i个位置上测得值的未知的权重；S0是预测的位置；N 是已知点（已测得值的点）的数目。

在距离权重倒数插值中，权重λi仅取决于距预测位置的距离。

然而，在克里格插值中，权重不仅建立在已知点和预测点位置间的距离的基础上，而且还要依据已知点的位置和已知点的值的整体的空间分布和排列。

应用权重的空间排列，空间自相关必须量化。

因此，运用普通克里格插值（Ordinary Kriging），权重λi取决于已知点的拟合模型、距预测位置的距离和预测点周围的已知点间的空间关系。

利用克里格方法进行预测，必须完成以下两个任务：（1）揭示相关性规则。

（2）进行预测。

要完成这两项任务，克里格插值方法通过以下两个步骤完成：（1）生成变异函数和协方差函数，用于估算单元值间的统计相关（也叫空间自相关），而变异函数和协方差函数也取决于自相关模型（拟合模型）。

（2）预测未知点的值。

因为前面已经说过的两个明确的任务，因此要用克里格方法对数据进行两次运算：第一次是估算这些数据的空间自相关而第二次是做出预测。

变异估计（Variography）变异估计就是拟合一个数学模型或空间模型，象已知的结构分析。

在已测点结构的空间建模中，首先得出经验半变异函数的曲线图，计算如下：半变异函数（距离h）＝ 0.5*均值[ (在i 位置的值－在j 位置的值)2 ]用于计算被距离h分隔的每一点对相对应的位置。

0x 克里格（Kringing ）插值法是建立在统计学理论基础上，实际上是利用区域化变量的原始数据和半方差数据的结构特征，对位采样点的区域化变量的取值进行线性最优无偏估计的一种方法，也就是根据待估样点有限领域内若干已经择定的测定的样点数据，在认真考虑了阳电的形状、大小和相互空间位置之间的关系，以及他们与待估样点见相互位置关系和编译函数提供的结构信息之后，对待估样点间相互位置关系的编译函数提供的结构信息之后，对待估样点值进行的一种线性最优无偏估计。

下图为运用克里格法计算未知点的值的一般步骤：其插值原理如下：设在某一研究内未知点0x 的属性为)(0x Z ，其周围相关范围内有n 个已知已测点),,2,1(n i x i ⋯=。

通过n 个测定值的线性组合求其估计值)(0x Z ：)()(10i n i i x Z x Z ∑==λ式中i λ为)(i x Z 位置有关的加权系数，并且∑==ni i 11λ克里格插值法是根据无偏估计和方差最小的要求来确定上式中的系数i λ。

1.构造半变异系数：设j x 和i x 的距离问为h 。

设n 个样点中mh 对样点的距离为h ，以他们的含量差)(-)(i j x Z x Z 构造的半变异函数为：2))()((21)(∑=--=h x x i j i j x Z x Z m h a 2.拟合得出变异系数：将n 个样点的含量带入公式，使用直线函数进行拟合3.构造矩阵和向量：求出任意两个已知点的半变异函数值，构造矩阵A:⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⋯⋯⋯⋯⋯⋯⋯⋯⋯=011110101021221112n n n n a a a a a a A 取任意一个已知点i x ，求出与未知点0x 的距离并代入求出该点与未知点0x 的半变异函数值0i a ，得到向量B:)1,,,,(02010n a a a B ⋯=方程AX=B 的姐的前n 个分量即为公式（）的权重系数i λ。

kriging插值法c代码克里格插值法（Kriging）是一种空间插值方法，用于根据已知数据点的位置和值来预测未知数据点的值。

以下是一个简单的克里格插值法的 C 语言代码示例：```c#include <stdio.h>#include <stdlib.h>#include <math.h>#define MAX_POINTS 100// 定义结构体来存储数据点的位置和值typedef struct {double x;double y;double z;} DataPoint;// 计算数据点之间的距离double calculateDistance(DataPoint p1, DataPoint p2) {double dx = p1.x - p2.x;double dy = p1.y - p2.y;return sqrt(dx * dx + dy * dy);}// 计算克里格插值的权重double calculateKrigingWeight(DataPoint *points, int nPoints, DataPoint queryPoint) {double sumWeight = 0.0;for (int i = 0; i < nPoints; i++) {DataPoint point = points[i];double distance = calculateDistance(queryPoint, point);double weight = 1.0 / (1.0 + distance * distance);sumWeight += weight;}return sumWeight / nPoints;}// 进行克里格插值double krigingInterpolation(DataPoint *points, int nPoints, DataPoint queryPoint) {double sumZ = 0.0;double sumWeight = 0.0;for (int i = 0; i < nPoints; i++) {DataPoint point = points[i];double distance = calculateDistance(queryPoint, point);double weight = 1.0 / (1.0 + distance * distance);sumZ += point.z * weight;sumWeight += weight;}return sumZ / sumWeight;}int main() {// 定义数据点的数量int nPoints = 5;// 分配内存来存储数据点DataPoint *points = (DataPoint *)malloc(nPoints * sizeof(DataPoint));// 输入数据点的位置和值points[0].x = 1.0; points[0].y = 2.0; points[0].z = 10.0;points[1].x = 2.0; points[1].y = 3.0; points[1].z = 15.0;points[2].x = 3.0; points[2].y = 4.0; points[2].z = 7.0;points[3].x = 4.0; points[3].y = 5.0; points[3].z = 12.0;points[4].x = 5.0; points[4].y = 6.0; points[4].z = 9.0;// 定义查询点的位置DataPoint queryPoint;queryPoint.x = 3.5; queryPoint.y = 4.5;// 计算克里格插值的权重double krigingWeight = calculateKrigingWeight(points, nPoints, queryPoint);// 进行克里格插值double interpolatedValue = krigingInterpolation(points, nPoints, queryPoint);// 输出结果printf("克里格插值的权重: %.2f\n", krigingWeight);printf("克里格插值的结果: %.2f\n", interpolatedValue);// 释放内存free(points);return 0;}```这段代码实现了一个简单的克里格插值法。

克里格估值方法(一)克里格估值方法详解什么是克里格估值法？克里格估值法（Kriging）是一种通过插值方法对未知地点进行估值的统计技术。

它将已知地点上的观测值用于预测未知地点上的数值，常用于地质、地理、环境等领域的研究。

克里格估值法通过建立空间相关性模型，可以提供对未知地点上现象的可信度估计。

克里格估值法的基本原理克里格估值法的基本原理是空间相关性。

其假设对空间上相邻点之间的值存在一定的相关性，且该相关性可通过距离进行量化。

基于该假设，克里格估值法可以通过已知点与未知点之间的空间距离进行权重的计算，进而进行预测。

克里格估值法的步骤1.数据获取：克里格估值法需要已知点的观测值作为输入，可以通过采集现有数据或者实地测量获得。

2.空间相关性分析：通过观测值之间的空间相关性判断模型类型，常用的模型包括球型模型、指数模型和高斯模型等。

3.参数估计：使用已知观测值中的半方差数据，通过最小二乘法或最大似然法对模型的空间相关参数进行估计。

4.半方差图绘制：通过绘制半方差图，可以了解观测值之间的空间相关性和变化趋势。

5.克里格估值：根据已知点的观测值和模型的参数，计算未知点上的估值。

常用的克里格估值方法包括简单克里格法、普通克里格法和泛克里格法等。

6.估值验证：通过验证估值和实际值之间的误差，评估克里格估值方法的精度和可靠性。

克里格估值法的优缺点克里格估值法作为一种插值方法具有以下优点： - 利用空间相关性进行预测，能够充分利用已知数据的信息； - 通过建立空间模型，可以对估值进行可靠的分析和解释； - 适用于各种数据类型和标度水平，可用于多种研究领域。

然而，克里格估值法也存在一些缺点： - 对观测值的空间相关性要求较高，如果空间相关性较弱，克里格估值的精度可能较低； - 克里格估值法对异常值敏感，对异常值进行处理是很重要的一步； - 克里格估值法无法考虑其他外部因素的影响，如地形、土壤等因素。

克里格估值法的应用领域克里格估值法广泛应用于地理信息系统（GIS）、环境调查和资源评价等领域，常见的应用包括： - 土壤污染程度评估； - 水资源管理及水质预测； - 土地利用规划和生态环境研究； - 地质勘探和矿产资源评估。

克立格插值（2）克立格插值（）空间信息统计深入 1.泛克立格方法的原理及应用 2.协同克立格方法的原理及应用1 泛克立格方法(universal kriging) ? 泛克立格法是用于解决非平稳区域化变量的插值方法? Kriging with a trend model y 漂移m(x) 实测值Z(x) x 1.1 与泛克立格方法有关的概念空间随机变量可以分解为：Z (x) = m (x) + R (x) 其中m(x)称为漂移，R(x)称其中m(x)称为漂移，R(x)称m(x)称为漂移为剩余漂移与剩余的平稳性分析剩余R(x)的数学期望不随空剩余R(x)的数学期望不随空R(x) 间位置的变化而变化，间位置的变化而变化，是平稳的区域化变量。

稳的区域化变量。

漂移m(x) m(x)随空间位置的不同漂移m(x)随空间位置的不同而变化，而变化，是非平稳的区域化变量。

变量。

非平稳条件下空间随机变量的数学期望m (x) = ∑ k l=0 al fl ( x) m ( x ) = a 0 + a1 x + a 2 y m ( x ) = a 0 + a 1 x + a 2 y + a 3 x 2 + a 4 xy + a 5 y 2 2.2 泛克立格方法可解决的问题? 求区域化变量的漂移求区域化变量的漂移m(x) ? 求x处的处的Z(x)值处的值? 求区域化变量的漂移系数(1) 求区域化变量的漂移求区域化变量的漂移m(x) m ( x) = ∑ ρα Z ( xα ) * n α =1 =1 k m ( x ) = ∑ al f l ( x ) l =0 无偏性条件的推导E*m ( x)+ = E*m( x)+ = m( x) * E*m ( x)+ = E*∑ ρα Z ( xα )+ = ∑ ρα E* Z ( xα )+ * n n α =1 α =1 = ∑ ρα E*m( xα )+ = ∑ ρα ∑ al f l ( xα ) α =1 k n n k α =1 l =0 = ∑ al ∑ ρα f l ( xα ) l =0 a =1 n 无偏条件E[m* ( x)] = m( x) 需要在任何条? 由于件下都成立，即：∑ a ∑ ρα f ( xα ) = ∑ a f ( x) l =0 l a =1 l l =0 l l k n k ∑ ρα f ( xα ) = f ( x) a =1 l l n 估计方差的推导Var * m ( x ) ? m * ( x )+ = E,m * ( x ) ? E* m * ( x )+-2 = E,∑ ρα Z ( xα ) ? ∑ ρα E* Z ( xα )+-2 α =1 n n n α =1 = E,∑ ρα * Z ( xα ) ? E ( Z ( xα ))+-2 α =1 n = E,∑ ρα * Z ( xα ) ? E ( Z ( xα ))+-, ∑ ρ β * Z ( x β ) ? E ( Z ( xβ ))+- α =1 n n β =1 = ∑ ∑ ρα ρ β E,Z ( xα ) ? E* Z ( xα )+-, Z ( xβ ) ? E* Z ( xβ )+- α =1 β =1 n n n = ∑ ∑ ρα ρ β Cov ( xα , xβ ) α =1 β =1 方差最小条件Q = ∑∑ ρα ρ βCov( xα , xβ ) ? 2∑ ?l *∑ ρα f l ( xα ) ? f l ( x)+ α =1 β =1 l =0 n n k n α =1 n k ?Q = 2∑ ρ β Cov( xα , xβ ) ? 2∑ ?l f l ( xα ) = 0 (α = 1,2, L n) ?ρα l =0 β =1 n ?Q = 2*∑ ρα f l ( xα ) ? f l ( x)+ = 0 (l = 0,1,2,L k ) ??l α =1 k ?n ?∑ ρ β Cov( xα , xβ ) ?∑ ?l f l ( xα ) = 0 (α = 1,2, L n) ? β =1 l =0 ? n ? ∑ ρα fl ( xα ) ? fl ( x) = 0 (l = 0,1,2,L k ) ? α =1 ? 估计方差δ *m( x)+ = ∑ ?l f l ( x) 2 k l =0 k 处的Z(x)值（2）求x处的）处的值Z UK ( x ) = ∑ λα Z ( xα ) * n α =1 处的Z(x)值（2）求x处的）处的值无偏条件* E[ Z UK ( x )] = E*∑ λα Z ( xα )+ = ∑ λα E* Z ( xα )+ α =1 α =1 n n = ∑ λα m( xα ) = ∑ λα ∑ al f l ( xa ) i =1 k n n k α =1 l =0 = ∑ al *∑ λα f l ( xa )+ l =0 n α =1 m( x) = ∑ al f l ( x) l =0 k 无偏条件∑ a *∑ λα f ( x )+ = ∑ a f ( x) α l =0 l =1 l a l =0 l l k n k ∑ λα f (x ) = f ( x) α =1 l a l n l = 0,1, L k 泛克立格估值的最优条件δ = E* Z 2 E n * UK ? Z ( x)+ = E*∑ λα Z ( xα ) ? Z ( x)+2 2 n α =1 = E,*∑ λα Z ( xα )+ ? 2∑ λα Z ( xα )Z ( x) + * Z ( x)+2 - 2 n α =1 n α =1 = ∑∑ λα λβ E* Z ( xα )Z ( xβ )+ ? 2∑ λα E* Z ( xα )Z ( x)+ + E* Z ( x)+2 α =1 β =1 n n n n α =1 = ∑∑ λα λβ ,C ( xα , xβ ) + E* Z ( xα )+E* Z ( xβ )+- ? 2∑ λα ,C ( xα , x) + E* Z ( xα )+E* Z ( x)+- α =1 β =1 α =1 n + ,C ( x, x) + E 2 * Z ( x)+- = ∑∑ λα λβ C ( xα , xβ ) ? 2∑ λα C ( xα , xβ ) + C ( x, x) + α =1 β =1 n n n n α =1 ,∑∑ λα λβ E* Z ( xα )+E* Z ( xβ )+ ? 2 E* Z ( x)+∑ λα E* Z ( xα )+ + E 2 * Z ( x)+- α =1 β =1 n n n α =1 = ∑∑ λα λβ C ( xα , xβ ) ? 2∑ λα C ( xα , x) + C ( x, x) α =1 β =1 α =1 n n 方差最小条件,∑∑ λα λβ E* Z ( xα )+E* Z ( xβ )+ ? 2 E* Z ( x)+∑ λα E* Z( xα )+ + E 2 * Z ( x)+- α =1 β =1 α =1 n n n = ,E*∑ λα Z ( xα )+- ? 2 E* Z ( x)+,∑ λα E* Z ( xα )+- + ,E* Z ( x)+-2 2 n n α =1 α =1 * * = ,E* ZUK ( x)+-2 ? 2 E* Z ( x)]E[ ZUK ( x)] + {E[ Z ( x)]}2 = {E[ Z ( x)]}2 ? 2 E[ Z ( x)]E[ Z ( x)] + {E[ Z ( x)]}2 =0 泛克立格方程组Q = ∑∑ λα λβ Cov( xα , xβ ) ? 2∑ λα C ( xα , x) + C ( x, x) ? 2∑ ?l * ∑ λα f l ( xα ) ? f l ( x)+ α =1 β =1 α =1 l =0 n n n k n α =1 n k ?Q = 2∑ λβ Cov( xα , xβ ) ? 2C ( xα , x) ? 2∑ ?l f l ( xα ) = 0 (α = 1,2,L n ) ?λα β =1 l =0 n ?Q = ?2* ∑ λαf l ( xα ) ? f l ( x)+ = 0 (l = 0,1,2,L k ) ??l α =1 k ?n ?∑ λβ C ( xα , xβ ) ?∑ ?l f l ( xα ) = C ( xα , x) (α = 1,2,L n) ? β =1 l =0 ? n ? ∑ λα fl ( xα ) ? f l ( x) = 0 (l = 0,1,2,L k ) ? α =1 ? 估计方差k n δ 2 UK = C ( x, x) + ∑ ?l f l ( x) ? ∑ λα C ( xα , x) l =0 α =1 2.3 非平稳区域化变量的变差函数γ R (h) = 1 2 E * R ( x ) + R ( x + h )+ 可以证明,Z(x)的变差函数在局部的范围内约等于剩余的变差函数2γ ( h ) = E,* Z ( x ) ? Z ( x + h )+ - 2 = E * R ( x ) + m ( x ) ? R ( x + h ) ? m( x + h )+ = 2γ R ( h ) + * m ( x ) ? m ( x + h )+ 2 2 2 = E,* R ( x ) ? R ( x + h )] + [ m ( x ) ? m ( x + h )]} 可以证明：E{[m( x) ? m( x + h)][ R ( x) ? R( x + h)]} = 0 当h足够小时：m( x) = m( x + h) R ( x) ? R ( x + h) = Z ( x) ? Z ( x + h) 泛克立格方法的应用泛克立格方法的应用2 协同克立格（Cokriging）协同克立格（）? 协同克立格是综合多个变量信息的估值方法? 协同克立格可以提高估值的准确性，协同克立格可以提高估值的准确性，并降低估计方差2.1 与协同克立格有关的概念? 互协方差函数C j,k (h) = E[Z j (x + h) ? Zk (x)]? mk mj ?互变差函数互变差函数γ j,k (h) = 1 E*Z j (x + h) ? Z j (x)+*Zk (x + h) ? Zk (x)+ 2 2.2 协同克立格方程组的条件Z V k0 = ∑ α∑ k =1 k K nk =1 λα Z α k k 无偏条件* E[ ZVk0 ? ZVk0 ] = E ( ZVk0 ) ? nk0 K ∑ λα α k0 nk0 =1 k0 E ( Zα k0 ) ? ∑ ∑ λα k E (Zα k ) k ≠ k 0 α k =1 K nk = mk0 *1 ? ∑ λα α k0 =1 k0 + ? ∑ mk0 ∑ λα k k ≠ k0 nk α k =1 n k0 ? ? ∑= 1λ α k 0 = 1 α k0 无偏条件? n ? k ? ∑= 1 λ α k = 0 k = 1, 2 .... K , k ≠ k 0 ?α k ? 协同克立格方程组的条件估计方差σ 2 E = E * Z Vk ? Z 0 * V k0 + 2 K nk = C k 0 , k 0 (V k 0 , V k 0 ) ? 2 ∑ ? k =1 ∑C α k =1 k 0 ,k (V k 0 , v α k ) ∑ ∑ α∑ β∑ λ α k =1 k ′=1 k K K nk n′k ′k =1 =1 k λ β ′C k ,k ′( vα , v β ′) k k k 协同克立格方程组∑ β∑ λ β k ′=1 k′K nk′=1 k′C k ′, k ( v β k ′, v α k ) ? ? k = C k 0, k (V k 0 , v α k ) α k = 1, 2 , L , n k ; k , k ′= 1, 2 ,..., K ∑ λα α k0 n k0 =1 k0 =1 = 0 ? k ≠ k0 K ∑ λα α k nk =1 k ∑ K k =1 nk + K 权系数，个个未知数∑ nk 权系数，K个u k =1 2.3 协同克立格应用应用条件? 主变量的信息不足，而次要变量的信息主变量的信息不足，充分? 主变量与次要变量存在一定的相关性? 估值时至少存在一个主变量采样点协同克立格应用实例本讲小结泛克立格方法的原理及应用协同克立格方法的原理及应用克里格法（Kriging）是地统计学的主要内容之一，从统计意义上说，是从变量相关性和变异性出发，在有限区域内对区域化变量的取值进行无偏、最优估计的一种方法；从插值角度讲是对空间分布的数据求线性最优、无偏内插估计一种方法。

克里格，或者说克里金插值Kriging。

法国krige名字来的。

特点是线性，无偏，方差小，适用于空间分析。

所以很适合地质学、气象学、地理学、制图学等。

相对于其他插值方法。

主要缺点：由于他要依次考虑（这也是克里格插值的一般顺序）计算影响范围，考虑各向异性否，选择变异函数模型，计算变异函数值，求解权重系数矩阵，拟合待估计点值，所以反映速度很慢。

（当然也看你算法设计和电脑反应速度了呵呵）。

而那些趋势面法，样条函数法等。

虽然较快，但是毕竟程度和适合用范围都大受限制。

具体对比如下：方法外推能力逼近程度运算能力适用范围距离反比加权法分布均匀时好差快分布均匀最近邻点插值法不高强很快分布均匀三角网线性插值高差慢分布均匀样条函数高强快分布密集时候克里金插值高强慢均可克里格插值又分为：简单，普通，块，对数，指示性，泛，离析克里金插值等。

克里金插值的变异函数球形模型，指数模型，高斯模型，纯块金模型，幂函数模型，迪维生模型等。

以下结合我的绘制等值线（等高线）的程序和高斯迭代解矩阵方程方法以及多元线性回归方法（此两方法实现另补充）说明克里格方法的实现：注：选择变异函数模型为球形模型，选择插值方法为普通克里金，我为了简化问题，考虑为各向同性，变差距离为固定。

int i,j,i0,i1,j0,j1,k,l,m,n,p,h;//循环变量double *r1Matrix;//系数矩阵double *r0Matrix;//已知向量double *langtaMatrix;//待求解向量double *x0;//已知点横坐标double *y0;//已知点纵坐标double * densgridz;//存储每次小方格内的已知值。

double densgridz0;//待求值int N1=0;//统计有多少个已知值double r[71],r0[71];int N[70];for(i=0;i<100;i++){for(j=0;j<100;j++){if(bdataprotected[i*100+j]) continue;//原值点不需要插值//1.遍历所有非保护网格。