基于无人机激光雷达技术的开采沉陷监测方法与参数反演

亓立壮， 安士凯， 田超， 周大伟*， 董祥

(1.中国矿业大学环境与测绘学院， 徐州 221116； 2.平安煤炭开采工程技术研究院有限责任公司， 淮南 232033； 3.黑龙江地理信息工程院， 哈尔滨 150081； 4.浙江中测新图地理信息技术有限公司， 湖州 313200)

煤炭资源的开采破坏了矿区的原地质结构，从而易引发地表沉陷、滑坡、泥石流、裂缝等地质灾害。开采沉陷监测就是要对开采引起的地表形变状况进行评估，进而找出沉陷变形规律，从而预测和预防沉陷灾害。目前开采沉陷应用最广泛的方法是概率积分模型，该模型需要事先以实测数据为依据获得预计参数。因此必须在煤矿开采区域建立观测站，收集开采引起的地表形变数据，用来反演概率积分预计参数。

目前开采沉陷监测手段主要有传统观测站测量、合成孔径雷达干涉测量(interferometric synthetic aperture radar，InSAR)、三维激光扫描以及无人机摄影测量等。传统开采沉陷监测手段包括水准仪导线测量、全站仪测量、全球导航卫星系统(global navigation satellite system，GNSS)测量等，虽然凭借其测量精度高、反映沉陷规律准确的优势暂时无法被完全取代，但是存在工作量大、效率低，监测范围局限等诸多问题[1]。INSAR是以雷达天线记录的回波信号为信息源，获取地球表面的三维地形、地表形变和地物特征变化等信息的测量技术，可以对矿区进行大面积、长时间的地表沉陷的监测，具有的全天候、大范围覆盖的监测优势，但由于其原理的局限性，存在大变形区域失相干的缺点，此外高分辨率高精度数据目前比较昂贵[2]。近年来，无人机(unmanned aerial vehical,UAV)摄影测量技术发展迅速，该技术具有易操作、飞行受地形因素影响小和监测周期短等优点，利用无人机影像可以生成矿区地表数字高程模型(digital elevation model，DEM)，两期DEM相减即可得到测区范围内的下沉盆地，但由于受相机分辨率的影响，以及畸变、标定误差和数据处理误差等因素的存在，导致其高程精度还处在分米级水平，而且相片均是地物的表面信息，当植被过于密集时就不能获取足够多的地面信息[3]。

三维激光扫描作为变形监测领域的新手段，被广泛应用于地形形变监测、沉陷滑坡变形监测、桥梁和土石坝变形监测以及基坑变形监测[4-6]等，数据量大且精度较高，但是传统地面三维激光扫描单站采集范围有限且存在地物遮挡，多站搬运步骤烦琐、效率低下[7]。随着无人机技术的迅速发展，特别是大疆等公司的不断推陈出新，无人机设备性能有了极大的提升，使得采用无人机作为飞行平台搭载三维激光扫描仪成为了现实。该方法进行开采沉陷监测具有快速灵活的特点，不需要固定测点，节约成本。由于激光可以穿透植被反射地表信息，并且获得点云密度大，因而该方法可以获得丰富的沉陷盆地观测数据。已得到实际应用，汤伏全等[8]利用低空无人机LiDAR对神榆矿区某工作面进行了开采沉陷监测，对无人机三维激光点云数据建模方法进行了研究；杨凡等[9]在利用无人机激光雷达进行植被高度监测时验证了该技术激光点云的精度可达厘米级；Ren等[10]对无人机平台、不同传感器及其应用领域等方面进行了探讨，构建了一个无人机激光雷达变形监测应用框架。

现以实际工程为例，采用无人机激光雷达技术进行矿区开采沉陷监测及预计参数反演，探讨无人机激光点云数据的处理方法。

1 研究区域概况及数据获取

1.1 研究区域概况

试验区位于内蒙古唐家会矿区，该矿区处在内蒙古自治区准格尔煤田东孔兑普查区的西南部，行政区划隶属准格尔旗大路镇管辖，具体位于准格尔旗人民政府所在地薛家湾镇西北约4 km处，其地理坐标为：111°10′27″E～111°14′34″E，39°52′45″N～ 39°57′22″N。测区地处鄂尔多斯高原，地形多种多样，属典型的侵蚀性丘陵地貌，区域内沟壑纵横，呈树枝状分布，地形起伏较大，植被较为稀疏，主要以低矮植物为主。唐家会煤矿该矿在卫星影像上的位置如图1所示。选择61202工作面上方地表作为观测地点，该工作面尺寸为240 m×1 960 m，平均倾角2°，平均采深538 m，平均采厚19 m，采用长壁后退式全垮综放采煤法开采，设计采放比1∶3。

图1 唐家会矿测区位置Fig.1 Location of survey area in Tangjiahui mine

1.2 UAV激光雷达地表沉陷监测的方法

1.2.1 UAV激光雷达地表沉陷监测原理

UAV激光雷达是以多旋翼无人机平台为载体，一体化集成高精度激光扫描仪、全球定位系统(global positioning system,GPS)、惯性测量单元(inertial measurement unit,IMU)等传感器，通过配备的数据处理和应用软件，快速生成数字地表模型(digital surface model,DSM)、DEM，具有重量轻、携带方便、成果处理效率高的优点。

在工作面开采上方地表划定适当范围的测区，在工作面推进的两个时刻分别对该测区进行观测，并利用采集到的激光点云数据建立两期数字高程模型DEM，这两期DEM相减即为这两个观测时间点之间由于工作面开采所产生的地表下沉盆地。基于概率积分法动态求参原理，利用全盆地的地表下沉数据，再结合工作面开采相关信息，即可获得该地区开采沉陷预计参数。

1.2.2 UAV激光雷达地表沉陷监测技术流程

UAV激光雷达地表沉陷监测流程主要为外业采集、求取下沉盆地、预计参数动态反演。其中外业采集的关键步骤为测区选择、航线规划、实地观测；求取下沉盆地需要先进行数据预处理，经过POS定位定姿系统(position and orientation system,POS)组合解算、数据融合解算之后生成具有较高精度三维坐标的地面激光点云，再经过分类滤波得到地面点云，最后构建两期DEM，DEM相减得到下沉盆地；预计参数动态反演主要是在全盆地均匀选点，根据概率积分原理进行预计参数求取。主要技术流程如图2所示。

图2 无人机激光雷达监测地表沉陷技术流程Fig.2 Technical process of UAV lidar monitoring surface subsidence

1.2.3 数据采集

采用ARS-450i激光测量系统，并利用大疆M600Pro无人机作为飞行平台。该设备集成了激光扫描仪、GPS、惯导等传感器，可同步获取三维激光点云和定位定姿数据，仪器参数如表1、表2所示。

表1 智喙ARS-450i激光测量系统主要参数Table 1 Main parameters of ARS-450i laser measurement system

表2 大疆M600Pro无人机飞行平台主要参数Table 2 Main parameters of DJI M600Pro UAV flight platform

首先根据测区现场实际情况以及无人机续航能力等因素，设计飞行方案。本实验测区面积2.36 km2，分6架次全覆盖扫描。该设备推荐相邻航线重叠率为30%，按照仪器参数设计飞行高度为70 m。分别于2020年8月和2021年3月对上述测区进行了两期数据采集。图3为操作无人机进行外业采集的现场照片，图4为外业采集航线设计示意图。

图3 无人机激光雷达数据采集Fig.3 UAV lidar data acquisition

图4 唐家会测区无人激光雷达航线示意图Fig.4 Route map of unmanned lidar in Tangjiahui survey area

现场实测前需要在测区附近的已知点位(或者现场设置木桩，并测量桩位三维坐标)架设GPS基站，该步骤的目的是采集GPS静态数据与无人机激光雷达采集的POS数据进行组合解算，得到有较高精度的定位定姿数据。确认基站接收数据后，先开启激光雷达静态对齐5 min，目的是对惯性导航系统进行初始位置对准。然后按照飞行方案进行飞行和数据采集，数据采集完毕飞机安全降落后，再次静态对齐5 min，确认数据无误后方可关闭基站。

为了验证无人机激光雷达数据精度，两期实验分别还在测区内均匀采集了若干检查点的三维坐标。检查点主要选在道路等交通便利、无植被的裸露地面，考虑到测区内地形复杂交通不便以及工作效率，实验利用唐家会矿区的连续运行卫星定位服务参考站(continuously operating reference stations,CORS)系统对检查点进行了实时动态(real-time kinematic,RTK)测量。

2 数据处理及沉陷盆地构建

2.1 激光点云数据处理

2.1.1 数据预处理

无人机激光雷达数据采集整理过后，并不能直接用于DEM构建，而是需要经过POS组合解算、数据融合之后，才能获得包含较高精度三维坐标信息的激光点云数据。具体步骤如下。

步骤1激光雷达外业采集的位姿数据包括机载GPS天线数据、惯性导航数据以及里程计数据，POS组合结算是指将上述三者与地面基站收集的同时段GPS静态数据进行联合结算，获得具有较高精度的定位定姿数据。解算方法本文采用了紧组合方式进行组合解算[11]，该方法是运用卡尔曼滤波原理解算信号残差，对GPS与IMU坐标速度之差进行校正，因其耦合程度较高目前得到了最广泛的应用。上述步骤主要在Inertial Explorer软件里实现。

步骤2数据融合解算是指利用步骤1得到的高精度的POS数据，对激光雷达采集的原始三维激光点云数据进行融合处理，得到实验所需坐标系下的具有较高精度三维坐标信息的激光点云数据。因为每个架次的飞行参数一致，所以采集到的数据都是在同一个坐标系下，可以自动拼接为一整个测区的完整点云数据。本实验利用检查点的RTK数据进行四参数求取和坐标转换，最终得到54坐标系下的测区激光点云。

2.1.2 地面点提取

监测区域位于荒漠化草原和草原之间的过渡地带，地面存在人工种植的苔草、沙棘、灌木等植被，此外还存在零星建筑物、车辆、电线杆等，这些会对DEM的构建造成影响，图5为两期原始点云数据。要想获得高精度的地面DEM，就必须对初始点云数据进行滤波处理仅留取地面点云。目前常用的滤波算法有基于布料模拟滤波算法、形态学滤波算法、基于不规则三角网滤波算法等[12-13]。Zhang等[14]提出的布料模拟激光点云滤波算法(cloth simulation filtering，CSF)是将原始点云数据看作一个整体，将其进行倒置，再通过模拟虚拟布料下降的物理过程实现地面点云的过滤。如图6所示，笔者在测区内选取尺寸为100 m的正方形区域作为实验数据，用手动滤波的方式制作标准地面点云，再用多种滤波算法分别进行滤波与之进行对比，具体比较结果如表3所示，在比较多种滤波算法结果后发现布料模拟滤波算法在测区内最为适用，故本实验激光点云数据处理均采用布料模拟算法进行滤波。最终得到两期地面激光点云(图7)。

图5 唐家会激光雷达两期原始激光雷达数据Fig.5 Original lidar data of two phases of Tangjiahui lidar

2.2 DEM构建及精度分析

利用两期地面激光点云数据，以1 m的间隔进行克里金插值，分别构建两个时间点的测区地表DEM。再利用测区内的检查点数据对两期DEM分别进行精度评定。两期DEM如图8所示，DEM精度评定结果如图9所示。通过计算得出第一期DEM高程中误差为0.034 m，第二期DEM高程中误差为0.037 m。

2.3 下沉盆地建立及分析

用第二期DEM减去第一期DEM，得到2020年8月15日—2021年3月13日唐家会测区61202工作面开采引起的地表下沉盆地，如图10所示。测区内最大下沉值约为6.9 m，最大下沉点位于工作面推进约475 m处，下沉盆地东北方向矸石山因堆积升高最大26 m，测区西北部挖土坑深13 m。

(1)从仪器自身角度来看：扫描仪测角、测距本身就存在一定的误差，GPS定位存在轨道误差、钟差、大气延迟、多路径效应等误差，IMU的陀螺仪和加速计也存在自身产生的偏差；机载GPS与IMU之间的相对位置关系也会影响测量系统的定姿定位，虽然厂商会提供相关的参数，但在实际使用中难免存在磕碰、松动等不利因素，导致传感器安置误差的产生。

(2)从数据采集的角度来看：内蒙古地区多风沙，虽然是在仪器允许条件下起飞，但是数据的采集质量依然受到天气条件的影响，这主要是受飞行平台的性能(如稳定性、抗干扰能力等)限制；航高也是影响数据精度的一个重要因素，且在一定范围内激光点云误差和航高呈线性关系，点云密度也会随着航高增加而降低。

图6 不同滤波方式误差分布Fig.6 Error distribution of different filtering methods

表3 不同滤波方式比较Table 3 Comparison of different filtering methods

图7 地面点云Fig.7 Ground point cloud

图8 唐家会激光雷达DEMFig.8 Tangjiahui lidar DEM

图9 唐家会激光雷达DEM精度统计Fig.9 DEM accuracy statistics of Tangjiahui lidar

图10 无人机激光雷达监测结果Fig.10 UAV lidar monitoring results

(3)从数据处理的角度来看：无人机激光雷达数据处理过程比较烦琐，尤其是滤波很难做到完美提取地面点，总是会存在噪声误差；利用地面点构建DEM进行插值时也难免会产生误差。

3 概率积分参数反演及分析

在煤矿开采监测领域中，为了研究复杂的沉降过程以及地表破坏程度，获取沉降值和沉降参数至关重要。预测开采沉陷也是核心研究内容之一[15-17]，对煤田生产和地面防护具有重要意义，可用于确定地层结构是否受采动影响及其相应程度。此外，预测开采沉陷是维修、加固、重建或采取地下措施的基础。预计参数是决定沉陷预测是否准确的关键因素，因此，求取准确的预计参数有重要意义。传统的常规求参方法需要设置观测线，但存在点位较少，需要获取最终下沉，且容易丢失破坏测点、野外劳动量大、监测周期长等诸多问题。无人机激光雷达无需获取完整采动过程的下沉值，只需对监测区域某一时间段或多个时间段内的两期或多期数据进行采集，得到高密度全盆地动态沉陷数据，然后通过建立动态求参模型进行参数求取，可以获取更加精确的沉陷参数。

3.1 概率积分法动态预计原理

在开采沉陷预计方法中，基于随机介质理论的概率积分法在中国的应用最为广泛[15]，该方法是同时考虑了地表沉降的空间位置和时间的相关性，即用概率积分法稳态预测模型乘以时间效应函数，以获取在煤矿开采下不同时间的地表动态下沉值；而概率积分法稳态预测模型是指工作面开采结束后所引起地表沉降的最终状态，可以用来说明影响范围内任意点的预测下沉，其单位开采所造成地表沉陷盆地的计算公式为

(1)

式(1)中：We(x)为x位置的地表下沉值；r为主要影响半径。

由于工作面开采和地表变形是在三维空间内的，则开采一个长为t宽为S的范围所造成的任一点A(x,y)的下沉可表示为

(2)

当工作继续推进，直至开采完毕时，由全部开采导致的A点下沉值的预测公式为

(3)

式(3)中：l为走向长度；L为倾向宽度；W0为达到充分开采的最大下沉值，W0=mqcosα，其中，q为下沉系数；m为采煤厚度；α为煤层倾角；主要影响半径r=H/tanβ，其中，tanβ为主要影响正切，H为开采深度。

因此，若要根据式(1)～式(3)推导出地表其他位置的下沉值，需要确定预测参数。

开采所造成的地表沉陷是一个连续的时空过程，与空间位置、时间、开采速度、覆岩性质等密切相关；其移动和变形可以看作连续的时空函数。动态预测模型也是建立沉陷与空间位置和时间的相关性。用Wt表示t时刻地表点的下沉值，则Wt=W0f(t)，其中f(t)为反映开采沉陷动态变化的时间效应函数。

Knothe时间函数[16]反映了任意时刻的下沉至与最终稳态下沉值之差和该时刻的下沉速度的比例关系，其最终形式为

f(t)=1-e-ct

(4)

式(4)中：c为时间常数。

3.2 基于动态预计的反演方法

运用概率积分法进行开采沉陷预计需要事先获取预计参数，这些参数可以根据实测数据运用数学方法进行反演估计[17]。目前参数反演方法主要有模矢法、线性最小二乘法、遗传算法、粒子群算法等。其中模矢法应用最为广泛，该方法的基本步骤如下。

步骤1构建目标函数(预计值与实际值之差的平方和)，给定一组参数初始值作为基点，代入计算函数值。

步骤2按照一定步长在参数初始值两侧取值，分别代入目标函数得到两组新函数值。

步骤3比较3个函数值大小，选择其中最小的那组参数作为矢点，然后沿坐标系搜索最佳点作为新的基点。

步骤4继续迭代直到求得满足要求的解，该步骤可以在开采沉陷预测预报系统中实现，该软件是由吴侃开发的一款基于概率积分原理的预计模拟软件，可以进行煤矿开采沉陷预计与参数反演。

3.3 求参结果及分析

以2.3节所述整个下沉盆地作为原始数据，运用本章节介绍的反演方法，进行唐家会测区工作面开采沉陷预计参数求取，主要是求取下沉系数q和主要影响正切tanβ。此外，为了对反演结果进行验证，收集唐家会矿区61101工作面的传统观测线数据，该数据是在工作面正上方沿走向和倾向布设两条观测线，采用GNSS技术，观测了该工作面开采前至稳沉的地表移动数据。以相同方法进行参数求取并进行对比，参数求取结果如表4所示。

两组结果对比，运用无人机激光雷达观测的全盆地数据求参结果和传统观测线求取结果差距较小，其中下沉系数相差13%，主要影响正切基本一致。对于下沉系数，传统观测站由于其线状的布站形式，极有可能捕捉不到最大下沉点，而无人机激光雷达则可以获取整个盆地的下沉数据，且2.3小结分析无人机激光雷达求下沉盆地精度约为0.05 m，仅占测区开采沉陷最大值6.9 m的0.7%，因此可以认为该方法反演下沉系数是可靠的；对于主要影响正切，虽然与传统观测线数据求参结果基本一致，但由于中国一般以下沉10 mm为沉陷盆地边界，无人机激光雷达受其精度限制无法捕捉盆地边界，故认为运用该数据求取主要影响正切理论上可能存在一定误差。文献[3]对全盆地求取预计参数的可靠性做了详细分析，当下沉数据中误差占最大下沉值的比例小于7%时，计算的下沉系数、主要影响正切是可靠的。

表4 预计参数Table 4 Estimated parameters

4 讨论

无人机激光雷达应用于开采沉陷监测其最大的优势就在于可以快速高效的矿区地表三维信息的获取和建模。具体表现在：①搭载的激光雷达可以获取较高精度的、高密度的矿区地物激光点云，由于激光具有穿透性，经过滤波处理即可获得大量的地面信息；②以无人机作为飞行平台可以节省大量的人力物力，相当于把测站设在了空中，克服了地面三维激光扫描仪视线遮挡的问题；无人机与三维激光扫描仪的结合可以发挥各自优势，扫描面积更大、获得的信息更全面，且数据处理的流程已经较为成熟，可以形成一整套无人机激光雷达数据采集、处理技术流程。中国煤炭开采正不断向大规模高强度发展，开采活动对地表的影响更加明显，与矿区居民和企业的安全发展密切相关，传统的监测手段存在获取信息量少、工作量大的问题，无法适应煤炭开采的发展趋势，在这种形势下，无人机激光雷达作就成为一种较为合理的技术选择。

分析无人机激光雷达的DEM建模平均高程中误差为3.5 cm，结合现有文献[4-9]可以确定运用该方法进行地表建模高程精度在5 cm之内。未作固定观测站来检核下沉盆地精度，利用误差传播定律得到下沉盆地精度为5.0 cm，事实上两期DEM相减可以消除部分系统误差，下沉盆地精度理论上应该更高。中国一般以下沉10 mm等值线作为开采沉陷盆地边界，因此无人机激光雷达的精度达不到监测沉陷盆地边界及其他小变形的要求，但是与该目标已经十分接近，并且由于其全盆地监测的优势，沉陷数据全面，在求取概率积分预计参数方面较为可靠。综合分析，无人机激光雷达可以在以下方面提高检测精度：一是进一步提高激光雷达的精密度，仪器本身为多传感器集成，激光扫描仪、惯导、GPS接收机均存在一定误差，三者的位置关系是否固定也影响系统定姿定位；二是进一步提高无人机性能，无人机的抗风能力和稳定性影响扫描姿态进而影响激光点云数据的精度；三是进一步改进数据处理方法，改进滤波算法以及DEM建模方法以获取更为精确的地面信息。

无人机除搭载激光雷达之外，还能搭载相机进行摄影测量。无人机摄影测量与无人机激光雷达相比各有优势与不足[18-21]，分析如下。

(1)精度方面，无人机摄影测量由于其相机性能、相片几何畸变、重叠度、飞行平台稳定性等因素的制约，目前还处于分米级的阶段，与无人机激光雷达还存在一定差距，但二者都能反映全盆地的沉陷信息，在进行预计参数反演尤其是下沉系数的求取时都具有相当高的可靠性。

(2)数据方面，无人机摄影测量获取的实时图像数据，首先由于相片本身是地物表面信息，因此当监测区域植被过于密集时，难以捕捉足够的地物信息，滤波过后存在较多空白区域；其次当地形起伏较大时，会造成连续的影像之间重叠度损失，实际工程常常被迫设置高重叠度从而增加工作任务，再者当大气质量较差或光线较暗时，相机难以获取高清晰度的影像信息，而激光雷达是主动发射激光接受反射信息，以其强大的穿透性可以穿过大气、植被直接反射地表信息，亦不受环境亮度的影响。

(3)工作效率方面，受点云密度的制约，搭载激光雷达的无人机往往比搭载相机的无人机飞行速度慢；激光雷达设备重量大，对飞机性能尤其是电池续航能力要求更高，实际应用中往往需要频繁起降更换电池；激光点云反射距离有限，使得激光雷达无人机飞行高度一般要低于摄影测量无人机；摄影测量需要在测区内布设若干像控点，并且随测区增大所需点数亦增大，激光扫描则不需要像控点的布设，因而节省了一部分人力物力。

综上所述，当需要快速大面积监测作业，且对精度要求较低时，选择无人机摄影测量更为合适，而测量范围相对较小，且要求厘米级精度时，无人机机光雷达则更具优势。随着无人机技术的继续发展和激光测量系统的进一步改进，无人机激光雷达在数据采集效率和精度上还有极大的进步空间，长远来看，该技术在开采沉陷监测领域更具发展优势和潜力。

5 结论

(1)目前主流的开采沉陷监测技术均存在着精度与效率两者不可兼得的问题，无人机激光雷达充分集合了三维激光扫描数据量大精度高与无人机操作方便工作迅速的优点，既能快速获取大范围的地表激光点云数据，又能有效节约时间与劳动成本，尤其适合地形复杂多样的西部矿区进行开采沉陷监测。

(2)对原始激光数据分别进行组合解算、数据融合、分类滤波之后，建立两期测区数字高程模型DEM，通过与检查点对比得到中误差分别为0.034 m和0.037 m；两期DEM相减得到观测时段内的开采沉陷盆地，最大下沉值约为6.9 m，根据误差传播定律得到下沉盆地中误差约为0.050 m。

(3)对下沉盆地进行均匀取点作为原始观测数据，采取基于动态预计的反演方法，求取该矿区概率积分预计参数，其中下沉系数主要影响正切比较可靠。收集同一矿区内临近工作面的传统观测线数据进行参数反演，得到的预计参数结果与无人机激光雷达全盆地求参基本一致。综上所述，利用无人机雷达进行开采沉降监测是可行的，利用无人机激光雷达监测所得完整沉陷盆地进行预计参数反演是可靠的。