• / 107
  • 下载费用:5 下载币  

2 kriging

关 键 词:
petrel 建模
资源描述:
克里金插值第二讲克里金方法( , 是以南非矿业工程师 克里格 )名字命名的一项实用空间估计技术,是 地质统计学 的重要组成部分,也是地质统计学的核心。地质统计学主要是为解决矿床储量计算和误差估计问题而发展起来的由法国巴黎国立高等矿业学院 G.马特隆教授于1962年所创立。H. S. 1947)1951) 克里金法,克立格法) :“根据样品空间位置不同、样品间相关程度的不同,对每个样品品位赋予不同的权,进行滑动加权平均,以估计中心块段平均品位”G. 962) 提出了 “地质统计学”概念(法文 表了专著 《 应用地质统计学论 》 。阐明了一整套区域化变量的理论,为地质统计学奠定了理论基础。区域化变量理论克里金估计随机模拟应用统计学方法研究金矿品位1977年我国开始引入克里金插值方法   井眼地震(普通克里金)(应用 随机函数 理论)不仅考虑待估点位置与已知数据位置的相互关系,而且还考虑变量的空间相关性。为一个实值变量,可根据概率分布取不同的值。每次取值(观测)结果 为随机变量 个实现。P一、随机变量与随机函数第一节 基本原理1. 随机变量连续变量:累积分布函数( ({( 条件累积分布函数( (|)({(|;( 离散变量(类型变量):)}(|)({(|;( Z (u)PP不同的取值方式: 估计( 拟( 续型地质变量构造深度砂体厚度有效厚度孔隙度渗透率含油饱和度离散型地质变量(范畴变量)砂体相流动单元隔夹层类型变量① 设 离散型随机变量 ξ的所有可能取值为… ,其相应的概率为P (ξ= k=1,2,…. 随机变量的特征值:(1)数学期望是随机变量 ξ的整体代表性特征数。则当级数 绝对收敛时,称此级数的和为 ξ的数学期望,记为 E(ξ),或 (ξ) = 11设 连续型随机变量 ξ的可能取值区间为 (-∞, +∞),p(x)为其概率密度函数,若无穷积分绝对收敛,则称它为 ξ的数学期望,记为 E(ξ)。 (E(ξ) =  (数学期望是随机变量的最基本的数字特征,相当于随机变量以其取值概率为权的加权平均数。从矩的角度说,数学期望是 ξ的一阶原点矩。对于一组样本:(1为随机变量 ξ的离散性特征数。若数学期望E[ξ)]2存在,则称它为 ξ的方差,记为 D(ξ),或 ),或 σξ2。σξ= 222 ])(E[-)(])()(  从矩的角度说,方差是 ξ的二阶中心矩。(2)方差其简算公式为D(ξ)=E(ξ2) –[E(ξ)]2D(ξ)= E[ξ)]2方差的平方根为标准差, 记为 σξ研究范围内的一组随机变量。}),({ 研究范围((|)(,,)({(|,,;,,( 1111 随机场:当随机函数依赖于多个自变量时,称为随机场。如具有三个自变量 (空间点的三个直角坐标 )的随机场2. 随机函数条件累积分布函数( 二个随机变量 ξ, η的协方差为二维随机变量 (ξ,η)的二阶混合中心矩 μ11,记为 , η),或 σξ,η。协方差 (,η) = σξ,η = E[ξ)][η)]其简算公式为,η) = E (ξη)) ·E(η)随机函数的特征值二、统计推断与平稳要求P任何统计推断( 学期望等)均要求重复取样。但在储层预测中,一个位置只能有一个样品。同一位置重复取样,得到 现实考虑邻近点,推断待估点空间一点处的观测值可解释为一个随机变量在该点处的一个随机实现。空间各点处随机变量的集合构成一个随机函数。区域化变量:能用其空间分布来表征一个自然现象的变量。(将空间位置作为随机函数的自变量)(可以应用随机函数理论解决插值和模拟问题)考虑邻近点,推断待估点,;,,(),,;,,( 1111 严格平稳);();( 对于单变量而言:可从研究区内所有数据的累积直方图推断而得(将邻近点当成重复取样点)太强的假设,不符合实际P当区域化变量 Z(u)满足下列二个条件时,则称其为二阶平稳或弱平稳:E[Z(u)] = E[Z(u+h)] = m(常数 )  x  绕 在整个研究区内有 Z(u)的数学期望存在,且等于常数,即:二阶平稳② 在整个研究区内, Z(u)的协方差函数存在且平稳(即只依赖于滞后 h,而与 , 即(u),Z(u+h)}= E[Z(u)Z(u+h)](u)]E[Z(u+h)]= E[Z(u)Z(u+h)]-㎡= C(h) 特殊地,当 h=0时,上式变为(u)]=C(0),即方差存在且为常数。协方差不依赖于空间绝对位置,而依赖于相对位置 ,即具有空间的平稳不变性。uu+h① 在整个研究区内有E[Z(u)-Z(u+h)] = 0 本征假设当区域化变量 Z(u)的增量 [Z(u)-Z(u+h)]满足下列二条件时,称其为满足本征假设或内蕴假设。可出现 E[Z(u)]不存在,但 E[Z(u)-Z(u+h)]存在并为零的情况(x)]可以变化 ,但 E[Z(u)-Z(u+h)]=0(比二阶平稳更弱的平稳假设)② 增量 [Z(u)-Z(u+h)]的方差函数 (变差函数 ,存在且平稳 (即不依赖于 u),即:(u)-Z(u+h)] = E[Z(u)-Z(u+h)]2-{E[Z(u)-Z(u+h)]}2= E[Z(u)-Z(u+h)]2 = 2γ(u,h) = 2γ(h), 相当于要求: Z(u)的变差函数存在且平稳。例:物理学上的著名的布朗运动是一种呈现出无限离散性的物理现象,其随机函数的理论模型就是维纳 程 (或随机游走过程 )。可出现协方差函数不存在,但变差函数存在的情况。既不能确定验前方差,也不能确定协方差函数。但是其增量却具有有限的方差:(x)-Z(x+h)] = 2 = A·|h| (其中, ,变差函数 = ·|h|,且随着 |h|线性地增大。2A)(h若区域化变量 Z(x)在整个区域内不满足二阶平稳 (或本征假设 ) ,但在有限大小的邻域内是二阶平稳 (或本征 )的,则称 Z(x)是准二阶平稳的 (或准本征的 )。准二阶平稳假设及准本征假设设 为区域上的一系列观测点,为相应的观测值。区域化变量在 处的值 可采用一个线性组合来估计:三、克里金估计基本思路,1     ,1 0x  0*            m 0*0 估计方差最小 被作为 选取的标准i 可知 为常数 ,有        0*11000可得到关系式 : 11 1)无偏条件Z*(在搜寻邻域内为常数,不同邻域可以有差别)      1,021200*   ( 2)估计方差最小                m i 00*00*用拉格朗日乘数法求条件极值Z*(     1 进一步推导,可得到 n+1阶的线性方程组 , 即克里金方程组当随机函数不满足二阶平稳 ,而满足内蕴(本征)假设时 ,可用变差函数来表示克里金方程组如下 :     1 Z*(小的估计方差 ,即克里金方差可用以下公式求解 :       00102  Z*(差函数 (或叫 变程方差函数 ,或 变异函数 )是地质统计学所特有的基本工具。它既能描述区域化变量的空间结构性变化,又能描述其随机性变化。跃迁现象1. 变差函数的概念与参数四、变差函数及其结构分析),( 设空间点 将区域化变量 Z(x)在 x, x+之差 的方差之半定义为 Z(x)在 为一维情况下的定义:(x)-Z(x+h)] E[Z(x)-Z(x+h)]2-{E[Z(x)-Z(x+h)]}2),( 21==21半变差函数 (或半变异函数 )在 二阶平稳假设,或作本征假设 ,此时:地质统计学中最常用的基本公式之一。E[Z(x)-Z(x+h)] = 0 h(x)-Z(x+h)] E[Z(x)-Z(x+h)]2-{E[Z(x)-Z(x+h)]}2),( 21==21E[Z(x)-Z(x+h)]2),( 1=则:)()0()( (二阶平稳假设条件下)变程 (: 指区域化变量在空间上具有相关性的范围。在变程范围之内,数据具有相关性;而在变程之外,数据之间互不相关,即在变程以外的观测值不对估计结果产生影响。具不同变程的克里金插值图象块金值 (:变差函数如果在原点间断,在地质统计学中称为“块金效应”,表现为在很短的距离内有较大的空间变异性,无论 个随机变量都不相关 。它可以由测量误差引起,也可以来自矿化现象的微观变异性。在数学上,块金值 果品位完全是典型的随机变量,则不论观测尺度大小,所得到的实验变差函数曲线总是接近于纯块金效应模型。当采样网格过大时,将掩盖小尺度的结构,而将采样尺度内的变化均视为块金常数。这种现象即为块金效应的尺度效应。块金效应的尺度效应1 21 11 3 333基台值 ( 代表变量在空间上的总变异性大小。即为变差函数在 块金值 高 高 为在取得有效数据的尺度上,可观测得到的变异性幅度大小。当块金值等于 0时,基台值即为拱高。= C(0) – C(h))(h几何各向异性: 变差函数在空间各个方向上的 变程不同 ,但 基台值不变 (即变化程度相等)。这种情况能用一个简单的几何坐标变换将各向异性结构变换为各向同性结构。带状各向异性: 不同方向的变差函数具有 不同的基台值 ,其中 变程可以不同,也可以相同 。这种情况不能通过坐标的线性变换转化为各向同性,因而结构套合是比较复杂的。地质变量相关性的各向异性1 21 11 3 333(2)2. 变差函数的理论模型设 Z(x)为满足本征假设的区域化变量,则常见的理论变差函数有以下几类:球状模型指数模型高斯模型幂函数模型空洞效应模型接近原点处,变差函数呈线性形状,在变程处达到基台值。原点处变差函数的切线在变程的 2/3处与基台值相交。 2123003球状模型:数模型:   x e x 变差函数渐近地逼近基台值。在实际变程处,变差函数为 模型在原点处为直线。高斯模型:   223e x 变差函数渐近地逼近基台值。在实际变程处,变差函数为 模型在原点处为抛物线。幂函数模型:   幂函数模型为一种无基台值的变差函数模型。这是一种特殊的模型。当 =1时,变差函数为一直线,即为线性模型,这一模型即为著名的 布朗运动(随机行走过程) 的变差函数模型;当 1时,变差函数为抛物线形状,为 分数布朗运动 (变差函数模型。布朗运动分数布朗运动分数布朗运动h 21  11       2c o se x 差函数并非单调增加,而显示出一定周期性的波动。模型可以有基台值,也可以无基台值;可以有块金值,也可以无块金值。空洞效应在地质上多沿垂向上出现,如富矿层与贫矿层互层、砂岩与泥岩频繁薄互层等等。( (h建相应的变差函数模型 , 以表征该变量的主要结构特征。(1)数据准备区域化变量的选取 、数据质量检查及校正 、数据的变换 (如对渗透率进行对数变换)、数据的统计 (如分相对储层参数计算平均值、方差,作直方图、相关散点图等)、丛聚数据的解串 等。3. 区域化变量的 结构分析(2)实验变差函数的计算实验变差函数是指应用观测值计算的变差函数。对于不同的滞后距 h,可算出相应的实验变差函数 。)(* h = N ( h )1h ) ]Z ( Z ( h )21一维实验变差函数的计算公式(i=1,…,N(h))[Z(Z(xi+h)]2的算术平均值一半即为一个 h,进行计算,得出各个 * h = N ( h )1h ) ]Z ( Z ( h )21h 3h 5h (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 )1h ) ]Z ( Z ( h )21721)1(*)2(*)3(*===[22+12+22+42+22+32+22] = 1442 = 12+32+22+22+12+12] = 1220 = 12+12+02+52+12] = 1028 = 1)分不同方向 ,进行 1增加垂向方向( 2)确定主变程方向次变程方向角度容限步长容限h 3h 5h 虑主变程方向的 走向、倾向和倾角)(3)理论变差函数的最优拟合与结构套合选择合适的理论变差函数模型,同时还需进行结构套合,从而得到一条反映不同层次(或不同空间规模)结构的、统一的、最优的变差函数曲线。球状模型指数模型高斯模型幂函数模型空洞效应模型复杂的区域化变量往往包含各种尺度上的多层次、多方向的变化性,反映在变差函数上即为多层次结构。将不同结构组合为统一结构的过程称为“结构套合”结构套合各层次套合例如,对于 200米宽的河道,在 h= 50分出来,但却无法区分 层理和矿物成分的变化性 (即无法找出更细微的结构来 ),它们在 50金效应”出现。若观测尺度为 500米,河道的变化也只能作为“块金效应”。1 21 11 3 333大尺度的变化性总是包含着小尺度的变化性,但却不能从大尺度的变化性中区分出小尺度的变化性。)()()()( 210  )(0 r = 。> 0,,0,00 1 r = 1 1131311 , ),2123(>2 r 2 2232322 , ),表微观变化性的变程极小的球状模型,可近似地看作纯块金效应型球状模型,没有块金常数,基台值为 程为 反映了小规模范围的变化球状模型,没有块金常数,基台值为 程较大,为 反映了大规模范围的变化可以用反映各种不同尺度变化性的多个变差函数之和来表示一个套合结构。(各层次理论模型可以不一样))(可以是不同模型的变差函数其中 21 则套合结构的表达式为)(r2210213232210133221122110a r , ),2123(0 ,r)(21)(230,03><< 。> 0,,0,00  1 1131311 , ),2123(> 2 2232322 , ),0 r)(1 r)(2 r===① 对于 几何各向异性 ,先根据异向比压缩距离轴,使之成为各向同性的模型;②对于 带状各向异性 ,运用模型叠加的方法加以处理。先用压缩距离轴的办法,使其变程变为相同,然后再把具有相同变程的两个球状模型叠加起来,构成一个新的球状模型各方向套合( 将各向异性套合为各向同性 ,以便于在克里金估计时,不同方向均可用统一的结构模型计算实际的变差函数值) (4)变差函数参数的最优性检验:变差函数是否符合实际,应该进行检验。一种实用的检验方法为 “ 交叉验证法 ”( 检验标准是在各实测点,根据周围点计算的 克里金估计值与该实测值的误差平方 平均最小。估计误差的平方 与 克里金估计方差 之比越接近 1,则说明变差函数与实际的符合程度越高。实际上,这种方法在检验变差函数的同时,也在检验所使用的克里金估计方法的适用性。Z*(、克里金插值中权系数的确定     1 (以普通克里金为例)协方差);解克里金方程组在结构分析的基础上设有一个油藏,在平面上 孔隙度值分别为 此估计 0设孔隙度 Z(x)是二阶平稳的。其在平面上的二维变差函数是一个各向同性的球状模型,其参数为:块金值 2,变程 a= 200,拱高 C= 20,即:实例200,h 22,200, ),( 2 0 0 )0,h 0,( h ) 33 41 普通克里金方程组的矩阵形式为 [K][ ]=[0 1 1 1 1 1 C C C C C C C C C C C C C[ K ],][444342413433323124232221141312114321,1  ][][][21      1 (求解)0 1 1 1 1 1 C C C C C C C C C C C C C[ K ],][444342413433323124232221141312114321)()0()( 求解 :  200,h 22,200, ),( 2 0 0 )0,h 0,( h ) 33? C(0) =σ2 = = 22,由于 C(h)=C(0) - (h)=22 - (h) ∴ 当 i≠(|=22 - (|.于是, 21=2- )250(,)200 250(21)200 25023(202{[22 3 ,0150(22 223113  0100(22 22024114  00100(22 223223  00150(22 224224  02 0 0(22 224334  0(2201  C, 5 0(2203  C将以上数值代入普通克里金方程组解的矩阵形式中,得 1 1 1 1 1 22 0 0 22 22 2214321经计算得 : ==234= ][][ 21 六、克里金插值中搜索策略搜索邻域注意 1:搜索邻域中的数据点才参加估计节省 意 2:参与计算的数据点不能太多,否则计算太慢一般软件中都 内置或可选 最大的 数据点数目(与待估点最近的数据点),如 10。注意 3:防止数据丛聚带来的数据代表性不强井眼 井眼垂向数据太密,若待估点与该井近,则可能忽视邻井数据八分搜寻 ,保证各象限均有代表数据若搜寻范围无数据,则应用边际概率。里金插值方法简单克里金 (普通克里金 (泛克里金 (协同克里金 ( 贝叶斯克里金( 示克里金 (一、简单克里金 ((有克里金估计都应用线性回归算法,形式为:)]()()[(][1*  求取权系数的克里金方程组的非平稳形式,2,1( ),,(),()(  求 (n+1)个 m(u), 求 (n+1)×(n+1) 个 C(u,u)二阶平稳假设E[Z(u)] = E[Z(u+h)] = m(常数 )C(u,u+h) = C(h)  )(1)()()( 简单克里金估计的平稳形式:)]()()[(][1*  E[Z(u)] = E[Z(u+h)] = m(常数 )应用条件:随机函数二阶平稳随机函数的期望值 知不能用于具有局部趋势的情况,2,1( ),()()(  简单克里金方程组的平稳形式:,2,1( ),,(),()(  C(u,u+h) = C(h) (   )(二、普通克里金 (     ,,1)()( 应用要求:随机函数二阶平稳或符合内蕴假设随机函数的期望值 通克里金相当于在每一个位置 u,重新估计 m。由于普通克里金估计常使用滑动数据邻域,相当于均值 Z*(u),此时,实际上是一种非平稳算法,对应于变化的均值和平稳的协方差。三、泛克里金 (    uu 非平稳随机函数的漂移函数 (简称为漂移或趋势)()()( 随机函数 = 趋势 + 残差区域化变量 Z(X)是非平稳的,即 E[Z(x)]=m(x) a (具有趋势的克里金m a fk ) ( )u u0用光滑的确定性函数来模拟,或用拟合方法趋势函数一维的线性趋势 m a a x( )u  0 1二维的二次趋势 : m a a x a y a x a y a )u      0 1 2 3 2 4 2 5R( )u 用均值为 0、协方差函数为 的平稳随机函数来模拟。)( )( ) ( ) ( )u u u    1泛克里金估计值 : 残差      ( )( )( ) ( ) ( ) ( ) ( ),, , ,( ) ( ) ( ), , , ,f f k Ku u u u u u uu u u      1 011 20 1  为权值( )KTk ( )u 是与 (K+1)个权值的限制条件相对应的 (K+1)个拉格朗日参数泛克里金方程组 )h 为残差协方差函数具有外部漂移的克里金(an ()()]([ 10 Z ( )( ) ( ) ( )u u u    1估计值当 K = 1时,线性趋势函数为m a fk ) ( )u u0趋势函数可理解为二级变量( 1)外部变量必须在空间光滑地变化,否则可能导致 2)在主变量的所有数据点 u处和待估计的位置 部变量都必须是已知的。n 1)(101)()()()(1)(,,1)()()()()()(   克里金 方程组: 可理解为地震数据(如深度)( K=0时,?)利用几个变量之间的空间相关性,对其中的一个或几个变量进行空间估计,尤其适用于被估计变量的观察数据较少的情况 。协同克里金估计值(初始变量和二级变量)四、 协同克里金( 处的估计值;*,1 ,1 ,1 及 m ,,1           01,,2,1,,,2,1,1102110111协同克里金方程组传统普通协克里金标准化普通协克里金)(110*YX  111  E[x(u)] E[y(u)]为协同克里金的简化形式,即如果二级变量密集取样时,只保留与估计点同位的二级变量。)()()()()(1  对应的协同克里金方程组只要求知道Z - 协方差函数 以及 协方差函数h))()(  )0(/)0()0( (同位两种数据的相关系数)(方差函数)同位协同克里金叶斯克里金( 1987)把线性贝叶斯理论用于克里金估计技术,提出了贝叶斯克里金估计技术。他构想了一个模型,把用于空间估计的数据分为两类:观察数据: 是指那些精度比较高,但数量比较少的数据猜测数据: 是指那些精度比较低,但分布广泛的数据在观测数据比较多的地方,估计结果主要受观测数据的影响;在观测数据比较少的地方,则主要受猜测数据的影响。显然,井数据和地震数据的关系符合贝叶斯估计中观测数据和猜测数据的关系。设 Z(x), x∈ A,是观察数据的区域化变量。设 M(x), x∈ A,是猜测数据的区域化变量。Z*(= E[Z( = M(E[M(x)] = μM(x) , x∈ AM(x)是对 Z(x)的一种猜测,误差为 设已得到 Z(x), x∈ 观察值{Z( i=1,2,…,N} 。定义一个新的随机函数:ZT(x)= Z(x)x), x∈ = Z(x), i=1,2,…,NZ(贝叶斯克里金估计量为* )()()()( 个观察值有(相当于误差 (误差的随机函数)基于无偏性和估计方差最小两个条件:         m 0*0        1,, 00|1|,2,1 Z(x' ,x" ) = Z| M(x' )+ M(x' ,x" )  将数据按照不同的门槛值编码为 1或 0的过程。对于模拟目标区内的每一类相,当它出现于某一位置时,指示变量为 1,否则为 0。A (100) B (010)A (100) C(001)类型变量的指示变换: 01u 属于范畴 A 其它六、 指示克里金 (982年由 A·G·尔奈耳 )教授提出(00111) (00001)(01111) (00011)首先将连续变量截断为类型变量,然后进行指示变换。如: z = 10, 15, 20, 25, 30    续变量的指示变换设沿空间某一方向,在间距为 对样品点处观测了 Z(Z(xα+h)的值 (α=1,2,
展开阅读全文
  石油文库所有资源均是用户自行上传分享,仅供网友学习交流,未经上传用户书面授权,请勿作他用。
0条评论

还可以输入200字符

暂无评论,赶快抢占沙发吧。

关于本文
本文标题:2 kriging
链接地址:http://www.oilwenku.com/p-20327.html
关于我们 - 网站声明 - 网站地图 - 资源地图 - 友情链接 - 网站客服客服 - 联系我们
copyright@ 2016-2020 石油文库网站版权所有
经营许可证编号:川B2-20120048,ICP备案号:蜀ICP备11026253号-10号
收起
展开