空天地一体化监测联合反演开采沉陷概率积分预计参数研究

胡东升，程小凯，张雅飞，李 涛，廉旭刚

(1.华阳新材料科技集团有限公司，山西 阳泉 045000；2.自然资源部第一大地测量队，陕西 西安 710054；3.太原理工大学 矿业工程学院，山西 太原 030024)

煤炭开采为社会带来经济效益的同时，其开采将导致地表沉降，给矿区以及周围环境造成一定的负面影响。通过对开采地区长期、动态的沉陷监测，可及时获取矿区地表形变数据，分析移动变化规律，科学地指导矿区的开采活动，从而最大限度地降低因矿区开采带来的损失。以往采用的水准测量和GPS(Global Positioning System，GPS)测量等传统的矿区形变监测手段，成熟度虽高，但工作强度大、耗费人力物力、易受天气影响、工作周期长、监测点易受损害、监测范围小，不能实现对矿区持续、大范围的动态监测[1]。

合成孔径雷达差分干涉测量(Differential Interferometric Synthetic Aperture Radar，DInSAR)作为新兴的对地观测技术，能够实现全天时、全天候对研究区进行观测，且监测精度高，覆盖范围广[2]。SAR卫星获取的数据，不仅包括强度信息，还包含相位信息[3]。杨氏双缝干涉实验为InSAR技术的出现奠定了理论基础。1989年，Gabriel[4]首次使用DInSAR技术将干涉相位中的形变相位和地形相位分离，开启了InSAR监测地表形变的新篇章。1993年，Gabriel等DInSAR技术对地表形变进行监测，证实了其精度可达厘米级甚至毫米级。1996年，Carnec C[5]对ERS-1影像数据进行处理分析，初步证明了DInSAR技术可用于矿区地表沉降监测。国内学者经过多年的努力，在该领域也取得了重大突破。吴立新[6]等对覆盖中国东部典型矿区的5景ERS1/2影像进行差分处理，验证了DInSAR技术用于矿区沉降监测的可行性；杨嘉威[7]等引入分布式散射体 (Distributed Scatterer，DS) 目标，通过DS-InSAR技术对矿区地表和铁路沿线的形变规律进行了分析。杨泽发、朱建军等人[8-10]分别基于单轨InSAR、基于InSAR时序形变以及联合InSAR技术和水准数据等手段对矿区开采引起的形变进行监测，并对其沉降规律进行分析，为该领域做出了许多贡献。廉旭刚[11]等结合Sentinel1A和1B数据提高DInSAR的时间分辨率，实现对大同矿区沉降的高精度监测。

无人机最早出现在1917年，是为满足军事上的需求。Przybilla和Wester Ebbinghaus[12]分别于1979年和1980年成功进行固定翼和旋翼无人机的航空摄影测试。Pawe[13]等利用无人机摄影测量手段监测到矿区的不连续形变，促进了无人机技术应用于矿区监测的发展。21世纪后，国内加强对无人机摄影测量技术的研究[14]。2003年，盛业华[15]等进行地表塌陷变形试验并采用数字摄影测量对其观测。Zhou[16]等利用无人机摄影测量手段对煤矿区地表沉降进行监测，并反演了开采沉降参数。邱亚辉[17]等人利用无人机航测技术获取露天采坑的DEM数据，并对多期DEM进行差值，得到研究区相关变化量。戴嵩[18]等使用无人机技术得到尾矿坝正射影像，监测出最大沉降范围在0.16m之内。廉旭刚[19]等通过无人机摄影测量技术实现了对阳煤矿区的地表沉陷监测，证明了该技术在大梯度形变区域监测精度高，监测微小形变的精度有待改善。

目前，针对无人机与InSAR联合监测地灾的研究很少。大多是使用二者中的其一手段进行观测，借助另一种手段对监测结果进行空间验证，未能利用二者监测精度不同的优势。考虑到InSAR和无人机技术各自的优缺点，本文将两种技术的优势进行互补，融合二者的监测值，并使用融合后的下沉曲线对一矿81403工作面的概率积分参数进行反演。

1 研究区概况

华阳集团一矿81403工作面采用放顶煤综合机械化采煤法，走向长为1345m，倾向长为226m，煤层的平均倾角是4°，平均煤厚7.24m，平均采深为446.8m。研究区位置关系如图1所示。依据《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规范》的有关要求，结合本区内观测地表情况，采用剖面线观测线的方式在81403工作面上方布置了三条观测线，分别为半条走向观测线A线和两条倾向观测线B线和C线。根据开采深度和设站目的，并根据实际情况，观测线上的测点按30m距离埋设。设计各观测线的总长度、分段长度以及点间距和点数见表1。地表移动观测站的首次观测日期为2019年7月24日，末次观测日期为2021年11月23日，共进行了23次地表移动观测。本次研究主要对走向A线33个有效监测点以及倾向C线17个有效监测点进行分析。

图1 研究区相对位置

2 数据采集及处理方法

2.1 数据采集

研究选用的SAR影像为Sentinel-1A卫星的数据，该卫星由欧空局研发，并于2014年4月发射。卫星上搭载C波段的传感器，波长为5.6cm，分辨率为5m×20m，是一个全天时、全天候工作的雷达成像系统。通过ASF(Alaska Satellite Facility)下载平台，覆盖一矿81403工作面的单视复数(Single Look Complex，SLC)、升轨轨道的SAR影像数据。选取的数据时间跨度为2020年3月22日至2022年1月11日，并且为保证差分干涉结果的质量，以卫星的重访周期12d为间隔进行两两差分干涉，共55景影像。使用SRTM(Shuttle Radar Topography Mission)的30m分辨率DEM(Digital Elevation Model，DEM)作为差分的外部DEM。

表1 观测线长度、点间距和点数明细

无人机数据使用飞马智能航测系统D2000采集，搭载D-LIDAR2000激光雷达模块进行LiDAR点云数据的采集，并且保证与SAR数据时间上的统一。该采集系统联合飞马系统高精度的组合导航产品，测量精度可达5cm(50m航高)，当矿区地表沉降量为分米级时，该手段可满足沉陷监测的需求。实地采集数据时，依据研究区的地形起伏等情况进行航线规划，重叠度不低于50%。使用变高航线，当相对地面的飞行高度为177m时，可以采集到密度为64点/m2的点云。数据采集前的航线规划如图2所示。

图2 机载LiDAR航线规划

2.2 数据处理

2.2.1 DInSAR数据处理

使用GAMMA软件的Diff&Geo (Differential Interferometry and Geocoding) 差分干涉与地理编码处理模块。对原始SAR数据进行多视处理，多视比为4∶1，以此抑制影像的斑点噪声。将主影像与副影像两景数据进行配准后，使用外部DEM进行二轨差分干涉。通过自适应滤波和去平地效应等步骤，将形变相位外的其余相位去除，再进行相位解缠以及地理编码，经雷达视线向(Line of Sight，LOS)形变转竖直向形变得到最终的沉降结果。二轨法DInSAR处理流程如图3所示。

图3 二轨法DInSAR处理流程

处理数据时将数据按时间序列组成54个干涉像对，以保证每个干涉像对的时间间隔最小，最后利用ArcGIS软件将差分结果进行叠加得到2020年3月22日至2022年1月11日的差分干涉图，如图4所示，各时间为叠加的最终时间。

对差分干涉结果图进行分析时，会发现DInSAR手段在监测大梯度变形时会产生失相干。DInSAR技术虽然可以监测到整体下沉盆地的范围，在沉降盆地边缘处表达能力强，但是无法监测到盆地中间的大梯度变形，由于沉降盆地中心沉降量级大、突变等原因，造成干涉相位的不连续，出现失相干及监测结果的错误，使得盆地中心的监测值与实际下沉情况不符合。

2.2.2 无人机机载 LiDAR处理

通过机载 LiDAR 获取原始数据的预处理，即对采集结果进行数据解算与转换，生成标准las点云数据文件。采用ATIN算法对点云数据进行点云滤波、反距离加权作插值处理，得到2020年3月26日与2022年1月16两期相应的DEM数据。使用ArcGIS将两期DEM数据进行相减，得到该时间段内的无人机监测的地表沉降，结果如图5所示。

图4 DInSAR差分结果

采用无人机监测开采沉陷能够较为全面地反映出矿区沉陷影响范围，在对大梯度形变区域进行监测时，测量精度高，可以监测出最大下沉。但监测精度达不到传统监测毫米级精度，难以对矿区边缘进行高精度监测。

图5 无人机首尾两期DEM相减

3 概率积分法参数反演

3.1 融合DInSAR与无人机技术监测值

DInSAR技术在大量级、大梯度区域的监测精度差，更适用于沉降盆地边缘沉降值的监测。无人机技术监测沉降盆地中心大变形时精度较高，在沉降盆地边缘监测精度低。因此把研究区的无人机监测结果与DInSAR差分干涉得到的沉降值进行融合，得到更加完整、准确的监测结果。

在实际的采煤沉陷研究中，一般取沉陷值为10mm的点作为边界点。然而，DInSAR技术的观测精度可以达到毫米级。因此，为了获取矿区更完整的高精度监测，使用DInSAR手段监测沉降盆地边缘的小变形，使用无人机技术监测沉降盆地中心的大梯度变形。

DInSAR干涉的相干性决定其相位解缠结果的准确性，因此使用DInSAR平均相干系数对DInSAR监测值进行筛选。Baran等通过经验统计等方法对DInSAR可检测的形变梯度与相干性的关系进行了深入研究，得到了以下数学模型[20]。

Dmax=dmax+0.002(γ-1)

(1)

式中，Dmax为最大可监测到的形变梯度；γ为差分干涉图的相干系数，该值介于0到1之间；dmax为DInSAR理论上可监测到的最大形变梯度，等于卫星传感器波长与2倍影像分辨率之比。

由式(1)可看出相干性增大，可监测的最大形变梯度也随之增大。实际中受许多种去相关因素的影响，DInSAR可监测的形变梯度一般远小于理论值。对于采用 C 波段、分辨率为20m的Sentinel-1A的影像，通过式(1)可知，若DInSAR相干系数小于0.3，则其最大可监测的形变梯度等于0。

基于以上分析，选取相干系数低于0.3并最靠近研究区工作面沉降盆地边缘的监测点，以该点作为DInSAR监测结果的筛选阈值。通过提取各监测点的相干系数，经对比发现，从A9点之后的A12监测点开始出现相干系数小于0.3的情况，且越靠近沉陷盆地中心，相干性系数越低甚至出现失相干现象。因此确定对于A9号监测点以及其之前的点采用DInSAR监测值，其余监测点采用无人机监测值。对于C9号与C24号之间的监测点，其监测值采用无人机监测值，其余采用DInSAR监测值。

3.2 概率积分预计参数反演

当前开采沉陷预计使用较为广泛的一种方法为概率积分法(Probability Integral Method，PIM)，由我国学者刘宝琛、廖国华等在随机介质理论的基础上发展而得[21]。概率积分法共有5个与沉陷相关的参数，各自为下沉系数q、主要影响角正切tanβ，水平移动系数b，拐点偏移距Si(i=1，2，3，4)和开采影响传播角θ0[22]。一般通过在工作面上方布置观测站来获取实测数据，以此数据为基础确定概率积分法最佳参数[23]。

本次研究使用InSAR技术与无人机摄影测量技术的融合监测值来反演概率积分参数。准备好预计所需的工作面及观测线数据，根据工作面实际情况设置初始概率积分参数，采用《山区煤矿开采地表移动变形预计系统(MMSPS)》计算软件进行计算。将计算结果与融合后的曲线进行对比，分析结果的差异以调整预计参数，再次进行模拟计算，重复该步骤直到取得最佳的模拟效果，该结果下的预计参数就是最优模拟参数。经多次模拟计算得到的概率积分参数结果见表2。

表2 模拟求参计算结果

由于采用InSAR与无人机技术监测开始于2020年3月，而现场实测开始于2019年7月24日，为了使时间统一，对实测数据进行了处理，将2022年1月的累计实测下沉减去2020年3月的实测累计沉降，得到2020年3月至2022年1月的实测下沉值。做差后的实测下沉曲线与融合曲线整体趋势相同。

A线和C线融合下沉曲线结果与模拟下沉曲线结果对比如图6所示。

图6 81403工作面A线和C线融合下沉曲线与模拟下沉曲线对比

使用融合DInSAR监测值与无人机监测值后的下沉曲线，对一矿81403工作面概率积分法参数进行反演，反演结果为：下沉率q=0.8，水平移动系数b=0.3，主要影响角正切tanβ=2.5，开采影响传播角θ0=85°，拐点偏移距S1=S2=S3=S4=-20m。通过对比，反演的参数结果与实测数据所拟合的参数结果基本接近，证明了获得的概率积分法参数的合理性和可靠性。经过概率积分法求参，只得到了非充分采动情况下的预计参数。依据《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规程》所述，要获得充分采动时的预计参数，需要对非充分采动条件下的预计参数进行修正，其中水平移动系数和开采影响传播角与回采尺寸之间的关系不明显，可以不予修正。按照规程中下沉系数、主要影响角正切、拐点偏移距等参数与回采尺寸间的关系，再结合研究区覆岩性质，对概率积分预计参数进行修正。阳泉矿区覆岩性质为中硬，修正后的参数：下沉系数q=0.82，水平移动系数b=0.3，主要影响角正切tanβ=2.6，开采影响传播角θ0=85°，拐点偏移距S1=S2=S3=S4=-33.3m。该研究获取的DInSAR与无人机监测值均为一维下沉值，反演出水平移动系数这个参数的参考意义不大。需要说明的是，因为一矿81403工作面的观测站于2019年7月24开始观测，而无人机与InSAR技术监测沉降是于2020年3月开始，此时该工作面已经开采了300多米，所以造成A线预计下沉曲线比无人机与DInSAR融合后的下沉曲线靠右。

4 结 论

1)通过对研究区的Sentinel-1A数据进行差分处理，结果表明：由于沉降量、突变等原因使得影像出现失相干的现象，对DInSAR监测结果影响较大，沉降边缘沉降量小的点，误差小；而在沉降中心沉降量大的点，DInSAR在该区域获取的沉降值仍十分微小，未能有效监测出正确的下沉值。符合DInSAR监测微小形变的特点。

2)由无人机机载LiDAR得到的监测值结果表明：采用无人机监测可以监测出最大下沉，较全面的反映出矿区沉陷影响范围。但是难以对矿区边缘进行高精度监测，边缘表达能力差，无法监测微小沉陷变形。

3)融合DInSAR和无人机两种非等精度的监测值，一定程度上可以弥补DInSAR手段在大量级、大梯度形变区失相干的问题以及无人机技术在监测边缘微小形变时精度低的不足。并利用融合后的下沉曲线对一矿81403工作面的概率积分参数进行反演。结果表明，反演的参数结果与实测参数结果基本接近。该方法对矿区地质灾害评估、预防矿山地质灾害具有一定的参考价值。