kriging(克里金方法,克里金插值)#高级教育

《kriging(克里金方法,克里金插值)#高级教育》由会员分享，可在线阅读，更多相关《kriging(克里金方法,克里金插值)#高级教育（107页珍藏版）》请在装配图网上搜索。

1、 第二讲第二讲 克里金方法（克里金方法（Kriging）,是以南非矿业是以南非矿业工程师工程师D.G.Krige(克里格克里格)名字命名的一项名字命名的一项实用空间估计技术，是实用空间估计技术，是地质统计学地质统计学 的重要的重要组成部分，也是地质统计学的核心。组成部分，也是地质统计学的核心。1优选内容主要是为解决矿床储量计算和误差估计问题而主要是为解决矿床储量计算和误差估计问题而发展起来的发展起来的 由法国巴黎国立高等矿业学院由法国巴黎国立高等矿业学院G马特隆教授于马特隆教授于1962年所创立。年所创立。2优选内容H.S.Sichel(1947)D.G.Krige(1951)Kriging法

2、法（克里金法，克立格（克里金法，克立格法）法）:“根据样品空间位置不同、样根据样品空间位置不同、样品间相关程度的不同，对每个样品品间相关程度的不同，对每个样品品位赋予不同的权，进行滑动加权品位赋予不同的权，进行滑动加权平均，以估计中心块段平均品位平均，以估计中心块段平均品位”G.Materon(1962)提出了提出了“地质统计学地质统计学”概念概念 (法文法文Geostatistique)发表了专著发表了专著应用地质统计学论应用地质统计学论。阐明了一整套区域化变量的理论，阐明了一整套区域化变量的理论，为地质统计学奠定了理论基础。为地质统计学奠定了理论基础。区域化变量理论区域化变量理论克里金估计

3、克里金估计随机模拟随机模拟应用统计学方法研究金矿品位应用统计学方法研究金矿品位1977年我国开始引入年我国开始引入 3优选内容 克里金插值方法克里金插值方法niiixzxz10*井眼地震（普通克里金）（应用（应用随机函数随机函数理论）理论）不仅考虑待估点位置与不仅考虑待估点位置与 已知数据位置的相互关已知数据位置的相互关 系，而且还考虑变量的系，而且还考虑变量的 空间相关性。空间相关性。4优选内容 为一个实值变量，可根据概率分布取不同的值。为一个实值变量，可根据概率分布取不同的值。每次取值（观测）结果每次取值（观测）结果z为一个确定的数值，称为为一个确定的数值，称为随机变量随机变量Z的的一个实

4、现。一个实现。P1.随机变量随机变量5优选内容连续变量：连续变量：累积分布函数（cdf）cumulative distribution function)(Pr);(zuZobzuF条件累积分布函数（ccdf）后验 conditional cumulative distribution function)(|)(Pr)(|;(nzuZobnzuF离散变量（类型变量）：离散变量（类型变量）：)(|)(Pr)(|;(nkuZobnkuFZ(u)PP不同的取值方式：估计（estimation）模拟（simulation）6优选内容连续型地质变量连续型地质变量构造深度构造深度砂体厚度砂体厚度有效厚度有

5、效厚度孔隙度孔隙度渗透率渗透率含油饱和度含油饱和度离散型地质变量离散型地质变量（范畴变量）（范畴变量）砂体砂体相相 流动单元流动单元隔夹层隔夹层断层断层类型变量类型变量7优选内容设设离散型随机变量离散型随机变量的所有可能取值为的所有可能取值为 x1，x2，其相应的概率为，其相应的概率为P(=xk)=pk,k=1,2,.随机变量的特征值：随机变量的特征值：(1)数学期望数学期望 是随机变量是随机变量的整体代表性特征数。的整体代表性特征数。则当级数 绝对收敛时，称此级数的和为的数学期望，记为E()，或E。E()=1kkpxk1kkpxk8优选内容设连续型随机变量的可能取值区间为(-，+)，p(x)

6、为其概率密度函数，若无穷积分 绝对收敛，则称它为的数学期望，记为E()。dxxxp)(E()=dxxxp)(数学期望是随机变量的最基本的数字特征，相当于随机变量以其取值概率为权的加权平均数。从矩的角度说，数学期望是的一阶原点矩。对于一组样本：对于一组样本：NzmNii)(19优选内容 为随机变量的离散性特征数。若数学期望E-E()2存在，则称它为的方差，记为D()，或Var()，或2。=222)(E-)()(E-E)(ED 从矩的角度说，方差是的二阶中心矩。(2)方差方差 其简算公式为 D()=E(2)E()2D()=E-E()2方差的平方根为标准差，记为 10优选内容 研究范围内的一组随机变

7、量。研究范围内的一组随机变量。),(研究范围uuZ)(uZ简记为)(|)(,)(Pr)(|,;,(1111nzuZzuZobnzzuuFKKKK 随机场：随机场：当随机函数依赖于多个当随机函数依赖于多个自变量时，称为随机场。自变量时，称为随机场。如具有三个自变量如具有三个自变量(空间空间点的三个直角坐标点的三个直角坐标)的随的随机场机场2.随机函数随机函数条件累积分布函数（ccdf）P11优选内容 二个随机变量二个随机变量，的协方差为二维随机变量的协方差为二维随机变量(，)的二阶混合中心矩的二阶混合中心矩11，记为，记为Cov(，)，或，或,。协方差协方差(Variance)：Cov(,)=,

8、=E-E()-E()其简算公式为其简算公式为 Cov(,)=E()-E()E()随机函数的特征值随机函数的特征值 12优选内容P 任何统计推断（cdf，数学期望等）均要求重复取样。但在储层预测中，一个位置只能有一个样品。同一位置重复取样，得到cdf，不现实13优选内容考虑邻近点，推断待估点 空间一点处的观测值可解释为一个随机变量在该点 处的一个随机实现。空间各点处随机变量的集合构成一个随机函数。区域化变量：能用其空间分布来表征一个自然现象的变量。（将空间位置作为随机函数的自变量）（可以应用随机函数理论解决插值和模拟问题）14优选内容考虑邻近点，推断待估点 -空间统计推断要求平稳假设),;,()

9、,;,(1111KKKKzzhuhuFzzuuF 严格平稳严格平稳);();(zhuFzuF对于单变量而言：可从研究区内所有数据的累积直方图推断而得 （将邻近点当成重复取样点）太强的假设，不符合实际P15优选内容 当区域化变量Z(u)满足下列二个条件时，则称其为二阶平稳或弱平稳：EZ(u)=EZ(u+h)=m(常数)xh 随机函数在空间上的变化没有明显趋势，随机函数在空间上的变化没有明显趋势，围绕围绕m值上下波动。值上下波动。在整个研究区内有在整个研究区内有Z(u)的数学期望存在，的数学期望存在，且等于常数，即：且等于常数，即：二阶平稳二阶平稳16优选内容 在整个研究区内，在整个研究区内，Z(

10、u)的协方差函数存在且平稳的协方差函数存在且平稳 (即只依赖于滞后即只依赖于滞后h，而与，而与u无关无关)，即即 CovZ(u),Z(u+h)=EZ(u)Z(u+h)-EZ(u)EZ(u+h)=EZ(u)Z(u+h)-=C(h)特殊地，当h=0时，上式变为VarZ(u)=C(0),即方差存在且为常数。协方差不依赖于空间绝对位置，而依赖于相对位置协方差不依赖于空间绝对位置，而依赖于相对位置，即具有空间的平稳不变性。即具有空间的平稳不变性。uu+h17优选内容 在整个研究区内有在整个研究区内有 EZ(u)-Z(u+h)=0 本征假设本征假设 当区域化变量Z(u)的增量Z(u)-Z(u+h)满足下列

11、二条件时，称其为满足本征假设或内蕴假设。可出现EZ(u)不存在，但EZ(u)-Z(u+h)存在并为零的情况存在并为零的情况 intrinsic hypotheseEZ(u)可以变化,但EZ(u)-Z(u+h)=0（比二阶平稳更弱的平稳假设）18优选内容 增量增量Z(u)-Z(u+h)的方差函数的方差函数(变差函数,Variogram)存在且平稳存在且平稳(即不依赖于即不依赖于u)，即：，即：VarZ(u)-Z(u+h)=EZ(u)-Z(u+h)2-EZ(u)-Z(u+h)2 =EZ(u)-Z(u+h)2 =2(u,h)=2(h),相当于要求：相当于要求：Z(u)的变差函数存在且平稳。的变差函数

12、存在且平稳。19优选内容例：物理学上的著名的布朗运动是一种呈现出无限离散性的物理现象，其随机函数的理论模型就是维纳-勒维(Wiener-Levy)过程(或随机游走过程)。布朗运动:可出现协方差函数不存在，但变差函数存在的情况。既不能确定验前方差，也不能确定协方差函数。但是其增量却具有有限的方差：VarZ(x)-Z(x+h)=2 =A|h|(其中，A是个常数)，变差函数=|h|，且随着|h|线性地增大。2A)(h20优选内容 若区域化变量若区域化变量Z(x)在整个区域内不满足二阶平在整个区域内不满足二阶平稳稳(或本征假设或本征假设)，但在有限大小的邻域内是二阶平，但在有限大小的邻域内是二阶平稳稳

13、(或本征或本征)的，则称的，则称Z(x)是准二阶平稳的是准二阶平稳的(或准本征或准本征的的)。准二阶平稳假设及准本征假设准二阶平稳假设及准本征假设21优选内容 设 为区域上的一系列观测点，为相应的观测值。区域化变量在 处的值 可采用一个线性组合来估计：nxx,1 nxzxz,10 x0*xzZ*(x0)niiixzxz10*min00*00*0 xZxZVarxZxZE无偏无偏最优最优无偏性和估计方差最小被作为 选取的标准 i-以普通克里金为例22优选内容从本征假设出发,可知 为常数,有 xZE 0*11000mmxZxZExZxZEniiniii可得到关系式:11nii（1）无偏条件）无偏条

14、件Z*(x0)（在搜寻邻域内为常数，不同邻域可以有差别）23优选内容njxZxZEnijj,1,021200*（2）估计方差最小）估计方差最小min200*200*00*xZxZExZxZExZxZE2k应用拉格朗日乘数法求条件极值Z*(x0)24优选内容niijniijinjxxCxxC1011,1进一步推导，可得到n+1阶的线性方程组,即克里金方程组 当随机函数不满足二阶平稳,而满足内蕴（本征）假设时,可用变差函数来表示克里金方程组如下:niijniijinjxxxx1011,1Z*(x0)25优选内容最小的估计方差,即克里金方差可用以下公式求解:niiikxxCxxC1000200102

15、xxxxniiikZ*(x0)26优选内容 变差函数变差函数(或叫或叫变程方差函数变程方差函数，或，或变异函数变异函数)是是地质统计学所特有的基本工具。它既能描述区域化地质统计学所特有的基本工具。它既能描述区域化变量的空间结构性变化，又能描述其随机性变化。变量的空间结构性变化，又能描述其随机性变化。跃迁现象1.变差函数的概念与参数变差函数的概念与参数 27优选内容),(hx 假设空间点假设空间点x只在一维的只在一维的x轴上变化，则将区域化轴上变化，则将区域化变量变量Z(x)在在x，x+h两点处的两点处的值之差值之差的方差之半定义的方差之半定义为为Z(x)在在x轴方向上的变差函数，记为轴方向上的

16、变差函数，记为一维情况下的定义：一维情况下的定义：VarZ(x)-Z(x+h)EZ(x)-Z(x+h)2-EZ(x)-Z(x+h)2),(hx21=21半变差函数(或半变异函数)28优选内容 在在二阶平稳假设，或作本征假设二阶平稳假设，或作本征假设，此时：，此时：地质统计学中最常用的基本公式之一。EZ(x)-Z(x+h)=0hVarZ(x)-Z(x+h)EZ(x)-Z(x+h)2-EZ(x)-Z(x+h)2),(hx21=21EZ(x)-Z(x+h)2 ),(hx21=则：29优选内容)()0()(hCCh（二阶平稳假设条件下边查函数与写防查的关系）30优选内容变程变程(Range)：指区域化

17、变量在空间上具有相关性的指区域化变量在空间上具有相关性的范围。在变程范围之内，数据具有相关性；而在变范围。在变程范围之内，数据具有相关性；而在变程之外，数据之间互不相关，即在变程以外的观测程之外，数据之间互不相关，即在变程以外的观测值不对估计结果产生影响。值不对估计结果产生影响。31优选内容具不同变程具不同变程的克里金插的克里金插值图象值图象32优选内容块金值块金值(Nugget)：变差函数如果在原点间断，在地质统计学中称：变差函数如果在原点间断，在地质统计学中称为为“块金效应块金效应”，表现为在很短的距离内有较大的空间变异性，表现为在很短的距离内有较大的空间变异性，无论无论h多小，两个随机变

18、量都不相关多小，两个随机变量都不相关。它可以由测量误差引起，。它可以由测量误差引起，也可以来自矿化现象的微观变异性。在数学上，块金值也可以来自矿化现象的微观变异性。在数学上，块金值c0相当于相当于变量纯随机性的部分。变量纯随机性的部分。33优选内容 如果品位完全是典型的随机变量，则不论如果品位完全是典型的随机变量，则不论观测尺度大小，所得到的实验变差函数曲线总观测尺度大小，所得到的实验变差函数曲线总是接近于纯块金效应模型。是接近于纯块金效应模型。当采样网格过大时，将掩盖小尺度的结构，当采样网格过大时，将掩盖小尺度的结构，而将采样尺度内的变化均视为块金常数。这种而将采样尺度内的变化均视为块金常数

19、。这种现象即为块金效应的尺度效应。现象即为块金效应的尺度效应。块金效应的尺度效应块金效应的尺度效应12111333334优选内容基台值基台值(Sill)：代表变量在空间上的总变异性大小。即为变代表变量在空间上的总变异性大小。即为变差函数在差函数在h大于变程时的值，为大于变程时的值，为块金值块金值c0和和拱高拱高cc之和。之和。拱高拱高为在取得有效数据的尺度上，可观测得到的变异性幅为在取得有效数据的尺度上，可观测得到的变异性幅度大小。当块金值等于度大小。当块金值等于0时，基台值即为拱高。时，基台值即为拱高。=C(0)C(h)(h35优选内容几何各向异性：几何各向异性：变差函数变差函数在空间各个方

20、向上的在空间各个方向上的变程变程不同不同，但，但基台值不变基台值不变（即（即变化程度相等）。这种情变化程度相等）。这种情况能用一个简单的几何坐况能用一个简单的几何坐标变换将各向异性结构变标变换将各向异性结构变换为各向同性结构。换为各向同性结构。带状各向异性：带状各向异性：不同方向不同方向的变差函数具有的变差函数具有不同的基不同的基台值台值，其中，其中变程可以不同，变程可以不同，也可以相同也可以相同。这种情况不。这种情况不能通过坐标的线性变换转能通过坐标的线性变换转化为各向同性，因而结构化为各向同性，因而结构套合是比较复杂的。套合是比较复杂的。地质变量相关性的各向异性地质变量相关性的各向异性12

21、1113333(2)36优选内容2.变差函数的理论模型变差函数的理论模型设Z(x)为满足本征假设的区域化变量，则常见的理论变差函数有以下几类：球状模型球状模型指数模型指数模型高斯模型高斯模型幂函数模型幂函数模型空洞效应模型空洞效应模型37优选内容 接近原点处，变差函接近原点处，变差函 数呈线性形状，在变数呈线性形状，在变 程处达到基台值。程处达到基台值。原点处变差函数的切原点处变差函数的切 线在变程的线在变程的2/3处与处与 基台值相交。基台值相交。ahcahahahchahSphch,2123003球状模型：球状模型：c为基台值，为基台值，a为变程，为变程，h为滞后距。为滞后距。38优选内容

22、指数模型：指数模型：ahcahExpch3exp1 变差函数渐近地逼近变差函数渐近地逼近 基台值。基台值。在实际变程处，变差在实际变程处，变差 函数为函数为0.95c。模型在原点处为直线。模型在原点处为直线。39优选内容高斯模型：高斯模型：223exp1ahch 变差函数渐近地逼近变差函数渐近地逼近 基台值。基台值。在实际变程处，变差函在实际变程处，变差函 数为数为0.95c。模型在原点处为抛物线。模型在原点处为抛物线。40优选内容幂函数模型：幂函数模型：hch.幂函数模型为一种无基幂函数模型为一种无基台值的变差函数模型。这台值的变差函数模型。这是一种特殊的模型。是一种特殊的模型。当当=1时，

23、变差函数为一时，变差函数为一直线，即为线性模型，这直线，即为线性模型，这一模型即为著名的一模型即为著名的布朗运布朗运动（随机行走过程）动（随机行走过程）的变的变差函数模型；差函数模型；当当 1时，变差函数为抛时，变差函数为抛物线形状，为物线形状，为分数布朗运分数布朗运动动(fBm)的变差函数模型。的变差函数模型。布朗运动布朗运动分数布朗运动分数布朗运动分数布朗运动分数布朗运动 h2111h41优选内容空洞效应模型空洞效应模型(Hole Effect)：2cosexp1.bhahch 变差函数并非单调增加，变差函数并非单调增加，而显示出一定周期性的而显示出一定周期性的 波动。波动。模型可以有基台

24、值，也模型可以有基台值，也 可以无基台值；可以有可以无基台值；可以有 块金值，也可以无块金块金值，也可以无块金 值。值。空洞效应在地质上多沿空洞效应在地质上多沿 垂向上出现，如富矿层垂向上出现，如富矿层 与贫矿层互层、砂岩与与贫矿层互层、砂岩与 泥岩频繁薄互层等等。泥岩频繁薄互层等等。（b为富矿化带重复距离）)(hh42优选内容 通过区域化变量的空间观测值来通过区域化变量的空间观测值来构建相应的变构建相应的变差函数模型差函数模型,以表征该变量的主要结构特征。以表征该变量的主要结构特征。(求变求变差差)(1)数据准备数据准备 区域化变量的选取区域化变量的选取、数据质量检查及校正数据质量检查及校正

25、、数据的变换数据的变换（如对渗透率进行对数变换）、（如对渗透率进行对数变换）、数据的统计数据的统计（如分相对储层参数计算平均值、（如分相对储层参数计算平均值、方差，作直方图、相关散点图等）、方差，作直方图、相关散点图等）、丛聚数据的解串丛聚数据的解串等。等。3.区域化变量的区域化变量的43优选内容(2)(2)实验变差函数的计算实验变差函数的计算 实验变差函数是指应用观测值计算的变差函实验变差函数是指应用观测值计算的变差函数。对于不同的滞后距数。对于不同的滞后距h h，可算出相应的实验变，可算出相应的实验变差函数差函数。)(*h=N(h)1i2iih)Z(x-)Z(xN(h)21一维实验变差函数

26、的计算公式(i=1,N(h)Z(xi)-Z(xi+h)2的算术平均值一半即为一个h的变差函数值44优选内容对不同的滞后h，进行计算，得出各个h的变差函数值)(*h=N(h)1i2iih)Z(x-)Z(xN(h)21h3h5hh45优选内容设Z(x)为一维区域化变量，满足本征假设，又已知Z(1)=2，Z(2)=4，Z(3)=3，Z(4)=1，Z(5)=5，Z(6)=3，Z(7)=6，Z(8)=4，)1(*)2(*)3(*例：例：试求：试求：)(*h=N(h)1i2iih)Z(x-)Z(xN(h)2146优选内容721)1(*)2(*)3(*=22+12+22+42+22+32+22=1442=3

27、.0062112+32+22+22+12+12=1220=1.6752112+12+02+52+12=1028=2.8047优选内容2D情况情况（1）分不同方向，进行1D变差函数计算3D情况情况：增加垂向方向（2）确定主变程方向 次变程方向角度容限步长容限h3h5hh四方向试算（考虑主变程方向的 走向、倾向和倾角）48优选内容(3)(3)理论变差函数的最优拟合与结构套合理论变差函数的最优拟合与结构套合 选择合适的理论变差函数模型，同时还需进选择合适的理论变差函数模型，同时还需进行结构套合，从而得到一条反映不同层次（或行结构套合，从而得到一条反映不同层次（或不同空间规模）结构的、统一的、最优的变

28、差不同空间规模）结构的、统一的、最优的变差函数曲线。函数曲线。球状模型球状模型指数模型指数模型高斯模型高斯模型幂函数模型幂函数模型空洞效应模型空洞效应模型49优选内容 复杂的区域化变量往往包含各种尺度上的多层次、多方向的变化性，反映在变差函数上即为多层次结构。将不同结构组合为统一结构的过程称为“结构套合”结构套合结构套合各层次套合各层次套合例如，对于200米宽的河道，在h50m的观测尺度上可以将其与河道间的变化性区分出来，但却无法区分层理和矿物成分的变化性(即无法找出更细微的结构来)，它们在50m尺度得到的结构上只能作为“块金效应”出现。若观测尺度为500米，河道的变化也只能作为“块金效应”。

29、121113333大尺度的变化性总是包含着小尺度的变化性，但却不能从大尺度的变化性中区分出小尺度的变化性。50优选内容)()()()(210rrrr)(0r=。0,0,00rCr)(1r=1 1131311 ar ,Car0 ),2123(ararC)(2r 2 2232322 ar ,Car0 ),r21-r23(aaC=代表微观变化性的变程极小的球状模型，可近似地看作纯块金效应型 球状模型，没有块金常数，基台值为C1，变程为a1，反映了小规模范围的变化 球状模型，没有块金常数，基台值为C2，变程较大，为a2，反映了大规模范围的变化 可以用反映各种不同尺度变化性的多个变差函数之和来表示一个套

30、合结构。（各层次理论模型可以不一样）)(ri可以是不同模型的变差函数51优选内容其中 21aa 则套合结构的表达式为)(r2210213232210133221122110a r ,CCCa ),2123(0 ,r)(21)(230,03raararCCCaraCaCraCaCCr=。0,0,00rCr 1 1131311 ar ,Car0 ),2123(ararC 2 2232322 ar ,Car0 ),r21-r23(aaC)(0r)(1r)(2r=52优选内容对于几何各向异性，先根据异向比压缩 距离轴，使之成为各向同性的模型；对于带状各向异性，运用模型叠加的方法加以处理。先用压缩距离轴

31、的办法，使其变程变为相同，然后再把具有相同变程的两个球状模型叠加起来，构成一个新的球状模型 各方向套合各方向套合（将各向异性套合为各向同性，以便于 在克里金估计时，不同方向均可用统一 的结构模型计算实际的变差函数值）53优选内容(4)(4)变差函数参数的最优性检验：变差函数参数的最优性检验：变差函数是否符合实际，应该进行检验。变差函数是否符合实际，应该进行检验。一种实用的检验方法为一种实用的检验方法为“交叉验证法交叉验证法”（Cross-validationCross-validation），检验标准是在各实测），检验标准是在各实测点，根据周围点计算的点，根据周围点计算的克里金估计值与该实测克

32、里金估计值与该实测值的误差平方值的误差平方平均最小。平均最小。估计误差的平方估计误差的平方与与克里金估计方差克里金估计方差之比越接之比越接近近1 1，则说明变差函数与实际的符合程度越高。，则说明变差函数与实际的符合程度越高。实际上，这种方法在检验变差函数的同时，也实际上，这种方法在检验变差函数的同时，也在检验所使用的克里金估计方法的适用性。在检验所使用的克里金估计方法的适用性。Z*(x0)54优选内容niijniijinjxxCxxC1011,1（以普通克里金为例）i求取变差函数（或协方差）；求取变差函数（或协方差）；解克里金方程组解克里金方程组55优选内容 设有一个油藏，在平面上S1，S2，

33、S3，S4处有四个井点，其孔隙度值分别为Z1，Z2，Z3，Z4。据此估计S0点处的孔隙度值Z0 设孔隙度Z(x)是二阶平稳的。其在平面上的二维变差函数是一个各向同性的球状模型，其参数为：块金值C02，变程a200，拱高C20，即：实例实例200,h 22,200,h0 ),(200)h21-200h2320(20,h 0,(h)3356优选内容Z0的估计量为 41iii*0ZZ普通克里金方程组的矩阵形式为 K =M2 0 1 1 1 1 1 C C C C 1 C C C C1 C C C C1 C C C C K,444342413433323124232221141312114321 ,1

34、 CCCCM04030201221MKniijniijinjxxCxxC1011,1(求解）57优选内容 0 1 1 1 1 1 C C C C 1 C C C C1 C C C C1 C C C C K,444342413433323124232221141312114321)()0()(hChC求解求解:Cij ,1 CCCCM040302012C11C12C01试求200,h 22,200,h0 ),(200)h21-200h2320(20,h 0,(h)33?58优选内容 C11=C22=C33=C44=C(0)=2 =C0+C=22,由于C(h)=C(0)-(h)=22-(h)当ij

35、时，Cij=C(|Si-Sj|)=22-(|Si-Sj|).于是，C12=C21=C04=22-)250(,84.9)200250(21)20025023(202223,22.1)50150(22223113CC,98.4)50100(2222024114CCC,32.2)100100(22223223CC,28.0)100150(22224224CC0)50200(22224334CC,66.12)50(2201C,72.1)150(2203C 59优选内容将以上数值代入普通克里金方程组解的矩阵形式中，得 19.841.724.9812.66 0 1 1 1 1 1 22 0 0.28 4.

36、98 1 0 22 2.32 1.221 0.28 2.32 22 9.841 4.98 1.22 9.84 22 14321经计算得:=0.5182,=0.0220,=0.0886,=0.3712。1234Z00.5182Z1+0.0220Z2+0.0886Z3+0.3712Z421MK60优选内容 搜索邻域搜索邻域注意注意1：搜索邻域中的数据点才参加估计节省CPU和内存局域平稳 搜索椭圆或椭球的选择方法与选择变差函数椭圆或椭球相同。61优选内容注意注意2：参与计算的数据点不能太多，否则计算太慢 一般软件中都内置或可选最大的 数据点数目（与待估点最近的数据点），如10。注意注意3：防止数据丛

37、聚带来的数据代表性不强井眼井眼垂向数据太密，若待估点与该井近，则可能忽视邻井数据八分搜寻，保证各象限均有代表数据62优选内容 若搜寻范围无数据，则应用边际概率。若搜寻范围无数据，则应用边际概率。63优选内容x0 简单克里金简单克里金(SK)(SK)普通克里金普通克里金(OK)(OK)泛克里金泛克里金(UK)(UK)协同克里金协同克里金(CK)(CK)贝叶斯克里金（贝叶斯克里金（BKBK）指示克里金指示克里金(IK)(IK)64优选内容所有克里金估计都应用线性回归算法，形式为：m为期望)()()(1*umuZumZnSK求取权系数的克里金方程组的非平稳形式nnuuCuuCu1),2,1(),()

38、,()(求(n+1)个m(u),求(n+1)(n+1)个C(u,u)65优选内容二阶平稳假设二阶平稳假设EZ(u)=EZ(u+h)=m(常数)C(u,u+h)=C(h)nnSKumuZuuZ11*)(1)()()(简单克里金估计的平稳形式：)()()(1*umuZumZnSKEZ(u)=EZ(u+h)=m(常数)66优选内容应用条件：应用条件：随机函数二阶平稳随机函数二阶平稳 随机函数的期望值 m为常数并已知已知不能用于具有局部趋势的情况nnuuCuuCu1),2,1(),()()(简单克里金方程组的平稳形式：nnuuCuuCu1),2,1(),(),()(C(u,u+h)=C(h)（C与位置

39、有关）（C与位置无关）67优选内容 nuZuuZ1*)(nnunuuCuuuCu111)(,1)()(68优选内容应用要求：应用要求：随机函数二阶平稳或符合内蕴假设随机函数二阶平稳或符合内蕴假设 随机函数的期望值 m在搜寻邻域内稳定但未知 协方差平稳 与简单克里金相比，普通克里金相当于在每一与简单克里金相比，普通克里金相当于在每一个位置个位置u，重新估计，重新估计 m。由于普通克里金估计常使用滑动数据邻域，由于普通克里金估计常使用滑动数据邻域，相当于均值相当于均值m随位置可变，即随位置可变，即Z*(u)，此时，实际，此时，实际上是一种非平稳算法，对应于变化的均值和平稳上是一种非平稳算法，对应于

40、变化的均值和平稳的协方差。的协方差。69优选内容 u uu umZE非平稳随机函数的漂移函数非平稳随机函数的漂移函数(drift),简称为漂移或趋势简称为漂移或趋势)()()(uuuRmZ随机函数随机函数=趋势趋势+残差残差区域化变量Z(X)是非平稳的，即EZ(x)=m(x)Kriging with a trend model (KT)具有趋势的克里金具有趋势的克里金70优选内容ma fkkkK()()uu0用光滑的确定性函数来模拟，或用拟合方法 趋势函数趋势函数一维的线性趋势 maa x()u 01二维的二次趋势:maa xa ya xa ya xy()u 0123242571优选内容R()

41、u用均值为0、协方差函数为 的平稳随机函数来模拟。CR()hZZKTKTn*()()()()uuu1泛克里金估计值泛克里金估计值:残差残差72优选内容()()()()()()(),()()(),KTRnkkkKRKTkknCfCnffkKuuuuuuuuuu1011201 为权值()KTk()u是与(K+1)个权值的限制条件相对应的(K+1)个拉格朗日参数泛克里金方程组泛克里金方程组 CR()h为残差协方差函数 73优选内容)()()(10uuuYaamZEZZKTKTn*()()()()uuu1估计值估计值 当K=1时，线性趋势函数为ma fkkkK()()uu0趋势函数可理解为二级变量74

42、优选内容 （1）外部变量必须在空间光滑）外部变量必须在空间光滑 地变化，否则可能导致地变化，否则可能导致KT 线性系统不稳定；线性系统不稳定；（2）在主变量的所有数据点）在主变量的所有数据点u 处和待估计的处和待估计的 位置位置u处，外部变量都必须是已知的。处，外部变量都必须是已知的。nKTnKTRnRKTYYnCYC1)(1)(101)()()()(1)(,1)()()()()()(uuuuuuuuuuuu克里金方程组：可理解为地震 数据（如深度）（K=0时，？）75优选内容 利用几个变量之间的空间相关性，对其中的一个或几利用几个变量之间的空间相关性，对其中的一个或几个变量进行空间估计，尤其

43、适用于被估计变量的观察数据个变量进行空间估计，尤其适用于被估计变量的观察数据较少的情况较少的情况。mjjjniiiyxZ110*协同克里金估计值（初始变量和二级变量）协同克里金估计值（初始变量和二级变量）-随机变量在位置0处的估计值；-初始变量的n个样本数据；-二级变量的m个样本数据；-需要确定的协同克里金加权系数。0*Znxx,1myy,1naa,1及m,176优选内容01,2,1,2,1,1102110111mjjniijmjjiinijiijmjjiinijiimjyxCyyCyxCnixxCxyCxxC协同克里金方程组协同克里金方程组传统普通协克里金传统普通协克里金77优选内容标准化普

44、通协克里金标准化普通协克里金)(110*YXmmyxZmjjjniii111mjjniimX=Ex(u)mY=Ey(u)78优选内容 为协同克里金的简化形式，即如果二级变量密为协同克里金的简化形式，即如果二级变量密集取样时，只保留与估计点同位的二级变量。集取样时，只保留与估计点同位的二级变量。)()()()()(1uYuuZuuZjniii对应的协同克里金方程组只要求知道Z-协方差函数以及Z-Y 互协方差函数CZ（h）)()(hChCZZY)0(/)0()0(ZYZYCCP（同位两种数据的 相关系数）（方差函数）同位协同克里金同位协同克里金 Collocated Cokriging79优选内容

45、H.Omre在（1987）把线性贝叶斯理论用于克里金估计技术，提出了贝叶斯克里金估计技术。他构想了一个模型，把用于空间估计的数据分为两类：观察数据：是指那些精度比较高，但数量比较少的数据 猜测数据：是指那些精度比较低，但分布广泛的数据 在观测数据比较多的地方，估计结果主要受观测数据的影响；在观测数据比较少的地方，则主要受猜测数据的影响。显然，井数据和地震数据的关系符合贝叶斯估计中观测数据和猜测数据的关系。80优选内容 设Z(x)，xA，是观察数据的区域化变量。设M(x)，xA，是猜测数据的区域化变量。Z*(x0)=EZ(x0)=a0+M(x0)EM(x)M(x)，xAM(x)是对Z(x)的一种

46、猜测，误差为a0 x0a0？81优选内容 设已得到设已得到Z(x)，xA的一组的一组(N个个)观察值观察值 Z(xi);i=1,2,N。定义一个新的随机函数：定义一个新的随机函数：ZT(x)Z(x)-M(x)，xAZT(xi)=Z(xi)-M(x)，i=1,2,NZ(x0)的贝叶斯克里金估计量为的贝叶斯克里金估计量为NiMiiBKxxZxZxZT100*0*)()()()(x0对这个N个观察值有（相当于误差a0）（误差的随机函数）82优选内容基于无偏性和估计方差最小两个条件：基于无偏性和估计方差最小两个条件：min00*00*0 xZxZVarxZxZE利用拉格朗日乘数法可得到贝叶斯克里金方程

47、组：利用拉格朗日乘数法可得到贝叶斯克里金方程组：1,00|1|jjiMiMZjjiMjiMZjaxxxxxxxxaNj,2,1 Z(x,x)=ZM(x-x)+M(x,x)83优选内容将数据按照不同的门槛值编码为将数据按照不同的门槛值编码为1或或0的过程。的过程。对于模拟目标区内的对于模拟目标区内的每一类相，当它出现于某每一类相，当它出现于某一位置时，指示变量为一位置时，指示变量为1，否则为否则为0。A(100)B(010)A(100)C(001)类型变量的指示变换：类型变量的指示变换：01ui 变量 u 属于范畴A 其它指示变换指示变换1982年由AGJournel(儒尔奈耳)教授提出 84优

48、选内容(00111)(00001)(01111)(00011)首先将连续变量截断首先将连续变量截断为类型变量，然后进为类型变量，然后进行指示变换。行指示变换。如：z=10，15，20，25,30 zxzzxzzui01;连续变量的指示变换连续变量的指示变换 85优选内容设沿空间某一方向，在间距为h的5对样品点处观测了Z(x)及Z(x+h)的值 (=1,2，5)。1 2 3 4 5 Z(x)0 2 3 6 9 Z(x+h)1 6 7 8 8 设指示XI(x;z)=设指示Y=I(x+h;z)=z)Z(,0z)Z(,1当当zh)Z(x,0zh)Z(x,1当当86优选内容假定只选定了5个门限值：0，2

49、，3，6，9 XI(x;z)Y=I(x+h;z)当 z=当 z=0 2 3 6 9 0 2 3 6 9 1 1 1 1 1 1 0 1 1 1 1 2 0 1 1 1 1 0 0 0 1 1 3 0 0 1 1 1 0 0 0 0 1 4 0 0 0 1 1 0 0 0 0 1 5 0 0 0 0 1 0 0 0 0 1 1 2 3 4 5 Z(x)0 2 3 6 9 Z(x+h)1 6 7 8 8 z)Z(,0z)Z(,1当当zh)Z(x,0zh)Z(x,1当当87优选内容 指示指示(函数函数)的数学期望的数学期望 当x固定时，若z再给定，则I(x；z)就是个随机变量，就有数学期望：EI(x

50、；z)1PI(x；z)=1+0PI(x；z)=0 =PI(x；z)=1=P =F(x；z)，zZ(x)z+在x点处区域化变量Z(x)的先验分布函数F(x；z)，z+就是x点处指示函数的数学期望 EI(x；z)i*(x；z)=),();(11zFzxinAAAnn88优选内容(1)z=3时，设指示随机变量X I (x;3)E(X)=EI(x;3)=51)(6.053)3;(51xmxiVar(X)=Var I(x;3)=mx(1-mx)=0.60.4=0.24=2xXI(x;z)Y=I(x+h;z)当 z=当 z=0 2 3 6 9 0 2 3 6 9 1 1 1 1 1 1 0 1 1 1 1

51、 2 0 1 1 1 1 0 0 0 1 1 3 0 0 1 1 1 0 0 0 0 1 4 0 0 0 1 1 0 0 0 0 1 5 0 0 0 0 1 0 0 0 0 1 试求指示随机变量的数学期望和方差：试求指示随机变量的数学期望和方差：89优选内容(2)z=6时，设指示随机变量YI(x+h;6)E(Y)=EI(x+h;6)=51)(4.052)6;(51YhmxiVar(Y)=Var I(x+h;6)=mY(1-mY)=0.40.6=0.24=2YXI(x;z)Y=I(x+h;z)当 z=当 z=0 2 3 6 9 0 2 3 6 9 1 1 1 1 1 1 0 1 1 1 1 2

52、0 1 1 1 1 0 0 0 1 1 3 0 0 1 1 1 0 0 0 0 1 4 0 0 0 1 1 0 0 0 0 1 5 0 0 0 0 1 0 0 0 0 1 90优选内容设设Z(X)是个一维区域化变量，在等间距的是个一维区域化变量，在等间距的10个点个点 处有处有10个个观测值：观测值：Z(1)=3，Z(2)=5，Z(3)=6，Z(4)=2，Z(5)=7，Z(6)=1，Z(7)=4，Z(8)=8，Z(9)=9，Z(10)=7，设门限值，设门限值Z分别等于分别等于2,3,4,5,6,7,8时，求指示变差函数的估计值时，求指示变差函数的估计值 (h;z)，h=1，2，3，4，5。计算

53、计算*I91优选内容)5;()5;(212hxIxIEI(h;5)=*I)(12)5;()5;()(21hnhxixihn(h;5)=首先，计算i(x;5)：i(x1;5)1,i(x2;5)1,i(x3;5)0,i(x4;5)=1,i(x5;5)0,i(x6;5)=1,i(x7;5)1,i(x8;5)0,i(x9;5)0,i(x10;5)0,*I;2778.0185)001011110(921222222222(1;5)=zxzzxzzxia01;92优选内容 1 2 3 4 5 6 7 8 9 10 i(x;2)i(x;3)i(x;4)i(x;5)i(x;6)i(x;7)i(x;8)h z

54、2 3 4 5 6 7 8 1 2 3 4 5 列表计算各指示值i(x;z)根据i(x;z)值算出不同z值的 (h;z)值*IZ(X)356271489793优选内容h z 2 3 4 5 6 7 8 1 0.2222 0.2778 0.2778 0.2778 0.1667 0.1111 0.1111 2 0.1250 0.1875 0.3125 0.2500 0.2500 0.1875 0.0625 3 0.2857 0.2143 0.2143 0.2857 0.2143 0.1429 0.0714 4 0.2500 0.3333 0.4167 0.3333 0.2500 0.1667 0.

55、0833 5 0.2000 0.1000 0.2000 0.1000 0.2000 0.2000 0.1000 根据i(x;z)值算出的不同z值的 (h;z)值*I 计算的各指示值i(x;z)1 2 3 4 5 6 7 8 9 10 i(x;2)0 0 0 1 0 1 0 0 0 0 i(x;3)1 0 0 1 0 1 0 0 0 0 i(x;4)1 0 0 1 0 1 1 0 0 0 i(x;5)1 1 0 1 0 1 1 0 0 0 i(x;6)1 1 1 1 0 1 1 0 0 0 i(x;7)1 1 1 1 1 1 1 0 0 1 i(x;8)1 1 1 1 1 1 1 1 0 1 9

56、4优选内容指示克里金指示克里金 线性估计量线性估计量 i*(x;z)=nzxizx1);();(其中，x(=1,2,，n)点处的Z(x)值已知，(x;z)(=1,2,，n)为IK权系数(00111)(00001)(01111)(00011)（普通指示克里金）95优选内容nnIIzxnzxxCzzxxCzx111);(.,2,1),;()();();(nnIIzxnzxxzzxxzx111);(.,2,1),;()();();(CI(h;z)为指示协方差，为指示变差函数);(zhI是拉格朗日乘数)(z有多少个门限值z，就有多少个IK方程组 IK方程组：方程组：96优选内容从IK方程组中解出),2

57、,1)(;(nzxi*(x;z)=求出i(x;z)的估计值 nzxizx1);();(i*(x;z)实际上是I(x;z)的条件数学期望EI(x;z)|(n)的估计值。EI(x;z)|(n)=0PI(x;z)=0|(n)+1PI(x;z)=1|(n)=PI(x;z)=1|(n)=PZ(x）z)|(n)=Fx;z|(n)为后验概率分布函数，或条件概率分布函数，等于在已知n个信息样品的条件下，概率P Z(x)z的大小。i*(x;z)=EI(x;z)|(n)即：i*(x;z)=EI(x;z)|(n)又：=PZ(x）z)|(n)97优选内容(00111)(00001)(01111)(00011)Z*(x

58、)=KkkxxZnkxP1*)()()()(各类的局部均值各类的后验概率估计的指示值98优选内容Sunbeson金矿，有3500个样品，边界品位值分别取z=0.009，0.015，0.020，0.033克/吨，用这4个边界品位可将金矿化分为5类：(1)(2)(3)(4)(5)0.033 应用IK法估计出以下的条件概率：P*x(1)|(n)=18.2%,P*x(2)|(n)=6.5%,P*x(3)|(n)=0%,P*x(4)|(n)=57.8%,P*x(5)|(n)=17.5%,各概率的总和为1 例：例：99优选内容各类的各类的(局部局部)均值为：均值为：Z(x)|x(1)*=0.004g/t,

59、Z(x)|x(2)*=0.0125g/t.Z(x)|x(3)*=0.020g/t,Z(x)|x(4)*=0.0276g/t.Z(x)|x(5)*=0.125g/t,按按IK法可得总平均品位的估计值为法可得总平均品位的估计值为 Z*(x)=51*kp x(k)(n)Z(x)x(k)*=0.0393683g/t.100优选内容A(100)B(010)A(100)C(001)离散变量的指示克里金估计离散变量的指示克里金估计分别对各类型进行指示克里分别对各类型进行指示克里金估计，得出各类型的概率金估计，得出各类型的概率估计。估计。00.10.20.30.40.50.60.70.8ABC相P（最大概率的

60、类型 即为待估点的类型）101优选内容 作为一种非参数统计方法，指示克里金估计作为一种非参数统计方法，指示克里金估计方法在处理特高值和特低值的分布方面，具有明方法在处理特高值和特低值的分布方面，具有明显的优势。显的优势。特色一：特色一：(00111)(00001)(01111)(00011)离散变量（类型变量）连续变量102优选内容 应用应用贝叶斯理论贝叶斯理论，可综合各种软信息，可综合各种软信息（与硬信息一起）进行指示克里金估计。（与硬信息一起）进行指示克里金估计。（如地震信息、试井解释、地质推理和解释等）硬信息：能准确标定所研究地质 变量取值的信息。软信息：不能准确标定所研究地质 变量取值

61、，而只能提供其 概率分布的信息。特色二：特色二：103优选内容00.20.40.60.81ABC相P 使用贝叶斯理论，将先验的局部累计分布使用贝叶斯理论，将先验的局部累计分布函数更新为后验的条件局部累计分布函数。函数更新为后验的条件局部累计分布函数。*)()(PrnnzxZob);();();();()(110zxizxzxizxzFnn整体整体先验概率先验概率硬指示数据硬指示数据软指示数据软指示数据A(100)B(010)A(100)C(001)104优选内容l l 估计的无偏性估计的无偏性 克里金方法的优点克里金方法的优点l l 反映了变量的空间结构性反映了变量的空间结构性l l 能得到估

62、计精度能得到估计精度 105优选内容 （1）克里金插值为局部估计方法，）克里金插值为局部估计方法，对估计值的对估计值的整体空间相关性考虑不够整体空间相关性考虑不够，它保证了数据的估计局它保证了数据的估计局部最优，却不能保证数据的总体最优部最优，却不能保证数据的总体最优，因为克里金，因为克里金估值的方差比原始数据的方差要小。因此，当井点估值的方差比原始数据的方差要小。因此，当井点较少且分布不均时可能会出现较大的估计误差，特较少且分布不均时可能会出现较大的估计误差，特别是在井点之外的无井区误差可能更大。别是在井点之外的无井区误差可能更大。克里金方法的局限性克里金方法的局限性106优选内容 （2）克里金插值法为光滑内插方法克里金插值法为光滑内插方法，为减，为减小估计方差而对真实观测数据的离散性进行了小估计方差而对真实观测数据的离散性进行了平滑处理，虽然可以得到由于光滑而更美观的平滑处理，虽然可以得到由于光滑而更美观的等值线图或三维图，但等值线图或三维图，但一些有意义的异常带也一些有意义的异常带也可能被光滑作用而可能被光滑作用而“光滑光滑”掉了掉了。所以，有时，。所以，有时，克里金方法被称为一种克里金方法被称为一种“移动光滑窗口移动光滑窗口”。（用于CCDF的求取，应用于随机建模）(3)107优选内容