地物三维变化检测方法、电子设备及存储介质

1.本发明属于测绘科学与技术领域，尤其涉及一种基于点云密度自适应局部欧式距离的地物三维变化检测方法、电子设备及存储介质。


背景技术：

2.由于点云存在高程方向的信息，能够反映地物几何上的属性，而地物的变化在几何上的反映往往不容易被误判，且更为强烈(地物的新建、拆除等)，因此基于点云的三维变化检测能够较好地获取地物的变化信息。同时，随着近年来传感器技术及无人机的发展，点云的精度和质量越来越高，所包含的信息越来越丰富，获取门槛、成本也有着大幅的降低，因此越来越多的研究者们开始进行三维变化检测相关的研究。
3.根据检测单元区分，三维变化检测主要可分为逐点和基于对象两类检测方法。逐点的三维变化检测通常利用逐点计算的高差、欧式距离或图割方法来提取候选变化区域，在此基础上，利用几何结构、纹理或颜色信息进行后精化处理，例如murakami h,nakagawa k,hasegawa h,et al.change detection of buildings using an airborne laser scanner[j].isprs journal of photogrammetry and remote sensing,1999,54(2-3):148-152.、jung f.detecting building changes from multitemporal aerial stereopairs[j].isprs journal of photogrammetry and remote sensing,2004,58(3-4):187-201.、彭代锋,张永军,熊小东.结合lidar点云和航空影像的建筑物三维变化检测[j].武汉大学学报信息科学版,2015,40(4):462-468.、以及du s,zhang y,qin r,et al.building change detection using old aerial images and new lidar data[j].remote sensing,2016,8(12):1030.。
[0004]
基于对象的三维变化检测方法通常需要对不同时相的数据进行分类，然后对比点云的分类结果从而提取变化区域，例如matikainen l,j,ahokas e,et al.automatic detection of buildings and changes in buildings for updating of maps[j].remote sensing,2010,2(5):1217-1248.、pang s,hu x,wang z,et al.object-based analysis of airborne lidar data for building change detection[j].remote sensing,2014,6(11):10733-10749.、以及qin r,huang x,gruen a,et al.object-based 3-d building change detection on multitemporal stereo images[j].ieee journal of selected topics in applied earth observations and remote sensing,2015,8(5):2125-2137.。
[0005]
逐点三维变化检测方法对点云的质量、配准都有着较高的要求，而且植被区域的季节性变化会导致误检测，使得这类方法受到一定的局限。与逐点的三维变化检测相比，基于对象的三维变化检测对点云的质量、配准要求相对较低，这类方法将分类与变化检测过程割裂，使得最终变化检测结果受到前序分类误差的影响，形成误差累积，检测精度依赖于点云分类的质量。
[0006]
总的来说，目前大部分三维变化检测方法在点云几何精度高、配准误差小且场景
中地物变化显著的情况下，变化检测的精度都能达到很高，但如果点云中存在地物变形，或是存在一定的配准误差，都会导致变化检测存在较多的误检、漏检地物，变化检测精度难以得到保证。


技术实现要素：

[0007]
本发明的目的在于提供一种地物三维变化检测方法、电子设备及存储介质，以解决在点云几何精度不高、配准误差较大时存在的误检、漏检问题。
[0008]
本发明是通过如下的技术方案来解决上述技术问题的：一种地物三维变化检测方法，包括以下步骤：
[0009]
步骤1：读取两个时相的点云，对其进行精配准，并以精配准后两个时相的点云中的其中一个作为参考时相点云，另一个作为待检测时相点云；
[0010]
步骤2：计算所述待检测时相点云、参考时相点云的法向量；
[0011]
步骤3：计算所述待检测时相点云中每个待检测点的法向量与z方向之间的夹角θ；
[0012]
若夹角θ小于角度阈值t
angle
，则搜索待检测点在参考时相点云中xy平面上的k个邻近点，并对该k个邻近点进行聚类，得到第一聚类结果，根据所述第一聚类结果选择合适的聚类簇计算局部欧式距离和局部密度；否则搜索待检测点在参考时相点云三维空间的k个邻近点，并计算局部欧式距离和局部密度；
[0013]
步骤4：根据所述局部密度改善局部欧式距离，得到密度自适应局部欧式距离；
[0014]
步骤5：将所述待检测时相点云中密度自适应局部欧式距离大于距离阈值td的待检测点提取到候选变化区域；
[0015]
步骤6：搜索所述候选变化区域内候选变化点qi在待检测时相点云中xy平面上的k个邻近点，并对该k个邻近点进行聚类，得到第二聚类结果，根据所述第二聚类结果对候选变化点qi进行移出或不移出候选变化区域的处理；
[0016]
步骤7：对经过所述步骤6处理后的候选变化区域内的点进行聚类，并根据每个聚类簇的平均密度计算该聚类簇对应的面积s；若面积s大于面积阈值t
area
，则认为该聚类簇发生变化；否则认为该聚类簇为误检测。
[0017]
优选地，所述步骤1中，采用icp算法对两个时相的点云进行精配准。
[0018]
进一步地，所述步骤2中，计算所述待检测时相点云、参考时相点云的法向量的具体实现过程为：
[0019]
建立所述待检测时相点云的三维kd树以及所述参考时相点云的三维kd树；
[0020]
根据所述三维kd树分别计算所述待检测时相点云、参考时相点云的法向量。
[0021]
进一步地，所述步骤3中，夹角θ的计算公式为：
[0022][0023]
其中，n
x
为待检测点的法向量在x方向上的值，ny为待检测点的法向量在y方向上的值。
[0024]
进一步地，所述步骤3中，选择合适的聚类簇计算局部欧式距离和局部密度的具体实现过程为：
[0025]
条件a：若所述第一聚类结果中聚类簇的数量为1，则将该k个邻近点放入参数计算点簇；
[0026]
条件b：若所述第一聚类结果中聚类簇的数量为2，且该k个邻近点中高程值最大的点的法向量与z方向之间的夹角小于角度阈值t
angle
，则将最远离检测点的聚类簇放入参数计算点簇；
[0027]
若不满足条件a和条件b，则将最近邻检测点的聚类簇放入参数计算点簇；
[0028]
根据所述参数计算点簇中的点计算局部欧式距离；
[0029]
计算所述参数计算点簇的中心，在参考时相点云中搜索kd个邻近点计算局部密度。
[0030]
优选地，所述局部欧式距离的计算公式为：
[0031][0032]
其中，d
loc
为局部欧式距离，kr为参数计算点簇中近邻点数量，为待检测点在参考时相点云的第j个近邻点，(xi,yi,zi)为待检测点的三维坐标，(xj,yj,zj)为点的三维坐标。
[0033]
优选地，所述局部密度的计算公式为：
[0034][0035]
其中，ρ
loc
为局部密度，kd为参考时相点云中近邻点的数量，r
max
为参数计算点簇的中心到kd个邻近点的最大距离。
[0036]
进一步地，所述步骤4中，密度自适应局部欧式距离的计算公式为：
[0037][0038]
其中，为密度自适应局部欧式距离，d
loc
为局部欧式距离，ρ
loc
为局部密度。
[0039]
进一步地，所述步骤6中，根据所述第二聚类结果对候选变化点qi进行移出或不移出候选变化区域的处理的具体实现过程为：
[0040]
若所述第二聚类结果中聚类簇的数量为1，则候选变化点qi不属于地物建模不好区域，不对候选变化点qi进行移出处理；
[0041]
若所述第二聚类结果中聚类簇的数量为2、3或4，则候选变化点qi属于地物建模不好区域，若候选变化点qi的第一顶面点存在且未发生变化，则将候选变化点qi移出候选变化区域；否则不对候选变化点qi进行处理；其中，候选变化点qi的第一顶面点是该k个邻近点中点的法向量与z方向之间的夹角θ小于角度阈值t
angle
且高程值最大的点；
[0042]
若所述第二聚类结果中聚类簇的数量大于4，则候选变化点qi属于地物建模不好区域，若候选变化点qi的第二顶面点发生变化，则不对候选变化点qi进行处理；否则将候选变化点qi移出候选变化区域；其中，候选变化点qi的第二顶面点是该k个邻近点中高程值最
大的点。
[0043]
进一步地，所述步骤7中，每个聚类簇的平均密度计算公式为：
[0044][0045]
其中，为当前聚类簇的平均密度，kc为当前聚类簇中点的数量，(xj,yj,zj)为当前聚类簇中第j个点的三维坐标，(xn,yn,zn)为与当前聚类簇中第j个点对应的最邻近点的三维坐标。
[0046]
进一步地，所述步骤7中，面积s的计算公式为：
[0047][0048]
其中，s为当前聚类簇的面积，为当前聚类簇的平均密度，kc为当前聚类簇中点的数量。
[0049]
本发明还提供一种电子设备，包括存储器和处理器，所述存储器上存储有能够在所述处理器上运行的计算机程序，所述处理器运行所述计算机程序时执行如上所述地物三维变化检测方法的步骤。
[0050]
本发明还提供一种计算机可读存储介质，所述计算机可读存储介质为非易失性存储介质或非瞬态存储介质，其上存储有计算机程序，所述计算机程序被处理器运行时执行如上所述地物三维变化检测方法的步骤。
[0051]
有益效果
[0052]
与现有技术相比，本发明的优点在于：
[0053]
本发明所提供的一种地物三维变化检测方法、电子设备及存储介质，在待检测点与三维空间单个近邻点之间欧式距离的基础上，引入待检测点与xy平面、三维空间上的k个邻近点之间的局部欧式距离，不仅减少了消失地物在边缘部分的误检测，而且大大提高了抗噪性和适应性；利用局部密度改善局部欧氏距离，减小了点云密度不一致以及配准误差导致计算得到的距离大于真实变化距离的影响；对候选变化区域xy平面上的邻近点进行聚类，根据聚类结果判断当前候选变化点是否属于地物建模不好区域，减少了点云变形造成的误检测；利用局部密度计算候选变化区域聚类簇的面积，剔除了面积较小的候选变化区域，提高了变化检测的质量。
[0054]
本发明所述方法具有较好的抗噪性能，更好的变化检测效果，能够大幅减少地物建模变形而造成的误检测；可以在点云几何精度不高、存在一定配准误差的情况下，较为完整、较高精度地提取点云中的变化区域。
附图说明
[0055]
为了更清楚地说明本发明的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一个实施例，对于本领域普通技术人员来说，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0056]
图1是本发明实施例中地物三维变化检测方法流程图；
[0057]
图2是本发明实施例中所使用的检测数据，其中，图(a)为2020年6月的建筑物点
云，图(b)为2020年9月的建筑物点云；
[0058]
图3是本发明实施例中三维变化检测结果，其中图(a)是以图2(b)为参考时相点云，图2(a)为待检测时相点云得到的变化检测结果；图(b)是以图2(a)为参考时相点云，图(b)为待检测时相点云得到的变化检测结果；
[0059]
图4是本发明实施例中所使用的人工标记的精度验证数据。
具体实施方式
[0060]
下面结合本发明实施例中的附图，对本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0061]
下面以具体地实施例对本技术的技术方案进行详细说明。下面这几个具体的实施例可以相互结合，对于相同或相似的概念或过程可能在某些实施例不再赘述。
[0062]
如图1所示，本发明实施例所提供的一种地物三维变化检测方法，包括以下步骤：
[0063]
步骤1：点云的读取与精配准
[0064]
读取两个时相的点云，对两个时相的点云进行精配准，并以精配准后两个时相的点云中的其中一个作为参考时相点云，另一个作为待检测时相点云。如图2所示，图(a)示出2020年6月的建筑物点云，图(b)示出2020年9月的建筑物点云，可以以图(a)为参考时相点云，图(b)为待检测时相点云，也可以以图(b)为参考时相点云，图(a)为待检测时相点云。
[0065]
本实施例，采用icp算法对两个时相的点云进行精配准，icp算法简单且易于实现，在初值较好的情况下，收敛速度和精度较好。采用icp算法对两个时相的点云进行精配准属于现有技术，可参考besl p j,mckay n d.method for registration of 3-d shapes[c]//sensor fusion iv:control paradigms and data structures.spie,1992,1611:586-606.。
[0066]
步骤2：检测时相点云、参考时相点云的法向量的计算，具体实现过程为：
[0067]
步骤2.1：建立待检测时相点云的三维kd树以及参考时相点云的三维kd树；
[0068]
步骤2.2：根据待检测时相点云的三维kd树计算待检测时相点云的法向量，根据参考时相点云的三维kd树计算参考时相点云的法向量。
[0069]
根据三维kd树计算时相点云的法向量为现有技术，可参考holz d,holzer s,rusu r b,et al.real-time plane segmentation using rgb-d cameras[c]//robot soccer world cup.springer,berlin,heidelberg,2011:306-317.。
[0070]
步骤3：局部欧式距离和局部密度的计算
[0071]
计算待检测时相点云中每个待检测点的法向量与z方向之间的夹角θ；若夹角θ小于角度阈值t
angle
，则搜索待检测点在参考时相点云中xy平面上的k个邻近点，并对该k个邻近点进行欧式聚类，得到第一聚类结果，根据第一聚类结果选择合适的聚类簇计算局部欧式距离和局部密度；否则搜索待检测点在参考时相点云三维空间的k个邻近点，并计算局部欧式距离和局部密度。
[0072]
本实施例中，待检测时相点云和参考时相点云的x、y坐标是基于wgs84利用utm投
影得到，z坐标即为wgs84坐标系的高程。
[0073]
本实施例中，待检测点的法向量与z方向之间的夹角θ的计算公式为：
[0074][0075]
其中，n
x
为待检测点的法向量在x方向上的值，ny为待检测点的法向量在y方向上的值。角度阈值t
angle
可以根据需求进行设置，在点云几何精度高，配准较好的情况下，如果对较小变化目标感兴趣，则角度阈值t
angle
可以设置小一点，否则设置大一点，本实施例中，角度阈值t
angle
设为20
°
。
[0076]
为了提高搜索速度，建立待检测时相点云和参考时相点云的二维kd树，在k个邻近点搜索时，基于参考时相点云的二维kd树能够快速搜索到待检测点在参考时相点云中xy平面上的k个邻近点。由于在参考时相点云xy平面搜索得到的邻近点会存在多个高程层次，因此需要对其进行欧式聚类；而在参考时相点云三维空间搜索得到的邻近点只有一个高程层次，无需进行欧式聚类。本实施例中，k为12。
[0077]
本实施例中，选择合适的聚类簇计算局部欧式距离和局部密度的具体实现过程为：
[0078]
步骤3.1：条件a：若第一聚类结果中聚类簇的数量为1，则将该k个邻近点放入参数计算点簇。
[0079]
步骤3.2：条件b：若第一聚类结果中聚类簇的数量为2，且该k个邻近点中高程值最大的点的法向量与z方向之间的夹角小于角度阈值t
angle
，则将最远离检测点的聚类簇放入参数计算点簇。
[0080]
步骤3.3：若不满足条件a和条件b，则将最近邻检测点的聚类簇放入参数计算点簇。
[0081]
步骤3.4：根据参数计算点簇中的点计算局部欧式距离，具体计算公式为：
[0082][0083]
其中，d
loc
为局部欧式距离，kr为参数计算点簇中近邻点数量，为待检测点在参考时相点云的第j个近邻点，(xi,yi,zi)为待检测点的三维坐标，(xj,yj,zj)为点的三维坐标。
[0084]
步骤3.5：计算参数计算点簇的中心，在参考时相点云中搜索kd个邻近点计算局部密度，具体计算公式为：
[0085][0086]
其中，ρ
loc
为局部密度，kd为参考时相点云中近邻点的数量，r
max
为参数计算点簇的中心到kd个邻近点的最大距离。
[0087]
步骤4：根据局部密度改善局部欧式距离，得到密度自适应局部欧式距离，具体计算公式为：
[0088][0089]
其中，为密度自适应局部欧式距离，d
loc
为局部欧式距离，ρ
loc
为局部密度。
[0090]
步骤5：将待检测时相点云中所有密度自适应局部欧式距离大于距离阈值td的待检测点提取到候选变化区域，密度自适应局部欧式距离小于等于距离阈值td的待检测点则认为未发生变化，不需要提取到候选变化区域中。
[0091]
本实施例中，距离阈值td设为1.2m。
[0092]
步骤6：搜索候选变化区域内候选变化点qi在待检测时相点云中xy平面上的k个邻近点，并对该k个邻近点进行欧式聚类，得到第二聚类结果，根据第二聚类结果对候选变化点qi进行移出或不移出候选变化区域的处理。
[0093]
本实施例中，根据第二聚类结果对候选变化点qi进行移出或不移出候选变化区域的处理的具体实现过程为：
[0094]
步骤6.1：当第二聚类结果中聚类簇的数量为1时，候选变化点qi不属于地物建模不好区域，不对候选变化点qi进行移出处理。
[0095]
步骤6.2：当第二聚类结果中聚类簇的数量为2、3或4时，候选变化点qi属于地物建模不好区域，若候选变化点qi的第一顶面点存在且未发生变化，则将候选变化点qi移出候选变化区域；否则将不对候选变化点qi进行处理；其中，候选变化点qi的第一顶面点是该k个邻近点中点的法向量与z方向之间的夹角θ小于角度阈值t
angle
且高程值最大的点。
[0096]
步骤6.3：当第二聚类结果中聚类簇的数量大于4时，候选变化点qi属于地物建模不好区域，若候选变化点qi的第二顶面点发生变化，则不对候选变化点qi进行处理；否则将候选变化点qi移出候选变化区域；其中，候选变化点qi的第二顶面点是该k个邻近点中高程值最大的点。
[0097]
步骤7：对经过步骤6处理后的候选变化区域内的点进行聚类，并根据每个聚类簇的平均密度计算该聚类簇对应的面积s；若面积s大于面积阈值t
area
，则认为该聚类簇真正发生变化；否则认为该聚类簇为误检测。
[0098]
本实施例中，每个聚类簇的平均密度计算公式为：
[0099][0100]
其中，为当前聚类簇的平均密度，kc为当前聚类簇中点的数量，(xj,yj,zj)为当前聚类簇中第j个点的三维坐标，(xn,yn,zn)为与当前聚类簇中第j个点对应的最邻近点的三维坐标。
[0101]
面积s的计算公式为：
[0102][0103]
其中，s为当前聚类簇的面积，为当前聚类簇的平均密度，kc为当前聚类簇中点的数量。本实施例中，面积阈值t
area
设为12平方米。
[0104]
以图2(b)为参考时相点云，图2(a)为待检测时相点云(第一种情况)，采用本发明方法进行三维变化检测得到的检测结果如图3(a)所示；以图2(a)为参考时相点云，图(b)为
待检测时相点云(第二种情况)，采用本发明方法进行三维变化检测得到的检测结果如图3(b)所示。经本发明所述方法进行检测后，检测结果中的每个待检测点均有一个标签(变或未变)，再将该检测结果与如图4所示人工标记的精度验证数据(人工标记每个待检测点的标签)进行比较，如果两者一致，则认为检测正确，否则认为检测错误，由此得到的检测精度如表1所示：
[0105]
表1本发明三维变化检测方法的检测精度
[0106]
类别精确度precision召回率recallf1-score第一种情况0.8470.8060.826第二种情况0.9160.8700.892
[0107]
采用传统的欧氏距离分割的三维变化检测方法进行检测，并与人工标记的精度验证数据比较，得到的检测精度如表2所示：
[0108]
表2传统的欧氏距离分割的三维变化检测方法的检测精度
[0109]
类别精确度precision召回率recallf1-score第一种情况0.3490.8420.493第二种情况0.8020.9430.867
[0110]
评价指标使用平均f1分数来评估检测性能，f1分数是每个点的精确度和召回率的调和平均值，公式定义如下：
[0111][0112]
其中，tp是本发明方法对待检测点检测的标签为变化，且人工对该点标记的标签也为变化的点的数量，tn是本发明方法对待检测点检测的标签为未变化，且人工对该点标记的标签也为未变化的点的数量，fp是本发明方法对待检测点检测的标签为未变化，但人工对该点标记的标签为变化的点的数量，fn是本发明方法对待检测点检测的标签为变化，但人工对该点标记的标签为未变化的点的数量。由表1和2可知，本发明的检测方法远优于传统检测方法，本发明具有更好的抗噪性能和更优的检测精度。
[0113]
以上所揭露的仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或变型，都应涵盖在本发明的保护范围之内。