一种基于点云的砖混建筑结构长度计算方法

1.本发明属于建筑结构长度计算技术领域，尤其涉及一种基于点云的砖混建筑结构长度计算方法。


背景技术：

2.砖混结构多用于乡村矮层建筑，该结构建筑几何结构较为简单，除门窗外，其墙体不含有用于装饰的附属结构，造价较为经济。但此类建筑物支持的使用年限不长，易受地表岩层移动的影响。且乡村为矿产资源开采的主要地区，该结构特点及采矿影响致使砖混建筑经常遭到破坏。因此，进行该类建筑的变形监测，以为建筑的健康状况及损坏等级评估提供依据，具有十分重要的意义。
3.针对砖混建筑的变形监测，传统监测方法以全站仪为主要设备，但该方法费时费力，且仅能获取少数离散点。随着测绘技术的发展，逐渐形成了以差分干涉测量(dinsar)、无人机摄影测量(uav photogrammetry)、和地面激光扫描(tls)技术为代表的新型监测方法。相比于dinsar、uav photogrammetry，tls可以高精度、海量点云的方式来表示物体表面信息，被广泛应用。
4.近年来，通过tls进行建筑物变形监测，可分为两类：a.结合建立的模型进行监测，b.基于点云直接进行监测。在a类监测中，juan et al.建立了建筑物信息模型(bim)，通过两个时间段的模型对比，计算了该期间的位移变化。raphael et al.同样建立了bim，最终计算了施加重力下的位移与应力变化。lian et al.建立了有限元模型(fem)，在煤层开采的条件下，计算了墙体主应力与剪切力变化。在b类监测中，batur et al.对点云进行配准、滤波等处理后，通过cloud-to-cloud的方式进行了位移计算。suchocki et al.将点云划分为条带状，通过optimum dataset方法减少数据量后，进行了墙体裂缝的检测。li et al.基于点云提取离散变形后，结合沉陷预计模型，计算了墙体下沉、水平移动、倾斜等变形。yang et al.在点云的基础上，结合alpha shape算法创建了三角网格向量模型，最终获取了墙体倾斜、裂缝深度。基于上述可以发现，a类监测在建立模型、施加力、对比模型时，需要借助多种商业软件。而b类监测是在点云层次上进行的，其监测方法及监测指标的选择更加灵活，此类监测似乎是一个发展趋势。
5.在中国进行建筑物监测，指标有下沉、水平移动、倾斜、曲率、水平变形等。其中，下沉、水平移动带来的影响主要为建筑物所在潜水位的变化，并不会造成结构破坏。在实际应用中，由于倾斜、曲率值较小，水平变形是造成结构破坏的主要因素。该值的计算通常需要在建筑物上布设人工标志。对于因煤矿开采引起的破坏，也可结合开采条件进行计算。但这两种计算方式获取水平变形均是少量的，且需将多期观测数据配准至同一坐标系下。增加了人力物力，且引入了误差。当水平变形造成的结构破坏达到一定程度时，最直接的现象在于墙体出现裂缝，反应在几何结构上，即为墙体结构长度发生了变化。因此，本发明直接选择结构长度做为监测指标，且考虑到基于点云直接进行监测的优势、及砖混建筑物几何结构特点，提出了基于点云的砖混建筑结构长度计算方法。


技术实现要素：

6.本发明的目的就在于为了解决上述问题而提供一种基于点云的砖混建筑结构长度计算方法。
7.本发明通过以下技术方案来实现上述目的：
8.一种基于点云的砖混建筑结构长度计算方法，包括以下步骤：
9.s1：获取建筑物墙体点云，对所述建筑物墙体点云进行预处理，获得处理后建筑物墙体点云；
10.s2：在所述处理后建筑物墙体点云内设置多个邻域，并对多个所述邻域内的点云进行过滤，得到有效邻域点，并计算出墙体点云厚度值；
11.s3：基于所述墙体点云厚度值确定第一距离阈值和第二距离阈值，基于所述第一距离阈值对所述处理后建筑物墙体点云进行建筑物点云分割，得到每个单面墙体的组成点云，基于所述第二距离阈值对每个所述单面墙体的组成点云进行墙体点云再分配，得到分配后墙体点云；
12.s4：采用改进后的基于法线的点云边界估计方法对所述分配后墙体点云进行墙体边界点估计，得到墙体边界点；
13.s5：对所述墙体边界点进行处理得到上边界线，在交点约束条件下采用选权迭代最小二乘法对所述上边界线进行拟合，得到拟合后上边界线；
14.s6：基于所述拟合后上边界线，在墙体上创建兴趣区域，根据所述兴趣区域内垂直边界点所在直线，计算得到墙体结构长度。
15.作为本发明的进一步优化方案，所述步骤s1中获取建筑物墙体点云，对所述建筑物墙体点云进行预处理，获得处理后建筑物墙体点云，具体过程如下：
16.对建筑物进行扫描，得到建筑物墙体点云，对所述建筑物墙体点云进行降噪、抽稀、地面滤波处理，得到处理后建筑物墙体点云。
17.作为本发明的进一步优化方案，所述步骤s2中在所述处理后建筑物墙体点云内设置多个邻域，并对多个所述邻域内的点云进行过滤，得到有效邻域点，并计算出墙体点云厚度值，具体过程如下：
18.在所述处理后建筑物墙体点云内随机选取n个点，每个点设定以r为半径划定邻域；
19.计算邻域内的点到该邻域所在平面的距离d，将所述距离d小于距离中误差的点作为有效邻域点，计算出每个邻域的点云厚度th
nei
，如下式所示：
[0020][0021]
其中，rmse1表示单个邻域的距离中误差；n、m分别为单个邻域中的原始点的数量和有效邻域点的数量；[d1,d2,
…
,di,
…
,dn]和[d1,d2,
…
,dj,
…
,dm]分别表示单个邻域中的原始点至该邻域所在平面距离的集合和单个邻域中的有效邻域点至该邻域所在平面距离的集合；d
mean
表示集合[d1,d2,
…
,dj,
…
,dm]中有效邻域点至邻域所在平面距离的均值；
[0022]
将每个邻域的所述点云厚度th
nei
进行均值计算，得到邻域点云厚度均值，取小于厚度中误差的邻域为有效邻域，计算得到所述墙体点云厚度值th
wall
，如下式所示：
[0023][0024]
其中，rmse2表示邻域点云厚度中误差，n、m分别为原始的邻域数量和有效邻域的数量；和和分别为原始的邻域的点云厚度集合和有效邻域的点云厚度集合；为集合中的有效邻域的点云厚度均值。
[0025]
作为本发明的进一步优化方案，所述步骤s3中基于所述墙体点云厚度值确定第一距离阈值和第二距离阈值，基于所述第一距离阈值对所述处理后建筑物墙体点云进行建筑物点云分割，得到每个单面墙体的组成点云，基于所述第二距离阈值对每个所述单面墙体的组成点云进行墙体点云再分配，得到分配后墙体点云，具体过程如下：
[0026]
基于所述墙体点云厚度值确定第一距离阈值和第二距离阈值，所述第一距离阈值大于所述墙体点云厚度值，所述第二距离阈值等于所述墙体点云厚度值；
[0027]
将所述处理后建筑物墙体点云投影至xoy面，基于第一距离阈值和随机抽样一致性算法对投影至xoy面的点云进行处理，得到单面墙体的组成点云；从所述处理后建筑物墙体点云中删除所述单面墙体的组成点云，完成数据集合的一次更新；重复随机抽样一致性算法处理步骤和数据集合更新步骤，直至获取每个单面墙体的组成点云；
[0028]
采用随机抽样一致性算法对每个所述单面墙体的组成点云进行处理，获得该单面墙体点云；对单面墙体的组成点云和单面墙体点云求差集，得到差集点云；计算差集点云中每个点至每个单面墙体的距离；取所述差集点云中的每一个到单面墙体的距离小于第二距离阈值的点，分配至该单面墙体的点云中，得到每个单面墙体的分配后墙体点云。
[0029]
作为本发明的进一步优化方案，所述步骤s4中采用改进后的基于法线的点云边界估计方法对所述分配后墙体点云进行墙体边界点估计，得到墙体边界点，具体过程如下：
[0030]
采用选权迭代最小二乘法对所述分配后墙体点云进行计算，得到平面的模型参数；基于平面的模型参数分别计算出分配后墙体点云与xoz面、xoy面间的旋转角，如下式所示：
[0031]
θ＝xaxb/|xa||xb|；
[0032]
其中，θ为旋转角；| |为模型参数的模长，其中xa、xb分别为预设的已知的平面a和平面b的模型参数；
[0033]
基于旋转角将分配后墙体点云旋转至与xoy面平行，得到旋转后的墙体点云；
[0034]
在旋转后的墙体点云中选择一点，将该点与其邻域点连线，计算相邻连线间的夹角，并按照角度大小排列；
[0035]
计算排列后的相邻的夹角之间的差值，以差值大于预设的角度阈值的点为墙体边界点。
[0036]
作为本发明的进一步优化方案，所述步骤s5中对所述墙体边界点进行处理得到上边界线，在交点约束条件下采用选权迭代最小二乘法对所述上边界线进行拟合，得到拟合后上边界线，具体过程如下：
[0037]
将所述墙体边界点分为内边界点和外边界点，得到外边界点所在直线的边界线，包括垂直边界线和水平边界线，将与屋顶毗邻的水平边界线作为上边界线；
[0038]
将所述上边界线垂直投影至xoy面，并采用选权迭代最小二乘法求解投影至xoy面上的上边界线的模型参数；
[0039]
基于xoy面上的上边界线的模型参数分别计算出上边界线所在墙体平面与xoz面、yoz面间的旋转角；
[0040]
将上边界线水平投影至旋转角较小的xoz面或yoz面上，采用选权迭代最小二乘法求解上边界线投影后的直线的模型参数，上边界线在两个面上的线a和线b分别由下列两个直线方程式表示：
[0041]
k1x+b1＝y＝(z-b2)/k2；
[0042]
k3x+b3＝y＝(z-b4)/k4；
[0043]
其中，k1、k2、b1、b2是线a在空间坐标系中的方程表达式参数，k3、k4、b3、b4是线b在空间坐标系中的方程表达式参数；
[0044]
对两个直线方程进行求解得到交点，设点a、点b为求解得到的交点，则拟合后上边界线如下式所示：
[0045][0046]
其中，(xa,ya,za)、(xb,yb,zb)分别为边界线交点a和点b的坐标；
[0047]
作为本发明的进一步优化方案，所述步骤s6中基于所述拟合后上边界线，在墙体上创建兴趣区域，根据所述兴趣区域内垂直边界点所在直线，计算得到墙体结构长度，具体过程如下：
[0048]
设定步长为h，将拟合后上边界线依次沿墙体平面向下平移，并以每次平移后的拟合后上边界线为中心线，将高度为2h的区域划定为兴趣区域；
[0049]
筛选出所述兴趣区域内位于该墙体本身的边界点和位于相邻墙体的边界点；
[0050]
将位于墙体本身的边界点和位于相邻墙体的边界点合并，并采用选权迭代最小二乘法计算出兴趣区域中的垂直边界线的模型参数，如下式所示：
[0051]
k5x+b5＝k6y+b6＝z；
[0052]
其中，k5、k6、b5、b6是垂直边界线的模型参数在空间坐标系中的方程表达式参数；
[0053]
设定平移次数为i，第i次平移后的拟合后上边界线如下式所示：
[0054]
(x-f)/g＝y＝(x-h-hi)/l；
[0055]
基于垂直边界线的模型参数和第i次平移后的拟合后上边界线，设定点a和点b是兴趣区域内两交点，则墙体长度为：
[0056][0057]
其中，l为当前平移位置的墙体结构长度；(xa,ya,za)、(xb,yb,zb)分别为兴趣区域内交点a和点b的坐标。
[0058]
本发明的有益效果在于：
[0059]
本发明通过选择结构长度做为监测指标，结合点云监测的优势和砖混建筑物的几何结构特点，可以十分精确地计算出砖混建筑结构长度，基于结构长度可更好的监测砖混建筑的变形情况。
附图说明
[0060]
图1是本发明的方法流程图；
[0061]
图2是本发明的建筑物及单个墙面点云图；
[0062]
图3是本发明的墙体边界点估计示意图；
[0063]
图4是本发明的边界线模型参数求解图；
[0064]
图5是本发明的墙体结构长度计算示意图；
[0065]
图6是本发明的实施例概况图；
[0066]
图7是本发明的变形监测流程图；
[0067]
图8是本发明的实施例中墙体点云厚度计算结果示意图；
[0068]
图9是本发明的实施例的流程中间结果展示图；
[0069]
图10是本发明的实施例中的结构长度变形情况图；
[0070]
图11是本发明的实施例中测量现场实拍图；
[0071]
图12是本发明的实施例中边界估计结果示意图；
[0072]
图13是本发明的边界点估计方法造成的结构长度误差情况示意图；
[0073]
图14是本发明的选权函数造成的结构长度误差情况示意图。
具体实施方式
[0074]
下面结合附图对本技术作进一步详细描述，有必要在此指出的是，以下具体实施方式只用于对本技术进行进一步的说明，不能理解为对本技术保护范围的限制，该领域的技术人员可以根据上述申请内容对本技术作出一些非本质的改进和调整。
[0075]
实施例1
[0076]
如图1-14所示，一种基于点云的砖混建筑结构长度计算方法，包括以下步骤：
[0077]
墙体点云厚度计算，计算点云厚度时，随机选择多个邻域，但仅保留部分邻域点、部分邻域，以达到稳健计算的目的；
[0078]
建筑物墙体拆分，基于墙体点云厚度，在确定的2个距离阈值后，分2步：建筑点云分割、墙体点云再分配实现了建筑物墙体拆分；
[0079]
墙体边界点估计，以单面墙体为基本单位，进行了墙体边界点估计，该估计方法由基于法线的点云边界估计方法变化得到，仅需计算一次法线，转换一次坐标系；
[0080]
上边界拟合，基于上边界点，将选权迭代最小二乘法引入，同时在相邻两上边界线相交于一点的约束下，进行该线拟合；
[0081]
结构长度计算，由上边界线，在墙体上创建了一可向下平移的兴趣区域，每次平移后，结合区域内垂直边界点所在直线，计算墙体结构长度。
[0082]
本发明在计算结构长度中多次采用了随机抽样一致性算法(ransac)、选权迭代最小二乘法(wils)，因此，将其作为基本原理进行描述。
[0083]
随机抽样一致性算法
[0084]
在数据集合可明确用某一模型描述，且不存在粗差值、异常值的情况下。最小二乘法(ls)求解模型参数可通过公式(1)描述。
[0085]
x＝(b
t
pb)-1bt
pl(1)
[0086]
其中，x为模型参数求解值矩阵；b为观测方程系数矩阵；l为观测值矩阵；p为观测
值权阵(pij＝0,i≠j)。由于不存在粗差值、异常值，p取单位矩阵。
[0087]
现实情况下，数据集合难以避免异常值的存在，且很难客观有效的进行正常值的分离。当确定描述数据的模型后，ransac可解决这一问题。进行分离时，其流程如下：
[0088]
step1，在数据集合中随机选择n个值，由公式(1)求解模型参数。
[0089]
step2，计算数据值与模型间的几何距离，并统计小于距离阈值的数量。
[0090]
step3，循环上述参数求解、数量统计。达到循环次数后，最大数量对应的小于距离阈值的各数据值，即为正常值。
[0091]
此外，ransac亦可进行模型参数求解，此时应以最大数量对应的模型参数为最终结果。
[0092]
选权迭代最小二乘法
[0093]
同样，数据集合也有部分粗差值的存在。且在求解模型参数时，ls、ransac均要受其影响。而wils可通过选权函数，改变粗差值在平差中的权重，以达到减弱粗差影响的目的。
[0094]
以igg选权函数为例，其表达式如公式(3)所示。通过wils求解模型参数，主要流程如下：
[0095]
step1，基于数据集合整体，由公式(1)求解模型参数。
[0096]
step2，由公式(2)计算各数据值求解残差。并根据公式(3)更新观测值权阵。
[0097]
step3，循环上述模型求解、权阵更新。达到循环次数后，以最终求解模型参数为结果。
[0098]
v＝bx-l(2)
[0099]
其中，v为各数据值的求解残差矩阵。
[0100][0101]
其中，vi为第i个数据值的求解残差；s1、s2是保证稳健性的阈值系数；σ0为求解残差的标准偏差。在这里，s1取1.0，s2取2.5。
[0102]
墙体点云厚度计算
[0103]
当点云可用平面模型描述时，点云厚度可通过点至平面距离表示。基于第一次观测的建筑物点云，随机选择了n个点，为每个点划定了以r为半径的邻域，用于计算墙体点云厚度。
[0104]
需要注意的是，上述每个邻域的点，并不全部参与至点云厚度计算。对于每个邻域，邻域点至邻域所在平面的距离被计算。考虑到点云噪声的存在，仅保留其距离小于距离中误差的点。在本实施例中，计算邻域内的点到该邻域所在平面的距离d，将所述距离d小于距离中误差的点作为有效邻域点，计算出每个邻域的点云厚度th
nei
，如下式所示：
[0105][0106]
其中，rmse1表示单个邻域的距离中误差；n、m分别为单个邻域中的原始点的数量和有效邻域点的数量；[d1,d2,
…
,di,
…
,dn]和[d1,d2,
…
,dj,
…
,dm]分别表示单个邻域中的
原始点至该邻域所在平面距离的集合和单个邻域中的有效邻域点至该邻域所在平面距离的集合；d
mean
表示集合[d1,d2,
…
,dj,
…
,dm]中有效邻域点至邻域所在平面距离的均值；
[0107]
将每个邻域的所述点云厚度th
nei
进行均值计算，得到邻域点云厚度均值，取小于厚度中误差的邻域为有效邻域，计算得到所述墙体点云厚度值th
wall
，如下式所示：
[0108][0109]
其中，rmse2表示邻域点云厚度中误差，n、m分别为原始的邻域数量和有效邻域的数量；和和分别为原始的邻域的点云厚度集合和有效邻域的点云厚度集合；为集合中的有效邻域的点云厚度均值。
[0110]
建筑物墙体拆分
[0111]
墙体是组成建筑物的基本单元，也是计算结构长度的基本单元。除建筑物墙体点云外，扫描仪同样会获取其扫描场内其他物体点云。并且，由于点云的无序性，各墙体点云无规则排列，共同组成建筑物点云。基于上述，建筑物墙体拆分分为建筑物点云分割何墙体点云再分配两个步骤：
[0112]
建筑物点云分割
[0113]
墙体在三维空间中近似垂直于xoy平面。以图2中四面墙体组成的建筑物为例，在此平面中，单面墙体的几何特征表现为直线形式，且建筑物的几何特征表现为多条直线闭合连接。基于这一特征，ransac被引入建筑物点云分割。
[0114]
具体来说，第一步，建筑物点云被投影至xoy面，以便建筑物点云被分割为墙体单元。第二步，进行ransac，最大数量对应的各点，即为单面墙体的组成点云。第三步，从建筑物点云中删除单面墙体的组成点云，完成数据集合的一次更新。第四步，重复ransac、集合更新，直至获取每一单面墙体的组成点云。
[0115]
其中，b、x、l构建如图2所示。为防止单面墙体投影在xoy面上的直线被分割为多条，距离阈值应大于墙体点云厚度。
[0116]
墙体点云再分配
[0117]
以大于点云厚度值为距离阈值，可避免单面墙体投影在xoy面上的直线被分割为多条。但同时，会造成建筑物点云的过度分割。
[0118]
按分割顺序，过度分割可分为三类：
[0119]
a.先被分割出的单面墙体，其组成点云包括该单面墙体的点云，但两侧存在相邻墙体点云。
[0120]
b.后被分割出的单面墙体，其组成点云仅包括该单面墙体点云，但其两侧点云缺失。
[0121]
c.一些墙体，其组成点包括该单面墙体点云，但一侧存在相邻墙体点云，另一侧点云缺失。因此，对于分割后的单面墙体组成点云，需进行再分配。
[0122]
第一步，针对每一单面墙体的组成点云，采用ransac，可获得该单面墙体点云。第二步，对单面墙体的组成点云、点云求差集。第三步，计算差集点云中，每一点至每一单面墙
体的距离。第四步，距离哪一单面墙体的距离小于距离阈值，则该点则被分配至该墙体的点云。
[0123]
其中，b、x、l构建如图3所示。此时，距离阈值应等于墙体点云厚度，并且同一点可被分配至多个单面墙体。
[0124]
墙体边界点估计
[0125]
目前，基于法线的点云边界估计方法已被广泛应用。在该方法中，首先需进行的是，以墙体点云中的点pi与邻域点构成的平面法线，建立三维坐标系。之后，在该坐标系的xoy面内，根据点pi与邻域点的几何关系，进行边界点与非边界点的判断。在该方法中，计算法线的邻域通常大于判断边界点的邻域，且邻域越大，边界估计效果越好，估计速度越慢。
[0126]
三维空间中，在两平面模型参数已知的情况下，可旋转一平面与另一平面平行。假定平面a、b的模型参数分别为xa、xb，则旋转角可如公式(4)表示。此时，可基于旋转角，将平面a旋转至于平面b平行。
[0127]
θ＝xaxb/|xa||xb|
ꢀꢀꢀ
(4)
[0128]
其中，θ为旋转角；||为模型参数的模长。
[0129]
考虑到单面墙体点云中所有点近似位于同一平面，且为减少给定参数，提高估计速度，本文对该估计方法进行了一些变化。主要表现为，将墙体点云转转至与xoy面平行后，进行边界估计。该旋转过程可分为两步：ⅰ.绕z轴旋转墙体点云与xoz面平行，ⅱ.绕x轴旋转墙体点云与xoy面平行，该过程可由公式(5)描述。
[0130][0131]
其中，(x,y,z)、(x’,y’,z’)分别为旋转前后的点坐标；θ1、θ2分别ⅰ、ⅱ步的旋转角，即与xoz面、xoy面的夹角。
[0132]
对于拆分后的单面墙体，本方法进行边界点估计的流程可如下所述：
[0133]
第一步，由wils求解平面模型参数。第二步，由公式(4)计算墙体点云与xoz、xoy面间的旋转角。第三步，由公式(5)旋转墙体点云与xoy面平行。第四步，如图3所示，在旋转后的墙体点云中，将点pi与其邻域点连线。第五步，计算相邻连接线间的夹角，并按照大小排列。第六步，计算相邻夹角差值，以差值大于角度阈值的点pi为边界点。第七步，对公式(5)进行逆变换，将墙体点云恢复至旋转前的姿态。
[0134]
上边界线拟合
[0135]
墙体边界点可分可分为两类：内边界点、外边界点。内边界点由门窗存在或点云空洞而产生，不参与至后续计算。边界线即为外边界点所在直线，亦可分为两类：垂直边界线，水平边界线。与屋顶毗邻的水平边界线，即上边界线，是进行建筑物结构长度计算的基础。
[0136]
如图4所示，观察上边界线几何特征可以发现，相邻两面墙体的上边界线具有同一交点。同时，上述估计的边界点中，粗差值必定存在，而ls、ransac均难以避免粗差值的影响。因此，本方法将wils引入，并在交点的约束下，进行上边界线拟合。
[0137]
第一步，上边界线被垂直投影至xoy面，并由wils求解直线的模型参数。第二步，由公式(4)计算上边界线所在墙体平面与xoz面、yoz面间的旋转角，即夹角。第三步，将上边界线水平投影至夹角较小的xoz面或yoz面后，同样由wils求解模型参数。其中，b、x、l构建如
图4所示。此时，上边界线a、b可分别由公式(6),公式(7)描述。在两直线方程已知的情况下，可求解其交点。假设点a、b是上边界线两交点，则该上边界线可由公式(8)描述。
[0138]
k1x+b1＝y＝(z-b2)/k2(6)
[0139]
k3x+b3＝y＝(z-b4)/k4(7)
[0140][0141]
其中，(xa,ya,za)、(xb,yb,zb)分别为边界线交点a、b坐标；
[0142]
g＝(xa-xb)/(ya-yb)为左侧式子分母；
[0143]
g＝xa+(xa-xb)/(ya-yb)为左侧式子分子符号之后的部分；
[0144]
h/l同理。
[0145]
结构长度计算
[0146]
为计算墙体结构长度，本方法将上边界线以h为步长，依次沿墙体平面向下平移。并以每次平移后的上边界线为中心线，定义以高度为2h的兴趣区域。
[0147]
如图5所示，在每个兴趣区域中，第一步，位于区域中的边界点被截取。由于进行边界点估计时，是以单面墙体为单位进行的，导致位于同一垂直边界上的边界点被分为两部分，分别位于墙体本身及其相邻墙体的边界点中。第二步，将上述被分为两部分的边界点合并，并由wils估计兴趣区域中的垂直边界线模型参数。其中，b、x、l构建如图4所示，该垂直边界线可由公式(9)描述。此时，第i次平移后的上边界线a可由公式(10)描述。假设点a和点b是某一兴趣区域内两交点，则墙体长度可由公式(11)计算。
[0148]
k5x+b5＝k6y+b6＝z
ꢀꢀꢀ
(9)
[0149]
(x-f)/g＝y＝(x-h-hi)/l
ꢀꢀꢀ
(10)
[0150][0151]
其中，l为当前平移位置的墙体结构长度；(xa,ya,za)、(xb,yb,zb)分别为兴趣区域内，点a和点b坐标。
[0152]
本方法研究案例所在地区为中国亳州市，在监测期间，共进行了两次观测，时间间隔约为11月，扫描仪参数及测量时间如图6(a)所示。对图6(b)中的点云进行降噪、抽稀、地面滤波等预处理后，根据图7，进行了该建筑物的结构长度计算，并进行变形分析。
[0153]
在本方法中，关于墙体点云厚度，随机选择了n＝40。对于r的取值，根据试错法的思想，首先在5-20cm的范围中，r以5cm为间隔依次增加。由图8(a)-(d)所示，各邻域、墙体点云厚度值随着r的增加而增加。但r达到15cm后，该变化并不明显。因此以10cm为间隔，使r继续增加了2次。此时，观察图8(a)-(f)可以发现，当r达到15cm后，被保留的邻域点云厚度没有明显变化，大部分不超过3cm。且由邻域点云厚度计算的墙体点云厚度同样没有明显变化，约为1.7cm。因此，本文以n＝40，r＝15cm计算的墙体点云厚度1.7cm作为确定距离阈值的基准。
[0154]
结构长度及变形结果分析
[0155]
在建筑物墙体拆分中，建筑物点云分割、墙体点云再分配中的距离阈值应分别大于、等于墙体点云厚度。基于上述中的计算结果，确定距离阈值分别为4.0cm、1.7cm，得到墙体拆分结果如图9a所示。由于墙体d附近有树木遮挡的缘故，采集获取的该墙体点云存在斑
块状漏洞，从而导致拆分后墙体d点云有一定的缺失。但相邻墙体已被完全拆分，且结构完整。之后，对于每一单面墙体点云，进行了边界点估计。在图9b所示的估计结果中，各墙体边界点组成的建筑物轮廓清晰。同样，因观测状况的影响，墙体d的边界点出现少许噪点，但仍可清晰区分内、外边界点，并不会对后续步骤造成影响。然后，基于上边界点，进行了上边界线拟合。为直观显示拟合结果，计算了拟合直线交点。观察图9c可以发现，上边界线交点与垂直边界点完全位于同一直线：垂直边界上，可反应拟合直线与实际情况相符合。在进行长度结构计算时，以h＝0.25m为兴趣区域高度，上边界线从上至下共平移20次，在每个兴趣区域中计算的交点及结构长度分别如图9d和表1所示。
[0156]
根据公式(15)，对相同平移次数的结构长度做差后，即可获取该平移位置的建筑物变形信息。为直观显示变形，其结果被绘制成图，如图10所示。
[0157]
le＝l
2-l1(15)
[0158]
其中，le为结构长度变形，l1为两次连续观测中，第一次观测时某平移位置的结构长度，l2为第二次观测时，相同平移位置的结构长度。
[0159]
表1结构长度计算结果
[0160][0161]
[0162]
观察表1和图10可以发现。对于每一墙体，两次观测间的结构长度的曲线形状具有较高的相似性。墙体a各高度处的变形分布在34-48mm。观测其变形趋势可以发现，以下移次数10为中心，即距屋顶2.5m处，变形呈越远离该处越大的趋势。结合现场实拍图图11a来看，墙体a的裂缝宽度从屋顶至地面呈逐渐减小趋势，但裂缝条数由1条变为2条。裂缝反应的结构长度变化与本文计算结果基本一致。墙体b各高度处的变形较小，位于0-12mm。在现场勘察中，并没有发现该墙体表面出现裂缝。墙体c各高度处的变形位于-18至12mm，变形略大于墙体b。如图11b所示，其墙体出现一条从上至下的裂纹，且墙体并未有裂缝出现。墙体d各高度处的变形分布在0-36mm。结合图6b可以发现，该面墙包含门窗等建筑物的主要采光结构，墙体结构较少，且裂缝出现在该墙体结构较少的区域。因此，该墙体各高度处的变形波动较小，没有出现与墙体a类似的变形波动趋势。
[0163]
边界点估计用时及结果对比
[0164]
针对墙体边界点估计，在基于法线的方法的基础上，进行了适应化。在本方法中，分别采用本文方法、基于法线的方法进行边界估计，并对比用时、结果。基于第一次观测中的墙体点云，获取估计时间和结果分别如表2和图12所示。
[0165]
表2边界点估计时间
[0166][0167][0168]
由表2可以看出，相对于基于法线的方法，本文方法所需时间约为其1/10。这是由于计算过程中，本文方法避免了计算每个点的邻域法线，同时也无需重新估算点云在每一个邻域坐标系下的坐标。仅需计算一次法线，转换一次坐标系。由图12可以看出，本文方法和基于法线的方法均可以成功估计墙体边界点。为进一步对比边界点估计结果，分别基于上述两种方法的估计结果，计算同位置结构长度的误差绝对值。根据图13所示结果，该误差绝对值均不超过8mm，大部分保持在0-5mm。由此可见，这两种方法估计的边界点具有较高的一致性，且本文方法所需时间远小于基于法线的方法。选权函数对结构长度计算的影响
[0169]
除本发明中的igg外，huber、hampel也是常用的选权函数。不同于将数据分为三段的igg选权函数，公式(16)所示的huber选权函数将数据分为两段，而公式(17)所示的hampel选权函数将数据分为四段。在本方法中，为验证不同选权函数的影响，基于第一次观测中的边界点，进行了结构长度的计算。将以huber、hampel为选权函数计算的结构长度，分别与以igg为选权函数计算的相对比后，获取同位置结构长度的误差绝对值。
[0170]
如图14所示结果，在huber与igg间的误差绝对值中，墙体a有1位置处约为6mm。在hampel与igg间的误差绝对值中，墙体b有1位置处约为7mm。其余各位置处的误差绝对值主要分布在0-5mm。这主要是由于上述三种选权函数的结构不同造成的，huber选权函数的第二段在igg选权函数中被划分为两部分，igg选权函数的第三段在hampel选权函数中被划分为两部分。但总体来说，上述三种选权函数对于结果的影响较小，均保持在5mm以内。
[0171][0172]
其中，s是保证稳健性的阈值系数，在本文中取1.5。
[0173][0174]
其中，s1、s2、s3是保证稳健性的阈值系数，在本文中分别取1.5、2.0、3.0。
[0175]
结构长度计算的准确度评价
[0176]
为评估准确度，在第二次观测当天，对建筑物墙体a进行了额外观测。同样基于墙体点云厚度，以h＝0.25m为兴趣区域高度，计算了墙体a的结构长度。并在获取两次观测间的变形后，计算变形绝对值与结构长度间的比值(即变形率)，其结果如表3所示。
[0177]
第二次观测和额外观测为同一天内的两次连续观测，其真实变形值应为0mm。因此，可将上表中变形值视为误差值，将变形率视为误差率，用于评价结果的准确性。从表3中可以看出，该两次连续观测间的变形值均分布在-6至6mm。这主要是由于扫描仪的观测误差、多站数据的拼接误差等偶然误差造成的。结合各位置的变形率不超过万分之5，可以表明本文方法获取的变形值具有较高的准确性，且可以达到mm级精度。
[0178]
表3准确性评价结果
[0179][0180][0181]
针对砖混建筑物，本方法提出了结构长度计算方法，并开展其在变形监测中的应用。经墙体点云厚度计算、建筑物墙体拆分、墙体边界点估计、上边界拟合及结构长度计算，
并结合案例研究及讨论，取得主要结论如下：
[0182]
(1)由本方法计算的结构长度而获取的变形情况，与墙体裂缝表现的变形情况基本吻合。
[0183]
(2)在边界点估计部分，本方法与基于法线的方法，均可以成功估计墙体边界点，且本方法用时仅为其1/10。在以估计结果为唯一变量时，计算的结构长度误差绝对值均不超过8mm，大部分保持在0-5mm。
[0184]
(3)在选权迭代最小二乘法中，以选权函数为唯一变量。将igg计算的结构长度为标准，huber、hampel计算结果与其之间的误差绝对值主要分布在0-5mm。
[0185]
(4)在假定两次间隔时间极短的连续观测间变形为0mm的情况下，本方法计算的结构长度误差分布在-6至6mm，误差率不超过万分之5，可以达到mm级精度。
[0186]
以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。

技术特征：
1.一种基于点云的砖混建筑结构长度计算方法，其特征在于，包括以下步骤：s1：获取建筑物墙体点云，对所述建筑物墙体点云进行预处理，获得处理后建筑物墙体点云；s2：在所述处理后建筑物墙体点云内设置多个邻域，并对多个所述邻域内的点云进行过滤，得到有效邻域点，并计算出墙体点云厚度值；s3：基于所述墙体点云厚度值确定第一距离阈值和第二距离阈值，基于所述第一距离阈值对所述处理后建筑物墙体点云进行建筑物点云分割，得到每个单面墙体的组成点云，基于所述第二距离阈值对每个所述单面墙体的组成点云进行墙体点云再分配，得到分配后墙体点云；s4：采用改进后的基于法线的点云边界估计方法对所述分配后墙体点云进行墙体边界点估计，得到墙体边界点；s5：对所述墙体边界点进行处理得到上边界线，在交点约束条件下采用选权迭代最小二乘法对所述上边界线进行拟合，得到拟合后上边界线；s6：基于所述拟合后上边界线，在墙体上创建兴趣区域，根据所述兴趣区域内垂直边界点所在直线，计算得到墙体结构长度。2.根据权利要求1所述的一种基于点云的砖混建筑结构长度计算方法，其特征在于：所述步骤s1中获取建筑物墙体点云，对所述建筑物墙体点云进行预处理，获得处理后建筑物墙体点云，具体过程如下：对建筑物进行扫描，得到建筑物墙体点云，对所述建筑物墙体点云进行降噪、抽稀、地面滤波处理，得到处理后建筑物墙体点云。3.根据权利要求1所述的一种基于点云的砖混建筑结构长度计算方法，其特征在于：所述步骤s2中在所述处理后建筑物墙体点云内设置多个邻域，并对多个所述邻域内的点云进行过滤，得到有效邻域点，并计算出墙体点云厚度值，具体过程如下：在所述处理后建筑物墙体点云内随机选取n个点，每个点设定以r为半径划定邻域；计算邻域内的点到该邻域所在平面的距离d，将所述距离d小于距离中误差的点作为有效邻域点，计算出每个邻域的点云厚度th
nei
，如下式所示：，如下式所示：其中，rmse1表示单个邻域的距离中误差；n、m分别为单个邻域中的原始点的数量和有效邻域点的数量；[d1,d2,
…
,d
i
,
…
,d
n
]和[d1,d2,
…
,d
j
,
…
,d
m
]分别表示单个邻域中的原始点至该邻域所在平面距离的集合和单个邻域中的有效邻域点至该邻域所在平面距离的集合；d
mean
表示集合[d1,d2,
…
,d
j
,
…
,d
m
]中有效邻域点至邻域所在平面距离的均值；将每个邻域的所述点云厚度th
nei
进行均值计算，得到邻域点云厚度均值，取小于厚度中误差的邻域为有效邻域，计算得到所述墙体点云厚度值th
wall
，如下式所示：，如下式所示：其中，rmse2表示邻域点云厚度中误差，n、m分别为原始的邻域数量和有效邻域的数量；
和和分别为原始的邻域的点云厚度集合和有效邻域的点云厚度集合；为集合中的有效邻域的点云厚度均值。4.根据权利要求1所述的一种基于点云的砖混建筑结构长度计算方法，其特征在于：所述步骤s3中基于所述墙体点云厚度值确定第一距离阈值和第二距离阈值，基于所述第一距离阈值对所述处理后建筑物墙体点云进行建筑物点云分割，得到每个单面墙体的组成点云，基于所述第二距离阈值对每个所述单面墙体的组成点云进行墙体点云再分配，得到分配后墙体点云，具体过程如下：基于所述墙体点云厚度值确定第一距离阈值和第二距离阈值，所述第一距离阈值大于所述墙体点云厚度值，所述第二距离阈值等于所述墙体点云厚度值；将所述处理后建筑物墙体点云投影至xoy面，基于第一距离阈值和随机抽样一致性算法对投影至xoy面的点云进行处理，得到单面墙体的组成点云；从所述处理后建筑物墙体点云中分离所述单面墙体的组成点云，完成数据集合的一次更新；重复随机抽样一致性算法处理步骤和数据集合更新步骤，直至获取每个单面墙体的组成点云；采用随机抽样一致性算法对每个所述单面墙体的组成点云进行处理，获得该单面墙体点云；对单面墙体的组成点云和单面墙体点云求差集，得到差集点云；计算差集点云中每个点至每个单面墙体的距离；取所述差集点云中的每一个到单面墙体的距离小于第二距离阈值的点，分配至该单面墙体的点云中，得到每个单面墙体的分配后墙体点云。5.根据权利要求1所述的一种基于点云的砖混建筑结构长度计算方法，其特征在于：所述步骤s4中采用改进后的基于法线的点云边界估计方法对所述分配后墙体点云进行墙体边界点估计，得到墙体边界点，具体过程如下：采用选权迭代最小二乘法对所述分配后墙体点云进行计算，得到平面的模型参数；基于平面的模型参数分别计算出分配后墙体点云与xoz面、xoy面间的旋转角，如下式所示：θ＝x
a
x
b
/|x
a
||x
b
|；其中，θ为旋转角；||为模型参数的模长，其中x
a
、x
b
分别为预设的已知的平面a和平面b的模型参数；基于旋转角将分配后墙体点云旋转至与xoy面平行，得到旋转后的墙体点云；在旋转后的墙体点云中选择一点，将该点与其邻域点连线，计算相邻连线间的夹角，并按照角度大小排列；计算排列后的相邻的夹角之间的差值，以差值大于预设的角度阈值的点为墙体边界点。6.根据权利要求5所述的一种基于点云的砖混建筑结构长度计算方法，其特征在于：所述步骤s5中对所述墙体边界点进行处理得到上边界线，在交点约束条件下采用选权迭代最小二乘法对所述上边界线进行拟合，得到拟合后上边界线，具体过程如下：将所述墙体边界点分为内边界点和外边界点，得到外边界点所在直线的边界线，包括垂直边界线和水平边界线，将与屋顶毗邻的水平边界线作为上边界线；将所述上边界线垂直投影至xoy面，并采用选权迭代最小二乘法求解投影至xoy面上的上边界线的模型参数；
基于xoy面上的上边界线的模型参数分别计算出上边界线所在墙体平面与xoz面、yoz面间的旋转角；将上边界线水平投影至旋转角较小的xoz面或yoz面上，采用选权迭代最小二乘法求解上边界线投影后的直线的模型参数，上边界线在两个面上的线a和线b分别由下列两个直线方程式表示：k1x+b1＝y＝(z-b2)/k2；k3x+b3＝y＝(z-b4)/k4；其中，k1、k2、b1、b2是线a在空间坐标系中的方程表达式参数，k3、k4、b3、b4是线b在空间坐标系中的方程表达式参数；对两个直线方程进行求解得到交点，设点a、点b为求解得到的交点，则拟合后上边界线如下式所示：其中，(x
a
,y
a
,z
a
)、(x
b
,y
b
,z
b
)分别为边界线交点a和点b的坐标。7.根据权利要求6所述的一种基于点云的砖混建筑结构长度计算方法，其特征在于：所述步骤s6中基于所述拟合后上边界线，在墙体上创建兴趣区域，根据所述兴趣区域内垂直边界点所在直线，计算得到墙体结构长度，具体过程如下：设定步长为h，将拟合后上边界线依次沿墙体平面向下平移，并以每次平移后的拟合后上边界线为中心线，将高度为2h的区域划定为兴趣区域；筛选出所述兴趣区域内位于该墙体本身的边界点和位于相邻墙体的边界点；将位于墙体本身的边界点和位于相邻墙体的边界点合并，并采用选权迭代最小二乘法计算出兴趣区域中的垂直边界线的模型参数，如下式所示：k5x+b5＝k6y+b6＝z；其中，k5、k6、b5、b6是垂直边界线的模型参数在空间坐标系中的方程表达式参数；设定平移次数为i，第i次平移后的拟合后上边界线如下式所示：(x-f)/g＝y＝(x-h-hi)/l；基于垂直边界线的模型参数和第i次平移后的拟合后上边界线，设定点a和点b是兴趣区域内两交点，则墙体长度为：其中，l为当前平移位置的墙体结构长度；(x
a
,y
a
,z
a
)、(x
b
,y
b
,z
b
)分别为兴趣区域内交点a和点b的坐标。

技术总结
本发明涉及一种基于点云的砖混建筑结构长度计算方法，包括以下步骤：墙体点云厚度计算、建筑物墙体拆分、墙体边界点估计、上边界拟合、结构长度计算。计算点云厚度时，随机选择多个邻域，但仅保留部分邻域点、部分邻域；基于墙体点云厚度，进行建筑点云分割、墙体点云再分配；以单面墙体为基本单位，进行了墙体边界点估计；基于上边界点，将选权迭代最小二乘法引入，进行该线拟合；由上边界线，在墙体上创建兴趣区域，结合区域内垂直边界点所在直线，计算墙体结构长度。本发明通过选择结构长度作为监测指标，结合点云监测的优势和砖混建筑物的几何结构特点，可以十分精确地计算出砖混建筑结构长度，基于结构长度可更好地监测砖混建筑的变形情况。变形情况。变形情况。


技术研发人员：王磊 张鲜妮 李靖宇 黄金中
受保护的技术使用者：安徽理工大学
技术研发日：2022.12.08
技术公布日：2023/7/11