基于时序InSAR的贵州喀斯特地区地表形变监测——以清镇市为例

王 硕，张正加，范 鹏,王猛猛

(1.自然资源部测绘发展研究中心，北京 100830；2.武昌理工学院 人工智能学院,湖北 武汉 430223；3.国网电力科学研究院武汉南瑞有限责任公司，湖北 武汉 430074)

0 引 言

地质灾害会对人民的生命安全和社会安全产生严重的威胁，严重阻碍社会的发展。我国自然灾害频发，是世界上地质灾害最严重、受威胁人口最多的国家之一，其中贵州省是我国受地质灾害影响最为严重的省份之一[1]。贵州的地质灾害主要以滑坡、泥石流、塌陷等地质灾害为主，对旅游景区和人口密集区危害极大。因此对贵州地区进行地表形变监测具有十分重要的意义。

虽然传统的地质灾害监测手段( GPS、水准测量等)能够有效地获取地面点形变信息，但其通常工作量大，且无法进行长时间、大范围观测。相对于传统地质灾害识别方法，遥感是一种现代高科技手段，具有客观、真实的特点，能为地质灾害调查监测提供重要的技术手段，能客观、真实、准确的提取地质灾害数据,为地质灾害的识别防治做出巨大贡献[2]。合成孔径雷达干涉测量(InSAR)技术作为雷达遥感的一个重要分支，具有较强的形变测量能力，其形变监测精度能够达到厘米级甚至毫米级[3-4]。InSAR技术具有监测范围广、精度高等优点，已经逐渐在地质灾害监测应用得到广泛认可[1,5]。为了克服传统InSAR技术固有的时间、大气去相干的影响[6]，研究学者提出时序InSAR方法，如永久散射体干涉测量(PSInSAR)[7]、小基线集(SBAS)[8]，并在城市沉降、地表形变、火山、地震等方面得到广泛应用[9-10]。近些年，时序InSAR技术已应用于山区地质灾害监测，特别是在四川、贵州、云南等地区得到应用[1,11-12]。

清镇市是贵阳市下辖县级市，位于贵州省中部，是贵州省地质灾害重灾区之一，目前清镇市地质灾害防治形势严峻，开展地质灾害监测是一项重要的工作[13]。因此，本文以贵阳清镇市为例，利用SBAS技术，使用26景Sentinel-1A数据，对贵州清镇市地区2018～2019年的地表形变监测，并识别出可能的地质灾害区域，并对其中形变严重地区进行详细分析。

1 清镇地理概况

清镇地处省会贵阳市以西，是贵阳市下辖的县级市，地处苗岭山脉北坡，乌江干流鸭池河东岸，东经106°07′～106°32′，北纬26°25′～26°56′之间[14]。东与白云区、乌当区毗邻，南与安顺市平坝区接壤，西为东风湖与织金县相连，北为鸭池河、猫跳河与黔西县、修文县隔岸相望，总面积1 383 km2。清镇地区植被覆盖茂密，覆盖面积达到了48%，如图1(a)所示。清镇市海拔变化从782 m至1 738 m，平均海拔1 200 m左右，如图1(b)所示。清镇市是贵州省地质灾害重灾区之一，高易发区主要分布在北部、西北部和中部的新店镇、流长乡、站街镇等乡。

图1 研究区位置范围

2 研究数据及处理

2.1 数据介绍

哨兵1号(Sentinel-1)卫星是欧洲航天局哥白尼计划(GMES)中的地球观测卫星，为双星系统，搭载C波段合成孔径雷达，可全天时、全天候不受天气条件影响地提供连续影像[15]。本文使用数据为Sentinel-1A卫星干涉宽幅模式的斜距单视复数产品(IW SLC)，图像的方位向和距离向分辨率为20 m和5 m，极化方式为VV。共收集了26景覆盖贵阳清镇市的哨兵数据，成像时间从2018年8月1日至2019年7月15日(表1)。清镇市的雷达图像如图2所示，在雷达图像中，水体表现为暗色，城市区域表现为亮色。本文收集了研究区的30 m分辨率的SRTM DEM数据，用于去除地形相位[16]。另外本文还收集了清镇市的Landsat-8影像，用于地质灾害识别与分析。

表1 SAR数据列表

图2 研究区SAR幅度图

2.2 SAR数据处理

研究区清镇市植被覆盖比较茂密，在干涉处理中容易引起时间去相干问题，导致反演结果精度较差。为了减小失相干影响，获得准确的地表沉降信息，本文使用小基线集(SBAS)的方法进行时序相干处理。本文的数据处理流程如图3所示，在本节重点介绍SBAS技术的关键步骤：干涉图生成、相干点选取和相位模型建立。

图3 数据处理流程图

2.2.1 干涉图生成

获取同一区域按照时间序列(t0,……,tN)排列的N+1幅SAR图像，选取其中的一幅影像作为主图像，并将其他SAR影像与主图像进行配准。根据干涉条件组合，得到M幅干涉图和相干图，以保证相干性的稳定。利用引入外部DEM，将干涉条纹图中的地形相位去掉得到差分干涉条纹图。为了抑制噪声的影响，使用自适应滤波方法对干涉相位进行滤波[17]。生成干涉条纹图以后，使用最小费用流的方法对差分干涉相位进行解缠[18]。

2.2.2 高相干点选取

选择适当的相干系数阈值，将N幅干涉图中，平均相干系数大于该阈值的点认为是相位稳定的点，即高相干点，如下式：

(1)

式中，ri表示单个像素在第i幅干涉图上的相干系数，thd_r是相干系数阈值。要注意，高阈值会舍弃许多有用的相干点，低的阈值又会带来许多噪声。在本文中，选择0.7作为相干点的选取阈值。

2.2.3 相位模型建立

假设从tA，tB两个时间获得的SAR图像产生第j幅干涉图，并假设tB>tA，去除平地和地形相位后，在x处的差分干涉相位[19]：

(2)

式中，d(tB,x)和d(tA,x)是相对于参考时间t1的LOS方向累计形变，d(t1,x)=0，φ(tB,x)和φ(tA,x)分别表示d(tB,x)和d(tA,x)引起的形变相位，λ是雷达波长，则利用如下线性模型估计N幅图像的形变：

AΦ=ΔΦ

(3)

式中，Φ表示待求点上的N个时刻的SAR图像上的未知形变相位组成的矩阵；ΔΦ为M幅差分干涉图上相位值组成的矩阵；系数矩阵A[M×N]每行对应于一幅干涉图，每列对应于一个时间上的SAR图像，主图像所在列为+1，辅图像所在列为-1，其余列为0。如果M>=N，且A的秩是N，则利用最小二乘法可得：

Φ=A#·ΔΦ,A#=(ATA)-1AT

(4)

上式给出了理想情况下的最优解。但实际上对于多个小基线集，ATA是一个奇异矩阵，如假设有L个子集，A的秩为N-L+1，方程组就会有无穷多解(设N≤M)。为了解决这个问题，SBAS算法采用矩阵的奇异值分解(SVD)技术从而求出方程(4)在最小范数意义上的最小二乘解。

将相位转化到平均相位速度：

(5)

代替式(2)中的相位，于是有：

∑(tk+1-tk)vk=Δφj,j=1,...,M

(6)

上述方程组也可表达为：Dv=ΔΦ，其中D是一个M×(N-1)矩阵，表示主辅图像之间干涉组合的时间基线，对其进行SVD分解可以得到速度矢量v的最小范数解。在线性形变估计的基础上，继续对残余相位进行合适的时间和空间滤波处理即能够分离出大气相位和非线性形变相位成分。

3 实验结果及分析

3.1 清镇市年平均沉降速率

对26景Sentinel-1A数据进行配准、预滤波等操作，设置相关的时间基线阈值(小于100 d)和空间基线阈值(小于60 m)，然后对生成的每个干涉图进行检查，剔除相位噪声较大的干涉对，最终得到78对干涉条纹图(图4)。然后根据计算得到的平均相干系数图，设置相关的阈值，得到研究区的高相干点。共选取了577 374个监测点，监测点密度达到了386/km2，选取的高相干点主要集中在城镇及周边区域，在山区选取的高相干点比较少。

图4 时间基线和空间基线分布

根据上述的方法和处理流程得到了清镇市201808～201907的年沉降速率结果，如图5所示。从图中可以看到，整体上清镇市的年沉降速率在-40～10 mm/year，大部分沉降速率在-10至10 mm/year区间，表明整体上比较稳定。沉降速率在空间上分布具有不均匀性，在清镇市的西部(犁倭镇)和中部(麦格苗族布依族乡)发现了沉降明显的区域，最大沉降速率达到了-40 mm/year。另外在卫城镇也监测到比较明显的沉降。

图5 清镇市201807～201907年平均沉降速率散点图

图6是根据沉降散点差值得到的清镇市沉降速率差值图，图中的线是沉降速率等值线。根据差值图，可以很明显地发现沉降区域分布的空间特征。沉降区域主要集中在清镇市的西北部和中部，在北部也零星分布了一些沉降区域，在这些区域应当加强地质灾害防治，避免在预计发生山体滑坡等地质灾害。

图6 清镇市沉降速率差值图

参考北京地区地面沉降危险性分级标准[20]，按照地面沉降速率大小将清镇分为地面沉降轻微区、一般区和较为严重区域，对应的沉降速率区间为：小于10 (mm·a-1),10～30(mm·a-1),大于30(mm·a-1)，统计结果如表2所示。从表2可见，沉降速率超过-30 (mm·a-1)的面积超过了0.42 km2，约占清镇市面积的0.03%。沉降速率在-30至-10mm/year的面积超过了54 km2，约占清镇市面积的3%,可以看到有明显沉降的区域面积比较大。

表2 沉降区间面积统计

3.2 区域分析

3.2.1 犁倭镇

本文提取的犁倭镇的平均沉降速率图如图7(a)所示，在犁倭镇的西南部分，出现了较为明显的沉降现象，最大的形变速率超多了-30 mm·a-1。在犁倭镇的东北则比较稳定，整个观测周期内的形变量小于5 mm·a-1。选取了3个点用于形变趋势分析，位置如图7(a)，时序形变如图7(c)。PS点1和点2的形变比较明显，并且形变趋势比较一致，观测周期内的形变分辨达到了33 mm和44 mm。从2018年8月至2019年5月形变速率比较缓慢，从2019年5月以后，形变速率明显加速。通常这种形变的加速分析可能与降雨量有关，特别是在山区[12]，这些区域需要注意，特别是在雨季，容易出现地质灾害隐患。PS3点在整个观测周期内则比较稳定。

(a)犁倭镇平均形变速率图；(b)犁倭镇光学图像；(c)选取3个点的累积形变量图7 犁倭镇沉降信息

3.2.2 麦格苗族布依族乡

麦格苗族布依族乡2018年8月至2019年7月的形变速率如图8(a)所示，可以看到在乡镇区域监测到了比较明显的沉降现象。我们选取了在形变区域选取了2个点(PS4和PS5)分析其时序形变规律。图8(c)展示了选取2个点的时间序列形变量，PS4和PS5两点的形变速率分别为25.4 mm/year和29.1 mm/year。需要主要的是，PS4和PS5两点的形变趋势比较一致，在观测周期内呈现线性变化，形变量分别达到了21 mm和31 mm。图8(d)是形变区域比较严重的光学放大图，我们可以看到在这些区域处于开发状态，地表植被遭到破坏，分析地表沉降可能与当地的矿产资源开采活动有关。

(a)麦格苗族布依族乡平均形变速率图；(b)麦格苗族布依族乡光学图像；(c)选取两个点的累积形变量；(d)形变区域的光学放大图图8 麦格苗族布依族乡沉降信息

3.2.3 卫城镇

图9(a)是卫城镇的平均沉降速率图，可以看到大部分区域都比较稳定。同样选取了PS6和PS7两个点对该区域进行分析。其中，PS6点为沉降比较明显的区域，PS7点位于稳定区域。可以发现，在PS6点的累积形变量达到了-20 mm，而PS7点的累积形变量在0～10 mm之间起伏。另外，从整体上看，围城镇的形变速率明显小于犁倭镇和麦格苗族布依族乡。

3.3 沉降因素分析

根据上述分析，可以看到清镇市的地表沉降主要集中在中部和西部。根据清镇市的地理环境，其沉降主要有3方面的因素：① 清镇市属于喀斯特岩溶地区，土地比较脆弱，土地容易石漠化[13-14],土地承载力较差，土层比较松散，在上层作荷载的作用下表面容易发生形变[12]；② 降雨的影响，雨水会进一步压实和冲刷松散土层，加速地表的形变。在犁倭镇PS1和PS2两点观测到了在雨季形变加速的现象；③ 人类活动的影响。清镇市有丰富的矿产资源，矿产资源的开采会引起比较明显的地表形变，如在麦格苗族布依族乡地区有丰富的矿产资源，当地的矿产开采活动一直在持续，一定程度上导致了该地区的沉降。

(a)卫城镇平均形变速率图；(b)卫城镇光学图像；(c)选取两个点的累积形变量图9 卫城镇沉降信息

4 结 语

本文利用SBAS-InSAR技术和26景Sentinel-1A数据，提取得到了贵阳清镇市2018年8月至2019年7月的地表沉降形变场。研究结果表明，研究区域内大部分地区表现较为稳定的形变，形变速率在-10 mm·a-1以内。沉降比较明显的地方主要集中在犁倭镇、麦格苗族布依族乡、卫城镇等3个乡镇地区，最大形变速率超过了-40 mm·a-1。沉降速率超过-10 mm·a-1的面积得到了54.7 km2,占全市面积的3%左右。在犁倭镇，发现了沉降信号在雨季有增加的趋势，分析这种形变增加趋势与降雨有关；另外发现了地质环境和人类活动有关等因素对形变的影响。

后期将继续利用SBAS技术对该地区的沉降持续监测，特别是地面沉降严重区域如卫城镇、犁倭镇、麦格苗族布依族乡，为清镇市的控沉规划的制定提供科学合理的建议。