GS_线性克立格估计

GS-线性克立格估计信息管理学院王玉兰E-mail:[e-mail],[e-mail]:84073385(o);13551261068克立格估计方法概述•克立格法是法国Mathron以南非矿山工程师D.G.Krige的名字命名的一类方法。简单地说是一种特殊的加权移动平均法。其基本是一种对空间分布数据求最优线性无偏内插估计量的方法。•关于估计方法：数学期望•线性估计•非线性估计克立格估计方法概述•常用的估计方法：取已知样本的加权平均值为估计值•假设：已知有n个中心点在xi、大小为vi、对应的区域化变量值为zi的已知样本点，•要求：估计以x为中心点大小为V的块段的区域化变量值zV*•求解方法：已知信息的加权和个样本的加权因子第izziniiiV:1*x1xnxjx2xix3x克立格估计方法概述•估计的一般要求：•1)无偏估计误差的期望为0；•2)方差最小误差平方的期望最小•Kriging（克立格估计法）：一种最佳线性无偏的局部估计方法(BLUE)。•实质：在考虑了信息样本的形状、大小及其与待估计块段之间的位置等几何特征及变量值的空间结构之后，为了达到线性、无偏和最小估计方差的估计，而对每个样本赋予一定的权系数，最后加权求和估计待估块段变量值的方法•关键：求权系数212*2*1*][][:VniiiVVEVViniiiVzzEzzEzzizz个样本的加权因子第克立格估计方法概述•克立格方法分类：•根据估计式形式、区域化变量平稳性、分布及涉及的区域化变量个数等有如下分类•线性克立格法•非线性克立格法(对数克立格、析取克立格等）•参数克立格法•非参数克立格法（指示、概率）•平稳克立格法（普通）•非平稳克立格法（泛克立格）•单变量•多变量个样本的加权因子第izziniiiV:1*x1xnxjx2xix3x普通克立格估计方法•假设所研究的区域化变量Z是满足二阶平稳或内蕴假设，即E[Z]=m,γ(h)或C(h)存在且平稳。•估计量为最佳线性无偏估计•1)无偏性：ijivvjijiVviiVVnininjjijiiiEnininjjijiiiEniiniiiVxdxxCdxvvvvCxdxxCdxVvvVCxdxxCdxVVVCvvVVvVvvCvVCVVCzz)(1),()(1),(,)(1),(),(),(),(2),(),(2),()21)121112111211*或普通克立格估计方法•2)最佳：估计方差最小•即在无偏条件下求权系数使估计方差达到最小•条件极值问题或约束优化问题•方法：拉格朗日乘数法（Lagrange)0/,,1,0/)1(2),(),(2),()21)11211121FniFFvvCvVCVVCiniiEnininjjijiiiEnii普通克立格估计方法•3)普通克立格方程（OrdinaryKrigingSystems)•即权系数满足的方程组1,,1),,(),(1,,1),,(),(0/,,1,0/)1(2111112njjinjjijnjjinjjijiniiEniVvvvorniVvCvvCFniFF普通克立格估计方法•4)普通克立格方差（KrigingVariance)),(),(),(),(:),(),(),(2),(),(2),(1,,1),,(),(1,,1),,(),(1212111211121111VVvVorvVCVVCvvVVvVvvCvVCVVCniVvvvorniVvCvvCniiiKniiiKnininjjijiiiEnininjjijiiiEnjjinjjijnjjinjjij有或由普通克立格估计方法•5)普通克立格法的矩阵表示形式•KΛ=M2•其中K称为普通克立格矩阵•解此方程可得所求权系数Λ2*1*1)1)*(1(2212122111:]111[2),(1]1),(),(),([,][1,,1),,(),(MKKvvCKKKKKVvCVvCVvCMniVvCvvCTnnnjinnTTnTnnjjinjjij则克立格方程为令普通克立格估计方法•普通克立格法的有关的问题KΛ=M2•1）关于方程唯一解的问题•当K严格正定时方程有唯一解，从而要求：•协方差函数为正C(h)0；vi与vi不重合•2）关于K与M2的问题•K只与已知样本点的位置有关，因此对利用同一批已知样本估计区域中的未知块段的变量值，只要计算一次矩阵K；•M2也只与已知样本点位置及待估计块段相对位置有关，而与区域化变量值无关。•因此如果区域化变量结构相同，数据构型相同则线性无偏估计的权系数肯定相同普通克立格估计方法•普通克立格法的有关的问题KΛ=M2•3）估计值的无偏内插性能•当V与vi重合时，zV*=zi•估计方差为0普通克立格估计方法•普通克立格法计算KΛ=M2•点克立格/•块段克立格•P1331234S0P0)(22)(,22)0(22))200(*24003(2020)(33hhCChhh普通克立格估计方法•普通克立格法权系数的特点•1）无偏性•2）估计方差最小•3）对称性：各向同性时几何形状、位置相对于V对称的样本权系数相同•4）减弱丛聚效应(declusteringeffect分解串珠效应)•5）屏蔽效应:块金效应大屏蔽效应小普通克立格估计方法•解决问题的过程•1）输入原始观测数据（区域化变量取样及研究区相关数据资料）•2）数据检查及预处理•3）直方图计算•4）计算变异函数•5）结构分析•6）克立格估计•7）绘图•8）编写出报告普通克立格估计方法•青海别勒滩卤水分布规律•渗透系数16610000166200001663000016640000166500001666000040700004075000408000040850004090000409500041000004105000411000041150004120000普通克立格估计方法•青海别勒滩卤水分布规律•渗透系数020004000600080001000012000140001600018000200002200024000LagDistance0100200300400500600700VariogramDirection:0.0Tolerance:22.5ColumnC020004000600080001000012000140001600018000200002200024000LagDistance0100200300400500600700VariogramDirection:45.0Tolerance:22.5ColumnC16610000166200001663000016640000166500001666000040700004075000408000040850004090000409500041000004105000411000041150004120000普通克立格估计方法•作业要求：•用普通克立格估计法分析新疆某地矿化情况•1）输入原始观测数据，任选一个区域化变量•2）数据检查及预处理•3）直方图计算•4）计算变异函数•5）结构分析•6）克立格估计•7）绘图•8）编写出报告泛克立格估计方法•问题的提出•当区域化变量Z不满足二阶平稳或内蕴假设，即γ(h)或C(h)存在且平稳的假设，而是存在漂移时。普通克立格法不适用，如何进行估计？•即当E[Z(x)]=m(x)如何对区域化变量做估计？•泛克立格法：在存在漂移和非平稳随机函数Z(x)的协方差和变异函数为已知的情况下，一种考虑有漂移的求无偏线性估计估计量的方法。泛克立格估计方法•实际存在的两种现象：•1）在某区域中，某区域化变量有明显的整体变化趋势，但在一个小的局部，又是相对平稳的。如中国大陆的海拔标高的东低西高的趋势，但有成都平原、华北平原等的存在）-----普通克立格法•2）从整体上看，某区域化变量是平稳的，但在一个局部范围内却呈现出明显的非平稳情况。如成都平原也有局部高低变化，即区域化变量Z不满足二阶平稳或内蕴假设，即γ(h)或C(h)存在且平稳的假设，而是存在漂移时。----泛克立格法泛克立格估计方法•有关概念•对非平稳的随机变量Z(x)•Z(x)=m(x)+R(x)•漂移(drift)：非平稳随机函数的数学期望•m(x)=E[Z(x)]•Z(x)的一阶矩，反映Z(x)的主要特征•涨落(fluctuation)：R(x)=Z(x)-m(x)•因为E[R(x)]=E[Z(x)-m(x)]=0，如果R(x)二阶平稳，则γR(h)=½E[R(x)-R(x+h)]2•R(x)有其自身特有结构的随机函数•剩余(residual)：ε(x)＝Z(x)－m*(x)•m*(x):Z(x)的无偏估计量泛克立格估计方法•漂移m(x)：Z(x)的规则而连续的变化•涨落R(x)：关于Z(x)的在m(x)附近的期望为０的有特定结构的变化•漂移m(x)的形式•由于漂移比较复杂，对整个研究区不能用简单的分析表达式模拟，一般用模型表示：•对给定的以点x0为中心的邻域V内，有任一x∈V，其漂移m(x)可用以下函数形式描述：klllxfaxfxm0)()()(已知函数（多项式）K+1个待定系数泛克立格法-一般假设和问题提出•一般假设•(1)假设Z(x)具有一非平稳的数学期望m(x)和非平稳的协方差C(x,y)–m(x)=E[Z(x)]–C(x,y)=E[(Z(x)-m(x))(Z(y)-m(y)]=E[Z(x)Z(y)]-m(x)m(y)•(2)假设Z(x)的增量[Z(x)-Z(y)]具有非平稳的数学期望和非平稳的方差函数–E[Z(x)-Z(y)]=[m(x)-m(y)]–1/2D2[Z(x)-Z(y)]=γ(x,y)•(3)假设Z(x)可分解为两部份：–Z(x)=m(x)+R(x)•(4)假设漂移的形式为：•(5)假设γ(x,y)或C(x,y)已知•(6)假设已知Z(x)的一个现实（有一组在研究区内的样本klllxfaxm0)()(泛克立格法-一般假设和问题提出•问题的提出•在上述假设下求：•(1)在x0点处Z(x)值的线性最优估计，或以x0为中心的邻域V上的平均值•(2)漂移m(x)的线性最优估计；或系数al的最优估计klllxfaxm0)()(泛克立格法-非平稳区域化变量的变异函数•涨落的协方差与变异函数•当Z(x)=m(x)+R(x)成立时，有Z(x)与R(x)具有相同的协方差函数或变异函数•只要能求出R(x)的变异函数即可得Z(x)的变异函数•但由于m(x)未知，不能由原始观测值Zi(x)直接求取)(),(:)(),()]()([))]()(())()([())]()(())()([()]()([),(),()]([)]([)]()([)]()([))]()())(()([(),(221221221221hhxxhxy，xRyxyRxRDymyZxmxZDymxmyZxZDyZxZDyxyxCyRExREyRxREyRxREymyZxmxZEyxCRzRzRz可得令设时满足二阶平稳或内蕴假当泛克立格法-非平稳区域化变量的变异函数•变异函数求取•?用趋势分析结果拟合漂移m(x)，用剩余计算变异函数•不行！！！！•变程范围