CN104463924A - 基于散乱点高程采样数据的数字高程地形模型生成方法 - Google Patents

Abstract

本发明公开了一种基于散乱高程采样数据的数字高程地形模型(DEM)生成方法。该方法从一组平面散乱采样点及其所对应的高程数据出发，首先确定DEM模型的采样区域Ω，通过为每个采样点定义邻近点集、支撑半径和支撑区域，进而为各采样点分别定义一个局部逼近函数和一个局部支撑函数；对于区域Ω内的任意一点，搜索支撑区域覆盖该点的各个采样点，并将它们的局部逼近函数用相应的局部支撑函数混合起来，定义覆盖区域Ω的全局高程函数；最后用全局高程函数计算平面上任意布局的DEM网格结点的高程值，生成所需的DEM模型。本发明方法能够过滤采样数据中的噪声，具有较好的内插能力，能够处理任意规模的散乱点采样数据。

Description

基于散乱点高程采样数据的数字高程地形模型生成方法

技术领域

本发明涉及一种基于散乱点高程采样数据的数字高程地形模型(DEM)生成方法，属于数字化地形建模技术领域，尤其涉及一种由散乱点高程采样数据构造数字高程地形模型(DEM)的方法。

背景技术

建立DEM的方法有多种：从数据源及采集方式上讲有：(1)直接由地形测量得到，例如用GPS、全站仪、野外测量等手段获取；(2)根据航空或航天影像，通过摄影测量途径获取，如立体坐标仪观测及空三加密法、解析测图、数字摄影测量等等；(3)从现有地形图上采集，如格网读点法、数字化仪手扶跟踪及扫描仪半自动采集，然后通过内插生成DEM等方法。DEM内插方法很多，主要有整体内插、分块内插和逐点内插三种。整体内插(参见：周兴华等,DEM内插方法与精度评定,测绘科学,2005,30(5):86-88)在整个区域内用一个数学函数式来表达地形曲面，该方法的拟合模型是由研究区内所有采样点的观测值(即采样点的高程值)建立的。该方法得到的整体内插函数保凸性较差，不易得到稳定的数值解，多项式系数的物理意义不明显，解算速度慢且对计算机容量要求较高，难以提供内插区域的局部地形特征。分块内插(参见：陈国良等,一种由等高线模型生成规则格网模型的算法,中国图象图形学报,2007,12(6):1110-1113)将采样点所在空间分成若干大小相同的块，对各分块使用不同的插值函数。其中，线性内插和双线性多项式内插是两种典型的分块内插方法。线性内插法算法简单、易于实现、运算速度快。但由于此内插法的原理采用的是选择邻近点内插，而不关心邻近点的特性，如果邻近点是极值点，那么内插出的结果与实际情况的出入会较大。双线性内插法得到的模型表面比较平滑并且精度较高，缺点是地形的特征点线面被网格平均化和平滑化，极值容易丢失。逐点内插(参见：李世平等,数字高程模型的建立与分析应用,辽宁工程技术大学学报(自然科学版),2008,5(27):31-33)以待插点为中心，定义一个局部函数去拟合周围的数据点，数据点的范围随待插位置的变化而变化，因此又称为移动拟合法。逐点内插法对未知点高程值的获得只考虑距离因素，不考虑距离之外的其它空间因素，也不关心未知点与周围其它点的拓扑关系和固有联系，该点的高程值过多地受其邻近点高程值的影响，因此采样点坐标或采样点所在的坐标系的微小变化都会使选点结果差别很大，结果可能造成数字高程模型表面的不连续。

目前，DEM模型的拓扑结构有规则网络结构和不规则三角网(Triangular IrregularNetwork，简称TIN)结构两种。常用的DEM生成方法是：先在采样平面上构造网格，然后再在所构造网格的基础上，通过线性和双线性内插建立DEM。用规则方格网高程数据记录地表起伏信息的优点是：采样点的位置信息可隐含，无需全部作为原始数据存储，后继的数据处理也比较方便，其缺点是数据采集量大，因为网格点不是特征点，一些微地形可能没有被记录。TIN结构数据的优点是：能以不同层次的分辨率来描述地表形态。与规则格网数据模型相比，TIN模型在某一特定分辨率下能用更少的空间和时间更精确地表示更加复杂的地形表面。

发明内容

本发明的目的是提供一种基于散乱点高程采样数据的数字高程地形模型(DEM)生成方法，该方法能够过滤采样数据中的噪声，具有较好的内插能力，能够处理任意规模的散乱点采样数据。

为实现以上目的，本发明采用如下技术方案：

一种基于散乱点高程采样数据的数字高程地形模型(DEM)生成方法，其具体步骤是：

a.输入平面采样点集合及相应的高程值集合：输入由n个平面采样点组成的采样点集合，以及各个采样点处的高程值采样数据；

b.确定数字高程地形模型(DEM)的采样区域：遍历平面采样点集内的各点，确定各采样点在平面上的分布范围；确定将要生成的数字地形模型(DEM)的采样区域；

c.确定任意采样点的邻点及支撑半径：对于每个采样点，计算与其距离最近的若干个邻点；计算采样点到其各邻点的最大距离，并以此距离来定义其支撑半径；

d.构造各采样点的局部逼近函数：对于每个采样点，由其自身及其邻点、以及相应的高程采样值，构造该采样点的局部高程逼近函数；

e.确定各个局部逼近函数的支撑函数：对于每个采样点，用该采样点的支撑半径为其定义支撑区域；然后构造定义于支撑区域上的支撑函数；

f.构造全局高程计算函数：根据支撑函数的局部覆盖特性，将相应局部逼近函数混合为一个全局高程计算函数，使得采样区域内的任一点均可由该函数计算出高程值；

g.构造DEM数据：将采样区域分割为适当的网格拓扑结构，通过全局高程函数计算网格各顶点处的高程值，生成数字化高程模型的DEM数据。

本发明的方法考虑到了更多的采样数据之间的联系，顾及到了距离、位置、范围等因素，同时可以有效地过滤采样高程数据中的噪声，既能生成规则网络结构的DEM模型，也能生成不规则三角网(TIN)结构的DEM模型。此外，该方法能够具有较好的内插能力，能够处理任意规模的散乱点采样数据。

附图说明

图1是本发明的方法流程图。

图2是采样点的邻点与支撑半径示意图，其中，采样点p₀共有8个邻点p_i(i＝1,2,Λ8)。

图3是支撑函数示例，其中，(a)非均匀衰减的类高斯函数、(b)均匀衰减的线性函数，随着自变量从0变化到R，函数值均由1衰减到0。

图4是支撑区域覆盖关系示意图，其中，p₁,p₂,p₃三个采样点的支撑区域覆盖点p。

图5是DEM的网格拓扑结构示意图，其中，(a)是规则四边形网格、(b)均匀三角网格、(c)具有细节层次的三角网格、(d)具有任意拓扑结构的三角网格。

具体实施方式

下面结合附图和实施例来进一步说明本发明。

如图1所示，基于散乱点高程采样数据的数字高程地形模型(DEM)生成方法，其具体步骤是：

a.输入平面采样点集合及相应的高程值集合：输入由n个平面采样点组成的采样点集合，以及各个采样点处的高程值采样数据。

b.确定数字高程地形模型(DEM)的采样区域：遍历平面采样点集内的各点，确定各采样点在平面上分布的范围，进而在该范围基础上确定将要生成的数字地形模型(DEM)的采样区域。

c.确定任意采样点的邻点及支撑半径：对于每个采样点，查找并记录与其距离最近的若干个邻点(不同采样点的邻近点数可以不同)；计算其到各个邻点的最大距离d，并以此距离进一步定义该采样点的支撑半径R，并使R≥d。如图2所示，采样点p₀的邻近点共有8个，即p_i(i＝1,2,Λ8)，其支撑半径取为采样点p₀到采样点p₆的距离。

d.构造各采样点的局部逼近函数：对于每个采样点，由其自身及其邻点、以及相应的高程采样值，构造该采样点的局部高程逼近函数，使该函数在上述各采样点处的值逼近相应的高程采样值。

e.确定各个局部逼近函数的支撑函数：对于每个采样点，用该采样点的支撑半径定义一个平面区域，称其为该采样点的支撑区域；然后构造定义于支撑区域上的支撑函数，支撑函数是一个径向单调非増函数，其值在该采样点处取为1，并沿径向外逐步衰减，到达支撑区域边界时取值为0。图3给出了两个支撑函数的示例，其中采样点到支撑区域边界的距离为R。

f.构造全局高程计算函数：根据支撑函数的局部覆盖特性，将相应局部逼近函数混合为一个全局高程计算函数，使得采样区域内的任一点均可由其计算出高程值。如图4所示，设p为DEM采样区域内的任一点，并且有且只有p₁、p₂、p₃三个采样点的支撑区域覆盖着点p，那么p点处的高程值由p₁、p₂、p₃的局部逼近函数与局部支撑函数共同确定。

g.构造DEM数据：根据需要，将采样区域分割为适当的网格拓扑结构(如图5所示，采样区域为一平面矩形区域；在区域上依据构造了规则四边形网格、均匀三角网格、具有细节层次的三角网格、以及具有任意拓扑结构的三角网格)，通过全局高程函数计算网格各顶点处的高程值，生成数字化高程模型的DEM数据。

实施例

本发明方法的输入和输出分别是：

输入：一组平面散乱采样点P＝{(x₁,y₁),(x₂,y₂),Λ,(x_n,y_n)}及其所对应的高程值H＝{h₁,h₂,Λ,h_n}，其中：采样点(x_i,y_i)处的高程值为h_i。

输出：在采样点分布的确定区域范围内，由采样点高程值所确定的数字化高程地形模型DEM。

具体的实施步骤为：

1.输入平面采样点集合及相应的高程值集合：输入由n个平面采样点组成的采样点集合P＝{p_i(x_i,y_i)∈R²,i＝1,Λ,n}，以及这n个采样点处的高程值集合H＝{h_i∈R,i＝1,Λ,n}，并令采样点p_i(x_i,y_i)处的高程值为h_i(i＝1,Λ,n)。

2.确定数字化地形DEM模型的采样区域：遍历平面点集P内的各个采样点，求出各点在橫坐标和纵坐标上分布范围的最大值和最小值，即

$a_1=\underset{i=1,\mathrm{\Λ},n}{\min }{x}_i$ 和 $a_2=\underset{i=1,\mathrm{\Λ},n}{\max }{x}_i---\left(1\right)$

$b_1=\underset{i=1,\mathrm{\Λ},n}{\min }{y}_i$ 和 $b_2=\underset{i=1,\mathrm{\Λ},n}{\max }{y}_i---\left(2\right)$

确定数字化地形DEM模型的采样区域为平面矩形Ω＝[a₁-Δ,a₂+Δ]×[b₁-Δ,b₂+Δ]。这里，取

$\mathrm{\Δ}=\frac{1}{20}\max \{a_2-a_1,b_2-b_1\}---\left(3\right)$

3.确定任意采样点的邻点及支撑半径：运用基于空间剖分的基本原理，将平面点集P的采样点组织在一棵Kd-树数据结构中，实现对P中采样点的快速查找(具体方法参见：朱林华等,一种节省内存的点云中K最近邻算法,兵工自动化,2008,27(7):23-26)。对于每个点p_i∈P，利用Kd-树结构查找与p_i距离最近的k个邻点(本实施例中统一取k＝15)。计算p_i到其k个邻近点的最大距离R_i：

$R_i=\underset{j=1,\mathrm{\Λ},k}{\max }\left|\right|p_i-p_j^i\left|\right|---\left(4\right)$

并定义R_i为p_i的支撑半径。

4.构造各采样点的局部逼近函数：对于P中每个点p_i∈P，由及p_i的k个邻点以及它们依次所对应的高程值及构造p_i的局部二次曲面逼近f_i(x,y)，具体方法为：

对于p_i的每个邻点(包括其自身)，构造向量

$b_j^i=\left[\begin{array}{cccccc}1& x_j^i& y_j^i& x_j^{i^2}& x_j^i{y}_j^i& y_j^{i^2}\end{array}\right]---\left(5\right)$

进而,构造矩阵

$M_j^i=b_j^{i^T}\·b_j^i=\left[\begin{array}{cccccc}1& x_j^i& y_j^i& x_j^{i^2}& x_j^i{y}_j^i& y_j^{i^2}\\ x_j^i& x_j^{i^2}& x_j^i{y}_j^i& x_j^{i^3}& x_j^{i^2}{y}_j^i& x_j^i{y}_j^{i^2}\\ y_j^i& x_j^i{y}_j^i& y_j^{i^2}& x_j^{i^2}{y}_j^i& x_j^i{y}_j^{i^2}& y_j^{i^3}\\ x_j^{i^2}& x_j^{i^3}& x_j^{i^2}{y}_j^i& x_j^{i^4}& x_j^{i^3}{y}_j^i& x_j^{i^2}{y}_j^{i^2}\\ x_j^i{y}_j^i& x_j^{i^2}{y}_j^i& x_j^i{y}_j^{i^2}& x_j^{i^3}{y}_j^i& x_j^{i^2}{y}_j^{i^2}& x_j^i{y}_j^{i^3}\\ y_j^{i^2}& x_j^i{y}_j^{i^2}& y_j^{i^3}& x_j^{i^2}{y}_j^{i^2}& x_j^i{y}_j^{i^3}& y_j^{i^4}\end{array}\right]---\left(6\right)$

及向量

$B_j^i=b_j^{i^T}\·h_j^i={\left[\begin{array}{cccccc}{h}_j^i& x_j^i{h}_j^i& y_j^i{h}_j^i& x_j^{i^2}{h}_j^i& x_j^i{y}_j^i{h}_j^i& y_j^{i^2}{h}_j^i\end{array}\right]}^T---\left(7\right)$

求解线性方程组：

$\left(\underset{j=0}{\overset{k}{\mathrm{\Σ}}}{M}_j^i\right)\·\left[\begin{array}{c}{c}_1^i\\ c_2^i\\ c_3^i\\ c_4^i\\ c_5^i\\ c_6^i\end{array}\right]=\left(\underset{j=0}{\overset{k}{\mathrm{\Σ}}}{B}_j^i\right)---\left(8\right)$

即

$\underset{j=0}{\overset{k}{\mathrm{\Σ}}}\left[\begin{array}{cccccc}1& x_j^i& y_j^i& x_j^{i^2}& x_j^i{y}_j^i& y_j^{i^2}\\ x_j^i& x_j^{i^2}& x_j^i{y}_j^i& x_j^{i^3}& x_j^{i^2}{y}_j^i& x_j^i{y}_j^{i^2}\\ y_j^i& x_j^i{y}_j^i& y_j^{i^2}& x_j^{i^2}{y}_j^i& x_j^i{y}_j^{i^2}& y_j^{i^3}\\ x_j^{i^2}& x_j^{i^3}& x_j^{i^2}{y}_j^i& x_j^{i^4}& x_j^{i^3}{y}_j^i& x_j^{i^2}{y}_j^{i^2}\\ x_j^i{y}_j^i& x_j^{i^2}{y}_j^i& x_j^i{y}_j^{i^2}& x_j^{i^3}{y}_j^i& x_j^{i^2}{y}_j^{i^2}& x_j^i{y}_j^{i^3}\\ y_j^{i^2}& x_j^i{y}_j^{i^2}& y_j^{i^3}& x_j^{i^2}{y}_j^{i^2}& x_j^i{y}_j^{i^3}& y_j^{i^4}\end{array}\right]\·\left[\begin{array}{c}{c}_1^i\\ c_2^i\\ c_3^i\\ c_4^i\\ c_5^i\\ c_6^i\end{array}\right]==\underset{j=0}{\overset{k}{\mathrm{\Σ}}}\left[\begin{array}{c}1\\ x_j^i\\ y_j^i\\ x_j^{i^2}\\ x_j^i{y}_j^i\\ y_j^{i^2}\end{array}\right]\·h_j^i---\left(9\right)$

可得则定义局部逼近函数为二元二次多项式函数：

$f_i(x,y)=c_1^i+c_2^{i}x+c_3^{i}y+c_4^i{x}^2+c_5^i\mathrm{xy}+c_6^i{y}^2---\left(10\right)$

即对于p_i近旁的任一点p(x,y)，可用函数f_i(x,y)为其估计一个高程值

5.确定各个局部逼近函数的支撑函数：对于每个采样点p_i∈P，令其支撑半径为R_i，本实施例定义以p_i为圆心、以R_i为半径的圆为点p_i的支撑区域；同时定义函数

$w_i\left(p\right)=\left\{\begin{array}{cc}{(1-\frac{\left|\right|p_i-p\left|\right|}{R_i})}^4\·(\frac{4\·\left|\right|p_i-p\left|\right|}{R_i}+1),& \left|\right|p_i-p\left|\right|\≤R_i\\ 0,& \left|\right|p_i-p\left|\right|>R_i\end{array}\right.---\left(11\right)$

为p_i的支撑函数。

6.构造全局高程计算函数：对于采样区域内的任意一点p∈Ω，搜索支撑域覆盖p的所有采样点，具体方法为：在Kd-树的辅助下，找到所有位于以p为圆心、以最大支撑半径为半径的圆内的所有采样点，再进一步求出支撑区域覆盖p的采样点。

不失一般性，设这些采样点为且它们对应的局部逼近函数和支撑函数分别为以及那么，定义p点处的全局高程计算函数为

$f^p\left(p\right)=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}(w_i^p\left(p\right)\·f_i^p\left(p\right))}{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{w}_i^p\left(p\right)}---\left(12\right)$

7.构造DEM数据：根据需要，将DEM模型采样区域Ω分割为适当的拓扑结构，如图5所示的规则四边形网格、均匀三角网格、具有细节层次的半规则细分三角网格、具有任意拓扑结构的三角网格等。平面网格上任一顶点的平面坐标都可以计算出来：令该网格顶点为p(x,y)∈Ω，由上步中的全局高程函数计算出p点处的高程值H_p＝f^p(p)，便可由各个p(x,y)及其H_p得到覆盖整个采样区域Ω的数字化高程模型的DEM数据，其定义了一块数字化地形的起伏情况。

Claims (1)

1.基于散乱点高程采样数据的数字高程地形模型生成方法，其特征在于，具体步骤是：

a.输入平面采样点集合及相应的高程值集合：输入由n个平面采样点组成的采样点集合，以及各个采样点处的高程值采样数据；

b.确定数字高程地形模型的采样区域：遍历平面采样点集内的各点，确定各采样点在平面上的分布范围；确定将要生成的数字地形模型的采样区域；

c.确定任意采样点的邻点及支撑半径：对于每个采样点，计算与其距离最近的若干个邻点；计算采样点到其各邻点的最大距离，并以此距离来定义其支撑半径；

d.构造各采样点的局部逼近函数：对于每个采样点，由其自身及其邻点、以及相应的高程采样值，构造该采样点的局部高程逼近函数；

e.确定各个局部逼近函数的支撑函数：对于每个采样点，用该采样点的支撑半径为其定义支撑区域；然后构造定义于支撑区域上的支撑函数；

f.构造全局高程计算函数：根据支撑函数的局部覆盖特性，将相应局部逼近函数混合为一个全局高程计算函数，使得采样区域内的任一点均可由该函数计算出高程值；

g.构造DEM数据：将采样区域分割为适当的网格拓扑结构，通过全局高程函数计算网格各顶点处的高程值，生成数字化高程模型的DEM数据。