23

 

§4.6 空间数据的插值方法

    在已观测点的区域内估算未观测点的数据的过程称为内插；在已观测点的区域外估算未观测点的数据的过程称为外推。空间数据的内插和外推在GIS中使用十分普遍。数字高程模型(DEM)，也称数字地形模型(DTM)，是一种对空间起伏变化的连续表示方法。

一、空间数据的插值 

    在已观测点的区域内估算未观测点的数据的过程称为内插；在已观测点的区域外估算未观测点的数据的过程称为外推。空间数据的内插和外推在GIS中使用十分普遍。一般情况下，空间位置越靠近的点越有可能获得与实际值相似的数据，而空间位置越远的点则获得与实际值相似的数据的可能性越小。下面介绍一些常用的内插方法:边界内插、趋势面分析、局部内插和移动平均法。局部内插包括线性内插、双线性多项式内插、和双三次多项式（样条函数）内插。

二、数字高程模型(DEM)的生成

    数字高程模型(DEM)，也称数字地形模型(DTM)，是一种对空间起伏变化的连续表示方法。由于DTM隐含有地形景观的意思，所以，常用DEM，以单纯表示高程。DEM有许多用途。DEM的表示方法有拟合法、等值线、格网DEM和不规则三角网DEM(TIN)。

三、格网DEM的建立

    格网DEM的数据可直接从解析测图仪获取，下面介绍的是如何由离散点来构建格网DEM的方法。

四、不规则三角网(DEM)的建立

    所谓建立不规则三角网DEM，就是由离散数据点构建三角网，即确定哪三个数据点构成一个三角形，也称为自动联接三角网。

Visualization of a Digital Elevation Model (Conrad, 1998).

一、 空间数据的插值

    用各种方法采集的空间数据往往是按用户自己的要求获取的采样观测值，亦既数据集合是由感兴趣的区域内的随机点或规则网点上的观测值组成的。但有时用户却需要获取未观测点上的数据，而已观测点上的数据的空间分布使我们有可能从已知点的数据推算出未知点的数据值。

    在已观测点的区域内估算未观测点的数据的过程称为内插；在已观测点的区域外估算未观测点的数据的过程称为外推。

    空间数据的内插和外推在GIS中使用十分普遍。一般情况下，空间位置越靠近的点越有可能获得与实际值相似的数据，而空间位置越远的点则获得与实际值相似的数据的可能性越小。下面介绍一些常用的内插方法。

1、边界内插

    使用边界内插法时，首先要假定任何重要的变化都发生在区域的边界上，边界内的变化则是均匀的、同质的。

图4-6-1

    边界内插的方法之一是泰森多边形法。泰森多边形法的基本原理是，未知点的最佳值由最邻近的观测值产生。如图4-6-1所示。                     

    泰森多边形的生成算法见§5.7。

2、趋势面分析

    趋势面分析是一种多项式回归分析技术。多项式回归的基本思想是用多项式表示线或面，按最小二乘法原理对数据点进行拟合，拟合时假定数据点的空间坐标X、Y为独立变量，而表示特征值的Z坐标为因变量。

    当数据为一维时，可用回归线近似表示为：

    其中，a₀、a₁为多项式的系数。当n个采样点方差和为最小时，则认为线性回归方程与被拟合曲线达到了最佳配准，如图4-6-2左图所示，即：

 

  图4-6-2

    当数据以更为复杂的方式变化时，如图4-6-2右图所示。在这种情况下，需要用到二次或高次多项式：

             （二次曲线）     在GIS中，数据往往是二维的，在这种情况下，需要用到二元二次或高次多项式：

   （二次曲面）     多项式的次数并非越高越好，超过3次的多元多项式往往会导致奇异解，因此，通常使用二次多项式。     趋势面是一种平滑函数，难以正好通过原始数据点，除非数据点数和多项式的系数的个数正好相同。这就是说，多重回归中的残差属正常分布的独立误差，而且趋势面拟合产生的偏差几乎都具有一定程度的空间非相关性。

3、局部内插

    在GIS中，实际的连续空间表面很难用一种数学多项式来描述，因此，往往使用局部内插技术，即利用局部范围内的已知采样点的数据内插出未知点的数据。常用的有线性内插、双线性多项式内插、双三次多项式（样条函数）内插。

(1)、线性内插

    线性内插的多项式函数为：

    只要将内插点周围的3个数据点的数据值带入多项式，即可解算出系数a₀、a₁、a₂ 。

(2)、双线性多项式内插

    双线性多项式内插的多项式函数为：

    只要将内插点周围的4个数据点的数据值带入多项式，即可解算出系数a₀、a₁、a₂、a₃ 。

图4-6-3

    如果数据是按正方形格网点布置的（如图4-6-3），则可用简单的公式即可计算出内存点的数据值。

    设正方形的四个角点为A、B、C、D，其相应的特征值为Z_A、Z_B、Z_C、Z_D，P点相对于A点的坐标为dX、dY，则插值点的特征值Z为：

(3)、双三次多项式（样条函数）内插

    双三次多项式是一种样条函数。样条函数是一种分段函数，对于n次多项式，在边界处其n-1阶导数连续。因此，样条函数每次只用少量的数据点，故内插速度很快；样条函数通过所有的数据点，故可用于精确的内插，可以保留微地貌特征；样条函数的n-1阶导数连续，故可用于平滑处理。

    双三次多项式内插的多项式函数为：

   

    将内插点周围的16个点的数据带入多项式，可计算出所有的系数。

4、移动平均法

    在未知点X处内插变量Z的值时，最常用的方法之一是在局部范围（或称窗口）内计算个数据点的平均值。既：

    对于二维平面的移动平均法也可用相同的公式，但位置X_i应被坐标矢量X_i代替。

    窗口的大小对内插的结果有决定性的影响。小窗口将增强近距离数据的影响；大窗口将增强远距离数据的影响，减小近距离数据的影响。

    当观测点的相互位置越近，其数据的相似性越强；当观测点的相互位置越远，其数据的相似性越低。因此，在应用移动平均法时，根据采样点到内插点的距离加权计算是很自然的。这就是加权移动平均法，即：

       其中，λ_i是采样点i对应的权值，常取的形式有：

     加权平均内插的结果随使用的函数及其参数、采样点的分布、窗口的大小等的不同而变化。通常使用的采样点数为6—8点。对于不规则分布的采样点需要不断地改变窗口的大小、形状和方向，以获取一定数量的采样点。

二、数字高程模型(DEM)的生成 

1、数字高程模型(DEM)的概念

    数字高程模型(DEM)，也称数字地形模型(DTM)，是一种对空间起伏变化的连续表示方法。由于DTM隐含有地形景观的意思，所以，常用DEM，以单纯表示高程。

    尽管DEM是为了模拟地面起伏而开始发展起来的，但也可以用于模拟其它二维表面的连续高度变化，如气温、降水量等。对于一些不具有三维空间连续分布特征的地理现象，如人口密度等，从宏观上讲，也可以用DEM来表示、分析和计算。

     DEM有许多用途，例如：在民用和军用的工程项目(如道路设计)中计算挖填土石方量；为武器精确制导进行地形匹配；为军事目的显示地形景观；进行越野通视情况分析；道路设计的路线选择、地址选择；不同地形的比较和统计分析；计算坡度和坡向，绘制坡度图、晕渲图等；用于地貌分析，计算浸蚀和径流等；与专题数据，如土壤等，进行组合分析；当用其它特征(如气温等)代替高程后，还可进行人口、地下水位等的分析。

http://www.pwtal.com/terrain.htm

2、DEM的表示方法

(1)拟合法

     拟合法是指用数学方法对表面进行拟合，主要利用连续的三维函数(如富立叶级数、高次多项式等)。但对于复杂的表面，进行整体的拟合是不可行的，所以，通常采用局部拟合法。

    局部拟合法将复杂表面分成正方形的小块，或面积大致相等的不规则形状的小块，用三维数学函数对每一小块进行拟合，由于在小块的边缘，表面的坡度不一定都是连续变化的，所以应使用加权函数来保证小块接边处的匹配。

    用拟合法表示DEM虽然在地形分析中用的不多，但在其它类型的机助设计系统(如飞机、汽车等的辅助设计)中应用广泛。

(2)等值线

    等值线是地图上表示DEM的最常用方法，但并不适用于坡度计算等地形分析工作，也不适用于制作晕渲图、立体图等。

(3)格网DEM

    格网DEM是DEM的最常用的形式，其数据的组织类似于图像栅格数据，只是每个像元的值是高程值。即格网DEM是一种高程矩阵(如图4—1)。其高程数据可直接由解析立体测图仪获取，也可由规则或不规则的离散数据内插产生。

    格网DEM的优点是：数据结构简单，便于管理；有利于地形分析，以及制作立体图。其缺点是：格网点高程的内插会损失精度；格网过大会损失地形的关键特征，如山峰、洼坑、山脊等；如不改变格网的大小，不能适用于起伏程度不同的地区；地形简单地区存在大量冗余数据。

100   110   120   140   110   105   90

120   115   130   135   120   110   100

135   120   120   130   130   120   110

145   130   115   120   120   115   118

150   140   135   130   135   120   110

145   135   150   140   138   125   120

表4-6-1

(4)不规则三角网DEM(TIN)

    不规则三角网DEM直接利用原始采样点进行地形表面的重建，由连续的相互联接的三角面组成(如图4—4)，三角面的形状和大小取决于不规则分布的观测点的密度和位置。不规则三角网DEM的优点是：能充分利用地貌的特征点、线，较好地表示复杂地形；可根据不同地形，选取合适的采样点数；进行地形分析和绘制立体图也很方便。其缺点是：由于数据结构复杂，因而不便于规范化管理，难以与矢量和栅格数据进行联合分析。

图4-6-4

    通常所说的DEM即指格网DEM和不规则三角网DEM，地形分析也基于此。

Department of Defense, Groundwater Modeling System

三、格网DEM的建立

    格网DEM的数据可直接从解析测图仪获取，下面介绍的是如何由离散点来构建格网DEM的方法。

离散点构格网DEM是在原始数据呈离散分布，或原有的格网DEM密度不够时需使用的方法。其基本思路是：选择一合理的数学模型，利用已知点上的信息求出函数的待定系数，然后求算规则格网点上的高程值。

    离散点构格网DEM所采用的是内插算法，插值的方法很多，如按距离加权法、多项式内插法、样条函数内插法、多面函数法等等。大量的实验证明，由于实际地形的非平稳性，不同的内插方法对DEM的精度并无显著影响，主要取决于原始采样点的密度和分布。简单而常用的为线性内插法和双线性多项式内插法。

    线性内插的数学模型为：

    双线性多项式内插的数学模型为：

    其中，x,y为平面坐标，Z为高程，a₁、a₂、a₃、a₄为待定系数。

    只要将与插值点距离最近的三个点(对线性内插)或四个点(对双线性多项式内插)的坐标值和高程值代入方程，即可解出全部系数，然后用插值点的坐标带入方程，即可计算出该点的高程值。

    距离加权法的数学模型为：

    为了计算点(x,y)的高程Z，可取该点周围的n个点(I=1,…,n)，按上述公式计算。Z_i为第I点的高程，P_i为第i点的权值，计算公式为：

           

    其中u是一个大于0的正数，通常取1或2，是i点(x_i,y_i)到格网点(x,y)的距离，即：

         

    在构建格网DEM的内插方法中，经常需要选取与插值点距离最近的若干个点。如图4-6-5左图，如果通过计算各离散点与待插值点的距离，然后选取距离最短的若干个点的方法来选取，尽管方法正确，但计算量大，影响插值的速度。因此，可以在插值点上建立一正方形的选取框，如图4-6-5右图，这样通过简单的坐标值比较就可找出落入框内的数据点。当落入框内的数据点较多时，可缩小框的尺寸；反之，增大框的尺寸。选取框尺寸的初始值可根据图幅中原始数据点的密度来确定。设图幅的面积为A，共有N个数据点，则每点的平均面积A₀为：

A₀＝A∕N

图4-6-5

    若需要选取插值点附近的K个点，则选取框的面积应为：A₁＝K·A₀ 。

        A₁^1/2就是该选取框的边长。在实际运算时，边长可再大一些，以尽量保证落入选取框中的数据点数大于或等于K，这样在缩小选取框时，只要对原框内数据点再判断即可。

四、不规则三角网(DEM)的建立

图4-6-6

        所谓建立不规则三角网DEM，就是由离散数据点构建三角网，如图4-6-6，即确定哪三个数据点构成一个三角形，也称为自动联接三角网。即对于平面上n个离散点，其平面坐标为(x_i，y_i)i＝1，2，…，n，将其中相近的三点构成最佳三角形，使每个离散点都成为三角形的顶点。自动联接三角网的结果为所有三角形的三个顶点的标号，如：

1,   2，  8

2,   8，  3

3,   8，  7

┇

    为了获得最佳三角形，在构三角网时，应尽可能使三角形的三内角均成锐角。其基本依据是三角形余弦定理(图4-6-7)：

cosC＝(a²＋b²－c²)／2ab

图4-6-7

    在已知A、B点，即已知c边时，要选取另一个三角形顶点C时，若C角最大，则要扩展的三角形的顶点C离扩展边c的距离最短。这样可保证由相邻最近的三点构成三角形，并且保证不让某个离散点在组成三角网时被漏掉。在自动连接三角网的软件设计时，要设L和K两个变量，L记录已形成的三角形数，K记录已扩展的三角形数。为了记录结果可设三个数组t₁[]、t₂[]、t₃[]，分别存放所构成的每个三角形的顶点编号。

    构建三角网的第一步是形成第一个三角形。可选择离散点中最相近的两个点作为第一个三角形的两个顶点，并置L＝1，把这两点的编号记录在t₁[L]和t₂[L]中。三角形的第三个顶点可按余弦定理，分别检查各个离散点，取与t₁[L]和t₂[L]连线最近的点，并把编号记录到t₃[L]中。

图4-6-8

    第一个三角形形成后，置K＝1，并以该三角形的三边向外扩展三角形。所谓扩展三角形是指从三角形的每条边向外再形成三角形。

    例如，如图4—45，对于t₁[K]、t₂[K]、t₃[K]构成的三角形，需对每一边向外扩展形成新的三角形。当从t₁[K]和t₂[K]组成的边向外扩展时，显然位于t₃[K]同侧的离散点应被排除在外。这可用直线判别正负区的原理来实现。直线方程的判别式为：

F(x,y)＝y－Ax－B

   其中：A＝(y₂－y₁)／(x₂－x₁)
         B＝(y₁x₂－y₂x₁)／(x₂－x₁) ₂
        F(x,y)＞0，点位于正区
        F(x,y)＝0，点位于线上
        F(x,y)＜0，点位于负区

    把t₃[K]的坐标代入判别式计算，记录下判别式的正负，则只有当离散点的判别式的值与t₃[K]的判别式的值符号相反时，才有可能成为被扩展的点。

    在这些可能被扩展的点中，找出对扩展边张角最大的点，就是要扩展的点。

    为了避免重复与交叉，还要进行一次检查判断，即要判断新的三角形的三条边是否已被已形成的三角形用过两次。若有一条边被用过两次，则此次扩展无效；否则此次扩展有效，且L加1。

    对三角形的第二条边和第三条边采用相同的方法扩展。三角形的三条边都扩展完后，令K加1，并继续扩展下一个三角形。

    当K＝L，且第K个三角形已扩展完毕时，表明对所有已形成的三角形都已进行了扩展，则构建三角网完成。