顾及形态精度的地形高程采样点加密方法与流程

本发明属于地形高程样点加密(即稀疏样点转换为密集样点)的技术领域，涉及线状要素高程采样点和区域高程采样点加密，特别是提出顾及地形要素形态约束的高程采样点加密策略。

背景技术：

数字高程模型(Digital Elevation Model，简称DEM)是国家基础地理数据的核心内容。目前，我国已初步完成1:100万、1:25万、1:5万和部分地区1:1万DEM数据的建设。这些数据在国民经济、国防建设及科学研究等方面发挥着了不可或缺的重要作用，其应用领域遍布测绘、交通、军事、水利、农业、环境、资源管理、规划与旅游等众多领域。构建DEM的首要一步是获取地形数据，主要是获取高程和平面位置数据，即地形高程采样点数据，简称高程采样点。DEM是实际地形表面的再现，其对地形表达的可信程度，很大程度上取决于高程采样点的分布与密度。然而，高程采样方法无法对地表面所有点进行观测，只能获取一定数量的高程采样点数据，这些采样点反映了地形的局部或者部分特征。在这种现实情况下，若要构建高质量的DEM数据，增加高程采样点密度，即进行高程采样点加密，实现高程采样点的全局控制与局部保真，便显得尤为重要。特别是，面对诸如城市地形等具有人工与自然结合、突变与渐变交错特征的地形建模时，高程样点加密需求更为迫切。

高程采样点加密是根据已知采样点估计出未知区域的地形高程采样点，通常分为两个步骤：(1)确定加密采样点的平面位置；(2)估计加密采样点的高程值。截止目前，鲜有直接关于高程采样点加密的研究报道，相关研究主要集中在空间内插方法方面。空间内插是根据分布在待内插点周围的样本点值求出其值。在这个过程中，待内插点的平面位置是给定的。现有空间内插方法主要可分为两大类：确定性插值法和地统计学插值法。确定性插值方法是以研究区域内部的相似性或者以平滑度为基础，使用一定的数学函数，由已知样点进行待内插点值的预测。这类插值方法主要包括有最近邻法、反距离加权法、多项式法(包括线性插值、双线性插值、三次卷积插值及趋势面插值等)及径向基函数法(包括样条函数法及高次曲面函数等)。确定性插值法可在样本点处的插值结果和原样本点实际值基本一致，且计算相对简单，但其忽略了样本点间存在的空间联系，插值结果往往受样本点的影响较大。地统计学插值方法通常指克里金插值法，是以空间自相关为基础，利用原始样本数据和半方差函数的结构性，对区域化变量的待内插点进行无偏最优估值。这类插值方法不仅考虑了待内插点与邻近样点数据的空间距离关系，还考虑了各参与预测的样本点之间的位置，充分利用了样本点数据的空间分布结构特征，其插值精度相对确定性插值法而言较高。然而，地统计学插值法计算步骤较为繁琐，计算速度较慢。

空间内插方法给高程样本加密提供了一定的参考基础，但在高程内插应用上其通常主要考虑高程数值精度，忽略了地表形态精度(即模拟地面形态与实际地面形态之间的差异程度)。更何况，现有空间内插方法无法直接应用于不同需求的高程采样点加密。目前，缺少面向线状要素及区域的高程采样点加密方法的研究，特别是在加密过程中，在确保高程数值精度一定的情况下，如何顾及地表形态特征是获取高保真高程采样点数据的关键。顾及形态精度的高程采样点加密方法的缺失，使得在高程采样点数据不足区域，特别是在人工与自然结合、突变与渐变交错特征的地形区域进行地形建模带来很大困难，在诸如水土保持与防治、城市规划及雨洪管理等工程应用上甚至会造成相关部门的决策失误，从而给国家和人民造成巨大的经济损失。

技术实现要素：

本发明针对上述问题，面向地形高程采样点数据，设计了一种顾及形态精度约束的高程采样点加密的技术方案。该技术方案涵盖了具有自适应地表形态特征的线状高程采样点加密和区域高程采样点加密两大关键技术环节，能够有效解决高精度地形建模中高程采样点不足而导致地形表达失真的问题。

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

顾及形态精度的地形高程采样点加密方法，包括如下步骤：

第一步，加密应用信息初始化：获取加密应用信息及地形采样数据信息，并进行判定，若判定为线状要素高程采样点，则执行第二步；若判定为区域高程采样点，则执行第三步；

第二步，对线状要素高程采样点进行加密：首先，根据高程采样点及线状地形要素数据，将高程采样点映射到线状地形要素数据；其次，依据加密步距，计算出加密点的平面位置信息；最后，依据地表形态精度保真方法，对加密点高程进行插值计算；

第三步，对区域高程采样点进行加密：首先，根据高程采样点及区域边界要素，采用泰森多边形邻近法，计算出加密点的平面位置信息；其次，依据地表形态保真方法，对加密点高程进行插值计算；最后，依据加密迭代次数，重复上述过程，完成区域高程采样点的加密处理；

第四步，生成高程采样点加密结果文件：回收加密过程中申请的内存，清空临时的外存数据；更新日志文件，统计并记录加密元数据信息。

所述第二步的具体计算过程如下：

步骤21.高程采样点映射到线状地形要素：基于高程采样点与线状地形要素的最邻近空间位置关系，将高程采样点添加为线状地形要素的结点，并将非高程采样点的线状地形要素结点的高程值修改为-9999.0；

步骤22.基于映射后的线状地形要素和加密步距，通过局部坐标转换和形态保真计算，得到加密点的平面位置(x，y)，同时将加密点的高程赋值为-9999.0；其中，形态保真计算包括：①利用加密点P的前两个高程点和后一个高程点，或者前一个高程点和后两个高程点，计算点P与前后高程点在形态上是否为共线特征；②若具有共线特征，则利用线性解算获取点P信息；③若不具有共线特征，则通过利用点P的前后四个高程点拟合的三次曲线解算出点P信息，同时依据点P的形态特征距离(指当前点与其前一个点与后一个点连线的距离)与到其前一个高程点及后一个高程点形态特征距离的关系(点P的形态特征距离大于其前后高程点的形态特征距离最大值时)，进行点P信息调整；

步骤23.基于高程采样点及加密高程点平面位置，通过三维坐标降维和形态保真计算，对加密点进行高程插值计算，得到加密点高程值；其中三位坐标降维的具体过程为：根据点P及其前两个、后两个点，将三维坐标转换为二维坐标，建立以五个点中第一个和最后一个点连线为X轴且第一个点为起点的二维坐标系，即各点以起点的累计平面距离为坐标x′值、以及z值为坐标y′值。

本发明的技术特点及有益效果：

(1)本发明提出的地表形态精度保真方法，兼顾了线状地形要素及区域地形要素，在高程采样点数值精度不变的情况下，融入了地形语义信息，最大程度挖掘和保持地表形态特征，该方法为地形数采集规范完善以及DEM形态精度研究提供了技术基础。

(2)本发明构建并实现的地形高程采样点加密方法，在顾及地表形态特征的基础上实现了稀疏样点向密集样点的转变，使得采样数据有限时真实地形得以最优数字化表达。同时，该方法涵盖了线状和区域地形要素，实现了高程采样方式的基本覆盖。该方法提高了地形建模数据精度，特别是诸如城市地形等复杂区域地形建模的数据精度，也为面向对象地形建模提供了技术参考。

(3)本发明提出的线状要素高程采样点加密模块与区域高程采样点加密模块可以独立使用，也可以集成使用。二者具有较高的可移植性，为不同领域的采集数据加工提供了技术参考。

附图说明

图1是本发明的初始化工作流程图；

图2是本发明高程采样点映射到线状地形要素示意图；

图3是本发明局部坐标转换示意图；

图4是本发明线状地形要素高程采样点加密形态保真计算示意图；

图5是本发明线状地形要素高程采样点加密起始区间及终止区间处理方法；

图6是本发明高程加密点平面位置计算示意图(1,2,…,6为计算序号)；

图7是本发明三维坐标降维示意图；

图8是本发明泰森多边形示意图；

图9是本发明区域地形要素加密点平面位置计算流程图；

图10是本发明区域地形要素加密形态保真计算示意图；

图11是本发明生成结果文件的基本流程图；

图12是本发明高程采样点加密技术的流程图；

图13为实施例中(a)城市道路高程采样点和(b)城市道路高程采样点加密结果；

图14为实施例中(a)城市地块高程采样点和(a)城市地块高程采样点加密结果；

图15为实施例中(a)加密前、(b)加密后城市道路和城市地块DEM构建结果。

具体实施方式

下面结合附图和具体实施例，对本发明作进一步详细说明。本发明中设计和实现的过程中主要类设计如表1。

表1主要类以及功能

第一步，加密应用信息初始化：本发明设计了针对不同格式地形采样数据加密的应用信息初始化操作(见图1)。通过该初始化工作流程，获取加密类型、加密步距、加密迭代次数及地形采样数据的各类参数，并预测判定加密后数据容量是否符合数据输出要求。

加密应用信息初始化的具体实现过程为：

(1.a)进入GDAL环境，获取加密应用信息及地形采样数据信息，包括有加密类型(线状要素加密或区域要素加密)、加密步距、加密迭代次数、地形高程采样点数据、线状地形要素数据、区域地形特征要素数据(含区域边界数据)以及高程值的有效范围等。

(1.b)预判加密后高程采样点数据容量。根据加密步距、加密迭代次数、地形高程采样点数据以及线状地形要素数据，预测加密后高程采样点的数量。若采样点数量超出7000万个，则需做出参数调整或者结束当前加密操作。预测公式如下：

式中，V为加密后的采样点数量，N、M分别表示线状要素数量及每条线状要素线段数，L_ij为第i条线状要素第j条线段的长度，l为加密步距，v为加密前采样点数量，n加密迭代次数，a、b、c为系数项，本发明分别取0.25、0.24、1.36。

(1.c)调用LogManagementClass模块，将加密元数据信息(表2所示)写入日志文件。

表2高程采样点加密元数据

注：斜体部分为加密应用信息初始化写入日志信息。

第二步，执行线状要素高程采样点加密模块。首先，根据高程采样点及线状地形要素数据，将高程采样点映射到线状地形要素数据；其次，依据加密步距，计算出加密点的平面位置(x，y)信息；最后，依据地表形态精度保真方法，对加密点高程进行插值计算。

线状要素高程采样点加密的具体实现过程是：

(2.a)高程采样点映射到线状地形要素。基于高程采样点与线状地形要素的最邻近空间位置关系，将高程采样点添加为线状地形要素的结点，并将非高程采样点的线状地形要素结点的高程值修改为-9999.0，以便区分高程采样点和非高程采样点(见图2)。

(2.b)基于映射后的线状地形要素和加密步距，计算得出各加密点的平面位置(x，y)，同时将加密点的高程赋值为-9999.0。

(1)局部坐标转换。如图3所示，由P₁(x₁，y₁)、P₂(x₂，y₂)、P₃(x₃，y₃)、P₄(x₄，y₄)四个点建立以P₁P₄为x轴、P₁原点的X'O'Y'局部坐标系，计算公式为：

x′＝(x-x₁)*cosα+(y-y₁)*sinα

y′＝(y-y₁)*cosα-(x-x₁)*sinα

式中，α为P₁P₄的倾斜角，x、y为XOY坐标下的坐标值。同理，依据P(x′，y′)、P₁(x₁，y₁)及P₄(x₄，y₄)，可反解算出P在原始XOY下的坐标值(x，y)。

(2)形态保真计算。如图4所示，已知P₁(x₁，y₁)、P₂(x₂，y₂)、P₃(x₃，y₃)、P₄(x₄，y₄)四个点以及P点的坐标x值，求算P点的坐标y值的过程为：①计算点P₂到线段P₁P₃的距离D₁和线段P₁P₃的长度L₁，以及点P₃到线段P₂P₄的距离D₂和线段P₂P₄的长度L₂，若D₁/L₁或者D₂/L₂小于一个很小阈值(通常可取万分之一米)，则P₁、P₂、P₃或P₂、P₃、P₄三点共线，P点的坐标y按照直线方程计算，否则执行下一步；②由P₁、P₂、P₃、P₄拟合一条三次曲线y＝ax³+bx²+cx+d，由拟合公式计算出P点的坐标y，同时计算出点P到线段P₂P₃的距离D。若D大于max(D₁，D₂)，则执行下一步；③将D修改为max(D₁，D₂)，并延点P垂直线段P₂P₃方向调整点P坐标，进而得到最终点P坐标。

(3)加密点平面位置生成。①如图5所示，对于起始区间P₁P₂及终止区间P_n-1P_n，加密点按照定比分段线性加密，相关计算公式如下：

L＝(x_e-x_s)²+(y_e-y_s)²

λ＝i*l/(L-i*l)

x＝(x_s+λ*x_e)/(1+λ)

y＝(y_s+λ*y_e)/(1+λ)

式中，(x_s，y_s)、(x_e，y_e)分别为区间起始点和终止点坐标，L为线段长度，l为加密步距，i表示加密点序号(i＝1,2,3…L/l)；②对于任意非起始/终止区间P_kP_k+1，则取P_k-1、P_k、P_k+1、P_k+2四个线状地形要素结点(可以为高程采样点或非高程采样点)。首先，采用局部坐标变换方法将四个已知点变换为X'O'Y'局部坐标系，构成初始已知点集。其次，依据P_k′、P_k+1′的坐标x值和加密步距，线性计算得出系列加密点P′的坐标x值。再者，如图6所示，遍历已知点集取4个相邻点，从加密点集找出位于第二、三已知点中间(近似)位置点P′，采用形态保真计算法得到P′(x′，y′)，并将其加入到已知点集。重复上述步骤，直到所有加密点P′的坐标y值都计算完成。最后，采用局部坐标变换方法将系列加密点P′转换为原始XOY下的坐标值(x，y)，且高程值取-9999.0。

(2.c)基于高程采样点及加密高程点平面位置，对加密点进行高程插值。

(1)三维坐标降维。如图7所示，已知P₁(x₁，y₁，z₁)、P₂(x₂，y₂，z₂)、P(x，y，z)、P₃(x₃，y₃，z₃)、P₄(x₄，y₄，z₄)五个点，其中P₁、P₂、P₃、P₄为高程采样点，P为高程加密点，通过降维计算，将三维坐标转换为二维坐标，即各点以P₁为起点的累计平面距离(XOY平面)为坐标x′值、z值为坐标y′值，计算公式为：

x₁′＝0

y₁′＝z₁ y′＝z y′₄＝z₄

y₂′＝z₂ y₃′＝z₃

式中，L_pipj为两点间的XOY平面距离。当二维恢复为三维时，由各点(x，y)与其对应的y′组合为(x，y，y′)即可。

(2)形态保真计算。计算过程同步骤2.b中形态保真计算。当将D修改为max(D₁，D₂)，并延点P垂直线段P₂P₃方向调整点P坐标时，这里只调整P的坐标y值。

(3)加密点高程插值。①当线状地形要素只有1个高程采样点P(x，y，z)时，所有加密点高程值修改为z；②当线状地形要素高程采样点大于1个且小于4个时，则依据距离权重对各加密点进行线性高程插值，相关计算公式如下：

当加密点P位于采样点P_s与P_e之外且位于P_s一侧

当加密点P位于采样点P_s与P_e的之间

式中，L_pipj为两点间的XOY平面距离；③当线状地形要素高程采样点大于4个时，对于高程采样点起始区间P₁P₂及高程采样点终止区间P_n-1P_n采用上述②方法处理。对于任意高程采样点非起始/终止区间P_kP_k+1，则取P_k-1、P_k、P_k+1、P_k+2四个高程采样点作为初始已知点集，P_k、P_k+1之间的加密点高程插值过程为：首先，查找出位于P_k、P_k+1之间的加密点集。最后，遍历已知点集取4个相邻点，从加密点集找出位于第二、三已知点中间(近似)位置点P，先采用三维坐标降维方法处理，再通过形态保真计算得到点P′的坐标y′值，并将其恢复为三维坐标P(x，y，y′)后加入到已知高程采样点集。重复上述步骤，直到所有加密点P的高程值都计算完成。

第三步，执行区域高程采样点加密模块。首先，根据高程采样点及区域边界要素，采用泰森多边形邻近法，计算出加密点的平面位置(x，y)信息；其次，依据地表形态保真方法，对加密点高程进行插值计算；最后，依据加密迭代次数，重复上述过程，完成区域高程采样点加密处理。

区域高程采样点加密的具体实现过程是：

(3.a)以高程采样点为初始高程点集，基于区域边界数据，计算得出各加密点的平面位置(x，y)，同时将加密点的高程赋值为-9999.0。

(1)泰森多边形生成。①由高程点集构建Delaunay三角网，并对高程点及三角形进行编号，从而记录高程点与三角形的对应关系；②对于每个高程点，按顺时针或逆时针为序找出以其为顶点的所有三角形，并计算这些三角形的外接圆圆心；③连接这些与高程点邻接三角形的外接圆圆心，即可生成该高程点的泰森多边形，如图8所示。

(2)加密点平面位置计算。如图9所示，①基于高程采样点和区域范围，采用泰森多边形生成法，获取区域高程采样点的泰森多边形；②根据泰森多边形及区域边界，提取非区域边界的泰森多边形边；③提取非区域边界的泰森多边形边中点，生成加密点平面位置(x，y)，并将其高程值赋值为-9999.0。

(3.b)加密点高程插值计算。基于初始高程点集及加密高程点平面位置，对加密点进行高程插值。

(1)形态保真计算。如图10所示，计算点P高程的过程为：①构建高程点Delaunay三角网。利用高程点集和区域地形线状特征要素(如陡坎线、地块边界线等)构建Delaunay三角网；②查找点P高程插值所需高程点。按照Delaunay三角网构建法则，利用点P编辑更新高程点Delaunay三角网，查找出与点P相连的高程点P₁、P₂、P₃。在查找过程中，地形单元边界线起到了阻碍作用，致使位于地形单元边界线P点相对一边的高程点无法被搜索到，从而保证了地形单元形态完整性和保真性；③插值计算点P高程。利用反距离权插值法，计算得到点P高程值，相关公式为：

式中，n为插值所需高程点个数，d_i为第i个高程点与点P的平面距离，z_i为第i个高程点的高程值，λ_i为第i个高程点的权重值。

(2)加密点高程插值。逐个遍历加密点集合，采用形态保真计算法计算得出加密点高程。若在形态保真计算过程中，无法查找到插值所需高程点，则剔除该加密点。

(3.c)迭代完成加密计算。将初始高程点集与加密点集合并为新的初始高程点集，重复上述步骤，直到达到指定加密迭代次数为止。

第四步，生成高程采样点加密结果文件。对加密计算过程中申请各种内存资源进行回收，清空临时的外存数据；更新日志文件，统计并记录加密元数据信息。生成结果文件的基本流程图如图11所示。

以上第一～第四步操作实现了线状地形要素和区域高程采样点加密技术，完成了顾及形态精度的地形高程采样点加密方法(如图12所示)。

实施例

下面以城市地形(道路和城市地块)高程采样点加密为例。

城市DEM是实现城市三维可视化表达、分析、规划的基础支撑。城市DEM通过有限的采样数据数字化表达地表形态，其构建核心是城市地表几何信息及语义信息即形态信息的获取，即城市地形采样数据的获取。在这些采样数据中，高程采样点是最为主要且最为直接的地形数据，对城市地形建模起着关键核心作用。城市地形被道路骨架网分割，并由城市道路和城市地块两大模块组成，其中城市道路被定义为次干路及以上等级的城市道路(依据《城市道路工程设计规范》)，城市地块是由城市道路围合而成的斑块区域。因此，城市地形建模可以分解为城市道路和城市地块地形建模两部分。然而，无论城市道路区域还是城市地块区域，从1:500大比例尺地形图中获取的城市地形高程点较为稀疏，需要对其进行地形高程点加密，以有效满足高精度城市地形建模的需求。

顾及形态精度的城市地形高程采样点加密过程为：

第一步、加密应用信息初始化。

(1)输入线状地形要素高程采样点加密参数：城市道路高程采样点、城市道路边界线、步距1m；

(2)输入区域高程采样点加密参数：城市地块高程采样点、城市地块边界线、陡坎线、等高线、地形单元边界线、迭代次数3；

(3)输入加密类型、高程值的有效范围、工作目录等参数；

(4)调用JudgeInfoClass模块，进行加密后数据容量进行预判；

(5)调用LogManagementClass模块，将相应的加密元数据写入日志文件。

第二步，城市道路高程采样点加密。调用LineDensifyClass模块，依据加密步距，对城市道路高程采样点进行加密处理。

第三步，城市地块高程采样点加密。调用AreaDensifyClass模块，依据加密迭代次数，对城市地块高程采样点进行加密处理。

第四步，生成城市道路和城市地块高程采样点加密结果。清理回收加密过程中申请的内存以及外存数据，并更新日志文件，记录相应加密元数据信息。

城市道路和城市地块高程采样点加密结果如图13、14所示。分别利用本发明的加密前高程采样点和加密后的高程采样点数据，以及相同的地形特征要素数据，构建了加密前后的城市道路DEM和城市地块DEM，如图15所示。从图15可看出，本发明有效解决了因高程采样点稀疏而引起地形建模形态失真问题。