CN111369436A - 一种顾及多元地形特征的机载LiDAR点云抽稀方法 - Google Patents

Abstract

本发明公开了一种顾及多元地形特征的机载LiDAR点云抽稀方法，属于测绘技术领域，其首先获取机载LiDAR点云数据，计算出激光点云坐标；然后对激光点进行局部地形的二次曲面LM(Levenberg‑Marquardt)算法拟合，实现高程标准差、坡度、高斯曲率和粗糙度四类地形特征参数的快速提取；最后基于主成分分析将四类地形特征参数线性构建为多元地形特征复杂度模型，依据获取的点云复杂度及抽稀准则，实现点云的高精度抽稀。本发明通过这种方法，实现了机载LiDAR点云的高精度、强鲁棒性抽稀，有效解决了海量点云在数据显示、处理、存储和传输过程中占用太多内存资源，使得处理器速度缓慢且运算效率不高的问题。

Description

一种顾及多元地形特征的机载LiDAR点云抽稀方法

技术领域

本发明属于测绘技术领域，具体涉及一种顾及多元地形特征的机载LiDAR点云抽稀方法。

背景技术

机载LiDAR系统，是一种集激光测距仪、数字航空摄影、全球定位系统(GPS)和惯性导航系统(INS)等多种尖端技术于一体的空间测量系统，是一种主动式的三维信息获取技术。该系统以飞机为载体，通过高速激光扫描测量的方式，可以高效率、高精度、自动化、直接性的获取地形三维数据信息。与传统的测量方法相比，机载LiDAR测量基本不受光照和天气条件的限制，并且具有周期短、速度快、成本低、效率高等优势，充分满足地形三维建模与更新的时效性要求，具有较高的研究价值和重要的现实意义。

机载LiDAR采集的点云数据往往是海量的，海量点云数据直接影响数据处理效率、数据存储速度，增加了数据交互的难度。机载LiDAR采集的点云数据是海量的，由于点云密度非常高，可达到每平方米十几个甚至上百个点，尤其对于大区域地形测量，所得数据量非常大。对于复杂地形区域，海量点云有利于地形的精确表达，但是对于地形平坦、地形变化趋势较小区域，海量点云存在大量冗余数据，这些冗余数据不仅会消耗大量的内存空间，还会增加点云数据后处理的难度。

点云抽稀是应用研究的重要环节。在具体实际应用研究中，很多工程不需要如此高密度的点云数据，庞大的数据量会加大应用研究的难度，同时也会影响工作效率。在满足工程实际应用精度的前提下，对海量点云进行合理抽稀，既可以保证工程应用精度，又可以去除冗余数据，有利于点云数据的存储及处理。因此，机载LiDAR点云抽稀算法研究对于资源勘探、城市规划、农业开发、水利工程、土地利用、环境监测、交通通讯、防震减灾及国家重点建设项目等方面，都具有重要的实际意义和良好的应用前景。

发明内容

针对现有技术中存在的上述技术问题，本发明提出了一种顾及多元地形特征的机载LiDAR点云抽稀方法，设计合理，解决了现有技术大都是基于单一地形特征实现的简单抽稀，不能综合全面地实现复杂地形的精确抽稀，该抽稀方法能够在充分保留地形地貌特征的前提下有效实现点云抽稀，对于不同复杂地形的机载LiDAR点云数据均具有较好的适用性及鲁棒性。

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

一种顾及多元地形特征的机载LiDAR点云抽稀方法，包括以下步骤：

步骤1：获取机载LiDAR点云数据，计算出激光点云坐标；

步骤2：对激光点进行局部地形的二次曲面拟合，采用LM算法迭代参数寻优，获得最优化结果下的地形拟合参数，并以地形拟合模型为基础，实现高程标准差、坡度、高斯曲率和粗糙度四类地形特征参数的提取；

步骤3：基于主成分分析，将四类地形特征变量因子线性构建为多元地形特征复杂度模型，依据获取的点云复杂度及抽稀准则，实现点云抽稀。

优选地，在步骤2中，高程标准差、坡度、高斯曲率和粗糙度四类地形特征参数的提取步骤如下：

步骤2.1：进行二次曲面LM算法拟合；

假设以r为半径，搜索给定半径区域内的所有激光点，由该区域内的所有激光点拟合形成的二次曲面模型如公式(1)所示：

f(x,y)＝ax²+by²+cxy+dx+ey+f (1)；

则二次曲面拟合的目标函数为：

其中，Q为目标函数；(x_j,y_j,z_j)为点集N_j(j＝1,2,…,k)内的激光点在局部坐标系中的地理坐标；a、b、c、d、e和f为二次曲面拟合参数；

步骤2.2：提取高程标准差参数；

假设点集N_j(j＝1,2,…,k)内的激光点的高程值为H_j，则高程标准差σ_H如公式(3)所示：

其中，H_j为激光点高程值，

为局部区域内激光点的高程均值，k为激光点个数；

步骤2.3：提取坡度参数；

坡度反映了地形的倾斜程度，是地形曲面在该点处的切平面与水平面的夹角，将坡度

分为x和y方向，则有

其中

p代表坡度

在x方向的高程差与水平距离的比值，q代表坡度

在y方向的高程差与水平距离的比值；

步骤2.4：提取高斯曲率参数；

激光点具有最大曲率C_max和最小曲率C_min，对应的曲率半径分别为R_max和R_min，R_max和R_min分别是公式(5)关于曲率半径R的两个根：

(rt-s²)R²+h[2pqs-(1+p²)t-(1+q²)r]R+h⁴＝0 (5)；

其中，

根据微分几何原理，高斯曲率C_G的计算如公式(6)所示：

步骤2.5：提取粗糙度参数；

地形粗糙度是地表的实际面积与投影面积之间的比值；由于地表的实际面积难以获取，故采用拟合的曲面面积代替地表实际面积进行计算；故地形粗糙度K_r由公式(7)和公式(8)计算：

其中，S_S为二次拟合曲面面积，S_E为曲面的投影平面面积，(x₁,y₁)、(x₂,y₂)、(x₃,y₃)、(x₄,y₄)分别是投影平面上点的水平坐标。

优选地，在步骤3中，多元地形特征复杂度模型构建及点云抽稀方法步骤如下：

步骤3.1：为实现不同地形特征参数间量纲的统一，需要对四类地形特征参数进行数据标准化处理，假设地形特征数据为：

标准化公式如(9)所示：

其中，x_ij ^*表示地形特征参数的标准化值；x_ij表示地形特征参数的原始值；

表示地形特征参数的均值，

n为样本数量；

步骤3.2：计算地形特征参数的协方差矩阵或相关系数矩阵；

由标准化后的矩阵向量得到各地形特征参数之间的相关系数矩阵R，作为主成分分析矩阵；r_wj(w，j＝1,2,3,4)为标准化后的各地形特征参数x_w和x_j的相关系数，且r_wj＝r_jw，根据公式(10)，计算各地形特征参数之间的相关系数矩阵R：

其中，

步骤3.3：解特征方程|λE-R|＝0，计算特征值λ_w及对应的特征向量e_w(w＝1,2,3,4)，使||e_w||＝1，即

e_wj表示地形特征参数j在第w个主成分中的系数值；将特征值按照大小顺序排列，λ₁≥λ₂≥λ₃≥λ₄≥0；根据公式(12)、(13)，计算主成分贡献率和累计贡献率：

步骤3.4：根据公式(14)，计算主成分载荷系数；

步骤3.5：根据累计贡献率的大小确定主成分的个数并分析主成分的含义；

步骤3.6：根据公式(15)，对保留的各主成分进行综合，计算各地形特征参数的综合系数f_j：

其中，f_j表示地形特征参数j的综合系数；m为保留的主成分个数；e_wj表示地形特征参数j在第w个主成分中的系数值；S_w表示第w个主成分的贡献率；

步骤3.7：构建多元地形特征复杂度模型，如公式(16)所示：

其中，T表示构建的多元地形特征复杂度模型；f_j为综合系数；t_j表示第j类地形特征因子；

步骤3.8：计算所有激光点的多元地形特征复杂度T，将其作为点云的属性计入点信息；依据抽稀的目标尺度，将地形划分为具有统一分辨率的规则格网，并将点云垂直投射至相对应的格网单元内；计算每个格网的区域地形特征复杂度T_m，T_m为该格网内所有激光点的地形复杂度均值；

步骤3.9：确定抽稀阈值；

选择典型的三种常见地形：平坦地形、一般复杂地形和特别复杂地形，分别计算出三种地形的平均地形特征复杂度：T₁、T₂和T₃，将min(T₁,T₂)和min(T₂,T₃)作为抽稀阈值；

步骤3.10：依据抽稀准则，完成LiDAR点云抽稀；

点云抽稀准则为：

当0<＝T_m<min(T₁,T₂)时，每个格网只保留一个特征点，该点为格网内高程均值点Z_m，其水平坐标为格网中心点坐标；

当min(T₁,T₂)<＝T_m<min(T₂,T₃)时，每个格网保留两个特征值点，分别为格网内高程最大值点Z_max和高程最小值点Z_min；

当T_m>＝min(T₂,T₃)时，每个格网保留三个特征值点，分别为Z_max、Z_min和Z_m点，遍历完成所有格网的判断与选择，实现LiDAR点云抽稀。

本发明所带来的有益技术效果：

本发明提出了顾及多元地形特征的机载LiDAR点云抽稀方法，与现有技术相比，本发明利用二次曲面LM算法对激光点进行局部地形拟合，实现了高程标准差、坡度、高斯曲率和粗糙度四类地形特征参数的快速提取；基于主成分分析算法将四类地形特征参数线性构建为多元地形特征复杂度模型，并依据抽稀准则实现了点云的高精度抽稀；有效解决了海量点云在数据显示、处理、存储和传输过程中占用太多内存资源，使得处理器速度缓慢且运算效率不高的问题。

附图说明

图1为本发明顾及多元地形特征的机载LiDAR点云抽稀方法的流程图。

图2为本发明中坡度、曲率和粗糙度等地形特征参数示意图。

图3为本发明中激光点云抽稀示意图。

图4为本发明中顾及多元地形特征的机载LiDAR点云抽稀方法的详细流程图。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

本发明提供了一种顾及多元地形特征的机载LiDAR点云抽稀方法，其流程如图1所示。

包括以下步骤：

步骤1：获取机载LiDAR点云数据，计算出激光点云坐标。

采用机载LiDAR系统获取激光点云数据；经过姿态改正、多源数据融合、坐标系转换后，得到有效激光点云坐标。

步骤2：二次曲面拟合，快速提取四类地形特征参数。

不同的地形特征参数反映了不同的地形变化程度，地形特征提取为点云抽稀提供了有力的数据支撑。为综合全面地实现复杂地形的精确抽稀，以二次曲面LM算法拟合模型为基础，实现高程标准差、坡度、高斯曲率和粗糙度四类地形特征参数的快速提取。其中坡度、曲率和粗糙度等地形特征参数见图2。

进一步的实施例中，步骤2具体包括如下步骤：

步骤2.1：进行二次曲面LM算法拟合；

假设以r为半径，搜索给定半径区域内的所有激光点，由该区域内的所有激光点拟合形成的二次曲面模型如公式(1)所示：

f(x,y)＝ax²+by²+cxy+dx+ey+f (1)；

则二次曲面拟合的目标函数为：

其中，Q为目标函数；(x_j,y_j,z_j)为点集N_j(j＝1,2,…,k)内的激光点在局部坐标系中的地理坐标；a、b、c、d、e和f为二次曲面拟合参数；

步骤2.2：提取高程标准差参数；

假设点集N_j(j＝1,2,…,k)内的激光点的高程值为H_j，则高程标准差σ_H如公式(3)所示：

其中，H_j为激光点高程值，

为局部区域内激光点的高程均值，k为激光点个数；

步骤2.3：提取坡度参数；

坡度反映了地形的倾斜程度，是地形曲面在该点处的切平面与水平面的夹角，将坡度

分为x和y方向，则有

其中

p代表坡度

在x方向的高程差与水平距离的比值，q代表坡度

在y方向的高程差与水平距离的比值；

步骤2.4：提取高斯曲率参数；

激光点具有最大曲率C_max和最小曲率C_min，对应的曲率半径分别为R_max和R_min，R_max和R_min分别是公式(5)关于曲率半径R的两个根：

(rt-s²)R²+h[2pqs-(1+p²)t-(1+q²)r]R+h⁴＝0 (5)；

其中，

根据微分几何原理，高斯曲率C_G的计算如公式(6)所示：

步骤2.5：提取粗糙度参数；

地形粗糙度是地表的实际面积与投影面积之间的比值；由于地表的实际面积难以获取，故采用拟合的曲面面积代替地表实际面积进行计算；故地形粗糙度K_r由公式(7)和公式(8)计算：

其中，S_S为二次拟合曲面面积，S_E为曲面的投影平面面积，(x₁,y₁)、(x₂,y₂)、(x₃,y₃)、(x₄,y₄)分别是投影平面上点的水平坐标。

具体实施时，基于LM算法的二次曲面拟合模型具有连续性和真实性，且稳定性强、效率高，为四类地形特征参数的快速提取提供了保障。

步骤3：构建多元地形特征复杂度模型，实现点云抽稀。

由于提取的四类地形特征参数均为定量值，且各特征参数间具有一定的相关性和相对独立性，故基于主成分分析将四类地形特征参数线性构建为多元地形特征复杂度模型，并设计了合理的抽稀准则，以此实现顾及多元地形特征的激光点云抽稀。激光点云抽稀示意图见图3，顾及多元地形特征的机载LiDAR点云抽稀方法的详细流程见图4。

进一步的实施例中，步骤3具体包括如下步骤：

步骤3.1：为实现不同地形特征参数间量纲的统一，需要对四类地形特征参数进行数据标准化处理，假设地形特征数据为：

标准化公式如(9)所示：

其中，x_ij ^*表示地形特征参数的标准化值；x_ij表示地形特征参数的原始值；

表示地形特征参数的均值，

n为样本数量；

步骤3.2：计算地形特征参数的协方差矩阵或相关系数矩阵；

由标准化后的矩阵向量得到各地形特征参数之间的相关系数矩阵R，作为主成分分析矩阵；r_wj(w，j＝1,2,3,4)为标准化后的各地形特征参数x_w和x_j的相关系数，且r_wj＝r_jw，根据公式(10)，计算各地形特征参数之间的相关系数矩阵R：

其中，

步骤3.3：解特征方程|λE-R|＝0，计算特征值λ_w及对应的特征向量e_w(w＝1,2,3,4)，使||e_w||＝1，即

e_wj表示地形特征参数j在第w个主成分中的系数值；将特征值按照大小顺序排列，λ₁≥λ₂≥λ₃≥λ₄≥0；根据公式(12)、(13)，计算主成分贡献率和累计贡献率：

步骤3.4：根据公式(14)，计算主成分载荷系数；

步骤3.5：根据累计贡献率的大小来确定主成分的个数并分析主成分的含义；

步骤3.6：根据公式(15)，对保留的各主成分进行综合，计算各地形特征参数的综合系数f_j：

其中，f_j表示地形特征参数j的综合系数；m为保留的主成分个数；e_wj表示地形特征参数j在第w个主成分中的系数值；S_w表示第w个主成分的贡献率。

步骤3.7：构建多元地形特征复杂度模型；

其中，T表示构建的多元地形特征复杂度模型；f_j为综合系数；t_j表示第j类地形特征因子。

步骤3.8：计算所有激光点的多元地形特征复杂度T，将其作为点云的属性计入点信息；依据抽稀的目标尺度，将地形划分为具有统一分辨率的规则格网，并将点云垂直投射至相对应的格网单元内；计算每个格网的区域地形特征复杂度T_m，T_m为该格网内所有激光点的地形复杂度均值；

步骤3.9：确定抽稀阈值；

选择典型的三种常见地形：平坦地形、一般复杂地形和特别复杂地形，分别计算出三种地形的平均地形特征复杂度为T₁、T₂和T₃，确定min(T₁,T₂)和min(T₂,T₃)作为抽稀阈值；

步骤3.10：依据抽稀准则，完成LiDAR点云抽稀；

设计点云抽稀准则为：当0<＝T_m<min(T₁,T₂)时，由于区域地形起伏较小，每个格网只保留一个特征点，该点为格网内高程均值点Z_m，其水平坐标为格网中心点坐标；当min(T₁,T₂)<＝T_m<min(T₂,T₃)时，区域地形起伏与褶皱程度较为明显，地形特征相对较大，故每个格网可保留两个特征值点，分别为格网内高程最大值点Z_max和高程最小值点Z_min；当T_m>＝min(T₂,T₃)时，区域地形特征特别明显，起伏与褶皱程度非常大，此时每个格网须保留三个特征值点，分别为Z_max、Z_min和Z_m点。遍历完成所有格网的判断与选择，实现LiDAR点云抽稀。

具体实施时，构建的多元地形特征复杂度模型可以综合评判地形特征，为点云抽稀提供有力的质量保障；通过合理设计抽稀阈值和抽稀准则，综合全面地实现了复杂地形的精确抽稀。该方法对于不同复杂地形的机载LiDAR点云数据均具有较好的适用性及鲁棒性。

综上所述，本发明提供了顾及多元地形特征的机载LiDAR点云抽稀方法，方法包括：在获得激光点坐标的前提下，进行局部地形的二次曲面LM算法拟合，实现了高程标准差、坡度、高斯曲率和粗糙度四类地形特征参数的快速提取；基于主成分分析将四类地形特征参数线性构建为多元地形特征复杂度模型，并依据抽稀准则实现点云的高精度抽稀。本发明有效解决了海量点云在数据显示、处理、存储和传输过程中占用太多内存资源，使得处理器速度缓慢且运算效率不高的问题，实现了机载LiDAR点云的高精度、强鲁棒性抽稀。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (3)

1.一种顾及多元地形特征的机载LiDAR点云抽稀方法，其特征在于：包括以下步骤：

步骤1：获取机载LiDAR点云数据，计算出激光点云坐标；

步骤2：对激光点进行局部地形的二次曲面拟合，采用LM算法迭代参数寻优，获得最优化结果下的地形拟合参数，并以地形拟合模型为基础，实现高程标准差、坡度、高斯曲率和粗糙度四类地形特征参数的提取；

步骤3：基于主成分分析，将四类地形特征变量因子线性构建为多元地形特征复杂度模型，依据获取的点云复杂度及抽稀准则，实现点云抽稀。

2.根据权利要求1所述的顾及多元地形特征的机载LiDAR点云抽稀方法，其特征在于：在步骤2中，高程标准差、坡度、高斯曲率和粗糙度四类地形特征参数的提取步骤如下：

步骤2.1：进行二次曲面LM算法拟合；

假设以r为半径，搜索给定半径区域内的所有激光点，由该区域内的所有激光点拟合形成的二次曲面模型如公式(1)所示：

f(x,y)＝ax²+by²+cxy+dx+ey+f (1)；

则二次曲面拟合的目标函数为：

其中，Q为目标函数；(x_j,y_j,z_j)为点集N_j(j＝1,2,…,k)内的激光点在局部坐标系中的地理坐标；a、b、c、d、e和f为二次曲面拟合参数；

步骤2.2：提取高程标准差参数；

假设点集N_j(j＝1,2,…,k)内的激光点的高程值为H_j，则高程标准差σ_H如公式(3)所示：

其中，H_j为激光点高程值，

为局部区域内激光点的高程均值，k为激光点个数；

步骤2.3：提取坡度参数；

坡度反映了地形的倾斜程度，是地形曲面在该点处的切平面与水平面的夹角，将坡度

分为x和y方向，则有

其中

p代表坡度

在x方向的高程差与水平距离的比值，q代表坡度

在y方向的高程差与水平距离的比值；

步骤2.4：提取高斯曲率参数；

激光点具有最大曲率C_max和最小曲率C_min，对应的曲率半径分别为R_max和R_min，R_max和R_min分别是公式(5)关于曲率半径R的两个根：

(rt-s²)R²+h[2pqs-(1+p²)t-(1+q²)r]R+h⁴＝0 (5)；

其中，

根据微分几何原理，高斯曲率C_G的计算如公式(6)所示：

步骤2.5：提取粗糙度参数；

地形粗糙度是地表的实际面积与投影面积之间的比值；由于地表的实际面积难以获取，故采用拟合的曲面面积代替地表实际面积进行计算；故地形粗糙度K_r由公式(7)和公式(8)计算：

其中，S_S为二次拟合曲面面积，S_E为曲面的投影平面面积，(x₁,y₁)、(x₂,y₂)、(x₃,y₃)、(x₄,y₄)分别是投影平面上点的水平坐标。

3.根据权利要求1所述的顾及多元地形特征的机载LiDAR点云抽稀方法，其特征在于：在步骤3中，多元地形特征复杂度模型构建及点云抽稀方法步骤如下：

步骤3.1：为实现不同地形特征参数间量纲的统一，需要对四类地形特征参数进行数据标准化处理，假设地形特征数据为：

标准化公式如(9)所示：

其中，x_ij ^*表示地形特征参数的标准化值；x_ij表示地形特征参数的原始值；

表示地形特征参数的均值，

n为样本数量；

步骤3.2：计算地形特征参数的协方差矩阵或相关系数矩阵；

由标准化后的矩阵向量得到各地形特征参数之间的相关系数矩阵R，作为主成分分析矩阵；r_wj(w，j＝1,2,3,4)为标准化后的各地形特征参数x_w和x_j的相关系数，且r_wj＝r_jw，根据公式(10)，计算各地形特征参数之间的相关系数矩阵R：

其中，

步骤3.3：解特征方程|λE-R|＝0，计算特征值λ_w及对应的特征向量e_w(w＝1,2,3,4)，使||e_w||＝1，即

e_wj表示地形特征参数j在第w个主成分中的系数值；将特征值按照大小顺序排列，λ₁≥λ₂≥λ₃≥λ₄≥0；根据公式(12)、(13)，计算主成分贡献率和累计贡献率：

步骤3.4：根据公式(14)，计算主成分载荷系数；

步骤3.5：根据累计贡献率的大小确定主成分的个数并分析主成分的含义；

步骤3.6：根据公式(15)，对保留的各主成分进行综合，计算各地形特征参数的综合系数f_j：

其中，f_j表示地形特征参数j的综合系数；m为保留的主成分个数；e_wj表示地形特征参数j在第w个主成分中的系数值；S_w表示第w个主成分的贡献率；

步骤3.7：构建多元地形特征复杂度模型，如公式(16)所示：

其中，T表示构建的多元地形特征复杂度模型；f_j为综合系数；t_j表示第j类地形特征因子；

步骤3.8：计算所有激光点的多元地形特征复杂度T，将其作为点云的属性计入点信息；依据抽稀的目标尺度，将地形划分为具有统一分辨率的规则格网，并将点云垂直投射至相对应的格网单元内；计算每个格网的区域地形特征复杂度T_m，T_m为该格网内所有激光点的地形复杂度均值；

步骤3.9：确定抽稀阈值；

选择典型的三种常见地形：平坦地形、一般复杂地形和特别复杂地形，分别计算出三种地形的平均地形特征复杂度：T₁、T₂和T₃，将min(T₁,T₂)和min(T₂,T₃)作为抽稀阈值；

步骤3.10：依据抽稀准则，完成LiDAR点云抽稀；

点云抽稀准则为：

当0<＝T_m<min(T₁,T₂)时，每个格网只保留一个特征点，该点为格网内高程均值点Z_m，其水平坐标为格网中心点坐标；

当min(T₁,T₂)<＝T_m<min(T₂,T₃)时，每个格网保留两个特征值点，分别为格网内高程最大值点Z_max和高程最小值点Z_min；

当T_m>＝min(T₂,T₃)时，每个格网保留三个特征值点，分别为Z_max、Z_min和Z_m点，遍历完成所有格网的判断与选择，实现LiDAR点云抽稀。

Images