CN109870106A - 一种基于无人机图片的建筑物体积测量方法 - Google Patents

Abstract

一种基于无人机图片的建筑物体积测量方法，利用无人机对建筑物进行环绕的图像信息采集，获取建筑物多幅图片序列；经双边滤波处理图像序列后，采用Harris角点检测算法获取图像序列的角点；先采用零均值互相关法进行角点粗匹配，建立一对多的对应关系；再采用松弛迭代法处理粗匹配角点，得到一对一精匹配的角点对；利用从运动中恢复结构算法，循环计算匹配点的空间坐标，并利用光束平差法对空间坐标进行优化，得到建筑物的三维点云；利用逐点插入法对三维点云进行Delaunay三角化处理；最后对构建的Delaunay三角网模型利用积分法计算建筑物的体积。

Description

一种基于无人机图片的建筑物体积测量方法

技术领域

本发明属于图像处理的领域，具体涉及一种基于无人机图片的建筑物体积测量方法。

背景技术

在工程领域中经常需要对建筑物体积进行测量预估，此类测量对象通常具有体积大、形状不规则等特点。因此直接对其进行体积测量难度较大，并且如果仅仅依靠人工估算，不但工作量大，而且测量速度慢，消耗时间长，甚至还可能带来误差大、不安全等一系列问题。

发明内容

本发明目的在于直接利用无人机获取的序列图像信息处理计算得到建筑物体积为工程建设提供手段，提出一种基于无人机图片的建筑物体积测量方法。

一种基于无人机图片的建筑物体积测量方法，包括以下步骤：

步骤1：利用无人机对建筑物进行环绕的图像信息采集，获取建筑物多幅图片序列；

步骤2：经双边滤波处理图像序列后，采用Harris角点检测算法获取图像序列的角点；

步骤3：先采用零均值互相关法进行角点粗匹配，建立一对多的对应关系；再采用松弛迭代法处理粗匹配角点，得到一对一精匹配的角点对；

步骤4：利用从运动中恢复结构算法，循环计算匹配点的空间坐标，并利用光束平差法对空间坐标进行优化，得到建筑物的三维点云；

步骤5：利用逐点插入法对三维点云进行Delaunay三角化处理；

步骤6：最后对构建的DelaunaV三角网模型利用积分法计算建筑物的体积。

进一步地，所述步骤2中，具体包括如下步骤：

步骤2-1：读取图像序列，采用双边滤波做图像去噪声，并对图像做灰度化处理，处理后的图像序列记为I_l(l＝1.2...，k)；

步骤2-2：计算图像I_l上每个像素点水平方向和垂直方向上的一阶导，即梯度I_x和I_y；

步骤2-3：设置矩阵M，

步骤2-4：平滑M得到M₁，方法如下：

式中为高斯模板；

步骤2-5：判断像素点是否为角点，具体地，首先计算响应函数W＝det M₁-k₁(traceM₁)²；其中，k₁取经验数，取值范围为0.04～0.06，其选值影响函数曲线形状的直观性，为便于计算设置为0.04；然后设置阈值T，其大小为所有像素点响应函数最大值的0.01倍，为便于计算设置为0.01，当W≥T时，该像素点则为角点；

步骤2-6：标记角点位置，输出角点。

进一步地，所述步骤3中，具体包括如下步骤：

步骤3-1：用零均值归一化互相关方法，依据角点的灰度值信息粗匹配角点，方法如下：

步骤3-1-1：获取图像I₁中的角点p(x，y)，其相同位置在图像I₂的搜寻区域尺寸大小(2l_x+1)(2l_y+1)下的角点为q(x′，y′)，l_x和l_y分别是的图像宽度和高度；

步骤3-1-2：获取图像中角点对应的灰度值H_l(x，y)；

步骤3-1-3：计算两个角点的互相关性：

其中，u_l(x，y)和σ_l(x，y)表示相关模板(2m+1)(2n+1)内灰度平均值和标准偏差，m＝n＝3表示角点局部窗口大小；

步骤3-1-4：根据两个角点局部区域内灰度值越相近，互相关值越近似为1，可以作为候选匹配点的要求，将阈值T₁设置为0.9，互相关值S_pq＞T₁的匹配点视为粗匹配点对，记为点集A；

步骤3-2：用松弛迭代法，依据匹配强度精匹配角点，其流程步骤如下：

步骤3-2-1：计算粗匹配点对的匹配强度：

其中，(p，q)和(f，g)表示粗匹配点对，N(p)和N(q)分别是以点p和q为中心、X为半径的邻域，f∈N(p)且g∈N(q)，δ表示相似度贡献度，d表示粗匹配点对平均距离，r表示粗匹配点对的相对距离偏差，r非常大时可以忽略匹配点(f，g)，ε_r表示相对距离偏差的阈值T₂，T₂的取值影响匹配的精度，设置T₂为0.3；

步骤3-2-2：判断粗匹配角点的唯一性，具体地，首先按照匹配强度大小从大到小对匹配点对排序，然后选取匹配强度最大和第二大的点对D₁和D₂，并计算选取匹配强度在D₁与D_p之间的点对为精匹配点对；

步骤3-2-3：删除点集A中正确的角点，从步骤3-2-1开始迭代剩下的点直到点集A为空；

步骤3-3：输出精匹配对的坐标点集合，记为点集B。

进一步地，所述步骤4中，具体包括如下步骤：

步骤4-1：通过投影关系计算建筑物图像上匹配点的空间坐标，步骤如下：

步骤4-1-1：对相机标定，获取相机的焦距f_d；

步骤4-1-2：计算基础矩阵F，具体地，首先获取点集B中图像I₁和I₂上的匹配点对m(x，y)和m′(x′，y′)，然后记m(x，y)在x，y轴方向上的距离尺寸为z_x和z_y，其对应像素坐标为(u，v)，记带入m和m′的齐次坐标计算m′^TFm＝0，得到基础矩阵F；

步骤4-1-3：计算本质矩阵E＝P_N ^TFP_N，其中

步骤4-1-4：对本质矩阵E奇异值分解得到旋转矩阵R和平移矩阵t；

步骤4-1-5：计算建筑物图像上匹配点的空间坐标M_s(X_i，Y_i，Z_i)，具体地，首先将M_s(X_i，Y_i，Z_i)投影到图像坐标中为q(x，y)，然后计算M_s＝P^-1q，其中P＝P_l[R|t]，P_l＝P_NP_M，

步骤4-1-6：循环加入图像I_l(l＝3，...，k)，直到所有图像加入完毕；

步骤4-2：用光束平差方法优化空间坐标M_s(X_i，Y_i，Z_i)，方法如下：

步骤4-2-1：计算空间坐标M_s(X_i，Y_i，Z_i)的重投影误差K，其由每幅图像所有距离误差的平方和相加，具体地，首先计算p_il，q_il和r_il：

然后计算

其中，I_il为可见标志值，第i点在第l幅图像中有像素坐标时，I_il＝1；第i点在第l幅图像中没有像素坐标时，I_il＝0；(x_il，y_il)是第i点在第l幅图像中的像素坐标，i＝1，...，N；

步骤4-2-2：依据非线性最小二乘法修正空间坐标M_s(X_i，Y_i，Z_i)以降低步骤4-2-1中的重投影误差K，直到误差函数K的值达到最小为止；

步骤4-3：输出误差函数K最小值对应的空间坐标M_s(X_i，Y_i，Z_i)的集合，记为点集G。

进一步地，所述步骤5中，具体包括如下步骤：

步骤5-1：读取建筑物的空间坐标点集G，构造一个包含所有数据点的超三角形C，将该超三角形C放入Delaunay三角形链表中，Delaunay三角形链表中任一三角形满足任意四点不共圆以及最小内角最大的要求，链表最终生成建筑物的网格模型；

步骤5-2：记超三角形C内某一正在处理的数据点为Q，将点Q插入Delaunay三角形链表中，步骤如下：

步骤5-2-1：在Delaunay三角形链表中找出其外接圆包含点Q的三角形，称为点Q的影响三角形；

步骤5-2-2：删除点Q所有影响三角形之间的公共边；

步骤5-2-3：连接点Q与所有影响三角形的顶点，依据任意四点不共圆和最小内角最大的要求进行筛选，将符合Delaunay三角形要求的新三角形放入Delaunay三角形链表中；

步骤5-3：循环加入超三角形C中新的数据点到Delaunay三角形链表中，执行步骤5-2，直到超三角形C中数据点全部加入Delaunay三角形链表为止，生成完整的建筑物点云三角网模型。

进一步地，所述步骤6中，具体包括如下步骤：

步骤6-1：读取步骤5生成的三角网模型；

步骤6-2：将模型顶部的投影分割成n个正方形，通过投影形成n个长方体，高为h_i；

步骤6-3：任取1个小长方体，假设采样间隔为d_l，其数值等于投影面的边长dx，即d_l＝dx＝dy；

步骤6-4：计算建筑物的体积

本发明的有益效果为：通过操作无人机，解决了当前靠人工估算而导致的工作量大，而且测量速度慢，消耗时间长，甚至还可能带来误差大、不安全等一系列问题。本发明使得建筑物体积测量更加快速方便，减少工作量和工作时间，同时本发明的测量方法精度高，也更加安全。

附图说明

图1为一种基于无人机图片的建筑物体积测量方法流程图。

图2为本发明所述建筑物图像Harris角点提取流程示意图。

图3为本发明所述角点粗匹配流程示意图。

图4为本发明所述角点精匹配流程示意图。

图5为本发明所述空间坐标计算流程示意图。

图6为本发明所述三维点云结果图。

图7为本发明所述Delaunay三角网流程示意图。

图8为本发明所述Delaunay三角网结构图。

图9为本发明所述体积计算投影示意图。

具体实施方式

下面结合说明书附图对本发明的技术方案做进一步的详细说明。

一种基于无人机图片的建筑物体积测量方法，包括以下步骤：

步骤1：利用无人机对建筑物进行环绕的图像信息采集，获取建筑物多幅图片序列。

步骤2：经双边滤波处理图像序列后，采用Harris角点检测算法获取图像序列的角点，流程如图2所示。

所述步骤2中，具体包括如下步骤：

步骤2-1：读取图像序列，采用双边滤波做图像去噪声，并对图像做灰度化处理，处理后的图像序列记为I_l(l＝1.2...，k)。

步骤2-2：计算图像I_l上每个像素点水平方向和垂直方向上的一阶导，即梯度I_x和I_y。

步骤2-3：设置矩阵M，

步骤2-4：平滑M得到M₁，方法如下：

式中为高斯模板。

步骤2-5：判断像素点是否为角点，具体地，首先计算响应函数W＝det M₁-k₁(traceM₁)²；其中，k₁取经验数，取值范围为0.04～0.06，其选值影响函数曲线形状的直观性，为便于计算设置为0.04；然后设置阈值T，其大小为所有像素点响应函数最大值的0.01倍，为便于计算设置为0.01，当W≥T时，该像素点则为角点。

步骤2-6：标记角点位置，输出角点。

步骤3：先采用零均值互相关法进行角点粗匹配，建立一对多的对应关系；再采用松弛迭代法处理粗匹配角点，得到一对一精匹配的角点对。

所述步骤3中，具体包括如下步骤：

步骤3-1：用零均值归一化互相关方法，依据角点的灰度值信息粗匹配角点，流程如图3所示，方法如下：

步骤3-1-1：获取图像I₁中的角点p(x，y)，其相同位置在图像I₂的搜寻区域尺寸大小(2l_x+1)(2l_y+1)下的角点为q(x′，y′)，l_x和l_y分别是的图像宽度和高度。

步骤3-1-2：获取图像中角点对应的灰度值H_l(x，y)。

步骤3-1-3：计算两个角点的互相关性：

其中，u_l(x，y)和σ_l(x，y)表示相关模板(2m+1)(2n+1)内灰度平均值和标准偏差，m＝n＝3表示角点局部窗口大小。

步骤3-1-4：根据两个角点局部区域内灰度值越相近，互相关值越近似为1，可以作为候选匹配点的要求，将阈值T₁设置为0.9，互相关值S_pq＞T₁的匹配点视为粗匹配点对，记为点集J。

步骤3-2：用松弛迭代法，依据匹配强度精匹配角点，流程如图4所示，步骤如下：

步骤3-2-1：计算粗匹配点对的匹配强度：

其中，(p，q)和(f，g)表示粗匹配点对，N(p)和N(q)分别是以点p和q为中心、X为半径的邻域，f∈N(p)且g∈N(q)，δ表示相似度贡献度，d表示粗匹配点对平均距离，r表示粗匹配点对的相对距离偏差，r非常大时可以忽略匹配点(f，g)，ε_r表示相对距离偏差的阈值T₂，T₂的取值影响匹配的精度，设置T₂为0.3。

步骤3-2-2：判断粗匹配角点的唯一性，具体地，首先按照匹配强度大小从大到小对匹配点对排序，然后选取匹配强度最大和第二大的点对D₁和D₂，并计算选取匹配强度在D₁与D_p之间的点对为精匹配点对。

步骤3-2-3：删除点集A中正确的角点，从步骤3-2-1开始迭代剩下的点直到点集A为空。

步骤3-3：输出精匹配对的坐标点集合，记为点集B。

步骤4：利用从运动中恢复结构算法，循环计算匹配点的空间坐标，并利用光束平差法对空间坐标进行优化，得到建筑物的三维点云，流程如图5所示。

所述步骤4中，具体包括如下步骤：

步骤4-1：通过投影关系计算建筑物图像上匹配点的空间坐标，步骤如下：

步骤4-1-1：对相机标定，获取相机的焦距f_d。

步骤4-1-2：计算基础矩阵F，具体地，首先获取点集B中图像I₁和I₂上的匹配点对m(x，y)和m′(x′，y′)，然后记m(x，y)在x，y轴方向上的距离尺寸为z_x和z_y，其对应像素坐标为(u，v)，记带入m和m′的齐次坐标计算m′^TFm＝0，得到基础矩阵F。

步骤4-1-3：计算本质矩阵E＝P_N ^TFP_N，其中

步骤4-1-4：对本质矩阵E奇异值分解得到旋转矩阵R和平移矩阵t。

步骤4-1-5：计算建筑物图像上匹配点的空间坐标M_s(X_i，Y_i，Z_i)，具体地，首先将M_s(X_i，Y_i，Z_i)投影到图像坐标中为q(x，y)，然后计算M_s＝P^-1q，其中P＝P_l[R|t]，P_l＝P_NP_M，

步骤4-1-6：循环加入图像I_l(l＝3，...，k)，直到所有图像加入完毕。

步骤4-2：用光束平差方法优化空间坐标M_s(X_i，Y_i，Z_i)，方法如下：

步骤4-2-1：计算空间坐标M_s(X_i，Y_i，Z_i)的重投影误差K，其由每幅图像所有距离误差的平方和相加，具体地，首先计算p_il，q_il和r_il：

然后计算

其中，I_il为可见标志值，第i点在第l幅图像中有像素坐标时，I_il＝1。第i点在第l幅图像中没有像素坐标时，I_il＝0；(x_il，y_il)是第i点在第l幅图像中的像素坐标，i＝1，...，N。

步骤4-2-2：依据非线性最小二乘法修正空间坐标M_s(X_i，Y_i，Z_i)以降低步骤4-2-1中的重投影误差K，直到误差函数K的值达到最小为止。

步骤4-3：输出误差函数K最小值对应的空间坐标M_s(X_i，Y_i，Z_i)的集合，记为点集G，具体结果如图6所示。

步骤5：利用逐点插入法对三维点云进行Delaunay三角化处理，其流程示意如图7所示。

所述步骤5中，具体包括如下步骤：

步骤5-1：读取建筑物的空间坐标点集G，构造一个包含所有数据点的超三角形C，将该超三角形C放入Delaunay三角形链表中，Delaunay三角形链表中任一三角形满足任意四点不共圆以及最小内角最大的要求，链表最终生成建筑物的网格模型。

步骤5-2：记超三角形C内某一正在处理的数据点为Q，将点Q插入Delaunay三角形链表中，步骤如下：

步骤5-2-1：在Delaunay三角形链表中找出其外接圆包含点Q的三角形，称为点Q的影响三角形。

步骤5-2-2：删除点Q所有影响三角形之间的公共边。

步骤5-2-3：连接点Q与所有影响三角形的顶点，依据任意四点不共圆和最小内角最大的要求进行筛选，将符合Delaunay三角形要求的新三角形放入Delaunay三角形链表中。

步骤5-3：循环加入超三角形C中新的数据点到Delaunay三角形链表中，执行步骤5-2，直到超三角形C中数据点全部加入Delaunay三角形链表为止，生成完整的建筑物点云三角网模型，具体结果如图8所示。

步骤6：最后对构建的Delaunay三角网模型利用积分法计算建筑物的体积。

所述步骤6中，具体包括如下步骤：

步骤6-1：读取步骤5生成的三角网模型。

步骤6-2：将模型顶部的投影分割成n个正方形，通过投影形成n个长方体，高为h_i。

步骤6-3：任取1个小长方体，假设采样间隔为d_l，其数值等于投影面的边长dx，即d_l＝dx＝dy，投影示意如图9所示。

步骤6-4：计算建筑物的体积

以上所述仅为本发明的较佳实施方式，本发明的保护范围并不以上述实施方式为限，但凡本领域普通技术人员根据本发明所揭示内容所作的等效修饰或变化，皆应纳入权利要求书中记载的保护范围内。

Claims (6)

1.一种基于无人机图片的建筑物体积测量方法，其特征在于：包括以下步骤：

步骤1：利用无人机对建筑物进行环绕的图像信息采集，获取建筑物多幅图片序列；

步骤2：经双边滤波处理图像序列后，采用Harris角点检测算法获取图像序列的角点；

步骤3：先采用零均值互相关法进行角点粗匹配，建立一对多的对应关系；再采用松弛迭代法处理粗匹配角点，得到一对一精匹配的角点对；

步骤4：利用从运动中恢复结构算法，循环计算匹配点的空间坐标，并利用光束平差法对空间坐标进行优化，得到建筑物的三维点云；

步骤5：利用逐点插入法对三维点云进行Delaunay三角化处理；

步骤6：最后对构建的Delaunay三角网模型利用积分法计算建筑物的体积。

2.根据权利要求1所述的一种基于无人机图片的建筑物体积测量方法，其特征在于：所述步骤2中，具体包括如下步骤：

步骤2-1：读取图像序列，采用双边滤波做图像去噪声，并对图像做灰度化处理，处理后的图像序列记为I_l(l＝l，2...，k)；

步骤2-2：计算图像I_l上每个像素点水平方向和垂直方向上的一阶导，即梯度I_x和I_y；

步骤2-3：设置矩阵M，

步骤2-4：平滑M得到M₁，方法如下：

式中为高斯模板；

步骤2-5：判断像素点是否为角点，具体地，首先计算响应函数W＝det M₁-k₁(trace M₁)²；其中，k₁取经验数，取值范围为0.04～0.06，其选值影响函数曲线形状的直观性，为便于计算设置为0.04；然后设置阈值T，其大小为所有像素点响应函数最大值的0.01倍，为便于计算设置为0.01，当W≥T时，该像素点则为角点；

步骤2-6：标记角点位置，输出角点。

3.根据权利要求1所述的一种基于无人机图片的建筑物体积测量方法，其特征在于：所述步骤3中，具体包括如下步骤：

步骤3-1：用零均值归一化互相关方法，依据角点的灰度值信息粗匹配角点，方法如下：

步骤3-1-1：获取图像I₁中的角点p(x，y)，其相同位置在图像I₂的搜寻区域尺寸大小(2l_x+1)(2l_y+1)下的角点为q(x′，y′)，l_x和l_y分别是的图像宽度和高度；

步骤3-1-2：获取图像中角点对应的灰度值H_l(x，y)；

步骤3-1-3：计算两个角点的互相关性：

其中，u_l(x，y)和σ_l(x，y)表示相关模板(2m+1)(2n+1)内灰度平均值和标准偏差，m＝n＝3表示角点局部窗口大小；

步骤3-1-4：根据两个角点局部区域内灰度值越相近，互相关值越近似为1，可以作为候选匹配点的要求，将阈值T₁设置为0.9，互相关值S_pq＞T₁的匹配点视为粗匹配点对，记为点集A；

步骤3-2：用松弛迭代法，依据匹配强度精匹配角点，其流程步骤如下：

步骤3-2-1：计算粗匹配点对的匹配强度：

其中，(p，q)和(f，g)表示粗匹配点对，N(p)和N(q)分别是以点p和q为中心、X为半径的邻域，f∈N(p)且g∈N(q)，δ表示相似度贡献度，d表示粗匹配点对平均距离，r表示粗匹配点对的相对距离偏差，r非常大时可以忽略匹配点(f，g)，ε_r表示相对距离偏差的阈值T₂，T₂的取值影响匹配的精度，设置T₂为0.3；

步骤3-2-2：判断粗匹配角点的唯一性，具体地，首先按照匹配强度大小从大到小对匹配点对排序，然后选取匹配强度最大和第二大的点对D₁和D₂，并计算选取匹配强度在D₁与D_p之间的点对为精匹配点对；

步骤3-2-3：删除点集A中正确的角点，从步骤3-2-1开始迭代剩下的点直到点集A为空；

步骤3-3：输出精匹配对的坐标点集合，记为点集B。

4.根据权利要求1所述的一种基于无人机图片的建筑物体积测量方法，其特征在于：所述步骤4中，具体包括如下步骤：

步骤4-1：通过投影关系计算建筑物图像上匹配点的空间坐标，步骤如下：

步骤4-1-1：对相机标定，获取相机的焦距f_d；

步骤4-1-2：计算基础矩阵F，具体地，首先获取点集B中图像I₁和I₂上的匹配点对m(x，y)和m′(x′，y′)，然后记m(x，y)在x，y轴方向上的距离尺寸为z_x和z_y，其对应像素坐标为(u，v)，记带入m和m′的齐次坐标计算m′^TFm＝0，得到基础矩阵F；

步骤4-1-3：计算本质矩阵E＝P_N ^TFP_N，其中

步骤4-1-4：对本质矩阵E奇异值分解得到旋转矩阵R和平移矩阵t；

步骤4-1-5：计算建筑物图像上匹配点的空间坐标M_s(X_i，Y_i，Z_i)，具体地，首先将M_s(X_i，Y_i，Z_i)投影到图像坐标中为q(x，y)，然后计算M_s＝P^-1q，其中P＝P_l[R|t]，P_l＝P_NP_M，

步骤4-1-6：循环加入图像I_l(l＝3，...，k)，直到所有图像加入完毕；

步骤4-2：用光束平差方法优化空间坐标M_s(X_i，Y_i，Z_i)，方法如下：

步骤4-2-1：计算空间坐标M_s(X_i，Y_i，Z_i)的重投影误差K，其由每幅图像所有距离误差的平方和相加，具体地，首先计算p_il，q_il和r_il：

然后计算

其中，I_il为可见标志值，第i点在第l幅图像中有像素坐标时，I_il＝1；第i点在第l幅图像中没有像素坐标时，I_il＝0；(x_il，y_il)是第i点在第l幅图像中的像素坐标，i＝1，...，N；

步骤4-2-2：依据非线性最小二乘法修正空间坐标M_s(X_i，Y_i，Z_i)以降低步骤4-2-1中的重投影误差K，直到误差函数K的值达到最小为止；

步骤4-3：输出误差函数K最小值对应的空间坐标M_s(X_i，X_i，Z_i)的集合，记为点集G。

5.根据权利要求1所述的一种基于无人机图片的建筑物体积测量方法，其特征在于：所述步骤5中，具体包括如下步骤：

步骤5-1：读取建筑物的空间坐标点集G，构造一个包含所有数据点的超三角形C，将该超三角形C放入Delaunay三角形链表中，Delaunay三角形链表中任一三角形满足任意四点不共圆以及最小内角最大的要求，链表最终生成建筑物的网格模型；

步骤5-2：记超三角形C内某一正在处理的数据点为Q，将点Q插入Delaunay三角形链表中，步骤如下：

步骤5-2-1：在Delaunay三角形链表中找出其外接圆包含点Q的三角形，称为点Q的影响三角形；

步骤5-2-2：删除点Q所有影响三角形之间的公共边；

步骤5-2-3：连接点Q与所有影响三角形的顶点，依据任意四点不共圆和最小内角最大的要求进行筛选，将符合Delaunay三角形要求的新三角形放入Delaunay三角形链表中；

步骤5-3：循环加入超三角形C中新的数据点到Delaunay三角形链表中，执行步骤5-2，直到超三角形C中数据点全部加入Delaunay三角形链表为止，生成完整的建筑物点云三角网模型。

6.根据权利要求1所述的一种基于无人机图片的建筑物体积测量方法，其特征在于：所述步骤6中，具体包括如下步骤：

步骤6-1：读取步骤5生成的三角网模型；

步骤6-2：将模型顶部的投影分割成n个正方形，通过投影形成n个长方体，高为h_i；

步骤6-3：任取1个小长方体，假设采样间隔为d_l，其数值等于投影面的边长d_x，即d_l＝dx＝dy；

步骤6-4：计算建筑物的体积