朱鹏程
(1. 陕西铁路工程职业技术学院,陕西 渭南 714000)
煤炭资源作为我国主要工业能源之一,在为经济建设做出重要贡献的同时也给环境带来了严重破坏,主要包括工业污染[1]、开采沉陷[2]、水土流失[3]等。由于井下采煤导致煤层上覆岩层发生断裂与位移,且随着开采范围的增大,岩层移动波及地表,引起地表位移,不可避免地对采动影响区域内的建(构)筑物、耕地、人类生产活动产生影响。为定量分析采动对地表的影响以及二者之间的对应关系,相关学者开展了大量研究,如李昱昊[4]等以无人机摄影测量技术为基础,利用BP 神经网络去除了沉降盆地内的噪点,进而提高了下沉盆地监测精度;廉旭刚[5]等研究了免像控无人机摄影测量技术监测开采沉陷的精度;杨绪霆[6]等通过地形跟随模式提高了无人机摄影测量方法监测矿区沉陷的精度。除无人机监测外,InSAR 技术也被广泛应用于地表移动监测,如谭志祥[7]、沙永莲[8]等采用SBAS-InSAR 技术实现了沉陷盆地的毫米级监测;张童康[9]等将InSAR 与支持向量机相结合,改进了开采沉陷预测模型。三维激光扫描技术也是常见的开采沉陷监测方法,如王磊[10]、贾秉松[11]等分别将其应用于矿区建筑物的变形监测和采动损害评估以及受采动影响的山区边坡沉陷监测。已有研究几乎均是通过某种技术手段监测地表某点的纵向位移,即单点的沉降量,但地表某点受采动影响的位移是一个三维移动的过程[12],既包含纵向位移,又包含横向位移(水平位移),二者叠加才是其真实的位移量。为了弥补已有监测方案的缺陷,更加全面地描述地表受采动影响的变化过程,本文通过无人机摄影测量技术与计算机视觉技术相结合的方法,实现了采动影响区域内地表的水平移动监测,可为井下开采、沉陷预计、人员安置等提供更全面的数据。
本文提出的无人机矿区地表水平移动监测方法的技术流程见图1,首先利用开采沉陷预计方法预计采动影响范围,并根据监测精度确定无人机作业的航高、重叠度、航线等参数;然后采集预计影响范围内不同时期的时序影像,通过Pix4Dmapper处理影像获得时序DOM;再利用Global Mapper 对DOM 进行网格化处理,并进行特征点匹配,计算特征点在不同时序影像上的位置;最后计算特征点的水平移动距离与方向,绘制水平移动矢量图。
图1 无人机矿区地表水平移动监测方法流程图
原始影像经Pix4Dmapper 处理后生成DOM,再利用Global Mapper 进行网格化处理,具体流程见图2,网格大小由沉陷预计的最大水平移动量ΔDmax确定,假设ΔDmax=500 mm,则网格边长为1.2ΔDmax即600 mm。划分网格大小的目的:①减小特征点匹配时的运算量,只需在对应网格内匹配同名特征点即可;②确定时序移动增量监测范围,随着采动的影响,地表发生水平移动与沉降后,其特征会与初期DOM 特征差异巨大,导致无法匹配,划定网格后只需匹配当前网格与上期对应网格,通过计算移动增量并累加增量的方法获得当期与初期之间的水平移动大小与方向。
图2 网格化处理
常见的特征点匹配算法包括Harris 角点检测算法、SIFT 算法、SURF 算法、FAST 角点检测算法、BRIEF 描述子以及基于FAST 角点检测算法和BRIEF描述子相结合的ORB算法等。大量学者已将上述算法应用于各场景中并做出了相关改进,如王丞[13]等利用Harris 角点检测算法实现了三维点云的自适应特征描述、提取与粗配准;张占平[14]等将Harris 算子与SIFT描述子相结合,并辅以无人机影像的POS数据与公开的SRTM 数据,提高了无人机倾斜摄影影像的匹配效率;徐启文[15]等利用一种具有动态阈值的改进SURF算法提升了图像拼接效率;陈伟[16]等结合FAST提取的特征点和SURF 算法的描述子实现了图像的快速拼接;陶卓[17]等利用稀疏光流法改进了传统ORB 算法,使算法的匹配精度更高、鲁棒性更强。本文对同一场景进行了常见特征检测与提取算法的测试,并统计了5种特征点提取算法的提取效果,结果见表1。通过对比提取效率、特征点提取数量以及特征点匹配准确性,最终确定本文的特征点提取与匹配采用经典SIFT算法。
表1 常见特征检测与提取算法比较
利用SIFT算子匹配两幅影像中同名特征点后,需将特征点像素坐标转换为地理坐标,本文采用GDAL空间数据转换坐标库进行坐标转换,转换公式为:
式中,trans(0)、trans(3)为影像左上角像素点左上角位置的地理坐标;trans(1)、trans(5)为影像上单个像素的宽度与高度;如若该影像是指北的,则trans(2)、trans(4)均为0;m、n为特征点在影像中的像素坐标;x、y为该特征点转换后的地理坐标。
将同名特征点的像素坐标转换为地理坐标后,便可计算该特征点的水平移动量和方向,假设同一网格匹配n对特征点,各对特征点水平移动量的计算公式为:
为使水平移动监测更加可视化,可计算水平移动方向,假定某点水平移动是由(x1'y1)指向(x2y2),且水平移动矢量与东方向形成的夹角为α(图3),则有:
图3 水平移动矢量示意图
由于矿区地物特征的特殊性和复杂性,采用SIFT算子提取的匹配特征点存在误差(图4a),为提高水平移动监测精度,本文提出了一种基于位移大小和方向的误匹配剔除方法。通常在同一网格内,各特征点间水平移动量差别不大,且水平移动方向具有一致性,根据这两个特性可剔除误匹配和不符合地表水平移动规律的匹配。首先计算同一网格内所有水平移动量的平均值ΔA,设置0.9ΔA~1.1ΔA为满足误匹配剔除的大小阈值范围;再计算同一网格内所有水平移动方向的平均值Δα,设置0.9ΔA~1.1ΔA为满足误匹配剔除的角度阈值范围;最后结合两个阈值范围,可基本剔除误匹配,剔除后效果见图4b。
图4 误匹配剔除前后对比
监测案例中采动区域为山东省济宁市南屯煤矿33上08工作面(图5),测区中心位于35°22′43″N、116°52′01″E,平均采深为300 m、走向长680 m、倾向长170 m、平均煤厚5.2 m。该工作面采动前已对地表村庄易地搬迁,但仍保留部分工业建筑,且该区域为高潜水位地区,采动后地下水会出露地表,因此对该区域的地表移动进行监测非常必要。
图5 监测区域
基于概率积分法的开采沉陷预计方法是目前较成熟的矿山开采沉陷预计手段,主要基于已有变形监测数据与对应工作面参数,根据随机介质理论,以统计学的视角将采动区域分解为无限个微小开采单元,进而反演计算工作面参数与变形监测数据之间的相互关系,主要通过下沉系数、水平移动系数、拐点偏移距、主要影响角正切以及开采影响传播角等预计参数进行描述和计算,并以此预计同地质采矿条件下的其他待采动工作面的变形情况[18]。由3308工作面参数与历史预计参数计算可知,该工作面预计最大下沉值为4 549 mm,东西最大水平移动为1 396 mm,南北最大水平移动为1 455 mm(图6)。
图6 33上08工作面开采地表下沉等值线
为监测3308工作面采动对地表的影响,检验本文方法的有效性,2021 年3—4 月对3308 工作面采动影响区域采集3 期无人机影像数据与像控点坐标数据。无人机采用大疆精灵RTK,作业高度为100 m,影像分辨率为2.73 cm/pixel,具体飞行数据见表2。像控点采用60 cm×60 cm 倒三角形像控点(图7),并采用GPS-RTK+DS3 级水准测量获得坐标与高程,其中平面误差为2 cm,高程误差为1 mm,共布设35 个像控点(图8),包括10 个像控点和25 个检查点。
表2 无人机飞行数据
图7 33上08工作面像控点布设与量测
图8 像控点布设位置分布图
利用本文方法对3 期数据进行处理,获得水平移动大小与方向图(图9),并将实测值与提取值进行比较,统计本文方法的水平移动提取误差(表3),首期水平移动由二期DOM 与初期DOM 比较获得,二期水平移动由三期DOM 与二期DOM 比较获得,经计算首期水平移动提取中误差为1.9 cm,二期水平移动提取中误差为3.3 cm,综合误差为3.8 cm。
表3 水平移动提取误差统计
图9 两期地表移动监测结果
本文基于无人机摄影测量技术提出了一种矿区地表水平移动监测方法,并通过实测案例检验了该方法的可行性。
1)地表水平移动作为受采动影响的地表变形的重要组成部分,前人的相关研究重点关注地表沉降监测,而对水平移动监测的研究存在缺失,本文是对相关研究的补充。
2)本文方法可基本完成地表水平移动的自动化监测,较少需要人为干预,且精度可满足矿区日常监测要求。在影像分辨率为2.73 cm/pixel 的前提下,水平移动提取精度为3.8 cm,且随着硬件设备性能的提升,该方法可获得更高精度的结果。
本文为受采动影响的地表水平移动监测提供了一种更全面、便利、可视化的监测手段,可广泛应用于采动区域的变形监测,为井下开采、耕地保护和搬迁移民等提供前瞻性和实时的数据支持,对于井下安全开采、采动区治理以及采动影响范围内的生命财产安全等具有重大意义。