地质数据处理_14－线性平稳地质统计学.ppt

《地质数据处理_14－线性平稳地质统计学.ppt》由会员分享，可在线阅读，更多相关《地质数据处理_14－线性平稳地质统计学.ppt（71页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、1,线性平稳地质统计学,洪金益中南大学地学院,地质数据处理基础14,2,第14章 线性平稳地质统计学,一、随机场与区域化变量1定义：以空间点x的三个直角坐标xu,xv,xw为自变量的随机场 Z(xu,xv,xw）=Z(x）称为一个区域化变量。 区域化变量具有两重性： 观测前，将Z(x）看作随机场；观测后，将Z(x）看作一个普通的三元实值函数。即空间点函数，一次观测后，就得到它的一个实观Z(x)。,第一节 区域化变量的理论,3,2.功能 能同时反映地质变量的结构性与随机性。 当空间点x固定后， Z(x）即为一个随机变量； x与x+h两点处的Z(x）具有某种程度的相关性（因随机场有相关函数R（x，

2、x+h）即为一个随机变量； 3.物理学或地质学特征 空间局限性；不同程度的连续性；不同类型的各向异性。,4,1. 协方差函数 若Z(x）是随机场，在空间两点x和x+h 处两个随机变量Z(x）和 Z（x+h）的二阶中心混合矩 称为随机场的Z(x）自协方差函数，简称协方差函数。一般地讲，它是依赖于点x和向量h 的函数。 特殊地：当h =0时， 就等于方差函数: 当其不依赖于x时简称方差，故有:,二、协方差函数与变差(异)函数,5,基本公式,6,在二维、三维情况下定义时，以一维变差函数为基础，需考虑各向异性，结构套合等问题。 当r（x，h）与x的取值无关时，r（x，h）只依赖与h（滞后、间隔、步长）

3、，则可将r（x，h）写成r（h），此时以h为横坐标，r（h）为纵坐标作出图形谓之变差(异)图。,7,问题：由数理统计知：要估计变差函数值 就要估计数学期望值 这必须有若干对Z（ x ）和Z（ x+h ）的值才可通过求 平均数的办法来估计上述数学期望。而这在实际地质，采矿工作中是不可实现的，因为不可能恰在空间同一点上重复直接取得二个样品。这就使统计陷入困境。需借助假设来解决。,三平稳假设与本征假设,两个重要的假设条件:,1. 平稳假设2. 本征假设,8,1 平稳假设 严格的平稳假设 区域化变量Z（x）的任意n维分布函数不因空间点 x发生位移h而改变。 即： 这种要求是Z（x）的各阶矩存在，且平稳

4、，这在实际中不能满足，且不好验证。所以实用上采用的只需一、二阶矩且平稳就够了。 二阶平稳（弱平稳）。,9, 二阶平稳假设满足下列两个条件1）整个研究区内，Z（x）的数学期望存在，且等于常数，2）整个研究区内，Z（x）的协方差函数存在且平稳（即只依赖于滞后h，而与x无关） 特殊地：当h=0时 =C（0）即方差存在且为常数。当上述条件仍不能满足时，条件进一步放宽，导致本征假设。,10,2. 本征假设 区域化变量Z（x）的增量Z（x）- Z（x+h）满足下列两个条件： 1） 在整个研究区内有: 2）增量Z（x） Z（x+h）的方差函数存在且平稳（不依赖于x）即： =2r(h), 即Z（x）的变差函数

5、存在且平稳。,11,3 .二阶平稳假设与本征假设的比较 总的结论：二阶平稳假设较强，本征假设较弱1） 由二阶平稳假设的第一个条件可推出本征假设条件一。 如：设 y为一服从柯西分布的随机变量，其概率密度为 则： ,不存在但： ，存在且为0,12,1） 二阶平稳假设的第二个条件可以推出本征假设条件之二 在二阶平稳假设满足时： 由二阶平稳假设条件之二 =C（0）， ，当h=o 故：同理有：而由h0 时的二阶平稳假设条件二有：则： 只要协方差函数存在，则C（0）存在，于是r(h)存在 ,13,4 准二阶平稳假设与准本征假设,区域化变量在整个区域内并不满足二阶平稳（或本征）假设而在有限的领域（如以X为中

6、心，X为半径的圆）内是二阶平稳（本征）的，则称区域化变量Z（X）是准二阶平稳（或准本征）的。 这才是在大多数情况下适用的，有了这一假设，我们便可根据N对z(x)和z(x+h)（i=1，2，n)的数值，通过求某种平均值的办法来估计变函数值了。,14,把x轴上相隔为h的N（h）对点xi 和xi+1 （i=1，2，N（h）处的N（h）对观测值Z(xi )和Z(xi+1 )（i=1，2，N（h）看成是Z（x）和Z（x+h）的N（h）对实现。于是一维实验变差函数 为： 例1：设Z（X）为一维区域化变量，满足本征上假设，又已知： Z（1）=2，Z（2）=4，Z（3）=3，Z（4）=1， Z（5）=5，Z（

7、6）=3，Z（7）=6 ，Z（8）=4,四. 实验变差函数的计算,15,五估计方差,定义：估计方差就是估计误差的方差，记作1 估计误差与估计方差 设平面上 处打一垂直钻孔，定义以下单元： 1） 块段V，其实平均品位 ，体积V 2） 若芯平均品位 ，钻于块段V的中心 若以 来估计 时，其估计误差为： 简记为一区域化变量在点 的实现，设R（x）是二阶平稳的（1）（2）,16,即有估计方差存在， ：表示平均估计误差的大小 ：表示估计误差对其分布中心 的离散程度按照统计上的参数估计，好的估计量应是无偏差和有效的，即方差最小，若 具有这两种性质，即有： 或 又若R（x）N（0， ），则由 和 构成 的置

8、信区间：95%的置信区间概率：（ -1.96 , +1.96 ）90%的置信区间概率：（ -1.64 , +1.64 ）,17,2 线性估计量的函数形式 n个样品品位值为,估计中心在x,体积为v的块段的平均品位ZV(x)。 一般情况下，估计量是诸的函数，设为： 我们总希望找到的这一函数f，使Z* 满足： 1）无偏: =0 2） 估计方差： = 最小 要计算上述两个期望值，需已知多维变量的联合分布，而我们仅掌握了其一个实现，因而无法求解问题。 采用一种限制，设该函数f为线性函数，采用线性估计量。求取方程的系数，即可得到该函数形式。 （ 为权系数）,18,3 估计方差的计算公式（算术平均构成） 设

9、点品位Z（x）是一个二阶平稳的区域化变量，具有期望值m，协方差（函数）C（h），变差函数 待估块段：中心点在x，体积为V。 信息样品：中心点在x，体积为v n个中心点在x的点集 被估计的平均品位： 只有一个信息样品时，中心在x ，体积v 无偏性： 又：,19,20,21,4 关于估计方差的几点说明： (1) 用域v 内的信息对块段V作估计求的 ，也称为v对V品位的外延方差，记 （v,V）。 (2) 不论v 和V是怎样的域， 的两个公式都完全适应V可以是不同的块段V1,V2， V= V1+V2 (3)可将(h)视为用Z（x+h）估计Z（x）的估计方差之半 (4) 的大小可以衡量估计量的优劣， 越

10、小，估计量越好， 的大小与许多因素有关： 待估块段V与信息样品v间的距离，距离越远，越差。 待估块段V的几何特征（大小、形状），块段越大，越好 信息样品v的几何特征，数量和空间排布，v大，样品多，相距远，越好。 变差函数的特征（矿化结构和空间连续性）。,22,5. 用加权平均作估计量的估计方差公式用 构成估计量来估计中心在x处块段v的平均品位值要满足无偏性时，就应有条件 。证明： 若要 0，即要而故：要 =m，则要 =1即：各权系数之和等于1,23,按照 的变差函数形式：式中： 表示当向量 的一端固定在点 ，另一端独立地扫描过体积V 时变差函数r（h）的平均值。要计算估计方差 ，首先要计算出变

11、差函数r(h)，然后再求r(h)的某种平均值。,24,六、承载效应和离散方差,1承载效应 承载对离散程度的影响，谓之承载效应。离散程度取决于两个因素： （1）变量（品位）的变化域V（离散域、开采面） （2）生产单元承载v（统计单元） P.Delfiner 的例子： 薄片面积V上按1818=324个规则网格测量，统计单元分别按：11、 22、 33、 66 时，则： 211=22.31、 222=6.42、 233=5.11、 266=1.99,25,2离散方差 设V表示以点x为中心的平面矩形开采面，将V分成N个形状相同，大小相等，分别以x i (i=1,2,N)为中心的生产单元（简记 v (x

12、 i )，均等于v，也为平面矩形） 于是： 设Z（y）为点y处的品位，则每个以点 x i 为中心的生产单元v (x i )的平均品位为: 以x为中心的开采面V的平均品位为：ZV(x) 是诸 Zv(xi) 的算术平均值,26,N个品位值Zv(xi) (i=1,2,N) 对它的平均值ZV(x) 的离散程度可用其方差来表示，即： 当x固点时（只要N确定，则xi 也确定）， Zv(xi) 和 ZV(x) 均为随机变量，故S2(x) 也是随机变量，可以讨论其数学期望。 在区域化变量z（y）（点品位）满足二阶平稳条件下， S2(x) 定义为在开采面V内N个生产单元v 的离散方差，记作：,27,3离散方差的

13、计算公式 根据前面估计方差的计算公式： 由z（y）满足二阶平稳假设知道：协方差c（h）是平稳的，不依赖于x或y，因此： 又按离散方差的公式： 导出：影响离散方差的因素： V的大小和形状； v的大小和形状； 变差函数r（h）,28,4 克里格关系式 若 ，则： G：矿床 V： 开采面 v：生产单元 证明:,29,线性平稳地质统计学的三大基本公式：,30,一、协方差数C（h）的性质（在二阶平稳假设下）,第二节 变差函数及结构分析,31,2变量函数 （h）的性质（Z（x）满足二阶平稳假设） （1） （0）=0 （2） （h）0 （3） （-h）= （h） （4） - （h）必须是条件非负定函数（即由

14、- （xi-xj）构成的 矩阵必须是条件非负定矩阵）。具体地，若 成立， 则 - （xi-xj） 为非负定阵。 （5） （）=C（0）,32,33,4. 交叉协方差函数和交叉变差函数的性质 （1）协同区域化 用一组K个相关的区域化变量Z1(x), Z2(x), Zk(x) 来表示的区域化谓之协同区域化 （2）在二阶平稳假设条件下，定义： EZk(x) =mk=常数， x，k=1,2, k 对每对区域化变量Zk(x) 和Zk(x) ，交叉协方差函数为： EZk(x+h). Zk(x)-mk mk=Ckk(h) x 对每对区域化变量Zk(x) 和Zk(x) ，交叉变差函数为： 1/2EZk(x+h

15、)- Zk(x)Zk(x+h)- Zk(x)= kk(h) x,34,（3）交叉协方差函数的性质 当k=k 时，交叉协方差函数（变差函数）变为协方差（变差）函数：Ckk(h)= Ck(h)， kk(h)= kk(h) x kk(h) 可以取负值，而k(h) 总是0， 负相关 交叉变差函数关于k和k 对称，关于h和（-h）对称 kk(h)= k k (h)， kk (-h)= kk (h) 交叉协方差函数： Ckk (h)= Ckk (h) 对k对称 Ckk (-h) Ckk (h) 对h不对称 在二阶平稳假设下，交叉协方差函数和交叉变差函数皆存在，且： 协同区域化中互相矢函数可定义为：在同一点

16、x处两个变量Zk(x)和Zk(x)之间点对点的互相关函数：,35,二、变差函数的功能1通过“变程”反映变量的影响范围2变差函数在原点处的性状反映了变量的空间连续性（1）抛物线型（或连续性） 高度连续性 当|h| 0时， (h) A|h|2 （A为常数）（2）线性型 平均连续性（均方意义下连续） 当|h| 0时， (h) A|h| （A为常数）（3）间断型 （或有“块金效应型”） 连续性很差（无平均 连续性）， (0) =0（4）随机型（“或纯块金效应型）（5）“过度型” 介于（1）和（4）之间3不同方面上的变差图反映矿化的各向异性。,36,37,设Z（x）是满足本征假设的区域化变量，它具有各向

17、同性的变差函数(h) ，则常见的变差函数 理论模型有：,三变异函数的理论模型,38,39,三种有基台值模型的比较,40,（4）幂函数模型： 实践上， 常采用线性模型： 注意 必须严格地小于2，因2，则（-r ）不再是条件非负定， r 就不能作为变差函数。（5）对数函数模型： 六十年代的 DeWijs 模型： 由于当 r 0 时，log r ，这与变差函数的性质不符合。因此，对数函数模型不能用来描述点承载的区域化变量。但却可以用来作为正则化变量的变差函数v(r)的模型。 如：对钻孔岩心样品以l=1 进行正则化后，点对数函数模型 (r)=ln(r) 变为正则化对数函数模型：,41,（6）纯块金效应

18、模型（7）空穴效应模型 当 (r) 并非单调递增，而显示出有一定的周期性的波动时，叫做 空穴效应（也叫孔穴效应） 常见的空穴效应模型公式： 其中：C0块金常数， C拱高， b高品位带的平均距离，a变程（指数模型）,0,1,1,0, (r), (r),r,r,12,=1,1,幂函数模型,对数函数模型,42,四、一个方向的套合结构,由于实际区域化变量并非是上述七种模型中的一种，而多数是多种结构的复合，即往往包含各种尺度上的多层次变化性。应由多种结构的变差函数来叠加，这谓之套合结构。大体上有以下几层结构和原因： （1）岩心采样率的波动，取样误差，以及在样品制备、分析和测定等过程中产生的变化性，反映在

19、变差函数上就是点承载（r0）一级结构； （2）由一种矿物成分转变为另一种矿物成分所引起的变化性，这在金矿、铀矿等品位变化剧烈的矿床上尤为明显（1cm 级）； （3）由岩矿层交替、或矿化透镜体和非矿围岩的交替引起的变化性，反映在变差函数上r 100 m一级的结构； （4）由区域构造运动，岩浆活动所造成的变化性，反映在变差函数上是r 100 km一级的结构 套合结构的表示：以反映各种不同尺度变化性的多个变差函数之和表示：(r)= 0(r)+ 1(r)+ 2(r) + + (r) + 其中，每个成分i(r) 可以是不同模型的变差函数。,43,五、不同方向上的结构套合,1各向异性的概念与种类 若Z(x

20、)的三维变差函数 ，则称变量Z(x)= Z(xu，xv，xw) 为各向同性的区域化变量，反之则为各向异性的。 （1）几何各向异性 当两个方向的变差函数具有相同的基台值C(设块金常数C0为0 )和不同的变程a1， a2 时，称这种各向异性为几何各向异性。（可经线性变换变为各向同性。） （2）带状各向异性： 凡不能通过坐标的线性变换化为各向同性的各向异性。即不同方向的变差函数(h) 都具有不同的基台值，而变程可以不同，也可以相同。,44,几何各向异性,不同变程、不同基台值,相同变程、不同基台值,带状各向异性,二维几何各向异性的方向变程图,45,4结构模型的一般表达式结构模型 (h) 总可看成由N个

21、向向同性结构i(|hi|) 套合而成，即：而i(|hi|) 则是经特定线性变换矩阵Ai 的坐标线性变换由某种各向异性(几何或带状的)结构转化而来的，这种线性变换将原坐标向量h变为新坐标向量hi,46,三、变差函数的拟合,.直观拟合球状模型（） 求变程a （见图）计算所用数据（用于计算实验变差函数*(h)的实验方差*2；）在实验变差函数图的纵坐标轴上过*点作一条平行于横坐标轴的直线；）以直线连接实验变差函数*(h)的头两至三个点，此直线与过*2点的直线相交，交点横坐标为 ，即假定其值为h0，则a =,（）求块金值C0：连接*(h) 的头两三个点的直线与纵坐标轴的交点的纵坐标为C0（见图）。若C0

22、a的情况都很简单，所以仅讨论0ha 的情况：令：y= (h)，x1=h，x2=h3，b0=C0， 则上式变为：y= b0+b2x+b2x2这样对球模型变差函数的拟合问题就变成了多元线性回归问题。,48,3. 用加权多项式回归法拟合二级套合球状模型,公式： (1) 将实验变差函数分两个区分另进行拟合先取前区0ha 来拟合：按前述方法，令：y= (h)，x1=h，x2=h3，b0=C0， (2)用后半区数据点aha 拟合第二段设： y= (h)，x1=h，x2=h3，b0=C0C，,49,七、结构分析的一般步骤,(一)选择区域化变量1选择：就研究目的而选取2样品承载：应具有一致的承载并加以注明；3

23、方法有具可加性：可加性：线性组合后仍保持相同的意义；适当性：能概括研究问题的主要特征均一性：注意观测尺度和空间域。,50,(二)审议数据(1)勘探方案是否合理？取样目的、设备、网格、方式等是否合适？取样方案在实施中有无遇到问题？(2) 数据是否有代表性，采样是否均匀、充分，数据在时间和空间上是否一致？岩心采样率多少，变化的大小，数据是否有系统误差，原因何在及处理方法，数据编录中有无问题等。,51,(三)定长岩心样品的改组(数据的正则化) 1. 样品不定长，但全部进行了分析，即将非等间格样品以长度加权平均分为等间格样品 2. 样品不定长，又未全部分析 这与选择性取样分析有关 : 模拟缺失分析地段

24、，根据地质情况赋以特定值 将变差函数计算限制在所分析的相带内 (四) 数据统计 1平均值、方差、相关系数矩阵 2直方图、相关散点图 3. 确定概率分布，确定奇异点，地质特征,52,53,54,3非列线又不等间距的数据构形 先将数据点合成角度组，将角度组内数据视为方向数据； 再将数据点按距离组kh(h)组合，凡落在角度范围(d)及距离范围kh(h)的数据均认为是与x0点在方向上相距为kh的数据点。,55,(六)、理论变差图的最优拟合。(七)、结构套合和比例效应的消除。(八)、变差函数参数的最优性检验 (1)“交叉验证”法(Jacknife法)：各测点上克里格估值与实测值误差平方最小(2)离散方差

25、检验法准则：克里格估值与测值标准差S*，好的变差图综合指标： I值越小，变差函数图参数确定越好。(九)对变差函数进行地质解释,56,第三节 整体估计,一、对 Z(V) 的整体估计 设 ，各钻孔位置xi 位于vi 的中心，有： （1）若各vi 均相等，则：Z(V) 的一个无偏估计量为： 因为： 而： 所以： （2）若各vi 不都相等时，则可取Z(V) 的估计量为：该估计量也是无偏的。 因为：,57,二、整体估计方差 v 表示全体信息样品承载的集合，于是： （1）由基本估计方差组成的法则 规则网格和不规则网格取样的情况 基本估计方差 用Z(xi) 估计V(vi) 时的估计方差。 当 Z(xi) 和

26、Z(xj) （ij） 相互独立时，可用 组成整体估计方差,58,规则网格取样的情况,（vi 均相等），设为v , 则： 由前已知， 所以： 由于当ij 时，Z(xi) 和 Z(xj) 独立，第二式为零。 所以：,59,不规则网格取样的情况,vi 都不相等，仍有 时： 而： 于是： 其中： 是以Z(xi) 估计Z(vi) 的估计方差。由于诸vi 都不相等，所以各估计方差也不一定相等。 特别地：当vi =v， i 时，有V=Nv，于是上式变为：,60,第四节 局部估计普通克里格法,1. 克里格法的定义 2. 克里格法的种类 3. 克里格法的使用信息和应用条件 4. 克里格方程组 5. 克里格方差,

27、61,1. 克里格法的定义 根据一个块段（盘区）内外的若干信息样品的某特征值（品位）数据，对该块段（盘区）品位（特征值）作出一种线性、无偏、最小估计方差的估计方法。2. 克里格法的种类 （1）普通克里格满足二阶平稳（或本征）假设的区域化变量 （2）泛克里格非平稳或具有漂移存在的区域化变量 （3）析取克里格计算可采储量时，需要采用非线性估计量时 （4）对数克里格区域化变量符合对数正态分布时 （5）随机克里格数据较少、分布不规则、对估计精度要求不高时 （6）指示克里格处理非参数（类型参数）的区域化变量时,62,63,用盘区（平均）品位y对样品（平均）品位x的回归直线，则不是y=x，而是一条斜率为（

28、m 时： 但由于m未知，以样品的平均值 代替 若取x=x1，则盘区估计值为： 设： 则： 这就是克里格法的雏形。各种克里格方法就是不同的权值估计法。,64,65,二、克里格方程组及其方差,66,2. 线性估计量的构造 Zi （i=1，2， ，n）是一组离散的信息样品数据，它们定义在点承载xi （i=1，2， ，n）上的；或是确定在以xi 点为中心的承载vi （i=1，2， ，n）上的平均值Zvi (xi) （简记Zi ）。且这n个承载vi （i=1，2， ，n）既不同于V，又各不相同。所采用的线性估计量为： 它是n个数值的线性组合。 克里格估值的原则：就是在保证这个估值ZV*是无偏的，且估计方

29、差最小的前提下，求出n个权系数i 。在这样的条件下求得的i 所构成的估计量ZV*称为ZV的克里格估计量，记为ZK* 。这时的估计方差称为克里格方差，记为K*。 当Z (x)的期望已知时：为简单克里格；未知时：为普通克里格,67,3. 简单克里格方程组和简单克里格方差（E(Z(x)=m 已知） 由于Z (x)的期望为已知，令：Y(x)=Z(x)-m 则：EY(x)=EZ(x)-m=EZ(x)-m=0 其协方差函数为： EY(x) Y( y) =C(x, y ) 对ZV的估计转化为对YV* 的估计，且有： 所用的估计量为： 只要求得YV的估计值YV* ，就能得到ZV的估计值ZV * 。,68,显然： YV*是YV的无偏估计量，且不需要任何条件： 为了求出i，使得 最小，需求 的公式： 所以： （1） 其中 表示协方差函数在待估域V上的平均值。,69,70,公式（1）与公式（4）中，所用的估计方差符号不一样，（1）式表示无偏估计量的估计方差，不能保证估计方差最小，故用记号 。公式（4）是在确保估计方差最小的前提下推导出来的，它是克里格方差，故记号为 。其中关键的区别在于i （i=1，2， ，n）在两个式中的意义不一样。 从克里格方程组解出i 后，即得到YV的简单克里格估计量： 所以：,71,结 束,