插值法及其应用毕业论文

地层面的拟合是多源地质建模中最为重要的步骤。无论是通过Delaunay细分方法增加节点还是通过网格分级增加的节点，都需要进一步求取其高程值。因此，必须借助插值方法来对有限的数据点信息所形成的初始地层面进行精确、光滑处理。插值是指在根据已知的数据计算未知值的过程，结果是形成一个连续分布的数据场(Spragueetal.，2005)。目前存在的插值方法众多，但适合三维数据场表达的插值方法还是比较有限的(胡小红等，2007)，常用的有距离加权反比法(Inverse-Distance Weighting，IDW)和普通Kriging法，其中距离加权反比法属于一种确定性差值方法，Lu et al.(2008)对该方法进行扩展，根据样本的数量和分布密度特征，使其能够根据样本的特征来确定其参数的取值；而Kriging法属于一种不确定性差值方法，Jessell(2001)对其进行了深入研究，基于该方法提出一种势场(potential-field)的插值方法，能够处理存在断层的不连续数据场。这里，仅对这两种方法进行讨论。

距离加权反比法

距离反比加权法是最常用的地质数据插值方法之一。它首先由气象学家及地质学工作者提出，后由进行改进，故该方法被称为Shepard方法。

距离反比加权法的基本思想是将插值函数f(P)定义为各数据点函数fk的加权平均，它认为与待插值点距离最近的若干个已知采样点对待插值点的值贡献最大，其贡献与距离的某次幂成反比。

距离反比加权法的基本原理可用下式表示：

数字地下空间与工程三维地质建模及应用研究

式中：f(P)是待插值点P的估计值；fi是第i(i=1，…，n)个已知采样点Pi的样本值；di是第i个样本点Pi与待插值点P的广义距离；v是距离的幂，它显著影响内插的结果，它的选择标准是最小平均绝对误差。相关研究结果表明，幂越高，内插结果越具有平滑的效果(Lu et al.，2008)。在幂指数为2时，不仅能得出较满意的内插结果，而且具有容易计算的优点。在实际应用中，一般用距离平方反比法来求待估算值。

传统的距离反比加权法是非常简单和自然的，但应用于实际的地质特征插值时，却存在以下明显的缺陷：

(1)当数据点的数目非常庞大时，f(P)的计算量将变得十分巨大，计算量的庞大甚至可能导致该方法变得无法实现。

(2)该方法只考虑了从Pi到P的距离，而没有考虑其方向。事实上，只考虑距离的大小是不充分的，有的已知离散点虽然距离待插值点较近，但它对待插值点的影响可能会被其他点屏蔽掉。

(3)在已知采样点Pi的邻域内，由于di≈0计算的误差将变得非常敏感，尤其是当两个项形式占优而又符号相反时，计算的误差更是如此。

因此，在实际应用时，必须对传统的距离反比加权法加以改进。考虑到地质特征数据的空间相关性，对距离反比加权法可附加地质体结构和影响距离两个限制条件，以提高空间几何或属性数据插值的合理性和精度。具体改进如下：

(1)地质体几何结构限制条件。在空间特征插值过程中，在选择影响待估计值的原始样本数据时只选取同一地层岩性地质内的样本；不是同一地层岩性地质体内的样本，即便距离很近，也不采用。即按照地质构造划分空间单元，并只在同一地层岩性地质体内选取样本点。

(2)邻近样本点的选择条件。选择待插值点的邻近点时，可考虑三个原则：一是距离原则，即根据地质属性数据的特征给出一距离r，在该距离之外的样本点对待插值点的估算无影响；二是点数原则，即给定一数据m，以距离待插值点最近m个样本点进行估算；三是利用Voronoi图求取待插值点的邻近点。

(3)建立数据点索引表，提高待插值点周围样本点的搜索效率，从而大大减少大数据量队的计算量。

普通Kriging插值方法

设研究区域为A，区域化变量(即欲研究的物理属性变量)为 {Z(x)∈A}，x表示空间位置。Z(x)在采样点xi(i=1，2，…，n)处的属性值(或称为区域化变量的一次实现)为Z(xi)(i=1，2，…，n)，则根据普通Kriging插值原理，未采样点x0处的属性值Z(x0)估计值是n个已知采样点属性值的加权和，即：

数字地下空间与工程三维地质建模及应用研究

其中，λi(i=1，2，…，n)为待求权系数。

假设区域化变量Z(x)在整个研究区域内满足二阶平稳假设：

(1)Z(x)的数学期望存在且等于常数：E[Z(x)]=m(常数)。

(2)Z(x)的协方差Cov(xi，xj)存在且只与两点之间的相对位置有关。或满足本征假设：

(3)E[Z(xi)-Z(xj)]=0。

(4)增量的方差存在且平稳：Var[Z(xi)-Z(xj)]=E[Z(xi)-Z(xj)]2。

依据无偏性要求：E[Z*(x0)]=E[Z(x0)]。

推导可得： 。

在无偏条件下使估计方差达到最小，即：

min{ Var[Z*(x0)-2μ( (λi-1))}，其中μ为拉格朗日乘子。

可求得求解权系数λii(=1，2，…，n)的方程组：

数字地下空间与工程三维地质建模及应用研究

求出诸权系数λii(=1，2，…，n)后，就可以求出采样点x0处的属性值Z*(x0)。

上述求解权系数λii(=1，2，…，n)的方程组中协方差Cov(xi，xj)若用变异函数γ(xi，xj)表示时，形式为：

数字地下空间与工程三维地质建模及应用研究

变异函数的定义为：

数字地下空间与工程三维地质建模及应用研究

由Kringing插值所得到的方差为：

数字地下空间与工程三维地质建模及应用研究

或

数字地下空间与工程三维地质建模及应用研究

插值法中的变异函数

变异函数是Kringing插值法插值的基础。插值中需要首先确定所研究的区域化变量的变异函数。假设研究的区域为A，区域A中有一区域化变量Z(x)，它在位置xi(i=1，2，…，N)上的一次采样为Z(xi)(i=1，2，…，N)，则Z(x)的变异函数的定义为：

数字地下空间与工程三维地质建模及应用研究

一个空间变量的空间变异性是指这个变量在空间中如何随着位置的不同而变化的性质。变异函数通过其自身的结构及其各项参数从不同的角度反映空间变异性，确定变异函数的过程就是一个对空间变异性进行结构分析的过程。

设h是一个模为r=|h|，方向为a的向量，如果存在着被向量h所隔开的Nh对观测数据点，则在a方向上相应于向量h的实验变异函数γ*(h)可表示为如下形式：

数字地下空间与工程三维地质建模及应用研究

其中，z(xi+h)和z(xi)分别位于点xi+h和xi(i=1，2，…，Nh)上的观测数据。

B.变异函数理论模型

当获取实验变异函数值后，需要先选择变异函数理论模型，然后对所选择的变异函数理论模型进行参数拟合，这一过程被称为“结构分析”。

变异函数理论模型参数一般包括：变程(range，一般用a表示)、基台(sill，一般用C(0)表示)、拱高(一般用C表示)、块金常数(Nugget，一般用C0表示)，如图所示。

图 理论变差函数

变程a表示了从空间相关性状态(|h|＜a)向不存在相关性状态(|h|＞a)转变的分界线；变异函数在原点处的间断性称为块金效应，相应的常数C0= (h)称为块金常数；基台C(0)具有协方差函数：C(h)=C(0)-γ(h)的二阶平稳区域化变量Z(x)的先验方差：Var{Z(x)}=C(0)=γ(h)；拱高C为变异函数中基台C(0)与块金常数C0之差：C=C(0)-C0。

C.变异函数理论模型分类

变异函数理论模型一般分为有基台值和无基台值两大类。有基台值的变异函数理论模型包括球状模型、指数模型、高斯模型等(图)。最常用的是球状模型。球状模型公式为：

数字地下空间与工程三维地质建模及应用研究

图 变异函数中的球状模型、指数模型、高斯模型

D.变异函数对空间变异性结构的反映

变异函数作为定量描述空间变异性的一种统计学工具，通过其自身的结构及其各项参数，从不同角度反映了空间变异性结构。利用变异函数可以对空间变量的连续性、相关性、变量的影响范围、尺度效应、原点处的间断性、各向异性等要素进行描述。

E.变异函数理论模型参数拟合

变异函数理论模型参数拟合就是利用原始采样点数据或实验变异函数取值对所选定的理论模型参数以特定的方法进行估计。拟合方法一般采用手工拟合法。

手工拟合就是依据实验变异函数的取值，一方面通过观察实验变异函数图；另一方面对所研究的区域化变量进行必要的分析，采用肉眼观察来确定变异函数模型参数，并对参数反复进行交叉验证，最终确定模型参数。其拟合的大致过程如下：

(1)首先对所研究的区域化变量进行必要的结构、背景等方面的分析，结合专家经验，确定变异函数理论模型。

(2)利用实验变异函数散点图确定变异函数参数中的块金常数、基台值、变程、各向异性角度以及各向异性比值。

(3)交叉验证。