基于多重分形理论的地表高程插值方法研究

基于多重分形理论的地表高程插值方法研究

摘要：从地表高程插值的研究现状出发，将多重分形理论引入Kriging和IDW两种常见插值方法，并对四种插值方法进行了比较分

析。结果表明，在相同的样本数据条件下，多重分形插值的精度要高于IDW插值方法和Kriging插值方法。通过对局部奇异性指数的计算，多重分形插值方法合理地强化了局部区域的估值结果，可以更好地表达整个地区高程的局部奇异性。

关键词：多重分形理论Krig ing IDW 高程插值

随着人们对空间高程数据的质量要求越来越高，高程插值被越来越广泛得应用在高程测量中，特别是由点源数据形成连续变化的曲面来表达高程值变化往往需要进行高程数据的插值或估计。目前地形高程插值的方法主要有反距离权重方法(IDW)、克里格插值方法(Kriging) 等［1,2］，这些方法为了体现某种趋势而采用加权平均估计，将实际高程值较大的地方消减，实际高程值较小的地方增大，对数值造成一定的光滑，此时高程的局部结构特征的表示与实际有较大差距。普通插值方法无法表达地表高程局部的复杂变异，无法对地表高程有较精确的预测。Mandelbrot上世纪60年代提出的分形理论以分维数、自相似性和幕函数等为工具，研究起伏的地形、复杂的水系等具有特征标度、极不规则但具有自相似性的复杂现象［3］。多重分形理论作为分形理论的进一步发展，不仅能够用来描述分形的复杂特征，还能特征描述分形本身的几何支撑度量。对于高程、水流等诸多非均匀的分形

现象，可以使用多重分形测度或维数的连续谱来表示［4］。目前多重

分形方法已应用到地球科学的诸多领域［5-8］，但应用于高程插值方面的研究还不多见。本文将多重分形理论引入到IDW和Kriging这两种常用的地形高程插值方法中，试图对高程这一非均匀分形现象进行更精确的预测。

1研究数据与方法

1.1研究数据

本研究所使用的数据为江苏省连云港市海州区内某地28.85km2范围内抽取的1085个点作为采样点（如图1所示），每个采样点使用全站仪测得每个点的高程值，其中高程最大值为187.109米，最小值为0.914米，平均值为14.629米。测量所得总体样本高程有一定差异，符合研究需要。将其中200个样点作为插值点用于插值研究（具体分布如图1中黑点），余下的885个样点作为检验点用作检验（具体分布如图1中白点）。

1.2研究方法

1.2.1克里格插值方法

克里格插值（Kriging）又称空间局部插值法，是以变异函数理论和

结构分析为基础，在一定的区域内对区域化变量进行无偏最优估计的一种方法，是地质学的主要方法之一。克里格法是根据待插值点与临

近实测高程点的空间位置，对待插值点的高程值进行线性无偏最优估计，通过生成一个关于高程的克里格插值图来表达研究区域的原始地形。在高程分析中，克里格插值方法是建立在平稳假设的基础上，这种假设在一定程度上要求所有数据值具有相同的变异性［2］。

本文采用普通克里格(O-Kriging)方法进行插值，O-Kriging是区域化变量的线性估计，它假设数据变化成正态分布，插值过程类似于加权移动，权重值的确定来自于空间数据分析［9,10］。

1.2.2 IDW插值方法

IDW(Inverse Distanee Weighted)是一种常用而简便的空间插值方法，它以插值点与样本点间的距离为权重进行加权平均，离插值点越近的样本点赋予的权重越大［11］。

反距离权重插值的一般公式为：

式中，是估计值，是第(=1,…个样本的实测值，是与测定点有关的权重，dij是第i个采样点与估计点之间的距离，r是距离的幕，用来控制权重值随距离变化的速度，当指数增加时距离远的观测点权重值则会下

降，它显著影响内插的结果。实际应用中的r取值范围一

般取1、2、3,其中2最为常用。

反距离加权法很简单，几乎任何GIS都能提供该方法。但是其缺陷也是明显的，它不能对内插的结果作精度评价，所得结果可能会出现很大的偏差，人为难以控制。它也没有考虑数据场在空间的分布情况，如果输入样点太少或不均匀，IDW则不能很好地表达期望的表面，影响内插精度。

123多重分形插值方法

多重分形(Multifractal)是指用多个维数来描述非均匀的复杂几何

体，以此来全面刻画其特征。多重分形也称为多标度分形，是定义在分形结构上的由多个不同标度指数的概率子集组成的非均匀分形测度的集合。在地形变化复杂的地区或者相对较大的区域，高程异常的

变化不规则，用规则的曲面取代实际上不规则的似大地水准面，必然会导致大地水准面的不规则变化，导致拟合后高程异常的精度降低、误差较大。多重分形可通过高程值计算出高程局部异常的奇异性指数，通过奇异性指数得到多重分形Kriging插值和多重分形IDW插值的结果。

局部奇异性分析方法实际上是将样点高程值段的密度在分形空间中进行度量，以确定分形密度和分形维数。奇异性指数是场值随量度范围大小的变化规律众多的空间插值基于对场值的某种滑动加权平均，公式如下：

式中Q(xO, £是围绕中心点xO、半径为£的小滑动窗口。？(||xO-x||)

是对在Q (xO, £中与中心点xO相隔||xO-x|距离的任意点x的加权函数。它往往与距离呈反相关。加权函数的选择不仅与距离有关，还与

场的空间性自相关以及处理目的有关。

成秋明［5］给出的多重分形方法将滑动平均关系表达为

这里a (X0是X0点处的局部奇异性指数。可以看出，以上表达式中不仅包含了空间相关性的成分，还有度量奇异性指数。如果a(XO)=2那么通过该方法所计算的加权平均值与通常的加权平均无异。然而，当处于高程差异较大地段而且局部地区具有奇异性，a (XO)<；2该方法所得的结果将高于通常的加权平均结果；相反，当高程值较低地段，a (xO)>2该方法所得的结果将低于通常的加权平均结果。由此可见，该方法有利于加强峰谷值，对于具有奇异性的空间结构来说，传统的插值方法不能很好的估计出峰谷值。指数a 对应分形空间维数，分形维数与正常的欧氏空间维数的差Aa二N(对于二维场)即可表示分形密度与正常密度的空间维数的差异，并通过ArcGIS软件计算得到奇异性指数。

1.2.4评价方法

对于四种不同预测方法的优劣可以采用插值数据集和验证数据集进行评价［12］。通过比较插值数据集相应点位高程的预测值和实际观测值可以评价预测结果的精度，而通过比较验证数据集中相应点位高程的预测值和实际观测值可以评价预测结果的准确性。通常使用均