CN112734929A - 基于网格细分算法的土石坝复杂土料场开挖方量计算方法 - Google Patents

Abstract

基于网格细分算法的土石坝复杂土料场开挖方量计算方法，包括以下步骤：扫描土料场开挖面的三维点云数据，并且通过点云配准算法实现不同扫描批次的点云的配准，以此获得开挖面原始不规则分布的点云。通过Doo‑Sabin细分算法处理得到的原始不规则分布的点云，使其形成比较规则的开挖面；然后，再进行Kriging插值计算得到参与开挖方量计算的空间点云。通过MATLAB编写网格细分算法，实现开挖方量的计算。本发明一种基于网格细分算法的土石坝复杂土料场开挖方量计算方法，能够对土料场开挖方量进行快速估算，该方法具备工作量小、测量相对安全、效率高。

Description

基于网格细分算法的土石坝复杂土料场开挖方量计算方法

技术领域

本发明属于土石坝施工技术领域，尤其涉及一种基于网格细分算法的土石坝复杂土料场开挖方量计算方法。

背景技术

近年来，心墙堆石坝由于具有基础适应性、能就地取材等诸多优点，在我国得到广泛使用，特别在我国在西南地区，海拔超过200m的特高坝中，心墙堆石坝数量众多，除已建成的糯扎渡、长河坝等工程外，两河口、双江口等工程也正在建设。心墙坝填筑中对土料要求较高，且需求量巨大，往往需要从多个料场开采。通常情况下，土料场形态各异，开采表面和底面形态极不均匀，传统的土方量计量方法主要是通过测量若干个典型断面的面积，再乘以断面间距进行累加得到，这种方法需要对多个断面数据进行测量，对于开挖形态复杂的断面，还需要加密测量断面，后处理工作量较大，且受断面选择的影响，往往存在一定的误差；对于分台阶开挖的料场，台阶较高，传统的测量方法还存在一定的困难和危险性。

发明内容

为解决上述技术问题，本发明提供一种基于网格细分算法的土石坝复杂土料场开挖方量计算方法，能够对土料场开挖方量进行快速估算，该方法具备计算方式简便、工作量小、计算精度可靠、效率高等优点。

本发明采取的技术方案为：

基于网格细分算法的土石坝复杂土料场开挖方量计算方法，包括以下步骤：

步骤1：扫描土料场开挖面的三维点云数据，并且通过点云配准算法实现不同扫描批次的点云的配准，以此获得开挖面原始不规则分布的点云。

步骤2：通过Doo-Sabin细分算法处理步骤1得到的原始不规则分布的点云，使其形成比较规则的开挖面；然后，再进行Kriging插值计算得到参与开挖方量计算的空间点云；

步骤3：通过MATLAB编写网格细分算法，实现开挖方量的计算。

所述步骤1中，点云配准算法用于三维激光扫描仪所扫描点云从相对坐标到大地坐标的转换，属于坐标点云预处理。

所述步骤2中，Doo-Sabin细分算法指的是附图1中点云处理中的部分。针对配准后的点云，即三维激光扫描仪所扫描后经过坐标转换至大地坐标下的点云，是原始的不规则分布的点云，并且基于开挖现场的情况，会存在开挖面局部凹凸不平的情况，从而影响计算精度，因此需要运用Doo-Sabin细分算法对开挖面原始不规则点云进行光滑处理。

所述步骤2中，Kriging插值计算指的是附图1中点云处理中的部分。针对光滑处理后的点云，并不能运用在网格细分算法的开挖方量估算中。由后述网格细分算法可知，参与开挖方量估算的点云是每个网格柱体的四个顶角坐标，因此，本发明引入Kriging插值算法用于网格单元柱体中顶角坐标的插值计算，以此获得用于实际开挖方量估算的计算点云。

所述步骤3中，网格细分算法指的是附图1中方量计算的部分，由上述可知，通过点云处理的过程获得了计算点云，下一步就是开挖方量估算。

本发明一种基于网格细分算法的土石坝复杂土料场开挖方量计算方法，技术效果如下：

1)本发明引入了Doo-Sabin细分算法，用于对原始不规则点云的光滑处理，经过一次细分将点云不仅保留了原始点云数据的规律性，同时保留了开挖面的完成性。将光滑处理后的点云应用到计算点云的插值计算，以此提高计算精度。通过实际验证，经过Doo-Sabin细分算法光滑处理后的点云参与开发方量计算能够使计算时间相应缩短以及计算精度有所提高。

2)本发明引入了Kriging插值计算，用于计算点云的插值计算，将复杂化的原始空间扫描点云经过插值后转化到空间柱体的顶点上，实现点云的简单化。同时，利用Kriging插值算法的空间自协方差最佳插值优点实现原始开挖面点云的插值，充分保留了原始开挖面的空间变化规律，友好的模拟了开挖面在空间内的三维形态变化。

3)本发明所提出的基于网格细分算法的土石坝复杂土料场开挖方量计算方法，通过两个开挖面点云进行相应的处理后即可进行一次开挖规划内的土料开挖方量计算。比起以往的计算方法，具有计算方式简便、计算结果可靠、效率高等优点。

附图说明

图1为本发明的计算流程图。

图2为本发明的实施例的开挖面点云数据图。

图3为本发明的计算基础参数含义示意图。

图4为本发明的边界内点云提取判断示意图。

图5为本发明的开挖方量估算流程示意图。

图6为本发明的开挖方量估算的网格划分示意图。

图7为本发明的开挖料场区域内三维模型图。

图8为本发明的实例计算的点云范围示意图。

图9(1)为本发明的实例计算的三维模型图；

图9(2)为图9(1)的俯视图；

图9(3)为图9(1)的侧视图。

图10(a)为本发明的Doo-Sabin细分算法拓扑规则示意图(原始网格)；

图10(b)为本发明的Doo-Sabin细分算法拓扑规则示意图(一次细分)；

图10(c)为本发明的Doo-Sabin细分算法拓扑规则示意图(F面)；

图10(d)为本发明的Doo-Sabin细分算法拓扑规则示意图(E面)；

图10(e)为本发明的Doo-Sabin细分算法拓扑规则示意图(V面)。

图11为本发明的相对坐标系与大地坐标系的三维基准关系示意图。

图12为本发明的坐标转换计算示意图。

具体实施方式

基于网格细分算法的土石坝复杂土料场开挖方量计算方法，包括以下步骤：

步骤1、点云获取：

如附图1所示的步骤流程示意图，首先需要通过三维激光扫描仪等设备扫描土料场开挖面的三维点云数据，并且通过点云配准算法实现不同扫描批次的点云的配准，以此获得开挖面原始不规则分布的点云，其点云如附图2所示。

所述步骤1的点云配准算法用于三维激光扫描仪所扫描点云从相对坐标到大地坐标的转换。土料场开挖面所扫描的点云数据通常需要多站多视角扫描来获得完整的开挖面点云数据，而各站的点云数据的坐标系统相互独立，后续处理需要将其转换到同一个坐标系下，这个转换过程称为点云配准。

本发明通过扫描数据特征及以往配准模型的研究，建立了基于布尔沙-沃尔夫模型的简易配准算法，其三维空间关系如图11所示：

为相对坐标系与大地坐标系的三维基准关系，xoy为相对坐标系，XOY为大地坐标系。相对坐标系与大地坐标系的Z轴方向必须同为向上，在实际测量中可通过三维激光扫描仪的控制使相对坐标系的Z轴方向向上。在经过坐标转换时需首先绕Z轴旋转，再进行X、Y、Z的平移即可完成坐标的配准。因此其转换模型可表示如下：

式(4)中，△X、△Y、△Z为平移参数，X'、Y'、Z'为相对坐标系点云的相对坐标，R(θ)为绕Z轴的旋转矩阵，表示如下：

如图12所示为坐标转换计算示意图:

xoy为相对坐标系，XOY为大地坐标系，坐标o点为相对原点(即测站点原点)。控制点为两坐标系统的公共点，公共点的在大地坐标系内的值可通过RTK或GPS测得。原点o、控制点1与控制点2的大地坐标分别记为(X₀,Y₀,Z₀)、(X₁,Y₁,Z₁)、(X₂,Y₂,Z₂)，原点、控制点1与控制点2的相对坐标为(0,0,0)、(x₁,y₁,z₁)、(x₂,y₂,z₂)。

则由几何关系可得如下方程式：

其中：△x＝x₁,△y＝y₁。

解上述方程组，求得△X,△Y如下：

分别求解△X、△Y、△Z，则控制点1大地坐标的X、Y、Z分别为：

以此类推，即可将所有扫描数据转换为大地世界坐标。

步骤2：点云处理：

通过Doo-Sabin细分算法处理步骤1得到的原始不规则分布的点云，使其形成比较规则的开挖面；然后，再进行Kriging插值计算得到参与开挖方量计算的空间点云。

所述步骤2的Doo-Sabin细分算法，旨在用于解决原始点云不规则分布导致开挖方量计算精度偏低的问题。通过Doo-Sabin细分算法对原始不规则点云进行一次细分光滑处理，不仅保留了原始点云数据的规律性，同时保留了开挖面的完成性。再将光滑处理后的点云应用到计算点云的插值计算，以此提高计算精度。具体实现如下：

对于初始控制网格M^k，其中任意一个多边形P^k的顶点V_i ^k(0≤i≤n)，细分后对应的多边形为P^k+1，每个顶点V_i ^k生成对应的新顶点V_i ^k+1：

其中：

其拓扑规则如附图10(a)～图10(e)所示，其拓扑步骤如下：

(1)对于初始控制网格M^k的每个顶点V_i ^k采用式(1)生成新的顶点V_i ^k+1；

(2)依次连接控制网格M^k的每个面中的所有新生成的顶点V_i ^k+1得到F面；

(3)依次连接控制网格M^k的每个边两边面所对应的新顶点V_i ^k+1得到E面；

(4)依次连接控制网格M^k的每个点新生成的新顶点V_i ^k+1得到V面。

所述步骤2的Kriging插值算法用于计算点云的插值计算，将复杂化的原始空间扫描点云经过插值后转化到空间柱体的顶点上，实现点云的简单化。Kriging插值算法所具有的空间自协方差最佳插值优点，能够适用于区域化变量存在空间相关性的情况。本发明所研究的开挖面三维点云数据经过预处理后，在空间范围内具有连续变化的规律，因此可以通过Kriging插值算法中的实验变差函数来模拟空间内的三维形态变化。基于此，本发明采用Kriging插值算法通过细分后的点云数据插值出给定平面坐标下的高程值。其具体描述如下：

设在开挖面的三维空间区域内有n个坐标点：

{(x₁,y₁,z₁),(x₂,y₂,z₂),(x₃,y₃,z₃),…,(x_n-1,y_n-1,z_n-1),(x_n,y_n,z_n)}

所组成的数据对,(x_i,y_i)为空间二维坐标，z_i看做样本点属性值，则对于空间内任意一点(x₀,y₀,z₀)，其估计属性值

可由以下线性组合公式进行估计：

式中：

——待求点属性估计值；

λ_i——权重系数，是用空间上所有已知点的数据加权求和来估计未知点的值，但是权重系数并非距离的倒数，而是能够满足点(x₀,y₀)处的估计值

与z₀的差最小的一套最优系数。

因此针对上述计算原理，针对开挖面待求点空间三维坐标的Kriging插值计算的步骤为：

①、对处理后的三维点云两两计算坐标点间距与半方差；

②、通过选用理论模型拟合距离与半方差的函数关系；

③、对于任意一点，自定义搜索半径及样本点数目，通过理论模型计算半方差；

④、通过OK插值算法(普通Kriging插值算法)的求解方程矩阵计算权重系数λ_i；

⑤、将权重系数λ_i代入公式(1)即可求得待求点坐标值。

步骤3：方量计算：

通过MATLAB编写网格细分算法，实现开挖方量的计算。步骤3包括以下两大步骤：

(1)计算点云数据集、以及基础参数的获取：

进行网格细分算法计算开挖方量之前，首先要获取两组计算点云数据集、以及一些基础参数，用于开挖方量计算，步骤如下：

步骤①、针对一次扫描的点云数据，第一步提取出X方向的最大值X_max，最小值X_min，Y方向的最大值Y_max，最小值Y_min；

步骤②、根据X方向的最大值X_max，最小值X_min，Y方向的最大值Y_max，最小值Y_min，确定X方向的最大跨度是D_x＝X_max-X_min，Y方向的最大跨度是D_y＝Y_max-Y_min；

步骤③、设最小插值密度为d_x、d_y，d_x、d_y分别表示建模时X方向、Y方向的网格间距；

步骤④、令

m、n分别表示X方向、Y方向的网格数量，此处的m、n是整数，考虑向下或向上取整即可。

步骤⑤、将三维点云数据按上述的方式，划分为m×n个网格，则每个网格均是d_x×d_y的矩形，得到一组新的平面二维点云数据集(X_i,Y_i)，点云数据集(X_i,Y_i)为m×n个网格的交点集；

步骤⑥、为了获得新的点云数据集(X_i,Y_i)中各点的Z坐标，将新的点云数据集(X_i,Y_i)代入原始点云数据中，通过Kriging插值法，插值计算出新的新的点云数据集(X_i,Y_i)中各点的Z坐标，得到了一组全新且规则的三维点云数据集(X_i,Y_i,Z_i)；

步骤⑦、将三维点云数据集(X_i,Y_i,Z_i)做网格填充处理，即可得到土料场开挖面曲面；

步骤⑧、将二次扫描获得的原始点云数据，同样按照上述步骤①-步骤⑦进行处理，获得第二组三维点云数据集(X'_i,Y'_i,Z'_i)。

经过上述步骤①-步骤⑧获得了两组计算点云数据集(X_i,Y_i,Z_i)、(X'_i,Y'_i,Z'_i)，分别表示一次开挖计划的两个断面点云，将两组计算点云数据集(X_i,Y_i,Z_i)、(X'_i,Y'_i,Z'_i)以及基础参数d_x、d_y、Z_i、Z'_i经过开挖方量计算，即可计算出一次开挖的开挖方量。网格细分算法所涉及的各参数具体表示如附图3所示。

(2)开挖方量计算方法：

由上述计算点云数据集的获取步骤，获得了两组新的计算点云数据集(X_i,Y_i,Z_i)和(X'_i,Y'_i,Z'_i)与基础计算参数d_x、d_y、Z_i、Z'_i。开挖方量计算时只需要再获得料场的开挖边界，即可进行计算，所述料场的开挖边界用闭合曲线L表示。

开挖方量计算方法步骤如下：

1、边界内点云提取：

计算时，需要逐一判断模型中的小网格的单元格子是否位于开挖边界以内，若网格部分位于开挖边界以内，则视为边界内网格，判断的方法是自网格形心点向任意方向引一条射线，判断射线与边界曲线的交点，若交点数为奇数，则该网格位于边界以内，否则，该网格位于边界以外。其具体步骤分述如下：

a：计算点云数据集(X_i,Y_i,Z_i)和(X'_i,Y'_i,Z'_i)与基础参数d_x、d_y、Z_i、Z'_i，以及料场的开挖边界数据获取；

b：在开挖面点云范围内，获得m×n个网格，以及每个网格的型心坐标：

(X_k,Y_k)(k＝1,2,3,…,m×n)；

c：对任意一型心坐标点(X_k,Y_k)，引一条任意方向的射线l，穿过开挖边界线L；其示意图如附图4所示。

d：判断射线l与开挖边界线L的交点个数，若交点个数为奇数，则该型心点位于开挖边界内，否则该型心点位于开挖边界外。

e：针对一次开挖计划的两个开挖断面点云，采用步骤c-步骤d重复执行，直到所有的点云判断完毕，得到最后位于开挖边界内的两组计算点云(x_i,y_i,z_i)和(x'_i,y'_i,z'_i)。

2、开挖方量估算：

经过上述边界内点云提取步骤，获得了两组开挖面内的两组计算点云(x_i,y_i,z_i)和(x'_i,y'_i,z'_i)，并且两组点云的数量是完全相同的，且各自组成的网格必须是一一对应的，即开挖体的计算模型，是由若干个d_x×d_y×h_i的棱体组成的，h_i为第i个棱柱体的高度，计算时，取该棱柱体的侧边线的高度的平均高度，作为棱柱体计算高度，即h_i＝z_i-z'_i，棱柱体的底面和顶面分别是由(x_i,y_i,z_i)和(x'_i,y'_i,z'_i)中一组对应的网格所围成的两个曲面，设计算模型中共有N个棱柱体，则土方量计算的公式为：

计算流程如附图5所示。土方量计算时，通过无限缩小X方向的最小插值密度d_x和Y方向的最小插值密度d_y，可以在一定程度上提高计算的精度，当d_x和d_y充分小时，计算模型中的棱柱体接近规则的平面柱体，减小了由h_i计算时带来的误差。计算示意图如附图6所示。

实施例：

以西南某水电工程一个土料场的开挖土方量计算为例：

该料场位于坝址上游约36km处，料场分布高程2740～2865m，料场前缘拔河高度约85m。料场内分布耕地、民居，属水库淹没区。其分区总面积约19.5m²，勘测设计总方量约681万m³。如附图7所示为该料场区域范围内根据原始地形等高线利用点云数据逆向建模所建立的料场区域内三维模型。其基本原理是通过三维激光扫描仪所获取的点云经过规则三角网格生成算法逆向生成开挖面的三角网，然后根据三角网建设置基准面生成三维模型。附图8所示为实例计算的点云范围示意图，开挖区位于料场的B1区。

根据上述规则三角网格生成算法，运用AutoCAD逆向建立一次开挖计划内开挖体的三维模型如附图9(1)、图9(2)、图9(3)所示。并用AutoCAD的体积查看功能查询得该开挖体的体积为29505.7m³。

以AutoCAD的体积查看功能查询得到的体积为参考，运用MATLAB编程实现本发明基于网格细分算法的复杂土料场开挖方量估算方法并按不同的网格密度进行开挖方量的估算。其计算平台为CPU：Inter(R)Core(TM)i5-4210M 2.60GHz,RAM 16G的Windows10计算机。计算结果如下表所示：

表1计算结果汇总表

由表1的计算结果可知，不同的网格密度会得到不同的计算结果，并且计算结果呈现网格密度越大，计算结果越精确的现象。从运行时间方面来看，网格密度的增加势必会增大程序的运行时间。表1中数据很明显的反映了这个结果。并且网格密度越大，时间达到两小时以上。但是从计算结果可看出，当网格密度大于2000后，相对误差百分比仅有很小的变化，因而可看出，过大的增加网格密度已无实际意义，反而会给工作带来巨大时间的消耗。

Claims (4)

1.基于网格细分算法的土石坝复杂土料场开挖方量计算方法，其特征在于包括以下步骤：

步骤1：扫描土料场开挖面的三维点云数据，并且通过点云配准算法实现不同扫描批次的点云的配准，以此获得开挖面原始不规则分布的点云；

步骤2：通过Doo-Sabin细分算法处理步骤1得到的原始不规则分布的点云，使其形成比较规则的开挖面；然后，再进行Kriging插值计算得到参与开挖方量计算的空间点云；

步骤3：通过MATLAB编写网格细分算法，实现开挖方量的计算。

2.根据权利要求1所述基于网格细分算法的土石坝复杂土料场开挖方量计算方法，其特征在于：所述步骤3包括计算点云数据集、以及基础参数的获取：

进行网格细分算法计算开挖方量之前，首先要获取两组计算点云数据集、以及一些基础参数，用于开挖方量计算，步骤如下：

步骤①、针对一次扫描的点云数据，第一步提取出X方向的最大值X_max，最小值X_min，Y方向的最大值Y_max，最小值Y_min；

步骤②、根据X方向的最大值X_max，最小值X_min，Y方向的最大值Y_max，最小值Y_min，确定X方向的最大跨度是D_x＝X_max-X_min，Y方向的最大跨度是D_y＝Y_max-Y_min；

步骤③、设最小插值密度为d_x、d_y，d_x、d_y分别表示建模时X方向、Y方向的网格间距；

步骤④、令

m、n分别表示X方向、Y方向的网格数量，此处的m、n是整数，考虑向下或向上取整即可；

步骤⑤、将三维点云数据按上述的方式，划分为m×n个网格，则每个网格均是d_x×d_y的矩形，得到一组新的平面二维点云数据集(X_i,Y_i)，点云数据集(X_i,Y_i)为m×n个网格的交点集；

步骤⑥、为了获得新的点云数据集(X_i,Y_i)中各点的Z坐标，将新的点云数据集(X_i,Y_i)代入原始点云数据中，通过Kriging插值法，插值计算出新的新的点云数据集(X_i,Y_i)中各点的Z坐标，得到了一组全新且规则的三维点云数据集(X_i,Y_i,Z_i)；

步骤⑦、将三维点云数据集(X_i,Y_i,Z_i)做网格填充处理，即可得到土料场开挖面曲面；

步骤⑧、将二次扫描获得的原始点云数据，同样按照上述步骤①-步骤⑦进行处理，获得第二组三维点云数据集(X'_i,Y_i’,Z'_i)；

经过上述步骤①-步骤⑧获得了两组计算点云数据集(X_i,Y_i,Z_i)、(X’_i,Y_i',Z’_i)，分别表示一次开挖计划的两个断面点云，将两组计算点云数据集(X_i,Y_i,Z_i)、(X’_i,Y_i',Z’_i)以及基础参数d_x、d_y、Z_i、Z’_i经过开挖方量计算，即可计算出一次开挖的开挖方量。

3.根据权利要求2所述基于网格细分算法的土石坝复杂土料场开挖方量计算方法，其特征在于：基于两组计算点云数据集(X_i,Y_i,Z_i)、(X’_i,Y_i',Z’_i)以及基础参数d_x、d_y、Z_i、Z’_i，开挖方量计算时只需要再获得料场的开挖边界，即可进行计算，所述料场的开挖边界用闭合曲线L表示。

4.根据权利要求2或3所述基于网格细分算法的土石坝复杂土料场开挖方量计算方法，其特征在于：所述步骤3中，开挖方量计算方法步骤如下：

1、边界内点云提取，其具体步骤分述如下：

a：计算点云数据集(X_i,Y_i,Z_i)和(X’_i,Y’_i,Z’_i)与基础参数d_x、d_y、Z_i、Z’_i，以及料场的开挖边界数据获取；

b：在开挖面点云范围内，获得m×n个网格，以及每个网格的型心坐标：

(X_k,Y_k)(k＝1,2,3,…,m×n)；

c：对任意一型心坐标点(X_k,Y_k)，引一条任意方向的射线l，穿过开挖边界线L；

d：判断射线l与开挖边界线L的交点个数，若交点个数为奇数，则该型心点位于开挖边界内，否则该型心点位于开挖边界外；

e：针对一次开挖计划的两个开挖断面点云，采用步骤c-步骤d重复执行，直到所有的点云判断完毕，得到最后位于开挖边界内的两组计算点云(x_i,y_i,z_i)和(x’_i,y’_i,z’_i)；

2、开挖方量计算：

经过上述边界内点云提取步骤，获得了两组开挖面内的两组计算点云(x_i,y_i,z_i)和(x’_i,y’_i,z’_i)，并且两组点云的数量是完全相同的，且各自组成的网格必须是一一对应的，即开挖体的计算模型，是由若干个d_x×d_y×h_i的棱体组成的，h_i为第i个棱柱体的高度，计算时，取该棱柱体的侧边线的高度的平均高度，作为棱柱体计算高度，即h_i＝z_i-z’_i，棱柱体的底面和顶面分别是由(x_i,y_i,z_i)和(x’_i,y’_i,z’_i)中一组对应的网格所围成的两个曲面，设计算模型中共有N个棱柱体，则土方量计算的公式为：

土方量计算时，通过无限缩小X方向的最小插值密度d_x和Y方向的最小插值密度d_y，可以在一定程度上提高计算的精度；当d_x和d_y充分小时，计算模型中的棱柱体接近规则的平面柱体，减小了由h_i计算时带来的误差。

Images