基于多重分形理论的地表高程插值方法研究摘要：从地表高程插值的研究现状出发，将多重分形理论引入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］。
通过比较插值数据集相应点位高程的预测值和实际观测值可以评价预测结果的精度，而通过比较验证数据集中相应点位高程的预测值和实际观测值可以评价预测结果的准确性。
通常使用均方根误差来评价预测精度和准确性。
RMSE值最小，插值精度越高；反之越低。
RMSE的计算方法如下:式中，表示第i点的实际测量值，z(xi)表示第i点的插值结果，m为插值的样点数。
2结论与分析2.1 Kriging插值结果与分析克里格插值结果如图2所示，该区域西边和南部以为平原地区为主，平均咼程在2m左右，东部和北部以山丘地为主平均咼程较咼，有地方高程为大于100m,中东部较高地区峰值达到150m,总体显示地形也是较为平滑，峰与谷的显示不明显。
图中A处为研究区域高程值最大的地方，也是最能说明本研究多重分形插值方法对于局部奇异性刻画是否成功之处，本文将这里作为典型来说明。
其真实情况为2个相邻近的山峰，而此插值方法得出的结果图为一条高程带，该插值对于局部的刻画不够精确，有削峰填谷的效应，将本来高程高的峰值降低，高程低的谷值增加，说明克里格插值方法有平滑作用。
2.2 IDW插值结果与分析IDW插值方法得出的结果如图3所示，与克里格插值结果没有太大的差异，总体数据值比较相似：西边和南部以小部分为平原地区为主，平均咼程在2m左右，东部和北部以山丘地为主平均咼程较咼，有地方咼程为大于100m,最咼中东部较咼地区峰值达到150m,咼峰和低谷处较为不明显，成片的地区较多，整体看起来较为平滑。
以A 处为例，插值结果与克里格几乎一致，说明IDW插值方法也是有平滑作用，插值结果不能很好描述局部奇异性较大的地方。
以上两种插值方法的插值结果总体来看并无差异，插值结果与实测值的比较直观上也无较大不同，两种插值方法的总体精度较高，但是对局部奇异性地区的刻画均有平滑作用。
2.3多重分形插值结果与分析多重分形克里格插值方法结果图，西边和南部以为平原地区为主，平均咼程在2m左右，东部和北部以山丘地为主平均咼程较咼，有地方高程为大于100m,与克里格插值所部同的是低谷的地方更低，峰处值更高，总体显示地形较普通克里格插值方法峰与谷的显示更为明显。
从图2A处和图4A处可以直观地看出多重分形克里格插值方法对于两个山峰的描述较为精确。
多重分形IDW插值方法结果（图5）所示，西边和南部以为平原地区为主，平均高程在2m左右，东部和北部以山丘地为主平均高程较高，有地方高程为大于100m,与IDW插值结果图比较可以看出多重分形插值IDW法峰值点较多，可以明显的看出奇异性明显的地区。
图4A处和图5A处直观看起来，二者皆减弱了原本插值的平滑作用，对于两个山峰的描述较Kriging与IDW插值更精确。
2.4不同插值方法的精度评价评价指标Kriging 多重分形KrigingIDW 多重分形IDWRMSE10.16810.15010.62310.607 从上表可以看出，Kriging 插值方法的RMSE值为10.168精确程度，要高于IDW插值方法的RMSE 值10.623,多重分形IDW插值方法和多重分形Kriging插值方法的精确程度要分别高于IDW插值方法和Kriging插值方法。
评价结果说明多重分形插值的精度要略高于IDW插值方法和Kriging插值方法。
3结论通过上述的比较分析可以得出，反距离权重方法、普通克里格方法的插值结果均能较好地反映出研究区域高程空间分布的实际情况，并且总体精度较高。
多重分形插值方法与普通插值方法结果的高程分布趋势基本一致，通过对局部奇异性指数的计算，合理地强化了局部区域的估值结果，可以更好地表达整个地区高程的局部奇异性。
参考文献[1] 吉玮，曹兆华，刘春华，等.神经网络方法在地表高程预测中的应用研究[J].资源开发与市场，2010, 26(9): 803-804.[2] 包世泰，廖衍旋，胡月明，等.基于Kriging的地形高程插值[J].地理与地理信息科学，2007, 23(3): 28-32.[3] Ma ndelbrot B B. How long is the coast of Britai n? Statistical self-similrity and fractional dimension[J]. Scienee, 1967, 156(3775): 636-638.[4] 刘鹏举，赵仁亮，朱金兆，等.保持地貌特征的数字高程模型生成方法研究[J].中国矿业大学学报，2006, 35(4): 523-526.[5] 成秋明.多重分形与地质统计学方法用于勘查地球化学异常空间结构和奇异性分析[J].地球科学(中国地质大学学报),2001, 26(2): 161-166.[6] 李锰,朱令人，龙海英.不同类型地貌的各向异性分形与多重分形特征研究[J].地球学报，2003, 24(3): 237-242.[7] 曹汉强，朱光喜，李旭涛，等.多重分形及其在地形特征分析中的应用[J].北京航空航天大学学报，2004, 30(12): 1182-1185.[8] 崔灵周，李占斌，郭彦彪，等.基于分形信息维数的流域地貌形态与侵蚀产沙关系[J]. 土壤学报，2007, 3(2): 197-203.[9] 袁贺，罗问，刘付程.广义回归神经网络残余Kriging方法预测地表高程[J].安徽大学学报(自然科学版),2010, 34(5): 21-26.[10] 汤国安，杨昕.ArcGIS地理信息系统空间分析实验教程[M].北京：科学出版社，2006.[11] 邓晓斌.基于ArcGIS两种空间插值方法的比较[J].地理空间信息，2008, 6(6): 85-87.[12] 刘付程，郭衍游，闫晓波.土壤属性空间预测的广义回归神经网络方法研究[J ].淮海工学院学报(自然科学版),2008, 17 (1): 68-71.。