一种基于稀疏编码的LiDAR点云数据修复方法与流程

本发明属于激光雷达三维点云数据修复领域，更进一步是一种基于稀疏编码理论的LiDAR(Light Detection And Ranging)点云数据修复方法。

背景技术：

随着三维激光扫描技术的日渐成熟和相关仪器设备的广泛普及，使得利用三维激光扫描技术对目标进行精确三维重建成为研究热点。在地形测绘方面，对于城市及危险地区的精细地形测绘，常规的方法往往显得无能为力，而利用机载LiDAR系统进行非接触式的测量能够直接获取高精度的三维数据，对物体没有限制，且具有实时性强、精度高、扫描速度快等特点，能够无接触地完成对复杂危险地区的精细地形测量。目前的LiDAR以单点扫描方式为主，通过离散点云数据内插后处理来拟合生成连续地形表面，现有的点云修补方法一般通过对点云数据进行“插值”实现，如克里金(Kriging)插值法、样条(Spline)插值法、离散平滑(Discrete Smooth Interpolation)插值法等，这种点云修补方法可能造成点云数据的失真及变形，使点云数据丢失特征；另外，地形地貌比较复杂，在利用LiDAR进行地面的三维点云数据获取的过程中，常因为高吸收性地物、航带无重叠度、地物遮挡等原因造成点云数据缺失，数据的缺失将无法生成完整的DSM(Digital Surface Model)和DEM(Digital Elevation Model)模型，影响后续的处理和应用。

技术实现要素：

针对现有技术的不足，本发明设计了一种全新的基于稀疏编码的LiDAR点云数据修复方法，主要是将LiDAR离散点云数据看作缺失了一部分数据的连续‘图像’，这样数据的处理过程就从矢量格式数据的处理变为了栅格格式数据的处理，使得利用稀疏编码强大的自适应学习能力来修复缺失数据成为可能，突破常规的离散点云修补方法多是在插值算法上来改进，容易丢失点云数据细节特征的不足，为LiDAR点云数据的处理提供了一种全新研究思路和方法。

为了解决上述技术问题，本发明提出的一种基于稀疏编码的LiDAR点云数据修复方法，包括如下步骤：

步骤一、利用ArcGIS软件对原始点云数据重采样，获取栅格化后的点云数据形成点云数据矩阵；

步骤二、对所述点云数据矩阵进行掩膜矩阵提取并将掩膜矩阵应用到上述点云数据矩阵，即将掩膜矩阵和点云数据矩阵进行点乘；

步骤三、通过稀疏编码的自适应学习过程，对步骤二得到的点云数据矩阵中缺失的数据进行修补，

式(1)中，λ为拉格朗日乘数，β为掩膜矩阵，y为上述点云数据矩阵，T_ij为二进制矩阵，为字典矩阵，为稀疏向量，为经修复的点云数据矩阵；具体步骤包括：

3-1)建立二进制矩阵T_ij，对输入的点云数据矩阵y及掩膜矩阵β在固定位置进行(6～10)×(6～10)矩阵提取并将提取出的(6～10)×(6～10)矩阵转换为向量的形式；

3-2)建立输出点云数据矩阵保存经修复的点云数据并将输入的点云数据矩阵y直接赋值给输出的点云数据矩阵完成初始化；

3-3)建立字典矩阵保存运算过程中产生的字典并利用外部数据集对字典矩阵初始化；

3-4)建立稀疏向量矩阵保存运算过程中产生的稀疏向量利用OMP算法求解最优化问题完成对稀疏向量矩阵的初始化：

3-5)更新字典矩阵

3-6)更新稀疏向量矩阵；

3-7)对步骤3-5)及步骤3-6)进行迭代，迭代次数为8～20；

步骤四、利用加权平均的方法对经过步骤三后的字典矩阵和稀疏向量矩阵进行运算，对点云数据矩阵中缺失的数据进行修复，输出经修复后的点云数据矩阵。

进一步讲，本发明步骤二中，所述掩膜矩阵的提取是指建立与输入点云数据矩阵大小相同的掩膜矩阵，并且对该掩膜矩阵中的数值二值化，若对应位置为原始数据，则填充数值1，否则填充数值0。

本发明步骤3-4)中利用OMP算法求解下述最优化问题：

式(2)中，为T_ij在输出点云数据矩阵某一固定位置提出的向量，为稀疏向量矩阵Coeffs与位置相对应的一列，为字典矩阵，T_ijβ∈R^n×1为T_ij在掩膜矩阵β的某一固定位置提出的向量，C是由所提取出的(6～10)×(6～10)矩阵大小决定的常数，代表点乘。

本发明步骤四中，加权平均公式为：

其中，λ是拉格朗日乘数，是输出的点云数据矩阵，y代表输入的点云数据矩阵，为字典矩阵，为稀疏向量。

与现有技术相比，本发明的有益效果是：

本发明方法解决了由于地形地貌比较复杂或者因为高吸收性地物、航带无重叠度、地物遮挡等原因造成点云数据缺失，可以较多的保留细节特征，生成完整的DSM(Digital Surface Model)和DEM(Digital Elevation Model)模型，更好地实现了后续的处理及应用。

附图说明

图1是本发明LiDAR点云修复方法流程图；

图2是本发明实施例2中预处理后的LiDAR三维点云图；

图3是本发明实施例2中经修复后的LiDAR三维点云图；

图4是与实施例2中所处理的区域对应的航拍影像图。

具体实施方式

下面结合附图和具体实施例对本发明技术方案作进一步详细描述，所描述的具体实施例仅对本发明进行解释说明，并不用以限制本发明。

实施例1：本发明提出的一种基于稀疏编码的LiDAR点云数据修复方法，如图1所示，包括如下步骤：

步骤一、利用ArcGIS软件对原始点云数据重采样，获取栅格化后的点云数据形成点云数据矩阵；根据机载雷达的扫描间隔，本实施例中，利用ArcGIS 10.2软件中的Conversion Tools对原始点云数据进行了矢量数据的栅格化。为了保留点云数据中的“孔洞”，所选用的空值填充方法为“None”。将栅格化后的点云数据称为输入点云数据矩阵y∈R^n×n。

步骤二、对所述点云数据矩阵进行掩膜矩阵提取，掩膜提取是指对输入点云数据矩阵y进行二值化，提取出对应的掩膜矩阵β∈R^n×n。首先构建与输入点云数据矩阵y大小相同的掩膜矩阵β，然后根据输入点云数据矩阵y的值对相应位置的掩膜矩阵β进行填充，如果对应点为原始点云数据，则填充数值1，否则填充数值0。将掩膜矩阵β应用到上述点云数据矩阵y，即将掩膜矩阵β和点云数据矩阵y进行点乘，使输入的点云数据矩阵y的孔洞值为0，减少后续运算量。

步骤三、主要作用是执行稀疏编码算法，即式(1)，对步骤二得到的点云数据进行修复，稀疏编码是自适应学习方法，需要进行一定的迭代步骤来使修复的效果达到最优，具体内容如下：

式(1)中，λ为拉格朗日乘数，β为掩膜矩阵，y为上述点云数据矩阵，T_ij为二进制矩阵，为字典矩阵，为稀疏向量，为经修复的点云数据矩阵；

3-1)建立二进制矩阵T_ij，对输入的点云数据矩阵y及掩膜矩阵β在固定位置[i,j]进行(6～10)×(6～10)矩阵提取并将提取出的(6～10)×(6～10)矩阵转换为向量的形式；其中，[i,j]代表矩阵的左上角坐标(在矩阵坐标[10,10]处提取出8×8的点云数据矩阵，将其转换为64×1的向量)。

3-2)建立输出点云数据矩阵保存经修复的点云数据并将输入的点云数据矩阵y直接赋值给输出的点云数据矩阵完成初始化；

3-3)建立字典矩阵的目的是为了存储在运算过程中产生的字典。字典的初始化过程通过导入外置的数据集来实现，经初始化后其中为字典矩阵的一列。

3-4)建立稀疏向量矩阵Coeffs保存运算过程中产生的稀疏向量稀疏向量的初始化过程通过OMP(orthogonal matching pursuit)算法求解最优化问题完成对稀疏向量矩阵的初始化：

式(2)中，为T_ij在输出点云数据矩阵某一固定位置提出的向量，为稀疏向量矩阵Coeffs与位置相对应的一列，为字典矩阵，T_ijβ∈R^n×1为T_ij在掩膜矩阵β的某一固定位置提出的向量，C是由所提取出的(6～10)×(6～10)矩阵大小决定的常数，代表点乘。

上述OMP的执行步骤如下：

输入：字典矩阵向量向量T_ijβ，稀疏度K；

输出：α_ij的K-稀疏的逼近

初始化：建立残差索引集t＝1；

循环执行步骤(1)到(5)；

(1)找出残差r与字典矩阵的列积中最大值所对应的脚标θ，即：

(2)更新索引集Λ_t＝Λ_t-1∪θ_t，记录找到的字典矩阵中重建原子集合

(3)由最小二乘得到

(4)更新残差

(5)判断是否满足t＞K,若满足，则停止迭代；若不满足，则执行步骤(1)。

通过上述公式以及固定字典矩阵就可以对稀疏向量矩阵Coeffs的每一列进行初始化。

3-5)更新字典矩阵

固定稀疏向量矩阵Coeffs，并且对字典矩阵的每一列执行以下步骤；

(1)挑选出稀疏向量矩阵Coeffs中不全为0的列集合所对应的小块集合，也就是说

式(3)中，代表稀疏向量矩阵Coeffs中列元素不全为0的列集合，ω_l代表与之对应的小块集合。

(2)对ω_l中的每一列按照公式(4)计算其误差向量

式(4)中，代表稀疏编码矩阵的列向量的第l个元素。

(3)建立误差矩阵其中E_l的每一列代表

(4)通过公式(5)更新和

式(4)是一个一阶近似方程，可以通过截断SVD(singular value decomposition)进行求解。式(5)中，是的简略表达，(E_l-dα^T)代表ωl的每一列的误差，β^l是一个大小与E_l相同的矩阵，并且β^l的每一列与在坐标[i,j]下对掩膜矩阵β进行移位提取并向量化后的值相同。

3-6)更新稀疏向量矩阵Coeffs，更新稀疏向量矩阵Coeffs的步骤与步骤3-4)相同，不同的是所用字典矩阵为经步骤3-5)更新后的字典矩阵。

3-7)对步骤3-5)及步骤3-6)进行迭代，迭代次数为8～20；本发明的迭代次数优选为10次；

步骤四、利用下述加权平均公式对经过步骤三后的字典矩阵和稀疏向量矩阵进行运算，经加权平均后的结果就是经修补后的点云数据矩阵。所述加权平均公式为：

式(6)中，λ是拉格朗日乘数，是输出的点云数据矩阵，y代表输入的点云数据矩阵，为字典矩阵，为稀疏向量。

实施例2：以从ISPRS(The International Society for Photogrammetry and Remote Sensing)测试项目中下载的点云为例进行点云修复，图2为经预处理后的点云数据，从图2中可以看出点云数据处于缺失状态，并且很难从图2中分辨出不同建筑物之间的高低。图3为按照上述实施例1的方法过程进行修复后的点云数据，图4为区域对应的航拍影像图，将图3与图4进行对比，可以明显的看出建筑物之间的高低，并且完整的保留了细节信息(如倾斜屋脊的高度变化信息、清晰的地物边界，树木、汽车等)。

尽管上面结合附图对本发明进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨的情况下，还可以做出很多变形，这些均属于本发明的保护之内。