时序InSAR的贵州省地质灾害监测
引用格式：高海英，赵争，章彭．时序ＩｎＳＡＲ 的贵州地质灾害监测 ［Ｊ］．测绘科学，２０２０， ４５（７）：９１—９９．
0 引言
贵州省是一个地质灾害发生尤其突出的省份，对旅游景区、重要工矿地区和人口集中区危害极大，滑坡、泥石流、采矿引起的塌陷以及尾矿堆的塌方是很多地区存在的主要安全隐患之一，造成严重经济损失和人员的伤亡；因此，对于这些危险区域的早期识别并持续形变监测具有十分重要的意义。遥感对地观测技术近年来不断发展，推动地质灾害隐患早期识别的工作进展，为防灾减灾做出巨大贡献[1]。合成孔径雷达干涉测量（interferometric synthetic aperture，InSAR）技术具有全天时、全天候、大范围、高分辨率对地观测能力，能够实现对微小地表形变的高精度探测[2-3]，近年来在业界大受好评，在地质灾害普查和监测中的巨大应用潜力在逐渐得到认知和肯定[4]。我国在雷达遥感地质灾害监测方面虽然起步较晚，但近年来经过众多学者的努力，得到了迅速发展。在2005年，文献[5]建立了GPS-CR（GPS-corner reflector）滑坡监测网，分析了InSAR观测技术与GPS观测技术在滑坡监测研究中的特点。同年，文献[6]对InSAR与GPS数据融合技术在矿区地表形变监测中应用的可行性进行了分析。2012年，文献[7]利用时间序列干涉图、相干图及强度图对美国加利福尼亚州北部和俄勒冈州南部的区域进行滑坡体制图，发现降水是影响当地滑坡形变的主要因素。2019年，文献[8]采用InSAR技术对甘肃黑方台地区2006年至2016年潜在的黄土滑坡开展了多时相编目、长时序监测以及失稳模式识别研究。地质灾害频频发生，滑坡、矿区塌陷、地表沉陷等发生原因并没有一个准确的标准，由于各个地区地形地貌、土壤结构、岩土类型、水文地质等条件的不同，造就了地质灾害发生原因的多可能性，本文以贵州省六盘水市水城县以及黔西南州贞丰县、安龙县为研究区，利用时序InSAR技术和全球降水观测计划（global precipitationmeasurement，GPM）中的全球降水数据以及实测数据，详细分析坡体失稳的形成原因。
对于山高谷深、绵延纵横的高原山地，具有植被覆盖率高、人工目标稀少等特点，采用C波段SAR数据和常规雷达差分干涉测量（differential InSAR，D-InSAR）方法受时空失相干影响严重，很难得到研究区域时间序列上可靠的监测结果[9]。应用传统的永久散射体干涉测量（persistent scatterer InSAR，PSI）技术也难以提取到足够数量和分布密度的稳定点目标，导致无法实现对地质灾害隐患的有效识别。本文利用小基线集（small baseline subsets，SBAS）InSAR方法，基于L波段ALOS-PALSAR 2数据，对贵州省黔西南州贞丰县和安龙县两个地区2017-2019年间的地质灾害隐患点进行早期识别，并结合实地调查数据、GPM全球降水数据分析其形变发展规律和空间分布特点以及形成原因。
1 实验区与实验数据
1.1 实验区
研究区位于贵州省黔西南州贞丰县和安龙县，地处云贵高原向广西低山丘陵过渡的斜坡地带，境内海拔高度最低324 m，最高2 572 m，相对高差2 248 m，境内地形起伏大、类型多，以山地为主，呈喀斯特地貌。研究区域如图1所示，黑色矩形为SAR影像覆盖范围，白色矩形为研究区的范围。两县矿产资源丰富，属金矿、煤矿等居多，研究区内包含两个地下矿，一个露天矿，随着矿山开采，以及每年汛期的到来，导致变形值不断增大，变形范围也随之扩大，随时会发生滑坡、泥石流以及崩塌等地质灾害，还伴随着地表沉降，对当地居民生命财产构成巨大威胁。最近几年逐渐利用多种监测手段对贵州省西南地区进行形变监测，包括外业实地调查、GPS测量等传统的监测方法，也包括InSAR技术进行短周期形变监测，相比前两种监测方法，InSAR技术不仅具有监测范围广、周期短的优势，而且也弥补了传统监测方法只能获取离散监测点的缺点。
图1 SAR数据覆盖范围
1.2 实验数据
选用的实验数据包括星载SAR数据、数字高程模型（digitalelevation model，DEM）、以及实地调查数据。星载SAR数据为12景L波段长波长的ALOS PALSAR-2降轨数据，表1列出了该数据的具体情况。本文收集的DEM数据为SRTM-1 DEM数据，其平面分辨率为30 m，即按经纬网中1″×1″格网大小进行采样，每个DEM子单元的数据幅宽为经纬度1°，3601×3601个像元。绝对水平精度为±20 m，绝对高程精度为±16 m。
本文还将引入全球降水计划GPM中的多星降水反演产品（integratedmulti-satellite retrievals for GPM，IMERG）。IMERG目前提供3套类型的卫星降水数据，分别是Early，Late和Final。本文选用的是Final产品，即称为IMERG-F产品，空间分辨率0.1°× 0.1°，时间分辨率为1个月，空间覆盖范围为（90°N～90°S），时间覆盖范围2000-06-01 ～2019-05-31[10-12]，从中获取了2017年4月至2019年5月的降水数据。以贵州省2017年9月1日的降水数据为例，如图红色代表降水量多，蓝色代表降水量少。研究人员需要准确、及时的降雨信息，以便更好的理解和建模何时何地可能发生严重洪灾、频繁的山体滑坡和破坏性干旱。相比上一代TRMM产品，GPM产品有着更大的覆盖范围和更高的时空分辨率，并增强了对微量和固态降水的探测能力，从而有效提高了探测精度[13]。IMERG-F产品在地形起伏较大区域具有较高精度，纳西效率系数为0.622，通过与该地区的历史天气情况对比，发现该降水数据与真实降水情况成正相关，满足本实验对降水数据精度的要求。GPM的全球降水数据有助于更好的准备和应对各种自然灾害[14-15]。该数据由美国国家航空航天局（National Aeronautics and Space Administration，NASA）提供，降水观测任务官方网站https://pmm.nasa.gov/precipitation-measurement-missions。
图3 提取贵州省2017-09-01至2017-09-30降水数据
2 SBAS技术原理与实验
2.1 SBAS技术原理
SBAS技术的基本思想是采用多主影像的方式，通过设置时间和空间基线阈值，组合出短时间和空间基线的干涉对。在很大程度上解决了由于空间基线过长而造成的失相干和大气效应等问题。SBAS技术的基本原理如下。
假设获取了监测区N+1幅雷达影像数据，首先选取一幅影像作为超级主影像，将其他N幅影像配准到超级主影像的成像空间。然后设置合适的时间、空间基线阈值，将符合阈值条件的干涉对进行差分干涉，生成M幅差分干涉图，M满足式（1）。
图4 时间序列InSAR技术处理流程
2.2 实验过程
利用ALOS PALSAR-2数据和SBAS技术，基于GAMMA软件进行贵州省贞丰县和安龙县的地质灾害监测具体实验过程包括数据读取、SAR影像配准、基线估计、D-InSAR处理、地理编码、获得形变速率图，下面将对关键技术进行具体说明：
1）SAR影像配准。设置2017-04-26为配准的主影像，其余11景影像均配准到该影像上，得到的影像在距离向和方位向上的偏移量均小于0.1个像元，达到了配准要求。
2）基线估计。为了得到较好的干涉图，实验中选取空间基线小于500 m的干涉对，共生成30对干涉对，最大空间基线437.919 4 m，干涉图组合如图5所示。
图5 ALOS PALSAR-2序列干涉图组合
3）常规D-InSAR处理过程。主要包括生成复数干涉图、生成SAR坐标系下的DEM，进行模拟地形相位、去除平地效应、自适应滤波，相位解缠使用最小费用流（minimum cost flow，MCF）方法，最后得到经过滤波、相位解缠的差分干涉图。
4）地理编码。将相干性图和差分干涉图从SAR坐标系下编码到国际通用地理坐标系GCS-WGS-1984下。图6列出了部分经过地理编码后得到的差分干涉图。
图6 地理坐标系下的差分干涉图
图7研究区雷达视线向平均形变速率图
3 实验结果分析
形变区主要分布在整景影像的西北方向，为了能够更加清晰的分析形变区在2017-04-26至2019-06-19时间段内的变化规律和分布状况，截取了形变区范围内的雷达视线向年平均形变量图，如图7所示。监测结果显示，共发现10处隐患点，有两处位于斜坡上，将重点分析，分别命名为P1和P2。如图8所示，是根据研究区的DEM数据提取的坡度图，P1点和P2点所在的斜坡坡度约在27.36°～50.88°之间，P1点最大坡度为40.59°，P2点附近最大坡度为44.64°，坡体较陡。另外有8处（P3-P10点）形变区域位于山顶平地处，极有可能是地下煤矿开采所导致的地表沉降，已形成明显的沉降漏斗，随时间的推移，在时间序列上沉降范围逐渐扩大，累积沉降面积达到0.372 km2。
图8 研究区坡度图
其中P1形变点位于龙山镇鑫发煤矿附近，该区域最大高程为1 815 m，最低高程为1 070 m，相对高差745 m，如图9（a）中光学影像所示，矿区山坡和山下均有居民点。该煤矿属地下煤矿，矿井主要开采煤层2层，自上而下编号为C3、C2煤层，煤层间距为80 m，目前开采C2煤层。通过对时序InSAR监测结果进行分析发现该区域最大形变速率约146 mm/a,靠近山体边缘处最大形变速率约120 mm/a。通过分析时间序列InSAR监测结果发现，2017-04-26至2018-06-20时间段内，靠近山体边缘处最大累积形变量约58 mm，而2017-04-26至2019-06-19时间段内，靠近山体边缘处最大累积形变量约254 mm，由此可以知道，在2018年6月20日这一时间节点山体发生剧烈形变，形变范围明显扩大，如图9（b）、图9（c）、图9（d）所示，形变方向及形变中心由山顶转移到山坡，已经构成地灾隐患。与实地调查数据相对照，该处存在崩塌隐患，发生形变的主要原因是人为的矿山开采，灾害孕育地层为基岩，险情等级已达到大型级别，严重威胁矿山厂区以及附近居民点。
图9 龙山镇鑫发煤矿累积形变和形变范围变化图
P2 形变点位于黔西南州贞丰县挽澜镇纳坎村某不稳定斜坡，该区域最大高程为1 518 m，最低高程为1 169 m，相对高差349 m。图10给出了利用ALOS PALSAR-2降轨数据集获取的P2点雷达视线向年平均形变速率图，该斜坡的形变方向是沿着坡顶向下发展到坡面，预测形变将逐渐发展到坡脚，最终可能发生滑坡。雨水容易使土壤饱和，强烈而长时间的降雨是最常见的滑坡触发因素[17]。雨水不断向坡体下渗，导致滑动面本身的重量加大，雨水长时间的浸泡使滑动面软化，逐渐降低了接触面的抗滑能力，当坡体承受能力达到极限时，将会发生滑坡。
图11 P2点隐患滑坡体时序形变与降水量关系
如图12所示，是P3-P10的形变速率图，通过提取下沉盆地边缘和中心沉降信息，得到了矿区在时间序列上的形变规律。观察并统计出了这8个沉降点的累积沉降面积和累积沉降量,如表2、表3所示，可以发现P3和P4点是累积沉降面积和累积沉降量最大的两个点，所以将以这两个点为例重点分析矿区沉降特点以及沉降变化规律。由统计数据可知P3和P4沉降点点累积沉降面积分别达到了0.136km2、0.374km2，累积沉降量分别达到了575.5mm、497mm，也就是说监测时间段内在水平方向和垂直方向都发生较大形变，与实地调查数据较为一致。根据地表沉陷盆地的特征可以知道，煤层埋深厚度越小，煤层开采厚度越大，对地表的影响越大。鑫发煤矿矿井开采区域在海拔+1414.32m以上，即上山开采，由此可知煤层埋深厚度约104m，属浅埋深，煤层厚度约3.22m，属中厚煤层。在这种情况下，地表极易受到地下开采的影响。图13显示的是进行地下煤层开采后，采空区覆岩移动分带情况，煤层开采形成采空区，放顶之后顶板垮落，形成垮落带，随着垮落带岩层的垮落松碎，填满采空区后，在地压作用下会产生大量裂缝，形成断裂带，但岩层整体状态基本不变，此时弯曲带会随着断裂带下沉而逐渐缓慢下沉，这样的沉降特点刚好符合时序InSAR监测到的监测结果，如图14所示。那么煤层埋深越浅，传递到地表的能量越大，对地表形变的影响也就越大，当能量全部传递到地表时，就会发生塌陷，造成附近房屋开裂、道路损坏等不可估量的后果。因此矿山开采所导致的地面沉降也是不容忽视的问题，采用InSAR技术监测矿区地表沉降问题既可以提前发现隐患点的具体位置，又可以为后期分析形变的发展规律提供有力依据。
表2 矿区沉降面积统计表
表3 矿区P3-P10点下沉盆地累积沉降量统计表
图14 P3和P4点时间序列沉降
4 结束语
本文以贵州省黔西南州贞丰县和安龙县为例，利用小基线集InSAR技术对该地区进行隐患点排查和形变监测。实验结果中共发现10处隐患点：j点P1和P2分别存在崩塌和滑坡的隐患，P1点位于陡坡上方，截止2019年6月山体边缘处累积形变量已达到254mm，且形变中心逐渐移向坡体下方，平均形变速率约120mm/a。P2点作为隐患滑坡点，利用全球降水数据分析它的变形规律，发现在监测日期内P2点日形变速率与月降水量有较高的相关性。说明该点形变受降雨影响严重。P1隐患崩塌点和P2隐患滑坡点的时间序列形变监测结果与实地调查数据表现一致。k统计了P3-P10沉降点的累积沉降量和累积沉降面积，并对典型的矿区沉降漏斗的形成原因进行分析，发现矿区开采造成了严重的地表沉降，其沉降量和沉降范围随着时间推移发生了巨大变化，P3和P4点是这8个点中变化最大的两个点。野外调查数据验证了P3、P4、P5、P6、P8点的时序InSAR监测结果确实发生了形变，由于没有收集到GPS或水准数据，不能进行定量验证。
通过对比和分析实地调查数据可以知道，InSAR技术对于隐患点的识别具有大范围、高精度的特点，可以识别出大部分形变的区域，并准确定位形变的位置，为提早预防提供数据支持。隐患滑坡点主要受强降雨影响，使用降水数据分析形变规律，达到了提前预警隐患点动态形变信息的目的，为防灾减灾提供依据。