兴仁县罗马山滑坡体形变监测及成因分析

章 彭，刘文明

(1.贵州省水利水电勘测设计研究院，贵阳 550002；2.中国测绘科学研究院 地球观测与时空信息科学国家测绘地理信息局重点实验室，北京 100830)

1 引言

作为一种常见的地质灾害，山体滑坡分布广泛，往往极大地危害附近居民的安全和财产。2010—2016年的统计结果表明，山体滑坡占中国所有地质灾害的70%以上(全国地质灾害通报 2010—2016年)，而我国地质灾害主要集中在南部及西南部地区。贵州省作为滑坡灾害的主要集中区(吕刚等 2016)，发生于2017年8月的纳雍县普洒镇塌方事件引起了国内外学者的广泛关注。

目前已有多种方法被联合用于滑坡体形变监测，如测地测量，基于全球定位系统的监测网络，以及遥感图像的解译。尽管这些方法获得了高度精确的结果，但由于人力和仪器成本，它们不容易获得高密度的测量点。作为有源微波传感器，合成孔径雷达(SAR)特别适用于贵州省这种多阴天和雨季地区。利用地形的先验知识，可以从在不同时间获取的两个图像中检测毫米变形。所谓的差分干涉SAR(D-InSAR)技术可以大规模地获得地表形变量，精度可达到厘米级。因此在地震，火山和采矿变形的监测等多种应用中得到了广泛的应用(Wang et al.2008；廖明生等 2013)。

本文选取贵州省黔西南布依族苗族自治州普安县罗马山作为研究区，共收集11景L波段的ALOS PALSAR2数据，生成10个干涉对，以形变时间序列的方式对研究区中的滑坡体进行形变监测分析。相对于X波段及C波段雷达信号，L波段具有较强的穿透能力，在黔西南地区这种多植被覆盖区域能够保证较小的时间去相干(田馨等 2013；Dong et al.2018；Zhao et al.2012)。文章最后结合该地区的降雨量信息和形变时间序列，发现降雨量于形变的高度相关性，这对滑坡灾害的预防具有重要意义(程海琴等 2014)。

2 研究区基本特征

图1 不稳定斜坡及矿区位置图

3 理论及方法

3.1 数据源简介

本文收集了2017年4月16日至2018年9月16日期间的11景ALOS PALSAR2雷达数据(见表1)，其工作模式为FBS(fine beam single polarization)，L波段(波长23.6 cm)，产品级别为Level1.1(单视复数影像数据)，HH极化，视角为32.8°，方位向像元大小为2.129 888 m，距离向像元大小为1.430 422 m，数据获取日期信息如表1所示。

表1 实验区ALOS PALSAR2数据

为了去除地形相位，在进行二轨差分法差分干涉测量时需要利用高精度DEM，本文采用美国地质调查局发布的30 m SRTM DEM数据。

3.2 二轨法D-InSAR技术

二轨法的基本思想是通过引入外部DEM去除干涉相位中的地形相位，实现过程中需要两景覆盖同一地区的SAR影像和辅助DEM(赵梦雪等 2017)。首先，需要将DEM与主影像进行配准处理；之后再利用DEM模拟成地形干涉条纹，得到地形相位；从主从影像生成的干涉相位中减去地形相位，即可得到形变相位；最后通过相高转换，实现相位到视线向形变量的转换。其公式可以表示为：

(1)

其数据处理的流程如图2所示。

图2 双轨法D-InSAR数据处理流程

(1)影像配准

目前，常用的配准方法有2种：①基于强度和相关系数法的配准方法其包括两景SAR影像之间距离向和方位向上初始偏移估计以及辅影像重采样到主影像，分为粗配准和精配准。精配准需要对距离向和方位向的偏移量多项式进行亚像元级的精确估计，并最终将配准精度限制在0.1个像元内。②基于查找表的配准方法。本文采用第一种配准方法。

(2)干涉图生成和相干性分析

主、辅影像精确配准后进行共轭相乘便得到了干涉图。在干涉图生成的同时，以计算相干系数作为评价干涉图质量的依据，由干涉系数产生的相干图直观地显示了干涉图的质量。相干系数的计算公式如下：

(2)

其中M，N分别为相干性估计是的窗口大小。相干系数值越大，形变监测值越精确。

(3)干涉基线估计

目前常用的基线估计方法有3种：①采用轨道信息估计基线；②采用快速傅里叶变换(Fast Fourier Transform，FFT)根据干涉图的条纹变化率估计基线；③基于地面控制点(Ground Control Point，GCP)估计基线。本文先根据方法1进行初始极限估计，然后使用方法2估计残余基线并将其加入初始估计基线中，以达到基线精化的目的。

(4)相位解缠

相位解缠是恢复雷达影像的原始相位值。目前常用的方法有2种：①枝切法；②最小费用流法(Minimum Cost Flow，MCF)。因为MCF基本不需要人工参与，本文选用后者。

(5)形变图生成及地理编码

地理编码是将相位解缠后的坐标转化为椭球坐标系下的坐标。这1步的关键是确定解缠后影像上每个像元对应的三维坐标。根据相位和形变的对应关系，从解缠后的形变相位转化为雷达视线向上的形变量。

4 实验结果

使用二轨法对收集的10景ALOS PALSAR2数据进行差分干涉处理，共生成9组干涉对，干涉对垂直基线及时间基线情况如表2所示。

表2 干涉对基本参数

续表2

编号主影像辅影像时间基线/d垂直基线/m52017102920171126288.0616201711262018021884-149.8647201802182018051384132.0938201805132018062442-178.9589201806242018080542212.464

根据第3.2节的处理流程对9个干涉对进行处理，最终得到9幅视线向形变量图，利用ENVI和ArcGIS等图像处理软件，将试验区中一处典型滑坡——罗马山滑坡制作成专题图，结果如图3所示。

图3 罗马山滑坡形变量图

需要说明的是，图中形变量的单位为m，而图中最大值(0.25)及最小值(-0.25)为人为设置，目的是为了统一图例。

同时，我们对罗马山滑坡体上一点的形变归算到型变速率，并进行了时间序列的统计。该点型变速率与降雨量叠加显示结果如图4所示。

图4 形变速率与降雨量关系图

5 结 论

从实验结果中可以看出，滑坡体的形变速率往往伴随着降雨量的增加而加快，形变速率与降雨量具有高度相关性。也就是说，强降雨是造成滑坡体形变的重要因素之一。同时，形变速率也是滑坡体稳定性评价的一个重要指标，滑坡发生前往往伴随着滑坡体的形变加速度，使用D-InSAR技术能够快速高效得获得滑坡体的形变量，对滑坡灾害的防治具有一定的应用价值。