本发明属于煤矿区地面沉降监测技术领域，涉及一种采煤沉陷区insar影像快速识别方法。
背景技术
目前中国80％的电能及70％的能源总量来自煤炭，每天需要消耗原煤600万吨。煤炭资源的大规模开采，破坏了覆岩层的稳定性，导致地表的崩塌、沉降和塌陷，严重威胁社会经济稳定和公民人身安全。快速准确监测采煤沉陷区的范围及沉降速率可为生态文明建设及区域生态修复提供支持。
我国传统的煤矿区地面沉降监测技术目前仍主要依靠野外实地大地水准测量、gps测量和全站仪测量。通过布设离散点来监测地表变化，其监测精度可以达到毫米级。然而，对于区域范围内较大空间尺度的微小形变监测，大范围布设有效监测点并不现实。此外，传统观测结果只能反映监测点自身的形变量，而不是空间整体的观测结果。
insar技术具有覆盖面广、监测周期短、全天时、全天候等优点，可有效监测地表形变。但旧有的insar技术受时空失相干和大气效应等因素的限制，在精确监测采煤沉陷区地表沉降方面仍有较多缺陷。小基线集技术(smallbaselinesubset，sbas)在很大程度上解决了传统insar技术由于空间基线过长造成的失相干和大气效应等问题，并且增加了时间采样频率。该技术对基于低分辨率、大尺度、时序形变监测能够取得很好的效果。
sentinel-1(哨兵1号)卫星是欧洲航天局哥白尼计划(gmes)中的地球观测卫星，由两颗卫星组成，载有c波段合成孔径雷达，可克服各种天气干扰提供全时连续影像。可用于海洋环境监控和北极海冰监测、海洋溢油和舰船监视、地表形变监测、水体管理和土壤保护、森林火灾监测等。sentinel-1数据是目前现势性较好的免费sar数据，影像分辨率高，重访周期仅6天，非常适合地表形变监测。
技术实现要素：
本发明的目的在于提出了一种采煤沉陷区insar影像快速识别方法。该方法结合采煤沉陷区危害大、监测迫切性强的情况，根据sentinel-1卫星影像及sbas的技术特点，运用sbas-insar技术监测研究区地表时间序列的形变动态演变。
其技术方案如下：
一种采煤沉陷区insar影像快速识别方法，包括以下步骤：
步骤1、数据预处理：遴选研究期间相同极化方式的16景影像数据，根据研究区矢量范围裁切影像得到案例研究区，获取研究区的多时相slc数据；下载最新的覆盖研究区的srtm-3version4dem影像数据，以辅助研究区地形信息的生成；
根据数据类型，设置流程处理的相关基础参数，其中azimuthlooks参数设置为1，rangelooks设为4，cartographicgridsize设为20m，并将range/azimuthwindowsizefine参数设置为40，重采样方法为四次立体卷积法；
步骤2、连接图生成：考虑到采煤沉陷区分布及影响范围的特殊性，将16景影像空间基线设为临界基线的0-30％，时间基线设为0-760d，采用德洛尼处理方法进行3d解缠，冗余度设置为高。本步骤同时生成相关工程文件auxiliary.xml，记录进行的每一步处理，以及数据索引，作为后面每一步骤的输入文件；
步骤3、干涉处理：去平、滤波和相位解缠，其他影像数据与超级主影像进行配准，为轨道精炼和重去平以及为sbas关键步骤反演做准备；
步骤4、轨道精炼和重去平：选择具有代表性的解缠及差分干涉后的结果，输入地表高程数据，在地形平坦、没有相位跃变和形变条纹的区域均匀选取典型的gcp控制点，采用三次轨道精炼多项式估算轨道精炼和相位偏移量，消除斜坡相位，并基于控制点对所有数据进行重去平；
步骤5、sbas形变估算：本步骤实质为估算研究区的形变速率和残余地形；输入步骤4生成的工程文件，选择线性处理模型并生成形变速率和残余地形相关文件；
步骤6、sbas形变估算精炼：去除大气干扰并计算研究时间序列上的位移；在步骤5得到的形变速率基础上，进行定制的大气滤波，低通滤波量级为1200m，高通滤波为365d，精炼方式为轨道精炼，以估算和去除大气相位，得到更加精确的研究区时间序列最终位移结果；
步骤7、地理编码：以dem数据为辅助，将步骤6生成的工程文件编码到gcs-wgs-84坐标系，输出结果为si_vel_geo图像，其中高度精度误差控制在5mm/a，变化速率精度误差控制在8mm/a；
步骤8、结果分析：分析取得的结果图像，图像内的像元值即为该区域的平均纵向位移速率，通过对图像进行密度分割彩色渲染，能发现地表沉降区域异质性情况；通过时间序列栅格分析工具，能定位到采煤塌陷明显的区域，绘制研究时段内该区域的形变图，通过野外实地验证，其精度能达毫米级。
进一步，步骤3生成的结果包括各个影像对去平和滤波后的干涉图(*_fint)、相干系数图(*_cc)、解缠结果(*_upha)和斜距下的强度数据(*_slant_pwr)；能查看解缠结果和相应的相干性图，检查有无相干性低和解缠结果较差的像对，利用连接图边界工具将不好的像对进行剔除。
本发明在保证识别结果精度的前提下提供了影像遴选的方法，提供了采煤沉陷区的快速识别的具体参数及处理流程。
本发明的有益效果为：
与传统采煤沉陷区范围提取方法相比，sbas-insar技术具有更广泛的监测范围，可以适应大区域长时间尺度的地表精细化形变监测。与差分干涉技术相比，sbas-insar技术监测地表形变的分辨率更高，精度可以达到毫米级别，可以得到更为精细的监测结果。与地下水过量开采影响城市地表沉降监测的相关研究相比，采煤沉陷区面积趋于斑块状破碎化分布，但其影响范围更广，往往呈现跨地市级别分布，治理难度更大。基于sbas-insar技术可以快速准确的监测采煤沉陷地表的范围及垂直变化速率，适合推广应用。
附图说明
图1是德洛尼3d图；
图2是相关性结果图，其中，图2a为(*_cc)；图2b为(*_fint)；图2c为(*_slant_pwr)；图2d为(*_upha)；
图3是研究区形变图示例；
图4是研究区形变时间序列图示例。
具体实施方式
下面结合附图和具体实施方式对本发明的技术方案作进一步详细地说明。
本实施例结合采煤沉陷区危害大、监测迫切性强的情况，根据sentinel-1卫星影像及sbas的技术特点，以山东省济宁市为例，选取2015年7月至2017年7月影像，运用sbas-insar技术监测研究区地表时间序列的形变动态演变。该技术方法的具体流程如下：
①数据预处理：遴选研究期间相同极化方式(vv方式)的16景影像数据，根据研究区矢量范围裁切影像得到案例研究区，获取研究区的多时相slc数据；下载最新的覆盖研究区的srtm-3version4dem影像数据，以辅助研究区地形信息的生成。
根据数据类型，设置流程处理的相关基础参数，其中azimuthlooks参数设置为1，rangelooks设为4，cartographicgridsize设为20m，并将range/azimuthwindowsizefine参数设置为40，重采样方法为四次立体卷积法。
②连接图生成：考虑到采煤沉陷区分布及影响范围的特殊性，本方法将16景影像空间基线设为临界基线的0-30％，时间基线设为0-100d，采用德洛尼处理方法进行3d解缠，冗余度设置为高。根据本步骤生成的影像配对连接表结果，研究区2015-2017年地表形变监测的超级主影像为2016-05-31获取的雷达影像(图1)。
③干涉处理：去平、滤波和相位解缠，其他影像数据与超级主影像进行配准，为轨道精炼和重去平以及为sbas关键步骤反演做准备。本步骤用时较长，需设置较多参数。生成的结果包括各个影像对去平和滤波后的干涉图(*_fint)、相干系数图(*_cc)、解缠结果(*_upha)和斜距下的强度数据(*_slant_pwr)。可以查看解缠结果和相应的相干性图，检查有无相干性低和解缠结果较差的像对，利用连接图边界工具将不好的像对进行剔除。图2为以2016-05-13影像为主影像，2015-07-30影像为从影像所生成的相关结果图。
④轨道精炼和重去平：选择具有代表性的解缠及差分干涉后的结果，输入地表高程数据，在地形平坦、没有相位跃变和形变条纹的区域均匀选取典型的gcp控制点，采用三次轨道精炼多项式估算轨道精炼和相位偏移量，消除斜坡相位，并基于控制点对所有数据进行重去平。
⑤sbas形变估算：本步骤实质为估算研究区的形变速率和残余地形。输入前一步生成的工程文件，选择线性处理模型并生成形变速率和残余地形相关文件。
⑥sbas形变估算精炼：其核心是去除大气干扰并计算研究时间序列上的位移。在前一步得到的形变速率基础上，进行定制的大气滤波，低通滤波量级为1200m，高通滤波为365d，精炼方式为轨道精炼，以估算和去除大气相位，得到更加精确的研究区时间序列最终位移结果。
⑦地理编码：以dem数据为辅助，将上一步生成的工程文件编码到gcs-wgs-84坐标系，输出结果为si_vel_geo图像，其中高度精度误差控制在5mm/a，变化速率精度误差控制在8mm/a。
⑧结果分析：如图3、图4所示，分析取得的结果图像，图像内的像元值即为该区域的平均纵向位移速率，通过对图像进行密度分割彩色渲染，可以发现案例研究区范围内的地表沉降区域异质性情况。此外，通过时间序列栅格分析工具，可以定位到采煤塌陷明显的区域，绘制研究时段内该区域的形变图，通过野外实地验证，其精度可达毫米级。
以上所述，仅为本发明较佳的具体实施方式，本发明的保护范围不限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可显而易见地得到的技术方案的简单变化或等效替换均落入本发明的保护范围内。