基于可信度信息的点云数据粗差定位和剔除方法与流程

本发明涉及数据误差处理，具体涉及一种基于可信度信息的点云数据粗差定位和剔除方法。

背景技术：

三维激光扫描测量技术是目前测绘工程中表示物体三维模型的常用方法之一，利用三维激光扫描测量技术可以扫描获取物体上大量的点云数据，包括物体的几何表面特征信息，诸如位置、法向量、颜色等，有助于便捷地重建出物体的数字模型，并可对其进行绘制、编辑等操作。通常在获取物体三维点云数据时，会因为设备本身的缺陷或者人为扰动使得扫描数据中不可避免地存在粗差数据。除此之外，无人机摄影测量技术获取的点云数据是由影像匹配获取的，由于影像获取的过程中，相机存在着镜头畸变，在软件处理中这种镜头畸变造成的误差并不能被完全消除，就会造成点云异位，也就造成粗差点的出现；除了镜头畸变之外，影像的质量及影像匹配误差的影响都会造成粗差点的出现。因此，不论哪种方式获取的高密度点云数据，都不可避免地存在粗差点。

如果三维点云数据点中带有粗差，就会影响模型参数化和重建的精度和效果，因此需要首先对点云数据进行粗差定位和剔除。现有的粗差探测和剔除技术进行点云数据的粗差处理，存在以下缺点：

(1)数据计算量巨大：针对散乱点云数据的粗差探测方法，由于散乱点云数据中只储存了被扫描物体表面的坐标信息，而没有明显的拓扑关系，为了建立点云数据中各点之间的空间拓扑与邻近关系，需要经过多次迭代才能实现，计算量非常大，效率低。

(2)粗差之间相互遮蔽：在对散乱点云数据进行邻域搜索时，需要首先设置点云数据邻域的大小，若点云数据邻域太大，会导致多个观测粗差或异常值处于相同的邻域内，这就会产生粗差或异常值之间的互相遮蔽。

(3)只能探测到较大的粗差：由于存在着粗差和异常值的相互遮蔽现象，致使某些数值较小的粗差根本无法探测出，并且异常值越大，对粗差的遮蔽越严重。

针对上述问题，现有技术中存在基于未确知滤波的点云数据粗差探测方法。根据点云数据粗差的特点，对趋势面滤波方法进行了改进，提出了基于未确知滤波的趋势面粗差探测方法。但在该技术中，只是考虑了粗差数据的偶然性和孤立性特点，并没有对点云数据的邻域搜索方法进行改进，也没有提出能够解决粗差遮蔽问题的邻域设置方式，因此，仍然存在计算量大、效率低，以及部分粗差之间相互遮蔽的缺陷。基于上述原因，该方法只适用于对数据量比较小、表面比较简单的点云数据进行粗差的定位和剔除。

技术实现要素：

发明目的：本申请的目的在于提供一种基于可信度信息的点云数据粗差定位和剔除方法，用于解决现有的点云数据处理技术中因为计算量大、粗差剔除伴随着特征信息丢失，从而导致粗差剔除不够快速、准确和适用性不强的问题。

技术方案：本发明一方面提供了一种基于可信度信息的点云数据粗差定位和剔除方法，包括：

(1)点云数据采集：通过测量设备采集待测物体的点云数据；

(2)点云数据预处理：对所测点云数据进行预处理，包括：点云数据降维、点云数据规则化；

(3)点云数据邻域确定和曲面拟合，包括：点云数据邻域自适应划分、邻域内点云数据曲面拟合；

(4)置信度计算，包括：根据迭代计算方法，计算出点云数据中所有点到其对应邻域内二次曲面的距离，选取点云数据中某点邻域中距离的平均值和方差，来计算该点的置信度，并据此计算出点云数据中其它所有点的置信度；

(5)粗差定位和剔除：根据点云数据中各点的置信度判断该点是否为粗差，若点云数据中某点判断为粗差，则以置信度作为权重，对该点所处邻域内的数据进行加权求和，得其邻域内点云数据的数学期望，以剔除定位到的粗差点。

进一步地，步骤(2)中，对点云数据进行预处理包括：

点云数据降维：利用坐标变换法，固定点云中的数据点，将其所处的坐标系进行平移或旋转，选取适当的坐标系表示点云数据的坐标，从而固定数据点坐标的其中一维，实现点云数据降维；

点云数据规则化：对点云数据按照行列进行规则化，形成按行列顺序排列的规则点云数据。

进一步地，步骤(3)中，

点云数据邻域自适应划分包括：将点云数据进行邻域划分，找到每个数据点p(i,j)所处的k×k最近邻域；其中i，j分别为该数据对应的行列号，且i＝1…m，j＝1…n，m为点云数据行总数，n为点云数据列总数，k为数据点p(i,j)所处邻域的大小，即邻域的行列总数；

邻域内点云数据曲面拟合包括：对点云数据中每个数据点所处的k×k邻域进行曲面拟合，即对该邻域内的数据利用最小二乘的方法构造二次曲面。

进一步地，在点云数据邻域自适应划分中，按照行列顺序确定邻域边界，根据以下原则进行邻域划分：

当i＝1时，取第1、2行作为点云数据的行邻域；当j＝1时，取第1、2列作为点云数据的列邻域；

当时，取第1至2i-1行作为点云数据的行邻域；当时，取第1至2j-1列作为点云数据的列邻域；

当时，取第至行作为点云数据的行邻域；当时，取第至列作为点云数据的列邻域；

当时，取第2i-m至m行作为点云数据的行邻域；当时，取第2j-n至n列作为点云数据的列邻域；

当i＝m时，取第m-1、m行作为点云数据的行邻域；当j＝n时，取第n-1、n列作为点云数据的列邻域。

进一步地，邻域内点云数据曲面拟合包括：

将数据点p(i,j)邻域内数据点lp(i,j)(r,s)(r＝1,2,…k；s＝1,2,…k)，在局部邻域内进行二次曲面拟合，拟合方程的一般形式为：

s(u,v)＝au²+buv+cv²+du+ev

拟合的目标函数为：

其中，u和v为曲面参数，a、b、c、d、e为曲面方程的系数，应用奇异值分解法可获得拟合曲面的最小二乘解。

进一步地，步骤(4)包括：

(41)求得数据点p(i,j)邻域内的每一个数据点沿法线方向到对应邻域内的二次曲面的距离dp(i,j)(r,s)(r＝1,2,…k；s＝1,2,…k)，计算得出邻域内数据点到邻域内二次曲面距离的平均值和均方差λ分别为：

其中，r,s分别为邻域的行列数；k为邻域的大小，即行列数；dp(i,j)(r,s)为邻域中点到二次曲面的距离；

(42)以数据点p(i,j)到拟合曲面的距离dp(i,j)和邻域内其它数据点到邻域曲面的距离dp(i,j)(r,s)(r＝1,2,…k；s＝1,2,…k)构造置信度邻域n，

n＝{d||dp(i,j)(r,s)-dp(i,j)|<δ,δ>0}(r＝1,2,…k；s＝1,2,…k)

取δ＝ωλ，即n＝{d||dp(i,j)(r,s)-dp(i,j)|＜ωλ}，其中ω为大于0的正整数，取值大小可根据点云数据的粗差剔除效果确定；若数据点p(i,j)的邻域集合中，各点到拟合曲面的距离处于闭区间[dp(i,j)-ωλdp(i,j)+ωλ]的个数为ξi,j，则数据点p(i,j)的置信度为：

其中，ξi,j表示p(i,j)点邻域中处于置信度区间数据点的个数；表示数据点p(i,j)是否为粗差的可信程度。

进一步地，步骤(5)中，设定粗差定位所需的置信度阈值，若数据点p(i,j)的置信度大于置信度阈值，则判定该数据点对应的观测值为非粗差值；否则，定位该数据点对应的观测值为粗差。

进一步地，对包括数据点p(i,j)的k×k邻域内所有的观测值求数学期望，即将其所处邻域内所有数据点的置信度作为权重进行加权求和，得出观测值的数学期望e(pi,j)，如下式：

其中，pi,j(r,s)(r＝1,2…k；s＝1,2…k)为p(i,j)点邻域中第r行和第s列的数据点；为该数据点对应的置信度。

另一方面，本发明还提供以一种计算机设备，包括存储器和处理器，存储器中存储计算机程序，计算机程序适于被处理器加载并执行基于可信度信息的点云数据粗差定位和剔除方法。

有益效果：与现有技术相比，本发明基于点云数据粗差探测中存在的计算复杂、效率低，存在粗差之间以及粗差和异常值之间相互遮蔽的问题，提高计算效率和粗差定位准确性。通过将点云数据进行降维和规则化处理，降低了计算的复杂程度；通过点云数据邻域的自适应设置，解决了点云数据的粗差的相互遮蔽问题；通过点云数据的可信度信息，对观测数据是否为粗差进行了准确判断，该发明不仅可以提高计算效率，而且可以有效地区分粗差和异常值，避免了粗差和异常值的混淆。

附图说明

图1为本申请的点云数据粗差定位和剔除方法流程图；

图2为本申请实施例中数据点到曲面距离计算示意图；

图3为工程实例一中待测量区域的三维扫描图；

图4为工程实例一中扫描测量的点云数据图；

图5为工程实例一中原始点云数据所在坐标系示意图；

图6为工程实例一中降维后点云数据所在坐标系示意图；

图7为工程实例一中规则化后的点云数据图；

图8为工程实例一中粗差定位的结果示意图；

图9为工程实例二中待测量区域的无人机摄影图；

图10为工程实例二中航测点云数据图；

图11为工程实例二中第一区域点云数据放大图；

图12为工程实例二中第二区域点云数据放大图；

图13为工程实例二中第三区域的点云数据放大图；

图14为工程实例二中摄影测量点云数据粗差定位结果；

图15为工程实例二中第一区域点云数据粗差定位结果放大图；

图16为工程实例二中第二区域点云数据粗差定位结果放大图；

图17为工程实例二中第三区域点云数据粗差定位结果放大图。

具体实施方式

下面结合附图和实施例对本发明做进一步描述：

本申请提供了一种基于可信度信息的点云数据粗差定位和剔除方法，如图1所示，包括：

s101点云数据采集：通过测量设备采集待测物体的点云数据。

具体地，在数据采集过程中可采用三维激光扫描仪进行扫描测量或采用无人机进行倾斜摄影测量。

s102点云数据预处理：对所测点云数据进行预处理，包括：点云数据降维、点云数据规则化。

具体地：

①点云数据降维：利用坐标变换法，固定点云数据中的数据点，将其所处的坐标系进行平移或旋转，选取适当的坐标系，以最简单的形式表示点云数据的坐标，例如，可以使空间中的平面与新坐标系的某个坐标轴垂直或者落在新坐标系的某个坐标平面上。如此，在新坐标系下，原先点的坐标将有一维为固定值，从而固定数据点坐标的其中一维，实现点云数据降维。

②点云数据规则化：对点云数据按照行列进行规则化，形成按行列顺序排列的规则点云数据。

s103点云数据邻域确定和曲面拟合，包括：点云数据邻域自适应划分、邻域内点云数据曲面拟合。

①点云数据邻域自适应划分：将点云数据进行邻域划分，找到每个数据点p(i,j)所处的k×k最近邻域；其中i，j分别为该数据对应的行列号，且i＝1…m，j＝1…n，m为点云数据行总数，n为点云数据列总数，k为数据点p(i,j)所处邻域的大小，即为邻域的行数和列数。在本实例中，根据以下原则进行邻域划分：

当i＝1时，取第1、2行作为点云数据的行邻域；当j＝1时，取第1、2列作为点云数据的列邻域；

当时，取第1至2i-1行作为点云数据的行邻域；当时，取第1至2j-1列作为点云数据的列邻域；

当时，取第至行作为点云数据的行邻域；当时，取第j-k2至列作为点云数据的列邻域；

当时，取第2i-m至m行作为点云数据的行邻域；当时，取第2j-n至n列作为点云数据的列邻域；

当i＝m时，取第m-1、m行作为点云数据的行邻域；当j＝n时，取第n-1、n列作为点云数据的列邻域。

②邻域内点云数据曲面拟合：对点云数据中每个数据点所处的k×k邻域进行曲面拟合，即对该邻域内的数据利用最小二乘的方法构造二次曲面。

将数据点p(i,j)邻域内数据点lp(i,j)(r,s)(r＝1,2,…k；s＝1,2,…k)，在局部邻域内进行二次曲面拟合，拟合方程的一般形式为：

s(u,v)＝au²+buv+cv²+du+ev

拟合的目标函数为：

其中，u和v为曲面参数，a、b、c、d、e为曲面方程的系数，应用奇异值分解法可获得拟合曲面的最小二乘解。

s104置信度计算，包括：根据迭代计算方法，计算出点云数据中所有点到其对应邻域内二次曲面的距离，选取点云数据中某点邻域中距离的平均值和方差，来计算该点的置信度，并据此计算出点云数据中其它所有点的置信度。

具体地，求得数据点pi,j邻域内的每一个数据点沿法线方向到对应邻域内的二次曲面距离dp(i,j)(r,s)(r＝1,2,…k；s＝1,2,…k)，计算得出邻域内数据点到邻域内二次曲面距离的平均值和均方差λ。

①算法原理：迭代法计算点到曲面的最小距离，如图2所示，参数曲面s(u,v)外一点q，利用网格对曲面进行划分，网格上距离q最近的点p0为初始点，曲面上q的投影点为pc，点q到曲面的最小距离可以表示为d，n是曲面在pc处的法矢量，于是有：

经过几次移动，当|δu|＜ε，|δv|＜ε时，就可以认为p0到达pc点了。该方法不仅可以提高计算速度，并且可以克服曲面边界处的迭代发散。

②求解δu、δv

用su、sv分别点乘式(1)两边，即

因为所以su·n·d＝sv·n·d＝0，设e＝su·su，f＝su·sv，g＝sv·sv，可得：

解上述方程组，可得δu、δv。得出δu、δv后，令u+δu→u；v+δv→v，经过迭代重复计算后，当|δu|＜ε，|δv|＜ε时，即可计算出q点到曲面的距离d。

③邻域内数据点到二次曲面距离的平均值由下式表示：

④邻域内数据点到二次曲面距离的均方差λ，由下式表示：

其中，r,s分别为邻域的行列数；k为邻域的大小，即行列总数；dp(i,j)(r,s)为邻域中点到二次曲面的距离。

s105粗差定位和剔除：根据点云数据中各点的置信度判断该点是否为粗差，若点云数据中某点判断为粗差，则以置信度作为权重，对该点所处邻域内的数据进行加权求和，得其邻域内点云数据的数学期望，以剔除定位到的粗差点。

具体地，

①点云数据的置信度邻域

假设点云数据表示形式如下：

points＝{p(i,j),i＝1,2…m；j＝1,2,…n}(6)

其中，points表示扫描测量的点云数据集合，m表示点云数据的行数，n表示点云数据的列数。p(i,j)＝(xi,j,yi,j,zi,j)表示点云数据的空间坐标。

以数据点p(i,j)到拟合曲面的距离dp(i,j)和邻域内其它数据点到邻域曲面的距离dp(i,j)(r,s)(r＝1,2,…k；s＝1,2,…k)构造置信度邻域n，

n＝{d||dp(i,j)(r,s)-dp(i,j)|<δ,δ>0}(r＝1,2,…k；s＝1,2,…k)(7)

取δ＝ωλ，即n＝{d||dp(i,j)(r,s)-dp(i,j)|＜ωλ}，其中ω为大于0的正整数，取值大小可根据点云数据的粗差剔除效果确定。

②点云数据置信度确定

若数据点p(i,j)的邻域集合中，各点到拟合曲面的距离处于闭区间[dp(i,j)-ωλdp(i,j)+ωλ]的个数为ξi,j，则数据点p(i,j)的置信度为：

上式中，ξi,j表示数据点p(i,j)邻域中处于置信度邻域的数据点的个数；表示p(i,j)是否为粗差的可信程度；若ξi,j数值较大，则认为p(i,j)的可信度就大，p(i,j)为非粗差测值的可能性就大，为粗差的可能性就小；反之当ξi,j数值较小，则认为p(i,j)的可信度就小，p(i,j)为非粗差测值的可能性就小，为粗差的可能性就大。

③粗差定位

根据式(8)对观测数据信息的描述，可以从每个测值相应的可信度反应其为粗差的可能性大小，从而实现粗差定位。

在本实施例中，设定粗差定位所需的置信度阈值，若数据点p(i,j)的置信度大于置信度阈值，则判定该数据点对应的观测值为非粗差值；否则，定位该数据点对应的观测值为粗差。

④粗差剔除

为了对定位到的粗差进行剔除或减弱粗差对观测值的影响，对包括数据点p(i,j)的k×k邻域内所有的观测值求数学期望，即将其所处邻域内所有数据点的置信度作为权重进行加权求和，得出各观测值的数学期望e(pi,j)，如下式：

其中，pi,j(r,s)(r＝1,2…k；s＝1,2…k)为p(i,j)点邻域中第r行和第s列的数据点；为该点对应的置信度。

如果点数据p(i,j)中带有粗差，其可信度就是零或接近于零，通过加权求和的方法可以对观测值中的粗差进行剔除。

另一方面，本发明还提供以一种计算机设备，包括存储器和处理器，存储器中存储计算机程序，计算机程序适于被处理器加载并执行上述基于可信度信息的点云数据粗差定位和剔除方法。

下面引用两个工程实例对本发明进行进一步说明：两实例对地面三维激光扫描点云和无人机倾斜摄影测量点云数据分别进行处理。实例一利用地面三维激光扫描仪对小范围的土方进行精细扫描，得到密集的点云数据。实例二利用无人机倾斜摄影测量技术获取大范围建筑场地点云数据。

工程实例一：

如图3所示，对一建筑物施工场地，采用美国faro公司focus^3d地面三维激光扫描仪对该原始场地进行扫描测量，扫描共布设了4个测站，并利用faroscene软件对点云数据进行了配准处理，得到的点云数据如图4所示。

由于场地上杂草和部分植被的存在使得扫描测量的点云数据存在大量的粗差，虽然手工已经剔除掉了部分较大的粗差，但仍然有少量的粗差存在点云数据当中。这些粗差和点云数据混杂在一起，并且测量地形高低起伏不平，很难利用人工的方式进行探测和剔除，利用本发明对点云数据中粗差进行探测，步骤如下：

(1)点云数据降维：如图5所示，由于地面三维激光扫描测量获取的点云数据，是以扫描仪中心为坐标原点的局部坐标系数据，不同测站获取的点云数据经过配准处理后会发生平移和旋转变化。因此，得到的原始点云整体数据在三维空间中呈倾斜状态，不能再保持水平或竖直状态，给点云的数据处理带来了困难。

为了降低数据处理的难度，将原始点云数据的坐标系变为以点云数据中心为坐标原点，以向上的方向为z轴方向，以点云数据的一条边为x轴的坐标系统。假设原始点云坐标系为o1-x1y1z1，变换后点云数据的坐标系为o2-x2y2z2，即：

o1＝[21.857437.357711.3693]

o2＝[38.005446.383711.9613]

由坐标系o-x1y1z1和o-x2y2z2的参数矩阵，计算出两个坐标系的旋转矩阵r和平移参数t：

t＝[16.14809.02600.5920]

利用计算的旋转和平移矩阵，将原始点云数据进行转换，变换后的点云数据坐标系图6所示。

(2)点云数据规则化：对测量的散乱点云数据进行内插，得到按行、列顺序排列的规则点云数据，如图7所示。

(3)点云数据邻域划分：根据点云数据密度，判断点云数据邻域大小k为15，进行邻域划分。以其中的p(54,56)点为例，按照邻域划分规则，其邻域的范围应为以p(54,56)为中心的15行和15列的点云数据。

(4)构造观测量：以p(54,56)点邻域内的225个点云数据进行曲面拟合，然后进行如下计算：

①点到曲面的距离：按照上述迭代计算方法的原理，计算出p(54,56)点到拟合曲面的距离为0.24m。

②邻域内距离平均值的计算：

③邻域内距离的均方差计算：

(5)可信度计算：以p(54,56)点的距离s54,56构造置信度邻域，并取ω＝2，即{d||d-d54,56|≤2λ}。该点的邻域距离集合中，处于闭区间[d54,56-2λd54,56+2λ]的个数为ξ54,56＝18，则p(54,56)点的可信度为：

(6)粗差判断：为了对点云数据的粗差进行判断，把置信度的参数阈值设置为0.2，当点云的置信度小于0.2时，则判断该观测值为粗差，由此可以判断p(54,56)为粗差点。

(7)粗差剔除：对p(54,56)邻域内的点云数据计算加权平均值，得到将粗差点p(54,56)剔除后的高程。

按照上述操作步骤，就可以定位到点云数据中的所有粗差，并对探测到的粗差进行剔除，如图8所示。

工程实例二：

利用大疆无人机对一大范围场地进行倾斜摄影测量，根据设计的航线飞行了6个架次，共获取了1300多张照片，选取图9中的三个区域进行测量。采用smart3d(contextcapture)软件对倾斜影像进行处理，得到该区域的密集点云数据，如图10所示，图11、图12和图13为第一、第二和第三区域放大后的点云数据。

按照本申请的点云数据粗差剔除的方法，对本测量区域点云数据进行粗差定位，具体的定位结果如图14所示，将第一、第二、第三区域的粗差定位结果放大，如图15、16、17所示。

由图8和15、16、17中的定位结果可知，本发明的方法不仅可以提高计算效率，而且可以准确地定位出点云数据中的粗差，较好地解决了粗差和异常值之间相互遮蔽的问题，即使对于大范围、高密度的复杂地形表面点云数据，也能较好地对粗差和地形变化特征进行区分，探测到数值较小的粗差。

本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd-rom、光学存储器等)上实施的计算机程序产品的形式。

本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。