一种机载雷达点云地形产品自动化生产方法及装置

技术领域

本发明属于图像信息处理领域，尤其涉及一种机载雷达点云地形产品自动化生产方法及装置。

背景技术

机载激光雷达(Light Detection and Ranging-LiDAR)是一种主动对地观测的新型技术，它克服了传统逐点测量的方式，不受天气和光照的影响，能够短时间快速、低成本的获取大范围地表场景的三维地理信息，在此基础上生产的地形产品在水文、气象、森林普查、电力巡检、地质、军事和国防建设中都有着广泛的应用。地形产品类型主要包括：数字表面模型(Digital Surface Model，DSM)，数字高程模型(Digital Elevation Model，DEM)，等高线等产品。

由于原始点云数据不包含分类信息，建筑物、植被、桥梁等地物点云，不能真实反应场景内的地形，所以从LiDAR点云数据中自动化识别地面点和非地面点是生产地形产品的前提，即点云滤波处理。虽然目前已经提出了很多地面点滤波算法，且很多算法在实际生产应用中也取得了比较理想的效果，但还在存一定的缺陷。首先，现有算法不能自适应处理复杂场景的数据，需要人为根据场景的地形特征设置不同的滤波参数，例布料模拟滤波算法(Cloth Simulation Filter，CSF)(已集成于CloudCompare)需要人为指定场景类型参数才能完成滤波，对于包含平原山地等地形的场景滤波效果不佳；渐进形态学滤波算法通过形态学腐蚀膨胀操作完成滤波，腐蚀膨胀的窗口大小需要根据地形和点密度手动设置；渐进三角网滤波算法(已集成于Terrasolid)通过对场景绘制规则格网，取最低点构建初始地形，逐步加密完成场景滤波，但初始格网大小的选择、沟壑、山脊等地形的处理效果欠佳。

数字表面模型DSM和数字高程模型DEM是基于初始点云和地面点生产的地形产品，通常对点云数据进行规则格网化来获取，对于地面点云数据而言，由于非地面点的缺失造成了地面点云中存在不同大小的空洞，常见的插值方法通常需要人为指定插值范围大小，需要一定的人工参与。

等高线地形产品的生产通常通过立体采集的方法，人为等间隔获取一些高程点通过插值构建等高线；常见的自动化提取算法也存在一定的不足，如等高线自相交、等高线锯齿状等问题，也需要进行后处理，一定程度降低了生产的效率。

发明内容

本发明旨在解决上述问题，提供一种机载雷达点云地形产品自动化生产方法及装置。

第一方面，本发明提供一种机载雷达点云地形产品自动化生产方法，包括DSM/DEM栅格地形产品生产、基于DEM栅格等高线生产和基于地面点等高线生产；

所述DSM/DEM栅格地形产品生产包括如下步骤：

S11、首先对场景全部/地面的点云数据进行处理；

S12、对处理后的点云数据构建三角网，依据需求判断是否进行体素化下采样，逐像素进行三角网定位；

S13、通过三角网顶点内插该像素位置的高程值对栅格进行赋值处理；通过建立三角网的方式进行插值生产无空洞的DSM/DEM数据；

所述基于DEM栅格等高线生产包括如下步骤：

S21、首先设定预设参数，依据输入的参数判断每条等高线的高程值，以指定的高程值对栅格数据进行二分类掩膜数据生产，高于该值则标记255，否则标记0；S22、其次采用线性内插的方式在栅格四边计算等高点并连线，判断等高线的类型；为了避免等高线穿过栅格角点，对所有栅格高程值增加0.00001米；

S23、然后采用等高点追踪方法，以任意一段等高线为起点，双向追踪处理，获取单条等高线；

S24、最后对等高线进行平滑处理；克服锯齿状的现象；

所述基于地面点等高线生产包括如下步骤：

S31、首先设定预设参数；

S32、其次建立测区地面点的三角网，判断三角网内是否存在等高线；如有，则采用线性内插的方式，在三角网三条边上确定等高点；

S33、然后采用等高点追踪方法，以任意一段等高线为起点，双向追踪处理，获取单条等高线；

S34、最后对等高线进行平滑处理，克服锯齿状的现象。

进一步，本发明所述机载雷达点云地形产品自动化生产方法，所述步骤S11具体包括如下步骤：

S111、噪点剔除；采用统计滤波的方式进行噪点剔除，建立点云数据Kdtree，采用统计K邻域内点云距离标准差的方式；初始点云中不免会存在噪点数据，通常认为最低的点为地面点，但实际数据中存在不少远低于地面的噪点数据，这些点通常表现为离散、高差异常；

S112、场景粗分类；选择场景内最大建筑物边长或场景内最大非地面点连通地物长度确定格网大小，选择格网内最低点为格网值，求解xy方向梯度，将梯度影像拉伸至0-255，并采用最大类间差法(OTSU)对场景进行前景和背景值分类，生产测区地形起伏的二值图，其中255为前景表示地形起伏大的区域，0为背景表示地形起伏较小的区域；

S113、渐进三角网滤波；采用多级格网划分处理获取初始地面点数据，并建立三角网，恢复初始地形，以迭代角、迭代距离、镜像判断为条件依次对新加入的点进行判断；如果是地面点，则更新三角网，否则跳过；添加新点的顺序是以初始地面点为中心，距离较近的优先判断；由于通常情况下，离初始地面点越近，越可能是地面点，添加新点的顺序会直接影响点云滤波的效果；

S114、地形二次恢复；在第一次滤波结果的基础上，以地面点为基础，以预设的半径范围，逐次向外扩散进行第二次地面点的判断，判断条件为点距离和生长角，实现地形的二次恢复；由于初始格网大小的不确定性以及三角网判断条件的缺陷，虽然引入镜像处理，但沟壑、陡坎、山脊等位置的地面点往往不满足三角网的判断条件，因此进行地形的二次恢复，该步骤可以恢复山脊被“削平”，陡坎边界缺点的现象；

S115、错分点剔除；以现有地面点为基础，以2-3倍平均点间隔为格网分辨率构建DEM，格网值选择格网内点云最低点，以预设高度(0.5米)构建缓冲层，将高于该缓冲层的错分数据进行剔除处理，得到地面点和非地面点，即处理后的点云数据。

进一步，本发明所述机载雷达点云地形产品自动化生产方法，步骤S113所述采用多级格网划分处理过程中，对山地区域进行二级格网细化，格网大小为初始格网的一半，对地形起伏大的区域(即步骤S112中标记为255的区域)的初始地面点进行补点操作，同时进行误差剔除处理。

进一步，本发明所述机载雷达点云地形产品自动化生产方法，所述步骤S113，步骤S114采用并行处理方式进行；设定单个处理单元的大小，对场景进行规则格网分块；通过采用并行处理方式可有效的提高大场景点云遍历效率，提高硬件内存的兼容范围。

进一步，本发明所述机载雷达点云地形产品自动化生产方法，所述格网分块在边界处构三角网的分块点云添加虚拟格网点。

进一步，本发明所述机载雷达点云地形产品自动化生产方法，所述预设参数包括等高线初始高程、终止高程、基本等高距、是否生产计曲线和助曲线成果类型。

第二方面，本发明提供一种机载雷达点云地形产品自动化生产装置，包括相电连接的处理器和存储器；所述存储器用于存储计算机程序；所述处理器执行前述计算机程序时，可实现前述第一方面所述的机载雷达点云地形产品自动化生产方法。

第三方面，本发明提供一种计算机可读存储介质，所述存储介质上存储有计算机程序；所述计算机程序被执行时，可实现前述第一方面所述的机载雷达点云地形产品自动化生产方法。

本发明所述机载雷达点云地形产品自动化生产方法及装置，采用三角网插值的策略，省略了插值范围参数，可以生产无空洞缺值的DSM/DEM地形产品；同时支持栅格和点云数据，克服多源数据的输入，并行处理每条等高线，并采用三次贝塞尔曲线对等高线进行平滑拟合，高效自动生产合理光滑的等高线产品；在点云数据处理时以渐进三角网滤波算法为基础，对地形粗分类，建立多级细化格网，适应复杂地形混合存在的场景；在初次虑结果基础上，采用点间隔及地形生长角为阈值对地形进行二次恢复，弥补山脊、陡坎地形缺失点的情况；算法对场景规则分块，采用并行处理策略，极大的提高了滤波的效率。

附图说明

图1为本发明实施例所述点云数据进行处理流程示意图；

图2为本发明实施例所述机载雷达点云地形产品自动化生产流程示意图；

图3为本发明实施例所述三维Kdtree示意图；

图4为本发明实施例所述点云分块和初始地面点多级格网细化示意图；

图5为本发明实施例所述多级格网细化加密点及虚拟点补充示意图；

图6为本发明实施例所述陡坎点镜像处理示意图；

图7为本发明实施例所述地形生长角定义示意图；

图8为本发明实施例所述二次地形恢复迭代外扩过程示意图；

图9为本发明实施例所述滤波处理结果进行对比示意图；

图10为本发明实施例所述掩膜图像中等高线可能存在情况示意图；

图11为本发明实施例所述三次贝塞尔曲线拟合示意图；

图12为本发明实施例所述等高线三次贝塞尔曲线平滑示意图；

图13为本发明实施例所述不规则三角网中等高线存在情况示意图；

图14为本发明实施例所述三角网等高线线段双向追踪示意图；

图15为本发明实施例所述地形产品成果展示示意图。

具体实施方式

下面通过附图及实施例对本发明所述机载雷达点云地形产品自动化生产方法及装置进行详细说明。

实施例一

本公开实施例公开一种机载雷达点云地形产品自动化生产方法中对场景全部/地面的点云数据进行处理的方法，如图1所示，具体包括如下过程：

(1)建立点云Kdtree：

点云数据区别于栅格数据，它是一种离散数据，无法像栅格一样通过下标进行邻域检索，故在处理前需要对场景点云建立Kdtree，加速邻域搜索速度。Kdtree是一种划分k维数据空间的数据结构，其本质也是一种二叉树，建立Kdtree是一个不断划分的过程，选择合适的轴当做分割面，随着深度的不断增加，直到建立所有点云的Kdtree，点云的检索则是通过根节点不断向下递归检索判断实现的。

(2)点云噪点剔除：

在本公开实施例中，点云噪点剔除采用点云及其邻域点的距离标准差阈值进行剔除，是一种全局滤波算法。如图3所示，建立点云Kdtree，通过Kdtree选取当前点的K邻域，计算当前点和所有邻域点的距离，统计全部点云的距离均值和标准差，采用以下公式计算距离阈值，大于该阈值即为噪点，需要剔除处理。

Dis

式中，Dis

(3)场景粗分类：

设定场景内最大建筑物或场景内最大连通地物的长度，对点云场景进行规则格网划分，格网高程值选择落入该格网的高程最低值；考虑到点云可能有空洞区域，故还需要对空值区域进行插值处理；分别在x,y方向进行梯度计算合并，获得xy方向的综合梯度图；进行灰度拉伸，将灰度值拉伸至0-255；采用最大类间差法(OTSU)对灰度图进行二分类，前景区域即为地形起伏明显区域，背景区域为起伏平缓区域，OTSU处理梯度图而不是深度图，其目的是为了识别地物变化明显(梯度大)的区域而非高程高的区域。

其中xy方向的梯度计算公式为：

式中，gradf(x,y)为xy方向梯度，

其中OTSU阈值自动确定方法可以通过下面公式实现。

D(T)＝P

图4为地形粗分类示意图，图4a为测区场景，主要包含平原、丘陵级沟壑几种地形，图4b为xy方向融合梯度图，可以发现高程变换平缓的区域灰度值较低，灰度值高的区域高程变换明显，图4c表示最大类间差法OTSU自动阈值分类的前景和背景区，坡度变化剧烈的区域都被标记。

(4)点云分块和初始地面点多级格网细化：

在获得OTSU二分类结果基础上，对前景(即灰度为255)进行格网细化加密，细化格网大小为初始格网的一半，格网值依旧为格网内点云最低值，新加入的点同样需要粗差点剔除工作；对点云进行规则格网分块处理，每块点云即为并行处理单元；考虑到测区和分块点云边界建网的问题，需要进行虚拟点云补充处理，补充点的高程选择离其最近点的高程。图5为多级格网细化加密点和前景区域补充的虚拟点示意图，相比于规则格网点，在地形起伏变化剧烈的区域补充加密点可以提高地形描述的准确性，由于该场景过小，暂不展示分块效果。

(5)渐进三角网滤波：

通过初始地面点建立不规则三角网，即TIN，计算场景内点云与初始地面点的距离并排序，从最小距离开始逐点进行判断，是否能够满足加密条件，加密条件为迭代角和迭代距离。首先寻找当前点在哪个三角网投影XOY平面的三角形中，依次计算当前点与三角网三个顶点的夹角和点到三角网平面的距离，并与迭代角和迭代距离阈值进行对比，小于该阈值则为地面点，并以该点更新地面三角网。考虑到建网和三角网遍历的效率问题，可以设定最小三角网的边长来控制三角网的更新。

渐进三角网滤波算法对陡坎处理加入了镜像处理的方式，如果三角网与XOY平面的夹角达到预设角度，则认为它为陡坎点，即下图的P

(6)二次地形恢复

在实际生产过程中，发现通过上述方法生产的结果中陡坎点、沟壑及山脊地形很多地面点都无法恢复，究其原因为这些点无法满足三角网判断的条件，故在第一次滤波的基础上对已有地形进行二次恢复。地形生长角的定义为，当前点与最近地面点R邻域范围拟合平面的夹角即图7表示的θ。具体方法是基于现有地面点以预设的半径范围进行迭代扩张，对扩张的点需要计算地形生长角，如果小于该阈值，则纳为地面点，如图8所示黄色区域为渐进三角网初次滤波的结果，蓝、红、绿和粉色为二次地形恢复过程中逐次迭代外扩的结果，陡坎边界区域缺失的点云已被恢复。

(7)错分点剔除

本实施例所述方法在对场景点云进行处理后，不可避免的会产生一些错分点，尤其是第二次地形恢复处理过程中，有个别非地面点也满足地形生长条件，故需要对错分点进行剔除。基于现有地面点进行DEM构建，格网值选择格网内点云最低值，在DEM基础上构建预设高度(0.5米)的缓冲层，高于该缓冲层的点云需要被归为非地面点。

(8)地面点滤波结束

通过上述描述的流程可以自适应、高效的完成复杂场景的地面点和非地面点的滤波处理，该数据为后续地形产品的生产提供了数据支持。

针对该测区数据，在本公开实施例中，选择PCM、CloudCompare两款软件中内置的滤波算法与本发明滤波处理结果进行对比，如图9所示，其中9a为三角网滤波结果示意图；图9b为形态学滤波结果示意图；图9c为布料模拟滤波结果示意图；图9d为本实施例所述方法处理结果示意图；红色表示非地面点，蓝色表示地面点，可以很直观的发现，本实施例所述方法滤波结果最好，在复杂地形情况下可以达到相对满意的处理结果，得到处理后的点云数据。

实施例二

在上述实施例一的基础上，本公开实施例公开一种机载雷达点云地形产品自动化生产方法，如图2所示具体包括如下过程：

(1)DSM和DEM生产

首先确定输出影像的分辨率，如不指定，则依据平均点间隔为分辨率输出，点云平均但间隔通过计算每个点与其最邻近点距离的均值确定；

依据场景点云数据，计算外包围盒范围，依此确定栅格影像的仿射变换矩阵，栅格产品以外包围盒左上角平面位置确定，影像宽高通过测区长宽和影像分辨率计算，初始影像栅格全部置-9999无效值；

依据分辨率大小确定是否需要对点云数据进行体素化下采样，即在xyz三个方向绘制规则格网，下采样间隔为影像分辨率大小，选择体素内距离体素中心最近的点为体素采样点，通过这种方法可以缩减构建三角网的数量，在不影响栅格精度基础上提高算法运行效率；

对下采样后的点云数据构建不规则三角网TIN，逐行列遍历像素，将像素坐标转换为地理坐标，确定该像素点存在于哪个三角网中，依据三角网三个顶点的空间坐标，拟合平面方程，通过平面坐标求解Z值坐标赋给栅格图像。其中，像素坐标转地理坐标、点是否在三角形(任意多边形)内部及空间平面方程的计算方式如下：

XGeo＝GeoTransform[0]+GeoTransform[1]*X

YGeo＝GeoTransform[3]+GeoTransform[4]*X

式中，地理坐标XGeo,YGEO分别采用X

t＝(P.x-B.x)*(A.y-B.y)-(A.x-B.x)*(P.y-B.y)

式中，A,B分别表示多边形前一个点和后一个点，.x,.y表示平面坐标，P表示当前点，迭代所有点，如果t<0个数为奇数则在任意多边形内部，反之不在。

Ax+By+Cz+D＝0(C≠0)

令：

采用最小二乘求解x:x＝(A

式中，(x

(2)等高线生产

本实施例为了兼容不同源数据数据的进行等高线生产，提供了基于DEM栅格的等高线生产策略和基于地面点云数据的等高线生产策略，能够满足数据的多元化输入，提高了算法的适应性。

基于DEM栅格的等高线生产过程包括：

首先确认等高线起始和终止高程，基本等高距和是否生产间曲线和助曲线等参数，计算场景内需要生产等高线的条数及对应高程，以每个高程的等高线为最小单元进行并行处理；

处理每条等高线策略，创建和栅格图相同大小的掩膜影像，初始化为0，遍历全部DEM栅格，如果栅格像素值大于等高线高程值，则该像素赋值255，生产对应等高线的掩膜影像；

图10表示所有掩膜影像中可能存在的所有情况，黑色值对应掩膜图像值为0，白色值对应掩膜图像为255，红色线表示等高线。逐像素进行判断，如果该像素属于下图中第一个或最后一个则不处理该像素，如果为其他剩余14种栅格，则在对应栅格种生产对应的等高线段，等高线段两端点必定与栅格四边交于两/四点(一/两条等高线)，端点的位置采用线性内插的方式，考虑到端点值有可能会处于栅格角点处的情况，故对每个栅格值进行增加0.00001米的操作。

在本公开实施例中等高线段追踪包括，将掩膜全部像素全部置为未标记状态，栅格内存在两段等高线的栅格特殊标记；以任意一段等高线为起点，记录其两个端点位置，从任意一端点开始寻找下一段等高线，将遍历过的栅格置于标记状态，直到找到所有等高线段，从另一端点开始反向寻找，以相同方法进行遍历，直到找到该条等高线的所有线段，并在端点处反向顺次连接；如果等高线两端端点距离相近，则对等高线进行闭合处理；反复操作，可以获得当前高程值的所有等高线；如果等高线过短，可以采用距离阈值对等高线进行删除，避免因为地形异常值造成多层等高线局部闭合现象的产生。

在本公开实施例中等高线平滑处理包括，基于栅格的生产的等高线其实是由多个等高线段顺次连接而成，等高线成果直观上会出现“锯齿状”的现象，为了去除这种现象，本发明采用了三次贝塞尔曲线对等高线进行平滑处理，生产的等高线更贴合实际的测区场景，值得一提的是，三次贝塞尔曲线是以四个等高点为控制点，以预设的步长(如0.1)重构并新增等高点并依次连接而成，一定程度上细化了等高点数据，对于分辨率较大的栅格影像也可以拟合出较为精细的等高线数据；

贝塞尔曲线实际上是以n个控制点来控制曲线的形状，通常用来绘制矢量曲线，贝塞尔曲线经过首尾两点，不经过中间点，设定合适步长，以图11所示四个等高点P

其中，三次贝塞尔曲线的计算方式如下：

式中，B(x,t),B(y,t),B(z,t)分别表示通过四个等高控制点P

图12表示三次贝塞尔曲线在平滑等高线中的实际应用，黑色维平滑前的等高线，红色表示平滑后的处理结果，可以很直观的发现，平滑后的等高线克服了“锯齿状”的现象，更美观，更符合实际地形的变化。

本公开实施例中对上述的等高线进行矢量写出，矢量类型为折线型，属性字段信息分别为Shape(矢量类型)，ID(等高线ID)，elev(高程字段)及Type(等高线类型：首曲线、助曲线、计曲线及间曲线)。

基于地面点云的等高线生产过程包括：

首先确认等高线起始和终止高程，基本等高距和是否生产间曲线和助曲线等参数，计算场景内需要生产等高线的条数及对应高程，以每个高程的等高线为最小单元进行并行处理；

对地面点云数据建立不规则三角网，区别于栅格影像，在三角网中，等高线只存在图13所示的三种情况，端点的位置采用线性内插的方式，考虑到端点值有可能会处于三角网顶点处的情况，故对每个三角网顶点高程增加0.00001米的操作，可以有效避免；分别在三角网中生产对应高程的等高线段；

如何判断三角网中是否存在等高线段，可以通过以下公式计算和判断：

如果n

式中，P

等高线段追踪和等高线平滑在前述基于DEM栅格等高线生成叙述中已详细描述，在此不再赘述，方法同上。图14表示不规则三角网中等高线段双向追踪示意图，红色为起始等高线段，黄色为左端点追踪的等高线，蓝色为右端点追踪的等高线段，故在构建三角网的过程中需要保存三角网和邻接三角网的关系；等高线矢量写出，写出内容和方法同上。

依据本实施例上述流程，对该测区进行多种地形产品高效自动化生产，部分地形产品成果如图15所示，其中图15a为测区DSM地形产品示意图，图15b为测区DEM地形产品示意图，图15c为测区等高线地形产品示意图，图15d为山体、断崖处等高线局部示意图；DSM/DEM分辨率为点云平均点间隔3米，等高线产品基本等高距10米，由于基于栅格和基于点云生产的等高线几近重合，故只展示一种产品。

从算法运行效率看，测区点云数量约96000+，硬件条件8核CPU、16G运存，点云滤波耗时2s，生产DSM/DEM耗时1s，生产等高线1s，算法运行效率比较高；从产品质量视角看，DSM/DEM产品能够真实的反映场景内的地物和地貌；等高线产品紧密贴合于点云表面，能准确反映测区的地形起伏变化，同时在平原、山地及断崖处都有较好的表现。

实施例三

本公开实施例公开一种机载雷达点云地形产品自动化生产装置装置，包括相电连接的处理器和存储器；所述存储器用于存储计算机程序；所述处理器执行前述计算机程序时，可实现前述实施例一、实施例二所述的机载雷达点云地形产品自动化生产方法，具体生产方法与前述实施例一、实施例二相同，再此不再赘述。

实施例四

本公开实施例公开一种计算机可读存储介质，所述存储介质上存储有计算机程序；所述计算机程序被执行时，可实现前述实施例一、实施例二所述的机载雷达点云地形产品自动化生产方法，具体生产方法与前述实施例一、实施例二相同，再此不再赘述。

本申请实施例所述计算机可以是通用计算机、专用计算机、计算机网络、或者其他可编程装置。所述计算机指令可以存储在计算机可读存储介质中，或者从一个计算机可读存储介质向另一个计算机可读存储介质传输。所述计算机可读存储介质可以是计算机能够读取的任何可用介质或者是包含一个或多个可用介质集成的服务器、数据中心等数据存储设备。所述可用介质可以是磁性介质，(例如，软盘、硬盘、磁带)、光介质(例如，数字通用光盘(DVD))或者半导体介质(例如，固态硬盘(SSD))等。计算机存储代码所形成软件可以位于随机存储器，闪存、只读存储器，可编程只读存储器或者电可擦写可编程存储器、寄存器等本领域成熟的存储介质中。

在本申请各实施例中的各功能模块可以集成在一个处理单元或模块中，也可以是各个模块单独物理存在，也可以两个或两个以上模块集成在一个单元或模块中。在上述实施例中，可以全部或部分地通过软件、硬件、固件或者其任意组合来实现。当使用软件实现时，可以全部或部分地以计算机程序产品的形式实现。所述计算机程序产品包括一个或多个计算机指令。在计算机上加载和执行所述计算机程序指令时，全部或部分地产生实现本申请实施例所述的流程或功能。

本发明所述机载雷达点云地形产品自动化生产方法及装置，能够自适应、高效、自动的完成复杂场景下机载LiDAR点云地面点滤波、DSM/DEM和等高线地形产品的生产。算法具有较强的数据兼容性、复杂的地形适配性及较少的参数自适应性，经过大量数据实验和工程化任务生产证明，该方法和装置有着较好的表现。