指示克立格估值计算的理论和方法

指示克立格估值计算的理论和方法

近十年来，随着人口的增加和经济的发展，许多地区的地下水开发利用程度日益严重，导致了一系列区域地质环境问题。地下水数值模拟作为一种常规研究工具,规模和范围也日趋扩大。此类区域性模型应用中的一个突出问题是勘探资料特别是含水层的渗透性参数资料分布不均或缺乏,在通常的地下水数值模拟中只能用简单的参数分区来描述渗透性的非均质性。由于分区过大而忽略分区内部参数的差异性,导致模拟误差的出现。寻求一种尽可能利用有限的勘探资料,对未知区域内的含水层参数进行合理估值的方法,是目前大区域地下水流模拟中的关键性问题。地质统计学中的克立格方法是线性无偏最优的估值方法。20世纪70年代末之后,克立格方法在水文地质学领域得到了广泛应用。Delhomme(1979)用克立格法估计法国下诺曼底含水层的导水系数并绘制了导水系数的最优等值线图;Dagan等(1982)用相同方法绘制了意大利威尼斯市地下水头分布的等水头线并以此作为研究地面沉降的依据;Neuman和Cliftun(1982)用克立格法给出了美国亚利桑拉Avra河谷含水层渗透系数的估计;Ahmed等(1987)也结合实例对使用克立格法求渗透系数与其他几种方法进行了比较。由于自然界大量的数据样本不服从正态分布,用高斯克立格法进行估值比较困难,为此Journel于1983年创立了指示克立格法,它不要求数据总体服从正态分布,且不受特异值的影响。指示克立格法自创立以来,在矿产(Fytas和Chaouai,1990)、土壤(Goovaerts,Webster和Dubois,1997)、土地管理(傅佩红,2004)、石油(Data-Gupta等,1999和Hohn等,1994)等领域有广泛应用。然而,使用指示克立格法对第四纪含水层渗透性空间分布进行估值国内外还鲜有报道。本文以华北某第四系地层为例,以指示变异函数为结构分析工具,运用指示克立格法对该区域含水层渗透性的空间分布进行三维估计,结果表明指示克立格法能很好地描述第四系含水层渗透性的空间分布规律。1拉格朗日函数中的克立格估计在研究含水层渗透性时,把含水层的渗透系数作为区域化变量,记作Z(u),它同时具有随机性和结构性的特点。建立变量Z(u)的模型时,设渗透系数在区域范围内的变化为二阶平稳过程,即其均值函数和协方差函数存在且平稳,则所有取样点所形成的分布可以代替总体分布,而区域内任一点处变量取值的概率分布服从这一总体。作所有取样点的累积频率曲线,将累积频率曲线由小到大截成K+1段,则有K个截断值(阈值)Zk,k=1,2,…,K,通常以其十分位数作为阈值(如果是类型变量,则以其类型作为阈值)。为求出位置u处的渗透系数Zu小于等于阈值Zk的概率,引进二值指示函数,其定义为Ι(u‚Ζk)={1当位置u处的渗透系数Ζu小于等于阈值Ζk时0其他(1)I(u‚Zk)=⎧⎩⎨⎪⎪1当位置u处的渗透系数Zu小于等于阈值Zk时0其他(1)指示变异函数可描述为γΙ(h)=1Ν(h)Ν(h)∑α=1[i(uα+h,Ζk)-i(uα,Ζk)]2(2)γI(h)=1N(h)∑α=1N(h)[i(uα+h,Zk)−i(uα,Zk)]2(2)它与协方差函数有如下关系:γI(h)=CI(0)-CI(h),CI(0)=Var[I(u,z)]为方差。在u处出现渗透系数Zu小于等于阈值Zk的数学期望是:E{Ι(u,Ζk)}=1⋅Ρrob{Ζ(u)≤Ζk}+0⋅Ρrob{Ζ(u)>Ζk}=Ρrob{Ζ(u)≤Ζk}=F(Ζk)(3)E{I(u,Zk)}=1⋅Prob{Z(u)≤Zk}+0⋅Prob{Z(u)>Zk}=Prob{Z(u)≤Zk}=F(Zk)(3)上式表明在u处出现渗透系数Zu小于等于阈值Zk的数学期望等于它出现的概率。所以,指示克立格法估计量能够用于估计特征出现的概率,这个估计值可通过实测渗透系数的线性函数表示:F(u,Ζk)*=E{Ι(u,Ζk)ㄧ(n)}*=Ρrob*{Ζ(u)≤Ζkㄧ(n)}=n∑α=1λα(u,Ζk)Ι(uα,Ζk)(4)F(u,Zk)∗=E{I(u,Zk)ㄧ(n)}∗=Prob∗{Z(u)≤Zkㄧ(n)}=∑α=1nλα(u,Zk)I(uα,Zk)(4)式中F(u,Zk)*为渗透系数在位置u处真实概率F(u,Zk)的克立格估计;uα是第α个实测渗透系数的位置,α=1,2,…,n;λα(u,Zk)是相应的截断值Zk的简单克立格法权重。要确定每个权,必须使估计方差达到最小,即:min{E[F(u,Ζk)*-F(u,Ζk)]2}min{E[F(u,Zk)∗−F(u,Zk)]2}且受无偏性条件n∑α=1λα(u,Ζk)=1的限制。构造拉格朗日函数L(λ1,…,λn,μ),μ为拉格朗日乘数。用L对每一个未知的权系数和拉格朗日乘数求偏导并使之为0,得正规方程组:{n∑β=1λβ(u,Ζk)CΙ(uβ-uα,Ζk)+μ=CΙ(u-uα,Ζk),α=1‚2‚⋯‚nn∑α=1λα(u,Ζk)=1(5)解此方程组可求得诸λα(u,Zk),代入(4)式即可得到概率F(u,Zk)的克立格估值,它就是待估点处阈值为Zk的条件累积分布概率。这里的CI(uα-uβ,Z)为指示协方差函数。因渗透系数为二阶平稳,指示协方差函数和变异函数与其位置无关,只与位置的增量h有关,故有:CΙ(uα-uβ,Ζ)=Cov{Ι(uα,Ζ),Ι(uβ,Ζ)}=CΙ(h),|α-β|=h(6)相应的最小估计方差为σ2(u)=C(0)-n∑β=1λβ(u,Ζk)CΙ(u-uβ,Ζk)(7)指示克立格法过程就是重复一系列K个截断值Zk,k=1,2,…,K的过程,对每一个指示函数都同样地使用上述的指示克立格算法。通过(4)式可以建立条件累积分布函数代表未采样值Z(u)的不确定性概率模型。估计累计分布函数(cdf)实质上是在每个估计位置或单元上综合所有的指示估计。图1是用多重指示克立格法估计出的条件cdf的示意图。应指出的是条件cdf的估计需要K个指示函数的K个变差函数模型。通过对区域化变量Z的多个截断值Zk处[i(u,Z)]*分别估值,得到待估点处的各渗透系数截断值所对应的条件累积分布概率[Prob{Z(u)≤Zk|(n)}]*。然后采用E型估计可最终得到待估点处的渗透系数估计值Z(u),即:[Ζ(u)]*=∫+∞-∞ΖdF(u,Ζ|(n))≈Κ+1∑k=1Ζk[F(u,Ζk)*-F(u,Ζk-1)](8)对每个待估点或单元分别采用以上步骤,就可得到所有单元上的渗透系数克立格估值。2研究区地质背景及含水层渗透率华北平原位于燕山山脉以南,淮河以北,太行山及豫西山地以东,东临渤海,由滦河、海河、黄河冲积而成,是一个地势坦荡、第四系深厚的全国最辽阔的冲洪积平原。华北平原自山前向滨海分为冲洪积扇区、洪泛平原区、滨海平原区,地层结构分区明显受物源区控制,平面上展示近山前是滚动组分占优势的沉积物,稍离山前地带即以跳跃组分为主的沉积物,再远离山前地带是跳跃、悬移物质组成的沉积物,其间以泥质为绝对优势的沉积物在洼地沉积下来。沉积物的粒度级配变化随离物源区渐远,水动力减弱而变细。表现在渗透性上由山前向滨海逐渐变低。垂向上,第四系地层结构与大型河流冲积扇发育历史、构造、气候背景等息息相关。平原沉积物主要来自河流搬运,出山河流可以有平行、汇聚、分散、曲线、杂乱、旋转、倒流等多种形式,造成河流摆动非常复杂,因而摆动形成的砂体在时空展布上的规律也非常复杂。本次选择华北平原中部钻孔控制程度较高的区域为研究区,其西面紧靠山前、东面临近滨海,范围为250km×200km,该范围内共有各类钻孔192个。为研究含水层的渗透性,将含水层的渗透率作为区域化变量。考虑到钻孔中水文地质钻孔较少,含水层渗透率的直接数据不多,这里以钻孔岩性为依据,参考当地岩性渗透性的经验值,并结合水文地质钻孔的数据给出各岩性的渗透系数。因克立格估值要求估计区域规则,且最好是区域化变量的最大变化方向与坐标轴方向一致,因此将坐标系统作一变换:先旋转坐标系使X轴与区域化变量的最大变化方向一致,然后平移坐标轴使区域的左下角与原点重合。为简化计算,X、Y轴的坐标值分别缩小1000倍,以km为单位。Z轴方向上,为消除地形起伏的影响,将高程转化为埋深。此次研究埋深为200m。钻孔揭露的不同岩层厚薄差异悬殊给变异函数的搜索造成困难,因此在钻孔方向上以10m为单位统一承载,以岩层厚度为权进行调和平均计算各承载上的渗透系数。3数据总体特征为右偏态分布,基于微织构在进行坐标转换和渗透系数数据整理后,首先进行传统统计学分析,作出渗透系数的直方图,考察其是否服从正态分布或对数正态分布,并求出各分位数。如果服从正态分布或对数正态分布,则可以进行高斯克立格估值,否则可用指示克立格进行估值,因为指示克立格法是非参数估值方法,它不要求研究的数据总体服从何种分布。其直方图如图2所示。本次分析共有3592个渗透系数数据,最大值为150m/d,最小值为0.001m/d,均值为1.969m/d,中位数为0.031m/d,数据总体呈右偏态分布,极端大值将算术平均数拉向极端大值一方,极端大值对均值贡献大。离散系数为5.701,说明变异值的离散程度极大,平均数代表性极差。数据总体不服从正态分布或对数正态分布。4实验变异函数拟合取中位数0.03,和1.0、10.0作为阈值,分X、Y、Z三个方向计算各阈值的指示变异函数。滞后距X、Y方向为10,Z方向为20。由于钻孔位置分布不规则,在搜索数据对时以滞后距的一半作为距离容限,以22.5°作为方向容限,计算实验变异函数。将所得到的实验变异函数以球状模型进行拟合,即:γ(h)=C0+γΙ(h)=C0+CΙ(3h2a-h32a3)(9)其中:C0为块金常数,CI为拱高,C0+CI为基台值,a为变程。限于篇幅只给出域值为1.0实验变异函数拟合图(图3)(实线为实验变异函数,虚线为拟合的理论变异函数)。三个阈值指示变异函数拟合的参数见表1。参数a(变程)反映区域化变量的相关程度,变程越大,渗透系数的相关性越好,越小则相关性越差。块金常数(C0)反映可测范围内原点处的性状,块金常数越大,渗透系数即使在很短矩离内其差异性也很大。由图表可以看出,不同阈值不同方向其参数各不相同。X轴方向变程约30km,Y轴方向约为20km,Z轴方向为30m,即X轴变程大于Y轴变程,当然远大于Z轴方向,反映出X轴方向渗透性的连续性较Y轴方向好,且远远大于Z轴方向。各个阈值各个方向的块金常数都较大,表明渗透性在相邻很短的距离内其变异性就很大。5渗透系数估计值的估计克立格估值为线性最优无偏估计,即估计值与真值的均值相同、方差最小,指示克立格法同样也遵循这一原则。将研究区250km×200km×200km的三维区域划分为50×40×20的网格,每个网格大小为5km×5km×10km,共有40000个矩形格子。对每个格子进行估值。根据各阈值指示变异函数拟合参数,由公式(9)得到各距离h处的γ(h)和CI(h),再代入公式(3)可得到[i(u,Z)]*的普通克立格估值,此为各阈值在每个格子上的条件累积分布概率[Prob{Z(u)≤Zk|(n)}]*。最后利用公式(7)进行E型估计即可得到各个格子的渗透系数估计值Z(u)。图4是利用斯坦福大学Gslib软件制作的研究区渗透性分布的水平和垂向两个方向的切面,它们分别是Z轴方向的第20个切面、Y轴方向的第20个切面。由于Z轴方向坐标单位比X轴、Y轴方向小3个数量级,垂直方向的切面图所反映的结构特征不甚明显。从切面图可以看出,水平方向上,近山前地区地层渗透性最大,约为80m/d,向滨海区逐渐过渡变小,约为0.1m/d,在东部有一渗透性异常小的区域约为0.01m/d,是为以泥质占绝对优势的沉积物在洼地沉积;在垂直方向上,渗透性变化不明显,浅部比深部略好,这与该区域的地层结构水平分布及其物性组成特征是一致的。克立格方法不仅能够给出待估点处的估计值,同时也能够给出估计方差,它是度量估计精度的一个很好的尺度。图5为渗透性估计方差等值线图,对估计方差进行标准化,使其在0~1范围之间,图中圆点为钻孔的位置。由图中可以看出,有钻孔控制的地方,其估计方差就小,估计精度高;钻孔密度小的地方,则估计方差相对就大。局部地方如图的西部估计方差较大,达到0.8以上,说明该区域地层结构复杂,渗透性变异性大,工程控制程度低。为使估计精度满足要求,应该在该区域增加工程控制。区域化变量的空间变异结构分析是否正确、结构模型是否适用于本区等需要经过交叉验证,只有经过验证证明其对本区渗透性的指示克立格法估计确实可靠且适用时才能采纳。以渗透性的空间结构模型参数为基础,对该区已知样点进行指示克立格法交叉验证,作差值频数分布直方图(图6)。由图6可以看出,差值基本服从均值为0的正态分布,变化范围在(-14.56,8.67)内,均值为0.0014;差值在(-1,1)区间内的约占总数的67.2%,标准差为1.89。上述统计结果表明,用该空间变异结构模型来对研究区渗透性进行指示克立格法估计,其结果是可信的。6含水层渗