西部矿区是我国主要煤炭生产基地，大规模、高强度地下采煤已导致地表大范围塌陷及各种衍生灾害[1-3]。在矿区开采沉陷的实地监测中，现有的GNSS、水准等大地测量方法及InSAR等遥感监测方法均存在一定的局限性。前者采用走向和倾向设站进行常规观测，难以获取整个地表移动盆地的变形特征，且野外工作量大[4];后者利用SAR影像进行干涉处理，易产生失相干现象，无法得到地表沉陷盆地内大梯度变形信息，一些学者采用偏移量追踪技术(offset tracing)虽然有助于解决上述问题，但所需的高分辨率SAR影像数据获取成本高，难以推广应用[5]。
近年来，低空无人机遥感技术已在矿区环境监测等领域开展多方面应用[6-10]，如将无人机遥感技术用于高潜水位矿区沉陷地受损作物植被指数反演及露天矿边坡的精细测绘与建模等。而无人机激光扫描(Light Detection and Ranging，LiDAR)作为一种低空激光遥感技术，通过对地面进行周期性扫描，可获取煤矿开采沉陷区高精度、高分辨率的点云数据，经过滤波和地面点插值形成DEM，通过对多期DEM进行叠加求差，可生成地表高程变化模型，即沉陷DEM。目前该技术在矿区沉陷监测中虽有一定应用，但尚处于探索阶段[11-13]。在激光点云滤波、插值、DEM叠加过程中，不可避免地存在噪声和建模误差，使得沉陷DEM包含显著的误差，一些学者通过各种滤波算法和模型去噪方法来改善上述沉陷模型的精度[12,14]，但仍存在以下主要问题:① 在复杂地理环境下无人机激光扫描系统获取的点云数据精度及其可靠性缺乏实测数据验证;② 在西部风积沙矿区植被覆盖较低和地形起伏条件下，针对现有主流的点云滤波和DEM插值算法的适用性及DEM建模的精度还缺乏定量评价;③ 对沉陷DEM进行去噪时，没有结合地表沉陷盆地本身的分布特征进行算法改进，造成模型去噪效果较差，导致地表沉陷建模的精度难以满足采煤沉陷监测的基本要求。
因此，笔者针对西部矿区地形起伏而植被较少的地理环境，利用机载LiDAR技术开展矿区沉陷监测建模实验研究，以榆神矿区金鸡滩煤矿某工作面地表为实验区，利用无人机LiDAR获取两期4组点云数据开展沉陷建模研究，对现有主流的点云滤波、插值算法的适应性进行对比选优，在获取初始沉陷DEM的基础上，结合沉陷盆地特征及模型误差分布特性，基于小波阈值去噪原理提出沉陷DEM的多尺度去噪优化方案，通过实测数据验证沉陷模型的去噪效果及实际精度。同时，结合沉陷模型边缘的倾斜变形特征，提出了沉陷模型边界的确定方法，为机载LiDAR在西部矿区采煤沉陷监测中的推广应用提供技术支撑。
1 实验区及数据概况
实验区位于陕北榆神矿区，地处毛乌素沙漠边缘，属于黄土高原的过渡地带，地形起伏较明显，植被覆盖度低且以低矮沙生植被为主[15]。实验区地表无水系通过，局部有水坑。井下工作面平均采高10.5 m，宽度300 m，总长度4 500 m，煤层倾角0.5°，平均埋深260 m。采用SZT-R250型无人机LiDAR系统对工作面开切眼一侧的沉陷区域进行两期地面扫描。外业扫描分别于工作面回采至450 m(2019-04-02)、980 m(2019-06-01)时进行。每次扫描矩形区域从开切眼一侧煤柱上方地表(距开切眼约200 m)沿工作面走向延伸1 200 m，沿工作面倾向边界向外侧各延伸250 m，去除边缘点云后扫描实验区长1 400 m、宽800 m。为了进行实验数据验证，在工作面地表布设常规观测站，采用GPS-RTK和四等水准测量分别进行测点平面坐标和高程监测。走向观测线沿工作面中央地表布设测点A1～A47，倾向观测线距开切眼510 m，布设测点B1～B49，与走向线正交于A45号点，测点间距15 m，2条观测线外共布设8个控制点K1～K8。工作面、地表观测站及实验区相对位置关系如图1所示。
图1 工作面、观测站及实验区相对位置关系
Fig.1 Diagram of the relative position of working face,observation station and experimental area
根据开采进度和常规监测结果，该工作面于2019-01-10开始回采。2019-01-31推进至120 m时，地表开始出现采动下沉;在3月26日工作面推进至385 m时，距开切眼150 m的地表测点A21达到最大下沉量4.849 m，已经接近充分下沉状态，但尚未稳定。两期扫描数据在工作面推进450 m(4月2日)至980 m(6月1日)期间进行外业采集，此阶段地表在测点A45两侧出现下沉平底区，基本达到充分采动状态，但由于放顶煤开采条件下采高变化较大，加上地表移动尚未稳定，两次扫描期间地表下沉盆地的平底区存在明显的高低起伏，平底区实测最大下沉量相差超过1 m。
每期激光扫描均重复一次，获得独立的两组数据。无人机飞行高度约70 m，使用GNSS和惯导系统INS进行定位，其中激光发射频率为100 kHz，角分辨率为0.001°，扫描数据精度为5 cm，在经过解算与转换处理后，生成标准点云数据文件。
两期扫描数据的点云平均密度约30点/m2，滤波后标准地面点云平均密度远超过10点/m2，点云平均间距小于0.5 m，按照0.5 m分辨率构建DEM时，可保证各栅格中有点云数据覆盖，因而点云数据满足0.5 m格网DEM的建模要求。
外业扫描时不同架次的坐标及高程基准与常规地表移动观测站的控制点保持一致。每一次作业的GNSS校准点为相同的实测控制点(图1)，保证了各批次点云数据之间统一基准。由于无人机惯导系统采用的INSCGCS2000为以ITRF97参考框架为基准2000历元下的坐标系统，不会随着时间而发生变化。因此只要扫描时基站已知点坐标、作业流程和航带路线保持不变，即可保证两期点云数据之间坐标基准的统一性。
为了减小各架次扫描点云的坐标误差，通过调整点云解算参数来实现多批次点云数据的配准。具体借助TerraSolid软件的TerraMatch模块完成，该模块通过处理GPS定位、INS定向以及扫描镜误差，解决激光扫描仪和惯性测量单元之间的比例因子和错位问题，实现各航带批次的点云配准。首先依次进行POS解算、点云解算，然后导入点云和航迹线数据，通过点云粗分类，在平地、建筑、斜坡处搜索连接线，即代表局部小范围内点云平均坐标位置的线段，通过各架次重叠部分的连接线来校准全区域点云数据。其中，平面区域的连接线主要用于各批次点云高程的平差配准，而建筑屋顶、斜坡顶的连接线主要用于点云平面坐标的平差配准，根据搜索到的连接线解算机载LiDAR扫描系统的航向角、横滚角和俯仰角、扫描镜误差修正值，并用于修正当前加载的所有连接线及点云数据，再用修正后的参数重新解算点云，迭代校正，以使各次数据间出现分层、偏移的连接线在x,y,z三个方向上的误差收敛接近为0。两期点云数据的平差配准同样参照上述原理，但由于存在沉陷变形区域，对第2期数据生成的连接线需进行编辑，仅保留稳定区部分，通过迭代实现两期点云数据的平差配准。
利用观测站控制点来验证各次点云配准效果，每次扫描前在各控制点上放置800 mm×800 mm的矩形反射平板，其中心及平板高程与控制点标志相同。上述反射平板上的点云不含非地面点噪声及复杂地形影响，将平板上点云高程的平均值作为该控制点的扫描高程，对配准后点云数据的高程进行精度验证，8个控制点的扫描高程误差均不超过0.05 m。
2 沉陷DEM获取及误差分析
2.1 点云插值及滤波算法对比优选
以所构建DEM的高程精度为出发点，在实验区均匀选取A,B,C,D四个200 m×200 m的实验区块，如图1所示，对比了专业化数字高程模型插值(ANUDEM)、反距离权重插值(IDW)、克里金插值(Kriging)、自然邻域插值(NN)、样条函数插值(Spline)等主流插值算法[16]的插值精度。通过手动编辑提取的标准地面点云进行验证，各算法误差均值、标准差、拟合优度对比如图2所示。
图2 不同插值算法对应DEM误差统计
Fig.2 DEM error statistics of different interpolation algorithms
经对比分析，专业化数字高程模型插值算法在实验区表现出最好的适用性，最小误差均值为0.045 m，误差标准差为0.046 m，最大拟合优度达99.95%。这与该算法的嵌套式多分辨率迭代处理相关，在不同分辨率下，拟合点处有就近点则采用就近点，没有就近点则通过高斯-塞德尔迭代逐次求出拟合值，该算法同时兼顾局部处理及全局处理，迭代次数的调整可优化插值结果的精度，多次迭代不会使其时间效率最高，但算法局部处理的高效性优化了其运行效率。
同时对比分析了基于高程阈值的滤波(ETEW)[17]、多尺度曲率滤波(MCC)[18]、基于坡度阈值的滤波(MLS)[19]、渐进形态学滤波(PM)[20]、三角网渐进加密滤波(PTDF)[21]等点云滤波算法对应的DEM误差。同样通过手动编辑提取的标准地面点云进行验证，利用各实验区块点云经各算法滤波、专业化数字高程模型插值生成的DEM误差指标对比如图3所示。
图3 不同滤波算法对应DEM误差统计
Fig.3 DEM error statistics of different filtering algorithms
结果表明，三角网渐进加密滤波在该区域滤波地面点插值的DEM精度最高，误差均值为0.007 m，标准差最小为0.013 m，拟合优度最大达99.998%。该算法滤波主要通过迭代距离及迭代角度两个参数协调控制加密三角网以选取地面点，能有效减少地形跳跃及去除研究区域的低矮植被。基于该算法构建的DEM误差均值和标准差最小，能有效去除研究区域的低矮植被。除了三角网渐进加密滤波算法，其他方法在阈值设置上缺少自适应性，多种算法的坡度阈值多为固定值，使得算法取得较差的滤波效果，George Vosselman等[19]也指出在地形复杂多变的区域，对整个实验区域设定统一的坡度阈值明显是不合理的，坡度阈值不应作为固定常量，而应随实际地形的变化而变化。
2.2 初始沉陷DEM获取
通过点云的滤波、插值及叠加处理，获取工作面开采时段的地表沉陷DEM。首先，对4月、6月两期数据(各取1组)进行滤波处理，采用三角网渐进加密滤波算法，结合ISPRS提出的点分类误差评定方法确定该算法的最优参数，滤波地面点平均密度约为10点/m2。然后，基于滤波结果进行专业化数字高程模型插值，插值结果像元大小最佳尺寸为0.5 m，最后将6月份与4月份两期DEM进行叠加，得到初始沉陷DEM。由于现有的滤波算法并不能完全准确的分离出所有的地面点，且点云经插值得到DEM又进一步引入了插值算法误差，导致所得到的初始沉陷不可避免地存在明显的误差。初始沉陷DEM三维视图如图4所示。图4中将高程方向拉伸50倍显示以突出模型误差的影响，主要误差区域用圆圈表示，其中黑色圈主要为空间位置偏移误差，绿色圈主要为未去除非地面点误差、蓝色圈主要为水面引起的插值误差。
图4 初始沉陷DEM的三维视图
Fig.4 3D view of initial subsidence DEM
2.3 误差特征及其处理
综上可知，利用对比确定的最优点云滤波和插值算法所构建的沉陷DEM仍然存在明显误差，分析原始点云数据及图4中的误差特征可知，沉陷DEM误差主要来源以下几方面:
(1)无人机LiDAR扫描点云的空间位置误差引起的沉陷模型误差。通过对每期两次独立扫描数据进行精度验证表明，实验条件下扫描点云的平面和高程相对误差不超过0.05 m，而沉陷建模是针对多期点云数据进行高程求差，如果基于同一位置的多期点云直接计算沉陷量(高差)时，其沉陷值的精度与目前GPS-RTK观测结果基本相当。但在利用点云内插生成DEM并进行叠加的过程中，由于扫描点云的密度所限，在地貌起伏较大的情况下，即使是微小的平面误差也可导致显著的DEM叠加误差。通过实验数据对比表明，在建筑物、人工堆挖地、陡坎陡坡边缘处因很小的平面误差可使差值DEM产生明显的误差[12]，主要分布如图4所示黑色圈注。因此，复杂地貌条件下点云本身的位置误差对沉陷DEM有较大影响。
为了减小这种点云平面位置误差造成的沉陷DEM误差，应尽量降低无人机飞行高度以增大扫描点云的密度，从而提高DEM分辨率，同时对两期点云构建的栅格DEM进行再次配准。针对第2期点云数据在沉陷区以外的稳定区域，选取多处高程变化显著的特征地貌单元(如建筑物、陡坎、台阶、小山包等)，利用TerraSolid软件分别沿坐标系x,y方向绘制点云剖面线。在小范围内及点云密度较大情况下，两期数据绘出的对应特征地貌单元剖面应具有相同的曲线形态，两曲线在x,y及高程z方向产生一定的偏差，利用其平均差值(dx，dy，dz)将第2期点云配准至第1期点云基准下。具体算法实现可借助TerraSolid软件的TerraMatch模块完成，也可采用与遥感图像或数码影像几何配准类似的方法实现。
2次配准的效果取决于所选地貌特征单元的高程复杂度及点云密度。由于激光扫描点云平面误差主要取决于无人机系统实时组合导航定位的稳定性，所采集的点云数据误差仍以随机误差为主，实验表明，上述点云相对配准的效果还有待进一步改善。应该指出，近年来随着无人机LiDAR后处理POS技术的发展，通过事后改正POS定位结果可显著提高系统的姿态精度和整组数据的一致性与数据平滑性，减小点云数据的随机误差，进一步通过多期点云的相对配准可减小点云的系统性误差，有效提高复杂地貌条件下所构建的沉陷DEM精度。
(2)滤波算法未去除的非地面点引起模型误差。由于点云滤波算法的局限性，一些植被、低矮地物等非地面点不能完全去除，加上不同扫描时期植被生长状况有所差异，未去除的植被点分布位置并不重叠，导致叠加后的沉陷DEM出现明显噪声，主要分布如图4所示绿色圈注。该类误差主要通过沉陷DEM的去噪算法进行处理。
(3)采用点云内插算法生成格网DEM引起的模型误差。因点云密度有限且多期点云位置分布有所不同，针对多期点云分别进行内插生成格网DEM时，会造成多次精度损失，而点云密度、格网尺寸和地形复杂程度均对点云内插精度产生显著影响。此类误差分布在整个区域，尤其在地形陡峭、植被覆盖大及地面点云密度较低的区域，所产生的模型误差较大。该类误差可通过增加点云密度和精度、改进点云滤波和插值算法及沉陷模型去噪等方法来减小。
(4)水域覆盖范围变化导致的沉陷模型误差。研究区域包含几处明显的水域，沉陷前、后水域范围变化明显，沉陷后的水域范围明显增大。由于激光雷达在水域范围不能获取点云数据，利用水域边缘点参与DEM建模时，当两期数据水域范围明显变化时，会导致叠生成的差值DEM出现较大误差，主要分布如图4所示蓝色圈注。该类误差可通过合理确定点云数据处理边界及优化点云插值算法来剔除。
上述4类误差因素对沉陷DEM模型误差的影响有所不同。在进行模型去噪时，可分别对每期数据生成的地形DEM去噪，也可直接对叠加后生成的沉陷DEM去噪。但地形DEM所反映的只是地貌形态特征，而直接对沉陷DEM去噪不仅能去除模型叠加导致的各项误差，还能利用沉陷盆地本身的分布特征和非沉陷稳定区DEM高差为零的先验知识，有利于优化去噪算法及及其参数选取。因此本文直接对叠加后的沉陷DEM进行去噪。
3 基于小波的沉陷DEM去噪
3.1 小波阈值去噪原理
由于沉陷DEM的误差多为突变误差，且服从高斯分布。而随着小波理论日益成熟及大量研究发现，其优良的时频局部化特征可以成功去除信号数据中的局部高频突变噪声[22-23]。利用小波(即有限长或快速衰减的振荡波形)的缩放和平移，实现对原始信号的良好匹配，再对基于小波展开匹配的系数做滤波处理，即可有效分离出信号的精细及粗糙成分。小波分析在图像去噪、分割及压缩等领域得到了广泛应用[24-26]，尤其在高斯噪声的滤除方面收到了很好的效果。基于此，本文选取小波去噪的方法实现矿区沉陷DEM的去噪优化。
小波去噪问题的本质是函数逼近问题，即如何通过小波基函数伸缩、平移获取的一系列函数空间，根据既定的衡量准则，确定对原始数据信号的最佳逼近，以区分有效信号和噪声信号。目前，小波去噪的基本方法有:模极大值重构去噪、基于小波系数的相关性去噪、小波变换阈值去噪、平移不变量小波去噪，本文选择应用较广的小波变换阈值去噪的方法对矿区沉陷DEM进行处理:首先，选择合适的小波基函数，确定分解层次，对原始图像进行小波变换，以分离出多个尺度层次的图像高频、低频信息，其次，设置阈值并依据阈值函数对变换后的小波系数进行阈值处理，将处理后的各层次高频及最后一个层次的低频信息进行小波逆变换以重构去噪后的图像。其中，分解层次为小波阈值去噪中的关键参数之一，不同信号或图像数据都存在一个去噪效果最好的分解层数，即最佳分解尺度。小波分解层数会对数据去噪的效果产生很大影响，当分解层数过多，且对各层小波空间的系数都进行阈值处理时，会损失大量原始数据的有效信息，同时增加运算量、降低处理效率;分解层数过少时，则数据噪声不能有效去除;阈值的选择同样是小波阈值去噪中的关键问题，阈值过大或过小都会影响去噪效果。如果阈值设置过小，处理结果中仍然会存在大量噪声;阈值设置过大又会造成原始数据有效信息的缺失，使获取的结果过于平滑。
3.2 矿区沉陷DEM去噪方案
针对沉陷DEM去噪，利用Matlab对小波阈值去噪的关键参数进行测试对比。其中，分解尺度层次及阈值通过两个向量来存储表示，n为尺度向量，p为阈值向量，为保证沉陷模型的平滑性，采用软阈值函数处理，依据分解层数较少时，噪声占比较大，阈值设为随分解层数的增加而减小。通过统一并变换尺度向量n及阈值向量p，以沉陷DEM处理后的三维可视化结果为辅助，把握非沉陷区的去噪标准，即非沉陷区没有明显突变且接近平地状态，测试对比Matlab中Haar,Daubechies,Biorthogonal,Coiflets,Symlets等各系列小波函数的去噪结果，选出去噪效果最好的小波函数。这里需要说明，非沉陷区的突变不是由开采沉陷引起的，而是由非地面点等其他非沉陷因素导致的，沉陷DEM中非沉陷区的沉陷值应接近为零，当这些突变噪声被消除后，可认为基本消除平滑了整个数据的非沉陷噪声。
由于实验区域两期扫描时间段的动态沉陷在倾向上不存在盆地平底区，沉陷盆地呈锥子型凹槽，若采用与非沉陷区相同的尺度及阈值处理，其底部起伏会被视为突变噪声造成过度平滑，故沉陷区与非沉陷区的小波去噪应取不同尺度的参数。因此，针对不存在超充分采动平底区域的沉陷DEM去噪处理方案为:将沉陷区与非沉陷区分别进行处理。先利用非沉陷区下沉值为0这一先验条件，针对全区域数据进行多层次小波分解;再对沉陷区采用低层次小波分解;然后将沉陷区处理结果与全区域处理结果做镶嵌处理，重叠区域取对应的最小下沉值。该方案中多层次小波分解以非沉陷区没有突变误差为基准，可去除整个区域由非开采沉陷因素引起的噪声和突变误差;低层次小波分解则保留了沉陷区真实下沉的细节信息，对沉陷区有一定的平滑去噪效果，但不会造成多层次分解后的过度平滑现象。沉陷DEM小波去噪处理流程如图5所示。
图5 矿区沉陷DEM小波去噪流程
Fig.5 Process of wavelet denoising of mining subsidence DEM
针对全区域的多层次分解采用bior5.5小波函数，进行10个层次的分解去噪，其中，各分解层次对应的小波系数阈值向量为:p=[25,23,20,18,15,13,10,8,5,3]，bior5.5小波具有双正交性、紧支撑性，支撑长度为11，消失矩为4阶，有利于有效信息的集中，较高阶的消失矩有利于数据去噪，且其非对称性适用于去除非沉陷因素引入的噪声信息[27]。
针对沉陷区域的低层次小波分解采用coif5小波函数，进行5个层次的分解去噪，各分解层次对应的小波系数阈值向量为:p=[10,8,6,4,2]，coif5小波具有双正交性、紧支撑性、近似对称性，其支撑长度为29，消失矩为10阶，经过实验验证，这一支撑长度有利于矿区下沉信息的提取，其高阶消失矩更利于沉陷区的去噪，小波的正交性则保证了对原始数据具有较好的平滑性能[27]。
3.3 矿区沉陷DEM去噪实验
采用上述优化方案进行小波阈值去噪。通过实验选择合适的小波函数、分解层次及小波系数阈值，并基于实测数据对比来验证其去噪效果。
外业采集的两期4组点云数据中，每期的2次扫描时间间隔在4 h以内，故同一期的两组数据之间可视为不存在沉陷量。实验首先将各组数据通过三角网渐进加密进行点云滤波、专业化数字高程模型作插值处理，获取4组DEM数据;再分别将4月、6月中第1组、第2组DEM对应进行叠加(求差)，得到两组含噪声的初始沉陷DEM。经小波测试分析，采用前述去噪方案用相同参数对两组沉陷DEM分别去噪，得到去噪结果1、结果2。
去噪后两组沉陷DEM应具有相同分布特征。将两者进行叠加分析，并与地表观测站的实测数据进行对比，以验证小波去噪效果。由于研究区4月、6月的水域覆盖范围存在变化，最大水域面长达200 m，宽50 m，鉴于大面积空白区点云插值得到的沉陷结果缺乏可靠性，实验中去除了覆水区域。
3.4 去噪结果分析与讨论
为了评价上述实验结果，通过两组沉陷DEM对比、与实测数据对比以及主断面处理前、后进行对比，来论证本文基于小波阈值的沉陷DEM去噪的实际效果。
3.4.1 去噪后两组沉陷DEM的相对精度
经过上述实验得到研究区的两组沉陷DEM去噪结果1与去噪结果2。显然，两组去噪结果之间差异性可以衡量沉陷DEM去噪结果的相对精度和可靠性。利用ArcGIS对格网下沉偏差进行统计表明，两组去噪结果之间的互差均值为0.034 m，互差标准差为0.031 m，均方根误差为0.046 m，各项误差指标均在5 cm以内，两个沉陷DEM去噪结果的差异极小。与激光点云本身的采集精度相比，经过去噪处理后沉陷DEM中下沉量的相对精度并未降低。图6展示了去噪结果1的三维视图。与图4的初始沉陷DEM三维视图对比，可从视觉上看出去噪效果非常明显。图6中沉陷DEM在非沉陷区域不再包含突变噪声信息，沉陷盆地以外的稳定区域内下沉值近乎为0;而沉陷区在去除突变噪声的同时，仍然较好地保留了地表沉陷盆地真实的细节特征，相比于常规的预计模型所描述的沉陷盆地，利用激光扫描点云构建的沉陷模型准确地反映了采煤沉陷盆地真实的非连续、非平缓特征。
图6 去噪结果的三维视图
Fig.6 3D view of the denoising result
3.4.2 去噪前、后沉陷DEM主断面下沉对比
主断面分析是开采沉陷分析的重要方式。为了更直观地反映本文方法进行沉陷模型去噪的效果，基于处理前、后的沉陷DEM分别沿下沉盆地走向和倾向主断面提取下沉剖面数据，以走向、倾向主断面交点为坐标原点，其下沉曲线如图7所示。
图7中，由于走向上存在4处覆水区域，导致下沉曲线不连续。由图7可知，去噪后取得了较好的平滑效果，也很好地保留了沉陷盆地的下沉细节特征，沿走向主断面上下沉盆地的充分采动区(盆地底部区域)在进行了多层次分解小波去噪及镶嵌处理后，仍然存在起伏波动情况，经实地测量验证了该曲线特征与实际相符。沿倾向主断面下沉盆地不存在平地区域，去噪后下沉曲线取得了较好的平滑效果并保留了沉陷底部有效的沉陷信息。
3.4.3 与实测数据对比
利用该工作面地表移动观测站的RTK及水准观测数据，将2组沉陷DEM去噪处理前、后的主断面下沉曲线与对应时间段两期实测数据获取的下沉曲线进行对比，如图8所示。
由图8可知，去噪后走向和倾向下沉曲线与实测曲线拟合效果更好，特别是盆地底部区域更符合实际情况。在与实测数据进行对比时发现，部分测点标志桩高出或低于地面0.05～0.10 m，激光扫描时受点云密度所限，可能并未扫描到标志桩而是周围的地面点，导致地面点云高程与实测的标志桩高程不相等。
进一步分析发现，由于实测数据是不同时期基于同一标志桩采集的三维坐标，在此期间标志桩产生下沉并伴随着发生平面位移，在两次扫描期间实测标志桩的最大平面位移接近2 m，而实测的下沉量是测桩原始高程与平面坐标变化后的测桩高程之差，实质上是不同平面位置的高程差值;而沉陷DEM下沉量是两期点云DEM中相同平面位置上两个栅格的高差值。当测桩处周围地形起伏大或测桩的平面位移较大时，会造成明显的对比误差。实验范围内的实地测桩数分别为走向线47点，倾向线49点，共计96点。其中有23个测桩高于或低于地面超过0.05 m，5个测桩位于陡坡处且水平位移量超过1 m，在进行实测数据对比时，舍弃上述测桩点，其余68个实测数据参与误差统计。
图7 两组沉陷DEM处理前后的主断面对比
Fig.7 Comparison of the main sections of two sets of subsidence DEM before and after denoising
图8 沉陷DEM1和DEM2去噪前、后及实测下沉曲线对比
Fig.8 Subsidence curves of the two subsidence models before and after denoising compared with the measured data
沉陷DEM1,2去噪前、后与实测数据间的误差均值、标准差及均方根误差统计结果见表1。由表1可知，采用本文去噪方案处理后的两组沉陷DEM误差指标均明显减小。其中，去噪结果1与RTK实测数据的均方根误差由去噪前的0.099 m降低到0.044 m;去噪结果2的均方根误差由去噪前的0.101 m降低为0.046 m，而且去噪结果1,2对应的各项误差指标都基本一致，也表明本文针对点云数据处理和沉陷模型去噪的结果具有较高的可靠性。
表1 沉陷DEM去噪前、后误差指标统计
Table 1 Error index statistics before and after denoising of the two subsidence DEMsm
项目误差均值处理前处理后标准差处理前处理后均方根误差处理前处理后沉陷DEM1 0.0730.0350.0980.0420.0990.044沉陷DEM2 0.0700.0360.0940.0430.1010.046
经统计，去噪后沉陷DEM误差小于10 mm的点数占比由处理前的8%提升至15%。对于沉陷DEM中达到毫米级精度的这部分栅格，可视为准确的沉陷信息，可进一步挖掘利用。
4 沉陷DEM边界提取
按照现行的开采沉陷规范要求，以地表下沉10 mm等值线作为开采沉陷边界，按照去噪后的沉陷DEM提取下沉10 mm的边界线，如图9所示中橙色等值线。由于沉陷DEM的下沉标准差为0.046 m，其误差明显大于10 mm，导致所得到的沉陷DEM在边缘区域下沉量出现正、负交替，故依据10 mm下沉量确定的沉陷盆地边界非常不规则，不符合开采沉陷的基本特征。若取下沉标准差的2倍作为极限误差，可认为沉陷DEM中下沉量超过-0.091 6 m等值线所圈定的范围属于确定的沉陷区，如图9所示中红色实线所示。图9中该等值线呈现明显的不规则形状并延伸到盆地外的稳定区域，这是由于沉陷DEM的不确定性误差所致。沉陷边界必然位于该下沉等值线以外的边缘区域。根据矿区地表沉陷盆地的三维空间特征，地表下沉盆地边缘的坡度由内向外呈现由大到小的确定性趋势，最终接近为0，理论上可取下沉坡度(实质为地表倾斜变形)为0°的等值线作为沉陷边界线，考虑到下沉量误差也会引起下沉坡度(倾斜)产生一定的误差，可根据沉陷DEM下沉标准差的大小，结合常规地表移动观测站中测点之间的距离，来确定下沉坡度(即倾斜变形)的误差值，以此作为判定沉陷边界的下沉坡度临界值。实验区条件下常规地表观测站测点的合理间距应为15 m，由DEM下沉标准差0.046 m计算出下沉坡度的临界值约为0.2°。
图9 基于下沉与坡度等值线提取的沉陷边界
Fig.9 Subsidence boundary extraction based on subsidence and slope contour lines
为了探讨这一判别方法的可行性，按照下沉坡度0.2°绘出对应的等值线，如图9黑色实线所示。该下沉坡度等值线不仅将确定的沉陷区域几乎包含在内，同时与10 mm下沉等值线的最内边界大部分相吻合，并具有较好的平滑性。利用走向和倾向线该时间段的实测高程值，确定主断面上地表下沉实测边界分别在A24,B10,B40附近，如图9所示。可见按本文方法提取的下沉边界与实测边界点基本一致。
应该指出，由于沉陷DEM边缘区域下沉误差具有随机性和正、负抵偿特性，可沿下沉盆地主断面方向选取若干相邻格网形成栅格组，计算栅格组的下沉均值，通过移动搜索出下沉均值接近为0或者10 mm所对应的栅格组，其中心位置便可视为沉陷边界的合理位置。同时，针对获取的多期沉陷DEM，将模型边缘非沉陷区内多期DEM叠加并取下沉均值，也能有效提高沉陷DEM边界的确定性。
5 结 论
(1)在西部榆神矿区地形起伏、植被较少的风积沙覆盖矿区，采用机载LiDAR系统进行地表沉陷监测时，按照现有的常规点云滤波、内插算法所生成沉陷DEM精度不足。实验条件下无人机激光扫描点云数据本身的测量精度可达0.03～0.05 m，而采用扫描数据按照现有技术流程生成沉陷DEM时，受点云平面误差、非地面点噪声、点云密度不足或复杂地貌环境带来的高程内插及DEM叠加误差等影响，所生成的沉陷DEM存在显著的噪声。实验结果表明，初始沉陷模型下沉量的均方根误差达0.1 m左右，在沉陷区和沉陷盆地外的稳定区域均存在明显的粗差，难以达到开采沉陷精准建模的基本要求。
(2)针对初始沉陷模型及其误差特征所提出的基于小波阈值的多尺度模型去噪方案，能够有效地去除沉陷模型噪声。针对沉陷盆地和非沉陷区域选用不同的小波参数去噪，不仅能有效地去除沉陷模型尤其是盆地边缘及非沉陷区域的随机噪声，也能很好地保留沉陷区尤其是盆地中央平底区下沉起伏变化的细节特征，展示出利用激光点云构建高分辨率矿区沉陷模型的优势。实验结果表明，经过本文方案去噪后沉陷DEM的各项精度指标都有明显改善，2组沉陷DEM的相对均方根误差为0.046 m，与常规实测沉陷量相比，沉陷DEM主断面上测点下沉的均方根误差为0.045 m，毫米级误差的栅格数量占比达15%。2组沉陷DEM各项误差指标十分接近，验证了点云数据的可靠性和本文去噪方案的有效性。
(3)根据矿区沉陷盆地三维特征，基于沉陷DEM误差所提出的沉陷边界提取方法具有一定的合理性。地表下沉盆地边缘区的坡度(倾斜变形)由内向外呈变小趋势，考虑到沉陷DEM误差带来的不确定性，利用下沉量的标准差给出了提取沉陷边界的临界下沉坡度值。本文实验条件下沉陷DEM边界提取的合理临界坡度为0.2°，所提取的沉陷边界在主断上与实测结果基本一致。
(4)上述实验获取的沉陷DEM总体标准差在50 mm以内，对于下沉量达到米级的榆神矿区变形监测而言，基本能够满足其精度要求，但在沉陷边缘区域的相对精度不足。基于实验结果表明，去噪后的沉陷DEM中有15%的栅格下沉值误差在10 mm以内，具有很好的精度和可靠性。通过初步实验发现，这些高精度下沉栅格主要分布在平坦光滑的地表处，通过计算格网高程起伏度并构建合理的算法，便能提取出这些高精度的下沉栅格，以作为沉陷边界精准提取和沉陷盆地二次内插建模的依据，笔者所在团队正在展开探索。应该指出，本文实验结果是在榆神矿区特定地理环境条件获取的，在地形起伏较大的黄土沟壑区及植被覆盖度较高的其他矿区，机载激光扫描监测地表沉陷的可行性有待进一步探索研究。
[1] 钱鸣高,许家林,王家臣.再论煤炭的科学开采[J].煤炭学报,2018,43(1):1-13.
QIAN Minggao,XU Jialin,WANG Jiachen.Further on the sustainable mining of coal[J].Journal of China Coal Society,2018,43(1):1-13.
[2] 袁亮.我国深部煤与瓦斯共采战略思考[J].煤炭学报,2016,41(1):1-6.
YUAN Liang.Strategic thinking of simultaneous exploitation of coal and gas in deep mining[J].Journal of China Coal Society,2016,41(1):1-6.
[3] FAN Limin,MA Xiongde.A review on investigation of water-pre-served coal mining in western China[J].International Journal of Coal Science and Technology,2018,5(4):411-416.
[4] 汪云甲.矿区生态扰动监测研究进展与展望[J].测绘学报,2017,46(10):507-518.
WANG Yunjia.Research progress and prospect on ecological disturbance monitoring in mining area[J].Acta Geodaetica et Cartographica Sinica,2017,46(10):507-518.
[5] 芦家欣,汤伏全,赵军仪,等.黄土矿区开采沉陷与地表损害研究述评[J].西安科技大学学报,2019,39(5):859-866.
LU Jiaxin,TANG Fuquan,ZHAO Junyi,et al.Review of study on mining subsidence and ground surface damage in loess mining area[J].Journal of Xi’an University of Science and Technology,2019,39(5):859-866.
[6] 肖武,陈佳乐,赵艳玲,等.利用无人机遥感反演高潜水位矿区沉陷地玉米叶绿素含量[J].煤炭学报,2019,44(1):295-306.
XIAO Wu,CHEN Jiale,ZHAO Yanling,et al.Identify maize chlorophyll impacted by coal mining subsidence in high groundwater table area based on UAV remote sensing[J].Journal of China Coal Society,2019,44(1):295-306.
[7] 胡晓,李新举.基于无人机的高潜水位煤矿区沉陷耕地提取方法比较[J].煤炭学报,2019,44(11):3547-3555.
HU Xiao,LI Xinju.Comparison the extraction methods of subsided cultivated land in high-groundwater-level coal mines based on unmanned aerial vehicle[J].Journal of China Coal Society,2019,44(11):3547-3555.
[8] 刘军,王鹤,王秋玲,等.无人机遥感技术在露天矿边坡测绘中的应用[J].红外与激光工程,2016,45(S1):118-121.
LIU Jun,WANG He,WANG Qiuling,et al.Application of UAV remote sensing technology in open-pit slop mapping[J].Infrared and Laser Engineering,2016,45(S1):118-121.
[9] 张智韬,王海峰,韩文霆,等.基于无人机多光谱遥感的土壤含水率反演研究[J].农业机械学报,2018,49(2):173-181.
ZHANG Zhitao,WANG Haifeng,HAN Wenting,et al.Inversion of soil moisture content based on multispectral remote sensing of UAVs[J].Transactions of the Chinese Society for Agriculture Machinery,2018,49(2):173-181.
[10] 赵长森,潘旭,杨胜天,等.低空遥感无人机影像反演河道流量[J].地理学报,2019,74(7):1392-1408.
ZHAO Changsen,PAN Xu,YANG Shengtian,et al.Measuring streamflow with low-altitude UAV imagery[J].Acta Geographica Sinica,2019,74(7):1392-1408.
[11] REN He,ZHAO Yanling,XIAO Wu,et al.A review of UAV monitoring in mining areas:Current status and future perspectives[J].International Journal of Coal Science and Technology,2019,6(3):320-333.
[12] 于海洋,杨礼,张春芳,等.基于LiDAR DEM不确定性分析的矿区沉陷信息提取[J].金属矿山,2017(10):1-7.
YU Haiyang,YANG Li,ZHANG Chunfang,et al.Mining subsidence information extraction based on uncertainty analysis of LiDAR DEM[J].Metal Mine,2017(10):1-7.
[13] NGUYEN QUOC Long,MICHA
M.Buczek,LA PHU Hien,et al.Accuracy assessment of mine walls’surface models derived from terrestrial laser scanning[J].International Journal of Coal Science and Technology,2018,5(3):328-338.
[14] 卢遥,余涛,卢小平,等.基于高差分析的点云数据提取矿区地表沉陷信息方法[J].测绘通报,2013(3):22-25.
LU Yao,YU Tao,LU Xiaoping,et al.Subsidence information extraction of mine surface with point cloud data based on elevation difference analysis[J].Bulletin of Surveying and mapping,2013(3):22-25.
[15] 侯恩科,车晓阳,冯洁,等.榆神府矿区含水层富水特征及保水采煤途径[J].煤炭学报,2019,44(3):813-820.
HOU Enke,CHE Xiaoyang,FENG Jie,et al.Abundance of aquifers in Yushenfu coal filed and the measures for water-preserved coal mining[J].Journal of China Coal Society,2019,44(3):813-820.
[16] 汤国安,杨昕,等.ArcGIS地理信息系统空间分析实验教程[M].北京:科学出版社,2012:289-295.
[17] ZHANG Keqi,WHITMAN Dean.Comparison of three algorithms for filtering airborne LiDAR data[J].Photogrammetric Engineering & Remote Sensing,2005,71(3):313-324.
[18] JEFFREY S Evans,ANDREW T Hudak.A multiscale curvature algorithm for classifying discrete return LiDAR in forested environments[J].IEEE Transactions on Geoscience and Remote Sensing,2007,45(4):1029-1038.
[19] GEORGE Vosselman.Slope based filtering of laser altimetry data[J].International Archives of Photogrammetry and Remote Sensing XXXIII(Part B3),2000,935-942.
[20] ZHANG Keqi,CHEN ShuChing,D Whitman,et al.A progressive morphological filter for removing nonground measurements from airborne LiDAR data[J].IEEE Transactions on Geoscience and Remote Sensing,2003,41(4):872-882.
[21] P Axelsson.DEM generation from laser scanner data using adaptive TIN models[J].International Archives of the Photogrammetry,Remote Sensing and Spatial Information Sciences,2000,33(B4/1):110-117.
[22] 马国兵,肖培如.基于小波的图像去噪研究综述[J].工业控制计算机,2013,26(5):91-95.
MA Guobing,XIAO Peiru.Study on Wavelet-based image denosing[J].Industrial Control Computer,2013,26(5):91-95.
[23] 曾艺辉,高鸣.基于Bayesian估计的小波自适应阈值方法对图像进行去噪处理的研究[J].生物医学工程研究,2018,37(4):410-413.
ZENG Yihui,GAO Ming.Wavelet adaptive threshold method based on Bayesian estimation for image denoising[J].Journal of Biomedical Engineering Research,2018,37(4):410-413.
[24] 杜春梅,冀志刚,张琛.基于小波阈值法的矿山遥感图像非局部均值去噪[J].金属矿山,2017(3):116-120.
DU Chunmei,JI Zhigang,ZHANG Chen.Non-local means filtering algorithm of mine remote sensing image based on wavelet thresholding method[J].Metal Mine,2017(3):116-120.
[25] 刘晓莉,任丽秋,李伟,等.阈值优化的遥感影像小波去噪[J].遥感信息,2016(2):109-113.
LIU Xiaoli,REN Liqiu,LI Wei,et al.Threshold optimized wavelet for remotely sensed image denoising[J].Remote Sensing Information,2016(2):109-113.
[26] 谢家林,李根强,谢家丽,等.改进阈值函数在图像去噪中的应用[J].空军工程大学学报:自然科学版,2016,17(1):72-76.
XIE Jialin,LI Genqiang,XIE Jiali,et al.Research on the application of the improved threshold function to image de-noising[J].Journal of Air Force Engineering University(Natural Science Edition),2016,17(1):72-76.
[27] 于万波.基于MATLAB的图像处理[M].北京:清华大学出版社,2011:131-152.
