无人机LiDAR在采空区沉陷监测中的应用

郑俊良,姚顽强,蔺小虎,张联队,张 咏,张 冬

(1.西安科技大学 测绘科学与技术学院,陕西 西安 710054;2.陕西彬长矿业集团有限公司,陕西 咸阳 712000;3.陕西延长石油靖边煤业有限公司,陕西 榆林 718500)

0 引 言

随着能源生产和消费革命的推进,煤炭在能源结构中的占比逐渐减少,但在一段时间内煤炭仍属于中国主体能源,具有重要战略地位[1-2]。西部煤矿地处黄土高原中上游,地形复杂多样,属于典型的西部黄土沟壑区,生态环境十分脆弱。机械化综采、垮落式顶板管理的采煤方法造成上覆岩层破裂、地表塌陷,引起地下水流失、土地沙漠化,破坏沉陷区生态环境、威胁地表建(构)筑物安全[3-5]。因此,针对西部煤矿所在的黄土沟壑区,开展可靠的矿区开采沉陷监测,准确掌握开采沉陷动态变化过程,对开采沉陷预计、地表建(构)筑物损坏评价、生态修复具有重要意义。

目前,开采沉陷监测方法包括传统预埋监测桩观测(水准测量、全站仪测量、GNSS测量等)、InSAR监测、地面三维激光测量、倾斜摄影测量和机载LiDAR测量。传统预埋监测桩观测方法局限于点位观测,具有较高的点位精度(毫米级、厘米级),但不能准确描述采空区地面整体的形变过程,且西部沟壑区地形结构复杂,存在工作量大和人员无法到达的问题[6];InSAR测量能够进行面状观测,且取得较高的测量精度(毫米级),但受限于观测周期,对于地表塌陷速度快、沉陷量大的区域,容易造成失相干现象[7-8];地面三维激光测量能够取得面状厘米级,甚至毫米级观测精度,但由于其视野范围、点云拼接误差等原因,在超长超宽工作面监测时存在工作量大、累积误差大等问题[9];近年来,随着影像密集特征点匹配算法的成熟,无人机倾斜摄影测量技术提高了摄影测量的观测精度,但在地表植被覆盖度高、地形起伏大的区域,易造成监测结果精度低,且精度受限于像控点密度[10-11];机载LiDAR测量具有精度高、观测范围大的优势,但有人机载LiDAR测量存在设备昂贵、天气要求高、空域申请困难等问题;随着无人机技术发展和激光扫描设备改进,无人机LiDAR测量技术逐渐成熟,在满足厘米级测量精度的同时能够极大地降低成本、提高效率,弥补有人机载LiDAR测量的不足[12-13]。无人机LiDAR通过对目标区域低空多时相扫描,可获取高密度点云数据,经点云滤波、插值和求差计算,可获取目标区域高程变化。但在激光点云滤波、插值和求差计算过程中不可避免地存在分类错误、噪声和误差,一些学者通过改进滤波算法和插值方法提高DEM精度,进而改善沉陷模型精度[14-17],但仍存在一定问题：①现有成熟点云滤波算法主要针对传统有人机载LiDAR点云数据,相比较而言,无人机LiDAR点云具有更高的点云密度和更复杂的噪声数据,如在山区应用效果较好的渐进三角网点云滤波算法(Progressire Trianglated Irregular Network,PTIN),在无人机LiDAR点云滤波时,由于高密度和多噪声的原因,易出现种子点错误问题,导致滤波错误;②无人机LiDAR获取的点云、数字高程模型(Digital Elevation Model,DEM)和沉陷模型之间的精度缺乏实测数据验证;在无人机LiDAR测量系统和环境条件相同情况下,沉陷模型精度是否优于原始点云和DEM精度缺乏试验对比;③传统LiDAR监测成果单一,对沉陷模型进一步研究较少,未进行深入的数据分析和挖掘。

1 研究区及数据获取

1.1 研究区概况

凉水井矿位于陕西省榆林市东北部、榆神府矿区东部,处于北纬38°47′～38°53′和东经110°14′～110°21′之间(图1)。井田地面标高1 100～1 326 m,地势西北高东南低。以凉水井煤矿431301工作面为研究对象,开采煤层为4-3煤层,该煤层平均埋深138 m,平均厚度1.3 m,采用长臂式综采、顶板自然塌陷方法开采;其中上覆4-2煤层已于2013年开采完成,地表于2014年达到稳定状态。

图1 研究对象位置Fig.1 Location of study object

基于《煤矿测量规程》(2013)设置地面监测点,结合临近矿区开采沉陷角量参数和现场实际情况,共布设3条监测线：走向观测线、倾向观测线和道路观测线,如图2所示。地面监测点通过埋设预制观测桩,采用四等水准测量,获取监测点多时相高程,实现地表沉陷监测。在工作面停采线附近,选择面积为0.4 km2的区域作为研究区,进行无人机LiDAR试验数据采集和分析。431301工作面回采期间共进行23次全面观测,根据开采沉陷一般规律,选取位于研究区域内的39个地面监测点作为无人机LiDAR参考数据,对无人机LiDAR监测成果进行精度分析,如图2中红色三角形所示。

图2 研究区与监测线布设Fig.2 Layout of study area and monitoring line

1.2 数据获取

无人机LiDAR型号为SZT-R250,集成了RIEGL miniVUX-1激光测距系统和高精度GNSS/IMU(Global Navigation Satellite System/Inertial Navigation System)系统,以及地面GNSS差分系统,扫描测程250 m。无人机有效飞行时间约20 min,试验飞行高度50 m,速度8 m/s,激光脉冲发射频率为100 kHz,扫描速度为100线/s。为了保证多时相数据精度一致,在每次数据采集时均采用相同的飞行和航测参数。在2020年11月和2021年5月共进行2次无人机LiDAR测量,数据采集范围覆盖该时间段地表下沉范围,航测时间与工作面回采位置、试验区域如图3所示。

图3 研究区与采空区位置关系Fig.3 Location relationship between study area and goaf

图3展示了井下巷道、工作面位置、采空区与地表航测范围位置关系。2020年11月23日和2021年5月12日分别进行无人机LiDAR航测,其对应的工作面回采位置如图3蓝色竖线所示。同时航测范围覆盖研究期间采煤导致的地表下沉影响范围,其中采空区1表示在首次数据采集时,已存在的采空区;采空区2表示第2次LiDAR数据采集时产生的采空区;图中黑色三角形表示地面实测监测点。

1.3 研究技术路线

无人机LiDAR开采沉陷监测主要步骤包括外业航测、LiDAR点云滤波、点云插值计算、DEM差值计算、沉陷模型获取及数据挖掘分析等。同时,以地面监测桩测量数据为依据进行LiDAR点云、DEM和沉陷模型精度验证,技术路线如图4所示。

图4 研究技术路线Fig.4 Research technology route

采用渐进三角网点云滤波算法进行点云滤波,该算法一般由局部最低点作为稀疏种子点,创建初始不规则三角网,通过激光点与三角网的距离和角度进行迭代分类。该算法适用性强、滤波效果较佳,但种子点的选择是渐进三角网滤波的关键。因此,结合矿区地表沉陷特征,对种子点的选取方法进行改进,确保种子点的准确性。点云插值计算参考文献[18]所述,选择最适用于黄土沟壑区的Kriging插值算法进行计算,获取地表DEM。点云、DEM和沉陷模型精度采用平均误差(Mean Error,ME)、平均绝对误差(Mean Absolute Error,MAE)、中误差(Root Mean Square Error,RMSE)验证,以水准测量结果作为参考值,计算方法如下

(1)

(2)

(3)

式中Rm为基于DEM提取的下沉值;Zm为地面监测桩观测下沉值;m为参考点数量。

2 LiDAR点云滤波与插值计算

2.1 常用滤波方法

无人机LiDAR能够获取海量地表点云数据,其中包括地面、植被、建筑物等,根据点云的几何特征、色彩、高程变化等信息提取地面点的过程即为点云滤波。点云滤波方法分为：基于坡度滤波(Slope Based Filtering,SBF)[19]、基于渐进三角网滤波[20]、形态学滤波(Morphological Filter,MF)[21]、布料模拟滤波(Cloth Simulation Filter,CSF)[22]等。SBF滤波算法采用点之间的坡度作为阈值参数进行滤波,该方法对于缓坡地形效果较好,但对于地形复杂区域,由于单一阈值的限定,滤波效果不足;MF滤波算法借鉴形态学图像处理方法进行点云滤波,剔除非地面点时不能有效保留地形信息,易造成地形平滑;CSF滤波算法采用点云翻转、模拟布料贴合倒置点云的方法,使用高程阈值分离地面点和非地面点,平地、缓坡和陡坡需采用不同参数,参数确定存在困难;PTIN滤波方法由Axelsson提出,根据种子点建立初始地形模型,对待定点进行距离和角度阈值判断,确定地面点或非地面点,逐步迭代完成点云滤波,该方法对于复杂的山区、森林地形具有较好的适应性,但对种子点的依赖度较高。

PTIN滤波方法具有较好的适用性,特别是对于山区、沟壑区效果较好,但其种子点的选取是滤波的关键[23]。特别是对于煤层开采导致地表出现阶梯裂缝区域种子点的选取,如种子点选取在裂缝或塌陷坑底部会造成局部滤波错误。因此,针对开采沉陷区域梯度变形和伪最低噪声点问题,对种子点选取方法进行改进。

2.2 PTIN种子点选取及滤波

针对矿区,特别是浅埋煤层、大采高开采造成的阶梯塌陷地表,无人机LiDAR渐进三角网点云滤波过程中种子点的选取至关重要。因此,针对矿区阶梯状地表,在渐进三角网滤波算法的基础上,进行种子点选取方法的改进,去除错误的种子点,提高种子点质量,计算过程如图5所示。

图5 改进的PTIN滤波流程Fig.5 Improved PTIN filter flow chart

格网划分主要是用规则的虚拟格网覆盖全部点云数据,将点云按照平面坐标分配到相应的虚拟网格内。种子点的选取是滤波的关键,采用改进的扩展局部最小值方法进行种子点选取,原理如图6所示。

图6 种子点选取计算示意Fig.6 Schematic diagram of seed point selection calculation

图7 局部点云滤波效果Fig.7 Partial point cloud filtering results

对2020年11月和2021年5月点云数据分别进行点云滤波,并对滤波结果进行统计,地面点云占比分别为93%和83%,由于不同时期地表植被覆盖度不同,导致了地面点云占比不同。

2.3 点云插值计算

插值算法优劣直接影响地面DEM精度,目前常用的插值算法包括反距离权重(Inverse Distance Weighted,IDW)、克里金(Kriging)、自然邻域(Natural Neighbor,NN)、样条函数(Spline)、不规则三角网(Triangulated Irregular Network,TIN)等[18]。IDW插值以采样点与插值点间的距离为权重计算加权平均值,离插值点越近的采样点权重越大,超过一定距离时权重忽略不计,适用于采样点均匀且密集的区域[24];Kriging插值与IDW类似,该方法采用半变异函数作为权重进行无偏最优估计,不仅考虑距离权重,还兼顾空间分布关系,适用于空间相关性较好的采样点,但容易出现边缘效应[25];自然邻域插值是基于Voronoi结构的插值计算,该方法对局部特征的继承性较好,但局部点缺失时容易造成失真[26];Spline以最小曲率面逼近各采样点,以获取最佳平滑曲面,顾及大范围的采样点,适用于复杂高曲率曲面,但平滑易造成地形失真[26];不规则三角网(TIN)利用采样点构建不规则三角形,能够在地形复杂区域保留地形特征细节,但其边缘梯度可能存在突然变化[25]。

文献[18]的研究成果表明,在黄土高原沟壑区,Kriging插值在无人机LiDAR点云插值计算中具有较好效果。因此,采用Kriging插值算法对地面点云进行插值计算,并将插值格网大小对结果精度的影响进行试验。采用DEM误差分布和误差值大小对比的方法进行结果分析,不同格网大小引起的DEM误差分布如图8所示。

图8 DEM误差分布Fig.8 DEM error distribution

同时,计算不同格网尺寸DEM的ME,MAE,RMSE,结果见表1。

表1 DEM误差Table 1 DEM errors

由表1可知,ME,MAE和RMSE在0.05～1 m区间内变化不大,当格网尺寸大于1 m时,误差急剧增大。同时,结合图8综合考虑,最终选择0.1 m的格网尺寸进行插值计算,获取较高精度的研究区DEM。

2.4 DEM差值计算

沉陷模型是由多时相地面DEM差值计算获得,与DEM结构和功能类似,能够反映一段时间内地表高程变化,即地表沉陷量,是一种高精度地表沉陷表示模型[13]。通过不同时相DEM之间的差值计算,可获取该段时间内地表的下沉变化。以首次观测时获取的DEM为基础,与开采结束后获取的DEM做减法操作,可得到该段时间内地表的沉陷模型,其三维效果如图9所示,由于地表人为破坏、建筑施工等,导致沉陷量部分异常值过大或过小。

图9 沉陷模型三维效果Fig.9 3D rendering of subsidence model

3 精度评定

基于地面监测桩实测数据对地面点云、DEM和沉陷模型进行精度评定。由于点云由一系列离散点构成,无法进行点对点差值直接计算,因此采用局部平均值的方法进行差值计算,计算监测桩附近局部点云的高程平均值,然后与实测数据对比,结果见表2,2020年11月LiDAR测量点云中误差为60.6 mm,2021年5月为59.9 mm。

表2 点云误差统计Table 2 Point cloud error statistics

其次,基于监测桩实测数据对DEM和沉陷模型误差进行计算,误差分布如图10所示。

图10 点云、DEM,沉陷模型误差分布Fig.10 Errors distribution of cloud,DEM,subsidence model

从图10可以看出,原始地面点云、DEM和沉陷模型三者中DEM的精度较差,沉陷模型精度最优。在无人机LiDAR航测过程中,多时相数据采集时使用相同航测参数,能够消除一定的系统误差,如激光扫描设备系统误差、飞行高度、地形变化等因素的影响,提高沉陷模型的精度。

基于沉陷模型提取地面监测桩的下沉量,与地表监测点监测值进行对比,绘制下沉折线图,如图11所示。

图11 下沉折线Fig.11 Subsidence lines

在2020年11月至2021年5月期间,地面水准测量获取的最大下沉值为1 872 mm,最小下沉值为0 mm,如图11中折线所示;基于沉陷模型获取的最大下沉值为1 860 mm;除A203监测点实测数据下沉异常,实测下沉折线图与基于沉陷模型获取的下沉折线图走势和拐点基本一致。相同位置实测下沉量与沉陷模型下沉量的误差绝对值最大为82 mm,最小为3 mm,平均为47 mm,中位数为42.5 mm,其中小于50 mm的占比60%,误差绝对值整体小于100 mm,误差分布相对集中。

4 基于沉陷模型的数据挖掘与分析

沉陷模型包含丰富的地表下沉信息,包括地表平面位置、任意位置下沉量、下沉变化等,基于沉陷模型进行数据挖掘分析,能够获取开采沉陷引起的下沉范围、边界角、下沉面积、不同程度损害区域面积等。

4.1 下沉量统计

为了准确获取地表不同位置的下沉量和不同下沉量的分布情况,对沉陷模型进行下沉量分级分类统计,结果如图12所示。

图12表示不同下沉量区域分布情况,下沉区域呈椭圆形分布,最大下沉区位于采空区中心,由中心逐渐向外扩散,下沉值逐渐减小,相应的下沉面积逐渐增大。图中圆形标注区域出现异常下沉,经实地调查发现是由于煤层开采期间,人为施工造成。

对下沉区域面积进行分级计算,以大于100 mm的下沉值作为开采沉陷引起的地表下沉,计算结果见表3。下沉值大于100 mm的下沉总面积为180 075 m2,采空区面积为49 789 m2。为准确确定开采沉陷引起的不同下沉量所占面积,暂定以下沉值大于100 mm的沉陷量作为煤层开采引起的地表下沉。结果表明,下沉值小于300 mm的面积占比超过一半以上,随着下沉量的增大,对应的下沉面积逐渐减小。

表3 下沉面积统计Table 3 Statistics of subsidence areas

4.2 下沉剖面信息提取

传统开采沉陷监测采用预埋监测桩的方式,监测煤层开采过程中地表变化,该方法获取的监测数据少、范围小。试验基于沉陷模型采用虚拟监测点和剖面线的形式,分别对沉陷区域进行走向和倾向方向下沉信息提取,虚拟监测点与剖面线位置如图13所示。

图13 虚拟监测点与剖面线位置Fig.13 Location map of virtual monitoring points and profile lines

虚拟监测点以15 m间隔在走向和倾向方向均匀分布,如图13中黑色三角形所示;剖面线分为走向剖面线和倾向剖面线,位置如图13中灰色实线所示。基于沉陷模型在虚拟监测点位提取地表下沉值,并进行折线图绘制,结果如图14所示。虚拟监测点最大下沉值为1 761 mm,走向下沉折线呈半盆地分布,倾向下沉折线图呈全盆地分布,符合开采沉陷下沉规律。

图14 虚拟监测点下沉折线Fig.14 Virtual monitoring point subsidence line

根据沉陷模型分辨率大小,下沉剖面线能够提取大量监测点数据,相比离散监测点具有更好的鲁棒性。因此,在沉陷区域设置一条长740 m的走向剖面线和长413 m的倾向剖面线,基于剖面线以沉陷模型分辨率0.1 m为间隔,提取下沉值。对提取结果进行散点绘制,并提取中心线,结果如图15所示。

图15 剖面下沉曲线Fig.15 Profile subsidence curves

基于沉陷模型提取的走向下沉曲线呈盆地形状,说明在走向方向达到超充分采动,符合开采沉陷一般规律;倾向方向下沉曲线呈山谷对称分布,其中存在部分异常值,经实地探查,造成该异常的主要原因是由于地表人为开挖施工,去除异常值,整体符合开采沉陷地表下沉特征。

4.3 沉陷面域分析

传统开采沉陷监测主要对固定点进行监测,基于沉陷模型不仅能够进行传统点监测、剖面线监测,获取下沉值和下沉折线图。同时也能够进行面域监测分析,进一步提取下沉等值线、下沉范围、地表损害区域图等。

等值线是一种由相同值连接构成的曲线,如等高线、等温线、等深线等。文中采用等值线表示地表下沉值,能够清晰地表示地表下沉情况。等值线的精度一般取决于等值线的间隔值,间隔值越小,图形越精细、等值线越密、表达的效果越好。但是等值线间隔值越小,越容易受异常值和离散数的影响,图形越复杂、异常数据影响越明显。因此,合适的等值线间隔能够完美表达地表下沉情况,同时能够降低异常值的影响。基于工作面下沉值大小和不同下沉值的分布区间,采用多组、分梯度的方法设置等值线间隔,分别选用30,50,100 mm和200 mm间隔值进行下沉等值线绘制,结果如图16所示。

从图16可以看出,不同等值距的下沉等值线所表达的特征详细度不同。较小的等值线区间值对应较详细的沉陷特征,同时也易受噪声和异常值的影响。通过综合比较发现,沉降间隔为100 mm的等值线既可以保留采煤沉陷特征细节,也能够降低噪声和异常值的影响。

传统开采沉陷以边界角确定煤层开采影响范围,但由于西部煤矿地形复杂,不同地形下沉特征有所区别,走向和倾向观测线获取的边界角难以准确描述所有地形条件下的影响边界。因此,基于沉陷模型进行开采边界的提取能够弥补传统测线的不足。但由于地表自然变化(如雨水冲刷、风沙侵蚀等)和无人机LiDAR测量精度的影响,获取10 mm下沉边界存在一定困难。因此,考虑自然因素影响,以100 mm下沉值作为煤层开采引起的下沉边界,提取地表下沉边界,如图17所示。

从图17可以看出,采空区A开采引起的沉陷区域面积140 271 m2,小于表4中计算的下沉面积180 075 m2,造成该结果的原因是由于研究区内存在人工开挖区域和远离采空区的下沉面积,该面积并不是由于煤层开采造成的地表下沉。下沉区域整体偏向采空区B,造成该结果的原因是在2020年11月无人机LiDAR航测前,采空区B已存在,其对地表造成的影响并未结束,因此后续地表下沉范围偏向采空区B。无人机LiDAR监测期间,地表下沉面积与采空区面积比例为2.8∶1,该比例能够为后续煤层开采地表影响范围的确定提供参考。

为了明晰煤层开采造成地表下沉的范围和严重程度,以不同下沉量所占下沉面积为参数,进行地表损害的分等定级,将地表下沉值进行分组,并确定其范围和对应面积,结果如图18所示,可以看出,不同沉降等级区域由内向外逐渐扩展,采空区从中心向外延伸。沉降可分为6级：100～300,300～600,600～900,900～1 200,1 200～1 500 mm和1 500 mm以上,相应的沉降区面积为65 101,29 083,16 354,13 178,11 123 m2和5 432 m2。沉降面积之比为46∶21∶12∶9∶8∶4,能为工作面开采过程中地表建(构)筑物的保护、移民搬迁提供依据。

图18 地表下沉量分级Fig.18 Classification of surface subsidence

5 结 论

1)在煤矿地形复杂区域,采用无人机LiDAR进行煤层开采沉陷监测,针对沉陷区地形特征,对PTIN滤波算法种子点的选取方法进行改进,剔除了最低点噪声,实现无人机LiDAR地面点云的稳健提取。采用Kriging插值算法对地面点云进行内插生成DEM,对格网尺寸进行试验,确定最优的格网尺寸大小(0.1 m),保证DEM误差分布和精度均处于最优,获取较高精度的地面DEM和沉陷模型。

2)无人机LiDAR测量的高程精度可达25～75 mm,受地面噪声、飞行条件、地形地貌的影响,DEM的精度可达25～100 mm,由于相似误差的抵偿,沉陷模型精度可达3～82 mm。因此,研究区基于无人机LiDAR获取的沉陷模型整体精度优于100 mm,平均误差小于47 mm。

3)根据沉陷模型下沉特征进行多方位分析,对沉陷模型进行挖掘和分析,获取地表下沉折线图、下沉等值线图、地表损害分布图、地表下沉面积与采空区比例关系等数据,能够为开采沉陷预计和地表建构筑物保护提供依据。