一种基于无人机影像序列的煤场采挖量计算方法与流程

本发明涉及一种体积测量方法，特别是涉及一种结合无人机倾斜摄影技术、测绘技术、计算机图形学和计算机数字几何处理技术的煤场采挖量体积计算方法。

背景技术：

大型露天煤场存煤量计算是煤场日常生产管理的重要环节，是企业测量验收工作中的重要指标，直接关乎企业的经济和效益。由于传统计算方法的缺点和不足，使得精确、快速的计算煤场工程量变得尤为重要。目前国内外主要采用土方计算原理，主要方法有：断面法、格网法、等高线法以及dtm三角网法。断面法是通过剖面分段化来计算采挖工程量。该方法操作繁琐、工作量大，主要使用于狭长地带等。格网法是通过设置格网间距和测量点位高程来计算采挖工程量。该方法工作量大，精度主要取决于地形环境，限制因素多。等高线法是通过利用等高线所围面积与相邻等高线高差的乘积得到采挖工程量，但等高线必须是闭合曲线，有时需要与断面法相结合，比较适合坡度变化较多的场合。目前，在实际生产中，dtm三角网法被广泛采用。该方法首先通过全站仪等设备测量煤场区域三维点数据，再通过实际采挖量生成空间地面模型，最后计算所有斜三棱柱体积得到采挖工程量。该方法操作相对复杂，自动化程度不高。

技术实现要素：

本发明的主要目的在于，针对现有技术中的不足，提供一种结合无人机倾斜摄影技术、测绘技术、计算机图形学和计算机数字几何处理技术的煤场采挖量体积计算方法，本文方法自动化程度高，能够自动定位到目标区域，并计算，使数据处理更加方便、快捷。

结合上述内容，本发明技术线路如下：

一种基于无人机影像序列的煤场采挖量计算方法，包括以下步骤：

1)煤场区域环境勘察；

采集好无人机起飞和降落点，并确定好航拍区域；

2)布设控制点；

在待测煤场周围布设控制点，采集这些点的三维坐标并做好记录；

3)无人机飞行作业；

采用四旋翼无人机进行航摄作业，采用倾斜摄影技术对煤场区域进行航拍，从而获取大量影像照片；

4)获取煤场开挖前三维点云数据；

在获取煤场倾斜数据之后，需要从图像序列中获取三维信息，此处采用运动过程恢复结构(structurefrommotion,sfm)重建出煤场三维结构模型，并得到点云数据。在重建时最终得到的是所有匹配的特征点的三维坐标，存在图像内容的特征点信息不足，并且只能检测到部分的特征点等原因，以至于恢复的三维场景点云信息较稀疏。此时通过多目立体视觉算法(patch-basedmultipleviewstereopsis,pmvs)把稀疏的三维点云转为致密的三维点云，最后实现三维重建；5)获取煤场开挖后三维点云数据；

与步骤四方法相同，先采用sfm算法得到稀疏的煤场三维点云数据，再用pmvs算法将稀疏的三维点云进行加密，从而得到致密的三维点云数据；

6)两期点云配准；

将煤场开挖前的点云数据和开挖后的点云数据进行配准，首先采用基于特征点的粗配准方法，通过布设在煤场周围的控制点来得到具体的坐标信息，进行坐标转换，选取同名标识点将前期点云数据转换到统一坐标系下，实现点云的粗配准。粗配准之后，前后两期点云数据已经统一在同一个坐标系下，此时采用迭代最近点匹配算法，icp算法(interativeclosestpoints,icp)对两期点云进行精配准。

7)开挖区域精确定位和点云模型边界确定；

煤场开挖过程中，同一个煤场存在多处作业的情况，对于整个煤场而言，采挖量的体积变化是指多个作业区域煤量变化的综合效应。此处首先对煤场不同采挖区域进行精确定位，从而对不同区域进行体积计算。再通过算法实现两期点云的边界确定；

8)点云网格化；

将得到的点云模型进行泊松重建，得到一个精确的封闭网格模型。

9)煤场开挖工程量计算；

对封闭的网格进行网格体积计算。

2、根据权利要求1所述的一种基于无人机影像序列的煤场采挖量计算方法，其特征在于：所述步骤3)中的无人机选取四旋翼无人机，搭载高清摄像头，调整好飞行高度。

3、根据权利要求1所述的一种基于无人机影像序列的煤场采挖量计算方法，其特征在于：所述步骤4)和步骤5)中sfm算法和pmvs算法，具体实现为，对于无人机采集的煤场影像序列，首先采用基于sfm的方法从原始影像序列中提取出特征点集，并进行特征点匹配，计算出模型的几何结构和相机的运动参数，最后重建出稀疏的点云三维结构。

影像数据经过结构和运动恢复之后，便得到了稀疏的三维点云，再采用多目立体视觉算法对稀疏点云进行表面重建。

算法首先进行特征点匹配，对于任意一副图像为ii，相机的光心是o(ii)，对于ii的特征点集合f,对于特征点f∈f，通过在其他的图像中找到对应的特征点f′，此时有特征点对(f,f′)，并将这个三维空间点赋值给c(p)，并得到与之相应的参数。

其次，进行面片扩展，其目的是为了得到更多的面片，即更加致密的点云。对于每个面片p，在每个可视的图像中，进行四领域扩展。当领域网格集合n(p)中某个网格已经有投影的面片，此时该网格就不需要扩展，即将其去除。

最后进行面片滤波，此处通过三种滤波器来进行滤波。首先进行可视一致性滤波，即：

其中，u(p)定义为与当前的面片p可见情况不一致的集合p′，p与p′不在同一个深度中，却投影到了同一个网格中。面片p如果发生异常，就会被滤掉。

第二个滤波器是在v^*(p)的集合中，利用每个摄像机的位置进行深度测试，观察是否有面片被遮挡，计算机没被遮挡的面片的数量，当小于规定阈值时，就将该面片删除。

最后进行几何一致性滤波，如面片p，当面片p的邻域的面片总和，与v(p)中面片p所在的网格和四邻近所包含的面片比例小于阈值时，就将该面片删除。

通过上述步骤，最后得到了致密的三维点云模型，并开始对这两期点云进行数字几何处理。

4、根据权利要求1所述的一种基于无人机影像序列的煤场采挖量计算方法，其特征在于：所述步骤6)中的前后两期点云配准，具体为，

首先进行点云粗配准，粗配准的目的在于缩小点云之间平移误差和旋转误差，为精配准做好准备。本文通过在煤场周围均匀布设控制点来获取点的坐标信息，进行坐标转换。以后期点云为基准，通过选取同名标识点将前期点云数据转换到统一坐标系下，实现点云的粗配准。

粗配准之后，前后两期点云数据已经统一在同一个坐标系下，此时需要采用精配准进行最终处理，以提高配准精度。此处采用迭代最近点匹配算法，icp算法(interativeclosestpoints,icp)对两期点云进行精配准。同时，煤场开挖时，使得前后两期点云产生形变，本文采用基于形态不变区的点云配准方法，对icp算法进行相应的调整以适用于本文场合。

icp算法是通过寻找两个点集的对应匹配关系，计算坐标转换参数，来完成匹配过程。通过确定对应关系点集并计算最优刚体变换的思路，使其满足收敛准则。由于外界环境影响，在煤场采挖过程中，前后两期点云形态可能发生了变化，在全局点云配准时，可能出现局部最优，从而影响配准的精度。因此，本文利用在两期点云数据保持不变的区域作为基准点，通过icp算法计算全局最优的转换参数，将这些转换参数对形态发生变化的区域进行刚性转换，从而实现点云精配准。5、根据权利要求1所述的一种基于无人机影像序列的煤场采挖量计算方法，其特征在于：所述步骤7)对煤场开挖区域精确定位和点云模型边界确定，具体为，煤场开挖过程中，同一个煤场存在多处作业的情况，对于整个煤场而言，采挖量的体积变化是指多个作业区域煤量变化的综合效应。故此处首先需要对煤场不同采挖区域进行精确定位，定位的同时，算法对两期点云重合之外部分形成的区域的边界进行寻找。算法流程如下，设前期点集为a，设p为点集a中的点，p∈a。后期点集为b，q为点集b中的点，q∈b。定义q^*为集合b中到点p距离最近的点。并有q^*＝argmin‖p-q‖，并设定一个阈值，当‖p-q^*‖≥δd时，把该点p存入点集a′中。同时，定义p^*为集合a中到点q距离最近的点。并有p^*＝argmin‖q-p‖，当‖q-p^*‖≥δd时，把该点q存入点集b′中。最终得到的点集a′、b′即为目标点集。

5-1)建立点云集合a、b、a′、b′，前期点云形成集合a，后期点云形成集合b，并建立两个新的点云集合a′、b′；

5-2)计算点间最小距离：在点集b中找到点q^*，使得点q^*满足q^*＝argmin‖p-q‖。其中，p为点集a中的点；q、q^*为点集b中的点；

5-3)寻找目标点集a′：在点集a、b中，当满足‖p-q^*‖≥δd时，且q^*＝argmin‖p-q‖，则把该点p存入目标点集a′中。其中，δd为距离约束阈值。

5-4)再次计算点间最小距离：在点集a中找到点p^*，使得点p^*满足p^*＝argmin‖q-p‖。

5-5)寻找目标点集b′：在点集a、b中，当满足‖q-p^*‖≥δd时，且p^*＝argmin‖q-p‖，则把该点q存入目标点集b′中。

5-6)重复步骤(2)、(3)、(4)、(5)：循环点集a、b中所哟点；

5-7)将得到的点集a′、b′进行合并，获取不同采挖区域的点云模型边界，并得到最终点云集。

6、根据权利要求1所述的一种基于无人机影像序列的煤场采挖量计算方法，其特征在于：所述步骤8)点云网格化，具体为，

将点集a′、b′合并之后，此时对合并后的点集从点云到网格进行poisson重建，得到一个精确的封闭的网格模型。泊松表面重建算法具体流程如下：

6-1)首先进行问题离散化，此处通过使用自适应的八叉树来表示隐式函数，并用样本点的位置采用一颗八叉树o，之后给八叉树o中的每一个节点o∈o附加一个函数f0。然后，定义一个基函数f：r³→r，使该函数空间可进行距离平移和尺度缩放。对于每一个节点o∈o，定义函数f0为节点中心o，以o的尺寸展开：

其中，o.c为节点o的中心，o.w为节点o的宽度。

此处函数f定义为一个带通滤波器与它自身的n次卷积：

f(x,y,z)≡(b(x)b(y)b(z))^*n

6-2)其次，定义向量场，为了满足节点的精度要求，采用三相线差值法把样本点分配到八个最邻近节点。

6-3)最后进行泊松问题求解，通过求解函数χ，使得相应方程达到最小值。最后进行等值面提取，选择逼近输入的采样点位置来提取等值面，有采样点的位置估计χ，之后使用平均值来提取等值面，最终得到网格模型。

与传统方法相比，实施本发明的一种结合无人机倾斜摄影技术、计算机图形学和计算机数字几何处理技术的煤场采挖量体积计算方法，具有以下有益效果：

本发明设计了一套完整的煤场体积测量方法，从煤场两期点云的数据获取到煤场体积开挖的计算，该方法快速、高效、精度高，能够满足煤场日常管理需求。

该方法操作方便，自动化程度高，不需要人为干预，从而减少了人为误差，使得煤场的生产作业更加高效。

本发明自主设计了一套经过配准后的两期点云中被开挖部分点云边界的确定的算法，该算法能够实现对开挖前、后的两期点云中，被开挖部分的点云进行精确定位和边界寻找，从而得到被开挖部分点云的精确模型，这是本设计发明的关键一步。

由于在日常作业时，同一个煤场中存在多处作业的情况发生，此时就会出现多处采挖区域，本发明能够同时精确的定位到不同的采挖区域，进而对这些区域进行采挖区域边界确定。

附图说明

图1为本发明一种基于无人机影像序列的大型露天煤场采挖量计算方法的流程图；

图2是本发明检测方法中三维体积计算的原理示意图；

图3是本发明检测方法中实验区及数据采集现场图；

图4为本发明检测方法中场景一下点云配准结果图；

图5是本发明检测方法中场景二下点云配准结果图；

图6是本发明检测方法中泊松重建之后得到的封闭网格模型图。

具体实施方式

下面结合说明书附图，对本发明作进一步的说明。

本发明提供一种基于无人机影像序列的大型露天煤场采挖量计算方法，如图1所示，包括以下步骤：

1)煤场区域环境勘察；

在对煤场环境进行勘察时，注意煤场附近有无高层建筑，是否会对无人机飞行造成影响，并确保不会发生信号干扰。在天气良好的情况下进行飞行，以及是否有煤场区域发生遮挡等。

2)布设控制点；

使用测绘专用设备rtk设备对煤场区域周围进行控制点测量，并有序做好坐标记录。控制点尽量呈矩形布设，选择地面特征点明显的位置进行布设，控制点数量尽量在四个及以上；

3)无人机飞行作业；

采用四旋翼无人机进行航摄作业，无人机选取四旋翼无人机，搭载高清摄像头，选取合适的飞行高度。首先在无人机地面站中规划好飞行区域，如果一个架次不能完成整个煤场的区域飞行，则需进行多架次飞行。划定好飞行区域之后，再进行飞行速度和飞行高度等的设定，对于航线而言，调整好无人机飞行的航向重叠率和纵向重叠率，且不同架次之间也要保证相应的航线重叠率。参数设置好之后，无人机开始自主飞行。对于部分区域，可能存在无人机自主飞行存在盲区，部分区域无法获取照片，此时就需要手动操作无人机进行飞行补拍。

4)获取煤场开挖前三维点云数据；

在获取煤场倾斜数据之后，利用sfm算法重建出煤场稀疏三维点云结构模型，再利用pmvs算法把稀疏的三维点云转为致密的三维点云，最后实现三维重建；

5)获取煤场开挖后三维点云数据；

与步骤四方法相同，先采用sfm算法得到稀疏的煤场三维点云数据，再用pmvs算法将稀疏的三维点云进行加密，从而得到致密的三维点云数据；

算法首先进行特征点匹配，对于任意一副图像为ii，相机的光心是o(ii)，对于ii的特征点集合f,对于特征点f∈f，通过在其他的图像中找到对应的特征点f′，此时有特征点对(f,f′)，并将这个三维空间点赋值给c(p)，并得到与之相应的参数。

其次，进行面片扩展，其目的是为了得到更多的面片，即更加致密的点云。对于每个面片p，在每个可视的图像中，进行四领域扩展。当领域网格集合n(p)中某个网格已经有投影的面片，此时该网格就不需要扩展，即将其去除。

最后进行面片滤波，此处通过三种滤波器来进行滤波。首先进行可视一致性滤波，即：

其中，u(p)定义为与当前的面片p可见情况不一致的集合p′，p与p′不在同一个深度中，却投影到了同一个网格中。面片p如果发生异常，就会被滤掉。

第二个滤波器是在v^*(p)的集合中，利用每个摄像机的位置进行深度测试，观察是否有面片被遮挡，计算机没被遮挡的面片的数量，当小于规定阈值时，就将该面片删除。

最后进行几何一致性滤波，如面片p，当面片p的邻域的面片总和，与v(p)中面片p所在的网格和四邻近所包含的面片比例小于阈值时，就将该面片删除。

通过上述步骤，最后得到了致密的三维点云模型，并开始对这两期点云进行数字几何处理。

6)两期点云配准；

将煤场开挖前的点云数据和开挖后的点云数据进行配准，首先采用基于特征点的粗配准方法，通过布设在煤场周围的控制点来得到具体的坐标信息，进行坐标转换，选取同名标识点将前期点云数据转换到统一坐标系下，实现点云的粗配准。在特征点粗配准之后，前后两期点云数据已经统一在同一个坐标系下，此时再采用icp算法进行精配。由于煤场在采挖前后两期点云形态发生变化，影响配准精度，此处利用在两期点云数据保持不变的区域作为基准点，通过icp算法计算全局最优的转换参数，再通过这些转换参数对形态发生变化的区域进行刚性转换，从而实现点云精配准。

7)开挖区域精确定位和点云模型边界确定；

煤场开挖过程中，同一个煤场存在多处作业的情况，对于整个煤场而言，采挖量的体积变化是指多个作业区域煤量变化的综合效应。此处首先对煤场不同采挖区域进行精确定位，从而对不同区域进行体积计算。再通过算法实现两期点云的边界确定，具体为：

煤场开挖过程中，同一个煤场存在多处作业的情况，对于整个煤场而言，采挖量的体积变化是指多个作业区域煤量变化的综合效应。故此处首先需要对煤场不同采挖区域进行精确定位，定位的同时，算法对两期点云重合之外部分形成的区域的边界进行寻找。算法流程如下，设前期点集为a，设p为点集a中的点，p∈a。后期点集为b，q为点集b中的点，q∈b。定义q^*为集合b中到点p距离最近的点。并有q^*＝argmin‖p-q‖，并设定一个阈值，当‖p-q^*‖≥δd时，把该点p存入点集a′中。同时，定义p^*为集合a中到点q距离最近的点。并有p^*＝argmin‖q-p‖，当‖q-p^*‖≥δd时，把该点q存入点集b′中。最终得到的点集a′、b′即为目标点集。7-1)建立点云集合a、b、a′、b′，前期点云形成集合a，后期点云形成集合b，并建立两个新的点云集合a′、b′；

7-2)计算点间最小距离：在点集b中找到点q^*，使得点q^*满足q^*＝argmin‖p-q‖。其中，p为点集a中的点；q、q^*为点集b中的点；

7-3)寻找目标点集a′：在点集a、b中，当满足‖p-q^*‖≥δd时，且q^*＝argmim‖p-q‖，则把该点p存入目标点集a′中。其中，δd为距离约束阈值。

7-4)再次计算点间最小距离：在点集a中找到点p^*，使得点p^*满足p^*＝argmin‖q-p‖。

7-5)寻找目标点集b′：在点集a、b中，当满足‖q-p^*‖≥δd时，且p^*＝argmin‖q-p‖，则把该点q存入目标点集b′中。

7-6)重复步骤(2)、(3)、(4)、(5)：循环点集a、b中所哟点；

7-7)将得到的点集a′、b′进行合并，获取不同采挖区域的点云模型边界，并得到最终点云集。

8)点云网格化；

将点集a′、b′合并之后，此时对合并后的点集从点云到网格进行poisson重建，得到一个精确的封闭的网格模型。泊松表面重建算法具体流程如下：

8-1)首先进行问题离散化，此处通过使用自适应的八叉树来表示隐式函数，并用样本点的位置采用一颗八叉树o，之后给八叉树o中的每一个节点o∈o附加一个函数f0。然后，定义一个基函数f：r³→r，使该函数空间可进行距离平移和尺度缩放。对于每一个节点o∈o，定义函数f0为节点中心o，以o的尺寸展开：

其中，o.c为节点o的中心，o.w为节点o的宽度。

此处函数f定义为一个带通滤波器与它自身的n次卷积：

f(x,y,z)≡(b(x)b(y)b(z))^*n

8-2)其次，定义向量场，为了满足节点的精度要求，采用三相线差值法把样本点分配到八个最邻近节点。

8-3)最后进行泊松问题求解，通过求解函数χ，使得相应方程达到最小值。最后进行等值面提取，选择逼近输入的采样点位置来提取等值面，有采样点的位置估计χ，之后使用平均值来提取等值面，最终得到网格模型。

9)煤场开挖工程量计算；

对封闭的网格进行网格体积计算。网格是由一组顶点和多边形表示。在计算体积之前，需要对模型进行预处理并确保所有多边形都是三角形。三角形的法线方向是通过顶点的顺序和右手定则来确定，对于两个邻近三角形，如果它们的公共边方向相反，则两个三角形的法线方向一致。在三维场景下，基本的计算单元是四面体，对于每个三角形，我们把它们的顶点分别与原点相连并形成一个四面体，如图2所示。在二维情况下，我们定义如下每个单元四面体的带符号体积：数值大小表示该四面体体积，符号是由原点是否在三角形防线方向的同一侧来决定的。如图2所示，三角形acb的法线为nacb，则四面体oacb的体积为：

因为原点o与法线nacb方向相反，所以四面体符号是正的。四面体符号也可以通过内积：来判断。

在实际应用中，我们只需要计算：

因此，煤场任一采挖区域的采挖量体积可表示为：

此处：i表示三角网格或者单元四面体的索引号；(xi1yi1zi1)、(xi2yi2zi2)和(xi3yi3zi3)为三角网格i的顶点坐标，因为有序排列所以三角网格i的法线方向与其他保持一致。三角网格模型的体积总是正值，最终结果可通过取v′total的绝对值来得到。获取网格模型体积之后，通过上述方法，我们可以依次计算出不同采挖区域的网格模型体积，把不同采挖区域网格体积进行求和，最后得到总的网格体积，所求体积即为煤场的采挖量。