技术摘要：
本发明公开了一种适用于机载LiDAR点云的地形自适应插值滤波方法。该方法包括低异常点去除、地面种子点选取、地面种子点核查以及地面点渐进选择等步骤。其中，在低异常点去除步骤中，本发明采用高程直方图实现低异常点自动剔除，以降低人工干预程度；在地面种子点选取步  全部
背景技术：
目前，机载激光雷达(LiDAR)技术以其主动性和高效性等优势已成为三维空间数 据采集的主要手段，广泛应用于数字高程模型生产、三维城市构建、森林资源调查以及山体 滑坡监测等领域。区别于直接接触式数据采集方式，原始机载LiDAR点云数据包括各种地面 点和非地面点，因此，原始机载LiDAR点云数据在使用前必须进行滤波处理。 过去二十年间，国内外众多研究人员提出了各种点云滤波算法。根据工作原理，现 有的滤波算法可分为坡度滤波、形态学滤波、插值滤波、分块滤波以及机器学习滤波等。 坡度滤波假设相邻地面点坡度较小，因此，如果两点间坡度超过设定阈值，则较高 的点为非地面点。然而，该滤波方法对坡度阈值非常敏感，特别在地形复杂区域容易将地形 特征点错分为非地面点。为此，相关人员提出了一系列自适应坡度滤波方法，确保坡度阈值 随地形复杂度自适应变化。由于坡度滤波仅利用相邻点之间的局部信息，该方法在地形断 裂区和陡坡区滤波误差较大。 形态学滤波主要借助形态学开运算对点云滤波，即以一定尺寸的结构元开运算栅 格化点云，如果处理前后栅格点高程差超过规定阈值，则对应点标记为非地面点。其中，结 构元尺寸显著影响滤波精度：该尺寸必须足够大以能够去除各种尺寸非地面物，同时，该尺 寸必须足够小以尽可能保留地形细节信息。为此，相关人员提出了一系列渐进形态学滤波 算法以剔除各种尺寸非地面物。然而，形态学滤波以栅格化点云为操作对象，因此在点云栅 格化过程中不可避免的造成精度损失。 插值滤波首先对地面种子点插值构建初始地面参考面，接着比较待滤波点与地面 参考面的距离，如果该距离小于设定阈值，则对应的点标记为地面点，最后用更新后的地面 种子点重复上述过程直至没有点标记为地面点。 作为插值滤波的典型代表，渐进三角网滤波因被纳入商业软件TerraSolid而被广 泛使用。然而，该滤波方法使用TIN插值构建地面参考面，因此在地形复杂区域滤波精度较 低。 为此，部分学者建议使用较高精度的非线性插值方法(如线性估计、薄板样条和克 里金等)代替TIN插值，并提出了多分辨分层插值滤波方法。该类滤波方法在计算过程中需 要构建一系列不同分辨率DEM作为地面参考面，并比较每个待分类点到对应DEM网格点的距 离。然而，当待分类点与DEM网格中心点不重合时，势必产生空间位置误差，而且该空间位置 误差与地形坡度呈正比，即坡度越大，误差越大。因此，插值滤波在陡坡区域滤波精度较低。 分块滤波以对点云分割后的块为基本处理单元，并利用块和块之间的拓扑关系、 块和地面参考面的高程差等对各个块滤波。鉴于点云块可捕捉到地形断裂点，理论上分块 滤波在地形断裂区具有较好的滤波效果。然而，分块滤波的优势是以对点云准确分块为前 提，不准确分块容易抵消分块滤波的高效性。而且，该类滤波在地形复杂区域由于各个块中 4 CN 111598780 A 说　明　书 2/12 页 点数较少，优势并不明显。 近年来，机器学习方法被广泛使用点云滤波。常用的机器学习算法包括人工神经 网络、条件随机场、支持向量机和深度学习等。然而，机器学习滤波方法需要大量的训练样 本，制作过程费时费力；而且，该类滤波具有计算量庞大以及人工参与程度高等缺陷。 根据上述分析，表1总结了上述五种滤波方法的优缺点。 表1各种机载LiDAR点云滤波方法优缺点比较 大量实验分析表明，现有的各种滤波方法在场景简单区可以取得较好的滤波精 度，但在场景复杂区域容易引起点云漏分或错分问题，降低了滤波精度，后续需要大量人工 编辑工作，费时费力，而且点云编辑精度完全取决于操作人员的熟练程度，主观性强。
技术实现要素：
本发明的目的在于提出一种适用于机载LiDAR点云的地形自适应插值滤波方法， 以提升各种场景下点云滤波的精度。本发明为了实现上述目的，采用如下技术方案： 一种适用于机载LiDAR点云的地形自适应插值滤波方法，包括如下步骤： I.利用高程直方图自动去除原始点云中的低异常点； II.选取地面种子点，具体过程如下： 首先将原始点云的研究区域覆盖一定尺寸的网格，并选择每个网格中的最低点作 为候选点；记录每个候选点的三维坐标，而且候选点按照栅格方式存储； 然后以栅格为操作对象，分别逐行和逐列选择初始地面种子点； 最后将行方式和列方式都选中的点标记为地面种子点； 其中，从栅格的任意一行或列栅格数据中选择初始地面种子点的步骤为： II.1.初始地面种子点选择； 利用一维移动窗口选择初始地面种子点；其中，选中的网格点含有高程值，未选中 的网格点标记为NaN；一维移动窗口的尺寸取值为最大建筑物的短边长，设为S； II.2.样条插值； 借助一维离散光滑样条法对初始地面种子点插值，估计其它未选中网格点的值； II.3.网格点高程差计算； 首先利用原始网格点高程减去插值后网格点高程，得到网格点高程差；然后将网 5 CN 111598780 A 说　明　书 3/12 页 格点高程差小于设定阈值eth的点标记为地面点； II.4.地面种子点更新； 利用步骤II.3中选中的地面点更新地面种子点； II.5.重复上述步骤II.2-步骤II.4，直至该行或列栅格数据没有地面点被选中； III.利用稳健曲面拟合方法对步骤II中选取的地面种子点进行核查，将非地面点 剔除； IV.对地面点进行渐进选择； 经过以上步骤I-步骤III中的低异常点剔除、地面种子点选择以及核查后，原始点 云被分为三部分，即非地面点集NonGP、地面点集GP以及候选点集CP； 其中，非地面点集NonGP中仅包括低异常点，地面点集GP中包括地面种子点，候选 点集CP中既包括地面点也包括非地面点； 下面从候选点集CP中选择地面点，候选点集CP中剩余的点即为非地面点集NonGP； 以薄板样条TPS作为插值函数，对地面点渐进选择的具体步骤为： IV.1.对于候选点集中任一待分类点P，利用kd树从地面点集GP选择点P的k个邻近 点； IV.2.以步骤IV.1中k个邻近点为已知点，利用TPS插值求算点P高程值； IV.3.对点P的观测值和估算值做差，求算点P的插值误差e； IV.4.如果插值误差e小于阈值t，则点P标记为地面点； 阈值t的计算公式为：t＝t0 c*slope； 其中，t0为初始误差阈值，c为坡度系数，slope为点P处地形坡度值； 利用薄板样条TPS基本方程求算slope的值，计算公式为： 其 中： 式中，fx表示点P处函数值在x方向偏导数，fy表示点P处函数值在y方向偏导数； αj表示第j个基函数系数； (x,y)表示点P的平面坐标，(xj,yj)表示点P的第j个邻近地面种子点平面坐标； rj表示点P与第j个邻近地面种子点的欧式距离，ε为大于0的常数； IV.5.重复上述步骤IV.1-步骤IV.4，直至候选点集CP中所有待分类点分类完毕； IV.6.将选中的地面点从候选点集CP中删除，利用选中的地面点更新地面点集GP； IV.7.重复上述步骤IV.1-步骤IV.6，直至候选点集CP中没有地面点被选中，则候 选点集CP中剩余的点标记为非地面点集NonGP。 优选地，步骤I具体为： I.1 .将原始点云的高程值分段，即将整个高程值的范围从左到右均匀分成m个间 隔Ci，i＝1,2,…,m，每个间隔的高程跨度L＝(Hmax-Hmin)/m； 其中，Hmax表示原始点云的高程最大值，Hmin表示原始点云的高程最小值； 6 CN 111598780 A 说　明　书 4/12 页 I.2.计算每个间隔Ci中包含的点数； I.3.如果第一个间隔C1中点数超过设定阈值n，则原始点云中认定没有低异常点； 否则，从间隔C1到间隔Ck中选择点数小于n的最右间隔，并将该最右间隔记为Cr； I.4.计算最右间隔Cr中包含的点中高程最大值MAX，并将该高程该最大值MAX作为 区分原始点云中正常点和低异常点的阈值； 将高程值大于MAX的点自动分类为正常点，否则自动分类为低异常点。 优选地，步骤II中，在原始点云研究区域覆盖的网格尺寸等于原始点云中点平均 距离； 其中，点平均距离是指所有相邻点距离的平均值。 优选地，步骤III中，利用稳健曲面拟合方法对地面种子点进行核查的具体过程如 下： 对于步骤II选取的地面种子点中任一点R，利用kd树选中点R周围的k个邻近点，并 对这些邻近点进行平面拟合，然后计算点R到该拟合平面的距离； 如果该距离大于设定阈值rth，则该点R为非地面点，应从地面种子点中剔除。 优选地，步骤III中，在进行平面拟合时，利用最小截断二乘法拟合切平面； 最小截断二乘法的目标函数为： 式中，(r2)(i)表示第i个邻近点到拟合平面的距离平方，且(r2)(1)≤(r2)(2)≤…≤ (r2)(k)，h＝k/2。 本发明具有如下优点： 如上所述，本发明提出了一种用于机载LiDAR点云的地形自适应插值滤波方法，该 方法利用高程直方图自动去除低异常点，以降低人工交互工作量；借助一维离散样条插值 方法选取尽量多的地面种子点，提升初始地面参考面精度；直接利用已选取的地面种子点 估算待分类点的高程，避免使用各种尺度DEM带来的空间位置误差影响；构建自适应滤波阈 值，以降低不同地形场景区点云错分和漏分的水平。本发明利于提升不同场景下点云滤波 的精度。 附图说明 图1为本发明实施例中适用于机载LiDAR点云的地形自适应插值滤波方法的流程 框图。 图2为本发明实施例中基于高程直方图的低异常点剔除示意图。 图3为本发明实施例中按照行方式选择地面种子点的流程示意图。 图4为本发明实施例中最小截断二乘法与传统最小二乘法拟合平面的结果比较示 意图。 图5为本发明实施例中对地面点渐进选择的流程示意图。 图6为空间位置误差对滤波结果影响示意图。 图7为高密度点云区正射影像图。 图8为高密度点云DEM山体阴影图。 图9为本发明方法和PMF、PTD、MHF滤波方法在样区3的DEM示意图。 7 CN 111598780 A 说　明　书 5/12 页