一种移动二维激光雷达的点云强度改正方法与流程

本发明涉及激光点云强度改正方法，尤其是一种地面移动二维激光雷达的点云强度改正方法。

背景技术：

激光强度表征为被扫描目标对发射激光光束后向散射回波的光功率。由于受到扫描仪特性、大气传输、探测器与放大电路噪声和扫描物体几何形状等因素响，激光强度值不能作为反射特性被直接提取，而是需要对其值进行改正，以消除各种因素的干扰。通过激光强度改正，扫描目标的属性信息被充分挖掘，可以开展利用目标特性进行分类的研究工作。

目前研究表明，改正后的激光在森林和地表资源勘探、历史文物修复、路面和隧道破损检测、冰川表面分类等领域有着重要的价值。激光强度信息常用来辅助提高地物点云分类精度。针对不同的扫描方式，现有的激光强度改正方法，主要分为以下两种：

(1)地面三维激光扫描(terrestriallaserscanning，tls)强度改正。地面三维激光扫描除了能获取高精度的三维空间信息外，还能记录目标的激光强度。由于地面三维激光雷达扫描距离相对较近，大气等环境因素可以忽略不计，激光强度改正主要考虑扫描距离和入射角因素。通过对实测数据的拟合，可以得到激光强度值与距离、入射角之间的关系。同时，通过反求运算，得到目标反射率。

(2)机载激光扫描(airbornelaserscanning，als)强度改正。机载激光扫描是依托全球定位系统、惯性测量装置、激光扫描仪获得高精度的三维点云数据，从激光强度在传播过程中的衰减机理出发，对激光传播距离、大气衰减、入射角因素进行改正，消除干扰因素对激光强度值的影响。

als和tls都采用成本较高的三维激光雷达，移动二维激光雷达成本低、数据冗余度小，可根据测量点的帧序号和帧内编号直接索引点云。目前还没有方法来消除距离和入射角对移动二维激光雷达强度的影响，因此，针对移动二维激光雷达强度改正方法有待深入研究。

技术实现要素：

本发明的目的是针对移动二维激光雷达获得的同质区域点云强度差异，提出了一种地面移动二维激光的点云强度改正方法。

本发明的技术方案是：

本发明提供一种移动二维激光雷达的点云强度改正方法，该方法包括以下步骤：

s1、利用二维激光雷达提取标准漫反射板在不同扫描距离、入射角度下的强度数据，分别建立扫描距离和强度值的函数，以及入射角度和强度值的函数；

s2、采用最小二乘法解算改正参数，根据拟合优度和均方根误差，确定模型阶数；

s3、采用移动二维激光扫描系统对实际场景进行移动扫描，获取实际场景的点云距离ρ、扫描角度θ，计算激光落脚点到激光原点的入射角α；

s4、建立基于扫描距离和入射角度对激光强度数据的改正模型，得到改正后的激光强度is。

进一步地，步骤s1具体为：

s1-1、设置入射角度、激光雷达的距离范围[rmin，rmax]以及激光雷达的距离间隔xm；采用激光雷达以同样的入射角度获取标准漫反射板在前述距离范围内，按照固定间隔调整距离所接收的各距离下的激光强度数据，建立扫描距离和强度值的拟合函数f2(r)；

其中：r表示激光雷达到标准漫反射板的距离；rsp表示拟合函数分段点，选取前述各距离下的激光强度最大值处的距离值，n和m分别表示f21(r)和f22(r)的多项式阶数，[a0，a1，…an]和[b0，b1，…bm]分别表示多项式的系数；

s1-2、设置激光雷达的距离、激光雷达的入射角度范围[0°，90°]以及标准漫反射板的旋转间隔；采用激光雷达以同样的扫描距离获取标准漫反射板在前述入射角度范围内，按照固定旋转间隔调整角度所接收的各入射角度下的激光强度数据，建立入射角度和强度值的拟合函数f3(cosα)：

f3(cosα)＝c0+c1cosα+…+ck(cosα)^k

其中：α表示激光雷达的入射角度；[c0，c1，…ck]表示f3(cosα)多项式系数，k表示多项式阶数。

进一步地，步骤s1-1中：对于激光雷达获取的激光强度数据，取一帧数据的中间点作为当前距离下的强度值；步骤s1-2中：对于激光雷达获取的激光强度数据，取一帧数据的中间点作为当前入射角下的强度值。

进一步地，步骤s1-2中：标准漫反射板的旋转间隔为10°。

进一步地，步骤2具体为：

s2-1、根据步骤1所建立的扫描距离和强度值的函数，以及入射角度和强度值的函数，通过最小二乘法计算，获取前述两个函数中的多项式系数[a0，a1，a2，…，an]^t、[b0，b1，b2，…，bm]^t、[c0，c1，c2，…，ck]^t；

s2-2、采用以下公式计算均方根误差(rmse)和拟合优度(r-square)，获取前述两个函数中的模型阶数：

其中：rmse1表示参与拟合的距离-强度数据点(r⁽ⁱ⁾，i2⁽ⁱ⁾)与拟合函数f2(r⁽ⁱ⁾)的均方根误差；i表示距离-强度数据点的编号；m表示距离-强度数据点的总数；f2(r⁽ⁱ⁾)表示编号为i的数据点的扫描距离和强度函数的拟合值，i2⁽ⁱ⁾表示编号为i的数据点对应距离下的激光雷达记录的原始强度值；

rmse2表示参与拟合的入射角-强度的数据点(α^j，i3^(j))与拟合函数f3(cosα^(j))的均方根误差；j表示入射角-强度数据点的编号，n表示入射角-强度数据点的总数；α表示激光入射角；f3(cosα^(j))表示编号为j的数据点的入射角和强度函数的拟合值，i3^(j)表示编号为j的数据点对应入射角下的激光雷达记录的原始强度值；

r-square1表示f2(r⁽ⁱ⁾)的拟合优度，r-square2表示f3(cosα^(j))的拟合优度；μ2表示原始强度i2⁽ⁱ⁾均值，μ3表示原始强度i3^(j)均值；

设定多项式阶数n、m和k的范围，n、m、k∈[0，7]，计算预设范围内所有系数对应的均方根误差(rmse)和拟合优度(r-square)，选择拟合优度接近1且均方根误差最小的阶数。

进一步地，步骤3具体为：

步骤31、以激光雷达初始位置为坐标原点o，建立直角坐标系o-xyz，x轴方向为激光雷达在滑台上的运动方向，y轴方向为激光雷达扫描深度方向，z轴方向为被扫描目标与地面垂直的高度方向，利用二维激光雷达采集第j帧内yz平面的极坐标数据其中n为激光雷达扫描点数，坐标变换公式：

其中，ρ(i′，j′)和θ(i′，j′)分别为激光雷达第j′帧内第i′个激光点的距离和扫描角度，x(i′，j′)表示第j′帧内第i′个激光点在x方向上的坐标，y(i′，j′)表示第j′帧内第i′个激光点在深度方向的坐标，z(i′，j′)表示第j′帧内第i′个激光落脚点在高度方向的坐标，δt表示激光雷达的扫描周期，v表示激光雷达移动速度v；

步骤32、计算点云的法向量，对于点云中每个点p(i′，j′)，通过最邻近算法建立邻域集，即获取最相近的k个相邻点，然后为这些点拟合一个最小二乘意义上的局部平面，即对于邻域中每个点qr，建立协方差矩阵m：

其中：k为最邻近点的个数，r为最邻近点的编号，是经过最邻近点的三维质心；

采用下述公式计算法向量；

其中：l为特征值的编号，取值1、2、3，λl为第l个特征值，为第l个特征值对应的特征向量，；通过主成分分析，即对协方差矩阵m进行特征值分解，将m矩阵中最小特征值λmin对应的特征向量作为拟合平面的法向量，这里记作

步骤33、采用下式得到激光落脚点p(i′，j′)((x(i′，j′)，y(i′，j′)，z(i′，j′))到激光原点o(x(i′，j′)，0，0)的距离r和入射角α；

式中，表示入射激光向量，表示法向量。

进一步地，步骤4具体为：建立距离改正模型如下：

建立改正模型，得到距离和入射角改正后的激光强度is：

其中；i′表示未改正的即待改正的距离所对应的原始强度值，rm表示需要改正到的标准距离，αm表示需要改正到的标准入射角，isr表示距离改正后的激光强度值，r表示激光落脚点p到激光原点o的距离，α表示入射角。

本发明的有益效果：

本发明的移动二维激光的点云强度改正方法，该方法从影响强度差异的因素入手，根据同类地物反射强度相同这一原则，利用二维激光雷达提取标准漫反射板在不同扫描距离、角度下的强度数据；然后采用移动二维激光扫描系统对实际场景进行移动扫描，获取实际场景的点云距离、扫描角度、强度等信息，再根据改正模型消除距离和入射角对强度值的影响，使改正后的同质区域强度趋于统一。

本发明的其它特征和优点将在随后具体实施方式部分予以详细说明。

附图说明

通过结合附图对本发明示例性实施方式进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施方式中，相同的参考标号通常代表相同部件。

图1是二维激光雷达数据采集坐标系

图2是扫描时刻与几何关系示意图

图3是距离与强度关系图

图4是入射角与强度关系图

图5是不同距离改正后强度图

图6是不同入射角改正后强度图

图7是扫描墙面示意图

图8是扫描墙面点云强度改正后伪彩色图、强度分布直方图。

图8(a)原始强度伪彩色图、图8(b)原始强度分布直方图、图8(c)入射角改正后强度伪彩色图、图8(d)入射角改正强度分布直图、图8(e)距离改正后强度伪彩色图、图8(f)距离改正强度分布直方图。

具体实施方式

下面将参照附图更详细地描述本发明的优选实施方式。虽然附图中显示了本发明的优选实施方式，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施方式所限制。

本发明提供一种移动二维激光雷达的点云强度改正方法，该方法包括以下步骤：

s1、利用二维激光雷达提取标准漫反射板在不同扫描距离、入射角度下的强度数据，分别建立扫描距离和强度值的函数，以及入射角度和强度值的函数，具体为：

s1-1、设置入射角度、激光雷达的距离范围[rmin，rmax]以及激光雷达的距离间隔xm；采用激光雷达以同样的入射角度获取标准漫反射板在前述距离范围内，按照固定间隔调整距离所接收的各距离下的激光强度数据，建立扫描距离和强度值的拟合函数f2(r)；

其中：r表示激光雷达到标准漫反射板的距离；rsp表示拟合函数分段点，选取前述各距离下的激光强度最大值处的距离值，n和m分别表示f21(r)和f22(r)的多项式阶数，[a0，a1，…an]和[b0，b1，…bm]分别表示多项式的系数；

s1-2、设置激光雷达的距离、激光雷达的入射角度范围[0°，90°]以及标准漫反射板的旋转间隔(间隔优选10°)；

采用激光雷达以同样的扫描距离获取标准漫反射板在前述入射角度范围内，按照固定旋转间隔调整角度所接收的各入射角度下的激光强度数据，建立入射角度和强度值的拟合函数f3(cosα)：

f3(cosα)＝c0+c1cosα+…+ck(cosα)^k

其中：α表示激光雷达的入射角度；[c0，c1，…ck]表示f3(cosα)多项式系数，k表示多项式阶数。

进一步地，步骤s1-1中：对于激光雷达获取的激光强度数据，取一帧数据的中间点作为当前距离下的强度值；步骤s1-2中：对于激光雷达获取的激光强度数据，取一帧数据的中间点作为当前入射角下的强度值。

s2、采用最小二乘法解算改正参数，根据拟合优度和均方根误差，确定模型阶数，具体为：

s2-1、根据步骤1所建立的扫描距离和强度值的函数，以及入射角度和强度值的函数，通过最小二乘法计算，获取前述两个函数中的多项式系数[a0，a1，a2，…，an]^t、[b0，b1，b2，…，bm]^t、[c0，c1，c2，…，ck]^t；

s2-2、采用以下公式计算均方根误差(rmse)和拟合优度(r-square)，获取前述两个函数中的模型阶数：

其中：rmse1表示参与拟合的距离-强度数据点(r⁽ⁱ⁾，i2⁽ⁱ⁾)与拟合函数f2(r⁽ⁱ⁾)的均方根误差；i表示距离-强度数据点的编号；m表示距离-强度数据点的总数；f2(r⁽ⁱ⁾)表示编号为i的数据点的扫描距离和强度函数的拟合值，i2⁽ⁱ⁾表示编号为i的数据点对应距离下的激光雷达记录的原始强度值；

rmse2表示参与拟合的入射角-强度的数据点(α^j，i3^(j))与拟合函数f3(cosα^(j))的均方根误差；j表示入射角-强度数据点的编号，n表示入射角-强度数据点的总数；α表示激光入射角；f3(cosα^(j))表示编号为j的数据点的入射角和强度函数的拟合值，i3^(j)表示编号为j的数据点对应入射角下的激光雷达记录的原始强度值；

r-square1表示f2(r⁽ⁱ⁾)的拟合优度，r-square2表示f3(cosα^(j))的拟合优度；μ2表示原始强度i2⁽ⁱ⁾均值，μ3表示原始强度i3^(j)均值；

设定多项式阶数n、m和k的范围，n、m、k∈[0，7]，计算预设范围内所有系数对应的均方根误差(rmse)和拟合优度(r-square)，选择拟合优度接近1且均方根误差最小的阶数。

s3、采用移动二维激光扫描系统对实际场景进行移动扫描，获取实际场景的点云距离ρ、扫描角度θ，计算激光落脚点到激光原点的入射角α；

步骤31、以激光雷达初始位置为坐标原点o，建立直角坐标系o-xyz，x轴方向为激光雷达在滑台上的运动方向，y轴方向为激光雷达扫描深度方向，z轴方向为被扫描目标与地面垂直的高度方向，利用二维激光雷达采集第j′帧内yz平面的极坐标数据其中n为激光雷达扫描点数，坐标变换公式：

其中，ρ(i′，j′)和θ(i′，j′)分别为激光雷达第j′帧内第i′个激光点的距离和扫描角度，x(i′，j′)表示第j′内第i′个激光点在x方向上的坐标，y(i′，j′)表示第j′帧内第i′个激光点在深度方向的坐标，z(i′，j′)表示第j′帧内第i′个激光落脚点在高度方向的坐标，δt表示激光雷达的扫描周期，v表示激光雷达移动速度v；

步骤32、计算点云的法向量，对于点云中每个点p(i′，j′)，通过最邻近算法建立邻域集，即获取最相近的k个相邻点，然后为这些点拟合一个最小二乘意义上的局部平面，即对于邻域中每个点qr，建立协方差矩阵m：

其中：k为最邻近点的个数，r为最邻近点的编号，是经过最邻近点的三维质心；

采用下述公式计算法向量

其中：l为特征值的个数，取值0、1、2，λmin为最小特征值，为最小特征值对应的第l个特征向量；通过主成分分析，即对协方差矩阵m进行特征值分解，将m矩阵中最小特征值对应的特征向量作为平面的法量；

步骤33、采用下式得到激光落脚点p(i′，j′)((x(i′，j′)，y(i′，j′)，z(i′，j′))到激光原点o(x(i′，j′)，0，0)的距离r和入射角α；

式中，表示入射激光向量，表示法向量。

s4、建立基于扫描距离和入射角度对激光强度数据的改正模型，得到改正后的激光强度is；

其中；i′表示未改正的即待改正的距离所对应的原始强度值，rm表示需要改正到的标准距离，αm表示需要改正到的标准入射角，is表示距离改正后的激光强度值，r表示激光落脚点p到激光原点o的距离，α表示入射角。

具体实施时：实验采用一种精确匀速同步滑台搭载二维激光雷达的强度校正装置，该装置主要有滑台移动控制单元和激光数据采集单元两部分组成，激光传感器与计算机通过usb串口通信。所选用的hokuyo公司生产的utm-30lx二维激光扫描仪，扫描范围[0°，270°]，最大测量距离30m，角度分辨率0.25°，扫描周期25ms。利用激光雷达获取标准漫反射板同质区域点云数据，采用反射率为50％的标准漫反射板为扫描对象，距离实验是将标准漫反射板正对于激光雷达，连续采集10帧，取每帧数据的中间点的强度数据并求其平均值，在近距离处设置18个站点，扫描距离范围是0.1m-1.8m，每隔0.1m设置一个站点对漫反射板进行扫描，在较远处设置21个站点，扫描范围是1.8m-14.4m，每隔0.6m设置一个站点进行扫描；入射角实验是在三个固定站点下分别完成的，每个站点分别距离漫反射板1.2m、1.8m、2.4m处，同样使漫反射板正对激光雷达，对于每个站点，在0°到80°的范围内，使用测角仪每隔10°水平旋转标准漫反射板并进行一次扫描，同样连续采集10帧数据，取每帧中间点的强度数据并求其平均值。在多项式拟合过程中，通过均方根误差(rmse)和确定系数(r-square)来评估多项式阶数对拟合结果的影响，f21(r)、f22(r)、f3(cosθ)的模型的确定系数分别为0.9913、0.9985、0.9979，非常接近于1，且均方根误差控制在10.2、5.3、6.4，因此选取的最佳模型阶数分别为n＝3、m＝5、k＝1。在对改正模型进行验证时，激光雷达安装在距离地面0.15m处，设计匀速同步滑台沿平行于墙面方向运动，移动速度为10mm/s，实验选取距离为3m处的墙面，选取的墙面点云个数为609248，待改正墙面的距离范围在[3m，4m]，入射角范围在[0°，60°]，实验发现改正后的强度值较为集中的分布在3000左右。

不同距离和角度下获取的墙面点云与移动二维激光雷达的几何关系如图2所示，用本发明阐述的方法得到改正前后墙面的强度变化如图8所示。(为便于显示，对点云强度进行彩色化处理)，发现改正前同一块墙面的不同区域分布的强度值十分散乱，改正后分布的强度值较集中，并呈现高斯分布。

以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。