基于无人机倾斜摄影的土石方工程量计量方法

2022-05-20覃毅媚

华东交通大学学报 2022年2期

杨轸，覃毅媚

（同济大学道路与交通工程教育部重点实验室，上海 201804）

工程建设中土石方是一个主体工程，无论是道路路基，还是房屋建筑尤其是公共建筑的施工，都离不开土石方的开挖与填筑，且具有工程量大、工期长等特点。土石方工程量的计算直接影响到工程造价和工程进度，以及施工方案的最优选择。传统的土石方工程量计算主要是通过单点测量进行离散点数据的外业采集后，建立模型进行计算，如等高线法、断面法和方格网法等[1-2]。但是由于地形条件变化复杂、土石方外形不规则，这些传统方法外业工作量大、数据采集困难，且很难准确地表征地貌真实的形态，对土石方的计算精度难以保证[3]。

随着数字地面模型技术的成熟和计算机软硬件的快速发展，国内外许多学者都对土石方量的三维计算进行了研究[4-9]。程建川等提出基于不规则三角网形式数字高程模型的体积计算方法[4]；Eliseev等使用Autodesk Civil 3D 建立地形表面模型用于土石方工程量的计算[8]；Tanoli 等基于参数化建模技术，通过比较两个三角网曲面来计算土方工程量[9]。

土石方工程量的三维计算，主要通过原始地面模型和施工后模型或设计模型的叠加运算来完成，地形模型的准确表达直接决定其计算精度[10]。在土石方工程量计算方面，目前建立三维模型的数据获取方式包括基于水准仪的断面高程测量，基于全站仪的特征点位测量，现有地图数字化，基于差分GPS的散点测量等[11-12]，这些方法通常难以同时兼顾效率和成本之间的关系。近年来发展起来的无人机航测[13]，尤其是倾斜摄影技术，具有操作灵活、数据采集快、成本低、分辨率高、数据呈现性强等优势，同时突破了传统航测从垂直角度拍摄的局限，能有效弥补以上方法的不足，成为快速获取地形全方位三维数据的有效方法[14]，可为施工过程中土石方量计算提供数据采集的手段[15-18]。

基于倾斜摄影进行土石方工程量计算需解决的关键问题有：①点云简化。倾斜摄影采集的地形点云数据具有超高密度的特点，掺杂着大量的冗余信息，极大地增加了地形表面重建的时间和空间成本，在地形关键特征信息不丢失的前提下，需要对原始点云数据进行简化。②表面重建。土石方量计算的准确性与重建的地形表面模型密切相关，为了保证体积计算的精度，地形表面重建结果不应出现缺损、空洞，且应保留更多的细节纹理特征，呈现良好的表面特性，需要选择合适的重建方法及参数。③体积计算。三维重建生成的施工前后两期地形模型为开放表面，不便直接进行体积计算，如何将其转换为封闭体模型求体积值是需要解决的问题。

针对上述问题，提出一种快速准确计算道路土石方工程量的具体方法，进行点云简化、表面重建、构造包围盒、布尔运算，最后计算体积差异量，使得开展土石方的挖填分析与计算直观有效。

1 数据采集

1.1 外业数据采集

1）航线设计。现场勘查确定测区范围后，进行航线设计，包括飞行姿态、飞行高度、飞行速度、重叠度、地面分辨率等飞行参数。需要根据低空数字航空摄影测量规范的相关规定，并结合测区的具体情况进行航线的设计，确保影像的重叠度和分辨率等符合作业要求。

2）像控点布设。像控点是空中三角测量过程中的关键内容，其精度和数量直接影响模型精度。像控点的数量和位置根据实际所需的精度和测量区域的大小均匀布设，并且通常设置在相对平坦的位置，在明显的地物拐角点或利用标靶板布设像控点，方便点位的查找和定位。

3）数据检查。无人机航摄飞行结束后，对外业影像数据和像控点的质量进行检查，如影像清晰度、像控点精度等，以便对不合格照片的区域进行及时补飞，对不满足精度要求的像控点进行重测。

1.2 内业数据处理

通过无人机搭载相机进行航摄飞行，倾斜摄影测量获得多角度的影像数据及其高度、姿态、经纬度等POS 信息，可通过相应的软件平台，导入外业像控点数据文件，进行空中三角测量。空中三角测量利用外业采集的影像、POS、像控点数据，解算出每张影像的位置和空中姿态即外方位元素，再结合多视影像密集匹配，可自动生成高密度三维点云模型。

2 点云简化

点云简化的核心思想就是在进行尽可能最大限度地简化处理之后还能表达三维重建所需的特征信息，并且有着较高的简化速度，由此选择合适的简化方法和简化率。

2.1 评价方法

目前，点云简化效果的评价一般是通过对重建表面进行视觉对比分析来实现的，主观性较大，当精度比较接近、肉眼难以判断时，就无法得出结论[19]。由于地形点云为面状点云，不能构成封闭的三维体模型，不适用于直接计算体积变化比的形状误差法[20]。考虑输入点云和基于其重建的表面模型位于同一个坐标系的相同位置，使用重建误差（reconstruction error）[21]评价法，通过输入点距重建表面的平均距离来评价重建结果的质量，重建误差越小，表示简化及重建效果越好。

2.2 简化方法

地形点云具有数据量庞大、纹理特征复杂的特点，在进行如图1 所示的简化效果的对比后，最终选取了对特征保留及表面重建更有利的加权局部最优投影（weighted locally optimal projection，WLOP）方法。可以看出，相较其他两种方法，WLOP 方法在保留必要特征信息的同时，进行了一定的平滑处理，以减少表面重建过程中空洞现象的发生。

图1 不同简化方法的效果对比Fig.1 Results using different simplified methods

WLOP 方法[22]是一种基于投影的简化算法，为每个点增加一个局部密度权重，然后对输入点云数据的子集进行不断迭代的投影，输出一个更依附于底层形状的新的点集，达到简化的目的，改善点云的不均匀采样和表面光顺度。对于一个原始点集P={pj}j∈J⊂R3，j=1，2， …，pj为P 中的每一个点。WLOP 尝试通过固定点迭代寻找一系列投影点X={xi}i∈I⊂R3，i=1，2，…，将X 投影到P 上，优化能量方程。如果将第k 次迭代中生成的投影点集记为Xk，k=0，1，…，则下一次迭代所得Xk+1需要最小化能量函数

2.3 最优简化率

不同的点云数据，其最优简化率不同。原始基坑开挖测区为200 m×200 m，约500 万个点，进行实验获取最优解的时间成本过高。从原始的测区点云数据中提取一小部分来进行试验，得出最优简化率，再推广到完整点云数据进行处理，以减少工作量。选取32 m×30 m 的区域，约8.5 万个点，包含起伏的地形及具有尖锐特征的地物，曲率跨度较大，可作为典型区域进行试验。

将不同简化率简化后的点云重建表面模型置于图2 中。图3 显示了不同简化率的简化时间和重建误差变化情况，结合图2 中重建模型的视觉纹理对比，可较为明显地看出，本组点云数据的最佳简化率为90%，既能保证一定的重建精度，又能减少时间成本。

图2 不同简化率下的点云重建效果对比Fig.2 Point cloud reconstruction results at different simplification rates

图3 不同简化率的简化时间和重建误差变化Fig.3 Simplification time and reconstruction error at different simplification rates

3 地形表面重建

3.1 重建方法

点云数据的三维重建主要有两种算法：基于组合结构和基于隐函数。针对地形点云数据生成的是开放而非封闭表面、大范围大面积的特点，较为适合的是基于组合结构的尺度空间表面重建和基于隐函数的泊松表面重建。

图4 对比了尺度空间重建和泊松表面重建方法的效果。可以看到，两种重建算法都没有空洞存在，可以进行体积计算；尺度空间算法虽大体上保留了地形地物的特征，但是在特征比较尖锐明显的前提下，而在平坦地面以及细节特征处未能体现出原始的纹理轮廓；相比之下，泊松重建采取隐函数拟合逼近的方式，重建出的模型是水密性的，对噪声点有很好的抵抗性，并且具有良好的表面和细节特性，这与地形表面点云的连续性是相适应的，对后续土石方体积的计算非常有利，因此选择泊松重建来完成倾斜摄影点云数据的地形表面重建。

图4 尺度空间重建和泊松重建效果对比Fig.4 Reconstruction results using different reconstruction methods

3.2 最大树深度

泊松表面重建算法[23]是一种隐函数方法。该算法将具有法向量的点云的表面重建转化为空间泊松问题，也就是泊松方程的求解过程。为了解泊松方程，要将其离散化，最简单的方法是用自适应的八叉树对隐函数进行表达。八叉树（Octree）是一种用于划分三维空间的树状数据结构，其根节点是一个封装了所有点的立方体包围框，然后以循环递归的方式将根节点进行立方体剖分，每次划分成8 个相同的子立方体，直至划分的子立方体达到既定的采样密度或最大递归深度时停止。当剖分深度为n时，该空间被划分为2n×2n×2n份。

八叉树划分的深度与重建模型表面网格的分辨率直接相关，树深度越大，重建表面模型的精度越高。取2.3 节同样部分点云进行实验对比与分析，寻找稳定的参数范围。如图5 所示，当深度大于8之后，重建细节的变化程度肉眼几乎无法分辨。图6和图7 显示了顶点面片数、表面重建误差和重建时间随不同深度的变化情况。

图5 不同重建深度的细节特征Fig.5 Detail features at different reconstruction depths

图6 不同深度的重建耗时和重建误差Fig.6 Reconstruction time and reconstruction error at different depths

图7 不同深度的重建耗时和顶点/面片数Fig.7 Reconstruction time and number of vertices/faces at different depths

可以看出在深度达到7 或8 之后，重建耗时仍在快速增长，但顶点面片的数量和重建误差达到稳定值。这是由于地形点云的分布具有非均匀性，而泊松重建采用的是自适应点云密度划分，就算还没达到最大树深度，对于点云密度较低的稀疏范围也不再进行剖分，以保证表面重建的合理性。这说明当深度达到8 的时候，这部分的点云数据已经剖分到最低采样密度了，已经足够符合我们的精度要求。

然而，把这个结论直接扩展到完整的点云，即采用深度为8 的泊松重建对完整地形点云数据进行重建时，重建表面明显丢失了地形地物的特征信息，图5 的细节特征在其上几乎没有显示。从原理上分析，八叉树是对空间进行划分，划分到同样的采样密度时，对于不同大小的空间，其划分深度应该不一样。根据八叉树的特点，只要求出封装了所有点的最小立体包围框的体积，就可以反求出原始点云的最大树深度。

CGAL（计算几何算法库）可求点集的最小包围球半径，根据最大树深度与空间划分次数的关系，得出完整测区点云的最大树深度计算公式

式中：d 为部分点云的最大树深度；D 为完整点云的最大树深度；r 为部分点云的最小包围球半径；R 为完整点云的最小包围球半径。

4 土石方工程量计量

获得施工前后的地形表面重建模型后，理论上土石方量可以由施工前后两个表面模型进行布尔运算得到的封闭体模型来计算。但是现有算法中，直接对两个开放的不规则网格表面模型做布尔运算非常不稳定[24]。考虑相对稳定地面与体的布尔运算，利用不规则的地形表面去切割一个规则的封闭体模型（即计算区域的包围盒），提取出土石方量的计算范围，完成面模型到体模型的转换，再进行体积计算。

4.1 包围盒

确定挖方或填方的计算范围，获得计算范围的顶点坐标后，将其构造成一个封闭的平行六面体模型。将计算范围内的土石方区域完整地包围在其中，称其为计算区域的包围盒，并利用右手定则手动设置其面片索引列表，构成一个OFF 格式的三维网格模型，以便布尔运算的进行。

4.2 布尔运算

泊松表面重建后获取的三维地形表面模型是以三角网格模型的形式进行存储的。将包围盒这样一个封闭的网格模型分别与施工前后的开放地形表面模型做布尔运算，得到两个以土石方计算区域为上表面的封闭体模型，分别计算两个体模型的体积，求差即可得到相应填挖量。

4.3 体积计算

目前流行的任意网格模型的体积计算方法中，拟蒙特卡罗法和切都只是对模型体积的估算，而投影法和行列式法虽然能比较准确地计算体积，但是其算法较为繁琐，对于编程的实现不够简单直观[25]。为了能够利用编程更快速准确地计算土石方量，根据三角网格模型的顶点坐标进行体积计算，进行算法的编写，实现土石方量的计算。

式中：i 为三角形的下标；（xi1，yi1，zi1），（xi2，yi2，zi2），（xi3，yi3，zi3）是三角形i 的三个顶点的坐标，它们的有序使得三角形i 的法线与其它三角形的法线方向一致。 Vsum的绝对值即为三维网格模型的体积。

5 工程实例应用与分析

5.1 前期准备

选取的基坑开挖区域整体范围为200 m×200 m。本文算法基于计算几何算法库（CGAL）和点云库（point cloud library，PCL）作为主要支持库，在VS2017 编译环境下由C++语言实现。

本次外业采集使用的无人机是上海瞰融信息技术发展有限公司的大疆精灵Phantom4 RTK，进行倾斜摄影测量后得到密集点云数据，施工前后的点云数量分别为5 071 981 个和4 807 052 个。

5.2 地形表面重建

以90%的简化率对点云进行简化后，进行泊松表面重建。 3.2 节r=16.459 m，其最大树深度d 取为8，而R=100.12 m，则由式（2）计算完整测区点云需要的泊松重建八叉树最大树深度得10.604 78 m，即当最大树深度取11 时，测区的完整点云能以最短的耗时达到最高的精度，重建结果如图8 所示。可以看出，中间区域施工前为堆土，施工后变为平坦地面，其挖方量即为本次工程实例的计算目标。

图8 施工前后重建结果对比图Fig.8 Reconstruction results before and after construction

放大局部细节如图9 所示，从视觉纹理上来说，深度为11 的泊松重建已经达到我们想要的视觉效果。从数学精度方面，首先取施工后的26 个测区野外像控点数据，计算其到施工后重建表面的平均距离即重建误差为0.06 m，对于200 m×200 m 的测区面积来说，6 cm 的重建误差是较低的。同时对三维模型进行精度校核，取10 个特征点的实测三维坐标与模型同名点三维坐标进行比较，结果如表1 所示（平面中误差为0.072 m，高程中误差为0.051 m）。由《三维地理信息模型数据产品规范》可知，采用该方法可满足1∶500 比例的三维精度要求，验证了式（2）的合理性，为土石方量的计算提供了良好的基础。

表1 施工后重建模型精度统计Tab.1 Accuracy statistics of the reconstructed model after construction m

图9 原始影像和泊松重建效果对比Fig.9 Original image and Poisson reconstruction result

5.3 土石方工程量计算

在施工后的正射影像中确定挖方的计算范围并构建包围盒如图10 所示。

图10 土石方计算范围及其包围盒Fig.10 Calculation range and bounding box

使用包围盒分别与施工前后的地形表面模型进行布尔运算，以地形表面对包围盒做剪切，布尔运算结果及其体积如图11 所示。两个体积值相减即为所求土石方挖方量：7 191 m3。

图11 布尔运算及体积计算结果Fig.11 Boolean operation and volume calculation results

图12 显示了从表面重建到体积计算流程总耗时和土石方量计算值随泊松重建不同最大树深度的变化情况。可以看出，当最大树深度到达10、11时，土石方计算量开始收敛，深度再增大时，耗时仍随之增加，但土石方量的计算结果处于一个稳定值，结合4.2 节D=10.604 78，验证了重建深度与体积计算结果之间的关系，以及论文提出的先提取部分点云进行重建分析，再将其参数变换推广到整体点云进行重建与计算的可行性。

图12 流程总耗时和土石方量计算值随最大树深度的变化情况Fig.12 Processing time and earthwork volume at different reconstruction depths

5.4 结果分析

5.2 节中重建地形表面的精度已经能够满足土石方工程量计算的要求，为了验证论文土石方计算结果的准确性，在重建模型上取点并利用传统的方格网法进行对比。方格网法将计算范围划分为多个方格，确定每个方格顶点的填挖高度，并根据填挖高度的加权平均值和方格面积计算土石方量。理论上，方格顶点间距越小，计算精度越高，但相应的计算量也越大。工程上采用全站仪进行外业高程点测量时，一般间距为10～40 m。为了得到更精确的结果进行对比，利用更高精度的3 m 方格网法，在4 652.3 m2的计算范围内，共取516 个点，计算结果与论文提出的计算方法偏差在5 %以内，满足实际工程现场中的土石方计算误差控制[26]，可用于工程项目土石方量计算。

以上分析结果显示，基于无人机倾斜摄影重建的地形表面满足1∶500 大比例尺三维地形模型的精度要求，同时计算得到的土石方工程量精度也达到一般实际工程的测量要求，验证了本文方法的可行性。对比传统方法，基于无人机倾斜摄影的土石方工程量计量方法一方面在外业测量时操作更灵活，不受地形限制，效率得到极大提升；另一方面能获取影像数据，具有可视化效果，整个施工过程中不同阶段的场地状态可三维动态呈现，并能精准地计算各个阶段的土石方工程量，解决实际工程建设中，审计人员往往无法参与项目实施的全过程，在工程竣工后无法再对土石方工程量与造价的真实性进行实测和复核的难题。