一种边坡点云曲面重建及变形信息提取方法与流程

本发明属于变形监测领域，具体涉及一种边坡点云曲面重建及变形信息提取方法。

背景技术：

20世纪90年代中期出现然后极为迅速地发展起来并且被人熟知的三维激光扫描技术，据有其他传统边坡变形监测方法无法比拟的优势；三维激光扫描技术获取高密度高精度的三维点云数据高效快速，对监测结果进行研究掌握其变形规律有很大帮助。在三维激光扫描应用于边坡变形监测的相关技术之中，散乱的边坡点云曲面重构技术和边坡点云变形信息提取方法选取是关键。

1971年，Hardy在研究曲面的二元函数中使用了MQ基函数来构造和拟合曲面，在他研究之后许多学者围绕着他的方法又研究出大量的成果，为后续被人广泛应用的的径向基函数的出现奠定了大量的理论研究基础；屈永华建立了一种适用于地质曲面重建的径向基函数(RBF)神经网络模型，并给出了具体算法和一般曲面插值方法相比，该模型对原始数据的分布和边界条件无特别的要求，因而适合于遭到破坏或勘探数据较少的地质曲面的重建，但这种进行的点云曲面重建的算法只适用于数据量小的散乱点云数据，曲面重建的精度比较低。

2010年赵小平等对边坡点云数据进行了处理，再通过扫描数据处理软件获得DEM数据边坡形态，他的研究中验证了点云预处理之后的精度，为边坡变形监测与灾害预报提供基础数据。2015年刘昌军等提出基于八叉树结构的点云与点云的精细直接比较算法，实现变形区域多个点云数据的自动比较，其研究运用的精细自动比较算法随点云的数据量增加，效率降低，不能快速高效的提取复杂大量数据的边坡点云变形信息。

技术实现要素：

针对现有技术的不足，本发明提出一种边坡点云曲面重建及变形信息提取方法，以达到提高提取变形信息效率和正确性目的。

一种边坡点云曲面重建及变形信息提取方法，包括以下步骤：

步骤1、通过三维激光扫描仪获取边坡三维点云数据，并通过全站仪测出每个反射片的中心坐标；

步骤2、根据实际需求间隔一段时间，再一次通过三维激光扫描仪获取边坡三维点云数据，并通过全站仪测出每个反射片的中心坐标；

步骤3、对三维激光扫描仪获取的边坡三维点云数据进行去噪和拼接预处理；

步骤4、基于最小包围盒的RBF点云曲面重建，构建边坡点云曲面；

步骤5、对三维激光扫描仪和全站仪测得的两次点云坐标进行比较，并根据预处理过的两期三维点云数据生成DEM模型，将两个DEM模型进行比较，获得边坡变形信息。

步骤1所述的通过三维激光扫描仪获取边坡三维点云数据，并通过全站仪测出每个反射片的中心坐标，具体为：三维激光扫描仪扫描时每站需布设有4～6个反射片。

步骤4所述的基于最小包围盒的RBF点云曲面重建，构建边坡点云曲面，具体步骤如下：

步骤4-1、输入预处理后的点云数据；

步骤4-2、用坐标系旋转确定最小包围盒，具体如下：

设三个坐标轴的初始转角均为0，按设定的角度步长，用三重循环分别使三个坐标由0°至360°进行旋转；在每一个旋转位置，确定包围盒的边界数据和体积，求出体积最小包围盒及相对应的转角；

步骤4-3、沿坐标系x轴方向分解包围盒，若点云数据集的点数大于设定的包围盒内最多点数，并且包围盒的边长度大于设定的最小盒子边长，则获取点云数据集所有点的重心；

步骤4-4、采用重心x坐标所在平面对点云数据集和包围盒进行切分，若切分后包围盒内的数据集小于设定的包围盒内最多点数，或者在x方向的包围盒边长小于等于设定的最小盒子边长，则执行步骤4-5，否则返回执行步骤4-4；

步骤4-5、对切分获得的相邻包围盒及相邻的数据集求并集，若某个包围盒内的数据大于设定的包围盒内最多点数，则在y轴和z轴方向上分割包围盒，直到包围盒内的数据均小于等于设定的包围盒内最多点数，且所有包围盒的并集包含该包围盒；

步骤4-6、采用径向基函数插值方法对分割后的包围盒的数据进行插值加密，根据插值后的数据构建边坡点云曲面。

步骤5所述的对三维激光扫描仪和全站仪测得的两次点云坐标进行比较，并根据预处理过的两期三维点云数据生成DEM模型，将两个DEM模型进行比较，获得边坡变形信息；

具体步骤如下：

步骤5-1、判断三维激光扫描仪得到的反射片的坐标值与全站仪得到的反射片的坐标值的差值是否在设定范围内，若是，则三维激光扫描仪获得的数据能满足边坡变形监测的要求，否则三维激光扫描仪重新采集数据；

步骤5-2、将预处理后的三维激光扫描仪两次采集的点云数据导入到Cyclone软件中；

步骤5-3、利用Cyclone软件的mesh模块生成两个点云的DEM模型；

步骤5-4、将两个点云的DEM模型进行插值运算，找到高程相同的点，连线生成高密度等高线图，并将其另存为dwf格式的文件；

步骤5-5、用CAD打开dwf格式的文件，并对变形信息进行提取和标注，得到用于分析变形的高密度等高线文件；

步骤5-6、以第一次扫描的DEM模型当做参考对象，第二次扫描的DEM模型当成结果对象，使用Geomagic Qualify软件对两个DEM模型进行比较，获得边坡变形信息。

本发明优点：

本专利针对边坡点云数据模型重建问题，采取了改进的RBF算法，通过对经过点云预处理之后的边坡点云数据进行了曲面建模，提高了RBF算法建模效率和精度；针对边坡点云信息提取问题，提出Cyclone结合CAD提取DEM与等高线法、三维质检软件法两种方法的可行性，并将三维质检软件法提取的变形数据与传统变形监测数据进行了对比，验证了方法的正确性。

附图说明

图1为本发明一种实施例的边坡点云曲面重建及变形信息提取方法流程图；

图2为本发明一种实施例的站位和控制点布局简略示意图；

图3为本发明一种实施例的最小包围盒的3种形式示意图，其中，图(a)为AABB包围盒示意图，图(b)为0BB包围盒示意图，图(c)为凸包围盒示意图；

图4为本发明一种实施例的确定最小包围盒的流程图；

图5为本发明一种实施例的点云边坡曲面与点云对比图，其中，图(a)是未被改进算法构建的点云边坡曲面与点云对比图，图(b)是改进算法构建的点云边坡曲面与点云对比图；

图6为本发明一种实施例的等高线文件示意图，其中，图(a)是第一次扫描的等高线，图(b)是第二次扫描的等高线；

图7为本发明一种实施例的边坡变形示意图。

具体实施方式

下面结合附图对本发明一种实施例做进一步说明。

本发明实施例中，边坡点云曲面重建及变形信息提取方法，如图1所示，包括以下步骤：

步骤1、通过三维激光扫描仪获取边坡三维点云数据，并通过全站仪测出每个反射片的中心坐标；

本发明实施例中，根据地形条件初步定出设站位置；本发明实施例中扫描对象为门头沟一个尾矿的边坡，环境较为复杂；尾矿底部为相对稳固的地面，可以将三维激光扫描仪架设在上面；为保证数据的完整性，同时为了避免数据重叠，总共设站6站，每一测站保证周边能观测到布设有4～6个反射体，并且在边坡上固定6个反射体(这六个反射体中心点坐标不进行点云拼接)，其位置必须稳固，不易移动；并用全站仪测出每一个测站每一个反射片的中心坐标。

步骤2、一个月后对尾矿边坡进行第二次扫描，再一次通过三维激光扫描仪获取边坡三维点云数据，并通过全站仪测出每个反射片的中心坐标；

步骤3、对三维激光扫描仪获取的边坡三维点云数据进行去噪和拼接预处理；

步骤3-1、由于尾矿边坡点云既包含坡上杂草、小灌木、电力设施和岩石早点，又包含三维激光扫描系统自身存在着一些系统误差，所以点云情况较复杂，本发明实施例中提出了基于特征区域分类的双边滤波点云去噪方法这一改进方法，对点云数据进行去噪；

具体如下：

双边滤波方法定义为p′＝p-αn，其中，p为原点云数据，p_j为其邻域，p′为去噪后的点云数据，α为双边滤波因子，n为法矢的方向；

其中，δ_c与δ_s为切平面上高斯滤波系数，W_c为光顺滤波权函数，W_s为特征保持权重函数；特征区域的选择，对大于阈值的点云网格，进行基于最小二乘的曲面拟合计算，将点云拟合到一个曲面上；

基于最小二乘的曲面拟合要求拟合曲面与实际值差值的平方和达到最小值，即：

z＝f(x，y) (4)

式中，E(f)为误差平方和，f(l_i)＝a_jb^j(l_i)表示待求的拟合曲面，z＝f(x，y)表示实际值，要使E(f)达到最小，则必须满足进一步计算可以得到曲面未知系数，即可以求出点云网格的拟合曲面。

点云特征区域就是曲面极值点所围成的区域，此时极值点z_i满足Q命题：在一定邻域内，若z_t≥z_i，则z_i为极大值，若z_t≤z_i，则z_i为极小值；根据求得的区域进行平滑去噪(z_t表示邻域点)。

步骤3-2、采用三点到三点的拼接方法，在扫描区域中设置控制点，使相邻测站点云数据中至少有三个同名控制点，如图2所示为站位和控制点布局简略图；通过控制点的强制附合，可以将相邻的点云数据统一到同一个坐标系下，对多个测站点云数据进行拼接时，需要将多个点云数据中的控制点组成一个闭合环，可以有效地防止拼接过程中坐标转换误差的积累；

步骤3-3、利用Geomagic 2013软件对点云进行等间距简化；

步骤4、基于最小包围盒的RBF点云曲面重建，构建边坡点云曲面；

具体步骤如下：

步骤4-1、输入预处理后的点云数据；

步骤4-2、用坐标系旋转确定最小包围盒；

本发明实施例中，如图3所示，考虑到紧密性和相交测试算法，本实施例中选择了凸包围盒作为最小包围盒，确定最小包围盒的具体步骤参见图4(包括：初始化、选择实体、绕z轴旋转、绕y轴旋转、绕x轴旋转、提取实体质量特征、提取包围盒边界数据、计算包围盒体积、判断包围盒体积是否小于0时的包围盒体积，若是，则把Vi记作最小包围盒体积，；否则，则把V0的值赋予Vi，判断x轴是否旋转360°，判断y轴是否旋转360°，判断z轴是否旋转360°，若是，输出数据，否则返回继续旋转)；

本发明实施例中，设三个坐标轴的转角α，β，γ分别为0，按一定的角度步长(3个轴的角度步长是一致的，只是表示的方式不同，本实例中角度步长取值5度，用三重循环分别使三个坐标进行0°到360°的旋转；在每一个旋转位置，确定包围盒的边界数据{X_i}和体积V_i，其中{X_i}＝{X_min，X_max，Y_min，Y_max，Z_min，Z_max}，其中，X_min表示包围盒的边界X轴最小值，X_max表示包围盒的边界X轴最大值，Y_min表示包围盒的边界Y轴最小值，Y_max表示包围盒的边界Y轴最大值，Z_min表示包围盒的边界Z轴最小值，Z_max表示包围盒的边界Z轴最大值，根据求得的体积V_i，确定最小包围盒及相对应的转角α₀，β₀，γ₀；

步骤4-3、沿坐标系x轴方向分解包围盒，若点云数据集的点数|S|大于设定的包围盒内最多点数N_max(本实例中取值15)，并且包围盒的边长度大于设定的最小盒子边长l(本实例中取值为4mm)，则获取点云数据集所有点的重心(x_c，y_c，z_c)；

步骤4-4、采用重心x坐标所在平面x＝x_c对点云数据集和包围盒进行切分，若切分后包围盒内的数据集小于设定的包围盒内最多点数N_max，或者在x方向的包围盒边长小于等于设定的最小盒子边长l，则执行步骤4-5，否则返回执行步骤4-4；

步骤4-5、对切分获得的相邻包围盒{B_i}及相邻的数据集{S_j}求并集，若某个包围盒内的数据大于设定的包围盒内最多点数N_max，则在y轴和z轴方向上分割包围盒，直到包围盒内的数据均小于等于设定的包围盒内最多点数N_max，且所有包围盒的并集包含该包围盒；

步骤4-6、采用径向基函数插值方法对分割后的包围盒的数据进行插值加密，根据插值后的数据构建边坡点云曲面；

本发明实施例中，对分割后的包围盒进行径向基函数插值，得到函数集{f_i}；其中基函数的选择对曲面构建影响很大，所以本发明实施例对基函数进行了改进，目前应用最好的基函数为：

φ(p_i，p_j)＝|p_i-p_j|³ (5)

其中，φ(p_i，p_j)表示函数，p_i、p_j表示邻域点；

给定数据点集{p_i}_i＝1，…n(n表示点的个数)和数集{d_i}，求f：R³→R(R表示实数集)，使得函数f(p_i)：

f(p_i)＝d_i，i＝1，…，n (6)

一般的径向基函数是基本函数Φ：R³×R³→R和多项式的线性组合：

P(x，y，z)是多项式，为了求解满足插值条件的权因子w_i，将p_i代入上式中经过计算得到下面关于函数与权因子的关系式，即本发明改进后的插值加密公式：

当点集很大的时候，方程组很庞大，计算机求解上述方程组不可能完成；为了获得更加鲁棒的结果，需要把内外两个等距面上的点作为插值点；这样就需要计算点云的法向量和距离；

本发明实施例中，对函数集{f_i}求和，构建边坡点云曲面，如图5所示，是未被改进算法(图5中图(a))和改进算法(图5中图(b))构建的点云边坡曲面与点云对比图；

步骤5、对三维激光扫描仪和全站仪测得的两次点云坐标进行比较，并根据预处理过的两期三维点云数据生成DEM模型，将两个DEM模型进行比较，获得边坡变形信息。

具体步骤如下：

步骤5-1、判断三维激光扫描仪得到的反射片的坐标值与全站仪得到的反射片的坐标值的差值是否在2mm内，若是，则三维激光扫描仪获得的数据能满足边坡变形监测的要求，否则三维激光扫描仪重新采集数据；

步骤5-2、将预处理后的三维激光扫描仪两次采集的点云数据导入到Cyclone软件中；

步骤5-3、利用Cyclone软件的mesh模块生成两个点云的DEM模型；

步骤5-4、将两个点云的DEM模型进行插值运算，找到高程相同的点，连线生成高密度等高线图，并将其另存为dwf格式的文件；

步骤5-5、用CAD打开dwf格式的文件，并对变形信息进行提取和标注，得到用于分析变形的高密度等高线文件；如图6所示，其中，图(a)是第一次扫描的等高线，图(b)是第二次扫描的等高线；

步骤5-6、以第一次扫描的DEM模型当做参考对象，第二次扫描的DEM模型当成结果对象，使用Geomagic Qualify软件对两个DEM模型进行比较，如图7所示，获得边坡变形信息。