基于空间分布权系数的邻域选点算法

基于空间分布权系数的邻域选点算法

khig法是1951年南非采矿工程师提出的一种计算采矿储量的方法。在计算krig值时，通常不需要参与整个变程范围内的原始数据，而只计算评估点的一些相邻点。在计算每个预测点之前，需要选择每个预测点，而不是选择所有采样点。边缘点的选拔是krig法计算的一个组成部分。在这项工作中，我们提出了一种基于空间分布系数的相邻区域选择算法。在同一网格内，基于样本点和估计点之间的距离。相邻点的空间分布系数算法不仅考虑了样本点和估计点之间的距离，还考虑了样本点之间的空间分布关系。这是一种动态和静态选择相邻点的算法。使用空间分布系数方法选择相邻点，其精度优于使用普通限制搜索法选择相邻点的精度。1共标区域划分邻域选点是空间插值计算一项必不可少的内容,目前主要算法有以下几种.1)固定距离法.以估值点为中心,以固定距离为半径,搜索所有落入圆内的点,该算法简单,但当样点分布不均匀的时候会出现选择的点数过多或过少等情况;2)固定数目法.预先固定邻域点的数目,选择距离估值点最近的样本点.该算法简单,克立格方程组维数固定,但样点分布不均的时候可能会导致外推;3)象限搜索法.把整个区域划分成若干象限,每个象限设定最多选点数量,如果少于最多选点数则全部采用,如果大于最多选点数则选择最近点,这种方法是目前广泛采用的一种方法,缺点是在同一象限选择邻域点只能根据离估值点的远近进行选取;4)Delaunay-固定距离邻域选择算法.对邻域点进行三角剖分,根据邻域点所在三角型的邻近阶次选择邻域点.利用Delaunay三角划分进行邻域点搜索可以避免外推,但算法复杂计算量大.2区域选点管理算法用于空间分布系数的区域选择2.1克立格方程方程设Z(x)是区域化变量,在整个区域范围内是2阶平稳或是本征的,Z(xi)(i=1,2,…,n)是一组离散的信息样品数据,根据已知样本点的值可以求得未知点的值,即Z(x)=∑i=1nλiZ(xi)(i=1,2,⋯,n)(1)Ζ(x)=∑i=1nλiΖ(xi)(i=1,2,⋯,n)(1)在保证Z(x)是线性无偏,且方差最小的情况下可得克立格方程组为⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪∑j=1nλjγ(xi,xj)+μ=γ(xi,x)(i=1,2,⋯n)∑i=1nλj=1(2){∑j=1nλjγ(xi,xj)+μ=γ(xi,x)(i=1,2,⋯n)∑i=1nλj=1(2)式中:λ为克立格方程系数,γ为变差函数.2.2比i计算a图1为邻域点空间位置关系图,图中样本点i,j,k到估值点o距离相等,但由于克立格本身的特点i,j,k点对o点的估值贡献不同,从图中可以看出i,j点的相关程度大于k,i点的相关程度,用jo在io上投影来表达相关程度,djcos(∠ioj)的值越大相关程度越大,即dj(1-cos(∠ioj))的值越小相关程度越大.因此i点与邻域内其余点的相关程度可以用式(3)表示.P¯¯¯i=∑j=1j≠indj(1−cos(∠ioj))(3)Ρ¯i=∑j=1j≠indj(1-cos(∠ioj))(3)为排除距离的影响可以将式(3)除以∑j=1j≠indj∑j=1j≠indj得到P¯¯¯¯¯¯i=∑j=1j≠indj(1−cos(∠ioj))/∑j=1j≠indj(4)Ρ¯¯i=∑j=1j≠indj(1-cos(∠ioj))/∑j=1j≠indj(4)P¯¯¯¯¯¯Ρ¯¯i的值越大表示i点与邻域内其余点的相关程度越低,反之相关程度越高.综合考虑样本点间相关程度和距离的综合影响可得样本点空间分布权函数为Pi=d−βi/∑j=1nd−βj+θP¯¯¯¯¯¯i(5)Ρi=di-β/∑j=1ndj-β+θΡ¯¯i(5)式中:β表示距离的幂,如β为1表示距离反比,β为2表示距离平方反比;θ表示距离与相关程度的影响比例,θ越大表示距离的影响越小.2.3当前限值的计算如果估值点邻域内任意1点与过多选出参与估值计算的邻域点进行权系数计算会淹没该点与周围邻近的几个点在估值计算中的权值比重,因此在计算权值时把选出参与估值计算的部分邻域点参与空间分布权系数计算,这样可以保证每个象限内所选择的邻域点和相邻象限之间所选择的邻域点之间有比较好的空间分布关系.已知β,θ,搜索半径h,象限数量P,每个象限最多取点数m,样本点集合Z(xi)(i=1,2,…,n),参与估值计算的邻域点集合SelectPt(x),把参与空间分布权值计算的点数记为js,当前象限已经被选为参与估值计算的邻域点个数记为jp.1)计算估值点与所有样本点的距离di(i=1,2,…,n),满足di≤h的点记为集合Z(xj)(j=1,2,…,k).2)集合Z(xj)(j=1,2,…,k)按照角度大小排序,根据象限数P,计算每个象限起始邻域点序号Start(l)(l=1,2,…,P)和终止邻域点序号End(l)(l=1,2,…,P).3)计算Z(x1)所在象限编号l,计算第l象限距离估值点距离最小的点,把距离最小点选入到结果集合SelectPt(x),令jp=1,js=1.4)判断当前象限点数是否小于m,如果小于m,把象限内的点依次选入到集合SelectPt(x),同时判断js是否小于m,如果小于m,js=js+1,转入步骤6;如果当前象限点数大于等于m,转入步骤5.5)如果jp大于等于m,转入步骤6;否则,判断js是否小于m,如果小于m,根据式(5)计算当前象限未被选入参与估值计算的邻域点与集合SelectPt(x)中点的权系数,选择权系数最大的点到集合SelectPt(x),js=js+1,jp=jp+1;如果js大于等于m,根据式(5)计算当前象限未被选入参与估值计算的邻域点与集合SelectPt(x)中后m个点的权系数,选择权系数最大的点到集合SelectPt(x),jp=jp+1,重复步骤5.6)如果象限编号l大于等于P,计算终止;否则,l=l+1,jp=0,重复步骤4和5进行计算.3实验分析3.1块金常数c0的表示笔者选取某金矿矿体的品位数据,共有20个钻孔,采集样品总数518个.通过对样品金品位计算,金品位数据满足对数正态分布,金品位对数分布直方图如图2所示.金品位变异函数在空间表现各向异性,变异函数采用球状模型如式(6)γ(h)=⎧⎩⎨⎪⎪⎪⎪0C0+C(32⋅ha−12⋅h3a3)C0+Ch=00<h≤ah>a(6)γ(h)={0h=0C0+C(32⋅ha-12⋅h3a3)0<h≤aC0+Ch>a(6)式中:C0表示块金常数;C表示基台值;a表示变程.通过对金品位数据的研究,这里分别用γ(hu),γ(hv),γ(hw)表示矿体走向、倾向、厚度方向的变异函数,走向方向Cu0=0.15,Cu=0.25,au=100,倾向方向Cv0=0.1,Cv=0.3,av=80,厚度方向Cw0=0.15,Cw=0.25,aw=50.三维变异函数采用式(7)进行变异函数套合计算.γ(h)=γ(h2u)+γ(h2v)+γ(h2w)−−−−−−−−−−−−−−−−−−√(7)γ(h)=γ(hu2)+γ(hv2)+γ(hw2)(7)3.2邻域点选择算法验证本文采用了交叉验证方法验证空间分布权系数邻域选点算法的优劣,具体做法是,对每个实测点均用其周围点上的值对该点进行估值.可以得到N个实测值Z和N个估计值Z*,再求其误差平方的平均值(Z∗−Z)2¯¯¯¯¯¯¯¯¯¯¯¯¯¯(Ζ*-Ζ)2¯,此均值的大小可以衡量算法的优劣.在邻域点的选择算法上,分别采用了空间分布权系数邻域选点算法与普通象限搜索算法对估值数据进行邻域选点.为验证空间分布权系数邻域选点算法对相关估值算法的实用性,估值方法在采用对数正态克立格法的同时,还采用了距离反比法和距离平方反比法.在式(5)中的β,θ的值均取1,象限P取4,搜索半径取50m,通过计算,对数正态克立格法交叉验证结果如表1,距离反比法估值的结果见表2,距离平方反比法估值的结果见表3.通过分析上述3个表的计算数据可以看出不论选取哪种估值方法,在象限最大选点数相同的条件下,采用空间分布权系数方法选择邻域点进行估值精度高于采用普通象限搜索法选择邻域点进