无人机DEM测点插值算法的矿区沉陷测绘技术

朱 涛

(河南省地质局矿产资源勘查中心，郑州 454000)

0 引言

DEM就是指数字高程模型，可以通过数字化模拟地形高程数据的方式，确定地面地形曲线的走向特征，简单来说，就是利用有序数值阵列来表示实体地面模型。分辨率是DEM模型刻画地形曲线的重要参考指标，其取值水平既影响数字化模拟手段的实施准确性，也影响有序数值阵列的实际排序形式[1]。无人机DEM测点插值算法是以DEM模型为基础提出的测量数据处理思想，对于所得地面地形数据可以进行按需排序，并可以在不改变地形曲线走向特征的前提下，得到不同的参考坐标取值结果，从而在一次计算的过程中，得到更多的数据样本坐标值。近年来，无人机DEM测点插值算法迅速发展，在工程建设、地质勘探等多个领域中都得到了广泛应用[2]。此外，为保证数据坐标取值能够准确表现出地面地形曲线的走向特征，还要求DEM模型不能在同一区域内对数值阵列向量进行重复取值。

矿区沉陷就是指由煤矿地下采矿行为引起的地表塌陷现象，矿产资源被开采出来之后，矿区周围岩体的力学平衡受到了破坏，而随着这种破坏作用的持续累积，矿区地表就会表现出区域性塌陷的情况。在实际应用过程中，如何根据矿区地表的实际沉陷程度，而制定具有针对性的测绘技术方案成为了一项亟待解决的难题。基于无人机激光雷达技术的沉陷监测方法根据矿区沉陷特征的三维点云排列形式，建立多期数字高程模型，再通过DEM指标连续相减的方式，确定核心沉陷区地表的下陷程度[3]。基于HOG特征的检测方法提取施工区域的边缘特征，并借助训练支持向量机原则，对相关特征参量进行分类，以便于实现对矿区沉陷程度的精准测量[4]。然而在上述两种方法作用下，所得沉陷深度测量结果与真实沉陷深度之间的差值水平并不满足实际应用需求，故而按需所制定测绘技术的适用性也就相对较为局限。为解决上述问题，针对无人机DEM测点插值算法的矿区沉陷测绘技术展开研究。

1 无人机DEM测点插值算法

无人机DEM测点插值算法是完善矿区沉陷测绘技术实施方案的基础，本章节将分别从DEM测点获取、分维值测定、插值拐点求解3个方面，分析该算法的应用流程。

1.1 DEM测点获取

实施无人机测绘时，所得DEM测点数据包括平面位置、高程两类信息参量[5]。其中，平面位置类的DEM测点数据对应矿区地表及外部岩体的基本地质特征，所得信息参量仅具有横轴、纵轴两个方向上的取值结果；高程类DEM测点数据对应矿区无人机航摄图像中的所有特征，所得信息参量具有横轴、纵轴、空间轴3个方向上的取值结果[6]。

对于平面位置信息的定义式公为：

(1)

oX、oY分别表示横轴、纵轴方向上的矿区地质特征，δX表示横轴方向上的测点参数，δY表示纵轴方向上的测点参数。

对于高程信息的定义式为：

(2)

oZ表示空间轴方向上的矿区地质特征，δZ表示空间轴方向上的测点参数。

联立公式(1)、公式(2)，可将无人机DEM测点获取表达式定义为：

(3)

ΔO表示无人机航摄图像中基本地质特征的单位累积量，χ表示数据样本测量参数。由于矿区的沉陷程度并不可控，所以oX、oY、oZ及δX、δY、δZ的取值存在相等的可能。

1.2 分维值测定

对于无人机DEM测点插值算法而言，分维值决定了测点参数之间的维数关系，且出于精确性考虑，分维值参数的取值应属于(0，1]的数值区间[7]。分维值等于“0”，表示所获取DEM测点信息完全不能描述矿区的沉陷特征，在实际应用过程中，即便是清晰度水平极低的无人机航摄图像中也不可能完全不存在可用信息，所以这种情况并不会出现。分维值等于“1”，表示所获取DEM测点信息能够100%精准描述出矿区的沉陷特征，对于清晰度水平极高的无人机航摄图像而言，这种表现情况是可能出现的，只是出现概率相对较小[8]。对于分维值测定表达式的求解满足公式(4)。

(4)

1.3 插值拐点

插值拐点是指矿区无人机航摄图像中DEM测点数据取值迅速发生改变的节点。所谓插值可以理解为安插在某一节点处的DEM测点数据，该数据样本虽然具有当前节点处矿区沉陷地质的全部特征，但由于其取值空间为无人机航摄图像，所以测点数据结果仅具有运算作用，并不能完全代替矿区沉陷特征[9]。图1反映了插值拐点的取样原则。

图1 插值拐点取样原则

在无人机航摄图像中，测绘方案实施方向并不一定与摄影方向保持一致，所以插值拐点之后DEM测点数据的取值既有增大也有缩小的可能(图1中1曲线表示DEM测点数据取值的增大变化状态，2曲线表示DEM测点数据取值的缩小变化状态)[10]。

设U0表示DEM测点数据的初始取值，Umax表示增大曲线上的DEM测点数据取值，Umin表示缩小曲线上的DEM测点数据取值，联立公式(4)，可将插值拐点定义式表示为：

(5)

定义插值拐点时，如果不参考分维值测定条件，则有可能导致测绘数据与矿区地表真实沉陷情况出现偏差。

2 矿区沉陷区域地理模型

在无人机DEM测点插值算法的作用下，构建矿区沉陷区域地理模型，还需根据像控点布设情况，定义空中三角，再联合相关数据参量，求解径向插值基函数表达式。

2.1 像控点布设

像控点布设就是在矿区沉陷航摄图像中，按照无人机DEM测点插值算法所定义的测绘节点，一般来说，为避免插值拐点取样结果无法真实反映出矿区地表的沉陷特征，在每一航摄方向上，都至少取两个像控点作为实际布设位置[11]。如图2为标准的像控点取样结果。

图2 像控点取样

对图2中的矿区沉陷航摄图像进行背景去除处理，得到如图3所示的像控点布设图像。

根据图3可知，对像控点进行布设时，应在矿区沉陷航摄图像中心定义一个标准控制区域，且该区域距离上、下两端边界的距离完全相等，距离左、右两端边界的距离也完全相等。航摄图像上、下端存在两组完全对称的像控点，且这些节点的布设位置非常靠近图像边界[12]。左、右两端像控点布设位置与标准控制区域的圆心等高，且这两个节点的存在状态也是完全对称的。

2.2 空中三角

空中三角是辅助矿区沉陷测绘技术实施的必要条件，可以对已布设像控点对象进行分别取样，再联合所有数值结果，定义一个完整的空间数据集合，由于任意3个像控点都可以连接成一个三角形，所以该空间数据集合被简称为“空中三角”[13]。简单来说，一个空中三角集合中所包含的全部数据样本都描述矿区地表的沉陷特征，且这些数据的获取完全依照无人机DEM测点插值算法。规定wX、wY、wZ分别表示横轴、纵轴、空间中方向上的地表沉陷数据取值，其定义式如下：

(6)

其中，γX、γY、γZ分别表示横轴、纵轴、空间中方向上的插值向量，RX、RY、RZ分别表示横轴、纵轴、空间中方向上的测绘数据样本。

在公式(6)的基础上，设φ表示矿区航测图像沉陷特征的全域测定参数，φ表示空间性标记参量，联立上述物理量，可将空中三角表达式定义为：

(7)

为保证数据样本空间的广域性，一般不会对同一类型像控点中的沉陷特征进行重复取样[14]。

2.3 径向插值基函数

(8)

如果测绘区域的划定范围较大，利用无人机DEM测点插值算法对沉陷数据进行取值时，极有可能出现DEM测点数据与插值点高程信息不匹配的问题，此时空中三角集合中会出现不满足径向插值基函数表达式的数据样本[16]。为避免上述情况的发生，在定义径向插值基函数表达式之前，应利用空中三角集合中的取样信息进行重复计算，直至将非合理信息参量完全去除。

3 测绘实施方法设计

测绘实施方法设计就是按照空间坐标系转换原则，对矢量数据进行叠加处理，再联合像片倾角与旋偏角，确定测量精度评价指标的取值范围，本章节将针对上述内容展开研究。

3.1 空间坐标系转换

空间坐标系转换是按照无人机DEM测点插值算法定义原则，确定在测绘矿区发生沉陷时，横向、纵向、空间向坐标轴之间的角度关系。转换之前，横向、纵向、空间向坐标轴之间的物理夹角均等于90°，这种定义模式只能适应微倾斜情况下的测绘行为，当无人机航摄方向处于水平或竖直状态时，至少有一个坐标轴在航摄图像中的投影为一个点，而二维图像中的点并不能描述出DEM测点数据的高程信息[17]。在转换后的空间坐标系中，横向、纵向、空间向坐标轴之间的夹角呈现出增大或缩小的数值变化状态，在面对水平或竖直方向的测绘射线时，坐标轴在航摄图像中的投影至少是具有一定程度的线段，而线段具有描述DEM测点数据高程信息的能力[18]。

对于空间坐标系转换原则的定义满足如下表达式：

(9)

3.2 矢量数据叠加

在实施测绘的过程中，矢量数据叠加就是将相似的DEM测点数据整合到一起，且处理后数据样本的高程信息并不会发生改变[19]。对于矢量数据的叠加处理主要遵循如图4所示的原则。

图4 矢量数据叠加原则

无人机DEM测点插值算法规定，在空中三角集合中提取出的沉陷深度测量值保持矢量状态。所谓矢量值是指同时具有方向特性与数值特性的数据样本，在空间坐标系中，与无人机航测方向保持一致的数据样本的取值方向为“+”，与无人机航测方向保持相反的数据样本的取值方向为“-”[20]。两个数值相同但取值方向相反的测绘数据样本必须处于同一沉陷区域内，但两个取值方向相同、数值不同的测绘数据样本有可能不属于同一沉陷区域。设l1、l2、…、ln表示n个不相等的测绘数据样本，κ1、κ2、…、κn分别表示与测绘数据样本匹配的叠加系数，联立公式(9)，可将矢量数据叠加原则表示为：

(10)

若存在数据样本与DEM测点信息不完全匹配的情况，则表示叠加处理后，会剩余取值方向与数值水平都不固定的数据参量。

3.3 像片倾角与旋偏角

像片倾角是指无人机航摄图像中测绘射线与沉陷深度所在地面水平线之间的物理夹角，如图5中的ω。旋偏角是指测绘射线全旋角与像片倾角之间的差值，在图5中可以表示为μ-ω。实施矿区沉陷测绘时，像片倾角的数值水平越大，就表示空间坐标系转换处理过程中，坐标轴之间的物理夹角越大；而旋偏角越大，则表示DEM测点数据高程信息与真实沉陷数据之间的差值水平越大[21]。为实现对矿区沉陷深度的精准测量，应同时控制像片倾角、旋偏角的实际取值。具体的像片倾角与旋偏角定义原则如图5所示。

图5 像片倾角与旋偏角定义

在无人机航摄图像中，全旋角数值一般等于180°，而像片倾角数值则小于90°，因此旋偏角取值属于90°-180°的数值区间[22]。

3.4 测量精度评价指标

测量精度评价指标影响无人机DEM测点插值算法对矿区沉陷程度的测绘处理能力，由于DEM测点数据高程信息始终保持定值状态，所以即便在测绘方案不符合实际应用需求的情况下，测量精度评价指标的计算数值也不会表现出波动变化行为[23-24]。如果沉陷区域面积保持不变，那么像片倾角与旋偏角数值也可以保持定值状态，当前情况下求解所得的测量精度评价指标也就能够符合测绘矿区的实际应用需求。对于测量精度评价指标的计算满足如下表达式：

(11)

4 实例分析

本次实验的主要目的是分析无人机DEM测点插值算法的矿区沉陷测绘技术、基于无人机激光雷达技术的沉陷监测方法、基于HOG特征的检测方法的实用性能力，在不考虑其他干扰条件的情况下，设计如下对比实验。

4.1 矿区情况

选择位于中海拔地区的某土石矿区作为实验环境，利用遥感技术对该区域的地质特征进行检测，详情如图6所示。

图6 矿区地质遥感图像

图6中部、下部两块小型区域的颜色明显比其他遥感区域的颜色更深，根据遥感成像特点可知，这两个区域为主要沉陷区，在上述两区域外侧划定一片区域作为本次测绘的具体实验区。为避免开采行为对矿区地质造成严重影响，实验过程中必须借助遥感技术对矿区沉陷程度进行实时监测。

4.2 沉陷测绘点取样

在实验区域内选择8个不重复、深度水平也不相同的测绘节点作为实验对象，并对每一节点处的地表沉陷深度进行测量，详情见表1。

表1 测绘点沉陷深度 m

分析表1可知，该实验区域内，地表沉陷程度的差异性较大，4号测绘点的最大沉陷深度与7号测绘点的最小沉陷深度之间的数值差达到了22.51 m。

4.3 实验结果分析

分别利用无人机DEM测点插值算法的矿区沉陷测绘技术、基于无人机激光雷达技术的沉陷监测方法、基于HOG特征的检测方法对8个测绘节点处的地表沉陷深度进行测量，其中第一组为实验组方法、第二组为对照A组方法、第三组为对照B组方法。

图7记录了实验组、对照组矿区地表沉陷深度的具体实验数值。

图7 矿区地表沉陷深度实验数值

分析图7可知，第7个测绘节点处，对照A组、对照B组沉陷深度测量值均为零，在精度分析方面的参考价值不大，因此在后续分析过程中，必须去除该测量结果。实验组方法所测得的矿区地表沉陷深度均值水平相对较高，对照A组方法次之，而对照B组测量值最小。

取第4测绘节点处的最大实验数值进行计算，可知沉陷深度测量结果与矿区真实沉陷深度之间的差值为：

(12)

其中，M1表示实验组差值，M2表示对照A组差值，M3表示对照B组差值，单位均为m。

根据公式(12)可知，实验组方法作用下，沉陷深度测量结果与矿区真实沉陷深度之间的差值最小，对照A组方法次之，对照B组方法最大。

由上述实验结果，可得出以下结论：

1)基于无人机激光雷达技术的沉陷监测方法的应用，对于沉陷深度测量结果与矿区真实沉陷深度之间差值水平的控制能力有限，不足以实现对矿区沉陷深度的精准测量；

2)相较于基于无人机激光雷达技术的沉陷监测方法，基于HOG特征的检测方法的测量能力相对较强，能够适当控制沉陷深度测量结果与矿区真实沉陷深度之间的差值水平，但依然达不到完善测绘技术实施方案的实际应用需求；

3)无人机DEM测点插值算法的矿区沉陷测绘技术的应用，可有效控制沉陷深度测量结果与矿区真实沉陷深度之间的差值水平，使其计算数值不超过1 m，不但实现了对矿区沉陷深度的精准测量，还可以对测绘技术具体实施方案进行完善，与无人机激光雷达技术的沉陷监测方法、基于无人机激光雷达技术的沉陷监测方法相比，更符合实际应用需求。

5 结束语

新型矿区沉陷测绘技术在无人机DEM测点插值算法的基础上，通过计算分维值测定结果的方式，确定插值拐点取样条件，又联合径向插值基函数表达式，定义空间坐标系转换原则，从而推导出测量精度评价指标的实际取值范围。实用性方面，无人机DEM测点插值算法的矿区沉陷测绘技术的作用下，沉陷深度测量结果与矿区真实沉陷深度之间的差值水平得到了有效控制，有利于实现对矿区沉陷深度的精准测量，对于完善测绘技术的具体实施方案，可以起到一定的促进性影响作用。