数字高程模型内插

空间插值常用于将离散点的测量数据加密和扩展，以便与其他空间现象的分布模式进行比较，它包括：空间内插算法：通过已知点的数据推求同一区域其他未知点数据的计算方法；空间外推算法：通过已知区域的数据，推求其他区域数据的方法。第五讲数字高程模型内插在以下几种情况下必须作空间插值：1）现有的离散点的分辨率，像元大小或方向与所要求的不符，需要重新插值。例如将一个扫描影像（航空像片、遥感影像）从一种分辨率或方向转换到另一种分辨率或方向的影像。2）现有的连续曲面的数据模型与所需的数据模型不符，需要重新插值。如将一个连续的曲面从一种空间切分方式变为另一种空间切分方式，从TIN到栅格、栅格到TIN或矢量多边形到栅格。3）现有的数据不能完全覆盖所要求的区域范围，需要插值。

空间插值一、DEM内插概念介绍定义：根据相邻若干参考点高程值求出待定点高程值；内插是DEM研究的核心问题，其渗透在DEM生产、质量控制、精度评定和分析应用等各个环节；其总的思想是：取地形表面一系列离散点构造连续函数来描述地形表面，利用该函数可以获取地表任一点的高程，实质是插值逼近或曲面拟合；

按内插点的分布范围，可将内插分为整体内插、分块内插和逐点内插三类，按二元函数逼近数学面和参考点的关系，内插方法又可分为纯二维内插（采样点个数=多项式系数）和曲面拟合内插；二维内插要求曲面通过内插范围的全部参考点，而曲面拟合则不要求曲面严格包括参考点，但要求遵从最小二乘法则；一、DEM内插概念介绍一、DEM内插概念介绍

DEM内插分类：分块内插插值：多项式内插、双线性内插拟合：多项式、样条函数内插、多层叠加面内插、最小二乘配置法内插逐点内插加权均值法、移动拟合法、最小二乘法

整体内插高次多项式二、整体内插整体内插：就是在整个区域用一个数学函数来表达地形曲面。整体内插函数通常是高次多项式，要求地形采样点的个数大于或等于多项式的系数数目。常用的二元多项式是：整体内插的缺点：尽管选择高次多项式能使函数面更接近实际地面，但计算量大，且函数不稳定，对于高次多项式而言，微小的测量误差扰动会造成多项式参数的很大变化，从而难以得到稳定解。具体包括：1）整体内插函数保凸性较差；2）不容易得到稳定的数值解；3）多项式系数物理意义不明显，这容易导致无意义的地形起伏现象；4）解算速度慢且对计算机容量要求较高；5）不能提供内插区域的局部地形特征。二、整体内插龙格现象二、整体内插优点：整个区域上的内插函数唯一；能得到全局光滑连续的DEM；可充分反映宏观地形特征；编程简单。整体内插函数常用来揭示整个区域内的地形宏观起伏态势。在DEM内插中，整体插值方法通常不直接用于空间插值，而是用来检测不同于总趋势的最大偏离部分，在去除了宏观地物特征后，可用剩余残差来进行局部插值。三、分块内插将地形区域按一定的方法进行分块，对每一块根据地形曲面特征单独进行曲面拟合和高程内插，称为DEM分块内插。区域分块简化了地形的曲面形态，使得每一块都可用不同的曲面进行表达。一般地，可按地形结构线或规则区域进行分块，而分块大小取决于地形的复杂程度、地形采样点的密度和分布。三、分块内插不同的分块单元可用不同的内插函数，常用的内插函数有：线性内插法双线性内插法Hardy多面函数法样条函数内插法Kriging法最小二乘拟合推估法有限元法1、线性内插法使用最靠近内插点的3个已知数据点确定一个平面,继而求出内插点的高程。其函数形式为：由于线性内插法属于平面拟合,此方法不适于丘陵地区和山区;1、线性内插法线性内插法当高程点近似一条线时，有可能难以得到稳定解，此时常采用双线性内插的方法；2、双线性多项式内插是使用最靠近内插点的4个已知数据点确定出一个函数，这样,由4个已知点构成的四边形内任一点的内插高程值就能唯一确定。其函数形式为：当把4个已知点坐标代入上式，就可以确定一个唯一的多项式：

A=X-1Z双线性多项式内插常用于GridDEM中。2、双线性多项式内插假设E为坐标原点，EF=a,EC=b。过A点作EC或FD的平行线，它们与EF、CD边分别交于G、H。

(1)利用E、F、C、D四点的坐标线性内插出点G和H的高程值：

(2)利用点G和H的高程值再次进行线性内插，得到：

综合以上两式得：

3、Hardy多面函数法由美国Hardy教授于1977年提出，其基本思想是在每个数据点上同各个已知点分别建立函数关系(这种函数称为核函数，它的图形是一个规则数学曲面)，将这些规则数学曲面按一定比例叠加起来,就可拟合出任何不规则的曲面,且能达到较好的拟合效果。3、Hardy多面函数法其函数表达式为：为了简单起见,核函数一般选用形式简单的对称函数，通常有如下几种形式:3、Hardy多面函数法对Hardy多面函数法(1)的解算，主要是确定几个待定系数ki(i=1…n)。设有m个已知三维坐标的数据点，当m>n时,由(1)式构成误差方程式：3、Hardy多面函数法

关于此方法的几点讨论：如果希望对地形增加各种约束和限制，可以设计某一函数将其增加到多面叠加的函数体内，如希望在内插中考虑地面坡度信息，则可以设计具有坡度特性的函数；在DEM中，如果数据点密度较小且数据点精度很高的情况下，可采用多面函数法，但在地表特征复杂的情况下，难以确定某一特定函数严格表示地形变化，故不易采用；由于该方法在大范围计算量大，选择函数困难，故应用较少；但也有人认为：Hardy多面函数法的计算简单、快捷，但是它要求参考点是地形表面的地形特征点,否则将导致计算失败；4、样条函数内插法该方法是将某一欲插值的区域分成若干块，对每一分块定义出一个不同的多项式曲面；为了保证各分块曲面之间的光滑性，必须保证所确定的ｎ次多项式曲面与相邻分块的边界上所有(n-1)次的导数都连续；这时的ｎ次多项式就称为样条函数；这种方法属于曲面拟合范畴，对于规则网格数据，由该法可对每个点的高程重新插值；4、样条函数内插法现以三次曲面法作为分块单元加以说明；任一矩形ABCD可构成双三次曲面方程：可根据ABCD四个角点高程在x和y方向上的一阶导数和二阶混合导数,求解这个方程式;4、样条函数内插法显然，在ABCD四个角点为已知的情况下,问题的关键就是设法求得三次曲面的一阶导数和二阶混合导数；可使用导数一阶差商中数求任一网格点Ａ(i,j)的导数值:4、样条函数内插法将系数代入方程，即可根据具体要求，在该矩形格内插值任一点的高程；根据以上分析求得的曲面在相邻边上的一阶导数是连续的，因此在整个区域的曲面连结是光滑的；样条函数的理论严密,但计算量大,并且在平滑之后有时会大幅度改变已知点的高程值；样条函数保留了微地形特征，拟和时仅使用较少的点，而且保证了块之间连续，意味着这种方法可修改曲面的某一分块而不必计算整个曲面；5、最小二乘配置法最小二乘配置法广泛用于测量学科。在测量中，测量值包含三部分：与某些参数有关的值，即趋势面；趋势面不能表达的值，即系统的信号部分；观测值的偶然误差，即随机噪声。最小二乘配置法E(信号)+E(随机噪声)=0，且COV(信号，随机噪声)=0使用最小二乘法求解趋势面参数。5、最小二乘配置法在数字地面模型内插中可以使用该方法：Hi

:参考点i的实测高程值；hi:i点投影到趋势面上的点的高程值；zi

:从趋势面起算的i点的高程值；si

:实际地面与趋势面在i点的高程差；ri

:i点的测量误差；5、最小二乘配置法目标：E(zi

)=E(si

)=E(ri

)=0，误差方程：Z=S+R=H-AW5、最小二乘配置法最小二乘配置的前提，是处理对象属于平稳随机过程，但地表起伏复杂，各种地貌形态未必都符合平稳随机过程的统计规律，趋势面起算高程未必仅和距离有关，所以前提条件不符合，难以保证内插质量；趋势面和协方差函数的参数确定，需要较长的迭代过程，当收敛速度较慢时，计算量较大。四、逐点内插法分块内插的分块范围一经确定，其形状大小和位置都保持不变，凡落在分块上的待插点都用该分块上的唯一数学模型面进行内插；而逐点内插是以待插点为中心，定义一个局部函数来拟合周围的数据点，其范围随待插点位置变化而变化，故又称为移动曲面法。四、逐点内插法逐点内插法主要有：移动拟合法加权平均法Voronoi图法顾及地貌特征的逐点内插1、移动拟合法对每个待插点，选其邻近的n个数据点拟合一个多项式曲面，如：移动曲面拟合法的关键是解决两方面的问题：如何确定待插点的最小邻域范围以保证有足够的参考点？如何确定各参考点的权重1、移动拟合法邻域范围的选择：基于点数选点：以待插点为中心，选取距其最近的n个点(n≥6)来拟合多项式曲面。1、移动拟合法邻域范围的选择：基于范围选点：选取以待插点为圆心半径为R的圆内的点。由于固定R值可能影响选点数目，可考虑采用动态圆半径的方法，其思想是基于数据点的平均密度：先确定圆内数据点的数目，如R=10；再解求圆的半径：πR2=10×(A/N)，N为总点数，A为总面积。1、移动拟合法邻域范围的选择：按方位取点：以待插点为中心，在不同方位上选取距离最近的固定数目的点：以待插点为中心将平面分隔成n个扇面；在每个扇面选取距离最近的k个点数。1、移动拟合法现象：由于两点相互位置越近，其相似性越强，距离越远，则相似性越小。采样点距离待插点的距离不同，对待插点的高程插值影响程度不同。因此，可以利用权来反映这种影响性质，权函数常使用距离来度量：2、加权平均法在移动拟合法中，往往要计算复杂的误差方程组，故在实际应用中经常使用经过简化的加权平均法。在解算待定点p的高程时，使用加权平均值代替误差方程：2、加权平均法移动拟合法选点的一个缺陷：以距离为基础进行选点和定义权重，难以很好地描述空间相邻性；3、Voronoi图法

Voronoi图将平面分成N个区域，每区域包括一个点，该区域是离该点最近的点的集合；Voronoi图对散点空间的剖分是唯一的，每个voronoi区域是一个凸多边形；利用Voronoi图可以找到最佳的邻近点，也可以非常方便地定出邻近点的权；3、Voronoi图法具体的解决方案是：设插入点x后其所在的voronoi多边形为Vx,设点x的相邻点集为p1,p2,p3,pn,设pi所在的Voronoi多边形为Vp，显然点p越接近x，其权越重，否则越轻；将邻接点Vp与Vx的相交区域记为Vi(I=1,…,n),则每个点的权重为Vi,设每点的高程为Hi，则插值点x的高程可由下式求出；4、考虑地貌特征的逐点内插在拟合曲面的插值过程中，由于使用光滑曲面表达地面，因而难以反映地性线（如山脊、山谷等），从而造成插值后的地形失真；解决方法是在插值前，判断拟合面内是否有地性线穿过，对含地性线的拟合面，按地性线将其分割，直到不含地性线为止，分割后的曲面如参考点个数不够，可扩展选点范围；4、考虑地貌特征的逐点内插如图，圆形曲面有山谷线穿过，内插点分布于山谷线两侧，如果不考虑山谷线的存在，前述拟合曲面方法均不能有效地逼近地表，但如采用加权平均内插，按距离反比取权，则与位于山脊另一侧的点反而取重权，造成插值地形失真，因此必须对山脊线两侧分割处理。五、等高线内插据规则格网DEM自动绘制等高线，主要包括以下两个步骤：利用DEM的矩形格网点的高程内插出格网边上的等高线点，并将这些点按顺序排列（即等高线跟踪）；利用这些排列后的等高线点的平面坐标进行插补，即进一步加密等高线点并绘制成光滑的曲线（即等高线光滑）；1、等高线跟踪等高线跟踪主要有两种方法：(1)对每条等高线内插边排序；(2)对同一高程等高线内插出所有点，再逐一排列每条等高线点；1、等高线跟踪

1）确定等高线高程按最低点高程与最高点高程计算最低等高线高程和最高等高线高程：由上式算得各等高线的高程为：由上式算得各等高线的高程为：1、等高线跟踪2）计算状态矩阵为了记录等高线通过DEM格网的情况，可设置两个状态矩阵H(K)与V(K)序列，分别表示等高线穿过DEM格网水平边与竖直边的状态：1、等高线跟踪格网（i,j）水平边有高程为Zk的等高线通过的条件是：等高线高程介于DEM格网水平边两端点高程之间，即：格网（i,j）竖直边有高程为Zk的等高线通过的条件:1、等高线跟踪状态矩阵H(K)与V(K)的元素为：1、等高线跟踪3）搜索等高线的起点与边界相交的等高线为开曲线，不与边界相交的等高线为闭曲线，通常首先跟踪开曲线，即沿DEM的四边搜索所有状态矩阵元素为1的格网边；处理完开曲线后，再处理闭曲线，可按先列后行或先行后列的顺序搜索DEM内部格网的水平边或垂直边；1、等高线跟踪4）内插等高线点

等高线点一般采用线性内插，格网（i,j）水平边与垂直边等高线插值点的坐标计算如右所示：1、等高线跟踪5）搜索下一个等高线点对任一格网，当已知进入边时，可按该格网状态矩阵中的值判断其离去边，方法是：1、等高线跟踪如格网中有两条高程为Zk的等高线通过，即该格网的四条边均有等高线通过，此时可按上述方法搜索，也可借助格网中心点高程Zc判断离去边，Zc

值取(Zi,j+Zi+1,j+Zi,j+1+Zi+1,j+1)/4，进入边仍为IN=1：当（

Zi