参数地质统计学方法综述

参数地质统计学方法综述

0不同属性的特值估计参数地质统计根据不同的研究对象提出了不同的研究方法。如在平稳条件下,研究已知或未知对象的数学期望时,可以使用适合于研究单个变量的简单克立格法和普通克立格法;在研究对象具漂移或趋势性变化时则使用具次级漂移的克立格法或具趋势的克立格法;若各区域化变量之间具有一定程度的相关性时,则应用各种协同克立格法。参数地质统计学以数据服从正态分布为基础,要求样品均匀分布且样本量充足。对于那些偏离正态分布的数据,常需要对数据进行各种处理来满足不同方法的要求。但是地质、物探化探数据的分布常出现一个长尾巴,这种分布特征影响了数据的整体分布规律及数据处理与估计的精度。引起这种非正态分布的因素主要有:①数据中有特异值,这种意义上的特异值有两种含义,可以指特低值,也可以指特高值。如在矿床储量计算时一般的特异值多为品位的特高值,表明矿化的富集作用,它对储量计算的影响很大。②数据由多个总体组成,对成矿来讲很可能是指成矿的多期多阶段或多要素控矿。这种数据一般不能满足正态假设的要求,用普通的克立格技术(即参数地质统计学方法)来研究那显然有些不稳健。为了解决这些问题,参数地质统计学家进行了许多卓有成效的研究,例如,若存在特异值时,根据数据的某种分布对特异值进行一定的处理(剔除或替代);利用平滑函数对数据进行平滑;限制特高品位值的外推;对数据分布具正偏或负偏的数据进行对数转换或某些线性转换等。若数据由多个总体混合而成,则对混合总体进行筛分。所有这些方法都是根据数据的分布特点进行变换,使数据的分布特征满足方法的要求。然而在实际应用中,被研究的母体常存在取样有限且分布不均的现象,特异值的存在使得变异函数不稳定,数据的分布假设无法检验等一系列的问题使得参数地质统计学的应用受到一定的限制,地质统计学者开始寻求一种更稳健的方法来满足研究的需要。AGJournel发展了无须对数据分布作任何假设的非参数地质统计学,提出了一些非参数地质统计学方法(如指示克立格法、概率克立格法、分位数克立格法及各种非参数统计量法)。单元指示克立格法(STIK)并没有充分利用包含在原始数据中的信息:在对连续变量的研究中,STIK只考虑原始数据值是否超过了临界值zl;在离散变量的研究中,如对岩石类型的空间变异性研究中,STIK只关心某类岩石是否出现或缺失,一次只保留一种岩石类型的信息。然而在一定的研究邻域里,实际研究对象中各区域化变量大都具有或多或少的联系,这些区域化变量之间表现出一种协同行为,绝对不相关的几个区域化变量是很少见的。如在油藏描述中,剩余油的分布就与孔隙度、渗透率及饱和度等都有关;在金矿床中对某点的金品位进行估计时,不但要考虑邻域内主要研究元素金的品位的影响,而且还要考虑研究邻域内与金矿床密切共生或伴生的其他元素品位(如铜品位)的影响,因为它们的存在或出现与主要研究变量具有某种必然联系。现代地质统计学考虑到了各区域化变量之间的相互关系,在时空域中的协同区域化理论的基础上,充分利用已有的各种信息,使得研究结果更合理、更稳健、更精确。本文是在空间指示克立格法(IK)和协同区域化理论研究的基础上,对时空多元指示克立格法进行研究。1关于法定改造目标的法—指示二阶距时空多元指示克立格法的研究是在对原始数据进行指示转换的基础上进行的,对不同类型的变量需要用不同的转换方法。原始数据Z[(t,x)a]服从随机分布,依此为基础进行转换得到的指示数据I[(t,x)a;zl]也为随机分布。如矿山储量计算中,金属(如Au)品位的原始数据服从随机分布,则在已知边界品位zl的条件下,金属品位经指示转换后的指示值I[(t,x)a;zl]也服从随机分布,且随机函数I[(t,x)a;zl]服从二项分布。二分指示随机函数I[(t,x)a;zl]的平稳期望是随机函数Z[(t,x)a]的累加分布函数(cdf—cumulativedistributionfunction)本身,其期望值为:E{Ι[(t,x)a;zl]}=1×Ρrob{Ζ[(t,x)a]≤zl}+0×Ρrob{Ζ[(t,x)a]＞zl}=Ρrob{Ζ[(t,x)a]≤zl}=F(zl)。(1)E{I[(t,x)a;zl]}=1×Prob{Z[(t,x)a]≤zl}+0×Prob{Z[(t,x)a]＞zl}=Prob{Z[(t,x)a]≤zl}=F(zl)。(1)式中的F(zl)表示品位分布函数在边界品位zl处的概率值。当I[(t,x)a;zl]和I[(t+τ,x+h)a;zl]为被矢量(τ,h)分隔的随机变量时,可定义如下二阶距:①非中心化协方差:ΚΙ(τ,h;zl)=Ρrob{Ζ[(t+τ,x+h)a]≤zl}和kΙ(τ,h;zl)=Ρrob{Ζ[(t,x)a]≤zl}‚(2)即非中心化协方差是2个点的品位值均小于或等于边界品位值的概率。②非中心化互协方差:ΚΙ(τ,h;zl,zl′)=Ρrob{Ζ[(t+τ,x+h)a]≤zl}‚ΚΙ(τ,h;zl,zl′)=Ρrob{Ζ[(t,x)a]≤zl′}‚(3)即非中心化互协方差是2个点的品位值均小于各自边界品位的概率。③中心化协方差:CI(τ,h;zl)=KI(τ,h;zl)-F2(zl)。(4)④方差:Var[Ι(t,x;zl)]=CΙ(0,0;zl)=F(zl)-F2(zl)。(5)⑤半变异函数:γΙ(τ,h;zl)=0.5×E{[Ι(t+τ,x+h;zl)-Ι(t,x;zl)]2}=CΙ(0,0;zl)-CΙ(τ,h;zl)=F(zl)-ΚΙ(τ,h;zl)。(6)⑥相关系数:γΙ(τ,h;zl)=CΙ(τ,h;zl)/CΙ(0,0;zl)=1-[γΙ(τ,h;zl)/CΙ(0,0;zl)]。(7)⑦指示半变异函数与指示协方差的关系:CΙ(τ,h;zl)=CΙ(0,0;zl)-γΙ(τ,h;zl)。(8)2样质结构参数的确定STCOIK充分利用研究邻域内所有具相关关系的区域化变量所能提供的信息,将IK技术与协同克立格技术在时空域中有机地结合起来进行研究。设有K个区域化变量,在某待估点(t,x)at有K个指示值{Ik[(t,x)atk;zk],k=1,2,…,K},对主变量k0的指示值的STCOIK估计表达式为:Ιk[(t,x)atk0;zk0]*SΤCΟΙΚ=Κ∑k=1Τ∑t=1Νtk∑atk=1λatk.Ιk[(t,x)atk;zk]‚(9)相应的协同克立格方程组需要K2个直接、互指示协方差矩阵:CΙ[(τ,h)atk;zk,zk′]=Cov{Ιk[(t,x)atk,zk],Ιk[(t+τ,x+h)atk′,zk′]}。(10)根据主变量与次级变量的不同作用,可以选用不同的协同克立格技术来进行各种估计:如主要考虑主变量的作用时用普通协同克立格法;当主变量与次级变量同样重要时,则可选用标准协同克立格法等。对于K较大时,直接推断和联合模拟K2个(互)协方差函数就有点不现实。一个解决方案是要求有一个先验二元分布模型;另一个方案是对那些不具互相关的指示变量值进行线性转换。采用先验二元分布模型相当于在参考一些或全部指示协方差模型的基础上模拟实际数据。用得最多的二元分布模型是对原始数据进行正态转换后的二元高斯模型:Ζ(W)→Y(W)=φ(Ζ(W)),(11)通过一些关系式,所有的指示(互)协方差都可由经正态转换后的Y协方差来决定。在析取克立格法(STDK)中的等因子模型对二元高斯模型进行了改进,这种改进体现在对原始数据Z(W)的Ψ(·)转换,这种非线性转换不同于正态转换φ(·)。这一改进除了引进Ψ(·)函数外,其他就没有什么与二元高斯模型有本质差别的地方。3stpck估计的基本原理对时空域中相距(τ0,h0)的指示协方差矩阵{CI[(τ0,h0)atk;zk,zk′],k,k′=1,2,…,K}利用克立格法求其主分量而形成STPCIK,通常将时空域中的空间距离h0和时间距离τ0取的很小或为0。K个主分量{Yk[(t,x)atk],k=1,2,…,K}是原来的K个指示变量的线性组合:Yk[(t,x)atk]=Κ∑k′=1ak,k′⋅Ιk[(t,x)atk′;zk′]‚(12)即Y=AT·I。其中K×K阶矩阵AT=[ak,k′]是正交矩阵A的转置矩阵,正交矩阵A是从协方差矩阵进行谱系分解而得到:CI[(τ0,h0)atk;zk,zk′]=A∧AT,(13)∧是具有K个指示协方差矩阵特征值的对角矩阵。通过这种构造,以时空矢量(τ0,h0)分隔的K个主分量Yk(t,x)相互正交:Cov{Yk[(t,x)atk],Yk′[(t+τ,x+h)atk′]}=0,∀τ,h,∀k≠k′。(14)通过协同克立格方程组来求得K个主分量Yk(t,x)的值,这大大减小了用单元克立格来进行估计和推断的花费。此外,当区域化变量呈连续分布时,将K个主分量按时空相关性以降序排列,高阶主分量与纯块金效应行为有关:Cov{Yk[(t,x)atk],Yk[(t+τ,x+h)atk]}≈0,h≠0,k接近Κ。(15)因而,对高阶主分量Yk(t,x)进行克立格估计实际就是在滑动研究邻域内对低阶主分量Yk[(t,x)atk]的简单平均。那么,原来协同克立格方程组中的K个指示值Ik[(t,x)atk;zk]就可以减少到按降序排列的前K′个,K′≤K。从指示克立格求出的主分量经线性转换就可以得到相应的对指示值的大致无偏的协同克立格方程组估计:Ik[(t,x)atk;zk]*SΤΡCΙΚ=Κ∑k′=1ak′,k⋅Yk′[(t,x)atk′],(16)即I=A·Y。现有的经验表明,对于连续的区域化变量,STPCIK的计算结果并没有明显好于直接的指示克立格计算,STPCIK主要用在评价不同分类变量的条件累加分布函数(ccdf—conditionalcumulativedistributionfunction)。在指示克立格中,对于多个变量信息的提取自然会用到协同指示克立格。STPCIK就是在协同指示克立格的基础上,结合主因子分析方法,对多个区域化变量的重要性进行评估,对区域化变量进行优化组合。4基于不同数据集的二次对接指示克立格法的灵活性以出现序次相关问题为代价,在有些特定的指示克立格研究中,有时可能会碰到一个或者一半甚至三分之二的ccdf函数值发生序次偏离的现象。通常,序次相关现象较少出现,平均的偏离概率一般只有百分之一左右。有2个原因会产生序次相关现象:1)负的指示克立格权值,解决办法是限定指示克立格方程的解非负;2)对于连续性变量,在某些分类中缺乏数据。实践表明,大多数序次相关问题都是与缺少数据有关:临界值zl-1和zl之间没有数据。在这种情况下,临界值Zl-1和zl的指示数据是一样的,而相应的指示变异函数模型很可能不同,因而二者的ccdf值就不同,序次相关问题就出现了。针对不同的数据集,有许多可以采取的措施来减少序次相关问题的出现:1)通过用一个在内连续变化的先验参数化核心函数F[(t,x)atk;z]来代替每个二元指示数据集{I[(t,x)atk;zl],l=1,2,…,K},可将原来离散分布的区域化变量z的值平滑成较连续的分布。由于核心函数F[(t,x)atk;z]围绕原始数据Z[(t,x)at]具有一定的改变速度,因而为了得到更连续的ccdf函数,便牺牲了估计精度。2)用逐渐(缓慢)增加或减小临界值的办法来定义或描述指示变异函数参数变化的连续性。用同一个具基本套合结构的变异函数来拟合所有的指示变异函数倒是一个好办法:γΙ((τ,h)at;zl)=C0(zl)+CΙ(zl)Sph(|h|/a(zl))。(17)其中Sph(·)是具变程a(zl)、单元基台的球状变异函数模型。绘出每个参数(如C0,Cl,a)与临界值zl的相关图,以保证参数值的平滑变化。如果平稳假设是合理的,则可将参数值平滑改变。此外,根据相关图,也可以对超过临界值zl的相应参数值进行内插和外推。选择稳健的中值指示变异函数将变成中值指示克立格法(STMIK),中值指示克立格法以欠灵活为代价戏剧性地减少了序次相关偏离的数目。3)将临界值取与研究