利用Stacking/SBAS技术在滇西北地区滑坡隐患的识别对比

2023-09-02董继红梁京涛

山东科技大学学报(自然科学版) 2023年4期

董继红,张肃,梁京涛,杨磊,赵聪

(1.四川省综合地质调查研究所稀有稀土战略资源评价与利用四川省重点实验室四川成都 610081;2.四川省智慧地质大数据有限公司,四川成都610081)

滑坡是全球最常见、分布最广泛的地质灾害类型之一,不仅造成人员伤亡和基础设施破坏,还会形成链式灾害,产生二次破坏[1]。据统计,1995—2014年,全球共发生3 876 处灾难性滑坡,造成16万余人死亡和1.1万余人受伤[2]。准确识别潜在滑坡并绘制滑坡隐患分布图是防灾减灾工作的重点内容[3]。预先获知滑坡隐患的分布位置、发育特征和变形趋势,对预测潜在风险区和防灾治理尤为重要。目前滑坡隐患识别手段多基于光学影像和现场调查等,工作强度大且效率低,亟需寻找一种高效、高精度技术手段对滑坡隐患进行识别和监测。

合成孔径雷达干涉测量(interferometric synthetic aperture radar,InSAR)作为一种非接触获取地表形变的技术[4],具有全天时、全天候、穿云透雾等优点,在滑坡隐患识别中被广泛应用[5]。随着InSAR技术的不断发展,形变监测从起初精度较低且提取信息单一的D-InSAR技术,发展到后来的时间序列InSAR(time series interferometric synthetic aperture radar,TS-InSAR)技术[6]。TS-InSAR技术能利用较多SAR影像,最大限度克服时空去相干、大气延迟和数字高程模型(digital elevation model,DEM)误差等的影响,拓展了InSAR技术的应用领域[7-8]。特别是Sandwell 等[9]提出的Stacking技术,能通过较少SAR数据生成的差分干涉图进行相位堆叠,快速获取地表形变结果而被广泛应用于滑坡隐患的早期识别[10]。Berardino等[11]提出的小基线集(small baselines subset,SBAS)技术通过将较短时空基线干涉对组合,在提高监测点密度的同时获取监测点形变量。目前上述两种技术均被广泛应用于滑坡隐患识别[12],如Liang等[13]以川西片区为例,开展大范围Stacking技术与SBAS技术滑坡隐患识别对比分析,何佳阳等[14]利用InSAR和SAR技术开展雅砻江沿线西昌区域滑坡隐患识别。上述文献对比了不同InSAR技术在滑坡隐患识别中的应用效果,但未分析隐患识别的差异原因,而且在西南高山峡谷区利用不同InSAR技术进行隐患识别和对比分析尤为重要。

本研究利用Stacking和SBAS技术处理Sentinel-1数据以获取地表形变速率,并结合SAR数据的可视性进行隐患识别。将不同技术的隐患识别效果进行对比,详细分析两种InSAR技术在滑坡隐患识别中的差异性和影响因素。

1 研究区域

研究区位于云南省丽江市与怒江、大理两州的交界处(图1),属云贵高原与横断山脉结合部位,地势西北高,东南低。金沙江、澜沧江和怒江流经该区域,为泛三江源区域,地貌复杂多样,属于低纬高原季风气候。受断裂、构造的影响,该区滑坡频繁发生[15],使得当地居民的生命财产和安全长期受到威胁。因此,对该区滑坡隐患识别进行系统分析,不仅对该区域具有重要意义,也能为地质环境相似的西南山区滑坡隐患识别、治理提供参考。

2 数据和采用的技术方法

2.1 数据

收集研究区2018年11月—2021年12月期间92期Sentinel-1卫星SAR影像数据,如图1所示,由于Sentinel-1数据幅宽较大,需根据研究区大小对SAR影像进行裁剪,图1为数据覆盖范围。表1为所采用SAR数据集的基本参数。收集Sentinel-1卫星对应的POD(precise orbit ephemerides)精密轨道星历数据,用于去除因轨道误差引起的系统误差。引入AW3D30 DSM数据去除InSAR干涉处理中的地形相位并辅助SAR影像地理编码,该数据也被用来计算研究区的可视性分布情况。

2.2 技术方法

本研究基于Stacking和SBAS技术获取地表形变结果,因SAR数据为侧视雷达成像,在地形起伏区域存在几何畸变现象,故首先采用改进后的R指数进行可视性计算,分析地形特征对结果的影响。

2.2.1R指数

SAR卫星采用侧视雷达成像,因此在地形起伏较大的区域容易形成几何畸变,不仅降低了SAR影像的质量,还会对隐患识别造成漏判、误判,因此有必要准确识别几何畸变区域,以提高InSAR识别隐患的准确率[15]。常见的几何畸变现象包括叠掩、阴影和透视收缩。

SAR影像的几何畸变与雷达卫星的入射角、方位角、地面坡度、坡向有关,本研究采用由Notti等[16]提出的R指数模型,以及Ren等[20]为识别较远被动叠掩区域而提出的改进R指数公式计算SAR影像几何畸变区域。

R=sin{θ+arctan[tan(α)×cos(φ-β)]}×Sh×La×Fa。

(1)

其中:θ为卫星的入射角;φ为卫星的方位角;α为地面坡度;β为坡向;Sh为山体阴影系数,阴影区域取0,其他区域取1.0;La为顶底倒置(Layover)系数,主动顶底倒置和被动顶底倒置区域的值为0,其他区域的值为1.0;Fa为较远被动叠掩系数,主动叠掩区域和被动叠掩区域的值为0,其他区域的值为1.0。Sh、La、Fa三个系数由ArcGIS计算获得,通过计算可视性分析,可直观了解不同地形特征对结果的影响。

2.2.2 Stacking技术

Stacking技术通过对多幅差分干涉图加权平均获取地表形变速率结果。该技术的前提是假设地表形变趋势为线性形变,通过相位堆叠有效抑制大气相位和DEM误差,从而提高形变信息的精度[17]。该技术通过处理少量SAR数据快速获取大范围滑坡隐患分布图,结果可靠,被广泛应用于滑坡隐患的早期识别[18]。

2.2.3 SBAS技术

SBAS技术基于配准后的多幅SAR影像,根据设置的时间基线阈值和空间基线阈值生成小基线集,进行差分干涉处理,经滤波降低相位噪声增加高相干点位,再对滤波后的干涉图进行相位解缠,利用奇异值分解(singular value decomposition,SVD)小基线集得到形变速率,最后对形变速率积分得到监测时间段内的累积形变量[19]。

3 数据处理过程

利用GAMMA软件对所获取的92景Sentinel-1A数据进行InSAR数据处理,具体处理流程如图2所示。主要包括以下5个步骤:

1) 数据预处理:选择位于中间的一期SAR影像作为主影像,根据研究区范围对其进行裁剪,完成主影像与DEM的配准、裁剪,得到SAR坐标系下的DEM数据;同时完成主影像与其他SAR数据配准;限制时间基线不超过60 d,垂直基线不大于±250 m,生成时空基线,共获取380对组合干涉对。

2) 干涉工作流处理:对上述组合干涉对进行主辅影像共轭相乘得到差分干涉图;采用精密轨道数据去除平地相位;采用DEM数据模拟并剔除地形相位;采用自适应滤波的方法对差分干涉图进行滤波处理,以消除或减弱噪声,生成相干系数图;根据R指数计算阴影叠掩区域,同时因水域没有有效干涉信息,对差分干涉图进行掩膜处理,去除阴影叠掩区域和水域。采用最小费用流(minimum costflow,MCF)法进行相位解缠,使用相干性掩模避开相干性较低、相位不可靠的区域,设定0.3为相干值的阈值,低于此相干值的区域设为空值并不参与解缠计算。

3) Stacking计算对相位解缠结果进行相位堆叠计算,获取形变速率结果,然后将形变转换到LOS向,借助DEM数据进行地理编码获取地理坐标系下形变结果。

4) 时间/空间域变形估算

该步骤是在上述第2)步完成之后进行,①相邻点间参数估计:将点目标连接构成不规则三角网,依据点间连接关系求解相邻点差分相位差;②残余高程计算和线性变形:依据基线组合,估算相邻点间的线性变形速率和DEM误差;③残余相位低通滤波:从差分干涉相位中减去步骤①中差分相位得到残余相位,对残余相位进行空间域低通滤波得到滤波后的残余相位;④奇异值分解处理:根据短基线像对组合关系,对上一步得到的滤波后残余相位进行奇异值分解(Singular Value Decomposition,SVD)处理,求解每个影像对应时刻的大气相位和非线性变形相位;⑤大气相位和非线性变形相位计算:对奇异值分解得到的大气相位和非线性变形相位进行空间域高通滤波,得到大气相位,并对滤波后的相位序列进行时域低通滤波,得到非线性变形相位。

5) 形变量计算

将上一步获取的非线性形变结果和线性形变结果相加,根据时间基线参数获得形变量结果,将结果转换到视线向形变,利用DEM数据进行基准修正和地理编码,获取地理坐标系下形变速率和时间序列结果。

4 结果分析

4.1 可视性分析

假设SAR卫星为太阳所在位置,利用表1中Sentinel-1数据参数,通过ArcGis软件对卫星高度角和方位角计算Sh、Fa和La系数值,然后通过式(1)计算得到最终的SAR地形可视性分布图(图3)。当R>sinθ,为好可视性区域;当0

(a)地形可视性分布;(b)几何畸变面积统计

4.2 地表形变结果分析

利用Stacking和SBAS技术分别对2018年11月—2021年12月的Sentinel-1数据进行处理,得到研究区雷达视线向形变速率(图4),正值表示地表形变靠近卫星飞行方向;负值表示地表形变背向卫星飞行方向。分别利用Stacking和SBAS技术获取的LOS向形变速率为-83～51和-79～29 mm/yr。对形变参考区域的平均值和标准差进行统计分析,确定以-10～10 mm/yr作为相对稳区,该阈值外的速率对应不同程度的形变,值越大表示形变越强烈,其中形变最大区域主要集中在河谷两侧。通过对比监测结果发现,SBAS技术获取的结果点位较稀疏,究其原因一方面受到几何畸变的影响,另一方面由某些区域相干性不连续导致。

图4 沿LOS向平均形变速率图

4.3 相关性分析

对两种不同InSAR技术获取的Stacking和SBAS结果进行相关性分析(图5),发现两者线性相关性系数达0.72,表明两种方法得到的结果具有较高的一致性。通过相关性分析,一方面验证两种技术获取结果的一致性,另一方面可验证内符合精度。

图5 Stacking与SBAS技术年平均形变速率相关性图

4.4 地质灾害隐患识别对比分析

为了对比Stacking和SBAS两种技术对地质灾害隐患的识别能力,借助该区域的DEM数据和光学影像数据,分别对结果进行解译和统计。Stacking和SBAS技术分别识别出隐患点45处和36处,结合野外验证得到滑坡隐患分布图(图6),图中绿色圆点为Stacking和SBAS技术共同识别的滑坡隐患,红色圆点为仅被Stacking技术识别出来的隐患点。表2为不同InSAR技术识别的滑坡隐患数目,可以看出,Stacking技术识别的滑坡隐患数目较多,但SBAS技术识别隐患的准确率更高。

图6 不同技术识别滑坡隐患分布图

为分析两种技术识别滑坡隐患差异的原因,对识别的32处隐患点的最大LOS向形变速率和R指数进行统计(表3)。由表3可以看出,能被两种技术共同识别的隐患点中,SBAS和Stacking技术能够识别的最小形变速率分别为14.8、15.6 mm/yr,且R指数均大于0.8,属于好可视性区域;仅被Stacking技术识别出的隐患点最小形变速率为14.2 mm/yr,R指数为0.52～0.63,属于透视收缩区域。

表3 不同InSAR技术识别滑坡特征统计表

为进一步分析两种技术识别的差异性,选择两处典型滑坡点进行分析:第一处滑坡点是被两种技术共同识别的中村滑坡;第二处滑坡点为仅被Stacking技术识别的营盘镇滑坡。

中村滑坡位于兰坪县,坡向朝东。从光学影像(图7(a))可以看出植被覆盖较差,滑坡前缘以澜沧江为界,由Stacking和SBAS结果可以看出滑坡体存在两处明显变形区域,即滑坡的中前部右侧和中后部左侧。对比两种技术可知,Stacking结果的有效监测点更密集,SBAS结果形变速率更大。为进一步分析不同InSAR技术的监测结果和R指数值,沿图7(b)中的剖面线提取形变曲线,如图8(a)所示,其中图7(c)中红色箭头标注区域对应图8(a)中灰色阴影区域,两种技术识别的变形区域一致,形变速率均超过-15 mm/yr,R指数均大于0.8。因此,能被两种技术共同识别的条件是形变速率大于-10 mm/yr且斜坡处于好可视性区域。另外,造成SBAS结果点位稀疏的主要原因是部分区域的相干性在时间域不连续。而SBAS技术获取的形变速率更大是因为Stacking技术是对形变相位的加权平均,适用于线性形变,因此监测过程会丢失非线性形变信息。

营盘镇滑坡位于兰坪县,坡向朝西,光学影像如图9(a)所示。SBAS结果无明显形变信息且监测点位稀疏,在R指数图中绝大部分区域为透视收缩区域。为深入分析不同InSAR技术和R指数对结果的影响,沿图9(b)中黑色剖线提取形变曲线,图9(c)中红色箭头区域对应图8(b)灰色阴影区域,两种技术在变形区域均出现明显的变形趋势,但SBAS结果中仅有4个像素值超过10 mm/yr,且最大值不超过12 mm/yr;在Stacking结果中阴影区域速率均超过10 mm/yr,部分区域变形速率超过12 mm/yr,R指数小于0.6,属于透视收缩区域。因此,该区域未被SBAS技术识别的原因是坡向朝西,在升轨数据中处于透视收缩区域,造成在地理编码之后监测点位稀疏,形变信息丧失;在Stacking结果中有明显形变信息的原因是两者求解形变速率的原理不同。但根据SAR数据成像几何关系可知该隐患点在降轨数据中应处于好可视性区域,可以弥补SBAS结果在该区域识别的劣势,增加SAR数据的可观测区域。

为证明降轨数据可以提高升轨数据的好可视性范围,对营盘镇滑坡的降轨Sentinel-1数据进行处理,并提取该滑坡体不同位置的3个时序监测点进行变形特点分析(见图10中的P1、P2、P3),对比图9(d)可以发现,降轨SBAS结果中变形特征明显,最大形变速率为-38 mm/yr,表明降轨数据可以提升升轨数据在透视收缩区域的监测效果。图11为营盘镇滑坡时间序列曲线图,发现该滑坡中前部(点P2、P3)变形强烈,且变形趋势接近,最大形变量为-83.3 mm,目前仍处于持续变形状态,存在较大隐患。