地学统计-第五章课件

第五章空间插值与克里格法华中农业大学杨勇8/5/20231第五章空间插值与克里格法华中农业大学杨勇7/31/20对自然各种属性的测量只能得到有限样本点的值，不可能对每个点都进行采样，然而我们总是想知道未测点的值，因此，就需要根据实测得到的离散数据，对未知点进行预测。8/5/20232华中农业大学资源与环境学院对自然各种属性的测量只能得到有限样本点的值，不可能对每个点都算法分析8/5/20233华中农业大学资源与环境学院算法分析7/31/20233华中农业大学资源与环境学院要想获得上图，就必须有各种空间分析手段支持，特别是各种空间插值方法。大多数空间插值（预测）方法，都可被看做是数据的加权平均数，因此，有下面的通用预测公式：如何分配权重是关键问题，多数空间插值方法只考虑到系统的或确定的变异，而没有考虑到其误差，地统计的克里格方法，不但能描述空间变异的分布特征，而且能够表达预测误差。8/5/20234华中农业大学资源与环境学院要想获得上图，就必须有各种空间分析手段支持，特别是各种空间插5.1主要空间插值法简介分类：确定性方法：基于实测数据的相似性程度或平滑程度，利用数学函数进行插值（如逆距离加权法）地统计方法：利用实测数据的统计特性来量化其空间自相关程度，生产插值面并评价预测的不确定性8/5/20235华中农业大学资源与环境学院5.1主要空间插值法简介分类：7/31/20235华中农业5.1主要空间插值法简介分类：整体插值法：利用整个实测数据集来预测局部插值法：在大面积的研究区域上选取较小的空间单元，利用预测点周围的临近样点来进行预测8/5/20236华中农业大学资源与环境学院5.1主要空间插值法简介分类：7/31/20236华中农业5.1主要空间插值法简介分类：精确性插值：实测点的预测值等于采得的实测值不精确插值：实测点的预测值不等于采得的值8/5/20237华中农业大学资源与环境学院5.1主要空间插值法简介分类：7/31/20237华中农业5.1.1空间整体插值法1、全局多项式插值法（趋势面分析法）：即用数学公式表达感兴趣区域上的一种渐变的趋势。平面：曲面：多项式中的参数系数往往用最小二乘法求解。但该方法是不精确的插值方法，很少有实测点刚好在生产的插值面上，而是或高或低于插值面，高低数值相加，之和近似为0。8/5/20238华中农业大学资源与环境学院5.1.1空间整体插值法1、全局多项式插值法（趋势面分析法5.1.1空间整体插值法全局多项式插值法的插值结果往往呈条带状（左图），适合于描述那些呈明显趋势分布的属性，不适合描述那些空间分布波动较大（较破碎，右图）的自然属性8/5/20239华中农业大学资源与环境学院5.1.1空间整体插值法全局多项式插值法的插值结果往往呈条5.1.1空间整体插值法2、变换函数插值法：根据一个或多个空间参量的经验方程进行整体空间插值，这种经验方程称为变换函数。即用与被预测属性相关的其他属性建立回归方程，进行空间预测：b0,b1,b2为回归系数，p1，p2为独立空间变量，z(x)为被预测属性8/5/202310华中农业大学资源与环境学院5.1.1空间整体插值法2、变换函数插值法：根据一个或多个5.1.2空间局部插值法1、泰森多边形插值：由一组连续多边形组成，多边形的边界是由相邻两点直线的垂直平分线组成。特性：（1）每个多边形内仅包含一个离散数据点。（2）在多边形内的任一点k(x，y)同Pi(xi，yi)之间距离总小于它同其它离散点Pj(xj，yj)之间距离。（3）泰森多边形的任意一个顶点必有三条边与它连接，这些边是相邻三个泰森多边形两两拼接的公共边。（4）泰森多边形的任意一个顶点周围存在三个离散点，将其连成三角形后其外接圆的圆心即为该顶点，该三角形称泰森三角形8/5/202311华中农业大学资源与环境学院5.1.2空间局部插值法1、泰森多边形插值：由一组连续多边5.1.2空间局部插值法各泰森多边形内的每一点属性均由各多边形内的已知点确定，若求数据域内任意一点数据属性Z（xi，yi），则需首先判断待求点所落入的多边形，然后再由控制该多边形的已知点Z（x，y）推算得到。8/5/202312华中农业大学资源与环境学院5.1.2空间局部插值法各泰森多边形内的每一点属性均由各多5.1.2空间局部插值法2、三角测量插值法：将采样点用直线与其相邻点连接成三角形，三角形内部包括任何样点，形成一个包括多个倾斜三角板的多面体（TIN）未测点只可能在三角形内或三角形边线上，利用线性插值即可求得缺点是每个预测值只是根据三个实测值得到，且有时会产生突变现象8/5/202313华中农业大学资源与环境学院5.1.2空间局部插值法2、三角测量插值法：将采样点用直线5.1.2空间局部插值法3、逆距离加权法（IDW）：利用被预测区域点周围的实测值来预测未采样点的值，实测点离预测点越近，则对插值的结果影响越大。其中p为实测值对预测值的影响级，若p=0，则每一个权重是一样的，预测值是所有实测值的平均值，当p增加时，相距较远的点的权重迅速减小，2最为常用。由于IDW方法只考虑距离进行权重分配，所以临近实测点的贡献往往很大，而造成空间分布的多点中心现象。8/5/202314华中农业大学资源与环境学院5.1.2空间局部插值法3、逆距离加权法（IDW）：利用被5.1.2空间局部插值法4、局部多项式插值法（移动内插法）：多项式插值法将整个区域考虑成一个平面或曲面，而局部多项式插值法是在划定的领域内（窗口内）用其中的实测数据来拟合不同次数的多项式。8/5/202315华中农业大学资源与环境学院5.1.2空间局部插值法4、局部多项式插值法（移动内插法）5.1.2空间局部插值法5、简单移动平均法：6、样条插值法：8/5/202316华中农业大学资源与环境学院5.1.2空间局部插值法5、简单移动平均法：6、样条插值法5.1.2空间局部插值法7、克里格方法：和IDW一样，也是一种局部估计的加权平均，但是它对各实测点权重的确定是通过半方差分析获取的，可分为线性克里格法和非线性克里格法。（1）普通克里格（6）概率克里格（2）简单克里格（7）贝叶斯克里格（3）泛克里格（8）普通协同克里格（4）指示克里格（5）析取克里格8/5/202317华中农业大学资源与环境学院5.1.2空间局部插值法7、克里格方法：和IDW一样，也是克里格法实质上是利用区域化变量的原始数据和变异函数的结构特点，对未采样点的区域化变量的取值进行线性无偏最优估计的一种方法，从数学角度讲就是一种对空间分布的数据求线性最优无偏内插估计量的一种方法。是根据待估样点有限领域内若干已测定的样点数据，在考虑样点形状、大小和空间相互位置关系，它们与待估样点相互空间位置关系，以及变异函数提供的结构信息之后，对该待估样点进行的一种线性无偏最优估计8/5/202318华中农业大学资源与环境学院克里格法实质上是利用区域化变量的原始数据和变异函数的结构特点5.2线性预测克里格法5.2.1普通克里格法：假定Z(x)是满足本证假设的一个随机过程，该随机过程有n个观测值z(xi)，要预测未采样点x0处的值，则线性预测值Z*(x0)可以表示如下：（5.2.1）Kriging是在使预测无偏并有最小方差的基础上，去确定最优的权重值，满足以下两个条件：（1）无偏性条件（5.2.2）（2）最优条件：（5.2.3）8/5/202319华中农业大学资源与环境学院5.2线性预测克里格法5.2.1普通克里格法：假定Z(x普通克里格法利用式5.2.1取代式5.2.2的左边部分，有：（5.2.4）根据本征假设，上式可以进一步表示为：要使上式成立，则必须满足以下条件：8/5/202320华中农业大学资源与环境学院普通克里格法利用式5.2.1取代式5.2.2的左边部分，有：本征假设条件1：条件2：r(h)称为半方差函数，也叫变异函数本征假设是地统计学中对随机函数的基本假设在研究区域内，区域化变量Z(x)的增量的数学期望对任意x和h存在且等于0在研究区域内，区域化变量的增量[Z(x)-Z(x+h)]的方差对任意x和h存在且平稳8/5/202321华中农业大学资源与环境学院本征假设条件1：r(h)称为半方差函数，也叫变异函数在研究区普通克里格法在本证假设条件下，公式（5.2.3）左边的式子可以表示为：根据方差最小原则，借助拉格朗日乘子，普通克里格的预测方程组为：预测方差为：8/5/202322华中农业大学资源与环境学院普通克里格法在本证假设条件下，公式（5.2.3）左边的式子可普通克里格法克里格公式也可以用矩阵的形式表示，对点状克里格，有：8/5/202323华中农业大学资源与环境学院普通克里格法克里格公式也可以用矩阵的形式表示，对点状克里格，普通克里格法实例8/5/202324华中农业大学资源与环境学院普通克里格法实例7/31/202324华中农业大学资源与普通克里格法实例8/5/202325华中农业大学资源与环境学院普通克里格法实例7/31/202325华中农业大学资源与8/5/202326华中农业大学资源与环境学院7/31/202326华中农业大学资源与环境学院8/5/202327华中农业大学资源与环境学院7/31/202327华中农业大学资源与环境学院8/5/202328华中农业大学资源与环境学院7/31/202328华中农业大学资源与环境学院

？8/5/202329华中农业大学资源与环境学院7/31/202329华中农业大学资源与环境学院对于上图，由于数据构型没有变，因此，新的未知点对克里格矩阵没有影响，对矩阵B有影响，因而值需要重新计算矩阵B即可，这样可以节省计算过程，同时也是为什么特别强调野外实际取样设计应尽可能采用规则格网数据结构的重要原因。8/5/202330华中农业大学资源与环境学院对于上图，由于数据构型没有变，因此，新的未知点对克里格矩阵没新的估计点x0的克里格方程组为：新的估计点的克里格估计值计算公式为：估计方差为：8/5/202331华中农业大学资源与环境学院新的估计点x0的克里格方程组为：7/31/202331华中农不规则网络取样数据实例在研究区域内，区域化变量Z(x)满足二阶平稳和本征假设，其协方差函数和变异函数存在，且均为各向同性的指数模型：8/5/202332华中农业大学资源与环境学院不规则网络取样数据实例在研究区域内，区域化变量Z(x)满足二区域中位置点x0，在其领域内有7个已知不规则分布的样点，每一个样点的坐标，样点的观测值与位置点x0的距离如表所示：8/5/202333华中农业大学资源与环境学院区域中位置点x0，在其领域内有7个已知不规则分布的样点，每一7个已知样点间的距离8/5/202334华中农业大学资源与环境学院7个已知样点间的距离7/31/202334华中农业大学资源8/5/202335华中农业大学资源与环境学院7/31/202335华中农业大学资源与环境学院未知点x0的估计值为：其克里格估计方差为：8/5/202336华中农业大学资源与环境学院未知点x0的估计值为：7/31/202336华中农业大学资块段普通克里格法（BlockKriging）如图，如果估计的不是x0点，而是以x0为中心快段的平均值，这时要用块段普通克里格法，由于与上例数据结构不变，因此，克里格矩阵不变，只是矩阵B变为8/5/202337华中农业大学资源与环境学院块段普通克里格法（BlockKriging）如图，如果估计其中V为块段面积，但积分经常很难解，因此，常常（1）将待估块段V离散化成若干点，计算已知样点与块段V内被离散化的若干个点之间的协方差函数C(h)或变异函数r(h)，最后求出平均数；（2）除了将待估块段V离散化外，如果已知样点也是一个小块段的支撑构成，将这块段也离散化成一只的若干点，计算点与点之间的协方差函数C(h)或变异函数r(h)x0x01x02x03x048/5/202338华中农业大学资源与环境学院其中x0x01x02x03x047/31/202338华中农离散化后，待估块段V中4个离散化点孩子见的平均协方差函数为：离散化是一个近似的方法，显然，离散化越多，得到的平均协方差函数值或变异函数值精度越高，但计算的工作量很大，因此，在离散化时，要考虑精度和工作量的平衡。8/5/202339华中农业大学资源与环境学院离散化后，待估块段V中4个离散化点孩子见的平均协方差函数为：有关克里格法的说明在克里格方程组中，克里格矩阵只取决于已知样点的相对集合特征，而与待估样点或块段无关，因此，只要两个数据构型相同时，其克里格矩阵也相同。计算过程中，只需求一次逆矩阵，就可以得到另一个接的列矩阵，如果两个待估块段或样点的几个特征也相同，那么所得克里格权重系数的解得列矩阵也一定相同，鉴于这方面的原因，空间取样设计时应尽量保证数据构形的系统性和规则性。这样，就可将单独一个克里格方案在整个估计中重复使用。8/5/202340华中农业大学资源与环境学院有关克里格法的说明在克里格方程组中，克里格矩阵只取决于已知样克里格估计的可靠性克里格估计量是一种线性、无偏最优估计量，而且可以给出估计的方差，它在已知样点上的估计等于样点值本身。因而它是真正的空间局部内插法，具有较高的可靠性。许多研究表明，克里格法的精度明显高于多边形法、三角形法、局部平均法和距离倒数平法加权法。8/5/202341华中农业大学资源与环境学院克里格估计的可靠性7/31/202341华中农业大学资源与克里格权重系数的对称性如图，每个样点的权重系数分别为k1,k2…k7，当区域化变量Z(X)具有各向同性结构时，在圈中系数之间有k2=k3,k4=k6.这主要是因为样点x2与x3，x4与x6对待估点x0几何位置是对称的，因而它们之间的克里格权重系数也具有对称性。8/5/202342华中农业大学资源与环境学院克里格权重系数的对称性如图，每个样点的权重系数分别为k1,k丛聚效应降低克里格权重系数在上图中，即使区域化变量Z(x)是各向同性结构的，x1与x5对待估样点x0几何对称，也不能存在k1=k5，因为在样点x5附近还存在样点x4,x6和x7，这三个样点与x5丛聚在一起，而x1是单独的一个样点，计算克里格权重系数的结果使k1>k5。x4,x5和x7的丛聚作用降低了x5对待估样点x0的影响。而x1不存在丛聚效应。因此，在克里格估计中，不会由于一些样点丛聚在一起而增大其权重系数，这也正是克里格法估计的优点。8/5/202343华中农业大学资源与环境学院丛聚效应降低克里格权重系数在上图中，即使区域化变量Z(x)是丛聚效应降低克里格权重系数8/5/202344华中农业大学资源与环境学院丛聚效应降低克里格权重系数7/31/202344华中农业大学屏蔽效应当块金值很小或不存在时，已知样点的克里格权重系数的大小受屏蔽效应影响，如图，已知样点x5虽然与样点x1到待估样点x0的距离相等，但是x1的克里格权重系数k1确大于x5的权重系数k5，这主要是因为样点x5受x4的屏蔽效应影响。8/5/202345华中农业大学资源与环境学院屏蔽效应当块金值很小或不存在时，已知样点的克里格权重系数的大屏蔽效应如图，待估点x0附近有x1,x2,……x12，共12个已知样点，由数据构形可知，x1,x2,x3,x4与x0的几何位置对称且相等，x5,x6…x12与x0的几何位置也对称且相等。内圈的每个已知样点的克里格权重系数是(1-k)/4，外圈的克里格权重系数是k/8，总和为1，因此，内圈的权重系数明显大于外圈样点的权重系数，这是由于内圈样点屏蔽了外圈样点的缘故。8/5/202346华中农业大学资源与环境学院屏蔽效应如图，待估点x0附近有x1,x2,……x12，共12屏蔽效应8/5/202347华中农业大学资源与环境学院屏蔽效应7/31/202347华中农业大学资源与环境学院屏蔽效应屏蔽效应还与块金常数有很大关系，当块金常数增大时，屏蔽效应减弱，当为纯块金效应时，所有样点之间相互独立，协方差函数为0，变异函数等于外延方差，即基台值，则待估样点x0与周围任何已知样点的克里格权重系数均相同，此时屏蔽效应消失，任何一点上的克里格线性无偏最优估计量都是所有样点的算术平均值。8/5/202348华中农业大学资源与环境学院屏蔽效应屏蔽效应还与块金常数有很大关系，当块金常数增大时，屏理论模型对克里格估值的影响数据由8个样点和1个待估点组成，结构分析后，采用球状模型作为理论模型样点位置、属性值及与待估点的距离8/5/202349华中农业大学资源与环境学院理论模型对克里格估值的影响数据由8个样点和1个待估点组成，结理论模型对克里格估值的影响球状模型：8/5/202350华中农业大学资源与环境学院理论模型对克里格估值的影响球状模型：7/31/202350华1）尺度对克里格估值的影响若令r2(h)=0.5r1(h)，球状模型尺度变化后，屏蔽效应减弱，方差较小，精度提高8/5/202351华中农业大学资源与环境学院1）尺度对克里格估值的影响若令r2(h)=0.5r1(h)，（1）尺度对克里格估值的影响r1(h)r2(h)8/5/202352华中农业大学资源与环境学院（1）尺度对克里格估值的影响r1(h)r2(h)7/31/2（2）变程对克里格估计的影响原球状模型变程为40，如果变为20，则新的球状模型为，则内圈样点的克里格权重系数增大，而外圈样点的克里格权重系数减小，增加了屏蔽效应，同时也增大了克里格估计方差，使估计的精度降低。8/5/202353华中农业大学资源与环境学院（2）变程对克里格估计的影响原球状模型变程为40，如果变为2（2）变程对克里格估计的影响8/5/202354华中农业大学资源与环境学院（2）变程对克里格估计的影响7/31/202354华中农业大（3）块金效应对克里格估值的影响由变异函数性质可知，，基台值相当样点之间自相关消失，相互独立时样本方差。对于一个已抽取的空间样本，其方差已确定，因此基台值=常数，当C0增加时，供高要降低，当C0增加到与样点独立时的方差时，C=0，此时基台值就是块金常数C0，因此，块金常数C0的增大，导致样点之间的相关性降低而独立性加大。8/5/202355华中农业大学资源与环境学院（3）块金效应对克里格估值的影响由变异函数性质可知，（3）块金效应对克里格估值的影响若把原来的球状模型块金值扩大一倍，即C0=0.34，基台值不变，则模型的曲线图：8/5/202356华中农业大学资源与环境学院（3）块金效应对克里格估值的影响若把原来的球状模型块金值扩大（3）块金效应对克里格估值的影响块金值变大后，降低了内圈样点的权重系数，增大了外圈样点的权重系数，而克里格方差也增大了近一倍，这说明块金效应增大可使屏蔽效应降低，显然，当C0=0.52时，8个样点的克里格权重系数均为0.125，此时的块金效应相当于纯块金效应。块金值不变时块金值增大1倍8/5/202357华中农业大学资源与环境学院（3）块金效应对克里格估值的影响块金值变大后，降低了内圈样点（4）理论模型的种类对克里格估计的影响如果换成线性有基台值模型：两种模型变异函数曲线（左）和协方差函数曲线（右）对比8/5/202