CN106338277A - 一种基于基线的建筑物变化检测方法 - Google Patents

Abstract

本发明公开了一种基于基线的建筑物变化检测方法，其主要步骤如下：点云数据采集；结构坐标系建立；坐标转换；数据分割，数据去噪；特征点选取；虚拟点提取；基线建立；基线投影；变化检测获取变形信息。与现有技术相比，本发明自动化程度高，整个变化检测过程无需配准，解决了传统变化检测方法配准误差对结果的影响，这是传统基于影像的变化检测方法以及以往基于LiDAR的变化检测方法所不能达到的。

Description

一种基于基线的建筑物变化检测方法

技术领域

本发明属于测绘科学或者结构工程领域，涉及一种基于基线的建筑物变化检测方法。

背景技术

在城市中，建筑物是人类活动的主要场所，其安全状况关系到人类日常生活及经济活动，因此，对建筑物进行变化检测具有十分重要的意义。因此，基于建筑物的变化检测和修复，尤其是由地震引起的建筑物变形与损坏，是近年来业内的研究重点。在建筑物结构的发展史上，砖石结构自古以来常被广泛用于建筑物的基础结构。用砖、石块、砌体及土坯等各种块体，以灰浆(砂浆、黏土浆等)砌筑而成的一种组合体成为砖石结构。砖石结构因造价低、耐火性、耐久性及施工简易被广泛使用，但是其砌体强度较低，抗震能力较差。因此，对砖石结构的建筑物进行变化检测的研究具有及其重要的现实意义。

众多学者已经提出了多种基于影像的变化检测方法，按照信息处理的层次，可以分为基于像素、基于特征和基于目标的变化检测。近年来，对于激光雷达(Light DetectionAnd Ranging)LiDAR数据的变化检测也逐渐开始在摄影测量与计算机视觉界得到研究。而以上传统变化检测方法均需配准不同时期的数据(点云，图像等)，实际上，配准的过程可能会引入小的误差，而经验表明即使是很小的误差也会导致错误的变化检测结果或者不可靠的变化检测结果。

因此，基于上述问题，本发明提供了一种基于基线的建筑物变化检测方法，该方法首次解决了配准误差对变化检测结果的影响。

发明内容

为解决上述问题，本发明的基于基线的建筑物变化检测方法所采用的技术方案为：

一种基于基线的建筑物变化检测方法，包括以下步骤：

1)采用激光扫描仪系统对同一建筑物进行两期扫描，获取建筑物表面点云数据；观测值为建筑物表面点的三维坐标和激光反射强度；

2)根据扫描场景内稳定区域选取稳定墙面点云数据，利用主成分分析方法对其进行降维分析，确定结构坐标系的坐标轴指向；将测站坐标系原点投影至固定墙面，计算平移向量；从而建立特殊的结构坐标系；

3)根据步骤2)获得的结构坐标系参数，计算原测站坐标系与结构坐标系之间的转换参数，将扫描的点云数据转化至结构坐标系；

4)将步骤3)获得的点云数据进行分割和去噪；

5)在扫描场景内选取具有代表意义的特征点；

6)分别提取虚拟点坐标来表示特征点；

7)利用步骤6)选取的特征点建立基线；

8)根据步骤7)建立的基线，将基线向量投影至结构坐标系；

9)根据步骤8)获得的两期相对应基线在结构坐标系上的投影分量进行变化检测。

步骤2)所述特殊的结构坐标系为一种根据该建筑物及周边环境而建立的坐标系，其中主成分分析方法及结构坐标系建立过程如下：

设扫描点X的三维坐标{X_i＝(x_i,yi,z_i)|i＝1,2,…,n}，构造相应的协方差矩阵C：

$C=\frac{1}{n}{M}^{T}M;$

其中， 为点集的重心坐标；()^T表示矩阵转置；

对协方差矩阵C进行主成分分析，可求得三个特征值：λ₁、λ₂和λ₃，对它们按降序排列，得到λ₁≥λ₂>λ₃>0；λ₃所对应的特征向量v₃，v₃为法向量，即为结构坐标系的X轴在测站坐标系下的单位向量；而结构坐标系的Z轴指向与测站坐标系一致，Y轴垂直于确定的XOZ平面，构成右手坐标系；

计算测站点坐标S(0,0,0)在固定墙面所在平面的投影S'(x_s,y_s,z_s)，将其作为结构坐标系的坐标原点，因此平移向量(Δx,Δy,Δz)＝(-x_s,-y_s,-z_s)；

在确立了平移参数及坐标轴旋转参数后，将两期配准后的点云数据旋转至结构坐标系。

步骤3)所述坐标转换的过程如下：

设矩阵A为A坐标系下的点云三维坐标，矩阵B为B坐标系下的点云三维坐标，A、B两坐标系的三维坐标转换方程如下所示：

${\left[\begin{array}{c}x\\ y\\ z\end{array}\right]}_B=\left[\begin{array}{c}\mathrm{\Δ}x\\ \mathrm{\Δ}y\\ \mathrm{\Δ}z\end{array}\right]+(1+k)R{\left[\begin{array}{c}x\\ y\\ z\end{array}\right]}_A$

式中：Δx、Δy和Δz分别表示坐标原点在X方向、Y方向和Z方向上的平移量，k为尺度因子，k＝0，R为A坐标系到B坐标系的旋转矩阵。

步骤6)所述虚拟点坐标的提取过程如下：

步骤1，根据获取点云数据的强度信息，采用K均值聚类算法对变化墙面的点云进行分类，分离得到砖块点云和砂浆点云；

步骤2，利用步骤1得到的砂浆点云，将点云坐标投影至Z方向和Y方向，定义固定窗口长度L_fix和移动窗口长度L_move，通过移动窗口，分别计算沿Z方向和Y方向点云线密度变化；

步骤3，根据步骤2得到的线密度变化，分别求取各砖块间横向和纵向分割线，计算每个砖块四个角点坐标，建立砖块模型；

步骤4，根据步骤3得到的砖块模型获取各砖块点云，计算各砖块中心；

步骤5，根据步骤4得到的两期相应的砖块中心三维坐标。

上述步骤1中基于强度信息的K均值聚类方法分离墙面砖块点云和砂浆点云的具体方法为：

采用聚类误差平方和函数E作为聚类准则函数，以点强度信息作为分类属性，

$E=\underset{i=1}{\overset{k}{\Σ}}\underset{j=1}{\overset{n_i}{\Σ}}\left|\right|x_{ij}-m_i|{|}^2,$

其中，x_ij是第i类第j个样本，m_i是第i类的聚类中心，n_i是第i类样本个数。

K均值聚类算法的实质就是通过反复迭代寻找k个最佳的聚类中心，k＝2，将全体n个样本点分配到离它最近的聚类中心，使得聚类误差的平方和E最小；其实现过程如下：

步骤一，随机指定k个聚类中心m_i(i＝1,2,…,k)；

步骤二，对每一个样本x_i找到离它最近的聚类中心，将其分配到该类；

步骤三，重新计算各簇新中心：

其中，i＝1,2,…,k；N_i是第i簇当前样本数；

步骤四，计算偏差，

步骤五，如果E值收敛，则返回m_i(i＝1,2,…,k)，算法终止；否则，重新执行步骤二至步骤五。

步骤2所提利用砂浆点云，基于窗口移动法计算点云线密度变化，其具体方法及说明如下：

设砂浆的平均宽度为L_mortar，该值为已知量,，定义固定窗口长度L_fix和移动窗口长度L_move，三者满足如下关系：

L_mortar≈L_fix+2L_move；

分别计算沿Z方向和Y方向移动窗口数目：

$\begin{array}{cc}{n}_y=\[\frac{y_{max}-y_{\min }-L_{fix}}{L_{move}}\]+1& n_z=\[\frac{z_{max}-z_{\min }-L_{fix}}{L_{move}}\]+1\end{array};$

其中[]为取整符号；y_max和y_min分别表示Y坐标的最大值和最小值；z_max和z_min分别表示Z坐标的最大值和最小值；

分别计算沿Z方向和Y方向各个窗口内的点数目：

n_zi(i＝1,2,…,n_y),n_yi(i＝1,2,…,n_z)；

分别计算沿Z方向和Y方向的点的线密度：

Density_z＝(n_z(i-1)+n_zi+n_z(i+1))/(3L_fix)(i＝2,3,…,(n_z-1))；

Density_y＝(n_y(i-1)+n_yi+n_y(i+1))/(3L_fix)(i＝2,3,…,(n_y-1))；

对于Z方向和Y方向各个窗口，计算其线密度变化率：

Grad(i,1)＝Density_y(i)-Density_y(i-1)

Grad(i,2)＝Density_y(i+1)-Density_y(i)。

利用线密度变化，计算各砖块间横向和纵向分割线，其具体方法及说明如下：

对于分析的Z方向或Y方向，点云线密度在垂直于该方向的砂浆接缝处明显高于其他地方，因此，选定该方向上窗口平均密度作为阈值，其中，ε_y和ε_z分别表示Y和Z方向上的阈值，则n_total为点云总数目，当窗口内点云密度大于阈值，且满足Grad(i,1)>0和Grad(i,2)<0时，此时所处固定窗口即为某砂浆接缝的范围；

通过砂浆接缝范围内的点云计算砂浆接缝中心线；

通过砖块四周接缝中心线计算砖块模型四个角点坐标，建立砖块模型。

所述基线向量投影方法如下：

假设为测期一基线向量，为测期二基线向量，那么为基线两期内变化，ΔX₁、ΔY₁、ΔZ₁和ΔX₂、ΔY₂、ΔZ₂分别为测期一和测期二中的对应基线投影；θ为基线旋转角。

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

本发明所提供的一种基线的建筑物变化检测方法，具有两方面的优点：一方面自动化程度高；另一方面变化检测过程无需配准，解决了传统变化检测方法配准误差对结果的影响，这是传统基于影像的变化检测方法以及以往基于LiDAR的变化检测方法所不能达到的。

附图说明

图1为本发明的基于基线的建筑物变化检测方法的流程图；

图2为本发明的测站坐标系示意图；

图3为本发明的结构坐标系示意图；

图4为本发明的砖块模型示意图；

图5为本发明的砖块中心提取实施实例；

图6为本发明基线向量投影示意图；

图7为B区域内各砖块中心与A区域第一个砖块中心构成基线的变化矢量图；

图8为C区域内各砖块中心与A区域第一个砖块中心构成基线的变化矢量图；

图9为D区域内各砖块中心与A区域第一个砖块中心构成基线的变化矢量图。

具体实施方式

下面结合附图和具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图1所示方法流程，一种基于基线的建筑物变化检测方法，包括如下步骤：

1)采用激光扫描仪系统对同一建筑物进行两期扫描，获取建筑物表面点云数据；观测值为建筑物表面点的三维坐标和激光反射强度。

2)根据扫描场景内稳定区域选取稳定墙面点云数据，利用主成分分析方法对其进行降维分析，确定结构坐标系的坐标轴指向；将测站坐标系原点投影至固定墙面，计算平移向量；从而建立特殊的结构坐标系；

3)根据步骤2)获得的结构坐标系参数，计算原测站坐标系与结构坐标系之间的转换参数，将扫描的点云数据转化至结构坐标系；

4)将步骤3)获得的点云数据进行分割和去噪；

5)在扫描场景内选取具有代表意义的特征点；

6)分别提取虚拟点坐标来表示不同类型的特征点；

7)利用步骤6)选取的特征点建立基线；

8)根据步骤7)建立的基线，将基线向量投影至结构坐标系；

9)根据步骤8)获得的两期相对应基线在结构坐标系上的投影分量进行变化检测。

实施例1

如图1、图2、图3和图4所示，根据图1的方法流程，以“某实验场砖石结构建筑在地震测试前后变化检测”为应用实例，对本发明作进一步阐述：

(1)利用Leica C10激光扫描仪系统对建筑物进行扫描，分别在地震测试前后进行扫描，获取两期建筑物表面激光点云数据。观测值包含两类：三维坐标，激光反射强度；

(2)根据扫描场景内稳定区域选取稳定墙面点云数据，利用主成分分析方法对其进行降维分析，获取稳定墙面法向量即为结构坐标系的X轴单位向量，其中测期一为[0.9348 0.3552 0.0015]；测期二为[0.9974 -0.0727 -0.0002]；利用测站坐标系原点，将其投影至稳定墙面，计算平移向量，其中测期一为[-8.2760 -3.1451 -0.0132]，测期二为[-8.8284 0.6432 0.0020]；根据上述参数，两期结构坐标系分别建立完毕；

(3)根据上述结构坐标系与测站坐标系间的转换关系，将扫描的点云数据转换至结构坐标系；

(4)综合考虑变化墙面各部位分布，以求尽可能覆盖整个墙面，选取变化墙面左上(设为A区域)、左下(设为C区域)、右上(设为B区域)及右下(设为D区域)四个研究对象，利用Cloudcompare软件对其进行去噪；

(5)选取砖块中心为特征点；

(6)采用如上所示方法对各个砖块中心进行提取；以左上区域为例，提取结果如图5所示；

(7)连接不同区域的不同砖块中心，构成基线，两期对应点连成的基线为对应基线，用于变化检测；

(8)按如图6所示方法，将两期对应基线向量投影至结构坐标系；

(9)通过比较两期基线投影后各分量得到各基线变化，以B区域、C区域和D区域内各砖块中心与A区域第一个砖块中心构成基线的变化检测结果为例，绘制其变化矢量图如图7-9所示。

Claims (9)

1.一种基于基线的建筑物变化检测方法，其特征在于：包括以下步骤：

1)采用激光扫描仪系统对同一建筑物进行两期扫描，获取建筑物表面点云数据；观测值为建筑物表面点的三维坐标和激光反射强度；

2)根据扫描场景内稳定区域选取稳定墙面点云数据，利用主成分分析方法对其进行降维分析，确定结构坐标系的坐标轴指向；将测站坐标系原点投影至稳定墙面，计算平移向量；从而建立特殊的结构坐标系；

3)根据步骤2)获得的结构坐标系参数，计算原测站坐标系与结构坐标系之间的转换参数，将扫描的点云数据转化至结构坐标系；

4)将步骤3)获得的点云数据进行分割和去噪；

5)在扫描场景内选取具有代表意义的特征点；

6)分别提取虚拟点坐标来表示特征点；

7)利用步骤6)选取的特征点建立基线；

8)根据步骤7)建立的基线，将基线向量投影至结构坐标系；

9)根据步骤8)获得的两期相对应基线在结构坐标系上的投影分量进行变化检测。

2.如权利要求1所述的基于基线的建筑物变化检测方法，其特征在于：步骤2)所述结构坐标系为一种根据该建筑物及周边环境而建立的坐标系，其中主成分分析方法及结构坐标系建立过程如下：

设扫描点X的三维坐标{X_i＝(x_i,y_i,z_i)|i＝1,2,…,n}，构造相应的协方差矩阵C：

$C=\frac{1}{n}{M}^{T}M;$

其中， 为点集的重心坐标；()^T表示矩阵转置；

对协方差矩阵C进行主成分分析，可求得三个特征值：λ₁、λ₂和λ₃，对它们按降序排列，得到λ₁≥λ₂>λ₃>0；λ₃所对应的特征向量v₃，v₃为法向量，即为结构坐标系的X轴在测站坐标系下的单位向量；而结构坐标系的Z轴指向与测站坐标系一致，Y轴垂直于确定的XOZ平面，构成右手坐标系；

计算测站点坐标S(0,0,0)在固定墙面所在平面的投影S'(x_s,y_s,z_s)，将其作为结构坐标系的坐标原点，因此平移向量(Δx,Δy,Δz)＝(-x_s,-y_s,-z_s)；

在确立了平移参数及坐标轴旋转参数后，将两期配准后的点云数据转换至结构坐标系。

3.如权利要求1所述的基于基线的建筑物变化检测方法，其特征在于：步骤3)所述坐标转换的过程如下：

设矩阵A为A坐标系下的点云三维坐标，矩阵B为B坐标系下的点云三维坐标，A、B两坐标系的三维坐标转换方程如下所示：

${\left[\begin{array}{c}x\\ y\\ z\end{array}\right]}_B=\left[\begin{array}{c}\mathrm{\&Delta;}x\\ \mathrm{\&Delta;}y\\ \mathrm{\&Delta;}z\end{array}\right]+(1+k)R{\left[\begin{array}{c}x\\ y\\ z\end{array}\right]}_A;$

式中：Δx、Δy和Δz分别表示坐标原点在X方向、Y方向和Z方向上的平移量，k为尺度因子，k＝0，R为A坐标系到B坐标系的旋转矩阵。

4.如权利要求1所述的基于基线的建筑物变化检测方法，其特征在于：步骤6)所述虚拟点坐标的提取过程如下：

步骤1，根据获取点云数据的强度信息，采用K均值聚类算法对变化墙面的点云进行分类，分离得到砖块点云和砂浆点云；

步骤2，利用步骤1得到的砂浆点云，将点云坐标投影至Z方向和Y方向，定义固定窗口长度L_fix和移动窗口长度L_move，通过移动窗口，分别计算沿Z方向和Y方向点云线密度变化；

步骤3，根据步骤2得到的线密度变化，分别求取各砖块间横向和纵向分割线，计算每个砖块四个角点坐标，建立砖块模型；

步骤4，根据步骤3得到的砖块模型获取各砖块点云，计算各砖块中心；

步骤5，根据步骤4得到的两期相应的砖块中心三维坐标。

5.如权利要求4所述的基于基线的建筑物变化检测方法，其特征在于：步骤1中基于强度信息的K均值聚类方法分离墙面砖块点云和砂浆点云的具体方法为：

采用聚类误差平方和函数E作为聚类准则函数，以点强度信息作为分类属性，

其中，x_ij是第i类第j个样本，m_i是第i类的聚类中心，n_i是第i类样本个数。

6.如权利要求4所述的基于基线的建筑物变化检测方法，其特征在于：K均值聚类算法的实质就是通过反复迭代寻找k个最佳的聚类中心，k＝2，将全体n个样本点分配到离它最近的聚类中心，使得聚类误差的平方和E最小；其实现过程如下：

步骤一，随机指定k个聚类中心m_i(i＝1,2,…,k)；

步骤二，对每一个样本x_i找到离它最近的聚类中心，将其分配到该类；

步骤三，重新计算各簇新中心：

其中，i＝1,2,…,k，N_i是第i簇当前样本数；

步骤四，计算偏差，

步骤五，如果E值收敛，则返回m_i(i＝1,2,…,k)，算法终止；否则，重新执行步骤二至步骤五。

7.如权利要求4所述的基于基线的建筑物变化检测方法，其特征在于：步骤2所提利用砂浆点云，基于窗口移动法计算点云线密度变化，其具体方法如下：

设砂浆的平均宽度为L_mortar，该值为已知量，定义固定窗口长度L_fix和移动窗口长度L_move，三者满足如下关系：

L_mortar≈L_fix+2L_move；

分别计算沿Z方向和Y方向移动窗口数目：

$\begin{array}{cc}{n}_y=\&lsqb;\frac{y_{max}-y_{\min }-L_{fix}}{L_{move}}\&rsqb;+1& n_z=\&lsqb;\frac{z_{max}-z_{\min }-L_{fix}}{L_{move}}\&rsqb;+1\end{array};$

其中，[]为取整符号；y_max和y_min分别表示Y坐标的最大值和最小值；

z_max和z_min分别表示Z坐标的最大值和最小值；

分别计算沿Z方向和Y方向各个窗口内的点数目：

n_zi(i＝1,2,…,n_y),n_yi(i＝1,2,…,n_z)；

分别计算沿Z方向和Y方向的点的线密度：

Density_z＝(n_z(i-1)+n_zi+n_z(i+1))/(3L_fix)(i＝2,3,…,(n_z-1))；

Density_y＝(n_y(i-1)+n_yi+n_y(i+1))/(3L_fix)(i＝2,3,…,(n_y-1))；

对于Z方向和Y方向各个窗口，计算其线密度变化率：

Grad(i,1)＝Density_y(i)-Density_y(i-1)

Grad(i,2)＝Density_y(i+1)-Density_y(i)。

8.如权利要求7所述的基于基线的建筑物变化检测方法，其特征在于：所述利用线密度变化，计算各砖块间横向和纵向分割线，其具体方法及说明如下：

对于分析的Z方向或Y方向，点云线密度在垂直于该方向的砂浆接缝处明显高于其他地方，因此，选定该方向上窗口平均密度作为阈值，其中，ε_y和ε_z分别表示Y和Z方向上的阈值，则n_total为点云总数目，当窗口内点云密度大于阈值，且满足Grad(i,1)>0和Grad(i,2)<0时，此时所处固定窗口即为某砂浆接缝的范围；

通过砂浆接缝范围内的点云计算砂浆接缝中心线；

通过砖块四周接缝中心线计算砖块模型四个角点坐标，建立砖块模型。

9.如权利要求1所述的基于基线的建筑物变化检测方法，其特征在于：步骤8)基线向量投影方法如下：

因三维向量是一个矢量，其可以看成以坐标原点为起点，带方向的线段，因此可以将基线一端平移至坐标原点；假设为测期一基线向量，为测期二基线向量，那么为基线两期内变化；ΔX₁、ΔY₁、ΔZ₁和ΔX₂、ΔY₂、ΔZ₂分别为测期一和测期二中的对应基线投影；θ为基线旋转角。